From 40f2dc51d68a82f117822d1cacf0f3c1ce6032b4 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Fri, 2 Oct 2020 10:11:26 +0200
Subject: [PATCH] improve output, and shower diagnostics: hunting all the
energy
---
Documentation/Examples/cascade_example.cc | 1 +
Documentation/Examples/em_shower.cc | 11 +++++---
Documentation/Examples/vertical_EAS.cc | 11 +++++---
.../ObservationPlane/ObservationPlane.cc | 25 ++++++++++++++++---
Processes/ObservationPlane/ObservationPlane.h | 7 ++++++
Processes/ParticleCut/ParticleCut.cc | 25 +++++++++++++------
Processes/ParticleCut/ParticleCut.h | 3 ++-
Processes/Proposal/ContinuousProcess.cc | 19 +++++++++++---
Processes/Proposal/ContinuousProcess.h | 8 ++++--
9 files changed, 88 insertions(+), 22 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 28c0315f7..ece23e9f5 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -162,4 +162,5 @@ int main() {
<< "relative difference (%): " << (Efinal / E0 - 1) * 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 5bc85b48c..e84931014 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -144,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;
@@ -157,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 be0749fed..5f7f7d343 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -207,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;
@@ -222,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 9907d3ef8..be331f223 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;
}
@@ -32,13 +35,16 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
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_) {
+ count_ground_++;
+ energy_ground_ += energy;
return process::EProcessReturn::eParticleAbsorbed;
} else {
return process::EProcessReturn::eOk;
@@ -58,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 e59addccc..ab3cd652a 100644
--- a/Processes/ParticleCut/ParticleCut.cc
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -94,15 +94,26 @@ namespace corsika::process {
energy_ = 0_GeV;
}
- void ParticleCut::ShowResults() {
+ void ParticleCut::ShowResults() const {
cout << " ******************************" << endl
- << " ParticleCut: " << endl
- << " energy in em. component (GeV): " << emEnergy_ / 1_GeV << endl
- << " no. of em. particles injected: " << emCount_ << endl
- << " energy in inv. component (GeV): " << invEnergy_ / 1_GeV << endl
- << " no. of inv. particles injected: " << invCount_ << endl
- << " energy below particle cut (GeV): " << energy_ / 1_GeV << endl
+ << " ParticleCut: " << endl;
+ if (discardEm_)
+ cout << " energy in em. component (GeV) : " << emEnergy_ / 1_GeV << endl
+ << " no. of em. particles removed : " << emCount_ << endl;
+ if (discardInv_)
+ cout << " energy in inv. component (GeV) : " << invEnergy_ / 1_GeV << endl
+ << " no. of inv. particles removed : " << invCount_ << endl;
+ cout << " energy below energy threshold (GeV): " << energy_ / 1_GeV << endl
<< " ******************************" << endl;
}
+
+ void ParticleCut::Reset() {
+ emEnergy_ = 0_GeV;
+ emCount_ = 0;
+ invEnergy_ = 0_GeV;
+ invCount_ = 0;
+ energy_ = 0_GeV;
+ }
+
} // namespace particle_cut
} // namespace corsika::process
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
index 0edab21e3..9d17c49dc 100644
--- a/Processes/ParticleCut/ParticleCut.h
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -38,7 +38,8 @@ namespace corsika::process {
bool ParticleIsEmParticle(particles::Code) const;
- void ShowResults();
+ void ShowResults() const;
+ void Reset();
units::si::HEPEnergyType GetECut() const { return eCut_; }
units::si::HEPEnergyType GetInvEnergy() const { return invEnergy_; }
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
--
GitLab