From 74fcacee4a190e398ba15aaef3152419e44b3bbf Mon Sep 17 00:00:00 2001
From: "theodore.efremov" <theodore.efremov@cea.fr>
Date: Tue, 5 Nov 2024 09:32:10 +0100
Subject: [PATCH] Added fit method and continued the analysis

---
 .../2024/macro/chio/CalibrationChio/ICCorr.C  | 147 ++++++++++++++----
 .../IC_Dampen_Recombination.cxx               |  98 ++++++++----
 .../CalibrationChio/IC_Dampen_Recombination.h |  10 ++
 .../macro/chio/CalibrationChio/SplineChioXY.C |  96 ++++++++----
 4 files changed, 258 insertions(+), 93 deletions(-)
 create mode 100644 Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.h

diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
index 80f5b4ab3..a12985b1e 100644
--- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
@@ -1,8 +1,4 @@
-#ifdef TICPHYSICS_H
-#define TICPHYSICS_H
-#include "TICPhysics.h"
-#endif
-
+#include <TICPhysics.h>
 #include <TCanvas.h>
 #include <TChain.h>
 #include <TF1.h>
@@ -14,6 +10,7 @@
 
 using namespace std;
 
+int GetNumberKey(TFile *infile ,const char* ClassName);
 //////////////////////////////////////////////////////////////////////////////////////////
 void ICCorr() {
   //===========================================================================================================
@@ -22,7 +19,6 @@ void ICCorr() {
   TChain* chain = new TChain("PhysicsTree");
   chain->Add("../../../root/analysis/VamosCalib241.root");
 
-
   TICPhysics* IC = new TICPhysics() ;
   double FF_IC_X, FF_IC_Y;
 
@@ -36,6 +32,7 @@ void ICCorr() {
 
   TFile* fHisto =  TFile::Open("Output/Histo.root");
   TFile* fSpline =  TFile::Open("Output/spline_Chio_2024.root");
+  TFile* fFit = TFile::Open("Output/FitIC.root");
 
   //Defining output histo
   vector<TH2F*> hIC,hICorr;
@@ -44,16 +41,31 @@ void ICCorr() {
 
   TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,100,1,1000,100,1);
   TH2F* hDE_Ecorr01 = new TH2F("DE_Ecorr01","DE_Ecorr01",1000,100,1,1000,100,1);
-  TH2F* hDEBis_Ecorr01 = new
-    TH2F("DEBis_Ecorr01","DEBis_Ecorr01",1000,0,20000,1000,0,25000);
+  TH2F* hDEBis_Ecorr_IC1 = new
+    TH2F("DEBis_Ecorr_IC1","DEBis_Ecorr_IC1",1000,0,20000,1000,0,25000);
   TH2F* hDE_E_IC1_Corr = new
     TH2F("hDE_E_IC1_Corr","hDE_E_IC1_Corr",1000,0,20000,1000,0,25000);
 
   TH2F* hDE_E_IC0_Corr = new
     TH2F("hDE_E_IC0_Corr","hDE_E_IC0_Corr",1000,0,20000,1000,0,25000);
- TH2F* hDEcut= 
+  TH2F* hDEcut= 
     new TH2F("DEcut","DEcut",1000,0,20000,1000,0,25000);
-    
+
+  TH2F* hDE_Ecorr_IC1 = new
+      TH2F("hDE_Ecorr_IC1","hDE_Ecorr_IC1",1000,0,20000,1000,0,25000);
+  
+  TH2F* hDE_Ecorr_IC0 = new
+      TH2F("hDE_Ecorr_IC0","hDE_Ecorr_IC0",1000,0,20000,1000,0,25000);
+
+  TH2F* hDE_Ecorr_IC0_IC1corr = new
+      TH2F("hDE_Ecorr_IC0_IC1corr","hDE_Ecorr_IC0_IC1corr",1000,0,20000,1000,0,25000);
+  //Fit
+  TH2F* hDE_E_ICFit_Corr = new
+      TH2F("hDE_E_ICFit_Corr","hDE_E_ICFit_Corr",1000,0,20000,1000,0,25000);
+
+  TH2F* hICorr_Y_Fit = new
+    TH2F("hICorr_Y_Fit","hICorr_Y_Fit",1000,-100,100,500,0,25000);
+
   hICorr.push_back(hICorr_X);
   hICorr.push_back(hICorr_Y);
 
@@ -61,12 +73,14 @@ void ICCorr() {
   hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_X_all"));
   hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_Y_all"));    
 
-
+  //===========================================================================================================
+  //                                            GETSPLINE
+  //===========================================================================================================
   // Get number of spline
   int SplineCount = 0 ;
   TIter next(fSpline->GetListOfKeys());
   TKey* key;
-  
+
   while ((key=(TKey*)next())){
     if (std::string(key->GetClassName()) == "TSpline3"){
       SplineCount ++;
@@ -78,13 +92,25 @@ void ICCorr() {
 
   vector<TSpline3*> spline;
   for (int i = 0; i < SplineCount  ; i++) {
-   cout << i << endl; 
+    cout << i << endl; 
     if (i<2) spline.push_back((TSpline3*)fSpline->Get(Form("fsplineIC01_%d",i)));
     else if ( i>=2 && i<4) spline.push_back((TSpline3*)fSpline->Get(Form("fsplineIC1_%d",i)));
 
     else if (i>=4 && i<6) spline.push_back((TSpline3 *)fSpline->Get(Form("fsplineIC0_%d",i)));
-  
+    else if (i==6) spline.push_back((TSpline3 *)fSpline->Get(Form("fsplineIC0_IC1corr")));
+
   } //End loop on histogram
+    //===========================================================================================================
+    //                                 get fit
+    //===========================================================================================================
+
+  int FitCount = GetNumberKey(fFit,"TF1"); 
+
+
+  vector<TF1*> Fit;
+  for (int i = 0; i < FitCount  ; i++) {
+    Fit.push_back((TF1*)fFit->Get(Form("Fit_Y_Expo_IC%d",i)));
+  } //End loop on Fit 
 
 
   //===========================================================================================================
@@ -118,8 +144,8 @@ void ICCorr() {
     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_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) ;
 
@@ -127,6 +153,11 @@ void ICCorr() {
     double IC0_CorrY = IC->fIC_raw[0] * spline[4]->Eval(0)/ spline[4]->Eval(FF_IC_Y);
     double IC0_CorrX = IC->fIC_raw[0] * spline[3]->Eval(0)/ spline[3]->Eval(FF_IC_X);
 
+    //Double correction
+    double IC0_CorrIC1_Y = IC->fIC_raw[1]/ IC_temp_CorrX * spline[6]->Eval(0) /
+      spline[6]->Eval(FF_IC_Y) ;
+    double IC0_Corr_IC1 = IC->fIC_raw[1] /IC0_CorrIC1_Y;
+
     double DE = 0 ,E =0 ,DEcorr=0 ,DEIC1corr =0, DE_IC0_Corr_Y=0;
 
     //===========================================================================================================
@@ -152,7 +183,7 @@ void ICCorr() {
     if ((IC0_CorrX < 1e100 && IC0_CorrX >-1000) || true){ 
       DE_IC0_Corr_Y += IC0_CorrX ;
     }
-       
+
     else DE_IC0_Corr_Y += IC->fIC_raw[0];
 
     DE += IC->fIC_raw[0] ;
@@ -176,21 +207,50 @@ void ICCorr() {
     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_Ecorr01->Fill(E,DE_Bis);
-    hDE_E->Fill(E,DE);
-    hDE_Ecorr01->Fill(E,DEcorr);
-    hDE_E_IC1_Corr->Fill(E,DEIC1corr);
- 
-   double DEcut;
+    double DEIC0corr = 0.5*(IC0_CorrY + IC1_CorrY +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+    double DEIC0_IC1corr = 0.5*(IC0_CorrIC1_Y + IC1_CorrY +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+
+    if (FF_IC_Y >-60 && FF_IC_Y <60){
+      hDEBis_Ecorr_IC1->Fill(E,DE_Bis);
+      hDE_Ecorr_IC1->Fill(E,DEIC1corr);
+      hDE_Ecorr_IC0->Fill(E,DEIC0corr);
+      hDE_Ecorr_IC0_IC1corr->Fill(E,DEIC0_IC1corr);
+
+      hDE_E->Fill(E,DE);
+      hDE_Ecorr01->Fill(E,DEcorr);
+      hDE_E_IC1_Corr->Fill(E,DEIC1corr);
+    }
+    double DEcut;
     if ((FF_IC_Y>-60 && FF_IC_Y<-55)|| (FF_IC_Y>55 && FF_IC_Y<60)){  DEcut=DE_Bis;    
-    hDEcut->Fill(E,DEcut);
+      hDEcut->Fill(E,DEcut);
     }
-  } //End loop on event
+    //===========================================================================================================
+    //
+    //===========================================================================================================
+    //===========================================================================================================
+    //                                            Fit Corr
+    //===========================================================================================================
 
-  //===========================================================================================================
-  //                                          Drawing histo
-  //===========================================================================================================
+    vector<double> IC_Corr_Fit(11);
+    double DEFit=0, EtotFit=0, EFit=0;
+    for (int seg = 0 ;  seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+      IC_Corr_Fit[seg] = IC->fIC_raw[seg]* Fit[seg]->Eval(0) /Fit[seg]->Eval(FF_IC_Y);
+    }
+
+    DEFit = (IC_Corr_Fit[1] + IC->fIC_raw[2]+ IC->fIC_raw[3])*0.5
+      + IC->fIC_raw[4];
+
+    EFit = IC->fIC_raw[5]+ IC->fIC_raw[6]+ IC->fIC_raw[7]+ IC->fIC_raw[8]+
+      IC->fIC_raw[9]+ IC->fIC_raw[10];
+
+    EtotFit =  DEFit + EFit ;
+
+    hICorr_Y_Fit->Fill(FF_IC_Y,IC_Corr_Fit[0]);
+    hDE_E_ICFit_Corr ->Fill(EFit,DEFit);
+  } //End loop on event
+    //===========================================================================================================
+    //                                          Drawing histo
+    //===========================================================================================================
 
   TCanvas* c1 = new TCanvas("c1","c1",1000,1000);
   c1->Divide(2,2);
@@ -206,12 +266,33 @@ void ICCorr() {
   TCanvas* c2 = new TCanvas("c2","c2",1500,1000);
   c2->Divide(4);
   c2->cd(1);
-  hDE_E->Draw();
+  gPad->SetLogz();
+  hDE_Ecorr_IC0_IC1corr->Draw();
   c2->cd(2);
-  hDE_Ecorr01->Draw();
+  gPad->SetLogz();
+  hDE_Ecorr_IC1->Draw();
   c2->cd(3);
-  hDEBis_Ecorr01->Draw();
+  gPad->SetLogz();
+  hDEBis_Ecorr_IC1->Draw();
   c2->cd(4);
-  hDEcut->Draw();
+  gPad->SetLogz();
+  hDE_E_ICFit_Corr->Draw();
 
 } // End spline chio XY
+
+
+int GetNumberKey(TFile *infile ,const char* ClassName){
+  // Get number of spline
+  int KeyCount = 0 ;
+  TIter next(infile->GetListOfKeys());
+  TKey* keyFit;
+
+  while ((keyFit=(TKey*)next())){
+    if (std::string(keyFit->GetClassName()) == ClassName){
+      KeyCount ++;
+    }
+  }
+
+  return KeyCount;
+
+}
diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx
index 05ce6484b..56e902673 100644
--- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx
@@ -1,33 +1,27 @@
-#include <TFile.h>
-#include <TICPhysics.h>
-#include <TChain.h>
-#include <TF1.h>
-#include <TH2.h>
-#include <Math/Functor.h> 
-#include <TTreeReaderValue.h>
-#include <vector>
+#include "IC_Dampen_Recombination.h"
 using namespace std;
 
 
-void HistoFillerIC();
+void HistoFillerIC(bool cut = false);
 
 void IC_Dampen_Recombination(bool Fill = false){
 
 
-  if (Fill==true) HistoFillerIC();
+  if (Fill==true) HistoFillerIC(false);
+
 
- 
   TFile *infile = TFile::Open("Output/HistoDampen.root"); 
+  TFile *OutFile = new TFile("Output/FitIC.root","recreate");
   vector<TH2F*> hICX, hICY ;
   vector<TH1F*> hICX_Profile, hICY_Profile ;
-  
-   for (int seg = 0 ; seg < 11; seg++ ){
-   
+
+  for (int seg = 0 ; seg < 11; seg++ ){
+
     hICX.push_back((TH2F*)infile->Get(Form("hChio_X_%d",seg)));
     hICY.push_back((TH2F*)infile->Get(Form("hChio_Y_%d",seg)));
     hICX_Profile.push_back((TH1F*)infile->Get(Form("hChio_X_%d_pfx",seg)));
     hICY_Profile.push_back((TH1F*)infile->Get(Form("hChio_Y_%d_pfx",seg)));
-  
+
   }
 
   //===========================================================================================================
@@ -37,27 +31,40 @@ void IC_Dampen_Recombination(bool Fill = false){
 
 
   vector<TF1*> Fit_Y(11);
+  double p0[11],p1[11],p2[11],p3[11], xmin[11];
+  
+  p0[0] = 2350  ; p1[0] = 0.00065  ; p3[0] = -54  ; xmin[0] = -60; 
+  p0[1] = 2350  ; p1[1] = 0.00065  ; p3[1] = -54  ; xmin[0] = -58; 
+  p0[2] = 2350  ; p1[2] = 0.00065  ; p3[2] = -54  ; xmin[0] = -60; 
+  p0[3] = 2350  ; p1[3] = 0.00065  ; p3[3] = -54  ; xmin[0] = -60; 
+  p0[4] = 2350  ; p1[4] = 0.00065  ; p3[4] = -54  ; xmin[0] = -60; 
+  p0[5] = 2350  ; p1[5] = 0.00065  ; p3[5] = -54  ; xmin[0] = -60; 
+  p0[6] = 2350  ; p1[6] = 0.00065  ; p3[6] = -54  ; xmin[0] = -60; 
+  p0[7] = 2350  ; p1[7] = 0.00065  ; p3[7] = -54  ; xmin[0] = -60; 
+  p0[8] = 2350  ; p1[8] = 0.00065  ; p3[8] = -54  ; xmin[0] = -60; 
+  p0[9] = 2350  ; p1[9] = 0.00065  ; p3[9] = -54  ; xmin[0] = -60; 
+  p0[10]= 2350  ; p1[10]= 0.00065  ; p3[10]= -54  ; xmin[0] = -60; 
+
+  OutFile->cd();
   for (int seg=0; seg<11 ; seg++){
 
     Fit_Y[seg] = new
       TF1(Form("Fit_Y_Expo_IC%d",seg),[](double *x, double*p) {
 
-          if (x[0]>-50 && x[0]<50) {
-            TF1::RejectPoint();
+          if (/*x[0]>-50 && x[0]<50 && x[0] > 60 && x[0] < -60*/ false) {
+          TF1::RejectPoint();
           }
           //Define here the function
-          return p[0] * exp(p[1]* (x[0]-p[3])) + p[2] ;
-      },-60,63,4);
+          return p[0] * exp(p[1]* (x[0])) ;
+          },-60,60,2);
 
 
-    Fit_Y[seg]->SetParameters(1,-10,4800,-60);
+    Fit_Y[seg]->SetParameters(p0[seg],p1[seg],p3[seg]);
 
 
     cout << "****************" << endl << "IC"<<seg << endl;
-    hICY_Profile[seg]->Fit(Fit_Y[seg] ,"Rnone");  vector<TH2F*> hICX, hICY ;
-    vector<TH1F*> hICX_Profile, hICY_Profile ;
-
-
+    hICY_Profile[seg]->Fit(Fit_Y[seg] ,"R"); 
+    Fit_Y[seg]->Write();
 
   }//end loop on histo
 
@@ -68,11 +75,11 @@ void IC_Dampen_Recombination(bool Fill = false){
   //===========================================================================================================
 
   vector<TCanvas*> cIC(11);
-  for (int seg=0; seg<1; seg++){
+  for (int seg=0; seg<2; seg++){
     cIC[seg] = new TCanvas
       (Form("cIC%d",seg),Form("cIC%d",seg),100,100,1000,600);
     cIC[seg]->Divide(2);
-    
+
     cIC[seg]->cd(1);
     hICY[seg]->Draw();
 
@@ -82,10 +89,12 @@ void IC_Dampen_Recombination(bool Fill = false){
 
 
   }
+
+  OutFile->Close();
 }
 
 
-void HistoFillerIC(){
+void HistoFillerIC(bool cut){
   //===========================================================================================================
   //                                              INPUT
   //===========================================================================================================
@@ -104,6 +113,10 @@ void HistoFillerIC(){
   chain->SetBranchAddress("IC", &IC);
 
 
+  TFile *fcut = new TFile("Output/CutZ.root","open") ;
+  TCutG *CutZ ;
+  CutZ = (TCutG*)fcut->Get("CutTry") ;
+
   //===========================================================================================================
   //                                        Declare and fill  histo
   //===========================================================================================================
@@ -115,8 +128,8 @@ void HistoFillerIC(){
 
   for (int seg = 0 ; seg < 11 ; seg++){
 
-    TH2F* hChio_X = new TH2F(Form("hChio_X_%d",seg) ,Form("hChio_X_%d",seg) ,1000 ,-800 ,800 ,500 ,2000,9500 );
-    TH2F* hChio_Y = new TH2F(Form("hChio_Y_%d",seg) ,Form("hChio_Y_%d",seg) ,1000 ,-160 ,160 ,500 ,2000,9500 );
+    TH2F* hChio_X = new TH2F(Form("hChio_X_%d",seg) ,Form("hChio_X_%d",seg) ,500 ,-800 ,800 ,500 ,2000,9500 );
+    TH2F* hChio_Y = new TH2F(Form("hChio_Y_%d",seg) ,Form("hChio_Y_%d",seg) ,500 ,-160 ,160 ,500 ,2000,9500 );
 
     TH1F* hProfile_X = (TH1F*)hChio_X->ProfileX();
     TH1F* hProfile_Y = (TH1F*)hChio_Y->ProfileX();
@@ -140,15 +153,32 @@ void HistoFillerIC(){
 
     chain->GetEntry(e);
 
-    for (int seg=0 ; seg<11 ; seg++ ){
 
-      hICX[seg]->Fill(FF_IC_X,IC->fIC_raw[seg]);
-      hICY[seg]->Fill(FF_IC_Y,IC->fIC_raw[seg]);
+    //if we want to do cut on DE_E we need to define those parameters
+    double DE = 0 ,E = 0, Etot = 0 ;
 
-    }//end loop seg
+    for (int seg = 0 ;  seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+      if (seg >4) E+= IC->fIC_raw[seg] ;
+    }
 
-  } //end loop entries
+    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];
 
+    Etot = E+DE_Bis;
+
+    //===========================================================================================================
+    //                                  Fill histo
+    //===========================================================================================================
+    if (cut == true && !(CutZ->IsInside(Etot,DE_Bis)) ){ ; } // end cut
+    else {
+      for (int seg=0 ; seg<11 ; seg++ ){
+
+        hICX[seg]->Fill(FF_IC_X,IC->fIC_raw[seg]);
+        hICY[seg]->Fill(FF_IC_Y,IC->fIC_raw[seg]);
+
+      }//end loop seg
+    }
+  } //end loop entries
 
   //ReProfile
 
diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.h b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.h
new file mode 100644
index 000000000..2a2b36437
--- /dev/null
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.h
@@ -0,0 +1,10 @@
+#include <TFile.h>
+#include <TICPhysics.h>
+#include <TChain.h>
+#include <TF1.h>
+#include <TH2.h>
+#include <Math/Functor.h> 
+#include <TTreeReaderValue.h>
+#include <vector>
+#include<TCutG.h>
+
diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C
index 925316654..7b9ec8f1e 100644
--- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C
@@ -9,13 +9,14 @@
 #include <TSpline.h>
 #include <fstream>
 #include <vector>
+#include <TCutG.h>
 
 using namespace std;
 
 TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double
     XMax, double YMin, double YMax);
-void HistoFiller();
-void HistoDrawer(vector<TH2F*> h);
+void HistoFiller(bool cut=true);
+void HistoDrawer(vector<TH2F*> h ,const char* CharName);
 
 //////////////////////////////////////////////////////////////////////////////////////////
 void SplineChioXY( bool Create = false) {
@@ -24,7 +25,7 @@ void SplineChioXY( bool Create = false) {
 
   if (Create == true){
     cout << "Creating Histo.root" << endl;
-    HistoFiller(); //Deactivated after one launch of the macro
+    HistoFiller(false); //Deactivated after one launch of the macro
     cout << "Histo.root created and loaded" << endl ;
   }
 
@@ -38,15 +39,26 @@ void SplineChioXY( bool Create = false) {
   }
 
 
-  vector<TH2F*> hIC;
+  vector<TH2F*> hIC,hICX,hICY;
   hIC.push_back((TH2F*)inFile->Get("hChio_IC0_IC1_X_all"));
   hIC.push_back((TH2F*)inFile->Get("hChio_IC0_IC1_Y_all"));
   hIC.push_back((TH2F*)inFile->Get("hChio_IC1_X"));
   hIC.push_back((TH2F*)inFile->Get("hChio_IC1_Y"));
   hIC.push_back((TH2F*)inFile->Get("hChio_IC0_X"));
-  hIC.push_back((TH2F*)inFile->Get("hChio_IC0_Y"));
+  hIC.push_back((TH2F*)inFile->Get("hChio_IC0_Y")); 
+  hIC.push_back((TH2F*)inFile->Get("hChio_IC0_IC1corr_Y_all"));
 
-  HistoDrawer(hIC);
+  hICX.push_back((TH2F*)inFile->Get("hChio_IC0_IC1_X_all"));
+  hICY.push_back((TH2F*)inFile->Get("hChio_IC0_IC1_Y_all"));
+  hICX.push_back((TH2F*)inFile->Get("hChio_IC1_X"));
+  hICY.push_back((TH2F*)inFile->Get("hChio_IC1_Y"));
+  hICX.push_back((TH2F*)inFile->Get("hChio_IC0_X"));
+  hICY.push_back((TH2F*)inFile->Get("hChio_IC0_Y")); 
+  hICY.push_back((TH2F*)inFile->Get("hChio_IC0_IC1corr_Y_all"));
+
+
+  HistoDrawer(hICX,"CanX");
+  HistoDrawer(hICY,"CanY");
 
 
   TFile* fspline = new TFile("Output/spline_Chio_2024.root", "recreate");
@@ -54,11 +66,12 @@ void SplineChioXY( bool Create = false) {
   int NHist = static_cast<int>(hIC.size());
   double XMin[NHist] , XMax[NHist], YMin[NHist] , YMax[NHist];
   XMin[0] = -520 ;   XMax[0] = 380 ; YMin[0] = 0 ; YMax[0] = 2 ; 
-  XMin[1] = -70  ;   XMax[1] = 79  ; YMin[1] = 0 ; YMax[1] = 2 ; 
+  XMin[1] = -60  ;   XMax[1] = 63  ; YMin[1] = 0 ; YMax[1] = 2 ; 
   XMin[2] = -520 ;   XMax[2] = 380 ; YMin[2] = 2000 ; YMax[2] = 9500 ; 
-  XMin[3] = -70  ;   XMax[3] = 79  ; YMin[3] = 2000 ; YMax[3] = 9500 ; 
+  XMin[3] = -60  ;   XMax[3] = 63  ; YMin[3] = 2000 ; YMax[3] = 9500 ; 
   XMin[4] = -520 ;   XMax[4] = 380 ; YMin[4] = 2000 ; YMax[4] = 9500 ; 
-  XMin[5] = -70  ;   XMax[5] = 79  ; YMin[5] = 2000 ; YMax[5] = 9500 ; 
+  XMin[5] = -60  ;   XMax[5] = 63  ; YMin[5] = 2000 ; YMax[5] = 9500 ; 
+  XMin[6] = -60  ;   XMax[6] = 63  ; YMin[6] = 0 ; YMax[6] = 2 ; 
 
 
 
@@ -67,19 +80,21 @@ void SplineChioXY( bool Create = false) {
     TSpline3* spline = MakeSpline(hIC[i] , i, XMin[i] , XMax[i], YMin[i], YMax[i]);
     if (i < 2) spline->SetName(Form("fsplineIC01_%d", i ));
     else if (i>=2 && i<4) spline->SetName(Form("fsplineIC1_%d",i));   
-    else if (i>=4) spline->SetName(Form("fsplineIC0_%d",i)); 
+    else if (i>=4 && i<6) spline->SetName(Form("fsplineIC0_%d",i)); 
+    else if (i==6) spline->SetName(Form("fsplineIC0_IC1corr")); 
     spline->Write();
   } //End loop on histogram
 
   fspline->Close();
 
+
 } // End spline chio XY
 
 ///////////////////////////////////////////////////////////////////////////////////////////
 /**
  * This function fill the histograms
  */
-void HistoFiller() {
+void HistoFiller(bool cut) {
 
   // Input and setters
   TChain* chain = new TChain("PhysicsTree");
@@ -95,6 +110,15 @@ void HistoFiller() {
   chain->SetBranchStatus("IC", true);
   chain->SetBranchAddress("IC", &IC);
 
+  //cut
+  TFile *fcut = new TFile("Output/CutZ.root","open") ;
+  TCutG *CutZ ;
+  CutZ = (TCutG*)fcut->Get("CutTry") ;
+
+  //spline to correct IC1 as a baseline of IC0
+  TFile* fSpline =  TFile::Open("Output/spline_Chio_2024.root");
+  TSpline3* SplineIC1 = (TSpline3 *)fSpline->Get(Form("fsplineIC1_3"));
+
   // Output of the spline
 
   double fIC0_CorrX;
@@ -115,7 +139,8 @@ void HistoFiller() {
   TH2F* hChio_IC0_Y = new TH2F("hChio_IC0_Y",
       "hChio_IC0_Y", 1000, -160, 160, 500, 2000, 9500);
 
-
+  TH2F* hChio_IC0_IC1corr_Y_all = new TH2F("hChio_IC0_IC1corr_Y_all",
+      "hChio_IC0_IC1corr_Y_all", 1000, -160, 160, 500, 0.8, 1.6);
 
   //===========================================================================================================
   //																	Beginning loop on entries
@@ -138,6 +163,22 @@ void HistoFiller() {
       fIC_raw[i] = IC->fIC_raw[i]; 
     }
 
+    //if we want to do cut on DE_E we need to define those parameters
+    double DE = 0 ,E = 0, Etot = 0 ;
+
+    for (int seg = 0 ;  seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+     if (seg >4) E+= IC->fIC_raw[seg] ;
+    }
+
+    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];
+
+    Etot = E+DE_Bis;
+
+
+    if (cut == true && !(CutZ->IsInside(Etot,DE_Bis)) ){; } // end cut
+    else {
+
     //if (IC->Eres>11000 && IC->Eres < 12000){
     hChio_IC0_IC1_X_all->Fill(FF_IC_X , fIC_raw[1] / fIC_raw[0] );
     hChio_IC0_IC1_Y_all->Fill(FF_IC_Y , fIC_raw[1] / fIC_raw[0] );
@@ -147,8 +188,16 @@ void HistoFiller() {
 
     hChio_IC0_X->Fill(FF_IC_X, fIC_raw[0]);
     hChio_IC0_Y->Fill(FF_IC_Y, fIC_raw[0]);
-    //}
-
+    
+    double IC1_Spline_Corr = IC->fIC_raw[1] * SplineIC1->Eval(0) / SplineIC1->Eval(FF_IC_Y);
+    hChio_IC0_IC1corr_Y_all->Fill(FF_IC_Y , IC1_Spline_Corr / fIC_raw[0] );
+
+    /*
+    cout << "*************************" << endl; 
+    cout << IC1_Spline_Corr/fIC_raw[0] << endl;
+    cout << FF_IC_Y << endl;
+    */
+    }
   } // end loop on entries
 
   //===========================================================================================================
@@ -162,6 +211,7 @@ void HistoFiller() {
   hChio_IC1_X->Write();
   hChio_IC0_X->Write();
   hChio_IC0_Y->Write();
+  hChio_IC0_IC1corr_Y_all->Write();
   outFile->Close();
 } // End HistoFiller
 
@@ -262,24 +312,18 @@ TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double
 
 } // end makespline
 
-void HistoDrawer(vector<TH2F*> h){
+void HistoDrawer(vector<TH2F*> h,const char* CharName){
 
   int NHisto = static_cast<int>(h.size());
 
-  TCanvas* CanX= new TCanvas("CanX","CanX",0,0,2000,1000);
-  TCanvas* CanY= new TCanvas("CanY","CanY",0,0,2000,1000);
+  TCanvas* Can= new TCanvas(CharName,CharName,0,0,2000,1000);
 
-  CanX->Divide(NHisto/2);
-  CanY->Divide(NHisto/2);
+  Can->Divide(NHisto);
 
-  for ( int i=0 ; i < NHisto ; i+=2 ){
+  for ( int i=0 ; i < NHisto ; i+=1 ){
     //All x are even and all y are odd
-    int CanvaNumber = (i+1)/2 + 1 ;
-    CanX->cd(CanvaNumber);
+    int CanvaNumber = (i) ;
+    Can->cd(CanvaNumber);
     h[i]->Draw("colz");
-
-    CanY->cd(CanvaNumber);
-    h[i+1]->Draw("colz");
-
   }
 }
-- 
GitLab