From b7c9ccec3306ae7d1964061dbb43fceca6f29ab0 Mon Sep 17 00:00:00 2001
From: "theodore.efremov" <theodore.efremov@cea.fr>
Date: Thu, 7 Nov 2024 16:07:36 +0100
Subject: [PATCH] P2P chio correction
---
.../2024/macro/chio/CalibrationChio/ICCorr.C | 338 +++++++++---------
.../chio/CalibrationChio/SplineChioP2P.C | 294 ++++++++++++---
2 files changed, 424 insertions(+), 208 deletions(-)
diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
index 660cd4177..6aa528b6f 100644
--- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
@@ -30,8 +30,8 @@ void ICCorr() {
chain->SetBranchAddress("IC", &IC);
- TFile* fHisto = TFile::Open("Output/SplineInputAll.root");
- TFile* fSpline = TFile::Open("Output/spline_ICALL_2024.root");
+ TFile* fHisto = TFile::Open("Output/HistoP2P.root");
+ TFile* fSpline = TFile::Open("Output/spline_P2P_2024.root");
TFile* fFit = TFile::Open("Output/FitIC.root");
//Defining output histo
@@ -45,172 +45,182 @@ void ICCorr() {
// Filling the input histo
int NSegment = 11 ;
- vector<TH2F*> hIC1_ICn, hICX,hICY ,hIC1_ICnX,hIC1_ICnY, hICcorr;
+ vector<TH2F*> hICseg_ICprev, hICX,hICY ,hICseg_ICprevX,hICseg_ICprevY, hICcorr;
for(int seg=0; seg<NSegment; seg++){
- hICX.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_X",seg)));
- hICY.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_Y",seg)));
-
- if (seg < NSegment-1){
- hIC1_ICn.push_back((TH2F*)fHisto->Get(Form("hIC1_IC%d_X",seg)));
- hIC1_ICn.push_back((TH2F*)fHisto->Get(Form("hIC1_IC%d_Y",seg)));
- hIC1_ICnX.push_back((TH2F*)fHisto->Get(Form("hIC1_IC%d_X",seg)));
- hIC1_ICnY.push_back((TH2F*)fHisto->Get(Form("hIC1_IC%d_Y",seg)));
- }
- hICcorr.push_back(new TH2F(Form("hICorr%d_Y",seg),Form("hICorr%d_Y",seg),1000,-100,100,500,2000,9500));
- }
-
- //===========================================================================================================
- // 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 ++;
- }
- }
-
-
- vector<TSpline3*> spline_X, spline_Y;
- for (int i = 0; i < SplineCount ; i++) {
- if (i == 0){
- spline_X.push_back((TSpline3 *)fSpline->Get("fsplineIC1_X"));
- }
-
- else if (i == 1){
- spline_Y.push_back((TSpline3 *)fSpline->Get("fsplineIC1_Y"));
- }
- else if (i>=2 && (i%2 == 0) && (i!=4) && false) {
- spline_X.push_back((TSpline3 *)fSpline->Get(Form("fsplineIC1_IC%d_X",int((i-2)/2))));
- }
- else if (i>=2 && !(i%2 == 0) && (i!=5)) {
- spline_Y.push_back((TSpline3 *)fSpline->Get(Form("fsplineIC1_IC%d_Y",int((i-2)/2))));
- }
- } //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
-
-
- //===========================================================================================================
- // Event Loop
- //===========================================================================================================
- int Nentries = chain->GetEntries();
-
- for (int e = 0; e < Nentries; e++) {
- // Printing status
- if (e % 100000 == 0){
- cout << "Please wait, we are " << e * 100. / Nentries << "% done ..\r" << flush;
- }
- chain->GetEntry(e);
-
-
- vector<double> ICcorr_Y(11), ICcorr_X(11); // the [0] element of spline is the
- // correction on the first
- // segment of ic
- for (int i = 0; i < ICcorr_Y.size() ; i++) {
- if (i == 0){
- ICcorr_Y.at(1) = IC->fIC_raw[1]*spline_Y[i]->Eval(0)/spline_Y[i]->Eval(FF_IC_Y);
- }
-
- else if (i == 1) {
- double temp = ICcorr_Y.at(1) / IC->fIC_raw[0] * spline_Y[i]->Eval(0)/ spline_Y[i]->Eval(FF_IC_Y);
- ICcorr_Y.at(0) = ICcorr_Y.at(1) / temp;
- if (!(ICcorr_Y.at(0)==ICcorr_Y.at(0))) ICcorr_Y.at(0) = 0;
- //cout << ICcorr_Y.at(0) << " " << IC->fIC_raw[0]<< endl;
- }
- else if (i > 1) {
- double temp = ICcorr_Y.at(1) / IC->fIC_raw[i] * spline_Y[i]->Eval(0)/ spline_Y[i]->Eval(FF_IC_Y);
- ICcorr_Y.at(i) = ICcorr_Y.at(1) / temp;
- if (!(ICcorr_Y.at(i)==ICcorr_Y.at(i))) ICcorr_Y.at(i) = 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(seg==1){
+ hICX.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_X",seg)));
+ hICY.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_Y",seg)));
+ 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)));
+ 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)));
+ hICcorr.push_back(new TH2F(Form("hICorr%d_Y",seg),Form("hICorr%d_Y",seg),1000,-100,100,500,0,2));
+ }
+ }
+
+ //===========================================================================================================
+ // 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 ++;
+ }
+ }
+
+
+ vector<TSpline3*> spline_X(11), spline_Y(11);
+ for (int i = 0; i < SplineCount ; i++) {
+ int seg = int((i-2)/2 +1 );
+
+ cout << " i " << i << " seg "<< seg << endl;
+ if (i < 2){
+ spline_X.at(0) = (TSpline3 *)fSpline->Get("fspline_IC1/IC0_X");
+ spline_Y.at(0) = (TSpline3 *)fSpline->Get("fspline_IC1/IC0_Y");
+ }
+
+ else if ( i>=2 && i<4){
+ spline_X.at(1) = (TSpline3 *)fSpline->Get("fsplineIC1_X");
+ spline_Y.at(1) = (TSpline3 *)fSpline->Get("fsplineIC1_Y");
+ }
+ else if (seg >=2 && (i%2 == 0) ) {
+ spline_X.at(seg) = (TSpline3 *)fSpline->Get(Form("fspline_IC%d/IC%d_X",seg,seg-1));
+ }
+ else if (seg >=2 && !(i%2 == 0) ) {
+ spline_Y.at(seg) = (TSpline3 *)fSpline->Get(Form("fspline_IC%d/IC%d_Y",seg,seg-1));
+ }
+ } //End loop on histogram
+ cout << spline_Y[0]->Eval(3)<< endl;
+ //===========================================================================================================
+ // 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
+
+
+ //===========================================================================================================
+ // Event Loop
+ //===========================================================================================================
+ int Nentries = chain->GetEntries();
+
+ for (int e = 0; e < Nentries; e++) {
+ // Printing status
+ if (e % 100000 == 0){
+ cout << "Please wait, we are " << e * 100. / Nentries << "% done ..\r" << flush;
+ }
+ chain->GetEntry(e);
+
+
+ vector<double> ICcorr_Y(11), ICcorr_X(11); // 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) {
+ cout << "bug " <<seg << endl;
+ 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();
} // 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/SplineChioP2P.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C
index 203af5bfc..e1be9a258 100644
--- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C
@@ -16,21 +16,24 @@ using namespace std;
TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double
XMax, double YMin, double YMax);
-vector<TH2F*>HistoFillerIC(int segment, vector<TSpline3*> SplineICprev);
+vector<TH2F*>HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_ICprev);
+vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Spline_All_ICprev);
void HistoDrawer(vector<TH2F*> h ,const char* CharName);
//////////////////////////////////////////////////////////////////////////////////////////
void SplineChioP2P() {
- int NSegment = 11 ;
+ int NSegment = 5 ;
double XMin[4] , XMax[4], YMin[4] , YMax[4];
XMin[0] = -520 ; XMax[0] = 380 ; YMin[0] = 2000 ; YMax[0] = 9500 ;
XMin[1] = -60 ; XMax[1] = 63 ; YMin[1] = 2000 ; YMax[1] = 9500 ;
XMin[2] = -520 ; XMax[2] = 380 ; YMin[2] = 0 ; YMax[2] = 2 ;
XMin[3] = -60 ; XMax[3] = 63 ; YMin[3] = 0 ; YMax[3] = 2 ;
-
+ //Reinitialising the outfile
+ TFile *OutFile = new TFile("Output/HistoP2P.root","recreate");
+ OutFile->Close();
//Initialising with IC1
int NCallSpline = 0;
int seg = 1;
@@ -38,19 +41,19 @@ void SplineChioP2P() {
vector<vector<TSpline3*>> SplineAllIC(11);
vector<TSpline3*> SplineICurrent;
- vector<vector<TH2F*>> hIC(11);
- hIC.at(1)= (HistoFillerIC(1,SplineICurrent));
+ 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(NCallSpline)->SetName("fsplineIC1_X");
+ 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(NCallSpline)->SetName("fsplineIC1_Y");
+ SplineICurrent.at(1)->SetName("fsplineIC1_Y");
NCallSpline+=1;
- SplineAllIC.at(1)= (SplineICurrent);
+ SplineAllIC.at(1)= (SplineICurrent); //
for (int seg = 0 ; seg < NSegment ; seg++){
@@ -60,8 +63,7 @@ void SplineChioP2P() {
else{
//Load Histo
- if (seg == 0) hIC.at(seg) = (HistoFillerIC(seg,SplineAllIC.at(1)));
- else if (seg != 0) hIC.at(seg) = (HistoFillerIC(seg,SplineAllIC.at(seg-1)));
+ hIC.at(seg) = (HistoFillerIC(seg,SplineAllIC));
//Reset spline holder
vector<TSpline3*> SplineICurrent;
@@ -74,20 +76,39 @@ void SplineChioP2P() {
// Y Spline
SplineICurrent.push_back(MakeSpline(hIC[seg][1] , NCallSpline, XMin[3] , XMax[3], YMin[3], YMax[3]));
- if (seg == 0) SplineICurrent.at(0)->SetName(Form("fspline_IC%d/IC%d_Y",seg+1,seg));
+ if (seg == 0) SplineICurrent.at(1)->SetName(Form("fspline_IC%d/IC%d_Y",seg+1,seg));
else if (seg != 0)SplineICurrent.at(1)->SetName(Form("fspline_IC%d/IC%d_Y",seg,seg-1));
NCallSpline+=1;
+ //push vector into memory
SplineAllIC.at(seg)= SplineICurrent;
}
}
-
+ //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();
}
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);
+
+ TCanvas *c1 = new TCanvas("c1","c1");
+ c1->Divide(NSegment);
+ for (int i=0 ; i<NSegment ; i++){
+ hIC_Corr_profile.at(i) = (TH1F*)hIC_Corr.at(i)->ProfileX();
+ c1->cd(i+1);
+ hIC_Corr_profile.at(i)->Draw();
+ }
+
+
} // End spline chio XY
///////////////////////////////////////////////////////////////////////////////////////////
@@ -95,7 +116,7 @@ void SplineChioP2P() {
* This function fill the histograms of ICseg/ICcorrSeg-1 expect for IC1 where
* it give IC1 and ICO where it give IC0/IC1corr
*/
-vector<TH2F*> HistoFillerIC(int segment, vector<TSpline3*> SplineICprev) {
+vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_ICprev) {
// Input and setters
TChain* chain = new TChain("PhysicsTree");
@@ -111,19 +132,26 @@ vector<TH2F*> HistoFillerIC(int segment, vector<TSpline3*> SplineICprev) {
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
- //spline to correct IC1 as a baseline of IC0
- TSpline3* SplineIC_X;
- TSpline3* SplineIC_Y;
+ // Get the number of previous spline to load
+ int NSpline = abs(segment-1);
+
+ //Filling all previous vector spline
+ vector<TSpline3*> SplineIC_X(NSpline+1);
+ vector<TSpline3*> SplineIC_Y(NSpline+1);
- if(!SplineICprev.empty()){
- SplineIC_X = SplineICprev.at(0);
- SplineIC_Y = SplineICprev.at(1);
+ //Fill all the previous spline , if the segment is 0 it'll fill spline 1 only
+ for (int pseg = 1 ; pseg<=NSpline; pseg++){
+ if(!Spline_All_ICprev.at(pseg).empty()){
+ SplineIC_X.at(pseg) = Spline_All_ICprev.at(pseg).at(0);
+ SplineIC_Y.at(pseg) = Spline_All_ICprev.at(pseg).at(1);
+ }
}
+
// Histograms
TH2F *hIC_X , *hIC_Y;
if (segment == 1 ) {
- hIC_Y = new TH2F(Form("hChio_IC%d__Y",segment),
+ hIC_Y = new TH2F(Form("hChio_IC%d_Y",segment),
Form("hChio_IC%d_Y",segment), 1000, -160, 160, 500, 2000, 9500);
hIC_X = new TH2F(Form("hChio_IC%d_X",segment),
Form("hChio_IC%d_X",segment), 1000, -800, 800, 500, 2000, 9500);
@@ -131,14 +159,14 @@ vector<TH2F*> HistoFillerIC(int segment, vector<TSpline3*> SplineICprev) {
}
else if (segment == 0){
- hIC_Y = new TH2F(Form("hChio_IC%d/IC%d__Y",segment+1,segment),
+ hIC_Y = new TH2F(Form("hChio_IC%d/IC%d_Y",segment+1,segment),
Form("hChio_IC%d/IC%d_Y",segment+1,segment), 1000, -160, 160, 500, 0, 2);
hIC_X = new TH2F(Form("hChio_IC%d/IC%d_X",segment+1,segment),
Form("hChio_IC%d/IC%d_X",segment+1,segment), 1000, -800, 800, 500, 0, 2);
}
else {
- hIC_Y = new TH2F(Form("hChio_IC%d/IC%d__Y",segment,segment-1),
+ hIC_Y = new TH2F(Form("hChio_IC%d/IC%d_Y",segment,segment-1),
Form("hChio_IC%d/IC%d_Y",segment,segment-1), 1000, -160, 160, 500, 0, 2);
hIC_X = new TH2F(Form("hChio_IC%d/IC%d_X",segment,segment-1),
Form("hChio_IC%d/IC%d_X",segment,segment-1), 1000, -800, 800, 500, 0, 2);
@@ -149,37 +177,66 @@ vector<TH2F*> HistoFillerIC(int segment, vector<TSpline3*> SplineICprev) {
// Beginning loop on entries
//===========================================================================================================
- int Nentries = chain->GetEntries();
-
+ // int Nentries = chain->GetEntries();
+ int Nentries = 100000;
+ 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 .. for segment " << segment << "\r" << flush;
- }
chain->GetEntry(e);
- //If the spline to correct IC1 doesnt exist we dont load it
- double IC_Spline_CorrX = IC->fIC_raw[segment];
- double IC_Spline_CorrY = IC->fIC_raw[segment];
-
- //Fetch corrected version of previous segment
- if ( !SplineICprev.empty() && segment !=0) {
- IC_Spline_CorrX = IC->fIC_raw[segment-2] * SplineIC_X->Eval(0) / SplineIC_X->Eval(FF_IC_Y);
- IC_Spline_CorrY = IC->fIC_raw[segment-2] * SplineIC_Y->Eval(0) / SplineIC_Y->Eval(FF_IC_Y);
- }
- else if (!SplineICprev.empty() && (segment == 0|| segment == 2)){
- IC_Spline_CorrX = IC->fIC_raw[1] * SplineIC_X->Eval(0) / SplineIC_X->Eval(FF_IC_Y);
- IC_Spline_CorrY = IC->fIC_raw[1] * SplineIC_Y->Eval(0) / SplineIC_Y->Eval(FF_IC_Y);
- }
+ vector<double> IC_Spline_CorrX(NSpline + 1),IC_Spline_CorrY(NSpline+1);
+
+ //Get the corrected value for all previous spline
+ for (int pseg = 1 ; pseg<=NSpline ; pseg++){
+ //If the spline doesnt exist define default value
+ 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];
+ }
+
+ //If it exist correct the current line
+ else if (!Spline_All_ICprev.at(pseg).empty()){
+ if(pseg ==1){
+ IC_Spline_CorrX.at(pseg) = IC->fIC_raw[pseg] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
+ IC_Spline_CorrY.at(pseg) = IC->fIC_raw[pseg] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
+ } //end ic 1
+ else {
+ IC_Spline_CorrX.at(pseg) = IC->fIC_raw[pseg] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
+ IC_Spline_CorrY.at(pseg) = IC->fIC_raw[pseg] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
+ } //end ic 2 to 10
+ } //End if non empty
+
+ } // End loop IC non 0
+
+ //===========================================================================================================
+ // Fill histo
+ //===========================================================================================================/
if (segment == 1){
hIC_Y->Fill(FF_IC_Y,IC->fIC_raw[segment]);
hIC_X->Fill(FF_IC_X,IC->fIC_raw[segment]);
}
+
+ else if (segment == 0){
+ hIC_Y->Fill(FF_IC_Y,IC_Spline_CorrY.at(1)/IC->fIC_raw[segment]);
+ hIC_X->Fill(FF_IC_X,IC_Spline_CorrX.at(1)/IC->fIC_raw[segment]);
+ }
+
+
else {
- hIC_Y->Fill(FF_IC_Y,IC->fIC_raw[segment]/IC_Spline_CorrY);
- hIC_X->Fill(FF_IC_X,IC->fIC_raw[segment]/IC_Spline_CorrX);
+ hIC_Y->Fill(FF_IC_Y,IC->fIC_raw[segment]/IC_Spline_CorrY.at(segment-1));
+ hIC_X->Fill(FF_IC_X,IC->fIC_raw[segment]/IC_Spline_CorrX.at(segment-1));
+ }
+
+ // Estimate time left every `updateInterval` iterations
+ 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 << "Treating segment : " << segment
+ << " | Estimated time left: " << timeLeft << " seconds" << "\r" << flush;
}
} // end loop on entries
@@ -189,9 +246,159 @@ vector<TH2F*> HistoFillerIC(int segment, vector<TSpline3*> SplineICprev) {
vector<TH2F*> hIC;
hIC.push_back(hIC_X);
hIC.push_back(hIC_Y);
+
+ TFile *OutFile = new TFile("Output/HistoP2P.root","UPDATE");
+ hIC_X->Write();
+ hIC_Y->Write();
+ OutFile->Close();
+
return hIC;
} // End HistoFiller
+vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Spline_All_ICprev){
+
+ // Input and setters
+ 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);
+
+ // Get the number of previous spline to load
+ segment = segment -1;
+ int NSpline = abs(segment);
+
+ //Filling all previous vector spline
+ vector<TSpline3*> SplineIC_X(NSpline+1);
+ vector<TSpline3*> SplineIC_Y(NSpline+1);
+
+ //Fill all the previous spline , if the segment is 0 it'll fill spline 1 only
+ for (int pseg = 0 ; pseg<=NSpline; pseg++){
+ if(!Spline_All_ICprev.at(pseg).empty()){
+ SplineIC_X.at(pseg) = Spline_All_ICprev.at(pseg).at(0);
+ SplineIC_Y.at(pseg) = Spline_All_ICprev.at(pseg).at(1);
+ }
+ }
+
+ // Histograms
+ vector<TH2F*> hIC_X(segment +1) , hIC_Y(segment+1);
+
+ for (int pseg = 0 ; pseg<=NSpline; pseg++){
+ if (pseg == 1 ) {
+ hIC_Y.at(pseg) = new TH2F(Form("hChiocorr_IC%d_Y",pseg),
+ Form("hChiocorr_IC%d_Y",pseg), 1000, -160, 160, 500, 2000, 9500);
+ hIC_X.at(pseg) = new TH2F(Form("hChiocorr_IC%d_X",pseg),
+ Form("hChiocorr_IC%d_X",pseg), 1000, -800, 800, 500, 2000, 9500);
+
+ }
+
+ else if (pseg == 0){
+ hIC_Y.at(pseg) = new TH2F(Form("hChiocorr_IC%d/IC%d_Y",pseg+1,pseg),
+ Form("hChiocorr_IC%d/IC%d_Y",pseg+1,pseg), 1000, -160, 160, 500, 0, 2);
+ hIC_X.at(pseg) = new TH2F(Form("hChiocorr_IC%d/IC%d_X",pseg+1,pseg),
+ Form("hChiocorr_IC%d/IC%d_X",pseg+1,pseg), 1000, -800, 800, 500, 0, 2);
+ }
+
+ else {
+ hIC_Y.at(pseg) = new TH2F(Form("hChiocorr_IC%d/IC%d_Y",pseg,pseg-1),
+ Form("hChiocorr_IC%d/IC%d_Y",pseg,pseg-1), 1000, -160, 160, 500, 0, 2);
+ hIC_X.at(pseg) = new TH2F(Form("hChiocorr_IC%d/IC%d_X",pseg,pseg-1),
+ Form("hChiocorr_IC%d/IC%d_X",pseg,pseg-1), 1000, -800, 800, 500, 0, 2);
+ }
+ }
+
+ //===========================================================================================================
+ // Beginning loop on entries
+ //===========================================================================================================
+
+ // int Nentries = chain->GetEntries();
+ int Nentries = 100000;
+ auto start = std::chrono::high_resolution_clock::now();
+ for (int e = 0; e < Nentries; e++) {
+
+ chain->GetEntry(e);
+
+ vector<double> IC_Spline_CorrX(NSpline + 1),IC_Spline_CorrY(NSpline+1);
+
+ //Get the corrected value for all previous spline
+ for (int pseg = 0 ; pseg<=NSpline ; pseg++){
+ //If the spline doesnt exist define default value
+ 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];
+ }
+
+ //If it exist correct the current line
+ else if (!Spline_All_ICprev.at(pseg).empty()){
+ if(pseg ==1){
+ IC_Spline_CorrX.at(pseg) = IC->fIC_raw[1] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
+ IC_Spline_CorrY.at(pseg) = IC->fIC_raw[1] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
+ } //end ic 1
+ else if (pseg !=0) {
+ IC_Spline_CorrX.at(pseg) = IC->fIC_raw[pseg] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
+ IC_Spline_CorrY.at(pseg) = IC->fIC_raw[pseg] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
+ } //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);
+ IC_Spline_CorrY.at(0) = IC->fIC_raw[0] / SplineIC_Y.at(0)->Eval(0) * SplineIC_Y.at(0)->Eval(FF_IC_Y);
+ }
+
+ for (int pseg = 0 ; pseg<=NSpline ; pseg++){
+ if (pseg == 1){
+ hIC_Y.at(pseg)->Fill(FF_IC_Y,IC_Spline_CorrY.at(1));
+ hIC_X.at(pseg)->Fill(FF_IC_X,IC_Spline_CorrY.at(1));
+ }
+
+ else if (pseg == 0){
+ hIC_Y.at(pseg)->Fill(FF_IC_Y,IC_Spline_CorrY.at(1)/IC_Spline_CorrY.at(pseg));
+ hIC_X.at(pseg)->Fill(FF_IC_X,IC_Spline_CorrY.at(1)/IC_Spline_CorrX.at(pseg));
+ }
+
+
+ else {
+ hIC_Y.at(pseg)->Fill(FF_IC_Y,IC_Spline_CorrY.at(pseg)/IC_Spline_CorrY.at(pseg-1));
+ hIC_X.at(pseg)->Fill(FF_IC_X,IC_Spline_CorrY.at(pseg)/IC_Spline_CorrX.at(pseg-1));
+ }
+ }
+ //===========================================================================================================
+ // Fill histo
+ //===========================================================================================================/
+
+
+ // Estimate time left every `updateInterval` iterations
+ 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 << "Treating segment : " << segment
+ << " | Estimated time left: " << timeLeft << " seconds" << "\r" << flush;
+ }
+ } // end loop on entries
+
+ //===========================================================================================================
+ // Output
+ //===========================================================================================================
+ vector<vector<TH2F*>> hIC;
+ hIC.push_back(hIC_X);
+ hIC.push_back(hIC_Y);
+
+ return hIC;
+} // End HistoFiller
+
+
/**
}
* Create the spline for one TH2F
@@ -244,7 +451,6 @@ TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double
LastBin = bin;
if (hProfile->GetBinLowEdge(bin)< XMax) break;
}
- cout << FirstBin << " " << LastBin << endl;
//===========================================================================================================
// Init Extended profile function
//===========================================================================================================
--
GitLab