#include "TICPhysics.h"
#include "TTimeData.h"
#include <TFPMWPhysics.h>
#include <TChain.h>
#include <TCutG.h>
#include <TF1.h>
#include <TFile.h>
#include <TH2.h>
#include <TSpectrum.h>
#include <TSpline.h>
#include <fstream>
using namespace std;


const int CutNumberLoader();
const int Ncuts = CutNumberLoader();
const int NSegment=11;
const int Nrun=8;
int RunNumber[Nrun]={245 , 246 ,247 ,248 ,249 ,250 ,251 ,252};
const char *CutPath = "../Step4/Output/CutAutoZ.root" ;

void MakeSpline();
void ApplySpline_perTCutG();
void ApplySpline();
vector<double> GetMinMaxX(TCutG* cut);
vector<double> GetMinMaxY(TCutG* cut) ;
TF1* Linear(double x1, double y1, double x2, double y2, const char* name = "f");
///////////////////////////////////////////////////////////////////////////
void SplineChio()
{
    //===========================================================================================================
    //                                      Loading Var
    //===========================================================================================================

    //**********************Cut**************************************
    TFile *fcut=new TFile(CutPath,"open");

    vector<TCutG*> cutZ(Ncuts);
    for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));

    //**********************DE_E****************************************************
    TChain *chain=new TChain("PhysicsTree");

    for(int i=0;i<Nrun;i++) 
    {
        chain->Add(Form("../../../../root/analysis/Run%d.root",RunNumber[i]));
        cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
    }

    TICPhysics* IC = new TICPhysics() ;
    TTimeData *Time = new TTimeData() ;


    chain->SetBranchStatus("IC", true);
    chain->SetBranchAddress("IC", &IC);

    double FF_IC_X;
    chain->SetBranchStatus("FF_IC_X", true);
    chain->SetBranchAddress("FF_IC_X", &FF_IC_X);


    // ****************** Load Time  **********************************

    chain->SetBranchStatus("Time", true);
    chain->SetBranchAddress("Time", &Time);

    //*****************************************************************

    TH2F *hChioDE_E_all = new TH2F("hChioDE_E","hChioDE_E",1000,0,36000,500,1000,26000);
    vector<TH2F*> hChioDE_E(Ncuts);
    for(int i=0;i<Ncuts;i++) hChioDE_E[i]=new TH2F(Form("hChioDE_E_%d",i+1),Form("hChioDE_E_%d",i+1),2000,0,36000,1250,1000,26000);

    auto start = std::chrono::high_resolution_clock::now();

    int Nentries=chain->GetEntries();
    //int Nentries=1000000;

    for(int e=0;e<Nentries;++e)
    {
        if (e % 100000  == 0 && e > 0 ) {
            auto now = std::chrono::high_resolution_clock::now();
            std::chrono::duration<double> elapsed = now - start;
            double avgTimePerIteration = elapsed.count() / e;
            double timeLeft = avgTimePerIteration * (Nentries - e);

            std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
        }

        chain->GetEntry(e);

        if(FF_IC_X>-530 ){
            hChioDE_E_all->Fill(IC->Eres,IC->DE);
        }

        //for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(Chio_E,FF_DE)) {hChioDE_E[i]->Fill(Chio_E,FF_DE); break;}
        for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(IC->Eres,IC->DE)) {hChioDE_E[i]->Fill(IC->Eres,IC->DE); break;}
    }//end loop event


    //**********************************Out***************************
    TFile *outFile=new TFile("histo/SingleZ_ChioDE_E.root","recreate");
    for(int i=0;i<Ncuts;i++) hChioDE_E[i]->Write();
    outFile->Close();

    TCanvas * can = new TCanvas(Form("ChioEbis_vs_ChioDE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            Form("ChioEbis_vs_ChiodE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            0,0,2000,1000);
    can->cd();
    hChioDE_E_all->Draw("colz");
    for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
}


///////////////////////////////////////////////////////////////////////////
void MakeSpline(){
    // *************************Cut***************************
    TFile *fcut=new TFile(CutPath, "open");
    vector<TCutG*> cutZ(Ncuts);
    // ********************** Load prev histo****************************
    TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root");

    int Z59POS = 28;
    vector<TSpline3*> gspline(Ncuts);  
    vector<TH2F*> h2 (Ncuts);
    vector<TH1F*> hProfile(Ncuts);
    vector<TH1F*> pfx(Ncuts); // extended hProfile 
    TFile* fspline=new TFile("Output/spline_Chio_2024.root","recreate");

    TCanvas * canfit = new TCanvas("canfit","canfit",0,0,2000,1500);

    for(int i=0;i<Ncuts;i++)
    {
        // Get the 2D plot in TCutG
        cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
        h2[i]=(TH2F*)inFile->Get(Form("hChioDE_E_%d",i+1));

        // Create the TProfile
        hProfile[i]=(TH1F*)h2[i]->ProfileX();
        hProfile[i]->SetLineWidth(2);
        hProfile[i]->SetDirectory(0);
        canfit->cd();
        if (i==0) {
            hProfile[i]->GetYaxis()->SetRangeUser(5000,26000);
            hProfile[i]->Draw();
        }
        else hProfile[i]->Draw("same");

        // Get the range of the TProfile
        int FirstBin, LastBin;
        double parpol3[4];
        for(int bin=1; bin<hProfile[i]->GetNbinsX(); bin++){
            FirstBin = bin;
            if (hProfile[i]->GetBinContent(bin)>6000) break;
        }
        double parpol1[2];
        for(int bin=hProfile[i]->GetNbinsX(); bin>1 ; bin--){
            LastBin = bin;
            if (hProfile[i]->GetBinContent(bin)>6000) break;
        }

        // Create the extended TProfile
        pfx[i] = new TH1F(Form("pfx_%02d",i+1),Form("pfx_%02d",i+1),hProfile[i]->GetNbinsX(),hProfile[i]->GetBinLowEdge(1),
                hProfile[i]->GetBinLowEdge(hProfile[i]->GetNbinsX())+hProfile[i]->GetBinWidth(hProfile[i]->GetNbinsX()));
        pfx[i]->SetLineColor(8);
        pfx[i]->SetDirectory(0);

        // find the function to extend the TProfile on the lower range
        TF1 * fitpol3 = new TF1(Form("FitPol3_pfx_%02d",i+1),"pol3",hProfile[i]->GetBinLowEdge(FirstBin),hProfile[i]->GetBinLowEdge(FirstBin+200));
        hProfile[i]->Fit(fitpol3,"R");
        fitpol3->GetParameters(parpol3);

        // find the function to extend the TProfile on the higher range
        TF1 * fitpol1 = new TF1(Form("FitPol1_pfx_%02d",i+1),"pol1",hProfile[i]->GetBinLowEdge(LastBin)-10000,hProfile[i]->GetBinLowEdge(LastBin));
        hProfile[i]->Fit(fitpol1,"R");
        fitpol1->GetParameters(parpol1);

        // fill the extended TProfile
        float newval,lastval,lasterr;
        for(int bin=1; bin<=FirstBin; bin++){
            newval=0;
            if(i < (Ncuts-4) && i!=0 && i!=1)
                for(int par=0; par<4; par++) newval += parpol3[par]*pow(hProfile[i]->GetBinCenter(bin),par);
            else if (i >= (Ncuts-4))
                newval = hProfile[i]->GetBinContent(FirstBin);
            else{
                newval=0;
                for(int par=0; par<2; par++) newval += parpol1[par]*pow(hProfile[i]->GetBinCenter(bin),par);
                pfx[i]->SetBinContent(bin,newval);
            }
            pfx[i]->SetBinContent(bin,newval);
        }
        for(int bin=FirstBin; bin<=LastBin; bin++){
            newval = hProfile[i]->GetBinContent(bin);
            if (newval!=0){
                pfx[i]->SetBinContent(bin,newval);
                pfx[i]->SetBinError(bin,hProfile[i]->GetBinError(bin));
                lastval = newval;
                lasterr = hProfile[i]->GetBinError(bin);
            }
            else{
                pfx[i]->SetBinContent(bin,lastval);
                pfx[i]->SetBinError(bin,lasterr);
            }
        }
        for(int bin=LastBin; bin<=hProfile[i]->GetNbinsX(); bin++){
            newval=0;
            for(int par=0; par<2; par++) newval += parpol1[par]*pow(hProfile[i]->GetBinCenter(bin),par);
            pfx[i]->SetBinContent(bin,newval);
        }
        pfx[i]->Draw("same");

        gspline[i]=new TSpline3(pfx[i]);
        gspline[i]->SetName(Form("fspline_%d",i+1));
        fspline->cd();
        gspline[i]->Write();

    }

    fspline->Close();

    TFile* CalOut=new TFile("../../../../Calibration/VAMOS/CHIO/Z_spline.root","recreate");
    for (TSpline3 * i : gspline){
        i->Write();
    } 
    CalOut->Close();

    TGraph * gr = new TGraph();
    gr->SetName("DEvsZ");

    TCanvas * can = new TCanvas("ChioEcorr_vs_Ebis_for_TSplines","ChioEcorr_vs_Ebis_for_TSplines",0,0,2500,2000);
    can->cd();
    ofstream ofile_par;
    string filename_par = "Output/mean_spline_par.dat";
    ofile_par.open(filename_par.c_str());
    cout << "Float_t FF_DEcorr[" << Ncuts << "] = { " << endl; 
    ofile_par << "Float_t FF_DEcorr[" << Ncuts << "] = { " << endl;



    vector<double> XEdgesMin = GetMinMaxX(cutZ.at(0));
    vector<double> XEdgesMax = GetMinMaxX(cutZ.at(Ncuts-1));
    vector<double> YEdgesMin = GetMinMaxY(cutZ.at(0));
    vector<double> YEdgesMax = GetMinMaxY(cutZ.at(Ncuts-1));

    double x1 = (XEdgesMin.at(0) + XEdgesMin.at(1)) /2;
    double x2 = (XEdgesMax.at(0) + XEdgesMax.at(1)) /2;
    double y1 = (YEdgesMin.at(0) + YEdgesMin.at(1)) /2;
    double y2 = (YEdgesMax.at(0) + YEdgesMax.at(1)) /2;

    TF1 *MeanPos = Linear(y1,x1,y2,x2,"MeanPos");

    for(int i=0;i<Ncuts;i++){
        can->cd();
        if(i==0) h2[i]->Draw("col");
        else     h2[i]->Draw("col,same");
        cout << h2[i]->GetMean(2) << ", " << endl;
        ofile_par << h2[i]->GetMean(2) << ", " << endl;
        //gr->SetPoint(i,sqrt(h2[i]->GetMean(2)),i+31);
        gr->SetPoint(i,gspline[i]->Eval(MeanPos->Eval((YEdgesMin.at(0)+YEdgesMin.at(1))/2)),i+Z59POS);
        //hProfile[i]->Draw("same");
        pfx[i]->Draw("same");
        gspline[i]->SetLineColor(kBlue);
        gspline[i]->SetLineWidth(3);
        gspline[i]->Draw("lsame");
    }

    TCanvas * cangr = new TCanvas("ToGetRoughCal","ToGetRoughCal",200,0,1200,1000);
    cangr->cd();
    gr->SetMarkerStyle(20);
    gr->Draw("AP");
    gr->Fit("pol5");
}


///////////////////////////////////////////////////////////////////////////
void ApplySpline_perTCutG()
{

    // *************************Cut***************************
    TFile *fcut=new TFile(CutPath);
    vector<TCutG*> cutZ(Ncuts);

    // ******************* TSpline DE*************************
    TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
    vector<TSpline3*> gspline(Ncuts);  

    for(int i=0;i<Ncuts;i++){ 
        cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
        gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
    }

    //**********************DE_E****************************************************
    TChain *chain=new TChain("PhysicsTree");


    for(int i=0;i<Nrun;i++) 
    {
        chain->Add(Form("../../../../root/analysis/Run%d.root",RunNumber[i]));
        cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
    }    
    TICPhysics* IC = new TICPhysics() ;

    chain->SetBranchStatus("IC", true);
    chain->SetBranchAddress("IC", &IC); 

    // variable after applying TSpline MUST FIND middles of each spline
    double FF_DEcorr;
    vector<double> FF_DEcorr0 = { 
        7375.73, 
        7982.83, 
        8419.21, 
        8801.52, 
        9227.89, 
        9667.77, 
        10163.1, 
        10680.5, 
        11222.2, 
        11758.3, 
        12287.6, 
        12802.2, 
        13324.2, 
        13828.4, 
        14325.2, 
        14773.7, 
        15207.3, 
        15674.6, 
        16098.9, 
        16548.5, 
        16926.9, 
        17384.1, 
        17778.3, 
        18196.6, 
        18597.6, 
        19043.9, 
        19437.8, 
        19889.4,}; 
    // 20301.1, 
    //20709.7, 
    //21112.6, 
    //21544, 
    //21953.4, 
    //22361.8};


TH2F *hChioDE_E_all  = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000);
TH2F *hChioDE_E      = new TH2F("hChioDE_E","hChioDE_E",1000,1000,41000,500,1000,26000);
TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);

int Nentries=chain->GetEntries();
//int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();

for(int e=0;e<Nentries;++e){
    if (e % 100000  == 0 && e > 0 ) {
        auto now = std::chrono::high_resolution_clock::now();
        std::chrono::duration<double> elapsed = now - start;
        double avgTimePerIteration = elapsed.count() / e;
        double timeLeft = avgTimePerIteration * (Nentries - e);

        std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
    }

    chain->GetEntry(e);

    hChioDE_E_all->Fill(IC->Eres,IC->DE); 
    for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(IC->Eres,IC->DE)) {
        hChioDE_E->Fill(IC->Eres,IC->DE); 
        FF_DEcorr = FF_DEcorr0[i] * IC->DE / gspline[i]->Eval(IC->Eres); 
        hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr); 
        break;
    }
}

TCanvas * can = new TCanvas(Form("ChioEbis_vs_ChioDE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
        Form("ChioEbis_vs_ChiodE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
        0,0,2000,1000);
can->Divide(1,3);
can->cd(1); gPad-> SetLogz(); hChioDE_E_all->Draw("colz");
can->cd(2); gPad-> SetLogz(); hChioDE_E->Draw("colz");
can->cd(3); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz");

}
void ApplySplineUnique()
{
    // *************************Cut***************************
    TFile *fcut=new TFile(CutPath, "open");
    vector<TCutG*> cutZ(Ncuts);

    // ******************* TSpline DE*************************
    TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
    vector<TSpline3*> gspline(Ncuts) ;  

    for(int i=0;i<Ncuts;i++){ 
        cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
        gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
    }

    //**********************DE_E****************************************************
    TChain *chain=new TChain("PhysicsTree");
    for(int i=0;i<Nrun;i++) 
    {
        chain->Add(Form("../../../../root/analysis/Run%d.root",RunNumber[i]));
        cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
    }   

    TICPhysics* IC = new TICPhysics() ;
    chain->SetBranchStatus("IC", true);
    chain->SetBranchAddress("IC", &IC); 


    // variable after applying TSpline
    double FF_DEcorr=0;
    double Zrough=0;
    double Zabitbetter=0;

    vector<double> FF_DEcorr0(Ncuts);
    for(int index=0; index<Ncuts; index++){  
        FF_DEcorr0[index] = gspline[index]->Eval(8500);
    }

    // histos
    TH2F *hChioDE_E_all  = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000);
    TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);
    TH1F *hChioDE_corr   = new TH1F("hChioDE_corr","hChioDE_corr",2500,1000,26000);
    TH2F *hChioZ_E_rough = new TH2F("hChioZ_E_rough","hChioZ_E_rough",1000,1000,41000,500,25,65); 
    TH1F *hChioZ_rough = new TH1F("hChioZ_rough","hChioZ_rough",2500,25,65); 
    TH1F *hChioZ_abitbetter = new TH1F("hChioZ_abitbetter","hChioZ_abitbetter",2500,25,65); 

    //int Nentries=1e6;
    int Nentries = chain->GetEntries();
    auto start = std::chrono::high_resolution_clock::now();

    double FF_Eres_prev = 0;
    for(int e=0;e<Nentries;e++)
    {
        if (e % 100000  == 0 && e > 0 ) {
            auto now = std::chrono::high_resolution_clock::now();
            std::chrono::duration<double> elapsed = now - start;
            double avgTimePerIteration = elapsed.count() / e;
            double timeLeft = avgTimePerIteration * (Nentries - e);

            std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
        }

        FF_DEcorr = -100;
        chain->GetEntry(e);

        if(IC->Eres==FF_Eres_prev) continue;
        FF_Eres_prev = IC->Eres;
        hChioDE_E_all->Fill(IC->Eres,IC->DE);

        float Eval_DEspline, DEspline0;
        int   index=0;

        Eval_DEspline = gspline[24]->Eval(IC->Eres);
        DEspline0 = FF_DEcorr0[24];

        FF_DEcorr = DEspline0 * IC->DE / Eval_DEspline ;
        hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr);
        //if(FF_DEcorr>15120 && FF_DEcorr<15130) {
        //  cout << e << " " << index << " " << IC->Eres << " " << FF_DEcorr << endl;
        //}

        // should be pol2 !!!
        //Zrough = -110.165 + 3.34475*sqrt(FF_DEcorr) - 0.0271123*FF_DEcorr + 8.60752e-05 * pow(sqrt(FF_DEcorr),3);
        Zrough = 16.8521 + 0.0017328*FF_DEcorr + 1.70774e-8*pow(FF_DEcorr,2);
        Zabitbetter = -127.117 + 3.83463*sqrt(FF_DEcorr) - 0.0317448 *FF_DEcorr + 0.000100428 * pow(sqrt(FF_DEcorr),3);
        hChioZ_E_rough->Fill(IC->Eres,Zrough);

        if(3000<IC->Eres && IC->Eres<25000){ 
            hChioDE_corr->Fill(FF_DEcorr);
            hChioZ_rough->Fill(Zrough);
            hChioZ_abitbetter->Fill(Zabitbetter);
        }
    }

    TCanvas * can = new TCanvas(Form("ChioEbis_vs_ChioDE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            Form("ChioEbis_vs_ChiodE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            0,0,2000,2000);
    can->Divide(1,3);
    can->cd(1); gPad-> SetLogz(); hChioDE_E_all->Draw("colz");
    for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
    can->cd(2); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz");
    can->cd(3); gPad-> SetLogz(); hChioZ_E_rough->Draw("colz");

    TCanvas * can1D = new TCanvas(Form("DE_and_Z_ifE_8000to25000_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            Form("DE_and_Z_ifE_8000to25000_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            0, 0, 2000,2000);
    can1D->Divide(1,3);
    can1D->cd(1);  hChioDE_corr->Draw(); 
    can1D->cd(2);  hChioZ_rough->Draw();
    can1D->cd(3);  hChioZ_abitbetter->Draw();

    TFile * fsave = new TFile(Form("Output/DEvsE_corr_run%04d_%04d.root",RunNumber[0],RunNumber[Nrun-1]),"recreate");  
    can->Write();
    hChioDE_E_all->Write();
    hChioDE_E_corr->Write();
    hChioZ_E_rough->Write();
    can1D->Write();
    hChioDE_corr->Write();
    hChioZ_rough->Write();
    hChioZ_abitbetter->Write();
    fsave->Close();
}


///////////////////////////////////////////////////////////////////////////
void ApplySpline()
{
    // *************************Cut***************************
    TFile *fcut=new TFile(CutPath, "open");
    vector<TCutG*> cutZ(Ncuts);

    // ******************* TSpline DE*************************
    TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
    vector<TSpline3*> gspline(Ncuts) ;  

    for(int i=0;i<Ncuts;i++){ 
        cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
        gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
    }

    //**********************DE_E****************************************************
    TChain *chain=new TChain("PhysicsTree");
    for(int i=0;i<Nrun;i++) 
    {
        chain->Add(Form("../../../../root/analysis/Run%d.root",RunNumber[i]));
        cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
    }   

    TICPhysics* IC = new TICPhysics() ;
    chain->SetBranchStatus("IC", true);
    chain->SetBranchAddress("IC", &IC); 

    double FF_IC_X;
    chain->SetBranchStatus("FF_IC_X", true);
    chain->SetBranchAddress("FF_IC_X", &FF_IC_X);


    // variable after applying TSpline
    double FF_DEcorr=0;
    double Zrough=0;
    double Zabitbetter=0;
    vector<double> FF_DEcorr0(Ncuts);

    vector<double> XEdgesMin = GetMinMaxX(cutZ.at(0));
    vector<double> XEdgesMax = GetMinMaxX(cutZ.at(Ncuts-1));
    vector<double> YEdgesMin = GetMinMaxY(cutZ.at(0));
    vector<double> YEdgesMax = GetMinMaxY(cutZ.at(Ncuts-1));

    double x1 = (XEdgesMin.at(0) + XEdgesMin.at(1)) /2;
    double x2 = (XEdgesMax.at(0) + XEdgesMax.at(1)) /2;
    double y1 = (YEdgesMin.at(0) + YEdgesMin.at(1)) /2;
    double y2 = (YEdgesMax.at(0) + YEdgesMax.at(1)) /2;

    TF1 *MeanPos = Linear(y1,x1,y2,x2,"MeanPos");

    for(int index=0; index<Ncuts; index++){  
        vector<double> YEdgesMinTemp = GetMinMaxY(cutZ.at(index));
        cout << " E " <<MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2) << endl;
        cout << "DE " << gspline[index]->Eval(MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2)) << endl;
        FF_DEcorr0[index] = gspline[index]->Eval(MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2));
    }
    // histos
    TH2F *hChioDE_E_all  = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000);
    TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);
    TH1F *hChioDE_corr   = new TH1F("hChioDE_corr","hChioDE_corr",1000,1000,26000);
    TH2F *hChioZ_E_rough = new TH2F("hChioZ_E_rough","hChioZ_E_rough",1000,1000,41000,500,25,65); 
    TH1F *hChioZ_rough = new TH1F("hChioZ_rough","hChioZ_rough",2500,25,65); 
    TH1F *hChioZ_abitbetter = new TH1F("hChioZ_abitbetter","hChioZ_abitbetter",2500,25,65); 

    std::ifstream file("Output/Calibration.txt");
    bool Loaded = true ;

    if (!file) {
        std::cerr << "Error opening file.\n";
        Loaded = false;
    }
    Double_t p[6]; // Assuming 5 values
    if (Loaded == true){
        for (int i = 0; i < 6; ++i) {
            file >> p[i];
        }
    }
    int Nentries=2e6;
    //int Nentries = chain->GetEntries();
    auto start = std::chrono::high_resolution_clock::now();

    double FF_Eres_prev = 0;
    for(int e=0;e<Nentries;e++)
    {
        if (e % 100000  == 0 && e > 0 ) {
            auto now = std::chrono::high_resolution_clock::now();
            std::chrono::duration<double> elapsed = now - start;
            double avgTimePerIteration = elapsed.count() / e;
            double timeLeft = avgTimePerIteration * (Nentries - e);

            std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
        }

        FF_DEcorr = -100;
        chain->GetEntry(e);

        if (IC->DE<3000) continue;
        if(IC->Eres==FF_Eres_prev) continue;
        if (FF_IC_X<-530) continue;
        FF_Eres_prev = IC->Eres;

        hChioDE_E_all->Fill(IC->Eres,IC->DE);

        float Eval_DEspline, DEspline0;
        int   index=0;
        for(int i=0; i<Ncuts; i++){
            Eval_DEspline = gspline[i]->Eval(IC->Eres); 
            if(IC->DE<Eval_DEspline) break;
            index = i;
        }

        Eval_DEspline = gspline[index]->Eval(IC->Eres);
        DEspline0 = FF_DEcorr0[index];
        float dmin, dsup;
        if( index < (Ncuts-1) && IC->DE > gspline[0]->Eval(IC->Eres)){
            dmin = abs(IC->DE - gspline[index]->Eval(IC->Eres));
            dsup = abs(gspline[index+1]->Eval(IC->Eres) - IC->DE);
            if(dmin<0) cout << "negative value of dmin = " << dmin << ", for index = " << index << endl;
            if(dsup<0) cout << "negative value of dsup = " << dsup << ", for index = " << index << endl;
            //if(dsup<dmin) {
            //  Eval_DEspline = gspline[index+1]->Eval(IC->Eres) ;
            //  DEspline0 = FF_DEcorr0[index] ;
            //}
            Eval_DEspline = dsup*gspline[index]->Eval(IC->Eres)/(dmin+dsup) + dmin*gspline[index+1]->Eval(IC->Eres)/(dmin+dsup) ;
            DEspline0 = dsup*FF_DEcorr0[index]/(dmin+dsup) + dmin*FF_DEcorr0[index+1]/(dmin+dsup) ;

            FF_DEcorr = DEspline0 * IC->DE / Eval_DEspline ;
            hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr);
            //if(FF_DEcorr>15120 && FF_DEcorr<15130) {
            //  cout << e << " " << index << " " << IC->Eres << " " << FF_DEcorr << endl;
            //}
        }

        // should be pol2 !!!
        //Zrough = -110.165 + 3.34475*sqrt(FF_DEcorr) - 0.0271123*FF_DEcorr + 8.60752e-05 * pow(sqrt(FF_DEcorr),3);

        if (Loaded == true){


            Zrough = p[0] + p[1]*FF_DEcorr +p[2]*pow(FF_DEcorr,2) + p[3]*pow(FF_DEcorr,3) +p[4]*pow(FF_DEcorr,4) + p[5]*pow(FF_DEcorr,5);
        }

        else {
            Zrough = -57.068 + 0.0275327*FF_DEcorr - 3.52303e-06*pow(FF_DEcorr,2) + 2.36906e-10*pow(FF_DEcorr,3) - 7.74553e-15*pow(FF_DEcorr,4) + 9.90229e-20*pow(FF_DEcorr,5);
        }
        //out << Zrough << endl;
        Zabitbetter = -127.117 + 3.83463*sqrt(FF_DEcorr) - 0.0317448 *FF_DEcorr + 0.000100428 * pow(sqrt(FF_DEcorr),3);
        hChioZ_E_rough->Fill(IC->Eres,Zrough);

        if(3000<IC->Eres && IC->Eres<25000){ 
            hChioDE_corr->Fill(FF_DEcorr);
            hChioZ_rough->Fill(Zrough);
            hChioZ_abitbetter->Fill(Zabitbetter);
        }
    }

    // ********************** Fit pol5 of each Z ***********************************
    //
    // PeakFinder
    TSpectrum *PeakFinder = new TSpectrum(100); 
    int nPeaks = PeakFinder->Search(hChioDE_corr,2,"",0.040); 
    Double_t * x = PeakFinder->GetPositionX();
    std::sort(x,x+nPeaks);

    TGraph * gr = new TGraph();
    for(int i=0;i<nPeaks;i++){
        cout << x[i] << endl;
        gr->SetPoint(i,x[i],i+32);
    }

    // Define and fit a function (linear in this case)
    TF1 *fitFunc = new TF1("fitFunc", "pol5");
    gr->Fit(fitFunc, "Q"); // Quiet mode

    // Open a text file to save fit results
    std::ofstream outFile("Output/Calibration.txt");
    std::ofstream outFileCal("../../../../Calibration/VAMOS/CHIO/Chio_Z_Calibration.cal");

    outFileCal << "IC_Z_CALIBRATION" << " ";
    if (outFile.is_open()) {
        for (int i = 0; i < fitFunc->GetNpar(); i++) {
            outFile <<  fitFunc->GetParameter(i) << " ";
            outFileCal <<  fitFunc->GetParameter(i) << " ";
        }
        outFile.close();
        outFileCal.close();
        std::cout << "Fit results saved to fit_results.txt\n";
    }

    std::ofstream CalFile("../../../../Calibration/VAMOS/CHIO/Chio_Z_Spline_Eval.txt");
    if (CalFile.is_open()) {
        for (int i = 0; i < cutZ.size(); i++) {
            vector<double> YEdgesMinTemp = GetMinMaxY(cutZ.at(i));
            CalFile << MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2) << " ";
        }
        CalFile.close();
        std::cout << "Eval position saved in calibration\n";
    }



    TCanvas * cangr = new TCanvas("Cal","Cal",200,0,1200,1000);
    cangr->cd();
    gr->SetMarkerStyle(20);
    gr->Draw("AP");
    gr->Fit("pol5");

    // ******************************************************************************
    TCanvas * can = new TCanvas(Form("ChioEbis_vs_ChioDE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            Form("ChioEbis_vs_ChiodE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            0,0,2000,2000);
    can->Divide(1,3);
    can->cd(1); gPad-> SetLogz(); hChioDE_E_all->Draw("colz");
    for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
    can->cd(2); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz");
    can->cd(3); gPad-> SetLogz(); hChioZ_E_rough->Draw("colz");

    TCanvas * can1D = new TCanvas(Form("DE_and_Z_ifE_8000to25000_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            Form("DE_and_Z_ifE_8000to25000_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
            0, 0, 2000,2000);
    can1D->Divide(1,3);
    can1D->cd(1);  hChioDE_corr->Draw(); 
    can1D->cd(2);  hChioZ_rough->Draw();
    can1D->cd(3);  hChioZ_abitbetter->Draw();

    TFile * fsave = new TFile(Form("Output/DEvsE_corr_run%04d_%04d.root",RunNumber[0],RunNumber[Nrun-1]),"recreate");  
    can->Write();
    hChioDE_E_all->Write();
    hChioDE_E_corr->Write();
    hChioZ_E_rough->Write();
    can1D->Write();
    hChioDE_corr->Write();
    hChioZ_rough->Write();
    hChioZ_abitbetter->Write();
    fsave->Close();
}

///////////////////////////////////////////////////////////////////////////
void Fit_DE_E_corr()
{

    TFile * f = new TFile(Form("Output/DEvsE_corr_run%04d_%04d.root",RunNumber[0],RunNumber[Nrun-1]),"read");
    f->ls();
    TH1F * h1DEcorr = (TH1F*)f->Get("hChioDE_corr");
    h1DEcorr->SetDirectory(0);
    f->Close();

    TCanvas * canfit = new TCanvas("DEcorrforfit","Decorrforfit",0,0,3000,1000);
    canfit->Divide(2,1);
    canfit->cd(1);   h1DEcorr->Draw();

    float InitMean[34] = {
        7.55521e+03,     7.82392e+03,    8.04283e+03,     8.36164e+03,   
        8.68128e+03,     9.06218e+03,    9.44613e+03,     9.82352e+03,   
        1.02433e+04,     1.06666e+04,    1.11215e+04,     1.15371e+04,   
        1.20246e+04,     1.24492e+04,    1.29190e+04,     1.34056e+04,   
        1.38367e+04,     1.42968e+04,    1.46816e+04,     1.51614e+04,   
        1.55115e+04,     1.58873e+04,    1.63049e+04,     1.66187e+04,   
        1.69730e+04,     1.73969e+04,    1.76965e+04,     1.80730e+04,   
        1.85234e+04,     1.88849e+04,    1.92127e+04,     1.96378e+04,   
        1.99820e+04,     2.02914e+04
    };

    float InitSigma[34] = {
        6.97215e+01,    7.76394e+01,    7.84506e+01,    8.23239e+01, 
        9.25943e+01,    1.06154e+02,    1.03566e+02,    1.02429e+02, 
        1.10066e+02,    1.15032e+02,    1.20449e+02,    1.15397e+02, 
        1.30322e+02,    1.29513e+02,    1.45757e+02,    1.53552e+02, 
        1.48061e+02,    1.49866e+02,    1.60311e+02,    1.56691e+02, 
        1.55461e+02,    1.78093e+02,    1.66829e+02,    1.67863e+02, 
        1.75509e+02,    1.75189e+02,    1.70589e+02,    1.76178e+02, 
        1.65192e+02,    1.87470e+02,    1.78153e+02,    1.58322e+02, 
        1.53678e+02,    1.97134e+02
    };

    char name[1000];
    TString totname;
    for(int i=0; i<34;i++){
        sprintf(name,"%i",i*3);
        totname += "gaus(" + TString(name) + ")+";
    }  
    totname.Remove(totname.Last('+'));
    TF1 * gtot = new TF1("MultiGaussianFit",totname,InitMean[0]-2*InitSigma[0],InitMean[33]+2*InitSigma[33]);

    for(int i=0; i<34;i++){
        gtot->SetParameter(i*3+1,InitMean[i]); 
        gtot->SetParameter(i*3+2,InitSigma[i]); 
    }  
    h1DEcorr->Fit(gtot,"R+");

    TGraph * gr = new TGraph();
    gr->SetName("DEfit_vs_Z");
    gr->SetMarkerStyle(20);
    for(int i=0; i<33; i++){
        gr->SetPoint(i,sqrt(gtot->GetParameter(i*3+1)),i+32); 
    }

    canfit->cd(2); gr->Draw("AP");

}


const int CutNumberLoader(){
    TFile *fcut=new TFile(CutPath,"open");

    int CutCount = 0 ;
    TIter next(fcut->GetListOfKeys());
    TKey* key;

    while ((key=(TKey*)next())){
        if (std::string(key->GetClassName()) == "TCutG"){
            CutCount ++;
        }
    }

    const int res = CutCount;
    return res;
}


std::vector<double> GetMinMaxX(TCutG* cut) {
    if (!cut) {
        std::cerr << "Invalid TCutG!" << std::endl;
        return {};
    }

    int nPoints = cut->GetN();
    double* xPoints = cut->GetX();

    // Find the minimum and maximum X values
    double minX = *std::min_element(xPoints, xPoints + nPoints);
    double maxX = *std::max_element(xPoints, xPoints + nPoints);

    return {minX, maxX};
}

std::vector<double> GetMinMaxY(TCutG* cut) {
    if (!cut) {
        std::cerr << "Invalid TCutG!" << std::endl;
        return {};
    }

    int nPoints = cut->GetN();
    double* yPoints = cut->GetY();

    // Find the minimum and maximum X values
    double minX = *std::min_element(yPoints, yPoints + nPoints);
    double maxX = *std::max_element(yPoints, yPoints + nPoints);

    return {minX, maxX};
}


TF1* Linear(double x1, double y1, double x2, double y2, const char* name ) {
    // Ensure the points are not the same
    if (x1 == x2) {
        throw std::invalid_argument("x1 and x2 must be different to define a function.");
    }

    // Calculate slope and intercept of the line
    double slope = (y2 - y1) / (x2 - x1);
    double intercept = y1 - slope * x1;

    // Define the function as a lambda
    auto func = [slope, intercept](double* x, double* p) {
        return slope * x[0] + intercept;
    };

    // Create and return the TF1
    TF1* f = new TF1(name, func, x1, x2, 0);
    return f;
}