Geant4_10
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 
54 {
55  theEnergies.clear();
56 
58 
59 #ifdef VERBOSE_ENERSPLIT
60  G4bool verbose = 1;
61  if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit()
62  << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size() << G4endl;
63 #endif
64  if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 0 ||
65  aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0) { // we are only counting dose deposit
66  return theEnergies.size();
67  }
68  if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
69  theEnergies.push_back(edep);
70  return theEnergies.size();
71  }
72 
73  if( !thePhantomParam ) GetPhantomParam(TRUE);
74 
75  if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
76 
77  //----- Distribute energy deposited in voxels
78  std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::Instance()->GetStepLengths();
79 
80  const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
81  G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
82  G4double kinEnergyPre = kinEnergyPreOrig;
83 
84  G4double stepLength = aStep->GetStepLength();
85  G4double slSum = 0.;
86  unsigned int ii;
87  for( ii = 0; ii < rnsl.size(); ii++ ){
88  G4double sl = rnsl[ii].second;
89  slSum += sl;
90 #ifdef VERBOSE_ENERSPLIT
91  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
92 #endif
93  }
94 
95 #ifdef VERBOSE_ENERSPLIT
96  if( verbose )
97  G4cout << "G4EnergySplitter RN: step length geom TOTAL " << slSum
98  << " true TOTAL " << stepLength
99  << " ratio " << stepLength/slSum
100  << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy()
101  << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName()
102  << " Number of geom steps " << rnsl.size() << G4endl;
103 #endif
104  //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
105  if( theNIterations == 0 ) {
106  for( ii = 0; ii < rnsl.size(); ii++ ){
107  G4double sl = rnsl[ii].second;
108  G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
109 #ifdef VERBOSE_ENERSPLIT
110  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii
111  << " edep " << edepStep << G4endl;
112 #endif
113 
114  theEnergies.push_back(edepStep);
115 
116  }
117  } else { // 1 or more iterations demanded
118 
119 #ifdef VERBOSE_ENERSPLIT
120  // print corrected energy at iteration 0
121  if(verbose) {
122  G4double slSum = 0.;
123  for( ii = 0; ii < rnsl.size(); ii++ ){
124  G4double sl = rnsl[ii].second;
125  slSum += sl;
126  }
127  for( ii = 0; ii < rnsl.size(); ii++ ){
128  G4cout << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
129  << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum
130  << G4endl;
131  }
132  }
133 #endif
134 
135  G4double slRatio = stepLength/slSum;
136 #ifdef VERBOSE_ENERSPLIT
137  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes RN: iter 0, step ratio " << slRatio << G4endl;
138 #endif
139 
140  //--- energy at each interaction
141  G4EmCalculator emcalc;
142  G4double totalELost = 0.;
143  std::vector<G4double> stepLengths;
144  for( int iiter = 1; iiter <= theNIterations; iiter++ ) {
145  //--- 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
146  if( iiter == 1 ) {
147  for( ii = 0; ii < rnsl.size(); ii++ ){
148  G4double sl = rnsl[ii].second;
149  stepLengths.push_back( sl * slRatio );
150 #ifdef VERBOSE_ENERSPLIT
151  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
152 #endif
153  }
154 
155  for( ii = 0; ii < rnsl.size(); ii++ ){
156  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
157  G4double dEdx = 0.;
158  if( kinEnergyPre > 0. ) { //t check this
159  dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
160  }
161  G4double elost = stepLengths[ii] * dEdx;
162 
163 #ifdef VERBOSE_ENERSPLIT
164  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost " << elost
165  << " energy at interaction " << kinEnergyPre
166  << " = stepLength " << stepLengths[ii]
167  << " * dEdx " << dEdx << G4endl;
168 #endif
169  kinEnergyPre -= elost;
170  theEnergies.push_back( elost );
171  totalELost += elost;
172  }
173 
174  } else{
175  //------ 2nd and other iterations
176  //----- Get step lengths corrected by changing geom2true correction
177  //-- Get ratios for each energy
178  slSum = 0.;
179  kinEnergyPre = kinEnergyPreOrig;
180  for( ii = 0; ii < rnsl.size(); ii++ ){
181  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
182  stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
183  kinEnergyPre -= theEnergies[ii];
184 
185 #ifdef VERBOSE_ENERSPLIT
186  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii
187  << " RN: iter" << iiter << " step length geom " << stepLengths[ii]
188  << " geom2true " << rnsl[ii].second / stepLengths[ii] << G4endl;
189 #endif
190 
191  slSum += stepLengths[ii];
192  }
193 
194  //Correct step lengths so that they sum the total step length
195  G4double slratio = aStep->GetStepLength()/slSum;
196 #ifdef VERBOSE_ENERSPLIT
197  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
198 #endif
199  for( ii = 0; ii < rnsl.size(); ii++ ){
200  stepLengths[ii] *= slratio;
201 #ifdef VERBOSE_ENERSPLIT
202  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
203 #endif
204  }
205 
206  //---- Recalculate energy lost with this new step lengths
207  kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
208  totalELost = 0.;
209  for( ii = 0; ii < rnsl.size(); ii++ ){
210  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
211  G4double dEdx = 0.;
212  if( kinEnergyPre > 0. ) {
213  dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
214  }
215  G4double elost = stepLengths[ii] * dEdx;
216 #ifdef VERBOSE_ENERSPLIT
217  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost
218  << " energy at interaction " << kinEnergyPre
219  << " = stepLength " << stepLengths[ii]
220  << " * dEdx " << dEdx << G4endl;
221 #endif
222  kinEnergyPre -= elost;
223  theEnergies[ii] = elost;
224  totalELost += elost;
225  }
226 
227  }
228 
229  //correct energies so that they reproduce the real step energy lost
230  G4double enerRatio = (edep/totalELost);
231 
232 #ifdef VERBOSE_ENERSPLIT
233  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
234 #endif
235 
236 #ifdef VERBOSE_ENERSPLIT
237  G4double elostTot = 0.;
238 #endif
239  for( ii = 0; ii < theEnergies.size(); ii++ ){
240  theEnergies[ii] *= enerRatio;
241 #ifdef VERBOSE_ENERSPLIT
242  elostTot += theEnergies[ii];
243  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii]
244  << " orig elost " << theEnergies[ii]/enerRatio
245  << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
246  << " energy after interaction " << kinEnergyPreOrig-elostTot
247  << G4endl;
248 #endif
249  }
250  }
251 
252  }
253 
254  return theEnergies.size();
255 }
256 
257 
258 //-----------------------------------------------------------------------
259 void G4EnergySplitter::GetPhantomParam(G4bool mustExist)
260 {
262  std::vector<G4VPhysicalVolume*>::iterator cite;
263  for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
264  // G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
265  if( IsPhantomVolume( *cite ) ) {
266  const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
267  G4VPVParameterisation* param = pvparam->GetParameterisation();
268  // if( static_cast<const G4PhantomParameterisation*>(param) ){
269  // if( static_cast<const G4PhantomParameterisation*>(param) ){
270  // G4cout << "G4PhantomParameterisation volume found " << (*cite)->GetName() << G4endl;
271  thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
272  }
273  }
274 
275  if( !thePhantomParam && mustExist )
276  G4Exception("G4EnergySplitter::GetPhantomParam",
277  "PhantomParamError", FatalException,
278  "No G4PhantomParameterisation found !");
279 }
280 
281 
282 //-----------------------------------------------------------------------
283 G4bool G4EnergySplitter::IsPhantomVolume( G4VPhysicalVolume* pv )
284 {
285  EAxis axis;
286  G4int nReplicas;
287  G4double width,offset;
288  G4bool consuming;
289  pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
290  EVolume type = (consuming) ? kReplica : kParameterised;
291  if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {
292  return TRUE;
293  } else {
294  return FALSE;
295  }
296 
297 }
298 
299 //-----------------------------------------------------------------------
301 {
302  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().begin())).first;
303 }
304 
305 //-----------------------------------------------------------------------
307 {
308  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().rbegin())).first;
309 }
310 
311 //-----------------------------------------------------------------------
312 void G4EnergySplitter::GetVoxelID( G4int stepNo, G4int& voxelID )
313 {
314  if( stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()) ) {
315  G4Exception("G4EnergySplitter::GetVoxelID",
316  "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
318  G4String("stepNo = " + G4UIcommand::ConvertToString(stepNo) + ", number of voxels = " + G4UIcommand::ConvertToString(G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size())) ).c_str());
319  }
320  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
321  advance( ite, stepNo );
322  voxelID = (*ite).first;
323 
324 }
325 
326 
327 //-----------------------------------------------------------------------
328 void G4EnergySplitter::GetStepLength( G4int stepNo, G4double& stepLength )
329 {
330  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
331  advance( ite, stepNo );
332  stepLength = (*ite).second;
333 }
G4ParticleDefinition * GetDefinition() const
G4double GetStepLength() const
const G4String & GetName() const
Definition: G4Material.hh:176
G4Material * GetMaterial() const
const std::vector< std::pair< G4int, G4double > > & GetStepLengths()
static G4String ConvertToString(G4bool boolVal)
Definition: G4UIcommand.cc:357
#define width
virtual ~G4EnergySplitter()
G4double TrueStepLength(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *part)
int G4int
Definition: G4Types.hh:78
Double_t edep
Definition: macro.C:13
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=0)
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
#define FALSE
Definition: globals.hh:52
TString part[npart]
Definition: Style.C:32
#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
G4int first
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:68
G4double GetKineticEnergy() const
double G4double
Definition: G4Types.hh:76
static G4RegularNavigationHelper * Instance()
G4Track * GetTrack() const
G4double GetPDGCharge() const
G4int SplitEnergyInVolumes(const G4Step *aStep)