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

[AlPhaPha] Optimized some macro

parent e7e0a4fa
No related branches found
No related tags found
No related merge requests found
Pipeline #390370 canceled
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/SplineChio.C
+ 10
− 8
View file @ 3ef27d99
......@@ -122,7 +122,7 @@ void MakeSpline(){
int Z59POS = 28;
vector<TSpline3*> gspline(Ncuts);
vector<TH2F*> h2 (Ncuts);
vector<TH1F*> hProfile(Ncuts);
vector<TProfile*> hProfile(Ncuts);
vector<TH1F*> pfx(Ncuts); // extended hProfile
TFile* fspline=new TFile("Output/spline_Chio_2024.root","recreate");
......@@ -135,7 +135,7 @@ void MakeSpline(){
h2[i]=(TH2F*)inFile->Get(Form("hChioDE_E_%d",i+1));
// Create the TProfile
hProfile[i]=(TH1F*)h2[i]->ProfileX();
hProfile[i]=h2[i]->ProfileX();
hProfile[i]->SetLineWidth(2);
hProfile[i]->SetDirectory(0);
canfit->cd();
......@@ -150,12 +150,12 @@ void MakeSpline(){
double parpol3[4];
for(int bin=1; bin<hProfile[i]->GetNbinsX(); bin++){
FirstBin = bin;
if (hProfile[i]->GetBinContent(bin)>6000) break;
if (hProfile[i]->GetBinEntries(bin)>10) break;
}
double parpol1[2];
for(int bin=hProfile[i]->GetNbinsX(); bin>1 ; bin--){
LastBin = bin;
if (hProfile[i]->GetBinContent(bin)>6000) break;
if (hProfile[i]->GetBinEntries(bin)>10) break;
}
// Create the extended TProfile
......@@ -582,8 +582,8 @@ void ApplySpline()
file >> p[i];
}
}
int Nentries=2e6;
//int Nentries = chain->GetEntries();
//int Nentries=2e6;
int Nentries = chain->GetEntries();
auto start = std::chrono::high_resolution_clock::now();
double FF_Eres_prev = 0;
......@@ -661,6 +661,7 @@ void ApplySpline()
}
}
/*
// ********************** Fit pol5 of each Z ***********************************
//
// PeakFinder
......@@ -693,7 +694,7 @@ void ApplySpline()
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++) {
......@@ -706,12 +707,13 @@ void ApplySpline()
/*
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]),
Form("ChioEbis_vs_ChiodE_run%04d_%04d",RunNumber[0],RunNumber[Nrun-1]),
......
Projects/AlPhaPha/2024/macro/Calibration/IC/Step6/CorrectionZ.C
+ 352
− 17
View file @ 3ef27d99
#include "TICPhysics.h"
#include "TFPMWPhysics.h"
#include <TChain.h>
#include <TCutG.h>
#include <TFileMerger.h>
#include <regex>
void HistoFiller(const char* Path);
void MakeSpline(TH2F *hIn , const char *outputPath , double yMin = 5000 , double yMax =26000, double iter=0 );
void ApplySpline(const char* Path);
void CopyHistogramWithTCutG(TH2F* h, TCutG* cut);
void CorrectionZ(const char* Path){
void CorrectionZ(const char* Path, bool Spline = 0, bool Cut = 0){
//===========================================================================================================
if( Path && Path[0] != '\0'){
HistoFiller(Path);
return;
}
// Fill corrected or uncorrected histo
//===========================================================================================================
if( Path && Path[0] != '\0'){
if (Spline == 0){
HistoFiller(Path);
return;
}
TFile *f = new TFile("Output/histo_merged.root","open");
else if (Spline == 1){
ApplySpline(Path);
return;
}
}
//===========================================================================================================
// Load filled histo
//===========================================================================================================
TH2F *hZvsTS = nullptr , *hZvsTheta = nullptr;
TFile *f ;
if (Spline == 0){
f = new TFile("Output/histo_merged.root","open");
}
else if (Spline == 1){
f = new TFile("Output/histo_corrected.root","open");
}
TH2F *hZvsTS = nullptr ;
TH2F *hZvsThetaf = nullptr, *hZvsThetafLin= nullptr , *hZvsThetaIn = nullptr, *hZvsThetaInLin=nullptr ;
f->GetObject("hZvsTS",hZvsTS);
f->GetObject("hZvsTheta",hZvsTheta);
f->GetObject("hZvsThetaf",hZvsThetaf);
f->GetObject("hZvsThetafLin",hZvsThetafLin);
f->GetObject("hZvsThetaIn",hZvsThetaIn);
f->GetObject("hZvsThetaInLin",hZvsThetaInLin);
//===========================================================================================================
if (Spline == 0){
double zMin = 5000 , zMax =26000 ;
//===========================================================================================================
// If there is a cut load it
//===========================================================================================================
TFile *fcut=new TFile("Output/CutZ.root", "open");
if (fcut && fcut->IsOpen()){
TCutG *cutZThetaf = (TCutG*)fcut->Get("Cutf");
TCutG *cutZThetaIn = (TCutG*)fcut->Get("CutIn");
if (cutZThetaf && Cut) CopyHistogramWithTCutG(hZvsThetafLin,cutZThetaf);
if (cutZThetaf && Cut) CopyHistogramWithTCutG(hZvsThetaf,cutZThetaf);
if (cutZThetaIn&& Cut) CopyHistogramWithTCutG(hZvsThetaIn,cutZThetaIn);
if (cutZThetaIn&& Cut) CopyHistogramWithTCutG(hZvsThetaInLin,cutZThetaIn);
}
//===========================================================================================================
MakeSpline(hZvsThetaf , "Output/SplineThetaf.root" , zMin ,zMax, 0);
MakeSpline(hZvsThetafLin , "Output/SplineThetafLin.root" , zMin ,zMax, 1);
MakeSpline(hZvsThetaIn , "Output/SplineThetaIn.root" , zMin ,zMax, 2);
MakeSpline(hZvsThetaInLin , "Output/SplineThetaInLin.root" , zMin ,zMax,3);
// Merge spline file
TFileMerger merger;
merger.SetFastMethod(kTRUE);
// Add file to be merged
merger.AddFile("Output/SplineThetaf.root" );
merger.AddFile("Output/SplineThetafLin.root" );
merger.AddFile("Output/SplineThetaIn.root" );
merger.AddFile("Output/SplineThetaInLin.root");
merger.OutputFile("Output/SplineTheta.root");
merger.Merge();
}
else if (Spline == 1){
// Save the spline to the Calibration
TFile *fSpline = new TFile("Output/SplineThetafLin.root","open");
TFile *fCal = new TFile("../../../../Calibration/VAMOS/CHIO/Spline_Z_Thetaf.root","recreate");
TSpline3* sThetafLin = (TSpline3*)fSpline->FindObjectAny(Form("fspline_%d",1));
sThetafLin->SetName("fspline_1");
sThetafLin->Write();
fSpline->Close();
fCal->Close();
}
//===========================================================================================================
TCanvas *c = new TCanvas("c","c");
c->Divide(2);
c->Divide(2,2);
c->cd(1);
hZvsTheta->Draw("colz");
hZvsThetaf->Draw("colz");
c->cd(2);
hZvsTS->Draw("colz");
hZvsThetafLin->Draw("colz");
c->cd(3);
hZvsThetaIn->Draw("colz");
c->cd(4);
hZvsThetaInLin->Draw("colz");
double y_min = 44.5; // Set your y-range minimum
double y_max = 45.5; // Set your y-range maximum
TProfile* profX = new TProfile("profX", "Profile in y-range", hZvsThetafLin->GetNbinsX(), hZvsThetafLin->GetXaxis()->GetXmin(), hZvsThetafLin->GetXaxis()->GetXmax());
for (int i = 1; i <= hZvsThetafLin->GetNbinsX(); ++i) {
for (int j = 1; j <= hZvsThetafLin->GetNbinsY(); ++j) {
double y_value = hZvsThetafLin->GetYaxis()->GetBinCenter(j);
if (y_value >= y_min && y_value <= y_max) {
int bin = hZvsThetafLin->GetBin(i, j);
double weight = hZvsThetafLin->GetBinContent(bin);
profX->Fill(hZvsThetafLin->GetXaxis()->GetBinCenter(i), y_value, weight);
}
}
}
TCanvas *c2 = new TCanvas("c2","c2");
profX->GetYaxis()->SetRangeUser(44.5,45.5);
profX->Draw();
} // end macro:
void HistoFiller(const char* Path){
......@@ -56,9 +151,16 @@ void HistoFiller(const char* Path){
//===========================================================================================================
//===========================================================================================================
// Histo initialisation + filling just to get the range
TH2F *hZvsTheta = new TH2F("hZvsTheta","hZvsTheta",1000,-0.5,0.5,1000,20,80);
TH2F *hZvsTS = new TH2F("hZvsTS","hZvsTS",2000,0,200,1000,20,80);
// Histo initialisation + filling just to get the range
//double zMin = 44.5 , zMax =45.5 ;
double zMin = 5000 , zMax = 26000 ;
double thetaMin = -0.3 , thetaMax= 0.3;
int Bin = 1000;
TH2F *hZvsThetaf = new TH2F("hZvsThetaf","hZvsThetaf",Bin,thetaMin,thetaMax,Bin,zMin,zMax);
TH2F *hZvsThetaIn = new TH2F("hZvsThetaIn","hZvsThetaIn",Bin,thetaMin,thetaMax , Bin,zMin,zMax);
TH2F *hZvsThetafLin = new TH2F("hZvsThetafLin","hZvsThetafLin",Bin,thetaMin,thetaMax,Bin,zMin,zMax);
TH2F *hZvsThetaInLin = new TH2F("hZvsThetaInLin","hZvsThetaInLin",Bin,thetaMin,thetaMax,Bin,zMin,zMax);
TH2F *hZvsTS = new TH2F("hZvsTS","hZvsTS",10000,0,1000,1000,zMin,zMax);
//===========================================================================================================
int Nentries = chain->GetEntries();
//int Nentries = 1e6;
......@@ -77,8 +179,11 @@ void HistoFiller(const char* Path){
chain->GetEntry(e);
if(FF_IC_X>-530 && VAMOS_TS_hour>0){
hZvsTheta->Fill(FPMW->Thetaf,IC->Chio_Z);
hZvsTS->Fill(VAMOS_TS_hour,IC->Chio_Z);
hZvsThetaf->Fill(FPMW->Thetaf,IC->Chio_ZRaw);
hZvsThetafLin->Fill(FPMW->ThetafLin,IC->Chio_ZRaw);
hZvsThetaIn->Fill(FPMW->Theta_in,IC->Chio_ZRaw);
hZvsThetaInLin->Fill(FPMW->Theta_in_Lin,IC->Chio_ZRaw);
hZvsTS->Fill(VAMOS_TS_hour,IC->Chio_ZRaw);
}
}
......@@ -103,7 +208,10 @@ void HistoFiller(const char* Path){
}
hZvsTS->Write();
hZvsTheta->Write();
hZvsThetaf->Write();
hZvsThetafLin->Write();
hZvsThetaInLin->Write();
hZvsThetaIn->Write();
// Don't forget to close the file
file->Close();
}
......@@ -112,3 +220,230 @@ void HistoFiller(const char* Path){
}
}
void MakeSpline(TH2F *hIn , const char *outputPath , double yMin , double yMax , double iter ){
//===========================================================================================================
TSpline3 *gspline;
TProfile *hProfile;
TH1F *pfx; // extended hProfile
TFile *fspline=new TFile(outputPath,"recreate");
//===========================================================================================================
TCanvas * canfit = new TCanvas(Form("canfit_%f",iter),Form("canfit_%f",iter),0,0,2000,1500);
// Create the TProfile
hProfile=hIn->ProfileX();
hProfile->SetLineWidth(2);
hProfile->SetDirectory(0);
canfit->cd();
hProfile->GetYaxis()->SetRangeUser(yMin,yMax);
hProfile->Draw();
//===========================================================================================================
// Get the range of the TProfile
int FirstBin, LastBin;
double parpol3[4];
for(int bin=1; bin<hProfile->GetNbinsX(); bin++){
FirstBin = bin;
if ((hProfile->GetBinEntries(bin)>300) &&(hProfile->GetBinEntries(bin+2)>200)) break;
}
double parpol2[2];
for(int bin=hProfile->GetNbinsX(); bin>1 ; bin--){
LastBin = bin;
if (hProfile->GetBinEntries(bin)>300) break;
}
cout << "FB " << FirstBin << " LB " << LastBin << endl;
//===========================================================================================================
// Create the extended TProfile
pfx = new TH1F("pfx","pfx",hProfile->GetNbinsX(),hProfile->GetBinLowEdge(1),
hProfile->GetBinLowEdge(hProfile->GetNbinsX())+hProfile->GetBinWidth(hProfile->GetNbinsX()));
pfx->SetLineColor(8);
pfx->SetDirectory(0);
// find the function to extend the TProfile on the lower range
TF1 * fitpol3 = new TF1(Form("FitPol3_pfx_%f",iter),"pol3",hProfile->GetBinLowEdge(FirstBin),hProfile->GetBinLowEdge(FirstBin+300));
hProfile->Fit(fitpol3,"R");
fitpol3->GetParameters(parpol3);
// find the function to extend the TProfile on the higher range
TF1 * fitpol2 = new TF1(Form("FitPol1_pfx_%f",iter),"pol2",hProfile->GetBinLowEdge(LastBin-50),hProfile->GetBinLowEdge(LastBin));
hProfile->Fit(fitpol2,"R");
fitpol2->GetParameters(parpol2);
//===========================================================================================================
// fill the extended TProfile
float newval,lastval,lasterr;
for(int bin=1; bin<=FirstBin; bin++){
newval=0;
for(int par=0; par<4; par++) newval += parpol3[par]*pow(hProfile->GetBinCenter(bin),par);
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);
}
}
for(int bin=LastBin; bin<=hProfile->GetNbinsX(); bin++){
newval=0;
for(int par=0; par<3; par++) newval += parpol2[par]*pow(hProfile->GetBinCenter(bin),par);
pfx->SetBinContent(bin,newval);
}
pfx->Draw("same");
gspline=new TSpline3(pfx);
int i = iter;
gspline->SetName(Form("fspline_%d",i));
fspline->cd();
gspline->Write();
fspline->Close();
}
void ApplySpline(const char* Path){
//===========================================================================================================
// Load Spline
//===========================================================================================================
// Load the merged file
TFile *fIn = new TFile("Output/SplineTheta.root","open");
TSpline3 *sThetaf , *sThetafLin , *sThetaIn , *sThetaInLin;
sThetaf = (TSpline3*)fIn->FindObjectAny(Form("fspline_%d",0));
sThetafLin = (TSpline3*)fIn->FindObjectAny(Form("fspline_%d",1));
sThetaIn = (TSpline3*)fIn->FindObjectAny(Form("fspline_%d",2));
sThetaInLin = (TSpline3*)fIn->FindObjectAny(Form("fspline_%d",3));
//===========================================================================================================
// Fill
//===========================================================================================================
//===========================================================================================================
// Import
//===========================================================================================================
TChain *chain = new TChain("PhysicsTree");
chain->Add (Path);
//===========================================================================================================
TICPhysics* IC = new TICPhysics() ;
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
TFPMWPhysics* FPMW = new TFPMWPhysics() ;
chain->SetBranchStatus("FPMW", true);
chain->SetBranchAddress("FPMW", &FPMW);
double FF_IC_X, VAMOS_TS_hour;
chain->SetBranchStatus("FF_IC_X", true);
chain->SetBranchAddress("FF_IC_X", &FF_IC_X);
chain->SetBranchStatus("VAMOS_TS_hour", true);
chain->SetBranchAddress("VAMOS_TS_hour", &VAMOS_TS_hour);
//===========================================================================================================
//===========================================================================================================
// Histo initialisation + filling just to get the range
double zMin = 5000 , zMax = 26000 ;
double thetaMin = -0.15 , thetaMax= 0.15;
int Bin = 1000;
TH2F *hZvsThetaf = new TH2F("hZvsThetaf","hZvsThetaf",Bin,thetaMin,thetaMax,Bin,zMin,zMax);
TH2F *hZvsThetaIn = new TH2F("hZvsThetaIn","hZvsThetaIn",Bin,thetaMin,thetaMax , Bin,zMin,zMax);
TH2F *hZvsThetafLin = new TH2F("hZvsThetafLin","hZvsThetafLin",Bin,thetaMin,thetaMax,Bin,zMin,zMax);
TH2F *hZvsThetaInLin = new TH2F("hZvsThetaInLin","hZvsThetaInLin",Bin,thetaMin,thetaMax,Bin,zMin,zMax);
TH2F *hZvsTS = new TH2F("hZvsTS","hZvsTS",10000,0,1000,1000,zMin,zMax);
//===========================================================================================================
int Nentries = chain->GetEntries();
//int Nentries = 1e6;
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);
if(FF_IC_X>-530 && VAMOS_TS_hour>0){
//cout << " Avant correction " << IC->Chio_ZRaw << endl;
//double newZ = sThetafLin->Eval(0)/sThetafLin->Eval(FPMW->ThetafLin)*IC->Chio_ZRaw;
double newZ = sThetafLin->Eval(0)/sThetafLin->Eval(FPMW->ThetafLin)*IC->Chio_ZRaw;
//cout << " Thetaf " << FPMW->Thetaf << " newZ " << newZ << endl;
hZvsThetaf->Fill(FPMW->Thetaf ,sThetaf->Eval(0)/sThetaf->Eval(FPMW->Thetaf)*IC->Chio_ZRaw);
hZvsThetafLin->Fill(FPMW->ThetafLin ,sThetafLin->Eval(0)/sThetafLin->Eval(FPMW->ThetafLin)*IC->Chio_ZRaw);
hZvsThetaIn->Fill(FPMW->Theta_in ,sThetaIn->Eval(0)/sThetaIn->Eval(FPMW->Theta_in)*IC->Chio_ZRaw);
hZvsThetaInLin->Fill(FPMW->Theta_in_Lin,sThetaInLin->Eval(0)/sThetaInLin->Eval(FPMW->Theta_in_Lin)*IC->Chio_ZRaw);
hZvsTS->Fill(VAMOS_TS_hour,IC->Chio_ZRaw);
}
}
//===========================================================================================================
// Save file
//===========================================================================================================
string strPath(Path);
regex pattern("Run(\\d{3})"); // Regex to match Run followed by exactly 3 digits
smatch matches;
if (regex_search(strPath, matches, pattern)) {
string number = matches[1]; // Extract the number
string fileName = "Output/histo_Corrected_" + number + ".root"; // Generate the filename
// Create the TFile
TFile* file = new TFile(fileName.c_str(), "RECREATE");
if (file->IsOpen()) {
cout << "File created: " << fileName << endl;
} else {
cout << "Failed to create the file." << endl;
}
hZvsTS->Write();
hZvsThetaf->Write();
hZvsThetafLin->Write();
hZvsThetaInLin->Write();
hZvsThetaIn->Write();
// Don't forget to close the file
file->Close();
}
else {
cout << "Pattern not found!" << endl;
}
}
void CopyHistogramWithTCutG(TH2F* h, TCutG* cut) {
TH2F* hnew = (TH2F*)h->Clone("h_cut");
hnew->Reset(); // Clear the content while keeping the binning
for (int i = 1; i <= h->GetNbinsX(); ++i) {
for (int j = 1; j <= h->GetNbinsY(); ++j) {
double x = h->GetXaxis()->GetBinCenter(i);
double y = h->GetYaxis()->GetBinCenter(j);
if (cut->IsInside(x, y)) { // Check if inside the graphical cut
hnew->SetBinContent(i, j, h->GetBinContent(i, j));
hnew->SetBinError(i, j, h->GetBinError(i, j));
}
}
}
h->Reset();
*h = *hnew;
hnew->Delete();
}
Projects/AlPhaPha/2024/macro/Calibration/IC/Step6/MultiRun.sh deleted 100644 → 0
+ 0
− 54
View file @ e7e0a4fa
#!/bin/bash
#launch job in parralel
MAX_JOBS=$(( $(nproc) -4 )) # Maximum parallel jobs
SEMAPHORE="job_semaphore"
# Create a named pipe (FIFO) for job control
mkfifo $SEMAPHORE
exec 3<> $SEMAPHORE
rm $SEMAPHORE
# Pre-fill the semaphore with tokens equal to MAX_JOBS
for ((i=0; i<MAX_JOBS; i++)); do
echo "" >&3
done
for run in {245..252}; do
read -u 3 # Wait for an available slot
{
root -l -q "./CorrectionZ.C(\"../../../../root/analysis/Run"$run"AoQ_test.root\")"
echo "" >&3 # Release the slot
} &
done
wait
exec 3>&- #close semaphore
# Define output file name
output_file="Output/histo_merged.root"
# Initialize an empty list of files
files=""
# Loop through the run numbers
for run in {245..252}; do
filename="Output/histo_${run}.root"
# Check if the file exists before adding it
if [[ -f "$filename" ]]; then
files+=" $filename"
else
echo "Warning: $filename not found, skipping."
fi
done
# Perform hadd only if there are files
if [[ -n "$files" ]]; then
echo "Merging files into $output_file..."
hadd "$output_file" $files
else
echo "No valid files found. Nothing to merge."
fi
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