From ac66c5ddea5c5fe36f8e8697ed87243fe18e431f Mon Sep 17 00:00:00 2001 From: ralfulrich <ralf.ulrich@kit.edu> Date: Sun, 11 Oct 2020 15:07:38 +0200 Subject: [PATCH] added SwitchProcessSequence, which is the much more general sibling of SwitchProcess (gone) --- .clang-format | 2 +- Documentation/Examples/CMakeLists.txt | 4 - Documentation/Examples/boundary_example.cc | 2 +- Documentation/Examples/cascade_example.cc | 4 +- .../Examples/cascade_proton_example.cc | 2 +- Documentation/Examples/em_shower.cc | 4 +- .../Examples/staticsequence_example.cc | 4 +- Documentation/Examples/vertical_EAS.cc | 48 ++- Environment/CMakeLists.txt | 2 +- Environment/ShowerAxis.cc | 17 +- Framework/Cascade/Cascade.h | 117 +++--- Framework/Cascade/testCascade.cc | 2 +- Framework/ProcessSequence/BaseProcess.h | 1 + Framework/ProcessSequence/CMakeLists.txt | 12 +- Framework/ProcessSequence/ContinuousProcess.h | 1 - Framework/ProcessSequence/DecayProcess.h | 15 +- .../ProcessSequence/InteractionProcess.h | 17 +- Framework/ProcessSequence/ProcessReturn.h | 12 + Framework/ProcessSequence/ProcessSequence.h | 310 ++++++++------- Framework/ProcessSequence/ProcessTraits.h | 33 ++ .../ProcessSequence/SecondariesProcess.h | 4 +- Framework/ProcessSequence/StackProcess.h | 4 +- .../ProcessSequence/SwitchProcessSequence.h | 365 ++++++++++++++++++ .../ProcessSequence/testProcessSequence.cc | 239 ++++++++++-- .../testSwitchProcessSequence.cc | 248 ++++++++++++ .../testSwitchProcessSequence.h | 12 - Framework/Random/CMakeLists.txt | 4 +- Framework/Random/RNGManager.cc | 17 +- Framework/Random/UniformRealDistribution.h | 1 - Framework/StackInterface/SecondaryView.h | 4 +- Framework/StackInterface/Stack.h | 4 +- .../StackInterface/StackIteratorInterface.h | 33 +- Processes/CMakeLists.txt | 1 - Processes/EnergyLoss/EnergyLoss.cc | 7 +- .../LongitudinalProfile.cc | 4 + .../ObservationPlane/ObservationPlane.cc | 21 +- Processes/ObservationPlane/ObservationPlane.h | 2 +- Processes/ParticleCut/ParticleCut.cc | 13 +- Processes/ParticleCut/ParticleCut.h | 7 +- Processes/Proposal/CMakeLists.txt | 1 + Processes/Proposal/ContinuousProcess.cc | 22 +- Processes/Pythia/Decay.h | 2 + Processes/Pythia/Interaction.cc | 2 +- Processes/Pythia/Interaction.h | 2 +- Processes/Sibyll/Decay.cc | 4 +- Processes/Sibyll/Interaction.cc | 326 ++++++++-------- Processes/Sibyll/Interaction.h | 2 +- Processes/StackInspector/StackInspector.cc | 6 +- Processes/StackInspector/StackInspector.h | 4 +- .../StackInspector/testStackInspector.cc | 3 +- Processes/UrQMD/UrQMD.cc | 29 +- 51 files changed, 1426 insertions(+), 576 deletions(-) create mode 100644 Framework/ProcessSequence/ProcessTraits.h create mode 100644 Framework/ProcessSequence/SwitchProcessSequence.h create mode 100644 Framework/ProcessSequence/testSwitchProcessSequence.cc rename Processes/SwitchProcess/testSwitchProcess.cc => Framework/ProcessSequence/testSwitchProcessSequence.h (94%) diff --git a/.clang-format b/.clang-format index 0ac489959..afc401b28 100644 --- a/.clang-format +++ b/.clang-format @@ -57,6 +57,7 @@ IncludeCategories: Priority: 2 - Regex: '.*' Priority: 3 +SortIncludes: false IndentCaseLabels: true IndentWidth: 2 IndentWrappedFunctionNames: false @@ -76,7 +77,6 @@ PenaltyExcessCharacter: 1000000 PenaltyReturnTypeOnItsOwnLine: 200 PointerAlignment: Left ReflowComments: true -SortIncludes: true SpaceAfterCStyleCast: false SpaceBeforeAssignmentOperators: true SpaceBeforeParens: ControlStatements diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 17c711669..aa6c72965 100644 --- a/Documentation/Examples/CMakeLists.txt +++ b/Documentation/Examples/CMakeLists.txt @@ -24,7 +24,6 @@ target_link_libraries (boundary_example ProcessParticleCut ProcessTrackingLine ProcessPythia8 - ProcessProposal CORSIKAprocesses CORSIKAparticles CORSIKAgeometry @@ -66,7 +65,6 @@ if (Pythia8_FOUND) ProcessSibyll ProcessPythia8 ProcessUrQMD - ProcessSwitch CORSIKAcascade ProcessCONEXSourceCut ProcessEnergyLoss @@ -95,7 +93,6 @@ if (Pythia8_FOUND) ProcessSibyll ProcessPythia8 ProcessUrQMD - ProcessSwitch CORSIKAcascade ProcessProposal ProcessPythia8 @@ -146,7 +143,6 @@ target_link_libraries (em_shower ProcessSibyll ProcessPythia8 ProcessUrQMD - ProcessSwitch CORSIKAcascade ProcessCONEXSourceCut ProcessEnergyLoss diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc index 40c60ebed..5d05e2bac 100644 --- a/Documentation/Examples/boundary_example.cc +++ b/Documentation/Examples/boundary_example.cc @@ -127,7 +127,7 @@ int main() { MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat"); // assemble all processes into an ordered process list - auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter; + auto sequence = sibyll % decay % cut % boundaryCrossing % trackWriter; // setup particle stack, and add primary particles setup::Stack stack; diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index fe8e05d87..ff0e419e0 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -147,8 +147,8 @@ int main() { process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()}; // assemble all processes into an ordered process list - auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut - << trackWriter; + auto sequence = stackInspect % sibyll % sibyllNuc % decay % eLoss % cut + % trackWriter; // define air shower object, run simulation cascade::Cascade EAS(env, tracking, sequence, stack); diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc index 6050f7fcf..b7d3970c3 100644 --- a/Documentation/Examples/cascade_proton_example.cc +++ b/Documentation/Examples/cascade_proton_example.cc @@ -136,7 +136,7 @@ int main() { // assemble all processes into an ordered process list // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter; - auto sequence = pythia << decay << cut << trackWriter << stackInspect; + auto sequence = pythia % decay % cut % trackWriter % stackInspect; // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() // << "\n"; diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc index 1c5c07cd4..208deaa93 100644 --- a/Documentation/Examples/em_shower.cc +++ b/Documentation/Examples/em_shower.cc @@ -144,8 +144,8 @@ int main(int argc, char** argv) { process::observation_plane::ObservationPlane observationLevel(obsPlane, "particles.dat"); - auto sequence = proposalCounted << em_continuous << longprof << cut << observationLevel - << trackWriter; + auto sequence = proposalCounted % em_continuous % longprof % cut % observationLevel + % trackWriter; // define air shower object, run simulation tracking_line::TrackingLine tracking; cascade::Cascade EAS(env, tracking, sequence, stack); diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc index ca75fe4b2..13dcd6de1 100644 --- a/Documentation/Examples/staticsequence_example.cc +++ b/Documentation/Examples/staticsequence_example.cc @@ -83,7 +83,7 @@ void modular() { Process3 m3; Process4 m4(0.9); - auto sequence = m1 << m2 << m3 << m4; + auto sequence = m1 % m2 % m3 % m4; DummyData particle; DummyTrajectory track; @@ -93,7 +93,7 @@ void modular() { for (int i = 0; i < nData; ++i) { // cout << p.p[i] << endl; - // assert(p.p[i] == n-i*100); + assert(particle.p[i] == n-i*100); } cout << " done (nothing...) " << endl; diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index 042e64bc2..ce28c92c9 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -22,6 +22,7 @@ #include <corsika/geometry/Sphere.h> #include <corsika/logging/Logging.h> #include <corsika/process/ProcessSequence.h> +#include <corsika/process/SwitchProcessSequence.h> #include <corsika/process/StackProcess.h> #include <corsika/process/energy_loss/EnergyLoss.h> #include <corsika/process/longitudinal_profile/LongitudinalProfile.h> @@ -34,7 +35,6 @@ #include <corsika/process/sibyll/Decay.h> #include <corsika/process/sibyll/Interaction.h> #include <corsika/process/sibyll/NuclearInteraction.h> -#include <corsika/process/switch_process/SwitchProcess.h> #include <corsika/process/tracking_line/TrackingLine.h> #include <corsika/process/urqmd/UrQMD.h> #include <corsika/random/RNGManager.h> @@ -60,6 +60,8 @@ using namespace corsika::environment; using namespace std; using namespace corsika::units::si; +using Particle = setup::Stack::StackIterator; + void registerRandomStreams(const int seed) { random::RNGManager::GetInstance().RegisterRandomStream("cascade"); random::RNGManager::GetInstance().RegisterRandomStream("qgsjet"); @@ -76,7 +78,7 @@ void registerRandomStreams(const int seed) { int main(int argc, char** argv) { - logging::SetLevel(logging::level::info); + logging::SetLevel(logging::level::trace); C8LOG_INFO("vertical_EAS"); @@ -149,23 +151,20 @@ int main(int argc, char** argv) { std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl; if (A != 1) { - stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, - corsika::stack::MomentumVector, geometry::Point, - units::si::TimeType, unsigned short, unsigned short>{ - beamCode, E0, plab, injectionPos, 0_ns, A, Z}); + stack.AddParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z)); } else { stack.AddParticle( - std::tuple<particles::Code, units::si::HEPEnergyType, - corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ - particles::Code::Proton, E0, plab, injectionPos, 0_ns}); + std::make_tuple(particles::Code::Proton, E0, plab, injectionPos, 0_ns)); } - std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02 + // we make the axis much longer than the inj-core distance since the + // profile will go beyond the core, depending on zenith angle + std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.5 << std::endl; environment::ShowerAxis const showerAxis{injectionPos, - (showerCore - injectionPos) * 1.02, env}; + (showerCore - injectionPos) * 1.5, env}; // setup processes, decays and interactions @@ -200,8 +199,7 @@ int main(int argc, char** argv) { decaySibyll.PrintDecayConfig(); - // PROPOSAL processs proposal{...}; - process::particle_cut::ParticleCut cut{60_GeV, false, true}; + process::particle_cut::ParticleCut cut{50_GeV, false, true}; process::proposal::Interaction proposal(env, cut.GetECut()); process::proposal::ContinuousProcess em_continuous(env, cut.GetECut()); process::interaction_counter::InteractionCounter proposalCounted(proposal); @@ -218,14 +216,22 @@ int main(int argc, char** argv) { process::interaction_counter::InteractionCounter urqmdCounted{urqmd}; // assemble all processes into an ordered process list - - auto sibyllSequence = sibyllNucCounted << sibyllCounted; - process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence, - 55_GeV); - auto decaySequence = decayPythia << decaySibyll; - - auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted - << em_continuous << cut << longprof << observationLevel; + struct EnergySwitch { + HEPEnergyType cutE_; + EnergySwitch(HEPEnergyType cutE) + : cutE_(cutE) {} + process::SwitchResult select(const Particle& p) { + if (p.GetEnergy() < cutE_) + return process::SwitchResult::First; + else + return process::SwitchResult::Second; + } + }; + auto hadronSequence = process::select(urqmdCounted, sibyllNucCounted % sibyllCounted, + EnergySwitch(55_GeV)); + auto decaySequence = decayPythia % decaySibyll; + auto sequence = hadronSequence + (reset_particle_mass + decaySequence + proposalCounted) * + em_continuous + cut + observationLevel + longprof; // define air shower object, run simulation tracking_line::TrackingLine tracking; diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt index 5536d2ea3..f0adcf74a 100644 --- a/Environment/CMakeLists.txt +++ b/Environment/CMakeLists.txt @@ -49,7 +49,7 @@ target_link_libraries ( CORSIKAgeometry CORSIKAparticles CORSIKAunits - CORSIKArandom + CORSIKAlogging ) target_include_directories ( diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc index fadc86ab3..ca5a4e45c 100644 --- a/Environment/ShowerAxis.cc +++ b/Environment/ShowerAxis.cc @@ -7,7 +7,9 @@ */ #include <corsika/environment/ShowerAxis.h> -#include <sstream> +#include <corsika/logging/Logging.h> + +#include <string> using namespace corsika::environment; using namespace corsika::units::si; @@ -20,19 +22,20 @@ GrammageType ShowerAxis::X(LengthType l) const { decltype(X_.size()) const upper = lower + 1; if (lower < 0) { + C8LOG_ERROR("cannot extrapolate to points behind point of injection l={} m", l/1_m); throw std::runtime_error("cannot extrapolate to points behind point of injection"); } - if (upper >= X_.size()) { - std::stringstream errormsg; - errormsg << "shower axis too short, cannot extrapolate (l / max_length_ = " - << l / max_length_ << ")"; - throw std::runtime_error(errormsg.str().c_str()); + if (upper >= X_.size()) { + const std::string err = fmt::format("shower axis too short, cannot extrapolate (l / max_length_ = {} )", + l / max_length_); + C8LOG_ERROR(err); + throw std::runtime_error(err.c_str()); } assert(0 <= lambda && lambda <= 1.); - std::cout << l << ": " << lower << " " << lambda << " " << upper << std::endl; + C8LOG_TRACE("ShowerAxis::X l={} m, lower={}, lambda={}, upper={}", l/1_m, lower, lambda, upper); // linear interpolation between X[lower] and X[upper] return X_[upper] * lambda + X_[lower] * (1 - lambda); diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index b520acddd..bf271fcb3 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -74,11 +74,11 @@ namespace corsika::cascade { private: // Data members - corsika::environment::Environment<MediumInterface> const& fEnvironment; - TTracking& fTracking; - TProcessList& fProcessSequence; - TStack& fStack; - corsika::random::RNG& fRNG = + corsika::environment::Environment<MediumInterface> const& environment_; + TTracking& tracking_; + TProcessList& process_sequence_; + TStack& stack_; + corsika::random::RNG& rng_ = corsika::random::RNGManager::GetInstance().GetRandomStream("cascade"); unsigned int count_ = 0; @@ -93,15 +93,15 @@ namespace corsika::cascade { */ Cascade(corsika::environment::Environment<MediumInterface> const& env, TTracking& tr, TProcessList& pl, TStack& stack) - : fEnvironment(env) - , fTracking(tr) - , fProcessSequence(pl) - , fStack(stack) + : environment_(env) + , tracking_(tr) + , process_sequence_(pl) + , stack_(stack) , count_(0) { C8LOG_INFO(c8_ascii_); -#ifdef WITH_HISTORY - C8LOG_INFO(" - With full cascade HISTORY."); -#endif + if constexpr (TStackView::has_event) { + C8LOG_INFO(" - With full cascade HISTORY."); + } } /** @@ -109,9 +109,9 @@ namespace corsika::cascade { * position */ void SetNodes() { - std::for_each(fStack.begin(), fStack.end(), [&](auto& p) { + std::for_each(stack_.begin(), stack_.end(), [&](auto& p) { auto const* numericalNode = - fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition()); + environment_.GetUniverse()->GetContainingNode(p.GetPosition()); p.SetNode(numericalNode); }); } @@ -123,18 +123,18 @@ namespace corsika::cascade { void Run() { SetNodes(); - while (!fStack.IsEmpty()) { - while (!fStack.IsEmpty()) { - C8LOG_TRACE(fmt::format("Stack: {}", fStack.as_string())); + while (!stack_.IsEmpty()) { + while (!stack_.IsEmpty()) { + C8LOG_TRACE("Stack: {}", stack_.as_string()); count_++; - auto pNext = fStack.GetNextParticle(); - C8LOG_DEBUG(fmt::format( + auto pNext = stack_.GetNextParticle(); + C8LOG_DEBUG( "============== next particle : count={}, pid={}, " ", stack entries={}" ", stack deleted={}", - count_, pNext.GetPID(), fStack.getEntries(), fStack.getDeleted())); + count_, pNext.GetPID(), stack_.getEntries(), stack_.getDeleted()); Step(pNext); - fProcessSequence.DoStack(fStack); + process_sequence_.DoStack(stack_); } // do cascade equations, which can put new particles on Stack, // thus, the double loop @@ -149,10 +149,10 @@ namespace corsika::cascade { */ void forceInteraction() { C8LOG_DEBUG("forced interaction!"); - auto vParticle = fStack.GetNextParticle(); + auto vParticle = stack_.GetNextParticle(); TStackView secondaries(vParticle); - interaction(vParticle, secondaries); - fProcessSequence.DoSecondaries(secondaries); + interaction(secondaries); + process_sequence_.DoSecondaries(secondaries); vParticle.Delete(); // todo: this should be reviewed, see below } @@ -172,16 +172,16 @@ namespace corsika::cascade { using namespace corsika::units::si; // determine geometric tracking - auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle); + auto [step, geomMaxLength, nextVol] = tracking_.GetTrack(vParticle); [[maybe_unused]] auto const& dummy_nextVol = nextVol; // determine combined total interaction length (inverse) InverseGrammageType const total_inv_lambda = - fProcessSequence.GetTotalInverseInteractionLength(vParticle); + process_sequence_.GetInverseInteractionLength(vParticle); // sample random exponential step length in grammage corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda); - GrammageType const next_interact = expDist(fRNG); + GrammageType const next_interact = expDist(rng_); C8LOG_DEBUG( "total_lambda={} g/cm2, " @@ -193,8 +193,8 @@ namespace corsika::cascade { // assert that particle stays outside void Universe if it has no // model properties set - assert(currentLogicalNode != &*fEnvironment.GetUniverse() || - fEnvironment.GetUniverse()->HasModelProperties()); + assert(currentLogicalNode != &*environment_.GetUniverse() || + environment_.GetUniverse()->HasModelProperties()); // convert next_step from grammage to length LengthType const distance_interact = @@ -202,16 +202,16 @@ namespace corsika::cascade { next_interact); // determine the maximum geometric step length from continuous processes - LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step); + LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step); C8LOG_DEBUG("distance_max={} m", distance_max / 1_m); // determine combined total inverse decay time InverseTimeType const total_inv_lifetime = - fProcessSequence.GetTotalInverseLifetime(vParticle); + process_sequence_.GetInverseLifetime(vParticle); // sample random exponential decay time corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime); - TimeType const next_decay = expDistDecay(fRNG); + TimeType const next_decay = expDistDecay(rng_); C8LOG_DEBUG( "total_lifetime={} s" ", next_decay={} s", @@ -236,12 +236,12 @@ namespace corsika::cascade { step.LimitEndTo(min_distance); // apply all continuous processes on particle + track - process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step); - - if (status == process::EProcessReturn::eParticleAbsorbed) { + if (process_sequence_.DoContinuous(vParticle, step) == + process::EProcessReturn::eParticleAbsorbed) { C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV", vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV); - vParticle.Delete(); + if (!vParticle.isDeleted()) + vParticle.Delete(); return; } @@ -271,10 +271,10 @@ namespace corsika::cascade { [[maybe_unused]] auto projectile = secondaries.GetProjectile(); if (min_distance == distance_interact) { - interaction(vParticle, secondaries); + interaction(secondaries); } else { assert(min_distance == distance_decay); - decay(vParticle, secondaries); + decay(secondaries); // make sure particle actually did decay if it should have done so if (secondaries.getSize() == 1 && projectile.GetPID() == secondaries.GetNextParticle().GetPID()) @@ -283,7 +283,7 @@ namespace corsika::cascade { particles::GetName(projectile.GetPID()))); } - fProcessSequence.DoSecondaries(secondaries); + process_sequence_.DoSecondaries(secondaries); vParticle.Delete(); } else { // step-length limitation within volume @@ -293,10 +293,9 @@ namespace corsika::cascade { [[maybe_unused]] auto const assertion = [&] { auto const* numericalNodeAfterStep = - fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition()); - C8LOG_TRACE(fmt::format( - "Geometry check: numericalNodeAfterStep={} currentLogicalNode={}", - fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode))); + environment_.GetUniverse()->GetContainingNode(vParticle.GetPosition()); + C8LOG_TRACE("Geometry check: numericalNodeAfterStep={} currentLogicalNode={}", + fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode)); return numericalNodeAfterStep == currentLogicalNode; }; @@ -306,41 +305,45 @@ namespace corsika::cascade { /* DoBoundary may delete the particle (or not) - small caveat: any changes to vParticle, or even the production + caveat: any changes to vParticle, or even the production of new secondaries is currently not passed to ParticleCut, thus, particles outside the desired phase space may be produced. + + todo: this must be fixed. */ - fProcessSequence.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol); + process_sequence_.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol); } } - auto decay(Particle& particle, TStackView& view) { + process::EProcessReturn decay(TStackView& view) { C8LOG_DEBUG("decay"); units::si::InverseTimeType const actual_decay_time = - fProcessSequence.GetTotalInverseLifetime(particle); + process_sequence_.GetInverseLifetime(view.parent()); random::UniformRealDistribution<units::si::InverseTimeType> uniDist( actual_decay_time); - const auto sample_process = uniDist(fRNG); - units::si::InverseTimeType inv_decay_count = units::si::InverseTimeType::zero(); - auto const returnCode = - fProcessSequence.SelectDecay(particle, view, sample_process, inv_decay_count); + const auto sample_process = uniDist(rng_); + auto const returnCode = process_sequence_.SelectDecay(view, sample_process); + if (returnCode != process::EProcessReturn::eDecayed) { + C8LOG_WARN("Particle did not decay!"); + } SetEventType(view, history::EventType::Decay); return returnCode; } - auto interaction(Particle& particle, TStackView& view) { + process::EProcessReturn interaction(TStackView& view) { C8LOG_DEBUG("collide"); units::si::InverseGrammageType const current_inv_length = - fProcessSequence.GetTotalInverseInteractionLength(particle); + process_sequence_.GetInverseInteractionLength(view.parent()); random::UniformRealDistribution<units::si::InverseGrammageType> uniDist( current_inv_length); - const auto sample_process = uniDist(fRNG); - auto inv_lambda_count = units::si::InverseGrammageType::zero(); - auto const returnCode = fProcessSequence.SelectInteraction( - particle, view, sample_process, inv_lambda_count); + const auto sample_process = uniDist(rng_); + auto const returnCode = process_sequence_.SelectInteraction(view, sample_process); + if (returnCode != process::EProcessReturn::eInteracted) { + C8LOG_WARN("Particle did not interace!"); + } SetEventType(view, history::EventType::Interaction); return returnCode; } diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 02f974f52..0c969d3d2 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -138,7 +138,7 @@ TEST_CASE("Cascade", "[Cascade]") { const HEPEnergyType Ecrit = 85_MeV; ProcessSplit split(X0); ProcessCut cut(Ecrit); - auto sequence = nullModel << stackInspect << split << cut; + auto sequence = nullModel % stackInspect % split % cut; TestCascadeStack stack; stack.Clear(); stack.AddParticle( diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h index 42934bbe9..c1a4e232e 100644 --- a/Framework/ProcessSequence/BaseProcess.h +++ b/Framework/ProcessSequence/BaseProcess.h @@ -21,6 +21,7 @@ namespace corsika::process { ProcessSequence. Both, the ProcessSequence and all its elements are of type BaseProcess<T> + \todo rename BaseProcess into just Process */ class _BaseProcess {}; diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt index 811e2c386..2c4a0fe5d 100644 --- a/Framework/ProcessSequence/CMakeLists.txt +++ b/Framework/ProcessSequence/CMakeLists.txt @@ -17,7 +17,9 @@ set ( StackProcess.h DecayProcess.h ProcessSequence.h + SwitchProcessSequence.h ProcessReturn.h + ProcessTraits.h ) CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAprocesssequence ${CORSIKAprocesssequence_NAMESPACE} ${CORSIKAprocesssequence_HEADERS}) @@ -27,6 +29,7 @@ target_link_libraries ( CORSIKAprocesssequence INTERFACE CORSIKAsetup + CORSIKAlogging ) target_include_directories ( @@ -53,8 +56,15 @@ target_link_libraries ( CORSIKA_ADD_TEST (testProcessSequence) target_link_libraries ( testProcessSequence - ProcessSwitch CORSIKAgeometry CORSIKAprocesssequence CORSIKAtesting ) + +# # CORSIKA_ADD_TEST (testSwitchProcessSequence) +# # target_link_libraries ( +# # testSwitchProcessSequence +# # CORSIKAgeometry +# # CORSIKAprocesssequence +# # CORSIKAtesting +# # ) diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h index 2c9e7bbc8..c4138a755 100644 --- a/Framework/ProcessSequence/ContinuousProcess.h +++ b/Framework/ProcessSequence/ContinuousProcess.h @@ -9,7 +9,6 @@ #pragma once #include <corsika/process/BaseProcess.h> -#include <corsika/process/ProcessReturn.h> // for convenience #include <corsika/units/PhysicalUnits.h> namespace corsika::process { diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h index eb0cf982f..d047d214e 100644 --- a/Framework/ProcessSequence/DecayProcess.h +++ b/Framework/ProcessSequence/DecayProcess.h @@ -9,8 +9,6 @@ #pragma once #include <corsika/process/BaseProcess.h> -#include <corsika/process/ProcessReturn.h> // for convenience -#include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> namespace corsika::process { @@ -35,12 +33,19 @@ namespace corsika::process { EProcessReturn DoDecay(TParticle&); template <typename TParticle> - corsika::units::si::TimeType GetLifetime(TParticle& p); + corsika::units::si::TimeType GetLifetime(const TParticle&); template <typename TParticle> - corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& vP) { - return 1. / GetRef().GetLifetime(vP); + corsika::units::si::InverseTimeType GetInverseLifetime(const TParticle& particle) { + return 1. / GetRef().GetLifetime(particle); } + + /* template <typename TParticle> + corsika::units::si::InverseTimeType GetInverseInteractionLength(TParticle&& particle) { + auto p = std::move(particle); + return 1. / GetRef().GetLifetime(p); + }*/ + }; } // namespace corsika::process diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h index cfecccc81..487668f5b 100644 --- a/Framework/ProcessSequence/InteractionProcess.h +++ b/Framework/ProcessSequence/InteractionProcess.h @@ -9,8 +9,6 @@ #pragma once #include <corsika/process/BaseProcess.h> -#include <corsika/process/ProcessReturn.h> // for convenience -#include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> namespace corsika::process { @@ -35,12 +33,19 @@ namespace corsika::process { EProcessReturn DoInteraction(TParticle&); template <typename TParticle> - corsika::units::si::GrammageType GetInteractionLength(TParticle& p); + corsika::units::si::GrammageType GetInteractionLength(const TParticle&); template <typename TParticle> - corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) { - return 1. / GetRef().GetInteractionLength(p); + corsika::units::si::InverseGrammageType GetInverseInteractionLength(const TParticle& particle) { + return 1. / GetRef().GetInteractionLength(particle); } - }; + + /* + template <typename TParticle> + corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle&& particle) { + auto p = std::move(particle); + return 1. / GetRef().GetInteractionLength(p); + }*/ +}; } // namespace corsika::process diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h index ae1e45b95..c4e7d5890 100644 --- a/Framework/ProcessSequence/ProcessReturn.h +++ b/Framework/ProcessSequence/ProcessReturn.h @@ -39,8 +39,20 @@ namespace corsika::process { return (static_cast<int>(a) & static_cast<int>(b)) != 0; } + inline bool isOk(const EProcessReturn a) { + return static_cast<int>(a & EProcessReturn::eOk); + } + inline bool isAbsorbed(const EProcessReturn a) { return static_cast<int>(a & EProcessReturn::eParticleAbsorbed); } + inline bool isDecayed(const EProcessReturn a) { + return static_cast<int>(a & EProcessReturn::eDecayed); + } + + inline bool isInteracted(const EProcessReturn a) { + return static_cast<int>(a & EProcessReturn::eInteracted); + } + } // namespace corsika::process diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index ea1254613..a959d40fd 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -9,6 +9,7 @@ #pragma once #include <corsika/process/BaseProcess.h> +#include <corsika/process/ProcessTraits.h> #include <corsika/process/BoundaryCrossingProcess.h> #include <corsika/process/ContinuousProcess.h> #include <corsika/process/DecayProcess.h> @@ -17,6 +18,7 @@ #include <corsika/process/SecondariesProcess.h> #include <corsika/process/StackProcess.h> #include <corsika/units/PhysicalUnits.h> +#include <corsika/logging/Logging.h> #include <cmath> #include <limits> @@ -30,99 +32,79 @@ namespace corsika::process { generate a new type based on template logic containing all the elements. + TProcess1 and TProcess2 must both be derived from BaseProcess, + and are both references if possible (lvalue), otherwise (rvalue) + they are just classes. This allows us to handle both, rvalue as + well as lvalue Processes in the ProcessSequence. + \comment Using CRTP pattern, https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern */ + template <typename TProcess1, typename TProcess2> + class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> { - // this is a marker to track which BaseProcess is also a ProcessSequence - template <typename TClass> - struct is_process_sequence : std::false_type {}; - - template <typename TClass> - bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value; - - // we also need a marker to identify SwitchProcess - namespace switch_process { - template <typename TProcess1, typename TProcess2> - class SwitchProcess; // fwd-decl. - } - - // to detect SwitchProcesses inside the ProcessSequence - template <typename TClass> - struct is_switch_process : std::false_type {}; - - template <typename TClass> - bool constexpr is_switch_process_v = is_switch_process<TClass>::value; - - template <typename Process1, typename Process2> - struct is_switch_process<switch_process::SwitchProcess<Process1, Process2>> - : std::true_type {}; - - /** - T1 and T2 are both references if possible (lvalue), otherwise - (rvalue) they are just classes. This allows us to handle both, - rvalue as well as lvalue Processes in the ProcessSequence. - */ - template <typename T1, typename T2> - class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2>> { + using TProcess1type = typename std::decay<TProcess1>::type; + using TProcess2type = typename std::decay<TProcess2>::type; - using T1type = typename std::decay<T1>::type; - using T2type = typename std::decay<T2>::type; + static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>; + static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>; - static bool constexpr t1ProcSeq = is_process_sequence_v<T1type>; - static bool constexpr t2ProcSeq = is_process_sequence_v<T2type>; + static bool constexpr t1SwitchProcSeq = is_switch_process_sequence_v<TProcess1type>; + static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<TProcess2type>; - static bool constexpr t1SwitchProc = is_switch_process_v<T1type>; - static bool constexpr t2SwitchProc = is_switch_process_v<T2type>; + protected: + TProcess1 A; // this is a reference, if possible + TProcess2 B; // this is a reference, if possible public: - T1 A; // this is a reference, if possible - T2 B; // this is a reference, if possible - - ProcessSequence(T1 in_A, T2 in_B) + ProcessSequence(TProcess1 in_A, TProcess2 in_B) : A(in_A) , B(in_B) {} template <typename Particle, typename VTNType> - EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from, + EProcessReturn DoBoundaryCrossing(Particle& particle, VTNType const& from, VTNType const& to) { EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T1type>, T1type> || + if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess1type>, + TProcess1type> || t1ProcSeq) { - ret |= A.DoBoundaryCrossing(p, from, to); + ret |= A.DoBoundaryCrossing(particle, from, to); } - if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T2type>, T2type> || + if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess2type>, + TProcess2type> || t2ProcSeq) { - ret |= B.DoBoundaryCrossing(p, from, to); + ret |= B.DoBoundaryCrossing(particle, from, to); } return ret; } template <typename TParticle, typename TTrack> - EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) { + inline EProcessReturn DoContinuous(TParticle& particle, TTrack& vT) { EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) { - if (!isAbsorbed(ret)) { ret |= A.DoContinuous(vP, vT); } + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, TProcess1type> || + t1ProcSeq) { + ret |= A.DoContinuous(particle, vT); } - if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) { - if (!isAbsorbed(ret)) { ret |= B.DoContinuous(vP, vT); } + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, TProcess2type> || + t2ProcSeq) { + if (!isAbsorbed(ret)) { ret |= B.DoContinuous(particle, vT); } } return ret; } template <typename TSecondaries> - EProcessReturn DoSecondaries(TSecondaries& vS) { - EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of_v<SecondariesProcess<T1type>, T1type> || t1ProcSeq) { - ret |= A.DoSecondaries(vS); + inline void DoSecondaries(TSecondaries& vS) { + if constexpr (std::is_base_of_v<SecondariesProcess<TProcess1type>, TProcess1type> || + t1ProcSeq) { + A.DoSecondaries(vS); } - if constexpr (std::is_base_of_v<SecondariesProcess<T2type>, T2type> || t2ProcSeq) { - ret |= B.DoSecondaries(vS); + if constexpr (std::is_base_of_v<SecondariesProcess<TProcess2type>, TProcess2type> || + t2ProcSeq) { + B.DoSecondaries(vS); } - return ret; } /** @@ -133,12 +115,14 @@ namespace corsika::process { tested if either A or B are StackProcess and if they are due for execution. */ - bool CheckStep() { + inline bool CheckStep() { bool ret = false; - if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) { + if constexpr (std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type> || + (t1ProcSeq && !t1SwitchProcSeq)) { ret |= A.CheckStep(); } - if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) { + if constexpr (std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type> || + (t2ProcSeq && !t2SwitchProcSeq)) { ret |= B.CheckStep(); } return ret; @@ -148,160 +132,161 @@ namespace corsika::process { Execute the StackProcess-es in the ProcessSequence */ template <typename TStack> - EProcessReturn DoStack(TStack& vS) { - EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) { - if (A.CheckStep()) { ret |= A.DoStack(vS); } + inline void DoStack(TStack& stack) { + if constexpr (std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type> || + (t1ProcSeq && !t1SwitchProcSeq)) { + if (A.CheckStep()) { A.DoStack(stack); } } - if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) { - if (B.CheckStep()) { ret |= B.DoStack(vS); } + if constexpr (std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type> || + (t2ProcSeq && !t2SwitchProcSeq)) { + if (B.CheckStep()) { B.DoStack(stack); } } - return ret; } template <typename TParticle, typename TTrack> - corsika::units::si::LengthType MaxStepLength(TParticle& vP, TTrack& vTrack) { + inline corsika::units::si::LengthType MaxStepLength(TParticle& particle, TTrack& vTrack) { corsika::units::si::LengthType max_length = // if no other process in the sequence implements it std::numeric_limits<double>::infinity() * corsika::units::si::meter; - if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) { - corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack); + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, TProcess1type> || + t1ProcSeq) { + corsika::units::si::LengthType const len = A.MaxStepLength(particle, vTrack); max_length = std::min(max_length, len); } - if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) { - corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack); + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, TProcess2type> || + t2ProcSeq) { + corsika::units::si::LengthType const len = B.MaxStepLength(particle, vTrack); max_length = std::min(max_length, len); } return max_length; } template <typename TParticle> - corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP) { - return 1. / GetInverseInteractionLength(vP); + inline corsika::units::si::GrammageType GetInteractionLength(TParticle&& particle) { + return 1. / GetInverseInteractionLength(particle); } template <typename TParticle> - corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength( - TParticle& vP) { - return GetInverseInteractionLength(vP); - } - - template <typename TParticle> - corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP) { + inline corsika::units::si::InverseGrammageType GetInverseInteractionLength( + TParticle&& particle) { using namespace corsika::units::si; - InverseGrammageType tot = 0 * meter * meter / gram; + InverseGrammageType tot = 0 * meter * meter / gram; // default value - if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type> || t1ProcSeq || - t1SwitchProc) { - tot += A.GetInverseInteractionLength(vP); + if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>, TProcess1type> || + t1ProcSeq) { + tot += A.GetInverseInteractionLength(particle); } - if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type> || t2ProcSeq || - t2SwitchProc) { - tot += B.GetInverseInteractionLength(vP); + if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>, TProcess2type> || + t2ProcSeq) { + tot += B.GetInverseInteractionLength(particle); } return tot; } - template <typename TParticle, typename TSecondaries> - EProcessReturn SelectInteraction( - TParticle& vP, TSecondaries& vS, - [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select, - corsika::units::si::InverseGrammageType& lambda_inv_count) { + template <typename TSecondaryView> + inline EProcessReturn SelectInteraction( + TSecondaryView& view, + [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_select, + [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_sum = + corsika::units::si::InverseGrammageType::zero()) { - if constexpr (t1ProcSeq || t1SwitchProc) { + // TODO: add check for lambda_inv_select>lambda_inv_tot + + if constexpr (t1ProcSeq) { // if A is a process sequence --> check inside const EProcessReturn ret = - A.SelectInteraction(vP, vS, lambda_select, lambda_inv_count); - // if A did succeed, stop routine + A.SelectInteraction(view, lambda_inv_select, lambda_inv_sum); + // if A did succeed, stop routine. Not checking other static branch B. if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type>) { + } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>, + TProcess1type>) { // if this is not a ContinuousProcess --> evaluate probability - lambda_inv_count += A.GetInverseInteractionLength(vP); + const auto particle = view.parent(); + lambda_inv_sum += A.GetInverseInteractionLength(particle); // check if we should execute THIS process and then EXIT - if (lambda_select < lambda_inv_count) { - A.DoInteraction(vS); + if (lambda_inv_select < lambda_inv_sum) { + A.DoInteraction(view); return EProcessReturn::eInteracted; } } // end branch A - if constexpr (t2ProcSeq || t2SwitchProc) { - // if A is a process sequence --> check inside - const EProcessReturn ret = - B.SelectInteraction(vP, vS, lambda_select, lambda_inv_count); - // if A did succeed, stop routine - if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type>) { + if constexpr (t2ProcSeq) { + // if B is a process sequence --> check inside + return B.SelectInteraction(view, lambda_inv_select, lambda_inv_sum); + } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>, + TProcess2type>) { // if this is not a ContinuousProcess --> evaluate probability - lambda_inv_count += B.GetInverseInteractionLength(vP); + lambda_inv_sum += B.GetInverseInteractionLength(view.parent()); // check if we should execute THIS process and then EXIT - if (lambda_select < lambda_inv_count) { - B.DoInteraction(vS); + if (lambda_inv_select < lambda_inv_sum) { + B.DoInteraction(view); return EProcessReturn::eInteracted; } - } // end branch A + } // end branch B return EProcessReturn::eOk; } template <typename TParticle> - corsika::units::si::TimeType GetTotalLifetime(TParticle& p) { - return 1. / GetInverseLifetime(p); - } - - template <typename TParticle> - corsika::units::si::InverseTimeType GetTotalInverseLifetime(TParticle& p) { - return GetInverseLifetime(p); + inline corsika::units::si::TimeType GetLifetime(TParticle& particle) { + return 1. / GetInverseLifetime(particle); } template <typename TParticle> - corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& p) { + inline corsika::units::si::InverseTimeType GetInverseLifetime(TParticle&& particle) { using namespace corsika::units::si; - corsika::units::si::InverseTimeType tot = 0 / second; + corsika::units::si::InverseTimeType tot = 0 / second; // default value - if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type> || t1ProcSeq) { - tot += A.GetInverseLifetime(p); + if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, TProcess1type> || + t1ProcSeq) { + tot += A.GetInverseLifetime(particle); } - if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type> || t2ProcSeq) { - tot += B.GetInverseLifetime(p); + if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, TProcess2type> || + t2ProcSeq) { + tot += B.GetInverseLifetime(particle); } return tot; } // select decay process - template <typename TParticle, typename TSecondaries> - EProcessReturn SelectDecay( - TParticle& vP, TSecondaries& vS, - [[maybe_unused]] corsika::units::si::InverseTimeType decay_select, - corsika::units::si::InverseTimeType& decay_inv_count) { + template <typename TSecondaryView> + inline EProcessReturn SelectDecay( + TSecondaryView& view, + [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_select, + [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_sum = + corsika::units::si::InverseTimeType::zero()) { + + // TODO: add check for decay_inv_select>decay_inv_tot + if constexpr (t1ProcSeq) { // if A is a process sequence --> check inside - const EProcessReturn ret = A.SelectDecay(vP, vS, decay_select, decay_inv_count); - // if A did succeed, stop routine + const EProcessReturn ret = A.SelectDecay(view, decay_inv_select, decay_inv_sum); + // if A did succeed, stop routine here (not checking other static branch B) if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type>) { + } else if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, + TProcess1type>) { // if this is not a ContinuousProcess --> evaluate probability - decay_inv_count += A.GetInverseLifetime(vP); + decay_inv_sum += A.GetInverseLifetime(view.parent()); // check if we should execute THIS process and then EXIT - if (decay_select < decay_inv_count) { // more pedagogical: rndm_select < - // decay_inv_count / decay_inv_tot - A.DoDecay(vS); + if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select < + // decay_inv_sum / decay_inv_tot + A.DoDecay(view); return EProcessReturn::eDecayed; } } // end branch A if constexpr (t2ProcSeq) { - // if A is a process sequence --> check inside - const EProcessReturn ret = B.SelectDecay(vP, vS, decay_select, decay_inv_count); - // if A did succeed, stop routine - if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type>) { + // if B is a process sequence --> check inside + return B.SelectDecay(view, decay_inv_select, decay_inv_sum); + } else if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, + TProcess2type>) { // if this is not a ContinuousProcess --> evaluate probability - decay_inv_count += B.GetInverseLifetime(vP); + decay_inv_sum += B.GetInverseLifetime(view.parent()); // check if we should execute THIS process and then EXIT - if (decay_select < decay_inv_count) { - B.DoDecay(vS); + if (decay_inv_select < decay_inv_sum) { + B.DoDecay(view); return EProcessReturn::eDecayed; } } // end branch B @@ -309,13 +294,14 @@ namespace corsika::process { } }; - // the << operator assembles many BaseProcess, ContinuousProcess, and + // the % operator assembles many BaseProcess, ContinuousProcess, and // Interaction/DecayProcess objects into a ProcessSequence, all combinatorics // must be allowed, this is why we define a macro to define all // combinations here: - // enable the << operator to construct ProcessSequence from two - // Processes, only if poth Processes derive from BaseProcesses + // enable the % operator to construct ProcessSequence from two + // Processes, only if both, Processes1 and Processes2, derive from + // BaseProcesses template <typename TProcess1, typename TProcess2> inline typename std::enable_if< @@ -324,11 +310,35 @@ namespace corsika::process { std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>, typename std::decay<TProcess2>::type>::value, ProcessSequence<TProcess1, TProcess2>>::type - operator<<(TProcess1&& vA, TProcess2&& vB) { + operator%(TProcess1&& vA, TProcess2&& vB) { + return ProcessSequence<TProcess1, TProcess2>(vA, vB); + } + + template <typename TProcess1, typename TProcess2> + inline typename std::enable_if< + std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>, + typename std::decay<TProcess1>::type>::value && + std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>, + typename std::decay<TProcess2>::type>::value, + ProcessSequence<TProcess1, TProcess2>>::type + operator+(TProcess1&& vA, TProcess2&& vB) { + C8LOG_TRACE("ProcessSequence (+) for\n\t\t {} and\n\t\t {}", typeid(vA).name(), typeid(vB).name()); + return ProcessSequence<TProcess1, TProcess2>(vA, vB); + } + + template <typename TProcess1, typename TProcess2> + inline typename std::enable_if< + std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>, + typename std::decay<TProcess1>::type>::value && + std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>, + typename std::decay<TProcess2>::type>::value, + ProcessSequence<TProcess1, TProcess2>>::type + operator*(TProcess1&& vA, TProcess2&& vB) { + C8LOG_TRACE("ProcessSequence (*) for\n\t\t {} and\n\t\t {}", typeid(vA).name(), typeid(vB).name()); return ProcessSequence<TProcess1, TProcess2>(vA, vB); } - /// marker to identify objectas ProcessSequence + /// traits marker to identify objectas ProcessSequence template <typename TProcess1, typename TProcess2> struct is_process_sequence<corsika::process::ProcessSequence<TProcess1, TProcess2>> : std::true_type {}; diff --git a/Framework/ProcessSequence/ProcessTraits.h b/Framework/ProcessSequence/ProcessTraits.h new file mode 100644 index 000000000..5c272d618 --- /dev/null +++ b/Framework/ProcessSequence/ProcessTraits.h @@ -0,0 +1,33 @@ +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +namespace corsika::process { + + /** + * A traits marker to track which BaseProcess is also a ProcessSequence + **/ + template <typename TClass> + struct is_process_sequence : std::false_type {}; + + template <typename TClass> + bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value; + + + /** + * A traits marker to identiy a BaseProcess that is also SwitchProcessesSequence + **/ + + template <typename TClass> + struct is_switch_process_sequence : std::false_type {}; + + template <typename TClass> + bool constexpr is_switch_process_sequence_v = is_switch_process_sequence<TClass>::value; + +} // namespace corsika::process diff --git a/Framework/ProcessSequence/SecondariesProcess.h b/Framework/ProcessSequence/SecondariesProcess.h index 355a465c7..a9bc3832a 100644 --- a/Framework/ProcessSequence/SecondariesProcess.h +++ b/Framework/ProcessSequence/SecondariesProcess.h @@ -9,8 +9,6 @@ #pragma once #include <corsika/process/BaseProcess.h> -#include <corsika/process/ProcessReturn.h> // for convenience -#include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> namespace corsika::process { @@ -30,7 +28,7 @@ namespace corsika::process { /// here starts the interface-definition part // -> enforce TDerived to implement DoSecondaries... template <typename TSecondaries> - inline EProcessReturn DoSecondaries(TSecondaries&); + inline void DoSecondaries(TSecondaries&); }; } // namespace corsika::process diff --git a/Framework/ProcessSequence/StackProcess.h b/Framework/ProcessSequence/StackProcess.h index 64a5dff97..ce5c0998d 100644 --- a/Framework/ProcessSequence/StackProcess.h +++ b/Framework/ProcessSequence/StackProcess.h @@ -9,8 +9,6 @@ #pragma once #include <corsika/process/BaseProcess.h> -#include <corsika/process/ProcessReturn.h> // for convenience -#include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> namespace corsika::process { @@ -36,7 +34,7 @@ namespace corsika::process { /// here starts the interface-definition part // -> enforce TDerived to implement DoStack... template <typename TStack> - inline EProcessReturn DoStack(TStack&); + inline void DoStack(TStack&); int GetStep() const { return fIStep; } bool CheckStep() { return !((++fIStep) % fNStep); } diff --git a/Framework/ProcessSequence/SwitchProcessSequence.h b/Framework/ProcessSequence/SwitchProcessSequence.h new file mode 100644 index 000000000..20fb1b0db --- /dev/null +++ b/Framework/ProcessSequence/SwitchProcessSequence.h @@ -0,0 +1,365 @@ +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/process/BaseProcess.h> +#include <corsika/process/ProcessTraits.h> +#include <corsika/process/BoundaryCrossingProcess.h> +#include <corsika/process/ContinuousProcess.h> +#include <corsika/process/DecayProcess.h> +#include <corsika/process/InteractionProcess.h> +#include <corsika/process/ProcessReturn.h> +#include <corsika/process/SecondariesProcess.h> +#include <corsika/process/StackProcess.h> +#include <corsika/units/PhysicalUnits.h> + +#include <cmath> +#include <limits> + +namespace corsika::process { + + // enum for the process switch selection: identify if First or + // Second process branch should be used. + enum class SwitchResult { First, Second }; + + /** + \class SwitchProcessSequence + + A compile time static list of processes that uses an internal + TSelect class to switch between different versions of processes + (or process sequence). + + TProcess1 and TProcess2 must be derived from BaseProcess and are + both references if possible (lvalue), otherwise (rvalue) they are + just classes. This allows us to handle both, rvalue as well as + lvalue Processes in the SwitchProcessSequence. + + TSelect has to implement a `select(const Particle&)` and has to + return either SwitchResult::First or SwitchResult::Second. Note: + TSelect may absolutely also use random numbers to sample between + its results. This can be used to achieve arbitrarily smooth + transition or mixtures of processes. + + Warning: do not put StackProcess into a SwitchProcessSequence + since this makes no sense. The StackProcess acts on an entire + particle stack and not on indiviidual particles. + + \comment See also class ProcessSequence + **/ + + template <typename TProcess1, typename TProcess2, typename TSelect> + class SwitchProcessSequence + : public BaseProcess<SwitchProcessSequence<TProcess1, TProcess2, TSelect>> { + + using TProcess1type = typename std::decay<TProcess1>::type; + using TProcess2type = typename std::decay<TProcess2>::type; + + static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>; + static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>; + + TSelect select_; // this is a reference, if possible + + public: + TProcess1 A; // this is a reference, if possible + TProcess2 B; // this is a reference, if possible + + SwitchProcessSequence(TProcess1 in_A, TProcess2 in_B, TSelect sel) + : select_(sel) + , A(in_A) + , B(in_B) {} + + template <typename Particle, typename VTNType> + EProcessReturn DoBoundaryCrossing(Particle& particle, VTNType const& from, + VTNType const& to) { + + switch (select_.select(particle)) { + case SwitchResult::First: { + if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess1type>, + TProcess1type> || + t1ProcSeq) { + return A.DoBoundaryCrossing(particle, from, to); + } + break; + } + case SwitchResult::Second: { + if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess2type>, + TProcess2type> || + t2ProcSeq) { + return B.DoBoundaryCrossing(particle, from, to); + } + break; + } + } + return EProcessReturn::eOk; + } + + template <typename TParticle, typename TTrack> + inline EProcessReturn DoContinuous(TParticle& particle, TTrack& vT) { + switch (select_.select(particle)) { + case SwitchResult::First: { + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, + TProcess1type> || + t1ProcSeq) { + return A.DoContinuous(particle, vT); + } + break; + } + case SwitchResult::Second: { + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, + TProcess2type> || + t2ProcSeq) { + return B.DoContinuous(particle, vT); + } + break; + } + } + return EProcessReturn::eOk; + } + + template <typename TSecondaries> + inline void DoSecondaries(TSecondaries& vS) { + const auto& particle = vS.parent(); + switch (select_.select(particle)) { + case SwitchResult::First: { + if constexpr (std::is_base_of_v<SecondariesProcess<TProcess1type>, + TProcess1type> || + t1ProcSeq) { + A.DoSecondaries(vS); + } + break; + } + case SwitchResult::Second: { + if constexpr (std::is_base_of_v<SecondariesProcess<TProcess2type>, + TProcess2type> || + t2ProcSeq) { + B.DoSecondaries(vS); + } + break; + } + } + } + + template <typename TParticle, typename TTrack> + inline corsika::units::si::LengthType MaxStepLength(TParticle& particle, + TTrack& vTrack) { + + switch (select_.select(particle)) { + case SwitchResult::First: { + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, + TProcess1type> || + t1ProcSeq) { + return A.MaxStepLength(particle, vTrack); + } + break; + } + case SwitchResult::Second: { + if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, + TProcess2type> || + t2ProcSeq) { + return B.MaxStepLength(particle, vTrack); + } + break; + } + } + + // if no other process in the sequence implements it + return std::numeric_limits<double>::infinity() * corsika::units::si::meter; + } + + template <typename TParticle> + inline corsika::units::si::GrammageType GetInteractionLength(TParticle&& particle) { + return 1. / GetInverseInteractionLength(particle); + } + + template <typename TParticle> + inline corsika::units::si::InverseGrammageType GetInverseInteractionLength( + TParticle&& particle) { + using namespace corsika::units::si; + + switch (select_.select(particle)) { + case SwitchResult::First: { + if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>, + TProcess1type> || + t1ProcSeq) { + return A.GetInverseInteractionLength(particle); + } + break; + } + case SwitchResult::Second: { + if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>, + TProcess2type> || + t2ProcSeq) { + return B.GetInverseInteractionLength(particle); + } + break; + } + } + return 0 * meter * meter / gram; // default value + } + + template <typename TSecondaryView> + inline EProcessReturn SelectInteraction( + TSecondaryView& view, + [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_select, + [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_sum = + corsika::units::si::InverseGrammageType::zero()) { + + switch (select_.select(view.parent())) { + case SwitchResult::First: { + if constexpr (t1ProcSeq) { + // if A is a process sequence --> check inside + const EProcessReturn ret = + A.SelectInteraction(view, lambda_inv_select, lambda_inv_sum); + // if A did succeed, stop routine. Not checking other static branch B. + if (ret != EProcessReturn::eOk) { return ret; } + } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>, + TProcess1type>) { + // if this is not a ContinuousProcess --> evaluate probability + lambda_inv_sum += A.GetInverseInteractionLength(view.parent()); + // check if we should execute THIS process and then EXIT + if (lambda_inv_select < lambda_inv_sum) { + A.DoInteraction(view); + return EProcessReturn::eInteracted; + } + } // end branch A + break; + } + + case SwitchResult::Second: { + + if constexpr (t2ProcSeq) { + // if B is a process sequence --> check inside + return B.SelectInteraction(view, lambda_inv_select, lambda_inv_sum); + } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>, + TProcess2type>) { + // if this is not a ContinuousProcess --> evaluate probability + lambda_inv_sum += B.GetInverseInteractionLength(view.parent()); + // check if we should execute THIS process and then EXIT + if (lambda_inv_select < lambda_inv_sum) { + B.DoInteraction(view); + return EProcessReturn::eInteracted; + } + } // end branch B + break; + } + } + return EProcessReturn::eOk; + } + + template <typename TParticle> + inline corsika::units::si::TimeType GetLifetime(TParticle&& particle) { + return 1. / GetInverseLifetime(particle); + } + + template <typename TParticle> + inline corsika::units::si::InverseTimeType GetInverseLifetime(TParticle&& particle) { + using namespace corsika::units::si; + + switch (select_.select(particle)) { + case SwitchResult::First: { + if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, TProcess1type> || + t1ProcSeq) { + return A.GetInverseLifetime(particle); + } + break; + } + + case SwitchResult::Second: { + if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, TProcess2type> || + t2ProcSeq) { + return B.GetInverseLifetime(particle); + } + break; + } + } + return 0 / second; // default value + } + + // select decay process + template <typename TSecondaryView> + inline EProcessReturn SelectDecay( + TSecondaryView& view, + [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_select, + [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_sum = + corsika::units::si::InverseTimeType::zero()) { + + switch (select_.select(view.parent())) { + case SwitchResult::First: { + if constexpr (t1ProcSeq) { + // if A is a process sequence --> check inside + const EProcessReturn ret = A.SelectDecay(view, decay_inv_select, decay_inv_sum); + // if A did succeed, stop routine here (not checking other static branch B) + if (ret != EProcessReturn::eOk) { return ret; } + } else if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, + TProcess1type>) { + // if this is not a ContinuousProcess --> evaluate probability + decay_inv_sum += A.GetInverseLifetime(view.parent()); + // check if we should execute THIS process and then EXIT + if (decay_inv_select < decay_inv_sum) { + // more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot + A.DoDecay(view); + return EProcessReturn::eDecayed; + } + } // end branch A + break; + } + + case SwitchResult::Second: { + + if constexpr (t2ProcSeq) { + // if B is a process sequence --> check inside + return B.SelectDecay(view, decay_inv_select, decay_inv_sum); + } else if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, + TProcess2type>) { + // if this is not a ContinuousProcess --> evaluate probability + decay_inv_sum += B.GetInverseLifetime(view.parent()); + // check if we should execute THIS process and then EXIT + if (decay_inv_select < decay_inv_sum) { + B.DoDecay(view); + return EProcessReturn::eDecayed; + } + } // end branch B + break; + } + } + return EProcessReturn::eOk; + } + }; + + // the method `select(proc1,proc1,selector)` assembles many + // BaseProcesses, and ProcessSequences into a SwitchProcessSequence, + // all combinatorics must be allowed, this is why we define a macro + // to define all combinations here: + + // Both, Processes1 and Processes2, must derive from BaseProcesses + + template <typename TProcess1, typename TProcess2, typename TSelect> + inline typename std::enable_if< + std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>, + typename std::decay<TProcess1>::type>::value && + std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>, + typename std::decay<TProcess2>::type>::value, + SwitchProcessSequence<TProcess1, TProcess2, TSelect>>::type + select(TProcess1&& vA, TProcess2&& vB, TSelect selector) { + return SwitchProcessSequence<TProcess1, TProcess2, TSelect>(vA, vB, selector); + } + + /// traits marker to identify objectas ProcessSequence + template <typename TProcess1, typename TProcess2, typename TSelect> + struct is_process_sequence< + corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>> + : std::true_type {}; + + /// traits marker to identify objectas SwitchProcessSequence + template <typename TProcess1, typename TProcess2, typename TSelect> + struct is_switch_process_sequence< + corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>> + : std::true_type {}; + +} // namespace corsika::process diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc index b51876886..8b9f0b581 100644 --- a/Framework/ProcessSequence/testProcessSequence.cc +++ b/Framework/ProcessSequence/testProcessSequence.cc @@ -11,9 +11,11 @@ #include <array> #include <iomanip> #include <iostream> +#include <typeinfo> #include <corsika/process/ProcessSequence.h> -#include <corsika/process/switch_process/SwitchProcess.h> +#include <corsika/process/SwitchProcessSequence.h> +#include <corsika/units/PhysicalUnits.h> using namespace corsika; using namespace corsika::units::si; @@ -22,7 +24,12 @@ using namespace std; static const int nData = 10; -int globalCount = 0; +int globalCount = 0; // simple counter + +int checkDecay = 0; // use this as a bit field +int checkInteract = 0; // use this as a bit field +int checkSec = 0; // use this as a bit field +int checkCont = 0; // use this as a bit field class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> { int fV = 0; @@ -38,7 +45,8 @@ public: template <typename D, typename T> inline EProcessReturn DoContinuous(D& d, T&) const { cout << "ContinuousProcess1::DoContinuous" << endl; - for (int i = 0; i < nData; ++i) d.p[i] += 0.933; + checkCont |= 1; + for (int i = 0; i < nData; ++i) d.data_[i] += 0.933; return EProcessReturn::eOk; } }; @@ -56,7 +64,27 @@ public: template <typename D, typename T> inline EProcessReturn DoContinuous(D& d, T&) const { cout << "ContinuousProcess2::DoContinuous" << endl; - for (int i = 0; i < nData; ++i) d.p[i] += 0.933; + checkCont |= 2; + for (int i = 0; i < nData; ++i) d.data_[i] += 0.933; + return EProcessReturn::eOk; + } +}; + +class ContinuousProcess3 : public ContinuousProcess<ContinuousProcess3> { + int fV = 0; + +public: + ContinuousProcess3(const int v) + : fV(v) { + cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl; + globalCount++; + } + + template <typename D, typename T> + inline EProcessReturn DoContinuous(D& d, T&) const { + cout << "ContinuousProcess3::DoContinuous" << endl; + checkCont |= 4; + for (int i = 0; i < nData; ++i) d.data_[i] += 0.933; return EProcessReturn::eOk; } }; @@ -69,12 +97,18 @@ public: globalCount++; } - template <typename D, typename S> - inline EProcessReturn DoInteraction(D& d, S&) const { - for (int i = 0; i < nData; ++i) d.p[i] += 1 + i; + template <typename TView> + inline EProcessReturn DoInteraction(TView& v) const { + checkInteract |= 1; + for (int i = 0; i < nData; ++i) v.parent().data_[i] += 1 + i; return EProcessReturn::eOk; } + template <typename TParticle> + corsika::units::si::GrammageType GetInteractionLength(TParticle&) const { + return 10_g / square(1_cm); + } + private: int fV; }; @@ -89,15 +123,16 @@ public: globalCount++; } - template <typename Particle> - inline EProcessReturn DoInteraction(Particle&) const { + template <typename TView> + inline EProcessReturn DoInteraction(TView&) const { + checkInteract |= 2; cout << "Process2::DoInteraction" << endl; return EProcessReturn::eOk; } template <typename Particle> GrammageType GetInteractionLength(Particle&) const { cout << "Process2::GetInteractionLength" << endl; - return 3_g / (1_cm * 1_cm); + return 20_g / (1_cm * 1_cm); } }; @@ -111,15 +146,16 @@ public: globalCount++; } - template <typename Particle> - inline EProcessReturn DoInteraction(Particle&) const { + template <typename TView> + inline EProcessReturn DoInteraction(TView&) const { + checkInteract |= 4; cout << "Process3::DoInteraction" << endl; return EProcessReturn::eOk; } template <typename Particle> GrammageType GetInteractionLength(Particle&) const { cout << "Process3::GetInteractionLength" << endl; - return 1_g / (1_cm * 1_cm); + return 30_g / (1_cm * 1_cm); } }; @@ -135,11 +171,14 @@ public: template <typename D, typename T> inline EProcessReturn DoContinuous(D& d, T&) const { - for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; } + std::cout << "Base::DoContinuous" << std::endl; + checkCont |= 8; + for (int i = 0; i < nData; ++i) { d.data_[i] /= 1.2; } return EProcessReturn::eOk; } - template <typename Particle> - EProcessReturn DoInteraction(Particle&) const { + template <typename TView> + EProcessReturn DoInteraction(TView&) const { + checkInteract |= 8; return EProcessReturn::eOk; } }; @@ -156,8 +195,28 @@ public: TimeType GetLifetime(Particle&) const { return 1_s; } + template <typename TView> + EProcessReturn DoDecay(TView&) const { + checkDecay |= 1; + return EProcessReturn::eOk; + } +}; + +class Decay2 : public DecayProcess<Decay2> { + +public: + Decay2(const int) { + cout << "Decay2()" << endl; + globalCount++; + } + template <typename Particle> - EProcessReturn DoDecay(Particle&) const { + TimeType GetLifetime(Particle&) const { + return 2_s; + } + template <typename TView> + EProcessReturn DoDecay(TView&) const { + checkDecay |= 2; return EProcessReturn::eOk; } }; @@ -173,15 +232,22 @@ public: fCount++; return EProcessReturn::eOk; } - int GetCount() const { return fCount; } + int GetCount() const { return fCount; } }; struct DummyStack {}; struct DummyData { - double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; }; struct DummyTrajectory {}; +struct DummyView { + DummyView(DummyData& p) + : p_(p) {} + DummyData& p_; + DummyData& parent() { return p_; } +}; + TEST_CASE("Process Sequence", "[Process Sequence]") { SECTION("Check construction") { @@ -195,7 +261,9 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { Process4 m4(3); CHECK(globalCount == 4); - [[maybe_unused]] auto sequence = m1 << m2 << m3 << m4; + auto sequence = m1 % m2 % m3 % m4; + CHECK(is_process_sequence_v<decltype(sequence)> == true); + CHECK(is_process_sequence_v<decltype(m2)> == false); } SECTION("interaction length") { @@ -206,10 +274,9 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { DummyData particle; - auto sequence2 = cp1 << m2 << m3; - GrammageType const tot = sequence2.GetTotalInteractionLength(particle); - InverseGrammageType const tot_inv = - sequence2.GetTotalInverseInteractionLength(particle); + auto sequence2 = cp1 % m2 % m3; + GrammageType const tot = sequence2.GetInteractionLength(particle); + InverseGrammageType const tot_inv = sequence2.GetInverseInteractionLength(particle); cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl; globalCount = 0; } @@ -223,9 +290,9 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { DummyData particle; - auto sequence2 = cp1 << m2 << m3 << d3; - TimeType const tot = sequence2.GetTotalLifetime(particle); - InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(particle); + auto sequence2 = cp1 % m2 % m3 % d3; + TimeType const tot = sequence2.GetLifetime(particle); + InverseTimeType const tot_inv = sequence2.GetInverseLifetime(particle); cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl; globalCount = 0; @@ -238,7 +305,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { Process2 m2(2); Process3 m3(3); - auto sequence2 = cp1 << m2 << m3 << cp2; + auto sequence2 = cp1 % m2 % m3 % cp2; DummyData particle; DummyTrajectory track; @@ -255,7 +322,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { cout << "Running loop with n=" << nLoop << endl; for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(particle, track); } for (int i = 0; i < nData; i++) { - cout << "data[" << i << "]=" << particle.p[i] << endl; + cout << "data_[" << i << "]=" << particle.data_[i] << endl; } cout << "done" << endl; } @@ -266,7 +333,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { Stack1 s1(1); Stack1 s2(2); - auto sequence = s1 << s2; + auto sequence = s1 % s2; DummyStack stack; @@ -278,14 +345,106 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { } } -/* - Note: there is a fine-grained dedicated test-suite for SwitchProcess - in Processes/SwitchProcess/testSwtichProcess - */ -TEST_CASE("SwitchProcess") { - globalCount = 0; - Process1 p1(0); - Process2 p2(1); - switch_process::SwitchProcess s(p1, p2, 10_GeV); - CHECK(is_switch_process_v<decltype(s)>); +TEST_CASE("Switch Process Sequence", "[Process Sequence]") { + + SECTION("Check construction") { + + struct TestSelect { + corsika::process::SwitchResult select(const DummyData& p) const { + std::cout << "TestSelect data=" << p.data_[0] << std::endl; + if (p.data_[0] > 0) return corsika::process::SwitchResult::First; + return corsika::process::SwitchResult::Second; + } + }; + TestSelect select; + + auto sequence1 = Process1(0) % ContinuousProcess2(0) % Decay1(0); + auto sequence2 = ContinuousProcess3(0) % Process2(0) % Decay2(0); + + auto sequence = ContinuousProcess1(0) % Process3(0) % + SwitchProcessSequence(sequence1, sequence2, select); + + auto sequence_alt = + (ContinuousProcess1(0) % Process3(0)) % + process::select(Process1(0) % ContinuousProcess2(0) % Decay1(0), + ContinuousProcess3(0) % Process2(0) % Decay2(0), select); + + // check that same process sequence can be build in different ways + CHECK(typeid(sequence) == typeid(sequence_alt)); + CHECK(is_process_sequence_v<decltype(sequence)> == true); + CHECK(is_process_sequence_v<decltype( + SwitchProcessSequence(sequence1, sequence2, select))> == true); + + DummyData particle; + DummyTrajectory track; + DummyView view(particle); + + checkDecay = 0; + checkInteract = 0; + checkSec = 0; + checkCont = 0; + particle.data_[0] = 100; // data positive + sequence.DoContinuous(particle, track); + CHECK(checkInteract == 0); + CHECK(checkDecay == 0); + CHECK(checkCont == 0b011); + CHECK(checkSec == 0); + + checkDecay = 0; + checkInteract = 0; + checkSec = 0; + checkCont = 0; + particle.data_[0] = -100; // data negative + sequence_alt.DoContinuous(particle, track); + CHECK(checkInteract == 0); + CHECK(checkDecay == 0); + CHECK(checkCont == 0b101); + CHECK(checkSec == 0); + + // 1/(30g/cm2) is Process3 + corsika::units::si::InverseGrammageType lambda_select = .9/30. * square(1_cm) / 1_g; + corsika::units::si::InverseTimeType time_select = 0.1 / second; + + checkDecay = 0; + checkInteract = 0; + checkSec = 0; + checkCont = 0; + particle.data_[0] = 100; // data positive + sequence.SelectInteraction(view, lambda_select); + sequence.SelectDecay(view, time_select); + CHECK(checkInteract == 0b100); // this is Process3 + CHECK(checkDecay == 0b001); // this is Decay1 + CHECK(checkCont == 0); + CHECK(checkSec == 0); + lambda_select = 1.01/30. * square(1_cm) / 1_g; + checkInteract = 0; + sequence.SelectInteraction(view, lambda_select); + CHECK(checkInteract == 0b001); // this is Process1 + + checkDecay = 0; + checkInteract = 0; + checkSec = 0; + checkCont = 0; + particle.data_[0] = -100; // data negative + sequence.SelectInteraction(view, lambda_select); + sequence.SelectDecay(view, time_select); + CHECK(checkInteract == 0b010); // this is Process2 + CHECK(checkDecay == 0b010); // this is Decay2 + CHECK(checkCont == 0); + CHECK(checkSec == 0); + + checkDecay = 0; + checkInteract = 0; + checkSec = 0; + checkCont = 0; + particle.data_[0] = -100; // data negative + sequence.DoSecondaries(view); + Stack1 stack(0); + sequence.DoStack(stack); + CHECK(checkInteract == 0); + CHECK(checkDecay == 0); + CHECK(checkCont == 0); + CHECK(checkSec == 0); + } } + diff --git a/Framework/ProcessSequence/testSwitchProcessSequence.cc b/Framework/ProcessSequence/testSwitchProcessSequence.cc new file mode 100644 index 000000000..304dc78a0 --- /dev/null +++ b/Framework/ProcessSequence/testSwitchProcessSequence.cc @@ -0,0 +1,248 @@ +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <corsika/process/switch_process/SwitchProcess.h> +#include <corsika/stack/SecondaryView.h> +#include <corsika/stack/Stack.h> +#include <corsika/units/PhysicalUnits.h> + +#include <catch2/catch.hpp> + +#include <algorithm> +#include <random> + +using namespace corsika; +using namespace corsika::process; +using namespace corsika::units::si; + +class TestStackData { + +public: + // these functions are needed for the Stack interface + void Clear() { fData.clear(); } + unsigned int GetSize() const { return fData.size(); } + unsigned int GetCapacity() const { return fData.size(); } + void Copy(int i1, int i2) { fData[i2] = fData[i1]; } + void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); } + + // custom data access function + void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; } + HEPEnergyType GetData(const unsigned int i) const { return fData[i]; } + + // these functions are also needed by the Stack interface + void IncrementSize() { fData.resize(fData.size() + 1); } + void DecrementSize() { + if (fData.size() > 0) { fData.pop_back(); } + } + + // custom private data section +private: + std::vector<HEPEnergyType> fData; +}; + +/** + * From static_cast of a StackIterator over entries in the + * TestStackData class you get and object of type + * TestParticleInterface defined here + * + * It provides Get/Set methods to read and write data to the "Data" + * storage of TestStackData obtained via + * "StackIteratorInterface::GetStackData()", given the index of the + * iterator "StackIteratorInterface::GetIndex()" + * + */ +template <typename StackIteratorInterface> +class TestParticleInterface + : public corsika::stack::ParticleBase<StackIteratorInterface> { + +public: + using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData; + using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex; + + /* + The SetParticleData methods are called for creating new entries + on the stack. You can specifiy various parametric versions to + perform this task: + */ + + // default version for particle-creation from input data + void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); } + void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/, + std::tuple<HEPEnergyType> v) { + SetEnergy(std::get<0>(v)); + } + + // here are the fundamental methods for access to TestStackData data + void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); } + HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); } +}; + +using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>; + +// see issue 161 +#if defined(__clang__) +using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>; +#elif defined(__GNUC__) || defined(__GNUG__) +using StackTestView = corsika::stack::MakeView<SimpleStack>::type; +#endif + +auto constexpr kgMSq = 1_kg / (1_m * 1_m); + +template <int N> +struct DummyProcess : InteractionProcess<DummyProcess<N>> { + + template <typename TParticle> + corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const { + return N * kgMSq; + } + + template <typename TSecondaries> + corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) { + // to figure out which process was selected in the end, we produce N + // secondaries for DummyProcess<N> + + for (int i = 0; i < N; ++i) { + vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N}); + } + + return EProcessReturn::eOk; + } +}; + +using DummyLowEnergyProcess = DummyProcess<1>; +using DummyHighEnergyProcess = DummyProcess<2>; +using DummyAdditionalProcess = DummyProcess<3>; + +TEST_CASE("SwitchProcess from InteractionProcess") { + DummyLowEnergyProcess low; + DummyHighEnergyProcess high; + DummyAdditionalProcess proc; + + switch_process::SwitchProcess switchProcess(low, high, 1_TeV); + auto seq = switchProcess << proc; + + SimpleStack stack; + + SECTION("low energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV}); + auto p = stack.GetNextParticle(); + + // low energy process returns 1 kg/m² + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1)); + REQUIRE(seq.GetInteractionLength(p) / kgMSq == Approx(3. / 4)); + } + + } + + SECTION("high energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV}); + auto p = stack.GetNextParticle(); + + // high energy process returns 2 kg/m² + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2)); + REQUIRE(seq.GetInteractionLength(p) / kgMSq == Approx(6. / 5)); + } + + // high energy process creates 2 secondaries + SECTION("SelectInteraction") { + typename SimpleStack::ParticleType theParticle = + stack.GetNextParticle(); // as in corsika::Cascade + StackTestView view(theParticle); + auto projectile = view.GetProjectile(); + + InverseGrammageType invLambda = 0 / kgMSq; + switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda); + + REQUIRE(view.getSize() == 2); + } + } +} + +TEST_CASE("SwitchProcess from ProcessSequence") { + DummyProcess<1> innerA; + DummyProcess<2> innerB; + DummyProcess<3> outer; + DummyProcess<4> additional; + + auto seq = innerA << innerB; + + switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV); + auto completeSeq = switchProcess << additional; + + SimpleStack stack; + + SECTION("low energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV}); + auto p = stack.GetNextParticle(); + + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3)); + REQUIRE(completeSeq.GetInteractionLength(p) / kgMSq == Approx(4. / 7)); + } + + SECTION("SelectInteraction") { + std::vector<int> numberOfSecondaries; + + for (int i = 0; i < 1000; ++i) { + typename SimpleStack::ParticleType theParticle = + stack.GetNextParticle(); // as in corsika::Cascade + StackTestView view(theParticle); + auto projectile = view.GetProjectile(); + + double r = i / 1000.; + InverseGrammageType invLambda = r * 7. / 4 / kgMSq; + + InverseGrammageType accumulator = 0 / kgMSq; + completeSeq.SelectInteraction(p, projectile, invLambda, accumulator); + + numberOfSecondaries.push_back(view.getSize()); + } + + auto const mean = + std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) / + numberOfSecondaries.size(); + REQUIRE(mean == Approx(12. / 7.).margin(0.01)); + } + } + + SECTION("high energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV}); + auto p = stack.GetNextParticle(); + + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3)); + REQUIRE(completeSeq.GetInteractionLength(p) / kgMSq == Approx(12. / 7.)); + } + + SECTION("SelectInteraction") { + std::vector<int> numberOfSecondaries; + + for (int i = 0; i < 1000; ++i) { + typename SimpleStack::ParticleType theParticle = + stack.GetNextParticle(); // as in corsika::Cascade + StackTestView view(theParticle); + auto projectile = view.GetProjectile(); + + double r = i / 1000.; + InverseGrammageType invLambda = r * 7. / 12. / kgMSq; + + InverseGrammageType accumulator = 0 / kgMSq; + completeSeq.SelectInteraction(p, projectile, invLambda, accumulator); + + numberOfSecondaries.push_back(view.getSize()); + } + + auto const mean = + std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) / + numberOfSecondaries.size(); + REQUIRE(mean == Approx(24. / 7.).margin(0.01)); + } + } +} diff --git a/Processes/SwitchProcess/testSwitchProcess.cc b/Framework/ProcessSequence/testSwitchProcessSequence.h similarity index 94% rename from Processes/SwitchProcess/testSwitchProcess.cc rename to Framework/ProcessSequence/testSwitchProcessSequence.h index 9ee148e18..aed29021e 100644 --- a/Processes/SwitchProcess/testSwitchProcess.cc +++ b/Framework/ProcessSequence/testSwitchProcessSequence.h @@ -138,18 +138,6 @@ TEST_CASE("SwitchProcess from InteractionProcess") { REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4)); } - // low energy process creates 1 secondary - //~ SECTION("SelectInteraction") { - //~ typename SimpleStack::ParticleType theParticle = - //~ stack.GetNextParticle(); // as in corsika::Cascade - //~ StackTestView view(theParticle); - //~ auto projectile = view.GetProjectile(); - - //~ InverseGrammageType invLambda = 0 / kgMSq; - //~ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda); - - //~ REQUIRE(view.GetSize() == 1); - //~ } } SECTION("high energy") { diff --git a/Framework/Random/CMakeLists.txt b/Framework/Random/CMakeLists.txt index f9c72763c..f909f77f1 100644 --- a/Framework/Random/CMakeLists.txt +++ b/Framework/Random/CMakeLists.txt @@ -20,9 +20,8 @@ CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKArandom ${CORSIKArandom_NAMESPACE} ${CO target_link_libraries ( CORSIKArandom - INTERFACE CORSIKAutilities - CORSIKAunits + CORSIKAlogging ) target_include_directories ( @@ -50,4 +49,5 @@ target_link_libraries ( testRandom CORSIKArandom CORSIKAtesting + CORSIKAunits ) diff --git a/Framework/Random/RNGManager.cc b/Framework/Random/RNGManager.cc index 238b38636..9aa094a4b 100644 --- a/Framework/Random/RNGManager.cc +++ b/Framework/Random/RNGManager.cc @@ -7,6 +7,8 @@ */ #include <corsika/random/RNGManager.h> +#include <corsika/logging/Logging.h> + #include <sstream> void corsika::random::RNGManager::RegisterRandomStream(std::string const& pStreamName) { @@ -42,14 +44,21 @@ std::stringstream corsika::random::RNGManager::dumpState() const { } void corsika::random::RNGManager::SeedAll(uint64_t vSeed) { - for (auto& entry : rngs) { entry.second.seed(vSeed++); } + for (auto& entry : rngs) { + auto seed = vSeed++; + C8LOG_TRACE("Random seed stream {} seed {}", entry.first, seed); + entry.second.seed(seed); + } } void corsika::random::RNGManager::SeedAll() { std::random_device rd; - + std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()}; for (auto& entry : rngs) { - std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()}; - entry.second.seed(sseq); + std::vector<std::uint32_t> seeds(1); + sseq.generate(seeds.begin(), seeds.end()); + std::uint32_t seed = seeds[0]; + C8LOG_TRACE("Random seed stream {} seed {}", entry.first, seed); + entry.second.seed(seed); } } diff --git a/Framework/Random/UniformRealDistribution.h b/Framework/Random/UniformRealDistribution.h index 7e092076d..0c63e8542 100644 --- a/Framework/Random/UniformRealDistribution.h +++ b/Framework/Random/UniformRealDistribution.h @@ -8,7 +8,6 @@ #pragma once -#include <corsika/units/PhysicalUnits.h> #include <random> namespace corsika::random { diff --git a/Framework/StackInterface/SecondaryView.h b/Framework/StackInterface/SecondaryView.h index dcdcab812..985be9953 100644 --- a/Framework/StackInterface/SecondaryView.h +++ b/Framework/StackInterface/SecondaryView.h @@ -168,8 +168,8 @@ namespace corsika::stack { * SecondaryView is derived from. This projectile should not be * used to modify the Stack! */ - ConstStackIteratorValue parent() const { - return ConstStackIteratorValue(inner_stack_, projectile_index_); + StackIteratorValue parent() const { // todo: check if this can't be ConstStackIteratorValue + return StackIteratorValue(inner_stack_, projectile_index_); } /** diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h index e166fccba..33d5554ae 100644 --- a/Framework/StackInterface/Stack.h +++ b/Framework/StackInterface/Stack.h @@ -301,9 +301,9 @@ namespace corsika::stack { /** * check if this particle was already deleted */ - bool isDeleted(const StackIterator& p) { return isDeleted(p.GetIndex()); } + bool isDeleted(const StackIterator& p) const { return isDeleted(p.GetIndex()); } bool isDeleted(const ConstStackIterator& p) const { return isDeleted(p.GetIndex()); } - bool isDeleted(const ParticleInterfaceType& p) { return isDeleted(p.GetIterator()); } + bool isDeleted(const ParticleInterfaceType& p) const { return isDeleted(p.GetIterator()); } /** * Function to ultimatively remove the last entry from the stack, diff --git a/Framework/StackInterface/StackIteratorInterface.h b/Framework/StackInterface/StackIteratorInterface.h index 2229c17bc..c23561dc6 100644 --- a/Framework/StackInterface/StackIteratorInterface.h +++ b/Framework/StackInterface/StackIteratorInterface.h @@ -64,7 +64,7 @@ namespace corsika::stack { For two examples see stack_example.cc, or the corsika::processes::sibyll::SibStack class - */ + **/ template <typename TStackData, template <typename> typename TParticleInterface, typename StackType = Stack<TStackData, TParticleInterface>> @@ -105,6 +105,10 @@ namespace corsika::stack { StackIteratorInterface() = delete; public: + StackIteratorInterface(StackIteratorInterface&& rhs) + : index_(std::move(rhs.index_)) + , data_(std::move(rhs.data_)) {} + StackIteratorInterface(StackIteratorInterface const& vR) : index_(vR.index_) , data_(vR.data_) {} @@ -118,7 +122,7 @@ namespace corsika::stack { /** iterator must always point to data, with an index: @param data reference to the stack [rw] @param index index on stack - */ + **/ StackIteratorInterface(StackType& data, const unsigned int index) : index_(index) , data_(&data) {} @@ -129,7 +133,7 @@ namespace corsika::stack { @param args variadic list of data to initialize stack entry, this must be consistent with the definition of the user-provided ParticleInterfaceType::SetParticleData(...) function - */ + **/ template <typename... Args> StackIteratorInterface(StackType& data, const unsigned int index, const Args... args) : index_(index) @@ -146,7 +150,7 @@ namespace corsika::stack { @param args variadic list of data to initialize stack entry, this must be consistent with the definition of the user-provided ParticleInterfaceType::SetParticleData(...) function - */ + **/ template <typename... Args> StackIteratorInterface(StackType& data, const unsigned int index, StackIteratorInterface& parent, const Args... args) @@ -160,7 +164,7 @@ namespace corsika::stack { public: /** @name Iterator interface @{ - */ + **/ StackIteratorInterface& operator++() { do { ++index_; @@ -195,14 +199,15 @@ namespace corsika::stack { /** * Convert iterator to value type, where value type is the user-provided particle * readout class - */ + **/ ParticleInterfaceType& operator*() { return static_cast<ParticleInterfaceType&>(*this); } + /** * Convert iterator to const value type, where value type is the user-provided * particle readout class - */ + **/ const ParticleInterfaceType& operator*() const { return static_cast<const ParticleInterfaceType&>(*this); } @@ -212,7 +217,7 @@ namespace corsika::stack { /** * @name Stack data access * @{ - */ + **/ /// Get current particle index inline unsigned int GetIndex() const { return index_; } /// Get current particle Stack object @@ -234,7 +239,7 @@ namespace corsika::stack { @class ConstStackIteratorInterface This is the iterator class for const-access to stack data - */ + **/ template <typename TStackData, template <typename> typename TParticleInterface, typename StackType = Stack<TStackData, TParticleInterface>> @@ -275,6 +280,10 @@ namespace corsika::stack { ConstStackIteratorInterface() = delete; public: + ConstStackIteratorInterface(ConstStackIteratorInterface&& rhs) + : index_(std::move(rhs.index_)) + , data_(std::move(rhs.data_)) {} + ConstStackIteratorInterface(const StackType& data, const unsigned int index) : index_(index) , data_(&data) {} @@ -290,13 +299,13 @@ namespace corsika::stack { \endverbatim See documentation of StackIteratorInterface for more details. - */ + **/ bool isDeleted() const { return GetStack().isDeleted(*this); } public: /** @name Iterator interface - */ + **/ ///@{ ConstStackIteratorInterface& operator++() { do { @@ -339,7 +348,7 @@ namespace corsika::stack { protected: /** @name Stack data access Only the const versions for read-only access - */ + **/ ///@{ inline unsigned int GetIndex() const { return index_; } inline const StackType& GetStack() const { return *data_; } diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index 5deaf153c..53a01ca25 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -31,7 +31,6 @@ add_subdirectory (ParticleCut) add_subdirectory (OnShellCheck) # meta-processes add_subdirectory (InteractionCounter) -add_subdirectory (SwitchProcess) ########################################## diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc index 69dcd7877..107fa813a 100644 --- a/Processes/EnergyLoss/EnergyLoss.cc +++ b/Processes/EnergyLoss/EnergyLoss.cc @@ -177,15 +177,10 @@ process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack co cout << "EnergyLoss dE=" << dE / 1_MeV << "MeV, " << " E=" << E / 1_GeV << "GeV, Ekin=" << Ekin / 1_GeV << ", Enew=" << Enew / 1_GeV << "GeV" << endl; - auto status = process::EProcessReturn::eOk; - if (E < emCut_) { - Enew = emCut_; - status = process::EProcessReturn::eParticleAbsorbed; - } p.SetEnergy(Enew); MomentumUpdate(p, Enew); FillProfile(t, dE); - return status; + return process::EProcessReturn::eOk; } LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle, diff --git a/Processes/LongitudinalProfile/LongitudinalProfile.cc b/Processes/LongitudinalProfile/LongitudinalProfile.cc index 0e174ee73..5c9956663 100644 --- a/Processes/LongitudinalProfile/LongitudinalProfile.cc +++ b/Processes/LongitudinalProfile/LongitudinalProfile.cc @@ -11,6 +11,8 @@ #include <corsika/particles/ParticleProperties.h> +#include <corsika/logging/Logging.h> + #include <corsika/setup/SetupStack.h> #include <corsika/setup/SetupTrajectory.h> @@ -39,6 +41,8 @@ corsika::process::EProcessReturn LongitudinalProfile::DoContinuous(Particle cons GrammageType const grammageStart = shower_axis_.projectedX(vTrack.GetPosition(0)); GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1)); + C8LOG_INFO("pos1={} m, pos2={}, X={} g/cm2", vTrack.GetPosition(0).GetCoordinates()/1_m, vTrack.GetPosition(1).GetCoordinates()/1_m, grammageStart/1_g*square(1_cm)); + const int binStart = std::ceil(grammageStart / dX_); const int binEnd = std::floor(grammageEnd / dX_); diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc index cb7ee95bd..ab3a80232 100644 --- a/Processes/ObservationPlane/ObservationPlane.cc +++ b/Processes/ObservationPlane/ObservationPlane.cc @@ -7,9 +7,9 @@ */ #include <corsika/process/observation_plane/ObservationPlane.h> +#include <corsika/logging/Logging.h> #include <fstream> -#include <iostream> using namespace corsika::process::observation_plane; using namespace corsika::units::si; @@ -25,7 +25,7 @@ ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane, } corsika::process::EProcessReturn ObservationPlane::DoContinuous( - setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) { + setup::Stack::ParticleType& particle, setup::Trajectory const& trajectory) { TimeType const timeOfIntersection = (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) / trajectory.GetV0().dot(plane_.GetNormal()); @@ -45,6 +45,7 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous( if (deleteOnHit_) { count_ground_++; energy_ground_ += energy; + particle.Delete(); return process::EProcessReturn::eParticleAbsorbed; } else { return process::EProcessReturn::eOk; @@ -62,15 +63,19 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&, } auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection); - return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001; + auto dist = (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001; + C8LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m); + return dist; } 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; + C8LOG_INFO( + " ******************************\n" + " ObservationPlane: \n" + " energy in ground (GeV) : {}\n" + " no. of particles in ground : {}\n" + " ******************************", + energy_ground_ / 1_GeV, count_ground_); } void ObservationPlane::Reset() { diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h index 7223defb1..86c5ba0f2 100644 --- a/Processes/ObservationPlane/ObservationPlane.h +++ b/Processes/ObservationPlane/ObservationPlane.h @@ -29,7 +29,7 @@ namespace corsika::process::observation_plane { ObservationPlane(geometry::Plane const&, std::string const&, bool = true); corsika::process::EProcessReturn DoContinuous( - corsika::setup::Stack::ParticleType const& vParticle, + corsika::setup::Stack::ParticleType& vParticle, corsika::setup::Trajectory const& vTrajectory); corsika::units::si::LengthType MaxStepLength( diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc index 72a97ee06..e5fb9d47e 100644 --- a/Processes/ParticleCut/ParticleCut.cc +++ b/Processes/ParticleCut/ParticleCut.cc @@ -27,9 +27,9 @@ namespace corsika::process { if (vP.GetPID() == particles::Code::Nucleus) { // calculate energy per nucleon auto const ElabNuc = energyLab / vP.GetNuclearA(); - return (ElabNuc <= fECut); + return (ElabNuc <= energy_cut_); } else { - return (energyLab <= fECut); + return (energyLab <= energy_cut_); } } @@ -87,13 +87,12 @@ namespace corsika::process { return false; // this particle will not be removed/cut } - EProcessReturn ParticleCut::DoSecondaries(corsika::setup::StackView& vS) { + void ParticleCut::DoSecondaries(corsika::setup::StackView& vS) { auto particle = vS.begin(); while (particle != vS.end()) { if (checkCutParticle(particle)) { particle.Delete(); } ++particle; // next entry in SecondaryView } - return EProcessReturn::eOk; } process::EProcessReturn ParticleCut::DoContinuous( @@ -102,13 +101,15 @@ namespace corsika::process { C8LOG_TRACE("ParticleCut::DoContinuous"); if (checkCutParticle(particle)) { C8LOG_TRACE("removing during continuous"); - return process::EProcessReturn::eParticleAbsorbed; + particle.Delete(); + // signal to upstream code that this particle was deleted + return process::EProcessReturn::eParticleAbsorbed; } return process::EProcessReturn::eOk; } ParticleCut::ParticleCut(const units::si::HEPEnergyType eCut, bool em, bool inv) - : fECut(eCut) + : energy_cut_(eCut) , bCutEm(em) , bCutInv(inv) { fEmEnergy = 0_GeV; diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h index 4705051d4..b0e897753 100644 --- a/Processes/ParticleCut/ParticleCut.h +++ b/Processes/ParticleCut/ParticleCut.h @@ -12,6 +12,7 @@ #include <corsika/process/ContinuousProcess.h> #include <corsika/process/SecondariesProcess.h> #include <corsika/setup/SetupStack.h> +#include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> namespace corsika::process { @@ -19,7 +20,7 @@ namespace corsika::process { class ParticleCut : public process::SecondariesProcess<ParticleCut>, public corsika::process::ContinuousProcess<ParticleCut> { - units::si::HEPEnergyType const fECut; + units::si::HEPEnergyType const energy_cut_; bool bCutEm; bool bCutInv; @@ -32,7 +33,7 @@ namespace corsika::process { public: ParticleCut(const units::si::HEPEnergyType eCut, bool em, bool inv); - EProcessReturn DoSecondaries(corsika::setup::StackView&); + void DoSecondaries(corsika::setup::StackView&); EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType& vParticle, corsika::setup::Trajectory const& vTrajectory); @@ -42,7 +43,7 @@ namespace corsika::process { return units::si::meter * std::numeric_limits<double>::infinity(); } - units::si::HEPEnergyType GetECut() const { return fECut; } + units::si::HEPEnergyType GetECut() const { return energy_cut_; } units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; } units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; } units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; } diff --git a/Processes/Proposal/CMakeLists.txt b/Processes/Proposal/CMakeLists.txt index cc0be5438..08a37a68b 100644 --- a/Processes/Proposal/CMakeLists.txt +++ b/Processes/Proposal/CMakeLists.txt @@ -22,6 +22,7 @@ TARGET_LINK_LIBRARIES ( CORSIKAunits CORSIKAthirdparty CORSIKAgeometry + CORSIKAlogging CORSIKAenvironment ProcessParticleCut PROPOSAL::PROPOSAL diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc index 59ec7e87f..d0e2056d5 100644 --- a/Processes/Proposal/ContinuousProcess.cc +++ b/Processes/Proposal/ContinuousProcess.cc @@ -15,6 +15,7 @@ #include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> #include <corsika/utl/COMBoost.h> +#include <corsika/logging/Logging.h> #include <iostream> @@ -106,21 +107,8 @@ namespace corsika::process::proposal { // if the particle has a charge take multiple scattering into account 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. - if (final_energy <= emCut_) { - vP.SetEnergy(emCut_); - vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm()); - return process::EProcessReturn::eParticleAbsorbed; - } - if (final_energy <= emCut_) { - vP.SetEnergy(emCut_); - vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm()); - return process::EProcessReturn::eParticleAbsorbed; - } vP.SetEnergy(final_energy); - vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm()); + vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm()); return process::EProcessReturn::eOk; } @@ -132,7 +120,7 @@ namespace corsika::process::proposal { 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 + // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a // hyper parameter which must be adjusted. // // in any case: never go below 0.99*emCut_ This needs to be @@ -152,7 +140,9 @@ namespace corsika::process::proposal { 1_g / square(1_cm); // return it in distance aequivalent - return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage); + auto dist = vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage); + C8LOG_TRACE("PROPOSAL::MaxStepLength X={} g/cm2, l={} m ", grammage/1_g*square(1_cm), dist/1_m); + return dist; } void ContinuousProcess::ShowResults() const { diff --git a/Processes/Pythia/Decay.h b/Processes/Pythia/Decay.h index fb4708df1..733b9eae1 100644 --- a/Processes/Pythia/Decay.h +++ b/Processes/Pythia/Decay.h @@ -11,6 +11,8 @@ #include <Pythia8/Pythia.h> #include <corsika/particles/ParticleProperties.h> #include <corsika/process/DecayProcess.h> +#include <corsika/units/PhysicalUnits.h> +#include <corsika/geometry/FourVector.h> namespace corsika::process { diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc index fe641a626..da9ef31e9 100644 --- a/Processes/Pythia/Interaction.cc +++ b/Processes/Pythia/Interaction.cc @@ -164,7 +164,7 @@ namespace corsika::process::pythia { } template <> - units::si::GrammageType Interaction::GetInteractionLength(Particle& p) { + units::si::GrammageType Interaction::GetInteractionLength(const Particle& p) { using namespace units; using namespace units::si; diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h index 4b06aeb6d..67a6ce423 100644 --- a/Processes/Pythia/Interaction.h +++ b/Processes/Pythia/Interaction.h @@ -47,7 +47,7 @@ namespace corsika::process::pythia { const corsika::units::si::HEPEnergyType CoMenergy); template <typename TParticle> - corsika::units::si::GrammageType GetInteractionLength(TParticle&); + corsika::units::si::GrammageType GetInteractionLength(const TParticle&); /** In this function PYTHIA is called to produce one event. The diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc index 62d0b4581..330dcfd27 100644 --- a/Processes/Sibyll/Decay.cc +++ b/Processes/Sibyll/Decay.cc @@ -23,7 +23,7 @@ using namespace corsika::setup; using SetupView = corsika::setup::StackView; using SetupProjectile = corsika::setup::StackView::ParticleType; -using SetupParticle = corsika::setup::Stack::ParticleType; +using Particle = corsika::setup::Stack::ParticleType; namespace corsika::process::sibyll { @@ -142,7 +142,7 @@ namespace corsika::process::sibyll { } template <> - units::si::TimeType Decay::GetLifetime(SetupParticle const& vP) { + units::si::TimeType Decay::GetLifetime(Particle const& vP) { using namespace units::si; const particles::Code pid = vP.GetPID(); diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index b028ca876..8ba3ab187 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/Processes/Sibyll/Interaction.cc @@ -25,7 +25,7 @@ using std::tuple; using namespace corsika; using namespace corsika::setup; -using SetupParticle = setup::Stack::StackIterator; +using Particle = setup::Stack::StackIterator; using SetupView = setup::StackView; using Track = Trajectory; @@ -81,8 +81,7 @@ namespace corsika::process::sibyll { } template <> - units::si::GrammageType Interaction::GetInteractionLength( - SetupParticle const& vP) const { + units::si::GrammageType Interaction::GetInteractionLength(Particle const& vP) const { using namespace units; using namespace units::si; @@ -172,187 +171,184 @@ namespace corsika::process::sibyll { auto const projectile = view.GetProjectile(); const auto corsikaBeamId = projectile.GetPID(); - C8LOG_DEBUG( - fmt::format("ProcessSibyll: " - "DoInteraction: {} interaction? ", - corsikaBeamId, process::sibyll::CanInteract(corsikaBeamId))); if (particles::IsNucleus(corsikaBeamId)) { // nuclei handled by different process, this should not happen throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!"); } - if (process::sibyll::CanInteract(corsikaBeamId)) { - // position and time of interaction, not used in Sibyll - Point const pOrig = projectile.GetPosition(); - TimeType const tOrig = projectile.GetTime(); + // position and time of interaction, not used in Sibyll + Point const pOrig = projectile.GetPosition(); + TimeType const tOrig = projectile.GetTime(); - // define projectile - HEPEnergyType const eProjectileLab = projectile.GetEnergy(); - auto const pProjectileLab = projectile.GetMomentum(); - const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem(); + // define projectile + HEPEnergyType const eProjectileLab = projectile.GetEnergy(); + auto const pProjectileLab = projectile.GetMomentum(); + const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem(); - // define target - // for Sibyll is always a single nucleon - // FOR NOW: target is always at rest - const auto eTargetLab = 0_GeV + constants::nucleonMass; - const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV); - const FourVector PtargLab(eTargetLab, pTargetLab); + C8LOG_DEBUG( + "ProcessSibyll: " + "DoInteraction: pid {} interaction ", + corsikaBeamId); - C8LOG_DEBUG( - fmt::format("Interaction: ebeam lab: {} GeV" - "Interaction: pbeam lab: {} GeV", - eProjectileLab / 1_GeV, pProjectileLab.GetComponents())); - C8LOG_DEBUG( - fmt::format("Interaction: etarget lab: {} GeV " - "Interaction: ptarget lab: {} GeV", - eTargetLab / 1_GeV, pTargetLab.GetComponents() / 1_GeV)); + // define target + // for Sibyll is always a single nucleon + // FOR NOW: target is always at rest + const auto eTargetLab = 0_GeV + constants::nucleonMass; + const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV); + const FourVector PtargLab(eTargetLab, pTargetLab); - const FourVector PprojLab(eProjectileLab, pProjectileLab); + C8LOG_DEBUG( + "Interaction: ebeam lab: {} GeV" + "Interaction: pbeam lab: {} GeV", + eProjectileLab / 1_GeV, pProjectileLab.GetComponents()); + C8LOG_DEBUG( + "Interaction: etarget lab: {} GeV " + "Interaction: ptarget lab: {} GeV", + eTargetLab / 1_GeV, pTargetLab.GetComponents() / 1_GeV); - // define target kinematics in lab frame - // define boost to and from CoM frame - // CoM frame definition in Sibyll projectile: +z - COMBoost const boost(PprojLab, constants::nucleonMass); - auto const& csPrime = boost.GetRotatedCS(); + const FourVector PprojLab(eProjectileLab, pProjectileLab); - // just for show: - // boost projecticle - auto const PprojCoM = boost.toCoM(PprojLab); + // define target kinematics in lab frame + // define boost to and from CoM frame + // CoM frame definition in Sibyll projectile: +z + COMBoost const boost(PprojLab, constants::nucleonMass); + auto const& csPrime = boost.GetRotatedCS(); - // boost target - auto const PtargCoM = boost.toCoM(PtargLab); + // just for show: + // boost projecticle + auto const PprojCoM = boost.toCoM(PprojLab); - C8LOG_DEBUG( - fmt::format("Interaction: ebeam CoM: {} GeV " - "Interaction: pbeam CoM: {} GeV ", - PprojCoM.GetTimeLikeComponent() / 1_GeV, - PprojCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV)); - C8LOG_DEBUG( - fmt::format("Interaction: etarget CoM: {} GeV " - "Interaction: ptarget CoM: {} GeV ", - PtargCoM.GetTimeLikeComponent() / 1_GeV, - PtargCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV)); - - C8LOG_DEBUG(fmt::format("Interaction: position of interaction: {} ", - pOrig.GetCoordinates())); - C8LOG_DEBUG(fmt::format("Interaction: time: {} ", tOrig)); - - HEPEnergyType Etot = eProjectileLab + eTargetLab; - MomentumVector Ptot = projectile.GetMomentum(); - // invariant mass, i.e. cm. energy - HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm()); - - // sample target mass number - auto const* currentNode = projectile.GetNode(); - auto const& mediumComposition = - currentNode->GetModelProperties().GetNuclearComposition(); - // get cross sections for target materials - /* - Here we read the cross section from the interaction model again, - should be passed from GetInteractionLength if possible - */ - //#warning reading interaction cross section again, should not be necessary - auto const& compVec = mediumComposition.GetComponents(); - std::vector<CrossSectionType> cross_section_of_components(compVec.size()); - - for (size_t i = 0; i < compVec.size(); ++i) { - auto const targetId = compVec[i]; - const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm); - [[maybe_unused]] const auto& dummy_sigEla = sigEla; - cross_section_of_components[i] = sigProd; - } + // boost target + auto const PtargCoM = boost.toCoM(PtargLab); - const auto targetCode = - mediumComposition.SampleTarget(cross_section_of_components, RNG_); - C8LOG_DEBUG(fmt::format("Interaction: target selected: {} ", targetCode)); - /* - FOR NOW: allow nuclei with A<18 or protons only. - when medium composition becomes more complex, approximations will have to be - allowed air in atmosphere also contains some Argon. - */ - int targetSibCode = -1; - if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode); - if (targetCode == particles::Proton::GetCode()) targetSibCode = 1; - C8LOG_DEBUG(fmt::format("Interaction: sibyll code: {}", targetSibCode)); - if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1) - throw std::runtime_error( - "Sibyll target outside range. Only nuclei with A<18 or protons are " - "allowed."); + C8LOG_DEBUG( + "Interaction: ebeam CoM: {} GeV " + "Interaction: pbeam CoM: {} GeV ", + PprojCoM.GetTimeLikeComponent() / 1_GeV, + PprojCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV); + C8LOG_DEBUG( + "Interaction: etarget CoM: {} GeV " + "Interaction: ptarget CoM: {} GeV ", + PtargCoM.GetTimeLikeComponent() / 1_GeV, + PtargCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV); + + C8LOG_DEBUG("Interaction: position of interaction: {} ", pOrig.GetCoordinates()); + C8LOG_DEBUG("Interaction: time: {} ", tOrig); + + HEPEnergyType Etot = eProjectileLab + eTargetLab; + MomentumVector Ptot = projectile.GetMomentum(); + // invariant mass, i.e. cm. energy + HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm()); + + // sample target mass number + auto const* currentNode = projectile.GetNode(); + auto const& mediumComposition = + currentNode->GetModelProperties().GetNuclearComposition(); + // get cross sections for target materials + /* + Here we read the cross section from the interaction model again, + should be passed from GetInteractionLength if possible + */ + //#warning reading interaction cross section again, should not be necessary + auto const& compVec = mediumComposition.GetComponents(); + std::vector<CrossSectionType> cross_section_of_components(compVec.size()); + + for (size_t i = 0; i < compVec.size(); ++i) { + auto const targetId = compVec[i]; + const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm); + [[maybe_unused]] const auto& dummy_sigEla = sigEla; + cross_section_of_components[i] = sigProd; + } - // beam id for sibyll - const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId); + const auto targetCode = + mediumComposition.SampleTarget(cross_section_of_components, RNG_); + C8LOG_DEBUG("Interaction: target selected: {} ", targetCode); + /* + FOR NOW: allow nuclei with A<18 or protons only. + when medium composition becomes more complex, approximations will have to be + allowed air in atmosphere also contains some Argon. + */ + int targetSibCode = -1; + if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode); + if (targetCode == particles::Proton::GetCode()) targetSibCode = 1; + C8LOG_DEBUG("Interaction: sibyll code: {}", targetSibCode); + if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1) + throw std::runtime_error( + "Sibyll target outside range. Only nuclei with A<18 or protons are " + "allowed."); + // beam id for sibyll + const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId); + + C8LOG_DEBUG( + "Interaction: " + " DoInteraction: E(GeV): {} " + " Ecm(GeV): {} ", + eProjectileLab / 1_GeV, Ecm / 1_GeV); + if (Ecm > GetMaxEnergyCoM()) + throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!"); + // FR: removed eProjectileLab < 8.5_GeV || + if (Ecm < GetMinEnergyCoM()) { C8LOG_DEBUG( - fmt::format("Interaction: " - " DoInteraction: E(GeV): {} " - " Ecm(GeV): {} ", - eProjectileLab / 1_GeV, Ecm / 1_GeV)); - if (Ecm > GetMaxEnergyCoM()) - throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!"); - // FR: removed eProjectileLab < 8.5_GeV || - if (Ecm < GetMinEnergyCoM()) { - C8LOG_DEBUG( - fmt::format("Interaction: " - " DoInteraction: should have dropped particle.. " - "THIS IS AN ERROR")); - throw std::runtime_error("energy too low for SIBYLL"); - } else { - count_++; - // Sibyll does not know about units.. - const double sqs = Ecm / 1_GeV; - // running sibyll, filling stack - sibyll_(kBeam, targetSibCode, sqs); - - // print final state - int print_unit = 6; - sib_list_(print_unit); - nucCount_ += get_nwounded() - 1; - - // add particles from sibyll to stack - // link to sibyll stack - SibStack ss; - - MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV}); - HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV; - for (auto& psib : ss) { - - // abort on particles that have decayed in Sibyll. Should not happen! - if (psib.HasDecayed()) - throw std::runtime_error("found particle that decayed in SIBYLL!"); - - // transform 4-momentum to lab. frame - // note that the momentum needs to be rotated back - auto const tmp = psib.GetMomentum().GetComponents(); - auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp); - HEPEnergyType const eCoM = psib.GetEnergy(); - auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM)); - auto const p3lab = Plab.GetSpaceLikeComponents(); - assert(p3lab.GetCoordinateSystem() == originalCS); // just to be sure! - - // add to corsika stack - auto pnew = view.AddSecondary( - make_tuple(process::sibyll::ConvertFromSibyll(psib.GetPID()), - Plab.GetTimeLikeComponent(), p3lab, pOrig, tOrig)); - - Plab_final += pnew.GetMomentum(); - Elab_final += pnew.GetEnergy(); - Ecm_final += psib.GetEnergy(); - } - C8LOG_DEBUG(fmt::format( - "conservation (all GeV):" - "Ecm_initial(per nucleon)={}, Ecm_final(per nucleon)={}, " - "Elab_initial={}, Elab_final={}, " - "diff (%)={}, " - "E in nucleons={}, " - "Plab_initial={}, " - "Plab_final={} ", - Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV, - Elab_final / 1_GeV, (Elab_final / Etot / get_nwounded() - 1) * 100, - constants::nucleonMass * get_nwounded() / 1_GeV, - (pProjectileLab / 1_GeV).GetComponents(), - (Plab_final / 1_GeV).GetComponents())); + "Interaction: " + " DoInteraction: should have dropped particle.. " + "THIS IS AN ERROR"); + throw std::runtime_error("energy too low for SIBYLL"); + } else { + count_++; + // Sibyll does not know about units.. + const double sqs = Ecm / 1_GeV; + // running sibyll, filling stack + sibyll_(kBeam, targetSibCode, sqs); + + // print final state + int print_unit = 6; + sib_list_(print_unit); + nucCount_ += get_nwounded() - 1; + + // add particles from sibyll to stack + // link to sibyll stack + SibStack ss; + + MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV}); + HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV; + for (auto& psib : ss) { + + // abort on particles that have decayed in Sibyll. Should not happen! + if (psib.HasDecayed()) + throw std::runtime_error("found particle that decayed in SIBYLL!"); + + // transform 4-momentum to lab. frame + // note that the momentum needs to be rotated back + auto const tmp = psib.GetMomentum().GetComponents(); + auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp); + HEPEnergyType const eCoM = psib.GetEnergy(); + auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM)); + auto const p3lab = Plab.GetSpaceLikeComponents(); + assert(p3lab.GetCoordinateSystem() == originalCS); // just to be sure! + + // add to corsika stack + auto pnew = view.AddSecondary( + make_tuple(process::sibyll::ConvertFromSibyll(psib.GetPID()), + Plab.GetTimeLikeComponent(), p3lab, pOrig, tOrig)); + + Plab_final += pnew.GetMomentum(); + Elab_final += pnew.GetEnergy(); + Ecm_final += psib.GetEnergy(); } + C8LOG_DEBUG( + "conservation (all GeV):" + "Ecm_initial(per nucleon)={}, Ecm_final(per nucleon)={}, " + "Elab_initial={}, Elab_final={}, " + "diff (%)={}, " + "E in nucleons={}, " + "Plab_initial={}, " + "Plab_final={} ", + Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV, + Elab_final / 1_GeV, (Elab_final / Etot / get_nwounded() - 1) * 100, + constants::nucleonMass * get_nwounded() / 1_GeV, + (pProjectileLab / 1_GeV).GetComponents(), (Plab_final / 1_GeV).GetComponents()); } return process::EProcessReturn::eOk; } diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h index cccc971bd..3bcf68807 100644 --- a/Processes/Sibyll/Interaction.h +++ b/Processes/Sibyll/Interaction.h @@ -45,7 +45,7 @@ namespace corsika::process::sibyll { const corsika::units::si::HEPEnergyType) const; template <typename TParticle> - corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const; + corsika::units::si::GrammageType GetInteractionLength(const TParticle&) const; /** In this function SIBYLL is called to produce one event. The diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc index da4a91baa..6bd1bceb3 100644 --- a/Processes/StackInspector/StackInspector.cc +++ b/Processes/StackInspector/StackInspector.cc @@ -40,7 +40,7 @@ template <typename TStack> StackInspector<TStack>::~StackInspector() {} template <typename TStack> -process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) { +void StackInspector<TStack>::DoStack(const TStack& vS) { [[maybe_unused]] int i = 0; HEPEnergyType Etot = 0_GeV; @@ -63,7 +63,7 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) { const std::chrono::duration<double> elapsed_seconds = now - StartTime_; std::time_t const now_time = std::chrono::system_clock::to_time_t(now); auto const dE = E0_ - Etot; - if (dE < dE_threshold_) return process::EProcessReturn::eOk; + if (dE < dE_threshold_) return; double const progress = dE / E0_; double const eta_seconds = elapsed_seconds.count() / progress; @@ -77,7 +77,7 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) { << ", nStep=" << GetStep() << ", stackEntries=" << vS.getEntries() << ", Estack=" << Etot / 1_GeV << " GeV" << ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << endl; - return process::EProcessReturn::eOk; + return; } #include <corsika/cascade/testCascade.h> diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h index 280ff15e6..a7222d37c 100644 --- a/Processes/StackInspector/StackInspector.h +++ b/Processes/StackInspector/StackInspector.h @@ -29,8 +29,8 @@ namespace corsika::process { StackInspector(const int vNStep, const bool vReportStack, const corsika::units::si::HEPEnergyType vE0); ~StackInspector(); - - EProcessReturn DoStack(const TStack&); + + void DoStack(const TStack&); /** * To set a new E0, for example when a new shower event is started diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc index 89edc2a44..d76d9fd09 100644 --- a/Processes/StackInspector/testStackInspector.cc +++ b/Processes/StackInspector/testStackInspector.cc @@ -46,7 +46,6 @@ TEST_CASE("StackInspector", "[processes]") { SECTION("interface") { StackInspector<TestCascadeStack> model(1, true, E0); - - [[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack); + model.DoStack(stack); } } diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc index 268142932..71bf1dd49 100644 --- a/Processes/UrQMD/UrQMD.cc +++ b/Processes/UrQMD/UrQMD.cc @@ -19,7 +19,6 @@ #include <cmath> #include <fstream> #include <functional> -#include <iostream> #include <random> #include <sstream> @@ -40,10 +39,6 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode, HEPEnergyType labEnergy) const { // translated to C++ from CORSIKA 7 subroutine cxtot_u - C8LOG_DEBUG("UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={}GeV", - particles::GetName(projectileCode), particles::GetName(targetCode), - labEnergy / 1_GeV); - auto const kinEnergy = labEnergy - particles::GetMass(projectileCode); assert(kinEnergy >= HEPEnergyType::zero()); @@ -91,10 +86,10 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode, case particles::Code::K0Bar: projectileIndex = 8; break; - default: - std::cout << "WARNING: UrQMD cross-section not tabulated for " << projectileCode - << std::endl; + default: { + C8LOG_WARN("WARNING: UrQMD cross-section not tabulated for {}", projectileCode); return CrossSectionType::zero(); + } } int targetIndex; @@ -120,6 +115,10 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode, xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i]; } + C8LOG_DEBUG("UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={}GeV, sigma={}", + particles::GetName(projectileCode), particles::GetName(targetCode), + labEnergy / 1_GeV, result); + return result; } @@ -248,14 +247,16 @@ GrammageType UrQMD::GetInteractionLength(SetupParticle const& particle) const { corsika::process::EProcessReturn UrQMD::DoInteraction(SetupView& view) { using namespace units::si; - auto const projectile = view.GetProjectile(); - + auto const projectile = view.GetProjectile(); + auto projectileCode = projectile.GetPID(); auto const projectileEnergyLab = projectile.GetEnergy(); auto const& projectileMomentumLab = projectile.GetMomentum(); auto const& projectilePosition = projectile.GetPosition(); auto const projectileTime = projectile.GetTime(); + C8LOG_DEBUG("UrQMD::DoInteraction pid={} E={} GeV", projectileCode, projectileEnergyLab/1_GeV); + // sample target particle auto const& mediumComposition = projectile.GetNode()->GetModelProperties().GetNuclearComposition(); @@ -348,14 +349,14 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupView& view) { auto const energy = sqrt(momentum.squaredNorm() + square(particles::GetMass(code))); momentum.rebase(originalCS); // transform back into standard lab frame - std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl; + C8LOG_DEBUG(" Secondary {} code {} p={} GeV", i, code, + momentum.GetComponents() / 1_GeV); view.AddSecondary( - std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point, - TimeType>{code, energy, momentum, projectilePosition, projectileTime}); + std::make_tuple(code, energy, momentum, projectilePosition, projectileTime)); } - std::cout << "UrQMD generated " << sys_.npart << " secondaries!" << std::endl; + C8LOG_DEBUG("UrQMD generated {} secondaries!", sys_.npart); return process::EProcessReturn::eOk; } -- GitLab