From de9bec2fbc258da016bcc8b82eb4999b6caf7024 Mon Sep 17 00:00:00 2001
From: Hugo Jacob <hugojacob57@gmail.com>
Date: Tue, 17 Oct 2023 15:01:19 +0200
Subject: [PATCH] Updating Analysis E805
---
NPLib/Detectors/CATS/TCATSPhysics.cxx | 11 +-------
NPLib/Detectors/CATS/TCATSPhysics.h | 14 +++++++---
NPLib/Detectors/MUST2/TMust2Physics.cxx | 1 -
NPLib/Utility/npreader.cxx | 17 ++-----------
Projects/E805/Analysis.cxx | 34 ++++++++++++++++++++-----
Projects/E805/Analysis.h | 12 ++++++---
6 files changed, 50 insertions(+), 39 deletions(-)
diff --git a/NPLib/Detectors/CATS/TCATSPhysics.cxx b/NPLib/Detectors/CATS/TCATSPhysics.cxx
index 9d3dc3d83..78ad6c461 100644
--- a/NPLib/Detectors/CATS/TCATSPhysics.cxx
+++ b/NPLib/Detectors/CATS/TCATSPhysics.cxx
@@ -115,7 +115,6 @@ void TCATSPhysics::BuildSimplePhysicalEvent(){
//////////////////////////////////////////////////////////////////////////////
void TCATSPhysics::BuildPhysicalEvent(){
-
// std::cout << "test 1" << std::endl;
if (NPOptionManager::getInstance()->IsReader() == true) {
m_EventData = &(**r_ReaderEventData);
@@ -143,7 +142,6 @@ void TCATSPhysics::BuildPhysicalEvent(){
DetectorHit.insert(m_PreTreatedData->GetCATSDetY(i));
}
// The number of CATS hit, i.e. the number of CATS that we are going to analyse
- sizeDet = DetectorHit.size();
@@ -185,7 +183,6 @@ void TCATSPhysics::BuildPhysicalEvent(){
// std::cout << "test 7" << std::endl;
double PosY = ReconstructionFunctionY[DetN](MaxQY[DetN],MapY[DetN], QSumY[DetN]);
//std::cout << "test " << std::reduce(QsumSample[DetN].begin(),QsumSample[DetN].end()) /(QsumSample[DetN]).size() << std::endl;
- //std::cout << "test Pos " << PosX << " " << PosY << std::endl;
StripNumberX.push_back(PosX);
StripNumberY.push_back(PosY);
DetNumber.push_back(DetN);
@@ -217,11 +214,6 @@ void TCATSPhysics::BuildPhysicalEvent(){
PositionY.push_back(py0+(py1-py0)*(PosY-sy0));
PositionZ.push_back(StripPositionZ[DetN]);
}
- else{
- PositionX.push_back(-1000);
- PositionY.push_back(-1000);
- PositionZ.push_back(-1000);
- }
}
// Sorting Positions depending on Z
@@ -262,7 +254,6 @@ void TCATSPhysics::BuildPhysicalEvent(){
BeamDirection = GetBeamDirection();
// Does not meet the conditions for target position and beam direction
-
return;
}
@@ -461,7 +452,7 @@ void TCATSPhysics::Clear(){
MapY.clear();
MaxQX.clear();
MaxQY.clear();
-
+
DetectorHit.clear();
DetectorHitX.clear();
}
diff --git a/NPLib/Detectors/CATS/TCATSPhysics.h b/NPLib/Detectors/CATS/TCATSPhysics.h
index d5aec479f..6e1e74046 100644
--- a/NPLib/Detectors/CATS/TCATSPhysics.h
+++ b/NPLib/Detectors/CATS/TCATSPhysics.h
@@ -88,7 +88,7 @@ class TCATSPhysics : public TObject, public NPL::VDetector, public TCATSPhysicsR
std::map<UShort_t,vector< vector<double> > > StripPositionX;//!
std::map<UShort_t,vector< vector<double> > > StripPositionY;//!
std::map<UShort_t,double> StripPositionZ;//!
- int m_NumberOfCATS;
+ int m_NumberOfCATS; //!
double m_TargetAngle; //!
double m_TargetThickness; //!
@@ -104,7 +104,7 @@ class TCATSPhysics : public TObject, public NPL::VDetector, public TCATSPhysicsR
public:
// Set the reconstruction Method used for the CATS plane
- void SetReconstructionMethod(unsigned int CATSNumber, string XorY, string MethodName);
+ void SetReconstructionMethod(unsigned int CATSNumber, string XorY, string MethodName); //!
private :
// Map of activated channel
@@ -116,8 +116,8 @@ class TCATSPhysics : public TObject, public NPL::VDetector, public TCATSPhysicsR
public: // Output data of interest
// for a CATS
- void SetTargetAngle(double m_TargetAngle) {m_TargetAngle = m_TargetAngle;}
- void SetTargetThickness(double m_TargetThickness) {m_TargetThickness = m_TargetThickness;}
+ void SetTargetAngle(double m_TargetAngle) {m_TargetAngle = m_TargetAngle;} //!
+ void SetTargetThickness(double m_TargetThickness) {m_TargetThickness = m_TargetThickness;} //!
// Remove bad channel, calibrate the data and apply threshold
@@ -214,6 +214,12 @@ class TCATSPhysics : public TObject, public NPL::VDetector, public TCATSPhysicsR
std::set<int> DetectorHit;//!
std::map<UShort_t, std::vector<double>> QsumSample;//!
+ //Debugging
+ unsigned long counter;//!
+ unsigned long long time_elapsed1;//!
+ unsigned long long time_elapsed2;//!
+ unsigned long long time_elapsed3;//!
+
private: // Spectra Class
TCATSSpectra* m_Spectra;//!
diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx
index 102c1e072..f9de9e77b 100644
--- a/NPLib/Detectors/MUST2/TMust2Physics.cxx
+++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx
@@ -326,7 +326,6 @@ void TMust2Physics::BuildSimplePhysicalEvent() { BuildPhysicalEvent(); }
///////////////////////////////////////////////////////////////////////////
void TMust2Physics::BuildPhysicalEvent() {
-
if (NPOptionManager::getInstance()->IsReader() == true) {
m_EventData = &(**r_ReaderEventData);
}
diff --git a/NPLib/Utility/npreader.cxx b/NPLib/Utility/npreader.cxx
index 670091b24..63bd65cf0 100644
--- a/NPLib/Utility/npreader.cxx
+++ b/NPLib/Utility/npreader.cxx
@@ -77,21 +77,14 @@ int main(int argc , char** argv){
// Instantiate the detector using a file
NPL::DetectorManager* myDetector = new NPL::DetectorManager();
- std::cout << "test npa1" << std::endl;
myDetector->ReadConfigurationFile(detectorfileName);
- std::cout << "test npa2" << std::endl;
myDetector->InitializeRootInput();
- std::cout << "test npa3" << std::endl;
myDetector->InitializeRootOutput();
- std::cout << "test npa4" << std::endl;
// Attempt to load an analysis
NPL::VAnalysis* UserAnalysis = NULL;
std::string libName = "./libNPAnalysis" + myOptionManager->GetSharedLibExtension();
- std::cout << "test npa5" << std::endl;
dlopen(libName.c_str(),RTLD_NOW | RTLD_GLOBAL);
- std::cout << "test npa5.1" << std::endl;
char* error = dlerror();
- std::cout << "test npa5.2" << std::endl;
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
std::cout << "Checking TreeReader adress: " << inputTreeReader << std::endl;
myDetector->SetTreeReader(inputTreeReader);
@@ -112,7 +105,6 @@ int main(int argc , char** argv){
}
}
- std::cout << "test npa6" << std::endl;
if(myOptionManager->GetOnline()){
// Request Detector manager to give the Spectra to the server
myDetector->SetSpectraServer();
@@ -121,7 +113,6 @@ int main(int argc , char** argv){
std::cout << std::endl << "///////// Starting Analysis ///////// "<< std::endl;
TChain* Chain = RootInput:: getInstance()->GetChain();
myOptionManager->GetNumberOfEntryToAnalyse();
- std::cout << "test npa7" << std::endl;
unsigned long first_entry = myOptionManager->GetFirstEntryToAnalyse(); // defaults to zero
unsigned long nentries = Chain->GetEntries();
@@ -144,14 +135,11 @@ int main(int argc , char** argv){
unsigned long new_nentries = 0 ;
int current_tree = 0 ;
int total_tree = Chain->GetNtrees();
- std::cout << "test npa8" << std::endl;
bool IsPhysics = myOptionManager->GetInputPhysicalTreeOption();
- std::cout << "test npa9" << std::endl;
if(UserAnalysis==NULL){
if(!IsPhysics){
- std::cout << "test npa branch 1" << std::endl;
while(inputTreeReader->Next()){
//for (unsigned long i = first_entry ; i < nentries + first_entry; i++) {
// // Get the raw Data
@@ -196,15 +184,14 @@ int main(int argc , char** argv){
else{
if(!IsPhysics){
- std::cout << "test npa branch 2" << std::endl;
while(inputTreeReader->Next()){
// Build the current event
if(UserAnalysis->UnallocateBeforeBuild()){
myDetector->BuildPhysicalEvent();
- // User Analysis
+ // User Analysis;
if(UserAnalysis->UnallocateBeforeTreat()){
- UserAnalysis->TreatEvent();
+ UserAnalysis->TreatEvent();
//Fill the tree
if(UserAnalysis->FillOutputCondition())
myDetector->FillOutputTree();
diff --git a/Projects/E805/Analysis.cxx b/Projects/E805/Analysis.cxx
index 0ed9df108..1ca337ef6 100755
--- a/Projects/E805/Analysis.cxx
+++ b/Projects/E805/Analysis.cxx
@@ -49,11 +49,25 @@ void Analysis::Init(){
string heavy_ejectile= NPL::ChangeNameToG4Standard(reaction->GetNucleus4()->GetName());
string light=NPL::ChangeNameToG4Standard(reaction->GetNucleus3()->GetName());
-
+ string Reaction_pd_s = "48Cr(p,d)47Cr@1620";
+ string Reaction_pt_s = "48Cr(p,t)46Cr@1620";
+ string Reaction_p3He_s = "48Cr(p,3He)46V@1620";
+ Reaction_pd = new Reaction(Reaction_pd_s);
+ Reaction_pt = new Reaction(Reaction_pt_s);
+ Reaction_p3He = new Reaction(Reaction_p3He_s);
//
ProtonSi = NPL::EnergyLoss(Path+ "proton_Si.G4table", "G4Table", 100);
+
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+ LightAl[ParticleType[i]] = NPL::EnergyLoss(Path+ParticleType[i]+"_Al.G4table","G4Table",100);
+ LightTarget[ParticleType[i]] = NPL::EnergyLoss(Path+ParticleType[i]+"_CH2.G4table","G4Table",100);
+ }
+ BeamTarget["48Cr"] = NPL::EnergyLoss(Path+"Cr48_CH2.G4table","G4Table",100);
+ Reaction_pd->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_pd->GetBeamEnergy(),TargetThickness*0.5,0));
+ Reaction_pt->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_pt->GetBeamEnergy(),TargetThickness*0.5,0));
+ Reaction_p3He->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_p3He->GetBeamEnergy(),TargetThickness*0.5,0));
Cal = CalibrationManager::getInstance();
}
///////////////////////////// Initialize some important parameters //////////////////////////////////
@@ -140,10 +154,12 @@ void Analysis::TreatEvent(){
ThetaNormalTarget = 0;
BeamImpact = TVector3(CATS->PositionOnTargetX,CATS->PositionOnTargetY,0);
- // BeamImpact = TVector3(0,0,0);
+ // std::cout << "Position On target : " << CATS->PositionOnTargetX << " " << CATS->PositionOnTargetY << std::endl;
+
+ //BeamImpact = TVector3(0,0,0);
BeamDirection = TVector3(CATS->PositionX[1] - CATS->PositionX[0],CATS->PositionY[1] - CATS->PositionY[0],CATS->PositionZ[1] - CATS->PositionZ[0]);
- // std::cout << CATS->PositionX[0] - CATS->PositionX[1] << " " << CATS->PositionY[0] - CATS->PositionY[1] << " " << CATS->PositionZ[0] - CATS->PositionZ[1] << std::endl;
+ // std::cout << "Position XY " << CATS->PositionX[1] - CATS->PositionX[0] << " " << CATS->PositionY[1] - CATS->PositionY[0] << " " << CATS->PositionZ[1] - CATS->PositionZ[0] << std::endl;
// BeamDirection = TVector3(0,0,1);
TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
@@ -172,6 +188,7 @@ void Analysis::TreatEvent(){
// The energy in CsI is calculate form dE/dx Table because
CsI_Energy[ParticleType[i]] = Cal->ApplyCalibration("MUST2/"+ParticleType[i]+"_T"+NPL::itoa(TelescopeNumber)+"_CsI"+NPL::itoa(CristalNb)+"_E",CsI_E_M2);
Energy[ParticleType[i]] = CsI_Energy[ParticleType[i]];
+ // Energy[ParticleType[i]] = LightAl[ParticleType[i]].EvaluateInitialEnergy(Energy[ParticleType[i]], 0.4 * micrometer, ThetaM2Surface);
//std::cout << ParticleType[i]+"MUST2/T"+NPL::itoa(TelescopeNumber)+"_CsI"+NPL::itoa(CristalNb)+"_E" << " " << Energy[ParticleType[i]] << "\n";
//Energy = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
//Energy+=Si_E_M2;
@@ -179,6 +196,11 @@ void Analysis::TreatEvent(){
// else
Energy[ParticleType[i]] += Si_E_M2;
+ if(Energy[ParticleType[i]] > 0)
+ Energy[ParticleType[i]] = LightTarget[ParticleType[i]].EvaluateInitialEnergy(Energy[ParticleType[i]] ,TargetThickness*0.5, ThetaNormalTarget);
+ else
+ Energy[ParticleType[i]] = -1000;
+ // std::cout << "DILO tata " << Energy[ParticleType[i]] << std::endl;
}
@@ -191,9 +213,9 @@ void Analysis::TreatEvent(){
// Part 3 : Excitation Energy Calculation
M2_Ex_p.push_back(reaction->ReconstructRelativistic( Energy["proton"] , M2_ThetaLab[countMust2] ));
- M2_Ex_d.push_back(reaction->ReconstructRelativistic( Energy["deuteron"] , M2_ThetaLab[countMust2] ));
- M2_Ex_t.push_back(reaction->ReconstructRelativistic( Energy["triton"] , M2_ThetaLab[countMust2] ));
- M2_Ex_a.push_back(reaction->ReconstructRelativistic( Energy["alpha"] , M2_ThetaLab[countMust2] ));
+ M2_Ex_d.push_back(Reaction_pd->ReconstructRelativistic( Energy["deuteron"] , M2_ThetaLab[countMust2] ));
+ 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] ));
M2_CsI_E_p.push_back(CsI_Energy["proton"]);
M2_CsI_E_d.push_back(CsI_Energy["deuteron"]);
diff --git a/Projects/E805/Analysis.h b/Projects/E805/Analysis.h
index ccaa73730..903c00996 100755
--- a/Projects/E805/Analysis.h
+++ b/Projects/E805/Analysis.h
@@ -280,6 +280,9 @@ class Analysis: public NPL::VAnalysis{
double CsI_E_M2 ;
std::vector<string> ParticleType{"proton","deuteron","triton","alpha"};
std::map<TString, double> Energy ;
+ std::map<TString, NPL::EnergyLoss> LightAl ;
+ std::map<TString, NPL::EnergyLoss> LightTarget ;
+ std::map<TString, NPL::EnergyLoss> BeamTarget ;
std::map<TString, double> CsI_Energy ;
double BeamEnergy;
double ThetaGDSurface ;
@@ -300,19 +303,22 @@ class Analysis: public NPL::VAnalysis{
+ NPL::Reaction* Reaction_pd;
+ NPL::Reaction* Reaction_pt;
+ NPL::Reaction* Reaction_p3He;
NPL::Reaction* reaction = new Reaction;
NPL::Particle* Co_57 = new Particle;
NPL::Particle* Ni_57 = new Particle;
NPL::EnergyLoss Beam_Target;
NPL::EnergyLoss Heavy_Target;
- NPL::EnergyLoss LightTarget;
+ // NPL::EnergyLoss LightTarget;
NPL::EnergyLoss ProtonSi;
std::vector<NPL::EnergyLoss> Heavy_IC_Gas;
std::vector<NPL::EnergyLoss> Heavy_IC_Windows;
std::vector<NPL::EnergyLoss> Heavy_IC_Mylar;
std::vector<NPL::EnergyLoss> Heavy_DC_Gas;
- NPL::EnergyLoss LightAl;
- NPL::EnergyLoss LightSi;
+ // NPL::EnergyLoss LightAl;
+ // NPL::EnergyLoss LightSi;
string BeamName = "48Cr";
--
GitLab