VertexBase.h

// -*- C++ -*-
//
// VertexBase.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2019 Peter Richardson, 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_VertexBase_H
#define ThePEG_VertexBase_H
//
// This is the declaration of the VertexBase class.
 
#include <ThePEG/Interface/Interfaced.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/Helicity/HelicityDefinitions.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "VertexBase.fh"
#include <array>
 
 
namespace ThePEG {
namespace Helicity {
 
namespace VertexType {
  typedef unsigned T;
  const T UNDEFINED = 0;
}
 
namespace ColourStructure {
  typedef unsigned T;
  const T UNDEFINED  = 0;
  const T SINGLET    = 1;
  const T SU3TFUND   = 2;
  const T SU3F       = 3;
  const T SU3T6      = 4;
  const T SU3K6      = 5;
  const T EPS        = 6;
  const T DELTA      = 7;
  const T SU3FF      = 8;
  const T SU3TTFUNDS = 9;
  const T SU3TTFUNDD = 10;
  const T SU3TT6     = 11;
  const T SU3I12I34  = 12;
  const T SU3I14I23  = 13;
  const T SU3T21T43  = 14;
  const T SU3T23T41  = 15;
}
 
namespace CouplingType {
  typedef unsigned T;
  const T UNDEFINED = 0;
  const T QED = 1;
  const T QCD = 2;
}
class VertexBase  : public Interfaced {
friend ostream & operator<<(ostream &, const VertexBase &);
 
public:
  
  VertexBase(VertexType::T name, bool kine=false);
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
  
  static void Init();
  
public:
 
  unsigned int size() const { return _particles.size(); }
 
public:
  bool isIncoming(tPDPtr p) const {
    return _inpart.find(p) != _inpart.end();
  }
 
  bool isOutgoing(tPDPtr p) const {
    return _outpart.find(p) != _outpart.end();
  }
 
  const set<tPDPtr> & incoming() const { return _inpart; }
 
  const set<tPDPtr> & outgoing() const { return _outpart; }
 
  Complex norm() const { return _norm; }
 
  vector<long> search(unsigned int ilist,long id) const;
 
  vector<tPDPtr> search(unsigned int ilist,tcPDPtr id) const;
 
  bool allowed(long id1, long id2, long id3, long id4 = 0) const;
 
  VertexType::T getName() const { return _theName; }
 
  unsigned int getNpoint() const { return _npoint; }
 
  int orderInGem() const { return couplingOrders_.at(CouplingType::QED); }
 
  int orderInGs() const { return couplingOrders_.at(CouplingType::QCD); }
 
  int orderInCoupling(CouplingType::T cType) const {
    if(couplingOrders_.find(cType) !=couplingOrders_.end())
      return couplingOrders_.at(cType);
    else
      return 0;
  }
 
  int orderInAllCouplings() const {
    int output(0);
    for(auto & p : couplingOrders_)
      output += p.second;
    return output;
  }
 
  ColourStructure::T colourStructure() const {return colourStructure_;}
 
public:
 
  double strongCoupling(Energy2 q2) const {
    if(_coupopt==0) {
      double val = 4.0*Constants::pi*generator()->standardModel()->alphaS(q2);
      assert(val>=0.);
      return sqrt(val);
    }
    else if(_coupopt==1)
      return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaS());
    else
      return _gs;
  }
 
  double electroMagneticCoupling(Energy2 q2) const {
    if(_coupopt==0)
      return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMME(q2));
    else if(_coupopt==1)
      return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMMZ());
    else
      return _ee;
  }
 
  double weakCoupling(Energy2 q2) const {
    if( _coupopt == 0 )
      return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMME(q2)/
                  generator()->standardModel()->sin2ThetaW());
    else if( _coupopt == 1 )
      return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMMZ()/
                  generator()->standardModel()->sin2ThetaW());
    else
      return _ee/_sw;
  }
 
  double sin2ThetaW() const {
    if( _coupopt == 0 || _coupopt  == 1)
      return generator()->standardModel()->sin2ThetaW();
    else
      return sqr(_sw);
  }
 
public:
 
  virtual void setCoupling(Energy2 q2,tcPDPtr part1,
                           tcPDPtr part2,tcPDPtr part3)=0;
 
  virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3,
                           tcPDPtr part4)=0;
 
protected:
 
  virtual void doinit();
 
  virtual void rebind(const TranslationMap & trans);
 
  virtual IVector getReferences();
 
protected:
  void addToList(const vector<long> & ids);
 
  void addToList(long ida, long idb, long idc, long idd = 0);
 
protected:
  void norm(const Complex & coup) { _norm = coup; }
 
  virtual Complex propagator(int iopt, Energy2 q2,tcPDPtr part,
                             complex<Energy> mass=-GeV,
                             complex<Energy> width=-GeV);
 
  Complex normPropagator(int iopt, Energy2 q2,tcPDPtr part,
                         complex<Energy> mass=-GeV, 
                         complex<Energy> width=-GeV) {
    return _norm*propagator(iopt,q2,part,mass,width);
  }
 
public:
 
  bool kinematics() const { return _calckinematics; }
 
  void kinematics(bool kine ) { _calckinematics=kine; }
 
  void calculateKinematics(const Lorentz5Momentum & p0,
                           const Lorentz5Momentum & p1,
                           const Lorentz5Momentum & p2) {
    _kine[0][0]=p0*p0;
    _kine[1][1]=p1*p1;
    _kine[2][2]=p2*p2;
    _kine[0][1]=p0*p1;_kine[1][0]=_kine[0][1];
    _kine[0][2]=p0*p2;_kine[2][0]=_kine[0][2];
    _kine[1][2]=p1*p2;_kine[2][1]=_kine[1][2];
  }
  
  void calculateKinematics(const Lorentz5Momentum & p0,
                           const Lorentz5Momentum & p1,
                           const Lorentz5Momentum & p2,
                           const Lorentz5Momentum & p3) {
    _kine[0][0]=p0*p0;
    _kine[1][1]=p1*p1;
    _kine[2][2]=p2*p2;
    _kine[3][3]=p3*p3;
    _kine[0][1]=p0*p1;_kine[1][0]=_kine[0][1];
    _kine[0][2]=p0*p2;_kine[2][0]=_kine[0][2];
    _kine[0][3]=p0*p3;_kine[3][0]=_kine[0][3];
    _kine[1][2]=p1*p2;_kine[2][1]=_kine[1][2];
    _kine[1][3]=p1*p3;_kine[3][1]=_kine[1][3];
    _kine[2][3]=p2*p3;_kine[3][2]=_kine[2][3];
  }
  
  void calculateKinematics(const vector<Lorentz5Momentum> & p) {
    for(size_t ix=0;ix<p.size();++ix) {
      for(size_t iy=0;iy<=ix;++ix) {
        _kine[ix][iy]=p[ix]*p[iy];
        _kine[iy][ix]=_kine[ix][iy];
      }
    }
  }
 
  Energy2 invariant(unsigned int ix ,unsigned int iy) const {
    assert ( ix < _npoint && iy < _npoint );
    return _kine[ix][iy];
  }
  
protected:
 
  void orderInGem(int order) { couplingOrders_[CouplingType::QED] = order; }
 
  void orderInGs (int order) { couplingOrders_[CouplingType::QCD] = order; }
 
  void orderInCoupling(CouplingType::T cType, int order) {
    couplingOrders_[cType] = order;
  }
 
  void colourStructure(ColourStructure::T structure) {
    colourStructure_ = structure;
  }
  
private:
  
  VertexBase & operator=(const VertexBase &) = delete;
  
private:
 
  vector<vector<PDPtr> > _particles;
 
  unsigned int _npoint;
 
  set<tPDPtr> _inpart;
 
  set<tPDPtr> _outpart;
 
  Complex _norm;
 
  bool _calckinematics;
 
  std::array<std::array<Energy2,5>,5> _kine;
 
  VertexType::T _theName;
 
  ColourStructure::T colourStructure_;
 
  map<CouplingType::T,int> couplingOrders_;
 
  unsigned int _coupopt;
 
  double _gs;
 
  double _ee;
 
  double _sw;
};
  
ostream & operator<<(ostream &, const VertexBase &);
  
}
}
 
#endif /* ThePEG_VertexBase_H */