diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 936ed40001c63b4c948883cf62360dc15f4f4c10..3a65b69183e5c8ecf94bb7850d57ce912541d5da 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -93,27 +93,6 @@ int main() { universe.AddChild(std::move(outerMedium)); - // setup processes, decays and interactions - tracking_line::TrackingLine tracking; - stack_inspector::StackInspector<setup::Stack> stackInspect(1, true); - - const std::vector<particles::Code> trackedHadrons = { - particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus, - particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short}; - - random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); - process::sibyll::Interaction sibyll; - process::sibyll::NuclearInteraction sibyllNuc(sibyll, env); - process::sibyll::Decay decay(trackedHadrons); - process::particle_cut::ParticleCut cut(20_GeV); - - process::track_writer::TrackWriter trackWriter("tracks.dat"); - process::energy_loss::EnergyLoss eLoss; - - // assemble all processes into an ordered process list - auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut - << trackWriter; - // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); @@ -148,6 +127,28 @@ int main() { beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ}); } + // setup processes, decays and interactions + tracking_line::TrackingLine tracking; + stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0); + + const std::vector<particles::Code> trackedHadrons = { + particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus, + particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short}; + + random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); + random::RNGManager::GetInstance().RegisterRandomStream("pythia"); + process::sibyll::Interaction sibyll; + process::sibyll::NuclearInteraction sibyllNuc(sibyll, env); + process::sibyll::Decay decay(trackedHadrons); + ProcessCut cut(20_GeV); + + process::track_writer::TrackWriter trackWriter("tracks.dat"); + process::energy_loss::EnergyLoss eLoss; + + // assemble all processes into an ordered process list + auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut + << trackWriter; + // define air shower object, run simulation cascade::Cascade EAS(env, tracking, sequence, stack); EAS.Init(); diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc index c917d9b4360c6663fd7c4ee6e63b8fed73217e93..6e84f5614e22431aee8864d003048901fffa9db3 100644 --- a/Documentation/Examples/cascade_proton_example.cc +++ b/Documentation/Examples/cascade_proton_example.cc @@ -11,6 +11,7 @@ #include <corsika/cascade/Cascade.h> #include <corsika/process/ProcessSequence.h> #include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h> +#include <corsika/process/stack_inspector/StackInspector.h> #include <corsika/process/tracking_line/TrackingLine.h> #include <corsika/setup/SetupStack.h> @@ -85,35 +86,6 @@ int main() { const CoordinateSystem& rootCS = env.GetCoordinateSystem(); - // setup processes, decays and interactions - tracking_line::TrackingLine tracking; - - const std::vector<particles::Code> trackedHadrons = { - particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus, - particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short}; - - random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); - random::RNGManager::GetInstance().RegisterRandomStream("pythia"); - // process::sibyll::Interaction sibyll(env); - process::pythia::Interaction pythia; - // process::sibyll::NuclearInteraction sibyllNuc(env, sibyll); - // process::sibyll::Decay decay(trackedHadrons); - process::pythia::Decay decay(trackedHadrons); - process::particle_cut::ParticleCut cut(20_GeV); - - // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); - // process::HadronicElasticModel::HadronicElasticInteraction - // hadronicElastic(env); - - process::track_writer::TrackWriter trackWriter("tracks.dat"); - - // assemble all processes into an ordered process list - // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter; - auto sequence = pythia << decay << cut << trackWriter; - - // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() - // << "\n"; - // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); @@ -145,6 +117,36 @@ int main() { beamCode, E0, plab, pos, 0_ns}); } + // setup processes, decays and interactions + tracking_line::TrackingLine tracking; + stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0); + + const std::vector<particles::Code> trackedHadrons = { + particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus, + particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short}; + + random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); + random::RNGManager::GetInstance().RegisterRandomStream("pythia"); + // process::sibyll::Interaction sibyll(env); + process::pythia::Interaction pythia; + // process::sibyll::NuclearInteraction sibyllNuc(env, sibyll); + // process::sibyll::Decay decay(trackedHadrons); + process::pythia::Decay decay(trackedHadrons); + ProcessCut cut(20_GeV); + + // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); + // process::HadronicElasticModel::HadronicElasticInteraction + // hadronicElastic(env); + + process::track_writer::TrackWriter trackWriter("tracks.dat"); + + // assemble all processes into an ordered process list + // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter; + auto sequence = pythia << decay << cut << trackWriter << stackInspect; + + // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() + // << "\n"; + // define air shower object, run simulation cascade::Cascade EAS(env, tracking, sequence, stack); EAS.Init(); diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index a133be5df380e8171fd5599e68347c76984d5ac7..fba55699b2f21618cd1e40f1d4a8146b14074900 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -12,6 +12,7 @@ #include <corsika/process/ProcessSequence.h> #include <corsika/process/energy_loss/EnergyLoss.h> #include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h> +#include <corsika/process/stack_inspector/StackInspector.h> #include <corsika/process/tracking_line/TrackingLine.h> #include <corsika/setup/SetupStack.h> @@ -86,32 +87,6 @@ int main() { const CoordinateSystem& rootCS = env.GetCoordinateSystem(); - // setup processes, decays and interactions - tracking_line::TrackingLine tracking; - - const std::vector<particles::Code> trackedHadrons = { - particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus, - particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short}; - - random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); - process::sibyll::Interaction sibyll; - process::sibyll::NuclearInteraction sibyllNuc(sibyll, env); - process::sibyll::Decay decay(trackedHadrons); - process::particle_cut::ParticleCut cut(20_GeV); - // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); - // process::HadronicElasticModel::HadronicElasticInteraction - // hadronicElastic(env); - - process::track_writer::TrackWriter trackWriter("tracks.dat"); - process::energy_loss::EnergyLoss eLoss; - - // assemble all processes into an ordered process list - // cut << trackWriter; - auto sequence = sibyll << sibyllNuc << decay << eLoss << cut; - - // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() - // << "\n"; - // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); @@ -146,6 +121,36 @@ int main() { beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ}); } + // setup processes, decays and interactions + tracking_line::TrackingLine tracking; + stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0); + + const std::vector<particles::Code> trackedHadrons = { + particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus, + particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short}; + + random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); + process::sibyll::Interaction sibyll; + process::sibyll::NuclearInteraction sibyllNuc(sibyll, env); + process::sibyll::Decay decay(trackedHadrons); + // random::RNGManager::GetInstance().RegisterRandomStream("pythia"); + // process::pythia::Decay decay(trackedHadrons); + ProcessCut cut(20_GeV); + + // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); + // process::HadronicElasticModel::HadronicElasticInteraction + // hadronicElastic(env); + + process::track_writer::TrackWriter trackWriter("tracks.dat"); + process::energy_loss::EnergyLoss eLoss; + + // assemble all processes into an ordered process list + // cut << trackWriter; + auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << stackInspect; + + // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() + // << "\n"; + // define air shower object, run simulation cascade::Cascade EAS(env, tracking, sequence, stack); EAS.Init(); diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 9c06c332e61eb25e934ef90713ec5fd179fa622c..fcabaabcfc672df70a45cd20239b4546331c4922 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -131,13 +131,16 @@ public: }; TEST_CASE("Cascade", "[Cascade]") { + + HEPEnergyType E0 = 100_GeV; + random::RNGManager& rmng = random::RNGManager::GetInstance(); rmng.RegisterRandomStream("cascade"); auto env = MakeDummyEnv(); tracking_line::TrackingLine tracking; - stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true); + stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0); null_model::NullModel nullModel; const GrammageType X0 = 20_g / square(1_cm); @@ -154,7 +157,6 @@ TEST_CASE("Cascade", "[Cascade]") { RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); stack.Clear(); - HEPEnergyType E0 = 100_GeV; stack.AddParticle( std::tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index 333a4f47d9678a0b7ee2d2bd61437240b71ec8fb..d45705900d2c3cbba97d4dc2fab9e5595f0e87b5 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -156,8 +156,8 @@ namespace corsika::process { } /** - Execute the StackProcess-es in the ProcessSequence - */ + Execute the StackProcess-es in the ProcessSequence + */ template <typename TStack> EProcessReturn DoStack(TStack& vS) { EProcessReturn ret = EProcessReturn::eOk; diff --git a/Framework/ProcessSequence/StackProcess.h b/Framework/ProcessSequence/StackProcess.h index 56cf4991c56add045b0a419262b0ec0a021a0da6..2b2c9fcd0effcb16bf59b3d51708f8e1991edc8f 100644 --- a/Framework/ProcessSequence/StackProcess.h +++ b/Framework/ProcessSequence/StackProcess.h @@ -42,6 +42,7 @@ namespace corsika::process { template <typename TStack> inline EProcessReturn DoStack(TStack&); + int GetStep() const { return fIStep; } bool CheckStep() { return !((++fIStep) % fNStep); } private: diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index e3057e2db42052314aea04b8e8c02f2aea963709..7ef912d1556dfb32f2e721abbc3cbe324bc33a37 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -1,17 +1,27 @@ +# general add_subdirectory (NullModel) +# tracking +add_subdirectory (TrackingLine) +# hadron interaction models add_subdirectory (Sibyll) if (PYTHIA8_FOUND) add_subdirectory (Pythia) endif (PYTHIA8_FOUND) -add_subdirectory (StackInspector) -add_subdirectory (TrackWriter) -add_subdirectory (TrackingLine) add_subdirectory (HadronicElasticModel) -add_subdirectory (EnergyLoss) add_subdirectory (UrQMD) add_subdirectory (SwitchProcess) add_subdirectory (ParticleCut) +# continuous physics +add_subdirectory (EnergyLoss) +add_subdirectory (TrackWriter) +# stack processes +add_subdirectory (StackInspector) +# secondaries process +# cuts, thinning, etc. + +########################################## +# add_custom_target(CORSIKAprocesses) add_library (CORSIKAprocesses INTERFACE) add_dependencies(CORSIKAprocesses ProcessNullModel) add_dependencies(CORSIKAprocesses ProcessSibyll) diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc index f9ad04daa9ce083c84c69c971b676a1c57298923..b982cdb7275d2526dbb4a8b91ff0dd952a3626f0 100644 --- a/Processes/StackInspector/StackInspector.cc +++ b/Processes/StackInspector/StackInspector.cc @@ -17,6 +17,8 @@ #include <corsika/setup/SetupTrajectory.h> +#include <chrono> +#include <iomanip> #include <iostream> #include <limits> using namespace std; @@ -27,41 +29,58 @@ using namespace corsika::units::si; using namespace corsika::process::stack_inspector; template <typename TStack> -StackInspector<TStack>::StackInspector(const int nStep, const bool aReport) - : StackProcess<StackInspector<TStack>>(nStep) - , fReport(aReport) - , fCountStep(0) {} +StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack, + const HEPEnergyType vE0) + : StackProcess<StackInspector<TStack>>(vNStep) + , fReportStack(vReportStack) + , fE0(vE0) + , fStartTime(std::chrono::system_clock::now()) {} template <typename TStack> StackInspector<TStack>::~StackInspector() {} template <typename TStack> -process::EProcessReturn StackInspector<TStack>::DoStack(TStack const& vS) { - if (!fReport) return process::EProcessReturn::eOk; +process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) { [[maybe_unused]] int i = 0; HEPEnergyType Etot = 0_GeV; - for (auto& iterP : vS) { + for (const auto& iterP : vS) { HEPEnergyType E = iterP.GetEnergy(); Etot += E; - geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance() - .GetRootCoordinateSystem(); // for printout - auto pos = iterP.GetPosition().GetCoordinates(rootCS); - cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30) - << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, " - << " pos=" << pos << " node = " << iterP.GetNode(); - if (iterP.GetPID() == Code::Nucleus) cout << " nuc_ref=" << iterP.GetNucleusRef(); - cout << endl; + if (fReportStack) { + geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance() + .GetRootCoordinateSystem(); // for printout + auto pos = iterP.GetPosition().GetCoordinates(rootCS); + cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30) + << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, " + << " pos=" << pos << " node = " << iterP.GetNode(); + if (iterP.GetPID() == Code::Nucleus) cout << " nuc_ref=" << iterP.GetNucleusRef(); + cout << endl; + } } - fCountStep++; - cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << vS.GetSize() - << " Estack=" << Etot / 1_GeV << " GeV" << endl; + + auto const now = std::chrono::system_clock::now(); + const std::chrono::duration<double> elapsed_seconds = now - fStartTime; + std::time_t const now_time = std::chrono::system_clock::to_time_t(now); + double const progress = (fE0 - Etot) / fE0; + double const eta_seconds = elapsed_seconds.count() / progress; + std::time_t const eta_time = std::chrono::system_clock::to_time_t( + fStartTime + std::chrono::seconds((int)eta_seconds)); + + cout << "StackInspector: " + << " time=" << std::put_time(std::localtime(&now_time), "%T") + << ", running=" << elapsed_seconds.count() << " seconds" + << " (" << setw(3) << int(progress * 100) << "%)" + << ", nStep=" << GetStep() << ", stackSize=" << vS.GetSize() + << ", Estack=" << Etot / 1_GeV << " GeV" + << ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << endl; return process::EProcessReturn::eOk; } template <typename TStack> void StackInspector<TStack>::Init() { - fCountStep = 0; + fReportStack = false; + fStartTime = std::chrono::system_clock::now(); } #include <corsika/cascade/testCascade.h> diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h index 24dbf3b9bc5127dc3497499acb7d4f4562e02476..28d6f04a12ce271e8920f88100848ffdbb93f727 100644 --- a/Processes/StackInspector/StackInspector.h +++ b/Processes/StackInspector/StackInspector.h @@ -15,6 +15,8 @@ #include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> +#include <chrono> + namespace corsika::process { namespace stack_inspector { @@ -24,16 +26,25 @@ namespace corsika::process { typedef typename TStack::ParticleType Particle; + using corsika::process::StackProcess<StackInspector<TStack>>::GetStep; + public: - StackInspector(const int nStep, const bool aReport); + StackInspector(const int vNStep, const bool vReportStack, + const corsika::units::si::HEPEnergyType vE0); ~StackInspector(); void Init(); - EProcessReturn DoStack(TStack const&); + EProcessReturn DoStack(const TStack&); + + /** + * To set a new E0, for example when a new shower event is started + */ + void SetE0(const corsika::units::si::HEPEnergyType vE0) { fE0 = vE0; } private: - bool fReport; - int fCountStep = 0; + bool fReportStack; + corsika::units::si::HEPEnergyType fE0; + decltype(std::chrono::system_clock::now()) fStartTime; }; } // namespace stack_inspector diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc index ba738ef8e39ce52e3150453526e5754ec95dda5f..59662d491ef24a4b871c041ccc903fd6cb4b759f 100644 --- a/Processes/StackInspector/testStackInspector.cc +++ b/Processes/StackInspector/testStackInspector.cc @@ -49,7 +49,7 @@ TEST_CASE("StackInspector", "[processes]") { SECTION("interface") { - StackInspector<TestCascadeStack> model(1, true); + StackInspector<TestCascadeStack> model(1, true, E0); model.Init(); [[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);