G4ScoreSplittingProcess.cc

Go to the documentation of this file.

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
//
//
 
#include "G4ScoreSplittingProcess.hh"
 
#include "G4EnergySplitter.hh"
#include "G4ParticleChange.hh"
#include "G4RegularNavigationHelper.hh"
#include "G4SDManager.hh"
#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "G4SystemOfUnits.hh"
#include "G4TouchableHistory.hh"
#include "G4TransportationManager.hh"
#include "G4VPhysicalVolume.hh"
#include "G4VSensitiveDetector.hh"
#include "G4VTouchable.hh"
#include "G4ios.hh"
 
//--------------------------------
// Constructor with name and type:
//--------------------------------
G4ScoreSplittingProcess::G4ScoreSplittingProcess(const G4String& processName, G4ProcessType theType)
  : G4VProcess(processName, theType)
{
  pParticleChange = &xParticleChange;
 
  fSplitStep = new G4Step();
  fSplitPreStepPoint = fSplitStep->GetPreStepPoint();
  fSplitPostStepPoint = fSplitStep->GetPostStepPoint();
 
  if (verboseLevel > 0) {
    G4cout << GetProcessName() << " is created " << G4endl;
  }
  fpEnergySplitter = new G4EnergySplitter();
}
 
// -----------
// Destructor:
// -----------
G4ScoreSplittingProcess::~G4ScoreSplittingProcess()
{
  delete fSplitStep;
  delete fpEnergySplitter;
}
 
//------------------------------------------------------
//
// StartTracking
//
//------------------------------------------------------
void G4ScoreSplittingProcess::StartTracking(G4Track* trk)
{
  // Setup initial touchables for the first step
  const G4Step* pStep = trk->GetStep();
 
  fOldTouchableH = trk->GetTouchableHandle();
  *fSplitPreStepPoint = *(pStep->GetPreStepPoint());  // Best to copy, so as to initialise
  fSplitPreStepPoint->SetTouchableHandle(fOldTouchableH);
  fNewTouchableH = fOldTouchableH;
  *fSplitPostStepPoint = *(pStep->GetPostStepPoint());  // Best to copy, so as to initialise
  fSplitPostStepPoint->SetTouchableHandle(fNewTouchableH);
 
  ///  Initialize
  fSplitPreStepPoint->SetStepStatus(fUndefined);
  fSplitPostStepPoint->SetStepStatus(fUndefined);
}
 
//----------------------------------------------------------
//
//  PostStepGetPhysicalInteractionLength()
//
//----------------------------------------------------------
G4double G4ScoreSplittingProcess::PostStepGetPhysicalInteractionLength(
  const G4Track& /*track*/, G4double /*previousStepSize*/, G4ForceCondition* condition)
{
  // This process must be invoked anyway to score the hit
  //   - to do the scoring if the current volume is a regular structure, or
  //   - else to toggle the flag so that the SteppingManager does the scoring.
  *condition = StronglyForced;
 
  // Future optimisation: check whether in regular structure.
  //  If it is in regular structure,  be StronglyForced
  //  If  not  in regular structure,
  //         ask to be called only if SteppingControl is AvoidHitInvocation
  //         in order to reset it to NormalCondition
 
  return DBL_MAX;
}
 
//------------------------------------
//
//             PostStepDoIt()
//
//------------------------------------
G4VParticleChange* G4ScoreSplittingProcess::PostStepDoIt(const G4Track& track, const G4Step& step)
{
  G4VPhysicalVolume* pCurrentVolume = track.GetVolume();
  G4LogicalVolume* pLogicalVolume = pCurrentVolume->GetLogicalVolume();
  G4VSensitiveDetector* ptrSD = pLogicalVolume->GetSensitiveDetector();
 
  pParticleChange->Initialize(track);
  if ((!pCurrentVolume->IsRegularStructure()) || (ptrSD == nullptr)
      || G4RegularNavigationHelper::Instance()->GetStepLengths().size() <= 1)
  {
    // Set the flag to make sure that Stepping Manager does the scoring
    pParticleChange->ProposeSteppingControl(NormalCondition);
  }
  else {
    G4ThreeVector preStepPosition, postStepPosition, direction, finalPostStepPosition;
    pParticleChange->ProposeSteppingControl(AvoidHitInvocation);
 
    G4double totalEnergyDeposit = step.GetTotalEnergyDeposit();
    G4StepStatus fullStepStatus = step.GetPostStepPoint()->GetStepStatus();
 
    CopyStepStart(step);
    fSplitPreStepPoint->SetSensitiveDetector(ptrSD);
    fOldTouchableH = fInitialTouchableH;
    fNewTouchableH = fOldTouchableH;
    *fSplitPostStepPoint = *(step.GetPreStepPoint());
 
    // Split the energy
    // ----------------
    G4int numberVoxelsInStep = fpEnergySplitter->SplitEnergyInVolumes(&step);
 
    preStepPosition = step.GetPreStepPoint()->GetPosition();
    finalPostStepPosition = step.GetPostStepPoint()->GetPosition();
    direction = (finalPostStepPosition - preStepPosition).unit();
 
    fFinalTouchableH = track.GetNextTouchableHandle();
 
    postStepPosition = preStepPosition;
    // Loop over the sub-parts of this step
    G4int iStep;
    for (iStep = 0; iStep < numberVoxelsInStep; iStep++) {
      G4int idVoxel = -1;  // Voxel ID
      G4double stepLength = 0.0, energyLoss = 0.0;
 
      *fSplitPreStepPoint = *fSplitPostStepPoint;
      fOldTouchableH = fNewTouchableH;
 
      preStepPosition = postStepPosition;
      fSplitPreStepPoint->SetPosition(preStepPosition);
      fSplitPreStepPoint->SetTouchableHandle(fOldTouchableH);
 
      fpEnergySplitter->GetLengthAndEnergyDeposited(iStep, idVoxel, stepLength, energyLoss);
 
      // Correct the material, so that the track->GetMaterial gives correct answer
      pLogicalVolume->SetMaterial(fpEnergySplitter->GetVoxelMaterial(iStep));  // idVoxel) );
 
      postStepPosition = preStepPosition + stepLength * direction;
      fSplitPostStepPoint->SetPosition(postStepPosition);
 
      // Load the Step with the new values
      fSplitStep->SetStepLength(stepLength);
      fSplitStep->SetTotalEnergyDeposit(energyLoss);
      if (iStep < numberVoxelsInStep - 1) {
        fSplitStep->GetPostStepPoint()->SetStepStatus(fGeomBoundary);
        G4int nextVoxelId = -1;
        fpEnergySplitter->GetVoxelID(iStep + 1, nextVoxelId);
 
        // Create new "next" touchable for each section ??
        G4VTouchable* fNewTouchablePtr = CreateTouchableForSubStep(nextVoxelId, postStepPosition);
        fNewTouchableH = G4TouchableHandle(fNewTouchablePtr);
        fSplitPostStepPoint->SetTouchableHandle(fNewTouchableH);
      }
      else {
        fSplitStep->GetPostStepPoint()->SetStepStatus(fullStepStatus);
        fSplitPostStepPoint->SetTouchableHandle(fFinalTouchableH);
      }
 
      // As first approximation, split the NIEL in the same fractions as the energy deposit
      G4double eLossFraction;
      eLossFraction = (totalEnergyDeposit > 0.0) ? energyLoss / totalEnergyDeposit : 1.0;
      fSplitStep->SetNonIonizingEnergyDeposit(step.GetNonIonizingEnergyDeposit() * eLossFraction);
 
      fSplitPostStepPoint->SetSensitiveDetector(ptrSD);
 
      // Call the Sensitive Detector
      ptrSD->Hit(fSplitStep);
 
      if (verboseLevel > 1) Verbose(step);
    }
  }
 
  // This must change the Stepping Control
  return pParticleChange;
}
 
G4TouchableHistory* G4ScoreSplittingProcess::CreateTouchableForSubStep(G4int newVoxelNum,
                                                                       G4ThreeVector)
{
  auto oldTouchableHistory = dynamic_cast<G4TouchableHistory*>(fOldTouchableH());
  G4TouchableHistory* ptrTouchableHistory =
    G4TransportationManager::GetTransportationManager()
      ->GetNavigatorForTracking()
      ->CreateTouchableHistory(oldTouchableHistory->GetHistory());
 
  // Change the history
  auto ptrNavHistory = const_cast<G4NavigationHistory*>(ptrTouchableHistory->GetHistory());
  G4VPhysicalVolume* curPhysicalVol = ptrNavHistory->GetTopVolume();
 
  EVolume curVolumeType = ptrNavHistory->GetTopVolumeType();
  if (curVolumeType == kParameterised) {
    ptrNavHistory->BackLevel();
    // G4VPVParameterised parameterisedPV= pNewMother
    G4VPVParameterisation* curParamstn = curPhysicalVol->GetParameterisation();
 
    // From G4ParameterisedNavigation::IdentifyAndPlaceSolid() inline method
    G4VSolid* sampleSolid = curParamstn->ComputeSolid(newVoxelNum, curPhysicalVol);
    sampleSolid->ComputeDimensions(curParamstn, newVoxelNum, curPhysicalVol);
    curParamstn->ComputeTransformation(newVoxelNum, curPhysicalVol);
 
    ptrNavHistory->NewLevel(curPhysicalVol, kParameterised, newVoxelNum);
  }
  else {
    G4cout << " Current volume type is not Parameterised. " << G4endl;
    G4Exception(
      "G4ScoreSplittingProcess::CreateTouchableForSubStep", "ErrorRegularParamaterisation",
      JustWarning,
      "Score Splitting Process is used for Regular Structure - but did not find one here.");
  }
  return ptrTouchableHistory;
}
 
void G4ScoreSplittingProcess::CopyStepStart(const G4Step& step)
{
  fSplitStep->SetTrack(step.GetTrack());
  fSplitStep->SetStepLength(step.GetStepLength());
  fSplitStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
  fSplitStep->SetNonIonizingEnergyDeposit(step.GetNonIonizingEnergyDeposit());
  fSplitStep->SetControlFlag(step.GetControlFlag());
 
  *fSplitPreStepPoint = *(step.GetPreStepPoint());
 
  fInitialTouchableH = (step.GetPreStepPoint())->GetTouchableHandle();
  fFinalTouchableH = (step.GetPostStepPoint())->GetTouchableHandle();
}
 
void G4ScoreSplittingProcess::Verbose(const G4Step& step) const
{
  G4cout << "In mass geometry ------------------------------------------------" << G4endl;
  G4cout << " StepLength : " << step.GetStepLength() / mm
         << "      TotalEnergyDeposit : " << step.GetTotalEnergyDeposit() / MeV << G4endl;
  G4cout << " PreStepPoint : " << step.GetPreStepPoint()->GetPhysicalVolume()->GetName() << " - ";
  if (step.GetPreStepPoint()->GetProcessDefinedStep() != nullptr) {
    G4cout << step.GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
  }
  else {
    G4cout << "NoProcessAssigned";
  }
  G4cout << G4endl;
  G4cout << "                " << step.GetPreStepPoint()->GetPosition() << G4endl;
  G4cout << " PostStepPoint : ";
  if (step.GetPostStepPoint()->GetPhysicalVolume() != nullptr) {
    G4cout << step.GetPostStepPoint()->GetPhysicalVolume()->GetName();
  }
  else {
    G4cout << "OutOfWorld";
  }
  G4cout << " - ";
  if (step.GetPostStepPoint()->GetProcessDefinedStep() != nullptr) {
    G4cout << step.GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
  }
  else {
    G4cout << "NoProcessAssigned";
  }
  G4cout << G4endl;
  G4cout << "                 " << step.GetPostStepPoint()->GetPosition() << G4endl;
 
  G4cout << "In ghost geometry ------------------------------------------------" << G4endl;
  G4cout << " StepLength : " << fSplitStep->GetStepLength() / mm
         << "      TotalEnergyDeposit : " << fSplitStep->GetTotalEnergyDeposit() / MeV << G4endl;
  G4cout << " PreStepPoint : " << fSplitStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()
         << " [" << fSplitStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber() << " ]"
         << " - ";
  if (fSplitStep->GetPreStepPoint()->GetProcessDefinedStep() != nullptr) {
    G4cout << fSplitStep->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
  }
  else {
    G4cout << "NoProcessAssigned";
  }
  G4cout << G4endl;
  G4cout << "                " << fSplitStep->GetPreStepPoint()->GetPosition() << G4endl;
  G4cout << " PostStepPoint : ";
  if (fSplitStep->GetPostStepPoint()->GetPhysicalVolume() != nullptr) {
    G4cout << fSplitStep->GetPostStepPoint()->GetPhysicalVolume()->GetName() << " ["
           << fSplitStep->GetPostStepPoint()->GetTouchable()->GetReplicaNumber() << " ]";
  }
  else {
    G4cout << "OutOfWorld";
  }
  G4cout << " - ";
  if (fSplitStep->GetPostStepPoint()->GetProcessDefinedStep() != nullptr) {
    G4cout << fSplitStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
  }
  else {
    G4cout << "NoProcessAssigned";
  }
  G4cout << G4endl;
  G4cout << "                 " << fSplitStep->GetPostStepPoint()->GetPosition()
         << " == " << fSplitStep->GetTrack()->GetMomentumDirection() << G4endl;
}
 
//----------------------------------------------------------
//
//  AtRestGetPhysicalInteractionLength()
//
//----------------------------------------------------------
G4double G4ScoreSplittingProcess::AtRestGetPhysicalInteractionLength(const G4Track& /*track*/,
                                                                     G4ForceCondition* condition)
{
  *condition = NotForced;  // Was Forced
  return DBL_MAX;
}
 
//---------------------------------------
//  AlongStepGetPhysicalInteractionLength
//---------------------------------------
G4double
G4ScoreSplittingProcess::AlongStepGetPhysicalInteractionLength(const G4Track&,  // track,
                                                               G4double,  // previousStepSize,
                                                               G4double,  // currentMinimumStep,
                                                               G4double&,  // proposedSafety,
                                                               G4GPILSelection* selection)
{
  *selection = NotCandidateForSelection;
  return DBL_MAX;
}
 
//------------------------------------
//             AlongStepDoIt()
//------------------------------------
 
G4VParticleChange* G4ScoreSplittingProcess::AlongStepDoIt(const G4Track& track, const G4Step&)
{
  // Dummy ParticleChange ie: does nothing
  // Expecting G4Transportation to move the track
  dummyParticleChange.Initialize(track);
  return &dummyParticleChange;
}
 
//------------------------------------
//             AtRestDoIt()
//------------------------------------
G4VParticleChange* G4ScoreSplittingProcess::AtRestDoIt(const G4Track& track, const G4Step&)
{
  pParticleChange->Initialize(track);
  return pParticleChange;
}