diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorr.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorr.C new file mode 100644 index 0000000000000000000000000000000000000000..cb5979369156add39a230845fa72dbc402ee5714 --- /dev/null +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorr.C @@ -0,0 +1,253 @@ +//#include <TICPhysics.h> +#include <TCanvas.h> +#include <TChain.h> +#include <TF1.h> +#include <TFile.h> +#include <TH2.h> +#include <TSpline.h> +#include <TCutG.h> +#include <fstream> +#include <vector> + +using namespace std; + +int GetNumberKey(TFile *infile ,const char* ClassName); + +TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double + XMax, double YMin, double YMax); +////////////////////////////////////////////////////////////////////////////////////////// +void DECorr(bool cut = false, bool spline = false) { + //=========================================================================================================== + // 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); + + TCutG *CutZ; + TCutG *CutZbis; + TFile *fCut; + if (cut) { + fCut = TFile::Open("Output/CutDeCorr.root"); + CutZbis = (TCutG*)fCut->Get("CutZbis") ; + CutZ = (TCutG*)fCut->Get("CutZbis") ; + } + + TFile *fSpline; + TSpline3 *splineDEcorr, *splineDEbis; + if(spline){ + fSpline= new TFile("Output/splineDE.root","open"); + splineDEcorr = (TSpline3*) fSpline->Get("SplineDe"); + splineDEbis = (TSpline3*) fSpline->Get("SplineDebis"); + } + + //Defining output histo + TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,100,1,1000,100,1); + TH2F* hDE_Ebis = new TH2F("DE_Ebis","DE_Ebis",1000,100,1,1000,100,1); + TH2F* hDE_Y = new TH2F("DE_Y","DE_Y",1000,-100,100,1000,0,22000); + TH2F* hDE_Y_bis = new TH2F("DE_Ybis","DE_Ybis",1000,-100,100,1000,0,22000); + TH2F* hDE_E_corr = new TH2F("DE_Ecorr","DE_Ecorr",1000,100,1,1000,0,22000); + TH2F* hDE_Ebis_corr = new TH2F("DE_Ebiscorr","DE_Ebiscorr",1000,100,1,1000,0,22000); + + + //=========================================================================================================== + // Event Loop + //=========================================================================================================== + 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); + + + double DE = 0 ,E =0 ,DE_Ycorr=0, DE_Ybis=0; + + for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ + if (seg >4){ + E += IC->fIC_raw[seg] ; + //Ecorr += ICcorr_Y.at(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]; + + DE = 0.5*(IC->fIC_raw[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; + + if (FF_IC_Y >-100 && FF_IC_Y <100){ + hDE_E->Fill(E,DE); + hDE_Ebis->Fill(E,DE_Bis); + + if (spline== true){ + DE_Ycorr = DE * 13750 / splineDEcorr->Eval(FF_IC_Y); + hDE_E_corr->Fill(E,DE_Ycorr); + + DE_Ybis = DE_Bis * 13750 / splineDEbis->Eval(FF_IC_Y); + hDE_Ebis_corr->Fill(E,DE_Ybis); + } + + if (cut == true && !(CutZ->IsInside(E,DE))) {;} + else { + hDE_Y->Fill(FF_IC_Y,DE); + }//end cut + + if (cut == true && !(CutZbis->IsInside(E,DE_Bis))) {;} + else { + hDE_Y_bis->Fill(FF_IC_Y,DE_Bis); + }//end cut + }//end cond Y + } //endl loop event + + // making spline + + TSpline3 *splineDE, *osplineDEbis; + if (cut == true){ + splineDE = MakeSpline(hDE_Y,0,-100,100,14300,15000); + osplineDEbis = MakeSpline(hDE_Y_bis,1,-100,100,14300,15000); + TFile *fSpline = new TFile("Output/splineDE.root","recreate"); + splineDE->SetName("SplineDe"); + splineDE->Write(); + osplineDEbis->SetName("SplineDebis"); + osplineDEbis->Write(); + fSpline->Close(); + fCut->cd(); + } + //=========================================================================================================== + // Drawing histo + //=========================================================================================================== + + TCanvas* c2 = new TCanvas("c2","c2",1500,1000); + c2->Divide(2 + 2*spline); + c2->cd(1); + gPad->SetLogz(); + hDE_E->Draw(); + c2->cd(2); + gPad->SetLogz(); + hDE_Ebis->Draw(); + if (spline){ + c2->cd(3); + gPad->SetLogz(); + hDE_E_corr->Draw(); + c2->cd(4); + gPad->SetLogz(); + hDE_Ebis_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; + +} + +TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double + XMax = 1000 , double YMin=0 , double YMax = 100) { + + TSpline3* gspline; + + TH1F* hProfile; + TH1F* pfx; // extended hProfile + + TCanvas* canfit = new TCanvas(Form("CanSpline_%02d",NCallee), + Form("canfit_%02d",NCallee), 0, 0, 1000, 500); + + + //=========================================================================================================== + // Init profile + //=========================================================================================================== + + // Create the TProfile + hProfile = (TH1F*)hInput->ProfileX(); + hProfile->SetLineWidth(2); + hProfile->SetDirectory(0); + canfit->cd(); + + hProfile->GetYaxis()->SetRangeUser(YMin, YMax); + hProfile->Draw(); + //=========================================================================================================== + // First and last bin to get to 6k + // event get promoted + //=========================================================================================================== + int FirstBin, LastBin; + double Treshold = 1; + for (int bin =1 ; bin<hProfile->GetNbinsX(); bin++) { + FirstBin = bin; + if (hProfile->GetBinLowEdge(bin)> XMin) break; + } + for (int bin = hProfile->GetNbinsX(); bin>1 ; bin--) { + LastBin = bin; + if (hProfile->GetBinLowEdge(bin)< XMax) break; + } + //=========================================================================================================== + // Init Extended profile function + //=========================================================================================================== + // Create the extended TProfile + pfx = new TH1F( + Form("pfx_%02d", Int_t(NCallee/2)), Form("pfx_%02d", Int_t(NCallee/2)), hProfile->GetNbinsX(), hProfile->GetBinLowEdge(1), + hProfile->GetBinLowEdge(hProfile->GetNbinsX()) + hProfile->GetBinWidth(hProfile->GetNbinsX())); + pfx->SetLineColor(8); + pfx->SetDirectory(0); + //=========================================================================================================== + // Fill extended profile + //=========================================================================================================== + + // fill the extended TProfile + float newval, lastval, lasterr; + for (int bin = 1; bin <= FirstBin; bin++) { + newval = 0; + pfx->SetBinContent(bin, newval); + } + + + for (int bin = FirstBin; bin <= LastBin; bin++) { + newval = hProfile->GetBinContent(bin); + if (newval != 0) { + pfx->SetBinContent(bin, newval); + pfx->SetBinError(bin, hProfile->GetBinError(bin)); + lastval = newval; + lasterr = hProfile->GetBinError(bin); + } + else { + pfx->SetBinContent(bin, lastval); + pfx->SetBinError(bin, lasterr); + } + } + pfx->Draw("same"); + + gspline = new TSpline3(pfx); + gspline->SetName(Form("fspline_%d", NCallee + 1)); + + return gspline; + +} // end makespline + diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorrParra.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorrParra.C new file mode 100644 index 0000000000000000000000000000000000000000..2e1be5cf48583c06b4d9f86ea1968da75cd2c05e --- /dev/null +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorrParra.C @@ -0,0 +1,253 @@ +//#include <TICPhysics.h> +#include <TCanvas.h> +#include <TChain.h> +#include <TF1.h> +#include <TFile.h> +#include <TH2.h> +#include <TSpline.h> +#include <TCutG.h> +#include <fstream> +#include <vector> + +using namespace std; + +int GetNumberKey(TFile *infile ,const char* ClassName); + +TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double + XMax, double YMin, double YMax); +////////////////////////////////////////////////////////////////////////////////////////// +void DECorrParra(bool cut = false, bool spline = false) { + //=========================================================================================================== + // 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); + + TCutG *CutZ; + TCutG *CutZbis; + TFile *fCut; + if (cut) { + fCut = TFile::Open("Output/CutDeCorr.root"); + CutZbis = (TCutG*)fCut->Get("CutZbis") ; + CutZ = (TCutG*)fCut->Get("CutZbis") ; + } + + TFile *fSpline; + TSpline3 *splineDEcorr, *splineDEbis; + if(spline){ + fSpline= new TFile("Output/splineDE.root","open"); + splineDEcorr = (TSpline3*) fSpline->Get("SplineDe"); + splineDEbis = (TSpline3*) fSpline->Get("SplineDebis"); + } + + //Defining output histo + TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,100,1,1000,100,1); + TH2F* hDE_Ebis = new TH2F("DE_Ebis","DE_Ebis",1000,100,1,1000,100,1); + TH2F* hDE_Y = new TH2F("DE_Y","DE_Y",1000,-100,100,1000,0,22000); + TH2F* hDE_Y_bis = new TH2F("DE_Ybis","DE_Ybis",1000,-100,100,1000,0,22000); + TH2F* hDE_E_corr = new TH2F("DE_Ecorr","DE_Ecorr",1000,100,1,1000,0,22000); + TH2F* hDE_Ebis_corr = new TH2F("DE_Ebiscorr","DE_Ebiscorr",1000,100,1,1000,0,22000); + + + //=========================================================================================================== + // Event Loop + //=========================================================================================================== + 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); + + + double DE = 0 ,E =0 ,DE_Ycorr=0, DE_Ybis=0; + + for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ + if (seg >4){ + E += IC->fIC_raw[seg] ; + //Ecorr += ICcorr_Y.at(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]; + + DE = 0.5*(IC->fIC_raw[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; + + if (FF_IC_Y >-100 && FF_IC_Y <100){ + hDE_E->Fill(E,DE); + hDE_Ebis->Fill(E,DE_Bis); + + if (spline== true){ + DE_Ycorr = DE * 13750 / splineDEcorr->Eval(FF_IC_Y); + hDE_E_corr->Fill(E,DE_Ycorr); + + DE_Ybis = DE_Bis * 13750 / splineDEbis->Eval(FF_IC_Y); + hDE_Ebis_corr->Fill(E,DE_Ybis); + } + + if (cut == true && !(CutZ->IsInside(E,DE))) {;} + else { + hDE_Y->Fill(FF_IC_Y,DE); + }//end cut + + if (cut == true && !(CutZbis->IsInside(E,DE_Bis))) {;} + else { + hDE_Y_bis->Fill(FF_IC_Y,DE_Bis); + }//end cut + }//end cond Y + } //endl loop event + + // making spline + + TSpline3 *splineDE, *osplineDEbis; + if (cut == true){ + splineDE = MakeSpline(hDE_Y,0,-100,100,14300,15000); + osplineDEbis = MakeSpline(hDE_Y_bis,1,-100,100,14300,15000); + TFile *fSpline = new TFile("Output/splineDE.root","recreate"); + splineDE->SetName("SplineDe"); + splineDE->Write(); + osplineDEbis->SetName("SplineDebis"); + osplineDEbis->Write(); + fSpline->Close(); + fCut->cd(); + } + //=========================================================================================================== + // Drawing histo + //=========================================================================================================== + + TCanvas* c2 = new TCanvas("c2","c2",1500,1000); + c2->Divide(2 + 2*spline); + c2->cd(1); + gPad->SetLogz(); + hDE_E->Draw(); + c2->cd(2); + gPad->SetLogz(); + hDE_Ebis->Draw(); + if (spline){ + c2->cd(3); + gPad->SetLogz(); + hDE_E_corr->Draw(); + c2->cd(4); + gPad->SetLogz(); + hDE_Ebis_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; + +} + +TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double + XMax = 1000 , double YMin=0 , double YMax = 100) { + + TSpline3* gspline; + + TH1F* hProfile; + TH1F* pfx; // extended hProfile + + TCanvas* canfit = new TCanvas(Form("CanSpline_%02d",NCallee), + Form("canfit_%02d",NCallee), 0, 0, 1000, 500); + + + //=========================================================================================================== + // Init profile + //=========================================================================================================== + + // Create the TProfile + hProfile = (TH1F*)hInput->ProfileX(); + hProfile->SetLineWidth(2); + hProfile->SetDirectory(0); + canfit->cd(); + + hProfile->GetYaxis()->SetRangeUser(YMin, YMax); + hProfile->Draw(); + //=========================================================================================================== + // First and last bin to get to 6k + // event get promoted + //=========================================================================================================== + int FirstBin, LastBin; + double Treshold = 1; + for (int bin =1 ; bin<hProfile->GetNbinsX(); bin++) { + FirstBin = bin; + if (hProfile->GetBinLowEdge(bin)> XMin) break; + } + for (int bin = hProfile->GetNbinsX(); bin>1 ; bin--) { + LastBin = bin; + if (hProfile->GetBinLowEdge(bin)< XMax) break; + } + //=========================================================================================================== + // Init Extended profile function + //=========================================================================================================== + // Create the extended TProfile + pfx = new TH1F( + Form("pfx_%02d", Int_t(NCallee/2)), Form("pfx_%02d", Int_t(NCallee/2)), hProfile->GetNbinsX(), hProfile->GetBinLowEdge(1), + hProfile->GetBinLowEdge(hProfile->GetNbinsX()) + hProfile->GetBinWidth(hProfile->GetNbinsX())); + pfx->SetLineColor(8); + pfx->SetDirectory(0); + //=========================================================================================================== + // Fill extended profile + //=========================================================================================================== + + // fill the extended TProfile + float newval, lastval, lasterr; + for (int bin = 1; bin <= FirstBin; bin++) { + newval = 0; + pfx->SetBinContent(bin, newval); + } + + + for (int bin = FirstBin; bin <= LastBin; bin++) { + newval = hProfile->GetBinContent(bin); + if (newval != 0) { + pfx->SetBinContent(bin, newval); + pfx->SetBinError(bin, hProfile->GetBinError(bin)); + lastval = newval; + lasterr = hProfile->GetBinError(bin); + } + else { + pfx->SetBinContent(bin, lastval); + pfx->SetBinError(bin, lasterr); + } + } + pfx->Draw("same"); + + gspline = new TSpline3(pfx); + gspline->SetName(Form("fspline_%d", NCallee + 1)); + + return gspline; + +} // end makespline + diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C index 7497d41b1a9ae4f73cd95da1a2734861866bfc37..015ccbf949f40a43ff7d629ba022ad1c8fdb0467 100644 --- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C @@ -1,4 +1,4 @@ -#include <TICPhysics.h> +//#include <TICPhysics.h> #include <TCanvas.h> #include <TChain.h> #include <TF1.h> @@ -54,13 +54,13 @@ void ICCorr() { hICcorr.push_back(new TH2F(Form("hICorr%d_Y",seg),Form("hICorr%d_Y",seg),1000,-100,100,500,2000,9500)); } else if (seg == 0){ - hICX.push_back((TH2F*)fHisto->Get(Form("hIC%d_IC%d_X",seg+1,seg))); - hICY.push_back((TH2F*)fHisto->Get(Form("hIC%d_IC%d_Y",seg+1,seg))); + hICX.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_IC%d_X",seg+1,seg))); + hICY.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_IC%d_Y",seg+1,seg))); hICcorr.push_back(new TH2F(Form("hICorr%d_Y",seg),Form("hICorr%d_Y",seg),1000,-100,100,500,0,2)); } else { - hICX.push_back((TH2F*)fHisto->Get(Form("hIC%d_IC%d_X",seg,seg-1))); - hICY.push_back((TH2F*)fHisto->Get(Form("hIC%d_IC%d_Y",seg,seg-1))); + hICX.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_IC%d_X",seg,seg-1))); + hICY.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_IC%d_Y",seg,seg-1))); hICcorr.push_back(new TH2F(Form("hICorr%d_Y",seg),Form("hICorr%d_Y",seg),1000,-100,100,500,0,2)); } } @@ -116,108 +116,119 @@ void ICCorr() { //=========================================================================================================== // Event Loop //=========================================================================================================== - int Nentries = chain->GetEntries(); - + int Nentries = chain->GetEntries(); + auto start = std::chrono::high_resolution_clock::now(); 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); + + 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: " << timeLeft << " seconds" << "\r" << flush; + } + + chain->GetEntry(e); - vector<double> ICcorr_Y(5), ICcorr_X(11); // the [0] element of spline is the + vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the // correction on the first // segment of ic for (int seg = 0; seg < ICcorr_Y.size() ; seg++) { - if (seg == 1){ - ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(0)->Eval(0)/spline_Y.at(0)->Eval(FF_IC_Y); - } - - if (seg == 0) { - double temp = ICcorr_Y.at(1) / IC->fIC_raw[0] * spline_Y.at(1)->Eval(0)/ spline_Y.at(1)->Eval(FF_IC_Y); - ICcorr_Y.at(seg) = ICcorr_Y.at(1) / temp; - //cout << ICcorr_Y.at(0) << " " << IC->fIC_raw[0]<< endl; - } - else if (seg > 1) { - double temp = 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; - } - if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0; - } //End loop on histogram - - - for (int i = 0; i < ICcorr_Y.size() ; i++) { - hICcorr.at(i)->Fill(FF_IC_Y,ICcorr_Y.at(i)); - } //End loop on histogram - - double DE = 0 ,E =0 ,DE_IC1_corrY =0, DE_ICALL_corrY=0, Ecorr=0; - - for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ - if (seg >4){ - E += IC->fIC_raw[seg] ; - //Ecorr += ICcorr_Y.at(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]; - - DE = 0.5*(IC->fIC_raw[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; - DE_IC1_corrY = 0.5*(ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; - DE_ICALL_corrY = 0.5*(ICcorr_Y[0]+ICcorr_Y[1] +ICcorr_Y[2]+ICcorr_Y[3]); - double Ealt = E + IC->fIC_raw[4]; - - if (FF_IC_Y >-60 && FF_IC_Y <60){ - hDEBis_Ecorr_IC1->Fill(E,DE_Bis); - hDE_E->Fill(E,DE); - hDE_Ecorr1->Fill(E,DE_IC1_corrY); - hDE_E_CorrAll->Fill(E,DE_ICALL_corrY); - } - } //endl loop event - //=========================================================================================================== - // Drawing histo - //=========================================================================================================== - - TCanvas* c1 = new TCanvas("c1","c1",10000,1000); - c1->Divide(11,2); - for (int i = 0 ; i<11 ;i++){ - c1->cd(i+1); - hICcorr.at(i)->ProfileX()->Draw(); - c1->cd(11+i+1); - hICY.at(i)->ProfileX()->Draw(); - } - - TCanvas* c2 = new TCanvas("c2","c2",1500,1000); - c2->Divide(4); - c2->cd(1); - gPad->SetLogz(); - hDE_E->Draw(); - c2->cd(2); - gPad->SetLogz(); - hDE_Ecorr1->Draw(); - c2->cd(3); - gPad->SetLogz(); - hDEBis_Ecorr_IC1->Draw(); - c2->cd(4); - gPad->SetLogz(); - hDE_E_CorrAll->Draw(); + if (spline_Y.at(seg)==0){ + ICcorr_Y.at(seg) = IC->fIC_raw[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); + } + + if (seg == 0) { + double temp = ICcorr_Y.at(1) / IC->fIC_raw[0] * spline_Y.at(0)->Eval(0)/ spline_Y.at(0)->Eval(FF_IC_Y); + ICcorr_Y.at(seg) = ICcorr_Y.at(1) / temp; + //cout << ICcorr_Y.at(0) << " " << IC->fIC_raw[0]<< endl; + } + else if (seg > 1) { + double temp = 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; + } + if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0; + } //End loop on histogram + } + + for (int i = 0; i < ICcorr_Y.size() ; i++) { + hICcorr.at(i)->Fill(FF_IC_Y,ICcorr_Y.at(i)); + } //End loop on histogram + + double DE = 0 ,E =0 ,DE_IC1_corrY =0, DE_ICALL_corrY=0, Ecorr=0; + + for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ + if (seg >4){ + E += IC->fIC_raw[seg] ; + //Ecorr += ICcorr_Y.at(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]; + + DE = 0.5*(IC->fIC_raw[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; + DE_IC1_corrY = 0.5*(ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; + DE_ICALL_corrY = 0.5*(ICcorr_Y[0]+ICcorr_Y[1] +ICcorr_Y[2]+ICcorr_Y[3])+IC->fIC_raw[4]; + double Ealt = E + IC->fIC_raw[4]; + + if (FF_IC_Y >-60 && FF_IC_Y <60){ + hDEBis_Ecorr_IC1->Fill(E,DE_Bis); + hDE_E->Fill(E,DE); + hDE_Ecorr1->Fill(E,DE_IC1_corrY); + hDE_E_CorrAll->Fill(E,DE_ICALL_corrY); + } + } //endl loop event + //=========================================================================================================== + // Drawing histo + //=========================================================================================================== + /* + TCanvas* c1 = new TCanvas("c1","c1",10000,1000); + c1->Divide(NSegment,2); + for (int i = 0 ; i<NSegment ;i++){ + c1->cd(i+1); + hICcorr.at(i)->ProfileX()->Draw(); + + c1->cd(NSegment+i+1); + hICY.at(i)->ProfileX()->Draw(); + } + */ + TCanvas* c2 = new TCanvas("c2","c2",1500,1000); + c2->Divide(4); + c2->cd(1); + gPad->SetLogz(); + hDE_E->Draw(); + c2->cd(2); + gPad->SetLogz(); + hDE_Ecorr1->Draw(); + c2->cd(3); + gPad->SetLogz(); + hDEBis_Ecorr_IC1->Draw(); + c2->cd(4); + gPad->SetLogz(); + hDE_E_CorrAll->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; + // 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.h b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.h index 2a2b364371d6e9393a586008c6eb68f4ea314afa..e4c7d79537f4062f13f402a2740f9327b86f6b7d 100644 --- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.h +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.h @@ -1,5 +1,5 @@ #include <TFile.h> -#include <TICPhysics.h> +//#include <TICPhysics.h> #include <TChain.h> #include <TF1.h> #include <TH2.h> diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ReadMe.md b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ReadMe.md index 74440f95c2f9bd869b0585ed48bedfdc95d196f0..1541c4dd52def5128368bc7f8696bc783bf7030e 100644 --- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ReadMe.md +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ReadMe.md @@ -15,3 +15,7 @@ of the electron by the grid. SplineXY.C serve the purpose of creating this spline and storing int a .root file. + +## Devel + +Apparently the main correction must come from the correction of the DE diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C index 6cb779075b87e972ce33f67beab260de85760d68..e906d72eadbb7532609a01b94ea7989c8c7a0fb5 100644 --- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C @@ -1,6 +1,6 @@ -//#ifdef TICPhysics +#ifdef TICPhysics #include <TICPhysics.h> -//#endif +#endif #include <TCanvas.h> #include <TChain.h> #include <TF1.h> @@ -39,21 +39,21 @@ void SplineChioP2P() { int seg = 1; vector<vector<TSpline3*>> SplineAllIC(11); - vector<TSpline3*> SplineICurrent; + vector<TSpline3*> SplineICurrent(2); vector<vector<TH2F*>> hIC(NSegment); hIC.at(1)= (HistoFillerIC(1,SplineAllIC)); - SplineICurrent.push_back(MakeSpline(hIC[seg][0] , NCallSpline, XMin[0] , XMax[0], YMin[0], YMax[0])); + SplineICurrent.at(0) = MakeSpline(hIC[seg][0] , NCallSpline, XMin[0] , XMax[0], YMin[0], YMax[0]); SplineICurrent.at(0)->SetName("fsplineIC1_X"); NCallSpline+=1; - - SplineICurrent.push_back(MakeSpline(hIC[seg][1] , NCallSpline, XMin[1] , XMax[1], YMin[1], YMax[1])); + SplineICurrent.at(1) = MakeSpline(hIC[seg][1] , NCallSpline, XMin[1] , XMax[1], YMin[1], YMax[1]); SplineICurrent.at(1)->SetName("fsplineIC1_Y"); NCallSpline+=1; SplineAllIC.at(1)= (SplineICurrent); // + //SplineAllIC.erase(SplineAllIC.begin() + 1); // for (int seg = 0 ; seg < NSegment ; seg++){ @@ -87,16 +87,17 @@ void SplineChioP2P() { //write spline TFile* fspline = new TFile("Output/spline_P2P_2024.root", "recreate"); for (int seg = 0 ; seg < NSegment ; seg++){ - SplineAllIC.at(seg)[0]->Write(); - SplineAllIC.at(seg)[1]->Write(); + if (!SplineAllIC.at(seg).empty()){ + SplineAllIC.at(seg)[0]->Write(); + SplineAllIC.at(seg)[1]->Write(); + } } fspline->Close(); - - + + vector<vector<TH2F*>> hIC_Corrall(NSegment); vector<TH1F*> hIC_Corr_profile(NSegment); vector<TH2F*> hIC_Corr(NSegment); - hIC_Corrall = HistoFillerICcorr(NSegment,SplineAllIC); hIC_Corr = hIC_Corrall.at(1); @@ -177,8 +178,8 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC // Beginning loop on entries //=========================================================================================================== - int Nentries = chain->GetEntries(); - //int Nentries = 100000; + //int Nentries = chain->GetEntries(); + int Nentries = 1000000; auto start = std::chrono::high_resolution_clock::now(); for (int e = 0; e < Nentries; e++) { @@ -318,8 +319,8 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp // Beginning loop on entries //=========================================================================================================== - int Nentries = chain->GetEntries(); - //int Nentries = 100000; + //int Nentries = chain->GetEntries(); + int Nentries = 1000000; auto start = std::chrono::high_resolution_clock::now(); for (int e = 0; e < Nentries; e++) { @@ -333,6 +334,7 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp if(Spline_All_ICprev.at(pseg).empty()){ IC_Spline_CorrX.at(pseg)= IC->fIC_raw[pseg]; IC_Spline_CorrY.at(pseg)= IC->fIC_raw[pseg]; + //cout << pseg << endl; } //If it exist correct the current line @@ -347,7 +349,7 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp } //end ic 2 to 10 } //End if non empty } // end loop - + //filling 0 if (!Spline_All_ICprev.at(0).empty()){ IC_Spline_CorrX.at(0) = IC->fIC_raw[0] / SplineIC_X.at(0)->Eval(0) * SplineIC_X.at(0)->Eval(FF_IC_X);