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: G4WentzelOKandVIxSection.cc,v 1.14 2010-11-13 19:08:27 vnivanch Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ 26 // 28 // 27 // ------------------------------------------- 29 // ------------------------------------------------------------------- 28 // 30 // 29 // GEANT4 Class file 31 // GEANT4 Class file 30 // 32 // 31 // 33 // 32 // File name: G4WentzelOKandVIxSection 34 // File name: G4WentzelOKandVIxSection 33 // 35 // 34 // Author: V.Ivanchenko 36 // Author: V.Ivanchenko 35 // 37 // 36 // Creation date: 09.04.2008 from G4MuMscModel 38 // Creation date: 09.04.2008 from G4MuMscModel 37 // 39 // 38 // Modifications: 40 // Modifications: 39 // 41 // >> 42 // >> 43 40 // ------------------------------------------- 44 // ------------------------------------------------------------------- 41 // 45 // 42 46 43 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 44 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 49 46 #include "G4WentzelOKandVIxSection.hh" 50 #include "G4WentzelOKandVIxSection.hh" 47 #include "G4ScreeningMottCrossSection.hh" << 48 #include "G4PhysicalConstants.hh" << 49 #include "G4SystemOfUnits.hh" << 50 #include "Randomize.hh" 51 #include "Randomize.hh" 51 #include "G4Electron.hh" 52 #include "G4Electron.hh" 52 #include "G4Positron.hh" 53 #include "G4Positron.hh" 53 #include "G4Proton.hh" 54 #include "G4Proton.hh" 54 #include "G4EmParameters.hh" << 55 #include "G4LossTableManager.hh" 55 #include "G4Log.hh" << 56 #include "G4Exp.hh" << 57 #include "G4AutoLock.hh" << 58 56 59 //....oooOO0OOooo........oooOO0OOooo........oo 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 60 58 61 G4double G4WentzelOKandVIxSection::ScreenRSqua << 59 G4double G4WentzelOKandVIxSection::ScreenRSquare[] = {0.0}; 62 G4double G4WentzelOKandVIxSection::ScreenRSqua << 60 G4double G4WentzelOKandVIxSection::FormFactor[] = {0.0}; 63 G4double G4WentzelOKandVIxSection::FormFactor[ << 64 << 65 namespace << 66 { << 67 G4Mutex theWOKVIMutex = G4MUTEX_INITIALIZER; << 68 } << 69 << 70 const G4double alpha2 = CLHEP::fine_structure_ << 71 const G4double factB1= 0.5*CLHEP::pi*CLHEP::fi << 72 const G4double numlimit = 0.1; << 73 const G4int nwarnlimit = 50; << 74 61 75 using namespace std; 62 using namespace std; 76 63 77 G4WentzelOKandVIxSection::G4WentzelOKandVIxSec << 64 G4WentzelOKandVIxSection::G4WentzelOKandVIxSection() : 78 temp(0.,0.,0.), << 65 numlimit(0.1), 79 isCombined(comb) << 66 nwarnings(0), >> 67 nwarnlimit(50), >> 68 alpha2(fine_structure_const*fine_structure_const) 80 { 69 { 81 fNistManager = G4NistManager::Instance(); 70 fNistManager = G4NistManager::Instance(); 82 fG4pow = G4Pow::GetInstance(); 71 fG4pow = G4Pow::GetInstance(); 83 << 84 theElectron = G4Electron::Electron(); 72 theElectron = G4Electron::Electron(); 85 thePositron = G4Positron::Positron(); 73 thePositron = G4Positron::Positron(); 86 theProton = G4Proton::Proton(); 74 theProton = G4Proton::Proton(); >> 75 lowEnergyLimit = 1.0*eV; >> 76 G4double p0 = electron_mass_c2*classic_electr_radius; >> 77 coeff = twopi*p0*p0; >> 78 particle = 0; >> 79 >> 80 // Thomas-Fermi screening radii >> 81 // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282 >> 82 >> 83 if(0.0 == ScreenRSquare[0]) { >> 84 G4double a0 = electron_mass_c2/0.88534; >> 85 G4double constn = 6.937e-6/(MeV*MeV); >> 86 >> 87 ScreenRSquare[0] = alpha2*a0*a0; >> 88 for(G4int j=1; j<100; ++j) { >> 89 G4double x = a0*fG4pow->Z13(j); >> 90 ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x; >> 91 //ScreenRSquare[j] = alpha2*x*x; >> 92 x = fNistManager->GetA27(j); >> 93 FormFactor[j] = constn*x*x; >> 94 } >> 95 } >> 96 currentMaterial = 0; >> 97 elecXSRatio = factB = factD = formfactA = screenZ = 0.0; >> 98 cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0; 87 99 88 G4double p0 = CLHEP::electron_mass_c2*CLHEP: << 100 factB1= 0.5*CLHEP::pi*fine_structure_const; 89 coeff = CLHEP::twopi*p0*p0; << 101 90 targetMass = CLHEP::proton_mass_c2; << 102 Initialise(theElectron, 1.0); >> 103 targetMass = proton_mass_c2; 91 } 104 } 92 105 93 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 94 107 95 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSe 108 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection() 96 { << 109 {} 97 delete fMottXSection; << 98 } << 99 110 100 //....oooOO0OOooo........oooOO0OOooo........oo 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 101 112 102 void G4WentzelOKandVIxSection::Initialise(cons 113 void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p, 103 G4do << 114 G4double CosThetaLim) 104 { 115 { 105 SetupParticle(p); 116 SetupParticle(p); 106 tkin = mom2 = momCM2 = 0.0; 117 tkin = mom2 = momCM2 = 0.0; 107 ecut = etag = DBL_MAX; 118 ecut = etag = DBL_MAX; 108 targetZ = 0; 119 targetZ = 0; 109 << 120 cosThetaMax = CosThetaLim; 110 // cosThetaMax is below 1.0 only when MSC is << 121 G4double a = 111 if(isCombined) { cosThetaMax = cosThetaLim; << 122 G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi; 112 G4EmParameters* param = G4EmParameters::Inst << 113 G4double a = param->FactorForAngleLimit()*CL << 114 factorA2 = 0.5*a*a; 123 factorA2 = 0.5*a*a; 115 currentMaterial = nullptr; << 124 currentMaterial = 0; 116 << 117 fNucFormfactor = param->NuclearFormfactorTyp << 118 if(0.0 == ScreenRSquare[0]) { InitialiseA(); << 119 << 120 // Mott corrections always added << 121 if((p == theElectron || p == thePositron) && << 122 fMottXSection = new G4ScreeningMottCrossSe << 123 fMottXSection->Initialise(p, 1.0); << 124 } << 125 /* << 126 G4cout << "G4WentzelOKandVIxSection::Initial << 127 << p->GetParticleName() << " cosThetaMax= " << 128 << " " << ScreenRSquare[0] << " coeff= " < << 129 */ << 130 } << 131 << 132 //....oooOO0OOooo........oooOO0OOooo........oo << 133 << 134 void G4WentzelOKandVIxSection::InitialiseA() << 135 { << 136 // Thomas-Fermi screening radii << 137 // Formfactors from A.V. Butkevich et al., N << 138 if(0.0 != ScreenRSquare[0]) { return; } << 139 G4AutoLock l(&theWOKVIMutex); << 140 if(0.0 == ScreenRSquare[0]) { << 141 const G4double invmev2 = 1./(CLHEP::MeV*CL << 142 G4double a0 = CLHEP::electron_mass_c2/0.88 << 143 G4double constn = 6.937e-6*invmev2; << 144 G4double fct = G4EmParameters::Instance()- << 145 << 146 G4double afact = 0.5*fct*alpha2*a0*a0; << 147 ScreenRSquare[0] = afact; << 148 ScreenRSquare[1] = afact; << 149 ScreenRSquareElec[1] = afact; << 150 FormFactor[1] = 3.097e-6*invmev2; << 151 << 152 for(G4int j=2; j<100; ++j) { << 153 G4double x = fG4pow->Z13(j); << 154 ScreenRSquare[j] = afact*(1 + G4Exp(-j*j << 155 ScreenRSquareElec[j] = afact*x*x; << 156 x = fNistManager->GetA27(j); << 157 FormFactor[j] = constn*x*x; << 158 } << 159 } << 160 l.unlock(); << 161 } 125 } 162 126 163 //....oooOO0OOooo........oooOO0OOooo........oo 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 164 128 165 void G4WentzelOKandVIxSection::SetupParticle(c 129 void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p) 166 { 130 { 167 particle = p; 131 particle = p; 168 mass = particle->GetPDGMass(); 132 mass = particle->GetPDGMass(); 169 spin = particle->GetPDGSpin(); 133 spin = particle->GetPDGSpin(); 170 if(0.0 != spin) { spin = 0.5; } 134 if(0.0 != spin) { spin = 0.5; } 171 G4double q = std::abs(particle->GetPDGCharge << 135 G4double q = std::fabs(particle->GetPDGCharge()/eplus); 172 chargeSquare = q*q; 136 chargeSquare = q*q; 173 charge3 = chargeSquare*q; 137 charge3 = chargeSquare*q; 174 tkin = 0.0; 138 tkin = 0.0; 175 currentMaterial = nullptr; << 139 currentMaterial = 0; 176 targetZ = 0; 140 targetZ = 0; 177 } 141 } 178 142 179 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 180 << 181 G4double << 182 G4WentzelOKandVIxSection::SetupKinematic(G4dou << 183 { << 184 if(ekin != tkin || mat != currentMaterial) { << 185 currentMaterial = mat; << 186 tkin = ekin; << 187 mom2 = tkin*(tkin + 2.0*mass); << 188 invbeta2 = 1.0 + mass*mass/mom2; << 189 factB = spin/invbeta2; << 190 cosTetMaxNuc = isCombined ? << 191 std::max(cosThetaMax, 1.-factorA2*mat->G << 192 : cosThetaMax; << 193 } << 194 return cosTetMaxNuc; << 195 } << 196 << 197 //....oooOO0OOooo........oooOO0OOooo........oo << 198 144 199 G4double 145 G4double 200 G4WentzelOKandVIxSection::SetupTarget(G4int Z, 146 G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut) 201 { 147 { 202 G4double cosTetMaxNuc2 = cosTetMaxNuc; 148 G4double cosTetMaxNuc2 = cosTetMaxNuc; 203 if(Z != targetZ || tkin != etag) { 149 if(Z != targetZ || tkin != etag) { 204 etag = tkin; 150 etag = tkin; 205 targetZ = std::min(Z, 99); << 151 targetZ = Z; 206 G4double massT = (1 == Z) ? CLHEP::proton_ << 152 if(targetZ > 99) { targetZ = 99; } 207 fNistManager->GetAtomicMassAmu(Z)*CLHEP: << 153 SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2); 208 SetTargetMass(massT); << 154 //G4double tmass2 = targetMass*targetMass; 209 << 155 //G4double etot = tkin + mass; 210 kinFactor = coeff*Z*chargeSquare*invbeta2/ << 156 //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass; 211 if(particle == theElectron && fMottXSectio << 157 //momCM2 = mom2*tmass2/invmass2; 212 fMottFactor = (1.0 + 2.0e-4*Z*Z); << 158 //gam0pcmp = (etot + targetMass)*targetMass/invmass2; 213 } << 159 //pcmp2 = tmass2/invmass2; 214 160 215 if(1 == Z) { << 161 kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2; 216 screenZ = ScreenRSquare[targetZ]/mom2; << 162 217 } else if(mass > MeV) { << 163 screenZ = ScreenRSquare[targetZ]/mom2; 218 screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z << 164 if(Z > 1) { 219 ScreenRSquare[targetZ]/mom2; << 220 } else { << 221 G4double tau = tkin/mass; 165 G4double tau = tkin/mass; 222 screenZ = std::min(Z*1.13,(1.13 +3.76*Z* << 166 screenZ *=std::min(Z*invbeta2, 223 *invbeta2*alpha2*std::sqrt(tau/(tau << 167 (1.13 +3.76*Z*Z*invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(Z))))); 224 ScreenRSquareElec[targetZ]/mom2; << 225 } 168 } 226 if(targetZ == 1 && particle == theProton & << 169 if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) { 227 cosTetMaxNuc2 = 0.0; 170 cosTetMaxNuc2 = 0.0; 228 } 171 } 229 formfactA = FormFactor[targetZ]*mom2; 172 formfactA = FormFactor[targetZ]*mom2; 230 173 >> 174 // allowing do not compute scattering off e- 231 cosTetMaxElec = 1.0; 175 cosTetMaxElec = 1.0; 232 ComputeMaxElectronScattering(cut); << 176 if(cut < DBL_MAX) { >> 177 /* >> 178 if(mass < MeV) { >> 179 if(cosTetMaxNuc < 1.0 && cosTetMaxNuc > 0.0 && tkin < 10*cut) { >> 180 cosTetMaxNuc2 *= 0.1*tkin/cut; >> 181 } >> 182 } >> 183 */ >> 184 ComputeMaxElectronScattering(cut); >> 185 } 233 } 186 } 234 //G4cout << "SetupTarget: Z= " << targetZ < << 235 // << " fMottFactor= " << fMottFactor << " << 236 return cosTetMaxNuc2; 187 return cosTetMaxNuc2; 237 } 188 } 238 189 239 //....oooOO0OOooo........oooOO0OOooo........oo 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 240 191 241 G4double 192 G4double 242 G4WentzelOKandVIxSection::ComputeTransportCros 193 G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax) 243 { 194 { >> 195 G4double xsec = 0.0; >> 196 if(cosTMax >= 1.0) { return xsec; } >> 197 244 G4double xSection = 0.0; 198 G4double xSection = 0.0; 245 if(cosTMax >= 1.0) { return xSection; } << 199 G4double x = 0; >> 200 G4double y = 0; >> 201 G4double x1= 0; >> 202 G4double x2= 0; >> 203 G4double xlog = 0.0; 246 204 247 G4double costm = std::max(cosTMax,cosTetMaxE 205 G4double costm = std::max(cosTMax,cosTetMaxElec); 248 G4double fb = screenZ*factB; 206 G4double fb = screenZ*factB; 249 207 250 // scattering off electrons 208 // scattering off electrons 251 if(costm < 1.0) { 209 if(costm < 1.0) { 252 G4double x = (1.0 - costm)/screenZ; << 210 x = (1.0 - costm)/screenZ; 253 if(x < numlimit) { 211 if(x < numlimit) { 254 G4double x2 = 0.5*x*x; << 212 x2 = 0.5*x*x; 255 xSection = x2*((1.0 - 1.3333333*x + 3*x2 << 213 y = x2*(1.0 - 1.3333333*x + 3*x2); >> 214 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); } 256 } else { 215 } else { 257 G4double x1= x/(1 + x); << 216 x1= x/(1 + x); 258 G4double xlog = G4Log(1.0 + x); << 217 xlog = log(1.0 + x); 259 xSection = xlog - x1 - fb*(x + x1 - 2*xl << 218 y = xlog - x1; >> 219 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); } 260 } 220 } 261 221 262 if(xSection < 0.0) { << 222 if(y < 0.0) { 263 ++nwarnings; 223 ++nwarnings; 264 if(nwarnings < nwarnlimit) { 224 if(nwarnings < nwarnlimit) { 265 G4cout << "G4WentzelOKandVIxSection::C << 225 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0" 266 << " scattering on e- <0" << 226 << G4endl; 267 << G4endl; << 227 G4cout << "y= " << y 268 G4cout << "cross= " << xSection << 228 << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 269 << " e(MeV)= " << tkin << " p(M << 229 << " Z= " << targetZ << " " 270 << " Z= " << targetZ << " " << 230 << particle->GetParticleName() << G4endl; 271 << particle->GetParticleName() << 231 G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 272 G4cout << " 1-costm= " << 1.0-costm << << 232 << " x= " << x << G4endl; 273 << " x= " << x << G4endl; << 274 } 233 } 275 xSection = 0.0; << 234 y = 0.0; 276 } 235 } >> 236 xSection = y; 277 } 237 } 278 /* << 238 /* 279 G4cout << "G4WentzelOKandVIxSection::Com << 239 G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 280 << " Z= " << targetZ << 240 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection 281 << " e(MeV)= " << tkin/MeV << " XSel= " << 241 << " cut(MeV)= " << ecut/MeV 282 << " zmaxE= " << (1.0 - cosTetMaxElec)/s << 242 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 283 << " zmaxN= " << (1.0 - cosThetaMax)/scr << 243 << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 284 << " 1-costm= " << 1.0 - cosThetaMax << << 244 << " 1-costm= " << 1.0 - cosThetaMax << G4endl; 285 */ 245 */ 286 // scattering off nucleus 246 // scattering off nucleus 287 if(cosTMax < 1.0) { 247 if(cosTMax < 1.0) { 288 G4double x = (1.0 - cosTMax)/screenZ; << 248 x = (1.0 - cosTMax)/screenZ; 289 G4double y; << 290 if(x < numlimit) { 249 if(x < numlimit) { 291 G4double x2 = 0.5*x*x; << 250 x2 = 0.5*x*x; 292 y = x2*((1.0 - 1.3333333*x + 3*x2) - fb* << 251 y = x2*(1.0 - 1.3333333*x + 3*x2); >> 252 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); } 293 } else { 253 } else { 294 G4double x1= x/(1 + x); << 254 x1= x/(1 + x); 295 G4double xlog = G4Log(1.0 + x); << 255 xlog = log(1.0 + x); 296 y = xlog - x1 - fb*(x + x1 - 2*xlog); << 256 y = xlog - x1; >> 257 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); } 297 } 258 } 298 259 299 if(y < 0.0) { 260 if(y < 0.0) { 300 ++nwarnings; 261 ++nwarnings; 301 if(nwarnings < nwarnlimit) { 262 if(nwarnings < nwarnlimit) { 302 G4cout << "G4WentzelOKandVIxSection::C << 263 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0" 303 << " scattering on nucleus <0" << 264 << G4endl; 304 << G4endl; << 265 G4cout << "y= " << y 305 G4cout << "y= " << y << 266 << " e(MeV)= " << tkin << " Z= " << targetZ << " " 306 << " e(MeV)= " << tkin << " Z= << 267 << particle->GetParticleName() << G4endl; 307 << particle->GetParticleName() << 268 G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ 308 G4cout << " formfactA= " << formfactA << 269 << " x= " << " x1= " << x1 <<G4endl; 309 << " x= " << x <<G4endl; << 310 } 270 } 311 y = 0.0; 271 y = 0.0; 312 } 272 } 313 xSection += y*targetZ; 273 xSection += y*targetZ; 314 } 274 } 315 xSection *= kinFactor; 275 xSection *= kinFactor; 316 << 276 /* 317 /* << 318 G4cout << "Z= " << targetZ << " XStot= " << 277 G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 319 << " screenZ= " << screenZ << " formF << 278 << " screenZ= " << screenZ << " formF= " << formfactA 320 << " for " << particle->GetParticleNa << 279 << " for " << particle->GetParticleName() 321 << " m= " << mass << " 1/v= " << sqrt(invbe << 280 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2) 322 << " p= " << sqrt(mom2) << 281 << " x= " << x 323 << " x= " << x << G4endl; << 282 << G4endl; 324 */ 283 */ 325 return xSection; 284 return xSection; 326 } 285 } 327 286 328 //....oooOO0OOooo........oooOO0OOooo........oo 287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 329 288 330 G4ThreeVector& << 289 G4ThreeVector 331 G4WentzelOKandVIxSection::SampleSingleScatteri 290 G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin, 332 << 291 G4double cosTMax, 333 << 292 G4double elecRatio) 334 { 293 { 335 temp.set(0.0,0.0,1.0); << 294 G4ThreeVector v(0.0,0.0,1.0); 336 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 337 295 338 G4double formf = formfactA; 296 G4double formf = formfactA; 339 G4double cost1 = cosTMin; 297 G4double cost1 = cosTMin; 340 G4double cost2 = cosTMax; 298 G4double cost2 = cosTMax; 341 if(elecRatio > 0.0) { 299 if(elecRatio > 0.0) { 342 if(rndmEngineMod->flat() <= elecRatio) { << 300 if(G4UniformRand() <= elecRatio) { 343 formf = 0.0; 301 formf = 0.0; 344 cost1 = std::max(cost1,cosTetMaxElec); 302 cost1 = std::max(cost1,cosTetMaxElec); 345 cost2 = std::max(cost2,cosTetMaxElec); 303 cost2 = std::max(cost2,cosTetMaxElec); 346 } 304 } 347 } 305 } 348 if(cost1 > cost2) { << 306 if(cost1 < cost2) { return v; } 349 G4double w1 = 1. - cost1; << 307 350 G4double w2 = 1. - cost2; << 308 G4double w1 = 1. - cost1 + screenZ; 351 G4double w3 = rndmEngineMod->flat()*(w2 - << 309 G4double w2 = 1. - cost2 + screenZ; 352 G4double z1 = ((w2 - w3)*screenZ + w1*w2)/ << 310 G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ; 353 G4double fm = 1.0; << 311 354 << 312 //G4double fm = 1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass); 355 if(fNucFormfactor == fExponentialNF) { << 313 G4double fm = 1.0 + formf*z1; 356 fm += formf*z1; << 314 //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm ); 357 fm = 1.0/(fm*fm); << 315 G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm ); 358 } else if(fNucFormfactor == fGaussianNF) { << 316 // "false" scattering 359 fm = G4Exp(-2*formf*z1); << 317 if( G4UniformRand() > grej ) { return v; } 360 } else if(fNucFormfactor == fFlatNF) { << 318 // } 361 static const G4double ccoef = 0.00508/CL << 319 G4double cost = 1.0 - z1; 362 G4double x = std::sqrt(2.*mom2*z1)*ccoef << 320 363 fm = FlatFormfactor(x); << 321 if(cost > 1.0) { cost = 1.0; } 364 fm *= FlatFormfactor(x*0.6*fG4pow->A13(f << 322 else if(cost < -1.0) { cost =-1.0; } 365 } << 323 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); 366 // G4cout << " fm=" << fm << " " << fMott << 324 //G4cout << "sint= " << sint << G4endl; 367 G4double grej; << 325 G4double phi = twopi*G4UniformRand(); 368 if(nullptr != fMottXSection) { << 326 G4double vx1 = sint*cos(phi); 369 fMottXSection->SetupKinematic(tkin, targ << 327 G4double vy1 = sint*sin(phi); 370 grej = fMottXSection->RatioMottRutherfor << 328 371 } else { << 329 // only direction is changed 372 grej = (1. - z1*factB + factB1*targetZ*s << 330 v.set(vx1,vy1,cost); 373 *fm/(1.0 + z1*factD); << 331 return v; 374 } << 375 if(fMottFactor*rndmEngineMod->flat() <= gr << 376 // exclude "false" scattering due to for << 377 G4double cost = 1.0 - z1; << 378 if(cost > 1.0) { cost = 1.0; } << 379 else if(cost < -1.0) { cost =-1.0; } << 380 G4double sint = sqrt((1.0 - cost)*(1.0 + << 381 //G4cout << "sint= " << sint << G4endl; << 382 G4double phi = twopi*rndmEngineMod->fla << 383 temp.set(sint*cos(phi),sint*sin(phi),cos << 384 } << 385 } << 386 return temp; << 387 } 332 } 388 333 389 //....oooOO0OOooo........oooOO0OOooo........oo 334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 390 335 391 void 336 void 392 G4WentzelOKandVIxSection::ComputeMaxElectronSc 337 G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy) 393 { 338 { >> 339 G4double tmax = tkin; 394 if(mass > MeV) { 340 if(mass > MeV) { 395 G4double ratio = electron_mass_c2/mass; 341 G4double ratio = electron_mass_c2/mass; 396 G4double tau = tkin/mass; 342 G4double tau = tkin/mass; 397 G4double tmax = 2.0*electron_mass_c2*tau*( << 343 tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/ 398 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*rat << 344 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); 399 cosTetMaxElec = 1.0 - std::min(cutEnergy, 345 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2; 400 } else { 346 } else { 401 347 402 G4double tmax = (particle == theElectron) << 348 if(particle == theElectron) { tmax *= 0.5; } 403 G4double t = std::min(cutEnergy, tmax); 349 G4double t = std::min(cutEnergy, tmax); 404 G4double mom21 = t*(t + 2.0*electron_mass_ 350 G4double mom21 = t*(t + 2.0*electron_mass_c2); 405 G4double t1 = tkin - t; 351 G4double t1 = tkin - t; 406 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax 352 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= " 407 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4end 353 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl; 408 if(t1 > 0.0) { 354 if(t1 > 0.0) { 409 G4double mom22 = t1*(t1 + 2.0*mass); 355 G4double mom22 = t1*(t1 + 2.0*mass); 410 G4double ctm = (mom2 + mom22 - mom21)*0. 356 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22); 411 if(ctm < 1.0) { cosTetMaxElec = ctm; } 357 if(ctm < 1.0) { cosTetMaxElec = ctm; } 412 if(particle == theElectron && cosTetMaxE << 358 //if(ctm < -1.0) { cosTetMaxElec = -1.0;} 413 cosTetMaxElec = 0.0; << 359 if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; } 414 } << 415 } 360 } 416 } 361 } 417 } << 418 << 419 //....oooOO0OOooo........oooOO0OOooo........oo << 420 << 421 G4double << 422 G4WentzelOKandVIxSection::ComputeSecondTranspo << 423 { << 424 return 0.0; << 425 } 362 } 426 363 427 //....oooOO0OOooo........oooOO0OOooo........oo 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 428 365