IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 3cc44628 authored by ralfulrich's avatar ralfulrich
Browse files

fixed clang history option

parent ca180706
No related branches found
No related tags found
No related merge requests found
Showing
with 130 additions and 62 deletions
......@@ -74,8 +74,6 @@ void registerRandomStreams(const int seed) {
int main(int argc, char** argv) {
logging::SetLevel(logging::level::debug);
C8LOG_INFO("vertical_EAS");
if (argc < 4) {
......
......@@ -17,6 +17,8 @@
#include <corsika/stack/SecondaryView.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/logging/Logging.h>
#include <corsika/setup/SetupTrajectory.h>
/* see Issue 161, we need to include SetupStack only because we need
......@@ -61,8 +63,9 @@ namespace corsika::cascade {
template <typename TTracking, typename TProcessList, typename TStack,
/*
TStackView is needed as template parameter because of issue 161 and the
inability of clang to understand "MakeView" so far.
TStackView is needed as explicit template parameter because
of issue 161 and the
inability of clang to understand "stack::MakeView" so far.
*/
typename TStackView = corsika::setup::StackView>
class Cascade {
......@@ -96,9 +99,8 @@ namespace corsika::cascade {
, fTracking(tr)
, fProcessSequence(pl)
, fStack(stack)
, energy_cut_(0 * corsika::units::si::electronvolt) {
std::cout << c8_ascii_ << std::endl;
, count_(0) {
C8LOG_INFO(c8_ascii_);
}
corsika::units::si::HEPEnergyType GetEnergyCut() const { return energy_cut_; }
......@@ -124,13 +126,17 @@ namespace corsika::cascade {
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
C8LOG_TRACE(fmt::format("Stack: {}", fStack.as_string()));
for (const auto& tmp_p : fStack)
std::cout << " test " << tmp_p.GetPID() << std::endl;
count_++;
auto pNext = fStack.GetNextParticle();
std::cout << "========= next: count=" << count_ << ", pid=" << pNext.GetPID()
<< ", stack entries=" << fStack.getEntries()
<< ", stack deleted=" << fStack.getDeleted() << std::endl;
C8LOG_DEBUG(fmt::format(
"============== next particle : count={}, pid={}, "
", stack entries={}"
", stack deleted={}",
count_, pNext.GetPID(), fStack.getEntries(), fStack.getDeleted()));
Step(pNext);
std::cout << "========= stack ============" << std::endl;
fProcessSequence.DoStack(fStack);
}
// do cascade equations, which can put new particles on Stack,
......@@ -145,7 +151,7 @@ namespace corsika::cascade {
* want to call forceInteraction() for the primary interaction.
*/
void forceInteraction() {
std::cout << "forced interaction!" << std::endl;
C8LOG_DEBUG("forced interaction!");
auto vParticle = fStack.GetNextParticle();
TStackView secondaries(vParticle);
interaction(vParticle, secondaries);
......@@ -180,8 +186,10 @@ namespace corsika::cascade {
corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
GrammageType const next_interact = expDist(fRNG);
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
C8LOG_DEBUG(
"total_lambda={} g/cm2, "
", next_interact={} g/cm2",
double((1. / total_inv_lambda) / 1_g * 1_cm * 1_cm), double(next_interact / 1_g * 1_cm * 1_cm));
auto const* currentLogicalNode = vParticle.GetNode();
......@@ -197,7 +205,7 @@ namespace corsika::cascade {
// determine the maximum geometric step length from continuous processes
LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
std::cout << "distance_max=" << distance_max << std::endl;
C8LOG_DEBUG("distance_max={} m", distance_max/1_m);
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime =
......@@ -206,8 +214,10 @@ namespace corsika::cascade {
// sample random exponential decay time
corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
TimeType const next_decay = expDistDecay(fRNG);
std::cout << "total_inv_lifetime=" << total_inv_lifetime
<< ", next_decay=" << next_decay << std::endl;
C8LOG_DEBUG(
"total_lifetime={} s"
", next_decay={} s",
(1/total_inv_lifetime)/1_s, next_decay/1_s);
// convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
......@@ -217,7 +227,7 @@ namespace corsika::cascade {
auto const min_distance =
std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
std::cout << " move particle by : " << min_distance << std::endl;
C8LOG_DEBUG("transport particle by : {} m", min_distance/1_m);
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
......@@ -231,15 +241,14 @@ namespace corsika::cascade {
process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
if (status == process::EProcessReturn::eParticleAbsorbed) {
std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
<< vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl;
energy_cut_ += vParticle.GetEnergy();
C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV);
vParticle.Delete();
return;
}
std::cout << "sth. happening before geometric limit ? "
<< ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
C8LOG_DEBUG("sth. happening before geometric limit ? {}",
((min_distance < geomMaxLength) ? "yes" : "no"));
if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
......@@ -257,7 +266,7 @@ namespace corsika::cascade {
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 (and not vParticle) for Interaction,
important to use projectle/view (and not vParticle) for Interaction,
and Decay!
*/
......@@ -277,15 +286,10 @@ namespace corsika::cascade {
}
fProcessSequence.DoSecondaries(secondaries);
vParticle.Delete(); // todo: this should be reviewed. Where
// exactly are particles best deleted, and
// where they should NOT be
// deleted... maybe Delete function should
// be "protected" and not accessible to physics
vParticle.Delete();
} else { // step-length limitation within volume
std::cout << "step-length limitation" << std::endl;
C8LOG_DEBUG("step-length limitation");
// no extra physics happens here. just proceed to next step.
}
......@@ -299,10 +303,8 @@ namespace corsika::cascade {
};
assert(assertion()); // numerical and logical nodes don't match
} else { // boundary crossing, step is limited by volume boundary
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
} else { // boundary crossing, step is limited by volume boundary
// C8LOG_DEBUG("boundary crossing! next node = {}", int(nextVol));
vParticle.SetNode(nextVol);
/*
DoBoundary may delete the particle (or not)
......@@ -316,7 +318,7 @@ namespace corsika::cascade {
}
auto decay(Particle& particle, TStackView& view) {
std::cout << "decay" << std::endl;
C8LOG_DEBUG("decay");
units::si::InverseTimeType const actual_decay_time =
fProcessSequence.GetTotalInverseLifetime(particle);
......@@ -329,7 +331,7 @@ namespace corsika::cascade {
}
auto interaction(Particle& particle, TStackView& view) {
std::cout << "collide" << std::endl;
C8LOG_DEBUG("collide");
units::si::InverseGrammageType const current_inv_length =
fProcessSequence.GetTotalInverseInteractionLength(particle);
......
......@@ -11,6 +11,7 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/logging/Logging.h>
namespace corsika::stack {
......@@ -37,10 +38,6 @@ namespace corsika::stack {
class CombinedParticleInterface
: public ParticleInterfaceB<ParticleInterfaceA<StackIterator>> {
// template<template <typename> typename _PI>
// template <typename StackDataType, template <typename> typename ParticleInterface>
// template<typename T1, template <typename> typename T2> friend class Stack<T1, T2>;
using PI_C =
CombinedParticleInterface<ParticleInterfaceA, ParticleInterfaceB, StackIterator>;
using PI_A = ParticleInterfaceA<StackIterator>;
......@@ -91,6 +88,11 @@ namespace corsika::stack {
PI_B::SetParticleData(static_cast<PI_B&>(p), vB);
}
///@}
std::string as_string() const {
return fmt::format("[[{}][{}]]",PI_A::as_string(), PI_B::as_string());
}
};
/**
......@@ -113,7 +115,7 @@ namespace corsika::stack {
unsigned int GetSize() const { return Stack1Impl::GetSize(); }
unsigned int GetCapacity() const { return Stack1Impl::GetCapacity(); }
/**
* Function to copy particle at location i1 in stack to i2
*/
......
......@@ -276,6 +276,7 @@ namespace corsika::stack {
*
*/
void Delete(StackIterator p) {
C8LOG_TRACE("SecondaryView::Delete");
if (IsEmpty()) { /*error*/
throw std::runtime_error("Stack, cannot delete entry since size is zero");
}
......@@ -323,6 +324,7 @@ namespace corsika::stack {
* the function will just return false and do nothing.
*/
bool purgeLastIfDeleted() {
C8LOG_TRACE("SecondaryView::purgeLastIfDeleted");
if (!isDeleted(getSize() - 1))
return false; // the last particle is not marked for deletion. Do nothing.
inner_stack_.purge(GetIndexFromIterator(getSize()));
......@@ -367,6 +369,7 @@ namespace corsika::stack {
* performed.
*/
unsigned int GetIndexFromIterator(const unsigned int vI) const {
// this is too much: C8LOG_TRACE("SecondaryView::GetIndexFromIterator({})={}", vI, (vI?indices_[vI-1]:projectile_index_));
if (vI == 0) return projectile_index_;
return indices_[vI - 1];
}
......
......@@ -10,10 +10,12 @@
#include <corsika/stack/StackIteratorInterface.h>
#include <corsika/utl/MetaProgramming.h>
#include <corsika/logging/Logging.h>
#include <stdexcept>
#include <type_traits>
#include <vector>
#include <string>
/**
All classes around management of particles on a stack.
......@@ -85,9 +87,9 @@ namespace corsika::stack {
* stacks, where the inner data container is always a reference
* and cannot be initialized here.
*/
template <typename... Args, typename _ = TStackData,
template <typename... TArgs, typename _ = TStackData,
typename = utl::disable_if<std::is_reference<_>>>
Stack(Args... args)
Stack(TArgs... args)
: data_(args...)
, deleted_(std::vector<bool>(data_.GetSize(), false))
, nDeleted_(0) {}
......@@ -125,7 +127,6 @@ namespace corsika::stack {
friend class StackIteratorInterface<StackDataValueType, MParticleInterface, Stack>;
friend class ConstStackIteratorInterface<StackDataValueType, MParticleInterface,
Stack>;
// friend class SecondaryView<StackDataValueType, MParticleInterface>;
template <typename T1, //=TStackData,
template <typename>
typename M1, //=MParticleInterface,
......@@ -145,8 +146,8 @@ namespace corsika::stack {
unsigned int getDeleted() const { return nDeleted_; }
unsigned int getEntries() const { return getSize() - getDeleted(); }
template <typename... Args>
void Clear(Args... args) {
template <typename... TArgs>
void Clear(TArgs... args) {
data_.Clear(args...);
deleted_ = std::vector<bool>(data_.GetSize(), false);
nDeleted_ = 0;
......@@ -195,7 +196,7 @@ namespace corsika::stack {
for (; i < getSize(); ++i) {
if (!deleted_[i]) break;
}
return ConstStackIterator(*this, 0);
return ConstStackIterator(*this, i);
}
ConstStackIterator cend() const { return ConstStackIterator(*this, getSize()); }
ConstStackIterator clast() const {
......@@ -215,8 +216,9 @@ namespace corsika::stack {
/**
* increase stack size, create new particle at end of stack
*/
template <typename... Args>
StackIterator AddParticle(const Args... v) {
template <typename... TArgs>
StackIterator AddParticle(const TArgs... v) {
C8LOG_TRACE("Stack::AddParticle");
data_.IncrementSize();
deleted_.push_back(false);
return StackIterator(*this, getSize() - 1, v...);
......@@ -230,8 +232,9 @@ namespace corsika::stack {
* This should only get internally called from a
* StackIterator::AddSecondary via ParticleBase
*/
template <typename... Args>
StackIterator AddSecondary(StackIterator& parent, const Args... v) {
template <typename... TArgs>
StackIterator AddSecondary(StackIterator& parent, const TArgs... v) {
C8LOG_TRACE("Stack::AddSecondary");
data_.IncrementSize();
deleted_.push_back(false);
return StackIterator(*this, getSize() - 1, parent, v...);
......@@ -239,27 +242,43 @@ namespace corsika::stack {
public:
void Swap(StackIterator a, StackIterator b) {
C8LOG_TRACE("Stack::Swap");
data_.Swap(a.GetIndex(), b.GetIndex());
std::swap(deleted_[a.GetIndex()], deleted_[b.GetIndex()]);
}
void Copy(StackIterator a, StackIterator b) {
C8LOG_TRACE("Stack::Copy");
data_.Copy(a.GetIndex(), b.GetIndex());
if (deleted_[b.GetIndex()] && !deleted_[a.GetIndex()]) nDeleted_--;
if (!deleted_[b.GetIndex()] && deleted_[a.GetIndex()]) nDeleted_++;
deleted_[b.GetIndex()] = deleted_[a.GetIndex()];
}
void Copy(ConstStackIterator a, StackIterator b) {
C8LOG_TRACE("Stack::Copy");
data_.Copy(a.GetIndex(), b.GetIndex());
if (deleted_[b.GetIndex()] && !deleted_[a.GetIndex()]) nDeleted_--;
if (!deleted_[b.GetIndex()] && deleted_[a.GetIndex()]) nDeleted_++;
deleted_[b.GetIndex()] = deleted_[a.GetIndex()];
}
std::string as_string() const {
std::string str(fmt::format("size {}, entries {}, deleted {} \n", getSize(), getEntries(), getDeleted()));
// we make our own begin/end since we want ALL entries
std::string new_line = " ";
for (unsigned int iPart = 0; iPart!=getSize(); ++iPart) {
ConstStackIterator itPart(*this, iPart);
str += fmt::format("{}{}{}", new_line, itPart.as_string(), (deleted_[itPart.GetIndex()]?" [deleted]":""));
new_line = "\n ";
}
return str;
}
/**
* delete this particle
*/
public:
void Delete(StackIterator p) {
C8LOG_TRACE("Stack::Delete");
if (IsEmpty()) { /*error*/
throw std::runtime_error("Stack, cannot delete entry since size is zero");
}
......@@ -294,6 +313,7 @@ namespace corsika::stack {
bool purgeLastIfDeleted() {
if (!deleted_.back())
return false; // the last particle is not marked for deletion. Do nothing.
C8LOG_TRACE("Stack::purgeLastIfDeleted: yes");
data_.DecrementSize();
nDeleted_--;
deleted_.pop_back();
......@@ -359,7 +379,10 @@ namespace corsika::stack {
* TStackData data_. By default (and in almost all cases) this
* should just be identiy. See class SecondaryView for an alternative implementation.
*/
unsigned int GetIndexFromIterator(const unsigned int vI) const { return vI; }
unsigned int GetIndexFromIterator(const unsigned int vI) const {
// this is too much: C8LOG_TRACE("Stack::GetIndexFromIterator({})={}", vI, vI);
return vI;
}
/**
* @name Return reference to TStackData object data_ for data access
......
......@@ -86,7 +86,7 @@ namespace corsika::setup {
typename corsika::setup::Stack::StackImpl,
// CHECK with CLANG: corsika::setup::Stack::MPIType>;
// corsika::setup::detail::StackWithGeometryInterface>;
corsika::setup::detail::StackWithHistoryInterface>;
corsika::setup::detail::StackWithHistoryInterface, StackViewProducer>;
#elif defined(__GNUC__) || defined(__GNUG__)
using TheStackView =
corsika::stack::MakeView<corsika::setup::Stack, StackViewProducer>::type;
......
......@@ -9,10 +9,12 @@
#pragma once
#include <corsika/particles/ParticleProperties.h>
#include <corsika/logging/Logging.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
#include <tuple>
#include <string>
namespace corsika::stack {
......@@ -46,6 +48,8 @@ namespace corsika::stack {
void SetParticleData(const std::tuple<NoData>& /*v*/) {}
void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
const std::tuple<NoData>& /*v*/) {}
std::string as_string() const { return "dummy-data"; }
};
/**
......
......@@ -10,6 +10,7 @@
#include <corsika/logging/Logging.h>
#include <corsika/stack/Stack.h>
#include <corsika/logging/Logging.h>
#include <tuple>
#include <utility>
......
......@@ -10,6 +10,7 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/history/SecondaryParticle.hpp>
#include <corsika/logging/Logging.h>
#include <iostream>
#include <memory>
......@@ -20,7 +21,7 @@ namespace corsika::history {
class Event;
using EventPtr = std::shared_ptr<history::Event>;
class Event {
size_t projectileIndex_ = 0; //!< index of projectile on stack
......@@ -28,20 +29,23 @@ namespace corsika::history {
EventPtr parent_event_;
std::optional<corsika::particles::Code>
targetCode_; // cannot be const, value set only after construction
targetCode_;
public:
Event() = default;
void setParentEvent(EventPtr const& evt) { parent_event_ = evt; }
EventPtr parentEvent() { return parent_event_; }
bool hasParentEvent() const { return bool(parent_event_); }
EventPtr& parentEvent() { return parent_event_; }
const EventPtr& parentEvent() const { return parent_event_; }
void setProjectileIndex(size_t i) { projectileIndex_ = i; }
size_t projectileIndex() const { return projectileIndex_; }
template <typename TStackIterator>
TStackIterator projectile(TStackIterator begin) {
// todo: change this
// MR: This is dangerous. You can pass any iterator though it must
// be stack.begin() to yield the correct projectile
......@@ -58,6 +62,10 @@ namespace corsika::history {
std::vector<SecondaryParticle>& secondaries() { return secondaries_; }
void setTargetCode(const particles::Code t) { targetCode_ = t; }
std::string as_string() const {
return fmt::format("hasParent={}, projIndex={}, Nsec={}", hasParentEvent(), projectileIndex_, secondaries_.size());
}
};
} // namespace corsika::history
......@@ -8,6 +8,7 @@
#include <corsika/logging/Logging.h>
#include <corsika/history/HistoryObservationPlane.hpp>
#include <corsika/logging/Logging.h>
#include <boost/histogram/ostream.hpp>
......
......@@ -43,7 +43,8 @@ namespace corsika::history {
// custom data access function
void SetEvent(const int i, EventPtr v) { historyData_[i].first = std::move(v); }
EventPtr GetEvent(const int i) const { return historyData_[i].first; }
const EventPtr& GetEvent(const int i) const { return historyData_[i].first; }
EventPtr& GetEvent(const int i) { return historyData_[i].first; }
void SetParentEventIndex(const int i, ParentEventIndex v) {
historyData_[i].second = std::move(v);
......@@ -108,6 +109,10 @@ namespace corsika::history {
int GetParentEventIndex() const {
return GetStackData().GetParentEventIndex(GetIndex());
}
std::string as_string() const {
return fmt::format("i_parent={}, [evt: {}]", GetParentEventIndex(), (bool(GetEvent())?GetEvent()->as_string():"n/a"));
}
};
template <typename T, typename TEvent>
......
......@@ -14,6 +14,8 @@
#include <corsika/stack/dummy/DummyStack.h>
#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
#include <corsika/logging/Logging.h>
#include <catch2/catch.hpp>
#include <iostream>
......
......@@ -18,6 +18,8 @@
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/logging/Logging.h>
#include <algorithm>
#include <tuple>
#include <vector>
......@@ -60,6 +62,7 @@ namespace corsika::stack {
protected:
using InnerParticleInterface<StackIteratorInterface>::GetStackData;
using InnerParticleInterface<StackIteratorInterface>::GetIndex;
using InnerParticleInterface<StackIteratorInterface>::as_string;
public:
void SetParticleData(
......@@ -141,6 +144,13 @@ namespace corsika::stack {
std::get<4>(v)});
}
std::string as_string() const {
return fmt::format(
"{}, nuc({})", InnerParticleInterface<StackIteratorInterface>::as_string(),
(isNucleus() ? fmt::format("A={}, Z={}", GetNuclearA(), GetNuclearZ())
: "n/a"));
}
/**
* @name individual setters
* @{
......@@ -184,6 +194,7 @@ namespace corsika::stack {
protected:
void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }
bool isNucleus() const { return GetStackData().isNucleus(GetIndex()); }
};
/**
......@@ -240,6 +251,8 @@ namespace corsika::stack {
throw std::runtime_error(err.str());
}
bool isNucleus(const unsigned int i) const { return fNucleusRef[i] >= 0; }
/**
* Function to copy particle at location i1 in stack to i2
*/
......
......@@ -12,6 +12,7 @@ target_link_libraries (
CORSIKAunits
CORSIKAparticles
CORSIKAgeometry
CORSIKAlogging
)
target_include_directories (
......@@ -33,8 +34,6 @@ install (
CORSIKA_ADD_TEST(testSuperStupidStack)
target_link_libraries (
testSuperStupidStack
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
SuperStupidStack
CORSIKAtesting
)
......@@ -16,6 +16,7 @@
#include <corsika/geometry/RootCoordinateSystem.h> // remove
#include <corsika/geometry/Vector.h>
#include <string>
#include <tuple>
#include <vector>
......@@ -40,6 +41,12 @@ namespace corsika::stack {
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
std::string as_string() const {
using namespace corsika::units::si;
return fmt::format("particle: i={}, PID={}, E={}GeV", GetIndex(),
particles::GetName(GetPID()), GetEnergy() / 1_GeV);
}
void SetParticleData(
const std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
MomentumVector, corsika::geometry::Point,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment