From 824e9fc37b47a986ec283c1f9ec4ffe6e9c9e9a6 Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@gmx.de>
Date: Mon, 2 Dec 2019 15:41:33 +0100
Subject: [PATCH] First very limited experiment of the electro magnetic process
with proposal.
The proposal library has to be installed manually and corresponding cmake variables have to be exported.
---
Documentation/Examples/CMakeLists.txt | 157 ++++++++++--------
Documentation/Examples/cascade_example.cc | 15 +-
Documentation/Examples/proposal_example.cc | 179 +++++++++++++++++++++
Processes/CMakeLists.txt | 5 +-
Processes/ParticleCut/ParticleCut.cc | 13 +-
Processes/ParticleCut/ParticleCut.h | 4 +-
Processes/Proposal/CMakeLists.txt | 46 ++++++
Processes/Proposal/ContinuousProcess.cc | 106 ++++++++++++
Processes/Proposal/ContinuousProcess.h | 128 +++++++++++++++
Processes/Proposal/Interaction.cc | 113 +++++++++++++
Processes/Proposal/Interaction.h | 107 ++++++++++++
11 files changed, 796 insertions(+), 77 deletions(-)
create mode 100644 Documentation/Examples/proposal_example.cc
create mode 100644 Processes/Proposal/CMakeLists.txt
create mode 100644 Processes/Proposal/ContinuousProcess.cc
create mode 100644 Processes/Proposal/ContinuousProcess.h
create mode 100644 Processes/Proposal/Interaction.cc
create mode 100644 Processes/Proposal/Interaction.h
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 1592372c2..1c4a2e76c 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -14,72 +14,26 @@ target_link_libraries (stack_example SuperStupidStack CORSIKAunits)
# address sanitizer is making this example too slow, so we only do "undefined"
CORSIKA_ADD_EXAMPLE (boundary_example)
target_link_libraries (boundary_example
- SuperStupidStack
- CORSIKAunits
- CORSIKAlogging
- CORSIKArandom
- ProcessSibyll
- CORSIKAcascade
- ProcessTrackWriter
- ProcessParticleCut
- ProcessTrackingLine
- CORSIKAprocesses
- CORSIKAparticles
- CORSIKAgeometry
- CORSIKAenvironment
- CORSIKAprocesssequence
- )
-
-CORSIKA_ADD_EXAMPLE (cascade_example)
-target_link_libraries (cascade_example
- SuperStupidStack
- CORSIKAunits
- CORSIKAlogging
- CORSIKArandom
- ProcessSibyll
- CORSIKAcascade
- ProcessEnergyLoss
- ProcessTrackWriter
- ProcessStackInspector
- ProcessTrackingLine
- ProcessParticleCut
- ProcessHadronicElasticModel
- ProcessStackInspector
- CORSIKAprocesses
- CORSIKAcascade
- CORSIKAparticles
- CORSIKAgeometry
- CORSIKAenvironment
- CORSIKAprocesssequence
- )
-
-if (Pythia8_FOUND)
- CORSIKA_ADD_EXAMPLE (cascade_proton_example)
- target_link_libraries (cascade_proton_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia8
+ ProcessProposal
CORSIKAcascade
- ProcessEnergyLoss
ProcessTrackWriter
- ProcessStackInspector
- ProcessTrackingLine
ProcessParticleCut
- ProcessHadronicElasticModel
- ProcessStackInspector
+ ProcessTrackingLine
CORSIKAprocesses
- CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
- CORSIKA_ADD_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
- target_link_libraries (vertical_EAS
+CORSIKA_ADD_EXAMPLE (cascade_example)
+target_link_libraries (cascade_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
@@ -91,14 +45,14 @@ if (Pythia8_FOUND)
CORSIKAcascade
ProcessCONEXSourceCut
ProcessEnergyLoss
- ProcessObservationPlane
- ProcessInteractionCounter
ProcessTrackWriter
+ ProcessStackInspector
ProcessTrackingLine
+ ProcessProposal
ProcessParticleCut
ProcessOnShellCheck
+ ProcessHadronicElasticModel
ProcessStackInspector
- ProcessLongitudinalProfile
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
@@ -106,21 +60,96 @@ if (Pythia8_FOUND)
CORSIKAenvironment
CORSIKAprocesssequence
)
+
+if (Pythia8_FOUND)
+ CORSIKA_ADD_EXAMPLE (cascade_proton_example)
+ target_link_libraries (cascade_proton_example
+ SuperStupidStack
+ CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ ProcessPythia
+ CORSIKAcascade
+ ProcessEnergyLoss
+ ProcessTrackWriter
+ ProcessStackInspector
+ ProcessTrackingLine
+ ProcessParticleCut
+ ProcessHadronicElasticModel
+ ProcessStackInspector
+ CORSIKAprocesses
+ CORSIKAcascade
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+
+ CORSIKA_ADD_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
+ target_link_libraries (vertical_EAS
+ SuperStupidStack
+ CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ ProcessPythia
+ ProcessUrQMD
+ ProcessSwitch
+ CORSIKAcascade
+ ProcessEnergyLoss
+ ProcessObservationPlane
+ ProcessInteractionCounter
+ ProcessTrackWriter
+ ProcessTrackingLine
+ ProcessParticleCut
+ ProcessStackInspector
+ ProcessLongitudinalProfile
+ CORSIKAprocesses
+ CORSIKAcascade
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
endif()
CORSIKA_ADD_EXAMPLE (stopping_power stopping_power)
target_link_libraries (stopping_power
- SuperStupidStack
- CORSIKAunits
- ProcessEnergyLoss
- CORSIKAparticles
- CORSIKAgeometry
- CORSIKAenvironment
- )
+ SuperStupidStack
+ CORSIKAunits
+ ProcessEnergyLoss
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ )
CORSIKA_ADD_EXAMPLE (staticsequence_example)
target_link_libraries (staticsequence_example
- CORSIKAprocesssequence
- CORSIKAunits
- CORSIKAgeometry
- CORSIKAlogging)
+ CORSIKAprocesssequence
+ CORSIKAunits
+ CORSIKAgeometry
+ CORSIKAlogging)
+
+CORSIKA_ADD_EXAMPLE (proposal_example RUN_OPTIONS 100.)
+target_link_libraries (proposal_example
+ SuperStupidStack
+ CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ CORSIKAcascade
+ ProcessObservationPlane
+ ProcessInteractionCounter
+ ProcessTrackWriter
+ ProcessProposal
+ ProcessTrackingLine
+ ProcessParticleCut
+ ProcessStackInspector
+ ProcessLongitudinalProfile
+ CORSIKAprocesses
+ CORSIKAcascade
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 93902e146..7aa0561c0 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -23,6 +23,8 @@
#include <corsika/geometry/Sphere.h>
+#include <corsika/process/proposal/Interaction.h>
+
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
@@ -135,31 +137,34 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+ random::RNGManager::GetInstance().RegisterRandomStream("proposal");
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay;
// cascade with only HE model ==> HE cut
process::particle_cut::ParticleCut cut(80_GeV);
+ process::proposal::Interaction proposal(env, cut);
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss{showerAxis};
// assemble all processes into an ordered process list
- auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
- << trackWriter;
+ auto sequence = stackInspect << sibyll << sibyllNuc << proposal << decay
+ /* << eLoss */
+ << cut << trackWriter;
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Run();
- eLoss.PrintProfile(); // print longitudinal profile
+ /* eLoss.PrintProfile(); // print longitudinal profile */
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
- cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
- << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+ /* cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl */
+ /* << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl; */
}
diff --git a/Documentation/Examples/proposal_example.cc b/Documentation/Examples/proposal_example.cc
new file mode 100644
index 000000000..3fb41d169
--- /dev/null
+++ b/Documentation/Examples/proposal_example.cc
@@ -0,0 +1,179 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/cascade/Cascade.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/ShowerAxis.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/StackProcess.h>
+#include <corsika/process/interaction_counter/InteractionCounter.h>
+#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
+#include <corsika/process/observation_plane/ObservationPlane.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/process/track_writer/TrackWriter.h>
+#include <corsika/process/proposal/Interaction.h>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+using namespace corsika::setup;
+using namespace corsika::geometry;
+using namespace corsika::environment;
+
+using namespace std;
+using namespace corsika::units::si;
+
+void registerRandomStreams() {
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("proposal");
+ // add PROPOSAL here (?)
+ random::RNGManager::GetInstance().SeedAll();
+}
+
+int main(int argc, char** argv) {
+ if (argc != 2) {
+ std::cerr << "usage: proposal_example <energy/GeV>" << std::endl;
+ return 1;
+ }
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ registerRandomStreams();
+
+ // setup environment, geometry
+ using EnvType = Environment<setup::IEnvironmentModel>;
+ EnvType env;
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ environment::LayeredSphericalAtmosphereBuilder builder{center};
+ builder.setNuclearComposition(
+ {{particles::Code::Nitrogen, particles::Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+ builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+ builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
+
+ builder.assemble(env);
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Electron;
+ auto const mass = particles::GetMass(beamCode);
+ const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
+ << endl;
+
+ auto const observationHeight = 1.4_km + builder.earthRadius;
+ auto const injectionHeight = 112.75_km + builder.earthRadius;
+ auto const t = -observationHeight * cos(thetaRad) +
+ sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
+ si::detail::static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore +
+ Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
+
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, injectionPos, 0_ns});
+
+ std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
+ << std::endl;
+
+ environment::ShowerAxis const showerAxis{injectionPos,
+ (showerCore - injectionPos) * 1.02, env};
+
+ // setup processes, decays and interactions
+
+ // PROPOSAL processs proposal{...};
+ PROPOSAL::InterpolationDef::path_to_tables = "~/.local/share/PROPOSAL/tables/";
+ PROPOSAL::InterpolationDef::path_to_tables_readonly = "~/.local/share/PROPOSAL/tables/";
+
+ process::particle_cut::ParticleCut cut(10_GeV);
+ process::proposal::Interaction proposal(env, cut);
+ process::interaction_counter::InteractionCounter proposalCounted(proposal);
+ process::track_writer::TrackWriter trackWriter("tracks.dat");
+
+ // energy cut; n.b. ParticleCut needs to be modified not to discard EM particles!
+ /* process::particle_cut::ParticleCut cut{60_GeV}; */
+
+ // long. profile; columns for gamma, e+, e- still need to be added
+ process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+ process::observation_plane::ObservationPlane observationLevel(obsPlane,
+ "particles.dat");
+
+ auto sequence = proposalCounted << longprof << proposal << cut << observationLevel << trackWriter;
+ // define air shower object, run simulation
+ tracking_line::TrackingLine tracking;
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+ EAS.Init();
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.SetNodes();
+ // EAS.forceInteraction();
+
+ EAS.Run();
+
+ cut.ShowResults();
+ const HEPEnergyType Efinal =
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+
+ auto const hists = proposalCounted.GetHistogram();
+ hists.saveLab("inthist_lab.txt");
+ hists.saveCMS("inthist_cms.txt");
+
+ longprof.save("longprof.txt");
+
+ std::ofstream finish("finished");
+ finish << "run completed without error" << std::endl;
+}
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 08c5dca48..c405641bb 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -1,7 +1,7 @@
# general
add_subdirectory (NullModel)
# tracking
-add_subdirectory (TrackingLine)
+add_subdirectory (TrackingLine)
# hadron interaction models
add_subdirectory (Sibyll)
add_subdirectory (QGSJetII)
@@ -13,6 +13,8 @@ if (CONEX_FOUND)
endif (CONEX_FOUND)
add_subdirectory (HadronicElasticModel)
add_subdirectory (UrQMD)
+add_subdirectory (SwitchProcess)
+add_subdirectory (Proposal)
# continuous physics
add_subdirectory (EnergyLoss)
@@ -34,6 +36,7 @@ add_subdirectory (SwitchProcess)
add_library (CORSIKAprocesses INTERFACE)
add_dependencies(CORSIKAprocesses ProcessNullModel)
add_dependencies(CORSIKAprocesses ProcessSibyll)
+add_dependencies(CORSIKAprocesses ProcessProposal)
if (Pythia8_FOUND)
add_dependencies(CORSIKAprocesses ProcessPythia)
endif (Pythia8_FOUND)
diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc
index bd5a5fb61..1a154b282 100644
--- a/Processes/ParticleCut/ParticleCut.cc
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -66,12 +66,13 @@ namespace corsika::process {
C8LOG_DEBUG(fmt::format("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV", pid,
energy / 1_GeV,
(fEmEnergy + fInvEnergy + fEnergy) / 1_GeV));
- if (ParticleIsEmParticle(pid)) {
- C8LOG_DEBUG("removing em. particle...");
- fEmEnergy += energy;
- fEmCount += 1;
- return true;
- } else if (ParticleIsInvisible(pid)) {
+ /* if (ParticleIsEmParticle(pid)) { */
+ /* C8LOG_DEBUG("removing em. particle..."); */
+ /* fEmEnergy += energy; */
+ /* fEmCount += 1; */
+ /* return true; */
+ /* } else */
+ if (ParticleIsInvisible(pid)) {
C8LOG_DEBUG("removing inv. particle...");
fInvEnergy += energy;
fInvCount += 1;
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
index 4c0f6371a..af73dbb47 100644
--- a/Processes/ParticleCut/ParticleCut.h
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -31,7 +31,8 @@ namespace corsika::process {
unsigned int fInvCount = 0;
public:
- ParticleCut(const units::si::HEPEnergyType vCut);
+ ParticleCut(const units::si::HEPEnergyType eCut)
+ : fECut(eCut) {}
EProcessReturn DoSecondaries(corsika::setup::StackView&);
@@ -44,6 +45,7 @@ namespace corsika::process {
void ShowResults();
+ units::si::HEPEnergyType GetECut() const { return fECut; }
units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
diff --git a/Processes/Proposal/CMakeLists.txt b/Processes/Proposal/CMakeLists.txt
new file mode 100644
index 000000000..8d2a19788
--- /dev/null
+++ b/Processes/Proposal/CMakeLists.txt
@@ -0,0 +1,46 @@
+# find_package(PROPOSAL REQUIRED ${LIB_INCLUDE})
+
+
+FILE (GLOB MODEL_SOURCES *.cc)
+FILE (GLOB MODEL_HEADERS *.h)
+SET (MODEL_NAMESPACE corsika/process/proposal)
+
+ADD_LIBRARY (ProcessProposal STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessProposal ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+SET_TARGET_PROPERTIES ( ProcessProposal PROPERTIES VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+ # PUBLIC_HEADER "${MODEL_HEADERS}"
+ )
+
+TARGET_LINK_LIBRARIES (
+ ProcessProposal
+ CORSIKAparticles
+ CORSIKAutilities
+ CORSIKAunits
+ CORSIKAthirdparty
+ CORSIKAgeometry
+ CORSIKAenvironment
+ ${PROPOSAL_LIBRARY}
+ )
+
+TARGET_INCLUDE_DIRECTORIES (
+ ProcessProposal
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+TARGET_INCLUDE_DIRECTORIES (
+ ProcessProposal
+ SYSTEM
+ PUBLIC ${PROPOSAL_INCLUDE_DIR}
+ INTERFACE ${PROPOSAL_INCLUDE_DIR}
+ )
+
+INSTALL (
+ TARGETS ProcessProposal
+ LIBRARY DESTINATION lib
+ ARCHIVE DESTINATION lib
+# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+ )
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
new file mode 100644
index 000000000..5914d5296
--- /dev/null
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -0,0 +1,106 @@
+#include <PROPOSAL/PROPOSAL.h>
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/proposal/ContinuousProcess.h>
+#include <corsika/process/proposal/Interaction.h>
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+
+using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupTrack = corsika::setup::Trajectory;
+
+namespace corsika::process::proposal {
+ using namespace corsika::setup;
+ using namespace corsika::environment;
+ using namespace corsika::units::si;
+
+ unordered_map<particles::Code, PROPOSAL::ParticleDef> ContinuousProcess::particles{
+ {particles::Code::Gamma, PROPOSAL::GammaDef()},
+ {particles::Code::Electron, PROPOSAL::EMinusDef()},
+ {particles::Code::Positron, PROPOSAL::EPlusDef()},
+ {particles::Code::MuMinus, PROPOSAL::MuMinusDef()},
+ {particles::Code::MuPlus, PROPOSAL::MuPlusDef()},
+ {particles::Code::TauPlus, PROPOSAL::TauPlusDef()},
+ {particles::Code::TauMinus, PROPOSAL::TauMinusDef()},
+ };
+
+ bool ContinuousProcess::CanInteract(particles::Code pcode) const noexcept {
+ auto search = particles.find(pcode);
+ if (search != particles.end()) return true;
+ return false;
+ }
+
+ template <>
+ ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env,
+ CORSIKA_ParticleCut const& _cut)
+ : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
+ false)) {
+ auto all_compositions = std::vector<const NuclearComposition*>();
+ _env.GetUniverse()->walk([&](auto& vtn) {
+ if (vtn.HasModelProperties())
+ all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
+ });
+ for (auto& ncarg : all_compositions) {
+ auto comp_vec = std::vector<PROPOSAL::Components::Component>();
+ auto frac_iter = ncarg->GetFractions().cbegin();
+ for (auto& pcode : ncarg->GetComponents()) {
+ comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
+ *frac_iter);
+ ++frac_iter;
+ }
+ media[ncarg] = PROPOSAL::Medium(
+ "Modified Air", 1., PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
+ PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
+ PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
+ }
+ }
+
+ HEPEnergyType ContinuousProcess::TotalEnergyLoss(SetupParticle const& vP,
+ GrammageType const vDX) {
+ auto calc_ptr = GetCalculator(vP);
+ auto upper_energy = calc_ptr->second->UpperLimitTrackIntegral(
+ vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm);
+ std::cout << "upper_energy: " << upper_energy << "MeV" << std::endl;
+ return upper_energy * 1_MeV;
+ }
+
+ void ContinuousProcess::Init() {}
+ template <>
+ EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP, SetupTrack const& vT) {
+ if (vP.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
+ std::cout << "DoContinuous..." << std::endl;
+ GrammageType const dX =
+ vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
+ HEPEnergyType dE = TotalEnergyLoss(vP, dX);
+ auto E = vP.GetEnergy();
+ const auto Ekin = E - vP.GetMass();
+ auto Enew = E + dE;
+ auto status = process::EProcessReturn::eOk;
+ if (-dE > Ekin) {
+ dE = -Ekin;
+ Enew = vP.GetMass();
+ status = process::EProcessReturn::eParticleAbsorbed;
+ }
+ vP.SetEnergy(Enew);
+ auto pnew = vP.GetMomentum();
+ vP.SetMomentum(pnew * Enew / pnew.GetNorm());
+ return status;
+ }
+
+ template <>
+ units::si::LengthType ContinuousProcess::MaxStepLength(SetupParticle const& vP,
+ SetupTrack const& vT) {
+ auto constexpr dX = 1_g / square(1_cm);
+ auto const dE = -TotalEnergyLoss(vP, dX); // dE > 0
+ auto const maxLoss = 0.01 * vP.GetEnergy();
+ auto const maxGrammage = maxLoss / dE * dX;
+
+ return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT,
+ maxGrammage) *
+ 1.0001; // to make sure particle gets absorbed when DoContinuous() is called
+ }
+
+} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
new file mode 100644
index 000000000..58b19ec29
--- /dev/null
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -0,0 +1,128 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _corsika_process_proposal_interaction_h_
+#define _corsika_process_proposal_interaction_h_
+
+#include <PROPOSAL/PROPOSAL.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <unordered_map>
+#include "PROPOSAL/PROPOSAL.h"
+
+using std::unordered_map;
+
+using namespace corsika::environment;
+
+using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
+
+namespace corsika::process::proposal {
+
+ class ContinuousProcess
+ : public corsika::process::ContinuousProcess<ContinuousProcess> {
+ private:
+ shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
+
+ static unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
+ unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
+
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("proposal");
+
+ bool CanInteract(particles::Code pcode) const noexcept;
+
+ struct interaction_hash {
+ size_t operator()(const std::pair<const NuclearComposition*, particles::Code>& p) const
+ {
+ auto hash1 = std::hash<const NuclearComposition*>{}(p.first);
+ auto hash2 = std::hash<particles::Code>{}(p.second);
+ return hash1 ^ hash2;
+ }
+ };
+
+ unordered_map<std::pair<const NuclearComposition*, particles::Code>, unique_ptr<PROPOSAL::Displacement>, interaction_hash> calc;
+
+ auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
+ auto medium = media.at(&comp);
+ if (corsika_code == particles::Code::Gamma) {
+ auto cross = GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true);
+ auto [insert_it, success] =
+ calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
+ if (corsika_code == particles::Code::Electron) {
+ auto cross = GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true);
+ auto [insert_it, success] =
+ calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
+ if (corsika_code == particles::Code::Positron) {
+ auto cross = GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true);
+ auto [insert_it, success] =
+ calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
+ if (corsika_code == particles::Code::MuMinus) {
+ auto cross = GetStdCrossSections(PROPOSAL::MuMinusDef(), media.at(&comp), cut, true);
+ auto [insert_it, success] =
+ calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
+ if (corsika_code == particles::Code::MuPlus) {
+ auto cross = GetStdCrossSections(PROPOSAL::MuPlusDef(), media.at(&comp), cut, true);
+ auto [insert_it, success] =
+ calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
+ if (corsika_code == particles::Code::TauMinus) {
+ auto cross = GetStdCrossSections(PROPOSAL::TauMinusDef(), media.at(&comp), cut, true);
+ auto [insert_it, success] =
+ calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
+ if (corsika_code == particles::Code::TauPlus) {
+ auto cross = GetStdCrossSections(PROPOSAL::TauPlusDef(), media.at(&comp), cut, true);
+ auto [insert_it, success] =
+ calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
+ throw std::runtime_error("PROPOSAL could not find corresponding builder");
+ }
+
+ template <typename Particle>
+ auto GetCalculator(Particle& vP) {
+ auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+ auto calc_it = calc.find(std::make_pair(&comp, vP.GetPID()));
+ if (calc_it != calc.end()) return calc_it;
+ return BuildCalculator(vP.GetPID(), comp);
+ }
+
+ units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
+ const units::si::GrammageType);
+
+ public:
+ template <typename TEnvironment>
+ ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+
+ void Init();
+
+ template <typename Particle, typename Track>
+ EProcessReturn DoContinuous(Particle&, Track const&) ;
+
+ template <typename Particle, typename Track>
+ units::si::LengthType MaxStepLength(Particle const& p, Track const& track) ;
+ };
+} // namespace corsika::process::proposal
+
+#endif
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
new file mode 100644
index 000000000..4394ae868
--- /dev/null
+++ b/Processes/Proposal/Interaction.cc
@@ -0,0 +1,113 @@
+
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/proposal/Interaction.h>
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+#include <limits>
+#include <memory>
+#include <random>
+#include <tuple>
+
+using Component_PROPOSAL = PROPOSAL::Components::Component;
+
+namespace corsika::process::proposal {
+ using namespace corsika::setup;
+ using namespace corsika::environment;
+ using namespace corsika::units::si;
+
+ std::unordered_map<particles::Code, PROPOSAL::ParticleDef> Interaction::particles{
+ {particles::Code::Gamma, PROPOSAL::GammaDef()},
+ {particles::Code::Electron, PROPOSAL::EMinusDef()},
+ {particles::Code::Positron, PROPOSAL::EPlusDef()},
+ {particles::Code::MuMinus, PROPOSAL::MuMinusDef()},
+ {particles::Code::MuPlus, PROPOSAL::MuPlusDef()},
+ {particles::Code::TauPlus, PROPOSAL::TauPlusDef()},
+ {particles::Code::TauMinus, PROPOSAL::TauMinusDef()},
+ };
+
+ bool Interaction::CanInteract(particles::Code pcode) const noexcept {
+ if (particles.find(pcode) != particles.end()) return true;
+ return false;
+ }
+
+ template <>
+ Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut const& _cut)
+ : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
+ false))
+ , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
+ auto all_compositions = std::vector<const NuclearComposition*>();
+ _env.GetUniverse()->walk([&](auto& vtn) {
+ if (vtn.HasModelProperties()) {
+ all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
+ }
+ });
+ for (auto& ncarg : all_compositions) {
+ auto comp_vec = std::vector<Component_PROPOSAL>();
+ auto frac_iter = ncarg->GetFractions().cbegin();
+ for (auto& pcode : ncarg->GetComponents()) {
+ comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
+ *frac_iter);
+ ++frac_iter;
+ }
+ media[ncarg] = PROPOSAL::Medium(
+ "Modified Air", 1., PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
+ PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
+ PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
+ }
+ }
+
+ template <>
+ corsika::process::EProcessReturn Interaction::DoInteraction(
+ setup::StackView::StackIterator& vP) {
+ if (CanInteract(vP.GetPID())) {
+ auto calc = GetCalculator(vP); // [CrossSections]
+ std::uniform_real_distribution<double> distr(0., 1.);
+ auto [type, comp_ptr, v] =
+ get<INTERACTION>(calc->second)
+ ->TypeInteraction(vP.GetEnergy() / 1_MeV, distr(fRNG));
+ auto rnd =
+ vector<double>(get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type));
+ for (auto& it : rnd) it = distr(fRNG);
+ auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
+ vP.GetPosition().GetY() / 1_cm,
+ vP.GetPosition().GetZ() / 1_cm);
+ auto d = vP.GetDirection().GetComponents();
+ auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
+ d.GetZ().magnitude());
+ auto loss = make_tuple(static_cast<int>(type), point, direction,
+ v * vP.GetEnergy() / 1_MeV, 0.);
+ auto sec = get<SECONDARIES>(calc->second)
+ ->CalculateSecondaries(vP.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd);
+ for (auto& s : sec) {
+ auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
+ auto vec = corsika::geometry::QuantityVector(
+ get<PROPOSAL::Loss::DIRECTION>(s).GetX() * E,
+ get<PROPOSAL::Loss::DIRECTION>(s).GetY() * E,
+ get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * E);
+ auto p = corsika::stack::MomentumVector(
+ corsika::geometry::RootCoordinateSystem::GetInstance()
+ .GetRootCoordinateSystem(),
+ vec);
+ auto sec_code = corsika::particles::ConvertFromPDG(
+ static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
+ vP.AddSecondary(make_tuple(sec_code, E, p, vP.GetPosition(), vP.GetTime()));
+ }
+ }
+ return process::EProcessReturn::eOk;
+ }
+
+ template <>
+ corsika::units::si::GrammageType Interaction::GetInteractionLength(
+ setup::Stack::StackIterator const& vP) {
+ if (CanInteract(vP.GetPID())) {
+ auto calc = GetCalculator(vP);
+ return get<INTERACTION>(calc->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
+ (1_cm * 1_cm);
+ }
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
new file mode 100644
index 000000000..54597a656
--- /dev/null
+++ b/Processes/Proposal/Interaction.h
@@ -0,0 +1,107 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _corsika_process_proposal_interaction_h_
+#define _corsika_process_proposalythia_interaction_h_
+
+#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <unordered_map>
+#include "PROPOSAL/PROPOSAL.h"
+
+using namespace corsika::environment;
+
+using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
+using std::make_pair;
+using std::make_tuple;
+
+namespace corsika::process::proposal {
+
+ class Interaction : public corsika::process::InteractionProcess<Interaction> {
+
+ shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
+
+ static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
+ std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
+
+ corsika::random::RNG& fRNG;
+
+ bool CanInteract(particles::Code pcode) const noexcept;
+
+ using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
+ unique_ptr<PROPOSAL::Interaction>>;
+ using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
+
+ struct interaction_hash {
+ size_t operator()(const calc_key_t& p) const {
+ return std::hash<const NuclearComposition*>{}(p.first) ^
+ std::hash<particles::Code>{}(p.second);
+ }
+ };
+
+ std::unordered_map<calc_key_t, calculator_t, interaction_hash> calculators;
+
+ template <typename Particle>
+ auto BuildCalculator(particles::Code corsika_code, Particle p_def,
+ NuclearComposition const& comp) {
+ auto cross = GetStdCrossSections(p_def, media.at(&comp), cut, true);
+ auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
+ auto [insert_it, success] = calculators.insert(
+ {make_pair(&comp, corsika_code),
+ make_tuple(PROPOSAL::make_secondaries(inter_types, p_def, media.at(&comp)),
+ PROPOSAL::make_interaction(cross, true))});
+ return insert_it;
+ }
+
+ auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
+ if (corsika_code == particles::Code::Gamma)
+ return BuildCalculator(particles::Code::Gamma, PROPOSAL::GammaDef(), comp);
+ if (corsika_code == particles::Code::Electron)
+ return BuildCalculator(particles::Code::Electron, PROPOSAL::EMinusDef(), comp);
+ if (corsika_code == particles::Code::Positron)
+ return BuildCalculator(particles::Code::Positron, PROPOSAL::EPlusDef(), comp);
+ if (corsika_code == particles::Code::MuMinus)
+ return BuildCalculator(particles::Code::MuMinus, PROPOSAL::MuMinusDef(), comp);
+ if (corsika_code == particles::Code::MuPlus)
+ return BuildCalculator(particles::Code::MuPlus, PROPOSAL::MuPlusDef(), comp);
+ if (corsika_code == particles::Code::TauMinus)
+ return BuildCalculator(particles::Code::TauMinus, PROPOSAL::TauMinusDef(), comp);
+ if (corsika_code == particles::Code::TauPlus)
+ return BuildCalculator(particles::Code::TauPlus, PROPOSAL::TauPlusDef(), comp);
+ throw std::runtime_error("PROPOSAL could not find corresponding builder");
+ }
+
+ enum { SECONDARIES, INTERACTION };
+ template <typename Particle>
+ auto GetCalculator(Particle& vP) {
+ auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+ auto calc_it = calculators.find(make_pair(&comp, vP.GetPID()));
+ if (calc_it != calculators.end()) return calc_it;
+ return BuildCalculator(vP.GetPID(), comp);
+ }
+
+ public:
+ template <typename TEnvironment>
+ Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+
+ void Init() {};
+
+ template <typename Particle>
+ corsika::process::EProcessReturn DoInteraction(Particle&);
+
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle const& p);
+ }; // namespace corsika::process::proposal
+} // namespace corsika::process::proposal
+#endif
--
GitLab