diff --git a/Geant4/G4LHCblists/G4LHCblists/G4EmStandardPhysics_option1NoApplyCuts.hh b/Geant4/G4LHCblists/G4LHCblists/G4EmStandardPhysics_option1NoApplyCuts.hh
index 113c45b24e088d820abba5955ccd0e327a8df41f..5ffc74344f7c41f11b2e3809fbe7aff574863462 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 61303945935597c30e90f1f92341e0ea3f4a3e0d..d55a0550343bdf2d62558cdc6dc958a1c00aa05a 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 d7820fb03d6caa7940547e582986783225e1b3ef..9eea2a83135d2e99c5bc9998b38890427da16298 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 01c7250bb5b8d5558800c8b8e3b31a694b030edf..60de452d2cf7bbdddb9916d24a32c3635951551c 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 2e5999ae7c6b8de367892ebc9700b4a040be6f1f..015f7144fc478cfc8143f72dd736f9a93d5184f4 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 b10fa70c2376d9a020b1d3dd8d95ef75f600523d..2fa4a1725c7923fdc03dd3390448c0f97d97c9df 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 dd1e07487ac40c0dc003124e881be39a48f24a1c..86e2511f37c07ac983a6035639a9b339247fc5de 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 a61926b7cf542ab65c41dfc7ee2e71a422639d6b..74e139ef29ac98e99866c7886ce181c4ec1164c3 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 aade440e3d26de57316ce79d352e09aa3f948829..369912b08c532c0a79b565c518bd01758bc0a136 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 05ca94a899c8368d23ec13e5ce1de45d242a3ba2..8705bc139d891a3b454840862e1bebe6d1f5419c 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 e1c816c18365873f45e22c0b81767795fa098292..768bace4e943825db10ce47314b4abc6b953ff54 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 572b272fd088f5cf9a114f4305b8a0b95ad1b381..931fcaaa386aca0d6c84fe5e72f5f2cdfd8a5580 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 b4ecbab389fa0b6573c82f8a03347d212850f6f0..cdd282be71c69e1077f09338149f281c4486abfd 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 5c7c34fe8540e69907484453796ea895e8c9d75d..5581ca76eb9d45acdf34b54b42bf9952006e6dc6 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 cafa24ce9dbc94df780dab1f7514fbe17e2f4692..9748795fa7950ccfaeb930d03467653bc5f3204b 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 0000000000000000000000000000000000000000..c86cea1997bb547714aab56aec58145617dd57ea
--- /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 34a1810229d4c35d2e3d8c08e9677b547fc13f23..df8c0c2c9b102d8b6d0096e3549ff209dfc69d4b 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
Binary files a/Geant4/G4examples/extended/electromagnetic/G4TestEm3/scripts/OUTPUT_3.000 and b/Geant4/G4examples/extended/electromagnetic/G4TestEm3/scripts/OUTPUT_3.000 differ
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 0000000000000000000000000000000000000000..20269737d5ad73fb1ef15a377d5da7b09c9ce242
--- /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 0000000000000000000000000000000000000000..6d1adb9a28d31479f6b2901b87cc0b97e3505ee1
--- /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 0000000000000000000000000000000000000000..db3bc550155706e46415813eddccf7483ae8cb21
--- /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 0000000000000000000000000000000000000000..f48c70b55eb81e81b4a5de1d6d74665f8473b721
--- /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 8d31eb18fda2df09c18ea693e4cb83051a6bbb72..464a6b549c04799c589fd042b4e9d523bd4ea645 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 2526f55e91fde9e6ab658ef55afd7da62644637d..eaabf906508aeca1823246ca0e7a7a171d06e6f7 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 bc3cd26a1c3ea2b1bfe8c5cffbc71f8970dcc061..5b1eadf273dc3bfbe49cd1c231a5fab54da8ea8e 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 0000000000000000000000000000000000000000..cf4b7cc05034621c0564c942f5e209e8d578dfbc
--- /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 0000000000000000000000000000000000000000..e633db1bf5265498355f5b200ba2ed6203699d0f
--- /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 0000000000000000000000000000000000000000..03e14623ce83d8dca6ead90eb104ce092aa6cd97
--- /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 0000000000000000000000000000000000000000..fae9cb79915c7a3eab7a1699230e37e8f083873f
--- /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 b91f83dd2bc1b028216cf664bd86be5ac02d5afc..21a973b7febe924653216ba4ba92de490733ebbb 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 8145f1c3692fdb129e6cadbfec96a701f929ac5d..5e9bf2bcced2954b0722791ffaf2868e7f9fce1d 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 2874c13cf195188d730bdb9a455c782f71d321de..ee0340997679be1e30528d341ed31a175ee78be6 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 9ac10a77543e64d6a3e6491ff92f692e7c993e7b..eb33c3b9aa89ac4df8701192042533f804716e80 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)