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

[AlPhaPha] Finished Step5

parent af5bcd6d
No related branches found
No related tags found
No related merge requests found
Pipeline #385644 canceled
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FillHistoCharge.C
+ 51
− 126
View file @ 4dfb2146
#include "TICPhysics.h"
#include <TChain.h>
#include <regex>
TChain* chain; TChain* chain;
TCutG* cut10Be[8];
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
void LoadCut(string part){ void FillHistoCharge(){
string path = "../Mass/cut/"+ part + "/cut10Be.root";
TFile* fcut = new TFile(path.c_str()); const char * Path = "../../../../root/analysis/Run247ZSpline.root";
for(int i=0; i<8; i++){ chain = new TChain("PhysicsTree");
cut10Be[i] = (TCutG*)fcut->Get(Form("cut10Be_det%i",i+1)); chain->Add(Path);
}
} TICPhysics* IC = new TICPhysics() ;
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
/////////////////////////////////////////////////////////////////
void FillHistoCharge(string part="part1"){
LoadCut(part);
chain = new TChain("PhysicsTree");
if(part=="part1"){
//chain->Add("../../root/analysis/run_28.root");
//chain->Add("../../root/analysis/run_29.root");
//chain->Add("../../root/analysis/run_30.root");
//chain->Add("../../root/analysis/run_32.root");
//chain->Add("../../root/analysis/run_36.root");
//chain->Add("../../root/analysis/run_37.root");
//chain->Add("../../root/analysis/run_39.root");
//chain->Add("../../root/analysis/run_41.root");
//chain->Add("../../root/analysis/run_42.root");
//chain->Add("../../root/analysis/run_43.root");
//chain->Add("../../root/analysis/run_44.root");
//chain->Add("../../root/analysis/run_45.root");
//chain->Add("../../root/analysis/run_46.root");
chain->Add("../../root/analysis/run_48.root");
//chain->Add("../../root/analysis/run_49.root");
//chain->Add("../../root/analysis/run_50.root");
//chain->Add("../../root/analysis/run_51.root");
//chain->Add("../../root/analysis/run_52.root");
//chain->Add("../../root/analysis/run_53.root");
}
else if(part=="part2"){
chain->Add("../../root/analysis/run_55.root");
chain->Add("../../root/analysis/run_56.root");
chain->Add("../../root/analysis/run_57.root");
chain->Add("../../root/analysis/run_58.root");
chain->Add("../../root/analysis/run_59.root");
chain->Add("../../root/analysis/run_60.root");
chain->Add("../../root/analysis/run_61.root");
chain->Add("../../root/analysis/run_62.root");
chain->Add("../../root/analysis/run_63.root");
chain->Add("../../root/analysis/run_64.root");
chain->Add("../../root/analysis/run_65.root");
chain->Add("../../root/analysis/run_66.root");
chain->Add("../../root/analysis/run_67.root");
chain->Add("../../root/analysis/run_68.root");
chain->Add("../../root/analysis/run_69.root");
}
double DeltaEcorr;
double Elab;
double Mass;
double Ex240Pu;
double VAMOS_TS_hour;
int Telescope;
int m_2alpha;
int FPMW_Section;
double FF_Z;
chain->SetBranchStatus("DeltaEcorr","true");
chain->SetBranchAddress("DeltaEcorr",&DeltaEcorr);
chain->SetBranchStatus("Elab","true");
chain->SetBranchAddress("Elab",&Elab);
chain->SetBranchStatus("FF_Mass13","true");
chain->SetBranchAddress("FF_Mass13",&Mass);
chain->SetBranchStatus("Ex240Pu","true");
chain->SetBranchAddress("Ex240Pu",&Ex240Pu);
chain->SetBranchStatus("VAMOS_TS_hour","true");
chain->SetBranchAddress("VAMOS_TS_hour",&VAMOS_TS_hour);
chain->SetBranchStatus("Telescope","true");
chain->SetBranchAddress("Telescope",&Telescope);
chain->SetBranchStatus("m_2alpha","true");
chain->SetBranchAddress("m_2alpha",&m_2alpha);
chain->SetBranchStatus("FPMW_Section","true");
chain->SetBranchAddress("FPMW_Section",&FPMW_Section);
chain->SetBranchStatus("FF_Z","true");
chain->SetBranchAddress("FF_Z",&FF_Z);
TH1F *hcharge[11];
for(int i=0; i<11; i++){
TString histo_name = Form("hcharge_%iMeV",i+6);
hcharge[i] = new TH1F(histo_name, histo_name,1000,20,80);
}
TH1F* hChargeAll = new TH1F("hcharge_all","hcharge_all",1000,20,80);
TH2D* hChargeVsEx = new TH2D("hcharge_vs_ex","hcharge_vs_ex",20,0,20,1000,20,80);
int nentries = chain->GetEntries();
cout << "**** Number of entries to be treated: " << nentries << " ****" << endl;
for(int i=0; i<nentries; i++){
chain->GetEntry(i);
if(i%100000==0){
cout << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush;
}
if(VAMOS_TS_hour>0){ TH1F *hChargeAll = new TH1F("hcharge_all","hcharge_all",2000,20,80);
hChargeAll->Fill(FF_Z); int Nentries = chain->GetEntries();
if(Telescope>0) { auto start = std::chrono::high_resolution_clock::now();
if(FPMW_Section>4 && cut10Be[Telescope-1]->IsInside(Elab,DeltaEcorr)){ for (int e = 0; e < Nentries; e++) {
hChargeVsEx->Fill(Ex240Pu,FF_Z); if (e % 100000 == 0 && e > 0 ) {
if(Ex240Pu>5.5 && Ex240Pu<6.5) hcharge[0]->Fill(FF_Z); auto now = std::chrono::high_resolution_clock::now();
else if(Ex240Pu>6.5 && Ex240Pu<7.5) hcharge[1]->Fill(FF_Z); std::chrono::duration<double> elapsed = now - start;
else if(Ex240Pu>7.5 && Ex240Pu<8.5) hcharge[2]->Fill(FF_Z); double avgTimePerIteration = elapsed.count() / e;
else if(Ex240Pu>8.5 && Ex240Pu<9.5) hcharge[3]->Fill(FF_Z); double timeLeft = avgTimePerIteration * (Nentries - e);
else if(Ex240Pu>9.5 && Ex240Pu<10.5) hcharge[4]->Fill(FF_Z);
else if(Ex240Pu>10.5 && Ex240Pu<11.5) hcharge[5]->Fill(FF_Z); std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
else if(Ex240Pu>11.5 && Ex240Pu<12.5) hcharge[6]->Fill(FF_Z);
else if(Ex240Pu>12.5 && Ex240Pu<13.5) hcharge[7]->Fill(FF_Z);
else if(Ex240Pu>13.5 && Ex240Pu<14.5) hcharge[8]->Fill(FF_Z);
else if(Ex240Pu>14.5 && Ex240Pu<15.5) hcharge[9]->Fill(FF_Z);
else if(Ex240Pu>15.5 && Ex240Pu<16.5) hcharge[10]->Fill(FF_Z);
} }
}
chain->GetEntry(e);
hChargeAll->Fill(IC->Chio_Z);
} }
}
TString filename = Form("histo_charge_%s.root",part.c_str());
TFile *ofile = new TFile(filename,"recreate");
for(int i=0; i<11; i++){
hcharge[i]->Write();
}
hChargeAll->Write();
hChargeVsEx->Write();
ofile->Close();
std::string strPath(Path);
std::regex pattern("Run(\\d{3})"); // Regex to match Run followed by exactly 3 digits
std::smatch matches;
if (std::regex_search(strPath, matches, pattern)) {
std::string number = matches[1]; // Extract the number
std::string fileName = "Output/histo_charge_" + number + ".root"; // Generate the filename
// Create the TFile
TFile* file = new TFile(fileName.c_str(), "RECREATE");
if (file->IsOpen()) {
std::cout << "File created: " << fileName << std::endl;
} else {
std::cout << "Failed to create the file." << std::endl;
}
hChargeAll->Write();
// Don't forget to close the file
file->Close();
} else {
std::cout << "Pattern not found!" << std::endl;
}
} }
Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FitCharge.C
+ 17
− 18
View file @ 4dfb2146
#include <TCanvas.h>
#include <TF1.h>
#include <TFile.h>
#include <TGraphErrors.h>
#include <TH1.h>
using namespace std;
void Fit(TH1F* hcharge, string Energy); void Fit(TH1F* hcharge, string Energy);
TFile* ofile; TFile* ofile;
int rebin=1; int rebin=1;
int RunNumber = 247 ;
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
void FitCharge(){ void FitCharge(){
TFile* ifile = new TFile("histo_charge_part1.root","read"); TFile* ifile = new TFile(Form("Output/histo_charge_%d.root",RunNumber),"read");
ofile = new TFile("histo_charge_fitted.root","recreate"); ofile = new TFile("Output/histo_charge_fitted.root","recreate");
TH1F* hall = (TH1F*)ifile->FindObjectAny("hcharge_all"); TH1F* hall = (TH1F*)ifile->FindObjectAny("hcharge_all");
Fit(hall,"AllEnergy"); Fit(hall,"AllEnergy");
/*for(int i=0; i<11; i++){
int Ex = i+6;
if(i>4) rebin=2;
else rebin=1;
string Energy = to_string(Ex) + "MeV";
TString histo_name = Form("hcharge_%iMeV",i+6);
TH1F* h1 = (TH1F*)ifile->FindObjectAny(histo_name);
Fit(h1,Energy);
}*/
ofile->Close(); ofile->Close();
} }
...@@ -31,7 +29,7 @@ void Fit(TH1F* hcharge,string Energy){ ...@@ -31,7 +29,7 @@ void Fit(TH1F* hcharge,string Energy){
hcharge->Rebin(rebin); hcharge->Rebin(rebin);
hcharge->Draw(); hcharge->Draw();
int Zmin = 32; int Zmin = 30;
int Zmax = 65; int Zmax = 65;
...@@ -40,11 +38,11 @@ void Fit(TH1F* hcharge,string Energy){ ...@@ -40,11 +38,11 @@ void Fit(TH1F* hcharge,string Energy){
int NumberOfZ = Zmax - Zmin; int NumberOfZ = Zmax - Zmin;
Double_t para[3*NumberOfZ]; Double_t para[3*NumberOfZ];
TF1* g[NumberOfZ]; vector<TF1*> g(NumberOfZ);
double Integral[NumberOfZ]; vector<double> Integral(NumberOfZ);
double Integral_err[NumberOfZ]; vector<double> Integral_err(NumberOfZ);
double total_integral = 0; double total_integral = 0;
double Yield[NumberOfZ]; vector<double> Yield(NumberOfZ);
for(int i=0; i<NumberOfZ; i++){ for(int i=0; i<NumberOfZ; i++){
g[i] = new TF1(Form("g%i",i),"gaus",Z-0.4,Z+0.4); g[i] = new TF1(Form("g%i",i),"gaus",Z-0.4,Z+0.4);
...@@ -124,6 +122,7 @@ void Fit(TH1F* hcharge,string Energy){ ...@@ -124,6 +122,7 @@ void Fit(TH1F* hcharge,string Energy){
//cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl; //cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl;
//cout << "Mean= " << mean << endl; //cout << "Mean= " << mean << endl;
cout << "Sigma= " << sigma << "+/-" << sigma_err << endl; cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;
cout << "Sigma/Z = " << 2.35*sigma /Z*100. << endl;
gsigma->SetPoint(i,Z,sigma); gsigma->SetPoint(i,Z,sigma);
gsigma->SetPointError(i,0,sigma_err); gsigma->SetPointError(i,0,sigma_err);
......
Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/SplineChio.C
+ 311
− 53
View file @ 4dfb2146
...@@ -6,6 +6,7 @@ ...@@ -6,6 +6,7 @@
#include <TF1.h> #include <TF1.h>
#include <TFile.h> #include <TFile.h>
#include <TH2.h> #include <TH2.h>
#include <TSpectrum.h>
#include <TSpline.h> #include <TSpline.h>
#include <fstream> #include <fstream>
using namespace std; using namespace std;
...@@ -14,13 +15,16 @@ using namespace std; ...@@ -14,13 +15,16 @@ using namespace std;
const int CutNumberLoader(); const int CutNumberLoader();
const int Ncuts = CutNumberLoader(); const int Ncuts = CutNumberLoader();
const int NSegment=11; const int NSegment=11;
const int Nrun=3; const int Nrun=1;
int RunNumber[Nrun]={246,247,248}; int RunNumber[Nrun]={247};
const char *CutPath = "../Step4/Output/CutAutoZ.root" ;
void MakeSpline(); void MakeSpline();
void ApplySpline_perTCutG(); void ApplySpline_perTCutG();
void ApplySpline(); void ApplySpline();
vector<double> GetMinMaxX(TCutG* cut);
vector<double> GetMinMaxY(TCutG* cut) ;
TF1* Linear(double x1, double y1, double x2, double y2, const char* name = "f");
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
void SplineChio() void SplineChio()
{ {
...@@ -29,7 +33,7 @@ void SplineChio() ...@@ -29,7 +33,7 @@ void SplineChio()
//=========================================================================================================== //===========================================================================================================
//**********************Cut************************************** //**********************Cut**************************************
TFile *fcut=new TFile("../Step4/Output/CutAutoZ.root","open"); TFile *fcut=new TFile(CutPath,"open");
vector<TCutG*> cutZ(Ncuts); vector<TCutG*> cutZ(Ncuts);
for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i)); for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
...@@ -103,7 +107,9 @@ void SplineChio() ...@@ -103,7 +107,9 @@ void SplineChio()
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
void MakeSpline(){ void MakeSpline(){
// *************************Cut***************************
TFile *fcut=new TFile(CutPath, "open");
vector<TCutG*> cutZ(Ncuts);
// ********************** Load prev histo**************************** // ********************** Load prev histo****************************
TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root"); TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root");
...@@ -119,6 +125,7 @@ void MakeSpline(){ ...@@ -119,6 +125,7 @@ void MakeSpline(){
for(int i=0;i<Ncuts;i++) for(int i=0;i<Ncuts;i++)
{ {
// Get the 2D plot in TCutG // Get the 2D plot in TCutG
cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
h2[i]=(TH2F*)inFile->Get(Form("hChioDE_E_%d",i+1)); h2[i]=(TH2F*)inFile->Get(Form("hChioDE_E_%d",i+1));
// Create the TProfile // Create the TProfile
...@@ -152,7 +159,7 @@ void MakeSpline(){ ...@@ -152,7 +159,7 @@ void MakeSpline(){
pfx[i]->SetDirectory(0); pfx[i]->SetDirectory(0);
// find the function to extend the TProfile on the lower range // find the function to extend the TProfile on the lower range
TF1 * fitpol3 = new TF1(Form("FitPol3_pfx_%02d",i+1),"pol3",hProfile[i]->GetBinLowEdge(FirstBin),hProfile[i]->GetBinLowEdge(FirstBin)+10000); TF1 * fitpol3 = new TF1(Form("FitPol3_pfx_%02d",i+1),"pol3",hProfile[i]->GetBinLowEdge(FirstBin),hProfile[i]->GetBinLowEdge(FirstBin)+hProfile[i]->GetBinLowEdge(FirstBin+200));
hProfile[i]->Fit(fitpol3,"R"); hProfile[i]->Fit(fitpol3,"R");
fitpol3->GetParameters(parpol3); fitpol3->GetParameters(parpol3);
...@@ -205,7 +212,13 @@ void MakeSpline(){ ...@@ -205,7 +212,13 @@ void MakeSpline(){
fspline->Close(); fspline->Close();
TGraph * gr = new TGraph(); TFile* CalOut=new TFile("../../../../Calibration/VAMOS/CHIO/Z_spline.root","recreate");
for (TSpline3 * i : gspline){
i->Write();
}
CalOut->Close();
TGraph * gr = new TGraph();
gr->SetName("DEvsZ"); gr->SetName("DEvsZ");
TCanvas * can = new TCanvas("ChioEcorr_vs_Ebis_for_TSplines","ChioEcorr_vs_Ebis_for_TSplines",0,0,2500,2000); TCanvas * can = new TCanvas("ChioEcorr_vs_Ebis_for_TSplines","ChioEcorr_vs_Ebis_for_TSplines",0,0,2500,2000);
...@@ -214,7 +227,22 @@ void MakeSpline(){ ...@@ -214,7 +227,22 @@ void MakeSpline(){
string filename_par = "Output/mean_spline_par.dat"; string filename_par = "Output/mean_spline_par.dat";
ofile_par.open(filename_par.c_str()); ofile_par.open(filename_par.c_str());
cout << "Float_t FF_DEcorr[" << Ncuts << "] = { " << endl; cout << "Float_t FF_DEcorr[" << Ncuts << "] = { " << endl;
ofile_par << "Float_t FF_DEcorr[" << Ncuts << "] = { " << endl; ofile_par << "Float_t FF_DEcorr[" << Ncuts << "] = { " << endl;
vector<double> XEdgesMin = GetMinMaxX(cutZ.at(0));
vector<double> XEdgesMax = GetMinMaxX(cutZ.at(Ncuts-1));
vector<double> YEdgesMin = GetMinMaxY(cutZ.at(0));
vector<double> YEdgesMax = GetMinMaxY(cutZ.at(Ncuts-1));
double x1 = (XEdgesMin.at(0) + XEdgesMin.at(1)) /2;
double x2 = (XEdgesMax.at(0) + XEdgesMax.at(1)) /2;
double y1 = (YEdgesMin.at(0) + YEdgesMin.at(1)) /2;
double y2 = (YEdgesMax.at(0) + YEdgesMax.at(1)) /2;
TF1 *MeanPos = Linear(y1,x1,y2,x2,"MeanPos");
for(int i=0;i<Ncuts;i++){ for(int i=0;i<Ncuts;i++){
can->cd(); can->cd();
if(i==0) h2[i]->Draw("col"); if(i==0) h2[i]->Draw("col");
...@@ -222,7 +250,7 @@ void MakeSpline(){ ...@@ -222,7 +250,7 @@ void MakeSpline(){
cout << h2[i]->GetMean(2) << ", " << endl; cout << h2[i]->GetMean(2) << ", " << endl;
ofile_par << h2[i]->GetMean(2) << ", " << endl; ofile_par << h2[i]->GetMean(2) << ", " << endl;
//gr->SetPoint(i,sqrt(h2[i]->GetMean(2)),i+31); //gr->SetPoint(i,sqrt(h2[i]->GetMean(2)),i+31);
gr->SetPoint(i,gspline[i]->Eval(8500),i+Z59POS); gr->SetPoint(i,gspline[i]->Eval(MeanPos->Eval((YEdgesMin.at(0)+YEdgesMin.at(1))/2)),i+Z59POS);
//hProfile[i]->Draw("same"); //hProfile[i]->Draw("same");
pfx[i]->Draw("same"); pfx[i]->Draw("same");
gspline[i]->SetLineColor(kBlue); gspline[i]->SetLineColor(kBlue);
...@@ -243,7 +271,7 @@ void ApplySpline_perTCutG() ...@@ -243,7 +271,7 @@ void ApplySpline_perTCutG()
{ {
// *************************Cut*************************** // *************************Cut***************************
TFile *fcut=new TFile("../Step4/Output/CutAutoZ.root"); TFile *fcut=new TFile(CutPath);
vector<TCutG*> cutZ(Ncuts); vector<TCutG*> cutZ(Ncuts);
// ******************* TSpline DE************************* // ******************* TSpline DE*************************
...@@ -346,12 +374,10 @@ can->cd(2); gPad-> SetLogz(); hChioDE_E->Draw("colz"); ...@@ -346,12 +374,10 @@ can->cd(2); gPad-> SetLogz(); hChioDE_E->Draw("colz");
can->cd(3); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz"); can->cd(3); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz");
} }
void ApplySplineUnique()
///////////////////////////////////////////////////////////////////////////
void ApplySpline()
{ {
// *************************Cut*************************** // *************************Cut***************************
TFile *fcut=new TFile("../Step4/Output/CutAutoZ.root", "open"); TFile *fcut=new TFile(CutPath, "open");
vector<TCutG*> cutZ(Ncuts); vector<TCutG*> cutZ(Ncuts);
// ******************* TSpline DE************************* // ******************* TSpline DE*************************
...@@ -385,41 +411,6 @@ void ApplySpline() ...@@ -385,41 +411,6 @@ void ApplySpline()
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);
} }
/*double FF_DEcorr0[Ncuts] = {
7375.73,
7982.83,
8419.21,
8801.52,
9227.89,
9667.77,
10163.1,
10680.5,
11222.2,
11758.3,
12287.6,
12802.2,
13324.2,
13828.4,
14325.2,
14773.7,
15207.3,
15674.6,
16098.9,
16548.5,
16926.9,
17384.1,
17778.3,
18196.6,
18597.6,
19043.9,
19437.8,
19889.4,
20301.1,
20709.7,
21112.6,
21544,
21953.4,
22361.8};*/
// histos // histos
TH2F *hChioDE_E_all = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000); TH2F *hChioDE_E_all = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000);
...@@ -433,6 +424,158 @@ void ApplySpline() ...@@ -433,6 +424,158 @@ void ApplySpline()
int Nentries = chain->GetEntries(); 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;
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;
}
FF_DEcorr = -100;
chain->GetEntry(e);
if(IC->Eres==FF_Eres_prev) continue;
FF_Eres_prev = IC->Eres;
hChioDE_E_all->Fill(IC->Eres,IC->DE);
float Eval_DEspline, DEspline0;
int index=0;
Eval_DEspline = gspline[24]->Eval(IC->Eres);
DEspline0 = FF_DEcorr0[24];
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 << " " << IC->Eres << " " << FF_DEcorr << endl;
//}
// should be pol2 !!!
//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(IC->Eres,Zrough);
if(3000<IC->Eres && IC->Eres<25000){
hChioDE_corr->Fill(FF_DEcorr);
hChioZ_rough->Fill(Zrough);
hChioZ_abitbetter->Fill(Zabitbetter);
}
}
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,2000);
can->Divide(1,3);
can->cd(1); gPad-> SetLogz(); hChioDE_E_all->Draw("colz");
for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
can->cd(2); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz");
can->cd(3); gPad-> SetLogz(); hChioZ_E_rough->Draw("colz");
TCanvas * can1D = new TCanvas(Form("DE_and_Z_ifE_8000to25000_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
Form("DE_and_Z_ifE_8000to25000_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
0, 0, 2000,2000);
can1D->Divide(1,3);
can1D->cd(1); hChioDE_corr->Draw();
can1D->cd(2); hChioZ_rough->Draw();
can1D->cd(3); hChioZ_abitbetter->Draw();
TFile * fsave = new TFile(Form("Output/DEvsE_corr_run%04d_%04d.root",RunNumber[0],RunNumber[Nrun-1]),"recreate");
can->Write();
hChioDE_E_all->Write();
hChioDE_E_corr->Write();
hChioZ_E_rough->Write();
can1D->Write();
hChioDE_corr->Write();
hChioZ_rough->Write();
hChioZ_abitbetter->Write();
fsave->Close();
}
///////////////////////////////////////////////////////////////////////////
void ApplySpline()
{
// *************************Cut***************************
TFile *fcut=new TFile(CutPath, "open");
vector<TCutG*> cutZ(Ncuts);
// ******************* TSpline DE*************************
TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
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));
}
//**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree");
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() ;
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
// variable after applying TSpline
double FF_DEcorr=0;
double Zrough=0;
double Zabitbetter=0;
vector<double> FF_DEcorr0(Ncuts);
vector<double> XEdgesMin = GetMinMaxX(cutZ.at(0));
vector<double> XEdgesMax = GetMinMaxX(cutZ.at(Ncuts-1));
vector<double> YEdgesMin = GetMinMaxY(cutZ.at(0));
vector<double> YEdgesMax = GetMinMaxY(cutZ.at(Ncuts-1));
double x1 = (XEdgesMin.at(0) + XEdgesMin.at(1)) /2;
double x2 = (XEdgesMax.at(0) + XEdgesMax.at(1)) /2;
double y1 = (YEdgesMin.at(0) + YEdgesMin.at(1)) /2;
double y2 = (YEdgesMax.at(0) + YEdgesMax.at(1)) /2;
TF1 *MeanPos = Linear(y1,x1,y2,x2,"MeanPos");
for(int index=0; index<Ncuts; index++){
vector<double> YEdgesMinTemp = GetMinMaxY(cutZ.at(index));
cout << " E " <<MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2) << endl;
cout << "DE " << gspline[index]->Eval(MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2)) << endl;
FF_DEcorr0[index] = gspline[index]->Eval(MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2));
}
// histos
TH2F *hChioDE_E_all = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000);
TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);
TH1F *hChioDE_corr = new TH1F("hChioDE_corr","hChioDE_corr",1000,1000,26000);
TH2F *hChioZ_E_rough = new TH2F("hChioZ_E_rough","hChioZ_E_rough",1000,1000,41000,500,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);
std::ifstream file("Output/Calibration.txt");
bool Loaded = true ;
if (!file) {
std::cerr << "Error opening file.\n";
Loaded = false;
}
Double_t p[6]; // Assuming 5 values
if (Loaded == true){
for (int i = 0; i < 6; ++i) {
file >> p[i];
}
}
int Nentries=1e6;
//int Nentries = chain->GetEntries();
auto start = std::chrono::high_resolution_clock::now();
double FF_Eres_prev = 0; double FF_Eres_prev = 0;
for(int e=0;e<Nentries;e++) for(int e=0;e<Nentries;e++)
{ {
...@@ -466,8 +609,8 @@ void ApplySpline() ...@@ -466,8 +609,8 @@ void ApplySpline()
DEspline0 = FF_DEcorr0[index]; DEspline0 = FF_DEcorr0[index];
float dmin, dsup; float dmin, dsup;
if( index < (Ncuts-1) && IC->DE > gspline[0]->Eval(IC->Eres)){ if( index < (Ncuts-1) && IC->DE > gspline[0]->Eval(IC->Eres)){
dmin = IC->DE-gspline[index]->Eval(IC->Eres); dmin = abs(IC->DE - gspline[index]->Eval(IC->Eres));
dsup = gspline[index+1]->Eval(IC->Eres)-IC->DE; dsup = abs(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) {
...@@ -486,7 +629,17 @@ void ApplySpline() ...@@ -486,7 +629,17 @@ void ApplySpline()
// should be pol2 !!! // should be pol2 !!!
//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);
if (Loaded == true){
Zrough = p[0] + p[1]*FF_DEcorr +p[2]*pow(FF_DEcorr,2) + p[3]*pow(FF_DEcorr,3) +p[4]*pow(FF_DEcorr,4) + p[5]*pow(FF_DEcorr,5);
}
else {
Zrough = -57.068 + 0.0275327*FF_DEcorr - 3.52303e-06*pow(FF_DEcorr,2) + 2.36906e-10*pow(FF_DEcorr,3) - 7.74553e-15*pow(FF_DEcorr,4) + 9.90229e-20*pow(FF_DEcorr,5);
}
//out << Zrough << endl;
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(IC->Eres,Zrough); hChioZ_E_rough->Fill(IC->Eres,Zrough);
...@@ -497,6 +650,58 @@ void ApplySpline() ...@@ -497,6 +650,58 @@ void ApplySpline()
} }
} }
// ********************** Fit pol5 of each Z ***********************************
//
// PeakFinder
TSpectrum *PeakFinder = new TSpectrum(100);
int nPeaks = PeakFinder->Search(hChioDE_corr,3,"",0.043);
Double_t * x = PeakFinder->GetPositionX();
std::sort(x,x+nPeaks);
TGraph * gr = new TGraph();
for(int i=0;i<nPeaks;i++){
cout << x[i] << endl;
gr->SetPoint(i,x[i],i+32);
}
// Define and fit a function (linear in this case)
TF1 *fitFunc = new TF1("fitFunc", "pol5");
gr->Fit(fitFunc, "Q"); // Quiet mode
// Open a text file to save fit results
std::ofstream outFile("Output/Calibration.txt");
std::ofstream outFileCal("../../../../Calibration/VAMOS/CHIO/Chio_Z_Calibration.cal");
outFileCal << "IC_Z_CALIBRATION" << " ";
if (outFile.is_open()) {
for (int i = 0; i < fitFunc->GetNpar(); i++) {
outFile << fitFunc->GetParameter(i) << " ";
outFileCal << fitFunc->GetParameter(i) << " ";
}
outFile.close();
outFileCal.close();
std::cout << "Fit results saved to fit_results.txt\n";
}
std::ofstream CalFile("../../../../Calibration/VAMOS/CHIO/Chio_Z_Spline_Eval.txt");
if (CalFile.is_open()) {
for (int i = 0; i < cutZ.size(); i++) {
vector<double> YEdgesMinTemp = GetMinMaxY(cutZ.at(i));
CalFile << MeanPos->Eval((YEdgesMinTemp.at(0)+YEdgesMinTemp.at(1))/2) << " ";
}
CalFile.close();
std::cout << "Eval position saved in calibration\n";
}
TCanvas * cangr = new TCanvas("Cal","Cal",200,0,1200,1000);
cangr->cd();
gr->SetMarkerStyle(20);
gr->Draw("AP");
gr->Fit("pol5");
// ******************************************************************************
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,2000); 0,0,2000,2000);
...@@ -590,8 +795,9 @@ void Fit_DE_E_corr() ...@@ -590,8 +795,9 @@ void Fit_DE_E_corr()
} }
const int CutNumberLoader(){ const int CutNumberLoader(){
TFile *fcut=new TFile("../Step4/Output/CutAutoZ.root","open"); TFile *fcut=new TFile(CutPath,"open");
int CutCount = 0 ; int CutCount = 0 ;
TIter next(fcut->GetListOfKeys()); TIter next(fcut->GetListOfKeys());
...@@ -608,3 +814,55 @@ const int CutNumberLoader(){ ...@@ -608,3 +814,55 @@ const int CutNumberLoader(){
} }
std::vector<double> GetMinMaxX(TCutG* cut) {
if (!cut) {
std::cerr << "Invalid TCutG!" << std::endl;
return {};
}
int nPoints = cut->GetN();
double* xPoints = cut->GetX();
// Find the minimum and maximum X values
double minX = *std::min_element(xPoints, xPoints + nPoints);
double maxX = *std::max_element(xPoints, xPoints + nPoints);
return {minX, maxX};
}
std::vector<double> GetMinMaxY(TCutG* cut) {
if (!cut) {
std::cerr << "Invalid TCutG!" << std::endl;
return {};
}
int nPoints = cut->GetN();
double* yPoints = cut->GetY();
// Find the minimum and maximum X values
double minX = *std::min_element(yPoints, yPoints + nPoints);
double maxX = *std::max_element(yPoints, yPoints + nPoints);
return {minX, maxX};
}
TF1* Linear(double x1, double y1, double x2, double y2, const char* name ) {
// Ensure the points are not the same
if (x1 == x2) {
throw std::invalid_argument("x1 and x2 must be different to define a function.");
}
// Calculate slope and intercept of the line
double slope = (y2 - y1) / (x2 - x1);
double intercept = y1 - slope * x1;
// Define the function as a lambda
auto func = [slope, intercept](double* x, double* p) {
return slope * x[0] + intercept;
};
// Create and return the TF1
TF1* f = new TF1(name, func, x1, x2, 0);
return f;
}
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