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