#include "CXGammaSearch.h"

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>

#include "TString.h"
#include "TSystem.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TMatrixD.h"
#include "TGListBox.h"
#include "TGNumberEntry.h"
#include "TGLabel.h"
#include "TGButton.h"
#include "TList.h"
#include "TGResourcePool.h"

#include "GwConfig.h"
#include "LevelScheme.h"
#include "GammaLink.h"

#include "CXMainWindow.h"
#include "CXNucleus.h"
#include "CXLevelSchemePlayer.h"
#include "CXSpreadIntensityMatrix.h"

using namespace std;

CXGammaSearch::CXGammaSearch(const TGWindow *p, const TGWindow *main, UInt_t w, UInt_t h, CXMainWindow *mwin) :
    TGTransientFrame(p, main, w, h, kHorizontalFrame),
    fMainWindow(mwin)
{
    fDataBaseFolder = Form("%s/DataBase/LevelScheme/ENSDF/",getenv("GWSYS"));

    fLevelScheme = new LevelScheme;
    fSp = new CXSpreadIntensityMatrix;

    TGCompositeFrame *Main = new TGCompositeFrame(this,600,100,kVerticalFrame);

    TGLayoutHints *GroupHints = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 5, 5, 0, 0);

    AddFrame(Main,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandY,10,10,10,10));

    //Gammas Energy

    TGGroupFrame *gFrame = new TGGroupFrame(Main, "Input Gammas rays", kVerticalFrame);
    gFrame->SetTextColor(CXblue);
    Main->AddFrame(gFrame, GroupHints);

    TGCompositeFrame *hFrame = new TGCompositeFrame(gFrame, 60, 20, kHorizontalFrame);
    hFrame->AddFrame(new TGLabel(hFrame, "Energy (keV)"),new TGLayoutHints(kLHintsTop | kLHintsLeft, 21, 0, 3, 0));
    hFrame->AddFrame(new TGLabel(hFrame, "Width (keV)"),new TGLayoutHints(kLHintsTop | kLHintsLeft, 35, 0, 3, 0));
    gFrame->AddFrame(hFrame,new TGLayoutHints(  kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));

    ///gammas
    for(int i=0 ; i<3 ; i++)
    {
        gFrame->AddFrame(hFrame = new TGCompositeFrame(gFrame, 60, 20, kHorizontalFrame),new TGLayoutHints(  kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
        hFrame->AddFrame(fCheckGammas[i] = new TGCheckButton(hFrame, "", 0), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,0,5,0,0));
        hFrame->AddFrame(fEnergies[i] = new TGNumberEntry(hFrame, 500.0, 8, 21,TGNumberFormat::kNESRealOne, TGNumberFormat::kNEANonNegative),new TGLayoutHints(  kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
        hFrame->AddFrame(fWidths[i] = new TGNumberEntry(hFrame, 1.0, 5, 21,TGNumberFormat::kNESRealOne, TGNumberFormat::kNEANonNegative), new TGLayoutHints(  kLHintsExpandX | kLHintsExpandY, 20, 0, 0, 0));
        fCheckGammas[i]->Connect("Clicked()", "CXGammaSearch", this, "HandleButtons()");
    }
    fCheckGammas[0]->SetState(kButtonDown);

    // Selected range

    gFrame = new TGGroupFrame(Main, "Selected range", kVerticalFrame);
    gFrame->SetTextColor(CXblue);
    Main->AddFrame(gFrame, GroupHints);

    hFrame = new TGCompositeFrame(gFrame, 60, 20, kHorizontalFrame);
    hFrame->AddFrame(fZRange[0] = new TGNumberEntry(hFrame, 1, 6,0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative),new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX,1,1,0,0));
    hFrame->AddFrame(new TGLabel(hFrame, " < Z < "),new TGLayoutHints(kLHintsCenterY, 0, 0, 0, 0));
    hFrame->AddFrame(fZRange[1] = new TGNumberEntry(hFrame, 100, 6,0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative),new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX,1,1,0,0));
    gFrame->AddFrame(hFrame,new TGLayoutHints(  kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));

    hFrame = new TGCompositeFrame(gFrame, 60, 20, kHorizontalFrame);
    hFrame->AddFrame(fNRange[0] = new TGNumberEntry(hFrame, 0, 6,0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative),new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX,1,1,0,0));
    hFrame->AddFrame(new TGLabel(hFrame, " < N < "),new TGLayoutHints(kLHintsCenterY, 0, 0, 0, 0));
    hFrame->AddFrame(fNRange[1] = new TGNumberEntry(hFrame, 100, 6,0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative),new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX,1,1,0,0));
    gFrame->AddFrame(hFrame,new TGLayoutHints(  kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));

    hFrame = new TGCompositeFrame(gFrame, 60, 20, kHorizontalFrame);
    hFrame->AddFrame(fARange[0] = new TGNumberEntry(hFrame, fZRange[0]->GetNumber()+fNRange[0]->GetNumber(), 6,0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative),new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX,1,1,0,0));
    hFrame->AddFrame(new TGLabel(hFrame, " < A < "),new TGLayoutHints(kLHintsCenterY, 0, 0, 0, 0));
    hFrame->AddFrame(fARange[1] = new TGNumberEntry(hFrame, fZRange[1]->GetNumber()+fNRange[1]->GetNumber(), 6,0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative),new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX,1,1,0,0));
    gFrame->AddFrame(hFrame,new TGLayoutHints(  kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));

    fZRange[0]->Connect("ValueSet(Long_t)", "CXGammaSearch", this, "HandleButtons()");
    fZRange[1]->Connect("ValueSet(Long_t)", "CXGammaSearch", this, "HandleButtons()");
    fNRange[0]->Connect("ValueSet(Long_t)", "CXGammaSearch", this, "HandleButtons()");
    fNRange[1]->Connect("ValueSet(Long_t)", "CXGammaSearch", this, "HandleButtons()");
    fARange[0]->Connect("ValueSet(Long_t)", "CXGammaSearch", this, "HandleButtons()");
    fARange[1]->Connect("ValueSet(Long_t)", "CXGammaSearch", this, "HandleButtons()");

    //Finder

    gFrame = new TGGroupFrame(Main, "Gamma search", kVerticalFrame);
    gFrame->SetTextColor(CXblue);
    Main->AddFrame(gFrame, GroupHints);

    fStartButton = new TGTextButton(gFrame, "Start");
    fStartButton->Connect("Clicked()", "CXGammaSearch", this, "SetCalMode()");
    gFrame->AddFrame(fStartButton,new TGLayoutHints(kLHintsCenterX | kLHintsCenterY | kLHintsExpandX,15,15,10,10));

    fNoCoincMode = new TGRadioButton(gFrame,"Only Gammas",M_WithOut_Cascade);
    fCoincMode = new TGRadioButton(gFrame,"In Coincidence",M_With_Cascade);
    fNoCoincMode->SetState(EButtonState::kButtonDown);
    fCurrentMode=M_WithOut_Cascade;

    fCoincMode->Connect("Clicked()","CXGammaSearch", this, "HandleButtons()");
    fNoCoincMode->Connect("Clicked()","CXGammaSearch", this, "HandleButtons()");

    gFrame->AddFrame(fNoCoincMode,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY,10,10,10,10));
    gFrame->AddFrame(fCoincMode,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY,10,10,10,10));


    //Results panel

    Main = new TGCompositeFrame(this,400,100,kVerticalFrame);
    TGLabel *lab = new TGLabel(Main,"Results");
    lab->SetTextColor(CXred);
    Main->AddFrame(lab,new TGLayoutHints(kLHintsTop | kLHintsLeft,0,0,0,0));

    fNNucAnalysed = new TGLabel(Main,Form("Analysed nuclei: %5.0d",0));
    fNNucAnalysed->SetTextColor(CXblue);
    Main->AddFrame(fNNucAnalysed,new TGLayoutHints(kLHintsTop | kLHintsLeft,0,0,0,0));
    fNGrayAnalysed = new TGLabel(Main,Form("Analysed gamma rays: %5.0d",0));
    fNGrayAnalysed->SetTextColor(CXblue);
    Main->AddFrame(fNGrayAnalysed,new TGLayoutHints(kLHintsTop | kLHintsLeft,0,0,0,0));

    fResultsBox = new TGListBox(Main);
    Main->AddFrame(fResultsBox, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,0,0,0,0));

    AddFrame(Main,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY,10,10,10,10));

    HandleButtons();

    SetCleanup(kDeepCleanup);
    SetWindowName("GammaSearch");
    CenterOnParent();
    MapSubwindows();
    Layout();
    MapWindow();
}

CXGammaSearch::~CXGammaSearch()
{
    fMainWindow->fGammaSearchWindow = nullptr;

    UnmapWindow();
    CloseWindow();
}

void CXGammaSearch::SetCalMode()
{
    switch (fCurrentMode)
    {
    case M_With_Cascade:
    {
        if(fNGammas<=1)
            cout<<"Soory, I cannot determine coincidence between one or less gamma ray..."<<endl;
        else if(fNGammas==2)
        {
            FindGammaRays(true);
            FindInDoubleCoincidence();
        }
        else
        {
            FindGammaRays(true);
            FindInTripleCoincidence();
        }
        break;
    }
    case M_WithOut_Cascade:
    {
        if(fNGammas>0)
            FindGammaRays();
        else
            cout<<"Soory, no requested energies..."<<endl;
        break;
    }
    }
}

void CXGammaSearch::FindGammaRays(Bool_t Bash)
{
    Int_t AnalysedNuclei=0;
    Int_t AnalysedGammaRays=0;

    fResultsBox->RemoveAll();

    Int_t ZMin = fZRange[0]->GetNumber();
    Int_t ZMax = fZRange[1]->GetNumber();
    Int_t AMin = fARange[0]->GetNumber();
    Int_t AMax = fARange[1]->GetNumber();

    Float_t Gates[fNGammas];
    Float_t Width[fNGammas];

    for(int ig=0 ; ig<fNGammas ; ig++)
    {
        Gates[ig] = fEnergies[ig]->GetNumber();
        Width[ig] = fWidths[ig]->GetNumber();
    }

    fListOfGoodGammas.clear();

    for(Int_t iz=ZMin ; iz<=ZMax ; iz++)
    {
        for(Int_t ia=AMin ; ia<=AMax ; ia++)
        {
            AnalysedNuclei++;

            CXNucleus Nuc(iz,ia);
            TString MajName = Nuc.GetSymbol();
            MajName.ToUpper();
            MajName.Append("NoVerbose");


            TString FileName = fDataBaseFolder + "/" + ((TString)Nuc.GetSymbol()) + ".ens";

            if(!gSystem->IsFileInIncludePath(FileName))
                continue;

            fLevelScheme->Clear();
            fLevelScheme->Import(FileName,MajName);

            if(fLevelScheme->GetLinks().GetSize() == 0)
                continue;

            vector<Gw::GammaLink *> GoodGammas;

            Bool_t gatetrig[fNGammas];
            memset(gatetrig,0,sizeof(fNGammas));

            for(Int_t i=0 ; i<fLevelScheme->GetLinks().GetSize() ; i++)
            {
                AnalysedGammaRays++;

                Gw::GammaLink *Link= (Gw::GammaLink*)fLevelScheme->GetLinks().At(i);

                Float_t E = Link->GetEnergy().GetValue();

                for(int ig=0 ; ig<fNGammas ; ig++)
                {
                    if(TMath::Abs(E-Gates[ig])<Width[ig])
                    {
                        gatetrig[ig] = true;
                        GoodGammas.push_back(Link);
                    }

                    if(gatetrig[ig])
                        continue;
                }
            }

            Int_t NGatesTrig=0;
            for(int ii=0 ; ii<fNGammas ; ii++)
                NGatesTrig += gatetrig[ii];

            if(NGatesTrig>=fNGammas)
            {
                //                cout<<NGatesTrig<<" "<<fNGammas<<endl;
                //                cout<<GoodGammas.size()<<endl;
                //                cout<<Form("%s (Z=%d, A=%d, N=%d)",Nuc.GetSymbol().Data(),Nuc.GetZ(),Nuc.GetA(),Nuc.GetN())<<endl;

                if(!Bash)
                    PrintInListBox(Form("%s (Z=%d, A=%d, N=%d)",Nuc.GetSymbol().Data(),Nuc.GetZ(),Nuc.GetA(),Nuc.GetN()),kInfo);

                std::vector< GammaTransition > avect;

                for(int ig=0 ; ig<GoodGammas.size() ; ig++)
                {
                    GammaTransition agammatrans;

                    GammaLink *Link = GoodGammas[ig];
                    NuclearLevel *NucLevI = (NuclearLevel*)Link->GetIL();
                    NuclearLevel *NucLevF = (NuclearLevel*)Link->GetFL();

                    Float_t Energy = Link->GetEnergy().GetValue();

                    Float_t ELevI = NucLevI->GetEnergy().GetValue();
                    TString spinI = GetSpinFromLev(NucLevI);
                    Float_t ELevF = NucLevF->GetEnergy().GetValue();
                    TString spinF = GetSpinFromLev(NucLevF);

                    Float_t LifeTime = NucLevI->GetT().GetValue();
                    TString LifeTimeText="";
                    if(LifeTime != 0) {
                        if(LifeTime<1e-12)
                            LifeTimeText = Form(" ; T1/2 = %.1f fs",LifeTime*1e15);
                        if(LifeTime>1e-12 && LifeTime<1e-9)
                            LifeTimeText = Form(" ; T1/2 = %.1f ps",LifeTime*1e12);
                        if(LifeTime>1e-9 && LifeTime<1e-6)
                            LifeTimeText = Form(" ; T1/2 = %.1f ns",LifeTime*1e9);
                        if(LifeTime>1e-6 && LifeTime<1e-3)
                            LifeTimeText = Form(" ; T1/2 = %.1f us",LifeTime*1e6);
                        if(LifeTime>1e-3 && LifeTime<1e-0)
                            LifeTimeText = Form(" ; T1/2 = %.1f ms",LifeTime*1e3);
                        if(LifeTime>1e-0)
                            LifeTimeText = Form(" ; T1/2 = %.1f s",LifeTime);
                    }

                    agammatrans.NucName = Nuc.GetSymbol();
                    agammatrans.EGamma = Energy;
                    agammatrans.EI = ELevI;
                    agammatrans.EF = ELevF;
                    agammatrans.SpinI = spinI;
                    agammatrans.SpinF = spinF;
                    agammatrans.LifeTime = LifeTimeText;

                    avect.push_back(agammatrans);

                    TString GammaTitle = Form(" => %6.1f keV : %5s (%6.1f keV) --> %5s (%6.1f keV)%s",Energy,spinI.Data(),ELevI,spinF.Data(),ELevF, LifeTimeText.Data());
                    if(!Bash) PrintInListBox(GammaTitle.Data(),kPrint);
                }

                fListOfGoodGammas.push_back(avect);

                fNNucAnalysed->SetText(Form("Analysed nuclei: %5.0d",AnalysedNuclei));
                fNGrayAnalysed->SetText(Form("Analysed gamma rays: %5.0d",AnalysedGammaRays));

                gSystem->ProcessEvents();
            }
        }
    }

    MapSubwindows();
    MapWindow();
    Layout();
}

void CXGammaSearch::FindInDoubleCoincidence(Bool_t Bash)
{
    //Pre-find nuclei with the right gammas

    std::vector<float> Tab_P;


    std::vector<int> Index_Nuclei;
    std::vector<string> Nuclei_ID;
    std::vector<Float_t> DeltaE;

    Int_t Number_Analysed_Gammas=0;
    Int_t Count=0;

    Int_t NumberOfNuclei = fListOfGoodGammas.size();

    for(int inuc=0 ; inuc<NumberOfNuclei ; inuc++)
    {
        fLevelScheme->Clear();

        vector < GammaTransition > avec = fListOfGoodGammas[inuc];

        if(avec.size()<2)
            continue;

        CXNucleus Nuc(((GammaTransition)avec[0]).NucName.Data());
        fLevelScheme = ImportRAWENSDF(Nuc.GetSymbol());

        /// taille du schema de niveau
        const Int_t Dim_LS = fLevelScheme->GetLinks().GetSize();

        Float_t Flag[Dim_LS];
        for(int ii=0 ; ii<Dim_LS ; ii++)
            Flag[ii] = -1;

        TMatrixD Mat_LS(Dim_LS,Dim_LS); // Raw matrix in TS Space
        TMatrixD Mat_P(Dim_LS,Dim_LS);  // Total probability matrix
        TMatrixD Mat_E(Dim_LS,1);       // Total probability matrix

        Mat_LS = Fill_TS_Matrix(fLevelScheme);
        Mat_P  = fSp->Get_Total_P(Mat_LS,Dim_LS);
        Mat_E  = Get_E_Matrix(fLevelScheme);

        Bool_t Det = fSp->IsDet(Mat_LS,Dim_LS);

        TString TransitionName[Dim_LS];
        Float_t ETrans[Dim_LS];
        Float_t EI[Dim_LS];
        Float_t EF[Dim_LS];

        //get the coordinates of the gammas
        for(Int_t ig=0 ; ig<Dim_LS ; ig++)
        {
            for(int j=0 ; j<avec.size() ; j++)
            {
                GammaTransition GT = avec[j];

                if(GT.EGamma == Mat_E[ig][0])
                {
                    Flag[ig] = Mat_E[ig][0];
                    TransitionName[ig] = Form(" %6.1f keV : %5s (%6.1f keV) --> %5s (%6.1f keV)%s ",GT.EGamma,GT.SpinI.Data(),GT.EI,GT.SpinF.Data(),GT.EF,GT.LifeTime.Data());
                    ETrans[ig] = GT.EGamma;
                    EI[ig] = GT.EI;
                    EF[ig] = GT.EF;
                }
            }
        }

        Bool_t gatetrig[fNGammas];
        Float_t Gates[fNGammas];
        Float_t Width[fNGammas];

        for(int ig=0 ; ig<fNGammas ; ig++) {
            Gates[ig] = fEnergies[ig]->GetNumber();
            Width[ig] = fWidths[ig]->GetNumber();
        }

        for(Int_t j=0 ; j<Dim_LS ; j++)
        {
            for(Int_t m=0 ; m<j ; m++)
            {
                memset(gatetrig,0,sizeof(fNGammas));
                Int_t NGatesTrig=0;

                for(int ig=0 ; ig<fNGammas ; ig++) {
                    if(TMath::Abs(ETrans[j]-Gates[ig])<Width[ig])
                        gatetrig[ig] = true;
                    if(TMath::Abs(ETrans[m]-Gates[ig])<Width[ig])
                        gatetrig[ig] = true;
                }

                for(int ii=0 ; ii<fNGammas ; ii++)
                    NGatesTrig += gatetrig[ii];

                if(NGatesTrig != 2)
                    continue;

                Float_t Prob = TMath::Abs(Mat_P[j][m]);
                if(Flag[j]>0 && Flag[m]>0 && Prob>0.01)
                {
                    //Save the probabilities
                    Tab_P.push_back(Prob);
                    Index_Nuclei.push_back(Count);

                    if(EI[j]>EI[m])
                        DeltaE.push_back(EF[j]-EI[m]);
                    else
                        DeltaE.push_back(EF[m]-EI[j]);

                    TString Name_Output = Form("%s# --> %s# --> %s, Prob = %1.3f", Nuc.GetSymbol().Data(),TransitionName[j].Data(),TransitionName[m].Data(),Prob);

                    if (Det==false)
                        Name_Output+="Warning Matrix Singular";

                    Nuclei_ID.push_back(Name_Output.Data());
                    Count++;
                }
            }
            Number_Analysed_Gammas++;
        }

        gSystem->ProcessEvents();
    }

    std::vector<int> Index_Sorted;
    Index_Sorted = Sort_Index(Index_Nuclei,Tab_P, DeltaE);

    if(!Bash)
    {
        for(Int_t j=(Int_t)Tab_P.size()-1 ; j>=0 ; j--)
        {
            TString Name = Nuclei_ID.at(Index_Sorted[j]);
            TObjArray *arr = Name.Tokenize("#");
            TString NucName = Form("%s : Energy Diff = %1.3f",arr->First()->GetName(),DeltaE[Index_Sorted[j]]);
            TString Transition1 = arr->At(1)->GetName();
            TString Transition2 = arr->At(2)->GetName();

            PrintInListBox(NucName,kInfo);
            PrintInListBox(Transition1,kPrint);
            PrintInListBox(Transition2,kPrint);
        }

        PrintInListBox(Form("Done... %d corresponding nuclei found",((Int_t)Tab_P.size())),kInfo);
    }
}

void CXGammaSearch::FindInTripleCoincidence()
{
    //Pre-find nuclei with the right gammas
    FindGammaRays(true);

    std::vector<float> Tab_P;

    std::vector<int> Index_Nuclei;
    std::vector<string> Nuclei_ID;
    std::vector<Float_t> DeltaE;

    Int_t Number_Analysed_Gammas=0;
    Int_t Count=0;

    Int_t NumberOfNuclei = fListOfGoodGammas.size();

    for(int inuc=0 ; inuc<NumberOfNuclei ; inuc++)
    {
        fLevelScheme->Clear();

        vector < GammaTransition > avec = fListOfGoodGammas[inuc];

        if(avec.size()<2)
            continue;

        CXNucleus Nuc(((GammaTransition)avec[0]).NucName.Data());
        fLevelScheme = ImportRAWENSDF(Nuc.GetSymbol());

        /// taille du schema de niveau
        const Int_t Dim_LS = fLevelScheme->GetLinks().GetSize();

        Float_t Flag[Dim_LS];
        for(int ii=0 ; ii<Dim_LS ; ii++)
            Flag[ii] = -1;

        TMatrixD Mat_LS(Dim_LS,Dim_LS); // Raw matrix in TS Space
        TMatrixD Mat_P(Dim_LS,Dim_LS);  // Total probability matrix
        TMatrixD Mat_E(Dim_LS,1);       // Total probability matrix

        Mat_LS = Fill_TS_Matrix(fLevelScheme);
        Mat_P  = fSp->Get_Total_P(Mat_LS,Dim_LS);
        Mat_E  = Get_E_Matrix(fLevelScheme);

        Bool_t Det = fSp->IsDet(Mat_LS,Dim_LS);

        TString TransitionName[Dim_LS];
        Float_t ETrans[Dim_LS];
        Float_t EI[Dim_LS];
        Float_t EF[Dim_LS];

        //get the coordinates of the gammas
        for(Int_t ig=0 ; ig<Dim_LS ; ig++)
        {
            for(int j=0 ; j<avec.size() ; j++)
            {
                GammaTransition GT = avec[j];

                if(GT.EGamma == Mat_E[ig][0])
                {
                    Flag[ig] = Mat_E[ig][0];
                    TransitionName[ig] = Form(" %6.1f keV : %5s (%6.1f keV) --> %5s (%6.1f keV)%s ",GT.EGamma,GT.SpinI.Data(),GT.EI,GT.SpinF.Data(),GT.EF,GT.LifeTime.Data());
                    ETrans[ig] = GT.EGamma;
                    EI[ig] = GT.EI;
                    EF[ig] = GT.EF;
                }
            }
        }

        Bool_t gatetrig[fNGammas];
        Float_t Gates[fNGammas];
        Float_t Width[fNGammas];

        for(int ig=0 ; ig<fNGammas ; ig++) {
            Gates[ig] = fEnergies[ig]->GetNumber();
            Width[ig] = fWidths[ig]->GetNumber();
        }

        for(Int_t j=0 ; j<Dim_LS ; j++)
        {
            for(Int_t m=0 ; m<j ; m++)
            {
                for(Int_t k=0 ; k<m ; k++)
                {
                    memset(gatetrig,0,sizeof(fNGammas));
                    Int_t NGatesTrig=0;

                    for(int ig=0 ; ig<fNGammas ; ig++) {
                        if(TMath::Abs(ETrans[j]-Gates[ig])<Width[ig])
                            gatetrig[ig] = true;
                        if(TMath::Abs(ETrans[m]-Gates[ig])<Width[ig])
                            gatetrig[ig] = true;
                        if(TMath::Abs(ETrans[k]-Gates[ig])<Width[ig])
                            gatetrig[ig] = true;
                    }

                    for(int ii=0 ; ii<fNGammas ; ii++)
                        NGatesTrig += gatetrig[ii];

                    if(NGatesTrig != 3)
                        continue;

                    Float_t Prob = TMath::Abs(Mat_P[j][m]*Mat_P[j][k]*Mat_P[m][k]);

                    if(Flag[j]>0 && Flag[m]>0 && Flag[k]>0 && Prob>0.01 )
                    {
                        //Save the probabilities
                        Tab_P.push_back(Prob);
                        Index_Nuclei.push_back(Count);

                        if(EI[j]>EI[m] && EI[j]>EI[k]) {
                            if(EI[m]>EI[k])
                                DeltaE.push_back(EF[j]-EI[m] + EF[m]-EI[k]);
                            else
                                DeltaE.push_back(EF[j]-EI[k] + EF[k]-EI[m]);
                        }
                        else if(EI[m]>EI[j] && EI[m]>EI[k]) {
                            if(EI[j]>EI[k])
                                DeltaE.push_back(EF[m]-EI[j] + EF[j]-EI[k]);
                            else
                                DeltaE.push_back(EF[m]-EI[k] + EF[k]-EI[j]);
                        }
                        else if(EI[k]>EI[j] && EI[k]>EI[m]) {
                            if(EI[j]>EI[m])
                                DeltaE.push_back(EF[k]-EI[j] + EF[j]-EI[m]);
                            else
                                DeltaE.push_back(EF[k]-EI[m] + EF[m]-EI[j]);
                        }

                        TString Name_Output = Form("%s# --> %s# --> %s# --> %s, Prob = %1.3f", Nuc.GetSymbol().Data(),TransitionName[j].Data(),TransitionName[m].Data(),TransitionName[k].Data(),Prob);

                        if (Det==false)
                            Name_Output+="Warning Matrix Singular";

                        Nuclei_ID.push_back(Name_Output.Data());
                        Count++;
                    }
                }
                Number_Analysed_Gammas++;
            }

            gSystem->ProcessEvents();
        }
    }

    std::vector<int> Index_Sorted;
    Index_Sorted = Sort_Index(Index_Nuclei,Tab_P, DeltaE);


    for(Int_t j=(Int_t)Tab_P.size()-1 ; j>=0 ; j--)
    {
        TString Name = Nuclei_ID.at(Index_Sorted[j]);
        TObjArray *arr = Name.Tokenize("#");
        TString NucName = Form("%s : Energy Diff = %1.3f",arr->First()->GetName(),DeltaE[Index_Sorted[j]]);
        TString Transition1 = arr->At(1)->GetName();
        TString Transition2 = arr->At(2)->GetName();
        TString Transition3 = arr->At(3)->GetName();

        PrintInListBox(NucName,kInfo);
        PrintInListBox(Transition1,kPrint);
        PrintInListBox(Transition2,kPrint);
        PrintInListBox(Transition3,kPrint);
    }

    PrintInListBox(Form("Done... %d corresponding nuclei found",((Int_t)Tab_P.size())),kInfo);
}

std::vector<int> CXGammaSearch::Sort_Index(std::vector<int> Index_Nuclei, std::vector<float> Tab_P, std::vector<float> Tab_DeltaE)
{
    std::vector<int> vecX;

    for(unsigned int i=1 ; i<Index_Nuclei.size() ; i++)
        for(unsigned int j=0 ; j<Index_Nuclei.size()-1 ; j++)
            if(Tab_P[j]>Tab_P[j+1]) {
                swap(Index_Nuclei[j],Index_Nuclei[j+1]);
                swap(Tab_P[j],Tab_P[j+1]);
                swap(Tab_DeltaE[j],Tab_DeltaE[j+1]);
            }

    for(unsigned int i=1 ; i<Index_Nuclei.size() ; i++)
        for(unsigned int j=0 ; j<Index_Nuclei.size()-1 ; j++) {
            if(TMath::Abs(Tab_P[j]-Tab_P[j+1])<1e-4) {
                if(Tab_DeltaE[j]<Tab_DeltaE[j+1]) {
                    swap(Index_Nuclei[j],Index_Nuclei[j+1]);
                    swap(Tab_P[j],Tab_P[j+1]);
                    swap(Tab_DeltaE[j],Tab_DeltaE[j+1]);
                }
            }
        }


    for(uint i=0 ; i<Tab_P.size() ; i++)
        vecX.push_back(Index_Nuclei[i]);

    return vecX ;
}

void CXGammaSearch::PrintInListBox(TString mess, Int_t Type)
{

#if (OS_TYPE == OS_LINUX)
    const TGFont *ufont;         // will reflect user font changes
    ufont = gClient->GetFont("-adobe-courier-medium-r-*-*-14-*-*-*-*-*-iso8859-1");
   // ufont = gClient->GetFont("-adobe-times-medium-r-*-*-12-*-*-*-*-*-iso8859-1");
    if (!ufont)
       ufont = fClient->GetResourcePool()->GetDefaultFont();

    TGGC   *uGC;           // will reflect user GC changes
    // graphics context changes
    GCValues_t val;
    val.fMask = kGCFont;
    val.fFont = ufont->GetFontHandle();
    uGC = gClient->GetGC(&val, kTRUE);

    TGTextLBEntry *entry = new TGTextLBEntry(fResultsBox->GetContainer(), new TGString(mess), fResultsBox->GetNumberOfEntries()+1, uGC->GetGC(), ufont->GetFontStruct());
#else
    TGTextLBEntry *entry = new TGTextLBEntry(fResultsBox->GetContainer(), new TGString(mess), fResultsBox->GetNumberOfEntries()+1);
#endif

    if(Type == kError)
        entry->SetBackgroundColor((Pixel_t)0xff0000);
    else if(Type == kInfo)
        entry->SetBackgroundColor((Pixel_t)0x87a7d2);
    else if(Type == kWarning)
        entry->SetBackgroundColor((Pixel_t)0xdfdf44);
    else if(Type == kPrint)
        entry->SetBackgroundColor((Pixel_t)0x90f269);

    fResultsBox->AddEntry((TGLBEntry *)entry, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX));
    fResultsBox->Layout();
}

void CXGammaSearch::HandleButtons()
{
    if(((TObject*)gTQSender)->InheritsFrom("TGNumberEntry"))
    {
        if(gTQSender == fNRange[0])
            fARange[0]->SetNumber(fZRange[0]->GetNumber()+fNRange[0]->GetNumber());
        if(gTQSender == fNRange[1])
            fARange[1]->SetNumber(fZRange[1]->GetNumber()+fNRange[1]->GetNumber());
        if(gTQSender == fARange[0])
            fNRange[0]->SetNumber(fARange[0]->GetNumber()-fZRange[0]->GetNumber());
        if(gTQSender == fARange[1])
            fNRange[1]->SetNumber(fARange[1]->GetNumber()-fZRange[1]->GetNumber());
        if(gTQSender == fZRange[0])
            fARange[0]->SetNumber(fZRange[0]->GetNumber()+fNRange[0]->GetNumber());
        if(gTQSender == fZRange[1])
            fARange[1]->SetNumber(fZRange[1]->GetNumber()+fNRange[1]->GetNumber());
    }

    if(((TObject*)gTQSender)->InheritsFrom("TGRadioButton"))
    {
        TGButton *btn = (TGButton *) gTQSender;
        Int_t id = btn->WidgetId();

        if(id==M_With_Cascade)
        {
            fNoCoincMode->SetState(kButtonUp);
            fCoincMode->SetState(kButtonDown);
            fCurrentMode=M_With_Cascade;
        }
        if(id ==M_WithOut_Cascade)
        {
            fCoincMode->SetState(kButtonUp);
            fNoCoincMode->SetState(kButtonDown);
            fCurrentMode=M_WithOut_Cascade;
        }
    }


    fNGammas = 0;
    for(int i=0 ; i<3 ; i++)
    {
        fEnergies[i]->SetState(fCheckGammas[i]->GetState());
        fWidths[i]->SetState(fCheckGammas[i]->GetState());

        if(fCheckGammas[i]->GetState() == kButtonDown)
            fNGammas++;
    }
}


LevelScheme *CXGammaSearch::ImportRAWENSDF(const char *NucSymbol)
{
    CXNucleus nuc(NucSymbol);
    TString ZMaj = ((TString)nuc.GetSymbol());
    ZMaj.ToUpper();
    ZMaj.Append("NoVerbose");

    TString FileName = fDataBaseFolder + "/" + ((TString)nuc.GetSymbol().Data()) + ".ens";

    if(gSystem->IsFileInIncludePath(FileName))
    {
        LevelScheme *lev  =  new LevelScheme;
        lev->Import(FileName.Data(),ZMaj.Data());

        if(lev->GetLinks().GetSize() == 0)
            return nullptr;
        else
            return lev;
    }
    else
        return nullptr;
}


TMatrixD CXGammaSearch::Fill_TS_Matrix(LevelScheme *lev)
{
    Int_t Size=lev->GetLinks().GetSize();

    TMatrixD Fill_M(Size,Size);

    for(Int_t i=0;i<(Int_t)Size;i++)
    {

        const GammaLink *Link=(GammaLink*)lev->GetLinks().At(i);
        Measure<Float_t> E = ((GammaLink*)Link)->GetEnergy();
        NuclearLevel *FL = (NuclearLevel*)Link->GetFL();
        Measure<Float_t> FLE = ((NuclearLevel*)FL)->GetEnergy();

        Float_t Br=0;

        for(Int_t j=0;j<(Int_t)Size;j++)//Get The branching ratio
        {

            const GammaLink *Link2=(GammaLink*)lev->GetLinks().At(j);
            Measure<Float_t> E2 = ((GammaLink*)Link2)->GetEnergy();
            NuclearLevel *IL = (NuclearLevel*)Link2->GetIL();
            Measure<Float_t> ILE = ((NuclearLevel*)IL)->GetEnergy();
            Measure<Float_t> S = ((GammaLink*)Link2)->GetStrength();

            Float_t toto = S.GetValue();
            if(toto<1e-6) toto = 1e-6;

            if(FLE.GetValue()==ILE.GetValue())
            {
                Br+=toto;
            }
        }

        for(Int_t j=0;j<(Int_t)Size;j++)//Fill the matrix
        {

            const GammaLink *Link2=(GammaLink*)lev->GetLinks().At(j);
            Measure<Float_t> E2 = ((GammaLink*)Link2)->GetEnergy();
            NuclearLevel *IL = (NuclearLevel*)Link2->GetIL();
            Measure<Float_t> ILE = ((NuclearLevel*)IL)->GetEnergy();
            Measure<Float_t> S = ((GammaLink*)Link2)->GetStrength();

            Float_t toto = S.GetValue();
            if(toto<1e-6) toto = 1e-6;

            if(FLE.GetValue()==ILE.GetValue())
            {
                if(!TMath::IsNaN(toto/Br))Fill_M[i][j]=toto/Br;
            }
        }
    }

    return Fill_M;
}

TMatrixD CXGammaSearch::Get_E_Matrix(LevelScheme *lev)
{

    const Int_t Size=lev->GetLinks().GetSize();

    TMatrixD Fill_E(Size,1);

    for(Int_t j=0;j<(Int_t)Size;j++)//Fill the matrix with the Energy
    {
        const GammaLink *Link2=(GammaLink*)lev->GetLinks().At(j);
        Measure<Float_t> E = ((GammaLink*)Link2)->GetEnergy();

        Fill_E[j][0]=E.GetValue();
    }

    return Fill_E;
}


TString CXGammaSearch::GetSpinFromLev(NuclearLevel *nuclev)
{
    Float_t spin = nuclev->GetSpin().Get();

    TString  PL = "+";
    if(nuclev->GetParity().IsParity(Gw::Parity::kMinus))
        PL = "-";

    TString SL = Form("%.1f",spin);

    TString Lev;

    if(SL.EndsWith(".5"))
        Lev = Form("%d/2%s",TMath::Nint(2*SL.Atof()),PL.Data());
    else
        Lev = Form("%d%s",SL.Atoi(),PL.Data());

    if(nuclev->GetSpin().IsData(InfoData::kUnknown))
        Lev="?";
    else if(nuclev->GetSpin().IsData(InfoData::kTentative))
        Lev.Append(")").Prepend("(");

    return Lev;
}


ClassImp(CXGammaSearch);