diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp index 450a50776ac2ba568857fb316d4e1b1333a680f5..159baed4f73255038b226236a674f0c916500570 100644 --- a/examples/hybrid_MC.cpp +++ b/examples/hybrid_MC.cpp @@ -82,6 +82,34 @@ void registerRandomStreams(uint64_t const seed) { RNGManager<>::getInstance().setSeed(seed); } +class TrackCheck : public ContinuousProcess<TrackCheck> { + +public: + TrackCheck(Plane const& plane) + : plane_(plane) {} + + template <typename TParticle, typename TTrack> + ProcessReturn doContinuous(TParticle const& particle, TTrack const&, bool const) { + auto const delta = plane_.getCenter() - particle.getPosition(); + auto const n = plane_.getNormal(); + auto const proj = n.dot(delta); + if (proj < -0_mm) { + CORSIKA_LOG_INFO("particle {} failes", particle.asString()); + throw std::runtime_error("particle below obs level"); + } + return ProcessReturn::Ok; + } + + template <typename TParticle, typename TTrack> + LengthType getMaxStepLength(TParticle const&, TTrack const&) const { + return std::numeric_limits<double>::infinity() * 1_m; + } + +private: + Plane plane_; +}; + + template <typename T> using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; @@ -151,7 +179,7 @@ int main(int argc, char** argv) { << ", norm = " << plab.getNorm() << endl; auto const observationHeight = 0_km + builder.getEarthRadius(); - auto const injectionHeight = 112.75_km + builder.getEarthRadius(); + auto const injectionHeight = 1_km + builder.getEarthRadius(); auto const t = -observationHeight * cos(thetaRad) + sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) + static_pow<2>(injectionHeight)); @@ -209,7 +237,7 @@ int main(int argc, char** argv) { decaySibyll.printDecayConfig(); - ParticleCut cut{60_GeV, false, true}; + ParticleCut cut{3_GeV, false, true}; BetheBlochPDG eLoss{showerAxis}; CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0, @@ -239,16 +267,16 @@ int main(int argc, char** argv) { auto hadronSequence = make_select(EnergySwitch(55_GeV), urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted)); auto decaySequence = make_sequence(decayPythia, decaySibyll); - auto sequence = make_sequence(hadronSequence, reset_particle_mass, decaySequence, eLoss, - cut, conex_model, longprof, observationLevel); + auto sequence = + make_sequence(TrackCheck(obsPlane), hadronSequence, reset_particle_mass, + decaySequence, eLoss, cut, conex_model, longprof, observationLevel); // define air shower object, run simulation setup::Tracking tracking; Cascade EAS(env, tracking, sequence, output, stack); // to fix the point of first interaction, uncomment the following two lines: - // EAS.SetNodes(); - // EAS.forceInteraction(); + EAS.forceInteraction(); output.startOfShower(); EAS.run();