Interaction.h

/**
 * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
 *
 * See file AUTHORS for a list of contributors.
 *
 * 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.
 */

#ifndef _corsika_process_sibyll_interaction_h_
#define _corsika_process_sibyll_interaction_h_

#include <corsika/process/InteractionProcess.h>

#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>

namespace corsika::process::sibyll {

  class Interaction : public corsika::process::InteractionProcess<Interaction> {

    int fCount = 0;
    int fNucCount = 0;

  public:
    Interaction(corsika::environment::Environment const& env)
        : fEnvironment(env) {}
    ~Interaction() {
      std::cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount
                << std::endl;
    }

    void Init() {

      using corsika::random::RNGManager;
      using std::cout;
      using std::endl;

      // initialize Sibyll
      sibyll_ini_();
    }

    std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
        const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
        const corsika::units::si::HEPEnergyType CoMenergy) {
      using namespace corsika::units::si;
      double sigProd, dummy, dum1, dum2, dum3, dum4;
      double dumdif[3];
      const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
      const double dEcm = CoMenergy / 1_GeV;
      if (corsika::particles::IsNucleus(TargetId)) {
        const int iTarget = corsika::particles::GetNucleusA(TargetId);
        if (iTarget > 18 || iTarget == 0)
          throw std::runtime_error(
              "Sibyll target outside range. Only nuclei with A<18 are allowed.");
        sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
        return std::make_tuple(sigProd * 1_mbarn, iTarget);
      } else if (TargetId == corsika::particles::Proton::GetCode()) {
        sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
        return std::make_tuple(sigProd * 1_mbarn, 1);
      } else {
        // no interaction in sibyll possible, return infinite cross section? or throw?
        sigProd = std::numeric_limits<double>::infinity();
        return std::make_tuple(sigProd * 1_mbarn, 0);
      }
    }

    template <typename Particle, typename Track>
    corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {

      using namespace corsika::units;
      using namespace corsika::units::si;
      using namespace corsika::geometry;
      using std::cout;
      using std::endl;

      // coordinate system, get global frame of reference
      CoordinateSystem& rootCS =
          RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();

      const particles::Code corsikaBeamId = p.GetPID();

      // beam particles for sibyll : 1, 2, 3 for p, pi, k
      // read from cross section code table
      const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);

      const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
                                              corsika::particles::Neutron::GetMass());

      // FOR NOW: assume target is at rest
      MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});

      // total momentum and energy
      HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
      MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
      pTotLab += p.GetMomentum();
      pTotLab += pTarget;
      auto const pTotLabNorm = pTotLab.norm();
      // calculate cm. energy
      const HEPEnergyType ECoM = sqrt(
          (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy

      std::cout << "Interaction: LambdaInt: \n"
                << " input energy: " << p.GetEnergy() / 1_GeV << endl
                << " beam can interact:" << kInteraction << endl
                << " beam pid:" << p.GetPID() << endl;

      if (kInteraction && Elab >= 8.5_GeV && ECoM >= 10_GeV) {

        // get target from environment
        /*
          the target should be defined by the Environment,
          ideally as full particle object so that the four momenta
          and the boosts can be defined..
        */
        const auto currentNode =
            fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
        const auto mediumComposition =
            currentNode->GetModelProperties().GetNuclearComposition();
        // determine average interaction length
        // weighted sum
        int i = -1;
        double avgTargetMassNumber = 0.;
        si::CrossSectionType weightedProdCrossSection = 0_mbarn;
        // get weights of components from environment/medium
        const auto w = mediumComposition.GetFractions();
        // loop over components in medium
        for (auto const targetId : mediumComposition.GetComponents()) {
          i++;
          cout << "Interaction: get interaction length for target: " << targetId << endl;

          auto const [productionCrossSection, numberOfNucleons] =
              GetCrossSection(corsikaBeamId, targetId, ECoM);

          std::cout << "Interaction: "
                    << " IntLength: sibyll return (mb): "
                    << productionCrossSection / 1_mbarn << std::endl;
          weightedProdCrossSection += w[i] * productionCrossSection;
          avgTargetMassNumber += w[i] * numberOfNucleons;
        }
        cout << "Interaction: "
             << "IntLength: weighted CrossSection (mb): "
             << weightedProdCrossSection / 1_mbarn << endl;

        // calculate interaction length in medium
        //#warning check interaction length units
        GrammageType const int_length =
            avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
        std::cout << "Interaction: "
                  << "interaction length (g/cm2): "
                  << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;

        return int_length;
      }

      return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
    }

    /**
       In this function SIBYLL is called to produce one event. The
       event is copied (and boosted) into the shower lab frame.
     */

    template <typename Particle, typename Stack>
    corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {

      using namespace corsika::units;
      using namespace corsika::utl;
      using namespace corsika::units::si;
      using namespace corsika::geometry;
      using std::cout;
      using std::endl;

      const auto corsikaBeamId = p.GetPID();
      cout << "ProcessSibyll: "
           << "DoInteraction: " << corsikaBeamId << " interaction? "
           << process::sibyll::CanInteract(corsikaBeamId) << endl;
      if (process::sibyll::CanInteract(corsikaBeamId)) {
        const CoordinateSystem& rootCS =
            RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();

        // position and time of interaction, not used in Sibyll
        Point pOrig = p.GetPosition();
        TimeType tOrig = p.GetTime();

        // define target
        // for Sibyll is always a single nucleon
        auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
                                             corsika::particles::Neutron::GetMass());
        // FOR NOW: target is always at rest
        const auto eTargetLab = 0_GeV + nucleon_mass;
        const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
        const FourVector PtargLab(eTargetLab, pTargetLab);

        // define projectile
        HEPEnergyType const eProjectileLab = p.GetEnergy();
        auto const pProjectileLab = p.GetMomentum();

        cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
             << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
             << endl;
        cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
             << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
             << endl;

        const FourVector PprojLab(eProjectileLab, pProjectileLab);

        // define target kinematics in lab frame
        // define boost to and from CoM frame
        // CoM frame definition in Sibyll projectile: +z
        COMBoost const boost(PprojLab, nucleon_mass);

        // just for show:
        // boost projecticle
        auto const PprojCoM = boost.toCoM(PprojLab);

        // boost target
        auto const PtargCoM = boost.toCoM(PtargLab);

        cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
             << endl
             << "Interaction: pbeam CoM: "
             << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
        cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
             << endl
             << "Interaction: ptarget CoM: "
             << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;

        cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
             << endl;
        cout << "Interaction: time: " << tOrig << endl;

        HEPEnergyType Etot = eProjectileLab + eTargetLab;
        MomentumVector Ptot = p.GetMomentum();
        // invariant mass, i.e. cm. energy
        HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());

        // sample target mass number
        const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
        const auto& mediumComposition =
            currentNode->GetModelProperties().GetNuclearComposition();
        // get cross sections for target materials
        /*
          Here we read the cross section from the interaction model again,
          should be passed from GetInteractionLength if possible
         */
        //#warning reading interaction cross section again, should not be necessary
        auto const& compVec = mediumComposition.GetComponents();
        std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());

        for (size_t i = 0; i < compVec.size(); ++i) {
          auto const targetId = compVec[i];
          const auto [sigProd, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
          cross_section_of_components[i] = sigProd;
          int ideleteme = nNuc;  // to avoid not used warning in array binding
          ideleteme = ideleteme; // to avoid not used warning in array binding
        }

        const auto targetCode = currentNode->GetModelProperties().SampleTarget(
            cross_section_of_components, fRNG);
        cout << "Interaction: target selected: " << targetCode << endl;
        /*
          FOR NOW: allow nuclei with A<18 or protons only.
          when medium composition becomes more complex, approximations will have to be
          allowed air in atmosphere also contains some Argon.
        */
        int targetSibCode = -1;
        if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
        if (targetCode == corsika::particles::Proton::GetCode()) targetSibCode = 1;
        cout << "Interaction: sibyll code: " << targetSibCode << endl;
        if (targetSibCode > 18 || targetSibCode < 1)
          throw std::runtime_error(
              "Sibyll target outside range. Only nuclei with A<18 or protons are "
              "allowed.");

        // beam id for sibyll
        const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);

        std::cout << "Interaction: "
                  << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
                  << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
        if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
          std::cout << "Interaction: "
                    << " DoInteraction: should have dropped particle.. "
                    << "THIS IS AN ERROR" << std::endl;
          // p.Delete(); delete later... different process
        } else {
          fCount++;
          // Sibyll does not know about units..
          const double sqs = Ecm / 1_GeV;
          // running sibyll, filling stack
          sibyll_(kBeam, targetSibCode, sqs);
          // running decays
          // setTrackedParticlesStable();
          decsib_();
          // print final state
          int print_unit = 6;
          sib_list_(print_unit);
          fNucCount += get_nwounded() - 1;

          // delete current particle
          p.Delete();

          // add particles from sibyll to stack
          // link to sibyll stack
          SibStack ss;

          MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
          HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
          for (auto& psib : ss) {

            // skip particles that have decayed in Sibyll
            if (psib.HasDecayed()) continue;

            // transform energy to lab. frame
            auto const pCoM = psib.GetMomentum();
            HEPEnergyType const eCoM = psib.GetEnergy();
            auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));

            // add to corsika stack
            auto pnew = s.NewParticle();
            pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
            pnew.SetEnergy(Plab.GetTimeLikeComponent());
            pnew.SetMomentum(Plab.GetSpaceLikeComponents());
            pnew.SetPosition(pOrig);
            pnew.SetTime(tOrig);

            Plab_final += pnew.GetMomentum();
            Elab_final += pnew.GetEnergy();
            Ecm_final += psib.GetEnergy();
          }
          std::cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV
                    << std::endl
                    << "Elab_final=" << Elab_final / 1_GeV
                    << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
                    << std::endl;
        }
      }
      return process::EProcessReturn::eOk;
    }

  private:
    corsika::environment::Environment const& fEnvironment;
    corsika::random::RNG& fRNG =
        corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
  };

} // namespace corsika::process::sibyll

#endif