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

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 17.5 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.60 2010/06/01 11:13:31 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-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// 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
42//              compute cross sections and sample scattering angle
43//
44//
45// Class Description:
46//
47// Implementation of the model of multiple scattering based on
48// G.Wentzel, Z. Phys. 40 (1927) 590.
49// H.W.Lewis, Phys Rev 78 (1950) 526.
50// J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
51// L.Urban, CERN-OPEN-2006-077.
52
53// -------------------------------------------------------------------
54//
55
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59#include "G4WentzelVIModel.hh"
60#include "Randomize.hh"
61#include "G4ParticleChangeForMSC.hh"
62#include "G4PhysicsTableHelper.hh"
63#include "G4ElementVector.hh"
64#include "G4ProductionCutsTable.hh"
65#include "G4LossTableManager.hh"
66#include "G4Pow.hh"
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70using namespace std;
71
72G4WentzelVIModel::G4WentzelVIModel(const G4String& nam) :
73  G4VMscModel(nam),
74  theLambdaTable(0),
75  numlimit(0.1),
76  currentCouple(0),
77  cosThetaMin(1.0),
78  isInitialized(false),
79  inside(false)
80{
81  invsqrt12 = 1./sqrt(12.);
82  tlimitminfix = 1.e-6*mm;
83  lowEnergyLimit = 1.0*eV;
84  particle = 0;
85  nelments = 5;
86  xsecn.resize(nelments);
87  prob.resize(nelments);
88  theManager = G4LossTableManager::Instance();
89  fG4pow = G4Pow::GetInstance();
90  wokvi = new G4WentzelOKandVIxSection();
91}
92
93//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94
95G4WentzelVIModel::~G4WentzelVIModel()
96{
97  delete wokvi;
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101
102void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
103                                  const G4DataVector& cuts)
104{
105  // reset parameters
106  SetupParticle(p);
107  currentRange = 0.0;
108  cosThetaMax = cos(PolarAngleLimit());
109  wokvi->Initialise(p, cosThetaMax);
110  /*
111  G4cout << "G4WentzelVIModel: factorA2(GeV^2) = " << factorA2/(GeV*GeV)
112         << "  1-cos(ThetaLimit)= " << 1 - cosThetaMax
113         << G4endl;
114  */
115  currentCuts = &cuts;
116
117  // set values of some data members
118  if(!isInitialized) {
119    isInitialized = true;
120    fParticleChange = GetParticleChangeForMSC();
121    InitialiseSafetyHelper();
122  }
123}
124
125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126
127G4double G4WentzelVIModel::ComputeCrossSectionPerAtom( 
128                             const G4ParticleDefinition* p,
129                             G4double kinEnergy,
130                             G4double Z, G4double,
131                             G4double cutEnergy, G4double)
132{
133  G4double xsec = 0.0;
134  if(p != particle) { SetupParticle(p); }
135  if(kinEnergy < lowEnergyLimit) { return xsec; }
136  DefineMaterial(CurrentCouple());
137  cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
138  if(cosTetMaxNuc < 1.0) {
139    cosTetMaxNuc = wokvi->SetupTarget(G4int(Z), cutEnergy);
140    xsec = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
141  /*     
142    G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy
143           << " 1-cosN= " << 1 - costm << " xsec(bn)= " << xsec/barn
144           << " " << particle->GetParticleName() << G4endl;
145  */
146  }
147  return xsec;
148}
149
150//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151
152G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
153                             const G4Track& track,
154                             G4PhysicsTable* theTable,
155                             G4double currentMinimalStep)
156{
157  G4double tlimit = currentMinimalStep;
158  const G4DynamicParticle* dp = track.GetDynamicParticle();
159  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
160  G4StepStatus stepStatus = sp->GetStepStatus();
161  //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
162  //     << stepStatus << G4endl;
163
164  // initialisation for 1st step 
165  if(stepStatus == fUndefined) {
166    inside = false;
167    SetupParticle(dp->GetDefinition());
168  }
169
170  // initialisation for each step, lambda may be computed from scratch
171  preKinEnergy  = dp->GetKineticEnergy();
172  DefineMaterial(track.GetMaterialCutsCouple());
173  theLambdaTable = theTable;
174  lambdaeff = GetLambda(preKinEnergy);
175  currentRange = 
176    theManager->GetRangeFromRestricteDEDX(particle,preKinEnergy,currentCouple);
177  cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
178
179  // extra check for abnormal situation
180  // this check needed to run MSC with eIoni and eBrem inactivated
181  if(tlimit > currentRange) { tlimit = currentRange; }
182
183  // stop here if small range particle
184  if(inside) { return tlimit; }
185
186  // pre step
187  G4double presafety = sp->GetSafety();
188
189  // compute presafety again if presafety <= 0 and no boundary
190  // i.e. when it is needed for optimization purposes
191  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
192    presafety = ComputeSafety(sp->GetPosition(), tlimit); 
193  }
194  /*
195  G4cout << "e(MeV)= " << preKinEnergy/MeV
196         << "  " << particle->GetParticleName()
197         << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
198         << " R(mm)= " <<currentRange/mm
199         << " L0(mm^-1)= " << lambdaeff*mm
200         <<G4endl;
201  */
202  // far from geometry boundary
203  if(currentRange < presafety) {
204    inside = true; 
205    return tlimit;
206  }
207
208  // natural limit for high energy
209  G4double rlimit = std::max(facrange*currentRange, 
210                             0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
211
212  // low-energy e-
213  if(cosThetaMax > cosTetMaxNuc) {
214    rlimit = std::min(rlimit, facsafety*presafety);
215  }
216   
217  // cut correction
218  G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
219  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= " << presafety
220  // << " 1-cosThetaMax= " <<1-cosThetaMax << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
221  // << G4endl;
222  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
223
224  if(rlimit < tlimit) { tlimit = rlimit; }
225
226  tlimit = std::max(tlimit, tlimitminfix);
227
228  // step limit in infinite media
229  tlimit = std::min(tlimit, 20*currentMaterial->GetRadlen());
230  /* 
231  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
232         << " L0= " << lambdaeff << " R= " << currentRange
233         << "tlimit= " << tlimit 
234         << " currentMinimalStep= " << currentMinimalStep << G4endl;
235  */
236  return tlimit;
237}
238
239//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240
241G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
242{
243  tPathLength  = truelength;
244  zPathLength  = tPathLength;
245
246  if(lambdaeff > 0.0) {
247    G4double tau = tPathLength/lambdaeff;
248    //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
249    //   << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
250    // small step
251    if(tau < numlimit) {
252      zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
253
254      // medium step
255    } else {
256      G4double e1 = 0.0;
257      if(currentRange > tPathLength) {
258        e1 = theManager->GetEnergy(particle,
259                                   currentRange-tPathLength,
260                                   currentCouple);
261      }
262      e1 = 0.5*(e1 + preKinEnergy);
263      cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
264      lambdaeff = GetLambda(e1);
265      zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
266    }
267  }
268  //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
269  return zPathLength;
270}
271
272//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
273
274G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
275{
276  // initialisation of single scattering x-section
277  xtsec = 0.0;
278
279  // pathalogical case
280  if(lambdaeff <= 0.0) { 
281    zPathLength = geomStepLength;
282    tPathLength = geomStepLength;
283    return tPathLength;
284  }
285
286  G4double tau = geomStepLength/lambdaeff;
287
288  // step defined by transportation
289  if(geomStepLength != zPathLength) { 
290
291    // step defined by transportation
292    zPathLength = geomStepLength;
293    tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); 
294
295    // energy correction for a big step
296    if(tau > numlimit) {
297      G4double e1 = 0.0;
298      if(currentRange > tPathLength) {
299        e1 = theManager->GetEnergy(particle,
300                                   currentRange-tPathLength,
301                                   currentCouple);
302      }
303      e1 = 0.5*(e1 + preKinEnergy);
304      cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
305      lambdaeff = GetLambda(e1);
306      tau = zPathLength/lambdaeff;
307     
308      if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); } 
309      else               { tPathLength = currentRange; }
310    }
311  }
312
313  // check of step length
314  // define threshold angle between single and multiple scattering
315  cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff;
316
317  // recompute transport cross section - do not change energy
318  // anymore - cannot be applied for big steps
319  if(cosThetaMin > cosTetMaxNuc) {
320
321    // new computation
322    G4double xsec = ComputeXSectionPerVolume();
323    //G4cout << "%%%% xsec= " << xsec << "  xtsec= " << xtsec << G4endl;
324    if(xtsec > 0.0) {
325      if(xsec > 0.0) { lambdaeff = 1./xsec; }
326      else           { lambdaeff = DBL_MAX; }
327
328      tau = zPathLength*xsec;
329      if(tau < numlimit) { tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); } 
330      else if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); } 
331      else               { tPathLength = currentRange; }
332    } 
333  } 
334
335  if(tPathLength > currentRange) { tPathLength = currentRange; }
336  if(tPathLength < zPathLength)  { tPathLength = zPathLength; }
337  /*   
338  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
339         <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
340  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
341         << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
342         << " e(MeV)= " << preKinEnergy/MeV << G4endl;
343  */
344  return tPathLength;
345}
346
347//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
348
349void G4WentzelVIModel::SampleScattering(const G4DynamicParticle* dynParticle,
350                                        G4double safety)
351{
352  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
353  //     << particle->GetParticleName() << G4endl;
354
355  // ignore scattering for zero step length and energy below the limit
356  if(dynParticle->GetKineticEnergy() < lowEnergyLimit || 
357     tPathLength <= DBL_MIN || lambdaeff <= DBL_MIN) 
358    { return; }
359 
360  G4double invlambda = 0.0;
361  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
362
363  // use average kinetic energy over the step
364  G4double cut = (*currentCuts)[currentMaterialIndex];
365  /*     
366  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
367         << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
368         << " x1= " <<  tPathLength*invlambda << " safety= " << safety << G4endl;
369  */
370
371  G4double length = tPathLength;
372  G4double lengthlim = tPathLength*1.e-6;
373
374  // step limit due msc
375  G4double x0 = length;
376  // large scattering angle case - two step approach
377  if(tPathLength*invlambda > 0.5 && length > tlimitminfix) { x0 *= 0.5; } 
378
379  // step limit due single scattering
380  G4double x1 = length;
381  if(xtsec > 0.0) { x1 = -log(G4UniformRand())/xtsec; }
382
383  const G4ElementVector* theElementVector = 
384    currentMaterial->GetElementVector();
385  G4int nelm = currentMaterial->GetNumberOfElements();
386
387  // geometry
388  G4double sint, cost, phi;
389  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
390  G4ThreeVector temp(0.0,0.0,1.0);
391
392  // current position and direction relative to the end point
393  // because of magnetic field geometry is computed relatively to the
394  // end point of the step
395  G4ThreeVector dir(0.0,0.0,1.0);
396  G4ThreeVector pos(0.0,0.0,-zPathLength);
397  G4double mscfac = zPathLength/tPathLength;
398
399  // start a loop
400  do {
401    G4double step = x0; 
402    G4bool singleScat = false;
403
404    // single scattering case
405    if(x1 < x0) { 
406      step = x1;
407      singleScat = true;
408    }
409
410    // new position
411    pos += step*mscfac*dir;
412
413    // added multiple scattering
414    G4double z; 
415    G4double tet2 = step*invlambda; 
416    do { z = -tet2*log(G4UniformRand()); } while (z >= 1.0); 
417
418    cost = 1.0 - 2.0*z;
419    sint = sqrt((1.0 - cost)*(1.0 + cost));
420    phi  = twopi*G4UniformRand();
421    G4double vx1 = sint*cos(phi);
422    G4double vy1 = sint*sin(phi);
423
424    // lateral displacement 
425    if (latDisplasment && safety > tlimitminfix) {
426      G4double rms = invsqrt12*sqrt(2.0*tet2);
427      G4double dx = step*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
428      G4double dy = step*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
429      G4double dz;
430      G4double d = (dx*dx + dy*dy)/(step*step);
431      if(d < numlimit)  { dz = -0.5*step*d*(1.0 + 0.25*d); }
432      else if(d < 1.0)  { dz = -step*(1.0 - sqrt(1.0 - d));}
433      else              { dx = dy = dz = 0.0; }
434
435      // change position
436      temp.set(dx,dy,dz);
437      temp.rotateUz(dir); 
438      pos += temp;
439    }
440
441    // direction is changed
442    temp.set(vx1,vy1,cost);
443    temp.rotateUz(dir);
444    dir = temp;
445
446    if(singleScat) {
447
448      // select element
449      G4int i = 0;
450      if(nelm > 1) {
451        G4double qsec = G4UniformRand()*xtsec;
452        for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
453        if(i >= nelm) { i = nelm - 1; }
454      }
455      G4double cosTetM = 
456        wokvi->SetupTarget(G4int((*theElementVector)[i]->GetZ()), cut);
457      temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
458      temp.rotateUz(dir);
459
460      // renew direction
461      dir = temp;
462
463      // new single scatetring
464      x1 = -log(G4UniformRand())/xtsec; 
465    }
466
467    // update step
468    length -= step;
469
470  } while (length > lengthlim);
471   
472  dir.rotateUz(oldDirection);
473  pos.rotateUz(oldDirection);
474
475  //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
476  // end of sampling -------------------------------
477
478  fParticleChange->ProposeMomentumDirection(dir);
479
480  // lateral displacement 
481  if (latDisplasment) {
482    G4double r = pos.mag();
483
484    /*   
485    G4cout << " r(mm)= " << r << " safety= " << safety
486           << " trueStep(mm)= " << tPathLength
487           << " geomStep(mm)= " << zPathLength
488           << G4endl;
489    */
490
491    if(r > tlimitminfix) {
492      pos /= r;
493      ComputeDisplacement(fParticleChange, pos, r, safety);
494    }
495  }
496  //G4cout << "G4WentzelVIModel::SampleScattering end" << G4endl;
497}
498
499//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
500
501G4double G4WentzelVIModel::ComputeXSectionPerVolume()
502{
503  // prepare recomputation of x-sections
504  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
505  const G4double* theAtomNumDensityVector = 
506    currentMaterial->GetVecNbOfAtomsPerVolume();
507  G4int nelm = currentMaterial->GetNumberOfElements();
508  if(nelm > nelments) {
509    nelments = nelm;
510    xsecn.resize(nelm);
511    prob.resize(nelm);
512  }
513  G4double cut = (*currentCuts)[currentMaterialIndex];
514  cosTetMaxNuc = wokvi->GetCosThetaNuc();
515
516  // check consistency
517  xtsec = 0.0;
518  if(cosTetMaxNuc > cosThetaMin) { return 0.0; }
519
520  // loop over elements
521  G4double xs = 0.0;
522  for (G4int i=0; i<nelm; ++i) {
523    G4double costm = 
524      wokvi->SetupTarget(G4int((*theElementVector)[i]->GetZ()), cut);
525    G4double density = theAtomNumDensityVector[i];
526
527    G4double esec = 0.0;
528    if(costm < cosThetaMin) { 
529
530      // recompute the transport x-section
531      xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosThetaMin);
532
533      // recompute the total x-section
534      G4double nsec = wokvi->ComputeNuclearCrossSection(cosThetaMin, costm);
535      esec = wokvi->ComputeElectronCrossSection(cosThetaMin, costm);
536      nsec += esec;
537      if(nsec > 0.0) { esec /= nsec; }
538      xtsec += nsec*density;
539    }
540    xsecn[i] = xtsec;
541    prob[i]  = esec;
542    //G4cout << i << "  xs= " << xs << " xtsec= " << xtsec << " 1-cosThetaMin= " << 1-cosThetaMin
543    //     << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
544  }
545 
546  //G4cout << "ComputeXS result:  xsec(1/mm)= " << xs
547  //     << " txsec(1/mm)= " << xtsec <<G4endl;
548  return xs;
549}
550
551//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.