diff --git a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C index 941c5e4be564c5f82947eae895a35a7f33b932d5..0fc833b1dd4f474186d1841d526474c0048d145f 100644 --- a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C +++ b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C @@ -8,19 +8,18 @@ #include <TH2.h> #include <TSpline.h> #include <fstream> -#include "../XYCalibration/ProfileEvaluator.h" using namespace std; -const int Ncuts=30; + +const int CutNumberLoader(); +const int Ncuts = CutNumberLoader(); const int NSegment=11; const int Nrun=1; -int RunNumber[Nrun]={241}; +int RunNumber[Nrun]={247}; void MakeSpline(); void ApplySpline_perTCutG(); void ApplySpline(); -vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline); -vector<double> TxtToVector(const char *Path); /////////////////////////////////////////////////////////////////////////// void SplineChio() @@ -31,76 +30,42 @@ void SplineChio() //**********************Cut************************************** TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root","open"); - TCutG *cutZ[Ncuts]; - for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i)); - - // ****************** Spline Y ********************************** - - TFile *fSplineDe; - TSpline3 *splineDE,*splineDE0, *splineDE02,*splineDE023,*splineDE0234,*splineDEbis; - fSplineDe= new TFile("../XYCalibration/Output/splineDE.root","open"); - splineDE = (TSpline3*) fSplineDe->Get("SplineDe"); - splineDE0 = (TSpline3*) fSplineDe->Get("SplineDe_0"); - splineDE02 = (TSpline3*) fSplineDe->Get("SplineDe_02"); - splineDE023 = (TSpline3*) fSplineDe->Get("SplineDe_023"); - splineDE0234 = (TSpline3*) fSplineDe->Get("SplineDe_0234"); - - - TFile *fSplineIC = new TFile("../XYCalibration/Output/spline_P2P_2024.root","open"); - - vector<vector<TSpline3*>> Spline = LoadSpline("../XYCalibration/Output/spline_P2P_2024.root"); - vector<TSpline3*> spline_X(11), spline_Y(11); - spline_X = Spline.at(0) ; spline_Y = Spline.at(1); - + vector<TCutG*> cutZ(Ncuts); + for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i)); //**********************DE_E**************************************************** TChain *chain=new TChain("PhysicsTree"); - chain->Add("../../../root/analysis/VamosCalib247.root"); + for(int i=0;i<Nrun;i++) + { + chain->Add(Form("../../../root/analysis/Run%d.root",RunNumber[i])); + cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl; + } TICPhysics* IC = new TICPhysics() ; TTimeData *Time = new TTimeData() ; - TFPMWPhysics *FPMW = new TFPMWPhysics(); - double FF_IC_X, FF_IC_Y, FF_V13; - double FF_DE, FF_Eres; chain->SetBranchStatus("IC", true); chain->SetBranchAddress("IC", &IC); - chain->SetBranchStatus("FF_IC_Y", true); - chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y); // ****************** Load Time ********************************** chain->SetBranchStatus("Time", true); chain->SetBranchAddress("Time", &Time); - vector<double> Toff13 , Toff14, Toff23, Toff24; - const char* Path13 = "../../mwpc/Toff/output/Toff13.txt"; - const char* Path14 = "../../mwpc/Toff/output/Toff14.txt"; - const char* Path23 = "../../mwpc/Toff/output/Toff23.txt"; - const char* Path24 = "../../mwpc/Toff/output/Toff24.txt"; - - Toff13 = TxtToVector(Path13); - Toff14 = TxtToVector(Path14); - Toff23 = TxtToVector(Path23); - Toff24 = TxtToVector(Path24); - - ProfileEvaluator Profile; - Profile.LoadProfile("../XYCalibration/Output/RatioProfile.root","ICOneZeroProfile"); - //***************************************************************** - TH2F *hChioDE_E_all = new TH2F("hChioDE_E","hChioDE_E",1000,1000,20000,500,1000,26000); - TH2F *hChioDE_E[Ncuts]; - for(int i=0;i<Ncuts;i++) hChioDE_E[i]=new TH2F(Form("hChioDE_E_%d",i+1),Form("hChioDE_E_%d",i+1),2000,1000,20000,1250,1000,26000); + TH2F *hChioDE_E_all = new TH2F("hChioDE_E","hChioDE_E",1000,1000,36000,500,1000,26000); + vector<TH2F*> hChioDE_E(Ncuts); + for(int i=0;i<Ncuts;i++) hChioDE_E[i]=new TH2F(Form("hChioDE_E_%d",i+1),Form("hChioDE_E_%d",i+1),2000,1000,36000,1250,1000,26000); auto start = std::chrono::high_resolution_clock::now(); - - //int Nentries=chain->GetEntries(); - int Nentries=1000000; - + + int Nentries=chain->GetEntries(); + //int Nentries=1000000; + for(int e=0;e<Nentries;++e) { if (e % 100000 == 0 && e > 0 ) { @@ -114,78 +79,11 @@ void SplineChio() chain->GetEntry(e); - vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the - vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first - // segment of ic - - - - if (Time->GetMWPC13Mult() ==1 && IC->fIC_TS.size()>=8){ //only mult 1 event - - //*************************************** GET TIME - UShort_t FPMW_Section = Time->GetSection_MWPC3(0); - double FF_DriftTime = 10* (IC->fIC_TS.at(0) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))) ; - - // ********************************* ALIGN IC*************************************************************** - for (int seg = 1; seg < IC->fIC_TS.size() ; seg++) { // loop on multiplicity of event - if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0 - - double FF_DriftTime_temp = 10* (IC->fIC_TS.at(seg) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))) ; - if (spline_Y.at(seg)==0){ - ICcorr_Y.at(seg) = IC->fIC_raw[seg]; - } - - else { - if (seg == 0) { - ICcorr_Y.at(seg) = ICcorr_Y.at(0) * spline_Y.at(0)->Eval(0)/spline_Y.at(0)->Eval(FF_DriftTime); - } - - else if (seg == 1){ - ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_DriftTime); - } - - else if (seg > 1) { - ICcorr_Y.at(seg) = IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_DriftTime); - } - if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0; - } //end if non empty - if (seg == 0) break; - }//endloop seg - - //********************* CORR IC 0 ****************************** - Double_t PolX = 1.37622; - - Double_t ICRatio = IC->fIC_raw[1]/IC->fIC_raw[0]; - Double_t ICRatioCorr = ICRatio * PolX / Profile.Evaluate(FF_IC_X,FF_DriftTime,false); - Double_t ICcorr_Y0 = IC->fIC_raw[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false); - - if ( ICRatioCorr<1.4 || ICRatioCorr >1.5){ - ICRatioCorr = ICRatio * PolX / Profile.Evaluate(FF_IC_X,FF_DriftTime,false); - ICcorr_Y0 = IC->fIC_raw[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false); - if (ICRatioCorr >100) { - ICcorr_Y0 = IC->fIC_raw[0]; - } - } - - // ************************** CORR Y GLOB ***************************** - - double FF_DE_precorr = 0.5*( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ ICcorr_Y.at(3)) + ICcorr_Y.at(4); - FF_DE = FF_DE_precorr * splineDE0234->Eval(0) / splineDE0234->Eval(FF_DriftTime); - FF_Eres = 0 ; - - for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ - if (seg >4){ - if (seg == 5) FF_Eres += double(IC->fIC_raw[seg]); - else FF_Eres += 2*double(IC->fIC_raw[seg]) ; - } - } - - hChioDE_E_all->Fill(FF_Eres,FF_DE); + hChioDE_E_all->Fill(IC->Eres,IC->DE); - //for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(Chio_E,FF_DE)) {hChioDE_E[i]->Fill(Chio_E,FF_DE); break;} - for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) {hChioDE_E[i]->Fill(FF_Eres,FF_DE); break;} - } // end max mult 1 + //for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(Chio_E,FF_DE)) {hChioDE_E[i]->Fill(Chio_E,FF_DE); break;} + for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(IC->Eres,IC->DE)) {hChioDE_E[i]->Fill(IC->Eres,IC->DE); break;} }//end loop event @@ -197,7 +95,7 @@ void SplineChio() TCanvas * can = new TCanvas(Form("ChioEbis_vs_ChioDE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]), Form("ChioEbis_vs_ChiodE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]), 0,0,2000,1000); - can->cd(); gPad-> SetLogz(); + can->cd(); hChioDE_E_all->Draw("colz"); for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same"); } @@ -205,14 +103,15 @@ void SplineChio() /////////////////////////////////////////////////////////////////////////// void MakeSpline(){ - TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root"); - TSpline3 *gspline[Ncuts]; + // ********************** Load prev histo**************************** + TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root"); - TH2F* h2 [Ncuts]; - TH1F* hProfile[Ncuts]; - TH1F* pfx[Ncuts]; // extended hProfile - TFile *fspline=new TFile("Output/spline_Chio_2024.root","recreate"); + vector<TSpline3*> gspline(Ncuts); + vector<TH2F*> h2 (Ncuts); + vector<TH1F*> hProfile(Ncuts); + vector<TH1F*> pfx(Ncuts); // extended hProfile + TFile* fspline=new TFile("Output/spline_Chio_2024.root","recreate"); TCanvas * canfit = new TCanvas("canfit","canfit",0,0,2000,1500); @@ -344,42 +243,34 @@ void ApplySpline_perTCutG() // *************************Cut*************************** TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root"); - TCutG *cutZ[Ncuts]; + vector<TCutG*> cutZ(Ncuts); // ******************* TSpline DE************************* TFile *fspline=new TFile("Output/spline_Chio_2024.root","read"); - TSpline3 *gspline[Ncuts]; + vector<TSpline3*> gspline(Ncuts); for(int i=0;i<Ncuts;i++){ cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i)); gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1)); } - // ****************** Spline Y ********************************** - TFile *fSplineIC = new TFile("../XYCalibration/Output/spline_P2P_2024.root","open"); - - - vector<vector<TSpline3*>> Spline = LoadSpline("../XYCalibration/Output/spline_P2P_2024.root"); - vector<TSpline3*> spline_X(11), spline_Y(11); - spline_X = Spline.at(0) ; spline_Y = Spline.at(1); - //**********************DE_E**************************************************** TChain *chain=new TChain("PhysicsTree"); - chain->Add("../../../root/analysis/VamosCalib241.root"); + for(int i=0;i<Nrun;i++) + { + chain->Add(Form("../../../root/analysis/Run%d.root",RunNumber[i])); + cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl; + } TICPhysics* IC = new TICPhysics() ; - double FF_IC_X, FF_IC_Y, FF_V13; - double FF_DE, FF_Eres; chain->SetBranchStatus("IC", true); chain->SetBranchAddress("IC", &IC); - chain->SetBranchStatus("FF_IC_Y", true); - chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y); - // variable after applying TSpline + // variable after applying TSpline MUST FIND middles of each spline double FF_DEcorr; - double FF_DEcorr0[Ncuts] = { + vector<double> FF_DEcorr0 = { 7375.73, 7982.83, 8419.21, @@ -421,6 +312,7 @@ TH2F *hChioDE_E = new TH2F("hChioDE_E","hChioDE_E",1000,1000,41000,500,1000 TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000); int Nentries=chain->GetEntries(); +//int Nentries = 1000000; auto start = std::chrono::high_resolution_clock::now(); for(int e=0;e<Nentries;++e){ @@ -435,55 +327,11 @@ for(int e=0;e<Nentries;++e){ chain->GetEntry(e); - - - vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the - vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first - // segment of ic - for (int seg = 1; seg < NSegment+1 ; seg++) { - - if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0 - - if (spline_Y.at(seg)==0){ - ICcorr_Y.at(seg) = IC->fIC_raw[seg]; - } - - else { - if (seg == 0) { - Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y); - ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg); - } - - else if (seg == 1){ - ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y); - } - - else if (seg > 1) { - Temp_Y.at(seg) = 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_Y.at(seg); - } - - if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0; - } //end if non empty - if (seg == 0) break; - }//endloop seg - - - FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; - FF_Eres = 0 ; - - for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ - if (seg >4){ - FF_Eres += IC->fIC_raw[seg] ; - //Ecorr += ICcorr_Y.at(seg) ; - } - } - - hChioDE_E_all->Fill(FF_Eres,FF_DE); - for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) { - hChioDE_E->Fill(FF_Eres,FF_DE); - FF_DEcorr = FF_DEcorr0[i] * FF_DE / gspline[i]->Eval(FF_Eres); - hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr); + hChioDE_E_all->Fill(IC->Eres,IC->DE); + for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(IC->Eres,IC->DE)) { + hChioDE_E->Fill(IC->Eres,IC->DE); + FF_DEcorr = FF_DEcorr0[i] * IC->DE / gspline[i]->Eval(IC->Eres); + hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr); break; } } @@ -502,39 +350,29 @@ can->cd(3); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz"); void ApplySpline() { // *************************Cut*************************** - TFile *fcut=new TFile("Cut/Cut_Z.root"); - TCutG *cutZ[Ncuts]; + TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root", "open"); + vector<TCutG*> cutZ(Ncuts); // ******************* TSpline DE************************* TFile *fspline=new TFile("Output/spline_Chio_2024.root","read"); - TSpline3 *gspline[Ncuts]; + vector<TSpline3*> gspline(Ncuts) ; for(int i=0;i<Ncuts;i++){ - cutZ[i]=(TCutG*)fcut->Get(Form("Z%i",i+1)); + cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i)); gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1)); } - // ****************** Spline Y ********************************** - TFile *fSplineIC = new TFile("../YCalibration/Output/spline_P2P_2024.root","open"); - - - vector<vector<TSpline3*>> Spline = LoadSpline("../YCalibration/Output/spline_P2P_2024.root"); - vector<TSpline3*> spline_X(11), spline_Y(11); - spline_X = Spline.at(0) ; spline_Y = Spline.at(1); - //**********************DE_E**************************************************** TChain *chain=new TChain("PhysicsTree"); - - chain->Add("../../../root/analysis/VamosCalib241.root"); + for(int i=0;i<Nrun;i++) + { + chain->Add(Form("../../../root/analysis/Run%d.root",RunNumber[i])); + cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl; + } TICPhysics* IC = new TICPhysics() ; - double FF_IC_X, FF_IC_Y, FF_V13; - double FF_DE, FF_Eres; - chain->SetBranchStatus("IC", true); chain->SetBranchAddress("IC", &IC); - chain->SetBranchStatus("FF_IC_Y", true); - chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y); // variable after applying TSpline @@ -542,7 +380,7 @@ void ApplySpline() double Zrough=0; double Zabitbetter=0; - double FF_DEcorr0[Ncuts]; + vector<double> FF_DEcorr0(Ncuts); for(int index=0; index<Ncuts; index++){ FF_DEcorr0[index] = gspline[index]->Eval(8500); } @@ -590,7 +428,8 @@ void ApplySpline() TH1F *hChioZ_rough = new TH1F("hChioZ_rough","hChioZ_rough",2500,25,65); TH1F *hChioZ_abitbetter = new TH1F("hChioZ_abitbetter","hChioZ_abitbetter",2500,25,65); - int Nentries=1e6;//chain->GetEntries(); + //int Nentries=1e6; + int Nentries = chain->GetEntries(); auto start = std::chrono::high_resolution_clock::now(); double FF_Eres_prev = 0; @@ -608,89 +447,39 @@ void ApplySpline() FF_DEcorr = -100; chain->GetEntry(e); - //=========================================================================================================== - // Splined Y - //=========================================================================================================== - vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the - vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first - // segment of ic - for (int seg = 1; seg < NSegment+1 ; seg++) { - - if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0 - - if (spline_Y.at(seg)==0){ - ICcorr_Y.at(seg) = IC->fIC_raw[seg]; - } - - else { - if (seg == 0) { - Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y); - ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg); - } - - else if (seg == 1){ - ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y); - } - - else if (seg > 1) { - Temp_Y.at(seg) = 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_Y.at(seg); - } - - if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0; - } //end if non empty - if (seg == 0) break; - }//endloop seg - - - FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4]; - FF_Eres = 0 ; + if (IC->DE<3000) continue; + if(IC->Eres==FF_Eres_prev) continue; + FF_Eres_prev = IC->Eres; - for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){ - if (seg >4){ - FF_Eres += IC->fIC_raw[seg] ; - //Ecorr += ICcorr_Y.at(seg) ; - } - } - - //=========================================================================================================== - // Spline DE - //===========================================================================================================& - - - if (FF_DE<3000) continue; - if(FF_Eres==FF_Eres_prev) continue; - FF_Eres_prev = FF_Eres; - - hChioDE_E_all->Fill(FF_Eres,FF_DE); + hChioDE_E_all->Fill(IC->Eres,IC->DE); float Eval_DEspline, DEspline0; int index=0; for(int i=0; i<Ncuts; i++){ - Eval_DEspline = gspline[i]->Eval(FF_Eres); - if(FF_DE<Eval_DEspline) break; + Eval_DEspline = gspline[i]->Eval(IC->Eres); + if(IC->DE<Eval_DEspline) break; index = i; } - Eval_DEspline = gspline[index]->Eval(FF_Eres); + Eval_DEspline = gspline[index]->Eval(IC->Eres); DEspline0 = FF_DEcorr0[index]; float dmin, dsup; - if( index < (Ncuts-1) && FF_DE > gspline[0]->Eval(FF_Eres)){ - dmin = FF_DE-gspline[index]->Eval(FF_Eres); - dsup = gspline[index+1]->Eval(FF_Eres)-FF_DE; + if( index < (Ncuts-1) && IC->DE > gspline[0]->Eval(IC->Eres)){ + dmin = IC->DE-gspline[index]->Eval(IC->Eres); + dsup = gspline[index+1]->Eval(IC->Eres)-IC->DE; if(dmin<0) cout << "negative value of dmin = " << dmin << ", for index = " << index << endl; if(dsup<0) cout << "negative value of dsup = " << dsup << ", for index = " << index << endl; //if(dsup<dmin) { - // Eval_DEspline = gspline[index+1]->Eval(FF_Eres) ; + // Eval_DEspline = gspline[index+1]->Eval(IC->Eres) ; // DEspline0 = FF_DEcorr0[index] ; //} - Eval_DEspline = dsup*gspline[index]->Eval(FF_Eres)/(dmin+dsup) + dmin*gspline[index+1]->Eval(FF_Eres)/(dmin+dsup) ; + Eval_DEspline = dsup*gspline[index]->Eval(IC->Eres)/(dmin+dsup) + dmin*gspline[index+1]->Eval(IC->Eres)/(dmin+dsup) ; DEspline0 = dsup*FF_DEcorr0[index]/(dmin+dsup) + dmin*FF_DEcorr0[index+1]/(dmin+dsup) ; - FF_DEcorr = DEspline0 * FF_DE / Eval_DEspline ; - hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr); + FF_DEcorr = DEspline0 * IC->DE / Eval_DEspline ; + hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr); //if(FF_DEcorr>15120 && FF_DEcorr<15130) { - // cout << e << " " << index << " " << FF_Eres << " " << FF_DEcorr << endl; + // cout << e << " " << index << " " << IC->Eres << " " << FF_DEcorr << endl; //} } @@ -698,9 +487,9 @@ void ApplySpline() //Zrough = -110.165 + 3.34475*sqrt(FF_DEcorr) - 0.0271123*FF_DEcorr + 8.60752e-05 * pow(sqrt(FF_DEcorr),3); Zrough = 16.8521 + 0.0017328*FF_DEcorr + 1.70774e-8*pow(FF_DEcorr,2); Zabitbetter = -127.117 + 3.83463*sqrt(FF_DEcorr) - 0.0317448 *FF_DEcorr + 0.000100428 * pow(sqrt(FF_DEcorr),3); - hChioZ_E_rough->Fill(FF_Eres,Zrough); + hChioZ_E_rough->Fill(IC->Eres,Zrough); - if(3000<FF_Eres && FF_Eres<25000){ + if(3000<IC->Eres && IC->Eres<25000){ hChioDE_corr->Fill(FF_DEcorr); hChioZ_rough->Fill(Zrough); hChioZ_abitbetter->Fill(Zabitbetter); @@ -800,66 +589,21 @@ void Fit_DE_E_corr() } +const int CutNumberLoader(){ + TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root","open"); -vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline){ - TFile *fSplineIC = new TFile(PathToSpline,"open"); - // Get number of spline - int SplineCount = 0 ; - TIter next(fSplineIC->GetListOfKeys()); + int CutCount = 0 ; + TIter next(fcut->GetListOfKeys()); TKey* key; while ((key=(TKey*)next())){ - if (std::string(key->GetClassName()) == "TSpline3"){ - SplineCount ++; + if (std::string(key->GetClassName()) == "TCutG"){ + CutCount ++; } } - vector<vector<TSpline3*>> Spline; - vector<TSpline3*> spline_X(11), spline_Y(11); - for (int i = 0; i < SplineCount ; i++) { - int seg = int((i-2)/2 +1 ); - if (i < 2){ - spline_X.at(0) = (TSpline3 *)fSplineIC->Get("fspline_IC0_X"); - spline_Y.at(0) = (TSpline3 *)fSplineIC->Get("fspline_IC0_Y"); - } - - else if ( i>=2 && i<4){ - spline_X.at(1) = (TSpline3 *)fSplineIC->Get("fspline_IC1_X"); - spline_Y.at(1) = (TSpline3 *)fSplineIC->Get("fspline_IC1_Y"); - } - else if (seg >=2 && (i%2 == 0) ) { - spline_X.at(seg) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_X",seg)); - } - else if (seg >=2 && !(i%2 == 0) ) { - spline_Y.at(seg) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_Y",seg)); - } - } //End loop on histogram - - Spline.push_back(spline_X); - Spline.push_back(spline_Y); - - return Spline; + const int res = CutCount; + return res; } - -vector<double> TxtToVector(const char *Path){ - string line; - vector<double> values; - ifstream file(Path); - - if (file.is_open()) { - while (std::getline(file, line)) { - try { - values.push_back(std::stod(line)); - } catch (const std::invalid_argument& e) { - std::cerr << "Invalid number in line: " << line << '\n'; - } - } - file.close(); - } else { - std::cerr << "Error opening file.\n"; - } - - return values; -}