diff --git a/bind/G4Calo.py b/bind/G4Calo.py index 09a8951..eb93c4d 100644 --- a/bind/G4Calo.py +++ b/bind/G4Calo.py @@ -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, diff --git a/src/RunAction.cc b/src/RunAction.cc index 160f566..b19eed9 100644 --- a/src/RunAction.cc +++ b/src/RunAction.cc @@ -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(); } diff --git a/src/SteppingAction.cc b/src/SteppingAction.cc index ec80e5d..1f84d52 100644 --- a/src/SteppingAction.cc +++ b/src/SteppingAction.cc @@ -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......