MEQCD2to2Fast.h

// -*- C++ -*-
//
// MEQCD2to2Fast.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MEQCD2to2Fast_H
#define HERWIG_MEQCD2to2Fast_H
//
// This is the declaration of the MEQCD2to2Fast class.
//
 
#include "Herwig/MatrixElement/HwMEBase.h"
#include "ThePEG/Repository/UseRandom.h"
 
namespace Herwig {
using namespace ThePEG;
using namespace ThePEG::Helicity;
 
class MEQCD2to2Fast: public HwMEBase {
 
public:
 
  MEQCD2to2Fast() :_maxflavour(5),_process(0),_strictFlavourScheme(false) {
    massOption(vector<unsigned int>(2,0));
  }
 
  virtual unsigned int orderInAlphaS() const;
 
  virtual unsigned int orderInAlphaEW() const;
 
  virtual double me2() const;
 
  virtual Energy2 scale() const;
 
  virtual void getDiagrams() const;
 
  virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
 
  virtual Selector<const ColourLines *>
  colourGeometries(tcDiagPtr diag) const;
 
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  double gg2ggME() const {
    Energy2 u(uHat()),t(tHat()),s(sHat());
    double output = 9./4.*(3.-t*u/s/s-s*u/t/t-s*t/u/u);
    double flow[3]={(1.-u*t/s/s-s*t/u/u+t*t/s/u),
                    (1.-t*u/s/s-s*u/t/t+u*u/s/t),
                    (1.-t*s/u/u-u*s/t/t+s*s/u/t)};
    _flow = 1+UseRandom::rnd3(flow[0],flow[1],flow[2]);
    double diag[3]={(sqr(u)+sqr(t))/sqr(s),
                    (sqr(s)+sqr(u))/sqr(t),
                    (sqr(s)+sqr(t))/sqr(u)};
    if(_flow==1)      diag[1]=0;
    else if(_flow==2) diag[2]=0;
    else if(_flow==3) diag[0]=0;
    _diagram=1+UseRandom::rnd3(diag[0],diag[1],diag[2]);
    return output;
  }
 
  double gg2qqbarME() const {
    Energy2 u(uHat()),t(tHat()),s(sHat());
    Energy4 u2(sqr(u)),t2(sqr(t)),s2(sqr(s));
    double output =(1./6./u/t-3./8./s2)*(t2+u2);
    double flow[2]={u2/(u2+t2),t2/(u2+t2)};
    _flow = 1+UseRandom::rnd2(flow[0],flow[1]);
    _diagram=3+_flow;
    return output;
  }
 
  double qqbar2ggME() const {
    Energy2 u(uHat()),t(tHat()),s(sHat());
    Energy4 s2(sqr(s)),u2(sqr(u)),t2(sqr(t));
    double output = 0.5*(32./27./u/t-8./3./s2)*(t2+u2);
    double flow[2] = {u2/(u2+t2),t2/(t2+u2)};
    _flow=1+UseRandom::rnd2(flow[0],flow[1]);
    _diagram=6+_flow;
    return output;
  }
 
  double qg2qgME() const {
    Energy2 u(uHat()),t(tHat()),s(sHat());
    Energy4 s2(sqr(s)),u2(sqr(u)),t2(sqr(t));
    double output = (-4./9./s/u+1./t2)*(s2+u2);
    double flow[2]={u2/(s2+u2),s2/(s2+u2)};
    _flow=1+UseRandom::rnd2(flow[0],flow[1]);
    _diagram=9+_flow;
    return output;
  }
  
  double qbarg2qbargME() const {
    // scale
    Energy2 u(uHat()),t(tHat()),s(sHat());
    Energy4 u2(sqr(u)),s2(sqr(s)); // t2(sqr(t))
    double flow[2]={u2/(s2+u2),s2/(s2+u2)};
    _flow=1+UseRandom::rnd2(flow[0],flow[1]);
    _diagram=12+_flow;
    return (-4./9./s/u+1./t/t)*(s*s+u*u);
  }
  
  double qq2qqME() const {
    Energy2 u(uHat()),t(tHat());
    Energy4 s2(sqr(sHat())),u2(sqr(u)),t2(sqr(t));
    double output;
    if(mePartonData()[0]->id()==mePartonData()[1]->id()) {
      output = 0.5*(4./9.*((s2+u2)/t2+(s2+t2)/u2)
                    -8./27.*s2/u/t);
      double flow[2]={(s2+u2)/t2,(s2+t2)/u2}; 
      _flow=1+UseRandom::rnd2(flow[0],flow[1]);
    }
    else {
      output = 4./9.*(s2+u2)/t2;
      _flow=2;
    }
    _diagram = 15+_flow;
    return output;
  }
 
  double qbarqbar2qbarqbarME() const {
    Energy2 u(uHat()),t(tHat());
    Energy4 u2(sqr(u)),t2(sqr(t)),s2(sqr(sHat()));
    double output;
    if(mePartonData()[0]->id()==mePartonData()[1]->id()) {
      output = 0.5*(4./9.*((s2+u2)/t2+(s2+t2)/u2)
                    -8./27.*s2/u/t);
      double flow[2]={(s2+u2)/t2,(s2+t2)/u2};
      _flow=1+UseRandom::rnd2(flow[0],flow[1]);
    }
    else {
      output = 4./9.*(s2+u2)/t2;
      _flow = 2;
    }
    _diagram = 17+_flow;
    // final part of colour and spin factors
    return output;
  }
  
  double qqbar2qqbarME() const {
    // type of process
    bool diagon[2]={mePartonData()[0]->id()== -mePartonData()[1]->id(),
                    mePartonData()[0]->id()==  mePartonData()[2]->id()};
    // scale
    Energy2 u(uHat()),t(tHat()),s(sHat());
    Energy4 s2(sqr(s)),t2(sqr(t)),u2(sqr(u));
    double output;
    if(diagon[0]&&diagon[1]) {
      output= (4./9.*((s2+u2)/t2+(u2+t2)/s2)
               -8./27.*u2/s/t);
      double flow[2]={(t2+u2)/s2,(s2+u2)/t2};
      _flow=1+UseRandom::rnd2(flow[0],flow[1]);
    }
    else if(diagon[0]) {
      output = (4./9.*(t2+u2)/s2);
      _flow=1;
    }
    else {
      output = (4./9.*(s2+u2)/t2);
      _flow=2;
    }
    _diagram=19+_flow;
    return output;
  }
  
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
protected:
 
  virtual void doinit();
 
 
private:
 
  MEQCD2to2Fast & operator=(const MEQCD2to2Fast &) = delete;
 
private:
 
  unsigned int _maxflavour;
 
  unsigned int _process;
 
  mutable unsigned int _flow;
 
  mutable unsigned int _diagram;
 
  bool _strictFlavourScheme;
 
};
 
}
 
#endif /* HERWIG_MEQCD2to2Fast_H */