Page MenuHomeHEPForge
Files F9501285

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc
@@ -1,310 +1,321 @@
// -*- C++ -*-
//
// FFMassiveInvertedTildeKinematics.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 FFMassiveInvertedTildeKinematics class.
//
#include "FFMassiveInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
using namespace Herwig;
-FFMassiveInvertedTildeKinematics::FFMassiveInvertedTildeKinematics() {}
+FFMassiveInvertedTildeKinematics::FFMassiveInvertedTildeKinematics()
+ : theFullMapping(false) {}
FFMassiveInvertedTildeKinematics::~FFMassiveInvertedTildeKinematics() {}
IBPtr FFMassiveInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FFMassiveInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FFMassiveInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
double mapping = 1.0;
vector<double> values;
pair<Energy,double> ptz = generatePtZ(mapping,r, &values);
if ( mapping == 0.0 ) {
jacobian(0.0);
return false;
}
// pt and zPrime = qi.nk / (qi+qj).nk are the generated variables
Energy pt = ptz.first;
Energy2 pt2 = sqr(pt);
double zPrime = ptz.second;
// Define the scale
Energy scale = (emitter+spectator).m();
Energy2 Qijk = sqr(scale);
// Most of the required values have been calculated in ptzAllowed
double y = values[0];
double z = values[1];
double xk = values[3];
double xij = values[4];
double suijk = values[5];
double suijk2 = sqr(values[5]);
// masses
double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
double mu2 = sqr( realEmissionData()->hardProcessMass() / scale );
double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
double Mui2 = sqr( bornEmitterData()->hardProcessMass() / scale );
double Muj2 = sqr( bornSpectatorData()->hardProcessMass() / scale );
// Construct reference momenta nk, nij, nt
Lorentz5Momentum nij = ( suijk2 / (suijk2-Mui2*Muj2) ) * (emitter - (Mui2/suijk)*spectator);
Lorentz5Momentum nk = ( suijk2 / (suijk2-Mui2*Muj2) ) * (spectator - (Muj2/suijk)*emitter);
// Following notation in notes, qt = sqrt(wt)*nt
double phi = 2.*Constants::pi*r[2];
Lorentz5Momentum qt = getKt(nij,nk,pt,phi);
// Construct qij, qk, qi and qj
Lorentz5Momentum qij = xij*nij + (Mui2/(xij*suijk))*nk;
Lorentz5Momentum spe = xk*nk + (Muj2/(xk*suijk))*nij;
Lorentz5Momentum em = zPrime*qij + ((pt2/Qijk + mui2 - zPrime*zPrime*Mui2)/(xij*suijk*zPrime))*nk + qt;
Lorentz5Momentum emm = (1.-zPrime)*qij + ((pt2/Qijk + mu2 - sqr(1.-zPrime)*Mui2)/(xij*suijk*(1.-zPrime)))*nk - qt;
em.setMass(realEmitterData()->hardProcessMass());
em.rescaleEnergy();
emm.setMass(realEmissionData()->hardProcessMass());
emm.rescaleEnergy();
spe.setMass(realSpectatorData()->hardProcessMass());
spe.rescaleEnergy();
// book
realEmitterMomentum() = em;
realEmissionMomentum() = emm;
realSpectatorMomentum() = spe;
// Store the jacobian
- mapping /= z*(1.-z);
- jacobian( mapping*(1.-y)*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)) *
- sqr(1.-mui2-mu2-muj2) / rootOfKallen(1.,Mui2,Muj2) );
+ // The jacobian here corresponds to dpt2 / sqr(lastscale) NOT dpt2 / pt2.
+ // This jacobian is the one-particle phase space
+ if ( true ) {
+ double bar = 1.-mui2-mu2-muj2;
+ mapping *= (1.-y)/(z*(1.-z));
+ jacobian( mapping*(sqr(lastScale())/sHat())*bar/rootOfKallen(1.,Mui2,Muj2) / (16.*sqr(Constants::pi)) );
+ }
- // Store the parameters
+ // Todo: Switch for the full jacobian including the z->zPrime jacobian
+ else if ( theFullMapping ) {
+ assert ( false );
+ }
+
+// Store the parameters
subtractionParameters().resize(3);
subtractionParameters()[0] = y;
subtractionParameters()[1] = z;
subtractionParameters()[2] = zPrime;
return true;
}
Energy FFMassiveInvertedTildeKinematics::lastPt() const {
Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m();
// masses
double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
double mu2 = sqr( realEmissionData()->hardProcessMass() / scale );
double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
double y = subtractionParameters()[0];
double zPrime = subtractionParameters()[2];
Energy ret = scale * sqrt( y * (1.-mui2-mu2-muj2) * zPrime*(1.-zPrime) - sqr(1.-zPrime)*mui2 - sqr(zPrime)*mu2 );
return ret;
}
double FFMassiveInvertedTildeKinematics::lastZ() const {
return subtractionParameters()[2];
}
Energy FFMassiveInvertedTildeKinematics::ptMax() const {
Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m();
// masses
double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
double mu2 = sqr( realEmissionData()->hardProcessMass() / scale );
double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
Energy ptmax = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
( 2.-2.*sqrt(muj2) ) * scale;
return ptmax > 0.*GeV ? ptmax : 0.*GeV;
}
// NOTE: bounds calculated at this step may be too loose
// These apply to zPrime, which is stored in lastZ()
pair<double,double> FFMassiveInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
if(pt>hardPt) return make_pair(0.5,0.5);
Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m();
// masses
double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
double mu2 = sqr( realEmissionData()->hardProcessMass() / scale );
double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
double zp = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) +
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
sqrt( 1.-sqr(pt/hardPt) ) ) /
( 2.*sqr(1.-sqrt(muj2)) );
double zm = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) -
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
sqrt( 1.-sqr(pt/hardPt) ) ) /
( 2.*sqr(1.-sqrt(muj2)) );
return make_pair(zm,zp);
}
bool FFMassiveInvertedTildeKinematics::ptzAllowed(pair<Energy,double> ptz, vector<double>* values) const {
Energy pt = ptz.first;
Energy2 pt2 = sqr(pt);
double zPrime = ptz.second;
// masses
double mui2 = sqr( realEmitterData()->hardProcessMass() / lastScale() );
double mu2 = sqr( realEmissionData()->hardProcessMass() / lastScale() );
double muj2 = sqr( realSpectatorData()->hardProcessMass() / lastScale() );
double Mui2 = sqr( bornEmitterData()->hardProcessMass() / lastScale() );
double Muj2 = sqr( bornSpectatorData()->hardProcessMass() / lastScale() );
// Calculate the scales that we need
Energy2 Qijk = sqr(lastScale());
double suijk = 0.5*( 1. - Mui2 - Muj2 + sqrt( sqr(1.-Mui2-Muj2) - 4.*Mui2*Muj2 ) );
double bar = 1.-mui2-mu2-muj2;
double y = ( pt2/Qijk + sqr(1.-zPrime)*mui2 + zPrime*zPrime*mu2 ) /
(zPrime*(1.-zPrime)*bar);
// Calculate A:=xij*w
double A = (1./(suijk*zPrime*(1.-zPrime))) * ( pt2/Qijk + zPrime*mu2 + (1.-zPrime)*mui2 - zPrime*(1.-zPrime)*Mui2 );
// Calculate the scaling factors, xk and xij
double lambdaK = 1. + (Muj2/suijk);
double lambdaIJ = 1. + (Mui2/suijk);
double xk = (1./(2.*lambdaK)) * ( (lambdaK + (Muj2/suijk)*lambdaIJ - A) + sqrt( sqr(lambdaK + (Muj2/suijk)*lambdaIJ - A) - 4.*lambdaK*lambdaIJ*Muj2/suijk) );
double xij = 1. - ( (Muj2/suijk) * (1.-xk) / xk );
// Calculate z
double z = ( (zPrime*xij*xk*suijk/2.) + (Muj2/ ( 2.*xk*xij*suijk*zPrime))*(pt2/Qijk + mui2) ) /
( (xij*xk*suijk/2.) + (Muj2/(2.*xk*xij))*(Mui2/suijk + A) );
// check (y,z) phasespace boundary
// TODO: is y boundary necessary?
double ym = 2.*sqrt(mui2)*sqrt(mu2)/bar;
double yp = 1. - 2.*sqrt(muj2)*(1.-sqrt(muj2))/bar;
// These limits apply to z, not zPrime
double zm = ( (2.*mui2+bar*y)*(1.-y) - sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
double zp = ( (2.*mui2+bar*y)*(1.-y) + sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
if ( y<ym || y>yp || z<zm || z>zp ) return false;
values->push_back(y);
values->push_back(z);
values->push_back(A);
values->push_back(xk);
values->push_back(xij);
values->push_back(suijk);
return true;
}
// This is used to generate pt and zPrime
pair<Energy,double> FFMassiveInvertedTildeKinematics::generatePtZ(double& jac, const double * r, vector<double> * values) const {
double kappaMin =
ptCut() != ZERO ?
sqr(ptCut()/ptMax()) :
sqr(0.1*GeV/GeV);
double kappa;
using namespace RandomHelpers;
if ( ptCut() > ZERO ) {
pair<double,double> kw =
generate(inverse(0.,kappaMin,1.),r[0]);
kappa = kw.first;
jac *= kw.second;
} else {
pair<double,double> kw =
generate((piecewise(),
flat(1e-4,kappaMin),
match(inverse(0.,kappaMin,1.))),r[0]);
kappa = kw.first;
jac *= kw.second;
}
Energy pt = sqrt(kappa)*ptMax();
pair<double,double> zLims = zBounds(pt);
pair<double,double> zw =
generate(inverse(0.,zLims.first,zLims.second)+
inverse(1.,zLims.first,zLims.second),r[1]);
double z = zw.first;
jac *= zw.second;
jac *= sqr(ptMax()/lastScale());
if( !ptzAllowed(make_pair(pt,z), values )) jac = 0.;
return make_pair(pt,z);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFMassiveInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FFMassiveInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FFMassiveInvertedTildeKinematics::Init() {
static ClassDocumentation<FFMassiveInvertedTildeKinematics> documentation
("FFMassiveInvertedTildeKinematics inverts the final-final tilde "
"kinematics.");
}
// *** 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<FFMassiveInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigFFMassiveInvertedTildeKinematics("Herwig::FFMassiveInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.h
@@ -1,178 +1,183 @@
// -*- C++ -*-
//
// FFMassiveInvertedTildeKinematics.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_FFMassiveInvertedTildeKinematics_H
#define HERWIG_FFMassiveInvertedTildeKinematics_H
//
// This is the declaration of the FFMassiveInvertedTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer, Stephen Webster
*
* \brief FFMassiveInvertedTildeKinematics inverts the final-final tilde
* kinematics.
*
*/
class FFMassiveInvertedTildeKinematics: public Herwig::InvertedTildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FFMassiveInvertedTildeKinematics();
/**
* The destructor.
*/
virtual ~FFMassiveInvertedTildeKinematics();
//@}
public:
/**
* Perform the mapping of the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the real
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap(const double *);
/**
* Return the pt associated to the last generated splitting.
*/
virtual Energy lastPt() const;
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() const;
/**
* Return the upper bound on pt
*/
virtual Energy ptMax() const;
/**
* Given a pt, return the boundaries on z
* Note that allowing parton masses these bounds may be too loose
*/
virtual pair<double,double> zBounds(Energy pt, Energy hardPt = ZERO) const;
/**
* For generated pt and z, check if this point is
* kinematically allowed
*/
/*virtual*/ bool ptzAllowed(pair<Energy,double> ptz, vector<double>* values ) const;
/**
* Generate pt and z
*/
virtual pair<Energy,double> generatePtZ(double& jac, const double * r, vector<double>* values) const;
public:
/**
* Triangular / Kallen function
*/
template <class T>
inline T rootOfKallen (T a, T b, T c) const {
return sqrt( a*a + b*b + c*c - 2.*( a*b+a*c+b*c ) ); }
// TODO: remove in both
/**
* stolen from FFMassiveKinematics.h
* Perform a rotation on both momenta such that the first one will
* point along the (positive) z axis. Rotate back to the original
* reference frame by applying rotateUz(returnedVector) to each momentum.
*/
ThreeVector<double> rotateToZ (Lorentz5Momentum& pTarget, Lorentz5Momentum& p1) {
ThreeVector<double> oldAxis = pTarget.vect().unit();
double ct = oldAxis.z(); double st = sqrt( 1.-sqr(ct) ); // cos,sin(theta)
double cp = oldAxis.x()/st; double sp = oldAxis.y()/st; // cos,sin(phi)
pTarget.setZ( pTarget.vect().mag() ); pTarget.setX( 0.*GeV ); pTarget.setY( 0.*GeV );
Lorentz5Momentum p1old = p1;
p1.setX( sp*p1old.x() - cp*p1old.y() );
p1.setY( ct*cp*p1old.x() + ct*sp*p1old.y() - st*p1old.z() );
p1.setZ( st*cp*p1old.x() + st*sp*p1old.y() + ct*p1old.z() );
return oldAxis;
}
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.
*/
FFMassiveInvertedTildeKinematics & operator=(const FFMassiveInvertedTildeKinematics &);
+ /**
+ * Option to use the full jacobian, including the z->zprime jacobian.
+ **/
+ bool theFullMapping;
+
};
}
#endif /* HERWIG_FFMassiveInvertedTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
@@ -1,134 +1,136 @@
// -*- C++ -*-
//
// FILightInvertedTildeKinematics.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 FILightInvertedTildeKinematics class.
//
#include "FILightInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FILightInvertedTildeKinematics::FILightInvertedTildeKinematics() {}
FILightInvertedTildeKinematics::~FILightInvertedTildeKinematics() {}
IBPtr FILightInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FILightInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FILightInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
double mapping = 1.0;
pair<Energy,double> ptz = generatePtZ(mapping,r);
if ( mapping == 0.0 ) {
jacobian(0.0);
return false;
}
Energy pt = ptz.first;
double z = ptz.second;
double y = sqr(pt/lastScale())/(z*(1.-z));
double x = 1./(1.+y);
if ( x < spectatorX() || x > 1. ||
z < 0. || z > 1. ) {
jacobian(0.0);
return false;
}
+ // This should (and does) have a factor of 1/x relative to
+ // the dipole shower jacobian.
mapping /= z*(1.-z);
jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
double phi = 2.*Constants::pi*r[2];
Lorentz5Momentum kt
= getKt(spectator,emitter,pt,phi,true);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = z;
realEmitterMomentum() = z*emitter + (1.-z)*((1.-x)/x)*spectator + kt;
realEmissionMomentum() = (1.-z)*emitter + z*((1.-x)/x)*spectator - kt;
realSpectatorMomentum() = (1./x)*spectator;
realEmitterMomentum().setMass(ZERO);
realEmitterMomentum().rescaleEnergy();
realEmissionMomentum().setMass(ZERO);
realEmissionMomentum().rescaleEnergy();
realSpectatorMomentum().setMass(ZERO);
realSpectatorMomentum().rescaleEnergy();
return true;
}
Energy FILightInvertedTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
return scale * sqrt(z*(1.-z)*(1.-x)/x);
}
double FILightInvertedTildeKinematics::lastZ() const {
return subtractionParameters()[1];
}
Energy FILightInvertedTildeKinematics::ptMax() const {
double x = spectatorX();
return sqrt((1.-x)/x)*lastScale()/2.;
}
pair<double,double> FILightInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
if(pt>hardPt) return make_pair(0.5,0.5);
double s = sqrt(1.-sqr(pt/hardPt));
return make_pair(0.5*(1.-s),0.5*(1.+s));
}
void FILightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FILightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FILightInvertedTildeKinematics::Init() {
static ClassDocumentation<FILightInvertedTildeKinematics> documentation
("FILightInvertedTildeKinematics inverts the final-initial tilde "
"kinematics.");
}
// *** 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<FILightInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigFILightInvertedTildeKinematics("Herwig::FILightInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc
@@ -1,179 +1,182 @@
// -*- C++ -*-
//
// FIMassiveInvertedTildeKinematics.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 FIMassiveInvertedTildeKinematics class.
//
#include "FIMassiveInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FIMassiveInvertedTildeKinematics::FIMassiveInvertedTildeKinematics() {}
FIMassiveInvertedTildeKinematics::~FIMassiveInvertedTildeKinematics() {}
IBPtr FIMassiveInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FIMassiveInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FIMassiveInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
double mapping = 1.0;
pair<Energy,double> ptz = generatePtZ(mapping,r);
if ( mapping == 0.0 ) {
jacobian(0.0);
return false;
}
Energy pt = ptz.first;
double z = ptz.second;
Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
Energy2 m2 = sqr(realEmissionData()->hardProcessMass());
Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
Energy2 scale=2.*emitter*spectator;
double y = (pt*pt+(1.-z)*mi2+z*m2-z*(1.-z)*Mi2) / (z*(1.-z)*scale);
double x = 1./(1.+y);
if ( x < spectatorX() ) {
jacobian(0.0);
return false;
}
- // no additional massive factors
+ // SW 05/12/2016: Checked this is correct
+ // This should appear to have a factor of 1/x relative
+ // to the dipole shower jacobian. It is cancelled by ratios of
+ // real and born cross sections in the units.
mapping /= z*(1.-z);
jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
double phi = 2.*Constants::pi*r[2];
Lorentz5Momentum kt = getKt(spectator,emitter,pt,phi,true);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = z;
realEmitterMomentum() = z*emitter +
(sqr(pt)+mi2-z*z*Mi2)/(z*scale)*spectator + kt;
realEmissionMomentum() = (1.-z)*emitter +
(pt*pt+m2-sqr(1.-z)*Mi2)/((1.-z)*scale)*spectator - kt;
realSpectatorMomentum() = (1.+y)*spectator;
double mui2 = x*mi2/scale;
double mu2 = x*m2/scale;
double Mui2 = x*Mi2/scale;
double xp = 1. + Mui2 - sqr(sqrt(mui2)+sqrt(mu2));
double zm = .5*( 1.-x+Mui2+mui2-mui2 -
sqrt( sqr(1.-x+Mui2-mui2-mu2)-4.*mui2*mu2 ) ) / (1.-x+Mui2);
double zp = .5*( 1.-x+Mui2+mui2-mui2 +
sqrt( sqr(1.-x+Mui2-mui2-mu2)-4.*mui2*mu2 ) ) / (1.-x+Mui2);
if ( x > xp || z < zm || z > zp ) {
jacobian(0.0);
return false;
}
realEmitterMomentum().setMass(sqrt(mi2));
realEmitterMomentum().rescaleEnergy();
realEmissionMomentum().setMass(sqrt(m2));
realEmissionMomentum().rescaleEnergy();
realSpectatorMomentum().setMass(ZERO);
realSpectatorMomentum().rescaleEnergy();
return true;
}
Energy FIMassiveInvertedTildeKinematics::lastPt() const {
Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
Energy2 m2 = sqr(realEmissionData()->hardProcessMass());
Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
Energy2 scale = Mi2 - (realEmitterMomentum()+realEmissionMomentum()-realSpectatorMomentum()).m2();
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
return sqrt( z*(1.-z)*(1.-x)/x*scale -
((1.-z)*mi2+z*m2-z*(1.-z)*Mi2) );
}
double FIMassiveInvertedTildeKinematics::lastZ() const {
return subtractionParameters()[1];
}
Energy FIMassiveInvertedTildeKinematics::ptMax() const {
Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
Energy2 m2 = sqr(realEmissionData()->hardProcessMass());
Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
double x = spectatorX();
// s^star/x
Energy2 scale=2.*bornEmitterMomentum()*bornSpectatorMomentum();
Energy2 s = scale * (1.-x)/x + Mi2;
Energy ptmax = .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
return ptmax > 0.*GeV ? ptmax : 0.*GeV;
}
pair<double,double> FIMassiveInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
if(pt>hardPt) return make_pair(0.5,0.5);
Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
Energy2 m2 = sqr(realEmissionData()->hardProcessMass());
Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
// s^star/x
Energy2 scale=2.*bornEmitterMomentum()*bornSpectatorMomentum();
Energy2 s = scale * (1.-spectatorX())/spectatorX() + Mi2;
double zm = .5*( 1.+(mi2-m2)/s - rootOfKallen(s/s,mi2/s,m2/s) *
sqrt( 1.-sqr(pt/hardPt) ) );
double zp = .5*( 1.+(mi2-m2)/s + rootOfKallen(s/s,mi2/s,m2/s) *
sqrt( 1.-sqr(pt/hardPt) ) );
return make_pair(zm, zp);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMassiveInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FIMassiveInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FIMassiveInvertedTildeKinematics::Init() {
static ClassDocumentation<FIMassiveInvertedTildeKinematics> documentation
("FIMassiveInvertedTildeKinematics inverts the final-initial tilde "
"kinematics.");
}
// *** 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<FIMassiveInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigFIMassiveInvertedTildeKinematics("Herwig::FIMassiveInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
@@ -1,149 +1,150 @@
// -*- C++ -*-
//
// IFLightInvertedTildeKinematics.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 IFLightInvertedTildeKinematics class.
//
#include "IFLightInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IFLightInvertedTildeKinematics::IFLightInvertedTildeKinematics() {}
IFLightInvertedTildeKinematics::~IFLightInvertedTildeKinematics() {}
IBPtr IFLightInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IFLightInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool IFLightInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
double mapping = 1.0;
pair<Energy,double> ptz = generatePtZ(mapping,r);
if ( mapping == 0.0 ) {
jacobian(0.0);
return false;
}
Energy pt = ptz.first;
double z = ptz.second;
double ratio = sqr(pt/lastScale());
double rho = 1. - 4.*ratio*z*(1.-z) / sqr(1. - z + ratio);
if ( rho < 0. ) {
jacobian(0.0);
return false;
}
double x = 0.5*(1./ratio)*(1.-z+ratio)*(1.-sqrt(rho));
double u = 0.5*(1./(1.-z))*(1.-z+ratio)*(1.-sqrt(rho));
if ( x < emitterX() || x > 1. ||
u < 0. || u > 1. ) {
jacobian(0.0);
return false;
}
+ // This jacobian is (1/x^2)*dx*du
mapping *= (1.-x)/((1.-z)*(z*(1.-z)+sqr(x-z)));
jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
double phi = 2.*Constants::pi*r[2];
Lorentz5Momentum kt = getKt(emitter,spectator,pt,phi,true);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = u;
realEmitterMomentum() = (1./x)*emitter;
realEmissionMomentum() = ((1.-x)*(1.-u)/x)*emitter + u*spectator + kt;
realSpectatorMomentum() = ((1.-x)*u/x)*emitter + (1.-u)*spectator - kt;
realEmitterMomentum().setMass(ZERO);
realEmitterMomentum().rescaleEnergy();
realEmissionMomentum().setMass(ZERO);
realEmissionMomentum().rescaleEnergy();
realSpectatorMomentum().setMass(ZERO);
realSpectatorMomentum().rescaleEnergy();
return true;
}
Energy IFLightInvertedTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
return scale * sqrt(u*(1.-u)*(1.-x)/x);
}
Energy IFLightInvertedTildeKinematics::ptMax() const {
double x = emitterX();
return sqrt((1.-x)/x)*lastScale()/2.;
}
pair<double,double> IFLightInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
if(pt>hardPt) return make_pair(0.5,0.5);
double s = sqrt(1.-sqr(pt/hardPt));
double x = emitterX();
return make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
}
double IFLightInvertedTildeKinematics::lastZ() const {
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
return 1. - (1.-x)*(1.-u);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IFLightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void IFLightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void IFLightInvertedTildeKinematics::Init() {
static ClassDocumentation<IFLightInvertedTildeKinematics> documentation
("IFLightInvertedTildeKinematics inverts the initial-final tilde "
"kinematics.");
}
// *** 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<IFLightInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigIFLightInvertedTildeKinematics("Herwig::IFLightInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc
@@ -1,153 +1,166 @@
// -*- C++ -*-
//
// IFMassiveInvertedTildeKinematics.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 IFMassiveInvertedTildeKinematics class.
//
#include "IFMassiveInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IFMassiveInvertedTildeKinematics::IFMassiveInvertedTildeKinematics() {}
IFMassiveInvertedTildeKinematics::~IFMassiveInvertedTildeKinematics() {}
IBPtr IFMassiveInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IFMassiveInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool IFMassiveInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
+ // Compute dipole scale
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
Energy2 scale = 2.*(spectator*emitter);
+ // Generate pt and z
double mapping = 1.0;
pair<Energy,double> ptz = generatePtZ(mapping,r);
if ( mapping == 0.0 ){
jacobian(0.0);
return false;
}
Energy pt = ptz.first;
double z = ptz.second;
- double ratio = sqr(pt)/scale;
+
+ // Compute x and u
+ double ratio = sqr(pt)/scale;
+ double muk2 = sqr(bornSpectatorData()->hardProcessMass())/scale;
+ double rho = 1. - 4.*ratio*(1.-muk2)*z*(1.-z)/sqr(1.-z+ratio);
- // x
- double u = ratio/(1.-z);
- double x = (z*(1.-z)-ratio)/(1.-z-ratio);
- double up = (1.-x) / (1.-x+(x*sqr(bornSpectatorData()->hardProcessMass())/scale));
+ double x = 0.5*((1.-z+ratio)/(ratio*(1.-muk2))) * (1. - sqrt(rho));
+ double u = x*ratio / (1.-z);
+
+ // Following Catani-Seymour paper
+ double muk2CS = x*muk2;
+ double up = (1.-x) /
+ ( 1.-x + muk2CS );
+
+ if ( x < emitterX() || x > 1. ||
+ u < 0. || u > up ) {
+ jacobian(0.0);
+ return false;
+ }
+
- if ( x < emitterX() || x > 1. || u > up) {
- jacobian(0.0);
- return false;
+ // Store x and u
+ subtractionParameters().resize(2);
+ subtractionParameters()[0] = x;
+ subtractionParameters()[1] = u;
+
+ // The jacobian here is the single particle phasespace
+ // saj*(1./x^2)*dx*du
+ // Note - lastScale() is not equal to scale!!!!!!!
+ double jac = abs( (1.+x*(muk2-1.))*(-u*(1.-u)/sqr(x)) - (1.+u*(muk2-1.))*((1.-2.*u)*(1.-x)/x - 2.*u*muk2) );
+ mapping /= x*x*jac;
+ jacobian( mapping*(sqr(lastScale())/sHat()) / (16.*sqr(Constants::pi)) );
+
+ // Compute the new momenta
+ double phi = 2.*Constants::pi*r[2];
+ Lorentz5Momentum kt = getKt(emitter,spectator,pt,phi,true);
+
+ realEmitterMomentum() = (1./x)*emitter;
+ realEmissionMomentum() = ((1.-x)*(1.-u)/x - 2.*u*muk2)*emitter + u*spectator + kt;
+ realSpectatorMomentum() = ((1.-x)*u/x + 2.*u*muk2)*emitter + (1.-u)*spectator - kt;
+
+ realEmitterMomentum().setMass(ZERO);
+ realEmitterMomentum().rescaleEnergy();
+ realEmissionMomentum().setMass(ZERO);
+ realEmissionMomentum().rescaleEnergy();
+ realSpectatorMomentum().setMass(bornSpectatorData()->hardProcessMass());
+ realSpectatorMomentum().rescaleEnergy();
+ return true;
+
}
- pt = sqrt(scale*u*(1.-u)*(1.-x));
- Energy magKt = sqrt(scale*u*(1.-u)*(1.-x)/x - sqr(u*bornSpectatorData()->hardProcessMass()));
-
- // TODO: why not mapping /= sqr(z*(1.-z)-ratio)/(1.-z) ? (see phd thesis, (5.74))
- mapping /= sqr(z*(1.-z)-ratio)/(1.-z-ratio);
- // mapping *= (1.-u)/(1.-2.*u+u*u*alpha)/x;
- jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
-
- double phi = 2.*Constants::pi*r[2];
- Lorentz5Momentum kt = getKt(emitter,spectator,magKt,phi,true);
- subtractionParameters().resize(2);
- subtractionParameters()[0] = x;
- subtractionParameters()[1] = u;
-
- realEmitterMomentum() = (1./x)*emitter;
- realEmissionMomentum() = (-kt*kt-u*u*sqr(bornSpectatorData()->hardProcessMass()))/(u*scale)*emitter +
- u*spectator - kt;
- realSpectatorMomentum() = (-kt*kt+sqr(bornSpectatorData()->hardProcessMass())*u*(2.-u))/((1.-u)*scale)*emitter +
- (1.-u)*spectator + kt;
-
- realEmitterMomentum().setMass(ZERO);
- realEmitterMomentum().rescaleEnergy();
- realEmissionMomentum().setMass(ZERO);
- realEmissionMomentum().rescaleEnergy();
- realSpectatorMomentum().setMass(bornSpectatorData()->hardProcessMass());
- realSpectatorMomentum().rescaleEnergy();
- return true;
-
-}
-
Energy IFMassiveInvertedTildeKinematics::lastPt() const {
- Energy2 scale = 2.*(bornEmitterMomentum()*bornSpectatorMomentum());
- double x = subtractionParameters()[0];
- double u = subtractionParameters()[1];
- // TODO: can't make sense of the following comment
- // there was no factor 1/x in massless case >> check
- // return scale * sqrt(u*(1.-u)*(1.-x)/x);
- return sqrt(scale*u*(1.-u)*(1.-x));
+ Energy2 scale = 2.*(bornEmitterMomentum()*bornSpectatorMomentum());
+ double muk2 = sqr(bornSpectatorData()->hardProcessMass())/scale;
+ double x = subtractionParameters()[0];
+ double u = subtractionParameters()[1];
+ return sqrt(scale * ( u*(1.-u)*(1.-x)/x - u*u*muk2 ));
}
double IFMassiveInvertedTildeKinematics::lastZ() const {
- double x = subtractionParameters()[0];
- double u = subtractionParameters()[1];
- return 1. - (1.-x)*(1.-u);
+ Energy2 scale = 2.*(bornEmitterMomentum()*bornSpectatorMomentum());
+ double muk2 = sqr(bornSpectatorData()->hardProcessMass())/scale;
+ double x = subtractionParameters()[0];
+ double u = subtractionParameters()[1];
+ return u + x + u*x*(muk2-1.);
}
Energy IFMassiveInvertedTildeKinematics::ptMax() const {
- Energy2 scale = 2.*(bornEmitterMomentum()*bornSpectatorMomentum());
- double x = emitterX();
- return sqrt(scale*(1.-x))/2.;
+ double xe = emitterX();
+ Energy2 scale = 2.*(bornEmitterMomentum()*bornSpectatorMomentum());
+ Energy2 A = scale*(1.-xe)/xe;
+ Energy2 mk2 = sqr(bornSpectatorData()->hardProcessMass());
+ return 0.5*A/sqrt(mk2+A);
}
pair<double,double> IFMassiveInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
if(pt>hardPt) return make_pair(0.5,0.5);
double s = sqrt(1.-sqr(pt/hardPt));
- double x = emitterX();
- return make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
+ double xe = emitterX();
+ return make_pair(0.5*(1.+xe-(1.-xe)*s),0.5*(1.+xe+(1.-xe)*s));
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IFMassiveInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void IFMassiveInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void IFMassiveInvertedTildeKinematics::Init() {
static ClassDocumentation<IFMassiveInvertedTildeKinematics> documentation
("IFMassiveInvertedTildeKinematics inverts the initial-final tilde "
"kinematics.");
}
// *** 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<IFMassiveInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigIFMassiveInvertedTildeKinematics("Herwig::IFMassiveInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.cc
@@ -1,123 +1,124 @@
// -*- C++ -*-
//
// IFMassiveTildeKinematics.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 IFMassiveTildeKinematics class.
//
#include "IFMassiveTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IFMassiveTildeKinematics::IFMassiveTildeKinematics() {}
IFMassiveTildeKinematics::~IFMassiveTildeKinematics() {}
IBPtr IFMassiveTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IFMassiveTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool IFMassiveTildeKinematics::doMap() {
Lorentz5Momentum emitter = realEmitterMomentum();
Lorentz5Momentum emission = realEmissionMomentum();
Lorentz5Momentum spectator = realSpectatorMomentum();
double x =
(- emission*spectator + emitter*spectator + emitter*emission) /
(emitter*emission + emitter*spectator);
double u = emitter*emission / (emitter*emission + emitter*spectator);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = u;
bornEmitterMomentum() = x*emitter;
bornSpectatorMomentum() = spectator + emission - (1.-x)*emitter;
bornEmitterMomentum().setMass(ZERO);
bornEmitterMomentum().rescaleEnergy();
bornSpectatorMomentum().setMass(bornSpectatorData()->hardProcessMass());
bornSpectatorMomentum().rescaleEnergy();
return true;
}
Energy IFMassiveTildeKinematics::lastPt() const {
-
- Energy scale = sqrt(2.*(realEmissionMomentum()*realEmitterMomentum()-realEmissionMomentum()*realSpectatorMomentum()+realEmitterMomentum()*realSpectatorMomentum()));
+ Energy2 scale = 2.*(realEmissionMomentum()*realEmitterMomentum()
+-realEmissionMomentum()*realSpectatorMomentum()
++realEmitterMomentum()*realSpectatorMomentum());
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
- return scale * sqrt(u*(1.-u)*(1.-x));
+ double muk2 = sqr(bornSpectatorData()->hardProcessMass())/scale;
+ return sqrt(scale * ( u*(1.-u)*(1.-x)/x - u*u*muk2 ));
+ }
+
+Energy IFMassiveTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
+ Energy2 scale = 2.*(emission*emitter-emission*spectator+emitter*spectator);
+ double x = 0.5*scale / (emitter*emission + emitter*spectator);
+ double u = emitter*emission / (emitter*emission + emitter*spectator);
+
+ double muk2 = sqr(spectator.mass())/scale;
+ return sqrt(scale * ( u*(1.-u)*(1.-x)/x - u*u*muk2 ));
+ }
+
+pair<double,double> IFMassiveTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
+ if(pt>hardPt) return make_pair(0.5,0.5);
+ double s = sqrt(1.-sqr(pt/hardPt));
+ double xe = emitterX();
+ return make_pair(0.5*(1.+xe-(1.-xe)*s),0.5*(1.+xe+(1.-xe)*s));
}
-Energy IFMassiveTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
-
- Energy scale = sqrt(2.*(emission*emitter-emission*spectator+emitter*spectator));
- double x =
- (- emission*spectator + emitter*spectator + emitter*emission) /
- (emitter*emission + emitter*spectator);
- double u = emitter*emission / (emitter*emission + emitter*spectator);
-
- return scale * sqrt(u*(1.-u)*(1.-x));
-}
-
-pair<double,double> IFMassiveTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
-
- if(pt>hardPt) return make_pair(0.5,0.5);
- double s = sqrt(1.-sqr(pt/hardPt));
- double x = emitterX();
- return make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
-}
-
-
-
-
double IFMassiveTildeKinematics::lastZ() const {
+Energy2 scale = 2.*(realEmissionMomentum()*realEmitterMomentum()
+-realEmissionMomentum()*realSpectatorMomentum()
++realEmitterMomentum()*realSpectatorMomentum());
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
- return 1. - (1.-x)*(1.-u);
-}
+ double muk2 = sqr(bornSpectatorData()->hardProcessMass())/scale;
+ return u + x - u*x*(1.-muk2);
+ }
+
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IFMassiveTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void IFMassiveTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void IFMassiveTildeKinematics::Init() {
static ClassDocumentation<IFMassiveTildeKinematics> documentation
("IFMassiveTildeKinematics implements the 'tilde' kinematics for "
"a initial-final subtraction dipole.");
}
// *** 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<IFMassiveTildeKinematics,TildeKinematics>
describeHerwigIFMassiveTildeKinematics("Herwig::IFMassiveTildeKinematics", "Herwig.so");

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:14 PM (2 h, 25 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4476299
Default Alt Text
(47 KB)

Event Timeline