diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 611606b56f9295476f8952b58b1251211e6d1bd3..7c3864b3a3ef036ec1b767c1d4256d66941fbeb3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -243,9 +243,10 @@ int main() {
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
corsika::process::sibyll::Interaction sibyll(env);
- corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
+ // corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
+ corsika::process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
corsika::process::sibyll::Decay decay;
- ProcessCut cut(200_GeV);
+ ProcessCut cut(20_GeV);
// corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// corsika::process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env);
@@ -282,6 +283,7 @@ int main() {
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
particle.SetEnergy(E0);
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index b403d09cfa79f5cdcc58cca8ff68ed54d1b65b23..f34c48926fdffcbfec6936dd6d1cbab324580408 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -31,7 +31,8 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
-
+ bool fInitialized = false;
+
public:
Interaction(corsika::environment::Environment const& env)
: fEnvironment(env) {}
@@ -47,32 +48,38 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize Sibyll
- sibyll_ini_();
+ if(!fInitialized){
+ sibyll_ini_();
+ fInitialized = true;
+ }
}
- std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
+ bool WasInitialized(){ return fInitialized;}
+
+ std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType, int> GetCrossSection(
const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy) {
using namespace corsika::units::si;
- double sigProd, dummy, dum1, dum2, dum3, dum4;
+ double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
const double dEcm = CoMenergy / 1_GeV;
+ int iTarget = -1;
if (corsika::particles::IsNucleus(TargetId)) {
- const int iTarget = corsika::particles::GetNucleusA(TargetId);
+ iTarget = corsika::particles::GetNucleusA(TargetId);
if (iTarget > 18 || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
- sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
- return std::make_tuple(sigProd * 1_mbarn, iTarget);
+ sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == corsika::particles::Proton::GetCode()) {
- sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
- return std::make_tuple(sigProd * 1_mbarn, 1);
+ sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+ iTarget = 1;
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
- return std::make_tuple(sigProd * 1_mbarn, 0);
+ sigEla = std::numeric_limits<double>::infinity();
}
+ return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
}
template <typename Particle, typename Track>
@@ -139,7 +146,7 @@ namespace corsika::process::sibyll {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
- auto const [productionCrossSection, numberOfNucleons] =
+ auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
std::cout << "Interaction: "
@@ -266,7 +273,7 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- const auto [sigProd, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ const auto [sigProd, sigEla, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
}
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 7ee1fb9201b3d244e1a9f4df62414d1e0842b76a..8ccce710b2c5a725815fb90e261f54a193e5a7ac 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -17,9 +17,6 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/process/sibyll/SibStack.h>
-#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/process/sibyll/nuclib.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
@@ -33,10 +30,12 @@ namespace corsika::process::sibyll {
int fNucCount = 0;
public:
- NuclearInteraction(corsika::environment::Environment const& env)
- : fEnvironment(env) {}
+ NuclearInteraction(corsika::environment::Environment const& env, corsika::process::sibyll::Interaction& hadint)
+ : fEnvironment(env),
+ fHadronicInteraction(hadint)
+ {}
~NuclearInteraction() {
- std::cout << "Sibyll::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
+ std::cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
<< std::endl;
}
@@ -47,45 +46,28 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize hadronic interaction module
- //sibyll_ini_();
+ // TODO: save to run multiple initializations?
+ if(!fHadronicInteraction.WasInitialized())
+ fHadronicInteraction.Init();
// initialize nuclib
+ // TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
}
- std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
+ // TODO: remove number of nucleons, avg target mass is available in environment
+ std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType, int> GetCrossSection(
const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy) {
using namespace corsika::units::si;
- double sigProd, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
+ double sigProd;
std::cout << "NuclearInteraction: GetCrossSection: called with: beamId= "<<BeamId
<< " TargetId= " << TargetId << " CoMenergy= " << CoMenergy / 1_GeV << std::endl;
if(!corsika::particles::IsNucleus(BeamId) ){
- // ask sibyll for hadron cross section
- // this is needed for sample target, which uses the nucleon cross section
- // TODO: generalize to parent hadronic interaction
- const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
- const double dEcm = CoMenergy / 1_GeV;
- // check target, hadron or nucleus?
- if (corsika::particles::IsNucleus(TargetId)) {
- const int iTarget = corsika::particles::GetNucleusA(TargetId);
- if (iTarget > 18 || iTarget == 0)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 are allowed.");
- sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
- return std::make_tuple(sigProd * 1_mbarn, iTarget);
- } else if (TargetId == corsika::particles::Proton::GetCode()) {
- sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
- return std::make_tuple(sigProd * 1_mbarn, 1);
- } else {
- throw std::runtime_error("GetCrossSection: no interaction in sibyll possible");
- // no interaction in sibyll possible, return infinite cross section? or throw?
- sigProd = std::numeric_limits<double>::infinity();
- return std::make_tuple(sigProd * 1_mbarn, 0);
- }
+ // ask hadronic interaction model for hadron cross section
+ return fHadronicInteraction.GetCrossSection( BeamId, TargetId, CoMenergy);
}
// use nuclib to calc. nuclear cross sections
@@ -101,7 +83,9 @@ namespace corsika::process::sibyll {
const double dEcm = CoMenergy / 1_GeV;
if(dEcm<10.)
throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
-
+
+ // TODO: these limitations are still sibyll specific.
+ // available target nuclei depends on the hadronic interaction model and the initialization
if (corsika::particles::IsNucleus(TargetId)) {
const int iTarget = corsika::particles::GetNucleusA(TargetId);
if (iTarget > 18 || iTarget == 0)
@@ -110,8 +94,9 @@ namespace corsika::process::sibyll {
std::cout<< "NuclearInteraction: calling signuc.." << std::endl;
signuc_(iBeam, dEcm, sigProd);
std::cout << "cross section: " << sigProd << std::endl;
- return std::make_tuple(sigProd * 1_mbarn, iTarget);
+ return std::make_tuple(sigProd * 1_mbarn, 0_mbarn, iTarget);
}
+ return std::make_tuple(std::numeric_limits<double>::infinity()* 1_mbarn, std::numeric_limits<double>::infinity()* 1_mbarn, 0);
}
template <typename Particle, typename Track>
@@ -134,7 +119,6 @@ namespace corsika::process::sibyll {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
- // beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
@@ -180,12 +164,12 @@ namespace corsika::process::sibyll {
for (auto const targetId : mediumComposition.GetComponents()) {
i++;
cout << "NuclearInteraction: get interaction length for target: " << targetId << endl;
-
- auto const [productionCrossSection, numberOfNucleons] =
+ // TODO: remove number of nucleons, avg target mass is available in environment
+ auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
-
+
std::cout << "NuclearInteraction: "
- << "IntLength: sibyll return (mb): "
+ << "IntLength: nuclib return (mb): "
<< productionCrossSection / 1_mbarn << std::endl;
weightedProdCrossSection += w[i] * productionCrossSection;
avgTargetMassNumber += w[i] * numberOfNucleons;
@@ -237,7 +221,7 @@ namespace corsika::process::sibyll {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- // position and time of interaction, not used in Sibyll
+ // position and time of interaction, not used in NUCLIB
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
@@ -272,7 +256,6 @@ namespace corsika::process::sibyll {
// define target
// always a nucleon
- // for Sibyll is always a single nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
// target is always at rest
@@ -326,7 +309,7 @@ namespace corsika::process::sibyll {
auto const targetId = compVec[i];
cout << "target component: " << targetId << endl;
cout << "beam id: " << targetId << endl;
- const auto [sigProd, nNuc] = GetCrossSection(beamId,targetId,EcmNN);
+ const auto [sigProd, sigEla, nNuc] = GetCrossSection(beamId,targetId,EcmNN);
cross_section_of_components[i] = sigProd;
}
@@ -341,7 +324,7 @@ namespace corsika::process::sibyll {
int kATarget = -1;
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
if (targetCode == corsika::particles::Proton::GetCode()) kATarget = 1;
- cout << "NuclearInteraction: sibyll target code: " << kATarget << endl;
+ cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
if (kATarget > 18 || kATarget < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
@@ -352,19 +335,12 @@ namespace corsika::process::sibyll {
cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "<< endl;
// get nucleon-nucleon cross section
// (needed to determine number of scatterings)
- //
- // TODO: this is an explicit dependence on sibyll
- // should be changed to a call to GetCrossSection() of the had. int. implementation
- // this also means GetCrossSection has to return the elastic cross section as well
- double sigProd, sigEla, dum1, dum2, dum3;
- double dumdif[3];
- const double sqsnuc = EcmNN / 1_GeV;
- const int iBeam = 1;
- // read hadron-proton cross section table (input: hadron id, CoMenergy, output: cross sections)
- sib_sigma_hp_(iBeam, sqsnuc, dum1, sigEla, sigProd, dumdif, dum2, dum3);
-
+ const auto protonId = corsika::particles::Proton::GetCode();
+ const auto [prodCrossSection, elaCrossSection, dum] = fHadronicInteraction.GetCrossSection(protonId,protonId,EcmNN);
+ const double sigProd = prodCrossSection / 1_mbarn;
+ const double sigEla = elaCrossSection / 1_mbarn;
// sample number of interactions (only input variables, output in common cnucms)
- // nuclear multiple scattering according to glauber (r.i.p.)
+ // nuclear multiple scattering according to glauber (r.i.p.)
int_nuc_( kATarget, kAProj, sigProd, sigEla);
cout << "number of nucleons in target : " << kATarget << endl
@@ -405,10 +381,6 @@ namespace corsika::process::sibyll {
<< " pz=" << fragments_.ppp[j][2]
<< endl;
- // bookeeping accross nucleon-nucleon interactions
- MomentumVector Plab_all(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_all = 0_GeV;
-
auto Nucleus = []( int iA ){
// int znew = iA / 2.15 + 0.7;
@@ -470,12 +442,11 @@ namespace corsika::process::sibyll {
pnew.SetMomentum(Plab.GetSpaceLikeComponents());
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
-
- Plab_all += Plab.GetSpaceLikeComponents();
- Elab_all += Plab.GetTimeLikeComponent();
}
-
+
// add elastic nucleons to corsika stack
+ // TODO: the elastic interaction could be external like the inelastic interaction,
+ // e.g. use existing ElasticModel
for(int j=0; j<nElasticNucleons; ++j){
auto pnew = s.NewParticle();
// TODO: sample proton or neutron
@@ -492,75 +463,39 @@ namespace corsika::process::sibyll {
pnew.SetMomentum(Plab.GetSpaceLikeComponents());
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
-
- Plab_all += Plab.GetSpaceLikeComponents();
- Elab_all += Plab.GetTimeLikeComponent();
}
- // add inelastic interactions
- for(int j=0; j<nInelNucleons; ++j){
-
- // create nucleon-nucleus inelastic interaction
- // TODO: switch to DoInteraction() of had. interaction implementation
- // for now use SIBYLL directly
- const int kBeamCode = process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode());
- cout << "creating interaction no. "<< j << endl;
- sibyll_( kBeamCode, kATarget, sqsnuc );
- decsib_();
- // print final state
- int print_unit = 6;
- sib_list_(print_unit);
-
- // add particles from sibyll to stack
- // link to sibyll stack
- SibStack ss;
-
- MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
- for (auto& psib : ss) {
-
- // skip particles that have decayed in Sibyll
- if (psib.HasDecayed()) continue;
-
- // // transform energy to lab. frame
- auto const pCoM = psib.GetMomentum();
- HEPEnergyType const eCoM = psib.GetEnergy();
- auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
-
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- pnew.SetEnergy(Plab.GetTimeLikeComponent());
- pnew.SetMomentum(Plab.GetSpaceLikeComponents());
- pnew.SetPosition(pOrig);
- pnew.SetTime(tOrig);
-
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
- Ecm_final += psib.GetEnergy();
-
- Plab_all += Plab.GetSpaceLikeComponents();
- Elab_all += Plab.GetTimeLikeComponent();
- }
- cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV
- << endl
- << "Elab_final=" << Elab_final / 1_GeV
- << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
- << endl;
+ // create dummy stack, to store nucleon projectile
+ Stack inelasticNucleons;
+ // add one nucleon per inelastic interaction
+ for(int j=0; j<nInelNucleons; ++j){
+
+ auto NucleonProjectile = inelasticNucleons.NewParticle();
+ // TODO: sample neutron or proton
+ auto pCode = corsika::particles::Proton::GetCode();
+ NucleonProjectile.SetPID(pCode);
+ NucleonProjectile.SetMomentum( PprojNucLab.GetSpaceLikeComponents() );
+ NucleonProjectile.SetEnergy( PprojNucLab.GetTimeLikeComponent() );
+ NucleonProjectile.SetPosition( pOrig );
+ NucleonProjectile.SetTime( tOrig );
+
}
- cout << "accross all nucleon interactions: Etot lab: " << Elab_all / 1_GeV << endl
- << "accross all nucleon interactions: Ptot lab: " << (Plab_all / 1_GeV).GetComponents()
- << endl;
- // delete current particle
+ // process stack of inelastic nucleons
+ for(auto& nucl : inelasticNucleons)
+ // create inelastic nucleon-nucleon interaction
+ fHadronicInteraction.DoInteraction( nucl, s);
+
+ // delete parent particle
p.Delete();
- //throw std::runtime_error(" stop here");
+ // throw std::runtime_error(" stop here");
return process::EProcessReturn::eOk;
}
private:
corsika::environment::Environment const& fEnvironment;
+ corsika::process::sibyll::Interaction& fHadronicInteraction;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
};
diff --git a/Processes/Sibyll/sibyll2.3c.f b/Processes/Sibyll/sibyll2.3c.f
index 7358ef9a78125d726da9155a054a790ef59b357a..e1f07491be0e017d84ac1f809ea00301615f4356 100644
--- a/Processes/Sibyll/sibyll2.3c.f
+++ b/Processes/Sibyll/sibyll2.3c.f
@@ -5516,7 +5516,7 @@ c external type declarations
c local type declarations
DOUBLE PRECISION SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO,
- & SIGPROD,SIGBDIF,S_RNDM,S,PF,PB,PD,P0,P1,P2,R
+ & SIGPROD,SIGBDIF,SIGELA,S_RNDM,S,PF,PB,PD,P0,P1,P2,R
DIMENSION SIGDIF(3)
INTEGER K
SAVE
@@ -5534,7 +5534,7 @@ c read hadron-nucleon cross section from table
c distinguish between nuclear cross sections..
IF(IAFLG.eq.0)THEN
c if target is nucleus calc. hadron-nucleus cross section (slow)
- CALL SIB_SIGMA_HNUC(L,IA,SQS,SIGprod,SIGbdif)
+ CALL SIB_SIGMA_HNUC(L,IA,SQS,SIGprod,SIGbdif,SIGela)
ELSE
c if target is air read hadron-air cross section from table
CALL SIB_SIGMA_HAIR(L,SQS,SIGprod,SIGbdif)
@@ -14010,7 +14010,7 @@ c internal
END
C=======================================================================
- SUBROUTINE SIB_SIGMA_HNUC (L,IAT,SQS,SIGprod,SIGbdif)
+ SUBROUTINE SIB_SIGMA_HNUC (L,IAT,SQS,SIGprod,SIGbdif,SIGela)
C-----------------------------------------------------------------------
C calculate Hadron-nucleus cross sections
@@ -14061,7 +14061,7 @@ C--------------------------------------------------------------------
+ SIGQSD,SIGPPT,SIGPPEL,SIGPPSD,ITG
c external
- DOUBLE PRECISION SQS,SIGPROD,SIGBDIF
+ DOUBLE PRECISION SQS,SIGPROD,SIGBDIF,SIGELA
INTEGER L,IAT
c internal
@@ -14092,6 +14092,8 @@ C particle production cross section
SIGprod = SIGT-SIGQE
C quasi elastic + elastic singl. diff cross section
SIGbdif = SIGQSD
+c elastic cross section
+ SIGela = SIGel
if(ndebug.gt.0)THEN
WRITE(LUN,'(1X,A,3F8.2)')
& 'SIB_SIGMA_HNUC: SIGprod, SIGbdif, ALAM:',
diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
index 5d3a5635b8eba54164297cc0bc31126f701c2e81..0b7f69a0da7c631f22725e28d0ef22e373aac4d5 100644
--- a/Processes/Sibyll/sibyll2.3c.h
+++ b/Processes/Sibyll/sibyll2.3c.h
@@ -95,7 +95,7 @@ void decsib_();
// interaction length
// double fpni_(double&, int&);
-void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&);
+ void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&, double&);
void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
double&);
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 79316ba8fee006bfba759a62fd297f6be560098c..e5aa388d3090ba22ab5d80a6bda088cf422de628 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -133,7 +133,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
setup::Stack stack;
auto particle = stack.NewParticle();
- NuclearInteraction model(env);
+ Interaction hmodel(env);
+ NuclearInteraction model(env, hmodel);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =