Page MenuHomeHEPForge
Files F7877743

No OneTemporary
Actions

View Options

diff --git a/DipoleShower/Kinematics/FFMassiveKinematics.cc b/DipoleShower/Kinematics/FFMassiveKinematics.cc
--- a/DipoleShower/Kinematics/FFMassiveKinematics.cc
+++ b/DipoleShower/Kinematics/FFMassiveKinematics.cc
@@ -1,350 +1,343 @@
// -*- C++ -*-
//
// FFMassiveKinematics.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 FFMassiveKinematics class.
//
#include "FFMassiveKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"
// TODO: remove after verification
// only for checking for NaN or inf
#include <gsl/gsl_math.h>
using namespace Herwig;
FFMassiveKinematics::FFMassiveKinematics()
: DipoleSplittingKinematics() {}
FFMassiveKinematics::~FFMassiveKinematics() {}
IBPtr FFMassiveKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FFMassiveKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> FFMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return make_pair(0.0,1.0);
}
pair<double,double> FFMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emissionData()->id() != ParticleID::g )
return make_pair(0.5*(1.-c),0.5*(1.+c));
double b = log((1.+c)/(1.-c));
return make_pair(-b,b);
}
return make_pair(-log(0.5*(1.+c)),-log(0.5*(1.-c)));
}
Energy FFMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
return (pEmitter+pSpectator).m();
}
Energy FFMassiveKinematics::ptMax(Energy dScale,
double, double,
const DipoleIndex& ind,
const DipoleSplittingKernel& split) const {
double mui = split.emitter(ind)->mass() / dScale;
double mu = split.emission(ind)->mass() / dScale;
double muj = split.spectator(ind)->mass() / dScale;
double mui2 = sqr( mui ), mu2 = sqr( mu ), muj2 = sqr( muj );
-
- if(mui>(1.+mu-muj)&&
- ( mui2*mui2 + mu2*mu2 + sqr(1.-muj)*sqr(1.-muj) - 2.*( mui2*mu2+mui2*sqr(1.-muj)+mu2*sqr(1.-muj) ))>0
- ) {
- return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.*sqrt(muj2) ) * dScale;
- }
- else
- return ZERO;
+ return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.*sqrt(muj2) ) * dScale;
}
Energy FFMassiveKinematics::QMax(Energy dScale,
double, double,
const DipoleIndex& ind,
const DipoleSplittingKernel&) const {
assert(false && "implementation missing");
double Muj = ind.spectatorData()->mass() / dScale;
return dScale * ( 1.-2.*Muj+sqr(Muj) );
}
Energy FFMassiveKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 pt2 = z*(1.-z)*sqr(scale) - (1-z)*sqr(mi) - z*sqr(m);
assert(pt2 >= ZERO);
return sqrt(pt2);
}
Energy FFMassiveKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 Q2 = (sqr(scale) + (1-z)*sqr(mi) + z*sqr(m))/(z*(1.-z));
return sqrt(Q2);
}
double FFMassiveKinematics::ptToRandom(Energy pt, Energy,
double,double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
// own kinematics close to Dinsdale,Ternick,Weinzierl
bool FFMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info,
const DipoleSplittingKernel&) {
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
if ( pt > info.hardPt() || pt < IRCutoff() ) {
jacobian(0.0);
return false;
}
double z;
double mapZJacobian;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emissionData()->id() != ParticleID::g ) {
z = xi;
mapZJacobian = 1.;
} else {
z = exp(xi)/(1.+exp(xi));
mapZJacobian = z*(1.-z);
}
} else {
z = 1.-exp(-xi);
mapZJacobian = 1.-z;
}
// masses
double mui2 = sqr( info.emitterData()->mass() / info.scale() );
double mu2 = sqr( info.emissionData()->mass() / info.scale() );
double muj2 = sqr( info.spectatorData()->mass() / info.scale() );
double Mui2 = 0.;
if ( info.emitterData()->id() + info.emissionData()->id() == 0 ) Mui2 = 0.; // gluon
else Mui2 = mui2; // (anti)quark
double Muj2 = muj2;
// new: 2011-08-31
// 2011-11-08: this does happen
if( sqrt(mui2)+sqrt(mu2)+sqrt(muj2) > 1. ){
jacobian(0.0);
return false;
}
double bar = 1.-mui2-mu2-muj2;
double y = ( sqr( pt / info.scale() ) + sqr(1.-z)*mui2 + z*z*mu2 ) /
(z*(1.-z)*bar);
// do this here for simplicity
Energy ptmax1 = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
( 2.-2.*sqrt(muj2) ) * info.scale();
Energy auxHardPt = ptmax1 > info.hardPt() ? info.hardPt() : ptmax1;
// 2011-11-09
assert(ptmax1>=info.hardPt());
// phasespace constraint to incorporate ptMax
double zp1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) +
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
sqrt( 1.-sqr(pt/auxHardPt) ) ) /
( 2.*sqr(1.-sqrt(muj2)) );
double zm1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) -
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
sqrt( 1.-sqr(pt/auxHardPt) ) ) /
( 2.*sqr(1.-sqrt(muj2)) );
if ( z > zp1 || z < zm1 ||
pt > auxHardPt) {
jacobian(0.0);
return false;
}
// kinematic phasespace boundaries for (y,z)
// same as in Dittmaier hep-ph/9904440v2 (equivalent to CS)
double ym = 2.*sqrt(mui2)*sqrt(mu2)/bar;
double yp = 1. - 2.*sqrt(muj2)*(1.-sqrt(muj2))/bar;
if ( y < ym || y > yp ) {
jacobian(0.0);
return false;
}
double zm = ( (2.*mui2+bar*y)*(1.-y) - sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
double zp = ( (2.*mui2+bar*y)*(1.-y) + sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
if ( z < zm || z > zp ) {
jacobian(0.0);
return false;
}
double phi = 2.*Constants::pi*rphi;
jacobian( 2. * mapZJacobian * (1.-y) *
log(0.5 * generator()->maximumCMEnergy()/IRCutoff()) *
bar / rootOfKallen(1.,Mui2,Muj2) );
lastPt(pt);
lastZ(z);
lastPhi(phi);
if ( theMCCheck )
theMCCheck->book(1.,1.,info.scale(),info.hardPt(),pt,z,jacobian());
return true;
}
// kinematics close to Dinsdale,Ternick,Weinzierl
// revised 2011-08-22
// revised 2011-11-06
void FFMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
double z = dInfo.lastZ();
Energy pt = dInfo.lastPt();
// masses
double mui2 = sqr( dInfo.emitterData()->mass() / dInfo.scale() );
double mu2 = sqr( dInfo.emissionData()->mass() / dInfo.scale() );
double muj2 = sqr( dInfo.spectatorData()->mass() / dInfo.scale() );
double y = ( sqr( pt / dInfo.scale() ) + sqr(1.-z)*mui2 + z*z*mu2 ) /
(z*(1.-z)*(1.-mui2-mu2-muj2));
Energy2 sbar = sqr(dInfo.scale()) *(1.-mui2-mu2-muj2);
// CMF: particle energies
Energy Ei = ( sbar*(1.-(1.-z)*(1.-y)) + 2.*sqr(dInfo.scale())*mui2 ) / (2.*dInfo.scale());
Energy E = ( sbar*(1.- z *(1.-y)) + 2.*sqr(dInfo.scale())*mu2 ) / (2.*dInfo.scale());
Energy Ej = ( sbar*(1.- y ) + 2.*sqr(dInfo.scale())*muj2 ) / (2.*dInfo.scale());
// CMF: momenta in z-direction (axis of pEmitter & pSpectator)
Energy qi3 = (2.*Ei*Ej-z*(1.-y)*sbar ) / 2./sqrt(Ej*Ej-sqr(dInfo.scale())*muj2);
Energy q3 = (2.*E *Ej-(1.-z)*(1.-y)*sbar) / 2./sqrt(Ej*Ej-sqr(dInfo.scale())*muj2);
Energy qj3 = sqrt( sqr(Ej) - sqr(dInfo.scale())*muj2 );
// get z axis in the dipole's CMF which is parallel to pSpectator
Boost toCMF = (pEmitter+pSpectator).findBoostToCM();
Lorentz5Momentum pjAux = pSpectator; pjAux.boost(toCMF);
ThreeVector<double> pjAxis = pjAux.vect().unit();
// set the momenta in this special reference frame
// note that pt might in some cases differ from the physical pt!
// phi is defined exactly as in getKt
Energy ptResc = sqrt( sqr(Ei)-sqr(dInfo.scale())*mui2-sqr(qi3) );
Lorentz5Momentum em ( ptResc*cos(dInfo.lastPhi()), -ptResc*sin(dInfo.lastPhi()), qi3, Ei );
Lorentz5Momentum emm ( -ptResc*cos(dInfo.lastPhi()), ptResc*sin(dInfo.lastPhi()), q3, E );
Lorentz5Momentum spe ( 0.*GeV, 0.*GeV, qj3, Ej );
// output the mismatch between pt and physical pt
// ofstream output1("ptDiffOnPtAxis-uub-m.dat",ofstream::app);
// ofstream output2("ptDiffOnCosAxis-uub-m.dat",ofstream::app);
// if( abs(dInfo.spectatorData()->id())==5 && dInfo.emitterData()->id()+dInfo.emissionData()->id()==0 &&
// abs(dInfo.emitterData()->id())==1 ) {
// output1 << pt/dInfo.scale() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
// output2 << em.vect().unit()*emm.vect().unit() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
// }
// output1.close(); output2.close();
// rotate back
em.rotateUz (pjAxis);
emm.rotateUz(pjAxis);
spe.rotateUz(pjAxis);
// boost back
em.boost (-toCMF);
emm.boost(-toCMF);
spe.boost(-toCMF);
// mass shells, rescale energy
em.setMass(dInfo.scale()*sqrt(mui2));
em.rescaleEnergy();
emm.setMass(dInfo.scale()*sqrt(mu2));
emm.rescaleEnergy();
spe.setMass(dInfo.scale()*sqrt(muj2));
spe.rescaleEnergy();
// book
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
// TODO: remove
// 2011-11-09: never occurred
if(em.t()/GeV>=0. && emm.t()/GeV>=0. && spe.t()/GeV>=0.);
else cout << "FFMassiveKinematics::generateKinematics momenta corrupt" << endl;
// 2011-11-03 LEP run with full masses:
// x,y,t no problem
// z order > 5.e-7 happend 41 times in 10000 runs
// maximum was 2.5e-6
// double order=2.4e-6;
// Lorentz5Momentum pDif=em+emm+spe-pEmitter-pSpectator, pSum=em+emm+spe+pEmitter+pSpectator;
// if(abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) <order || (pDif-pSum).x()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: x " <<
// abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) << endl <<
// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
// if(abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) <order || (pDif-pSum).y()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: y " <<
// abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) << endl;
// if(abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) <order | (pDif-pSum).z()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: z " <<
// abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) << endl <<
// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
// if(abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) <order || (pDif-pSum).t()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: t " <<
// // abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) << endl;
// cout << endl;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFMassiveKinematics::persistentOutput(PersistentOStream & ) const {
}
void FFMassiveKinematics::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FFMassiveKinematics> FFMassiveKinematics::initFFMassiveKinematics;
// Definition of the static class description member.
void FFMassiveKinematics::Init() {
static ClassDocumentation<FFMassiveKinematics> documentation
("FFMassiveKinematics implements massive splittings "
"off a final-final dipole.");
}
diff --git a/DipoleShower/Kinematics/FFMassiveKinematics.h b/DipoleShower/Kinematics/FFMassiveKinematics.h
--- a/DipoleShower/Kinematics/FFMassiveKinematics.h
+++ b/DipoleShower/Kinematics/FFMassiveKinematics.h
@@ -1,261 +1,262 @@
// -*- C++ -*-
//
// FFMassiveKinematics.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_FFMassiveKinematics_H
#define HERWIG_FFMassiveKinematics_H
//
// This is the declaration of the FFMassiveKinematics class.
//
#include "DipoleSplittingKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief FFMassiveKinematics implements massive splittings
* off a final-final dipole.
*
*/
class FFMassiveKinematics: public DipoleSplittingKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FFMassiveKinematics();
/**
* The destructor.
*/
virtual ~FFMassiveKinematics();
//@}
public:
/**
* Return the boundaries in between the evolution
* variable random number is to be sampled; the lower
* cuoff is assumed to correspond to the infrared cutoff.
*/
virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries in between the momentum
* fraction random number is to be sampled.
*/
virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries on the momentum fraction
*/
virtual pair<double,double> zBoundaries(Energy,
const DipoleSplittingInfo&,
const DipoleSplittingKernel&) const {
return make_pair(0.0,1.0);
}
/**
* Return the dipole scale associated to the
* given pair of emitter and spectator. This
* should be the invariant mass or absolute value
* final/final or initial/initial and the absolute
* value of the momentum transfer for intial/final or
* final/initial dipoles.
*/
virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy dScale,
double emX, double specX,
const DipoleIndex& ind,
const DipoleSplittingKernel& split) const;
/**
* Return the maximum virtuality for the given dipole scale.
*/
virtual Energy QMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel& split) const;
/**
* Return the pt given a virtuality.
*/
virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the virtuality given a pt.
*/
virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the random number associated to
* the given pt.
*/
virtual double ptToRandom(Energy pt, Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel& split) const;
/**
* Generate splitting variables given three random numbers
* and the momentum fractions of the emitter and spectator.
* Return true on success.
*/
virtual bool generateSplitting(double kappa, double xi, double phi,
DipoleSplittingInfo& dIndex,
const DipoleSplittingKernel& split);
/**
* Generate the full kinematics given emitter and
* spectator momentum and a previously completeted
* DipoleSplittingInfo object.
*/
virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo);
public:
/**
* Triangular / Kallen function
*/
template <class T>
inline double rootOfKallen (T a, T b, T c) const {
- return sqrt( a*a + b*b + c*c - 2.*( a*b+a*c+b*c ) ); }
+ double sres=a*a + b*b + c*c - 2.*( a*b+a*c+b*c );
+ return sres>0.?sqrt( sres ):0.; }
/**
* Perform a rotation on both momenta such that the first one will
* point along the (positive) z axis. Rotate back to the original
* reference frame by applying rotateUz(returnedVector) to each momentum.
*/
ThreeVector<double> rotateToZ (Lorentz5Momentum& pTarget, Lorentz5Momentum& p1){
ThreeVector<double> oldAxis = pTarget.vect().unit();
double ct = oldAxis.z(); double st = sqrt( 1.-sqr(ct) ); // cos,sin(theta)
double cp = oldAxis.x()/st; double sp = oldAxis.y()/st; // cos,sin(phi)
pTarget.setZ( pTarget.vect().mag() ); pTarget.setX( 0.*GeV ); pTarget.setY( 0.*GeV );
Lorentz5Momentum p1old = p1;
p1.setX( sp*p1old.x() - cp*p1old.y() );
p1.setY( ct*cp*p1old.x() + ct*sp*p1old.y() - st*p1old.z() );
p1.setZ( st*cp*p1old.x() + st*sp*p1old.y() + ct*p1old.z() );
return oldAxis;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<FFMassiveKinematics> initFFMassiveKinematics;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FFMassiveKinematics & operator=(const FFMassiveKinematics &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of FFMassiveKinematics. */
template <>
struct BaseClassTrait<Herwig::FFMassiveKinematics,1> {
/** Typedef of the first base class of FFMassiveKinematics. */
typedef Herwig::DipoleSplittingKinematics NthBase;
};
/** This template specialization informs ThePEG about the name of
* the FFMassiveKinematics class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::FFMassiveKinematics>
: public ClassTraitsBase<Herwig::FFMassiveKinematics> {
/** Return a platform-independent class name */
static string className() { return "Herwig::FFMassiveKinematics"; }
/**
* The name of a file containing the dynamic library where the class
* FFMassiveKinematics is implemented. It may also include several, space-separated,
* libraries if the class FFMassiveKinematics 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_FFMassiveKinematics_H */
diff --git a/DipoleShower/Kinematics/FIMassiveKinematics.cc b/DipoleShower/Kinematics/FIMassiveKinematics.cc
--- a/DipoleShower/Kinematics/FIMassiveKinematics.cc
+++ b/DipoleShower/Kinematics/FIMassiveKinematics.cc
@@ -1,282 +1,279 @@
// -*- C++ -*-
//
// FIMassiveKinematics.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 FIMassiveKinematics class.
//
#include "FIMassiveKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"
using namespace Herwig;
FIMassiveKinematics::FIMassiveKinematics()
: DipoleSplittingKinematics() {}
FIMassiveKinematics::~FIMassiveKinematics() {}
IBPtr FIMassiveKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FIMassiveKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> FIMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return make_pair(0.0,1.0);
}
pair<double,double> FIMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emissionData()->id() != ParticleID::g )
return make_pair(0.5*(1.-c),0.5*(1.+c));
double b = log((1.+c)/(1.-c));
return make_pair(-b,b);
}
return make_pair(-log(0.5*(1.+c)),-log(0.5*(1.-c)));
}
// sbar
Energy FIMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
return sqrt(2.*(pEmitter*pSpectator));
}
Energy FIMassiveKinematics::ptMax(Energy dScale,
double, double specX,
const DipoleIndex& ind,
const DipoleSplittingKernel& split) const {
Energy mi = split.emitter(ind)->mass(), m = split.emission(ind)->mass();
Energy2 mi2 = sqr(mi), m2 = sqr(m);
// Energy2 Mi2 = split.emitter(int)->id() + split.emission(int)->id() == 0 ?
// 0.*GeV2 : mi2;
Energy2 Mi2 = mi2 == m2 ? 0.*GeV2 : mi2;
// s^star/x
Energy2 s = sqr(dScale) * (1.-specX)/specX + Mi2;
Energy roots = sqrt(s);
- if(roots>mi+m)
- return .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
- else
- return ZERO;
+ return .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
}
// what is this? in FF it is given by y+*dScale = sqrt( 2qi*q / bar )->max
Energy FIMassiveKinematics::QMax(Energy dScale,
double, double specX,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
assert(false && "implementation missing");
cout << "FIMassiveKinematics::QMax called.\n" << flush;
// this is sqrt( 2qi*q ) -> max;
return dScale * sqrt((1.-specX)/specX);
}
Energy FIMassiveKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 pt2 = z*(1.-z)*sqr(scale) - (1-z)*sqr(mi) - z*sqr(m);
assert(pt2 >= ZERO);
return sqrt(pt2);
}
Energy FIMassiveKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 Q2 = (sqr(scale) + (1-z)*sqr(mi) + z*sqr(m))/(z*(1.-z));
return sqrt(Q2);
}
double FIMassiveKinematics::ptToRandom(Energy pt, Energy,
double,double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
bool FIMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info,
const DipoleSplittingKernel&) {
if ( info.spectatorX() < xMin() ) {
jacobian(0.0);
return false;
}
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
if ( pt > info.hardPt() || pt < IRCutoff() ) {
jacobian(0.0);
return false;
}
double z;
double mapZJacobian;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emissionData()->id() != ParticleID::g ) {
z = xi;
mapZJacobian = 1.;
} else {
z = exp(xi)/(1.+exp(xi));
mapZJacobian = z*(1.-z);
}
} else {
z = 1.-exp(-xi);
mapZJacobian = 1.-z;
}
// double s = z*(1.-z);
// double xs = info.spectatorX();
// double x = 1. / ( 1. + sqr(pt/info.scale()) / s );
// double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
// double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));
Energy2 mi2 = sqr(info.emitterData()->mass());
Energy2 m2 = sqr(info.emissionData()->mass());
Energy2 Mi2 = info.emitterData()->id()+info.emissionData()->id() == 0 ?
0.*GeV2 : mi2;
// s^star/x
Energy2 s = sqr(info.scale()) * (1.-info.spectatorX())/info.spectatorX() + Mi2;
double xs = info.spectatorX();
double x = 1. / ( 1. +
( sqr(pt) + (1.-z)*mi2 + z*m2 - z*(1.-z)*Mi2 ) /
( z*(1.-z)*s ) );
double zm1 = .5*( 1.+(mi2-m2)/s - rootOfKallen(s/s,mi2/s,m2/s) *
sqrt( 1.-sqr(pt/info.hardPt()) ) );
double zp1 = .5*( 1.+(mi2-m2)/s + rootOfKallen(s/s,mi2/s,m2/s) *
sqrt( 1.-sqr(pt/info.hardPt()) ) );
if ( // pt < IRCutoff() ||
// pt > info.hardPt() ||
z > zp1 || z < zm1 ||
x < xs ) {
jacobian(0.0);
return false;
}
// additional purely kinematic constraints
double mui2 = x*mi2/sqr(info.scale());
double mu2 = x*m2/sqr(info.scale());
double Mui2 = x*Mi2/sqr(info.scale());
double xp = 1. + Mui2 - sqr(sqrt(mui2)+sqrt(mu2));
double root = sqr(1.-x+Mui2-mui2-mu2)-4.*mui2*mu2;
if( root < 0. && root>-1e-10 )
root = 0.;
else if (root <0. ) {
jacobian(0.0);
return false;
}
root = sqrt(root);
double zm = .5*( 1.-x+Mui2+mui2-mui2 - root ) / (1.-x+Mui2);
double zp = .5*( 1.-x+Mui2+mui2-mui2 + root ) / (1.-x+Mui2);
if (x > xp ||
z > zp || z < zm ) {
jacobian(0.0);
return false;
}
double phi = 2.*Constants::pi*rphi;
jacobian( 2. * mapZJacobian * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
lastPt(pt);
lastZ(z);
lastPhi(phi);
lastSpectatorZ(x);
if ( theMCCheck )
theMCCheck->book(1.,info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());
return true;
}
void FIMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
Energy pt = dInfo.lastPt();
double z = dInfo.lastZ();
Lorentz5Momentum kt =
getKt (pSpectator, pEmitter, pt, dInfo.lastPhi(),true);
Energy2 mi2 = sqr(dInfo.emitterData()->mass());
Energy2 m2 = sqr(dInfo.emissionData()->mass());
Energy2 Mi2 = dInfo.emitterData()->id() + dInfo.emissionData()->id() == 0 ?
0.*GeV2 : mi2;
double xInv = ( 1. +
(pt*pt+(1.-z)*mi2+z*m2-z*(1.-z)*Mi2) /
(z*(1.-z)*sqr(dInfo.scale())) );
Lorentz5Momentum em = z*pEmitter +
(sqr(pt)+mi2-z*z*Mi2)/(z*sqr(dInfo.scale()))*pSpectator + kt;
em.setMass(sqrt(mi2));
em.rescaleEnergy();
Lorentz5Momentum emm = (1.-z)*pEmitter +
(pt*pt+m2-sqr(1.-z)*Mi2)/((1.-z)*sqr(dInfo.scale()))*pSpectator - kt;
emm.setMass(sqrt(m2));
emm.rescaleEnergy();
Lorentz5Momentum spe = xInv*pSpectator;
spe.setMass(ZERO);
spe.rescaleEnergy();
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMassiveKinematics::persistentOutput(PersistentOStream & ) const {
}
void FIMassiveKinematics::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FIMassiveKinematics> FIMassiveKinematics::initFIMassiveKinematics;
// Definition of the static class description member.
void FIMassiveKinematics::Init() {
static ClassDocumentation<FIMassiveKinematics> documentation
("FIMassiveKinematics implements massless splittings "
"off a final-initial dipole.");
}
diff --git a/DipoleShower/Kinematics/FIMassiveKinematics.h b/DipoleShower/Kinematics/FIMassiveKinematics.h
--- a/DipoleShower/Kinematics/FIMassiveKinematics.h
+++ b/DipoleShower/Kinematics/FIMassiveKinematics.h
@@ -1,244 +1,246 @@
// -*- C++ -*-
//
// FILightKinematics.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_FILightKinematics_H
#define HERWIG_FILightKinematics_H
//
// This is the declaration of the FILightKinematics class.
//
#include "DipoleSplittingKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer, Martin Stoll
*
* \brief FIMassiveKinematics implements massless splittings
* off a final-initial dipole.
*
*/
class FIMassiveKinematics: public DipoleSplittingKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FIMassiveKinematics();
/**
* The destructor.
*/
virtual ~FIMassiveKinematics();
//@}
public:
/**
* Return the boundaries in between the evolution
* variable random number is to be sampled; the lower
* cuoff is assumed to correspond to the infrared cutoff.
*/
virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries in between the momentum
* fraction random number is to be sampled.
*/
virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries on the momentum fraction
*/
virtual pair<double,double> zBoundaries(Energy,
const DipoleSplittingInfo&,
const DipoleSplittingKernel&) const {
return make_pair(0.0,1.0);
}
/**
* Return the dipole scale associated to the
* given pair of emitter and spectator. This
* should be the invariant mass or absolute value
* final/final or initial/initial and the absolute
* value of the momentum transfer for intial/final or
* final/initial dipoles.
*/
virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel&) const;
/**
* Return the maximum virtuality for the given dipole scale.
*/
virtual Energy QMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel&) const;
/**
* Return the pt given a virtuality.
*/
virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the virtuality given a pt.
*/
virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the random number associated to
* the given pt.
*/
virtual double ptToRandom(Energy pt, Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel&) const;
/**
* Generate splitting variables given three random numbers
* and the momentum fractions of the emitter and spectator.
* Return true on success.
*/
virtual bool generateSplitting(double kappa, double xi, double phi,
DipoleSplittingInfo& dIndex,
const DipoleSplittingKernel&);
/**
* Generate the full kinematics given emitter and
* spectator momentum and a previously completeted
* DipoleSplittingInfo object.
*/
virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo);
public:
/**
* Triangular / Kallen function
*/
template <class T>
inline double rootOfKallen (T a, T b, T c) const {
- return sqrt( a*a + b*b + c*c - 2.*( a*b+a*c+b*c ) ); }
+ double sres=a*a + b*b + c*c - 2.*( a*b+a*c+b*c );
+ return sres>0.?sqrt( sres ):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 static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<FIMassiveKinematics> initFIMassiveKinematics;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FIMassiveKinematics & operator=(const FIMassiveKinematics &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of FIMassiveKinematics. */
template <>
struct BaseClassTrait<Herwig::FIMassiveKinematics,1> {
/** Typedef of the first base class of FIMassiveKinematics. */
typedef Herwig::DipoleSplittingKinematics NthBase;
};
/** This template specialization informs ThePEG about the name of
* the FIMassiveKinematics class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::FIMassiveKinematics>
: public ClassTraitsBase<Herwig::FIMassiveKinematics> {
/** Return a platform-independent class name */
static string className() { return "Herwig::FIMassiveKinematics"; }
/**
* The name of a file containing the dynamic library where the class
* FIMassiveKinematics is implemented. It may also include several, space-separated,
* libraries if the class FIMassiveKinematics 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_FIMassiveKinematics_H */
diff --git a/MatrixElement/Matchbox/Matching/DipoleMatching.cc b/MatrixElement/Matchbox/Matching/DipoleMatching.cc
--- a/MatrixElement/Matchbox/Matching/DipoleMatching.cc
+++ b/MatrixElement/Matchbox/Matching/DipoleMatching.cc
@@ -1,129 +1,139 @@
// -*- C++ -*-
//
// DipoleMatching.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 DipoleMatching class.
//
#include "DipoleMatching.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
using namespace Herwig;
DipoleMatching::DipoleMatching() {}
DipoleMatching::~DipoleMatching() {}
IBPtr DipoleMatching::clone() const {
return new_ptr(*this);
}
IBPtr DipoleMatching::fullclone() const {
return new_ptr(*this);
}
CrossSection DipoleMatching::dSigHatDR() const {
double xme2 = 0.;
pair<int,int> ij(dipole()->bornEmitter(),
dipole()->bornSpectator());
double ccme2 =
dipole()->underlyingBornME()->largeNColourCorrelatedME2(ij,theLargeNBasis);
+ if(ccme2==0.)return 0.*nanobarn;
+
+ double lnme2=dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
+ if(lnme2==0){
+ generator()->log() <<"\nQTildeMatching: ";
+ generator()->log() <<"\n LargeNME2 is ZERO, while largeNColourCorrelatedME2 is not ZERO." ;
+ generator()->log() <<"\n This is too seriuos.\n" ;
+ generator()->log() << Exception::runerror;
+ }
+
+
ccme2 *=
- dipole()->underlyingBornME()->me2() /
- dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
+ dipole()->underlyingBornME()->me2() /lnme2;
xme2 = dipole()->me2Avg(ccme2);
xme2 /= dipole()->underlyingBornME()->lastXComb().lastAlphaS();
double bornPDF = bornPDFWeight(dipole()->underlyingBornME()->lastScale());
if ( bornPDF == 0.0 )
return ZERO;
xme2 *= bornPDF;
if ( profileScales() )
xme2 *= profileScales()->hardScaleProfile(dipole()->showerHardScale(),dipole()->lastPt());
CrossSection res =
sqr(hbarc) *
realXComb()->jacobian() *
subtractionScaleWeight() *
xme2 /
(2. * realXComb()->lastSHat());
return res;
}
double DipoleMatching::me2() const {
throw Exception() << "DipoleMatching::me2(): Not intented to use. Disable the ShowerApproximationGenerator."
<< Exception::runerror;
return 0.;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleMatching::persistentOutput(PersistentOStream & os) const {
os << theShowerHandler;
}
void DipoleMatching::persistentInput(PersistentIStream & is, int) {
is >> theShowerHandler;
}
void DipoleMatching::doinit() {
ShowerApproximation::doinit();
if ( theShowerHandler ) {
if ( theShowerHandler->scaleFactorOption() < 2 ) {
hardScaleFactor(theShowerHandler->hardScaleFactor());
factorizationScaleFactor(theShowerHandler->factorizationScaleFactor());
renormalizationScaleFactor(theShowerHandler->renormalizationScaleFactor());
}
profileScales(theShowerHandler->profileScales());
restrictPhasespace(theShowerHandler->restrictPhasespace());
hardScaleIsMuF(theShowerHandler->hardScaleIsMuF());
}
}
// *** 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<DipoleMatching,Herwig::ShowerApproximation>
describeHerwigDipoleMatching("Herwig::DipoleMatching", "HwDipoleMatching.so HwShower.so");
void DipoleMatching::Init() {
static ClassDocumentation<DipoleMatching> documentation
("DipoleMatching implements NLO matching with the dipole shower.");
static Reference<DipoleMatching,ShowerHandler> interfaceShowerHandler
("ShowerHandler",
"",
&DipoleMatching::theShowerHandler, false, false, true, true, false);
}
diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.cc b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
--- a/MatrixElement/Matchbox/Matching/QTildeMatching.cc
+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
@@ -1,501 +1,511 @@
// -*- C++ -*-
//
// QTildeMatching.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 QTildeMatching class.
//
#include "QTildeMatching.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
using namespace Herwig;
QTildeMatching::QTildeMatching()
: theCorrectForXZMismatch(true) {}
QTildeMatching::~QTildeMatching() {}
IBPtr QTildeMatching::clone() const {
return new_ptr(*this);
}
IBPtr QTildeMatching::fullclone() const {
return new_ptr(*this);
}
void QTildeMatching::checkCutoff() {
if ( showerTildeKinematics() ) {
showerTildeKinematics()->
prepare(realCXComb(),bornCXComb());
showerTildeKinematics()->dipole(dipole());
showerTildeKinematics()->getShowerVariables();
}
}
void QTildeMatching::getShowerVariables() {
// already filled from checkCutoff in this case
if ( showerTildeKinematics() )
return;
// get the shower variables
calculateShowerVariables();
// check for the cutoff
dipole()->isAboveCutoff(isAboveCutoff());
// get the hard scale
dipole()->showerHardScale(hardScale());
// check for phase space
dipole()->isInShowerPhasespace(isInShowerPhasespace());
}
bool QTildeMatching::isInShowerPhasespace() const {
assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) &&
"implementation only provided for default and pt cutoff");
Energy qtildeHard = ZERO;
Energy qtilde = dipole()->showerScale();
assert(!dipole()->showerParameters().empty());
double z = dipole()->showerParameters()[0];
// FF
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
qtildeHard =
theQTildeFinder->
calculateFinalFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
bornCXComb()->meMomenta()[dipole()->bornSpectator()],
bornCXComb()->mePartonData()[dipole()->bornEmitter()]->iColour() == PDT::Colour3).first;
}
// FI
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
qtildeHard =
theQTildeFinder->
calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornSpectator()],
bornCXComb()->meMomenta()[dipole()->bornEmitter()],false).second;
}
// IF
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
qtildeHard =
theQTildeFinder->
calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
bornCXComb()->meMomenta()[dipole()->bornSpectator()],false).first;
if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
return false;
}
// II
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
qtildeHard =
theQTildeFinder->
calculateInitialInitialScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
bornCXComb()->meMomenta()[dipole()->bornSpectator()]).first;
if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
return false;
}
Energy Qg = theQTildeSudakov->kinScale();
Energy2 pt2 = ZERO;
if ( dipole()->bornEmitter() > 1 ) {
Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu);
else
pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg);
}
if ( dipole()->bornEmitter() < 2 ) {
pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg);
}
if ( pt2 < theQTildeSudakov->pT2min() )
return false;
bool hardVeto = restrictPhasespace() && sqrt(pt2) >= dipole()->showerHardScale();
return qtilde <= qtildeHard && !hardVeto;
}
bool QTildeMatching::isAboveCutoff() const {
assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) &&
"implementation only provided for default and pt cutoff");
Energy qtilde = dipole()->showerScale();
assert(!dipole()->showerParameters().empty());
double z = dipole()->showerParameters()[0];
Energy Qg = theQTildeSudakov->kinScale();
if ( dipole()->bornEmitter() > 1 ) {
Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
return sqr(z*(1.-z)*qtilde) - sqr(mu) >= theQTildeSudakov->pT2min();
else
return sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg) >= theQTildeSudakov->pT2min();
}
if ( dipole()->bornEmitter() < 2 ) {
return
sqr((1.-z)*qtilde) - z*sqr(Qg) >= theQTildeSudakov->pT2min();
}
return false;
}
CrossSection QTildeMatching::dSigHatDR() const {
assert(!dipole()->showerParameters().empty());
pair<Energy2,double> vars =
make_pair(sqr(dipole()->showerScale()),
dipole()->showerParameters()[0]);
pair<int,int> ij(dipole()->bornEmitter(),
dipole()->bornSpectator());
double ccme2 =
dipole()->underlyingBornME()->largeNColourCorrelatedME2(ij,theLargeNBasis);
-
+
+ if(ccme2==0.)return 0.*nanobarn;
+
+ double lnme2=dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
+ if(lnme2==0){
+ generator()->log() <<"\nQTildeMatching: ";
+ generator()->log() <<"\n largeNME2 is ZERO, while largeNColourCorrelatedME2 is not ZERO." ;
+ generator()->log() <<"\n This is too seriuos.\n" ;
+ generator()->log() << Exception::runerror;
+ }
+
ccme2 *=
- dipole()->underlyingBornME()->me2() /
- dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
+ dipole()->underlyingBornME()->me2() /lnme2;
+
Energy2 prop = ZERO;
if ( dipole()->bornEmitter() > 1 ) {
prop =
(realCXComb()->meMomenta()[dipole()->realEmitter()] +
realCXComb()->meMomenta()[dipole()->realEmission()]).m2()
- bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2();
} else {
prop =
2.*vars.second*(realCXComb()->meMomenta()[dipole()->realEmitter()]*
realCXComb()->meMomenta()[dipole()->realEmission()]);
}
// note alphas included downstream from subtractionScaleWeight()
double xme2 = -8.*Constants::pi*ccme2*splitFn(vars)*realXComb()->lastSHat()/prop;
xme2 *=
pow(realCXComb()->lastSHat() / bornCXComb()->lastSHat(),
bornCXComb()->mePartonData().size()-4.);
double bornPDF = bornPDFWeight(dipole()->underlyingBornME()->lastScale());
if ( bornPDF == 0.0 )
return ZERO;
xme2 *= bornPDF;
xme2 *= dipole()->realEmissionME()->finalStateSymmetry() /
dipole()->underlyingBornME()->finalStateSymmetry();
// take care of mismatch between z and x as we are approaching the
// hard phase space boundary
// TODO get rid of this useless scale option business and simplify PDF handling in here
if ( dipole()->bornEmitter() < 2 && theCorrectForXZMismatch ) {
Energy2 emissionScale = ZERO;
if ( emissionScaleInSubtraction() == showerScale ) {
emissionScale = showerFactorizationScale();
} else if ( emissionScaleInSubtraction() == realScale ) {
emissionScale = dipole()->realEmissionME()->lastScale();
} else if ( emissionScaleInSubtraction() == bornScale ) {
emissionScale = dipole()->underlyingBornME()->lastScale();
}
double xzMismatch =
dipole()->subtractionParameters()[0] / dipole()->showerParameters()[0];
double realCorrectedPDF =
dipole()->bornEmitter() == 0 ?
dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,
xzMismatch) :
dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,
xzMismatch);
double realPDF =
dipole()->bornEmitter() == 0 ?
dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,1.0) :
dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,1.0);
if ( realPDF == 0.0 || realCorrectedPDF == 0.0 )
return ZERO;
xme2 *= realCorrectedPDF / realPDF;
}
Energy qtilde = sqrt(vars.first);
double z = vars.second;
Energy2 pt2 = ZERO;
Energy Qg = theQTildeSudakov->kinScale();
if ( dipole()->bornEmitter() > 1 ) {
Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu);
else
pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg);
}
if ( dipole()->bornEmitter() < 2 ) {
pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg);
}
assert(pt2 >= ZERO);
if ( profileScales() )
xme2 *= profileScales()->hardScaleProfile(dipole()->showerHardScale(),sqrt(pt2));
CrossSection res =
sqr(hbarc) *
realXComb()->jacobian() *
subtractionScaleWeight() *
xme2 /
(2. * realXComb()->lastSHat());
return res;
}
double QTildeMatching::me2() const {
throw Exception() << "QTildeMatching::me2(): Not intented to use. Disable the ShowerApproximationGenerator."
<< Exception::runerror;
return 0.;
}
void QTildeMatching::calculateShowerVariables() const {
Lorentz5Momentum n;
Energy2 Q2 = ZERO;
const Lorentz5Momentum& pb = bornCXComb()->meMomenta()[dipole()->bornEmitter()];
const Lorentz5Momentum& pc = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
if ( dipole()->bornEmitter() > 1 ) {
Q2 = (pb+pc).m2();
} else {
Q2 = -(pb-pc).m2();
}
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
double b = sqr(bornCXComb()->meMomenta()[dipole()->bornEmitter()].m())/Q2;
double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
double lambda = sqrt(1.+sqr(b)+sqr(c)-2.*b-2.*c-2.*b*c);
n = (1.-0.5*(1.-b+c-lambda))*pc - 0.5*(1.-b+c-lambda)*pb;
}
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
}
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
n = (1.+c)*pc - c*pb;
}
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
}
// the light-cone condition is numerically not very stable, so we
// explicitly push it on the light-cone here
n.setMass(ZERO);
n.rescaleEnergy();
double z = 0.0;
if ( dipole()->bornEmitter() > 1 ) {
z = 1. -
(n*realCXComb()->meMomenta()[dipole()->realEmission()])/
(n*bornCXComb()->meMomenta()[dipole()->bornEmitter()]);
} else {
z = 1. -
(n*realCXComb()->meMomenta()[dipole()->realEmission()])/
(n*realCXComb()->meMomenta()[dipole()->realEmitter()]);
}
Energy2 qtilde2 = ZERO;
Energy2 q2 = ZERO;
if ( dipole()->bornEmitter() > 1 ) {
q2 =
(realCXComb()->meMomenta()[dipole()->realEmitter()] + realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
qtilde2 = (q2 - bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(z*(1.-z));
} else {
q2 =
-(realCXComb()->meMomenta()[dipole()->realEmitter()] - realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
qtilde2 = (q2 + bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(1.-z);
}
assert(qtilde2 >= ZERO && z >= 0.0 && z <= 1.0);
dipole()->showerScale(sqrt(qtilde2));
dipole()->showerParameters().resize(1);
dipole()->showerParameters()[0] = z;
}
double QTildeMatching::splitFn(const pair<Energy2,double>& vars) const {
const Energy2& qtilde2 = vars.first;
const double& z = vars.second;
double Nc = SM().Nc();
// final state branching
if ( dipole()->bornEmitter() > 1 ) {
// final state quark quark branching
if ( abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 7 ) {
Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
return
((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z);
}
// final state gluon branching
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) {
if ( realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
// ATTENTION the factor 2 here is intentional as it cancels to the 1/2
// stemming from the large-N colour correlator
return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z));
}
if ( abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
Energy m = realCXComb()->mePartonData()[dipole()->realEmission()]->hardProcessMass();
return (1./2.)*(1.-2.*z*(1.-z)+2.*sqr(m)/(z*(1.-z)*qtilde2));
}
}
// final state squark branching
if ((abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 1000000 &&
abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 1000007) ||
(abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 2000000 &&
abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 2000007)){
Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
return ((sqr(Nc)-1.)/Nc)*(z-sqr(m)/(z*qtilde2))/(1.-z);
}
// final state gluino branching
if (bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == 1000021){
Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
return Nc*(1.+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z);
}
}
// initial state branching
if ( dipole()->bornEmitter() < 2 ) {
// g/g
if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
// see above for factor of 2
return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z));
}
// q/q
if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
return
((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z))/(1.-z);
}
// g/q
if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
return (1./2.)*(1.-2.*z*(1.-z));
}
// q/g
if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
return
((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(1.-z))/z;
}
}
return 0.0;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void QTildeMatching::doinit() {
assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
theShowerHandler->init();
theQTildeFinder->init();
theQTildeSudakov->init();
if ( theShowerHandler->scaleFactorOption() < 2 ) {
hardScaleFactor(theShowerHandler->hardScaleFactor());
factorizationScaleFactor(theShowerHandler->factorizationScaleFactor());
renormalizationScaleFactor(theShowerHandler->renormalizationScaleFactor());
}
profileScales(theShowerHandler->profileScales());
restrictPhasespace(theShowerHandler->restrictPhasespace());
hardScaleIsMuF(theShowerHandler->hardScaleIsMuF());
ShowerApproximation::doinit();
}
void QTildeMatching::doinitrun() {
assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
theShowerHandler->initrun();
theQTildeFinder->initrun();
theQTildeSudakov->initrun();
ShowerApproximation::doinitrun();
}
void QTildeMatching::persistentOutput(PersistentOStream & os) const {
os << theQTildeFinder << theQTildeSudakov
<< theShowerHandler << theCorrectForXZMismatch;
}
void QTildeMatching::persistentInput(PersistentIStream & is, int) {
is >> theQTildeFinder >> theQTildeSudakov
>> theShowerHandler >> theCorrectForXZMismatch;
}
// *** 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<QTildeMatching,Herwig::ShowerApproximation>
describeHerwigQTildeMatching("Herwig::QTildeMatching", "HwShower.so HwQTildeMatching.so");
void QTildeMatching::Init() {
static ClassDocumentation<QTildeMatching> documentation
("QTildeMatching implements NLO matching with the default shower.");
static Reference<QTildeMatching,QTildeFinder> interfaceQTildeFinder
("QTildeFinder",
"Set the partner finder to calculate hard scales.",
&QTildeMatching::theQTildeFinder, false, false, true, false, false);
static Reference<QTildeMatching,QTildeSudakov> interfaceQTildeSudakov
("QTildeSudakov",
"Set the partner finder to calculate hard scales.",
&QTildeMatching::theQTildeSudakov, false, false, true, false, false);
static Reference<QTildeMatching,ShowerHandler> interfaceShowerHandler
("ShowerHandler",
"",
&QTildeMatching::theShowerHandler, false, false, true, true, false);
static Switch<QTildeMatching,bool> interfaceCorrectForXZMismatch
("CorrectForXZMismatch",
"Correct for x/z mismatch near hard phase space boundary.",
&QTildeMatching::theCorrectForXZMismatch, true, false, false);
static SwitchOption interfaceCorrectForXZMismatchYes
(interfaceCorrectForXZMismatch,
"Yes",
"Include the correction factor.",
true);
static SwitchOption interfaceCorrectForXZMismatchNo
(interfaceCorrectForXZMismatch,
"No",
"Do not include the correction factor.",
false);
}
diff --git a/Shower/Couplings/ShowerAlphaQCD.cc b/Shower/Couplings/ShowerAlphaQCD.cc
--- a/Shower/Couplings/ShowerAlphaQCD.cc
+++ b/Shower/Couplings/ShowerAlphaQCD.cc
@@ -1,377 +1,372 @@
// -*- C++ -*-
//
// ShowerAlphaQCD.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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 ShowerAlphaQCD class.
//
#include "ShowerAlphaQCD.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Config/Constants.h"
using namespace Herwig;
DescribeClass<ShowerAlphaQCD,ShowerAlpha>
describeShowerAlphaQCD("Herwig::ShowerAlphaQCD","HwShower.so");
IBPtr ShowerAlphaQCD::clone() const {
return new_ptr(*this);
}
IBPtr ShowerAlphaQCD::fullclone() const {
return new_ptr(*this);
}
void ShowerAlphaQCD::persistentOutput(PersistentOStream & os) const {
os << _asType << _asMaxNP << ounit(_qmin,GeV) << _nloop << _lambdaopt << _thresopt
<< ounit(_lambdain,GeV) << _alphain << _inopt
<< _tolerance << _maxtry << _alphamin
<< ounit(_thresholds,GeV) << ounit(_lambda,GeV);
}
void ShowerAlphaQCD::persistentInput(PersistentIStream & is, int) {
is >> _asType >> _asMaxNP >> iunit(_qmin,GeV) >> _nloop >> _lambdaopt >> _thresopt
>> iunit(_lambdain,GeV) >> _alphain >> _inopt
>> _tolerance >> _maxtry >> _alphamin
>> iunit(_thresholds,GeV) >> iunit(_lambda,GeV);
}
void ShowerAlphaQCD::Init() {
static ClassDocumentation<ShowerAlphaQCD> documentation
("This (concrete) class describes the QCD alpha running.");
static Switch<ShowerAlphaQCD, int> intAsType
("NPAlphaS",
"Behaviour of AlphaS in the NP region",
&ShowerAlphaQCD::_asType, 1, false, false);
static SwitchOption intAsTypeZero
(intAsType, "Zero","zero below Q_min", 1);
static SwitchOption intAsTypeConst
(intAsType, "Const","const as(qmin) below Q_min", 2);
static SwitchOption intAsTypeLin
(intAsType, "Linear","growing linearly below Q_min", 3);
static SwitchOption intAsTypeQuad
(intAsType, "Quadratic","growing quadratically below Q_min", 4);
static SwitchOption intAsTypeExx1
(intAsType, "Exx1", "quadratic from AlphaMaxNP down to as(Q_min)", 5);
static SwitchOption intAsTypeExx2
(intAsType, "Exx2", "const = AlphaMaxNP below Q_min", 6);
// default such that as(qmin) = 1 in the current parametrization.
// min = Lambda3
static Parameter<ShowerAlphaQCD,Energy> intQmin
("Qmin", "Q < Qmin is treated with NP parametrization as of (unit [GeV])",
&ShowerAlphaQCD::_qmin, GeV, 0.630882*GeV, 0.330445*GeV,
100.0*GeV,false,false,false);
static Parameter<ShowerAlphaQCD,double> interfaceAlphaMaxNP
("AlphaMaxNP",
"Max value of alpha in NP region, only relevant if NPAlphaS = 5,6",
&ShowerAlphaQCD::_asMaxNP, 1.0, 0., 100.0,
false, false, Interface::limited);
static Parameter<ShowerAlphaQCD,unsigned int> interfaceNumberOfLoops
("NumberOfLoops",
"The number of loops to use in the alpha_S calculation",
&ShowerAlphaQCD::_nloop, 3, 1, 3,
false, false, Interface::limited);
static Switch<ShowerAlphaQCD,bool> interfaceLambdaOption
("LambdaOption",
"Option for the calculation of the Lambda used in the simulation from the input"
" Lambda_MSbar",
&ShowerAlphaQCD::_lambdaopt, false, false, false);
static SwitchOption interfaceLambdaOptionfalse
(interfaceLambdaOption,
"Same",
"Use the same value",
false);
static SwitchOption interfaceLambdaOptionConvert
(interfaceLambdaOption,
"Convert",
"Use the conversion to the Herwig scheme from NPB349, 635",
true);
static Parameter<ShowerAlphaQCD,Energy> interfaceLambdaQCD
("LambdaQCD",
"Input value of Lambda_MSBar",
&ShowerAlphaQCD::_lambdain, MeV, 0.208364*GeV, 100.0*MeV, 500.0*MeV,
false, false, Interface::limited);
static Parameter<ShowerAlphaQCD,double> interfaceAlphaMZ
("AlphaMZ",
"The input value of the strong coupling at the Z mass ",
&ShowerAlphaQCD::_alphain, 0.118, 0.1, 0.2,
false, false, Interface::limited);
static Switch<ShowerAlphaQCD,bool> interfaceInputOption
("InputOption",
"Option for inputing the initial value of the coupling",
&ShowerAlphaQCD::_inopt, true, false, false);
static SwitchOption interfaceInputOptionAlphaMZ
(interfaceInputOption,
"AlphaMZ",
"Use the value of alpha at MZ to calculate the coupling",
true);
static SwitchOption interfaceInputOptionLambdaQCD
(interfaceInputOption,
"LambdaQCD",
"Use the input value of Lambda to calculate the coupling",
false);
static Parameter<ShowerAlphaQCD,double> interfaceTolerance
("Tolerance",
"The tolerance for discontinuities in alphaS at thresholds.",
&ShowerAlphaQCD::_tolerance, 1e-10, 1e-20, 1e-4,
false, false, Interface::limited);
static Parameter<ShowerAlphaQCD,unsigned int> interfaceMaximumIterations
("MaximumIterations",
"The maximum number of iterations for the Newton-Raphson method to converge.",
&ShowerAlphaQCD::_maxtry, 100, 10, 1000,
false, false, Interface::limited);
static Switch<ShowerAlphaQCD,bool> interfaceThresholdOption
("ThresholdOption",
"Whether to use the consistuent or normal masses for the thresholds",
&ShowerAlphaQCD::_thresopt, true, false, false);
static SwitchOption interfaceThresholdOptionCurrent
(interfaceThresholdOption,
"Current",
"Use the current masses",
true);
static SwitchOption interfaceThresholdOptionConstituent
(interfaceThresholdOption,
"Constituent",
"Use the constitent masses.",
false);
static Command<ShowerAlphaQCD> interfaceValue
("Value",
"",
&ShowerAlphaQCD::value, false);
}
void ShowerAlphaQCD::doinit() {
ShowerAlpha::doinit();
// calculate the value of 5-flavour lambda
// evaluate the initial
// value of Lambda from alphas if needed using Newton-Raphson
- if(_inopt)
- {_lambda[2]=computeLambda(getParticleData(ParticleID::Z0)->mass(),_alphain,5);}
+ if(_inopt) {
+ _lambda[2]=computeLambda(getParticleData(ParticleID::Z0)->mass(),_alphain,5);
+ }
// otherwise it was an input parameter
else{_lambda[2]=_lambdain;}
// convert lambda to the Monte Carlo scheme if needed
using Constants::pi;
- if(_lambdaopt){_lambda[2] *=exp(0.5*(67.-3.*sqr(pi)-50./3.)/23.)/sqrt(2.);}
+ if(_lambdaopt) _lambda[2] *=exp(0.5*(67.-3.*sqr(pi)-50./3.)/23.)/sqrt(2.);
// compute the threshold matching
// top threshold
for(int ix=1;ix<4;++ix) {
if(_thresopt)
_thresholds[ix]=getParticleData(ix+3)->mass();
else
_thresholds[ix]=getParticleData(ix+3)->constituentMass();
}
// compute 6 flavour lambda by matching at top mass using Newton Raphson
_lambda[3]=computeLambda(_thresholds[3],alphaS(_thresholds[3],_lambda[2],5),6);
// bottom threshold
// compute 4 flavour lambda by matching at bottom mass using Newton Raphson
_lambda[1]=computeLambda(_thresholds[2],alphaS(_thresholds[2],_lambda[2],5),4);
// charm threshold
// compute 3 flavour lambda by matching at charm mass using Newton Raphson
_lambda[0]=computeLambda(_thresholds[1],alphaS(_thresholds[1],_lambda[1],4),3);
// final threshold is qmin
_thresholds[0]=_qmin;
// compute the maximum value of as
if ( _asType < 5 ) _alphamin = value(sqr(_qmin)+1.0e-8*sqr(MeV)); // approx as = 1
else _alphamin = max(_asMaxNP, value(sqr(_qmin)+1.0e-8*sqr(MeV)));
// check consistency lambda_3 < qmin
if(_lambda[0]>_qmin)
Throw<InitException>() << "The value of Qmin is less than Lambda_3 in"
<< " ShowerAlphaQCD::doinit " << Exception::abortnow;
-
}
double ShowerAlphaQCD::value(const Energy2 scale) const {
pair<short,Energy> nflam;
Energy q = sqrt(scale);
double val(0.);
// special handling if the scale is less than Qmin
if (q < _qmin) {
nflam = getLamNfTwoLoop(_qmin);
double val0 = alphaS(_qmin, nflam.second, nflam.first);
switch (_asType) {
case 1:
// flat, zero; the default type with no NP effects.
val = 0.;
break;
case 2:
// flat, non-zero alpha_s = alpha_s(q2min).
val = val0;
break;
case 3:
// linear in q
val = val0*q/_qmin;
break;
case 4:
// quadratic in q
val = val0*sqr(q/_qmin);
break;
case 5:
// quadratic in q, starting off at asMaxNP, ending on as(qmin)
val = (val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP;
break;
case 6:
// just asMaxNP and constant
val = _asMaxNP;
break;
}
} else {
// the 'ordinary' case
nflam = getLamNfTwoLoop(q);
val = alphaS(q, nflam.second, nflam.first);
}
return scaleFactor() * val;
}
double ShowerAlphaQCD::overestimateValue() const {
return scaleFactor() * _alphamin;
}
double ShowerAlphaQCD::ratio(const Energy2 scale) const {
pair<short,Energy> nflam;
Energy q = sqrt(scale);
double val(0.);
// special handling if the scale is less than Qmin
if (q < _qmin) {
nflam = getLamNfTwoLoop(_qmin);
double val0 = alphaS(_qmin, nflam.second, nflam.first);
switch (_asType) {
case 1:
// flat, zero; the default type with no NP effects.
val = 0.;
break;
case 2:
// flat, non-zero alpha_s = alpha_s(q2min).
val = val0;
break;
case 3:
// linear in q
val = val0*q/_qmin;
break;
case 4:
// quadratic in q
val = val0*sqr(q/_qmin);
break;
case 5:
// quadratic in q, starting off at asMaxNP, ending on as(qmin)
val = (val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP;
break;
case 6:
// just asMaxNP and constant
val = _asMaxNP;
break;
}
} else {
// the 'ordinary' case
nflam = getLamNfTwoLoop(q);
val = alphaS(q, nflam.second, nflam.first);
}
// denominator
return val/_alphamin;
}
string ShowerAlphaQCD::value (string scale) {
istringstream readscale(scale);
double inScale; readscale >> inScale;
Energy theScale = inScale * GeV;
initialize();
ostringstream showvalue ("");
showvalue << "alpha_s (" << theScale/GeV << " GeV) = "
<< value (sqr(theScale));
return showvalue.str();
}
pair<short, Energy> ShowerAlphaQCD::getLamNfTwoLoop(Energy q) const {
short nf = 6;
// get lambda and nf according to the thresholds
if (q < _thresholds[1]) nf = 3;
else if (q < _thresholds[2]) nf = 4;
else if (q < _thresholds[3]) nf = 5;
return pair<short,Energy>(nf, _lambda[nf-3]);
}
Energy ShowerAlphaQCD::computeLambda(Energy match,
double alpha,
unsigned int nflav) const {
Energy lamtest=200.0*MeV;
double xtest;
unsigned int ntry=0;
- do
- {
- ++ntry;
- xtest=log(sqr(match/lamtest));
- xtest+= (alpha-alphaS(match,lamtest,nflav))/derivativealphaS(match,lamtest,nflav);
- lamtest=match/exp(0.5*xtest);
- }
+ do {
+ ++ntry;
+ xtest = log(sqr(match/lamtest));
+ xtest += (alpha-alphaS(match,lamtest,nflav))/derivativealphaS(match,lamtest,nflav);
+ Energy newLambda = match/exp(0.5*xtest);
+ lamtest = newLambda<match ? newLambda : 0.5*(lamtest+match);
+ }
while(abs(alpha-alphaS(match,lamtest,nflav)) > _tolerance && ntry < _maxtry);
return lamtest;
}
-double ShowerAlphaQCD::derivativealphaS(Energy q, Energy lam, int nf) const
-{
+double ShowerAlphaQCD::derivativealphaS(Energy q, Energy lam, int nf) const {
using Constants::pi;
double lx = log(sqr(q/lam));
double b0 = 11. - 2./3.*nf;
double b1 = 51. - 19./3.*nf;
double b2 = 2857. - 5033./9.*nf + 325./27.*sqr(nf);
if(_nloop==1)
return -4.*pi/(b0*sqr(lx));
else if(_nloop==2)
- return -4.*pi/(b0*sqr(lx))*(1.-2.*b1/sqr(b0)/lx*(1.-2.*log(lx)));
+ return -4.*pi/(b0*sqr(lx))*(1.+2.*b1/sqr(b0)/lx*(1.-2.*log(lx)));
else
- return -4.*pi/(b0*sqr(lx))*(1.
- - 2.*b1/sqr(b0)/lx*(1.-2.*log(lx))
- + 4.*sqr(b1)/(sqr(sqr(b0))*sqr(lx))*
- (1.
- - 2.*log(lx)
- + 3.*(sqr(log(lx) - 0.5)+b2*b0/(8.*sqr(b1))-1.25)
- )
- );
+ return -4.*pi/(b0*sqr(lx))*
+ (1. + 2.*b1/sqr(b0)/lx*(1.-2.*log(lx))
+ + 4.*sqr(b1)/(sqr(sqr(b0))*sqr(lx))*(1. - 2.*log(lx)
+ + 3.*(sqr(log(lx) - 0.5)+b2*b0/(8.*sqr(b1))-1.25)));
}
double ShowerAlphaQCD::alphaS(Energy q, Energy lam, int nf) const {
using Constants::pi;
double lx(log(sqr(q/lam)));
double b0 = 11. - 2./3.*nf;
double b1 = 51. - 19./3.*nf;
double b2 = 2857. - 5033./9.*nf + 325./27.*sqr(nf);
// one loop
if(_nloop==1)
{return 4.*pi/(b0*lx);}
// two loop
else if(_nloop==2) {
return 4.*pi/(b0*lx)*(1.-2.*b1/sqr(b0)*log(lx)/lx);
}
// three loop
else
{return 4.*pi/(b0*lx)*(1.-2.*b1/sqr(b0)*log(lx)/lx +
4.*sqr(b1)/(sqr(sqr(b0))*sqr(lx))*
(sqr(log(lx) - 0.5) + b2*b0/(8.*sqr(b1)) - 5./4.));}
}
diff --git a/Tests/Rivet/Templates/Hadron-Dipole.in b/Tests/Rivet/Templates/Hadron-Dipole.in
--- a/Tests/Rivet/Templates/Hadron-Dipole.in
+++ b/Tests/Rivet/Templates/Hadron-Dipole.in
@@ -1,42 +1,55 @@
##################################################
## Collider type
##################################################
read Matchbox/PPCollider.in
##################################################
## Matrix element library selection
##################################################
read Matchbox/MadGraph-OpenLoops.in
##################################################
## Matching and shower selection
## Please also see flavour scheme settings
## towards the end of the input file.
##################################################
read Matchbox/MCatNLO-DipoleShower.in
##################################################
## PDF choice
##################################################
read Matchbox/FiveFlavourScheme.in
## required for dipole shower in five flavour scheme
# read Matchbox/FiveFlavourNoBMassScheme.in
read Matchbox/MMHT2014.in
##################################################
# Create the Herwig analysis
##################################################
create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
set /Herwig/Analysis/RivetAnalysis:Debug No
##################################################
## Save the generator
##################################################
do /Herwig/MatrixElements/Matchbox/Factory:ProductionMode
## Select the process
cd /Herwig/MatrixElements/Matchbox
+
+do Factory:StartParticleGroup pnob
+insert Factory:ParticleGroup 0 /Herwig/Particles/c
+insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/s
+insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/d
+insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/u
+insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
+insert Factory:ParticleGroup 0 /Herwig/Particles/g
+do Factory:EndParticleGroup
+
${process}
read ${parameterFile}
cd /Herwig/Generators
saverun ${runname} EventGenerator
diff --git a/Tests/python/make_input_files.py b/Tests/python/make_input_files.py
--- a/Tests/python/make_input_files.py
+++ b/Tests/python/make_input_files.py
@@ -1,494 +1,495 @@
#! /usr/bin/env python
import logging,sys,os
from string import strip, Template
import sys
if sys.version_info[:3] < (2,4,0):
print "rivet scripts require Python version >= 2.4.0... exiting"
sys.exit(1)
if __name__ == "__main__":
import logging
from optparse import OptionParser, OptionGroup
parser = OptionParser(usage="%prog name [...]")
(opts, args) = parser.parse_args()
## Check args
if len(args) != 1:
logging.error("Must specify at least input file")
sys.exit(1)
name=args[0]
collider=""
# select the template to load
# collider
if(name.find("BFactory")==0) :
collider="BFactory"
elif(name.find("LEP")==0) :
collider="LEP"
elif(name.find("DIS")==0) :
collider="DIS"
elif(name.find("TVT")==0) :
collider="TVT"
elif(name.find("LHC")==0) :
collider="LHC"
simulation=""
if(name.find("Matchbox")>0) :
simulation="Matchbox"
elif(name.find("Dipole")>0) :
simulation="Dipole"
elif(name.find("Powheg")>0) :
simulation="Powheg"
if(collider=="") :
logging.error("Can\'t find collider")
sys.exit(1)
# find the template
if(simulation=="") :
if(collider.find("TVT")>=0 or collider.find("LHC")>=0) :
templateName="Hadron.in"
elif(collider.find("BFactory")<0) :
templateName= "%s.in" % (collider)
else :
templateName= "LEP.in"
else :
if(collider.find("TVT")>=0 or collider.find("LHC")>=0) :
templateName= "Hadron-%s.in" % (simulation)
elif(collider.find("BFactory")<0) :
templateName= "%s-%s.in" % (collider,simulation)
else :
templateName= "LEP-%s.in" % (simulation)
with open(os.path.join("Rivet/Templates",templateName), 'r') as f:
templateText = f.read()
template = Template( templateText )
# work out the name of the parameter file
if(simulation=="") :
istart = 1
else :
istart = 2
nameSplit=name.split("-")
parameterName=nameSplit[istart]
for i in range(istart+1,len(nameSplit)) :
parameterName += "-%s" % nameSplit[i]
# work out the process and parameters
process=""
bscheme="read Matchbox/FiveFlavourScheme.in"
# Bfactory
if(collider=="BFactory") :
if(simulation=="") :
process = "set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4"
if(parameterName=="10.58") :
process += "\ncreate Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so\nset /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)\nset /Herwig/MatrixElements/MEUpsilon:Coupling 0.0004151809\ninsert /Herwig/MatrixElements/SimpleEE:MatrixElements 0 /Herwig/MatrixElements/MEUpsilon"
elif(simulation=="Powheg") :
process = "set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 4"
if(parameterName=="10.58") :
process += "\ncreate Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so\nset /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)\nset /Herwig/MatrixElements/MEUpsilon:Coupling 0.0004151809\ninsert /Herwig/MatrixElements/SimpleEE:MatrixElements 0 /Herwig/MatrixElements/MEUpsilon"
elif(simulation=="Matchbox" or simulation == "Dipole" ) :
process = "do Factory:Process e- e+ -> u ubar\ndo Factory:Process e- e+ -> d dbar\ndo Factory:Process e- e+ -> c cbar\ndo Factory:Process e- e+ -> s sbar"
if(parameterName=="10.58") :
process += "\ninsert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers 0 /Herwig/MatrixElements/SimpleEE\ncreate Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so\nset /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)\nset /Herwig/MatrixElements/MEUpsilon:Coupling 0.0004151809\ninsert /Herwig/MatrixElements/SimpleEE:MatrixElements 0 /Herwig/MatrixElements/MEUpsilon\n"
# DIS
elif(collider=="DIS") :
if(simulation=="") :
if(parameterName.find("NoME")>=0) :
process = "set /Herwig/Shower/Evolver:MECorrMode 0"
parameterName=parameterName.replace("NoME-","")
else :
process = ""
elif(simulation=="Powheg") :
process = ""
elif(simulation=="Matchbox" or simulation == "Dipole" ) :
if(parameterName.find("e-")>=0) :
process="do Factory:Process e- p -> e- j"
else :
process="do Factory:Process e+ p -> e+ j"
# LEP
elif(collider=="LEP") :
+ print parameterName
if(simulation=="") :
process=""
- if(parameterName=="-10") :
+ if(parameterName=="10") :
process="set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4"
elif(simulation=="Powheg") :
process=""
- if(parameterName=="-10") :
+ if(parameterName=="10") :
process="set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 4"
elif(simulation=="Matchbox" or simulation == "Dipole" ) :
- if(parameterName=="-10") :
+ if(parameterName=="10") :
process="do Factory:Process e- e+ -> u ubar\ndo Factory:Process e- e+ -> d dbar\ndo Factory:Process e- e+ -> c cbar\ndo Factory:Process e- e+ -> s sbar"
else :
process="do Factory:Process e- e+ -> j j"
# TVT
elif(collider=="TVT") :
process="set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/pbar-\n"
if(parameterName.find("Run-II")>=0) :
process+="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 1960.0\n"
elif(parameterName.find("Run-I")>=0) :
process+="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 1800.0\n"
if(simulation=="") :
if(parameterName.find("Run-I-WZ")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\ninsert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
elif(parameterName.find("Run-I-W")>=0 or parameterName.find("Run-II-W")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\n"
elif(parameterName.find("Run-I-Z")>=0 or parameterName.find("Run-II-Z-e")>=0) :
process +="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
elif(parameterName.find("Run-II-Z-mu")>=0) :
process +="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
elif(simulation=="Powheg") :
if(parameterName.find("Run-I-WZ")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\ninsert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
elif(parameterName.find("Run-I-W")>=0 or parameterName.find("Run-II-W")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\n"
elif(parameterName.find("Run-I-Z")>=0 or parameterName.find("Run-II-Z-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
elif(parameterName.find("Run-II-Z-mu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Muon\n"
elif(simulation=="Matchbox" or simulation == "Dipole" ) :
if(parameterName.find("Run-I-WZ")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar e+ e-\ndo Factory:Process p pbar e+ nu\ndo Factory:Process p pbar e- nu\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV\n"
elif(parameterName.find("Run-I-W")>=0 or parameterName.find("Run-II-W")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar e+ nu\ndo Factory:Process p pbar e- nu\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
elif(parameterName.find("Run-I-Z")>=0 or parameterName.find("Run-II-Z-e")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV\n"
elif(parameterName.find("Run-II-Z-mu")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar mu+ mu-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV\n"
# LHC
elif(collider=="LHC") :
if(parameterName.find("-7-")>=0) :
process="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 7000.0\n"
elif(parameterName.find("-8-")>=0) :
process="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 8000.0\n"
else :
process="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 7000.0\n"
if(simulation=="") :
if(parameterName.find("WZ")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process WZ\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+;\n"
process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_ebar,e-; /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("WW-emu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process WW\n"
process+="set /Herwig/Particles/W+:Synchronized 0\n"
process+="set /Herwig/Particles/W-:Synchronized 0\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+;\n"
process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("WW-ll")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process WW\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+; /Herwig/Particles/W+/W+->nu_tau,tau+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("ZZ-ll")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process ZZ\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("ZZ-lv")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process ZZ\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+; /Herwig/Particles/Z0/Z0->nu_e,nu_ebar; /Herwig/Particles/Z0/Z0->nu_mu,nu_mubar; /Herwig/Particles/Z0/Z0->nu_tau,nu_taubar;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("LHC-W-Z-e")>0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\n"
elif(parameterName.find("LHC-W-Z-mu")>0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Muon\n"
elif(parameterName.find("W-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\n"
elif(parameterName.find("W-mu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Muon\n"
elif(parameterName.find("Z-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
elif(parameterName.find("Z-mu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
elif(parameterName.find("Z-LowMass-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.*GeV\nset /Herwig/Cuts/MassCut:MinM 20.*GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
elif(parameterName.find("Z-MedMass-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
process+="set /Herwig/Cuts/QCDCuts:MHatMin 40.*GeV\nset /Herwig/Cuts/MassCut:MinM 40.*GeV\nset /Herwig/Cuts/MassCut:MaxM 130.*GeV\n"
elif(parameterName.find("Z-LowMass-mu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.*GeV\nset /Herwig/Cuts/MassCut:MinM 10.*GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
elif(parameterName.find("W-Jet")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEWJet\nset MEWJet:WDecay Electron\n"
if(parameterName.find("W-Jet-1-e")>=0) :
process+="set /Herwig/Cuts/WBosonKtCut:MinKT 100.0*GeV\n"
parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
elif(parameterName.find("W-Jet-2-e")>=0) :
process+="set /Herwig/Cuts/WBosonKtCut:MinKT 190.0*GeV\n"
parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
elif(parameterName.find("W-Jet-3-e")>=0) :
process+="set /Herwig/Cuts/WBosonKtCut:MinKT 270.0*GeV\n"
parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
elif(parameterName.find("Z-Jet")>=0) :
process+="insert SimpleQCD:MatrixElements[0] MEZJet\nset MEZJet:ZDecay Electron\n"
if(parameterName.find("Z-Jet-1-e")>=0) :
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 100.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
elif(parameterName.find("Z-Jet-2-e")>=0) :
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 190.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
elif(parameterName.find("Z-Jet-3-e")>=0) :
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 270.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
else :
logging.error(" Process %s not supported for internal matrix elements" % name)
sys.exit(1)
elif(simulation=="Powheg") :
if(parameterName.find("WZ")>=0) :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process WZ\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+;\n"
process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_ebar,e-; /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("WW-emu")>=0) :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process WW\n"
process+="set /Herwig/Particles/W+:Synchronized 0\n"
process+="set /Herwig/Particles/W-:Synchronized 0\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+;\n"
process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("WW-ll")>=0) :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process WW\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+; /Herwig/Particles/W+/W+->nu_tau,tau+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("ZZ-ll")>=0) :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process ZZ\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("ZZ-lv")>=0) :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process ZZ\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+; /Herwig/Particles/Z0/Z0->nu_e,nu_ebar; /Herwig/Particles/Z0/Z0->nu_mu,nu_mubar; /Herwig/Particles/Z0/Z0->nu_tau,nu_taubar;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("LHC-W-Z-e")>0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\n"
elif(parameterName.find("LHC-W-Z-mu")>0) :
process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Muon\n"
elif(parameterName.find("W-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\n"
elif(parameterName.find("W-mu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Muon\n"
elif(parameterName.find("Z-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
elif(parameterName.find("Z-mu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Muon\n"
elif(parameterName.find("Z-LowMass-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.*GeV\nset /Herwig/Cuts/MassCut:MinM 20.*GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
elif(parameterName.find("Z-MedMass-e")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
process+="set /Herwig/Cuts/QCDCuts:MHatMin 40.*GeV\nset /Herwig/Cuts/MassCut:MinM 40.*GeV\nset /Herwig/Cuts/MassCut:MaxM 130.*GeV\n"
elif(parameterName.find("Z-LowMass-mu")>=0) :
process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Muon\n"
process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.*GeV\nset /Herwig/Cuts/MassCut:MinM 10.*GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
else :
logging.error(" Process %s not supported for internal POWHEG matrix elements" % name)
sys.exit(1)
elif(simulation=="Matchbox" or simulation == "Dipole" ) :
if(parameterName.find("WZ")>=0) :
process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p W+ Z0\ndo Factory:Process p p W- Z0\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 171.6*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n\n"
process+="set /Herwig/Cuts/MassCut:MinM 66*GeV\nset /Herwig/Cuts/MassCut:MaxM 116*GeV\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+;\n"
process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_ebar,e-; /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("WW-emu")>=0) :
process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process pnob pnob W+ W-\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="set /Herwig/Particles/W+:Synchronized 0\n"
process+="set /Herwig/Particles/W-:Synchronized 0\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+;\n"
process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("WW-ll")>=0) :
process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process pnob pnob W+ W-\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+; /Herwig/Particles/W+/W+->nu_tau,tau+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("ZZ-ll")>=0) :
process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p Z0 Z0\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("ZZ-lv")>=0) :
process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p Z0 Z0\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+; /Herwig/Particles/Z0/Z0->nu_e,nu_ebar; /Herwig/Particles/Z0/Z0->nu_mu,nu_mubar; /Herwig/Particles/Z0/Z0->nu_tau,nu_taubar;\n"
process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
elif(parameterName.find("LHC-W-Z-e")>0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\n"
process+="do Factory:Process p p e+ e-\ndo Factory:Process p p e+ nu\ndo Factory:Process p p e- nu\n"
process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV\n"
elif(parameterName.find("LHC-W-Z-mu")>0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\n"
process+="do Factory:Process p p mu+ mu-\ndo Factory:Process p p mu+ nu\ndo Factory:Process p p mu- nu\n"
process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV\n"
elif(parameterName.find("W-e")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ nu\ndo Factory:Process p p e- nu\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
elif(parameterName.find("W-mu")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p mu+ nu\ndo Factory:Process p p mu- nu\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
elif(parameterName.find("Z-e")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
elif(parameterName.find("Z-mu")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p mu+ mu-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV\n"
elif(parameterName.find("Z-LowMass-e")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 20*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 70*GeV\n"
elif(parameterName.find("Z-MedMass-e")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 40*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 130*GeV\n"
elif(parameterName.find("Z-LowMass-mu")>=0) :
process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p mu+ mu-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 20*GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 70*GeV\n"
elif(parameterName.find("W-Jet")>=0) :
process+="set Factory:OrderInAlphaS 1\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ nu j\ndo Factory:Process p p e- nu j\n\n"
process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/HTScale\n"
process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
if(parameterName.find("W-Jet-1-e")>=0) :
process+="set /Herwig/Cuts/FirstJet:PtMin 100.*GeV\n"
parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
elif(parameterName.find("W-Jet-2-e")>=0) :
process+="set /Herwig/Cuts/FirstJet:PtMin 190.0*GeV\n"
parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
elif(parameterName.find("W-Jet-3-e")>=0) :
process+="set /Herwig/Cuts/FirstJet:PtMin 270.0*GeV\n"
parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
elif(parameterName.find("Z-Jet")>=0) :
process+="set Factory:OrderInAlphaS 1\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e- j\n"
process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/HTScale\n"
process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
if(parameterName.find("Z-Jet-1-e")>=0) :
process+="set /Herwig/Cuts/FirstJet:PtMin 100.*GeV\n"
parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
elif(parameterName.find("Z-Jet-2-e")>=0) :
process+="set /Herwig/Cuts/FirstJet:PtMin 190.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
elif(parameterName.find("Z-Jet-3-e")>=0) :
process+="set /Herwig/Cuts/FirstJet:PtMin 270.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
elif(parameterName.find("Z-bb")>=0) :
bscheme="read Matchbox/FourFlavourScheme.in"
process+="set /Herwig/Particles/b:HardProcessMass 4.2*GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2*GeV\n"
process+="set Factory:OrderInAlphaS 2\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e- b bbar\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 91.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="set /Herwig/Cuts/MassCut:MinM 66*GeV\nset /Herwig/Cuts/MassCut:MaxM 116*GeV\n"
process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/SecondJet\n"
process+="set /Herwig/Cuts/FirstJet:PtMin 18.*GeV\n"
process+="set /Herwig/Cuts/SecondJet:PtMin 15.*GeV\n"
elif(parameterName.find("Z-b")>=0) :
process+="set Factory:OrderInAlphaS 1\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e- b\ndo Factory:Process p p e+ e- bbar\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 91.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="set /Herwig/Cuts/MassCut:MinM 66*GeV\nset /Herwig/Cuts/MassCut:MaxM 116*GeV\n"
process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
process+="set /Herwig/Cuts/FirstJet:PtMin 15.*GeV\n"
elif(parameterName.find("W-b")>=0) :
bscheme="read Matchbox/FourFlavourScheme.in"
process += "set /Herwig/Particles/b:HardProcessMass 4.2*GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2*GeV\n"
process += "set Factory:OrderInAlphaS 2\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ nu b bbar\ndo Factory:Process p p e- nu b bbar\n"
process += "set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 80.4*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
process+="set /Herwig/Cuts/FirstJet:PtMin 30.*GeV\n"
else :
logging.error(" Process %s not supported for Matchbox matrix elements" % name)
sys.exit(1)
# write the file
if(simulation=="Matchbox" or simulation =="Dipole") :
with open(os.path.join("Rivet",name+".in") ,'w') as f:
f.write( template.substitute({ 'process' : process,
'runname' : name,
'bscheme' : bscheme,
'parameterFile' : os.path.join(collider,collider+"-"+parameterName+".in") }))
else :
with open(os.path.join("Rivet",name+".in") ,'w') as f:
f.write( template.substitute({ 'process' : process,
'runname' : name,
'parameterFile' : os.path.join(collider,collider+"-"+parameterName+".in") }))

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:21 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805128
Default Alt Text
(114 KB)

Event Timeline