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

Continuing analysis, tried IC1 not corrected and only correcting the DE

parent 3a794d7a
No related branches found
No related tags found
1 merge request!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...
Pipeline #366030 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Showing
with 637 additions and 114 deletions
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorr.C 0 → 100644
+ 253
− 0
View file @ 1719af04
//#include <TICPhysics.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TF1.h>
#include <TFile.h>
#include <TH2.h>
#include <TSpline.h>
#include <TCutG.h>
#include <fstream>
#include <vector>
using namespace std;
int GetNumberKey(TFile *infile ,const char* ClassName);
TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double
XMax, double YMin, double YMax);
//////////////////////////////////////////////////////////////////////////////////////////
void DECorr(bool cut = false, bool spline = false) {
//===========================================================================================================
// Loading var
//===========================================================================================================
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);
TCutG *CutZ;
TCutG *CutZbis;
TFile *fCut;
if (cut) {
fCut = TFile::Open("Output/CutDeCorr.root");
CutZbis = (TCutG*)fCut->Get("CutZbis") ;
CutZ = (TCutG*)fCut->Get("CutZbis") ;
}
TFile *fSpline;
TSpline3 *splineDEcorr, *splineDEbis;
if(spline){
fSpline= new TFile("Output/splineDE.root","open");
splineDEcorr = (TSpline3*) fSpline->Get("SplineDe");
splineDEbis = (TSpline3*) fSpline->Get("SplineDebis");
}
//Defining output histo
TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,100,1,1000,100,1);
TH2F* hDE_Ebis = new TH2F("DE_Ebis","DE_Ebis",1000,100,1,1000,100,1);
TH2F* hDE_Y = new TH2F("DE_Y","DE_Y",1000,-100,100,1000,0,22000);
TH2F* hDE_Y_bis = new TH2F("DE_Ybis","DE_Ybis",1000,-100,100,1000,0,22000);
TH2F* hDE_E_corr = new TH2F("DE_Ecorr","DE_Ecorr",1000,100,1,1000,0,22000);
TH2F* hDE_Ebis_corr = new TH2F("DE_Ebiscorr","DE_Ebiscorr",1000,100,1,1000,0,22000);
//===========================================================================================================
// Event Loop
//===========================================================================================================
int Nentries = chain->GetEntries();
//int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();
for (int e = 0; e < Nentries; e++) {
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 << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
}
chain->GetEntry(e);
double DE = 0 ,E =0 ,DE_Ycorr=0, DE_Ybis=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];
if (FF_IC_Y >-100 && FF_IC_Y <100){
hDE_E->Fill(E,DE);
hDE_Ebis->Fill(E,DE_Bis);
if (spline== true){
DE_Ycorr = DE * 13750 / splineDEcorr->Eval(FF_IC_Y);
hDE_E_corr->Fill(E,DE_Ycorr);
DE_Ybis = DE_Bis * 13750 / splineDEbis->Eval(FF_IC_Y);
hDE_Ebis_corr->Fill(E,DE_Ybis);
}
if (cut == true && !(CutZ->IsInside(E,DE))) {;}
else {
hDE_Y->Fill(FF_IC_Y,DE);
}//end cut
if (cut == true && !(CutZbis->IsInside(E,DE_Bis))) {;}
else {
hDE_Y_bis->Fill(FF_IC_Y,DE_Bis);
}//end cut
}//end cond Y
} //endl loop event
// making spline
TSpline3 *splineDE, *osplineDEbis;
if (cut == true){
splineDE = MakeSpline(hDE_Y,0,-100,100,14300,15000);
osplineDEbis = MakeSpline(hDE_Y_bis,1,-100,100,14300,15000);
TFile *fSpline = new TFile("Output/splineDE.root","recreate");
splineDE->SetName("SplineDe");
splineDE->Write();
osplineDEbis->SetName("SplineDebis");
osplineDEbis->Write();
fSpline->Close();
fCut->cd();
}
//===========================================================================================================
// Drawing histo
//===========================================================================================================
TCanvas* c2 = new TCanvas("c2","c2",1500,1000);
c2->Divide(2 + 2*spline);
c2->cd(1);
gPad->SetLogz();
hDE_E->Draw();
c2->cd(2);
gPad->SetLogz();
hDE_Ebis->Draw();
if (spline){
c2->cd(3);
gPad->SetLogz();
hDE_E_corr->Draw();
c2->cd(4);
gPad->SetLogz();
hDE_Ebis_corr->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;
}
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;
}
//===========================================================================================================
// 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
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/DECorrParra.C 0 → 100644
+ 253
− 0
View file @ 1719af04
//#include <TICPhysics.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TF1.h>
#include <TFile.h>
#include <TH2.h>
#include <TSpline.h>
#include <TCutG.h>
#include <fstream>
#include <vector>
using namespace std;
int GetNumberKey(TFile *infile ,const char* ClassName);
TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double
XMax, double YMin, double YMax);
//////////////////////////////////////////////////////////////////////////////////////////
void DECorrParra(bool cut = false, bool spline = false) {
//===========================================================================================================
// Loading var
//===========================================================================================================
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);
TCutG *CutZ;
TCutG *CutZbis;
TFile *fCut;
if (cut) {
fCut = TFile::Open("Output/CutDeCorr.root");
CutZbis = (TCutG*)fCut->Get("CutZbis") ;
CutZ = (TCutG*)fCut->Get("CutZbis") ;
}
TFile *fSpline;
TSpline3 *splineDEcorr, *splineDEbis;
if(spline){
fSpline= new TFile("Output/splineDE.root","open");
splineDEcorr = (TSpline3*) fSpline->Get("SplineDe");
splineDEbis = (TSpline3*) fSpline->Get("SplineDebis");
}
//Defining output histo
TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,100,1,1000,100,1);
TH2F* hDE_Ebis = new TH2F("DE_Ebis","DE_Ebis",1000,100,1,1000,100,1);
TH2F* hDE_Y = new TH2F("DE_Y","DE_Y",1000,-100,100,1000,0,22000);
TH2F* hDE_Y_bis = new TH2F("DE_Ybis","DE_Ybis",1000,-100,100,1000,0,22000);
TH2F* hDE_E_corr = new TH2F("DE_Ecorr","DE_Ecorr",1000,100,1,1000,0,22000);
TH2F* hDE_Ebis_corr = new TH2F("DE_Ebiscorr","DE_Ebiscorr",1000,100,1,1000,0,22000);
//===========================================================================================================
// Event Loop
//===========================================================================================================
int Nentries = chain->GetEntries();
//int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();
for (int e = 0; e < Nentries; e++) {
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 << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
}
chain->GetEntry(e);
double DE = 0 ,E =0 ,DE_Ycorr=0, DE_Ybis=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];
if (FF_IC_Y >-100 && FF_IC_Y <100){
hDE_E->Fill(E,DE);
hDE_Ebis->Fill(E,DE_Bis);
if (spline== true){
DE_Ycorr = DE * 13750 / splineDEcorr->Eval(FF_IC_Y);
hDE_E_corr->Fill(E,DE_Ycorr);
DE_Ybis = DE_Bis * 13750 / splineDEbis->Eval(FF_IC_Y);
hDE_Ebis_corr->Fill(E,DE_Ybis);
}
if (cut == true && !(CutZ->IsInside(E,DE))) {;}
else {
hDE_Y->Fill(FF_IC_Y,DE);
}//end cut
if (cut == true && !(CutZbis->IsInside(E,DE_Bis))) {;}
else {
hDE_Y_bis->Fill(FF_IC_Y,DE_Bis);
}//end cut
}//end cond Y
} //endl loop event
// making spline
TSpline3 *splineDE, *osplineDEbis;
if (cut == true){
splineDE = MakeSpline(hDE_Y,0,-100,100,14300,15000);
osplineDEbis = MakeSpline(hDE_Y_bis,1,-100,100,14300,15000);
TFile *fSpline = new TFile("Output/splineDE.root","recreate");
splineDE->SetName("SplineDe");
splineDE->Write();
osplineDEbis->SetName("SplineDebis");
osplineDEbis->Write();
fSpline->Close();
fCut->cd();
}
//===========================================================================================================
// Drawing histo
//===========================================================================================================
TCanvas* c2 = new TCanvas("c2","c2",1500,1000);
c2->Divide(2 + 2*spline);
c2->cd(1);
gPad->SetLogz();
hDE_E->Draw();
c2->cd(2);
gPad->SetLogz();
hDE_Ebis->Draw();
if (spline){
c2->cd(3);
gPad->SetLogz();
hDE_E_corr->Draw();
c2->cd(4);
gPad->SetLogz();
hDE_Ebis_corr->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;
}
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;
}
//===========================================================================================================
// 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
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
+ 108
− 97
View file @ 1719af04
#include <TICPhysics.h>
//#include <TICPhysics.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TF1.h>
......@@ -54,13 +54,13 @@ void ICCorr() {
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)));
hICX.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_IC%d_X",seg+1,seg)));
hICY.push_back((TH2F*)fHisto->Get(Form("hChio_IC%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)));
hICX.push_back((TH2F*)fHisto->Get(Form("hChio_IC%d_IC%d_X",seg,seg-1)));
hICY.push_back((TH2F*)fHisto->Get(Form("hChio_IC%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));
}
}
......@@ -116,108 +116,119 @@ void ICCorr() {
//===========================================================================================================
// Event Loop
//===========================================================================================================
int Nentries = chain->GetEntries();
int Nentries = chain->GetEntries();
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 ..\r" << flush;
}
chain->GetEntry(e);
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 << " | Estimated time left: " << timeLeft << " seconds" << "\r" << flush;
}
chain->GetEntry(e);
vector<double> ICcorr_Y(5), ICcorr_X(11); // the [0] element of spline is the
vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // 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) {
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();
if (spline_Y.at(seg)==0){
ICcorr_Y.at(seg) = IC->fIC_raw[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);
}
if (seg == 0) {
double temp = ICcorr_Y.at(1) / IC->fIC_raw[0] * spline_Y.at(0)->Eval(0)/ spline_Y.at(0)->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) {
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])+IC->fIC_raw[4];
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(NSegment,2);
for (int i = 0 ; i<NSegment ;i++){
c1->cd(i+1);
hICcorr.at(i)->ProfileX()->Draw();
c1->cd(NSegment+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/IC_Dampen_Recombination.h
+ 1
− 1
View file @ 1719af04
#include <TFile.h>
#include <TICPhysics.h>
//#include <TICPhysics.h>
#include <TChain.h>
#include <TF1.h>
#include <TH2.h>
......
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ReadMe.md
+ 4
− 0
View file @ 1719af04
......@@ -15,3 +15,7 @@ of the electron by the grid.
SplineXY.C serve the purpose of creating this spline and storing int a .root
file.
## Devel
Apparently the main correction must come from the correction of the DE
Projects/AlPhaPha/2024/macro/chio/CalibrationChio/SplineChioP2P.C
+ 18
− 16
View file @ 1719af04
//#ifdef TICPhysics
#ifdef TICPhysics
#include <TICPhysics.h>
//#endif
#endif
#include <TCanvas.h>
#include <TChain.h>
#include <TF1.h>
......@@ -39,21 +39,21 @@ void SplineChioP2P() {
int seg = 1;
vector<vector<TSpline3*>> SplineAllIC(11);
vector<TSpline3*> SplineICurrent;
vector<TSpline3*> SplineICurrent(2);
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(0) = MakeSpline(hIC[seg][0] , NCallSpline, XMin[0] , XMax[0], YMin[0], YMax[0]);
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(1) = MakeSpline(hIC[seg][1] , NCallSpline, XMin[1] , XMax[1], YMin[1], YMax[1]);
SplineICurrent.at(1)->SetName("fsplineIC1_Y");
NCallSpline+=1;
SplineAllIC.at(1)= (SplineICurrent); //
//SplineAllIC.erase(SplineAllIC.begin() + 1); //
for (int seg = 0 ; seg < NSegment ; seg++){
......@@ -87,16 +87,17 @@ void SplineChioP2P() {
//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();
if (!SplineAllIC.at(seg).empty()){
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);
......@@ -177,8 +178,8 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
// Beginning loop on entries
//===========================================================================================================
int Nentries = chain->GetEntries();
//int Nentries = 100000;
//int Nentries = chain->GetEntries();
int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();
for (int e = 0; e < Nentries; e++) {
......@@ -318,8 +319,8 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
// Beginning loop on entries
//===========================================================================================================
int Nentries = chain->GetEntries();
//int Nentries = 100000;
//int Nentries = chain->GetEntries();
int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();
for (int e = 0; e < Nentries; e++) {
......@@ -333,6 +334,7 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
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];
//cout << pseg << endl;
}
//If it exist correct the current line
......@@ -347,7 +349,7 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
} //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);
......
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