G4EnergySplitter Class Reference

#include <G4EnergySplitter.hh>

Collaboration diagram for G4EnergySplitter:

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)

Private Member Functions
void	GetStepLength (G4int stepNo, G4double &stepLength)
void	GetPhantomParam (G4bool mustExist)
G4bool	IsPhantomVolume (G4VPhysicalVolume *pv)

Private Attributes
G4EnergyLossForExtrapolator *	theElossExt
G4int	theNIterations
std::vector< G4double >	theEnergies
G4PhantomParameterisation *	thePhantomParam

Detailed Description

Definition at line 46 of file G4EnergySplitter.hh.

Constructor & Destructor Documentation

G4EnergySplitter()

G4EnergySplitter::G4EnergySplitter

(

)

Definition at line 43 of file G4EnergySplitter.cc.

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

~G4EnergySplitter()

G4EnergySplitter::~G4EnergySplitter ( )

virtual

Definition at line 50 of file G4EnergySplitter.cc.

{
    delete theElossExt;
}

Member Function Documentation

GetFirstVoxelID()

void G4EnergySplitter::GetFirstVoxelID

(

G4int &

voxelID

)

Definition at line 308 of file G4EnergySplitter.cc.

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

Here is the call graph for this function:

GetLastVoxelID()

void G4EnergySplitter::GetLastVoxelID

(

G4int &

voxelID

)

Definition at line 302 of file G4EnergySplitter.cc.

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

Here is the call graph for this function:

GetLengthAndEnergyDeposited()

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

inline

Here is the caller graph for this function:

GetLengthAndInitialEnergy()

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

inline

GetPhantomParam()

void G4EnergySplitter::GetPhantomParam ( G4bool  mustExist )

private

Definition at line 261 of file G4EnergySplitter.cc.

{
  G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
  std::vector<G4VPhysicalVolume*>::iterator cite;
  for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
    //    G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
    if( IsPhantomVolume( *cite ) ) {
      const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
      G4VPVParameterisation* param = pvparam->GetParameterisation();
      //    if( static_cast<const G4PhantomParameterisation*>(param) ){
      //    if( static_cast<const G4PhantomParameterisation*>(param) ){
      //      G4cout << "G4PhantomParameterisation volume found  " << (*cite)->GetName() << G4endl;
      thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
    }
  }
  
  if( !thePhantomParam && mustExist )
    G4Exception("G4EnergySplitter::GetPhantomParam",
                "PhantomParamError", FatalException,
                "No G4PhantomParameterisation found !");
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetStepLength()

void G4EnergySplitter::GetStepLength ( G4int  stepNo, G4double &  stepLength  )

private

Definition at line 330 of file G4EnergySplitter.cc.

{
  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
  advance( ite, stepNo );
  stepLength = (*ite).second;
}

Here is the call graph for this function:

GetVoxelID()

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

Definition at line 314 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());
  }
  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
  advance( ite, stepNo );
  voxelID = (*ite).first;

}

Here is the call graph for this function:

Here is the caller graph for this function:

GetVoxelIDAndLength()

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

inline

GetVoxelMaterial()

G4Material* G4EnergySplitter::GetVoxelMaterial ( G4int  stepNo )

inline

Here is the caller graph for this function:

IsPhantomVolume()

G4bool G4EnergySplitter::IsPhantomVolume ( G4VPhysicalVolume *  pv )

private

Definition at line 285 of file G4EnergySplitter.cc.

{
  EAxis axis;
  G4int nReplicas;
  G4double width,offset;
  G4bool consuming;
  pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
  EVolume type = (consuming) ? kReplica : kParameterised;
  if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {  
    return TRUE;
  } else {
    return FALSE;
  }

} 

Here is the call graph for this function:

Here is the caller graph for this function:

SetNIterations()

void G4EnergySplitter::SetNIterations ( G4int  niter )

inline

SplitEnergyInVolumes()

G4int G4EnergySplitter::SplitEnergyInVolumes

(

const G4Step *

aStep

)

Definition at line 55 of file G4EnergySplitter.cc.

{
  theEnergies.clear();

  G4double edep = aStep->GetTotalEnergyDeposit();
 
#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().size() == 0 ||
      aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0)  { // we are only counting dose deposit
    return theEnergies.size();
  }
  if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
    theEnergies.push_back(edep);
    return theEnergies.size();
  }

  if( !thePhantomParam ) GetPhantomParam(TRUE);
  
  if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
  
  //----- 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( int 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 theEnergies.size();
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

theElossExt

G4EnergyLossForExtrapolator* G4EnergySplitter::theElossExt

private

Definition at line 71 of file G4EnergySplitter.hh.

theEnergies

std::vector<G4double> G4EnergySplitter::theEnergies

private

Definition at line 75 of file G4EnergySplitter.hh.

theNIterations

G4int G4EnergySplitter::theNIterations

private

Definition at line 73 of file G4EnergySplitter.hh.

thePhantomParam

G4PhantomParameterisation* G4EnergySplitter::thePhantomParam

private

Definition at line 77 of file G4EnergySplitter.hh.

---

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

• G4EnergySplitter.hh
• G4EnergySplitter.cc