testCascade.cc

#include <corsika/cascade/Cascade.h>
#include <corsika/geometry/LineTrajectory.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/stack/super_stupid/SuperStupidStack.h>

#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
                          // cpp file
#include <catch2/catch.hpp>

using namespace corsika::process;
using namespace corsika::units;

#include <iostream>
using namespace std;

static int fCount = 0;

class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
  ProcessSplit() {}

  template <typename Particle>
  double MinStepLength(Particle&) const {
    return 0;
  }

  template <typename Particle, typename Trajectory, typename Stack>
  void DoContinuous(Particle& p, Trajectory& t, Stack& s) const {}

  template <typename Particle, typename Stack>
  void DoDiscrete(Particle& p, Stack& s) const {
    EnergyType E = p.GetEnergy();
    if (E < 85_MeV) {
      p.Delete();
      fCount++;
    } else {
      p.SetEnergy(E / 2);
      s.NewParticle().SetEnergy(E / 2);
    }    
  }

  void Init() { fCount = 0; }

  int GetCount() { return fCount; }

private:
};

class ProcessReport : public corsika::process::BaseProcess<ProcessReport> {
  bool fReport = false;

public:
  ProcessReport(bool v)
      : fReport(v) {}

  template <typename Particle>
  double MinStepLength(Particle&) const {
    return 0;
  }

  template <typename Particle, typename Trajectory, typename Stack>
  void DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
    static int countStep = 0;
    if (!fReport) return;
    //std::cout << "generation  " << countStep << std::endl;
    int i = 0;
    EnergyType Etot = 0_GeV;
    for (auto& iterP : s) {
      EnergyType E = iterP.GetEnergy();
      Etot += E;
      /*      std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()
                << ", E=" << double(E / 1_GeV) << " GeV "
                << " | " << std::endl;
      */
    }
    countStep++;
    //cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl;
    cout << countStep << " " << s.GetSize() << " " << Etot/1_GeV << " " << fCount << endl;
  }

  template <typename Particle, typename Stack>
  void DoDiscrete(Particle& p, Stack& s) const {}
  void Init() {}
};

TEST_CASE("Cascade", "[Cascade]") {

  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 < 0; ++i) {
      stack.Clear();
      auto particle = stack.NewParticle();
      EnergyType E0 = 100_GeV * pow(10, i);
      particle.SetEnergy(E0);
      EAS.Init();
      EAS.Run();
       
      //cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
    }
  }
}