G4EnergySplitter Class Reference

#include <G4EnergySplitter.hh>

Public Member Functions
	G4EnergySplitter ()
virtual	~G4EnergySplitter ()
G4int	SplitEnergyInVolumes (const G4Step *aStep)
void	GetLastVoxelID (G4int &voxelID)
void	GetFirstVoxelID (G4int &voxelID)
void	GetVoxelID (G4int stepNo, G4int &voxelID)
void	GetVoxelIDAndLength (G4int stepNo, G4int &voxelID, G4double &stepLength)
void	GetLengthAndEnergyDeposited (G4int stepNo, G4int &voxelID, G4double &stepLength, G4double &energyLoss)
void	GetLengthAndInitialEnergy (G4double &preStepEnergy, G4int stepNo, G4int &voxelID, G4double &stepLength, G4double &initialEnergy)
void	SetNIterations (G4int niter)
G4Material *	GetVoxelMaterial (G4int stepNo)

Detailed Description

Definition at line 49 of file G4EnergySplitter.hh.

Constructor & Destructor Documentation

G4EnergySplitter()

G4EnergySplitter::G4EnergySplitter

(

)

Definition at line 44 of file G4EnergySplitter.cc.

{
  theElossExt = new G4EnergyLossForExtrapolator(0);
  thePhantomParam = nullptr;
  theNIterations = 2;
}

~G4EnergySplitter()

G4EnergySplitter::~G4EnergySplitter ( )

virtual

Definition at line 51 of file G4EnergySplitter.cc.

{
  delete theElossExt;
}

Member Function Documentation

GetFirstVoxelID()

void G4EnergySplitter::GetFirstVoxelID

(

G4int &

voxelID

)

Definition at line 317 of file G4EnergySplitter.cc.

{
  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().crbegin())).first;
}

GetLastVoxelID()

void G4EnergySplitter::GetLastVoxelID

(

G4int &

voxelID

)

Definition at line 311 of file G4EnergySplitter.cc.

{
  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().cbegin())).first;
}

GetLengthAndEnergyDeposited()

void G4EnergySplitter::GetLengthAndEnergyDeposited ( G4int stepNo, G4int & voxelID, G4double & stepLength, G4double & energyLoss )

inline

GetLengthAndInitialEnergy()

void G4EnergySplitter::GetLengthAndInitialEnergy ( G4double & preStepEnergy, G4int stepNo, G4int & voxelID, G4double & stepLength, G4double & initialEnergy )

inline

GetVoxelID()

void G4EnergySplitter::GetVoxelID	(	G4int	stepNo,
		G4int &	voxelID )

Definition at line 323 of file G4EnergySplitter.cc.

{
  if (stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()))
  {
    G4Exception("G4EnergySplitter::GetVoxelID",
                "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
                FatalErrorInArgument,
                G4String("stepNo = " + G4UIcommand::ConvertToString(stepNo)
                         + ", number of voxels = "
                         + G4UIcommand::ConvertToString(
                           G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size())))
                  .c_str());
  }
  auto ite = G4RegularNavigationHelper::Instance()->GetStepLengths().cbegin();
  advance(ite, stepNo);
  voxelID = (*ite).first;
}

GetVoxelIDAndLength()

void G4EnergySplitter::GetVoxelIDAndLength ( G4int stepNo, G4int & voxelID, G4double & stepLength )

inline

GetVoxelMaterial()

G4Material * G4EnergySplitter::GetVoxelMaterial ( G4int stepNo )

inline

SetNIterations()

void G4EnergySplitter::SetNIterations ( G4int niter )

inline

SplitEnergyInVolumes()

G4int G4EnergySplitter::SplitEnergyInVolumes

(

const G4Step *

aStep

)

Definition at line 56 of file G4EnergySplitter.cc.

{
  theEnergies.clear();
 
  if (aStep == nullptr) return false;  // it is 0 when called by GmScoringMgr after last event
 
  G4double edep = aStep->GetTotalEnergyDeposit();
 
  if( edep == 0. ) {
    return 0; 
  }
  
#ifdef VERBOSE_ENERSPLIT
  G4bool verbose = 1;
  if (verbose)
    G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit()
           << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size()
           << G4endl;
#endif
  if (G4RegularNavigationHelper::Instance()->GetStepLengths().empty()
      || aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0)
  {  // we are only counting dose deposit
    return (G4int)theEnergies.size();
  }
  if (G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1) {
    theEnergies.push_back(edep);
    return (G4int)theEnergies.size();
  }
 
  if (thePhantomParam == nullptr) GetPhantomParam(true);
 
  //----- Distribute energy deposited in voxels
  std::vector<std::pair<G4int, G4double>> rnsl =
    G4RegularNavigationHelper::Instance()->GetStepLengths();
 
  const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
  G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
  G4double kinEnergyPre = kinEnergyPreOrig;
 
  G4double stepLength = aStep->GetStepLength();
  G4double slSum = 0.;
  unsigned int ii;
  for (ii = 0; ii < rnsl.size(); ++ii) {
    G4double sl = rnsl[ii].second;
    slSum += sl;
#ifdef VERBOSE_ENERSPLIT
    if (verbose)
      G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter1 step length geom "
             << sl << G4endl;
#endif
  }
 
#ifdef VERBOSE_ENERSPLIT
  if (verbose)
    G4cout << "G4EnergySplitter RN:  step length geom TOTAL " << slSum << " true TOTAL "
           << stepLength << " ratio " << stepLength / slSum << " Energy "
           << aStep->GetPreStepPoint()->GetKineticEnergy() << " Material "
           << aStep->GetPreStepPoint()->GetMaterial()->GetName() << " Number of geom steps "
           << rnsl.size() << G4endl;
#endif
  //----- No iterations to correct elost and msc => distribute energy deposited according to
  // geometrical step length in each voxel
  if (theNIterations == 0) {
    for (ii = 0; ii < rnsl.size(); ++ii) {
      G4double sl = rnsl[ii].second;
      G4double edepStep = edep * sl / slSum;  // divide edep along steps, proportional to step
                                              // length
#ifdef VERBOSE_ENERSPLIT
      if (verbose)
        G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " edep " << edepStep << G4endl;
#endif
 
      theEnergies.push_back(edepStep);
    }
  }
  else {  //  1 or more iterations demanded
 
#ifdef VERBOSE_ENERSPLIT
    // print corrected energy at iteration 0
    if (verbose) {
      G4double slSum = 0.;
      for (ii = 0; ii < rnsl.size(); ++ii) {
        G4double sl = rnsl[ii].second;
        slSum += sl;
      }
      for (ii = 0; ii < rnsl.size(); ii++) {
        G4cout << "G4EnergySplitter::SplitEnergyInVolumes " << ii
               << " RN: iter0 corrected energy lost " << edep * rnsl[ii].second / slSum << G4endl;
      }
    }
#endif
 
    G4double slRatio = stepLength / slSum;
#ifdef VERBOSE_ENERSPLIT
    if (verbose)
      G4cout << "G4EnergySplitter::SplitEnergyInVolumes  RN: iter 0, step ratio " << slRatio
             << G4endl;
#endif
 
    //--- energy at each interaction
    G4EmCalculator emcalc;
    G4double totalELost = 0.;
    std::vector<G4double> stepLengths;
    for (G4int iiter = 1; iiter <= theNIterations; ++iiter) {
      //--- iter1: distribute true step length in each voxel: geom SL in each voxel is multiplied by
      // a constant so that the sum gives the total true step length
      if (iiter == 1) {
        for (ii = 0; ii < rnsl.size(); ++ii) {
          G4double sl = rnsl[ii].second;
          stepLengths.push_back(sl * slRatio);
#ifdef VERBOSE_ENERSPLIT
          if (verbose)
            G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter
                   << " corrected step length " << sl * slRatio << G4endl;
#endif
        }
 
        for (ii = 0; ii < rnsl.size(); ++ii) {
          const G4Material* mate = thePhantomParam->GetMaterial(rnsl[ii].first);
          G4double dEdx = 0.;
          if (kinEnergyPre > 0.) {  // t check this
            dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
          }
          G4double elost = stepLengths[ii] * dEdx;
 
#ifdef VERBOSE_ENERSPLIT
          if (verbose)
            G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter1 energy lost "
                   << elost << " energy at interaction " << kinEnergyPre << " = stepLength "
                   << stepLengths[ii] << " * dEdx " << dEdx << G4endl;
#endif
          kinEnergyPre -= elost;
          theEnergies.push_back(elost);
          totalELost += elost;
        }
      }
      else {
        //------ 2nd and other iterations
        //----- Get step lengths corrected by changing geom2true correction
        //-- Get ratios for each energy
        slSum = 0.;
        kinEnergyPre = kinEnergyPreOrig;
        for (ii = 0; ii < rnsl.size(); ++ii) {
          const G4Material* mate = thePhantomParam->GetMaterial(rnsl[ii].first);
          stepLengths[ii] = theElossExt->TrueStepLength(kinEnergyPre, rnsl[ii].second, mate, part);
          kinEnergyPre -= theEnergies[ii];
 
#ifdef VERBOSE_ENERSPLIT
          if (verbose)
            G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter
                   << " step length geom " << stepLengths[ii] << " geom2true "
                   << rnsl[ii].second / stepLengths[ii] << G4endl;
#endif
 
          slSum += stepLengths[ii];
        }
 
        // Correct step lengths so that they sum the total step length
        G4double slratio = aStep->GetStepLength() / slSum;
#ifdef VERBOSE_ENERSPLIT
        if (verbose)
          G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter
                 << " step ratio " << slRatio << G4endl;
#endif
        for (ii = 0; ii < rnsl.size(); ++ii) {
          stepLengths[ii] *= slratio;
#ifdef VERBOSE_ENERSPLIT
          if (verbose)
            G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter
                   << " corrected step length " << stepLengths[ii] << G4endl;
#endif
        }
 
        //---- Recalculate energy lost with this new step lengths
        kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
        totalELost = 0.;
        for (ii = 0; ii < rnsl.size(); ++ii) {
          const G4Material* mate = thePhantomParam->GetMaterial(rnsl[ii].first);
          G4double dEdx = 0.;
          if (kinEnergyPre > 0.) {
            dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
          }
          G4double elost = stepLengths[ii] * dEdx;
#ifdef VERBOSE_ENERSPLIT
          if (verbose)
            G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter
                   << " energy lost " << elost << " energy at interaction " << kinEnergyPre
                   << " = stepLength " << stepLengths[ii] << " * dEdx " << dEdx << G4endl;
#endif
          kinEnergyPre -= elost;
          theEnergies[ii] = elost;
          totalELost += elost;
        }
      }
 
      // correct energies so that they reproduce the real step energy lost
      G4double enerRatio = (edep / totalELost);
 
#ifdef VERBOSE_ENERSPLIT
      if (verbose)
        G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter
               << " energy ratio " << enerRatio << G4endl;
#endif
 
#ifdef VERBOSE_ENERSPLIT
      G4double elostTot = 0.;
#endif
      for (ii = 0; ii < theEnergies.size(); ++ii) {
        theEnergies[ii] *= enerRatio;
#ifdef VERBOSE_ENERSPLIT
        elostTot += theEnergies[ii];
        if (verbose)
          G4cout << "G4EnergySplitter::SplitEnergyInVolumes " << ii << " RN: iter" << iiter
                 << " corrected energy lost " << theEnergies[ii] << " orig elost "
                 << theEnergies[ii] / enerRatio << " energy before interaction "
                 << kinEnergyPreOrig - elostTot + theEnergies[ii] << " energy after interaction "
                 << kinEnergyPreOrig - elostTot << G4endl;
#endif
      }
    }
  }
 
  return (G4int)theEnergies.size();
}

---

The documentation for this class was generated from the following files:

• G4EnergySplitter.hh
• G4EnergySplitter.cc