From c96767ee8747aa3d7bfbc9da155ef1ebaaa76722 Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@tu-dortmund.de>
Date: Wed, 30 Sep 2020 15:48:56 +0000
Subject: [PATCH] add further doc and some minor fixes in particle cut
---
Documentation/Examples/vertical_EAS.cc | 11 -----
Processes/ParticleCut/ParticleCut.cc | 2 +-
Processes/ParticleCut/ParticleCut.h | 4 +-
Processes/Proposal/ContinuousProcess.cc | 66 +++++++++++++++++++------
Processes/Proposal/ContinuousProcess.h | 9 +---
Processes/Proposal/Interaction.cc | 20 +++++++-
Processes/Proposal/Interaction.h | 26 ++++++++--
7 files changed, 99 insertions(+), 39 deletions(-)
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 26b259c84..738bd73e4 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -195,7 +195,6 @@ int main(int argc, char** argv) {
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
- process::energy_loss::EnergyLoss eLoss(showerAxis);
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
@@ -205,10 +204,6 @@ int main(int argc, char** argv) {
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
- /* process::conex_source_cut::CONEXSourceCut conexSource( */
- /* center, showerAxis, t, injectionHeight, E0, */
- /* particles::GetPDG(particles::Code::Proton)); */
-
// assemble all processes into an ordered process list
auto sibyllSequence = sibyllNucCounted << sibyllCounted;
@@ -216,8 +211,6 @@ int main(int argc, char** argv) {
55_GeV);
auto decaySequence = decayPythia << decaySibyll;
- //auto sequence = switchProcess << reset_particle_mass << decaySequence << conexSource
- // << longprof << eLoss << cut << observationLevel;
auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << cut << em_continuous
<< longprof << observationLevel;
@@ -231,16 +224,12 @@ int main(int argc, char** argv) {
EAS.Run();
- eLoss.PrintProfile(); // print longitudinal profile
- /* conexSource.SolveCE(); */
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;
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc
index defc68bad..e59addccc 100644
--- a/Processes/ParticleCut/ParticleCut.cc
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -84,7 +84,7 @@ namespace corsika::process {
ParticleCut::ParticleCut(const units::si::HEPEnergyType eCut, bool discardEm, bool discardInv)
: eCut_(eCut)
- , discardEm_(cutEm)
+ , discardEm_(discardEm)
, discardInv_(discardInv) {
emEnergy_ = 0_GeV;
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
index 68639471a..0edab21e3 100644
--- a/Processes/ParticleCut/ParticleCut.h
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -18,8 +18,8 @@ namespace corsika::process {
class ParticleCut : public process::SecondariesProcess<ParticleCut> {
units::si::HEPEnergyType const eCut_;
- bool cutEm_;
- bool cutInv_;
+ bool discardEm_;
+ bool discardInv_;
units::si::HEPEnergyType energy_ = 0 * units::si::electronvolt;
units::si::HEPEnergyType emEnergy_ = 0 * units::si::electronvolt;
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 49427f03e..eb5f8d3f7 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -23,10 +23,19 @@ namespace corsika::process::proposal {
void ContinuousProcess::BuildCalculator(particles::Code code,
environment::NuclearComposition const& comp) {
+ // search crosssection builder for given particle
auto p_cross = cross.find(code);
if (p_cross == cross.end())
throw std::runtime_error("PROPOSAL could not find corresponding builder");
+
+ // interpolate the crosssection for given media and energy cut. These may
+ // take some minutes if you have to build the tables and cannot read the
+ // from disk
auto c = p_cross->second(media.at(&comp), cut);
+
+ // Build displacement integral and scattering object and interpolate them too and
+ // saved in the calc map by a key build out of a hash of composed of the component and
+ // particle code.
auto disp = PROPOSAL::make_displacement(c, true);
auto scatter = PROPOSAL::make_scattering("highland", particle[code], media.at(&comp));
calc[std::make_pair(&comp, code)] =
@@ -38,31 +47,35 @@ namespace corsika::process::proposal {
particle_cut::ParticleCut& _cut)
: ProposalProcessBase(_env, _cut) {}
- template <>
- HEPEnergyType ContinuousProcess::TotalEnergyLoss(setup::Stack::ParticleType const& vP,
- GrammageType const& vDX) {
- auto c = GetCalculator(vP, calc);
- return vP.GetEnergy() - get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
- vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm) *
- 1_MeV;
- }
-
template <>
void ContinuousProcess::Scatter(setup::Stack::ParticleType& vP,
HEPEnergyType const& loss,
GrammageType const& grammage) {
+ // Get or build corresponding calculators
auto c = GetCalculator(vP, calc);
+
+ // Cast corsika vector to proposal vector
auto d = vP.GetDirection().GetComponents();
auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
d.GetZ().magnitude());
- auto E_f = vP.GetEnergy() - loss; // final energy
+
+ auto E_f = vP.GetEnergy() - loss;
+
+ // draw random numbers required for scattering process
std::uniform_real_distribution<double> distr(0., 1.);
auto rnd = array<double, 4>();
for (auto& it : rnd) it = distr(fRNG);
+
+ // calculate deflection based on particle energy, loss
auto [mean_dir, final_dir] = get<SCATTERING>(c->second)->Scatter(
grammage / 1_g * square(1_cm), vP.GetEnergy() / 1_MeV, E_f / 1_MeV, direction,
rnd);
+
+ // TODO: neglect mean direction deflection because Trajectory is a const ref
(void)mean_dir;
+
+ // update particle direction after continuous loss caused by multiple
+ // scattering
auto vec = corsika::geometry::QuantityVector(
final_dir.GetX() * E_f, final_dir.GetY() * E_f, final_dir.GetZ() * E_f);
vP.SetMomentum(corsika::stack::MomentumVector(
@@ -74,11 +87,26 @@ namespace corsika::process::proposal {
EProcessReturn ContinuousProcess::DoContinuous(setup::Stack::ParticleType& vP,
setup::Trajectory const& vT) {
if (!CanInteract(vP.GetPID())) return process::EProcessReturn::eOk;
+
+ // calculate passed grammage
auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
- auto energy_loss = TotalEnergyLoss(vP, dX);
- if (vP.GetChargeNumber() != 0) Scatter(vP, energy_loss, dX);
- vP.SetEnergy(vP.GetEnergy() - energy_loss);
- if (vP.GetEnergy() <= cut.GetECut()) return process::EProcessReturn::eParticleAbsorbed;
+
+ // Get or build corresponding track integral calculator and solve the
+ // integral
+ auto c = GetCalculator(vP, calc);
+ auto final_energy = get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
+ vP.GetEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
+ 1_MeV;
+
+ // if the particle has a charge take multiple scattering into account
+ if (vP.GetChargeNumber() != 0 && vP.GetEnergy() > final_energy)
+ Scatter(vP, vP.GetEnergy() - final_energy, dX);
+
+ // Update the energy and absorbe the particle if it's below the energy
+ // threshold, because it will no longer propagated.
+ vP.SetEnergy(final_energy);
+ if (std::fabs((vP.GetEnergy() - cut.GetECut() )/ cut.GetECut()) < 0.01)
+ return process::EProcessReturn::eParticleAbsorbed;
vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
return process::EProcessReturn::eOk;
}
@@ -88,12 +116,22 @@ namespace corsika::process::proposal {
setup::Stack::ParticleType const& vP, setup::Trajectory const& vT) {
if (!CanInteract(vP.GetPID()))
return units::si::meter * std::numeric_limits<double>::infinity();
+
+ // Limit the step size of a conitnous loss. The maximal continuous loss seems to be a
+ // hyper parameter which must be adjusted.
auto energy_lim = 0.9 * vP.GetEnergy();
+
+ // If energy lim is below the cut it couldn't be descriped by these process. So the
+ // energy_lim is replaced by the cut.
if (cut.GetECut() > energy_lim) energy_lim = cut.GetECut();
+
+ // solving the track integral for giving energy lim
auto c = GetCalculator(vP, calc);
auto grammage = get<DISPLACEMENT>(c->second)->SolveTrackIntegral(
vP.GetEnergy() / 1_MeV, energy_lim / 1_MeV) *
1_g / square(1_cm);
+
+ // return it in distance aequivalent
return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage);
}
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index c50ac8d54..9449d9c74 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -43,19 +43,14 @@ namespace corsika::process::proposal {
void BuildCalculator(particles::Code, environment::NuclearComposition const&) final;
public:
+
+ //!
//! Produces the continuous loss calculator for leptons based on nuclear
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
ContinuousProcess(TEnvironment const&, particle_cut::ParticleCut&);
- //!
- //! Calculates the energy of the particle which has passed the given
- //! grammage.
- //!
- template <typename Particle>
- corsika::units::si::HEPEnergyType TotalEnergyLoss(
- Particle const&, corsika::units::si::GrammageType const&);
//!
//! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index 2ac0e7b70..b88653360 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -31,10 +31,20 @@ namespace corsika::process::proposal {
void Interaction::BuildCalculator(particles::Code code,
environment::NuclearComposition const& comp) {
+ // search crosssection builder for given particle
auto p_cross = cross.find(code);
if (p_cross == cross.end())
throw std::runtime_error("PROPOSAL could not find corresponding builder");
+
+ // interpolate the crosssection for given media and energy cut. These may
+ // take some minutes if you have to build the tables and cannot read the
+ // from disk
auto c = p_cross->second(media.at(&comp), cut);
+
+ // Look which interactions take place and build the corresponding
+ // interaction and secondarie builder. The interaction integral will
+ // interpolated too and saved in the calc map by a key build out of a hash
+ // of composed of the component and particle code.
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
calc[std::make_pair(&comp, code)] = std::make_tuple(
PROPOSAL::make_secondaries(inter_types, particle[code], media.at(&comp)),
@@ -45,11 +55,19 @@ namespace corsika::process::proposal {
corsika::process::EProcessReturn Interaction::DoInteraction(
setup::StackView::StackIterator& vP) {
if (CanInteract(vP.GetPID())) {
- auto c = GetCalculator(vP, calc); // [CrossSections]
+ // Get or build corresponding calculators
+ auto c = GetCalculator(vP, 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<INTERACTION>(c->second)->Rates(vP.GetEnergy() / 1_MeV);
auto [type, comp_ptr, v] = get<INTERACTION>(c->second)->SampleLoss(
vP.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<SECONDARIES>(c->second)->RequiredRandomNumbers(type));
for (auto& it : rnd) it = distr(fRNG);
auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index 990cf1115..c96cf00df 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -23,24 +23,44 @@ namespace corsika::process::proposal {
using namespace corsika::units::si;
+ //!
+ //! Electro-magnetic and gamma stochastic losses produced by proposal. It makes
+ //! use of interpolation tables which are runtime intensive calculation, but can be
+ //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
+ //!
class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
+ enum { SECONDARIES, INTERACTION };
using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
unique_ptr<PROPOSAL::Interaction>>;
- std::unordered_map<calc_key_t, calculator_t, hash> calc;
+ std::unordered_map<calc_key_t, calculator_t, hash>
+ calc; //!< Stores the secondaries and interaction calculators.
+ //!
+ //! Build the secondaries and interaction calculators and add it to calc.
+ //!
void BuildCalculator(particles::Code, environment::NuclearComposition const&) final;
- enum { SECONDARIES, INTERACTION };
-
public:
+ //!
+ //! Produces the stoachastic loss calculator for leptons based on nuclear
+ //! compositions and stochastic description limited by the particle cut.
+ //!
template <typename TEnvironment>
Interaction(TEnvironment const& env, particle_cut::ParticleCut& cut);
+ //!
+ //! Calculate the rates for the different targets and interactions. Sample a
+ //! 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&);
+ //!
+ //! Calculates the mean free path length
+ //!
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const& p);
};
--
GitLab