Page MenuHomeHEPForge
Files F8310135

No OneTemporary
Actions

View Options

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,525 +1,525 @@
// -*- C++ -*-
//
// QTildeMatching.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// 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 < max(theQTildeSudakov->pT2min(),sqr(safeCut()) ))
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) >=
max(theQTildeSudakov->pT2min(),sqr(safeCut()));
else
return sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg) >=
max(theQTildeSudakov->pT2min(),sqr(safeCut()));
}
if ( dipole()->bornEmitter() < 2 ) {
return
sqr((1.-z)*qtilde) - z*sqr(Qg) >=
max(theQTildeSudakov->pT2min(),sqr(safeCut()));
}
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() /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()]);
}
// allow small violations (numerical inaccuracies)
if ( z <= 0 && z >= -1e-6 ) {
z = std::numeric_limits<double>::epsilon();
} else if ( z >= 1 && z <= 1+1e-6 ) {
z = 1-std::numeric_limits<double>::epsilon();
}
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);
}
if ( qtilde2 < ZERO ) {
qtilde2 = ZERO;
}
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();
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
+ static Reference<QTildeMatching,PartnerFinder> interfaceQTildeFinder
("QTildeFinder",
"Set the partner finder to calculate hard scales.",
&QTildeMatching::theQTildeFinder, false, false, true, false, false);
interfaceQTildeFinder.rank(-1);
static Reference<QTildeMatching,QTildeSudakov> interfaceQTildeSudakov
("QTildeSudakov",
"Set the partner finder to calculate hard scales.",
&QTildeMatching::theQTildeSudakov, false, false, true, false, false);
interfaceQTildeSudakov.rank(-1);
static Reference<QTildeMatching,ShowerHandler> interfaceShowerHandler
("ShowerHandler",
"The QTilde shower handler to use.",
&QTildeMatching::theShowerHandler, false, false, true, true, false);
interfaceShowerHandler.rank(-1);
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);
interfaceCorrectForXZMismatch.rank(-1);
}
diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.h b/MatrixElement/Matchbox/Matching/QTildeMatching.h
--- a/MatrixElement/Matchbox/Matching/QTildeMatching.h
+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.h
@@ -1,201 +1,201 @@
// -*- C++ -*-
//
// QTildeMatching.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_QTildeMatching_H
#define Herwig_QTildeMatching_H
//
// This is the declaration of the QTildeMatching class.
//
#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
#include "Herwig/Shower/ShowerHandler.h"
-#include "Herwig/Shower/QTilde/Default/QTildeFinder.h"
+#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
#include "Herwig/Shower/QTilde/Default/QTildeSudakov.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief QTildeMatching implements NLO matching with the default shower.
*
*/
class QTildeMatching: public Herwig::ShowerApproximation {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
QTildeMatching();
/**
* The destructor.
*/
virtual ~QTildeMatching();
//@}
public:
/**
* Return the shower approximation to the real emission cross
* section for the given pair of Born and real emission
* configurations.
*/
virtual CrossSection dSigHatDR() const;
/**
* Return the shower approximation splitting kernel for the given
* pair of Born and real emission configurations in units of the
* Born center of mass energy squared, and including a weight to
* project onto the splitting given by the dipole used.
*/
virtual double me2() const;
/**
* Determine if the configuration is below or above the cutoff.
*/
virtual void checkCutoff();
/**
* Determine all kinematic variables which are not provided by the
* dipole kinematics; store all shower variables in the respective
* dipole object for later use.
*/
virtual void getShowerVariables();
protected:
/**
* Return true, if the shower was able to generate an emission
* leading from the given Born to the given real emission process.
*/
virtual bool isInShowerPhasespace() const;
/**
* Return true, if the shower emission leading from the given Born
* to the given real emission process would have been generated
* above the shower's infrared cutoff.
*/
virtual bool isAboveCutoff() const;
/**
* Calculate qtilde^2 and z for the splitting considered
*/
void calculateShowerVariables() const;
/**
* Return the splitting function as a function of the kinematic
* variables
*/
double splitFn(const pair<Energy2,double>&) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
QTildeMatching & operator=(const QTildeMatching &);
/**
* The shower handler to be used
*/
Ptr<ShowerHandler>::ptr theShowerHandler;
/**
* The qtilde partner finder for calculating the hard scales
*/
- Ptr<QTildeFinder>::ptr theQTildeFinder;
+ Ptr<PartnerFinder>::ptr theQTildeFinder;
/**
* The qtilde Sudakov to access the cutoff
*/
Ptr<QTildeSudakov>::ptr theQTildeSudakov;
/**
* True, if PDF weight should be corrected for z/x mismatch at the
* hard phase space boundary
*/
bool theCorrectForXZMismatch;
};
}
#endif /* Herwig_QTildeMatching_H */
diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc
--- a/Shower/QTilde/Base/PartnerFinder.cc
+++ b/Shower/QTilde/Base/PartnerFinder.cc
@@ -1,536 +1,757 @@
// -*- C++ -*-
//
// PartnerFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PartnerFinder class.
//
#include "PartnerFinder.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
-DescribeAbstractClass<PartnerFinder,Interfaced>
+DescribeClass<PartnerFinder,Interfaced>
describePartnerFinder ("Herwig::PartnerFinder","HwShower.so");
// some useful functions to avoid using #define
namespace {
// return bool if final-state particle
inline bool FS(const tShowerParticlePtr a) {
return a->isFinalState();
}
// return colour line pointer
inline Ptr<ThePEG::ColourLine>::transient_pointer
CL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->colourLines()[index]);
}
// return colour line size
inline unsigned int
CLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->colourLines().size();
}
inline Ptr<ThePEG::ColourLine>::transient_pointer
ACL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->antiColourLines()[index]);
}
inline unsigned int
ACLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->antiColourLines().size();
}
}
void PartnerFinder::persistentOutput(PersistentOStream & os) const {
- os << partnerMethod_ << QEDPartner_ << scaleChoice_;
+ os << partnerMethod_ << QEDPartner_ << scaleChoice_
+ << _finalFinalConditions << _initialFinalDecayConditions
+ << _initialInitialConditions;
}
void PartnerFinder::persistentInput(PersistentIStream & is, int) {
- is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
+ is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_
+ >> _finalFinalConditions >> _initialFinalDecayConditions
+ >>_initialInitialConditions;
}
void PartnerFinder::Init() {
static ClassDocumentation<PartnerFinder> documentation
("This class is responsible for finding the partners for each interaction types ",
"and within the evolution scale range specified by the ShowerVariables ",
"then to determine the initial evolution scales for each pair of partners.");
static Switch<PartnerFinder,int> interfacePartnerMethod
("PartnerMethod",
"Choice of partner finding method for gluon evolution.",
&PartnerFinder::partnerMethod_, 0, false, false);
static SwitchOption interfacePartnerMethodRandom
(interfacePartnerMethod,
"Random",
"Choose partners of a gluon randomly.",
0);
static SwitchOption interfacePartnerMethodMaximum
(interfacePartnerMethod,
"Maximum",
"Choose partner of gluon with largest angle.",
1);
static Switch<PartnerFinder,int> interfaceQEDPartner
("QEDPartner",
"Control of which particles to use as the partner for QED radiation",
&PartnerFinder::QEDPartner_, 0, false, false);
static SwitchOption interfaceQEDPartnerAll
(interfaceQEDPartner,
"All",
"Consider all possible choices which give a positive contribution"
" in the soft limit.",
0);
static SwitchOption interfaceQEDPartnerIIandFF
(interfaceQEDPartner,
"IIandFF",
"Only allow initial-initial or final-final combinations",
1);
static SwitchOption interfaceQEDPartnerIF
(interfaceQEDPartner,
"IF",
"Only allow initial-final combinations",
2);
static Switch<PartnerFinder,int> interfaceScaleChoice
("ScaleChoice",
"The choice of the evolution scales",
&PartnerFinder::scaleChoice_, 0, false, false);
static SwitchOption interfaceScaleChoicePartner
(interfaceScaleChoice,
"Partner",
"Scale of all interactions is that of the evolution partner",
0);
static SwitchOption interfaceScaleChoiceDifferent
(interfaceScaleChoice,
"Different",
"Allow each interaction to have different scales",
1);
+
+ static Switch<PartnerFinder,unsigned int> interfaceFinalFinalConditions
+ ("FinalFinalConditions",
+ "The initial conditions for the shower of a final-final colour connection",
+ &PartnerFinder::_finalFinalConditions, 0, false, false);
+ static SwitchOption interfaceFinalFinalConditionsSymmetric
+ (interfaceFinalFinalConditions,
+ "Symmetric",
+ "The symmetric choice",
+ 0);
+ static SwitchOption interfaceFinalFinalConditionsColoured
+ (interfaceFinalFinalConditions,
+ "Coloured",
+ "Maximal radiation from the coloured particle",
+ 1);
+ static SwitchOption interfaceFinalFinalConditionsAntiColoured
+ (interfaceFinalFinalConditions,
+ "AntiColoured",
+ "Maximal emission from the anticoloured particle",
+ 2);
+ static SwitchOption interfaceFinalFinalConditionsRandom
+ (interfaceFinalFinalConditions,
+ "Random",
+ "Randomly selected maximal emission from one of the particles",
+ 3);
+
+ static Switch<PartnerFinder,unsigned int> interfaceInitialFinalDecayConditions
+ ("InitialFinalDecayConditions",
+ "The initial conditions for the shower of an initial-final"
+ " decay colour connection.",
+ &PartnerFinder::_initialFinalDecayConditions, 0, false, false);
+ static SwitchOption interfaceInitialFinalDecayConditionsSymmetric
+ (interfaceInitialFinalDecayConditions,
+ "Symmetric",
+ "The symmetric choice",
+ 0);
+ static SwitchOption interfaceInitialFinalDecayConditionsMaximal
+ (interfaceInitialFinalDecayConditions,
+ "Maximal",
+ "Maximal radiation from the decay product",
+ 1);
+ static SwitchOption interfaceInitialFinalDecayConditionsSmooth
+ (interfaceInitialFinalDecayConditions,
+ "Smooth",
+ "Smooth matching in the soft limit",
+ 2);
+
+ static Switch<PartnerFinder,unsigned int> interfaceInitialInitialConditions
+ ("InitialInitialConditions",
+ "The initial conditions for the shower of an initial-initial"
+ " colour connection.",
+ &PartnerFinder::_initialInitialConditions, 0, false, false);
+ static SwitchOption interfaceInitialInitialConditionsSymmetric
+ (interfaceInitialInitialConditions,
+ "Symmetric",
+ "The symmetric choice",
+ 0);
+ static SwitchOption interfaceInitialInitialConditionsMaximiseB
+ (interfaceInitialInitialConditions,
+ "MaximiseB",
+ "Maximal radiation from parton b",
+ 1);
+ static SwitchOption interfaceInitialInitialConditionsMaximiseC
+ (interfaceInitialInitialConditions,
+ "MaximiseC",
+ "Maximal radiation from parton c",
+ 2);
}
void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
ShowerInteraction type,
const bool setPartners) {
// clear the existing partners
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) (*cit)->clearPartners();
// set them
if(type==ShowerInteraction::QCD) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
}
else if(type==ShowerInteraction::QED) {
setInitialQEDEvolutionScales(particles,isDecayCase,setPartners);
}
else {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
setInitialQEDEvolutionScales(particles,isDecayCase,false);
}
// print out for debugging
if(Debug::level>=10) {
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
generator()->log() << "Particle: " << **cit << "\n";
if(!(**cit).partner()) continue;
generator()->log() << "Primary partner: " << *(**cit).partner() << "\n";
for(vector<ShowerParticle::EvolutionPartner>::const_iterator it= (**cit).partners().begin();
it!=(**cit).partners().end();++it) {
generator()->log() << static_cast<long>(it->type) << " "
<< it->weight << " "
<< it->scale/GeV << " "
<< *(it->partner)
<< "\n";
}
}
generator()->log() << flush;
}
}
void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// set the partners and the scales
if(setPartners) {
// Loop over particles and consider only coloured particles which don't
// have already their colour partner fixed and that don't have children
// (the latter requirement is relaxed in the case isDecayCase is true).
// Build a map which has as key one of these particles (i.e. a pointer
// to a ShowerParticle object) and as a corresponding value the vector
// of all its possible *normal* candidate colour partners, defined as follows:
// --- have colour, and no children (this is not required in the case
// isDecayCase is true);
// --- if both are initial (incoming) state particles, then the (non-null) colourLine()
// of one of them must match the (non-null) antiColourLine() of the other.
// --- if one is an initial (incoming) state particle and the other is
// a final (outgoing) state particle, then both must have the
// same (non-null) colourLine() or the same (non-null) antiColourLine();
// Notice that this definition exclude the special case of baryon-violating
// processes (as in R-parity Susy), which will show up as particles
// without candidate colour partners, and that we will be treated a part later
// (this means that no modifications of the following loop is needed!)
ShowerParticleVector::const_iterator cit, cjt;
for(cit = particles.begin(); cit != particles.end(); ++cit) {
// Skip colourless particles
if(!(*cit)->data().coloured()) continue;
// find the partners
vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners =
findQCDPartners(*cit,particles);
// must have a partner
if(partners.empty()) {
throw Exception() << "`Failed to make colour connections in "
<< "PartnerFinder::setQCDInitialEvolutionScales"
<< (**cit)
<< Exception::eventerror;
}
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase));
}
// In the case of more than one candidate colour partners,
// there are now two approaches to choosing the partner. The
// first method is based on two assumptions:
// 1) the choice of which is the colour partner is done
// *randomly* between the available candidates.
// 2) the choice of which is the colour partner is done
// *independently* from each particle: in other words,
// if for a particle "i" its selected colour partner is
// the particle "j", then the colour partner of "j"
// does not have to be necessarily "i".
// The second method always chooses the furthest partner
// for hard gluons and gluinos.
// store the choice
int position(-1);
// random choice
if( partnerMethod_ == 0 ) {
// random choice of partner
position = UseRandom::irnd(partners.size());
}
// take the one with largest angle
else if (partnerMethod_ == 1 ) {
if ((*cit)->perturbative() == 1 &&
(*cit)->dataPtr()->iColour()==PDT::Colour8 ) {
assert(partners.size()==2);
// Determine largest angle
double maxAngle(0.);
for(unsigned int ix=0;ix<partners.size();++ix) {
double angle = (*cit)->momentum().vect().
angle(partners[ix].second->momentum().vect());
if(angle>maxAngle) {
maxAngle = angle;
position = ix;
}
}
}
else
position = UseRandom::irnd(partners.size());
}
else
assert(false);
// set the evolution partner
(*cit)->partner(partners[position].second);
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
1.,partners[ix].first,
scales[ix].first));
}
// set scales for all interactions to that of the partner, default
Energy scale = scales[position].first;
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
(**cit).scales().QCD_c =
(**cit).scales().QCD_c_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
(**cit).scales().QCD_ac =
(**cit).scales().QCD_ac_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else
assert(false);
}
}
}
// primary partner set, set the others and do the scale
else {
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
// Skip colourless particles
if(!(*cit)->data().coloured()) continue;
// find the partners
vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners =
findQCDPartners(*cit,particles);
// must have a partner
if(partners.empty()) {
throw Exception() << "`Failed to make colour connections in "
<< "PartnerFinder::setQCDInitialEvolutionScales"
<< (**cit)
<< Exception::eventerror;
}
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
int position(-1);
for(unsigned int ix=0;ix< partners.size();++ix) {
if(partners[ix].second) position = ix;
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase));
}
assert(position>=0);
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
1.,partners[ix].first,
scales[ix].first));
}
// set scales for all interactions to that of the partner, default
Energy scale = scales[position].first;
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
(**cit).scales().QCD_c =
(**cit).scales().QCD_c_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
(**cit).scales().QCD_ac =
(**cit).scales().QCD_ac_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else {
assert(false);
}
}
}
}
}
void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// loop over all the particles
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
// not charged or photon continue
if(!(**cit).dataPtr()->charged()) continue;
// find the potential partners
vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(*cit,particles,isDecayCase);
if(partners.empty()) {
throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
}
// calculate the probabilities
double prob(0.);
for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
// normalise
for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
// set the partner if required
int position(-1);
// use QCD partner if set
if(!setPartners&&(*cit)->partner()) {
for(unsigned int ix=0;ix<partners.size();++ix) {
if((*cit)->partner()==partners[ix].second) {
position = ix;
break;
}
}
}
// set the partner
if(setPartners||!(*cit)->partner()||position<0) {
prob = UseRandom::rnd();
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first>prob) {
position = ix;
break;
}
prob -= partners[ix].first;
}
if(position>=0&&(setPartners||!(*cit)->partner())) {
(*cit)->partner(partners[position].second);
}
}
// must have a partner
if(position<0) throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase));
}
// store all the possible partners
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
partners[ix].first,
ShowerPartnerType::QED,
scales[ix].first));
}
// set scales
(**cit).scales().QED = scales[position].first;
(**cit).scales().QED_noAO = scales[position].first;
}
}
pair<Energy,Energy> PartnerFinder::
calculateInitialEvolutionScales(const ShowerPPair &particlePair,
const bool isDecayCase) {
bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
- if(FS1 && FS2)
- return calculateFinalFinalScales(particlePair);
+ if(FS1 && FS2) {
+ bool colouredFirst =
+ particlePair.first->colourLine()&&
+ particlePair.first->colourLine()==particlePair.second->antiColourLine();
+ return calculateFinalFinalScales(particlePair.first->momentum(),
+ particlePair.second->momentum(),
+ colouredFirst);
+ }
else if(FS1 && !FS2) {
ShowerPPair a(particlePair.second, particlePair.first);
- pair<Energy,Energy> rval = calculateInitialFinalScales(a,isDecayCase);
+ pair<Energy,Energy> rval = calculateInitialFinalScales(a.first->momentum(),
+ a.second->momentum(),
+ isDecayCase);
return pair<Energy,Energy>(rval.second,rval.first);
}
else if(!FS1 &&FS2)
- return calculateInitialFinalScales(particlePair,isDecayCase);
+ return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
else
- return calculateInitialInitialScales(particlePair);
+ return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
}
vector< pair<ShowerPartnerType, tShowerParticlePtr> >
PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles) {
vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
ShowerParticleVector::const_iterator cjt;
for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
if(!(*cjt)->data().coloured() || particle==*cjt) continue;
// one initial-state and one final-state particle
if(FS(particle) != FS(*cjt)) {
// loop over all the colours of both particles
for(unsigned int ix=0; ix<CLSIZE(particle); ++ix) {
for(unsigned int jx=0; jx<CLSIZE(*cjt); ++jx) {
if((CL(particle,ix) && CL(particle,ix)==CL(*cjt,jx))) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
//loop over all the anti-colours of both particles
for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix) {
for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx) {
if((ACL(particle,ix) && ACL(particle,ix)==ACL(*cjt,jx))) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
}
// two initial-state or two final-state particles
else {
//loop over the colours of the first particle and the anti-colours of the other
for(unsigned int ix=0; ix<CLSIZE(particle); ++ix){
for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx){
if(CL(particle,ix) && CL(particle,ix)==ACL(*cjt,jx)) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
//loop over the anti-colours of the first particle and the colours of the other
for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix){
for(unsigned int jx=0; jx<CLSIZE(*cjt); jx++){
if(ACL(particle,ix) && ACL(particle,ix)==CL(*cjt,jx)) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
}
}
// if we haven't found any partners look for RPV
if (partners.empty()) {
// special for RPV
tColinePtr col = CL(particle);
if(FS(particle)&&col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))||
(!FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second ))) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
else if(col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))||
(!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
col = ACL(particle);
if(FS(particle)&&col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))||
(!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second ))) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
else if(col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))||
(!FS(*cjt) && (ACL(*cjt) == cpair.first ||ACL(*cjt) == cpair.second))) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
}
// return the partners
return partners;
}
vector< pair<double, tShowerParticlePtr> >
PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles,
const bool isDecayCase) {
vector< pair<double, tShowerParticlePtr> > partners;
ShowerParticleVector::const_iterator cjt;
double pcharge = particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
vector< pair<double, tShowerParticlePtr> > photons;
for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
if(particle == *cjt) continue;
if((**cjt).id()==ParticleID::gamma) photons.push_back(make_pair(1.,*cjt));
if(!(*cjt)->data().charged() ) continue;
double charge = pcharge*double((*cjt)->data().iCharge());
if( FS(particle) != FS(*cjt) ) charge *=-1.;
if( QEDPartner_ != 0 && !isDecayCase ) {
// only include II and FF as requested
if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) )
continue;
// only include IF is requested
else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) )
continue;
}
if(particle->id()==ParticleID::gamma) charge = -abs(charge);
// only keep positive dipoles
if(charge<0.) partners.push_back(make_pair(-charge,*cjt));
}
if(particle->id()==ParticleID::gamma&& partners.empty()) {
return photons;
}
return partners;
}
+
+pair<Energy,Energy>
+PartnerFinder::calculateFinalFinalScales(
+ const Lorentz5Momentum & p1,
+ const Lorentz5Momentum & p2,
+ bool colouredFirst)
+{
+ static const double eps=1e-7;
+ // Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
+ // ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
+ // find momenta in rest frame of system
+ // calculate quantities for the scales
+ Energy2 Q2 = (p1+p2).m2();
+ double b = p1.mass2()/Q2;
+ double c = p2.mass2()/Q2;
+ if(b<0.) {
+ if(b<-eps) {
+ throw Exception() << "Negative Mass squared b = " << b
+ << "in PartnerFinder::calculateFinalFinalScales()"
+ << Exception::eventerror;
+ }
+ b = 0.;
+ }
+ if(c<0.) {
+ if(c<-eps) {
+ throw Exception() << "Negative Mass squared c = " << c
+ << "in PartnerFinder::calculateFinalFinalScales()"
+ << Exception::eventerror;
+ }
+ c = 0.;
+ }
+ // KMH & PR - 16 May 2008 - swapped lambda calculation from
+ // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
+ // which should be identical for p1 & p2 onshell in their COM
+ // but in the inverse construction for the Nason method, this
+ // was not the case, leading to misuse.
+ double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
+ *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b)));
+ // symmetric case
+ unsigned int iopt=finalFinalConditions();
+ Energy firstQ,secondQ;
+ if(iopt==0) {
+ firstQ = sqrt(0.5*Q2*(1.+b-c+lam));
+ secondQ = sqrt(0.5*Q2*(1.-b+c+lam));
+ }
+ // assymetric choice
+ else {
+ double kappab,kappac;
+ // calculate kappa with coloured line getting maximum
+ if((iopt==1&&colouredFirst)|| // first particle coloured+maximal for coloured
+ (iopt==2&&!colouredFirst)|| // first particle anticoloured+maximal for acoloured
+ (iopt==3&&UseRandom::rndbool(0.5))) { // random choice
+ kappab=4.*(1.-2.*sqrt(c)-b+c);
+ kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b);
+ }
+ else {
+ kappac=4.*(1.-2.*sqrt(b)-c+b);
+ kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c);
+ }
+ // calculate the scales
+ firstQ = sqrt(Q2*kappab);
+ secondQ = sqrt(Q2*kappac);
+ }
+ return pair<Energy,Energy>(firstQ, secondQ);
+}
+
+
+pair<Energy,Energy>
+PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
+ const bool isDecayCase) {
+ if(!isDecayCase) {
+ // In this case from JHEP 12(2003)045 we find the conditions
+ // ktilde_b = (1+c) and ktilde_c = (1+2c)
+ // We also find that c = m_c^2/Q^2. The process is a+b->c where
+ // particle a is not colour connected (considered as a colour singlet).
+ // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
+ // q_c = sqrt(Q^2+2 m_c^2)
+ // We also assume that the first particle in the pair is the initial
+ // state particle and the second is the final state one c
+ Energy2 mc2 = sqr(pc.mass());
+ Energy2 Q2 = -(pb-pc).m2();
+ return pair<Energy,Energy>(sqrt(Q2+mc2), sqrt(Q2+2*mc2));
+ }
+ else {
+ // In this case from JHEP 12(2003)045 we find, for the decay
+ // process b->c+a(neutral), the condition
+ // (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda).
+ // We also assume that the first particle in the pair is the initial
+ // state particle (b) and the second is the final state one (c).
+ // - We find maximal phase space coverage through emissions from
+ // c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
+ // - We find the most 'symmetric' way to populate the phase space
+ // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
+ // - We find the most 'smooth' way to populate the phase space
+ // occurs for...
+ Energy2 mb2(sqr(pb.mass()));
+ double a=(pb-pc).m2()/mb2;
+ double c=sqr(pc.mass())/mb2;
+ double lambda = 1. + a*a + c*c - 2.*a - 2.*c - 2.*a*c;
+ lambda = sqrt(max(lambda,0.));
+ double PROD = 0.25*sqr(1. - a + c + lambda);
+ double ktilde_b, ktilde_c,cosi(0.);
+ switch(initialFinalDecayConditions()) {
+ case 0: // the 'symmetric' choice
+ ktilde_c = 0.5*(1-a+c+lambda) + c ;
+ ktilde_b = 1.+PROD/(ktilde_c-c) ;
+ break;
+ case 1: // the 'maximal' choice
+ ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
+ ktilde_b = 1.+PROD/(ktilde_c-c) ;
+ break;
+ case 2: // the 'smooth' choice
+ // c is a problem if very small here use 1GeV as minimum
+ c = max(c,1.*GeV2/mb2);
+ cosi = (sqr(1-sqrt(c))-a)/lambda;
+ ktilde_b = 2.0/(1.0-cosi);
+ ktilde_c = (1.0-a+c+lambda)*(1.0+c-a-lambda*cosi)/(2.0*(1.0+cosi));
+ break;
+ default:
+ throw Exception() << "Invalid option for decay shower's phase space"
+ << " PartnerFinder::calculateInitialFinalScales"
+ << Exception::abortnow;
+ }
+ return pair<Energy,Energy>(sqrt(mb2*ktilde_b),sqrt(mb2*ktilde_c));
+ }
+}
+
+pair<Energy,Energy>
+PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
+ // This case is quite simple. From JHEP 12(2003)045 we find the condition
+ // that ktilde_b = ktilde_c = 1. In this case we have the process
+ // b+c->a so we need merely boost to the CM frame of the two incoming
+ // particles and then qtilde is equal to the energy in that frame
+ Energy Q = sqrt((p1+p2).m2());
+ if(_initialInitialConditions==1) {
+ return pair<Energy,Energy>(sqrt(2.0)*Q,sqrt(0.5)*Q);
+ } else if(_initialInitialConditions==2) {
+ return pair<Energy,Energy>(sqrt(0.5)*Q,sqrt(2.0)*Q);
+ } else {
+ return pair<Energy,Energy>(Q,Q);
+ }
+}
diff --git a/Shower/QTilde/Base/PartnerFinder.h b/Shower/QTilde/Base/PartnerFinder.h
--- a/Shower/QTilde/Base/PartnerFinder.h
+++ b/Shower/QTilde/Base/PartnerFinder.h
@@ -1,208 +1,297 @@
// -*- C++ -*-
//
// PartnerFinder.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_PartnerFinder_H
#define HERWIG_PartnerFinder_H
//
// This is the declaration of the PartnerFinder class.
//
#include "Herwig/Shower/QTilde/ShowerConfig.h"
#include "ThePEG/Interface/Interfaced.h"
#include "PartnerFinder.fh"
namespace Herwig {
using namespace ThePEG;
/**
* typedef of a pair of particle for calculating the evolution scales
*/
typedef pair<tShowerParticlePtr,tShowerParticlePtr> ShowerPPair;
/** \ingroup Shower
*
* This class is responsible of two related tasks:
* - it finds the partners
* - for each pair of partners (and interaction therefore)
* it sets the initial evolution scales of both of them.
*
* In general the finding of the partners is performed by this class but
* the calculation of the initial evolution scales should be implemented
* for different shower evolution models in classes inheriting from this one.
* Notice that a given particle has, in general, a different partner
* for each different interaction; however, given a partner, its
* initial evolution scale, Q, is purely a kinematical relationship
* between the pair, without dependence on the dynamics (i.e. type of interaction).
*
* @see \ref PartnerFinderInterfaces "The interfaces"
* defined for PartnerFinder.
*/
class PartnerFinder: public Interfaced {
public:
/**
* The default constructor.
*/
- PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0) {}
+ PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0),
+ _finalFinalConditions(0),
+ _initialFinalDecayConditions(0),
+ _initialInitialConditions(0) {}
/**
* Given in input a collection of particles (ShowerParticle objects),
* each of these methods set the initial evolution scales of those particles,
* between the ones given in input, that do not have yet their
* evolution scale set.
* The input collection of particles can be either the full collection of
* showering particles (kept in the main class ShowerHandler,
* in the case isDecayCase is false, or simply, in the case isDecayCase
* is true, the decaying particle and its decay products.
* The methods returns true, unless something wrong (inconsistencies,
* or undefined values) happens.
*
* These methods are virtual but in most cases inheriting classes should not
* need to overide them as they simply find the relevant partner and call
* one of the calculateScale members to calculate the scale.
*/
//@{
/**
* Set the initial scales
* @param particles The particles to be considered
* @param isDecayCase Whether or not this is a decay
* @param setPartners Whether to set the colour partners or just the scales
*/
virtual void setInitialEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
ShowerInteraction,
const bool setPartners=true);
//@}
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
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:
/**
+ * Access function for the initial conditions for the shower
+ */
+ //@{
+ /**
+ * Initial conditions for the shower of a final-final colour connection
+ * - 0 is the symmetric choice
+ * - 1 is maximal emmision from the coloured particle
+ * - 2 is maximal emmision from the anticoloured particle
+ * - 3 is randomly selected maximal emmision
+ */
+ unsigned int finalFinalConditions() const
+ {return _finalFinalConditions;}
+
+ /**
+ * Initial conditions for the shower of an initial-final decay colour connection
+ * - 0 is the symmetric choice
+ * - 1 is maximal emission from the decay product
+ * - 2 is the smooth choice
+ */
+ unsigned int initialFinalDecayConditions() const
+ {return _initialFinalDecayConditions;}
+ //@}
+
+protected:
+
+ /**
* Members to set the scales for different interactions
*/
//@{
/**
* Set initial scales for a QCD interaction
*/
virtual void setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners=true);
/**
* Set initial scales for a QED interaction
*/
virtual void setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners=true);
//@}
/**
* Find the QCD partners
* @param particle The particle to find the partners for
* @param particles The full set of particles to search
*/
vector< pair<ShowerPartnerType, tShowerParticlePtr> >
findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles);
/**
* Find the QED partners
* @param particle The particle to find the partners for
* @param particles The full set of particles to search
*/
vector< pair<double, tShowerParticlePtr> >
findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles,
const bool isDecayCase);
+public:
/**
* Given a pair of particles, supposedly partners w.r.t. an interaction,
* this method returns their initial evolution scales as a pair.
* If something wrong happens, it returns the null (ZERO,ZERO) pair.
* This method is used by the above setXXXInitialEvolutionScales
* methods.
* These methods must be overiden in inheriting classes
*/
//@{
/**
* General method to calculate the initial evolution scales
*/
- virtual pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
- const bool isDecayCase);
+ pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
+ const bool isDecayCase);
/**
- * Calculate the initial evolution scales for two final-state particles
+ * Calculate the initial evolution scales given momenta
*/
- virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &)=0;
+ pair<Energy,Energy> calculateFinalFinalScales(const Lorentz5Momentum & p1,
+ const Lorentz5Momentum & p2,
+ bool colouredfirst);
/**
- * Calculate the initial evolution scales for two initial-state particles
+ * Calculate the initial evolution scales given momenta
*/
- virtual pair<Energy,Energy> calculateInitialInitialScales(const ShowerPPair &)=0;
+ pair<Energy,Energy> calculateInitialInitialScales(const Lorentz5Momentum& p1,
+ const Lorentz5Momentum& p2);
/**
- * Calculate the initial evolution scales for one initial
- * and one final-state particles
+ * Calculate the initial evolution scales given momenta
*/
- virtual pair<Energy,Energy> calculateInitialFinalScales(const ShowerPPair &,
- const bool isDecayCase)=0;
+ pair<Energy,Energy> calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
+ const bool isDecayCase);
+
+
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
PartnerFinder & operator=(const PartnerFinder &);
private:
/**
* Method for choosing colour partner
*/
int partnerMethod_;
/**
* Choice for the QED radiation partner
*/
int QEDPartner_;
/**
* Choice of the scale
*/
int scaleChoice_;
+
+ /**
+ * Flags controlling the initial conditions for the shower
+ */
+ //@{
+ /**
+ * Initial conditions for the shower with a final-final colour
+ * connection
+ */
+ unsigned int _finalFinalConditions;
+
+ /**
+ * Initial conditions for the shower with an initial-final decay colour
+ * connection. This is done according to the top decay colour
+ * connection calculation in JHEP12(2003)_045. The options act as follows:
+ * 0: This is the default 'symmetric' choice which more or less divides
+ * the phase space evenly between the parent and its charged child.
+ * 1: This 'maximal' choice maximises the phase space available for
+ * gluons emitted from the charged child.
+ * 2: This (experimental) 'smooth' choice does not suffer from
+ * a discontinuity at the boundary between the region populated by
+ * emissions from the charged child and the region populated by emissions
+ * from the parent. This does, however, mean that the phase space
+ * available for emissions from the charged child is fairly minimal.
+ */
+ unsigned int _initialFinalDecayConditions;
+
+ /**
+ * Initial conditions for the shower with an initial-initial colour
+ * connection. This is done according to the colour connection
+ * calculation in JHEP12(2003)_045. The options act as follows:
+ * 0: This is the default 'symmetric' choice which more or less divides
+ * the phase space evenly between the two incoming partons.
+ * 1: This increases the phase space for emission from "parton b".
+ * 2: This increases the phase space for emission from "parton c".
+ */
+ unsigned int _initialInitialConditions;
+ //@}
};
}
#endif /* HERWIG_PartnerFinder_H */
diff --git a/Shower/QTilde/Default/QTildeFinder.cc b/Shower/QTilde/Default/QTildeFinder.cc
deleted file mode 100644
--- a/Shower/QTilde/Default/QTildeFinder.cc
+++ /dev/null
@@ -1,273 +0,0 @@
-// -*- C++ -*-
-//
-// QTildeFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
-// Copyright (C) 2002-2017 The Herwig Collaboration
-//
-// Herwig is licenced under version 3 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
-//
-// This is the implementation of the non-inlined, non-templated member
-// functions of the QTildeFinder class.
-//
-
-#include "QTildeFinder.h"
-#include "ThePEG/Interface/ClassDocumentation.h"
-#include "ThePEG/Persistency/PersistentOStream.h"
-#include "ThePEG/Persistency/PersistentIStream.h"
-#include "ThePEG/Interface/Switch.h"
-#include "ThePEG/Repository/EventGenerator.h"
-#include "ThePEG/EventRecord/Event.h"
-#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
-#include "ThePEG/Repository/UseRandom.h"
-#include "ThePEG/Utilities/DescribeClass.h"
-
-using namespace Herwig;
-
-DescribeClass<QTildeFinder,Herwig::PartnerFinder>
-describeQTildeFinder ("Herwig::QTildeFinder","HwShower.so");
-
-void QTildeFinder::persistentOutput(PersistentOStream & os) const {
- os << _finalFinalConditions << _initialFinalDecayConditions
- << _initialInitialConditions;
-}
-
-void QTildeFinder::persistentInput(PersistentIStream & is, int) {
- is >> _finalFinalConditions >> _initialFinalDecayConditions
- >>_initialInitialConditions;
-}
-
-void QTildeFinder::Init() {
-
- static ClassDocumentation<QTildeFinder> documentation
- ("This class is responsible for finding the partners for each interaction types ",
- "and within the evolution scale range specified by the ShowerVariables ",
- "then to determine the initial evolution scales for each pair of partners.");
-
- static Switch<QTildeFinder,unsigned int> interfaceFinalFinalConditions
- ("FinalFinalConditions",
- "The initial conditions for the shower of a final-final colour connection",
- &QTildeFinder::_finalFinalConditions, 0, false, false);
- static SwitchOption interfaceFinalFinalConditionsSymmetric
- (interfaceFinalFinalConditions,
- "Symmetric",
- "The symmetric choice",
- 0);
- static SwitchOption interfaceFinalFinalConditionsColoured
- (interfaceFinalFinalConditions,
- "Coloured",
- "Maximal radiation from the coloured particle",
- 1);
- static SwitchOption interfaceFinalFinalConditionsAntiColoured
- (interfaceFinalFinalConditions,
- "AntiColoured",
- "Maximal emission from the anticoloured particle",
- 2);
- static SwitchOption interfaceFinalFinalConditionsRandom
- (interfaceFinalFinalConditions,
- "Random",
- "Randomly selected maximal emission from one of the particles",
- 3);
-
- static Switch<QTildeFinder,unsigned int> interfaceInitialFinalDecayConditions
- ("InitialFinalDecayConditions",
- "The initial conditions for the shower of an initial-final"
- " decay colour connection.",
- &QTildeFinder::_initialFinalDecayConditions, 0, false, false);
- static SwitchOption interfaceInitialFinalDecayConditionsSymmetric
- (interfaceInitialFinalDecayConditions,
- "Symmetric",
- "The symmetric choice",
- 0);
- static SwitchOption interfaceInitialFinalDecayConditionsMaximal
- (interfaceInitialFinalDecayConditions,
- "Maximal",
- "Maximal radiation from the decay product",
- 1);
- static SwitchOption interfaceInitialFinalDecayConditionsSmooth
- (interfaceInitialFinalDecayConditions,
- "Smooth",
- "Smooth matching in the soft limit",
- 2);
-
- static Switch<QTildeFinder,unsigned int> interfaceInitialInitialConditions
- ("InitialInitialConditions",
- "The initial conditions for the shower of an initial-initial"
- " colour connection.",
- &QTildeFinder::_initialInitialConditions, 0, false, false);
- static SwitchOption interfaceInitialInitialConditionsSymmetric
- (interfaceInitialInitialConditions,
- "Symmetric",
- "The symmetric choice",
- 0);
- static SwitchOption interfaceInitialInitialConditionsMaximiseB
- (interfaceInitialInitialConditions,
- "MaximiseB",
- "Maximal radiation from parton b",
- 1);
- static SwitchOption interfaceInitialInitialConditionsMaximiseC
- (interfaceInitialInitialConditions,
- "MaximiseC",
- "Maximal radiation from parton c",
- 2);
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateInitialFinalScales(const ShowerPPair &ppair, const bool isDecayCase) {
- return
- calculateInitialFinalScales(ppair.first->momentum(),ppair.second->momentum(),isDecayCase);
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
- const bool isDecayCase) {
- if(!isDecayCase) {
- // In this case from JHEP 12(2003)045 we find the conditions
- // ktilde_b = (1+c) and ktilde_c = (1+2c)
- // We also find that c = m_c^2/Q^2. The process is a+b->c where
- // particle a is not colour connected (considered as a colour singlet).
- // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
- // q_c = sqrt(Q^2+2 m_c^2)
- // We also assume that the first particle in the pair is the initial
- // state particle and the second is the final state one c
- Energy2 mc2 = sqr(pc.mass());
- Energy2 Q2 = -(pb-pc).m2();
- return pair<Energy,Energy>(sqrt(Q2+mc2), sqrt(Q2+2*mc2));
- }
- else {
- // In this case from JHEP 12(2003)045 we find, for the decay
- // process b->c+a(neutral), the condition
- // (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda).
- // We also assume that the first particle in the pair is the initial
- // state particle (b) and the second is the final state one (c).
- // - We find maximal phase space coverage through emissions from
- // c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
- // - We find the most 'symmetric' way to populate the phase space
- // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
- // - We find the most 'smooth' way to populate the phase space
- // occurs for...
- Energy2 mb2(sqr(pb.mass()));
- double a=(pb-pc).m2()/mb2;
- double c=sqr(pc.mass())/mb2;
- double lambda = 1. + a*a + c*c - 2.*a - 2.*c - 2.*a*c;
- lambda = sqrt(max(lambda,0.));
- double PROD = 0.25*sqr(1. - a + c + lambda);
- double ktilde_b, ktilde_c,cosi(0.);
- switch(initialFinalDecayConditions()) {
- case 0: // the 'symmetric' choice
- ktilde_c = 0.5*(1-a+c+lambda) + c ;
- ktilde_b = 1.+PROD/(ktilde_c-c) ;
- break;
- case 1: // the 'maximal' choice
- ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
- ktilde_b = 1.+PROD/(ktilde_c-c) ;
- break;
- case 2: // the 'smooth' choice
- // c is a problem if very small here use 1GeV as minimum
- c = max(c,1.*GeV2/mb2);
- cosi = (sqr(1-sqrt(c))-a)/lambda;
- ktilde_b = 2.0/(1.0-cosi);
- ktilde_c = (1.0-a+c+lambda)*(1.0+c-a-lambda*cosi)/(2.0*(1.0+cosi));
- break;
- default:
- throw Exception() << "Invalid option for decay shower's phase space"
- << " QTildeFinder::calculateInitialFinalScales"
- << Exception::abortnow;
- }
- return pair<Energy,Energy>(sqrt(mb2*ktilde_b),sqrt(mb2*ktilde_c));
- }
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateInitialInitialScales(const ShowerPPair &ppair) {
- return
- calculateInitialInitialScales(ppair.first->momentum(),
- ppair.second->momentum());
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
- // This case is quite simple. From JHEP 12(2003)045 we find the condition
- // that ktilde_b = ktilde_c = 1. In this case we have the process
- // b+c->a so we need merely boost to the CM frame of the two incoming
- // particles and then qtilde is equal to the energy in that frame
- Energy Q = sqrt((p1+p2).m2());
- if(_initialInitialConditions==1) {
- return pair<Energy,Energy>(sqrt(2.0)*Q,sqrt(0.5)*Q);
- } else if(_initialInitialConditions==2) {
- return pair<Energy,Energy>(sqrt(0.5)*Q,sqrt(2.0)*Q);
- } else {
- return pair<Energy,Energy>(Q,Q);
- }
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateFinalFinalScales(const ShowerPPair & pp) {
- bool colouredFirst =
- pp.first->colourLine()&&
- pp.first->colourLine()==pp.second->antiColourLine();
- return calculateFinalFinalScales(pp.first->momentum(),pp.second->momentum(),
- colouredFirst);
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateFinalFinalScales(Lorentz5Momentum p1, Lorentz5Momentum p2,
- bool colouredFirst) {
- static const double eps=1e-7;
- // Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
- // ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
- // find momenta in rest frame of system
- // calculate quantities for the scales
- Energy2 Q2 = (p1+p2).m2();
- double b = p1.mass2()/Q2;
- double c = p2.mass2()/Q2;
- if(b<0.) {
- if(b<-eps) {
- throw Exception() << "Negative Mass squared b = " << b
- << "in QTildeFinder::calculateFinalFinalScales()"
- << Exception::eventerror;
- }
- b = 0.;
- }
- if(c<0.) {
- if(c<-eps) {
- throw Exception() << "Negative Mass squared c = " << c
- << "in QTildeFinder::calculateFinalFinalScales()"
- << Exception::eventerror;
- }
- c = 0.;
- }
- // KMH & PR - 16 May 2008 - swapped lambda calculation from
- // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
- // which should be identical for p1 & p2 onshell in their COM
- // but in the inverse construction for the Nason method, this
- // was not the case, leading to misuse.
- double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
- *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b)));
- // symmetric case
- unsigned int iopt=finalFinalConditions();
- Energy firstQ,secondQ;
- if(iopt==0) {
- firstQ = sqrt(0.5*Q2*(1.+b-c+lam));
- secondQ = sqrt(0.5*Q2*(1.-b+c+lam));
- }
- // assymetric choice
- else {
- double kappab,kappac;
- // calculate kappa with coloured line getting maximum
- if((iopt==1&&colouredFirst)|| // first particle coloured+maximal for coloured
- (iopt==2&&!colouredFirst)|| // first particle anticoloured+maximal for acoloured
- (iopt==3&&UseRandom::rndbool(0.5))) { // random choice
- kappab=4.*(1.-2.*sqrt(c)-b+c);
- kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b);
- }
- else {
- kappac=4.*(1.-2.*sqrt(b)-c+b);
- kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c);
- }
- // calculate the scales
- firstQ = sqrt(Q2*kappab);
- secondQ = sqrt(Q2*kappac);
- }
- return pair<Energy,Energy>(firstQ, secondQ);
-}
diff --git a/Shower/QTilde/Default/QTildeFinder.h b/Shower/QTilde/Default/QTildeFinder.h
deleted file mode 100644
--- a/Shower/QTilde/Default/QTildeFinder.h
+++ /dev/null
@@ -1,206 +0,0 @@
-// -*- C++ -*-
-//
-// QTildeFinder.h is a part of Herwig - A multi-purpose Monte Carlo event generator
-// Copyright (C) 2002-2017 The Herwig Collaboration
-//
-// Herwig is licenced under version 3 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
-#ifndef HERWIG_QTildeFinder_H
-#define HERWIG_QTildeFinder_H
-//
-// This is the declaration of the QTildeFinder class.
-//
-
-#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
-#include "Herwig/Shower/QTilde/ShowerConfig.h"
-#include "ThePEG/Interface/Interfaced.h"
-
-namespace Herwig {
-using namespace ThePEG;
-
-/** \ingroup Shower
- *
- * The QTildeFinder class is responsible for finding the partners and
- * setting the initial evolution scales for the shower evolution described
- * in JHEP 0312:045,2003.
- *
- * @see \ref QTildeFinderInterfaces "The interfaces"
- * defined for QTildeFinder.
- */
-class QTildeFinder: public PartnerFinder {
-
-public:
-
- /**
- * The default constructor.
- */
- QTildeFinder() : _finalFinalConditions(0),
- _initialFinalDecayConditions(0),
- _initialInitialConditions(0) {}
-
- /** @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();
-
-public:
-
- /**
- * Calculate the initial evolution scales given momenta
- */
- pair<Energy,Energy> calculateFinalFinalScales(Lorentz5Momentum p1, Lorentz5Momentum p2,
- bool colouredfirst);
-
- /**
- * Calculate the initial evolution scales given momenta
- */
- pair<Energy,Energy> calculateInitialInitialScales(const Lorentz5Momentum& p1,
- const Lorentz5Momentum& p2);
-
- /**
- * Calculate the initial evolution scales given momenta
- */
- pair<Energy,Energy> calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
- const bool isDecayCase);
-
-protected:
-
- /**
- * Given a pair of particles, supposedly partners w.r.t. an interaction,
- * this method returns their initial evolution scales as a pair.
- * If something wrong happens, it returns the null (ZERO,ZERO) pair.
- * This method is used by the above setXXXInitialEvolutionScales
- * methods.
- */
- //@{
- /**
- * Calculate the initial evolution scales for two final-state particles
- */
- virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &);
-
- /**
- * Calculate the initial evolution scales for two initial-state particles
- */
- virtual pair<Energy,Energy> calculateInitialInitialScales(const ShowerPPair &);
-
- /**
- * Calculate the initial evolution scales for one initial
- * and one final-state particles
- */
- virtual pair<Energy,Energy> calculateInitialFinalScales(const ShowerPPair &,
- const bool isDecayCase);
- //@}
-
- /**
- * Access function for the initial conditions for the shower
- */
- //@{
- /**
- * Initial conditions for the shower of a final-final colour connection
- * - 0 is the symmetric choice
- * - 1 is maximal emmision from the coloured particle
- * - 2 is maximal emmision from the anticoloured particle
- * - 3 is randomly selected maximal emmision
- */
- unsigned int finalFinalConditions() const
- {return _finalFinalConditions;}
-
- /**
- * Initial conditions for the shower of an initial-final decay colour connection
- * - 0 is the symmetric choice
- * - 1 is maximal emission from the decay product
- * - 2 is the smooth choice
- */
- unsigned int initialFinalDecayConditions() const
- {return _initialFinalDecayConditions;}
- //@}
-
-protected:
-
- /** @name Clone Methods. */
- //@{
- /**
- * Make a simple clone of this object.
- * @return a pointer to the new object.
- */
- virtual IBPtr clone() const {return new_ptr(*this);}
-
- /** 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 {return new_ptr(*this);}
- //@}
-
-private:
-
- /**
- * The assignment operator is private and must never be called.
- * In fact, it should not even be implemented.
- */
- QTildeFinder & operator=(const QTildeFinder &);
-
-private:
-
- /**
- * Flags controlling the initial conditions for the shower
- */
- //@{
- /**
- * Initial conditions for the shower with a final-final colour
- * connection
- */
- unsigned int _finalFinalConditions;
-
- /**
- * Initial conditions for the shower with an initial-final decay colour
- * connection. This is done according to the top decay colour
- * connection calculation in JHEP12(2003)_045. The options act as follows:
- * 0: This is the default 'symmetric' choice which more or less divides
- * the phase space evenly between the parent and its charged child.
- * 1: This 'maximal' choice maximises the phase space available for
- * gluons emitted from the charged child.
- * 2: This (experimental) 'smooth' choice does not suffer from
- * a discontinuity at the boundary between the region populated by
- * emissions from the charged child and the region populated by emissions
- * from the parent. This does, however, mean that the phase space
- * available for emissions from the charged child is fairly minimal.
- */
- unsigned int _initialFinalDecayConditions;
-
- /**
- * Initial conditions for the shower with an initial-initial colour
- * connection. This is done according to the colour connection
- * calculation in JHEP12(2003)_045. The options act as follows:
- * 0: This is the default 'symmetric' choice which more or less divides
- * the phase space evenly between the two incoming partons.
- * 1: This increases the phase space for emission from "parton b".
- * 2: This increases the phase space for emission from "parton c".
- */
- unsigned int _initialInitialConditions;
- //@}
-};
-
-}
-
-#endif /* HERWIG_QTildeFinder_H */
diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am
--- a/Shower/QTilde/Makefile.am
+++ b/Shower/QTilde/Makefile.am
@@ -1,43 +1,42 @@
SUBDIRS = Matching
pkglib_LTLIBRARIES = HwShower.la
HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 25:0:0
HwShower_la_SOURCES = \
Couplings/ShowerAlphaQCD.h Couplings/ShowerAlphaQCD.cc \
Couplings/ShowerAlphaQED.h Couplings/ShowerAlphaQED.cc\
QTildeShowerHandler.h QTildeShowerHandler.fh QTildeShowerHandler.cc \
SplittingFunctions/HalfHalfOneSplitFn.h SplittingFunctions/HalfHalfOneSplitFn.cc\
SplittingFunctions/OneOneOneSplitFn.h SplittingFunctions/OneOneOneSplitFn.cc\
SplittingFunctions/OneOneOneMassiveSplitFn.h SplittingFunctions/OneOneOneMassiveSplitFn.cc\
SplittingFunctions/ZeroZeroOneSplitFn.h SplittingFunctions/ZeroZeroOneSplitFn.cc\
SplittingFunctions/OneHalfHalfSplitFn.h SplittingFunctions/OneHalfHalfSplitFn.cc\
SplittingFunctions/HalfOneHalfSplitFn.h SplittingFunctions/HalfOneHalfSplitFn.cc\
SplittingFunctions/CMWOneOneOneSplitFn.h SplittingFunctions/CMWOneOneOneSplitFn.cc\
SplittingFunctions/CMWHalfHalfOneSplitFn.h SplittingFunctions/CMWHalfHalfOneSplitFn.cc\
Default/QTildeSudakov.cc Default/QTildeSudakov.h\
Default/Decay_QTildeShowerKinematics1to2.cc \
Default/Decay_QTildeShowerKinematics1to2.h \
Default/IS_QTildeShowerKinematics1to2.cc Default/IS_QTildeShowerKinematics1to2.h \
Default/FS_QTildeShowerKinematics1to2.cc Default/FS_QTildeShowerKinematics1to2.h \
-Default/QTildeFinder.cc Default/QTildeFinder.h\
Default/QTildeReconstructor.cc Default/QTildeReconstructor.h Default/QTildeReconstructor.tcc \
Base/KinematicsReconstructor.cc \
Base/KinematicsReconstructor.h \
Base/KinematicsReconstructor.fh \
Base/HardTree.cc Base/HardTree.h Base/HardTree.fh \
Base/HardBranching.h Base/HardBranching.fh Base/HardBranching.cc\
Base/PartnerFinder.h Base/PartnerFinder.fh Base/PartnerFinder.cc \
Base/ShowerVeto.h Base/ShowerVeto.fh Base/ShowerVeto.cc \
Base/FullShowerVeto.h Base/FullShowerVeto.fh Base/FullShowerVeto.cc \
SplittingFunctions/SplittingGenerator.cc SplittingFunctions/SplittingGenerator.h\
SplittingFunctions/SplittingGenerator.fh \
Base/ShowerTree.h Base/ShowerTree.fh Base/ShowerTree.cc \
ShowerConfig.h ShowerConfig.cc \
Base/Branching.h \
Base/ShowerParticle.cc Base/ShowerParticle.fh Base/ShowerParticle.h \
Base/ShowerKinematics.fh Base/ShowerKinematics.h Base/ShowerKinematics.cc \
Base/ShowerBasis.fh Base/ShowerBasis.h Base/ShowerBasis.cc \
Base/ShowerProgenitor.fh Base/ShowerProgenitor.h \
Base/SudakovFormFactor.cc Base/SudakovFormFactor.h Base/SudakovFormFactor.fh \
SplittingFunctions/SplittingFunction.h SplittingFunctions/SplittingFunction.fh \
SplittingFunctions/SplittingFunction.cc \
Base/ShowerVertex.cc Base/ShowerVertex.fh Base/ShowerVertex.h
diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in
--- a/src/defaults/Shower.in
+++ b/src/defaults/Shower.in
@@ -1,310 +1,310 @@
# -*- ThePEG-repository -*-
############################################################
# Setup of default parton shower
#
# Useful switches for users are marked near the top of
# this file.
#
# Don't edit this file directly, but reset the switches
# in your own input files!
############################################################
library HwMPI.so
library HwShower.so
library HwMatching.so
mkdir /Herwig/Shower
cd /Herwig/Shower
create Herwig::QTildeShowerHandler ShowerHandler
newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
# use LO PDFs for Shower, can be changed later
newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
#####################################
# initial setup, don't change these!
#####################################
create Herwig::SplittingGenerator SplittingGenerator
create Herwig::ShowerAlphaQCD AlphaQCD
create Herwig::ShowerAlphaQED AlphaQED
-create Herwig::QTildeFinder PartnerFinder
+create Herwig::PartnerFinder PartnerFinder
newdef PartnerFinder:PartnerMethod 1
newdef PartnerFinder:ScaleChoice 1
create Herwig::QTildeReconstructor KinematicsReconstructor
newdef KinematicsReconstructor:ReconstructionOption Colour3
newdef KinematicsReconstructor:InitialStateReconOption SofterFraction
newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost
newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD
newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED
newdef ShowerHandler:PartnerFinder PartnerFinder
newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor
newdef ShowerHandler:SplittingGenerator SplittingGenerator
newdef ShowerHandler:SpinCorrelations Yes
newdef ShowerHandler:SoftCorrelations Singular
##################################################################
# Intrinsic pT
#
# Recommended:
# 1.9 GeV for Tevatron W/Z production.
# 2.1 GeV for LHC W/Z production at 10 TeV
# 2.2 GeV for LHC W/Z production at 14 TeV
#
# Set all parameters to 0 to disable
##################################################################
newdef ShowerHandler:IntrinsicPtGaussian 1.3*GeV
newdef ShowerHandler:IntrinsicPtBeta 0
newdef ShowerHandler:IntrinsicPtGamma 0*GeV
newdef ShowerHandler:IntrinsicPtIptmax 0*GeV
#############################################################
# Set up truncated shower handler.
#############################################################
create Herwig::PowhegShowerHandler PowhegShowerHandler
set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PartnerFinder PartnerFinder
newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor
newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator
newdef PowhegShowerHandler:Interactions QCDandQED
newdef PowhegShowerHandler:SpinCorrelations Yes
newdef PowhegShowerHandler:SoftCorrelations Singular
newdef PowhegShowerHandler:IntrinsicPtGaussian 1.3*GeV
newdef PowhegShowerHandler:IntrinsicPtBeta 0
newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV
newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV
newdef PowhegShowerHandler:ReconstructionOption OffShell5
#############################################################
# End of interesting user servicable section.
#
# Anything that follows below should only be touched if you
# know what you're doing.
#
# Really.
#############################################################
#
# a few default values
newdef ShowerHandler:MECorrMode 1
newdef ShowerHandler:ReconstructionOption OffShell5
newdef AlphaQCD:ScaleFactor 1.0
newdef AlphaQCD:NPAlphaS 2
newdef AlphaQCD:Qmin 0.935
newdef AlphaQCD:NumberOfLoops 2
newdef AlphaQCD:InputOption 1
newdef AlphaQCD:AlphaMZ 0.126234
#
#
# Lets set up all the splittings
create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn
set QtoQGammaSplitFn:InteractionType QED
set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral
set QtoQGammaSplitFn:AngularOrdered Yes
create Herwig::HalfHalfOneSplitFn QtoQGSplitFn
newdef QtoQGSplitFn:InteractionType QCD
newdef QtoQGSplitFn:ColourStructure TripletTripletOctet
set QtoQGSplitFn:AngularOrdered Yes
create Herwig::OneOneOneSplitFn GtoGGSplitFn
newdef GtoGGSplitFn:InteractionType QCD
newdef GtoGGSplitFn:ColourStructure OctetOctetOctet
set GtoGGSplitFn:AngularOrdered Yes
create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn
newdef WtoWGammaSplitFn:InteractionType QED
newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral
set WtoWGammaSplitFn:AngularOrdered Yes
create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn
newdef GtoQQbarSplitFn:InteractionType QCD
newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet
set GtoQQbarSplitFn:AngularOrdered Yes
create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn
newdef GammatoQQbarSplitFn:InteractionType QED
newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged
set GammatoQQbarSplitFn:AngularOrdered Yes
create Herwig::HalfOneHalfSplitFn QtoGQSplitFn
newdef QtoGQSplitFn:InteractionType QCD
newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet
set QtoGQSplitFn:AngularOrdered Yes
create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn
newdef QtoGammaQSplitFn:InteractionType QED
newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged
set QtoGammaQSplitFn:AngularOrdered Yes
#
# Now the Sudakovs
create Herwig::QTildeSudakov SudakovCommon
newdef SudakovCommon:Alpha AlphaQCD
newdef SudakovCommon:cutoffKinScale 0.0*GeV
newdef SudakovCommon:PDFmax 1.0
newdef SudakovCommon:CutOffOption pT
newdef SudakovCommon:pTmin 1.222798*GeV
cp SudakovCommon QtoQGSudakov
newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn
newdef QtoQGSudakov:PDFmax 1.9
cp SudakovCommon QtoQGammaSudakov
set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn
set QtoQGammaSudakov:Alpha AlphaQED
set QtoQGammaSudakov:PDFmax 1.9
cp QtoQGammaSudakov LtoLGammaSudakov
# Technical parameter to stop evolution.
set LtoLGammaSudakov:pTmin 0.000001
cp SudakovCommon GtoGGSudakov
newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn
newdef GtoGGSudakov:PDFmax 2.0
cp SudakovCommon WtoWGammaSudakov
newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn
set WtoWGammaSudakov:Alpha AlphaQED
cp SudakovCommon GtoQQbarSudakov
newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtoQQbarSudakov:PDFmax 120.0
cp SudakovCommon GammatoQQbarSudakov
newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn
set GammatoQQbarSudakov:Alpha AlphaQED
newdef GammatoQQbarSudakov:PDFmax 120.0
cp SudakovCommon GtobbbarSudakov
newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtobbbarSudakov:PDFmax 40000.0
cp SudakovCommon GtoccbarSudakov
newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtoccbarSudakov:PDFmax 2000.0
cp SudakovCommon QtoGQSudakov
newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn
cp SudakovCommon QtoGammaQSudakov
newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn
set QtoGammaQSudakov:Alpha AlphaQED
cp SudakovCommon utoGuSudakov
newdef utoGuSudakov:SplittingFunction QtoGQSplitFn
newdef utoGuSudakov:PDFFactor OverOneMinusZ
newdef utoGuSudakov:PDFmax 5.0
cp SudakovCommon dtoGdSudakov
newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn
newdef dtoGdSudakov:PDFFactor OverOneMinusZ
#
# Now add the final splittings
#
do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov
#
do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov
#
do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov
do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov
do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov
do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov
do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov
do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov
#
do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov
#
do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov
#
# Now lets add the initial splittings. Remember the form a->b,c; means
# that the current particle b is given and we backward branch to new
# particle a which is initial state and new particle c which is final state
#
do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov
#
do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov
do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov
#
do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov
#
do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov
do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov
do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov
do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov
do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov
#
do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:27 PM (9 h, 25 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023638
Default Alt Text
(99 KB)

Event Timeline