From b6d17c302ea04efcc1e66e8807d5332a345db285 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 6 Oct 2020 12:36:14 +0200
Subject: [PATCH] final adaptions of proposal and history
---
Documentation/Examples/em_shower.cc | 9 ++--
Documentation/Examples/vertical_EAS.cc | 3 +-
Framework/Cascade/Cascade.h | 2 -
Framework/Cascade/testCascade.h | 1 +
.../StackInterface/StackIteratorInterface.h | 2 -
Processes/ParticleCut/ParticleCut.h | 3 --
Processes/ParticleCut/testParticleCut.cc | 8 ++--
Processes/Proposal/Interaction.cc | 43 ++++++++++---------
Processes/Proposal/Interaction.h | 4 +-
.../testGeometryNodeStackExtension.cc | 8 ++--
Stack/History/testHistoryStack.cc | 2 +-
Stack/History/testHistoryView.cc | 33 --------------
do-copyright.py | 2 +-
13 files changed, 41 insertions(+), 79 deletions(-)
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index 5840bc597..1c5c07cd4 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -15,7 +15,6 @@
#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>
@@ -28,6 +27,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/process/interaction_counter/InteractionCounter.hpp>
#include <iomanip>
#include <iostream>
@@ -107,8 +107,8 @@ int main(int argc, char** argv) {
auto const observationHeight = 1.4_km + builder.getEarthRadius();
auto const injectionHeight = 112.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
- sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
- si::detail::static_pow<2>(injectionHeight));
+ sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+ static_pow<2>(injectionHeight));
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
showerCore +
@@ -159,10 +159,9 @@ int main(int argc, char** argv) {
cut.ShowResults();
em_continuous.ShowResults();
observationLevel.ShowResults();
- cout << "Cascade energy cut: " << EAS.GetEnergyCut() / 1_GeV << " GeV" << endl;
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.GetEnergyLost() +
- observationLevel.GetEnergyGround() + EAS.GetEnergyCut();
+ observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index cfdc75edd..5ceb0de05 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -240,10 +240,9 @@ int main(int argc, char** argv) {
cut.ShowResults();
em_continuous.ShowResults();
observationLevel.ShowResults();
- cout << "Cascade energy cut: " << EAS.GetEnergyCut() / 1_GeV << " GeV" << endl;
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.GetEnergyLost() +
- observationLevel.GetEnergyGround() + EAS.GetEnergyCut();
+ observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 5e2e0f0d8..b520acddd 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -104,8 +104,6 @@ namespace corsika::cascade {
#endif
}
- corsika::units::si::HEPEnergyType GetEnergyCut() const { return energy_cut_; }
-
/**
* set the nodes for all particles on the stack according to their numerical
* position
diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h
index fa756414a..55740c560 100644
--- a/Framework/Cascade/testCascade.h
+++ b/Framework/Cascade/testCascade.h
@@ -11,6 +11,7 @@
#include <corsika/environment/Environment.h>
#include <corsika/stack/CombinedStack.h>
+#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/node/GeometryNodeStackExtension.h>
#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
diff --git a/Framework/StackInterface/StackIteratorInterface.h b/Framework/StackInterface/StackIteratorInterface.h
index a7a65bd14..2229c17bc 100644
--- a/Framework/StackInterface/StackIteratorInterface.h
+++ b/Framework/StackInterface/StackIteratorInterface.h
@@ -330,8 +330,6 @@ namespace corsika::stack {
StackType>& rhs) const {
return index_ != rhs.index_;
}
- bool operator==(const StackIteratorInterface<TStackData, TParticleInterface, StackType>& rhs) const { return index_ == rhs.index_; }
- bool operator!=(const StackIteratorInterface<TStackData, TParticleInterface, StackType>& rhs) const { return index_ != rhs.index_; }
const ParticleInterfaceType& operator*() const {
return static_cast<const ParticleInterfaceType&>(*this);
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
index 57315bf4e..4705051d4 100644
--- a/Processes/ParticleCut/ParticleCut.h
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -42,7 +42,6 @@ namespace corsika::process {
return units::si::meter * std::numeric_limits<double>::infinity();
}
-
units::si::HEPEnergyType GetECut() const { return fECut; }
units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
@@ -62,8 +61,6 @@ namespace corsika::process {
bool ParticleIsEmParticle(particles::Code) const;
bool ParticleIsInvisible(particles::Code) const;
-
-
};
} // namespace particle_cut
} // namespace corsika::process
diff --git a/Processes/ParticleCut/testParticleCut.cc b/Processes/ParticleCut/testParticleCut.cc
index 05805cd92..2092636fa 100644
--- a/Processes/ParticleCut/testParticleCut.cc
+++ b/Processes/ParticleCut/testParticleCut.cc
@@ -55,7 +55,7 @@ TEST_CASE("ParticleCut", "[processes]") {
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
// view on secondary particles
- setup::StackView view{particle};
+ setup::StackView view(particle);
// ref. to primary particle through the secondary view.
// only this way the secondary view is populated
auto projectile = view.GetProjectile();
@@ -70,7 +70,7 @@ TEST_CASE("ParticleCut", "[processes]") {
cut.DoSecondaries(view);
- CHECK(view.GetSize() == 9);
+ CHECK(view.getEntries() == 9);
}
SECTION("cut on particle type: em") {
@@ -85,7 +85,7 @@ TEST_CASE("ParticleCut", "[processes]") {
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
// view on secondary particles
- corsika::stack::SecondaryView view(particle);
+ corsika::setup::StackView view(particle);
// ref. to primary particle through the secondary view.
// only this way the secondary view is populated
auto projectile = view.GetProjectile();
@@ -100,7 +100,7 @@ TEST_CASE("ParticleCut", "[processes]") {
cut.DoSecondaries(view);
- CHECK(view.GetSize() == 9);
+ CHECK(view.getEntries() == 9);
}
SECTION("cut low energy") {
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index 064e146ea..b2417b426 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -49,48 +49,51 @@ namespace corsika::process::proposal {
}
template <>
- corsika::process::EProcessReturn Interaction::DoInteraction(
- setup::StackView::StackIterator& vP) {
+ corsika::process::EProcessReturn Interaction::DoInteraction(setup::StackView& view) {
using namespace corsika::units::si; // required for operator::_MeV
- if (CanInteract(vP.GetPID())) {
+ auto const projectile = view.GetProjectile();
+
+ if (CanInteract(projectile.GetPID())) {
// Get or build corresponding calculators
- auto c = GetCalculator(vP, calc);
+ auto c = GetCalculator(projectile, calc);
// Get the rates of the interaction types for every component.
std::uniform_real_distribution<double> distr(0., 1.);
// sample a interaction-type, loss and component
- auto rates = get<eINTERACTION>(c->second)->Rates(vP.GetEnergy() / 1_MeV);
+ auto rates = get<eINTERACTION>(c->second)->Rates(projectile.GetEnergy() / 1_MeV);
auto [type, comp_ptr, v] = get<eINTERACTION>(c->second)->SampleLoss(
- vP.GetEnergy() / 1_MeV, rates, distr(fRNG));
+ projectile.GetEnergy() / 1_MeV, rates, distr(fRNG));
// Read how much random numbers are required to calculate the secondaries.
// Calculate the secondaries and deploy them on the corsika stack.
auto rnd =
vector<double>(get<eSECONDARIES>(c->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 vP_dir = vP.GetDirection();
- auto d = vP_dir.GetComponents();
+ auto point = PROPOSAL::Vector3D(projectile.GetPosition().GetX() / 1_cm,
+ projectile.GetPosition().GetY() / 1_cm,
+ projectile.GetPosition().GetZ() / 1_cm);
+ auto projectile_dir = projectile.GetDirection();
+ auto d = projectile_dir.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.);
+ v * projectile.GetEnergy() / 1_MeV, 0.);
auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries(
- vP.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd);
+ projectile.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(vP_dir.GetCoordinateSystem(), vec);
+ auto p =
+ corsika::stack::MomentumVector(projectile_dir.GetCoordinateSystem(), 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()));
+ view.AddSecondary(
+ make_tuple(sec_code, E, p, projectile.GetPosition(), projectile.GetTime()));
}
}
return process::EProcessReturn::eOk;
@@ -98,13 +101,13 @@ namespace corsika::process::proposal {
template <>
corsika::units::si::GrammageType Interaction::GetInteractionLength(
- setup::Stack::StackIterator const& vP) {
+ setup::Stack::StackIterator const& projectile) {
using namespace corsika::units::si; // required for operator::_MeV
- if (CanInteract(vP.GetPID())) {
- auto c = GetCalculator(vP, calc);
- return get<eINTERACTION>(c->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
- (1_cm * 1_cm);
+ if (CanInteract(projectile.GetPID())) {
+ auto c = GetCalculator(projectile, calc);
+ return get<eINTERACTION>(c->second)->MeanFreePath(projectile.GetEnergy() / 1_MeV) *
+ 1_g / (1_cm * 1_cm);
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index b506dc739..7617e8c27 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -52,8 +52,8 @@ namespace corsika::process::proposal {
//! pair of interaction-type, component and rate, followed by sampling a loss and
//! produce the corresponding secondaries and store them on the particle stack.
//!
- template <typename Particle>
- corsika::process::EProcessReturn DoInteraction(Particle&);
+ template <typename TSecondaryView>
+ corsika::process::EProcessReturn DoInteraction(TSecondaryView&);
//!
//! Calculates the mean free path length
diff --git a/Stack/GeometryNodeStackExtension/testGeometryNodeStackExtension.cc b/Stack/GeometryNodeStackExtension/testGeometryNodeStackExtension.cc
index be14441ec..b566c2c8e 100644
--- a/Stack/GeometryNodeStackExtension/testGeometryNodeStackExtension.cc
+++ b/Stack/GeometryNodeStackExtension/testGeometryNodeStackExtension.cc
@@ -51,7 +51,7 @@ TEST_CASE("GeometryNodeStackExtension", "[stack]") {
TestStack s;
s.AddParticle(std::make_tuple(noData), std::tuple<const int*>{&data});
- CHECK(s.GetSize() == 1);
+ CHECK(s.getEntries() == 1);
}
SECTION("write/read node") {
@@ -60,7 +60,7 @@ TEST_CASE("GeometryNodeStackExtension", "[stack]") {
TestStack s;
auto p = s.AddParticle(std::make_tuple(noData));
p.SetNode(&data);
- CHECK(s.GetSize() == 1);
+ CHECK(s.getEntries() == 1);
const auto pout = s.GetNextParticle();
CHECK(*(pout.GetNode()) == 15);
@@ -77,7 +77,7 @@ TEST_CASE("GeometryNodeStackExtension", "[stack]") {
p.SetNode(&data);
}
- CHECK(s.GetSize() == 99);
+ CHECK(s.getEntries() == 99);
double v = 0;
for (int i = 0; i < 99; ++i) {
auto p = s.GetNextParticle();
@@ -85,6 +85,6 @@ TEST_CASE("GeometryNodeStackExtension", "[stack]") {
p.Delete();
}
CHECK(v == 99 * data);
- CHECK(s.GetSize() == 0);
+ CHECK(s.getEntries() == 0);
}
}
diff --git a/Stack/History/testHistoryStack.cc b/Stack/History/testHistoryStack.cc
index 648ba224b..c068bdb69 100644
--- a/Stack/History/testHistoryStack.cc
+++ b/Stack/History/testHistoryStack.cc
@@ -58,7 +58,7 @@ TEST_CASE("HistoryStackExtension", "[stack]") {
auto p = s.AddParticle(std::tuple<dummy::NoData>{noData});
SECTION("add lone particle") {
- CHECK(s.GetSize() == 1);
+ CHECK(s.getEntries() == 1);
EvtPtr evt = p.GetEvent();
CHECK(evt == nullptr);
diff --git a/Stack/History/testHistoryView.cc b/Stack/History/testHistoryView.cc
index de0f87d77..588537198 100644
--- a/Stack/History/testHistoryView.cc
+++ b/Stack/History/testHistoryView.cc
@@ -165,39 +165,6 @@ TEST_CASE("HistoryStackExtension", "[stack]") {
CHECK(count_generations(sec.GetEvent().get()) == 3);
}
-
- /*
- Now, let's perform some history reading and checking based on p3
-
- p3 should have 20 secondaries, and a projectile (with 15
- secondaries), with another projectil (with 10 secondaries), with
- antother projectil (5 secondaries), with NO parent
-
-
- {
- auto test_ev3 = p3.GetEvent();
- auto test_sec3 = test_ev3->secondaries();
- CHECK(test_sec3.size() == 20);
-
- auto test_proj3 = test_ev3->projectile(s.begin());
- CHECK(test_proj3.GetEvent() == ev3);
- auto test_ev2 = test_ev3->parentEvent();
- auto test_sec2 = test_ev2->secondaries();
- CHECK(test_sec2.size() == 15);
-
- auto test_proj2 = test_ev2->projectile(s.begin());
- CHECK(test_proj2.GetEvent() == ev2);
- auto test_ev1 = test_ev2->parentEvent();
- auto test_sec1 = test_ev1->secondaries();
- CHECK(test_sec1.size() == 10);
-
- auto test_proj1 = test_ev1->projectile(s.begin());
- CHECK(test_proj1.GetEvent() == ev1);
-
- CHECK(test_proj1.GetEvent()->parentEvent() == ev0);
- CHECK(test_proj1.GetEvent()->parentEvent()->parentEvent() == nullptr);
- }
- */
}
SECTION("also test projectile access") {
diff --git a/do-copyright.py b/do-copyright.py
index 4c2b29cc7..5079a3205 100755
--- a/do-copyright.py
+++ b/do-copyright.py
@@ -21,7 +21,7 @@ Debug settings are 0: nothing, 1: checking, 2: filesystem
"""
Debug = 0
-excludeDirs = ["ThirdParty", "git", "build", "install", "PROPOSAL"]
+excludeDirs = ["ThirdParty", "git", "build", "install", "PROPOSAL", "inexlib_rio"]
excludeFiles = ['PhysicalConstants.h','CorsikaFenvOSX.cc', 'sgn.h']
extensions = [".cc", ".h", ".test"]
--
GitLab