Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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)
 
G4MaterialGetVoxelMaterial (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.

44 {
45  theElossExt = new G4EnergyLossForExtrapolator(0);
46  thePhantomParam = 0;
47  theNIterations = 2;
48 }
G4EnergySplitter::~G4EnergySplitter ( )
virtual

Definition at line 50 of file G4EnergySplitter.cc.

51 {
52  delete theElossExt;
53 }

Member Function Documentation

void G4EnergySplitter::GetFirstVoxelID ( G4int voxelID)

Definition at line 308 of file G4EnergySplitter.cc.

309 {
310  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().rbegin())).first;
311 }
const std::vector< std::pair< G4int, G4double > > & GetStepLengths()
static G4RegularNavigationHelper * Instance()

Here is the call graph for this function:

void G4EnergySplitter::GetLastVoxelID ( G4int voxelID)

Definition at line 302 of file G4EnergySplitter.cc.

303 {
304  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().begin())).first;
305 }
const std::vector< std::pair< G4int, G4double > > & GetStepLengths()
static G4RegularNavigationHelper * Instance()

Here is the call graph for this function:

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

Here is the caller graph for this function:

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

Definition at line 314 of file G4EnergySplitter.cc.

315 {
316  if( stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()) ) {
317  G4Exception("G4EnergySplitter::GetVoxelID",
318  "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
320  G4String("stepNo = " + G4UIcommand::ConvertToString(stepNo) + ", number of voxels = " + G4UIcommand::ConvertToString(G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size())) ).c_str());
321  }
322  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
323  advance( ite, stepNo );
324  voxelID = (*ite).first;
325 
326 }
const std::vector< std::pair< G4int, G4double > > & GetStepLengths()
static G4String ConvertToString(G4bool boolVal)
Definition: G4UIcommand.cc:372
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4RegularNavigationHelper * Instance()

Here is the call graph for this function:

Here is the caller graph for this function:

void G4EnergySplitter::GetVoxelIDAndLength ( G4int  stepNo,
G4int voxelID,
G4double stepLength 
)
inline
G4Material* G4EnergySplitter::GetVoxelMaterial ( G4int  stepNo)
inline

Here is the caller graph for this function:

void G4EnergySplitter::SetNIterations ( G4int  niter)
inline
G4int G4EnergySplitter::SplitEnergyInVolumes ( const G4Step aStep)

Definition at line 55 of file G4EnergySplitter.cc.

56 {
57  theEnergies.clear();
58 
59  G4double edep = aStep->GetTotalEnergyDeposit();
60 
61 #ifdef VERBOSE_ENERSPLIT
62  G4bool verbose = 1;
63  if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit()
64  << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size() << G4endl;
65 #endif
66  if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 0 ||
67  aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0) { // we are only counting dose deposit
68  return theEnergies.size();
69  }
70  if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
71  theEnergies.push_back(edep);
72  return theEnergies.size();
73  }
74 
75  if( !thePhantomParam ) GetPhantomParam(TRUE);
76 
77  if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
78 
79  //----- Distribute energy deposited in voxels
80  std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::Instance()->GetStepLengths();
81 
82  const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
83  G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
84  G4double kinEnergyPre = kinEnergyPreOrig;
85 
86  G4double stepLength = aStep->GetStepLength();
87  G4double slSum = 0.;
88  unsigned int ii;
89  for( ii = 0; ii < rnsl.size(); ii++ ){
90  G4double sl = rnsl[ii].second;
91  slSum += sl;
92 #ifdef VERBOSE_ENERSPLIT
93  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
94 #endif
95  }
96 
97 #ifdef VERBOSE_ENERSPLIT
98  if( verbose )
99  G4cout << "G4EnergySplitter RN: step length geom TOTAL " << slSum
100  << " true TOTAL " << stepLength
101  << " ratio " << stepLength/slSum
102  << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy()
103  << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName()
104  << " Number of geom steps " << rnsl.size() << G4endl;
105 #endif
106  //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
107  if( theNIterations == 0 ) {
108  for( ii = 0; ii < rnsl.size(); ii++ ){
109  G4double sl = rnsl[ii].second;
110  G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
111 #ifdef VERBOSE_ENERSPLIT
112  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii
113  << " edep " << edepStep << G4endl;
114 #endif
115 
116  theEnergies.push_back(edepStep);
117 
118  }
119  } else { // 1 or more iterations demanded
120 
121 #ifdef VERBOSE_ENERSPLIT
122  // print corrected energy at iteration 0
123  if(verbose) {
124  G4double slSum = 0.;
125  for( ii = 0; ii < rnsl.size(); ii++ ){
126  G4double sl = rnsl[ii].second;
127  slSum += sl;
128  }
129  for( ii = 0; ii < rnsl.size(); ii++ ){
130  G4cout << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
131  << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum
132  << G4endl;
133  }
134  }
135 #endif
136 
137  G4double slRatio = stepLength/slSum;
138 #ifdef VERBOSE_ENERSPLIT
139  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes RN: iter 0, step ratio " << slRatio << G4endl;
140 #endif
141 
142  //--- energy at each interaction
143  G4EmCalculator emcalc;
144  G4double totalELost = 0.;
145  std::vector<G4double> stepLengths;
146  for( int iiter = 1; iiter <= theNIterations; iiter++ ) {
147  //--- 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
148  if( iiter == 1 ) {
149  for( ii = 0; ii < rnsl.size(); ii++ ){
150  G4double sl = rnsl[ii].second;
151  stepLengths.push_back( sl * slRatio );
152 #ifdef VERBOSE_ENERSPLIT
153  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
154 #endif
155  }
156 
157  for( ii = 0; ii < rnsl.size(); ii++ ){
158  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
159  G4double dEdx = 0.;
160  if( kinEnergyPre > 0. ) { //t check this
161  dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
162  }
163  G4double elost = stepLengths[ii] * dEdx;
164 
165 #ifdef VERBOSE_ENERSPLIT
166  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost " << elost
167  << " energy at interaction " << kinEnergyPre
168  << " = stepLength " << stepLengths[ii]
169  << " * dEdx " << dEdx << G4endl;
170 #endif
171  kinEnergyPre -= elost;
172  theEnergies.push_back( elost );
173  totalELost += elost;
174  }
175 
176  } else{
177  //------ 2nd and other iterations
178  //----- Get step lengths corrected by changing geom2true correction
179  //-- Get ratios for each energy
180  slSum = 0.;
181  kinEnergyPre = kinEnergyPreOrig;
182  for( ii = 0; ii < rnsl.size(); ii++ ){
183  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
184  stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
185  kinEnergyPre -= theEnergies[ii];
186 
187 #ifdef VERBOSE_ENERSPLIT
188  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii
189  << " RN: iter" << iiter << " step length geom " << stepLengths[ii]
190  << " geom2true " << rnsl[ii].second / stepLengths[ii] << G4endl;
191 #endif
192 
193  slSum += stepLengths[ii];
194  }
195 
196  //Correct step lengths so that they sum the total step length
197  G4double slratio = aStep->GetStepLength()/slSum;
198 #ifdef VERBOSE_ENERSPLIT
199  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
200 #endif
201  for( ii = 0; ii < rnsl.size(); ii++ ){
202  stepLengths[ii] *= slratio;
203 #ifdef VERBOSE_ENERSPLIT
204  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
205 #endif
206  }
207 
208  //---- Recalculate energy lost with this new step lengths
209  kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
210  totalELost = 0.;
211  for( ii = 0; ii < rnsl.size(); ii++ ){
212  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
213  G4double dEdx = 0.;
214  if( kinEnergyPre > 0. ) {
215  dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
216  }
217  G4double elost = stepLengths[ii] * dEdx;
218 #ifdef VERBOSE_ENERSPLIT
219  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost
220  << " energy at interaction " << kinEnergyPre
221  << " = stepLength " << stepLengths[ii]
222  << " * dEdx " << dEdx << G4endl;
223 #endif
224  kinEnergyPre -= elost;
225  theEnergies[ii] = elost;
226  totalELost += elost;
227  }
228 
229  }
230 
231  //correct energies so that they reproduce the real step energy lost
232  G4double enerRatio = (edep/totalELost);
233 
234 #ifdef VERBOSE_ENERSPLIT
235  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
236 #endif
237 
238 #ifdef VERBOSE_ENERSPLIT
239  G4double elostTot = 0.;
240 #endif
241  for( ii = 0; ii < theEnergies.size(); ii++ ){
242  theEnergies[ii] *= enerRatio;
243 #ifdef VERBOSE_ENERSPLIT
244  elostTot += theEnergies[ii];
245  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii]
246  << " orig elost " << theEnergies[ii]/enerRatio
247  << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
248  << " energy after interaction " << kinEnergyPreOrig-elostTot
249  << G4endl;
250 #endif
251  }
252  }
253 
254  }
255 
256  return theEnergies.size();
257 }
G4ParticleDefinition * GetDefinition() const
G4double GetStepLength() const
const G4String & GetName() const
Definition: G4Material.hh:178
G4Material * GetMaterial() const
const std::vector< std::pair< G4int, G4double > > & GetStepLengths()
static constexpr double second
Definition: G4SIunits.hh:157
G4double TrueStepLength(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *part)
G4StepPoint * GetPreStepPoint() const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
#define FALSE
Definition: globals.hh:52
#define TRUE
Definition: globals.hh:55
G4double GetTotalEnergyDeposit() const
G4Material * GetMaterial(size_t nx, size_t ny, size_t nz) const
#define G4endl
Definition: G4ios.hh:61
G4double GetKineticEnergy() const
double G4double
Definition: G4Types.hh:76
static G4RegularNavigationHelper * Instance()
G4Track * GetTrack() const
G4double GetPDGCharge() const

Here is the call graph for this function:

Here is the caller graph for this function:


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