From 367eb78b666f5cb15ca7e5d1e5b32d081f7d8695 Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@tu-dortmund.de>
Date: Thu, 1 Oct 2020 08:42:05 +0000
Subject: [PATCH] Now stores interpolation tables in $CORSIKA_DATA/PROPOSAL.
Remove proposal from cascade example, improve error message. Improve output,
and shower diagnostics: hunting all the energy Produce REAL error message.
---
Documentation/Examples/cascade_example.cc | 17 ++++-------
Documentation/Examples/em_shower.cc | 14 +++++----
Documentation/Examples/vertical_EAS.cc | 14 +++++----
.../ObservationPlane/ObservationPlane.cc | 30 ++++++++++++++++---
Processes/ObservationPlane/ObservationPlane.h | 7 +++++
Processes/ParticleCut/ParticleCut.cc | 9 ++++++
Processes/ParticleCut/ParticleCut.h | 3 +-
Processes/Proposal/ContinuousProcess.cc | 19 ++++++++++--
Processes/Proposal/ContinuousProcess.h | 8 +++--
Processes/Proposal/ProposalProcessBase.cc | 13 ++++++++
10 files changed, 101 insertions(+), 33 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index c988b4f5d..ece23e9f5 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -23,8 +23,6 @@
#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>
@@ -137,35 +135,32 @@ 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, true, true);
- // 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 /* << proposal*/
- << decay
- /* << eLoss */
- << cut << trackWriter;
+ auto sequence = stackInspect << sibyll << sibyllNuc << 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;
+ cut.Reset();
}
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index d83cdacf5..e84931014 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -130,9 +130,6 @@ int main(int argc, char** argv) {
// 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, false, true);
process::proposal::Interaction proposal(env, cut);
process::proposal::ContinuousProcess em_continuous(env, cut);
@@ -147,7 +144,7 @@ int main(int argc, char** argv) {
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
- auto sequence = proposalCounted << cut << em_continuous << longprof << observationLevel
+ auto sequence = proposalCounted << em_continuous << longprof << cut << observationLevel
<< trackWriter;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
@@ -160,11 +157,16 @@ int main(int argc, char** argv) {
EAS.Run();
cut.ShowResults();
+ em_continuous.ShowResults();
+ observationLevel.ShowResults();
const HEPEnergyType Efinal =
- cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-
+ observationLevel.Reset();
+ cut.Reset();
+ em_continuous.Reset();
+
auto const hists = proposalCounted.GetHistogram();
hists.saveLab("inthist_lab.txt");
hists.saveCMS("inthist_cms.txt");
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 9c0a72a87..5f7f7d343 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -184,9 +184,6 @@ int main(int argc, char** argv) {
decaySibyll.PrintDecayConfig();
// 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{60_GeV, false, true};
process::proposal::Interaction proposal(env, cut);
process::proposal::ContinuousProcess em_continuous(env, cut);
@@ -210,8 +207,8 @@ int main(int argc, char** argv) {
55_GeV);
auto decaySequence = decayPythia << decaySibyll;
- auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << cut << em_continuous
- << longprof << observationLevel;
+ auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << em_continuous
+ << cut << longprof << observationLevel;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
@@ -225,10 +222,15 @@ int main(int argc, char** argv) {
cut.ShowResults();
+ em_continuous.ShowResults();
+ observationLevel.ShowResults();
const HEPEnergyType Efinal =
- cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ observationLevel.Reset();
+ cut.Reset();
+ em_continuous.Reset();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index 2f58e96aa..ce2a20b7d 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -9,6 +9,7 @@
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <fstream>
+#include <iostream>
using namespace corsika::process::observation_plane;
using namespace corsika::units::si;
@@ -17,7 +18,9 @@ ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
std::string const& filename, bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
- , deleteOnHit_(deleteOnHit) {
+ , deleteOnHit_(deleteOnHit)
+ , energy_ground_(0_GeV)
+ , count_ground_(0) {
outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
}
@@ -33,13 +36,19 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
return process::EProcessReturn::eOk;
}
+ const auto energy = particle.GetEnergy();
outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
- << particle.GetEnergy() * (1 / 1_eV) << ' '
+ << energy / 1_eV << ' '
<< (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
<< std::endl;
- if (deleteOnHit_) { return process::EProcessReturn::eParticleAbsorbed; }
- return process::EProcessReturn::eOk;
+ if (deleteOnHit_) {
+ count_ground_++;
+ energy_ground_ += energy;
+ return process::EProcessReturn::eParticleAbsorbed;
+ } else {
+ return process::EProcessReturn::eOk;
+ }
}
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
@@ -55,3 +64,16 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
}
+
+void ObservationPlane::ShowResults() const {
+ std::cout << " ******************************" << std::endl
+ << " ObservationPlane: " << std::endl;
+ std::cout << " energy in ground (GeV) : " << energy_ground_ / 1_GeV << std::endl
+ << " no. of particles in ground : " << count_ground_ << std::endl
+ << " ******************************" << std::endl;
+}
+
+void ObservationPlane::Reset() {
+ energy_ground_ = 0_GeV;
+ count_ground_ = 0;
+}
diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h
index 9b1ecf613..96d2e77d4 100644
--- a/Processes/ObservationPlane/ObservationPlane.h
+++ b/Processes/ObservationPlane/ObservationPlane.h
@@ -36,9 +36,16 @@ namespace corsika::process::observation_plane {
corsika::setup::Stack::ParticleType const&,
corsika::setup::Trajectory const& vTrajectory);
+ void ShowResults() const;
+ void Reset();
+ corsika::units::si::HEPEnergyType GetEnergyGround() const { return energy_ground_; }
+
private:
geometry::Plane const plane_;
std::ofstream outputStream_;
bool const deleteOnHit_;
+
+ units::si::HEPEnergyType energy_ground_ = 0 * units::si::electronvolt;
+ unsigned int count_ground_ = 0;
};
} // namespace corsika::process::observation_plane
diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc
index 5fdfbc085..3fa22d9f7 100644
--- a/Processes/ParticleCut/ParticleCut.cc
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -131,5 +131,14 @@ namespace corsika::process {
" ******************************",
fEmEnergy / 1_GeV, uiEmCount, fInvEnergy / 1_GeV, uiInvCount, fEnergy / 1_GeV));
}
+
+ void ParticleCut::Reset() {
+ fEmEnergy = 0_GeV;
+ uiEmCount = 0;
+ fInvEnergy = 0_GeV;
+ uiInvCount = 0;
+ fEnergy = 0_GeV;
+ }
+
} // namespace particle_cut
} // namespace corsika::process
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
index a7f3898b1..21a6181f0 100644
--- a/Processes/ParticleCut/ParticleCut.h
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -52,7 +52,8 @@ namespace corsika::process {
bool ParticleIsEmParticle(particles::Code) const;
- void ShowResults();
+ void ShowResults() const;
+ void Reset();
units::si::HEPEnergyType GetECut() const { return fECut; }
units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 61a64c0bb..ebd0ba087 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -17,6 +17,8 @@
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
+#include <iostream>
+
namespace corsika::process::proposal {
using namespace corsika::units::si;
@@ -97,10 +99,11 @@ namespace corsika::process::proposal {
auto final_energy = get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
vP.GetEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
1_MeV;
-
+ auto dE = vP.GetEnergy() - final_energy;
+ energy_lost_ += dE;
+
// if the particle has a charge take multiple scattering into account
- if (vP.GetChargeNumber() != 0)
- Scatter(vP, vP.GetEnergy() - final_energy, dX);
+ if (vP.GetChargeNumber() != 0) Scatter(vP, dE, dX);
// Update the energy and absorbe the particle if it's below the energy
// threshold, because it will no longer propagated.
@@ -131,4 +134,14 @@ namespace corsika::process::proposal {
return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage);
}
+ void ContinuousProcess::ShowResults() const {
+ std::cout << " ******************************" << std::endl
+ << " PROCESS::ContinuousProcess: " << std::endl;
+ std::cout << " energy lost dE (GeV) : " << energy_lost_ / 1_GeV << std::endl;
+ }
+
+ void ContinuousProcess::Reset() {
+ energy_lost_ = 0_GeV;
+ }
+
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index 9449d9c74..6a3dd177d 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -20,8 +20,6 @@
namespace corsika::process::proposal {
- using namespace corsika::units::si;
-
//!
//! Electro-magnetic and gamma continous losses produced by proposal. It makes
//! use of interpolation tables which are runtime intensive calculation, but can be
@@ -37,6 +35,8 @@ namespace corsika::process::proposal {
std::unordered_map<calc_key_t, calc_t, hash>
calc; //!< Stores the displacement and scattering calculators.
+ units::si::HEPEnergyType energy_lost_ = 0 * units::si::electronvolt;
+
//!
//! Build the displacement and scattering calculators and add it to calc.
//!
@@ -74,5 +74,9 @@ namespace corsika::process::proposal {
//!
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle const&, Track const&);
+
+ void ShowResults() const;
+ void Reset();
+ corsika::units::si::HEPEnergyType GetEnergyLost() const { return energy_lost_; }
};
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/ProposalProcessBase.cc b/Processes/Proposal/ProposalProcessBase.cc
index 815f032b5..b24143570 100644
--- a/Processes/Proposal/ProposalProcessBase.cc
+++ b/Processes/Proposal/ProposalProcessBase.cc
@@ -14,6 +14,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
+#include <cstdlib>
#include <limits>
#include <memory>
#include <random>
@@ -50,6 +51,18 @@ namespace corsika::process::proposal {
PROPOSAL::InterpolationDef::order_of_interpolation = 2;
PROPOSAL::InterpolationDef::nodes_cross_section = 100;
PROPOSAL::InterpolationDef::nodes_propagate = 1000;
+
+ //! If corsika data exist store interpolation tables to the corresponding
+ //! path, otherwise interpolation tables would only stored in main memory if
+ //! no explicit intrpolation def is specified.
+ if (auto data_path = std::getenv("CORSIKA_DATA")) {
+ PROPOSAL::InterpolationDef::path_to_tables = std::string(data_path) + "/PROPOSAL";
+ } else {
+ throw std::runtime_error(
+ "It is not recommended to run PROPOSAL without its tables in "
+ "$CORSIKA_DATA/PROPOSAL. This would be extremely slow. Please provide the "
+ "table directory. ");
+ }
}
size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept {
--
GitLab