Skip to content
Snippets Groups Projects
Commit cceb0f13 authored by Theodore Efremov's avatar Theodore Efremov :hibiscus:
Browse files

Added iterative spline maker for chio

parent 9e36743d
No related branches found
No related tags found
1 merge request!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...
Pipeline #363883 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
+ 84
− 166
View file @ cceb0f13
...@@ -30,48 +30,34 @@ void ICCorr() { ...@@ -30,48 +30,34 @@ void ICCorr() {
chain->SetBranchAddress("IC", &IC); chain->SetBranchAddress("IC", &IC);
TFile* fHisto = TFile::Open("Output/Histo.root"); TFile* fHisto = TFile::Open("Output/SplineInputAll.root");
TFile* fSpline = TFile::Open("Output/spline_Chio_2024.root"); TFile* fSpline = TFile::Open("Output/spline_ICALL_2024.root");
TFile* fFit = TFile::Open("Output/FitIC.root"); TFile* fFit = TFile::Open("Output/FitIC.root");
//Defining output histo //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_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* hDEBis_Ecorr_IC1 = new
TH2F("DEBis_Ecorr_IC1","DEBis_Ecorr_IC1",1000,0,20000,1000,0,25000); TH2F("DEBis_Ecorr_IC1","DEBis_Ecorr_IC1",1000,0,20000,1000,0,25000);
TH2F* hDE_E_IC1_Corr = new TH2F* hDE_E_CorrAll = new
TH2F("hDE_E_IC1_Corr","hDE_E_IC1_Corr",1000,0,20000,1000,0,25000); TH2F("hDE_E_CorrAll","hDE_E_CorrAll",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* 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 // Filling the input histo
hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_X_all")); int NSegment = 11 ;
hIC.push_back((TH2F*)fHisto->Get("hChio_IC0_IC1_Y_all")); 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 // GETSPLINE
...@@ -88,21 +74,26 @@ void ICCorr() { ...@@ -88,21 +74,26 @@ void ICCorr() {
} }
vector<TSpline3*> spline_X, spline_Y;
vector<TSpline3*> spline;
for (int i = 0; i < SplineCount ; i++) { for (int i = 0; i < SplineCount ; i++) {
cout << i << endl; if (i == 0){
if (i<2) spline.push_back((TSpline3*)fSpline->Get(Form("fsplineIC01_%d",i))); spline_X.push_back((TSpline3 *)fSpline->Get("fsplineIC1_X"));
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")));
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 } //End loop on histogram
//===========================================================================================================
// get fit //===========================================================================================================
//=========================================================================================================== // get fit
//===========================================================================================================
int FitCount = GetNumberKey(fFit,"TF1"); int FitCount = GetNumberKey(fFit,"TF1");
...@@ -126,158 +117,85 @@ void ICCorr() { ...@@ -126,158 +117,85 @@ void ICCorr() {
chain->GetEntry(e); chain->GetEntry(e);
// Correction 1on0 vector<double> ICcorr_Y(11), ICcorr_X(11); // the [0] element of spline is the
double IC_CorrX = IC->fIC_raw[1]/ IC->fIC_raw[0] * spline[0]->Eval(0) / // correction on the first
spline[0]->Eval(FF_IC_X) ; // segment of ic
for (int i = 0; i < ICcorr_Y.size() ; i++) {
double IC_temp_CorrX = IC->fIC_raw[1] / IC_CorrX ; if (i == 0){
ICcorr_Y.at(1) = IC->fIC_raw[1]*spline_Y[i]->Eval(0)/spline_Y[i]->Eval(FF_IC_Y);
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);
//Correction 1 else if (i == 1) {
double IC0_Corr = IC->fIC_raw[1] /IC_CorrY; double temp = ICcorr_Y.at(1) / IC->fIC_raw[0] * spline_Y[i]->Eval(0)/ spline_Y[i]->Eval(FF_IC_Y);
if (IC0_Corr != IC0_Corr ) IC0_Corr = 0; 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 for (int i = 0; i < ICcorr_Y.size() ; i++) {
double IC0_CorrY = IC->fIC_raw[0] * spline[4]->Eval(0)/ spline[4]->Eval(FF_IC_Y); hICcorr.at(i)->Fill(FF_IC_Y,ICcorr_Y.at(i));
double IC0_CorrX = IC->fIC_raw[0] * spline[3]->Eval(0)/ spline[3]->Eval(FF_IC_X); } //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++ ){ 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] ; if (seg >4){
else if (seg >4) E+= IC->fIC_raw[seg] ; E += IC->fIC_raw[seg] ;
} Ecorr += ICcorr_Y.at(seg) ;
}
//===========================================================================================================
// Filling DE with ICO
//===========================================================================================================
DEcorr += DE ;
DE_IC0_Corr_Y += DE;
if ((IC0_Corr < 1e100 && IC0_Corr >-1000) || true){
DEcorr += IC0_Corr ;
} }
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 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*(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]; DE = 0.5*(IC->fIC_raw[1] +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]; DE_IC1_corrY = 0.5*(ICcorr_Y[1] +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_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){ if (FF_IC_Y >-60 && FF_IC_Y <60){
hDEBis_Ecorr_IC1->Fill(E,DE_Bis); 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_E->Fill(E,DE);
hDE_Ecorr01->Fill(E,DEcorr); hDE_Ecorr1->Fill(E,DE_IC1_corrY);
hDE_E_IC1_Corr->Fill(E,DEIC1corr); hDE_E_CorrAll->Fill(E,DE_ICALL_corrY);
} }
double DEcut; } //endl loop event
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
//=========================================================================================================== //===========================================================================================================
// Drawing histo // Drawing histo
//=========================================================================================================== //===========================================================================================================
TCanvas* c1 = new TCanvas("c1","c1",1000,1000); TCanvas* c1 = new TCanvas("c1","c1",10000,1000);
c1->Divide(2,2); c1->Divide(11,2);
c1->cd(1); for (int i = 0 ; i<11 ;i++){
hICorr[0]->Draw(); c1->cd(i+1);
c1->cd(2); hICcorr.at(i)->ProfileX()->Draw();
hIC[0]->Draw();
c1->cd(3); c1->cd(11+i+1);
hICorr[1]->Draw(); hICY.at(i)->ProfileX()->Draw();
c1->cd(4); }
hIC[1]->Draw();
TCanvas* c2 = new TCanvas("c2","c2",1500,1000); TCanvas* c2 = new TCanvas("c2","c2",1500,1000);
c2->Divide(4); c2->Divide(4);
c2->cd(1); c2->cd(1);
gPad->SetLogz(); gPad->SetLogz();
hDE_Ecorr_IC0_IC1corr->Draw(); hDE_E->Draw();
c2->cd(2); c2->cd(2);
gPad->SetLogz(); gPad->SetLogz();
hDE_Ecorr_IC1->Draw(); hDE_Ecorr1->Draw();
c2->cd(3); c2->cd(3);
gPad->SetLogz(); gPad->SetLogz();
hDEBis_Ecorr_IC1->Draw(); hDEBis_Ecorr_IC1->Draw();
c2->cd(4); c2->cd(4);
gPad->SetLogz(); gPad->SetLogz();
hDE_E_ICFit_Corr->Draw(); hDE_E_CorrAll->Draw();
} // End spline chio XY } // End spline chio XY
......
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioAllXY.C
+ 36
− 30
View file @ cceb0f13
...@@ -21,7 +21,6 @@ void HistoDrawer(vector<TH2F*> h ,const char* CharName); ...@@ -21,7 +21,6 @@ void HistoDrawer(vector<TH2F*> h ,const char* CharName);
////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////
void SplineChioAllXY( bool Create = false , bool Cut = false) { void SplineChioAllXY( bool Create = false , bool Cut = false) {
TFile* inFile = TFile::Open("Output/SplineInputAll.root");
if (Create == true){ if (Create == true){
cout << "Creating SplineInputAll.root" << endl; cout << "Creating SplineInputAll.root" << endl;
...@@ -29,8 +28,9 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) { ...@@ -29,8 +28,9 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) {
cout << "SplineInputAll.root created and loaded" << endl ; 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; cout << "Histogram loaded" << endl;
} }
...@@ -38,6 +38,7 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) { ...@@ -38,6 +38,7 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) {
cerr << "Make your mind" << endl; cerr << "Make your mind" << endl;
} }
int NSegment = 11 ; int NSegment = 11 ;
vector<TH2F*> hIC,hICX,hICY; vector<TH2F*> hIC,hICX,hICY;
hIC.push_back((TH2F*)inFile->Get("hChio_IC1_X")); hIC.push_back((TH2F*)inFile->Get("hChio_IC1_X"));
...@@ -68,19 +69,23 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) { ...@@ -68,19 +69,23 @@ void SplineChioAllXY( bool Create = false , bool Cut = false) {
for (int i = 0; i < NHist; i++) { for (int i = 0; i < NHist; i++) {
TSpline3 *spline; 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 = 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) { 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]); 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)) { 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]); YMax[3]);
spline->SetName(Form("fsplineIC1_IC%d_Y",i)); spline->SetName(Form("fsplineIC1_IC%d_Y",int((i-2)/2)));
} }
spline->Write(); spline->Write();
} //End loop on histogram } //End loop on histogram
...@@ -112,7 +117,7 @@ void HistoFiller(bool cut) { ...@@ -112,7 +117,7 @@ void HistoFiller(bool cut) {
//spline to correct IC1 as a baseline of IC0 //spline to correct IC1 as a baseline of IC0
TFile* fSpline = TFile::Open("Output/spline_ICALL_2024.root"); 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 //cut
...@@ -121,20 +126,20 @@ void HistoFiller(bool cut) { ...@@ -121,20 +126,20 @@ void HistoFiller(bool cut) {
CutZ = (TCutG*)fcut->Get("CutTry") ; CutZ = (TCutG*)fcut->Get("CutTry") ;
// Histograms // 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; 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++){ for (int seg = 0 ; seg< NSegment; seg++){
hIC1_ICn_X.at(seg) = new TH2F(Form("hIC1_IC%d_X",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); 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), 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); 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) { ...@@ -186,28 +191,30 @@ void HistoFiller(bool cut) {
if (fSpline && !fSpline->IsZombie()) { if (fSpline && !fSpline->IsZombie()) {
IC1_Spline_Corr = IC->fIC_raw[1] * SplineIC1->Eval(0) / SplineIC1->Eval(FF_IC_Y); 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 }//loop on segment
} // end else } // 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"); TFile* outFile = new TFile("Output/SplineInputAll.root", "recreate");
hChio_IC1_Y->Write();
hChio_IC1_X->Write();
for(int seg=0; seg<NSegment; seg++){ for(int seg=0; seg<NSegment; seg++){
hIC1_ICn_X.at(seg)->Write(); hIC1_ICn_X.at(seg)->Write();
hIC1_ICn_Y.at(seg)->Write(); hIC1_ICn_Y.at(seg)->Write();
hChio_ICn_Y.at(seg)->Write();
hChio_ICn_X.at(seg)->Write();
}//loop on segment }//loop on segment
outFile->Close(); outFile->Close();
} // End HistoFiller } // End HistoFiller
...@@ -270,11 +277,10 @@ TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double ...@@ -270,11 +277,10 @@ TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double
//=========================================================================================================== //===========================================================================================================
// Create the extended TProfile // Create the extended TProfile
pfx = new TH1F( 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())); hProfile->GetBinLowEdge(hProfile->GetNbinsX()) + hProfile->GetBinWidth(hProfile->GetNbinsX()));
pfx->SetLineColor(8); pfx->SetLineColor(8);
pfx->SetDirectory(0); pfx->SetDirectory(0);
//=========================================================================================================== //===========================================================================================================
// Fill extended profile // Fill extended profile
//=========================================================================================================== //===========================================================================================================
......
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C 0 → 100644
+ 305
− 0
View file @ cceb0f13
//#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");
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment