TestSpline.C

#include "TICPhysics.h"
#include <TCanvas.h>
#include <TChain.h>
#include <TF1.h>
#include <TFile.h>
#include <TH2.h>
#include <TSpline.h>
#include <fstream>
#include <vector>

using namespace std;

//////////////////////////////////////////////////////////////////////////////////////////
void TestSpline() {
  //===========================================================================================================
  //                           Loading var
  //===========================================================================================================
  TChain* chain = new TChain("PhysicsTree");
  chain->Add("../../../root/analysis/VamosCalib241.root");
  TICPhysics* IC = new TICPhysics() ;

  double FF_IC_X, FF_IC_Y;

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


  TFile* fHisto =  TFile::Open("Output/Histo.root");
  TFile* fSpline =  TFile::Open("Output/spline_Chio_2024.root");

  vector<TH2F*> hIC,hICorr;
  TH2F* hICorr_X = new TH2F("hICorr_X","hICorr_X",1000,-600,400,500,1,1.6);
  TH2F* hICorr_Y = new TH2F("hICorr_Y","hICorr_Y",1000,-100,100,500,1,1.6);
  TH2F* hDE_Ecorr = new TH2F("DE_Ecorr","DE_Ecorr",1000,100,1,1000,100,1);
  TH2F* hDEBis_Ecorr = new
    TH2F("DEBis_Ecorr","DEBis_Ecorr",1000,0,20000,1000,0,25000);
  TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,100,1,1000,100,1);
  TH2F* hDE_E_IC1_Corr = new
    TH2F("hDE_E_IC1_Corr","hDE_E_IC1_Corr",1000,0,20000,1000,0,25000);

  hICorr.push_back(hICorr_X);
  hICorr.push_back(hICorr_Y);
  hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_X_all"));
  hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_Y_all"));    

  vector<TSpline3*> spline;
  for (int i = 0; i < 4; i++) {
    
    if (i<2) spline.push_back((TSpline3*)fSpline->Get(Form("fsplineIC0_%d",i)));
    else if ( i>=2 && i<4) spline.push_back((TSpline3*)fSpline->Get(Form("fsplineIC1_%d",i)));
  } //End loop on histogram


  //===========================================================================================================
  //                                        Event Loop
  //===========================================================================================================
  int Nentries = chain->GetEntries();

  for (int e = 0; e < Nentries; e++) {
    // Printing status
    if (e % 100000 == 0){
      cout <<  "Please wait, we are  " << e * 100. / Nentries << "% done ..\r" << flush;
    }
    chain->GetEntry(e);

    double IC_CorrX = IC->fIC_raw[1]/ IC->fIC_raw[0]  * spline[0]->Eval(0) /
      spline[0]->Eval(FF_IC_X) ;

    double IC0_CorrX = IC->fIC_raw[1] / IC_CorrX ;

    double IC_CorrY = IC->fIC_raw[1]/ IC0_CorrX * spline[1]->Eval(0) /
      spline[1]->Eval(FF_IC_Y) ;

    hICorr[0]->Fill(FF_IC_X,IC_CorrX);
    hICorr[1]->Fill(FF_IC_Y,IC_CorrY);


    double IC0_Corr = IC->fIC_raw[1] /IC_CorrY;
    if (IC0_Corr != IC0_Corr ) IC0_Corr = 0;


    double IC1_CorrY = IC->fIC_raw[1] * spline[3]->Eval(0) /
      spline[3]->Eval(FF_IC_Y) ;
    double IC1_CorrX = IC->fIC_raw[1] * spline[2]->Eval(0) /
      spline[2]->Eval(FF_IC_X) ;


    double DE = 0 ,E =0 ,DEcorr=0 ,DEIC1corr =0;

    //===========================================================================================================
    //                                Fill De from seg 1 to 3
    //===========================================================================================================
    for (int seg = 0 ;  seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
      if(seg > 1 && seg < 4) DE += IC->fIC_raw[seg] ;
      else if (seg >4) E+= IC->fIC_raw[seg] ;
    }



    //===========================================================================================================
    //                                      Filling DE with ICO
    //===========================================================================================================
    DEcorr += DE ;

    if ((IC0_Corr < 1e100 && IC0_Corr >-1000) || true){ 
      DEcorr += IC0_Corr ;
    }
    else DEcorr += IC->fIC_raw[0];
    DE += IC->fIC_raw[0] ;

    DE = DE *0.5 + IC->fIC_raw[4];
    DEcorr = DEcorr *0.5 + IC->fIC_raw[4];

    /* 
       cout << IC0_Corr << "  " << IC->fIC_raw[0] << endl; 
       cout << DEcorr << "  " << DE << endl;
       cout << "*************************" << endl;
       */

    //===========================================================================================================
    //                            Chio online
    //===========================================================================================================

    double IC1corr = (IC->fIC_raw[1]*(1-0.000686068*FF_IC_Y))*(1-4.88238e-05*FF_IC_Y+7.40395e-06*FF_IC_Y*FF_IC_Y);
    double DE_Bis = 0.5*(IC1corr+IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
    //double DE_Bis = 0.5*(IC->fIC_raw[1]+IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
    DEIC1corr = 0.5*(IC1_CorrY +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];

    hDEBis_Ecorr->Fill(E,DE_Bis);
    hDE_E->Fill(E,DE);
    hDE_Ecorr->Fill(E,DEcorr);
    hDE_E_IC1_Corr->Fill(E,DEIC1corr);
  } //End loop on event

  //===========================================================================================================
  //                                          Drawing histo
  //===========================================================================================================

  TCanvas* c1 = new TCanvas("c1","c1",1000,1000);
  c1->Divide(2,2);
  c1->cd(1);
  hICorr[0]->Draw();
  c1->cd(2);
  hIC[0]->Draw();
  c1->cd(3);
  hICorr[1]->Draw();
  c1->cd(4);
  hIC[1]->Draw();

  TCanvas* c2 = new TCanvas("c2","c2",1500,1000);
  c2->Divide(4);
  c2->cd(1);
  hDE_E->Draw();
  c2->cd(2);
  hDE_Ecorr->Draw();
  c2->cd(3);
  hDEBis_Ecorr->Draw();
  c2->cd(4);
  hDE_E_IC1_Corr->Draw();

} // End spline chio XY