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