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