FlatInvertiblePhasespace.h

// -*- C++ -*-
//
// FlatInvertiblePhasespace.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_FlatInvertiblePhasespace_H
#define Herwig_FlatInvertiblePhasespace_H
//
// This is the declaration of the FlatInvertiblePhasespace class.
//
 
#include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class FlatInvertiblePhasespace: public MatchboxPhasespace {
 
public:
 
  virtual double generateTwoToNKinematics(const double*,
                                          vector<Lorentz5Momentum>& momenta);
 
  virtual int nDimPhasespace(int nFinal) const {
    if ( nFinal == 1 )
      return 1;
    return 3*nFinal - 4;
  }
 
public:
 
  virtual bool isInvertible() const { return true; }
 
  virtual double invertTwoToNKinematics(const vector<Lorentz5Momentum>& momenta,
                                        double* r) const {
    return invertKinematics(momenta,(momenta[0]+momenta[1]).m(),r);
  }
 
private:
 
  double bisect(double v, double n, 
                double target = -16., double maxLevel = 80.) const;
 
  double rho(Energy M, Energy N, Energy m) const {
    return sqrt((sqr(M)-sqr(N+m))*(sqr(M)-sqr(N-m)))/(8.*sqr(M));
  }
 
  double generateIntermediates(vector<Energy>& K,
                               const double* r) const;
 
  double invertIntermediates(const vector<Energy>& K,
                             double* r) const;
 
  double generateIntermediates(vector<Energy>& M,
                               const vector<Energy>& m,
                               const double* r) const;
 
  double invertIntermediates(const vector<Energy>& M,
                             const vector<Energy>& m,
                             double* r) const;
 
  double generateKinematics(vector<Lorentz5Momentum>& P,
                            Energy Ecm,
                            const double* r) const;
 
  double invertKinematics(const vector<Lorentz5Momentum>& P,
                          Energy Ecm,
                          double* r) const;
 
  long double flatWeights(int n) const;
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
 
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
 
private:
 
  FlatInvertiblePhasespace & operator=(const FlatInvertiblePhasespace &) = delete;
 
};
 
}
 
#endif /* Herwig_FlatInvertiblePhasespace_H */