diff --git a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/ReadMe.md b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/ReadMe.md
new file mode 100644
index 0000000000000000000000000000000000000000..33bb0e4de760864f6193d2d8abe39236fb2b99c5
--- /dev/null
+++ b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/ReadMe.md
@@ -0,0 +1,5 @@
+# Requirement
+
+You need to fill Cut folder with the cut of each charge
+You can do this simply by using the DECorr.C macro in the YCalibration folder
+The name of the cut should be Z1->ZN and you should space them a little
diff --git a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
index cdf4136a2d57b976b5c983118e742b7a00162b3f..eb7c1011f258abdd4000e661ba873e344470ab96 100644
--- a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
+++ b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
@@ -1,196 +1,274 @@
+#include "TICPhysics.h"
+#include <TChain.h>
+#include <TCutG.h>
+#include <TF1.h>
+#include <TFile.h>
+#include <TH2.h>
+#include <TSpline.h>
#include <fstream>
using namespace std;
-const int Ncuts=34;
+const int Ncuts=28;
+const int NSegment=11;
const int Nrun=1;
-int RunNumber[Nrun]={48};
+int RunNumber[Nrun]={241};
void MakeSpline();
void ApplySpline_perTCutG();
void ApplySpline();
+vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline);
///////////////////////////////////////////////////////////////////////////
void SplineChio()
{
+ //===========================================================================================================
+ // Loading Var
+ //===========================================================================================================
+
+ //**********************Cut**************************************
+ TFile *fcut=new TFile("Cut/Cut_Z.root","open");
+ TCutG *cutZ[Ncuts];
+ for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Z%i",i+1));
+
+ // ****************** Spline Y **********************************
+ TFile *fSplineIC = new TFile("../YCalibration/Output/spline_P2P_2024.root","open");
+
+
+ vector<vector<TSpline3*>> Spline = LoadSpline("../YCalibration/Output/spline_P2P_2024.root");
+ vector<TSpline3*> spline_X(11), spline_Y(11);
+ spline_X = Spline.at(0) ; spline_Y = Spline.at(1);
+
+ //**********************DE_E****************************************************
+ TChain *chain=new TChain("PhysicsTree");
+
+ chain->Add("../../../root/analysis/VamosCalib241.root");
+
+ TICPhysics* IC = new TICPhysics() ;
+ double FF_IC_X, FF_IC_Y, FF_V13;
+ double FF_DE, FF_Eres;
+
+ chain->SetBranchStatus("IC", true);
+ chain->SetBranchAddress("IC", &IC);
+ chain->SetBranchStatus("FF_IC_Y", true);
+ chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y);
+
+
+ TH2F *hChioDE_E_all = new TH2F("hChioDE_E","hChioDE_E",1000,1000,20000,500,1000,26000);
+ 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,1000,20000,1250,1000,26000);
+
+ auto start = std::chrono::high_resolution_clock::now();
+ int Nentries=chain->GetEntries();
+ 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);
+
+ vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the
+ vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first
+ // segment of ic
+ for (int seg = 1; seg < NSegment+1 ; seg++) {
+
+ if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0
+
+ if (spline_Y.at(seg)==0){
+ ICcorr_Y.at(seg) = IC->fIC_raw[seg];
+ }
+
+ else {
+ if (seg == 0) {
+ Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg);
+ }
+
+ else if (seg == 1){
+ ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y);
+ }
+
+ else if (seg > 1) {
+ Temp_Y.at(seg) = IC->fIC_raw[seg]/ICcorr_Y.at(seg-1) * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(seg-1) * Temp_Y.at(seg);
+ }
+
+ if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
+ } //end if non empty
+ if (seg == 0) break;
+ }//endloop seg
+
- TFile *fcut=new TFile("cut_chio.root");
- TCutG *cutZ[Ncuts];
- for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("cut%i",i+31));
-
- double FF_DE, FF_Eres;
- TChain *chain=new TChain("PhysicsTree");
- for(int i=0;i<Nrun;i++)
- {
- chain->Add(Form("../../root/analysis/run_%d.root",RunNumber[i]));
- cout << "RunNumber=" << RunNumber[i] << endl;
- }
- chain->SetBranchStatus("FF_DE",true);
- chain->SetBranchAddress("FF_DE",&FF_DE);
- chain->SetBranchStatus("FF_Eres",true);
- chain->SetBranchAddress("FF_Eres",&FF_Eres);
-
- TH2F *hChioDE_E_all = new TH2F("hChioDE_E","hChioDE_E",1000,1000,20000,500,1000,26000);
- 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,1000,20000,1250,1000,26000);
-
- int Nentries=chain->GetEntries();
- for(int e=0;e<Nentries;++e)
- {
- if(e%100000==0) cout << e*100./Nentries << "% \r" << flush;
-
- chain->GetEntry(e);
- hChioDE_E_all->Fill(FF_Eres,FF_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(FF_Eres,FF_DE)) {hChioDE_E[i]->Fill(FF_Eres,FF_DE); break;}
- }
-
- 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(); gPad-> SetLogz();
- hChioDE_E_all->Draw("colz");
- for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
+ FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+ FF_Eres = 0 ;
+
+ for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+ if (seg >4){
+ FF_Eres += IC->fIC_raw[seg] ;
+ //Ecorr += ICcorr_Y.at(seg) ;
+ }
+ }
+
+
+ hChioDE_E_all->Fill(FF_Eres,FF_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(FF_Eres,FF_DE)) {hChioDE_E[i]->Fill(FF_Eres,FF_DE); break;}
+ }
+
+
+ //**********************************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(); gPad-> SetLogz();
+ hChioDE_E_all->Draw("colz");
+ for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
}
///////////////////////////////////////////////////////////////////////////
void MakeSpline()
{
- TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root");
-
- TSpline3 *gspline[Ncuts];
-
- TH2F* h2 [Ncuts];
- TH1F* hProfile[Ncuts];
- TH1F* pfx[Ncuts]; // extended hProfile
- TFile *fspline=new TFile("spline_Chio_2023.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
- 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;
- }
+ TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root");
+
+ TSpline3 *gspline[Ncuts];
+
+ TH2F* h2 [Ncuts];
+ TH1F* hProfile[Ncuts];
+ 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
+ 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)+10000);
+ 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();
- // 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)+10000);
- 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();
- fspline->Close();
+ TGraph * gr = new TGraph();
+ gr->SetName("DEvsZ");
- 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 = "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;
- for(int i=0;i<Ncuts;i++){
+ TCanvas * can = new TCanvas("ChioEcorr_vs_Ebis_for_TSplines","ChioEcorr_vs_Ebis_for_TSplines",0,0,2500,2000);
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(8500),i+31);
- //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("pol2");
+ 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;
+ 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(8500),i+31);
+ //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("pol2");
}
@@ -198,331 +276,501 @@ void MakeSpline()
void ApplySpline_perTCutG()
{
- // TCutG from Lucas
- TFile *fcut=new TFile("cut_chio.root");
- TCutG *cutZ[Ncuts];
-
- // TSpline
- TFile *fspline=new TFile("spline_Chio_2023.root","read");
- TSpline3 *gspline[Ncuts];
-
- for(int i=0;i<Ncuts;i++){
- cutZ[i]=(TCutG*)fcut->Get(Form("cut%i",i+31));
- gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
- }
-
- // input data
- TChain *chain=new TChain("PhysicsTree");
- for(int i=0;i<Nrun;i++)
- {
- chain->Add(Form("../root/analysis/run_%d.root",RunNumber[i]));
- cout << "RunNumber=" << RunNumber[i] << endl;
- }
- double FF_DE, FF_Eres;
- chain->SetBranchStatus("FF_DE",true);
- chain->SetBranchAddress("FF_DE",&FF_DE);
- chain->SetBranchStatus("FF_Eres",true);
- chain->SetBranchAddress("FF_Eres",&FF_Eres);
-
- // variable after applying TSpline
- double FF_DEcorr;
- double FF_DEcorr0[Ncuts] = {
- 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();
-
- for(int e=0;e<Nentries;++e)
- {
- if(e%100000==0) cout << e*100./Nentries << "% \r" << flush;
-
- chain->GetEntry(e);
- hChioDE_E_all->Fill(FF_Eres,FF_DE);
- for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) {
- hChioDE_E->Fill(FF_Eres,FF_DE);
- FF_DEcorr = FF_DEcorr0[i] * FF_DE / gspline[i]->Eval(FF_Eres);
- hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr);
- break;
+ // *************************Cut***************************
+ TFile *fcut=new TFile("Cut/Cut_Z.root");
+ TCutG *cutZ[Ncuts];
+
+ // ******************* TSpline DE*************************
+ TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
+ TSpline3 *gspline[Ncuts];
+
+ for(int i=0;i<Ncuts;i++){
+ cutZ[i]=(TCutG*)fcut->Get(Form("Z%i",i+1));
+ gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
+ }
+
+ // ****************** Spline Y **********************************
+ TFile *fSplineIC = new TFile("../YCalibration/Output/spline_P2P_2024.root","open");
+
+
+ vector<vector<TSpline3*>> Spline = LoadSpline("../YCalibration/Output/spline_P2P_2024.root");
+ vector<TSpline3*> spline_X(11), spline_Y(11);
+ spline_X = Spline.at(0) ; spline_Y = Spline.at(1);
+
+ //**********************DE_E****************************************************
+ TChain *chain=new TChain("PhysicsTree");
+
+ chain->Add("../../../root/analysis/VamosCalib241.root");
+
+ TICPhysics* IC = new TICPhysics() ;
+ double FF_IC_X, FF_IC_Y, FF_V13;
+ double FF_DE, FF_Eres;
+
+ chain->SetBranchStatus("IC", true);
+ chain->SetBranchAddress("IC", &IC);
+ chain->SetBranchStatus("FF_IC_Y", true);
+ chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y);
+
+ // variable after applying TSpline
+ double FF_DEcorr;
+ double FF_DEcorr0[Ncuts] = {
+ 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();
+ 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);
+
+
+
+ vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the
+ vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first
+ // segment of ic
+ for (int seg = 1; seg < NSegment+1 ; seg++) {
+
+ if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0
+
+ if (spline_Y.at(seg)==0){
+ ICcorr_Y.at(seg) = IC->fIC_raw[seg];
+ }
+
+ else {
+ if (seg == 0) {
+ Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg);
+ }
+
+ else if (seg == 1){
+ ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y);
+ }
+
+ else if (seg > 1) {
+ Temp_Y.at(seg) = IC->fIC_raw[seg]/ICcorr_Y.at(seg-1) * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(seg-1) * Temp_Y.at(seg);
+ }
+
+ if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
+ } //end if non empty
+ if (seg == 0) break;
+ }//endloop seg
+
+
+ FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+ FF_Eres = 0 ;
+
+ for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+ if (seg >4){
+ FF_Eres += IC->fIC_raw[seg] ;
+ //Ecorr += ICcorr_Y.at(seg) ;
+ }
+ }
+
+ hChioDE_E_all->Fill(FF_Eres,FF_DE);
+ for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) {
+ hChioDE_E->Fill(FF_Eres,FF_DE);
+ FF_DEcorr = FF_DEcorr0[i] * FF_DE / gspline[i]->Eval(FF_Eres);
+ hChioDE_E_corr->Fill(FF_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");
+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 ApplySpline()
{
- // TCutG from Lucas
- TFile *fcut=new TFile("cut_chio.root");
- TCutG *cutZ[Ncuts];
-
- // TSpline
- TFile *fspline=new TFile("spline_Chio_2023.root","read");
- TSpline3 *gspline[Ncuts];
-
- for(int i=0;i<Ncuts;i++){
- cutZ[i]=(TCutG*)fcut->Get(Form("cut%d",i+31));
- gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
- }
-
- // input data
- TChain *chain=new TChain("PhysicsTree");
- for(int i=0;i<Nrun;i++)
- {
- chain->Add(Form("../../root/analysis/run_%d.root",RunNumber[i]));
- cout << "RunNumber=" << RunNumber[i] << endl;
- }
- double FF_DE, FF_Eres;
- chain->SetBranchStatus("FF_DE",true);
- chain->SetBranchAddress("FF_DE",&FF_DE);
- chain->SetBranchStatus("FF_Eres",true);
- chain->SetBranchAddress("FF_Eres",&FF_Eres);
-
- // variable after applying TSpline
- double FF_DEcorr=0;
- double Zrough=0;
- double Zabitbetter=0;
-
- double FF_DEcorr0[Ncuts];
- for(int index=0; index<Ncuts; index++){
- FF_DEcorr0[index] = gspline[index]->Eval(8500);
- }
- /*double FF_DEcorr0[Ncuts] = {
- 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};*/
-
- // 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;//chain->GetEntries();
-
- double FF_Eres_prev = 0;
- for(int e=0;e<Nentries;e++)
- {
- if(e%100000==0) cout << e*100./Nentries << "% \r" << flush;
+ // *************************Cut***************************
+ TFile *fcut=new TFile("Cut/Cut_Z.root");
+ TCutG *cutZ[Ncuts];
-
- FF_DEcorr = -100;
- chain->GetEntry(e);
- if (FF_DE<3000) continue;
- if(FF_Eres==FF_Eres_prev) continue;
- FF_Eres_prev = FF_Eres;
-
- hChioDE_E_all->Fill(FF_Eres,FF_DE);
-
- float Eval_DEspline, DEspline0;
- int index=0;
- for(int i=0; i<Ncuts; i++){
- Eval_DEspline = gspline[i]->Eval(FF_Eres);
- if(FF_DE<Eval_DEspline) break;
- index = i;
- }
+ // ******************* TSpline DE*************************
+ TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
+ TSpline3 *gspline[Ncuts];
- Eval_DEspline = gspline[index]->Eval(FF_Eres);
- DEspline0 = FF_DEcorr0[index];
- float dmin, dsup;
- if( index < (Ncuts-1) && FF_DE > gspline[0]->Eval(FF_Eres)){
- dmin = FF_DE-gspline[index]->Eval(FF_Eres);
- dsup = gspline[index+1]->Eval(FF_Eres)-FF_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(FF_Eres) ;
- // DEspline0 = FF_DEcorr0[index] ;
- //}
- Eval_DEspline = dsup*gspline[index]->Eval(FF_Eres)/(dmin+dsup) + dmin*gspline[index+1]->Eval(FF_Eres)/(dmin+dsup) ;
- DEspline0 = dsup*FF_DEcorr0[index]/(dmin+dsup) + dmin*FF_DEcorr0[index+1]/(dmin+dsup) ;
-
- FF_DEcorr = DEspline0 * FF_DE / Eval_DEspline ;
- hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr);
- //if(FF_DEcorr>15120 && FF_DEcorr<15130) {
- // cout << e << " " << index << " " << FF_Eres << " " << FF_DEcorr << endl;
- //}
+ for(int i=0;i<Ncuts;i++){
+ cutZ[i]=(TCutG*)fcut->Get(Form("Z%i",i+1));
+ gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
}
- // 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(FF_Eres,Zrough);
+ // ****************** Spline Y **********************************
+ TFile *fSplineIC = new TFile("../YCalibration/Output/spline_P2P_2024.root","open");
+
+
+ vector<vector<TSpline3*>> Spline = LoadSpline("../YCalibration/Output/spline_P2P_2024.root");
+ vector<TSpline3*> spline_X(11), spline_Y(11);
+ spline_X = Spline.at(0) ; spline_Y = Spline.at(1);
+
+ //**********************DE_E****************************************************
+ TChain *chain=new TChain("PhysicsTree");
+
+ chain->Add("../../../root/analysis/VamosCalib241.root");
+
+ TICPhysics* IC = new TICPhysics() ;
+ double FF_IC_X, FF_IC_Y, FF_V13;
+ double FF_DE, FF_Eres;
+
+ chain->SetBranchStatus("IC", true);
+ chain->SetBranchAddress("IC", &IC);
+ chain->SetBranchStatus("FF_IC_Y", true);
+ chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y);
+
+
+ // variable after applying TSpline
+ double FF_DEcorr=0;
+ double Zrough=0;
+ double Zabitbetter=0;
- if(3000<FF_Eres && FF_Eres<25000){
- hChioDE_corr->Fill(FF_DEcorr);
- hChioZ_rough->Fill(Zrough);
- hChioZ_abitbetter->Fill(Zabitbetter);
+ double FF_DEcorr0[Ncuts];
+ for(int index=0; index<Ncuts; index++){
+ FF_DEcorr0[index] = gspline[index]->Eval(8500);
}
- }
-
- 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("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();
+ /*double FF_DEcorr0[Ncuts] = {
+ 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};*/
+
+ // 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;//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);
+
+ //===========================================================================================================
+ // Splined Y
+ //===========================================================================================================
+ vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the
+ vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first
+ // segment of ic
+ for (int seg = 1; seg < NSegment+1 ; seg++) {
+
+ if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0
+
+ if (spline_Y.at(seg)==0){
+ ICcorr_Y.at(seg) = IC->fIC_raw[seg];
+ }
+
+ else {
+ if (seg == 0) {
+ Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg);
+ }
+
+ else if (seg == 1){
+ ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y);
+ }
+
+ else if (seg > 1) {
+ Temp_Y.at(seg) = IC->fIC_raw[seg]/ICcorr_Y.at(seg-1) * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(seg-1) * Temp_Y.at(seg);
+ }
+
+ if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
+ } //end if non empty
+ if (seg == 0) break;
+ }//endloop seg
+
+
+ FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+ FF_Eres = 0 ;
+
+ for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+ if (seg >4){
+ FF_Eres += IC->fIC_raw[seg] ;
+ //Ecorr += ICcorr_Y.at(seg) ;
+ }
+ }
+
+ //===========================================================================================================
+ // Spline DE
+ //===========================================================================================================&
+
+
+ if (FF_DE<3000) continue;
+ if(FF_Eres==FF_Eres_prev) continue;
+ FF_Eres_prev = FF_Eres;
+
+ hChioDE_E_all->Fill(FF_Eres,FF_DE);
+
+ float Eval_DEspline, DEspline0;
+ int index=0;
+ for(int i=0; i<Ncuts; i++){
+ Eval_DEspline = gspline[i]->Eval(FF_Eres);
+ if(FF_DE<Eval_DEspline) break;
+ index = i;
+ }
+
+ Eval_DEspline = gspline[index]->Eval(FF_Eres);
+ DEspline0 = FF_DEcorr0[index];
+ float dmin, dsup;
+ if( index < (Ncuts-1) && FF_DE > gspline[0]->Eval(FF_Eres)){
+ dmin = FF_DE-gspline[index]->Eval(FF_Eres);
+ dsup = gspline[index+1]->Eval(FF_Eres)-FF_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(FF_Eres) ;
+ // DEspline0 = FF_DEcorr0[index] ;
+ //}
+ Eval_DEspline = dsup*gspline[index]->Eval(FF_Eres)/(dmin+dsup) + dmin*gspline[index+1]->Eval(FF_Eres)/(dmin+dsup) ;
+ DEspline0 = dsup*FF_DEcorr0[index]/(dmin+dsup) + dmin*FF_DEcorr0[index+1]/(dmin+dsup) ;
+
+ FF_DEcorr = DEspline0 * FF_DE / Eval_DEspline ;
+ hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr);
+ //if(FF_DEcorr>15120 && FF_DEcorr<15130) {
+ // cout << e << " " << index << " " << FF_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(FF_Eres,Zrough);
+
+ if(3000<FF_Eres && FF_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 Fit_DE_E_corr()
{
- TFile * f = new TFile(Form("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");
+ 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");
+
+}
+
+
+vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline){
+ TFile *fSplineIC = new TFile(PathToSpline,"open");
+ // Get number of spline
+ int SplineCount = 0 ;
+ TIter next(fSplineIC->GetListOfKeys());
+ TKey* key;
+
+ while ((key=(TKey*)next())){
+ if (std::string(key->GetClassName()) == "TSpline3"){
+ SplineCount ++;
+ }
+ }
+ vector<vector<TSpline3*>> Spline;
+ vector<TSpline3*> spline_X(11), spline_Y(11);
+ for (int i = 0; i < SplineCount ; i++) {
+ int seg = int((i-2)/2 +1 );
+ if (i < 2){
+ spline_X.at(0) = (TSpline3 *)fSplineIC->Get("fspline_IC0_X");
+ spline_Y.at(0) = (TSpline3 *)fSplineIC->Get("fspline_IC0_Y");
+ }
+
+ else if ( i>=2 && i<4){
+ spline_X.at(1) = (TSpline3 *)fSplineIC->Get("fspline_IC1_X");
+ spline_Y.at(1) = (TSpline3 *)fSplineIC->Get("fspline_IC1_Y");
+ }
+ else if (seg >=2 && (i%2 == 0) ) {
+ spline_X.at(seg) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_X",seg));
+ }
+ else if (seg >=2 && !(i%2 == 0) ) {
+ spline_Y.at(seg) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_Y",seg));
+ }
+ } //End loop on histogram
+
+ Spline.push_back(spline_X);
+ Spline.push_back(spline_Y);
+
+ return Spline;
}