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

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

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