Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // Author: Luciano Pandola 27 // Author: Luciano Pandola 28 // 28 // 29 // History: 29 // History: 30 // -------- 30 // -------- 31 // 15 Feb 2010 L Pandola Implementation 31 // 15 Feb 2010 L Pandola Implementation 32 // 18 Mar 2010 L Pandola Removed GetAtomsPe 32 // 18 Mar 2010 L Pandola Removed GetAtomsPerMolecule(), now demanded 33 // to G4PenelopeOsc 33 // to G4PenelopeOscillatorManager 34 // 01 Feb 2011 L Pandola Suppress fake ener 34 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger is 35 // active. 35 // active. 36 // Make sure that flu 36 // Make sure that fluorescence/Auger is generated only 37 // if above thresho 37 // if above threshold 38 // 24 May 2011 L Pandola Renamed (make v200 38 // 24 May 2011 L Pandola Renamed (make v2008 as default Penelope) 39 // 10 Jun 2011 L Pandola Migrate atomic dee 39 // 10 Jun 2011 L Pandola Migrate atomic deexcitation interface 40 // 09 Oct 2013 L Pandola Migration to MT 40 // 09 Oct 2013 L Pandola Migration to MT 41 // 25 Jul 2023 D Iuso Fix for possible i << 42 // 41 // 43 #include "G4PenelopeComptonModel.hh" 42 #include "G4PenelopeComptonModel.hh" 44 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh" 45 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh" 46 #include "G4ParticleDefinition.hh" 45 #include "G4ParticleDefinition.hh" 47 #include "G4MaterialCutsCouple.hh" 46 #include "G4MaterialCutsCouple.hh" 48 #include "G4DynamicParticle.hh" 47 #include "G4DynamicParticle.hh" 49 #include "G4VEMDataSet.hh" 48 #include "G4VEMDataSet.hh" 50 #include "G4PhysicsTable.hh" 49 #include "G4PhysicsTable.hh" 51 #include "G4PhysicsLogVector.hh" 50 #include "G4PhysicsLogVector.hh" 52 #include "G4AtomicTransitionManager.hh" 51 #include "G4AtomicTransitionManager.hh" 53 #include "G4AtomicShell.hh" 52 #include "G4AtomicShell.hh" 54 #include "G4Gamma.hh" 53 #include "G4Gamma.hh" 55 #include "G4Electron.hh" 54 #include "G4Electron.hh" 56 #include "G4PenelopeOscillatorManager.hh" 55 #include "G4PenelopeOscillatorManager.hh" 57 #include "G4PenelopeOscillator.hh" 56 #include "G4PenelopeOscillator.hh" 58 #include "G4LossTableManager.hh" 57 #include "G4LossTableManager.hh" 59 #include "G4Exp.hh" 58 #include "G4Exp.hh" 60 59 61 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 62 61 >> 62 63 G4PenelopeComptonModel::G4PenelopeComptonModel 63 G4PenelopeComptonModel::G4PenelopeComptonModel(const G4ParticleDefinition* part, 64 const G4String& nam) 64 const G4String& nam) 65 :G4VEmModel(nam),fParticleChange(nullptr),fP << 65 :G4VEmModel(nam),fParticleChange(0),fParticle(0), 66 fAtomDeexcitation(nullptr), << 66 isInitialised(false),fAtomDeexcitation(0), 67 fOscManager(nullptr),fIsInitialised(false) << 67 oscManager(0) 68 { 68 { 69 fIntrinsicLowEnergyLimit = 100.0*eV; 69 fIntrinsicLowEnergyLimit = 100.0*eV; 70 fIntrinsicHighEnergyLimit = 100.0*GeV; 70 fIntrinsicHighEnergyLimit = 100.0*GeV; 71 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 71 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 72 // 72 // 73 fOscManager = G4PenelopeOscillatorManager::G << 73 oscManager = G4PenelopeOscillatorManager::GetOscillatorManager(); 74 74 75 if (part) 75 if (part) 76 SetParticle(part); 76 SetParticle(part); 77 77 78 fVerboseLevel= 0; << 78 verboseLevel= 0; 79 // Verbosity scale: 79 // Verbosity scale: 80 // 0 = nothing 80 // 0 = nothing 81 // 1 = warning for energy non-conservation 81 // 1 = warning for energy non-conservation 82 // 2 = details of energy budget 82 // 2 = details of energy budget 83 // 3 = calculation of cross sections, file o 83 // 3 = calculation of cross sections, file openings, sampling of atoms 84 // 4 = entering in methods 84 // 4 = entering in methods 85 85 86 //Mark this model as "applicable" for atomic 86 //Mark this model as "applicable" for atomic deexcitation 87 SetDeexcitationFlag(true); 87 SetDeexcitationFlag(true); 88 88 89 fTransitionManager = G4AtomicTransitionManag 89 fTransitionManager = G4AtomicTransitionManager::Instance(); 90 } 90 } 91 91 92 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 93 93 94 G4PenelopeComptonModel::~G4PenelopeComptonMode 94 G4PenelopeComptonModel::~G4PenelopeComptonModel() 95 {;} 95 {;} 96 96 97 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 98 99 void G4PenelopeComptonModel::Initialise(const 99 void G4PenelopeComptonModel::Initialise(const G4ParticleDefinition* part, 100 const G4DataVector&) 100 const G4DataVector&) 101 { 101 { 102 if (fVerboseLevel > 3) << 102 if (verboseLevel > 3) 103 G4cout << "Calling G4PenelopeComptonModel: 103 G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl; 104 104 105 fAtomDeexcitation = G4LossTableManager::Inst 105 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 106 //Issue warning if the AtomicDeexcitation ha 106 //Issue warning if the AtomicDeexcitation has not been declared 107 if (!fAtomDeexcitation) 107 if (!fAtomDeexcitation) 108 { 108 { 109 G4cout << G4endl; 109 G4cout << G4endl; 110 G4cout << "WARNING from G4PenelopeCompto 110 G4cout << "WARNING from G4PenelopeComptonModel " << G4endl; 111 G4cout << "Atomic de-excitation module i 111 G4cout << "Atomic de-excitation module is not instantiated, so there will not be "; 112 G4cout << "any fluorescence/Auger emissi 112 G4cout << "any fluorescence/Auger emission." << G4endl; 113 G4cout << "Please make sure this is inte 113 G4cout << "Please make sure this is intended" << G4endl; 114 } 114 } 115 115 116 SetParticle(part); 116 SetParticle(part); 117 117 118 if (IsMaster() && part == fParticle) 118 if (IsMaster() && part == fParticle) 119 { 119 { 120 120 121 if (fVerboseLevel > 0) << 121 if (verboseLevel > 0) 122 { 122 { 123 G4cout << "Penelope Compton model v2008 is 123 G4cout << "Penelope Compton model v2008 is initialized " << G4endl 124 << "Energy range: " 124 << "Energy range: " 125 << LowEnergyLimit() / keV << " keV - " 125 << LowEnergyLimit() / keV << " keV - " 126 << HighEnergyLimit() / GeV << " GeV"; 126 << HighEnergyLimit() / GeV << " GeV"; 127 } 127 } 128 //Issue a warning, if the model is going 128 //Issue a warning, if the model is going to be used down to a 129 //energy which is outside the validity o 129 //energy which is outside the validity of the model itself 130 if (LowEnergyLimit() < fIntrinsicLowEner 130 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit) 131 { 131 { 132 G4ExceptionDescription ed; 132 G4ExceptionDescription ed; 133 ed << "Using the Penelope Compton model ou 133 ed << "Using the Penelope Compton model outside its intrinsic validity range. " 134 << G4endl; 134 << G4endl; 135 ed << "-> LowEnergyLimit() in process = " 135 ed << "-> LowEnergyLimit() in process = " << LowEnergyLimit()/keV << "keV " << G4endl; 136 ed << "-> Instrinsic low-energy limit = " 136 ed << "-> Instrinsic low-energy limit = " << fIntrinsicLowEnergyLimit/keV << "keV " 137 << G4endl; 137 << G4endl; 138 ed << "Result of the simulation have to be 138 ed << "Result of the simulation have to be taken with care" << G4endl; 139 G4Exception("G4PenelopeComptonModel::Initi 139 G4Exception("G4PenelopeComptonModel::Initialise()", 140 "em2100",JustWarning,ed); 140 "em2100",JustWarning,ed); 141 } 141 } 142 } 142 } 143 143 144 if(fIsInitialised) return; << 144 if(isInitialised) return; 145 fParticleChange = GetParticleChangeForGamma( 145 fParticleChange = GetParticleChangeForGamma(); 146 fIsInitialised = true; << 146 isInitialised = true; 147 147 148 } 148 } 149 149 150 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 151 151 152 void G4PenelopeComptonModel::InitialiseLocal(c 152 void G4PenelopeComptonModel::InitialiseLocal(const G4ParticleDefinition* part, 153 153 G4VEmModel *masterModel) 154 { 154 { 155 if (fVerboseLevel > 3) << 155 if (verboseLevel > 3) 156 G4cout << "Calling G4PenelopeComptonModel 156 G4cout << "Calling G4PenelopeComptonModel::InitialiseLocal()" << G4endl; >> 157 157 // 158 // 158 //Check that particle matches: one might hav 159 //Check that particle matches: one might have multiple master models (e.g. 159 //for e+ and e-). 160 //for e+ and e-). 160 // 161 // 161 if (part == fParticle) 162 if (part == fParticle) 162 { 163 { 163 //Get the const table pointers from the 164 //Get the const table pointers from the master to the workers 164 const G4PenelopeComptonModel* theModel = 165 const G4PenelopeComptonModel* theModel = 165 static_cast<G4PenelopeComptonModel*> ( 166 static_cast<G4PenelopeComptonModel*> (masterModel); 166 167 167 //Same verbosity for all workers, as the 168 //Same verbosity for all workers, as the master 168 fVerboseLevel = theModel->fVerboseLevel; << 169 verboseLevel = theModel->verboseLevel; 169 } 170 } >> 171 170 return; 172 return; 171 } 173 } 172 174 173 175 174 //....oooOO0OOooo........oooOO0OOooo........oo 176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 175 177 176 G4double G4PenelopeComptonModel::CrossSectionP 178 G4double G4PenelopeComptonModel::CrossSectionPerVolume(const G4Material* material, 177 179 const G4ParticleDefinition* p, 178 180 G4double energy, 179 181 G4double, 180 182 G4double) 181 { 183 { 182 // Penelope model v2008 to calculate the Com 184 // Penelope model v2008 to calculate the Compton scattering cross section: 183 // D. Brusa et al., Nucl. Instrum. Meth. A 3 185 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167 184 // 186 // 185 // The cross section for Compton scattering 187 // The cross section for Compton scattering is calculated according to the Klein-Nishina 186 // formula for energy > 5 MeV. 188 // formula for energy > 5 MeV. 187 // For E < 5 MeV it is used a parametrizatio 189 // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega, 188 // which is integrated numerically in cos(th 190 // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection(). 189 // The parametrization includes the J(p) 191 // The parametrization includes the J(p) 190 // distribution profiles for the atomic shel 192 // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations 191 // from F. Biggs et al., At. Data Nucl. Data 193 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201 192 // 194 // 193 if (fVerboseLevel > 3) << 195 if (verboseLevel > 3) 194 G4cout << "Calling CrossSectionPerVolume() 196 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeComptonModel" << G4endl; 195 SetupForMaterial(p, material, energy); 197 SetupForMaterial(p, material, energy); 196 198 >> 199 197 G4double cs = 0; 200 G4double cs = 0; 198 //Force null cross-section if below the low- 201 //Force null cross-section if below the low-energy edge of the table 199 if (energy < LowEnergyLimit()) 202 if (energy < LowEnergyLimit()) 200 return cs; 203 return cs; 201 204 202 //Retrieve the oscillator table for this mat 205 //Retrieve the oscillator table for this material 203 G4PenelopeOscillatorTable* theTable = fOscMa << 206 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); 204 207 205 if (energy < 5*MeV) //explicit calculation f 208 if (energy < 5*MeV) //explicit calculation for E < 5 MeV 206 { 209 { 207 size_t numberOfOscillators = theTable->s 210 size_t numberOfOscillators = theTable->size(); 208 for (size_t i=0;i<numberOfOscillators;i+ 211 for (size_t i=0;i<numberOfOscillators;i++) 209 { 212 { 210 G4PenelopeOscillator* theOsc = (*theTable) 213 G4PenelopeOscillator* theOsc = (*theTable)[i]; 211 //sum contributions coming from each oscil 214 //sum contributions coming from each oscillator 212 cs += OscillatorTotalCrossSection(energy,t 215 cs += OscillatorTotalCrossSection(energy,theOsc); 213 } 216 } 214 } 217 } 215 else //use Klein-Nishina for E>5 MeV 218 else //use Klein-Nishina for E>5 MeV 216 cs = KleinNishinaCrossSection(energy,mater 219 cs = KleinNishinaCrossSection(energy,material); 217 220 218 //cross sections are in units of pi*classic_ 221 //cross sections are in units of pi*classic_electr_radius^2 219 cs *= pi*classic_electr_radius*classic_elect 222 cs *= pi*classic_electr_radius*classic_electr_radius; 220 223 221 //Now, cs is the cross section *per molecule 224 //Now, cs is the cross section *per molecule*, let's calculate the 222 //cross section per volume 225 //cross section per volume >> 226 223 G4double atomDensity = material->GetTotNbOfA 227 G4double atomDensity = material->GetTotNbOfAtomsPerVolume(); 224 G4double atPerMol = fOscManager->GetAtomsPe << 228 G4double atPerMol = oscManager->GetAtomsPerMolecule(material); 225 229 226 if (fVerboseLevel > 3) << 230 if (verboseLevel > 3) 227 G4cout << "Material " << material->GetName 231 G4cout << "Material " << material->GetName() << " has " << atPerMol << 228 "atoms per molecule" << G4endl; 232 "atoms per molecule" << G4endl; 229 233 230 G4double moleculeDensity = 0.; 234 G4double moleculeDensity = 0.; 231 235 232 if (atPerMol) 236 if (atPerMol) 233 moleculeDensity = atomDensity/atPerMol; 237 moleculeDensity = atomDensity/atPerMol; 234 238 235 G4double csvolume = cs*moleculeDensity; 239 G4double csvolume = cs*moleculeDensity; 236 240 237 if (fVerboseLevel > 2) << 241 if (verboseLevel > 2) 238 G4cout << "Compton mean free path at " << 242 G4cout << "Compton mean free path at " << energy/keV << " keV for material " << 239 material->GetName() << " = " << (1 243 material->GetName() << " = " << (1./csvolume)/mm << " mm" << G4endl; 240 return csvolume; 244 return csvolume; 241 } 245 } 242 246 243 //....oooOO0OOooo........oooOO0OOooo........oo 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 244 248 245 //This is a dummy method. Never inkoved by the 249 //This is a dummy method. Never inkoved by the tracking, it just issues 246 //a warning if one tries to get Cross Sections 250 //a warning if one tries to get Cross Sections per Atom via the 247 //G4EmCalculator. 251 //G4EmCalculator. 248 G4double G4PenelopeComptonModel::ComputeCrossS 252 G4double G4PenelopeComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 249 253 G4double, 250 254 G4double, 251 255 G4double, 252 256 G4double, 253 257 G4double) 254 { 258 { 255 G4cout << "*** G4PenelopeComptonModel -- WAR 259 G4cout << "*** G4PenelopeComptonModel -- WARNING ***" << G4endl; 256 G4cout << "Penelope Compton model v2008 does 260 G4cout << "Penelope Compton model v2008 does not calculate cross section _per atom_ " << G4endl; 257 G4cout << "so the result is always zero. For 261 G4cout << "so the result is always zero. For physics values, please invoke " << G4endl; 258 G4cout << "GetCrossSectionPerVolume() or Get 262 G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl; 259 return 0; 263 return 0; 260 } 264 } 261 265 262 //....oooOO0OOooo........oooOO0OOooo........oo 266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 263 267 264 void G4PenelopeComptonModel::SampleSecondaries 268 void G4PenelopeComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 265 const G4MaterialCutsCouple* c 269 const G4MaterialCutsCouple* couple, 266 const G4DynamicParticle* aDyna 270 const G4DynamicParticle* aDynamicGamma, 267 G4double, 271 G4double, 268 G4double) 272 G4double) 269 { 273 { >> 274 270 // Penelope model v2008 to sample the Compto 275 // Penelope model v2008 to sample the Compton scattering final state. 271 // D. Brusa et al., Nucl. Instrum. Meth. A 3 276 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167 272 // The model determines also the original sh 277 // The model determines also the original shell from which the electron is expelled, 273 // in order to produce fluorescence de-excit 278 // in order to produce fluorescence de-excitation (from G4DeexcitationManager) 274 // 279 // 275 // The final state for Compton scattering is 280 // The final state for Compton scattering is calculated according to the Klein-Nishina 276 // formula for energy > 5 MeV. In this case, 281 // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and 277 // one can assume that the target electron i 282 // one can assume that the target electron is at rest. 278 // For E < 5 MeV it is used the parametrizat 283 // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega, 279 // to sample the scattering angle and the en 284 // to sample the scattering angle and the energy of the emerging electron, which is 280 // G4PenelopeComptonModel::DifferentialCross 285 // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is 281 // used to sample cos(theta). The efficiency 286 // used to sample cos(theta). The efficiency increases monotonically with photon energy and is 282 // nearly independent on the Z; typical valu 287 // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV, 283 // respectively. 288 // respectively. 284 // The parametrization includes the J(p) dis 289 // The parametrization includes the J(p) distribution profiles for the atomic shells, that are 285 // tabulated 290 // tabulated 286 // from Hartree-Fock calculations from F. Bi 291 // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201. 287 // Doppler broadening is included. 292 // Doppler broadening is included. 288 // 293 // 289 294 290 if (fVerboseLevel > 3) << 295 if (verboseLevel > 3) 291 G4cout << "Calling SampleSecondaries() of 296 G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl; 292 297 293 G4double photonEnergy0 = aDynamicGamma->GetK 298 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); 294 299 295 // do nothing below the threshold 300 // do nothing below the threshold 296 // should never get here because the XS is z 301 // should never get here because the XS is zero below the limit 297 if(photonEnergy0 < LowEnergyLimit()) 302 if(photonEnergy0 < LowEnergyLimit()) 298 return; 303 return; 299 304 300 G4ParticleMomentum photonDirection0 = aDynam 305 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); 301 const G4Material* material = couple->GetMate 306 const G4Material* material = couple->GetMaterial(); 302 307 303 G4PenelopeOscillatorTable* theTable = fOscMa << 308 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); 304 309 305 const G4int nmax = 64; 310 const G4int nmax = 64; 306 G4double rn[nmax]={0.0}; 311 G4double rn[nmax]={0.0}; 307 G4double pac[nmax]={0.0}; 312 G4double pac[nmax]={0.0}; 308 313 309 G4double S=0.0; 314 G4double S=0.0; 310 G4double epsilon = 0.0; 315 G4double epsilon = 0.0; 311 G4double cosTheta = 1.0; 316 G4double cosTheta = 1.0; 312 G4double hartreeFunc = 0.0; 317 G4double hartreeFunc = 0.0; 313 G4double oscStren = 0.0; 318 G4double oscStren = 0.0; 314 size_t numberOfOscillators = theTable->size( 319 size_t numberOfOscillators = theTable->size(); 315 size_t targetOscillator = 0; 320 size_t targetOscillator = 0; 316 G4double ionEnergy = 0.0*eV; 321 G4double ionEnergy = 0.0*eV; 317 322 318 G4double ek = photonEnergy0/electron_mass_c2 323 G4double ek = photonEnergy0/electron_mass_c2; 319 G4double ek2 = 2.*ek+1.0; 324 G4double ek2 = 2.*ek+1.0; 320 G4double eks = ek*ek; 325 G4double eks = ek*ek; 321 G4double ek1 = eks-ek2-1.0; 326 G4double ek1 = eks-ek2-1.0; 322 327 323 G4double taumin = 1.0/ek2; 328 G4double taumin = 1.0/ek2; 324 //This is meant to fix a possible (rare) flo << 329 G4double a1 = std::log(ek2); 325 //causing an infinite loop. The maximum of t << 326 //be represented (i.e. ~ 1. - 1e-16). Fix by << 327 static G4double taumax = std::nexttoward(1.0 << 328 if (fVerboseLevel > 3) << 329 G4cout << "G4PenelopeComptonModel: maximum << 330 //To here. << 331 G4double a1 = G4Log(ek2); << 332 G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2); 330 G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2); 333 331 334 G4double TST = 0; 332 G4double TST = 0; 335 G4double tau = 0.; 333 G4double tau = 0.; 336 334 337 //If the incoming photon is above 5 MeV, the 335 //If the incoming photon is above 5 MeV, the quicker approach based on the 338 //pure Klein-Nishina formula is used 336 //pure Klein-Nishina formula is used 339 if (photonEnergy0 > 5*MeV) 337 if (photonEnergy0 > 5*MeV) 340 { 338 { 341 do{ 339 do{ 342 do{ 340 do{ 343 if ((a2*G4UniformRand()) < a1) 341 if ((a2*G4UniformRand()) < a1) 344 tau = std::pow(taumin,G4UniformRand()); 342 tau = std::pow(taumin,G4UniformRand()); 345 else 343 else 346 tau = std::sqrt(1.0+G4UniformRand()*(tau 344 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0)); 347 //rejection function 345 //rejection function 348 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(e 346 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau)); 349 }while (G4UniformRand()> TST); 347 }while (G4UniformRand()> TST); 350 if (tau > taumax) tau = taumax; //prevent FP << 351 epsilon=tau; 348 epsilon=tau; 352 cosTheta = 1.0 - (1.0-tau)/(ek*tau); 349 cosTheta = 1.0 - (1.0-tau)/(ek*tau); 353 350 354 //Target shell electrons 351 //Target shell electrons 355 TST = fOscManager->GetTotalZ(material)*G4Uni << 352 TST = oscManager->GetTotalZ(material)*G4UniformRand(); 356 targetOscillator = numberOfOscillators-1; // 353 targetOscillator = numberOfOscillators-1; //last level 357 S=0.0; 354 S=0.0; 358 G4bool levelFound = false; 355 G4bool levelFound = false; 359 for (size_t j=0;j<numberOfOscillators && !le 356 for (size_t j=0;j<numberOfOscillators && !levelFound; j++) 360 { 357 { 361 S += (*theTable)[j]->GetOscillatorStreng 358 S += (*theTable)[j]->GetOscillatorStrength(); 362 if (S > TST) 359 if (S > TST) 363 { 360 { 364 targetOscillator = j; 361 targetOscillator = j; 365 levelFound = true; 362 levelFound = true; 366 } 363 } 367 } 364 } 368 //check whether the level is valid 365 //check whether the level is valid 369 ionEnergy = (*theTable)[targetOscillator]->G 366 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy(); 370 }while((epsilon*photonEnergy0-photonEner 367 }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0); 371 } 368 } 372 else //photonEnergy0 < 5 MeV 369 else //photonEnergy0 < 5 MeV 373 { 370 { 374 //Incoherent scattering function for the 371 //Incoherent scattering function for theta=PI 375 G4double s0=0.0; 372 G4double s0=0.0; 376 G4double pzomc=0.0; 373 G4double pzomc=0.0; 377 G4double rni=0.0; 374 G4double rni=0.0; 378 G4double aux=0.0; 375 G4double aux=0.0; 379 for (size_t i=0;i<numberOfOscillators;i+ 376 for (size_t i=0;i<numberOfOscillators;i++) 380 { 377 { 381 ionEnergy = (*theTable)[i]->GetIonisationE 378 ionEnergy = (*theTable)[i]->GetIonisationEnergy(); 382 if (photonEnergy0 > ionEnergy) 379 if (photonEnergy0 > ionEnergy) 383 { 380 { 384 G4double aux2 = photonEnergy0*(photonE 381 G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0; 385 hartreeFunc = (*theTable)[i]->GetHartr 382 hartreeFunc = (*theTable)[i]->GetHartreeFactor(); 386 oscStren = (*theTable)[i]->GetOscillat 383 oscStren = (*theTable)[i]->GetOscillatorStrength(); 387 pzomc = hartreeFunc*(aux2-electron_mas 384 pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/ 388 (electron_mass_c2*std::sqrt(2.0*aux2+ionEn 385 (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy)); 389 if (pzomc > 0) 386 if (pzomc > 0) 390 rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+st 387 rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)* 391 (std::sqrt(0.5)+std::sqrt(2.0)* 388 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc)); 392 else 389 else 393 rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::s 390 rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)* 394 (std::sqrt(0.5)-std::sqrt(2.0)*pzom 391 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc)); 395 s0 += oscStren*rni; 392 s0 += oscStren*rni; 396 } 393 } 397 } 394 } 398 //Sampling tau 395 //Sampling tau 399 G4double cdt1 = 0.; 396 G4double cdt1 = 0.; 400 do 397 do 401 { 398 { 402 if ((G4UniformRand()*a2) < a1) 399 if ((G4UniformRand()*a2) < a1) 403 tau = std::pow(taumin,G4UniformRand()); 400 tau = std::pow(taumin,G4UniformRand()); 404 else 401 else 405 tau = std::sqrt(1.0+G4UniformRand()*(tau 402 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0)); 406 if (tau > taumax) tau = taumax; //prevent << 407 cdt1 = (1.0-tau)/(ek*tau); 403 cdt1 = (1.0-tau)/(ek*tau); 408 //Incoherent scattering function 404 //Incoherent scattering function 409 S = 0.; 405 S = 0.; 410 for (size_t i=0;i<numberOfOscillators;i++) 406 for (size_t i=0;i<numberOfOscillators;i++) 411 { 407 { 412 ionEnergy = (*theTable)[i]->GetIonisat 408 ionEnergy = (*theTable)[i]->GetIonisationEnergy(); 413 if (photonEnergy0 > ionEnergy) //sum o 409 if (photonEnergy0 > ionEnergy) //sum only on excitable levels 414 { 410 { 415 aux = photonEnergy0*(photonEnergy0-ionEn 411 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1; 416 hartreeFunc = (*theTable)[i]->GetHartree 412 hartreeFunc = (*theTable)[i]->GetHartreeFactor(); 417 oscStren = (*theTable)[i]->GetOscillator 413 oscStren = (*theTable)[i]->GetOscillatorStrength(); 418 pzomc = hartreeFunc*(aux-electron_mass_c 414 pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/ 419 (electron_mass_c2*std::sqrt(2.0*aux+io 415 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy)); 420 if (pzomc > 0) 416 if (pzomc > 0) 421 rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0 417 rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)* 422 (std::sqrt(0.5)+std::sqrt(2.0)* 418 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc)); 423 else 419 else 424 rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)- 420 rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)* 425 (std::sqrt(0.5)-std::sqrt(2.0)*pzom 421 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc)); 426 S += oscStren*rn[i]; 422 S += oscStren*rn[i]; 427 pac[i] = S; 423 pac[i] = S; 428 } 424 } 429 else 425 else 430 pac[i] = S-1e-6; 426 pac[i] = S-1e-6; 431 } 427 } 432 //Rejection function 428 //Rejection function 433 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/ 429 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau)); 434 }while ((G4UniformRand()*s0) > TST); 430 }while ((G4UniformRand()*s0) > TST); 435 431 436 cosTheta = 1.0 - cdt1; 432 cosTheta = 1.0 - cdt1; 437 G4double fpzmax=0.0,fpz=0.0; 433 G4double fpzmax=0.0,fpz=0.0; 438 G4double A=0.0; 434 G4double A=0.0; 439 //Target electron shell 435 //Target electron shell 440 do 436 do 441 { 437 { 442 do 438 do 443 { 439 { 444 TST = S*G4UniformRand(); 440 TST = S*G4UniformRand(); 445 targetOscillator = numberOfOscillators 441 targetOscillator = numberOfOscillators-1; //last level 446 G4bool levelFound = false; 442 G4bool levelFound = false; 447 for (size_t i=0;i<numberOfOscillators 443 for (size_t i=0;i<numberOfOscillators && !levelFound;i++) 448 { 444 { 449 if (pac[i]>TST) 445 if (pac[i]>TST) 450 { 446 { 451 targetOscillator = i; 447 targetOscillator = i; 452 levelFound = true; 448 levelFound = true; 453 } 449 } 454 } 450 } 455 A = G4UniformRand()*rn[targetOscillato 451 A = G4UniformRand()*rn[targetOscillator]; 456 hartreeFunc = (*theTable)[targetOscill 452 hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor(); 457 oscStren = (*theTable)[targetOscillato 453 oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength(); 458 if (A < 0.5) 454 if (A < 0.5) 459 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-G4Lo << 455 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/ 460 (std::sqrt(2.0)*hartreeFunc); 456 (std::sqrt(2.0)*hartreeFunc); 461 else 457 else 462 pzomc = (std::sqrt(0.5-G4Log(2.0-2.0*A))-s << 458 pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/ 463 (std::sqrt(2.0)*hartreeFunc); 459 (std::sqrt(2.0)*hartreeFunc); 464 } while (pzomc < -1); 460 } while (pzomc < -1); 465 461 466 // F(EP) rejection 462 // F(EP) rejection 467 G4double XQC = 1.0+tau*(tau-2.0*cosTheta); 463 G4double XQC = 1.0+tau*(tau-2.0*cosTheta); 468 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau 464 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC); 469 if (AF > 0) 465 if (AF > 0) 470 fpzmax = 1.0+AF*0.2; 466 fpzmax = 1.0+AF*0.2; 471 else 467 else 472 fpzmax = 1.0-AF*0.2; 468 fpzmax = 1.0-AF*0.2; 473 fpz = 1.0+AF*std::max(std::min(pzomc,0.2), 469 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2); 474 }while ((fpzmax*G4UniformRand())>fpz); 470 }while ((fpzmax*G4UniformRand())>fpz); 475 471 476 //Energy of the scattered photon 472 //Energy of the scattered photon 477 G4double T = pzomc*pzomc; 473 G4double T = pzomc*pzomc; 478 G4double b1 = 1.0-T*tau*tau; 474 G4double b1 = 1.0-T*tau*tau; 479 G4double b2 = 1.0-T*tau*cosTheta; 475 G4double b2 = 1.0-T*tau*cosTheta; 480 if (pzomc > 0.0) 476 if (pzomc > 0.0) 481 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2 477 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T)))); 482 else 478 else 483 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2 479 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T)))); 484 } //energy < 5 MeV 480 } //energy < 5 MeV 485 481 486 //Ok, the kinematics has been calculated. 482 //Ok, the kinematics has been calculated. 487 G4double sinTheta = std::sqrt(1-cosTheta*cos 483 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta); 488 G4double phi = twopi * G4UniformRand() ; 484 G4double phi = twopi * G4UniformRand() ; 489 G4double dirx = sinTheta * std::cos(phi); 485 G4double dirx = sinTheta * std::cos(phi); 490 G4double diry = sinTheta * std::sin(phi); 486 G4double diry = sinTheta * std::sin(phi); 491 G4double dirz = cosTheta ; 487 G4double dirz = cosTheta ; 492 488 493 // Update G4VParticleChange for the scattere 489 // Update G4VParticleChange for the scattered photon 494 G4ThreeVector photonDirection1(dirx,diry,dir 490 G4ThreeVector photonDirection1(dirx,diry,dirz); 495 photonDirection1.rotateUz(photonDirection0); 491 photonDirection1.rotateUz(photonDirection0); 496 fParticleChange->ProposeMomentumDirection(ph 492 fParticleChange->ProposeMomentumDirection(photonDirection1) ; 497 493 498 G4double photonEnergy1 = epsilon * photonEne 494 G4double photonEnergy1 = epsilon * photonEnergy0; 499 495 500 if (photonEnergy1 > 0.) 496 if (photonEnergy1 > 0.) 501 fParticleChange->SetProposedKineticEnergy( 497 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ; 502 else 498 else 503 { 499 { 504 fParticleChange->SetProposedKineticEnergy( 500 fParticleChange->SetProposedKineticEnergy(0.) ; 505 fParticleChange->ProposeTrackStatus(fStopA 501 fParticleChange->ProposeTrackStatus(fStopAndKill); 506 } 502 } 507 503 508 //Create scattered electron 504 //Create scattered electron 509 G4double diffEnergy = photonEnergy0*(1-epsil 505 G4double diffEnergy = photonEnergy0*(1-epsilon); 510 ionEnergy = (*theTable)[targetOscillator]->G 506 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy(); 511 507 512 G4double Q2 = 508 G4double Q2 = 513 photonEnergy0*photonEnergy0+photonEnergy1* 509 photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta); 514 G4double cosThetaE = 0.; //scattering angle 510 G4double cosThetaE = 0.; //scattering angle for the electron 515 511 516 if (Q2 > 1.0e-12) 512 if (Q2 > 1.0e-12) 517 cosThetaE = (photonEnergy0-photonEnergy1*c 513 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2); 518 else 514 else 519 cosThetaE = 1.0; 515 cosThetaE = 1.0; 520 G4double sinThetaE = std::sqrt(1-cosThetaE*c 516 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE); 521 517 522 //Now, try to handle fluorescence 518 //Now, try to handle fluorescence 523 //Notice: merged levels are indicated with Z 519 //Notice: merged levels are indicated with Z=0 and flag=30 524 G4int shFlag = (*theTable)[targetOscillator] 520 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag(); 525 G4int Z = (G4int) (*theTable)[targetOscillat 521 G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ(); 526 522 527 //initialize here, then check photons create 523 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e- 528 G4double bindingEnergy = 0.*eV; 524 G4double bindingEnergy = 0.*eV; 529 const G4AtomicShell* shell = 0; 525 const G4AtomicShell* shell = 0; 530 526 531 //Real level 527 //Real level 532 if (Z > 0 && shFlag<30) 528 if (Z > 0 && shFlag<30) 533 { 529 { 534 shell = fTransitionManager->Shell(Z,shFl 530 shell = fTransitionManager->Shell(Z,shFlag-1); 535 bindingEnergy = shell->BindingEnergy(); 531 bindingEnergy = shell->BindingEnergy(); 536 } 532 } 537 533 538 G4double ionEnergyInPenelopeDatabase = ionEn 534 G4double ionEnergyInPenelopeDatabase = ionEnergy; 539 //protection against energy non-conservation 535 //protection against energy non-conservation 540 ionEnergy = std::max(bindingEnergy,ionEnergy 536 ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase); 541 537 542 //subtract the excitation energy. If not emi 538 //subtract the excitation energy. If not emitted by fluorescence 543 //the ionization energy is deposited as loca 539 //the ionization energy is deposited as local energy deposition 544 G4double eKineticEnergy = diffEnergy - ionEn 540 G4double eKineticEnergy = diffEnergy - ionEnergy; 545 G4double localEnergyDeposit = ionEnergy; 541 G4double localEnergyDeposit = ionEnergy; 546 G4double energyInFluorescence = 0.; //testin 542 G4double energyInFluorescence = 0.; //testing purposes only 547 G4double energyInAuger = 0; //testing purpos 543 G4double energyInAuger = 0; //testing purposes 548 544 549 if (eKineticEnergy < 0) 545 if (eKineticEnergy < 0) 550 { 546 { 551 //It means that there was some problem/m 547 //It means that there was some problem/mismatch between the two databases. 552 //Try to make it work 548 //Try to make it work 553 //In this case available Energy (diffEne 549 //In this case available Energy (diffEnergy) < ionEnergy 554 //Full residual energy is deposited loca 550 //Full residual energy is deposited locally 555 localEnergyDeposit = diffEnergy; 551 localEnergyDeposit = diffEnergy; 556 eKineticEnergy = 0.0; 552 eKineticEnergy = 0.0; 557 } 553 } 558 554 559 //the local energy deposit is what remains: 555 //the local energy deposit is what remains: part of this may be spent for fluorescence. 560 //Notice: shell might be nullptr (invalid!) 556 //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected 561 //Now, take care of fluorescence, if require 557 //Now, take care of fluorescence, if required 562 if (fAtomDeexcitation && shell) 558 if (fAtomDeexcitation && shell) 563 { 559 { 564 G4int index = couple->GetIndex(); 560 G4int index = couple->GetIndex(); 565 if (fAtomDeexcitation->CheckDeexcitation 561 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) 566 { 562 { 567 size_t nBefore = fvect->size(); 563 size_t nBefore = fvect->size(); 568 fAtomDeexcitation->GenerateParticles(fvect 564 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index); 569 size_t nAfter = fvect->size(); 565 size_t nAfter = fvect->size(); 570 566 571 if (nAfter > nBefore) //actual production 567 if (nAfter > nBefore) //actual production of fluorescence 572 { 568 { 573 for (size_t j=nBefore;j<nAfter;j++) // 569 for (size_t j=nBefore;j<nAfter;j++) //loop on products 574 { 570 { 575 G4double itsEnergy = ((*fvect)[j])->GetK 571 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy(); 576 if (itsEnergy < localEnergyDeposit) // v 572 if (itsEnergy < localEnergyDeposit) // valid secondary, generate it 577 { 573 { 578 localEnergyDeposit -= itsEnergy; 574 localEnergyDeposit -= itsEnergy; 579 if (((*fvect)[j])->GetParticleDefini 575 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition()) 580 energyInFluorescence += itsEnergy; 576 energyInFluorescence += itsEnergy; 581 else if (((*fvect)[j])->GetParticleD 577 else if (((*fvect)[j])->GetParticleDefinition() == 582 G4Electron::Definition()) 578 G4Electron::Definition()) 583 energyInAuger += itsEnergy; 579 energyInAuger += itsEnergy; 584 } 580 } 585 else //invalid secondary: takes more tha 581 else //invalid secondary: takes more than the available energy: delete it 586 { 582 { 587 delete (*fvect)[j]; 583 delete (*fvect)[j]; 588 (*fvect)[j] = nullptr; 584 (*fvect)[j] = nullptr; 589 } 585 } 590 } 586 } 591 } 587 } 592 588 593 } 589 } 594 } 590 } 595 591 596 //Always produce explicitly the electron << 592 >> 593 /* >> 594 if(DeexcitationFlag() && Z > 5 && shellId>0) { >> 595 >> 596 const G4ProductionCutsTable* theCoupleTable= >> 597 G4ProductionCutsTable::GetProductionCutsTable(); >> 598 >> 599 size_t index = couple->GetIndex(); >> 600 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index]; >> 601 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index]; >> 602 >> 603 // Generation of fluorescence >> 604 // Data in EADL are available only for Z > 5 >> 605 // Protection to avoid generating photons in the unphysical case of >> 606 // shell binding energy > photon energy >> 607 if (localEnergyDeposit > cutg || localEnergyDeposit > cute) >> 608 { >> 609 G4DynamicParticle* aPhoton; >> 610 deexcitationManager.SetCutForSecondaryPhotons(cutg); >> 611 deexcitationManager.SetCutForAugerElectrons(cute); >> 612 >> 613 photonVector = deexcitationManager.GenerateParticles(Z,shellId); >> 614 if(photonVector) >> 615 { >> 616 size_t nPhotons = photonVector->size(); >> 617 for (size_t k=0; k<nPhotons; k++) >> 618 { >> 619 aPhoton = (*photonVector)[k]; >> 620 if (aPhoton) >> 621 { >> 622 G4double itsEnergy = aPhoton->GetKineticEnergy(); >> 623 G4bool keepIt = false; >> 624 if (itsEnergy <= localEnergyDeposit) >> 625 { >> 626 //check if good! >> 627 if(aPhoton->GetDefinition() == G4Gamma::Gamma() >> 628 && itsEnergy >= cutg) >> 629 { >> 630 keepIt = true; >> 631 energyInFluorescence += itsEnergy; >> 632 } >> 633 if (aPhoton->GetDefinition() == G4Electron::Electron() && >> 634 itsEnergy >= cute) >> 635 { >> 636 energyInAuger += itsEnergy; >> 637 keepIt = true; >> 638 } >> 639 } >> 640 //good secondary, register it >> 641 if (keepIt) >> 642 { >> 643 localEnergyDeposit -= itsEnergy; >> 644 fvect->push_back(aPhoton); >> 645 } >> 646 else >> 647 { >> 648 delete aPhoton; >> 649 (*photonVector)[k] = 0; >> 650 } >> 651 } >> 652 } >> 653 delete photonVector; >> 654 } >> 655 } >> 656 } >> 657 */ >> 658 >> 659 >> 660 //Always produce explicitely the electron 597 G4DynamicParticle* electron = 0; 661 G4DynamicParticle* electron = 0; 598 662 599 G4double xEl = sinThetaE * std::cos(phi+pi); 663 G4double xEl = sinThetaE * std::cos(phi+pi); 600 G4double yEl = sinThetaE * std::sin(phi+pi); 664 G4double yEl = sinThetaE * std::sin(phi+pi); 601 G4double zEl = cosThetaE; 665 G4double zEl = cosThetaE; 602 G4ThreeVector eDirection(xEl,yEl,zEl); //ele 666 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction 603 eDirection.rotateUz(photonDirection0); 667 eDirection.rotateUz(photonDirection0); 604 electron = new G4DynamicParticle (G4Electron 668 electron = new G4DynamicParticle (G4Electron::Electron(), 605 eDirection,eKineticEnergy) ; 669 eDirection,eKineticEnergy) ; 606 fvect->push_back(electron); 670 fvect->push_back(electron); 607 671 >> 672 608 if (localEnergyDeposit < 0) //Should not be: 673 if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning) 609 { 674 { 610 G4Exception("G4PenelopeComptonModel::Sam 675 G4Exception("G4PenelopeComptonModel::SampleSecondaries()", 611 "em2099",JustWarning,"WARNING: Negative 676 "em2099",JustWarning,"WARNING: Negative local energy deposit"); 612 localEnergyDeposit=0.; 677 localEnergyDeposit=0.; 613 } 678 } 614 fParticleChange->ProposeLocalEnergyDeposit(l 679 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit); 615 680 616 G4double electronEnergy = 0.; 681 G4double electronEnergy = 0.; 617 if (electron) 682 if (electron) 618 electronEnergy = eKineticEnergy; 683 electronEnergy = eKineticEnergy; 619 if (fVerboseLevel > 1) << 684 if (verboseLevel > 1) 620 { 685 { 621 G4cout << "----------------------------- 686 G4cout << "-----------------------------------------------------------" << G4endl; 622 G4cout << "Energy balance from G4Penelop 687 G4cout << "Energy balance from G4PenelopeCompton" << G4endl; 623 G4cout << "Incoming photon energy: " << 688 G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl; 624 G4cout << "----------------------------- 689 G4cout << "-----------------------------------------------------------" << G4endl; 625 G4cout << "Scattered photon: " << photon 690 G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl; 626 G4cout << "Scattered electron " << elect 691 G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl; 627 if (energyInFluorescence) 692 if (energyInFluorescence) 628 G4cout << "Fluorescence x-rays: " << energyI 693 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl; 629 if (energyInAuger) 694 if (energyInAuger) 630 G4cout << "Auger electrons: " << energyInAug 695 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl; 631 G4cout << "Local energy deposit " << loc 696 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; 632 G4cout << "Total final state: " << (phot 697 G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+ 633 localEnergyDeposit+energyInAuger)/ 698 localEnergyDeposit+energyInAuger)/keV << 634 " keV" << G4endl; 699 " keV" << G4endl; 635 G4cout << "----------------------------- 700 G4cout << "-----------------------------------------------------------" << G4endl; 636 } 701 } 637 if (fVerboseLevel > 0) << 702 if (verboseLevel > 0) 638 { 703 { 639 G4double energyDiff = std::fabs(photonEn 704 G4double energyDiff = std::fabs(photonEnergy1+ 640 electronEnergy+energyInFluoresce 705 electronEnergy+energyInFluorescence+ 641 localEnergyDeposit+energyInAuger 706 localEnergyDeposit+energyInAuger-photonEnergy0); 642 if (energyDiff > 0.05*keV) 707 if (energyDiff > 0.05*keV) 643 G4cout << "Warning from G4PenelopeCompton: p 708 G4cout << "Warning from G4PenelopeCompton: problem with energy conservation: " << 644 (photonEnergy1+electronEnergy+energyInFluo 709 (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/keV << 645 " keV (final) vs. " << 710 " keV (final) vs. " << 646 photonEnergy0/keV << " keV (initial)" << G 711 photonEnergy0/keV << " keV (initial)" << G4endl; 647 } 712 } 648 } 713 } 649 714 650 //....oooOO0OOooo........oooOO0OOooo........oo 715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 651 716 652 G4double G4PenelopeComptonModel::DifferentialC 717 G4double G4PenelopeComptonModel::DifferentialCrossSection(G4double cosTheta,G4double energy, 653 G4PenelopeOscillator* osc) 718 G4PenelopeOscillator* osc) 654 { 719 { 655 // 720 // 656 // Penelope model v2008. Single differential 721 // Penelope model v2008. Single differential cross section *per electron* 657 // for photon Compton scattering by 722 // for photon Compton scattering by 658 // electrons in the given atomic oscillator, 723 // electrons in the given atomic oscillator, differential in the direction of the 659 // scattering photon. This is in units of pi 724 // scattering photon. This is in units of pi*classic_electr_radius**2 660 // 725 // 661 // D. Brusa et al., Nucl. Instrum. Meth. A 3 726 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167 662 // The parametrization includes the J(p) dis 727 // The parametrization includes the J(p) distribution profiles for the atomic shells, 663 // that are tabulated from Hartree-Fock calc 728 // that are tabulated from Hartree-Fock calculations 664 // from F. Biggs et al., At. Data Nucl. Data 729 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201 665 // 730 // 666 G4double ionEnergy = osc->GetIonisationEnerg 731 G4double ionEnergy = osc->GetIonisationEnergy(); 667 G4double harFunc = osc->GetHartreeFactor(); 732 G4double harFunc = osc->GetHartreeFactor(); 668 733 669 static const G4double k2 = std::sqrt(2.); 734 static const G4double k2 = std::sqrt(2.); 670 static const G4double k1 = 1./k2; 735 static const G4double k1 = 1./k2; 671 736 672 if (energy < ionEnergy) 737 if (energy < ionEnergy) 673 return 0; 738 return 0; 674 739 675 //energy of the Compton line 740 //energy of the Compton line 676 G4double cdt1 = 1.0-cosTheta; 741 G4double cdt1 = 1.0-cosTheta; 677 G4double EOEC = 1.0+(energy/electron_mass_c2 742 G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1; 678 G4double ECOE = 1.0/EOEC; 743 G4double ECOE = 1.0/EOEC; 679 744 680 //Incoherent scattering function (analytical 745 //Incoherent scattering function (analytical profile) 681 G4double aux = energy*(energy-ionEnergy)*cdt 746 G4double aux = energy*(energy-ionEnergy)*cdt1; 682 G4double Pzimax = 747 G4double Pzimax = 683 (aux - electron_mass_c2*ionEnergy)/(electr 748 (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy)); 684 G4double sia = 0.0; 749 G4double sia = 0.0; 685 G4double x = harFunc*Pzimax; 750 G4double x = harFunc*Pzimax; 686 if (x > 0) 751 if (x > 0) 687 sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x 752 sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x)); 688 else 753 else 689 sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x)); 754 sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x)); 690 755 691 //1st order correction, integral of Pz times 756 //1st order correction, integral of Pz times the Compton profile. 692 //Calculated approximately using a free-elec 757 //Calculated approximately using a free-electron gas profile 693 G4double pf = 3.0/(4.0*harFunc); 758 G4double pf = 3.0/(4.0*harFunc); 694 if (std::fabs(Pzimax) < pf) 759 if (std::fabs(Pzimax) < pf) 695 { 760 { 696 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE* 761 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta; 697 G4double p2 = Pzimax*Pzimax; 762 G4double p2 = Pzimax*Pzimax; 698 G4double dspz = std::sqrt(QCOE2)* 763 G4double dspz = std::sqrt(QCOE2)* 699 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc 764 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc 700 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf)); 765 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf)); 701 sia += std::max(dspz,-1.0*sia); 766 sia += std::max(dspz,-1.0*sia); 702 } 767 } 703 768 704 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosThe 769 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta; 705 770 706 //Differential cross section (per electron, 771 //Differential cross section (per electron, in units of pi*classic_electr_radius**2) 707 G4double diffCS = ECOE*ECOE*XKN*sia; 772 G4double diffCS = ECOE*ECOE*XKN*sia; 708 773 709 return diffCS; 774 return diffCS; 710 } 775 } 711 776 712 //....oooOO0OOooo........oooOO0OOooo........oo 777 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 713 778 714 G4double G4PenelopeComptonModel::OscillatorTot 779 G4double G4PenelopeComptonModel::OscillatorTotalCrossSection(G4double energy,G4PenelopeOscillator* osc) 715 { 780 { 716 //Total cross section (integrated) for the g 781 //Total cross section (integrated) for the given oscillator in units of 717 //pi*classic_electr_radius^2 782 //pi*classic_electr_radius^2 718 783 719 //Integrate differential cross section for e 784 //Integrate differential cross section for each oscillator 720 G4double stre = osc->GetOscillatorStrength() 785 G4double stre = osc->GetOscillatorStrength(); 721 786 722 // here one uses the using the 20-point 787 // here one uses the using the 20-point 723 // Gauss quadrature method with an adaptive 788 // Gauss quadrature method with an adaptive bipartition scheme 724 const G4int npoints=10; 789 const G4int npoints=10; 725 const G4int ncallsmax=20000; 790 const G4int ncallsmax=20000; 726 const G4int nst=256; 791 const G4int nst=256; 727 static const G4double Abscissas[10] = {7.652 792 static const G4double Abscissas[10] = {7.652651133497334e-02,2.2778585114164508e-01,3.7370608871541956e-01, 728 5.1086700195082710e-01,6.36053680726 793 5.1086700195082710e-01,6.3605368072651503e-01,7.4633190646015079e-01, 729 8.3911697182221882e-01,9.12234428251 794 8.3911697182221882e-01,9.1223442825132591e-01,9.6397192727791379e-01, 730 9.9312859918509492e-01}; 795 9.9312859918509492e-01}; 731 static const G4double Weights[10] = {1.52753 796 static const G4double Weights[10] = {1.5275338713072585e-01,1.4917298647260375e-01,1.4209610931838205e-01, 732 1.3168863844917663e-01,1.1819453196151 797 1.3168863844917663e-01,1.1819453196151842e-01,1.0193011981724044e-01, 733 8.3276741576704749e-02,6.2672048334109 798 8.3276741576704749e-02,6.2672048334109064e-02,4.0601429800386941e-02, 734 1.7614007139152118e-02}; 799 1.7614007139152118e-02}; 735 800 736 G4double MaxError = 1e-5; 801 G4double MaxError = 1e-5; 737 //Error control 802 //Error control 738 G4double Ctol = std::min(std::max(MaxError,1 803 G4double Ctol = std::min(std::max(MaxError,1e-13),1e-02); 739 G4double Ptol = 0.01*Ctol; 804 G4double Ptol = 0.01*Ctol; 740 G4double Err=1e35; 805 G4double Err=1e35; 741 806 742 //Gauss integration from -1 to 1 807 //Gauss integration from -1 to 1 743 G4double LowPoint = -1.0; 808 G4double LowPoint = -1.0; 744 G4double HighPoint = 1.0; 809 G4double HighPoint = 1.0; 745 810 746 G4double h=HighPoint-LowPoint; 811 G4double h=HighPoint-LowPoint; 747 G4double sumga=0.0; 812 G4double sumga=0.0; 748 G4double a=0.5*(HighPoint-LowPoint); 813 G4double a=0.5*(HighPoint-LowPoint); 749 G4double b=0.5*(HighPoint+LowPoint); 814 G4double b=0.5*(HighPoint+LowPoint); 750 G4double c=a*Abscissas[0]; 815 G4double c=a*Abscissas[0]; 751 G4double d= Weights[0]* 816 G4double d= Weights[0]* 752 (DifferentialCrossSection(b+c,energy,osc)+ 817 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); 753 for (G4int i=2;i<=npoints;i++) 818 for (G4int i=2;i<=npoints;i++) 754 { 819 { 755 c=a*Abscissas[i-1]; 820 c=a*Abscissas[i-1]; 756 d += Weights[i-1]* 821 d += Weights[i-1]* 757 (DifferentialCrossSection(b+c,energy,osc)+Di 822 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); 758 } 823 } 759 G4int icall = 2*npoints; 824 G4int icall = 2*npoints; 760 G4int LH=1; 825 G4int LH=1; 761 G4double S[nst],x[nst],sn[nst],xrn[nst]; 826 G4double S[nst],x[nst],sn[nst],xrn[nst]; 762 S[0]=d*a; 827 S[0]=d*a; 763 x[0]=LowPoint; 828 x[0]=LowPoint; 764 829 765 G4bool loopAgain = true; 830 G4bool loopAgain = true; 766 831 767 //Adaptive bipartition scheme 832 //Adaptive bipartition scheme 768 do{ 833 do{ 769 G4double h0=h; 834 G4double h0=h; 770 h=0.5*h; //bipartition 835 h=0.5*h; //bipartition 771 G4double sumr=0; 836 G4double sumr=0; 772 G4int LHN=0; 837 G4int LHN=0; 773 G4double si,xa,xb,xc; 838 G4double si,xa,xb,xc; 774 for (G4int i=1;i<=LH;i++){ 839 for (G4int i=1;i<=LH;i++){ 775 si=S[i-1]; 840 si=S[i-1]; 776 xa=x[i-1]; 841 xa=x[i-1]; 777 xb=xa+h; 842 xb=xa+h; 778 xc=xa+h0; 843 xc=xa+h0; 779 a=0.5*(xb-xa); 844 a=0.5*(xb-xa); 780 b=0.5*(xb+xa); 845 b=0.5*(xb+xa); 781 c=a*Abscissas[0]; 846 c=a*Abscissas[0]; 782 G4double dLocal = Weights[0]* 847 G4double dLocal = Weights[0]* 783 (DifferentialCrossSection(b+c,energy,osc)+Di 848 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); 784 849 785 for (G4int j=1;j<npoints;j++) 850 for (G4int j=1;j<npoints;j++) 786 { 851 { 787 c=a*Abscissas[j]; 852 c=a*Abscissas[j]; 788 dLocal += Weights[j]* 853 dLocal += Weights[j]* 789 (DifferentialCrossSection(b+c,energy,osc 854 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); 790 } 855 } 791 G4double s1=dLocal*a; 856 G4double s1=dLocal*a; 792 a=0.5*(xc-xb); 857 a=0.5*(xc-xb); 793 b=0.5*(xc+xb); 858 b=0.5*(xc+xb); 794 c=a*Abscissas[0]; 859 c=a*Abscissas[0]; 795 dLocal=Weights[0]* 860 dLocal=Weights[0]* 796 (DifferentialCrossSection(b+c,energy,osc)+Di 861 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); 797 862 798 for (G4int j=1;j<npoints;j++) 863 for (G4int j=1;j<npoints;j++) 799 { 864 { 800 c=a*Abscissas[j]; 865 c=a*Abscissas[j]; 801 dLocal += Weights[j]* 866 dLocal += Weights[j]* 802 (DifferentialCrossSection(b+c,energy,osc 867 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc)); 803 } 868 } 804 G4double s2=dLocal*a; 869 G4double s2=dLocal*a; 805 icall=icall+4*npoints; 870 icall=icall+4*npoints; 806 G4double s12=s1+s2; 871 G4double s12=s1+s2; 807 if (std::abs(s12-si)<std::max(Ptol*std:: 872 if (std::abs(s12-si)<std::max(Ptol*std::abs(s12),1e-35)) 808 sumga += s12; 873 sumga += s12; 809 else 874 else 810 { 875 { 811 sumr += s12; 876 sumr += s12; 812 LHN += 2; 877 LHN += 2; 813 sn[LHN-1]=s2; 878 sn[LHN-1]=s2; 814 xrn[LHN-1]=xb; 879 xrn[LHN-1]=xb; 815 sn[LHN-2]=s1; 880 sn[LHN-2]=s1; 816 xrn[LHN-2]=xa; 881 xrn[LHN-2]=xa; 817 } 882 } 818 883 819 if (icall>ncallsmax || LHN>nst) 884 if (icall>ncallsmax || LHN>nst) 820 { 885 { 821 G4cout << "G4PenelopeComptonModel: " << G4 886 G4cout << "G4PenelopeComptonModel: " << G4endl; 822 G4cout << "LowPoint: " << LowPoint << ", H 887 G4cout << "LowPoint: " << LowPoint << ", High Point: " << HighPoint << G4endl; 823 G4cout << "Tolerance: " << MaxError << G4e 888 G4cout << "Tolerance: " << MaxError << G4endl; 824 G4cout << "Calls: " << icall << ", Integra 889 G4cout << "Calls: " << icall << ", Integral: " << sumga << ", Error: " << Err << G4endl; 825 G4cout << "Number of open subintervals: " 890 G4cout << "Number of open subintervals: " << LHN << G4endl; 826 G4cout << "WARNING: the required accuracy 891 G4cout << "WARNING: the required accuracy has not been attained" << G4endl; 827 loopAgain = false; 892 loopAgain = false; 828 } 893 } 829 } 894 } 830 Err=std::abs(sumr)/std::max(std::abs(sumr+ 895 Err=std::abs(sumr)/std::max(std::abs(sumr+sumga),1e-35); 831 if (Err < Ctol || LHN == 0) 896 if (Err < Ctol || LHN == 0) 832 loopAgain = false; //end of cycle 897 loopAgain = false; //end of cycle 833 LH=LHN; 898 LH=LHN; 834 for (G4int i=0;i<LH;i++) 899 for (G4int i=0;i<LH;i++) 835 { 900 { 836 S[i]=sn[i]; 901 S[i]=sn[i]; 837 x[i]=xrn[i]; 902 x[i]=xrn[i]; 838 } 903 } 839 }while(Ctol < 1.0 && loopAgain); 904 }while(Ctol < 1.0 && loopAgain); 840 905 >> 906 841 G4double xs = stre*sumga; 907 G4double xs = stre*sumga; 842 908 843 return xs; 909 return xs; 844 } 910 } 845 911 846 //....oooOO0OOooo........oooOO0OOooo........oo 912 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 847 913 848 G4double G4PenelopeComptonModel::KleinNishinaC 914 G4double G4PenelopeComptonModel::KleinNishinaCrossSection(G4double energy, 849 const G4Material* material) 915 const G4Material* material) 850 { 916 { 851 // use Klein-Nishina formula 917 // use Klein-Nishina formula 852 // total cross section in units of pi*classi 918 // total cross section in units of pi*classic_electr_radius^2 >> 919 853 G4double cs = 0; 920 G4double cs = 0; 854 921 855 G4double ek =energy/electron_mass_c2; 922 G4double ek =energy/electron_mass_c2; 856 G4double eks = ek*ek; 923 G4double eks = ek*ek; 857 G4double ek2 = 1.0+ek+ek; 924 G4double ek2 = 1.0+ek+ek; 858 G4double ek1 = eks-ek2-1.0; 925 G4double ek1 = eks-ek2-1.0; 859 926 860 G4double t0 = 1.0/ek2; 927 G4double t0 = 1.0/ek2; 861 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*G4Lo << 928 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*std::log(t0)-(1.0/t0); 862 929 863 G4PenelopeOscillatorTable* theTable = fOscMa << 930 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); 864 931 865 for (size_t i=0;i<theTable->size();i++) 932 for (size_t i=0;i<theTable->size();i++) 866 { 933 { 867 G4PenelopeOscillator* theOsc = (*theTabl 934 G4PenelopeOscillator* theOsc = (*theTable)[i]; 868 G4double ionEnergy = theOsc->GetIonisati 935 G4double ionEnergy = theOsc->GetIonisationEnergy(); 869 G4double tau=(energy-ionEnergy)/energy; 936 G4double tau=(energy-ionEnergy)/energy; 870 if (tau > t0) 937 if (tau > t0) 871 { 938 { 872 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1 << 939 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*std::log(tau)-(1.0/tau); 873 G4double stre = theOsc->GetOscillatorStren 940 G4double stre = theOsc->GetOscillatorStrength(); 874 941 875 cs += stre*(csu-csl); 942 cs += stre*(csu-csl); 876 } 943 } 877 } 944 } >> 945 878 cs /= (ek*eks); 946 cs /= (ek*eks); 879 947 880 return cs; 948 return cs; 881 949 882 } 950 } 883 951 884 //....oooOO0OOooo........oooOO0OOooo........oo 952 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 885 953 886 void G4PenelopeComptonModel::SetParticle(const 954 void G4PenelopeComptonModel::SetParticle(const G4ParticleDefinition* p) 887 { 955 { 888 if(!fParticle) { 956 if(!fParticle) { 889 fParticle = p; 957 fParticle = p; 890 } 958 } 891 } 959 } 892 960