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

P2P chio correction

parent 62df19d6
No related branches found
No related tags found
1 merge request!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...
Pipeline #364228 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
+ 174
− 164
View file @ b7c9ccec
......@@ -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;
}
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C
+ 250
− 44
View file @ b7c9ccec
......@@ -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
//===========================================================================================================
......
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