source: trunk/source/processes/electromagnetic/standard/src/G4WentzelOKandVIxSection.cc @ 1348

Last change on this file since 1348 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 11.0 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26// $Id: G4WentzelOKandVIxSection.cc,v 1.12 2010/11/04 17:30:32 vnivanch Exp $
27// GEANT4 tag $Name: emstand-V09-03-25 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:   G4WentzelOKandVIxSection
35//
36// Author:      V.Ivanchenko
37//
38// Creation date: 09.04.2008 from G4MuMscModel
39//
40// Modifications:
41//
42//
43
44// -------------------------------------------------------------------
45//
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49
50#include "G4WentzelOKandVIxSection.hh"
51#include "Randomize.hh"
52#include "G4Electron.hh"
53#include "G4Positron.hh"
54#include "G4Proton.hh"
55#include "G4LossTableManager.hh"
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59G4double G4WentzelOKandVIxSection::ScreenRSquare[] = {0.0};
60G4double G4WentzelOKandVIxSection::FormFactor[]    = {0.0};
61
62using namespace std;
63
64G4WentzelOKandVIxSection::G4WentzelOKandVIxSection() :
65  numlimit(0.1),
66  nwarnings(0),
67  nwarnlimit(50),
68  alpha2(fine_structure_const*fine_structure_const)
69{
70  fNistManager = G4NistManager::Instance();
71  fG4pow = G4Pow::GetInstance();
72  theElectron = G4Electron::Electron();
73  thePositron = G4Positron::Positron();
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] = alpha2*x*x;
91      x = fNistManager->GetA27(j);
92      FormFactor[j] = constn*x*x;
93    } 
94  }
95  currentMaterial = 0;
96  elecXSRatio = factB = formfactA = screenZ = 0.0;
97  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = 1.0;
98
99  Initialise(theElectron, 1.0);
100  SetTargetMass(proton_mass_c2);
101}
102
103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104
105G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection()
106{}
107
108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109
110void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p, 
111                                          G4double CosThetaLim)
112{
113  SetupParticle(p);
114  tkin = mom2 = 0.0;
115  ecut = etag = DBL_MAX;
116  targetZ = 0;
117  cosThetaMax = CosThetaLim;
118  G4double a = 
119    G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
120  factorA2 = 0.5*a*a;
121  currentMaterial = 0;
122}
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p)
127{
128  particle = p;
129  mass = particle->GetPDGMass();
130  spin = particle->GetPDGSpin();
131  if(0.0 != spin) { spin = 0.5; }
132  G4double q = std::fabs(particle->GetPDGCharge()/eplus);
133  chargeSquare = q*q;
134  charge3 = chargeSquare*q;
135  tkin = 0.0;
136  currentMaterial = 0;
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140 
141G4double
142G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut)
143{
144  G4double cosTetMaxNuc2 = cosTetMaxNuc;
145  if(Z != targetZ || tkin != etag) {
146    etag    = tkin; 
147    targetZ = Z;
148    if(targetZ > 99) { targetZ = 99; }
149    SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*amu_c2);
150    kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2;
151
152    screenZ = ScreenRSquare[targetZ]/mom2;
153    if(Z > 1) {
154      G4double tau = tkin/mass;
155      screenZ *=std::min(Z*invbeta2,
156        (1.13 +3.76*Z*Z*invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(Z)))));
157    }
158    if(targetZ == 1 && cosTetMaxNuc < 0.0 && particle == theProton) {
159      cosTetMaxNuc = 0.0;
160    }
161    formfactA = FormFactor[targetZ]*mom2;
162
163    // allowing do not compute scattering off e-
164    cosTetMaxElec = 1.0;
165    if(cut < DBL_MAX) { 
166      if(mass < MeV) { 
167        if(cosTetMaxNuc < 1.0 && cosTetMaxNuc > 0.0 && tkin < 10*cut) { 
168          cosTetMaxNuc2 *= 0.1*tkin/cut;
169        }
170      }
171      ComputeMaxElectronScattering(cut); 
172    }
173  }
174  return cosTetMaxNuc2;
175} 
176
177//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178
179G4double
180G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
181{
182  G4double xsec = 0.0;
183  if(cosTMax >= 1.0) { return xsec; }
184 
185  G4double xSection = 0.0;
186  G4double x = 0; 
187  G4double y = 0;
188  G4double x1= 0;
189  G4double x2= 0;
190  G4double xlog = 0.0;
191
192  G4double costm = std::max(cosTMax,cosTetMaxElec); 
193  G4double fb = screenZ*factB;
194
195  // scattering off electrons
196  if(costm < 1.0) {
197    x = (1.0 - costm)/screenZ;
198    x1= x/(1 + x);
199    if(x < numlimit) { 
200      x2 = 0.5*x*x;
201      y  = x2*(1.0 - 1.3333333*x + 3*x2); 
202    } else { 
203      xlog = log(1.0 + x); 
204      y = xlog - x1; 
205    }
206
207    if(0.0 < factB) {
208      if(x < numlimit) { y -= fb*x2*x*(0.6666667 - x); } 
209      else             { y -= fb*(x + x1 - 2*xlog); }
210    }
211
212    if(y < 0.0) {
213      ++nwarnings;
214      if(nwarnings < nwarnlimit) {
215        G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
216               << G4endl;
217        G4cout << "y= " << y
218               << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
219               << " Z= " << targetZ << "  " 
220               << particle->GetParticleName() << G4endl;
221        G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
222               << " x= " << x << G4endl;
223      }
224      y = 0.0;
225    }
226    xSection = y;
227  }
228  /*
229  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
230         << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
231         << " cut(MeV)= " << ecut/MeV 
232         << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
233         << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
234         << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
235  */
236  // scattering off nucleus
237  if(cosTMax < 1.0) {
238    x = (1.0 - cosTMax)/screenZ;
239    x1= x/(1 + x);
240    if(x < numlimit) { 
241      x2 = 0.5*x*x;
242      y  = x2*(1.0 - 1.3333333*x + 3*x2); 
243    } else { 
244      xlog = log(1.0 + x); 
245      y = xlog - x1; 
246    }
247
248    if(0.0 < factB) {
249      if(x < numlimit) { y -= fb*x2*x*(0.6666667 - x); } 
250      else             { y -= fb*(x + x1 - 2*xlog); }
251    }
252    if(y < 0.0) {
253      ++nwarnings;
254      if(nwarnings < nwarnlimit) {
255        G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
256               << G4endl;
257        G4cout << "y= " << y
258               << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
259               << particle->GetParticleName() << G4endl;
260        G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
261               << " x= " << " x1= " << x1 <<G4endl;
262      }
263      y = 0.0;
264    }
265    xSection += y*targetZ; 
266  }
267  xSection *= kinFactor;
268  /*
269  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
270         << " screenZ= " << screenZ << " formF= " << formfactA
271         << " for " << particle->GetParticleName()
272         << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
273         << " x= " << x
274         << G4endl;
275  */
276  return xSection; 
277}
278
279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
280
281G4ThreeVector
282G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin,
283                                                 G4double cosTMax,
284                                                 G4double elecRatio)
285{
286  G4ThreeVector v(0.0,0.0,1.0);
287 
288  G4double formf = formfactA;
289  G4double cost1 = cosTMin;
290  G4double cost2 = cosTMax;
291  if(elecRatio > 0.0) {
292    if(G4UniformRand() <= elecRatio) {
293      formf = 0.0;
294      cost1 = std::max(cost1,cosTetMaxElec);
295      cost2 = std::max(cost2,cosTetMaxElec);
296    }
297  }
298  if(cost1 < cost2) { return v; }
299
300  G4double w1 = 1. - cost1 + screenZ;
301  G4double w2 = 1. - cost2 + screenZ;
302  G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
303  if(factB > 0.0 || formf > 0.0 || factD > 0.01) {
304    G4double grej = (1. - z1*factB)/
305      ( (1.0 + z1*factD)*(1.0 + formf*z1)*(1.0 + formf*z1) );
306    if( G4UniformRand() > grej ) { return v; }
307  } 
308  G4double cost = 1.0 - z1;
309  if(cost > 1.0)       { cost = 1.0; }
310  else if(cost < -1.0) { cost =-1.0; }
311  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
312  //G4cout << "sint= " << sint << G4endl;
313  G4double phi  = twopi*G4UniformRand();
314  G4double vx1 = sint*cos(phi);
315  G4double vy1 = sint*sin(phi);
316
317  // only direction is changed
318  v.set(vx1,vy1,cost);
319  return v;
320}
321
322//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
323
324void 
325G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy)
326{
327  G4double tmax = tkin;
328  if(mass > MeV) {
329    G4double ratio = electron_mass_c2/mass;
330    G4double tau = tkin/mass;
331    tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
332      (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); 
333    cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
334  } else {
335
336    if(particle == theElectron) { tmax *= 0.5; }
337    G4double t = std::min(cutEnergy, tmax);
338    G4double mom21 = t*(t + 2.0*electron_mass_c2);
339    G4double t1 = tkin - t;
340    //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
341    //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
342    if(t1 > 0.0) {
343      G4double mom22 = t1*(t1 + 2.0*mass);
344      G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
345      if(ctm <  1.0) { cosTetMaxElec = ctm; }
346      //if(ctm < -1.0) { cosTetMaxElec = -1.0;}
347      if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; }
348    }
349  }
350}
351
352//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.