Track the occurrence of gamma hits in ChiNu. Volume name from ProcessScorers stored.

NPLib/Detectors/ChiNu/TChiNuData.cxx

@@ -50,6 +50,9 @@ void TChiNuData::Clear() {
   // Time
   fChiNu_T_DetectorNbr.clear();
   fChiNu_Time.clear();
+  // PArticle type
+  fChiNu_IG_DetectorNbr.clear();
+  fChiNu_IsGamma.clear();
 }
 
 

NPLib/Detectors/ChiNu/TChiNuData.h

@@ -42,6 +42,10 @@ class TChiNuData : public TObject {
     vector<UShort_t>   fChiNu_T_DetectorNbr;
     vector<Double_t>   fChiNu_Time;
 
+    // Particle Type
+    vector<UShort_t>   fChiNu_IG_DetectorNbr;
+    vector<Bool_t>     fChiNu_IsGamma;
+
 
   //////////////////////////////////////////////////////////////
   // Constructor and destructor
@@ -77,6 +81,12 @@ class TChiNuData : public TObject {
       fChiNu_Time.push_back(Time);
     };//!
 
+    // Particle Type
+    inline void SetIsGamma(const UShort_t& DetNbr,const Bool_t& Type)	{
+      fChiNu_IG_DetectorNbr.push_back(DetNbr);     
+      fChiNu_IsGamma.push_back(Type);
+    };//!
+
 
     //////////////////////    GETTERS    ////////////////////////
     // Energy
@@ -95,6 +105,14 @@ class TChiNuData : public TObject {
     inline Double_t Get_Time(const unsigned int &i) const 
       {return fChiNu_Time[i];}//!
 
+    // Particle Type
+    inline UShort_t GetMultIsGamma() const
+      {return fChiNu_IG_DetectorNbr.size();}
+    inline UShort_t GetIG_DetectorNbr(const unsigned int &i) const 
+      {return fChiNu_IG_DetectorNbr[i];}//!
+    inline Double_t Get_IsGamma(const unsigned int &i) const 
+      {return fChiNu_IsGamma[i];}//!
+
 
   //////////////////////////////////////////////////////////////
   // Required for ROOT dictionnary

NPLib/Detectors/ChiNu/TChiNuPhysics.cxx

@@ -92,6 +92,7 @@ void TChiNuPhysics::BuildPhysicalEvent() {
         DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e));
         Energy.push_back(m_PreTreatedData->Get_Energy(e));
         Time.push_back(m_PreTreatedData->Get_Time(t));
+        IsGamma.push_back(m_PreTreatedData->Get_IsGamma(e));
       }
     }
   }
@@ -125,6 +126,14 @@ void TChiNuPhysics::PreTreat() {
     Double_t Time= Cal->ApplyCalibration("ChiNu/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i));
     m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time);
   }
+
+  // IsGamma 
+  mysize = m_EventData->GetMultIsGamma();
+  for (UShort_t i = 0; i < mysize; ++i) {
+    Bool_t IsGamma= m_EventData->Get_IsGamma(i);
+    m_PreTreatedData->SetIsGamma(m_EventData->GetIG_DetectorNbr(i), IsGamma);
+  }
+ 
 }
 
 
@@ -198,6 +207,7 @@ void TChiNuPhysics::Clear() {
   DetectorNumber.clear();
   Energy.clear();
   Time.clear();
+  IsGamma.clear();
 }
 
 

NPLib/Detectors/ChiNu/TChiNuPhysics.h

@@ -65,6 +65,7 @@ class TChiNuPhysics : public TObject, public NPL::VDetector {
     vector<int>      DetectorNumber;
     vector<double>   Energy;
     vector<double>   Time;
+    vector<bool>     IsGamma;
 
   /// A usefull method to bundle all operation to add a detector
   void AddDetector(TVector3 POS); 

NPSimulation/Detectors/ChiNu/ChiNu.cc

@@ -23,6 +23,7 @@
 #include <sstream>
 #include <cmath>
 #include <limits>
+#include <regex>
 //G4 Geometry object
 #include "G4Tubs.hh"
 #include "G4Box.hh"
@@ -45,6 +46,7 @@
 // NPTool header
 #include "ChiNu.hh"
 #include "CalorimeterScorers.hh"
+#include "ProcessScorers.hh"
 #include "InteractionScorers.hh"
 #include "RootOutput.h"
 #include "MaterialManager.hh"
@@ -318,31 +320,89 @@ void ChiNu::InitializeRootOutput(){
 void ChiNu::ReadSensitive(const G4Event* ){
   m_Event->Clear();
 
+
+  ///////////
+  // Process scorer
+  ProcessScorers::PS_Process* Process_scorer = (ProcessScorers::PS_Process*) m_ChiNuScorer->GetPrimitive(2);
+  unsigned int ProcessMult = Process_scorer->GetMult();
+  vector<int> IsGamma;
+
+  int m_VID=-1;
+  int g_hit=0;
+
+  for(unsigned int i = 0 ; i < ProcessMult ; i++){
+    string particle_name = Process_scorer->GetParticleName(i);
+    string volume_name   = Process_scorer->GetVolumeName(i);
+    string process_name   = Process_scorer->GetProcessName(i);
+
+
+    regex pattern("impr_(\\d+)_logic");
+    smatch match;
+    int VolumeNbr;
+
+    if(regex_search( volume_name,match,pattern))
+      VolumeNbr = stoi(match[1]);
+    else VolumeNbr = -1;
+
+    if(VolumeNbr!=m_VID && particle_name=="gamma")
+      {// if the first particle in the volume is a gamma, it will give a wrong time
+	// cout << "!!!!!!!! First particle is gamma !!!!!!!!!!" << endl;
+	// cout << "i=" << i << " particle_name=" << particle_name << endl;
+	// cout << "volume_name=" << volume_name << endl;
+	// cout << "process_name=" << process_name << endl;
+	g_hit=1;
+      }
+    else if(VolumeNbr!=m_VID)
+      g_hit=0;
+    else if(g_hit==1){
+      // Only if an electron ionizes the liquid scintillator will there be a measured energy deposit and associated time.
+      // cout << "//////////" << endl;
+      // cout << "i=" << i << " particle_name=" << particle_name << endl;
+      // cout << "volume_name=" << volume_name << endl;
+      // cout << "process_name=" << process_name << endl;
+      if(particle_name=="e-")
+	{
+	  g_hit++;
+	  IsGamma.push_back(m_VID);
+	}
+    }
+
+    m_VID=VolumeNbr;
+  }
+
   ///////////
   // Calorimeter scorer
   CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_ChiNuScorer->GetPrimitive(0);
 
-  unsigned int size = Scorer->GetMult(); 
+  unsigned int size = Scorer->GetMult();
+  vector<int> DetID;
   for(unsigned int i = 0 ; i < size ; i++){
     vector<unsigned int> level = Scorer->GetLevel(i);
     int DetectorNbr = level[0];
-    
+    DetID.push_back(DetectorNbr);
     // Light output equation:
     // L(Ep)=A*Ep - B*(1-exp(-C*pow(Ep,D))) from F. Pino et al. Applied Radiation and Isotopes 89 (2014) 79-84.
     double A=0.62, B=1.3, C=0.39, D=0.97;
     double Ep=RandGauss::shoot(Scorer->GetEnergy(i),ChiNu_NS::ResoEnergyOffset+ChiNu_NS::ResoEnergySlope*Scorer->GetEnergy(i));
     double Energy = A*Ep - B*(1-exp(-C*pow(Ep,D)));
     ////////////////////////
-    
-    //    Energy = RandGauss::shoot(Energy,ChiNu_NS::ResoEnergyOffset+ChiNu_NS::ResoEnergySlope*Energy);
+ 
     
     if(Energy>ChiNu_NS::EnergyThreshold){
       //      double Time = RandGauss::shoot(Scorer->GetTime(i),ChiNu_NS::ResoTime); // One time resolution for all detectors
       double Time = RandGauss::shoot(Scorer->GetTime(i),m_ResoTime.at(DetectorNbr-1)); // Time resolution that can be set different for each detector
       m_Event->SetEnergy(DetectorNbr,Energy);
-      m_Event->SetTime(DetectorNbr,Time); 
+      m_Event->SetTime(DetectorNbr,Time);
+
+      // Check if the detector was first lighted by a gamma ray, which would have given a wrong time.
+      auto it = find(IsGamma.begin(), IsGamma.end(), DetectorNbr);
+      if(it != IsGamma.end())
+	m_Event->SetIsGamma(DetectorNbr,true);
+      else
+	m_Event->SetIsGamma(DetectorNbr,false);
     }
   }
+  IsGamma.clear();
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -359,9 +419,11 @@ void ChiNu::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_ChiNuScorer->RegisterPrimitive(Calorimeter);
   m_ChiNuScorer->RegisterPrimitive(Interaction);
+  m_ChiNuScorer->RegisterPrimitive(Process);
   G4SDManager::GetSDMpointer()->AddNewDetector(m_ChiNuScorer) ;
 }
 

NPSimulation/Scorers/ProcessScorers.cc

@@ -29,6 +29,7 @@ PS_Process::PS_Process(G4String name, G4int depth) : G4VPrimitiveScorer(name, de
   // m_NestingLevel = NestingLevel;
   auto tree = RootOutput::getInstance()->GetTree();
   tree->Branch("ProcScor_Particle", &t_particlename);
+  tree->Branch("ProcScor_Volume", &t_volumename);
   tree->Branch("ProcScor_Process", &t_processname);
   tree->Branch("ProcScor_Time", &t_processtime);
   tree->Branch("ProcScor_Efficiency", &t_efficiency);
@@ -71,6 +72,7 @@ G4bool PS_Process::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
     //step_length = aStep->GetTrack()->GetStepLength();
     track_kineE = aStep->GetTrack()->GetKineticEnergy();
     t_particlename.push_back(particlename);
+    t_volumename.push_back(volumename);
     t_kinetic_energy.push_back(track_kineE);
 
     //Looking after the neutrons 
@@ -143,6 +145,7 @@ void PS_Process::EndOfEvent(G4HCofThisEvent*) {}
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 void PS_Process::clear() {
   t_particlename.clear();
+  t_volumename.clear();
   t_processname.clear();
   t_processtime.clear();
   t_gamma_energy.clear();

NPSimulation/Scorers/ProcessScorers.hh

@@ -59,6 +59,7 @@ namespace ProcessScorers {
       //double t_Time;
       //vector<unsigned int> t_Level;
       vector<string> t_particlename;
+      vector<string> t_volumename;
       vector<string> t_processname;
       vector<double> t_processtime;
       vector<double> t_gamma_energy;
@@ -74,6 +75,7 @@ namespace ProcessScorers {
       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 string GetVolumeName(const unsigned int& i) {return t_volumename[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;};