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

Timestamp:

Apr 6, 2009, 12:30:29 PM (16 years ago)

Author:

garnier

Message:

update processes

File:

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

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4BGGPionElasticXS.cc,v 1.1 2007/03/13 15:19:30 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4BGGPionElasticXS.cc,v 1.3 2008/12/01 16:50:23 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -46 +46 @@
 #include "G4GlauberGribovCrossSection.hh"
 #include "G4UPiNuclearCrossSection.hh"
+#include "G4HadronNucleonXsc.hh"
 #include "G4PionPlus.hh"
 #include "G4PionMinus.hh"
@@ -52 +53 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4BGGPionElasticXS::G4BGGPionElasticXS(const G4ParticleDefinition* p) 
+G4BGGPionElasticXS::G4BGGPionElasticXS(const G4ParticleDefinition*) 
 {
   verboseLevel = 0;
-  thEnergy     = 100.*GeV;
-  if(p == G4PionPlus::PionPlus() || p == G4PionMinus::PionMinus()) {
-    fPion = new G4UPiNuclearCrossSection();
-    fGlauber = new G4GlauberGribovCrossSection();
-    particle = p;
-    Initialise();
-  } else {
-    fPion = 0;
-    fGlauber = 0;
-    particle = 0;
-    if(p) G4cout << "### G4BGGPionElasticXS WARNING: is not applicable to " 
-                 << p->GetParticleName()
-                 << G4endl;
-    else  G4cout << "### G4BGGPionElasticXS WARNING: particle is not defined " 
-                 << G4endl;
-  }
+  fGlauberEnergy = 91.*GeV;
+  fLowEnergy = 20.*MeV;
+  fPion = 0;
+  fGlauber = 0;
+  fHadron  = 0;
+  particle = 0;
+  isPiplus = false;
+  isInitialized = false;
 }
 
@@ -79 +72 @@
   delete fGlauber;
   delete fPion;
+  delete fHadron;
 }
 
@@ -90 +84 @@
   G4double cross = 0.0;
   G4double ekin = dp->GetKineticEnergy();
-  G4int iz = G4int(Z + 0.5);
+  G4int iz = G4int(Z);
   if(iz > 92) iz = 92;
 
-  if(ekin > thEnergy) {
-    cross = theFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
+  if(ekin <= fLowEnergy) {
+    cross = theCoulombFac[iz];
+    if(isPiplus) { cross *= CoulombFactor(ekin, A); } 
+
+  } else if(iz == 1) {
+    if( A < 1.5) {
+      //fHadron->GetHadronNucleonXscPDG(dp, G4Proton::Proton());
+      //fHadron->GetHadronNucleonXscEL(dp, G4Proton::Proton());
+      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
+      //fHadron->GetHadronNucleonXscVU(dp, G4Proton::Proton());
+      //fHadron->GetHadronNucleonXscMK(dp, G4Proton::Proton());
+      cross = fHadron->GetElasticHadronNucleonXsc();
+    } else {
+      fHadron->GetHadronNucleonXscNS(dp, G4Proton::Proton());
+      cross = fHadron->GetElasticHadronNucleonXsc();
+      fHadron->GetHadronNucleonXscNS(dp, G4Neutron::Neutron());
+      cross += fHadron->GetElasticHadronNucleonXsc();
+    }
+  } else if(ekin > fGlauberEnergy) {
+    cross = theGlauberFac[iz]*fGlauber->GetElasticGlauberGribov(dp, Z, A);
   } else {
     cross = fPion->GetElasticCrossSection(dp, Z, A);
@@ -112 +124 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4BGGPionElasticXS::BuildPhysicsTable(const G4ParticleDefinition&)
-{
+void G4BGGPionElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
+{
+  if(&p == G4PionPlus::PionPlus() || &p == G4PionMinus::PionMinus()) {
+    particle = &p;
+    Initialise();
+  } else {
+    G4cout << "### G4BGGPionElasticXS WARNING: is not applicable to " 
+           << p.GetParticleName()
+           << G4endl;
+  }
 }
 
@@ -127 +147 @@
 void G4BGGPionElasticXS::Initialise() 
 {
+  if(isInitialized) return;
+  isInitialized = true;
+
+  fPion = new G4UPiNuclearCrossSection();
+  fGlauber = new G4GlauberGribovCrossSection();
+  fHadron  = new G4HadronNucleonXsc();
+  fPion->BuildPhysicsTable(*particle);
+  fGlauber->BuildPhysicsTable(*particle);
+  if(particle == G4PionPlus::PionPlus()) isPiplus = true;
+
   G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
   G4ThreeVector mom(0.0,0.0,1.0);
-  G4DynamicParticle dp(part, mom, thEnergy);
+  G4DynamicParticle dp(part, mom, fGlauberEnergy);
 
   G4NistManager* nist = G4NistManager::Instance();
-  G4double A = nist->GetAtomicMassAmu(2);
-
-  G4double csup, csdn;
+
+  G4double A, csup, csdn;
 
   if(verboseLevel > 0) G4cout << "### G4BGGPionElasticXS::Initialise for "
@@ -147 +176 @@
     csdn = fPion->GetElasticCrossSection(&dp, Z, A);
 
-    theFac[iz] = csdn/csup;
+    theGlauberFac[iz] = csdn/csup;
     if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A 
-                                << " factor= " << theFac[iz] << G4endl; 
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-
+                                << " factor= " << theGlauberFac[iz] << G4endl; 
+  }
+  dp.SetKineticEnergy(fLowEnergy);
+  fHadron->GetHadronNucleonXscNS(&dp, G4Proton::Proton());
+  theCoulombFac[1] = fHadron->GetElasticHadronNucleonXsc();
+  if(isPiplus) { theCoulombFac[1] /= CoulombFactor(fLowEnergy,1.0); }
+
+  for(G4int iz=2; iz<93; iz++) {
+
+    G4double Z = G4double(iz);
+    A = nist->GetAtomicMassAmu(iz);
+
+    theCoulombFac[iz] = fPion->GetElasticCrossSection(&dp, Z, A);
+    if(isPiplus) { theCoulombFac[iz] /= CoulombFactor(fLowEnergy,A); }
+    if(verboseLevel > 0) G4cout << "Z= " << Z <<  "  A= " << A 
+                                << " factor= " << theCoulombFac[iz] << G4endl; 
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4BGGPionElasticXS::CoulombFactor(G4double kinEnergy, G4double A)
+{
+  G4double res= 0.0;
+  if(kinEnergy <= DBL_MIN) return res;
+  else if(A < 1.5) return kinEnergy*kinEnergy;
+  
+  G4double elog = std::log10(kinEnergy/GeV);
+
+  // from G4ProtonInelasticCrossSection
+  G4double f1 = 8.0  - 8.0/A - 0.008*A;
+  G4double f2 = 2.34 - 5.4/A - 0.0028*A;
+
+  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))));
+  return res;  
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+