Skip to content

Commit

Permalink
Add (leakage) spectrum analyzer
Browse files Browse the repository at this point in the history
Add a portable particle spectrum analyzer singleton that can be used
within the simulation without affecting it. Its usage is managed by a
dedicated compiler definition. Used to study leaking particles, but can
be used for different studies as well.
  • Loading branch information
lopezzot committed Aug 28, 2023
1 parent 0aa1923 commit 87daaa3
Show file tree
Hide file tree
Showing 7 changed files with 228 additions and 5 deletions.
11 changes: 11 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,17 @@ if(WITH_ATLTileCalTB_PulseOutput)
add_compile_definitions(ATLTileCalTB_PulseOutput)
endif()

#----------------------------------------------------------------------------
#Option to enable leakage particle spectrum analysis
#
option(WITH_LEAKAGEANALYSIS "enable SpectrumAnalysis on leakage" OFF)
if(WITH_LEAKAGEANALYSIS)
add_compile_definitions(ATLTileCalTB_LEAKANALYSIS)
message(STATUS "WITH_LEAKAGEANALYSIS = ON : building with ATLTileCalTB_LEAKANALYSIS compiler definition.")
else()
message(STATUS "WITH_LEAKAGEANALYSIS = OFF : building without ATLTileCalTB_LEAKANALYSIS compiler definition.")
endif()

#----------------------------------------------------------------------------
# Option to disable noise
#
Expand Down
3 changes: 2 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -287,7 +287,8 @@ Custom options:
- `WITH_ATLTileCalTB_NoNoise`: if set to `ON`, the simulation will not put electronic noise on the
signal (per cell) and disable the 2 sigma noise cut. Only relevant for noise calibration.
- `WITH_GEANT4_UIVIS`: if set to `ON` (default), build with UI and visualization drivers.
- `G4_USE_FLUKA`: if set to `ON` build against the Fluka.Cern interface (default `OFF`).
- `G4_USE_FLUKA`: if set to `ON` build against the Fluka.Cern interface (default `OFF`).
- `WITH_LEAKAGEANALYSIS`: if set to `ON` build with leakage spectrum analyzer (default `OFF`).

Relevant built-in options:
- `CMAKE_BUILD_TYPE`: set to `Debug` for debugging and to `Release` for production (faster).
Expand Down
90 changes: 90 additions & 0 deletions include/SpectrumAnalyzer.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
//**************************************************
// \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., othersScore = 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;
G4double othersScore;

// 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

//**************************************************
13 changes: 11 additions & 2 deletions src/ATLTileCalTBEventAction.cc
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,9 @@
#include "ATLTileCalTBGeometry.hh"
#include "ATLTileCalTBConstants.hh"
#include "ATLTileCalTBPrimaryGenAction.hh"
#ifdef ATLTileCalTB_LEAKANALYSIS
#include "SpectrumAnalyzer.hh"
#endif

//Includers from Geant4
//
Expand Down Expand Up @@ -67,7 +70,10 @@ void ATLTileCalTBEventAction::BeginOfEventAction([[maybe_unused]] const G4Event*
pulse_event_path = std::filesystem::path("ATLTileCalTBpulse_Run" + std::to_string(runNumber) + "/Ev" + std::to_string(eventNumber));
std::filesystem::create_directory(pulse_event_path);
#endif


#ifdef ATLTileCalTB_LEAKANALYSIS
SpectrumAnalyzer::GetInstance()->ClearEventFields();
#endif
}

//GetHitsCollection method()
Expand Down Expand Up @@ -200,7 +206,10 @@ void ATLTileCalTBEventAction::EndOfEventAction( const G4Event* event ) {
analysisManager->FillNtupleFColumn(7, fPrimaryGenAction->GetParticlenGun()->GetParticleEnergy());

analysisManager->AddNtupleRow();


#ifdef ATLTileCalTB_LEAKANALYSIS
SpectrumAnalyzer::GetInstance()->FillEventFields();
#endif
}

//**************************************************
13 changes: 11 additions & 2 deletions src/ATLTileCalTBRunAction.cc
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,9 @@
//
#include "ATLTileCalTBRunAction.hh"
#include "ATLTileCalTBEventAction.hh"
#ifdef ATLTileCalTB_LEAKANALYSIS
#include "SpectrumAnalyzer.hh"
#endif

//Includers from Geant4
//
Expand Down Expand Up @@ -73,7 +76,10 @@ ATLTileCalTBRunAction::ATLTileCalTBRunAction( ATLTileCalTBEventAction* eventActi
analysisManager->CreateNtupleIColumn("PDGID");
analysisManager->CreateNtupleFColumn("EBeam");
analysisManager->FinishNtuple();


#ifdef ATLTileCalTB_LEAKANALYSIS
SpectrumAnalyzer::GetInstance()->CreateNtupleAndScorer("ke");
#endif
}

ATLTileCalTBRunAction::~ATLTileCalTBRunAction() {
Expand Down Expand Up @@ -109,7 +115,10 @@ void ATLTileCalTBRunAction::EndOfRunAction(const G4Run* /*run*/) {
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->Write();
analysisManager->CloseFile();


#ifdef ATLTileCalTB_LEAKANALYSIS
SpectrumAnalyzer::GetInstance()->ClearNtupleID();
#endif
}

//**************************************************
6 changes: 6 additions & 0 deletions src/ATLTileCalTBStepAction.cc
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@
//Includers from project files
//
#include "ATLTileCalTBStepAction.hh"
#ifdef ATLTileCalTB_LEAKANALYSIS
#include "SpectrumAnalyzer.hh"
#endif

//Constructor and de-constructor
//
Expand All @@ -27,6 +30,9 @@ void ATLTileCalTBStepAction::UserSteppingAction( const G4Step* aStep ) {
//
if ( !aStep->GetTrack()->GetNextVolume() ){
fEventAction->Add( 0, aStep->GetTrack()->GetKineticEnergy() );
#ifdef ATLTileCalTB_LEAKANALYSIS
SpectrumAnalyzer::GetInstance()->Analyze(aStep);
#endif
}

if ( aStep->GetTrack()->GetTouchableHandle()->GetVolume()->GetName() != "CALO::CALO" ||
Expand Down
97 changes: 97 additions & 0 deletions src/SpectrumAnalyzer.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
//**************************************************
// \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->CreateNtupleDColumn("othersScore");
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->FillNtupleDColumn(ntupleID, 5, othersScore);
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 {
othersScore += val;
}

#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
}

//**************************************************

0 comments on commit 87daaa3

Please sign in to comment.