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