diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 2ef7cdcf0d03a973e5e57099d0fa11f48eabd58f..54910cf9c890cbfd5d95d62b0b8fcbe64d11c157 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -247,28 +247,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"); - 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; - // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); @@ -303,6 +281,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 e5b4b3456c42c8e03f28847e41d0babdee9ab464..9bc0636d242da9990b55f1ada2bf0fbe6808f794 100644 --- a/Documentation/Examples/cascade_proton_example.cc +++ b/Documentation/Examples/cascade_proton_example.cc @@ -12,6 +12,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> @@ -237,35 +238,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); - 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; - - // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() - // << "\n"; - // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); @@ -297,6 +269,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 d029a35cf6c9eb0d1af938c5522d2ee68a4bf75a..e0d5d8d085cb0173778cff86c0af490cfaa99944 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -13,6 +13,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> @@ -241,35 +242,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); - // 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; - - // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() - // << "\n"; - // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); @@ -304,6 +276,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 89fd65ef42150149a9c0ffbd7183fed865ba659d..4813d64432dff8bfdf0eaa98d2cf0154f6ac2af9 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -132,13 +132,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); @@ -155,7 +158,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 af55e7efd0f633122efe9d008d6d5829f4677117..0972fc591a0fc003d8eb9962a429fd3a8bbdad9e 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -141,8 +141,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 1458fa75ca5938cdd42d50d0b4d22af46f641302..7c81405918e0aaccc94a2e7e8b4bccccd94ce994 100644 --- a/Framework/ProcessSequence/StackProcess.h +++ b/Framework/ProcessSequence/StackProcess.h @@ -43,6 +43,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 2835a7043cc80101ef22d1301af981b5735de6bb..df9d2e64afb6694c117e24741e327d6db837e7cb 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -1,17 +1,24 @@ - +# 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) +# continuous physics +add_subdirectory (EnergyLoss) +add_subdirectory (TrackWriter) +# stack processes +add_subdirectory (StackInspector) +# secondaries process +# cuts, thinning, etc. -#add_custom_target(CORSIKAprocesses) +########################################## +# 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 e56503378ed5abfc15578c474ed0afd4d54d6412..9fbdc1cce885745262c3ad3be4df68637af27e5c 100644 --- a/Processes/StackInspector/StackInspector.cc +++ b/Processes/StackInspector/StackInspector.cc @@ -18,6 +18,8 @@ #include <corsika/setup/SetupTrajectory.h> +#include <chrono> +#include <iomanip> #include <iostream> #include <limits> using namespace std; @@ -28,41 +30,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 62f36eeebae774f74d26953697707e698c75da6e..761a7c26d01392b29e70f66382090de0d3fae571 100644 --- a/Processes/StackInspector/StackInspector.h +++ b/Processes/StackInspector/StackInspector.h @@ -16,6 +16,8 @@ #include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> +#include <chrono> + namespace corsika::process { namespace stack_inspector { @@ -25,16 +27,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 832bea00a8074214ef9f3a235a6a574f7286cf02..b575aa24c45901bd004ba75ab2df308f26fe4372 100644 --- a/Processes/StackInspector/testStackInspector.cc +++ b/Processes/StackInspector/testStackInspector.cc @@ -50,7 +50,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);