physics/src/G4LowEnergyGammaConversion.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // --------------------------------------------------------------------
0027 ///
0028 //
0029 //
0030 // --------------------------------------------------------------
0031 //
0032 // Author: A. Forti
0033 //         Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
0034 //
0035 // History:
0036 // --------
0037 // 02/03/1999 A. Forti 1st implementation
0038 // 14.03.2000 Veronique Lefebure;
0039 // Change initialisation of lowestEnergyLimit from 1.22 to 1.022.
0040 // Note that the hard coded value 1.022 should be used instead of
0041 // 2*electron_mass_c2 in order to agree with the value of the data bank EPDL97
0042 // 24.04.01 V.Ivanchenko remove RogueWave
0043 // 27.07.01 F.Longo correct bug in energy distribution
0044 // 21.01.03 V.Ivanchenko Cut per region
0045 // 25.03.03 F.Longo fix in angular distribution of e+/e-
0046 // 24.04.03 V.Ivanchenko - Cut per region mfpt
0047 //
0048 // --------------------------------------------------------------
0049
0050 #include "G4LowEnergyGammaConversion.hh"
0051
0052 #include "Randomize.hh"
0053 #include "G4PhysicalConstants.hh"
0054 #include "G4SystemOfUnits.hh"
0055 #include "G4ParticleDefinition.hh"
0056 #include "G4Track.hh"
0057 #include "G4Step.hh"
0058 #include "G4ForceCondition.hh"
0059 #include "G4Gamma.hh"
0060 #include "G4Electron.hh"
0061 #include "G4DynamicParticle.hh"
0062 #include "G4VParticleChange.hh"
0063 #include "G4ThreeVector.hh"
0064 #include "G4Positron.hh"
0065 #include "G4IonisParamElm.hh"
0066 #include "G4Material.hh"
0067 #include "G4RDVCrossSectionHandler.hh"
0068 #include "G4RDCrossSectionHandler.hh"
0069 #include "G4RDVEMDataSet.hh"
0070 #include "G4RDVDataSetAlgorithm.hh"
0071 #include "G4RDLogLogInterpolation.hh"
0072 #include "G4RDVRangeTest.hh"
0073 #include "G4RDRangeTest.hh"
0074 #include "G4MaterialCutsCouple.hh"
0075
0076 G4LowEnergyGammaConversion::G4LowEnergyGammaConversion(const G4String& processName)
0077   : G4VDiscreteProcess(processName),
0078     lowEnergyLimit(1.022000*MeV),
0079     highEnergyLimit(100*GeV),
0080     intrinsicLowEnergyLimit(1.022000*MeV),
0081     intrinsicHighEnergyLimit(100*GeV),
0082     smallEnergy(2.*MeV)
0083
0084 {
0085   if (lowEnergyLimit < intrinsicLowEnergyLimit ||
0086       highEnergyLimit > intrinsicHighEnergyLimit)
0087     {
0088       G4Exception("G4LowEnergyGammaConversion::G4LowEnergyGammaConversion()",
0089                   "OutOfRange", FatalException,
0090                   "Energy limit outside intrinsic process validity range!");
0091     }
0092
0093   // The following pointer is owned by G4DataHandler
0094
0095   crossSectionHandler = new G4RDCrossSectionHandler();
0096   crossSectionHandler->Initialise(0,1.0220*MeV,100.*GeV,400);
0097   meanFreePathTable = 0;
0098   rangeTest = new G4RDRangeTest;
0099
0100    if (verboseLevel > 0)
0101      {
0102        G4cout << GetProcessName() << " is created " << G4endl
0103           << "Energy range: "
0104           << lowEnergyLimit / MeV << " MeV - "
0105           << highEnergyLimit / GeV << " GeV"
0106           << G4endl;
0107      }
0108 }
0109
0110 G4LowEnergyGammaConversion::~G4LowEnergyGammaConversion()
0111 {
0112   delete meanFreePathTable;
0113   delete crossSectionHandler;
0114   delete rangeTest;
0115 }
0116
0117 void G4LowEnergyGammaConversion::BuildPhysicsTable(const G4ParticleDefinition& )
0118 {
0119
0120   crossSectionHandler->Clear();
0121   G4String crossSectionFile = "pair/pp-cs-";
0122   crossSectionHandler->LoadData(crossSectionFile);
0123
0124   delete meanFreePathTable;
0125   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
0126 }
0127
0128 G4VParticleChange* G4LowEnergyGammaConversion::PostStepDoIt(const G4Track& aTrack,
0129                                 const G4Step& aStep)
0130 {
0131 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
0132 // cross sections with Coulomb correction. A modified version of the random
0133 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
0134
0135 // Note 1 : Effects due to the breakdown of the Born approximation at low
0136 // energy are ignored.
0137 // Note 2 : The differential cross section implicitly takes account of
0138 // pair creation in both nuclear and atomic electron fields. However triplet
0139 // prodution is not generated.
0140
0141   aParticleChange.Initialize(aTrack);
0142
0143   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
0144
0145   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
0146   G4double photonEnergy = incidentPhoton->GetKineticEnergy();
0147   G4ParticleMomentum photonDirection = incidentPhoton->GetMomentumDirection();
0148
0149   G4double epsilon ;
0150   G4double epsilon0 = electron_mass_c2 / photonEnergy ;
0151
0152   // Do it fast if photon energy < 2. MeV
0153   if (photonEnergy < smallEnergy )
0154     {
0155       epsilon = epsilon0 + (0.5 - epsilon0) * G4UniformRand();
0156     }
0157   else
0158     {
0159       // Select randomly one element in the current material
0160       const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
0161
0162       if (element == 0)
0163     {
0164       G4cout << "G4LowEnergyGammaConversion::PostStepDoIt - element = 0" << G4endl;
0165     }
0166       G4IonisParamElm* ionisation = element->GetIonisation();
0167        if (ionisation == 0)
0168     {
0169       G4cout << "G4LowEnergyGammaConversion::PostStepDoIt - ionisation = 0" << G4endl;
0170     }
0171
0172       // Extract Coulomb factor for this Element
0173       G4double fZ = 8. * (ionisation->GetlogZ3());
0174       if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
0175
0176       // Limits of the screening variable
0177       G4double screenFactor = 136. * epsilon0 / (element->GetIonisation()->GetZ3()) ;
0178       G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
0179       G4double screenMin = std::min(4.*screenFactor,screenMax) ;
0180
0181       // Limits of the energy sampling
0182       G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
0183       G4double epsilonMin = std::max(epsilon0,epsilon1);
0184       G4double epsilonRange = 0.5 - epsilonMin ;
0185
0186       // Sample the energy rate of the created electron (or positron)
0187       G4double screen;
0188       G4double gReject ;
0189
0190       G4double f10 = ScreenFunction1(screenMin) - fZ;
0191       G4double f20 = ScreenFunction2(screenMin) - fZ;
0192       G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
0193       G4double normF2 = std::max(1.5 * f20,0.);
0194
0195       do {
0196     if (normF1 / (normF1 + normF2) > G4UniformRand() )
0197       {
0198         epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.3333) ;
0199         screen = screenFactor / (epsilon * (1. - epsilon));
0200         gReject = (ScreenFunction1(screen) - fZ) / f10 ;
0201       }
0202     else
0203       {
0204         epsilon = epsilonMin + epsilonRange * G4UniformRand();
0205         screen = screenFactor / (epsilon * (1 - epsilon));
0206         gReject = (ScreenFunction2(screen) - fZ) / f20 ;
0207       }
0208       } while ( gReject < G4UniformRand() );
0209
0210     }   //  End of epsilon sampling
0211
0212   // Fix charges randomly
0213
0214   G4double electronTotEnergy;
0215   G4double positronTotEnergy;
0216
0217   if (CLHEP::RandBit::shootBit())
0218     {
0219       electronTotEnergy = (1. - epsilon) * photonEnergy;
0220       positronTotEnergy = epsilon * photonEnergy;
0221     }
0222   else
0223     {
0224       positronTotEnergy = (1. - epsilon) * photonEnergy;
0225       electronTotEnergy = epsilon * photonEnergy;
0226     }
0227
0228   // Scattered electron (positron) angles. ( Z - axis along the parent photon)
0229   // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
0230   // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
0231
0232   G4double u;
0233   const G4double a1 = 0.625;
0234   G4double a2 = 3. * a1;
0235   //  G4double d = 27. ;
0236
0237   //  if (9. / (9. + d) > G4UniformRand())
0238   if (0.25 > G4UniformRand())
0239     {
0240       u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
0241     }
0242   else
0243     {
0244       u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
0245     }
0246
0247   G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
0248   G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
0249   G4double phi  = twopi * G4UniformRand();
0250
0251   G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
0252   G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
0253
0254
0255   // Kinematics of the created pair:
0256   // the electron and positron are assumed to have a symetric angular
0257   // distribution with respect to the Z axis along the parent photon
0258
0259   G4double localEnergyDeposit = 0. ;
0260
0261   aParticleChange.SetNumberOfSecondaries(2) ;
0262   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
0263
0264   // Generate the electron only if with large enough range w.r.t. cuts and safety
0265
0266   G4double safety = aStep.GetPostStepPoint()->GetSafety();
0267
0268   if (rangeTest->Escape(G4Electron::Electron(),couple,electronKineEnergy,safety))
0269     {
0270       G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
0271       electronDirection.rotateUz(photonDirection);
0272
0273       G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
0274                                 electronDirection,
0275                                 electronKineEnergy);
0276       aParticleChange.AddSecondary(particle1) ;
0277     }
0278   else
0279     {
0280       localEnergyDeposit += electronKineEnergy ;
0281     }
0282
0283   // The e+ is always created (even with kinetic energy = 0) for further annihilation
0284   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
0285
0286   // Is the local energy deposit correct, if the positron is always created?
0287   if (! (rangeTest->Escape(G4Positron::Positron(),couple,positronKineEnergy,safety)))
0288     {
0289       localEnergyDeposit += positronKineEnergy ;
0290       positronKineEnergy = 0. ;
0291     }
0292
0293   G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
0294   positronDirection.rotateUz(photonDirection);
0295
0296   // Create G4DynamicParticle object for the particle2
0297   G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
0298                                positronDirection, positronKineEnergy);
0299   aParticleChange.AddSecondary(particle2) ;
0300
0301   aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit) ;
0302
0303   // Kill the incident photon
0304   aParticleChange.ProposeMomentumDirection(0.,0.,0.) ;
0305   aParticleChange.ProposeEnergy(0.) ;
0306   aParticleChange.ProposeTrackStatus(fStopAndKill) ;
0307
0308   //  Reset NbOfInteractionLengthLeft and return aParticleChange
0309   return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
0310 }
0311
0312 G4bool G4LowEnergyGammaConversion::IsApplicable(const G4ParticleDefinition& particle)
0313 {
0314   return ( &particle == G4Gamma::Gamma() );
0315 }
0316
0317 G4double G4LowEnergyGammaConversion::GetMeanFreePath(const G4Track& track,
0318                              G4double, // previousStepSize
0319                              G4ForceCondition*)
0320 {
0321   const G4DynamicParticle* photon = track.GetDynamicParticle();
0322   G4double energy = photon->GetKineticEnergy();
0323   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
0324   size_t materialIndex = couple->GetIndex();
0325
0326   G4double meanFreePath;
0327   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
0328   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
0329   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
0330   return meanFreePath;
0331 }
0332
0333 G4double G4LowEnergyGammaConversion::ScreenFunction1(G4double screenVariable)
0334 {
0335   // Compute the value of the screening function 3*phi1 - phi2
0336
0337   G4double value;
0338
0339   if (screenVariable > 1.)
0340     value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
0341   else
0342     value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
0343
0344   return value;
0345 }
0346
0347 G4double G4LowEnergyGammaConversion::ScreenFunction2(G4double screenVariable)
0348 {
0349   // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
0350
0351   G4double value;
0352
0353   if (screenVariable > 1.)
0354     value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
0355   else
0356     value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
0357
0358   return value;
0359 }