source: trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 9.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: G4UniversalFluctuation.cc,v 1.16 2008/10/22 16:04:33 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4UniversalFluctuation
35//
36// Author: Vladimir Ivanchenko
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 in G4UniversalFluctuations
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//
60
61//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63
64#include "G4UniversalFluctuation.hh"
65#include "Randomize.hh"
66#include "G4Poisson.hh"
67#include "G4Step.hh"
68#include "G4Material.hh"
69#include "G4DynamicParticle.hh"
70#include "G4ParticleDefinition.hh"
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74using namespace std;
75
76G4UniversalFluctuation::G4UniversalFluctuation(const G4String& nam)
77 :G4VEmFluctuationModel(nam),
78 particle(0),
79 minNumberInteractionsBohr(10.0),
80 theBohrBeta2(50.0*keV/proton_mass_c2),
81 minLoss(10.*eV),
82 nmaxCont1(4.),
83 nmaxCont2(16.)
84{
85 lastMaterial = 0;
86 facwidth = 1.000/keV;
87 oldloss = 0.;
88 samestep = 0.;
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92
93G4UniversalFluctuation::~G4UniversalFluctuation()
94{}
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97
98void G4UniversalFluctuation::InitialiseMe(const G4ParticleDefinition* part)
99{
100 particle = part;
101 particleMass = part->GetPDGMass();
102 G4double q = part->GetPDGCharge()/eplus;
103 chargeSquare = q*q;
104}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108G4double G4UniversalFluctuation::SampleFluctuations(const G4Material* material,
109 const G4DynamicParticle* dp,
110 G4double& tmax,
111 G4double& length,
112 G4double& meanLoss)
113{
114// Calculate actual loss from the mean loss.
115// The model used to get the fluctuations is essentially the same
116// as in Glandz in Geant3 (Cern program library W5013, phys332).
117// L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
118
119 // shortcut for very very small loss (out of validity of the model)
120 //
121 if (meanLoss < minLoss)
122 {
123 oldloss = meanLoss;
124 return meanLoss;
125 }
126
127 if(!particle) InitialiseMe(dp->GetDefinition());
128
129 G4double tau = dp->GetKineticEnergy()/particleMass;
130 G4double gam = tau + 1.0;
131 G4double gam2 = gam*gam;
132 G4double beta2 = tau*(tau + 2.0)/gam2;
133
134 G4double loss(0.), siga(0.);
135
136 // Gaussian regime
137 // for heavy particles only and conditions
138 // for Gauusian fluct. has been changed
139 //
140 if ((particleMass > electron_mass_c2) &&
141 (meanLoss >= minNumberInteractionsBohr*tmax))
142 {
143 G4double massrate = electron_mass_c2/particleMass ;
144 G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
145 (1.+massrate*(2.*gam+massrate)) ;
146 if (tmaxkine <= 2.*tmax)
147 {
148 electronDensity = material->GetElectronDensity();
149 siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
150 * electronDensity * chargeSquare;
151 siga = sqrt(siga);
152 G4double twomeanLoss = meanLoss + meanLoss;
153 if (twomeanLoss < siga) {
154 G4double x;
155 do {
156 loss = twomeanLoss*G4UniformRand();
157 x = (loss - meanLoss)/siga;
158 } while (1.0 - 0.5*x*x < G4UniformRand());
159 } else {
160 do {
161 loss = G4RandGauss::shoot(meanLoss,siga);
162 } while (loss < 0. || loss > twomeanLoss);
163 }
164 return loss;
165 }
166 }
167
168 // Glandz regime : initialisation
169 //
170 if (material != lastMaterial) {
171 f1Fluct = material->GetIonisation()->GetF1fluct();
172 f2Fluct = material->GetIonisation()->GetF2fluct();
173 e1Fluct = material->GetIonisation()->GetEnergy1fluct();
174 e2Fluct = material->GetIonisation()->GetEnergy2fluct();
175 e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
176 e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
177 ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy();
178 ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
179 e0 = material->GetIonisation()->GetEnergy0fluct();
180 lastMaterial = material;
181 }
182
183 // very small step or low-density material
184 if(tmax <= e0) return meanLoss;
185
186 G4double a1 = 0. , a2 = 0., a3 = 0. ;
187
188 // correction to get better width even using stepmax
189 if(abs(meanLoss- oldloss) < 1.*eV)
190 samestep += 1;
191 else
192 samestep = 1.;
193 oldloss = meanLoss;
194 G4double width = 1.+samestep*facwidth*meanLoss;
195 if(width > 4.50) width = 4.50;
196 e1 = width*e1Fluct;
197 e2 = width*e2Fluct;
198
199 // cut and material dependent rate
200 G4double rate = 1.0;
201 if(tmax > ipotFluct) {
202 G4double w2 = log(2.*electron_mass_c2*beta2*gam2)-beta2;
203
204 if(w2 > ipotLogFluct && w2 > e2LogFluct) {
205
206 rate = 0.03+0.23*log(log(tmax/ipotFluct));
207 G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
208 a1 = C*f1Fluct*(w2-e1LogFluct)/e1;
209 a2 = C*f2Fluct*(w2-e2LogFluct)/e2;
210 }
211 }
212
213 G4double w1 = tmax/e0;
214 if(tmax > e0)
215 a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
216
217 //'nearly' Gaussian fluctuation if a1>nmaxCont2&&a2>nmaxCont2&&a3>nmaxCont2
218 G4double emean = 0.;
219 G4double sig2e = 0., sige = 0.;
220 G4double p1 = 0., p2 = 0., p3 = 0.;
221
222 // excitation of type 1
223 if(a1 > nmaxCont2)
224 {
225 emean += a1*e1;
226 sig2e += a1*e1*e1;
227 }
228 else if(a1 > 0.)
229 {
230 p1 = G4double(G4Poisson(a1));
231 loss += p1*e1;
232 if(p1 > 0.)
233 loss += (1.-2.*G4UniformRand())*e1;
234 }
235
236 // excitation of type 2
237 if(a2 > nmaxCont2)
238 {
239 emean += a2*e2;
240 sig2e += a2*e2*e2;
241 }
242 else if(a2 > 0.)
243 {
244 p2 = G4double(G4Poisson(a2));
245 loss += p2*e2;
246 if(p2 > 0.)
247 loss += (1.-2.*G4UniformRand())*e2;
248 }
249
250 // ionisation
251 G4double lossc = 0.;
252 if(a3 > 0.)
253 {
254 p3 = a3;
255 G4double alfa = 1.;
256 if(a3 > nmaxCont2)
257 {
258 alfa = w1*(nmaxCont2+a3)/(w1*nmaxCont2+a3);
259 G4double alfa1 = alfa*log(alfa)/(alfa-1.);
260 G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
261 emean += namean*e0*alfa1;
262 sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
263 p3 = a3-namean;
264 }
265
266 G4double w2 = alfa*e0;
267 G4double w = (tmax-w2)/tmax;
268 G4int nb = G4Poisson(p3);
269 if(nb > 0)
270 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
271 }
272
273 if(emean > 0.)
274 {
275 sige = sqrt(sig2e);
276 loss += max(0.,G4RandGauss::shoot(emean,sige));
277 }
278
279 loss += lossc;
280
281 return loss;
282
283}
284
285//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
286
287
288G4double G4UniversalFluctuation::Dispersion(
289 const G4Material* material,
290 const G4DynamicParticle* dp,
291 G4double& tmax,
292 G4double& length)
293{
294 if(!particle) InitialiseMe(dp->GetDefinition());
295
296 electronDensity = material->GetElectronDensity();
297
298 G4double gam = (dp->GetKineticEnergy())/particleMass + 1.0;
299 G4double beta2 = 1.0 - 1.0/(gam*gam);
300
301 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
302 * electronDensity * chargeSquare;
303
304 return siga;
305}
306
307//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.