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