IAP GITLAB

Commit d3e58cb7 authored by Ralf Ulrich's avatar Ralf Ulrich

fixed, improvements and switch off debugging output

parent c4b8af2e
Pipeline #3537 passed with stages
in 37 minutes and 48 seconds
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
......
...@@ -121,11 +121,17 @@ namespace corsika { ...@@ -121,11 +121,17 @@ namespace corsika {
ContinuousProcessStepLength const continuousMaxStep = ContinuousProcessStepLength const continuousMaxStep =
sequence_.getMaxStepLength(vParticle, step); sequence_.getMaxStepLength(vParticle, step);
LengthType const continuous_max_dist = continuousMaxStep; LengthType const continuous_max_dist = continuousMaxStep;
ContinuousProcessIndex const limitingId = continuousMaxStep; ContinuousProcessIndex limitingId;
// take minimum of geometry, interaction, decay for next step // take minimum of geometry, interaction, decay for next step
auto const min_distance = LengthType const min_discrete = std::min(distance_interact, distance_decay);
std::min({distance_interact, distance_decay, continuous_max_dist, geomMaxLength}); LengthType const min_non_continuous = std::min(min_discrete, geomMaxLength);
LengthType const min_distance = std::min(min_non_continuous, continuous_max_dist);
if (continuous_max_dist < min_non_continuous) {
limitingId =
continuousMaxStep; // the current step IS limited by known continuous process
}
CORSIKA_LOG_DEBUG( CORSIKA_LOG_DEBUG(
"transport particle by : {} m " "transport particle by : {} m "
...@@ -157,51 +163,43 @@ namespace corsika { ...@@ -157,51 +163,43 @@ namespace corsika {
} }
return; return;
} }
if (continuous_max_dist < min_non_continuous) {
return; // there is nothing further further
}
CORSIKA_LOG_DEBUG("sth. happening before geometric limit ? {}", CORSIKA_LOG_DEBUG("sth. happening before geometric limit ? {}",
((min_distance < geomMaxLength) ? "yes" : "no")); ((min_distance < geomMaxLength) ? "yes" : "no"));
if (min_distance < geomMaxLength) { // interaction to happen within geometric limit if (geomMaxLength < min_discrete) {
// geometric / tracking limit
// check whether decay or interaction limits this step the if (nextVol != currentLogicalNode) {
// outcome of decay or interaction MAY be a) new particles in // boundary crossing, step is limited by volume boundary
// secondaries, b) the projectile particle deleted (or
// changed)
TStackView secondaries(vParticle); CORSIKA_LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
if (min_distance < continuous_max_dist) { if (nextVol == environment_.getUniverse().get()) {
CORSIKA_LOG_DEBUG(
"particle left physics world, is now in unknown space -> delete");
vParticle.erase();
}
vParticle.setNode(nextVol);
/* /*
Create SecondaryView object on Stack. The data container doBoundary may delete the particle (or not)
remains untouched and identical, and 'projectil' is identical
to 'vParticle' above this line. However,
projectil.AddSecondaries populate the SecondaryView, which can
then be used afterwards for further processing. Thus: it is
important to use projectle/view (and not vParticle) for Interaction,
and Decay!
*/
[[maybe_unused]] auto projectile = secondaries.getProjectile(); caveat: any changes to vParticle, or even the production
of new secondaries is currently not passed to ParticleCut,
if (distance_interact < distance_decay) { thus, particles outside the desired phase space may be produced.
interaction(secondaries);
} else {
decay(secondaries);
// make sure particle actually did decay if it should have done so
if (secondaries.getSize() == 1 &&
projectile.getPID() == secondaries.getNextParticle().getPID())
throw std::runtime_error(fmt::format("Particle {} decays into itself!",
get_name(projectile.getPID())));
}
sequence_.doSecondaries(secondaries); \todo: this must be fixed.
vParticle.erase(); */
} else { // step-length limitation within volume
CORSIKA_LOG_DEBUG("step-length limitation"); sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
// no further physics happens here. just proceed to next step. return;
} }
CORSIKA_LOG_DEBUG("step limit reached. nothing further happens.");
[[maybe_unused]] auto const assertion = [&] { [[maybe_unused]] auto const assertion = [&] {
auto const* numericalNodeAfterStep = auto const* numericalNodeAfterStep =
environment_.getUniverse()->getContainingNode(vParticle.getPosition()); environment_.getUniverse()->getContainingNode(vParticle.getPosition());
...@@ -214,31 +212,42 @@ namespace corsika { ...@@ -214,31 +212,42 @@ namespace corsika {
assert(assertion()); // numerical and logical nodes should match, since assert(assertion()); // numerical and logical nodes should match, since
// we did not cross any volume boundary // we did not cross any volume boundary
} else { // boundary crossing, step is limited by volume boundary // step length limit
return;
if (nextVol != currentLogicalNode) { }
CORSIKA_LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
if (nextVol == environment_.getUniverse().get()) {
CORSIKA_LOG_DEBUG(
"particle left physics world, is now in unknown space -> delete");
vParticle.erase();
}
vParticle.setNode(nextVol);
/*
doBoundary may delete the particle (or not)
caveat: any changes to vParticle, or even the production // interaction or decay to happen in this step
of new secondaries is currently not passed to ParticleCut, // the outcome of decay or interaction MAY be a) new particles in
thus, particles outside the desired phase space may be produced. // secondaries, b) the projectile particle deleted (or
// changed)
\todo: this must be fixed. TStackView secondaries(vParticle);
*/
sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol); /*
} Create SecondaryView object on Stack. The data container
remains untouched and identical, and 'projectil' is identical
to 'vParticle' above this line. However,
projectil.AddSecondaries populate the SecondaryView, which can
then be used afterwards for further processing. Thus: it is
important to use projectle/view (and not vParticle) for Interaction,
and Decay!
*/
[[maybe_unused]] auto projectile = secondaries.getProjectile();
if (distance_interact < distance_decay) {
interaction(secondaries);
} else {
decay(secondaries);
// make sure particle actually did decay if it should have done so
if (secondaries.getSize() == 1 &&
projectile.getPID() == secondaries.getNextParticle().getPID())
throw std::runtime_error(fmt::format("Particle {} decays into itself!",
get_name(projectile.getPID())));
} }
sequence_.doSecondaries(secondaries);
vParticle.erase();
} }
template <typename TTracking, typename TProcessList, typename TStack, template <typename TTracking, typename TProcessList, typename TStack,
......
...@@ -35,11 +35,11 @@ namespace corsika { ...@@ -35,11 +35,11 @@ namespace corsika {
GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0)); GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1)); GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
CORSIKA_LOG_INFO("pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2", CORSIKA_LOG_TRACE("pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2",
vTrack.getPosition(0).getCoordinates() / 1_m, vTrack.getPosition(0).getCoordinates() / 1_m,
vTrack.getPosition(1).getCoordinates() / 1_m, vTrack.getPosition(1).getCoordinates() / 1_m,
grammageStart / 1_g * square(1_cm), grammageStart / 1_g * square(1_cm),
grammageEnd / 1_g * square(1_cm)); grammageEnd / 1_g * square(1_cm));
// Note: particle may go also "upward", thus, grammageEnd<grammageStart // Note: particle may go also "upward", thus, grammageEnd<grammageStart
const int binStart = std::ceil(grammageStart / dX_); const int binStart = std::ceil(grammageStart / dX_);
......
...@@ -24,8 +24,8 @@ ...@@ -24,8 +24,8 @@
namespace corsika { namespace corsika {
template <typename TStack> template <typename TStack>
inline StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack, inline StackInspector<TStack>::StackInspector(int const vNStep, bool const vReportStack,
const HEPEnergyType vE0) HEPEnergyType const vE0)
: StackProcess<StackInspector<TStack>>(vNStep) : StackProcess<StackInspector<TStack>>(vNStep)
, ReportStack_(vReportStack) , ReportStack_(vReportStack)
, E0_(vE0) , E0_(vE0)
...@@ -35,12 +35,12 @@ namespace corsika { ...@@ -35,12 +35,12 @@ namespace corsika {
inline StackInspector<TStack>::~StackInspector() {} inline StackInspector<TStack>::~StackInspector() {}
template <typename TStack> template <typename TStack>
inline void StackInspector<TStack>::doStack(const TStack& vS) { inline void StackInspector<TStack>::doStack(TStack const& vS) {
[[maybe_unused]] int i = 0; [[maybe_unused]] int i = 0;
HEPEnergyType Etot = 0_GeV; HEPEnergyType Etot = 0_GeV;
for (const auto& iterP : vS) { for (auto const& iterP : vS) {
HEPEnergyType E = iterP.getEnergy(); HEPEnergyType E = iterP.getEnergy();
Etot += E; Etot += E;
if (ReportStack_) { if (ReportStack_) {
...@@ -60,7 +60,7 @@ namespace corsika { ...@@ -60,7 +60,7 @@ namespace corsika {
} }
auto const now = std::chrono::system_clock::now(); auto const now = std::chrono::system_clock::now();
const std::chrono::duration<double> elapsed_seconds = now - StartTime_; std::chrono::duration<double> const elapsed_seconds = now - StartTime_;
std::time_t const now_time = std::chrono::system_clock::to_time_t(now); std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
auto const dE = E0_ - Etot; auto const dE = E0_ - Etot;
if (dE < dE_threshold_) return; if (dE < dE_threshold_) return;
......
...@@ -82,7 +82,7 @@ namespace corsika::pythia8 { ...@@ -82,7 +82,7 @@ namespace corsika::pythia8 {
inline void Decay::setHandleDecay(Code const vParticleCode) { inline void Decay::setHandleDecay(Code const vParticleCode) {
handleAllDecays_ = false; handleAllDecays_ = false;
CORSIKA_LOG_INFO("Pythia::Decay: set to handle decay of {} ", vParticleCode); CORSIKA_LOG_DEBUG("Pythia::Decay: set to handle decay of {} ", vParticleCode);
if (Decay::canHandleDecay(vParticleCode)) if (Decay::canHandleDecay(vParticleCode))
handledDecays_.insert(vParticleCode); handledDecays_.insert(vParticleCode);
else else
...@@ -106,12 +106,12 @@ namespace corsika::pythia8 { ...@@ -106,12 +106,12 @@ namespace corsika::pythia8 {
} }
inline void Decay::setUnstable(Code const pCode) { inline void Decay::setUnstable(Code const pCode) {
CORSIKA_LOG_INFO("Pythia::Decay: setting {} unstable..", pCode); CORSIKA_LOG_DEBUG("Pythia::Decay: setting {} unstable..", pCode);
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true); Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
} }
inline void Decay::setStable(Code const pCode) { inline void Decay::setStable(Code const pCode) {
CORSIKA_LOG_INFO("Pythia::Decay: setting {} stable..", pCode); CORSIKA_LOG_DEBUG("Pythia::Decay: setting {} stable..", pCode);
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false); Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
} }
...@@ -122,8 +122,8 @@ namespace corsika::pythia8 { ...@@ -122,8 +122,8 @@ namespace corsika::pythia8 {
inline bool Decay::canDecay(Code const pCode) { inline bool Decay::canDecay(Code const pCode) {
bool const ans = bool const ans =
Pythia8::Pythia::particleData.canDecay(static_cast<int>(get_PDG(pCode))); Pythia8::Pythia::particleData.canDecay(static_cast<int>(get_PDG(pCode)));
CORSIKA_LOG_INFO("Pythia::Decay: checking if particle: {} can decay in PYTHIA? {} ", CORSIKA_LOG_DEBUG("Pythia::Decay: checking if particle: {} can decay in PYTHIA? {} ",
pCode, ans); pCode, ans);
return ans; return ans;
} }
...@@ -153,13 +153,13 @@ namespace corsika::pythia8 { ...@@ -153,13 +153,13 @@ namespace corsika::pythia8 {
TimeType const t0 = get_lifetime(pid); TimeType const t0 = get_lifetime(pid);
auto const lifetime = gamma * t0; auto const lifetime = gamma * t0;
CORSIKA_LOG_INFO("Pythia::Decay: code: {}", particle.getPID()); CORSIKA_LOG_TRACE("Pythia::Decay: code: {}", particle.getPID());
CORSIKA_LOG_INFO("Pythia::Decay: MinStep: t0: {}", t0); CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: t0: {}", t0);
CORSIKA_LOG_INFO("Pythia::Decay: MinStep: energy: {} GeV", E / 1_GeV); CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: energy: {} GeV", E / 1_GeV);
CORSIKA_LOG_INFO("Pythia::Decay: momentum: {} GeV", CORSIKA_LOG_TRACE("Pythia::Decay: momentum: {} GeV",
particle.getMomentum().getComponents() / 1_GeV); particle.getMomentum().getComponents() / 1_GeV);
CORSIKA_LOG_INFO("Pythia::Decay: MinStep: gamma: {}", gamma); CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: gamma: {}", gamma);
CORSIKA_LOG_INFO("Pythia::Decay: MinStep: tau: {} ", lifetime); CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: tau: {} ", lifetime);
return lifetime; return lifetime;
} else } else
...@@ -209,7 +209,7 @@ namespace corsika::pythia8 { ...@@ -209,7 +209,7 @@ namespace corsika::pythia8 {
if (!Pythia8::Pythia::next()) if (!Pythia8::Pythia::next())
throw std::runtime_error("Pythia::Decay: decay failed!"); throw std::runtime_error("Pythia::Decay: decay failed!");
else else
CORSIKA_LOG_INFO("Pythia::Decay: particles after decay: {} ", event.size()); CORSIKA_LOG_DEBUG("Pythia::Decay: particles after decay: {} ", event.size());
if (print_listing_) { if (print_listing_) {
// list final state // list final state
...@@ -227,9 +227,10 @@ namespace corsika::pythia8 { ...@@ -227,9 +227,10 @@ namespace corsika::pythia8 {
FourVector const fourMomRest{Erest, pRest}; FourVector const fourMomRest{Erest, pRest};
auto const fourMomLab = boost.fromCoM(fourMomRest); auto const fourMomLab = boost.fromCoM(fourMomRest);
CORSIKA_LOG_INFO("particle: id={} momentum={} energy={} ", pyId, CORSIKA_LOG_TRACE(
fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV, "particle: id={} momentum={} energy={} ", pyId,
fourMomLab.getTimeLikeComponent()); fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
fourMomLab.getTimeLikeComponent());
view.addSecondary(std::make_tuple(pyId, fourMomLab.getTimeLikeComponent(), view.addSecondary(std::make_tuple(pyId, fourMomLab.getTimeLikeComponent(),
fourMomLab.getSpaceLikeComponents(), decayPoint, fourMomLab.getSpaceLikeComponents(), decayPoint,
......
...@@ -18,6 +18,8 @@ namespace corsika { ...@@ -18,6 +18,8 @@ namespace corsika {
class ContinuousProcessIndex { class ContinuousProcessIndex {
public: public:
ContinuousProcessIndex()
: id_(-1) {} // default
ContinuousProcessIndex(int const id) ContinuousProcessIndex(int const id)
: id_(id) {} : id_(id) {}
void setIndex(int const id) { id_ = id; } void setIndex(int const id) { id_ = id; }
......
...@@ -10,6 +10,14 @@ import os ...@@ -10,6 +10,14 @@ import os
import sys import sys
import re import re
do_progress = False
try:
from progress.bar import ChargingBar
do_progress = True
except ImportError as e:
do_progress = False
# no progress bar
parser = argparse.ArgumentParser(description=__doc__) parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--apply', action="store_true", parser.add_argument('--apply', action="store_true",
help="Apply clang-format to files which need changes.") help="Apply clang-format to files which need changes.")
...@@ -85,17 +93,26 @@ version = subp.check_output(cmd.split() + ["--version"]).decode("utf-8") ...@@ -85,17 +93,26 @@ version = subp.check_output(cmd.split() + ["--version"]).decode("utf-8")
print (version) print (version)
print ("Note: the clang-format version has an impact on the result. Make sure you are consistent with current CI. Consider \'--docker\' option.") print ("Note: the clang-format version has an impact on the result. Make sure you are consistent with current CI. Consider \'--docker\' option.")
bar = None
if do_progress:
bar = ChargingBar('Processing', max=len(filelist))
if args.apply: if args.apply:
for filename in filelist: for filename in filelist:
if bar: bar.next()
subp.check_call(cmd.split() + ["-i", filename]) subp.check_call(cmd.split() + ["-i", filename])
if bar: bar.finish()
else: else:
# only print files which need formatting # only print files which need formatting
files_need_formatting = 0 files_need_formatting = 0
for filename in filelist: for filename in filelist:
if bar: bar.next()
a = open(filename, "rb").read() a = open(filename, "rb").read()
b = subp.check_output(cmd.split() + [filename]) b = subp.check_output(cmd.split() + [filename])
if a != b: if a != b:
files_need_formatting += 1 files_need_formatting += 1
print(filename) print(filename)
if bar: bar.finish()
sys.exit(1 if files_need_formatting > 0 else 0) sys.exit(1 if files_need_formatting > 0 else 0)
...@@ -6,8 +6,6 @@ ...@@ -6,8 +6,6 @@
* the license. * the license.
*/ */
#define DEBUG 1
/* clang-format off */ /* clang-format off */
// InteractionCounter used boost/histogram, which // InteractionCounter used boost/histogram, which
// fails if boost/type_traits have been included before. Thus, we have // fails if boost/type_traits have been included before. Thus, we have
...@@ -96,7 +94,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; ...@@ -96,7 +94,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) { int main(int argc, char** argv) {
corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v"); corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
logging::set_level(logging::level::trace); logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("vertical_EAS"); CORSIKA_LOG_INFO("vertical_EAS");
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment