| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
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 // 7 //
>> 8 // $Id: G4PhotoElectricEffect.cc,v 1.12 2001/02/22 16:05:43 maire Exp $
>> 9 // GEANT4 tag $Name: geant4-03-01 $
26 // 10 //
27 // << 11 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
28 //------------------ G4PhotoElectricEffect phy <<
29 // by Michel Maire, 24 May 1 <<
30 // <<
31 // Modified by Michel Maire and Vladimir Ivanc <<
32 // <<
33 // ------------------------------------------- <<
34 12
35 #include "G4PhotoElectricEffect.hh" << 13 // 12-06-96, Added SelectRandomAtom() method, by M.Maire
36 #include "G4SystemOfUnits.hh" << 14 // 21-06-96, SetCuts implementation, M.Maire
37 #include "G4PEEffectFluoModel.hh" << 15 // 17-09-96, PartialSumSigma(i)
38 #include "G4Electron.hh" << 16 // split of ComputeBindingEnergy, M.Maire
39 #include "G4EmParameters.hh" << 17 // 08-01-97, crossection table + meanfreepath table, M.Maire
>> 18 // 13-03-97, adapted for the new physics scheme, M.Maire
>> 19 // 28-03-97, protection in BuildPhysicsTable, M.Maire
>> 20 // 04-06-98, in DoIt, secondary production condition: range>G4std::min(threshold,safety)
>> 21 // 13-08-98, new methods SetBining() PrintInfo()
>> 22 // 17-11-98, use table of Atomic shells in PostStepDoIt
>> 23 // 06-01-99, use Sandia crossSection below 50 keV, V.Grichine mma
>> 24 // 20-05-99, protection against very low energy photons ,L.Urban
>> 25 // 08-06-99, removed this above protection from the DoIt. mma
>> 26 // 21-06-00, in DoIt, killing photon: aParticleChange.SetEnergyChange(0.); mma
>> 27 // 22-06-00, in DoIt, absorbe very low energy photon (back to 20-05-99); mma
>> 28 // 22-02-01, back to 08-06-99 after correc in SandiaTable (materials-V03-00-05)
>> 29 //
>> 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 32
41 //....oooOO0OOooo........oooOO0OOooo........oo << 33 #include "G4PhotoElectricEffect.hh"
>> 34 #include "G4EnergyLossTables.hh"
>> 35 #include "G4UnitsTable.hh"
42 36
43 using namespace std; << 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 38
>> 39 // constructor
>> 40
>> 41 G4PhotoElectricEffect::G4PhotoElectricEffect(const G4String& processName)
>> 42 : G4VDiscreteProcess (processName), // initialization
>> 43 theCrossSectionTable(NULL),
>> 44 theMeanFreePathTable(NULL),
>> 45 LowestEnergyLimit (50*keV),
>> 46 HighestEnergyLimit(50*MeV),
>> 47 NumbBinTable(100)
>> 48 { }
44 49
45 G4PhotoElectricEffect::G4PhotoElectricEffect(c << 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 G4ProcessType type):G4VEmProcess (processNam << 51
>> 52 // destructor
>> 53
>> 54 G4PhotoElectricEffect::~G4PhotoElectricEffect()
47 { 55 {
48 SetBuildTableFlag(false); << 56 if (theCrossSectionTable) {
49 SetSecondaryParticle(G4Electron::Electron()) << 57 theCrossSectionTable->clearAndDestroy();
50 SetProcessSubType(fPhotoElectricEffect); << 58 delete theCrossSectionTable;
51 SetMinKinEnergyPrim(200*keV); << 59 }
>> 60
>> 61 if (theMeanFreePathTable) {
>> 62 theMeanFreePathTable->clearAndDestroy();
>> 63 delete theMeanFreePathTable;
>> 64 }
52 } 65 }
53 66
54 //....oooOO0OOooo........oooOO0OOooo........oo << 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 68
56 G4PhotoElectricEffect::~G4PhotoElectricEffect( << 69 void G4PhotoElectricEffect::SetPhysicsTableBining(G4double lowE, G4double highE, G4int nBins)
>> 70 {
>> 71 LowestEnergyLimit = lowE; HighestEnergyLimit = highE; NumbBinTable = nBins;
>> 72 }
57 73
58 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 75
>> 76 void G4PhotoElectricEffect::BuildPhysicsTable(const G4ParticleDefinition& PhotonType)
59 77
60 G4bool G4PhotoElectricEffect::IsApplicable(con << 78 // Build microscopic cross section table and mean free path table
61 { 79 {
62 return (&p == G4Gamma::Gamma()); << 80 G4double LowEdgeEnergy, Value;
63 } << 81 G4PhysicsLogVector* ptrVector;
>> 82
>> 83 // Build microscopic cross section tables for the Photo Electric Effect
64 84
65 //....oooOO0OOooo........oooOO0OOooo........oo << 85 if (theCrossSectionTable) {
>> 86 theCrossSectionTable->clearAndDestroy(); delete theCrossSectionTable; }
66 87
67 void G4PhotoElectricEffect::InitialiseProcess( << 88 theCrossSectionTable = new G4PhysicsTable( G4Element::GetNumberOfElements()) ;
>> 89 const G4ElementTable* theElementTable = G4Element::GetElementTable() ;
>> 90 G4double AtomicNumber;
>> 91 G4int J;
>> 92
>> 93 for ( J=0 ; J < G4Element::GetNumberOfElements(); J++ )
>> 94 {
>> 95 //create physics vector then fill it ....
>> 96 ptrVector = new G4PhysicsLogVector(LowestEnergyLimit, HighestEnergyLimit,
>> 97 NumbBinTable ) ;
>> 98 AtomicNumber = (*theElementTable)(J)->GetZ();
>> 99
>> 100 for ( G4int i = 0 ; i < NumbBinTable ; i++ )
>> 101 {
>> 102 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
>> 103 Value = ComputeCrossSectionPerAtom( LowEdgeEnergy, AtomicNumber);
>> 104 ptrVector->PutValue( i , Value ) ;
>> 105 }
>> 106
>> 107 theCrossSectionTable->insertAt( J , ptrVector ) ;
>> 108
>> 109 }
>> 110
>> 111 // Build mean free path table for the Photo Electric Effect
>> 112
>> 113 if (theMeanFreePathTable) {
>> 114 theMeanFreePathTable->clearAndDestroy(); delete theMeanFreePathTable; }
>> 115
>> 116 theMeanFreePathTable = new G4PhysicsTable( G4Material::GetNumberOfMaterials() ) ;
>> 117 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
>> 118 G4Material* material;
>> 119
>> 120 for ( J=0 ; J < G4Material::GetNumberOfMaterials(); J++ )
>> 121 {
>> 122 //create physics vector then fill it ....
>> 123 ptrVector = new G4PhysicsLogVector(LowestEnergyLimit, HighestEnergyLimit,
>> 124 NumbBinTable ) ;
>> 125 material = (*theMaterialTable)(J);
>> 126
>> 127 for ( G4int i = 0 ; i < NumbBinTable ; i++ )
>> 128 {
>> 129 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
>> 130 Value = ComputeMeanFreePath( LowEdgeEnergy, material);
>> 131 ptrVector->PutValue( i , Value ) ;
>> 132 }
>> 133
>> 134 theMeanFreePathTable->insertAt( J , ptrVector ) ;
>> 135
>> 136 }
>> 137
>> 138 PrintInfoDefinition();
>> 139 }
>> 140
>> 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 142
>> 143 G4double G4PhotoElectricEffect::ComputeCrossSectionPerAtom (G4double PhotonEnergy,
>> 144 G4double AtomicNumber)
>> 145
>> 146 // Calculates the microscopic cross section in GEANT4 internal units.
>> 147 // A parametrized formula from L. Urban is used to estimate the total cross section.
>> 148 // It gives a good description of the elements : 5 < Atomic Number < 100 and
>> 149 // from 10 keV to 50 MeV.
>> 150
68 { 151 {
69 if(!isInitialised) { << 152 G4double CrossSection = 0.0 ;
70 isInitialised = true; << 153 if ( AtomicNumber < 1. ) return CrossSection;
71 if(nullptr == EmModel()) { SetEmModel(new << 154 if ( PhotonEnergy > 50.*MeV ) return CrossSection;
72 G4EmParameters* param = G4EmParameters::In << 155
73 EmModel()->SetLowEnergyLimit(param->MinKin << 156 static const G4double
74 EmModel()->SetHighEnergyLimit(param->MaxKi << 157 p1K =-8.8893e+2*nanobarn, p2K = 2.4394 *nanobarn, p3K = 2.8835e+2*nanobarn,
75 AddEmModel(1, EmModel()); << 158 p4K = 1.2133e+1*nanobarn, p5K =-3.1104e+2*nanobarn, p6K =-1.7284e-1*nanobarn,
76 } << 159 p7K = 1.4400e+1*nanobarn, p8K = 6.8357e+1*nanobarn, p9K = 7.3945e-4*nanobarn,
>> 160 p10K=-4.8149e-2*nanobarn, p11K= 5.5823e-1*nanobarn, p12K=-1.0089e-1*nanobarn;
>> 161 static const G4double
>> 162 p1L1=-1.0927e+3*nanobarn, p2L1=-9.7897e-1*nanobarn, p3L1= 1.2854e+2*nanobarn;
>> 163 static const G4double
>> 164 p1L2=-4.5803e+3*nanobarn, p2L2= 1.6858e-3*nanobarn, p3L2= 1.2013e+2*nanobarn;
>> 165 static const G4double
>> 166 p1M = 1.6924e+1*nanobarn;
>> 167
>> 168 const G4double pwZ = 3.845 , pwE = 2.975 ;
>> 169
>> 170 G4double Z = AtomicNumber, Z2 = Z*Z, Z3 = Z*Z*Z;
>> 171 G4double Em = PhotonEnergy/electron_mass_c2, Em2 = Em*Em, Em3 = Em*Em*Em;
>> 172
>> 173 CrossSection = pow(Z,pwZ)/pow(Em,pwE);
>> 174
>> 175 if (PhotonEnergy > ComputeKBindingEnergy(Z) ) {
>> 176 CrossSection *= (p1K/Z + p2K/Em + p3K + p4K*Z + p5K*Em
>> 177 + p6K*Z2 + p7K *Z *Em + p8K *Em2
>> 178 + p9K*Z3 + p10K*Z2*Em + p11K*Z*Em2 + p12K*Em3);
>> 179 if (CrossSection < 0.) CrossSection = 0. ;
>> 180 }
>> 181
>> 182 else if (PhotonEnergy > ComputeL1BindingEnergy(Z) ) {
>> 183 CrossSection *= (p1L1/Z + p2L1/Em + p3L1 );
>> 184 if (CrossSection < 0.) CrossSection = 0. ;
>> 185 }
>> 186
>> 187 else if (PhotonEnergy > ComputeL2BindingEnergy(Z) ) {
>> 188 CrossSection *= (p1L2/Z + p2L2/Em + p3L2 );
>> 189 if (CrossSection < 0.) CrossSection = 0. ;
>> 190 }
>> 191
>> 192 else CrossSection *= p1M;
>> 193
>> 194 return CrossSection;
>> 195 }
>> 196
>> 197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 198
>> 199 G4double G4PhotoElectricEffect::ComputeSandiaCrossSection(G4double PhotonEnergy,
>> 200 G4double AtomicNumber)
>> 201 {
>> 202 G4double energy2 = PhotonEnergy*PhotonEnergy, energy3 = PhotonEnergy*energy2,
>> 203 energy4 = energy2*energy2;
>> 204
>> 205 G4double* SandiaCof
>> 206 = G4SandiaTable::GetSandiaCofPerAtom((int)AtomicNumber,PhotonEnergy);
>> 207
>> 208 return SandiaCof[0]/PhotonEnergy + SandiaCof[1]/energy2 +
>> 209 SandiaCof[2]/energy3 + SandiaCof[3]/energy4;
>> 210 }
>> 211
>> 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 213
>> 214 G4VParticleChange* G4PhotoElectricEffect::PostStepDoIt(const G4Track& aTrack,
>> 215 const G4Step& aStep)
>> 216 //
>> 217 // Generate an electron resulting of a photo electric effect.
>> 218 // The incident photon disappear.
>> 219 // GEANT4 internal units
>> 220 //
>> 221
>> 222 { aParticleChange.Initialize(aTrack);
>> 223 G4Material* aMaterial = aTrack.GetMaterial();
>> 224
>> 225 const G4DynamicParticle* aDynamicPhoton = aTrack.GetDynamicParticle();
>> 226
>> 227 G4double PhotonEnergy = aDynamicPhoton->GetKineticEnergy();
>> 228 G4ParticleMomentum PhotonDirection = aDynamicPhoton->GetMomentumDirection();
>> 229
>> 230 // select randomly one element constituing the material.
>> 231 G4Element* anElement = SelectRandomAtom(aDynamicPhoton, aMaterial);
>> 232
>> 233 //
>> 234 // Photo electron
>> 235 //
>> 236
>> 237 G4int NbOfShells = anElement->GetNbOfAtomicShells();
>> 238 G4int i=0;
>> 239 while ((i<NbOfShells)&&(PhotonEnergy<anElement->GetAtomicShell(i))) i++;
>> 240
>> 241 if (i==NbOfShells) return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
>> 242
>> 243 G4double ElecKineEnergy = PhotonEnergy - anElement->GetAtomicShell(i);
>> 244 if ((G4EnergyLossTables::GetRange(G4Electron::Electron(),
>> 245 ElecKineEnergy,aMaterial)>aStep.GetPostStepPoint()->GetSafety())
>> 246 ||
>> 247 (ElecKineEnergy >
>> 248 (G4Electron::Electron()->GetCutsInEnergy())[aMaterial->GetIndex()]))
>> 249 {
>> 250 // the electron is created in the direction of the incident photon ...
>> 251 G4DynamicParticle* aElectron= new G4DynamicParticle (G4Electron::Electron(),
>> 252 PhotonDirection, ElecKineEnergy) ;
>> 253 aParticleChange.SetNumberOfSecondaries(1) ;
>> 254 aParticleChange.AddSecondary( aElectron ) ;
>> 255 }
>> 256 else
>> 257 {
>> 258 ElecKineEnergy = 0. ;
>> 259 aParticleChange.SetNumberOfSecondaries(0) ;
>> 260 }
>> 261
>> 262 //
>> 263 // Kill the incident photon
>> 264 //
>> 265 aParticleChange.SetLocalEnergyDeposit(PhotonEnergy-ElecKineEnergy);
>> 266 aParticleChange.SetEnergyChange(0.);
>> 267 aParticleChange.SetStatusChange(fStopAndKill);
>> 268
>> 269 // Reset NbOfInteractionLengthLeft and return aParticleChange
>> 270 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
77 } 271 }
78 272
79 //....oooOO0OOooo........oooOO0OOooo........oo << 273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 274
81 void G4PhotoElectricEffect::ProcessDescription << 275 G4Element*
>> 276 G4PhotoElectricEffect::SelectRandomAtom(const G4DynamicParticle* aDynamicPhoton,
>> 277 G4Material* aMaterial)
82 { 278 {
83 out << " Photoelectric effect"; << 279 // select randomly 1 element within the material
84 G4VEmProcess::ProcessDescription(out); << 280
>> 281 const G4int NumberOfElements = aMaterial->GetNumberOfElements();
>> 282 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 283 if (NumberOfElements == 1) return (*theElementVector)(0);
>> 284
>> 285 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
>> 286
>> 287 G4double PartialSumSigma = 0. ;
>> 288 G4double rval = G4UniformRand();
>> 289
>> 290 for ( G4int elm=0 ; elm < NumberOfElements ; elm++ )
>> 291 { PartialSumSigma += NbOfAtomsPerVolume[elm] *
>> 292 GetCrossSectionPerAtom(aDynamicPhoton,
>> 293 (*theElementVector)(elm));
>> 294 if (rval <= PartialSumSigma*MeanFreePath) return ((*theElementVector)(elm));
>> 295 }
>> 296 return ((*theElementVector)(NumberOfElements-1));
85 } 297 }
86 298
87 //....oooOO0OOooo........oooOO0OOooo........oo << 299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 300
>> 301 void G4PhotoElectricEffect::PrintInfoDefinition()
>> 302 {
>> 303 G4String comments = "Total cross sections from a parametrisation. ";
>> 304 comments += "Good description from 10 KeV to 50 MeV for all Z";
>> 305 comments += "\n Sandia crossSection below 50 KeV";
>> 306
>> 307 G4cout << G4endl << GetProcessName() << ": " << comments
>> 308 << "\n PhysicsTables from " << G4BestUnit(LowestEnergyLimit,"Energy")
>> 309 << " to " << G4BestUnit(HighestEnergyLimit,"Energy")
>> 310 << " in " << NumbBinTable << " bins. \n";
>> 311 }
>> 312
>> 313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 314