source: trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc@ 1005

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

fichier oublies
File size: 14.4 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: G4IonFluctuations.cc,v 1.25 2009/02/19 19:17:50 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4IonFluctuation
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// 23-05-03 Add control on parthalogical cases (V.Ivanchenko)
46// 16-10-03 Changed interface to Initialisation (V.Ivanchenko)
47// 27-09-07 Use FermiEnergy from material, add cut dependence (V.Ivanchenko)
48// 01-02-08 Add protection for small energies and optimise the code (V.Ivanchenko)
49// 01-06-08 Added initialisation of effective charge prestep (V.Ivanchenko)
50//
51// Class Description:
52//
53// -------------------------------------------------------------------
54//
55
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58
59#include "G4IonFluctuations.hh"
60#include "Randomize.hh"
61#include "G4Poisson.hh"
62#include "G4Material.hh"
63#include "G4DynamicParticle.hh"
64
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
67using namespace std;
68
69G4IonFluctuations::G4IonFluctuations(const G4String& nam)
70 : G4VEmFluctuationModel(nam),
71 particle(0),
72 particleMass(proton_mass_c2),
73 charge(1.0),
74 chargeSquare(1.0),
75 effChargeSquare(1.0),
76 parameter(10.0*CLHEP::MeV/CLHEP::proton_mass_c2),
77 minNumberInteractionsBohr(0.0),
78 theBohrBeta2(50.0*keV/CLHEP::proton_mass_c2),
79 minFraction(0.2),
80 xmin(0.2),
81 minLoss(0.001*eV)
82{}
83
84//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85
86G4IonFluctuations::~G4IonFluctuations()
87{}
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91void G4IonFluctuations::InitialiseMe(const G4ParticleDefinition* part)
92{
93 particle = part;
94 particleMass = part->GetPDGMass();
95 charge = part->GetPDGCharge()/eplus;
96 chargeSquare = charge*charge;
97 effChargeSquare= chargeSquare;
98 uniFluct.InitialiseMe(part);
99}
100
101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102
103G4double G4IonFluctuations::SampleFluctuations(const G4Material* material,
104 const G4DynamicParticle* dp,
105 G4double& tmax,
106 G4double& length,
107 G4double& meanLoss)
108{
109 // G4cout << "### meanLoss= " << meanLoss << G4endl;
110 if(meanLoss <= minLoss) return meanLoss;
111
112 //G4cout << "G4IonFluctuations::SampleFluctuations E(MeV)= " << dp->GetKineticEnergy()
113 // << " Elim(MeV)= " << parameter*charge*particleMass << G4endl;
114
115 // Vavilov fluctuations
116 if(dp->GetKineticEnergy() > parameter*charge*particleMass) {
117 return uniFluct.SampleFluctuations(material,dp,tmax,length,meanLoss);
118 }
119
120 G4double siga = Dispersion(material,dp,tmax,length);
121 G4double loss = meanLoss;
122
123 G4double navr = minNumberInteractionsBohr;
124 navr = meanLoss*meanLoss/siga;
125 //G4cout << "### siga= " << sqrt(siga) << " navr= " << navr << G4endl;
126
127 // Gaussian fluctuation
128 if (navr >= minNumberInteractionsBohr) {
129
130 // Increase fluctuations for big fractional energy loss
131 //G4cout << "siga= " << siga << G4endl;
132 if ( meanLoss > minFraction*kineticEnergy ) {
133 G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0;
134 G4double b2 = 1.0 - 1.0/(gam*gam);
135 if(b2 < xmin*beta2) b2 = xmin*beta2;
136 G4double x = b2/beta2;
137 G4double x3 = 1.0/(x*x*x);
138 siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) );
139 }
140 // G4cout << "siga= " << siga << G4endl;
141 siga = sqrt(siga);
142 G4double lossmax = meanLoss+meanLoss;
143
144 if(siga > 5.0*meanLoss) {
145 loss = lossmax*G4UniformRand();
146 } else {
147 do {
148 loss = G4RandGauss::shoot(meanLoss,siga);
149 } while (0.0 > loss || loss > lossmax);
150 }
151 // Poisson fluctuations
152 } else {
153
154 G4double n = (G4double)(G4Poisson(navr));
155 loss = meanLoss*n/navr;
156 }
157
158 //G4cout << "meanLoss= " << meanLoss << " loss= " << loss << G4endl;
159 return loss;
160}
161
162//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163
164G4double G4IonFluctuations::Dispersion(const G4Material* material,
165 const G4DynamicParticle* dp,
166 G4double& tmax,
167 G4double& length)
168{
169 kineticEnergy = dp->GetKineticEnergy();
170 G4double etot = kineticEnergy + particleMass;
171 beta2 = kineticEnergy*(kineticEnergy + 2.*particleMass)/(etot*etot);
172
173 G4double electronDensity = material->GetElectronDensity();
174
175 /*
176 G4cout << "e= " << kineticEnergy << " m= " << particleMass
177 << " tmax= " << tmax << " l= " << length
178 << " q^2= " << effChargeSquare << " beta2=" << beta2<< G4endl;
179 */
180 G4double siga = (1. - beta2*0.5)*tmax*length*electronDensity*
181 twopi_mc2_rcl2*chargeSquare/beta2;
182
183 // Low velocity - additional ion charge fluctuations according to
184 // Q.Yang et al., NIM B61(1991)149-155.
185 //G4cout << "sigE= " << sqrt(siga) << " charge= " << charge <<G4endl;
186
187 G4double Z = electronDensity/material->GetTotNbOfAtomsPerVolume();
188
189 G4double fac = Factor(material, Z);
190
191 // heavy ion correction
192 G4double f1 = 1.065e-4*chargeSquare;
193 if(beta2 > theBohrBeta2) f1/= beta2;
194 else f1/= theBohrBeta2;
195 if(f1 > 2.5) f1 = 2.5;
196 fac *= (1.0 + f1);
197
198 // taking into account the cut
199 if(fac > 1.0) {
200 siga *= (1.0 + (fac - 1.0)*2.0*electron_mass_c2*beta2/(tmax*(1.0 - beta2)));
201 }
202 //G4cout << "siga(keV)= " << sqrt(siga)/keV << " fac= " << fac
203 // << " f1= " << f1 << G4endl;
204
205 return siga;
206}
207
208//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209
210G4double G4IonFluctuations::Factor(const G4Material* material, G4double Z)
211{
212 // The aproximation of energy loss fluctuations
213 // Q.Yang et al., NIM B61(1991)149-155.
214
215 // Reduced energy in MeV/AMU
216 G4double energy = kineticEnergy *amu_c2/(particleMass*MeV) ;
217
218 // simple approximation for higher beta2
219 G4double s1 = RelativisticFactor(material, Z);
220
221 // tabulation for lower beta2
222 if( beta2 < 3.0*theBohrBeta2*Z ) {
223
224 static G4double a[96][4] = {
225 {-0.3291, -0.8312, 0.2460, -1.0220},
226 {-0.5615, -0.5898, 0.5205, -0.7258},
227 {-0.5280, -0.4981, 0.5519, -0.5865},
228 {-0.5125, -0.4625, 0.5660, -0.5190},
229 {-0.5127, -0.8595, 0.5626, -0.8721},
230 {-0.5174, -1.1930, 0.5565, -1.1980},
231 {-0.5179, -1.1850, 0.5560, -1.2070},
232 {-0.5209, -0.9355, 0.5590, -1.0250},
233 {-0.5255, -0.7766, 0.5720, -0.9412},
234
235 {-0.5776, -0.6665, 0.6598, -0.8484},
236 {-0.6013, -0.6045, 0.7321, -0.7671},
237 {-0.5781, -0.5518, 0.7605, -0.6919},
238 {-0.5587, -0.4981, 0.7835, -0.6195},
239 {-0.5466, -0.4656, 0.7978, -0.5771},
240 {-0.5406, -0.4690, 0.8031, -0.5718},
241 {-0.5391, -0.5061, 0.8024, -0.5974},
242 {-0.5380, -0.6483, 0.7962, -0.6970},
243 {-0.5355, -0.7722, 0.7962, -0.7839},
244 {-0.5329, -0.7720, 0.7988, -0.7846},
245
246 {-0.5335, -0.7671, 0.7984, -0.7933},
247 {-0.5324, -0.7612, 0.7998, -0.8031},
248 {-0.5305, -0.7300, 0.8031, -0.7990},
249 {-0.5307, -0.7178, 0.8049, -0.8216},
250 {-0.5248, -0.6621, 0.8165, -0.7919},
251 {-0.5180, -0.6502, 0.8266, -0.7986},
252 {-0.5084, -0.6408, 0.8396, -0.8048},
253 {-0.4967, -0.6331, 0.8549, -0.8093},
254 {-0.4861, -0.6508, 0.8712, -0.8432},
255 {-0.4700, -0.6186, 0.8961, -0.8132},
256
257 {-0.4545, -0.5720, 0.9227, -0.7710},
258 {-0.4404, -0.5226, 0.9481, -0.7254},
259 {-0.4288, -0.4778, 0.9701, -0.6850},
260 {-0.4199, -0.4425, 0.9874, -0.6539},
261 {-0.4131, -0.4188, 0.9998, -0.6332},
262 {-0.4089, -0.4057, 1.0070, -0.6218},
263 {-0.4039, -0.3913, 1.0150, -0.6107},
264 {-0.3987, -0.3698, 1.0240, -0.5938},
265 {-0.3977, -0.3608, 1.0260, -0.5852},
266 {-0.3972, -0.3600, 1.0260, -0.5842},
267
268 {-0.3985, -0.3803, 1.0200, -0.6013},
269 {-0.3985, -0.3979, 1.0150, -0.6168},
270 {-0.3968, -0.3990, 1.0160, -0.6195},
271 {-0.3971, -0.4432, 1.0050, -0.6591},
272 {-0.3944, -0.4665, 1.0010, -0.6825},
273 {-0.3924, -0.5109, 0.9921, -0.7235},
274 {-0.3882, -0.5158, 0.9947, -0.7343},
275 {-0.3838, -0.5125, 0.9999, -0.7370},
276 {-0.3786, -0.4976, 1.0090, -0.7310},
277 {-0.3741, -0.4738, 1.0200, -0.7155},
278
279 {-0.3969, -0.4496, 1.0320, -0.6982},
280 {-0.3663, -0.4297, 1.0430, -0.6828},
281 {-0.3630, -0.4120, 1.0530, -0.6689},
282 {-0.3597, -0.3964, 1.0620, -0.6564},
283 {-0.3555, -0.3809, 1.0720, -0.6454},
284 {-0.3525, -0.3607, 1.0820, -0.6289},
285 {-0.3505, -0.3465, 1.0900, -0.6171},
286 {-0.3397, -0.3570, 1.1020, -0.6384},
287 {-0.3314, -0.3552, 1.1130, -0.6441},
288 {-0.3235, -0.3531, 1.1230, -0.6498},
289
290 {-0.3150, -0.3483, 1.1360, -0.6539},
291 {-0.3060, -0.3441, 1.1490, -0.6593},
292 {-0.2968, -0.3396, 1.1630, -0.6649},
293 {-0.2935, -0.3225, 1.1760, -0.6527},
294 {-0.2797, -0.3262, 1.1940, -0.6722},
295 {-0.2704, -0.3202, 1.2100, -0.6770},
296 {-0.2815, -0.3227, 1.2480, -0.6775},
297 {-0.2880, -0.3245, 1.2810, -0.6801},
298 {-0.3034, -0.3263, 1.3270, -0.6778},
299 {-0.2936, -0.3215, 1.3430, -0.6835},
300
301 {-0.3282, -0.3200, 1.3980, -0.6650},
302 {-0.3260, -0.3070, 1.4090, -0.6552},
303 {-0.3511, -0.3074, 1.4470, -0.6442},
304 {-0.3501, -0.3064, 1.4500, -0.6442},
305 {-0.3490, -0.3027, 1.4550, -0.6418},
306 {-0.3487, -0.3048, 1.4570, -0.6447},
307 {-0.3478, -0.3074, 1.4600, -0.6483},
308 {-0.3501, -0.3283, 1.4540, -0.6669},
309 {-0.3494, -0.3373, 1.4550, -0.6765},
310 {-0.3485, -0.3373, 1.4570, -0.6774},
311
312 {-0.3462, -0.3300, 1.4630, -0.6728},
313 {-0.3462, -0.3225, 1.4690, -0.6662},
314 {-0.3453, -0.3094, 1.4790, -0.6553},
315 {-0.3844, -0.3134, 1.5240, -0.6412},
316 {-0.3848, -0.3018, 1.5310, -0.6303},
317 {-0.3862, -0.2955, 1.5360, -0.6237},
318 {-0.4262, -0.2991, 1.5860, -0.6115},
319 {-0.4278, -0.2910, 1.5900, -0.6029},
320 {-0.4303, -0.2817, 1.5940, -0.5927},
321 {-0.4315, -0.2719, 1.6010, -0.5829},
322
323 {-0.4359, -0.2914, 1.6050, -0.6010},
324 {-0.4365, -0.2982, 1.6080, -0.6080},
325 {-0.4253, -0.3037, 1.6120, -0.6150},
326 {-0.4335, -0.3245, 1.6160, -0.6377},
327 {-0.4307, -0.3292, 1.6210, -0.6447},
328 {-0.4284, -0.3204, 1.6290, -0.6380},
329 {-0.4227, -0.3217, 1.6360, -0.6438}
330 } ;
331
332 G4int iz = G4int(Z) - 2;
333 if( 0 > iz ) iz = 0;
334 else if(95 < iz ) iz = 95;
335
336 G4double ss = 1.0 + a[iz][0]*pow(energy,a[iz][1])+
337 + a[iz][2]*pow(energy,a[iz][3]);
338
339 // protection for the validity range for low beta
340 G4double slim = 0.001;
341 if(ss < slim) s1 = 1.0/slim;
342 // for high value of beta
343 else if(s1*ss < 1.0) s1 = 1.0/ss;
344 }
345
346 G4int i = 0 ;
347 G4double factor = 1.0 ;
348
349 // The index of set of parameters i = 0 for protons(hadrons) in gases
350 // 1 for protons(hadrons) in solids
351 // 2 for ions in atomic gases
352 // 3 for ions in molecular gases
353 // 4 for ions in solids
354 static G4double b[5][4] = {
355 {0.1014, 0.3700, 0.9642, 3.987},
356 {0.1955, 0.6941, 2.522, 1.040},
357 {0.05058, 0.08975, 0.1419, 10.80},
358 {0.05009, 0.08660, 0.2751, 3.787},
359 {0.01273, 0.03458, 0.3951, 3.812}
360 } ;
361
362 // protons (hadrons)
363 if(1.5 > charge) {
364 if( kStateGas != material->GetState() ) i = 1 ;
365
366 // ions
367 } else {
368
369 factor = charge * pow(charge/Z, 0.33333333);
370
371 if( kStateGas == material->GetState() ) {
372 energy /= (charge * sqrt(charge)) ;
373
374 if(1 == (material->GetNumberOfElements())) {
375 i = 2 ;
376 } else {
377 i = 3 ;
378 }
379
380 } else {
381 energy /= (charge * sqrt(charge*Z)) ;
382 i = 4 ;
383 }
384 }
385
386 G4double x = b[i][2];
387 G4double y = energy * b[i][3];
388 if(y <= 0.2) x *= (y*(1.0 - 0.5*y));
389 else x *= (1.0 - exp(-y));
390
391 y = energy - b[i][1];
392
393 G4double s2 = factor * x * b[i][0] / (y*y + x*x);
394 /*
395 G4cout << "s1= " << s1 << " s2= " << s2 << " q^2= " << effChargeSquare
396 << " e= " << energy << G4endl;
397 */
398 return s1*effChargeSquare/chargeSquare + s2;
399}
400
401//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
402
403G4double G4IonFluctuations::RelativisticFactor(const G4Material* mat,
404 G4double Z)
405{
406 G4double eF = mat->GetIonisation()->GetFermiEnergy();
407 G4double I = mat->GetIonisation()->GetMeanExcitationEnergy();
408
409 // H.Geissel et al. NIM B, 195 (2002) 3.
410 G4double bF2= 2.0*eF/electron_mass_c2;
411 G4double f = 0.4*(1.0 - beta2)/((1.0 - 0.5*beta2)*Z);
412 if(beta2 > bF2) f *= log(2.0*electron_mass_c2*beta2/I)*bF2/beta2;
413 else f *= log(4.0*eF/I);
414
415 // G4cout << "f= " << f << " beta2= " << beta2
416 // << " bf2= " << bF2 << " q^2= " << chargeSquare << G4endl;
417
418 return 1.0 + f;
419}
420
421//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
422
423void G4IonFluctuations::SetParticleAndCharge(const G4ParticleDefinition* part,
424 G4double q2)
425{
426 if(part != particle) {
427 particle = part;
428 particleMass = part->GetPDGMass();
429 charge = part->GetPDGCharge()/eplus;
430 chargeSquare = charge*charge;
431 }
432 effChargeSquare = q2;
433 uniFluct.SetParticleAndCharge(part, q2);
434}
435
436//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.