From bb8337012004828f9f57d5028174a9f481d1c550 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sun, 30 Dec 2018 19:45:08 +0000
Subject: [PATCH] moved target sampling from sibyll process to medium
---
Documentation/Examples/cascade_example.cc | 6 ++--
Environment/CMakeLists.txt | 1 +
Environment/HomogeneousMedium.h | 26 +++++++++++++++-
Environment/IMediumModel.h | 2 ++
Processes/Sibyll/Interaction.h | 36 ++++++-----------------
5 files changed, 41 insertions(+), 30 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 61d157663..5fc6942dc 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -239,7 +239,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const hep::EnergyType E0 = 100_GeV;
+ const hep::EnergyType E0 = 100_TeV;
double theta = 0.;
double phi = 0.;
{
@@ -269,6 +269,8 @@ int main() {
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
+ const hep::EnergyType 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 c4083f6f3..a476b89e8 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 cc612d044..4c9a03f36 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 ddae05642..24393e805 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/Interaction.h b/Processes/Sibyll/Interaction.h
index fc3590d74..60553ef42 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -248,37 +248,19 @@ namespace corsika::process::sibyll {
// sample target mass number
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
- // get cross sections and nucleon number for target materials
- std::vector<si::CrossSectionType> cross_section_components;
+ // 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_components.push_back( sigProd );
+ cross_section_of_components.push_back( sigProd );
}
- // this routine could be moved to the environment as GetTarget( std::vector<si::CrossSectionType> )
- const auto sample_target = [](const corsika::environment::NuclearComposition& comp, const std::vector<si::CrossSectionType> & sigma_comp )
- {
- int i=-1;
- si::CrossSectionType total_weighted_sigma = 0._mbarn;
- std::vector<float> fractions;
- for(auto w: comp.GetFractions() ){
- i++;
- cout << "fraction: " << w << endl;
- total_weighted_sigma += w * sigma_comp[i];
- fractions.push_back( w * sigma_comp[i] / 1_mbarn );
- }
-
- for(auto f: fractions){
- f = f / ( total_weighted_sigma / 1_mbarn );
- cout << "reweighted fraction: " << f << 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 comp.GetComponents()[ichannel];
- };
- const auto targetCode = sample_target( mediumComposition, cross_section_components );
+ 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.
--
GitLab