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