MEHiggsPairOL.h

// -*- C++ -*-
//
// MEHiggsPair.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2009-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_MEHiggsPairOL_H
#define HERWIG_MEHiggsPairOL_H
//
// This is the declaration of the MEHiggsPairOL class.
//
// The implementation of this process is based upon hep-ph/0112161 by G.F. Giudice, R. Rattazzi, J.D. Wells.
 
#include "Herwig/MatrixElement/HwMEBase.h"
#include "ThePEG/Repository/UseRandom.h"
#include "Herwig/PDT/GenericMassGenerator.h"
#include "Herwig/Utilities/Kinematics.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
#include "HiggsPair.h"
 
 
namespace Herwig {
using namespace ThePEG;
 
 
/* 
 * Interface to external FORTRAN hpair functions (T.Plehn, M.Spira & P.Zerwas, arXiv:hep-ph/9603205)
 * original Program found at http://people.web.psi.ch/spira/hpair/ 
 */
extern "C" { 
  complex<double> eta_(complex<double>* C1, complex<double>* C2);
  complex<double> d04_(double* P1, double* P2, double* P3, double* P4, double* P12,double* P23, double* M1, double* M2, double* M3, double* M4);
  complex<double> c03_(double* P1,double* P2,double* P3,double* M1,double* M2,double* M3); 
  complex<double> cspen_(double* Z);
  complex<double> sqe_(double* A, double* B, double* C);
  complex<double> etas_(complex<double>*Y,complex<double>* R,complex<double>*RS);
  void formfac_(double* AMQ, double* S, double* T, double* U, double* M1, double* M2, complex<double>* C0AB,complex<double>* C0AC,complex<double>* C0AD,complex<double>* C0BC,complex<double>* C0BD,complex<double>* C0CD,complex<double>* D0ABC,complex<double>* D0BAC,complex<double>* D0ACB);
}
 
extern "C" {
  void ol_setparameter_int(const char* param, int val);
  void ol_setparameter_double(const char* param, double val);
  void ol_setparameter_string(const char* param, char* val);
  int ol_register_process(const char* process, int amptype);
  int ol_n_external(int id);
  void ol_phase_space_point(int id, double sqrt_s, double* pp);
  void ol_start();
  void ol_finish();
  void ol_evaluate_tree(int id, double* pp, double* m2_tree);
  void ol_evaluate_loop(int id, double* pp, double* m2_tree, double* m2_loop, double* acc);
  void ol_evaluate_loop2(int id, double* pp, double* m2_loop2, double* acc);
}
 
extern "C" { 
  extern struct{
      complex<double> A1;
      complex<double> A2;
      complex<double> H1;
      complex<double> H2;
      complex<double> Z1;
      complex<double> Z2; 
      complex<double> AA1;
      complex<double> AA2;
      complex<double> HH1;
      complex<double> HH2; 
      complex<double> AH1;
      complex<double> AH2;
  } form_;
}   
 
class MEHiggsPairOL: public HwMEBase {
 
 
public:
 
  MEHiggsPairOL();
 
  virtual unsigned int orderInAlphaS() const { return 0; }
 
  virtual unsigned int orderInAlphaEW() const { return 0; }
 
  virtual Energy2 scale() const;
 
  virtual void getDiagrams() const;
 
  virtual bool generateKinematics(const double * r);
 
 
  virtual CrossSection dSigHatDR() const;
 
 
  
 
  virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
 
  virtual Selector<const ColourLines *>
  colourGeometries(tcDiagPtr diag) const;
 
  virtual double me2() const;
 
 
  virtual void setKinematics();
 
  virtual int nDim() const;
 
 
protected:
 
  virtual void doinit();
 
  virtual void doinitrun();
 
 
  virtual void rebind(const TranslationMap & trans);
 
  virtual IVector getReferences();
 
  PDPtr higgs() const { return _higgs; }
 
  void higgs(PDPtr in) {_higgs =in;}
 
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
 
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
protected:
 
private:
  static ClassDescription<MEHiggsPairOL> initMEHiggsPairOL;
 
  MEHiggsPairOL & operator=(const MEHiggsPairOL &) = delete;
 
 
  int _id, _idtri, _idbox, _idint;
 
  PDPtr _higgs;
 
  /* 
   * The top mass
   */
  Energy _topmass;
 
  /* 
   * The bottom mass
   */
 Energy _bottommass;
 
  /*
   * The Z boson mass
   */ 
  Energy _zmass;
  
 
  /*                                           
   * Higgs boson mass(es) 
   */
  Energy _m1, _m2;
 
  /*                                           
   * Heavy H mass
   */
  Energy _heavyHmass;
 
 
  /*                                           
   * Heavy H width
   */
  Energy _heavyHwidth;
 
 
 
 
private:
  tcSMPtr _theModel;
 
  GenericMassGeneratorPtr _hmass;
 
  double _selfcoupling;
 
 
  unsigned int _process;
 
  Energy _mh;
 
 
  Energy _wh;
 
  /*
   * Fix alphaS
   */ 
  unsigned int _fixedalphaS;
 
 /*
   * Value of AlphaS if fixed using option 2 above
   */ 
  double _alphasfixedvalue;
 
  /*
   * Fix scale of whole process
   */ 
  unsigned int _fixedscale;
  
  /*
   * Scale to use for alpha_S if fixed
   */
  Energy _alphascale;
 
  /*
   * Choose between the OpenLoops and HPAIR implementations
   */ 
  unsigned int _implementation;
 
  /*
   * include the widths
   */
  unsigned int _includeWidths;
 
  /*
   * include the b quark loops
   */
  unsigned int _includeBquark;
 
 
  /* 
   * Base scale to use if chosen to be fixed
   */ 
  Energy _basescale;
 
  /* 
   * scale multiplier 
   */ 
  double _scalemultiplier;
 
 
  /* 
   * Masses of fermions 
   */
  double Mass_E;
  double Mass_M;
  double Mass_L;
  double Mass_T;
  double Mass_U;
  double Mass_C;
  double Mass_D;
  double Mass_S;
  double Mass_B;
  double Mass_Z;
  double Mass_W;
  double Mass_H;
  double Width_C;
  double Width_B;
  double Width_T;
  double Width_W;
  double Width_Z;
  double Width_H;
 
 
  double getMass_E() const { return Mass_E; }
  double getMass_M() const { return Mass_M; }
  double getMass_L() const { return Mass_L; }
 
  double getMass_T() const { return Mass_T; }
  double getMass_U() const { return Mass_U; }  
  double getMass_C() const { return Mass_C; }
  double getMass_D() const { return Mass_D; }
  double getMass_S() const { return Mass_S; }
  double getMass_B() const { return Mass_B; }
 
  /*
   * Masses of bosons
   */
  double getMass_Z() const { return Mass_Z; }
  double getMass_W() const { return Mass_W; }
  double getMass_H() const { return Mass_H; }
 
  /*
   * Widths of particles
   */ 
  double getWidth_C() const { return Width_C; }
  double getWidth_B() const { return Width_B; }
  double getWidth_T() const { return Width_T; }
  double getWidth_W() const { return Width_W; }
  double getWidth_Z() const { return Width_Z; }
  double getWidth_H() const { return Width_H; }
 
 
  /*
   * QCD and QED couplings 
   */
  double Coupl_Alpha_QED; 
  double Coupl_Alpha_QCD; 
 
  
  /* Scalar integral initialization from hpair.f (T.Plehn, M.Spira & P.Zerwas, arXiv:hep-ph/9603205)
     Program found at http://people.web.psi.ch/spira/hpair/ */
  
  virtual vector<Complex> iniscal(double AMQ, double S, double T,double U, double M1, double M2) const;
  
  /* Matrix element calculation from hpair */
 
  virtual double MATRIX(double S, double T,double U, double M1, double M2) const;
 
 
 /* OpenLoops initialization */
  virtual void InitOpenLoops();
 
};
 
}
 
#include "ThePEG/Utilities/ClassTraits.h"
 
namespace ThePEG {
 
template <>
struct BaseClassTrait<Herwig::MEHiggsPairOL,1> {
  typedef Herwig::HwMEBase NthBase;
 
};
 
template <>
struct ClassTraits<Herwig::MEHiggsPairOL>
  : public ClassTraitsBase<Herwig::MEHiggsPairOL> {
  static string className() { return "Herwig::MEHiggsPairOL"; }
  static string library() { return "HiggsPair.so"; }
};
 
}
 
 
#endif /* HERWIG_MEHiggsPairOL_H */