Geant4  10.02.p01
G4UniversalFluctuation.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: G4UniversalFluctuation.cc 91726 2015-08-03 15:41:36Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4UniversalFluctuation
34 //
35 // Author: Laszlo Urban
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 // 16-10-03 Changed interface to Initialisation (V.Ivanchenko)
45 // 07-11-03 Fix problem of rounding of double in G4UniversalFluctuations
46 // 06-02-04 Add control on big sigma > 2*meanLoss (V.Ivanchenko)
47 // 26-04-04 Comment out the case of very small step (V.Ivanchenko)
48 // 07-02-05 define problim = 5.e-3 (mma)
49 // 03-05-05 conditions of Gaussian fluctuation changed (bugfix)
50 // + smearing for very small loss (L.Urban)
51 // 03-10-05 energy dependent rate -> cut dependence of the
52 // distribution is much weaker (L.Urban)
53 // 17-10-05 correction for very small loss (L.Urban)
54 // 20-03-07 'GLANDZ' part rewritten completely, no 'very small loss'
55 // regime any more (L.Urban)
56 // 03-04-07 correction to get better width of eloss distr.(L.Urban)
57 // 13-07-07 add protection for very small step or low-density material (VI)
58 // 19-03-09 new width correction (does not depend on previous steps) (L.Urban)
59 // 20-03-09 modification in the width correction (L.Urban)
60 // 14-06-10 fixed tail distribution - do not use uniform function (L.Urban)
61 // 08-08-10 width correction algorithm has bee modified -->
62 // better results for thin targets (L.Urban)
63 // 06-02-11 correction for very small losses (L.Urban)
64 //
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 
70 #include "G4PhysicalConstants.hh"
71 #include "G4SystemOfUnits.hh"
72 #include "Randomize.hh"
73 #include "G4Poisson.hh"
74 #include "G4Step.hh"
75 #include "G4Material.hh"
76 #include "G4MaterialCutsCouple.hh"
77 #include "G4DynamicParticle.hh"
78 #include "G4ParticleDefinition.hh"
79 #include "G4Log.hh"
80 #include "G4Exp.hh"
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 
84 using namespace std;
85 
88  particle(0),
89  minNumberInteractionsBohr(10.0),
90  theBohrBeta2(50.0*keV/proton_mass_c2),
91  minLoss(10.*eV),
92  nmaxCont(16.),
93  rate(0.55),
94  fw(4.)
95 {
96  lastMaterial = 0;
97 
100  = e1 = e2 = 0;
102  sizearray = 30;
103  rndmarray = new G4double[30];
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 
109 {
110  delete [] rndmarray;
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
116 {
117  particle = part;
118  particleMass = part->GetPDGMass();
119  G4double q = part->GetPDGCharge()/eplus;
120 
121  // Derived quantities
123  m_massrate = electron_mass_c2 * m_Inv_particleMass ;
124  chargeSquare = q*q;
125 }
126 
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 
129 G4double
131  const G4DynamicParticle* dp,
132  G4double tmax,
133  G4double length,
134  G4double averageLoss)
135 {
136  // Calculate actual loss from the mean loss.
137  // The model used to get the fluctuations is essentially the same
138  // as in Glandz in Geant3 (Cern program library W5013, phys332).
139  // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
140 
141  // shortcut for very small loss or from a step nearly equal to the range
142  // (out of validity of the model)
143  //
144  G4double meanLoss = averageLoss;
145  G4double tkin = dp->GetKineticEnergy();
146  //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
147  if (meanLoss < minLoss) { return meanLoss; }
148 
149  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
150 
151  CLHEP::HepRandomEngine* rndmEngineF = G4Random::getTheEngine();
152 
153  G4double tau = tkin * m_Inv_particleMass;
154  G4double gam = tau + 1.0;
155  G4double gam2 = gam*gam;
156  G4double beta2 = tau*(tau + 2.0)/gam2;
157 
158  G4double loss(0.), siga(0.);
159 
160  const G4Material* material = couple->GetMaterial();
161 
162  // Gaussian regime
163  // for heavy particles only and conditions
164  // for Gauusian fluct. has been changed
165  //
166  if ((particleMass > electron_mass_c2) &&
167  (meanLoss >= minNumberInteractionsBohr*tmax))
168  {
169  G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
170  (1.+m_massrate*(2.*gam+m_massrate)) ;
171  if (tmaxkine <= 2.*tmax)
172  {
173  electronDensity = material->GetElectronDensity();
174  siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
176 
177  G4double sn = meanLoss/siga;
178 
179  // thick target case
180  if (sn >= 2.0) {
181 
182  G4double twomeanLoss = meanLoss + meanLoss;
183  do {
184  loss = G4RandGauss::shoot(rndmEngineF,meanLoss,siga);
185  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
186  } while (0.0 > loss || twomeanLoss < loss);
187 
188  // Gamma distribution
189  } else {
190 
191  G4double neff = sn*sn;
192  loss = meanLoss*G4RandGamma::shoot(rndmEngineF,neff,1.0)/neff;
193  }
194  //G4cout << "Gauss: " << loss << G4endl;
195  return loss;
196  }
197  }
198 
199  // Glandz regime : initialisation
200  //
201  if (material != lastMaterial) {
202  f1Fluct = material->GetIonisation()->GetF1fluct();
203  f2Fluct = material->GetIonisation()->GetF2fluct();
204  e1Fluct = material->GetIonisation()->GetEnergy1fluct();
205  e2Fluct = material->GetIonisation()->GetEnergy2fluct();
210  e0 = material->GetIonisation()->GetEnergy0fluct();
211  esmall = 0.5*sqrt(e0*ipotFluct);
212  lastMaterial = material;
213  }
214 
215  // very small step or low-density material
216  if(tmax <= e0) { return meanLoss; }
217 
218  G4double losstot = 0.;
219  G4int nstep = 1;
220  if(meanLoss < 25.*ipotFluct)
221  {
222  if(rndmEngineF->flat()*ipotFluct< 0.04*meanLoss)
223  { nstep = 1; }
224  else
225  {
226  nstep = 2;
227  meanLoss *= 0.5;
228  }
229  }
230 
231  for (G4int istep=0; istep < nstep; ++istep) {
232 
233  loss = 0.;
234 
235  G4double a1 = 0. , a2 = 0., a3 = 0. ;
236 
237  if(tmax > ipotFluct) {
238  G4double w2 = G4Log(2.*electron_mass_c2*beta2*gam2)-beta2;
239 
240  if(w2 > ipotLogFluct) {
241  G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
242  a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
243  if(w2 > e2LogFluct) {
244  a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
245  }
246  if(a1 < nmaxCont) {
247  //small energy loss
248  G4double sa1 = sqrt(a1);
249  if(rndmEngineF->flat() < G4Exp(-sa1))
250  {
251  e1 = esmall;
252  a1 = meanLoss*(1.-rate)/e1;
253  a2 = 0.;
254  e2 = e2Fluct;
255  }
256  else
257  {
258  a1 = sa1 ;
259  e1 = sa1*e1Fluct;
260  e2 = e2Fluct;
261  }
262 
263  } else {
264  //not small energy loss
265  //correction to get better fwhm value
266  a1 /= fw;
267  e1 = fw*e1Fluct;
268  e2 = e2Fluct;
269  }
270  }
271  }
272 
273  G4double w1 = tmax/e0;
274  if(tmax > e0) {
275  a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*G4Log(w1));
276  if(a1+a2 <= 0.) {
277  a3 /= rate;
278  }
279  }
280  //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont
281  G4double emean = 0.;
282  G4double sig2e = 0., sige = 0.;
283  G4double p1 = 0., p2 = 0., p3 = 0.;
284 
285  // excitation of type 1
286  if(a1 > nmaxCont)
287  {
288  emean += a1*e1;
289  sig2e += a1*e1*e1;
290  }
291  else if(a1 > 0.)
292  {
293  p1 = G4double(G4Poisson(a1));
294  loss += p1*e1;
295  if(p1 > 0.) {
296  loss += (1.-2.*rndmEngineF->flat())*e1;
297  }
298  }
299 
300  // excitation of type 2
301  if(a2 > nmaxCont)
302  {
303  emean += a2*e2;
304  sig2e += a2*e2*e2;
305  }
306  else if(a2 > 0.)
307  {
308  p2 = G4double(G4Poisson(a2));
309  loss += p2*e2;
310  if(p2 > 0.)
311  loss += (1.-2.*rndmEngineF->flat())*e2;
312  }
313  if(emean > 0.)
314  {
315  sige = sqrt(sig2e);
316  loss += max(0.,G4RandGauss::shoot(rndmEngineF,emean,sige));
317  }
318 
319  // ionisation
320  G4double lossc = 0.;
321  if(a3 > 0.) {
322  emean = 0.;
323  sig2e = 0.;
324  sige = 0.;
325  p3 = a3;
326  G4double alfa = 1.;
327  if(a3 > nmaxCont)
328  {
329  alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
330  G4double alfa1 = alfa*G4Log(alfa)/(alfa-1.);
331  G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
332  emean += namean*e0*alfa1;
333  sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
334  p3 = a3-namean;
335  }
336 
337  G4double w2 = alfa*e0;
338  G4double w = (tmax-w2)/tmax;
339  const G4int nb = G4Poisson(p3);
340  if(nb > 0) {
341  if(nb > sizearray) {
342  sizearray = nb;
343  delete [] rndmarray;
344  rndmarray = new G4double[nb];
345  }
346  rndmEngineF->flatArray(nb, rndmarray);
347  for (G4int k=0; k<nb; ++k) { lossc += w2/(1.-w*rndmarray[k]); }
348  }
349 
350  if(emean > 0.)
351  {
352  sige = sqrt(sig2e);
353  lossc += max(0.,G4RandGauss::shoot(rndmEngineF,emean,sige));
354  }
355  }
356 
357  loss += lossc;
358 
359  losstot += loss;
360  }
361  //G4cout << "Vavilov: " << losstot << " Nstep= " << nstep << G4endl;
362 
363  return losstot;
364 
365 }
366 
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
368 
369 
371  const G4Material* material,
372  const G4DynamicParticle* dp,
373  G4double tmax,
374  G4double length)
375 {
376  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
377 
378  electronDensity = material->GetElectronDensity();
379 
381  G4double beta2 = 1.0 - 1.0/(gam*gam);
382 
383  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
385 
386  return siga;
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390 
391 void
393  G4double q2)
394 {
395  if(part != particle) {
396  particle = part;
397  particleMass = part->GetPDGMass();
398 
399  // Derived quantities
400  if( particleMass != 0.0 ){
402  m_massrate = electron_mass_c2 * m_Inv_particleMass ;
403  }else{
406  }
407  }
408  chargeSquare = q2;
409 }
410 
411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ThreeVector shoot(const G4int Ap, const G4int Af)
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetKineticEnergy() const
G4double GetEnergy2fluct() const
static const G4double a1
const G4Material * lastMaterial
const G4double w[NPOINTSGL]
G4double GetLogEnergy2fluct() const
G4UniversalFluctuation(const G4String &nam="UniFluc")
G4ParticleDefinition * GetDefinition() const
double C(double temp)
G4double GetLogMeanExcEnergy() const
int G4int
Definition: G4Types.hh:78
virtual G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double, G4double)
G4double GetEnergy0fluct() const
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
G4double GetElectronDensity() const
Definition: G4Material.hh:217
virtual void InitialiseMe(const G4ParticleDefinition *)
const G4ParticleDefinition * particle
static const G4double a3
virtual G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double, G4double)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const double eV
Definition: G4SIunits.hh:212
G4double GetLogEnergy1fluct() const
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetMeanExcitationEnergy() const
G4double GetF2fluct() const
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:196
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
G4double GetF1fluct() const
static const G4double a2
const G4Material * GetMaterial() const
G4double GetEnergy1fluct() const