From c46231e2318da4dbec174e778102f93d5f6218bf Mon Sep 17 00:00:00 2001
From: Morfouace <pierre.morfouace@gmail.com>
Date: Thu, 22 Dec 2022 10:49:21 +0100
Subject: [PATCH] Updating ProcessScorer for Vendeta

---
 Inputs/EventGenerator/neutron.source      |  2 +-
 NPSimulation/Detectors/Vendeta/Vendeta.cc | 29 +++++++-
 NPSimulation/Scorers/ProcessScorers.cc    | 86 ++++++++++++++---------
 NPSimulation/Scorers/ProcessScorers.hh    |  2 +
 Projects/Vendeta_sim/run.mac              |  2 +-
 5 files changed, 83 insertions(+), 38 deletions(-)

diff --git a/Inputs/EventGenerator/neutron.source b/Inputs/EventGenerator/neutron.source
index 4f7489c76..775ce6b40 100644
--- a/Inputs/EventGenerator/neutron.source
+++ b/Inputs/EventGenerator/neutron.source
@@ -13,7 +13,7 @@ Isotropic
  %EnergyDistribution= 0.619676*TMath::SinH(sqrt(1.07777*x))*exp(-0.847212*x)
  %EnergyDistribution= 1.5*TMath::SinH(sqrt(1.3*x))*exp(-0.89*x)
  HalfOpenAngleMin= 0
- HalfOpenAngleMax= 180
+ HalfOpenAngleMax= 2
  x0= 0 
  y0= 0 
  z0= 0 mm 
diff --git a/NPSimulation/Detectors/Vendeta/Vendeta.cc b/NPSimulation/Detectors/Vendeta/Vendeta.cc
index acbb987d5..af97c3a19 100644
--- a/NPSimulation/Detectors/Vendeta/Vendeta.cc
+++ b/NPSimulation/Detectors/Vendeta/Vendeta.cc
@@ -40,6 +40,7 @@
 
 // NPTool header
 #include "Vendeta.hh"
+#include "ProcessScorers.hh"
 #include "CalorimeterScorers.hh"
 #include "InteractionScorers.hh"
 #include "RootOutput.h"
@@ -75,7 +76,7 @@ Vendeta::Vendeta(){
   m_VendetaDetector = 0;
   m_SensitiveCell = 0;
   m_MecanicalStructure = 0;  
-  m_Build_MecanicalStructure = 1;
+  m_Build_MecanicalStructure = 0;
 
   // RGB Color + Transparency
   m_VisAl      = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5));   
@@ -351,7 +352,6 @@ void Vendeta::ReadSensitive(const G4Event* ){
       m_Event->SetHGQ2(Energy);
       m_Event->SetHGTime(Time); 
       m_Event->SetHGQmax(0); 
-      m_Event->SetHGIsSat(0); 
 
       // Filling LG
       m_Event->SetLGDetectorNbr(DetectorNbr);
@@ -359,12 +359,33 @@ void Vendeta::ReadSensitive(const G4Event* ){
       m_Event->SetLGQ2(Energy);
       m_Event->SetLGTime(Time); 
       m_Event->SetLGQmax(0); 
-      m_Event->SetLGIsSat(0); 
 
     }
   }
+
+  ///////////
+  // Process scorer
+  ProcessScorers::PS_Process* Process_scorer = (ProcessScorers::PS_Process*) m_VendetaScorer->GetPrimitive(2);
+  unsigned int ProcessMult = Process_scorer->GetMult();
+  if(ProcessMult>0){
+    string particle_name = Process_scorer->GetParticleName(0);
+    if(particle_name=="gamma"){
+      m_Event->SetHGIsSat(1); 
+      m_Event->SetLGIsSat(1); 
+    }
+
+    if(particle_name=="neutron"){
+      m_Event->SetHGIsSat(0); 
+      m_Event->SetLGIsSat(0); 
+    }
+ 
+  }
+  //for(unsigned int i=0; i<ProcessMult; i++){
+  //  string particle_name = Process_scorer->GetParticleName(i);
+  //}
 }
 
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 ////////////////////////////////////////////////////////////////   
 void Vendeta::InitializeScorers() { 
@@ -379,9 +400,11 @@ void Vendeta::InitializeScorers() {
   vector<int> level; level.push_back(0);
   G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter",level, 0) ;
   G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ;
+  G4VPrimitiveScorer* Process= new ProcessScorers::PS_Process("Process", 0) ;
   //and register it to the multifunctionnal detector
   m_VendetaScorer->RegisterPrimitive(Calorimeter);
   m_VendetaScorer->RegisterPrimitive(Interaction);
+  m_VendetaScorer->RegisterPrimitive(Process);
   G4SDManager::GetSDMpointer()->AddNewDetector(m_VendetaScorer) ;
 }
 
diff --git a/NPSimulation/Scorers/ProcessScorers.cc b/NPSimulation/Scorers/ProcessScorers.cc
index 301442aa4..334a86a9d 100644
--- a/NPSimulation/Scorers/ProcessScorers.cc
+++ b/NPSimulation/Scorers/ProcessScorers.cc
@@ -27,8 +27,8 @@ 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("Process", &t_processname);
+  //auto tree = RootOutput::getInstance()->GetTree();
+  //tree->Branch("ProcessBranch", &t_processname);
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -37,13 +37,32 @@ PS_Process::~PS_Process() {}
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 G4bool PS_Process::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
   // Contain Process Name
-  G4String processname;
+  G4String process_name;
+  G4String particle_name;
+  G4String volume_name;
+  //G4int trackID;
+  //G4int parentID;
   if (aStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL) {
-    processname = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
-
-    t_processname.push_back(processname);
+    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;
@@ -51,32 +70,32 @@ G4bool PS_Process::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
   //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);
-      }
-      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());
-        }
-      }
-    }
+  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);
+  }
+  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;
@@ -90,6 +109,7 @@ void PS_Process::EndOfEvent(G4HCofThisEvent*) {}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 void PS_Process::clear() {
+  t_particlename.clear();
   t_processname.clear();
   t_processtime.clear();
   t_gamma_energy.clear();
diff --git a/NPSimulation/Scorers/ProcessScorers.hh b/NPSimulation/Scorers/ProcessScorers.hh
index c0de3b6b5..7382a9550 100644
--- a/NPSimulation/Scorers/ProcessScorers.hh
+++ b/NPSimulation/Scorers/ProcessScorers.hh
@@ -59,6 +59,7 @@ namespace ProcessScorers {
       //double t_Time;
       //vector<unsigned int> t_Level;
       vector<G4String> t_processname;
+      vector<G4String> t_particlename;
       vector<double> t_processtime;
       vector<double> t_gamma_energy;
       vector<double> t_proton_energy;
@@ -69,6 +70,7 @@ namespace ProcessScorers {
     public:
       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> GetProtonEnergy() {return t_proton_energy;};
diff --git a/Projects/Vendeta_sim/run.mac b/Projects/Vendeta_sim/run.mac
index 325d88eef..b379ed9bc 100644
--- a/Projects/Vendeta_sim/run.mac
+++ b/Projects/Vendeta_sim/run.mac
@@ -1 +1 @@
-/run/beamOn 1000000
+/run/beamOn 10000
-- 
GitLab