Page MenuHomeHEPForge
Files F7878347

No OneTemporary
Actions

View Options

diff --git a/Makefile.am b/Makefile.am
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,29 +1,29 @@
-SUBDIRS = include \
+SUBDIRS = include MatrixElement Shower \
Utilities PDT Decay PDF Models \
- Shower Hadronization MatrixElement \
+ Hadronization \
UnderlyingEvent Analysis Looptools Sampling \
lib src Doc Contrib Tests
EXTRA_DIST = GUIDELINES
## DISTCHECK_CONFIGURE_FLAGS = --enable-debug --with-thepeg=$(THEPEGPATH) --with-gsl=$(GSLPATH)
ACLOCAL_AMFLAGS = -I m4
DISTCLEANFILES = config.herwig
libclean:
find . -name '*.la' -print0 | xargs -0 rm -rf
cd lib && $(MAKE) $(AM_MAKEFLAGS) clean
cd src && $(MAKE) $(AM_MAKEFLAGS) clean
tests:
cd Tests && $(MAKE) $(AM_MAKEFLAGS) tests
## ThePEG registration
unregister:
cd src && $(MAKE) $(AM_MAKEFLAGS) unregister
register:
cd src && $(MAKE) $(AM_MAKEFLAGS) register
diff --git a/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.cc
@@ -1,98 +1,114 @@
// -*- C++ -*-
//
// FFLightTildeKinematics.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 FFLightTildeKinematics class.
//
#include "FFLightTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FFLightTildeKinematics::FFLightTildeKinematics() {}
FFLightTildeKinematics::~FFLightTildeKinematics() {}
IBPtr FFLightTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FFLightTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FFLightTildeKinematics::doMap() {
Lorentz5Momentum emitter = realEmitterMomentum();
Lorentz5Momentum emission = realEmissionMomentum();
Lorentz5Momentum spectator = realSpectatorMomentum();
double y = emission*emitter / (emission*emitter + emission*spectator + emitter*spectator);
double z = emitter*spectator / (emitter*spectator + emission*spectator);
subtractionParameters().resize(2);
subtractionParameters()[0] = y;
subtractionParameters()[1] = z;
bornEmitterMomentum() = emitter+emission-(y/(1.-y))*spectator;
bornSpectatorMomentum() = spectator/(1.-y);
bornEmitterMomentum().setMass(ZERO);
bornEmitterMomentum().rescaleEnergy();
bornSpectatorMomentum().setMass(ZERO);
bornSpectatorMomentum().rescaleEnergy();
return true;
}
Energy FFLightTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
return scale * sqrt(y*z*(1.-z));
}
+
+Energy FFLightTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
+ Energy scale = (emitter+emission+spectator).m();
+ double y = emission*emitter/(emission*emitter + emission*spectator + emitter*spectator);
+ double z = emitter*spectator / (emitter*spectator + emission*spectator);
+ Energy ret = scale * sqrt( y * z*(1.-z) );
+ return ret;
+}
+
+pair<double,double> FFLightTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
+ 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));
+}
+
+
double FFLightTildeKinematics::lastZ() const {
return subtractionParameters()[1];
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFLightTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FFLightTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FFLightTildeKinematics::Init() {
static ClassDocumentation<FFLightTildeKinematics> documentation
("FFLightTildeKinematics implements the 'tilde' kinematics for "
"a final-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<FFLightTildeKinematics,TildeKinematics>
describeHerwigFFLightTildeKinematics("Herwig::FFLightTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h
@@ -1,134 +1,144 @@
// -*- C++ -*-
//
// FFLightTildeKinematics.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_FFLightTildeKinematics_H
#define HERWIG_FFLightTildeKinematics_H
//
// This is the declaration of the FFLightTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief FFLightTildeKinematics implements the 'tilde' kinematics for
* a final-final subtraction dipole.
*
*/
class FFLightTildeKinematics: public TildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FFLightTildeKinematics();
/**
* The destructor.
*/
virtual ~FFLightTildeKinematics();
//@}
public:
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap();
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const;
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const ;
/**
+ * Given a pt, return the boundaries on z
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const;
+
+ /**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() const;
/*
* True if phase space point is above the alpha cut for this dipole.
*/
bool aboveAlpha() const {return dipole()->alpha()<subtractionParameters()[0];}
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.
*/
FFLightTildeKinematics & operator=(const FFLightTildeKinematics &);
};
}
#endif /* HERWIG_FFLightTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.cc
@@ -1,122 +1,167 @@
// -*- C++ -*-
//
// FFMassiveTildeKinematics.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 FFMassiveTildeKinematics class.
//
#include "FFMassiveTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FFMassiveTildeKinematics::FFMassiveTildeKinematics() {}
FFMassiveTildeKinematics::~FFMassiveTildeKinematics() {}
IBPtr FFMassiveTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FFMassiveTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FFMassiveTildeKinematics::doMap() {
Lorentz5Momentum emitter = realEmitterMomentum();
Lorentz5Momentum emission = realEmissionMomentum();
Lorentz5Momentum spectator = realSpectatorMomentum();
double y = emission*emitter / (emission*emitter + emission*spectator + emitter*spectator);
double z = emitter*spectator / (emitter*spectator + emission*spectator);
subtractionParameters().resize(2);
subtractionParameters()[0] = y;
subtractionParameters()[1] = z;
Lorentz5Momentum pTot = emitter+emission+spectator;
Energy scale = pTot.m();
// 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 );
double bar = 1.-mui2-mu2-muj2;
// from Catani,Seymour,Dittmaier,Trocsanyi (CSm!=0) (5.9)
// has the right massless limit
bornSpectatorMomentum() =
rootOfKallen( 1., Mui2, Muj2 ) / rootOfKallen( 1., mui2+mu2+y*bar, muj2 ) *
( spectator - (pTot*spectator)/sqr(scale)*pTot ) +
0.5*( 1.+Muj2-Mui2 ) * pTot;
bornEmitterMomentum() = pTot - bornSpectatorMomentum();
bornEmitterMomentum().setMass( sqrt(Mui2)*scale );
bornEmitterMomentum().rescaleEnergy();
bornSpectatorMomentum().setMass( sqrt(Muj2)*scale );
bornSpectatorMomentum().rescaleEnergy();
return true;
}
Energy FFMassiveTildeKinematics::lastPt() const {
Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m();
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
// masses
double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
double mu2 = sqr( realEmissionData()->hardProcessMass() / scale );
double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
Energy ret = scale * sqrt( y * (1.-mui2-mu2-muj2) * z*(1.-z) - sqr(1.-z)*mui2 - sqr(z)*mu2 );
return ret;
}
+
+Energy FFMassiveTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
+
+
+ Energy scale = (emitter+emission+spectator).m();
+
+ double y = emission*emitter/(emission*emitter + emission*spectator + emitter*spectator);
+ double z = emitter*spectator / (emitter*spectator + emission*spectator);
+
+ // masses
+ double mui2 = sqr( emitter.m() / scale );
+ double mu2 = sqr( emission.m() / scale );
+ double muj2 = sqr( spectator.m() / scale );
+
+ Energy ret = scale * sqrt( y * (1.-mui2-mu2-muj2) * z*(1.-z) - sqr(1.-z)*mui2 - sqr(z)*mu2 );
+
+ return ret;
+}
+
+
+ // NOTE: bounds calculated at this step may be too loose
+pair<double,double> FFMassiveTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
+
+ 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);
+}
+
+
+
double FFMassiveTildeKinematics::lastZ() const {
return subtractionParameters()[1];
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFMassiveTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FFMassiveTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FFMassiveTildeKinematics::Init() {
static ClassDocumentation<FFMassiveTildeKinematics> documentation
("FFMassiveTildeKinematics implements the 'tilde' kinematics for "
"a final-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<FFMassiveTildeKinematics,TildeKinematics>
describeHerwigFFMassiveTildeKinematics("Herwig::FFMassiveTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h
@@ -1,137 +1,148 @@
// -*- C++ -*-
//
// FFMassiveTildeKinematics.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_FFMassiveTildeKinematics_H
#define HERWIG_FFMassiveTildeKinematics_H
//
// This is the declaration of the FFMassiveTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief FFMassiveTildeKinematics implements the 'tilde' kinematics for
* a final-final subtraction dipole.
*
*/
class FFMassiveTildeKinematics: public TildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FFMassiveTildeKinematics();
/**
* The destructor.
*/
virtual ~FFMassiveTildeKinematics();
//@}
public:
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap();
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const;
+
+ /**
+ * Given a pt, return the boundaries on z
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const;
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const ;
+
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() 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 ) ); }
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.
*/
FFMassiveTildeKinematics & operator=(const FFMassiveTildeKinematics &);
};
}
#endif /* HERWIG_FFMassiveTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h
@@ -1,138 +1,141 @@
// -*- C++ -*-
//
// FILightInvertedTildeKinematics.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_FILightInvertedTildeKinematics_H
#define HERWIG_FILightInvertedTildeKinematics_H
//
// This is the declaration of the FILightInvertedTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief FILightInvertedTildeKinematics inverts the final-final tilde
* kinematics.
*
*/
class FILightInvertedTildeKinematics: public Herwig::InvertedTildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FILightInvertedTildeKinematics();
/**
* The destructor.
*/
virtual ~FILightInvertedTildeKinematics();
//@}
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
*/
virtual pair<double,double> zBounds(Energy pt, Energy hardPt = ZERO) 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.
*/
FILightInvertedTildeKinematics & operator=(const FILightInvertedTildeKinematics &);
};
}
#endif /* HERWIG_FILightInvertedTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.cc
@@ -1,99 +1,121 @@
// -*- C++ -*-
//
// FILightTildeKinematics.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 FILightTildeKinematics class.
//
#include "FILightTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FILightTildeKinematics::FILightTildeKinematics() {}
FILightTildeKinematics::~FILightTildeKinematics() {}
IBPtr FILightTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FILightTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FILightTildeKinematics::doMap() {
Lorentz5Momentum emitter = realEmitterMomentum();
Lorentz5Momentum emission = realEmissionMomentum();
Lorentz5Momentum spectator = realSpectatorMomentum();
double x =
(- emission*emitter + emission*spectator + emitter*spectator) /
(emitter*spectator + emission*spectator);
double z = emitter*spectator / (emitter*spectator + emission*spectator);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = z;
bornEmitterMomentum() = emitter+emission-(1.-x)*spectator;
bornSpectatorMomentum() = x*spectator;
bornEmitterMomentum().setMass(ZERO);
bornEmitterMomentum().rescaleEnergy();
bornSpectatorMomentum().setMass(ZERO);
bornSpectatorMomentum().rescaleEnergy();
return true;
}
Energy FILightTildeKinematics::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);
}
+Energy FILightTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
+
+
+ Energy2 scale = - (emitter+emission-spectator).m2();
+
+
+ double x =
+ (- emission*emitter + emission*spectator + emitter*spectator) /
+ (emitter*spectator + emission*spectator);
+ double z = emitter*spectator / (emitter*spectator + emission*spectator);
+
+ return sqrt( z*(1.-z)*(1.-x)/x*scale ) ;
+}
+
+pair<double,double> FILightTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
+ 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));
+}
+
+
+
double FILightTildeKinematics::lastZ() const {
return subtractionParameters()[1];
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FILightTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FILightTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FILightTildeKinematics::Init() {
static ClassDocumentation<FILightTildeKinematics> documentation
("FILightTildeKinematics implements the 'tilde' kinematics for "
"a final-initial 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<FILightTildeKinematics,TildeKinematics>
describeHerwigFILightTildeKinematics("Herwig::FILightTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h
@@ -1,133 +1,144 @@
// -*- C++ -*-
//
// FILightTildeKinematics.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_FILightTildeKinematics_H
#define HERWIG_FILightTildeKinematics_H
//
// This is the declaration of the FILightTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief FILightTildeKinematics implements the 'tilde' kinematics for
* a final-initial subtraction dipole.
*
*/
class FILightTildeKinematics: public TildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FILightTildeKinematics();
/**
* The destructor.
*/
virtual ~FILightTildeKinematics();
//@}
public:
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap();
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const;
+
+ /**
+ * Given a pt, return the boundaries on z
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const;
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const ;
+
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() const;
/*
* True if phase space point is above the alpha cut for this dipole.
*/
bool aboveAlpha() const {return dipole()->alpha()<1.-subtractionParameters()[0];}
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.
*/
FILightTildeKinematics & operator=(const FILightTildeKinematics &);
};
}
#endif /* HERWIG_FILightTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.cc
@@ -1,109 +1,149 @@
// -*- C++ -*-
//
// FIMassiveTildeKinematics.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 FIMassiveTildeKinematics class.
//
#include "FIMassiveTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FIMassiveTildeKinematics::FIMassiveTildeKinematics() {}
FIMassiveTildeKinematics::~FIMassiveTildeKinematics() {}
IBPtr FIMassiveTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FIMassiveTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FIMassiveTildeKinematics::doMap() {
Lorentz5Momentum emitter = realEmitterMomentum();
Lorentz5Momentum emission = realEmissionMomentum();
Lorentz5Momentum spectator = realSpectatorMomentum();
Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
Energy2 m2 = sqr(realEmissionData()->hardProcessMass());
Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
double x =
(- emission*emitter + emission*spectator + emitter*spectator +
0.5*(Mi2-mi2-m2)) /
(emitter*spectator + emission*spectator);
double z = emitter*spectator / (emitter*spectator + emission*spectator);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = z;
bornEmitterMomentum() = emitter+emission-(1.-x)*spectator;
bornSpectatorMomentum() = x*spectator;
bornEmitterMomentum().setMass(bornEmitterData()->hardProcessMass());
bornEmitterMomentum().rescaleEnergy();
bornSpectatorMomentum().setMass(bornSpectatorData()->hardProcessMass());
bornSpectatorMomentum().rescaleEnergy();
return true;
}
Energy FIMassiveTildeKinematics::lastPt() const {
Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
Energy2 m2 = sqr(realEmissionData()->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) );
}
+
+
+Energy FIMassiveTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
+ // g->QQ or Q -> Qg
+ Energy2 Mi2 = emitter.m()==emission.m()?0.*GeV2:max(emitter.m2(),emission.m2());
+ Energy2 mi2 = emitter.m2();
+ Energy2 m2 = emission.m2();
+
+ Energy2 scale = Mi2 - (emitter+emission-spectator).m2();
+
+ double x =
+ (- emission*emitter + emission*spectator + emitter*spectator +
+ 0.5*(Mi2-mi2-m2)) /
+ (emitter*spectator + emission*spectator);
+ double z = emitter*spectator / (emitter*spectator + emission*spectator);
+
+ return sqrt( z*(1.-z)*(1.-x)/x*scale -
+ ((1.-z)*mi2+z*m2-z*(1.-z)*Mi2) );
+}
+
+
+pair<double,double> FIMassiveTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
+ 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);
+
+}
+
+
+
double FIMassiveTildeKinematics::lastZ() const {
return subtractionParameters()[1];
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMassiveTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FIMassiveTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FIMassiveTildeKinematics::Init() {
static ClassDocumentation<FIMassiveTildeKinematics> documentation
("FIMassiveTildeKinematics implements the 'tilde' kinematics for "
"a final-initial 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<FIMassiveTildeKinematics,TildeKinematics>
describeHerwigFIMassiveTildeKinematics("Herwig::FIMassiveTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h
@@ -1,128 +1,144 @@
// -*- C++ -*-
//
// FIMassiveTildeKinematics.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_FIMassiveTildeKinematics_H
#define HERWIG_FIMassiveTildeKinematics_H
//
// This is the declaration of the FIMassiveTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer, Martin Stoll
*
* \brief FIMassiveTildeKinematics implements the 'tilde' kinematics for
* a final-initial subtraction dipole.
*
*/
class FIMassiveTildeKinematics: public TildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FIMassiveTildeKinematics();
/**
* The destructor.
*/
virtual ~FIMassiveTildeKinematics();
//@}
public:
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap();
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const;
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const ;
+
+ /**
+ * Given a pt, return the boundaries on z
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const;
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() 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 ) ); }
/** @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.
*/
FIMassiveTildeKinematics & operator=(const FIMassiveTildeKinematics &);
};
}
#endif /* HERWIG_FIMassiveTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.cc
@@ -1,101 +1,125 @@
// -*- C++ -*-
//
// IFLightTildeKinematics.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 IFLightTildeKinematics class.
//
#include "IFLightTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IFLightTildeKinematics::IFLightTildeKinematics() {}
IFLightTildeKinematics::~IFLightTildeKinematics() {}
IBPtr IFLightTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IFLightTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool IFLightTildeKinematics::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(ZERO);
bornSpectatorMomentum().rescaleEnergy();
return true;
}
Energy IFLightTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
+ assert(abs(scale * sqrt(u*(1.-u)*(1.-x)/x)-realEmissionMomentum().perp())<1.*GeV);
return scale * sqrt(u*(1.-u)*(1.-x)/x);
}
+Energy IFLightTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
+
+ double x =
+ (- emission*spectator + emitter*spectator + emitter*emission) /
+ (emitter*emission + emitter*spectator);
+ double u = emitter*emission / (emitter*emission + emitter*spectator);
+
+
+ Energy scale = sqrt(2.*((x*emitter)*(spectator + emission - (1.-x)*emitter)));
+ assert(scale * sqrt(u*(1.-u)*(1.-x))-emission.perp()<1*GeV);
+
+ return emission.perp();
+}
+
+pair<double,double> IFLightTildeKinematics::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 IFLightTildeKinematics::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 IFLightTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void IFLightTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void IFLightTildeKinematics::Init() {
static ClassDocumentation<IFLightTildeKinematics> documentation
("IFLightTildeKinematics 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<IFLightTildeKinematics,TildeKinematics>
describeHerwigIFLightTildeKinematics("Herwig::IFLightTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h
@@ -1,135 +1,146 @@
// -*- C++ -*-
//
// IFLightTildeKinematics.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_IFLightTildeKinematics_H
#define HERWIG_IFLightTildeKinematics_H
//
// This is the declaration of the IFLightTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief IFLightTildeKinematics implements the 'tilde' kinematics for
* a initial-final subtraction dipole.
*
*/
class IFLightTildeKinematics: public TildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
IFLightTildeKinematics();
/**
* The destructor.
*/
virtual ~IFLightTildeKinematics();
//@}
public:
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap();
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const;
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const ;
+
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() const;
+
+ /**
+ * Given a pt, return the boundaries on z
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const;
/*
* True if phase space point is above the alpha cut for this dipole.
*/
bool aboveAlpha() const {return dipole()->alpha()<subtractionParameters()[1];}
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.
*/
IFLightTildeKinematics & operator=(const IFLightTildeKinematics &);
};
}
#endif /* HERWIG_IFLightTildeKinematics_H */
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,101 +1,125 @@
// -*- 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()));
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
return scale * sqrt(u*(1.-u)*(1.-x));
}
+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);
+
+ assert(scale * sqrt(u*(1.-u)*(1.-x))-emission.perp()<1*GeV);
+
+ return emission.perp();
+}
+
+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 {
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 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");
diff --git a/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h
@@ -1,128 +1,138 @@
// -*- C++ -*-
//
// IFMassiveTildeKinematics.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_IFMassiveTildeKinematics_H
#define HERWIG_IFMassiveTildeKinematics_H
//
// This is the declaration of the IFMassiveTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer, Martin Stoll
*
* \brief IFMassiveTildeKinematics implements the 'tilde' kinematics for
* a initial-final subtraction dipole.
*
*/
class IFMassiveTildeKinematics: public TildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
IFMassiveTildeKinematics();
/**
* The destructor.
*/
virtual ~IFMassiveTildeKinematics();
//@}
public:
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap();
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const;
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const ;
+
+ /**
+ * Given a pt, return the boundaries on z
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const;
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() 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.
*/
IFMassiveTildeKinematics & operator=(const IFMassiveTildeKinematics &);
};
}
#endif /* HERWIG_IFMassiveTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc
@@ -1,113 +1,125 @@
// -*- C++ -*-
//
// IILightTildeKinematics.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 IILightTildeKinematics class.
//
#include "IILightTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IILightTildeKinematics::IILightTildeKinematics() {}
IILightTildeKinematics::~IILightTildeKinematics() {}
IBPtr IILightTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IILightTildeKinematics::fullclone() const {
return new_ptr(*this);
}
Lorentz5Momentum IILightTildeKinematics::transform(const Lorentz5Momentum& k) const {
LorentzMomentum res =
k - 2.*((k*(K+Ktilde)/(K+Ktilde).m2())*(K+Ktilde)-((k*K)/(K.m2()))*Ktilde);
return res;
}
bool IILightTildeKinematics::doMap() {
Lorentz5Momentum emitter = realEmitterMomentum();
Lorentz5Momentum emission = realEmissionMomentum();
Lorentz5Momentum spectator = realSpectatorMomentum();
double x = (emitter*spectator - emitter*emission - spectator*emission)/(emitter*spectator);
double v = (emitter*emission)/(emitter*spectator);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = v;
bornEmitterMomentum() = x * emitter;
bornSpectatorMomentum() = spectator;
bornEmitterMomentum().setMass(ZERO);
bornEmitterMomentum().rescaleEnergy();
bornSpectatorMomentum().setMass(ZERO);
bornSpectatorMomentum().rescaleEnergy();
K = emitter + spectator - emission;
Ktilde = x * emitter + spectator;
return true;
}
Energy IILightTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double v = subtractionParameters()[1];
return scale * sqrt(v*(1.-x-v)/x);
}
+
+Energy IILightTildeKinematics::lastPt(Lorentz5Momentum emitter,Lorentz5Momentum emission,Lorentz5Momentum spectator)const {
+ return emission.perp();
+}
+
+pair<double,double> IILightTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
+ if(pt>hardPt) return make_pair(0.5,0.5);
+ double root = (1.-emitterX())*sqrt(1.-sqr(pt/hardPt));
+ return make_pair(0.5*( 1.+emitterX() - root),0.5*( 1.+emitterX() + root));
+}
+
+
double IILightTildeKinematics::lastZ() const {
double x = subtractionParameters()[0];
double v = subtractionParameters()[1];
return x + v;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IILightTildeKinematics::persistentOutput(PersistentOStream & os) const {
os << ounit(K,GeV) << ounit(Ktilde,GeV);
}
void IILightTildeKinematics::persistentInput(PersistentIStream & is, int) {
is >> iunit(K,GeV) >> iunit(Ktilde,GeV);
}
void IILightTildeKinematics::Init() {
static ClassDocumentation<IILightTildeKinematics> documentation
("IILightTildeKinematics implements the 'tilde' kinematics for "
"a initial-initial 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<IILightTildeKinematics,TildeKinematics>
describeHerwigIILightTildeKinematics("Herwig::IILightTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h
@@ -1,160 +1,170 @@
// -*- C++ -*-
//
// IILightTildeKinematics.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_IILightTildeKinematics_H
#define HERWIG_IILightTildeKinematics_H
//
// This is the declaration of the IILightTildeKinematics class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \breif IILightTildeKinematics implements the 'tilde' kinematics for
* a initial-initial subtraction dipole.
*
*/
class IILightTildeKinematics: public TildeKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
IILightTildeKinematics();
/**
* The destructor.
*/
virtual ~IILightTildeKinematics();
//@}
public:
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap();
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const;
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const ;
+
+ /**
+ * Given a pt, return the boundaries on z
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const;
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() const;
/**
* Return true, if this TildeKinematics object needs to transform
* all other particles in the process except the emitter and spectator
*/
virtual bool doesTransform() const { return true; }
/**
* If this TildeKinematics object needs to transform all other particles
* in the process except the emitter and spectator, return the transformed
* momentum.
*/
virtual Lorentz5Momentum transform(const Lorentz5Momentum& p) const;
/*
* True if phase space point is above the alpha cut for this dipole.
*/
bool aboveAlpha() const {return dipole()->alpha()<subtractionParameters()[1];}
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 K momentum used to transform the final state.
*/
Lorentz5Momentum K;
/**
* The Ktilde momentum used to transform the final state.
*/
Lorentz5Momentum Ktilde;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
IILightTildeKinematics & operator=(const IILightTildeKinematics &);
};
}
#endif /* HERWIG_IILightTildeKinematics_H */
diff --git a/MatrixElement/Matchbox/Phasespace/TildeKinematics.h b/MatrixElement/Matchbox/Phasespace/TildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/TildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/TildeKinematics.h
@@ -1,348 +1,384 @@
// -*- C++ -*-
//
// TildeKinematics.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_TildeKinematics_H
#define HERWIG_TildeKinematics_H
//
// This is the declaration of the TildeKinematics class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief TildeKinematics is the base class for the 'tilde'
* kinematics being used for subtraction terms in the
* formalism of Catani and Seymour.
*
*/
class TildeKinematics: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
TildeKinematics();
/**
* The destructor.
*/
virtual ~TildeKinematics();
//@}
public:
/**
* Clone this object
*/
Ptr<TildeKinematics>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<TildeKinematics>::ptr>(clone());
}
/** @name Access to kinematic quantities. */
//@{
/**
* Return the momentum of the emitter in the real emission process
*/
const Lorentz5Momentum& realEmitterMomentum() const {
return theRealXComb->meMomenta()[theDipole->realEmitter()];
}
/**
* Return the momentum of the emission in the real emission process
*/
const Lorentz5Momentum& realEmissionMomentum() const {
return theRealXComb->meMomenta()[theDipole->realEmission()];
}
/**
* Return the momentum of the spectator in the real emission process
*/
const Lorentz5Momentum& realSpectatorMomentum() const {
return theRealXComb->meMomenta()[theDipole->realSpectator()];
}
/**
* Return the momentum of the emitter in the underlying Born process
*/
const Lorentz5Momentum& bornEmitterMomentum() const { return theBornEmitterMomentum; }
/**
* Return the momentum of the spectator in the underlying Born process
*/
const Lorentz5Momentum& bornSpectatorMomentum() const { return theBornSpectatorMomentum; }
/**
* Return the vector of dimensionless variables calculated
*/
const vector<double>& subtractionParameters() const { return theDipole->subtractionParameters(); }
/**
* Return true, if this TildeKinematics object needs to transform
* all other particles in the process except the emitter and spectator
*/
virtual bool doesTransform() const { return false; }
/**
* If this TildeKinematics object needs to transform all other particles
* in the process except the emitter and spectator, return the transformed
* momentum.
*/
virtual Lorentz5Momentum transform(const Lorentz5Momentum& p) const { return p; }
//@}
/**
* If this tilde kinematics is implementing a mapping different from
* the baseline dipole mapping, determine the relevant shower
* parameters and check for phase space boundaries. Note that real
* emission kinematics only are available at this stage.
*/
virtual void getShowerVariables() const {}
/**
* If this tilde kinematics is implementing a mapping different from
* the baseline dipole mapping, return the ratio of phase space
* factorization Jacobians for this and the nominal dipole
* mapping. This is used for matching subtractions.
*/
virtual double jacobianRatio() const { return 1.; }
public:
/** @name Access to process data. */
//@{
/**
* Prepare given a dipole, and XCombs describing the real emission
* and underlying Born processes, respectively.
*/
void prepare(tcStdXCombPtr newRealXComb,
tcStdXCombPtr newBornXComb) {
theRealXComb = newRealXComb; theBornXComb = newBornXComb;
}
/**
* Set the current dipole
*/
void dipole(Ptr<SubtractionDipole>::tptr dip) { theDipole = dip; }
/**
* Return the current dipole
*/
Ptr<SubtractionDipole>::tptr dipole() { return theDipole; }
/**
* Return the current dipole
*/
Ptr<SubtractionDipole>::tcptr dipole() const { return theDipole; }
/**
* Perform the mapping to the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the tilde
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap() = 0;
/**
* Return the pt associated to the last merged splitting.
*/
virtual Energy lastPt() const = 0;
+
+
+ /**
+ * Return the pt associated to emitter emission and sppectator momentum.
+ */
+ virtual Energy lastPt(Lorentz5Momentum,Lorentz5Momentum,Lorentz5Momentum) const =0 ;
+
+ /**
+ * Given a pt and a hard pt, return the boundaries on z;
+ */
+ virtual pair<double,double> zBounds(Energy pt, Energy hardPt ) const = 0;
+
+
/**
* Return the momentum fraction associated to the last splitting.
*/
virtual double lastZ() const = 0;
/**
* Return the relevant dipole scale
*/
virtual Energy lastScale() const;
virtual bool aboveAlpha() const {
cerr<<"only implemented for light kinematics";
assert(false);
return false;
}
/**
* Return the particle type of the emitter in the real emission process
*/
cPDPtr realEmitterData() const {
return
(theDipole && theRealXComb) ?
theRealXComb->mePartonData()[theDipole->realEmitter()] :
cPDPtr();
}
/**
* Return the particle type of the emission in the real emission process
*/
cPDPtr realEmissionData() const {
return
(theDipole && theRealXComb) ?
theRealXComb->mePartonData()[theDipole->realEmission()] :
cPDPtr();
}
/**
* Return the particle type of the spectator in the real emission process
*/
cPDPtr realSpectatorData() const {
return
(theDipole && theRealXComb) ?
theRealXComb->mePartonData()[theDipole->realSpectator()] :
cPDPtr();
}
/**
* Return the particle type of the emitter in the underlying Born process
*/
cPDPtr bornEmitterData() const {
return
(theDipole && theBornXComb) ?
theBornXComb->mePartonData()[theDipole->bornEmitter()] :
cPDPtr();
}
/**
* Return the particle type of the spectator in the underlying Born process
*/
cPDPtr bornSpectatorData() const {
return
(theDipole && theBornXComb) ?
theBornXComb->mePartonData()[theDipole->bornSpectator()] :
cPDPtr();
}
//@}
protected:
/**
* Access the momentum of the emitter in the underlying Born process
*/
Lorentz5Momentum& bornEmitterMomentum() { return theBornEmitterMomentum; }
/**
* Access the momentum of the spectator in the underlying Born process
*/
Lorentz5Momentum& bornSpectatorMomentum() { return theBornSpectatorMomentum; }
/**
* Access the vector of dimensionless variables calculated
*/
vector<double>& subtractionParameters() { return theDipole->subtractionParameters(); }
+
+ /**
+ * Return the momentum fraction of the emitter
+ */
+ double emitterX() const {
+ return
+ theDipole->bornEmitter() == 0 ?
+ theBornXComb->lastX1() :
+ theBornXComb->lastX2();
+ }
+
+
+ /**
+ * Return the momentum fraction of the spectator
+ */
+ double spectatorX() const {
+ return
+ theDipole->bornSpectator() == 0 ?
+ theBornXComb->lastX1() :
+ theBornXComb->lastX2();
+ }
+
+
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();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans);
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
//@}
private:
/**
* The last dipole this TildeKinematics has been selected for
*/
Ptr<SubtractionDipole>::tptr theDipole;
/**
* The XComb object describing the real emission process
*/
tcStdXCombPtr theRealXComb;
/**
* The XComb object describing the underlying Born process
*/
tcStdXCombPtr theBornXComb;
/**
* The momentum of the emitter in the underlying Born process
*/
Lorentz5Momentum theBornEmitterMomentum;
/**
* The momentum of the spectator in the underlying Born process
*/
Lorentz5Momentum theBornSpectatorMomentum;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
TildeKinematics & operator=(const TildeKinematics &);
};
}
#endif /* HERWIG_TildeKinematics_H */
diff --git a/Shower/Dipole/Base/DipoleSplittingGenerator.cc b/Shower/Dipole/Base/DipoleSplittingGenerator.cc
--- a/Shower/Dipole/Base/DipoleSplittingGenerator.cc
+++ b/Shower/Dipole/Base/DipoleSplittingGenerator.cc
@@ -1,731 +1,731 @@
// -*- C++ -*-
//
// DipoleSplittingGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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 DipoleSplittingGenerator class.
//
#include <config.h>
#include "DipoleSplittingGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
using namespace Herwig;
DipoleSplittingGenerator::DipoleSplittingGenerator()
: HandlerBase(),
theExponentialGenerator(0), prepared(false), presampling(false),
theDoCompensate(false), theSplittingWeight(1.) {
if ( ShowerHandler::currentHandler() )
setGenerator(ShowerHandler::currentHandler()->generator());
}
DipoleSplittingGenerator::~DipoleSplittingGenerator() {
if ( theExponentialGenerator ) {
delete theExponentialGenerator;
theExponentialGenerator = 0;
}
}
IBPtr DipoleSplittingGenerator::clone() const {
return new_ptr(*this);
}
IBPtr DipoleSplittingGenerator::fullclone() const {
return new_ptr(*this);
}
void DipoleSplittingGenerator::wrap(Ptr<DipoleSplittingGenerator>::ptr other) {
assert(!prepared);
theOtherGenerator = other;
}
void DipoleSplittingGenerator::resetVariations() {
for ( map<string,double>::iterator w = currentWeights.begin();
w != currentWeights.end(); ++w )
w->second = 1.;
}
void DipoleSplittingGenerator::veto(const vector<double>&, double p, double r) {
double factor = 1.;
if ( splittingReweight() ) {
if ( ( ShowerHandler::currentHandler()->firstInteraction() && splittingReweight()->firstInteraction() ) ||
( !ShowerHandler::currentHandler()->firstInteraction() && splittingReweight()->secondaryInteractions() ) ) {
factor = splittingReweight()->evaluate(generatedSplitting);
theSplittingWeight *= (r-factor*p)/(r-p);
}
}
splittingKernel()->veto(generatedSplitting, factor*p, r, currentWeights);
}
void DipoleSplittingGenerator::accept(const vector<double>&, double p, double r) {
double factor = 1.;
if ( splittingReweight() ) {
if ( ( ShowerHandler::currentHandler()->firstInteraction() && splittingReweight()->firstInteraction() ) ||
( !ShowerHandler::currentHandler()->firstInteraction() && splittingReweight()->secondaryInteractions() ) ) {
factor = splittingReweight()->evaluate(generatedSplitting);
theSplittingWeight *= factor;
}
}
splittingKernel()->accept(generatedSplitting, factor*p, r, currentWeights);
}
void DipoleSplittingGenerator::prepare(const DipoleSplittingInfo& sp) {
generatedSplitting = sp;
generatedSplitting.splittingKinematics(splittingKernel()->splittingKinematics());
generatedSplitting.splittingParameters().resize(splittingKernel()->nDimAdditional());
if ( wrapping() ) {
generatedSplitting.emitterData(theSplittingKernel->emitter(generatedSplitting.index()));
generatedSplitting.spectatorData(theSplittingKernel->spectator(generatedSplitting.index()));
generatedSplitting.emissionData(theSplittingKernel->emission(generatedSplitting.index()));
parameters.resize(theOtherGenerator->nDim());
prepared = true;
return;
}
generatedSplitting.emitterData(splittingKernel()->emitter(generatedSplitting.index()));
generatedSplitting.spectatorData(splittingKernel()->spectator(generatedSplitting.index()));
generatedSplitting.emissionData(splittingKernel()->emission(generatedSplitting.index()));
presampledSplitting = generatedSplitting;
prepared = true;
parameters.resize(nDim());
theExponentialGenerator =
new exsample::exponential_generator<DipoleSplittingGenerator,UseRandom>();
theExponentialGenerator->sampling_parameters().maxtry = maxtry();
theExponentialGenerator->sampling_parameters().presampling_points = presamplingPoints();
theExponentialGenerator->sampling_parameters().freeze_grid = freezeGrid();
theExponentialGenerator->detuning(detuning());
theExponentialGenerator->docompensate(theDoCompensate);
theExponentialGenerator->function(this);
theExponentialGenerator->initialize();
}
void DipoleSplittingGenerator::fixParameters(const DipoleSplittingInfo& sp,
Energy optHardPt) {
assert(generator());
assert(!presampling);
assert(prepared);
assert(sp.index() == generatedSplitting.index());
generatedSplitting.scale(sp.scale());
parameters[3] = sp.scale()/generator()->maximumCMEnergy();
generatedSplitting.hardPt(sp.hardPt());
parameters[0] = splittingKinematics()->ptToRandom(optHardPt == ZERO ?
generatedSplitting.hardPt() :
min(generatedSplitting.hardPt(),optHardPt),
sp.scale(),
sp.emitterX(), sp.spectatorX(),
generatedSplitting.index(),
*splittingKernel());
size_t shift = 4;
if ( generatedSplitting.index().emitterPDF().pdf() &&
generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.emitterX(sp.emitterX());
generatedSplitting.spectatorX(sp.spectatorX());
parameters[4] = sp.emitterX();
parameters[5] = sp.spectatorX();
shift += 2;
}
if ( generatedSplitting.index().emitterPDF().pdf() &&
!generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.emitterX(sp.emitterX());
parameters[4] = sp.emitterX();
++shift;
}
if ( !generatedSplitting.index().emitterPDF().pdf() &&
generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.spectatorX(sp.spectatorX());
parameters[4] = sp.spectatorX();
++shift;
}
if ( splittingKernel()->nDimAdditional() )
copy(sp.lastSplittingParameters().begin(),sp.lastSplittingParameters().end(),parameters.begin()+shift);
if ( sp.emitter() )
generatedSplitting.emitter(sp.emitter());
if ( sp.spectator() )
generatedSplitting.spectator(sp.spectator());
}
int DipoleSplittingGenerator::nDim() const {
assert(!wrapping());
assert(prepared);
int ret = 4; // 0 pt, 1 z, 2 phi, 3 scale, 4/5 xs + parameters
if ( generatedSplitting.index().emitterPDF().pdf() ) {
++ret;
}
if ( generatedSplitting.index().spectatorPDF().pdf() ) {
++ret;
}
ret += splittingKernel()->nDimAdditional();
return ret;
}
const vector<bool>& DipoleSplittingGenerator::sampleFlags() {
assert(!wrapping());
if ( !theFlags.empty() )
return theFlags;
theFlags.resize(nDim(),false);
theFlags[0] = true; theFlags[1] = true; theFlags[2] = true; // 0 pt, 1 z, 2 phi
return theFlags;
}
const pair<vector<double>,vector<double> >& DipoleSplittingGenerator::support() {
assert(!wrapping());
if ( !theSupport.first.empty() )
return theSupport;
vector<double> lower(nDim(),0.);
vector<double> upper(nDim(),1.);
pair<double,double> kSupport =
generatedSplitting.splittingKinematics()->kappaSupport(generatedSplitting);
pair<double,double> xSupport =
generatedSplitting.splittingKinematics()->xiSupport(generatedSplitting);
lower[0] = kSupport.first;
lower[1] = xSupport.first;
upper[0] = kSupport.second;
upper[1] = xSupport.second;
theSupport.first = lower;
theSupport.second = upper;
return theSupport;
}
void DipoleSplittingGenerator::startPresampling() {
assert(!wrapping());
splittingKernel()->startPresampling(generatedSplitting.index());
presampling = true;
}
void DipoleSplittingGenerator::stopPresampling() {
assert(!wrapping());
splittingKernel()->stopPresampling(generatedSplitting.index());
presampling = false;
}
bool DipoleSplittingGenerator::haveOverestimate() const {
assert(!wrapping());
assert(prepared);
return
generatedSplitting.splittingKinematics()->haveOverestimate() &&
splittingKernel()->haveOverestimate(generatedSplitting);
}
bool DipoleSplittingGenerator::overestimate(const vector<double>& point) {
assert(!wrapping());
assert(prepared);
assert(!presampling);
assert(haveOverestimate());
if ( ! generatedSplitting.splittingKinematics()->generateSplitting(point[0],point[1],point[2],
generatedSplitting,
*splittingKernel()) )
return 0.;
generatedSplitting.splittingKinematics()->prepareSplitting(generatedSplitting);
return
( generatedSplitting.splittingKinematics()->jacobianOverestimate() *
splittingKernel()->overestimate(generatedSplitting) );
}
double DipoleSplittingGenerator::invertOverestimateIntegral(double value) const {
assert(!wrapping());
assert(prepared);
assert(!presampling);
assert(haveOverestimate());
return
splittingKernel()->invertOverestimateIntegral(generatedSplitting,value);
}
//TODO: make some proper flag..
double DipoleSplittingGenerator::evaluate(const vector<double>& point,Energy fixed) {
assert(!wrapping());
assert(prepared);
assert(generator());
DipoleSplittingInfo& split =
( !presampling ? generatedSplitting : presampledSplitting );
if(fixed>0.*GeV){
split.fixedScale(fixed);
}else{
split.fixedScale(-1.*GeV);
}
split.continuesEvolving();
size_t shift = 4;
if ( presampling ) {
split.scale(point[3] * generator()->maximumCMEnergy());
if ( split.index().emitterPDF().pdf() &&
split.index().spectatorPDF().pdf() ) {
split.emitterX(point[4]);
split.spectatorX(point[5]);
shift += 2;
}
if ( split.index().emitterPDF().pdf() &&
!split.index().spectatorPDF().pdf() ) {
split.emitterX(point[4]);
++shift;
}
if ( !split.index().emitterPDF().pdf() &&
split.index().spectatorPDF().pdf() ) {
split.spectatorX(point[4]);
++shift;
}
if ( splittingKernel()->nDimAdditional() )
copy(point.begin()+shift,point.end(),split.splittingParameters().begin());
split.hardPt(split.splittingKinematics()->ptMax(split.scale(),
split.emitterX(),
split.spectatorX(),
split.index(),
*splittingKernel()));
}
if ( ! split.splittingKinematics()->generateSplitting(point[0],point[1],point[2],split,*splittingKernel()) ) {
split.lastValue(0.);
return 0.;
}
split.splittingKinematics()->prepareSplitting(split);
if ( split.stoppedEvolving() ) {
split.lastValue(0.);
return 0.;
}
if ( !presampling )
splittingKernel()->clearAlphaPDFCache();
double kernel = splittingKernel()->evaluate(split);
double jac = split.splittingKinematics()->jacobian();
// multiply in the profile scales when relevant
assert(ShowerHandler::currentHandler());
if ( ShowerHandler::currentHandler()->firstInteraction() &&
ShowerHandler::currentHandler()->profileScales() &&
!presampling ) {
Energy hard = ShowerHandler::currentHandler()->hardScale();
if ( hard > ZERO )
kernel *= ShowerHandler::currentHandler()->profileScales()->hardScaleProfile(hard,split.lastPt());
}
split.lastValue( abs(jac) * kernel );
if ( ! isfinite(split.lastValue()) ) {
generator()->log() << "DipoleSplittingGenerator:evaluate(): problematic splitting kernel encountered for "
<< splittingKernel()->name() << "\n" << flush;
split.lastValue(0.0);
}
if ( kernel < 0. )
return 0.;
return split.lastValue();
}
void DipoleSplittingGenerator::doGenerate(map<string,double>& variations,
Energy optCutoff) {
assert(!wrapping());
double res = 0.;
Energy startPt = generatedSplitting.hardPt();
double optKappaCutoff = 0.0;
if ( optCutoff > splittingKinematics()->IRCutoff() ) {
optKappaCutoff = splittingKinematics()->ptToRandom(optCutoff,
generatedSplitting.scale(),
generatedSplitting.emitterX(),
generatedSplitting.spectatorX(),
generatedSplitting.index(),
*splittingKernel());
}
resetVariations();
theSplittingWeight = 1.;
while (true) {
try {
if ( optKappaCutoff == 0.0 ) {
res = theExponentialGenerator->generate();
} else {
res = theExponentialGenerator->generate(optKappaCutoff);
}
} catch (exsample::exponential_regenerate&) {
resetVariations();
theSplittingWeight = 1.;
generatedSplitting.hardPt(startPt);
continue;
} catch (exsample::hit_and_miss_maxtry&) {
throw DipoleShowerHandler::RedoShower();
} catch (exsample::selection_maxtry&) {
throw DipoleShowerHandler::RedoShower();
}
break;
}
for ( map<string,double>::const_iterator w = currentWeights.begin();
w != currentWeights.end(); ++w ) {
map<string,double>::iterator v = variations.find(w->first);
if ( v != variations.end() )
v->second *= w->second;
else
variations[w->first] = w->second;
}
if ( res == 0. ) {
generatedSplitting.lastPt(0.0*GeV);
generatedSplitting.didStopEvolving();
} else {
generatedSplitting.continuesEvolving();
if ( theMCCheck )
theMCCheck->book(generatedSplitting.emitterX(),
generatedSplitting.spectatorX(),
generatedSplitting.scale(),
startPt,
generatedSplitting.lastPt(),
generatedSplitting.lastZ(),
1.);
}
}
Energy DipoleSplittingGenerator::generate(const DipoleSplittingInfo& split,
map<string,double>& variations,
Energy optHardPt,
Energy optCutoff) {
fixParameters(split,optHardPt);
if ( wrapping() ) {
return theOtherGenerator->generateWrapped(generatedSplitting,variations,optHardPt,optCutoff);
}
doGenerate(variations,optCutoff);
return generatedSplitting.lastPt();
}
double DipoleSplittingGenerator::unlops(const DipoleSplittingInfo& split,Energy down,Energy fixedScale){
fixParameters(split);
if ( wrapping() ) {
return theOtherGenerator->wrappedUnlops( generatedSplitting, down,fixedScale);
}
return dounlops( split, down,fixedScale);
}
double DipoleSplittingGenerator::sudakov(const DipoleSplittingInfo& split,Energy down){
fixParameters(split);
if ( wrapping() ) {
return theOtherGenerator->wrappedSudakov( generatedSplitting, down);
}
return dosudakov( split, down);
}
double DipoleSplittingGenerator::dounlops(const DipoleSplittingInfo& ,Energy down,Energy fixedScale){
assert(down > splittingKinematics()->IRCutoff());
double optKappaCutoff = splittingKinematics()->ptToRandom(down,
generatedSplitting.scale(),
generatedSplitting.emitterX(),
generatedSplitting.spectatorX(),
generatedSplitting.index(),
*splittingKernel());
vector<double> RN;
RN.resize(parameters.size());
double res=0.;
double resq=0.;
double varx=10.;
int k=0;
while ((k<500.||varx>0.1)&&k<5000){
k+=1.;
RN[0]= optKappaCutoff+(1-optKappaCutoff)*UseRandom::rnd();
for (size_t rn=1;rn< parameters.size();rn++)RN[rn]=UseRandom::rnd();
double tmp=(1-optKappaCutoff)*evaluate(RN,fixedScale);
res+= tmp;
resq+=pow(tmp,2.);
if(k%50==0.){
varx=sqrt((resq/pow(1.*k,2)-pow(res,2)/pow(1.*k,3)));
}
}
return -res/(1.0*k);
}
double DipoleSplittingGenerator::dosudakov(const DipoleSplittingInfo& ,Energy down){
double optKappaCutoff = splittingKinematics()->ptToRandom(down,
generatedSplitting.scale(),
generatedSplitting.emitterX(),
generatedSplitting.spectatorX(),
generatedSplitting.index(),
*splittingKernel());
vector<double> RN;
RN.resize(parameters.size());
double res=0.;
double resq=0.;
double var=10.;
double varx=10.;
int k=0;
- while (((k<100.||var>0.05)&&k<50000)){
+ while (((k<100.||var>0.1)&&k<50000)){
k+=1.;
RN[0]= optKappaCutoff+(1-optKappaCutoff)*UseRandom::rnd();
for (size_t rn=1;rn< parameters.size();rn++)RN[rn]=UseRandom::rnd();
double tmp=(1-optKappaCutoff)*evaluate(RN);
res+= tmp;
resq+=pow(tmp,2.);
varx=sqrt((resq/pow(1.*k,2)-pow(res,2)/pow(1.*k,3)));
if(k%50==0.)var= exp(-(res)/(1.0*k)+varx)-exp(-(res)/(1.0*k)-varx);
}
return exp(-res/(1.0*k));
}
double DipoleSplittingGenerator::wrappedUnlops(DipoleSplittingInfo& split,
Energy down,Energy fixedScale) {
assert(!wrapping());
DipoleSplittingInfo backup = generatedSplitting;
generatedSplitting = split;
fixParameters(split);
double res=dounlops( split, down,fixedScale);
split = generatedSplitting;
generatedSplitting = backup;
return res;
}
double DipoleSplittingGenerator::wrappedSudakov(DipoleSplittingInfo& split,
Energy down) {
assert(!wrapping());
DipoleSplittingInfo backup = generatedSplitting;
generatedSplitting = split;
fixParameters(split);
double res=dosudakov( split, down);
split = generatedSplitting;
generatedSplitting = backup;
return res;
}
Energy DipoleSplittingGenerator::generateWrapped(DipoleSplittingInfo& split,
map<string,double>& variations,
Energy optHardPt,
Energy optCutoff) {
assert(!wrapping());
DipoleSplittingInfo backup = generatedSplitting;
generatedSplitting = split;
fixParameters(split,optHardPt);
try {
doGenerate(variations,optCutoff);
} catch (...) {
split = generatedSplitting;
generatedSplitting = backup;
throw;
}
Energy pt = generatedSplitting.lastPt();
split = generatedSplitting;
generatedSplitting = backup;
return pt;
}
void DipoleSplittingGenerator::completeSplitting(DipoleSplittingInfo& sp) const {
pair<bool,bool> conf = sp.configuration();
sp = generatedSplitting;
sp.configuration(conf);
}
Ptr<DipoleSplittingKernel>::tptr DipoleSplittingGenerator::splittingKernel() const {
if ( wrapping() )
return theOtherGenerator->splittingKernel();
return theSplittingKernel;
}
Ptr<DipoleSplittingReweight>::tptr DipoleSplittingGenerator::splittingReweight() const {
if ( wrapping() )
return theOtherGenerator->splittingReweight();
return theSplittingReweight;
}
Ptr<DipoleSplittingKinematics>::tptr DipoleSplittingGenerator::splittingKinematics() const {
if ( wrapping() )
return theOtherGenerator->splittingKinematics();
return theSplittingKernel->splittingKinematics();
}
void DipoleSplittingGenerator::splittingKernel(Ptr<DipoleSplittingKernel>::tptr sp) {
theSplittingKernel = sp;
if ( theSplittingKernel->mcCheck() )
theMCCheck = theSplittingKernel->mcCheck();
}
void DipoleSplittingGenerator::splittingReweight(Ptr<DipoleSplittingReweight>::tptr sp) {
theSplittingReweight = sp;
}
void DipoleSplittingGenerator::debugGenerator(ostream& os) const {
os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";
os << " generating splittings using\n"
<< " splittingKernel = " << splittingKernel()->name()
<< " splittingKinematics = " << generatedSplitting.splittingKinematics()->name() << "\n"
<< " to sample splittings of type:\n";
os << generatedSplitting;
os << "--------------------------------------------------------------------------------\n";
}
void DipoleSplittingGenerator::debugLastEvent(ostream& os) const {
os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";
os << " last generated event:\n";
os << generatedSplitting;
os << "--------------------------------------------------------------------------------\n";
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleSplittingGenerator::persistentOutput(PersistentOStream & os) const {
os << theOtherGenerator << theSplittingKernel << theSplittingReweight << theMCCheck << theDoCompensate;
}
void DipoleSplittingGenerator::persistentInput(PersistentIStream & is, int) {
is >> theOtherGenerator >> theSplittingKernel >> theSplittingReweight >> theMCCheck >> theDoCompensate;
}
ClassDescription<DipoleSplittingGenerator> DipoleSplittingGenerator::initDipoleSplittingGenerator;
// Definition of the static class description member.
void DipoleSplittingGenerator::Init() {
static ClassDocumentation<DipoleSplittingGenerator> documentation
("DipoleSplittingGenerator is used by the dipole shower "
"to sample splittings from a given dipole splitting kernel.");
static Reference<DipoleSplittingGenerator,DipoleSplittingKernel> interfaceSplittingKernel
("SplittingKernel",
"Set the splitting kernel to sample from.",
&DipoleSplittingGenerator::theSplittingKernel, false, false, true, false, false);
static Reference<DipoleSplittingGenerator,DipoleSplittingReweight> interfaceSplittingReweight
("SplittingReweight",
"Set the splitting reweight.",
&DipoleSplittingGenerator::theSplittingReweight, false, false, true, true, false);
static Reference<DipoleSplittingGenerator,DipoleMCCheck> interfaceMCCheck
("MCCheck",
"[debug option] MCCheck",
&DipoleSplittingGenerator::theMCCheck, false, false, true, true, false);
interfaceMCCheck.rank(-1);
}
diff --git a/Shower/Dipole/DipoleShowerHandler.h b/Shower/Dipole/DipoleShowerHandler.h
--- a/Shower/Dipole/DipoleShowerHandler.h
+++ b/Shower/Dipole/DipoleShowerHandler.h
@@ -1,525 +1,530 @@
// -*- C++ -*-
//
// DipoleShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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_DipoleShowerHandler_H
#define HERWIG_DipoleShowerHandler_H
//
// This is the declaration of the DipoleShowerHandler class.
//
#include "Herwig/Shower/ShowerHandler.h"
#include "Herwig/Shower/Dipole/DipoleShowerHandler.fh"
#include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h"
#include "Herwig/Shower/Dipole/Base/DipoleSplittingReweight.h"
#include "Herwig/Shower/Dipole/Kernels/DipoleSplittingKernel.h"
#include "Herwig/Shower/Dipole/Base/DipoleSplittingGenerator.h"
#include "Herwig/Shower/Dipole/Base/DipoleEventRecord.h"
#include "Herwig/Shower/Dipole/Base/DipoleEvolutionOrdering.h"
#include "Herwig/Shower/Dipole/Base/DipoleEventReweight.h"
#include "Herwig/Shower/Dipole/Utility/ConstituentReshuffler.h"
#include "Herwig/Shower/Dipole/Utility/IntrinsicPtGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.fh"
#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief The DipoleShowerHandler class manages the showering using
* the dipole shower algorithm.
*
* @see \ref DipoleShowerHandlerInterfaces "The interfaces"
* defined for DipoleShowerHandler.
*/
class DipoleShowerHandler: public ShowerHandler {
friend class Merger;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleShowerHandler();
/**
* The destructor.
*/
virtual ~DipoleShowerHandler();
//@}
public:
/**
* Indicate a problem in the shower.
*/
struct RedoShower {};
/**
* Insert an additional splitting kernel.
*/
void addSplitting(Ptr<DipoleSplittingKernel>::ptr sp) {
kernels.push_back(sp);
}
/**
* Reset the alpha_s for all splitting kernels.
*/
void resetAlphaS(Ptr<AlphaSBase>::tptr);
/**
* Reset the splitting reweight for all splitting kernels.
*/
void resetReweight(Ptr<DipoleSplittingReweight>::tptr);
/**
* Return true, if the shower handler can generate a truncated
* shower for POWHEG style events generated using Matchbox
*/
virtual bool canHandleMatchboxTrunc() const { return false; }
/**
* Return true, if this cascade handler will perform reshuffling from hard
* process masses.
*/
virtual bool isReshuffling() const { return false; }
/**
* Return the relevant hard scale to be used in the profile scales
*/
virtual Energy hardScale() const {
return muPt;
}
/**
* Calculate the alpha_s value the shower uses for Q.
*/
double as(Energy Q)const{return theGlobalAlphaS->value(sqr(Q));}
+ /**
+ *
+ */
+
+ double Nf(Energy Q)const{return theGlobalAlphaS->Nf(sqr(Q));}
protected:
typedef multimap<DipoleIndex,Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap;
/**
* The main method which manages the showering of a subprocess.
*/
virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb) {
return cascade(sub,xcomb,ZERO,ZERO);
}
/**
* The main method which manages the showering of a subprocess.
*/
tPPair cascade(tSubProPtr sub, XCPtr xcomb,
Energy optHardPt, Energy optCutoff);
/**
* Implementation of the merging algorithm.
**/
double reweightCKKW(int, int);
/**
* Build splitting generators for the given
* dipole index.
*/
void getGenerators(const DipoleIndex&,
Ptr<DipoleSplittingReweight>::tptr rw =
Ptr<DipoleSplittingReweight>::tptr());
/**
* Setup the hard scales.
*/
void hardScales(Energy2 scale);
/**
* Return the evolution ordering
*/
Ptr<DipoleEvolutionOrdering>::tptr evolutionOrdering() const { return theEvolutionOrdering; }
/**
* Reshuffle to constituent mass shells
*/
void constituentReshuffle();
/**
* Access the generator map
*/
GeneratorMap& generators() { return theGenerators; }
/**
* Access the event record
*/
DipoleEventRecord& eventRecord() { return theEventRecord; }
/**
* Return the event record
*/
const DipoleEventRecord& eventRecord() const { return theEventRecord; }
/**
* Return the splitting kernels.
*/
const vector<Ptr<DipoleSplittingKernel>::ptr>& splittingKernels() const {
return kernels;
}
/**
* Realign the event such as to have the incoming partons along thre
* beam axes.
*/
bool realign();
protected:
/**
* Perform the cascade.
*/
void doCascade(unsigned int& emDone,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Set the number of emissions
**/
void setNEmissions(unsigned int n){nEmissions=n;}
/**
* Get the winning splitting for the
* given dipole and configuration.
*/
Energy getWinner(DipoleSplittingInfo& winner,
const Dipole& dip,
pair<bool,bool> conf,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Get the winning splitting for the
* given dipole and configuration.
*/
Energy getWinner(SubleadingSplittingInfo& winner,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Get the winning splitting for the
* given dipole and configuration.
*/
Energy getWinner(DipoleSplittingInfo& winner,
const DipoleIndex& index,
double emitterX, double spectatorX,
pair<bool,bool> conf,
tPPtr emitter, tPPtr spectator,
Energy startScale,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
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).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The splitting kernels to be used.
*/
vector<Ptr<DipoleSplittingKernel>::ptr> kernels;
/**
* The evolution ordering considered
*/
Ptr<DipoleEvolutionOrdering>::ptr theEvolutionOrdering;
/**
* The ConstituentReshuffler to be used
*/
Ptr<ConstituentReshuffler>::ptr constituentReshuffler;
/**
* The intrinsic pt generator to be used.
*/
Ptr<IntrinsicPtGenerator>::ptr intrinsicPtGenerator;
/**
* A global alpha_s to be used for all splitting kernels.
*/
Ptr<AlphaSBase>::ptr theGlobalAlphaS;
/**
* Apply chain ordering to events from matrix
* element corrections.
*/
bool chainOrderVetoScales;
/**
* Limit the number of emissions.
* Limit applied if > 0.
*/
unsigned int nEmissions;
/**
* Discard events which did not radiate.
*/
bool discardNoEmissions;
/**
* Perform the first MC@NLO emission only.
*/
bool firstMCatNLOEmission;
/**
* The realignment scheme
*/
int realignmentScheme;
private:
/**
* The verbosity level.
* 0 - print no info
* 1 - print diagnostic information on setting up
* splitting generators etc.
* 2 - print detailed event information for up to
* printEvent events.
* 3 - print dipole chains after each splitting.
*/
int verbosity;
/**
* See verbosity.
*/
int printEvent;
private:
/**
* The splitting generators indexed by the dipole
* indices they can work on.
*/
GeneratorMap theGenerators;
/**
* The evnt record used.
*/
DipoleEventRecord theEventRecord;
/**
* The number of shoer tries so far.
*/
unsigned int nTries;
/**
* Whether or not we did radiate anything
*/
bool didRadiate;
/**
* Whether or not we did realign the event
*/
bool didRealign;
private:
/**
* A freezing value for the renormalization scale
*/
Energy theRenormalizationScaleFreeze;
/**
* A freezing value for the factorization scale
*/
Energy theFactorizationScaleFreeze;
/**
* The matching subtraction, if appropriate
*/
Ptr<ShowerApproximation>::tptr theShowerApproximation;
/**
* True, if sampler should apply compensation
*/
bool theDoCompensate;
/**
* Return the number of accepted points after which the grid should
* be frozen
*/
unsigned long theFreezeGrid;
/**
* The detuning factor applied to the sampling overestimate kernel
*/
double theDetuning;
/**
* A pointer to the dipole event reweight object
*/
Ptr<DipoleEventReweight>::ptr theEventReweight;
/**
* A pointer to a global dipole splitting reweight
*/
Ptr<DipoleSplittingReweight>::ptr theSplittingReweight;
/**
* True if no warnings have been issued yet
*/
static bool firstWarn;
/**
* The shower starting scale for the last event encountered
*/
Energy maxPt;
/**
* The shower hard scale for the last event encountered
*/
Energy muPt;
Ptr<MergerBase>::ptr theMergingHelper;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<DipoleShowerHandler> initDipoleShowerHandler;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleShowerHandler & operator=(const DipoleShowerHandler &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleShowerHandler. */
template <>
struct BaseClassTrait<Herwig::DipoleShowerHandler,1> {
/** Typedef of the first base class of DipoleShowerHandler. */
typedef Herwig::ShowerHandler NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleShowerHandler class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleShowerHandler>
: public ClassTraitsBase<Herwig::DipoleShowerHandler> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleShowerHandler"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleShowerHandler is implemented. It may also include several, space-separated,
* libraries if the class DipoleShowerHandler depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleShowerHandler_H */
diff --git a/Shower/Dipole/Kernels/DipoleSplittingKernel.cc b/Shower/Dipole/Kernels/DipoleSplittingKernel.cc
--- a/Shower/Dipole/Kernels/DipoleSplittingKernel.cc
+++ b/Shower/Dipole/Kernels/DipoleSplittingKernel.cc
@@ -1,395 +1,395 @@
// -*- C++ -*-
//
// DipoleSplittingKernel.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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 DipoleSplittingKernel class.
//
#include "DipoleSplittingKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/ShowerHandler.h"
using namespace Herwig;
DipoleSplittingKernel::DipoleSplittingKernel()
: HandlerBase(), theScreeningScale(0.0*GeV),
thePresamplingPoints(2000), theMaxtry(100000),
theFreezeGrid(500000),
theDetuning(1.0),
theStrictLargeN(false),
theFactorizationScaleFactor(1.0),
theRenormalizationScaleFactor(1.0),
theRenormalizationScaleFreeze(1.*GeV),
theFactorizationScaleFreeze(1.*GeV),
theCMWScheme(false),
theVirtualitySplittingScale(false),
presampling(false) {}
DipoleSplittingKernel::~DipoleSplittingKernel() {}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleSplittingKernel::persistentOutput(PersistentOStream & os) const {
os << theAlphaS << ounit(theScreeningScale,GeV) << theSplittingKinematics << thePDFRatio
<< thePresamplingPoints << theMaxtry << theFreezeGrid << theDetuning
<< theFlavour << theMCCheck << theStrictLargeN
<< theFactorizationScaleFactor
<< theRenormalizationScaleFactor
<< ounit(theRenormalizationScaleFreeze,GeV)
<< ounit(theFactorizationScaleFreeze,GeV)
<< theVirtualitySplittingScale<<theCMWScheme;
}
void DipoleSplittingKernel::persistentInput(PersistentIStream & is, int) {
is >> theAlphaS >> iunit(theScreeningScale,GeV) >> theSplittingKinematics >> thePDFRatio
>> thePresamplingPoints >> theMaxtry >> theFreezeGrid >> theDetuning
>> theFlavour >> theMCCheck >> theStrictLargeN
>> theFactorizationScaleFactor
>> theRenormalizationScaleFactor
>> iunit(theRenormalizationScaleFreeze,GeV)
>> iunit(theFactorizationScaleFreeze,GeV)
>> theVirtualitySplittingScale>>theCMWScheme;
}
double DipoleSplittingKernel::alphaPDF(const DipoleSplittingInfo& split,
Energy optScale,
double rScaleFactor,
double fScaleFactor) const {
Energy pt = optScale == ZERO ? split.lastPt() : optScale;
Energy2 scale = ZERO;
if ( !virtualitySplittingScale() ) {
scale = sqr(pt) + sqr(theScreeningScale);
} else {
scale = sqr(splittingKinematics()->QFromPt(pt,split)) + sqr(theScreeningScale);
}
if(split.fixedScale()>0.*GeV){
scale=sqr(split.fixedScale());
}
Energy2 rScale = sqr(theRenormalizationScaleFactor*rScaleFactor)*scale;
rScale = rScale > sqr(renormalizationScaleFreeze()) ? rScale : sqr(renormalizationScaleFreeze());
Energy2 fScale = sqr(theFactorizationScaleFactor*fScaleFactor)*scale;
fScale = fScale > sqr(factorizationScaleFreeze()) ? fScale : sqr(factorizationScaleFreeze());
double alphas = 1.0;
double pdf = 1.0;
// check if we are potentially reweighting and cache evaluations
bool evaluatePDF = true;
bool evaluateAlphaS = true;
bool variations =
!ShowerHandler::currentHandler()->showerVariations().empty() &&
!presampling;
if ( variations ) {
/*
cerr << "--------------------------------------------------------------------------------\n"
<< flush;
cerr << "variations ... central scale/GeV = "
<< sqrt(scale/GeV2) << " r = "
<< rScaleFactor << " f = " << fScaleFactor << "\n"
<< flush;
*/
map<double,double>::const_iterator pit = thePDFCache.find(fScaleFactor);
evaluatePDF = (pit == thePDFCache.end());
if ( !evaluatePDF ) {
//cerr << "PDF is cached: ";
pdf = pit->second;
//cerr << pdf << "\n" << flush;
}
map<double,double>::const_iterator ait = theAlphaSCache.find(rScaleFactor);
evaluateAlphaS = (ait == theAlphaSCache.end());
if ( !evaluateAlphaS ) {
//cerr << "alphas is cached: ";
alphas = ait->second;
//cerr << alphas << "\n" << flush;
}
}
if ( evaluateAlphaS )
alphas = alphaS()->value(rScale);
if ( evaluatePDF ) {
if ( split.index().initialStateEmitter() ) {
assert(pdfRatio());
pdf *=
split.lastEmitterZ() *
(*pdfRatio())(split.index().emitterPDF(), fScale,
split.index().emitterData(),split.emitterData(),
split.emitterX(),split.lastEmitterZ());
}
if ( split.index().initialStateSpectator() ) {
assert(pdfRatio());
pdf *=
split.lastSpectatorZ() *
(*pdfRatio())(split.index().spectatorPDF(), fScale,
split.index().spectatorData(),split.spectatorData(),
split.spectatorX(),split.lastSpectatorZ());
}
}
if ( evaluatePDF && variations ) {
//cerr << "caching PDF = " << pdf << "\n" << flush;
thePDFCache[fScaleFactor] = pdf;
}
if ( evaluateAlphaS && variations ) {
//cerr << "caching alphas = " << alphas << "\n" << flush;
theAlphaSCache[rScaleFactor] = alphas;
}
double ret = pdf;
if(split.fixedScale()<0.*GeV){
- ret *= alphas / (2.*Constants::pi)*(theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*5.)*alphaS()->value(rScale)/2./Constants::pi)):1.);
+ ret *= alphas / (2.*Constants::pi)*(theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*alphaS()->Nf(rScale))*alphaS()->value(rScale)/2./Constants::pi)):1.);
}else{
ret *=1.;
}
if ( ret < 0. )
ret = 0.;
return ret;
}
void DipoleSplittingKernel::accept(const DipoleSplittingInfo& split,
double, double,
map<string,double>& weights) const {
if ( ShowerHandler::currentHandler()->showerVariations().empty() )
return;
double reference = alphaPDF(split);
assert(reference > 0.);
for ( map<string,ShowerVariation>::const_iterator var =
ShowerHandler::currentHandler()->showerVariations().begin();
var != ShowerHandler::currentHandler()->showerVariations().end(); ++var ) {
if ( ( ShowerHandler::currentHandler()->firstInteraction() && var->second.firstInteraction ) ||
( !ShowerHandler::currentHandler()->firstInteraction() && var->second.secondaryInteractions ) ) {
/*
cerr << "reweighting in accept for: "
<< var->first
<< " using "
<< var->second.renormalizationScaleFactor << " "
<< var->second.factorizationScaleFactor << " "
<< var->second.firstInteraction << " "
<< var->second.secondaryInteractions << "\n" << flush;
*/
double varied = alphaPDF(split,ZERO,
var->second.renormalizationScaleFactor,
var->second.factorizationScaleFactor);
if ( varied != reference ) {
map<string,double>::iterator wi = weights.find(var->first);
if ( wi != weights.end() )
wi->second *= varied/reference;
else
weights[var->first] = varied/reference;
}
}
}
}
void DipoleSplittingKernel::veto(const DipoleSplittingInfo& split,
double p, double r,
map<string,double>& weights) const {
if ( ShowerHandler::currentHandler()->showerVariations().empty() )
return;
double reference = alphaPDF(split);
// this is dangerous, but we have no other choice currently -- need to
// carefully check for the effects; the assumption is that if the central
// one ius zero, then so will be the variations.
if ( reference == 0.0 )
return;
for ( map<string,ShowerVariation>::const_iterator var =
ShowerHandler::currentHandler()->showerVariations().begin();
var != ShowerHandler::currentHandler()->showerVariations().end(); ++var ) {
if ( ( ShowerHandler::currentHandler()->firstInteraction() && var->second.firstInteraction ) ||
( !ShowerHandler::currentHandler()->firstInteraction() && var->second.secondaryInteractions ) ) {
/*
cerr << "reweighting in veto for: "
<< var->first
<< " using "
<< var->second.renormalizationScaleFactor << " "
<< var->second.factorizationScaleFactor << " "
<< var->second.firstInteraction << " "
<< var->second.secondaryInteractions << "\n" << flush;
*/
double varied = alphaPDF(split,ZERO,
var->second.renormalizationScaleFactor,
var->second.factorizationScaleFactor);
if ( varied != reference ) {
map<string,double>::iterator wi = weights.find(var->first);
if ( wi != weights.end() )
wi->second *= (r - varied*p/reference) / (r-p);
else
weights[var->first] = (r - varied*p/reference) / (r-p);
}
}
}
}
AbstractClassDescription<DipoleSplittingKernel> DipoleSplittingKernel::initDipoleSplittingKernel;
// Definition of the static class description member.
void DipoleSplittingKernel::Init() {
static ClassDocumentation<DipoleSplittingKernel> documentation
("DipoleSplittingKernel is the base class for all kernels "
"used within the dipole shower.");
static Reference<DipoleSplittingKernel,AlphaSBase> interfaceAlphaS
("AlphaS",
"The strong coupling to be used by this splitting kernel.",
&DipoleSplittingKernel::theAlphaS, false, false, true, true, false);
static Parameter<DipoleSplittingKernel,Energy> interfaceScreeningScale
("ScreeningScale",
"A colour screening scale",
&DipoleSplittingKernel::theScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Reference<DipoleSplittingKernel,DipoleSplittingKinematics> interfaceSplittingKinematics
("SplittingKinematics",
"The splitting kinematics to be used by this splitting kernel.",
&DipoleSplittingKernel::theSplittingKinematics, false, false, true, false, false);
static Reference<DipoleSplittingKernel,PDFRatio> interfacePDFRatio
("PDFRatio",
"Set the optional PDF ratio object to evaluate this kernel",
&DipoleSplittingKernel::thePDFRatio, false, false, true, true, false);
static Parameter<DipoleSplittingKernel,unsigned long> interfacePresamplingPoints
("PresamplingPoints",
"The number of points used to presample this kernel.",
&DipoleSplittingKernel::thePresamplingPoints, 2000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,unsigned long> interfaceMaxtry
("Maxtry",
"The maximum number of attempts to generate a splitting.",
&DipoleSplittingKernel::theMaxtry, 10000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,unsigned long> interfaceFreezeGrid
("FreezeGrid",
"",
&DipoleSplittingKernel::theFreezeGrid, 500000, 1, 0,
false, false, Interface::lowerlim);
static Reference<DipoleSplittingKernel,ParticleData> interfaceFlavour
("Flavour",
"Set the flavour to be produced if ambiguous.",
&DipoleSplittingKernel::theFlavour, false, false, true, true, false);
static Reference<DipoleSplittingKernel,DipoleMCCheck> interfaceMCCheck
("MCCheck",
"[debug option] MCCheck",
&DipoleSplittingKernel::theMCCheck, false, false, true, true, false);
interfaceMCCheck.rank(-1);
static Switch<DipoleSplittingKernel,bool> interfaceStrictLargeN
("StrictLargeN",
"Work in a strict large-N limit.",
&DipoleSplittingKernel::theStrictLargeN, false, false, false);
static SwitchOption interfaceStrictLargeNOn
(interfaceStrictLargeN,
"On",
"Replace C_F -> C_A/2 where present",
true);
static SwitchOption interfaceStrictLargeNOff
(interfaceStrictLargeN,
"Off",
"Keep C_F=4/3",
false);
interfaceStrictLargeN.rank(-2);
static Switch<DipoleSplittingKernel,bool> interfaceCMWScheme
("CMWScheme",
"Add the CMW Scheme related Kg expression to the splitting",
&DipoleSplittingKernel::theCMWScheme, false, false, false);
static SwitchOption interfaceCMWSchemeOn
(interfaceCMWScheme,"On","", true);
static SwitchOption interfaceCMWSchemeOff
(interfaceCMWScheme,"Off","",false);
static Parameter<DipoleSplittingKernel,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&DipoleSplittingKernel::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
interfaceFactorizationScaleFactor.rank(-2);
static Parameter<DipoleSplittingKernel,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&DipoleSplittingKernel::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
interfaceRenormalizationScaleFactor.rank(-2);
static Parameter<DipoleSplittingKernel,Energy> interfaceRenormalizationScaleFreeze
("RenormalizationScaleFreeze",
"The freezing scale for the renormalization scale.",
&DipoleSplittingKernel::theRenormalizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,Energy> interfaceFactorizationScaleFreeze
("FactorizationScaleFreeze",
"The freezing scale for the factorization scale.",
&DipoleSplittingKernel::theFactorizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Switch<DipoleSplittingKernel,bool> interfaceVirtualitySplittingScale
("VirtualitySplittingScale",
"Use the virtuality as the splitting scale.",
&DipoleSplittingKernel::theVirtualitySplittingScale, false, false, false);
static SwitchOption interfaceVirtualitySplittingScaleYes
(interfaceVirtualitySplittingScale,
"Yes",
"Use vrituality.",
true);
static SwitchOption interfaceVirtualitySplittingScaleNo
(interfaceVirtualitySplittingScale,
"No",
"Use transverse momentum.",
false);
static Parameter<DipoleSplittingKernel,double> interfaceDetuning
("Detuning",
"A value to detune the overestimate kernel.",
&DipoleSplittingKernel::theDetuning, 1.0, 1.0, 0,
false, false, Interface::lowerlim);
}
diff --git a/Shower/Dipole/Merging/Merger.cc b/Shower/Dipole/Merging/Merger.cc
--- a/Shower/Dipole/Merging/Merger.cc
+++ b/Shower/Dipole/Merging/Merger.cc
@@ -1,1516 +1,1552 @@
// -*- C++ -*-
//
// Merger.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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 Merger class.
//
#include "Merger.h"
#include "Node.h"
#include "MFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Config/Constants.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "fastjet/ClusterSequence.hh"
using namespace Herwig;
Merger::Merger()
: MergerBase() {
- theNf=5;
+
minusL=false;
Unlopsweights=true;
theCMWScheme=true;
theSmearing=0.;
theDipoleShowerHandler=
Ptr<DipoleShowerHandler>::ptr();
+ FFLTK=new_ptr(FFLightTildeKinematics());
+ FILTK=new_ptr(FILightTildeKinematics());
+ IFLTK=new_ptr(IFLightTildeKinematics());
+ IILTK=new_ptr(IILightTildeKinematics());
+ FFMTK=new_ptr(FFMassiveTildeKinematics());
+ FIMTK=new_ptr(FIMassiveTildeKinematics());
+ IFMTK=new_ptr(IFMassiveTildeKinematics());
+
isUnitarized=true;
isNLOUnitarized=true;
defMERegionByJetAlg=false;
theOpenInitialStateZ=false;
theChooseHistory=8;
xiRenME=1.;
xiFacME=1.;
xiRenSh=1.;
xiFacSh=1.;
xiQSh=1.;
theGamma=1.;
ee_ycut=1000000;
pp_dcut=1000000;
therenormscale=0.*GeV;
theIRSafePT=1.*GeV;
theMergePt=3.*GeV;
theCentralMergePt=3.*GeV;
}
Merger::~Merger() {}
IBPtr Merger::clone() const {
return new_ptr(*this);
}
IBPtr Merger::fullclone() const {
return new_ptr(*this);
}
CrossSection Merger::MergingDSigDRBornGamma(NPtr Node){
double weightCL=0.;
weight=1.;
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), theMergePt))weight*= 0.;
}
NPtr Born= Node-> getHistory(true,xiQSh);
NPtr CLNode;
NPtr BornCL;
if( !Node->children().empty()){
if (UseRandom::rnd()<.5){
CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
weight*=inhist?(-2.*int(Node->children().size())):0.;
projected=true;Node->nodeME()->projectorStage(1);
weightCL=2.*int(Node->children().size());
BornCL=CLNode-> getHistory(false,xiQSh);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->legsize()-Node->nodeME()->projectorStage()))
return ZERO;
if(!Node->allAbove(mergePt()-0.1*GeV))weight=0.;
if(CLNode&&!CLNode->allAbove(mergePt()-0.1*GeV))weightCL=0.;
if (!Born->xcomb()->willPassCuts()){
if (Born->parent()) {
//cout<<"\n parent not pass";
}
return ZERO;
}
CrossSection res=ZERO;
- bool maxMulti=Node->legsize() == maxLegsLO();
+ bool maxMulti=Node->legsize() == int(maxLegsLO());
if(weight!=0.){
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node);
if (!fillProjector(projectedscale))weight=0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.&&weightCL==0.)return ZERO;
res+=weight*TreedSigDR(startscale,Node,(!maxMulti&&!projected)?theGamma:1.);
}
if(CLNode&&theGamma!=1.){
Energy startscale=CKKW_StartScale(BornCL);
Energy projectedscale=startscale;
fillHistory( startscale, BornCL, CLNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weightCL*=history.back().weight*alphaReweight()*pdfReweight();
CrossSection diff=ZERO;
Node->nodeME()->factory()->setAlphaParameter(1.);
diff-=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME));
Node->nodeME()->factory()->setAlphaParameter(theGamma);
string alp=(CLNode->dipol()->aboveAlpha()?"Above":"Below");
diff+=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME));
Node->nodeME()->factory()->setAlphaParameter(1.);
res+=diff;
}
return res;
}
CrossSection Merger::MergingDSigDRBornStandard(NPtr Node){
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), theMergePt))return ZERO;
}
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
return ZERO;
assert(Node->allAbove(mergePt()-0.1*GeV));
if (!Born->xcomb()->willPassCuts()){
return ZERO;
}
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
double w1=history.back().weight;
double w2=alphaReweight();
double w3=pdfReweight();
if(std::isnan(w1)){cout<<"\nhistory weight is nan";w1=0.;};
if(std::isnan(w2)){cout<<"\nalphaReweight weight is nan";w1=0.;};
if(std::isnan(w3)){cout<<"\npdfReweight weight is nan";w1=0.;};
weight*=w1*w2*w3;
if(weight==0.)return ZERO;
bool maxMulti=Node->xcomb()->meMomenta().size() == maxLegsLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection w4=TreedSigDR(startscale,Node,1.);
if(std::isnan(double(w4/nanobarn))){cout<<"\nTreedSigDR is nan";w1=0.;};
return weight*w4;
}
CrossSection Merger::MergingDSigDRVirtualStandard(NPtr Node ){
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node);
//if (history.size()==1&&Node->children().size()!=0) {
// cout<<"\n1-->"<<startscale/GeV<<" "<<weight;
//}
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=Node->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection matrixElement=LoopdSigDR(startscale,Node);
- assert(false);
CrossSection Bornweight=Node->nodeME()->dSigHatDRB();
CrossSection unlopsweight =(-sumpdfReweightUnlops()
-sumalphaReweightUnlops()
-sumfillHistoryUnlops())
*Bornweight
*SM().alphaS()/(2.*ThePEG::Constants::pi);
return weight*
as(startscale*xiRenSh)/SM().alphaS()*
(matrixElement+unlopsweight);
}
CrossSection Merger::MergingDSigDRRealStandard(NPtr Node){
bool allAbove=Node->allAbove(mergePt());
if(!Node->allAbove(theIRSafePT))return ZERO;
if (allAbove)return MergingDSigDRRealAllAbove(Node);
if (UseRandom::rnd()<.5)
return 2.*MergingDSigDRRealBelowSubReal( Node );
return 2.*MergingDSigDRRealBelowSubInt( Node);
}
CrossSection Merger::MergingDSigDRRealAllAbove(NPtr Node){
NPtr Born= Node-> getHistory(true,xiQSh);
NPtr CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
if (!CLNode->allAbove(mergePt()))return ZERO;
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection res= weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*
((inhist?TreedSigDR(startscale,Node):ZERO)
+CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME)));
// cout<<"\nCLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn "<<CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
return res;
}
CrossSection Merger::MergingDSigDRRealBelowSubReal(NPtr Node){
NPtrVec children=Node->children();
Selector<NPtr> HistNodeSel;
Energy minScale=generator()->maximumCMEnergy();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
if ((*child)->dipol()->lastPt()<minScale) {
minScale=(*child)->dipol()->lastPt();
HistNodeSel.insert(1.,*child);
}
}
NPtr HistNode=HistNodeSel.select(UseRandom::rnd());
NPtr Born=HistNode-> getHistory(false,xiQSh);
if( Born!= HistNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(1);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(0);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, HistNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=HistNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection sumPS=ZERO;
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
if ((*child)->allAbove(mergePt()))
sumPS-=(*child)->calcPs(startscale*xiFacME);
}
CrossSection me=TreedSigDR(startscale,Node);
//cout<<"\nSubReal "<<Node->miniPt()/GeV<<" "<<me<<" "<<sumPS<<" "<<me/sumPS;
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(me-sumPS);
}
CrossSection Merger::MergingDSigDRRealBelowSubInt(NPtr Node){
NPtr CLNode= Node->randomChild();
NPtr Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
if (!CLNode->allAbove(mergePt()))return ZERO;
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
pair<CrossSection,CrossSection> DipAndPs=
CLNode->calcDipandPS(startscale*xiFacME);
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*(DipAndPs.second-DipAndPs.first);
}
CrossSection Merger::TreedSigDR(Energy startscale,NPtr Node,double diffAlpha){
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
CrossSection res=DeepHead->nodeME()->dSigHatDRB();
if (diffAlpha!=1.) {
res+=DeepHead->nodeME()->dSigHatDRAlphaDiff(diffAlpha);
}
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
CrossSection Merger::LoopdSigDR(Energy startscale,NPtr Node){
// The deephead should be calculated here.
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
DeepHead->nodeME()->doOneLoopNoBorn();
CrossSection res=DeepHead->nodeME()->dSigHatDRV();
DeepHead->nodeME()->noOneLoopNoBorn();
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
bool Merger::fillProjector(Energy& prerunning){
// in the shower handler the scale is multiplied by xiQSh
// so here we need to devide this
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
if(history.begin()->node->deepHead()->nodeME()->projectorStage() == 0){
prerunning=(history.size()==1?1.:(1./xiQShfactor))*history.back().scale;
return true;
}
for (History::iterator it=history.begin();it!=history.end();it++){
if (projectorStage((*it).node)&&history.begin()->node->deepHead()->nodeME()->projectorStage() != 0){
history.begin()->node->deepHead()->xcomb()->lastProjector((*it).node->xcomb());
prerunning=(it==history.begin()?1.:(1./xiQShfactor))*(*it).scale;
return true;
}
}
return false;
}
double Merger::pdfReweight(){
double res=1.;
double facfactor=(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
for(int side=0;side!=2;side++){
if(history[0].node->xcomb()->mePartonData()[side]->coloured()){
History::iterator it=history.begin();
for (;it+1!=history.end();it++){
res *= pdfratio((*it).node, facfactor*(*it).scale,xiFacSh*((*it).node->dipol()->lastPt()), side);
facfactor=xiFacSh;
}
res*=pdfratio(history.back().node,facfactor*history.back().scale ,history[0].scale*xiFacME, side);
}
}
return res;
}
double Merger::alphaReweight(){
double res=1.;
Energy Q_R=(history[0].node->legsize()==N0()?xiRenME:xiRenSh)*history[0].scale;
res *= pow(as(Q_R) / SM().alphaS(), history[0].node->nodeME()->orderInAlphaS());
res *= pow(history[0].node->deepHead()->xcomb()->eventHandler().SM().alphaEMPtr()->value(history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED())/ SM().alphaEMMZ(), history[0].node->nodeME()->orderInAlphaEW());
if (!(history[0].node->children().empty())){
- res *=pow((theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*5.)*as(Q_R))/2./Constants::pi):1.),int(history[0].node->legsize()-N0()));
+ res *=pow((theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(Q_R))*as(Q_R))/2./Constants::pi):1.),int(history[0].node->legsize()-N0()));
}
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
if ((*it).node->parent()){
Energy q_i=xiRenSh* (*it).node->dipol()->lastPt();
res *= as(q_i)/ SM().alphaS()
- *(theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*5.)*as(q_i))/2./Constants::pi):1.);
+ *(theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(q_i))*as(q_i))/2./Constants::pi):1.);
}
}
return res;
}
void Merger::fillHistory(Energy scale, NPtr Begin, NPtr EndNode){
history.clear();
double sudakov0_n=1.;
history.push_back(HistoryStep(Begin,sudakov0_n,scale));
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
scale*=xiQShfactor;
if (Begin->parent()||!isUnitarized) {
while (Begin->parent() && (Begin != EndNode)) {
if (!dosudakov(Begin,scale, Begin->dipol()->lastPt(),sudakov0_n)){
history.push_back(HistoryStep(Begin->parent(),0.,scale));
}
scale=Begin->dipol()->lastPt();
history.push_back(HistoryStep(Begin->parent(),sudakov0_n,Begin->dipol()->lastPt()));
Begin = Begin->parent();
}
Energy notunirunning=scale;
if (!isUnitarized&&N()+N0() > int(Begin->deepHead()->legsize())) {
if (!dosudakov(Begin,notunirunning,mergePt(),sudakov0_n)){
history.back().weight=0.;
}else{
history.back().weight=sudakov0_n;
}
}
}
if( history.size()==1)scale/=xiQSh;
}
double Merger::sumpdfReweightUnlops(){
double res=0.;
Energy beam1Scale=history[0].scale*(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
Energy beam2Scale=history[0].scale*(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
History::iterator it=history.begin();
for (;it!=history.end();it++){
History::iterator ittmp=it;
ittmp++;
if((*it).node->xcomb()->mePartonData()[0]->coloured()&&(*it).node->nodeME()->lastX1()>0.99)return 0.;
if((*it).node->xcomb()->mePartonData()[1]->coloured()&&(*it).node->nodeME()->lastX2()>0.99)return 0.;
if (ittmp!=history.end()){
if((*it).node->nodeME()->lastX1()<0.00001)return 0.;
if((*it).node->nodeME()->lastX2()<0.00001)return 0.;
if ((*it).node->dipol()->bornEmitter() == 0 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
- theNf,
+ Nf(history[0].scale),
history[0].scale);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornEmitter() == 1 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
- theNf,
+ Nf(history[0].scale),
history[0].scale);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 0 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
- theNf,
+ Nf(history[0].scale),
history[0].scale);
//pdfratio((*it).node, beam1Scale, sqrt((*it).node->dipol()->lastPt()), 1);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 1 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
- theNf,
+ Nf(history[0].scale),
history[0].scale);
//pdfratio((*it).node, beam2Scale , sqrt((*it).node->dipol()->lastPt()), 2);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
}
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[0]->coloured()){
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().first->dataPtr(),
history.back().node->nodeME()->lastPartons().first->dataPtr(),
history.back().node->xcomb()->partonBins().first->pdf(),
beam1Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX1(),
- theNf,
+ Nf(history[0].scale),
history[0].scale);
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[1]->coloured()) {
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().second->dataPtr(),
history.back().node->nodeME()->lastPartons().second->dataPtr(),
history.back().node->xcomb()->partonBins().second->pdf(),
beam2Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX2(),
- theNf,
+ Nf(history[0].scale),
history[0].scale);
//history[0].node->deepHead()->nodeME()->pdf2(sqr(beam2Scale))/history[0].node->deepHead()->nodeME()->pdf2(sqr(history[0].scale));
}
return res;
}
double Merger::pdfUnlops(tcPDPtr particle,tcPDPtr parton,tcPDFPtr pdf,Energy running, Energy next,double x,int nlp,Energy fixedScale) {
//copied from PKOperator
double res=0.;
int number=10;//*int((max(running/next,next/running)));
for (int nr=0;nr<number;nr++){
double restmp=0;
using namespace RandomHelpers;
double r = UseRandom::rnd();
double eps = 1e-3;
pair<double,double> zw =
generate((piecewise(),
flat(0.0,x),
match(inverse(0.0,x,1.0) + inverse(1.0+eps,x,1.0))),r);
double z = zw.first;
double mapz = zw.second;
double PDFxparton=pdf->xfx(particle,parton,sqr(fixedScale),x)/x;
double CA = SM().Nc();
double CF = (SM().Nc()*SM().Nc()-1.0)/(2.*SM().Nc());
if (abs(parton->id()) < 7) {
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
double PDFxByzparton=pdf->xfx(particle,parton,sqr(fixedScale),x/z)*z/x;
assert(abs(parton->id()) < 7);
restmp += CF*(3./2.+2.*log(1.-x)) * PDFxparton;
if ( z > x ) {
restmp+= 0.5 * ( sqr(z) + sqr(1.-z) ) * PDFxByzgluon / z;
restmp += CF*2.*(PDFxByzparton - z*PDFxparton)/(z*(1.-z));
restmp -= CF*PDFxByzparton * (1.+z)/z;
}
}else{
assert(parton->id() == ParticleID::g);
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
if ( z > x ){
double factor = CF * ( 1. + sqr(1.-z) ) / sqr(z);
for ( int f = -nlp; f <= nlp; ++f ) {
if ( f == 0 )
continue;
restmp += pdf->xfx(particle,getParticleData(f),sqr(fixedScale),x/z)*z/x*factor;
}
}
- restmp += ( (11./6.) * CA - (1./3.)*theNf + 2.*CA*log(1.-x) ) *PDFxparton;
+ restmp += ( (11./6.) * CA - (1./3.)*Nf(history[0].scale) + 2.*CA*log(1.-x) ) *PDFxparton;
if ( z > x ) {
restmp += 2. * CA * ( PDFxByzgluon - z*PDFxparton ) / (z*(1.-z));
restmp += 2.* CA *( (1.-z)/z - 1. + z*(1.-z) ) * PDFxByzgluon / z;
}
}
if (PDFxparton<1e-8)restmp= 0.;
res+=1*restmp*log(sqr(running/next))/PDFxparton*mapz;
}
return res/number;
}
double Merger::sumalphaReweightUnlops(){
double res=0.;
if (!(history[0].node->children().empty())){
res +=alphasUnlops(history[0].scale*xiRenSh,
history[0].scale*xiRenME)*int(history[0].node->legsize()-N0());
assert(int(history[0].node->legsize()-N0()>0));
}
// dsig is calculated with fixed alpha_s
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
assert((*it).node->parent());
res +=alphasUnlops((*it).node->dipol()->lastPt()*xiRenSh ,history[0].scale);
}
return res;
}
double Merger::sumfillHistoryUnlops(){
double res=0.;
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
for (History::iterator it = history.begin(); (it+1) != history.end();it++){
doUNLOPS((*it).node,(it == history.begin()?xiQShfactor:1.)*(*it).scale , (*it).node->dipol()->lastPt() , history[0].scale, res);
}
return res;
}
Ptr<MFactory>::ptr Merger::treefactory(){return theTreeFactory;}
void Merger::doinit(){
assert(DSH()->hardScaleIsMuF());
}
CrossSection Merger::MergingDSigDR() {
history.clear();
- theNf=5;//TODO
if (theFirstNodeMap.find(theCurrentME)==theFirstNodeMap.end()) {
cout<<"\nnot in map:"<<theFirstNodeMap.size();
}
NPtr Node = theFirstNodeMap[theCurrentME]; assert(Node);
xiRenME=theCurrentME->renormalizationScaleFactor();
xiFacME=theCurrentME->factorizationScaleFactor();
xiRenSh=DSH()->renormalizationScaleFactor();
xiFacSh=DSH()->factorizationScaleFactor();
xiQSh=DSH()->hardScaleFactor();
DSH()->setCurrentHandler();
DSH()->eventHandler(generator()->eventHandler());
assert(DSH()->hardScaleIsMuF());
CrossSection res=ZERO;
if(Node->deepHead()->subtractedReal()){
res=MergingDSigDRRealStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(Node->deepHead()->virtualContribution()){
res=MergingDSigDRVirtualStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(theGamma!=1.){
res=MergingDSigDRBornGamma(Node);
theCurrentMaxLegs=maxLegsLO();
}else{
res=MergingDSigDRBornStandard(Node);
theCurrentMaxLegs=maxLegsLO();
}
//cout<<"\nrunning "<<Node->runningPt()/GeV<<" "<<Node->legsize();
theCurrentME->lastXCombPtr()->lastCentralScale(sqr(Node->runningPt()));
theCurrentME->lastXCombPtr()->lastShowerScale(sqr(Node->runningPt()));
if(theCurrentME->lastXCombPtr()->lastProjector()){
theCurrentME->lastXCombPtr()->lastProjector()->lastCentralScale(sqr(Node->runningPt()));
theCurrentME->lastXCombPtr()->lastProjector()->lastShowerScale(sqr(Node->runningPt()));
}
renormscale(0.0*GeV);
if (res == ZERO){
history.clear();
return res;
}
cleanup(Node);
DSH()->eventRecord().clear();
theCurrentME->lastXCombPtr()->subProcess(SubProPtr());
history.clear();
if(std::isnan(double(res/nanobarn))|| !std::isfinite(double(res/nanobarn))){
cout<<"Warning merger weight is "<<res/nanobarn<<" -> setting to 0";
return ZERO;
}
return res;
}
//----------------------------------Reviewed--------------------------------------//
void Merger::CKKW_PrepareSudakov(NPtr Born,Energy running){
//cleanup(Born);
tSubProPtr sub = Born->xcomb()->construct();
DSH()->resetPDFs(make_pair(Born->xcomb()->partonBins().first->pdf(),
Born->xcomb()->partonBins().second->pdf()),
Born->xcomb()->partonBins());
DSH()->setCurrentHandler();
DSH()->currentHandler()->generator()->currentEventHandler(Born->deepHead()->xcomb()->eventHandlerPtr());
DSH()->currentHandler()->remnantDecayer()->setHadronContent(Born->deepHead()->xcomb()->lastParticles());
DSH()->eventRecord().clear();
DSH()->eventRecord().prepare(sub, dynamic_ptr_cast<tStdXCombPtr>(Born->xcomb()), DSH()->pdfs(), Born->deepHead()->xcomb()->lastParticles());
DSH()->hardScales(sqr(running));
}
Energy Merger::CKKW_StartScale(NPtr Born){
Energy res=generator()->maximumCMEnergy();;
if(!Born->children().empty()){
for (int i=0;i<Born->legsize();i++){
if (!Born->nodeME()->mePartonData()[i]->coloured())continue;
for (int j=i+1;j<Born->legsize();j++){
if (i==j||!Born->nodeME()->mePartonData()[j]->coloured())continue;
res= min(res,sqrt(2.*Born->nodeME()->lastMEMomenta()[i]*Born->nodeME()->lastMEMomenta()[j]));
}
}
}else{
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res= sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale());
}
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res=max(res, sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale()));
return res;
}
double Merger::alphasUnlops( Energy next,Energy fixedScale) {
- double betaZero = (11./6.)*SM().Nc() - (1./3.)*theNf;
+ double betaZero = (11./6.)*SM().Nc() - (1./3.)*Nf(history[0].scale);
return (betaZero*log(sqr(fixedScale/next)))+
- (theCMWScheme?(((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*theNf))):0.);
+ (theCMWScheme?(((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(history[0].scale)))):0.);
}
double Merger::pdfratio(NPtr Born,Energy nominator_scale, Energy denominator_scale,int side){
StdXCombPtr bXc = Born->xcomb();
assert (bXc->mePartonData()[side]->coloured());
double from,to;
if (side==0){
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf1(sqr(denominator_scale));
to =Born->nodeME()->pdf1(sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
else{
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf2(sqr(denominator_scale));
to=Born->nodeME()->pdf2( sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
return to/from;
}
bool Merger::dosudakov(NPtr Born,Energy running, Energy next, double& sudakov0_n) {
CKKW_PrepareSudakov(Born, running);
for ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
sudakov0_n*=singlesudakov( dip, next,running,make_pair(true,false) );
sudakov0_n*=singlesudakov( dip, next,running,make_pair(false,true) );
if (sudakov0_n==0.0){
cleanup(Born);
return false;
}
}
}
cleanup(Born);
return true;
}
bool Merger::doUNLOPS(NPtr Born,Energy running, Energy next,Energy fixedScale, double& UNLOPS) {
CKKW_PrepareSudakov(Born, running);
for ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(true,false) );;
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(false,true) );
}
}
cleanup(Born);
return true;
}
bool Merger::projectorStage(NPtr Born){
return (Born->deepHead()->nodeME()->projectorStage() ==
int((Born->deepHead()->legsize()
- Born->legsize())));
}
void Merger::cleanup(NPtr Born) {
DSH()->eventRecord().clear();
if(!Born->xcomb()->subProcess())return;
ParticleVector vecfirst = Born->xcomb()->subProcess()->incoming().first->children();
for ( ParticleVector::const_iterator it = vecfirst.begin() ; it != vecfirst.end() ; it++ )
Born->xcomb()->subProcess()->incoming().first->abandonChild(*it);
ParticleVector vecsecond = Born->xcomb()->subProcess()->incoming().second->children();
for ( ParticleVector::const_iterator it = vecsecond.begin() ; it != vecsecond.end() ; it++ )
Born->xcomb()->subProcess()->incoming().second->abandonChild(*it);
Born->xcomb()->subProcess(SubProPtr());
}
double Merger::singlesudakov(list<Dipole>::iterator dip ,Energy next,Energy running,pair<bool,bool> conf ){
double res=1.;
tPPtr emitter = dip->emitter(conf);
tPPtr spectator = dip->spectator(conf);
DipoleSplittingInfo candidate((*dip).index(conf),conf,(*dip).emitterX(conf),(*dip).spectatorX(conf),emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=(*gen).second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res*=gen->second->sudakov(candidate,next);
}
}
return res;
}
double Merger::singleUNLOPS(list<Dipole>::iterator dip ,Energy next,Energy running,Energy fixedScale,pair<bool,bool> conf ){
double res=0.;
tPPtr emitter = dip->emitter(conf);
tPPtr spectator = dip->spectator(conf);
DipoleSplittingInfo candidate((*dip).index(conf),conf,(*dip).emitterX(conf),(*dip).spectatorX(conf),emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=gen->second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res+=gen->second->unlops(candidate,next,fixedScale);
}
}
return res;
}
void Merger::firstNodeMap(Ptr<MatchboxMEBase>::ptr a,NPtr b){theFirstNodeMap.insert(make_pair(a,b));}
map<Ptr<MatchboxMEBase>::ptr,NPtr> Merger::firstNodeMap(){return theFirstNodeMap;}
void Merger::setXComb(Ptr<MatchboxMEBase>::ptr me,tStdXCombPtr xc,int st){
theFirstNodeMap[me]->setXComb(xc, st);
}
void Merger::setKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->setKinematics();
}
void Merger::clearKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->clearKinematics();
}
bool Merger::generateKinematics(Ptr<MatchboxMEBase>::ptr me,const double * r){
theFirstNodeMap[me]->firstgenerateKinematics(r, 0,me->lastXCombPtr()->lastSHat());
if (theFirstNodeMap[me]->cutStage()==0 ){
bool inAlphaPS=false;
NPtrVec children=theFirstNodeMap[me]->children();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
treefactory()->setAlphaParameter(theGamma);
- inAlphaPS|=(*child)->dipol()->aboveAlpha();
+ inAlphaPS|=theGamma!=1.?(*child)->dipol()->aboveAlpha():false;
treefactory()->setAlphaParameter(1.);
}
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first&&!inAlphaPS)return false;
}
}
if (theFirstNodeMap[me]->cutStage()==1 ){
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first && !it->second.second)return false;
}
}
return true;
}
bool Merger::calculateInNode() const{
return theCalculateInNode;
}
void Merger::fillProjectors(Ptr<MatchboxMEBase>::ptr me){
for (unsigned int i = 0; i < (theFirstNodeMap[me]->Projector()).size(); ++i) {
me->lastXCombPtr()->projectors().insert(
(theFirstNodeMap[me]->Projector())[i].first,
(theFirstNodeMap[me]->Projector())[i].second->xcomb());
}
}
pair<bool,bool> Merger::clusterSafe(Ptr<MatchboxMEBase>::ptr me,int emit,int emis,int spec){
return theFirstNodeMap[me]->clusterSafe().find(make_pair(make_pair(emit,emis),spec))->second;
}
bool Merger::matrixElementRegion(PVector particles,Energy winnerScale,Energy cutscale){
//cout<<"\nparticles s"<<particles.size()<<" "<<particles[0]<<" "<<particles[1]<<flush;
/*
if (defMERegionByJetAlg && !particles[0]->coloured()&& !particles[1]->coloured()) {
assert(false);
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
input_particles.push_back(fastjet::PseudoJet(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[em]->momentum().e()/GeV));
}
fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets_ycut(ee_ycut);
return n==exclusive_jets.size();
}
if (defMERegionByJetAlg ) {
assert(false);
size_t noncol=0;
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
if (particles[em]->coloured())
input_particles.push_back(fastjet::PseudoJet(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[em]->momentum().e()/GeV));
else
noncol++;
}
double Rparam = 1.0;
fastjet::Strategy strategy = fastjet::Best;
fastjet::RecombinationScheme recomb_scheme = fastjet::E_scheme;
fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, recomb_scheme, strategy);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2-noncol;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets(pp_dcut);
// cout<<"\nn= "<<n<<" jets= "<<exclusive_jets.size();
//for (size_t i=0; i<exclusive_jets.size(); i++) {
//cout<<"\nn= "<<n<<" jetspt= "<<exclusive_jets[i].perp();
//}
return n==exclusive_jets.size();
}
*/
//assert(false);
Energy ptx=1000000.*GeV;
bool foundwinnerpt=false;
//FF
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (em==spe||emm==spe) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
- Energy scale = sqrt(2.*(emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom));
- double y = emissionmom*emittermom / (emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom);
- double z = emittermom*spectatormom / (emittermom*spectatormom + emissionmom*spectatormom);
- if (abs(scale * sqrt(y*z*(1.-z))-winnerScale)<0.0001*GeV) {
+ Energy pt=0*GeV;
+ if (emittermom.m()==0*GeV&&emissionmom.m()==0*GeV&&spectatormom.m()==0*GeV) {
+ pt=FFLTK->lastPt(emittermom,emissionmom,spectatormom);
+ }else{
+ pt=FFMTK->lastPt(emittermom,emissionmom,spectatormom);
+ }
+
+
+ if (abs(pt-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
// if(scale * sqrt(y*z*(1.-z))<optVeto&&winnerScale>optVeto)cout<<"\nFF "<<(scale * sqrt(y*z*(1.-z))/GeV);
- ptx =min(ptx,scale * sqrt(y*z*(1.-z)));
+ ptx =min(ptx,pt);
}
}
}
//FI
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
- Lorentz5Momentum Emsp=emittermom+emissionmom-spectatormom;
- Energy scale = sqrt(-1.*Emsp*Emsp);
- double x = (- emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom) /(emittermom*spectatormom + emissionmom*spectatormom);
- double z = emittermom*spectatormom / (emittermom*spectatormom + emissionmom*spectatormom);
- // if(scale * sqrt(z*(1.-z)*(1.-x)/x)<optVeto&&winnerScale>optVeto)cout<<"\nFI "<<(scale * sqrt(z*(1.-z)*(1.-x)/x)/GeV);
+ Energy pt=0*GeV;
+ if (emittermom.m()==0*GeV&&emissionmom.m()==0*GeV&&spectatormom.m()==0*GeV) {
+ pt=FILTK->lastPt(emittermom,emissionmom,spectatormom);
+ }else{
+ pt=FIMTK->lastPt(emittermom,emissionmom,spectatormom);
+ }
- if (abs(scale *sqrt(z*(1.-z)*(1.-x)/x)-winnerScale)<0.0001*GeV) {
+
+ if (abs(pt-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
- if((z*(1.-z)*(1.-x)/x)>0.)
- ptx =min(ptx,scale * sqrt(z*(1.-z)*(1.-x)/x));
+ if(pt>0.*GeV)
+ ptx =min(ptx,pt);
}
}
}
//IF
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (emm==spe) continue;
if (!(particles[em]->id()>6|| particles[em]->id()==particles[emm]->id() ||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
- Lorentz5Momentum Emsp=-emittermom+emissionmom+spectatormom;
- Energy scale = sqrt(-1.*Emsp*Emsp);
- double x = (-emissionmom*spectatormom + emittermom*spectatormom + emittermom*emissionmom) /(emittermom*emissionmom + emittermom*spectatormom);
- double u = emittermom*emissionmom / (emittermom*emissionmom + emittermom*spectatormom);
- if (abs(scale *sqrt(u*(1.-u)*(1.-x)/x)-winnerScale)<0.0001*GeV) {
+
+ Energy pt=0*GeV;
+ if (emittermom.m()==0*GeV&&emissionmom.m()==0*GeV&&spectatormom.m()==0*GeV) {
+ pt=IFLTK->lastPt(emittermom,emissionmom,spectatormom);
+ }else{
+ pt=IFMTK->lastPt(emittermom,emissionmom,spectatormom);
+ }
+
+
+ if (abs(pt-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
- if((u*(1.-u)*(1.-x)/x)>0.)
- ptx =min(ptx,scale * sqrt(u*(1.-u)*(1.-x)/x));
+ ptx =min(ptx,pt);
}
}
}
//II
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==spe) continue;
if (!(particles[em]->id()>6||
particles[em]->id()==particles[emm]->id() ||
particles[emm]->id()>6))continue;
//assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
- Lorentz5Momentum Emsp=emittermom-emissionmom+spectatormom;
- Energy scale = sqrt(Emsp*Emsp);
- double x = (emittermom*spectatormom - emittermom*emissionmom - spectatormom*emissionmom)/(emittermom*spectatormom);
- double v = (emittermom*emissionmom)/(emittermom*spectatormom);
- if (abs(scale * sqrt(v*(1.-x-v)/x)-winnerScale)<0.0001*GeV) {
+ Energy pt=IILTK->lastPt(emittermom,emissionmom,spectatormom);
+ if (abs(pt-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
- if ((v*(1.-x-v)/x)>0.)
- ptx =min(ptx, scale * sqrt(v*(1.-x-v)/x));
+ ptx =min(ptx, pt);
}
}
}
// cout<<"\n"<<cutscale/GeV<< " "<<foundwinnerpt<<" "<<ptx/GeV;
assert(foundwinnerpt);
return (ptx>cutscale) ;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Merger::persistentOutput(PersistentOStream & os) const {
os << minusL<< Unlopsweights<<
theCMWScheme<< projected<<
isUnitarized<< isNLOUnitarized<<
defMERegionByJetAlg<<theOpenInitialStateZ<<
theChooseHistory<<
theN0<< theOnlyN<<
theN<< theM<<
xiRenME<< xiFacME<<
xiRenSh<< xiFacSh<<
- xiQSh<< theNf<<
- weight<<weightCB<<theGamma<<ee_ycut<<pp_dcut<<theSmearing<<ounit(therenormscale, GeV)<<ounit(theIRSafePT, GeV)<<ounit(theMergePt, GeV)<<ounit(theCentralMergePt, GeV)<<theMergingJetFinder<<theLargeNBasis<<theDipoleShowerHandler<< theTreeFactory<<theFirstNodeMap ;
+ xiQSh<< weight<<weightCB<<theGamma<<ee_ycut<<pp_dcut<<theSmearing<<ounit(therenormscale, GeV)<<ounit(theIRSafePT, GeV)<<ounit(theMergePt, GeV)<<ounit(theCentralMergePt, GeV)<<theMergingJetFinder
+
+ <<theLargeNBasis
+
+
+ << FFLTK
+ << FILTK
+ << IFLTK
+ << IILTK
+ << FFMTK
+ << FIMTK
+ << IFMTK
+
+ <<theDipoleShowerHandler<< theTreeFactory<<theFirstNodeMap ;
}
void Merger::persistentInput(PersistentIStream & is, int) {
is >> minusL>> Unlopsweights>>
theCMWScheme>> projected>>
isUnitarized>> isNLOUnitarized>>
defMERegionByJetAlg>>theOpenInitialStateZ>>
theChooseHistory>>
theN0>> theOnlyN>>
theN>> theM>>
xiRenME>> xiFacME>>
xiRenSh>> xiFacSh>>
- xiQSh>> theNf>>
- weight>>weightCB>>theGamma>>ee_ycut>>pp_dcut>>theSmearing>>iunit(therenormscale, GeV)>>iunit(theIRSafePT, GeV)>>iunit(theMergePt, GeV)>>iunit(theCentralMergePt, GeV)>>theMergingJetFinder>>theLargeNBasis>>theDipoleShowerHandler>> theTreeFactory>>theFirstNodeMap ;
+ xiQSh>>
+ weight>>weightCB>>theGamma>>ee_ycut>>pp_dcut>>theSmearing>>iunit(therenormscale, GeV)>>iunit(theIRSafePT, GeV)>>iunit(theMergePt, GeV)>>iunit(theCentralMergePt, GeV)>>theMergingJetFinder>>theLargeNBasis>>
+
+
+ FFLTK
+ >> FILTK
+ >> IFLTK
+ >> IILTK
+ >> FFMTK
+ >> FIMTK
+ >> IFMTK>>
+
+ theDipoleShowerHandler>> theTreeFactory>>theFirstNodeMap ;
}
// *** 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<Merger, Herwig::MergerBase>
describeHerwigMerger("Herwig::Merger", "HwDipoleShower.so");
//ClassDescription<Merger> Merger::initMerger;
// Definition of the static class description member.
void Merger::Init() {
static ClassDocumentation<Merger> documentation
("Merger.");
static Reference<Merger,DipoleShowerHandler> interfaceShowerHandler
("DipoleShowerHandler",
"",
&Merger::theDipoleShowerHandler, false, false, true, true, false);
static Switch<Merger,bool>
interfaceminusL ("minusL","",&Merger::minusL, false, false, false);
static SwitchOption interfaceminusLYes
(interfaceminusL,"Yes","",true);
static SwitchOption interfaceminusLNo
(interfaceminusL,"No","",false);
static Switch<Merger,bool>
interfaceUnlopsweights ("Unlopsweights","",&Merger::Unlopsweights, false, false, false);
static SwitchOption interfaceUnlopsweightsYes
(interfaceUnlopsweights,"Yes","",true);
static SwitchOption interfaceUnlopsweightsNo
(interfaceUnlopsweights,"No","",false);
static Switch<Merger,bool>
interfacetheCMWScheme ("CMWScheme","",&Merger::theCMWScheme, false, false, false);
static SwitchOption interfacetheCMWSchemeYes
(interfacetheCMWScheme,"Yes","",true);
static SwitchOption interfacetheCMWSchemeNo
(interfacetheCMWScheme,"No","",false);
static Parameter<Merger,Energy> interfaceMergerScale
("MergingScale",
"The MergingScale.",
&Merger::theCentralMergePt, GeV, 20.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Reference<Merger,MFactory> interfaceMergerHelper
("MFactory",
"",
&Merger::theTreeFactory, false, false, true, true, false);
static Parameter<Merger,double> interfacedcut
("pp_dcut",
"The MergingScale.",
&Merger::pp_dcut, 5.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<Merger,double> interfaceycut
("ee_ycut",
"The MergingScale.",
&Merger::ee_ycut, 0.2, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<Merger,double> interfacegamma
("gamma",
"gamma parameter.",
&Merger::theGamma, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Reference<Merger,JetFinder> interfaceMergingJetFinder
("MergingJetFinder",
"A reference to the jet finder used in Merging.",
&Merger::theMergingJetFinder, false, false, true, false, false);
static Reference<Merger,ColourBasis> interfaceLargeNBasis
("LargeNBasis",
"Set the large-N colour basis implementation.",
&Merger::theLargeNBasis, false, false, true, true, false);
static Switch<Merger,bool>
interfacedefMERegionByJetAlg
("MERegionByJetAlg","",&Merger::defMERegionByJetAlg, false, false, false);
static SwitchOption interfacedefMERegionByJetAlgYes
(interfacedefMERegionByJetAlg,"Yes","",true);
static SwitchOption interfacedefMERegionByJetAlgNo
(interfacedefMERegionByJetAlg,"No","",false);
static Switch<Merger,bool>
interfaceOpenInitialSateZ
("OpenInitialStateZ","",&Merger::theOpenInitialStateZ, false, false, false);
static SwitchOption interfaceOpenInitialSateZYes
(interfaceOpenInitialSateZ,"Yes","",true);
static SwitchOption interfaceOpenInitialSateZNo
(interfaceOpenInitialSateZ,"No","",false);
static Parameter<Merger, Energy>
interfaceIRSafePT
("IRSafePT", "Set the pt for which a matrixelement should be treated as IR-safe.",
&Merger::theIRSafePT,
GeV, 0.0 * GeV, ZERO, Constants::MaxEnergy, true, false, Interface::limited);
interfaceIRSafePT.setHasDefault(false);
static Parameter<Merger, double> interfacemergePtsmearing("MergingScaleSmearing", "Set the percentage the merging pt should be smeared.",
&Merger::theSmearing, 0., 0.,
0.0, 0.5, true, false, Interface::limited);
static Parameter<Merger, int> interfacechooseHistory("chooseHistory", "different ways of choosing the history", &Merger::theChooseHistory, 3, -1, 0,
false, false, Interface::lowerlim);
static Parameter<Merger, int> interfaceaddLOLegs("addLOLegs", "Set the number of additional jets to consider.", &Merger::theN, 0, 0,
0, false, false, Interface::lowerlim);
static Parameter<Merger, int> interfaceaddNLOLegs("addNLOLegs",
"Set the number of virtual corrections to consider. -1 is default for no virtual correction.", &Merger::theM, -1, -1, 0, false, false,
Interface::lowerlim);
}
diff --git a/Shower/Dipole/Merging/Merger.h b/Shower/Dipole/Merging/Merger.h
--- a/Shower/Dipole/Merging/Merger.h
+++ b/Shower/Dipole/Merging/Merger.h
@@ -1,367 +1,388 @@
// -*- C++ -*-
//
// Merger.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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_Merger_H
#define HERWIG_Merger_H
//
// This is the declaration of the Merger class.
//
#include "MFactory.fh"
#include "Node.fh"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Vectors/SpinOneLorentzRotation.h"
#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
#include "Herwig/Shower/Dipole/Base/DipoleSplittingGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
+#include "Herwig/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h"
+#include "Herwig/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
+#include "Herwig/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h"
+#include "Herwig/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h"
+#include "Herwig/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
+#include "Herwig/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
+#include "Herwig/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h"
+
#include "ThePEG/Cuts/JetFinder.h"
#include "ThePEG/Cuts/Cuts.h"
namespace Herwig {
using namespace ThePEG;
typedef Ptr<Node>::ptr NPtr;
typedef vector<Ptr<Node>::ptr> NPtrVec;
struct HistoryStep {
HistoryStep(){}
HistoryStep(NPtr cn,double w ,Energy sc){
node=cn;
weight=w;
scale=sc;
}
NPtr node;
double weight;
Energy scale;
};
typedef vector< HistoryStep > History;
typedef multimap<DipoleIndex,Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap2;
/**
* \ingroup DipoleShower
* \author Johannes Bellm
*
* \brief Merger handles the Merger ....... //TODO .
*
* @see \ref MergerInterfaces "The interfaces"
* defined for Merger.
*/
class Merger: public MergerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
Merger();
/**
* The destructor.
*/
virtual ~Merger();
//@}
public:
double singlesudakov(list<Dipole>::iterator,Energy,Energy,pair<bool,bool>);
bool dosudakov(NPtr Born,Energy running, Energy next, double& sudakov0_n);
//bool dosudakovold(NPtr, Energy, Energy, double&);
//bool sudakov(NPtr Born, Energy running, Energy next);
void cleanup(NPtr);
bool projectorStage(NPtr);
Energy CKKW_StartScale(NPtr);
void CKKW_PrepareSudakov(NPtr,Energy);
CrossSection TreedSigDR(Energy startscale,NPtr,double diffalpha=1.);
CrossSection LoopdSigDR(Energy startscale,NPtr);
bool fillProjector(Energy&);
void fillHistory(Energy, NPtr, NPtr );
double sumpdfReweightUnlops();
double sumalphaReweightUnlops();
double pdfratio(NPtr,Energy,Energy,int);
double pdfReweight();
double alphaReweight();
double sumfillHistoryUnlops();
double alphasUnlops( Energy next,Energy fixedScale);
double pdfUnlops(tcPDPtr,tcPDPtr,tcPDFPtr,Energy,Energy,double,int,Energy);
double singleUNLOPS(list<Dipole>::iterator,Energy,Energy,Energy,pair<bool,bool>);
bool doUNLOPS(NPtr Born,Energy running, Energy next,Energy fixedScale, double& UNLOPS);
void setME(Ptr<MatchboxMEBase>::ptr me){theCurrentME=me;}
CrossSection MergingDSigDR() ;
CrossSection MergingDSigDRBornStandard(NPtr Node);
CrossSection MergingDSigDRBornGamma(NPtr Node);
CrossSection MergingDSigDRVirtualStandard(NPtr Node);
CrossSection MergingDSigDRRealStandard(NPtr Node);
CrossSection MergingDSigDRRealAllAbove(NPtr Node);
CrossSection MergingDSigDRRealBelowSubReal(NPtr Node);
CrossSection MergingDSigDRRealBelowSubInt(NPtr Node);
double as(Energy q){return DSH()->as(q);}
Ptr<DipoleShowerHandler>::ptr DSH(){return theDipoleShowerHandler;}
+ Ptr<DipoleShowerHandler>::ptr DSH()const{return theDipoleShowerHandler;}
size_t maxLegs() const {return theCurrentMaxLegs;}
size_t maxLegsLO() const {return N0()+N();}
- size_t maxLegsNLO()const {return N0()+M()+1;}
+ size_t maxLegsNLO()const {return N0()+M();}
+ double Nf(Energy scale)const{return DSH()->Nf(scale);}
+
Ptr<MFactory>::ptr treefactory();
map<Ptr<MatchboxMEBase>::ptr,NPtr> firstNodeMap();
void setXComb(Ptr<MatchboxMEBase>::ptr,tStdXCombPtr,int);
void setKinematics(Ptr<MatchboxMEBase>::ptr);
void clearKinematics(Ptr<MatchboxMEBase>::ptr);
bool generateKinematics(Ptr<MatchboxMEBase>::ptr,const double *);
bool calculateInNode() const;
void fillProjectors(Ptr<MatchboxMEBase>::ptr);
pair<bool,bool> clusterSafe(Ptr<MatchboxMEBase>::ptr,int,int,int);
bool theCalculateInNode;
bool matrixElementRegion(PVector particles,Energy winnerScale=0.*GeV,Energy cutscale=0.*GeV);
void renormscale(Energy x) {
therenormscale = x;
}
Energy renormscale() {
return therenormscale;
}
Energy mergingScale()const{
return theMergePt;
}
Energy mergePt() {
return theMergePt;
}
void mergePt(Energy x) {
theMergePt = x;
}
Energy centralMergePt() {
return theCentralMergePt;
}
void centralMergePt(Energy x) {
theCentralMergePt = x;
}
void firstNodeMap(Ptr<MatchboxMEBase>::ptr,NPtr);
void smeareMergePt(){theMergePt=centralMergePt()*(1.+0.*(-1.+2.*UseRandom::rnd())*smear());}
double gamma()const{return theGamma;}
double smear()const{return theSmearing;}
bool MERegionByJetAlg(){return defMERegionByJetAlg;}
Ptr<ColourBasis>::ptr largeNBasis(){return theLargeNBasis;}
void largeNBasis(Ptr<ColourBasis>::ptr x){theLargeNBasis=x;}
int M()const{return theM;}
int N()const{return theN;}
int N0()const{return theN0;}
void N0(int n){ theN0=n;}
bool openInitialStateZ(){return theOpenInitialStateZ;}
private:
/**
* Variations
*/
bool minusL;
bool Unlopsweights;
bool theCMWScheme;
bool projected;
bool isUnitarized;
bool isNLOUnitarized;
bool defMERegionByJetAlg;
bool theOpenInitialStateZ;
int theChooseHistory;
int theN0;
int theOnlyN;
int theN;
int theM;
size_t theCurrentMaxLegs;
double xiRenME;
double xiFacME;
double xiRenSh;
double xiFacSh;
double xiQSh;
double theNf;
double weight,weightCB;
double theGamma;
double ee_ycut;
double pp_dcut;
double theSmearing;
Energy therenormscale;
/**
* If any clustering below the CutStage has a lower pT than the MergePt
* the phase space point has to be thrown away.
*/
Energy theIRSafePT;
Energy theMergePt;
Energy theCentralMergePt;
History history;
Ptr<JetFinder>::ptr theMergingJetFinder;
Ptr<ColourBasis>::ptr theLargeNBasis;
Ptr<MatchboxMEBase>::ptr theCurrentME;
+ Ptr<FFLightTildeKinematics>::ptr FFLTK;
+ Ptr<FILightTildeKinematics>::ptr FILTK;
+ Ptr<IFLightTildeKinematics>::ptr IFLTK;
+ Ptr<IILightTildeKinematics>::ptr IILTK;
+ Ptr<FFMassiveTildeKinematics>::ptr FFMTK;
+ Ptr<FIMassiveTildeKinematics>::ptr FIMTK;
+ Ptr<IFMassiveTildeKinematics>::ptr IFMTK;
+
+
+
/**
* The mean of the Gaussian distribution for
* the intrinsic pt of sea partons.
*/
Ptr<DipoleShowerHandler>::ptr theDipoleShowerHandler;
Ptr<MFactory>::ptr theTreeFactory;
map<Ptr<MatchboxMEBase>::ptr,NPtr> theFirstNodeMap;
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
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.
*/
Merger & operator=(const Merger &);
};
}
#endif /* HERWIG_Merger_H */
diff --git a/Shower/Dipole/Merging/Node.cc b/Shower/Dipole/Merging/Node.cc
--- a/Shower/Dipole/Merging/Node.cc
+++ b/Shower/Dipole/Merging/Node.cc
@@ -1,660 +1,622 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Node class.
//
#include "Node.h"
#include "MFactory.h"
#include "Merger.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
Node::Node() {
}
bool NodeDebug=false;
Node::Node(Ptr<MatchboxMEBase>::ptr nodeME, int deepprostage, int cutstage) {
nodeME->maxMultCKKW(1);
nodeME->minMultCKKW(0);
theDeepHead = this;
thenodeMEPtr = nodeME;
theDeepProStage = deepprostage;
theCutStage = cutstage;
theSubtractedReal=false;
theVirtualContribution=false;
}
Node::Node(Ptr<Node>::ptr deephead, Ptr<Node>::ptr head, Ptr<SubtractionDipole>::ptr dipol, Ptr<MatchboxMEBase>::ptr nodeME,
int deepprostage, int cutstage) {
theDeepHead = deephead;
theparent = head;
thedipol = dipol;
thenodeMEPtr = nodeME;
theDeepProStage = deepprostage;
theCutStage = cutstage;
theSubtractedReal=false;
theVirtualContribution=false;
}
Node::~Node() { }
Ptr<SubtractionDipole>::ptr Node::dipol() const {
return thedipol;
}
/** returns the matrix element pointer */
const Ptr<MatchboxMEBase>::ptr Node::nodeME() const {
return thenodeMEPtr;
}
/** access the matrix element pointer */
Ptr<MatchboxMEBase>::ptr Node::nodeME() {
return thenodeMEPtr;
}
int Node::legsize() const {return nodeME()->legsize();}
Ptr<Node>::ptr Node::randomChild() {
return thechildren[(int)(UseRandom::rnd() * thechildren.size())];
}
Energy Node::miniPt() const{
Energy res=1000000000*GeV;
for (vector<Ptr<Node>::ptr>::const_iterator it = thechildren.begin(); it != thechildren.end(); it++) {
res=min(res,(*it)->dipol()->lastPt());
}
return res;
}
bool Node::allAbove(Energy pt){
for (vector<Ptr<Node>::ptr>::iterator it = thechildren.begin(); it != thechildren.end(); it++) {
if((*it)->dipol()->lastPt()<pt)return false;
}
return true;
}
bool Node::isInHistoryOf(Ptr<Node>::ptr other){
while (other->parent()) {
if(other==this)return true;
other=other->parent();
}
return false;
}
void Node::flushCaches() {
this->theProjectors.clear();
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->underlyingBornME()->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r = thechildren[i]->dipol()->reweights().begin() ; r != thechildren[i]->dipol()->reweights().end() ;
++r ) {
(**r).flushCaches();
}
if ( thechildren[i]->xcomb() ) thechildren[i]->xcomb()->clean();
thechildren[i]->nodeME()->flushCaches();
thechildren[i]->flushCaches();
}
}
void Node::setKinematics() {
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
thechildren[i]->dipol()->setKinematics();
thechildren[i]->nodeME()->setKinematics();
thechildren[i]->setKinematics();
}
}
void Node::clearKinematics() {
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
thechildren[i]->nodeME()->clearKinematics();
thechildren[i]->dipol()->clearKinematics();
thechildren[i]->clearKinematics();
}
}
bool Node::generateKinematics(const double *r, int stage, Energy2 shat) {
isOrdered = false;
isthissafe=true;
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
if ( !thechildren[i]->dipol()->generateKinematics(r) ) cout << "stop";
if ( dipol()->lastPt() < thechildren[i]->dipol()->lastPt() && parent()->ordered() ) {
isOrdered = true;
deepHead()->setOrderedSteps(stage + 1);
}
thechildren[i]->generateKinematics(r, stage + 1, thechildren[i]->xcomb()->lastSHat());
isthissafe = (isthissafe && thechildren[i]->dipol()->lastPt() >=deepHead()->MH()->mergePt());
}
return isthissafe;
}
void Node::firstgenerateKinematics(const double *r, int stage, Energy2 shat) {
flushCaches();
//Set here the new merge Pt for the next phase space point.( Smearing!!!)
MH()->smeareMergePt();
isOrdered = true;
// if we consider the hard scale to play a role in the steps we must change here to 0.
setOrderedSteps(1);
this->orderedSteps();
clustersafer.clear();
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
bool ifirst = true;
bool isecond = true;
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
if ( !thechildren[i]->dipol()->generateKinematics(r) ) cout << "stop";
isecond = thechildren[i]->generateKinematics(r, stage + 1, thechildren[i]->xcomb()->lastSHat());
ifirst = (thechildren[i]->dipol()->lastPt() >= deepHead()->MH()->mergePt());
pair<pair<int, int>, int> EmitEmisSpec = make_pair(make_pair(thechildren[i]->dipol()->realEmitter(), thechildren[i]->dipol()->realEmission()),
thechildren[i]->dipol()->realSpectator());
clustersafer.insert(make_pair(EmitEmisSpec, make_pair(ifirst, isecond)));
}
}
void Node::setXComb(tStdXCombPtr xc, int proStage) {
if ( !parent() ) this->xcomb(xc);
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
if ( !thechildren[i]->xcomb() ) {
thechildren[i]->xcomb(thechildren[i]->dipol()->makeBornXComb(xc));
assert(thechildren[i]->xcomb());
thechildren[i]->xcomb()->head(xc);
if ( !thechildren[i]->dipol()->lastXCombPtr() ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
}
thechildren[i]->setXComb(thechildren[i]->xcomb(), (proStage - 1));
} else {
if ( !(thechildren[i]->dipol()->lastXCombPtr()->lastScale() == thechildren[i]->xcomb()->lastScale()) ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
}
if ( thechildren[i]->xcomb()->head() != xc ) thechildren[i]->xcomb()->head(xc);
thechildren[i]->setXComb(thechildren[i]->xcomb(), (proStage - 1));
}
}
}
void Node::birth(vector<Ptr<MatchboxMEBase>::ptr> vec) {
vector<Ptr<SubtractionDipole>::ptr> dipoles =
nodeME()->getDipoles(DipoleRepository::dipoles(
nodeME()->factory()->dipoleSet()), vec,true);
for ( unsigned int j = 0 ; j < dipoles.size() ; ++j ) {
dipoles[j]->doSubtraction();
Ptr<Node>::ptr node = new_ptr(Node(theDeepHead,this, dipoles[j],
dipoles[j]->underlyingBornME(),
theDeepHead->DeepProStage(),
theDeepHead->cutStage()));
thechildren.push_back(node);
}
}
vector<Ptr<Node>::ptr> Node::getNextOrderedNodes(bool normal,double hardScaleFactor) {
vector<Ptr<Node>::ptr> temp = children();
vector<Ptr<Node>::ptr> res;
for ( vector<Ptr<Node>::ptr>::const_iterator it = temp.begin() ; it != temp.end() ; ++it ) {
if(deepHead()->MH()->mergePt()>(*it)->dipol()->lastPt()){
res.clear();
// if any of the nodes is below the merging scale return empty vector
return res;
continue;
}
if (parent()&& normal){
if ( (*it)->dipol()->lastPt() < dipol()->lastPt() ){
continue;
}
}
if ( (*it)->children().size() != 0 ){
vector<Ptr<Node>::ptr> tempdown = (*it)->children();
for ( vector<Ptr<Node>::ptr>::const_iterator itd = tempdown.begin() ; itd != tempdown.end() ; ++itd ) {
if( (*itd)->dipol()->lastPt() > (*it)->dipol()->lastPt()&&(*it)->inShowerPS((*itd)->dipol()->lastPt()) ){
res.push_back(*it);
break;
}
}
}
else {
(*it)->nodeME()->factory()->scaleChoice()->setXComb((*it)->xcomb());
if ( sqr(hardScaleFactor)*(*it)->nodeME()->factory()->scaleChoice()->renormalizationScale()
>= (*it)->dipol()->lastPt() * (*it)->dipol()->lastPt()&&
(*it)->inShowerPS(hardScaleFactor*sqrt((*it)->nodeME()->factory()->scaleChoice()->renormalizationScale()))) {
res.push_back(*it);
}
}
}
return res;
}
bool Node::inShowerPS(Energy hardpT){
assert(deepHead()->MH()->largeNBasis());
+ double z_=dipol()->lastZ();
+ // II
+ if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()<2&&deepHead()->MH()->openInitialStateZ()) return true;
+ // IF
+ if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()>=2&&deepHead()->MH()->openInitialStateZ())
+ return true;
- double x=0.;
- double z_=dipol()->lastZ();
+ pair<double,double> zbounds=
+ dipol()->tildeKinematics()->zBounds(dipol()->lastPt(),hardpT);
- // if (dipol()->lastPt()>50*GeV) {
- // return false;
- // }
-
- string type;
- // II
- if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()<2){
- type="II";
- x =dipol()->bornEmitter()==0?xcomb()->lastX1():xcomb()->lastX2();
- double ratio = sqr(dipol()->lastPt()/dipol()->lastDipoleScale());
- double x__ = z_*(1.-z_)/(1.-z_+ratio);
- double v = ratio*z_ /(1.-z_+ratio);
- if (dipol()->lastPt()>(1.-x) * dipol()->lastDipoleScale()/ (2.*sqrt(x)))return false;
- assert(v< 1.-x__&&x > 0. && x < 1. && v > 0.);
- if (deepHead()->MH()->openInitialStateZ()) {
- return true;
- }
- }
- // IF
- if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()>=2){
- type="IF";
- x =dipol()->bornEmitter()==0?xcomb()->lastX1():xcomb()->lastX2();
- if (dipol()->lastPt()>dipol()->lastDipoleScale()* sqrt((1.- x)/x) /2.)return false;
- if (deepHead()->MH()->openInitialStateZ()) {
- return true;
- }
- }
- // FI
- if( dipol()->bornEmitter()>=2&&dipol()->bornSpectator()<2){
- type="FI";
- double lastx=dipol()->bornSpectator()==0?xcomb()->lastX1():xcomb()->lastX2();
- if (dipol()->lastPt()>dipol()->lastDipoleScale()*sqrt((1.-lastx)/lastx) /2.)return false;
- }
- // FF
- if( dipol()->bornEmitter()>=2&&dipol()->bornSpectator()>=2){
- type="FF";
- if (dipol()->lastPt()>dipol()->lastDipoleScale()/2.)return false;
- }
-
- double kappa=sqr(dipol()->lastPt()/hardpT);
-
- double s = sqrt(1.-kappa);
- pair<double,double> zbounds= make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
-
return (zbounds.first<z_&&z_<zbounds.second);
}
Ptr<Node>::ptr Node::getHistory(bool normal,double hardScaleFactor) {
Ptr<Node>::ptr res=this;
vector<Ptr<Node>::ptr> temp = getNextOrderedNodes(normal,hardScaleFactor);
Energy minpt=100000.*GeV;
Selector<Ptr<Node>::ptr> subprosel;
while (temp.size()!=0){
minpt=100000.*GeV;
subprosel.clear();
for (vector<Ptr<Node>::ptr>::iterator it=temp.begin();it!=temp.end();it++){
if( (*it)->dipol()->underlyingBornME()->largeNColourCorrelatedME2(
make_pair((*it)->dipol()->bornEmitter(),(*it)->dipol()->bornSpectator()),
deepHead()->MH()->largeNBasis())!=0.
){
double weight=abs((*it)->dipol()->dSigHatDR()/nanobarn);
if(weight!=0.){
subprosel.insert(weight , (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
/*
if((*it)->nodeME()->dSigHatDR()/nanobarn!=0.){
subprosel.insert((abs((*it)->dipol()->dSigHatDR() /
(*it)->nodeME()->dSigHatDR()*deepHead()->MH()->as((*it)->dipol()->lastPt()))), (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
*/
//TODO choosehistories
}
}
if (subprosel.empty())
return res;
res = subprosel.select(UseRandom::rnd());
temp = res->getNextOrderedNodes(true,hardScaleFactor);
}
return res;
}
bool Node::headContribution(double hardScaleFactor){
bool allabove=true;
vector<Ptr<Node>::ptr> temp2 = children();
for (vector<Ptr<Node>::ptr>::iterator it = temp2.begin(); it != temp2.end(); it++) {
allabove&=(*it)->dipol()->lastPt()>deepHead()->MH()->mergePt();
}
if(allabove){
Ptr<Node>::ptr tmpBorn = getHistory(true,hardScaleFactor);
if(!tmpBorn->parent())return false;
}
return true;
}
bool Node::DipolesAboveMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
double Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
if ((Di!=0.)&&(*it)->xcomb()->willPassCuts()) {
assert((*it)->dipol()->clustersafe());
minpt=min(minpt,(*it)->dipol()->lastPt());
assert((*it)->dipol()->clustersafe());
sum+=Di;
first_subpro.insert(1., (*it));
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
pair<CrossSection,CrossSection> Node::calcDipandPS(Energy scale){
return dipol()->dipandPs(sqr(scale),deepHead()->MH()->largeNBasis());
}
CrossSection Node::calcPs(Energy scale){
return dipol()->ps(sqr(scale),deepHead()->MH()->largeNBasis());
}
/*
bool Node::diffPsDipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
bool isInSomePS=false;
Energy maxx=0.*GeV;
Energy minn=100000*GeV;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
isInSomePS|=(*it)->inShowerPS((*it2)->dipol()->lastPt());
maxx=max(maxx,(*it)->dipol()->lastPt());
minn=min(minn,(*it)->dipol()->lastPt());
}
double Di=0.;
if(isInSomePS||(tmp2.empty())){
Di=-1.* (*it)->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->MH()->largeNBasis())/nanobarn;
}else{
Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
}
if ((Di!=0.)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<MH()->mergePt())) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++) assert(((*it2)->dipol()->lastPt()>deepHead()->MH()->mergePt()));
if (Di!=0.) {
first_subpro.insert(1., (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
vector<Ptr<Node>::ptr> tmp2=selectedNode->children();
bool isInSomePS=false;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
isInSomePS|=selectedNode->inShowerPS((*it2)->dipol()->lastPt());
}
if((isInSomePS||(tmp2.empty()))){//selectedNode->dipol()->lastPt()<MH()->mergePt()&&
sum=-1.* selectedNode->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->MH()->largeNBasis())/nanobarn;
}else{
sum=-1.* selectedNode->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
}
return true;
}
return false;
}
*/
/*
bool Node::psBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
double Di=-1.* (*it)->dipol()->ps(sqr(10.*GeV),deepHead()->MH()->largeNBasis())/nanobarn;
if ((Di!=0)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<MH()->mergePt())) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
bool isInSomePS=false;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
assert(((*it2)->dipol()->lastPt()>deepHead()->MH()->mergePt())||((*it)->dipol()->lastPt()>deepHead()->MH()->mergePt()));
isInSomePS|=(*it)->inShowerPS((*it2)->dipol()->lastPt());
}
if(!isInSomePS&&!(tmp2.empty()))continue;
sum+=Di;
if (Di!=0) {
first_subpro.insert(1., (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
*/
/*
bool Node::dipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
if (true
//(*it)->dipol()->clustersafe() &&
// (*it)->xcomb()->willPassCuts()
) {
bool calcdip=true;
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
assert(((*it2)->dipol()->lastPt()>deepHead()->MH()->mergePt())||((*it)->dipol()->lastPt()>deepHead()->MH()->mergePt()));
if(((*it2)->dipol()->lastPt()<deepHead()->MH()->mergePt())){
if((*it)->dipol()->lastPt()<deepHead()->MH()->mergePt()){
sum=0;
return false;
}
calcdip=false;
}
}
if(calcdip){
double Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
sum+=Di;
minpt=min(minpt,(*it)->dipol()->lastPt());
if (Di!=0) {
first_subpro.insert(1., (*it));
}
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
*/
IBPtr Node::clone() const {
return new_ptr(*this);
}
IBPtr Node::fullclone() const {
return new_ptr(*this);
}
void Node::persistentOutput(PersistentOStream & os) const {
os <<theheadxcomb<<
thexcomb<<
thenodeMEPtr<<
thedipol<<
thechildren<<
theparent<<
theProjectors<<
theDeepHead<<
theCutStage<<
isthissafe<<
theDeepProStage<<
clustersafer<<
isOrdered<<
theOrderdSteps<<
theNeedFullOrderedHistory<<
theSudakovSteps<<
ounit(theVetoPt,GeV)<<
ounit(theRunningPt,GeV)<<
theSubtractedReal<<
theVirtualContribution<<
theMergingHelper;
}
void Node::persistentInput(PersistentIStream & is, int) {
is >>theheadxcomb>>
thexcomb>>
thenodeMEPtr>>
thedipol>>
thechildren>>
theparent>>
theProjectors>>
theDeepHead>>
theCutStage>>
isthissafe>>
theDeepProStage>>
clustersafer>>
isOrdered>>
theOrderdSteps>>
theNeedFullOrderedHistory>>
theSudakovSteps>>
iunit(theVetoPt,GeV)>>
iunit(theRunningPt,GeV)>>
theSubtractedReal>>
theVirtualContribution>>
theMergingHelper;
}
// *** 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<Node,Interfaced> describeHerwigNode("Herwig::Node", "HwDipoleShower.so");
void Node::Init() {
static ClassDocumentation<Node> documentation("There is no documentation for the Node class");
}
diff --git a/src/Merging/2LO-1NLO-LHC-Z.in b/src/Merging/2LO-1NLO-LHC-Z.in
--- a/src/Merging/2LO-1NLO-LHC-Z.in
+++ b/src/Merging/2LO-1NLO-LHC-Z.in
@@ -1,28 +1,32 @@
read Merging/MergingLHC.in
cd /Herwig/Merging
set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
set PPMFactory:RenormalizationScaleFactor 1.
set PPMFactory:FactorizationScaleFactor 1.
set PPMFactory:onlyabove 0
set Merger:chooseHistory 8
set Merger:addNLOLegs 0
set Merger:addLOLegs 1
set Merger:MergingScale 10.*GeV
set Merger:MergingScaleSmearing 0.1
set Merger:gamma 1.
cd /Herwig/Merging
set PPMFactory:OrderInAlphaS 0
set PPMFactory:OrderInAlphaEW 2
do PPMFactory:Process p p e+ e-
+cd /Herwig/MatrixElements/Matchbox/Utility
+insert DiagramGenerator:ExcludeInternal 0 /Herwig/Particles/gamma
+
+
cd /Herwig/Generators
saverun 2LO-1NLO-LHC-Z LHCGenerator
diff --git a/src/Merging/MergingLEP.in b/src/Merging/MergingLEP.in
--- a/src/Merging/MergingLEP.in
+++ b/src/Merging/MergingLEP.in
@@ -1,263 +1,315 @@
set /Herwig/Particles/d:NominalMass 0*GeV
set /Herwig/Particles/dbar:NominalMass 0*GeV
set /Herwig/Particles/u:NominalMass 0*GeV
set /Herwig/Particles/ubar:NominalMass 0*GeV
set /Herwig/Particles/s:NominalMass 0*GeV
set /Herwig/Particles/sbar:NominalMass 0*GeV
set /Herwig/Particles/c:NominalMass 0*GeV
set /Herwig/Particles/cbar:NominalMass 0*GeV
-set /Herwig/Particles/b:NominalMass 0*GeV
-set /Herwig/Particles/bbar:NominalMass 0*GeV
+#set /Herwig/Particles/b:NominalMass 0*GeV
+#set /Herwig/Particles/bbar:NominalMass 0*GeV
#set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
#set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
-
-set /Herwig/DipoleShower/Kinematics/FFLightKinematics:IRCutoff 0.78
set /Herwig/Couplings/NLOAlphaS:min_active_flavours 4
set /Herwig/Couplings/NLOAlphaS:input_alpha_s 0.118
set /Herwig/Model:QCD/RunningAlphaS /Herwig/Couplings/NLOAlphaS
set /Herwig/DipoleShower/DipoleShowerHandler:GlobalAlphaS /Herwig/Couplings/NLOAlphaS
cd /Herwig/Merging
create Herwig::Merger Merger
set Merger:DipoleShowerHandler /Herwig/DipoleShower/DipoleShowerHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper Merger
cd /Herwig/DipoleShower/Kernels
-clear /Herwig/DipoleShower/DipoleShowerHandler:Kernels
-insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ggxDipoleKernel
-insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFqx2qgxDipoleKernel
-insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ddxDipoleKernel
-insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2uuxDipoleKernel
-insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ccxDipoleKernel
-insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ssxDipoleKernel
-insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2bbxDipoleKernel
+#clear /Herwig/DipoleShower/DipoleShowerHandler:Kernels
+#insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ggxDipoleKernel
+#insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFqx2qgxDipoleKernel
+#insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ddxDipoleKernel
+#insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2uuxDipoleKernel
+#insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ccxDipoleKernel
+#insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ssxDipoleKernel
+#insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2bbxDipoleKernel
cd /Herwig/EventHandlers
set /Herwig/Generators/LEPGenerator:EventHandler:LuminosityFunction:Energy 91.2*GeV
cd /Herwig/Merging
insert /Herwig/Generators/LEPGenerator:EventHandler:SubProcessHandlers[0] EEMFactory
set EEMFactory:OrderInAlphaS 0
set EEMFactory:OrderInAlphaEW 2
do EEMFactory:Process e- e+ j j
create Herwig::SimpleColourBasis LargeNColourBasis
set LargeNColourBasis:LargeN On
set Merger:LargeNBasis LargeNColourBasis
cd /Herwig/Generators
set LEPGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
set LEPGenerator:EventHandler:CascadeHandler:MPIHandler NULL
set LEPGenerator:EventHandler:MultipleInteractionHandler NULL
set /Herwig/Samplers/Sampler:Verbose On
set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
set /Herwig/Generators/LEPGenerator:EventHandler:CollisionCuts Off
cd /Herwig/MatrixElements/Matchbox/Phasespace
set TreePhasespace:M0 0.001*GeV
set TreePhasespace:MC 0.00001*GeV
set /Herwig/Merging/EEMFactory:Phasespace TreePhasespace
cd /Herwig/EventHandlers
set LEPHandler:Sampler /Herwig/Samplers/Sampler
cd /Herwig/MatrixElements/Matchbox
cd Utility
insert DiagramGenerator:ExcludeInternal 0 /Herwig/Particles/gamma
cd /Herwig/Samplers
create Herwig::MonacoSampler Monaco
set /Herwig/Samplers/Sampler:BinSampler Monaco
set Monaco:EnhancementFactor 1.2
set Monaco:InitialPoints 1000
set Monaco:LuminosityMapperBins 0
set Monaco:NIterations 4
set Monaco:RemapChannelDimension Yes
set Monaco:RemapperMinSelection 0.0001
set Monaco:RemapperPoints 1000
set Monaco:UseAllIterations No
set Sampler:UpdateAfter 3000
set Sampler:AddUpSamplers Off
set Sampler:GlobalMaximumWeight Off
set Sampler:FlatSubprocesses Off
set Sampler:MinSelection 0.00001
set Sampler:AlmostUnweighted On
set Sampler:BinSampler:Kappa 1.
set Sampler:RunCombinationData Off
set Sampler:Verbose On
set /Herwig/EventHandlers/LEPHandler:Sampler /Herwig/Samplers/Sampler
cd /Herwig/MatrixElements/Matchbox/Scales/
set /Herwig/Merging/MScale:ScaleChoice FixedScale
set FixedScale:FixedScale 91.2*GeV
cd /Herwig/Generators
#set LEPGenerator:EventHandler:DecayHandler NULL
#set LEPGenerator:EventHandler:HadronizationHandler NULL
cd /Herwig/Analysis
clear /Herwig/Generators/LEPGenerator:AnalysisHandlers
create ThePEG::RivetAnalysis RivetAnalysis RivetAnalysis.so
insert /Herwig/Generators/LEPGenerator:AnalysisHandlers 0 RivetAnalysis
cd /Herwig/Generators
set /Herwig/Generators/LEPGenerator:IntermediateOutput Yes
set /Herwig/Generators/LEPGenerator:EventHandler /Herwig/EventHandlers/LEPHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MaxPtIsMuF Yes
cd /Herwig/Cuts
set EECuts:Fuzzy FuzzyTheta
set MatchboxJetMatcher:Factory /Herwig/Merging/EEMFactory
cd /Herwig/Merging
set EEMFactory:MergingHelper Merger
set Merger:MFactory EEMFactory
set MScale:MergingHelper Merger
set Merger:MergingJetFinder /Herwig/Cuts/JetFinder
cd /Herwig/Generators
set LEPGenerator:DebugLevel 1
set LEPGenerator:PrintEvent 10
set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper:minusL No
### From Tests/Rivet/LEP/LEP-91.in
##################################################
# LEP physics parameters (override defaults)
##################################################
##################################################
# select the analyses
##################################################
# Validated
##################################################
# ALEPH charged particle multiplicity
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1991_S2435284
# ALEPH main LEP I QCD summary paper
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1996_S3486095
# ALEPH D*
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1999_S4193598
# OPAL charged hadron analysis
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1994_S2927284
# OPAL Delta++ analysis
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1995_S3198391
# OPAL J/Psi analysis analysis
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1996_S3257789
# ALEPH eta/omega analysis
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2002_S4823664
# OPAL K*0 analysis
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_S3608263
# OPAL flavour specific charged multiplicities etc
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3780481
# OPAL f_0,f_2 and phi production
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3702294
# OPAL gamma,pi0,eta,eta',rho+/-,a0+/-
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3749908
# OPAL K0
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2000_S4418603
# SLD flavour specific charged multiplicities etc
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_1996_S3398250
# SLD flavour specific charged multiplicities etc
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_1999_S3743934
# SLD flavour specific charged multiplicities etc
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_2004_S5693039
# OPAL event shapes and multiplicities at different energies
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2004_S6132243
# ALEPH jet and event shapes at many energies
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2004_S5765862
# OPAL/JADE jet rates at many energies
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 JADE_OPAL_2000_S4300807
# DELPHI strange baryon production
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1995_S3137023
# DELPHI f_0, rho_0 and f_2 production
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1999_S3960137
# OPAL strange baryon production
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_S3396100
# DELPHI tuning paper
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1996_S3430090
# DELPHI b quark
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2002_069_CONF_603
-insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2011_I890503
+#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2011_I890503
# ALEPH b quark
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2001_S4656318
# SLD b quark
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_2002_S4869273
# PDG hadron multiplicities and ratios
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 PDG_HADRON_MULTIPLICITIES
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 PDG_HADRON_MULTIPLICITIES_RATIOS
##################################################
# unvalidated
##################################################
# OPAL 4 jet angles
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2001_S4553896
cd /Herwig/DipoleShower/Kernels
set FFgx2ggxDipoleKernel:CMWScheme On
set FFqx2qgxDipoleKernel:CMWScheme On
set FFgx2ddxDipoleKernel:CMWScheme On
set FFgx2uuxDipoleKernel:CMWScheme On
set FFgx2ccxDipoleKernel:CMWScheme On
set FFgx2ssxDipoleKernel:CMWScheme On
set FFgx2bbxDipoleKernel:CMWScheme On
set FFMgx2ggxDipoleKernel:CMWScheme On
set FFMdx2dgxDipoleKernel:CMWScheme On
set FFMux2ugxDipoleKernel:CMWScheme On
set FFMcx2cgxDipoleKernel:CMWScheme On
set FFMsx2sgxDipoleKernel:CMWScheme On
set FFMbx2bgxDipoleKernel:CMWScheme On
set FFMtx2tgxDipoleKernel:CMWScheme On
set FFMgx2ddxDipoleKernel:CMWScheme On
set FFMgx2uuxDipoleKernel:CMWScheme On
set FFMgx2ccxDipoleKernel:CMWScheme On
set FFMgx2ssxDipoleKernel:CMWScheme On
-set FFMgx2bbxDipoleKernel:CMWScheme On
+set FFMgx2bbxDipoleKernel:CMWScheme Off
+#set FFMbx2bgxDipoleKernel:RenormalizationScaleFreeze 8.4*GeV
+set FFMgx2bbxDipoleKernel:RenormalizationScaleFreeze 8.4*GeV
+set FFMbx2bgxDipoleKernel:VirtualitySplittingScale No
+set FFMgx2bbxDipoleKernel:VirtualitySplittingScale Yes
+
+cd /Herwig/Hadronization
+
+#set /Herwig/DipoleShower/Kinematics/FFLightKinematics:IRCutoff 0.78
+#set /Herwig/DipoleShower/Kinematics/FFMassiveKinematics:IRCutoff 0.1018
+
+#set ClusterFissioner:ClMaxLight 3.30
+#set ClusterFissioner:ClPowLight 2.50
+#set ClusterFissioner:PSplitLight 1.29
+#set ClusterDecayer:ClDirLight 1
+#set ClusterDecayer:ClSmrLight 3.118342
+
+#set ClusterFissioner:ClMaxCharm 3.11*GeV
+#set ClusterFissioner:ClPowCharm 1.62
+#set ClusterFissioner:PSplitCharm 2.54
+#set ClusterDecayer:ClDirCharm 1
+#set ClusterDecayer:ClSmrCharm 0.
+#set HadronSelector:SingleHadronLimitCharm 0.0
+
+#set ClusterFissioner:ClMaxBottom 4.30*GeV
+get ClusterFissioner:ClMaxBottom
+#set ClusterFissioner:ClPowBottom 3.518
+get ClusterFissioner:ClPowBottom
+#set ClusterFissioner:PSplitBottom 1.5
+#set ClusterDecayer:ClDirBottom 1
+#set ClusterDecayer:ClSmrBottom .1
+#set HadronSelector:SingleHadronLimitBottom 0.1
+#set /Herwig/Particles/b:NominalMass 4.0
+#set /Herwig/Particles/b:ConstituentMass 4.3
+#set /Herwig/Particles/bbar:ConstituentMass 4.3
+#set /Herwig/Particles/g:ConstituentMass .7
+get /Herwig/Particles/u:ConstituentMass
+get /Herwig/Particles/d:ConstituentMass
+get /Herwig/Particles/s:ConstituentMass
+
+
+#set HadronSelector:PwtUquark 1.0
+#set HadronSelector:PwtDquark 1.0
+#set HadronSelector:PwtSquark 1.09
+#set HadronSelector:PwtCquark 1.0
+#set HadronSelector:PwtBquark .1
+#set HadronSelector:PwtDIquark 0.66
+#set HadronSelector:SngWt 1.0
+#set HadronSelector:DecWt 1.0
+
+#set /Herwig/Couplings/NLOAlphaS:input_alpha_s 0.118
+#set /Herwig/Couplings/NLOAlphaS:min_active_flavours 3
+#fulldescribe /Herwig/Couplings/NLOAlphaS:exact_evaluation exact
+
+#set /Herwig/Couplings/NLOAlphaS:max_active_flavours 6
diff --git a/src/Merging/MergingLHC.in b/src/Merging/MergingLHC.in
--- a/src/Merging/MergingLHC.in
+++ b/src/Merging/MergingLHC.in
@@ -1,376 +1,376 @@
cd /Herwig/Partons
create ThePEG::LHAPDF PDFSet ThePEGLHAPDF.so
set PDFSet:RemnantHandler HadronRemnants
set /Herwig/Partons/PDFSet:PDFName MMHT2014nlo68cl
set /Herwig/Particles/p+:PDF PDFSet
set /Herwig/Particles/pbar-:PDF PDFSet
set HardLOPDF:PDFName MMHT2014nlo68cl
set HardNLOPDF:PDFName MMHT2014nlo68cl
set ShowerLOPDF:PDFName MMHT2014nlo68cl
set ShowerNLOPDF:PDFName MMHT2014nlo68cl
#set MPIPDF:PDFName MMHT2014nlo68cl
#set RemnantPDF:PDFName MMHT2014nlo68cl
-set /Herwig/DipoleShower/IntrinsicPtGenerator:ValenceIntrinsicPtScale 2.0*GeV
-set /Herwig/DipoleShower/IntrinsicPtGenerator:SeaIntrinsicPtScale 2.0*GeV
+set /Herwig/DipoleShower/IntrinsicPtGenerator:ValenceIntrinsicPtScale 2.4*GeV
+set /Herwig/DipoleShower/IntrinsicPtGenerator:SeaIntrinsicPtScale 2.4*GeV
set /Herwig/DipoleShower/Kinematics/FFLightKinematics:IRCutoff 0.78
set /Herwig/DipoleShower/Kinematics/FILightKinematics:IRCutoff 0.78
set /Herwig/DipoleShower/Kinematics/IFLightKinematics:IRCutoff 0.78
set /Herwig/DipoleShower/Kinematics/IILightKinematics:IRCutoff 0.78
set /Herwig/Model:QCD/RunningAlphaS /Herwig/Couplings/NLOAlphaS
set /Herwig/DipoleShower/DipoleShowerHandler:GlobalAlphaS /Herwig/Couplings/NLOAlphaS
set /Herwig/DipoleShower/DipoleShowerHandler:GlobalAlphaS:input_alpha_s 0.118
set /Herwig/Couplings/NLOAlphaS:input_alpha_s 0.118
read Matchbox/FiveFlavourNoBMassScheme.in
cd /Herwig/Particles
set d:NominalMass 0*GeV
set dbar:NominalMass 0*GeV
set u:NominalMass 0*GeV
set ubar:NominalMass 0*GeV
set s:NominalMass 0*GeV
set sbar:NominalMass 0*GeV
set c:NominalMass 0*GeV
set cbar:NominalMass 0*GeV
set c:HardProcessMass 0*GeV
set cbar:HardProcessMass 0*GeV
set b:HardProcessMass 0*GeV
set bbar:HardProcessMass 0*GeV
set e+:HardProcessMass 0*GeV
set e-:HardProcessMass 0*GeV
set mu+:HardProcessMass 0*GeV
set mu-:HardProcessMass 0*GeV
set nu_e:HardProcessMass 0*GeV
set nu_ebar:HardProcessMass 0*GeV
set nu_mu:HardProcessMass 0*GeV
set nu_mubar:HardProcessMass 0*GeV
set nu_tau:HardProcessMass 0*GeV
set nu_taubar:HardProcessMass 0*GeV
cd /Herwig/Cuts
set MassCut:MinM 66*GeV
set MassCut:MaxM 116*GeV
set /Herwig/Generators/LHCGenerator:EventHandler:CollisionCuts Off
cd /Herwig/MatrixElements/Matchbox/Phasespace
set TreePhasespace:M0 0.01*GeV
set TreePhasespace:MC 0.0001*GeV
set /Herwig/Generators/LHCGenerator:MaxErrors 10000
cd /Herwig/EventHandlers
set /Herwig/EventHandlers/LHCHandler:LuminosityFunction:Energy 7000.0*GeV
set /Herwig/Generators/LHCGenerator:NumberOfEvents 100000000
cd /Herwig/Merging
insert /Herwig/Generators/LHCGenerator:EventHandler:SubProcessHandlers[0] PPMFactory
set PPMFactory:QuarkFlavourDiagonal Yes
set /Herwig/Vertices/FFWMatchboxVertex:Diagonal Yes
clear PPMFactory:ParticleGroup p
do PPMFactory:StartParticleGroup p
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/g
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/u
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/ubar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/d
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/dbar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/s
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/sbar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/c
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/cbar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/b
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/bbar
do PPMFactory:EndParticleGroup
clear PPMFactory:ParticleGroup j
do PPMFactory:StartParticleGroup j
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/g
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/u
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/ubar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/d
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/dbar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/s
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/sbar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/c
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/cbar
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/b
insert PPMFactory:ParticleGroup 0 /Herwig/Particles/bbar
do PPMFactory:EndParticleGroup
cd /Herwig/Merging
create Herwig::Merger Merger
create Herwig::MergingReweight MPreWeight HwDipoleShower.so
#insert PPMFactory:Preweighters 0 MPreWeight
set MPreWeight:HTPower 0
set MPreWeight:MaxPTPower 0
set MPreWeight:OnlyColoured False
create Herwig::SimpleColourBasis LargeNColourBasis
set LargeNColourBasis:LargeN On
set Merger:LargeNBasis LargeNColourBasis
set Merger:DipoleShowerHandler /Herwig/DipoleShower/DipoleShowerHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper Merger
cd /Herwig/Generators
set LHCGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
-#set LHCGenerator:EventHandler:CascadeHandler:MPIHandler NULL
-#set LHCGenerator:EventHandler:MultipleInteractionHandler NULL
-#set LHCGenerator:EventHandler:DecayHandler NULL
-#set LHCGenerator:EventHandler:HadronizationHandler NULL
+set LHCGenerator:EventHandler:CascadeHandler:MPIHandler NULL
+set LHCGenerator:EventHandler:MultipleInteractionHandler NULL
+set LHCGenerator:EventHandler:DecayHandler NULL
+set LHCGenerator:EventHandler:HadronizationHandler NULL
#erase /Herwig/Generators/LHCGenerator:EventHandler:PostSubProcessHandlers[0]
create Herwig::MonacoSampler /Herwig/Samplers/Monaco
set /Herwig/Samplers/Sampler:BinSampler /Herwig/Samplers/Monaco
cd /Herwig/MatrixElements/Matchbox/Phasespace
set TreePhasespace:M0 0.01*GeV
set TreePhasespace:MC 0.0001*GeV
set /Herwig/Merging/PPMFactory:Phasespace TreePhasespace
cd /Herwig/EventHandlers
set LHCHandler:Sampler /Herwig/Samplers/Sampler
cd /Herwig/Samplers
set Monaco:EnhancementFactor 1.2
set Monaco:InitialPoints 1000
set Monaco:LuminosityMapperBins 8
set Monaco:NIterations 4
set Monaco:RemapChannelDimension Yes
set Monaco:RemapperMinSelection 0.0001
set Monaco:RemapperPoints 1000
set Monaco:UseAllIterations No
set Sampler:UpdateAfter 1000
set Sampler:AddUpSamplers Off
set Sampler:GlobalMaximumWeight Off
set Sampler:FlatSubprocesses Off
set Sampler:MinSelection 0.000001
set Sampler:AlmostUnweighted On
set Sampler:BinSampler:Kappa 1.
set Sampler:RunCombinationData Off
set Sampler:Verbose On
cd /Herwig/MatrixElements/Matchbox/Amplitudes
#clear /Herwig/Merging/PPMFactory:Amplitudes
#set GenericProcesses:TreeLevelAmplitude MadGraph
#set GenericProcesses:OneLoopAmplitude OpenLoops
cd /Herwig/Cuts
set MassCut:MinM 66*GeV
set MassCut:MaxM 116*GeV
set /Herwig/Generators/LHCGenerator:MaxErrors 10000
cd /Herwig/EventHandlers
set /Herwig/EventHandlers/LHCHandler:LuminosityFunction:Energy 7000.0*GeV
cd /Herwig/MatrixElements/Matchbox/Scales/
set /Herwig/Merging/MScale:ScaleChoice LeptonPairMassScale
cd /Herwig/Analysis
clear /Herwig/Generators/LHCGenerator:AnalysisHandlers
create ThePEG::RivetAnalysis RivetAnalysis RivetAnalysis.so
insert /Herwig/Generators/LHCGenerator:AnalysisHandlers 0 RivetAnalysis
cd /Herwig/Generators
set /Herwig/Generators/LHCGenerator:IntermediateOutput Yes
set /Herwig/Generators/LHCGenerator:EventHandler /Herwig/EventHandlers/LHCHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MaxPtIsMuF Yes
cd /Herwig/Cuts
set QCDCuts:Fuzzy FuzzyTheta
set MatchboxJetMatcher:Factory /Herwig/Merging/PPMFactory
cd /Herwig/Merging
set PPMFactory:MergingHelper Merger
set Merger:MFactory PPMFactory
set MScale:MergingHelper Merger
set Merger:MergingJetFinder /Herwig/Cuts/JetFinder
cd /Herwig/Generators
set LHCGenerator:DebugLevel 1
set LHCGenerator:PrintEvent 10
set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper:minusL No
##################################################
# select the analyses
##################################################
# General analysis
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 MC_ZINC
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 MC_ZJETS
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 MC_ZKTSPLITTINGS
# ATLAS pT
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2011_S9131140
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2014_I1300647
# ATLAS Z+jets
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2011_I945498
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2013_I1230812
# ATLAS phi*
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1204784
# CMS Z + b-hadron
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2013_I1256943
# CMS Z pt and y
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2012_I941555
# ATLAS Z + bjets
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2014_I1306294
# ATLAS Z
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2011_I928289_Z
# CMS Z AFB
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2013_I1122847
# CMS Z+jets
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2015_I1310737
# ATLAS event shapes in Z events
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2016_I1424838
# ATLAS forwrd backward
insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2015_I1351916_EL
cd /Herwig/DipoleShower/Kernels
set FFgx2ggxDipoleKernel:CMWScheme On
set FFqx2qgxDipoleKernel:CMWScheme On
set FFgx2ddxDipoleKernel:CMWScheme On
set FFgx2uuxDipoleKernel:CMWScheme On
set FFgx2ccxDipoleKernel:CMWScheme On
set FFgx2ssxDipoleKernel:CMWScheme On
set FFgx2bbxDipoleKernel:CMWScheme On
set FFMgx2ggxDipoleKernel:CMWScheme On
set FFMdx2dgxDipoleKernel:CMWScheme On
set FFMux2ugxDipoleKernel:CMWScheme On
set FFMcx2cgxDipoleKernel:CMWScheme On
set FFMsx2sgxDipoleKernel:CMWScheme On
set FFMbx2bgxDipoleKernel:CMWScheme On
set FFMtx2tgxDipoleKernel:CMWScheme On
set FFMgx2ddxDipoleKernel:CMWScheme On
set FFMgx2uuxDipoleKernel:CMWScheme On
set FFMgx2ccxDipoleKernel:CMWScheme On
set FFMgx2ssxDipoleKernel:CMWScheme On
set FFMgx2bbxDipoleKernel:CMWScheme On
set FIgx2ggxDipoleKernel:CMWScheme On
set FIqx2qgxDipoleKernel:CMWScheme On
set FIgx2ddxDipoleKernel:CMWScheme On
set FIgx2uuxDipoleKernel:CMWScheme On
set FIgx2ccxDipoleKernel:CMWScheme On
set FIgx2ssxDipoleKernel:CMWScheme On
set FIgx2bbxDipoleKernel:CMWScheme On
set FIMdx2dgxDipoleKernel:CMWScheme On
set FIMux2ugxDipoleKernel:CMWScheme On
set FIMcx2cgxDipoleKernel:CMWScheme On
set FIMsx2sgxDipoleKernel:CMWScheme On
set FIMbx2bgxDipoleKernel:CMWScheme On
set FIMtx2tgxDipoleKernel:CMWScheme On
set FIMgx2ddxDipoleKernel:CMWScheme On
set FIMgx2uuxDipoleKernel:CMWScheme On
set FIMgx2ccxDipoleKernel:CMWScheme On
set FIMgx2ssxDipoleKernel:CMWScheme On
set FIMgx2bbxDipoleKernel:CMWScheme On
#set FIMgx2ttxDipoleKernel:CMWScheme On
set IFgx2ggxDipoleKernel:CMWScheme On
set IFqx2qgxDipoleKernel:CMWScheme On
set IFqx2gqxDipoleKernel:CMWScheme On
set IFgx2ddbarxDipoleKernel:CMWScheme On
set IFgx2dbardxDipoleKernel:CMWScheme On
set IFgx2uubarxDipoleKernel:CMWScheme On
set IFgx2ubaruxDipoleKernel:CMWScheme On
set IFgx2ccbarxDipoleKernel:CMWScheme On
set IFgx2cbarcxDipoleKernel:CMWScheme On
set IFgx2ssbarxDipoleKernel:CMWScheme On
set IFgx2sbarsxDipoleKernel:CMWScheme On
set IFMgx2ggxDipoleKernel:CMWScheme On
set IFMqx2qgxDipoleKernel:CMWScheme On
set IFMqx2gqxDipoleKernel:CMWScheme On
set IFMgx2ddbarxDipoleKernel:CMWScheme On
set IFMgx2dbardxDipoleKernel:CMWScheme On
set IFMgx2uubarxDipoleKernel:CMWScheme On
set IFMgx2ubaruxDipoleKernel:CMWScheme On
set IFMgx2ccbarxDipoleKernel:CMWScheme On
set IFMgx2cbarcxDipoleKernel:CMWScheme On
set IFMgx2ssbarxDipoleKernel:CMWScheme On
set IFMgx2sbarsxDipoleKernel:CMWScheme On
set IIgx2ggxDipoleKernel:CMWScheme On
set IIqx2qgxDipoleKernel:CMWScheme On
set IIqx2gqxDipoleKernel:CMWScheme On
set IIgx2ddbarxDipoleKernel:CMWScheme On
set IIgx2dbardxDipoleKernel:CMWScheme On
set IIgx2uubarxDipoleKernel:CMWScheme On
set IIgx2ubaruxDipoleKernel:CMWScheme On
set IIgx2ccbarxDipoleKernel:CMWScheme On
set IIgx2cbarcxDipoleKernel:CMWScheme On
set IIgx2ssbarxDipoleKernel:CMWScheme On
set IIgx2sbarsxDipoleKernel:CMWScheme On

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:55 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3800223
Default Alt Text
(235 KB)

Event Timeline