diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 0caa88ae6fd495258c0c1c473f2e29cf685c0cc4..ee8e7db37d5d01af89f70e72e9e20ac2bdb97da3 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -133,6 +133,21 @@ namespace corsika::cascade { } } + /** + * Force an interaction of the top particle of the stack at its current position. + * Note that SetNodes() or an equivalent procedure needs to be called first if you + * want to call forceInteraction() for the primary interaction. + */ + void forceInteraction() { + std::cout << "forced interaction!" << std::endl; + auto vParticle = fStack.GetNextParticle(); + TStackView secondaries(vParticle); + auto projectile = secondaries.GetProjectile(); + interaction(vParticle, projectile); + fProcessSequence.DoSecondaries(secondaries); + vParticle.Delete(); // todo: this should be reviewed, see below + } + private: /** * The Step function is executed for each particle from the @@ -243,28 +258,10 @@ namespace corsika::cascade { [[maybe_unused]] auto projectile = secondaries.GetProjectile(); if (min_distance == distance_interact) { - std::cout << "collide" << std::endl; - - InverseGrammageType const current_inv_length = - fProcessSequence.GetTotalInverseInteractionLength(vParticle); - - random::UniformRealDistribution<InverseGrammageType> uniDist( - current_inv_length); - const auto sample_process = uniDist(fRNG); - InverseGrammageType inv_lambda_count = 0. * meter * meter / gram; - fProcessSequence.SelectInteraction(vParticle, projectile, sample_process, - inv_lambda_count); + interaction(vParticle, projectile); } else { assert(min_distance == distance_decay); - std::cout << "decay" << std::endl; - InverseTimeType const actual_decay_time = - fProcessSequence.GetTotalInverseLifetime(vParticle); - - random::UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time); - const auto sample_process = uniDist(fRNG); - InverseTimeType inv_decay_count = 0 / second; - fProcessSequence.SelectDecay(vParticle, projectile, sample_process, - inv_decay_count); + decay(vParticle, projectile); // make sure particle actually did decay if it should have done so if (secondaries.GetSize() == 1 && projectile.GetPID() == secondaries.GetNextParticle().GetPID()) @@ -299,6 +296,35 @@ namespace corsika::cascade { } } + auto decay(Particle& particle, + decltype(std::declval<TStackView>().GetProjectile()) projectile) { + std::cout << "decay" << std::endl; + units::si::InverseTimeType const actual_decay_time = + fProcessSequence.GetTotalInverseLifetime(particle); + + 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(); + return fProcessSequence.SelectDecay(particle, projectile, sample_process, + inv_decay_count); + } + + auto interaction(Particle& particle, + decltype(std::declval<TStackView>().GetProjectile()) projectile) { + std::cout << "collide" << std::endl; + + units::si::InverseGrammageType const current_inv_length = + fProcessSequence.GetTotalInverseInteractionLength(particle); + + random::UniformRealDistribution<units::si::InverseGrammageType> uniDist( + current_inv_length); + const auto sample_process = uniDist(fRNG); + auto inv_lambda_count = units::si::InverseGrammageType::zero(); + return fProcessSequence.SelectInteraction(particle, projectile, sample_process, + inv_lambda_count); + } + private: corsika::environment::Environment<MediumInterface> const& fEnvironment; TTracking& fTracking; diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 7992c32afcd348ef890dcaf98c27a23fa104e6fc..05e8898b6c214413c0591ed00ea7cfc0f18c09a7 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -136,6 +136,7 @@ TEST_CASE("Cascade", "[Cascade]") { rmng.RegisterRandomStream("cascade"); auto env = MakeDummyEnv(); + auto const& rootCS = env.GetCoordinateSystem(); tracking_line::TrackingLine tracking; stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0); @@ -147,13 +148,6 @@ TEST_CASE("Cascade", "[Cascade]") { ProcessCut cut(Ecrit); auto sequence = nullModel << stackInspect << split << cut; TestCascadeStack stack; - - cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack, - TestCascadeStackView> - EAS(env, tracking, sequence, stack); - CoordinateSystem const& rootCS = - RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); - stack.Clear(); stack.AddParticle( std::tuple<particles::Code, units::si::HEPEnergyType, @@ -161,10 +155,25 @@ TEST_CASE("Cascade", "[Cascade]") { particles::Code::Electron, E0, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}), Point(rootCS, {0_m, 0_m, 10_km}), 0_ns}); + + cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack, + TestCascadeStackView> + EAS(env, tracking, sequence, stack); + EAS.Init(); - EAS.Run(); - CHECK(cut.GetCount() == 2048); - CHECK(cut.GetCalls() == 2047); - CHECK(split.GetCalls() == 2047); + SECTION("full cascade") { + EAS.Run(); + + CHECK(cut.GetCount() == 2048); + CHECK(cut.GetCalls() == 2047); + CHECK(split.GetCalls() == 2047); + } + + SECTION("forced interaction") { + EAS.SetNodes(); + EAS.forceInteraction(); + CHECK(stack.GetSize() == 2); + CHECK(split.GetCalls() == 1); + } }