source: trunk/source/processes/hadronic/models/photolepton_hadron/muon_nuclear/src/G4MuNuclearInteraction.cc

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

update ti head
File size: 20.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: G4MuNuclearInteraction.cc,v 1.14 2010/09/08 08:59:29 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-ref-09 $
28//
29// $Id:
30// --------------------------------------------------------------
31//      GEANT 4 class implementation file
32//
33//      History: first implementation, based on object model of
34//      2nd December 1995, G.Cosmo
35//      -------- G4MuNuclearInteraction physics process ---------
36//                by Laszlo Urban, May 1998
37//      added simple model for hadronic vertex, J.P. Wellisch, November 1998
38// --------------------------------------------------------------
39// 26/10/1998: new corr.s from R.Kokoulin + cleanup , L.Urban
40// 23/01/2009  V.Ivanchenko Add deregistration
41//
42
43#include "G4MuNuclearInteraction.hh"
44#include "G4UnitsTable.hh"
45#include "G4HadronicProcessStore.hh"
46
47// static members ........
48G4int G4MuNuclearInteraction::nzdat =  5 ;
49G4double G4MuNuclearInteraction::zdat[]={1.,4.,13.,29.,92.};
50G4double G4MuNuclearInteraction::adat[]={1.01,9.01,26.98,63.55,238.03};
51G4int G4MuNuclearInteraction::ntdat = 8 ;
52G4double G4MuNuclearInteraction::tdat[]={1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,
53                                                             1.e9,1.e10};
54G4int G4MuNuclearInteraction::NBIN = 1000 ; 
55G4double G4MuNuclearInteraction::ya[1001]={0.};
56G4double G4MuNuclearInteraction::proba[5][8][1001]={{{0.}}};
57 
58G4MuNuclearInteraction::G4MuNuclearInteraction(const G4String& processName)
59  : G4VDiscreteProcess(processName), 
60    theMeanFreePathTable(0),
61    theCrossSectionTable(0),
62    LowestKineticEnergy (1.*GeV),
63    HighestKineticEnergy (1000000.*TeV),
64    TotBin(100),
65    CutFixed ( 0.200*GeV),
66    GramPerMole(g/mole),
67    theMuonMinus ( G4MuonMinus::MuonMinus() ),
68    theMuonPlus ( G4MuonPlus::MuonPlus() ),
69    thePionZero (G4PionZero::PionZero() )
70{ 
71  SetProcessSubType(fHadronInelastic);
72  G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
73}
74 
75G4MuNuclearInteraction::~G4MuNuclearInteraction()
76{
77  G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
78
79  if (theMeanFreePathTable) {
80    theMeanFreePathTable->clearAndDestroy();
81    delete theMeanFreePathTable;
82  }
83  if (theCrossSectionTable) {
84    theCrossSectionTable->clearAndDestroy();
85    delete theCrossSectionTable;
86  }
87
88  if (&PartialSumSigma) {
89    PartialSumSigma.clearAndDestroy();
90  }
91}
92 
93void G4MuNuclearInteraction::SetPhysicsTableBining(G4double lowE,
94                                                   G4double highE, G4int nBins)
95{
96  LowestKineticEnergy = lowE; HighestKineticEnergy = highE ; TotBin = nBins ;
97}
98
99
100G4bool G4MuNuclearInteraction::IsApplicable(const G4ParticleDefinition& particle)
101{
102   return(   (&particle == theMuonMinus)||(&particle == theMuonPlus)) ;
103}
104
105void G4MuNuclearInteraction::PreparePhysicsTable(
106                                    const G4ParticleDefinition& aParticleType)
107{
108  G4HadronicProcessStore::Instance()
109    ->RegisterParticleForExtraProcess(this, &aParticleType);
110}
111
112void G4MuNuclearInteraction::BuildPhysicsTable(
113                                    const G4ParticleDefinition& aParticleType)
114{
115  G4HadronicProcessStore::Instance()->PrintInfo(&aParticleType);
116
117  G4double LowEdgeEnergy , Value;
118  G4PhysicsLogVector* ptrVector;
119   
120  if (theCrossSectionTable) {
121    theCrossSectionTable->clearAndDestroy() ;
122    delete theCrossSectionTable ;
123  }
124
125  // make tables for the sampling at initialization
126  if (theMeanFreePathTable == 0) MakeSamplingTables(&aParticleType);
127
128  theCrossSectionTable = new G4PhysicsTable (G4Element::GetNumberOfElements()); 
129  const G4ElementTable* theElementTable = G4Element::GetElementTable() ;
130  G4double AtomicNumber,AtomicWeight ;
131
132  for (size_t J=0; J < G4Element::GetNumberOfElements(); J++ )
133  {
134    ptrVector = new G4PhysicsLogVector(LowestKineticEnergy,
135                                       HighestKineticEnergy,TotBin) ;
136    AtomicNumber = (*theElementTable )[J]->GetZ() ;
137    AtomicWeight = (*theElementTable )[J]->GetA() ;
138
139    for ( G4int i = 0 ; i < TotBin ; i++)
140    {
141      LowEdgeEnergy = ptrVector->GetLowEdgeEnergy(i) ;
142      Value = ComputeMicroscopicCrossSection(&aParticleType,
143                                             LowEdgeEnergy,
144                                             AtomicNumber,AtomicWeight) ;
145      ptrVector->PutValue(i,Value) ;
146    }
147
148    theCrossSectionTable->insertAt( J , ptrVector ) ;
149  }
150
151  G4double FixedEnergy = (LowestKineticEnergy + HighestKineticEnergy)/2. ;
152  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
153  if (theMeanFreePathTable) {
154    theMeanFreePathTable->clearAndDestroy();
155    delete theMeanFreePathTable;
156  }
157  theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
158
159  PartialSumSigma.clearAndDestroy();
160  PartialSumSigma.resize(G4Material::GetNumberOfMaterials());
161
162
163  for (size_t K=0 ; K < G4Material::GetNumberOfMaterials(); K++ ) 
164  { 
165    ptrVector = new G4PhysicsLogVector(LowestKineticEnergy,
166                                       HighestKineticEnergy,
167                                       TotBin ) ;
168
169    const G4Material* material= (*theMaterialTable)[K];
170
171    for ( G4int i = 0 ; i < TotBin ; i++ )     
172    {
173      LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
174      Value = ComputeMeanFreePath( &aParticleType, LowEdgeEnergy,
175                                         material ); 
176       
177      ptrVector->PutValue( i , Value ) ;
178    }
179
180    theMeanFreePathTable->insertAt( K , ptrVector );
181
182    // Compute the PartialSumSigma table at a given fixed energy
183    ComputePartialSumSigma( &aParticleType, FixedEnergy, material);       
184  }
185
186  if (&aParticleType == theMuonPlus) PrintInfoDefinition();
187}
188
189void G4MuNuclearInteraction::ComputePartialSumSigma(
190                                   const G4ParticleDefinition* ParticleType,
191                                   G4double KineticEnergy,
192                                   const G4Material* aMaterial)
193
194// Build the table of cross section per element. The table is built for MATERIALS.
195// This table is used by DoIt to select randomly an element in the material.
196{
197   G4int Imate = aMaterial->GetIndex();
198   G4int NbOfElements = aMaterial->GetNumberOfElements();
199   const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 
200   const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
201
202   PartialSumSigma[Imate] = new G4DataVector();
203
204   G4double SIGMA = 0. ;
205
206   for ( G4int Ielem=0 ; Ielem < NbOfElements ; Ielem++ )
207     {             
208       SIGMA += theAtomNumDensityVector[Ielem] * 
209         ComputeMicroscopicCrossSection( ParticleType, KineticEnergy,
210                                         (*theElementVector)[Ielem]->GetZ(),
211                                         (*theElementVector)[Ielem]->GetA()) ;
212       PartialSumSigma[Imate]->push_back(SIGMA);
213     }
214}
215
216G4double G4MuNuclearInteraction::ComputeMicroscopicCrossSection(
217                                const G4ParticleDefinition* ParticleType,
218                                G4double KineticEnergy,
219                                G4double AtomicNumber,G4double AtomicWeight)
220{
221  static const G4double
222  xgi[] ={ 0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801 };
223  static const G4double
224  wgi[] ={ 0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506 };
225  static const G4double ak1=6.9 ;
226  static const G4double ak2=1.0 ;
227
228  G4double Mass,epmin,epmax,epln,ep,aaa,bbb,hhh,x ;
229  G4int kkk ;
230
231  Mass = ParticleType->GetPDGMass() ;
232
233  G4double CrossSection = 0.0 ;
234
235  if ( AtomicNumber < 1. ) return CrossSection;
236  if ( KineticEnergy <= CutFixed  ) return CrossSection; 
237
238  epmin = CutFixed ;
239  epmax = KineticEnergy + Mass - 0.5*proton_mass_c2 ;
240  if (epmax <= epmin ) return CrossSection; // NaN bug correction
241
242  aaa = std::log(epmin) ;
243  bbb = std::log(epmax) ;
244  kkk = G4int((bbb-aaa)/ak1 +ak2) ;
245  hhh = (bbb-aaa)/kkk ;
246
247  for (G4int l=0 ; l<kkk; l++)
248  {
249    x = aaa + hhh*l ;
250    for (G4int ll=0; ll<8; ll++)
251    {
252      epln=x+xgi[ll]*hhh ;
253      ep = std::exp(epln) ;
254      CrossSection += ep*wgi[ll]* ComputeDMicroscopicCrossSection(ParticleType,
255                                                KineticEnergy, 
256                                                AtomicNumber,AtomicWeight,     
257                                                ep) ;             
258    }
259  }
260  CrossSection *= hhh ;
261
262  if (CrossSection < 0.) CrossSection = 0.;
263
264  return CrossSection;
265}
266
267G4double G4MuNuclearInteraction::ComputeDMicroscopicCrossSection(
268                                 const G4ParticleDefinition* ParticleType,
269                                 G4double KineticEnergy,
270                                 G4double, G4double AtomicWeight,
271                                 G4double epsilon)
272 // Calculates the differential (D) microscopic cross section
273 //   using the cross section formula of R.P. Kokoulin (18/01/98)
274{
275  const G4double alam2 = 0.400*GeV*GeV ;
276  const G4double alam  = 0.632456*GeV ;
277  const G4double coeffn = fine_structure_const/pi ;   
278
279  G4double ep,a,aeff,sigph,v,v1,v2,mass2,up,down ;
280  G4double ParticleMass = ParticleType->GetPDGMass() ;
281  G4double TotalEnergy = KineticEnergy + ParticleMass ;
282
283  G4double DCrossSection = 0. ;
284
285  if((epsilon >= TotalEnergy - 0.5*proton_mass_c2) 
286  ||
287     (epsilon <= CutFixed))
288       return DCrossSection ;
289
290  ep = epsilon/GeV ;
291  a = AtomicWeight/(g/mole) ;
292
293  aeff = 0.22*a+0.78*std::exp(0.89*std::log(a)) ;                   //shadowing
294  sigph = (49.2+11.1*std::log(ep)+151.8/std::sqrt(ep))*microbarn ; //!!!!!!!!!!!
295 
296  v=epsilon/TotalEnergy ;
297  v1 = 1.-v ;
298  v2 = v*v ;
299  mass2 = ParticleMass*ParticleMass ;
300
301  up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1)) ;
302  down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam) ;
303
304  DCrossSection = coeffn*aeff*sigph/epsilon*
305                  (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*std::log(up/down)) ;
306
307  if( DCrossSection < 0.) 
308      DCrossSection = 0. ; 
309
310  return DCrossSection ;
311}
312
313void G4MuNuclearInteraction::MakeSamplingTables(
314                                     const G4ParticleDefinition* ParticleType)
315{
316 
317  G4int nbin;
318  G4double AtomicNumber,AtomicWeight,KineticEnergy,
319           TotalEnergy,Maxep ; 
320
321  G4double ParticleMass = ParticleType->GetPDGMass() ;
322
323  for (G4int iz=0; iz<nzdat; iz++)
324  {
325    AtomicNumber = zdat[iz];
326    AtomicWeight = adat[iz]*GramPerMole ; 
327
328    for (G4int it=0; it<ntdat; it++)
329    {
330      KineticEnergy = tdat[it];
331      TotalEnergy = KineticEnergy + ParticleMass;
332      Maxep = TotalEnergy - 0.5*proton_mass_c2 ;
333
334      G4double CrossSection = 0.0 ;
335
336      G4double c,y,ymin,ymax,dy,yy,dx,x,ep ;
337
338      // calculate the differential cross section
339      // numerical integration in   
340      //  log ...............
341      c = std::log(Maxep/CutFixed) ;
342      ymin = -5. ;
343      ymax = 0. ;
344      dy = (ymax-ymin)/NBIN ; 
345
346      nbin=-1;             
347
348      y = ymin - 0.5*dy ;
349      yy = ymin - dy ;
350      for (G4int i=0 ; i<NBIN; i++)
351      {
352        y += dy ;
353        x = std::exp(y) ;
354        yy += dy ;
355        dx = std::exp(yy+dy)-std::exp(yy) ;
356     
357        ep = CutFixed*std::exp(c*x) ;
358
359        CrossSection += ep*dx*ComputeDMicroscopicCrossSection(ParticleType,
360                                                 KineticEnergy,AtomicNumber,
361                                                 AtomicWeight,ep) ;
362        if(nbin<NBIN)
363        {
364          nbin += 1 ;
365          ya[nbin]=y ;
366          proba[iz][it][nbin] = CrossSection ;
367        }
368      }
369      ya[NBIN]=0. ;
370
371      if(CrossSection > 0.) 
372      {
373        for(G4int ib=0; ib<=nbin; ib++)
374        {
375          proba[iz][it][ib] /= CrossSection ;
376        }
377      }
378    }
379  }
380} 
381
382G4VParticleChange* G4MuNuclearInteraction::PostStepDoIt(
383                                                  const G4Track& trackData,
384                                                  const G4Step& stepData)
385{
386   static const G4double Mass=theMuonPlus->GetPDGMass() ;
387   static const G4double m0=0.2*GeV ;
388
389   aParticleChange.Initialize(trackData);
390   G4Material* aMaterial=trackData.GetMaterial() ;
391
392   const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle(); 
393
394   G4double           KineticEnergy     = aDynamicParticle->GetKineticEnergy();
395   G4ParticleMomentum ParticleDirection = 
396                                      aDynamicParticle->GetMomentumDirection();
397
398   // limits of the energy sampling
399   G4double TotalEnergy = KineticEnergy + Mass ;
400   G4double epmin = CutFixed ;
401   G4double epmax = TotalEnergy - 0.5*proton_mass_c2 ;
402   // check against insufficient energy
403   if (epmax <= epmin )   
404   {
405     aParticleChange.ProposeMomentumDirection( ParticleDirection );
406     aParticleChange.ProposeEnergy( KineticEnergy );
407     aParticleChange.ProposeLocalEnergyDeposit (0.); 
408     aParticleChange.SetNumberOfSecondaries(0);
409     return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
410   }
411
412   // select randomly one element constituing the material 
413   G4Element* anElement = SelectRandomAtom(aMaterial);
414
415   // sample  energy of the secondary ( pi0)
416   //  sampling using tables
417   G4double ep,x,y ;
418   G4int iy ;
419
420   // select sampling table ;
421   G4double lnZ = std::log(anElement->GetZ()) ;
422   G4double delmin = 1.e10 ;
423   G4double del ;
424   G4int izz=0,itt=0,NBINminus1 ;
425   NBINminus1 = NBIN-1 ;
426   for (G4int iz=0; iz<nzdat; iz++)
427   {
428     del = std::abs(lnZ-std::log(zdat[iz])) ;
429     if(del<delmin)
430     {
431        delmin=del ;
432        izz=iz ;
433     }
434   }
435   delmin = 1.e10 ;
436   for (G4int it=0; it<ntdat; it++)
437   {
438     del = std::abs(std::log(KineticEnergy)-std::log(tdat[it])) ;
439     if(del<delmin)
440     {
441       delmin=del;
442       itt=it ;
443     }
444   }
445
446  //sample energy transfer according to the sampling table
447
448   G4double r = G4UniformRand() ;
449 
450   iy = -1 ;
451   do {
452        iy += 1 ;
453      } while (((proba[izz][itt][iy]) < r)&&(iy < NBINminus1)) ;
454
455   //sampling is Done uniformly in y in the bin
456   if( iy < NBIN )
457     y = ya[iy] + G4UniformRand() * ( ya[iy+1] - ya[iy] ) ;
458   else
459     y = ya[iy] ;
460
461   x = std::exp(y) ;
462   ep = epmin*std::exp(x*std::log(epmax/epmin)) ;                             
463
464   // sample scattering angle of mu, but first t should be sampled.
465   G4double yy = ep/TotalEnergy ;
466   G4double tmin=Mass*Mass*yy*yy/(1.-yy) ;
467   G4double tmax=2.*proton_mass_c2*ep ;
468   G4double t1,t2,t,w1,w2,w3,y1,y2,y3,rej ;
469   if(m0<ep)
470   {
471     t1=m0*m0;
472     t2=ep*ep;
473   }
474   else
475   {
476     t1=ep*ep;
477     t2=m0*m0;
478   }
479
480   w1=tmax*t1;
481   w2=tmax+t1 ;
482   w3=tmax*(tmin+t1)/(tmin*w2);
483   y1=1.-yy;
484   y2=0.5*yy*yy;
485   y3=y1+y2;
486
487   // now the sampling of t
488   G4int ntry=0;
489   do
490   {
491     ntry += 1 ;
492     t=w1/(w2*std::exp(G4UniformRand()*std::log(w3))-tmax) ;
493     rej = (1.-t/tmax)*(y1*(1.-tmin/t)+y2)/(y3*(1.-t/t2)); 
494   } while (G4UniformRand() > rej) ;
495
496   // compute angle from t
497   G4double sinth2,theta ; //  sinth2 = std::sin(theta/2)*std::sin(theta/2) !
498   sinth2 = 0.5*(t-tmin)/(2.*(TotalEnergy*(TotalEnergy-ep)-Mass*Mass)-tmin) ;
499   theta = std::acos(1.-2.*sinth2) ;
500   
501   G4double phi=twopi*G4UniformRand() ;
502   G4double sinth=std::sin(theta) ;
503   G4double dirx=sinth*std::cos(phi) , diry=sinth*std::sin(phi) , dirz=std::cos(theta);
504
505   G4ThreeVector finalDirection(dirx,diry,dirz) ;
506   finalDirection.rotateUz(ParticleDirection) ;
507
508   G4double NewKinEnergy = KineticEnergy - ep ;
509   G4double finalMomentum=std::sqrt(NewKinEnergy*
510                       (NewKinEnergy+2.*Mass)) ;
511
512   G4double Ef=NewKinEnergy+Mass ;
513   G4double initMomentum=std::sqrt(KineticEnergy*(TotalEnergy+Mass)) ;
514
515   // G4double Q2=2.*(TotalEnergy*Ef-initMomentum*finalMomentum*std::cos(theta)-Mass*Mass) ;
516
517   aParticleChange.ProposeMomentumDirection( finalDirection );
518   aParticleChange.ProposeEnergy( NewKinEnergy );
519
520   G4LorentzVector primaryMomentum(initMomentum*ParticleDirection, TotalEnergy);
521   G4LorentzVector fsMomentum(finalMomentum*finalDirection, Ef);
522   G4LorentzVector momentumTransfere = primaryMomentum-fsMomentum;
523
524   G4DynamicParticle* aGamma = 
525           new G4DynamicParticle(G4Gamma::Gamma(), momentumTransfere);
526   G4Track gammaTrack(aGamma, trackData.GetGlobalTime(), trackData.GetPosition() );
527   gammaTrack.SetStep(trackData.GetStep());
528   G4Nucleus theTarget(aMaterial);
529
530   G4VParticleChange* aHadronicFS = theHadronicVertex.ApplyYourself(theTarget, gammaTrack);
531   //delete aGamma;
532
533   G4int numSecondaries = aHadronicFS->GetNumberOfSecondaries();
534   aParticleChange.SetNumberOfSecondaries(numSecondaries);
535
536   //   G4ParticleMomentum secondaryMomentum = G4ThreeVector(0.,0.,0.);
537   for(G4int iSec=0; iSec<numSecondaries; iSec++) 
538   {
539     //secondaryMomentum
540     //  = secondaryMomentum + aHadronicFS->GetSecondary(iSec)->GetMomentum();
541     aParticleChange.AddSecondary(aHadronicFS->GetSecondary(iSec));
542   }
543   aHadronicFS->Clear();
544
545   return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
546}
547
548G4Element* G4MuNuclearInteraction::SelectRandomAtom(G4Material* aMaterial) const
549{
550  // select randomly 1 element within the material
551
552  G4int Index = aMaterial->GetIndex();
553  G4int NumberOfElements = aMaterial->GetNumberOfElements();
554  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
555  if(1 == NumberOfElements) return ((*theElementVector)[0]);
556
557  G4double rval = G4UniformRand()*((*PartialSumSigma[Index])[NumberOfElements-1]);
558  for ( G4int i=0; i < NumberOfElements; i++ ) {
559    if (rval <= (*PartialSumSigma[Index])[i]) return ((*theElementVector)[i]);
560  }
561  G4cout << "G4MuNuclearInteraction WARNING !!! no element selected for '"
562         << aMaterial->GetName()
563         << " 1st element returned." << G4endl;
564  return ((*theElementVector)[0]);
565}
566
567void G4MuNuclearInteraction::PrintInfoDefinition()
568{
569  G4String comments = "cross sections from R. Kokoulin \n ";
570           comments += "         Good description up to 1000 PeV.";
571
572  G4cout << G4endl << GetProcessName() << ":  " << comments
573         << "\n    PhysicsTables from " << G4BestUnit(LowestKineticEnergy,
574                                                     "Energy")
575         << " to " << G4BestUnit(HighestKineticEnergy,"Energy")
576         << " in " << TotBin << " bins. \n";
577
578  G4cout << G4endl;
579}
580
581G4double G4MuNuclearInteraction::GetMeanFreePath(const G4Track& trackData,
582                                                 G4double,
583                                                 G4ForceCondition* condition)
584{
585   const G4DynamicParticle* aDynamicParticle;
586   G4Material* aMaterial;
587   G4double MeanFreePath;
588   G4bool isOutRange ;
589 
590   *condition = NotForced ;
591
592   aDynamicParticle = trackData.GetDynamicParticle();
593   aMaterial = trackData.GetMaterial();
594
595   G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
596
597   if (KineticEnergy <  LowestKineticEnergy)
598     MeanFreePath = DBL_MAX ;
599   else {
600     if (KineticEnergy > HighestKineticEnergy) 
601                                 KineticEnergy = 0.99*HighestKineticEnergy ;
602     MeanFreePath = (*theMeanFreePathTable)(aMaterial->GetIndex())->
603                                 GetValue( KineticEnergy, isOutRange );
604   }
605   return MeanFreePath; 
606} 
607
608G4double G4MuNuclearInteraction::ComputeMeanFreePath(
609                                     const G4ParticleDefinition* ParticleType,
610                                     G4double KineticEnergy,
611                                     const G4Material* aMaterial)
612{
613  const G4ElementVector* theElementVector = aMaterial->GetElementVector() ;
614  const G4double* theAtomNumDensityVector = 
615                                      aMaterial->GetAtomicNumDensityVector();
616
617  G4double SIGMA = 0 ;
618
619  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ )
620  {             
621    SIGMA += theAtomNumDensityVector[i] * 
622             ComputeMicroscopicCrossSection( ParticleType, KineticEnergy,
623                                          (*theElementVector)[i]->GetZ(),
624                                          (*theElementVector)[i]->GetA()) ;
625  }       
626
627  return SIGMA<=0.0 ? DBL_MAX : 1./SIGMA ;
628}
629
630
Note: See TracBrowser for help on using the repository browser.