testInteractionCounter.cpp

/*
 * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
 *
 * This software is distributed under the terms of the GNU General Public
 * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
 * the license.
 */

#include <corsika/modules/InteractionCounter.hpp>

#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>

#include <corsika/setup/SetupStack.hpp>

#include <catch2/catch.hpp>

#include <numeric>

using namespace corsika;

const std::string refDataDir = std::string(REFDATADIR); // from cmake

struct DummyProcess {
  template <typename TParticle>
  GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
    return 100_g / 1_cm / 1_cm;
  }

  template <typename TParticle>
  auto DoInteraction([[maybe_unused]] TParticle& projectile) {
    return nullptr;
  }
};

TEST_CASE("InteractionCounter", "[process]") {

  logging::set_level(logging::level::info);
  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");

  DummyProcess d;
  InteractionCounter countedProcess(d);

  SECTION("GetInteractionLength") {
    CHECK(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
  }

  auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
  [[maybe_unused]] auto& env_dummy = env;

  SECTION("DoInteraction nucleus") {
    unsigned short constexpr A = 14, Z = 7;
    auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
                                                             Z, 105_TeV, nodePtr, *csPtr);
    CHECK(stackPtr->getEntries() == 1);
    CHECK(secViewPtr->getEntries() == 0);

    auto const ret = countedProcess.DoInteraction(*secViewPtr);
    CHECK(ret == nullptr);

    auto const& h = countedProcess.GetHistogram().labHist();
    CHECK(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
    CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);

    auto const& h2 = countedProcess.GetHistogram().CMSHist();
    CHECK(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
    CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);

    countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
                                          utl::SaveMode::overwrite);
    countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz",
                                          utl::SaveMode::overwrite);
  }

  SECTION("DoInteraction Lambda") {
    auto constexpr code = particles::Code::Lambda0;
    auto [stackPtr, secViewPtr] =
        setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
    CHECK(stackPtr->getEntries() == 1);
    CHECK(secViewPtr->getEntries() == 0);

    auto const ret = countedProcess.DoInteraction(*secViewPtr);
    CHECK(ret == nullptr);

    auto const& h = countedProcess.GetHistogram().labHist();
    CHECK(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
    CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);

    auto const& h2 = countedProcess.GetHistogram().CMSHist();
    CHECK(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
    CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
  }
}

#include <algorithm>
#include <iterator>
#include <string>
#include <fstream>

TEST_CASE("InteractionCounterOutput", "[output validation]") {

  logging::set_level(logging::level::info);
  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
  
  auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
                       "testInteractionCounter_file2");

  SECTION(std::string("check saved data, ") + file + ".npz") {

    std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;

    // compare to binary reference data
    std::ifstream file1(file + ".npz");
    std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");

    std::istreambuf_iterator<char> begin1(file1);
    std::istreambuf_iterator<char> begin1ref(file1ref);

    std::istreambuf_iterator<char> end;

    while (begin1 != end && begin1ref != end) {
      CHECK(*begin1 == *begin1ref);
      ++begin1;
      ++begin1ref;
    }
    CHECK(begin1 == end);
    CHECK(begin1ref == end);
    file1.close();
    file1ref.close();
  }
}