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