MatchboxPhasespace.h

// -*- C++ -*-
//
// MatchboxPhasespace.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_MatchboxPhasespace_H
#define HERWIG_MatchboxPhasespace_H
//
// This is the declaration of the MatchboxPhasespace class.
//
 
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig/MatrixElement/Matchbox/Utility/ProcessData.fh"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "Herwig/MatrixElement/Matchbox/Phasespace/PhasespaceCouplings.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
struct StreamingRnd {
 
  const double* numbers;
 
  size_t nRnd;
 
  StreamingRnd()
    : numbers(0), nRnd(0) {}
 
  explicit StreamingRnd(const double* newNumbers,
                        size_t n)
    : numbers(newNumbers), nRnd(n) {}
 
  inline double operator()() {
    assert(numbers && nRnd > 0);
    const double ret = numbers[0];
    ++numbers; --nRnd;
    return ret;
  }
 
};
 
class MatchboxPhasespace: 
    public HandlerBase, 
    public LastXCombInfo<StandardXComb>,
    public LastMatchboxXCombInfo {
 
public:
 
  MatchboxPhasespace();
 
public:
 
  virtual void setXComb(tStdXCombPtr xc) { 
    theLastXComb = xc;
    lastMatchboxXComb(xc);
  }
 
  Ptr<MatchboxFactory>::tcptr factory() const;
 
  Ptr<ProcessData>::tptr processData() const;
 
  virtual double generateKinematics(const double* r,
                                    vector<Lorentz5Momentum>& momenta);
 
  virtual double generateTwoToNKinematics(const double*,
                                          vector<Lorentz5Momentum>& momenta) = 0;
 
  virtual double generateTwoToOneKinematics(const double*,
                                            vector<Lorentz5Momentum>& momenta);
 
  virtual int nDim(const cPDVector&) const;
 
  virtual int nDimPhasespace(int nFinal) const = 0;
 
  virtual bool haveX1X2() const { return false; }
 
  virtual bool wantCMS() const { return true; }
 
  bool useMassGenerators() const { return theUseMassGenerators; }
 
  virtual Selector<MEBase::DiagramIndex> selectDiagrams(const MEBase::DiagramVector&) const;
 
  pair<double,Lorentz5Momentum> timeLikeWeight(const Tree2toNDiagram& diag,
                                               int branch, double flatCut) const;
 
  double spaceLikeWeight(const Tree2toNDiagram& diag,
                         const Lorentz5Momentum& incoming,
                         int branch, double flatCut) const;
 
  double diagramWeight(const Tree2toNDiagram& diag) const {
    assert( !diagramWeights().empty() );
    return diagramWeights().find(diag.id())->second;
  }
 
  void fillDiagramWeights(double flatCut = 0.0);
 
  void clearDiagramWeights() {
    diagramWeights().clear();
  }
 
  Ptr<MatchboxPhasespace>::ptr cloneMe() const {
    return dynamic_ptr_cast<Ptr<MatchboxPhasespace>::ptr>(clone());
  }
 
  virtual void cloneDependencies(const std::string& prefix = "");
 
public:
 
  virtual bool isInvertible() const { return false; }
 
  virtual double invertKinematics(const vector<Lorentz5Momentum>& momenta,
                                  double* r) const;
 
  virtual double invertTwoToNKinematics(const vector<Lorentz5Momentum>&,
                                        double*) const {
    return 0.;
  }
 
  virtual double invertTwoToOneKinematics(const vector<Lorentz5Momentum>&, double*) const;
 
public:
 
  void singularLimit(size_t i, size_t j) {
    if ( i > j )
      swap(i,j);
    singularLimits().insert(make_pair(i,j));
  }
 
  const pair<size_t,size_t>& lastSingularIndices() const {
    assert(lastSingularLimit() != singularLimits().end());
    return *lastSingularLimit();
  }
 
  bool matchConstraints(const vector<Lorentz5Momentum>& momenta);
 
protected:
 
  void setCoupling(long a, long b, long c,
                   double coupling, bool includeCrossings = true);
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
public:
 
  static void Init();
 
 
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
 
private:
 
  Energy singularCutoff;
 
  bool theUseMassGenerators;
 
  Ptr<PhasespaceCouplings>::ptr theCouplings;
 
  string doSetCoupling(string);
 
  string doSetPhysicalCoupling(string);
 
  int theLoopParticleIdMin;
 
  int theLoopParticleIdMax;
 
  MatchboxPhasespace & operator=(const MatchboxPhasespace &) = delete;
 
};
 
}
 
#endif /* HERWIG_MatchboxPhasespace_H */