minicalosim/src/SteppingAction.cc

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
/// \file B4/B4a/src/SteppingAction.cc
/// \brief Implementation of the B4a::SteppingAction class

#include "SteppingAction.hh"
#include "EventAction.hh"
#include "DetectorConstruction.hh"

#include "G4Step.hh"
#include "G4RunManager.hh"
#include "G4AnalysisManager.hh"
#include "G4Track.hh"
#include "G4ParticleDefinition.hh"
#include "G4TransportationManager.hh"
#include "G4FieldManager.hh"
#include "G4Field.hh"
#include "G4SystemOfUnits.hh"

using namespace B4;

namespace B4a
{

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

SteppingAction::SteppingAction(EventAction* eventAction)
  : fEventAction(eventAction)
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void SteppingAction::UserSteppingAction(const G4Step* step)
{
// Collect energy and track length step by step

  // get volume of the current step
  auto volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();

  // energy deposit
  auto edep = step->GetTotalEnergyDeposit();

  //return;
  
  //find the layer through the GeometryDescriptor that corresponds to the volume
  auto cw = DetectorConstruction::getDetectorConstruction()->getGeometryDescriptor();
  if(cw == nullptr){
    G4cout << "GeometryDescriptor not found" << G4endl;
    throw std::runtime_error("GeometryDescriptor not found");
  }

  
  auto sensor = DetectorConstruction::getDetectorConstruction()->getGeometryDescriptor()->getSensorByVolume(volume);
  if(sensor != nullptr){
    //G4cout << "Sensor found in " << volume->GetName() << " with copy number " << volume->GetCopyNo()<< G4endl; //DEBUG
    sensor->energy += edep;
  }

  // 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();

  // Write one Spawning row per secondary born in this step (track IDs pre-assigned)
  const auto* secondaries = step->GetSecondaryInCurrentStep();
  if (secondaries) {
    for (const G4Track* sec : *secondaries) {
      analysisManager->FillNtupleIColumn(2, 0, evtId);
      analysisManager->FillNtupleIColumn(2, 1, track->GetTrackID());
      analysisManager->FillNtupleIColumn(2, 2, track->GetCurrentStepNumber());
      analysisManager->FillNtupleIColumn(2, 3, sec->GetTrackID());
      analysisManager->AddNtupleRow(2);
    }
  }

  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());

  // A: energy deposited in medium at this step (≠ pre_E - post_E when secondaries are created)
  analysisManager->FillNtupleDColumn(1, 12, edep);

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

  // E: layer index (-1 if outside all layers) and material name at pre-step point
  int layerId = -1;
  const auto& layers = cw->getLayers();
  for (int i = 0; i < (int)layers.size(); i++) {
    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);

  // 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());

  // 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;
  const auto* fm = G4TransportationManager::GetTransportationManager()->GetFieldManager();
  if (fm && fm->DoesFieldExist()) {
    const G4Field* field = fm->GetDetectorField();
    if (field) {
      const G4double point[4] = {pre->GetPosition().x(), pre->GetPosition().y(),
                                  pre->GetPosition().z(), pre->GetGlobalTime()};
      G4double fieldVal[6] = {0,0,0,0,0,0};
      field->GetFieldValue(point, fieldVal);
      Bx = fieldVal[0] / tesla;
      By = fieldVal[1] / tesla;
      Bz = fieldVal[2] / tesla;
      Ex = fieldVal[3] / (volt/m);
      Ey = fieldVal[4] / (volt/m);
      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->AddNtupleRow(1);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

}