diff --git a/NPLib/Detectors/Epic/TEpicData.cxx b/NPLib/Detectors/Epic/TEpicData.cxx
index 955f35ec264c0f785dbe84c534f1d358ac7b2366..d754eb4fa5af4b0c7847458b1598a56bc49fded4 100644
--- a/NPLib/Detectors/Epic/TEpicData.cxx
+++ b/NPLib/Detectors/Epic/TEpicData.cxx
@@ -50,8 +50,8 @@ void TEpicData::Dump() const {
cout << "MultAnodes = : " << mysize << endl;
for (size_t i = 0 ; i < mysize ; i++){
- cout << "AnodeNbr: " << GetAnodeNbr(i)
- << " Q1: " << GetQ1(i)
- << " Time: " << GetTime(i) << endl;
+ cout << "==================================" << endl;
+ cout << "TRACKID #" << GetTrackID(i) << ": Particle |" << GetParticleName(i) << "|"
+ << " Q1: " << GetQ1(i) << " Time: " << GetTime(i) << endl;
}
}
diff --git a/NPLib/Detectors/Epic/TEpicData.h b/NPLib/Detectors/Epic/TEpicData.h
index 4d612196b47ce90ec0ccd436ebd8e7b6af418f8d..f1489b5dbd47060f90e8f66e8d4a24a5c0990c1f 100644
--- a/NPLib/Detectors/Epic/TEpicData.h
+++ b/NPLib/Detectors/Epic/TEpicData.h
@@ -33,23 +33,24 @@ using namespace std;
class TEpicData : public TObject {
public:
+ // per TrackID
struct EpicAnodeData {
+ int TrackID;
+ string ParticleName;
UShort_t AnodeNbr;
Double_t Q1;
Double_t Time;
// vector of {Zstep, DEstep, DTstep}
- std::vector<string> ParticleName;
- std::vector<int> TrackID;
std::vector<Double_t> PosZ;
std::vector<Double_t> DE;
- std::vector<Double_t> DT; // DT=Ti-T0, with T0=Tfission or Talpha
+ std::vector<Double_t> DT; // DT=Ti-T0, with T0=Tfission or Tfirst_alpha
// Getters
+ int GetTrackID() const { return TrackID; }
+ string GetParticleName() const { return ParticleName; }
UShort_t GetAnodeNbr() const { return AnodeNbr; }
Double_t GetQ1() const { return Q1; }
Double_t GetTime() const { return Time; }
- const std::vector<string>& GetParticleName() const { return ParticleName; }
- const std::vector<int>& GetTrackID() const { return TrackID; }
const std::vector<Double_t>& GetPosZ() const { return PosZ; }
const std::vector<Double_t>& GetDE() const { return DE; }
const std::vector<Double_t>& GetDT() const { return DT; }
@@ -83,21 +84,19 @@ class TEpicData : public TObject {
// add //! to avoid ROOT creating dictionnary for the methods
public:
////////////////////// SETTERS ////////////////////////
- void Set(const UShort_t& n, const Double_t& q1, const Double_t& t,
- const std::vector<string>& name, const std::vector<int>& tid,
- const std::vector<double>& z, const std::vector<double>& de,
- const std::vector<double>& dt) {
- fEpic_Data.push_back({n, q1, t, name, tid, z, de, dt});
+ void Set(const int& tid, const string& name, const UShort_t& n, const Double_t& q1, const Double_t& t,
+ const std::vector<double>& z, const std::vector<double>& de, const std::vector<double>& dt) {
+ fEpic_Data.push_back({tid, name, n, q1, t, z, de, dt});
}
const EpicAnodeData& operator[](const unsigned int& i) const {return fEpic_Data[i];}
//////////////////////// GETTERS ////////////////////////
- inline UShort_t GetMultiplicity() const {return fEpic_Data.size();}
+ inline UShort_t GetMultiplicity() const {return fEpic_Data.size();}
+ int GetTrackID(const unsigned int& i) const { return fEpic_Data[i].TrackID; }
+ string GetParticleName(const unsigned int& i) const { return fEpic_Data[i].ParticleName; }
UShort_t GetAnodeNbr(const unsigned short& i) const { return fEpic_Data[i].AnodeNbr; }
Double_t GetQ1(const unsigned int& i) const { return fEpic_Data[i].Q1; }
Double_t GetTime(const unsigned int& i) const { return fEpic_Data[i].Time; }
- const std::vector<string>& GetParticleName(const unsigned int& i) const { return fEpic_Data[i].ParticleName; }
- const std::vector<int>& GetTrackID(const unsigned int& i) const { return fEpic_Data[i].TrackID; }
const std::vector<Double_t>& GetPosZ(const unsigned int& i) const { return fEpic_Data[i].PosZ; }
const std::vector<Double_t>& GetDE(const unsigned int& i) const { return fEpic_Data[i].DE; }
const std::vector<Double_t>& GetDT(const unsigned int& i) const { return fEpic_Data[i].DT; }
diff --git a/NPLib/Detectors/Epic/TEpicPhysics.cxx b/NPLib/Detectors/Epic/TEpicPhysics.cxx
index 1a3015acd3656060aaa23698f90de659caf70fb6..77a85451a86c21ec1a22febc0b2693ff3403d7d0 100644
--- a/NPLib/Detectors/Epic/TEpicPhysics.cxx
+++ b/NPLib/Detectors/Epic/TEpicPhysics.cxx
@@ -113,9 +113,9 @@ void TEpicPhysics::PreTreat() {
double TimeOffset = 0;
TimeOffset = Cal->GetValue("Epic/ANODE"+NPL::itoa(AnodeNumber)+"_TIMEOFFSET",0);
double Time = m_EventData->GetTime(i);
- m_PreTreatedData->Set(AnodeNumber, Q1, Time,
+ m_PreTreatedData->Set(m_EventData->GetTrackID(i),
m_EventData->GetParticleName(i),
- m_EventData->GetTrackID(i),
+ AnodeNumber, Q1, Time,
m_EventData->GetPosZ(i),
m_EventData->GetDE(i),
m_EventData->GetDT(i));
@@ -194,12 +194,7 @@ void TEpicPhysics::ReadAnalysisConfig() {
void TEpicPhysics::Clear() {
AnodeNumber.clear();
Q1.clear();
- Q2.clear();
- Qmax.clear();
Time.clear();
- Time_HF.clear();
- ToF.clear();
- isFakeFission.clear();
DT_FC.clear();
}
diff --git a/NPLib/Detectors/Epic/TEpicPhysics.h b/NPLib/Detectors/Epic/TEpicPhysics.h
index 6223dc2942ddcf45182de5446c0a9d94157a1043..8733238496be6f7ba2ecf7d4b07beb6b650a1c46 100644
--- a/NPLib/Detectors/Epic/TEpicPhysics.h
+++ b/NPLib/Detectors/Epic/TEpicPhysics.h
@@ -64,12 +64,7 @@ class TEpicPhysics : public TObject, public NPL::VDetector {
public:
vector<int> AnodeNumber;
vector<double> Q1;
- vector<double> Q2;
- vector<double> Qmax;
vector<double> Time;
- vector<bool> isFakeFission;
- vector<double> Time_HF;
- vector<double> ToF;
vector<double> DT_FC;
/// A usefull method to bundle all operation to add a detector
diff --git a/NPSimulation/Detectors/Epic/Epic.cc b/NPSimulation/Detectors/Epic/Epic.cc
index 539aa377a9f799c3231f0b4d170eee34dc7c5d4c..f3673dff1d39bb5f224bb965c0fea3b435298052 100644
--- a/NPSimulation/Detectors/Epic/Epic.cc
+++ b/NPSimulation/Detectors/Epic/Epic.cc
@@ -544,7 +544,7 @@ void Epic::ReadSensitive(const G4Event* ){
vector<unsigned int> level = Scorer->GetLevel(i);
//double Time = RandGauss::shoot(Scorer->GetTime(i),Epic_NS::ResoTime);
//double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Epic_NS::ResoEnergy);
- double Energy = Scorer->GetEnergy(i); // Attention Energy of all particles (FF, alpha, e-, ...)
+ //double Energy = Scorer->GetEnergy(i); // Attention Energy of all particles (FF, alpha, e-, ...)
vector<string> step_name = Scorer->GetParticleName(i);
vector<int> step_trackID = Scorer->GetTrackID(i);
vector<double> step_posZ = Scorer->GetStepPosZ(i);
@@ -556,56 +556,75 @@ void Epic::ReadSensitive(const G4Event* ){
double posZ_anode = (1-m_nA)*(m_Distance_AK*mm + m_Thickness_K*mm) - std::trunc(0.5*m_nA)*(m_InterDistance_KK*mm + thickness_A);
posZ_anode += m_mapping_A[Anode] * (2. * (m_Distance_AK*mm + m_Thickness_K*mm) + m_InterDistance_KK*mm + thickness_A) ;
- //// FIXME : need a fixed binning in dz of 25*um
- //// FIXME : one input per TrackID -> change the TEpicData TrackID and name should not be a vector anymore
- //vector<string> name;
- //vector<int> trackID;
- //vector<double> dz;
- //vector<double> de;
- //vector<double> dt;
-
- //double influence_pertrackID = 0 ;
- //bool end_of_new_trackID = false ;
- //cout << "step_name.size() = " << step_name.size() << endl;
- //for(int j=0; j<step_name.size(); j++){
- // cout << "j = " << j << " : TrackID = " << step_trackID.at(j) << " [previous is " << previous_trackID << "]. Particle = " << step_name.at(j) << " Actual m_Event mult : " << m_Event->GetMultiplicity() << endl;
- // if(step_name.at(j)!="e-" || step_name.at(j)=="anti_nu_e"){
- // continue;
- // }
- // else if(step_name.at(j)=="alpha"){ // FF or alpha // to be fixed !
- //
- // cout << "j = " << j << " : TrackID = " << step_trackID.at(j) << " [previous is " << previous_trackID << "]. Particle = " << step_name.at(j) << " Actual m_Event mult : " << m_Event->GetMultiplicity() << endl;
- //
- // if(step_trackID.at(j) != previous_trackID){
- // //if(m_Event->GetMultiplicity()>0) end_of_new_trackID = true;
- // if(previous_trackID != -1) end_of_new_trackID = true;
- // previous_trackID = step_trackID.at(j) ;
- // }
- //
- // if(end_of_new_trackID == true) {
- // m_Event->Set(m_mapping_A[Anode], // set anode number
- // influence_pertrackID, // set Q1 FF or alpha only
- // step_time.at(j), // set Time
- // name,trackID, dz, de, dt);
- // name.clear();
- // trackID.clear();
- // dz.clear();
- // de.clear();
- // dt.clear();
- // end_of_new_trackID = false;
- // }
-
- // double dz_anode = TMath::Abs(step_posZ.at(j) - posZ_anode) - 0.5*thickness_A ; // need to remove half thickness of Anode
- // influence_pertrackID += (step_DE.at(j) * dz_anode) / m_Distance_AK*mm;
- //
- // name.push_back(step_name.at(j));
- // trackID.push_back(step_trackID.at(j));
- // dz.push_back(dz_anode);
- // de.push_back(step_DE.at(j));
- // dt.push_back(step_time.at(j));
- // }// end of else
- //}// end of loop over j
+ //// FIXME : need a fixed binning in dz of 25*um and in time
+ //// FIXME : one TrackID per input for TEpicData
+ //// FIXME : fill another event class FasterData with only the vectors sorted in time
+
+
+ // === Epic Data
+ int trackID = -1;
+ string name_pertrackID = "";
+ double influence_pertrackID = 0;
+ double time_pertrackID = -1;
+ vector<double> dz;
+ vector<double> de;
+ vector<double> dt;
+
+ bool end_of_new_trackID = false ;
+ cout << "step_name.size() = " << step_name.size() << endl;
+ for(int j=0; j<step_name.size(); j++){
+
+ //if(step_name.at(j)!="e-" && step_name.at(j)!="anti_nu_e"){
+ // cout << " ======================================================================== " << endl;
+ // cout << "STEP #" << j << " Particle |" << step_name.at(j) << "|" << endl;
+ // cout << " " << " TrackID " << step_trackID.at(j) << " [previous " << previous_trackID
+ // << "] Actual m_Event mult : " << m_Event->GetMultiplicity()
+ // << " GetStepTime = " << step_time.at(j) << endl;
+ //}
+ if(step_trackID.at(j) != previous_trackID){
+ //if(m_Event->GetMultiplicity()>0) end_of_new_trackID = true;
+ if(previous_trackID != -1) end_of_new_trackID = true;
+ else {
+ trackID = step_trackID.at(j);
+ name_pertrackID = step_name.at(j);
+ time_pertrackID = step_time.at(j);
+ }
+ previous_trackID = step_trackID.at(j) ;
+ }
+
+ if(end_of_new_trackID == true || j==step_name.size()-1) {
+ if(name_pertrackID!="e-" && name_pertrackID!="anti_nu_e"){
+ //cout << "TEpicData m_Event->Set() : TrackID #" << trackID << " |" << name_pertrackID << "| influence = " << influence_pertrackID << " and time at " << time_pertrackID << endl;
+ m_Event->Set(trackID,
+ name_pertrackID,
+ m_mapping_A[Anode], // set anode number
+ influence_pertrackID, // set Q1 FF or alpha only
+ time_pertrackID, // set Time
+ dz, de, dt);
+ dz.clear();
+ de.clear();
+ dt.clear();
+ }
+ trackID = step_trackID.at(j);
+ name_pertrackID = step_name.at(j);
+ influence_pertrackID = 0;
+ time_pertrackID = step_time.at(j);
+ end_of_new_trackID = false;
+ }
+
+ double dz_anode = TMath::Abs(step_posZ.at(j) - posZ_anode) - 0.5*thickness_A ; // need to remove half thickness of Anode
+ influence_pertrackID += (step_DE.at(j) * dz_anode) / m_Distance_AK*mm;
+ dz.push_back(dz_anode);
+ de.push_back(step_DE.at(j));
+ dt.push_back(step_time.at(j));
+ }// end of loop over j
}// end of loop over i
+
+ //// FIXME : fill the FasterData
+ if(m_Event->GetMultiplicity()>0){
+ // fill clean-up data to FasterData
+ }
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 941c5e4be564c5f82947eae895a35a7f33b932d5..0fc833b1dd4f474186d1841d526474c0048d145f 100644
--- a/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
+++ b/Projects/AlPhaPha/2024/macro/chio/DE_E_Spline/SplineChio.C
@@ -8,19 +8,18 @@
#include <TH2.h>
#include <TSpline.h>
#include <fstream>
-#include "../XYCalibration/ProfileEvaluator.h"
using namespace std;
-const int Ncuts=30;
+
+const int CutNumberLoader();
+const int Ncuts = CutNumberLoader();
const int NSegment=11;
const int Nrun=1;
-int RunNumber[Nrun]={241};
+int RunNumber[Nrun]={247};
void MakeSpline();
void ApplySpline_perTCutG();
void ApplySpline();
-vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline);
-vector<double> TxtToVector(const char *Path);
///////////////////////////////////////////////////////////////////////////
void SplineChio()
@@ -31,76 +30,42 @@ void SplineChio()
//**********************Cut**************************************
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("Cut %i",i));
-
- // ****************** Spline Y **********************************
-
- 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");
-
-
- 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);
-
+ vector<TCutG*> cutZ(Ncuts);
+ for(int i=0;i<Ncuts;i++) cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
//**********************DE_E****************************************************
TChain *chain=new TChain("PhysicsTree");
- chain->Add("../../../root/analysis/VamosCalib247.root");
+ for(int i=0;i<Nrun;i++)
+ {
+ chain->Add(Form("../../../root/analysis/Run%d.root",RunNumber[i]));
+ cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
+ }
TICPhysics* IC = new TICPhysics() ;
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);
+ TH2F *hChioDE_E_all = new TH2F("hChioDE_E","hChioDE_E",1000,1000,36000,500,1000,26000);
+ vector<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,36000,1250,1000,26000);
auto start = std::chrono::high_resolution_clock::now();
-
- //int Nentries=chain->GetEntries();
- int Nentries=1000000;
-
+
+ int Nentries=chain->GetEntries();
+ //int Nentries=1000000;
+
for(int e=0;e<Nentries;++e)
{
if (e % 100000 == 0 && e > 0 ) {
@@ -114,78 +79,11 @@ void SplineChio()
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
-
-
-
- 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))) ;
-
- // ********************************* 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
-
- 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];
- }
-
- 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 *****************************
-
- 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){
- 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(IC->Eres,IC->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;}
- } // end max mult 1
+ //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(IC->Eres,IC->DE)) {hChioDE_E[i]->Fill(IC->Eres,IC->DE); break;}
}//end loop event
@@ -197,7 +95,7 @@ void SplineChio()
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);
- can->cd(); gPad-> SetLogz();
+ can->cd();
hChioDE_E_all->Draw("colz");
for(int i = 0 ; i < Ncuts ; i++) cutZ[i]->Draw("same");
}
@@ -205,14 +103,15 @@ void SplineChio()
///////////////////////////////////////////////////////////////////////////
void MakeSpline(){
- TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root");
- TSpline3 *gspline[Ncuts];
+ // ********************** Load prev histo****************************
+ TFile *inFile=new TFile("histo/SingleZ_ChioDE_E.root");
- TH2F* h2 [Ncuts];
- TH1F* hProfile[Ncuts];
- TH1F* pfx[Ncuts]; // extended hProfile
- TFile *fspline=new TFile("Output/spline_Chio_2024.root","recreate");
+ vector<TSpline3*> gspline(Ncuts);
+ vector<TH2F*> h2 (Ncuts);
+ vector<TH1F*> hProfile(Ncuts);
+ vector<TH1F*> pfx(Ncuts); // extended hProfile
+ TFile* fspline=new TFile("Output/spline_Chio_2024.root","recreate");
TCanvas * canfit = new TCanvas("canfit","canfit",0,0,2000,1500);
@@ -344,42 +243,34 @@ void ApplySpline_perTCutG()
// *************************Cut***************************
TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root");
- TCutG *cutZ[Ncuts];
+ vector<TCutG*> cutZ(Ncuts);
// ******************* TSpline DE*************************
TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
- TSpline3 *gspline[Ncuts];
+ vector<TSpline3*> gspline(Ncuts);
for(int i=0;i<Ncuts;i++){
cutZ[i]=(TCutG*)fcut->Get(Form("Cut %i",i));
gspline[i] = (TSpline3*)fspline->FindObjectAny(Form("fspline_%d",i+1));
}
- // ****************** Spline Y **********************************
- 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");
+ for(int i=0;i<Nrun;i++)
+ {
+ chain->Add(Form("../../../root/analysis/Run%d.root",RunNumber[i]));
+ cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
+ }
TICPhysics* IC = new TICPhysics() ;
- 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);
- // variable after applying TSpline
+ // variable after applying TSpline MUST FIND middles of each spline
double FF_DEcorr;
- double FF_DEcorr0[Ncuts] = {
+ vector<double> FF_DEcorr0 = {
7375.73,
7982.83,
8419.21,
@@ -421,6 +312,7 @@ TH2F *hChioDE_E = new TH2F("hChioDE_E","hChioDE_E",1000,1000,41000,500,1000
TH2F *hChioDE_E_corr = new TH2F("hChioDE_E_corr","hChioDE_E_corr",1000,1000,41000,500,1000,26000);
int Nentries=chain->GetEntries();
+//int Nentries = 1000000;
auto start = std::chrono::high_resolution_clock::now();
for(int e=0;e<Nentries;++e){
@@ -435,55 +327,11 @@ for(int e=0;e<Nentries;++e){
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
-
- 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);
- }
-
- 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 > 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) ;
- }
- }
-
- 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);
+ hChioDE_E_all->Fill(IC->Eres,IC->DE);
+ for(int i=0;i<Ncuts;i++) if(cutZ[i]->IsInside(IC->Eres,IC->DE)) {
+ hChioDE_E->Fill(IC->Eres,IC->DE);
+ FF_DEcorr = FF_DEcorr0[i] * IC->DE / gspline[i]->Eval(IC->Eres);
+ hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr);
break;
}
}
@@ -502,39 +350,29 @@ can->cd(3); gPad-> SetLogz(); hChioDE_E_corr->Draw("colz");
void ApplySpline()
{
// *************************Cut***************************
- TFile *fcut=new TFile("Cut/Cut_Z.root");
- TCutG *cutZ[Ncuts];
+ TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root", "open");
+ vector<TCutG*> cutZ(Ncuts);
// ******************* TSpline DE*************************
TFile *fspline=new TFile("Output/spline_Chio_2024.root","read");
- TSpline3 *gspline[Ncuts];
+ vector<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");
-
-
- vector<vector<TSpline3*>> Spline = LoadSpline("../YCalibration/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");
+ for(int i=0;i<Nrun;i++)
+ {
+ chain->Add(Form("../../../root/analysis/Run%d.root",RunNumber[i]));
+ cout << "Run number " << RunNumber[i] << " loaded sucessfully !" << endl;
+ }
TICPhysics* IC = new TICPhysics() ;
- 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);
// variable after applying TSpline
@@ -542,7 +380,7 @@ void ApplySpline()
double Zrough=0;
double Zabitbetter=0;
- double FF_DEcorr0[Ncuts];
+ vector<double> FF_DEcorr0(Ncuts);
for(int index=0; index<Ncuts; index++){
FF_DEcorr0[index] = gspline[index]->Eval(8500);
}
@@ -590,7 +428,8 @@ void ApplySpline()
TH1F *hChioZ_rough = new TH1F("hChioZ_rough","hChioZ_rough",2500,25,65);
TH1F *hChioZ_abitbetter = new TH1F("hChioZ_abitbetter","hChioZ_abitbetter",2500,25,65);
- int Nentries=1e6;//chain->GetEntries();
+ //int Nentries=1e6;
+ int Nentries = chain->GetEntries();
auto start = std::chrono::high_resolution_clock::now();
double FF_Eres_prev = 0;
@@ -608,89 +447,39 @@ void ApplySpline()
FF_DEcorr = -100;
chain->GetEntry(e);
- //===========================================================================================================
- // Splined Y
- //===========================================================================================================
- 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);
- }
-
- 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 > 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 ;
+ if (IC->DE<3000) continue;
+ if(IC->Eres==FF_Eres_prev) continue;
+ FF_Eres_prev = IC->Eres;
- 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) ;
- }
- }
-
- //===========================================================================================================
- // Spline DE
- //===========================================================================================================&
-
-
- if (FF_DE<3000) continue;
- if(FF_Eres==FF_Eres_prev) continue;
- FF_Eres_prev = FF_Eres;
-
- hChioDE_E_all->Fill(FF_Eres,FF_DE);
+ hChioDE_E_all->Fill(IC->Eres,IC->DE);
float Eval_DEspline, DEspline0;
int index=0;
for(int i=0; i<Ncuts; i++){
- Eval_DEspline = gspline[i]->Eval(FF_Eres);
- if(FF_DE<Eval_DEspline) break;
+ Eval_DEspline = gspline[i]->Eval(IC->Eres);
+ if(IC->DE<Eval_DEspline) break;
index = i;
}
- Eval_DEspline = gspline[index]->Eval(FF_Eres);
+ Eval_DEspline = gspline[index]->Eval(IC->Eres);
DEspline0 = FF_DEcorr0[index];
float dmin, dsup;
- if( index < (Ncuts-1) && FF_DE > gspline[0]->Eval(FF_Eres)){
- dmin = FF_DE-gspline[index]->Eval(FF_Eres);
- dsup = gspline[index+1]->Eval(FF_Eres)-FF_DE;
+ if( index < (Ncuts-1) && IC->DE > gspline[0]->Eval(IC->Eres)){
+ dmin = IC->DE-gspline[index]->Eval(IC->Eres);
+ dsup = gspline[index+1]->Eval(IC->Eres)-IC->DE;
if(dmin<0) cout << "negative value of dmin = " << dmin << ", for index = " << index << endl;
if(dsup<0) cout << "negative value of dsup = " << dsup << ", for index = " << index << endl;
//if(dsup<dmin) {
- // Eval_DEspline = gspline[index+1]->Eval(FF_Eres) ;
+ // Eval_DEspline = gspline[index+1]->Eval(IC->Eres) ;
// DEspline0 = FF_DEcorr0[index] ;
//}
- Eval_DEspline = dsup*gspline[index]->Eval(FF_Eres)/(dmin+dsup) + dmin*gspline[index+1]->Eval(FF_Eres)/(dmin+dsup) ;
+ Eval_DEspline = dsup*gspline[index]->Eval(IC->Eres)/(dmin+dsup) + dmin*gspline[index+1]->Eval(IC->Eres)/(dmin+dsup) ;
DEspline0 = dsup*FF_DEcorr0[index]/(dmin+dsup) + dmin*FF_DEcorr0[index+1]/(dmin+dsup) ;
- FF_DEcorr = DEspline0 * FF_DE / Eval_DEspline ;
- hChioDE_E_corr->Fill(FF_Eres,FF_DEcorr);
+ FF_DEcorr = DEspline0 * IC->DE / Eval_DEspline ;
+ hChioDE_E_corr->Fill(IC->Eres,FF_DEcorr);
//if(FF_DEcorr>15120 && FF_DEcorr<15130) {
- // cout << e << " " << index << " " << FF_Eres << " " << FF_DEcorr << endl;
+ // cout << e << " " << index << " " << IC->Eres << " " << FF_DEcorr << endl;
//}
}
@@ -698,9 +487,9 @@ void ApplySpline()
//Zrough = -110.165 + 3.34475*sqrt(FF_DEcorr) - 0.0271123*FF_DEcorr + 8.60752e-05 * pow(sqrt(FF_DEcorr),3);
Zrough = 16.8521 + 0.0017328*FF_DEcorr + 1.70774e-8*pow(FF_DEcorr,2);
Zabitbetter = -127.117 + 3.83463*sqrt(FF_DEcorr) - 0.0317448 *FF_DEcorr + 0.000100428 * pow(sqrt(FF_DEcorr),3);
- hChioZ_E_rough->Fill(FF_Eres,Zrough);
+ hChioZ_E_rough->Fill(IC->Eres,Zrough);
- if(3000<FF_Eres && FF_Eres<25000){
+ if(3000<IC->Eres && IC->Eres<25000){
hChioDE_corr->Fill(FF_DEcorr);
hChioZ_rough->Fill(Zrough);
hChioZ_abitbetter->Fill(Zabitbetter);
@@ -800,66 +589,21 @@ void Fit_DE_E_corr()
}
+const int CutNumberLoader(){
+ TFile *fcut=new TFile("../XYCalibration/Output/CutAutoZ.root","open");
-vector<vector<TSpline3*>> LoadSpline(const char* PathToSpline){
- TFile *fSplineIC = new TFile(PathToSpline,"open");
- // Get number of spline
- int SplineCount = 0 ;
- TIter next(fSplineIC->GetListOfKeys());
+ int CutCount = 0 ;
+ TIter next(fcut->GetListOfKeys());
TKey* key;
while ((key=(TKey*)next())){
- if (std::string(key->GetClassName()) == "TSpline3"){
- SplineCount ++;
+ if (std::string(key->GetClassName()) == "TCutG"){
+ CutCount ++;
}
}
- vector<vector<TSpline3*>> Spline;
- 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));
- }
- } //End loop on histogram
-
- Spline.push_back(spline_X);
- Spline.push_back(spline_Y);
-
- return Spline;
+ const int res = CutCount;
+ return res;
}
-
-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/Epic_sim_proto/Pu240.mac b/Projects/Epic_sim_proto/Pu240.mac
deleted file mode 100644
index 9a01f41a956e6169ef83b2c81c72e656d5408046..0000000000000000000000000000000000000000
--- a/Projects/Epic_sim_proto/Pu240.mac
+++ /dev/null
@@ -1 +0,0 @@
-/process/had/rdm/thresholdForVeryLongDecayTime 1.e+4 y
diff --git a/Projects/Epic_sim_proto/U238.mac b/Projects/Epic_sim_proto/U238.mac
deleted file mode 100644
index f89f13d8e44d1e42f9ddc5b01fd59b0ec957b8f3..0000000000000000000000000000000000000000
--- a/Projects/Epic_sim_proto/U238.mac
+++ /dev/null
@@ -1 +0,0 @@
-/process/had/rdm/thresholdForVeryLongDecayTime 1.e+20 y
diff --git a/Projects/nptuto_gdr/Analysis.cxx b/Projects/nptuto_gdr/Analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c71f7c79de96f97faecbbd846234cdea4bed1972
--- /dev/null
+++ b/Projects/nptuto_gdr/Analysis.cxx
@@ -0,0 +1,242 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include <iostream>
+using namespace std;
+
+#include "Analysis.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "NPFunction.h"
+#include "NPOptionManager.h"
+// #include <unistd.h>
+////////////////////////////////////////////////////////////////////////////////
+Analysis::Analysis() {}
+////////////////////////////////////////////////////////////////////////////////
+Analysis::~Analysis() {}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init() {
+ if (NPOptionManager::getInstance()->HasDefinition("simulation")) {
+ cout << "Considering input data as simulation" << endl;
+ simulation = true;
+ }
+ else {
+ cout << "Considering input data as real" << endl;
+ simulation = false;
+ }
+ simulation = true;
+
+ // initialize input and output branches
+ if (simulation) {
+ Initial = new TInitialConditions();
+ ReactionConditions = new TReactionConditions();
+ }
+
+ InitOutputBranch();
+ InitInputBranch();
+ // get MUST2 objects
+ M2 = (TMust2Physics*)m_DetectorManager->GetDetector("M2Telescope");
+
+ // get reaction information
+ reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+ OriginalBeamEnergy = reaction.GetBeamEnergy();
+ // target thickness
+ TargetThickness = m_DetectorManager->GetTargetThickness();
+ string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+
+ // energy losses
+ string light = NPL::ChangeNameToG4Standard(reaction.GetNucleus3()->GetName());
+ string beam = NPL::ChangeNameToG4Standard(reaction.GetNucleus1()->GetName());
+ cout << light << " " << beam << " " << TargetMaterial << " " << TargetThickness << endl;
+ LightTarget = NPL::EnergyLoss(light + "_" + TargetMaterial + ".G4table", "G4Table", 100);
+ LightAl = NPL::EnergyLoss(light + "_Al.G4table", "G4Table", 100);
+ LightSi = NPL::EnergyLoss(light + "_Si.G4table", "G4Table", 100);
+ BeamTarget = NPL::EnergyLoss(beam + "_" + TargetMaterial + ".G4table", "G4Table", 100);
+
+ FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, TargetThickness * 0.5, 0);
+ // FinalBeamEnergy = OriginalBeamEnergy;
+ // cout << "Original beam energy: " << OriginalBeamEnergy << " MeV Mid-target beam energy: " << FinalBeamEnergy
+ // << "MeV " << endl;
+ reaction.SetBeamEnergy(FinalBeamEnergy);
+
+ // initialize various parameters
+ Rand = TRandom3();
+ DetectorNumber = 0;
+ ThetaNormalTarget = 0;
+ ThetaM2Surface = 0;
+ Si_E_M2 = 0;
+ CsI_E_M2 = 0;
+ Energy = 0;
+ ThetaGDSurface = 0;
+ X = 0;
+ Y = 0;
+ Z = 0;
+ dE = 0;
+ BeamDirection = TVector3(0, 0, 1);
+ nbTrack = 0;
+
+ sleep(5);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::TreatEvent() {
+ // Reinitiate calculated variable
+ ReInitValue();
+ double XTarget, YTarget;
+ TVector3 BeamDirection;
+ XTarget = 0;
+ YTarget = 0;
+ BeamDirection = TVector3(0, 0, 1);
+ BeamImpact = TVector3(XTarget, YTarget, 0);
+
+ // determine beam energy for a randomized interaction point in target
+ // 1% FWHM randominastion (E/100)/2.35
+ // reaction.SetBeamEnergy(Rand.Gaus(BeamEnergy, BeamEnergy * 1. / 100. / 2.35));
+ // reaction.SetBeamEnergy(BeamEnergy);
+
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////// LOOP on MUST2 ////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////
+ for (unsigned int countMust2 = 0; countMust2 < M2->Si_E.size(); countMust2++) {
+ /************************************************/
+ // Part 0 : Get the usefull Data
+ // MUST2
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
+ TVector3 HitDirection = M2->GetPositionOfInteraction(countMust2) - BeamImpact;
+ ThetaLab = HitDirection.Angle(BeamDirection);
+
+ X = M2->GetPositionOfInteraction(countMust2).X();
+ Y = M2->GetPositionOfInteraction(countMust2).Y();
+ Z = M2->GetPositionOfInteraction(countMust2).Z();
+
+ ThetaM2Surface = HitDirection.Angle(-M2->GetTelescopeNormal(countMust2));
+ ThetaNormalTarget = HitDirection.Angle(TVector3(0, 0, 1));
+
+ /************************************************/
+
+ /************************************************/
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2 = M2->CsI_E[countMust2];
+
+ // if CsI
+ if (CsI_E_M2 > 0) {
+ // The energy in CsI is calculate form dE/dx Table because
+ Energy = CsI_E_M2;
+ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
+ Energy += Si_E_M2;
+ }
+ else {
+ Energy = Si_E_M2;
+ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
+ }
+
+ // Evaluate energy using the thickness
+ ELab = Energy;
+
+ // Target Correction
+ ELab = LightTarget.EvaluateInitialEnergy(ELab, TargetThickness * 0.5, ThetaNormalTarget);
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ // Ex = reaction.ReconstructRelativistic(ELab, ThetaLab);
+ // reaction.SetBeamEnergy(Initial->GetIncidentFinalKineticEnergy());
+ Ex = reaction.ReconstructRelativistic(ELab, ThetaLab);
+ // ExNoBeam = reaction.ReconstructRelativistic(ELab, ThetaLab);
+ // reaction.SetBeamEnergy(FinalBeamEnergy);
+ // ExNoProton = reaction.ReconstructRelativistic(ReactionConditions->GetKineticEnergy(0),
+ // ReactionConditions->GetParticleDirection(0).Angle(TVector3(0, 0,
+ // 1)));
+ ThetaLab = ThetaLab / deg;
+
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = reaction.EnergyLabToThetaCM(ELab, ThetaLab) / deg;
+ /************************************************/
+ } // end loop MUST2
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::End() {}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex, "Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab, "ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab, "ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM/D");
+ RootOutput::getInstance()->GetTree()->Branch("Run", &Run, "Run/I");
+ RootOutput::getInstance()->GetTree()->Branch("X", &X, "X/D");
+ RootOutput::getInstance()->GetTree()->Branch("Y", &Y, "Y/D");
+ RootOutput::getInstance()->GetTree()->Branch("Z", &Z, "Z/D");
+ RootOutput::getInstance()->GetTree()->Branch("dE", &dE, "dE/D");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch() {
+ // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
+ if (!simulation) {
+ }
+ else {
+ RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions", true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fIC_*", true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions", &Initial);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("ReactionConditions", true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*", true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &ReactionConditions);
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue() {
+ Ex = -1000;
+ ELab = -1000;
+ ThetaLab = -1000;
+ ThetaCM = -1000;
+ X = -1000;
+ Y = -1000;
+ Z = -1000;
+ dE = -1000;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the AnalysisFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VAnalysis* Analysis::Construct() { return (NPL::VAnalysis*)new Analysis(); }
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C" {
+class proxy_analysis {
+ public:
+ proxy_analysis() { NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct); }
+};
+
+proxy_analysis p_analysis;
+}
diff --git a/Projects/nptuto_gdr/Analysis.h b/Projects/nptuto_gdr/Analysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..f314fb353d1faae7ca8ba4c63f9536243880dd03
--- /dev/null
+++ b/Projects/nptuto_gdr/Analysis.h
@@ -0,0 +1,114 @@
+#ifndef Analysis_h
+#define Analysis_h
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include "NPEnergyLoss.h"
+#include "NPReaction.h"
+#include "NPVAnalysis.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TCATSPhysics.h"
+#include "TInitialConditions.h"
+#include "TMust2Physics.h"
+#include "TReactionConditions.h"
+#include <TMath.h>
+#include <TRandom3.h>
+#include <TVector3.h>
+
+class Analysis : public NPL::VAnalysis {
+ public:
+ Analysis();
+ ~Analysis();
+
+ public:
+ void Init();
+ void TreatEvent();
+ void End();
+
+ void InitOutputBranch();
+ void InitInputBranch();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
+
+ private:
+ double Ex;
+ double ELab;
+ double ThetaLab;
+ double ThetaCM;
+
+ double TargetThickness;
+ NPL::Reaction reaction;
+ // Energy loss table: the G4Table are generated by the simulation
+ NPL::EnergyLoss LightTarget;
+ NPL::EnergyLoss LightAl;
+ NPL::EnergyLoss LightSi;
+ NPL::EnergyLoss BeamTarget;
+
+ // Beam Energy
+ double OriginalBeamEnergy; // AMEV
+ double FinalBeamEnergy;
+
+ // intermediate variable
+ TVector3 BeamDirection;
+ TVector3 BeamImpact;
+ TRandom3 Rand;
+ int Run;
+ int DetectorNumber;
+ double ThetaNormalTarget;
+ double ThetaM2Surface;
+ double Si_E_M2;
+ double CsI_E_M2;
+ double Energy;
+ double BeamEnergy;
+
+ double ThetaGDSurface;
+ double X;
+ double Y;
+ double Z;
+ // Vamos Branches
+ unsigned long long int LTS;
+
+ // Agata branches
+ double agata_zShift;
+ unsigned long long int TStrack;
+ int nbHits;
+ int nbTrack;
+ float* trackE = new float(100);
+ float* trackX1 = new float(100);
+ float* trackY1 = new float(100);
+ float* trackZ1 = new float(100);
+ float* trackT = new float(100);
+ int* trackCrystalID = new int(100);
+ int nbCores;
+ int* coreId = new int(100);
+ ULong64_t* coreTS = new ULong64_t(100);
+ float* coreE0 = new float(100);
+ //
+ double dE;
+ double dTheta;
+ // Branches and detectors
+ TMust2Physics* M2;
+ bool simulation;
+ TInitialConditions* Initial;
+ TReactionConditions* ReactionConditions;
+};
+#endif
diff --git a/Projects/nptuto_gdr/CMakeLists.txt b/Projects/nptuto_gdr/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..22c74affdfc45019bdda2594f8439c52d4ab97ec
--- /dev/null
+++ b/Projects/nptuto_gdr/CMakeLists.txt
@@ -0,0 +1,5 @@
+cmake_minimum_required (VERSION 2.8)
+# Setting the policy to match Cmake version
+cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
+# include the default NPAnalysis cmake file
+include("../../NPLib/ressources/CMake/NPAnalysis.cmake")
diff --git a/Projects/nptuto_gdr/configs/ConfigMust2.dat b/Projects/nptuto_gdr/configs/ConfigMust2.dat
new file mode 100644
index 0000000000000000000000000000000000000000..216b16b66a665ababbd3ffe819626fbe65c5e89e
--- /dev/null
+++ b/Projects/nptuto_gdr/configs/ConfigMust2.dat
@@ -0,0 +1,5 @@
+ConfigMust2
+ CSI_SIZE 40
+ CAL_PIXEL
+ PIXEL_SIZE 8
+ MAX_STRIP_MULTIPLICITY 10
diff --git a/Projects/nptuto_gdr/detector/MUGAST_LISE.detector b/Projects/nptuto_gdr/detector/MUGAST_LISE.detector
new file mode 100644
index 0000000000000000000000000000000000000000..f03211b192acbea5191614709a79dcb4569db21b
--- /dev/null
+++ b/Projects/nptuto_gdr/detector/MUGAST_LISE.detector
@@ -0,0 +1,56 @@
+%%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+Target
+ THICKNESS= 50 micrometer
+ ANGLE= 0 deg
+ RADIUS= 25 mm
+ MATERIAL= CH2
+ X= 0
+ Y= 0
+ Z= 0
+
+%%%%%%% Telescope 1 %%%%%%%
+M2Telescope
+ X128_Y128= 115.88 9.61 154.54 mm
+ X128_Y1= 104.8 101.89 125.09 mm
+ X1_Y1= 14.55 102.4 160.63 mm
+ X1_Y128= 25.63 10.12 190.08 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 2 %%%%%%%
+M2Telescope
+ X128_Y128= -11.23 102.42 160.87 mm
+ X128_Y1= -101.39 102.39 124.37 mm
+ X1_Y1= -113.17 10.36 153.56 mm
+ X1_Y128= -23.03 10.38 190.05 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 3 %%%%%%%
+M2Telescope
+ X128_Y128= -113.28 -12.52 153.32 mm
+ X128_Y1= -101.58 -104.77 124.82 mm
+ X1_Y1= -11.39 -104.58 161.48 mm
+ X1_Y128= -23.1 -12.34 189.98 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X128_Y128= 13.82 -104.92 160.72 mm
+ X128_Y1= 104.3 -104.95 125.08 mm
+ X1_Y1= 115.75 -12.73 153.76 mm
+ X1_Y128= 25.23 -12.65 189.43 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/nptuto_gdr/reactions/48Crpd.reaction b/Projects/nptuto_gdr/reactions/48Crpd.reaction
new file mode 100644
index 0000000000000000000000000000000000000000..ef5c42de50074fa865e5825f3a9630e28eb0722c
--- /dev/null
+++ b/Projects/nptuto_gdr/reactions/48Crpd.reaction
@@ -0,0 +1,32 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 48Cr
+ Energy= 1440 MeV
+ SigmaEnergy= 40 MeV
+ ExcitationEnergy= 0 MeV
+ SigmaThetaX= 0 mrad
+ SigmaPhiY= 0 mrad
+ SigmaX= 1 mm
+ SigmaY= 1 mm
+ MeanThetaX= 0 mrad
+ MeanPhiY= 0 mrad
+ MeanX= 0 mm
+ MeanY= 0 mm
+ %EnergyProfilePath= eLise.txt EL
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 48Cr
+ Target= 1H
+ Light= 2H
+ Heavy= 47Cr
+ ExcitationEnergyLight= 0.0 MeV
+ ExcitationEnergyHeavy= 0.0 MeV
+ CrossSectionPath= flat.txt
+ ShootLight= 1
+ ShootHeavy= 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+