Geant4 Cross Reference |
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 // 27 // Author: A. Forti 28 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch) 29 // 30 // History: 31 // -------- 32 // Added Livermore data table construction methods A. Forti 33 // Modified BuildMeanFreePath to read new data tables A. Forti 34 // Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti 35 // Added SelectRandomAtom A. Forti 36 // Added map of the elements A. Forti 37 // 24.04.2001 V.Ivanchenko - Remove RogueWave 38 // 06.08.2001 MGP - Revised according to a design iteration 39 // 22.01.2003 V.Ivanchenko - Cut per region 40 // 10.03.2003 V.Ivanchenko - Remove CutPerMaterial warning 41 // 24.04.2003 V.Ivanchenko - Cut per region mfpt 42 // 43 // ------------------------------------------------------------------- 44 45 #include "G4LowEnergyCompton.hh" 46 #include "Randomize.hh" 47 #include "G4PhysicalConstants.hh" 48 #include "G4SystemOfUnits.hh" 49 #include "G4ParticleDefinition.hh" 50 #include "G4Track.hh" 51 #include "G4Step.hh" 52 #include "G4ForceCondition.hh" 53 #include "G4Gamma.hh" 54 #include "G4Electron.hh" 55 #include "G4DynamicParticle.hh" 56 #include "G4VParticleChange.hh" 57 #include "G4ThreeVector.hh" 58 #include "G4EnergyLossTables.hh" 59 #include "G4RDVCrossSectionHandler.hh" 60 #include "G4RDCrossSectionHandler.hh" 61 #include "G4RDVEMDataSet.hh" 62 #include "G4RDCompositeEMDataSet.hh" 63 #include "G4RDVDataSetAlgorithm.hh" 64 #include "G4RDLogLogInterpolation.hh" 65 #include "G4RDVRangeTest.hh" 66 #include "G4RDRangeTest.hh" 67 #include "G4RDRangeNoTest.hh" 68 #include "G4MaterialCutsCouple.hh" 69 70 71 G4LowEnergyCompton::G4LowEnergyCompton(const G4String& processName) 72 : G4VDiscreteProcess(processName), 73 lowEnergyLimit(250*eV), 74 highEnergyLimit(100*GeV), 75 intrinsicLowEnergyLimit(10*eV), 76 intrinsicHighEnergyLimit(100*GeV) 77 { 78 if (lowEnergyLimit < intrinsicLowEnergyLimit || 79 highEnergyLimit > intrinsicHighEnergyLimit) 80 { 81 G4Exception("G4LowEnergyCompton::G4LowEnergyCompton()", 82 "OutOfRange", FatalException, 83 "Energy outside intrinsic process validity range!"); 84 } 85 86 crossSectionHandler = new G4RDCrossSectionHandler; 87 88 G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation; 89 G4String scatterFile = "comp/ce-sf-"; 90 scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation, 1., 1.); 91 scatterFunctionData->LoadData(scatterFile); 92 93 meanFreePathTable = 0; 94 95 rangeTest = new G4RDRangeNoTest; 96 97 // For Doppler broadening 98 shellData.SetOccupancyData(); 99 100 if (verboseLevel > 0) 101 { 102 G4cout << GetProcessName() << " is created " << G4endl 103 << "Energy range: " 104 << lowEnergyLimit / keV << " keV - " 105 << highEnergyLimit / GeV << " GeV" 106 << G4endl; 107 } 108 } 109 110 G4LowEnergyCompton::~G4LowEnergyCompton() 111 { 112 delete meanFreePathTable; 113 delete crossSectionHandler; 114 delete scatterFunctionData; 115 delete rangeTest; 116 } 117 118 void G4LowEnergyCompton::BuildPhysicsTable(const G4ParticleDefinition& ) 119 { 120 121 crossSectionHandler->Clear(); 122 G4String crossSectionFile = "comp/ce-cs-"; 123 crossSectionHandler->LoadData(crossSectionFile); 124 125 delete meanFreePathTable; 126 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials(); 127 128 // For Doppler broadening 129 G4String file = "/doppler/shell-doppler"; 130 shellData.LoadData(file); 131 } 132 133 G4VParticleChange* G4LowEnergyCompton::PostStepDoIt(const G4Track& aTrack, 134 const G4Step& aStep) 135 { 136 // The scattered gamma energy is sampled according to Klein - Nishina formula. 137 // then accepted or rejected depending on the Scattering Function multiplied 138 // by factor from Klein - Nishina formula. 139 // Expression of the angular distribution as Klein Nishina 140 // angular and energy distribution and Scattering fuctions is taken from 141 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth. 142 // Phys. Res. B 101 (1995). Method of sampling with form factors is different 143 // data are interpolated while in the article they are fitted. 144 // Reference to the article is from J. Stepanek New Photon, Positron 145 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10 146 // TeV (draft). 147 // The random number techniques of Butcher & Messel are used 148 // (Nucl Phys 20(1960),15). 149 150 aParticleChange.Initialize(aTrack); 151 152 // Dynamic particle quantities 153 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle(); 154 G4double photonEnergy0 = incidentPhoton->GetKineticEnergy(); 155 156 if (photonEnergy0 <= lowEnergyLimit) 157 { 158 aParticleChange.ProposeTrackStatus(fStopAndKill); 159 aParticleChange.ProposeEnergy(0.); 160 aParticleChange.ProposeLocalEnergyDeposit(photonEnergy0); 161 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); 162 } 163 164 G4double e0m = photonEnergy0 / electron_mass_c2 ; 165 G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection(); 166 G4double epsilon0 = 1. / (1. + 2. * e0m); 167 G4double epsilon0Sq = epsilon0 * epsilon0; 168 G4double alpha1 = -std::log(epsilon0); 169 G4double alpha2 = 0.5 * (1. - epsilon0Sq); 170 G4double wlPhoton = h_Planck*c_light/photonEnergy0; 171 172 // Select randomly one element in the current material 173 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple(); 174 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0); 175 176 // Sample the energy of the scattered photon 177 G4double epsilon; 178 G4double epsilonSq; 179 G4double oneCosT; 180 G4double sinT2; 181 G4double gReject; 182 do 183 { 184 if ( alpha1/(alpha1+alpha2) > G4UniformRand()) 185 { 186 epsilon = std::exp(-alpha1 * G4UniformRand()); // std::pow(epsilon0,G4UniformRand()) 187 epsilonSq = epsilon * epsilon; 188 } 189 else 190 { 191 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand(); 192 epsilon = std::sqrt(epsilonSq); 193 } 194 195 oneCosT = (1. - epsilon) / ( epsilon * e0m); 196 sinT2 = oneCosT * (2. - oneCosT); 197 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm); 198 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1); 199 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction; 200 201 } while(gReject < G4UniformRand()*Z); 202 203 G4double cosTheta = 1. - oneCosT; 204 G4double sinTheta = std::sqrt (sinT2); 205 G4double phi = twopi * G4UniformRand() ; 206 G4double dirX = sinTheta * std::cos(phi); 207 G4double dirY = sinTheta * std::sin(phi); 208 G4double dirZ = cosTheta ; 209 210 // Doppler broadening - Method based on: 211 // Y. Namito, S. Ban and H. Hirayama, 212 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 213 // NIM A 349, pp. 489-494, 1994 214 215 // Maximum number of sampling iterations 216 G4int maxDopplerIterations = 1000; 217 G4double bindingE = 0.; 218 G4double photonEoriginal = epsilon * photonEnergy0; 219 G4double photonE = -1.; 220 G4int iteration = 0; 221 G4double eMax = photonEnergy0; 222 do 223 { 224 iteration++; 225 // Select shell based on shell occupancy 226 G4int shell = shellData.SelectRandomShell(Z); 227 bindingE = shellData.BindingEnergy(Z,shell); 228 229 eMax = photonEnergy0 - bindingE; 230 231 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units) 232 G4double pSample = profileData.RandomSelectMomentum(Z,shell); 233 // Rescale from atomic units 234 G4double pDoppler = pSample * fine_structure_const; 235 G4double pDoppler2 = pDoppler * pDoppler; 236 G4double var2 = 1. + oneCosT * e0m; 237 G4double var3 = var2*var2 - pDoppler2; 238 G4double var4 = var2 - pDoppler2 * cosTheta; 239 G4double var = var4*var4 - var3 + pDoppler2 * var3; 240 if (var > 0.) 241 { 242 G4double varSqrt = std::sqrt(var); 243 G4double scale = photonEnergy0 / var3; 244 // Random select either root 245 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale; 246 else photonE = (var4 + varSqrt) * scale; 247 } 248 else 249 { 250 photonE = -1.; 251 } 252 } while ( iteration <= maxDopplerIterations && 253 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) ); 254 255 // End of recalculation of photon energy with Doppler broadening 256 // Revert to original if maximum number of iterations threshold has been reached 257 if (iteration >= maxDopplerIterations) 258 { 259 photonE = photonEoriginal; 260 bindingE = 0.; 261 } 262 263 // Update G4VParticleChange for the scattered photon 264 265 G4ThreeVector photonDirection1(dirX,dirY,dirZ); 266 photonDirection1.rotateUz(photonDirection0); 267 aParticleChange.ProposeMomentumDirection(photonDirection1); 268 269 G4double photonEnergy1 = photonE; 270 //G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl; 271 272 if (photonEnergy1 > 0.) 273 { 274 aParticleChange.ProposeEnergy(photonEnergy1) ; 275 } 276 else 277 { 278 aParticleChange.ProposeEnergy(0.) ; 279 aParticleChange.ProposeTrackStatus(fStopAndKill); 280 } 281 282 // Kinematics of the scattered electron 283 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE; 284 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2; 285 286 G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2; 287 G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2; 288 G4double sinThetaE = -1.; 289 G4double cosThetaE = 0.; 290 if (electronP2 > 0.) 291 { 292 cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2); 293 sinThetaE = -1. * std::sqrt(1. - cosThetaE * cosThetaE); 294 } 295 296 G4double eDirX = sinThetaE * std::cos(phi); 297 G4double eDirY = sinThetaE * std::sin(phi); 298 G4double eDirZ = cosThetaE; 299 300 // Generate the electron only if with large enough range w.r.t. cuts and safety 301 302 G4double safety = aStep.GetPostStepPoint()->GetSafety(); 303 304 if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety)) 305 { 306 G4ThreeVector eDirection(eDirX,eDirY,eDirZ); 307 eDirection.rotateUz(photonDirection0); 308 309 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),eDirection,eKineticEnergy) ; 310 aParticleChange.SetNumberOfSecondaries(1); 311 aParticleChange.AddSecondary(electron); 312 // Binding energy deposited locally 313 aParticleChange.ProposeLocalEnergyDeposit(bindingE); 314 } 315 else 316 { 317 aParticleChange.SetNumberOfSecondaries(0); 318 aParticleChange.ProposeLocalEnergyDeposit(eKineticEnergy + bindingE); 319 } 320 321 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep); 322 } 323 324 G4bool G4LowEnergyCompton::IsApplicable(const G4ParticleDefinition& particle) 325 { 326 return ( &particle == G4Gamma::Gamma() ); 327 } 328 329 G4double G4LowEnergyCompton::GetMeanFreePath(const G4Track& track, 330 G4double, // previousStepSize 331 G4ForceCondition*) 332 { 333 const G4DynamicParticle* photon = track.GetDynamicParticle(); 334 G4double energy = photon->GetKineticEnergy(); 335 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 336 size_t materialIndex = couple->GetIndex(); 337 338 G4double meanFreePath; 339 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex); 340 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX; 341 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex); 342 return meanFreePath; 343 } 344