Page MenuHomeHEPForge

No OneTemporary

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

Mime Type
text/x-diff
Expires
Sat, Dec 21, 3:24 PM (1 d, 40 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023276
Default Alt Text
(55 KB)

Event Timeline