From e457caabdc75c30aeda47b93e66b70602da8d9f9 Mon Sep 17 00:00:00 2001
From: "theodore.efremov" <theodore.efremov@cea.fr>
Date: Thu, 24 Oct 2024 16:01:01 +0200
Subject: [PATCH] Updated Spline XY code

---
 .../macro/chio/CalibrationChio/SplineChioXY.C | 146 ++++++++----------
 1 file changed, 67 insertions(+), 79 deletions(-)

diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C
index c5e8d8af8..24e712216 100644
--- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioXY.C
@@ -10,16 +10,33 @@
 
 using namespace std;
 
-TSpline3* MakeSpline(TH2F* hInput);
+TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee);
 void ApplySplinePerCut();
 void ApplySpline();
 void HistoFiller();
 
 //////////////////////////////////////////////////////////////////////////////////////////
 void SplineChioXY() {
-  HistoFiller();
 
-  TFile* inFile = new TFile("Output/Spline.root");
+  TFile* inFile =  TFile::Open("Output/Histo.root");
+  bool Create = false;
+
+  if (Create == true){
+    cout << "Creating Histo.root" << endl;
+    HistoFiller(); //Deactivated after one launch of the macro
+    cout << "Histo.root created and loaded" << endl ;
+    }
+
+
+  else if (inFile && !inFile->IsZombie()) {
+    cout << "Histogram loaded" << endl;
+  }
+
+  else {
+   cerr << "Make your mind" << endl;
+  }
+
+
   vector<TH2F*> hIC;
   hIC.push_back((TH2F*)inFile->Get("hChio_IC0_IC1_X_all"));
   hIC.push_back((TH2F*)inFile->Get("hChio_IC0_IC1_Y_all"));
@@ -28,14 +45,13 @@ void SplineChioXY() {
 
   for (int i = 0; i < hIC.size(); i++) {
 
-    TSpline3* spline = MakeSpline(hIC[0]);
-    spline->SetName(Form("fspline_%d", i + 1));
+    TSpline3* spline = MakeSpline(hIC[i] , i);
+    spline->SetName(Form("fspline_%d", i ));
     spline->Write();
   } //End loop on histogram
 
   fspline->Close();
 
-
   /*
   //===========================================================================================================
   //                                    Plot Histograms
@@ -52,20 +68,20 @@ void SplineChioXY() {
   cout << "Float_t F[" << Ncuts << "] = { " << endl;
   ofile_par << "Float_t FF_DEcorr[" << Ncuts << "] = { " << endl;
   for (int i = 0; i < Ncuts; i++) {
-    can->cd();
-    if (i == 0)
-      h2->Draw("col");
-    else
-      h2->Draw("col,same");
-    cout << h2->GetMean(2) << ", " << endl;
-    ofile_par << h2->GetMean(2) << ", " << endl;
-    // gr->SetPoint(i,sqrt(h2->GetMean(2)),i+31);
-    gr->SetPoint(i, gspline->Eval(8500), i + 31);
-    // hProfile->Draw("same");
-    pfx->Draw("same");
-    gspline->SetLineColor(kBlue);
-    gspline->SetLineWidth(3);
-    gspline->Draw("lsame");
+  can->cd();
+  if (i == 0)
+  h2->Draw("col");
+  else
+  h2->Draw("col,same");
+  cout << h2->GetMean(2) << ", " << endl;
+  ofile_par << h2->GetMean(2) << ", " << endl;
+  // gr->SetPoint(i,sqrt(h2->GetMean(2)),i+31);
+  gr->SetPoint(i, gspline->Eval(8500), i + 31);
+  // hProfile->Draw("same");
+  pfx->Draw("same");
+  gspline->SetLineColor(kBlue);
+  gspline->SetLineWidth(3);
+  gspline->Draw("lsame");
   }
 
   TCanvas* cangr = new TCanvas("ToGetRoughCal", "ToGetRoughCal", 200, 0, 1200, 1000);
@@ -73,7 +89,7 @@ void SplineChioXY() {
   gr->SetMarkerStyle(20);
   gr->Draw("AP");
   gr->Fit("pol2");
- */
+  */
 } // End spline chio XY
 
 ///////////////////////////////////////////////////////////////////////////////////////////
@@ -103,8 +119,10 @@ void HistoFiller() {
 
   // Histograms
 
-  TH2F* hChio_IC0_IC1_X_all = new TH2F("hChio_IC0_IC1_X_all", "hChio_IC0_IC1_X_all", 1000, -660, 400, 500, 0, 2);
-  TH2F* hChio_IC0_IC1_Y_all = new TH2F("hChio_IC0_IC1_Y_all", "hChio_IC0_IC1_Y_all", 1000, -70, 70, 500, 0, 2);
+  TH2F* hChio_IC0_IC1_X_all = new TH2F("hChio_IC0_IC1_X_all",
+      "hChio_IC0_IC1_X_all", 1000, -600, 400, 500, 0, 2);
+  TH2F* hChio_IC0_IC1_Y_all = new TH2F("hChio_IC0_IC1_Y_all",
+      "hChio_IC0_IC1_Y_all", 1000, -100, 100, 500, 0, 2);
 
   //===========================================================================================================
   //																	Beginning loop on entries
@@ -116,7 +134,7 @@ void HistoFiller() {
 
     // Printing status
     if (e % 100000 == 0){
-      cout << e * 100. / Nentries << "% \r" << flush;
+      cout <<  "Please wait, we are  " << e * 100. / Nentries << "% done ..\r" << flush;
     }
     chain->GetEntry(e);
 
@@ -126,7 +144,7 @@ void HistoFiller() {
     for (int i=0 ; i<11 ; i++){
       fIC_raw[i] = IC->fIC_raw[i]; 
     }
-       
+
     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] );
 
@@ -146,13 +164,14 @@ void HistoFiller() {
   canY->cd();
   hChio_IC0_IC1_Y_all->Draw("colz");
 
-  TFile* outFile = new TFile("Output/Spline.root", "recreate");
+  TFile* outFile = new TFile("Output/Histo.root", "recreate");
   hChio_IC0_IC1_X_all->Write();
   hChio_IC0_IC1_Y_all->Write();
   outFile->Close();
 } // End HistoFiller
 
 /**
+  }
  *  Create the spline for one TH2F
  *
  *  Input :
@@ -165,14 +184,16 @@ void HistoFiller() {
  *
  */
 
-TSpline3* MakeSpline(TH2F* hInput) {
+TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee) {
 
   TSpline3* gspline;
 
   TH1F* hProfile;
   TH1F* pfx; // extended hProfile
 
-  TCanvas* canfit = new TCanvas("canfit", "canfit", 0, 0, 2000, 1500);
+  TCanvas* canfit = new TCanvas(Form("CanSpline_%02d",NCallee),
+      Form("canfit_%02d",NCallee), 0, 0, 1000, 500);
+  
 
   //===========================================================================================================
   //                                            Init profile
@@ -186,70 +207,43 @@ TSpline3* MakeSpline(TH2F* hInput) {
 
   hProfile->GetYaxis()->SetRangeUser(0, 2);
   hProfile->Draw();
-
-  // Get the range of the TProfile
+  //===========================================================================================================
+  //                                          First and last bin to get to 6k
+  //                                          event get promoted
+  //===========================================================================================================
   int FirstBin, LastBin;
-  Double_t parpol3[4];
-
-  for (int bin = 1; bin < hProfile->GetNbinsX(); bin++) {
+  double Treshold = 1;
+  for (int bin =1 ; bin<hProfile->GetNbinsX(); bin++) {
     FirstBin = bin;
-    if (hProfile->GetBinContent(bin) > 6000)
-      break;
+    if (hProfile->GetBinContent(bin)>Treshold) break;
   }
-
-
-  Double_t* parpol1;
-
-  for (int bin = hProfile->GetNbinsX(); bin > 1; bin--) {
+  for (int bin = hProfile->GetNbinsX(); bin>1 ; bin--) {
     LastBin = bin;
-    if (hProfile->GetBinContent(bin) > 6000)
-      break;
+    if (hProfile->GetBinContent(bin)>Treshold) break;
   }
-
+  cout << FirstBin << " " << LastBin << endl;
   //===========================================================================================================
-  //                                          Init Extended profile
+  //                                          Init Extended profile function
   //===========================================================================================================
   // Create the extended TProfile
-  int i=0;
   pfx = new TH1F(
-      Form("pfx_%02d", i + 1), Form("pfx_%02d", i + 1), hProfile->GetNbinsX(), hProfile->GetBinLowEdge(1),
+      Form("pfx_%02d", NCallee + 1), Form("pfx_%02d", NCallee + 1), hProfile->GetNbinsX(), hProfile->GetBinLowEdge(1),
       hProfile->GetBinLowEdge(hProfile->GetNbinsX()) + hProfile->GetBinWidth(hProfile->GetNbinsX()));
   pfx->SetLineColor(8);
   pfx->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->GetBinLowEdge(FirstBin),
-      hProfile->GetBinLowEdge(FirstBin) + 10000);
-  hProfile->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->GetBinLowEdge(LastBin) - 10000,
-      hProfile->GetBinLowEdge(LastBin));
-  hProfile->Fit(fitpol1, "R");
-  fitpol1->GetParameters(parpol1);
+  //===========================================================================================================
+  //                                           Fill extended profile
+  //===========================================================================================================
 
   // fill the extended TProfile
   float newval, lastval, lasterr;
   for (int bin = 1; bin <= FirstBin; bin++) {
     newval = 0;
-    if ( i != 0 && i != 1)
-      for (int par = 0; par < 4; par++)
-        newval += parpol3[par] * pow(hProfile->GetBinCenter(bin), par);
-
-    /*
-       else if (i >= (Ncuts - 4))
-       newval = hProfile->GetBinContent(FirstBin);
-       */
-
-    else {
-      newval = 0;
-      for (int par = 0; par < 2; par++)
-        newval += parpol1[par] * pow(hProfile->GetBinCenter(bin), par);
-      pfx->SetBinContent(bin, newval);
-    }
     pfx->SetBinContent(bin, newval);
   }
+
+
   for (int bin = FirstBin; bin <= LastBin; bin++) {
     newval = hProfile->GetBinContent(bin);
     if (newval != 0) {
@@ -263,16 +257,10 @@ TSpline3* MakeSpline(TH2F* hInput) {
       pfx->SetBinError(bin, lasterr);
     }
   }
-  for (int bin = LastBin; bin <= hProfile->GetNbinsX(); bin++) {
-    newval = 0;
-    for (int par = 0; par < 2; par++)
-      newval += parpol1[par] * pow(hProfile->GetBinCenter(bin), par);
-    pfx->SetBinContent(bin, newval);
-  }
   pfx->Draw("same");
 
   gspline = new TSpline3(pfx);
-  gspline->SetName(Form("fspline_%d", i + 1));
+  gspline->SetName(Form("fspline_%d", NCallee + 1));
 
   return gspline;
 
-- 
GitLab