Page MenuHomeHEPForge
Files F8723768

No OneTemporary
Actions

View Options

diff --git a/Shower/Dipole/Kernels/FFMgx2ggxDipoleKernel.cc b/Shower/Dipole/Kernels/FFMgx2ggxDipoleKernel.cc
--- a/Shower/Dipole/Kernels/FFMgx2ggxDipoleKernel.cc
+++ b/Shower/Dipole/Kernels/FFMgx2ggxDipoleKernel.cc
@@ -1,146 +1,142 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FFMgx2ggxDipoleKernel class.
//
#include "FFMgx2ggxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FFMgx2ggxDipoleKernel::FFMgx2ggxDipoleKernel()
: DipoleSplittingKernel(){}
FFMgx2ggxDipoleKernel::~FFMgx2ggxDipoleKernel() {}
IBPtr FFMgx2ggxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr FFMgx2ggxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool FFMgx2ggxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
useThisKernel() &&
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() != ZERO &&
!ind.initialStateEmitter() && !ind.initialStateSpectator() &&
!ind.incomingDecayEmitter() && !ind.incomingDecaySpectator();
}
bool FFMgx2ggxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emitter(b)->id() == ParticleID::g &&
sk.emission(b)->id() == ParticleID::g &&
abs(spectator(a)->mass()) == abs(sk.spectator(b)->mass());
}
tcPDPtr FFMgx2ggxDipoleKernel::emitter(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr FFMgx2ggxDipoleKernel::emission(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr FFMgx2ggxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double FFMgx2ggxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
// These are the physical variables as used in the
// standard form of the kernel (i.e. do not redefine variables or kernel)
const double z = split.lastZ();
const Energy pt = split.lastPt();
// Need zPrime to calculate y,
// TODO: Should just store y in the dipole splitting info everywhere anyway!!!
// The only value stored in dInfo.lastSplittingParameters() should be zPrime
const double zPrime = split.lastSplittingParameters()[0];
// Construct mass squared variables
// Construct mass squared variables
double muj2 = sqr(split.spectatorMass() / split.scale());
double bar = 1. - muj2;
// Calculate y
double y = sqr(pt)/sqr(split.scale()) / (bar*zPrime*(1.-zPrime));
double vijk = sqrt( sqr(2.*muj2+bar*(1.-y))-4.*muj2 ) / (bar*(1.-y));
double viji = 1.;
double zp = 0.5*(1.+viji*vijk);
double zm = 0.5*(1.-viji*vijk);
// how to choose kappa?
double kappa = 0.;
double S1=1./(1.-z*(1.-y));
double S2=1./(1.-(1.-z)*(1.-y));
double NS=(z*(1.-z)-(1.-kappa)*zp*zm-2.)/vijk;
if( theAsymmetryOption == 0 ){
ret *= 3.*( S1 + 0.5 * NS);
}else if ( theAsymmetryOption == 1 ){
ret *= 3.*z*( S1 +S2 + NS );
}else{
ret *= 3.*0.5*( S1 + S2 + NS );
}
-
- ret *= 3.*(1./(1.-z*(1.-y))+ 0.5*(z*(1.-z)-(1.-kappa)*zp*zm-2.)/vijk);
- return ret > 0. ? ret : 0.;
-
+ return ret > 0. ? ret : 0.;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFMgx2ggxDipoleKernel::persistentOutput(PersistentOStream & os) const {
os << theAsymmetryOption;
}
void FFMgx2ggxDipoleKernel::persistentInput(PersistentIStream & is, int) {
is >> theAsymmetryOption;
}
ClassDescription<FFMgx2ggxDipoleKernel> FFMgx2ggxDipoleKernel::initFFMgx2ggxDipoleKernel;
// Definition of the static class description member.
void FFMgx2ggxDipoleKernel::Init() {
static ClassDocumentation<FFMgx2ggxDipoleKernel> documentation
("FFMgx2ggxDipoleKernel");
static Parameter<FFMgx2ggxDipoleKernel,int> interfacetheAsymmetryOption
("AsymmetryOption",
"The asymmetry option for final state gluon spliitings.",
&FFMgx2ggxDipoleKernel::theAsymmetryOption, 1, 0, 0,
false, false, Interface::lowerlim);
}
-
diff --git a/Shower/Dipole/Kernels/FFgx2ggxDipoleKernel.cc b/Shower/Dipole/Kernels/FFgx2ggxDipoleKernel.cc
--- a/Shower/Dipole/Kernels/FFgx2ggxDipoleKernel.cc
+++ b/Shower/Dipole/Kernels/FFgx2ggxDipoleKernel.cc
@@ -1,115 +1,113 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FFgx2ggxDipoleKernel class.
//
#include "FFgx2ggxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FFgx2ggxDipoleKernel::FFgx2ggxDipoleKernel()
: DipoleSplittingKernel(){}
FFgx2ggxDipoleKernel::~FFgx2ggxDipoleKernel() {}
IBPtr FFgx2ggxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr FFgx2ggxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool FFgx2ggxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
useThisKernel() &&
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() == ZERO &&
!ind.initialStateEmitter() && !ind.initialStateSpectator();
}
bool FFgx2ggxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emitter(b)->id() == ParticleID::g &&
sk.emission(b)->id() == ParticleID::g;
}
tcPDPtr FFgx2ggxDipoleKernel::emitter(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr FFgx2ggxDipoleKernel::emission(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr FFgx2ggxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double FFgx2ggxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
double z = split.lastZ();
double y = sqr(split.lastPt() / split.scale()) / (z*(1.-z));
double S1=1./(1.-z*(1.-y));
double S2=1./(1.-(1.-z)*(1.-y));
double NS=(-2 + z*(1.-z));
if( theAsymmetryOption == 0 ){
ret *= 3.*( S1 + 0.5 * NS);
}else if ( theAsymmetryOption == 1 ){
ret *= 3.*z*( S1 +S2 + NS );
}else{
ret *= 3.*0.5*( S1 + S2 + NS );
}
return ret > 0. ? ret : 0.;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFgx2ggxDipoleKernel::persistentOutput(PersistentOStream & os) const {
-
os << theAsymmetryOption;
}
void FFgx2ggxDipoleKernel::persistentInput(PersistentIStream & is, int) {
is >> theAsymmetryOption;
}
ClassDescription<FFgx2ggxDipoleKernel> FFgx2ggxDipoleKernel::initFFgx2ggxDipoleKernel;
// Definition of the static class description member.
void FFgx2ggxDipoleKernel::Init() {
static ClassDocumentation<FFgx2ggxDipoleKernel> documentation
("FFgx2ggxDipoleKernel");
static Parameter<FFgx2ggxDipoleKernel,int> interfacetheAsymmetryOption
("AsymmetryOption",
"The asymmetry option for final state gluon spliitings.",
&FFgx2ggxDipoleKernel::theAsymmetryOption, 1, 0, 0,
false, false, Interface::lowerlim);
}
-
diff --git a/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.cc b/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.cc
--- a/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.cc
+++ b/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.cc
@@ -1,155 +1,159 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FIMDecaygx2ggxDipoleKernel class.
//
#include "FIMDecaygx2ggxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FIMDecaygx2ggxDipoleKernel::FIMDecaygx2ggxDipoleKernel()
- : DipoleSplittingKernel(),theSymmetryFactor(0.5){}
+ : DipoleSplittingKernel(){}
FIMDecaygx2ggxDipoleKernel::~FIMDecaygx2ggxDipoleKernel() {}
IBPtr FIMDecaygx2ggxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr FIMDecaygx2ggxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool FIMDecaygx2ggxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
useThisKernel() &&
ind.incomingDecaySpectator() && !ind.incomingDecayEmitter() &&
ind.emitterData()->id() == ParticleID::g &&
!(ind.spectatorData()->mass() == ZERO) &&
// Initial state here refers to the entire event
!ind.initialStateEmitter() && !ind.initialStateSpectator();
}
bool FIMDecaygx2ggxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emission(b)->id() == ParticleID::g &&
sk.emitter(b)->id() == ParticleID::g &&
abs(sk.spectator(b)->mass()) == abs(spectator(a)->mass());
}
tcPDPtr FIMDecaygx2ggxDipoleKernel::emitter(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr FIMDecaygx2ggxDipoleKernel::emission(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr FIMDecaygx2ggxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double FIMDecaygx2ggxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
// These are the physical variables as used in the standard form of the kernel (i.e. do not redefine variables or kernel)
double z = split.lastZ();
Energy pt = split.lastPt();
// Need zPrime to calculate y,
// TODO: Should just store y in the dipole splitting info everywhere anyway!!!
// The only value stored in dInfo.lastSplittingParameters() should be zPrime
//assert(split.lastSplittingParameters().size() == 1 );
double zPrime = split.lastSplittingParameters()[0];
// Construct mass squared variables
- double mua2 = sqr( split.spectatorData()->mass() / split.scale() );
+ double mua2 = sqr( split.spectatorMass() / split.scale() );
// Recoil system mass
double muj2 = sqr(split.recoilMass() / split.scale());
double bar = 1. - muj2;
// Calculate y
double y = (sqr(pt)/sqr(split.scale())) / (bar*zPrime*(1.-zPrime));
if( sqr(2.*muj2+bar*(1.-y))-4.*muj2 < 0. ){
generator()->logWarning( Exception()
<< "error in FIMDecaygx2ggxDipoleKernel::evaluate -- " <<
"muj2 " << muj2 << " y " << y << Exception::warning );
return 0.0;
}
double vijk = sqrt( sqr(2.*muj2+bar*(1.-y))-4.*muj2 ) / (bar*(1.-y));
double viji = 1.;
double vbar = 1.;
double zp = 0.5*(1.+viji*vijk);
double zm = 0.5*(1.-viji*vijk);
// how to choose kappa?
double kappa = 0.;
- ret *=
- theSymmetryFactor
- * 3.*( (2.*y + 1.)/((1.+y)-z*(1.-y)) + (2.*y + 1.)/((1.+y)-(1.-z)*(1.-y))
- + (1./vijk)*( z*(1.-z) - (1.-kappa)*zp*zm - 2. ) )
- +
- (!strictLargeN() ? 4./3. : 3./2.)
- * (
- y/(1.-z*(1.-y)) * ( 2.*(2.*y + 1.)/((1.+y)-z*(1.-y))
- - (vbar/vijk)*(2. + 2.*mua2/((1.-z*(1.-y))*bar)) )
- +
- y/(1.-(1.-z)*(1.-y)) * ( 2.*(2.*y + 1.)/((1.+y)-(1.-z)*(1.-y))
- - (vbar/vijk)*(2. + 2.*mua2/((1.-(1.-z)*(1.-y))*bar)) )
- );
-
+ double S1 = 0.5*3.*(2.*y + 1.)/((1.+y)-z*(1.-y)) +
+ (!strictLargeN() ? 4./3. : 3./2.)*
+ y/(1.-z*(1.-y)) * ( 2.*(2.*y + 1.)/((1.+y)-z*(1.-y))
+ - (vbar/vijk)*(2. + 2.*mua2/((1.-z*(1.-y))*bar)) );
+ double S2 = 0.5*3.*(2.*y + 1.)/((1.+y)-(1.-z)*(1.-y)) +
+ (!strictLargeN() ? 4./3. : 3./2.)*
+ y/(1.-(1.-z)*(1.-y)) * ( 2.*(2.*y + 1.)/((1.+y)-(1.-z)*(1.-y))
+ - (vbar/vijk)*(2. + 2.*mua2/((1.-(1.-z)*(1.-y))*bar)) );
+ double NS = 0.5*3.*(z*(1.-z)-(1.-kappa)*zp*zm - 2.)/vijk;
+
+
+ if( theAsymmetryOption == 0 ){
+ ret *= 2.*S1 + NS;
+ }else if ( theAsymmetryOption == 1 ){
+ ret *= 2.*z*( S1 + S2 + NS );
+ }else{
+ ret *= S1 + S2 + NS;
+ }
+
return ret > 0. ? ret : 0.;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMDecaygx2ggxDipoleKernel::persistentOutput(PersistentOStream & os) const {
-
- os<<theSymmetryFactor;
+ os<<theAsymmetryOption;
}
void FIMDecaygx2ggxDipoleKernel::persistentInput(PersistentIStream & is, int) {
-
- is>>theSymmetryFactor;
+ is>>theAsymmetryOption;
}
ClassDescription<FIMDecaygx2ggxDipoleKernel> FIMDecaygx2ggxDipoleKernel::initFIMDecaygx2ggxDipoleKernel;
// Definition of the static class description member.
void FIMDecaygx2ggxDipoleKernel::Init() {
static ClassDocumentation<FIMDecaygx2ggxDipoleKernel> documentation
("FIMDecaygx2ggxDipoleKernel");
- static Parameter<FIMDecaygx2ggxDipoleKernel,double> interfaceSymmetryFactor
- ("SymmetryFactor",
- "The symmetry factor for final state gluon splittings.",
- &FIMDecaygx2ggxDipoleKernel::theSymmetryFactor, 1.0, 0.0, 0,
+ static Parameter<FIMDecaygx2ggxDipoleKernel,int> interfacetheAsymmetryOption
+ ("AsymmetryOption",
+ "The asymmetry option for final state gluon spliitings.",
+ &FIMDecaygx2ggxDipoleKernel::theAsymmetryOption, 0, 0, 0,
false, false, Interface::lowerlim);
+
}
diff --git a/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.h b/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.h
--- a/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.h
+++ b/Shower/Dipole/Kernels/FIMDecaygx2ggxDipoleKernel.h
@@ -1,188 +1,186 @@
// -*- C++ -*-
#ifndef HERWIG_FIMDecaygx2ggxDipoleKernel_H
#define HERWIG_FIMDecaygx2ggxDipoleKernel_H
//
// This is the declaration of the FIMDecaygx2ggxDipoleKernel class.
//
#include "DipoleSplittingKernel.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Stephen Webster
*
* \brief FIMDecaygx2ggxDipoleKernel implements the g -> gg
* splitting off a final-initial decay dipole and includes the
* contribution from the splitting of the intial / decay particle
*
*/
class FIMDecaygx2ggxDipoleKernel: public DipoleSplittingKernel {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FIMDecaygx2ggxDipoleKernel();
/**
* The destructor.
*/
virtual ~FIMDecaygx2ggxDipoleKernel();
//@}
public:
/**
* Return true, if this splitting kernel
* applies to the given dipole index.
*/
virtual bool canHandle(const DipoleIndex&) const;
/**
* Return true, if this splitting kernel is
* the same for the given index a, as the given
* splitting kernel for index b.
*/
virtual bool canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const;
/**
* Return the emitter data after splitting, given
* a dipole index.
*/
virtual tcPDPtr emitter(const DipoleIndex&) const;
/**
* Return the emission data after splitting, given
* a dipole index.
*/
virtual tcPDPtr emission(const DipoleIndex&) const;
/**
* Return the spectator data after splitting, given
* a dipole index.
*/
virtual tcPDPtr spectator(const DipoleIndex&) const;
/**
* Evaluate this splitting kernel for the given
* dipole splitting.
*/
virtual double evaluate(const DipoleSplittingInfo&) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
- * Symmetry factor for final state gluon splittings (should be 1/2).
- */
-
- double theSymmetryFactor;
-
- /**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<FIMDecaygx2ggxDipoleKernel> initFIMDecaygx2ggxDipoleKernel;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FIMDecaygx2ggxDipoleKernel & operator=(const FIMDecaygx2ggxDipoleKernel &);
+ /**
+ * Asymmetry option for final state gluon splittings.
+ */
+ int theAsymmetryOption=0;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of FIMDecaygx2ggxDipoleKernel. */
template <>
struct BaseClassTrait<Herwig::FIMDecaygx2ggxDipoleKernel,1> {
/** Typedef of the first base class of FIMDecaygx2ggxDipoleKernel. */
typedef Herwig::DipoleSplittingKernel NthBase;
};
/** This template specialization informs ThePEG about the name of
* the FIMDecaygx2ggxDipoleKernel class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::FIMDecaygx2ggxDipoleKernel>
: public ClassTraitsBase<Herwig::FIMDecaygx2ggxDipoleKernel> {
/** Return a platform-independent class name */
static string className() { return "Herwig::FIMDecaygx2ggxDipoleKernel"; }
/**
* The name of a file containing the dynamic library where the class
* FIMDecaygx2ggxDipoleKernel is implemented. It may also include several, space-separated,
* libraries if the class FIMDecaygx2ggxDipoleKernel depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_FIMDecaygx2ggxDipoleKernel_H */
diff --git a/Shower/Dipole/Kernels/FIMDecaygx2qqxDipoleKernel.cc b/Shower/Dipole/Kernels/FIMDecaygx2qqxDipoleKernel.cc
--- a/Shower/Dipole/Kernels/FIMDecaygx2qqxDipoleKernel.cc
+++ b/Shower/Dipole/Kernels/FIMDecaygx2qqxDipoleKernel.cc
@@ -1,142 +1,142 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FIMDecaygx2qqxDipoleKernel class.
//
#include "FIMDecaygx2qqxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FIMDecaygx2qqxDipoleKernel::FIMDecaygx2qqxDipoleKernel()
: DipoleSplittingKernel() {}
FIMDecaygx2qqxDipoleKernel::~FIMDecaygx2qqxDipoleKernel() {}
IBPtr FIMDecaygx2qqxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr FIMDecaygx2qqxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool FIMDecaygx2qqxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
useThisKernel() &&
ind.incomingDecaySpectator() && !ind.incomingDecayEmitter() &&
ind.emitterData()->id() == ParticleID::g &&
!(ind.spectatorData()->mass() == ZERO) &&
// Initial state here refers to the entire event
!ind.initialStateEmitter() && !ind.initialStateSpectator();
}
bool FIMDecaygx2qqxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emitter(b)->id() + sk.emission(b)->id() == 0 &&
abs(sk.emitter(b)->id()) < 6 &&
emitter(a)->id() == sk.emitter(b)->id() &&
abs(sk.spectator(b)->mass()) == abs(spectator(a)->mass());
}
tcPDPtr FIMDecaygx2qqxDipoleKernel::emitter(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6);
return flavour();
}
tcPDPtr FIMDecaygx2qqxDipoleKernel::emission(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6);
return flavour()->CC();
}
tcPDPtr FIMDecaygx2qqxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double FIMDecaygx2qqxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
// These are the physical variables as used in the standard form of the kernel (i.e. do not redefine variables or kernel)
double z = split.lastZ();
Energy pt = split.lastPt();
// Need zPrime to calculate y,
// TODO: Should just store y in the dipole splitting info everywhere anyway!!!
// The only value stored in dInfo.lastSplittingParameters() should be zPrime
//assert(split.lastSplittingParameters().size() == 1 );
double zPrime = split.lastSplittingParameters()[0];
// Construct mass squared variables
- double mua2 = sqr( split.spectatorData()->mass() / split.scale() );
double mui2 = sqr(split.emitterData()->mass() / split.scale());
double mu2 = mui2;
+ //double mua2 = sqr( split.spectatorMass() / split.scale() );
// Recoil system mass
double muj2 = sqr(split.recoilMass() / split.scale());
double bar = 1. - mui2 - mu2 - muj2;
// Calculate y
double y = (sqr(pt)/sqr(split.scale()) + sqr(1.-zPrime)*mui2 + sqr(zPrime)*mu2) / (bar*zPrime*(1.-zPrime));
if( sqr(2.*muj2+bar*(1.-y))-4.*muj2 < 0. ){
generator()->logWarning( Exception()
<< "error in FIMDecaygx2qqxDipoleKernel::evaluate -- " <<
"muj2 " << muj2 << " mui2 " << mui2 << " y " << y << Exception::warning );
return 0.0;
}
double vijk = sqrt( sqr(2.*muj2+bar*(1.-y))-4.*muj2 ) / (bar*(1.-y));
double viji = sqrt( sqr(bar*y)-4.*sqr(mui2) ) / (bar*y+2.*mui2);
- double vbar = sqrt( 1.+sqr(mui2)+sqr(muj2)-2.*(mui2+muj2+mui2*muj2) ) / bar;
+ //double vbar = sqrt( 1.+sqr(mui2)+sqr(muj2)-2.*(mui2+muj2+mui2*muj2) ) / bar;
double zp = 0.5*(1.+viji*vijk);
double zm = 0.5*(1.-viji*vijk);
// how to choose kappa?
double kappa = 0.;
ret *= 0.25 / vijk * ( 1. - 2.*( z*(1.-z) - (1.-kappa)*zp*zm - kappa*mui2/(2*mui2+bar*y) ) );
return ret > 0. ? ret : 0.;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMDecaygx2qqxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void FIMDecaygx2qqxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FIMDecaygx2qqxDipoleKernel> FIMDecaygx2qqxDipoleKernel::initFIMDecaygx2qqxDipoleKernel;
// Definition of the static class description member.
void FIMDecaygx2qqxDipoleKernel::Init() {
static ClassDocumentation<FIMDecaygx2qqxDipoleKernel> documentation
("FIMDecaygx2qqxDipoleKernel");
}
diff --git a/Shower/Dipole/Kernels/FIMDecayqx2qgxDipoleKernel.cc b/Shower/Dipole/Kernels/FIMDecayqx2qgxDipoleKernel.cc
--- a/Shower/Dipole/Kernels/FIMDecayqx2qgxDipoleKernel.cc
+++ b/Shower/Dipole/Kernels/FIMDecayqx2qgxDipoleKernel.cc
@@ -1,141 +1,144 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FIMDecayqx2qgxDipoleKernel class.
//
#include "FIMDecayqx2qgxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FIMDecayqx2qgxDipoleKernel::FIMDecayqx2qgxDipoleKernel() : DipoleSplittingKernel() {}
FIMDecayqx2qgxDipoleKernel::~FIMDecayqx2qgxDipoleKernel() {}
IBPtr FIMDecayqx2qgxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr FIMDecayqx2qgxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool FIMDecayqx2qgxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
useThisKernel() &&
ind.incomingDecaySpectator() && !ind.incomingDecayEmitter() &&
abs(ind.emitterData()->id()) < 7 &&
// This line matches to the kernel declared in a .in file for the given emitter flavour
abs(ind.emitterData()->id()) == abs(flavour()->id()) &&
!(ind.spectatorData()->mass() == ZERO) &&
// Initial state here refers to the entire event
!ind.initialStateEmitter() && !ind.initialStateSpectator();
}
bool FIMDecayqx2qgxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emission(b)->id() == ParticleID::g &&
abs(sk.emitter(b)->id()) < 7 &&
abs(sk.emitter(b)->mass()) == abs(emitter(a)->mass()) &&
abs(sk.spectator(b)->mass()) == abs(spectator(a)->mass());
}
tcPDPtr FIMDecayqx2qgxDipoleKernel::emitter(const DipoleIndex& ind) const {
assert(flavour());
assert(abs(flavour()->id()) < 7);
return ind.emitterData()->id() > 0 ?
(tcPDPtr) flavour() : (tcPDPtr) flavour()->CC();
}
tcPDPtr FIMDecayqx2qgxDipoleKernel::emission(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr FIMDecayqx2qgxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double FIMDecayqx2qgxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
// These are the physical variables as used in the standard form of the kernel (i.e. do not redefine variables or kernel)
double z = split.lastZ();
Energy pt = split.lastPt();
// Need zPrime to calculate y,
// TODO: Should just store y in the dipole splitting info everywhere anyway!!!
// The only value stored in dInfo.lastSplittingParameters() should be zPrime
//assert(split.lastSplittingParameters().size() == 1 );
double zPrime = split.lastSplittingParameters()[0];
// Construct mass squared variables
- double mui2 = sqr(split.emitterData()->mass() / split.scale());
+ // Note for q->qg can use the emitterMass
+ // (i.e. mass of emitter before splitting = mass of emitter after)
+ double mui2 = sqr(split.emitterMass() / split.scale());
// Recoil system mass
double muj2 = sqr(split.recoilMass() / split.scale());
- double mua2 = sqr( split.spectatorData()->mass() / split.scale() );
+ // This should be equal to one
+ double mua2 = sqr( split.spectatorMass() / split.scale() );
double bar = 1. - mui2 - muj2;
// Calculate y
double y = (sqr(pt)/sqr(split.scale()) + sqr(1.-zPrime)*mui2) / (bar*zPrime*(1.-zPrime));
if( sqr(2.*muj2+bar*(1.-y))-4.*muj2 < 0. ){
generator()->logWarning( Exception()
<< "error in FIMDecayqx2qgxDipoleKernel::evaluate -- " <<
"muj2 " << muj2 << " mui2 " << mui2 << " y " << y << Exception::warning );
return 0.0;
}
double vijk = sqrt( sqr(2.*muj2 + bar*(1.-y))-4.*muj2 ) / (bar*(1.-y));
double vbar = sqrt( 1.+sqr(mui2)+sqr(muj2)-2.*(mui2+muj2+mui2*muj2) ) / bar;
ret *=
(!strictLargeN() ? 4./3. : 3./2.)
* ( ( 2.*(2.*mui2/bar + 2.*y + 1.)/((1.+y)-z*(1.-y))
- (vbar/vijk)*((1.+z) + 2.*mui2/(y*bar)) )
+ y/(1.-z*(1.-y)) * ( 2.*(2.*mui2/bar + 2.*y + 1.)/((1.+y)-z*(1.-y))
- (vbar/vijk)*(2. + 2.*mua2/((1.-z*(1.-y))*bar)) ) );
return ret > 0. ? ret : 0.;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMDecayqx2qgxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void FIMDecayqx2qgxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FIMDecayqx2qgxDipoleKernel> FIMDecayqx2qgxDipoleKernel::initFIMDecayqx2qgxDipoleKernel;
// Definition of the static class description member.
void FIMDecayqx2qgxDipoleKernel::Init() {
static ClassDocumentation<FIMDecayqx2qgxDipoleKernel> documentation
("FIMDecayqx2qgxDipoleKernel");
}
diff --git a/Shower/Dipole/Kinematics/FIMassiveDecayKinematics.cc b/Shower/Dipole/Kinematics/FIMassiveDecayKinematics.cc
--- a/Shower/Dipole/Kinematics/FIMassiveDecayKinematics.cc
+++ b/Shower/Dipole/Kinematics/FIMassiveDecayKinematics.cc
@@ -1,504 +1,512 @@
// -*- C++ -*-
//
// FIMassiveDecayKinematics.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 FIMassiveDecayKinematics class.
//
#include "FIMassiveDecayKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h"
#include "Herwig/Shower/Dipole/Kernels/DipoleSplittingKernel.h"
#include "ThePEG/Interface/Switch.h"
using namespace Herwig;
FIMassiveDecayKinematics::FIMassiveDecayKinematics()
: DipoleSplittingKinematics(), theFullJacobian(true) {}
FIMassiveDecayKinematics::~FIMassiveDecayKinematics() {}
IBPtr FIMassiveDecayKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FIMassiveDecayKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> FIMassiveDecayKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return {0.0,1.0};
}
pair<double,double> FIMassiveDecayKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emissionData()->id() != ParticleID::g ){
return {0.5*(1.-c),0.5*(1.+c)};
}
double b = log((1.+c)/(1.-c));
return {-b,b};
}
return {-log(0.5*(1.+c)),-log(0.5*(1.-c))};
}
Energy FIMassiveDecayKinematics::dipoleScale(const Lorentz5Momentum&,
const Lorentz5Momentum& pSpectator) const {
return pSpectator.m();
}
Energy FIMassiveDecayKinematics::recoilMassKin(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
Lorentz5Momentum pk = pSpectator - pEmitter;
Energy pkmass = pk.m();
return pkmass;
}
Energy FIMassiveDecayKinematics::ptMax(Energy dScale,
double, double,
const DipoleSplittingInfo& dInfo,
const DipoleSplittingKernel& split) const {
DipoleIndex ind = dInfo.index();
double mui2 = 0.;
// g->gg and g->qqbar
if ( abs(split.emitter(ind)->id()) == abs(split.emission(ind)->id()) ) {
mui2 = sqr(split.emitter(ind)->mass() / dScale);
}
// Otherwise have X->Xg (should work for SUSY)
else {
mui2 = sqr(dInfo.emitterMass()/dScale);
}
double mu2 = sqr(split.emission(ind)->mass() / dScale);
// Mass of recoil system
double muj = dInfo.recoilMass() / dScale;
return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.*muj ) * dScale;
}
Energy FIMassiveDecayKinematics::ptMax(Energy dScale,
double, double,
const DipoleIndex& ind,
const DipoleSplittingKernel& split,
tPPtr emitter, tPPtr spectator) const {
double mui2 = 0.;
// g->gg and g->qqbar
if ( abs(split.emitter(ind)->id()) == abs(split.emission(ind)->id()) ) {
mui2 = sqr(split.emitter(ind)->mass() / dScale);
}
// Otherwise have X->Xg (should work for SUSY)
else {
mui2 = sqr(emitter->mass()/dScale);
}
double mu2 = sqr(split.emission(ind)->mass() / dScale);
// Mass of recoil system
double muj = recoilMassKin(emitter->momentum(),spectator->momentum()) / dScale;
return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.*muj ) * dScale;
}
Energy FIMassiveDecayKinematics::QMax(Energy dScale,
double, double,
const DipoleSplittingInfo& dInfo,
const DipoleSplittingKernel&) const {
assert(false && "implementation missing");
// Mass of recoil system
double Muj2 = sqr( dInfo.recoilMass() / dScale );
double Muj = sqrt( Muj2 );
return dScale * ( 1.-2.*Muj+Muj2 );
}
// The name of this function is misleading
// scale here is defined as sqr(scale) = sqr(qi+qj)
// Here, scale is Q
Energy FIMassiveDecayKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
double zPrime = split.lastSplittingParameters()[0];
// masses
Energy2 mi2 = ZERO;
// g->gg and g->qqbar
if ( abs(split.emitterData()->id()) == abs(split.emissionData()->id()) ) {
mi2 = sqr( split.emitterData()->mass());
}
// Otherwise have X->Xg (should work for SUSY)
else {
mi2 = sqr(split.emitterMass());
}
Energy2 m2 = sqr(split.emissionData()->mass());
Energy2 pt2 = zPrime*(1.-zPrime)*sqr(scale) - (1-zPrime)*mi2 - zPrime*m2;
assert(pt2 >= ZERO);
return sqrt(pt2);
}
// This is simply the inverse of PtFromQ
Energy FIMassiveDecayKinematics::QFromPt(Energy pt, const DipoleSplittingInfo& split) const {
double zPrime = split.lastSplittingParameters()[0];
// masses
Energy2 mi2 = ZERO;
// g->gg and g->qqbar
if ( abs(split.emitterData()->id()) == abs(split.emissionData()->id()) ) {
mi2 = sqr( split.emitterData()->mass());
}
// Otherwise have X->Xg (should work for SUSY)
else {
mi2 = sqr(split.emitterMass());
}
Energy2 m2 = sqr(split.emissionData()->mass());
Energy2 Q2 = (sqr(pt) + (1-zPrime)*mi2 + zPrime*m2)/(zPrime*(1.-zPrime));
return sqrt(Q2);
}
double FIMassiveDecayKinematics::ptToRandom(Energy pt, Energy,
double,double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
bool FIMassiveDecayKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info,
const DipoleSplittingKernel&) {
// scaled masses
double Mui2 = sqr(info.emitterMass() / info.scale());
double Muj2 = sqr(info.recoilMass() / info.scale());
double muj2 = Muj2;
double mui2 = 0.;
// g->gg and g->qqbar
if ( abs(info.emitterData()->id()) == abs(info.emissionData()->id()) ) {
mui2 = sqr(info.emitterData()->mass()/info.scale());
}
// Otherwise have X->Xg (should work for SUSY)
else {
mui2 = Mui2;
}
double mu2 = sqr(info.emissionData()->mass()/info.scale() );
// To solve issue with scale during presampling
// need to enforce that Qijk-mij2-mk2 = 2*pij.pk > 0,
// so combine checks by comparing against square root.
if ( 1.-Mui2-Muj2 < sqrt(4.*Mui2*Muj2) ) {
jacobian(0.0);
return false;
}
Energy2 Qijk = sqr(info.scale());
double suijk = 0.5*( 1. - Mui2 - Muj2 + sqrt( sqr(1.-Mui2-Muj2) - 4.*Mui2*Muj2 ) );
// Calculate pt
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
Energy2 pt2 = sqr(pt);
if ( pt > info.hardPt() || pt < IRCutoff() ) {
jacobian(0.0);
return false;
}
// Generate zPrime (i.e. the new definition of z specific to massive FF and decays)
double zPrime;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emissionData()->id() != ParticleID::g ) {
zPrime = xi;
}
else {
zPrime = exp(xi)/(1.+exp(xi));
}
}
else {
zPrime = 1.-exp(-xi);
}
// new: 2011-08-31
// 2011-11-08: this does happen
if( sqrt(mui2)+sqrt(mu2)+sqrt(muj2) > 1. ){
jacobian(0.0);
return false;
}
// Have derived and checked the equations for zp1 and zm1, these apply to zPrime
// phasespace constraint to incorporate ptMax
- Energy hard=info.hardPt();
+ Energy hard = info.hardPt();
+
+ if(openZBoundaries()>0){
+ // From ptMax(..)
+ hard = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
+ ( 2.-2.*sqrt(muj2) ) * info.scale();
+ assert(pt<=hard);
+ }
+
double ptRatio = sqrt(1.-sqr(pt/hard));
double zp1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) +
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) * ptRatio) /
( 2.*sqr(1.-sqrt(muj2)) );
double zm1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) -
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) * ptRatio) /
( 2.*sqr(1.-sqrt(muj2)) );
if ( zPrime > zp1 || zPrime < zm1 ) {
jacobian(0.0);
return false;
}
// Calculate A:=xij*w
double A = (1./(suijk*zPrime*(1.-zPrime))) * ( pt2/Qijk + zPrime*mu2 + (1.-zPrime)*mui2 - zPrime*(1.-zPrime)*Mui2 );
// Calculate y from A (can also write explicitly in terms of qt, zPrime and masses however we need A anyway)
double bar = 1.-mui2-mu2-muj2;
double y = (1./bar) * (A*suijk + Mui2 - mui2 - mu2 );
// kinematic phasespace boundaries for y
// same as in Dittmaier hep-ph/9904440v2 (equivalent to CS)
double ym = 2.*sqrt(mui2)*sqrt(mu2)/bar;
double yp = 1. - 2.*sqrt(muj2)*(1.-sqrt(muj2))/bar;
if ( y < ym || y > yp ) {
jacobian(0.0);
return false;
}
// Calculate xk and xij
double lambdaK = 1. + (Muj2/suijk);
double lambdaIJ = 1. + (Mui2/suijk);
double xk = (1./(2.*lambdaK)) * ( (lambdaK + (Muj2/suijk)*lambdaIJ - A) + sqrt( sqr(lambdaK + (Muj2/suijk)*lambdaIJ - A) - 4.*lambdaK*lambdaIJ*Muj2/suijk) );
double xij = 1. - ( (Muj2/suijk) * (1.-xk) / xk );
// Transform to standard z definition as used in the kernels (i.e. that used in CS and standard sudakov parametrisations)
double z = ( (zPrime*xij*xk*suijk/2.) + (Muj2/ ( 2.*xk*xij*suijk*zPrime))*(pt2/Qijk + mui2) ) /
( (xij*xk*suijk/2.) + (Muj2/(2.*xk*xij))*(Mui2/suijk + A) );
// These apply to z, not zPrime
double zm = ( (2.*mui2+bar*y)*(1.-y) - sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
double zp = ( (2.*mui2+bar*y)*(1.-y) + sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
if ( z < zm || z > zp ) {
jacobian(0.0);
return false;
}
double phi = 2.*Constants::pi*rphi;
double mapZJacobian;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emissionData()->id() != ParticleID::g ) {
mapZJacobian = 1.;
}
else {
mapZJacobian = z*(1.-z);
}
}
else {
mapZJacobian = 1.-z;
}
// Compute and store the jacobian
double jac = 0.0;
double propCntrb = 1./ ( 1. + (mui2+mu2-Mui2)/(bar*y) );
// TODO - SW - Tidy up notation everywhere alongside writing for the manual
// The full jacobian including the z->zprime jacobian
if ( theFullJacobian ) {
// Sort out variables
Energy2 sbar = Qijk*bar;
Energy2 sijk = Qijk*suijk;
Energy2 mi2 = Qijk*mui2;
Energy2 mj2 = Qijk*mu2;
Energy2 mk2 = Qijk*muj2;
// Compute dy/dzPrime and pt2* dy/dpt2
double dyBydzPrime = (1./sbar) * ( -pt2*(1.-2.*zPrime)/sqr(zPrime*(1.-zPrime)) - mi2/sqr(zPrime) + mj2/sqr(1.-zPrime) );
InvEnergy2 dyBydpt2 = 1./(sbar*zPrime*(1.-zPrime));
// Compute dA/dzPrime and dA/dpt2
double dABydzPrime = (sbar/sijk) * dyBydzPrime;
InvEnergy2 dABydpt2 = (sbar/sijk) * dyBydpt2;
// Compute dxk/dzPrime, dxk/dpt2, dxij/dzPrime and dxij/dpt2
double factor = (0.5/lambdaK) * (-1. - (1./sqrt( sqr(lambdaK + (mk2/sijk)*lambdaIJ - A) - 4.*lambdaK*lambdaIJ*mk2/sijk)) * (lambdaK + (mk2/sijk)*lambdaIJ - A));
double dxkBydzPrime = factor * dABydzPrime;
InvEnergy2 dxkBydpt2 = factor * dABydpt2;
double dxijBydzPrime = (mk2/sijk) * (1./sqr(xk)) * dxkBydzPrime;
InvEnergy2 dxijBydpt2 = (mk2/sijk) * (1./sqr(xk)) * dxkBydpt2;
Energy2 dqiDotqkBydzPrime = xij*xk*0.5*sijk + zPrime*dxijBydzPrime*xk*0.5*sijk + zPrime*xij*dxkBydzPrime*0.5*sijk
+ 0.5*(mk2/sijk)*(pt2 + mi2) * (-1./(xk*xij*sqr(zPrime)) - dxkBydzPrime/(zPrime*xij*sqr(xk)) - dxijBydzPrime/(zPrime*xk*sqr(xij)));
double dqiDotqkBydpt2 = dxijBydpt2*zPrime*xk*0.5*sijk + zPrime*xij*dxkBydpt2*0.5*sijk
+ (0.5*mk2/sijk) * (1./(zPrime*xk*xij)) * (1. + (pt2+mi2)*(-dxkBydpt2/xk - dxijBydpt2/xij) );
// Compute dzBydzPrime and dzBydpt2
Energy2 qiDotqk = (zPrime*xij*xk*sijk*0.5) + (mk2/ ( 2.*xk*xij*sijk*zPrime))*(pt2 + mi2);
double dzBydzPrime = (1./sbar) * ( 2.*qiDotqk*dyBydzPrime/sqr(1.-y) + (1./(1.-y)) * 2.*dqiDotqkBydzPrime );
InvEnergy2 dzBydpt2 = (1./sbar) * ( 2.*qiDotqk*dyBydpt2/sqr(1.-y) + (1./(1.-y)) * 2.*dqiDotqkBydpt2 );
double pt2Jac = pt2*abs(dzBydpt2*dyBydzPrime - dzBydzPrime*dyBydpt2);
// Include the other terms and calculate the jacobian
jac = propCntrb * bar / rootOfKallen(1.,Mui2,Muj2) * (1.-y)/y * pt2Jac;
}
else {
double jacPt2 = 1. / ( 1. + sqr(1.-zPrime)*Qijk*mui2/pt2 + zPrime*zPrime*Qijk*mu2/pt2 );
jac = propCntrb * bar / rootOfKallen(1.,Mui2,Muj2) * (1.-y) * jacPt2;
}
jacobian(jac * mapZJacobian * 2. * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()) );
// Record the physical variables, as used by the CS kernel definitions
lastPt(pt);
lastZ(z);
lastPhi(phi);
// Record zPrime for use in kinematics generation
splittingParameters().clear();
splittingParameters().push_back(zPrime);
if ( theMCCheck ) {
theMCCheck->book(1.,1.,info.scale(),info.hardPt(),pt,z,jacobian());
}
return true;
}
void FIMassiveDecayKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
// The only value stored in dInfo.lastSplittingParameters() should be zPrime
assert(dInfo.lastSplittingParameters().size() == 1 );
double zPrime = dInfo.lastSplittingParameters()[0];
Energy pt = dInfo.lastPt();
Energy2 pt2 = sqr(pt);
// Momentum of the recoil system
Lorentz5Momentum pk = pSpectator - pEmitter;
Lorentz5Momentum pij = pEmitter;
// scaled masses
double Mui2 = sqr(dInfo.emitterMass() / dInfo.scale());
double Muk2 = sqr(dInfo.recoilMass() / dInfo.scale());
//double muk2 = Muk2;
Energy mi = ZERO;
// g->gg and g->qqbar
if ( abs(dInfo.emitterData()->id()) == abs(dInfo.emissionData()->id()) ) {
mi = dInfo.emitterData()->mass();
}
// Otherwise have X->Xg (should work for SUSY)
else {
mi = dInfo.emitterMass();
}
double mui2 = sqr(mi/dInfo.scale());
double mu2 = sqr(dInfo.emissionData()->mass()/dInfo.scale() );
Energy2 Qijk = sqr(dInfo.scale());
double suijk = 0.5*( 1. - Mui2 - Muk2 + sqrt( sqr(1.-Mui2-Muk2) - 4.*Mui2*Muk2 ) );
double suijk2 = sqr(suijk);
// Calculate A:=xij*w
double A = (1./(suijk*zPrime*(1.-zPrime))) * ( pt2/Qijk + zPrime*mu2 + (1.-zPrime)*mui2 - zPrime*(1.-zPrime)*Mui2 );
// Calculate the scaling factors, xk and xij
double lambdaK = 1. + (Muk2/suijk);
double lambdaIJ = 1. + (Mui2/suijk);
double xk = (1./(2.*lambdaK)) * ( (lambdaK + (Muk2/suijk)*lambdaIJ - A) + sqrt( sqr(lambdaK + (Muk2/suijk)*lambdaIJ - A) - 4.*lambdaK*lambdaIJ*Muk2/suijk) );
double xij = 1. - ( (Muk2/suijk) * (1.-xk) / xk );
// Construct reference momenta nk, nij, nt
Lorentz5Momentum nij = ( suijk2 / (suijk2-Mui2*Muk2) ) * (pij - (Mui2/suijk)*pk);
Lorentz5Momentum nk = ( suijk2 / (suijk2-Mui2*Muk2) ) * (pk - (Muk2/suijk)*pij);
// Following notation in notes, qt = sqrt(wt)*nt
Lorentz5Momentum qt = getKt(nij,nk,pt,dInfo.lastPhi());
// Construct qij, qk, qi and qj
Lorentz5Momentum qij = xij*nij + (Mui2/(xij*suijk))*nk;
Lorentz5Momentum qk = xk*nk + (Muk2/(xk*suijk))*nij;
// No need to actually calculate nt and wt:
Lorentz5Momentum qi = zPrime*qij + ((pt2/Qijk + mui2 - zPrime*zPrime*Mui2)/(xij*suijk*zPrime))*nk + qt;
Lorentz5Momentum qj = (1.-zPrime)*qij + ((pt2/Qijk + mu2 - sqr(1.-zPrime)*Mui2)/(xij*suijk*(1.-zPrime)))*nk - qt;
qi.setMass(mi);
qi.rescaleEnergy();
qj.setMass(dInfo.emissionData()->mass());
qj.rescaleEnergy();
emitterMomentum(qi);
emissionMomentum(qj);
spectatorMomentum(pSpectator);
// Required for absorbing recoil in DipoleEventRecord::update
splitRecoilMomentum(qk);
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMassiveDecayKinematics::persistentOutput(PersistentOStream & os) const {
os << theFullJacobian;
}
void FIMassiveDecayKinematics::persistentInput(PersistentIStream & is, int) {
is >> theFullJacobian;
}
ClassDescription<FIMassiveDecayKinematics> FIMassiveDecayKinematics::initFIMassiveDecayKinematics;
// Definition of the static class description member.
void FIMassiveDecayKinematics::Init() {
static ClassDocumentation<FIMassiveDecayKinematics> documentation
("FIMassiveDecayKinematics implements implements massive splittings "
"off a final-initial decay dipole.");
static Switch<FIMassiveDecayKinematics,bool> interfaceFullJacobian
("FullJacobian",
"Use the full jacobian expression for the FI Decay kinematics.",
&FIMassiveDecayKinematics::theFullJacobian, true, false, false);
static SwitchOption interfaceFullJacobianYes
(interfaceFullJacobian,
"Yes",
"Use the full jacobian.",
true);
static SwitchOption interfaceFullJacobianNo
(interfaceFullJacobian,
"No",
"Do not use the full jacobian.",
false);
interfaceFullJacobian.rank(-1);
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:18 PM (1 d, 1 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242461
Default Alt Text
(46 KB)

Event Timeline