Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EnergySplitter.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 #include "G4EnergySplitter.hh"
27 #include "G4VSolid.hh"
28 #include "G4UnitsTable.hh"
31 #include "G4EmCalculator.hh"
32 #include "G4PhysicalVolumeStore.hh"
33 #include "G4Step.hh"
34 #include "G4PVParameterised.hh"
35 
37 // (Description)
38 //
39 // Created:
40 //
42 
44 {
45  theElossExt = new G4EnergyLossForExtrapolator(0);
46  thePhantomParam = 0;
47  theNIterations = 2;
48 }
49 
51 {
52  delete theElossExt;
53 }
54 
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 }
258 
259 
260 //-----------------------------------------------------------------------
261 void G4EnergySplitter::GetPhantomParam(G4bool mustExist)
262 {
264  std::vector<G4VPhysicalVolume*>::iterator cite;
265  for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
266  // G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
267  if( IsPhantomVolume( *cite ) ) {
268  const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
269  G4VPVParameterisation* param = pvparam->GetParameterisation();
270  // if( static_cast<const G4PhantomParameterisation*>(param) ){
271  // if( static_cast<const G4PhantomParameterisation*>(param) ){
272  // G4cout << "G4PhantomParameterisation volume found " << (*cite)->GetName() << G4endl;
273  thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
274  }
275  }
276 
277  if( !thePhantomParam && mustExist )
278  G4Exception("G4EnergySplitter::GetPhantomParam",
279  "PhantomParamError", FatalException,
280  "No G4PhantomParameterisation found !");
281 }
282 
283 
284 //-----------------------------------------------------------------------
285 G4bool G4EnergySplitter::IsPhantomVolume( G4VPhysicalVolume* pv )
286 {
287  EAxis axis;
288  G4int nReplicas;
289  G4double width,offset;
290  G4bool consuming;
291  pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
292  EVolume type = (consuming) ? kReplica : kParameterised;
293  if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {
294  return TRUE;
295  } else {
296  return FALSE;
297  }
298 
299 }
300 
301 //-----------------------------------------------------------------------
303 {
304  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().begin())).first;
305 }
306 
307 //-----------------------------------------------------------------------
309 {
310  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().rbegin())).first;
311 }
312 
313 //-----------------------------------------------------------------------
314 void G4EnergySplitter::GetVoxelID( G4int stepNo, G4int& voxelID )
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 }
327 
328 
329 //-----------------------------------------------------------------------
330 void G4EnergySplitter::GetStepLength( G4int stepNo, G4double& stepLength )
331 {
332  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
333  advance( ite, stepNo );
334  stepLength = (*ite).second;
335 }
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 G4String ConvertToString(G4bool boolVal)
Definition: G4UIcommand.cc:372
#define width
static constexpr double second
Definition: G4SIunits.hh:157
virtual ~G4EnergySplitter()
G4double TrueStepLength(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *part)
int G4int
Definition: G4Types.hh:78
static G4PhysicalVolumeStore * GetInstance()
void GetFirstVoxelID(G4int &voxelID)
G4StepPoint * GetPreStepPoint() const
G4VPVParameterisation * GetParameterisation() const
G4GLOB_DLL std::ostream G4cout
void GetLastVoxelID(G4int &voxelID)
bool G4bool
Definition: G4Types.hh:79
virtual G4int GetRegularStructureId() const =0
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
#define FALSE
Definition: globals.hh:52
#define TRUE
Definition: globals.hh:55
Definition: G4Step.hh:76
G4double GetTotalEnergyDeposit() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4Material * GetMaterial(size_t nx, size_t ny, size_t nz) const
void GetVoxelID(G4int stepNo, G4int &voxelID)
EAxis
Definition: geomdefs.hh:54
virtual void GetReplicationData(EAxis &axis, G4int &nReplicas, G4double &width, G4double &offset, G4bool &consuming) const =0
#define G4endl
Definition: G4ios.hh:61
EVolume
Definition: geomdefs.hh:69
G4double GetKineticEnergy() const
double G4double
Definition: G4Types.hh:76
static G4RegularNavigationHelper * Instance()
G4Track * GetTrack() const
G4double GetPDGCharge() const
G4int SplitEnergyInVolumes(const G4Step *aStep)