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: G4WentzelVIRelXSection.cc 85306 2014-10-27 14:17:47Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class file 30 // GEANT4 Class file 30 // 31 // 31 // 32 // 32 // File name: G4WentzelVIRelXSection 33 // File name: G4WentzelVIRelXSection 33 // 34 // 34 // Author: V.Ivanchenko 35 // Author: V.Ivanchenko 35 // 36 // 36 // Creation date: 08.06.2012 from G4WentzelOKa 37 // Creation date: 08.06.2012 from G4WentzelOKandVIxSection 37 // 38 // 38 // Modifications: 39 // Modifications: 39 // 40 // 40 // 41 // 41 42 42 // ------------------------------------------- 43 // ------------------------------------------------------------------- 43 // 44 // 44 45 45 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 46 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 47 48 48 #include "G4WentzelVIRelXSection.hh" 49 #include "G4WentzelVIRelXSection.hh" >> 50 #include "G4PhysicalConstants.hh" >> 51 #include "G4SystemOfUnits.hh" >> 52 #include "Randomize.hh" >> 53 #include "G4Electron.hh" >> 54 #include "G4Positron.hh" >> 55 #include "G4Proton.hh" >> 56 #include "G4EmParameters.hh" >> 57 #include "G4Log.hh" 49 58 50 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 51 60 52 G4WentzelVIRelXSection::G4WentzelVIRelXSection << 61 G4double G4WentzelVIRelXSection::ScreenRSquare[] = {0.0}; >> 62 G4double G4WentzelVIRelXSection::FormFactor[] = {0.0}; >> 63 >> 64 using namespace std; >> 65 >> 66 G4WentzelVIRelXSection::G4WentzelVIRelXSection(G4bool combined) : >> 67 numlimit(0.1), >> 68 nwarnings(0), >> 69 nwarnlimit(50), >> 70 isCombined(combined), >> 71 alpha2(fine_structure_const*fine_structure_const) >> 72 { >> 73 fNistManager = G4NistManager::Instance(); >> 74 fG4pow = G4Pow::GetInstance(); >> 75 theElectron = G4Electron::Electron(); >> 76 thePositron = G4Positron::Positron(); >> 77 theProton = G4Proton::Proton(); >> 78 lowEnergyLimit = 1.0*eV; >> 79 G4double p0 = electron_mass_c2*classic_electr_radius; >> 80 coeff = twopi*p0*p0; >> 81 particle = 0; >> 82 >> 83 // Thomas-Fermi screening radii >> 84 // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282 >> 85 >> 86 if(0.0 == ScreenRSquare[0]) { >> 87 G4double a0 = electron_mass_c2/0.88534; >> 88 G4double constn = 6.937e-6/(MeV*MeV); >> 89 >> 90 ScreenRSquare[0] = alpha2*a0*a0; >> 91 for(G4int j=1; j<100; ++j) { >> 92 G4double x = a0*fG4pow->Z13(j); >> 93 //ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x; >> 94 ScreenRSquare[j] = 0.5*alpha2*x*x; >> 95 x = fNistManager->GetA27(j); >> 96 FormFactor[j] = constn*x*x; >> 97 } >> 98 } >> 99 currentMaterial = 0; >> 100 elecXSRatio = factB = factD = formfactA = screenZ = 0.0; >> 101 cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0; >> 102 cosThetaMax = 1.0; >> 103 >> 104 factB1= 0.5*CLHEP::pi*fine_structure_const; >> 105 >> 106 Initialise(theElectron, 1.0); >> 107 targetMass = proton_mass_c2; >> 108 } 53 109 54 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 55 111 56 G4WentzelVIRelXSection::~G4WentzelVIRelXSectio << 112 G4WentzelVIRelXSection::~G4WentzelVIRelXSection() >> 113 {} 57 114 58 //....oooOO0OOooo........oooOO0OOooo........oo 115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 59 116 60 G4double G4WentzelVIRelXSection::SetupKinemati << 117 void G4WentzelVIRelXSection::Initialise(const G4ParticleDefinition* p, 61 cons << 118 G4double cosThetaLim) 62 { 119 { 63 if(kinEnergy != tkin || mat != currentMateri << 120 SetupParticle(p); >> 121 tkin = mom2 = momCM2 = 0.0; >> 122 ecut = etag = DBL_MAX; >> 123 targetZ = 0; >> 124 >> 125 // cosThetaMax is below 1.0 only when MSC is combined with SS >> 126 if(isCombined) { cosThetaMax = cosThetaLim; } >> 127 >> 128 G4double a = G4EmParameters::Instance()->FactorForAngleLimit() >> 129 *CLHEP::hbarc/CLHEP::fermi; >> 130 factorA2 = 0.5*a*a; >> 131 currentMaterial = 0; >> 132 } 64 133 65 currentMaterial = mat; << 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 66 tkin = kinEnergy; << 135 67 G4double momLab2 = tkin*(tkin + 2.0*mass) << 136 void G4WentzelVIRelXSection::SetupParticle(const G4ParticleDefinition* p) 68 << 137 { 69 G4double etot = tkin + mass; << 138 particle = p; 70 G4double ptot = std::sqrt(momLab2); << 139 mass = particle->GetPDGMass(); 71 G4double m12 = mass*mass; << 140 spin = particle->GetPDGSpin(); 72 << 141 if(0.0 != spin) { spin = 0.5; } 73 // relativistic reduced mass from publucat << 142 G4double q = std::fabs(particle->GetPDGCharge()/eplus); 74 // A.P. Martynenko, R.N. Faustov, Teoret. << 143 chargeSquare = q*q; 75 << 144 charge3 = chargeSquare*q; 76 //incident particle & target nucleus << 145 tkin = 0.0; 77 G4double Ecm = std::sqrt(m12 + targetMass* << 146 currentMaterial = 0; 78 G4double mu_rel = mass*targetMass/Ecm; << 147 targetZ = 0; 79 G4double momCM = ptot*targetMass/Ecm; << 80 // relative system << 81 mom2 = momCM*momCM; << 82 invbeta2 = 1.0 + mu_rel*mu_rel/mom2; << 83 << 84 factB = spin/invbeta2; << 85 factD = std::sqrt(mom2)/targetMass; << 86 cosTetMaxNuc = isCombined ? << 87 std::max(cosThetaMax, 1.-factorA2*mat->G << 88 : cosThetaMax; << 89 } << 90 return cosTetMaxNuc; << 91 } 148 } 92 149 >> 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 151 >> 152 G4double >> 153 G4WentzelVIRelXSection::SetupTarget(G4int Z, G4double cut) >> 154 { >> 155 G4double cosTetMaxNuc2 = cosTetMaxNuc; >> 156 if(Z != targetZ || tkin != etag) { >> 157 etag = tkin; >> 158 targetZ = Z; >> 159 if(targetZ > 99) { targetZ = 99; } >> 160 SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2); >> 161 //G4double tmass2 = targetMass*targetMass; >> 162 //G4double etot = tkin + mass; >> 163 //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass; >> 164 //momCM2 = mom2*tmass2/invmass2; >> 165 //gam0pcmp = (etot + targetMass)*targetMass/invmass2; >> 166 //pcmp2 = tmass2/invmass2; >> 167 >> 168 kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2; >> 169 >> 170 screenZ = ScreenRSquare[targetZ]/mom2; >> 171 if(Z > 1) { >> 172 screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare); >> 173 /* >> 174 if(mass > MeV) { >> 175 screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare); >> 176 } else { >> 177 G4double tau = tkin/mass; >> 178 // screenZ *= std::min(Z*invbeta2, >> 179 screenZ *= std::min(Z*1.13, >> 180 (1.13 +3.76*Z*Z*invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(Z))))); >> 181 } >> 182 */ >> 183 } >> 184 if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) { >> 185 cosTetMaxNuc2 = 0.0; >> 186 } >> 187 formfactA = FormFactor[targetZ]*mom2; >> 188 >> 189 cosTetMaxElec = 1.0; >> 190 ComputeMaxElectronScattering(cut); >> 191 } >> 192 return cosTetMaxNuc2; >> 193 } >> 194 93 //....oooOO0OOooo........oooOO0OOooo........oo 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 94 196 >> 197 G4double >> 198 G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax) >> 199 { >> 200 G4double xsec = 0.0; >> 201 if(cosTMax >= 1.0) { return xsec; } >> 202 >> 203 G4double xSection = 0.0; >> 204 G4double x = 0; >> 205 G4double y = 0; >> 206 G4double x1= 0; >> 207 G4double x2= 0; >> 208 G4double xlog = 0.0; >> 209 >> 210 G4double costm = std::max(cosTMax,cosTetMaxElec); >> 211 G4double fb = screenZ*factB; >> 212 >> 213 // scattering off electrons >> 214 if(costm < 1.0) { >> 215 x = (1.0 - costm)/screenZ; >> 216 if(x < numlimit) { >> 217 x2 = 0.5*x*x; >> 218 y = x2*(1.0 - 1.3333333*x + 3*x2); >> 219 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); } >> 220 } else { >> 221 x1= x/(1 + x); >> 222 xlog = G4Log(1.0 + x); >> 223 y = xlog - x1; >> 224 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); } >> 225 } >> 226 >> 227 if(y < 0.0) { >> 228 ++nwarnings; >> 229 if(nwarnings < nwarnlimit) { >> 230 G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0" >> 231 << G4endl; >> 232 G4cout << "y= " << y >> 233 << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) >> 234 << " Z= " << targetZ << " " >> 235 << particle->GetParticleName() << G4endl; >> 236 G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ >> 237 << " x= " << x << G4endl; >> 238 } >> 239 y = 0.0; >> 240 } >> 241 xSection = y; >> 242 } >> 243 /* >> 244 G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ >> 245 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection >> 246 << " cut(MeV)= " << ecut/MeV >> 247 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ >> 248 << " zmaxN= " << (1.0 - cosThetaMax)/screenZ >> 249 << " 1-costm= " << 1.0 - cosThetaMax << G4endl; >> 250 */ >> 251 // scattering off nucleus >> 252 if(cosTMax < 1.0) { >> 253 x = (1.0 - cosTMax)/screenZ; >> 254 if(x < numlimit) { >> 255 x2 = 0.5*x*x; >> 256 y = x2*(1.0 - 1.3333333*x + 3*x2); >> 257 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); } >> 258 } else { >> 259 x1= x/(1 + x); >> 260 xlog = G4Log(1.0 + x); >> 261 y = xlog - x1; >> 262 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); } >> 263 } >> 264 >> 265 if(y < 0.0) { >> 266 ++nwarnings; >> 267 if(nwarnings < nwarnlimit) { >> 268 G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0" >> 269 << G4endl; >> 270 G4cout << "y= " << y >> 271 << " e(MeV)= " << tkin << " Z= " << targetZ << " " >> 272 << particle->GetParticleName() << G4endl; >> 273 G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ >> 274 << " x= " << " x1= " << x1 <<G4endl; >> 275 } >> 276 y = 0.0; >> 277 } >> 278 xSection += y*targetZ; >> 279 } >> 280 xSection *= kinFactor; >> 281 /* >> 282 G4cout << "Z= " << targetZ << " XStot= " << xSection/barn >> 283 << " screenZ= " << screenZ << " formF= " << formfactA >> 284 << " for " << particle->GetParticleName() >> 285 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2) >> 286 << " x= " << x >> 287 << G4endl; >> 288 */ >> 289 return xSection; >> 290 } >> 291 >> 292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 293 >> 294 G4ThreeVector >> 295 G4WentzelVIRelXSection::SampleSingleScattering(G4double cosTMin, >> 296 G4double cosTMax, >> 297 G4double elecRatio) >> 298 { >> 299 G4ThreeVector v(0.0,0.0,1.0); >> 300 >> 301 G4double formf = formfactA; >> 302 G4double cost1 = cosTMin; >> 303 G4double cost2 = cosTMax; >> 304 if(elecRatio > 0.0) { >> 305 if(G4UniformRand() <= elecRatio) { >> 306 formf = 0.0; >> 307 cost1 = std::max(cost1,cosTetMaxElec); >> 308 cost2 = std::max(cost2,cosTetMaxElec); >> 309 } >> 310 } >> 311 if(cost1 < cost2) { return v; } >> 312 >> 313 G4double w1 = 1. - cost1 + screenZ; >> 314 G4double w2 = 1. - cost2 + screenZ; >> 315 G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ; >> 316 >> 317 //G4double fm = 1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass); >> 318 G4double fm = 1.0 + formf*z1; >> 319 //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm ); >> 320 G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm ); >> 321 // "false" scattering >> 322 if( G4UniformRand() > grej ) { return v; } >> 323 // } >> 324 G4double cost = 1.0 - z1; >> 325 >> 326 if(cost > 1.0) { cost = 1.0; } >> 327 else if(cost < -1.0) { cost =-1.0; } >> 328 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); >> 329 //G4cout << "sint= " << sint << G4endl; >> 330 G4double phi = twopi*G4UniformRand(); >> 331 G4double vx1 = sint*cos(phi); >> 332 G4double vy1 = sint*sin(phi); >> 333 >> 334 // only direction is changed >> 335 v.set(vx1,vy1,cost); >> 336 return v; >> 337 } >> 338 >> 339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 340 >> 341 void >> 342 G4WentzelVIRelXSection::ComputeMaxElectronScattering(G4double cutEnergy) >> 343 { >> 344 if(mass > MeV) { >> 345 G4double ratio = electron_mass_c2/mass; >> 346 G4double tau = tkin/mass; >> 347 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/ >> 348 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); >> 349 //tmax = std::min(tmax, targetZ*targetZ*10*eV); >> 350 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2; >> 351 } else { >> 352 >> 353 G4double tmax = tkin; >> 354 if(particle == theElectron) { tmax *= 0.5; } >> 355 //tmax = std::min(tmax, targetZ*targetZ*10*eV); >> 356 G4double t = std::min(cutEnergy, tmax); >> 357 G4double mom21 = t*(t + 2.0*electron_mass_c2); >> 358 G4double t1 = tkin - t; >> 359 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= " >> 360 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl; >> 361 if(t1 > 0.0) { >> 362 G4double mom22 = t1*(t1 + 2.0*mass); >> 363 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22); >> 364 if(ctm < 1.0) { cosTetMaxElec = ctm; } >> 365 if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; } >> 366 } >> 367 } >> 368 } >> 369 >> 370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 95 371