| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 G4EnergyLossForExtrapolato
47 thePhantomParam = nullptr;
48 theNIterations = 2;
49 }
50
51 G4EnergySplitter::~G4EnergySplitter()
52 {
53 delete theElossExt;
54 }
55
56 G4int G4EnergySplitter::SplitEnergyInVolumes(c
57 {
58 theEnergies.clear();
59
60 if (aStep == nullptr) return false; // it i
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::SplitEnergyIn
72 << " Nsteps " << G4RegularNavigatio
73 << G4endl;
74 #endif
75 if (G4RegularNavigationHelper::Instance()->G
76 || aStep->GetTrack()->GetDefinition()->G
77 { // we are only counting dose deposit
78 return (G4int)theEnergies.size();
79 }
80 if (G4RegularNavigationHelper::Instance()->G
81 theEnergies.push_back(edep);
82 return (G4int)theEnergies.size();
83 }
84
85 if (thePhantomParam == nullptr) GetPhantomPa
86
87 //----- Distribute energy deposited in voxel
88 std::vector<std::pair<G4int, G4double>> rnsl
89 G4RegularNavigationHelper::Instance()->Get
90
91 const G4ParticleDefinition* part = aStep->Ge
92 G4double kinEnergyPreOrig = aStep->GetPreSte
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::SplitEnergy
104 << sl << G4endl;
105 #endif
106 }
107
108 #ifdef VERBOSE_ENERSPLIT
109 if (verbose)
110 G4cout << "G4EnergySplitter RN: step leng
111 << stepLength << " ratio " << stepL
112 << aStep->GetPreStepPoint()->GetKin
113 << aStep->GetPreStepPoint()->GetMat
114 << rnsl.size() << G4endl;
115 #endif
116 //----- No iterations to correct elost and m
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;
122
123 #ifdef VERBOSE_ENERSPLIT
124 if (verbose)
125 G4cout << "G4EnergySplitter::SplitEner
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::SplitEner
143 << " RN: iter0 corrected energy
144 }
145 }
146 #endif
147
148 G4double slRatio = stepLength / slSum;
149 #ifdef VERBOSE_ENERSPLIT
150 if (verbose)
151 G4cout << "G4EnergySplitter::SplitEnergy
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 <= theNIterati
160 //--- iter1: distribute true step length
161 // a constant so that the sum gives the
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::Split
169 << " corrected step length
170 #endif
171 }
172
173 for (ii = 0; ii < rnsl.size(); ++ii) {
174 const G4Material* mate = thePhantomP
175 G4double dEdx = 0.;
176 if (kinEnergyPre > 0.) { // t check
177 dEdx = emcalc.GetDEDX(kinEnergyPre
178 }
179 G4double elost = stepLengths[ii] * d
180
181 #ifdef VERBOSE_ENERSPLIT
182 if (verbose)
183 G4cout << "G4EnergySplitter::Split
184 << elost << " energy at int
185 << stepLengths[ii] << " * d
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
195 //-- Get ratios for each energy
196 slSum = 0.;
197 kinEnergyPre = kinEnergyPreOrig;
198 for (ii = 0; ii < rnsl.size(); ++ii) {
199 const G4Material* mate = thePhantomP
200 stepLengths[ii] = theElossExt->TrueS
201 kinEnergyPre -= theEnergies[ii];
202
203 #ifdef VERBOSE_ENERSPLIT
204 if (verbose)
205 G4cout << "G4EnergySplitter::Split
206 << " step length geom " <<
207 << rnsl[ii].second / stepLe
208 #endif
209
210 slSum += stepLengths[ii];
211 }
212
213 // Correct step lengths so that they s
214 G4double slratio = aStep->GetStepLengt
215 #ifdef VERBOSE_ENERSPLIT
216 if (verbose)
217 G4cout << "G4EnergySplitter::SplitEn
218 << " step ratio " << slRatio
219 #endif
220 for (ii = 0; ii < rnsl.size(); ++ii) {
221 stepLengths[ii] *= slratio;
222 #ifdef VERBOSE_ENERSPLIT
223 if (verbose)
224 G4cout << "G4EnergySplitter::Split
225 << " corrected step length
226 #endif
227 }
228
229 //---- Recalculate energy lost with th
230 kinEnergyPre = aStep->GetPreStepPoint(
231 totalELost = 0.;
232 for (ii = 0; ii < rnsl.size(); ++ii) {
233 const G4Material* mate = thePhantomP
234 G4double dEdx = 0.;
235 if (kinEnergyPre > 0.) {
236 dEdx = emcalc.GetDEDX(kinEnergyPre
237 }
238 G4double elost = stepLengths[ii] * d
239 #ifdef VERBOSE_ENERSPLIT
240 if (verbose)
241 G4cout << "G4EnergySplitter::Split
242 << " energy lost " << elost
243 << " = stepLength " << step
244 #endif
245 kinEnergyPre -= elost;
246 theEnergies[ii] = elost;
247 totalELost += elost;
248 }
249 }
250
251 // correct energies so that they reprodu
252 G4double enerRatio = (edep / totalELost)
253
254 #ifdef VERBOSE_ENERSPLIT
255 if (verbose)
256 G4cout << "G4EnergySplitter::SplitEner
257 << " energy ratio " << enerRati
258 #endif
259
260 #ifdef VERBOSE_ENERSPLIT
261 G4double elostTot = 0.;
262 #endif
263 for (ii = 0; ii < theEnergies.size(); ++
264 theEnergies[ii] *= enerRatio;
265 #ifdef VERBOSE_ENERSPLIT
266 elostTot += theEnergies[ii];
267 if (verbose)
268 G4cout << "G4EnergySplitter::SplitEn
269 << " corrected energy lost "
270 << theEnergies[ii] / enerRati
271 << kinEnergyPreOrig - elostTo
272 << kinEnergyPreOrig - elostTo
273 #endif
274 }
275 }
276 }
277
278 return (G4int)theEnergies.size();
279 }
280
281 //--------------------------------------------
282 void G4EnergySplitter::GetPhantomParam(G4bool
283 {
284 G4PhysicalVolumeStore* pvs = G4PhysicalVolum
285 for (const auto pv : *pvs) {
286 if (IsPhantomVolume(pv)) {
287 const auto pvparam = static_cast<const G
288 G4VPVParameterisation* param = pvparam->
289 thePhantomParam = static_cast<G4PhantomP
290 }
291 }
292
293 if ((thePhantomParam == nullptr) && mustExis
294 G4Exception("G4EnergySplitter::GetPhantomP
295 "No G4PhantomParameterisation
296 }
297
298 //--------------------------------------------
299 G4bool G4EnergySplitter::IsPhantomVolume(G4VPh
300 {
301 EAxis axis;
302 G4int nReplicas;
303 G4double width, offset;
304 G4bool consuming;
305 pv->GetReplicationData(axis, nReplicas, widt
306 EVolume type = (consuming) ? kReplica : kPar
307 return type == kParameterised && pv->GetRegu
308 }
309
310 //--------------------------------------------
311 void G4EnergySplitter::GetLastVoxelID(G4int& v
312 {
313 voxelID = (*(G4RegularNavigationHelper::Inst
314 }
315
316 //--------------------------------------------
317 void G4EnergySplitter::GetFirstVoxelID(G4int&
318 {
319 voxelID = (*(G4RegularNavigationHelper::Inst
320 }
321
322 //--------------------------------------------
323 void G4EnergySplitter::GetVoxelID(G4int stepNo
324 {
325 if (stepNo < 0 || stepNo >= G4int(G4RegularN
326 {
327 G4Exception("G4EnergySplitter::GetVoxelID"
328 "Invalid stepNo, smaller than
329 FatalErrorInArgument,
330 G4String("stepNo = " + G4UIcom
331 + ", number of voxels
332 + G4UIcommand::Conver
333 G4int(G4RegularNavi
334 .c_str());
335 }
336 auto ite = G4RegularNavigationHelper::Instan
337 advance(ite, stepNo);
338 voxelID = (*ite).first;
339 }
340
341 //--------------------------------------------
342 void G4EnergySplitter::GetStepLength(G4int ste
343 {
344 auto ite = G4RegularNavigationHelper::Instan
345 advance(ite, stepNo);
346 stepLength = (*ite).second;
347 }
348