HeavyIons.h

// HeavyIons.h is a part of the PYTHIA event generator.
// Copyright (C) 2020 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This file contains the definition of the HeavyIons class which
// Provides Pythia with infrastructure to combine several nucleon
// collisions into a single heavy ion collision event. This file also
// includes the definition of the Angantyr class which implements the
// default model for heavy ion collisions in Pythia.

#ifndef Pythia8_HeavyIons_H
#define Pythia8_HeavyIons_H

#include "Pythia8/HIUserHooks.h"
#include "Pythia8/PhysicsBase.h"

namespace Pythia8 {

class Pythia;

//==========================================================================


class HeavyIons : public PhysicsBase {

public:

  HeavyIons(Pythia & mainPythiaIn)
    : mainPythiaPtr(&mainPythiaIn), HIHooksPtr(0),
      pythia(vector<Pythia*>(1, &mainPythiaIn)) {
  }

  virtual ~HeavyIons() {}

  virtual bool init() = 0;

  virtual bool next() = 0;

  static void addSpecialSettings(Settings & settings);

  static bool isHeavyIon(Settings & settings);

  bool setHIUserHooksPtr(HIUserHooksPtr userHooksPtrIn) {
    HIHooksPtr = userHooksPtrIn; return true;
  }

protected:

  void sumUpMessages(Info & in, string tag, const Info * other);

  void updateInfo();

  void clearProcessLevel(Pythia & pyt);

  static void setupSpecials(Settings & settings, string match);

  static void setupSpecials(Pythia & p, string match);

  Pythia * mainPythiaPtr;

  SigmaTotal sigtot;

  HIUserHooksPtr HIHooksPtr;

  vector<Pythia *> pythia;

  vector<string> pythiaNames;

  vector<Info*> info;

  struct InfoGrabber : public UserHooks {

    Info * getInfo() {
      return infoPtr;
    }

  };

public:

  HIInfo hiInfo;

  // Print out statistics.
  virtual void stat();

};

//==========================================================================


class Angantyr: public HeavyIons {

public:

  enum PythiaObject {
    HADRON = 0,   
    MBIAS = 1,    
    SASD = 2,    
    SIGPP = 3,    
    SIGPN = 4,    
    SIGNP = 5,    
    SIGNN = 6,    
    ALL = 7       
  };

public:

  Angantyr(Pythia & mainPythiaIn);

  virtual ~Angantyr();

  virtual bool init() override;

  virtual bool next() override;

  bool setUserHooksPtr(PythiaObject sel, UserHooksPtr userHooksPtrIn);

  vector<Nucleon>::iterator projBegin() {
   return projectile.begin();
  }
  vector<Nucleon>::iterator targBegin() {
   return target.begin();
  }
  vector<Nucleon>::iterator projEnd() {
   return projectile.end();
  }
  vector<Nucleon>::iterator targEnd() {
   return target.end();
  }

protected:

  virtual void onInitInfoPtr() override {
    registerSubObject(sigtot); }

  bool init(PythiaObject sel, string name, int n = 0);

  EventInfo mkEventInfo(Pythia &, Info &, const SubCollision * coll = 0);

  EventInfo getSignal(const SubCollision & coll);
  EventInfo getND() { return getMBIAS(0, 101); }
  EventInfo getND(const SubCollision & coll) { return getMBIAS(&coll, 101); }
  EventInfo getEl(const SubCollision & coll) { return getMBIAS(&coll, 102); }
  EventInfo getSDP(const SubCollision & coll) { return getMBIAS(&coll, 103); }
  EventInfo getSDT(const SubCollision & coll) { return getMBIAS(&coll, 104); }
  EventInfo getDD(const SubCollision & coll) { return getMBIAS(&coll, 105); }
  EventInfo getCD(const SubCollision & coll) { return getMBIAS(&coll, 106); }
  EventInfo getSDabsP(const SubCollision & coll)
    { return getSASD(&coll, 103); }
  EventInfo getSDabsT(const SubCollision & coll)
    { return getSASD(&coll, 104); }
  EventInfo getMBIAS(const SubCollision * coll, int procid);
  EventInfo getSASD(const SubCollision * coll, int procid);
  bool genAbs(const multiset<SubCollision> & coll,
              list<EventInfo> & subevents);
  void addSASD(const multiset<SubCollision> & coll);
  bool addDD(const multiset<SubCollision> & coll, list<EventInfo> & subevents);
  bool addSD(const multiset<SubCollision> & coll, list<EventInfo> & subevents);
  void addSDsecond(const multiset<SubCollision> & coll);
  bool addCD(const multiset<SubCollision> & coll, list<EventInfo> & subevents);
  void addCDsecond(const multiset<SubCollision> & coll);
  bool addEL(const multiset<SubCollision> & coll, list<EventInfo> & subevents);
  void addELsecond(const multiset<SubCollision> & coll);
  bool buildEvent(list<EventInfo> & subevents,
                  const vector<Nucleon> & proj,
                  const vector<Nucleon> & targ);

  bool setupFullCollision(EventInfo & ei, const SubCollision & coll,
                          Nucleon::Status ptype, Nucleon::Status ttype);
  bool isRemnant(const EventInfo & ei, int i, int past = 1 ) const {
    int statNow = ei.event[i].status()*past;
    if ( statNow == 63 ) return true;
    if ( statNow > 70 && statNow < 80 )
      return isRemnant(ei, ei.event[i].mother1(), -1);
    return false;
  }
  bool fixIsoSpin(EventInfo & ei);
  EventInfo & shiftEvent(EventInfo & ei);
  static int getBeam(Event & ev, int i);

  bool nextSASD(int proc);

  bool addNucleonExcitation(EventInfo & orig, EventInfo & add,
                            bool colConnect = false);

  // Find the recoilers in the current event to conserve energy and
  // momentum in addNucleonExcitation.
  vector<int> findRecoilers(const Event & e, bool tside, int beam, int end,
                            const Vec4 & pdiff, const Vec4 & pbeam);

  void addSubEvent(Event & evnt, Event & sub);
  static void addJunctions(Event & evnt, Event & sub, int coloff);

  bool addNucleusRemnants(const vector<Nucleon> & proj,
                          const vector<Nucleon> & targ);

public:

  static bool
  getTransforms(Vec4 p1, Vec4 p2, const Vec4 & p1p,
                pair<RotBstMatrix,RotBstMatrix> & R12, int, int);
  static double mT2(const Vec4 & p) { return p.pPos()*p.pNeg(); }
  static double mT(const Vec4 & p) { return sqrt(max(mT2(p), 0.0)); }

private:

  // Private UserHooks class to select a specific process.
  struct ProcessSelectorHook: public UserHooks {

    ProcessSelectorHook(): proc(0), b(-1.0) {}

    // Yes we can veto event after process-level selection.
    virtual bool canVetoProcessLevel() {
      return true;
    }

    // Veto any unwanted process.
    virtual bool doVetoProcessLevel(Event&) {
      return proc > 0 && infoPtr->code() != proc;
    }

    // Can set the overall impact parameter for the MPI treatment.
    virtual bool canSetImpactParameter() const {
      return b >= 0.0;
    }

    // Set the overall impact parameter for the MPI treatment.
    virtual double doSetImpactParameter() {
      return b;
    }

    // The wanted process;
    int proc;

    // The selected b-value.
    double b;

  };

  // Holder class to temporarily select a specific process
  struct HoldProcess {

    // Set the given process for the given hook object.
    HoldProcess(shared_ptr<ProcessSelectorHook> hook, int proc,
      double b = -1.0) : saveHook(hook), saveProc(hook->proc), saveB(hook->b) {
      hook->proc = proc;
      hook->b = b;
    }

    // Reset the process of the hook object given in the constructor.
    ~HoldProcess() {
      if ( saveHook ) {
        saveHook->proc = saveProc;
        saveHook->b = saveB;
      }
    }

    // The hook object.
    shared_ptr<ProcessSelectorHook> saveHook;

    // The previous process of the hook object.
    int saveProc;

    // The previous b-value of the hook object.
    double saveB;

  };

  // The process selector for standard minimum bias processes.
  shared_ptr<ProcessSelectorHook> selectMB;

  // The process selector for the SASD object.
  shared_ptr<ProcessSelectorHook> selectSASD;

private:

  static const int MAXTRY = 999;
  static const int MAXEVSAVE = 999;

  vector<Nucleon> projectile;
  vector<Nucleon> target;

  multiset<SubCollision> subColls;

  bool hasSignal;

  ImpactParameterGenerator * bGenPtr;

  NucleusModel * projPtr;
  NucleusModel * targPtr;

  SubCollisionModel * collPtr;

  int recoilerMode;

  int bMode;

public:

  class Redirect {

  public:

    Redirect(ostream & os, ostream & ros)
      : osSave(&os), bufferSave(os.rdbuf()) {
      osSave->rdbuf(ros.rdbuf());
    }

    ~Redirect() {
      osSave->rdbuf(bufferSave);
    }

    ostream * osSave;

    std::streambuf * bufferSave;

  };

};

//==========================================================================

} // end namespace Pythia8

#endif // Pythia8_HeavyIons_H