Changeset 1055 for trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc

Timestamp:

May 28, 2009, 4:26:57 PM (15 years ago)

Author:

garnier

Message:

maj sur la beta de geant 4.9.3

File:

• trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc (modified) (9 diffs)

trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4eCoulombScatteringModel.cc,v 1.59 2008/10/22 18:39:29 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eCoulombScatteringModel.cc,v 1.69 2009/05/10 16:09:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -66 +66 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4double G4eCoulombScatteringModel::ScreenRSquare[] = {0.0};
+G4double G4eCoulombScatteringModel::FormFactor[]    = {0.0};
+
 using namespace std;
 
@@ -84 +87 @@
   currentMaterial = 0; 
   currentElement  = 0;
-  a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
+  lowEnergyLimit = keV;
   G4double p0 = electron_mass_c2*classic_electr_radius;
   coeff  = twopi*p0*p0;
-  constn = 6.937e-6/(MeV*MeV);
   tkin = targetZ = mom2 = DBL_MIN;
   elecXSection = nucXSection = 0.0;
-  recoilThreshold = DBL_MAX;
+  recoilThreshold = 100.*keV;
   ecut = DBL_MAX;
   particle = 0;
   currentCouple = 0;
-  for(size_t j=0; j<100; j++) {
-    FF[j] = 0.0;
-  } 
+
+  // Thomas-Fermi screening radii
+  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
+
+  if(0.0 == ScreenRSquare[0]) {
+    G4double a0 = electron_mass_c2/0.88534; 
+    G4double constn = 6.937e-6/(MeV*MeV);
+
+    ScreenRSquare[0] = alpha2*a0*a0;
+    for(G4int j=1; j<100; j++) {
+      G4double x = a0*fNistManager->GetZ13(j);
+      ScreenRSquare[j] = alpha2*x*x;
+      x = fNistManager->GetA27(j); 
+      FormFactor[j] = constn*x*x;
+    } 
+  }
 }
 
@@ -121 +136 @@
   if(!isInitialised) {
     isInitialised = true;
-
-    if(pParticleChange)
-      fParticleChange = 
-        reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-    else
-      fParticleChange = new G4ParticleChangeForGamma();
+    fParticleChange = GetParticleChangeForGamma();
   }
   if(mass < GeV && particle->GetParticleType() != "nucleus") {
@@ -157 +167 @@
       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
       //G4cout << "ctm= " << ctm << G4endl;
-      if(ctm < 1.0) cosTetMaxElec = ctm;
+      if(ctm <  1.0) cosTetMaxElec = ctm;
+      if(ctm < -1.0) cosTetMaxElec = -1.0;
     }
   }
@@ -196 +207 @@
 {
   // This method needs initialisation before be called
-
-  G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2; 
+  G4double fac = coeff*targetZ*chargeSquare*kinFactor;
   elecXSection = 0.0;
   nucXSection  = 0.0;
@@ -216 +226 @@
     G4double s  = screenZ*formfactA;
     G4double z1 = 1.0 - cosTetMaxNuc2 + screenZ;
-    G4double d  = (1.0 - s)/formfactA;
+    G4double s1 = 1.0 - s;
+    G4double d  = s1/formfactA;
     //G4cout <<"x1= "<<x1<<" z1= " <<z1<<" s= "<<s << " d= " <<d <<G4endl;
     if(d < 0.2*x1) {
@@ -226 +237 @@
       G4double x2 = x1 + d;
       G4double z2 = z1 + d;
-      x = (1.0 + 2.0*s)*((cosTetMinNuc - cosTetMaxNuc2)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
-                         2.0*log(z1*x2/(z2*x1))/d);
+      x = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
     }
     nucXSection += fac*targetZ*x;
   }
-
   //G4cout<<" cross(bn)= "<<nucXSection/barn<<" xsElec(bn)= "<<elecXSection/barn
   //    << " Asc= " << screenZ << G4endl; 
@@ -283 +292 @@
     G4int ia = SelectIsotopeNumber(currentElement);
     G4double Trec = z1*mom2/(amu_c2*G4double(ia));
-    G4double th = 
-      std::min(recoilThreshold,
-               targetZ*currentElement->GetIonisation()->GetMeanExcitationEnergy());
-
-    if(Trec > th) {
-      G4int iz = G4int(targetZ);
+
+    if(Trec > recoilThreshold) {
       G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
       Trec = z1*mom2/ion->GetPDGMass();