diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 49e2c11b77e2704c84297600a21f5e1933594227..1b428066812c63c2a1cad080e0010446fb69ac08 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -211,13 +211,15 @@ int main() {
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
+ // fraction of oxygen
+ const float fox = 0.20946;
using MyHomogeneousModel =
corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
corsika::environment::NuclearComposition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Oxygen},
- std::vector<float>{1.}));
+ std::vector<corsika::particles::Code>{corsika::particles::Code::Nitrogen,corsika::particles::Code::Oxygen},
+ std::vector<float>{(float)1.-fox, fox}));
universe.AddChild(std::move(theMedium));
@@ -226,7 +228,7 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
- corsika::process::sibyll::Interaction sibyll;
+ corsika::process::sibyll::Interaction sibyll(env);
corsika::process::sibyll::Decay decay;
ProcessCut cut(8_GeV);
corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
@@ -270,6 +272,8 @@ int main() {
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
+ const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): "
- << (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy()) / 1_GeV << endl;
+ << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1. ) * 100 << endl;
}
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index c4083f6f391babff17c107475b287eef9fbac812..a476b89e8a7ada3f41152a9c0166aafa185a4d3f 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -22,6 +22,7 @@ target_link_libraries (
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
+ CORSIKArandom
)
target_include_directories (
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index cc612d044cb8c1d3f2d96ec081894806e6ab6dad..4c9a03f366a52c68b96e3ffe90bca8421b8f7810 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -17,6 +17,7 @@
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/random/RNGManager.h>
/**
* a homogeneous medium
@@ -41,6 +42,29 @@ namespace corsika::environment {
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+ corsika::particles::Code const& GetTarget( std::vector<corsika::units::si::CrossSectionType> &sigma) const override {
+ using namespace corsika::units::si;
+ int i=-1;
+ corsika::units::si::CrossSectionType total_weighted_sigma = 0._mbarn;
+ std::vector<float> fractions;
+ for(auto w: fNuclComp.GetFractions() ){
+ i++;
+ std::cout << "HomogeneousMedium: fraction: " << w << std::endl;
+ total_weighted_sigma += w * sigma[i];
+ fractions.push_back( w * sigma[i] / 1_mbarn );
+ }
+
+ for(auto f: fractions){
+ f = f / ( total_weighted_sigma / 1_mbarn );
+ std::cout << "HomogeneousMedium: reweighted fraction: " << f << std::endl;
+ }
+ std::discrete_distribution channelDist( fractions.begin(), fractions.end() );
+ static corsika::random::RNG& kRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ const int ichannel = channelDist(kRNG);
+ return fNuclComp.GetComponents()[ichannel];
+ }
+
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::LengthType pTo) const override {
@@ -52,7 +76,7 @@ namespace corsika::environment {
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
- }
+ }
};
} // namespace corsika::environment
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index ddae0564200a1a97428a0791f8333feea8c8932b..24393e805eb46a90febe4ac568c2949334f88bfe 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -38,6 +38,8 @@ namespace corsika::environment {
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
+
+ virtual corsika::particles::Code const& GetTarget( std::vector<corsika::units::si::CrossSectionType> &) const = 0;
};
} // namespace corsika::environment
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index ab8f9b3e094a42b00a63899442774d6b42c1b1f7..bf89db756b5aaca2be901596c1015b5d724245e7 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -69,6 +69,7 @@ target_link_libraries (
CORSIKAunits
CORSIKAthirdparty
CORSIKAgeometry
+ CORSIKAenvironment
)
target_include_directories (
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fac0dd578b77cfea776b1b942b6a42b9e899f809..b8709216a75bd7a827fea1b063520fa615805553 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -21,6 +21,8 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
namespace corsika::process::sibyll {
@@ -28,9 +30,8 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
-
public:
- Interaction() {}
+ Interaction(corsika::environment::Environment const& env) : fEnvironment(env) { }
~Interaction() {
std::cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount
<< std::endl;
@@ -44,11 +45,38 @@ namespace corsika::process::sibyll {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm");
-
+
// initialize Sibyll
sibyll_ini_();
}
+ std::tuple<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];
+ const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
+ const double dEcm = CoMenergy / 1_GeV;
+ 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
+ // 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 );
+ }
+ }
+
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) {
@@ -66,60 +94,65 @@ namespace corsika::process::sibyll {
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
- const int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- // target nuclei: A < 18
- // FOR NOW: assume target is oxygen
- const int kTarget = corsika::particles::Oxygen::GetNucleusA();
-
+
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
- HEPEnergyType const Elab = p.GetEnergy();
- HEPEnergyType const Etot = Elab + nucleon_mass;
- MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ corsika::particles::Neutron::GetMass());
+
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+
+ // total momentum and energy
+ HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+ MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
Ptot += p.GetMomentum();
Ptot += pTarget;
// calculate cm. energy
- const HEPEnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
- const double Ecm = sqs / 1_GeV;
+ const HEPEnergyType ECoM = sqrt(Elab * Elab - Ptot.squaredNorm());
std::cout << "Interaction: LambdaInt: \n"
- << " input energy: " << Elab / 1_GeV << endl
- << " beam can interact:" << kBeam << endl
- << " beam XS code:" << kBeam << endl
- << " beam pid:" << p.GetPID() << endl
- << " target mass number:" << kTarget << std::endl;
-
- if (kInteraction && Elab >= 8.5_GeV && sqs >= 10_GeV) {
-
- double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
-
- if (kTarget == 1)
- sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
- else
- sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
-
- std::cout << "Interaction: "
- << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- si::CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "Interaction: "
- << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
-
- std::cout << "Interaction: "
- << "nucleon mass " << nucleon_mass << std::endl;
- // calculate interaction length in medium
- GrammageType int_length = kTarget * corsika::units::constants::u / sig;
- std::cout << "Interaction: "
- << "interaction length (g/cm2): " << int_length << std::endl;
+ << " input energy: " << p.GetEnergy() / 1_GeV << endl
+ << " beam can interact:" << kInteraction << endl
+ << " beam pid:" << p.GetPID() << endl;
+
+ if (kInteraction && Elab >= 8.5_GeV && ECoM >= 10_GeV) {
+
+ // get target from environment
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+ */
+ const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
+ // determine average interaction length
+ // weighted sum
+ int i =-1;
+ double avgTargetMassNumber = 0.;
+ si::CrossSectionType weightedProdCrossSection = 0_mbarn;
+ // get weights of components from environment/medium
+ const auto w = mediumComposition.GetFractions();
+ // loop over components in medium
+ for (auto targetId: mediumComposition.GetComponents() ){
+ i++;
+ cout << "Interaction: get interaction length for target: " << targetId << endl;
+
+ auto const [ productionCrossSection, numberOfNucleons ] = GetCrossSection( corsikaBeamId, targetId, ECoM);
+
+ std::cout << "Interaction: "
+ << " IntLength: sibyll return (mb): " << productionCrossSection / 1_mbarn << std::endl;
+ weightedProdCrossSection += w[i] * productionCrossSection;
+ avgTargetMassNumber += w[i] * numberOfNucleons;
+ }
+ cout << "Interaction: "
+ << "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mbarn
+ << endl;
+
+ // calculate interaction length in medium
+#warning check interaction length units
+ GrammageType int_length = avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
+ std::cout << "Interaction: "
+ << "interaction length (g/cm2): " << int_length / ( 0.001_kg ) * 1_cm * 1_cm << std::endl;
return int_length;
}
@@ -137,10 +170,11 @@ namespace corsika::process::sibyll {
using std::cout;
using std::endl;
+ const auto corsikaBeamId = p.GetPID();
cout << "ProcessSibyll: "
- << "DoInteraction: " << p.GetPID() << " interaction? "
- << process::sibyll::CanInteract(p.GetPID()) << endl;
- if (process::sibyll::CanInteract(p.GetPID())) {
+ << "DoInteraction: " << corsikaBeamId << " interaction? "
+ << process::sibyll::CanInteract( corsikaBeamId ) << endl;
+ if (process::sibyll::CanInteract(corsikaBeamId) ) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -171,18 +205,8 @@ namespace corsika::process::sibyll {
// double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
const CoordinateSystem tempCS = rootCS.rotate(zAxis, phi);
const CoordinateSystem sibyllCS = tempCS.rotate(yAxis, theta);
-
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
-
- here we need: GetTargetMassNumber() or GetTargetPID()??
- GetTargetMomentum() (zero in EAS)
- */
- // FOR NOW: set target to oxygen
- const int kTarget = corsika::particles::Oxygen::
- GetNucleusA(); // env.GetTargetParticle().GetPID();
+
+
// FOR NOW: target is always at rest
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
@@ -213,6 +237,39 @@ namespace corsika::process::sibyll {
MomentumVector Ptot = p.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+
+
+ // sample target mass number
+ const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
+ const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
+ // get cross sections for target materials
+ /*
+ Here we read the cross section from the interaction model again,
+ should be passed from GetInteractionLength if possible
+ */
+#warning reading interaction cross section again, should not be necessary
+ std::vector<si::CrossSectionType> cross_section_of_components;
+ for(auto targetId: mediumComposition.GetComponents() ){
+ const auto [sigProd, nNuc] = GetCrossSection( corsikaBeamId, targetId, Ecm);
+ cross_section_of_components.push_back( sigProd );
+ }
+
+ const auto targetCode = currentNode->GetModelProperties().GetTarget( cross_section_of_components);
+ cout << "Interaction: target selected: " << targetCode << endl;
+ /*
+ FOR NOW: allow nuclei with A<18 or protons only.
+ when medium composition becomes more complex, approximations will have to be allowed
+ air in atmosphere also contains some Argon.
+ */
+ int targetSibCode = -1;
+ if( IsNucleus(targetCode))
+ targetSibCode = GetNucleusA( targetCode );
+ if( targetCode == corsika::particles::Proton::GetCode())
+ targetSibCode = 1;
+ cout << "Interaction: sibyll code: " << targetSibCode << endl;
+ if(targetSibCode>18||targetSibCode<1)
+ throw std::runtime_error("Sibyll target outside range. Only nuclei with A<18 or protons are allowed.");
+
/*
get transformation between Stack-frame and SibStack-frame
for EAS Stack-frame is lab. frame, could be different for CRMC-mode
@@ -246,7 +303,7 @@ namespace corsika::process::sibyll {
cout << "Interaction: new boost: ebeam com: " << eProjectileCoM / 1_GeV << endl
<< "Interaction: new boost: pbeam com: " << pProjectileCoM / 1_GeV << endl;
- int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
std::cout << "Interaction: "
<< " DoInteraction: E(GeV):" << E / 1_GeV
@@ -260,7 +317,7 @@ namespace corsika::process::sibyll {
// Sibyll does not know about units..
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
- sibyll_(kBeam, kTarget, sqs);
+ sibyll_(kBeam, targetSibCode, sqs);
// running decays
// setTrackedParticlesStable();
decsib_();
@@ -334,6 +391,10 @@ namespace corsika::process::sibyll {
}
return process::EProcessReturn::eOk;
}
+
+ private:
+ corsika::environment::Environment const& fEnvironment;
+
};
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 3d7d425d9c63aaacf5ffcef5723dfac730283b7e..eb3259e792aa86899f6569b7a1247dce8acd4cf6 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -74,25 +74,48 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/particles/ParticleProperties.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
using namespace corsika::units::si;
using namespace corsika::units;
TEST_CASE("SibyllInterface", "[processes]") {
- auto const& cs =
- geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ // setup environment, geometry
+ corsika::environment::Environment env;
+ auto& universe = *(env.GetUniverse());
+
+ auto theMedium = corsika::environment::Environment::CreateNode<geometry::Sphere>(
+ geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel =
+ corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ corsika::environment::NuclearComposition(
+ std::vector<corsika::particles::Code>{corsika::particles::Code::Oxygen},
+ std::vector<float>{1.}));
+
+ universe.AddChild(std::move(theMedium));
+
+ const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
cs, 0_m / second, 0_m / second, 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
+
SECTION("InteractionInterface") {
setup::Stack stack;
auto particle = stack.NewParticle();
- Interaction model;
+ Interaction model(env);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =