| 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 //
27 //
28 // -------------------------------------------
29 // first version created by P.Urban , 06/04/19
30 // modifications + "precise" functions added b
31 // modifications , TOF functions , 26/10/98, L
32 // cache mechanism in order to gain time, 11/0
33 // bug fixed , 12/04/99 , L.Urban
34 // 10.11.99: moved from RWT hash dictionary to
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 // -------------------------------------------
42
43 #include "G4EnergyLossTables.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4MaterialCutsCouple.hh"
46 #include "G4RegionStore.hh"
47 #include "G4LossTableManager.hh"
48
49
50 //....oooOO0OOooo........oooOO0OOooo........oo
51
52 G4EnergyLossTablesHelper *G4EnergyLossTables::
53 G4EnergyLossTablesHelper *G4EnergyLossTables::
54 G4ParticleDefinition* G4EnergyLossTables::last
55 G4double G4EnergyLossTables::QQPositron = 1.0;
56 G4double G4EnergyLossTables::Chargesquare ;
57 G4int G4EnergyLossTables::oldIndex = -1 ;
58 G4double G4EnergyLossTables::rmin = 0. ;
59 G4double G4EnergyLossTables::rmax = 0. ;
60 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
90 //....oooOO0OOooo........oooOO0OOooo........oo
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
102 //....oooOO0OOooo........oooOO0OOooo........oo
103
104 void G4EnergyLossTables::Register(
105 const G4ParticleDefinition* p,
106 const G4PhysicsTable* tDEDX,
107 const G4PhysicsTable* tRange,
108 const G4PhysicsTable* tInverseRange,
109 const G4PhysicsTable* tLabTime,
110 const G4PhysicsTable* tProperTime,
111 G4double lowestKineticEnergy,
112 G4double highestKineticEnergy,
113 G4double massRatio,
114 G4int NumberOfBins)
115 {
116 if (!dict) dict = new G4EnergyLossTables::he
117 if (!null_loss) null_loss = new G4EnergyLoss
118 if (!t) t = new G4EnergyLossTablesHelper;
119
120 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX,
121 tLabTime,tProperTime,lowes
122 highestKineticEnergy, massRatio,Number
123
124 *t = GetTables(p) ; // important for cach
125 lastParticle = (G4ParticleDefinition*) p ;
126 Chargesquare = (p->GetPDGCharge())*(p->GetPD
127 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 }
647
648 //....oooOO0OOooo........oooOO0OOooo........oo
649
650 G4double G4EnergyLossTables::GetPreciseDEDX(
651 const G4ParticleDefinition *aParticle,
652 G4double KineticEnergy,
653 const G4Material *aMaterial)
654 {
655 if (!t) t = new G4EnergyLossTablesHelper;
656
657 CPRWarning();
658 if( aParticle != (const G4ParticleDefinition
659 {
660 *t= GetTables(aParticle);
661 lastParticle = (G4ParticleDefinition*) aPa
662 Chargesquare = (aParticle->GetPDGCharge())
663 (aParticle->GetPDGCharge())
664 QQPositron ;
665 oldIndex = -1 ;
666 }
667 const G4PhysicsTable* dEdxTable= t->theDEDX
668 if (!dEdxTable) {
669 ParticleHaveNoLoss(aParticle,"dEdx");
670 return 0.0;
671 }
672
673 G4int materialIndex = (G4int)aMaterial->GetI
674 G4double scaledKineticEnergy = KineticEnergy
675 G4double dEdx;
676 G4bool isOut;
677
678 if (scaledKineticEnergy<t->theLowestKineticE
679
680 dEdx = std::sqrt(scaledKineticEnergy/t->t
681 *(*dEdxTable)(materialIndex)->GetV
682 t->theLowestKineticEnergy,isOut)
683
684 } else if (scaledKineticEnergy>t->theHighest
685
686 dEdx = (*dEdxTable)(materialIndex)->GetVa
687 t->theHighestKineticEnergy,isOut);
688
689 } else {
690
691 dEdx = (*dEdxTable)(materialIndex)->GetV
692 scaledKineticEnergy,
693
694 }
695
696 return dEdx*Chargesquare;
697 }
698
699 //....oooOO0OOooo........oooOO0OOooo........oo
700
701 G4double G4EnergyLossTables::GetPreciseRangeF
702 const G4ParticleDefinition *aParticle,
703 G4double KineticEnergy,
704 const G4Material *aMaterial)
705 {
706 if (!t) t = new G4EnergyLossTablesHelper;
707
708 CPRWarning();
709 if( aParticle != (const G4ParticleDefinition
710 {
711 *t= GetTables(aParticle);
712 lastParticle = (G4ParticleDefinition*) aPa
713 Chargesquare = (aParticle->GetPDGCharge())
714 (aParticle->GetPDGCharge())
715 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 }
724 G4int materialIndex = (G4int)aMaterial->GetI
725
726 G4double Thighr = t->theHighestKineticEnergy
727 (*rangeTable)(materialIndex
728 GetLowEdgeEnergy(1) ;
729
730 G4double scaledKineticEnergy = KineticEnergy
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
1013 GetLowEdgeEnergy(1) ;
1014
1015 G4double scaledKineticEnergy = KineticEnerg
1016 G4double Range;
1017 G4bool isOut;
1018
1019 if (scaledKineticEnergy<t->theLowestKinetic
1020
1021 Range = std::sqrt(scaledKineticEnergy/t->
1022 (*rangeTable)(materialIndex)->Get
1023 t->theLowestKineticEnergy,isOut
1024
1025 } else if (scaledKineticEnergy>Thighr) {
1026
1027 Range = (*rangeTable)(materialIndex)->Get
1028 Thighr,isOut)+
1029 (scaledKineticEnergy-Thighr)/
1030 (*dEdxTable)(materialIndex)->GetV
1031 Thighr,isOut);
1032
1033 } else {
1034
1035 Range = (*rangeTable)(materialIndex)->Ge
1036 scaledKineticEnergy,is
1037
1038 }
1039
1040 return Range/(Chargesquare*t->theMassRatio)
1041 }
1042
1043 //....oooOO0OOooo........oooOO0OOooo........o
1044
1045 void G4EnergyLossTables::CPRWarning()
1046 {
1047 if (let_counter < let_max_num_warnings) {
1048
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
1057 } else if (let_counter == let_max_num_warni
1058
1059 G4cout << "##### G4EnergyLossTable WARNIN
1060 let_counter++;
1061 }
1062 }
1063
1064 //....oooOO0OOooo........oooOO0OOooo........o
1065
1066 void
1067 G4EnergyLossTables::ParticleHaveNoLoss(const
1068 const G4String& /*q*/)
1069 {
1070 }
1071
1072 //....oooOO0OOooo........oooOO0OOooo........o
1073