GeneralTwoBodyDecayer.h

// -*- C++ -*-
//
// GeneralTwoBodyDecayer.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_GeneralTwoBodyDecayer_H
#define HERWIG_GeneralTwoBodyDecayer_H
//
// This is the declaration of the GeneralTwoBodyDecayer class.
//
 
#include "Herwig/Decay/PerturbativeDecayer.h"
#include "Herwig/Decay/PhaseSpaceMode.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "GeneralTwoBodyDecayer.fh"
 
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
 
class GeneralTwoBodyDecayer: public PerturbativeDecayer {
 
public:
  
  typedef pair<tcPDPtr, Energy> PMPair;
 
public:
 
  GeneralTwoBodyDecayer() : maxWeight_(1.), colour_(1,DVector(1,1.))
  {}
 
 
  virtual ParticleVector decay(const Particle & parent,
                               const tPDVector & children) const;
 
  virtual int modeNumber(bool & cc, tcPDPtr parent,const tPDVector & children) const;
  
  virtual Energy partialWidth(PMPair inpart, PMPair outa, 
                              PMPair outb) const;
 
  virtual bool twoBodyMEcode(const DecayMode & dm, int & mecode,
                             double & coupling) const;
 
  virtual double brat(const DecayMode & dm, const Particle & p,
                      double oldbrat) const;
 
  virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
                            vector<VertexBasePtr>,
                            map<ShowerInteraction,VertexBasePtr> &,
                            const vector<map<ShowerInteraction,VertexBasePtr> > &,
                            map<ShowerInteraction,VertexBasePtr>) =0;
 
protected:
  
  void setWeight(double wgt) { maxWeight_ = wgt; }
 
  void colourConnections(const Particle & parent, 
                         const ParticleVector & out) const;
 
  double colourFactor(tcPDPtr in, tcPDPtr out1, tcPDPtr out2) const;
 
  double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
                            const ParticleVector & decay3, MEOption meopt,
                            ShowerInteraction inter);
 
  void decayInfo(PDPtr incoming, PDPair outgoing);
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual void doinit();
 
  virtual void doinitrun();
 
protected:
 
  typedef vector<pair<int,double > > CFlowPairVec;
  typedef vector<CFlowPairVec> CFlow;
 
  const vector<DVector> & getColourFactors(const Particle & inpart, 
                                           const ParticleVector & decay, 
                                           unsigned int & nflow); 
 
  const CFlow & colourFlows(const Particle & inpart,
                            const ParticleVector & decay);
 
  virtual double threeBodyME(const int , const Particle & inpart,
                             const ParticleVector & decay,
                             ShowerInteraction inter, MEOption meopt);
 
 
private:
 
  GeneralTwoBodyDecayer & operator=(const GeneralTwoBodyDecayer &) = delete;
 
private:
 
  PDPtr incoming_;
 
  vector<PDPtr> outgoing_;
 
  double maxWeight_;
 
  vector<DVector> colour_;
 
};
 
 
template<typename T, typename Cmp, typename A>
inline PersistentOStream & operator<<(PersistentOStream & os,
                                      const map<ShowerInteraction,T,Cmp,A> & m) {
  os << m.size();
  if(m.find(ShowerInteraction::QCD)!=m.end()) {
    os << 0 << m.at(ShowerInteraction::QCD);
  }
  if(m.find(ShowerInteraction::QED)!=m.end()) {
    os << 1 << m.at(ShowerInteraction::QED);
  }
  return os;
}
 
template <typename T, typename Cmp, typename A>
inline PersistentIStream & operator>>(PersistentIStream & is, map<ShowerInteraction,T,Cmp,A> & m) {
  m.clear();
  long size;
  int k;
  is >> size;
  while ( size-- && is ) {
    is >> k;
    if(k==0)
      is >> m[ShowerInteraction::QCD];
    else if(k==1)
      is >> m[ShowerInteraction::QED];
    else
      assert(false);
  }
  return is;
}
}
 
#endif /* HERWIG_GeneralTwoBodyDecayer_H */