Changeset 1340 for trunk/source/processes/hadronic/cross_sections/src/G4BGGPionElasticXS.cc

Timestamp:

Nov 5, 2010, 3:45:55 PM (14 years ago)

Author:

garnier

Message:

update ti head

File:

• trunk/source/processes/hadronic/cross_sections/src/G4BGGPionElasticXS.cc (modified) (7 diffs)

trunk/source/processes/hadronic/cross_sections/src/G4BGGPionElasticXS.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4BGGPionElasticXS.cc,v 1.7 2009/11/19 19:40:45 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4BGGPionElasticXS.cc,v 1.8 2010/10/12 06:03:37 dennis Exp $
+// GEANT4 tag $Name: hadr-cross-V09-03-12 $
 //
 // -------------------------------------------------------------------
@@ -77 +77 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4BGGPionElasticXS::GetIsoZACrossSection(const G4DynamicParticle* dp, 
-                                                    G4double Z,
-                                                    G4double A, 
-                                                    G4double)
+G4double
+G4BGGPionElasticXS::GetZandACrossSection(const G4DynamicParticle* dp, 
+                                         G4int Z, G4int A, G4double)
 {
   G4double cross = 0.0;
   G4double ekin = dp->GetKineticEnergy();
-  G4int iz = G4int(Z);
-  if(iz > 92) iz = 92;
+//  G4int iz = G4int(Z);
+  if(Z > 92) Z = 92;
 
   if(ekin <= fLowEnergy) {
-    cross = theCoulombFac[iz];
+    cross = theCoulombFac[Z];
     if(isPiplus) { cross *= CoulombFactor(ekin, A); } 
 
-  } else if(iz == 1) {
-    if( A < 1.5) {
+  } else if(Z == 1) {
+    if( A < 2) {
       fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
       cross = fHadron->GetElasticHadronNucleonXsc();
@@ -102 +101 @@
     }
   } else if(ekin > fGlauberEnergy) {
-    cross = theGlauberFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
+    cross = theGlauberFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
   } else {
     cross = fPion->GetElasticCrossSection(dp, Z, A);
@@ -159 +158 @@
   G4NistManager* nist = G4NistManager::Instance();
 
-  G4double A, csup, csdn;
+  G4double csup, csdn;
+  G4int A;
 
   if(verboseLevel > 0) G4cout << "### G4BGGPionElasticXS::Initialise for "
@@ -167 +167 @@
 
     G4double Z = G4double(iz);
-    A = nist->GetAtomicMassAmu(iz);
-
-    csup = fGlauber->GetElasticGlauberGribov(&dp, Z, A);
-    csdn = fPion->GetElasticCrossSection(&dp, Z, A);
+    A = G4lrint(nist->GetAtomicMassAmu(iz));
+
+    csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A);
+    csdn = fPion->GetElasticCrossSection(&dp, iz, A);
 
     theGlauberFac[iz] = csdn/csup;
@@ -179 +179 @@
   fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
   theCoulombFac[1] = fHadron->GetElasticHadronNucleonXsc();
-  if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1.0); }
+  if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1); }
 
   for(G4int iz=2; iz<93; iz++) {
 
     G4double Z = G4double(iz);
-    A = nist->GetAtomicMassAmu(iz);
-
-    theCoulombFac[iz] = fPion->GetElasticCrossSection(&dp, Z, A);
+    A = G4lrint(nist->GetAtomicMassAmu(iz));
+
+    theCoulombFac[iz] = fPion->GetElasticCrossSection(&dp, iz, A);
     if(isPiplus) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy,A); }
     if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A 
@@ -193 +193 @@
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4BGGPionElasticXS::CoulombFactor(G4double kinEnergy, G4double A)
+
+G4double G4BGGPionElasticXS::CoulombFactor(G4double kinEnergy, G4int A)
 {
   G4double res= 0.0;
   if(kinEnergy <= DBL_MIN) return res;
-  else if(A < 1.5) return kinEnergy*kinEnergy;
+  else if(A < 2) return kinEnergy*kinEnergy;
   
   G4double elog = std::log10(kinEnergy/GeV);
+  G4double aa = A;
 
   // from G4ProtonInelasticCrossSection
-  G4double f1 = 8.0  - 8.0/A - 0.008*A;
-  G4double f2 = 2.34 - 5.4/A - 0.0028*A;
+  G4double f1 = 8.0  - 8.0/aa - 0.008*aa;
+  G4double f2 = 2.34 - 5.4/aa - 0.0028*aa;
 
   res = 1.0/(1.0 + std::exp(-f1*(elog + f2)));
  
-  f1 = 5.6 - 0.016*A;
-  f2 = 1.37 + 1.37/A;
-  res *= ( 1.0 + (0.8 + 18./A - 0.002*A)/(1.0 + std::exp(f1*(elog + f2))));
+  f1 = 5.6 - 0.016*aa;
+  f2 = 1.37 + 1.37/aa;
+  res *= ( 1.0 + (0.8 + 18./aa - 0.002*aa)/(1.0 + std::exp(f1*(elog + f2))));
   return res;  
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-