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

[AlPhaPha] 1D plot of charge extracted

parent eef51f28
No related branches found
No related tags found
No related merge requests found
Pipeline #384547 passed
...@@ -8,19 +8,18 @@ ...@@ -8,19 +8,18 @@
#include <TH2.h> #include <TH2.h>
#include <TSpline.h> #include <TSpline.h>
#include <fstream> #include <fstream>
#include "../XYCalibration/ProfileEvaluator.h"
using namespace std; using namespace std;
const int Ncuts=30;
const int CutNumberLoader();
const int Ncuts = CutNumberLoader();
const int NSegment=11; const int NSegment=11;
const int Nrun=1; const int Nrun=1;
int RunNumber[Nrun]={241}; int RunNumber[Nrun]={247};
void MakeSpline(); void MakeSpline();
void ApplySpline_perTCutG(); void ApplySpline_perTCutG();
void ApplySpline(); void ApplySpline();
vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline);
vector<double> TxtToVector(const char *Path);
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
void SplineChio() void SplineChio()
...@@ -31,76 +30,42 @@ void SplineChio() ...@@ -31,76 +30,42 @@ void SplineChio()
//**********************Cut************************************** //**********************Cut**************************************
TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root","open"); 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**************************************************** //**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree"); 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() ; TICPhysics* IC = new TICPhysics() ;
TTimeData *Time = new TTimeData() ; 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->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC); chain->SetBranchAddress("IC", &IC);
chain->SetBranchStatus("FF_IC_Y", true);
chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y);
// ****************** Load Time ********************************** // ****************** Load Time **********************************
chain->SetBranchStatus("Time", true); chain->SetBranchStatus("Time", true);
chain->SetBranchAddress("Time", &Time); 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_all = new TH2F("hChioDE_E","hChioDE_E",1000,1000,36000,500,1000,26000);
TH2F *hChioDE_E[Ncuts]; 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,20000,1250,1000,26000); 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(); auto start = std::chrono::high_resolution_clock::now();
//int Nentries=chain->GetEntries(); int Nentries=chain->GetEntries();
int Nentries=1000000; //int Nentries=1000000;
for(int e=0;e<Nentries;++e) for(int e=0;e<Nentries;++e)
{ {
if (e % 100000 == 0 && e > 0 ) { if (e % 100000 == 0 && e > 0 ) {
...@@ -114,78 +79,11 @@ void SplineChio() ...@@ -114,78 +79,11 @@ void SplineChio()
chain->GetEntry(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
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(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;} 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 max mult 1
}//end loop event }//end loop event
...@@ -197,7 +95,7 @@ void SplineChio() ...@@ -197,7 +95,7 @@ void SplineChio()
TCanvas * can = new TCanvas(Form("ChioEbis_vs_ChioDE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]), 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]), Form("ChioEbis_vs_ChiodE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
0,0,2000,1000); 0,0,2000,1000);
can->cd(); gPad-> SetLogz(); can->cd();
hChioDE_E_all->Draw("colz"); hChioDE_E_all->Draw("colz");
for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same"); for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
} }
...@@ -205,14 +103,15 @@ void SplineChio() ...@@ -205,14 +103,15 @@ void SplineChio()
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
void MakeSpline(){ 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]; vector<TSpline3*> gspline(Ncuts);
TH1F* hProfile[Ncuts]; vector<TH2F*> h2 (Ncuts);
TH1F* pfx[Ncuts]; // extended hProfile vector<TH1F*> hProfile(Ncuts);
TFile *fspline=new TFile("Output/spline_Chio_2024.root","recreate"); 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); TCanvas * canfit = new TCanvas("canfit","canfit",0,0,2000,1500);
...@@ -344,42 +243,34 @@ void ApplySpline_perTCutG() ...@@ -344,42 +243,34 @@ void ApplySpline_perTCutG()
// *************************Cut*************************** // *************************Cut***************************
TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root"); TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root");
TCutG *cutZ[Ncuts]; vector<TCutG*> cutZ(Ncuts);
// ******************* TSpline DE************************* // ******************* TSpline DE*************************
TFile *fspline=new TFile("Output/spline_Chio_2024.root","read"); TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
TSpline3 *gspline[Ncuts]; vector<TSpline3*> gspline(Ncuts);
for(int i=0;i<Ncuts;i++){ for(int i=0;i<Ncuts;i++){
cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i)); cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1)); 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**************************************************** //**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree"); 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() ; TICPhysics* IC = new TICPhysics() ;
double FF_IC_X, FF_IC_Y, FF_V13;
double FF_DE, FF_Eres;
chain->SetBranchStatus("IC", true); chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC); 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_DEcorr;
double FF_DEcorr0[Ncuts] = { vector<double> FF_DEcorr0 = {
7375.73, 7375.73,
7982.83, 7982.83,
8419.21, 8419.21,
...@@ -421,6 +312,7 @@ TH2F *hChioDE_E = new TH2F("hChioDE_E","hChioDE_E",1000,1000,41000,500,1000 ...@@ -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); TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);
int Nentries=chain->GetEntries(); int Nentries=chain->GetEntries();
//int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now(); auto start = std::chrono::high_resolution_clock::now();
for(int e=0;e<Nentries;++e){ for(int e=0;e<Nentries;++e){
...@@ -435,55 +327,11 @@ for(int e=0;e<Nentries;++e){ ...@@ -435,55 +327,11 @@ for(int e=0;e<Nentries;++e){
chain->GetEntry(e); chain->GetEntry(e);
hChioDE_E_all->Fill(IC->Eres,IC->DE);
for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(IC->Eres,IC->DE)) {
vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the hChioDE_E->Fill(IC->Eres,IC->DE);
vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first FF_DEcorr = FF_DEcorr0[i] * IC->DE / gspline[i]->Eval(IC->Eres);
// segment of ic hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr);
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);
break; break;
} }
} }
...@@ -502,39 +350,29 @@ can->cd(3); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz"); ...@@ -502,39 +350,29 @@ can->cd(3); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz");
void ApplySpline() void ApplySpline()
{ {
// *************************Cut*************************** // *************************Cut***************************
TFile *fcut=new TFile("Cut/Cut_Z.root"); TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root", "open");
TCutG *cutZ[Ncuts]; vector<TCutG*> cutZ(Ncuts);
// ******************* TSpline DE************************* // ******************* TSpline DE*************************
TFile *fspline=new TFile("Output/spline_Chio_2024.root","read"); TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
TSpline3 *gspline[Ncuts]; vector<TSpline3*> gspline(Ncuts) ;
for(int i=0;i<Ncuts;i++){ 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)); 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**************************************************** //**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree"); TChain *chain=new TChain("PhysicsTree");
for(int i=0;i<Nrun;i++)
chain->Add("../../../root/analysis/VamosCalib241.root"); {
chain->Add(Form("../../../root/analysis/Run%d.root",RunNumber[i]));
cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
}
TICPhysics* IC = new TICPhysics() ; TICPhysics* IC = new TICPhysics() ;
double FF_IC_X, FF_IC_Y, FF_V13;
double FF_DE, FF_Eres;
chain->SetBranchStatus("IC", true); chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC); 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
...@@ -542,7 +380,7 @@ void ApplySpline() ...@@ -542,7 +380,7 @@ void ApplySpline()
double Zrough=0; double Zrough=0;
double Zabitbetter=0; double Zabitbetter=0;
double FF_DEcorr0[Ncuts]; vector<double> FF_DEcorr0(Ncuts);
for(int index=0; index<Ncuts; index++){ for(int index=0; index<Ncuts; index++){
FF_DEcorr0[index] = gspline[index]->Eval(8500); FF_DEcorr0[index] = gspline[index]->Eval(8500);
} }
...@@ -590,7 +428,8 @@ void ApplySpline() ...@@ -590,7 +428,8 @@ void ApplySpline()
TH1F *hChioZ_rough = new TH1F("hChioZ_rough","hChioZ_rough",2500,25,65); TH1F *hChioZ_rough = new TH1F("hChioZ_rough","hChioZ_rough",2500,25,65);
TH1F *hChioZ_abitbetter = new TH1F("hChioZ_abitbetter","hChioZ_abitbetter",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(); auto start = std::chrono::high_resolution_clock::now();
double FF_Eres_prev = 0; double FF_Eres_prev = 0;
...@@ -608,89 +447,39 @@ void ApplySpline() ...@@ -608,89 +447,39 @@ void ApplySpline()
FF_DEcorr = -100; FF_DEcorr = -100;
chain->GetEntry(e); chain->GetEntry(e);
//=========================================================================================================== if (IC->DE<3000) continue;
// Splined Y if(IC->Eres==FF_Eres_prev) continue;
//=========================================================================================================== FF_Eres_prev = IC->Eres;
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++ ){ hChioDE_E_all->Fill(IC->Eres,IC->DE);
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);
float Eval_DEspline, DEspline0; float Eval_DEspline, DEspline0;
int index=0; int index=0;
for(int i=0; i<Ncuts; i++){ for(int i=0; i<Ncuts; i++){
Eval_DEspline = gspline[i]->Eval(FF_Eres); Eval_DEspline = gspline[i]->Eval(IC->Eres);
if(FF_DE<Eval_DEspline) break; if(IC->DE<Eval_DEspline) break;
index = i; index = i;
} }
Eval_DEspline = gspline[index]->Eval(FF_Eres); Eval_DEspline = gspline[index]->Eval(IC->Eres);
DEspline0 = FF_DEcorr0[index]; DEspline0 = FF_DEcorr0[index];
float dmin, dsup; float dmin, dsup;
if( index < (Ncuts-1) && FF_DE > gspline[0]->Eval(FF_Eres)){ if( index < (Ncuts-1) && IC->DE > gspline[0]->Eval(IC->Eres)){
dmin = FF_DE-gspline[index]->Eval(FF_Eres); dmin = IC->DE-gspline[index]->Eval(IC->Eres);
dsup = gspline[index+1]->Eval(FF_Eres)-FF_DE; dsup = gspline[index+1]->Eval(IC->Eres)-IC->DE;
if(dmin<0) cout << "negative value of dmin = " << dmin << ", for index = " << index << endl; 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<0) cout << "negative value of dsup = " << dsup << ", for index = " << index << endl;
//if(dsup<dmin) { //if(dsup<dmin) {
// Eval_DEspline = gspline[index+1]->Eval(FF_Eres) ; // Eval_DEspline = gspline[index+1]->Eval(IC->Eres) ;
// DEspline0 = FF_DEcorr0[index] ; // 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) ; DEspline0 = dsup*FF_DEcorr0[index]/(dmin+dsup) + dmin*FF_DEcorr0[index+1]/(dmin+dsup) ;
FF_DEcorr = DEspline0 * FF_DE / Eval_DEspline ; FF_DEcorr = DEspline0 * IC->DE / Eval_DEspline ;
hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr); hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr);
//if(FF_DEcorr>15120 && FF_DEcorr<15130) { //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() ...@@ -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 = -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); 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); 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); hChioDE_corr->Fill(FF_DEcorr);
hChioZ_rough->Fill(Zrough); hChioZ_rough->Fill(Zrough);
hChioZ_abitbetter->Fill(Zabitbetter); hChioZ_abitbetter->Fill(Zabitbetter);
...@@ -800,66 +589,21 @@ void Fit_DE_E_corr() ...@@ -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){ int CutCount = 0 ;
TFile *fSplineIC = new TFile(PathToSpline,"open"); TIter next(fcut->GetListOfKeys());
// Get number of spline
int SplineCount = 0 ;
TIter next(fSplineIC->GetListOfKeys());
TKey* key; TKey* key;
while ((key=(TKey*)next())){ while ((key=(TKey*)next())){
if (std::string(key->GetClassName()) == "TSpline3"){ if (std::string(key->GetClassName()) == "TCutG"){
SplineCount ++; CutCount ++;
} }
} }
vector<vector<TSpline3*>> Spline; const int res = CutCount;
vector<TSpline3*> spline_X(11), spline_Y(11); return res;
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;
} }
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;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment