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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

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