MEDiffraction.h

// -*- C++ -*-
#ifndef HERWIG_MEDiffraction_H
#define HERWIG_MEDiffraction_H
//
// This is the declaration of the MEDiffraction class.
//
 
#include "Herwig/MatrixElement/HwMEBase.h"
#include "Herwig/Shower/UEBase.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class MEDiffraction: public HwMEBase {
 
public:
 
  MEDiffraction();
 
  virtual unsigned int orderInAlphaS() const;
 
  virtual unsigned int orderInAlphaEW() const;
 
  virtual double me2() const;
 
  double correctionweight() const;
 
  virtual Energy2 scale() const;
 
  virtual void setKinematics();
 
  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;
 
  /* virtual bool wantCMS() const { return false; } */
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual void doinit();
 
  virtual void doinitrun();
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
private:
 
  /* The matrix element squared */
  double theme2;
 
  /* Use only delta as excited state */
  bool deltaOnly;
 
  /* Only produce single-diffractive events */
  // Default = off
  bool singleDiffractionOnly = 0;
 
  /* Direction of the excited proton */
  unsigned int diffDirection;
 
  /* Number of clusters the dissociated proton decays into */
  unsigned int dissociationDecay;
 
  /* The mass of the consitutent quark */
  Energy mq() const {return Energy(0.325*GeV);}
 
  /* The mass of the constituent diquark */
  Energy mqq() const {return Energy(0.650*GeV);}
 
  /* The proton-pomeron slope */
  double theprotonPomeronSlope;
 
  /* The soft pomeron intercept */
  double thesoftPomeronIntercept;
 
  /* The soft pomeron slope */
  double thesoftPomeronSlope;
 
 
  pair<pair<Energy2,Energy2>,Energy2> diffractiveMassAndMomentumTransfer() const;
 
 
  Energy2 randomM2() const;
 
  Energy2 randomt(Energy2 M2) const;
 
  Energy2 doublediffrandomt(Energy2 M12, Energy2 M22) const;
 
 
  pair<Lorentz5Momentum,Lorentz5Momentum> twoBodyDecayMomenta(Lorentz5Momentum pp) const;
 
  InvEnergy2 protonPomeronSlope() const;
 
  double softPomeronIntercept() const;
 
  //M12 and M22 are masses squared of
  //outgoing particles
 
  Energy2 tminfun(Energy2 s, Energy2 M12, Energy2 M22) const;
 
  Energy2 tmaxfun(Energy2 s , Energy2 M12, Energy2 M22) const;
 
  //lowest possible mass given the constituent masses of quark and diquark
  Energy2 M2min() const{return sqr(getParticleData(2212)->mass()+mq()+mqq());}
 
  Energy2 M2max() const{
    return sqr(generator()->maximumCMEnergy()-getParticleData(2212)->mass());
  }//TODO:modify to get proper parameters
 
  InvEnergy2 softPomeronSlope() const;
 
 
 
  /* Kallen function */
  template<typename A, typename B, typename C>
  auto kallen(A a, B b, C c) const -> decltype(a*a)
  {
    return a*a + b*b + c*c - 2.0*(a*b + b*c + c*a);
  }
 
 
 
  MEDiffraction & operator=(const MEDiffraction &) = delete;
 
  bool isInRunPhase;
 
 
  /* The proton mass */
  Energy theProtonMass;
 
  UEBasePtr MPIHandler_;
 
};
 
}
 
#endif /* HERWIG_MEDiffraction_H */