physics/src/G4LowEnergyCompton.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 // Author: A. Forti
0028 //         Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
0029 //
0030 // History:
0031 // --------
0032 // Added Livermore data table construction methods A. Forti
0033 // Modified BuildMeanFreePath to read new data tables A. Forti
0034 // Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti
0035 // Added SelectRandomAtom A. Forti
0036 // Added map of the elements A. Forti
0037 // 24.04.2001 V.Ivanchenko - Remove RogueWave
0038 // 06.08.2001 MGP          - Revised according to a design iteration
0039 // 22.01.2003 V.Ivanchenko - Cut per region
0040 // 10.03.2003 V.Ivanchenko - Remove CutPerMaterial warning
0041 // 24.04.2003 V.Ivanchenko - Cut per region mfpt
0042 //
0043 // -------------------------------------------------------------------
0044
0045 #include "G4LowEnergyCompton.hh"
0046 #include "Randomize.hh"
0047 #include "G4PhysicalConstants.hh"
0048 #include "G4SystemOfUnits.hh"
0049 #include "G4ParticleDefinition.hh"
0050 #include "G4Track.hh"
0051 #include "G4Step.hh"
0052 #include "G4ForceCondition.hh"
0053 #include "G4Gamma.hh"
0054 #include "G4Electron.hh"
0055 #include "G4DynamicParticle.hh"
0056 #include "G4VParticleChange.hh"
0057 #include "G4ThreeVector.hh"
0058 #include "G4EnergyLossTables.hh"
0059 #include "G4RDVCrossSectionHandler.hh"
0060 #include "G4RDCrossSectionHandler.hh"
0061 #include "G4RDVEMDataSet.hh"
0062 #include "G4RDCompositeEMDataSet.hh"
0063 #include "G4RDVDataSetAlgorithm.hh"
0064 #include "G4RDLogLogInterpolation.hh"
0065 #include "G4RDVRangeTest.hh"
0066 #include "G4RDRangeTest.hh"
0067 #include "G4RDRangeNoTest.hh"
0068 #include "G4MaterialCutsCouple.hh"
0069
0070
0071 G4LowEnergyCompton::G4LowEnergyCompton(const G4String& processName)
0072   : G4VDiscreteProcess(processName),
0073     lowEnergyLimit(250*eV),
0074     highEnergyLimit(100*GeV),
0075     intrinsicLowEnergyLimit(10*eV),
0076     intrinsicHighEnergyLimit(100*GeV)
0077 {
0078   if (lowEnergyLimit < intrinsicLowEnergyLimit ||
0079       highEnergyLimit > intrinsicHighEnergyLimit)
0080     {
0081       G4Exception("G4LowEnergyCompton::G4LowEnergyCompton()",
0082                   "OutOfRange", FatalException,
0083                   "Energy outside intrinsic process validity range!");
0084     }
0085
0086   crossSectionHandler = new G4RDCrossSectionHandler;
0087
0088   G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
0089   G4String scatterFile = "comp/ce-sf-";
0090   scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation, 1., 1.);
0091   scatterFunctionData->LoadData(scatterFile);
0092
0093   meanFreePathTable = 0;
0094
0095   rangeTest = new G4RDRangeNoTest;
0096
0097   // For Doppler broadening
0098   shellData.SetOccupancyData();
0099
0100    if (verboseLevel > 0)
0101      {
0102        G4cout << GetProcessName() << " is created " << G4endl
0103           << "Energy range: "
0104           << lowEnergyLimit / keV << " keV - "
0105           << highEnergyLimit / GeV << " GeV"
0106           << G4endl;
0107      }
0108 }
0109
0110 G4LowEnergyCompton::~G4LowEnergyCompton()
0111 {
0112   delete meanFreePathTable;
0113   delete crossSectionHandler;
0114   delete scatterFunctionData;
0115   delete rangeTest;
0116 }
0117
0118 void G4LowEnergyCompton::BuildPhysicsTable(const G4ParticleDefinition& )
0119 {
0120
0121   crossSectionHandler->Clear();
0122   G4String crossSectionFile = "comp/ce-cs-";
0123   crossSectionHandler->LoadData(crossSectionFile);
0124
0125   delete meanFreePathTable;
0126   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
0127
0128   // For Doppler broadening
0129   G4String file = "/doppler/shell-doppler";
0130   shellData.LoadData(file);
0131 }
0132
0133 G4VParticleChange* G4LowEnergyCompton::PostStepDoIt(const G4Track& aTrack,
0134                             const G4Step& aStep)
0135 {
0136   // The scattered gamma energy is sampled according to Klein - Nishina formula.
0137   // then accepted or rejected depending on the Scattering Function multiplied
0138   // by factor from Klein - Nishina formula.
0139   // Expression of the angular distribution as Klein Nishina
0140   // angular and energy distribution and Scattering fuctions is taken from
0141   // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
0142   // Phys. Res. B 101 (1995). Method of sampling with form factors is different
0143   // data are interpolated while in the article they are fitted.
0144   // Reference to the article is from J. Stepanek New Photon, Positron
0145   // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
0146   // TeV (draft).
0147   // The random number techniques of Butcher & Messel are used
0148   // (Nucl Phys 20(1960),15).
0149
0150   aParticleChange.Initialize(aTrack);
0151
0152   // Dynamic particle quantities
0153   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
0154   G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
0155
0156   if (photonEnergy0 <= lowEnergyLimit)
0157     {
0158       aParticleChange.ProposeTrackStatus(fStopAndKill);
0159       aParticleChange.ProposeEnergy(0.);
0160       aParticleChange.ProposeLocalEnergyDeposit(photonEnergy0);
0161       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
0162     }
0163
0164   G4double e0m = photonEnergy0 / electron_mass_c2 ;
0165   G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
0166   G4double epsilon0 = 1. / (1. + 2. * e0m);
0167   G4double epsilon0Sq = epsilon0 * epsilon0;
0168   G4double alpha1 = -std::log(epsilon0);
0169   G4double alpha2 = 0.5 * (1. - epsilon0Sq);
0170   G4double wlPhoton = h_Planck*c_light/photonEnergy0;
0171
0172   // Select randomly one element in the current material
0173   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
0174   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
0175
0176   // Sample the energy of the scattered photon
0177   G4double epsilon;
0178   G4double epsilonSq;
0179   G4double oneCosT;
0180   G4double sinT2;
0181   G4double gReject;
0182   do
0183     {
0184       if ( alpha1/(alpha1+alpha2) > G4UniformRand())
0185     {
0186       epsilon = std::exp(-alpha1 * G4UniformRand());  // std::pow(epsilon0,G4UniformRand())
0187       epsilonSq = epsilon * epsilon;
0188     }
0189       else
0190     {
0191       epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
0192       epsilon = std::sqrt(epsilonSq);
0193     }
0194
0195       oneCosT = (1. - epsilon) / ( epsilon * e0m);
0196       sinT2 = oneCosT * (2. - oneCosT);
0197       G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
0198       G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
0199       gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
0200
0201     }  while(gReject < G4UniformRand()*Z);
0202
0203   G4double cosTheta = 1. - oneCosT;
0204   G4double sinTheta = std::sqrt (sinT2);
0205   G4double phi = twopi * G4UniformRand() ;
0206   G4double dirX = sinTheta * std::cos(phi);
0207   G4double dirY = sinTheta * std::sin(phi);
0208   G4double dirZ = cosTheta ;
0209
0210   // Doppler broadening -  Method based on:
0211   // Y. Namito, S. Ban and H. Hirayama,
0212   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
0213   // NIM A 349, pp. 489-494, 1994
0214
0215   // Maximum number of sampling iterations
0216   G4int maxDopplerIterations = 1000;
0217   G4double bindingE = 0.;
0218   G4double photonEoriginal = epsilon * photonEnergy0;
0219   G4double photonE = -1.;
0220   G4int iteration = 0;
0221   G4double eMax = photonEnergy0;
0222   do
0223     {
0224       iteration++;
0225       // Select shell based on shell occupancy
0226       G4int shell = shellData.SelectRandomShell(Z);
0227       bindingE = shellData.BindingEnergy(Z,shell);
0228
0229       eMax = photonEnergy0 - bindingE;
0230
0231       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
0232       G4double pSample = profileData.RandomSelectMomentum(Z,shell);
0233       // Rescale from atomic units
0234       G4double pDoppler = pSample * fine_structure_const;
0235       G4double pDoppler2 = pDoppler * pDoppler;
0236       G4double var2 = 1. + oneCosT * e0m;
0237       G4double var3 = var2*var2 - pDoppler2;
0238       G4double var4 = var2 - pDoppler2 * cosTheta;
0239       G4double var = var4*var4 - var3 + pDoppler2 * var3;
0240       if (var > 0.)
0241     {
0242       G4double varSqrt = std::sqrt(var);
0243       G4double scale = photonEnergy0 / var3;
0244           // Random select either root
0245       if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
0246       else photonE = (var4 + varSqrt) * scale;
0247     }
0248       else
0249     {
0250       photonE = -1.;
0251     }
0252    } while ( iteration <= maxDopplerIterations &&
0253          (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
0254
0255   // End of recalculation of photon energy with Doppler broadening
0256   // Revert to original if maximum number of iterations threshold has been reached
0257   if (iteration >= maxDopplerIterations)
0258     {
0259       photonE = photonEoriginal;
0260       bindingE = 0.;
0261     }
0262
0263   // Update G4VParticleChange for the scattered photon
0264
0265   G4ThreeVector photonDirection1(dirX,dirY,dirZ);
0266   photonDirection1.rotateUz(photonDirection0);
0267   aParticleChange.ProposeMomentumDirection(photonDirection1);
0268
0269   G4double photonEnergy1 = photonE;
0270   //G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
0271
0272   if (photonEnergy1 > 0.)
0273     {
0274       aParticleChange.ProposeEnergy(photonEnergy1) ;
0275     }
0276   else
0277     {
0278       aParticleChange.ProposeEnergy(0.) ;
0279       aParticleChange.ProposeTrackStatus(fStopAndKill);
0280     }
0281
0282   // Kinematics of the scattered electron
0283   G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
0284   G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
0285
0286   G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2;
0287   G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2;
0288   G4double sinThetaE = -1.;
0289   G4double cosThetaE = 0.;
0290   if (electronP2 > 0.)
0291     {
0292       cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2);
0293       sinThetaE = -1. * std::sqrt(1. - cosThetaE * cosThetaE);
0294     }
0295
0296   G4double eDirX = sinThetaE * std::cos(phi);
0297   G4double eDirY = sinThetaE * std::sin(phi);
0298   G4double eDirZ = cosThetaE;
0299
0300   // Generate the electron only if with large enough range w.r.t. cuts and safety
0301
0302   G4double safety = aStep.GetPostStepPoint()->GetSafety();
0303
0304   if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
0305     {
0306       G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
0307       eDirection.rotateUz(photonDirection0);
0308
0309       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),eDirection,eKineticEnergy) ;
0310       aParticleChange.SetNumberOfSecondaries(1);
0311       aParticleChange.AddSecondary(electron);
0312       // Binding energy deposited locally
0313       aParticleChange.ProposeLocalEnergyDeposit(bindingE);
0314     }
0315   else
0316     {
0317       aParticleChange.SetNumberOfSecondaries(0);
0318       aParticleChange.ProposeLocalEnergyDeposit(eKineticEnergy + bindingE);
0319     }
0320
0321   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
0322 }
0323
0324 G4bool G4LowEnergyCompton::IsApplicable(const G4ParticleDefinition& particle)
0325 {
0326   return ( &particle == G4Gamma::Gamma() );
0327 }
0328
0329 G4double G4LowEnergyCompton::GetMeanFreePath(const G4Track& track,
0330                          G4double, // previousStepSize
0331                          G4ForceCondition*)
0332 {
0333   const G4DynamicParticle* photon = track.GetDynamicParticle();
0334   G4double energy = photon->GetKineticEnergy();
0335   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
0336   size_t materialIndex = couple->GetIndex();
0337
0338   G4double meanFreePath;
0339   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
0340   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
0341   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
0342   return meanFreePath;
0343 }