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 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4KleinNishinaModel 32 // File name: G4KleinNishinaModel 33 // 33 // 34 // Author: Vladimir Ivanchenko on base 34 // Author: Vladimir Ivanchenko on base of G4KleinNishinaCompton 35 // 35 // 36 // Creation date: 13.06.2010 36 // Creation date: 13.06.2010 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 39 // 40 // Class Description: 40 // Class Description: 41 // 41 // 42 // ------------------------------------------- 42 // ------------------------------------------------------------------- 43 // 43 // 44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 46 47 #include "G4KleinNishinaModel.hh" 47 #include "G4KleinNishinaModel.hh" 48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4Electron.hh" 50 #include "G4Electron.hh" 51 #include "G4Gamma.hh" 51 #include "G4Gamma.hh" 52 #include "Randomize.hh" 52 #include "Randomize.hh" 53 #include "G4RandomDirection.hh" 53 #include "G4RandomDirection.hh" 54 #include "G4DataVector.hh" 54 #include "G4DataVector.hh" 55 #include "G4ParticleChangeForGamma.hh" 55 #include "G4ParticleChangeForGamma.hh" 56 #include "G4VAtomDeexcitation.hh" 56 #include "G4VAtomDeexcitation.hh" 57 #include "G4AtomicShells.hh" 57 #include "G4AtomicShells.hh" 58 #include "G4LossTableManager.hh" 58 #include "G4LossTableManager.hh" 59 #include "G4Log.hh" 59 #include "G4Log.hh" 60 #include "G4Exp.hh" 60 #include "G4Exp.hh" 61 61 62 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 63 63 64 using namespace std; 64 using namespace std; 65 65 66 G4KleinNishinaModel::G4KleinNishinaModel(const 66 G4KleinNishinaModel::G4KleinNishinaModel(const G4String& nam) 67 : G4VEmModel(nam), 67 : G4VEmModel(nam), 68 lv1(0.,0.,0.,0.), 68 lv1(0.,0.,0.,0.), 69 lv2(0.,0.,0.,0.), 69 lv2(0.,0.,0.,0.), 70 bst(0.,0.,0.) 70 bst(0.,0.,0.) 71 { 71 { 72 theGamma = G4Gamma::Gamma(); 72 theGamma = G4Gamma::Gamma(); 73 theElectron = G4Electron::Electron(); 73 theElectron = G4Electron::Electron(); 74 lowestSecondaryEnergy = 10*eV; 74 lowestSecondaryEnergy = 10*eV; 75 limitFactor = 4; 75 limitFactor = 4; 76 fProbabilities.resize(9,0.0); 76 fProbabilities.resize(9,0.0); 77 SetDeexcitationFlag(true); 77 SetDeexcitationFlag(true); 78 fParticleChange = nullptr; 78 fParticleChange = nullptr; 79 fAtomDeexcitation = nullptr; 79 fAtomDeexcitation = nullptr; 80 } 80 } 81 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 83 83 84 G4KleinNishinaModel::~G4KleinNishinaModel() = << 84 G4KleinNishinaModel::~G4KleinNishinaModel() >> 85 {} 85 86 86 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 87 88 88 void G4KleinNishinaModel::Initialise(const G4P 89 void G4KleinNishinaModel::Initialise(const G4ParticleDefinition* p, 89 const G4D 90 const G4DataVector& cuts) 90 { 91 { 91 fAtomDeexcitation = G4LossTableManager::Inst 92 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 92 if(IsMaster()) { InitialiseElementSelectors( 93 if(IsMaster()) { InitialiseElementSelectors(p, cuts); } 93 if(nullptr == fParticleChange) { 94 if(nullptr == fParticleChange) { 94 fParticleChange = GetParticleChangeForGamm 95 fParticleChange = GetParticleChangeForGamma(); 95 } 96 } 96 } 97 } 97 98 98 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 99 100 100 void G4KleinNishinaModel::InitialiseLocal(cons 101 void G4KleinNishinaModel::InitialiseLocal(const G4ParticleDefinition*, 101 G4VE 102 G4VEmModel* masterModel) 102 { 103 { 103 SetElementSelectors(masterModel->GetElementS 104 SetElementSelectors(masterModel->GetElementSelectors()); 104 } 105 } 105 106 106 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 108 108 G4double 109 G4double 109 G4KleinNishinaModel::ComputeCrossSectionPerAto 110 G4KleinNishinaModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 110 111 G4double gammaEnergy, 111 112 G4double Z, G4double, 112 113 G4double, G4double) 113 { 114 { 114 G4double xSection = 0.0 ; 115 G4double xSection = 0.0 ; 115 if (gammaEnergy <= LowEnergyLimit()) { retur 116 if (gammaEnergy <= LowEnergyLimit()) { return xSection; } 116 117 117 static const G4double a = 20.0 , b = 230.0 , 118 static const G4double a = 20.0 , b = 230.0 , c = 440.0; 118 119 119 static const G4double 120 static const G4double 120 d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLH 121 d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLHEP::barn, 121 d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLH 122 d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLHEP::barn, 122 e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLH 123 e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLHEP::barn, 123 e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLH 124 e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLHEP::barn, 124 f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLH 125 f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLHEP::barn, 125 f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLH 126 f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLHEP::barn; 126 127 127 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = 128 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z), 128 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = 129 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z); 129 130 130 G4double T0 = 15.0*keV; 131 G4double T0 = 15.0*keV; 131 if (Z < 1.5) { T0 = 40.0*keV; } 132 if (Z < 1.5) { T0 = 40.0*keV; } 132 133 133 G4double X = max(gammaEnergy, T0) / electr 134 G4double X = max(gammaEnergy, T0) / electron_mass_c2; 134 xSection = p1Z*G4Log(1.+2.*X)/X 135 xSection = p1Z*G4Log(1.+2.*X)/X 135 + (p2Z + p3Z*X + p4Z*X*X)/(1. + 136 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X); 136 137 137 // modification for low energy. (special ca 138 // modification for low energy. (special case for Hydrogen) 138 static const G4double dT0 = keV; 139 static const G4double dT0 = keV; 139 if (gammaEnergy < T0) { 140 if (gammaEnergy < T0) { 140 X = (T0+dT0) / electron_mass_c2 ; 141 X = (T0+dT0) / electron_mass_c2 ; 141 G4double sigma = p1Z*G4Log(1.+2*X)/X 142 G4double sigma = p1Z*G4Log(1.+2*X)/X 142 + (p2Z + p3Z*X + p4Z*X*X)/ 143 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X); 143 G4double c1 = -T0*(sigma-xSection)/(xSec 144 G4double c1 = -T0*(sigma-xSection)/(xSection*dT0); 144 G4double c2 = 0.150; 145 G4double c2 = 0.150; 145 if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); 146 if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); } 146 G4double y = G4Log(gammaEnergy/T0); 147 G4double y = G4Log(gammaEnergy/T0); 147 xSection *= G4Exp(-y*(c1+c2*y)); 148 xSection *= G4Exp(-y*(c1+c2*y)); 148 } 149 } 149 150 150 if(xSection < 0.0) { xSection = 0.0; } 151 if(xSection < 0.0) { xSection = 0.0; } 151 // G4cout << "e= " << GammaEnergy << " Z= " 152 // G4cout << "e= " << GammaEnergy << " Z= " << Z 152 // << " cross= " << xSection << G4endl; 153 // << " cross= " << xSection << G4endl; 153 return xSection; 154 return xSection; 154 } 155 } 155 156 156 //....oooOO0OOooo........oooOO0OOooo........oo 157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 158 158 void G4KleinNishinaModel::SampleSecondaries( 159 void G4KleinNishinaModel::SampleSecondaries( 159 std::vector<G4Dyn 160 std::vector<G4DynamicParticle*>* fvect, 160 const G4MaterialC 161 const G4MaterialCutsCouple* couple, 161 const G4DynamicPa 162 const G4DynamicParticle* aDynamicGamma, 162 G4double, 163 G4double, 163 G4double) 164 G4double) 164 { 165 { 165 // primary gamma 166 // primary gamma 166 G4double energy = aDynamicGamma->GetKineticE 167 G4double energy = aDynamicGamma->GetKineticEnergy(); 167 168 168 // do nothing below the threshold 169 // do nothing below the threshold 169 if(energy <= LowEnergyLimit()) { return; } 170 if(energy <= LowEnergyLimit()) { return; } 170 171 171 G4ThreeVector direction = aDynamicGamma->Get 172 G4ThreeVector direction = aDynamicGamma->GetMomentumDirection(); 172 173 173 // select atom 174 // select atom 174 const G4Element* elm = SelectRandomAtom(coup 175 const G4Element* elm = SelectRandomAtom(couple, theGamma, energy); 175 176 176 // select shell first 177 // select shell first 177 G4int nShells = elm->GetNbOfAtomicShells(); 178 G4int nShells = elm->GetNbOfAtomicShells(); 178 if(nShells > (G4int)fProbabilities.size()) { 179 if(nShells > (G4int)fProbabilities.size()) { fProbabilities.resize(nShells); } 179 G4double totprob = 0.0; 180 G4double totprob = 0.0; 180 G4int i; 181 G4int i; 181 for(i=0; i<nShells; ++i) { 182 for(i=0; i<nShells; ++i) { 182 //G4double bindingEnergy = elm->GetAtomicS 183 //G4double bindingEnergy = elm->GetAtomicShell(i); 183 totprob += elm->GetNbOfShellElectrons(i); 184 totprob += elm->GetNbOfShellElectrons(i); 184 //totprob += elm->GetNbOfShellElectrons(i) 185 //totprob += elm->GetNbOfShellElectrons(i)/(bindingEnergy*bindingEnergy); 185 fProbabilities[i] = totprob; 186 fProbabilities[i] = totprob; 186 } 187 } 187 188 188 // Loop on sampling 189 // Loop on sampling 189 static const G4int nlooplim = 1000; 190 static const G4int nlooplim = 1000; 190 G4int nloop = 0; 191 G4int nloop = 0; 191 192 192 G4double bindingEnergy, ePotEnergy, eKinEner 193 G4double bindingEnergy, ePotEnergy, eKinEnergy; 193 G4double gamEnergy0, gamEnergy1; 194 G4double gamEnergy0, gamEnergy1; 194 195 195 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 196 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine(); 196 G4double rndm[4]; 197 G4double rndm[4]; 197 198 198 do { 199 do { 199 ++nloop; 200 ++nloop; 200 201 201 // 4 random numbers to select e- 202 // 4 random numbers to select e- 202 rndmEngineMod->flatArray(4, rndm); 203 rndmEngineMod->flatArray(4, rndm); 203 G4double xprob = totprob*rndm[0]; 204 G4double xprob = totprob*rndm[0]; 204 205 205 // select shell 206 // select shell 206 for(i=0; i<nShells; ++i) { if(xprob <= fPr 207 for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) { break; } } 207 208 208 bindingEnergy = elm->GetAtomicShell(i); 209 bindingEnergy = elm->GetAtomicShell(i); 209 lv1.set(0.0,0.0,energy,energy); 210 lv1.set(0.0,0.0,energy,energy); 210 /* 211 /* 211 G4cout << "nShells= " << nShells << " i= " 212 G4cout << "nShells= " << nShells << " i= " << i 212 << " Egamma= " << energy << " Ebind= " 213 << " Egamma= " << energy << " Ebind= " << bindingEnergy 213 << G4endl; 214 << G4endl; 214 */ 215 */ 215 // for rest frame of the electron 216 // for rest frame of the electron 216 G4double x = -G4Log(rndm[1]); 217 G4double x = -G4Log(rndm[1]); 217 eKinEnergy = bindingEnergy*x; 218 eKinEnergy = bindingEnergy*x; 218 ePotEnergy = bindingEnergy*(1.0 + x); 219 ePotEnergy = bindingEnergy*(1.0 + x); 219 220 220 // for rest frame of the electron 221 // for rest frame of the electron 221 G4double eTotMomentum = sqrt(eKinEnergy*(e 222 G4double eTotMomentum = sqrt(eKinEnergy*(eKinEnergy + 2*electron_mass_c2)); 222 G4double phi = rndm[2]*twopi; 223 G4double phi = rndm[2]*twopi; 223 G4double costet = 2*rndm[3] - 1; 224 G4double costet = 2*rndm[3] - 1; 224 G4double sintet = sqrt((1 - costet)*(1 + c 225 G4double sintet = sqrt((1 - costet)*(1 + costet)); 225 lv2.set(eTotMomentum*sintet*cos(phi),eTotM 226 lv2.set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi), 226 eTotMomentum*costet,eKinEnergy + e 227 eTotMomentum*costet,eKinEnergy + electron_mass_c2); 227 bst = lv2.boostVector(); 228 bst = lv2.boostVector(); 228 lv1.boost(-bst); 229 lv1.boost(-bst); 229 230 230 gamEnergy0 = lv1.e(); 231 gamEnergy0 = lv1.e(); 231 232 232 // In the rest frame of the electron 233 // In the rest frame of the electron 233 // The scattered gamma energy is sampled a 234 // The scattered gamma energy is sampled according to Klein-Nishina formula 234 // The random number techniques of Butcher 235 // The random number techniques of Butcher & Messel are used 235 // (Nuc Phys 20(1960),15). 236 // (Nuc Phys 20(1960),15). 236 G4double E0_m = gamEnergy0/electron_mass_c 237 G4double E0_m = gamEnergy0/electron_mass_c2; 237 238 238 //G4cout << "Nloop= "<< nloop << " Ecm(keV 239 //G4cout << "Nloop= "<< nloop << " Ecm(keV)= " << gamEnergy0/keV << G4endl; 239 // 240 // 240 // sample the energy rate of the scattered 241 // sample the energy rate of the scattered gamma 241 // 242 // 242 243 243 G4double epsilon, epsilonsq, onecost, sint 244 G4double epsilon, epsilonsq, onecost, sint2, greject ; 244 245 245 G4double eps0 = 1./(1 + 2*E0_m); 246 G4double eps0 = 1./(1 + 2*E0_m); 246 G4double epsilon0sq = eps0*eps0; 247 G4double epsilon0sq = eps0*eps0; 247 G4double alpha1 = - G4Log(eps0); 248 G4double alpha1 = - G4Log(eps0); 248 G4double alpha2 = alpha1 + 0.5*(1 - ep 249 G4double alpha2 = alpha1 + 0.5*(1 - epsilon0sq); 249 250 250 do { 251 do { 251 ++nloop; 252 ++nloop; 252 // false interaction if too many iterati 253 // false interaction if too many iterations 253 if(nloop > nlooplim) { return; } 254 if(nloop > nlooplim) { return; } 254 255 255 // 3 random numbers to sample scattering 256 // 3 random numbers to sample scattering 256 rndmEngineMod->flatArray(3, rndm); 257 rndmEngineMod->flatArray(3, rndm); 257 258 258 if ( alpha1 > alpha2*rndm[0] ) { 259 if ( alpha1 > alpha2*rndm[0] ) { 259 epsilon = G4Exp(-alpha1*rndm[1]); 260 epsilon = G4Exp(-alpha1*rndm[1]); // epsilon0**r 260 epsilonsq = epsilon*epsilon; 261 epsilonsq = epsilon*epsilon; 261 262 262 } else { 263 } else { 263 epsilonsq = epsilon0sq + (1.- epsilon0 264 epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndm[1]; 264 epsilon = sqrt(epsilonsq); 265 epsilon = sqrt(epsilonsq); 265 } 266 } 266 267 267 onecost = (1.- epsilon)/(epsilon*E0_m); 268 onecost = (1.- epsilon)/(epsilon*E0_m); 268 sint2 = onecost*(2.-onecost); 269 sint2 = onecost*(2.-onecost); 269 greject = 1. - epsilon*sint2/(1.+ epsilo 270 greject = 1. - epsilon*sint2/(1.+ epsilonsq); 270 271 271 // Loop checking, 03-Aug-2015, Vladimir 272 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 272 } while (greject < rndm[2]); 273 } while (greject < rndm[2]); 273 gamEnergy1 = epsilon*gamEnergy0; 274 gamEnergy1 = epsilon*gamEnergy0; 274 275 275 // before scattering total 4-momentum in e 276 // before scattering total 4-momentum in e- system 276 lv2.set(0.0,0.0,0.0,electron_mass_c2); 277 lv2.set(0.0,0.0,0.0,electron_mass_c2); 277 lv2 += lv1; 278 lv2 += lv1; 278 279 279 // 280 // 280 // scattered gamma angles. ( Z - axis alon 281 // scattered gamma angles. ( Z - axis along the parent gamma) 281 // 282 // 282 if(sint2 < 0.0) { sint2 = 0.0; } 283 if(sint2 < 0.0) { sint2 = 0.0; } 283 costet = 1. - onecost; 284 costet = 1. - onecost; 284 sintet = sqrt(sint2); 285 sintet = sqrt(sint2); 285 phi = twopi * rndmEngineMod->flat(); 286 phi = twopi * rndmEngineMod->flat(); 286 287 287 // e- recoil 288 // e- recoil 288 // 289 // 289 // in rest frame of the electron 290 // in rest frame of the electron 290 G4ThreeVector gamDir = lv1.vect().unit(); 291 G4ThreeVector gamDir = lv1.vect().unit(); 291 G4ThreeVector v = G4ThreeVector(sintet*cos 292 G4ThreeVector v = G4ThreeVector(sintet*cos(phi),sintet*sin(phi),costet); 292 v.rotateUz(gamDir); 293 v.rotateUz(gamDir); 293 lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(), 294 lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(),gamEnergy1*v.z(),gamEnergy1); 294 lv2 -= lv1; 295 lv2 -= lv1; 295 //G4cout<<"Egam(keV)= " << lv1.e()/keV 296 //G4cout<<"Egam(keV)= " << lv1.e()/keV 296 // <<" Ee(keV)= " << (lv2.e()-ele 297 // <<" Ee(keV)= " << (lv2.e()-electron_mass_c2)/keV << G4endl; 297 lv2.boost(bst); 298 lv2.boost(bst); 298 eKinEnergy = lv2.e() - electron_mass_c2 - 299 eKinEnergy = lv2.e() - electron_mass_c2 - ePotEnergy; 299 //G4cout << "Nloop= " << nloop << " eKinEn 300 //G4cout << "Nloop= " << nloop << " eKinEnergy= " << eKinEnergy << G4endl; 300 301 301 // Loop checking, 03-Aug-2015, Vladimir Iv 302 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 302 } while ( eKinEnergy < 0.0 ); 303 } while ( eKinEnergy < 0.0 ); 303 304 304 // 305 // 305 // update G4VParticleChange for the scattere 306 // update G4VParticleChange for the scattered gamma 306 // 307 // 307 308 308 lv1.boost(bst); 309 lv1.boost(bst); 309 gamEnergy1 = lv1.e(); 310 gamEnergy1 = lv1.e(); 310 if(gamEnergy1 > lowestSecondaryEnergy) { 311 if(gamEnergy1 > lowestSecondaryEnergy) { 311 G4ThreeVector gamDirection1 = lv1.vect().u 312 G4ThreeVector gamDirection1 = lv1.vect().unit(); 312 gamDirection1.rotateUz(direction); 313 gamDirection1.rotateUz(direction); 313 fParticleChange->ProposeMomentumDirection( 314 fParticleChange->ProposeMomentumDirection(gamDirection1); 314 } else { 315 } else { 315 fParticleChange->ProposeTrackStatus(fStopA 316 fParticleChange->ProposeTrackStatus(fStopAndKill); 316 gamEnergy1 = 0.0; 317 gamEnergy1 = 0.0; 317 } 318 } 318 fParticleChange->SetProposedKineticEnergy(ga 319 fParticleChange->SetProposedKineticEnergy(gamEnergy1); 319 320 320 // 321 // 321 // kinematic of the scattered electron 322 // kinematic of the scattered electron 322 // 323 // 323 324 324 if(eKinEnergy > lowestSecondaryEnergy) { 325 if(eKinEnergy > lowestSecondaryEnergy) { 325 G4ThreeVector eDirection = lv2.vect().unit 326 G4ThreeVector eDirection = lv2.vect().unit(); 326 eDirection.rotateUz(direction); 327 eDirection.rotateUz(direction); 327 auto dp = new G4DynamicParticle(theElectro << 328 G4DynamicParticle* dp = >> 329 new G4DynamicParticle(theElectron,eDirection,eKinEnergy); 328 fvect->push_back(dp); 330 fvect->push_back(dp); 329 } else { eKinEnergy = 0.0; } 331 } else { eKinEnergy = 0.0; } 330 332 331 G4double edep = energy - gamEnergy1 - eKinEn 333 G4double edep = energy - gamEnergy1 - eKinEnergy; 332 G4double esec = 0.0; 334 G4double esec = 0.0; 333 335 334 // sample deexcitation 336 // sample deexcitation 335 // 337 // 336 if(nullptr != fAtomDeexcitation) { << 338 if(fAtomDeexcitation) { 337 G4int index = couple->GetIndex(); 339 G4int index = couple->GetIndex(); 338 if(fAtomDeexcitation->CheckDeexcitationAct 340 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) { 339 G4int Z = elm->GetZasInt(); 341 G4int Z = elm->GetZasInt(); 340 auto as = (G4AtomicShellEnumerator)(i); << 342 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i); 341 const G4AtomicShell* shell = fAtomDeexci 343 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 342 G4int nbefore = (G4int)fvect->size(); << 344 G4int nbefore = fvect->size(); 343 fAtomDeexcitation->GenerateParticles(fve 345 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index); 344 G4int nafter = (G4int)fvect->size(); << 346 G4int nafter = fvect->size(); 345 //G4cout << "N1= " << nbefore << " N2= 347 //G4cout << "N1= " << nbefore << " N2= " << nafter << G4endl; 346 for (G4int j=nbefore; j<nafter; ++j) { 348 for (G4int j=nbefore; j<nafter; ++j) { 347 G4double e = ((*fvect)[j])->GetKinetic 349 G4double e = ((*fvect)[j])->GetKineticEnergy(); 348 if(esec + e > edep) { 350 if(esec + e > edep) { 349 // correct energy in order to have e 351 // correct energy in order to have energy balance 350 e = edep - esec; 352 e = edep - esec; 351 ((*fvect)[j])->SetKineticEnergy(e); 353 ((*fvect)[j])->SetKineticEnergy(e); 352 esec += e; 354 esec += e; 353 /* 355 /* 354 G4cout << "### G4KleinNishinaModel 356 G4cout << "### G4KleinNishinaModel Edep(eV)= " << edep/eV 355 << " Esec(eV)= " << esec/eV 357 << " Esec(eV)= " << esec/eV 356 << " E["<< j << "](eV)= " < 358 << " E["<< j << "](eV)= " << e/eV 357 << " N= " << nafter 359 << " N= " << nafter 358 << " Z= " << Z << " shell= 360 << " Z= " << Z << " shell= " << i 359 << " Ebind(keV)= " << bind 361 << " Ebind(keV)= " << bindingEnergy/keV 360 << " Eshell(keV)= " << she 362 << " Eshell(keV)= " << shell->BindingEnergy()/keV 361 << G4endl; 363 << G4endl; 362 */ 364 */ 363 // delete the rest of secondaries (s 365 // delete the rest of secondaries (should not happens) 364 for (G4int jj=nafter-1; jj>j; --jj) 366 for (G4int jj=nafter-1; jj>j; --jj) { 365 delete (*fvect)[jj]; 367 delete (*fvect)[jj]; 366 fvect->pop_back(); 368 fvect->pop_back(); 367 } 369 } 368 break; 370 break; 369 } 371 } 370 esec += e; 372 esec += e; 371 } 373 } 372 edep -= esec; 374 edep -= esec; 373 } 375 } 374 } 376 } 375 if(std::abs(energy - gamEnergy1 - eKinEnergy 377 if(std::abs(energy - gamEnergy1 - eKinEnergy - esec - edep) > eV) { 376 G4cout << "### G4KleinNishinaModel dE(eV)= 378 G4cout << "### G4KleinNishinaModel dE(eV)= " 377 << (energy - gamEnergy1 - eKinEnerg 379 << (energy - gamEnergy1 - eKinEnergy - esec - edep)/eV 378 << " shell= " << i 380 << " shell= " << i 379 << " E(keV)= " << energy/keV 381 << " E(keV)= " << energy/keV 380 << " Ebind(keV)= " << bindingEnerg 382 << " Ebind(keV)= " << bindingEnergy/keV 381 << " Eg(keV)= " << gamEnergy1/keV 383 << " Eg(keV)= " << gamEnergy1/keV 382 << " Ee(keV)= " << eKinEnergy/keV 384 << " Ee(keV)= " << eKinEnergy/keV 383 << " Esec(keV)= " << esec/keV 385 << " Esec(keV)= " << esec/keV 384 << " Edep(keV)= " << edep/keV 386 << " Edep(keV)= " << edep/keV 385 << G4endl; 387 << G4endl; 386 } 388 } 387 // energy balance 389 // energy balance 388 if(edep > 0.0) { 390 if(edep > 0.0) { 389 fParticleChange->ProposeLocalEnergyDeposit 391 fParticleChange->ProposeLocalEnergyDeposit(edep); 390 } 392 } 391 } 393 } 392 394 393 //....oooOO0OOooo........oooOO0OOooo........oo 395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 394 396 395 397