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: G4KleinNishinaCompton.cc 93362 2015-10-19 13:45:19Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class file 30 // GEANT4 Class file 30 // 31 // 31 // 32 // 32 // File name: G4KleinNishinaCompton 33 // File name: G4KleinNishinaCompton 33 // 34 // 34 // Author: Vladimir Ivanchenko on base 35 // Author: Vladimir Ivanchenko on base of Michel Maire code 35 // 36 // 36 // Creation date: 15.03.2005 37 // Creation date: 15.03.2005 37 // 38 // 38 // Modifications: 39 // Modifications: 39 // 18-04-05 Use G4ParticleChangeForGamma (V.Iv 40 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivantchenko) 40 // 27-03-06 Remove upper limit of cross sectio 41 // 27-03-06 Remove upper limit of cross section (V.Ivantchenko) 41 // 42 // 42 // Class Description: 43 // Class Description: 43 // 44 // 44 // ------------------------------------------- 45 // ------------------------------------------------------------------- 45 // 46 // 46 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 48 49 49 #include "G4KleinNishinaCompton.hh" 50 #include "G4KleinNishinaCompton.hh" 50 #include "G4PhysicalConstants.hh" 51 #include "G4PhysicalConstants.hh" 51 #include "G4SystemOfUnits.hh" 52 #include "G4SystemOfUnits.hh" 52 #include "G4Electron.hh" 53 #include "G4Electron.hh" 53 #include "G4Gamma.hh" 54 #include "G4Gamma.hh" 54 #include "Randomize.hh" 55 #include "Randomize.hh" 55 #include "G4DataVector.hh" 56 #include "G4DataVector.hh" 56 #include "G4ParticleChangeForGamma.hh" 57 #include "G4ParticleChangeForGamma.hh" 57 #include "G4Log.hh" 58 #include "G4Log.hh" 58 #include "G4Exp.hh" 59 #include "G4Exp.hh" 59 60 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 61 62 62 using namespace std; 63 using namespace std; 63 64 >> 65 static const G4double >> 66 d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLHEP::barn, >> 67 d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLHEP::barn, >> 68 e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLHEP::barn, >> 69 e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLHEP::barn, >> 70 f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLHEP::barn, >> 71 f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLHEP::barn; >> 72 static const G4int nlooplim = 1000; >> 73 64 G4KleinNishinaCompton::G4KleinNishinaCompton(c 74 G4KleinNishinaCompton::G4KleinNishinaCompton(const G4ParticleDefinition*, 65 c 75 const G4String& nam) 66 : G4VEmModel(nam) 76 : G4VEmModel(nam) 67 { 77 { 68 theGamma = G4Gamma::Gamma(); 78 theGamma = G4Gamma::Gamma(); 69 theElectron = G4Electron::Electron(); 79 theElectron = G4Electron::Electron(); 70 lowestSecondaryEnergy = 100.0*eV; 80 lowestSecondaryEnergy = 100.0*eV; 71 fParticleChange = nullptr; 81 fParticleChange = nullptr; 72 } 82 } 73 83 74 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 85 76 G4KleinNishinaCompton::~G4KleinNishinaCompton( << 86 G4KleinNishinaCompton::~G4KleinNishinaCompton() >> 87 {} 77 88 78 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 79 90 80 void G4KleinNishinaCompton::Initialise(const G 91 void G4KleinNishinaCompton::Initialise(const G4ParticleDefinition* p, 81 const G 92 const G4DataVector& cuts) 82 { 93 { 83 if(IsMaster()) { InitialiseElementSelectors( 94 if(IsMaster()) { InitialiseElementSelectors(p, cuts); } 84 if(nullptr == fParticleChange) { 95 if(nullptr == fParticleChange) { 85 fParticleChange = GetParticleChangeForGamm 96 fParticleChange = GetParticleChangeForGamma(); 86 } 97 } 87 } 98 } 88 99 89 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 90 101 91 void G4KleinNishinaCompton::InitialiseLocal(co 102 void G4KleinNishinaCompton::InitialiseLocal(const G4ParticleDefinition*, 92 G4 103 G4VEmModel* masterModel) 93 { 104 { 94 SetElementSelectors(masterModel->GetElementS 105 SetElementSelectors(masterModel->GetElementSelectors()); 95 } 106 } 96 107 97 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 109 99 G4double G4KleinNishinaCompton::ComputeCrossSe 110 G4double G4KleinNishinaCompton::ComputeCrossSectionPerAtom( 100 const G 111 const G4ParticleDefinition*, 101 G 112 G4double GammaEnergy, 102 G 113 G4double Z, G4double, 103 G 114 G4double, G4double) 104 { 115 { 105 G4double xSection = 0.0 ; 116 G4double xSection = 0.0 ; 106 if (GammaEnergy <= LowEnergyLimit()) { retur 117 if (GammaEnergy <= LowEnergyLimit()) { return xSection; } 107 118 108 static const G4double a = 20.0 , b = 230.0 , 119 static const G4double a = 20.0 , b = 230.0 , c = 440.0; 109 << 110 static const G4double << 111 d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLH << 112 d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLH << 113 e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLH << 114 e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLH << 115 f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLH << 116 f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLH << 117 120 118 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = 121 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z), 119 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = 122 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z); 120 123 121 G4double T0 = 15.0*keV; 124 G4double T0 = 15.0*keV; 122 if (Z < 1.5) { T0 = 40.0*keV; } 125 if (Z < 1.5) { T0 = 40.0*keV; } 123 126 124 G4double X = max(GammaEnergy, T0) / electr 127 G4double X = max(GammaEnergy, T0) / electron_mass_c2; 125 xSection = p1Z*G4Log(1.+2.*X)/X 128 xSection = p1Z*G4Log(1.+2.*X)/X 126 + (p2Z + p3Z*X + p4Z*X*X)/(1. + 129 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X); 127 130 128 // modification for low energy. (special ca 131 // modification for low energy. (special case for Hydrogen) 129 if (GammaEnergy < T0) { 132 if (GammaEnergy < T0) { 130 static const G4double dT0 = keV; 133 static const G4double dT0 = keV; 131 X = (T0+dT0) / electron_mass_c2 ; 134 X = (T0+dT0) / electron_mass_c2 ; 132 G4double sigma = p1Z*G4Log(1.+2*X)/X 135 G4double sigma = p1Z*G4Log(1.+2*X)/X 133 + (p2Z + p3Z*X + p4Z*X*X)/ 136 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X); 134 G4double c1 = -T0*(sigma-xSection)/(xSec 137 G4double c1 = -T0*(sigma-xSection)/(xSection*dT0); 135 G4double c2 = 0.150; 138 G4double c2 = 0.150; 136 if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); 139 if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); } 137 G4double y = G4Log(GammaEnergy/T0); 140 G4double y = G4Log(GammaEnergy/T0); 138 xSection *= G4Exp(-y*(c1+c2*y)); 141 xSection *= G4Exp(-y*(c1+c2*y)); 139 } 142 } 140 // G4cout<<"e= "<< GammaEnergy<<" Z= "<<Z<<" 143 // G4cout<<"e= "<< GammaEnergy<<" Z= "<<Z<<" cross= " << xSection << G4endl; 141 return xSection; 144 return xSection; 142 } 145 } 143 146 144 //....oooOO0OOooo........oooOO0OOooo........oo 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 145 148 146 void G4KleinNishinaCompton::SampleSecondaries( 149 void G4KleinNishinaCompton::SampleSecondaries( 147 std::vector<G4Dyna 150 std::vector<G4DynamicParticle*>* fvect, 148 const G4MaterialCu 151 const G4MaterialCutsCouple*, 149 const G4DynamicPar 152 const G4DynamicParticle* aDynamicGamma, 150 G4double, 153 G4double, 151 G4double) 154 G4double) 152 { 155 { 153 // The scattered gamma energy is sampled acc 156 // The scattered gamma energy is sampled according to Klein - Nishina formula. 154 // The random number techniques of Butcher & 157 // The random number techniques of Butcher & Messel are used 155 // (Nuc Phys 20(1960),15). 158 // (Nuc Phys 20(1960),15). 156 // Note : Effects due to binding of atomic e 159 // Note : Effects due to binding of atomic electrons are negliged. 157 160 158 G4double gamEnergy0 = aDynamicGamma->GetKine 161 G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy(); 159 162 160 // do nothing below the threshold 163 // do nothing below the threshold 161 if(gamEnergy0 <= LowEnergyLimit()) { return; 164 if(gamEnergy0 <= LowEnergyLimit()) { return; } 162 165 163 G4double E0_m = gamEnergy0 / electron_mass_c 166 G4double E0_m = gamEnergy0 / electron_mass_c2 ; 164 167 165 G4ThreeVector gamDirection0 = aDynamicGamma- 168 G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection(); 166 169 167 // 170 // 168 // sample the energy rate of the scattered g 171 // sample the energy rate of the scattered gamma 169 // 172 // 170 173 171 G4double epsilon, epsilonsq, onecost, sint2, 174 G4double epsilon, epsilonsq, onecost, sint2, greject ; 172 175 173 G4double eps0 = 1./(1. + 2.*E0_m); 176 G4double eps0 = 1./(1. + 2.*E0_m); 174 G4double epsilon0sq = eps0*eps0; 177 G4double epsilon0sq = eps0*eps0; 175 G4double alpha1 = - G4Log(eps0); 178 G4double alpha1 = - G4Log(eps0); 176 G4double alpha2 = alpha1 + 0.5*(1.- epsi 179 G4double alpha2 = alpha1 + 0.5*(1.- epsilon0sq); 177 180 178 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 181 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine(); 179 G4double rndm[3]; 182 G4double rndm[3]; 180 183 181 static const G4int nlooplim = 1000; << 182 G4int nloop = 0; 184 G4int nloop = 0; 183 do { 185 do { 184 ++nloop; 186 ++nloop; 185 // false interaction if too many iteration 187 // false interaction if too many iterations 186 if(nloop > nlooplim) { return; } 188 if(nloop > nlooplim) { return; } 187 189 188 // 3 random numbers to sample scattering 190 // 3 random numbers to sample scattering 189 rndmEngineMod->flatArray(3, rndm); 191 rndmEngineMod->flatArray(3, rndm); 190 192 191 if ( alpha1 > alpha2*rndm[0] ) { 193 if ( alpha1 > alpha2*rndm[0] ) { 192 epsilon = G4Exp(-alpha1*rndm[1]); // 194 epsilon = G4Exp(-alpha1*rndm[1]); // eps0**r 193 epsilonsq = epsilon*epsilon; 195 epsilonsq = epsilon*epsilon; 194 196 195 } else { 197 } else { 196 epsilonsq = epsilon0sq + (1.- epsilon0sq 198 epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndm[1]; 197 epsilon = sqrt(epsilonsq); 199 epsilon = sqrt(epsilonsq); 198 }; 200 }; 199 201 200 onecost = (1.- epsilon)/(epsilon*E0_m); 202 onecost = (1.- epsilon)/(epsilon*E0_m); 201 sint2 = onecost*(2.-onecost); 203 sint2 = onecost*(2.-onecost); 202 greject = 1. - epsilon*sint2/(1.+ epsilons 204 greject = 1. - epsilon*sint2/(1.+ epsilonsq); 203 205 204 // Loop checking, 03-Aug-2015, Vladimir Iv 206 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 205 } while (greject < rndm[2]); 207 } while (greject < rndm[2]); 206 208 207 // 209 // 208 // scattered gamma angles. ( Z - axis along 210 // scattered gamma angles. ( Z - axis along the parent gamma) 209 // 211 // 210 212 211 if(sint2 < 0.0) { sint2 = 0.0; } 213 if(sint2 < 0.0) { sint2 = 0.0; } 212 G4double cosTeta = 1. - onecost; 214 G4double cosTeta = 1. - onecost; 213 G4double sinTeta = sqrt (sint2); 215 G4double sinTeta = sqrt (sint2); 214 G4double Phi = twopi * rndmEngineMod->fl 216 G4double Phi = twopi * rndmEngineMod->flat(); 215 217 216 // 218 // 217 // update G4VParticleChange for the scattere 219 // update G4VParticleChange for the scattered gamma 218 // 220 // 219 221 220 G4ThreeVector gamDirection1(sinTeta*cos(Phi) 222 G4ThreeVector gamDirection1(sinTeta*cos(Phi), sinTeta*sin(Phi), cosTeta); 221 gamDirection1.rotateUz(gamDirection0); 223 gamDirection1.rotateUz(gamDirection0); 222 G4double gamEnergy1 = epsilon*gamEnergy0; 224 G4double gamEnergy1 = epsilon*gamEnergy0; 223 G4double edep = 0.0; 225 G4double edep = 0.0; 224 if(gamEnergy1 > lowestSecondaryEnergy) { 226 if(gamEnergy1 > lowestSecondaryEnergy) { 225 fParticleChange->ProposeMomentumDirection( 227 fParticleChange->ProposeMomentumDirection(gamDirection1); 226 fParticleChange->SetProposedKineticEnergy( 228 fParticleChange->SetProposedKineticEnergy(gamEnergy1); 227 } else { 229 } else { 228 fParticleChange->ProposeTrackStatus(fStopA 230 fParticleChange->ProposeTrackStatus(fStopAndKill); 229 fParticleChange->SetProposedKineticEnergy( 231 fParticleChange->SetProposedKineticEnergy(0.0); 230 edep = gamEnergy1; 232 edep = gamEnergy1; 231 } 233 } 232 234 233 // 235 // 234 // kinematic of the scattered electron 236 // kinematic of the scattered electron 235 // 237 // 236 238 237 G4double eKinEnergy = gamEnergy0 - gamEnergy 239 G4double eKinEnergy = gamEnergy0 - gamEnergy1; 238 240 239 if(eKinEnergy > lowestSecondaryEnergy) { 241 if(eKinEnergy > lowestSecondaryEnergy) { 240 G4ThreeVector eDirection = gamEnergy0*gamD 242 G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1; 241 eDirection = eDirection.unit(); 243 eDirection = eDirection.unit(); 242 244 243 // create G4DynamicParticle object for the 245 // create G4DynamicParticle object for the electron. 244 auto dp = new G4DynamicParticle(theElectro << 246 G4DynamicParticle* dp = new G4DynamicParticle(theElectron,eDirection,eKinEnergy); 245 fvect->push_back(dp); 247 fvect->push_back(dp); 246 } else { 248 } else { 247 edep += eKinEnergy; 249 edep += eKinEnergy; 248 } 250 } 249 // energy balance 251 // energy balance 250 if(edep > 0.0) { 252 if(edep > 0.0) { 251 fParticleChange->ProposeLocalEnergyDeposit 253 fParticleChange->ProposeLocalEnergyDeposit(edep); 252 } 254 } 253 } 255 } 254 256 255 //....oooOO0OOooo........oooOO0OOooo........oo 257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 256 258 257 259 258 260