Changeset 1196 for trunk/source/processes/electromagnetic/highenergy/src/G4AnnihiToMuPair.cc

Timestamp:

Nov 25, 2009, 5:13:58 PM (15 years ago)

Author:

garnier

Message:

update CVS release candidate geant4.9.3.01

File:

• trunk/source/processes/electromagnetic/highenergy/src/G4AnnihiToMuPair.cc (modified) (7 diffs)

trunk/source/processes/electromagnetic/highenergy/src/G4AnnihiToMuPair.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4AnnihiToMuPair.cc,v 1.5 2008/10/16 14:29:48 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4AnnihiToMuPair.cc,v 1.6 2009/11/09 18:24:07 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //         ------------ G4AnnihiToMuPair physics process ------
@@ -65 +65 @@
  HighestEnergyLimit = 1000*TeV;
  
+ CurrentSigma = 0.0;
  CrossSecFactor = 1.;
  SetProcessSubType(6);
@@ -88 +89 @@
 //here no tables, just calling PrintInfoDefinition
 {
-   PrintInfoDefinition();
+  CurrentSigma = 0.0;
+  PrintInfoDefinition();
 }
 
@@ -117 +119 @@
   G4double SigmaEl = Sig0*xi*(1.+xi/2.)*sqrt(1.-xi); // per electron
   CrossSection = SigmaEl*Z;         // number of electrons per atom
-  CrossSection *= CrossSecFactor;   //increase the CrossSection by  (default 1)
   return CrossSection;
 }
@@ -123 +124 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack,
-                                           G4double, G4ForceCondition*)
-
-// returns the positron mean free path in GEANT4 internal units
-
-{
-  const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle();
-  G4double PositronEnergy = aDynamicPositron->GetKineticEnergy()
-                                              +electron_mass_c2;
-  G4Material* aMaterial = aTrack.GetMaterial();
+G4double G4AnnihiToMuPair::CrossSectionPerVolume(G4double PositronEnergy, 
+                                                 const G4Material* aMaterial)
+{
   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
 
-  G4double SIGMA = 0 ;
-
-  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ )
+  G4double SIGMA = 0.0;
+
+  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; ++i )
   {
     G4double AtomicZ = (*theElementVector)[i]->GetZ();
@@ -144 +138 @@
       ComputeCrossSectionPerAtom(PositronEnergy,AtomicZ);
   }
-  return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
+  return SIGMA;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack,
+                                           G4double, G4ForceCondition*)
+
+// returns the positron mean free path in GEANT4 internal units
+
+{
+  const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle();
+  G4double PositronEnergy = aDynamicPositron->GetKineticEnergy()
+                                              +electron_mass_c2;
+  G4Material* aMaterial = aTrack.GetMaterial();
+  CurrentSigma = CrossSectionPerVolume(PositronEnergy, aMaterial);
+
+  // increase the CrossSection by CrossSecFactor (default 1)
+  G4double mfp = DBL_MAX;
+  if(CurrentSigma > DBL_MIN) mfp = 1.0/(CurrentSigma*CrossSecFactor);
+
+  return mfp;
 }
 
@@ -162 +177 @@
   // current Positron energy and direction, return if energy too low
   const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle();
-  G4double Epos = aDynamicPositron->GetKineticEnergy()+Mele;
-
- if (Epos < LowestEnergyLimit)
-  { G4cout 
-        << "error in G4AnnihiToMuPair::PostStepDoIt called with energy below"
-           " threshold Epos= "
-        << Epos << G4endl;       // shoud never happen
-        G4Exception(10);
+  G4double Epos = aDynamicPositron->GetKineticEnergy() + Mele; 
+
+  // test of cross section
+  if(CurrentSigma*G4UniformRand() > 
+     CrossSectionPerVolume(Epos, aTrack.GetMaterial())) 
+    {
+      return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
+    }
+
+  if (Epos < LowestEnergyLimit) {
+     return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
   }
-
-  if (Epos < LowestEnergyLimit)
-     return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
 
   G4ParticleMomentum PositronDirection =