Particle.h

// -*- C++ -*-
//
// Particle.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_Particle_H
#define ThePEG_Particle_H
// This is the decalaration of the Particle class.
 
#include "EventConfig.h"
#include "ThePEG/Vectors/Lorentz5Vector.h"
#include "ThePEG/Vectors/LorentzRotation.h"
#include "ThePEG/Utilities/ClassDescription.h"
#include "ThePEG/EventRecord/MultiColour.h"
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ThePEG/PDT/ParticleData.h"
 
namespace ThePEG {
 
class Particle: public EventRecordBase {
 
public:
 
  friend class Event;
  friend class Collision;
  friend class Step;
  friend class SubProcess;
  friend class ParticleData;
 
  struct ParticleRep;
 
public:
 
  Particle(tcEventPDPtr newData) : theData(newData), theRep(0), theStatus(0) {}
 
  Particle(const Particle &);
 
  virtual ~Particle();
 
  bool decayed() const {
    return hasRep() && !rep().theChildren.empty();
  }
 
  const ParticleVector & children() const {
    static const ParticleVector null;
    return hasRep() ? rep().theChildren : null;
  }
 
  void addChild(tPPtr c) {
      rep().theChildren.push_back(c);
      (c->rep()).theParents.push_back(this);
  }
 
  void abandonChild(tPPtr child) {
    removeChild(child);
    child->removeParent(this);
  }
 
  const tParticleVector & parents() const {
    static const tParticleVector null;
    return hasRep() ? rep().theParents : null;
  }
 
  tParticleSet siblings() const;
 
  void undecay() {
    if ( hasRep() ) {
      rep().theChildren.clear();
      rep().theNext = tPPtr();
    }
  }
 
  void decayMode(tDMPtr dm) { rep().theDecayMode = dm; }
 
  tDMPtr decayMode() const {
    return hasRep() ? rep().theDecayMode : tDMPtr();
  }
 
  tPPtr next() const {
    return hasRep() ? rep().theNext : PPtr();
  }
 
  tPPtr previous() const {
    return hasRep() ? rep().thePrevious : tPPtr();
  }
 
  tcPPtr original() const {
    return previous() ? tcPPtr(previous()->original()) : tcPPtr(this);
  }
 
  tPPtr original() {
    return previous() ? previous()->original() : tPPtr(this);
  }
 
 
  tcPPtr final() const {
    return next() ? tcPPtr(next()->final()) : tcPPtr(this);
  }
 
  tPPtr final() {
    return next() ? next()->final() : tPPtr(this);
  }
 
 
  tStepPtr birthStep() const { 
    return hasRep() ? rep().theBirthStep : tStepPtr();
  }
 
  int number() const { 
    return hasRep() ? rep().theNumber : 0; 
  }
 
  int status() const { return theStatus; }
 
  void status(int n) { theStatus = n; }
 
  const ParticleDataClass & data() const { return *theData; }
 
  tcEventPDPtr dataPtr() const { return theData; }
 
  const string & PDGName() const { return data().PDGName(); }
 
  long id() const { return data().id(); }
 
  const Lorentz5Momentum & momentum() const { return theMomentum; }
 
  void set3Momentum(const Momentum3 & p) {
    theMomentum.setVect(p);
    theMomentum.rescaleEnergy();
  }
 
  void setMomentum(const LorentzMomentum & p) {
    theMomentum = p;
  }
 
  void set5Momentum(const Lorentz5Momentum & p) {
    theMomentum = p;
  }
 
  Energy mass() const { return momentum().mass(); }
 
  Energy nominalMass() const { return data().mass(); }
 
  Energy2 scale() const { 
    return hasRep() ? rep().theScale : -1.0*GeV2;
  }
 
  void scale(Energy2 q2) { rep().theScale = q2; }
 
  Energy2 vetoScale() const { 
    return hasRep() ? rep().theVetoScale : -1.0*GeV2;
  }
 
  void vetoScale(Energy2 q2) { rep().theVetoScale = q2; }
 
  Energy2 mt2() const { return sqr(momentum().t()) - sqr(momentum().z()); }
 
  Energy mt() const { return sqrt(mt2()); }
 
  Energy2 perpmass2() const { return momentum().perp2() + momentum().mass2(); }
 
  Energy perpmass() const { return sqrt(perpmass2()); }
 
  double rapidity() const {
    return ( Pplus() > ZERO && Pminus() > ZERO )?
      0.5*log(Pplus()/Pminus()) : Constants::MaxFloat;
  }
 
  double eta() const {
    Energy rho = momentum().rho();
    return rho > abs(momentum().z())?
      0.5*log((rho+momentum().z())/(rho-momentum().z())) : Constants::MaxFloat;
  }
 
  Energy Pplus() const { return momentum().plus(); }
  Energy Pminus() const { return momentum().minus(); }
 
  const LorentzPoint & vertex() const {
    static const LorentzPoint null;
    return hasRep() ? rep().theVertex : null;
  }
 
  LorentzPoint labVertex() const;
 
  LorentzPoint decayVertex() const { 
    return vertex() + lifeLength();
  }
 
  LorentzPoint labDecayVertex() const {
    return labVertex() + lifeLength();
  }
 
  const Lorentz5Distance & lifeLength() const {
    static const Lorentz5Distance null;
    return hasRep() ? rep().theLifeLength : null;
  }
 
  void setVertex(const LorentzPoint & p) {
    rep().theVertex = p;
  }
 
  void setLabVertex(const LorentzPoint &);
 
  void setLifeLength(const Distance & d) {
    rep().theLifeLength.setVect(d);
    rep().theLifeLength.rescaleEnergy();
  }
 
  void setLifeLength(const LorentzDistance & d) {
    rep().theLifeLength = d;
  }
 
  void setLifeLength(const Lorentz5Distance & d) {
    rep().theLifeLength = d;
  }
 
  Time lifeTime() const { return lifeLength().m(); }
 
 
  void transform(const LorentzRotation & r);
 
  void boost(double bx, double by, double bz) {
    transform(LorentzRotation(Boost(bx, by, bz)));
  }
 
  void boost(const Boost & b) { transform(LorentzRotation(b)); }
 
  void rotateX(double a);
 
  void rotateY(double a);
 
  void rotateZ(double a);
 
  void rotate(double a, const Axis & axis);
 
  void mirror() { theMomentum.setZ(-theMomentum.z()); }
 
  void deepTransform(const LorentzRotation & r);
 
  void deepBoost(double bx, double by, double bz) {
    deepTransform(LorentzRotation(Boost(bx, by, bz)));
  }
 
  void deepBoost(const Boost & b) { deepTransform(LorentzRotation(b)); }
 
  void deepRotateX(double a);
 
  void deepRotateY(double a);
 
  void deepRotateZ(double a);
 
  void deepRotate(double a, const Axis & axis);
 
 
  double massError() const { return theMomentum.massError(); }
 
  double energyError() const { return theMomentum.energyError(); }
 
  double rhoError() const { return theMomentum.rhoError(); }
 
  void rescaleEnergy() { theMomentum.rescaleEnergy(); }
 
  void rescaleRho() { theMomentum.rescaleRho(); }
 
  void rescaleMass() { theMomentum.rescaleMass(); }
 
  bool hasColourInfo() const {
    return hasRep() && rep().theColourInfo;
  }
 
  tColinePtr antiColourLine() const {
    return hasColourInfo() ? colourInfo()->antiColourLine() : tColinePtr();
  }
 
  tColinePtr colourLine(bool anti = false) const {
    if ( anti ) return antiColourLine();
    return hasColourInfo() ? colourInfo()->colourLine() : tColinePtr();
  }
 
  bool hasColourLine(tcColinePtr line, bool anti = false) const {
    return hasColourInfo() ? colourInfo()->hasColourLine(line, anti) : false;
  }
 
  bool hasAntiColourLine(tcColinePtr line) const {
    return hasColourLine(line, true);
  }
 
  bool coloured() const { return data().coloured(); }
 
  bool hasColour(bool anti = false) const { return data().hasColour(anti); }
 
  bool hasAntiColour() const { return data().hasAntiColour(); }
 
  tcCBPtr colourInfo() const {
    return hasRep() ? rep().theColourInfo : CBPtr();
  }
 
  tCBPtr colourInfo() {
    if ( !rep().theColourInfo ) {
      switch(theData->iColour()) {
      case PDT::Colour6: 
      case PDT::Colour6bar:
        rep().theColourInfo = new_ptr(MultiColour());
        break;
      default:
        rep().theColourInfo = new_ptr(ColourBase());
      }
    }
    return rep().theColourInfo;
  }
 
  void colourInfo(tCBPtr c) {
    rep().theColourInfo = c;
  }
 
  template <typename Iterator>
  typename std::iterator_traits<Iterator>::value_type
  colourNeighbour(Iterator first, Iterator last, bool anti = false) const;
 
  template <typename Iterator>
  typename std::iterator_traits<Iterator>::value_type
  antiColourNeighbour(Iterator first, Iterator last) const {
    return colourNeighbour(first, last, true);
  }
 
  void colourNeighbour(tPPtr, bool anti = false);
 
  void antiColourNeighbour(tPPtr p) { colourNeighbour(p, true); }
 
  void antiColourConnect(tPPtr neighbour) {
    colourConnect(neighbour, true);
  }
 
  void colourConnect(tPPtr neighbour, bool anti = false) {
    colourNeighbour(neighbour, anti);
  }
 
  tPPtr incomingColour(bool anti = false) const;
 
  tPPtr incomingAntiColour() const { return incomingColour(true); }
 
  void incomingColour(tPPtr p, bool anti = false) { p->outgoingColour(this, anti); }
 
  void incomingAntiColour(tPPtr p) { p->outgoingColour(this, true); }
 
  tPPtr outgoingColour(bool anti = false) const;
  tPPtr outgoingAntiColour() const { return outgoingColour(true); }
 
  void outgoingColour(tPPtr, bool anti = false);
 
  void outgoingAntiColour(tPPtr p) { outgoingColour(p, true); }
 
  void colourFlow(tPPtr child, bool anti = false) {
    outgoingColour(child, anti);
  }
 
  void antiColourFlow(tPPtr child) { colourFlow(child, true); }
 
  void resetColour() {
    if ( hasColourInfo() ) rep().theColourInfo = CBPtr();
  }
 
 
  tcSpinPtr spinInfo() const { 
    return hasRep() ? rep().theSpinInfo : SpinPtr(); 
  }
 
  tSpinPtr spinInfo() {
    return hasRep() ? rep().theSpinInfo : SpinPtr(); 
  }
 
  void spinInfo(tSpinPtr s) { rep().theSpinInfo = s; }
 
  const EIVector & getInfo() const {
    static const EIVector null;
    return hasRep() ? rep().theExtraInfo : null;
  }
 
  EIVector & getInfo() { return rep().theExtraInfo; }
 
public:
 
  bool hasRep() const { return theRep; }
 
  void initFull();
 
 
public:
 
  void persistentOutput(PersistentOStream &) const;
 
  void persistentInput(PersistentIStream &, int);
 
 
  ostream & print(ostream & os, tcStepPtr step = tcStepPtr()) const;
 
  template <typename Iterator>
  static void PrintParticles(ostream & os, Iterator first, Iterator last,
                             tcStepPtr step = tcStepPtr());
 
  template <typename Cont>
  static inline void PrintParticles(ostream & os, const Cont & c,
                             tcStepPtr step = tcStepPtr()) {
    PrintParticles(os, c.begin(), c.end(), step);
  }
   
  static void Init();
 
  static string outputFormat;
 
private:
 
  virtual PPtr clone() const;
 
  virtual void rebind(const EventTranslationMap &);
 
  void number(int n) { rep().theNumber = n; }
 
  void removeChild(tPPtr c) {
    if ( hasRep() )
      rep().theChildren.erase(remove(rep().theChildren.begin(),
                                     rep().theChildren.end(), c),
                              rep().theChildren.end());
  }
 
  void removeParent(tPPtr p) {
    if ( hasRep() )
      rep().theParents.erase(remove(rep().theParents.begin(),
                                    rep().theParents.end(), p),
                             rep().theParents.end());
  }
 
  void mass(Energy m) { theMomentum.setMass(m); }
 
  void lifeTime(Length t) { rep().theLifeLength.setTau(t); }
 
  ParticleRep & rep() {
    if ( !hasRep() ) initFull();
    return *theRep;
  }
 
  const ParticleRep & rep() const {
    static const ParticleRep null;
    return hasRep() ? *theRep : null;
  }
 
  cEventPDPtr theData;
 
  Lorentz5Momentum theMomentum;
 
  ParticleRep * theRep;
 
  int theStatus;
 
public:
 
  struct ParticleRep {
 
    ParticleRep() : theScale(-1.0*GeV2), theVetoScale(-1.0*GeV2), theNumber(0) {}
 
    ParticleRep(const ParticleRep &);
 
    tParticleVector theParents;
 
    ParticleVector theChildren;
 
    tPPtr thePrevious;
 
    PPtr theNext;
 
    tDMPtr theDecayMode;
 
    tStepPtr theBirthStep;
 
    LorentzPoint theVertex;
 
    Lorentz5Distance theLifeLength;
 
    Energy2 theScale;
 
    Energy2 theVetoScale;
 
    int theNumber;
 
    CBPtr theColourInfo;
 
    SpinPtr theSpinInfo;
 
    EIVector theExtraInfo;
 
  };
 
public:
 
  virtual void debugme() const;
 
protected:
 
  Particle() : theRep(0), theStatus(0) {}
 
  friend struct ClassTraits<Particle>;
 
private:
 
  Particle & operator=(const Particle &) = delete;
 
  static ClassDescription<Particle> initParticle;
 
};
 
ostream & operator<<(ostream &, const Particle &);
 
 
template <>
struct BaseClassTrait<Particle,1>: public ClassTraitsType {
  typedef EventRecordBase NthBase;
};
 
template <>
struct ClassTraits<Particle>: public ClassTraitsBase<Particle> {
  static string className() { return "ThePEG::Particle"; }
  static TPtr create() { return TPtr::Create(Particle()); }
};
 
}
 
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "Particle.tcc"
#endif
 
#endif /* ThePEG_Particle_H */