| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4EnergyLossTables.cc,v 1.14.2.1 2001/06/28 19:12:48 gunter Exp $
>> 25 // GEANT4 tag $Name: $
27 // 26 //
28 // ------------------------------------------- << 27 // -------------------------------------------------------------------
29 // first version created by P.Urban , 06/04/19 28 // first version created by P.Urban , 06/04/1998
30 // modifications + "precise" functions added b 29 // modifications + "precise" functions added by L.Urban , 27/05/98
31 // modifications , TOF functions , 26/10/98, L 30 // modifications , TOF functions , 26/10/98, L.Urban
32 // cache mechanism in order to gain time, 11/0 31 // cache mechanism in order to gain time, 11/02/99, L.Urban
33 // bug fixed , 12/04/99 , L.Urban 32 // bug fixed , 12/04/99 , L.Urban
34 // 10.11.99: moved from RWT hash dictionary to << 33 // 10/11/99: moved from RWT hash dictionary to STL map, G.Barrand, M.Maire
35 // 27.09.01 L.Urban , bug fixed (negative ener <<
36 // 26.10.01 all static functions moved from .i <<
37 // 15.01.03 Add interfaces required for "cut p <<
38 // 12.03.03 Add warnings to obsolete interface <<
39 // 10.04.03 Add call to G4LossTableManager is <<
40 // <<
41 // ------------------------------------------- 34 // -------------------------------------------------------------------
42 35
43 #include "G4EnergyLossTables.hh" 36 #include "G4EnergyLossTables.hh"
44 #include "G4SystemOfUnits.hh" <<
45 #include "G4MaterialCutsCouple.hh" <<
46 #include "G4RegionStore.hh" <<
47 #include "G4LossTableManager.hh" <<
48 <<
49 37
50 //....oooOO0OOooo........oooOO0OOooo........oo 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 39
52 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 40 G4EnergyLossTablesHelper G4EnergyLossTables::t ;
53 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 41 const G4ParticleDefinition* G4EnergyLossTables::lastParticle = NULL ;
54 G4ParticleDefinition* G4EnergyLossTables::last << 42 G4double G4EnergyLossTables::QQPositron = eplus*eplus ;
55 G4double G4EnergyLossTables::QQPositron = 1.0; <<
56 G4double G4EnergyLossTables::Chargesquare ; 43 G4double G4EnergyLossTables::Chargesquare ;
57 G4int G4EnergyLossTables::oldIndex = -1 ; 44 G4int G4EnergyLossTables::oldIndex = -1 ;
58 G4double G4EnergyLossTables::rmin = 0. ; 45 G4double G4EnergyLossTables::rmin = 0. ;
59 G4double G4EnergyLossTables::rmax = 0. ; 46 G4double G4EnergyLossTables::rmax = 0. ;
60 G4double G4EnergyLossTables::Thigh = 0. ; 47 G4double G4EnergyLossTables::Thigh = 0. ;
61 G4int G4EnergyLossTables::let_counter = 0; <<
62 G4int G4EnergyLossTables::let_max_num_warni <<
63 G4bool G4EnergyLossTables::first_loss = true <<
64 <<
65 G4EnergyLossTables::helper_map *G4EnergyLossTa <<
66 <<
67 //....oooOO0OOooo........oooOO0OOooo........oo <<
68 <<
69 G4EnergyLossTablesHelper::G4EnergyLossTablesHe <<
70 const G4PhysicsTable* aDEDXTable, <<
71 const G4PhysicsTable* aRangeTable, <<
72 const G4PhysicsTable* anInverseRangeTable, <<
73 const G4PhysicsTable* aLabTimeTable, <<
74 const G4PhysicsTable* aProperTimeTable, <<
75 G4double aLowestKineticEnergy, <<
76 G4double aHighestKineticEnergy, <<
77 G4double aMassRatio, <<
78 G4int aNumberOfBins) <<
79 : <<
80 theDEDXTable(aDEDXTable), theRangeTable(aRan <<
81 theInverseRangeTable(anInverseRangeTable), <<
82 theLabTimeTable(aLabTimeTable), <<
83 theProperTimeTable(aProperTimeTable), <<
84 theLowestKineticEnergy(aLowestKineticEnergy) <<
85 theHighestKineticEnergy(aHighestKineticEnerg <<
86 theMassRatio(aMassRatio), <<
87 theNumberOfBins(aNumberOfBins) <<
88 { } <<
89 48
90 //....oooOO0OOooo........oooOO0OOooo........oo << 49 G4EnergyLossTables::helper_map G4EnergyLossTables::dict;
91 <<
92 G4EnergyLossTablesHelper::G4EnergyLossTablesHe <<
93 { <<
94 theLowestKineticEnergy = 0.0; <<
95 theHighestKineticEnergy= 0.0; <<
96 theMassRatio = 0.0; <<
97 theNumberOfBins = 0; <<
98 theDEDXTable = theRangeTable = theInverseRan <<
99 = theProperTimeTable = nullptr; <<
100 } <<
101 50
102 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 52
104 void G4EnergyLossTables::Register( 53 void G4EnergyLossTables::Register(
105 const G4ParticleDefinition* p, 54 const G4ParticleDefinition* p,
106 const G4PhysicsTable* tDEDX, 55 const G4PhysicsTable* tDEDX,
107 const G4PhysicsTable* tRange, 56 const G4PhysicsTable* tRange,
108 const G4PhysicsTable* tInverseRange, 57 const G4PhysicsTable* tInverseRange,
109 const G4PhysicsTable* tLabTime, 58 const G4PhysicsTable* tLabTime,
110 const G4PhysicsTable* tProperTime, 59 const G4PhysicsTable* tProperTime,
111 G4double lowestKineticEnergy, 60 G4double lowestKineticEnergy,
112 G4double highestKineticEnergy, 61 G4double highestKineticEnergy,
113 G4double massRatio, 62 G4double massRatio,
114 G4int NumberOfBins) 63 G4int NumberOfBins)
115 { 64 {
116 if (!dict) dict = new G4EnergyLossTables::he << 65 dict[p]= G4EnergyLossTablesHelper(tDEDX, tRange,tInverseRange,
117 if (!null_loss) null_loss = new G4EnergyLoss <<
118 if (!t) t = new G4EnergyLossTablesHelper; <<
119 <<
120 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, <<
121 tLabTime,tProperTime,lowes 66 tLabTime,tProperTime,lowestKineticEnergy,
122 highestKineticEnergy, massRatio,Number 67 highestKineticEnergy, massRatio,NumberOfBins);
123 << 68
124 *t = GetTables(p) ; // important for cach << 69 t = GetTables(p) ; // important for cache !!!!!
125 lastParticle = (G4ParticleDefinition*) p ; << 70 lastParticle = p ;
126 Chargesquare = (p->GetPDGCharge())*(p->GetPD 71 Chargesquare = (p->GetPDGCharge())*(p->GetPDGCharge())/
127 QQPositron ; 72 QQPositron ;
128 if (first_loss ) { <<
129 *null_loss = G4EnergyLossTablesHelper( <<
130 nullptr, nullptr, nullptr, nu <<
131 first_loss = false; <<
132 } <<
133 } <<
134 <<
135 //....oooOO0OOooo........oooOO0OOooo........oo <<
136 <<
137 const G4PhysicsTable* G4EnergyLossTables::GetD <<
138 const G4ParticleDefinition* p) <<
139 { <<
140 if (!dict) dict = new G4EnergyLossTables::he <<
141 helper_map::iterator it; <<
142 if((it=dict->find(p))==dict->end()) return n <<
143 return (*it).second.theDEDXTable; <<
144 } <<
145 <<
146 //....oooOO0OOooo........oooOO0OOooo........oo <<
147 <<
148 const G4PhysicsTable* G4EnergyLossTables::GetR <<
149 const G4ParticleDefinition* p) <<
150 { <<
151 if (!dict) dict = new G4EnergyLossTables::he <<
152 helper_map::iterator it; <<
153 if((it=dict->find(p))==dict->end()) return n <<
154 return (*it).second.theRangeTable; <<
155 } <<
156 <<
157 //....oooOO0OOooo........oooOO0OOooo........oo <<
158 <<
159 const G4PhysicsTable* G4EnergyLossTables::GetI <<
160 const G4ParticleDefinition* p) <<
161 { <<
162 if (!dict) dict = new G4EnergyLossTables::he <<
163 helper_map::iterator it; <<
164 if((it=dict->find(p))==dict->end()) return n <<
165 return (*it).second.theInverseRangeTable; <<
166 } <<
167 <<
168 //....oooOO0OOooo........oooOO0OOooo........oo <<
169 <<
170 const G4PhysicsTable* G4EnergyLossTables::GetL <<
171 const G4ParticleDefinition* p) <<
172 { <<
173 if (!dict) dict = new G4EnergyLossTables::he <<
174 helper_map::iterator it; <<
175 if((it=dict->find(p))==dict->end()) return n <<
176 return (*it).second.theLabTimeTable; <<
177 } <<
178 <<
179 //....oooOO0OOooo........oooOO0OOooo........oo <<
180 <<
181 const G4PhysicsTable* G4EnergyLossTables::GetP <<
182 const G4ParticleDefinition* p) <<
183 { <<
184 if (!dict) dict = new G4EnergyLossTables::he <<
185 helper_map::iterator it; <<
186 if((it=dict->find(p))==dict->end()) return n <<
187 return (*it).second.theProperTimeTable; <<
188 } <<
189 <<
190 //....oooOO0OOooo........oooOO0OOooo........oo <<
191 <<
192 G4EnergyLossTablesHelper G4EnergyLossTables::G <<
193 const G4ParticleDefinition* p) <<
194 { <<
195 if (!dict) dict = new G4EnergyLossTables::he <<
196 if (!null_loss) null_loss = new G4EnergyLoss <<
197 <<
198 helper_map::iterator it; <<
199 if ((it=dict->find(p))==dict->end()) { <<
200 return *null_loss; <<
201 } <<
202 return (*it).second; <<
203 } <<
204 <<
205 //....oooOO0OOooo........oooOO0OOooo........oo <<
206 <<
207 G4double G4EnergyLossTables::GetDEDX( <<
208 const G4ParticleDefinition *aParticle, <<
209 G4double KineticEnergy, <<
210 const G4Material *aMaterial) <<
211 { <<
212 if (!t) t = new G4EnergyLossTablesHelper; <<
213 <<
214 CPRWarning(); <<
215 if(aParticle != (const G4ParticleDefinition* <<
216 { <<
217 *t= GetTables(aParticle); <<
218 lastParticle = (G4ParticleDefinition*) aPa <<
219 Chargesquare = (aParticle->GetPDGCharge()) <<
220 (aParticle->GetPDGCharge()) <<
221 QQPositron ; <<
222 oldIndex = -1 ; <<
223 } <<
224 const G4PhysicsTable* dEdxTable= t->theDEDX <<
225 if (!dEdxTable) { <<
226 ParticleHaveNoLoss(aParticle,"dEdx"); <<
227 return 0.0; <<
228 } <<
229 <<
230 G4int materialIndex = (G4int)aMaterial->GetI <<
231 G4double scaledKineticEnergy = KineticEnergy <<
232 G4double dEdx; <<
233 G4bool isOut; <<
234 <<
235 if (scaledKineticEnergy<t->theLowestKineticE <<
236 <<
237 dEdx =(*dEdxTable)(materialIndex)->GetVal <<
238 t->theLowestKineticEnergy,isOut) <<
239 *std::sqrt(scaledKineticEnergy/t->t <<
240 <<
241 } else if (scaledKineticEnergy>t->theHighest <<
242 <<
243 dEdx = (*dEdxTable)(materialIndex)->GetVa <<
244 t->theHighestKineticEnergy,isOut); <<
245 <<
246 } else { <<
247 <<
248 dEdx = (*dEdxTable)(materialIndex)->GetVal <<
249 scaledKineticEnergy,isOut); <<
250 <<
251 } <<
252 <<
253 return dEdx*Chargesquare; <<
254 } <<
255 <<
256 //....oooOO0OOooo........oooOO0OOooo........oo <<
257 <<
258 G4double G4EnergyLossTables::GetLabTime( <<
259 const G4ParticleDefinition *aParticle, <<
260 G4double KineticEnergy, <<
261 const G4Material *aMaterial) <<
262 { <<
263 if (!t) t = new G4EnergyLossTablesHelper; <<
264 <<
265 CPRWarning(); <<
266 if(aParticle != (const G4ParticleDefinition* <<
267 { <<
268 *t= GetTables(aParticle); <<
269 lastParticle = (G4ParticleDefinition*) aPa <<
270 oldIndex = -1 ; <<
271 } <<
272 const G4PhysicsTable* labtimeTable= t->theLa <<
273 if (!labtimeTable) { <<
274 ParticleHaveNoLoss(aParticle,"LabTime"); <<
275 return 0.0; <<
276 } <<
277 <<
278 const G4double parlowen=0.4 , ppar=0.5-parlo <<
279 G4int materialIndex = (G4int)aMaterial->GetI <<
280 G4double scaledKineticEnergy = KineticEnergy <<
281 G4double time; <<
282 G4bool isOut; <<
283 <<
284 if (scaledKineticEnergy<t->theLowestKineticE <<
285 <<
286 time = std::exp(ppar*std::log(scaledKinet <<
287 (*labtimeTable)(materialIndex)->Ge <<
288 t->theLowestKineticEnergy,isOut) <<
289 <<
290 <<
291 } else if (scaledKineticEnergy>t->theHighest <<
292 <<
293 time = (*labtimeTable)(materialIndex)->Ge <<
294 t->theHighestKineticEnergy,isOut <<
295 <<
296 } else { <<
297 <<
298 time = (*labtimeTable)(materialIndex)->Get <<
299 scaledKineticEnergy,isOut); <<
300 <<
301 } <<
302 <<
303 return time/t->theMassRatio ; <<
304 } <<
305 <<
306 //....oooOO0OOooo........oooOO0OOooo........oo <<
307 <<
308 G4double G4EnergyLossTables::GetDeltaLabTime( <<
309 const G4ParticleDefinition *aParticle, <<
310 G4double KineticEnergyStart, <<
311 G4double KineticEnergyEnd, <<
312 const G4Material *aMaterial) <<
313 { <<
314 if (!t) t = new G4EnergyLossTablesHelper; <<
315 <<
316 CPRWarning(); <<
317 if(aParticle != (const G4ParticleDefinition* <<
318 { <<
319 *t= GetTables(aParticle); <<
320 lastParticle = (G4ParticleDefinition*) aPa <<
321 oldIndex = -1 ; <<
322 } <<
323 const G4PhysicsTable* labtimeTable= t->theLa <<
324 if (!labtimeTable) { <<
325 ParticleHaveNoLoss(aParticle,"LabTime"); <<
326 return 0.0; <<
327 } <<
328 <<
329 const G4double parlowen=0.4 , ppar=0.5-parlo <<
330 const G4double dToverT = 0.05 , facT = 1. -d <<
331 G4double timestart,timeend,deltatime,dTT; <<
332 G4bool isOut; <<
333 <<
334 G4int materialIndex = (G4int)aMaterial->GetI <<
335 G4double scaledKineticEnergy = KineticEnergy <<
336 <<
337 if (scaledKineticEnergy<t->theLowestKineticE <<
338 <<
339 timestart = std::exp(ppar*std::log(scaled <<
340 (*labtimeTable)(materialIndex) <<
341 t->theLowestKineticEnergy,isOu <<
342 <<
343 <<
344 } else if (scaledKineticEnergy>t->theHighest <<
345 <<
346 timestart = (*labtimeTable)(materialIndex <<
347 t->theHighestKineticEnergy,isO <<
348 <<
349 } else { <<
350 <<
351 timestart = (*labtimeTable)(materialIndex) <<
352 scaledKineticEnergy,isOut); <<
353 <<
354 } <<
355 <<
356 dTT = (KineticEnergyStart - KineticEnergyEnd <<
357 <<
358 if( dTT < dToverT ) <<
359 scaledKineticEnergy = facT*KineticEnergySt <<
360 else <<
361 scaledKineticEnergy = KineticEnergyEnd*t-> <<
362 <<
363 if (scaledKineticEnergy<t->theLowestKineticE <<
364 <<
365 timeend = std::exp(ppar*std::log(scaledKi <<
366 (*labtimeTable)(materialIndex) <<
367 t->theLowestKineticEnergy,isOu <<
368 <<
369 <<
370 } else if (scaledKineticEnergy>t->theHighest <<
371 <<
372 timeend = (*labtimeTable)(materialIndex)- <<
373 t->theHighestKineticEnergy,isO <<
374 <<
375 } else { <<
376 <<
377 timeend = (*labtimeTable)(materialIndex)-> <<
378 scaledKineticEnergy,isOut); <<
379 <<
380 } <<
381 <<
382 deltatime = timestart - timeend ; <<
383 <<
384 if( dTT < dToverT ) <<
385 deltatime *= dTT/dToverT; <<
386 <<
387 return deltatime/t->theMassRatio ; <<
388 } <<
389 <<
390 //....oooOO0OOooo........oooOO0OOooo........oo <<
391 <<
392 G4double G4EnergyLossTables::GetProperTime( <<
393 const G4ParticleDefinition *aParticle, <<
394 G4double KineticEnergy, <<
395 const G4Material *aMaterial) <<
396 { <<
397 if (!t) t = new G4EnergyLossTablesHelper; <<
398 <<
399 CPRWarning(); <<
400 if(aParticle != (const G4ParticleDefinition* <<
401 { <<
402 *t= GetTables(aParticle); <<
403 lastParticle = (G4ParticleDefinition*) aPa <<
404 oldIndex = -1 ; <<
405 } <<
406 const G4PhysicsTable* propertimeTable= t->th <<
407 if (!propertimeTable) { <<
408 ParticleHaveNoLoss(aParticle,"ProperTime") <<
409 return 0.0; <<
410 } <<
411 <<
412 const G4double parlowen=0.4 , ppar=0.5-parlo <<
413 G4int materialIndex = (G4int)aMaterial->GetI <<
414 G4double scaledKineticEnergy = KineticEnergy <<
415 G4double time; <<
416 G4bool isOut; <<
417 <<
418 if (scaledKineticEnergy<t->theLowestKineticE <<
419 <<
420 time = std::exp(ppar*std::log(scaledKinet <<
421 (*propertimeTable)(materialIndex)- <<
422 t->theLowestKineticEnergy,isOut) <<
423 <<
424 <<
425 } else if (scaledKineticEnergy>t->theHighest <<
426 <<
427 time = (*propertimeTable)(materialIndex)- <<
428 t->theHighestKineticEnergy,isOut <<
429 <<
430 } else { <<
431 <<
432 time = (*propertimeTable)(materialIndex)-> <<
433 scaledKineticEnergy,isOut); <<
434 <<
435 } <<
436 <<
437 return time/t->theMassRatio ; <<
438 } <<
439 <<
440 //....oooOO0OOooo........oooOO0OOooo........oo <<
441 <<
442 G4double G4EnergyLossTables::GetDeltaProperTim <<
443 const G4ParticleDefinition *aParticle, <<
444 G4double KineticEnergyStart, <<
445 G4double KineticEnergyEnd, <<
446 const G4Material *aMaterial) <<
447 { <<
448 if (!t) t = new G4EnergyLossTablesHelper; <<
449 <<
450 CPRWarning(); <<
451 if(aParticle != (const G4ParticleDefinition* <<
452 { <<
453 *t= GetTables(aParticle); <<
454 lastParticle = (G4ParticleDefinition*) aPa <<
455 oldIndex = -1 ; <<
456 } <<
457 const G4PhysicsTable* propertimeTable= t->th <<
458 if (!propertimeTable) { <<
459 ParticleHaveNoLoss(aParticle,"ProperTime") <<
460 return 0.0; <<
461 } <<
462 <<
463 const G4double parlowen=0.4 , ppar=0.5-parlo <<
464 const G4double dToverT = 0.05 , facT = 1. -d <<
465 G4double timestart,timeend,deltatime,dTT; <<
466 G4bool isOut; <<
467 <<
468 G4int materialIndex = (G4int)aMaterial->GetI <<
469 G4double scaledKineticEnergy = KineticEnergy <<
470 <<
471 if (scaledKineticEnergy<t->theLowestKineticE <<
472 <<
473 timestart = std::exp(ppar*std::log(scaled <<
474 (*propertimeTable)(materialInd <<
475 t->theLowestKineticEnergy,isOu <<
476 <<
477 <<
478 } else if (scaledKineticEnergy>t->theHighest <<
479 <<
480 timestart = (*propertimeTable)(materialIn <<
481 t->theHighestKineticEnergy,isO <<
482 <<
483 } else { <<
484 <<
485 timestart = (*propertimeTable)(materialInd <<
486 scaledKineticEnergy,isOut); <<
487 <<
488 } <<
489 <<
490 dTT = (KineticEnergyStart - KineticEnergyEnd <<
491 <<
492 if( dTT < dToverT ) <<
493 scaledKineticEnergy = facT*KineticEnergySt <<
494 else <<
495 scaledKineticEnergy = KineticEnergyEnd*t-> <<
496 <<
497 if (scaledKineticEnergy<t->theLowestKineticE <<
498 <<
499 timeend = std::exp(ppar*std::log(scaledKi <<
500 (*propertimeTable)(materialInd <<
501 t->theLowestKineticEnergy,isOu <<
502 <<
503 <<
504 } else if (scaledKineticEnergy>t->theHighest <<
505 <<
506 timeend = (*propertimeTable)(materialInde <<
507 t->theHighestKineticEnergy,isO <<
508 <<
509 } else { <<
510 <<
511 timeend = (*propertimeTable)(materialIndex <<
512 scaledKineticEnergy,isOut); <<
513 <<
514 } <<
515 <<
516 deltatime = timestart - timeend ; <<
517 <<
518 if( dTT < dToverT ) <<
519 deltatime *= dTT/dToverT ; <<
520 <<
521 return deltatime/t->theMassRatio ; <<
522 } <<
523 <<
524 //....oooOO0OOooo........oooOO0OOooo........oo <<
525 <<
526 G4double G4EnergyLossTables::GetRange( <<
527 const G4ParticleDefinition *aParticle, <<
528 G4double KineticEnergy, <<
529 const G4Material *aMaterial) <<
530 { <<
531 if (!t) t = new G4EnergyLossTablesHelper; <<
532 <<
533 CPRWarning(); <<
534 if(aParticle != (const G4ParticleDefinition* <<
535 { <<
536 *t= GetTables(aParticle); <<
537 lastParticle = (G4ParticleDefinition*) aPa <<
538 Chargesquare = (aParticle->GetPDGCharge()) <<
539 (aParticle->GetPDGCharge()) <<
540 QQPositron ; <<
541 oldIndex = -1 ; <<
542 } <<
543 const G4PhysicsTable* rangeTable= t->theRang <<
544 const G4PhysicsTable* dEdxTable= t->theDEDX <<
545 if (!rangeTable) { <<
546 ParticleHaveNoLoss(aParticle,"Range"); <<
547 return 0.0; <<
548 } <<
549 <<
550 G4int materialIndex = (G4int)aMaterial->GetI <<
551 G4double scaledKineticEnergy = KineticEnergy <<
552 G4double Range; <<
553 G4bool isOut; <<
554 <<
555 if (scaledKineticEnergy<t->theLowestKineticE <<
556 <<
557 Range = std::sqrt(scaledKineticEnergy/t->t <<
558 (*rangeTable)(materialIndex)->GetV <<
559 t->theLowestKineticEnergy,isOut) <<
560 <<
561 } else if (scaledKineticEnergy>t->theHighest <<
562 <<
563 Range = (*rangeTable)(materialIndex)->GetV <<
564 t->theHighestKineticEnergy,isOut)+ <<
565 (scaledKineticEnergy-t->theHighest <<
566 (*dEdxTable)(materialIndex)->GetVa <<
567 t->theHighestKineticEnergy,isOut <<
568 <<
569 } else { <<
570 <<
571 Range = (*rangeTable)(materialIndex)->GetV <<
572 scaledKineticEnergy,isOut); <<
573 <<
574 } <<
575 <<
576 return Range/(Chargesquare*t->theMassRatio); <<
577 } <<
578 <<
579 //....oooOO0OOooo........oooOO0OOooo........oo <<
580 <<
581 G4double G4EnergyLossTables::GetPreciseEnergyF <<
582 const G4P <<
583 G4d <<
584 const G4M <<
585 // it returns the value of the kinetic energy <<
586 { <<
587 if (!t) t = new G4EnergyLossTablesHelper; <<
588 <<
589 CPRWarning(); <<
590 if( aParticle != (const G4ParticleDefinition <<
591 { <<
592 *t= GetTables(aParticle); <<
593 lastParticle = (G4ParticleDefinition*) aPa <<
594 Chargesquare = (aParticle->GetPDGCharge()) <<
595 (aParticle->GetPDGCharge()) <<
596 QQPositron ; <<
597 oldIndex = -1 ; <<
598 } <<
599 const G4PhysicsTable* dEdxTable= t->theDEDX <<
600 const G4PhysicsTable* inverseRangeTable= t- <<
601 if (!inverseRangeTable) { <<
602 ParticleHaveNoLoss(aParticle,"InverseRange <<
603 return 0.0; <<
604 } <<
605 <<
606 G4double scaledrange,scaledKineticEnergy ; <<
607 G4bool isOut ; <<
608 <<
609 G4int materialIndex = (G4int)aMaterial->GetI <<
610 <<
611 if(materialIndex != oldIndex) <<
612 { <<
613 oldIndex = materialIndex ; <<
614 rmin = (*inverseRangeTable)(materialIndex) <<
615 GetLowEdgeEnergy <<
616 rmax = (*inverseRangeTable)(materialIndex) <<
617 GetLowEdgeEnergy(t->theNumb <<
618 Thigh = (*inverseRangeTable)(materialIndex <<
619 GetValue(rmax,is <<
620 } <<
621 <<
622 scaledrange = range*Chargesquare*t->theMassR <<
623 <<
624 if(scaledrange < rmin) <<
625 { <<
626 scaledKineticEnergy = t->theLowestKineticE <<
627 scaledrange*scaledrange/(rm <<
628 } <<
629 else <<
630 { <<
631 if(scaledrange < rmax) <<
632 { <<
633 scaledKineticEnergy = (*inverseRangeTabl <<
634 GetValue( scaled <<
635 } <<
636 else <<
637 { <<
638 scaledKineticEnergy = Thigh + <<
639 (scaledrange-rmax)* <<
640 (*dEdxTable)(materialInd <<
641 GetValue(Thig <<
642 } <<
643 } <<
644 <<
645 return scaledKineticEnergy/t->theMassRatio ; <<
646 } 73 }
647 74
648 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
649 76
650 G4double G4EnergyLossTables::GetPreciseDEDX( 77 G4double G4EnergyLossTables::GetPreciseDEDX(
651 const G4ParticleDefinition *aParticle, 78 const G4ParticleDefinition *aParticle,
652 G4double KineticEnergy, 79 G4double KineticEnergy,
653 const G4Material *aMaterial) 80 const G4Material *aMaterial)
654 { 81 {
655 if (!t) t = new G4EnergyLossTablesHelper; << 82 if( aParticle != lastParticle)
656 <<
657 CPRWarning(); <<
658 if( aParticle != (const G4ParticleDefinition <<
659 { 83 {
660 *t= GetTables(aParticle); << 84 t= GetTables(aParticle);
661 lastParticle = (G4ParticleDefinition*) aPa << 85 lastParticle = aParticle;
662 Chargesquare = (aParticle->GetPDGCharge()) 86 Chargesquare = (aParticle->GetPDGCharge())*
663 (aParticle->GetPDGCharge()) 87 (aParticle->GetPDGCharge())/
664 QQPositron ; 88 QQPositron ;
665 oldIndex = -1 ; <<
666 } <<
667 const G4PhysicsTable* dEdxTable= t->theDEDX <<
668 if (!dEdxTable) { <<
669 ParticleHaveNoLoss(aParticle,"dEdx"); <<
670 return 0.0; <<
671 } 89 }
>> 90 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
672 91
673 G4int materialIndex = (G4int)aMaterial->GetI << 92 G4int materialIndex = aMaterial->GetIndex();
674 G4double scaledKineticEnergy = KineticEnergy << 93 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
675 G4double dEdx; 94 G4double dEdx;
676 G4bool isOut; 95 G4bool isOut;
677 96
678 if (scaledKineticEnergy<t->theLowestKineticE << 97 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
679 98
680 dEdx = std::sqrt(scaledKineticEnergy/t->t << 99 dEdx = sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)
681 *(*dEdxTable)(materialIndex)->GetV 100 *(*dEdxTable)(materialIndex)->GetValue(
682 t->theLowestKineticEnergy,isOut) << 101 t.theLowestKineticEnergy,isOut);
683 102
684 } else if (scaledKineticEnergy>t->theHighest << 103 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
685 104
686 dEdx = (*dEdxTable)(materialIndex)->GetVa 105 dEdx = (*dEdxTable)(materialIndex)->GetValue(
687 t->theHighestKineticEnergy,isOut); << 106 t.theHighestKineticEnergy,isOut);
688 107
689 } else { 108 } else {
690 << 109
691 dEdx = (*dEdxTable)(materialIndex)->GetV 110 dEdx = (*dEdxTable)(materialIndex)->GetValue(
692 scaledKineticEnergy, 111 scaledKineticEnergy,isOut) ;
693 112
694 } 113 }
695 114
696 return dEdx*Chargesquare; 115 return dEdx*Chargesquare;
697 } 116 }
698 117
699 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
700 119
701 G4double G4EnergyLossTables::GetPreciseRangeF 120 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy(
702 const G4ParticleDefinition *aParticle, 121 const G4ParticleDefinition *aParticle,
703 G4double KineticEnergy, 122 G4double KineticEnergy,
704 const G4Material *aMaterial) 123 const G4Material *aMaterial)
705 { 124 {
706 if (!t) t = new G4EnergyLossTablesHelper; << 125 if( aParticle != lastParticle)
707 <<
708 CPRWarning(); <<
709 if( aParticle != (const G4ParticleDefinition <<
710 { 126 {
711 *t= GetTables(aParticle); << 127 t= GetTables(aParticle);
712 lastParticle = (G4ParticleDefinition*) aPa << 128 lastParticle = aParticle;
713 Chargesquare = (aParticle->GetPDGCharge()) 129 Chargesquare = (aParticle->GetPDGCharge())*
714 (aParticle->GetPDGCharge()) 130 (aParticle->GetPDGCharge())/
715 QQPositron ; 131 QQPositron ;
716 oldIndex = -1 ; <<
717 } <<
718 const G4PhysicsTable* rangeTable= t->theRang <<
719 const G4PhysicsTable* dEdxTable= t->theDEDX <<
720 if (!rangeTable) { <<
721 ParticleHaveNoLoss(aParticle,"Range"); <<
722 return 0.0; <<
723 } 132 }
724 G4int materialIndex = (G4int)aMaterial->GetI << 133 const G4PhysicsTable* rangeTable= t.theRangeTable;
>> 134 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
725 135
726 G4double Thighr = t->theHighestKineticEnergy << 136 G4int materialIndex = aMaterial->GetIndex();
727 (*rangeTable)(materialIndex <<
728 GetLowEdgeEnergy(1) ; <<
729 137
730 G4double scaledKineticEnergy = KineticEnergy << 138 G4double Thighr = t.theHighestKineticEnergy*t.theLowestKineticEnergy/
731 G4double Range; <<
732 G4bool isOut; <<
733 <<
734 if (scaledKineticEnergy<t->theLowestKineticE <<
735 <<
736 Range = std::sqrt(scaledKineticEnergy/t->t <<
737 (*rangeTable)(materialIndex)->GetV <<
738 t->theLowestKineticEnergy,isOut) <<
739 <<
740 } else if (scaledKineticEnergy>Thighr) { <<
741 <<
742 Range = (*rangeTable)(materialIndex)->GetV <<
743 Thighr,isOut)+ <<
744 (scaledKineticEnergy-Thighr)/ <<
745 (*dEdxTable)(materialIndex)->GetVa <<
746 Thighr,isOut); <<
747 <<
748 } else { <<
749 <<
750 Range = (*rangeTable)(materialIndex)->Get <<
751 scaledKineticEnergy,isO <<
752 <<
753 } <<
754 <<
755 return Range/(Chargesquare*t->theMassRatio); <<
756 } <<
757 <<
758 //....oooOO0OOooo........oooOO0OOooo........oo <<
759 <<
760 G4double G4EnergyLossTables::GetDEDX( <<
761 const G4ParticleDefinition *aParticle, <<
762 G4double KineticEnergy, <<
763 const G4MaterialCutsCouple *couple, <<
764 G4bool check) <<
765 { <<
766 if (!t) t = new G4EnergyLossTablesHelper; <<
767 <<
768 if(aParticle != (const G4ParticleDefinition* <<
769 { <<
770 *t= GetTables(aParticle); <<
771 lastParticle = (G4ParticleDefinition*) aPa <<
772 Chargesquare = (aParticle->GetPDGCharge()) <<
773 (aParticle->GetPDGCharge()) <<
774 QQPositron ; <<
775 oldIndex = -1 ; <<
776 } <<
777 const G4PhysicsTable* dEdxTable= t->theDEDX <<
778 <<
779 if (!dEdxTable ) { <<
780 if (check) return G4LossTableManager::Inst <<
781 else ParticleHaveNoLoss(aParticle, " <<
782 return 0.0; <<
783 } <<
784 <<
785 G4int materialIndex = couple->GetIndex(); <<
786 G4double scaledKineticEnergy = KineticEnergy <<
787 G4double dEdx; <<
788 G4bool isOut; <<
789 <<
790 if (scaledKineticEnergy<t->theLowestKineticE <<
791 <<
792 dEdx =(*dEdxTable)(materialIndex)->GetVal <<
793 t->theLowestKineticEnergy,isOut) <<
794 *std::sqrt(scaledKineticEnergy/t->t <<
795 <<
796 } else if (scaledKineticEnergy>t->theHighest <<
797 <<
798 dEdx = (*dEdxTable)(materialIndex)->GetVa <<
799 t->theHighestKineticEnergy,isOut); <<
800 <<
801 } else { <<
802 <<
803 dEdx = (*dEdxTable)(materialIndex)->GetVal <<
804 scaledKineticEnergy,isOut); <<
805 <<
806 } <<
807 <<
808 return dEdx*Chargesquare; <<
809 } <<
810 <<
811 //....oooOO0OOooo........oooOO0OOooo........oo <<
812 <<
813 G4double G4EnergyLossTables::GetRange( <<
814 const G4ParticleDefinition *aParticle, <<
815 G4double KineticEnergy, <<
816 const G4MaterialCutsCouple *couple, <<
817 G4bool check) <<
818 { <<
819 if (!t) t = new G4EnergyLossTablesHelper; <<
820 <<
821 if(aParticle != (const G4ParticleDefinition* <<
822 { <<
823 *t= GetTables(aParticle); <<
824 lastParticle = (G4ParticleDefinition*) aPa <<
825 Chargesquare = (aParticle->GetPDGCharge()) <<
826 (aParticle->GetPDGCharge()) <<
827 QQPositron ; <<
828 oldIndex = -1 ; <<
829 } <<
830 const G4PhysicsTable* rangeTable= t->theRang <<
831 const G4PhysicsTable* dEdxTable= t->theDEDX <<
832 if (!rangeTable) { <<
833 if(check) return G4LossTableManager::Insta <<
834 else return DBL_MAX; <<
835 //ParticleHaveNoLoss(aParticle,"Range"); <<
836 } <<
837 <<
838 G4int materialIndex = couple->GetIndex(); <<
839 G4double scaledKineticEnergy = KineticEnergy <<
840 G4double Range; <<
841 G4bool isOut; <<
842 <<
843 if (scaledKineticEnergy<t->theLowestKineticE <<
844 <<
845 Range = std::sqrt(scaledKineticEnergy/t->t <<
846 (*rangeTable)(materialIndex)->GetV <<
847 t->theLowestKineticEnergy,isOut) <<
848 <<
849 } else if (scaledKineticEnergy>t->theHighest <<
850 <<
851 Range = (*rangeTable)(materialIndex)->GetV <<
852 t->theHighestKineticEnergy,isOut)+ <<
853 (scaledKineticEnergy-t->theHighest <<
854 (*dEdxTable)(materialIndex)->GetVa <<
855 t->theHighestKineticEnergy,isOut <<
856 <<
857 } else { <<
858 <<
859 Range = (*rangeTable)(materialIndex)->GetV <<
860 scaledKineticEnergy,isOut); <<
861 <<
862 } <<
863 <<
864 return Range/(Chargesquare*t->theMassRatio); <<
865 } <<
866 <<
867 //....oooOO0OOooo........oooOO0OOooo........oo <<
868 <<
869 G4double G4EnergyLossTables::GetPreciseEnergyF <<
870 const G4P <<
871 G4d <<
872 const G4M <<
873 G4bool check) <<
874 // it returns the value of the kinetic energy <<
875 { <<
876 if (!t) t = new G4EnergyLossTablesHelper; <<
877 <<
878 if( aParticle != (const G4ParticleDefinition <<
879 { <<
880 *t= GetTables(aParticle); <<
881 lastParticle = (G4ParticleDefinition*) aPa <<
882 Chargesquare = (aParticle->GetPDGCharge()) <<
883 (aParticle->GetPDGCharge()) <<
884 QQPositron ; <<
885 oldIndex = -1 ; <<
886 } <<
887 const G4PhysicsTable* dEdxTable= t->theDEDX <<
888 const G4PhysicsTable* inverseRangeTable= t- <<
889 <<
890 if (!inverseRangeTable) { <<
891 if(check) return G4LossTableManager::Insta <<
892 else return DBL_MAX; <<
893 // else ParticleHaveNoLoss(aPartic <<
894 } <<
895 <<
896 G4double scaledrange,scaledKineticEnergy ; <<
897 G4bool isOut ; <<
898 <<
899 G4int materialIndex = couple->GetIndex() ; <<
900 <<
901 if(materialIndex != oldIndex) <<
902 { <<
903 oldIndex = materialIndex ; <<
904 rmin = (*inverseRangeTable)(materialIndex) <<
905 GetLowEdgeEnergy <<
906 rmax = (*inverseRangeTable)(materialIndex) <<
907 GetLowEdgeEnergy(t->theNumb <<
908 Thigh = (*inverseRangeTable)(materialIndex <<
909 GetValue(rmax,is <<
910 } <<
911 <<
912 scaledrange = range*Chargesquare*t->theMassR <<
913 <<
914 if(scaledrange < rmin) <<
915 { <<
916 scaledKineticEnergy = t->theLowestKineticE <<
917 scaledrange*scaledrange/(rm <<
918 } <<
919 else <<
920 { <<
921 if(scaledrange < rmax) <<
922 { <<
923 scaledKineticEnergy = (*inverseRangeTabl <<
924 GetValue( scaled <<
925 } <<
926 else <<
927 { <<
928 scaledKineticEnergy = Thigh + <<
929 (scaledrange-rmax)* <<
930 (*dEdxTable)(materialInd <<
931 GetValue(Thig <<
932 } <<
933 } <<
934 <<
935 return scaledKineticEnergy/t->theMassRatio ; <<
936 } <<
937 <<
938 //....oooOO0OOooo........oooOO0OOooo........oo <<
939 <<
940 G4double G4EnergyLossTables::GetPreciseDEDX( <<
941 const G4ParticleDefinition *aParticle, <<
942 G4double KineticEnergy, <<
943 const G4MaterialCutsCouple *couple) <<
944 { <<
945 if (!t) t = new G4EnergyLossTablesHelper; <<
946 <<
947 if( aParticle != (const G4ParticleDefinition <<
948 { <<
949 *t= GetTables(aParticle); <<
950 lastParticle = (G4ParticleDefinition*) aPa <<
951 Chargesquare = (aParticle->GetPDGCharge()) <<
952 (aParticle->GetPDGCharge()) <<
953 QQPositron ; <<
954 oldIndex = -1 ; <<
955 } <<
956 const G4PhysicsTable* dEdxTable= t->theDEDX <<
957 if ( !dEdxTable ) <<
958 return G4LossTableManager::Instance()->Get <<
959 <<
960 G4int materialIndex = couple->GetIndex(); <<
961 G4double scaledKineticEnergy = KineticEnergy <<
962 G4double dEdx; <<
963 G4bool isOut; <<
964 <<
965 if (scaledKineticEnergy<t->theLowestKineticE <<
966 <<
967 dEdx = std::sqrt(scaledKineticEnergy/t->t <<
968 *(*dEdxTable)(materialIndex)->GetV <<
969 t->theLowestKineticEnergy,isOut) <<
970 <<
971 } else if (scaledKineticEnergy>t->theHighest <<
972 <<
973 dEdx = (*dEdxTable)(materialIndex)->GetVa <<
974 t->theHighestKineticEnergy,isOut); <<
975 <<
976 } else { <<
977 <<
978 dEdx = (*dEdxTable)(materialIndex)->GetV <<
979 scaledKineticEnergy, <<
980 <<
981 } <<
982 <<
983 return dEdx*Chargesquare; <<
984 } <<
985 <<
986 //....oooOO0OOooo........oooOO0OOooo........oo <<
987 <<
988 G4double G4EnergyLossTables::GetPreciseRangeFr <<
989 const G4ParticleDefinition *aParticle, <<
990 G4double KineticEnergy, <<
991 const G4MaterialCutsCouple *couple) <<
992 { <<
993 if (!t) t = new G4EnergyLossTablesHelper; <<
994 <<
995 if( aParticle != (const G4ParticleDefinition <<
996 { <<
997 *t= GetTables(aParticle); <<
998 lastParticle = (G4ParticleDefinition*) aPa <<
999 Chargesquare = (aParticle->GetPDGCharge()) <<
1000 (aParticle->GetPDGCharge() <<
1001 QQPositron ; <<
1002 oldIndex = -1 ; <<
1003 } <<
1004 const G4PhysicsTable* rangeTable= t->theRan <<
1005 const G4PhysicsTable* dEdxTable= t->theDED <<
1006 if ( !dEdxTable || !rangeTable) <<
1007 return G4LossTableManager::Instance()->Ge <<
1008 <<
1009 G4int materialIndex = couple->GetIndex(); <<
1010 <<
1011 G4double Thighr = t->theHighestKineticEnerg <<
1012 (*rangeTable)(materialInde 139 (*rangeTable)(materialIndex)->
1013 GetLowEdgeEnergy(1) ; 140 GetLowEdgeEnergy(1) ;
1014 141
1015 G4double scaledKineticEnergy = KineticEnerg << 142 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
1016 G4double Range; 143 G4double Range;
1017 G4bool isOut; 144 G4bool isOut;
1018 145
1019 if (scaledKineticEnergy<t->theLowestKinetic << 146 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
1020 147
1021 Range = std::sqrt(scaledKineticEnergy/t-> << 148 Range = sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)*
1022 (*rangeTable)(materialIndex)->Get 149 (*rangeTable)(materialIndex)->GetValue(
1023 t->theLowestKineticEnergy,isOut << 150 t.theLowestKineticEnergy,isOut);
1024 151
1025 } else if (scaledKineticEnergy>Thighr) { 152 } else if (scaledKineticEnergy>Thighr) {
1026 153
1027 Range = (*rangeTable)(materialIndex)->Get 154 Range = (*rangeTable)(materialIndex)->GetValue(
1028 Thighr,isOut)+ 155 Thighr,isOut)+
1029 (scaledKineticEnergy-Thighr)/ 156 (scaledKineticEnergy-Thighr)/
1030 (*dEdxTable)(materialIndex)->GetV 157 (*dEdxTable)(materialIndex)->GetValue(
1031 Thighr,isOut); 158 Thighr,isOut);
1032 159
1033 } else { 160 } else {
1034 << 161
1035 Range = (*rangeTable)(materialIndex)->Ge 162 Range = (*rangeTable)(materialIndex)->GetValue(
1036 scaledKineticEnergy,is 163 scaledKineticEnergy,isOut) ;
1037 164
1038 } 165 }
1039 166
1040 return Range/(Chargesquare*t->theMassRatio) << 167 return Range/(Chargesquare*t.theMassRatio);
1041 } 168 }
1042 169
1043 //....oooOO0OOooo........oooOO0OOooo........o 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1044 171
1045 void G4EnergyLossTables::CPRWarning() <<
1046 { <<
1047 if (let_counter < let_max_num_warnings) { <<
1048 172
1049 G4cout << G4endl; <<
1050 G4cout << "##### G4EnergyLossTable WARNIN <<
1051 G4cout << "##### RESULTS ARE NOT GARANTEE <<
1052 G4cout << "##### Please, substitute G4Mat <<
1053 G4cout << "##### Obsolete interface will <<
1054 G4cout << G4endl; <<
1055 let_counter++; <<
1056 173
1057 } else if (let_counter == let_max_num_warni <<
1058 <<
1059 G4cout << "##### G4EnergyLossTable WARNIN <<
1060 let_counter++; <<
1061 } <<
1062 } <<
1063 174
1064 //....oooOO0OOooo........oooOO0OOooo........o <<
1065 <<
1066 void <<
1067 G4EnergyLossTables::ParticleHaveNoLoss(const <<
1068 const G4String& /*q*/) <<
1069 { <<
1070 } <<
1071 <<
1072 //....oooOO0OOooo........oooOO0OOooo........o <<
1073 175