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