Geant4  10.02.p02
G4IonFluctuations.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 // $Id: G4IonFluctuations.cc 93567 2015-10-26 14:51:41Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4IonFluctuation
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 03.01.2002
38 //
39 // Modifications:
40 //
41 // 28-12-02 add method Dispersion (V.Ivanchenko)
42 // 07-02-03 change signature (V.Ivanchenko)
43 // 13-02-03 Add name (V.Ivanchenko)
44 // 23-05-03 Add control on parthalogical cases (V.Ivanchenko)
45 // 16-10-03 Changed interface to Initialisation (V.Ivanchenko)
46 // 27-09-07 Use FermiEnergy from material, add cut dependence (V.Ivanchenko)
47 // 01-02-08 Add protection for small energies and optimise the code (V.Ivanchenko)
48 // 01-06-08 Added initialisation of effective charge prestep (V.Ivanchenko)
49 //
50 // Class Description:
51 //
52 // -------------------------------------------------------------------
53 //
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57 
58 #include "G4IonFluctuations.hh"
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "Randomize.hh"
62 #include "G4Poisson.hh"
63 #include "G4MaterialCutsCouple.hh"
64 #include "G4Material.hh"
65 #include "G4DynamicParticle.hh"
66 #include "G4Pow.hh"
67 #include "G4Log.hh"
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
71 using namespace std;
72 
74  : G4VEmFluctuationModel(nam),
75  particle(0),
76  particleMass(proton_mass_c2),
77  charge(1.0),
78  chargeSquare(1.0),
79  effChargeSquare(1.0),
80  parameter(10.0*CLHEP::MeV/CLHEP::proton_mass_c2),
81  minNumberInteractionsBohr(0.0),
82  theBohrBeta2(50.0*keV/CLHEP::proton_mass_c2),
83  minFraction(0.2),
84  xmin(0.2),
85  minLoss(0.001*eV)
86 {
87  kineticEnergy = 0.0;
88  beta2 = 0.0;
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 
95 {}
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
100 {
101  particle = part;
102  particleMass = part->GetPDGMass();
103  charge = part->GetPDGCharge()/eplus;
106  uniFluct.InitialiseMe(part);
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
111 G4double
113  const G4DynamicParticle* dp,
114  G4double tmax,
115  G4double length,
116  G4double meanLoss)
117 {
118  // G4cout << "### meanLoss= " << meanLoss << G4endl;
119  if(meanLoss <= minLoss) return meanLoss;
120 
121  //G4cout << "G4IonFluctuations::SampleFluctuations E(MeV)= "
122  // << dp->GetKineticEnergy()
123  // << " Elim(MeV)= " << parameter*charge*particleMass << G4endl;
124 
125  // Vavilov fluctuations
127  return uniFluct.SampleFluctuations(couple,dp,tmax,length,meanLoss);
128  }
129 
130  const G4Material* material = couple->GetMaterial();
131  G4double siga = Dispersion(material,dp,tmax,length);
132  G4double loss = meanLoss;
133 
134  //G4cout << "### siga= " << sqrt(siga) << " navr= " << navr << G4endl;
135 
136  // Gaussian fluctuation
137 
138  // Increase fluctuations for big fractional energy loss
139  //G4cout << "siga= " << siga << G4endl;
140  if ( meanLoss > minFraction*kineticEnergy ) {
141  G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0;
142  G4double b2 = 1.0 - 1.0/(gam*gam);
143  if(b2 < xmin*beta2) b2 = xmin*beta2;
144  G4double x = b2/beta2;
145  G4double x3 = 1.0/(x*x*x);
146  siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) );
147  }
148  siga = sqrt(siga);
149  G4double sn = meanLoss/siga;
150  G4double twomeanLoss = meanLoss + meanLoss;
151  // G4cout << "siga= " << siga << " sn= " << sn << G4endl;
152 
153  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
154  // thick target case
155  if (sn >= 2.0) {
156 
157  do {
158  loss = G4RandGauss::shoot(rndmEngine,meanLoss,siga);
159  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
160  } while (0.0 > loss || twomeanLoss < loss);
161 
162  // Gamma distribution
163  } else if(sn > 0.1) {
164 
165  G4double neff = sn*sn;
166  loss = meanLoss*G4RandGamma::shoot(rndmEngine,neff,1.0)/neff;
167 
168  // uniform distribution for very small steps
169  } else {
170  loss = twomeanLoss*rndmEngine->flat();
171  }
172 
173  //G4cout << "meanLoss= " << meanLoss << " loss= " << loss << G4endl;
174  return loss;
175 }
176 
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178 
180  const G4DynamicParticle* dp,
181  G4double tmax,
182  G4double length)
183 {
186  beta2 = kineticEnergy*(kineticEnergy + 2.*particleMass)/(etot*etot);
187 
188  G4double electronDensity = material->GetElectronDensity();
189 
190  /*
191  G4cout << "e= " << kineticEnergy << " m= " << particleMass
192  << " tmax= " << tmax << " l= " << length
193  << " q^2= " << effChargeSquare << " beta2=" << beta2<< G4endl;
194  */
195  G4double siga = (1. - beta2*0.5)*tmax*length*electronDensity*
196  twopi_mc2_rcl2*chargeSquare/beta2;
197 
198  // Low velocity - additional ion charge fluctuations according to
199  // Q.Yang et al., NIM B61(1991)149-155.
200  //G4cout << "sigE= " << sqrt(siga) << " charge= " << charge <<G4endl;
201 
202  G4double Z = material->GetIonisation()->GetZeffective();
203 
204  G4double fac = Factor(material, Z);
205 
206  // heavy ion correction
207  // G4double f1 = 1.065e-4*chargeSquare;
208  // if(beta2 > theBohrBeta2) f1/= beta2;
209  // else f1/= theBohrBeta2;
210  // if(f1 > 2.5) f1 = 2.5;
211  // fac *= (1.0 + f1);
212 
213  // taking into account the cut
214  G4double fac_cut = 1.0 + (fac - 1.0)*2.0*electron_mass_c2*beta2
215  /(tmax*(1.0 - beta2));
216  if(fac_cut > 0.01 && fac > 0.01) {
217  siga *= fac_cut;
218  }
219 
220  //G4cout << "siga(keV)= " << sqrt(siga)/keV << " fac= " << fac
221  // << " f1= " << f1 << G4endl;
222 
223  return siga;
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227 
229 {
230  // The aproximation of energy loss fluctuations
231  // Q.Yang et al., NIM B61(1991)149-155.
232 
233  // Reduced energy in MeV/AMU
235 
236  // simple approximation for higher beta2
237  G4double s1 = RelativisticFactor(material, Z);
238 
239  // tabulation for lower beta2
240  if( beta2 < 3.0*theBohrBeta2*Z ) {
241 
242  static const G4double a[96][4] = {
243  {-0.3291, -0.8312, 0.2460, -1.0220},
244  {-0.5615, -0.5898, 0.5205, -0.7258},
245  {-0.5280, -0.4981, 0.5519, -0.5865},
246  {-0.5125, -0.4625, 0.5660, -0.5190},
247  {-0.5127, -0.8595, 0.5626, -0.8721},
248  {-0.5174, -1.1930, 0.5565, -1.1980},
249  {-0.5179, -1.1850, 0.5560, -1.2070},
250  {-0.5209, -0.9355, 0.5590, -1.0250},
251  {-0.5255, -0.7766, 0.5720, -0.9412},
252 
253  {-0.5776, -0.6665, 0.6598, -0.8484},
254  {-0.6013, -0.6045, 0.7321, -0.7671},
255  {-0.5781, -0.5518, 0.7605, -0.6919},
256  {-0.5587, -0.4981, 0.7835, -0.6195},
257  {-0.5466, -0.4656, 0.7978, -0.5771},
258  {-0.5406, -0.4690, 0.8031, -0.5718},
259  {-0.5391, -0.5061, 0.8024, -0.5974},
260  {-0.5380, -0.6483, 0.7962, -0.6970},
261  {-0.5355, -0.7722, 0.7962, -0.7839},
262  {-0.5329, -0.7720, 0.7988, -0.7846},
263 
264  {-0.5335, -0.7671, 0.7984, -0.7933},
265  {-0.5324, -0.7612, 0.7998, -0.8031},
266  {-0.5305, -0.7300, 0.8031, -0.7990},
267  {-0.5307, -0.7178, 0.8049, -0.8216},
268  {-0.5248, -0.6621, 0.8165, -0.7919},
269  {-0.5180, -0.6502, 0.8266, -0.7986},
270  {-0.5084, -0.6408, 0.8396, -0.8048},
271  {-0.4967, -0.6331, 0.8549, -0.8093},
272  {-0.4861, -0.6508, 0.8712, -0.8432},
273  {-0.4700, -0.6186, 0.8961, -0.8132},
274 
275  {-0.4545, -0.5720, 0.9227, -0.7710},
276  {-0.4404, -0.5226, 0.9481, -0.7254},
277  {-0.4288, -0.4778, 0.9701, -0.6850},
278  {-0.4199, -0.4425, 0.9874, -0.6539},
279  {-0.4131, -0.4188, 0.9998, -0.6332},
280  {-0.4089, -0.4057, 1.0070, -0.6218},
281  {-0.4039, -0.3913, 1.0150, -0.6107},
282  {-0.3987, -0.3698, 1.0240, -0.5938},
283  {-0.3977, -0.3608, 1.0260, -0.5852},
284  {-0.3972, -0.3600, 1.0260, -0.5842},
285 
286  {-0.3985, -0.3803, 1.0200, -0.6013},
287  {-0.3985, -0.3979, 1.0150, -0.6168},
288  {-0.3968, -0.3990, 1.0160, -0.6195},
289  {-0.3971, -0.4432, 1.0050, -0.6591},
290  {-0.3944, -0.4665, 1.0010, -0.6825},
291  {-0.3924, -0.5109, 0.9921, -0.7235},
292  {-0.3882, -0.5158, 0.9947, -0.7343},
293  {-0.3838, -0.5125, 0.9999, -0.7370},
294  {-0.3786, -0.4976, 1.0090, -0.7310},
295  {-0.3741, -0.4738, 1.0200, -0.7155},
296 
297  {-0.3969, -0.4496, 1.0320, -0.6982},
298  {-0.3663, -0.4297, 1.0430, -0.6828},
299  {-0.3630, -0.4120, 1.0530, -0.6689},
300  {-0.3597, -0.3964, 1.0620, -0.6564},
301  {-0.3555, -0.3809, 1.0720, -0.6454},
302  {-0.3525, -0.3607, 1.0820, -0.6289},
303  {-0.3505, -0.3465, 1.0900, -0.6171},
304  {-0.3397, -0.3570, 1.1020, -0.6384},
305  {-0.3314, -0.3552, 1.1130, -0.6441},
306  {-0.3235, -0.3531, 1.1230, -0.6498},
307 
308  {-0.3150, -0.3483, 1.1360, -0.6539},
309  {-0.3060, -0.3441, 1.1490, -0.6593},
310  {-0.2968, -0.3396, 1.1630, -0.6649},
311  {-0.2935, -0.3225, 1.1760, -0.6527},
312  {-0.2797, -0.3262, 1.1940, -0.6722},
313  {-0.2704, -0.3202, 1.2100, -0.6770},
314  {-0.2815, -0.3227, 1.2480, -0.6775},
315  {-0.2880, -0.3245, 1.2810, -0.6801},
316  {-0.3034, -0.3263, 1.3270, -0.6778},
317  {-0.2936, -0.3215, 1.3430, -0.6835},
318 
319  {-0.3282, -0.3200, 1.3980, -0.6650},
320  {-0.3260, -0.3070, 1.4090, -0.6552},
321  {-0.3511, -0.3074, 1.4470, -0.6442},
322  {-0.3501, -0.3064, 1.4500, -0.6442},
323  {-0.3490, -0.3027, 1.4550, -0.6418},
324  {-0.3487, -0.3048, 1.4570, -0.6447},
325  {-0.3478, -0.3074, 1.4600, -0.6483},
326  {-0.3501, -0.3283, 1.4540, -0.6669},
327  {-0.3494, -0.3373, 1.4550, -0.6765},
328  {-0.3485, -0.3373, 1.4570, -0.6774},
329 
330  {-0.3462, -0.3300, 1.4630, -0.6728},
331  {-0.3462, -0.3225, 1.4690, -0.6662},
332  {-0.3453, -0.3094, 1.4790, -0.6553},
333  {-0.3844, -0.3134, 1.5240, -0.6412},
334  {-0.3848, -0.3018, 1.5310, -0.6303},
335  {-0.3862, -0.2955, 1.5360, -0.6237},
336  {-0.4262, -0.2991, 1.5860, -0.6115},
337  {-0.4278, -0.2910, 1.5900, -0.6029},
338  {-0.4303, -0.2817, 1.5940, -0.5927},
339  {-0.4315, -0.2719, 1.6010, -0.5829},
340 
341  {-0.4359, -0.2914, 1.6050, -0.6010},
342  {-0.4365, -0.2982, 1.6080, -0.6080},
343  {-0.4253, -0.3037, 1.6120, -0.6150},
344  {-0.4335, -0.3245, 1.6160, -0.6377},
345  {-0.4307, -0.3292, 1.6210, -0.6447},
346  {-0.4284, -0.3204, 1.6290, -0.6380},
347  {-0.4227, -0.3217, 1.6360, -0.6438}
348  } ;
349 
350  G4int iz = G4lrint(Z) - 2;
351  if( 0 > iz ) { iz = 0; }
352  else if(95 < iz ) { iz = 95; }
353 
354  // G4double ss = 1.0 + a[iz][0]*pow(energy,a[iz][1])+
355  // + a[iz][2]*pow(energy,a[iz][3]);
356  G4double ss = 1.0 + a[iz][0]*g4pow->powA(energy,a[iz][1])+
357  + a[iz][2]*g4pow->powA(energy,a[iz][3]);
358 
359  // protection for the validity range for low beta
360  static const G4double slim = 0.001;
361  if(ss < slim) { s1 = 1.0/slim; }
362  // for high value of beta
363  else if(s1*ss < 1.0) { s1 = 1.0/ss; }
364  }
365  G4int i = 0 ;
366  G4double factor = 1.0 ;
367 
368  // The index of set of parameters i = 0 for protons(hadrons) in gases
369  // 1 for protons(hadrons) in solids
370  // 2 for ions in atomic gases
371  // 3 for ions in molecular gases
372  // 4 for ions in solids
373  static const G4double b[5][4] = {
374  {0.1014, 0.3700, 0.9642, 3.987},
375  {0.1955, 0.6941, 2.522, 1.040},
376  {0.05058, 0.08975, 0.1419, 10.80},
377  {0.05009, 0.08660, 0.2751, 3.787},
378  {0.01273, 0.03458, 0.3951, 3.812}
379  } ;
380 
381  // protons (hadrons)
382  if(1.5 > charge) {
383  if( kStateGas != material->GetState() ) i = 1 ;
384 
385  // ions
386  } else {
387 
388  // factor = charge * pow(charge/Z, 0.33333333);
389  factor = charge * g4pow->A13(charge/Z);
390 
391  if( kStateGas == material->GetState() ) {
392  energy /= (charge * sqrt(charge)) ;
393 
394  if(1 == (material->GetNumberOfElements())) {
395  i = 2 ;
396  } else {
397  i = 3 ;
398  }
399 
400  } else {
401  energy /= (charge * sqrt(charge*Z)) ;
402  i = 4 ;
403  }
404  }
405 
406  G4double x = b[i][2];
407  G4double y = energy * b[i][3];
408  if(y <= 0.2) x *= (y*(1.0 - 0.5*y));
409  else x *= (1.0 - g4pow->expA(-y));
410  // else x *= (1.0 - exp(-y));
411 
412  y = energy - b[i][1];
413 
414  G4double s2 = factor * x * b[i][0] / (y*y + x*x);
415  /*
416  G4cout << "s1= " << s1 << " s2= " << s2 << " q^2= " << effChargeSquare
417  << " e= " << energy << G4endl;
418  */
419  return s1*effChargeSquare/chargeSquare + s2;
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
423 
425  G4double Z)
426 {
427  G4double eF = mat->GetIonisation()->GetFermiEnergy();
429 
430  // H.Geissel et al. NIM B, 195 (2002) 3.
431  G4double bF2= 2.0*eF/electron_mass_c2;
432  G4double f = 0.4*(1.0 - beta2)/((1.0 - 0.5*beta2)*Z);
433  if(beta2 > bF2) f *= G4Log(2.0*electron_mass_c2*beta2/I)*bF2/beta2;
434  else f *= G4Log(4.0*eF/I);
435 
436  // G4cout << "f= " << f << " beta2= " << beta2
437  // << " bf2= " << bF2 << " q^2= " << chargeSquare << G4endl;
438 
439  return 1.0 + f;
440 }
441 
442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
443 
445  G4double q2)
446 {
447  if(part != particle) {
448  particle = part;
449  particleMass = part->GetPDGMass();
450  charge = part->GetPDGCharge()/eplus;
452  }
453  effChargeSquare = q2;
454  uniFluct.SetParticleAndCharge(part, q2);
455 }
456 
457 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
ThreeVector shoot(const G4int Ap, const G4int Af)
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
static const double MeV
Definition: G4SIunits.hh:211
G4double GetKineticEnergy() const
G4double expA(G4double A) const
Definition: G4Pow.hh:235
G4IonFluctuations(const G4String &nam="IonFluc")
void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmax, G4double length, G4double meanLoss)
G4double a
Definition: TRTMaterials.hh:39
void InitialiseMe(const G4ParticleDefinition *)
int G4int
Definition: G4Types.hh:78
G4double GetZeffective() const
G4double RelativisticFactor(const G4Material *, G4double Zeff)
G4double GetFermiEnergy() const
G4UniversalFluctuation uniFluct
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double tmax, G4double length)
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4double iz
Definition: TRTMaterials.hh:39
static const G4double b2
virtual void InitialiseMe(const G4ParticleDefinition *)
virtual ~G4IonFluctuations()
virtual G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double, G4double)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
static const double eV
Definition: G4SIunits.hh:212
static const G4double factor
G4double GetPDGMass() const
G4double A13(G4double A) const
Definition: G4Pow.hh:132
int G4lrint(double ad)
Definition: templates.hh:163
G4double energy(const ThreeVector &p, const G4double m)
const G4double x[NPOINTSGL]
G4double Factor(const G4Material *, G4double Zeff)
G4double GetMeanExcitationEnergy() const
const G4ParticleDefinition * particle
static const G4double fac
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double keV
Definition: G4SIunits.hh:213
G4State GetState() const
Definition: G4Material.hh:181
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:196
G4double GetPDGCharge() const
const G4Material * GetMaterial() const