DecayIntegrator.h

// -*- C++ -*-
//
// DecayIntegrator.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_DecayIntegrator_H
#define Herwig_DecayIntegrator_H
//
// This is the declaration of the DecayIntegrator class.
//
 
#include "DecayIntegrator.fh"
#include "HwDecayerBase.h"
#include "PhaseSpaceMode.fh"
#include "Herwig/PDT/WidthCalculatorBase.fh"
#include "Radiation/DecayRadiationGenerator.h"
#include <Herwig/Decay/DecayVertex.h>
 
namespace Herwig {
 
using namespace ThePEG;
 
class DecayIntegrator: public HwDecayerBase {
 
public:
  
  friend class PhaseSpaceMode;
 
  enum MEOption {Initialize,Calculate,Terminate};
 
public:
 
  DecayIntegrator() : nIter_(10), nPoint_(10000), nTry_(500),
                       generateInter_(false), iMode_(-1),
                       realME_(false), virtualME_(false), eps_(ZERO), warnings_(false)
  {}
 
public:
  
  virtual bool accept(tcPDPtr parent, const tPDVector & children) const {
    bool cc;
    return modeNumber(cc,parent,children)>=0;
  }
  
  virtual ParticleVector decay(const Particle & parent,
                               const tPDVector & children) const;
  
  virtual int modeNumber(bool & cc, tcPDPtr parent, 
                         const tPDVector & children) const = 0;
 
  int imode() const {return iMode_;}
 
  void addMode(PhaseSpaceModePtr mode) const;
 
  virtual double me2(const int ichan, const Particle & part,
                     const tPDVector & outgoing,
                     const vector<Lorentz5Momentum> & momenta,
                     MEOption meopt) const = 0;
 
  virtual void constructSpinInfo(const Particle & part,
                                 ParticleVector outgoing) const = 0;
  
  virtual void dataBaseOutput(ofstream & os,bool header) const;
  
  void setPartialWidth(const DecayMode & dm, int imode);
  virtual bool twoBodyMEcode(const DecayMode &, int & mecode,
                             double & coupling) const {
    coupling = 1.;
    mecode   = -1;
    return false;
  }                                  
  
  virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
  
  virtual double threeBodyMatrixElement(const int imode,  const Energy2 q2,
                                        const Energy2 s3, const Energy2 s2, 
                                        const Energy2 s1, const Energy  m1, 
                                        const Energy  m2, const Energy  m3) const;
  
  virtual InvEnergy threeBodydGammads(const int imode, const Energy2 q2,
                                      const Energy2 s,
                                      const Energy m1, const Energy m2, 
                                      const Energy m3) const;
 
  int findMode(const DecayMode & dm);
  
public:
 
  ParticleVector generatePhotons(const Particle & p,ParticleVector children) {
    return photonGen_->generatePhotons(p,children,this);
  }
 
  bool canGeneratePhotons() {return photonGen_;}
 
  virtual double oneLoopVirtualME(unsigned int imode,
                                  const Particle & part, 
                                  const ParticleVector & products);
 
  bool hasOneLoopME() {return virtualME_;}
  
  virtual InvEnergy2 realEmissionME(unsigned int imode,
                                    const Particle & part, 
                                    ParticleVector & products,
                                    unsigned int iemitter,
                                    double ctheta, double stheta,
                                    const LorentzRotation & rot1,
                                    const LorentzRotation & rot2);
 
  bool hasRealEmissionME() {return realME_;}
  
public:
 
  friend ostream & operator<<(ostream & os, const DecayIntegrator & decay);
 
  
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
  
protected:
 
  virtual void doinitrun();
 
protected:
 
  ParticleVector generate(bool inter,bool cc, const unsigned int & imode,
                          const Particle & inpart) const;
 
  void imode(int in) { iMode_ = in;}
  
  void ME(DecayMEPtr in) const { matrixElement_ = in;}
  
  DecayMEPtr ME() const {return matrixElement_;}
   
  void resetIntermediate(tcPDPtr part, Energy mass, Energy width);
 
  Energy initializePhaseSpaceMode(unsigned int imode,bool init, bool onShell=false) const;
  
protected:
  
  void generateIntermediates(bool in) {generateInter_=in;}
 
  bool generateIntermediates() const {return generateInter_;}
  
  void hasOneLoopME(bool in) {virtualME_=in;}
 
  void hasRealEmissionME(bool in) {realME_=in;}
 
  void epsilonPS(Energy in) {eps_=in;}
 
  void clearModes() {modes_.clear();}
 
protected:
  
  unsigned int numberModes() const {return modes_.size();}
 
  tPhaseSpaceModePtr mode(unsigned int ix) {
    return modes_[ix];
  }
  tcPhaseSpaceModePtr mode(unsigned int ix) const {
    return modes_[ix];
  }
  
public:
 
bool warnings() const {
  return warnings_;
}
 
private:
 
  DecayIntegrator & operator=(const DecayIntegrator &) = delete;
 
  unsigned int nIter_;
 
  unsigned int nPoint_;
 
  unsigned int nTry_;
  
  mutable vector<PhaseSpaceModePtr> modes_;
 
  
  bool generateInter_;
  
  DecayRadiationGeneratorPtr photonGen_;
 
  mutable int iMode_;
 
  mutable DecayMEPtr matrixElement_;
 
  bool realME_;
 
  bool virtualME_;
 
  Energy eps_;
 
  bool warnings_;
 
protected:
 
  class DecayIntegratorError: public Exception {};
 
};
  
ostream & operator<<(ostream &, const DecayIntegrator &);
 
}
 
#endif /* Herwig_DecayIntegrator_H */