Getting and saving photon energies

[hship2]

Hello! I am modifying example B4a to calculate the energy of photons that pass through a piece of aerogel and then touch a tally surface.

I already have the tally surface counting the number of photons, and I just need to add in a function such as GetKineticEnergy to save each particle’s energy/momentum.

However, since the built-in functionality is counting up energy deposition and step length by just summing each quantity in aggregate, I don’t know how to make sure the kinetic energy stays as a separate quantity. In the end, I want to have a Root histogram that shows how many photons were collected at each energy quanta.

I’m happy to share code snippets if that would be helpful, but so far I haven’t modified the functions much beyond the usual B4a example. I appreciate any help. Thanks!

Geant4 Version: 11.2.1
Operating System: MacOS
CMake Version: 3.29.4

---

[svanlang]

I do something similar in my simulations using a vector. If you’re interested I list some modifications you can try.

ActionInitialization:

  EventAction* eventAction = new EventAction();
  SetUserAction(eventAction);
  RunAction* runAction = new RunAction(eventAction);
  SetUserAction(runAction);

RunAction:

RunAction::RunAction(EventAction* ev) : fEvAction(ev){ ... }

analysisManager->CreateNtupleDColumn("EnergySpectrum", fEvAction->GetEmissionSpectrum());

EventAction:

#include <vector>
std::vector<G4double> vES;
std::vector<G4double>&     GetEmissionSpectrum(){return vES;}
Add2ES(G4double E){vES.push_back(E);}

SteppingAction:

if (...) {
  fEvAction->Add2ES(step->GetTrack()->GetKineticEnergy()/keV);
}

[maire1]

I do not know exactly how you count photons which touch the tally (please post you code, if possible).
At the same place, you can add something like :
G4double ekin = step->GetTrack()->GetKineticEnergy();
and fill a 1D histogram with ekin.

[hship2]

Filling the histogram with ekin worked, thanks for the help!

I tried implementing the vector solution, but it was giving me an error that I hadn’t initialized RunAction properly, that it couldn’t convert EventAction* to RunAction type, and that RunAction wasn’t supposed to take an argument. I think as a result of those errors it was also complaining about the

RunAction::RunAction(EventAction* ev) : fEvAction(ev){ ... }

line, because when I put that in the RunAction file the compiler stopped recognizing analysisManager and quit when more than 20 errors had been generated.

Anyway, now in SteppingAction I have

// photon counter
    auto photcheck = step->GetTrack()->GetDefinition();
    auto photcount = 0.;
    
    if ( photcheck == G4OpticalPhoton::OpticalPhotonDefinition() ) {
         photcount = photcount + 1.;
         G4double ekin = step->GetTrack()->GetKineticEnergy();
        auto analysisManager = G4AnalysisManager::Instance();
        analysisManager->FillH1(7,ekin);
        analysisManager->FillNtupleDColumn(7,ekin);
...
if ( volume == fDetConstruction->GetStepPV() ) {
      fEventAction->AddStep(edep,stepLength, photcount);

This counts the photons and saves their energies in a histogram (Step is the name of the tally volume).
For any other newbies out there, the function AddStep is defined in EventAction as

public:
    void AddStep(G4double de, G4double dl, G4double dc);
    
private:
    G4double  fEnergyStep = 0.;
    G4double  fTrackLStep = 0.;
    G4double  fPCountStep = 0.;
...
inline void EventAction::AddStep(G4double de, G4double dl, G4double dc) {
    fEnergyStep += de;
    fTrackLStep += dl;
    fPCountStep += dc;
}