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

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

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