00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "G4PenelopeIonisationModel.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4ParticleDefinition.hh"
00050 #include "G4MaterialCutsCouple.hh"
00051 #include "G4ProductionCutsTable.hh"
00052 #include "G4DynamicParticle.hh"
00053 #include "G4AtomicTransitionManager.hh"
00054 #include "G4AtomicShell.hh"
00055 #include "G4Gamma.hh"
00056 #include "G4Electron.hh"
00057 #include "G4Positron.hh"
00058 #include "G4AtomicDeexcitation.hh"
00059 #include "G4PenelopeOscillatorManager.hh"
00060 #include "G4PenelopeOscillator.hh"
00061 #include "G4PenelopeCrossSection.hh"
00062 #include "G4PhysicsFreeVector.hh"
00063 #include "G4PhysicsLogVector.hh"
00064 #include "G4LossTableManager.hh"
00065 #include "G4PenelopeIonisationXSHandler.hh"
00066
00067
00068
00069
00070 G4PenelopeIonisationModel::G4PenelopeIonisationModel(const G4ParticleDefinition*,
00071 const G4String& nam)
00072 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
00073 fAtomDeexcitation(0),kineticEnergy1(0.*eV),
00074 cosThetaPrimary(1.0),energySecondary(0.*eV),
00075 cosThetaSecondary(0.0),targetOscillator(-1),
00076 theCrossSectionHandler(0)
00077 {
00078 fIntrinsicLowEnergyLimit = 100.0*eV;
00079 fIntrinsicHighEnergyLimit = 100.0*GeV;
00080
00081 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00082 nBins = 200;
00083
00084 oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
00085
00086 verboseLevel= 0;
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096 SetDeexcitationFlag(true);
00097 }
00098
00099
00100
00101 G4PenelopeIonisationModel::~G4PenelopeIonisationModel()
00102 {
00103 if (theCrossSectionHandler)
00104 delete theCrossSectionHandler;
00105 }
00106
00107
00108
00109 void G4PenelopeIonisationModel::Initialise(const G4ParticleDefinition*,
00110 const G4DataVector&)
00111 {
00112 if (verboseLevel > 3)
00113 G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl;
00114
00115 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00116
00117 if (!fAtomDeexcitation)
00118 {
00119 G4cout << G4endl;
00120 G4cout << "WARNING from G4PenelopeIonisationModel " << G4endl;
00121 G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
00122 G4cout << "any fluorescence/Auger emission." << G4endl;
00123 G4cout << "Please make sure this is intended" << G4endl;
00124 }
00125
00126
00127 nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
00128 nBins = std::max(nBins,(size_t)100);
00129
00130
00131 if (theCrossSectionHandler)
00132 {
00133 delete theCrossSectionHandler;
00134 theCrossSectionHandler = 0;
00135 }
00136 theCrossSectionHandler = new G4PenelopeIonisationXSHandler(nBins);
00137 theCrossSectionHandler->SetVerboseLevel(verboseLevel);
00138
00139 if (verboseLevel > 2) {
00140 G4cout << "Penelope Ionisation model v2008 is initialized " << G4endl
00141 << "Energy range: "
00142 << LowEnergyLimit() / keV << " keV - "
00143 << HighEnergyLimit() / GeV << " GeV. Using "
00144 << nBins << " bins."
00145 << G4endl;
00146 }
00147
00148 if(isInitialised) return;
00149 fParticleChange = GetParticleChangeForLoss();
00150 isInitialised = true;
00151 }
00152
00153
00154
00155 G4double G4PenelopeIonisationModel::CrossSectionPerVolume(const G4Material* material,
00156 const G4ParticleDefinition*
00157 theParticle,
00158 G4double energy,
00159 G4double cutEnergy,
00160 G4double)
00161 {
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 if (verboseLevel > 3)
00179 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" << G4endl;
00180
00181 SetupForMaterial(theParticle, material, energy);
00182
00183 G4double totalCross = 0.0;
00184 G4double crossPerMolecule = 0.;
00185
00186 G4PenelopeCrossSection* theXS =
00187 theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
00188 material,
00189 cutEnergy);
00190
00191 if (theXS)
00192 crossPerMolecule = theXS->GetHardCrossSection(energy);
00193
00194 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
00195 G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
00196
00197 if (verboseLevel > 3)
00198 G4cout << "Material " << material->GetName() << " has " << atPerMol <<
00199 "atoms per molecule" << G4endl;
00200
00201 G4double moleculeDensity = 0.;
00202 if (atPerMol)
00203 moleculeDensity = atomDensity/atPerMol;
00204
00205 G4double crossPerVolume = crossPerMolecule*moleculeDensity;
00206
00207 if (verboseLevel > 2)
00208 {
00209 G4cout << "G4PenelopeIonisationModel " << G4endl;
00210 G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " <<
00211 energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
00212 if (theXS)
00213 totalCross = (theXS->GetTotalCrossSection(energy))*moleculeDensity;
00214 G4cout << "Total free path for ionisation (no threshold) at " <<
00215 energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl;
00216 }
00217 return crossPerVolume;
00218 }
00219
00220
00221
00222
00223
00224
00225 G4double G4PenelopeIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00226 G4double,
00227 G4double,
00228 G4double,
00229 G4double,
00230 G4double)
00231 {
00232 G4cout << "*** G4PenelopeIonisationModel -- WARNING ***" << G4endl;
00233 G4cout << "Penelope Ionisation model v2008 does not calculate cross section _per atom_ " << G4endl;
00234 G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
00235 G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
00236 return 0;
00237 }
00238
00239
00240
00241 G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
00242 const G4ParticleDefinition* theParticle,
00243 G4double kineticEnergy,
00244 G4double cutEnergy)
00245 {
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261 if (verboseLevel > 3)
00262 G4cout << "Calling ComputeDEDX() of G4PenelopeIonisationModel" << G4endl;
00263
00264
00265 G4PenelopeCrossSection* theXS =
00266 theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,material,
00267 cutEnergy);
00268 G4double sPowerPerMolecule = 0.0;
00269 if (theXS)
00270 sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
00271
00272
00273 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
00274 G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
00275
00276 G4double moleculeDensity = 0.;
00277 if (atPerMol)
00278 moleculeDensity = atomDensity/atPerMol;
00279
00280 G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
00281
00282 if (verboseLevel > 2)
00283 {
00284 G4cout << "G4PenelopeIonisationModel " << G4endl;
00285 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
00286 kineticEnergy/keV << " keV = " <<
00287 sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
00288 }
00289 return sPowerPerVolume;
00290 }
00291
00292
00293
00294 G4double G4PenelopeIonisationModel::MinEnergyCut(const G4ParticleDefinition*,
00295 const G4MaterialCutsCouple*)
00296 {
00297 return fIntrinsicLowEnergyLimit;
00298 }
00299
00300
00301
00302 void G4PenelopeIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00303 const G4MaterialCutsCouple* couple,
00304 const G4DynamicParticle* aDynamicParticle,
00305 G4double cutE, G4double)
00306 {
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 if (verboseLevel > 3)
00340 G4cout << "Calling SamplingSecondaries() of G4PenelopeIonisationModel" << G4endl;
00341
00342 G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy();
00343 const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
00344
00345 if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
00346 {
00347 fParticleChange->SetProposedKineticEnergy(0.);
00348 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0);
00349 return ;
00350 }
00351
00352 const G4Material* material = couple->GetMaterial();
00353 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(material);
00354
00355 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
00356
00357
00358
00359 kineticEnergy1=kineticEnergy0;
00360 cosThetaPrimary=1.0;
00361 energySecondary=0.0;
00362 cosThetaSecondary=1.0;
00363 targetOscillator = -1;
00364
00365 if (theParticle == G4Electron::Electron())
00366 SampleFinalStateElectron(material,cutE,kineticEnergy0);
00367 else if (theParticle == G4Positron::Positron())
00368 SampleFinalStatePositron(material,cutE,kineticEnergy0);
00369 else
00370 {
00371 G4ExceptionDescription ed;
00372 ed << "Invalid particle " << theParticle->GetParticleName() << G4endl;
00373 G4Exception("G4PenelopeIonisationModel::SamplingSecondaries()",
00374 "em0001",FatalException,ed);
00375
00376 }
00377 if (energySecondary == 0) return;
00378
00379 if (verboseLevel > 3)
00380 {
00381 G4cout << "G4PenelopeIonisationModel::SamplingSecondaries() for " <<
00382 theParticle->GetParticleName() << G4endl;
00383 G4cout << "Final eKin = " << kineticEnergy1 << " keV" << G4endl;
00384 G4cout << "Final cosTheta = " << cosThetaPrimary << G4endl;
00385 G4cout << "Delta-ray eKin = " << energySecondary << " keV" << G4endl;
00386 G4cout << "Delta-ray cosTheta = " << cosThetaSecondary << G4endl;
00387 G4cout << "Oscillator: " << targetOscillator << G4endl;
00388 }
00389
00390
00391 G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
00392 G4double phiPrimary = twopi * G4UniformRand();
00393 G4double dirx = sint * std::cos(phiPrimary);
00394 G4double diry = sint * std::sin(phiPrimary);
00395 G4double dirz = cosThetaPrimary;
00396
00397 G4ThreeVector electronDirection1(dirx,diry,dirz);
00398 electronDirection1.rotateUz(particleDirection0);
00399
00400 if (kineticEnergy1 > 0)
00401 {
00402 fParticleChange->ProposeMomentumDirection(electronDirection1);
00403 fParticleChange->SetProposedKineticEnergy(kineticEnergy1);
00404 }
00405 else
00406 fParticleChange->SetProposedKineticEnergy(0.);
00407
00408
00409
00410 G4double ionEnergyInPenelopeDatabase =
00411 (*theTable)[targetOscillator]->GetIonisationEnergy();
00412
00413
00414
00415 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
00416 G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();
00417
00418
00419 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
00420 G4double bindingEnergy = 0.*eV;
00421
00422
00423 const G4AtomicShell* shell = 0;
00424
00425 if (Z > 0 && shFlag<30)
00426 {
00427 shell = transitionManager->Shell(Z,shFlag-1);
00428 bindingEnergy = shell->BindingEnergy();
00429
00430 }
00431
00432
00433
00434
00435 energySecondary += ionEnergyInPenelopeDatabase-bindingEnergy;
00436
00437 G4double localEnergyDeposit = bindingEnergy;
00438
00439 G4double energyInFluorescence = 0;
00440 G4double energyInAuger = 0;
00441
00442 if (energySecondary < 0)
00443 {
00444
00445
00446
00447
00448 localEnergyDeposit += energySecondary;
00449 energySecondary = 0.0;
00450 }
00451
00452
00453
00454 if (fAtomDeexcitation && shell)
00455 {
00456 G4int index = couple->GetIndex();
00457 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
00458 {
00459 size_t nBefore = fvect->size();
00460 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
00461 size_t nAfter = fvect->size();
00462
00463 if (nAfter > nBefore)
00464 {
00465 for (size_t j=nBefore;j<nAfter;j++)
00466 {
00467 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
00468 localEnergyDeposit -= itsEnergy;
00469 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
00470 energyInFluorescence += itsEnergy;
00471 else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
00472 energyInAuger += itsEnergy;
00473 }
00474 }
00475 }
00476 }
00477
00478
00479 if (energySecondary > cutE)
00480 {
00481 G4DynamicParticle* electron = 0;
00482 G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary);
00483 G4double phiEl = phiPrimary+pi;
00484 G4double xEl = sinThetaE * std::cos(phiEl);
00485 G4double yEl = sinThetaE * std::sin(phiEl);
00486 G4double zEl = cosThetaSecondary;
00487 G4ThreeVector eDirection(xEl,yEl,zEl);
00488 eDirection.rotateUz(particleDirection0);
00489 electron = new G4DynamicParticle (G4Electron::Electron(),
00490 eDirection,energySecondary) ;
00491 fvect->push_back(electron);
00492 }
00493 else
00494 {
00495 localEnergyDeposit += energySecondary;
00496 energySecondary = 0;
00497 }
00498
00499 if (localEnergyDeposit < 0)
00500 {
00501 G4cout << "WARNING-"
00502 << "G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
00503 << G4endl;
00504 localEnergyDeposit=0.;
00505 }
00506 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
00507
00508 if (verboseLevel > 1)
00509 {
00510 G4cout << "-----------------------------------------------------------" << G4endl;
00511 G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
00512 G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
00513 G4cout << "-----------------------------------------------------------" << G4endl;
00514 G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
00515 G4cout << "Delta ray " << energySecondary/keV << " keV" << G4endl;
00516 if (energyInFluorescence)
00517 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
00518 if (energyInAuger)
00519 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
00520 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
00521 G4cout << "Total final state: " << (energySecondary+energyInFluorescence+kineticEnergy1+
00522 localEnergyDeposit+energyInAuger)/keV <<
00523 " keV" << G4endl;
00524 G4cout << "-----------------------------------------------------------" << G4endl;
00525 }
00526
00527 if (verboseLevel > 0)
00528 {
00529 G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+
00530 localEnergyDeposit+energyInAuger-kineticEnergy0);
00531 if (energyDiff > 0.05*keV)
00532 G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " <<
00533 (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit+energyInAuger)/keV <<
00534 " keV (final) vs. " <<
00535 kineticEnergy0/keV << " keV (initial)" << G4endl;
00536 }
00537
00538 }
00539
00540
00541 void G4PenelopeIonisationModel::SampleFinalStateElectron(const G4Material* mat,
00542 G4double cutEnergy,
00543 G4double kineticEnergy)
00544 {
00545
00546
00547
00548
00549
00550
00551
00552
00553 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
00554 size_t numberOfOscillators = theTable->size();
00555 G4PenelopeCrossSection* theXS =
00556 theCrossSectionHandler->GetCrossSectionTableForCouple(G4Electron::Electron(),mat,
00557 cutEnergy);
00558 G4double delta = theCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);
00559
00560
00561 G4double TST = G4UniformRand();
00562 targetOscillator = numberOfOscillators-1;
00563 G4double XSsum = 0.;
00564
00565 for (size_t i=0;i<numberOfOscillators-1;i++)
00566 {
00567
00568
00569
00570
00571
00572
00573 XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
00574
00575 if (XSsum > TST)
00576 {
00577 targetOscillator = (G4int) i;
00578 break;
00579 }
00580 }
00581
00582
00583 if (verboseLevel > 3)
00584 {
00585 G4cout << "SampleFinalStateElectron: sampled oscillator #" << targetOscillator << "." << G4endl;
00586 G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV <<
00587 " eV " << G4endl;
00588 G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV << " eV "
00589 << G4endl;
00590 }
00591
00592 G4double rb = kineticEnergy + 2.0*electron_mass_c2;
00593 G4double gam = 1.0+kineticEnergy/electron_mass_c2;
00594 G4double gam2 = gam*gam;
00595 G4double beta2 = (gam2-1.0)/gam2;
00596 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
00597
00598
00599 G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
00600 G4double invResEne = 1.0/resEne;
00601 G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
00602 G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
00603 G4double XHDL = 0.;
00604 G4double XHDT = 0.;
00605 G4double QM = 0.;
00606 G4double cps = 0.;
00607 G4double cp = 0.;
00608
00609
00610 if (resEne > cutEnergy && resEne < kineticEnergy)
00611 {
00612 cps = kineticEnergy*rb;
00613 cp = std::sqrt(cps);
00614 G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
00615 if (resEne > 1.0e-6*kineticEnergy)
00616 {
00617 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
00618 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
00619 }
00620 else
00621 {
00622 QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
00623 QM *= (1.0-QM*0.5/electron_mass_c2);
00624 }
00625 if (QM < cutoffEne)
00626 {
00627 XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
00628 *invResEne;
00629 XHDT = XHDT0*invResEne;
00630 }
00631 else
00632 {
00633 QM = cutoffEne;
00634 XHDL = 0.;
00635 XHDT = 0.;
00636 }
00637 }
00638 else
00639 {
00640 QM = cutoffEne;
00641 cps = 0.;
00642 cp = 0.;
00643 XHDL = 0.;
00644 XHDT = 0.;
00645 }
00646
00647
00648 G4double EE = kineticEnergy + ionEne;
00649 G4double wmaxc = 0.5*EE;
00650 G4double wcl = std::max(cutEnergy,cutoffEne);
00651 G4double rcl = wcl/EE;
00652 G4double XHC = 0.;
00653 if (wcl < wmaxc)
00654 {
00655 G4double rl1 = 1.0-rcl;
00656 G4double rrl1 = 1.0/rl1;
00657 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
00658 (1.0-amol)*std::log(rcl*rrl1))/EE;
00659 }
00660
00661
00662 G4double XHTOT = XHC + XHDL + XHDT;
00663
00664
00665 if (XHTOT < 1.e-14*barn)
00666 {
00667 kineticEnergy1=kineticEnergy;
00668 cosThetaPrimary=1.0;
00669 energySecondary=0.0;
00670 cosThetaSecondary=1.0;
00671 targetOscillator = numberOfOscillators-1;
00672 return;
00673 }
00674
00675
00676 TST = XHTOT*G4UniformRand();
00677
00678
00679
00680 G4double TS1 = XHC;
00681
00682 if (TST < TS1)
00683 {
00684 G4double A = 5.0*amol;
00685 G4double ARCL = A*0.5*rcl;
00686 G4double rk=0.;
00687 G4bool loopAgain = false;
00688 do
00689 {
00690 loopAgain = false;
00691 G4double fb = (1.0+ARCL)*G4UniformRand();
00692 if (fb < 1)
00693 rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
00694 else
00695 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
00696 G4double rk2 = rk*rk;
00697 G4double rkf = rk/(1.0-rk);
00698 G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
00699 if (G4UniformRand()*(1.0+A*rk2) > phi)
00700 loopAgain = true;
00701 }while(loopAgain);
00702
00703 G4double deltaE = rk*EE;
00704
00705 kineticEnergy1 = kineticEnergy - deltaE;
00706 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
00707
00708 energySecondary = deltaE - ionEne;
00709 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
00710 if (verboseLevel > 3)
00711 G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl;
00712 return;
00713 }
00714
00715
00716 TS1 += XHDL;
00717 G4double deltaE = resEne;
00718 kineticEnergy1 = kineticEnergy - deltaE;
00719
00720 if (TST < TS1)
00721 {
00722 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
00723 G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
00724 - (QS*0.5/electron_mass_c2));
00725 G4double QTREV = Q*(Q+2.0*electron_mass_c2);
00726 G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
00727 cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
00728 if (cosThetaPrimary > 1.)
00729 cosThetaPrimary = 1.0;
00730
00731 energySecondary = deltaE - ionEne;
00732 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
00733 if (cosThetaSecondary > 1.0)
00734 cosThetaSecondary = 1.0;
00735 if (verboseLevel > 3)
00736 G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl;
00737 return;
00738 }
00739
00740
00741 cosThetaPrimary = 1.0;
00742
00743 energySecondary = deltaE - ionEne;
00744 cosThetaSecondary = 0.5;
00745 if (verboseLevel > 3)
00746 G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl;
00747
00748 return;
00749 }
00750
00751
00752
00753
00754 void G4PenelopeIonisationModel::SampleFinalStatePositron(const G4Material* mat,
00755 G4double cutEnergy,
00756 G4double kineticEnergy)
00757 {
00758
00759
00760
00761
00762
00763
00764
00765
00766 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
00767 size_t numberOfOscillators = theTable->size();
00768 G4PenelopeCrossSection* theXS = theCrossSectionHandler->GetCrossSectionTableForCouple(G4Positron::Positron(),mat,
00769 cutEnergy);
00770 G4double delta = theCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);
00771
00772
00773 G4double TST = G4UniformRand();
00774 targetOscillator = numberOfOscillators-1;
00775 G4double XSsum = 0.;
00776 for (size_t i=0;i<numberOfOscillators-1;i++)
00777 {
00778 XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
00779 if (XSsum > TST)
00780 {
00781 targetOscillator = (G4int) i;
00782 break;
00783 }
00784 }
00785
00786 if (verboseLevel > 3)
00787 {
00788 G4cout << "SampleFinalStatePositron: sampled oscillator #" << targetOscillator << "." << G4endl;
00789 G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV
00790 << " eV " << G4endl;
00791 G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV
00792 << " eV " << G4endl;
00793 }
00794
00795
00796 G4double rb = kineticEnergy + 2.0*electron_mass_c2;
00797 G4double gam = 1.0+kineticEnergy/electron_mass_c2;
00798 G4double gam2 = gam*gam;
00799 G4double beta2 = (gam2-1.0)/gam2;
00800 G4double g12 = (gam+1.0)*(gam+1.0);
00801 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
00802
00803 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
00804 G4double bha2 = amol*(3.0+1.0/g12);
00805 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
00806 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
00807
00808
00809
00810
00811 G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
00812 G4double invResEne = 1.0/resEne;
00813 G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
00814 G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
00815
00816 G4double XHDL = 0.;
00817 G4double XHDT = 0.;
00818 G4double QM = 0.;
00819 G4double cps = 0.;
00820 G4double cp = 0.;
00821
00822
00823 if (resEne > cutEnergy && resEne < kineticEnergy)
00824 {
00825 cps = kineticEnergy*rb;
00826 cp = std::sqrt(cps);
00827 G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
00828 if (resEne > 1.0e-6*kineticEnergy)
00829 {
00830 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
00831 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
00832 }
00833 else
00834 {
00835 QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
00836 QM *= (1.0-QM*0.5/electron_mass_c2);
00837 }
00838 if (QM < cutoffEne)
00839 {
00840 XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
00841 *invResEne;
00842 XHDT = XHDT0*invResEne;
00843 }
00844 else
00845 {
00846 QM = cutoffEne;
00847 XHDL = 0.;
00848 XHDT = 0.;
00849 }
00850 }
00851 else
00852 {
00853 QM = cutoffEne;
00854 cps = 0.;
00855 cp = 0.;
00856 XHDL = 0.;
00857 XHDT = 0.;
00858 }
00859
00860 G4double wmaxc = kineticEnergy;
00861 G4double wcl = std::max(cutEnergy,cutoffEne);
00862 G4double rcl = wcl/kineticEnergy;
00863 G4double XHC = 0.;
00864 if (wcl < wmaxc)
00865 {
00866 G4double rl1 = 1.0-rcl;
00867 XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
00868 + (bha3/2.0)*(rcl*rcl-1.0)
00869 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
00870 }
00871
00872
00873 G4double XHTOT = XHC + XHDL + XHDT;
00874
00875
00876 if (XHTOT < 1.e-14*barn)
00877 {
00878 kineticEnergy1=kineticEnergy;
00879 cosThetaPrimary=1.0;
00880 energySecondary=0.0;
00881 cosThetaSecondary=1.0;
00882 targetOscillator = numberOfOscillators-1;
00883 return;
00884 }
00885
00886
00887 TST = XHTOT*G4UniformRand();
00888
00889
00890 G4double TS1 = XHC;
00891 if (TST < TS1)
00892 {
00893 G4double rl1 = 1.0-rcl;
00894 G4double rk=0.;
00895 G4bool loopAgain = false;
00896 do
00897 {
00898 loopAgain = false;
00899 rk = rcl/(1.0-G4UniformRand()*rl1);
00900 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
00901 if (G4UniformRand() > phi)
00902 loopAgain = true;
00903 }while(loopAgain);
00904
00905 G4double deltaE = rk*kineticEnergy;
00906 kineticEnergy1 = kineticEnergy - deltaE;
00907 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
00908
00909 energySecondary = deltaE - ionEne;
00910 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
00911 if (verboseLevel > 3)
00912 G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl;
00913 return;
00914 }
00915
00916
00917 TS1 += XHDL;
00918 G4double deltaE = resEne;
00919 kineticEnergy1 = kineticEnergy - deltaE;
00920 if (TST < TS1)
00921 {
00922 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
00923 G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
00924 - (QS*0.5/electron_mass_c2));
00925 G4double QTREV = Q*(Q+2.0*electron_mass_c2);
00926 G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
00927 cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
00928 if (cosThetaPrimary > 1.)
00929 cosThetaPrimary = 1.0;
00930
00931 energySecondary = deltaE - ionEne;
00932 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
00933 if (cosThetaSecondary > 1.0)
00934 cosThetaSecondary = 1.0;
00935 if (verboseLevel > 3)
00936 G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl;
00937 return;
00938 }
00939
00940
00941 cosThetaPrimary = 1.0;
00942
00943 energySecondary = deltaE - ionEne;
00944 cosThetaSecondary = 0.5;
00945
00946 if (verboseLevel > 3)
00947 G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl;
00948
00949 return;
00950 }
00951