diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 8c29a77987316d40f1aebf6df392f1abf581e446..6fe5e1700e603b682bf9ab1152ead901eeda8d4a 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -27,8 +27,20 @@ namespace corsika::cascade { void Run() { while (!fStack.IsEmpty()) { while (!fStack.IsEmpty()) { - Particle& p = *fStack.GetNextParticle(); - Step(p); + //Particle& p = *fStack.GetNextParticle(); + EnergyType Emin; + typename Stack::StackIterator pMin(fStack, 0); + bool first = true; + for (typename Stack::StackIterator ip = fStack.begin(); ip!=fStack.end(); ++ip) + { + if (first || ip.GetEnergy()<Emin) { + first = false; + pMin = ip; + Emin = pMin.GetEnergy(); + } + } + + Step(pMin); } // do cascade equations, which can put new particles on Stack, // thus, the double loop diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 3f9f1b0fcc61ae6764adb13424fc447dfc2a5d7e..7e0bd76ba6e94655dc8c3a4c83c06c1acd87c7a5 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -13,6 +13,8 @@ using namespace corsika::units; #include <iostream> using namespace std; +static int fCount = 0; + class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> { public: ProcessSplit() {} @@ -28,13 +30,13 @@ public: template <typename Particle, typename Stack> void DoDiscrete(Particle& p, Stack& s) const { EnergyType E = p.GetEnergy(); - if (E < 1_GeV) { + if (E < 85_MeV) { p.Delete(); fCount++; } else { p.SetEnergy(E / 2); s.NewParticle().SetEnergy(E / 2); - } + } } void Init() { fCount = 0; } @@ -42,7 +44,6 @@ public: int GetCount() { return fCount; } private: - mutable int fCount = 0; }; class ProcessReport : public corsika::process::BaseProcess<ProcessReport> { @@ -59,17 +60,22 @@ public: template <typename Particle, typename Trajectory, typename Stack> void DoContinuous(Particle& p, Trajectory& t, Stack& s) const { + static int countStep = 0; if (!fReport) return; - static int fCount = 0; - std::cout << "generation " << fCount << std::endl; + //std::cout << "generation " << countStep << std::endl; int i = 0; + EnergyType Etot = 0_GeV; for (auto& iterP : s) { EnergyType E = iterP.GetEnergy(); - std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID() + Etot += E; + /* std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID() << ", E=" << double(E / 1_GeV) << " GeV " << " | " << std::endl; + */ } - fCount++; + countStep++; + //cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl; + cout << countStep << " " << s.GetSize() << " " << Etot/1_GeV << " " << fCount << endl; } template <typename Particle, typename Stack> @@ -79,25 +85,32 @@ public: TEST_CASE("Cascade", "[Cascade]") { - ProcessReport p0(false); + ProcessReport p0(true); ProcessSplit p1; const auto sequence = p0 + p1; corsika::stack::super_stupid::SuperStupidStack stack; - + corsika::cascade::Cascade<corsika::geometry::LineTrajectory, decltype(sequence), decltype(stack)> EAS(sequence, stack); + stack.Clear(); + auto particle = stack.NewParticle(); + EnergyType E0 = 100_GeV; + particle.SetEnergy(E0); + EAS.Init(); + EAS.Run(); + SECTION("sectionTwo") { - for (int i = 0; i < 5; ++i) { + for (int i = 0; i < 0; ++i) { stack.Clear(); auto particle = stack.NewParticle(); EnergyType E0 = 100_GeV * pow(10, i); particle.SetEnergy(E0); EAS.Init(); EAS.Run(); - - cout << "E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl; + + //cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl; } } } diff --git a/Framework/StackInterface/StackIterator.h b/Framework/StackInterface/StackIterator.h index 8137e9b53b8268679cdcb936c477cba914cb9c8e..a6ca8bed67c4d25fbbd97abab1bda1bb3157e9ec 100644 --- a/Framework/StackInterface/StackIterator.h +++ b/Framework/StackInterface/StackIterator.h @@ -63,9 +63,12 @@ namespace corsika::stack { StackIteratorInterface(const StackIteratorInterface& mit) : fData(mit.fData) , fIndex(mit.fIndex) {} + + public: StackIteratorInterface& operator=(const StackIteratorInterface& mit) { fData = mit.fData; fIndex = mit.fIndex; + return *this; } public: