Changeset 1340 for trunk/source/processes/hadronic/cross_sections/src/G4BGGNucleonInelasticXS.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/G4BGGNucleonInelasticXS.cc (modified) (8 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4BGGNucleonInelasticXS.cc,v 1.7 2009/11/19 19:40:45 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4BGGNucleonInelasticXS.cc,v 1.9 2010/10/12 06:16:19 dennis Exp $
+// GEANT4 tag $Name: hadr-cross-V09-03-12 $
 //
 // -------------------------------------------------------------------
@@ -75 +75 @@
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4BGGNucleonInelasticXS::GetIsoZACrossSection(const G4DynamicParticle* dp, 
-                                                       G4double Z,
-                                                       G4double A, 
-                                                       G4double)
+
+G4double
+G4BGGNucleonInelasticXS::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]*CoulombFactor(ekin, A);
-  } else if(iz == 1) {
-    if( A < 1.5) {
+    cross = theCoulombFac[Z]*CoulombFactor(ekin, A);
+  } else if(Z == 1) {
+    if( A < 2) {
       fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
       cross = fHadron->GetInelasticHadronNucleonXsc();
@@ -100 +98 @@
     }
   } else if(ekin > fGlauberEnergy) {
-    cross = theGlauberFac[iz]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
+    cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, A);
   } else {
-    cross = fNucleon->GetIsoZACrossSection(dp, Z, A);
+    cross = fNucleon->GetZandACrossSection(dp, Z, A);
   }
 
@@ -156 +154 @@
 
   G4NistManager* nist = G4NistManager::Instance();
-  G4double A = nist->GetAtomicMassAmu(2);
+  G4int A = G4lrint(nist->GetAtomicMassAmu(2));
 
   G4double csup, csdn;
@@ -166 +164 @@
 
     G4double Z = G4double(iz);
-    A = nist->GetAtomicMassAmu(iz);
-
-    csup = fGlauber->GetInelasticGlauberGribov(&dp, Z, A);
-    csdn = fNucleon->GetIsoZACrossSection(&dp, Z, A);
+    A = G4lrint(nist->GetAtomicMassAmu(iz));
+
+    csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
+//    csdn = fNucleon->GetIsoZACrossSection(&dp, Z, A);
+    csdn = fNucleon->GetZandACrossSection(&dp, iz, A);
 
     theGlauberFac[iz] = csdn/csup;
@@ -178 +177 @@
   fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
   theCoulombFac[1] = 
-    fHadron->GetInelasticHadronNucleonXsc()/CoulombFactor(fLowEnergy,1.0);
+    fHadron->GetInelasticHadronNucleonXsc()/CoulombFactor(fLowEnergy,1);
      
   for(G4int iz=2; iz<93; iz++) {
 
     G4double Z = G4double(iz);
-    A = nist->GetAtomicMassAmu(iz);
+    A = G4lrint(nist->GetAtomicMassAmu(iz));
 
     theCoulombFac[iz] = 
-      fNucleon->GetIsoZACrossSection(&dp, Z, A)/CoulombFactor(fLowEnergy,A);
+      fNucleon->GetZandACrossSection(&dp, iz, A)/CoulombFactor(fLowEnergy,A);
 
     if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A 
@@ -193 +192 @@
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4double A)
+
+G4double G4BGGNucleonInelasticXS::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
   if(isProton) {
-    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))));
   } else {
 
-    G4double p3 = 0.6 + 13./A - 0.0005*A;
-    G4double p4 = 7.2449 - 0.018242*A;
-    G4double p5 = 1.36 + 1.8/A + 0.0005*A;
-    G4double p6 = 1. + 200./A + 0.02*A;
-    G4double p7 = 3.0 - (A-70.)*(A-200.)/11000.;
+    G4double p3 = 0.6 + 13./aa - 0.0005*aa;
+    G4double p4 = 7.2449 - 0.018242*aa;
+    G4double p5 = 1.36 + 1.8/aa + 0.0005*aa;
+    G4double p6 = 1. + 200./aa + 0.02*aa;
+    G4double p7 = 3.0 - (aa-70.)*(aa-200.)/11000.;
 
     res = (1.+p3/(1. + std::exp(p4*(elog+p5))))/(1.+std::exp(-p6*(elog+p7)));
@@ -227 +226 @@
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-