Merge pull request #27 from lopezzot/lp-addleakagecounter

Add (leakage) particle spectrum analyzer

CMakeLists.txt

@@ -22,6 +22,16 @@ else()
     find_package(Geant4 REQUIRED)
 endif()
 
+#Option to enable leakage particle spectrum analysis
+#
+option(WITH_LEAKAGEANALYSIS "enable SpectrumAnalysis on leakage" OFF)
+if(WITH_LEAKAGEANALYSIS)
+    add_compile_definitions(ATLHECTB_LEAKANALYSIS)
+    message(STATUS "WITH_LEAKAGEANALYSIS = ON : building with ATLHECTB_LEAKANALYSIS compiler definition.")
+else()
+    message(STATUS "WITH_LEAKAGEANALYSIS = OFF : building without ATLHECTB_LEAKANALYSIS compiler definition.")
+endif()
+
 #Setup Geant4 include directories and compile definitions
 #
 include(${Geant4_USE_FILE})

README.md

@@ -234,6 +234,9 @@ Parser options
    * -t integer: pass number of threads for multi-thread execution (example -t 3, default t=2)
    * -pl Physics_List: select Geant4 physics list (example -pl FTFP_BERT)
    * It is possible to select alternative FTF tunings with PL_tuneID (example -pl FTFP_BERT_tune0) [only for Geant4-11.1.0 or higher]
+
+Custom CMake options
+   * `WITH_LEAKAGEANALYSIS`: if set to `ON` include the usage of the `SpectrumAnalysis` singleton to study the particle leakage (default `OFF`)
 5. (optional) It is possible to install the executable in `bin` under `CMAKE_INSTALL_PREFIX`
    ```sh
    cmake -DCMAKE_INSTALL_PREFIX=/absolute-path-to-installdir/ -DGeant4_DIR=/absolute_path_to/geant4.10.07_p01-install/lib/Geant4-10.7.1/ relative_path_to/ATLHECTB/

include/SpectrumAnalyzer.hh

@@ -0,0 +1,89 @@
+//**************************************************
+// \file SpectrumAnalyzer.hh
+// \brief: Declaration of SpectrumAnalyzer class
+// \author: Lorenzo Pezzotti (CERN EP-SFT-sim)
+//          @lopezzot
+// \start date: 28 August 2023
+//**************************************************
+
+// A portable Geant4-based particle spectrum analyzer
+// to be used within a Geant4 simulation without affecting it.
+// Instead of coding it in the simulation, create a singleton
+// and manage its usage with (#ifdef) compiler definition.
+
+#ifndef SpectrumAnalyzer_h
+#  define SpectrumAnalyzer_h
+
+// Includers from Geant4
+//
+#  include "G4Step.hh"
+#  include "G4ThreadLocalSingleton.hh"
+
+// Includers from C++
+//
+#  include <functional>
+
+class SpectrumAnalyzer
+{
+    friend class G4ThreadLocalSingleton<SpectrumAnalyzer>;
+
+  public:
+    // Return pointer to class instance
+    static SpectrumAnalyzer* GetInstance()
+    {
+      static G4ThreadLocalSingleton<SpectrumAnalyzer> instance{};
+      return instance.Instance();
+    }
+
+    // Methods
+    //
+    // Run-wise methods
+    void CreateNtupleAndScorer(const G4String scName = "te");
+    inline void ClearNtupleID() { ntupleID = 99; }
+    // Event-wise methods
+    inline void ClearEventFields()
+    {
+      neutronScore = 0.;
+      protonScore = 0., pionScore = 0., gammaScore = 0., electronScore = 0.;
+    }
+    void FillEventFields() const;
+    // Step-wise methods
+    void Analyze(const G4Step* step);
+
+  private:
+    // Members
+    //
+    // Run-wise members
+    G4int ntupleID;
+    std::function<G4double(const G4Step* step)> scorer;
+    G4String scorerName{};
+    // Event-wise members
+    G4double neutronScore;
+    G4double protonScore;
+    G4double pionScore;
+    G4double gammaScore;
+    G4double electronScore;
+
+    // Scoring quantities
+    inline static G4double GetMomentum(const G4Step* step)
+    {
+      return step->GetTrack()->GetMomentum().mag();
+    };
+    inline static G4double GetKE(const G4Step* step)
+    {
+      return step->GetTrack()->GetKineticEnergy();
+    };
+    inline static G4double GetTE(const G4Step* step) { return step->GetTrack()->GetTotalEnergy(); };
+
+  private:
+    // Private constructor
+    SpectrumAnalyzer() = default;
+
+  public:
+    SpectrumAnalyzer(SpectrumAnalyzer const&) = delete;
+    void operator=(SpectrumAnalyzer const&) = delete;
+};
+
+#endif  // SpectrumAnalyzer_h
+
+//**************************************************

src/ATLHECTBEventAction.cc

@@ -10,6 +10,7 @@
 #include "ATLHECTBEventAction.hh"
 
 #include "ATLHECTBRunAction.hh"
+#include "SpectrumAnalyzer.hh"
 
 // Includers from Geant4
 //
@@ -114,6 +115,10 @@ void ATLHECTBEventAction::BeginOfEventAction(const G4Event*)
   for (unsigned int i = 0; i < 20; i++) {
     M3L4BirkeLayer.push_back(0.);
   }
+
+#ifdef ATLHECTB_LEAKANALYSIS
+  SpectrumAnalyzer::GetInstance()->ClearEventFields();
+#endif
 }
 
 void ATLHECTBEventAction::EndOfEventAction(const G4Event*)
@@ -133,6 +138,10 @@ void ATLHECTBEventAction::EndOfEventAction(const G4Event*)
   analysisManager->FillNtupleDColumn(5, elAr);
   analysisManager->FillNtupleDColumn(6, BirkelAr);
   analysisManager->AddNtupleRow();
+
+#ifdef ATLHECTB_LEAKANALYSIS
+  SpectrumAnalyzer::GetInstance()->FillEventFields();
+#endif
 }
 
 //**************************************************

src/ATLHECTBRunAction.cc

@@ -10,6 +10,7 @@
 #include "ATLHECTBRunAction.hh"
 
 #include "ATLHECTBEventAction.hh"
+#include "SpectrumAnalyzer.hh"
 
 // Includers from Geant4
 //
@@ -59,6 +60,10 @@ ATLHECTBRunAction::ATLHECTBRunAction(ATLHECTBEventAction* eventAction)
   analysisManager->CreateNtupleDColumn("M2L4BirkeLayer", fEventAction->GetM2L4BirkeLayer());
   analysisManager->CreateNtupleDColumn("M3L4BirkeLayer", fEventAction->GetM3L4BirkeLayer());
   analysisManager->FinishNtuple();
+
+#ifdef ATLHECTB_LEAKANALYSIS
+  SpectrumAnalyzer::GetInstance()->CreateNtupleAndScorer("ke");
+#endif
 }
 
 // Define deconstructor
@@ -91,6 +96,10 @@ void ATLHECTBRunAction::EndOfRunAction(const G4Run*)
 
   analysisManager->Write();
   analysisManager->CloseFile();
+
+#ifdef ATLHECTB_LEAKANALYSIS
+  SpectrumAnalyzer::GetInstance()->ClearNtupleID();
+#endif
 }
 
 //**************************************************

src/ATLHECTBSteppingAction.cc

@@ -11,6 +11,7 @@
 
 #include "ATLHECTBDetectorConstruction.hh"
 #include "ATLHECTBEventAction.hh"
+#include "SpectrumAnalyzer.hh"
 
 // Includers from Geant4
 //
@@ -49,6 +50,9 @@ void ATLHECTBSteppingAction::UserSteppingAction(const G4Step* step)
   //
   if (!step->GetTrack()->GetNextVolume()) {
     fEventAction->Addeleak(step->GetTrack()->GetKineticEnergy());
+#ifdef ATLHECTB_LEAKANALYSIS
+    SpectrumAnalyzer::GetInstance()->Analyze(step);
+#endif
   }
 
   // Collect energy deposited in test beam prototype

src/SpectrumAnalyzer.cc

@@ -0,0 +1,94 @@
+//**************************************************
+// \file SpectrumAnalyzer.cc
+// \brief: Definition of SpectrumAnalyzer class
+// \author: Lorenzo Pezzotti (CERN EP-SFT-sim)
+//          @lopezzot
+// \start date: 28 August 2023
+//**************************************************
+
+// Includers from project files
+//
+#include "SpectrumAnalyzer.hh"
+
+// Includers from Geant4
+//
+#include "G4Version.hh"
+#if G4VERSION_NUMBER < 1100
+#  include "g4root.hh"
+#else
+#  include "G4AnalysisManager.hh"
+#endif
+#include "G4Step.hh"
+
+// #define DEBUG
+
+void SpectrumAnalyzer::CreateNtupleAndScorer(const G4String scName)
+{
+  auto AM = G4AnalysisManager::Instance();
+
+  ntupleID = AM->CreateNtuple("Spectrum", "Spectrum");
+  AM->CreateNtupleDColumn("neutronScore");
+  AM->CreateNtupleDColumn("protonScore");
+  AM->CreateNtupleDColumn("pionScore");
+  AM->CreateNtupleDColumn("gammaScore");
+  AM->CreateNtupleDColumn("electronScore");
+  AM->FinishNtuple();
+
+  // Define scorer type
+  scorerName = scName;
+  if (scorerName == "te") {
+    scorer = GetTE;
+  }
+  if (scorerName == "momentum") {
+    scorer = GetMomentum;
+  }
+  else if (scorerName == "ke") {
+    scorer = GetKE;
+  }
+  else {
+    scorer = GetTE;
+  }  // default case
+}
+
+void SpectrumAnalyzer::FillEventFields() const
+{
+  auto AM = G4AnalysisManager::Instance();
+  AM->FillNtupleDColumn(ntupleID, 0, neutronScore);
+  AM->FillNtupleDColumn(ntupleID, 1, protonScore);
+  AM->FillNtupleDColumn(ntupleID, 2, pionScore);
+  AM->FillNtupleDColumn(ntupleID, 3, gammaScore);
+  AM->FillNtupleDColumn(ntupleID, 4, electronScore);
+  AM->AddNtupleRow(ntupleID);
+}
+
+void SpectrumAnalyzer::Analyze(const G4Step* step)
+{
+  auto PDGID = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
+  auto val = scorer(step);
+  if (PDGID == 2112 || PDGID == -2112) {
+    neutronScore += val;
+  }
+  else if (PDGID == 2212 || PDGID == 2212) {
+    protonScore += val;
+  }
+  else if (PDGID == 211 || PDGID == 211) {
+    pionScore += val;
+  }
+  else if (PDGID == 22) {
+    gammaScore += val;
+  }
+  else if (PDGID == -11 || PDGID == 11) {
+    electronScore += val;
+  }
+  else {
+  }
+
+#ifdef DEBUG
+  G4cout << "-->SpectrumAnalyzer::Analyze, scorer name " << scorerName << " " << PDGID << " "
+         << step->GetTrack()->GetParticleDefinition()->GetParticleName() << " Total Energy "
+         << GetTE(step) << " Momentum " << GetMomentum(step) << " Kinetic Energy " << GetKE(step)
+         << G4endl;
+#endif
+}
+
+//**************************************************