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

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 10.8 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.10 2010/06/01 13:34:21 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
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}
96
97//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98
99G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection()
100{}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103
104void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p, 
105                                          G4double CosThetaLim)
106{
107  SetupParticle(p);
108  tkin = mom2 = 0.0;
109  ecut = etag = DBL_MAX;
110  targetZ = 0;
111  cosThetaMax = CosThetaLim;
112  G4double a = 
113    G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
114  factorA2 = 0.5*a*a;
115  currentMaterial = 0;
116}
117
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119
120void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p)
121{
122  particle = p;
123  mass = particle->GetPDGMass();
124  spin = particle->GetPDGSpin();
125  if(0.0 != spin) { spin = 0.5; }
126  G4double q = std::fabs(particle->GetPDGCharge()/eplus);
127  chargeSquare = q*q;
128  charge3 = chargeSquare*q;
129  tkin = 0.0;
130  currentMaterial = 0;
131}
132
133//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134 
135G4double
136G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut)
137{
138  G4double cosTetMaxNuc2 = cosTetMaxNuc;
139  if(Z != targetZ || tkin != etag) {
140    etag    = tkin; 
141    targetZ = Z;
142    if(targetZ > 99) { targetZ = 99; }
143    SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*amu_c2);
144    kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2;
145
146    screenZ = ScreenRSquare[targetZ]/mom2;
147    if(Z > 1) {
148      G4double tau = tkin/mass;
149      screenZ *=std::min(Z*invbeta2,
150        (1.13 +3.76*Z*Z*invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(Z)))));
151    }
152    if(targetZ == 1 && cosTetMaxNuc < 0.0 && particle == theProton) {
153      cosTetMaxNuc = 0.0;
154    }
155    formfactA = FormFactor[targetZ]*mom2;
156
157    // allowing do not compute scattering off e-
158    cosTetMaxElec = 1.0;
159    if(cut < DBL_MAX) { 
160      if(mass < MeV) { 
161        if(cosTetMaxNuc < 1.0 && cosTetMaxNuc > 0.0 && tkin < 10*cut) { 
162          cosTetMaxNuc2 *= 0.1*tkin/cut;
163        }
164      }
165      ComputeMaxElectronScattering(cut); 
166    }
167  }
168  return cosTetMaxNuc2;
169} 
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172
173G4double
174G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
175{
176  G4double xsec = 0.0;
177  if(cosTMax >= 1.0) { return xsec; }
178 
179  G4double xSection = 0.0;
180  G4double x = 0; 
181  G4double y = 0;
182  G4double x1= 0;
183  G4double x2= 0;
184  G4double xlog = 0.0;
185
186  G4double costm = std::max(cosTMax,cosTetMaxElec); 
187  G4double fb = screenZ*factB;
188
189  // scattering off electrons
190  if(costm < 1.0) {
191    x = (1.0 - costm)/screenZ;
192    x1= x/(1 + x);
193    if(x < numlimit) { 
194      x2 = 0.5*x*x;
195      y  = x2*(1.0 - 1.3333333*x + 3*x2); 
196    } else { 
197      xlog = log(1.0 + x); 
198      y = xlog - x1; 
199    }
200
201    if(0.0 < factB) {
202      if(x < numlimit) { y -= fb*x2*x*(0.6666667 - x); } 
203      else             { y -= fb*(x + x1 - 2*xlog); }
204    }
205
206    if(y < 0.0) {
207      ++nwarnings;
208      if(nwarnings < nwarnlimit) {
209        G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
210               << G4endl;
211        G4cout << "y= " << y
212               << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
213               << " Z= " << targetZ << "  " 
214               << particle->GetParticleName() << G4endl;
215        G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
216               << " x= " << x << G4endl;
217      }
218      y = 0.0;
219    }
220    xSection = y;
221  }
222  /*
223  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
224         << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
225         << " cut(MeV)= " << ecut/MeV 
226         << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
227         << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
228         << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
229  */
230  // scattering off nucleus
231  if(cosTMax < 1.0) {
232    x = (1.0 - cosTMax)/screenZ;
233    x1= x/(1 + x);
234    if(x < numlimit) { 
235      x2 = 0.5*x*x;
236      y  = x2*(1.0 - 1.3333333*x + 3*x2); 
237    } else { 
238      xlog = log(1.0 + x); 
239      y = xlog - x1; 
240    }
241
242    if(0.0 < factB) {
243      if(x < numlimit) { y -= fb*x2*x*(0.6666667 - x); } 
244      else             { y -= fb*(x + x1 - 2*xlog); }
245    }
246    if(y < 0.0) {
247      ++nwarnings;
248      if(nwarnings < nwarnlimit) {
249        G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
250               << G4endl;
251        G4cout << "y= " << y
252               << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
253               << particle->GetParticleName() << G4endl;
254        G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
255               << " x= " << " x1= " << x1 <<G4endl;
256      }
257      y = 0.0;
258    }
259    xSection += y*targetZ; 
260  }
261  xSection *= kinFactor;
262  /*
263  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
264         << " screenZ= " << screenZ << " formF= " << formfactA
265         << " for " << particle->GetParticleName()
266         << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
267         << " x= " << x
268         << G4endl;
269  */
270  return xSection; 
271}
272
273//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274
275G4ThreeVector
276G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin,
277                                                 G4double cosTMax,
278                                                 G4double elecRatio)
279{
280  G4ThreeVector v(0.0,0.0,1.0);
281 
282  G4double formf = formfactA;
283  G4double cost1 = cosTMin;
284  G4double cost2 = cosTMax;
285  if(elecRatio > 0.0) {
286    if(G4UniformRand() <= elecRatio) {
287      formf = 0.0;
288      cost1 = std::max(cost1,cosTetMaxElec);
289      cost2 = std::max(cost2,cosTetMaxElec);
290    }
291  }
292  if(cost1 < cost2) { return v; }
293
294  G4double w1 = 1. - cost1 + screenZ;
295  G4double w2 = 1. - cost2 + screenZ;
296  G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
297  if(factB > 0.0 || formf > 0.0 || factD > 0.01) {
298    G4double grej = (1. - z1*factB)/
299      ( (1.0 + z1*factD)*(1.0 + formf*z1)*(1.0 + formf*z1) );
300    if( G4UniformRand() > grej ) { return v; }
301  } 
302  G4double cost = 1.0 - z1;
303  if(cost > 1.0)       { cost = 1.0; }
304  else if(cost < -1.0) { cost =-1.0; }
305  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
306  //G4cout << "sint= " << sint << G4endl;
307  G4double phi  = twopi*G4UniformRand();
308  G4double vx1 = sint*cos(phi);
309  G4double vy1 = sint*sin(phi);
310
311  // only direction is changed
312  v.set(vx1,vy1,cost);
313  return v;
314}
315
316//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
317
318void 
319G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy)
320{
321  G4double tmax = tkin;
322  if(mass > MeV) {
323    G4double ratio = electron_mass_c2/mass;
324    G4double tau = tkin/mass;
325    tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
326      (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); 
327    cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
328  } else {
329
330    if(particle == theElectron) { tmax *= 0.5; }
331    G4double t = std::min(cutEnergy, tmax);
332    G4double mom21 = t*(t + 2.0*electron_mass_c2);
333    G4double t1 = tkin - t;
334    //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
335    //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
336    if(t1 > 0.0) {
337      G4double mom22 = t1*(t1 + 2.0*mass);
338      G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
339      if(ctm <  1.0) { cosTetMaxElec = ctm; }
340      //if(ctm < -1.0) { cosTetMaxElec = -1.0;}
341      if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; }
342    }
343  }
344}
345
346//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.