#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) |
Definition at line 46 of file G4EnergySplitter.hh.
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] |
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 }
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 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 }