From 509ad94acb8acc83da4766ad3cab5452e9708229 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Mon, 27 Mar 2023 11:11:22 +0200
Subject: [PATCH] * Importing modification to Process Scorer made by Louis
 Lemair for nebula

---
 NPSimulation/Scorers/ProcessScorers.cc | 152 +++++++++++++++----------
 NPSimulation/Scorers/ProcessScorers.hh |  17 +--
 2 files changed, 103 insertions(+), 66 deletions(-)

diff --git a/NPSimulation/Scorers/ProcessScorers.cc b/NPSimulation/Scorers/ProcessScorers.cc
index 334a86a9d..8d8c4fd73 100644
--- a/NPSimulation/Scorers/ProcessScorers.cc
+++ b/NPSimulation/Scorers/ProcessScorers.cc
@@ -27,8 +27,14 @@ using namespace ProcessScorers;
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 PS_Process::PS_Process(G4String name, G4int depth) : G4VPrimitiveScorer(name, depth) {
   // m_NestingLevel = NestingLevel;
-  //auto tree = RootOutput::getInstance()->GetTree();
-  //tree->Branch("ProcessBranch", &t_processname);
+  auto tree = RootOutput::getInstance()->GetTree();
+  tree->Branch("ProcScor_Particle", &t_particlename);
+  tree->Branch("ProcScor_Process", &t_processname);
+  tree->Branch("ProcScor_Time", &t_processtime);
+  tree->Branch("ProcScor_Efficiency", &t_efficiency);
+  tree->Branch("ProcScor_El1_InEl_2", &t_FC_process);
+  tree->Branch("ProcScor_Pos_Z", &t_pos_Z);
+  tree->Branch("ProcScor_KinEn", &t_kinetic_energy);
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -37,66 +43,93 @@ PS_Process::~PS_Process() {}
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 G4bool PS_Process::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
   // Contain Process Name
-  G4String process_name;
-  G4String particle_name;
-  G4String volume_name;
-  //G4int trackID;
-  //G4int parentID;
+  G4String processname;
+  double pos_z;
   if (aStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL) {
-    process_name  = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
-    particle_name = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
-    volume_name   = aStep->GetTrack()->GetVolume()->GetName();
-    //parentID      = aStep->GetTrack()->GetParentID();
-    //trackID       = aStep->GetTrack()->GetTrackID();
-
-    t_processname.push_back(process_name);
-    t_particlename.push_back(particle_name);
-    t_processtime.push_back(aStep->GetPreStepPoint()->GetGlobalTime());
- 
-    /*if(particle_name=="gamma"){
-      cout << "**************" << endl;
-      cout << "size= " << t_particlename.size() << endl;
-      cout << "process name= " << process_name << endl;
-      cout << "particle name= " << particle_name << endl;
-      cout << "volume name= " << volume_name << endl;
-      cout << "parent ID= " << parentID << endl;
-      cout << "track ID= " << trackID << endl;
-    }*/
- }
-
-  /*G4int trackID;
-  //G4int parentID;
+    processname = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
+    t_processname.push_back(processname);
+    if(processname!="Transportation"){
+      pos_z = aStep->GetPostStepPoint()->GetPosition()[2];
+      t_pos_Z.push_back(pos_z);
+    }
+    
+    //cout << "Process :" << processname << endl;
+    t_processtime.push_back(aStep->GetPostStepPoint()->GetGlobalTime());
+  }
+
+  G4int trackID;
+  G4int parentID;
   G4String volumename, particlename;
   //G4double step_length;
   G4double track_kineE;
   if(aStep->GetTrack()->GetNextVolume() != 0){
-  trackID = aStep->GetTrack()->GetTrackID();
-  //parentID = aStep->GetTrack()->GetParentID();
-  volumename = aStep->GetTrack()->GetVolume()->GetName();
-  particlename = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
-  //step_length = aStep->GetTrack()->GetStepLength();
-  track_kineE = aStep->GetTrack()->GetKineticEnergy();
-
-  if(volumename=="av_1_impr_1_gas_box_logic_pv_0"){
-  if(processname=="hadElastic")t_FC_process.push_back(1);
-  if(processname=="neutronInelastic")t_FC_process.push_back(2);
-  }
-
-  if(HasBeenTracked[trackID]==0){
-  if(particlename=="gamma"){
-  HasBeenTracked[trackID]=1;
-  //cout << trackID  << " " << track_kineE << endl;
-  t_gamma_energy.push_back(track_kineE);
+    trackID = aStep->GetTrack()->GetTrackID();
+    parentID = aStep->GetTrack()->GetParentID();
+    volumename = aStep->GetTrack()->GetVolume()->GetName();
+    particlename = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
+    //cout << "Particule :" << particlename << endl;
+    //step_length = aStep->GetTrack()->GetStepLength();
+    track_kineE = aStep->GetTrack()->GetKineticEnergy();
+    t_particlename.push_back(particlename);
+    t_kinetic_energy.push_back(track_kineE);
+
+    //Looking after the neutrons 
+    if(particlename=="neutron"){
+
+      //If it's an original neutron
+      if(parentID==0){
+        //Checking what reactions it prefers between El and InEl 
+        if(processname=="hadElastic")t_FC_process.push_back(5);
+        if(processname=="neutronInelastic")t_FC_process.push_back(6);
+        //Checking wall efficiency by discriminating whether it has interactied or not 
+        if(HasBeenTracked.find(trackID)==HasBeenTracked.end()){//first interaction
+          if(processname!="Transportation"){
+            // # SPECIAL NEBULA PROJECT #
+            if(volumename=="NebulaModule" || volumename=="NebulaVeto"){
+              //if the hit happens in Nebula, it counts as efficiency
+              HasBeenTracked.insert(trackID);
+              t_efficiency.push_back(1);
+            }
+            else{
+              //the purpose of this is for dicrimination, 0 being "didn't interact", 1 "interacted", 5 "interacted elsewhere"
+              HasBeenTracked.insert(trackID);
+              t_efficiency.push_back(5);
+            }
+            //Checking if first interaction is El or InEl
+            if(processname=="hadElastic")t_FC_process.push_back(1);
+            if(processname=="neutronInelastic")t_FC_process.push_back(2);
+          }
+        }
+      }
+      //Otherwise it's a neutron coming from secondary reactions
+      else{
+        if(processname=="hadElastic")t_FC_process.push_back(9);
+        if(processname=="neutronInelastic")t_FC_process.push_back(10);
+      }
+    }//endif neutron
+
+
+    if(HasBeenTracked.find(trackID)==HasBeenTracked.end()){
+      if(particlename=="gamma"){
+        HasBeenTracked.insert(trackID);
+        //cout << trackID  << " " << track_kineE << endl;
+        t_gamma_energy.push_back(track_kineE);
+      }
+      if(particlename=="proton"){
+        HasBeenTracked.insert(trackID);
+        if(track_kineE>0.1){
+          t_proton_energy.push_back(track_kineE);
+          t_proton_time.push_back(aStep->GetPreStepPoint()->GetGlobalTime());
+        }
+      }
+    }//endif first interaction (other than neutron)
+  }//endif track unended
+  else{
+    trackID = aStep->GetTrack()->GetTrackID();
+    if(HasBeenTracked.find(trackID)==HasBeenTracked.end()){
+      t_efficiency.push_back(0);
+    }
   }
-  if(particlename=="proton"){
-  HasBeenTracked[trackID]=1;
-  if(track_kineE>0.1){
-  t_proton_energy.push_back(track_kineE);
-  t_proton_time.push_back(aStep->GetPreStepPoint()->GetGlobalTime());
-  }
-  }
-  }
-  }*/
 
   return TRUE;
 }
@@ -116,10 +149,11 @@ void PS_Process::clear() {
   t_proton_energy.clear();
   t_proton_time.clear();
   t_FC_process.clear();
+  t_efficiency.clear();
+  t_pos_Z.clear();
+  t_kinetic_energy.clear();
 
-  for (unsigned int i = 0; i < 100; i++) {
-    HasBeenTracked[i] = 0;
-  }
+  HasBeenTracked.clear();
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Scorers/ProcessScorers.hh b/NPSimulation/Scorers/ProcessScorers.hh
index 7382a9550..af460d038 100644
--- a/NPSimulation/Scorers/ProcessScorers.hh
+++ b/NPSimulation/Scorers/ProcessScorers.hh
@@ -49,30 +49,33 @@ namespace ProcessScorers {
       void DrawAll();
       void PrintAll();
 
-    private: // How much level of volume nesting should be considered
+    private: 
+      // How much level of volume nesting should be considered
       // Give the list of the nesting level at which the copy number should be return.
       // 0 is the lowest level possible (the actual volume copy number in which the interaction happen)
 
     private: 
-      //CalorimeterDataVector m_Data;
-      //double t_Energy;
       //double t_Time;
       //vector<unsigned int> t_Level;
-      vector<G4String> t_processname;
-      vector<G4String> t_particlename;
+      vector<string> t_particlename;
+      vector<string> t_processname;
       vector<double> t_processtime;
       vector<double> t_gamma_energy;
       vector<double> t_proton_energy;
       vector<double> t_proton_time;
       vector<int> t_FC_process;
+      vector<int> t_efficiency;
+      vector<double> t_pos_Z;
+      vector<double> t_kinetic_energy;
 
-      int HasBeenTracked[100];
+      std::set<int> HasBeenTracked; //map<ID,HasBeenTrackedStatus>
     public:
-      inline unsigned int  GetMult() {return t_processname.size();};
+      inline unsigned int GetMult() {return t_processname.size();};
       inline string GetProcessName(const unsigned int& i) {return t_processname[i];};
       inline string GetParticleName(const unsigned int& i) {return t_particlename[i];};
       inline double GetProcessTime(const unsigned int& i) {return t_processtime[i];};
       inline vector<double> GetGammaEnergy() {return t_gamma_energy;};
+      inline vector<double> GetKinEnergy() {return t_kinetic_energy;};
       inline vector<double> GetProtonEnergy() {return t_proton_energy;};
       inline vector<double> GetProtonTime() {return t_proton_time;};
       inline vector<int> GetFCProcess() {return t_FC_process;};
-- 
GitLab