HIUserHooks.h

// HIUserHooks.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 HIUserHooks class and a
// set of other classes used inside Pythia8 to model collisions
// involving heavy ions.
// Nucleon: represents a proton or a neutron inside a necleus.
// SubCollision: a collision between a projectile and a target Nucleon.
// NucleusModel: models the Nucleon distribution in a nucleus.
// WoodsSaxonModel: NucleusModel implementing a simple Woods-Saxon.
// GLISSANDOModel: NucleusModel implementing the GLISSANDO prescription.
// ImpactParameterGenerator: distributes nuclei in impact parameter space.
// SubCollisionModel: Models the collision probabilities of nucleons.
// NaiveSubCollisionModel: A very simple SubCollisionModel.
// DoubleStrikman: A more advanced SubCollisionModel.
// EventInfo: stores full nucleon-nucleon events with corresponding Info.
// HIInfo: info about a Heavy Ion run and its produced events.
// HIUserHooks: User hooks for HeavyIons models.

#ifndef Pythia8_HIUserHooks_H
#define Pythia8_HIUserHooks_H

#include "Pythia8/Pythia.h"
#include <limits>

namespace Pythia8 {

class Pythia;

class EventInfo;

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

// The HIUnits namespace defines the unitsystem used by the heavy ion
// machinery in Pythia8. In particular all lengths are in femtometer
// and cross sections are therefore in squared femtometer.
namespace HIUnits {

// Lengths
const double femtometer = 1.0;
const double millimeter = 1.0e12;

// Cross sections
const double femtometer2 = 1.0;
const double millibarn = 0.1;
const double nanobarn = 1.0e-7;

}

using namespace HIUnits;

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


class Nucleon {

  friend class SubCollisionModel;

public:

  enum Status {
    UNWOUNDED = 0,  
    ELASTIC = 1,    
    DIFF = 2,       
    ABS = 3         
  };

  typedef vector<double> State;

  Nucleon(int idIn = 0, int indexIn = 0, const Vec4 & pos = Vec4())
    : idSave(idIn), indexSave(indexIn), nPosSave(pos), bPosSave(pos),
      statusSave(UNWOUNDED),eventp(0), isDone(0) {}


  int id() const { return idSave; }

  int index() const { return indexSave; }

  const Vec4 & nPos() const { return nPosSave; }

  const Vec4 & bPos() const { return bPosSave; }

  const Vec4 & bShift(const Vec4 & bvec) { return bPosSave += bvec; }

  Status status() const { return statusSave; }

  bool done() const { return isDone; }

  EventInfo * event() const { return eventp; }

  const State & state() const { return stateSave; }

  const State & altState(int i = 0) {
    static State nullstate;
    return i < int(altStatesSave.size())? altStatesSave[i]: nullstate;
  }


  void status(Status s) { statusSave = s; }

  void state(State s) { stateSave = s; }

  void addAltState(State s) { altStatesSave.push_back(s); }

  void select(EventInfo & evp, Status s) {
    eventp = &evp;
    isDone = true;
    status(s);
  }

  void select() { isDone = true; }

  void debug();

private:

  int idSave;

  int indexSave;

  Vec4 nPosSave;

  Vec4 bPosSave;

  Status statusSave;

  State stateSave;

  vector<State> altStatesSave;

  EventInfo * eventp;

  bool isDone;

  void reset() {
    statusSave = UNWOUNDED;
    altStatesSave.clear();
    bPosSave = nPosSave;
    isDone = false;
    eventp = 0;
  }

};

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


class SubCollision {

public:

  enum Type {
    NONE,       
    ELASTIC,    
    SDEP,       
    SDET,       
    DDE,        
    CDE,        
    ABS         
  };

  SubCollision(Nucleon & projIn, Nucleon & targIn,
               double bIn, double bpIn, Type typeIn)
    : proj(&projIn), targ(&targIn), b(bIn), bp(bpIn), type(typeIn) {}

  SubCollision()
    : proj(0), targ(0), b(0.0), bp(0.0), type(NONE) {}

  // Used to order sub-collisions in a set.
  bool operator< (const SubCollision & s) const { return b < s.b; }

  // Return 0 if neither proj or target are neutrons, 1 if target is
  // neutron, 2 if projectile is neutron, and 3 if both are neutrons.
  int nucleons() const {
    return ( abs(targ->id()) == 2112? 1: 0 ) +
           ( abs(proj->id()) == 2112? 2: 0 );
  }

  // The projectile nucleon.
  Nucleon * proj;

  Nucleon * targ;

  double b;

  double bp;

  mutable Type type;

};

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


class NucleusModel {

public:

  NucleusModel()
    : idSave(0), ISave(0), ASave(0), ZSave(0), LSave(0), RSave(0.0),
      settingsPtr(0), particleDataPtr(0), rndPtr(0) {}

  virtual ~NucleusModel() {}

  void initPtr(int idIn, Settings & settingsIn,
               ParticleData & particleDataIn, Rndm & rndIn);
  virtual bool init();

  virtual Particle produceIon(bool istarg);

  virtual vector<Nucleon> generate() const = 0;

  int id() const { return idSave; }
  int I() const { return ISave; }
  int A() const { return ASave; }
  int Z() const { return ZSave; }
  int L() const { return LSave; }
  double R() const { return RSave; }

protected:

  int idSave;

  int ISave, ASave, ZSave, LSave;

  double RSave;

  Settings * settingsPtr;
  ParticleData * particleDataPtr;
  Rndm * rndPtr;

};

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


class WoodsSaxonModel: public NucleusModel {

public:

  WoodsSaxonModel(): aSave(0.0), intlo(0.0),
                    inthi0(0.0), inthi1(0.0), inthi2(0.0) {}

  double a() const { return aSave; }

protected:

  Vec4 generateNucleon() const;

  virtual ~WoodsSaxonModel() {}

  virtual bool init() {
    intlo = R()*R()*R()/3.0;
    inthi0 = a()*R()*R();
    inthi1 = 2.0*a()*a()*R();
    inthi2 = 2.0*a()*a()*a();
    return NucleusModel::init();
  }

protected:

  double aSave;

private:

  double intlo, inthi0, inthi1, inthi2;

};


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


class GLISSANDOModel: public WoodsSaxonModel {

public:

  GLISSANDOModel(): RhSave(0.0), gaussHardCore(false) {}

  virtual ~GLISSANDOModel() {}

  bool init();

  virtual vector<Nucleon> generate() const;

  double Rh() const { return RhSave; }

  double RhGauss() const { return RhSave*abs(rndPtr->gauss()); };

private:

  double RhSave;

  bool gaussHardCore;

};

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

class SubCollisionModel;
class NucleusModel;


class ImpactParameterGenerator {

public:

  ImpactParameterGenerator()
    : widthSave(0.0), collPtr(0), projPtr(0), targPtr(0),
      settingsPtr(0), rndPtr(0) {}

  virtual ~ImpactParameterGenerator() {}

  virtual bool init();
  void initPtr(SubCollisionModel & collIn,
               NucleusModel & projIn,
               NucleusModel & targIn,
               Settings & settingsIn,
               Rndm & rndIn);

  virtual Vec4 generate(double & weight) const;

  void width(double widthIn) { widthSave = widthIn; }

  double width() const { return widthSave; }

private:

  double widthSave;

protected:

  SubCollisionModel * collPtr;
  NucleusModel * projPtr;
  NucleusModel * targPtr;
  Settings * settingsPtr;
  Rndm * rndPtr;

};


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


class SubCollisionModel {

public:

  struct SigEst {
    vector<double> sig;

    vector<double> dsig2;

    vector<bool> fsig;

    double avNDb, davNDb2;

    SigEst(): sig(8, 0.0), dsig2(8, 0.0), fsig(8, false),
              avNDb(0.0), davNDb2(0.0) {}

  };

public:

  SubCollisionModel(): sigTarg(8, 0.0), sigErr(8, 0.05), NInt(100000),
    NGen(20), NPop(20), sigFuzz(0.2), fitPrint(true), avNDb(1.0*femtometer),
    projPtr(), targPtr(), sigTotPtr(), settingsPtr(), infoPtr(), rndPtr() {}

  virtual ~SubCollisionModel() {}

  virtual bool init();

  void initPtr(NucleusModel & projIn, NucleusModel & targIn,
               SigmaTotal & sigTotIn, Settings & settingsIn,
               Info & infoIn, Rndm & rndIn) {
    projPtr = &projIn;
    targPtr = &targIn;
    sigTotPtr = &sigTotIn;
    settingsPtr = &settingsIn;
    infoPtr = &infoIn;
    rndPtr = & rndIn;
  }

  virtual multiset<SubCollision> getCollisions(vector<Nucleon> & proj,
                                               vector<Nucleon> & targ,
                                               const Vec4 & bvec,
                                               double & T) = 0;


  double sigTot() const {
    return sigTarg[0];
  }

  double sigEl() const { return sigTarg[6]; }

  double sigCDE() const { return sigTarg[5]; }

  double sigSDE() const { return sigTarg[3] + sigTarg[4]; }

  double sigSDEP() const { return sigTarg[3]; }

  double sigSDET() const { return sigTarg[4]; }

  double sigDDE() const { return sigTarg[2]; }

  double sigND() const { return sigTarg[1]; }

  double bSlope() const { return sigTarg[7]; }

  virtual SigEst getSig() const {
    return SigEst();
  }

  double avNDB() const {
    return avNDb;
  }

  double Chi2(const SigEst & sigs, int npar) const;

  virtual bool evolve();

  virtual void setParm(const vector<double> &) {}

  virtual vector<double> getParm() const {
    return vector<double>();
  }
  virtual vector<double> minParm() const {
    return vector<double>();
  }
  virtual vector<double> maxParm() const {
    return vector<double>();
  }

private:

  vector<double> sigTarg, sigErr;

protected:

  int NInt, NGen, NPop;
  double sigFuzz;
  bool fitPrint;

  double avNDb;

  NucleusModel * projPtr;
  NucleusModel * targPtr;
  SigmaTotal * sigTotPtr;
  Settings * settingsPtr;
  Info * infoPtr;
  Rndm * rndPtr;

};


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


class BlackSubCollisionModel: public SubCollisionModel {

public:

  BlackSubCollisionModel() {}

  virtual ~BlackSubCollisionModel() {}

  virtual multiset<SubCollision>
  getCollisions(vector<Nucleon> & proj, vector<Nucleon> & targ,
                const Vec4 & bvec, double & T);

};

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


class NaiveSubCollisionModel: public SubCollisionModel {

public:

  NaiveSubCollisionModel() {}

  virtual ~NaiveSubCollisionModel() {}

  virtual multiset<SubCollision>
  getCollisions(vector<Nucleon> & proj, vector<Nucleon> & targ,
                const Vec4 & bvec, double & T);

};

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


class DoubleStrikman: public SubCollisionModel {

public:

  DoubleStrikman(int modein = 0)
    : r0(0.0), k0(4.0), sigd(75.0), alpha(0.5), opacityMode(modein) {}

  virtual ~DoubleStrikman() {}

  virtual multiset<SubCollision>
  getCollisions(vector<Nucleon> & proj, vector<Nucleon> & targ,
                const Vec4 & bvec, double & T);

  double gamma() const;

  double opacity(double sig) const {
    // *** THINK *** maybe sig/sigd?
    sig /= sigd;
    if ( opacityMode == 1 ) sig = 1.0/sig;
    return sig > std::numeric_limits<double>::epsilon()?
      pow(-expm1(-1.0/sig), alpha): 1.0;
  }

  double Tpt(const Nucleon::State & p,
             const Nucleon::State & t, double b) const {
    double sig = M_PI*pow2(p[0] + t[0]);
    double grey = opacity(sig);
    return sig/grey > b*b*2.0*M_PI? grey: 0.0;
  }

  SigEst getSig() const;

  virtual void setParm(const vector<double> &);

  virtual vector<double> getParm() const;
  virtual vector<double> minParm() const;
  virtual vector<double> maxParm() const;

  // Helper functions
  static void shuffle(double PND1, double PND2,
                      double & PW1, double & PW2);
  static void shuffel(double & PEL11, double P11,
                      double P12, double P21, double P22);
  static double PNW(double PWp, double PWt, double PND) {
    return ( 1.0 - PWp <= 0.0 || 1.0 - PWt <= 0.0 )?
      0.0: (1.0 - PWp)*(1.0 - PWt)/(1.0 - PND);
  }

protected:

  double r0;

  double k0;

  double sigd;

  double alpha;

  int opacityMode;

};


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


class MultiRadial: public SubCollisionModel {

public:

  MultiRadial(int NrIn = 0)
    : Nr(max(1, NrIn)) {
    dR = T0 = c = phi = vector<double>(Nr, 0.0);
  }

  virtual ~MultiRadial() {}

  virtual multiset<SubCollision>
  getCollisions(vector<Nucleon> & proj, vector<Nucleon> & targ,
                const Vec4 & bvec, double & T);

  double Tpt(const Nucleon::State & p,
             const Nucleon::State & t, double b) const {
    return b < p[0] + t[0]? p[1]*t[1]: 0.0;
  }

  SigEst getSig() const;

  virtual void setParm(const vector<double> &);

  virtual vector<double> getParm() const;
  virtual vector<double> minParm() const;
  virtual vector<double> maxParm() const;

protected:

  // Set the probabilities according to the angle parameters.
  void setProbs();

  int choose() const;

  int Nr;


  vector<double> c;

  vector<double> dR;

  vector<double> T0;

  vector<double> phi;

};


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

// Class for storing Events and Info objects.

class EventInfo {

public:

  EventInfo(): code(0), ordering(-1.0), coll(0), ok(false) {}

  Event event;

  // The corresponding Info object.
  Info info;

  int code;

  double ordering;
  bool operator<(const EventInfo & ei) const {
    return ordering < ei.ordering;
  }

  const SubCollision * coll;

  bool ok;

  map<Nucleon *, pair<int,int> > projs, targs;

};

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


class HIInfo {

public:

  friend class HeavyIons;
  friend class Angantyr;

  HIInfo()
    : idProjSave(0), idTargSave(0), bSave(0.0), NSave(0), NAccSave(0),
      sigmaTotSave(0.0), sigmaNDSave(0.0), sigErr2TotSave(0.0),
      sigErr2NDSave(0.0), weightSave(0.0), weightSumSave(0.0),
      nCollSave(10, 0), nProjSave(10, 0), nTargSave(10, 0), nFailSave(0),
      subColsPtr(NULL) {}

  double b() const {
    return bSave;
  }

  double sigmaTot() const {
    return sigmaTotSave/millibarn;
  }

  double sigmaTotErr() const {
    return sqrt(sigErr2TotSave/max(1.0, double(NSave)))/millibarn;
  }

  double sigmaND() const {
    return sigmaNDSave/millibarn;
  }

  double sigmaNDErr() const {
    return sqrt(sigErr2NDSave/max(1.0, double(NSave)));
  }

  long nAttempts() const {
    return NSave;
  }

  long nAccepted() const {
    return NAccSave;
  }

  int nCollTot() const { return nCollSave[0]; }

  int nCollND() const { return nCollSave[1]; }

  int nCollNDTot() const { return nProjSave[1] + nTargSave[1] - nCollSave[1]; }

  int nCollSDP() const { return nCollSave[2]; }

  int nCollSDT() const { return nCollSave[3]; }

  int nCollDD() const { return nCollSave[4]; }

  int nCollCD() const { return nCollSave[5]; }

  int nCollEL() const { return nCollSave[6]; }

  int nPartProj() const { return nProjSave[0]; }

  int nAbsProj() const { return nProjSave[1]; }

  int nDiffProj() const { return nProjSave[2]; }

  int nElProj() const { return nProjSave[3]; }

  int nPartTarg() const { return nTargSave[0]; }

  int nAbsTarg() const { return nTargSave[1]; }

  int nDiffTarg() const { return nTargSave[2]; }

  int nElTarg() const { return nTargSave[3]; }

  double weight() const { return weightSave; }

  double weightSum() const { return weightSumSave; }

  int nFail() const {
    return nFailSave;
  }

  void failedExcitation() {
    ++nFailSave;
  }

private:

  void addAttempt(double T, double bin, double bweight);

  void reweight(double w) {
    weightSave *= w;
  }

  // Select the primary process.
  void select(Info & in) {
    primInfo = in;
    primInfo.hiInfo = this;
  }

  // Accept an event and update statistics in info.
  void accept();

  void reject() {}

  int addSubCollision(const SubCollision & c);

  int addProjectileNucleon(const Nucleon & n);
  int addTargetNucleon(const Nucleon & n);


  int idProjSave, idTargSave;

  double bSave;

  long NSave, NAccSave;
  double sigmaTotSave, sigmaNDSave, sigErr2TotSave, sigErr2NDSave;
  double weightSave, weightSumSave;

  vector<int> nCollSave, nProjSave, nTargSave;

  // Map of primary processes and the number of events and the sum of
  // weights.
  map<int,double> sumPrimW, sumPrimW2;
  map<int,int> NPrim;
  map<int,string> NamePrim;

  // The info object of the primary process.
  Info primInfo;

  // Number of failed nucleon excitations.
  int nFailSave;


public:
  // Access to subcollision to be extracted by the user.
  multiset<SubCollision>* subCollisionsPtr() { return subColsPtr; }

  void subCollisionsPtr(multiset<SubCollision> * sPtrIn) {
    subColsPtr = sPtrIn; }

private:

  // Full information about the Glauber calculation, consisting of
  // all subcollisions.
  multiset<SubCollision>* subColsPtr;

};

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


class HIUserHooks {

public:

  HIUserHooks(): idProjSave(0), idTargSave(0) {}

  virtual ~HIUserHooks() {}

  virtual void init(int idProjIn, int idTargIn) {
    idProjSave = idProjIn;
    idTargSave = idTargIn;
  }

  virtual bool hasImpactParameterGenerator() const { return false; }
  virtual ImpactParameterGenerator * impactParameterGenerator() const {
    return 0; }

  virtual bool hasProjectileModel() const { return false; }
  virtual NucleusModel * projectileModel() const { return 0; }
  virtual bool hasTargetModel() const { return false; }
  virtual NucleusModel * targetModel() const { return 0; }

  virtual bool hasSubCollisionModel() { return false; }
  virtual SubCollisionModel * subCollisionModel() { return 0; }

  virtual bool hasEventOrdering() const { return false; }
  virtual double eventOrdering(const Event &, const Info &) { return -1; }

  virtual bool canFixIsoSpin() const { return false; }
  virtual bool fixIsoSpin(EventInfo &) { return false; }

  virtual bool canShiftEvent() const { return false; }
  virtual EventInfo & shiftEvent(EventInfo & ei) const { return ei; }

  bool canAddNucleonExcitation() const { return false; }
  bool addNucleonExcitation(EventInfo &, EventInfo &, bool) const {
    return false; }

  virtual bool canForceHadronLevel() const { return false; }
  virtual bool forceHadronLevel(Pythia &) { return false; }

  virtual bool canFindRecoilers() const { return false; }
  virtual vector<int>
  findRecoilers(const Event &, bool /* tside */, int /* beam */, int /* end */,
               const Vec4 & /* pdiff */, const Vec4 & /* pbeam */) const {
    return vector<int>();
  }

protected:

  int idProjSave, idTargSave;

};

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

} // end namespace Pythia8

#endif // Pythia8_HIUserHooks_H