Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // 23 // >> 24 // $Id: G4PhotoElectricEffect.cc,v 1.27 2002/05/02 11:37:22 maire Exp $ >> 25 // GEANT4 tag $Name: geant4-04-01 $ 27 // 26 // 28 //------------------ G4PhotoElectricEffect phy << 27 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 29 // by Michel Maire, 24 May 1 << 28 30 // << 29 // 12-06-96, Added SelectRandomAtom() method, by M.Maire 31 // Modified by Michel Maire and Vladimir Ivanc << 30 // 21-06-96, SetCuts implementation, M.Maire 32 // << 31 // 17-09-96, PartialSumSigma(i) 33 // ------------------------------------------- << 32 // split of ComputeBindingEnergy, M.Maire >> 33 // 08-01-97, crossection table + meanfreepath table, M.Maire >> 34 // 13-03-97, adapted for the new physics scheme, M.Maire >> 35 // 28-03-97, protection in BuildPhysicsTable, M.Maire >> 36 // 04-06-98, in DoIt, secondary production condition: >> 37 // range > G4std::min(threshold,safety) >> 38 // 13-08-98, new methods SetBining() PrintInfo() >> 39 // 17-11-98, use table of Atomic shells in PostStepDoIt >> 40 // 06-01-99, use Sandia crossSection below 50 keV, V.Grichine mma >> 41 // 20-05-99, protection against very low energy photons ,L.Urban >> 42 // 08-06-99, removed this above protection from the DoIt. mma >> 43 // 21-06-00, in DoIt, killing photon: aParticleChange.SetEnergyChange(0.); mma >> 44 // 22-06-00, in DoIt, absorbe very low energy photon (back to 20-05-99); mma >> 45 // 22-02-01, back to 08-06-99 after correc in SandiaTable (materials-V03-00-05) >> 46 // 28-05-01, V.Ivanchenko minor changes to provide ANSI -wall compilation >> 47 // 13-07-01, DoIt: suppression of production cut of the electron (mma) >> 48 // 06-08-01, new methods Store/Retrieve PhysicsTable (mma) >> 49 // 06-08-01, BuildThePhysicsTable() called from constructor (mma) >> 50 // 17-09-01, migration of Materials to pure STL (mma) >> 51 // 20-09-01, DoIt: fminimalEnergy of generated electron = 1*eV (mma) >> 52 // 01-10-01, come back to BuildPhysicsTable(const G4ParticleDefinition&) >> 53 // 10-01-02, moved few function from icc to cc >> 54 // 17-04-02, Keep only Sandia crossSections. Remove BuildPhysicsTables. >> 55 // Simplify public interface (mma) >> 56 // 29-04-02, Generate theta angle of the photoelectron from Sauter-Gavrila >> 57 // distribution (mma) >> 58 // >> 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 34 61 35 #include "G4PhotoElectricEffect.hh" 62 #include "G4PhotoElectricEffect.hh" 36 #include "G4SystemOfUnits.hh" << 63 #include "G4UnitsTable.hh" 37 #include "G4PEEffectFluoModel.hh" << 38 #include "G4Electron.hh" << 39 #include "G4EmParameters.hh" << 40 64 41 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 66 >> 67 // constructor >> 68 >> 69 G4PhotoElectricEffect::G4PhotoElectricEffect(const G4String& processName) >> 70 : G4VDiscreteProcess (processName), >> 71 fminimalEnergy(1*eV) >> 72 >> 73 { PrintInfoDefinition();} 42 74 43 using namespace std; << 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 76 >> 77 // destructor >> 78 >> 79 G4PhotoElectricEffect::~G4PhotoElectricEffect() >> 80 { } >> 81 >> 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 44 83 45 G4PhotoElectricEffect::G4PhotoElectricEffect(c << 84 G4double G4PhotoElectricEffect::ComputeCrossSectionPerAtom(G4double GammaEnergy, 46 G4ProcessType type):G4VEmProcess (processNam << 85 G4double AtomicNumber) >> 86 >> 87 // returns the photoElectric cross Section in GEANT4 internal units 47 { 88 { 48 SetBuildTableFlag(false); << 89 G4double* SandiaCof 49 SetSecondaryParticle(G4Electron::Electron()) << 90 = G4SandiaTable::GetSandiaCofPerAtom((int)AtomicNumber,GammaEnergy); 50 SetProcessSubType(fPhotoElectricEffect); << 91 51 SetMinKinEnergyPrim(200*keV); << 92 G4double energy2 = GammaEnergy*GammaEnergy, energy3 = GammaEnergy*energy2, >> 93 energy4 = energy2*energy2; >> 94 >> 95 return SandiaCof[0]/GammaEnergy + SandiaCof[1]/energy2 + >> 96 SandiaCof[2]/energy3 + SandiaCof[3]/energy4; >> 97 52 } 98 } 53 99 54 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 55 101 56 G4PhotoElectricEffect::~G4PhotoElectricEffect( << 102 G4double G4PhotoElectricEffect::ComputeMeanFreePath(G4double GammaEnergy, 57 << 103 G4Material* aMaterial) 58 //....oooOO0OOooo........oooOO0OOooo........oo << 59 104 60 G4bool G4PhotoElectricEffect::IsApplicable(con << 105 // returns the gamma mean free path in GEANT4 internal units 61 { 106 { 62 return (&p == G4Gamma::Gamma()); << 107 G4double* SandiaCof = aMaterial->GetSandiaTable() >> 108 ->GetSandiaCofForMaterial(GammaEnergy); >> 109 >> 110 G4double energy2 = GammaEnergy*GammaEnergy, energy3 = GammaEnergy*energy2, >> 111 energy4 = energy2*energy2; >> 112 >> 113 >> 114 G4double SIGMA = SandiaCof[0]/GammaEnergy + SandiaCof[1]/energy2 + >> 115 SandiaCof[2]/energy3 + SandiaCof[3]/energy4; >> 116 >> 117 return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX; 63 } 118 } 64 119 65 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 121 >> 122 G4VParticleChange* G4PhotoElectricEffect::PostStepDoIt(const G4Track& aTrack, >> 123 const G4Step& aStep) >> 124 // >> 125 // Generate an electron resulting of a photo electric effect. >> 126 // The incident photon disappear. >> 127 // GEANT4 internal units >> 128 // >> 129 >> 130 { aParticleChange.Initialize(aTrack); >> 131 G4Material* aMaterial = aTrack.GetMaterial(); >> 132 >> 133 const G4DynamicParticle* aDynamicPhoton = aTrack.GetDynamicParticle(); >> 134 >> 135 G4double PhotonEnergy = aDynamicPhoton->GetKineticEnergy(); >> 136 G4ParticleMomentum PhotonDirection = aDynamicPhoton->GetMomentumDirection(); >> 137 >> 138 // select randomly one element constituing the material. >> 139 G4Element* anElement = SelectRandomAtom(aDynamicPhoton, aMaterial); >> 140 >> 141 // >> 142 // Photo electron >> 143 // >> 144 >> 145 G4int NbOfShells = anElement->GetNbOfAtomicShells(); >> 146 G4int i=0; >> 147 while ((i<NbOfShells)&&(PhotonEnergy<anElement->GetAtomicShell(i))) i++; >> 148 >> 149 if (i==NbOfShells) return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); >> 150 >> 151 G4double ElecKineEnergy = PhotonEnergy - anElement->GetAtomicShell(i); >> 152 >> 153 if (ElecKineEnergy > fminimalEnergy) >> 154 { >> 155 // direction of the photo electron >> 156 // >> 157 G4double cosTeta = ElecThetaDistribution(ElecKineEnergy); >> 158 G4double sinTeta = sqrt(1.-cosTeta*cosTeta); >> 159 G4double Phi = twopi * G4UniformRand(); >> 160 G4double dirx = sinTeta*cos(Phi),diry = sinTeta*sin(Phi),dirz = cosTeta; >> 161 G4ThreeVector ElecDirection(dirx,diry,dirz); >> 162 ElecDirection.rotateUz(PhotonDirection); >> 163 // >> 164 G4DynamicParticle* aElectron = new G4DynamicParticle ( >> 165 G4Electron::Electron(),ElecDirection, ElecKineEnergy); >> 166 aParticleChange.SetNumberOfSecondaries(1); >> 167 aParticleChange.AddSecondary(aElectron); >> 168 } >> 169 else >> 170 { >> 171 ElecKineEnergy = 0.; >> 172 aParticleChange.SetNumberOfSecondaries(0); >> 173 } >> 174 >> 175 // >> 176 // Kill the incident photon >> 177 // >> 178 aParticleChange.SetLocalEnergyDeposit(PhotonEnergy-ElecKineEnergy); >> 179 aParticleChange.SetEnergyChange(0.); >> 180 aParticleChange.SetStatusChange(fStopAndKill); 66 181 67 void G4PhotoElectricEffect::InitialiseProcess( << 182 // Reset NbOfInteractionLengthLeft and return aParticleChange 68 { << 183 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 69 if(!isInitialised) { << 70 isInitialised = true; << 71 if(nullptr == EmModel()) { SetEmModel(new << 72 G4EmParameters* param = G4EmParameters::In << 73 EmModel()->SetLowEnergyLimit(param->MinKin << 74 EmModel()->SetHighEnergyLimit(param->MaxKi << 75 AddEmModel(1, EmModel()); << 76 } << 77 } 184 } 78 185 79 //....oooOO0OOooo........oooOO0OOooo........oo 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 80 187 81 void G4PhotoElectricEffect::ProcessDescription << 188 G4Element* G4PhotoElectricEffect::SelectRandomAtom( >> 189 const G4DynamicParticle* aDynamicPhoton, >> 190 G4Material* aMaterial) 82 { 191 { 83 out << " Photoelectric effect"; << 192 // select randomly 1 element within the material 84 G4VEmProcess::ProcessDescription(out); << 193 >> 194 const G4int NumberOfElements = aMaterial->GetNumberOfElements(); >> 195 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); >> 196 if (NumberOfElements == 1) return (*theElementVector)[0]; >> 197 >> 198 G4double GammaEnergy = aDynamicPhoton->GetKineticEnergy(); >> 199 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume(); >> 200 >> 201 G4double PartialSumSigma = 0. ; >> 202 G4double rval = G4UniformRand(); >> 203 >> 204 for ( G4int elm=0 ; elm < NumberOfElements ; elm++ ) >> 205 {PartialSumSigma += NbOfAtomsPerVolume[elm] * >> 206 ComputeCrossSectionPerAtom(GammaEnergy, >> 207 (*theElementVector)[elm]->GetZ()); >> 208 if (rval<=PartialSumSigma*MeanFreePath) return ((*theElementVector)[elm]); >> 209 } >> 210 return ((*theElementVector)[NumberOfElements-1]); 85 } 211 } >> 212 >> 213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 214 >> 215 G4double G4PhotoElectricEffect::ElecThetaDistribution(G4double kineEnergy) >> 216 { >> 217 // Compute Theta distribution of the emitted electron, with respect to the >> 218 // incident Gamma. >> 219 // The Sauter-Gavrila distribution for the K-shell is used. >> 220 // >> 221 G4double gamma = 1. + kineEnergy/electron_mass_c2; >> 222 G4double beta = sqrt(gamma*gamma-1.)/gamma; >> 223 G4double b = 0.5*gamma*(gamma-1.)*(gamma-2); >> 224 >> 225 G4double rndm,costeta,term,greject,grejsup; >> 226 if (gamma < 2.) grejsup = gamma*gamma*(1.+b-beta*b); >> 227 else grejsup = gamma*gamma*(1.+b+beta*b); >> 228 >> 229 do { rndm = 1.-2*G4UniformRand(); >> 230 costeta = (rndm+beta)/(rndm*beta+1.); >> 231 term = 1.-beta*costeta; >> 232 greject = (1.-costeta*costeta)*(1.+b*term)/(term*term); >> 233 } while(greject < G4UniformRand()*grejsup); >> 234 >> 235 return costeta; >> 236 >> 237 } >> 238 >> 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 240 >> 241 void G4PhotoElectricEffect::PrintInfoDefinition() >> 242 { >> 243 G4String comments = "Total cross sections from Sandia parametrisation. "; >> 244 >> 245 G4cout << G4endl << GetProcessName() << ": " << comments << G4endl; >> 246 } 86 247 87 //....oooOO0OOooo........oooOO0OOooo........oo 248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 88 249