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