Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8309446
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
55 KB
Subscribers
None
View Options
diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am
--- a/Shower/QTilde/Makefile.am
+++ b/Shower/QTilde/Makefile.am
@@ -1,47 +1,48 @@
SUBDIRS = Matching
pkglib_LTLIBRARIES = HwShower.la
HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 26: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/HalfHalfOneEWSplitFn.h SplittingFunctions/HalfHalfOneEWSplitFn.cc\
SplittingFunctions/HalfHalfZeroEWSplitFn.h SplittingFunctions/HalfHalfZeroEWSplitFn.cc\
+SplittingFunctions/OneOneOneEWSplitFn.h SplittingFunctions/OneOneOneEWSplitFn.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\
Kinematics/Decay_QTildeShowerKinematics1to2.cc \
Kinematics/Decay_QTildeShowerKinematics1to2.h \
Kinematics/IS_QTildeShowerKinematics1to2.cc Kinematics/IS_QTildeShowerKinematics1to2.h \
Kinematics/FS_QTildeShowerKinematics1to2.cc Kinematics/FS_QTildeShowerKinematics1to2.h \
Kinematics/KinematicsReconstructor.cc \
Kinematics/KinematicsReconstructor.tcc \
Kinematics/KinematicsReconstructor.h \
Kinematics/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 \
Kinematics/ShowerKinematics.fh Kinematics/ShowerKinematics.h Kinematics/ShowerKinematics.cc \
Kinematics/ShowerBasis.fh Kinematics/ShowerBasis.h Kinematics/ShowerBasis.cc \
Base/ShowerProgenitor.fh Base/ShowerProgenitor.h \
SplittingFunctions/SudakovFormFactor.cc SplittingFunctions/SudakovFormFactor.h SplittingFunctions/SudakovFormFactor.fh \
SplittingFunctions/SudakovCutOff.cc SplittingFunctions/SudakovCutOff.h SplittingFunctions/SudakovCutOff.fh \
SplittingFunctions/PTCutOff.cc SplittingFunctions/PTCutOff.h \
SplittingFunctions/MassCutOff.cc SplittingFunctions/MassCutOff.h \
SplittingFunctions/VariableMassCutOff.cc SplittingFunctions/VariableMassCutOff.h \
SplittingFunctions/SplittingFunction.h SplittingFunctions/SplittingFunction.fh \
SplittingFunctions/SplittingFunction.cc \
Base/ShowerVertex.cc Base/ShowerVertex.fh Base/ShowerVertex.h
diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfOneEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/HalfHalfOneEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/HalfHalfOneEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/HalfHalfOneEWSplitFn.cc
@@ -1,216 +1,216 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HalfHalfOneEWSplitFn class.
//
#include "HalfHalfOneEWSplitFn.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/ParticleData.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
using namespace Herwig;
IBPtr HalfHalfOneEWSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr HalfHalfOneEWSplitFn::fullclone() const {
return new_ptr(*this);
}
void HalfHalfOneEWSplitFn::persistentOutput(PersistentOStream & os) const {
os << gZ_ << gWL_;
}
void HalfHalfOneEWSplitFn::persistentInput(PersistentIStream & is, int) {
is >> gZ_ >> gWL_;
}
-// The following static variable is needed for the type description system in ThePEG.
+// The following static variable is needed for the type description system in ThePEG.
DescribeClass<HalfHalfOneEWSplitFn,SplittingFunction>
describeHerwigHalfHalfOneEWSplitFn("Herwig::HalfHalfOneEWSplitFn", "HwShower.so");
void HalfHalfOneEWSplitFn::Init() {
static ClassDocumentation<HalfHalfOneEWSplitFn> documentation
("The HalfHalfOneEWSplitFn class implements the splitting q->qWand q->qZ");
}
void HalfHalfOneEWSplitFn::doinit() {
SplittingFunction::doinit();
tcSMPtr sm = generator()->standardModel();
double sw2 = sm->sin2ThetaW();
// left-handled W coupling
gWL_ = 1./sqrt(2.*sw2);
// Z couplings
double fact = 0.25/sqrt(sw2*(1.-sw2));
for(int ix=1;ix<4;++ix) {
gZ_[2*ix-1] = make_pair(fact*(sm->vd() + sm->ad()),
fact*(sm->vd() - sm->ad() ));
gZ_[2*ix ] = make_pair(fact*(sm->vu() + sm->au() ),
fact*(sm->vu() - sm->au() ));
gZ_[2*ix+9 ] = make_pair(fact*(sm->ve() + sm->ae() ),
fact*(sm->ve() - sm->ae() ));
gZ_[2*ix+10] = make_pair(fact*(sm->vnu() + sm->anu()),
fact*(sm->vnu() - sm->anu()));
}
}
void HalfHalfOneEWSplitFn::getCouplings(double & gL, double & gR, const IdList & ids) const {
if(ids[2]->id()==ParticleID::Z0) {
map<long,pair<double,double> >::const_iterator it = gZ_.find(abs(ids[0]->id()));
assert(it!=gZ_.end());
gL = it->second.first ;
gR = it->second.second;
}
else if(abs(ids[2]->id())==ParticleID::Wplus) {
gL = gWL_;
}
else
assert(false);
if(ids[0]->id()<0) swap(gL,gR);
}
double HalfHalfOneEWSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
double gL(0.),gR(0.);
getCouplings(gL,gR,ids);
double val = (1. + sqr(z))/(1.-z);
if(mass) {
- Energy m = ids[2]->mass();
+ Energy m = ids[2]->mass();
val -= sqr(m)/t;
}
- val *= (sqr(gL)*norm(rho(0,0))+sqr(gR)*norm(rho(1,1)));
+ val *= (sqr(gL)*abs(rho(0,0))+sqr(gR)*abs(rho(1,1)));
return colourFactor(ids)*val;
}
double HalfHalfOneEWSplitFn::overestimateP(const double z,
const IdList & ids) const {
double gL(0.),gR(0.);
getCouplings(gL,gR,ids);
- return 2.*max(sqr(gL),sqr(gR))*colourFactor(ids)/(1.-z);
+ return 2.*max(sqr(gL),sqr(gR))*colourFactor(ids)/(1.-z);
}
double HalfHalfOneEWSplitFn::ratioP(const double z, const Energy2 t,
const IdList & ids, const bool mass,
const RhoDMatrix & rho) const {
double gL(0.),gR(0.);
getCouplings(gL,gR,ids);
double val = 1. + sqr(z);
if(mass) {
- Energy m = ids[2]->mass();
+ Energy m = ids[2]->mass();
val -= (1.-z)*sqr(m)/t;
}
val *= (sqr(gL)*abs(rho(0,0))+sqr(gR)*abs(rho(1,1)))/max(sqr(gL),sqr(gR));
return 0.5*val;
-}
+}
double HalfHalfOneEWSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
double gL(0.),gR(0.);
getCouplings(gL,gR,ids);
double pre = colourFactor(ids)*max(sqr(gL),sqr(gR));
switch (PDFfactor) {
case 0:
return -2.*pre*Math::log1m(z);
case 1:
return 2.*pre*log(z/(1.-z));
case 2:
return 2.*pre/(1.-z);
case 3:
default:
throw Exception() << "HalfHalfOneEWSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
- }
+ }
}
double HalfHalfOneEWSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
double gL(0.),gR(0.);
getCouplings(gL,gR,ids);
double pre = colourFactor(ids)*max(sqr(gL),sqr(gR));
switch (PDFfactor) {
case 0:
- return 1. - exp(- 0.5*r/pre);
+ return 1. - exp(- 0.5*r/pre);
case 1:
return 1./(1.-exp(-0.5*r/pre));
case 2:
return 1.-2.*pre/r;
case 3:
default:
throw Exception() << "HalfHalfOneEWSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
- }
+ }
}
bool HalfHalfOneEWSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3) return false;
if(ids[2]->id()==ParticleID::Z0) {
- if(ids[0]->id()==ids[1]->id() &&
+ if(ids[0]->id()==ids[1]->id() &&
((ids[0]->id()>=1 && ids[0]->id()<=6) || (ids[0]->id()>=11&&ids[0]->id()<=16) )) return true;
}
else if(abs(ids[2]->id())==ParticleID::Wplus) {
if(!((ids[0]->id()>=1 && ids[0]->id()<=6) || (ids[0]->id()>=11&&ids[0]->id()<=16) )) return false;
if(!((ids[1]->id()>=1 && ids[1]->id()<=6) || (ids[1]->id()>=11&&ids[1]->id()<=16) )) return false;
if(ids[0]->id()+1!=ids[1]->id() && ids[0]->id()-1!=ids[1]->id()) return false;
int out = ids[1]->iCharge()+ids[2]->iCharge();
if(ids[0]->iCharge()==out) return true;
}
return false;
}
-vector<pair<int, Complex> >
+vector<pair<int, Complex> >
HalfHalfOneEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
-vector<pair<int, Complex> >
+vector<pair<int, Complex> >
HalfHalfOneEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
-DecayMEPtr HalfHalfOneEWSplitFn::matrixElement(const double z, const Energy2 t,
+DecayMEPtr HalfHalfOneEWSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool) {
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1)));
Energy m = ids[2]->mass();
double gL(0.),gR(0.);
getCouplings(gL,gR,ids);
double mt = m/sqrt(t);
double root = sqrt(1.-sqr(m)/t/(1-z));
double romz = sqrt(1.-z);
double rz = sqrt(z);
double r2 = sqrt(2.);
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
(*kernal)(0,0,0) = -phase*root*gL/romz;
(*kernal)(1,1,2) = cphase*root*gR/romz;
(*kernal)(0,0,2) = cphase*z*root*gL/romz;
(*kernal)(1,1,0) = -phase*z*root*gR/romz;
// long terms
(*kernal)(1,1,1) =-gR*mt*r2*rz/(1-z);
(*kernal)(0,0,1) =-gL*mt*r2*rz/(1-z);
// +- -+ terms zero due quark mass
for(unsigned int ix=0;ix<3;++ix) {
(*kernal)(1,0,ix) = 0.;
(*kernal)(0,1,ix) = 0.;
}
// return the answer
return kernal;
}
diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
@@ -1,248 +1,249 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HalfHalfZeroEWSplitFn class.
//
#include "HalfHalfZeroEWSplitFn.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/ParticleData.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "Herwig/Models/StandardModel/SMFFHVertex.h"
using namespace Herwig;
IBPtr HalfHalfZeroEWSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr HalfHalfZeroEWSplitFn::fullclone() const {
return new_ptr(*this);
}
void HalfHalfZeroEWSplitFn::persistentOutput(PersistentOStream & os) const {
os << ghqq_ << _theSM;
}
void HalfHalfZeroEWSplitFn::persistentInput(PersistentIStream & is, int) {
is >> ghqq_ >> _theSM;
}
// The following static variable is needed for the type description system in ThePEG.
DescribeClass<HalfHalfZeroEWSplitFn,SplittingFunction>
describeHerwigHalfHalfZeroEWSplitFn("Herwig::HalfHalfZeroEWSplitFn", "HwShower.so");
void HalfHalfZeroEWSplitFn::Init() {
static ClassDocumentation<HalfHalfZeroEWSplitFn> documentation
("The HalfHalfZeroEWSplitFn class implements the splitting q->qH");
}
void HalfHalfZeroEWSplitFn::doinit() {
SplittingFunction::doinit();
tcSMPtr sm = generator()->standardModel();
double sw2 = sm->sin2ThetaW();
ghqq_ = 1./sqrt(4.*sw2);
_theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
}
void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids) const {
if(abs(ids[2]->id())==ParticleID::h0) {
//get quark masses
Energy mq;
if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
else if(abs(ids[0]->id())==ParticleID::t)
mq = getParticleData(ParticleID::t)->mass();
Energy mW = getParticleData(ParticleID::Wplus)->mass();
gH = ghqq_*(mq/mW);
}
else
assert(false);
}
void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids, const Energy2 t) const {
if(abs(ids[2]->id())==ParticleID::h0) {
//get quark masses
Energy mq;
if(abs(ids[0]->id())==ParticleID::c)
mq = _theSM->mass(t,getParticleData(ParticleID::c));
else if(abs(ids[0]->id())==ParticleID::b)
mq = _theSM->mass(t,getParticleData(ParticleID::b));
else if(abs(ids[0]->id())==ParticleID::t)
mq = _theSM->mass(t,getParticleData(ParticleID::t));
Energy mW = getParticleData(ParticleID::Wplus)->mass();
//Energy mW = _theSM->mass(t,getParticleData(ParticleID::Wplus));
gH = ghqq_*(mq/mW);
}
else
assert(false);
}
double HalfHalfZeroEWSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
double gH(0.);
getCouplings(gH,ids,t);
double val = (1.-z);
Energy mq, mH;
//get masses
if(mass) {
mq = ids[0]->mass();
mH = ids[2]->mass();
}
else { // to assure the particle mass in non-zero
if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
else if(abs(ids[0]->id())==ParticleID::t)
mq = getParticleData(ParticleID::t)->mass();
mH = getParticleData(ParticleID::h0)->mass();
}
val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z);
val *= sqr(gH);
return colourFactor(ids)*val;
}
double HalfHalfZeroEWSplitFn::overestimateP(const double z,
const IdList & ids) const {
double gH(0.);
getCouplings(gH,ids);
return sqr(gH)*colourFactor(ids)*(1.-z);
}
double HalfHalfZeroEWSplitFn::ratioP(const double z, const Energy2 t,
const IdList & ids, const bool mass,
const RhoDMatrix & rho) const {
double gH(0.);
getCouplings(gH,ids,t);
double val = 1.;
Energy mq, mH;
if(mass) {
mq = ids[0]->mass();
mH = ids[2]->mass();
}
else { // to assure the particle mass in non-zero
if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
else if(abs(ids[0]->id())==ParticleID::t)
mq = getParticleData(ParticleID::t)->mass();
mH = getParticleData(ParticleID::h0)->mass();
}
val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z);
return val;
}
double HalfHalfZeroEWSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
double gH(0.);
getCouplings(gH,ids);
double pre = colourFactor(ids)*sqr(gH);
switch (PDFfactor) {
case 0: //OverP
return pre*(z-sqr(z)/2.);
case 1: //OverP/z
return pre*(log(z)-z);
case 2: //OverP/(1-z)
return pre*z;
case 3: //OverP/[z(1-z)]
return pre*log(z);
default:
throw Exception() << "HalfHalfZeroEWSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double HalfHalfZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
double gH(0.);
getCouplings(gH,ids);
double pre = colourFactor(ids)*sqr(gH);
switch (PDFfactor) {
case 0:
return 1. - sqrt(1. - 2.*r/pre);
case 1: //OverP/z
case 2: //OverP/(1-z)
return r/pre;
case 3: //OverP/[z(1-z)]
return exp(r/pre);
default:
throw Exception() << "HalfHalfZeroEWSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool HalfHalfZeroEWSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3) return false;
if(ids[2]->id()==ParticleID::h0) {
if(ids[0]->id()==ids[1]->id() && (ids[0]->id()==4 || ids[0]->id()==5 || ids[0]->id()==6))
return true;
}
return false;
}
vector<pair<int, Complex> >
HalfHalfZeroEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
vector<pair<int, Complex> >
HalfHalfZeroEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
DecayMEPtr HalfHalfZeroEWSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool) {
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0)));
//get masses
Energy mq, mH;
if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
else if(abs(ids[0]->id())==ParticleID::t)
mq = getParticleData(ParticleID::t)->mass();
mH = getParticleData(ParticleID::h0)->mass();
double gH(0.);
- getCouplings(gH,ids,t);
- double mqt = mq/sqrt(t);
- double mHt = mH/sqrt(t);
+ Energy2 tC = t/(z*(1-z));
+ getCouplings(gH,ids,tC);
+ double mqt = mq/sqrt(tC);
+ double mHt = mH/sqrt(tC);
double num1 = gH*(1.+z)*mqt;
double num2 = gH*sqrt(-sqr(mqt)*(1.-z) - sqr(mHt)*z + z*(1.-z)*(sqr(mqt)+z*(1.-z))); //watch this
double dnum = sqrt(2.)*sqrt((1.-z)*sqr(z));
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
(*kernal)(0,0,0) = num1/dnum;
(*kernal)(0,1,0) = cphase*num2/dnum;
(*kernal)(1,0,0) = -phase*num2/dnum;
(*kernal)(1,1,0) = num1/dnum;
cerr << "testing kernal";
for(unsigned int ix=0;ix<2;++ix) {
for(unsigned int iy=0;iy<2;++iy)
cerr << (*kernal)(ix,iy,0) << " ";
cerr << "\n";
}
// return the answer
return kernal;
}
diff --git a/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.cc
new file mode 100644
--- /dev/null
+++ b/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.cc
@@ -0,0 +1,302 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the OneOneOneEWSplitFn class.
+//
+
+#include "OneOneOneEWSplitFn.h"
+#include "ThePEG/StandardModel/StandardModelBase.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/PDT/ParticleData.h"
+#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
+#include "Herwig/Models/StandardModel/SMFFHVertex.h"
+
+using namespace Herwig;
+
+IBPtr OneOneOneEWSplitFn::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr OneOneOneEWSplitFn::fullclone() const {
+ return new_ptr(*this);
+}
+
+void OneOneOneEWSplitFn::persistentOutput(PersistentOStream & os) const {
+ os << gWWG_ << gWWZ_ << _theSM;
+}
+
+void OneOneOneEWSplitFn::persistentInput(PersistentIStream & is, int) {
+ is >> gWWG_ >> gWWZ_ >> _theSM;
+}
+
+// The following static variable is needed for the type description system in ThePEG.
+DescribeClass<OneOneOneEWSplitFn,SplittingFunction>
+describeHerwigOneOneOneEWSplitFn("Herwig::OneOneOneEWSplitFn", "HwShower.so");
+
+
+void OneOneOneEWSplitFn::Init() {
+
+ static ClassDocumentation<OneOneOneEWSplitFn> documentation
+ ("The OneOneOneEWSplitFn class implements the splitting W->WG, W->WZ and Z->ZZ");
+
+}
+
+
+void OneOneOneEWSplitFn::doinit() {
+ SplittingFunction::doinit();
+ tcSMPtr sm = generator()->standardModel();
+ double sw2 = sm->sin2ThetaW();
+ // WWZ coupling
+ gWWZ_ = 1.;
+ // WWG coupling
+ gWWG_ = sqrt((1.-sw2)/sw2);
+ // to employ running masses, wherever needed
+ _theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
+}
+
+
+void OneOneOneEWSplitFn::getCouplings(double & gvvv, const IdList & ids) const {
+ // G > WW
+ if(ids[0]->id()==ParticleID::gamma && abs(ids[1]->id())==ParticleID::Wplus
+ && abs(ids[2]->id())==ParticleID::Wplus){
+ gvvv = gWWG_;
+ }
+ // Z > WW
+ else if(ids[0]->id()==ParticleID::Z0 && abs(ids[1]->id())==ParticleID::Wplus
+ && abs(ids[2]->id())==ParticleID::Wplus){
+ gvvv = gWWZ_;
+ }
+ // W > WG
+ else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[2]->id()==ParticleID::gamma){
+ gvvv = gWWG_;
+ }
+ // W > WZ
+ else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[2]->id()==ParticleID::Z0){
+ gvvv = gWWZ_;
+ }
+ else
+ assert(false);
+}
+
+
+double OneOneOneEWSplitFn::P(const double z, const Energy2 t,
+ const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
+ double gvvv(0.);
+ getCouplings(gvvv,ids);
+ double val(0.);
+ val = ((2.*sqr(1.-(1.-z)*z))/((1.-z)*z))*abs(rho(-1,-1)+rho(1,1));
+ if(mass) {
+ double mWt2 = sqr(getParticleData(ParticleID::Wplus)->mass())/t;
+ double mZt2 = sqr(getParticleData(ParticleID::Z0)->mass())/t;
+ // G > WW
+ if(ids[0]->id()==ParticleID::gamma && abs(ids[1]->id())==ParticleID::Wplus
+ && abs(ids[2]->id())==ParticleID::Wplus){
+ val += (-2.*mWt2*(2.-(1.-z)*z)*abs(rho(-1,-1)+rho(1,1)))/((1.-z)*z);
+ }
+ // Z > WW
+ else if(ids[0]->id()==ParticleID::Z0 && abs(ids[1]->id())==ParticleID::Wplus
+ && abs(ids[2]->id())==ParticleID::Wplus){
+ val += (-2.*mWt2*(2.-(1.-z)*z)*abs(rho(-1,-1)+rho(1,1)))/((1.-z)*z);
+ val += (-2.*mZt2*(z*(-2.+z*(8.+z*(-13.-2.*(-4.+z)*z)))*abs(rho(0,0))
+ - (1.-z)*sqr(1.-(1.-z)*z)*abs(rho(-1,-1)+rho(1,1))))/(sqr(1.-z)*z);
+ }
+ // W > WG
+ else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[2]->id()==ParticleID::gamma){
+ val += 4.*mWt2*sqr(1.-z)*abs(rho(0,0));
+ val -= 2.*mWt2*(1.+sqr(1.-z))*abs(rho(-1,-1) + rho(1,1));
+ }
+ // W > WZ
+ else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[2]->id()==ParticleID::Z0){
+ val += (2.*mZt2*(pow(z,3)*abs(rho(0,0)) + (1.-z)*(-1.+sqr(1.-z)*z)
+ * abs(rho(-1,-1)+rho(1,1))))/(sqr(1.-z)*z);
+ val += -2.*mWt2*(-2.*sqr(1.-z)*abs(rho(0,0))+(2.+(-2.+z)*z)
+ * abs(rho(-1,-1)+rho(1,1)));
+ }
+ }
+ return sqr(gvvv)*val;
+}
+
+
+double OneOneOneEWSplitFn::overestimateP(const double z,
+ const IdList & ids) const {
+ double gvvv(0.);
+ getCouplings(gvvv,ids);
+ return sqr(gvvv)*(2.*sqr(gvvv)/(z*(1.-z)));
+}
+
+
+double OneOneOneEWSplitFn::ratioP(const double z, const Energy2 t,
+ const IdList & ids, const bool mass,
+ const RhoDMatrix & rho) const {
+ double gvvv(0.);
+ getCouplings(gvvv,ids);
+ double val = sqr(1.-(1.-z)*z)*abs(rho(-1,-1)+rho(1,1));
+ if(mass) {
+ double mWt2 = sqr(getParticleData(ParticleID::Wplus)->mass())/t;
+ double mZt2 = sqr(getParticleData(ParticleID::Z0)->mass())/t;
+ // G > WW
+ if(ids[0]->id()==ParticleID::gamma && abs(ids[1]->id())==ParticleID::Wplus
+ && abs(ids[2]->id())==ParticleID::Wplus) {
+ val += (-2.*mWt2*(2.-(1.-z)*z)*abs(rho(-1,-1)+rho(1,1)))/((1.-z)*z);
+ }
+ // Z > WW
+ else if(ids[0]->id()==ParticleID::Z0 && abs(ids[1]->id())==ParticleID::Wplus
+ && abs(ids[2]->id())==ParticleID::Wplus){
+ val += (-2.*mWt2*(2.-(1.-z)*z)*abs(rho(-1,-1)+rho(1,1)))/((1.-z)*z);
+ val += (-2.*mZt2*(z*(-2.+z*(8.+z*(-13.-2.*(-4.+z)*z)))*abs(rho(0,0))
+ - (1.-z)*sqr(1.-(1.-z)*z)*abs(rho(-1,-1)+rho(1,1))))/(sqr(1.-z)*z);
+ }
+ // W > WG
+ else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[2]->id()==ParticleID::gamma) {
+ val += 4.*mWt2*sqr(1.-z)*abs(rho(0,0));
+ val -= 2.*mWt2*(1.+sqr(1.-z))*abs(rho(-1,-1) + rho(1,1));
+ }
+ // W > WZ
+ else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[2]->id()==ParticleID::Z0) {
+ val += (2.*mZt2*(pow(z,3)*abs(rho(0,0)) + (1.-z)*(-1.+sqr(1.-z)*z)
+ * abs(rho(-1,-1)+rho(1,1))))/(sqr(1.-z)*z);
+ val += -2.*mWt2*(-2.*sqr(1.-z)*abs(rho(0,0))+(2.+(-2.+z)*z)
+ * abs(rho(-1,-1)+rho(1,1)));
+ }
+ }
+ return sqr(gvvv)*val;
+}
+
+
+double OneOneOneEWSplitFn::integOverP(const double z,
+ const IdList & ids,
+ unsigned int PDFfactor) const {
+ double gvvv(0.);
+ getCouplings(gvvv,ids);
+ double pre = sqr(gvvv);
+ switch (PDFfactor) {
+ case 0:
+ return 2.*pre*(log(z)-log(1.-z));
+ case 1:
+ //return -2.*pre*(1./z+log(1.-z)-log(z));
+ case 2:
+ //return 2.*pre*(2.*log(z)+(2.*z-1.)/(z*(1.-z))-2.*log(1.-z));
+ case 3:
+ //return 2.*pre*(1./(1.-z)-1./z-2.*log(1.-z)+2.*log(z));
+ default:
+ throw Exception() << "OneOneOneEWSplitFn::integOverP() invalid PDFfactor = "
+ << PDFfactor << Exception::runerror;
+ }
+}
+
+double OneOneOneEWSplitFn::invIntegOverP(const double r, const IdList & ids,
+ unsigned int PDFfactor) const {
+ double gvvv(0.);
+ getCouplings(gvvv,ids);
+ double pre = sqr(gvvv);
+ switch (PDFfactor) {
+ case 0:
+ return exp(0.5*r/pre)/(1.+exp(0.5*r/pre));
+ case 1:
+ case 2:
+ case 3:
+ default:
+ throw Exception() << "OneOneOneEWSplitFn::invIntegOverP() invalid PDFfactor = "
+ << PDFfactor << Exception::runerror;
+ }
+}
+
+bool OneOneOneEWSplitFn::accept(const IdList &ids) const {
+ if(ids.size()!=3) return false;
+ if(ids[0]->id()==ParticleID::gamma && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[1]->id()==-ids[2]->id())
+ return true;
+
+ if(ids[0]->id()==ParticleID::Z0 && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[1]->id()==-ids[2]->id())
+ return true;
+
+ if(abs(ids[0]->id())==ParticleID::Wplus) {
+ if(ids[1]->id()==ids[0]->id() && ids[2]->id()==ParticleID::gamma)
+ return true;
+ if(ids[1]->id()==ids[0]->id() && ids[2]->id()==ParticleID::Z0)
+ return true;
+ }
+ return false;
+}
+
+
+vector<pair<int, Complex> >
+OneOneOneEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
+ const RhoDMatrix &) {
+ // no dependence on the spin density matrix, dependence on off-diagonal terms cancels
+ // and rest = splitting function for Tr(rho)=1 as required by defn
+ return vector<pair<int, Complex> >(1,make_pair(0,1.));
+}
+
+
+vector<pair<int, Complex> >
+OneOneOneEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
+ const RhoDMatrix &) {
+ // no dependence on the spin density matrix, dependence on off-diagonal terms cancels
+ // and rest = splitting function for Tr(rho)=1 as required by defn
+ return vector<pair<int, Complex> >(1,make_pair(0,1.));
+}
+
+
+DecayMEPtr OneOneOneEWSplitFn::matrixElement(const double z, const Energy2 t,
+ const IdList & ids, const double phi,
+ bool) {
+ // calculate the kernal
+ DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1)));
+ double gvvv(0.);
+ getCouplings(gvvv,ids);
+ // defining dummies
+ double m0t = ids[0]->mass()/sqrt(t);
+ double m1t = ids[1]->mass()/sqrt(t);
+ double m2t = ids[2]->mass()/sqrt(t);
+ Complex phase = exp(Complex(0.,1.)*phi);
+ Complex cphase = conj(phase);
+
+ double z1_z = z*(1.-z);
+ double sqrtmass = sqrt(sqr(m0t)-sqr(m1t)/z-sqr(m2t)/(1.-z)+1.);
+ double r2 = sqrt(2.);
+ // assign kernel
+ (*kernal)(0,0,0) = gvvv*(phase/sqrt(z1_z))*sqrtmass;
+ (*kernal)(0,0,1) = gvvv*r2*m2t*(z/(1.-z)); //2>4
+ (*kernal)(0,0,2) = -gvvv*cphase*sqrt(z/(1.-z))*sqrtmass;
+ (*kernal)(0,1,0) = -gvvv*r2*m1t*(1.-z)/z; //2>4
+ (*kernal)(0,1,1) = 0.;
+ (*kernal)(0,1,2) = 0.;
+ (*kernal)(0,2,0) = -gvvv*(1.-z)*cphase*sqrt((1.-z)/z)*sqrtmass;
+ (*kernal)(0,2,1) = 0.;
+ (*kernal)(0,2,2) = 0.;
+
+ (*kernal)(1,0,0) = 0.;
+ (*kernal)(1,0,1) = 0.; //2>4
+ (*kernal)(1,0,2) = -gvvv*r2*m0t*(1.-z); //2>4
+ (*kernal)(1,1,0) = gvvv*phase*(m0t/m1t)*sqrt(z/(1.-z)); //221>421
+ (*kernal)(1,1,1) = gvvv*r2*(m0t*m2t/m1t)*(z/(1.-z)); //222>424
+ (*kernal)(1,1,2) = -gvvv*cphase*(m0t/m1t)*sqrt(z/(1.-z)); //223>423
+ (*kernal)(1,2,0) = -gvvv*r2*m0t*(1.-z); //2>4
+ (*kernal)(1,2,1) = 0.; //2>4
+ (*kernal)(1,2,2) = 0.; //2>4
+
+ (*kernal)(2,0,0) = 0.;
+ (*kernal)(2,0,1) = 0.;
+ (*kernal)(2,0,2) = gvvv*(1.-z)*phase*sqrt((1.-z)/z)*sqrtmass;
+ (*kernal)(2,1,0) = 0.;
+ (*kernal)(2,1,1) = 0.; //2>4
+ (*kernal)(2,1,2) = -gvvv*r2*m1t*((1.-z)/z);//2>4
+ (*kernal)(2,2,0) = gvvv*phase*sqrt(z/(1.-z))*sqrtmass;
+ (*kernal)(2,2,1) = gvvv*r2*m2t*(z/(1.-z)); //2>4
+ (*kernal)(2,2,2) = -gvvv*(cphase/sqrt(z1_z))*sqrtmass;
+
+ // return the answer
+ return kernal;
+}
diff --git a/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.h b/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.h
new file mode 100644
--- /dev/null
+++ b/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.h
@@ -0,0 +1,217 @@
+// -*- C++ -*-
+#ifndef Herwig_OneOneOneEWSplitFn_H
+#define Herwig_OneOneOneEWSplitFn_H
+//
+// This is the declaration of the OneOneOneEWSplitFn class.
+//
+
+#include "SplittingFunction.h"
+#include "Herwig/Models/StandardModel/StandardModel.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The OneOneOneEWSplitFn class implements the splitting function for
+ * \f$\1\to q\1 1\f$ where the spin-1 particles are massive electroweak gauge bosons.
+ *
+ * @see \ref OneOneOneEWSplitFnInterfaces "The interfaces"
+ * defined for OneOneOneEWSplitFn.
+ */
+class OneOneOneEWSplitFn: public SplittingFunction {
+
+public:
+
+ /**
+ * Concrete implementation of the method to determine whether this splitting
+ * function can be used for a given set of particles.
+ * @param ids The PDG codes for the particles in the splitting.
+ */
+ virtual bool accept(const IdList & ids) const;
+
+ /**
+ * Methods to return the splitting function.
+ */
+ //@{
+ /**
+ * The concrete implementation of the splitting function, \f$P(z,t)\f$.
+ * @param z The energy fraction.
+ * @param t The scale.
+ * @param ids The PDG codes for the particles in the splitting.
+ * @param mass Whether or not to include the mass dependent terms
+ * @param rho The spin density matrix
+ */
+ virtual double P(const double z, const Energy2 t, const IdList & ids,
+ const bool mass, const RhoDMatrix & rho) const;
+
+ /**
+ * The concrete implementation of the overestimate of the splitting function,
+ * \f$P_{\rm over}\f$.
+ * @param z The energy fraction.
+ * @param ids The PDG codes for the particles in the splitting.
+ */
+ virtual double overestimateP(const double z, const IdList & ids) const;
+
+ /**
+ * The concrete implementation of the
+ * the ratio of the splitting function to the overestimate, i.e.
+ * \f$P(z,t)/P_{\rm over}(z)\f$.
+ * @param z The energy fraction.
+ * @param t The scale.
+ * @param ids The PDG codes for the particles in the splitting.
+ * @param mass Whether or not to include the mass dependent terms
+ * @param rho The spin density matrix
+ */
+ virtual double ratioP(const double z, const Energy2 t, const IdList & ids,
+ const bool mass, const RhoDMatrix & rho) const;
+
+ /**
+ * The concrete implementation of the indefinite integral of the
+ * overestimated splitting function, \f$P_{\rm over}\f$.
+ * @param z The energy fraction.
+ * @param ids The PDG codes for the particles in the splitting.
+ * @param PDFfactor Which additional factor to include for the PDF
+ * 0 is no additional factor,
+ * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
+ */
+ virtual double integOverP(const double z, const IdList & ids,
+ unsigned int PDFfactor=0) const;
+
+ /**
+ * The concrete implementation of the inverse of the indefinite integral.
+ * @param r Value of the splitting function to be inverted
+ * @param ids The PDG codes for the particles in the splitting.
+ * @param PDFfactor Which additional factor to include for the PDF
+ * 0 is no additional factor,
+ * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
+ */
+ virtual double invIntegOverP(const double r, const IdList & ids,
+ unsigned int PDFfactor=0) const;
+ //@}
+
+ /**
+ * Method to calculate the azimuthal angle
+ * @param z The energy fraction
+ * @param t The scale \f$t=2p_j\cdot p_k\f$.
+ * @param ids The PDG codes for the particles in the splitting.
+ * @param The azimuthal angle, \f$\phi\f$.
+ * @return The weight
+ */
+ virtual vector<pair<int,Complex> >
+ generatePhiForward(const double z, const Energy2 t, const IdList & ids,
+ const RhoDMatrix &);
+
+ /**
+ * Method to calculate the azimuthal angle for backward evolution
+ * @param z The energy fraction
+ * @param t The scale \f$t=2p_j\cdot p_k\f$.
+ * @param ids The PDG codes for the particles in the splitting.
+ * @param The azimuthal angle, \f$\phi\f$.
+ * @return The weight
+ */
+ virtual vector<pair<int,Complex> >
+ generatePhiBackward(const double z, const Energy2 t, const IdList & ids,
+ const RhoDMatrix &);
+
+ /**
+ * Calculate the matrix element for the splitting
+ * @param z The energy fraction
+ * @param t The scale \f$t=2p_j\cdot p_k\f$.
+ * @param ids The PDG codes for the particles in the splitting.
+ * @param The azimuthal angle, \f$\phi\f$.
+ */
+ virtual DecayMEPtr matrixElement(const double z, const Energy2 t,
+ const IdList & ids, const double phi, bool timeLike);
+
+protected:
+
+ /**
+ * Get the couplings
+ */
+ void getCouplings(double & gvvv, const IdList & ids) 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;
+ //@}
+
+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();
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ OneOneOneEWSplitFn & operator=(const OneOneOneEWSplitFn &) = delete;
+
+private:
+
+ /**
+ * W^{\pm} -> W^{\pm} G couplings
+ */
+ double gWWG_;
+
+ /**
+ * W^{\pm} -> W^{\pm} Z couplings
+ */
+ double gWWZ_;
+
+ /**
+ * Pointer to the SM object.
+ */
+ tcHwSMPtr _theSM;
+};
+
+}
+
+#endif /* Herwig_OneOneOneEWSplitFn_H */
diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in
--- a/src/defaults/Shower.in
+++ b/src/defaults/Shower.in
@@ -1,380 +1,398 @@
# -*- ThePEG-repository -*-
############################################################
# Setup of default parton shower
#
-# Useful switches for users are marked near the top of
+# 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
set AlphaQED:CouplingSource Thompson
create Herwig::ShowerAlphaQED AlphaEW
set AlphaEW:CouplingSource MZ
create Herwig::PartnerFinder PartnerFinder
newdef PartnerFinder:PartnerMethod 1
newdef PartnerFinder:ScaleChoice 1
create Herwig::KinematicsReconstructor KinematicsReconstructor
newdef KinematicsReconstructor:ReconstructionOption Colour3
newdef KinematicsReconstructor:InitialStateReconOption SofterFraction
newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost
newdef KinematicsReconstructor:FinalFinalWeight Yes
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:Interactions QEDQCD
newdef ShowerHandler:SpinCorrelations Yes
newdef ShowerHandler:SoftCorrelations Singular
##################################################################
# Intrinsic pT
#
-# Recommended:
+# 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 QEDQCD
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:EvolutionScheme DotProduct
#############################################################
# End of interesting user servicable section.
#
-# Anything that follows below should only be touched if you
-# know what you're doing.
+# 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:EvolutionScheme DotProduct
newdef AlphaQCD:ScaleFactor 1.0
newdef AlphaQCD:NPAlphaS 2
newdef AlphaQCD:Qmin 0.935
newdef AlphaQCD:NumberOfLoops 2
newdef AlphaQCD:AlphaIn 0.1186
#
#
# Lets set up all the splittings
create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn
set QtoQGammaSplitFn:InteractionType QED
set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral
set QtoQGammaSplitFn:AngularOrdered Yes
set QtoQGammaSplitFn:StrictAO Yes
create Herwig::HalfHalfOneSplitFn QtoQGSplitFn
newdef QtoQGSplitFn:InteractionType QCD
newdef QtoQGSplitFn:ColourStructure TripletTripletOctet
set QtoQGSplitFn:AngularOrdered Yes
set QtoQGSplitFn:StrictAO Yes
create Herwig::OneOneOneSplitFn GtoGGSplitFn
newdef GtoGGSplitFn:InteractionType QCD
newdef GtoGGSplitFn:ColourStructure OctetOctetOctet
set GtoGGSplitFn:AngularOrdered Yes
set GtoGGSplitFn:StrictAO Yes
create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn
newdef WtoWGammaSplitFn:InteractionType QED
newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral
set WtoWGammaSplitFn:AngularOrdered Yes
set WtoWGammaSplitFn:StrictAO Yes
create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn
newdef GtoQQbarSplitFn:InteractionType QCD
newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet
set GtoQQbarSplitFn:AngularOrdered Yes
set GtoQQbarSplitFn:StrictAO Yes
create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn
newdef GammatoQQbarSplitFn:InteractionType QED
newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged
set GammatoQQbarSplitFn:AngularOrdered Yes
set GammatoQQbarSplitFn:StrictAO Yes
create Herwig::HalfOneHalfSplitFn QtoGQSplitFn
newdef QtoGQSplitFn:InteractionType QCD
newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet
set QtoGQSplitFn:AngularOrdered Yes
set QtoGQSplitFn:StrictAO Yes
create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn
newdef QtoGammaQSplitFn:InteractionType QED
newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged
set QtoGammaQSplitFn:AngularOrdered Yes
set QtoGammaQSplitFn:StrictAO Yes
create Herwig::HalfHalfOneEWSplitFn QtoQWZSplitFn
newdef QtoQWZSplitFn:InteractionType EW
newdef QtoQWZSplitFn:ColourStructure EW
create Herwig::HalfHalfZeroEWSplitFn QtoQHSplitFn
newdef QtoQHSplitFn:InteractionType EW
newdef QtoQHSplitFn:ColourStructure EW
+create Herwig::OneOneOneEWSplitFn VtoVVSplitFn
+newdef VtoVVSplitFn:InteractionType EW
+newdef VtoVVSplitFn:ColourStructure EW
+
#
# Now the Sudakovs
create Herwig::PTCutOff PTCutOff
newdef PTCutOff:pTmin 0.958*GeV
create Herwig::SudakovFormFactor SudakovCommon
newdef SudakovCommon:Alpha AlphaQCD
newdef SudakovCommon:Cutoff PTCutOff
newdef SudakovCommon:PDFmax 1.0
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
cp PTCutOff LtoLGammaPTCutOff
# Technical parameter to stop evolution.
set LtoLGammaPTCutOff:pTmin 0.000001
set LtoLGammaSudakov:Cutoff LtoLGammaPTCutOff
cp SudakovCommon QtoQWZSudakov
set QtoQWZSudakov:SplittingFunction QtoQWZSplitFn
set QtoQWZSudakov:Alpha AlphaEW
set QtoQWZSudakov:PDFmax 1.9
cp QtoQWZSudakov LtoLWZSudakov
cp SudakovCommon QtoQHSudakov
set QtoQHSudakov:SplittingFunction QtoQHSplitFn
set QtoQHSudakov:Alpha AlphaEW
set QtoQHSudakov:PDFmax 1.9
cp QtoQHSudakov LtoLHSudakov
+cp SudakovCommon VtoVVSudakov
+set VtoVVSudakov:SplittingFunction VtoVVSplitFn
+set VtoVVSudakov:Alpha AlphaEW
+set VtoVVSudakov:PDFmax 1.9
+
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
#
# Electroweak
#
do SplittingGenerator:AddFinalSplitting u->u,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting d->d,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting s->s,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting c->c,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting b->b,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting t->t,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting u->u,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting d->d,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting s->s,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting c->c,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting b->b,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting t->t,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting u->d,W+; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting c->s,W+; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting d->u,W-; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting s->c,W-; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting u->d,W+; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting c->s,W+; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting d->u,W-; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting s->c,W-; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting c->c,h0; QtoQHSudakov
do SplittingGenerator:AddFinalSplitting b->b,h0; QtoQHSudakov
do SplittingGenerator:AddFinalSplitting t->t,h0; QtoQHSudakov
do SplittingGenerator:AddInitialSplitting c->c,h0; QtoQHSudakov
do SplittingGenerator:AddInitialSplitting b->b,h0; QtoQHSudakov
do SplittingGenerator:AddInitialSplitting t->t,h0; QtoQHSudakov
+
+do SplittingGenerator:AddFinalSplitting gamma->W+,W-; VtoVVSudakov
+do SplittingGenerator:AddFinalSplitting gamma->W-,W+; VtoVVSudakov
+do SplittingGenerator:AddFinalSplitting Z0->W+,W-; VtoVVSudakov
+do SplittingGenerator:AddFinalSplitting Z0->W-,W+; VtoVVSudakov
+do SplittingGenerator:AddFinalSplitting W+->W+,gamma; VtoVVSudakov
+do SplittingGenerator:AddFinalSplitting W-->W-,gamma; VtoVVSudakov
+do SplittingGenerator:AddFinalSplitting W+->W+,Z0; VtoVVSudakov
+do SplittingGenerator:AddFinalSplitting W-->W-,Z0; VtoVVSudakov
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 3:24 PM (1 d, 6 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023276
Default Alt Text
(55 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment