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