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/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/utl/COMBoost.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.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;

      corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
      rmng.RegisterRandomStream("s_rndm");
     
      // 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::hep::EnergyType& 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);