source: trunk/source/processes/electromagnetic/muons/src/G4MuMscModel.cc @ 961

Last change on this file since 961 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 18.0 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: G4MuMscModel.cc,v 1.6 2007/11/11 17:40:48 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:   G4MuMscModel
35//
36// Author:      Laszlo Mu
37//
38// Creation date: 03.03.2001
39//
40// Modifications:
41//
42// 27-03-03 Move model part from G4MultipleScattering80 (V.Ivanchenko)
43//
44
45// Class Description:
46//
47// Implementation of the model of multiple scattering based on
48// H.W.Lewis Phys Rev 78 (1950) 526 and others
49
50// -------------------------------------------------------------------
51//
52
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57#include "G4MuMscModel.hh"
58#include "Randomize.hh"
59#include "G4Electron.hh"
60#include "G4LossTableManager.hh"
61#include "G4ParticleChangeForMSC.hh"
62#include "G4TransportationManager.hh"
63#include "G4SafetyHelper.hh"
64#include "G4eCoulombScatteringModel.hh"
65#include "G4PhysicsTableHelper.hh"
66#include "G4ElementVector.hh"
67#include "G4ProductionCutsTable.hh"
68#include "G4PhysicsLogVector.hh"
69//#include "G4Poisson.hh"
70
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72
73using namespace std;
74
75G4MuMscModel::G4MuMscModel(G4double frange, 
76                           G4double thetaMax, 
77                           G4double tMax, 
78                           const G4String& nam)
79  : G4eCoulombScatteringModel(0.0,thetaMax,false,tMax,nam),
80    theLambdaTable(0),
81    theLambda2Table(0),
82    dtrl(0.05),
83    facrange(frange),
84    thetaLimit(thetaMax),
85    numlimit(0.2),
86    lowBinEnergy(keV),
87    highBinEnergy(PeV),
88    nbins(60),
89    nwarnings(0),
90    nwarnlimit(50),
91    currentCouple(0),
92    isInitialized(false),
93    buildTables(true),
94    newrun(true),
95    inside(false)
96{
97  invsqrt12 = 1./sqrt(12.);
98  tlimitminfix = 1.e-6*mm;
99  theManager = G4LossTableManager::Instance(); 
100}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103
104G4MuMscModel::~G4MuMscModel()
105{}
106
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108
109void G4MuMscModel::Initialise(const G4ParticleDefinition* p,
110                              const G4DataVector& cuts)
111{
112  SetupParticle(p);
113  newrun = true;
114  xSection = currentRange = targetZ = ecut = tkin = 0.0;
115  // set values of some data members
116  if(!isInitialized) {
117    isInitialized = true;
118    if(p->GetParticleName() == "GenericIon") buildTables = false;
119
120    if (pParticleChange)
121      fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
122    else
123      fParticleChange = new G4ParticleChangeForMSC();
124
125    safetyHelper = G4TransportationManager::GetTransportationManager()
126      ->GetSafetyHelper();
127    safetyHelper->InitialiseHelper();
128  }
129  G4eCoulombScatteringModel::Initialise(p, cuts);
130  currentCuts = &cuts;
131  if(buildTables)
132    theLambda2Table = G4PhysicsTableHelper::PreparePhysicsTable(theLambda2Table);
133}
134
135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136
137void G4MuMscModel::BuildTables()
138{
139  //G4cout << "G4MuMscModel::BuildTables flags newrun= " << newrun
140  //     << "  buildTables= " << buildTables << G4endl;
141  newrun = false;
142  if(!buildTables) return;
143
144  // Access to materials
145  const G4ProductionCutsTable* theCoupleTable=
146        G4ProductionCutsTable::GetProductionCutsTable();
147  size_t numOfCouples = theCoupleTable->GetTableSize();
148  G4double e, s, cut;
149
150  for(size_t i=0; i<numOfCouples; i++) {
151
152    if (theLambda2Table->GetFlag(i)) {
153
154      // create physics vector and fill it
155      DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i));
156      cut = (*currentCuts)[currentMaterialIndex];
157      G4PhysicsVector* aVector =
158        new G4PhysicsLogVector(lowBinEnergy, highBinEnergy, nbins);
159      for(G4int j=0; j<nbins; j++) {
160        e = aVector->GetLowEdgeEnergy(j);
161        s = ComputeLambda2(e, cut);
162        //G4cout << j << "  " << currentCouple->GetMaterial()->GetName()
163        //       << "  e(MeV)= " << e << " cut(MeV)= " << cut
164        //       << " L2= " << s << G4endl; 
165        aVector->PutValue(j, s);
166      }
167     
168      G4PhysicsTableHelper::SetPhysicsVector(theLambda2Table, i, aVector);
169    }
170  }
171}
172
173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174
175G4double G4MuMscModel::ComputeCrossSectionPerAtom( 
176                             const G4ParticleDefinition* p,
177                             G4double kinEnergy,
178                             G4double Z, G4double A,
179                             G4double cutEnergy, G4double)
180{
181  if(p == particle && kinEnergy == tkin && Z == targetZ &&
182     cutEnergy == ecut) return xSection;
183  ecut = cutEnergy;
184  xSection = 0.0;
185  SetupParticle(p);
186  G4double ekin = std::max(keV, kinEnergy);
187  SetupTarget(Z, A, ekin);
188
189  G4double tmax = tkin;
190  if(p == theElectron) tmax *= 0.5;
191  else if(p != thePositron) {
192    G4double ratio = electron_mass_c2/mass;
193    tmax = 2.0*mom2/
194      (electron_mass_c2*(1.0 + ratio*(tkin/mass + 1.0) + ratio*ratio));
195  }
196  G4double t = std::min(cutEnergy, tmax);
197  G4double mom21 = t*(t + 2.0*electron_mass_c2);
198  t = tkin - t;
199  G4double mom22 = t*(t + 2.0*mass);
200  cosTetMaxElec = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
201  if(cosTetMaxElec < cosTetMaxNuc) cosTetMaxElec = cosTetMaxNuc;
202
203  if(cosTetMaxElec < 1.0) {
204    G4double x2 = screenZ/(1.0 - cosTetMaxElec + screenZ);
205    xSection += (x2 - 1.0 - log(x2))/Z;
206  }
207  //  G4cout << "cut= " << ecut << " e= " << tkin << " croosE= "
208  //  << xSection/barn << G4endl;
209
210  if(cosTetMaxNuc < 1.0) {
211    G4double x1 = screenZ*formfactA;
212    G4double x2 = 1.0 - cosTetMaxNuc + screenZ;
213    G4double x3 = 1.0 - x1; 
214    G4double x4 = 1.0/(formfactA*x2 + x3);
215    G4double x5 = screenZ/x2;
216    xSection += ((1.0 - 2.0*x1/x3)*log(x4/x5) - 1.0 + 
217                 x5 - (1.0 - 4.0*x1)*(1.0 - x4))/(x3*x3); 
218  }
219  xSection *= coeff*Z*Z*chargeSquare*invbeta2/mom2; 
220  //  G4cout << " croosE= " << xSection/barn << " screenZ= "
221  //     << screenZ << " formF= " << formfactA << G4endl;
222  return xSection; 
223}
224
225//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
226
227G4double G4MuMscModel::ComputeLambda2(G4double kinEnergy,
228                                      G4double cutEnergy)
229{
230  G4double res = 0.0;
231  SetupParticle(particle);
232  G4double ekin = std::max(keV, kinEnergy);
233
234  const G4Material* mat = currentCouple->GetMaterial();
235  const G4ElementVector* theElementVector = mat->GetElementVector();
236  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
237  size_t nelm = mat->GetNumberOfElements();
238
239  SetupKinematic(ekin);
240
241  G4double tmax = tkin;
242  if(particle == theElectron) tmax *= 0.5;
243  else if(particle != thePositron) {
244    G4double ratio = electron_mass_c2/mass;
245    tmax = 2.0*mom2/
246      (electron_mass_c2*(1.0 + ratio*(tkin/mass + 1.0) + ratio*ratio));
247  }
248  G4double t = std::min(cutEnergy, tmax);
249  G4double mom21 = t*(t + 2.0*electron_mass_c2);
250  t = tkin - t;
251  G4double mom22 = t*(t + 2.0*mass);
252  cosTetMaxElec = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
253  if(cosTetMaxElec < 0.0) cosTetMaxElec = 0.0;
254
255  G4double x, x1, x2, y;
256
257  for (size_t i=0; i<nelm; i++) {
258    const G4Element* elm = (*theElementVector)[i];
259    G4double Z = elm->GetZ();
260    SetupTarget(Z, elm->GetN(), tkin);
261    G4double s = 0.0; 
262    G4double costm = cosTetMaxElec;
263    if(costm < cosTetMaxNuc) costm = cosTetMaxNuc;
264    if(costm < 1.0) {
265      x = 1.0 - costm + screenZ;
266      y = (x - screenZ*(screenZ/x + 2.0*log(x/screenZ)))/Z;
267      if(y < 0.0) {
268        nwarnings++;
269        if(nwarnings < nwarnlimit) 
270          G4cout << "Electron scattering <0 for L2 " << y << G4endl;
271        y = 0.0;
272      }
273      s += y;
274    }
275    //  G4cout << "cut= " << cut << " e= " << tkin << " croosE= "
276    //  << xSection/barn << G4endl;
277
278    // limit main integral because of nuclear size effect
279
280    if(cosTetMaxNuc < 1.0) {
281      x1 = screenZ*formfactA;
282      x2 = 1.0 - cosTetMaxNuc + screenZ;
283      G4double x3 = 1.0 - x1; 
284      G4double f  = 1.0/formfactA;
285      G4double d  = f - screenZ;
286      G4double x4 = f/(x2 + d);
287      G4double x5 = screenZ/x2;
288      y = (screenZ*(1.0 - x5) + (d*d - screenZ*(2.0*d - 3.0*screenZ))*(1.0 - x4)/f -
289           2.0*screenZ*f*log(x4/x5)/d)/(x3*x3); 
290      if(y < 0.0) {
291        nwarnings++;
292        if(nwarnings < nwarnlimit) 
293          G4cout << "Nuclear scattering <0 for L2 " << y << G4endl;
294        y = 0.0;
295      }
296      s += y;
297    }
298
299    res += Z*Z*s*theAtomNumDensityVector[i]; 
300  }
301  res *= 0.25*coeff*chargeSquare*invbeta2/mom2; 
302  //  G4cout << " croosE= " << xSection/barn << " screenZ= "
303  //     << screenZ << " formF= " << formfactA << G4endl;
304  return res; 
305}
306
307//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308
309G4double G4MuMscModel::ComputeTruePathLengthLimit(
310                             const G4Track& track,
311                             G4PhysicsTable* theTable,
312                             G4double currentMinimalStep)
313{
314  G4double tlimit = currentMinimalStep;
315  const G4DynamicParticle* dp = track.GetDynamicParticle();
316
317  // initialisation for 1st step 
318  if(track.GetCurrentStepNumber() == 1) {
319    inside = false;
320    SetupParticle(dp->GetDefinition());
321    theLambdaTable = theTable;
322    if(newrun && buildTables) BuildTables();
323  }
324
325  // initialisation for each step 
326  preKinEnergy = dp->GetKineticEnergy();
327  DefineMaterial(track.GetMaterialCutsCouple());
328  lambda0 = GetLambda(preKinEnergy);
329  currentRange = 
330    theManager->GetRangeFromRestricteDEDX(particle,preKinEnergy,currentCouple);
331
332  // extra check for abnormal situation
333  // this check needed to run MSC with eIoni and eBrem inactivated
334  if(tlimit > currentRange) tlimit = currentRange;
335
336  // stop here if small range particle
337  if(inside) return tlimit;   
338
339  // pre step
340  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
341  G4StepStatus stepStatus = sp->GetStepStatus();
342  G4double presafety = sp->GetSafety();
343
344  // compute presafety again if presafety <= 0 and no boundary
345  // i.e. when it is needed for optimization purposes
346  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) 
347    presafety = safetyHelper->ComputeSafety(sp->GetPosition()); 
348
349  //  G4cout << "G4MuMscModel::ComputeTruePathLengthLimit tlimit= "
350  //     <<tlimit<<" safety= " << presafety
351  //     << " range= " <<currentRange<<G4endl;
352
353  // far from geometry boundary
354  if(currentRange < presafety) {
355    inside = true; 
356   
357    // limit mean scattering angle
358  } else {
359    tlimit = std::min(facrange*lambda0, tlimit);
360  }
361  /*
362  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
363         << " L0= " << lambda0 << " R= " << currentRange
364         << "tlimit= " << tlimit 
365         << " currentMinimalStep= " << currentMinimalStep << G4endl;
366  */
367  return tlimit;
368}
369
370//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
371
372G4double G4MuMscModel::ComputeGeomPathLength(G4double truelength)
373{
374  tPathLength  = truelength;
375  zPathLength  = tPathLength;
376
377  G4double tau = tPathLength/lambda0;
378  lambdaeff    = lambda0;
379  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
380  //     << " lambda0= " << lambda0 << " tau= " << tau << G4endl;
381  // small step
382  if(tau < numlimit) {
383    par1 = -1. ; 
384    par2 = par3 = 0. ; 
385    zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
386
387    // medium step
388  } else if(tPathLength < currentRange*dtrl) {
389    zPathLength = lambda0*(1.0 - exp(-tau));
390
391  } else if(tkin < mass) {
392
393    par1 = 1./currentRange;
394    par2 = 1./(par1*lambda0);
395    par3 = 1.+ par2;
396    lambdaeff = 1.0/(par1*par3);
397    G4double x = tPathLength/currentRange;
398    G4double x1;
399    if(x < numlimit) x1 = x*(1.0  - 0.5*x + x*x/3.0);
400    else             x1 = log(1.0 - x); 
401    zPathLength = lambdaeff*(1.-exp(par3*x1));
402
403  } else {
404
405    G4double T1 = theManager->GetEnergy(particle,
406                                        currentRange-tPathLength,
407                                        currentCouple);
408    G4double lambda1 = GetLambda(T1);
409
410    par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
411    par2 = 1./(par1*lambda0) ;
412    par3 = 1.+ par2 ;
413    lambdaeff = 1.0/(par1*par3);
414    zPathLength = lambdaeff*(1.-exp(par3*log(lambda1/lambda0)));
415  }
416
417  //  if(zPathLength > lambda0) zPathLength = lambda0;
418  if(zPathLength > tPathLength) zPathLength = tPathLength;
419
420  return zPathLength;
421}
422
423//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
424
425G4double G4MuMscModel::ComputeTrueStepLength(G4double geomStepLength)
426{
427  // step defined other than transportation
428  if(geomStepLength == zPathLength) return tPathLength;
429
430  tPathLength  = geomStepLength;
431  zPathLength  = geomStepLength;
432  G4double tau = geomStepLength/lambda0;
433  if(tau < numlimit) {
434    tPathLength *= (1.0 + 0.5*tau - tau*tau/3.0); 
435 
436  } else if(par1 <  0.) {
437    tPathLength = -lambda0*log(1.0 - tau); 
438
439  } else {
440    G4double x = par1*par3*geomStepLength;
441    if(x < numlimit) 
442      tPathLength = (1.- exp(- x*(1.- 0.5*x + x*x/3.0)/par3))/par1 ;
443    else if (x < 1.0) 
444      tPathLength = (1.-exp(log(1.- x)/par3))/par1;
445    else 
446      tPathLength = currentRange;
447  }
448  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
449
450  return tPathLength;
451}
452
453//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
454
455void G4MuMscModel::SampleScattering(const G4DynamicParticle* dynParticle,
456                                    G4double safety)
457{
458  G4double kinEnergy = dynParticle->GetKineticEnergy();
459  if(kinEnergy == 0.0) return;
460  G4double x1  = 0.5*tPathLength/lambdaeff;
461 
462  /*
463  G4cout << "G4MuMscModel::SampleScattering t(mm)= " << tPathLength
464         << " 1/lambdaeff= " << 1.0/lambdaeff
465         << " matIdx= " << currentMaterialIndex << G4endl;
466  */
467  /*
468  G4double y1  = 1.0 - x1;
469  G4double x2  = tPathLength*GetLambda2(0.5*(preKinEnergy + kinEnergy));
470  G4double x3  = (x2 - x1*x1)/(x1*y1);
471  if(x3 <= 0.0 || x3 >= 0.33) {
472    nwarnings++;
473    if(nwarnings < nwarnlimit)
474      G4cout << "G4MuMscModel::SampleScattering: ePre(MeV)= " << preKinEnergy/MeV
475             << " ePost(MeV)= " << kinEnergy/MeV
476             << " <x>= " << x1 << " sqrt(<x^2>)= " << sqrt(x2)
477             << " x3= " << x3
478             << G4endl;
479    x3 = std::min(1.0/y1,0.16666);
480  }
481  G4double x4 = 0.25*(3.0*x3 + sqrt(x3*(x3 + 8.0)))/(1.0 - x3);
482  */
483
484  G4double x  = G4UniformRand();
485  G4double z;
486 
487  //if(x < y1) z = x1*pow(x/y1,x4);
488  //else       z = 1.0 - y1*pow((1.0 - x)/x1,x4);
489 
490  z = -x1*log(x);
491
492  G4double cost = 1.0 - 2.0*z;
493  if(cost < -1.0) cost = -1.0;
494  else if(cost > 1.0) cost = 1.0;
495  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
496
497  G4double phi  = twopi*G4UniformRand();
498
499  G4double dirx = sint*cos(phi);
500  G4double diry = sint*sin(phi);
501
502  //  G4cout << "G4MuMscModel::SampleSecondaries: tstep(mm)= " << truestep/mm
503  //     << " lambdaeff= " << lambdaeff
504  //     << " rms= " << rms << G4endl;
505
506  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
507  G4ThreeVector newDirection(dirx,diry,cost);
508  newDirection.rotateUz(oldDirection);
509  fParticleChange->ProposeMomentumDirection(newDirection);
510
511  if (latDisplasment && safety > tlimitminfix) {
512    G4double rms= sqrt(2.0*x1);
513    G4double rx = zPathLength*(0.5*dirx + invsqrt12*G4RandGauss::shoot(0.0,rms));
514    G4double ry = zPathLength*(0.5*diry + invsqrt12*G4RandGauss::shoot(0.0,rms));
515    G4double r  = sqrt(rx*rx + ry*ry);
516/*
517    G4cout << "G4MuMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
518           << " sinTheta= " << sth << " r(mm)= " << r
519           << " trueStep(mm)= " << truestep
520           << " geomStep(mm)= " << zPathLength
521           << G4endl;
522*/
523   
524    G4ThreeVector latDirection(rx,ry,0.0);
525    latDirection.rotateUz(oldDirection);
526
527    G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
528    G4double fac = 1.;
529    if(r >  safety) {
530      //  ******* so safety is computed at boundary too ************
531      G4double newsafety = safetyHelper->ComputeSafety(Position);
532      if(r > newsafety)
533        fac = newsafety/r ;
534    } 
535
536    if(fac > 0.) {
537      // compute new endpoint of the Step
538      G4ThreeVector newPosition = Position+fac*r*latDirection;
539
540      // definitely not on boundary
541      if(1. == fac) {
542        safetyHelper->ReLocateWithinVolume(newPosition);
543           
544      } else {
545        // check safety after displacement
546        G4double postsafety = safetyHelper->ComputeSafety(newPosition);
547
548        // displacement to boundary
549        if(postsafety <= 0.0) {
550          safetyHelper->Locate(newPosition, newDirection);
551
552          // not on the boundary
553        } else { 
554              safetyHelper->ReLocateWithinVolume(newPosition);
555        }
556      }
557      fParticleChange->ProposePosition(newPosition);
558    } 
559  }
560}
561
562//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
563
564void G4MuMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
565                                     const G4MaterialCutsCouple*,
566                                     const G4DynamicParticle*,
567                                     G4double,
568                                     G4double)
569{}
570
571//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.