From 98f7a4ecae988b8c66e02d66d1c929f7856a49b1 Mon Sep 17 00:00:00 2001 From: ralfulrich <ralf.ulrich@kit.edu> Date: Tue, 2 Apr 2019 11:30:45 +0200 Subject: [PATCH] style guide --- Framework/Cascade/CMakeLists.txt | 3 +- Framework/Cascade/Cascade.h | 99 ++++++++++++++++++++------------ Framework/Cascade/testCascade.cc | 14 +++-- 3 files changed, 71 insertions(+), 45 deletions(-) diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt index 64669224..9937ac16 100644 --- a/Framework/Cascade/CMakeLists.txt +++ b/Framework/Cascade/CMakeLists.txt @@ -42,8 +42,9 @@ target_link_libraries ( CORSIKArandom ProcessSibyll CORSIKAcascade - ProcessStackInspector +# ProcessStackInspector ProcessTrackingLine + ProcessNullModel CORSIKAstackinterface CORSIKAprocesses CORSIKAparticles diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index c55ac0b4..524b99ed 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -16,9 +16,19 @@ #include <corsika/process/ProcessReturn.h> #include <corsika/random/RNGManager.h> #include <corsika/random/UniformRealDistribution.h> -#include <corsika/setup/SetupTrajectory.h> +#include <corsika/stack/SecondaryView.h> #include <corsika/units/PhysicalUnits.h> +#include <corsika/setup/SetupTrajectory.h> + +/* see Issue 161, we need to include SetupStack only because we need + to globally define StackView. This is clearly not nice and should + be changed, when possible. It might be that StackView needs to be + templated in Cascade, but this would be even worse... so we don't + do that until it is really needed. + */ +#include <corsika/setup/SetupStack.h> + #include <cmath> #include <iostream> @@ -35,25 +45,24 @@ namespace corsika::cascade { * it very versatile. Via the template arguments physics models are * plugged into the cascade simulation. * - * <b>Tracking</b> must be a class according to the + * <b>TTracking</b> must be a class according to the * TrackingInterface providing the functions: <code>void * Init();</code> and <code>auto GetTrack(Particle const& p)</auto>, * where the latter has a return type of <code> * geometry::Trajectory<corsika::geometry::Line or Helix> </code> * - * <b>ProcessList</b> must be a ProcessSequence. + * <b>TProcessList</b> must be a ProcessSequence. * TimeOfIntersection(corsika::geometry::Line const& line, * * <b>Stack</b> is the storage object for particle data, i.e. with * Particle class type <code>Stack::ParticleType</code> * * - */ - template <typename Tracking, typename ProcessList, typename Stack> + template <typename TTracking, typename TProcessList, typename TStack> class Cascade { - using Particle = typename Stack::ParticleType; + using Particle = typename TStack::ParticleType; // we only want fully configured objects Cascade() = delete; @@ -63,8 +72,8 @@ namespace corsika::cascade { * Cascade class cannot be default constructed, but needs a valid * list of physics processes for configuration at construct time. */ - Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl, - Stack& stack) + Cascade(corsika::environment::Environment const& env, TTracking& tr, TProcessList& pl, + TStack& stack) : fEnvironment(env) , fTracking(tr) , fProcessSequence(pl) @@ -107,15 +116,16 @@ namespace corsika::cascade { * New particles produced in one step are subject to further * processing, e.g. thinning, etc. */ - void Step(Particle& particle) { + void Step(Particle& vParticle) { + using namespace corsika; using namespace corsika::units::si; // determine geometric tracking - corsika::setup::Trajectory step = fTracking.GetTrack(particle); + setup::Trajectory step = fTracking.GetTrack(vParticle); // determine combined total interaction length (inverse) InverseGrammageType const total_inv_lambda = - fProcessSequence.GetTotalInverseInteractionLength(particle, step); + fProcessSequence.GetTotalInverseInteractionLength(vParticle, step); // sample random exponential step length in grammage std::exponential_distribution expDist(total_inv_lambda * (1_g / (1_m * 1_m))); @@ -126,7 +136,7 @@ namespace corsika::cascade { // convert next_step from grammage to length auto const* currentNode = - fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition()); + fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition()); if (currentNode == &*fEnvironment.GetUniverse()) { throw std::runtime_error("particle entered void universe"); @@ -136,12 +146,12 @@ namespace corsika::cascade { currentNode->GetModelProperties().ArclengthFromGrammage(step, next_interact); // determine the maximum geometric step length - LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step); + LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step); std::cout << "distance_max=" << distance_max << std::endl; // determine combined total inverse decay time InverseTimeType const total_inv_lifetime = - fProcessSequence.GetTotalInverseLifetime(particle); + fProcessSequence.GetTotalInverseLifetime(vParticle); // sample random exponential decay time std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s); @@ -150,9 +160,8 @@ namespace corsika::cascade { << ", next_decay=" << next_decay << std::endl; // convert next_decay from time to length [m] - LengthType const distance_decay = next_decay * particle.GetMomentum().norm() / - particle.GetEnergy() * - corsika::units::constants::c; + LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() / + vParticle.GetEnergy() * units::constants::c; // take minimum of geometry, interaction, decay for next step auto const min_distance = @@ -161,33 +170,45 @@ namespace corsika::cascade { std::cout << " move particle by : " << min_distance << std::endl; // here the particle is actually moved along the trajectory to new position: - // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step); - particle.SetPosition(step.PositionFromArclength(min_distance)); + // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step); + vParticle.SetPosition(step.PositionFromArclength(min_distance)); // .... also update time, momentum, direction, ... + vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c); step.LimitEndTo(min_distance); // particle.GetNode(); // previous VolumeNode - particle.SetNode( + vParticle.SetNode( currentNode); // NOTE @Max : here we need to distinguish: IF particle step is // limited by tracking (via fTracking.GetTrack()), THEN we need // to check/update VolumeNodes. In all other cases it is // guaranteed that we are still in the same volume // apply all continuous processes on particle + track - corsika::process::EProcessReturn status = - fProcessSequence.DoContinuous(particle, step, fStack); + process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step); - if (status == corsika::process::EProcessReturn::eParticleAbsorbed) { - std::cout << "Cascade: delete absorbed particle " << particle.GetPID() << " " - << particle.GetEnergy() / 1_GeV << "GeV" << std::endl; - particle.Delete(); + if (status == process::EProcessReturn::eParticleAbsorbed) { + std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " " + << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl; + vParticle.Delete(); return; } std::cout << "sth. happening before geometric limit ? " << ((min_distance < distance_max) ? "yes" : "no") << std::endl; + /* + Create SecondaryView object on Stack. The data container + remains untouched and identical, and 'projectil' is identical + to 'vParticle' above this line. However, + projectil.AddSecondaries populate the SecondaryView, which can + then be used afterwards for further processing. Thus: it is + important to use projectle (and not vParticle) for Interaction, + and Decay! + */ + setup::StackView secondaries(vParticle); + [[maybe_unused]] auto projectile = secondaries.GetProjectile(); + if (min_distance < distance_max) { // interaction to happen within geometric limit // check whether decay or interaction limits this step @@ -195,33 +216,35 @@ namespace corsika::cascade { std::cout << "collide" << std::endl; InverseGrammageType const actual_inv_length = - fProcessSequence.GetTotalInverseInteractionLength(particle, step); + fProcessSequence.GetTotalInverseInteractionLength(vParticle, step); - corsika::random::UniformRealDistribution<InverseGrammageType> uniDist( - actual_inv_length); + random::UniformRealDistribution<InverseGrammageType> uniDist(actual_inv_length); const auto sample_process = uniDist(fRNG); InverseGrammageType inv_lambda_count = 0. * meter * meter / gram; - fProcessSequence.SelectInteraction(particle, step, fStack, sample_process, + fProcessSequence.SelectInteraction(vParticle, projectile, step, sample_process, inv_lambda_count); } else { std::cout << "decay" << std::endl; InverseTimeType const actual_decay_time = - fProcessSequence.GetTotalInverseLifetime(particle); + fProcessSequence.GetTotalInverseLifetime(vParticle); - corsika::random::UniformRealDistribution<InverseTimeType> uniDist( - actual_decay_time); + random::UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time); const auto sample_process = uniDist(fRNG); InverseTimeType inv_decay_count = 0 / second; - fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count); + fProcessSequence.SelectDecay(vParticle, projectile, sample_process, + inv_decay_count); } + + fProcessSequence.DoSecondaries(secondaries); + vParticle.Delete(); // last thing in Step function } } private: - corsika::environment::Environment const& fEnvironment; - Tracking& fTracking; - ProcessList& fProcessSequence; - Stack& fStack; + environment::Environment const& fEnvironment; + TTracking& fTracking; + TProcessList& fProcessSequence; + TStack& fStack; corsika::random::RNG& fRNG = corsika::random::RNGManager::GetInstance().GetRandomStream("cascade"); }; diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 1e8a7a4d..e88131e7 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -16,6 +16,7 @@ #include <corsika/cascade/Cascade.h> #include <corsika/process/ProcessSequence.h> +#include <corsika/process/null_model/NullModel.h> #include <corsika/process/stack_inspector/StackInspector.h> #include <corsika/process/tracking_line/TrackingLine.h> @@ -75,13 +76,13 @@ public: ProcessSplit(HEPEnergyType e) : fEcrit(e) {} - template <typename Particle, typename T> - LengthType MaxStepLength(Particle&, T&) const { + template <typename Particle, typename Track> + LengthType MaxStepLength(Particle&, Track&) const { return 1_m; } - template <typename Particle, typename T, typename Stack> - EProcessReturn DoContinuous(Particle& p, T&, Stack&) { + template <typename Particle, typename Track> + EProcessReturn DoContinuous(Particle& p, Track&) { fCalls++; HEPEnergyType E = p.GetEnergy(); if (E < fEcrit) { @@ -115,11 +116,12 @@ TEST_CASE("Cascade", "[Cascade]") { auto env = MakeDummyEnv(); tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env); - stack_inspector::StackInspector<setup::Stack> p0(true); + // stack_inspector::StackInspector<setup::Stack> stackInspect(true); + null_model::NullModel nullModel; const HEPEnergyType Ecrit = 85_MeV; ProcessSplit p1(Ecrit); - auto sequence = p0 << p1; + auto sequence = nullModel << p1; setup::Stack stack; cascade::Cascade EAS(env, tracking, sequence, stack); -- GitLab