[1316] | 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 | // |
---|
| 26 | // $Id: G4ICRU73QOModel.cc,v 1.3 2010/06/04 09:08:09 vnivanch Exp $ |
---|
[1337] | 27 | // GEANT4 tag $Name: geant4-09-04-beta-01 $ |
---|
[1316] | 28 | // |
---|
| 29 | // ------------------------------------------------------------------- |
---|
| 30 | // |
---|
| 31 | // GEANT4 Class file |
---|
| 32 | // |
---|
| 33 | // |
---|
| 34 | // File name: G4ICRU73QOModel |
---|
| 35 | // |
---|
| 36 | // Author: Alexander Bagulya |
---|
| 37 | // |
---|
| 38 | // Creation date: 21.05.2010 |
---|
| 39 | // |
---|
| 40 | // Modifications: |
---|
| 41 | // |
---|
| 42 | // |
---|
| 43 | // ------------------------------------------------------------------- |
---|
| 44 | // |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 48 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 49 | |
---|
| 50 | #include "G4ICRU73QOModel.hh" |
---|
| 51 | #include "Randomize.hh" |
---|
| 52 | #include "G4Electron.hh" |
---|
| 53 | #include "G4ParticleChangeForLoss.hh" |
---|
| 54 | #include "G4LossTableManager.hh" |
---|
| 55 | #include "G4AntiProton.hh" |
---|
| 56 | |
---|
| 57 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 58 | |
---|
| 59 | using namespace std; |
---|
| 60 | |
---|
| 61 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 62 | |
---|
| 63 | G4ICRU73QOModel::G4ICRU73QOModel(const G4ParticleDefinition* p, const G4String& nam) |
---|
| 64 | : G4VEmModel(nam), |
---|
| 65 | particle(0), |
---|
| 66 | isInitialised(false) |
---|
| 67 | { |
---|
| 68 | if(p) SetParticle(p); |
---|
| 69 | SetHighEnergyLimit(10.0*MeV); |
---|
| 70 | |
---|
| 71 | lowestKinEnergy = 10.0*keV; |
---|
| 72 | |
---|
| 73 | sizeL0 = 67; |
---|
| 74 | sizeL1 = 22; |
---|
| 75 | sizeL2 = 14; |
---|
| 76 | |
---|
| 77 | theElectron = G4Electron::Electron(); |
---|
| 78 | |
---|
| 79 | for (G4int i = 0; i < 100; ++i) |
---|
| 80 | { |
---|
| 81 | indexZ[i] = -1; |
---|
| 82 | } |
---|
| 83 | for(G4int i = 0; i < NQOELEM; ++i) |
---|
| 84 | { |
---|
| 85 | if(ZElementAvailable[i] > 0) { |
---|
| 86 | indexZ[ZElementAvailable[i]] = i; |
---|
| 87 | } |
---|
| 88 | } |
---|
| 89 | } |
---|
| 90 | |
---|
| 91 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 92 | |
---|
| 93 | G4ICRU73QOModel::~G4ICRU73QOModel() |
---|
| 94 | {} |
---|
| 95 | |
---|
| 96 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 97 | |
---|
| 98 | G4double G4ICRU73QOModel::MinEnergyCut(const G4ParticleDefinition*, |
---|
| 99 | const G4MaterialCutsCouple* ) |
---|
| 100 | { |
---|
| 101 | // return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy(); |
---|
| 102 | return 100*keV; |
---|
| 103 | } |
---|
| 104 | |
---|
| 105 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 106 | |
---|
| 107 | void G4ICRU73QOModel::Initialise(const G4ParticleDefinition* p, |
---|
| 108 | const G4DataVector&) |
---|
| 109 | { |
---|
| 110 | if(p != particle) SetParticle(p); |
---|
| 111 | |
---|
| 112 | // always false before the run |
---|
| 113 | SetDeexcitationFlag(false); |
---|
| 114 | |
---|
| 115 | if(!isInitialised) { |
---|
| 116 | isInitialised = true; |
---|
| 117 | |
---|
| 118 | G4String pname = particle->GetParticleName(); |
---|
| 119 | fParticleChange = GetParticleChangeForLoss(); |
---|
| 120 | const G4MaterialTable* mtab = G4Material::GetMaterialTable(); |
---|
| 121 | denEffData = (*mtab)[0]->GetIonisation()->GetDensityEffectData(); |
---|
| 122 | } |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 126 | |
---|
| 127 | G4double G4ICRU73QOModel::ComputeCrossSectionPerElectron( |
---|
| 128 | const G4ParticleDefinition* p, |
---|
| 129 | G4double kineticEnergy, |
---|
| 130 | G4double cutEnergy, |
---|
| 131 | G4double maxKinEnergy) |
---|
| 132 | { |
---|
| 133 | G4double cross = 0.0; |
---|
| 134 | G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); |
---|
| 135 | G4double maxEnergy = std::min(tmax,maxKinEnergy); |
---|
| 136 | if(cutEnergy < maxEnergy) { |
---|
| 137 | |
---|
| 138 | G4double energy = kineticEnergy + mass; |
---|
| 139 | G4double energy2 = energy*energy; |
---|
| 140 | G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2; |
---|
| 141 | cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax; |
---|
| 142 | |
---|
| 143 | cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2; |
---|
| 144 | } |
---|
| 145 | |
---|
| 146 | return cross; |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 150 | |
---|
| 151 | G4double G4ICRU73QOModel::ComputeCrossSectionPerAtom( |
---|
| 152 | const G4ParticleDefinition* p, |
---|
| 153 | G4double kineticEnergy, |
---|
| 154 | G4double Z, G4double, |
---|
| 155 | G4double cutEnergy, |
---|
| 156 | G4double maxEnergy) |
---|
| 157 | { |
---|
| 158 | G4double cross = Z*ComputeCrossSectionPerElectron |
---|
| 159 | (p,kineticEnergy,cutEnergy,maxEnergy); |
---|
| 160 | return cross; |
---|
| 161 | } |
---|
| 162 | |
---|
| 163 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 164 | |
---|
| 165 | G4double G4ICRU73QOModel::CrossSectionPerVolume( |
---|
| 166 | const G4Material* material, |
---|
| 167 | const G4ParticleDefinition* p, |
---|
| 168 | G4double kineticEnergy, |
---|
| 169 | G4double cutEnergy, |
---|
| 170 | G4double maxEnergy) |
---|
| 171 | { |
---|
| 172 | G4double eDensity = material->GetElectronDensity(); |
---|
| 173 | G4double cross = eDensity*ComputeCrossSectionPerElectron |
---|
| 174 | (p,kineticEnergy,cutEnergy,maxEnergy); |
---|
| 175 | return cross; |
---|
| 176 | } |
---|
| 177 | |
---|
| 178 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 179 | |
---|
| 180 | G4double G4ICRU73QOModel::ComputeDEDXPerVolume(const G4Material* material, |
---|
| 181 | const G4ParticleDefinition* p, |
---|
| 182 | G4double kineticEnergy, |
---|
| 183 | G4double cutEnergy) |
---|
| 184 | { |
---|
| 185 | SetParticle(p); |
---|
| 186 | G4double tmax = MaxSecondaryEnergy(p, kineticEnergy); |
---|
| 187 | G4double tkin = kineticEnergy/massRate; |
---|
| 188 | G4double dedx = 0.0; |
---|
| 189 | if(tkin > lowestKinEnergy) { dedx = DEDX(material, tkin); } |
---|
| 190 | else { dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); } |
---|
| 191 | |
---|
| 192 | if (cutEnergy < tmax) { |
---|
| 193 | |
---|
| 194 | G4double tau = kineticEnergy/mass; |
---|
| 195 | G4double gam = tau + 1.0; |
---|
| 196 | G4double bg2 = tau * (tau+2.0); |
---|
| 197 | G4double beta2 = bg2/(gam*gam); |
---|
| 198 | G4double x = cutEnergy/tmax; |
---|
| 199 | |
---|
| 200 | dedx += chargeSquare*( log(x) + (1.0 - x)*beta2 ) * twopi_mc2_rcl2 |
---|
| 201 | * material->GetElectronDensity()/beta2; |
---|
| 202 | } |
---|
| 203 | |
---|
| 204 | return dedx; |
---|
| 205 | } |
---|
| 206 | |
---|
| 207 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 208 | |
---|
| 209 | G4double G4ICRU73QOModel::DEDX(const G4Material* material, |
---|
| 210 | G4double kineticEnergy) |
---|
| 211 | { |
---|
| 212 | G4double eloss = 0.0; |
---|
| 213 | const G4int numberOfElements = material->GetNumberOfElements(); |
---|
| 214 | const G4double* theAtomicNumDensityVector = |
---|
| 215 | material->GetAtomicNumDensityVector(); |
---|
| 216 | |
---|
| 217 | // Bragg's rule calculation |
---|
| 218 | const G4ElementVector* theElementVector = |
---|
| 219 | material->GetElementVector() ; |
---|
| 220 | |
---|
| 221 | // loop for the elements in the material |
---|
| 222 | for (G4int i=0; i<numberOfElements; ++i) |
---|
| 223 | { |
---|
| 224 | const G4Element* element = (*theElementVector)[i] ; |
---|
| 225 | eloss += DEDXPerElement(G4int(element->GetZ()), kineticEnergy) |
---|
| 226 | * theAtomicNumDensityVector[i] * G4int(element->GetZ()); |
---|
| 227 | } |
---|
| 228 | return eloss; |
---|
| 229 | } |
---|
| 230 | |
---|
| 231 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 232 | |
---|
| 233 | G4double G4ICRU73QOModel::DEDXPerElement(G4int AtomicNumber, |
---|
| 234 | G4double kineticEnergy) |
---|
| 235 | { |
---|
| 236 | G4int Z = AtomicNumber; |
---|
| 237 | if(Z > 97) { Z = 97; } |
---|
| 238 | G4int nbOfShells = GetNumberOfShells(Z); |
---|
| 239 | if(nbOfShells < 1) nbOfShells = 1; |
---|
| 240 | |
---|
| 241 | G4double v = CLHEP::c_light * std::sqrt( 2.0*kineticEnergy/proton_mass_c2 ); |
---|
| 242 | |
---|
| 243 | G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/v; |
---|
| 244 | |
---|
| 245 | G4double tau = kineticEnergy/proton_mass_c2; |
---|
| 246 | G4double gam = tau + 1.0; |
---|
| 247 | G4double bg2 = tau * (tau+2.0); |
---|
| 248 | G4double beta2 = bg2/(gam*gam); |
---|
| 249 | |
---|
| 250 | G4double l0Term = 0, l1Term = 0, l2Term = 0; |
---|
| 251 | |
---|
| 252 | for (G4int nos = 0; nos < nbOfShells; ++nos){ |
---|
| 253 | |
---|
| 254 | G4double NormalizedEnergy = (2.0*CLHEP::electron_mass_c2*beta2) / |
---|
| 255 | GetShellEnergy(Z,nos); |
---|
| 256 | |
---|
| 257 | G4double shStrength = GetShellStrength(Z,nos); |
---|
| 258 | |
---|
| 259 | G4double l0 = GetL0(NormalizedEnergy); |
---|
| 260 | l0Term += shStrength * l0; |
---|
| 261 | |
---|
| 262 | G4double l1 = GetL1(NormalizedEnergy); |
---|
| 263 | l1Term += shStrength * l1; |
---|
| 264 | |
---|
| 265 | G4double l2 = GetL2(NormalizedEnergy); |
---|
| 266 | l2Term += shStrength * l2; |
---|
| 267 | |
---|
| 268 | } |
---|
| 269 | G4double dedx = 2*CLHEP::twopi_mc2_rcl2*chargeSquare*factorBethe[Z]* |
---|
| 270 | (l0Term + q*fBetheVelocity*l1Term |
---|
| 271 | + chargeSquare*fBetheVelocity*fBetheVelocity*l2Term)/beta2; |
---|
| 272 | return dedx; |
---|
| 273 | } |
---|
| 274 | |
---|
| 275 | |
---|
| 276 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 277 | |
---|
| 278 | G4double G4ICRU73QOModel::GetOscillatorEnergy(G4int Z, |
---|
| 279 | G4int nbOfTheShell) const |
---|
| 280 | { |
---|
| 281 | G4int idx = denEffData->GetElementIndex(Z, kStateUndefined); |
---|
| 282 | if(idx == -1) { idx = denEffData->GetElementIndex(Z-1, kStateUndefined); } |
---|
| 283 | G4double PlasmaEnergy = denEffData->GetPlasmaEnergy(idx); |
---|
| 284 | |
---|
| 285 | G4double PlasmaEnergy2 = PlasmaEnergy * PlasmaEnergy; |
---|
| 286 | |
---|
| 287 | G4double plasmonTerm = 0.66667 * G4AtomicShells::GetNumberOfElectrons(Z,nbOfTheShell) |
---|
| 288 | * PlasmaEnergy2 / (Z*Z) ; |
---|
| 289 | |
---|
| 290 | G4double ionTerm = std::exp(0.5) * (G4AtomicShells::GetBindingEnergy(Z,nbOfTheShell)) ; |
---|
| 291 | G4double ionTerm2 = ionTerm*ionTerm ; |
---|
| 292 | |
---|
| 293 | G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm ); |
---|
| 294 | |
---|
| 295 | return oscShellEnergy; |
---|
| 296 | } |
---|
| 297 | |
---|
| 298 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 299 | |
---|
| 300 | G4double G4ICRU73QOModel::GetL0(G4double normEnergy) const |
---|
| 301 | { |
---|
| 302 | G4int n; |
---|
| 303 | |
---|
| 304 | for(n = 0; n < sizeL0; n++) { |
---|
| 305 | if( normEnergy < L0[n][0] ) break; |
---|
| 306 | } |
---|
| 307 | if(0 == n) n = 1 ; |
---|
| 308 | if(n >= sizeL0) n = sizeL0 - 1 ; |
---|
| 309 | |
---|
| 310 | G4double l0 = L0[n][1]; |
---|
| 311 | G4double l0p = L0[n-1][1]; |
---|
| 312 | G4double bethe = l0p + (l0 - l0p) * ( normEnergy - L0[n-1][0]) / |
---|
| 313 | (L0[n][0] - L0[n-1][0]); |
---|
| 314 | |
---|
| 315 | return bethe ; |
---|
| 316 | } |
---|
| 317 | |
---|
| 318 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 319 | |
---|
| 320 | G4double G4ICRU73QOModel::GetL1(G4double normEnergy) const |
---|
| 321 | { |
---|
| 322 | G4int n; |
---|
| 323 | |
---|
| 324 | for(n = 0; n < sizeL1; n++) { |
---|
| 325 | if( normEnergy < L1[n][0] ) break; |
---|
| 326 | } |
---|
| 327 | if(0 == n) n = 1 ; |
---|
| 328 | if(n >= sizeL1) n = sizeL1 - 1 ; |
---|
| 329 | |
---|
| 330 | G4double l1 = L1[n][1]; |
---|
| 331 | G4double l1p = L1[n-1][1]; |
---|
| 332 | G4double barkas= l1p + (l1 - l1p) * ( normEnergy - L1[n-1][0]) / |
---|
| 333 | (L1[n][0] - L1[n-1][0]); |
---|
| 334 | |
---|
| 335 | return barkas; |
---|
| 336 | } |
---|
| 337 | |
---|
| 338 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 339 | |
---|
| 340 | G4double G4ICRU73QOModel::GetL2(G4double normEnergy) const |
---|
| 341 | { |
---|
| 342 | G4int n; |
---|
| 343 | for(n = 0; n < sizeL2; n++) { |
---|
| 344 | if( normEnergy < L2[n][0] ) break; |
---|
| 345 | } |
---|
| 346 | if(0 == n) n = 1 ; |
---|
| 347 | if(n >= sizeL2) n = sizeL2 - 1 ; |
---|
| 348 | |
---|
| 349 | G4double l2 = L2[n][1]; |
---|
| 350 | G4double l2p = L2[n-1][1]; |
---|
| 351 | G4double bloch = l2p + (l2 - l2p) * ( normEnergy - L2[n-1][0]) / |
---|
| 352 | (L2[n][0] - L2[n-1][0]); |
---|
| 353 | |
---|
| 354 | return bloch; |
---|
| 355 | } |
---|
| 356 | |
---|
| 357 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 358 | |
---|
| 359 | void G4ICRU73QOModel::CorrectionsAlongStep(const G4MaterialCutsCouple*, |
---|
| 360 | const G4DynamicParticle*, |
---|
| 361 | G4double&, |
---|
| 362 | G4double&, |
---|
| 363 | G4double) |
---|
| 364 | {} |
---|
| 365 | |
---|
| 366 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 367 | |
---|
| 368 | void G4ICRU73QOModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, |
---|
| 369 | const G4MaterialCutsCouple*, |
---|
| 370 | const G4DynamicParticle* dp, |
---|
| 371 | G4double xmin, |
---|
| 372 | G4double maxEnergy) |
---|
| 373 | { |
---|
| 374 | G4double tmax = MaxSecondaryKinEnergy(dp); |
---|
| 375 | G4double xmax = std::min(tmax, maxEnergy); |
---|
| 376 | if(xmin >= xmax) return; |
---|
| 377 | |
---|
| 378 | G4double kineticEnergy = dp->GetKineticEnergy(); |
---|
| 379 | G4double energy = kineticEnergy + mass; |
---|
| 380 | G4double energy2 = energy*energy; |
---|
| 381 | G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2; |
---|
| 382 | G4double grej = 1.0; |
---|
| 383 | G4double deltaKinEnergy, f; |
---|
| 384 | |
---|
| 385 | G4ThreeVector direction = dp->GetMomentumDirection(); |
---|
| 386 | |
---|
| 387 | // sampling follows ... |
---|
| 388 | do { |
---|
| 389 | G4double q = G4UniformRand(); |
---|
| 390 | deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q); |
---|
| 391 | |
---|
| 392 | f = 1.0 - beta2*deltaKinEnergy/tmax; |
---|
| 393 | |
---|
| 394 | if(f > grej) { |
---|
| 395 | G4cout << "G4ICRU73QOModel::SampleSecondary Warning! " |
---|
| 396 | << "Majorant " << grej << " < " |
---|
| 397 | << f << " for e= " << deltaKinEnergy |
---|
| 398 | << G4endl; |
---|
| 399 | } |
---|
| 400 | |
---|
| 401 | } while( grej*G4UniformRand() >= f ); |
---|
| 402 | |
---|
| 403 | G4double deltaMomentum = |
---|
| 404 | sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2)); |
---|
| 405 | G4double totMomentum = energy*sqrt(beta2); |
---|
| 406 | G4double cost = deltaKinEnergy * (energy + electron_mass_c2) / |
---|
| 407 | (deltaMomentum * totMomentum); |
---|
| 408 | if(cost > 1.0) cost = 1.0; |
---|
| 409 | G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); |
---|
| 410 | |
---|
| 411 | G4double phi = twopi * G4UniformRand() ; |
---|
| 412 | |
---|
| 413 | G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ; |
---|
| 414 | deltaDirection.rotateUz(direction); |
---|
| 415 | |
---|
| 416 | // Change kinematics of primary particle |
---|
| 417 | kineticEnergy -= deltaKinEnergy; |
---|
| 418 | G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum; |
---|
| 419 | finalP = finalP.unit(); |
---|
| 420 | |
---|
| 421 | fParticleChange->SetProposedKineticEnergy(kineticEnergy); |
---|
| 422 | fParticleChange->SetProposedMomentumDirection(finalP); |
---|
| 423 | |
---|
| 424 | // create G4DynamicParticle object for delta ray |
---|
| 425 | G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection, |
---|
| 426 | deltaKinEnergy); |
---|
| 427 | |
---|
| 428 | vdp->push_back(delta); |
---|
| 429 | } |
---|
| 430 | |
---|
| 431 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 432 | |
---|
| 433 | G4double G4ICRU73QOModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd, |
---|
| 434 | G4double kinEnergy) |
---|
| 435 | { |
---|
| 436 | if(pd != particle) SetParticle(pd); |
---|
| 437 | G4double tau = kinEnergy/mass; |
---|
| 438 | G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) / |
---|
| 439 | (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); |
---|
| 440 | return tmax; |
---|
| 441 | } |
---|
| 442 | |
---|
| 443 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 444 | |
---|
| 445 | const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] = {1,2,4,6,7,8,10,13,14,-18, |
---|
| 446 | 22,26,28,29,32,36,42,47, |
---|
| 447 | 50,54,73,74,78,79,82,92}; |
---|
| 448 | |
---|
| 449 | const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] = {1,1,2,3,3,3,3,4,5,4, |
---|
| 450 | 5,5,5,5,6,4,6,6, |
---|
| 451 | 7,6,6,8,7,7,9,9}; |
---|
| 452 | |
---|
| 453 | const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] = {0,1,2,4,7,10,13,16,20,25, |
---|
| 454 | 29,34,39,44,49,55,59,65, |
---|
| 455 | 71,78,84,90,98,105,112,121}; |
---|
| 456 | |
---|
| 457 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 458 | |
---|
| 459 | // SubShellOccupation = Z * ShellStrength |
---|
| 460 | const G4double G4ICRU73QOModel::SubShellOccupation[NQODATA] = |
---|
| 461 | { |
---|
| 462 | 1.000, // H 0 |
---|
| 463 | 2.000, // He 1 |
---|
| 464 | 1.930, 2.070, // Be 2-3 |
---|
| 465 | 1.992, 1.841, 2.167, // C 4-6 |
---|
| 466 | 1.741, 1.680, 3.579, // N 7-9 |
---|
| 467 | 1.802, 1.849, 4.349, // O 10-12 |
---|
| 468 | 1.788, 2.028, 6.184, // Ne 13-15 |
---|
| 469 | 1.623, 2.147, 6.259, 2.971, // Al 16-19 |
---|
| 470 | 1.631, 2.094, 6.588, 2.041, 1.646, // Si 20-24 |
---|
| 471 | 1.535, 8.655, 1.706, 6.104, // Ar 25-28 |
---|
| 472 | 1.581, 8.358, 8.183, 2.000, 1.878, // Ti 29-33 |
---|
| 473 | 1.516, 8.325, 8.461, 6.579, 1.119, // Fe 34-38 |
---|
| 474 | 1.422, 7.81, 8.385, 8.216, 2.167, // Ni 39-43 |
---|
| 475 | 1.458, 8.049, 8.79, 9.695, 1.008, // Cu 44-48 |
---|
| 476 | 1.442, 7.791, 7.837, 10.122, 2.463, 2.345, // Ge 49-54 |
---|
| 477 | 1.645, 7.765, 19.192, 7.398, // Kr 55-58 |
---|
| 478 | 1.313, 6.409, 19.229, 8.633, 5.036, 1.380, // Mo 59-64 |
---|
| 479 | 1.295, 6.219, 18.751, 8.748, 10.184, 1.803, // Ag 65-70 |
---|
| 480 | 1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948, // Sn 71-77 |
---|
| 481 | 1.563, 6.312, 21.868, 5.762, 11.245, 7.250, // Xe 78-83 |
---|
| 482 | 0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284, // Ta 84-89 |
---|
| 483 | 1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108, // W 90-97 |
---|
| 484 | 1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025, // Pt 98-104 |
---|
| 485 | 1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322, // Au 105-111 |
---|
| 486 | 2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000, // Pb 112-120 |
---|
| 487 | 2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000 // U 121-129 |
---|
| 488 | }; |
---|
| 489 | |
---|
| 490 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 491 | |
---|
| 492 | // ShellEnergy in eV |
---|
| 493 | const G4double G4ICRU73QOModel::ShellEnergy[NQODATA] = |
---|
| 494 | { |
---|
| 495 | 19.2, // H |
---|
| 496 | 41.8, // He |
---|
| 497 | 209.11, 21.68, // Be |
---|
| 498 | 486.2, 60.95, 23.43, // C |
---|
| 499 | 732.61, 100.646, 23.550, // N |
---|
| 500 | 965.1, 129.85, 31.60, // O |
---|
| 501 | 1525.9, 234.9, 56.18, // Ne |
---|
| 502 | 2701, 476.5, 150.42, 16.89, // Al |
---|
| 503 | 3206.1, 586.4, 186.8, 23.52, 14.91, // Si |
---|
| 504 | 5551.6, 472.43, 124.85, 22.332, // Ar |
---|
| 505 | 8554.6, 850.58, 93.47, 39.19, 19.46, // Ti |
---|
| 506 | 12254.7, 1279.29, 200.35, 49.19, 17.66, // Fe |
---|
| 507 | 14346.9, 1532.28, 262.71, 74.37, 23.03, // Ni |
---|
| 508 | 15438.5, 1667.96, 294.1, 70.69, 16.447, // Cu |
---|
| 509 | 19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95, // Ge |
---|
| 510 | 24643, 2906.4, 366.85, 22.24, // Kr |
---|
| 511 | 34394, 4365.3, 589.36, 129.42, 35.59, 18.42, // Mo |
---|
| 512 | 43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63, // Ag |
---|
| 513 | 49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54, // Sn |
---|
| 514 | 58987, 8159, 1296.6, 356.75, 101.03, 16.52, // Xe |
---|
| 515 | 88926, 18012, 3210, 575, 108.7, 30.8, // Ta |
---|
| 516 | 115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25, // W |
---|
| 517 | 128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04, // Pt |
---|
| 518 | 131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575, // Au |
---|
| 519 | 154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17, // Pb |
---|
| 520 | 167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43 // U |
---|
| 521 | }; |
---|
| 522 | |
---|
| 523 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 524 | |
---|
| 525 | // Data for L0 from: Sigmund P., Haagerup U. Phys. Rev. A34 (1986) 892-910 |
---|
| 526 | const G4double G4ICRU73QOModel::L0[67][2] = |
---|
| 527 | { |
---|
| 528 | {0.00, 0.000001}, |
---|
| 529 | {0.10, 0.000001}, |
---|
| 530 | {0.12, 0.00001}, |
---|
| 531 | {0.14, 0.00005}, |
---|
| 532 | {0.16, 0.00014}, |
---|
| 533 | {0.18, 0.00030}, |
---|
| 534 | {0.20, 0.00057}, |
---|
| 535 | {0.25, 0.00189}, |
---|
| 536 | {0.30, 0.00429}, |
---|
| 537 | {0.35, 0.00784}, |
---|
| 538 | {0.40, 0.01248}, |
---|
| 539 | {0.45, 0.01811}, |
---|
| 540 | {0.50, 0.02462}, |
---|
| 541 | {0.60, 0.03980}, |
---|
| 542 | {0.70, 0.05731}, |
---|
| 543 | {0.80, 0.07662}, |
---|
| 544 | {0.90, 0.09733}, |
---|
| 545 | {1.00, 0.11916}, |
---|
| 546 | {1.20, 0.16532}, |
---|
| 547 | {1.40, 0.21376}, |
---|
| 548 | {1.60, 0.26362}, |
---|
| 549 | {1.80, 0.31428}, |
---|
| 550 | {2.00, 0.36532}, |
---|
| 551 | {2.50, 0.49272}, |
---|
| 552 | {3.00, 0.61765}, |
---|
| 553 | {3.50, 0.73863}, |
---|
| 554 | {4.00, 0.85496}, |
---|
| 555 | {4.50, 0.96634}, |
---|
| 556 | {5.00, 1.07272}, |
---|
| 557 | {6.00, 1.27086}, |
---|
| 558 | {7.00, 1.45075}, |
---|
| 559 | {8.00, 1.61412}, |
---|
| 560 | {9.00, 1.76277}, |
---|
| 561 | {10.00, 1.89836}, |
---|
| 562 | {12.00, 2.13625}, |
---|
| 563 | {14.00, 2.33787}, |
---|
| 564 | {16.00, 2.51093}, |
---|
| 565 | {18.00, 2.66134}, |
---|
| 566 | {20.00, 2.79358}, |
---|
| 567 | {25.00, 3.06539}, |
---|
| 568 | {30.00, 3.27902}, |
---|
| 569 | {35.00, 3.45430}, |
---|
| 570 | {40.00, 3.60281}, |
---|
| 571 | {45.00, 3.73167}, |
---|
| 572 | {50.00, 3.84555}, |
---|
| 573 | {60.00, 4.04011}, |
---|
| 574 | {70.00, 4.20264}, |
---|
| 575 | {80.00, 4.34229}, |
---|
| 576 | {90.00, 4.46474}, |
---|
| 577 | {100.00, 4.57378}, |
---|
| 578 | {120.00, 4.76155}, |
---|
| 579 | {140.00, 4.91953}, |
---|
| 580 | {160.00, 5.05590}, |
---|
| 581 | {180.00, 5.17588}, |
---|
| 582 | {200.00, 5.28299}, |
---|
| 583 | {250.00, 5.50925}, |
---|
| 584 | {300.00, 5.69364}, |
---|
| 585 | {350.00, 5.84926}, |
---|
| 586 | {400.00, 5.98388}, |
---|
| 587 | {450.00, 6.10252}, |
---|
| 588 | {500.00, 6.20856}, |
---|
| 589 | {600.00, 6.39189}, |
---|
| 590 | {700.00, 6.54677}, |
---|
| 591 | {800.00, 6.68084}, |
---|
| 592 | {900.00, 6.79905}, |
---|
| 593 | {1000.00, 6.90474} |
---|
| 594 | }; |
---|
| 595 | |
---|
| 596 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 597 | |
---|
| 598 | // Data for L1 from: Mikkelsen H.H., Sigmund P. Phys. Rev. A40 (1989) 101-116 |
---|
| 599 | const G4double G4ICRU73QOModel::L1[22][2] = |
---|
| 600 | { |
---|
| 601 | {0.00, -0.000001}, |
---|
| 602 | {0.10, -0.00001}, |
---|
| 603 | {0.20, -0.00049}, |
---|
| 604 | {0.30, -0.00084}, |
---|
| 605 | {0.40, 0.00085}, |
---|
| 606 | {0.50, 0.00519}, |
---|
| 607 | {0.60, 0.01198}, |
---|
| 608 | {0.70, 0.02074}, |
---|
| 609 | {0.80, 0.03133}, |
---|
| 610 | {0.90, 0.04369}, |
---|
| 611 | {1.00, 0.06035}, |
---|
| 612 | {2.00, 0.24023}, |
---|
| 613 | {3.00, 0.44284}, |
---|
| 614 | {4.00, 0.62012}, |
---|
| 615 | {5.00, 0.77031}, |
---|
| 616 | {6.00, 0.90390}, |
---|
| 617 | {7.00, 1.02705}, |
---|
| 618 | {8.00, 1.10867}, |
---|
| 619 | {9.00, 1.17546}, |
---|
| 620 | {10.00, 1.21599}, |
---|
| 621 | {15.00, 1.24349}, |
---|
| 622 | {20.00, 1.16752} |
---|
| 623 | }; |
---|
| 624 | |
---|
| 625 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 626 | |
---|
| 627 | // Data for L2 from: Mikkelsen H.H. Nucl. Instr. Meth. B58 (1991) 136-148 |
---|
| 628 | const G4double G4ICRU73QOModel::L2[14][2] = |
---|
| 629 | { |
---|
| 630 | {0.00, 0.000001}, |
---|
| 631 | {0.10, 0.00001}, |
---|
| 632 | {0.20, 0.00000}, |
---|
| 633 | {0.40, -0.00120}, |
---|
| 634 | {0.60, -0.00036}, |
---|
| 635 | {0.80, 0.00372}, |
---|
| 636 | {1.00, 0.01298}, |
---|
| 637 | {2.00, 0.08296}, |
---|
| 638 | {4.00, 0.21953}, |
---|
| 639 | {6.00, 0.23903}, |
---|
| 640 | {8.00, 0.20893}, |
---|
| 641 | {10.00, 0.10879}, |
---|
| 642 | {20.00, -0.88409}, |
---|
| 643 | {40.00, -1.13902} |
---|
| 644 | }; |
---|
| 645 | |
---|
| 646 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 647 | |
---|
| 648 | // Correction obtained by V.Ivanchenko using G4BetheBlochModel |
---|
| 649 | const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0, |
---|
| 650 | 0.9637, 0.9877, 0.9467, 0.9856, 0.9066, 0.9886, 0.9256, 0.9822, 0.8729, 0.9965, // 1 - 10 |
---|
| 651 | 0.8564, 0.8879, 0.9577, 0.9606, 0.9104, 0.8962, 0.9316, 0.9378, 0.9189, 0.9543, // 11 - 20 |
---|
| 652 | 0.9114, 0.9846, 0.8497, 0.8311, 0.8081, 0.985, 0.76, 0.9696, 0.9964, 0.7361, // 21 - 30 |
---|
| 653 | 0.7787, 1.066, 0.8547, 0.8759, 0.9135, 1.042, 0.9256, 0.9555, 0.965, 0.9619, // 31 - 40 |
---|
| 654 | 0.9361, 0.9078, 0.9291, 0.9128, 0.9115, 0.876, 0.8701, 0.7955, 0.8197, 0.8551, // 41 - 50 |
---|
| 655 | 0.8318, 0.8646, 0.8802, 0.8871, 0.8881, 0.9103, 0.9113, 0.8879, 0.8588, 0.8279, // 51 - 60 |
---|
| 656 | 0.7953, 0.7685, 0.7433, 0.7219, 0.7003, 0.6732, 0.6485, 0.63, 0.6232, 0.611, // 61 - 70 |
---|
| 657 | 0.6385, 0.6642, 0.8452, 0.8127, 0.6982, 0.6971, 0.6889, 0.8416, 0.8341, 0.6956, // 71 - 80 |
---|
| 658 | 0.7308, 0.7777, 0.7692, 0.7924, 0.8054, 0.8203, 0.8386, 0.8618, 0.8654, 0.8558, // 81 - 90 |
---|
| 659 | 0.8323, 0.8362, 0.8044, 0.7811, 0.7685, 0.7714, 0.759, 0.7506}; |
---|