| 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 // -------------------------------------------
30 // GEANT 4 class implementation file
31 //
32 // History: first implementation, based on ob
33 // 2nd December 1995, G.Cosmo
34 // -------------------------------------------
35 //
36 // Modifications:
37 // 20/09/00 update fluctuations V.Ivanchenko
38 // 22/11/00 minor fix in fluctuations V.Ivanch
39 // 10/05/01 V.Ivanchenko Clean up againist Li
40 // 22/05/01 V.Ivanchenko Update range calcula
41 // 23/11/01 V.Ivanchenko Move static member-f
42 // 22/01/03 V.Ivanchenko Cut per region
43 // 11/02/03 V.Ivanchenko Add limits to fluctu
44 // 24/04/03 V.Ivanchenko Fix the problem of t
45 //
46 // -------------------------------------------
47
48 #include "G4RDVeLowEnergyLoss.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4ProductionCutsTable.hh"
52
53 G4double G4RDVeLowEnergyLoss::ParticleMass
54 G4double G4RDVeLowEnergyLoss::taulow
55 G4double G4RDVeLowEnergyLoss::tauhigh
56 G4double G4RDVeLowEnergyLoss::ltaulow
57 G4double G4RDVeLowEnergyLoss::ltauhigh
58
59
60 G4bool G4RDVeLowEnergyLoss::rndmStepFlag
61 G4bool G4RDVeLowEnergyLoss::EnlossFlucFl
62 G4double G4RDVeLowEnergyLoss::dRoverRange
63 G4double G4RDVeLowEnergyLoss::finalRange
64 G4double G4RDVeLowEnergyLoss::c1lim = dRov
65 G4double G4RDVeLowEnergyLoss::c2lim = 2.*(
66 G4double G4RDVeLowEnergyLoss::c3lim = -(1.
67
68
69 //
70
71 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()
72 :G4VContinuousDiscreteProce
73 lastMaterial(0),
74 nmaxCont1(4),
75 nmaxCont2(16)
76 {
77 G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowE
78 FatalException, "Default constru
79 }
80
81 //
82
83 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const
84 : G4VContinuousDiscreteProce
85 lastMaterial(0),
86 nmaxCont1(4),
87 nmaxCont2(16)
88 {
89 }
90
91 //
92
93 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss()
94 {
95 }
96
97 //
98
99 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDV
100 : G4VContinuousDiscreteProce
101 lastMaterial(0),
102 nmaxCont1(4),
103 nmaxCont2(16)
104 {
105 }
106
107 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool v
108 {
109 rndmStepFlag = value;
110 }
111
112 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool
113 {
114 EnlossFlucFlag = value;
115 }
116
117 void G4RDVeLowEnergyLoss::SetStepFunction (G4d
118 {
119 dRoverRange = c1;
120 finalRange = c2;
121 c1lim=dRoverRange;
122 c2lim=2.*(1-dRoverRange)*finalRange;
123 c3lim=-(1.-dRoverRange)*finalRange*finalRan
124 }
125
126 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang
127 G4PhysicsTable* theRangeTable,
128 G4double lowestKineticEnergy,
129 G4double highestKineticEnergy,
130 G4int TotBin)
131 // Build range table from the energy loss tabl
132 {
133
134 G4int numOfCouples = theDEDXTable->length()
135
136 if(theRangeTable)
137 { theRangeTable->clearAndDestroy();
138 delete theRangeTable; }
139 theRangeTable = new G4PhysicsTable(numOfCou
140
141 // loop for materials
142
143 for (G4int J=0; J<numOfCouples; J++)
144 {
145 G4PhysicsLogVector* aVector;
146 aVector = new G4PhysicsLogVector(lowestKi
147 highestKineticEn
148 BuildRangeVector(theDEDXTable,lowestKinet
149 TotBin,J,aVector);
150 theRangeTable->insert(aVector);
151 }
152 return theRangeTable ;
153 }
154
155 //
156
157 void G4RDVeLowEnergyLoss::BuildRangeVector(G4P
158 G4dou
159 G4dou
160 G4int
161 G4int
162 G4Phy
163
164 // create range vector for a material
165 {
166 G4bool isOut;
167 G4PhysicsVector* physicsVector= (*theDEDXTab
168 G4double energy1 = lowestKineticEnergy;
169 G4double dedx = physicsVector->GetValue(e
170 G4double range = 0.5*energy1/dedx;
171 rangeVector->PutValue(0,range);
172 G4int n = 100;
173 G4double del = 1.0/(G4double)n ;
174
175 for (G4int j=1; j<TotBin; j++) {
176
177 G4double energy2 = rangeVector->GetLowEdge
178 G4double de = (energy2 - energy1) * del ;
179 G4double dedx1 = dedx ;
180
181 for (G4int i=1; i<n; i++) {
182 G4double energy = energy1 + i*de ;
183 G4double dedx2 = physicsVector->GetValu
184 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2)
185 dedx1 = dedx2;
186 }
187 rangeVector->PutValue(j,range);
188 dedx = dedx1 ;
189 energy1 = energy2 ;
190 }
191 }
192
193 //
194
195 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4Ph
196 G4int nbin)
197 // num. integration, linear binning
198 {
199 G4double dtau,Value,taui,ti,lossi,ci;
200 G4bool isOut;
201 dtau = (tauhigh-taulow)/nbin;
202 Value = 0.;
203
204 for (G4int i=0; i<=nbin; i++)
205 {
206 taui = taulow + dtau*i ;
207 ti = ParticleMass*taui;
208 lossi = physicsVector->GetValue(ti,isOut);
209 if(i==0)
210 ci=0.5;
211 else
212 {
213 if(i<nbin)
214 ci=1.;
215 else
216 ci=0.5;
217 }
218 Value += ci/lossi;
219 }
220 Value *= ParticleMass*dtau;
221 return Value;
222 }
223
224 //
225
226 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4Ph
227 G4int nbin)
228 // num. integration, logarithmic binning
229 {
230 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci
231 G4bool isOut;
232 ltt = ltauhigh-ltaulow;
233 dltau = ltt/nbin;
234 Value = 0.;
235
236 for (G4int i=0; i<=nbin; i++)
237 {
238 ui = ltaulow+dltau*i;
239 taui = std::exp(ui);
240 ti = ParticleMass*taui;
241 lossi = physicsVector->GetValue(ti,isOut);
242 if(i==0)
243 ci=0.5;
244 else
245 {
246 if(i<nbin)
247 ci=1.;
248 else
249 ci=0.5;
250 }
251 Value += ci*taui/lossi;
252 }
253 Value *= ParticleMass*dltau;
254 return Value;
255 }
256
257
258 //
259
260 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabT
261 G4PhysicsTable* theLabTimeTab
262 G4double lowestKineticEnergy,
263 G4double highestKineticEnergy
264
265 {
266
267 G4int numOfCouples = G4ProductionCutsTable::
268
269 if(theLabTimeTable)
270 { theLabTimeTable->clearAndDestroy();
271 delete theLabTimeTable; }
272 theLabTimeTable = new G4PhysicsTable(numOfCo
273
274
275 for (G4int J=0; J<numOfCouples; J++)
276 {
277 G4PhysicsLogVector* aVector;
278
279 aVector = new G4PhysicsLogVector(lowestKin
280 highestKineticEner
281
282 BuildLabTimeVector(theDEDXTable,
283 lowestKineticEnergy,highestKinet
284 theLabTimeTable->insert(aVector);
285
286
287 }
288 return theLabTimeTable ;
289 }
290
291 //
292
293 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProp
294 G4PhysicsTable* theProperTimeTab
295 G4double lowestKineticEnergy,
296 G4double highestKineticEnergy,G4
297
298 {
299
300 G4int numOfCouples = G4ProductionCutsTable::
301
302 if(theProperTimeTable)
303 { theProperTimeTable->clearAndDestroy();
304 delete theProperTimeTable; }
305 theProperTimeTable = new G4PhysicsTable(numO
306
307
308 for (G4int J=0; J<numOfCouples; J++)
309 {
310 G4PhysicsLogVector* aVector;
311
312 aVector = new G4PhysicsLogVector(lowestKin
313 highestKineticEner
314
315 BuildProperTimeVector(theDEDXTable,
316 lowestKineticEnergy,highestKinet
317 theProperTimeTable->insert(aVector);
318
319
320 }
321 return theProperTimeTable ;
322 }
323
324 //
325
326 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G
327 G4double, // lowestKineticEnergy,
328 G4double, // highestKineticEnergy
329 G4i
330 G4int materialIndex, G4PhysicsLog
331 // create lab time vector for a material
332 {
333
334 G4int nbin=100;
335 G4bool isOut;
336 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p
337 G4double losslim,clim,taulim,timelim,
338 LowEdgeEnergy,tau,Value ;
339
340 G4PhysicsVector* physicsVector= (*theDEDXTab
341
342 // low energy part first...
343 losslim = physicsVector->GetValue(tlim,isOut
344 taulim=tlim/ParticleMass ;
345 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh
346
347 G4int i=-1;
348 G4double oldValue = 0. ;
349 G4double tauold ;
350 do
351 {
352 i += 1 ;
353 LowEdgeEnergy = timeVector->GetLowEdgeEner
354 tau = LowEdgeEnergy/ParticleMass ;
355 if ( tau <= taulim )
356 {
357 Value = clim*std::exp(ppar*std::log(tau/
358 }
359 else
360 {
361 timelim=clim ;
362 ltaulow = std::log(taulim);
363 ltauhigh = std::log(tau);
364 Value = timelim+LabTimeIntLog(physicsVec
365 }
366 timeVector->PutValue(i,Value);
367 oldValue = Value ;
368 tauold = tau ;
369 } while (tau<=taulim) ;
370 i += 1 ;
371 for (G4int j=i; j<TotBin; j++)
372 {
373 LowEdgeEnergy = timeVector->GetLowEdgeEner
374 tau = LowEdgeEnergy/ParticleMass ;
375 ltaulow = std::log(tauold);
376 ltauhigh = std::log(tau);
377 Value = oldValue+LabTimeIntLog(physicsVect
378 timeVector->PutValue(j,Value);
379 oldValue = Value ;
380 tauold = tau ;
381 }
382 }
383
384 //
385
386 void G4RDVeLowEnergyLoss::BuildProperTimeVecto
387 G4double, // lowestKineticEner
388 G4double, // highestKineticEne
389
390 G4int materialIndex, G4Physics
391 // create proper time vector for a material
392 {
393 G4int nbin=100;
394 G4bool isOut;
395 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p
396 G4double losslim,clim,taulim,timelim,
397 LowEdgeEnergy,tau,Value ;
398
399 G4PhysicsVector* physicsVector= (*theDEDXTab
400 //const G4MaterialTable* theMaterialTable =
401
402 // low energy part first...
403 losslim = physicsVector->GetValue(tlim,isOut
404 taulim=tlim/ParticleMass ;
405 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh
406
407 G4int i=-1;
408 G4double oldValue = 0. ;
409 G4double tauold ;
410 do
411 {
412 i += 1 ;
413 LowEdgeEnergy = timeVector->GetLowEdgeEner
414 tau = LowEdgeEnergy/ParticleMass ;
415 if ( tau <= taulim )
416 {
417 Value = clim*std::exp(ppar*std::log(tau/
418 }
419 else
420 {
421 timelim=clim ;
422 ltaulow = std::log(taulim);
423 ltauhigh = std::log(tau);
424 Value = timelim+ProperTimeIntLog(physics
425 }
426 timeVector->PutValue(i,Value);
427 oldValue = Value ;
428 tauold = tau ;
429 } while (tau<=taulim) ;
430 i += 1 ;
431 for (G4int j=i; j<TotBin; j++)
432 {
433 LowEdgeEnergy = timeVector->GetLowEdgeEner
434 tau = LowEdgeEnergy/ParticleMass ;
435 ltaulow = std::log(tauold);
436 ltauhigh = std::log(tau);
437 Value = oldValue+ProperTimeIntLog(physicsV
438 timeVector->PutValue(j,Value);
439 oldValue = Value ;
440 tauold = tau ;
441 }
442 }
443
444 //
445
446 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4
447 G4int nbin)
448 // num. integration, logarithmic binning
449 {
450 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci
451 G4bool isOut;
452 ltt = ltauhigh-ltaulow;
453 dltau = ltt/nbin;
454 Value = 0.;
455
456 for (G4int i=0; i<=nbin; i++)
457 {
458 ui = ltaulow+dltau*i;
459 taui = std::exp(ui);
460 ti = ParticleMass*taui;
461 lossi = physicsVector->GetValue(ti,isOut);
462 if(i==0)
463 ci=0.5;
464 else
465 {
466 if(i<nbin)
467 ci=1.;
468 else
469 ci=0.5;
470 }
471 Value += ci*taui*(ti+ParticleMass)/(std::s
472 }
473 Value *= ParticleMass*dltau/c_light;
474 return Value;
475 }
476
477 //
478
479 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog
480 G4int nbin)
481 // num. integration, logarithmic binning
482 {
483 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci
484 G4bool isOut;
485 ltt = ltauhigh-ltaulow;
486 dltau = ltt/nbin;
487 Value = 0.;
488
489 for (G4int i=0; i<=nbin; i++)
490 {
491 ui = ltaulow+dltau*i;
492 taui = std::exp(ui);
493 ti = ParticleMass*taui;
494 lossi = physicsVector->GetValue(ti,isOut);
495 if(i==0)
496 ci=0.5;
497 else
498 {
499 if(i<nbin)
500 ci=1.;
501 else
502 ci=0.5;
503 }
504 Value += ci*taui*ParticleMass/(std::sqrt(t
505 }
506 Value *= ParticleMass*dltau/c_light;
507 return Value;
508 }
509
510 //
511
512 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInve
513 G4PhysicsTable*,
514 G4PhysicsTable*,
515 G4PhysicsTable*,
516 G4PhysicsTable* theInverseRang
517 G4double, // lowestKineticEner
518 G4double, // highestKineticEne
519 G4int ) // nbins
520 // Build inverse table of the range table
521 {
522 G4bool b;
523
524 G4int numOfCouples = G4ProductionCutsTable::
525
526 if(theInverseRangeTable)
527 { theInverseRangeTable->clearAndDestroy();
528 delete theInverseRangeTable; }
529 theInverseRangeTable = new G4PhysicsTable(
530
531 // loop for materials
532 for (G4int i=0; i<numOfCouples; i++)
533 {
534
535 G4PhysicsVector* pv = (*theRangeTable)[i];
536 size_t nbins = pv->GetVectorLength();
537 G4double elow = pv->GetLowEdgeEnergy(0);
538 G4double ehigh = pv->GetLowEdgeEnergy(nbin
539 G4double rlow = pv->GetValue(elow, b);
540 G4double rhigh = pv->GetValue(ehigh, b);
541
542 rhigh *= std::exp(std::log(rhigh/rlow)/((G
543
544 G4PhysicsLogVector* v = new G4PhysicsLogVe
545
546 v->PutValue(0,elow);
547 G4double energy1 = elow;
548 G4double range1 = rlow;
549 G4double energy2 = elow;
550 G4double range2 = rlow;
551 size_t ilow = 0;
552 size_t ihigh;
553
554 for (size_t j=1; j<nbins; j++) {
555
556 G4double range = v->GetLowEdgeEnergy(j);
557
558 for (ihigh=ilow+1; ihigh<nbins; ihigh++)
559 energy2 = pv->GetLowEdgeEnergy(ihigh);
560 range2 = pv->GetValue(energy2, b);
561 if(range2 >= range || ihigh == nbins-1
562 ilow = ihigh - 1;
563 energy1 = pv->GetLowEdgeEnergy(ilow)
564 range1 = pv->GetValue(energy1, b);
565 break;
566 }
567 }
568
569 G4double e = std::log(energy1) + std::lo
570
571 v->PutValue(j,std::exp(e));
572 }
573 theInverseRangeTable->insert(v);
574
575 }
576 return theInverseRangeTable ;
577 }
578
579 //
580
581 void G4RDVeLowEnergyLoss::InvertRangeVector(G4
582 G4PhysicsTable* theRangeCoeffATabl
583 G4PhysicsTable* theRangeCoeffBTabl
584 G4PhysicsTable* theRangeCoeffCTabl
585 G4double lowestKineticEnergy,
586 G4double highestKineticEnergy, G4i
587 G4int materialIndex, G4PhysicsLog
588 // invert range vector for a material
589 {
590 G4double LowEdgeRange,A,B,C,discr,KineticEne
591 G4double RTable = std::exp(std::log(highestK
592 G4double Tbin = lowestKineticEnergy/RTable ;
593 G4double rangebin = 0.0 ;
594 G4int binnumber = -1 ;
595 G4bool isOut ;
596
597 //loop for range values
598 for( G4int i=0; i<TotBin; i++)
599 {
600 LowEdgeRange = aVector->GetLowEdgeEnergy(i
601 if( rangebin < LowEdgeRange )
602 {
603 do
604 {
605 binnumber += 1 ;
606 Tbin *= RTable ;
607 rangebin = (*theRangeTable)(materialIn
608 }
609 while ((rangebin < LowEdgeRange) && (bin
610 }
611
612 if(binnumber == 0)
613 KineticEnergy = lowestKineticEnergy ;
614 else if(binnumber == TotBin-1)
615 KineticEnergy = highestKineticEnergy ;
616 else
617 {
618 A = (*(*theRangeCoeffATable)(materialInd
619 B = (*(*theRangeCoeffBTable)(materialInd
620 C = (*(*theRangeCoeffCTable)(materialInd
621 if(A==0.)
622 KineticEnergy = (LowEdgeRange -C )/B
623 else
624 {
625 discr = B*B - 4.*A*(C-LowEdgeRange);
626 discr = discr>0. ? std::sqrt(discr) :
627 KineticEnergy = 0.5*(discr-B)/A ;
628 }
629 }
630
631 aVector->PutValue(i,KineticEnergy) ;
632 }
633 }
634
635 //
636
637 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang
638 G4PhysicsTable* theRangeCoeffAT
639 G4double lowestKineticEnergy,
640 G4double highestKineticEnergy,
641 // Build tables of coefficients for the energy
642 // create table for coefficients "A"
643 {
644
645 G4int numOfCouples = G4ProductionCutsTable::
646
647 if(theRangeCoeffATable)
648 { theRangeCoeffATable->clearAndDestroy();
649 delete theRangeCoeffATable; }
650 theRangeCoeffATable = new G4PhysicsTable(num
651
652 G4double RTable = std::exp(std::log(highestK
653 G4double R2 = RTable*RTable ;
654 G4double R1 = RTable+1.;
655 G4double w = R1*(RTable-1.)*(RTable-1.);
656 G4double w1 = RTable/w , w2 = -RTable*R1/w ,
657 G4double Ti , Tim , Tip , Ri , Rim , Rip , V
658 G4bool isOut;
659
660 // loop for materials
661 for (G4int J=0; J<numOfCouples; J++)
662 {
663 G4int binmax=TotBin ;
664 G4PhysicsLinearVector* aVector =
665 new G4PhysicsLinear
666 Ti = lowestKineticEnergy ;
667 G4PhysicsVector* rangeVector= (*theRangeTa
668
669 for ( G4int i=0; i<TotBin; i++)
670 {
671 Ri = rangeVector->GetValue(Ti,isOut) ;
672 if ( i==0 )
673 Rim = 0. ;
674 else
675 {
676 Tim = Ti/RTable ;
677 Rim = rangeVector->GetValue(Tim,isOut)
678 }
679 if ( i==(TotBin-1))
680 Rip = Ri ;
681 else
682 {
683 Tip = Ti*RTable ;
684 Rip = rangeVector->GetValue(Tip,isOut)
685 }
686 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti
687
688 aVector->PutValue(i,Value);
689 Ti = RTable*Ti ;
690 }
691
692 theRangeCoeffATable->insert(aVector);
693 }
694 return theRangeCoeffATable ;
695 }
696
697 //
698
699 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang
700 G4PhysicsTable* theRangeCoeffBT
701 G4double lowestKineticEnergy,
702 G4double highestKineticEnergy,
703 // Build tables of coefficients for the energy
704 // create table for coefficients "B"
705 {
706
707 G4int numOfCouples = G4ProductionCutsTable::
708
709 if(theRangeCoeffBTable)
710 { theRangeCoeffBTable->clearAndDestroy();
711 delete theRangeCoeffBTable; }
712 theRangeCoeffBTable = new G4PhysicsTable(num
713
714 G4double RTable = std::exp(std::log(highestK
715 G4double R2 = RTable*RTable ;
716 G4double R1 = RTable+1.;
717 G4double w = R1*(RTable-1.)*(RTable-1.);
718 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3
719 G4double Ti , Tim , Tip , Ri , Rim , Rip , V
720 G4bool isOut;
721
722 // loop for materials
723 for (G4int J=0; J<numOfCouples; J++)
724 {
725 G4int binmax=TotBin ;
726 G4PhysicsLinearVector* aVector =
727 new G4PhysicsLinearVec
728 Ti = lowestKineticEnergy ;
729 G4PhysicsVector* rangeVector= (*theRangeTa
730
731 for ( G4int i=0; i<TotBin; i++)
732 {
733 Ri = rangeVector->GetValue(Ti,isOut) ;
734 if ( i==0 )
735 Rim = 0. ;
736 else
737 {
738 Tim = Ti/RTable ;
739 Rim = rangeVector->GetValue(Tim,isOut)
740 }
741 if ( i==(TotBin-1))
742 Rip = Ri ;
743 else
744 {
745 Tip = Ti*RTable ;
746 Rip = rangeVector->GetValue(Tip,isOut)
747 }
748 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
749
750 aVector->PutValue(i,Value);
751 Ti = RTable*Ti ;
752 }
753 theRangeCoeffBTable->insert(aVector);
754 }
755 return theRangeCoeffBTable ;
756 }
757
758 //
759
760 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang
761 G4PhysicsTable* theRangeCoeffCT
762 G4double lowestKineticEnergy,
763 G4double highestKineticEnergy,
764 // Build tables of coefficients for the energy
765 // create table for coefficients "C"
766 {
767
768 G4int numOfCouples = G4ProductionCutsTable::
769
770 if(theRangeCoeffCTable)
771 { theRangeCoeffCTable->clearAndDestroy();
772 delete theRangeCoeffCTable; }
773 theRangeCoeffCTable = new G4PhysicsTable(num
774
775 G4double RTable = std::exp(std::log(highestK
776 G4double R2 = RTable*RTable ;
777 G4double R1 = RTable+1.;
778 G4double w = R1*(RTable-1.)*(RTable-1.);
779 G4double w1 = 1./w , w2 = -RTable*R1/w , w3
780 G4double Ti , Tim , Tip , Ri , Rim , Rip , V
781 G4bool isOut;
782
783 // loop for materials
784 for (G4int J=0; J<numOfCouples; J++)
785 {
786 G4int binmax=TotBin ;
787 G4PhysicsLinearVector* aVector =
788 new G4PhysicsLinearVecto
789 Ti = lowestKineticEnergy ;
790 G4PhysicsVector* rangeVector= (*theRangeTa
791
792 for ( G4int i=0; i<TotBin; i++)
793 {
794 Ri = rangeVector->GetValue(Ti,isOut) ;
795 if ( i==0 )
796 Rim = 0. ;
797 else
798 {
799 Tim = Ti/RTable ;
800 Rim = rangeVector->GetValue(Tim,isOut)
801 }
802 if ( i==(TotBin-1))
803 Rip = Ri ;
804 else
805 {
806 Tip = Ti*RTable ;
807 Rip = rangeVector->GetValue(Tip,isOut)
808 }
809 Value = w1*Rip + w2*Ri + w3*Rim ;
810
811 aVector->PutValue(i,Value);
812 Ti = RTable*Ti ;
813 }
814 theRangeCoeffCTable->insert(aVector);
815 }
816 return theRangeCoeffCTable ;
817 }
818
819 //
820
821 G4double G4RDVeLowEnergyLoss::GetLossWithFluct
822 c
823 G4double MeanLoss,
824 G4double step)
825 // calculate actual loss from the mean loss
826 // The model used to get the fluctuation is e
827 {
828 static const G4double minLoss = 1.*eV ;
829 static const G4double probLim = 0.01 ;
830 static const G4double sumaLim = -std::log(p
831 static const G4double alim=10.;
832 static const G4double kappa = 10. ;
833 static const G4double factor = twopi_mc2_rc
834 const G4Material* aMaterial = couple->GetMat
835
836 // check if the material has changed ( cache
837
838 if (aMaterial != lastMaterial)
839 {
840 lastMaterial = aMaterial;
841 imat = couple->GetIndex();
842 f1Fluct = aMaterial->GetIonisation(
843 f2Fluct = aMaterial->GetIonisation(
844 e1Fluct = aMaterial->GetIonisation(
845 e2Fluct = aMaterial->GetIonisation(
846 e1LogFluct = aMaterial->GetIonisation(
847 e2LogFluct = aMaterial->GetIonisation(
848 rateFluct = aMaterial->GetIonisation(
849 ipotFluct = aMaterial->GetIonisation(
850 ipotLogFluct = aMaterial->GetIonisation(
851 }
852 G4double threshold,w1,w2,C,
853 beta2,suma,e0,loss,lossc,w;
854 G4double a1,a2,a3;
855 G4int p1,p2,p3;
856 G4int nb;
857 G4double Corrfac, na,alfa,rfac,namean,sa,alf
858 // G4double dp1;
859 G4double dp3;
860 G4double siga ;
861
862 // shortcut for very very small loss
863 if(MeanLoss < minLoss) return MeanLoss ;
864
865 // get particle data
866 G4double Tkin = aParticle->GetKineticEnerg
867
868 // G4cout << "MGP -- Fluc Tkin " << Tkin/ke
869
870 threshold = (*((G4ProductionCutsTable::GetPr
871 ->GetEnergyCutsVector(1)))[ima
872 G4double rmass = electron_mass_c2/ParticleMa
873 G4double tau = Tkin/ParticleMass, tau1 = t
874 G4double Tm = 2.*electron_mass_c2*tau2/(1
875
876 // G4cout << "MGP Particle mass " << Particl
877
878 if(Tm > threshold) Tm = threshold;
879 beta2 = tau2/(tau1*tau1);
880
881 // Gaussian fluctuation ?
882 if(MeanLoss >= kappa*Tm || MeanLoss <= kappa
883 {
884 G4double electronDensity = aMaterial->GetE
885 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
886 factor*electronDensity/beta2)
887 do {
888 loss = G4RandGauss::shoot(MeanLoss,siga)
889 } while (loss < 0. || loss > 2.0*MeanLoss)
890 return loss ;
891 }
892
893 w1 = Tm/ipotFluct;
894 w2 = std::log(2.*electron_mass_c2*tau2);
895
896 C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct
897
898 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct
899 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct
900 a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipot
901
902 suma = a1+a2+a3;
903
904 loss = 0. ;
905
906 if(suma < sumaLim) // very small
907 {
908 e0 = aMaterial->GetIonisation()->GetEner
909 // G4cout << "MGP e0 = " << e0/keV << G4
910
911 if(Tm == ipotFluct)
912 {
913 a3 = MeanLoss/e0;
914
915 if(a3>alim)
916 {
917 siga=std::sqrt(a3) ;
918 p3 = std::max(0,G4int(G4RandGauss::sho
919 }
920 else p3 = G4Poisson(a3);
921
922 loss = p3*e0 ;
923
924 if(p3 > 0) loss += (1.-2.*G4UniformRand())
925 // G4cout << "MGP very small step " << los
926 }
927 else
928 {
929 // G4cout << "MGP old Tm = " << Tm << "
930 Tm = Tm-ipotFluct+e0 ;
931
932 // MGP ---- workaround to avoid log argume
933 if (Tm <= 0.)
934 {
935 loss = MeanLoss;
936 p3 = 0;
937 // G4cout << "MGP correction loss = Me
938 }
939 else
940 {
941 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(
942
943 // G4cout << "MGP new Tm = " << Tm <<
944
945 if(a3>alim)
946 {
947 siga=std::sqrt(a3) ;
948 p3 = std::max(0,G4int(G4RandGauss::shoot
949 }
950 else
951 p3 = G4Poisson(a3);
952 //G4cout << "MGP p3 " << p3 << G4endl;
953
954 }
955
956 if(p3 > 0)
957 {
958 w = (Tm-e0)/Tm ;
959 if(p3 > nmaxCont2)
960 {
961 // G4cout << "MGP dp3 " << dp3 << " p3 "
962 dp3 = G4double(p3) ;
963 Corrfac = dp3/G4double(nmaxCont2) ;
964 p3 = nmaxCont2 ;
965 }
966 else
967 Corrfac = 1. ;
968
969 for(G4int i=0; i<p3; i++) loss += 1./(
970 loss *= e0*Corrfac ;
971 // G4cout << "MGP Corrfac = " << Corrf
972 }
973 }
974 }
975
976 else // not so
977 {
978 // excitation type 1
979 if(a1>alim)
980 {
981 siga=std::sqrt(a1) ;
982 p1 = std::max(0,int(G4RandGauss::shoot
983 }
984 else
985 p1 = G4Poisson(a1);
986
987 // excitation type 2
988 if(a2>alim)
989 {
990 siga=std::sqrt(a2) ;
991 p2 = std::max(0,int(G4RandGauss::shoot
992 }
993 else
994 p2 = G4Poisson(a2);
995
996 loss = p1*e1Fluct+p2*e2Fluct;
997
998 // smearing to avoid unphysical peaks
999 if(p2 > 0)
1000 loss += (1.-2.*G4UniformRand())*e2Flu
1001 else if (loss>0.)
1002 loss += (1.-2.*G4UniformRand())*e1Flu
1003
1004 // ionisation .........................
1005 if(a3 > 0.)
1006 {
1007 if(a3>alim)
1008 {
1009 siga=std::sqrt(a3) ;
1010 p3 = std::max(0,int(G4RandGauss::shoo
1011 }
1012 else
1013 p3 = G4Poisson(a3);
1014
1015 lossc = 0.;
1016 if(p3 > 0)
1017 {
1018 na = 0.;
1019 alfa = 1.;
1020 if (p3 > nmaxCont2)
1021 {
1022 dp3 = G4double(p3);
1023 rfac = dp3/(G4double(nmaxCont2)+d
1024 namean = G4double(p3)*rfac;
1025 sa = G4double(nmaxCont1)*rfac;
1026 na = G4RandGauss::shoot(namean,
1027 if (na > 0.)
1028 {
1029 alfa = w1*G4double(nmaxCont2+p3)/
1030 alfa1 = alfa*std::log(alfa)/(alfa-
1031 ea = na*ipotFluct*alfa1;
1032 sea = ipotFluct*std::sqrt(na*(al
1033 lossc += G4RandGauss::shoot(ea,sea)
1034 }
1035 }
1036
1037 nb = G4int(G4double(p3)-na);
1038 if (nb > 0)
1039 {
1040 w2 = alfa*ipotFluct;
1041 w = (Tm-w2)/Tm;
1042 for (G4int k=0; k<nb; k++) lossc += w2/
1043 }
1044 }
1045
1046 loss += lossc;
1047 }
1048 }
1049
1050 return loss ;
1051 }
1052
1053 //
1054
1055