SMZDecayer.h

// -*- C++ -*-
//
// SMZDecayer.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_SMZDecayer_H
#define HERWIG_SMZDecayer_H
//
// This is the declaration of the SMZDecayer class.
//
#include "Herwig/Decay/PerturbativeDecayer.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "Herwig/Decay/PhaseSpaceMode.h"
 
namespace Herwig {
using namespace ThePEG;
using namespace ThePEG::Helicity;
 
class SMZDecayer: public PerturbativeDecayer {
 
public:
 
  SMZDecayer();
 
  virtual bool hasMECorrection() {return true;}
 
  virtual void initializeMECorrection(RealEmissionProcessPtr , double & ,
                                      double & );
 
  virtual bool softMatrixElementVeto(PPtr parent,
                                     PPtr progenitor,
                                     const bool & fs,
                                     const Energy & highestpT,
                                     const tcPDVector & ids,
                                     const double & z,
                                     const Energy & scale,
                                     const Energy & pT);
    
  virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
 
public:
 
  virtual int modeNumber(bool & cc, tcPDPtr parent, 
                         const tPDVector & children) const;
 
  virtual ParticleVector decay(const Particle & parent,const tPDVector & children) const;
 
  double me2(const int ichan,const Particle & part,
             const tPDVector & outgoing,
             const vector<Lorentz5Momentum> & momenta,
             MEOption meopt) const;
 
  virtual void constructSpinInfo(const Particle & part,
                                 ParticleVector outgoing) const;
 
  virtual void dataBaseOutput(ofstream & os,bool header) const;
 
  virtual double oneLoopVirtualME(unsigned int imode,
                                  const Particle & part, 
                                  const ParticleVector & products);
  
  virtual InvEnergy2 realEmissionME(unsigned int imode,
                                    const Particle & part, 
                                    ParticleVector & products,
                                    unsigned int iemitter,
                                    double ctheta, double stheta,
                                    const LorentzRotation & rot1,
                                    const LorentzRotation & rot2);
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
  
  static void Init();
 
protected:
  
  virtual IBPtr clone() const {return new_ptr(*this);}
 
  virtual IBPtr fullclone() const {return new_ptr(*this);}
  
protected:
  
  virtual void doinit();
 
  virtual void doinitrun();
 
protected:
 
  void setRho(double);
 
  void setKtildeSymm();
 
  void setKtilde2();
 
  double getZfromX(double, double);
 
  double getKfromX(double, double);
 
  void getXXbar(double kt, double z, double & x, double & xbar);
 
  double qWeight(double x, double xbar); 
 
  double qbarWeight(double x, double xbar);
 
  double qWeightX(Energy qtilde, double z);
 
  double qbarWeightX(Energy qtilde, double z);
 
  double u(double);
 
  double MEV(double, double);
 
  double MEA(double, double);
 
  double PS(double x1, double x2);
  
protected:
 
  double calculateRealEmission(double x1, double x2, 
                               tcPDVector outgoing,
                               vector<Lorentz5Momentum> momenta,
                               double phi,
                               bool subtract,
                               int emitter) const;
 
  double calculateRealEmission(double x1, double x2, 
                               vector<PPtr> hardProcess,
                               double phi,
                               bool subtract) const;
 
  bool checkZMomenta(double x1, double x2, double x3, double y, Energy pT) const;
 
  InvEnergy calculateJacobian(double x1, double x2, Energy pT) const;
 
  double meRatio(tcPDVector partons, 
                 vector<Lorentz5Momentum> momenta,
                 unsigned int iemitter,bool subtract) const;
 
  virtual double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
                                    const ParticleVector & decay3, MEOption meopt,
                                    ShowerInteraction inter);
  
  double loME(const tcPDVector & partons, 
              const vector<Lorentz5Momentum> & momenta) const;
 
  InvEnergy2 realME(const tcPDVector & partons, 
                    const vector<Lorentz5Momentum> & momenta,
                    ShowerInteraction inter) const;
 
  bool getEvent(vector<PPtr> hardProcess);
 
private:
 
  SMZDecayer & operator=(const SMZDecayer &) = delete;
 
private:
 
 
 
  FFVVertexPtr FFZVertex_;
  
  AbstractFFVVertexPtr FFPVertex_;
 
  AbstractFFVVertexPtr FFGVertex_;
 
  vector<double> quarkWeight_;
 
  vector<double> leptonWeight_;
 
  mutable RhoDMatrix _rho;
 
  mutable vector<VectorWaveFunction> _vectors;
 
  mutable vector<SpinorWaveFunction> _wave;
 
  mutable vector<SpinorBarWaveFunction> _wavebar;
 
private:
 
  Energy d_Q_;
 
  Energy d_m_;
 
  double d_rho_;
 
  double d_v_;
 
  double d_kt1_;
 
  double d_kt2_;
 
  static const double EPS_;
 
private:
 
  double CF_;
 
  mutable Energy mZ_;
 
  mutable double mu_;
 
  mutable double mu2_;
 
  mutable double aS_;
 
  mutable Energy2 scale_;
 
  tcPDPtr gluon_;
 
  tcPDVector partons_;
 
  vector<Lorentz5Momentum> quark_;
 
  Lorentz5Momentum gauge_;
 
  PPtr zboson_;
 
  Energy2 mz2_;
 
  bool NLO_;
};
 
}
 
 
#endif /* HERWIG_SMZDecayer_H */