Page MenuHomeHEPForge
Files F8725592

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc
@@ -0,0 +1,419 @@
+// -*- C++ -*-
+//
+// FlatInvertiblePhasespace.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
+// Copyright (C) 2002-2012 The Herwig Collaboration
+//
+// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the FlatInvertiblePhasespace class.
+//
+
+#include "FlatInvertiblePhasespace.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "Herwig++/Utilities/GSLBisection.h"
+#include "ThePEG/Cuts/Cuts.h"
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+
+FlatInvertiblePhasespace::FlatInvertiblePhasespace() {}
+
+FlatInvertiblePhasespace::~FlatInvertiblePhasespace() {}
+
+IBPtr FlatInvertiblePhasespace::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr FlatInvertiblePhasespace::fullclone() const {
+ return new_ptr(*this);
+}
+
+void FlatInvertiblePhasespace::prepare(tStdXCombPtr xc, bool) {
+ theLastXComb = xc;
+}
+
+double FlatInvertiblePhasespace::bisect(double v, double n,
+ double target, double maxLevel) const {
+
+ if ( v != 0.0 && v != 1.0 ) {
+
+ double level = 0;
+ double left = 0;
+ double right = 1;
+
+ double checkV = -1.;
+ double u = -1;
+
+ while ( level < maxLevel ) {
+
+ u = (left+right)*pow(0.5,level+1.);
+ checkV =
+ pow(u,n+1.)*(n+2.-(n+1.)*u);
+
+ if ( log10(abs(1.-checkV/v)) <= target )
+ break;
+
+ left *= 2.;
+ right *= 2.;
+
+ if ( v <= checkV ) {
+ right -= 1.;
+ ++level;
+ }
+
+ if ( v > checkV ) {
+ left += 1.;
+ ++level;
+ }
+
+ }
+
+ return u;
+
+ }
+
+ return v;
+
+}
+
+static double flatWeights[7] = {
+
+ -1.,-1.,
+ 0.039788735772973833942,
+ 0.00012598255637968550463,
+ 1.3296564302788840628E-7,
+ 7.0167897579949011130E-11,
+ 2.2217170114046130768E-14
+
+};
+
+double FlatInvertiblePhasespace::generateIntermediates(vector<Energy>& K,
+ const double* r) const {
+
+ size_t n = K.size() + 1;
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ double u = bisect(r[i-2],n-1-i);
+ K[i-1] = sqrt(u*sqr(K[i-2]));
+ }
+
+ return flatWeights[n];
+
+}
+
+double FlatInvertiblePhasespace::invertIntermediates(const vector<Energy>& K,
+ double* r) const {
+
+ size_t n = K.size() + 1;
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ double u = sqr(K[i-1]/K[i-2]);
+ r[i-2] = (n+1-i)*pow(u,(double)(n-i)) - (n-i)*pow(u,(double)(n+1-i));
+ }
+
+ return flatWeights[n];
+
+}
+
+double FlatInvertiblePhasespace::generateIntermediates(vector<Energy>& M,
+ const vector<Energy>& m,
+ const double* r) const {
+
+ size_t n = M.size() + 1;
+
+ vector<Energy> K = M;
+ for ( size_t i = 1; i <= n; ++i )
+ K[0] -= m[i-1];
+
+ double w0 = generateIntermediates(K,r);
+
+ M = K;
+ for ( size_t i = 1; i <= n-1; ++i ) {
+ for ( size_t k = i; k <= n; ++k )
+ M[i-1] += m[k-1];
+ }
+
+ double weight = 8.*w0*rho(M[n-2],m[n-1],m[n-2]);
+
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ weight *=
+ (rho(M[i-2],M[i-1],m[i-2])/rho(K[i-2],K[i-1],ZERO)) * (M[i-1]/K[i-1]);
+ }
+
+ weight *= pow(K[0]/M[0],2.*n-4.);
+
+ return weight;
+
+}
+
+double FlatInvertiblePhasespace::invertIntermediates(const vector<Energy>& M,
+ const vector<Energy>& m,
+ double* r) const {
+
+ size_t n = M.size() + 1;
+
+ vector<Energy> K = M;
+ for ( size_t i = 1; i <= n-1; ++i ) {
+ for ( size_t k = i; k <= n; ++k )
+ K[i-1] -= m[k-1];
+ }
+
+ double w0 = invertIntermediates(K,r);
+
+ double weight = 8.*w0*rho(M[n-2],m[n-1],m[n-2]);
+
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ weight *=
+ (rho(M[i-2],M[i-1],m[i-2])/rho(K[i-2],K[i-1],ZERO)) * (M[i-1]/K[i-1]);
+ }
+
+ weight *= pow(K[0]/M[0],2.*n-4.);
+
+ return weight;
+
+}
+
+
+double FlatInvertiblePhasespace::generateKinematics(vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ const double* r) const {
+
+ vector<Energy> m;
+ for ( vector<Lorentz5Momentum>::const_iterator p =
+ P.begin(); p != P.end(); ++p )
+ m.push_back(p->mass());
+
+ size_t n = P.size();
+ vector<Energy> M(n-1);
+ M[0] = Ecm;
+
+ double weight = generateIntermediates(M,m,r);
+
+ M.push_back(m.back());
+
+ Lorentz5Momentum Q(M[0]);
+ Lorentz5Momentum nextQ;
+
+ for ( size_t i = 2; i <= n; ++i ) {
+
+ Energy q = 4.*M[i-2]*rho(M[i-2],M[i-1],m[i-2]);
+
+ double c = 2.*r[n-6+2*i]-1.;
+ double s = sqrt(1.-sqr(c));
+ double phi = 2.*Constants::pi*r[n-5+2*i];
+ double cphi = cos(phi);
+ double sphi = sqrt(1.-sqr(cphi));
+ if ( phi > Constants::pi )
+ sphi = -sphi;
+
+ P[i-2].setX(q*cphi*s);
+ P[i-2].setY(q*sphi*s);
+ P[i-2].setZ(q*c);
+ P[i-2].rescaleEnergy();
+ P[i-2].boost(Q.boostVector());
+ P[i-2].rescaleEnergy();
+
+ nextQ = Q - P[i-2];
+ nextQ.setMass(M[i-1]);
+ nextQ.rescaleEnergy();
+
+ Q = nextQ;
+
+ }
+
+ P.back() = Q;
+
+ return weight;
+
+}
+
+double FlatInvertiblePhasespace::invertKinematics(const vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ double* r) const {
+
+ vector<Energy> m;
+ for ( vector<Lorentz5Momentum>::const_iterator p =
+ P.begin(); p != P.end(); ++p )
+ m.push_back(p->mass());
+
+ size_t n = P.size();
+ vector<Energy> M(n-1);
+ M[0] = Ecm;
+
+ vector<Lorentz5Momentum> Q(n-1);
+ Q[0] = Lorentz5Momentum(M[0]);
+
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ for ( size_t k = i; k <= n; ++k )
+ Q[i-1] += P[k-1];
+ M[i-1] = Q[i-1].m();
+ }
+
+ double weight = invertIntermediates(M,m,r);
+
+ for ( size_t i = 2; i <= n; ++i ) {
+ Lorentz5Momentum p = P[i-2];
+ p.boost(-Q[i-2].boostVector());
+ r[n-6+2*i] = (p.cosTheta()+1.)/2.;
+ double phi = p.phi();
+ if ( phi < 0. )
+ phi = 2.*Constants::pi + phi;
+ r[n-5+2*i] = phi/(2.*Constants::pi);
+ }
+
+ return weight;
+
+}
+
+
+double FlatInvertiblePhasespace::generateKinematics(const double* r,
+ vector<Lorentz5Momentum>& momenta) {
+
+ cPDVector::const_iterator pd = mePartonData().begin();
+ vector<Lorentz5Momentum>::iterator p = momenta.begin();
+ for ( ; pd != mePartonData().end(); ++pd, ++p )
+ p->setMass((**pd).mass());
+
+ Energy ECMMin = ZERO;
+ pd = mePartonData().begin() + 2;
+ for ( ; pd != mePartonData().end(); ++pd )
+ ECMMin += (**pd).mass();
+
+ double weight = 1.;
+
+ double tauMin = sqr(ECMMin)/lastXCombPtr()->lastS();
+ double tau = 1.;
+
+ if ( momenta.size() > 3 ) {
+ tau = tauMin + r[1]*(1.-tauMin);
+ } else {
+ tau = tauMin;
+ }
+
+ double ltau = log(tau)/2.;
+ weight *= -2.*ltau;
+
+ if ( momenta.size() > 3 ) {
+ weight *= 1.-tauMin;
+ } else {
+ weight *= 2.*Constants::pi;
+ }
+
+ double y = ltau - 2.*r[0]*ltau;
+ double x1 = sqrt(tau)*exp(y);
+ double x2 = sqrt(tau)*exp(-y);
+
+ ThreeVector<Energy> p1 =
+ x1*(lastXCombPtr()->lastParticles().first->momentum().vect());
+
+ ThreeVector<Energy> p2 =
+ x2*(lastXCombPtr()->lastParticles().second->momentum().vect());
+
+ momenta[0] = Lorentz5Momentum(momenta[0].mass(),p1);
+ momenta[1] = Lorentz5Momentum(momenta[1].mass(),p2);
+
+ lastXCombPtr()->lastX1X2(make_pair(x1,x2));
+ lastXCombPtr()->lastSHat((momenta[0]+momenta[1]).m2());
+
+ if ( momenta.size() == 3 ) {
+ ThreeVector<Energy> q = p1 + p2;
+ momenta[2] = Lorentz5Momentum(momenta[2].mass(),q);
+ return weight;
+ }
+
+ vector<Lorentz5Momentum> cmMomenta;
+ copy(momenta.begin()+2,momenta.end(),back_inserter(cmMomenta));
+
+ weight *= generateKinematics(cmMomenta,sqrt(lastXCombPtr()->lastSHat()),r+2);
+
+ Boost toLab = (momenta[0]+momenta[1]).boostVector();
+
+ vector<Lorentz5Momentum>::iterator pcm = cmMomenta.begin();
+ vector<Lorentz5Momentum>::iterator plab = momenta.begin() + 2;
+
+ for ( ; pcm != cmMomenta.end(); ++pcm, ++plab )
+ *plab = pcm->boost(toLab);
+
+ return weight;
+
+}
+
+double FlatInvertiblePhasespace::invertKinematics(const vector<Lorentz5Momentum>& momenta,
+ double* r) const {
+
+ Energy ECMMin = ZERO;
+ cPDVector::const_iterator pd = mePartonData().begin() + 2;
+ for ( ; pd != mePartonData().end(); ++pd )
+ ECMMin += (**pd).mass();
+
+ double weight = 1.;
+
+ double tauMin = sqr(ECMMin)/lastXCombPtr()->lastS();
+ double tau = lastXCombPtr()->lastTau();
+
+ if ( momenta.size() > 3 ) {
+ r[1] = (tau - tauMin)/(1.-tauMin);
+ }
+
+ double ltau = log(tau)/2.;
+ weight *= -2.*ltau;
+
+ if ( momenta.size() > 3 ) {
+ weight *= 1.-tauMin;
+ } else {
+ weight *= 2.*Constants::pi;
+ }
+
+ double y = lastXCombPtr()->lastY();
+ r[0] = (ltau - y)/(2.*ltau);
+
+ if ( momenta.size() == 3 )
+ return weight;
+
+ vector<Lorentz5Momentum> cmMomenta;
+ copy(momenta.begin()+2,momenta.end(),back_inserter(cmMomenta));
+
+ Boost toCM = -(momenta[0]+momenta[1]).boostVector();
+
+ for ( vector<Lorentz5Momentum>::iterator p = cmMomenta.begin();
+ p != cmMomenta.end(); ++p )
+ p->boost(toCM);
+
+ weight *= invertKinematics(cmMomenta,sqrt(lastXCombPtr()->lastSHat()),r+2);
+
+ return weight;
+
+}
+
+// If needed, insert default implementations of virtual function defined
+// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
+
+
+void FlatInvertiblePhasespace::persistentOutput(PersistentOStream &) const {}
+
+void FlatInvertiblePhasespace::persistentInput(PersistentIStream &, int) {}
+
+
+// *** Attention *** The following static variable is needed for the type
+// description system in ThePEG. Please check that the template arguments
+// are correct (the class and its base class), and that the constructor
+// arguments are correct (the class name and the name of the dynamically
+// loadable library where the class implementation can be found).
+DescribeClass<FlatInvertiblePhasespace,MatchboxPhasespace>
+ describeHerwigFlatInvertiblePhasespace("Herwig::FlatInvertiblePhasespace", "HwMatchbox.so");
+
+void FlatInvertiblePhasespace::Init() {
+
+ static ClassDocumentation<FlatInvertiblePhasespace> documentation
+ ("FlatInvertiblePhasespace implements flat, invertible phase space generation.");
+
+}
+
diff --git a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h
@@ -0,0 +1,209 @@
+// -*- C++ -*-
+//
+// FlatInvertiblePhasespace.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
+// Copyright (C) 2002-2012 The Herwig Collaboration
+//
+// Herwig++ is licenced under version 2 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;
+
+/**
+ * \ingroup Matchbox
+ * \author Simon Platzer
+ *
+ * \brief FlatInvertiblePhasespace implements flat, invertible phase space generation.
+ *
+ */
+class FlatInvertiblePhasespace: public MatchboxPhasespace {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ FlatInvertiblePhasespace();
+
+ /**
+ * The destructor.
+ */
+ virtual ~FlatInvertiblePhasespace();
+ //@}
+
+public:
+
+ /**
+ * Prepare a phase space generator for the given xcomb object.
+ */
+ virtual void prepare(tStdXCombPtr, bool verbose = false);
+
+ /**
+ * Generate a phase space point and return its weight.
+ */
+ virtual double generateKinematics(const double*,
+ vector<Lorentz5Momentum>& momenta);
+
+ /**
+ * Return the number of random numbers required to produce a given
+ * multiplicity final state.
+ */
+ virtual int nDim(int nFinal) const {
+ if ( nFinal == 1 )
+ return 1;
+ return 3*nFinal - 2;
+ }
+
+ /**
+ * Return true, if this phasespace generator will generate incoming
+ * partons itself.
+ */
+ virtual bool haveX1X2() const { return true; }
+
+ /**
+ * Return true, if this phase space generator expects
+ * the incoming partons in their center-of-mass system
+ */
+ virtual bool wantCMS() const { return false; }
+
+public:
+
+ /**
+ * Return true, if this phase space generator is invertible
+ */
+ virtual bool isInvertible() const { return true; }
+
+ /**
+ * Invert the given phase space point to the random numbers which
+ * would have generated it.
+ */
+ virtual double invertKinematics(const vector<Lorentz5Momentum>&,
+ double*) const;
+
+private:
+
+ /**
+ * Solve v = (n+2) * u^(n+1) - (n+1) * u^(n+2) for u
+ */
+ double bisect(double v, double n,
+ double target = -16., double maxLevel = 80.) const;
+
+ /**
+ * Return rho
+ */
+ double rho(Energy M, Energy N, Energy m) const {
+ return sqrt((sqr(M)-sqr(N+m))*(sqr(M)-sqr(N-m)))/(8.*sqr(M));
+ }
+
+ /**
+ * Generate intermediate masses for a massless final state
+ */
+ double generateIntermediates(vector<Energy>& K,
+ const double* r) const;
+
+ /**
+ * Invert intermediate masses for a massless final state
+ */
+ double invertIntermediates(const vector<Energy>& K,
+ double* r) const;
+
+ /**
+ * Generate intermediate masses for a massive final state
+ */
+ double generateIntermediates(vector<Energy>& M,
+ const vector<Energy>& m,
+ const double* r) const;
+
+ /**
+ * Invert intermediate masses for a massive final state
+ */
+ double invertIntermediates(const vector<Energy>& M,
+ const vector<Energy>& m,
+ double* r) const;
+
+ /**
+ * Generate momenta in the CMS
+ */
+ double generateKinematics(vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ const double* r) const;
+
+ /**
+ * Invert momenta in the CMS
+ */
+ double invertKinematics(const vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ double* r) const;
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr fullclone() const;
+ //@}
+
+
+// If needed, insert declarations of virtual function defined in the
+// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
+
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ FlatInvertiblePhasespace & operator=(const FlatInvertiblePhasespace &);
+
+};
+
+}
+
+#endif /* Herwig_FlatInvertiblePhasespace_H */
diff --git a/MatrixElement/Matchbox/Phasespace/Makefile.am b/MatrixElement/Matchbox/Phasespace/Makefile.am
--- a/MatrixElement/Matchbox/Phasespace/Makefile.am
+++ b/MatrixElement/Matchbox/Phasespace/Makefile.am
@@ -1,50 +1,52 @@
if WANT_DIPOLE
noinst_LTLIBRARIES = libHwMatchboxPhasespace.la
endif
libHwMatchboxPhasespace_la_SOURCES = \
FFLightInvertedTildeKinematics.cc \
FFLightInvertedTildeKinematics.h \
FFLightTildeKinematics.cc \
FFLightTildeKinematics.h \
FFMassiveInvertedTildeKinematics.cc \
FFMassiveInvertedTildeKinematics.h \
FFMassiveTildeKinematics.cc \
FFMassiveTildeKinematics.h \
FILightInvertedTildeKinematics.cc \
FILightInvertedTildeKinematics.h \
FILightTildeKinematics.cc \
FILightTildeKinematics.h \
FIMassiveInvertedTildeKinematics.cc \
FIMassiveInvertedTildeKinematics.h \
FIMassiveTildeKinematics.cc \
FIMassiveTildeKinematics.h \
IFLightInvertedTildeKinematics.cc \
IFLightInvertedTildeKinematics.h \
IFLightTildeKinematics.cc \
IFLightTildeKinematics.h \
IFMassiveInvertedTildeKinematics.cc \
IFMassiveInvertedTildeKinematics.h \
IFMassiveTildeKinematics.cc \
IFMassiveTildeKinematics.h \
IILightInvertedTildeKinematics.cc \
IILightInvertedTildeKinematics.h \
IILightTildeKinematics.cc \
IILightTildeKinematics.h \
InvertedTildeKinematics.cc \
InvertedTildeKinematics.h \
MatchboxPhasespace.cc \
MatchboxPhasespace.h \
MatchboxRambo.cc \
MatchboxRambo.h \
PhasespaceHelpers.cc \
PhasespaceHelpers.h \
RandomHelpers.h \
TildeKinematics.cc \
TildeKinematics.h \
TreePhasespace.cc \
TreePhasespace.h \
TreePhasespaceChannels.cc \
TreePhasespaceChannels.h \
MatchboxReference.h \
-MatchboxReference.cc
+MatchboxReference.cc \
+FlatInvertiblePhasespace.h \
+FlatInvertiblePhasespace.cc
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
--- a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
@@ -1,260 +1,276 @@
// -*- C++ -*-
//
// MatchboxPhasespace.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Wrap around a vector of random numbers to behave as a stream
* of those.
*/
struct StreamingRnd {
/**
* The random numbers
*/
const double* numbers;
/**
* The number of random numbers available.
*/
size_t nRnd;
/**
* Default constructor.
*/
StreamingRnd()
: numbers(0), nRnd(0) {}
/**
* Construct from random numbers.
*/
explicit StreamingRnd(const double* newNumbers,
size_t n)
: numbers(newNumbers), nRnd(n) {}
/**
* Return next random number
*/
inline double operator()() {
assert(numbers && nRnd > 0);
const double ret = numbers[0];
++numbers; --nRnd;
return ret;
}
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxPhasespace defines an abstract interface to a phase
* space generator.
*
*/
class MatchboxPhasespace:
public HandlerBase, public LastXCombInfo<StandardXComb> {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxPhasespace();
/**
* The destructor.
*/
virtual ~MatchboxPhasespace();
//@}
public:
/**
* Prepare a phase space generator for the given xcomb object.
*/
virtual void prepare(tStdXCombPtr, bool verbose = false) = 0;
/**
* Generate a phase space point and return its weight.
*/
virtual double generateKinematics(const double*,
vector<Lorentz5Momentum>& momenta) = 0;
/**
* Return the number of random numbers required to produce a given
* multiplicity final state.
*/
virtual int nDim(int nFinal) const = 0;
/**
* Return true, if this phasespace generator will generate incoming
* partons itself.
*/
virtual bool haveX1X2() const { return false; }
/**
* Return true, if this phase space generator expects
* the incoming partons in their center-of-mass system
*/
virtual bool wantCMS() const { return true; }
/**
* Fill a diagram selector for the last phase space point.
*/
virtual Selector<MEBase::DiagramIndex> selectDiagrams(const MEBase::DiagramVector&) const;
/**
* Return the momentum and weight appropriate to the given timelike
* branch of the diagram.
*/
pair<double,Lorentz5Momentum> timeLikeWeight(const Tree2toNDiagram& diag,
int branch, double flatCut) const;
/**
* Return the weight appropriate to the given spacelike branch of
* the diagram.
*/
double spaceLikeWeight(const Tree2toNDiagram& diag,
const Lorentz5Momentum& incoming,
int branch, double flatCut) const;
/**
* Return the weight appropriate to the given diagram.
*/
double diagramWeight(const Tree2toNDiagram& diag) const {
assert( !diagramWeights.empty() );
return diagramWeights.find(diag.id())->second;
}
/**
* Fill the diagram weights.
*/
void fillDiagramWeights(double flatCut = 0.0);
/**
* Clone this phase space generator.
*/
Ptr<MatchboxPhasespace>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxPhasespace>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
virtual void cloneDependencies(const std::string& prefix = "");
public:
/**
+ * Return true, if this phase space generator is invertible
+ */
+ virtual bool isInvertible() const { return false; }
+
+ /**
+ * Invert the given phase space point to the random numbers which
+ * would have generated it.
+ */
+ virtual double invertKinematics(const vector<Lorentz5Momentum>&,
+ double*) const {
+ return 0.;
+ }
+
+public:
+
+ /**
* Limit phasespace generation to a given collinear or soft limit.
*/
void singularLimit(size_t i, size_t j) {
if ( i > j )
swap(i,j);
singularLimits.insert(make_pair(i,j));
}
/**
* Return the last matched singular limit.
*/
const pair<size_t,size_t>& lastSingularLimit() const {
assert(theLastSingularLimit != singularLimits.end());
return *theLastSingularLimit;
}
/**
* Return true, if constraints on phasespace generation have been met.
*/
bool matchConstraints(const vector<Lorentz5Momentum>& momenta);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
public:
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The diagram weights indexed by diagram id.
*/
map<int,double> diagramWeights;
/**
* If not empty, the entries here serve to limit phasespace
* generation to the corresponding collinear limits, or soft limits
* if both pair entries are equal.
*/
set<pair<size_t,size_t> > singularLimits;
/**
* The last matched singular limit.
*/
set<pair<size_t,size_t> >::const_iterator theLastSingularLimit;
/**
* A cutoff below which a region is considered singular.
*/
Energy singularCutoff;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxPhasespace & operator=(const MatchboxPhasespace &);
};
}
#endif /* HERWIG_MatchboxPhasespace_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 1:59 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243581
Default Alt Text
(24 KB)

Event Timeline