From 68fda5ba5060863be2d2c1d649a0255be91ee441 Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@tu-dortmund.de>
Date: Thu, 9 Jul 2020 11:12:02 +0000
Subject: [PATCH] First electromagnetic shower cascade is completly possible.
There still have to be added multiple scattering effects and some deflections effects.
The directions and asymetrie of interaction products will be sampled.
There seems to be an issue with the particle cut or a implementation failure on proposal side that if a particle energy is below the cut energy, the particle is simply deleted without taken the particle cut relative difference into account.
---
Documentation/Examples/proposal_example.cc | 9 +-
Processes/Proposal/ContinuousProcess.cc | 62 +++++-------
Processes/Proposal/ContinuousProcess.h | 106 +++++++++------------
Processes/Proposal/Interaction.cc | 5 +-
Processes/Proposal/Interaction.h | 23 ++---
5 files changed, 90 insertions(+), 115 deletions(-)
diff --git a/Documentation/Examples/proposal_example.cc b/Documentation/Examples/proposal_example.cc
index 3fb41d169..de16fc65d 100644
--- a/Documentation/Examples/proposal_example.cc
+++ b/Documentation/Examples/proposal_example.cc
@@ -29,6 +29,7 @@
#include <corsika/utl/CorsikaFenv.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/proposal/Interaction.h>
+#include <corsika/process/proposal/ContinuousProcess.h>
#include <iomanip>
#include <iostream>
@@ -137,6 +138,7 @@ int main(int argc, char** argv) {
process::particle_cut::ParticleCut cut(10_GeV);
process::proposal::Interaction proposal(env, cut);
+ process::proposal::ContinuousProcess em_continuous(env, cut);
process::interaction_counter::InteractionCounter proposalCounted(proposal);
process::track_writer::TrackWriter trackWriter("tracks.dat");
@@ -150,7 +152,12 @@ int main(int argc, char** argv) {
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
- auto sequence = proposalCounted << longprof << proposal << cut << observationLevel << trackWriter;
+ auto sequence = proposalCounted
+ << cut
+ << em_continuous
+ << longprof
+ << observationLevel
+ << trackWriter;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 5914d5296..d222a80b1 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -1,6 +1,7 @@
#include <PROPOSAL/PROPOSAL.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
#include <corsika/setup/SetupEnvironment.h>
@@ -28,20 +29,20 @@ namespace corsika::process::proposal {
};
bool ContinuousProcess::CanInteract(particles::Code pcode) const noexcept {
- auto search = particles.find(pcode);
- if (search != particles.end()) return true;
+ if (particles.find(pcode) != 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)) {
+ CORSIKA_ParticleCut& _cut)
+ : cut(_cut)
+ , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
auto all_compositions = std::vector<const NuclearComposition*>();
_env.GetUniverse()->walk([&](auto& vtn) {
- if (vtn.HasModelProperties())
+ if (vtn.HasModelProperties()) {
all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
+ }
});
for (auto& ncarg : all_compositions) {
auto comp_vec = std::vector<PROPOSAL::Components::Component>();
@@ -60,47 +61,32 @@ namespace corsika::process::proposal {
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;
+ std::cout << "distance: " << vDX / 1_g * 1_cm * 1_cm << std::endl;
+ return vP.GetEnergy() - GetCalculator(vP)->second->UpperLimitTrackIntegral(
+ vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm) *
+ 1_MeV;
}
- void ContinuousProcess::Init() {}
template <>
- EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP, SetupTrack const& vT) {
+ 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;
+ auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
+ vP.SetEnergy(vP.GetEnergy() - TotalEnergyLoss(vP, dX));
+ if (vP.GetEnergy() < cut.GetECut()) return process::EProcessReturn::eParticleAbsorbed;
+ vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
+ return process::EProcessReturn::eOk;
}
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
+ auto energy_lim = 0.9 * vP.GetEnergy() / 1_MeV;
+ auto grammage = GetCalculator(vP)->second->SolveTrackIntegral(vP.GetEnergy() / 1_MeV,
+ energy_lim) *
+ 1_g / square(1_cm);
+ return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage) *
+ 1.0001;
}
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index 58b19ec29..e0cfe6e26 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -24,6 +24,7 @@
using std::unordered_map;
using namespace corsika::environment;
+using namespace corsika::units::si;
using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
@@ -31,75 +32,56 @@ namespace corsika::process::proposal {
class ContinuousProcess
: public corsika::process::ContinuousProcess<ContinuousProcess> {
- private:
- shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
-
+ CORSIKA_ParticleCut& cut;
+ corsika::random::RNG& fRNG;
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;
- }
+ using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
+
+ struct disp_hash {
+ size_t operator()(const calc_key_t& p) const {
+ return std::hash<const NuclearComposition*>{}(p.first) ^
+ std::hash<particles::Code>{}(p.second);
+ }
};
- unordered_map<std::pair<const NuclearComposition*, particles::Code>, unique_ptr<PROPOSAL::Displacement>, interaction_hash> calc;
+ unordered_map<calc_key_t, unique_ptr<PROPOSAL::Displacement>, disp_hash> calc;
+
+ template <typename Particle>
+ auto BuildCalculator(particles::Code code, Particle p_def,
+ NuclearComposition const& comp) {
+ auto medium = media.at(&comp);
+ auto cross = GetStdCrossSections(
+ p_def, media.at(&comp),
+ make_shared<const PROPOSAL::EnergyCutSettings>(cut.GetECut() / 1_MeV, 1, false),
+ true);
+ auto [insert_it, success] = calc.insert(
+ {std::make_pair(&comp, code), PROPOSAL::make_displacement(cross, true)});
+ return insert_it;
+ }
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");
+ 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");
}
-
+
template <typename Particle>
auto GetCalculator(Particle& vP) {
auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
@@ -113,15 +95,15 @@ namespace corsika::process::proposal {
public:
template <typename TEnvironment>
- ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+ ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut& cut);
- void Init();
+ void Init(){};
template <typename Particle, typename Track>
- EProcessReturn DoContinuous(Particle&, Track const&) ;
+ EProcessReturn DoContinuous(Particle&, Track const&);
template <typename Particle, typename Track>
- units::si::LengthType MaxStepLength(Particle const& p, Track const& track) ;
+ units::si::LengthType MaxStepLength(Particle const& p, Track const& track);
};
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index 4394ae868..60ae3c81f 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -35,9 +35,8 @@ namespace corsika::process::proposal {
}
template <>
- Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut const& _cut)
- : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
- false))
+ Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut& _cut)
+ : cut(_cut)
, fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
auto all_compositions = std::vector<const NuclearComposition*>();
_env.GetUniverse()->walk([&](auto& vtn) {
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index 54597a656..6d0ad2f37 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -21,6 +21,7 @@
#include "PROPOSAL/PROPOSAL.h"
using namespace corsika::environment;
+using namespace corsika::units::si;
using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
using std::make_pair;
@@ -29,14 +30,11 @@ using std::make_tuple;
namespace corsika::process::proposal {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
-
- shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
-
+ CORSIKA_ParticleCut& cut;
+ corsika::random::RNG& fRNG;
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>,
@@ -53,12 +51,15 @@ namespace corsika::process::proposal {
std::unordered_map<calc_key_t, calculator_t, interaction_hash> calculators;
template <typename Particle>
- auto BuildCalculator(particles::Code corsika_code, Particle p_def,
+ auto BuildCalculator(particles::Code code, Particle p_def,
NuclearComposition const& comp) {
- auto cross = GetStdCrossSections(p_def, media.at(&comp), cut, true);
+ auto cross = GetStdCrossSections(
+ p_def, media.at(&comp),
+ make_shared<const PROPOSAL::EnergyCutSettings>(cut.GetECut() / 1_MeV, 1, false),
+ true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
- {make_pair(&comp, corsika_code),
+ {make_pair(&comp, code),
make_tuple(PROPOSAL::make_secondaries(inter_types, p_def, media.at(&comp)),
PROPOSAL::make_interaction(cross, true))});
return insert_it;
@@ -93,15 +94,15 @@ namespace corsika::process::proposal {
public:
template <typename TEnvironment>
- Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+ Interaction(TEnvironment const& env, CORSIKA_ParticleCut& cut);
- void Init() {};
+ 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