add per-step kinematics ntuple ("Steps") to ROOT output

Adds a second ntuple alongside the existing "Hits" ntuple with one row
per Geant4 step. Records event_id, track_id, step_no, PDG code, and
pre/post position + kinetic energy for every step in the simulation.

Also extends run_batch() with return_steps=False; when True, returns
(hits_mf, steps_mf) instead of just hits_mf.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
2026-06-08 11:37:39 +02:00
parent 91020714e7
commit 19b1ef125b
3 changed files with 53 additions and 11 deletions
+13 -4
View File
@@ -334,7 +334,8 @@ def run_batch(
maxEnergy_GeV: float = -1.0,
filename: str = "",
no_mp: bool = False,
manual_seed : int = -1):
manual_seed : int = -1,
return_steps: bool = False):
"""
Run a full Geant4 simulation batch with automatic parallelisation.
@@ -361,12 +362,17 @@ no_mp : bool, optional
Disable multiprocessing (run single-threaded) for debugging.
manual_seed : int, optional
Optional random seed. If < 0, seeds are generated automatically.
return_steps : bool, optional
If True, also return a MiniFrame of step-level data (one row per
Geant4 step across all events). Returns a tuple (hits, steps).
Returns
-------
MiniFrame
Table-like structure with one row per event, containing both scalar
values and fixed-size sensor arrays.
MiniFrame or tuple[MiniFrame, MiniFrame]
If return_steps is False (default): event-level MiniFrame with one
row per event. If return_steps is True: (hits, steps) where steps
has one row per Geant4 step with columns event_id, track_id,
step_no, pdg, pre_x/y/z, pre_E, post_x/y/z, post_E.
Notes
-----
@@ -463,6 +469,9 @@ Example
if os.path.exists(f):
os.remove(f)
print('G4Calo: concatenation finished after {:.2f} seconds'.format(sw.elapsed()))
if return_steps:
steps_df = _assemble_results_to_mini_df(rp, tree_name="Steps")
return df, steps_df
return df
def _fill_event(gd : GeometryDescriptor,
+15
View File
@@ -80,7 +80,22 @@ RunAction::RunAction()
analysisManager->CreateNtupleIColumn("sensor_layer",hitLayer);
analysisManager->CreateNtupleIColumn("sensor_copy_number",hitCopyNumber);
analysisManager->FinishNtuple();
// Steps ntuple (one row per step, ntuple id=1)
analysisManager->CreateNtuple("Steps", "Steps");
analysisManager->CreateNtupleIColumn("event_id"); // col 0
analysisManager->CreateNtupleIColumn("track_id"); // col 1
analysisManager->CreateNtupleIColumn("step_no"); // col 2
analysisManager->CreateNtupleIColumn("pdg"); // col 3
analysisManager->CreateNtupleDColumn("pre_x"); // col 4
analysisManager->CreateNtupleDColumn("pre_y"); // col 5
analysisManager->CreateNtupleDColumn("pre_z"); // col 6
analysisManager->CreateNtupleDColumn("pre_E"); // col 7
analysisManager->CreateNtupleDColumn("post_x"); // col 8
analysisManager->CreateNtupleDColumn("post_y"); // col 9
analysisManager->CreateNtupleDColumn("post_z"); // col 10
analysisManager->CreateNtupleDColumn("post_E"); // col 11
analysisManager->FinishNtuple();
}
+25 -7
View File
@@ -33,6 +33,9 @@
#include "G4Step.hh"
#include "G4RunManager.hh"
#include "G4AnalysisManager.hh"
#include "G4Track.hh"
#include "G4ParticleDefinition.hh"
using namespace B4;
@@ -68,16 +71,31 @@ void SteppingAction::UserSteppingAction(const G4Step* step)
auto sensor = DetectorConstruction::getDetectorConstruction()->getGeometryDescriptor()->getSensorByVolume(volume);
if(sensor == nullptr){
//can simply be not an active volume, so no need to throw an error
//G4cout << "Sensor not found in" << volume->GetName() << G4endl; //DEBUG
return;
if(sensor != nullptr){
//G4cout << "Sensor found in " << volume->GetName() << " with copy number " << volume->GetCopyNo()<< G4endl; //DEBUG
sensor->energy += edep;
}
//G4cout << "Sensor found in " << volume->GetName() << " with copy number " << volume->GetCopyNo()<< G4endl; //DEBUG
sensor->energy += edep;
//find sensor by volume copy number
// Record step kinematics for every step
auto analysisManager = G4AnalysisManager::Instance();
auto pre = step->GetPreStepPoint();
auto post = step->GetPostStepPoint();
auto track = step->GetTrack();
int evtId = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
analysisManager->FillNtupleIColumn(1, 0, evtId);
analysisManager->FillNtupleIColumn(1, 1, track->GetTrackID());
analysisManager->FillNtupleIColumn(1, 2, track->GetCurrentStepNumber());
analysisManager->FillNtupleIColumn(1, 3, track->GetDefinition()->GetPDGEncoding());
analysisManager->FillNtupleDColumn(1, 4, pre->GetPosition().x());
analysisManager->FillNtupleDColumn(1, 5, pre->GetPosition().y());
analysisManager->FillNtupleDColumn(1, 6, pre->GetPosition().z());
analysisManager->FillNtupleDColumn(1, 7, pre->GetKineticEnergy());
analysisManager->FillNtupleDColumn(1, 8, post->GetPosition().x());
analysisManager->FillNtupleDColumn(1, 9, post->GetPosition().y());
analysisManager->FillNtupleDColumn(1, 10, post->GetPosition().z());
analysisManager->FillNtupleDColumn(1, 11, post->GetKineticEnergy());
analysisManager->AddNtupleRow(1);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......