Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // Author: A. Forti 28 // Maria Grazia Pia (Maria.Grazia.Pia@ 29 // 30 // History: 31 // -------- 32 // Added Livermore data table construction met 33 // Modified BuildMeanFreePath to read new data 34 // Modified PostStepDoIt to insert sampling wi 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 39 // 22.01.2003 V.Ivanchenko - Cut per region 40 // 10.03.2003 V.Ivanchenko - Remove CutPerMate 41 // 24.04.2003 V.Ivanchenko - Cut per region mf 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 G 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 > intrinsicHighEnergyLim 80 { 81 G4Exception("G4LowEnergyCompton::G4LowEn 82 "OutOfRange", FatalException 83 "Energy outside intrinsic pr 84 } 85 86 crossSectionHandler = new G4RDCrossSectionHa 87 88 G4RDVDataSetAlgorithm* scatterInterpolation 89 G4String scatterFile = "comp/ce-sf-"; 90 scatterFunctionData = new G4RDCompositeEMDat 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 crea 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(con 119 { 120 121 crossSectionHandler->Clear(); 122 G4String crossSectionFile = "comp/ce-cs-"; 123 crossSectionHandler->LoadData(crossSectionFi 124 125 delete meanFreePathTable; 126 meanFreePathTable = crossSectionHandler->Bui 127 128 // For Doppler broadening 129 G4String file = "/doppler/shell-doppler"; 130 shellData.LoadData(file); 131 } 132 133 G4VParticleChange* G4LowEnergyCompton::PostSte 134 const G4Step& aStep) 135 { 136 // The scattered gamma energy is sampled acc 137 // then accepted or rejected depending on th 138 // by factor from Klein - Nishina formula. 139 // Expression of the angular distribution as 140 // angular and energy distribution and Scatt 141 // D. E. Cullen "A simple model of photon tr 142 // Phys. Res. B 101 (1995). Method of sampli 143 // data are interpolated while in the articl 144 // Reference to the article is from J. Stepa 145 // and Electron Interaction Data for GEANT i 146 // TeV (draft). 147 // The random number techniques of Butcher & 148 // (Nucl Phys 20(1960),15). 149 150 aParticleChange.Initialize(aTrack); 151 152 // Dynamic particle quantities 153 const G4DynamicParticle* incidentPhoton = aT 154 G4double photonEnergy0 = incidentPhoton->Get 155 156 if (photonEnergy0 <= lowEnergyLimit) 157 { 158 aParticleChange.ProposeTrackStatus(fStop 159 aParticleChange.ProposeEnergy(0.); 160 aParticleChange.ProposeLocalEnergyDeposi 161 return G4VDiscreteProcess::PostStepDoIt( 162 } 163 164 G4double e0m = photonEnergy0 / electron_mass 165 G4ParticleMomentum photonDirection0 = incide 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/photonE 171 172 // Select randomly one element in the curren 173 const G4MaterialCutsCouple* couple = aTrack. 174 G4int Z = crossSectionHandler->SelectRandomA 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) > G4UniformR 185 { 186 epsilon = std::exp(-alpha1 * G4UniformRand 187 epsilonSq = epsilon * epsilon; 188 } 189 else 190 { 191 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) 192 epsilon = std::sqrt(epsilonSq); 193 } 194 195 oneCosT = (1. - epsilon) / ( epsilon * e 196 sinT2 = oneCosT * (2. - oneCosT); 197 G4double x = std::sqrt(oneCosT/2.) / (wl 198 G4double scatteringFunction = scatterFun 199 gReject = (1. - epsilon * sinT2 / (1. + 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 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 * photonE 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.SelectRandomShel 227 bindingE = shellData.BindingEnergy(Z,she 228 229 eMax = photonEnergy0 - bindingE; 230 231 // Randomly sample bound electron moment 232 G4double pSample = profileData.RandomSel 233 // Rescale from atomic units 234 G4double pDoppler = pSample * fine_struc 235 G4double pDoppler2 = pDoppler * pDoppler 236 G4double var2 = 1. + oneCosT * e0m; 237 G4double var3 = var2*var2 - pDoppler2; 238 G4double var4 = var2 - pDoppler2 * cosTh 239 G4double var = var4*var4 - var3 + pDoppl 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 246 else photonE = (var4 + varSqrt) * scale; 247 } 248 else 249 { 250 photonE = -1.; 251 } 252 } while ( iteration <= maxDopplerIterations 253 (photonE < 0. || photonE > eMax || phot 254 255 // End of recalculation of photon energy wit 256 // Revert to original if maximum number of i 257 if (iteration >= maxDopplerIterations) 258 { 259 photonE = photonEoriginal; 260 bindingE = 0.; 261 } 262 263 // Update G4VParticleChange for the scattere 264 265 G4ThreeVector photonDirection1(dirX,dirY,dir 266 photonDirection1.rotateUz(photonDirection0); 267 aParticleChange.ProposeMomentumDirection(pho 268 269 G4double photonEnergy1 = photonE; 270 //G4cout << "--> PHOTONENERGY1 = " << photon 271 272 if (photonEnergy1 > 0.) 273 { 274 aParticleChange.ProposeEnergy(photonEner 275 } 276 else 277 { 278 aParticleChange.ProposeEnergy(0.) ; 279 aParticleChange.ProposeTrackStatus(fStop 280 } 281 282 // Kinematics of the scattered electron 283 G4double eKineticEnergy = photonEnergy0 - ph 284 G4double eTotalEnergy = eKineticEnergy + ele 285 286 G4double electronE = photonEnergy0 * (1. - e 287 G4double electronP2 = electronE*electronE - 288 G4double sinThetaE = -1.; 289 G4double cosThetaE = 0.; 290 if (electronP2 > 0.) 291 { 292 cosThetaE = (eTotalEnergy + photonEnergy 293 sinThetaE = -1. * std::sqrt(1. - cosThet 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 301 302 G4double safety = aStep.GetPostStepPoint()-> 303 304 if (rangeTest->Escape(G4Electron::Electron() 305 { 306 G4ThreeVector eDirection(eDirX,eDirY,eDi 307 eDirection.rotateUz(photonDirection0); 308 309 G4DynamicParticle* electron = new G4Dyna 310 aParticleChange.SetNumberOfSecondaries(1 311 aParticleChange.AddSecondary(electron); 312 // Binding energy deposited locally 313 aParticleChange.ProposeLocalEnergyDeposi 314 } 315 else 316 { 317 aParticleChange.SetNumberOfSecondaries(0 318 aParticleChange.ProposeLocalEnergyDeposi 319 } 320 321 return G4VDiscreteProcess::PostStepDoIt( aTr 322 } 323 324 G4bool G4LowEnergyCompton::IsApplicable(const 325 { 326 return ( &particle == G4Gamma::Gamma() ); 327 } 328 329 G4double G4LowEnergyCompton::GetMeanFreePath(c 330 G4double, // previousStepSize 331 G4ForceCondition*) 332 { 333 const G4DynamicParticle* photon = track.GetD 334 G4double energy = photon->GetKineticEnergy() 335 const G4MaterialCutsCouple* couple = track.G 336 size_t materialIndex = couple->GetIndex(); 337 338 G4double meanFreePath; 339 if (energy > highEnergyLimit) meanFreePath = 340 else if (energy < lowEnergyLimit) meanFreePa 341 else meanFreePath = meanFreePathTable->FindV 342 return meanFreePath; 343 } 344