From ac2294504c392577c20aa5a04aa75bdf30d4f682 Mon Sep 17 00:00:00 2001
From: Marco Clemencic <marco.clemencic@cern.ch>
Date: Wed, 9 Sep 2015 17:36:09 +0200
Subject: [PATCH] v96r4p1
---
.../G4EmStandardPhysics_option1NoApplyCuts.hh | 5 +-
Geant4/G4LHCblists/cmt/requirements | 2 +-
Geant4/G4LHCblists/doc/release.notes | 13 +-
.../G4EmStandardPhysics_option1NoApplyCuts.cc | 182 +-
Geant4/G4config/cmt/copy_examples_source.csh | 6 +
Geant4/G4config/cmt/copy_source.csh | 2 +-
Geant4/G4config/cmt/requirements | 4 +-
Geant4/G4config/doc/release.notes | 7 +
Geant4/G4examples/cmt/requirements | 2 +-
Geant4/G4examples/doc/release.notes | 3 +
.../G4TestEm3/Test/opt1NoApplyCuts/Plot.C | 58 +-
.../Test/opt1NoApplyCuts/batchrun.sh | 8 +-
.../Test/opt1NoApplyCuts/opt1noapplycuts.mac | 83 +-
.../G4TestEm3/Test/opt1NoApplyCuts/run.sh | 6 +-
.../G4TestEm3/Test/opt2/opt2.mac | 2 +-
.../G4TestEm3/cmt/copyPatchedSource.py | 95 +
.../G4TestEm3/doc/release.notes | 21 +-
.../G4TestEm3/scripts/OUTPUT_3.000 | Bin 166 -> 166 bytes
.../G4TestEm3/srcnew/HistoManager.cc | 107 ++
.../G4TestEm3/srcnew/PhysicsList.cc | 357 ++++
.../G4TestEm3/srcnew/RunAction.cc | 422 ++++
.../G4TestEm3/srcnew/SteppingAction.cc | 174 ++
Geant4/G4processes/cmt/copyPatchedSource.py | 11 +-
Geant4/G4processes/cmt/requirements | 2 +-
Geant4/G4processes/doc/release.notes | 6 +
.../G4ComponentGGHadronNucleusXsc.cc | 1701 +++++++++++++++++
.../G4ComponentGGHadronNucleusXsc.hh | 278 +++
.../cross_sections/G4HadronNucleonXsc.cc | 1676 ++++++++++++++++
.../cross_sections/G4HadronNucleonXsc.hh | 162 ++
Geant4Sys/cmt/requirements | 10 +-
Geant4Sys/doc/release.notes | 15 +
cmt/project.cmt | 2 +-
toolchain.cmake | 2 +-
33 files changed, 5335 insertions(+), 89 deletions(-)
create mode 100755 Geant4/G4examples/extended/electromagnetic/G4TestEm3/cmt/copyPatchedSource.py
create mode 100644 Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/HistoManager.cc
create mode 100644 Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/PhysicsList.cc
create mode 100644 Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/RunAction.cc
create mode 100644 Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/SteppingAction.cc
create mode 100644 Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.cc
create mode 100644 Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.hh
create mode 100644 Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.cc
create mode 100644 Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.hh
diff --git a/Geant4/G4LHCblists/G4LHCblists/G4EmStandardPhysics_option1NoApplyCuts.hh b/Geant4/G4LHCblists/G4LHCblists/G4EmStandardPhysics_option1NoApplyCuts.hh
index 113c45b24..5ffc74344 100644
--- a/Geant4/G4LHCblists/G4LHCblists/G4EmStandardPhysics_option1NoApplyCuts.hh
+++ b/Geant4/G4LHCblists/G4LHCblists/G4EmStandardPhysics_option1NoApplyCuts.hh
@@ -23,8 +23,7 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
-// $Id: G4EmStandardPhysics_option1NoApplyCuts.hh,v 1.2 2010-06-02 17:21:29 vnivanch Exp $
-// GEANT4 tag $Name: not supported by cvs2svn $
+// $Id$
//
//---------------------------------------------------------------------------
//
@@ -36,8 +35,6 @@
// 05.12.2005 V.Ivanchenko add controlled verbosity
// 13.11.2006 V.Ivanchenko set default msc step limit false
// 15.05.2007 V.Ivanchenko rename to _option1
-// 12.07.2011 G.Corti rename as _option1NoApplyCuts to get back to G49.2
-// functionality
//
//----------------------------------------------------------------------------
//
diff --git a/Geant4/G4LHCblists/cmt/requirements b/Geant4/G4LHCblists/cmt/requirements
index 613039459..d55a05503 100755
--- a/Geant4/G4LHCblists/cmt/requirements
+++ b/Geant4/G4LHCblists/cmt/requirements
@@ -3,7 +3,7 @@
# Maintainer : Gloria CORTI
#============================================================================
package G4LHCblists
-version v3r0
+version v3r2
# Structure, i.e. directories to process.
#============================================================================
diff --git a/Geant4/G4LHCblists/doc/release.notes b/Geant4/G4LHCblists/doc/release.notes
index d7820fb03..9eea2a831 100755
--- a/Geant4/G4LHCblists/doc/release.notes
+++ b/Geant4/G4LHCblists/doc/release.notes
@@ -2,7 +2,18 @@
! Package : Geant4/G4LHCblists
! Responsible : Gloria CORTI
! Purpose : Private LHCb physics lists
-!-----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
+
+!======================== G4LHCblists v3r2 2015-04-30 ========================
+! 2015-04-30 - Timothy Williams
+ - Replaced G4EmStandardPhysics_option1NoApplyCuts.cc and .hh with v9.6.p04
+of G4EmStandardPhysics_option1.cc and .hh but with the EmProcessOption of "SetApplyCuts"
+removed. Previous version was v9.4 of G4EmStandardPhysics_option1.cc with "SetApplyCuts" removed.
+
+!======================== G4LHCblists v3r1 2015-03-07 ========================
+! 2015-03-07 - Timothy Williams
+ - Modified G4EmStandardPhysics_option1NoApplyCuts.cc to use UrbanMsc93 and
+ WentzelVI models for multiple scattering as suggested by V.Ivanchenko.
!======================== G4LHCblists v3r0 2013-09-16 ========================
! 2013-09-16 - Gloria Corti
diff --git a/Geant4/G4LHCblists/src/G4EmStandardPhysics_option1NoApplyCuts.cc b/Geant4/G4LHCblists/src/G4EmStandardPhysics_option1NoApplyCuts.cc
index 01c7250bb..60de452d2 100644
--- a/Geant4/G4LHCblists/src/G4EmStandardPhysics_option1NoApplyCuts.cc
+++ b/Geant4/G4LHCblists/src/G4EmStandardPhysics_option1NoApplyCuts.cc
@@ -23,8 +23,7 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
-// $Id: G4EmStandardPhysics_option1NoApplyCuts.cc,v 1.22 2010-12-19 18:11:05 vnivanch Exp $
-// GEANT4 tag $Name: not supported by cvs2svn $
+// $Id$
//
//---------------------------------------------------------------------------
//
@@ -42,15 +41,14 @@
// 15.05.2007 V.Ivanchenko rename to _option1
// 21.04.2008 V.Ivanchenko add long-lived D and B mesons
// 19.11.2010 V.Ivanchenko added WentzelVI model for muons;
-// 12.07.2011 G.Corti disable ApplyCuts option to reproduce 9.2 and rename
-// _option1NoApplyCuts
+// disable ApplyCut option for compatibility with 9.2
//
//----------------------------------------------------------------------------
//
#include "G4EmStandardPhysics_option1NoApplyCuts.hh"
+#include "G4SystemOfUnits.hh"
#include "G4ParticleDefinition.hh"
-#include "G4ProcessManager.hh"
#include "G4LossTableManager.hh"
#include "G4EmProcessOptions.hh"
@@ -61,12 +59,16 @@
#include "G4eMultipleScattering.hh"
#include "G4MuMultipleScattering.hh"
#include "G4hMultipleScattering.hh"
+#include "G4eCoulombScatteringModel.hh"
#include "G4CoulombScattering.hh"
#include "G4WentzelVIModel.hh"
+#include "G4UrbanMscModel93.hh"
+#include "G4UrbanMscModel95.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
+#include "G4UAtomicDeexcitation.hh"
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
@@ -74,6 +76,11 @@
#include "G4hBremsstrahlung.hh"
#include "G4hPairProduction.hh"
+#include "G4MuBremsstrahlungModel.hh"
+#include "G4MuPairProductionModel.hh"
+#include "G4hBremsstrahlungModel.hh"
+#include "G4hPairProductionModel.hh"
+
#include "G4hIonisation.hh"
#include "G4ionIonisation.hh"
#include "G4alphaIonisation.hh"
@@ -95,12 +102,22 @@
#include "G4Alpha.hh"
#include "G4GenericIon.hh"
+#include "G4PhysicsListHelper.hh"
+#include "G4BuilderType.hh"
+
+// factory
+#include "G4PhysicsConstructorFactory.hh"
+//
+G4_DECLARE_PHYSCONSTR_FACTORY(G4EmStandardPhysics_option1NoApplyCuts);
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4EmStandardPhysics_option1NoApplyCuts::G4EmStandardPhysics_option1NoApplyCuts(G4int ver)
: G4VPhysicsConstructor("G4EmStandard_opt1"), verbose(ver)
{
G4LossTableManager::Instance();
+ SetPhysicsType(bElectromagnetic);
+ // std::cout<<"This is the v9.6 option 1 physics list with ApplyCuts option commented out"<<std::endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -109,6 +126,7 @@ G4EmStandardPhysics_option1NoApplyCuts::G4EmStandardPhysics_option1NoApplyCuts(G
: G4VPhysicsConstructor("G4EmStandard_opt1"), verbose(ver)
{
G4LossTableManager::Instance();
+ SetPhysicsType(bElectromagnetic);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -151,12 +169,36 @@ void G4EmStandardPhysics_option1NoApplyCuts::ConstructParticle()
void G4EmStandardPhysics_option1NoApplyCuts::ConstructProcess()
{
- // Add standard EM Processes
+ G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+
+ // muon & hadron bremsstrahlung and pair production
+ G4MuBremsstrahlung* mub = new G4MuBremsstrahlung();
+ G4MuPairProduction* mup = new G4MuPairProduction();
+ G4hBremsstrahlung* pib = new G4hBremsstrahlung();
+ G4hPairProduction* pip = new G4hPairProduction();
+ G4hBremsstrahlung* kb = new G4hBremsstrahlung();
+ G4hPairProduction* kp = new G4hPairProduction();
+ G4hBremsstrahlung* pb = new G4hBremsstrahlung();
+ G4hPairProduction* pp = new G4hPairProduction();
+
+ // muon & hadron multiple scattering
+ G4MuMultipleScattering* mumsc = new G4MuMultipleScattering();
+ mumsc->AddEmModel(0, new G4WentzelVIModel());
+ G4MuMultipleScattering* pimsc = new G4MuMultipleScattering();
+ pimsc->AddEmModel(0, new G4WentzelVIModel());
+ G4MuMultipleScattering* kmsc = new G4MuMultipleScattering();
+ kmsc->AddEmModel(0, new G4WentzelVIModel());
+ G4MuMultipleScattering* pmsc = new G4MuMultipleScattering();
+ pmsc->AddEmModel(0, new G4WentzelVIModel());
+ G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc");
+
+ // high energy limit for e+- scattering models and bremsstrahlung
+ G4double highEnergyLimit = 100*MeV;
+ // Add standard EM Processes
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
- G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if(verbose > 1)
G4cout << "### " << GetPhysicsName() << " instantiates for "
@@ -164,63 +206,111 @@ void G4EmStandardPhysics_option1NoApplyCuts::ConstructProcess()
if (particleName == "gamma") {
- pmanager->AddDiscreteProcess(new G4PhotoElectricEffect);
- pmanager->AddDiscreteProcess(new G4ComptonScattering);
- pmanager->AddDiscreteProcess(new G4GammaConversion);
+ ph->RegisterProcess(new G4PhotoElectricEffect(), particle);
+ ph->RegisterProcess(new G4ComptonScattering(), particle);
+ ph->RegisterProcess(new G4GammaConversion(), particle);
} else if (particleName == "e-") {
G4eIonisation* eioni = new G4eIonisation();
eioni->SetStepFunction(0.8, 1.0*mm);
+
G4eMultipleScattering* msc = new G4eMultipleScattering;
msc->SetStepLimitType(fMinimal);
- pmanager->AddProcess(msc, -1, 1, 1);
- pmanager->AddProcess(eioni, -1, 2, 2);
- pmanager->AddProcess(new G4eBremsstrahlung, -1,-3, 3);
+ G4UrbanMscModel93* msc1 = new G4UrbanMscModel93();
+ G4WentzelVIModel* msc2 = new G4WentzelVIModel();
+ msc1->SetHighEnergyLimit(highEnergyLimit);
+ msc2->SetLowEnergyLimit(highEnergyLimit);
+ msc->AddEmModel(0, msc1);
+ msc->AddEmModel(0, msc2);
+
+ G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
+ G4CoulombScattering* ss = new G4CoulombScattering();
+ ss->SetEmModel(ssm, 1);
+ ss->SetMinKinEnergy(highEnergyLimit);
+ ssm->SetLowEnergyLimit(highEnergyLimit);
+ ssm->SetActivationLowEnergyLimit(highEnergyLimit);
+
+ ph->RegisterProcess(msc, particle);
+ ph->RegisterProcess(eioni, particle);
+ ph->RegisterProcess(new G4eBremsstrahlung(), particle);
+ ph->RegisterProcess(ss, particle);
} else if (particleName == "e+") {
G4eIonisation* eioni = new G4eIonisation();
eioni->SetStepFunction(0.8, 1.0*mm);
+
G4eMultipleScattering* msc = new G4eMultipleScattering;
msc->SetStepLimitType(fMinimal);
- pmanager->AddProcess(msc, -1, 1, 1);
- pmanager->AddProcess(eioni, -1, 2, 2);
- pmanager->AddProcess(new G4eBremsstrahlung, -1,-3, 3);
- pmanager->AddProcess(new G4eplusAnnihilation, 0,-1, 4);
+ G4UrbanMscModel93* msc1 = new G4UrbanMscModel93();
+ G4WentzelVIModel* msc2 = new G4WentzelVIModel();
+ msc1->SetHighEnergyLimit(highEnergyLimit);
+ msc2->SetLowEnergyLimit(highEnergyLimit);
+ msc->AddEmModel(0, msc1);
+ msc->AddEmModel(0, msc2);
+
+ G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
+ G4CoulombScattering* ss = new G4CoulombScattering();
+ ss->SetEmModel(ssm, 1);
+ ss->SetMinKinEnergy(highEnergyLimit);
+ ssm->SetLowEnergyLimit(highEnergyLimit);
+ ssm->SetActivationLowEnergyLimit(highEnergyLimit);
+
+ ph->RegisterProcess(msc, particle);
+ ph->RegisterProcess(eioni, particle);
+ ph->RegisterProcess(new G4eBremsstrahlung(), particle);
+ ph->RegisterProcess(new G4eplusAnnihilation(), particle);
+ ph->RegisterProcess(ss, particle);
} else if (particleName == "mu+" ||
- particleName == "mu-" ) {
+ particleName == "mu-" ) {
- G4MuMultipleScattering* msc = new G4MuMultipleScattering();
- msc->AddEmModel(0, new G4WentzelVIModel());
- pmanager->AddProcess(msc, -1, 1, 1);
- pmanager->AddProcess(new G4MuIonisation, -1, 2, 2);
- pmanager->AddProcess(new G4MuBremsstrahlung, -1,-3, 3);
- pmanager->AddProcess(new G4MuPairProduction, -1,-4, 4);
- pmanager->AddDiscreteProcess(new G4CoulombScattering());
+ ph->RegisterProcess(mumsc, particle);
+ ph->RegisterProcess(new G4MuIonisation(), particle);
+ ph->RegisterProcess(mub, particle);
+ ph->RegisterProcess(mup, particle);
+ ph->RegisterProcess(new G4CoulombScattering(), particle);
- } else if (particleName == "alpha" ||
- particleName == "He3") {
+ } else if (particleName == "alpha" ||
+ particleName == "He3" ) {
- pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
- pmanager->AddProcess(new G4ionIonisation, -1, 2, 2);
+ //ph->RegisterProcess(hmsc, particle);
+ ph->RegisterProcess(new G4hMultipleScattering(), particle);
+ ph->RegisterProcess(new G4ionIonisation(), particle);
} else if (particleName == "GenericIon") {
- pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
- pmanager->AddProcess(new G4ionIonisation, -1, 2, 2);
+ ph->RegisterProcess(hmsc, particle);
+ ph->RegisterProcess(new G4ionIonisation(), particle);
+
+ } else if (particleName == "pi+" ||
+ particleName == "pi-" ) {
+
+ //G4hMultipleScattering* pimsc = new G4hMultipleScattering();
+ ph->RegisterProcess(pimsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
+ ph->RegisterProcess(pib, particle);
+ ph->RegisterProcess(pip, particle);
- } else if (particleName == "pi+" ||
- particleName == "pi-" ||
- particleName == "kaon+" ||
- particleName == "kaon-" ||
- particleName == "proton" ) {
+ } else if (particleName == "kaon+" ||
+ particleName == "kaon-" ) {
- pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
- pmanager->AddProcess(new G4hIonisation, -1, 2, 2);
- pmanager->AddProcess(new G4hBremsstrahlung, -1,-3, 3);
- pmanager->AddProcess(new G4hPairProduction, -1,-4, 4);
+ //G4hMultipleScattering* kmsc = new G4hMultipleScattering();
+ ph->RegisterProcess(kmsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
+ ph->RegisterProcess(kb, particle);
+ ph->RegisterProcess(kp, particle);
+
+ // } else if (particleName == "proton" ) {
+ } else if (particleName == "proton" ||
+ particleName == "anti_proton") {
+
+ //G4hMultipleScattering* pmsc = new G4hMultipleScattering();
+ ph->RegisterProcess(pmsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
+ ph->RegisterProcess(pb, particle);
+ ph->RegisterProcess(pp, particle);
} else if (particleName == "B+" ||
particleName == "B-" ||
@@ -233,7 +323,6 @@ void G4EmStandardPhysics_option1NoApplyCuts::ConstructProcess()
particleName == "anti_deuteron" ||
particleName == "anti_lambda_c+" ||
particleName == "anti_omega-" ||
- particleName == "anti_proton" ||
particleName == "anti_sigma_c+" ||
particleName == "anti_sigma_c++" ||
particleName == "anti_sigma+" ||
@@ -254,14 +343,19 @@ void G4EmStandardPhysics_option1NoApplyCuts::ConstructProcess()
particleName == "xi_c+" ||
particleName == "xi-" ) {
- pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
- pmanager->AddProcess(new G4hIonisation, -1, 2, 2);
+ ph->RegisterProcess(hmsc, particle);
+ ph->RegisterProcess(new G4hIonisation(), particle);
}
}
G4EmProcessOptions opt;
opt.SetVerbose(verbose);
- opt.SetPolarAngleLimit(0.2);
+ opt.SetPolarAngleLimit(CLHEP::pi);
//opt.SetApplyCuts(true);
+
+ // Deexcitation
+ //
+ G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
+ G4LossTableManager::Instance()->SetAtomDeexcitation(de);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/Geant4/G4config/cmt/copy_examples_source.csh b/Geant4/G4config/cmt/copy_examples_source.csh
index 2e5999ae7..015f7144f 100755
--- a/Geant4/G4config/cmt/copy_examples_source.csh
+++ b/Geant4/G4config/cmt/copy_examples_source.csh
@@ -47,3 +47,9 @@ else
endif
cd $here
+
+if( ${pack} == 'TestEm3' ) then
+ ./copyPatchedSource.py
+echo ' executed copyPatchedSource.py '
+endif
+
diff --git a/Geant4/G4config/cmt/copy_source.csh b/Geant4/G4config/cmt/copy_source.csh
index b10fa70c2..2fa4a1725 100755
--- a/Geant4/G4config/cmt/copy_source.csh
+++ b/Geant4/G4config/cmt/copy_source.csh
@@ -55,7 +55,7 @@ cd $here
# NKW 22/02/2014
#Not pretty but will work.
-if ( $pack == "processes" ) then
+if ( $pack == "processes" || $pack == "physics_lists" ) then
echo "Copying patched sources"
./copyPatchedSource.py
endif
diff --git a/Geant4/G4config/cmt/requirements b/Geant4/G4config/cmt/requirements
index dd1e07487..86e2511f3 100755
--- a/Geant4/G4config/cmt/requirements
+++ b/Geant4/G4config/cmt/requirements
@@ -1,5 +1,5 @@
package G4config
-version v96r1p0
+version v96r3p0
branches cmt doc
@@ -19,7 +19,7 @@ use Geant4Files v96r* -no_auto_imports
#==============================================================================
set G4_native_version "9.6.p04" \
override-geant4-version "${G4_NATIVE_VERSION}"
-set G4VERS v96r1p0
+set G4VERS v96r3p0
# =============================================================================
# set Geant4 environment variables
diff --git a/Geant4/G4config/doc/release.notes b/Geant4/G4config/doc/release.notes
index a61926b7c..74e139ef2 100755
--- a/Geant4/G4config/doc/release.notes
+++ b/Geant4/G4config/doc/release.notes
@@ -3,6 +3,13 @@
! Responsible : Gloria Corti/ Nigel Watson
! Purpose : Configuration package for Geant4 build
!-----------------------------------------------------------------------------
+! ======================= G4config v96r3p0 2015-05-11 ========================
+! 2015-05-11 - Timothy Williams
+ - Added execution of srcnew copying to copy_examples_source.csh
+
+! ======================= G4config v96r2p0 2015-04-24 ========================
+! 2015-04-24 - Nige Watson
+ - Updated copy scripts to copy new source for G4physics_list
! ======================= G4config v96r1p0 2015-03-06 ========================
! 2015-03-06 - Nige Watson
diff --git a/Geant4/G4examples/cmt/requirements b/Geant4/G4examples/cmt/requirements
index aade440e3..369912b08 100644
--- a/Geant4/G4examples/cmt/requirements
+++ b/Geant4/G4examples/cmt/requirements
@@ -3,7 +3,7 @@
# Maintainer : James Mccarthy
#============================================================================
package G4examples
-version v5r0
+version v5r1
#============================================================================
# Structure, i.e. directories to process.
diff --git a/Geant4/G4examples/doc/release.notes b/Geant4/G4examples/doc/release.notes
index 05ca94a89..8705bc139 100644
--- a/Geant4/G4examples/doc/release.notes
+++ b/Geant4/G4examples/doc/release.notes
@@ -3,6 +3,9 @@
! Responsible : James Mccarthy
! Purpose : Examples packages provided my Geant4, used as standalone tests
!-----------------------------------------------------------------------------
+!===================== Geant4/G4examples v5r1 2015-05-11 =====================
+! 2015-05-11- Tim Williams
+ -Added copy patchedsource script and srcnew to TestEm3 to copy private source files into place upon build.
!===================== Geant4/G4examples v5r0 2015-03-06 =====================
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/Plot.C b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/Plot.C
index e1c816c18..768bace4e 100644
--- a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/Plot.C
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/Plot.C
@@ -23,11 +23,23 @@ struct M{
int rebin;
};
+struct SeperateProfiles
+{
+ TH1D* PhotonsLead;
+ TH1D* ElectronsLead;
+ TH1D* PositronsLead;
+ TH1D* PhotonsScintillator;
+ TH1D* ElectronsScintillator;
+ TH1D* PositronsScintillator;
+};
+vector<SeperateProfiles> Profiles(15);
vector<M> Control(15);
double nEvents=10000;
int main()
{
+
+
TCanvas *c[15];
ofstream outt;
outt.open ("results.txt");
@@ -109,7 +121,33 @@ int main()
hnew[i]= (TH1D*)h2[i]->Clone("hnew");
hnew[i]->Rebin(Control[i].rebin);
cout<<"nbins after cloning= "<<hnew[i]->GetNbinsX()<<endl;
-
+
+ Profiles[i].PhotonsLead=(TH1D*)In->Get("23");
+ Profiles[i].ElectronsLead=(TH1D*)In->Get("24");
+ Profiles[i].PositronsLead=(TH1D*)In->Get("25");
+ Profiles[i].PhotonsScintillator=(TH1D*)In->Get("26");
+ Profiles[i].ElectronsScintillator=(TH1D*)In->Get("27");
+ Profiles[i].PositronsScintillator=(TH1D*)In->Get("28");
+ Char_t Phl[50];
+ Char_t El[50];
+ Char_t Pol[50];
+ Char_t Phs[50];
+ Char_t Es[50];
+ Char_t Pos[50];
+ sprintf(Phl,"v9.6.p04 Photons in Lead Shower Profile at %.2f GeV",Control[i].Energy);
+ sprintf(El,"v9.6.p04 Electrons in Lead Shower Profile at %.2f GeV",Control[i].Energy);
+ sprintf(Pol,"v9.6.p04 Positrons in Lead Shower Profile at %.2f GeV",Control[i].Energy);
+ sprintf(Phs,"v9.6.p04 Photons in Scintillator Shower Profile at %.2f GeV",Control[i].Energy);
+ sprintf(Es,"v9.6.p04 Electrons in Scintillator Shower Profile at %.2f GeV",Control[i].Energy);
+ sprintf(Pos,"v9.6.p04 Positrons in Scintillator Shower Profile at %.2f GeV;",Control[i].Energy);
+
+ Profiles[i].PhotonsLead->SetTitle(Phl);
+ Profiles[i].ElectronsLead->SetTitle(El);
+ Profiles[i].PositronsLead->SetTitle(Pol);
+ Profiles[i].PhotonsScintillator->SetTitle(Phs);
+ Profiles[i].ElectronsScintillator->SetTitle(Es);
+ Profiles[i].PositronsScintillator->SetTitle(Pos);
+
//set titles of histograms
Char_t titleActive[50];
Char_t titlePassive[50];
@@ -145,6 +183,7 @@ int main()
}
+
double SigmaOverE[15]; //sigma over E
double Xaxis[15]; //1 / sqrt(E)
double SigmaOverEError[15];
@@ -171,6 +210,17 @@ CovElement[k])/AMean[k]));
TFile* out= new TFile("Save.root","RECREATE");
TCanvas*CStraightRes = new TCanvas("CStraightRes","CStraightRes",1800,1000);
+ for(int l=1;l<14;++l)
+ {
+ Profiles[l].PhotonsLead->Write();
+ Profiles[l].ElectronsLead->Write();
+ Profiles[l].PositronsLead->Write();
+ Profiles[l].PhotonsScintillator->Write();
+ Profiles[l].ElectronsScintillator->Write();
+ Profiles[l].PositronsScintillator->Write();
+ }
+
+
//------------1/sqrt(E)------------------------------------------------------------------------
TF1* fitstraight= new TF1("fitstraight","sqrt(([0]*x)**2+(([1])**2))",0,2);
fitstraight->SetParName(0,"a");
@@ -260,10 +310,7 @@ CovElement[k])/AMean[k]));
double LeadSFError[15];
double ScintSFError[15];
double Energ[15];
- double EnergyEr[15];
-
-
- for(int n=1;n<14;++n){
+ double EnergyEr[15];for(int n=1;n<14;++n){
int k=n-1;
Double_t xmin=l2[n]->GetXaxis()->GetXmin();
Double_t xmax=l2[n]->GetXaxis()->GetXmax();
@@ -369,6 +416,7 @@ CovElement[k])/AMean[k]));
// hnew[j]->Write();
}
+
for(int E=0;E<13;++E)
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/batchrun.sh b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/batchrun.sh
index 572b272fd..931fcaaa3 100755
--- a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/batchrun.sh
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/batchrun.sh
@@ -1,13 +1,11 @@
source /afs/cern.ch/lhcb/software/releases/LBSCRIPTS/prod/InstallArea/scripts/LbLogin.sh
-SetupProject ROOT
+SetupProject Geant4 HEAD --nightly lhcb-gauss-dev --user-area /afs/cern.ch/work/t/tiwillia/private/9.6.p04cmtuser
-SetupProject Geant4 HEAD --nightly lhcb-gauss-dev
+export TESTAREAROOT="/afs/cern.ch/work/t/tiwillia/private/9.6.p04cmtuser/Geant4_HEAD/GEANT4/GEANT4_HEAD/Geant4/G4examples/extended/electromagnetic/G4TestEm3"
-export TESTAREAROOT="/afs/cern.ch/user/t/tiwillia/cmtuser/GEANT4/GEANT4_HEAD/Geant4/G4examples/extended/electromagnetic/G4TestEm3"
-
-export PLIST="Test/opt1NoApplyCuts"
+export PLIST="Test/opt1noapplycuts"
$TESTAREAROOT/x86_64-slc6-gcc48-opt/testEm3.exe $TESTAREAROOT/$PLIST/opt1noapplycuts.mac #generated the data if run with sim option
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/opt1noapplycuts.mac b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/opt1noapplycuts.mac
index b4ecbab38..cdd282be7 100644
--- a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/opt1noapplycuts.mac
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/opt1noapplycuts.mac
@@ -13,7 +13,7 @@
/testem/det/setAbsor 2 Scintillator 4 mm
/testem/det/setSizeYZ 20 cm
#
-/testem/phys/addPhysics emstandard_opt1nocuts
+/testem/phys/addPhysics emstandard_opt1nocuts
#
#/gun/particle e-
#/gun/energy 100 GeV
@@ -43,6 +43,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -57,6 +63,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -71,6 +83,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -85,6 +103,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -99,6 +123,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -113,7 +143,13 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
-/run/beamOn 10000
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
+/run/beamOn 10000
#
/gun/particle e-
/gun/energy 4.94 GeV
@@ -127,6 +163,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -141,6 +183,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -155,6 +203,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -169,6 +223,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -183,6 +243,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -197,6 +263,7 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+
/run/beamOn 10000
#
/gun/particle e-
@@ -211,6 +278,12 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
#
/gun/particle e-
@@ -225,4 +298,10 @@
/analysis/h1/set 12 66 0. 66. none #long. profile in absor2
/analysis/h1/set 21 102 0. 102. none #energy flow
/analysis/h1/set 22 102 0. 102. none #lateral energy leakage
+/analysis/h1/setActivation 23 true
+/analysis/h1/setActivation 24 true
+/analysis/h1/setActivation 25 true
+/analysis/h1/setActivation 26 true
+/analysis/h1/setActivation 27 true
+/analysis/h1/setActivation 28 true
/run/beamOn 10000
\ No newline at end of file
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/run.sh b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/run.sh
index 5c7c34fe8..5581ca76e 100755
--- a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/run.sh
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt1NoApplyCuts/run.sh
@@ -1,7 +1,11 @@
+#run with . run.sh to keep environment for SetupProject
+
if [ "$1" = "sim" ]
then
- ../../x86_64-slc6-gcc48-opt/testEm3.exe opt1noapplycuts.mac #generated the data if run with sim option
+ ../../x86_64-slc6-gcc48-opt/testEm3.exe opt1noapplycuts.mac
fi
+#There is a bug with ROOT 6 opening Geant4 created histograms, have to use old root for now
+SetupProject ROOT --version 72a
g++ -c `root-config --cflags` Plot.C
g++ -o Plot `root-config --glibs` Plot.o
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt2/opt2.mac b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt2/opt2.mac
index cafa24ce9..9748795fa 100644
--- a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt2/opt2.mac
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/Test/opt2/opt2.mac
@@ -13,7 +13,7 @@
/testem/det/setAbsor 2 Scintillator 4 mm
/testem/det/setSizeYZ 20 cm
#
-/testem/phys/addPhysics emstandard_opt2
+/testem/phys/addPhysics emstandard_opt1nocuts
#
#/gun/particle e-
#/gun/energy 100 GeV
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/cmt/copyPatchedSource.py b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/cmt/copyPatchedSource.py
new file mode 100755
index 000000000..c86cea199
--- /dev/null
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/cmt/copyPatchedSource.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+#
+# copyPatchedSource.py
+# Author Nigel Watson 14 Feb 2012
+#
+# Copy private G4 headers from package to install area
+#
+import os, sys, fnmatch, shutil
+
+def main():
+# Dir in which we keep updated source/headers.
+ # SRCNEW_Dir="../srcnew"
+ #SRCNEW_cross_sections_Dir = "../srcnew/cross_sections"
+#11/2014 SRCNEW_diffraction_Dir="../srcnew/diffraction"
+#11/2014 SRCNEW_management_Dir="../srcnew/management"
+
+
+# Dir for original source
+#11/2014 SRC_Dir="../hadronic/models/chiral_inv_phase_space/interface/src"
+# Dir for original FTF source to be replaced
+#11/2014 SRC_diffraction_Dir="../hadronic/models/parton_string/diffraction/src"
+#11/2014 SRC_management_Dir="../hadronic/models/parton_string/management/src"
+ #SRC_cross_sections_Dir="../"
+# Can we get cmt macro values in .python?
+# Expanded by cmt to be $(GEANT4_home)/$(GEANT4_installarea_prefix)/include
+# INSTALLAREA_project=os.environ['GEANT4_install_include']
+# INSTALLAREA_package=os.environ['G4PROCESSESROOT']+'/G4processes'
+ #INSTALLAREA_project = "../../../InstallArea/include"
+ #INSTALLAREA_package = "../G4processes"
+
+#$(G4processes_root)/G4processes
+
+# Find any files in the srcnew dir (CHIPS cross-section fixes)
+#11/2014 for newfile in os.listdir(SRCNEW_Dir):
+#11/2014 if fnmatch.fnmatch(newfile, '*.hh'):
+#11/2014#Replace old headers in install areas with new.
+#11/2014 fname = os.path.join(SRCNEW_Dir, newfile)
+#11/2014 shutil.copy2(fname,INSTALLAREA_project)
+#11/2014 shutil.copy2(fname,INSTALLAREA_package)
+#11/2014# Replace old .cc with new.
+#11/2014 if fnmatch.fnmatch(newfile, '*.cc'):
+#11/2014 fname = os.path.join(SRCNEW_Dir, newfile)
+#11/2014 shutil.copy2(fname,SRC_Dir)
+
+
+#11/2014# Find any files in the srcnew dir (FTF mass problem fixes)
+#11/2014 for newfile in os.listdir(SRCNEW_diffraction_Dir):
+#11/2014# Only .cc to replace for FTF fix.
+#11/2014# if fnmatch.fnmatch(newfile, '*.hh'):
+#11/2014#Replace old headers in install areas with new.
+#11/2014# fname = os.path.join(SRCNEW_Dir, newfile)
+#11/2014# shutil.copy2(fname,INSTALLAREA_project)
+#11/2014# shutil.copy2(fname,INSTALLAREA_package)
+#11/2014# Replace old .cc with new.
+#11/2014 if fnmatch.fnmatch(newfile, '*.cc'):
+#11/2014 fname = os.path.join(SRCNEW_diffraction_Dir, newfile)
+#11/2014 shutil.copy2(fname,SRC_diffraction_Dir)
+
+#11/2014# Find any files in the srcnew dir (string max. retries problem fixes)
+#11/2014 for newfile in os.listdir(SRCNEW_management_Dir):
+#11/2014# Only .cc to replace for max. retries fix.
+#11/2014# if fnmatch.fnmatch(newfile, '*.hh'):
+#11/2014#Replace old headers in install areas with new.
+#11/2014# fname = os.path.join(SRCNEW_Dir, newfile)
+#11/2014# shutil.copy2(fname,INSTALLAREA_project)
+#11/2014# shutil.copy2(fname,INSTALLAREA_package)
+#11/2014# Replace old .cc with new.
+#11/2014 if fnmatch.fnmatch(newfile, '*.cc'):
+#11/2014 fname = os.path.join(SRCNEW_management_Dir, newfile)
+#11/2014 shutil.copy2(fname,SRC_management_Dir)
+
+# Find any files in the srcnew dir (fixes missing assert in LHCb CLHEP version)
+ #for newfile in os.listdir(SRCNEW_cross_sections_Dir):
+# Only .cc to replace for CLHEP fix.
+# if fnmatch.fnmatch(newfile, '*.hh'):
+#Replace old headers in install areas with new.
+# fname = os.path.join(SRCNEW_Dir, newfile)
+# shutil.copy2(fname,INSTALLAREA_project)
+# shutil.copy2(fname,INSTALLAREA_package)
+# Replace old .cc with new.
+ # if fnmatch.fnmatch(newfile, '*.cc'):
+ # fname = os.path.join(SRCNEW_cross_sections_Dir, newfile)
+ # shutil.copy2(fname,SRC_cross_sections_Dir)
+
+#========================================TESTEM3 PRIVATE SRC COPYING========================================
+ SRCNEW_DIR="../srcnew/"
+ SRC_DIR="../src/"
+
+ for privatesource in os.listdir(SRCNEW_DIR):
+ if privatesource!='.svn':
+ fname=os.path.join(SRCNEW_DIR,privatesource)
+ shutil.copy2(fname,SRC_DIR)
+
+if __name__ == "__main__":
+ main()
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/doc/release.notes b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/doc/release.notes
index 34a181022..df8c0c2c9 100644
--- a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/doc/release.notes
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/doc/release.notes
@@ -3,16 +3,23 @@
! Responsible : Timothy Williams
! Purpose : Standalone Geant4 example used to test EM shower simulation
!-----------------------------------------------------------------------------
-
+================================================================================
+! 2015-05-07 - Timothy Williams
+-Added srcnew diractory to store private source files. These get copied over with copyPatchedsource.py in cmt after sources are copied from LCG area.
+================================================================================
+! 2014-12-07- Timothy Williams
+Updated Analysis code Plot.C, added more information and better formatting to plots, removed all binary files, added run.sh script to run entire analysis. Run as "run.sh sim" to run geant4 and create output.
+================================================================================
+! 2014-11-20 - Timothy Williams
+-added README to Test/opt2 exaplining analysis for LHCbPR.2
+================================================================================
+! 2014-11-19 - Timothy Williams
+-added Simulated data and analysis code in directory Test/opt2. This will be used for LHCbPR.
+================================================================================
! 2014-11-13 - Timothy Williams
-Initial Version.
-Full desription of example available at:- http://geant4.web.cern.ch/geant4/UserDocumentation/Doxygen/examples_doc/html/ExampleTestEm3.html
-2014-11-19 - Timothy Williams
--added Simulated data and analysis code in directory Test/opt2. This will be used for LHCbPR.
-2014-11-20 - Timothy Williams
--added README to Test/opt2 exaplining analysis for LHCbPR.2
-2015-12-07- Timothy Williams
-Updated Analysis code Plot.C, added more information and better formatting to plots, removed all binary files, added run.sh script to run entire analysis. Run as "run.sh sim" to run geant4 and create output.
\ No newline at end of file
+
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/scripts/OUTPUT_3.000 b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/scripts/OUTPUT_3.000
index f8fff504f9e33673993832f264fb3db8723679a1..b9f997dd90a6e5c19fa6984ab39f4e972bf32181 100644
GIT binary patch
delta 18
XcmZ3+xQuZ^9Y?@q@%JD+u}uL0LkS0N
delta 18
XcmZ3+xQuZ^9f#lTg6kkWu}uL0OCtz;
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/HistoManager.cc b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/HistoManager.cc
new file mode 100644
index 000000000..20269737d
--- /dev/null
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/HistoManager.cc
@@ -0,0 +1,107 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+// $Id$
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "HistoManager.hh"
+#include "G4UnitsTable.hh"
+#include "DetectorConstruction.hh"
+#include<map>
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+HistoManager::HistoManager()
+ : fFileName("testem3")
+{
+ Book();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+HistoManager::~HistoManager()
+{
+ delete G4AnalysisManager::Instance();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void HistoManager::Book()
+{
+ // Create or get analysis manager
+ // The choice of analysis technology is done via selection of a namespace
+ // in HistoManager.hh
+ G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
+ analysisManager->SetFileName(fFileName);
+ analysisManager->SetVerboseLevel(1);
+ analysisManager->SetActivation(true); // enable inactivation of histograms
+
+ // Define histograms start values
+
+ const G4String id[] = { "0", "1", "2", "3", "4", "5", "6", "7", "8", "9",
+ "10","11","12","13","14","15","16","17","18","19",
+ "20","21","22"};
+ G4String title;
+ typedef std::map<G4String,G4String> HistPropertyCont;
+ HistPropertyCont Properties;
+
+
+ // define new histogram properties
+ Properties.insert(std::make_pair("23","Photons in Lead Shower Profile"));
+ Properties.insert(std::make_pair("24","Electrons in Lead Shower Profile"));
+ Properties.insert(std::make_pair("25","Positrons in Lead Shower Profile"));
+ Properties.insert(std::make_pair("26","Photons in Scintillator Shower Profile"));
+ Properties.insert(std::make_pair("27","Electrons in Scintillator Shower Profile"));
+ Properties.insert(std::make_pair("28","Positrons in Scintillator Shower Profile"));
+ // Default values (to be reset via /analysis/h1/set command)
+ G4int nbins = 100;
+ G4double vmin = 0.;
+ G4double vmax = 100.;
+
+ // Create all histograms as inactivated
+ // as we have not yet set nbins, vmin, vmax
+ for (G4int k=0; k<MaxHisto; k++) {
+ if (k < MaxAbsor) title = "Edep in absorber " + id[k];
+ if (k > MaxAbsor) title = "Edep longit. profile (MeV/event) in absorber "
+ + id[k-MaxAbsor];
+ if (k == 2*MaxAbsor+1) title = "energy flow (MeV/event)";
+ if (k == 2*MaxAbsor+2) title = "lateral energy leak (MeV/event)";
+ G4int ih = analysisManager->CreateH1(id[k], title, nbins, vmin, vmax);
+ std::cout<<"histo k number= "<<k<<std::endl;
+ analysisManager->SetActivation(G4VAnalysisManager::kH1, ih, false);
+ }
+
+ //G4int i=analysisManager->CreateH1("23","A Test Histogram",66,1,67);
+ // analysisManager->SetActivation(G4VAnalysisManager::kH1,i,false);
+
+ HistPropertyCont::iterator h=Properties.begin();
+ for(;h!=Properties.end();++h){
+ G4int i= analysisManager->CreateH1((*h).first,(*h).second,66,1,67);
+ analysisManager->SetActivation(G4VAnalysisManager::kH1,i,false);
+ }
+}
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/PhysicsList.cc b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/PhysicsList.cc
new file mode 100644
index 000000000..6d1adb9a2
--- /dev/null
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/PhysicsList.cc
@@ -0,0 +1,357 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+/// \file electromagnetic/TestEm3/src/PhysicsList.cc
+/// \brief Implementation of the PhysicsList class
+//
+// $Id$
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "PhysicsList.hh"
+#include "PhysicsListMessenger.hh"
+
+#include "PhysListEmStandard.hh"
+#include "PhysListEmStandardWVI.hh"
+
+#include "G4EmStandardPhysics.hh"
+#include "G4EmStandardPhysics_option1.hh"
+#include "G4EmStandardPhysics_option2.hh"
+#include "G4EmStandardPhysics_option3.hh"
+#include "G4EmStandardPhysics_option4.hh"
+#include "G4EmStandardPhysics_option1NoApplyCuts.hh"
+#include "G4EmLivermorePhysics.hh"
+#include "G4EmPenelopePhysics.hh"
+
+#include "G4UnitsTable.hh"
+#include "G4SystemOfUnits.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//a comment
+PhysicsList::PhysicsList() : G4VModularPhysicsList()
+{
+ fCurrentDefaultCut = 1.0*mm;
+ fCutForGamma = fCurrentDefaultCut;
+ fCutForElectron = fCurrentDefaultCut;
+ fCutForPositron = fCurrentDefaultCut;
+
+ fMessenger = new PhysicsListMessenger(this);
+
+ SetVerboseLevel(1);
+
+ // EM physics
+ fEmName = G4String("local");
+ fEmPhysicsList = new PhysListEmStandard(fEmName);
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+PhysicsList::~PhysicsList()
+{
+ delete fMessenger;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Bosons
+#include "G4ChargedGeantino.hh"
+#include "G4Geantino.hh"
+#include "G4Gamma.hh"
+#include "G4OpticalPhoton.hh"
+
+// leptons
+#include "G4MuonPlus.hh"
+#include "G4MuonMinus.hh"
+#include "G4NeutrinoMu.hh"
+#include "G4AntiNeutrinoMu.hh"
+
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+#include "G4NeutrinoE.hh"
+#include "G4AntiNeutrinoE.hh"
+
+// Mesons
+#include "G4PionPlus.hh"
+#include "G4PionMinus.hh"
+#include "G4PionZero.hh"
+#include "G4Eta.hh"
+#include "G4EtaPrime.hh"
+
+#include "G4KaonPlus.hh"
+#include "G4KaonMinus.hh"
+#include "G4KaonZero.hh"
+#include "G4AntiKaonZero.hh"
+#include "G4KaonZeroLong.hh"
+#include "G4KaonZeroShort.hh"
+
+// Baryons
+#include "G4Proton.hh"
+#include "G4AntiProton.hh"
+#include "G4Neutron.hh"
+#include "G4AntiNeutron.hh"
+
+// Nuclei
+#include "G4Deuteron.hh"
+#include "G4Triton.hh"
+#include "G4Alpha.hh"
+#include "G4GenericIon.hh"
+
+void PhysicsList::ConstructParticle()
+{
+// pseudo-particles
+ G4Geantino::GeantinoDefinition();
+ G4ChargedGeantino::ChargedGeantinoDefinition();
+
+// gamma
+ G4Gamma::GammaDefinition();
+
+// optical photon
+ G4OpticalPhoton::OpticalPhotonDefinition();
+
+// leptons
+ G4Electron::ElectronDefinition();
+ G4Positron::PositronDefinition();
+ G4MuonPlus::MuonPlusDefinition();
+ G4MuonMinus::MuonMinusDefinition();
+
+ G4NeutrinoE::NeutrinoEDefinition();
+ G4AntiNeutrinoE::AntiNeutrinoEDefinition();
+ G4NeutrinoMu::NeutrinoMuDefinition();
+ G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
+
+// mesons
+ G4PionPlus::PionPlusDefinition();
+ G4PionMinus::PionMinusDefinition();
+ G4PionZero::PionZeroDefinition();
+ G4Eta::EtaDefinition();
+ G4EtaPrime::EtaPrimeDefinition();
+ G4KaonPlus::KaonPlusDefinition();
+ G4KaonMinus::KaonMinusDefinition();
+ G4KaonZero::KaonZeroDefinition();
+ G4AntiKaonZero::AntiKaonZeroDefinition();
+ G4KaonZeroLong::KaonZeroLongDefinition();
+ G4KaonZeroShort::KaonZeroShortDefinition();
+
+// barions
+ G4Proton::ProtonDefinition();
+ G4AntiProton::AntiProtonDefinition();
+ G4Neutron::NeutronDefinition();
+ G4AntiNeutron::AntiNeutronDefinition();
+
+// ions
+ G4Deuteron::DeuteronDefinition();
+ G4Triton::TritonDefinition();
+ G4Alpha::AlphaDefinition();
+ G4GenericIon::GenericIonDefinition();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4ProcessManager.hh"
+#include "G4Decay.hh"
+
+void PhysicsList::ConstructProcess()
+{
+ AddTransportation();
+
+ // electromagnetic Physics List
+ //
+ fEmPhysicsList->ConstructProcess();
+
+ // Add Decay Process
+ //
+ G4Decay* fDecayProcess = new G4Decay();
+
+ theParticleIterator->reset();
+ while( (*theParticleIterator)() ){
+ G4ParticleDefinition* particle = theParticleIterator->value();
+ G4ProcessManager* pmanager = particle->GetProcessManager();
+
+ if (fDecayProcess->IsApplicable(*particle)) {
+
+ pmanager ->AddProcess(fDecayProcess);
+
+ // set ordering for PostStepDoIt and AtRestDoIt
+ pmanager ->SetProcessOrdering(fDecayProcess, idxPostStep);
+ pmanager ->SetProcessOrdering(fDecayProcess, idxAtRest);
+
+ }
+ }
+
+ // stepLimitation (as a full process)
+ //
+ AddStepMax();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void PhysicsList::AddPhysicsList(const G4String& name)
+{
+ if (verboseLevel>1) {
+ G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
+ }
+
+ if (name == fEmName) return;
+
+ if (name == "local") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new PhysListEmStandard(name);
+
+ } else if (name == "emstandard_opt0") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmStandardPhysics();
+
+ } else if (name == "emstandard_opt1") {
+ std::cout<<"opt1 physics list found"<<std::endl;
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmStandardPhysics_option1();
+
+ } else if (name == "emstandard_opt2") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmStandardPhysics_option2();
+
+ } else if (name == "emstandard_opt3") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmStandardPhysics_option3();
+
+ } else if (name == "emstandard_opt4") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmStandardPhysics_option4();
+
+ } else if (name == "emstandard_opt1nocuts")
+ {
+ fEmName = name;
+ std::cout<<"Passed if statement in Physicslist.cc"<<std::endl;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmStandardPhysics_option1NoApplyCuts();
+ }
+ else if (name == "standardWVI") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new PhysListEmStandardWVI(name);
+
+ } else if (name == "emlivermore") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmLivermorePhysics();
+
+ } else if (name == "empenelope") {
+
+ fEmName = name;
+ delete fEmPhysicsList;
+ fEmPhysicsList = new G4EmPenelopePhysics();
+
+ } else {
+
+ G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
+ << " is not defined"
+ << G4endl;
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "StepMax.hh"
+
+void PhysicsList::AddStepMax()
+{
+ // Step limitation seen as a process
+ fStepMaxProcess = new StepMax();
+
+ theParticleIterator->reset();
+ while ((*theParticleIterator)()){
+ G4ParticleDefinition* particle = theParticleIterator->value();
+ G4ProcessManager* pmanager = particle->GetProcessManager();
+
+ if (fStepMaxProcess->IsApplicable(*particle))
+ {
+ pmanager ->AddDiscreteProcess(fStepMaxProcess);
+ }
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void PhysicsList::SetCuts()
+{
+ if (verboseLevel >0) {
+ G4cout << "PhysicsList::SetCuts:";
+ G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
+ }
+
+ // set cut values for gamma at first and for e- second and next for e+,
+ // because some processes for e+/e- need cut values for gamma
+ SetCutValue(fCutForGamma, "gamma");
+ SetCutValue(fCutForElectron, "e-");
+ SetCutValue(fCutForPositron, "e+");
+
+ // Cut for proton not used in EM processes except single scattering
+ // so electron cut is used in this example
+ SetCutValue(fCutForElectron, "proton");
+
+ if (verboseLevel>0) DumpCutValuesTable();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void PhysicsList::SetCutForGamma(G4double cut)
+{
+ fCutForGamma = cut;
+ SetParticleCuts(fCutForGamma, G4Gamma::Gamma());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void PhysicsList::SetCutForElectron(G4double cut)
+{
+ fCutForElectron = cut;
+ SetParticleCuts(fCutForElectron, G4Electron::Electron());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void PhysicsList::SetCutForPositron(G4double cut)
+{
+ fCutForPositron = cut;
+ SetParticleCuts(fCutForPositron, G4Positron::Positron());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/RunAction.cc b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/RunAction.cc
new file mode 100644
index 000000000..db3bc5501
--- /dev/null
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/RunAction.cc
@@ -0,0 +1,422 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+/// \file electromagnetic/TestEm3/src/RunAction.cc
+/// \brief Implementation of the RunAction class
+//
+// $Id$
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "RunAction.hh"
+
+#include "PrimaryGeneratorAction.hh"
+#include "RunActionMessenger.hh"
+#include "HistoManager.hh"
+#include "EmAcceptance.hh"
+
+#include "G4Run.hh"
+#include "G4RunManager.hh"
+
+#include "G4ParticleTable.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4Track.hh"
+#include "G4Gamma.hh"
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+#include "G4ProductionCutsTable.hh"
+#include "G4LossTableManager.hh"
+
+#include "G4UnitsTable.hh"
+#include "G4SystemOfUnits.hh"
+
+#include "Randomize.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
+:fDetector(det), fPrimary(prim), fRunMessenger(0), fHistoManager(0)
+{
+ fRunMessenger = new RunActionMessenger(this);
+ fHistoManager = new HistoManager();
+ fApplyLimit = false;
+
+ fChargedStep = fNeutralStep = 0.0;
+
+ for (G4int k=0; k<MaxAbsor; k++) { fEdeptrue[k] = fRmstrue[k] = 1.;
+ fLimittrue[k] = DBL_MAX;
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+RunAction::~RunAction()
+{
+ delete fRunMessenger;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void RunAction::BeginOfRunAction(const G4Run* aRun)
+{
+ G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
+
+ // save Rndm status
+ //
+ ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
+ CLHEP::HepRandom::showEngineStatus();
+
+ //initialize cumulative quantities
+ //
+ for (G4int k=0; k<MaxAbsor; k++) {
+ fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
+ fEnergyDeposit[k].clear();
+ }
+
+ fChargedStep = fNeutralStep = 0.0;
+
+ fN_gamma = 0;
+ fN_elec = 0;
+ fN_pos = 0;
+
+ //initialize Eflow
+ //
+ G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
+ fEnergyFlow.resize(nbPlanes);
+ fLateralEleak.resize(nbPlanes);
+ for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }
+
+ //histograms
+ //
+ G4AnalysisManager* analysis = G4AnalysisManager::Instance();
+ if (analysis->IsActive()) analysis->OpenFile();
+
+ //example of print dEdx tables
+ //
+ ////PrintDedxTables();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void RunAction::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
+{
+ //accumulate statistic with restriction
+ //
+ if(fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
+ fSumEAbs[kAbs] += EAbs; fSum2EAbs[kAbs] += EAbs*EAbs;
+ fSumLAbs[kAbs] += LAbs; fSum2LAbs[kAbs] += LAbs*LAbs;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
+void RunAction::EndOfRunAction(const G4Run* aRun)
+{
+ G4int nEvt = aRun->GetNumberOfEvent();
+ G4double norm = G4double(nEvt);
+ if(norm > 0) norm = 1./norm;
+ G4double qnorm = std::sqrt(norm);
+
+ fChargedStep *= norm;
+ fNeutralStep *= norm;
+
+ //compute and print statistic
+ //
+ G4double beamEnergy = fPrimary->GetParticleGun()->GetParticleEnergy();
+ G4double sqbeam = std::sqrt(beamEnergy/GeV);
+
+ G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
+ G4double MeanLAbs,MeanLAbs2,rmsLAbs;
+
+ std::ios::fmtflags mode = G4cout.flags();
+ G4int prec = G4cout.precision(2);
+ G4cout << "\n------------------------------------------------------------\n";
+ G4cout << std::setw(14) << "material"
+ << std::setw(17) << "Edep RMS"
+ << std::setw(33) << "sqrt(E0(GeV))*rmsE/Emean"
+ << std::setw(23) << "total tracklen \n \n";
+
+ for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++)
+ {
+ MeanEAbs = fSumEAbs[k]*norm;
+ MeanEAbs2 = fSum2EAbs[k]*norm;
+ rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
+ //G4cout << "k= " << k << " RMS= " << rmsEAbs
+ // << " fApplyLimit: " << fApplyLimit << G4endl;
+ if(fApplyLimit) {
+ G4int nn = 0;
+ G4double sume = 0.0;
+ G4double sume2 = 0.0;
+ // compute trancated means
+ G4double lim = rmsEAbs * 2.5;
+ for(G4int i=0; i<nEvt; i++) {
+ G4double e = (fEnergyDeposit[k])[i];
+ if(std::abs(e - MeanEAbs) < lim) {
+ sume += e;
+ sume2 += e*e;
+ nn++;
+ }
+ }
+ G4double norm1 = G4double(nn);
+ if(norm1 > 0.0) norm1 = 1.0/norm1;
+ MeanEAbs = sume*norm1;
+ MeanEAbs2 = sume2*norm1;
+ rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
+ }
+
+ resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
+ rmsres = resolution*qnorm;
+
+ // Save mean and RMS
+ fSumEAbs[k] = MeanEAbs;
+ fSum2EAbs[k] = rmsEAbs;
+
+ MeanLAbs = fSumLAbs[k]*norm;
+ MeanLAbs2 = fSum2LAbs[k]*norm;
+ rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
+
+ //print
+ //
+ G4cout
+ << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
+ << std::setprecision(5)
+ << std::setw(6) << G4BestUnit(MeanEAbs,"Energy") << " : "
+ << std::setprecision(4)
+ << std::setw(5) << G4BestUnit( rmsEAbs,"Energy")
+ << std::setw(10) << resolution << " +- "
+ << std::setw(5) << rmsres << " %"
+ << std::setprecision(3)
+ << std::setw(10) << G4BestUnit(MeanLAbs,"Length") << " +- "
+ << std::setw(4) << G4BestUnit( rmsLAbs,"Length")
+ << G4endl;
+ }
+ G4cout << "\n------------------------------------------------------------\n";
+
+ G4cout << " Beam particle "
+ << fPrimary->GetParticleGun()->
+ GetParticleDefinition()->GetParticleName()
+ << " E = " << G4BestUnit(beamEnergy,"Energy") << G4endl;
+ G4cout << " Mean number of gamma " << (G4double)fN_gamma*norm << G4endl;
+ G4cout << " Mean number of e- " << (G4double)fN_elec*norm << G4endl;
+ G4cout << " Mean number of e+ " << (G4double)fN_pos*norm << G4endl;
+ G4cout << std::setprecision(6)
+ << " Mean number of charged steps " << fChargedStep << G4endl;
+ G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl;
+ G4cout << "------------------------------------------------------------\n";
+
+ //Energy flow
+ //
+ G4AnalysisManager* analysis = G4AnalysisManager::Instance();
+ G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor());
+ for (G4int Id=1; Id<=Idmax+1; Id++) {
+ analysis->FillH1(2*MaxAbsor+1, (G4double)Id, fEnergyFlow[Id]);
+ analysis->FillH1(2*MaxAbsor+2, (G4double)Id, fLateralEleak[Id]);
+ }
+
+ //Energy deposit from energy flow balance
+ //
+ G4double EdepTot[MaxAbsor];
+ for (G4int k=0; k<MaxAbsor; k++) EdepTot[k] = 0.;
+
+ G4int nbOfAbsor = fDetector->GetNbOfAbsor();
+ for (G4int Id=1; Id<=Idmax; Id++) {
+ G4int iAbsor = Id%nbOfAbsor; if (iAbsor==0) iAbsor = nbOfAbsor;
+ EdepTot [iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id+1] - fLateralEleak[Id]);
+ }
+
+ G4cout << std::setprecision(3)
+ << "\n Energy deposition from Energy flow balance : \n"
+ << std::setw(10) << " material \t Total Edep \n \n";
+ G4cout.precision(6);
+
+ for (G4int k=1; k<=nbOfAbsor; k++) {
+ EdepTot [k] *= norm;
+ G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
+ << "\t " << G4BestUnit(EdepTot [k],"Energy") << "\n";
+ }
+
+ G4cout << "\n------------------------------------------------------------\n"
+ << G4endl;
+
+ G4cout.setf(mode,std::ios::floatfield);
+ G4cout.precision(prec);
+
+ // Acceptance
+ EmAcceptance acc;
+ G4bool isStarted = false;
+ for (G4int j=1; j<=fDetector->GetNbOfAbsor(); j++) {
+ if (fLimittrue[j] < DBL_MAX) {
+ if (!isStarted) {
+ acc.BeginOfAcceptance("Sampling Calorimeter",nEvt);
+ isStarted = true;
+ }
+ MeanEAbs = fSumEAbs[j];
+ rmsEAbs = fSum2EAbs[j];
+ G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
+ acc.EmAcceptanceGauss("Edep"+mat, nEvt, MeanEAbs,
+ fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
+ acc.EmAcceptanceGauss("Erms"+mat, nEvt, rmsEAbs,
+ fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
+ }
+ }
+ if(isStarted) acc.EndOfAcceptance();
+
+ //normalize histograms
+ //
+ for (G4int ih = MaxAbsor+1; ih < MaxHisto; ih++) {
+ analysis->ScaleH1(ih,norm/MeV);
+ std::cout<<"MeV= "<<MeV<<std::endl;
+ }
+
+ for (G4int i= 23;i <=28;++i)
+ {
+ analysis->ScaleH1(i,norm/MeV);
+ }
+
+
+ //save histograms
+ if (analysis->IsActive()) {
+ std::cout<<"got to before write"<<std::endl;
+ analysis->Write();
+ analysis->CloseFile();
+ std::cout<<"got to after close"<<std::endl;
+ }
+
+ // show Rndm status
+ CLHEP::HepRandom::showEngineStatus();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void RunAction::PrintDedxTables()
+{
+ //Print dE/dx tables with binning identical to the Geant3 JMATE bank.
+ //The printout is readable as Geant3 ffread data cards (by the program g4mat).
+ //
+ const G4double tkmin=10*keV, tkmax=10*TeV;
+ const G4int nbin=90;
+ G4double tk[nbin];
+
+ const G4int ncolumn = 5;
+
+ //compute the kinetic energies
+ //
+ const G4double dp = std::log10(tkmax/tkmin)/nbin;
+ const G4double dt = std::pow(10.,dp);
+ tk[0] = tkmin;
+ for (G4int i=1; i<nbin; ++i) tk[i] = tk[i-1]*dt;
+
+ //print the kinetic energies
+ //
+ std::ios::fmtflags mode = G4cout.flags();
+ G4cout.setf(std::ios::fixed,std::ios::floatfield);
+ G4int prec = G4cout.precision(3);
+
+ G4cout << "\n kinetic energies \n ";
+ for (G4int j=0; j<nbin; ++j) {
+ G4cout << G4BestUnit(tk[j],"Energy") << "\t";
+ if ((j+1)%ncolumn == 0) G4cout << "\n ";
+ }
+ G4cout << G4endl;
+
+ //print the dE/dx tables
+ //
+ G4cout.setf(std::ios::scientific,std::ios::floatfield);
+
+ G4ParticleDefinition*
+ part = fPrimary->GetParticleGun()->GetParticleDefinition();
+
+ G4ProductionCutsTable* theCoupleTable =
+ G4ProductionCutsTable::GetProductionCutsTable();
+ size_t numOfCouples = theCoupleTable->GetTableSize();
+ const G4MaterialCutsCouple* couple = 0;
+
+ for (G4int iab=1;iab <= fDetector->GetNbOfAbsor(); iab++)
+ {
+ G4Material* mat = fDetector->GetAbsorMaterial(iab);
+ G4int index = 0;
+ for (size_t i=0; i<numOfCouples; i++) {
+ couple = theCoupleTable->GetMaterialCutsCouple(i);
+ if (couple->GetMaterial() == mat) {index = i; break;}
+ }
+ G4cout << "\nLIST";
+ G4cout << "\nC \nC dE/dx (MeV/cm) for " << part->GetParticleName()
+ << " in " << mat ->GetName() << "\nC";
+ G4cout << "\nKINE (" << part->GetParticleName() << ")";
+ G4cout << "\nMATE (" << mat ->GetName() << ")";
+ G4cout.precision(2);
+ G4cout << "\nERAN " << tkmin/GeV << " (ekmin)\t"
+ << tkmax/GeV << " (ekmax)\t"
+ << nbin << " (nekbin)";
+ G4double cutgam =
+ (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
+ if (cutgam < tkmin) cutgam = tkmin;
+ if (cutgam > tkmax) cutgam = tkmax;
+ G4double cutele =
+ (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
+ if (cutele < tkmin) cutele = tkmin;
+ if (cutele > tkmax) cutele = tkmax;
+ G4cout << "\nCUTS " << cutgam/GeV << " (cutgam)\t"
+ << cutele/GeV << " (cutele)";
+
+ G4cout.precision(6);
+ G4cout << "\nG4VAL \n ";
+ for (G4int l=0;l<nbin; ++l)
+ {
+ G4double dedx = G4LossTableManager::Instance()
+ ->GetDEDX(part,tk[l],couple);
+ G4cout << dedx/(MeV/cm) << "\t";
+ if ((l+1)%ncolumn == 0) G4cout << "\n ";
+ }
+ G4cout << G4endl;
+ }
+
+ G4cout.precision(prec);
+ G4cout.setf(mode,std::ios::floatfield);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void RunAction::AddSecondaryTrack(const G4Track* track)
+{
+ const G4ParticleDefinition* d = track->GetDefinition();
+ if(d == G4Gamma::Gamma()) { ++fN_gamma; }
+ else if (d == G4Electron::Electron()) { ++fN_elec; }
+ else if (d == G4Positron::Positron()) { ++fN_pos; }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void RunAction::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
+{
+ if (i>=0 && i<MaxAbsor) {
+ fEdeptrue [i] = edep;
+ fRmstrue [i] = rms;
+ fLimittrue[i] = lim;
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/SteppingAction.cc b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/SteppingAction.cc
new file mode 100644
index 000000000..f48c70b55
--- /dev/null
+++ b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/srcnew/SteppingAction.cc
@@ -0,0 +1,174 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+/// \file electromagnetic/TestEm3/src/SteppingAction.cc
+/// \brief Implementation of the SteppingAction class
+//
+// $Id$
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "SteppingAction.hh"
+
+#include "DetectorConstruction.hh"
+#include "RunAction.hh"
+#include "EventAction.hh"
+#include "HistoManager.hh"
+
+#include "G4Step.hh"
+#include "G4Positron.hh"
+#include "G4RunManager.hh"
+#include "G4PhysicalConstants.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+SteppingAction::SteppingAction(DetectorConstruction* det, RunAction* run,
+ EventAction* evt)
+:G4UserSteppingAction(),fDetector(det),fRunAct(run),fEventAct(evt)
+{ }
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+SteppingAction::~SteppingAction()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void SteppingAction::UserSteppingAction(const G4Step* aStep)
+{
+ //track informations
+ const G4StepPoint* prePoint = aStep->GetPreStepPoint();
+ const G4StepPoint* endPoint = aStep->GetPostStepPoint();
+ const G4ParticleDefinition* particle = aStep->GetTrack()->GetDefinition();
+
+ //if World, return
+ //
+ G4VPhysicalVolume* volume = prePoint->GetTouchableHandle()->GetVolume();
+ //if sum of absorbers do not fill exactly a layer: check material, not volume.
+ G4Material* mat = volume->GetLogicalVolume()->GetMaterial();
+ if (mat == fDetector->GetWorldMaterial()) return;
+
+ //here we are in an absorber. Locate it
+ //
+ G4int absorNum = prePoint->GetTouchableHandle()->GetCopyNumber(0);
+ G4int layerNum = prePoint->GetTouchableHandle()->GetCopyNumber(1);
+
+ // collect energy deposit taking into account track weight
+ G4double edep = aStep->GetTotalEnergyDeposit()*aStep->GetTrack()->GetWeight();
+
+ // collect step length of charged particles
+ G4double stepl = 0.;
+ if (particle->GetPDGCharge() != 0.) {
+ stepl = aStep->GetStepLength();
+ fRunAct->AddChargedStep();
+ } else { fRunAct->AddNeutralStep(); }
+
+ // G4cout << "Nabs= " << absorNum << " edep(keV)= " << edep << G4endl;
+
+ // sum up per event
+ fEventAct->SumEnergy(absorNum,edep,stepl);
+ G4int Encoding=particle->GetPDGEncoding();
+ G4AnalysisManager* analysisManager=G4AnalysisManager::Instance();
+
+ //longitudinal profile of edep per absorber
+ if (edep>0.) {
+ G4AnalysisManager::Instance()->FillH1(MaxAbsor+absorNum,G4double(layerNum+1), edep);
+ if(absorNum==1){
+ if(Encoding==11){
+ analysisManager->FillH1(24,G4double(layerNum+1),edep);
+ }
+ else if(Encoding==-11){
+ analysisManager->FillH1(25,G4double(layerNum+1),edep);
+ }
+ else if(Encoding==22){
+ analysisManager->FillH1(23,G4double(layerNum+1),edep);
+ }
+ }
+ if(absorNum==2){
+ if(Encoding==11){
+ analysisManager->FillH1(27,G4double(layerNum+1),edep);
+ }
+ else if(Encoding==-11){
+ analysisManager->FillH1(28,G4double(layerNum+1),edep);
+ }
+ else if(Encoding==22){
+ analysisManager->FillH1(26,G4double(layerNum+1),edep);
+ }
+ }
+ if((Encoding!=11)&&(Encoding!=22)&&(Encoding!=-11)){
+ std::cout<<"Different particle created with id: "<<Encoding<<std::endl;
+ }
+ }
+
+ //energy flow
+ //
+ // unique identificator of layer+absorber
+ G4int Idnow = (fDetector->GetNbOfAbsor())*layerNum + absorNum;
+ G4int plane;
+ //
+ //leaving the absorber ?
+ if (endPoint->GetStepStatus() == fGeomBoundary) {
+ G4ThreeVector position = endPoint->GetPosition();
+ G4ThreeVector direction = endPoint->GetMomentumDirection();
+ G4double sizeYZ = 0.5*fDetector->GetCalorSizeYZ();
+ G4double Eflow = endPoint->GetKineticEnergy();
+ if (particle == G4Positron::Positron()) Eflow += 2*electron_mass_c2;
+ if ((std::abs(position.y()) >= sizeYZ) || (std::abs(position.z()) >= sizeYZ))
+ fRunAct->SumLateralEleak(Idnow, Eflow);
+ else if (direction.x() >= 0.) fRunAct->SumEnergyFlow(plane=Idnow+1, Eflow);
+ else fRunAct->SumEnergyFlow(plane=Idnow, -Eflow);
+ }
+
+//// example of Birk attenuation
+///G4double destep = aStep->GetTotalEnergyDeposit();
+///G4double response = BirksAttenuation(aStep);
+///G4cout << " Destep: " << destep/keV << " keV"
+/// << " response after Birks: " << response/keV << " keV" << G4endl;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double SteppingAction::BirksAttenuation(const G4Step* aStep)
+{
+ //Example of Birk attenuation law in organic scintillators.
+ //adapted from Geant3 PHYS337. See MIN 80 (1970) 239-244
+ //
+ G4Material* material = aStep->GetTrack()->GetMaterial();
+ G4double birk1 = material->GetIonisation()->GetBirksConstant();
+ G4double destep = aStep->GetTotalEnergyDeposit();
+ G4double stepl = aStep->GetStepLength();
+ G4double charge = aStep->GetTrack()->GetDefinition()->GetPDGCharge();
+ //
+ G4double response = destep;
+ if (birk1*destep*stepl*charge != 0.)
+ {
+ response = destep/(1. + birk1*destep/stepl);
+ }
+ return response;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
diff --git a/Geant4/G4processes/cmt/copyPatchedSource.py b/Geant4/G4processes/cmt/copyPatchedSource.py
index 8d31eb18f..464a6b549 100755
--- a/Geant4/G4processes/cmt/copyPatchedSource.py
+++ b/Geant4/G4processes/cmt/copyPatchedSource.py
@@ -25,11 +25,10 @@ def main():
# Expanded by cmt to be $(GEANT4_home)/$(GEANT4_installarea_prefix)/include
# INSTALLAREA_project=os.environ['GEANT4_install_include']
# INSTALLAREA_package=os.environ['G4PROCESSESROOT']+'/G4processes'
- INSTALLAREA_project = "../../../InstallArea/include"
+ INSTALLAREA_project = "../../../InstallArea/" + os.environ['CMTCONFIG'] + "/include"
INSTALLAREA_package = "../G4processes"
#$(G4processes_root)/G4processes
-
# Find any files in the srcnew dir (CHIPS cross-section fixes)
#11/2014 for newfile in os.listdir(SRCNEW_Dir):
#11/2014 if fnmatch.fnmatch(newfile, '*.hh'):
@@ -72,11 +71,11 @@ def main():
# Find any files in the srcnew dir (fixes missing assert in LHCb CLHEP version)
for newfile in os.listdir(SRCNEW_cross_sections_Dir):
# Only .cc to replace for CLHEP fix.
-# if fnmatch.fnmatch(newfile, '*.hh'):
+ if fnmatch.fnmatch(newfile, '*.hh'):
#Replace old headers in install areas with new.
-# fname = os.path.join(SRCNEW_Dir, newfile)
-# shutil.copy2(fname,INSTALLAREA_project)
-# shutil.copy2(fname,INSTALLAREA_package)
+ fname = os.path.join(SRCNEW_cross_sections_Dir, newfile) #TW changed SRCNEW_Dir for SRCNEW_cross_sections_Dir
+ shutil.copy2(fname,INSTALLAREA_project)
+ shutil.copy2(fname,INSTALLAREA_package)
# Replace old .cc with new.
if fnmatch.fnmatch(newfile, '*.cc'):
fname = os.path.join(SRCNEW_cross_sections_Dir, newfile)
diff --git a/Geant4/G4processes/cmt/requirements b/Geant4/G4processes/cmt/requirements
index 2526f55e9..eaabf9065 100755
--- a/Geant4/G4processes/cmt/requirements
+++ b/Geant4/G4processes/cmt/requirements
@@ -3,7 +3,7 @@
# Maintainer : Gloria CORTI
#============================================================================
package G4processes
-version v9r0
+version v9r1
# Structure, i.e. directories to process.
#============================================================================
diff --git a/Geant4/G4processes/doc/release.notes b/Geant4/G4processes/doc/release.notes
index bc3cd26a1..5b1eadf27 100755
--- a/Geant4/G4processes/doc/release.notes
+++ b/Geant4/G4processes/doc/release.notes
@@ -3,9 +3,15 @@
! Responsible : Gloria CORTI, Nigel Watson
! Purpose :
!-----------------------------------------------------------------------------
+
+! ======================= G4processes v9r1 2015-04-24 ========================
+! 2015-04-24 Nigel Watson
+ - Include new alternative sources from Witek for kaon cross-sections
+
! ======================= G4processes v9r0 2015-03-06 ========================
! 2015-03-06 Nigel Watson
- Ready for G4 9.6p4 releease
+
! 2014-11-05 - Nigel Watson/Tim Williams
- Updated copyPatchedSource.py for 9.6.p03 (only single source to fix)
diff --git a/Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.cc b/Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.cc
new file mode 100644
index 000000000..cf4b7cc05
--- /dev/null
+++ b/Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.cc
@@ -0,0 +1,1701 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+// author: V. Grichine
+//
+// 25.04.12 V. Grichine - first implementation
+
+#include "G4ComponentGGHadronNucleusXsc.hh"
+
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4DynamicParticle.hh"
+#include "G4HadronNucleonXsc.hh"
+
+
+//////////////////////////////////////////////////////////////////////////////
+//
+
+G4ComponentGGHadronNucleusXsc::G4ComponentGGHadronNucleusXsc()
+ : G4VComponentCrossSection("Glauber-Gribov"),
+// fUpperLimit(100000*GeV),
+ fLowerLimit(10.*MeV),// fLowerLimit(3*GeV),
+ fRadiusConst(1.08*fermi), // 1.1, 1.3 ?
+ fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
+ fDiffractionXsc(0.0)
+// , fHadronNucleonXsc(0.0)
+{
+ theGamma = G4Gamma::Gamma();
+ theProton = G4Proton::Proton();
+ theNeutron = G4Neutron::Neutron();
+ theAProton = G4AntiProton::AntiProton();
+ theANeutron = G4AntiNeutron::AntiNeutron();
+ thePiPlus = G4PionPlus::PionPlus();
+ thePiMinus = G4PionMinus::PionMinus();
+ thePiZero = G4PionZero::PionZero();
+ theKPlus = G4KaonPlus::KaonPlus();
+ theKMinus = G4KaonMinus::KaonMinus();
+ theK0S = G4KaonZeroShort::KaonZeroShort();
+ theK0L = G4KaonZeroLong::KaonZeroLong();
+ theL = G4Lambda::Lambda();
+ theAntiL = G4AntiLambda::AntiLambda();
+ theSPlus = G4SigmaPlus::SigmaPlus();
+ theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
+ theSMinus = G4SigmaMinus::SigmaMinus();
+ theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
+ theS0 = G4SigmaZero::SigmaZero();
+ theAS0 = G4AntiSigmaZero::AntiSigmaZero();
+ theXiMinus = G4XiMinus::XiMinus();
+ theXi0 = G4XiZero::XiZero();
+ theAXiMinus = G4AntiXiMinus::AntiXiMinus();
+ theAXi0 = G4AntiXiZero::AntiXiZero();
+ theOmega = G4OmegaMinus::OmegaMinus();
+ theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
+ theD = G4Deuteron::Deuteron();
+ theT = G4Triton::Triton();
+ theA = G4Alpha::Alpha();
+ theHe3 = G4He3::He3();
+
+ hnXsc = new G4HadronNucleonXsc();
+}
+
+///////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4ComponentGGHadronNucleusXsc::~G4ComponentGGHadronNucleusXsc()
+{
+ if (hnXsc) delete hnXsc;
+}
+
+////////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetTotalIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, A);
+ delete aDP;
+
+ return fTotalXsc;
+}
+
+//////////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetTotalElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
+ delete aDP;
+
+ return fTotalXsc;
+}
+
+////////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetInelasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, A);
+ delete aDP;
+
+ return fInelasticXsc;
+}
+
+////////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetProductionIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, A);
+ delete aDP;
+
+ return fProductionXsc;
+}
+
+/////////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetInelasticElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
+ delete aDP;
+
+ return fInelasticXsc;
+}
+
+/////////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetProductionElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
+ delete aDP;
+
+ return fProductionXsc;
+}
+
+//////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetElasticElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
+ delete aDP;
+
+ return fElasticXsc;
+}
+
+///////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::GetElasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, A);
+ delete aDP;
+
+ return fElasticXsc;
+}
+
+////////////////////////////////////////////////////////////////
+
+G4double G4ComponentGGHadronNucleusXsc::ComputeQuasiElasticRatio(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A)
+{
+ G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
+ kinEnergy);
+ fTotalXsc = GetIsoCrossSection(aDP, Z, A);
+ delete aDP;
+ G4double ratio = 0.;
+
+ if(fInelasticXsc > 0.)
+ {
+ ratio = (fInelasticXsc - fProductionXsc)/fInelasticXsc;
+ if(ratio < 0.) ratio = 0.;
+ }
+ return ratio;
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////
+
+G4bool
+G4ComponentGGHadronNucleusXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
+ G4int Z, G4int /*A*/,
+ const G4Element*,
+ const G4Material*)
+{
+ G4bool applicable = false;
+ // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
+ G4double kineticEnergy = aDP->GetKineticEnergy();
+
+ const G4ParticleDefinition* theParticle = aDP->GetDefinition();
+
+ if ( ( kineticEnergy >= fLowerLimit &&
+ // Z > 1 && // >= He
+ ( theParticle == theAProton ||
+ theParticle == theGamma ||
+ theParticle == theSMinus ||
+ theParticle == theProton ||
+ theParticle == theNeutron ||
+ theParticle == thePiPlus ||
+ theParticle == thePiMinus ) &&
+ Z >= 1 && // >= H for kaons
+ (
+ theParticle == theKPlus ||
+ theParticle == theKMinus ||
+ theParticle == theK0L ||
+ theParticle == theK0S
+ ) ) ) applicable = true;
+
+ return applicable;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////
+//
+// Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
+// Glauber model with Gribov correction calculated in the dipole approximation on
+// light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model.
+// [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetIsoCrossSection(const G4DynamicParticle* aParticle,
+ G4int Z, G4int A,
+ const G4Isotope*,
+ const G4Element*,
+ const G4Material*)
+{
+ G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
+ G4double hpInXsc(0.), hnInXsc(0.);
+ G4double R = GetNucleusRadius(A);
+
+ G4int N = A - Z; // number of neutrons
+ if (N < 0) N = 0;
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ if( theParticle == theProton ||
+ theParticle == theNeutron ||
+ theParticle == thePiPlus ||
+ theParticle == thePiMinus )
+ {
+ // sigma = GetHadronNucleonXscNS(aParticle, A, Z);
+
+ sigma = Z*hnXsc->GetHadronNucleonXscNS(aParticle, theProton);
+
+ hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
+
+ sigma += N*hnXsc->GetHadronNucleonXscNS(aParticle, theNeutron);
+
+ hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
+
+ cofInelastic = 2.4;
+ cofTotal = 2.0;
+ }
+ else if( theParticle == theKPlus ||
+ theParticle == theKMinus ||
+ theParticle == theK0S ||
+ theParticle == theK0L )
+ {
+ // sigma = GetKaonNucleonXscVector(aParticle, A, Z);
+
+ sigma = Z*hnXsc->GetKaonNucleonXscGG(aParticle, theProton);
+
+ hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
+
+ sigma += N*hnXsc->GetKaonNucleonXscGG(aParticle, theNeutron);
+
+ hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
+
+ cofInelastic = 2.2;
+ cofTotal = 2.0;
+ R = 1.3*fermi;
+ R *= std::pow(G4double(A), 0.3333);
+ }
+ else
+ {
+ sigma = GetHadronNucleonXscNS(aParticle, A, Z);
+ cofInelastic = 2.2;
+ cofTotal = 2.0;
+ }
+ // cofInelastic = 2.0;
+
+ if( A > 1 )
+ {
+ nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
+ ratio = sigma/nucleusSquare;
+
+ xsection = nucleusSquare*std::log( 1. + ratio );
+
+ xsection *= GetParticleBarCorTot(theParticle, Z);
+
+ fTotalXsc = xsection;
+
+
+
+ fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
+
+ fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);
+
+ fElasticXsc = fTotalXsc - fInelasticXsc;
+
+ if(fElasticXsc < 0.) fElasticXsc = 0.;
+
+ G4double difratio = ratio/(1.+ratio);
+
+ fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
+
+
+ // sigma = GetHNinelasticXsc(aParticle, A, Z);
+
+ sigma = Z*hpInXsc + N*hnInXsc;
+
+ ratio = sigma/nucleusSquare;
+
+ fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
+
+ fProductionXsc *= GetParticleBarCorIn(theParticle, Z);
+
+ if (fElasticXsc < 0.) fElasticXsc = 0.;
+ }
+ else // H
+ {
+ fTotalXsc = sigma;
+ xsection = sigma;
+
+ fInelasticXsc = hnXsc->GetInelasticHadronNucleonXsc();
+
+ if ( theParticle != theAProton )
+ {
+ fElasticXsc = hnXsc->GetElasticHadronNucleonXsc();
+
+ // sigma = GetHNinelasticXsc(aParticle, A, Z);
+ // fInelasticXsc = sigma;
+ // fElasticXsc = fTotalXsc - fInelasticXsc;
+ }
+ else if( theParticle == theKPlus ||
+ theParticle == theKMinus ||
+ theParticle == theK0S ||
+ theParticle == theK0L )
+ {
+ fInelasticXsc = hpInXsc;
+ fElasticXsc = fTotalXsc - fInelasticXsc;
+ }
+ else
+ {
+ fInelasticXsc = hpInXsc;
+ fElasticXsc = fTotalXsc - fInelasticXsc;
+ }
+ if (fElasticXsc < 0.) fElasticXsc = 0.;
+
+ }
+ return xsection;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//
+// Return single-diffraction/inelastic cross-section ratio
+
+G4double G4ComponentGGHadronNucleusXsc::
+GetRatioSD(const G4DynamicParticle* aParticle, G4int A, G4int Z)
+{
+ G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
+ G4double R = GetNucleusRadius(A);
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ if( theParticle == theProton ||
+ theParticle == theNeutron ||
+ theParticle == thePiPlus ||
+ theParticle == thePiMinus )
+ {
+ sigma = GetHadronNucleonXscNS(aParticle, A, Z);
+ cofInelastic = 2.4;
+ cofTotal = 2.0;
+ }
+ else
+ {
+ sigma = GetHadronNucleonXscNS(aParticle, A, Z);
+ cofInelastic = 2.2;
+ cofTotal = 2.0;
+ }
+ nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
+ ratio = sigma/nucleusSquare;
+
+ fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
+
+ G4double difratio = ratio/(1.+ratio);
+
+ fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
+
+ if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
+ else ratio = 0.;
+
+ return ratio;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//
+// Return suasi-elastic/inelastic cross-section ratio
+
+G4double G4ComponentGGHadronNucleusXsc::
+GetRatioQE(const G4DynamicParticle* aParticle, G4int A, G4int Z)
+{
+ G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
+ G4double R = GetNucleusRadius(A);
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ if( theParticle == theProton ||
+ theParticle == theNeutron ||
+ theParticle == thePiPlus ||
+ theParticle == thePiMinus )
+ {
+ sigma = GetHadronNucleonXscNS(aParticle, A, Z);
+ cofInelastic = 2.4;
+ cofTotal = 2.0;
+ }
+ else
+ {
+ sigma = GetHadronNucleonXscNS(aParticle, A, Z);
+ cofInelastic = 2.2;
+ cofTotal = 2.0;
+ }
+ nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
+ ratio = sigma/nucleusSquare;
+
+ fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
+
+ sigma = GetHNinelasticXsc(aParticle, A, Z);
+ ratio = sigma/nucleusSquare;
+
+ fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
+
+ if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
+ else ratio = 0.;
+ if ( ratio < 0. ) ratio = 0.;
+
+ return ratio;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon Xsc according to differnt parametrisations:
+// [2] E. Levin, hep-ph/9710546
+// [3] U. Dersch, et al, hep-ex/9910052
+// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
+ const G4Element* anElement)
+{
+ G4int At = G4lrint(anElement->GetN()); // number of nucleons
+ G4int Zt = G4lrint(anElement->GetZ()); // number of protons
+
+ return GetHadronNucleonXsc(aParticle, At, Zt);
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon Xsc according to differnt parametrisations:
+// [2] E. Levin, hep-ph/9710546
+// [3] U. Dersch, et al, hep-ex/9910052
+// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
+ G4int At, G4int /*Zt*/)
+{
+ G4double xsection;
+
+ //G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
+
+ G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
+
+ G4double proj_mass = aParticle->GetMass();
+ G4double proj_momentum = aParticle->GetMomentum().mag();
+ G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
+
+ sMand /= GeV*GeV; // in GeV for parametrisation
+ proj_momentum /= GeV;
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ G4double aa = At;
+
+ if(theParticle == theGamma)
+ {
+ xsection = aa*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theNeutron) // as proton ???
+ {
+ xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theProton)
+ {
+ xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
+ // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
+ // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
+ }
+ else if(theParticle == theAProton)
+ {
+ xsection = aa*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == thePiPlus)
+ {
+ xsection = aa*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == thePiMinus)
+ {
+ // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
+ xsection = aa*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theKPlus)
+ {
+ xsection = aa*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theKMinus)
+ {
+ xsection = aa*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
+ }
+ else // as proton ???
+ {
+ xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
+ }
+ xsection *= millibarn;
+ return xsection;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
+// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
+ const G4Element* anElement)
+{
+ G4int At = G4lrint(anElement->GetN()); // number of nucleons
+ G4int Zt = G4lrint(anElement->GetZ()); // number of protons
+
+ return GetHadronNucleonXscPDG(aParticle, At, Zt);
+}
+
+
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
+// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
+// At = number of nucleons, Zt = number of protons
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
+ G4int At, G4int Zt)
+{
+ G4double xsection;
+
+ G4int Nt = At-Zt; // number of neutrons
+ if (Nt < 0) Nt = 0;
+
+ G4double zz = Zt;
+ G4double aa = At;
+ G4double nn = Nt;
+
+ G4double targ_mass = G4ParticleTable::GetParticleTable()->
+ GetIonTable()->GetIonMass(Zt, At);
+
+ targ_mass = 0.939*GeV; // ~mean neutron and proton ???
+
+ G4double proj_mass = aParticle->GetMass();
+ G4double proj_momentum = aParticle->GetMomentum().mag();
+
+ G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
+
+ sMand /= GeV*GeV; // in GeV for parametrisation
+
+ // General PDG fit constants
+
+ G4double s0 = 5.38*5.38; // in Gev^2
+ G4double eta1 = 0.458;
+ G4double eta2 = 0.458;
+ G4double B = 0.308;
+
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+
+ if(theParticle == theNeutron) // proton-neutron fit
+ {
+ xsection = zz*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
+ xsection += nn*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
+ }
+ else if(theParticle == theProton)
+ {
+
+ xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
+
+ xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theAProton)
+ {
+ xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
+
+ xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == thePiPlus)
+ {
+ xsection = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
+ + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == thePiMinus)
+ {
+ xsection = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
+ + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theKPlus || theParticle == theK0L )
+ {
+ xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
+ + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
+
+ xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
+ + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theKMinus || theParticle == theK0S )
+ {
+ xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
+ + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
+
+ xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
+ + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theSMinus)
+ {
+ xsection = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
+ - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theGamma) // modify later on
+ {
+ xsection = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
+ + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
+
+ }
+ else // as proton ???
+ {
+ xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
+
+ xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
+ }
+ xsection *= millibarn; // parametrised in mb
+ return xsection;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
+// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
+ const G4Element* anElement)
+{
+ G4int At = G4lrint(anElement->GetN()); // number of nucleons
+ G4int Zt = G4lrint(anElement->GetZ()); // number of protons
+
+ return GetHadronNucleonXscNS(aParticle, At, Zt);
+}
+
+
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
+// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
+ G4int At, G4int Zt)
+{
+ G4double xsection(0);
+ // G4double Delta; DHW 19 May 2011: variable set but not used
+ G4double A0, B0;
+ G4double hpXscv(0);
+ G4double hnXscv(0);
+
+ G4int Nt = At-Zt; // number of neutrons
+ if (Nt < 0) Nt = 0;
+
+ G4double aa = At;
+ G4double zz = Zt;
+ G4double nn = Nt;
+
+ G4double targ_mass = G4ParticleTable::GetParticleTable()->
+ GetIonTable()->GetIonMass(Zt, At);
+
+ targ_mass = 0.939*GeV; // ~mean neutron and proton ???
+
+ G4double proj_mass = aParticle->GetMass();
+ G4double proj_energy = aParticle->GetTotalEnergy();
+ G4double proj_momentum = aParticle->GetMomentum().mag();
+
+ G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
+
+ sMand /= GeV*GeV; // in GeV for parametrisation
+ proj_momentum /= GeV;
+ proj_energy /= GeV;
+ proj_mass /= GeV;
+
+ // General PDG fit constants
+
+ G4double s0 = 5.38*5.38; // in Gev^2
+ G4double eta1 = 0.458;
+ G4double eta2 = 0.458;
+ G4double B = 0.308;
+
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+
+ if(theParticle == theNeutron)
+ {
+ if( proj_momentum >= 373.)
+ {
+ return GetHadronNucleonXscPDG(aParticle,At,Zt);
+ }
+ else if( proj_momentum >= 10.)
+ // if( proj_momentum >= 2.)
+ {
+ // Delta = 1.; // DHW 19 May 2011: variable set but not used
+ // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
+
+ if(proj_momentum >= 10.)
+ {
+ B0 = 7.5;
+ A0 = 100. - B0*std::log(3.0e7);
+
+ xsection = A0 + B0*std::log(proj_energy) - 11
+ + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
+ 0.93827*0.93827,-0.165); // mb
+ }
+ xsection *= zz + nn;
+ }
+ else
+ {
+ // nn to be pp
+
+ if( proj_momentum < 0.73 )
+ {
+ hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
+ }
+ else if( proj_momentum < 1.05 )
+ {
+ hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
+ (std::log(proj_momentum/0.73));
+ }
+ else // if( proj_momentum < 10. )
+ {
+ hnXscv = 39.0+
+ 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
+ }
+ // pn to be np
+
+ if( proj_momentum < 0.8 )
+ {
+ hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
+ }
+ else if( proj_momentum < 1.4 )
+ {
+ hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
+ }
+ else // if( proj_momentum < 10. )
+ {
+ hpXscv = 33.3+
+ 20.8*(std::pow(proj_momentum,2.0)-1.35)/
+ (std::pow(proj_momentum,2.50)+0.95);
+ }
+ xsection = hpXscv*zz + hnXscv*nn;
+ }
+ }
+ else if(theParticle == theProton)
+ {
+ if( proj_momentum >= 373.)
+ {
+ return GetHadronNucleonXscPDG(aParticle,At,Zt);
+ }
+ else if( proj_momentum >= 10.)
+ // if( proj_momentum >= 2.)
+ {
+ // Delta = 1.; DHW 19 May 2011: variable set but not used
+ // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
+
+ if(proj_momentum >= 10.)
+ {
+ B0 = 7.5;
+ A0 = 100. - B0*std::log(3.0e7);
+
+ xsection = A0 + B0*std::log(proj_energy) - 11
+ + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
+ 0.93827*0.93827,-0.165); // mb
+ }
+ xsection *= zz + nn;
+ }
+ else
+ {
+ // pp
+
+ if( proj_momentum < 0.73 )
+ {
+ hpXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
+ }
+ else if( proj_momentum < 1.05 )
+ {
+ hpXscv = 23 + 40*(std::log(proj_momentum/0.73))*
+ (std::log(proj_momentum/0.73));
+ }
+ else // if( proj_momentum < 10. )
+ {
+ hpXscv = 39.0+
+ 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
+ }
+ // pn to be np
+
+ if( proj_momentum < 0.8 )
+ {
+ hnXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
+ }
+ else if( proj_momentum < 1.4 )
+ {
+ hnXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
+ }
+ else // if( proj_momentum < 10. )
+ {
+ hnXscv = 33.3+
+ 20.8*(std::pow(proj_momentum,2.0)-1.35)/
+ (std::pow(proj_momentum,2.50)+0.95);
+ }
+ xsection = hpXscv*zz + hnXscv*nn;
+ // xsection = hpXscv*(Zt + Nt);
+ // xsection = hnXscv*(Zt + Nt);
+ }
+ // xsection *= 0.95;
+ }
+ else if( theParticle == theAProton )
+ {
+ // xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ // + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
+
+ // xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ // + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
+
+ G4double logP = std::log(proj_momentum);
+
+ if( proj_momentum <= 1.0 )
+ {
+ xsection = zz*(65.55 + 53.84/(proj_momentum+1.e-6) );
+ }
+ else
+ {
+ xsection = zz*( 41.1 + 77.2*std::pow( proj_momentum, -0.68)
+ + 0.293*logP*logP - 1.82*logP );
+ }
+ if ( nn > 0.)
+ {
+ xsection += nn*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
+ }
+ else // H
+ {
+ fInelasticXsc = 38.0 + 38.0*std::pow( proj_momentum, -0.96)
+ - 0.169*logP*logP;
+ fInelasticXsc *= millibarn;
+ }
+ }
+ else if( theParticle == thePiPlus )
+ {
+ if(proj_momentum < 0.4)
+ {
+ G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
+ hpXscv = Ex3+20.0;
+ }
+ else if( proj_momentum < 1.15 )
+ {
+ G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
+ hpXscv = Ex4+14.0;
+ }
+ else if(proj_momentum < 3.5)
+ {
+ G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
+ hpXscv = Ex1+Ex2+27.5;
+ }
+ else // if(proj_momentum > 3.5) // mb
+ {
+ hpXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
+ }
+ // pi+n = pi-p??
+
+ if(proj_momentum < 0.37)
+ {
+ hnXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
+ }
+ else if(proj_momentum<0.65)
+ {
+ hnXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
+ }
+ else if(proj_momentum<1.3)
+ {
+ hnXscv = 36.1+
+ 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
+ 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
+ }
+ else if(proj_momentum<3.0)
+ {
+ hnXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
+ 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
+ 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
+ }
+ else // mb
+ {
+ hnXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
+ }
+ xsection = hpXscv*zz + hnXscv*nn;
+ }
+ else if(theParticle == thePiMinus)
+ {
+ // pi-n = pi+p??
+
+ if(proj_momentum < 0.4)
+ {
+ G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
+ hnXscv = Ex3+20.0;
+ }
+ else if(proj_momentum < 1.15)
+ {
+ G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
+ hnXscv = Ex4+14.0;
+ }
+ else if(proj_momentum < 3.5)
+ {
+ G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
+ hnXscv = Ex1+Ex2+27.5;
+ }
+ else // if(proj_momentum > 3.5) // mb
+ {
+ hnXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
+ }
+ // pi-p
+
+ if(proj_momentum < 0.37)
+ {
+ hpXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
+ }
+ else if(proj_momentum<0.65)
+ {
+ hpXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
+ }
+ else if(proj_momentum<1.3)
+ {
+ hpXscv = 36.1+
+ 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
+ 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
+ }
+ else if(proj_momentum<3.0)
+ {
+ hpXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
+ 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
+ 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
+ }
+ else // mb
+ {
+ hpXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
+ }
+ xsection = hpXscv*zz + hnXscv*nn;
+ }
+ else if(theParticle == theKPlus)
+ {
+ xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
+ + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
+
+ xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
+ + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theKMinus)
+ {
+ xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
+ + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
+
+ xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
+ + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theSMinus)
+ {
+ xsection = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
+ - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theGamma) // modify later on
+ {
+ xsection = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
+ + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
+
+ }
+ else // as proton ???
+ {
+ xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
+
+ xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
+ }
+ xsection *= millibarn; // parametrised in mb
+ return xsection;
+}
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector(const G4DynamicParticle* aParticle,
+ G4int At, G4int Zt)
+{
+ G4double Tkin, logTkin, xsc, xscP, xscN;
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ G4int Nt = At-Zt; // number of neutrons
+ if (Nt < 0) Nt = 0;
+
+ Tkin = aParticle->GetKineticEnergy(); // Tkin in MeV
+
+ if( Tkin > 70*GeV ) return GetHadronNucleonXscPDG(aParticle,At,Zt);
+
+ logTkin = std::log(Tkin); // Tkin in MeV!!!
+
+ if( theParticle == theKPlus )
+ {
+ xscP = hnXsc->GetKpProtonTotXscVector(logTkin);
+ xscN = hnXsc->GetKpNeutronTotXscVector(logTkin);
+ }
+ else if( theParticle == theKMinus )
+ {
+ xscP = hnXsc->GetKmProtonTotXscVector(logTkin);
+ xscN = hnXsc->GetKmNeutronTotXscVector(logTkin);
+ }
+ else // K-zero as half of K+ and K-
+ {
+ xscP = (hnXsc->GetKpProtonTotXscVector(logTkin)+hnXsc->GetKmProtonTotXscVector(logTkin))*0.5;
+ xscN = (hnXsc->GetKpNeutronTotXscVector(logTkin)+hnXsc->GetKmNeutronTotXscVector(logTkin))*0.5;
+ }
+ xsc = xscP*Zt + xscN*Nt;
+ return xsc;
+}
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon inelastic cross-section based on proper parametrisation
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
+ const G4Element* anElement)
+{
+ G4int At = G4lrint(anElement->GetN()); // number of nucleons
+ G4int Zt = G4lrint(anElement->GetZ()); // number of protons
+
+ return GetHNinelasticXsc(aParticle, At, Zt);
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
+ G4int At, G4int Zt)
+{
+ const G4ParticleDefinition* hadron = aParticle->GetDefinition();
+ G4double sumInelastic;
+ G4int Nt = At - Zt;
+ if(Nt < 0) Nt = 0;
+
+ if( hadron == theKPlus )
+ {
+ sumInelastic = GetHNinelasticXscVU(aParticle, At, Zt);
+ }
+ else
+ {
+ //sumInelastic = Zt*GetHadronNucleonXscMK(aParticle, theProton);
+ // sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);
+ sumInelastic = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
+ sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);
+ }
+ return sumInelastic;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
+ G4int At, G4int Zt)
+{
+ G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
+ G4int absPDGcode = std::abs(PDGcode);
+
+ G4double Elab = aParticle->GetTotalEnergy();
+ // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
+ G4double Plab = aParticle->GetMomentum().mag();
+ // std::sqrt(Elab * Elab - 0.88);
+
+ Elab /= GeV;
+ Plab /= GeV;
+
+ G4double LogPlab = std::log( Plab );
+ G4double sqrLogPlab = LogPlab * LogPlab;
+
+ //G4cout<<"Plab = "<<Plab<<G4endl;
+
+ G4double NumberOfTargetProtons = G4double(Zt);
+ G4double NumberOfTargetNucleons = G4double(At);
+ G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
+
+ if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;
+
+ G4double Xtotal, Xelastic, Xinelastic;
+
+ if( absPDGcode > 1000 ) //------Projectile is baryon --------
+ {
+ G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
+ 0.522*sqrLogPlab - 4.51*LogPlab;
+
+ G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
+ 0.513*sqrLogPlab - 4.27*LogPlab;
+
+ G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
+ 0.169*sqrLogPlab - 1.85*LogPlab;
+
+ G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
+ 0.169*sqrLogPlab - 1.85*LogPlab;
+
+ Xtotal = (NumberOfTargetProtons * XtotPP +
+ NumberOfTargetNeutrons * XtotPN);
+
+ Xelastic = (NumberOfTargetProtons * XelPP +
+ NumberOfTargetNeutrons * XelPN);
+ }
+ else if( PDGcode == 211 ) //------Projectile is PionPlus -------
+ {
+ G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) +
+ 0.19 *sqrLogPlab - 0.0 *LogPlab;
+
+ G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) +
+ 0.456*sqrLogPlab - 4.03*LogPlab;
+
+ G4double XelPiP = 0.0 + 11.4*std::pow(Plab,-0.40) +
+ 0.079*sqrLogPlab - 0.0 *LogPlab;
+
+ G4double XelPiN = 1.76 + 11.2*std::pow(Plab,-0.64) +
+ 0.043*sqrLogPlab - 0.0 *LogPlab;
+
+ Xtotal = ( NumberOfTargetProtons * XtotPiP +
+ NumberOfTargetNeutrons * XtotPiN );
+
+ Xelastic = ( NumberOfTargetProtons * XelPiP +
+ NumberOfTargetNeutrons * XelPiN );
+ }
+ else if( PDGcode == -211 ) //------Projectile is PionMinus -------
+ {
+ G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) +
+ 0.456*sqrLogPlab - 4.03*LogPlab;
+
+ G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) +
+ 0.19 *sqrLogPlab - 0.0 *LogPlab;
+
+ G4double XelPiP = 1.76 + 11.2*std::pow(Plab,-0.64) +
+ 0.043*sqrLogPlab - 0.0 *LogPlab;
+
+ G4double XelPiN = 0.0 + 11.4*std::pow(Plab,-0.40) +
+ 0.079*sqrLogPlab - 0.0 *LogPlab;
+
+ Xtotal = ( NumberOfTargetProtons * XtotPiP +
+ NumberOfTargetNeutrons * XtotPiN );
+
+ Xelastic = ( NumberOfTargetProtons * XelPiP +
+ NumberOfTargetNeutrons * XelPiN );
+ }
+ else if( PDGcode == 111 ) //------Projectile is PionZero -------
+ {
+ G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) +
+ 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
+ 33.0 + 14.0 *std::pow(Plab,-1.36) +
+ 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
+
+ G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) +
+ 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
+ 16.4 + 19.3 *std::pow(Plab,-0.42) +
+ 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
+
+ G4double XelPiP =( 0.0 + 11.4*std::pow(Plab,-0.40) +
+ 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
+ 1.76 + 11.2*std::pow(Plab,-0.64) +
+ 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
+
+ G4double XelPiN =( 1.76 + 11.2*std::pow(Plab,-0.64) +
+ 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
+ 0.0 + 11.4*std::pow(Plab,-0.40) +
+ 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
+
+ Xtotal = ( NumberOfTargetProtons * XtotPiP +
+ NumberOfTargetNeutrons * XtotPiN );
+
+ Xelastic = ( NumberOfTargetProtons * XelPiP +
+ NumberOfTargetNeutrons * XelPiN );
+ }
+ else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
+ {
+ G4double XtotKP = 18.1 + 0. *std::pow(Plab, 0. ) +
+ 0.26 *sqrLogPlab - 1.0 *LogPlab;
+ G4double XtotKN = 18.7 + 0. *std::pow(Plab, 0. ) +
+ 0.21 *sqrLogPlab - 0.89*LogPlab;
+
+ G4double XelKP = 5.0 + 8.1*std::pow(Plab,-1.8 ) +
+ 0.16 *sqrLogPlab - 1.3 *LogPlab;
+
+ G4double XelKN = 7.3 + 0. *std::pow(Plab,-0. ) +
+ 0.29 *sqrLogPlab - 2.4 *LogPlab;
+
+ Xtotal = ( NumberOfTargetProtons * XtotKP +
+ NumberOfTargetNeutrons * XtotKN );
+
+ Xelastic = ( NumberOfTargetProtons * XelKP +
+ NumberOfTargetNeutrons * XelKN );
+ }
+ else if( PDGcode ==-321 ) //------Projectile is KaonMinus ------
+ {
+ G4double XtotKP = 32.1 + 0. *std::pow(Plab, 0. ) +
+ 0.66 *sqrLogPlab - 5.6 *LogPlab;
+ G4double XtotKN = 25.2 + 0. *std::pow(Plab, 0. ) +
+ 0.38 *sqrLogPlab - 2.9 *LogPlab;
+
+ G4double XelKP = 7.3 + 0. *std::pow(Plab,-0. ) +
+ 0.29 *sqrLogPlab - 2.4 *LogPlab;
+
+ G4double XelKN = 5.0 + 8.1*std::pow(Plab,-1.8 ) +
+ 0.16 *sqrLogPlab - 1.3 *LogPlab;
+
+ Xtotal = ( NumberOfTargetProtons * XtotKP +
+ NumberOfTargetNeutrons * XtotKN );
+
+ Xelastic = ( NumberOfTargetProtons * XelKP +
+ NumberOfTargetNeutrons * XelKN );
+ }
+ else if( PDGcode == 311 ) //------Projectile is KaonZero ------
+ {
+ G4double XtotKP = ( 18.1 + 0. *std::pow(Plab, 0. ) +
+ 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
+ 32.1 + 0. *std::pow(Plab, 0. ) +
+ 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
+
+ G4double XtotKN = ( 18.7 + 0. *std::pow(Plab, 0. ) +
+ 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
+ 25.2 + 0. *std::pow(Plab, 0. ) +
+ 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
+
+ G4double XelKP = ( 5.0 + 8.1*std::pow(Plab,-1.8 )
+ + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
+ 7.3 + 0. *std::pow(Plab,-0. ) +
+ 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
+
+ G4double XelKN = ( 7.3 + 0. *std::pow(Plab,-0. ) +
+ 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
+ 5.0 + 8.1*std::pow(Plab,-1.8 ) +
+ 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
+
+ Xtotal = ( NumberOfTargetProtons * XtotKP +
+ NumberOfTargetNeutrons * XtotKN );
+
+ Xelastic = ( NumberOfTargetProtons * XelKP +
+ NumberOfTargetNeutrons * XelKN );
+ }
+ else //------Projectile is undefined, Nucleon assumed
+ {
+ G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
+ 0.522*sqrLogPlab - 4.51*LogPlab;
+
+ G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
+ 0.513*sqrLogPlab - 4.27*LogPlab;
+
+ G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
+ 0.169*sqrLogPlab - 1.85*LogPlab;
+ G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
+ 0.169*sqrLogPlab - 1.85*LogPlab;
+
+ Xtotal = ( NumberOfTargetProtons * XtotPP +
+ NumberOfTargetNeutrons * XtotPN );
+
+ Xelastic = ( NumberOfTargetProtons * XelPP +
+ NumberOfTargetNeutrons * XelPN );
+ }
+ Xinelastic = Xtotal - Xelastic;
+
+ if( Xinelastic < 0.) Xinelastic = 0.;
+
+ return Xinelastic*= millibarn;
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetNucleusRadius(const G4DynamicParticle* ,
+ const G4Element* anElement)
+{
+ G4int At = G4lrint(anElement->GetN());
+ G4double oneThird = 1.0/3.0;
+ G4double cubicrAt = std::pow(G4double(At), oneThird);
+
+ G4double R; // = fRadiusConst*cubicrAt;
+ /*
+ G4double tmp = std::pow( cubicrAt-1., 3.);
+ tmp += At;
+ tmp *= 0.5;
+
+ if (At > 20.) // 20.
+ {
+ R = fRadiusConst*std::pow (tmp, oneThird);
+ }
+ else
+ {
+ R = fRadiusConst*cubicrAt;
+ }
+ */
+
+ R = fRadiusConst*cubicrAt;
+
+ G4double meanA = 21.;
+
+ G4double tauA1 = 40.;
+ G4double tauA2 = 10.;
+ G4double tauA3 = 5.;
+
+ G4double a1 = 0.85;
+ G4double b1 = 1. - a1;
+
+ G4double b2 = 0.3;
+ G4double b3 = 4.;
+
+ if (At > 20) // 20.
+ {
+ R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
+ }
+ else if (At > 3)
+ {
+ R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
+ }
+ else
+ {
+ R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
+ }
+ return R;
+
+}
+////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double
+G4ComponentGGHadronNucleusXsc::GetNucleusRadius(G4int At)
+{
+ G4double oneThird = 1.0/3.0;
+ G4double cubicrAt = std::pow(G4double(At), oneThird);
+
+ G4double R; // = fRadiusConst*cubicrAt;
+
+ /*
+ G4double tmp = std::pow( cubicrAt-1., 3.);
+ tmp += At;
+ tmp *= 0.5;
+
+ if (At > 20.)
+ {
+ R = fRadiusConst*std::pow (tmp, oneThird);
+ }
+ else
+ {
+ R = fRadiusConst*cubicrAt;
+ }
+ */
+
+ R = fRadiusConst*cubicrAt;
+
+ G4double meanA = 20.;
+ G4double tauA = 20.;
+
+ if (At > 20) // 20.
+ {
+ R *= ( 0.8 + 0.2*std::exp( -(G4double(At) - meanA)/tauA) );
+ }
+ else
+ {
+ R *= ( 1.0 + 0.1*( 1. - std::exp( (G4double(At) - meanA)/tauA) ) );
+ }
+
+ return R;
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double G4ComponentGGHadronNucleusXsc::CalculateEcmValue( const G4double mp ,
+ const G4double mt ,
+ const G4double Plab )
+{
+ G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
+ G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
+ // G4double Pcm = Plab * mt / Ecm;
+ // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
+
+ return Ecm ; // KEcm;
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double G4ComponentGGHadronNucleusXsc::CalcMandelstamS( const G4double mp ,
+ const G4double mt ,
+ const G4double Plab )
+{
+ G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
+ G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
+
+ return sMand;
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+void G4ComponentGGHadronNucleusXsc::CrossSectionDescription(std::ostream& outFile) const
+{
+ outFile << "G4ComponentGGHadronNucleusXsc calculates total, inelastic and\n"
+ << "elastic cross sections for hadron-nucleus cross sections using\n"
+ << "the Glauber model with Gribov corrections. It is valid for all\n"
+ << "targets except hydrogen, and for incident p, pbar, n, sigma-,\n"
+ << "pi+, pi-, K+, K- and gammas with energies above 3 GeV. This is\n"
+ << "a cross section component which is to be used to build a cross\n"
+ << "data set.\n";
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// Correction arrays for GG <-> Bar changea at ~ 90 GeV
+
+const G4double G4ComponentGGHadronNucleusXsc::fNeutronBarCorrectionTot[93] = {
+
+ 1.0, 1.0, 1.42517e+00, // 1.118517e+00,
+1.082002e+00, 1.116171e+00, 1.078747e+00, 1.061315e+00,
+1.058205e+00, 1.082663e+00, 1.068500e+00, 1.076912e+00, 1.083475e+00, 1.079117e+00,
+1.071856e+00, 1.071990e+00, 1.073774e+00, 1.079356e+00, 1.081314e+00, 1.082056e+00,
+1.090772e+00, 1.096776e+00, 1.095828e+00, 1.097678e+00, 1.099157e+00, 1.103677e+00,
+1.105132e+00, 1.109806e+00, 1.110816e+00, 1.117378e+00, 1.115165e+00, 1.115710e+00,
+1.111855e+00, 1.110482e+00, 1.110112e+00, 1.106676e+00, 1.108706e+00, 1.105549e+00,
+1.106318e+00, 1.106242e+00, 1.107672e+00, 1.107342e+00, 1.108119e+00, 1.106655e+00,
+1.102588e+00, 1.096657e+00, 1.092920e+00, 1.086629e+00, 1.083592e+00, 1.076030e+00,
+1.083777e+00, 1.089460e+00, 1.086545e+00, 1.079924e+00, 1.082218e+00, 1.077798e+00,
+1.077062e+00, 1.072825e+00, 1.072241e+00, 1.072104e+00, 1.072490e+00, 1.069829e+00,
+1.070398e+00, 1.065458e+00, 1.064968e+00, 1.060524e+00, 1.060048e+00, 1.057620e+00,
+1.056428e+00, 1.055366e+00, 1.055017e+00, 1.052304e+00, 1.051767e+00, 1.049728e+00,
+1.048745e+00, 1.047399e+00, 1.045876e+00, 1.042972e+00, 1.041824e+00, 1.039993e+00,
+1.039021e+00, 1.036627e+00, 1.034176e+00, 1.032526e+00, 1.033633e+00, 1.036107e+00,
+1.037803e+00, 1.031266e+00, 1.032991e+00, 1.033284e+00, 1.035015e+00, 1.033945e+00,
+1.037075e+00, 1.034721e+00
+
+};
+
+const G4double G4ComponentGGHadronNucleusXsc::fNeutronBarCorrectionIn[93] = {
+
+1.0, 1.0, 1.167421e+00, 1.156250e+00, 1.205364e+00, 1.154225e+00, 1.120391e+00, // 6
+1.124632e+00, 1.129460e+00, 1.107863e+00, 1.102152e+00, 1.104593e+00, 1.100285e+00, // 12
+1.098450e+00, 1.092677e+00, 1.101124e+00, 1.106461e+00, 1.115049e+00, 1.123903e+00, // 18
+1.126661e+00, 1.131259e+00, 1.133949e+00, 1.134185e+00, 1.133767e+00, 1.132813e+00, // 24
+1.131515e+00, 1.144338e+00, // 1.130338e+00,
+1.134171e+00, 1.139206e+00, 1.148474e+00, // 1.141474e+00,
+1.142189e+00,
+1.140725e+00, 1.140100e+00, 1.139848e+00, 1.137674e+00, 1.138645e+00, 1.136339e+00,
+1.136439e+00, 1.135946e+00, 1.136431e+00, 1.135702e+00, 1.135703e+00, 1.134113e+00,
+1.131935e+00, 1.128381e+00, 1.126373e+00, 1.122453e+00, 1.120908e+00, 1.115953e+00,
+1.115947e+00, 1.114426e+00, 1.111749e+00, 1.106207e+00, 1.107494e+00, 1.103622e+00,
+1.102576e+00, 1.098816e+00, 1.097889e+00, 1.097306e+00, 1.097130e+00, 1.094578e+00,
+1.094552e+00, 1.090222e+00, 1.089358e+00, 1.085409e+00, 1.084560e+00, 1.082182e+00,
+1.080773e+00, 1.079464e+00, 1.078724e+00, 1.076121e+00, 1.075235e+00, 1.073159e+00,
+1.071920e+00, 1.070395e+00, 1.069503e+00, 1.067525e+00, 1.066919e+00, 1.065779e+00,
+1.065319e+00, 1.063730e+00, 1.062092e+00, 1.061085e+00, 1.059908e+00, 1.059815e+00,
+1.059109e+00, 1.051920e+00, 1.051258e+00, 1.049473e+00, 1.048823e+00, 1.045984e+00,
+1.046435e+00, 1.042614e+00
+
+};
+
+const G4double G4ComponentGGHadronNucleusXsc::fProtonBarCorrectionTot[93] = {
+
+1.0, 1.0,
+1.118515e+00, 1.082000e+00, 1.116169e+00, 1.078745e+00, 1.061313e+00, 1.058203e+00,
+1.082661e+00, 1.068498e+00, 1.076910e+00, 1.083474e+00, 1.079115e+00, 1.071854e+00,
+1.071988e+00, 1.073772e+00, 1.079355e+00, 1.081312e+00, 1.082054e+00, 1.090770e+00,
+1.096774e+00, 1.095827e+00, 1.097677e+00, 1.099156e+00, 1.103676e+00, 1.105130e+00,
+1.109805e+00, 1.110814e+00, 1.117377e+00, 1.115163e+00, 1.115708e+00, 1.111853e+00,
+1.110480e+00, 1.110111e+00, 1.106674e+00, 1.108705e+00, 1.105548e+00, 1.106317e+00,
+1.106241e+00, 1.107671e+00, 1.107341e+00, 1.108118e+00, 1.106654e+00, 1.102586e+00,
+1.096655e+00, 1.092918e+00, 1.086628e+00, 1.083590e+00, 1.076028e+00, 1.083776e+00,
+1.089458e+00, 1.086543e+00, 1.079923e+00, 1.082216e+00, 1.077797e+00, 1.077061e+00,
+1.072824e+00, 1.072239e+00, 1.072103e+00, 1.072488e+00, 1.069828e+00, 1.070396e+00,
+1.065456e+00, 1.064966e+00, 1.060523e+00, 1.060047e+00, 1.057618e+00, 1.056427e+00,
+1.055365e+00, 1.055016e+00, 1.052303e+00, 1.051766e+00, 1.049727e+00, 1.048743e+00,
+1.047397e+00, 1.045875e+00, 1.042971e+00, 1.041823e+00, 1.039992e+00, 1.039019e+00,
+1.036626e+00, 1.034175e+00, 1.032525e+00, 1.033632e+00, 1.036106e+00, 1.037802e+00,
+1.031265e+00, 1.032990e+00, 1.033283e+00, 1.035014e+00, 1.033944e+00, 1.037074e+00,
+1.034720e+00
+
+};
+
+const G4double G4ComponentGGHadronNucleusXsc::fProtonBarCorrectionIn[93] = {
+
+1.0, 1.0,
+1.147419e+00, // 1.167419e+00,
+1.156248e+00, 1.205362e+00, 1.154224e+00, 1.120390e+00, 1.124630e+00, // 7
+1.129459e+00, 1.107861e+00, 1.102151e+00, 1.104591e+00, 1.100284e+00, 1.098449e+00, // 13
+1.092675e+00, 1.101122e+00, 1.106460e+00, 1.115048e+00, 1.123902e+00, 1.126659e+00, // 19
+1.131258e+00, 1.133948e+00, 1.134183e+00, 1.133766e+00, 1.132812e+00, 1.131514e+00, // 25
+1.140337e+00, // 1.130337e+00,
+
+1.134170e+00, 1.139205e+00, 1.151472e+00, // 1.141472e+00,
+1.142188e+00, 1.140724e+00,
+1.140099e+00, 1.139847e+00, 1.137672e+00, 1.138644e+00, 1.136338e+00, 1.136438e+00,
+1.135945e+00, 1.136429e+00, 1.135701e+00, 1.135702e+00, 1.134112e+00, 1.131934e+00,
+1.128380e+00, 1.126371e+00, 1.122452e+00, 1.120907e+00, 1.115952e+00, 1.115946e+00,
+1.114425e+00, 1.111748e+00, 1.106205e+00, 1.107493e+00, 1.103621e+00, 1.102575e+00,
+1.098815e+00, 1.097888e+00, 1.097305e+00, 1.097129e+00, 1.094577e+00, 1.094551e+00,
+1.090221e+00, 1.089357e+00, 1.085408e+00, 1.084559e+00, 1.082181e+00, 1.080772e+00,
+1.079463e+00, 1.078723e+00, 1.076120e+00, 1.075234e+00, 1.073158e+00, 1.071919e+00,
+1.070394e+00, 1.069502e+00, 1.067524e+00, 1.066918e+00, 1.065778e+00, 1.065318e+00,
+1.063729e+00, 1.062091e+00, 1.061084e+00, 1.059907e+00, 1.059814e+00, 1.059108e+00,
+1.051919e+00, 1.051257e+00, 1.049472e+00, 1.048822e+00, 1.045983e+00, 1.046434e+00,
+1.042613e+00
+
+};
+
+
+const G4double G4ComponentGGHadronNucleusXsc::fPionPlusBarCorrectionTot[93] = {
+
+1.0, 1.0,
+1.075927e+00, 1.074407e+00, 1.126098e+00, 1.100127e+00, 1.089742e+00, 1.083536e+00,
+1.089988e+00, 1.103566e+00, 1.096922e+00, 1.126573e+00, 1.132734e+00, 1.136512e+00,
+1.136629e+00, 1.133086e+00, 1.132428e+00, 1.129299e+00, 1.125622e+00, 1.126992e+00,
+1.127840e+00, 1.162670e+00, 1.160392e+00, 1.157864e+00, 1.157227e+00, 1.154627e+00,
+1.192555e+00, 1.197243e+00, 1.197911e+00, 1.200326e+00, 1.220053e+00, 1.215019e+00,
+1.211703e+00, 1.209080e+00, 1.204248e+00, 1.203328e+00, 1.198671e+00, 1.196840e+00,
+1.194392e+00, 1.193037e+00, 1.190408e+00, 1.188583e+00, 1.206127e+00, 1.210028e+00,
+1.206434e+00, 1.204456e+00, 1.200547e+00, 1.199058e+00, 1.200174e+00, 1.200276e+00,
+1.198912e+00, 1.213048e+00, 1.207160e+00, 1.208020e+00, 1.203814e+00, 1.202380e+00,
+1.198306e+00, 1.197002e+00, 1.196027e+00, 1.195449e+00, 1.192563e+00, 1.192135e+00,
+1.187556e+00, 1.186308e+00, 1.182124e+00, 1.180900e+00, 1.178224e+00, 1.176471e+00,
+1.174811e+00, 1.173702e+00, 1.170827e+00, 1.169581e+00, 1.167205e+00, 1.165626e+00,
+1.180244e+00, 1.177626e+00, 1.175121e+00, 1.173903e+00, 1.172192e+00, 1.171128e+00,
+1.168997e+00, 1.166826e+00, 1.164130e+00, 1.165412e+00, 1.165504e+00, 1.165020e+00,
+1.158462e+00, 1.158014e+00, 1.156519e+00, 1.156081e+00, 1.153602e+00, 1.154190e+00,
+1.152974e+00
+
+};
+
+const G4double G4ComponentGGHadronNucleusXsc::fPionPlusBarCorrectionIn[93] = {
+
+1.0, 1.0,
+1.140246e+00, 1.097872e+00, 1.104301e+00, 1.068722e+00, 1.056495e+00, 1.062622e+00, // 7
+1.047987e+00, 1.037032e+00, 1.035686e+00, 1.042870e+00, 1.052222e+00, 1.075100e+00, // 13
+1.084480e+00, 1.078286e+00, 1.081488e+00, 1.089713e+00, 1.099105e+00, 1.098003e+00, // 19
+1.102175e+00, 1.117707e+00, 1.121734e+00, 1.125229e+00, 1.126457e+00, 1.128905e+00, // 25
+1.163312e+00, 1.126263e+00, 1.126459e+00, 1.135191e+00, 1.116986e+00, 1.117184e+00, // 31
+1.117037e+00, 1.116777e+00, 1.115858e+00, 1.115745e+00, 1.114489e+00, 1.113993e+00, // 37
+1.113226e+00, 1.112818e+00, 1.111890e+00, 1.111238e+00, 1.111209e+00, 1.111775e+00, // 43
+1.110256e+00, 1.109414e+00, 1.107647e+00, 1.106980e+00, 1.106096e+00, 1.107331e+00, // 49
+1.107849e+00, 1.106407e+00, 1.103426e+00, 1.103896e+00, 1.101756e+00, 1.101031e+00, // 55
+1.098915e+00, 1.098260e+00, 1.097768e+00, 1.097487e+00, 1.095964e+00, 1.095773e+00, // 61
+1.093348e+00, 1.092687e+00, 1.090465e+00, 1.089821e+00, 1.088394e+00, 1.087462e+00, // 67
+1.086571e+00, 1.085997e+00, 1.084451e+00, 1.083798e+00, 1.082513e+00, 1.081670e+00, // 73
+1.080735e+00, 1.075659e+00, 1.074341e+00, 1.073689e+00, 1.072787e+00, 1.072237e+00, // 79
+1.071107e+00, 1.069955e+00, 1.074856e+00, 1.065873e+00, 1.065938e+00, 1.065694e+00,
+1.062192e+00, 1.061967e+00, 1.061180e+00, 1.060960e+00, 1.059646e+00, 1.059975e+00,
+1.059658e+00
+
+};
+
+
+const G4double G4ComponentGGHadronNucleusXsc::fPionMinusBarCorrectionTot[93] = {
+
+1.0, 1.0,
+1.3956e+00, 1.077959e+00, 1.129145e+00, 1.102088e+00, 1.089765e+00, 1.083542e+00, // 7
+1.089995e+00, 1.104895e+00, 1.097154e+00, 1.127663e+00, 1.133063e+00, 1.137425e+00, // 13
+1.136724e+00, 1.133859e+00, 1.132498e+00, 1.130276e+00, 1.127896e+00, 1.127656e+00, // 19
+1.127905e+00, 1.164210e+00, 1.162259e+00, 1.160075e+00, 1.158978e+00, 1.156649e+00, // 25
+1.194157e+00, 1.199177e+00, 1.198983e+00, 1.202325e+00, 1.221967e+00, 1.217548e+00,
+1.214389e+00, 1.211760e+00, 1.207335e+00, 1.206081e+00, 1.201766e+00, 1.199779e+00,
+1.197283e+00, 1.195706e+00, 1.193071e+00, 1.191115e+00, 1.208838e+00, 1.212681e+00,
+1.209235e+00, 1.207163e+00, 1.203451e+00, 1.201807e+00, 1.203283e+00, 1.203388e+00,
+1.202244e+00, 1.216509e+00, 1.211066e+00, 1.211504e+00, 1.207539e+00, 1.205991e+00,
+1.202143e+00, 1.200724e+00, 1.199595e+00, 1.198815e+00, 1.196025e+00, 1.195390e+00,
+1.191137e+00, 1.189791e+00, 1.185888e+00, 1.184575e+00, 1.181996e+00, 1.180229e+00,
+1.178545e+00, 1.177355e+00, 1.174616e+00, 1.173312e+00, 1.171016e+00, 1.169424e+00,
+1.184120e+00, 1.181478e+00, 1.179085e+00, 1.177817e+00, 1.176124e+00, 1.175003e+00,
+1.172947e+00, 1.170858e+00, 1.168170e+00, 1.169397e+00, 1.169304e+00, 1.168706e+00,
+1.162774e+00, 1.162217e+00, 1.160740e+00, 1.160196e+00, 1.157857e+00, 1.158220e+00,
+1.157267e+00
+};
+
+
+const G4double G4ComponentGGHadronNucleusXsc::fPionMinusBarCorrectionIn[93] = {
+
+1.0, 1.0,
+1.463e+00, 1.100898e+00, 1.106773e+00, 1.070289e+00, 1.040514e+00, 1.062628e+00, // 7
+1.047992e+00, 1.038041e+00, 1.035862e+00, 1.043679e+00, 1.052466e+00, 1.065780e+00, // 13
+1.070551e+00, 1.078869e+00, 1.081541e+00, 1.090455e+00, 1.100847e+00, 1.098511e+00, // 19
+1.102226e+00, 1.118865e+00, 1.123143e+00, 1.126904e+00, 1.127785e+00, 1.130444e+00, // 25
+1.148502e+00, 1.127678e+00, 1.127244e+00, 1.123634e+00, 1.118347e+00, 1.118988e+00,
+1.118957e+00, 1.118696e+00, 1.118074e+00, 1.117722e+00, 1.116717e+00, 1.116111e+00,
+1.115311e+00, 1.114745e+00, 1.113814e+00, 1.113069e+00, 1.113141e+00, 1.113660e+00,
+1.112249e+00, 1.111343e+00, 1.109718e+00, 1.108942e+00, 1.108310e+00, 1.109549e+00,
+1.110227e+00, 1.108846e+00, 1.106183e+00, 1.106354e+00, 1.104388e+00, 1.103583e+00,
+1.101632e+00, 1.100896e+00, 1.100296e+00, 1.099873e+00, 1.098420e+00, 1.098082e+00,
+1.095892e+00, 1.095162e+00, 1.093144e+00, 1.092438e+00, 1.091083e+00, 1.090142e+00,
+1.089236e+00, 1.088604e+00, 1.087159e+00, 1.086465e+00, 1.085239e+00, 1.084388e+00,
+1.083473e+00, 1.078373e+00, 1.077136e+00, 1.076450e+00, 1.075561e+00, 1.074973e+00,
+1.073898e+00, 1.072806e+00, 1.067706e+00, 1.068684e+00, 1.068618e+00, 1.068294e+00,
+1.065241e+00, 1.064939e+00, 1.064166e+00, 1.063872e+00, 1.062659e+00, 1.062828e+00,
+1.062699e+00
+
+};
+
+
+//
+//
+///////////////////////////////////////////////////////////////////////////////////////
diff --git a/Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.hh b/Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.hh
new file mode 100644
index 000000000..e633db1bf
--- /dev/null
+++ b/Geant4/G4processes/srcnew/cross_sections/G4ComponentGGHadronNucleusXsc.hh
@@ -0,0 +1,278 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+// Calculation of the total, elastic and inelastic cross-sections
+// based on parametrisations of (proton, pion, kaon, photon) nucleon
+// cross-sections and the hadron-nucleous cross-section model in
+// the framework of Glauber-Gribov approach
+//
+//
+//
+//
+//
+// 25.04.12 V. Grichine - first implementation based on G4GlauberGribovCrossSection old interface
+//
+//
+
+#ifndef G4ComponentGGHadronNucleusXsc_h
+#define G4ComponentGGHadronNucleusXsc_h 1
+
+#include "globals.hh"
+#include "G4Proton.hh"
+#include "G4Nucleus.hh"
+
+#include "G4VComponentCrossSection.hh"
+
+class G4ParticleDefinition;
+class G4HadronNucleonXsc;
+
+class G4ComponentGGHadronNucleusXsc : public G4VComponentCrossSection
+{
+public:
+
+ G4ComponentGGHadronNucleusXsc ();
+ virtual ~G4ComponentGGHadronNucleusXsc ();
+
+
+ // virtual interface methods
+
+virtual
+ G4double GetTotalIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A);
+
+virtual
+ G4double GetTotalElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A);
+
+virtual
+ G4double GetInelasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A);
+
+virtual
+ G4double GetProductionIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A);
+
+virtual
+ G4double GetInelasticElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A);
+virtual
+ G4double GetProductionElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A);
+
+virtual
+ G4double GetElasticElementCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4double A);
+
+virtual
+ G4double GetElasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A);
+
+virtual
+ G4double ComputeQuasiElasticRatio(const G4ParticleDefinition* aParticle,
+ G4double kinEnergy,
+ G4int Z, G4int A);
+
+
+
+ // virtual
+ G4bool IsIsoApplicable(const G4DynamicParticle* aDP, G4int Z, G4int A,
+ const G4Element* elm = 0,
+ const G4Material* mat = 0);
+
+ // virtual
+ G4double GetIsoCrossSection(const G4DynamicParticle*, G4int Z, G4int A,
+ const G4Isotope* iso = 0,
+ const G4Element* elm = 0,
+ const G4Material* mat = 0);
+
+ G4double GetRatioSD(const G4DynamicParticle*, G4int At, G4int Zt);
+ G4double GetRatioQE(const G4DynamicParticle*, G4int At, G4int Zt);
+
+ G4double GetHadronNucleonXsc(const G4DynamicParticle*, const G4Element*);
+ G4double GetHadronNucleonXsc(const G4DynamicParticle*, G4int At, G4int Zt);
+
+ G4double GetHadronNucleonXscPDG(const G4DynamicParticle*, const G4Element*);
+ G4double GetHadronNucleonXscPDG(const G4DynamicParticle*, G4int At, G4int Zt);
+ G4double GetHadronNucleonXscNS(const G4DynamicParticle*, const G4Element*);
+ G4double GetHadronNucleonXscNS(const G4DynamicParticle*, G4int At, G4int Zt);
+ G4double GetKaonNucleonXscVector(const G4DynamicParticle*, G4int At, G4int Zt);
+
+ G4double GetHNinelasticXsc(const G4DynamicParticle*, const G4Element*);
+ G4double GetHNinelasticXsc(const G4DynamicParticle*, G4int At, G4int Zt);
+ G4double GetHNinelasticXscVU(const G4DynamicParticle*, G4int At, G4int Zt);
+
+ G4double CalculateEcmValue ( const G4double , const G4double , const G4double );
+
+ G4double CalcMandelstamS( const G4double , const G4double , const G4double );
+
+ G4double GetNucleusRadius(const G4DynamicParticle*, const G4Element*);
+ G4double GetNucleusRadius(G4int At);
+
+ virtual void CrossSectionDescription(std::ostream&) const;
+
+ inline G4double GetElasticGlauberGribov(const G4DynamicParticle*, G4int Z, G4int A);
+ inline G4double GetInelasticGlauberGribov(const G4DynamicParticle*, G4int Z, G4int A);
+
+ inline G4double GetTotalGlauberGribovXsc() { return fTotalXsc; };
+ inline G4double GetElasticGlauberGribovXsc() { return fElasticXsc; };
+ inline G4double GetInelasticGlauberGribovXsc(){ return fInelasticXsc; };
+ inline G4double GetProductionGlauberGribovXsc(){ return fProductionXsc; };
+ inline G4double GetDiffractionGlauberGribovXsc(){ return fDiffractionXsc; };
+ inline G4double GetRadiusConst() { return fRadiusConst; };
+
+ inline G4double GetParticleBarCorTot(const G4ParticleDefinition* theParticle, G4int Z);
+ inline G4double GetParticleBarCorIn(const G4ParticleDefinition* theParticle, G4int Z);
+
+ inline void SetEnergyLowerLimit(G4double E ){fLowerLimit=E;};
+
+private:
+
+// const G4double fUpperLimit;
+ G4double fLowerLimit;
+ const G4double fRadiusConst;
+
+ static const G4double fNeutronBarCorrectionTot[93];
+ static const G4double fNeutronBarCorrectionIn[93];
+
+ static const G4double fProtonBarCorrectionTot[93];
+ static const G4double fProtonBarCorrectionIn[93];
+
+ static const G4double fPionPlusBarCorrectionTot[93];
+ static const G4double fPionPlusBarCorrectionIn[93];
+
+ static const G4double fPionMinusBarCorrectionTot[93];
+ static const G4double fPionMinusBarCorrectionIn[93];
+
+ G4double fTotalXsc, fElasticXsc, fInelasticXsc, fProductionXsc, fDiffractionXsc;
+// G4double fHadronNucleonXsc;
+
+ G4ParticleDefinition* theGamma;
+ G4ParticleDefinition* theProton;
+ G4ParticleDefinition* theNeutron;
+ G4ParticleDefinition* theAProton;
+ G4ParticleDefinition* theANeutron;
+ G4ParticleDefinition* thePiPlus;
+ G4ParticleDefinition* thePiMinus;
+ G4ParticleDefinition* thePiZero;
+ G4ParticleDefinition* theKPlus;
+ G4ParticleDefinition* theKMinus;
+ G4ParticleDefinition* theK0S;
+ G4ParticleDefinition* theK0L;
+ G4ParticleDefinition* theL;
+ G4ParticleDefinition* theAntiL;
+ G4ParticleDefinition* theSPlus;
+ G4ParticleDefinition* theASPlus;
+ G4ParticleDefinition* theSMinus;
+ G4ParticleDefinition* theASMinus;
+ G4ParticleDefinition* theS0;
+ G4ParticleDefinition* theAS0;
+ G4ParticleDefinition* theXiMinus;
+ G4ParticleDefinition* theXi0;
+ G4ParticleDefinition* theAXiMinus;
+ G4ParticleDefinition* theAXi0;
+ G4ParticleDefinition* theOmega;
+ G4ParticleDefinition* theAOmega;
+ G4ParticleDefinition* theD;
+ G4ParticleDefinition* theT;
+ G4ParticleDefinition* theA;
+ G4ParticleDefinition* theHe3;
+
+ G4HadronNucleonXsc* hnXsc;
+
+};
+
+////////////////////////////////////////////////////////////////
+//
+// Inlines
+
+inline
+G4double
+G4ComponentGGHadronNucleusXsc::GetElasticGlauberGribov(const G4DynamicParticle* dp,
+ G4int Z, G4int A)
+{
+ GetIsoCrossSection(dp, Z, A);
+ return fElasticXsc;
+}
+
+/////////////////////////////////////////////////////////////////
+
+inline
+G4double
+G4ComponentGGHadronNucleusXsc::GetInelasticGlauberGribov(const G4DynamicParticle* dp,
+ G4int Z, G4int A)
+{
+ GetIsoCrossSection(dp, Z, A);
+ return fInelasticXsc;
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// return correction at Tkin = 90*GeV GG -> Barashenkov tot xsc, when it
+// is available, else return 1.0
+
+
+inline G4double G4ComponentGGHadronNucleusXsc::GetParticleBarCorTot(
+ const G4ParticleDefinition* theParticle, G4int Z)
+{
+ if(Z >= 2 && Z <= 92)
+ {
+ if( theParticle == theProton ) return fProtonBarCorrectionTot[Z];
+ else if( theParticle == theNeutron) return fNeutronBarCorrectionTot[Z];
+ else if( theParticle == thePiPlus ) return fPionPlusBarCorrectionTot[Z];
+ else if( theParticle == thePiMinus) return fPionMinusBarCorrectionTot[Z];
+ else return 1.0;
+ }
+ else return 1.0;
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// return correction at Tkin = 90*GeV GG -> Barashenkov in xsc, when it
+// is available, else return 1.0
+
+
+inline G4double G4ComponentGGHadronNucleusXsc::GetParticleBarCorIn(
+ const G4ParticleDefinition* theParticle, G4int Z)
+{
+ if(Z >= 2 && Z <= 92)
+ {
+ if( theParticle == theProton ) return fProtonBarCorrectionIn[Z];
+ else if( theParticle == theNeutron) return fNeutronBarCorrectionIn[Z];
+ else if( theParticle == thePiPlus ) return fPionPlusBarCorrectionIn[Z];
+ else if( theParticle == thePiMinus) return fPionMinusBarCorrectionIn[Z];
+ else return 1.0;
+ }
+ else return 1.0;
+}
+
+#endif
diff --git a/Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.cc b/Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.cc
new file mode 100644
index 000000000..03e14623c
--- /dev/null
+++ b/Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.cc
@@ -0,0 +1,1676 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+// 14.03.07 V. Grichine - first implementation
+//
+
+#include "G4HadronNucleonXsc.hh"
+
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4HadTmpUtil.hh"
+
+
+G4HadronNucleonXsc::G4HadronNucleonXsc()
+:
+// fUpperLimit( 10000 * GeV ),
+ fLowerLimit( 0.03 * MeV ),
+ fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0)
+// , fHadronNucleonXsc(0.0)
+{
+ theGamma = G4Gamma::Gamma();
+ theProton = G4Proton::Proton();
+ theNeutron = G4Neutron::Neutron();
+ theAProton = G4AntiProton::AntiProton();
+ theANeutron = G4AntiNeutron::AntiNeutron();
+ thePiPlus = G4PionPlus::PionPlus();
+ thePiMinus = G4PionMinus::PionMinus();
+ thePiZero = G4PionZero::PionZero();
+ theKPlus = G4KaonPlus::KaonPlus();
+ theKMinus = G4KaonMinus::KaonMinus();
+ theK0S = G4KaonZeroShort::KaonZeroShort();
+ theK0L = G4KaonZeroLong::KaonZeroLong();
+ theL = G4Lambda::Lambda();
+ theAntiL = G4AntiLambda::AntiLambda();
+ theSPlus = G4SigmaPlus::SigmaPlus();
+ theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
+ theSMinus = G4SigmaMinus::SigmaMinus();
+ theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
+ theS0 = G4SigmaZero::SigmaZero();
+ theAS0 = G4AntiSigmaZero::AntiSigmaZero();
+ theXiMinus = G4XiMinus::XiMinus();
+ theXi0 = G4XiZero::XiZero();
+ theAXiMinus = G4AntiXiMinus::AntiXiMinus();
+ theAXi0 = G4AntiXiZero::AntiXiZero();
+ theOmega = G4OmegaMinus::OmegaMinus();
+ theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
+ theD = G4Deuteron::Deuteron();
+ theT = G4Triton::Triton();
+ theA = G4Alpha::Alpha();
+ theHe3 = G4He3::He3();
+
+ InitialiseKaonNucleonTotXsc();
+}
+
+
+G4HadronNucleonXsc::~G4HadronNucleonXsc()
+{}
+
+void G4HadronNucleonXsc::CrossSectionDescription(std::ostream& outFile) const
+{
+ outFile << "G4HadronNucleonXsc calculates the total, inelastic and elastic\n"
+ << "hadron-nucleon cross sections using several different\n"
+ << "parameterizations within the Glauber-Gribov approach. It is\n"
+ << "valid for all incident gammas and long-lived hadrons at\n"
+ << "energies above 30 keV. This is a cross section component which\n"
+ << "is to be used to build a cross section data set.\n";
+}
+
+G4bool
+G4HadronNucleonXsc::IsApplicable(const G4DynamicParticle* aDP,
+ const G4Element* anElement)
+{
+ G4int Z = G4lrint(anElement->GetZ());
+ G4int A = G4lrint(anElement->GetN());
+ return IsIsoApplicable(aDP, Z, A);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////
+//
+
+G4bool
+G4HadronNucleonXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
+ G4int Z, G4int)
+{
+ G4bool applicable = false;
+ // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
+ G4double kineticEnergy = aDP->GetKineticEnergy();
+
+ const G4ParticleDefinition* theParticle = aDP->GetDefinition();
+
+ if ( ( kineticEnergy >= fLowerLimit &&
+ Z > 1 && // >= He
+ ( theParticle == theAProton ||
+ theParticle == theGamma ||
+ theParticle == theKPlus ||
+ theParticle == theKMinus ||
+ theParticle == theSMinus) ) ||
+
+ ( kineticEnergy >= 0.1*fLowerLimit &&
+ Z > 1 && // >= He
+ ( theParticle == theProton ||
+ theParticle == theNeutron ||
+ theParticle == thePiPlus ||
+ theParticle == thePiMinus ) ) ) applicable = true;
+
+ return applicable;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon Xsc according to differnt parametrisations:
+// [2] E. Levin, hep-ph/9710546
+// [3] U. Dersch, et al, hep-ex/9910052
+// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
+
+G4double
+G4HadronNucleonXsc::GetHadronNucleonXscEL(const G4DynamicParticle* aParticle,
+ const G4ParticleDefinition* nucleon )
+{
+ G4double xsection;
+
+
+ G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
+
+ G4double proj_mass = aParticle->GetMass();
+ G4double proj_momentum = aParticle->GetMomentum().mag();
+ G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
+
+ sMand /= GeV*GeV; // in GeV for parametrisation
+ proj_momentum /= GeV;
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
+
+
+ if(theParticle == theGamma && pORn )
+ {
+ xsection = (0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theNeutron && pORn ) // as proton ???
+ {
+ xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theProton && pORn )
+ {
+ xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
+
+ // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
+ // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
+ }
+ else if(theParticle == theAProton && pORn )
+ {
+ xsection = ( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == thePiPlus && pORn )
+ {
+ xsection = (13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == thePiMinus && pORn )
+ {
+ // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
+ xsection = (13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theKPlus && pORn )
+ {
+ xsection = (11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
+ }
+ else if(theParticle == theKMinus && pORn )
+ {
+ xsection = (11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
+ }
+ else // as proton ???
+ {
+ xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
+ }
+ xsection *= millibarn;
+
+ fTotalXsc = xsection;
+ fInelasticXsc = 0.83*xsection;
+ fElasticXsc = fTotalXsc - fInelasticXsc;
+ if (fElasticXsc < 0.)fElasticXsc = 0.;
+
+ return xsection;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
+// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
+// At = number of nucleons, Zt = number of protons
+
+G4double
+G4HadronNucleonXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
+ const G4ParticleDefinition* nucleon )
+{
+ G4double xsection(0);
+ G4int Zt=1, Nt=1, At=1;
+
+ G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
+
+ G4double proj_mass = aParticle->GetMass();
+ G4double proj_momentum = aParticle->GetMomentum().mag();
+
+ G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
+
+ sMand /= GeV*GeV; // in GeV for parametrisation
+
+ // General PDG fit constants
+
+ G4double s0 = 5.38*5.38; // in Gev^2
+ G4double eta1 = 0.458;
+ G4double eta2 = 0.458;
+ G4double B = 0.308;
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
+ G4bool proton = (nucleon == theProton);
+ G4bool neutron = (nucleon == theNeutron);
+
+ if(theParticle == theNeutron) // proton-neutron fit
+ {
+ if ( proton )
+ {
+ xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));// on p
+ }
+ if ( neutron )
+ {
+ xsection = Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // on n pp for nn
+ }
+ }
+ else if(theParticle == theProton)
+ {
+ if ( proton )
+ {
+ xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
+ }
+ if ( neutron )
+ {
+ xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
+ }
+ }
+ else if(theParticle == theAProton)
+ {
+ if ( proton )
+ {
+ xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
+ }
+ if ( neutron )
+ {
+ xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
+ }
+ }
+ else if(theParticle == thePiPlus && pORn )
+ {
+ xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
+ + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == thePiMinus && pORn )
+ {
+ xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
+ + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
+ }
+ else if(theParticle == theKPlus)
+ {
+ if ( proton )
+ {
+ xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
+ + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
+ }
+ if ( neutron )
+ {
+ xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
+ + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
+ }
+ }
+ else if(theParticle == theKMinus)
+ {
+ if ( proton )
+ {
+ xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
+ + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
+ }
+ if ( neutron )
+ {
+ xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
+ + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2) );
+ }
+ }
+ else if(theParticle == theSMinus && pORn )
+ {
+ xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
+ - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2) );
+ }
+ else if(theParticle == theGamma && pORn ) // modify later on
+ {
+ xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
+ + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2) );
+
+ }
+ else // as proton ???
+ {
+ if ( proton )
+ {
+ xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2) );
+ }
+ if ( neutron )
+ {
+ xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
+ }
+ }
+ xsection *= millibarn; // parametrised in mb
+
+ fTotalXsc = xsection;
+ fInelasticXsc = 0.75*xsection;
+ fElasticXsc = fTotalXsc - fInelasticXsc;
+ if (fElasticXsc < 0.) fElasticXsc = 0.;
+
+ return xsection;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
+// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
+
+G4double
+G4HadronNucleonXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
+ const G4ParticleDefinition* nucleon )
+{
+ G4double xsection(0);
+
+ G4double A0, B0;
+ G4double hpXsc(0);
+ G4double hnXsc(0);
+
+
+ G4double tM = 0.939*GeV; // ~mean neutron and proton ???
+
+ G4double pM = aParticle->GetMass();
+ G4double pE = aParticle->GetTotalEnergy();
+ G4double pLab = aParticle->GetMomentum().mag();
+
+ G4double sMand = CalcMandelstamS ( pM , tM , pLab );
+
+ sMand /= GeV*GeV; // in GeV for parametrisation
+ pLab /= GeV;
+ pE /= GeV;
+ pM /= GeV;
+
+ G4double logP = std::log(pLab);
+
+
+ // General PDG fit constants
+
+ G4double s0 = 5.38*5.38; // in Gev^2
+ G4double eta1 = 0.458;
+ G4double eta2 = 0.458;
+ G4double B = 0.308;
+ G4double minLogP = 3.5; // min of (lnP-minLogP)^2
+ G4double cofLogE = .0557; // elastic (lnP-minLogP)^2
+ G4double cofLogT = .3; // total (lnP-minLogP)^2
+ G4double pMin = .1; // fast LE calculation
+ G4double pMax = 1000.; // fast HE calculation
+
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
+ G4bool proton = (nucleon == theProton);
+ G4bool neutron = (nucleon == theNeutron);
+
+ if( theParticle == theNeutron && pORn )
+ {
+ if( pLab >= 373.)
+ {
+ xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
+
+ fElasticXsc = 6.5 + 0.308*std::pow(std::log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
+
+ fTotalXsc = xsection;
+
+ }
+ else if( pLab >= 100.)
+ {
+ B0 = 7.5;
+ A0 = 100. - B0*std::log(3.0e7);
+
+ xsection = A0 + B0*std::log(pE) - 11
+ // + 103*std::pow(2*0.93827*pE + pM*pM+0.93827*0.93827,-0.165); // mb
+ + 103*std::pow(sMand,-0.165); // mb
+
+ fElasticXsc = 5.53 + 0.308*std::pow(std::log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
+
+ fTotalXsc = xsection;
+ }
+ else if( pLab >= 10.)
+ {
+ B0 = 7.5;
+ A0 = 100. - B0*std::log(3.0e7);
+
+ xsection = A0 + B0*std::log(pE) - 11
+ + 103*std::pow(2*0.93827*pE + pM*pM+
+ 0.93827*0.93827,-0.165); // mb
+ fTotalXsc = xsection;
+ fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
+ }
+ else // pLab < 10 GeV/c
+ {
+ if( neutron ) // nn to be pp
+ {
+ if( pLab < 0.4 )
+ {
+ hnXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
+ fElasticXsc = hnXsc;
+ }
+ else if( pLab < 0.73 )
+ {
+ hnXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
+ fElasticXsc = hnXsc;
+ }
+ else if( pLab < 1.05 )
+ {
+ hnXsc = 23 + 40*(std::log(pLab/0.73))*
+ (std::log(pLab/0.73));
+ fElasticXsc = 23 + 20*(std::log(pLab/0.73))*
+ (std::log(pLab/0.73));
+ }
+ else // 1.05 - 10 GeV/c
+ {
+ hnXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
+
+ fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
+ }
+ fTotalXsc = hnXsc;
+ }
+ if( proton ) // pn to be np
+ {
+ if( pLab < 0.02 )
+ {
+ hpXsc = 4100+30*std::pow(std::log(1.3/pLab),3.6); // was as pLab < 0.8
+ fElasticXsc = hpXsc;
+ }
+ else if( pLab < 0.8 )
+ {
+ hpXsc = 33+30*std::pow(std::log(pLab/1.3),4.0);
+ fElasticXsc = hpXsc;
+ }
+ else if( pLab < 1.05 )
+ {
+ hpXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
+ fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
+ }
+ else if( pLab < 1.4 )
+ {
+ hpXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
+ fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
+ }
+ else // 1.4 < pLab < 10. )
+ {
+ hpXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
+
+ fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
+ }
+ fTotalXsc = hpXsc;
+ }
+ }
+ }
+ else if( theParticle == theProton && pORn ) ////// proton //////////////////////////////////////////////
+ {
+ if( pLab >= 373.) // pdg due to TOTEM data
+ {
+ xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
+
+ fElasticXsc = 6.5 + 0.308*std::pow(std::log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
+
+ fTotalXsc = xsection;
+ }
+ else if( pLab >= 100.)
+ {
+ B0 = 7.5;
+ A0 = 100. - B0*std::log(3.0e7);
+
+ xsection = A0 + B0*std::log(pE) - 11 + 103*std::pow(sMand,-0.165); // mb
+
+ fElasticXsc = 5.53 + 0.308*std::pow(std::log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
+
+ fTotalXsc = xsection;
+ }
+ else if( pLab >= 10.)
+ {
+ B0 = 7.5;
+ A0 = 100. - B0*std::log(3.0e7);
+
+ xsection = A0 + B0*std::log(pE) - 11 + 103*std::pow(sMand,-0.165); // mb
+
+ fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
+
+ fTotalXsc = xsection;
+ }
+ else
+ {
+ // pp
+
+ if( proton )
+ {
+ if( pLab < 0.4 )
+ {
+ hpXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
+ fElasticXsc = hpXsc;
+ }
+ else if( pLab < 0.73 )
+ {
+ hpXsc = 23 + 50*( std::pow( std::log(0.73/pLab), 3.5 ) );
+ fElasticXsc = hpXsc;
+ }
+ else if( pLab < 1.05 )
+ {
+ hpXsc = 23 + 40*(std::log(pLab/0.73))*
+ (std::log(pLab/0.73));
+ fElasticXsc = 23 + 20*(std::log(pLab/0.73))*
+ (std::log(pLab/0.73));
+ }
+ else // 1.05 - 10 GeV/c
+ {
+ hpXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
+
+ fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
+ }
+ fTotalXsc = hpXsc;
+ }
+ if( neutron ) // pn to be np
+ {
+ if( pLab < 0.02 )
+ {
+ hnXsc = 4100+30*std::pow(std::log(1.3/pLab),3.6); // was as pLab < 0.8
+ fElasticXsc = hnXsc;
+ }
+ else if( pLab < 0.8 )
+ {
+ hnXsc = 33+30*std::pow(std::log(pLab/1.3),4.0);
+ fElasticXsc = hnXsc;
+ }
+ else if( pLab < 1.05 )
+ {
+ hnXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
+ fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
+ }
+ else if( pLab < 1.4 )
+ {
+ hnXsc = 33+30*std::pow(std::log(pLab/0.95),2.0);
+ fElasticXsc = 6 + 52/( std::log(0.511/pLab)*std::log(0.511/pLab) + 1.6 );
+ }
+ else // 1.4 < pLab < 10. )
+ {
+ hnXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
+
+ fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
+ }
+ fTotalXsc = hnXsc;
+ }
+ }
+ }
+ else if( theParticle == theAProton && pORn ) /////////////////// p_bar ///////////////////////////
+ {
+ if( proton )
+ {
+ xsection = 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2);
+ }
+ if( neutron ) // ???
+ {
+ xsection = 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2);
+ }
+ fTotalXsc = xsection;
+ }
+ else if( theParticle == thePiPlus && pORn ) // pi+ /////////////////////////////////////////////
+ {
+ if( proton ) // pi+ p
+ {
+ if( pLab < 0.28 )
+ {
+ hpXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
+ fElasticXsc = hpXsc;
+ }
+ else if( pLab < 0.4 )
+ {
+ hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
+ fElasticXsc = hpXsc;
+ }
+ else if( pLab < 0.68 )
+ {
+ hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
+ fElasticXsc = hpXsc;
+ }
+ else if( pLab < 0.85 )
+ {
+ G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
+ hpXsc = Ex4 + 14.9;
+ fElasticXsc = hpXsc*std::exp(-3.*(pLab - 0.68));
+ }
+ else if( pLab < 1.15 )
+ {
+ G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
+ hpXsc = Ex4 + 14.9;
+
+ fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
+ }
+ else if( pLab < 1.4) // ns original
+ {
+ G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
+ hpXsc = Ex1 + Ex2 + 27.5;
+ fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
+ }
+ else if( pLab < 2.0 ) // ns original
+ {
+ G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
+ hpXsc = Ex1 + Ex2 + 27.5;
+ fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
+ }
+ else if( pLab < 3.5 ) // ns original
+ {
+ G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
+ hpXsc = Ex1 + Ex2 + 27.5;
+ fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
+ }
+ else if( pLab < 200. ) // my
+ {
+ hpXsc = 10.6 + 2.*std::log(pE) + 25*std::pow(pE, -0.43 ); // ns original
+ // hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
+ fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
+ }
+ else // pLab > 100 // my
+ {
+ hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
+ fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
+ }
+ fTotalXsc = hpXsc;
+ }
+ if( neutron ) // pi+ n = pi- p??
+ {
+ if( pLab < 0.28 )
+ {
+ hnXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
+ fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
+ }
+ else if( pLab < 0.395676 ) // first peak
+ {
+ hnXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
+ fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
+ }
+ else if( pLab < 0.5 )
+ {
+ hnXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
+ fElasticXsc = 0.37*hnXsc;
+ }
+ else if( pLab < 0.65 )
+ {
+ hnXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
+ fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
+ }
+ else if( pLab < 0.72 )
+ {
+ hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
+ }
+ else if( pLab < 0.88 )
+ {
+ hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
+ }
+ else if( pLab < 1.03 )
+ {
+ hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
+ }
+ else if( pLab < 1.15 )
+ {
+ hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
+ }
+ else if( pLab < 1.3 )
+ {
+ hnXsc = 36.1 + 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 3. + 13./pLab;
+ }
+ else if( pLab < 2.6 ) // < 3.0) // ns original
+ {
+ hnXsc = 36.1 + 0.079-4.313*std::log(pLab)+
+ 3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
+ 1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
+ fElasticXsc = 3. + 13./pLab;
+ }
+ else if( pLab < 20. ) // < 3.0) // ns original
+ {
+ hnXsc = 36.1 + 0.079 - 4.313*std::log(pLab)+
+ 3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
+ 1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
+ fElasticXsc = 3. + 13./pLab;
+ }
+ else // mb
+ {
+ hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
+ fElasticXsc = 3. + 13./pLab;
+ }
+ fTotalXsc = hnXsc;
+ }
+ }
+ else if( theParticle == thePiMinus && pORn ) /// pi- ////////////////////////////////////////////
+ {
+ if( neutron ) // pi- n = pi+ p??
+ {
+ if( pLab < 0.28 )
+ {
+ hnXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
+ fElasticXsc = hnXsc;
+ }
+ else if( pLab < 0.4 )
+ {
+ hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
+ fElasticXsc = hnXsc;
+ }
+ else if( pLab < 0.68 )
+ {
+ hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
+ fElasticXsc = hnXsc;
+ }
+ else if( pLab < 0.85 )
+ {
+ G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
+ hnXsc = Ex4 + 14.9;
+ fElasticXsc = hnXsc*std::exp(-3.*(pLab - 0.68));
+ }
+ else if( pLab < 1.15 )
+ {
+ G4double Ex4 = 88*(std::log(pLab/0.77))*(std::log(pLab/0.77));
+ hnXsc = Ex4 + 14.9;
+
+ fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
+ }
+ else if( pLab < 1.4) // ns original
+ {
+ G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
+ hnXsc = Ex1 + Ex2 + 27.5;
+ fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
+ }
+ else if( pLab < 2.0 ) // ns original
+ {
+ G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
+ hnXsc = Ex1 + Ex2 + 27.5;
+ fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
+ }
+ else if( pLab < 3.5 ) // ns original
+ {
+ G4double Ex1 = 3.2*std::exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
+ G4double Ex2 = 12*std::exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
+ hnXsc = Ex1 + Ex2 + 27.5;
+ fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
+ }
+ else if( pLab < 200. ) // my
+ {
+ hnXsc = 10.6 + 2.*std::log(pE) + 25*std::pow(pE, -0.43 ); // ns original
+ fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
+ }
+ else // pLab > 100 // my
+ {
+ hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
+ fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
+ }
+ fTotalXsc = hnXsc;
+ }
+ if( proton ) // pi- p
+ {
+ if( pLab < 0.28 )
+ {
+ hpXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
+ fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
+ }
+ else if( pLab < 0.395676 ) // first peak
+ {
+ hpXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
+ fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
+ }
+ else if( pLab < 0.5 )
+ {
+ hpXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
+ fElasticXsc = 0.37*hpXsc;
+ }
+ else if( pLab < 0.65 )
+ {
+ hpXsc = 26 + 110*(std::log(pLab/0.48))*(std::log(pLab/0.48));
+ fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
+ }
+ else if( pLab < 0.72 )
+ {
+ hpXsc = 36.1+
+ 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
+ }
+ else if( pLab < 0.88 )
+ {
+ hpXsc = 36.1+
+ 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
+ }
+ else if( pLab < 1.03 )
+ {
+ hpXsc = 36.1+
+ 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
+ }
+ else if( pLab < 1.15 )
+ {
+ hpXsc = 36.1+
+ 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
+ }
+ else if( pLab < 1.3 )
+ {
+ hpXsc = 36.1+
+ 10*std::exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
+ 24*std::exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
+ fElasticXsc = 3. + 13./pLab;
+ }
+ else if( pLab < 2.6 ) // < 3.0) // ns original
+ {
+ hpXsc = 36.1+0.079-4.313*std::log(pLab)+
+ 3*std::exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
+ 1.5*std::exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
+ fElasticXsc = 3. +13./pLab; // *std::log(pLab*6.79);
+ }
+ else // mb
+ {
+ hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
+ fElasticXsc = 3. + 13./pLab;
+ }
+ fTotalXsc = hpXsc;
+ }
+ }
+ else if( (theParticle == theKMinus || theParticle == theK0S) && proton ) // Kmp/K0p /////////////////////////////////
+ {
+ if( pLab < pMin)
+ {
+ G4double psp = pLab*std::sqrt(pLab);
+ fElasticXsc = 5.2/psp;
+ fTotalXsc = 14./psp;
+ }
+ else if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
+ }
+ else
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ G4double sp = std::sqrt(pLab);
+ G4double psp = pLab*sp;
+ G4double p2 = pLab*pLab;
+ G4double p4 = p2*p2;
+ G4double lm = pLab - .39;
+ G4double md = lm*lm + .000356;
+
+ G4double lh1 = pLab - 0.78;
+ G4double hd1 = lh1*lh1 + .00166;
+
+ G4double lh = pLab - 1.01;
+ G4double hd = lh*lh + .011;
+
+ G4double lh2 = pLab - 1.63;
+ G4double hd2 = lh2*lh2 + .007;
+
+ fElasticXsc = 5.2/psp + (1.1*cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4)
+ + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd; // small peaks were added
+
+ fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
+ + .006/md + 0.01/hd1+ 0.02/hd2 + .20/hd ;
+ }
+ }
+ else if( (theParticle == theKMinus || theParticle == theK0S) && neutron ) // Kmn/K0n //////////////////////////////
+ {
+ if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
+ }
+ else
+ {
+
+ G4double lh = pLab - 0.98;
+ G4double hd = lh*lh + .021;
+
+ G4double LogPlab = std::log( pLab );
+ G4double sqrLogPlab = LogPlab * LogPlab;
+
+ fElasticXsc = // 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md
+ 5.0 + 8.1*std::pow(pLab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab + .15/hd;
+ fTotalXsc = // 14./psp +
+ // (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
+ 25.2 + 0. *std::pow(pLab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab
+ // + .006/md + 0.01/hd1+ 0.02/hd2
+ + 0.30/hd ;
+ }
+ }
+ else if( (theParticle == theKPlus || theParticle == theK0L) && proton ) // Kpp/aKp ////////////////////////
+ {
+ if( pLab < pMin )
+ {
+ G4double lr = pLab - .38;
+ G4double lm = pLab - 1.;
+ G4double md = lm*lm + .392;
+ fElasticXsc = .7/(lr*lr + .076) + 2./md;
+ fTotalXsc = .7/(lr*lr + .076) + 2.6/md;
+ }
+ else if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = cofLogT*ld2 + 19.2;
+ }
+ else
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ G4double lr = pLab - .38;
+ G4double LE = .7/(lr*lr + .076);
+ G4double sp = std::sqrt(pLab);
+ G4double p2 = pLab*pLab;
+ G4double p4 = p2*p2;
+ G4double lm = pLab - 1.;
+ G4double md = lm*lm + .392;
+ fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
+ fTotalXsc = LE + (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 2.6/md;
+ }
+ }
+ else if( (theParticle == theKPlus || theParticle == theK0L) && neutron ) // Kpn/aKn ///////////////////////
+ {
+ if( pLab < pMin )
+ {
+ G4double lm = pLab - 0.94;
+ G4double md = lm*lm + .392;
+ fElasticXsc = 2./md;
+ fTotalXsc = 4.6/md;
+ }
+ else if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = cofLogT*ld2 + 19.2;
+ }
+ else
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ G4double sp = std::sqrt(pLab);
+ G4double p2 = pLab*pLab;
+ G4double p4 = p2*p2;
+ G4double lm = pLab - 0.94;
+ G4double md = lm*lm + .392;
+ fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
+ fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 4.6/md;
+ }
+ }
+ else if( theParticle == theSMinus && pORn )
+ {
+ xsection = 35.20 + B*std::pow(std::log(sMand/s0),2.)
+ - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2);
+ }
+ else if( theParticle == theGamma && pORn ) // modify later on
+ {
+ xsection = 0.0 + B*std::pow(std::log(sMand/s0),2.)
+ + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2);
+ fTotalXsc = xsection;
+ }
+ else // other then p,n,pi+,pi-,K+,K- as proton ???
+ {
+ if( proton )
+ {
+ xsection = 35.45 + B*std::pow(std::log(sMand/s0),2.)
+ + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2);
+ }
+ if( neutron )
+ {
+ xsection += 35.80 + B*std::pow(std::log(sMand/s0),2.)
+ + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2);
+ }
+ fTotalXsc = xsection;
+ }
+ fTotalXsc *= millibarn; // parametrised in mb
+ fElasticXsc *= millibarn; // parametrised in mb
+
+ if( proton && aParticle->GetDefinition()->GetPDGCharge() > 0. )
+ {
+ G4double cB = GetCoulombBarrier(aParticle, nucleon);
+ fTotalXsc *= cB;
+ fElasticXsc *= cB;
+ }
+ fInelasticXsc = fTotalXsc - fElasticXsc;
+ if( fInelasticXsc < 0. ) fInelasticXsc = 0.;
+
+ // G4cout<<fTotalXsc/millibarn<<"; "<<fElasticXsc/millibarn<<"; "<<fInelasticXsc/millibarn<<G4endl;
+
+ return fTotalXsc;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns kaon-nucleon cross-section based on smoothed NS for GG model
+
+G4double
+G4HadronNucleonXsc::GetKaonNucleonXscGG(const G4DynamicParticle* aParticle,
+ const G4ParticleDefinition* nucleon )
+{
+ G4double pLab = aParticle->GetMomentum().mag();
+
+ pLab /= GeV;
+ G4double LogPlab = std::log( pLab );
+ G4double sqrLogPlab = LogPlab * LogPlab;
+
+ G4double minLogP = 3.5; // min of (lnP-minLogP)^2
+ G4double cofLogE = .0557; // elastic (lnP-minLogP)^2
+ G4double cofLogT = .3; // total (lnP-minLogP)^2
+ G4double pMin = .1; // fast LE calculation
+ G4double pMax = 1000.; // fast HE calculation
+
+ const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
+
+ G4bool proton = (nucleon == theProton);
+ G4bool neutron = (nucleon == theNeutron);
+
+ if( (theParticle == theKMinus || theParticle == theK0S) && proton ) // (K-,K0)on p ////////////////////////////
+ {
+
+ if( pLab < pMin)
+ {
+ G4double psp = pLab*std::sqrt(pLab);
+ fElasticXsc = 5.2/psp;
+ fTotalXsc = 14./psp;
+ }
+ else if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
+ }
+ else
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ G4double sp = std::sqrt(pLab);
+ G4double psp = pLab*sp;
+ G4double p2 = pLab*pLab;
+ G4double p4 = p2*p2;
+
+ G4double lh = pLab - 0.98;
+ G4double hd = lh*lh + .045;
+
+
+ fElasticXsc = 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) // + .004/md
+ + .15/hd;
+ fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
+ // + .006/md + 0.01/hd1 + 0.02/hd2
+ + .60/hd;
+ }
+ }
+ else if( (theParticle == theKMinus || theParticle == theK0S) && neutron ) // Kmn/K0n /////////////////////////////
+ {
+ if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
+ }
+ else
+ {
+
+ G4double lh = pLab - 0.98;
+ G4double hd = lh*lh + .045;
+
+ fElasticXsc = // 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md
+ 5.0 + 8.1*std::pow(pLab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab + .15/hd;
+ fTotalXsc = // 14./psp +
+ // (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
+ 25.2 + 0. *std::pow(pLab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab
+ // + .006/md + 0.01/hd1+ 0.02/hd2
+ + 0.60/hd ;
+ }
+ }
+ else if( (theParticle == theKPlus || theParticle == theK0L) && proton ) // Kpp/aKp //////////////////////
+ {
+ if( pLab < pMin )
+ {
+ G4double lr = pLab - .38;
+ G4double lm = pLab - 1.;
+ G4double md = lm*lm + .392;
+ fElasticXsc = .7/(lr*lr + .076) + 2./md;
+ fTotalXsc = // .7/(lr*lr + .076) +
+ 2.6/md;
+ }
+ else if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = cofLogT*ld2 + 19.2;
+ }
+ else
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ G4double lr = pLab - .38;
+ G4double LE = .7/(lr*lr + .076);
+ G4double sp = std::sqrt(pLab);
+ G4double p2 = pLab*pLab;
+ G4double p4 = p2*p2;
+ G4double lm = pLab - 0.8;
+ G4double md = lm*lm + .652;
+ fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
+ fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 7.6/md; // + LE;
+ }
+ }
+ else if( (theParticle == theKPlus || theParticle == theK0L) && neutron ) // Kpn/aKn //////////////////////////////////
+ {
+ if( pLab < pMin )
+ {
+ G4double lm = pLab - 0.94;
+ G4double md = lm*lm + .392;
+ fElasticXsc = 2./md;
+ fTotalXsc = 4.6/md;
+ }
+ else if( pLab > pMax )
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ fElasticXsc = cofLogE*ld2 + 2.23;
+ fTotalXsc = cofLogT*ld2 + 19.2;
+ }
+ else
+ {
+ G4double ld = std::log(pLab) - minLogP;
+ G4double ld2 = ld*ld;
+ G4double sp = std::sqrt(pLab);
+ G4double p2 = pLab*pLab;
+ G4double p4 = p2*p2;
+ G4double lm = pLab - 0.8;
+ G4double md = lm*lm + .652;
+ fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
+ fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 7.6/md;
+ }
+ }
+ fTotalXsc *= millibarn; // parametrised in mb
+ fElasticXsc *= millibarn; // parametrised in mb
+
+ if( proton && aParticle->GetDefinition()->GetPDGCharge() > 0. )
+ {
+ G4double cB = GetCoulombBarrier(aParticle, nucleon);
+ fTotalXsc *= cB;
+ fElasticXsc *= cB;
+ }
+ fInelasticXsc = fTotalXsc - fElasticXsc;
+ if( fInelasticXsc < 0. ) fInelasticXsc = 0.;
+
+ // G4cout<<fTotalXsc/millibarn<<"; "<<fElasticXsc/millibarn<<"; "<<fInelasticXsc/millibarn<<G4endl;
+
+ return fTotalXsc;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+//
+// Returns hadron-nucleon cross-section based on V. Uzjinsky parametrisation of
+// data from G4FTFCrossSection class
+
+G4double
+G4HadronNucleonXsc::GetHadronNucleonXscVU(const G4DynamicParticle* aParticle,
+ const G4ParticleDefinition* nucleon )
+{
+ G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
+ G4int absPDGcode = std::abs(PDGcode);
+ G4double Elab = aParticle->GetTotalEnergy();
+ // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
+ G4double Plab = aParticle->GetMomentum().mag();
+ // std::sqrt(Elab * Elab - 0.88);
+
+ Elab /= GeV;
+ Plab /= GeV;
+
+ G4double LogPlab = std::log( Plab );
+ G4double sqrLogPlab = LogPlab * LogPlab;
+
+ G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
+ G4bool proton = (nucleon == theProton);
+ G4bool neutron = (nucleon == theNeutron);
+
+
+ if( absPDGcode > 1000 && pORn ) //------Projectile is baryon -
+ {
+ if(proton)
+ {
+ fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
+ fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
+ }
+ if(neutron)
+ {
+ fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
+ fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
+ }
+ }
+ else if( PDGcode == 211 && pORn ) //------Projectile is PionPlus ----
+ {
+ if(proton)
+ {
+ fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
+ fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
+ }
+ if(neutron)
+ {
+ fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
+ fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
+ }
+ }
+ else if( PDGcode == -211 && pORn ) //------Projectile is PionMinus ----
+ {
+ if(proton)
+ {
+ fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
+ fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
+ }
+ if(neutron)
+ {
+ fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
+ fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
+ }
+ }
+ else if( PDGcode == 111 && pORn ) //------Projectile is PionZero --
+ {
+ if(proton)
+ {
+ fTotalXsc = (16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
+ 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
+
+ fElasticXsc = ( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
+ 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
+
+ }
+ if(neutron)
+ {
+ fTotalXsc = (33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
+ 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
+ fElasticXsc = ( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
+ 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
+ }
+ }
+ else if( PDGcode == 321 && pORn ) //------Projectile is KaonPlus --
+ {
+ if(proton)
+ {
+ fTotalXsc = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
+ fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
+ }
+ if(neutron)
+ {
+ fTotalXsc = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
+ fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
+ }
+ }
+ else if( PDGcode ==-321 && pORn ) //------Projectile is KaonMinus ----
+ {
+ if(proton)
+ {
+ fTotalXsc = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66*sqrLogPlab - 5.6*LogPlab;
+ fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29*sqrLogPlab - 2.4*LogPlab;
+ }
+ if(neutron)
+ {
+ fTotalXsc = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab;
+ fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab;
+ }
+ }
+ else if( PDGcode == 311 && pORn ) //------Projectile is KaonZero -----
+ {
+ if(proton)
+ {
+ fTotalXsc = ( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
+ 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
+ fElasticXsc = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
+ 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
+ }
+ if(neutron)
+ {
+ fTotalXsc = ( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
+ 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
+ fElasticXsc = ( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
+ 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
+ }
+ }
+ else //------Projectile is undefined, Nucleon assumed
+ {
+ if(proton)
+ {
+ fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
+ fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
+ }
+ if(neutron)
+ {
+ fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
+ fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
+ }
+ }
+ fTotalXsc *= millibarn;
+ fElasticXsc *= millibarn;
+ fInelasticXsc = fTotalXsc - fElasticXsc;
+ if (fInelasticXsc < 0.) fInelasticXsc = 0.;
+
+ return fTotalXsc;
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double G4HadronNucleonXsc::CalculateEcmValue( const G4double mp ,
+ const G4double mt ,
+ const G4double Plab )
+{
+ G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
+ G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
+ // G4double Pcm = Plab * mt / Ecm;
+ // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
+
+ return Ecm ; // KEcm;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double G4HadronNucleonXsc::CalcMandelstamS( const G4double mp ,
+ const G4double mt ,
+ const G4double Plab )
+{
+ G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
+ G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
+
+ return sMand;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double G4HadronNucleonXsc::GetCoulombBarrier(const G4DynamicParticle* aParticle,
+ const G4ParticleDefinition* nucleon )
+{
+ G4double ratio;
+
+ G4double tR = 0.895*fermi, pR;
+
+ if ( aParticle->GetDefinition() == theProton ) pR = 0.895*fermi;
+ else if( aParticle->GetDefinition() == thePiPlus ) pR = 0.663*fermi;
+ else if( aParticle->GetDefinition() == theKPlus ) pR = 0.340*fermi;
+ else pR = 0.500*fermi;
+
+ G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
+ G4double tZ = nucleon->GetPDGCharge();
+
+ G4double pTkin = aParticle->GetKineticEnergy();
+
+ G4double pM = aParticle->GetDefinition()->GetPDGMass();
+ G4double tM = nucleon->GetPDGMass();
+
+ G4double pElab = pTkin + pM;
+
+ G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
+
+ G4double totTcm = totEcm - pM -tM;
+
+ G4double bC = fine_structure_const*hbarc*pZ*tZ;
+ bC /= pR + tR;
+ bC /= 2.; // 4., 2. parametrisation cof ??? vmg
+
+ // G4cout<<"pTkin = "<<pTkin/GeV<<"; pPlab = "
+ // <<pPlab/GeV<<"; bC = "<<bC/GeV<<"; pTcm = "<<pTcm/GeV<<G4endl;
+
+ if( totTcm <= bC ) ratio = 0.;
+ else ratio = 1. - bC/totTcm;
+
+ // if(ratio < DBL_MIN) ratio = DBL_MIN;
+ if( ratio < 0.) ratio = 0.;
+
+ // G4cout <<"ratio = "<<ratio<<G4endl;
+ return ratio;
+}
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////////
+//
+// Initialaise K(p,m)-(p,n) total cross section vectors
+
+
+void G4HadronNucleonXsc::InitialiseKaonNucleonTotXsc()
+{
+ G4int i = 0, iMax;
+ G4double tmpxsc[106];
+
+ // Kp-proton tot xsc
+
+ iMax = 66;
+ fKpProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKpProtonTotTkin[0], fKpProtonTotTkin[iMax-1]);
+ fKpProtonTotXscVector.SetSpline(true);
+
+ for( i = 0; i < iMax; i++)
+ {
+ tmpxsc[i] = 0.;
+
+ if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpProtonTotXsc[i];
+ else tmpxsc[i] = 0.5*(fKpProtonTotXsc[i-1]+fKpProtonTotXsc[i+1]);
+
+ fKpProtonTotXscVector.PutValues(size_t(i), fKpProtonTotTkin[i], tmpxsc[i]*millibarn);
+ }
+
+ // Kp-neutron tot xsc
+
+ iMax = 75;
+ fKpNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKpNeutronTotTkin[0], fKpNeutronTotTkin[iMax-1]);
+ fKpNeutronTotXscVector.SetSpline(true);
+
+ for( i = 0; i < iMax; i++)
+ {
+ tmpxsc[i] = 0.;
+ if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpNeutronTotXsc[i];
+ else tmpxsc[i] = 0.5*(fKpNeutronTotXsc[i-1]+fKpNeutronTotXsc[i+1]);
+
+ fKpNeutronTotXscVector.PutValues(size_t(i), fKpNeutronTotTkin[i], tmpxsc[i]*millibarn);
+ }
+
+ // Km-proton tot xsc
+
+ iMax = 106;
+ fKmProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKmProtonTotTkin[0], fKmProtonTotTkin[iMax-1]);
+ fKmProtonTotXscVector.SetSpline(true);
+
+ for( i = 0; i < iMax; i++)
+ {
+ tmpxsc[i] = 0.;
+
+ if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmProtonTotXsc[i];
+ else tmpxsc[i] = 0.5*(fKmProtonTotXsc[i-1]+fKmProtonTotXsc[i+1]);
+
+ fKmProtonTotXscVector.PutValues(size_t(i), fKmProtonTotTkin[i], tmpxsc[i]*millibarn);
+ }
+
+ // Km-neutron tot xsc
+
+ iMax = 68;
+ fKmNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKmNeutronTotTkin[0], fKmNeutronTotTkin[iMax-1]);
+ fKmNeutronTotXscVector.SetSpline(true);
+
+ for( i = 0; i < iMax; i++)
+ {
+ tmpxsc[i] = 0.;
+ if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmNeutronTotXsc[i];
+ else tmpxsc[i] = 0.5*(fKmNeutronTotXsc[i-1]+fKmNeutronTotXsc[i+1]);
+
+ fKmNeutronTotXscVector.PutValues(size_t(i), fKmNeutronTotTkin[i], tmpxsc[i]*millibarn);
+ }
+}
+
+///////////////////////////////////////////////////////
+//
+// K-nucleon tot xsc (mb) fit data, std::log(Tkin(MeV))
+
+const G4double G4HadronNucleonXsc::fKpProtonTotXsc[66] = {
+0.000000e+00, 1.592400e-01, 3.184700e-01, 7.961800e-01, 1.433120e+00, 2.070060e+00,
+2.866240e+00, 3.582800e+00, 4.378980e+00, 5.015920e+00, 5.573250e+00, 6.449040e+00,
+7.404460e+00, 8.200640e+00, 8.837580e+00, 9.713380e+00, 1.027070e+01, 1.090764e+01,
+1.130573e+01, 1.170382e+01, 1.242038e+01, 1.281847e+01, 1.321656e+01, 1.337580e+01,
+1.345541e+01, 1.329618e+01, 1.265924e+01, 1.242038e+01, 1.250000e+01, 1.305732e+01,
+1.369427e+01, 1.425159e+01, 1.544586e+01, 1.648089e+01, 1.751592e+01, 1.791401e+01,
+1.791401e+01, 1.775478e+01, 1.751592e+01, 1.735669e+01, 1.719745e+01, 1.711783e+01,
+1.703822e+01, 1.695860e+01, 1.695860e+01, 1.695860e+01, 1.695860e+01, 1.687898e+01,
+1.687898e+01, 1.703822e+01, 1.719745e+01, 1.735669e+01, 1.751592e+01, 1.767516e+01,
+1.783439e+01, 1.799363e+01, 1.815287e+01, 1.839172e+01, 1.855095e+01, 1.871019e+01,
+1.886943e+01, 1.918790e+01, 1.942675e+01, 1.966561e+01, 2.006369e+01, 2.054140e+01
+}; // 66
+
+
+const G4double G4HadronNucleonXsc::fKpProtonTotTkin[66] = {
+2.100000e+00, 2.180770e+00, 2.261540e+00, 2.396150e+00, 2.476920e+00, 2.557690e+00,
+2.557690e+00, 2.584620e+00, 2.638460e+00, 2.665380e+00, 2.719230e+00, 2.746150e+00,
+2.800000e+00, 2.853850e+00, 2.934620e+00, 3.042310e+00, 3.150000e+00, 3.311540e+00,
+3.446150e+00, 3.607690e+00, 3.930770e+00, 4.226920e+00, 4.361540e+00, 4.846150e+00,
+4.980770e+00, 5.088460e+00, 5.465380e+00, 5.653850e+00, 5.950000e+00, 6.084620e+00,
+6.246150e+00, 6.300000e+00, 6.380770e+00, 6.515380e+00, 6.730770e+00, 6.838460e+00,
+7.000000e+00, 7.161540e+00, 7.323080e+00, 7.457690e+00, 7.619230e+00, 7.780770e+00,
+7.915380e+00, 8.130770e+00, 8.265380e+00, 8.453850e+00, 8.642310e+00, 8.803850e+00,
+9.019230e+00, 9.234620e+00, 9.530770e+00, 9.773080e+00, 1.001538e+01, 1.017692e+01,
+1.033846e+01, 1.058077e+01, 1.082308e+01, 1.098462e+01, 1.114615e+01, 1.138846e+01,
+1.160385e+01, 1.173846e+01, 1.192692e+01, 1.216923e+01, 1.238461e+01, 1.257308e+01
+}; // 66
+
+const G4double G4HadronNucleonXsc::fKpNeutronTotXsc[75] = {
+3.980900e-01, 3.184700e-01, 3.184700e-01, 3.980900e-01, 3.980900e-01, 3.980900e-01,
+3.980900e-01, 3.980900e-01, 3.980900e-01, 4.777100e-01, 3.980900e-01, 3.980900e-01,
+4.777100e-01, 5.573200e-01, 1.035030e+00, 1.512740e+00, 2.149680e+00, 2.786620e+00,
+3.503180e+00, 4.219750e+00, 5.015920e+00, 5.652870e+00, 6.289810e+00, 7.245220e+00,
+8.121020e+00, 8.837580e+00, 9.633760e+00, 1.042994e+01, 1.114650e+01, 1.194268e+01,
+1.265924e+01, 1.329618e+01, 1.393312e+01, 1.449045e+01, 1.496815e+01, 1.552548e+01,
+1.592357e+01, 1.664013e+01, 1.727707e+01, 1.783439e+01, 1.831210e+01, 1.902866e+01,
+1.902866e+01, 1.878981e+01, 1.847134e+01, 1.831210e+01, 1.807325e+01, 1.791401e+01,
+1.783439e+01, 1.767516e+01, 1.759554e+01, 1.743631e+01, 1.743631e+01, 1.751592e+01,
+1.743631e+01, 1.735669e+01, 1.751592e+01, 1.759554e+01, 1.767516e+01, 1.783439e+01,
+1.783439e+01, 1.791401e+01, 1.815287e+01, 1.823248e+01, 1.847134e+01, 1.878981e+01,
+1.894905e+01, 1.902866e+01, 1.934713e+01, 1.966561e+01, 1.990446e+01, 2.014331e+01,
+2.030255e+01, 2.046178e+01, 2.085987e+01
+}; // 75
+
+const G4double G4HadronNucleonXsc::fKpNeutronTotTkin[75] = {
+2.692000e-02, 1.615400e-01, 2.961500e-01, 4.576900e-01, 6.461500e-01, 7.538500e-01,
+8.884600e-01, 1.103850e+00, 1.211540e+00, 1.400000e+00, 1.561540e+00, 1.776920e+00,
+1.992310e+00, 2.126920e+00, 2.342310e+00, 2.423080e+00, 2.557690e+00, 2.692310e+00,
+2.800000e+00, 2.988460e+00, 3.203850e+00, 3.365380e+00, 3.500000e+00, 3.688460e+00,
+3.850000e+00, 4.011540e+00, 4.173080e+00, 4.415380e+00, 4.630770e+00, 4.873080e+00,
+5.061540e+00, 5.276920e+00, 5.492310e+00, 5.707690e+00, 5.896150e+00, 6.030770e+00,
+6.138460e+00, 6.219230e+00, 6.273080e+00, 6.326920e+00, 6.407690e+00, 6.650000e+00,
+6.784620e+00, 7.026920e+00, 7.242310e+00, 7.350000e+00, 7.484620e+00, 7.619230e+00,
+7.807690e+00, 7.915380e+00, 8.050000e+00, 8.211540e+00, 8.453850e+00, 8.588460e+00,
+8.830770e+00, 9.073080e+00, 9.288460e+00, 9.476920e+00, 9.665380e+00, 9.826920e+00,
+1.004231e+01, 1.031154e+01, 1.052692e+01, 1.071538e+01, 1.095769e+01, 1.120000e+01,
+1.138846e+01, 1.155000e+01, 1.176538e+01, 1.190000e+01, 1.214231e+01, 1.222308e+01,
+1.238461e+01, 1.246538e+01, 1.265385e+01
+}; // 75
+
+const G4double G4HadronNucleonXsc::fKmProtonTotXsc[106] = {
+1.136585e+02, 9.749129e+01, 9.275262e+01, 8.885017e+01, 8.334146e+01, 7.955401e+01,
+7.504530e+01, 7.153658e+01, 6.858537e+01, 6.674913e+01, 6.525784e+01, 6.448781e+01,
+6.360279e+01, 6.255401e+01, 6.127526e+01, 6.032404e+01, 5.997910e+01, 5.443554e+01,
+5.376307e+01, 5.236934e+01, 5.113937e+01, 5.090941e+01, 4.967944e+01, 4.844948e+01,
+4.705575e+01, 4.638327e+01, 4.571080e+01, 4.475958e+01, 4.397213e+01, 4.257840e+01,
+4.102090e+01, 4.090592e+01, 3.906969e+01, 3.839721e+01, 3.756097e+01, 3.644599e+01,
+3.560976e+01, 3.533101e+01, 3.533101e+01, 3.644599e+01, 3.811847e+01, 3.839721e+01,
+3.979094e+01, 4.090592e+01, 4.257840e+01, 4.341463e+01, 4.425087e+01, 4.564460e+01,
+4.759582e+01, 4.703833e+01, 4.843206e+01, 4.787457e+01, 4.452962e+01, 4.202090e+01,
+4.034843e+01, 3.839721e+01, 3.616725e+01, 3.365854e+01, 3.170732e+01, 3.087108e+01,
+3.170732e+01, 3.254355e+01, 3.310104e+01, 3.254355e+01, 3.142857e+01, 3.059233e+01,
+2.947735e+01, 2.891986e+01, 2.836237e+01, 2.752613e+01, 2.696864e+01, 2.641115e+01,
+2.501742e+01, 2.473868e+01, 2.418118e+01, 2.362369e+01, 2.334495e+01, 2.278746e+01,
+2.250871e+01, 2.222997e+01, 2.167247e+01, 2.139373e+01, 2.139373e+01, 2.139373e+01,
+2.111498e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01,
+2.083624e+01, 2.055749e+01, 2.055749e+01, 2.055749e+01, 2.027875e+01, 2.000000e+01,
+2.055749e+01, 2.027875e+01, 2.083624e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01,
+2.083624e+01, 2.083624e+01, 2.139373e+01, 2.139373e+01
+}; // 106
+
+const G4double G4HadronNucleonXsc::fKmProtonTotTkin[106] = {
+4.017980e+00, 4.125840e+00, 4.179780e+00, 4.251690e+00, 4.287640e+00, 4.341570e+00,
+4.395510e+00, 4.467420e+00, 4.503370e+00, 4.575280e+00, 4.683150e+00, 4.737080e+00,
+4.773030e+00, 4.826970e+00, 4.880900e+00, 4.916850e+00, 4.952810e+00, 4.988760e+00,
+4.988760e+00, 5.006740e+00, 5.006740e+00, 5.042700e+00, 5.078650e+00, 5.114610e+00,
+5.132580e+00, 5.150560e+00, 5.186520e+00, 5.204490e+00, 5.276400e+00, 5.348310e+00,
+5.366290e+00, 5.384270e+00, 5.456180e+00, 5.564040e+00, 5.600000e+00, 5.671910e+00,
+5.743820e+00, 5.833710e+00, 5.905620e+00, 5.977530e+00, 6.085390e+00, 6.085390e+00,
+6.157300e+00, 6.175280e+00, 6.211240e+00, 6.229210e+00, 6.247190e+00, 6.337080e+00,
+6.391010e+00, 6.516850e+00, 6.462920e+00, 6.498880e+00, 6.570790e+00, 6.606740e+00,
+6.660670e+00, 6.678650e+00, 6.696630e+00, 6.732580e+00, 6.804490e+00, 6.876400e+00,
+6.948310e+00, 7.020220e+00, 7.074160e+00, 7.182020e+00, 7.235960e+00, 7.289890e+00,
+7.397750e+00, 7.523600e+00, 7.631460e+00, 7.757300e+00, 7.901120e+00, 8.062920e+00,
+8.260670e+00, 8.386520e+00, 8.530340e+00, 8.674160e+00, 8.817980e+00, 8.943820e+00,
+9.087640e+00, 9.267420e+00, 9.429210e+00, 9.573030e+00, 9.698880e+00, 9.896630e+00,
+1.002247e+01, 1.016629e+01, 1.031011e+01, 1.048989e+01, 1.063371e+01, 1.077753e+01,
+1.095730e+01, 1.108315e+01, 1.120899e+01, 1.135281e+01, 1.149663e+01, 1.162247e+01,
+1.174831e+01, 1.187416e+01, 1.200000e+01, 1.212584e+01, 1.221573e+01, 1.234157e+01,
+1.239551e+01, 1.250337e+01, 1.261124e+01, 1.273708e+01
+}; // 106
+
+const G4double G4HadronNucleonXsc::fKmNeutronTotXsc[68] = {
+2.621810e+01, 2.741123e+01, 2.868413e+01, 2.963889e+01, 3.067343e+01, 3.178759e+01,
+3.282148e+01, 3.417466e+01, 3.536778e+01, 3.552620e+01, 3.544576e+01, 3.496756e+01,
+3.433030e+01, 3.401166e+01, 3.313537e+01, 3.257772e+01, 3.178105e+01, 3.138264e+01,
+3.074553e+01, 2.970952e+01, 2.891301e+01, 2.827542e+01, 2.787700e+01, 2.715978e+01,
+2.660181e+01, 2.612394e+01, 2.564574e+01, 2.516721e+01, 2.421098e+01, 2.365235e+01,
+2.317366e+01, 2.261437e+01, 2.237389e+01, 2.205427e+01, 2.181395e+01, 2.165357e+01,
+2.149335e+01, 2.133297e+01, 2.109232e+01, 2.093128e+01, 2.069030e+01, 2.052992e+01,
+2.028927e+01, 2.012824e+01, 1.996737e+01, 1.996590e+01, 1.988530e+01, 1.964432e+01,
+1.948361e+01, 1.940236e+01, 1.940040e+01, 1.931882e+01, 1.947593e+01, 1.947429e+01,
+1.939320e+01, 1.939157e+01, 1.946922e+01, 1.962715e+01, 1.970481e+01, 1.970301e+01,
+1.993958e+01, 2.009669e+01, 2.025380e+01, 2.033178e+01, 2.049003e+01, 2.064747e+01,
+2.080540e+01, 2.096333e+01
+}; // 68
+
+const G4double G4HadronNucleonXsc::fKmNeutronTotTkin[68] = {
+5.708500e+00, 5.809560e+00, 5.896270e+00, 5.954120e+00, 5.997630e+00, 6.041160e+00,
+6.142160e+00, 6.171410e+00, 6.272470e+00, 6.344390e+00, 6.416230e+00, 6.459180e+00,
+6.487690e+00, 6.501940e+00, 6.544740e+00, 6.573280e+00, 6.616110e+00, 6.644710e+00,
+6.658840e+00, 6.744700e+00, 6.773150e+00, 6.830410e+00, 6.859010e+00, 6.916240e+00,
+6.973530e+00, 6.987730e+00, 7.030670e+00, 7.102360e+00, 7.173880e+00, 7.288660e+00,
+7.374720e+00, 7.547000e+00, 7.690650e+00, 7.791150e+00, 7.920420e+00, 8.020980e+00,
+8.107160e+00, 8.207720e+00, 8.365740e+00, 8.523790e+00, 8.710560e+00, 8.811110e+00,
+8.969140e+00, 9.127190e+00, 9.270860e+00, 9.400230e+00, 9.486440e+00, 9.673210e+00,
+9.802510e+00, 9.946220e+00, 1.011870e+01, 1.029116e+01, 1.047808e+01, 1.062181e+01,
+1.075114e+01, 1.089488e+01, 1.106739e+01, 1.118244e+01, 1.135496e+01, 1.151307e+01,
+1.171439e+01, 1.190130e+01, 1.208822e+01, 1.223199e+01, 1.231829e+01, 1.247646e+01,
+1.259150e+01, 1.270655e+01
+}; // 68
+
+
+
+
+
+//
+//
+///////////////////////////////////////////////////////////////////////////////////////
diff --git a/Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.hh b/Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.hh
new file mode 100644
index 000000000..fae9cb799
--- /dev/null
+++ b/Geant4/G4processes/srcnew/cross_sections/G4HadronNucleonXsc.hh
@@ -0,0 +1,162 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+// Calculation of the total, elastic and inelastic cross-sections
+// based on parametrisations of (proton, pion, kaon, photon) nucleon
+// cross-sections and the hadron-nucleous cross-section model in
+// the framework of Glauber-Gribov approach
+//
+//
+//
+// 14.03.07 V. Grichine - first implementation
+// 04.11.11 V. Grichine - update for kaon-(p,n) xsc, vector spline
+// 21.02.12 V. Grichine - update for pion-(p,n) xsc, NS fit++, vector spline
+//
+//
+
+#ifndef G4HadronNucleonXsc_h
+#define G4HadronNucleonXsc_h
+
+#include "globals.hh"
+#include "G4Proton.hh"
+#include "G4Nucleus.hh"
+#include "G4LPhysicsFreeVector.hh"
+
+
+class G4ParticleDefinition;
+class G4DynamicParticle;
+
+class G4HadronNucleonXsc
+{
+public:
+
+ G4HadronNucleonXsc ();
+ virtual ~G4HadronNucleonXsc ();
+
+ virtual
+ G4bool IsApplicable(const G4DynamicParticle* aDP, const G4Element*);
+
+ virtual
+ G4bool IsIsoApplicable(const G4DynamicParticle* aDP, G4int Z, G4int A);
+
+ virtual
+ void DumpPhysicsTable(const G4ParticleDefinition&)
+ {G4cout << "G4HadronNucleonXsc: uses parametrisation"<<G4endl;}
+
+ void CrossSectionDescription(std::ostream&) const;
+
+ // Xsc parametrisations
+
+ G4double GetHadronNucleonXscEL(const G4DynamicParticle*, const G4ParticleDefinition*);
+
+ G4double GetHadronNucleonXscPDG(const G4DynamicParticle*, const G4ParticleDefinition*);
+
+ G4double GetHadronNucleonXscNS(const G4DynamicParticle*, const G4ParticleDefinition*);
+
+ G4double GetKaonNucleonXscGG(const G4DynamicParticle*, const G4ParticleDefinition*);
+
+ G4double GetHadronNucleonXscVU(const G4DynamicParticle*, const G4ParticleDefinition*);
+
+ // kinematics and set/get
+
+ G4double CalculateEcmValue ( const G4double , const G4double , const G4double );
+
+ G4double CalcMandelstamS( const G4double , const G4double , const G4double );
+
+ G4double GetCoulombBarrier(const G4DynamicParticle* aParticle, const G4ParticleDefinition* nucleon );
+
+ G4double GetTotalHadronNucleonXsc() { return fTotalXsc; };
+ G4double GetElasticHadronNucleonXsc() { return fElasticXsc; };
+ G4double GetInelasticHadronNucleonXsc(){ return fInelasticXsc; };
+
+ void InitialiseKaonNucleonTotXsc();
+
+ G4double GetKpProtonTotXscVector(G4double logEnergy){ return fKpProtonTotXscVector.Value(logEnergy); };
+ G4double GetKpNeutronTotXscVector(G4double logEnergy){ return fKpNeutronTotXscVector.Value(logEnergy); };
+ G4double GetKmProtonTotXscVector(G4double logEnergy){ return fKmProtonTotXscVector.Value(logEnergy); };
+ G4double GetKmNeutronTotXscVector(G4double logEnergy){ return fKmNeutronTotXscVector.Value(logEnergy); };
+
+private:
+
+// const G4double fUpperLimit;
+ const G4double fLowerLimit;
+
+ G4double fTotalXsc, fElasticXsc, fInelasticXsc;
+// G4double fHadronNucleonXsc;
+
+ // K-nucleon tot xsc (mb) fit data, std::log(Tkin(MeV))
+
+ static const G4double fKpProtonTotXsc[66];
+ static const G4double fKpProtonTotTkin[66];
+
+ static const G4double fKpNeutronTotXsc[75];
+ static const G4double fKpNeutronTotTkin[75];
+
+ static const G4double fKmProtonTotXsc[106];
+ static const G4double fKmProtonTotTkin[106];
+
+ static const G4double fKmNeutronTotXsc[68];
+ static const G4double fKmNeutronTotTkin[68];
+
+ G4LPhysicsFreeVector fKpProtonTotXscVector;
+ G4LPhysicsFreeVector fKpNeutronTotXscVector;
+ G4LPhysicsFreeVector fKmProtonTotXscVector;
+ G4LPhysicsFreeVector fKmNeutronTotXscVector;
+
+ G4ParticleDefinition* theGamma;
+ G4ParticleDefinition* theProton;
+ G4ParticleDefinition* theNeutron;
+ G4ParticleDefinition* theAProton;
+ G4ParticleDefinition* theANeutron;
+ G4ParticleDefinition* thePiPlus;
+ G4ParticleDefinition* thePiMinus;
+ G4ParticleDefinition* thePiZero;
+ G4ParticleDefinition* theKPlus;
+ G4ParticleDefinition* theKMinus;
+ G4ParticleDefinition* theK0S;
+ G4ParticleDefinition* theK0L;
+ G4ParticleDefinition* theL;
+ G4ParticleDefinition* theAntiL;
+ G4ParticleDefinition* theSPlus;
+ G4ParticleDefinition* theASPlus;
+ G4ParticleDefinition* theSMinus;
+ G4ParticleDefinition* theASMinus;
+ G4ParticleDefinition* theS0;
+ G4ParticleDefinition* theAS0;
+ G4ParticleDefinition* theXiMinus;
+ G4ParticleDefinition* theXi0;
+ G4ParticleDefinition* theAXiMinus;
+ G4ParticleDefinition* theAXi0;
+ G4ParticleDefinition* theOmega;
+ G4ParticleDefinition* theAOmega;
+ G4ParticleDefinition* theD;
+ G4ParticleDefinition* theT;
+ G4ParticleDefinition* theA;
+ G4ParticleDefinition* theHe3;
+
+};
+
+
+#endif
diff --git a/Geant4Sys/cmt/requirements b/Geant4Sys/cmt/requirements
index b91f83dd2..21a973b7f 100755
--- a/Geant4Sys/cmt/requirements
+++ b/Geant4Sys/cmt/requirements
@@ -1,9 +1,9 @@
package Geant4Sys
-version v96r4
+version v96r4p1
branches cmt doc
-use G4config v96r1p0 Geant4
+use G4config v96r3p0 Geant4
# =============================================================================
# =========== global libraries ================================================
@@ -13,7 +13,7 @@ use G4intercoms v4r1p1 Geant4
use G4particles v6r2p1 Geant4
use G4track v5r1p1 Geant4
use G4geometry v6r2p1 Geant4
-use G4processes v9r0 Geant4
+use G4processes v9r1 Geant4
use G4physics_lists v1r3 Geant4
use G4tracking v6r1p1 Geant4
use G4global v5r1p1 Geant4
@@ -46,11 +46,11 @@ use G4UIbasic v4r4 Geant4
use G4UIGAG v4r3 Geant4
# LHCb additional package (eg. for MC11 G4LHCblists)
-use G4LHCblists v3r0 Geant4
+use G4LHCblists v3r2 Geant4
# G4 examples package - for checks
use G4analysis v1r0 Geant4
-use G4examples v5r0 Geant4
+use G4examples v5r1 Geant4
#=============================================================================
# ======== the end ============================================================
diff --git a/Geant4Sys/doc/release.notes b/Geant4Sys/doc/release.notes
index 8145f1c36..5e9bf2bcc 100755
--- a/Geant4Sys/doc/release.notes
+++ b/Geant4Sys/doc/release.notes
@@ -3,6 +3,21 @@ Package : Geant4Sys
Package manager(s) : Gloria Corti, Hubert DeGaudenzi, Nigel Watson
Purpose : LHCb build using cmt of Geant4
+!=============================================================================
+!</PRE><H1><A NAME=v96r5>2015-05-11 Geant4Sys v96r5</A></H1><PRE>
+
+! ======================= G4config v96r4p1 2015-05-11 ========================
+ - Updates include:
+ - Geant4/G4config v96r3p0
+ - Geant4/G4processes v9r1 - new kaon cross-sections
+ - Geant4/G4examples v5r1
+ - Geant4/G4LHCblists v3r2
+
+
+!=============================================================================
+!</PRE><H1><A NAME=v95r2p8>2015-02-04 Geant4Sys v95r2p8</A></H1><PRE>
+
+
!=============================================================================
!</PRE><H1><A NAME=v96r4>2015-03-06 Geant4Sys v96r4</A></H1><PRE>
diff --git a/cmt/project.cmt b/cmt/project.cmt
index 2874c13cf..ee0340997 100755
--- a/cmt/project.cmt
+++ b/cmt/project.cmt
@@ -1,6 +1,6 @@
project GEANT4
-use GAUDI GAUDI_v25r6p1
+use GAUDI GAUDI_v26r1
use PARAM
build_strategy with_installarea
diff --git a/toolchain.cmake b/toolchain.cmake
index 9ac10a775..eb33c3b9a 100644
--- a/toolchain.cmake
+++ b/toolchain.cmake
@@ -1,5 +1,5 @@
# Special wrapper to load the declared version of the heptools toolchain.
-set(heptools_version 71root6)
+set(heptools_version 74root6)
find_file(use_heptools_module UseHEPTools.cmake
HINTS ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
--
GitLab