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

Last change on this file since 1183 was 1055, checked in by garnier, 15 years ago

maj sur la beta de geant 4.9.3
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.32 2009/05/10 16:09:29 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
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 = DBL_MIN;
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  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
161  SetupKinematic(ekin, cutEnergy);
162  SetupTarget(Z, ekin);
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  if(kinEnergy <= DBL_MIN || tPathLength <= DBL_MIN) return;
395
396  G4double ekin = preKinEnergy;
397  if(ekin - kinEnergy > ekin*dtrl) {
398    ekin = 0.5*(preKinEnergy + kinEnergy);
399    lambdaeff = GetLambda(ekin);
400  } 
401 
402  G4double x1 = 0.5*tPathLength/lambdaeff;
403  G4double cut= (*currentCuts)[currentMaterialIndex];
404  /* 
405  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy<<" Eeff(MeV)= "<<ekin/MeV
406         << " L0= " << lambda0 << " Leff= " << lambdaeff
407         << " x1= " << x1 << " safety= " << safety << G4endl;
408  */
409
410  G4double xsec = 0.0;
411  G4bool largeAng = false;
412
413  // large scattering angle case
414  if(x1 > 0.5) {
415    x1 *= 0.5;
416    largeAng = true;
417
418    // normal case
419  } else {
420
421    // define threshold angle between single and multiple scattering
422    cosThetaMin = 1.0 - 3.0*x1;
423
424    // for low-energy e-,e+ no limit
425    ekin = std::max(ekin, lowEnergyLimit);
426    SetupKinematic(ekin, cut);
427 
428    // recompute transport cross section
429    if(cosThetaMin > cosTetMaxNuc) {
430
431      xsec = ComputeXSectionPerVolume();
432
433      if(xtsec > DBL_MIN) x1 = 0.5*tPathLength*xtsec;
434      else                x1 = 0.0;
435
436      /*     
437        G4cout << "cosTetMaxNuc= " << cosTetMaxNuc
438        << " cosThetaMin= " << cosThetaMin
439        << " cosThetaMax= " << cosThetaMax
440        << " cosTetMaxElec2= " << cosTetMaxElec2 << G4endl;
441        G4cout << "Recomputed xsec(1/mm)= " << xsec << " x1= " << x1 << G4endl;
442      */
443    }
444  }
445
446  // result of central part sampling
447  G4double z; 
448  do {
449    z = -x1*log(G4UniformRand());
450  } while (z > 1.0); 
451
452  // cost is sampled ------------------------------
453  G4double cost = 1.0 - 2.0*z;
454  if(cost < -1.0) cost = -1.0;
455  else if(cost > 1.0) cost = 1.0;
456  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
457
458  G4double phi  = twopi*G4UniformRand();
459
460  G4double dirx = sint*cos(phi);
461  G4double diry = sint*sin(phi);
462 
463  //G4cout << "G4WentzelVIModel: step(mm)= " << tPathLength/mm
464  //     << " sint= " << sint << " cost= " << cost<< G4endl;
465 
466  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
467  G4ThreeVector newDirection(dirx,diry,cost);
468  G4ThreeVector temp(0.0,0.0,1.0);
469  G4ThreeVector pos(0.0,0.0,-zPathLength);
470  G4ThreeVector dir(0.0,0.0,1.0);
471  G4bool isscat = false;
472
473  // sample MSC scattering for large angle
474  // extra central scattering for holf step
475  if(largeAng) {
476    isscat = true;
477    pos.setZ(-0.5*zPathLength);
478    do {
479      z = -x1*log(G4UniformRand());
480    } while (z > 1.0); 
481    cost = 1.0 - 2.0*z;
482    if(std::abs(cost) > 1.0) cost = 1.0;
483
484    sint = sqrt((1.0 - cost)*(1.0 + cost));
485    phi  = twopi*G4UniformRand();
486
487    // position and direction for secondary scattering
488    dir.set(sint*cos(phi),sint*sin(phi),cost);
489    pos += 0.5*dir*zPathLength;
490    x1 *= 2.0;
491  }
492
493  // sample Reserford scattering for large angle
494  if(xsec > DBL_MIN) {
495    G4double t = tPathLength;
496    G4int nelm = currentMaterial->GetNumberOfElements();
497    const G4ElementVector* theElementVector = 
498      currentMaterial->GetElementVector();
499    do{
500      G4double x  = -log(G4UniformRand())/xsec;     
501      pos += dir*(zPathLength*std::min(x,t)/tPathLength);
502      t -= x;
503      if(t > 0.0) {
504        G4double zz1 = 1.0;
505        G4double qsec = G4UniformRand()*xsec;
506
507        // scattering off nucleus
508        G4int i = 0;
509        if(nelm > 1) {
510          for (; i<nelm; i++) {if(xsecn[i] >= qsec) break;}
511          if(i >= nelm) i = nelm - 1;
512        }
513        SetupTarget((*theElementVector)[i]->GetZ(), tkin);
514        G4double formf = formfactA;
515        G4double costm = cosTetMaxNuc2;
516        if(prob[i] > 0.0) {
517          if(G4UniformRand() <= prob[i]) {
518            formf = 0.0;
519            costm = cosTetMaxElec2;
520          }
521        }
522        if(cosThetaMin > costm) {
523
524          G4double w1 = 1. - cosThetaMin + screenZ;
525          G4double w2 = 1. - costm + screenZ;
526          G4double w3 = cosThetaMin - costm;
527          G4double grej, zz; 
528          do {
529            zz = w1*w2/(w1 + G4UniformRand()*w3) - screenZ;
530            grej = 1.0/(1.0 + formf*zz);
531          } while ( G4UniformRand() > grej*grej ); 
532          if(zz < 0.0) zz = 0.0;
533          else if(zz > 2.0) zz = 2.0;
534          zz1 = 1.0 - zz;
535        }
536        if(zz1 < 1.0) {
537          isscat = true;
538          //G4cout << "Reserford zz1= " << zz1 << " t= " << t << G4endl;
539          sint = sqrt((1.0 - zz1)*(1.0 + zz1));
540          //G4cout << "sint= " << sint << G4endl;
541          phi  = twopi*G4UniformRand();
542          G4double vx1 = sint*cos(phi);
543          G4double vy1 = sint*sin(phi);
544          temp.set(vx1,vy1,zz1);
545          temp.rotateUz(dir);
546          dir = temp;
547        }
548      }
549    } while (t > 0.0); 
550  }
551  if(isscat) newDirection.rotateUz(dir);
552  newDirection.rotateUz(oldDirection);
553
554  //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
555  // end of sampling -------------------------------
556
557  fParticleChange->ProposeMomentumDirection(newDirection);
558
559  if (latDisplasment && safety > tlimitminfix) {
560    G4double rms = invsqrt12*sqrt(2.0*x1);
561    G4double dx = zPathLength*(0.5*dirx + rms*G4RandGauss::shoot(0.0,1.0));
562    G4double dy = zPathLength*(0.5*diry + rms*G4RandGauss::shoot(0.0,1.0));
563    G4double dz;
564    G4double d = (dx*dx + dy*dy)/(zPathLength*zPathLength);
565    if(d < numlimit)  dz = -0.5*zPathLength*d*(1.0 + 0.25*d);
566    else if(d < 1.0)  dz = -zPathLength*(1.0 - sqrt(1.0 - d));
567    else {
568      dx = dy = dz = 0.0;
569    }
570
571    temp.set(dx,dy,dz);
572    if(isscat) temp.rotateUz(dir);
573    pos += temp;
574   
575    pos.rotateUz(oldDirection);
576
577    G4double r = pos.mag();
578
579    /*   
580    G4cout << " r(mm)= " << r << " safety= " << safety
581           << " trueStep(mm)= " << tPathLength
582           << " geomStep(mm)= " << zPathLength
583           << G4endl;
584    */
585
586    if(r > tlimitminfix) {
587      pos /= r;
588      ComputeDisplacement(fParticleChange, pos, r, safety);
589    }
590  }
591  //G4cout << "G4WentzelVIModel::SampleScattering end" << G4endl;
592}
593
594//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
595
596G4double G4WentzelVIModel::ComputeXSectionPerVolume()
597{
598  const G4ElementVector* theElementVector = 
599    currentMaterial->GetElementVector();
600  const G4double* theAtomNumDensityVector = 
601    currentMaterial->GetVecNbOfAtomsPerVolume();
602  G4int nelm = currentMaterial->GetNumberOfElements();
603  if(nelm > nelments) {
604    nelments = nelm;
605    xsecn.resize(nelments);
606    prob.resize(nelments);
607  }
608
609  xtsec = 0.0;
610  G4double xs = 0.0;
611
612  for (G4int i=0; i<nelm; i++) {
613    SetupTarget((*theElementVector)[i]->GetZ(), tkin);
614    G4double density = theAtomNumDensityVector[i];
615    G4double cosnm = cosTetMaxNuc2;
616    G4double cosem = cosTetMaxElec2;
617
618    // recompute the angular limit
619    cosTetMaxNuc2  = std::max(cosnm,cosThetaMin); 
620    cosTetMaxElec2 = std::max(cosem,cosThetaMin); 
621    xtsec += ComputeTransportXSectionPerAtom()*density;
622    // return limit back
623    cosTetMaxElec2 = cosem;
624    cosTetMaxNuc2  = cosnm;
625
626    G4double esec = 0.0;
627    G4double nsec = 0.0;
628    G4double x1 = 1.0 - cosThetaMin + screenZ;
629    G4double f  = kinFactor*density; 
630
631    // scattering off electrons
632    if(cosThetaMin > cosem) {
633      esec = f*(cosThetaMin - cosem)/(x1*(1.0 - cosem + screenZ));
634    }
635
636    // scattering off nucleaus
637    if(cosThetaMin > cosnm) {
638
639      // Reserford part
640      G4double s  = screenZ*formfactA;
641      G4double z1 = 1.0 - cosnm + screenZ;
642      G4double s1 = 1.0 - s;
643      G4double d  = s1/formfactA;
644
645      // check numerical limit
646      if(d < numlimit*x1) {
647        G4double x2 = x1*x1;
648        G4double z2 = z1*z1;
649        nsec = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/
650          (3.0*formfactA*formfactA);
651      } else {
652        G4double x2 = x1 + d;
653        G4double z2 = z1 + d;
654        nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
655      }
656      nsec *= f*targetZ;
657    }
658    nsec += esec;
659    if(nsec > 0.0) esec /= nsec;
660    xs += nsec;
661    xsecn[i] = xs;
662    prob[i]  = esec;
663    //G4cout << i << "  xs= " << xs << " cosThetaMin= " << cosThetaMin
664    //     << " costm= " << costm << G4endl;
665  }
666 
667  //G4cout << "ComputeXS result:  xsec(1/mm)= " << xs
668  //<< " txsec(1/mm)= " << xtsec <<G4endl;
669  return xs;
670}
671
672//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
673
674/*
675G4double G4MuMscModel::ComputeXSectionPerVolume()
676{
677  const G4ElementVector* theElementVector =
678    currentMaterial->GetElementVector();
679  const G4double* theAtomNumDensityVector =
680    currentMaterial->GetVecNbOfAtomsPerVolume();
681  size_t nelm = currentMaterial->GetNumberOfElements();
682
683  xsece1 = 0.0;
684  xsece2 = 0.0;
685  xsecn2 = 0.0;
686  zcorr  = 0.0;
687
688  G4double fac = coeff*chargeSquare*invbeta2/mom2;
689
690  for (size_t i=0; i<nelm; i++) {
691    const G4Element* elm = (*theElementVector)[i];
692    G4double Z = elm->GetZ();
693    SetupTarget(Z, tkin);
694    G4double den = fac*theAtomNumDensityVector[i]*Z;
695
696    G4double x  = 1.0 - cosThetaMin;
697    G4double x1 = x + screenZ;
698    G4double x2 = 1.0/(x1*x1);
699    G4double x3 = 1.0 + x*formfactA;
700   
701    //G4cout << "x= " << x << " den= " << den << " cosE= " << cosTetMaxElec << G4endl;
702    //G4cout << "cosThtaMin= " << cosThetaMin << G4endl;
703    //G4cout << "cosTetMaxNuc= " << cosTetMaxNuc << " q2Limit= " << q2Limit << G4endl;
704   
705    // scattering off electrons
706    if(cosTetMaxElec < cosThetaMin) {
707
708      // flat part
709      G4double s = den*x2*x;
710      xsece1 += s;
711      zcorr  += 0.5*x*s;
712
713      // Reserford part
714      G4double z1 = 1.0 - cosTetMaxElec + screenZ;
715      G4double z2 = (cosThetaMin - cosTetMaxElec)/x1;
716      if(z2 < 0.2) s = z2*(x - 0.5*z2*(x - screenZ))/x1;
717      else         s = log(1.0 + z2)  - screenZ*z2/z1;
718      xsece2  += den*z2/z1;
719      zcorr   += den*s;
720    }
721    den *= Z;
722
723    //G4cout << "Z= " << Z<< " cosL= " << cosTetMaxNuc << " cosMin= " << cosThetaMin << G4endl;
724    // scattering off nucleaus
725    if(cosTetMaxNuc < cosThetaMin) {
726
727      // flat part
728      G4double s = den*x2*x/(x3*x3);
729      xsece1 += s;
730      zcorr  += 0.5*x*s;
731
732      // Reserford part
733      s  = screenZ*formfactA;
734      G4double w  = 1.0 + 2.0*s;
735      G4double z1 = 1.0 - cosTetMaxNuc + screenZ;
736      G4double d  = (1.0 - s)/formfactA;
737      G4double x4 = x1 + d;
738      G4double z4 = z1 + d;
739      G4double t1 = 1.0/(x1*z1);
740      G4double t4 = 1.0/(x4*z4);
741      G4double w1 = cosThetaMin - cosTetMaxNuc;
742      G4double w2 = log(z1*x4/(x1*z4));
743
744      den *= w;     
745      xsecn2  += den*(w1*(t1 + t4) - 2.0*w2/d);
746      zcorr   += den*(w*w2 - w1*(screenZ*t1 + t4/formfactA));
747    }
748    xsece[i] = xsece2;
749    xsecn[i] = xsecn2;
750    //    G4cout << i << "  xsece2= " << xsece2 << "  xsecn2= " << xsecn2 << G4endl;
751  }
752  G4double xsec = xsece1 + xsece2 + xsecn2;
753 
754    //G4cout << "xsece1= " << xsece1 << "  xsece2= " << xsece2
755    //<< "  xsecn2= " << xsecn2
756        // << " zsec= " << zcorr*0.5*tPathLength << G4endl;
757  zcorr *= 0.5*tPathLength;
758
759  return xsec;
760}
761*/
762
763//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
764
765void G4WentzelVIModel::ComputeMaxElectronScattering(G4double cutEnergy)
766{
767  ecut = cutEnergy;
768  G4double tmax = tkin;
769  cosTetMaxElec = 1.0;
770  if(mass > MeV) {
771    G4double ratio = electron_mass_c2/mass;
772    G4double tau = tkin/mass;
773    tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
774      (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); 
775    cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
776  } else {
777
778    if(particle == theElectron) tmax *= 0.5;
779    G4double t = std::min(cutEnergy, tmax);
780    G4double mom21 = t*(t + 2.0*electron_mass_c2);
781    G4double t1 = tkin - t;
782    //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
783    //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
784    if(t1 > 0.0) {
785      G4double mom22 = t1*(t1 + 2.0*mass);
786      G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
787      if(ctm <  1.0) cosTetMaxElec = ctm;
788      if(ctm < -1.0) cosTetMaxElec = -1.0;
789    }
790  }
791  if(cosTetMaxElec < cosTetMaxNuc) cosTetMaxElec = cosTetMaxNuc;
792}
793
794//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.