| 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 79188 2014-02-20 09:22:48Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // -------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class file 30 // GEANT4 Class file
30 // 31 //
31 // 32 //
32 // File name: G4UniversalFluctuation 33 // File name: G4UniversalFluctuation
33 // 34 //
34 // Author: V. Ivanchenko for Laszlo Urb << 35 // Author: Laszlo Urban
35 // 36 //
36 // Creation date: 03.01.2002 37 // Creation date: 03.01.2002
37 // 38 //
38 // Modifications: 39 // Modifications:
39 // 40 //
>> 41 // 28-12-02 add method Dispersion (V.Ivanchenko)
>> 42 // 07-02-03 change signature (V.Ivanchenko)
>> 43 // 13-02-03 Add name (V.Ivanchenko)
>> 44 // 16-10-03 Changed interface to Initialisation (V.Ivanchenko)
>> 45 // 07-11-03 Fix problem of rounding of double in G4UniversalFluctuations
>> 46 // 06-02-04 Add control on big sigma > 2*meanLoss (V.Ivanchenko)
>> 47 // 26-04-04 Comment out the case of very small step (V.Ivanchenko)
>> 48 // 07-02-05 define problim = 5.e-3 (mma)
>> 49 // 03-05-05 conditions of Gaussian fluctuation changed (bugfix)
>> 50 // + smearing for very small loss (L.Urban)
>> 51 // 03-10-05 energy dependent rate -> cut dependence of the
>> 52 // distribution is much weaker (L.Urban)
>> 53 // 17-10-05 correction for very small loss (L.Urban)
>> 54 // 20-03-07 'GLANDZ' part rewritten completely, no 'very small loss'
>> 55 // regime any more (L.Urban)
>> 56 // 03-04-07 correction to get better width of eloss distr.(L.Urban)
>> 57 // 13-07-07 add protection for very small step or low-density material (VI)
>> 58 // 19-03-09 new width correction (does not depend on previous steps) (L.Urban)
>> 59 // 20-03-09 modification in the width correction (L.Urban)
>> 60 // 14-06-10 fixed tail distribution - do not use uniform function (L.Urban)
>> 61 // 08-08-10 width correction algorithm has bee modified -->
>> 62 // better results for thin targets (L.Urban)
>> 63 // 06-02-11 correction for very small losses (L.Urban)
40 // 64 //
41 65
42 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 68
45 #include "G4UniversalFluctuation.hh" 69 #include "G4UniversalFluctuation.hh"
46 #include "G4PhysicalConstants.hh" 70 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh" 71 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh" 72 #include "Randomize.hh"
49 #include "G4Poisson.hh" 73 #include "G4Poisson.hh"
>> 74 #include "G4Step.hh"
50 #include "G4Material.hh" 75 #include "G4Material.hh"
51 #include "G4MaterialCutsCouple.hh" 76 #include "G4MaterialCutsCouple.hh"
52 #include "G4DynamicParticle.hh" 77 #include "G4DynamicParticle.hh"
53 #include "G4ParticleDefinition.hh" 78 #include "G4ParticleDefinition.hh"
54 #include "G4Log.hh" 79 #include "G4Log.hh"
>> 80 #include "G4Exp.hh"
55 81
56 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 83
>> 84 using namespace std;
>> 85
58 G4UniversalFluctuation::G4UniversalFluctuation 86 G4UniversalFluctuation::G4UniversalFluctuation(const G4String& nam)
59 :G4VEmFluctuationModel(nam), 87 :G4VEmFluctuationModel(nam),
60 minLoss(10.*CLHEP::eV) << 88 particle(0),
>> 89 minNumberInteractionsBohr(10.0),
>> 90 theBohrBeta2(50.0*keV/proton_mass_c2),
>> 91 minLoss(10.*eV),
>> 92 nmaxCont(16.),
>> 93 rate(0.55),
>> 94 fw(4.)
61 { 95 {
62 rndmarray = new G4double[sizearray]; << 96 lastMaterial = 0;
>> 97
>> 98 particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct
>> 99 = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall
>> 100 = e1 = e2 = 0;
>> 101 m_Inv_particleMass = m_massrate = DBL_MAX;
63 } 102 }
64 103
65 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 105
67 G4UniversalFluctuation::~G4UniversalFluctuatio 106 G4UniversalFluctuation::~G4UniversalFluctuation()
68 { << 107 {}
69 delete [] rndmarray; <<
70 } <<
71 108
72 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 110
74 void G4UniversalFluctuation::InitialiseMe(cons 111 void G4UniversalFluctuation::InitialiseMe(const G4ParticleDefinition* part)
75 { 112 {
76 particle = part; << 113 particle = part;
77 particleMass = part->GetPDGMass(); << 114 particleMass = part->GetPDGMass();
78 const G4double q = part->GetPDGCharge()/CLHE << 115 G4double q = part->GetPDGCharge()/eplus;
79 116
80 // Derived quantities 117 // Derived quantities
81 m_Inv_particleMass = 1.0 / particleMass; 118 m_Inv_particleMass = 1.0 / particleMass;
82 m_massrate = CLHEP::electron_mass_c2 * m_Inv << 119 m_massrate = electron_mass_c2 * m_Inv_particleMass ;
83 chargeSquare = q*q; << 120 chargeSquare = q*q;
84 } 121 }
85 122
86 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 124
88 G4double 125 G4double
89 G4UniversalFluctuation::SampleFluctuations(con 126 G4UniversalFluctuation::SampleFluctuations(const G4MaterialCutsCouple* couple,
90 con << 127 const G4DynamicParticle* dp,
91 con << 128 G4double tmax,
92 con << 129 G4double length,
93 con << 130 G4double averageLoss)
94 con <<
95 { 131 {
96 // Calculate actual loss from the mean loss. 132 // Calculate actual loss from the mean loss.
97 // The model used to get the fluctuations is 133 // The model used to get the fluctuations is essentially the same
98 // as in Glandz in Geant3 (Cern program libr 134 // as in Glandz in Geant3 (Cern program library W5013, phys332).
99 // L. Urban et al. NIM A362, p.416 (1995) an 135 // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
100 136
101 // shortcut for very small loss or from a st 137 // shortcut for very small loss or from a step nearly equal to the range
102 // (out of validity of the model) 138 // (out of validity of the model)
103 // 139 //
104 if (averageLoss < minLoss) { return averageL << 140 G4double meanLoss = averageLoss;
105 meanLoss = averageLoss; << 141 G4double tkin = dp->GetKineticEnergy();
106 const G4double tkin = dp->GetKineticEnergy( <<
107 //G4cout<< "Emean= "<< meanLoss<< " tmax= "< 142 //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
>> 143 if (meanLoss < minLoss) { return meanLoss; }
108 144
109 if(dp->GetDefinition() != particle) { Initia 145 if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
110 << 146
111 CLHEP::HepRandomEngine* rndmEngineF = G4Rand << 147 G4double tau = tkin * m_Inv_particleMass;
112 << 148 G4double gam = tau + 1.0;
113 const G4double gam = tkin * m_Inv_particle << 149 G4double gam2 = gam*gam;
114 const G4double gam2 = gam*gam; << 150 G4double beta2 = tau*(tau + 2.0)/gam2;
115 const G4double beta = dp->GetBeta(); <<
116 const G4double beta2 = beta*beta; <<
117 151
118 G4double loss(0.), siga(0.); 152 G4double loss(0.), siga(0.);
119 153
120 const G4Material* material = couple->GetMate 154 const G4Material* material = couple->GetMaterial();
121 155
122 // Gaussian regime 156 // Gaussian regime
123 // for heavy particles only and conditions 157 // for heavy particles only and conditions
124 // for Gauusian fluct. has been changed 158 // for Gauusian fluct. has been changed
125 // 159 //
126 if (particleMass > CLHEP::electron_mass_c2 & << 160 if ((particleMass > electron_mass_c2) &&
127 meanLoss >= minNumberInteractionsBohr*tc << 161 (meanLoss >= minNumberInteractionsBohr*tmax))
>> 162 {
>> 163 G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
>> 164 (1.+m_massrate*(2.*gam+m_massrate)) ;
>> 165 if (tmaxkine <= 2.*tmax)
>> 166 {
>> 167 electronDensity = material->GetElectronDensity();
>> 168 siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
>> 169 * electronDensity * chargeSquare);
128 170
129 siga = std::sqrt((tmax/beta2 - 0.5*tcut)*C << 171 G4double sn = meanLoss/siga;
130 length*chargeSquare*mate <<
131 const G4double sn = meanLoss/siga; <<
132 172
133 // thick target case << 173 // thick target case
134 if (sn >= 2.0) { << 174 if (sn >= 2.0) {
135 175
136 const G4double twomeanLoss = meanLoss + << 176 G4double twomeanLoss = meanLoss + meanLoss;
137 do { << 177 do {
138 loss = G4RandGauss::shoot(rndmEngineF, meanL << 178 loss = G4RandGauss::shoot(meanLoss,siga);
139 // Loop checking, 03-Aug-2015, Vladimir Ivan << 179 } while (0.0 > loss || twomeanLoss < loss);
140 } while (0.0 > loss || twomeanLoss < lo <<
141 180
142 // Gamma distribution << 181 // Gamma distribution
143 } else { << 182 } else {
144 183
145 const G4double neff = sn*sn; << 184 G4double neff = sn*sn;
146 loss = meanLoss*G4RandGamma::shoot(rndmE << 185 loss = meanLoss*G4RandGamma::shoot(neff,1.0)/neff;
>> 186 }
>> 187 //G4cout << "Gauss: " << loss << G4endl;
>> 188 return loss;
147 } 189 }
148 //G4cout << "Gauss: " << loss << G4endl; <<
149 return loss; <<
150 } 190 }
151 191
152 auto ioni = material->GetIonisation(); << 192 // Glandz regime : initialisation
153 e0 = ioni->GetEnergy0fluct(); << 193 //
>> 194 if (material != lastMaterial) {
>> 195 f1Fluct = material->GetIonisation()->GetF1fluct();
>> 196 f2Fluct = material->GetIonisation()->GetF2fluct();
>> 197 e1Fluct = material->GetIonisation()->GetEnergy1fluct();
>> 198 e2Fluct = material->GetIonisation()->GetEnergy2fluct();
>> 199 e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
>> 200 e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
>> 201 ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy();
>> 202 ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
>> 203 e0 = material->GetIonisation()->GetEnergy0fluct();
>> 204 esmall = 0.5*sqrt(e0*ipotFluct);
>> 205 lastMaterial = material;
>> 206 }
154 207
155 // very small step or low-density material 208 // very small step or low-density material
156 if(tcut <= e0) { return meanLoss; } << 209 if(tmax <= e0) { return meanLoss; }
157 210
158 ipotFluct = ioni->GetMeanExcitationEnergy(); << 211 G4double losstot = 0.;
159 ipotLogFluct = ioni->GetLogMeanExcEnergy(); << 212 G4int nstep = 1;
>> 213 if(meanLoss < 25.*ipotFluct)
>> 214 {
>> 215 if(G4UniformRand()*ipotFluct< 0.04*meanLoss)
>> 216 { nstep = 1; }
>> 217 else
>> 218 {
>> 219 nstep = 2;
>> 220 meanLoss *= 0.5;
>> 221 }
>> 222 }
160 223
161 // width correction for small cuts << 224 for (G4int istep=0; istep < nstep; ++istep) {
162 const G4double scaling = std::min(1.+0.5*CLH << 225
163 meanLoss /= scaling; << 226 loss = 0.;
164 << 227
165 w2 = (tcut > ipotFluct) ? << 228 G4double a1 = 0. , a2 = 0., a3 = 0. ;
166 G4Log(2.*CLHEP::electron_mass_c2*beta2*gam << 229
167 return SampleGlandz(rndmEngineF, material, t << 230 if(tmax > ipotFluct) {
168 } << 231 G4double w2 = G4Log(2.*electron_mass_c2*beta2*gam2)-beta2;
>> 232
>> 233 if(w2 > ipotLogFluct) {
>> 234 G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
>> 235 a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
>> 236 if(w2 > e2LogFluct) {
>> 237 a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
>> 238 }
>> 239 if(a1 < nmaxCont) {
>> 240 //small energy loss
>> 241 G4double sa1 = sqrt(a1);
>> 242 if(G4UniformRand() < G4Exp(-sa1))
>> 243 {
>> 244 e1 = esmall;
>> 245 a1 = meanLoss*(1.-rate)/e1;
>> 246 a2 = 0.;
>> 247 e2 = e2Fluct;
>> 248 }
>> 249 else
>> 250 {
>> 251 a1 = sa1 ;
>> 252 e1 = sa1*e1Fluct;
>> 253 e2 = e2Fluct;
>> 254 }
>> 255
>> 256 } else {
>> 257 //not small energy loss
>> 258 //correction to get better fwhm value
>> 259 a1 /= fw;
>> 260 e1 = fw*e1Fluct;
>> 261 e2 = e2Fluct;
>> 262 }
>> 263 }
>> 264 }
169 265
170 //....oooOO0OOooo........oooOO0OOooo........oo << 266 G4double w1 = tmax/e0;
>> 267 if(tmax > e0) {
>> 268 a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*G4Log(w1));
>> 269 if(a1+a2 <= 0.) {
>> 270 a3 /= rate;
>> 271 }
>> 272 }
>> 273 //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont
>> 274 G4double emean = 0.;
>> 275 G4double sig2e = 0., sige = 0.;
>> 276 G4double p1 = 0., p2 = 0., p3 = 0.;
>> 277
>> 278 // excitation of type 1
>> 279 if(a1 > nmaxCont)
>> 280 {
>> 281 emean += a1*e1;
>> 282 sig2e += a1*e1*e1;
>> 283 }
>> 284 else if(a1 > 0.)
>> 285 {
>> 286 p1 = G4double(G4Poisson(a1));
>> 287 loss += p1*e1;
>> 288 if(p1 > 0.) {
>> 289 loss += (1.-2.*G4UniformRand())*e1;
>> 290 }
>> 291 }
171 292
172 G4double << 293 // excitation of type 2
173 G4UniversalFluctuation::SampleGlandz(CLHEP::He << 294 if(a2 > nmaxCont)
174 const G4M << 295 {
175 const G4d << 296 emean += a2*e2;
176 { << 297 sig2e += a2*e2*e2;
177 G4double a1(0.0), a3(0.0); << 298 }
178 G4double loss = 0.0; << 299 else if(a2 > 0.)
179 G4double e1 = ipotFluct; << 300 {
180 << 301 p2 = G4double(G4Poisson(a2));
181 if(tcut > e1) { << 302 loss += p2*e2;
182 a1 = meanLoss*(1.-rate)/e1; << 303 if(p2 > 0.)
183 if(a1 < a0) { << 304 loss += (1.-2.*G4UniformRand())*e2;
184 const G4double fwnow = 0.1+(fw-0.1)*std: << 305 }
185 a1 /= fwnow; << 306 if(emean > 0.)
186 e1 *= fwnow; << 307 {
187 } else { << 308 sige = sqrt(sig2e);
188 a1 /= fw; << 309 loss += max(0.,G4RandGauss::shoot(emean,sige));
189 e1 *= fw; << 310 }
190 } <<
191 } <<
192 311
193 const G4double w1 = tcut/e0; << 312 // ionisation
194 a3 = rate*meanLoss*(tcut - e0)/(e0*tcut*G4Lo << 313 G4double lossc = 0.;
195 if(a1 <= 0.) { a3 /= rate; } << 314 if(a3 > 0.) {
196 << 315 emean = 0.;
197 //'nearly' Gaussian fluctuation if a1>nmaxCo << 316 sig2e = 0.;
198 G4double emean = 0.; << 317 sige = 0.;
199 G4double sig2e = 0.; << 318 p3 = a3;
200 << 319 G4double alfa = 1.;
201 // excitation of type 1 << 320 if(a3 > nmaxCont)
202 if(a1 > 0.0) { AddExcitation(rndmEngineF, a1 << 321 {
203 << 322 alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
204 if(sig2e > 0.0) { SampleGauss(rndmEngineF, e << 323 G4double alfa1 = alfa*G4Log(alfa)/(alfa-1.);
205 << 324 G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
206 // ionisation << 325 emean += namean*e0*alfa1;
207 if(a3 > 0.) { << 326 sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
208 emean = 0.; << 327 p3 = a3-namean;
209 sig2e = 0.; << 328 }
210 G4double p3 = a3; << 329
211 G4double alfa = 1.; << 330 G4double w2 = alfa*e0;
212 if(a3 > nmaxCont) { << 331 G4double w = (tmax-w2)/tmax;
213 alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3) << 332 G4int nb = G4Poisson(p3);
214 const G4double alfa1 = alfa*G4Log(alfa) << 333 if(nb > 0) {
215 const G4double namean = a3*w1*(alfa-1.)/ << 334 for (G4int k=0; k<nb; k++) { lossc += w2/(1.-w*G4UniformRand()); }
216 emean += namean*e0*alfa1; << 335 }
217 sig2e += e0*e0*namean*(alfa-alfa1*alfa1) << 336
218 p3 = a3 - namean; << 337 if(emean > 0.)
219 } << 338 {
220 << 339 sige = sqrt(sig2e);
221 const G4double w3 = alfa*e0; << 340 lossc += max(0.,G4RandGauss::shoot(emean,sige));
222 if(tcut > w3) { <<
223 const G4double w = (tcut-w3)/tcut; <<
224 const G4int nnb = (G4int)G4Poisson(p3); <<
225 if(nnb > 0) { <<
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 } 341 }
234 } 342 }
235 if(sig2e > 0.0) { SampleGauss(rndmEngineF, << 343
>> 344 loss += lossc;
>> 345
>> 346 losstot += loss;
236 } 347 }
237 //G4cout << "### loss=" << loss << G4endl; << 348 //G4cout << "Vavilov: " << losstot << " Nstep= " << nstep << G4endl;
238 return loss; << 349
>> 350 return losstot;
>> 351
239 } 352 }
240 353
241 //....oooOO0OOooo........oooOO0OOooo........oo 354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 355
243 356
244 G4double G4UniversalFluctuation::Dispersion( 357 G4double G4UniversalFluctuation::Dispersion(
245 const G4Material* ma 358 const G4Material* material,
246 const G4DynamicParti 359 const G4DynamicParticle* dp,
247 const G4double tcut, << 360 G4double tmax,
248 const G4double tmax, << 361 G4double length)
249 const G4double lengt <<
250 { 362 {
251 if(dp->GetDefinition() != particle) { Initia 363 if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
252 const G4double beta = dp->GetBeta(); << 364
253 return (tmax/(beta*beta) - 0.5*tcut) * CLHEP << 365 electronDensity = material->GetElectronDensity();
254 * material->GetElectronDensity() * chargeS << 366
>> 367 G4double gam = (dp->GetKineticEnergy())*m_Inv_particleMass + 1.0;
>> 368 G4double beta2 = 1.0 - 1.0/(gam*gam);
>> 369
>> 370 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
>> 371 * electronDensity * chargeSquare;
>> 372
>> 373 return siga;
255 } 374 }
256 375
257 //....oooOO0OOooo........oooOO0OOooo........oo 376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
258 377
259 void 378 void
260 G4UniversalFluctuation::SetParticleAndCharge(c 379 G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,
261 G << 380 G4double q2)
262 { 381 {
263 if(part != particle) { 382 if(part != particle) {
264 particle = part; << 383 particle = part;
265 particleMass = part->GetPDGMass(); << 384 particleMass = part->GetPDGMass();
266 385
267 // Derived quantities 386 // Derived quantities
268 m_Inv_particleMass = 1.0 / particleMass; << 387 if( particleMass != 0.0 ){
269 m_massrate = CLHEP::electron_mass_c2 * m_I << 388 m_Inv_particleMass = 1.0 / particleMass;
>> 389 m_massrate = electron_mass_c2 * m_Inv_particleMass ;
>> 390 }else{
>> 391 m_Inv_particleMass = DBL_MAX;
>> 392 m_massrate = DBL_MAX;
>> 393 }
270 } 394 }
271 chargeSquare = q2; 395 chargeSquare = q2;
272 } 396 }
273 397
274 //....oooOO0OOooo........oooOO0OOooo........oo 398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275 399