| 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: G4UniversalFluctuation.cc,v 1.16 2008/10/22 16:04:33 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-01 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class file 31 // GEANT4 Class file
30 // 32 //
31 // 33 //
32 // File name: G4UniversalFluctuation 34 // File name: G4UniversalFluctuation
33 // 35 //
34 // Author: V. Ivanchenko for Laszlo Urb << 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 // << 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 //
41 60
42 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 63
45 #include "G4UniversalFluctuation.hh" 64 #include "G4UniversalFluctuation.hh"
46 #include "G4PhysicalConstants.hh" <<
47 #include "G4SystemOfUnits.hh" <<
48 #include "Randomize.hh" 65 #include "Randomize.hh"
49 #include "G4Poisson.hh" 66 #include "G4Poisson.hh"
>> 67 #include "G4Step.hh"
50 #include "G4Material.hh" 68 #include "G4Material.hh"
51 #include "G4MaterialCutsCouple.hh" <<
52 #include "G4DynamicParticle.hh" 69 #include "G4DynamicParticle.hh"
53 #include "G4ParticleDefinition.hh" 70 #include "G4ParticleDefinition.hh"
54 #include "G4Log.hh" <<
55 71
56 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 73
>> 74 using namespace std;
>> 75
58 G4UniversalFluctuation::G4UniversalFluctuation 76 G4UniversalFluctuation::G4UniversalFluctuation(const G4String& nam)
59 :G4VEmFluctuationModel(nam), 77 :G4VEmFluctuationModel(nam),
60 minLoss(10.*CLHEP::eV) << 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.)
61 { 84 {
62 rndmarray = new G4double[sizearray]; << 85 lastMaterial = 0;
>> 86 facwidth = 1.000/keV;
>> 87 oldloss = 0.;
>> 88 samestep = 0.;
63 } 89 }
64 90
65 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 92
67 G4UniversalFluctuation::~G4UniversalFluctuatio 93 G4UniversalFluctuation::~G4UniversalFluctuation()
68 { << 94 {}
69 delete [] rndmarray; <<
70 } <<
71 95
72 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 97
74 void G4UniversalFluctuation::InitialiseMe(cons 98 void G4UniversalFluctuation::InitialiseMe(const G4ParticleDefinition* part)
75 { 99 {
76 particle = part; << 100 particle = part;
77 particleMass = part->GetPDGMass(); << 101 particleMass = part->GetPDGMass();
78 const G4double q = part->GetPDGCharge()/CLHE << 102 G4double q = part->GetPDGCharge()/eplus;
79 << 103 chargeSquare = q*q;
80 // Derived quantities <<
81 m_Inv_particleMass = 1.0 / particleMass; <<
82 m_massrate = CLHEP::electron_mass_c2 * m_Inv <<
83 chargeSquare = q*q; <<
84 } 104 }
85 105
86 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 107
88 G4double << 108 G4double G4UniversalFluctuation::SampleFluctuations(const G4Material* material,
89 G4UniversalFluctuation::SampleFluctuations(con << 109 const G4DynamicParticle* dp,
90 con << 110 G4double& tmax,
91 con << 111 G4double& length,
92 con << 112 G4double& meanLoss)
93 con <<
94 con <<
95 { 113 {
96 // Calculate actual loss from the mean loss. << 114 // Calculate actual loss from the mean loss.
97 // The model used to get the fluctuations is << 115 // The model used to get the fluctuations is essentially the same
98 // as in Glandz in Geant3 (Cern program libr << 116 // as in Glandz in Geant3 (Cern program library W5013, phys332).
99 // L. Urban et al. NIM A362, p.416 (1995) an << 117 // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
100 118
101 // shortcut for very small loss or from a st << 119 // shortcut for very very small loss (out of validity of the model)
102 // (out of validity of the model) <<
103 // 120 //
104 if (averageLoss < minLoss) { return averageL << 121 if (meanLoss < minLoss)
105 meanLoss = averageLoss; << 122 {
106 const G4double tkin = dp->GetKineticEnergy( << 123 oldloss = meanLoss;
107 //G4cout<< "Emean= "<< meanLoss<< " tmax= "< << 124 return meanLoss;
108 << 125 }
109 if(dp->GetDefinition() != particle) { Initia <<
110 <<
111 CLHEP::HepRandomEngine* rndmEngineF = G4Rand <<
112 <<
113 const G4double gam = tkin * m_Inv_particle <<
114 const G4double gam2 = gam*gam; <<
115 const G4double beta = dp->GetBeta(); <<
116 const G4double beta2 = beta*beta; <<
117 126
118 G4double loss(0.), siga(0.); << 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;
119 133
120 const G4Material* material = couple->GetMate << 134 G4double loss(0.), siga(0.);
121 135
122 // Gaussian regime 136 // Gaussian regime
123 // for heavy particles only and conditions 137 // for heavy particles only and conditions
124 // for Gauusian fluct. has been changed 138 // for Gauusian fluct. has been changed
125 // 139 //
126 if (particleMass > CLHEP::electron_mass_c2 & << 140 if ((particleMass > electron_mass_c2) &&
127 meanLoss >= minNumberInteractionsBohr*tc << 141 (meanLoss >= minNumberInteractionsBohr*tmax))
128 << 142 {
129 siga = std::sqrt((tmax/beta2 - 0.5*tcut)*C << 143 G4double massrate = electron_mass_c2/particleMass ;
130 length*chargeSquare*mate << 144 G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
131 const G4double sn = meanLoss/siga; << 145 (1.+massrate*(2.*gam+massrate)) ;
132 << 146 if (tmaxkine <= 2.*tmax)
133 // thick target case << 147 {
134 if (sn >= 2.0) { << 148 electronDensity = material->GetElectronDensity();
135 << 149 siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
136 const G4double twomeanLoss = meanLoss + << 150 * electronDensity * chargeSquare;
137 do { << 151 siga = sqrt(siga);
138 loss = G4RandGauss::shoot(rndmEngineF, meanL << 152 G4double twomeanLoss = meanLoss + meanLoss;
139 // Loop checking, 03-Aug-2015, Vladimir Ivan << 153 if (twomeanLoss < siga) {
140 } while (0.0 > loss || twomeanLoss < lo << 154 G4double x;
141 << 155 do {
142 // Gamma distribution << 156 loss = twomeanLoss*G4UniformRand();
143 } else { << 157 x = (loss - meanLoss)/siga;
144 << 158 } while (1.0 - 0.5*x*x < G4UniformRand());
145 const G4double neff = sn*sn; << 159 } else {
146 loss = meanLoss*G4RandGamma::shoot(rndmE << 160 do {
>> 161 loss = G4RandGauss::shoot(meanLoss,siga);
>> 162 } while (loss < 0. || loss > twomeanLoss);
>> 163 }
>> 164 return loss;
147 } 165 }
148 //G4cout << "Gauss: " << loss << G4endl; <<
149 return loss; <<
150 } 166 }
151 167
152 auto ioni = material->GetIonisation(); << 168 // Glandz regime : initialisation
153 e0 = ioni->GetEnergy0fluct(); << 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 }
154 182
155 // very small step or low-density material 183 // very small step or low-density material
156 if(tcut <= e0) { return meanLoss; } << 184 if(tmax <= e0) return meanLoss;
157 185
158 ipotFluct = ioni->GetMeanExcitationEnergy(); << 186 G4double a1 = 0. , a2 = 0., a3 = 0. ;
159 ipotLogFluct = ioni->GetLogMeanExcEnergy(); <<
160 187
161 // width correction for small cuts << 188 // correction to get better width even using stepmax
162 const G4double scaling = std::min(1.+0.5*CLH << 189 if(abs(meanLoss- oldloss) < 1.*eV)
163 meanLoss /= scaling; << 190 samestep += 1;
164 << 191 else
165 w2 = (tcut > ipotFluct) ? << 192 samestep = 1.;
166 G4Log(2.*CLHEP::electron_mass_c2*beta2*gam << 193 oldloss = meanLoss;
167 return SampleGlandz(rndmEngineF, material, t << 194 G4double width = 1.+samestep*facwidth*meanLoss;
168 } << 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 }
169 212
170 //....oooOO0OOooo........oooOO0OOooo........oo << 213 G4double w1 = tmax/e0;
>> 214 if(tmax > e0)
>> 215 a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
171 216
172 G4double << 217 //'nearly' Gaussian fluctuation if a1>nmaxCont2&&a2>nmaxCont2&&a3>nmaxCont2
173 G4UniversalFluctuation::SampleGlandz(CLHEP::He <<
174 const G4M <<
175 const G4d <<
176 { <<
177 G4double a1(0.0), a3(0.0); <<
178 G4double loss = 0.0; <<
179 G4double e1 = ipotFluct; <<
180 <<
181 if(tcut > e1) { <<
182 a1 = meanLoss*(1.-rate)/e1; <<
183 if(a1 < a0) { <<
184 const G4double fwnow = 0.1+(fw-0.1)*std: <<
185 a1 /= fwnow; <<
186 e1 *= fwnow; <<
187 } else { <<
188 a1 /= fw; <<
189 e1 *= fw; <<
190 } <<
191 } <<
192 <<
193 const G4double w1 = tcut/e0; <<
194 a3 = rate*meanLoss*(tcut - e0)/(e0*tcut*G4Lo <<
195 if(a1 <= 0.) { a3 /= rate; } <<
196 <<
197 //'nearly' Gaussian fluctuation if a1>nmaxCo <<
198 G4double emean = 0.; 218 G4double emean = 0.;
199 G4double sig2e = 0.; << 219 G4double sig2e = 0., sige = 0.;
200 << 220 G4double p1 = 0., p2 = 0., p3 = 0.;
>> 221
201 // excitation of type 1 222 // excitation of type 1
202 if(a1 > 0.0) { AddExcitation(rndmEngineF, a1 << 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 }
203 235
204 if(sig2e > 0.0) { SampleGauss(rndmEngineF, e << 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 }
205 249
206 // ionisation 250 // ionisation
207 if(a3 > 0.) { << 251 G4double lossc = 0.;
208 emean = 0.; << 252 if(a3 > 0.)
209 sig2e = 0.; << 253 {
210 G4double p3 = a3; << 254 p3 = a3;
211 G4double alfa = 1.; 255 G4double alfa = 1.;
212 if(a3 > nmaxCont) { << 256 if(a3 > nmaxCont2)
213 alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3) << 257 {
214 const G4double alfa1 = alfa*G4Log(alfa) << 258 alfa = w1*(nmaxCont2+a3)/(w1*nmaxCont2+a3);
215 const G4double namean = a3*w1*(alfa-1.)/ << 259 G4double alfa1 = alfa*log(alfa)/(alfa-1.);
216 emean += namean*e0*alfa1; << 260 G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
217 sig2e += e0*e0*namean*(alfa-alfa1*alfa1) << 261 emean += namean*e0*alfa1;
218 p3 = a3 - namean; << 262 sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
>> 263 p3 = a3-namean;
219 } 264 }
220 265
221 const G4double w3 = alfa*e0; << 266 G4double w2 = alfa*e0;
222 if(tcut > w3) { << 267 G4double w = (tmax-w2)/tmax;
223 const G4double w = (tcut-w3)/tcut; << 268 G4int nb = G4Poisson(p3);
224 const G4int nnb = (G4int)G4Poisson(p3); << 269 if(nb > 0)
225 if(nnb > 0) { << 270 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
226 if(nnb > sizearray) { <<
227 sizearray = nnb; <<
228 delete [] rndmarray; <<
229 rndmarray = new G4double[nnb]; <<
230 } <<
231 rndmEngineF->flatArray(nnb, rndmarray) <<
232 for (G4int k=0; k<nnb; ++k) { loss += <<
233 } <<
234 } <<
235 if(sig2e > 0.0) { SampleGauss(rndmEngineF, <<
236 } 271 }
237 //G4cout << "### loss=" << loss << G4endl; << 272
>> 273 if(emean > 0.)
>> 274 {
>> 275 sige = sqrt(sig2e);
>> 276 loss += max(0.,G4RandGauss::shoot(emean,sige));
>> 277 }
>> 278
>> 279 loss += lossc;
>> 280
238 return loss; 281 return loss;
>> 282
239 } 283 }
240 284
241 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 286
243 287
244 G4double G4UniversalFluctuation::Dispersion( 288 G4double G4UniversalFluctuation::Dispersion(
245 const G4Material* ma 289 const G4Material* material,
246 const G4DynamicParti 290 const G4DynamicParticle* dp,
247 const G4double tcut, << 291 G4double& tmax,
248 const G4double tmax, << 292 G4double& length)
249 const G4double lengt <<
250 { 293 {
251 if(dp->GetDefinition() != particle) { Initia << 294 if(!particle) InitialiseMe(dp->GetDefinition());
252 const G4double beta = dp->GetBeta(); <<
253 return (tmax/(beta*beta) - 0.5*tcut) * CLHEP <<
254 * material->GetElectronDensity() * chargeS <<
255 } <<
256 295
257 //....oooOO0OOooo........oooOO0OOooo........oo << 296 electronDensity = material->GetElectronDensity();
258 297
259 void << 298 G4double gam = (dp->GetKineticEnergy())/particleMass + 1.0;
260 G4UniversalFluctuation::SetParticleAndCharge(c << 299 G4double beta2 = 1.0 - 1.0/(gam*gam);
261 G << 300
262 { << 301 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
263 if(part != particle) { << 302 * electronDensity * chargeSquare;
264 particle = part; << 303
265 particleMass = part->GetPDGMass(); << 304 return siga;
266 <<
267 // Derived quantities <<
268 m_Inv_particleMass = 1.0 / particleMass; <<
269 m_massrate = CLHEP::electron_mass_c2 * m_I <<
270 } <<
271 chargeSquare = q2; <<
272 } 305 }
273 306
274 //....oooOO0OOooo........oooOO0OOooo........oo 307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275 308