diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C index a12985b1e5dff7387e40138c37acb5115eeefbd1..660cd41778b524c99476e2f8aafc8c96e8a4801e 100644 --- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C @@ -30,48 +30,34 @@ void ICCorr() { chain->SetBranchAddress("IC", &IC); - TFile* fHisto = TFile::Open("Output/Histo.root"); - TFile* fSpline = TFile::Open("Output/spline_Chio_2024.root"); + TFile* fHisto = TFile::Open("Output/SplineInputAll.root"); + TFile* fSpline = TFile::Open("Output/spline_ICALL_2024.root"); TFile* fFit = TFile::Open("Output/FitIC.root"); //Defining output histo - vector<TH2F*> hIC,hICorr; - TH2F* hICorr_X = new TH2F("hICorr_X","hICorr_X",1000,-600,400,500,1,1.6); - TH2F* hICorr_Y = new TH2F("hICorr_Y","hICorr_Y",1000,-100,100,500,1,1.6); - TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,100,1,1000,100,1); - TH2F* hDE_Ecorr01 = new TH2F("DE_Ecorr01","DE_Ecorr01",1000,100,1,1000,100,1); + TH2F* hDE_Ecorr1 = new TH2F("DE_Ecorr1","DE_Ecorr1",1000,100,1,1000,100,1); TH2F* hDEBis_Ecorr_IC1 = new TH2F("DEBis_Ecorr_IC1","DEBis_Ecorr_IC1",1000,0,20000,1000,0,25000); - TH2F* hDE_E_IC1_Corr = new - TH2F("hDE_E_IC1_Corr","hDE_E_IC1_Corr",1000,0,20000,1000,0,25000); - - TH2F* hDE_E_IC0_Corr = new - TH2F("hDE_E_IC0_Corr","hDE_E_IC0_Corr",1000,0,20000,1000,0,25000); - TH2F* hDEcut= - new TH2F("DEcut","DEcut",1000,0,20000,1000,0,25000); - - TH2F* hDE_Ecorr_IC1 = new - TH2F("hDE_Ecorr_IC1","hDE_Ecorr_IC1",1000,0,20000,1000,0,25000); - - TH2F* hDE_Ecorr_IC0 = new - TH2F("hDE_Ecorr_IC0","hDE_Ecorr_IC0",1000,0,20000,1000,0,25000); - - TH2F* hDE_Ecorr_IC0_IC1corr = new - TH2F("hDE_Ecorr_IC0_IC1corr","hDE_Ecorr_IC0_IC1corr",1000,0,20000,1000,0,25000); - //Fit - TH2F* hDE_E_ICFit_Corr = new - TH2F("hDE_E_ICFit_Corr","hDE_E_ICFit_Corr",1000,0,20000,1000,0,25000); + TH2F* hDE_E_CorrAll = new + TH2F("hDE_E_CorrAll","hDE_E_CorrAll",1000,0,20000,1000,0,25000); - TH2F* hICorr_Y_Fit = new - TH2F("hICorr_Y_Fit","hICorr_Y_Fit",1000,-100,100,500,0,25000); - - hICorr.push_back(hICorr_X); - hICorr.push_back(hICorr_Y); // Filling the input histo - hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_X_all")); - hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_Y_all")); + int NSegment = 11 ; + vector<TH2F*> hIC1_ICn, hICX,hICY ,hIC1_ICnX,hIC1_ICnY, 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 @@ -88,21 +74,26 @@ void ICCorr() { } - - - vector<TSpline3*> spline; + vector<TSpline3*> spline_X, spline_Y; for (int i = 0; i < SplineCount ; i++) { - cout << i << endl; - if (i<2) spline.push_back((TSpline3*)fSpline->Get(Form("fsplineIC01_%d",i))); - else if ( i>=2 && i<4) spline.push_back((TSpline3*)fSpline->Get(Form("fsplineIC1_%d",i))); - - else if (i>=4 && i<6) spline.push_back((TSpline3 *)fSpline->Get(Form("fsplineIC0_%d",i))); - else if (i==6) spline.push_back((TSpline3 *)fSpline->Get(Form("fsplineIC0_IC1corr"))); + 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 - //=========================================================================================================== + + //=========================================================================================================== + // get fit + //=========================================================================================================== int FitCount = GetNumberKey(fFit,"TF1"); @@ -126,158 +117,85 @@ void ICCorr() { chain->GetEntry(e); - // Correction 1on0 - double IC_CorrX = IC->fIC_raw[1]/ IC->fIC_raw[0] * spline[0]->Eval(0) / - spline[0]->Eval(FF_IC_X) ; - - double IC_temp_CorrX = IC->fIC_raw[1] / IC_CorrX ; - - double IC_CorrY = IC->fIC_raw[1]/ IC_temp_CorrX * spline[1]->Eval(0) / - spline[1]->Eval(FF_IC_Y) ; - - hICorr[0]->Fill(FF_IC_X,IC_CorrX); - hICorr[1]->Fill(FF_IC_Y,IC_CorrY); - + 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); + } - //Correction 1 - double IC0_Corr = IC->fIC_raw[1] /IC_CorrY; - if (IC0_Corr != IC0_Corr ) IC0_Corr = 0; + 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 - double IC1_CorrY = IC->fIC_raw[1] + (spline[3]->Eval(0) - - spline[3]->Eval(FF_IC_Y)) ; - double IC1_CorrX = IC->fIC_raw[1] * spline[2]->Eval(0) / - spline[2]->Eval(FF_IC_X) ; - //Correction 0 - double IC0_CorrY = IC->fIC_raw[0] * spline[4]->Eval(0)/ spline[4]->Eval(FF_IC_Y); - double IC0_CorrX = IC->fIC_raw[0] * spline[3]->Eval(0)/ spline[3]->Eval(FF_IC_X); + 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; - //Double correction - double IC0_CorrIC1_Y = IC->fIC_raw[1]/ IC_temp_CorrX * spline[6]->Eval(0) / - spline[6]->Eval(FF_IC_Y) ; - double IC0_Corr_IC1 = IC->fIC_raw[1] /IC0_CorrIC1_Y; - - double DE = 0 ,E =0 ,DEcorr=0 ,DEIC1corr =0, DE_IC0_Corr_Y=0; - - //=========================================================================================================== - // Fill De from seg 1 to 3 - //=========================================================================================================== for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ - if(seg > 1 && seg < 4) DE += IC->fIC_raw[seg] ; - else if (seg >4) E+= IC->fIC_raw[seg] ; - } - - - - //=========================================================================================================== - // Filling DE with ICO - //=========================================================================================================== - DEcorr += DE ; - DE_IC0_Corr_Y += DE; - - if ((IC0_Corr < 1e100 && IC0_Corr >-1000) || true){ - DEcorr += IC0_Corr ; + if (seg >4){ + E += IC->fIC_raw[seg] ; + Ecorr += ICcorr_Y.at(seg) ; + } } - else DEcorr += IC->fIC_raw[0]; - if ((IC0_CorrX < 1e100 && IC0_CorrX >-1000) || true){ - DE_IC0_Corr_Y += IC0_CorrX ; - } - - else DE_IC0_Corr_Y += IC->fIC_raw[0]; - - DE += IC->fIC_raw[0] ; - - DE = DE *0.5 + IC->fIC_raw[4]; - DEcorr = DEcorr *0.5 + IC->fIC_raw[4]; - DE_IC0_Corr_Y = DE_IC0_Corr_Y *0.5 + IC->fIC_raw[4]; - - /* - cout << IC0_Corr << " " << IC->fIC_raw[0] << endl; - cout << DEcorr << " " << DE << endl; - cout << "*************************" << endl; - */ - - //=========================================================================================================== - // Chio online - //=========================================================================================================== - 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]; - //double DE_Bis = 0.5*(IC->fIC_raw[1]+IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; - DEIC1corr = 0.5*(IC1_CorrY +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; - double DEIC0corr = 0.5*(IC0_CorrY + IC1_CorrY +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; - double DEIC0_IC1corr = 0.5*(IC0_CorrIC1_Y + IC1_CorrY +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; + + 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_Ecorr_IC1->Fill(E,DEIC1corr); - hDE_Ecorr_IC0->Fill(E,DEIC0corr); - hDE_Ecorr_IC0_IC1corr->Fill(E,DEIC0_IC1corr); - hDE_E->Fill(E,DE); - hDE_Ecorr01->Fill(E,DEcorr); - hDE_E_IC1_Corr->Fill(E,DEIC1corr); + hDE_Ecorr1->Fill(E,DE_IC1_corrY); + hDE_E_CorrAll->Fill(E,DE_ICALL_corrY); } - double DEcut; - if ((FF_IC_Y>-60 && FF_IC_Y<-55)|| (FF_IC_Y>55 && FF_IC_Y<60)){ DEcut=DE_Bis; - hDEcut->Fill(E,DEcut); - } - //=========================================================================================================== - // - //=========================================================================================================== - //=========================================================================================================== - // Fit Corr - //=========================================================================================================== - - vector<double> IC_Corr_Fit(11); - double DEFit=0, EtotFit=0, EFit=0; - for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ - IC_Corr_Fit[seg] = IC->fIC_raw[seg]* Fit[seg]->Eval(0) /Fit[seg]->Eval(FF_IC_Y); - } - - DEFit = (IC_Corr_Fit[1] + IC->fIC_raw[2]+ IC->fIC_raw[3])*0.5 - + IC->fIC_raw[4]; - - EFit = IC->fIC_raw[5]+ IC->fIC_raw[6]+ IC->fIC_raw[7]+ IC->fIC_raw[8]+ - IC->fIC_raw[9]+ IC->fIC_raw[10]; - - EtotFit = DEFit + EFit ; - - hICorr_Y_Fit->Fill(FF_IC_Y,IC_Corr_Fit[0]); - hDE_E_ICFit_Corr ->Fill(EFit,DEFit); - } //End loop on event + } //endl loop event //=========================================================================================================== // Drawing histo //=========================================================================================================== - TCanvas* c1 = new TCanvas("c1","c1",1000,1000); - c1->Divide(2,2); - c1->cd(1); - hICorr[0]->Draw(); - c1->cd(2); - hIC[0]->Draw(); - c1->cd(3); - hICorr[1]->Draw(); - c1->cd(4); - hIC[1]->Draw(); + 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_Ecorr_IC0_IC1corr->Draw(); + hDE_E->Draw(); c2->cd(2); gPad->SetLogz(); - hDE_Ecorr_IC1->Draw(); + hDE_Ecorr1->Draw(); c2->cd(3); gPad->SetLogz(); hDEBis_Ecorr_IC1->Draw(); c2->cd(4); gPad->SetLogz(); - hDE_E_ICFit_Corr->Draw(); - + hDE_E_CorrAll->Draw(); } // End spline chio XY diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioAllXY.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioAllXY.C index 6a00498d37e24fb98f9f8379dda8ebf759d8437e..fe9f9cbbf07197a01ddb2cc2d0a612e49d7e6e16 100644 --- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioAllXY.C +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioAllXY.C @@ -21,7 +21,6 @@ void HistoDrawer(vector<TH2F*> h ,const char* CharName); ////////////////////////////////////////////////////////////////////////////////////////// void SplineChioAllXY( bool Create = false , bool Cut = false) { - TFile* inFile = TFile::Open("Output/SplineInputAll.root"); if (Create == true){ cout << "Creating SplineInputAll.root" << endl; @@ -29,8 +28,9 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) { cout << "SplineInputAll.root created and loaded" << endl ; } + TFile* inFile = TFile::Open("Output/SplineInputAll.root"); - else if (inFile && !inFile->IsZombie()) { + if (inFile && !inFile->IsZombie()) { cout << "Histogram loaded" << endl; } @@ -38,6 +38,7 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) { cerr << "Make your mind" << endl; } + int NSegment = 11 ; vector<TH2F*> hIC,hICX,hICY; hIC.push_back((TH2F*)inFile->Get("hChio_IC1_X")); @@ -68,19 +69,23 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) { for (int i = 0; i < NHist; i++) { TSpline3 *spline; - if (i < 2){ + if (i == 1){ + spline = MakeSpline(hIC[i] , i, XMin[i] , XMax[i], YMin[i], YMax[i]); + spline->SetName("fsplineIC1_Y"); + } + else if (i == 0){ spline = MakeSpline(hIC[i] , i, XMin[i] , XMax[i], YMin[i], YMax[i]); - spline->SetName("fsplineIC1"); + spline->SetName("fsplineIC1_X"); } else if (i>=2 && (i%2 == 0) && (i!=4) && false) { - TSpline3* spline = MakeSpline(hIC[i] , i, XMin[2] , XMax[2], YMin[2], + spline = MakeSpline(hIC[i] , i, XMin[2] , XMax[2], YMin[2], YMax[2]); - spline->SetName(Form("fsplineIC1_IC%d_X",i)); + spline->SetName(Form("fsplineIC1_IC%d_X",int((i-2)/2))); } else if (i>=2 && !(i%2 == 0) && (i!=5)) { - TSpline3* spline = MakeSpline(hIC[i] , i, XMin[3] , XMax[3], YMin[3], + spline = MakeSpline(hIC[i] , i, XMin[3] , XMax[3], YMin[3], YMax[3]); - spline->SetName(Form("fsplineIC1_IC%d_Y",i)); + spline->SetName(Form("fsplineIC1_IC%d_Y",int((i-2)/2))); } spline->Write(); } //End loop on histogram @@ -112,7 +117,7 @@ void HistoFiller(bool cut) { //spline to correct IC1 as a baseline of IC0 TFile* fSpline = TFile::Open("Output/spline_ICALL_2024.root"); - TSpline3* SplineIC1 = (TSpline3 *)fSpline->Get(Form("fsplineIC1")); + TSpline3* SplineIC1 = (TSpline3 *)fSpline->Get(Form("fsplineIC1_Y")); //cut @@ -121,20 +126,20 @@ void HistoFiller(bool cut) { CutZ = (TCutG*)fcut->Get("CutTry") ; // Histograms - - TH2F* hChio_IC1_Y = new TH2F("hChio_IC1_Y", - "hChio_IC1_Y", 1000, -160, 160, 500, 2000, 9500); - TH2F* hChio_IC1_X = new TH2F("hChio_IC1_X", - "hChio_IC1_X", 1000, -800, 800, 500, 2000, 9500); - double NSegment = 11; - vector<TH2F*> hIC1_ICn_X(NSegment), hIC1_ICn_Y(NSegment); + vector<TH2F*> hIC1_ICn_X(NSegment), hIC1_ICn_Y(NSegment), hChio_ICn_Y(NSegment) , hChio_ICn_X(NSegment) ; for (int seg = 0 ; seg< NSegment; seg++){ hIC1_ICn_X.at(seg) = new TH2F(Form("hIC1_IC%d_X",seg), Form("hIC1_IC%d_X",seg), 1000, -800, 800, 500, 0, 2); hIC1_ICn_Y.at(seg) = new TH2F(Form("hIC1_IC%d_Y",seg), Form("hIC1_IC%d_Y",seg), 1000, -160, 160, 500, 0, 2); + + hChio_ICn_Y.at(seg) = new TH2F(Form("hChio_IC%d_Y",seg), + Form("hChio_IC%d_Y",seg), 1000, -160, 160, 500, 2000, 9500); + hChio_ICn_X.at(seg) = new TH2F(Form("hChio_IC%d_X",seg), + Form("hChio_IC%d_X",seg), 1000, -800, 800, 500, 2000, 9500); + } @@ -186,28 +191,30 @@ void HistoFiller(bool cut) { if (fSpline && !fSpline->IsZombie()) { IC1_Spline_Corr = IC->fIC_raw[1] * SplineIC1->Eval(0) / SplineIC1->Eval(FF_IC_Y); } + for(int seg=0; seg<NSegment; seg++){ + hIC1_ICn_X.at(seg)->Fill(FF_IC_X,IC1_Spline_Corr/fIC_raw[seg]); + hIC1_ICn_Y.at(seg)->Fill(FF_IC_Y,IC1_Spline_Corr/fIC_raw[seg]); + + hChio_ICn_X.at(seg)->Fill(FF_IC_X, fIC_raw[seg]); + hChio_ICn_Y.at(seg)->Fill(FF_IC_Y, fIC_raw[seg]); - hChio_IC1_X->Fill(FF_IC_X, fIC_raw[1]); - hChio_IC1_Y->Fill(FF_IC_Y, fIC_raw[1]); - for(int seg=0; seg<NSegment; seg++){ - hIC1_ICn_X.at(seg)->Fill(FF_IC_X,IC1_Spline_Corr/fIC_raw[seg]); - hIC1_ICn_Y.at(seg)->Fill(FF_IC_Y,IC1_Spline_Corr/fIC_raw[seg]); }//loop on segment } // end else - } // end loop on entries + } // end loop on entries - //=========================================================================================================== - // Output Canvas and files - //=========================================================================================================== + //=========================================================================================================== + // Output Canvas and files + //=========================================================================================================== TFile* outFile = new TFile("Output/SplineInputAll.root", "recreate"); - hChio_IC1_Y->Write(); - hChio_IC1_X->Write(); - for(int seg=0; seg<NSegment; seg++){ hIC1_ICn_X.at(seg)->Write(); hIC1_ICn_Y.at(seg)->Write(); + hChio_ICn_Y.at(seg)->Write(); + hChio_ICn_X.at(seg)->Write(); + + }//loop on segment outFile->Close(); } // End HistoFiller @@ -270,11 +277,10 @@ TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double //=========================================================================================================== // Create the extended TProfile pfx = new TH1F( - Form("pfx_%02d", NCallee + 1), Form("pfx_%02d", NCallee + 1), hProfile->GetNbinsX(), hProfile->GetBinLowEdge(1), + 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 //=========================================================================================================== diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C new file mode 100644 index 0000000000000000000000000000000000000000..203af5bfc95fae1c02b346bf1e9d7b9d60f2f322 --- /dev/null +++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C @@ -0,0 +1,305 @@ +//#ifdef TICPhysics +#include <TICPhysics.h> +//#endif +#include <TCanvas.h> +#include <TChain.h> +#include <TF1.h> +#include <TFile.h> +#include <TH2.h> +#include <TSpline.h> +#include <fstream> +#include <vector> +#include <TCutG.h> + +using namespace std; + +TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double + XMax, double YMin, double YMax); + +vector<TH2F*>HistoFillerIC(int segment, vector<TSpline3*> SplineICprev); + +void HistoDrawer(vector<TH2F*> h ,const char* CharName); + +////////////////////////////////////////////////////////////////////////////////////////// +void SplineChioP2P() { + + int NSegment = 11 ; + 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 ; + + + //Initialising with IC1 + int NCallSpline = 0; + int seg = 1; + + vector<vector<TSpline3*>> SplineAllIC(11); + vector<TSpline3*> SplineICurrent; + + vector<vector<TH2F*>> hIC(11); + hIC.at(1)= (HistoFillerIC(1,SplineICurrent)); + + SplineICurrent.push_back(MakeSpline(hIC[seg][0] , NCallSpline, XMin[0] , XMax[0], YMin[0], YMax[0])); + SplineICurrent.at(NCallSpline)->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"); + NCallSpline+=1; + + SplineAllIC.at(1)= (SplineICurrent); + + + for (int seg = 0 ; seg < NSegment ; seg++){ + if(seg ==1) { + } + + 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))); + + //Reset spline holder + vector<TSpline3*> SplineICurrent; + + // X Spline + SplineICurrent.push_back(MakeSpline(hIC[seg][0] , NCallSpline, XMin[2] , XMax[2], YMin[2], YMax[2])); + if (seg == 0) SplineICurrent.at(0)->SetName(Form("fspline_IC%d/IC%d_X",seg+1,seg)); + else if (seg != 0) SplineICurrent.at(0)->SetName(Form("fspline_IC%d/IC%d_X",seg,seg-1)); + NCallSpline+=1; + + // 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)); + else if (seg != 0)SplineICurrent.at(1)->SetName(Form("fspline_IC%d/IC%d_Y",seg,seg-1)); + NCallSpline+=1; + + SplineAllIC.at(seg)= SplineICurrent; + } + } + + 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(); +} // End spline chio XY + +/////////////////////////////////////////////////////////////////////////////////////////// +/** + * 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) { + + // 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); + + //spline to correct IC1 as a baseline of IC0 + TSpline3* SplineIC_X; + TSpline3* SplineIC_Y; + + if(!SplineICprev.empty()){ + SplineIC_X = SplineICprev.at(0); + SplineIC_Y = SplineICprev.at(1); + } + // Histograms + TH2F *hIC_X , *hIC_Y; + + if (segment == 1 ) { + 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); + + } + + else if (segment == 0){ + 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), + 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); + } + + + //=========================================================================================================== + // Beginning loop on entries + //=========================================================================================================== + + 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 .. 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); + } + + if (segment == 1){ + hIC_Y->Fill(FF_IC_Y,IC->fIC_raw[segment]); + hIC_X->Fill(FF_IC_X,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); + } + } // end loop on entries + + //=========================================================================================================== + // Output + //=========================================================================================================== + vector<TH2F*> hIC; + hIC.push_back(hIC_X); + hIC.push_back(hIC_Y); + return hIC; +} // End HistoFiller + +/** + } + * Create the spline for one TH2F + * + * Input : + * Pointer to a TH2F + * + * + * Output: + * + * TSpline3 : the output spline + * + */ + +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; + } + cout << FirstBin << " " << LastBin << endl; + //=========================================================================================================== + // 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 + +void HistoDrawer(vector<TH2F*> h,const char* CharName){ + + int NHisto = static_cast<int>(h.size()); + + TCanvas* Can= new TCanvas(CharName,CharName,0,0,2000,1000); + + Can->Divide(NHisto); + + for ( int i=0 ; i < NHisto ; i+=1 ){ + //All x are even and all y are odd + int CanvaNumber = (i) ; + Can->cd(CanvaNumber); + h[i]->Draw("colz"); + } +}