diff --git a/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Chio_Z_Calibration.cal b/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Chio_Z_Calibration.cal
index d1e6d52bcf6ef4c414ffb70580310fae394e9920..4fa19c2c756121469e3682499d144f52b5f7958b 100644
--- a/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Chio_Z_Calibration.cal
+++ b/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Chio_Z_Calibration.cal
@@ -1 +1 @@
-IC_Z_CALIBRATION -8.06274 0.00930079 -8.74778e-07 4.91672e-11 -1.24593e-15 1.10243e-20
\ No newline at end of file
+IC_Z_CALIBRATION -32.1669 0.0175119 -2.00422e-06 1.24893e-10 -3.72936e-15 4.29223e-20
\ No newline at end of file
diff --git a/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Z_spline.root b/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Z_spline.root
index de8f9034b29ee7d36cb2289c30245aeff3c173f0..dc383c18a3778acc00c0de64c7f8aea895c0c341 100644
Binary files a/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Z_spline.root and b/Projects/AlPhaPha/2024/Calibration/VAMOS/CHIO/Z_spline.root differ
diff --git a/Projects/AlPhaPha/2024/Snakefile b/Projects/AlPhaPha/2024/Snakefile
index c2b42c1666d746178afc26f844ea7121daa18566..c16e34a1fc805673b3359de151ffab358385a6de 100644
--- a/Projects/AlPhaPha/2024/Snakefile
+++ b/Projects/AlPhaPha/2024/Snakefile
@@ -3,7 +3,7 @@ import subprocess
# Lire le répertoire d'entrée depuis les arguments de configuration
#input_directory = config["folder"]
-input_directory = os.getcwd() + "/../DataMacro/output/run_247"
+input_directory = os.getcwd() + "/../DataMacro/output/run_252"
origin = []
# Iterate over files in input_directory
for filename in os.listdir(input_directory):
@@ -12,14 +12,14 @@ for filename in os.listdir(input_directory):
origin.append(filename)
# Définir le répertoire de sortie pour les fichiers convertis
-phy_directory = os.getcwd() + "/../DataMacro/output/analysis/run_247"
+phy_directory = os.getcwd() + "/../DataMacro/output/analysis/run_252"
#phy_directory = "./"
# define target files directory
analysedfile = []
for inputfile in origin:
#analysedfile.append("/home/morfouacep/Physics/NPTool/nptool/Projects/ana_e850/root/analysis/"+inputfile.replace("_raw_","_"))
- analysedfile.append( os.getcwd() + "/../DataMacro/output/analysis/run_247/"+inputfile)
+ analysedfile.append( os.getcwd() + "/../DataMacro/output/analysis/run_252/"+inputfile)
## batch rules
rule all:
diff --git a/Projects/AlPhaPha/2024/convert_snakemake_generic.sh b/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
index 352444b7dd95ed9aeffc48e11c29e62ecc74d455..0d258992645730e45df4f8cfae6503c061b75e5c 100755
--- a/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
+++ b/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
@@ -1,7 +1,8 @@
#!/bin/bash
#Declare run number
-RUN_NUMBER="247"
+#RUN_NUMBER="246"
+RUN_NUMBER="$1"
# Check if Snakefile exists
if [[ ! -f "Snakefile" ]]; then
@@ -19,7 +20,7 @@ snakemake --cores 30 --forceall --keep-incomplete --keep-going --rerun-incomplet
echo "- snakemake executed successfully!"
echo "- Merging file..."
-OName="\"Run${RUN_NUMBER}ZSpline\""
+OName="\"Run${RUN_NUMBER}\""
OPATH="\"root/analysis\""
Path="\"../DataMacro/output/analysis/run_${RUN_NUMBER}/run_raw_${RUN_NUMBER}_\""
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/ReadMe.md b/Projects/AlPhaPha/2024/macro/Calibration/IC/ReadMe.md
index 12fd80c691be4e977900573264d31a72d0d12e24..213558496a04b07803ce30b96754cf597a664af0 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/ReadMe.md
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/ReadMe.md
@@ -44,6 +44,8 @@ Those spline are the one used to flatten the distribution and are the calibratio
**If you want to improve the Z resolution** : The IC 2-3 are well behaved but the 4 has a factor 2 in the analysis.
Therefore one should be very careful about the range in which the spline is made for IC4 !!
This is a point to be improved !
+The range in Y should be custom for each section of the chio. One should take time to code something to make this range finding
+automatic, for instance project the ratio in a histogram then fitting it by a gaussian and take the range as mean +- 3 * sigma.
# Step 3 : Correction of DE in Y
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/AlignIC.C b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/AlignIC.C
index cb9376140c3db7577f59f867a425a9f2da2f3469..812e7203637c81f313f039fab32d8b5410b6ff9c 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/AlignIC.C
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/AlignIC.C
@@ -28,10 +28,11 @@ void HistoDrawer(vector<TH2F*> h ,const char* CharName);
// Range of histo
-double XMin[5]= {-520,-1000,-520,-1000,-1000};
-double XMax[5]= {380,1000,380,1000,1000};
-double YMin[5]= {2000,1000,0,0,2};
-double YMax[5]= {9500,10000,3,3,4};
+double XMin[5]= {-520,-1000,-520,-800,-800};
+double XMax[5]= {380,1000,380,800,800};
+double YMin[5]= {2000,1000,0,0.8,2.};
+double YMax[5]= {9500,10000,3,1.5,2.4};
+double CondX[2]= {-200,200};
//////////////////////////////////////////////////////////////////////////////////////////
void AlignIC() {
@@ -77,7 +78,7 @@ void AlignIC() {
vector<TSpline3*> SplineICurrent;
// X Spline
- SplineICurrent.push_back(MakeSpline(hIC[seg][0] , NCallSpline, XMin[2] , XMax[2], YMin[2], YMax[2]));
+ SplineICurrent.push_back(MakeSpline(hIC[seg][0] , NCallSpline, XMin[2] , XMax[2], YMin[4], YMax[4]));
SplineICurrent.at(0)->SetName(Form("fspline_IC%d_X",seg));
NCallSpline+=1;
@@ -96,11 +97,12 @@ void AlignIC() {
//Load Histo
hIC.at(seg) = (HistoFillerIC(seg,SplineAllIC));
+ cout << hIC.at(seg).at(0)->GetYaxis()->GetXmin() << endl;
//Reset spline holder
vector<TSpline3*> SplineICurrent;
// X Spline
- SplineICurrent.push_back(MakeSpline(hIC[seg][0] , NCallSpline, XMin[2] , XMax[2], YMin[2], YMax[2]));
+ SplineICurrent.push_back(MakeSpline(hIC[seg][0] , NCallSpline, XMin[2] , XMax[2], YMin[3], YMax[3]));
SplineICurrent.at(0)->SetName(Form("fspline_IC%d_X",seg));
NCallSpline+=1;
@@ -159,7 +161,9 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
// Input and setters
TChain* chain = new TChain("PhysicsTree");
+ chain->Add("../../../../root/analysis/VamosCalib246.root");
chain->Add("../../../../root/analysis/VamosCalib247.root");
+ chain->Add("../../../../root/analysis/VamosCalib248.root");
TICPhysics* IC = new TICPhysics();
TTimeData *Time = new TTimeData();
@@ -222,12 +226,12 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
else {
hIC_Y = new TH2F(Form("hChio_IC%d_IC%d_Y",segment,segment-1),
- Form("hChio_IC%d_IC%d_Y",segment,segment-1), 1000, -1000, 1000, 500, 0, 4);
+ Form("hChio_IC%d_IC%d_Y",segment,segment-1), 1000, -1000, 1000, 500, YMin[3], YMax[3]);
hIC_X = new TH2F(Form("hChio_IC%d_IC%d_X",segment,segment-1),
- Form("hChio_IC%d_IC%d_X",segment,segment-1), 1000, -800, 800, 500, 0, 4);
+ Form("hChio_IC%d_IC%d_X",segment,segment-1), 1000, -800, 800, 500, YMin[3], YMax[3]);
}
-
+ hIC_Y->GetYaxis()->SetCanExtend(kFALSE);
//===========================================================================================================
// Beginning loop on entries
//===========================================================================================================
@@ -242,7 +246,7 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
vector<double> IC_Spline_CorrX(NSpline + 1),IC_Spline_CorrY(NSpline+1);
- if (Time->GetMWPC13Mult() ==1 && IC->fIC_TS.size()>= (segment+1)){ //only mult 1 event
+ if (Time->GetMWPC13Mult() ==1 && IC->fIC_TS.size()>= (segment+1) && FF_IC_X >-530 ){ //only mult 1 event
UShort_t FPMW_Section = Time->GetSection_MWPC3(0);
@@ -299,6 +303,7 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
} // end loop on entries
+ cout << endl;
//===========================================================================================================
// Output
//===========================================================================================================
@@ -380,9 +385,9 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
else {
hIC_Y.at(pseg) = new TH2F(Form("hChiocorr_IC%d_IC%d_Y",pseg,pseg-1),
- Form("hChiocorr_IC%d_IC%d_Y",pseg,pseg-1), 1000, -1000, 1000, 500, 0, 4);
+ Form("hChiocorr_IC%d_IC%d_Y",pseg,pseg-1), 1000, -1000, 1000, 500, YMin[3], YMax[3]);
hIC_X.at(pseg) = new TH2F(Form("hChiocorr_IC%d_IC%d_X",pseg,pseg-1),
- Form("hChiocorr_IC%d_IC%d_X",pseg,pseg-1), 1000, -800, 800, 500, 0, 4);
+ Form("hChiocorr_IC%d_IC%d_X",pseg,pseg-1), 1000, -800, 800, 500, YMin[3], YMax[3]);
}
}
@@ -454,6 +459,7 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
<< " | Estimated time left: " << timeLeft << " seconds" << "\r" << flush;
}
} // end loop on entries
+
} //end only m2
//===========================================================================================================
// Output
@@ -462,6 +468,7 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
hIC.push_back(hIC_X);
hIC.push_back(hIC_Y);
+ cout << endl;
return hIC;
} // End HistoFiller
@@ -497,104 +504,136 @@ TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin = 0 , double
//===========================================================================================================
// Create the TProfile
+ //hInput->GetYaxis()->SetRangeUser(YMin+0.01,YMax-0.01);
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;
+ int FirstBin, LastBin;
+ double Threshold = hInput->GetEntries() / hInput->GetNbinsX() / 10;
+ for (int bin =1 ; bin<hInput->GetNbinsX(); bin++) {
+
+ int nEvents = 0;
+ for (int ybin = 1; ybin <= hInput->GetNbinsY(); ybin++) {
+ nEvents += hInput->GetBinContent(bin, ybin);
+ }
+ FirstBin = bin;
+ if (nEvents> Threshold && hProfile->GetXaxis()->GetBinLowEdge(bin) > XMin ) break;
+ }
+
+ for (int bin = hInput->GetNbinsX(); bin>1 ; bin--) {
+ int nEvents = 0;
+ for (int ybin = 1; ybin <= hInput->GetNbinsY(); ybin++) {
+ nEvents += hInput->GetBinContent(bin, ybin);
+ }
+ LastBin = bin;
+ if (nEvents > Threshold && hProfile->GetXaxis()->GetBinUpEdge(bin) < XMax ) break;
+ }
+
+ cout << "BinEdges " << FirstBin << " " << LastBin << endl;
+ //===========================================================================================================
+ // 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);
+
+
+ double parpol3[4];
+ double parpol1[2];
+ // find the function to extend the TProfile on the lower range
+ TF1 * fitpol3 = new TF1(Form("FitPol3_pfx_%02d",Int_t(NCallee/2)),"pol3",hProfile->GetBinLowEdge(FirstBin),hProfile->GetBinLowEdge(LastBin));
+ hProfile->Fit(fitpol3,"RQ");
+ fitpol3->GetParameters(parpol3);
+
+ // find the function to extend the TProfile on the higher range
+ TF1 * fitpol1 = new TF1(Form("FitPol1_pfx_%02d",Int_t(NCallee/2)),"pol1",hProfile->GetBinLowEdge(FirstBin),hProfile->GetBinLowEdge(LastBin));
+ hProfile->Fit(fitpol1,"RQ");
+ fitpol1->GetParameters(parpol1);
+
+
+ //===========================================================================================================
+ // Fill extended profile
+ //===========================================================================================================
+
+ // 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<4; par++) newval += parpol3[par]*pow(hProfile->GetBinCenter(bin),par);
+ pfx->SetBinContent(bin,newval);
+ }
+ pfx->Draw("same");
+
+ gspline = new TSpline3(pfx);
+ gspline->SetName(Form("fspline_%d", NCallee + 1));
+
+ return gspline;
} // end makespline
void HistoDrawer(vector<TH2F*> h,const char* CharName){
- int NHisto = static_cast<int>(h.size());
+ int NHisto = static_cast<int>(h.size());
- TCanvas* Can= new TCanvas(CharName,CharName,0,0,2000,1000);
+ TCanvas* Can= new TCanvas(CharName,CharName,0,0,2000,1000);
- Can->Divide(NHisto);
+ Can->Divide(NHisto);
- for ( int i=0 ; i < NHisto ; i+=1 ){
- //All x are even and all y are odd
- int CanvaNumber = (i) ;
- Can->cd(CanvaNumber);
- h[i]->Draw("colz");
- }
+ for ( int i=0 ; i < NHisto ; i+=1 ){
+ //All x are even and all y are odd
+ int CanvaNumber = (i) ;
+ Can->cd(CanvaNumber);
+ h[i]->Draw("colz");
+ }
}
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;
+ 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;
}
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/IC0_XYCorrection.C b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/IC0_XYCorrection.C
index 8e50b0faaa9311bb897bd1ef3eafcaf011960bd5..7a3b3e499cba8245144653f0db520c3e971ac436 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/IC0_XYCorrection.C
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step2/IC0_XYCorrection.C
@@ -17,6 +17,8 @@ void IC0_XYCorrection(){
//===========================================================================================================
TChain* chain = new TChain("PhysicsTree");
chain->Add("../../../../root/analysis/VamosCalib247.root");
+ chain->Add("../../../../root/analysis/VamosCalib246.root");
+ chain->Add("../../../../root/analysis/VamosCalib248.root");
TTimeData *Time = new TTimeData();
TICPhysics* IC = new TICPhysics() ;
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step3/DEYcorrection.C b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step3/DEYcorrection.C
index 051cd11702e97e2f7e9e78d6c96f80499ff3ddd1..8bbb57d38ccfd7b6e19ae9c8d791f0fd2956ff60 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step3/DEYcorrection.C
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step3/DEYcorrection.C
@@ -27,6 +27,8 @@ void DEYcorrection(bool cut = false, bool spline = false) {
//===========================================================================================================
TChain* chain = new TChain("PhysicsTree");
chain->Add("../../../../root/analysis/VamosCalib247.root");
+ chain->Add("../../../../root/analysis/VamosCalib246.root");
+ chain->Add("../../../../root/analysis/VamosCalib248.root");
//chain->Add("../../../root/analysis/Run246.root");
TICPhysics* IC = new TICPhysics() ;
@@ -247,8 +249,8 @@ void DEYcorrection(bool cut = false, bool spline = false) {
for (int seg = 0 ; seg < sizeof(IC->fIC_PID)/sizeof(IC->fIC_PID[0]) ; seg++ ){
if (seg >4){
- if (seg == 5) E += double(IC->fIC_PID[seg]);
- else E +=double(IC->fIC_PID[seg]) ;
+ //E +=double(IC->fIC_PID[seg]) ;
+ E += ICcorr_Y.at(seg);
}
}
@@ -383,7 +385,6 @@ void DEYcorrection(bool cut = false, bool spline = false) {
hDE02_Y->Fill(FF_DriftTime,DE02);
hDE023_Y->Fill(FF_DriftTime,DE023);
hDE0234_Y->Fill(FF_DriftTime,DE0234);
-
}
}
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/CutAutoDEE_FineTuning.C b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/CutAutoDEE_FineTuning.C
new file mode 100644
index 0000000000000000000000000000000000000000..f8652bd7a7e080d7e205cdb68ec8405cf1430caa
--- /dev/null
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/CutAutoDEE_FineTuning.C
@@ -0,0 +1,698 @@
+#include "CutAutoDEE.h"
+#include <numeric>
+#include <set>
+#include <iostream>
+#include <unordered_map>
+
+// Important variable
+#define STEP 2; // Step to parse through the histogram in X
+#define BINSIZE 10; // Size of the bin in X on which to projectY the plot
+#define THRESHOLDY 75.0;
+#define THRESHOLDX 200.0;
+
+//Main
+void CutAutoDEE(){
+
+ // Retrieve initial histogram
+ TFile *inFile = new TFile("Output/DE_E_merged.root","read");
+ TH2F *hDE_E = (TH2F*)inFile->Get("DE_E");
+ // Inverse data to use tspectrum
+ for (int xBin=1; xBin<= hDE_E->GetNbinsX();xBin++){
+ for (int yBin=1; yBin<= hDE_E->GetNbinsY();yBin++){
+ hDE_E->SetBinContent(xBin,yBin , -hDE_E->GetBinContent(xBin,yBin) );
+ }
+ }
+
+ //Store result of peak finding
+ //This is not a good way to do this : We should normalize the structure of the code by
+ //directly using a PeakData class to store each of the peak.
+ vector<int> *vPeaks = new vector<int>() ;
+ vector<TSpectrum*> *vPeakFinder = new vector<TSpectrum*>();
+ vector<vector<Double_t>> *LinePos = new vector<vector<Double_t>> ;
+
+ // Tspectrum used to find the peak : If you have problem with peak finding
+ // adjust the setting in this function
+ // Step and Binsize are used in this function
+ Double_t sigma = 3;
+ Double_t threshold = 0.0001;
+ PeakFinder(vPeakFinder,vPeaks,hDE_E,LinePos,sigma,threshold);
+
+
+ // Link the peak, this is where the most critical variable are.
+ // Specifically the maximum distance between each peak and the minimum number of peak in a line
+ vector<TCutG*> CutLine = PeakLinker(LinePos);
+
+ // Inverse data for plotting
+ for (int xBin=1; xBin<= hDE_E->GetNbinsX();xBin++){
+ for (int yBin=1; yBin<= hDE_E->GetNbinsY();yBin++){
+ hDE_E->SetBinContent(xBin,yBin , -hDE_E->GetBinContent(xBin,yBin));
+ }
+ }
+
+ //Plot
+ gStyle->SetPalette(kRainBow);
+ TCanvas* c2 = new TCanvas("c2","c2");
+ hDE_E->Draw("colz");
+
+ TCanvas* c3 = new TCanvas("c3","c3");
+ hDE_E->Draw("colz");
+ for (auto cut : CutLine){
+ cut->Draw("SAME");
+ }
+
+ //Save the cut !
+ TFile *ofile = new TFile("Output/CutAutoZ.root","RECREATE");
+ for (auto cut : CutLine){
+ cut->Write();
+ }
+ ofile->Close();
+}
+
+/**
+ * Find the peak of an histogram by subdivinding the histogram on it's x axis.
+ * Store the number of peak and their position in different vector
+ *
+ * input and Output:
+ * vector<TSpectrum*> *sPeak : Vector of TSpectrum to output
+ * vector<int> *vPeak : Vector of number of peak for each bin range
+ * vector<vector<Double_t>> *LinePos : Vector of position of each peak for each bin range
+ */
+void PeakFinder(vector<TSpectrum*> *sPeak , vector<int> *vPeak, TH2F *h, vector<vector<Double_t>> *LinePos, Double_t sigma, Double_t threshold){
+
+ // Initialising
+ int nbinsX = h->GetNbinsX();
+ int iter = 0;
+ int stepParsing = STEP;
+ int SizeBin = BINSIZE;
+
+ // Looping on the histogram x axis.
+ for (int i = 1; i <= nbinsX; i += stepParsing) {
+
+ // Limit in binnumber
+ iter +=1;
+ int x1 = i;
+ int x2 = std::min(i + SizeBin - 1, nbinsX);
+
+ // Projection on Y between the binnumber
+ TString projName = Form("projY_%d_to_%d", x1, x2);
+ TH1D* projY = h->ProjectionY(projName, x1, x2);
+ projY->SetTitle(Form("Projection Y: bins %d to %d", x1, x2));
+
+ // PeakFinder
+ TSpectrum *PeakFinder = new TSpectrum(600);
+ int nPeaks = PeakFinder->Search(projY,sigma,"",threshold);
+
+ //Store npeaks and tspectrum
+ vPeak->push_back(nPeaks); // number of peak for this bin range
+ sPeak->push_back(PeakFinder);
+
+ //Save the position of each peak
+ Double_t * x = PeakFinder->GetPositionX();
+
+ for(int peaks = 0 ; peaks < nPeaks ; peaks++){
+ vector<Double_t> Position(2);
+ Position.at(0) = (double(i) + double(SizeBin)/2.0) * h->GetXaxis()->GetBinWidth(4) ; //Retrieve E position as the center
+ //of the bin
+ Position.at(1) = PeakFinder->GetPositionX()[peaks]; // Position in DE
+ LinePos->push_back(Position);
+ }
+
+
+ //Test
+ if (true && iter ==50){
+ // Draw the histogram
+ TCanvas* c = new TCanvas(Form("c%d",i), Form("Histogram with E = %f",(h->GetXaxis()->GetBinCenter(i))), 800, 600);
+ projY->Draw();
+ }
+
+ }
+}
+
+/**
+ * Link the different point to create TCUT from the peak found previously
+ * To do so a max range is given and then all point separated by a distance inferior to it are linked.
+ *
+ * Input :
+ * vec<vec<Double>> *LinePos : pointer to a vector of length number of peak holding position in x and y for each peak
+ * Output:
+ * TCUTG : a cut for each charge in vamos
+ */
+vector<TCutG*> PeakLinker(vector<vector<Double_t>> *LinePos){
+
+ //Retrieve position of peaks
+ vector<Double_t> X , Y ;
+ for (int ipeak = 0; ipeak<LinePos->size() ; ipeak ++){
+ X.push_back(LinePos->at(ipeak).at(0));
+ Y.push_back(LinePos->at(ipeak).at(1));
+ }
+
+ //===========================================================================================================
+ // 1. Find the threshold between two point
+ //===========================================================================================================
+
+ // Find the distance between two point in X.
+ //Sort the vec
+ vector<Double_t> XSub(X);
+ sort(XSub.begin(),XSub.end());
+
+ //Take only the unique element
+ set<Double_t> SetX(XSub.begin() , XSub.end());
+ vector<Double_t> XThresh ;
+ for (Double_t elem : SetX) {
+ XThresh.push_back(elem);
+ }
+
+ //Find the average dist in X
+ vector<Double_t> DiffX = diff(XThresh);
+ Double_t ThresholdXestim ;
+ Double_t ThresholdX ;
+ for (int i =0 ; i<DiffX.size() ; i++) { ThresholdXestim += DiffX.at(i);}
+ ThresholdXestim = ThresholdXestim / DiffX.size();
+
+
+ // In Y the threshold is not easily findable, However it's critical that the threshold
+ // is superior to the average distance in X. If not line can't be found by construction
+ // The threshold should also be inferior to the distance between each line in Y
+ // To do so adjust the bin STEP at the beginning of the code.
+
+ //Very important set
+ ThresholdX = THRESHOLDX;
+ Double_t ThresholdY = THRESHOLDY;
+ cout << "Threshold X : " << ThresholdX << endl;
+ cout << "Threshold X approx : " << ThresholdXestim << endl;
+ cout << "N of X " << SetX.size() << endl;
+ cout << "ThresholdY " << ThresholdY << endl;
+
+ //===========================================================================================================
+ // 2.LINK POINT
+ //===========================================================================================================
+
+ // Now that the threshold is found we will link the point by comparing each point
+ // Notice that if two point have the same X they can't be linked
+
+ // Create a vector holding the id of each peak
+ int ID = 0 ; // Initial Id of the first line
+ vector<int> vID(LinePos->size(), -1000);
+
+ // We want to loop on the UNIDENTIFIED peak.
+ for (int ipeak = 0 ; ipeak < LinePos->size() ; ipeak ++){
+
+ //Here do ID manipulation
+ //Create a set of ID without the unidentified peak
+ set<int> SetId(vID.begin() , vID.end());
+ if (SetId.size() > 1 ) {SetId.erase(-1000) ; ID = *prev(SetId.end()) + 1 ;}
+
+ //If the current peak is identified set newID to it's ID
+ if (vID.at(ipeak) != -1000) ID = vID.at(ipeak); // If the point is identified propagate it!
+
+ Double_t XPosPeak = LinePos->at(ipeak).at(0);
+ Double_t YPosPeak = LinePos->at(ipeak).at(1);
+
+
+ // Parse through all the next peak
+ // This could be done faster by sorting the list by X and Y before and stoping after a threshold
+ // For now it'll do just fine
+ for (int ipeakCompare = ipeak + 1 ; ipeakCompare < LinePos->size() ; ipeakCompare ++){
+
+ if ((vID.at(ipeakCompare) != -1000) || (XPosPeak == LinePos->at(ipeakCompare).at(0))) continue; // If the point is identified skip it !
+
+ //Calculate the distance between the two peak
+ Double_t XPosPeakComp = LinePos->at(ipeakCompare).at(0);
+ Double_t YPosPeakComp = LinePos->at(ipeakCompare).at(1);
+ Double_t dist = sqrt( pow(abs(XPosPeak-XPosPeakComp),2) + pow(abs(YPosPeak-YPosPeakComp),2) );
+ Double_t distX = sqrt( pow(abs(XPosPeak-XPosPeakComp),2) );
+ Double_t distY = sqrt( pow(abs(YPosPeak-YPosPeakComp),2) );
+
+ //If the distance is inferior to the threshold then set both ID to the new one !
+ if ((distX < ThresholdX )&& (distY <ThresholdY )) {
+ vID.at(ipeak) = ID ;
+ vID.at(ipeakCompare) = ID ;
+ } // end threshold if
+ } // end compare
+ } // end loop peak
+
+ //Check if needed
+ /*
+ for (int i=0 ; i<vID.size() ; i++){
+ if (vID.at(i)!=-1000){
+ cout << vID.at(i) << " , ";
+ }
+ }
+ */
+
+
+
+
+ // Here we flush all the id that have less element than a threshold value
+ //===========================================================================================================
+ // 3. Flush line with not enough value
+ //===========================================================================================================
+
+ //This is a pretty sensitive variable
+ int MinElem = SetX.size()/4;
+ cout << "Min Elem " << MinElem << endl ;
+
+
+ set<int> SetIdTemp(vID.begin() , vID.end());
+ for (int Id : SetIdTemp){
+ int NId = 0 ;
+ for (int ipeak = 0 ; ipeak < vID.size() ; ipeak ++) {
+ if (vID.at(ipeak) == Id) NId += 1 ;
+ }
+ if ( NId < MinElem ) {
+
+ for (int ipeak = 0 ; ipeak < vID.size() ; ipeak ++) {
+ if (vID.at(ipeak) == Id) vID.at(ipeak) = -1000 ;
+ }
+
+ }
+ }
+
+ //===========================================================================================================
+ // 4. TCutG Creator
+ //===========================================================================================================
+
+ // We now have identified each peak to a id that correspond to a line in the DE E
+ // Now we need to link the ADJACENT line and then transform those point into tcutg
+
+ set<int> SetId(vID.begin() , vID.end()); // This is all the unique element of vID
+ vector<int> IndexPeak(LinePos->size());
+ for (int ipeak = 0 ; ipeak < LinePos->size() ; ipeak ++){
+ IndexPeak.at(ipeak) = ipeak ;
+ }
+
+ //===========================================================================================================
+ // 4.1 Sort by id
+ //===========================================================================================================
+ // The next part of the code SHOULD sort by meanY but does not.
+ // It therefore only serve as a convoluted way to sort by ID.
+ // We need to sort by id either way so no need to rewrite i
+
+ vector<Double_t> MeanYId(SetId.size()) ;
+ vector<Double_t> NelemYId(SetId.size()) ;
+
+ //Loop on the element of the set
+ for (int Id = 0 ; Id < SetId.size() ; Id ++){
+ Double_t Sum = 0 , Npeak = 1;
+ auto itSet = SetId.begin();
+ advance(itSet,Id);
+
+ //Loop on all the peak
+ for (int ipeak = 0 ; ipeak < LinePos->size() ; ipeak ++){
+
+ // If the ID of the peak correspond to one of the set then increment his Y and NPeak value
+ if (vID.at(ipeak) == *itSet )
+ {
+ Npeak += 1 ;
+ Sum += LinePos->at(ipeak).at(1);
+ } // End if Id
+ } // Loop on Peak
+
+ NelemYId.at(Id) = Npeak ;
+ MeanYId.at(Id) = Sum /Npeak ;
+ // cout << Id << " " << MeanYId.at(Id) << endl;
+ } // Loop on Id
+
+ //Now we will sort the MeanY vector and return a vector holding the ID sorted
+ vector<size_t> IndexMeanY(MeanYId.size());
+ iota(IndexMeanY.begin(),IndexMeanY.end(),0);
+ sort(IndexMeanY.begin(),IndexMeanY.end() , [&](size_t i, size_t j) {
+ return MeanYId.at(i) < MeanYId.at(j);
+ });
+
+ // Now we have to sort the SetId according to this index ordering
+ // Since set can't be sorted we need to copy it to a vector
+ vector<int> vSetId(SetId.begin(), SetId.end());
+ vector<int> SortedSetId(SetId.begin(), SetId.end());
+ for (int Index = 0 ; Index < IndexMeanY.size() ; Index ++){ SortedSetId.at(IndexMeanY.at(Index)) = vSetId.at(Index) ;}
+
+ /*
+ cout << endl << " SORTED" <<endl ;
+ for (int i = 0 ; i<IndexMeanY.size() ; i++){
+ if (IndexMeanY.at(i) == 49 || IndexMeanY.at(i) == 32){
+ cout << NelemYId.at(IndexMeanY.at(i)) << endl;
+ }
+ cout << IndexMeanY.at(i) << " " << SortedSetId.at(i) << " " << MeanYId.at(IndexMeanY.at(i)) << endl;
+ }
+
+ cout << endl << " *********************" <<endl ;
+ */
+
+ //===========================================================================================================
+ // 4.2 Recompacting data into a struct
+ //===========================================================================================================
+
+
+ // We now have the different lines sorted by Id we just need to link them together by group of two.
+ // to do so we will sort each id by x and then link them
+
+ // 1st step is to sort the original IndexPeak according to the MeanY sort
+ // To do so we need to sort the vID vector according to meanY and apply this transformation to IndexPeak too
+
+ // vID as the ID of each peak -> correspondance of peak to a line
+ // IndexPeak -> index of each of these peak of size number of peak
+ // SortedSetId -> Set of sorted ID of peak according to mean y
+
+ //Create a map for the sort of ID
+ unordered_map<int,int> order_map ;
+ for (size_t i = 0 ; i < SortedSetId.size() ; i++ ) {
+ order_map[SortedSetId[i]] = i ;
+ }
+
+ // Combine Index and ID into a single vector of pair
+ vector<pair<int,int>> PeakId ;
+ for(size_t i = 0 ; i < IndexPeak.size() ; i++){
+ PeakId.emplace_back(IndexPeak.at(i), vID.at(i));
+ }
+
+ // Sort the PeakId on the custom order of IDs
+ sort(PeakId.begin(), PeakId.end(), [&order_map](const auto& a, const auto& b){
+ return order_map[a.second] < order_map[b.second];
+ } );
+
+ //2nd step is to take the original Position vector and sort it thanks to PeakId
+
+ vector<vector<Double_t>> LinePosSorted = *LinePos ;
+
+ //Swap LinePosSorted elemn according to Index of PeakID
+ for (int i = 0 ; i < PeakId.size() ; i++ ){
+ LinePosSorted.at(i) = LinePos->at(PeakId.at(i).first);
+ }
+
+ // 3rd step is to compile x y id and index in one structure
+
+ if (PeakId.size() != LinePos->size() ){
+ cerr << "Vector not same size" << endl;
+ }
+
+ vector<PeakData> PeakStruct ;
+ for (size_t i=0 ; i<LinePosSorted.size() ; i++){
+ PeakData TempPeak = {LinePosSorted.at(i).at(0), LinePosSorted.at(i).at(1) , PeakId.at(i).second , PeakId.at(i).first};
+ PeakStruct.push_back(TempPeak);
+ }
+
+ //===========================================================================================================
+ // 4.3 Cleanupi and sortY
+ //===========================================================================================================
+ // We drop all non clustered element and sort the whole struct by the mean value of Y for their ID.
+
+ // 4th step is to drop all element = -1000
+ //
+
+ auto groupBegin = find_if(PeakStruct.begin(),PeakStruct.end(),[&](const PeakData& elem){
+ return elem.LineId == -1000 ;
+ });
+
+ auto groupEnd = find_if(groupBegin,PeakStruct.end(),[&](const PeakData& elem){
+ return elem.LineId != -1000 ;
+ });
+
+ PeakStruct.erase(groupBegin,groupEnd);
+
+ //4.1 step we will resort by mean y because it didnt work before , probably due to data structure
+
+ auto ItMeanY = PeakStruct.begin();
+ vector<int> SortedId ;
+ vector<Indexed2D> MeanY ;
+ int index = 0;
+
+ while (ItMeanY != PeakStruct.end()){
+
+ //Get position of end of the current id
+ auto IdEnd = find_if(ItMeanY, PeakStruct.end() , [&](const PeakData& elem){
+ return elem.LineId != ItMeanY->LineId;
+ });
+
+ // Fill MeanY of the current Id
+
+ Double_t Sum = 0;
+ Double_t Npts = 0;
+
+ for (auto i = ItMeanY ; i!=IdEnd ; i++){
+ Sum += i->Y;
+ Npts += 1;
+ }
+
+ MeanY.push_back({Sum/Npts ,ItMeanY->LineId , index });
+
+ index +=1;
+ ItMeanY = IdEnd;
+ }
+
+ // We now have a all the mean Y of each ID and their index
+ // Let's now sort the id from this meanY
+
+ sort(MeanY.begin(),MeanY.end(),[&](Indexed2D a , Indexed2D b){
+ return a.x < b.x;
+ });
+
+ unordered_map<int,size_t> ordered_map ;
+ for (size_t i=0 ; i<MeanY.size() ; i++){
+ ordered_map[MeanY.at(i).id]=i;
+ }
+
+
+ sort(PeakStruct.begin(),PeakStruct.end(), [&ordered_map](const PeakData& a , const PeakData &b) {
+ return ordered_map[a.LineId]<ordered_map[b.LineId];
+ });
+
+ //===========================================================================================================
+ // 4.4 Sort by X alternating order
+ //===========================================================================================================
+ // We sort by X in alternating order because a TCutG need the point to be aligned to link them together
+
+ // Iterate and sort groups with the same LineId
+ auto it = PeakStruct.begin();// create an iterator that point to the elements of the PeakData
+ bool Oddness = true; // to know in which way sort
+ while (it != PeakStruct.end()){
+
+ // Find range where id is the same
+ // This return the first iterator where the condition is fulfilled
+ auto groupEnd = find_if(it,PeakStruct.end(),[lineId = it->LineId](const PeakData& elem){
+ return elem.LineId != lineId ;
+ });
+
+ // Now we sort this whole group by x
+
+ sort(it,groupEnd,[Oddness](const PeakData& a , const PeakData& b){
+ if (Oddness) return a.X < b.X ;
+ if (!Oddness) return a.X > b.X ;
+ else return a.X < b.X;
+ });
+
+ // Move to next group
+ it = groupEnd;
+ Oddness = not(Oddness) ;
+ }
+
+ // works !
+ /*
+ for (int i : SortedSetId) {cout << i << endl;}
+
+ int tempId= -1000;
+ for (PeakData i : PeakStruct) {
+ if (tempId != i.LineId ){
+ tempId = i.LineId;
+ cout << " Id " << i.LineId << " Y " << i.Y << endl;
+ }
+ }
+ */
+ //===========================================================================================================
+ // 4.5 Create TCut
+ //===========================================================================================================
+ // We need to make sublist for each 2 lines
+
+ // We will iterate through our data and create substructure for each adjacent pair of LineID
+ // and store them in a vector of TCutG
+
+ vector<TCutG*> Res ;
+
+ it = PeakStruct.begin();
+ int iteration = 0;
+
+
+ while(it != PeakStruct.end()){
+
+ // 1st line
+ int LineIdBottom = it->LineId;
+ auto BottomLine = find_if(it,PeakStruct.end(), [&] (const PeakData& elem){
+ return elem.LineId != LineIdBottom ;
+ }) ;
+
+ //2nd line
+ int LineIdTop = BottomLine->LineId;
+ auto TopLine = find_if(BottomLine,PeakStruct.end(), [&] (const PeakData& elem){
+ return elem.LineId != LineIdTop ;
+ }) ;
+
+ //cout << LineIdBottom << " " << LineIdTop << endl;
+
+ // first we need to find the shortest line
+ int SizeBottom = distance(it,BottomLine);
+ int SizeTop = distance(BottomLine,TopLine);
+
+ //Subvector with the longest line first then the sortest
+ vector<PeakData> CutVec;
+ vector<PeakData> LongVec;
+ vector<PeakData> ShortVec;
+
+ if(SizeBottom >= SizeTop ){
+ CutVec.insert(CutVec.end(),it,TopLine);
+ LongVec.insert(LongVec.end(),it,BottomLine);
+ ShortVec.insert(ShortVec.end(),BottomLine,TopLine);
+ }
+
+ else {
+ CutVec.insert(CutVec.end(),BottomLine,TopLine);
+ CutVec.insert(CutVec.end(),it,BottomLine);
+ LongVec.insert(LongVec.end(),BottomLine,TopLine);
+ ShortVec.insert(ShortVec.end(),it,BottomLine);
+ }
+
+
+
+ //nbr of point
+ int Npoints = CutVec.size();
+ //Close the loop by finding the closest point to the last one
+ int closestIndexLast = 0;
+ int closestIndexBegin = LongVec.size()-1;
+ int SizeArray = 0;
+ if (ShortVec.size() == 0) {
+ closestIndexLast = CutVec.size();
+ SizeArray = CutVec.size();
+ }
+ else{
+ PeakData lastPoint = ShortVec.back();
+ double minDistance = distancePeak(lastPoint,LongVec.at(0));
+ for (size_t i = 1 ; i < LongVec.size(); i++){
+ double dist = distancePeak(lastPoint,LongVec.at(i));
+ if(dist<minDistance){
+ minDistance = dist;
+ closestIndexLast = i;
+ }
+ }
+ PeakData firstPoint = ShortVec.at(0);
+ double minDistanceFirst = distancePeak(firstPoint,LongVec.at(0));
+ for (size_t i = 1 ; i < LongVec.size(); i++){
+ double dist = distancePeak(firstPoint,LongVec.at(i));
+ if(dist<minDistanceFirst){
+ minDistanceFirst = dist;
+ closestIndexBegin = i;
+ }
+ }
+ SizeArray = ShortVec.size() + (closestIndexBegin - closestIndexLast);
+ }
+ //Create array to hold position
+ double *x = new double[SizeArray + 1 ];
+ double *y = new double[SizeArray + 1 ];
+
+ if (SizeTop == 0){
+
+ for(int i = 0 ; i<CutVec.size() ; i++) {
+ x[i] = CutVec.at(i).X ;
+ y[i] = CutVec.at(i).Y ;
+ }
+ x[SizeArray ] = CutVec.at(0).X;
+ y[SizeArray ] = CutVec.at(0).Y;
+
+ }
+
+ else {
+ for(int i = 0 ; i<ShortVec.size() ; i++) {
+ x[i] = ShortVec.at(i).X ;
+ y[i] = ShortVec.at(i).Y ;
+ }
+
+ for(int i = closestIndexLast ; i<=closestIndexBegin ; i++) {
+ x[i + ShortVec.size() - closestIndexLast] = LongVec.at(i).X ;
+ y[i + ShortVec.size() - closestIndexLast ] = LongVec.at(i).Y ;
+ }
+
+
+
+ x[SizeArray] = ShortVec.at(0).X;
+ y[SizeArray] = ShortVec.at(0).Y;
+
+ }
+ //Saving and moving on
+
+ if (Npoints >10 && (TopLine != PeakStruct.end())){ //Dont link last line to itself
+ TCutG* Cut = new TCutG(Form("Cut %d",iteration) , SizeArray+1, x , y);
+ Cut->SetLineWidth(2);
+ Cut->SetLineColor(kRed);
+ Res.push_back(Cut);
+ }
+ iteration += 1 ;
+ it = BottomLine ;
+ } // end while
+
+ return Res;
+}
+
+
+vector<Double_t> diff(const vector<Double_t>& input){
+ vector<Double_t> result;
+ if (input.size() < 2) return result; // Return empty vector if size < 2
+
+ for (size_t i = 1; i < input.size(); ++i) {
+ result.push_back(input[i] - input[i - 1]);
+ }
+ return result;
+}
+
+void SortAndRenameTCutG() {
+ TIter next(gROOT->GetListOfSpecials());
+ std::vector<std::pair<int, TCutG*>> cuts;
+
+ TObject *obj;
+ while ((obj = next())) {
+ TCutG *cut = dynamic_cast<TCutG*>(obj);
+ if (cut) {
+ std::string name = cut->GetName();
+ if (name.rfind("Cut", 0) == 0) { // Check if name starts with "Cut"
+ try {
+ int num = std::stoi(name.substr(3));
+ cuts.emplace_back(num, cut);
+ } catch (...) {
+ continue; // Ignore malformed names
+ }
+ }
+ }
+ }
+
+ // Sort by extracted number
+ std::sort(cuts.begin(), cuts.end(), [](const auto &a, const auto &b) {
+ return a.first < b.first;
+ });
+
+ // Rename sequentially
+ for (size_t i = 0; i < cuts.size(); ++i) {
+ cuts[i].second->SetName(Form("Cut %zu", i));
+ }
+
+ std::cout << "TCutG objects renamed from original numbering to 0:N.\n";
+}
+void SaveTCutGToFile() {
+ TFile file("Output/CutAutoZ.root", "RECREATE");
+ if (!file.IsOpen()) {
+ std::cerr << "Error: Could not open Output/CutAutoZ.root for writing.\n";
+ return;
+ }
+
+ TIter next(gROOT->GetListOfSpecials());
+ TObject *obj;
+ while ((obj = next())) {
+ TCutG *cut = dynamic_cast<TCutG*>(obj);
+ if (cut) {
+ cut->Write();
+ }
+ }
+
+ file.Close();
+ std::cout << "All TCutG objects saved to Output/CutAutoZ.root.\n";
+}
+
+double distancePeak(const PeakData& p1, const PeakData& p2) {
+ return abs(p1.X - p2.X);
+}
+
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/CutAutoDEE_FineTuning.h b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/CutAutoDEE_FineTuning.h
new file mode 100644
index 0000000000000000000000000000000000000000..effc07e99a8b0b397326bc1846f69059dee66611
--- /dev/null
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/CutAutoDEE_FineTuning.h
@@ -0,0 +1,36 @@
+#include <TCanvas.h>
+#include <TCutG.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TH2.h>
+#include <TSpectrum.h>
+#include <vector>
+#include <TStyle.h>
+
+using namespace std;
+
+void PeakFinder(vector<TSpectrum*> *sPeak , vector<int> *vPeak, TH2F *h , vector<vector<Double_t>> *LinePos, Double_t sigma, Double_t threshold);
+vector<TCutG*> PeakLinker(vector<vector<Double_t>> *LinePos);
+
+struct PeakData{
+ Double_t X ;
+ Double_t Y ;
+ int LineId ;
+ int index ;
+
+} ;
+
+struct Indexed1D{
+ Double_t x ;
+ int index ;
+} ;
+
+struct Indexed2D{
+ Double_t x ;
+ int id ;
+ int index ;
+} ;
+
+vector<Double_t> diff(const vector<Double_t>& input);
+
+double distancePeak(const PeakData& p1, const PeakData& p2) ;
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/HistoFillerDEE.cpp b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/HistoFillerDEE.cpp
index 8a0f3e42010befef729749a14665f55399dfc105..18b87f6ee464a7daea23c0fbb718d5c421d9d659 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/HistoFillerDEE.cpp
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step4/HistoFillerDEE.cpp
@@ -17,7 +17,11 @@ void HistoFillerDEE(const char* Path){
TICPhysics *IC = new TICPhysics();
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
-
+
+ double FF_IC_X;
+ chain->SetBranchStatus("FF_IC_X", true);
+ chain->SetBranchAddress("FF_IC_X", &FF_IC_X);
+
TH2F *hDE_E = new TH2F("DE_E" , "DE_E", 1000,0,35000,1000,0,26000);
//===========================================================================================================
// Fill
@@ -37,17 +41,20 @@ void HistoFillerDEE(const char* Path){
} // end timer
chain->GetEntry(e);
- hDE_E->Fill(IC->Eres,IC->DE);
+
+ if(FF_IC_X>-530 ){
+ hDE_E->Fill(IC->Eres,IC->DE);
+ }
}
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/DE_E_" + number + ".root"; // Generate the filename
-
+
// Create the TFile
TFile* file = new TFile(fileName.c_str(), "RECREATE");
if (file->IsOpen()) {
@@ -55,7 +62,7 @@ void HistoFillerDEE(const char* Path){
} else {
std::cout << "Failed to create the file." << std::endl;
}
-
+
hDE_E->Write();
// Don't forget to close the file
file->Close();
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FillHistoCharge.C b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FillHistoCharge.C
index baadd87dfd8310164ed00290a03934c7f6ce9731..15c8c5f8baaabbe8d560ab00f12934abd849d526 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FillHistoCharge.C
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FillHistoCharge.C
@@ -15,7 +15,9 @@ void FillHistoCharge(){
chain->SetBranchAddress("IC", &IC);
-
+ double FF_IC_X;
+ chain->SetBranchStatus("FF_IC_X", true);
+ chain->SetBranchAddress("FF_IC_X", &FF_IC_X);
TH1F *hChargeAll = new TH1F("hcharge_all","hcharge_all",2000,20,80);
@@ -34,9 +36,11 @@ void FillHistoCharge(){
}
chain->GetEntry(e);
- hChargeAll->Fill(IC->Chio_Z);
- }
+ if(FF_IC_X>-530 ){
+ hChargeAll->Fill(IC->Chio_Z);
+ }
+ }
std::string strPath(Path);
std::regex pattern("Run(\\d{3})"); // Regex to match Run followed by exactly 3 digits
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FitCharge.C b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FitCharge.C
index 056531cf5841dacea24835756387e1fae80e3c3b..199123850da42e276926de7886dcdc0380f2b820 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FitCharge.C
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/FitCharge.C
@@ -26,10 +26,11 @@ void FitCharge(){
//////////////////////////////////////////////////////
void Fit(TH1F* hcharge,string Energy){
TGraphErrors* gsigma = new TGraphErrors();
+ TGraphErrors* gsigmaZ = new TGraphErrors();
hcharge->Rebin(rebin);
hcharge->Draw();
- int Zmin = 30;
+ int Zmin = 29;
int Zmax = 65;
@@ -93,6 +94,7 @@ void Fit(TH1F* hcharge,string Energy){
hcharge->Fit(total,"RBq");
Z = Zmin;
+ double MeanSigmaZ = 0 ;
//TF1* one_gauss = new TF1("one_gaus","gaus");
for(int i=0; i<NumberOfZ; i++){
@@ -123,10 +125,13 @@ void Fit(TH1F* hcharge,string Energy){
//cout << "Mean= " << mean << endl;
cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;
cout << "Sigma/Z = " << 2.35*sigma /Z*100. << endl;
+
+ if (2.35*sigma/Z*100 > 1 && 2.35*sigma/Z*100 < 3){MeanSigmaZ += 2.35*sigma/Z*100;}
gsigma->SetPoint(i,Z,sigma);
gsigma->SetPointError(i,0,sigma_err);
+ gsigmaZ->SetPoint(i,Z,2.35*sigma/Z*100);
Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
Integral_err[i] = (12.5/rebin)*sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
@@ -137,6 +142,8 @@ void Fit(TH1F* hcharge,string Energy){
Z++;
}
+ MeanSigmaZ = MeanSigmaZ/NumberOfZ;
+ cout << "********************" << endl << "Mean Sigma/Z " << MeanSigmaZ << endl;
double histo_integral = hcharge->Integral();
hcharge->Write();
total->Write();
@@ -164,4 +171,13 @@ void Fit(TH1F* hcharge,string Energy){
csig->cd();
gsigma->Draw("ap");
+ TCanvas* csigZ = new TCanvas("csigZ","csigZ",800,800);
+ gsigmaZ->SetMarkerStyle(8);
+ gsigmaZ->GetYaxis()->SetTitle("2.35*#sigma/Z*100");
+ gsigmaZ->GetXaxis()->SetTitle("Z");
+ csigZ->cd();
+ gsigmaZ->Draw("ap");
+ total->Draw("lsame");
+
+
}
diff --git a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/SplineChio.C b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/SplineChio.C
index 1781a9cdd04a5236011edf1822f854c69310aa61..f9dd91e5590972646bf06d8049d4874ffb145601 100644
--- a/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/SplineChio.C
+++ b/Projects/AlPhaPha/2024/macro/Calibration/IC/Step5/SplineChio.C
@@ -43,7 +43,7 @@ void SplineChio()
for(int i=0;i<Nrun;i++)
{
- chain->Add(Form("../../../../root/analysis/Run%d.root",RunNumber[i]));
+ chain->Add(Form("../../../../root/analysis/Run%dEcorr.root",RunNumber[i]));
cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
}
@@ -159,7 +159,7 @@ void MakeSpline(){
pfx[i]->SetDirectory(0);
// 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)+hProfile[i]->GetBinLowEdge(FirstBin+200));
+ TF1 * fitpol3 = new TF1(Form("FitPol3_pfx_%02d",i+1),"pol3",hProfile[i]->GetBinLowEdge(FirstBin),hProfile[i]->GetBinLowEdge(FirstBin+200));
hProfile[i]->Fit(fitpol3,"R");
fitpol3->GetParameters(parpol3);
@@ -572,7 +572,7 @@ void ApplySpline()
file >> p[i];
}
}
- int Nentries=1e6;
+ int Nentries=2e6;
//int Nentries = chain->GetEntries();
auto start = std::chrono::high_resolution_clock::now();
@@ -654,7 +654,7 @@ void ApplySpline()
//
// PeakFinder
TSpectrum *PeakFinder = new TSpectrum(100);
- int nPeaks = PeakFinder->Search(hChioDE_corr,3,"",0.043);
+ int nPeaks = PeakFinder->Search(hChioDE_corr,2,"",0.040);
Double_t * x = PeakFinder->GetPositionX();
std::sort(x,x+nPeaks);