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 46 of file G4EnergySplitter.hh.

Constructor & Destructor Documentation

G4EnergySplitter::G4EnergySplitter (  ) 

Definition at line 43 of file G4EnergySplitter.cc.

00044 {
00045    theElossExt = new G4EnergyLossForExtrapolator(0);
00046    thePhantomParam = 0;
00047    theNIterations = 2;
00048 }

G4EnergySplitter::~G4EnergySplitter (  )  [virtual]

Definition at line 50 of file G4EnergySplitter.cc.

00051 {;}

Member Function Documentation

void G4EnergySplitter::GetFirstVoxelID ( G4int &  voxelID  ) 

Definition at line 306 of file G4EnergySplitter.cc.

References G4RegularNavigationHelper::GetStepLengths(), and G4RegularNavigationHelper::Instance().

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

void G4EnergySplitter::GetLastVoxelID ( G4int &  voxelID  ) 

Definition at line 300 of file G4EnergySplitter.cc.

References G4RegularNavigationHelper::GetStepLengths(), and G4RegularNavigationHelper::Instance().

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

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

Definition at line 37 of file G4EnergySplitter.icc.

References GetVoxelIDAndLength().

Referenced by G4ScoreSplittingProcess::PostStepDoIt().

00038 {
00039   GetVoxelIDAndLength( stepNo, voxelID, stepLength );
00040 
00041   energyLoss = theEnergies[stepNo];
00042 }

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

Definition at line 45 of file G4EnergySplitter.icc.

References GetVoxelIDAndLength().

00046 {
00047   GetVoxelIDAndLength( stepNo, voxelID, stepLength );
00048 
00049   initialEnergy = preStepEnergy;
00050   for( G4int ii = 0; ii < stepNo; ii++ ){
00051     initialEnergy -= theEnergies[stepNo];
00052   }
00053 }

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

Definition at line 312 of file G4EnergySplitter.cc.

References G4UIcommand::ConvertToString(), FatalErrorInArgument, G4Exception(), G4RegularNavigationHelper::GetStepLengths(), and G4RegularNavigationHelper::Instance().

Referenced by GetVoxelIDAndLength(), GetVoxelMaterial(), and G4ScoreSplittingProcess::PostStepDoIt().

00313 {
00314   if( stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()) ) {
00315   G4Exception("G4EnergySplitter::GetVoxelID",
00316               "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
00317               FatalErrorInArgument,
00318               G4String("stepNo = " + G4UIcommand::ConvertToString(stepNo) + ", number of voxels = " + G4UIcommand::ConvertToString(G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size())) ).c_str());
00319   }
00320   std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
00321   advance( ite, stepNo );
00322   voxelID = (*ite).first;
00323 
00324 }

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

Definition at line 30 of file G4EnergySplitter.icc.

References GetVoxelID().

Referenced by GetLengthAndEnergyDeposited(), and GetLengthAndInitialEnergy().

00031 {
00032   GetVoxelID( stepNo, voxelID );
00033   GetStepLength( stepNo, stepLength);
00034 }

G4Material * G4EnergySplitter::GetVoxelMaterial ( G4int  stepNo  )  [inline]

Definition at line 62 of file G4EnergySplitter.icc.

References G4PhantomParameterisation::GetMaterial(), GetVoxelID(), and TRUE.

Referenced by G4ScoreSplittingProcess::PostStepDoIt().

00063 {
00064   if( !thePhantomParam ) GetPhantomParam(TRUE);
00065   G4int voxelID;
00066   GetVoxelID( stepNo, voxelID );
00067   return thePhantomParam->GetMaterial( voxelID );
00068 }

void G4EnergySplitter::SetNIterations ( G4int  niter  )  [inline]

Definition at line 56 of file G4EnergySplitter.icc.

00057 {
00058   theNIterations = niter;
00059 }

G4int G4EnergySplitter::SplitEnergyInVolumes ( const G4Step *  aStep  ) 

Definition at line 53 of file G4EnergySplitter.cc.

References FALSE, G4cout, G4endl, G4EmCalculator::GetDEDX(), G4Track::GetDefinition(), G4StepPoint::GetKineticEnergy(), G4PhantomParameterisation::GetMaterial(), G4StepPoint::GetMaterial(), G4Material::GetName(), G4ParticleDefinition::GetPDGCharge(), G4Step::GetPreStepPoint(), G4Step::GetStepLength(), G4RegularNavigationHelper::GetStepLengths(), G4Step::GetTotalEnergyDeposit(), G4Step::GetTrack(), G4RegularNavigationHelper::Instance(), TRUE, and G4EnergyLossForExtrapolator::TrueStepLength().

Referenced by G4ScoreSplittingProcess::PostStepDoIt().

00054 {
00055   theEnergies.clear();
00056 
00057   G4double edep = aStep->GetTotalEnergyDeposit();
00058  
00059 #ifdef VERBOSE_ENERSPLIT
00060   G4bool verbose = 1;
00061   if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit() 
00062                        << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size() << G4endl;
00063 #endif    
00064   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 0 ||
00065       aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0)  { // we are only counting dose deposit
00066     return theEnergies.size();
00067   }
00068   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
00069     theEnergies.push_back(edep);
00070     return theEnergies.size();
00071   }
00072 
00073   if( !thePhantomParam ) GetPhantomParam(TRUE);
00074   
00075   if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
00076   
00077   //----- Distribute energy deposited in voxels 
00078   std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::Instance()->GetStepLengths(); 
00079 
00080   const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
00081   G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
00082   G4double kinEnergyPre = kinEnergyPreOrig;
00083   
00084   G4double stepLength = aStep->GetStepLength();
00085   G4double slSum = 0.;
00086   unsigned int ii;
00087   for( ii = 0; ii < rnsl.size(); ii++ ){
00088     G4double sl = rnsl[ii].second;
00089     slSum += sl;
00090 #ifdef VERBOSE_ENERSPLIT
00091     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
00092 #endif
00093   }
00094   
00095 #ifdef VERBOSE_ENERSPLIT
00096   if( verbose )
00097     G4cout << "G4EnergySplitter RN:  step length geom TOTAL " << slSum 
00098            << " true TOTAL " << stepLength 
00099            << " ratio " << stepLength/slSum 
00100            << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy() 
00101            << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName() 
00102            << " Number of geom steps " << rnsl.size() << G4endl;
00103 #endif
00104   //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
00105   if( theNIterations == 0 ) { 
00106     for( ii = 0; ii < rnsl.size(); ii++ ){
00107       G4double sl = rnsl[ii].second;
00108       G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
00109 #ifdef VERBOSE_ENERSPLIT
00110       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii 
00111                           << " edep " << edepStep << G4endl;
00112 #endif
00113         
00114       theEnergies.push_back(edepStep);
00115         
00116     }
00117   } else { //  1 or more iterations demanded
00118       
00119 #ifdef VERBOSE_ENERSPLIT
00120     // print corrected energy at iteration 0 
00121     if(verbose)  {
00122       G4double slSum = 0.;
00123       for( ii = 0; ii < rnsl.size(); ii++ ){
00124         G4double sl = rnsl[ii].second;
00125         slSum += sl;
00126       }
00127       for( ii = 0; ii < rnsl.size(); ii++ ){
00128         G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
00129                 << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum  
00130                 << G4endl;
00131       }
00132     }
00133 #endif
00134 
00135     G4double slRatio = stepLength/slSum;
00136 #ifdef VERBOSE_ENERSPLIT
00137     if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes  RN: iter 0, step ratio " << slRatio << G4endl;
00138 #endif
00139       
00140     //--- energy at each interaction
00141     G4EmCalculator emcalc;
00142     G4double totalELost = 0.;
00143     std::vector<G4double> stepLengths;
00144     for( int iiter = 1; iiter <= theNIterations; iiter++ ) {
00145       //--- 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
00146       if( iiter == 1 ) {
00147         for( ii = 0; ii < rnsl.size(); ii++ ){
00148           G4double sl = rnsl[ii].second;
00149           stepLengths.push_back( sl * slRatio );
00150 #ifdef VERBOSE_ENERSPLIT
00151           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
00152 #endif
00153         }
00154         
00155         for( ii = 0; ii < rnsl.size(); ii++ ){
00156           const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
00157           G4double dEdx = 0.;
00158           if( kinEnergyPre > 0. ) {  //t check this 
00159             dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
00160           }
00161           G4double elost = stepLengths[ii] * dEdx;
00162           
00163 #ifdef VERBOSE_ENERSPLIT
00164           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost "  << elost 
00165                               << " energy at interaction " << kinEnergyPre 
00166                               << " = stepLength " << stepLengths[ii] 
00167                               << " * dEdx " << dEdx << G4endl;
00168 #endif
00169           kinEnergyPre -= elost;
00170           theEnergies.push_back( elost );
00171           totalELost += elost;
00172         }
00173         
00174       } else{
00175         //------ 2nd and other iterations
00176         //----- Get step lengths corrected by changing geom2true correction
00177         //-- Get ratios for each energy 
00178         slSum = 0.;
00179         kinEnergyPre = kinEnergyPreOrig;
00180         for( ii = 0; ii < rnsl.size(); ii++ ){
00181           const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
00182           stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
00183           kinEnergyPre -= theEnergies[ii];
00184           
00185 #ifdef VERBOSE_ENERSPLIT
00186           if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii 
00187                              << " RN: iter" << iiter << " step length geom " << stepLengths[ii] 
00188                              << " geom2true " << rnsl[ii].second / stepLengths[ii]  << G4endl;
00189 #endif
00190           
00191           slSum += stepLengths[ii];
00192         }
00193         
00194         //Correct step lengths so that they sum the total step length
00195         G4double slratio = aStep->GetStepLength()/slSum;
00196 #ifdef VERBOSE_ENERSPLIT
00197         if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
00198 #endif
00199         for( ii = 0; ii < rnsl.size(); ii++ ){
00200           stepLengths[ii] *= slratio;
00201 #ifdef VERBOSE_ENERSPLIT
00202           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
00203 #endif
00204         }
00205         
00206         //---- Recalculate energy lost with this new step lengths
00207         kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
00208         totalELost = 0.;
00209         for( ii = 0; ii < rnsl.size(); ii++ ){
00210           const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
00211           G4double dEdx = 0.;
00212           if( kinEnergyPre > 0. ) {
00213             dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
00214           }
00215           G4double elost = stepLengths[ii] * dEdx;
00216 #ifdef VERBOSE_ENERSPLIT
00217           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost 
00218                               << " energy at interaction " << kinEnergyPre 
00219                               << " = stepLength " << stepLengths[ii] 
00220                               << " * dEdx " << dEdx << G4endl;
00221 #endif
00222           kinEnergyPre -= elost;
00223           theEnergies[ii] = elost;
00224           totalELost += elost;
00225         }
00226         
00227       }
00228       
00229       //correct energies so that they reproduce the real step energy lost
00230       G4double enerRatio = (edep/totalELost);
00231       
00232 #ifdef VERBOSE_ENERSPLIT
00233       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
00234 #endif
00235         
00236 #ifdef VERBOSE_ENERSPLIT
00237       G4double elostTot = 0.; 
00238 #endif
00239       for( ii = 0; ii < theEnergies.size(); ii++ ){
00240         theEnergies[ii] *= enerRatio;
00241 #ifdef VERBOSE_ENERSPLIT
00242         elostTot += theEnergies[ii];
00243         if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii] 
00244                             << " orig elost " << theEnergies[ii]/enerRatio 
00245                             << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
00246                             << " energy after interaction " << kinEnergyPreOrig-elostTot
00247                             << G4endl;
00248 #endif
00249       }
00250     }
00251     
00252   }
00253   
00254   return theEnergies.size();
00255 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:54 2013 for Geant4 by  doxygen 1.4.7