source: trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation93.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 13 years ago

update to last version 4.9.4
File size: 10.7 KB
Line 
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: G4UniversalFluctuation93.cc,v 1.2 2010/10/26 10:06:12 vnivanch Exp $
27// GEANT4 tag $Name: emstand-V09-03-24 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4UniversalFluctuation93
35//
36// Author:        V.Ivanchenko make a class with the Laszlo Urban model
37//
38// Creation date: 03.01.2002
39//
40// Modifications:
41//
42// 28-12-02 add method Dispersion (V.Ivanchenko)
43// 07-02-03 change signature (V.Ivanchenko)
44// 13-02-03 Add name (V.Ivanchenko)
45// 16-10-03 Changed interface to Initialisation (V.Ivanchenko)
46// 07-11-03 Fix problem of rounding of double
47// 06-02-04 Add control on big sigma > 2*meanLoss (V.Ivanchenko)
48// 26-04-04 Comment out the case of very small step (V.Ivanchenko)
49// 07-02-05 define problim = 5.e-3 (mma)
50// 03-05-05 conditions of Gaussian fluctuation changed (bugfix)
51//          + smearing for very small loss (L.Urban)
52// 03-10-05 energy dependent rate -> cut dependence of the
53//          distribution is much weaker (L.Urban)
54// 17-10-05 correction for very small loss (L.Urban)
55// 20-03-07 'GLANDZ' part rewritten completely, no 'very small loss'
56//          regime any more (L.Urban)
57// 03-04-07 correction to get better width of eloss distr.(L.Urban)
58// 13-07-07 add protection for very small step or low-density material (VI)
59// 19-03-09 new width correction (does not depend on previous steps) (L.Urban)
60// 20-03-09 modification in the width correction (L.Urban)
61// 14-06-10 saved version of 9.3 model with the name G4UniversalFluctuation93
62//
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66
67#include "G4UniversalFluctuation93.hh"
68#include "Randomize.hh"
69#include "G4Poisson.hh"
70#include "G4Step.hh"
71#include "G4Material.hh"
72#include "G4DynamicParticle.hh"
73#include "G4ParticleDefinition.hh"
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76
77using namespace std;
78
79G4UniversalFluctuation93::G4UniversalFluctuation93(const G4String& nam)
80 :G4VEmFluctuationModel(nam),
81  particle(0),
82  minNumberInteractionsBohr(10.0),
83  theBohrBeta2(50.0*keV/proton_mass_c2),
84  minLoss(10.*eV),
85  nmaxCont1(4.),
86  nmaxCont2(16.)
87{
88  lastMaterial = 0;
89
90  particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct
91    = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 
92    = e1 = e2 = 0;
93}
94
95//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96
97G4UniversalFluctuation93::~G4UniversalFluctuation93()
98{}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101
102void G4UniversalFluctuation93::InitialiseMe(const G4ParticleDefinition* part)
103{
104  particle       = part;
105  particleMass   = part->GetPDGMass();
106  G4double q     = part->GetPDGCharge()/eplus;
107  chargeSquare   = q*q;
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111
112G4double
113G4UniversalFluctuation93::SampleFluctuations(const G4Material* material,
114                                             const G4DynamicParticle* dp,
115                                             G4double& tmax,
116                                             G4double& length,
117                                             G4double& meanLoss)
118{
119// Calculate actual loss from the mean loss.
120// The model used to get the fluctuations is essentially the same
121// as in Glandz in Geant3 (Cern program library W5013, phys332).
122// L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
123
124  // shortcut for very very small loss (out of validity of the model)
125  //
126  if (meanLoss < minLoss)
127    return meanLoss;
128
129  if(!particle) InitialiseMe(dp->GetDefinition());
130
131  G4double tau   = dp->GetKineticEnergy()/particleMass;
132  G4double gam   = tau + 1.0;
133  G4double gam2  = gam*gam;
134  G4double beta2 = tau*(tau + 2.0)/gam2;
135
136  G4double loss(0.), siga(0.);
137 
138  // Gaussian regime
139  // for heavy particles only and conditions
140  // for Gauusian fluct. has been changed
141  //
142  if ((particleMass > electron_mass_c2) &&
143      (meanLoss >= minNumberInteractionsBohr*tmax))
144  {
145    G4double massrate = electron_mass_c2/particleMass ;
146    G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
147                        (1.+massrate*(2.*gam+massrate)) ;
148    if (tmaxkine <= 2.*tmax)   
149    {
150      electronDensity = material->GetElectronDensity();
151      siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
152                                * electronDensity * chargeSquare;
153      siga = sqrt(siga);
154      G4double twomeanLoss = meanLoss + meanLoss;
155      if (twomeanLoss < siga) {
156        G4double x;
157        do {
158          loss = twomeanLoss*G4UniformRand();
159          x = (loss - meanLoss)/siga;
160        } while (1.0 - 0.5*x*x < G4UniformRand());
161      } else {
162        do {
163          loss = G4RandGauss::shoot(meanLoss,siga);
164        } while (loss < 0. || loss > twomeanLoss);
165      }
166      return loss;
167    }
168  }
169
170  // Glandz regime : initialisation
171  //
172  if (material != lastMaterial) {
173    f1Fluct      = material->GetIonisation()->GetF1fluct();
174    f2Fluct      = material->GetIonisation()->GetF2fluct();
175    e1Fluct      = material->GetIonisation()->GetEnergy1fluct();
176    e2Fluct      = material->GetIonisation()->GetEnergy2fluct();
177    e1LogFluct   = material->GetIonisation()->GetLogEnergy1fluct();
178    e2LogFluct   = material->GetIonisation()->GetLogEnergy2fluct();
179    ipotFluct    = material->GetIonisation()->GetMeanExcitationEnergy();
180    ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
181    e0 = material->GetIonisation()->GetEnergy0fluct();
182    lastMaterial = material;
183 
184    // modification of some model parameters
185    // (this part should go to materials later)
186    G4double p = 1.40;
187    f2Fluct *= p;
188    f1Fluct = 1.-f2Fluct;
189    G4double q = 1.00;
190    e2Fluct *= q;
191    e2LogFluct = log(e2Fluct);
192    e1LogFluct = (ipotLogFluct-f2Fluct*e2LogFluct)/f1Fluct;
193    e1Fluct = exp(e1LogFluct);
194  }
195
196  // very small step or low-density material
197  if(tmax <= e0) return meanLoss;
198
199  G4double a1 = 0. , a2 = 0., a3 = 0. ;
200
201  // cut and material dependent rate
202  G4double rate = 1.0;
203  if(tmax > ipotFluct) {
204    G4double w2 = log(2.*electron_mass_c2*beta2*gam2)-beta2;
205
206    if(w2 > ipotLogFluct && w2 > e2LogFluct) {
207
208      rate = 0.03+0.23*log(log(tmax/ipotFluct));
209      G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
210      a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
211      a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
212      // correction in order to get better FWHM values
213      // ( scale parameters a1 and e1)
214      G4double width = 1.;
215      if(meanLoss > 10.*e1Fluct)
216      {
217        width = 3.1623/sqrt(meanLoss/e1Fluct);
218        if(width < a2/a1)
219        width = a2/a1;
220      } 
221      a1 *= width;
222      e1 = e1Fluct/width;
223      e2 = e2Fluct;
224    }
225  }
226
227  G4double w1 = tmax/e0;
228  if(tmax > e0) 
229    a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
230
231  //'nearly' Gaussian fluctuation if a1>nmaxCont2&&a2>nmaxCont2&&a3>nmaxCont2 
232  G4double emean = 0.;
233  G4double sig2e = 0., sige = 0.;
234  G4double p1 = 0., p2 = 0., p3 = 0.;
235 
236  // excitation of type 1
237  if(a1 > nmaxCont2)
238  {
239    emean += a1*e1;
240    sig2e += a1*e1*e1;
241  }
242  else if(a1 > 0.)
243  {
244    p1 = G4double(G4Poisson(a1));
245    loss += p1*e1;
246    if(p1 > 0.) 
247      loss += (1.-2.*G4UniformRand())*e1;
248  }
249
250  // excitation of type 2
251  if(a2 > nmaxCont2)
252  {
253    emean += a2*e2;
254    sig2e += a2*e2*e2;
255  }
256  else if(a2 > 0.)
257  {
258    p2 = G4double(G4Poisson(a2));
259    loss += p2*e2;
260    if(p2 > 0.) 
261      loss += (1.-2.*G4UniformRand())*e2;
262  }
263
264  // ionisation
265  G4double lossc = 0.;
266  if(a3 > 0.)
267  {
268    p3 = a3;
269    G4double alfa = 1.;
270    if(a3 > nmaxCont2)
271    {
272       alfa            = w1*(nmaxCont2+a3)/(w1*nmaxCont2+a3);
273       G4double alfa1  = alfa*log(alfa)/(alfa-1.);
274       G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
275       emean          += namean*e0*alfa1;
276       sig2e          += e0*e0*namean*(alfa-alfa1*alfa1);
277       p3              = a3-namean;
278    }
279
280    G4double w2 = alfa*e0;
281    G4double w  = (tmax-w2)/tmax;
282    G4int nb = G4Poisson(p3);
283    if(nb > 0)
284      for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
285  }
286
287  if(emean > 0.)
288  {
289    sige   = sqrt(sig2e);
290    loss += max(0.,G4RandGauss::shoot(emean,sige));
291  }
292
293  loss += lossc;
294
295  return loss;
296
297}
298
299//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
300
301
302G4double G4UniversalFluctuation93::Dispersion(
303                          const G4Material* material,
304                          const G4DynamicParticle* dp,
305                                G4double& tmax,
306                                G4double& length)
307{
308  if(!particle) InitialiseMe(dp->GetDefinition());
309
310  electronDensity = material->GetElectronDensity();
311
312  G4double gam   = (dp->GetKineticEnergy())/particleMass + 1.0;
313  G4double beta2 = 1.0 - 1.0/(gam*gam);
314
315  G4double siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
316                 * electronDensity * chargeSquare;
317
318  return siga;
319}
320
321//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
322
323void G4UniversalFluctuation93::SetParticleAndCharge(
324          const G4ParticleDefinition* part, G4double q2)
325{
326  if(part != particle) {
327    particle       = part;
328    particleMass   = part->GetPDGMass();
329  }
330  chargeSquare = q2;
331}
332
333//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.