From 6cc34af28c9b8874e468d58415ae17989a292e1a Mon Sep 17 00:00:00 2001
From: Hugo Jacob <hugojacob57@gmail.com>
Date: Fri, 26 Jan 2024 15:14:03 +0100
Subject: [PATCH] Updating E805 Analysis
---
NPLib/Detectors/Exogam/TExogamPhysics.cxx | 33 +++++++++++++------
NPLib/Detectors/Exogam/TExogamPhysics.h | 8 +++--
NPLib/Physics/NPReaction.cxx | 9 ++++++
NPLib/Physics/NPReaction.h | 3 ++
Projects/E805/Analysis.cxx | 39 +++++++++++++++++++----
Projects/E805/Analysis.h | 7 ++++
6 files changed, 80 insertions(+), 19 deletions(-)
diff --git a/NPLib/Detectors/Exogam/TExogamPhysics.cxx b/NPLib/Detectors/Exogam/TExogamPhysics.cxx
index ed0f8f2a2..523dbfd04 100644
--- a/NPLib/Detectors/Exogam/TExogamPhysics.cxx
+++ b/NPLib/Detectors/Exogam/TExogamPhysics.cxx
@@ -180,20 +180,31 @@ void TExogamPhysics::BuildPhysicalEvent() {
// Adding all AddBack (AB) related stuff
E_AB.push_back(E_AddBack);
- FlangeN_AB.push_back(flange_nbr);
- Size_AB.push_back((*it).second.size());
+ FlangeN_ABD.push_back(flange_nbr);
+ Size_ABD.push_back((*it).second.size());
// Adding these parameters for Doppler correction purposes (D)
CrystalN_ABD.push_back(crystal_nbr);
int MaxOuterId = GetMaxOuter(Id_Max);
OuterN_ABD.push_back(GetMaxOuter(Id_Max));
-
- // If a max Outer is found, Do doppler correction, else push_back -1000;
- double EnergyDoppler = -1000;
+
if(MaxOuterId > -1){
- EnergyDoppler= GetDoppler(E_AddBack, flange_nbr, crystal_nbr, MaxOuterId);
+ Exogam_struc = Ask_For_Angles(flange_nbr, ComputeMeanFreePath(E_AddBack));
+ double Theta_seg = Exogam_struc.Theta_Crystal_Seg[crystal_nbr][MaxOuterId];
+ double Phi_seg = Exogam_struc.Phi_Crystal_Seg[crystal_nbr][MaxOuterId];
+ Theta_D.push_back(Theta_seg);
+ Phi_D.push_back(Phi_seg);
}
- E_ABD.push_back(EnergyDoppler);
+ else{
+ Theta_D.push_back(-1000);
+ Phi_D.push_back(-1000);
+ }
+ // If a max Outer is found, Do doppler correction, else push_back -1000;
+ //double EnergyDoppler = -1000;
+ //if(MaxOuterId > -1){
+ // EnergyDoppler= GetDoppler(E_AddBack, flange_nbr, crystal_nbr, MaxOuterId);
+ //}
+ //E_ABD.push_back(EnergyDoppler);
}
}
@@ -286,11 +297,13 @@ void TExogamPhysics::Clear() {
CrystalN.clear();
// E_Doppler.clear();
E_AB.clear();
- FlangeN_AB.clear();
- Size_AB.clear();
+ FlangeN_ABD.clear();
+ Size_ABD.clear();
CrystalN_ABD.clear();
OuterN_ABD.clear();
- E_ABD.clear();
+ Theta_D.clear();
+ Phi_D.clear();
+ // E_ABD.clear();
//
// ECC_CloverNumber.clear();
diff --git a/NPLib/Detectors/Exogam/TExogamPhysics.h b/NPLib/Detectors/Exogam/TExogamPhysics.h
index a923e2cf4..3e42c7672 100644
--- a/NPLib/Detectors/Exogam/TExogamPhysics.h
+++ b/NPLib/Detectors/Exogam/TExogamPhysics.h
@@ -75,11 +75,13 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
// Previous treatment + Add_Back (size of vectors are not the same because of AB !)
std::vector<double> E_AB;
- std::vector<unsigned int> FlangeN_AB;
- std::vector<unsigned int> Size_AB;
+ std::vector<unsigned int> FlangeN_ABD;
+ std::vector<unsigned int> Size_ABD;
std::vector<unsigned int> CrystalN_ABD;
std::vector<int> OuterN_ABD;
- std::vector<double> E_ABD;
+ std::vector<double> Theta_D;
+ std::vector<double> Phi_D;
+ // std::vector<double> E_ABD;
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index 46aefd569..e79397f08 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -309,6 +309,15 @@ double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab, doub
return Eex;
}
+TLorentzVector Reaction::LorentzAfterReaction(double EnergyLab, double ThetaLab, double PhiLab){
+ double E3 = m3 + EnergyLab;
+ double p_Lab_3 = sqrt(E3 * E3 - m3 * m3);
+ fEnergyImpulsionLab_3 =
+ TLorentzVector(p_Lab_3 * sin(ThetaLab) * cos(PhiLab), p_Lab_3 * sin(PhiLab), p_Lab_3 * cos(ThetaLab), E3);
+ fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
+ return fEnergyImpulsionLab_4;
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Return ThetaCM
double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab) {
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 455aadcc6..71388278f 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -243,6 +243,9 @@ namespace NPL {
// Return Excitation Energy
double ReconstructRelativistic(double EnergyLab, double ThetaLab, double PhiLab = 0);
+
+ // Return Lorentz Vector of Heavy Recoil after reaction
+ TLorentzVector LorentzAfterReaction(double EnergyLab, double ThetaLab, double PhiLab = 0);
// Return ThetaCM
// EnergyLab: energy measured in the laboratory frame
diff --git a/Projects/E805/Analysis.cxx b/Projects/E805/Analysis.cxx
index 8fd78683a..3a757df1b 100755
--- a/Projects/E805/Analysis.cxx
+++ b/Projects/E805/Analysis.cxx
@@ -104,11 +104,11 @@ void Analysis::TreatEvent(){
TreatCATS();
if(bCATS){
TreatMUST2();
+ TreatEXO();
}
//if(bCATS){
// TreatZDD();
// TreatTAC();
- // TreatEXO();
//}
/*//for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
//Si_E_M2 = M2->Si_E[countMust2];
@@ -166,8 +166,8 @@ void Analysis::TreatMUST2(){
ThetaM2Surface = 0;
ThetaNormalTarget = 0;
- BeamImpact = TVector3(0,0,0);
- BeamDirection = TVector3(0,0,1);
+ //BeamImpact = TVector3(0,0,0);
+ //BeamDirection = TVector3(0,0,1);
TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact;
M2_ThetaLab.push_back(HitDirection.Angle( BeamDirection ));
//std::cout << BeamImpact.X() << " " << BeamImpact.Y() << " " << BeamImpact.Z() << std::endl;
@@ -217,6 +217,14 @@ void Analysis::TreatMUST2(){
M2_Ex_t.push_back(Reaction_pt->ReconstructRelativistic( Energy["triton"] , M2_ThetaLab[countMust2] ));
M2_Ex_a.push_back(Reaction_p3He->ReconstructRelativistic( Energy["alpha"] , M2_ThetaLab[countMust2] ));
+ TLorentzVector PHeavy_pd = Reaction_pd->LorentzAfterReaction(Energy["deuteron"] , M2_ThetaLab[countMust2]);
+ TLorentzVector PHeavy_pt = Reaction_pd->LorentzAfterReaction(Energy["triton"] , M2_ThetaLab[countMust2]);
+ TLorentzVector PHeavy_p3He = Reaction_pd->LorentzAfterReaction(Energy["alpha"] , M2_ThetaLab[countMust2]);
+ Beta_pd.push_back(PHeavy_pd.Beta());
+ Beta_pt.push_back(PHeavy_pt.Beta());
+ Beta_p3He.push_back(PHeavy_p3He.Beta());
+
+
M2_CsI_E_p.push_back(CsI_Energy["proton"]);
M2_CsI_E_d.push_back(CsI_Energy["deuteron"]);
M2_CsI_E_t.push_back(CsI_Energy["triton"]);
@@ -232,6 +240,15 @@ void Analysis::TreatMUST2(){
}
void Analysis::TreatEXO(){
+ int EXO_AB_size = Exogam->E_AB.size();
+ for(unsigned int countExo = 0 ; countExo < EXO_AB_size; countExo++){
+ // Doing Doppler correction only if one reaction occurs
+ if(Beta_pd.size() == 1){
+ EXO_Doppler_pd.push_back(Doppler_Correction(Exogam->Theta_D[countExo], Exogam->Phi_D[countExo], 0,0,Beta_pd[0],Exogam->E_AB[countExo]));
+ EXO_Doppler_pt.push_back(Doppler_Correction(Exogam->Theta_D[countExo], Exogam->Phi_D[countExo], 0,0,Beta_pt[0],Exogam->E_AB[countExo]));
+ EXO_Doppler_p3He.push_back(Doppler_Correction(Exogam->Theta_D[countExo], Exogam->Phi_D[countExo], 0,0,Beta_p3He[0],Exogam->E_AB[countExo]));
+ }
+ }
// }
/*
if(M2->Si_E.size() > 0){
@@ -299,9 +316,16 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("M2_X",&M2_X);
RootOutput::getInstance()->GetTree()->Branch("M2_Y",&M2_Y);
RootOutput::getInstance()->GetTree()->Branch("M2_Z",&M2_Z);
- RootOutput::getInstance()->GetTree()->Branch("M2_dE",&M2_dE);
- RootOutput::getInstance()->GetTree()->Branch("CsI_E_M2",&CsI_E_M2);
- // RootOutput::getInstance()->GetTree()->Branch("M2_ECsI_from_deltaE",&M2_ECsI_from_deltaE);
+ // RootOutput::getInstance()->GetTree()->Branch("M2_dE",&M2_dE);
+ // RootOutput::getInstance()->GetTree()->Branch("CsI_E_M2",&CsI_E_M2);
+
+ RootOutput::getInstance()->GetTree()->Branch("EXO_Doppler_pd",&EXO_Doppler_pd);
+ RootOutput::getInstance()->GetTree()->Branch("EXO_Doppler_pt",&EXO_Doppler_pt);
+ RootOutput::getInstance()->GetTree()->Branch("EXO_Doppler_p3He",&EXO_Doppler_p3He);
+
+ RootOutput::getInstance()->GetTree()->Branch("Beta_pd",&Beta_pd);
+ RootOutput::getInstance()->GetTree()->Branch("Beta_p3He",&Beta_p3He);
+ RootOutput::getInstance()->GetTree()->Branch("Beta_pt",&Beta_pt);
/*
RootOutput:: getInstance()->GetTree()->Branch("TAC_CATS_PL",&TAC_CATS_PL,"TAC_CATS_PL/s");
RootOutput:: getInstance()->GetTree()->Branch("TAC_CATS_PLTS",&TAC_CATS_PLTS,"TAC_CATS_PLTS/l");
@@ -466,6 +490,9 @@ void Analysis::ReInit(){
M2_Z.clear();
M2_dE.clear();
+ EXO_Doppler_p3He.clear();
+ EXO_Doppler_pt.clear();
+ EXO_Doppler_pd.clear();
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/E805/Analysis.h b/Projects/E805/Analysis.h
index c08cf7c7f..271d55d2a 100755
--- a/Projects/E805/Analysis.h
+++ b/Projects/E805/Analysis.h
@@ -118,6 +118,13 @@ class Analysis: public NPL::VAnalysis{
std::vector<double> M2_Y;
std::vector<double> M2_Z;
std::vector<double> M2_dE;
+
+ std::vector<double> EXO_Doppler_pd;
+ std::vector<double> EXO_Doppler_pt;
+ std::vector<double> EXO_Doppler_p3He;
+ std::vector<double> Beta_pd;
+ std::vector<double> Beta_pt;
+ std::vector<double> Beta_p3He;
double OriginalBeamEnergy ; // AMEV
--
GitLab