diff --git a/Projects/AlPhaPha/2024/Analysis.cxx b/Projects/AlPhaPha/2024/Analysis.cxx
index a27c146286e260235b19e09ce647ea3654b3d146..e3120a941be4299a3701e405a48b859bc11c60f6 100644
--- a/Projects/AlPhaPha/2024/Analysis.cxx
+++ b/Projects/AlPhaPha/2024/Analysis.cxx
@@ -466,8 +466,32 @@ void Analysis::VamosAnalysis(){
UShort_t FPMWPat = MTOF_FP0_T0VN[0]; //20 section
FPMW_Section = FPMWPat;
+
+ //Toff
+ vector<double> Toff13 , Toff14, Toff23, Toff24;
+ const char* Path13 = "macro/mwpc/Toff/output/Toff13.txt";
+ const char* Path14 = "macro/mwpc/Toff/output/Toff14.txt";
+ const char* Path23 = "macro/mwpc/Toff/output/Toff23.txt";
+ const char* Path24 = "macro/mwpc/Toff/output/Toff24.txt";
+
+ Toff13 = TxtToVector(Path13);
+ Toff14 = TxtToVector(Path14);
+ Toff23 = TxtToVector(Path23);
+ Toff24 = TxtToVector(Path24);
+
+ for (int i = 0 ; i<Toff13.size() ; i++){
+ Time->SetToff_DT13(Toff13.at(i));
+ Time->SetToff_DT14(Toff14.at(i));
+ }
+ FF_IC_Y = FPMW->Yf + (1442.6+6774.4-7600)*tan(FPMW->Phif);
+ FF_IC_X = FPMW->Xf + (1442.6+6774.4-7600)*tan(FPMW->Thetaf);
+
IC->SetFPMWSection(FPMW_Section);
+ IC->SetTTimeData(Time);
+ IC->SetX(FF_IC_X);
+ IC->SetY(FF_IC_Y);
IC->BuildSimplePhysicalEvent();
+
double Theta = -1000;
@@ -484,29 +508,19 @@ void Analysis::VamosAnalysis(){
FF_Y3 = FPMW->PositionY[2];
//FF_IC_Y = FPMW->Yf + (1442.6+6774.4-7600)*tan(FPMW->Phif/1000);
- FF_IC_Y = FPMW->Yf + (1442.6+6774.4-7600)*tan(FPMW->Phif);
- FF_IC_X = FPMW->Xf + (1442.6+6774.4-7600)*tan(FPMW->Thetaf);
+ //FF_IC_Y = FPMW->Yf + (1442.6+6774.4-7600)*tan(FPMW->Phif);
+ //FF_IC_X = FPMW->Xf + (1442.6+6774.4-7600)*tan(FPMW->Thetaf);
// T13 //
double path1 = FPMW->GetDetectorPositionZ(0)/10./cos(FPMW->Theta_in)/cos(FPMW->Phi_in);
double path2 = (FPMW->GetDetectorPositionZ(2)-7600)/10./cos(FPMW->Thetaf);
- //Toff
- vector<double> Toff13 , Toff14, Toff23, Toff24;
- const char* Path13 = "macro/mwpc/Toff/output/Toff13.txt";
- const char* Path14 = "macro/mwpc/Toff/output/Toff14.txt";
- const char* Path23 = "macro/mwpc/Toff/output/Toff23.txt";
- const char* Path24 = "macro/mwpc/Toff/output/Toff24.txt";
- Toff13 = TxtToVector(Path13);
- Toff14 = TxtToVector(Path14);
- Toff23 = TxtToVector(Path23);
- Toff24 = TxtToVector(Path24);
double Toff[20] = {0, 588.0, 588.5, 587.95, 588.04, 587.72, 587.92, 587.9, 587.9, 588.66, 588.80, 588.67, 588.64, 588.75, 588.47, 588.65, 588.65, 588.67, 589.05, 590.3};
// To know time of electron drift in a section of the chio
for (int seg = 0 ; seg < IC->fIC_TS.size() ; seg++){
- FF_DriftTime.push_back(10* (IC->fIC_TS.at(seg) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))));
+ FF_DriftTime.push_back(10* (IC->fIC_TS.at(seg) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))));
}
diff --git a/Projects/AlPhaPha/2024/Snakefile b/Projects/AlPhaPha/2024/Snakefile
index 374a37903e8b95991d85c8577260d2d86afebc5d..6146edc8fefe5bc9fa4e482df0901725fe3ca33b 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_246"
+input_directory = os.getcwd() + "/../DataMacro/output/run_248"
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_246"
+phy_directory = os.getcwd() + "/../DataMacro/output/analysis/run_248"
#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_246/"+inputfile)
+ analysedfile.append( os.getcwd() + "/../DataMacro/output/analysis/run_248/"+inputfile)
## batch rules
rule all:
diff --git a/Projects/AlPhaPha/2024/configs/ConfigIC.dat b/Projects/AlPhaPha/2024/configs/ConfigIC.dat
index b81f9e5915c43cb423b6d3c3ca800d20f3446860..c8b2ded39f487f32dd86573a8596ddb5104d221a 100644
--- a/Projects/AlPhaPha/2024/configs/ConfigIC.dat
+++ b/Projects/AlPhaPha/2024/configs/ConfigIC.dat
@@ -1,3 +1,6 @@
ConfigIC
%LOAD_Z_SPLINE ./Calibration/VAMOS/CHIO/Z_spline.root
DATA_YEAR 2024
+ LOAD_DE_SPLINE ./Calibration/VAMOS/CHIO/Spline_DE.root
+ LOAD_Y_SPLINE ./Calibration/VAMOS/CHIO/Spline_Y.root
+ LOAD_XY0_PROFILE ./Calibration/VAMOS/CHIO/ICXY_Profile.root
diff --git a/Projects/AlPhaPha/2024/convert_snakemake_generic.sh b/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
index 526de6830dad5f03f9665d95bbb4ff429da0db62..354e11b5761d2323f353fca7cde853728273ec28 100755
--- a/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
+++ b/Projects/AlPhaPha/2024/convert_snakemake_generic.sh
@@ -1,6 +1,26 @@
#!/bin/bash
+
+#Declare run number
+RUN_NUMBER="248"
+
+# Check if Snakefile exists
+if [[ ! -f "Snakefile" ]]; then
+ echo "Error: Snakefile not found!"
+ exit 1
+fi
+
+# Replace all occurrences of run_<anything> with run_244
+sed -i "s/run_[a-zA-Z0-9]\{3\}/run_${RUN_NUMBER}/g" Snakefile
+echo "All occurrences of 'run_*' have been replaced with 'run_${RUN_NUMBER}' in Snakefile."
+
echo "- executing snakemake file for npanalysis..."
snakemake --cores 30 --forceall --keep-incomplete --keep-going --rerun-incomplete
+
echo "- snakemake executed successfully!"
echo "- Merging file..."
-root -q '../DataMacro/Merger.C(30,"root/analysis","Run246","../DataMacro/output/analysis/run_246/run_raw_246_")'
+
+OName="\"Run${RUN_NUMBER}\""
+OPATH="\"root/analysis\""
+Path="\"../DataMacro/output/analysis/run_${RUN_NUMBER}/run_raw_${RUN_NUMBER}_\""
+
+root -l -q "../DataMacro/Merger.C(30,${OPATH},${OName},${Path})"
diff --git a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/ReadMe.md b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/ReadMe.md
index 3e6f98ff227ed13a32aec8477ee7090fd282e46a..7486114a8d2b4ee6f0e4740343915ed62d052fd5 100644
--- a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/ReadMe.md
+++ b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/ReadMe.md
@@ -2,6 +2,7 @@
You need to fill Cut folder with the cut of each charge
-You can do this simply by using the DECorr.C macro in the YCalibration folder
+You can do this simply by using the DECorr.C macro in the XYCalibration folder
+You can also use the AutoCut Code from XYCalibration folder
-The name of the cut should be Z1->ZN and you should space them a little
+Be careful with cut name
diff --git a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
index eb7c1011f258abdd4000e661ba873e344470ab96..941c5e4be564c5f82947eae895a35a7f33b932d5 100644
--- a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
+++ b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
@@ -1,4 +1,6 @@
#include "TICPhysics.h"
+#include "TTimeData.h"
+#include <TFPMWPhysics.h>
#include <TChain.h>
#include <TCutG.h>
#include <TF1.h>
@@ -6,10 +8,10 @@
#include <TH2.h>
#include <TSpline.h>
#include <fstream>
-
+#include "../XYCalibration/ProfileEvaluator.h"
using namespace std;
-const int Ncuts=28;
+const int Ncuts=30;
const int NSegment=11;
const int Nrun=1;
int RunNumber[Nrun]={241};
@@ -18,6 +20,7 @@ void MakeSpline();
void ApplySpline_perTCutG();
void ApplySpline();
vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline);
+vector<double> TxtToVector(const char *Path);
///////////////////////////////////////////////////////////////////////////
void SplineChio()
@@ -25,41 +28,79 @@ void SplineChio()
//===========================================================================================================
// Loading Var
//===========================================================================================================
-
+
//**********************Cut**************************************
- TFile *fcut=new TFile("Cut/Cut_Z.root","open");
+ TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root","open");
TCutG *cutZ[Ncuts];
- for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Z%i",i+1));
+ for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
// ****************** Spline Y **********************************
- TFile *fSplineIC = new TFile("../YCalibration/Output/spline_P2P_2024.root","open");
+ TFile *fSplineDe;
+ TSpline3 *splineDE,*splineDE0, *splineDE02,*splineDE023,*splineDE0234,*splineDEbis;
+ fSplineDe= new TFile("../XYCalibration/Output/splineDE.root","open");
+ splineDE = (TSpline3*) fSplineDe->Get("SplineDe");
+ splineDE0 = (TSpline3*) fSplineDe->Get("SplineDe_0");
+ splineDE02 = (TSpline3*) fSplineDe->Get("SplineDe_02");
+ splineDE023 = (TSpline3*) fSplineDe->Get("SplineDe_023");
+ splineDE0234 = (TSpline3*) fSplineDe->Get("SplineDe_0234");
- vector<vector<TSpline3*>> Spline = LoadSpline("../YCalibration/Output/spline_P2P_2024.root");
+
+ TFile *fSplineIC = new TFile("../XYCalibration/Output/spline_P2P_2024.root","open");
+
+ vector<vector<TSpline3*>> Spline = LoadSpline("../XYCalibration/Output/spline_P2P_2024.root");
vector<TSpline3*> spline_X(11), spline_Y(11);
spline_X = Spline.at(0) ; spline_Y = Spline.at(1);
+
+
//**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree");
-
- chain->Add("../../../root/analysis/VamosCalib241.root");
-
+
+ chain->Add("../../../root/analysis/VamosCalib247.root");
+
TICPhysics* IC = new TICPhysics() ;
+ TTimeData *Time = new TTimeData() ;
+ TFPMWPhysics *FPMW = new TFPMWPhysics();
+
double FF_IC_X, FF_IC_Y, FF_V13;
double FF_DE, FF_Eres;
-
+
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
chain->SetBranchStatus("FF_IC_Y", true);
chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y);
+ // ****************** Load Time **********************************
+
+ chain->SetBranchStatus("Time", true);
+ chain->SetBranchAddress("Time", &Time);
+
+ vector<double> Toff13 , Toff14, Toff23, Toff24;
+ const char* Path13 = "../../mwpc/Toff/output/Toff13.txt";
+ const char* Path14 = "../../mwpc/Toff/output/Toff14.txt";
+ const char* Path23 = "../../mwpc/Toff/output/Toff23.txt";
+ const char* Path24 = "../../mwpc/Toff/output/Toff24.txt";
+
+ Toff13 = TxtToVector(Path13);
+ Toff14 = TxtToVector(Path14);
+ Toff23 = TxtToVector(Path23);
+ Toff24 = TxtToVector(Path24);
+
+ ProfileEvaluator Profile;
+ Profile.LoadProfile("../XYCalibration/Output/RatioProfile.root","ICOneZeroProfile");
+
+ //*****************************************************************
TH2F *hChioDE_E_all = new TH2F("hChioDE_E","hChioDE_E",1000,1000,20000,500,1000,26000);
TH2F *hChioDE_E[Ncuts];
for(int i=0;i<Ncuts;i++) hChioDE_E[i]=new TH2F(Form("hChioDE_E_%d",i+1),Form("hChioDE_E_%d",i+1),2000,1000,20000,1250,1000,26000);
auto start = std::chrono::high_resolution_clock::now();
- int Nentries=chain->GetEntries();
+
+ //int Nentries=chain->GetEntries();
+ int Nentries=1000000;
+
for(int e=0;e<Nentries;++e)
{
if (e % 100000 == 0 && e > 0 ) {
@@ -76,51 +117,76 @@ void SplineChio()
vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the
vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first
// segment of ic
- for (int seg = 1; seg < NSegment+1 ; seg++) {
- if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0
- if (spline_Y.at(seg)==0){
- ICcorr_Y.at(seg) = IC->fIC_raw[seg];
- }
- else {
- if (seg == 0) {
- Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
- ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg);
- }
+ if (Time->GetMWPC13Mult() ==1 && IC->fIC_TS.size()>=8){ //only mult 1 event
+
+ //*************************************** GET TIME
+ UShort_t FPMW_Section = Time->GetSection_MWPC3(0);
+ double FF_DriftTime = 10* (IC->fIC_TS.at(0) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))) ;
- else if (seg == 1){
- ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y);
- }
+ // ********************************* ALIGN IC***************************************************************
+ for (int seg = 1; seg < IC->fIC_TS.size() ; seg++) { // loop on multiplicity of event
+ if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0
- else if (seg > 1) {
- Temp_Y.at(seg) = IC->fIC_raw[seg]/ICcorr_Y.at(seg-1) * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
- ICcorr_Y.at(seg) = ICcorr_Y.at(seg-1) * Temp_Y.at(seg);
+ double FF_DriftTime_temp = 10* (IC->fIC_TS.at(seg) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))) ;
+ if (spline_Y.at(seg)==0){
+ ICcorr_Y.at(seg) = IC->fIC_raw[seg];
}
- if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
- } //end if non empty
- if (seg == 0) break;
- }//endloop seg
+ else {
+ if (seg == 0) {
+ ICcorr_Y.at(seg) = ICcorr_Y.at(0) * spline_Y.at(0)->Eval(0)/spline_Y.at(0)->Eval(FF_DriftTime);
+ }
+
+ else if (seg == 1){
+ ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_DriftTime);
+ }
+
+ else if (seg > 1) {
+ ICcorr_Y.at(seg) = IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_DriftTime);
+ }
+ if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
+ } //end if non empty
+ if (seg == 0) break;
+ }//endloop seg
+
+ //********************* CORR IC 0 ******************************
+ Double_t PolX = 1.37622;
+
+ Double_t ICRatio = IC->fIC_raw[1]/IC->fIC_raw[0];
+ Double_t ICRatioCorr = ICRatio * PolX / Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
+ Double_t ICcorr_Y0 = IC->fIC_raw[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
+
+ if ( ICRatioCorr<1.4 || ICRatioCorr >1.5){
+ ICRatioCorr = ICRatio * PolX / Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
+ ICcorr_Y0 = IC->fIC_raw[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
+ if (ICRatioCorr >100) {
+ ICcorr_Y0 = IC->fIC_raw[0];
+ }
+ }
+ // ************************** CORR Y GLOB *****************************
- FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
- FF_Eres = 0 ;
+ double FF_DE_precorr = 0.5*( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ ICcorr_Y.at(3)) + ICcorr_Y.at(4);
+ FF_DE = FF_DE_precorr * splineDE0234->Eval(0) / splineDE0234->Eval(FF_DriftTime);
+ FF_Eres = 0 ;
- for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
- if (seg >4){
- FF_Eres += IC->fIC_raw[seg] ;
- //Ecorr += ICcorr_Y.at(seg) ;
+ for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+ if (seg >4){
+ if (seg == 5) FF_Eres += double(IC->fIC_raw[seg]);
+ else FF_Eres += 2*double(IC->fIC_raw[seg]) ;
+ }
}
- }
- hChioDE_E_all->Fill(FF_Eres,FF_DE);
+ hChioDE_E_all->Fill(FF_Eres,FF_DE);
- //for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(Chio_E,FF_DE)) {hChioDE_E[i]->Fill(Chio_E,FF_DE); break;}
- for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) {hChioDE_E[i]->Fill(FF_Eres,FF_DE); break;}
- }
+ //for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(Chio_E,FF_DE)) {hChioDE_E[i]->Fill(Chio_E,FF_DE); break;}
+ for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) {hChioDE_E[i]->Fill(FF_Eres,FF_DE); break;}
+ } // end max mult 1
+ }//end loop event
//**********************************Out***************************
@@ -136,9 +202,9 @@ void SplineChio()
for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
}
+
///////////////////////////////////////////////////////////////////////////
-void MakeSpline()
-{
+void MakeSpline(){
TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root");
TSpline3 *gspline[Ncuts];
@@ -277,7 +343,7 @@ void ApplySpline_perTCutG()
{
// *************************Cut***************************
- TFile *fcut=new TFile("Cut/Cut_Z.root");
+ TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root");
TCutG *cutZ[Ncuts];
// ******************* TSpline DE*************************
@@ -285,27 +351,27 @@ void ApplySpline_perTCutG()
TSpline3 *gspline[Ncuts];
for(int i=0;i<Ncuts;i++){
- cutZ[i]=(TCutG*)fcut->Get(Form("Z%i",i+1));
+ cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
}
// ****************** Spline Y **********************************
- TFile *fSplineIC = new TFile("../YCalibration/Output/spline_P2P_2024.root","open");
+ TFile *fSplineIC = new TFile("../XYCalibration/Output/spline_P2P_2024.root","open");
- vector<vector<TSpline3*>> Spline = LoadSpline("../YCalibration/Output/spline_P2P_2024.root");
+ vector<vector<TSpline3*>> Spline = LoadSpline("../XYCalibration/Output/spline_P2P_2024.root");
vector<TSpline3*> spline_X(11), spline_Y(11);
spline_X = Spline.at(0) ; spline_Y = Spline.at(1);
//**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree");
-
+
chain->Add("../../../root/analysis/VamosCalib241.root");
-
+
TICPhysics* IC = new TICPhysics() ;
double FF_IC_X, FF_IC_Y, FF_V13;
double FF_DE, FF_Eres;
-
+
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
chain->SetBranchStatus("FF_IC_Y", true);
@@ -342,86 +408,86 @@ void ApplySpline_perTCutG()
19043.9,
19437.8,
19889.4,};
- // 20301.1,
- //20709.7,
- //21112.6,
- //21544,
- //21953.4,
- //22361.8};
+ // 20301.1,
+ //20709.7,
+ //21112.6,
+ //21544,
+ //21953.4,
+ //22361.8};
- TH2F *hChioDE_E_all = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000);
- TH2F *hChioDE_E = new TH2F("hChioDE_E","hChioDE_E",1000,1000,41000,500,1000,26000);
- TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);
+TH2F *hChioDE_E_all = new TH2F("hChioDE_E_all","hChioDE_E_all",1000,1000,41000,500,1000,26000);
+TH2F *hChioDE_E = new TH2F("hChioDE_E","hChioDE_E",1000,1000,41000,500,1000,26000);
+TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);
- int Nentries=chain->GetEntries();
- 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);
+int Nentries=chain->GetEntries();
+auto start = std::chrono::high_resolution_clock::now();
- std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
- }
+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);
- chain->GetEntry(e);
+ std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
+ }
+ chain->GetEntry(e);
- vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the
- vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first
- // segment of ic
- for (int seg = 1; seg < NSegment+1 ; seg++) {
- if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0
+ vector<double> ICcorr_Y(NSegment), ICcorr_X(NSegment); // the [0] element of spline is the
+ vector<double> Temp_X(NSegment) , Temp_Y(NSegment); // correction on the first
+ // segment of ic
+ for (int seg = 1; seg < NSegment+1 ; seg++) {
- if (spline_Y.at(seg)==0){
- ICcorr_Y.at(seg) = IC->fIC_raw[seg];
- }
+ if (seg == NSegment) seg = 0; //from 1to NSeg finishing with 0
- else {
- if (seg == 0) {
- Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
- ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg);
- }
+ if (spline_Y.at(seg)==0){
+ ICcorr_Y.at(seg) = IC->fIC_raw[seg];
+ }
- else if (seg == 1){
- ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y);
- }
+ else {
+ if (seg == 0) {
+ Temp_Y.at(seg) = ICcorr_Y.at(1) / IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(1) / Temp_Y.at(seg);
+ }
- else if (seg > 1) {
- Temp_Y.at(seg) = IC->fIC_raw[seg]/ICcorr_Y.at(seg-1) * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
- ICcorr_Y.at(seg) = ICcorr_Y.at(seg-1) * Temp_Y.at(seg);
- }
+ else if (seg == 1){
+ ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_IC_Y);
+ }
- if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
- } //end if non empty
- if (seg == 0) break;
- }//endloop seg
+ else if (seg > 1) {
+ Temp_Y.at(seg) = IC->fIC_raw[seg]/ICcorr_Y.at(seg-1) * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_IC_Y);
+ ICcorr_Y.at(seg) = ICcorr_Y.at(seg-1) * Temp_Y.at(seg);
+ }
+ if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
+ } //end if non empty
+ if (seg == 0) break;
+ }//endloop seg
- FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
- FF_Eres = 0 ;
- for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
- if (seg >4){
- FF_Eres += IC->fIC_raw[seg] ;
- //Ecorr += ICcorr_Y.at(seg) ;
- }
- }
+ FF_DE = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+ FF_Eres = 0 ;
- hChioDE_E_all->Fill(FF_Eres,FF_DE);
- for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) {
- hChioDE_E->Fill(FF_Eres,FF_DE);
- FF_DEcorr = FF_DEcorr0[i] * FF_DE / gspline[i]->Eval(FF_Eres);
- hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr);
- break;
+ for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+ if (seg >4){
+ FF_Eres += IC->fIC_raw[seg] ;
+ //Ecorr += ICcorr_Y.at(seg) ;
}
}
+ hChioDE_E_all->Fill(FF_Eres,FF_DE);
+ for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(FF_Eres,FF_DE)) {
+ hChioDE_E->Fill(FF_Eres,FF_DE);
+ FF_DEcorr = FF_DEcorr0[i] * FF_DE / gspline[i]->Eval(FF_Eres);
+ hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr);
+ break;
+ }
+}
+
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,1000);
@@ -458,13 +524,13 @@ void ApplySpline()
//**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree");
-
+
chain->Add("../../../root/analysis/VamosCalib241.root");
-
+
TICPhysics* IC = new TICPhysics() ;
double FF_IC_X, FF_IC_Y, FF_V13;
double FF_DE, FF_Eres;
-
+
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
chain->SetBranchStatus("FF_IC_Y", true);
@@ -774,3 +840,26 @@ vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline){
return Spline;
}
+
+
+
+vector<double> TxtToVector(const char *Path){
+ string line;
+ vector<double> values;
+ ifstream file(Path);
+
+ if (file.is_open()) {
+ while (std::getline(file, line)) {
+ try {
+ values.push_back(std::stod(line));
+ } catch (const std::invalid_argument& e) {
+ std::cerr << "Invalid number in line: " << line << '\n';
+ }
+ }
+ file.close();
+ } else {
+ std::cerr << "Error opening file.\n";
+ }
+
+ return values;
+}
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
index 4225a9fd4bbfc0e62eff7e64bb9614221384114f..dcc3c68990a3365895f68954008e20b065bc4ef1 100644
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
@@ -13,8 +13,8 @@
void CutAutoDEE(){
// Retrieve initial histogram
- TFile *inFile = new TFile("Output/DE_E.root","open");
- TH2F *hDE_E = (TH2F*)inFile->Get("DE0234_Ecorr");
+ TFile *inFile = new TFile("Output/DE_E_merged.root","open");
+ TH2F *hDE_E = (TH2F*)inFile->Get("DE_E");
// Inverse data to use tspectrum
for (int xBin=1; xBin<= hDE_E->GetNbinsX();xBin++){
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
index 461cd1e6158683dc4571d3bebfe75565f6ec6d51..cdf08fc2bfd5b0635987c8c67efb88e481a6af0a 100644
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
@@ -11,7 +11,7 @@
#include <TStyle.h>
#include <fstream>
#include <vector>
-#include "ProfileEvaluator.h"
+#include "TProfileEvaluator.h"
using namespace std;
vector<double> TxtToVector(const char *Path);
@@ -27,7 +27,7 @@ void DECorr(bool cut = false, bool spline = false) {
//===========================================================================================================
TChain* chain = new TChain("PhysicsTree");
chain->Add("../../../root/analysis/VamosCalib247.root");
- chain->Add("../../../root/analysis/Run246.root");
+ //chain->Add("../../../root/analysis/Run246.root");
TICPhysics* IC = new TICPhysics() ;
TTimeData *Time = new TTimeData() ;
@@ -105,27 +105,13 @@ void DECorr(bool cut = false, bool spline = false) {
SplineCount ++;
}
}
-
vector<TSpline3*> spline_X(11), spline_Y(11);
- for (int i = 0; i < SplineCount ; i++) {
- int seg = int((i-2)/2 +1 );
- if (i < 2){
- spline_X.at(0) = (TSpline3 *)fSplineIC->Get("fspline_IC0_X");
- spline_Y.at(0) = (TSpline3 *)fSplineIC->Get("fspline_IC0_Y");
- }
-
- else if ( i>=2 && i<4){
- spline_X.at(1) = (TSpline3 *)fSplineIC->Get("fspline_IC1_X");
- spline_Y.at(1) = (TSpline3 *)fSplineIC->Get("fspline_IC1_Y");
- }
- else if (seg >=2 && (i%2 == 0) ) {
- spline_X.at(seg) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_X",seg));
- }
- else if (seg >=2 && !(i%2 == 0) ) {
- spline_Y.at(seg) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_Y",seg));
- }
+ for (int i = 0; i < 11 ; i++) {
+ spline_X.at(i) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_X",i));
+ spline_Y.at(i) = (TSpline3 *)fSplineIC->Get(Form("fspline_IC%d_Y",i));
} //End loop on histogram
-
+ //
+
double XmaxDE = 26000, XMaxE = 35000;
//Defining output histo
TH2F* hDE_E = new TH2F("DE_E","DE_E",1000,0,XMaxE,1000,0,XmaxDE);
@@ -188,9 +174,10 @@ void DECorr(bool cut = false, bool spline = false) {
//===========================================================================================================
int Nentries = chain->GetEntries();
//int Nentries = 1000000;
+ //int Nentries = 22000;
auto start = std::chrono::high_resolution_clock::now();
int NSegment = 11;
-
+ ofstream debug("Output/Debug.txt");
//save data for clustering
for (int e = 0; e < Nentries; e++) {
@@ -209,7 +196,7 @@ void DECorr(bool cut = false, bool spline = false) {
- if (Time->GetMWPC13Mult() ==1 && IC->fIC_TS.size()>=8){ //only mult 1 event
+ if (Time->GetMWPC13Mult() ==1 && IC->fIC_TS.size()>=8 ){ //only mult 1 event
UShort_t FPMW_Section = Time->GetSection_MWPC3(0);
double FF_DriftTime = 10* (IC->fIC_TS.at(0) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))) ;
@@ -218,7 +205,7 @@ void DECorr(bool cut = false, bool spline = false) {
double FF_DriftTime_temp = 10* (IC->fIC_TS.at(seg) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))) ;
if (spline_Y.at(seg)==0){
- ICcorr_Y.at(seg) = IC->fIC_raw[seg];
+ ICcorr_Y.at(seg) = IC->fIC_PID[seg];
}
else {
@@ -227,11 +214,11 @@ void DECorr(bool cut = false, bool spline = false) {
}
else if (seg == 1){
- ICcorr_Y.at(seg) = IC->fIC_raw[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_DriftTime);
+ ICcorr_Y.at(seg) = IC->fIC_PID[1]*spline_Y.at(1)->Eval(0)/spline_Y.at(1)->Eval(FF_DriftTime);
}
else if (seg > 1) {
- ICcorr_Y.at(seg) = IC->fIC_raw[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_DriftTime);
+ ICcorr_Y.at(seg) = IC->fIC_PID[seg] * spline_Y.at(seg)->Eval(0)/ spline_Y.at(seg)->Eval(FF_DriftTime);
}
if (!(ICcorr_Y.at(seg)==ICcorr_Y.at(seg))) ICcorr_Y.at(seg) = 0;
} //end if non empty
@@ -240,45 +227,45 @@ void DECorr(bool cut = false, bool spline = false) {
- Double_t PolX = -3.10262 * 1e-7 *FF_IC_X*FF_IC_X - 7.02108 * 1e-5 * FF_IC_X + 1.37622;
- //Double_t PolX = 1.37622;
+ //Double_t PolX = -3.10262 * 1e-7 *FF_IC_X*FF_IC_X - 7.02108 * 1e-5 * FF_IC_X + 1.37622;
+ Double_t PolX = 1.37622;
- Double_t ICRatio = IC->fIC_raw[1]/IC->fIC_raw[0];
+ Double_t ICRatio = IC->fIC_PID[1]/IC->fIC_PID[0];
Double_t ICRatioCorr = ICRatio * PolX / Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
- Double_t ICcorr_Y0 = IC->fIC_raw[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
+ Double_t ICcorr_Y0 = IC->fIC_PID[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
if ( ICRatioCorr<1.4 || ICRatioCorr >1.5){
ICRatioCorr = ICRatio * PolX / Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
- ICcorr_Y0 = IC->fIC_raw[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
+ ICcorr_Y0 = IC->fIC_PID[0] / PolX * Profile.Evaluate(FF_IC_X,FF_DriftTime,false);
if (ICRatioCorr >100) {
- ICcorr_Y0 = IC->fIC_raw[0];
+ ICcorr_Y0 = IC->fIC_PID[0];
}
}
double DE = 0 , DE_splined=0 ,DE_corr=0, DE_Ybis=0;
double E =0, E_splined=0;
- for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+ 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_raw[seg]);
- else E += 2*double(IC->fIC_raw[seg]) ;
+ if (seg == 5) E += double(IC->fIC_PID[seg]);
+ else E +=double(IC->fIC_PID[seg]) ;
}
}
// ********************* DE setter *********************
- double IC1corr = (IC->fIC_raw[1]*(1-0.000686068*FF_IC_Y))*(1-4.88238e-05*FF_IC_Y+7.40395e-06*FF_IC_Y*FF_IC_Y);
- double DE_Bis = 0.5*(IC1corr+IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+ double IC1corr = (IC->fIC_PID[1]*(1-0.000686068*FF_IC_Y))*(1-4.88238e-05*FF_IC_Y+7.40395e-06*FF_IC_Y*FF_IC_Y);
+ double DE_Bis = (IC1corr+IC->fIC_PID[2]+IC->fIC_PID[3])+IC->fIC_PID[4];
- DE = 0.5* ( IC->fIC_raw[0]+ double(IC->fIC_raw[1]) +double(IC->fIC_raw[2])+ double(IC->fIC_raw[3])) + double(IC->fIC_raw[4]);
- DE_splined = 0.5* ( ICcorr_Y0+ ICcorr_Y.at(1) +IC->fIC_raw[2]+ IC->fIC_raw[3]) + IC->fIC_raw[4];
- double DE0 = 0.5* ( ICcorr_Y0 + ICcorr_Y.at(1) +IC->fIC_raw[2] + IC->fIC_raw[3]) + IC->fIC_raw[4];
- double DE02 = 0.5* ( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ IC->fIC_raw[3]) + IC->fIC_raw[4];
- double DE023 = 0.5* ( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ ICcorr_Y.at(3)) + IC->fIC_raw[4];
- double DE0234 = 0.5*( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ ICcorr_Y.at(3)) + ICcorr_Y.at(4);
- double DE1234 = 0.5*( double(IC->fIC_raw[1]) +double(IC->fIC_raw[2])+ double(IC->fIC_raw[3])) + double(IC->fIC_raw[4]);
- //DE_splined = 0.5*( ICcorr_Y[1] +IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+ DE = ( IC->fIC_PID[0]+ double(IC->fIC_PID[1]) +double(IC->fIC_PID[2])+ double(IC->fIC_PID[3])) + double(IC->fIC_PID[4]);
+ DE_splined = ( ICcorr_Y0+ ICcorr_Y.at(1) +IC->fIC_PID[2]+ IC->fIC_PID[3]) + IC->fIC_PID[4];
+ double DE0 = ( ICcorr_Y0 + ICcorr_Y.at(1) +IC->fIC_PID[2] + IC->fIC_PID[3]) + IC->fIC_PID[4];
+ double DE02 = ( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ IC->fIC_PID[3]) + IC->fIC_PID[4];
+ double DE023 = ( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ ICcorr_Y.at(3)) + IC->fIC_PID[4];
+ double DE0234 = ( ICcorr_Y0+ ICcorr_Y.at(1) +ICcorr_Y.at(2)+ ICcorr_Y.at(3)) + ICcorr_Y.at(4);
+ double DE1234 = ( double(IC->fIC_PID[1]) +double(IC->fIC_PID[2])+ double(IC->fIC_PID[3])) + double(IC->fIC_PID[4]);
+ //DE_splined = 0.5*( ICcorr_Y[1] +IC->fIC_PID[2]+IC->fIC_PID[3])+IC->fIC_PID[4];
//************* Cut ******************
@@ -290,7 +277,7 @@ void DECorr(bool cut = false, bool spline = false) {
//***********Fill histo*************
- if (FF_DriftTime >Ymin && FF_DriftTime <Ymax && (FF_Theta >0)){
+ if (/*FF_DriftTime >Ymin && FF_DriftTime <Ymax && (FF_Theta >0)*/ true){
//*************Init DE_E*****************
hDE_E -> Fill(E,DE);
@@ -312,6 +299,13 @@ void DECorr(bool cut = false, bool spline = false) {
double DE023_corr = DE023 * splineDE023->Eval(0) / splineDE023->Eval(FF_DriftTime);
double DE0234_corr = DE0234 * splineDE0234->Eval(0) / splineDE0234->Eval(FF_DriftTime);
+ for(int i =0 ; i<5 ; i++){
+ if (i==0) debug << ICcorr_Y0 << " ";
+ else debug << ICcorr_Y.at(i) << " " ;
+ }
+ debug << IC->fIC_PID[4];
+ debug << endl;
+
hDE_E_corr->Fill(E,DE_corr);
hDE0_E_corr->Fill(E,DE0_corr);
hDE02_E_corr->Fill(E,DE02_corr);
@@ -355,7 +349,7 @@ void DECorr(bool cut = false, bool spline = false) {
hDE_Y_splined->Fill(FF_DriftTime,DE_splined);
hDE_Y_corr->Fill(FF_DriftTime,DE_corr);
- hIC_Y->Fill(FF_DriftTime,IC->fIC_raw[0]);
+ hIC_Y->Fill(FF_DriftTime,IC->fIC_PID[0]);
hIC_Y_splined->Fill(FF_DriftTime,ICcorr_Y0);
@@ -380,7 +374,7 @@ void DECorr(bool cut = false, bool spline = false) {
hDE_Y_splined->Fill(FF_DriftTime,DE_splined);
hDE_Y_corr->Fill(FF_DriftTime,DE_corr);
- hIC_Y->Fill(FF_DriftTime,IC->fIC_raw[0]);
+ hIC_Y->Fill(FF_DriftTime,IC->fIC_PID[0]);
hIC_Y_splined->Fill(FF_DriftTime,ICcorr_Y0);
if (CutV ) {
@@ -431,11 +425,16 @@ void DECorr(bool cut = false, bool spline = false) {
splineDe0234->SetName("SplineDe_0234");
splineDe0234->Write();
-
-
osplineDEbis->SetName("SplineDebis");
osplineDEbis->Write();
fSplineDE->Close();
+
+ //output calibration
+ TFile *fSplineDEY = new TFile("../../../Calibration/VAMOS/CHIO/Spline_DE.root","recreate");
+ splineDe0234->SetName("fspline_1");
+ splineDe0234->Write();
+ fSplineDEY->Close();
+
//===========================================================================================================
// Drawing histo
//===========================================================================================================
@@ -443,54 +442,54 @@ void DECorr(bool cut = false, bool spline = false) {
gStyle->SetPalette(kRainBow);
// =========== Poster plot=============================
- TCanvas* c8 = new TCanvas("DE Before correction","DE Before correction",4000,4000);
-
- hDE_E->GetXaxis()->SetTitle("E (A.U)");
- hDE_E->GetXaxis()->SetRangeUser(3000,35000);
- hDE_E->GetXaxis()->CenterTitle();
-
- hDE_E->GetYaxis()->SetTitle("#Delta E (A.U)");
- hDE_E->GetYaxis()->SetRangeUser(6000,26000);
- hDE_E->GetYaxis()->CenterTitle();
-
- hDE_E->GetZaxis()->SetRangeUser(0,140);
- hDE_E->GetZaxis()->SetTitle("Number of events");
- hDE_E->GetZaxis()->CenterTitle();
-
- hDE_E->SetTitle("PID VAMOS before correction");
- hDE_E->SetStats(0);
- hDE_E->Draw("COLZ");
-
- c8->SetCanvasSize(4000,3000);
- c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.png");
- c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.pdf");
- c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.jpg");
-
- TCanvas* c9 = new TCanvas("DE after correction","DE after correction",1000,1000);
-
- hDE0234_E_corr->GetXaxis()->SetTitle("E (A.U)");
- hDE0234_E_corr->GetXaxis()->SetRangeUser(3000,35000);
- hDE0234_E_corr->GetXaxis()->CenterTitle();
-
- hDE0234_E_corr->GetYaxis()->SetTitle("#Delta E (A.U)");
- hDE0234_E_corr->GetYaxis()->SetRangeUser(6000,26000);
- hDE0234_E_corr->GetYaxis()->CenterTitle();
-
- hDE0234_E_corr->GetZaxis()->SetRangeUser(0,140);
- hDE0234_E_corr->GetZaxis()->SetTitle("Number of events");
- hDE0234_E_corr->GetZaxis()->CenterTitle();
-
- hDE0234_E_corr->SetTitle("PID VAMOS after correction");
- hDE0234_E_corr->SetStats(0);
- hDE0234_E_corr->Draw("COLZ");
-
- c9->SetCanvasSize(4000,3000);
- c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.png");
- c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.pdf");
- c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.jpg");
+ /*
+ TCanvas* c8 = new TCanvas("DE Before correction","DE Before correction",4000,4000);
+ hDE_E->GetXaxis()->SetTitle("E (A.U)");
+ hDE_E->GetXaxis()->SetRangeUser(3000,35000);
+ hDE_E->GetXaxis()->CenterTitle();
+
+ hDE_E->GetYaxis()->SetTitle("#Delta E (A.U)");
+ hDE_E->GetYaxis()->SetRangeUser(6000,26000);
+ hDE_E->GetYaxis()->CenterTitle();
+
+ hDE_E->GetZaxis()->SetRangeUser(0,140);
+ hDE_E->GetZaxis()->SetTitle("Number of events");
+ hDE_E->GetZaxis()->CenterTitle();
+
+ hDE_E->SetTitle("PID VAMOS before correction");
+ hDE_E->SetStats(0);
+ hDE_E->Draw("COLZ");
+
+ c8->SetCanvasSize(4000,3000);
+ c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.png");
+ c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.pdf");
+ c8->SaveAs("~/These/Plot/Poster/VamosPIDBeforeCorr.jpg");
+
+ TCanvas* c9 = new TCanvas("DE after correction","DE after correction",1000,1000);
+
+ hDE0234_E_corr->GetXaxis()->SetTitle("E (A.U)");
+ hDE0234_E_corr->GetXaxis()->SetRangeUser(3000,35000);
+ hDE0234_E_corr->GetXaxis()->CenterTitle();
+
+ hDE0234_E_corr->GetYaxis()->SetTitle("#Delta E (A.U)");
+ hDE0234_E_corr->GetYaxis()->SetRangeUser(6000,26000);
+ hDE0234_E_corr->GetYaxis()->CenterTitle();
+
+ hDE0234_E_corr->GetZaxis()->SetRangeUser(0,140);
+ hDE0234_E_corr->GetZaxis()->SetTitle("Number of events");
+ hDE0234_E_corr->GetZaxis()->CenterTitle();
+
+ hDE0234_E_corr->SetTitle("PID VAMOS after correction");
+ hDE0234_E_corr->SetStats(0);
+ hDE0234_E_corr->Draw("COLZ");
+
+ c9->SetCanvasSize(4000,3000);
+ c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.png");
+ c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.pdf");
+ c9->SaveAs("~/These/Plot/Poster/VamosPIDAfterCorr.jpg");
//gStyle->SetOptStat(1);
// =========== Poster plot=============================
-
+ */
TCanvas* c7 = new TCanvas("correction succesive","correction succesive",1500,1000);
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/ProfileEvaluator.h b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/ProfileEvaluator.h
deleted file mode 100644
index 0994ba971a8eb0dc2bb445420355b322baf2af80..0000000000000000000000000000000000000000
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/ProfileEvaluator.h
+++ /dev/null
@@ -1,127 +0,0 @@
-#ifndef PROFILE_EVAL
-#define PROFILE_EVAL
-
-#include <TFile.h>
-#include <TProfile2D.h>
-#include <iostream>
-using namespace std;
-
-class ProfileEvaluator {
- private:
- TFile* fFile;
- TProfile2D* fProfile;
-
- public:
-
- //***********************************************************//
-
- ProfileEvaluator() : fFile(nullptr), fProfile(nullptr) {}
-
- //***********************************************************//
-
- bool LoadProfile(const char* filename, const char* profname) {
- fFile = TFile::Open(filename, "READ");
- if (!fFile || fFile->IsZombie()) return false;
-
- fProfile = (TProfile2D*)fFile->Get(profname)->Clone("Profile");
- fProfile->SetDirectory(0);
- return fProfile != nullptr;
- }
-
- //***********************************************************//
-
- Double_t Evaluate(Double_t x, Double_t y, Bool_t useInterpolation = true) {
- if (!fProfile) return 0.0;
-
- Int_t binx = fProfile->GetXaxis()->FindBin(x);
- Int_t biny = fProfile->GetYaxis()->FindBin(y);
-
- if (!useInterpolation) return fProfile->GetBinContent(binx, biny);
-
- //Def Bin dist
- Double_t BinDistX = fProfile->GetXaxis()->GetBinCenter(2) - fProfile->GetXaxis()->GetBinCenter(1);
- Double_t BinDistY = fProfile->GetYaxis()->GetBinCenter(2) - fProfile->GetYaxis()->GetBinCenter(1);
-
- // Normalisation var
- Double_t norma = 0;
-
- // Fetch Bin Center for surrounding 3 by 3 bin and value
- vector<Double_t> Cx(3) , Cy(3) ;
- vector<vector<Double_t>> Q(3, vector<Double_t>(3)), W(3,vector<Double_t>(3));
-
- for (int i = -1 ; i<=1 ; i++){
- // Loading nearby bins
- int TempX = binx +i ;
- int PosX = i + 1 ;
-
- // Edge cases handling
- if (TempX < 1) TempX = 1 ;
- if (TempX > fProfile->GetNbinsX()) TempX = fProfile->GetNbinsX() ;
-
- // Loading center
- Cx.at(PosX) = fProfile->GetXaxis()->GetBinCenter(TempX);
- //Normat that we use POSX here
-
- // Loading value and weight
- for (int j = -1 ; j<=1 ; j++){
- // Need to translate to acess the vector
- int PosY = j + 1 ;
- int TempY = biny + j ;
-
- if (TempY < 1) TempY = 1 ;
- if (TempY > fProfile->GetNbinsY()) TempY = fProfile->GetNbinsY() ;
-
- Cy.at(PosY) = fProfile->GetYaxis()->GetBinCenter(TempY);
-
- // Loading value
- Q.at(PosX).at(PosY) = fProfile->GetBinContent(TempX,TempY);
-
- // Calculating weight
- Double_t Wx = 0 , Wy = 0 ;
- Wx = 1 - ( abs( x - Cx.at(PosX)) / (BinDistX) );
- Wy = 1 - ( abs( y - Cy.at(PosY)) / (BinDistY) );
-
- //handling max range
- if(Wx < 0 || Wy < 0) W.at(PosX).at(PosY) = 0;
- else W.at(PosX).at(PosY) = Wx * Wy;
- // Verifiy that sum of the weight is 1
- norma += W.at(PosX).at(PosY) ;
- }
- }
-
- // now just calculate the return value
-
- Double_t InterPol = 0 ;
-
- for (int PosX = 0 ; PosX <= 2 ; PosX ++ ){
- for (int PosY = 0 ; PosY <= 2 ; PosY ++ ){
- InterPol += Q.at(PosX).at(PosY) * W.at(PosX).at(PosY);
- }
- }
- return InterPol;
- }
-
- //***********************************************************//
-
- void PrintInfo() {
- if (!fProfile) return;
- printf("Profile: %s\n", fProfile->GetName());
- printf("X bins: %d, range: [%.2f, %.2f]\n",
- fProfile->GetNbinsX(),
- fProfile->GetXaxis()->GetXmin(),
- fProfile->GetXaxis()->GetXmax());
- printf("Y bins: %d, range: [%.2f, %.2f]\n",
- fProfile->GetNbinsY(),
- fProfile->GetYaxis()->GetXmin(),
- fProfile->GetYaxis()->GetXmax());
- }
-
- //***********************************************************//
-
- TProfile2D* GetProfile() { return fProfile; }
-
- //***********************************************************//
- ~ProfileEvaluator() {
- }
-};
-#endif
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/SplineChioRobust.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/SplineChioRobust.C
index 331abe5bdb4df46ef97850f2474af9504197dbbb..07bd7b2a34fdeaa44fcbf56908e08ef277beca8c 100644
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/SplineChioRobust.C
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/SplineChioRobust.C
@@ -33,8 +33,8 @@ void SplineChioRobust() {
double XMin[4] , XMax[4], YMin[4] , YMax[4];
XMin[0] = -520 ; XMax[0] = 380 ; YMin[0] = 2000 ; YMax[0] = 9500 ;
XMin[1] = -1000 ; XMax[1] = 1000 ; YMin[1] = 1000 ; YMax[1] = 10000 ;
- XMin[2] = -520 ; XMax[2] = 380 ; YMin[2] = 0 ; YMax[2] = 1.3 ;
- XMin[3] = -1000 ; XMax[3] = 1000 ; YMin[3] = 0 ; YMax[3] = 1.3 ;
+ XMin[2] = -520 ; XMax[2] = 380 ; YMin[2] = 0 ; YMax[2] = 3 ;
+ XMin[3] = -1000 ; XMax[3] = 1000 ; YMin[3] = 0 ; YMax[3] = 3 ;
//Reinitialising the outfile
TFile *OutFile = new TFile("Output/HistoP2P.root","recreate");
@@ -67,6 +67,29 @@ void SplineChioRobust() {
if(seg ==1) {
}
+ else if (seg == 4){
+
+ //Load Histo
+ hIC.at(seg) = (HistoFillerIC(seg,SplineAllIC));
+
+ //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.at(0)->SetName(Form("fspline_IC%d_X",seg));
+ NCallSpline+=1;
+
+ // Y Spline
+ SplineICurrent.push_back(MakeSpline(hIC[seg][1] , NCallSpline, XMin[3] , XMax[3],0,2.5));
+ SplineICurrent.at(1)->SetName(Form("fspline_IC%d_Y",seg));
+ NCallSpline+=1;
+
+ //push vector into memory
+ SplineAllIC.at(seg)= SplineICurrent;
+
+ }
+
else{
//Load Histo
@@ -99,6 +122,15 @@ void SplineChioRobust() {
}
fspline->Close();
+ //Write for analysis
+ TFile* ospline = new TFile("../../../Calibration/VAMOS/CHIO/Spline_Y.root", "recreate");
+ for (int seg = 2 ; seg < NSegment ; seg++){
+ if (!SplineAllIC.at(seg).empty()){
+ SplineAllIC.at(seg)[1]->SetName(Form("fspline_%d",seg-1));
+ SplineAllIC.at(seg)[1]->Write();
+ }
+ }
+ fspline->Close();
vector<vector<TH2F*>> hIC_Corrall(NSegment);
vector<TH1F*> hIC_Corr_profile(NSegment);
@@ -179,11 +211,19 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
}
+ else if (segment==4) {
+ 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, 2, 2.5);
+ 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, 2, 2.5);
+ }
+
+
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, 2);
+ Form("hChio_IC%d_IC%d_Y",segment,segment-1), 1000, -1000, 1000, 500, 0, 4);
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, 2);
+ Form("hChio_IC%d_IC%d_X",segment,segment-1), 1000, -800, 800, 500, 0, 4);
}
@@ -191,8 +231,8 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
// Beginning loop on entries
//===========================================================================================================
- int Nentries = chain->GetEntries();
- //int Nentries = 1000000;
+ //int Nentries = chain->GetEntries();
+ int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();
for (int e = 0; e < Nentries; e++) {
@@ -200,7 +240,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
UShort_t FPMW_Section = Time->GetSection_MWPC3(0);
@@ -211,14 +251,14 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
for (int pseg = 1 ; pseg<=NSpline ; pseg++){
//If the spline doesnt exist define default value
if(Spline_All_ICprev.at(pseg).empty()){
- IC_Spline_CorrX.at(pseg)= IC->fIC_raw[pseg];
- IC_Spline_CorrY.at(pseg)= IC->fIC_raw[pseg];
+ IC_Spline_CorrX.at(pseg)= IC->fIC_PID[pseg];
+ IC_Spline_CorrY.at(pseg)= IC->fIC_PID[pseg];
}
//If it exist correct the current line
else if (!Spline_All_ICprev.at(pseg).empty()){
- IC_Spline_CorrX.at(pseg) = IC->fIC_raw[pseg] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
- IC_Spline_CorrY.at(pseg) = IC->fIC_raw[pseg] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
+ IC_Spline_CorrX.at(pseg) = IC->fIC_PID[pseg] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
+ IC_Spline_CorrY.at(pseg) = IC->fIC_PID[pseg] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
} //End if non empty
} // End loop IC non 0
@@ -230,19 +270,18 @@ vector<TH2F*> HistoFillerIC(int segment, vector<vector<TSpline3*>> Spline_All_IC
//===========================================================================================================/
if (segment == 1){
- hIC_Y->Fill(FF_IC_Y,IC->fIC_raw[segment]);
- hIC_X->Fill(FF_IC_X,IC->fIC_raw[segment]);
+ hIC_Y->Fill(FF_IC_Y,IC->fIC_PID[segment]);
+ hIC_X->Fill(FF_IC_X,IC->fIC_PID[segment]);
}
else if (segment == 0){
- hIC_Y->Fill(FF_IC_Y,IC->fIC_raw[segment]/IC_Spline_CorrY.at(1));
- hIC_X->Fill(FF_IC_X,IC->fIC_raw[segment]/IC_Spline_CorrX.at(1));
+ hIC_Y->Fill(FF_IC_Y,IC->fIC_PID[segment]/IC_Spline_CorrY.at(1));
+ hIC_X->Fill(FF_IC_X,IC->fIC_PID[segment]/IC_Spline_CorrX.at(1));
}
-
else {
- hIC_Y->Fill(FF_IC_Y,IC->fIC_raw[segment]/IC_Spline_CorrY.at(segment-1));
- hIC_X->Fill(FF_IC_X,IC->fIC_raw[segment]/IC_Spline_CorrX.at(segment-1));
+ hIC_Y->Fill(FF_IC_Y,IC->fIC_PID[segment]/IC_Spline_CorrY.at(segment-1));
+ hIC_X->Fill(FF_IC_X,IC->fIC_PID[segment]/IC_Spline_CorrX.at(segment-1));
}
} //end only m2
@@ -332,9 +371,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, 2);
+ Form("hChiocorr_IC%d_IC%d_Y",pseg,pseg-1), 1000, -1000, 1000, 500, 0, 4);
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, 2);
+ Form("hChiocorr_IC%d_IC%d_X",pseg,pseg-1), 1000, -800, 800, 500, 0, 4);
}
}
@@ -342,8 +381,8 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
// Beginning loop on entries
//===========================================================================================================
- int Nentries = chain->GetEntries();
- //int Nentries = 1000000;
+ //int Nentries = chain->GetEntries();
+ int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();
for (int e = 0; e < Nentries; e++) {
@@ -361,15 +400,15 @@ vector<vector<TH2F*>>HistoFillerICcorr(int segment, vector<vector<TSpline3*>> Sp
for (int pseg = 0 ; pseg<=NSpline ; pseg++){
//If the spline doesnt exist define default value
if(Spline_All_ICprev.at(pseg).empty()){
- IC_Spline_CorrX.at(pseg)= IC->fIC_raw[pseg];
- IC_Spline_CorrY.at(pseg)= IC->fIC_raw[pseg];
+ IC_Spline_CorrX.at(pseg)= IC->fIC_PID[pseg];
+ IC_Spline_CorrY.at(pseg)= IC->fIC_PID[pseg];
//cout << pseg << endl;
}
//If it exist correct the current line
else if (!Spline_All_ICprev.at(pseg).empty()){
- IC_Spline_CorrX.at(pseg) = IC->fIC_raw[pseg] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
- IC_Spline_CorrY.at(pseg) = IC->fIC_raw[pseg] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
+ IC_Spline_CorrX.at(pseg) = IC->fIC_PID[pseg] * SplineIC_X.at(pseg)->Eval(0) / SplineIC_X.at(pseg)->Eval(FF_IC_X);
+ IC_Spline_CorrY.at(pseg) = IC->fIC_PID[pseg] * SplineIC_Y.at(pseg)->Eval(0) / SplineIC_Y.at(pseg)->Eval(FF_IC_Y);
} //End if non empty
} // end loop
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/VerifCorrelation.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/VerifCorrelation.C
index ba867ff95f7dae4cab22dacff3760196281bd812..0c9a87d410524d5b2658315059e63c9ea72b4fe2 100644
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/VerifCorrelation.C
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/VerifCorrelation.C
@@ -49,7 +49,7 @@ void VerifCorrelation(){
TH3F *ICXY = new TH3F("ICXY","ICXY",100,-550,400,100,-1000,1000,100,0,2);
TH2F *ICX = new TH2F("ICX","ICX",1000,-500,400,1000,0,2);
- TProfile2D *Pxyz = new TProfile2D("ICOneZeroProfile","ICOneZeroProfile",100,-550.0,400.0,100,-900.0,500.0);
+ TProfile2D *Pxyz = new TProfile2D("ICOneZeroProfile","ICOneZeroProfile",100,-550.0,400.0,100,-800.0,480.0);
int Nentries = chain->GetEntries();
//int Nentries = 1000000;
@@ -68,14 +68,14 @@ void VerifCorrelation(){
chain->GetEntry(e);
- if(Time->GetMWPC13Mult()==1 && IC->fIC_raw[8] >0){
+ if(Time->GetMWPC13Mult()==1 && IC->fIC_PID[8] >0){
UShort_t FPMW_Section = Time->GetSection_MWPC3(0);
FF_DriftTime = 10* (IC->fIC_TS.at(0) - Time->GetTS_MWPC13(0)) - ((Time->GetTime_MWPC13(0)+Toff13.at(FPMW_Section))) ;
- Double_t ICRatio = IC->fIC_raw[1]/IC->fIC_raw[0];
+ Double_t ICRatio = IC->fIC_PID[1]/IC->fIC_PID[0];
- if(IC->fIC_raw[8] >0 && ICRatio <= 3 && ICRatio >0.1 ) {
+ if(IC->fIC_PID[8] >0 && ICRatio <= 3 && ICRatio >0.1 ) {
ICXY->Fill(FF_IC_X,FF_DriftTime,ICRatio);
Pxyz->Fill(FF_IC_X,FF_DriftTime,ICRatio);
@@ -134,6 +134,10 @@ void VerifCorrelation(){
TFile *ofile = new TFile("Output/RatioProfile.root","recreate");
Pxyz->Write();
+
+
+ TFile *oconf = new TFile("../../../Calibration/VAMOS/CHIO/ICXY_Profile.root","recreate");
+ Pxyz->Write();
}
vector<double> TxtToVector(const char *Path){