source: trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc @ 1250

Last change on this file since 1250 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
File size: 24.4 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: G4WentzelVIModel.cc,v 1.37 2009/10/28 10:14:13 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:   G4WentzelVIModel
35//
36// Author:      V.Ivanchenko
37//
38// Creation date: 09.04.2008 from G4MuMscModel
39//
40// Modifications:
41//
42//
43// Class Description:
44//
45// Implementation of the model of multiple scattering based on
46// G.Wentzel, Z. Phys. 40 (1927) 590.
47// H.W.Lewis, Phys Rev 78 (1950) 526.
48// J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
49// L.Urban, CERN-OPEN-2006-077.
50
51// -------------------------------------------------------------------
52//
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57#include "G4WentzelVIModel.hh"
58#include "Randomize.hh"
59#include "G4LossTableManager.hh"
60#include "G4ParticleChangeForMSC.hh"
61#include "G4PhysicsTableHelper.hh"
62#include "G4ElementVector.hh"
63#include "G4ProductionCutsTable.hh"
64#include "G4PhysicsLogVector.hh"
65#include "G4Electron.hh"
66#include "G4Positron.hh"
67#include "G4Proton.hh"
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70
71G4double G4WentzelVIModel::ScreenRSquare[] = {0.0};
72G4double G4WentzelVIModel::FormFactor[]    = {0.0};
73
74using namespace std;
75
76G4WentzelVIModel::G4WentzelVIModel(const G4String& nam) :
77  G4VMscModel(nam),
78  theLambdaTable(0),
79  theLambda2Table(0),
80  numlimit(0.2),
81  nbins(60),
82  nwarnings(0),
83  nwarnlimit(50),
84  currentCouple(0),
85  cosThetaMin(1.0),
86  q2Limit(TeV*TeV),
87  alpha2(fine_structure_const*fine_structure_const),
88  isInitialized(false),
89  inside(false)
90{
91  invsqrt12 = 1./sqrt(12.);
92  tlimitminfix = 1.e-6*mm;
93  theManager = G4LossTableManager::Instance(); 
94  fNistManager = G4NistManager::Instance();
95  theElectron = G4Electron::Electron();
96  thePositron = G4Positron::Positron();
97  theProton   = G4Proton::Proton();
98  lowEnergyLimit = 0.1*keV;
99  G4double p0 = electron_mass_c2*classic_electr_radius;
100  coeff  = twopi*p0*p0;
101  tkin = targetZ = mom2 = DBL_MIN;
102  ecut = etag = DBL_MAX;
103  particle = 0;
104  nelments = 5;
105  xsecn.resize(nelments);
106  prob.resize(nelments);
107
108  // Thomas-Fermi screening radii
109  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
110
111  if(0.0 == ScreenRSquare[0]) {
112    G4double a0 = electron_mass_c2/0.88534; 
113    G4double constn = 6.937e-6/(MeV*MeV);
114
115    ScreenRSquare[0] = alpha2*a0*a0;
116    for(G4int j=1; j<100; j++) {
117      G4double x = a0*fNistManager->GetZ13(j);
118      ScreenRSquare[j] = alpha2*x*x;
119      x = fNistManager->GetA27(j); 
120      FormFactor[j] = constn*x*x;
121    } 
122  }
123}
124
125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126
127G4WentzelVIModel::~G4WentzelVIModel()
128{}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131
132void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
133                                  const G4DataVector& cuts)
134{
135  // reset parameters
136  SetupParticle(p);
137  tkin = targetZ = mom2 = 0.0;
138  ecut = etag = DBL_MAX;
139  currentRange = 0.0;
140  cosThetaMax = cos(PolarAngleLimit());
141  currentCuts = &cuts;
142
143  // set values of some data members
144  if(!isInitialized) {
145    isInitialized = true;
146    fParticleChange = GetParticleChangeForMSC();
147    InitialiseSafetyHelper();
148  }
149}
150
151//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152
153G4double G4WentzelVIModel::ComputeCrossSectionPerAtom( 
154                             const G4ParticleDefinition* p,
155                             G4double kinEnergy,
156                             G4double Z, G4double,
157                             G4double cutEnergy, G4double)
158{
159  SetupParticle(p);
160  if(kinEnergy < lowEnergyLimit) return 0.0;
161  SetupKinematic(kinEnergy, cutEnergy);
162  SetupTarget(Z, kinEnergy);
163  G4double xsec = ComputeTransportXSectionPerAtom();
164  /*   
165  G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2
166         << " cosN= " << cosTetMaxNuc2 << " xsec(bn)= " << xsec/barn
167         << " " << particle->GetParticleName() << G4endl;
168  */
169  return xsec;
170}
171
172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173
174G4double G4WentzelVIModel::ComputeTransportXSectionPerAtom()
175{
176  G4double xSection = 0.0;
177  G4double x, y, x1, x2, x3, x4;
178
179  // scattering off electrons
180  if(cosTetMaxElec2 < 1.0) {
181    x = (1.0 - cosTetMaxElec2)/screenZ;
182    if(x < numlimit) y = 0.5*x*x*(1.0 - 1.3333333*x + 1.5*x*x); 
183    else             y = log(1.0 + x) - x/(1.0 + x);
184    if(y < 0.0) {
185      nwarnings++;
186      if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) {
187        G4cout << "Electron scattering <0 for L1 " << y
188               << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
189               << " Z= " << targetZ << "  " 
190               << particle->GetParticleName() << G4endl;
191        G4cout << " z= " << 1.0-cosTetMaxElec2 << " screenZ= " << screenZ
192               << " x= " << x << G4endl;
193      }
194      y = 0.0;
195    }
196    xSection = y;
197  }
198  /* 
199  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
200         << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
201         << " cut(MeV)= " << ecut/MeV 
202         << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
203         << " zmaxN= " << (1.0 - cosTetMaxNuc2)/screenZ
204         << " costm= " << cosTetMaxNuc2 << G4endl;
205  */
206  // scattering off nucleus
207  if(cosTetMaxNuc2 < 1.0) {
208    x  = 1.0 - cosTetMaxNuc2;
209    x1 = screenZ*formfactA;
210    x2 = 1.0 - x1; 
211    x3 = x/screenZ;
212    x4 = formfactA*x;
213    // low-energy limit
214    if(x3 < numlimit && x1 < numlimit) {
215      y = 0.5*x3*x3*(1.0 - 1.3333333*x3 + 1.5*x3*x3 - 1.5*x1
216                     + 3.0*x1*x1 + 2.666666*x3*x1)/(x2*x2*x2);
217      // high energy limit
218    } else if(x2 <= 0.0) {
219      x4 = x1*(1.0 + x3);
220      y  = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4);
221      // middle energy
222    } else {
223      y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4)) 
224           - x3/(1. + x3) - x4/(1. + x4))/(x2*x2); 
225    }
226    //G4cout << "y= " << y << " x1= " <<x1<<"  x2= " <<x2
227    //     <<"  x3= "<<x3<<"  x4= " << x4<<G4endl;
228    if(y < 0.0) {
229      nwarnings++;
230      if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) { 
231        G4cout << "Nuclear scattering <0 for L1 " << y
232               << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
233               << particle->GetParticleName() << G4endl;
234        G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
235               << " x= " << " x1= " << x1 << " x2= " << x2
236               << " x3= " << x3 << " x4= " << x4 <<G4endl;
237      }
238      y = 0.0;
239    }
240    xSection += y*targetZ; 
241  }
242  xSection *= kinFactor;
243  /*
244  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
245         << " screenZ= " << screenZ << " formF= " << formfactA
246         << " for " << particle->GetParticleName()
247         << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
248         << " x= " << x
249         << G4endl;
250  */
251  return xSection; 
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255
256G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
257                             const G4Track& track,
258                             G4PhysicsTable* theTable,
259                             G4double currentMinimalStep)
260{
261  G4double tlimit = currentMinimalStep;
262  const G4DynamicParticle* dp = track.GetDynamicParticle();
263  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
264  G4StepStatus stepStatus = sp->GetStepStatus();
265
266  // initialisation for 1st step 
267  if(stepStatus == fUndefined) {
268    inside = false;
269    SetupParticle(dp->GetDefinition());
270    theLambdaTable = theTable;
271  }
272
273  // initialisation for each step, lambda may be computed from scratch
274  preKinEnergy  = dp->GetKineticEnergy();
275  DefineMaterial(track.GetMaterialCutsCouple());
276  lambda0 = GetLambda(preKinEnergy);
277  currentRange = 
278    theManager->GetRangeFromRestricteDEDX(particle,preKinEnergy,currentCouple);
279
280  // extra check for abnormal situation
281  // this check needed to run MSC with eIoni and eBrem inactivated
282  if(tlimit > currentRange) tlimit = currentRange;
283
284  // stop here if small range particle
285  if(inside) return tlimit;   
286
287  // pre step
288  G4double presafety = sp->GetSafety();
289
290  // compute presafety again if presafety <= 0 and no boundary
291  // i.e. when it is needed for optimization purposes
292  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) 
293    presafety = ComputeSafety(sp->GetPosition(), tlimit); 
294  /*
295  G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= "
296         <<tlimit<<" safety= " << presafety
297         << " range= " <<currentRange<<G4endl;
298  */
299  // far from geometry boundary
300  if(currentRange < presafety) {
301    inside = true; 
302   
303    // limit mean scattering angle
304  } else {
305    G4double rlimit = facrange*lambda0;
306    G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
307    if(rcut > rlimit) rlimit = std::pow(rcut*rcut*rlimit,0.33333333);
308    rlimit = std::min(rlimit, facgeom*currentMaterial->GetRadlen());
309    if(rlimit < tlimit) tlimit = rlimit;
310  }
311  /*
312  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
313         << " L0= " << lambda0 << " R= " << currentRange
314         << "tlimit= " << tlimit 
315         << " currentMinimalStep= " << currentMinimalStep << G4endl;
316  */
317  return tlimit;
318}
319
320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
321
322G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
323{
324  tPathLength  = truelength;
325  zPathLength  = tPathLength;
326  lambdaeff    = lambda0;
327
328  if(lambda0 > 0.0) {
329    G4double tau = tPathLength/lambda0;
330    //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
331    //   << " lambda0= " << lambda0 << " tau= " << tau << G4endl;
332    // small step
333    if(tau < numlimit) {
334      zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
335
336      // medium step
337    } else {
338      //      zPathLength = lambda0*(1.0 - exp(-tPathLength/lambda0));
339      G4double e1 = 0.0;
340      if(currentRange > tPathLength) {
341        e1 = theManager->GetEnergy(particle,
342                                   currentRange-tPathLength,
343                                   currentCouple);
344      }
345      lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
346      zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
347    }
348  }
349  //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
350  return zPathLength;
351}
352
353//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
354
355G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
356{
357  // step defined other than transportation
358  if(geomStepLength == zPathLength) return tPathLength;
359
360  // step defined by transportation
361  tPathLength  = geomStepLength;
362  zPathLength  = geomStepLength;
363  G4double tau = zPathLength/lambdaeff;
364  tPathLength *= (1.0 + 0.5*tau + tau*tau/3.0); 
365
366  if(tau > numlimit) {
367    G4double e1 = 0.0;
368    if(currentRange > tPathLength) {
369      e1 = theManager->GetEnergy(particle,
370                                 currentRange-tPathLength,
371                                 currentCouple);
372    }
373    lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
374    tau = zPathLength/lambdaeff;
375
376    if(tau < 0.999999) tPathLength = -lambdaeff*log(1.0 - tau); 
377    else               tPathLength = currentRange;
378
379    if(tPathLength < zPathLength) tPathLength = zPathLength;
380  }
381  if(tPathLength > currentRange) tPathLength = currentRange;
382  //G4cout<<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
383  return tPathLength;
384}
385
386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
387
388void G4WentzelVIModel::SampleScattering(const G4DynamicParticle* dynParticle,
389                                        G4double safety)
390{
391  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
392  //     << particle->GetParticleName() << G4endl;
393  G4double kinEnergy = dynParticle->GetKineticEnergy();
394
395  // ignore scattering for zero step length and enegy below the limit
396  if(kinEnergy < lowEnergyLimit || tPathLength <= DBL_MIN) return;
397
398  G4double ekin = preKinEnergy;
399  if(ekin - kinEnergy > ekin*dtrl) {
400    ekin = 0.5*(preKinEnergy + kinEnergy);
401    lambdaeff = GetLambda(ekin);
402  } 
403 
404  G4double x1 = 0.5*tPathLength/lambdaeff;
405  G4double cut= (*currentCuts)[currentMaterialIndex];
406  /* 
407  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy<<" Eeff(MeV)= "<<ekin/MeV
408         << " L0= " << lambda0 << " Leff= " << lambdaeff
409         << " x1= " << x1 << " safety= " << safety << G4endl;
410  */
411
412  G4double xsec = 0.0;
413  G4bool largeAng = false;
414
415  // large scattering angle case
416  if(x1 > 0.5) {
417    x1 *= 0.5;
418    largeAng = true;
419
420    // normal case
421  } else {
422
423    // define threshold angle between single and multiple scattering
424    cosThetaMin = 1.0 - 3.0*x1;
425
426    // for low-energy e-,e+ no limit
427    SetupKinematic(ekin, cut);
428 
429    // recompute transport cross section
430    if(cosThetaMin > cosTetMaxNuc) {
431
432      xsec = ComputeXSectionPerVolume();
433
434      if(xtsec > 0.0) x1 = 0.5*tPathLength*xtsec;
435      else            x1 = 0.0;
436
437      /*     
438        G4cout << "cosTetMaxNuc= " << cosTetMaxNuc
439        << " cosThetaMin= " << cosThetaMin
440        << " cosThetaMax= " << cosThetaMax
441        << " cosTetMaxElec2= " << cosTetMaxElec2 << G4endl;
442        G4cout << "Recomputed xsec(1/mm)= " << xsec << " x1= " << x1 << G4endl;
443      */
444    }
445  }
446
447  // result of central part sampling
448  G4double z; 
449  do {
450    z = -x1*log(G4UniformRand());
451  } while (z > 1.0); 
452
453  // cost is sampled ------------------------------
454  G4double cost = 1.0 - 2.0*z;
455  //  if(cost < -1.0) cost = -1.0;
456  // else if(cost > 1.0) cost = 1.0;
457  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
458
459  G4double phi  = twopi*G4UniformRand();
460
461  G4double dirx = sint*cos(phi);
462  G4double diry = sint*sin(phi);
463 
464  //G4cout << "G4WentzelVIModel: step(mm)= " << tPathLength/mm
465  //     << " sint= " << sint << " cost= " << cost<< G4endl;
466 
467  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
468  G4ThreeVector newDirection(dirx,diry,cost);
469  G4ThreeVector temp(0.0,0.0,1.0);
470  G4ThreeVector pos(0.0,0.0,-zPathLength);
471  G4ThreeVector dir(0.0,0.0,1.0);
472  G4bool isscat = false;
473
474  // sample MSC scattering for large angle
475  // extra central scattering for half step
476  if(largeAng) {
477    isscat = true;
478    pos.setZ(-0.5*zPathLength);
479    do {
480      z = -x1*log(G4UniformRand());
481    } while (z > 1.0); 
482    cost = 1.0 - 2.0*z;
483    //if(std::fabs(cost) > 1.0) cost = 1.0;
484
485    sint = sqrt((1.0 - cost)*(1.0 + cost));
486    phi  = twopi*G4UniformRand();
487
488    // position and direction for secondary scattering
489    dir.set(sint*cos(phi),sint*sin(phi),cost);
490    pos += 0.5*dir*zPathLength;
491    x1 *= 2.0;
492  }
493
494  // sample Rutherford scattering for large angle
495  if(xsec > DBL_MIN) {
496    G4double t = tPathLength;
497    G4int nelm = currentMaterial->GetNumberOfElements();
498    const G4ElementVector* theElementVector = 
499      currentMaterial->GetElementVector();
500    do{
501      G4double x  = -log(G4UniformRand())/xsec;     
502      pos += dir*(zPathLength*std::min(x,t)/tPathLength);
503      t -= x;
504      if(t > 0.0) {
505        G4double zz1 = 1.0;
506        G4double qsec = G4UniformRand()*xsec;
507
508        // scattering off nucleus
509        G4int i = 0;
510        if(nelm > 1) {
511          for (; i<nelm; i++) {if(xsecn[i] >= qsec) break;}
512          if(i >= nelm) i = nelm - 1;
513        }
514        SetupTarget((*theElementVector)[i]->GetZ(), tkin);
515        G4double formf = formfactA;
516        G4double costm = cosTetMaxNuc2;
517        if(prob[i] > 0.0) {
518          if(G4UniformRand() <= prob[i]) {
519            formf = 0.0;
520            costm = cosTetMaxElec2;
521          }
522        }
523        if(cosThetaMin > costm) {
524
525          G4double w1 = 1. - cosThetaMin + screenZ;
526          G4double w2 = 1. - costm + screenZ;
527          G4double w3 = cosThetaMin - costm;
528          G4double grej, zz; 
529          do {
530            zz = w1*w2/(w1 + G4UniformRand()*w3) - screenZ;
531            grej = 1.0/(1.0 + formf*zz);
532          } while ( G4UniformRand() > grej*grej ); 
533          if(zz < 0.0) zz = 0.0;
534          else if(zz > 2.0) zz = 2.0;
535          zz1 = 1.0 - zz;
536        }
537        if(zz1 < 1.0) {
538          isscat = true;
539          //G4cout << "Rutherford zz1= " << zz1 << " t= " << t << G4endl;
540          sint = sqrt((1.0 - zz1)*(1.0 + zz1));
541          //G4cout << "sint= " << sint << G4endl;
542          phi  = twopi*G4UniformRand();
543          G4double vx1 = sint*cos(phi);
544          G4double vy1 = sint*sin(phi);
545          temp.set(vx1,vy1,zz1);
546          temp.rotateUz(dir);
547          dir = temp;
548        }
549      }
550    } while (t > 0.0); 
551  }
552  if(isscat) newDirection.rotateUz(dir);
553  newDirection.rotateUz(oldDirection);
554
555  //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
556  // end of sampling -------------------------------
557
558  fParticleChange->ProposeMomentumDirection(newDirection);
559
560  if (latDisplasment && safety > tlimitminfix) {
561    G4double rms = invsqrt12*sqrt(2.0*x1);
562    G4double dx = zPathLength*(0.5*dirx + rms*G4RandGauss::shoot(0.0,1.0));
563    G4double dy = zPathLength*(0.5*diry + rms*G4RandGauss::shoot(0.0,1.0));
564    G4double dz;
565    G4double d = (dx*dx + dy*dy)/(zPathLength*zPathLength);
566    if(d < numlimit)  dz = -0.5*zPathLength*d*(1.0 + 0.25*d);
567    else if(d < 1.0)  dz = -zPathLength*(1.0 - sqrt(1.0 - d));
568    else {
569      dx = dy = dz = 0.0;
570    }
571
572    temp.set(dx,dy,dz);
573    if(isscat) temp.rotateUz(dir);
574    pos += temp;
575   
576    pos.rotateUz(oldDirection);
577
578    G4double r = pos.mag();
579
580    /*   
581    G4cout << " r(mm)= " << r << " safety= " << safety
582           << " trueStep(mm)= " << tPathLength
583           << " geomStep(mm)= " << zPathLength
584           << G4endl;
585    */
586
587    if(r > tlimitminfix) {
588      pos /= r;
589      ComputeDisplacement(fParticleChange, pos, r, safety);
590    }
591  }
592  //G4cout << "G4WentzelVIModel::SampleScattering end" << G4endl;
593}
594
595//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
596
597G4double G4WentzelVIModel::ComputeXSectionPerVolume()
598{
599  const G4ElementVector* theElementVector = 
600    currentMaterial->GetElementVector();
601  const G4double* theAtomNumDensityVector = 
602    currentMaterial->GetVecNbOfAtomsPerVolume();
603  G4int nelm = currentMaterial->GetNumberOfElements();
604  if(nelm > nelments) {
605    nelments = nelm;
606    xsecn.resize(nelments);
607    prob.resize(nelments);
608  }
609
610  xtsec = 0.0;
611  G4double xs = 0.0;
612
613  for (G4int i=0; i<nelm; i++) {
614    SetupTarget((*theElementVector)[i]->GetZ(), tkin);
615    G4double density = theAtomNumDensityVector[i];
616    G4double cosnm = cosTetMaxNuc2;
617    G4double cosem = cosTetMaxElec2;
618
619    // recompute the angular limit
620    cosTetMaxNuc2  = std::max(cosnm,cosThetaMin); 
621    cosTetMaxElec2 = std::max(cosem,cosThetaMin); 
622    xtsec += ComputeTransportXSectionPerAtom()*density;
623    // return limit back
624    cosTetMaxElec2 = cosem;
625    cosTetMaxNuc2  = cosnm;
626
627    G4double esec = 0.0;
628    G4double nsec = 0.0;
629    G4double x1 = 1.0 - cosThetaMin + screenZ;
630    G4double f  = kinFactor*density; 
631
632    // scattering off electrons
633    if(cosThetaMin > cosem) {
634      esec = f*(cosThetaMin - cosem)/(x1*(1.0 - cosem + screenZ));
635    }
636
637    // scattering off nucleaus
638    if(cosThetaMin > cosnm) {
639
640      // Rutherford part
641      G4double s  = screenZ*formfactA;
642      G4double z1 = 1.0 - cosnm + screenZ;
643      G4double s1 = 1.0 - s;
644      G4double d  = s1/formfactA;
645
646      // check numerical limit
647      if(d < numlimit*x1) {
648        G4double x2 = x1*x1;
649        G4double z2 = z1*z1;
650        nsec = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/
651          (3.0*formfactA*formfactA);
652      } else {
653        G4double x2 = x1 + d;
654        G4double z2 = z1 + d;
655        nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
656      }
657      nsec *= f*targetZ;
658    }
659    nsec += esec;
660    if(nsec > 0.0) esec /= nsec;
661    xs += nsec;
662    xsecn[i] = xs;
663    prob[i]  = esec;
664    //G4cout << i << "  xs= " << xs << " cosThetaMin= " << cosThetaMin
665    //     << " costm= " << costm << G4endl;
666  }
667 
668  //G4cout << "ComputeXS result:  xsec(1/mm)= " << xs
669  //<< " txsec(1/mm)= " << xtsec <<G4endl;
670  return xs;
671}
672
673//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
674
675/*
676G4double G4MuMscModel::ComputeXSectionPerVolume()
677{
678  const G4ElementVector* theElementVector =
679    currentMaterial->GetElementVector();
680  const G4double* theAtomNumDensityVector =
681    currentMaterial->GetVecNbOfAtomsPerVolume();
682  size_t nelm = currentMaterial->GetNumberOfElements();
683
684  xsece1 = 0.0;
685  xsece2 = 0.0;
686  xsecn2 = 0.0;
687  zcorr  = 0.0;
688
689  G4double fac = coeff*chargeSquare*invbeta2/mom2;
690
691  for (size_t i=0; i<nelm; i++) {
692    const G4Element* elm = (*theElementVector)[i];
693    G4double Z = elm->GetZ();
694    SetupTarget(Z, tkin);
695    G4double den = fac*theAtomNumDensityVector[i]*Z;
696
697    G4double x  = 1.0 - cosThetaMin;
698    G4double x1 = x + screenZ;
699    G4double x2 = 1.0/(x1*x1);
700    G4double x3 = 1.0 + x*formfactA;
701   
702    //G4cout << "x= " << x << " den= " << den << " cosE= " << cosTetMaxElec << G4endl;
703    //G4cout << "cosThtaMin= " << cosThetaMin << G4endl;
704    //G4cout << "cosTetMaxNuc= " << cosTetMaxNuc << " q2Limit= " << q2Limit << G4endl;
705   
706    // scattering off electrons
707    if(cosTetMaxElec < cosThetaMin) {
708
709      // flat part
710      G4double s = den*x2*x;
711      xsece1 += s;
712      zcorr  += 0.5*x*s;
713
714      // Rutherford part
715      G4double z1 = 1.0 - cosTetMaxElec + screenZ;
716      G4double z2 = (cosThetaMin - cosTetMaxElec)/x1;
717      if(z2 < 0.2) s = z2*(x - 0.5*z2*(x - screenZ))/x1;
718      else         s = log(1.0 + z2)  - screenZ*z2/z1;
719      xsece2  += den*z2/z1;
720      zcorr   += den*s;
721    }
722    den *= Z;
723
724    //G4cout << "Z= " << Z<< " cosL= " << cosTetMaxNuc << " cosMin= " << cosThetaMin << G4endl;
725    // scattering off nucleaus
726    if(cosTetMaxNuc < cosThetaMin) {
727
728      // flat part
729      G4double s = den*x2*x/(x3*x3);
730      xsece1 += s;
731      zcorr  += 0.5*x*s;
732
733      // Rutherford part
734      s  = screenZ*formfactA;
735      G4double w  = 1.0 + 2.0*s;
736      G4double z1 = 1.0 - cosTetMaxNuc + screenZ;
737      G4double d  = (1.0 - s)/formfactA;
738      G4double x4 = x1 + d;
739      G4double z4 = z1 + d;
740      G4double t1 = 1.0/(x1*z1);
741      G4double t4 = 1.0/(x4*z4);
742      G4double w1 = cosThetaMin - cosTetMaxNuc;
743      G4double w2 = log(z1*x4/(x1*z4));
744
745      den *= w;     
746      xsecn2  += den*(w1*(t1 + t4) - 2.0*w2/d);
747      zcorr   += den*(w*w2 - w1*(screenZ*t1 + t4/formfactA));
748    }
749    xsece[i] = xsece2;
750    xsecn[i] = xsecn2;
751    //    G4cout << i << "  xsece2= " << xsece2 << "  xsecn2= " << xsecn2 << G4endl;
752  }
753  G4double xsec = xsece1 + xsece2 + xsecn2;
754 
755    //G4cout << "xsece1= " << xsece1 << "  xsece2= " << xsece2
756    //<< "  xsecn2= " << xsecn2
757        // << " zsec= " << zcorr*0.5*tPathLength << G4endl;
758  zcorr *= 0.5*tPathLength;
759
760  return xsec;
761}
762*/
763
764//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
765
766void G4WentzelVIModel::ComputeMaxElectronScattering(G4double cutEnergy)
767{
768  ecut = cutEnergy;
769  G4double tmax = tkin;
770  cosTetMaxElec = 1.0;
771  if(mass > MeV) {
772    G4double ratio = electron_mass_c2/mass;
773    G4double tau = tkin/mass;
774    tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
775      (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); 
776    cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
777  } else {
778
779    if(particle == theElectron) tmax *= 0.5;
780    G4double t = std::min(cutEnergy, tmax);
781    G4double mom21 = t*(t + 2.0*electron_mass_c2);
782    G4double t1 = tkin - t;
783    //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
784    //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
785    if(t1 > 0.0) {
786      G4double mom22 = t1*(t1 + 2.0*mass);
787      G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
788      if(ctm <  1.0) cosTetMaxElec = ctm;
789      if(ctm < -1.0) cosTetMaxElec = -1.0;
790    }
791  }
792  if(cosTetMaxElec < cosTetMaxNuc) cosTetMaxElec = cosTetMaxNuc;
793}
794
795//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.