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

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

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