Decay.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 _include_corsika_process_sibyll_decay_h_
#define _include_corsika_process_sibyll_decay_h_

#include <corsika/process/DecayProcess.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>

#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>

#include <corsika/particles/ParticleProperties.h>

#include <fenv.h>

namespace corsika::process {

  namespace sibyll {

    class Decay : public corsika::process::DecayProcess<Decay> {
      int fCount = 0;

    public:
      Decay() {}
      ~Decay() { std::cout << "Sibyll::Decay n=" << fCount << std::endl; }
      void Init() {
        setHadronsUnstable();
        setTrackedParticlesStable();
      }

      void setTrackedParticlesStable() {
        /*
           Sibyll is hadronic generator
           only hadrons decay
         */
        // set particles unstable
        setHadronsUnstable();
        // make tracked particles stable
        std::cout << "Interaction: setting tracked hadrons stable.." << std::endl;
        const std::vector<corsika::particles::Code> particleList = {
            particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
            particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};

        for (auto p : particleList) {
          // set particle stable by setting table value negative
          const int sibid = process::sibyll::ConvertToSibyllRaw(p);
          s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
        }
      }

      void setUnstable(const corsika::particles::Code pCode) {
        int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
        s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
      }

      void setStable(const corsika::particles::Code pCode) {
        int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
        s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
      }

      void setAllStable() {
        // name? also makes EM particles stable

        using std::cout;
        using std::endl;

        std::cout << "Decay: setting all particles stable.." << std::endl;

        // loop over all particles in sibyll
        // should be changed to loop over human readable list
        // i.e. corsika::particles::ListOfParticles()
        for (auto& p : corsika2sibyll) {
          // std::cout << (int)p << std::endl;
          const int sibCode = (int)p;
          // skip unknown and antiparticles
          if (sibCode < 1) continue;
          s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
        }
      }

      void setHadronsUnstable() {

        using std::cout;
        using std::endl;
        // using namespace corsika::io;
        using namespace corsika::units::si;

        // name? also makes EM particles stable

        // loop over all particles in sibyll
        // should be changed to loop over human readable list
        // i.e. corsika::particles::ListOfParticles()
        std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
        // make ALL particles unstable, then set EM stable
        for (int sibCode : corsika2sibyll) {
          // std::cout << (int)p << std::endl;
          // skip unknown and antiparticles
          if (sibCode < 1) continue;
          // std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
          // static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
          s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
          // std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
          // std::endl;
        }
        // set Leptons and Proton and Neutron stable
        // use stack to loop over particles
        const std::vector<corsika::particles::Code> particleList = {
            corsika::particles::Code::Proton,   corsika::particles::Code::Neutron,
            corsika::particles::Code::Electron, corsika::particles::Code::Positron,
            corsika::particles::Code::NuE,      corsika::particles::Code::NuEBar,
            corsika::particles::Code::MuMinus,  corsika::particles::Code::MuPlus,
            corsika::particles::Code::NuMu,     corsika::particles::Code::NuMuBar};

        for (auto p : particleList) {
          // set particle stable by setting table value negative
          //	cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
          //     << " stable in Sibyll .." << endl;
          const int sibid = process::sibyll::ConvertToSibyllRaw(p);
          s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
        }
      }

      template <typename Particle>
      corsika::units::si::TimeType GetLifetime(Particle& p) {
        using std::cout;
        using std::endl;
        using namespace corsika::units::si;

        corsika::units::hep::EnergyType E = p.GetEnergy();
        corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());

        const double gamma = E / m;

        const corsika::units::si::TimeType t0 =
            corsika::particles::GetLifetime(p.GetPID());
        auto const lifetime = gamma * t0;

        const auto mkin =
            (E * E - p.GetMomentum().squaredNorm()); // delta_mass(p.GetMomentum(), E, m);
        cout << "Decay: code: " << p.GetPID() << endl;
        cout << "Decay: MinStep: t0: " << t0 << endl;
        cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
        cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV"
             << endl;
        cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV
             << " " << m / 1_GeV * m / 1_GeV << endl;
        // cout << "Decay: sib mass: " << s_mass1_.am2[
        // process::sibyll::ConvertToSibyllRaw(p.GetPID()) ] << endl;
        auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
        cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
        cout << "Decay: MinStep: gamma: " << gamma << endl;
        cout << "Decay: MinStep: tau: " << lifetime << endl;

        return lifetime;
      }

      template <typename Particle, typename Stack>
      void DoDecay(Particle& p, Stack& s) {
        using corsika::geometry::Point;
        using namespace corsika::units::si;

	// TODO: this should be done in a central, common place. Not here..
#ifndef CORSIKA_OSX
        feenableexcept(FE_INVALID);
#endif

        fCount++;
        SibStack ss;
        ss.Clear();
        // copy particle to sibyll stack
        auto pin = ss.NewParticle();
        const corsika::particles::Code pCode = p.GetPID();
        pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
        pin.SetEnergy(p.GetEnergy());
        pin.SetMomentum(p.GetMomentum());
        // setting particle mass with Corsika values, may be inconsistent with sibyll
        // internal values
	// TODO: #warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
        pin.SetMass(corsika::particles::GetMass(pCode));
        // remember position
        Point decayPoint = p.GetPosition();
        TimeType t0 = p.GetTime();
        // remove original particle from corsika stack
        p.Delete();
        // set all particles/hadrons unstable
        // setHadronsUnstable();
        setUnstable(pCode);
        // call sibyll decay
        std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
        decsib_();
        // reset to stable
        setStable(pCode);
        // print output
        int print_unit = 6;
        sib_list_(print_unit);

        // copy particles from sibyll stack to corsika
        for (auto& psib : ss) {
          // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
          if (psib.HasDecayed()) continue;
          // add to corsika stack
          auto pnew = s.NewParticle();
          pnew.SetEnergy(psib.GetEnergy());
          pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
          pnew.SetMomentum(psib.GetMomentum());
          pnew.SetPosition(decayPoint);
          pnew.SetTime(t0);
        }
        // empty sibyll stack
        ss.Clear();

	// TODO: this should be done in a central, common place. Not here..
#ifndef CORSIKA_OSX
        fedisableexcept(FE_INVALID);
#endif
      }
    };
  } // namespace sibyll
} // namespace corsika::process

#endif