MEPP2GammaGammaPowheg.h

// -*- C++ -*-
#ifndef HERWIG_MEPP2GammaGammaPowheg_H
#define HERWIG_MEPP2GammaGammaPowheg_H
//
// This is the declaration of the MEPP2GammaGammaPowheg class.
//
 
#include "Herwig/MatrixElement/HwMEBase.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
#include "Herwig/Shower/ShowerAlpha.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class MEPP2GammaGammaPowheg: public HwMEBase {
 
  enum DipoleType {           IIQCD1=2,           IIQCD2=4, 
                    IFQED1=5, FIQED1=6, IFQED2=7, FIQED2=8 };
 
  enum RadiationType {Subtraction,Hard,Shower};
 
public:
 
  MEPP2GammaGammaPowheg();
 
  virtual POWHEGType hasPOWHEGCorrection() {return ISR;}
 
  virtual RealEmissionProcessPtr generateHardest(RealEmissionProcessPtr,
                                                 ShowerInteraction);
 
public:
 
  virtual unsigned int orderInAlphaS() const;
 
  virtual unsigned int orderInAlphaEW() const;
 
  virtual double me2() const;
 
  virtual Energy2 scale() const;
 
  virtual int nDim() const;
 
  virtual bool generateKinematics(const double * r);
 
  virtual CrossSection dSigHatDR() 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:
 
  virtual double NLOWeight() const;
 
  double loGammaGammaME(const cPDVector & particles,
                        const vector<Lorentz5Momentum> & momenta,
                        bool first=false) const;
 
  double loGammaqME(const cPDVector & particles,
                const vector<Lorentz5Momentum> & momenta,
                bool first=false) const;
 
  double loGammaqbarME(const cPDVector & particles,
                       const vector<Lorentz5Momentum> & momenta,
                       bool first=false) const;
  
  double loGammagME(const cPDVector & particles,
                    const vector<Lorentz5Momentum> & momenta,
                    bool first=false) const;
 
  InvEnergy2 realGammaGammagME(const cPDVector & particles,
                               const vector<Lorentz5Momentum> & momenta,
                               DipoleType dipole, RadiationType rad,
                               bool first=false) const;
  
  InvEnergy2 realGammaGammaqME(const cPDVector & particles,
                               const vector<Lorentz5Momentum> & momenta,
                               DipoleType dipole, RadiationType rad,
                               bool first=false) const;
  
  InvEnergy2 realGammaGammaqbarME(const cPDVector & particles,
                                  const vector<Lorentz5Momentum> & momenta,
                                  DipoleType dipole, RadiationType rad,
                                  bool first=false) const;
  
  double subtractedVirtual() const;
 
  double subtractedReal(pair<double,double> x, double z,
                        double zJac, double oldqPDF, double newqPDF,
                        double newgPDF,bool order) const;
 
  double collinearQuark(double x, Energy2 mu2, double jac, double z,
                        double oldPDF, double newPDF) const;
 
  double collinearGluon(Energy2 mu2, double jac, double z,
                        double oldPDF, double newPDF) const;
 
  double realME(const cPDVector & particles,
                const vector<Lorentz5Momentum> & momenta) const;
 
  RealEmissionProcessPtr hardQCDEmission(RealEmissionProcessPtr,
                                         ParticleVector,
                                         pair<double,double>);
 
  RealEmissionProcessPtr hardQEDEmission(RealEmissionProcessPtr,
                                         ParticleVector,
                                         pair<double,double>);
 
  pair<double,double> supressionFunction(Energy pT,Energy scale) const {
    if(supressionScale_==0) scale = lambda_;
    Energy2 scale2 = sqr(scale), pT2 = sqr(pT);
    switch( supressionFunction_ ) {
    case 0:
      return make_pair(1.,0.);
    case 1:
      if(pT < scale ) return make_pair(1.,0.);
      else            return make_pair(0.,1.);
    case 2:
      return make_pair(scale2/(pT2+scale2),pT2/(pT2+scale2));
    default:
      assert(false);
      return make_pair(0.,0.);
    }
  }
 
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
protected:
 
  virtual void doinit();
 
private:
 
  MEPP2GammaGammaPowheg & operator=(const MEPP2GammaGammaPowheg &) = delete;
 
private:
 
  AbstractFFVVertexPtr FFPvertex_;
 
  AbstractFFVVertexPtr FFGvertex_;
 
  mutable double zTilde_;
 
  mutable double vTilde_;
 
  mutable double phi_;
 
  unsigned int contrib_;
 
  double power_;
 
  tcPDPtr gluon_;
 
  unsigned int process_;
 
  unsigned int threeBodyProcess_;
 
  int maxflavour_;
 
  mutable double CFfact_;
 
  mutable double TRfact_;
 
  mutable double alphaS_;
 
  bool fixedAlphaS_;
 
  mutable double alphaEM_;
 
  mutable double loME_;
 
  mutable ProductionMatrixElement me_;
 
  mutable DipoleType dipole_;
 
  unsigned int supressionFunction_;
 
  unsigned int supressionScale_;
 
  Energy lambda_;
 
 
  bool quarkplus_;
 
  vector<tcBeamPtr> beams_;
  
  vector<tcPDPtr> partons_;
 
  double preQCDqqbarq_;
  double preQCDqqbarqbar_;
 
  double preQCDqg_;
 
  double preQCDgqbar_;
 
  double preQEDqqbarq_;
  double preQEDqqbarqbar_;
  double preQEDqgq_;
  double preQEDgqbarqbar_;
 
  vector<double> prefactor_;
 
  Energy minpT_;
 
  ShowerAlphaPtr alphaQCD_;
 
  ShowerAlphaPtr alphaQED_;
 
  unsigned int scaleChoice_;
 
  double scalePreFactor_;
 
};
 
}
 
#endif /* HERWIG_MEPP2GammaGammaPowheg_H */