| 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 // GEANT 4 class implementation file
30 //
31 // History: based on object model of
32 // 2nd December 1995, G.Cosmo
33 // ---------- G4hEnergyLoss physics proce
34 // by Laszlo Urban, 30 May 1997
35 //
36 // *******************************************
37 // It is the first implementation of the NEW U
38 // It calculates the energy loss of charged ha
39 // *******************************************
40 //
41 // 7/10/98 bug fixes + some cleanup , L.Urb
42 // 22/10/98 cleanup , L.Urban
43 // 07/12/98 works for ions as well+ bug corr
44 // 02/02/99 several bugs fixed, L.Urban
45 // 31/03/00 rename to lowenergy as G4hRDEner
46 // 05/11/00 new method to calculate particle
47 // 10/05/01 V.Ivanchenko Clean up againist L
48 // 23/11/01 V.Ivanchenko Move static member-
49 // 28/10/02 V.Ivanchenko Optimal binning for
50 // 21/01/03 V.Ivanchenko Cut per region
51 // 23/01/03 V.Ivanchenko Fix in table build
52
53 // 31 Jul 2008 MGP Short term supply of en
54 // former G4hLowEnergyLoss
55 // To be replaced by rewor
56 // issues properly
57
58 // 25 Nov 2010 MGP Added some protections
59 // The whole energy loss d
60
61 // 20 Jun 2011 MGP Corrected some compilat
62
63 // -------------------------------------------
64
65 #include "G4hRDEnergyLoss.hh"
66 #include "G4PhysicalConstants.hh"
67 #include "G4SystemOfUnits.hh"
68 #include "G4EnergyLossTables.hh"
69 #include "G4Poisson.hh"
70 #include "G4ProductionCutsTable.hh"
71
72
73 // Initialisation of static members **********
74 // contributing processes : ion.loss ->NumberO
75 // to 1 . YOU DO NOT HAVE TO CHANGE this var
76 // You have to change NumberOfProcesses
77 // if you invent a new process contributing to
78 // NumberOfProcesses should be 2 in this cas
79 // or for debugging purposes.
80 // The NumberOfProcesses data member can be c
81 // functions Get/Set/Plus/MinusNumberOfProces
82
83
84 // The following vectors have a fixed dimensio
85 // filled in with more than 100 elements will
86 G4ThreadLocal G4int G4hRDEnergyLoss::NumberOfP
87
88 G4ThreadLocal G4int G4hRDEnergyLoss
89 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss
90
91 G4ThreadLocal G4int G4hRDEnergyLoss
92 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss
93
94 G4ThreadLocal G4int G4hRDEnergyLoss
95 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss
96
97 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
98 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
99 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
100 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
101 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
102 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
103 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
104 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
105 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
106 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
107
108 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
109 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
110 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
111 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
112 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
113 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
114
115 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
116 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
117 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
118 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
119 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
120
121 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
122 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
123 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss:
124
125 //const G4Proton* G4hRDEnergyLoss::theProton=G
126 //const G4AntiProton* G4hRDEnergyLoss::theAnti
127
128 G4ThreadLocal G4double G4hRDEnergyLoss::Partic
129 G4ThreadLocal G4double G4hRDEnergyLoss::ptable
130 G4ThreadLocal G4double G4hRDEnergyLoss::pbarta
131
132 G4ThreadLocal G4double G4hRDEnergyLoss::Mass,
133 G4hRDEnergyLoss::taulow
134 G4hRDEnergyLoss::tauhig
135 G4hRDEnergyLoss::ltaulo
136 G4hRDEnergyLoss::ltauhi
137
138 G4ThreadLocal G4double G4hRDEnergyLoss::dRover
139 G4ThreadLocal G4double G4hRDEnergyLoss::finalR
140
141 G4ThreadLocal G4double G4hRDEnergyLoss::c1lim
142 G4ThreadLocal G4double G4hRDEnergyLoss::c2lim
143 // 2.*(1.-(0.20))*(200.*micrometer)
144 G4ThreadLocal G4double G4hRDEnergyLoss::c3lim
145 // -(1.-(0.20))*(200.*micrometer)*(200.*micro
146
147 G4ThreadLocal G4double G4hRDEnergyLoss::Charge
148
149 G4ThreadLocal G4bool G4hRDEnergyLoss::rndmSt
150 G4ThreadLocal G4bool G4hRDEnergyLoss::Enloss
151
152 G4ThreadLocal G4double G4hRDEnergyLoss::Lowest
153 G4ThreadLocal G4double G4hRDEnergyLoss::Highes
154 G4ThreadLocal G4int G4hRDEnergyLoss::TotBin
155 G4ThreadLocal G4double G4hRDEnergyLoss::RTable
156
157 //--------------------------------------------
158
159 // constructor and destructor
160
161 G4hRDEnergyLoss::G4hRDEnergyLoss(const G4Strin
162 : G4VContinuousDiscreteProcess (processName)
163 MaxExcitationNumber (1.e6),
164 probLimFluct (0.01),
165 nmaxDirectFluct (100),
166 nmaxCont1(4),
167 nmaxCont2(16),
168 theLossTable(0),
169 linLossLimit(0.05),
170 MinKineticEnergy(0.0)
171 {
172 if (!RecorderOfpbarProcess) RecorderOfpbarPr
173 if (!RecorderOfpProcess) RecorderOfpProcess
174 if (!RecorderOfProcess) RecorderOfProcess =
175 }
176
177 //--------------------------------------------
178
179 G4hRDEnergyLoss::~G4hRDEnergyLoss()
180 {
181 if(theLossTable)
182 {
183 theLossTable->clearAndDestroy();
184 delete theLossTable;
185 }
186 }
187
188 //--------------------------------------------
189
190 G4int G4hRDEnergyLoss::GetNumberOfProcesses()
191 {
192 return NumberOfProcesses;
193 }
194
195 //--------------------------------------------
196
197 void G4hRDEnergyLoss::SetNumberOfProcesses(G4i
198 {
199 NumberOfProcesses=number;
200 }
201
202 //--------------------------------------------
203
204 void G4hRDEnergyLoss::PlusNumberOfProcesses()
205 {
206 NumberOfProcesses++;
207 }
208
209 //--------------------------------------------
210
211 void G4hRDEnergyLoss::MinusNumberOfProcesses()
212 {
213 NumberOfProcesses--;
214 }
215
216 //--------------------------------------------
217
218 void G4hRDEnergyLoss::SetdRoverRange(G4double
219 {
220 dRoverRange = value;
221 }
222
223 //--------------------------------------------
224
225 void G4hRDEnergyLoss::SetRndmStep (G4bool va
226 {
227 rndmStepFlag = value;
228 }
229
230 //--------------------------------------------
231
232 void G4hRDEnergyLoss::SetEnlossFluc (G4bool va
233 {
234 EnlossFlucFlag = value;
235 }
236
237 //--------------------------------------------
238
239 void G4hRDEnergyLoss::SetStepFunction (G4doubl
240 {
241 dRoverRange = c1;
242 finalRange = c2;
243 c1lim=dRoverRange;
244 c2lim=2.*(1-dRoverRange)*finalRange;
245 c3lim=-(1.-dRoverRange)*finalRange*finalRang
246 }
247
248 //--------------------------------------------
249
250 void G4hRDEnergyLoss::BuildDEDXTable(const G4P
251 {
252 // calculate data members TotBin,LOGRTable,
253
254 if (!RecorderOfpbarProcess) RecorderOfpbarPr
255 if (!RecorderOfpProcess) RecorderOfpProcess
256 if (!RecorderOfProcess) RecorderOfProcess =
257
258 const G4ProductionCutsTable* theCoupleTable=
259 G4ProductionCutsTable::GetProductionCutsTa
260 std::size_t numOfCouples = theCoupleTable->G
261
262 // create/fill proton or antiproton tables d
263 Charge = aParticleType.GetPDGCharge()/eplus;
264 ParticleMass = aParticleType.GetPDGMass() ;
265
266 if (Charge>0.) {theDEDXTable= theDEDXpTable;
267 else {theDEDXTable= theDEDXpbarTab
268
269 if( ((Charge>0.) && (theDEDXTable==0)) ||
270 ((Charge<0.) && (theDEDXTable==0))
271 )
272 {
273
274 // Build energy loss table as a sum of t
275 // different processes.
276 if( Charge >0.)
277 {
278 RecorderOfProcess=RecorderOfpProcess;
279 CounterOfProcess=CounterOfpProcess;
280
281 if(CounterOfProcess == NumberOfProcesses)
282 {
283 if(theDEDXpTable)
284 { theDEDXpTable->clearAndDestroy();
285 delete theDEDXpTable; }
286 theDEDXpTable = new G4PhysicsTable(num
287 theDEDXTable = theDEDXpTable;
288 }
289 }
290 else
291 {
292 RecorderOfProcess=RecorderOfpbarProcess;
293 CounterOfProcess=CounterOfpbarProcess;
294
295 if(CounterOfProcess == NumberOfProcesses)
296 {
297 if(theDEDXpbarTable)
298 { theDEDXpbarTable->clearAndDestroy();
299 delete theDEDXpbarTable; }
300 theDEDXpbarTable = new G4PhysicsTable(
301 theDEDXTable = theDEDXpbarTable;
302 }
303 }
304
305 if(CounterOfProcess == NumberOfProcesses
306 {
307 // loop for materials
308 G4double LowEdgeEnergy , Value ;
309 G4bool isOutRange ;
310 G4PhysicsTable* pointer ;
311
312 for (std::size_t J=0; J<numOfCouples; ++J)
313 {
314 // create physics vector and fill it
315 G4PhysicsLogVector* aVector =
316 new G4PhysicsLogVector(LowestK
317 Highest
318
319 // loop for the kinetic energy
320 for (G4int i=0; i<TotBin; i++)
321 {
322 LowEdgeEnergy = aVector->GetLowEdgeEnerg
323 Value = 0. ;
324
325 // loop for the contributing processes
326 for (G4int process=0; process < NumberOf
327 {
328 pointer= RecorderOfProcess[process];
329 Value += (*pointer)[J]->
330 GetValue(LowEdgeEnergy,isOutRange) ;
331 }
332
333 aVector->PutValue(i,Value) ;
334 }
335
336 theDEDXTable->insert(aVector) ;
337 }
338
339 // reset counter to zero ................
340 if( Charge >0.)
341 CounterOfpProcess=0 ;
342 else
343 CounterOfpbarProcess=0 ;
344
345 // Build range table
346 BuildRangeTable( aParticleType);
347
348 // Build lab/proper time tables
349 BuildTimeTables( aParticleType) ;
350
351 // Build coeff tables for the energy loss
352 BuildRangeCoeffATable( aParticleType);
353 BuildRangeCoeffBTable( aParticleType);
354 BuildRangeCoeffCTable( aParticleType);
355
356 // invert the range table
357
358 BuildInverseRangeTable(aParticleType);
359 }
360 }
361 // make the energy loss and the range table
362
363 G4EnergyLossTables::Register(&aParticleType,
364 (Charge>0)?
365 theDEDXpTable: theDEDXpbarTable,
366 (Charge>0)?
367 theRangepTable: theRangepbarTable
368 (Charge>0)?
369 theInverseRangepTable: theInverse
370 (Charge>0)?
371 theLabTimepTable: theLabTimepbarT
372 (Charge>0)?
373 theProperTimepTable: theProperTim
374 LowestKineticEnergy, HighestKinet
375 proton_mass_c2/aParticleType.GetP
376 TotBin);
377
378 }
379
380 //--------------------------------------------
381
382 void G4hRDEnergyLoss::BuildRangeTable(const G4
383 {
384 // Build range table from the energy loss ta
385
386 Mass = aParticleType.GetPDGMass();
387
388 const G4ProductionCutsTable* theCoupleTable=
389 G4ProductionCutsTable::GetProductionCutsTa
390 G4int numOfCouples = (G4int)theCoupleTable->
391
392 if( Charge >0.)
393 {
394 if(theRangepTable)
395 { theRangepTable->clearAndDestroy();
396 delete theRangepTable; }
397 theRangepTable = new G4PhysicsTable(numO
398 theRangeTable = theRangepTable ;
399 }
400 else
401 {
402 if(theRangepbarTable)
403 { theRangepbarTable->clearAndDestroy();
404 delete theRangepbarTable; }
405 theRangepbarTable = new G4PhysicsTable(n
406 theRangeTable = theRangepbarTable ;
407 }
408
409 // loop for materials
410
411 for (G4int J=0; J<numOfCouples; ++J)
412 {
413 G4PhysicsLogVector* aVector;
414 aVector = new G4PhysicsLogVector(LowestK
415 HighestKineticEnergy,TotBin);
416
417 BuildRangeVector(J, aVector);
418 theRangeTable->insert(aVector);
419 }
420 }
421
422 //--------------------------------------------
423
424 void G4hRDEnergyLoss::BuildTimeTables(const G4
425 {
426 const G4ProductionCutsTable* theCoupleTable=
427 G4ProductionCutsTable::GetProductionCutsTa
428 G4int numOfCouples = (G4int)theCoupleTable->
429
430 if(&aParticleType == G4Proton::Proton())
431 {
432 if(theLabTimepTable)
433 { theLabTimepTable->clearAndDestroy();
434 delete theLabTimepTable; }
435 theLabTimepTable = new G4PhysicsTable(nu
436 theLabTimeTable = theLabTimepTable ;
437
438 if(theProperTimepTable)
439 { theProperTimepTable->clearAndDestroy();
440 delete theProperTimepTable; }
441 theProperTimepTable = new G4PhysicsTable
442 theProperTimeTable = theProperTimepTable
443 }
444
445 if(&aParticleType == G4AntiProton::AntiProto
446 {
447 if(theLabTimepbarTable)
448 { theLabTimepbarTable->clearAndDestroy();
449 delete theLabTimepbarTable; }
450 theLabTimepbarTable = new G4PhysicsTable
451 theLabTimeTable = theLabTimepbarTable ;
452
453 if(theProperTimepbarTable)
454 { theProperTimepbarTable->clearAndDestroy();
455 delete theProperTimepbarTable; }
456 theProperTimepbarTable = new G4PhysicsTa
457 theProperTimeTable = theProperTimepbarTa
458 }
459
460 for (G4int J=0; J<numOfCouples; ++J)
461 {
462 G4PhysicsLogVector* aVector;
463 G4PhysicsLogVector* bVector;
464
465 aVector = new G4PhysicsLogVector(LowestK
466 HighestKineticEnergy,TotBin);
467
468 BuildLabTimeVector(J, aVector);
469 theLabTimeTable->insert(aVector);
470
471 bVector = new G4PhysicsLogVector(LowestK
472 HighestKineticEnergy,TotBin);
473
474 BuildProperTimeVector(J, bVector);
475 theProperTimeTable->insert(bVector);
476 }
477 }
478
479 //--------------------------------------------
480
481 void G4hRDEnergyLoss::BuildRangeVector(G4int m
482 G4PhysicsLogVector* rangeVector
483 {
484 // create range vector for a material
485
486 G4bool isOut;
487 G4PhysicsVector* physicsVector= (*theDEDXTab
488 G4double energy1 = rangeVector->GetLowEdgeEn
489 G4double dedx = physicsVector->GetValue(e
490 G4double range = 0.5*energy1/dedx;
491 rangeVector->PutValue(0,range);
492 G4int n = 100;
493 G4double del = 1.0/(G4double)n ;
494
495 for (G4int j=1; j<TotBin; j++) {
496
497 G4double energy2 = rangeVector->GetLowEdge
498 G4double de = (energy2 - energy1) * del ;
499 G4double dedx1 = dedx ;
500
501 for (G4int i=1; i<n; i++) {
502 G4double energy = energy1 + i*de ;
503 G4double dedx2 = physicsVector->GetValu
504 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2)
505 dedx1 = dedx2;
506 }
507 rangeVector->PutValue(j,range);
508 dedx = dedx1 ;
509 energy1 = energy2 ;
510 }
511 }
512
513 //--------------------------------------------
514
515 void G4hRDEnergyLoss::BuildLabTimeVector(G4int
516 G4PhysicsLogVector* timeVector)
517 {
518 // create lab time vector for a material
519
520 G4int nbin=100;
521 G4bool isOut;
522 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p
523 //G4double losslim,clim,taulim,timelim,ltaul
524 G4double losslim,clim,taulim,timelim,
525 LowEdgeEnergy,tau,Value ;
526
527 G4PhysicsVector* physicsVector= (*theDEDXTab
528
529 // low energy part first...
530 losslim = physicsVector->GetValue(tlim,isOut
531 taulim=tlim/ParticleMass ;
532 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh
533 //ltaulim = std::log(taulim);
534 //ltaumax = std::log(HighestKineticEnergy/Pa
535
536 G4int i=-1;
537 G4double oldValue = 0. ;
538 G4double tauold ;
539 do
540 {
541 i += 1 ;
542 LowEdgeEnergy = timeVector->GetLowEdgeEn
543 tau = LowEdgeEnergy/ParticleMass ;
544 if ( tau <= taulim )
545 {
546 Value = clim*std::exp(ppar*std::log(tau/ta
547 }
548 else
549 {
550 timelim=clim ;
551 ltaulow = std::log(taulim);
552 ltauhigh = std::log(tau);
553 Value = timelim+LabTimeIntLog(physicsVecto
554 }
555 timeVector->PutValue(i,Value);
556 oldValue = Value ;
557 tauold = tau ;
558 } while (tau<=taulim) ;
559
560 i += 1 ;
561 for (G4int j=i; j<TotBin; j++)
562 {
563 LowEdgeEnergy = timeVector->GetLowEdgeEn
564 tau = LowEdgeEnergy/ParticleMass ;
565 ltaulow = std::log(tauold);
566 ltauhigh = std::log(tau);
567 Value = oldValue+LabTimeIntLog(physicsVe
568 timeVector->PutValue(j,Value);
569 oldValue = Value ;
570 tauold = tau ;
571 }
572 }
573
574 //--------------------------------------------
575
576 void G4hRDEnergyLoss::BuildProperTimeVector(G4
577 G4PhysicsLogVector* timeVector)
578 {
579 // create proper time vector for a material
580
581 G4int nbin=100;
582 G4bool isOut;
583 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p
584 //G4double losslim,clim,taulim,timelim,ltaul
585 G4double losslim,clim,taulim,timelim,
586 LowEdgeEnergy,tau,Value ;
587
588 G4PhysicsVector* physicsVector= (*theDEDXTab
589
590 // low energy part first...
591 losslim = physicsVector->GetValue(tlim,isOut
592 taulim=tlim/ParticleMass ;
593 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh
594 //ltaulim = std::log(taulim);
595 //ltaumax = std::log(HighestKineticEnergy/Pa
596
597 G4int i=-1;
598 G4double oldValue = 0. ;
599 G4double tauold ;
600 do
601 {
602 i += 1 ;
603 LowEdgeEnergy = timeVector->GetLowEdgeEn
604 tau = LowEdgeEnergy/ParticleMass ;
605 if ( tau <= taulim )
606 {
607 Value = clim*std::exp(ppar*std::log(tau/ta
608 }
609 else
610 {
611 timelim=clim ;
612 ltaulow = std::log(taulim);
613 ltauhigh = std::log(tau);
614 Value = timelim+ProperTimeIntLog(physicsVe
615 }
616 timeVector->PutValue(i,Value);
617 oldValue = Value ;
618 tauold = tau ;
619 } while (tau<=taulim) ;
620
621 i += 1 ;
622 for (G4int j=i; j<TotBin; j++)
623 {
624 LowEdgeEnergy = timeVector->GetLowEdgeEn
625 tau = LowEdgeEnergy/ParticleMass ;
626 ltaulow = std::log(tauold);
627 ltauhigh = std::log(tau);
628 Value = oldValue+ProperTimeIntLog(physic
629 timeVector->PutValue(j,Value);
630 oldValue = Value ;
631 tauold = tau ;
632 }
633 }
634
635 //--------------------------------------------
636
637 G4double G4hRDEnergyLoss::RangeIntLin(G4Physic
638 G4int nbin)
639 {
640 // num. integration, linear binning
641
642 G4double dtau,Value,taui,ti,lossi,ci;
643 G4bool isOut;
644 dtau = (tauhigh-taulow)/nbin;
645 Value = 0.;
646
647 for (G4int i=0; i<=nbin; i++)
648 {
649 taui = taulow + dtau*i ;
650 ti = Mass*taui;
651 lossi = physicsVector->GetValue(ti,isOut
652 if(i==0)
653 ci=0.5;
654 else
655 {
656 if(i<nbin)
657 ci=1.;
658 else
659 ci=0.5;
660 }
661 Value += ci/lossi;
662 }
663 Value *= Mass*dtau;
664 return Value;
665 }
666
667 //--------------------------------------------
668
669 G4double G4hRDEnergyLoss::RangeIntLog(G4Physic
670 G4int nbin)
671 {
672 // num. integration, logarithmic binning
673
674 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci
675 G4bool isOut;
676 ltt = ltauhigh-ltaulow;
677 dltau = ltt/nbin;
678 Value = 0.;
679
680 for (G4int i=0; i<=nbin; i++)
681 {
682 ui = ltaulow+dltau*i;
683 taui = std::exp(ui);
684 ti = Mass*taui;
685 lossi = physicsVector->GetValue(ti,isOut
686 if(i==0)
687 ci=0.5;
688 else
689 {
690 if(i<nbin)
691 ci=1.;
692 else
693 ci=0.5;
694 }
695 Value += ci*taui/lossi;
696 }
697 Value *= Mass*dltau;
698 return Value;
699 }
700
701 //--------------------------------------------
702
703 G4double G4hRDEnergyLoss::LabTimeIntLog(G4Phys
704 G4int nbin)
705 {
706 // num. integration, logarithmic binning
707
708 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci
709 G4bool isOut;
710 ltt = ltauhigh-ltaulow;
711 dltau = ltt/nbin;
712 Value = 0.;
713
714 for (G4int i=0; i<=nbin; i++)
715 {
716 ui = ltaulow+dltau*i;
717 taui = std::exp(ui);
718 ti = ParticleMass*taui;
719 lossi = physicsVector->GetValue(ti,isOut
720 if(i==0)
721 ci=0.5;
722 else
723 {
724 if(i<nbin)
725 ci=1.;
726 else
727 ci=0.5;
728 }
729 Value += ci*taui*(ti+ParticleMass)/(std:
730 }
731 Value *= ParticleMass*dltau/c_light;
732 return Value;
733 }
734
735 //--------------------------------------------
736
737 G4double G4hRDEnergyLoss::ProperTimeIntLog(G4P
738 G4int nbin)
739 {
740 // num. integration, logarithmic binning
741
742 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci
743 G4bool isOut;
744 ltt = ltauhigh-ltaulow;
745 dltau = ltt/nbin;
746 Value = 0.;
747
748 for (G4int i=0; i<=nbin; i++)
749 {
750 ui = ltaulow+dltau*i;
751 taui = std::exp(ui);
752 ti = ParticleMass*taui;
753 lossi = physicsVector->GetValue(ti,isOut
754 if(i==0)
755 ci=0.5;
756 else
757 {
758 if(i<nbin)
759 ci=1.;
760 else
761 ci=0.5;
762 }
763 Value += ci*taui*ParticleMass/(std::sqrt
764 }
765 Value *= ParticleMass*dltau/c_light;
766 return Value;
767 }
768
769 //--------------------------------------------
770
771 void G4hRDEnergyLoss::BuildRangeCoeffATable( c
772 {
773 // Build tables of coefficients for the ener
774 // create table for coefficients "A"
775
776 G4int numOfCouples = (G4int)G4ProductionCuts
777
778 if(Charge>0.)
779 {
780 if(thepRangeCoeffATable)
781 { thepRangeCoeffATable->clearAndDestroy();
782 delete thepRangeCoeffATable; }
783 thepRangeCoeffATable = new G4PhysicsTabl
784 theRangeCoeffATable = thepRangeCoeffATab
785 theRangeTable = theRangepTable ;
786 }
787 else
788 {
789 if(thepbarRangeCoeffATable)
790 { thepbarRangeCoeffATable->clearAndDestroy()
791 delete thepbarRangeCoeffATable; }
792 thepbarRangeCoeffATable = new G4PhysicsT
793 theRangeCoeffATable = thepbarRangeCoeffA
794 theRangeTable = theRangepbarTable ;
795 }
796
797 G4double R2 = RTable*RTable ;
798 G4double R1 = RTable+1.;
799 G4double w = R1*(RTable-1.)*(RTable-1.);
800 G4double w1 = RTable/w , w2 = -RTable*R1/w ,
801 G4double Ti , Tim , Tip , Ri , Rim , Rip , V
802 G4bool isOut;
803
804 // loop for materials
805 for (G4int J=0; J<numOfCouples; J++)
806 {
807 G4int binmax=TotBin ;
808 G4PhysicsLinearVector* aVector =
809 new G4PhysicsLinearVector(0.,binmax, TotBin)
810 Ti = LowestKineticEnergy ;
811 if (Ti < DBL_MIN) Ti = 1.e-8;
812 G4PhysicsVector* rangeVector= (*theRange
813
814 for ( G4int i=0; i<TotBin; i++)
815 {
816 Ri = rangeVector->GetValue(Ti,isOut) ;
817 if (Ti < DBL_MIN) Ti = 1.e-8;
818 if ( i==0 )
819 Rim = 0. ;
820 else
821 {
822 // ---- MGP ---- Modified to avoid a f
823 // The correction is just a temporary
824 // Original: Tim = Ti/RTable results
825 if (RTable != 0.)
826 {
827 Tim = Ti/RTable ;
828 }
829 else
830 {
831 Tim = 0.;
832 }
833 Rim = rangeVector->GetValue(Tim,isOut)
834 }
835 if ( i==(TotBin-1))
836 Rip = Ri ;
837 else
838 {
839 Tip = Ti*RTable ;
840 Rip = rangeVector->GetValue(Tip,isOut)
841 }
842 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti)
843
844 aVector->PutValue(i,Value);
845 Ti = RTable*Ti ;
846 }
847
848 theRangeCoeffATable->insert(aVector);
849 }
850 }
851
852 //--------------------------------------------
853
854 void G4hRDEnergyLoss::BuildRangeCoeffBTable( c
855 {
856 // Build tables of coefficients for the ener
857 // create table for coefficients "B"
858
859 G4int numOfCouples =
860 (G4int)G4ProductionCutsTable::GetProdu
861
862 if(Charge>0.)
863 {
864 if(thepRangeCoeffBTable)
865 { thepRangeCoeffBTable->clearAndDestroy();
866 delete thepRangeCoeffBTable; }
867 thepRangeCoeffBTable = new G4PhysicsTabl
868 theRangeCoeffBTable = thepRangeCoeffBTab
869 theRangeTable = theRangepTable ;
870 }
871 else
872 {
873 if(thepbarRangeCoeffBTable)
874 { thepbarRangeCoeffBTable->clearAndDestroy()
875 delete thepbarRangeCoeffBTable; }
876 thepbarRangeCoeffBTable = new G4PhysicsT
877 theRangeCoeffBTable = thepbarRangeCoeffB
878 theRangeTable = theRangepbarTable ;
879 }
880
881 G4double R2 = RTable*RTable ;
882 G4double R1 = RTable+1.;
883 G4double w = R1*(RTable-1.)*(RTable-1.);
884 if (w < DBL_MIN) w = DBL_MIN;
885 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3
886 G4double Ti , Tim , Tip , Ri , Rim , Rip , V
887 G4bool isOut;
888
889 // loop for materials
890 for (G4int J=0; J<numOfCouples; J++)
891 {
892 G4int binmax=TotBin ;
893 G4PhysicsLinearVector* aVector =
894 new G4PhysicsLinearVector(0.,binmax, TotBin)
895 Ti = LowestKineticEnergy ;
896 if (Ti < DBL_MIN) Ti = 1.e-8;
897 G4PhysicsVector* rangeVector= (*theRange
898
899 for ( G4int i=0; i<TotBin; i++)
900 {
901 Ri = rangeVector->GetValue(Ti,isOut) ;
902 if (Ti < DBL_MIN) Ti = 1.e-8;
903 if ( i==0 )
904 Rim = 0. ;
905 else
906 {
907 if (RTable < DBL_MIN) RTable = DBL_MIN
908 Tim = Ti/RTable ;
909 Rim = rangeVector->GetValue(Tim,isOut)
910 }
911 if ( i==(TotBin-1))
912 Rip = Ri ;
913 else
914 {
915 Tip = Ti*RTable ;
916 Rip = rangeVector->GetValue(Tip,isOut)
917 }
918 if (Ti < DBL_MIN) Ti = DBL_MIN;
919 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
920
921 aVector->PutValue(i,Value);
922 Ti = RTable*Ti ;
923 }
924 theRangeCoeffBTable->insert(aVector);
925 }
926 }
927
928 //--------------------------------------------
929
930 void G4hRDEnergyLoss::BuildRangeCoeffCTable( c
931 {
932 // Build tables of coefficients for the ener
933 // create table for coefficients "C"
934
935 G4int numOfCouples =
936 (G4int)G4ProductionCutsTable::GetProdu
937
938 if(Charge>0.)
939 {
940 if(thepRangeCoeffCTable)
941 { thepRangeCoeffCTable->clearAndDestroy();
942 delete thepRangeCoeffCTable; }
943 thepRangeCoeffCTable = new G4PhysicsTabl
944 theRangeCoeffCTable = thepRangeCoeffCTab
945 theRangeTable = theRangepTable ;
946 }
947 else
948 {
949 if(thepbarRangeCoeffCTable)
950 { thepbarRangeCoeffCTable->clearAndDestroy()
951 delete thepbarRangeCoeffCTable; }
952 thepbarRangeCoeffCTable = new G4PhysicsT
953 theRangeCoeffCTable = thepbarRangeCoeffC
954 theRangeTable = theRangepbarTable ;
955 }
956
957 G4double R2 = RTable*RTable ;
958 G4double R1 = RTable+1.;
959 G4double w = R1*(RTable-1.)*(RTable-1.);
960 G4double w1 = 1./w , w2 = -RTable*R1/w , w3
961 G4double Ti , Tim , Tip , Ri , Rim , Rip , V
962 G4bool isOut;
963
964 // loop for materials
965 for (G4int J=0; J<numOfCouples; J++)
966 {
967 G4int binmax=TotBin ;
968 G4PhysicsLinearVector* aVector =
969 new G4PhysicsLinearVector(0.,binmax, TotBin)
970 Ti = LowestKineticEnergy ;
971 G4PhysicsVector* rangeVector= (*theRange
972
973 for ( G4int i=0; i<TotBin; i++)
974 {
975 Ri = rangeVector->GetValue(Ti,isOut) ;
976 if ( i==0 )
977 Rim = 0. ;
978 else
979 {
980 Tim = Ti/RTable ;
981 Rim = rangeVector->GetValue(Tim,isOut)
982 }
983 if ( i==(TotBin-1))
984 Rip = Ri ;
985 else
986 {
987 Tip = Ti*RTable ;
988 Rip = rangeVector->GetValue(Tip,isOut)
989 }
990 Value = w1*Rip + w2*Ri + w3*Rim ;
991
992 aVector->PutValue(i,Value);
993 Ti = RTable*Ti ;
994 }
995 theRangeCoeffCTable->insert(aVector);
996 }
997 }
998
999 //--------------------------------------------
1000
1001 void G4hRDEnergyLoss::
1002 BuildInverseRangeTable(const G4ParticleDefini
1003 {
1004 // Build inverse table of the range table
1005
1006 G4bool b;
1007
1008 const G4ProductionCutsTable* theCoupleTable
1009 G4ProductionCutsTable::GetProductionCutsT
1010 std::size_t numOfCouples = theCoupleTable->
1011
1012 if(&aParticleType == G4Proton::Proton())
1013 {
1014 if(theInverseRangepTable)
1015 { theInverseRangepTable->clearAndDestroy();
1016 delete theInverseRangepTable; }
1017 theInverseRangepTable = new G4PhysicsTa
1018 theInverseRangeTable = theInverseRangep
1019 theRangeTable = theRangepTable ;
1020 theDEDXTable = theDEDXpTable ;
1021 theRangeCoeffATable = thepRangeCoeffATa
1022 theRangeCoeffBTable = thepRangeCoeffBTa
1023 theRangeCoeffCTable = thepRangeCoeffCTa
1024 }
1025
1026 if(&aParticleType == G4AntiProton::AntiProt
1027 {
1028 if(theInverseRangepbarTable)
1029 { theInverseRangepbarTable->clearAndDestroy
1030 delete theInverseRangepbarTable; }
1031 theInverseRangepbarTable = new G4Physic
1032 theInverseRangeTable = theInverseRangep
1033 theRangeTable = theRangepbarTable ;
1034 theDEDXTable = theDEDXpbarTable ;
1035 theRangeCoeffATable = thepbarRangeCoeff
1036 theRangeCoeffBTable = thepbarRangeCoeff
1037 theRangeCoeffCTable = thepbarRangeCoeff
1038 }
1039
1040 // loop for materials
1041 for (std::size_t i=0; i<numOfCouples; ++i)
1042 {
1043
1044 G4PhysicsVector* pv = (*theRangeTable)[
1045 std::size_t nbins = pv->GetVectorLeng
1046 G4double elow = pv->GetLowEdgeEnergy(0
1047 G4double ehigh = pv->GetLowEdgeEnergy(n
1048 G4double rlow = pv->GetValue(elow, b);
1049 G4double rhigh = pv->GetValue(ehigh, b)
1050
1051 if (rlow <DBL_MIN) rlow = 1.e-8;
1052 if (rhigh > 1.e16) rhigh = 1.e16;
1053 if (rhigh < 1.e-8) rhigh =1.e-8;
1054 G4double tmpTrick = rhigh/rlow;
1055
1056 //std::cout << nbins << ", elow " << el
1057 // << ", rlow " << rlow << ", rhigh
1058 // << ", trick " << tmpT
1059
1060 if (tmpTrick <= 0. || tmpTrick < DBL_MI
1061 if (tmpTrick > 1.e16) tmpTrick = 1.e16;
1062
1063 rhigh *= std::exp(std::log(tmpTrick)/((
1064
1065 G4PhysicsLogVector* v = new G4PhysicsLo
1066
1067 v->PutValue(0,elow);
1068 G4double energy1 = elow;
1069 G4double range1 = rlow;
1070 G4double energy2 = elow;
1071 G4double range2 = rlow;
1072 std::size_t ilow = 0;
1073 std::size_t ihigh;
1074
1075 for (std::size_t j=1; j<nbins; ++j) {
1076
1077 G4double range = v->GetLowEdgeEnergy(j);
1078
1079 for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
1080 energy2 = pv->GetLowEdgeEnergy(ihigh);
1081 range2 = pv->GetValue(energy2, b);
1082 if(range2 >= range || ihigh == nbins-1) {
1083 ilow = ihigh - 1;
1084 energy1 = pv->GetLowEdgeEnergy(ilow);
1085 range1 = pv->GetValue(energy1, b);
1086 break;
1087 }
1088 }
1089
1090 G4double e = std::log(energy1) + std::log(e
1091
1092 v->PutValue(j,std::exp(e));
1093 }
1094 theInverseRangeTable->insert(v);
1095
1096 }
1097 }
1098
1099 //-------------------------------------------
1100
1101 void G4hRDEnergyLoss::InvertRangeVector(G4int
1102 G4PhysicsLogVector* aVector)
1103 {
1104 // invert range vector for a material
1105
1106 G4double LowEdgeRange,A,B,C,discr,KineticEn
1107 G4double Tbin = LowestKineticEnergy/RTable
1108 G4double rangebin = 0.0 ;
1109 G4int binnumber = -1 ;
1110 G4bool isOut ;
1111
1112
1113 //loop for range values
1114 for( G4int i=0; i<TotBin; i++)
1115 {
1116 LowEdgeRange = aVector->GetLowEdgeEnerg
1117
1118 if( rangebin < LowEdgeRange )
1119 {
1120 do
1121 {
1122 binnumber += 1 ;
1123 Tbin *= RTable ;
1124 rangebin = (*theRangeTable)(materialI
1125 }
1126 while ((rangebin < LowEdgeRange) && (binn
1127 }
1128
1129 if(binnumber == 0)
1130 KineticEnergy = LowestKineticEnergy ;
1131 else if(binnumber == TotBin-1)
1132 KineticEnergy = HighestKineticEnergy ;
1133 else
1134 {
1135 A = (*(*theRangeCoeffATable)(materialInde
1136 B = (*(*theRangeCoeffBTable)(materialInde
1137 C = (*(*theRangeCoeffCTable)(materialInde
1138 if(A==0.)
1139 KineticEnergy = (LowEdgeRange -C )/B ;
1140 else
1141 {
1142 discr = B*B - 4.*A*(C-LowEdgeRange);
1143 discr = discr>0. ? std::sqrt(discr) :
1144 KineticEnergy = 0.5*(discr-B)/A ;
1145 }
1146 }
1147
1148 aVector->PutValue(i,KineticEnergy) ;
1149 }
1150 }
1151
1152 //-------------------------------------------
1153
1154 G4bool G4hRDEnergyLoss::CutsWhereModified()
1155 {
1156 G4bool wasModified = false;
1157 const G4ProductionCutsTable* theCoupleTable
1158 G4ProductionCutsTable::GetProductionCutsT
1159 G4int numOfCouples = (G4int)theCoupleTable-
1160
1161 for (G4int j=0; j<numOfCouples; ++j){
1162 if (theCoupleTable->GetMaterialCutsCouple
1163 wasModified = true;
1164 break;
1165 }
1166 }
1167 return wasModified;
1168 }
1169
1170 //-------------------------------------------
1171 //--- --- --- --- --- --- --- --- --- --- ---
1172
1173 //G4bool G4hRDEnergyLoss::IsApplicable(const
1174 //
1175 //{
1176 // return((particle.GetPDGCharge()!= 0.) &&
1177 //}
1178
1179 //G4double G4hRDEnergyLoss::GetContinuousStep
1180 //
1181 //
1182 //
1183 //
1184 //{
1185 //
1186 // G4double Step =
1187 // GetConstraints(track.GetDynamicParticle
1188 //
1189 // if((Step>0.0)&&(Step<currentMinimumStep))
1190 // currentMinimumStep = Step ;
1191 //
1192 // return Step ;
1193 //}
1194