diff --git a/include/RunAction.hh b/include/RunAction.hh index 854bb55..20a1e8d 100644 --- a/include/RunAction.hh +++ b/include/RunAction.hh @@ -83,6 +83,8 @@ class RunAction : public G4UserRunAction mutable std::vector hitLayer; mutable std::vector hitCopyNumber; + mutable std::vector stepChildIds; + mutable G4String filename; }; diff --git a/include/SteppingAction.hh b/include/SteppingAction.hh index 85cdd0f..6cf86cb 100644 --- a/include/SteppingAction.hh +++ b/include/SteppingAction.hh @@ -33,6 +33,7 @@ #include "G4UserSteppingAction.hh" #include "GeometryDescriptor.hh" #include "DetectorConstruction.hh" +#include "RunAction.hh" namespace B4a @@ -54,8 +55,11 @@ public: void UserSteppingAction(const G4Step* step) override; + void setRunAction(const B4::RunAction* runAction) { fRunAction = runAction; } + private: EventAction* fEventAction = nullptr; + const B4::RunAction* fRunAction = nullptr; }; } diff --git a/src/ActionInitialization.cc b/src/ActionInitialization.cc index d862832..5383cea 100644 --- a/src/ActionInitialization.cc +++ b/src/ActionInitialization.cc @@ -70,6 +70,7 @@ void ActionInitialization::Build() const eventAction->setPrimaryGeneratorAction(gen); SetUserAction(eventAction); auto steppingAction = new SteppingAction(eventAction); + steppingAction->setRunAction(runact); SetUserAction(steppingAction); } diff --git a/src/RunAction.cc b/src/RunAction.cc index d9d9afa..2d7507f 100644 --- a/src/RunAction.cc +++ b/src/RunAction.cc @@ -98,19 +98,19 @@ RunAction::RunAction() analysisManager->CreateNtupleDColumn("post_E"); // col 11 analysisManager->CreateNtupleDColumn("edep"); // col 12 [MeV] analysisManager->CreateNtupleDColumn("step_length"); // col 13 [mm] - analysisManager->CreateNtupleIColumn("parent_id"); // col 14 - analysisManager->CreateNtupleSColumn("process"); // col 15 - analysisManager->CreateNtupleIColumn("layer_id"); // col 16 (-1 = outside all layers) - analysisManager->CreateNtupleSColumn("material"); // col 17 - analysisManager->CreateNtupleDColumn("pre_dx"); // col 18 - analysisManager->CreateNtupleDColumn("pre_dy"); // col 19 - analysisManager->CreateNtupleDColumn("pre_dz"); // col 20 - analysisManager->CreateNtupleDColumn("Bx"); // col 21 [T] - analysisManager->CreateNtupleDColumn("By"); // col 22 [T] - analysisManager->CreateNtupleDColumn("Bz"); // col 23 [T] - analysisManager->CreateNtupleDColumn("Ex"); // col 24 [V/m] - analysisManager->CreateNtupleDColumn("Ey"); // col 25 [V/m] - analysisManager->CreateNtupleDColumn("Ez"); // col 26 [V/m] + analysisManager->CreateNtupleSColumn("process"); // col 14 + analysisManager->CreateNtupleIColumn("layer_id"); // col 15 (-1 = outside all layers) + analysisManager->CreateNtupleSColumn("material"); // col 16 + analysisManager->CreateNtupleDColumn("pre_dx"); // col 17 + analysisManager->CreateNtupleDColumn("pre_dy"); // col 18 + analysisManager->CreateNtupleDColumn("pre_dz"); // col 19 + analysisManager->CreateNtupleDColumn("Bx"); // col 20 [T] + analysisManager->CreateNtupleDColumn("By"); // col 21 [T] + analysisManager->CreateNtupleDColumn("Bz"); // col 22 [T] + analysisManager->CreateNtupleDColumn("Ex"); // col 23 [V/m] + analysisManager->CreateNtupleDColumn("Ey"); // col 24 [V/m] + analysisManager->CreateNtupleDColumn("Ez"); // col 25 [V/m] + analysisManager->CreateNtupleIColumn("child_track_ids", stepChildIds); // col 26 (variable-length array) analysisManager->FinishNtuple(); // Spawning ntuple (one row per secondary born, ntuple id=2) diff --git a/src/SteppingAction.cc b/src/SteppingAction.cc index 4a6ac67..cb4146f 100644 --- a/src/SteppingAction.cc +++ b/src/SteppingAction.cc @@ -87,10 +87,13 @@ void SteppingAction::UserSteppingAction(const G4Step* step) auto track = step->GetTrack(); int evtId = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID(); - // Write one Spawning row per secondary born in this step (track IDs pre-assigned) + // Write one Spawning row per secondary born in this step (track IDs pre-assigned); + // also populate stepChildIds for the child_track_ids vector column in Steps. + fRunAction->stepChildIds.clear(); const auto* secondaries = step->GetSecondaryInCurrentStep(); if (secondaries) { for (const G4Track* sec : *secondaries) { + fRunAction->stepChildIds.push_back(sec->GetTrackID()); analysisManager->FillNtupleIColumn(2, 0, evtId); analysisManager->FillNtupleIColumn(2, 1, track->GetTrackID()); analysisManager->FillNtupleIColumn(2, 2, track->GetCurrentStepNumber()); @@ -118,14 +121,11 @@ void SteppingAction::UserSteppingAction(const G4Step* step) // B: geometric step length analysisManager->FillNtupleDColumn(1, 13, step->GetStepLength()); - // C: parent track ID (0 for primary) - analysisManager->FillNtupleIColumn(1, 14, track->GetParentID()); - // D: process that ended this step G4String processName = ""; const auto* postProc = post->GetProcessDefinedStep(); if (postProc) processName = postProc->GetProcessName(); - analysisManager->FillNtupleSColumn(1, 15, processName); + analysisManager->FillNtupleSColumn(1, 14, processName); // E: layer index (-1 if outside all layers) and material name at pre-step point int layerId = -1; @@ -134,14 +134,14 @@ void SteppingAction::UserSteppingAction(const G4Step* step) if (layers[i].physicalVolume == volume) { layerId = i; break; } } G4String materialName = pre->GetMaterial() ? pre->GetMaterial()->GetName() : ""; - analysisManager->FillNtupleIColumn(1, 16, layerId); - analysisManager->FillNtupleSColumn(1, 17, materialName); + analysisManager->FillNtupleIColumn(1, 15, layerId); + analysisManager->FillNtupleSColumn(1, 16, materialName); // F: pre-step momentum direction (unit vector) const auto dir = pre->GetMomentumDirection(); - analysisManager->FillNtupleDColumn(1, 18, dir.x()); - analysisManager->FillNtupleDColumn(1, 19, dir.y()); - analysisManager->FillNtupleDColumn(1, 20, dir.z()); + analysisManager->FillNtupleDColumn(1, 17, dir.x()); + analysisManager->FillNtupleDColumn(1, 18, dir.y()); + analysisManager->FillNtupleDColumn(1, 19, dir.z()); // Field: B [T] and E [V/m] at pre-step position; zero if no field is registered G4double Bx=0, By=0, Bz=0, Ex=0, Ey=0, Ez=0; @@ -161,12 +161,13 @@ void SteppingAction::UserSteppingAction(const G4Step* step) Ez = fieldVal[5] / (volt/m); } } - analysisManager->FillNtupleDColumn(1, 21, Bx); - analysisManager->FillNtupleDColumn(1, 22, By); - analysisManager->FillNtupleDColumn(1, 23, Bz); - analysisManager->FillNtupleDColumn(1, 24, Ex); - analysisManager->FillNtupleDColumn(1, 25, Ey); - analysisManager->FillNtupleDColumn(1, 26, Ez); + analysisManager->FillNtupleDColumn(1, 20, Bx); + analysisManager->FillNtupleDColumn(1, 21, By); + analysisManager->FillNtupleDColumn(1, 22, Bz); + analysisManager->FillNtupleDColumn(1, 23, Ex); + analysisManager->FillNtupleDColumn(1, 24, Ey); + analysisManager->FillNtupleDColumn(1, 25, Ez); + // col 26 child_track_ids: vector column, auto-read from fRunAction->stepChildIds at AddNtupleRow analysisManager->AddNtupleRow(1); }