Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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$
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 
69 #include "G4UniversalFluctuation.hh"
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 "G4DynamicParticle.hh"
77 #include "G4ParticleDefinition.hh"
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 
81 using namespace std;
82 
83 G4UniversalFluctuation::G4UniversalFluctuation(const G4String& nam)
84  :G4VEmFluctuationModel(nam),
85  particle(0),
86  minNumberInteractionsBohr(10.0),
87  theBohrBeta2(50.0*keV/proton_mass_c2),
88  minLoss(10.*eV),
89  nmaxCont(16.),
90  rate(0.55),
91  fw(4.)
92 {
93  lastMaterial = 0;
94 
95  particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct
96  = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall
97  = e1 = e2 = 0;
98 
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 
103 G4UniversalFluctuation::~G4UniversalFluctuation()
104 {}
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 
108 void G4UniversalFluctuation::InitialiseMe(const G4ParticleDefinition* part)
109 {
110  particle = part;
111  particleMass = part->GetPDGMass();
112  G4double q = part->GetPDGCharge()/eplus;
113  chargeSquare = q*q;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 G4double G4UniversalFluctuation::SampleFluctuations(const G4Material* material,
119  const G4DynamicParticle* dp,
120  G4double& tmax,
121  G4double& length,
122  G4double& averageLoss)
123 {
124  // Calculate actual loss from the mean loss.
125  // The model used to get the fluctuations is essentially the same
126  // as in Glandz in Geant3 (Cern program library W5013, phys332).
127  // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
128 
129  // shortcut for very small loss or from a step nearly equal to the range
130  // (out of validity of the model)
131  //
132  G4double meanLoss = averageLoss;
133  G4double tkin = dp->GetKineticEnergy();
134  //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
135  if (meanLoss < minLoss) { return meanLoss; }
136 
137  if(!particle) { InitialiseMe(dp->GetDefinition()); }
138 
139  G4double tau = tkin/particleMass;
140  G4double gam = tau + 1.0;
141  G4double gam2 = gam*gam;
142  G4double beta2 = tau*(tau + 2.0)/gam2;
143 
144  G4double loss(0.), siga(0.);
145 
146  // Gaussian regime
147  // for heavy particles only and conditions
148  // for Gauusian fluct. has been changed
149  //
150  if ((particleMass > electron_mass_c2) &&
151  (meanLoss >= minNumberInteractionsBohr*tmax))
152  {
153  G4double massrate = electron_mass_c2/particleMass ;
154  G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
155  (1.+massrate*(2.*gam+massrate)) ;
156  if (tmaxkine <= 2.*tmax)
157  {
158  electronDensity = material->GetElectronDensity();
159  siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
160  * electronDensity * chargeSquare);
161 
162 
163  G4double sn = meanLoss/siga;
164 
165  // thick target case
166  if (sn >= 2.0) {
167 
168  G4double twomeanLoss = meanLoss + meanLoss;
169  do {
170  loss = G4RandGauss::shoot(meanLoss,siga);
171  } while (0.0 > loss || twomeanLoss < loss);
172 
173  // Gamma distribution
174  } else {
175 
176  G4double neff = sn*sn;
177  loss = meanLoss*CLHEP::RandGamma::shoot(neff,1.0)/neff;
178  }
179  //G4cout << "Gauss: " << loss << G4endl;
180  return loss;
181  }
182  }
183 
184  // Glandz regime : initialisation
185  //
186  if (material != lastMaterial) {
187  f1Fluct = material->GetIonisation()->GetF1fluct();
188  f2Fluct = material->GetIonisation()->GetF2fluct();
189  e1Fluct = material->GetIonisation()->GetEnergy1fluct();
190  e2Fluct = material->GetIonisation()->GetEnergy2fluct();
191  e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
192  e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
193  ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy();
194  ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
195  e0 = material->GetIonisation()->GetEnergy0fluct();
196  esmall = 0.5*sqrt(e0*ipotFluct);
197  lastMaterial = material;
198  }
199 
200  // very small step or low-density material
201  if(tmax <= e0) { return meanLoss; }
202 
203  G4double losstot = 0.;
204  G4int nstep = 1;
205  if(meanLoss < 25.*ipotFluct)
206  {
207  if(G4UniformRand() < 0.04*meanLoss/ipotFluct)
208  { nstep = 1; }
209  else
210  {
211  nstep = 2;
212  meanLoss *= 0.5;
213  }
214  }
215 
216  for (G4int istep=0; istep < nstep; ++istep) {
217 
218  loss = 0.;
219 
220  G4double a1 = 0. , a2 = 0., a3 = 0. ;
221 
222  if(tmax > ipotFluct) {
223  G4double w2 = log(2.*electron_mass_c2*beta2*gam2)-beta2;
224 
225  if(w2 > ipotLogFluct) {
226  G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
227  a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
228  if(w2 > e2LogFluct) {
229  a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
230  }
231  if(a1 < nmaxCont) {
232  //small energy loss
233  G4double sa1 = sqrt(a1);
234  if(G4UniformRand() < exp(-sa1))
235  {
236  e1 = esmall;
237  a1 = meanLoss*(1.-rate)/e1;
238  a2 = 0.;
239  e2 = e2Fluct;
240  }
241  else
242  {
243  a1 = sa1 ;
244  e1 = sa1*e1Fluct;
245  e2 = e2Fluct;
246  }
247 
248  } else {
249  //not small energy loss
250  //correction to get better fwhm value
251  a1 /= fw;
252  e1 = fw*e1Fluct;
253  e2 = e2Fluct;
254  }
255  }
256  }
257 
258  G4double w1 = tmax/e0;
259  if(tmax > e0) {
260  a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
261  }
262  //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont
263  G4double emean = 0.;
264  G4double sig2e = 0., sige = 0.;
265  G4double p1 = 0., p2 = 0., p3 = 0.;
266 
267  // excitation of type 1
268  if(a1 > nmaxCont)
269  {
270  emean += a1*e1;
271  sig2e += a1*e1*e1;
272  }
273  else if(a1 > 0.)
274  {
275  p1 = G4double(G4Poisson(a1));
276  loss += p1*e1;
277  if(p1 > 0.) {
278  loss += (1.-2.*G4UniformRand())*e1;
279  }
280  }
281 
282 
283  // excitation of type 2
284  if(a2 > nmaxCont)
285  {
286  emean += a2*e2;
287  sig2e += a2*e2*e2;
288  }
289  else if(a2 > 0.)
290  {
291  p2 = G4double(G4Poisson(a2));
292  loss += p2*e2;
293  if(p2 > 0.)
294  loss += (1.-2.*G4UniformRand())*e2;
295  }
296 
297  if(emean > 0.)
298  {
299  sige = sqrt(sig2e);
300  loss += max(0.,G4RandGauss::shoot(emean,sige));
301  }
302 
303  // ionisation
304  G4double lossc = 0.;
305  if(a3 > 0.) {
306  emean = 0.;
307  sig2e = 0.;
308  sige = 0.;
309  p3 = a3;
310  G4double alfa = 1.;
311  if(a3 > nmaxCont)
312  {
313  alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
314  G4double alfa1 = alfa*log(alfa)/(alfa-1.);
315  G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
316  emean += namean*e0*alfa1;
317  sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
318  p3 = a3-namean;
319  }
320 
321  G4double w2 = alfa*e0;
322  G4double w = (tmax-w2)/tmax;
323  G4int nb = G4Poisson(p3);
324  if(nb > 0) {
325  for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
326  }
327  }
328 
329  if(emean > 0.)
330  {
331  sige = sqrt(sig2e);
332  lossc += max(0.,G4RandGauss::shoot(emean,sige));
333  }
334 
335  loss += lossc;
336 
337  losstot += loss;
338  }
339  //G4cout << "Vavilov: " << losstot << " Nstep= " << nstep << G4endl;
340 
341  return losstot;
342 
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
346 
347 
348 G4double G4UniversalFluctuation::Dispersion(
349  const G4Material* material,
350  const G4DynamicParticle* dp,
351  G4double& tmax,
352  G4double& length)
353 {
354  if(!particle) { InitialiseMe(dp->GetDefinition()); }
355 
356  electronDensity = material->GetElectronDensity();
357 
358  G4double gam = (dp->GetKineticEnergy())/particleMass + 1.0;
359  G4double beta2 = 1.0 - 1.0/(gam*gam);
360 
361  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
362  * electronDensity * chargeSquare;
363 
364  return siga;
365 }
366 
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
368 
369 void
370 G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,
371  G4double q2)
372 {
373  if(part != particle) {
374  particle = part;
375  particleMass = part->GetPDGMass();
376  }
377  chargeSquare = q2;
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4