Page MenuHomeHEPForge
Files F8309674

No OneTemporary
Actions

View Options

diff --git a/Decay/Perturbative/SMWDecayer.cc b/Decay/Perturbative/SMWDecayer.cc
--- a/Decay/Perturbative/SMWDecayer.cc
+++ b/Decay/Perturbative/SMWDecayer.cc
@@ -1,929 +1,924 @@
// -*- C++ -*-
//
// SMWDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SMWDecayer class.
//
#include "SMWDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Decay/DecayVertex.h"
#include "ThePEG/Helicity/VectorSpinInfo.h"
#include "ThePEG/Helicity/FermionSpinInfo.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "Herwig/Shower/Core/Base/ShowerProgenitor.h"
#include "Herwig/Shower/Core/Base/ShowerParticle.h"
#include "Herwig/Shower/Core/Base/Branching.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include <numeric>
using namespace Herwig;
using namespace ThePEG::Helicity;
const double SMWDecayer::EPS_=0.00000001;
SMWDecayer::SMWDecayer()
- : quarkWeight_(6,0.), leptonWeight_(3,0.), CF_(4./3.), pTmin_(1.*GeV), NLO_(false) {
+ : quarkWeight_(6,0.), leptonWeight_(3,0.), CF_(4./3.),
+ NLO_(false) {
quarkWeight_[0] = 1.01596;
quarkWeight_[1] = 0.0537308;
quarkWeight_[2] = 0.0538085;
quarkWeight_[3] = 1.01377;
quarkWeight_[4] = 1.45763e-05;
quarkWeight_[5] = 0.0018143;
leptonWeight_[0] = 0.356594;
leptonWeight_[1] = 0.356593;
leptonWeight_[2] = 0.356333;
// intermediates
generateIntermediates(false);
}
void SMWDecayer::doinit() {
PerturbativeDecayer::doinit();
// get the vertices from the Standard Model object
tcHwSMPtr hwsm=dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException() << "Must have Herwig StandardModel object in"
<< "SMWDecayer::doinit()"
<< Exception::runerror;
FFWVertex_ = hwsm->vertexFFW();
FFGVertex_ = hwsm->vertexFFG();
// make sure they are initialized
FFGVertex_->init();
FFWVertex_->init();
// now set up the decay modes
DecayPhaseSpaceModePtr mode;
tPDVector extpart(3);
vector<double> wgt(0);
// W modes
extpart[0]=getParticleData(ParticleID::Wplus);
// loop for the quarks
unsigned int iz=0;
for(int ix=1;ix<6;ix+=2) {
for(int iy=2;iy<6;iy+=2) {
// check that the combination of particles is allowed
if(!FFWVertex_->allowed(-ix,iy,ParticleID::Wminus))
throw InitException() << "SMWDecayer::doinit() the W vertex"
<< "cannot handle all the quark modes"
<< Exception::abortnow;
extpart[1] = getParticleData(-ix);
extpart[2] = getParticleData( iy);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
addMode(mode,quarkWeight_[iz],wgt);
++iz;
}
}
// loop for the leptons
for(int ix=11;ix<17;ix+=2) {
// check that the combination of particles is allowed
// if(!FFWVertex_->allowed(-ix,ix+1,ParticleID::Wminus))
// throw InitException() << "SMWDecayer::doinit() the W vertex"
// << "cannot handle all the lepton modes"
// << Exception::abortnow;
extpart[1] = getParticleData(-ix);
extpart[2] = getParticleData(ix+1);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
addMode(mode,leptonWeight_[(ix-11)/2],wgt);
}
gluon_ = getParticleData(ParticleID::g);
}
int SMWDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
int imode(-1);
if(children.size()!=2) return imode;
int id0=parent->id();
tPDVector::const_iterator pit = children.begin();
int id1=(**pit).id();
++pit;
int id2=(**pit).id();
if(abs(id0)!=ParticleID::Wplus) return imode;
int idd(0),idu(0);
if(abs(id1)%2==1&&abs(id2)%2==0) {
idd=abs(id1);
idu=abs(id2);
}
else if(abs(id1)%2==0&&abs(id2)%2==1) {
idd=abs(id2);
idu=abs(id1);
}
if(idd==0&&idu==0) {
return imode;
}
else if(idd<=5) {
imode=idd+idu/2-2;
}
else {
imode=(idd-1)/2+1;
}
cc= (id0==ParticleID::Wminus);
return imode;
}
void SMWDecayer::persistentOutput(PersistentOStream & os) const {
os << FFWVertex_ << quarkWeight_ << leptonWeight_
- << FFGVertex_ << gluon_ << ounit( pTmin_, GeV ) << NLO_;
+ << FFGVertex_ << gluon_ << NLO_;
}
void SMWDecayer::persistentInput(PersistentIStream & is, int) {
is >> FFWVertex_ >> quarkWeight_ >> leptonWeight_
- >> FFGVertex_ >> gluon_ >> iunit( pTmin_, GeV ) >> NLO_;
+ >> FFGVertex_ >> gluon_ >> NLO_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SMWDecayer,PerturbativeDecayer>
describeHerwigSMWDecayer("Herwig::SMWDecayer", "HwPerturbativeDecay.so");
void SMWDecayer::Init() {
static ClassDocumentation<SMWDecayer> documentation
("The SMWDecayer class is the implementation of the decay"
" of the W boson to the Standard Model fermions.");
static ParVector<SMWDecayer,double> interfaceWquarkMax
("QuarkMax",
"The maximum weight for the decay of the W to quarks",
&SMWDecayer::quarkWeight_,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<SMWDecayer,double> interfaceWleptonMax
("LeptonMax",
"The maximum weight for the decay of the W to leptons",
&SMWDecayer::leptonWeight_,
0, 0, 0, -10000, 10000, false, false, true);
- static Parameter<SMWDecayer, Energy> interfacePtMin
- ("minpT",
- "The pt cut on hardest emision generation",
- &SMWDecayer::pTmin_, GeV, 1.*GeV, 0*GeV, 100000.0*GeV,
- false, false, Interface::limited);
-
static Switch<SMWDecayer,bool> interfaceNLO
("NLO",
"Whether to return the LO or NLO result",
&SMWDecayer::NLO_, false, false, false);
static SwitchOption interfaceNLOLO
(interfaceNLO,
"No",
"Leading-order result",
false);
static SwitchOption interfaceNLONLO
(interfaceNLO,
"Yes",
"NLO result",
true);
}
// return the matrix element squared
double SMWDecayer::me2(const int, const Particle & part,
const ParticleVector & decay,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin1Half,PDT::Spin1Half)));
int iferm(1),ianti(0);
if(decay[0]->id()>0) swap(iferm,ianti);
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(vectors_,rho_,
const_ptr_cast<tPPtr>(&part),
incoming,false);
}
if(meopt==Terminate) {
VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast<tPPtr>(&part),
incoming,true,false);
SpinorBarWaveFunction::
constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
SpinorWaveFunction::
constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
return 0.;
}
SpinorBarWaveFunction::
calculateWaveFunctions(wavebar_,decay[iferm],outgoing);
SpinorWaveFunction::
calculateWaveFunctions(wave_ ,decay[ianti],outgoing);
// compute the matrix element
Energy2 scale(sqr(part.mass()));
for(unsigned int ifm=0;ifm<2;++ifm) {
for(unsigned int ia=0;ia<2;++ia) {
for(unsigned int vhel=0;vhel<3;++vhel) {
if(iferm>ianti) (*ME())(vhel,ia,ifm)=
FFWVertex_->evaluate(scale,wave_[ia],wavebar_[ifm],vectors_[vhel]);
else (*ME())(vhel,ifm,ia)=
FFWVertex_->evaluate(scale,wave_[ia],wavebar_[ifm],vectors_[vhel]);
}
}
}
double output=(ME()->contract(rho_)).real()*UnitRemoval::E2/scale;
if(abs(decay[0]->id())<=6) output*=3.;
if(decay[0]->hasColour()) decay[0]->antiColourNeighbour(decay[1]);
else if(decay[1]->hasColour()) decay[1]->antiColourNeighbour(decay[0]);
// leading-order result
if(!NLO_) return output;
// check decay products coloured, otherwise return
if(!decay[0]->dataPtr()->coloured()) return output;
// inital masses, couplings etc
// W mass
mW_ = part.mass();
// strong coupling
aS_ = SM().alphaS(sqr(mW_));
// reduced mass
double mu1 = (decay[0]->dataPtr()->mass())/mW_;
double mu2 = (decay[1]->dataPtr()->mass())/mW_;
// scale
scale_ = sqr(mW_);
// now for the nlo loop correction
double virt = CF_*aS_/Constants::pi;
// now for the real correction
double realFact=0.;
for(int iemit=0;iemit<2;++iemit) {
double phi = UseRandom::rnd()*Constants::twopi;
// set the emitter and the spectator
double muj = iemit==0 ? mu1 : mu2;
double muk = iemit==0 ? mu2 : mu1;
double muj2 = sqr(muj);
double muk2 = sqr(muk);
// calculate y
double yminus = 0.;
double yplus = 1.-2.*muk*(1.-muk)/(1.-muj2-muk2);
double y = yminus + UseRandom::rnd()*(yplus-yminus);
double v = sqrt(sqr(2.*muk2 + (1.-muj2-muk2)*(1.-y))-4.*muk2)
/(1.-muj2-muk2)/(1.-y);
double zplus = (1.+v)*(1.-muj2-muk2)*y/2./(muj2+(1.-muj2-muk2)*y);
double zminus = (1.-v)*(1.-muj2-muk2)*y/2./(muj2+(1.-muj2-muk2)*y);
double z = zminus + UseRandom::rnd()*(zplus-zminus);
double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
// calculate x1,x2,x3,xT
double x2 = 1.-y*(1.-muj2-muk2)-muj2+muk2;
double x1 = 1.+muj2-muk2-z*(x2-2.*muk2);
// copy the particle objects over for calculateRealEmission
vector<PPtr> hardProcess(3);
hardProcess[0] = const_ptr_cast<PPtr>(&part);
hardProcess[1] = decay[0];
hardProcess[2] = decay[1];
realFact += 0.25*jac*sqr(1.-muj2-muk2)/
sqrt((1.-sqr(muj-muk))*(1.-sqr(muj+muk)))/Constants::twopi
*2.*CF_*aS_*calculateRealEmission(x1, x2, hardProcess, phi,
muj, muk, iemit, true);
}
// the born + virtual + real
output *= (1. + virt + realFact);
return output;
}
void SMWDecayer::doinitrun() {
PerturbativeDecayer::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<numberModes();++ix) {
if(ix<6) quarkWeight_ [ix]=mode(ix)->maxWeight();
else leptonWeight_[ix-6]=mode(ix)->maxWeight();
}
}
}
void SMWDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
for(unsigned int ix=0;ix<quarkWeight_.size();++ix) {
output << "newdef " << name() << ":QuarkMax " << ix << " "
<< quarkWeight_[ix] << "\n";
}
for(unsigned int ix=0;ix<leptonWeight_.size();++ix) {
output << "newdef " << name() << ":LeptonMax " << ix << " "
<< leptonWeight_[ix] << "\n";
}
// parameters for the PerturbativeDecayer base class
PerturbativeDecayer::dataBaseOutput(output,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
void SMWDecayer::
initializeMECorrection(RealEmissionProcessPtr born, double & initial,
double & final) {
// get the quark and antiquark
ParticleVector qq;
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix)
qq.push_back(born->bornOutgoing()[ix]);
// ensure quark first
if(qq[0]->id()<0) swap(qq[0],qq[1]);
// centre of mass energy
d_Q_ = (qq[0]->momentum() + qq[1]->momentum()).m();
// quark mass
d_m_ = 0.5*(qq[0]->momentum().m()+qq[1]->momentum().m());
// set the other parameters
setRho(sqr(d_m_/d_Q_));
setKtildeSymm();
// otherwise can do it
initial=1.;
final =1.;
}
RealEmissionProcessPtr SMWDecayer::
applyHardMatrixElementCorrection(RealEmissionProcessPtr born) {
// get the quark and antiquark
ParticleVector qq;
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix)
qq.push_back(born->bornOutgoing()[ix]);
if(!qq[0]->dataPtr()->coloured()) return RealEmissionProcessPtr();
// ensure quark first
bool order = qq[0]->id()<0;
if(order) swap(qq[0],qq[1]);
// get the momenta
vector<Lorentz5Momentum> newfs = applyHard(qq);
// return if no emission
if(newfs.size()!=3) return RealEmissionProcessPtr();
// perform final check to ensure energy greater than constituent mass
for (int i=0; i<2; i++) {
if (newfs[i].e() < qq[i]->data().constituentMass()) return RealEmissionProcessPtr();
}
if (newfs[2].e() < getParticleData(ParticleID::g)->constituentMass())
return RealEmissionProcessPtr();
// set masses
for (int i=0; i<2; i++) newfs[i].setMass(qq[i]->mass());
newfs[2].setMass(ZERO);
// decide which particle emits
bool firstEmits=
newfs[2].vect().perp2(newfs[0].vect())<
newfs[2].vect().perp2(newfs[1].vect());
// create the new quark, antiquark and gluon
PPtr newg = getParticleData(ParticleID::g)->produceParticle(newfs[2]);
PPtr newq = qq[0]->dataPtr()->produceParticle(newfs[0]);
PPtr newa = qq[1]->dataPtr()->produceParticle(newfs[1]);
// create the output real emission process
for(unsigned int ix=0;ix<born->bornIncoming().size();++ix) {
born->incoming().push_back(born->bornIncoming()[ix]);
}
if(!order) {
born->outgoing().push_back(newq);
born->outgoing().push_back(newa);
born->outgoing().push_back(newg);
}
else {
born->outgoing().push_back(newa);
born->outgoing().push_back(newq);
born->outgoing().push_back(newg);
firstEmits = !firstEmits;
}
// make colour connections
newg->colourNeighbour(newq);
newa->colourNeighbour(newg);
if(firstEmits) {
born->emitter(1);
born->spectator(2);
}
else {
born->emitter(2);
born->spectator(1);
}
born->emitted(3);
born->interaction(ShowerInteraction::QCD);
return born;
}
vector<Lorentz5Momentum> SMWDecayer::
applyHard(const ParticleVector &p) {
double x, xbar;
vector<Lorentz5Momentum> fs;
// return if no emission
if (getHard(x, xbar) < UseRandom::rnd() || p.size() != 2) return fs;
// centre of mass energy
Lorentz5Momentum pcm = p[0]->momentum() + p[1]->momentum();
// momenta of quark,antiquark and gluon
Lorentz5Momentum pq, pa, pg;
if (p[0]->id() > 0) {
pq = p[0]->momentum();
pa = p[1]->momentum();
} else {
pa = p[0]->momentum();
pq = p[1]->momentum();
}
// boost to boson rest frame
Boost beta = (pcm.findBoostToCM());
pq.boost(beta);
pa.boost(beta);
// return if fails ?????
double xg = 2.-x-xbar;
if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return fs;
Axis u1, u2, u3;
// moduli of momenta in units of Q and cos theta
// stick to q direction?
// p1 is the one that is kept, p2 is the other fermion, p3 the gluon.
Energy e1, e2, e3;
Energy pp1, pp2, pp3;
bool keepq = true;
if (UseRandom::rnd() > sqr(x)/(sqr(x)+sqr(xbar)))
keepq = false;
if (keepq) {
pp1 = d_Q_*sqrt(sqr(x)-4.*d_rho_)/2.;
pp2 = d_Q_*sqrt(sqr(xbar)-4.*d_rho_)/2.;
e1 = d_Q_*x/2.;
e2 = d_Q_*xbar/2.;
u1 = pq.vect().unit();
} else {
pp2 = d_Q_*sqrt(sqr(x)-4.*d_rho_)/2.;
pp1 = d_Q_*sqrt(sqr(xbar)-4.*d_rho_)/2.;
e2 = d_Q_*x/2.;
e1 = d_Q_*xbar/2.;
u1 = pa.vect().unit();
}
pp3 = d_Q_*xg/2.;
e3 = pp3;
u2 = u1.orthogonal();
u2 /= u2.mag();
u3 = u1.cross(u2);
u3 /= u3.mag();
double ct2=-2., ct3=-2.;
if (pp1 == ZERO || pp2 == ZERO || pp3 == ZERO) {
bool touched = false;
if (pp1 == ZERO) {
ct2 = 1;
ct3 = -1;
touched = true;
}
if (pp2 == ZERO || pp3 == ZERO) {
ct2 = 1;
ct3 = 1;
touched = true;
}
if (!touched)
throw Exception() << "SMWDecayer::applyHard()"
<< " did not set ct2/3"
<< Exception::abortnow;
} else {
ct3 = (sqr(pp1)+sqr(pp3)-sqr(pp2))/(2.*pp1*pp3);
ct2 = (sqr(pp1)+sqr(pp2)-sqr(pp3))/(2.*pp1*pp2);
}
double phi = Constants::twopi*UseRandom::rnd();
double cphi = cos(phi);
double sphi = sin(phi);
double st2 = sqrt(1.-sqr(ct2));
double st3 = sqrt(1.-sqr(ct3));
ThreeVector<Energy> pv1, pv2, pv3;
pv1 = pp1*u1;
pv2 = -ct2*pp2*u1 + st2*cphi*pp2*u2 + st2*sphi*pp2*u3;
pv3 = -ct3*pp3*u1 - st3*cphi*pp3*u2 - st3*sphi*pp3*u3;
if (keepq) {
pq = Lorentz5Momentum(pv1, e1);
pa = Lorentz5Momentum(pv2, e2);
} else {
pa = Lorentz5Momentum(pv1, e1);
pq = Lorentz5Momentum(pv2, e2);
}
pg = Lorentz5Momentum(pv3, e3);
pq.boost(-beta);
pa.boost(-beta);
pg.boost(-beta);
fs.push_back(pq);
fs.push_back(pa);
fs.push_back(pg);
return fs;
}
double SMWDecayer::getHard(double &x1, double &x2) {
double w = 0.0;
double y1 = UseRandom::rnd(),y2 = UseRandom::rnd();
// simply double MC efficiency
// -> weight has to be divided by two (Jacobian)
if (y1 + y2 > 1) {
y1 = 1.-y1;
y2 = 1.-y2;
}
bool inSoft = false;
if (y1 < 0.25) {
if (y2 < 0.25) {
inSoft = true;
if (y1 < y2) {
y1 = 0.25-y1;
y2 = y1*(1.5 - 2.*y2);
} else {
y2 = 0.25 - y2;
y1 = y2*(1.5 - 2.*y1);
}
} else {
if (y2 < y1 + 2.*sqr(y1)) return w;
}
} else {
if (y2 < 0.25) {
if (y1 < y2 + 2.*sqr(y2)) return w;
}
}
// inside PS?
x1 = 1.-y1;
x2 = 1.-y2;
if(y1*y2*(1.-y1-y2) < d_rho_*sqr(y1+y2)) return w;
double k1 = getKfromX(x1, x2);
double k2 = getKfromX(x2, x1);
// Is it in the quark emission zone?
if (k1 < d_kt1_) return 0.0;
// No...is it in the anti-quark emission zone?
if (k2 < d_kt2_) return 0.0;
// Point is in dead zone: compute q qbar g weight
w = MEV(x1, x2);
// for axial:
// w = MEA(x1, x2);
// Reweight soft region
if (inSoft) {
if (y1 < y2) w *= 2.*y1;
else w *= 2.*y2;
}
// alpha and colour factors
Energy2 pt2 = sqr(d_Q_)*(1.-x1)*(1.-x2);
w *= 1./3./Constants::pi*coupling()->value(pt2);
return w;
}
bool SMWDecayer::
softMatrixElementVeto(ShowerProgenitorPtr initial,ShowerParticlePtr parent,Branching br) {
// check we should be applying the veto
if(parent->id()!=initial->progenitor()->id()||
br.ids[0]!=br.ids[1]||
br.ids[2]->id()!=ParticleID::g) return false;
// calculate pt
double d_z = br.kinematics->z();
Energy d_qt = br.kinematics->scale();
Energy2 d_m2 = parent->momentum().m2();
Energy2 pPerp2 = sqr(d_z*d_qt) - d_m2;
if(pPerp2<ZERO) {
parent->vetoEmission(br.type,br.kinematics->scale());
return true;
}
Energy pPerp = (1.-d_z)*sqrt(pPerp2);
// if not hardest so far don't apply veto
if(pPerp<initial->highestpT()) return false;
// calculate the weight
double weight = 0.;
if(parent->id()>0) weight = qWeightX(d_qt, d_z);
else weight = qbarWeightX(d_qt, d_z);
// compute veto from weight
bool veto = !UseRandom::rndbool(weight);
// if vetoing reset the scale
if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
// return the veto
return veto;
}
void SMWDecayer::setRho(double r)
{
d_rho_ = r;
d_v_ = sqrt(1.-4.*d_rho_);
}
void SMWDecayer::setKtildeSymm() {
d_kt1_ = (1. + sqrt(1. - 4.*d_rho_))/2.;
setKtilde2();
}
void SMWDecayer::setKtilde2() {
double num = d_rho_ * d_kt1_ + 0.25 * d_v_ *(1.+d_v_)*(1.+d_v_);
double den = d_kt1_ - d_rho_;
d_kt2_ = num/den;
}
double SMWDecayer::getZfromX(double x1, double x2) {
double uval = u(x2);
double num = x1 - (2. - x2)*uval;
double den = sqrt(x2*x2 - 4.*d_rho_);
return uval + num/den;
}
double SMWDecayer::getKfromX(double x1, double x2) {
double zval = getZfromX(x1, x2);
return (1.-x2)/(zval*(1.-zval));
}
double SMWDecayer::MEV(double x1, double x2) {
// Vector part
double num = (x1+2.*d_rho_)*(x1+2.*d_rho_) + (x2+2.*d_rho_)*(x2+2.*d_rho_)
- 8.*d_rho_*(1.+2.*d_rho_);
double den = (1.+2.*d_rho_)*(1.-x1)*(1.-x2);
return (num/den - 2.*d_rho_/((1.-x1)*(1.-x1))
- 2*d_rho_/((1.-x2)*(1.-x2)))/d_v_;
}
double SMWDecayer::MEA(double x1, double x2) {
// Axial part
double num = (x1+2.*d_rho_)*(x1+2.*d_rho_) + (x2+2.*d_rho_)*(x2+2.*d_rho_)
+ 2.*d_rho_*((5.-x1-x2)*(5.-x1-x2) - 19.0 + 4*d_rho_);
double den = d_v_*d_v_*(1.-x1)*(1.-x2);
return (num/den - 2.*d_rho_/((1.-x1)*(1.-x1))
- 2*d_rho_/((1.-x2)*(1.-x2)))/d_v_;
}
double SMWDecayer::u(double x2) {
return 0.5*(1. + d_rho_/(1.-x2+d_rho_));
}
void SMWDecayer::
getXXbar(double kti, double z, double &x, double &xbar) {
double w = sqr(d_v_) + kti*(-1. + z)*z*(2. + kti*(-1. + z)*z);
if (w < 0) {
x = -1.;
xbar = -1;
} else {
x = (1. + sqr(d_v_)*(-1. + z) + sqr(kti*(-1. + z))*z*z*z
+ z*sqrt(w)
- kti*(-1. + z)*z*(2. + z*(-2 + sqrt(w))))/
(1. - kti*(-1. + z)*z + sqrt(w));
xbar = 1. + kti*(-1. + z)*z;
}
}
double SMWDecayer::qWeight(double x, double xbar) {
double rval;
double xg = 2. - xbar - x;
// always return one in the soft gluon region
if(xg < EPS_) return 1.0;
// check it is in the phase space
if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return 0.0;
double k1 = getKfromX(x, xbar);
double k2 = getKfromX(xbar, x);
// Is it in the quark emission zone?
if(k1 < d_kt1_) {
rval = MEV(x, xbar)/PS(x, xbar);
// is it also in the anti-quark emission zone?
if(k2 < d_kt2_) rval *= 0.5;
return rval;
}
return 1.0;
}
double SMWDecayer::qbarWeight(double x, double xbar) {
double rval;
double xg = 2. - xbar - x;
// always return one in the soft gluon region
if(xg < EPS_) return 1.0;
// check it is in the phase space
if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return 0.0;
double k1 = getKfromX(x, xbar);
double k2 = getKfromX(xbar, x);
// Is it in the antiquark emission zone?
if(k2 < d_kt2_) {
rval = MEV(x, xbar)/PS(xbar, x);
// is it also in the quark emission zone?
if(k1 < d_kt1_) rval *= 0.5;
return rval;
}
return 1.0;
}
double SMWDecayer::qWeightX(Energy qtilde, double z) {
double x, xb;
getXXbar(sqr(qtilde/d_Q_), z, x, xb);
// if exceptionally out of phase space, leave this emission, as there
// is no good interpretation for the soft ME correction.
if (x < 0 || xb < 0) return 1.0;
return qWeight(x, xb);
}
double SMWDecayer::qbarWeightX(Energy qtilde, double z) {
double x, xb;
getXXbar(sqr(qtilde/d_Q_), z, xb, x);
// see above in qWeightX.
if (x < 0 || xb < 0) return 1.0;
return qbarWeight(x, xb);
}
double SMWDecayer::PS(double x, double xbar) {
double u = 0.5*(1. + d_rho_ / (1.-xbar+d_rho_));
double z = u + (x - (2.-xbar)*u)/sqrt(xbar*xbar - 4.*d_rho_);
double brack = (1.+z*z)/(1.-z)- 2.*d_rho_/(1-xbar);
// interesting: the splitting function without the subtraction
// term. Actually gives a much worse approximation in the collinear
// limit. double brack = (1.+z*z)/(1.-z);
double den = (1.-xbar)*sqrt(xbar*xbar - 4.*d_rho_);
return brack/den;
}
double SMWDecayer::matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
const ParticleVector & decay3, MEOption) {
// extract partons and LO momentas
vector<cPDPtr> partons(1,inpart.dataPtr());
vector<Lorentz5Momentum> lomom(1,inpart.momentum());
for(unsigned int ix=0;ix<2;++ix) {
partons.push_back(decay2[ix]->dataPtr());
lomom.push_back(decay2[ix]->momentum());
}
vector<Lorentz5Momentum> realmom(1,inpart.momentum());
for(unsigned int ix=0;ix<3;++ix) {
if(ix==2) partons.push_back(decay3[ix]->dataPtr());
realmom.push_back(decay3[ix]->momentum());
}
if(partons[0]->id()<0) {
swap(partons[1],partons[2]);
swap(lomom[1],lomom[2]);
swap(realmom[1],realmom[2]);
}
double lome = loME(partons,lomom);
InvEnergy2 reme = realME(partons,realmom);
double ratio = CF_*reme/lome*sqr(inpart.mass());
return ratio;
}
double SMWDecayer::meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter, bool subtract) const {
Lorentz5Momentum q = momenta[1]+momenta[2]+momenta[3];
Energy2 Q2=q.m2();
Energy2 lambda = sqrt((Q2-sqr(momenta[1].mass()+momenta[2].mass()))*
(Q2-sqr(momenta[1].mass()-momenta[2].mass())));
InvEnergy2 D[2];
double lome(0.);
for(unsigned int iemit=0;iemit<2;++iemit) {
unsigned int ispect = iemit==0 ? 1 : 0;
Energy2 pipj = momenta[3 ] * momenta[1+iemit ];
Energy2 pipk = momenta[3 ] * momenta[1+ispect];
Energy2 pjpk = momenta[1+iemit] * momenta[1+ispect];
double y = pipj/(pipj+pipk+pjpk);
double z = pipk/( pipk+pjpk);
Energy mij = sqrt(2.*pipj+sqr(momenta[1+iemit].mass()));
Energy2 lamB = sqrt((Q2-sqr(mij+momenta[1+ispect].mass()))*
(Q2-sqr(mij-momenta[1+ispect].mass())));
Energy2 Qpk = q*momenta[1+ispect];
Lorentz5Momentum pkt =
lambda/lamB*(momenta[1+ispect]-Qpk/Q2*q)
+0.5/Q2*(Q2+sqr(momenta[1+ispect].mass())-sqr(momenta[1+ispect].mass()))*q;
Lorentz5Momentum pijt =
q-pkt;
double muj = momenta[1+iemit ].mass()/sqrt(Q2);
double muk = momenta[1+ispect].mass()/sqrt(Q2);
double vt = sqrt((1.-sqr(muj+muk))*(1.-sqr(muj-muk)))/(1.-sqr(muj)-sqr(muk));
double v = sqrt(sqr(2.*sqr(muk)+(1.-sqr(muj)-sqr(muk))*(1.-y))-4.*sqr(muk))
/(1.-y)/(1.-sqr(muj)-sqr(muk));
// dipole term
D[iemit] = 0.5/pipj*(2./(1.-(1.-z)*(1.-y))
-vt/v*(2.-z+sqr(momenta[1+iemit].mass())/pipj));
// matrix element
vector<Lorentz5Momentum> lomom(3);
lomom[0] = momenta[0];
if(iemit==0) {
lomom[1] = pijt;
lomom[2] = pkt ;
}
else {
lomom[2] = pijt;
lomom[1] = pkt ;
}
if(iemit==0) lome = loME(partons,lomom);
}
InvEnergy2 ratio = realME(partons,momenta)/lome*abs(D[iemitter])
/(abs(D[0])+abs(D[1]));
if(subtract)
return Q2*(ratio-2.*D[iemitter]);
else
return Q2*ratio;
}
double SMWDecayer::loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
VectorWaveFunction win (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
}
for(unsigned int ix=0;ix<3;++ix){
win.reset(ix);
vin.push_back(win);
}
// temporary storage of the different diagrams
// sum over helicities to get the matrix element
double total(0.);
for(unsigned int inhel=0;inhel<3;++inhel) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
Complex diag1 = FFWVertex()->evaluate(scale_,aout[outhel2],fout[outhel1],vin[inhel]);
total += norm(diag1);
}
}
}
// return the answer
return total;
}
InvEnergy2 SMWDecayer::realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
vector<VectorWaveFunction> gout;
VectorWaveFunction win (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
VectorWaveFunction gluon(momenta[3],partons[3],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
gluon.reset(2*ix);
gout.push_back(gluon);
}
for(unsigned int ix=0;ix<3;++ix){
win.reset(ix);
vin.push_back(win);
}
vector<Complex> diag(2,0.);
double total(0.);
for(unsigned int inhel1=0;inhel1<3;++inhel1) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
for(unsigned int outhel3=0;outhel3<2;++outhel3) {
SpinorBarWaveFunction off1 =
FFGVertex()->evaluate(scale_,3,partons[1],fout[outhel1],gout[outhel3]);
diag[0] = FFWVertex()->evaluate(scale_,aout[outhel2],off1,vin[inhel1]);
SpinorWaveFunction off2 =
FFGVertex()->evaluate(scale_,3,partons[2],aout[outhel2],gout[outhel3]);
diag[1] = FFWVertex()->evaluate(scale_,off2,fout[outhel1],vin[inhel1]);
// sum of diagrams
Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.));
// me2
total += norm(sum);
}
}
}
}
// divide out the coupling
total /= norm(FFGVertex()->norm());
// return the total
return total*UnitRemoval::InvE2;
}
double SMWDecayer::calculateRealEmission(double x1, double x2,
vector<PPtr> hardProcess,
double phi, double muj,
double muk, int iemit,
bool subtract) const {
// make partons data object for meRatio
vector<cPDPtr> partons (3);
for(int ix=0; ix<3; ++ix)
partons[ix] = hardProcess[ix]->dataPtr();
partons.push_back(gluon_);
// calculate x3
double x3 = 2.-x1-x2;
double xT = sqrt(max(0.,sqr(x3)-0.25*sqr(sqr(x2)+sqr(x3)-sqr(x1)-4.*sqr(muk)+4.*sqr(muj))
/(sqr(x2)-4.*sqr(muk))));
// calculate the momenta
Energy M = mW_;
Lorentz5Momentum pspect(ZERO,ZERO,-0.5*M*sqrt(max(sqr(x2)-4.*sqr(muk),0.)),
0.5*M*x2,M*muk);
Lorentz5Momentum pemit (-0.5*M*xT*cos(phi),-0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x1)-sqr(xT)-4.*sqr(muj),0.)),
0.5*M*x1,M*muj);
Lorentz5Momentum pgluon(0.5*M*xT*cos(phi), 0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x3)-sqr(xT),0.)),0.5*M*x3,ZERO);
if(abs(pspect.z()+pemit.z()-pgluon.z())/M<1e-6)
pgluon.setZ(-pgluon.z());
else if(abs(pspect.z()-pemit.z()+pgluon.z())/M<1e-6)
pemit .setZ(- pemit.z());
// boost and rotate momenta
LorentzRotation eventFrame( ( hardProcess[1]->momentum() +
hardProcess[2]->momentum() ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*hardProcess[iemit+1]->momentum();
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() );
eventFrame.invert();
vector<Lorentz5Momentum> momenta(3);
momenta[0] = hardProcess[0]->momentum();
if(iemit==0) {
momenta[2] = eventFrame*pspect;
momenta[1] = eventFrame*pemit ;
}
else {
momenta[1] = eventFrame*pspect;
momenta[2] = eventFrame*pemit ;
}
momenta.push_back(eventFrame*pgluon);
// calculate the weight
double realwgt(0.);
if(1.-x1>1e-5 && 1.-x2>1e-5)
realwgt = meRatio(partons,momenta,iemit,subtract);
return realwgt;
}
diff --git a/Decay/Perturbative/SMWDecayer.h b/Decay/Perturbative/SMWDecayer.h
--- a/Decay/Perturbative/SMWDecayer.h
+++ b/Decay/Perturbative/SMWDecayer.h
@@ -1,520 +1,511 @@
// -*- C++ -*-
//
// SMWDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SMWDecayer_H
#define HERWIG_SMWDecayer_H
//
// This is the declaration of the SMWDecayer class.
//
#include "Herwig/Decay/PerturbativeDecayer.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "Herwig/Decay/DecayPhaseSpaceMode.h"
-#include "Herwig/Shower/Core/Couplings/ShowerAlpha.fh"
namespace Herwig {
using namespace ThePEG;
using namespace ThePEG::Helicity;
/** \ingroup Decay
*
* The <code>SMWDecayer</code> is designed to perform the decay of the
* W boson to the Standard Model fermions, including the first order
* electroweak corrections.
*
* @see PerturbativeDecayer
*
*/
class SMWDecayer: public PerturbativeDecayer {
public:
/**
* Default constructor.
*/
SMWDecayer();
public:
/**
* Virtual members to be overridden by inheriting classes
* which implement hard corrections
*/
//@{
/**
* Has an old fashioned ME correction
*/
virtual bool hasMECorrection() {return true;}
/**
* Initialize the ME correction
*/
virtual void initializeMECorrection(RealEmissionProcessPtr , double & ,
double & );
/**
* Apply the hard matrix element correction to a given hard process or decay
*/
virtual RealEmissionProcessPtr applyHardMatrixElementCorrection(RealEmissionProcessPtr);
/**
* Apply the soft matrix element correction
* @param initial The particle from the hard process which started the
* shower
* @param parent The initial particle in the current branching
* @param br The branching struct
* @return If true the emission should be vetoed
*/
virtual bool softMatrixElementVeto(ShowerProgenitorPtr initial,
ShowerParticlePtr parent,Branching br);
/**
* Virtual members to be overridden by inheriting classes
* which implement hard corrections
*/
//@{
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
//@}
public:
/**
* Which of the possible decays is required
* @param cc Is this mode the charge conjugate
* @param parent The decaying particle
* @param children The decay products
*/
virtual int modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const;
/**
* Return the matrix element squared for a given mode and phase-space channel.
* @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
* @param decay The particles produced in the decay.
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
virtual double me2(const int ichan, const Particle & part,
const ParticleVector & decay,MEOption meopt) const;
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
*/
virtual void dataBaseOutput(ofstream & os,bool header) 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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
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 {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving and
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/**
* Apply the hard matrix element
*/
vector<Lorentz5Momentum> applyHard(const ParticleVector &p);
/**
* Get the weight for hard emission
*/
double getHard(double &, double &);
/**
* Set the \f$\rho\f$ parameter
*/
void setRho(double);
/**
* Set the \f$\tilde{\kappa}\f$ parameters symmetrically
*/
void setKtildeSymm();
/**
* Set second \f$\tilde{\kappa}\f$, given the first.
*/
void setKtilde2();
/**
* Translate the variables from \f$x_q,x_{\bar{q}}\f$ to \f$\tilde{\kappa},z\f$
*/
//@{
/**
* Calculate \f$z\f$.
*/
double getZfromX(double, double);
/**
* Calculate \f$\tilde{\kappa}\f$.
*/
double getKfromX(double, double);
//@}
/**
* Calculate \f$x_{q},x_{\bar{q}}\f$ from \f$\tilde{\kappa},z\f$.
* @param kt \f$\tilde{\kappa}\f$
* @param z \f$z\f$
* @param x \f$x_{q}\f$
* @param xbar \f$x_{\bar{q}}\f$
*/
void getXXbar(double kt, double z, double & x, double & xbar);
/**
* Soft weight
*/
//@{
/**
* Soft quark weight calculated from \f$x_{q},x_{\bar{q}}\f$
* @param x \f$x_{q}\f$
* @param xbar \f$x_{\bar{q}}\f$
*/
double qWeight(double x, double xbar);
/**
* Soft antiquark weight calculated from \f$x_{q},x_{\bar{q}}\f$
* @param x \f$x_{q}\f$
* @param xbar \f$x_{\bar{q}}\f$
*/
double qbarWeight(double x, double xbar);
/**
* Soft quark weight calculated from \f$\tilde{q},z\f$
* @param qtilde \f$\tilde{q}\f$
* @param z \f$z\f$
*/
double qWeightX(Energy qtilde, double z);
/**
* Soft antiquark weight calculated from \f$\tilde{q},z\f$
* @param qtilde \f$\tilde{q}\f$
* @param z \f$z\f$
*/
double qbarWeightX(Energy qtilde, double z);
//@}
/**
* ????
*/
double u(double);
/**
* Vector and axial vector parts of the matrix element
*/
//@{
/**
* Vector part of the matrix element
*/
double MEV(double, double);
/**
* Axial vector part of the matrix element
*/
double MEA(double, double);
/**
* The matrix element, given \f$x_1\f$, \f$x_2\f$.
* @param x1 \f$x_1\f$
* @param x2 \f$x_2\f$
*/
double PS(double x1, double x2);
//@}
protected:
/**
* Pointer to the fermion-antifermion W vertex
*/
AbstractFFVVertexPtr FFWVertex() const {return FFWVertex_;}
/**
* Pointer to the fermion-antifermion G vertex
*/
AbstractFFVVertexPtr FFGVertex() const {return FFGVertex_;}
/**
* Real emission term, for use in generating the hardest emission
*/
double calculateRealEmission(double x1, double x2,
vector<PPtr> hardProcess,
double phi, double muj, double muk,
int iemit, bool subtract) const;
/**
* Calculate the ratio between NLO & LO ME
*/
double meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter,bool subtract) const;
/**
* Calculate matrix element ratio R/B
*/
virtual double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
const ParticleVector & decay3, MEOption meopt);
/**
* Calculate the LO ME
*/
double loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const;
/**
* Calculate the NLO real emission piece of ME
*/
InvEnergy2 realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const;
private:
/**
* Private and non-existent assignment operator.
*/
SMWDecayer & operator=(const SMWDecayer &);
private:
/**
* Pointer to the fermion-antifermion W vertex
*/
AbstractFFVVertexPtr FFWVertex_;
/**
* Pointer to the fermion-antifermion G vertex
*/
AbstractFFVVertexPtr FFGVertex_;
/**
* maximum weights for the different integrations
*/
//@{
/**
* Weights for the W to quarks decays.
*/
vector<double> quarkWeight_;
/**
* Weights for the W to leptons decays.
*/
vector<double> leptonWeight_;
//@}
/**
* Spin density matrix for the decay
*/
mutable RhoDMatrix rho_;
/**
* Polarization vectors for the decay
*/
mutable vector<VectorWaveFunction> vectors_;
/**
* Spinors for the decay
*/
mutable vector<SpinorWaveFunction> wave_;
/**
* Barred spinors for the decay
*/
mutable vector<SpinorBarWaveFunction> wavebar_;
private:
/**
* CM energy
*/
Energy d_Q_;
/**
* Quark mass
*/
Energy d_m_;
/**
* The rho parameter
*/
double d_rho_;
/**
* The v parameter
*/
double d_v_;
/**
* The initial kappa-tilde values for radiation from the quark
*/
double d_kt1_;
/**
* The initial kappa-tilde values for radiation from the antiquark
*/
double d_kt2_;
/**
* Cut-off parameter
*/
static const double EPS_;
private:
/**
* The colour factor
*/
double CF_;
/**
* The W mass
*/
mutable Energy mW_;
// TODO: delete this
mutable double mu_;
/**
* The reduced mass of particle 1
*/
mutable double mu1_;
/**
* The reduced mass of particle 1 squared
*/
mutable double mu12_;
/**
* The reduceed mass of particle 2
*/
mutable double mu2_;
/**
* The reduceed mass of particle 2 squared
*/
mutable double mu22_;
/**
* The strong coupling
*/
mutable double aS_;
/**
* The scale
*/
mutable Energy2 scale_;
/**
* Stuff for the POWHEG correction
*/
//@{
/**
* ParticleData object for the gluon
*/
tcPDPtr gluon_;
/**
- * The cut off on pt, assuming massless quarks.
- */
- Energy pTmin_;
-
- // radiative variables (pt,y)
- Energy pT_;
-
- /**
* The ParticleData objects for the fermions
*/
vector<tcPDPtr> partons_;
/**
* The fermion momenta
*/
vector<Lorentz5Momentum> quark_;
/**
* The momentum of the radiated gauge boson
*/
Lorentz5Momentum gauge_;
/**
* The W boson
*/
PPtr wboson_;
/**
* W mass squared
*/
Energy2 mw2_;
//@}
/**
* Whether ro return the LO or NLO result
*/
bool NLO_;
};
}
#endif /* HERWIG_SMWDecayer_H */
diff --git a/Decay/Perturbative/SMZDecayer.cc b/Decay/Perturbative/SMZDecayer.cc
--- a/Decay/Perturbative/SMZDecayer.cc
+++ b/Decay/Perturbative/SMZDecayer.cc
@@ -1,1588 +1,1331 @@
// -*- C++ -*-
//
// SMZDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SMZDecayer class.
//
#include "SMZDecayer.h"
#include "Herwig/Utilities/Maths.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Decay/DecayVertex.h"
#include "ThePEG/Helicity/VectorSpinInfo.h"
#include "ThePEG/Helicity/FermionSpinInfo.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "Herwig/Shower/Core/Base/ShowerProgenitor.h"
#include "Herwig/Shower/Core/Base/ShowerParticle.h"
#include "Herwig/Shower/Core/Base/Branching.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include <numeric>
using namespace Herwig;
using namespace ThePEG::Helicity;
const double SMZDecayer::EPS_=0.00000001;
SMZDecayer::SMZDecayer()
- : quarkWeight_(5,0.), leptonWeight_(6,0.), CF_(4./3.), pTmin_(1.*GeV),
+ : quarkWeight_(5,0.), leptonWeight_(6,0.), CF_(4./3.),
NLO_(false) {
quarkWeight_[0] = 0.488029;
quarkWeight_[1] = 0.378461;
quarkWeight_[2] = 0.488019;
quarkWeight_[3] = 0.378027;
quarkWeight_[4] = 0.483207;
leptonWeight_[0] = 0.110709;
leptonWeight_[1] = 0.220276;
leptonWeight_[2] = 0.110708;
leptonWeight_[3] = 0.220276;
leptonWeight_[4] = 0.110458;
leptonWeight_[5] = 0.220276;
// intermediates
generateIntermediates(false);
// QED corrections
hasRealEmissionME(true);
hasOneLoopME(true);
}
void SMZDecayer::doinit() {
PerturbativeDecayer::doinit();
// get the vertices from the Standard Model object
tcHwSMPtr hwsm=dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException() << "Must have Herwig StandardModel object in"
<< "SMZDecayer::doinit()"
<< Exception::runerror;
// cast the vertices
FFZVertex_ = dynamic_ptr_cast<FFVVertexPtr>(hwsm->vertexFFZ());
FFZVertex_->init();
FFGVertex_ = hwsm->vertexFFG();
FFGVertex_->init();
FFPVertex_ = hwsm->vertexFFP();
FFPVertex_->init();
gluon_ = getParticleData(ParticleID::g);
// now set up the decay modes
DecayPhaseSpaceModePtr mode;
tPDVector extpart(3);
vector<double> wgt(0);
// the Z decay modes
extpart[0]=getParticleData(ParticleID::Z0);
// loop over the quarks and the leptons
for(int istep=0;istep<11;istep+=10) {
for(int ix=1;ix<7;++ix) {
int iy=istep+ix;
if(iy==6) continue;
extpart[1] = getParticleData(-iy);
extpart[2] = getParticleData( iy);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
if(iy<=6) addMode(mode, quarkWeight_.at(ix-1),wgt);
else addMode(mode,leptonWeight_.at(iy-11),wgt);
}
}
}
int SMZDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
int imode(-1);
if(children.size()!=2) return imode;
int id0=parent->id();
tPDVector::const_iterator pit = children.begin();
int id1=(**pit).id();
++pit;
int id2=(**pit).id();
// Z to quarks or leptons
cc =false;
if(id0!=ParticleID::Z0) return imode;
if(abs(id1)<6&&id1==-id2) {
imode=abs(id1)-1;
}
else if(abs(id1)>=11&&abs(id1)<=16&&id1==-id2) {
imode=abs(id1)-6;
}
cc = false;
return imode;
}
void SMZDecayer::persistentOutput(PersistentOStream & os) const {
os << FFZVertex_ << FFPVertex_ << FFGVertex_
- << quarkWeight_ << leptonWeight_ << alpha_ << NLO_
- << gluon_ << ounit( pTmin_, GeV );
+ << quarkWeight_ << leptonWeight_ << NLO_
+ << gluon_;
}
void SMZDecayer::persistentInput(PersistentIStream & is, int) {
is >> FFZVertex_ >> FFPVertex_ >> FFGVertex_
- >> quarkWeight_ >> leptonWeight_ >> alpha_ >> NLO_
- >> gluon_ >> iunit( pTmin_, GeV );
+ >> quarkWeight_ >> leptonWeight_ >> NLO_
+ >> gluon_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SMZDecayer,PerturbativeDecayer>
describeHerwigSMZDecayer("Herwig::SMZDecayer", "HwPerturbativeDecay.so");
void SMZDecayer::Init() {
static ClassDocumentation<SMZDecayer> documentation
("The SMZDecayer class is the implementation of the decay"
" Z boson to the Standard Model fermions.");
static ParVector<SMZDecayer,double> interfaceZquarkMax
("QuarkMax",
"The maximum weight for the decay of the Z to quarks",
&SMZDecayer::quarkWeight_,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<SMZDecayer,double> interfaceZleptonMax
("LeptonMax",
"The maximum weight for the decay of the Z to leptons",
&SMZDecayer::leptonWeight_,
0, 0, 0, -10000, 10000, false, false, true);
- static Reference<SMZDecayer,ShowerAlpha> interfaceCoupling
- ("Coupling",
- "Pointer to the object to calculate the coupling for the correction",
- &SMZDecayer::alpha_, false, false, true, false, false);
-
- static Parameter<SMZDecayer, Energy> interfacePtMin
- ("minpT",
- "The pt cut on hardest emision generation",
- &SMZDecayer::pTmin_, GeV, 1.*GeV, 0*GeV, 100000.0*GeV,
- false, false, Interface::limited);
-
static Switch<SMZDecayer,bool> interfaceNLO
("NLO",
"Whether to return the LO or NLO result",
&SMZDecayer::NLO_, false, false, false);
static SwitchOption interfaceNLOLO
(interfaceNLO,
"No",
"Leading-order result",
false);
static SwitchOption interfaceNLONLO
(interfaceNLO,
"Yes",
"NLO result",
true);
}
// return the matrix element squared
double SMZDecayer::me2(const int, const Particle & part,
const ParticleVector & decay,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin1Half,PDT::Spin1Half)));
int iferm(1),ianti(0);
if(decay[0]->id()>0) swap(iferm,ianti);
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(_vectors,_rho,
const_ptr_cast<tPPtr>(&part),
incoming,false);
}
if(meopt==Terminate) {
VectorWaveFunction::constructSpinInfo(_vectors,const_ptr_cast<tPPtr>(&part),
incoming,true,false);
SpinorBarWaveFunction::
constructSpinInfo(_wavebar,decay[iferm],outgoing,true);
SpinorWaveFunction::
constructSpinInfo(_wave ,decay[ianti],outgoing,true);
return 0.;
}
SpinorBarWaveFunction::
calculateWaveFunctions(_wavebar,decay[iferm],outgoing);
SpinorWaveFunction::
calculateWaveFunctions(_wave ,decay[ianti],outgoing);
// compute the matrix element
Energy2 scale(sqr(part.mass()));
unsigned int ifm,ia,vhel;
for(ifm=0;ifm<2;++ifm) {
for(ia=0;ia<2;++ia) {
for(vhel=0;vhel<3;++vhel) {
if(iferm>ianti) (*ME())(vhel,ia,ifm)=
FFZVertex_->evaluate(scale,_wave[ia],_wavebar[ifm],_vectors[vhel]);
else (*ME())(vhel,ifm,ia)=
FFZVertex_->evaluate(scale,_wave[ia],_wavebar[ifm],_vectors[vhel]);
}
}
}
double output=(ME()->contract(_rho)).real()*UnitRemoval::E2/scale;
if(abs(decay[0]->id())<=6) output*=3.;
if(decay[0]->hasColour()) decay[0]->antiColourNeighbour(decay[1]);
else if(decay[1]->hasColour()) decay[1]->antiColourNeighbour(decay[0]);
// if LO return
if(!NLO_) return output; // check decay products coloured, otherwise return
if(!decay[0]->dataPtr()->coloured()) return output;
// inital masses, couplings etc
// fermion mass
Energy particleMass = decay[0]->dataPtr()->mass();
// Z mass
mZ_ = part.mass();
// strong coupling
aS_ = SM().alphaS(sqr(mZ_));
// reduced mass
mu_ = particleMass/mZ_;
mu2_ = sqr(mu_);
// scale
scale_ = sqr(mZ_);
// compute the spinors
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
vector<VectorWaveFunction> vin;
SpinorBarWaveFunction qkout(decay[0]->momentum(),decay[0]->dataPtr(),outgoing);
SpinorWaveFunction qbout(decay[1]->momentum(),decay[1]->dataPtr(),outgoing);
VectorWaveFunction zin (part.momentum() ,part.dataPtr() ,incoming);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
}
for(unsigned int ix=0;ix<3;++ix){
zin.reset(ix);
vin.push_back(zin);
}
// temporary storage of the different diagrams
// sum over helicities to get the matrix element
double total=0.;
if(mu_!=0.) {
LorentzPolarizationVector momDiff =
(decay[0]->momentum()-decay[1]->momentum())/2./
(decay[0]->momentum().mass()+decay[1]->momentum().mass());
// scalars
Complex scalar1 = zin.wave().dot(momDiff);
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
for(unsigned int inhel=0;inhel<3;++inhel) {
// first the LO bit
Complex diag1 = FFZVertex_->evaluate(scale_,aout[outhel2],
fout[outhel1],vin[inhel]);
// extra stuff for NLO
LorentzPolarizationVector left =
aout[outhel2].wave().leftCurrent(fout[outhel1].wave());
LorentzPolarizationVector right =
aout[outhel2].wave().rightCurrent(fout[outhel1].wave());
Complex scalar =
aout[outhel2].wave().scalar(fout[outhel1].wave());
// nlo specific pieces
Complex diag3 =
Complex(0.,1.)*FFZVertex_->norm()*
(FFZVertex_->right()*( left.dot(zin.wave())) +
FFZVertex_-> left()*(right.dot(zin.wave())) -
( FFZVertex_-> left()+FFZVertex_->right())*scalar1*scalar);
// nlo piece
total += real(diag1*conj(diag3) + diag3*conj(diag1));
}
}
}
// rescale
total *= UnitRemoval::E2/scale_;
}
else {
total = ZERO;
}
// now for the NLO bit
double mu4 = sqr(mu2_);
double lmu = mu_!=0. ? log(mu_) : 0.;
double v = sqrt(1.-4.*mu2_),v2(sqr(v));
double omv = 4.*mu2_/(1.+v);
double f1,f2,fNS,VNS;
double r = omv/(1.+v);
double lr = mu_!=0. ? log(r) : 0.;
// normal form
if(mu_>1e-4) {
f1 = CF_*aS_/Constants::pi*
( +1. + 3.*log(0.5*(1.+v)) - 1.5*log(0.5*(1.+v2)) + sqr(Constants::pi)/6.
- 0.5*sqr(lr) - (1.+v2)/v*(lr*log(1.+v2) + sqr(Constants::pi)/12.
-0.5*log(4.*mu2_)*lr + 0.25*sqr(lr)));
fNS = -0.5*(1.+2.*v2)*lr/v + 1.5*lr - 2./3.*sqr(Constants::pi) + 0.5*sqr(lr)
+ (1.+v2)/v*(Herwig::Math::ReLi2(r) + sqr(Constants::pi)/3. - 0.25*sqr(lr)
+ lr*log((2.*v/ (1.+v))));
VNS = 1.5*log(0.5*(1.+v2))
+ 0.5*(1.+v2)/v*( 2.*lr*log(2.*(1.+v2)/sqr(1.+v))
+ 2.*Herwig::Math::ReLi2(sqr(r))
- 2.*Herwig::Math::ReLi2(2.*v/(1.+v)) - sqr(Constants::pi)/6.)
+ log(1.-mu_) - 2.*log(1.-2.*mu_) - 4.*mu2_/(1.+v2)*log(mu_/(1.-mu_))
- mu_/(1.-mu_)
+ 4.*(2.*mu2_-mu_)/(1.+v2) + 0.5*sqr(Constants::pi);
f2 = CF_*aS_/Constants::pi*mu2_*lr/v;
}
// small mass limit
else {
f1 = -CF_*aS_/Constants::pi/6.*
( - 6. - 24.*lmu*mu2_ - 15.*mu4 - 12.*mu4*lmu - 24.*mu4*sqr(lmu)
+ 2.*mu4*sqr(Constants::pi) - 12.*mu2_*mu4 - 96.*mu2_*mu4*sqr(lmu)
+ 8.*mu2_*mu4*sqr(Constants::pi) - 80.*mu2_*mu4*lmu);
fNS = - mu2_/18.*( + 36.*lmu - 36. - 45.*mu2_ + 216.*lmu*mu2_ - 24.*mu2_*sqr(Constants::pi)
+ 72.*mu2_*sqr(lmu) - 22.*mu4 + 1032.*mu4 * lmu
- 96.*mu4*sqr(Constants::pi) + 288.*mu4*sqr(lmu));
VNS = - mu2_/1260.*(-6930. + 7560.*lmu + 2520.*mu_ - 16695.*mu2_
+ 1260.*mu2_*sqr(Constants::pi)
+ 12600.*lmu*mu2_ + 1344.*mu_*mu2_ - 52780.*mu4 + 36960.*mu4*lmu
+ 5040.*mu4*sqr(Constants::pi) - 12216.*mu_*mu4);
f2 = CF_*aS_*mu2_/Constants::pi*( 2.*lmu + 4.*mu2_*lmu + 2.*mu2_ + 12.*mu4*lmu + 7.*mu4);
}
// add up bits for f1
f1 += CF_*aS_/Constants::pi*(fNS+VNS);
- // now for the real correction
- double phi = UseRandom::rnd()*Constants::twopi;
- // calculate y
- double yminus = 0.;
- double yplus = 1.-2.*mu_*(1.-mu_)/(1.-2*mu2_);
- double y = yminus + UseRandom::rnd()*(yplus-yminus);
- // calculate z
- double v1 = sqrt(sqr(2.*mu2_+(1.-2.*mu2_)*(1.-y))-4.*mu2_)/(1.-2.*mu2_)/(1.-y);
- double zplus = (1.+v1)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
- double zminus = (1.-v1)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
- double z = zminus + UseRandom::rnd()*(zplus-zminus);
- double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
- // calculate x1,x2
- double x2 = 1. - y*(1.-2.*mu2_);
- double x1 = 1. - z*(x2-2.*mu2_);
- // copy the particle objects over for calculateRealEmission
- vector<PPtr> hardProcess(3);
- hardProcess[0] = const_ptr_cast<PPtr>(&part);
- hardProcess[1] = decay[0];
- hardProcess[2] = decay[1];
- // total real emission contribution
- double realFact = 0.25*jac*sqr(1.-2.*mu2_)/
+ double realFact(0.);
+ for(int iemit=0;iemit<2;++iemit) {
+ // now for the real correction
+ double phi = UseRandom::rnd()*Constants::twopi;
+ // calculate y
+ double yminus = 0.;
+ double yplus = 1.-2.*mu_*(1.-mu_)/(1.-2*mu2_);
+ double y = yminus + UseRandom::rnd()*(yplus-yminus);
+ // calculate z
+ double v1 = sqrt(sqr(2.*mu2_+(1.-2.*mu2_)*(1.-y))-4.*mu2_)/(1.-2.*mu2_)/(1.-y);
+ double zplus = (1.+v1)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
+ double zminus = (1.-v1)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
+ double z = zminus + UseRandom::rnd()*(zplus-zminus);
+ double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
+ // calculate x1,x2
+ double x2 = 1. - y*(1.-2.*mu2_);
+ double x1 = 1. - z*(x2-2.*mu2_);
+ // copy the particle objects over for calculateRealEmission
+ vector<PPtr> hardProcess(3);
+ hardProcess[0] = const_ptr_cast<PPtr>(&part);
+ hardProcess[1] = decay[0];
+ hardProcess[2] = decay[1];
+ // total real emission contribution
+ realFact += 0.25*jac*sqr(1.-2.*mu2_)/
sqrt(1.-4.*mu2_)/Constants::twopi
- *2.*CF_*aS_*calculateRealEmission(x1, x2, hardProcess, phi, true);
+ *2.*CF_*aS_*calculateRealEmission(x1, x2, hardProcess, phi,
+ iemit, true);
+ }
// the born + virtual + real
output = output*(1. + f1 + realFact) + f2*total;
return output;
}
void SMZDecayer::doinitrun() {
PerturbativeDecayer::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<numberModes();++ix) {
if(ix<5) quarkWeight_ [ix ]=mode(ix)->maxWeight();
else if(ix<11) leptonWeight_[ix-5 ]=mode(ix)->maxWeight();
}
}
}
void SMZDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
for(unsigned int ix=0;ix<quarkWeight_.size();++ix) {
output << "newdef " << name() << ":QuarkMax " << ix << " "
<< quarkWeight_[ix] << "\n";
}
for(unsigned int ix=0;ix<leptonWeight_.size();++ix) {
output << "newdef " << name() << ":LeptonMax " << ix << " "
<< leptonWeight_[ix] << "\n";
}
// parameters for the PerturbativeDecayer base class
PerturbativeDecayer::dataBaseOutput(output,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
InvEnergy2 SMZDecayer::
realEmissionME(unsigned int,const Particle &parent,
ParticleVector &children,
unsigned int iemitter,
double ctheta, double stheta,
const LorentzRotation & rot1,
const LorentzRotation & rot2) {
// check the number of products and parent
assert(children.size()==3 && parent.id()==ParticleID::Z0);
// the electric charge
double e = sqrt(SM().alphaEM()*4.*Constants::pi);
// azimuth of the photon
double phi = children[2]->momentum().phi();
// wavefunctions for the decaying particle in the rotated dipole frame
vector<VectorWaveFunction> vec1 = _vectors;
for(unsigned int ix=0;ix<vec1.size();++ix) {
vec1[ix].transform(rot1);
vec1[ix].transform(rot2);
}
// wavefunctions for the decaying particle in the rotated rest frame
vector<VectorWaveFunction> vec2 = _vectors;
for(unsigned int ix=0;ix<vec1.size();++ix) {
vec2[ix].transform(rot2);
}
// find the outgoing particle and antiparticle
unsigned int iferm(0),ianti(1);
if(children[iferm]->id()<0) swap(iferm,ianti);
// wavefunctions for the particles before the radiation
// wavefunctions for the outgoing fermion
SpinorBarWaveFunction wavebartemp;
Lorentz5Momentum ptemp = - _wavebar[0].momentum();
ptemp *= rot2;
if(ptemp.perp()/ptemp.e()<1e-10) {
ptemp.setX(ZERO);
ptemp.setY(ZERO);
}
wavebartemp = SpinorBarWaveFunction(ptemp,_wavebar[0].particle(),outgoing);
// wavefunctions for the outgoing antifermion
SpinorWaveFunction wavetemp;
ptemp = - _wave[0].momentum();
ptemp *= rot2;
if(ptemp.perp()/ptemp.e()<1e-10) {
ptemp.setX(ZERO);
ptemp.setY(ZERO);
}
wavetemp = SpinorWaveFunction(ptemp,_wave[0].particle(),outgoing);
// loop over helicities
vector<SpinorWaveFunction> wf_old;
vector<SpinorBarWaveFunction> wfb_old;
for(unsigned int ihel=0;ihel<2;++ihel) {
wavetemp.reset(ihel);
wf_old.push_back(wavetemp);
wavebartemp.reset(ihel);
wfb_old.push_back(wavebartemp);
}
// calculate the wave functions for the new fermions
// ensure the momenta have pT=0
for(unsigned int ix=0;ix<2;++ix) {
Lorentz5Momentum ptemp = children[ix]->momentum();
if(ptemp.perp()/ptemp.e()<1e-10) {
ptemp.setX(ZERO);
ptemp.setY(ZERO);
children[ix]->set5Momentum(ptemp);
}
}
// calculate the wavefunctions
vector<SpinorBarWaveFunction> wfb;
SpinorBarWaveFunction::calculateWaveFunctions(wfb,children[iferm],outgoing);
vector<SpinorWaveFunction> wf;
SpinorWaveFunction::calculateWaveFunctions (wf ,children[ianti],outgoing);
// wave functions for the photons
vector<VectorWaveFunction> photon;
VectorWaveFunction::calculateWaveFunctions(photon,children[2],outgoing,true);
// loop to calculate the matrix elements
Complex lome[3][2][2],diffme[3][2][2][2],summe[3][2][2][2];
Energy2 scale(sqr(parent.mass()));
Complex diff[2]={0.,0.};
Complex sum [2]={0.,0.};
for(unsigned int ifm=0;ifm<2;++ifm) {
for(unsigned int ia=0;ia<2;++ia) {
for(unsigned int vhel=0;vhel<3;++vhel) {
// calculation of the leading-order matrix element
Complex loamp = FFZVertex_->evaluate(scale,wf_old[ia],
wfb_old[ifm],vec2[vhel]);
Complex lotemp = FFZVertex_->evaluate(scale,wf[ia],
wfb[ifm],vec1[vhel]);
if(iferm>ianti) lome[vhel][ia][ifm] = loamp;
else lome[vhel][ifm][ia] = loamp;
// photon loop for the real emmision terms
for(unsigned int phel=0;phel<2;++phel) {
// radiation from the antifermion
// normal case with small angle treatment
if(children[2 ]->momentum().z()/
children[iferm]->momentum().z()>=ZERO && iemitter == iferm ) {
Complex dipole = e*double(children[iferm]->dataPtr()->iCharge())/3.*
UnitRemoval::E*loamp*
(children[iferm]->momentum()*photon[2*phel].wave())/
(children[iferm]->momentum()*children[2]->momentum());
// sum and difference
SpinorBarWaveFunction foff =
FFPVertex_->evaluateSmall(ZERO,3,children[iferm]->dataPtr()->CC(),
wfb[ifm],photon[2*phel],
ifm,2*phel,ctheta,phi,stheta,false);
diff[0] = FFZVertex_->evaluate(scale,wf[ia],foff,vec1[vhel]) +
e*double(children[iferm]->dataPtr()->iCharge())/3.*
UnitRemoval::E*(lotemp-loamp)*
(children[iferm]->momentum()*photon[2*phel].wave())/
(children[iferm]->momentum()*children[2]->momentum());
sum [0] = diff[0]+2.*dipole;
}
// special if fermion backwards
else {
SpinorBarWaveFunction foff =
FFPVertex_->evaluate(ZERO,3,children[iferm]->dataPtr()->CC(),
wfb[ifm],photon[2*phel]);
Complex diag =
FFZVertex_->evaluate(scale,wf[ia],foff,vec1[vhel]);
Complex dipole = e*double(children[iferm]->dataPtr()->iCharge())/3.*
UnitRemoval::E*loamp*
(children[iferm]->momentum()*photon[2*phel].wave())/
(children[iferm]->momentum()*children[2]->momentum());
diff[0] = diag-dipole;
sum [0] = diag+dipole;
}
// radiation from the anti fermion
// small angle case in general
if(children[2 ]->momentum().z()/
children[ianti]->momentum().z()>=ZERO && iemitter == ianti ) {
Complex dipole = e*double(children[ianti]->dataPtr()->iCharge())/3.*
UnitRemoval::E*loamp*
(children[ianti]->momentum()*photon[2*phel].wave())/
(children[ianti]->momentum()*children[2]->momentum());
// sum and difference
SpinorWaveFunction foff =
FFPVertex_->evaluateSmall(ZERO,3,children[ianti]->dataPtr()->CC(),
wf[ia],photon[2*phel],
ia,2*phel,ctheta,phi,stheta,false);
diff[1] = FFZVertex_->evaluate(scale,foff ,wfb[ifm],vec1[vhel]) +
e*double(children[ianti]->dataPtr()->iCharge())/3.*
UnitRemoval::E*(lotemp-loamp)*
(children[ianti]->momentum()*photon[2*phel].wave())/
(children[ianti]->momentum()*children[2]->momentum());
sum [1] = diff[1]+2.*dipole;
}
// special if fermion backwards after radiation
else {
SpinorWaveFunction foff =
FFPVertex_->evaluate(ZERO,3,children[ianti]->dataPtr()->CC(),
wf[ia],photon[2*phel]);
Complex diag =
FFZVertex_->evaluate(scale,foff ,wfb[ifm],vec1[vhel]);
Complex dipole = e*double(children[ianti]->dataPtr()->iCharge())/3.*
UnitRemoval::E*loamp*
(children[ianti]->momentum()*photon[2*phel].wave())/
(children[ianti]->momentum()*children[2]->momentum());
// sum and difference
diff[1] = diag - dipole;
sum [1] = diag + dipole;
}
// add to me
if(iferm>ianti) {
diffme[vhel][ia][ifm][phel] = diff[0] + diff[1];
summe [vhel][ia][ifm][phel] = sum[0] + sum[1] ;
}
else {
diffme [vhel][ifm][ia][phel] = diff[0] + diff[1];
summe [vhel][ifm][ia][phel] = sum[0] + sum[1] ;
}
}
}
}
}
// cerr << parent << "\n";
// for(unsigned int ix=0;ix<children.size();++ix) {
// cerr << *children[ix] << "\n";
// }
// _rho = RhoDMatrix(PDT::Spin1);
Complex lo(0.),difference(0.);
for(unsigned int vhel1=0;vhel1<3;++vhel1) {
for(unsigned int vhel2=0;vhel2<3;++vhel2) {
for(unsigned int ifm=0;ifm<2;++ifm) {
for(unsigned int ia=0;ia<2;++ia) {
lo += _rho(vhel1,vhel2)*lome[vhel1][ifm][ia]*conj(lome[vhel2][ifm][ia]);
for(unsigned int phel=0;phel<2;++phel) {
difference +=
_rho(vhel1,vhel2)*diffme[vhel1][ifm][ia][phel]*conj(summe[vhel2][ifm][ia][phel]);
}
}
}
}
}
// // analytic result
// double iCharge = children[0]->dataPtr()->iCharge()*
// children[1]->dataPtr()->iCharge()/9.;
// Energy2 ubar = 2.*children[0]->momentum()*children[2]->momentum();
// Energy2 tbar = 2.*children[1]->momentum()*children[2]->momentum();
// double mu2 = sqr(children[1]->mass()/parent.mass());
// double gL = (FFZVertex_->left() *FFZVertex_->norm()).real();
// double gR = (FFZVertex_->right()*FFZVertex_->norm()).real();
// Energy2 den = sqr(parent.mass())*(((sqr(gL)+sqr(gR))*(1-mu2)+6.*mu2*gL*gR));
// InvEnergy2 anal = -iCharge*( 2.*(ubar/tbar+tbar/ubar)/sqr(parent.mass())+
// 4.*mu2/den*((sqr(gL)+sqr(gR))*(1+ubar/tbar+tbar/ubar)
// -2.*gL*gR*(1.+2.*(ubar/tbar+tbar/ubar))));
// cerr << "testing ratio " << parent.PDGName()
// << " " << difference.real()/sqr(e)/lo.real()*UnitRemoval::InvE2/(anal) << "\n"
// << stheta << " " << ctheta << "\n";
return difference.real()/sqr(e)/lo.real()*UnitRemoval::InvE2;
}
double SMZDecayer::oneLoopVirtualME(unsigned int,
const Particle & parent,
const ParticleVector & children) {
assert(children.size()==2);
// velocities of the particles
double beta = sqrt(1.-4.*sqr(children[0]->mass()/parent.mass()));
double opb = 1.+beta;
double omb = 4.*sqr(children[0]->mass()/parent.mass())/opb;
// couplings
double gL = (FFZVertex_->left() *FFZVertex_->norm()).real();
double gR = (FFZVertex_->right()*FFZVertex_->norm()).real();
double gA = 0.5*(gL-gR);
double gV = 0.5*(gL+gR);
// correction terms
double ln = log(omb/opb);
double f1 = 1. + ln*beta;
double fA = 1. + ln/beta;
InvEnergy f2 = 0.5*sqrt(omb*opb)/parent.mass()/beta*ln;
// momentum difference for the loop
Lorentz5Momentum q = children[0]->momentum()-children[1]->momentum();
if(children[0]->id()<0) q *= -1.;
// spinors
vector<LorentzSpinor <SqrtEnergy> > sp;
vector<LorentzSpinorBar<SqrtEnergy> > sbar;
for(unsigned int ix=0;ix<2;++ix) {
sp .push_back( _wave[ix].dimensionedWave());
sbar.push_back(_wavebar[ix].dimensionedWave());
}
// polarization vectors
vector<LorentzPolarizationVector> pol;
for(unsigned int ix=0;ix<3;++ix)
pol.push_back(_vectors[ix].wave());
// matrix elements
complex<Energy> lome[3][2][2],loopme[3][2][2];
for(unsigned int vhel=0;vhel<3;++vhel) {
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
complex<Energy> vector =
sp[ihel1].generalCurrent(sbar[ihel2], 1.,1.).dot(pol[vhel]);
complex<Energy> axial =
sp[ihel1].generalCurrent(sbar[ihel2],-1.,1.).dot(pol[vhel]);
complex<Energy2> scalar =
sp[ihel1].scalar(sbar[ihel2])*(q*pol[vhel]);
lome [vhel][ihel1][ihel2] = gV* vector-gA* axial;
loopme[vhel][ihel1][ihel2] = gV*f1*vector-gA*fA*axial+scalar*f2*gV;
}
}
}
// sum sums
complex<Energy2> den(ZERO),num(ZERO);
for(unsigned int vhel1=0;vhel1<3;++vhel1) {
for(unsigned int vhel2=0;vhel2<3;++vhel2) {
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
num += _rho(vhel1,vhel2)*
( lome[vhel1][ihel1][ihel2]*conj(loopme[vhel2][ihel1][ihel2])+
loopme[vhel1][ihel1][ihel2]*conj( lome[vhel2][ihel1][ihel2]));
den += _rho(vhel1,vhel2)*
lome[vhel1][ihel1][ihel2]*conj(lome[vhel2][ihel1][ihel2]);
}
}
}
}
// prefactor
double iCharge = children[0]->dataPtr()->iCharge()*
children[1]->dataPtr()->iCharge()/9.;
double pre = 0.5*SM().alphaEM()*iCharge/Constants::pi;
// output
return pre*num.real()/den.real();
}
void SMZDecayer::
initializeMECorrection(RealEmissionProcessPtr born, double & initial,
double & final) {
// get the quark and antiquark
ParticleVector qq;
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix)
qq.push_back(born->bornOutgoing()[ix]);
// ensure quark first
if(qq[0]->id()<0) swap(qq[0],qq[1]);
// centre of mass energy
d_Q_ = (qq[0]->momentum() + qq[1]->momentum()).m();
// quark mass
d_m_ = 0.5*(qq[0]->momentum().m()+qq[1]->momentum().m());
// set the other parameters
setRho(sqr(d_m_/d_Q_));
setKtildeSymm();
// otherwise can do it
initial=1.;
final =1.;
}
RealEmissionProcessPtr SMZDecayer::
applyHardMatrixElementCorrection(RealEmissionProcessPtr born) {
// get the quark and antiquark
ParticleVector qq;
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix)
qq.push_back(born->bornOutgoing()[ix]);
if(!qq[0]->dataPtr()->coloured()) return RealEmissionProcessPtr();
// ensure quark first
bool order = qq[0]->id()<0;
if(order) swap(qq[0],qq[1]);
// get the momenta
vector<Lorentz5Momentum> newfs = applyHard(qq);
// return if no emission
if(newfs.size()!=3) return RealEmissionProcessPtr();
// perform final check to ensure energy greater than constituent mass
for (int i=0; i<2; i++) {
if (newfs[i].e() < qq[i]->data().constituentMass()) return RealEmissionProcessPtr();
}
if (newfs[2].e() < getParticleData(ParticleID::g)->constituentMass())
return RealEmissionProcessPtr();
// set masses
for (int i=0; i<2; i++) newfs[i].setMass(qq[i]->mass());
newfs[2].setMass(ZERO);
// decide which particle emits
bool firstEmits=
newfs[2].vect().perp2(newfs[0].vect())<
newfs[2].vect().perp2(newfs[1].vect());
// create the new quark, antiquark and gluon
PPtr newg = getParticleData(ParticleID::g)->produceParticle(newfs[2]);
PPtr newq = qq[0]->dataPtr()->produceParticle(newfs[0]);
PPtr newa = qq[1]->dataPtr()->produceParticle(newfs[1]);
// create the output real emission process
for(unsigned int ix=0;ix<born->bornIncoming().size();++ix) {
born->incoming().push_back(born->bornIncoming()[ix]);
}
if(!order) {
born->outgoing().push_back(newq);
born->outgoing().push_back(newa);
born->outgoing().push_back(newg);
}
else {
born->outgoing().push_back(newa);
born->outgoing().push_back(newq);
born->outgoing().push_back(newg);
firstEmits = !firstEmits;
}
// make colour connections
newg->colourNeighbour(newq);
newa->colourNeighbour(newg);
if(firstEmits) {
born->emitter(1);
born->spectator(2);
}
else {
born->emitter(2);
born->spectator(1);
}
born->emitted(3);
born->interaction(ShowerInteraction::QCD);
return born;
}
vector<Lorentz5Momentum> SMZDecayer::
applyHard(const ParticleVector &p) {
double x, xbar;
vector<Lorentz5Momentum> fs;
// return if no emission
if (getHard(x, xbar) < UseRandom::rnd() || p.size() != 2) return fs;
// centre of mass energy
Lorentz5Momentum pcm = p[0]->momentum() + p[1]->momentum();
// momenta of quark,antiquark and gluon
Lorentz5Momentum pq, pa, pg;
if (p[0]->id() > 0) {
pq = p[0]->momentum();
pa = p[1]->momentum();
} else {
pa = p[0]->momentum();
pq = p[1]->momentum();
}
// boost to boson rest frame
Boost beta = (pcm.findBoostToCM());
pq.boost(beta);
pa.boost(beta);
// return if fails ?????
double xg = 2.-x-xbar;
if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return fs;
Axis u1, u2, u3;
// moduli of momenta in units of Q and cos theta
// stick to q direction?
// p1 is the one that is kept, p2 is the other fermion, p3 the gluon.
Energy e1, e2, e3;
Energy pp1, pp2, pp3;
bool keepq = true;
if (UseRandom::rnd() > sqr(x)/(sqr(x)+sqr(xbar)))
keepq = false;
if (keepq) {
pp1 = d_Q_*sqrt(sqr(x)-4.*d_rho_)/2.;
pp2 = d_Q_*sqrt(sqr(xbar)-4.*d_rho_)/2.;
e1 = d_Q_*x/2.;
e2 = d_Q_*xbar/2.;
u1 = pq.vect().unit();
} else {
pp2 = d_Q_*sqrt(sqr(x)-4.*d_rho_)/2.;
pp1 = d_Q_*sqrt(sqr(xbar)-4.*d_rho_)/2.;
e2 = d_Q_*x/2.;
e1 = d_Q_*xbar/2.;
u1 = pa.vect().unit();
}
pp3 = d_Q_*xg/2.;
e3 = pp3;
u2 = u1.orthogonal();
u2 /= u2.mag();
u3 = u1.cross(u2);
u3 /= u3.mag();
double ct2=-2., ct3=-2.;
if (pp1 == ZERO || pp2 == ZERO || pp3 == ZERO) {
bool touched = false;
if (pp1 == ZERO) {
ct2 = 1;
ct3 = -1;
touched = true;
}
if (pp2 == ZERO || pp3 == ZERO) {
ct2 = 1;
ct3 = 1;
touched = true;
}
if (!touched)
throw Exception() << "SMZDecayer::applyHard()"
<< " did not set ct2/3"
<< Exception::abortnow;
} else {
ct3 = (sqr(pp1)+sqr(pp3)-sqr(pp2))/(2.*pp1*pp3);
ct2 = (sqr(pp1)+sqr(pp2)-sqr(pp3))/(2.*pp1*pp2);
}
double phi = Constants::twopi*UseRandom::rnd();
double cphi = cos(phi);
double sphi = sin(phi);
double st2 = sqrt(1.-sqr(ct2));
double st3 = sqrt(1.-sqr(ct3));
ThreeVector<Energy> pv1, pv2, pv3;
pv1 = pp1*u1;
pv2 = -ct2*pp2*u1 + st2*cphi*pp2*u2 + st2*sphi*pp2*u3;
pv3 = -ct3*pp3*u1 - st3*cphi*pp3*u2 - st3*sphi*pp3*u3;
if (keepq) {
pq = Lorentz5Momentum(pv1, e1);
pa = Lorentz5Momentum(pv2, e2);
} else {
pa = Lorentz5Momentum(pv1, e1);
pq = Lorentz5Momentum(pv2, e2);
}
pg = Lorentz5Momentum(pv3, e3);
pq.boost(-beta);
pa.boost(-beta);
pg.boost(-beta);
fs.push_back(pq);
fs.push_back(pa);
fs.push_back(pg);
return fs;
}
double SMZDecayer::getHard(double &x1, double &x2) {
double w = 0.0;
double y1 = UseRandom::rnd(),y2 = UseRandom::rnd();
// simply double MC efficiency
// -> weight has to be divided by two (Jacobian)
if (y1 + y2 > 1) {
y1 = 1.-y1;
y2 = 1.-y2;
}
bool inSoft = false;
if (y1 < 0.25) {
if (y2 < 0.25) {
inSoft = true;
if (y1 < y2) {
y1 = 0.25-y1;
y2 = y1*(1.5 - 2.*y2);
} else {
y2 = 0.25 - y2;
y1 = y2*(1.5 - 2.*y1);
}
} else {
if (y2 < y1 + 2.*sqr(y1)) return w;
}
} else {
if (y2 < 0.25) {
if (y1 < y2 + 2.*sqr(y2)) return w;
}
}
// inside PS?
x1 = 1.-y1;
x2 = 1.-y2;
if(y1*y2*(1.-y1-y2) < d_rho_*sqr(y1+y2)) return w;
double k1 = getKfromX(x1, x2);
double k2 = getKfromX(x2, x1);
// Is it in the quark emission zone?
if (k1 < d_kt1_) return 0.0;
// No...is it in the anti-quark emission zone?
if (k2 < d_kt2_) return 0.0;
// Point is in dead zone: compute q qbar g weight
w = MEV(x1, x2);
// for axial:
// w = MEA(x1, x2);
// Reweight soft region
if (inSoft) {
if (y1 < y2) w *= 2.*y1;
else w *= 2.*y2;
}
// alpha and colour factors
Energy2 pt2 = sqr(d_Q_)*(1.-x1)*(1.-x2);
- w *= 1./3./Constants::pi*alpha_->value(pt2);
+ w *= 1./3./Constants::pi*coupling()->value(pt2);
return w;
}
bool SMZDecayer::
softMatrixElementVeto(ShowerProgenitorPtr initial,ShowerParticlePtr parent,Branching br) {
// check we should be applying the veto
if(parent->id()!=initial->progenitor()->id()||
br.ids[0]->id()!=br.ids[1]->id()||
br.ids[2]->id()!=ParticleID::g) return false;
// calculate pt
double d_z = br.kinematics->z();
Energy d_qt = br.kinematics->scale();
Energy2 d_m2 = parent->momentum().m2();
Energy pPerp = (1.-d_z)*sqrt( sqr(d_z*d_qt) - d_m2);
// if not hardest so far don't apply veto
if(pPerp<initial->highestpT()) return false;
// calculate the weight
double weight = 0.;
if(parent->id()>0) weight = qWeightX(d_qt, d_z);
else weight = qbarWeightX(d_qt, d_z);
// compute veto from weight
bool veto = !UseRandom::rndbool(weight);
// if vetoing reset the scale
if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
// return the veto
return veto;
}
void SMZDecayer::setRho(double r)
{
d_rho_ = r;
d_v_ = sqrt(1.-4.*d_rho_);
}
void SMZDecayer::setKtildeSymm() {
d_kt1_ = (1. + sqrt(1. - 4.*d_rho_))/2.;
setKtilde2();
}
void SMZDecayer::setKtilde2() {
double num = d_rho_ * d_kt1_ + 0.25 * d_v_ *(1.+d_v_)*(1.+d_v_);
double den = d_kt1_ - d_rho_;
d_kt2_ = num/den;
}
double SMZDecayer::getZfromX(double x1, double x2) {
double uval = u(x2);
double num = x1 - (2. - x2)*uval;
double den = sqrt(x2*x2 - 4.*d_rho_);
return uval + num/den;
}
double SMZDecayer::getKfromX(double x1, double x2) {
double zval = getZfromX(x1, x2);
return (1.-x2)/(zval*(1.-zval));
}
double SMZDecayer::MEV(double x1, double x2) {
// Vector part
double num = (x1+2.*d_rho_)*(x1+2.*d_rho_) + (x2+2.*d_rho_)*(x2+2.*d_rho_)
- 8.*d_rho_*(1.+2.*d_rho_);
double den = (1.+2.*d_rho_)*(1.-x1)*(1.-x2);
return (num/den - 2.*d_rho_/((1.-x1)*(1.-x1))
- 2*d_rho_/((1.-x2)*(1.-x2)))/d_v_;
}
double SMZDecayer::MEA(double x1, double x2) {
// Axial part
double num = (x1+2.*d_rho_)*(x1+2.*d_rho_) + (x2+2.*d_rho_)*(x2+2.*d_rho_)
+ 2.*d_rho_*((5.-x1-x2)*(5.-x1-x2) - 19.0 + 4*d_rho_);
double den = d_v_*d_v_*(1.-x1)*(1.-x2);
return (num/den - 2.*d_rho_/((1.-x1)*(1.-x1))
- 2*d_rho_/((1.-x2)*(1.-x2)))/d_v_;
}
double SMZDecayer::u(double x2) {
return 0.5*(1. + d_rho_/(1.-x2+d_rho_));
}
void SMZDecayer::
getXXbar(double kti, double z, double &x, double &xbar) {
double w = sqr(d_v_) + kti*(-1. + z)*z*(2. + kti*(-1. + z)*z);
if (w < 0) {
x = -1.;
xbar = -1;
} else {
x = (1. + sqr(d_v_)*(-1. + z) + sqr(kti*(-1. + z))*z*z*z
+ z*sqrt(w)
- kti*(-1. + z)*z*(2. + z*(-2 + sqrt(w))))/
(1. - kti*(-1. + z)*z + sqrt(w));
xbar = 1. + kti*(-1. + z)*z;
}
}
double SMZDecayer::qWeight(double x, double xbar) {
double rval;
double xg = 2. - xbar - x;
// always return one in the soft gluon region
if(xg < EPS_) return 1.0;
// check it is in the phase space
if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return 0.0;
double k1 = getKfromX(x, xbar);
double k2 = getKfromX(xbar, x);
// Is it in the quark emission zone?
if(k1 < d_kt1_) {
rval = MEV(x, xbar)/PS(x, xbar);
// is it also in the anti-quark emission zone?
if(k2 < d_kt2_) rval *= 0.5;
return rval;
}
return 1.0;
}
double SMZDecayer::qbarWeight(double x, double xbar) {
double rval;
double xg = 2. - xbar - x;
// always return one in the soft gluon region
if(xg < EPS_) return 1.0;
// check it is in the phase space
if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return 0.0;
double k1 = getKfromX(x, xbar);
double k2 = getKfromX(xbar, x);
// Is it in the antiquark emission zone?
if(k2 < d_kt2_) {
rval = MEV(x, xbar)/PS(xbar, x);
// is it also in the quark emission zone?
if(k1 < d_kt1_) rval *= 0.5;
return rval;
}
return 1.0;
}
double SMZDecayer::qWeightX(Energy qtilde, double z) {
double x, xb;
getXXbar(sqr(qtilde/d_Q_), z, x, xb);
// if exceptionally out of phase space, leave this emission, as there
// is no good interpretation for the soft ME correction.
if (x < 0 || xb < 0) return 1.0;
return qWeight(x, xb);
}
double SMZDecayer::qbarWeightX(Energy qtilde, double z) {
double x, xb;
getXXbar(sqr(qtilde/d_Q_), z, xb, x);
// see above in qWeightX.
if (x < 0 || xb < 0) return 1.0;
return qbarWeight(x, xb);
}
double SMZDecayer::PS(double x, double xbar) {
double u = 0.5*(1. + d_rho_ / (1.-xbar+d_rho_));
double z = u + (x - (2.-xbar)*u)/sqrt(xbar*xbar - 4.*d_rho_);
double brack = (1.+z*z)/(1.-z)- 2.*d_rho_/(1-xbar);
// interesting: the splitting function without the subtraction
// term. Actually gives a much worse approximation in the collinear
// limit. double brack = (1.+z*z)/(1.-z);
double den = (1.-xbar)*sqrt(xbar*xbar - 4.*d_rho_);
return brack/den;
}
-
-RealEmissionProcessPtr SMZDecayer::
-generateHardest(RealEmissionProcessPtr born) {
- assert(born->bornOutgoing().size()==2);
- // check coloured
- if(!born->bornOutgoing()[0]->dataPtr()->coloured()) return RealEmissionProcessPtr();
- // extract required info
- partons_.resize(2);
- quark_.resize(2);
- vector<PPtr> hardProcess;
- zboson_ = born->bornIncoming()[0];
- hardProcess.push_back(zboson_);
- for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
- partons_[ix] = born->bornOutgoing()[ix]->dataPtr();
- quark_[ix] = born->bornOutgoing()[ix]->momentum();
- quark_[ix].setMass(partons_[ix]->mass());
- hardProcess.push_back(born->bornOutgoing()[ix]);
+double SMZDecayer::matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
+ const ParticleVector & decay3, MEOption) {
+ // extract partons and LO momentas
+ vector<cPDPtr> partons(1,inpart.dataPtr());
+ vector<Lorentz5Momentum> lomom(1,inpart.momentum());
+ for(unsigned int ix=0;ix<2;++ix) {
+ partons.push_back(decay2[ix]->dataPtr());
+ lomom.push_back(decay2[ix]->momentum());
}
- bool order = partons_[0]->id()<0;
- if(order) {
- swap(partons_[0] ,partons_[1] );
- swap(quark_[0] ,quark_[1] );
- swap(hardProcess[1],hardProcess[2]);
+ vector<Lorentz5Momentum> realmom(1,inpart.momentum());
+ for(unsigned int ix=0;ix<3;++ix) {
+ if(ix==2) partons.push_back(decay3[ix]->dataPtr());
+ realmom.push_back(decay3[ix]->momentum());
}
- gauge_.setMass(0.*MeV);
- // Get the Z boson mass.
- mz2_ = (quark_[0] + quark_[1]).m2();
- // Generate emission and set _quark[0,1] and _gauge to be the
- // momenta of q, qbar and g after the hardest emission:
- if(!getEvent(hardProcess)) {
- born->pT()[ShowerInteraction::QCD] = pTmin_;
- return born;
+ if(partons[0]->id()<0) {
+ swap(partons[1],partons[2]);
+ swap(lomom[1],lomom[2]);
+ swap(realmom[1],realmom[2]);
}
- // Ensure the energies are greater than the constituent masses:
- for (int i=0; i<2; i++) {
- if (quark_[i].e() < partons_[i]->constituentMass()) return RealEmissionProcessPtr();
- }
- if (gauge_.e() < gluon_ ->constituentMass()) return RealEmissionProcessPtr();
- // set masses
- quark_[0].setMass( partons_[0]->mass() );
- quark_[1].setMass( partons_[1]->mass() );
- gauge_ .setMass( ZERO );
- // assign the emitter based on evolution scales
- unsigned int iemitter = quark_[0]*gauge_ > quark_[1]*gauge_ ? 2 : 1;
- unsigned int ispectator = iemitter==1 ? 1 : 2;
- // create new partices and insert
- PPtr zboson = zboson_->dataPtr()->produceParticle(zboson_->momentum());
- born->incoming().push_back(zboson);
- PPtr newq = partons_[0]->produceParticle(quark_[0]);
- PPtr newa = partons_[1]->produceParticle(quark_[1]);
- PPtr newg = gluon_->produceParticle(gauge_);
- // make colour connections
- newg->colourNeighbour(newq);
- newa->colourNeighbour(newg);
- // insert in output structure
- if(!order) {
- born->outgoing().push_back(newq);
- born->outgoing().push_back(newa);
- }
- else {
- born->outgoing().push_back(newa);
- born->outgoing().push_back(newq);
- swap(iemitter,ispectator);
- }
- born->outgoing().push_back(newg);
- born->emitter (iemitter );
- born->spectator(ispectator);
- born->emitted (3);
- born->pT()[ShowerInteraction::QCD] = pT_;
- // return process
- born->interaction(ShowerInteraction::QCD);
- return born;
+ double lome = loME(partons,lomom);
+ InvEnergy2 reme = realME(partons,realmom);
+ double ratio = CF_*reme/lome*sqr(inpart.mass());
+ return ratio;
}
double SMZDecayer::meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter, bool subtract) const {
Lorentz5Momentum q = momenta[1]+momenta[2]+momenta[3];
Energy2 Q2=q.m2();
Energy2 lambda = sqrt((Q2-sqr(momenta[1].mass()+momenta[2].mass()))*
(Q2-sqr(momenta[1].mass()-momenta[2].mass())));
InvEnergy2 D[2];
double lome[2];
for(unsigned int iemit=0;iemit<2;++iemit) {
unsigned int ispect = iemit==0 ? 1 : 0;
Energy2 pipj = momenta[3 ] * momenta[1+iemit ];
Energy2 pipk = momenta[3 ] * momenta[1+ispect];
Energy2 pjpk = momenta[1+iemit] * momenta[1+ispect];
double y = pipj/(pipj+pipk+pjpk);
double z = pipk/( pipk+pjpk);
Energy mij = sqrt(2.*pipj+sqr(momenta[1+iemit].mass()));
Energy2 lamB = sqrt((Q2-sqr(mij+momenta[1+ispect].mass()))*
(Q2-sqr(mij-momenta[1+ispect].mass())));
Energy2 Qpk = q*momenta[1+ispect];
Lorentz5Momentum pkt =
lambda/lamB*(momenta[1+ispect]-Qpk/Q2*q)
+0.5/Q2*(Q2+sqr(momenta[1+ispect].mass())-sqr(momenta[1+ispect].mass()))*q;
Lorentz5Momentum pijt =
q-pkt;
double muj = momenta[1+iemit ].mass()/sqrt(Q2);
double muk = momenta[1+ispect].mass()/sqrt(Q2);
double vt = sqrt((1.-sqr(muj+muk))*(1.-sqr(muj-muk)))/(1.-sqr(muj)-sqr(muk));
double v = sqrt(sqr(2.*sqr(muk)+(1.-sqr(muj)-sqr(muk))*(1.-y))-4.*sqr(muk))
/(1.-y)/(1.-sqr(muj)-sqr(muk));
// dipole term
D[iemit] = 0.5/pipj*(2./(1.-(1.-z)*(1.-y))
-vt/v*(2.-z+sqr(momenta[1+iemit].mass())/pipj));
// matrix element
vector<Lorentz5Momentum> lomom(3);
lomom[0] = momenta[0];
if(iemit==0) {
lomom[1] = pijt;
lomom[2] = pkt ;
}
else {
lomom[2] = pijt;
lomom[1] = pkt ;
}
lome[iemit] = loME(partons,lomom);
}
InvEnergy2 ratio = realME(partons,momenta)*abs(D[iemitter])
/(abs(D[0]*lome[0])+abs(D[1]*lome[1]));
if(subtract)
return Q2*(ratio-2.*D[iemitter]);
else
return Q2*ratio;
}
double SMZDecayer::loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
VectorWaveFunction zin (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
}
for(unsigned int ix=0;ix<3;++ix){
zin.reset(ix);
vin.push_back(zin);
}
// temporary storage of the different diagrams
// sum over helicities to get the matrix element
double total(0.);
for(unsigned int inhel=0;inhel<3;++inhel) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
Complex diag1 = FFZVertex_->evaluate(scale_,aout[outhel2],fout[outhel1],vin[inhel]);
total += norm(diag1);
}
}
}
// return the answer
return total;
}
InvEnergy2 SMZDecayer::realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
vector<VectorWaveFunction> gout;
VectorWaveFunction zin (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
VectorWaveFunction gluon(momenta[3],partons[3],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
gluon.reset(2*ix);
gout.push_back(gluon);
}
for(unsigned int ix=0;ix<3;++ix){
zin.reset(ix);
vin.push_back(zin);
}
vector<Complex> diag(2,0.);
double total(0.);
for(unsigned int inhel1=0;inhel1<3;++inhel1) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
for(unsigned int outhel3=0;outhel3<2;++outhel3) {
SpinorBarWaveFunction off1 =
FFGVertex_->evaluate(scale_,3,partons[1],fout[outhel1],gout[outhel3]);
diag[0] = FFZVertex_->evaluate(scale_,aout[outhel2],off1,vin[inhel1]);
SpinorWaveFunction off2 =
FFGVertex_->evaluate(scale_,3,partons[2],aout[outhel2],gout[outhel3]);
diag[1] = FFZVertex_->evaluate(scale_,off2,fout[outhel1],vin[inhel1]);
// sum of diagrams
Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.));
// me2
total += norm(sum);
}
}
}
}
// divide out the coupling
total /= norm(FFGVertex_->norm());
// return the total
return total*UnitRemoval::InvE2;
}
-bool SMZDecayer::getEvent(vector<PPtr> hardProcess) {
- Energy particleMass = ZERO;
- for(unsigned int ix=0;ix<hardProcess.size();++ix) {
- if(hardProcess[ix]->id()==ParticleID::Z0) {
- mZ_ = hardProcess[ix]->mass();
- }
- else {
- particleMass = hardProcess[ix]->mass();
- }
- }
- // reduced mass
- mu_ = particleMass/mZ_;
- mu2_ = sqr(mu_);
- // scale
- scale_ = sqr(mZ_);
- // max pT
- Energy pTmax = 0.5*sqrt(mz2_);
- if(pTmax<pTmin_) return false;
- // Define over valued y_max & y_min according to the associated pt_min cut.
- double ymax = acosh(pTmax/pTmin_);
- double ymin = -ymax;
- // pt of the emmission
- pT_ = pTmax;
- // prefactor
- double overEst = 4.;
- double prefactor = overEst*alphaS()->overestimateValue()*CF_*
- (ymax-ymin)/Constants::twopi;
- // loop to generate the pt and rapidity
- bool reject;
- //arrays to hold the temporary probabilities whilst the for loop progresses
- double probTemp[2][2]={{0.,0.},{0.,0.}};
- probTemp[0][0]=probTemp[0][1]=probTemp[1][0]=probTemp[1][1]=0.;
- double x1Solution[2][2] = {{0.,0.},{0.,0.}};
- double x2Solution[2][2] = {{0.,0.},{0.,0.}};
- double x3Solution[2] = {0.,0.};
- Energy pT[2] = {pTmax,pTmax};
- double yTemp[2] = {0.,0.};
- double phi = 0.;
- // do the competition
- for(int i=0; i<2; i++) {
- do {
- //generation of phi
- phi = UseRandom::rnd() * Constants::twopi;
- // reject the emission
- reject = true;
- // generate pt
- pT[i] *= pow(UseRandom::rnd(),1./prefactor);
- Energy2 pT2 = sqr(pT[i]);
- if(pT[i]<pTmin_) {
- pT[i] = -GeV;
- break;
- }
- // generate y
- yTemp[i] = ymin + UseRandom::rnd()*(ymax-ymin);
- //generate x3 & x1 from pT & y
- double x1Plus = 1.;
- double x1Minus = 2.*mu_;
- x3Solution[i] = 2.*pT[i]*cosh(yTemp[i])/mZ_;
- // prefactor
- double weightPrefactor = 0.5/sqrt(1.-4.*mu2_)/overEst;
- // calculate x1 & x2 solutions
- Energy4 discrim2 = (sqr(x3Solution[i]*mZ_) - 4.*pT2)*
- (mz2_*(x3Solution[i]-1.)*(4.*mu2_+x3Solution[i]-1.)-4.*mu2_*pT2);
- //check discriminant2 is > 0
- if( discrim2 < ZERO) continue;
- Energy2 discriminant = sqrt(discrim2);
- Energy2 fact1 = 3.*mz2_*x3Solution[i]-2.*mz2_+2.*pT2*x3Solution[i]
- -4.*pT2-mz2_*sqr(x3Solution[i]);
- Energy2 fact2 = 2.*mz2_*(x3Solution[i]-1.)-2.*pT2;
- // two solns for x1
- x1Solution[i][0] = (fact1 + discriminant)/fact2;
- x1Solution[i][1] = (fact1 - discriminant)/fact2;
-
- bool found = false;
- for(unsigned int j=0;j<2;++j) {
- x2Solution[i][0] = 2.-x3Solution[i]-x1Solution[i][0];
- x2Solution[i][1] = 2.-x3Solution[i]-x1Solution[i][1];
- // set limits on x2
- double root = max(0.,sqr(x1Solution[i][j])-4.*mu2_);
- root = sqrt(root);
- double x2Plus = 1.-0.5*(1.-x1Solution[i][j])/(1.-x1Solution[i][j]+mu2_)
- *(x1Solution[i][j]-2.*mu2_-root);
- double x2Minus = 1.-0.5*(1.-x1Solution[i][j])/(1.-x1Solution[i][j]+mu2_)
- *(x1Solution[i][j]-2.*mu2_+root);
- if(x1Solution[i][j]>=x1Minus && x1Solution[i][j]<=x1Plus &&
- x2Solution[i][j]>=x2Minus && x2Solution[i][j]<=x2Plus &&
- checkZMomenta(x1Solution[i][j], x2Solution[i][j], x3Solution[i], yTemp[i], pT[i])) {
- probTemp[i][j] = weightPrefactor*pT[i]*
- calculateJacobian(x1Solution[i][j], x2Solution[i][j], pT[i])*
- calculateRealEmission(x1Solution[i][j], x2Solution[i][j],
- hardProcess, phi, false, i);
- found = true;
- }
- else {
- probTemp[i][j] = 0.;
- }
- }
- if(!found) continue;
- // alpha S piece
- double wgt = (probTemp[i][0]+probTemp[i][1])*alphaS()->ratio(sqr(pT[i]));
- // matrix element weight
- reject = UseRandom::rnd()>wgt;
- }
- while(reject);
- } //end of emitter for loop
- // no emission
- if(pT[0]<ZERO&&pT[1]<ZERO) return false;
- //pick the spectator and x1 x2 values
- double x1,x2,y;
- //particle 1 emits, particle 2 spectates
- unsigned int iemit=0;
- if(pT[0]>pT[1]){
- pT_ = pT[0];
- y=yTemp[0];
- if(probTemp[0][0]>UseRandom::rnd()*(probTemp[0][0]+probTemp[0][1])) {
- x1 = x1Solution[0][0];
- x2 = x2Solution[0][0];
- }
- else {
- x1 = x1Solution[0][1];
- x2 = x2Solution[0][1];
- }
- }
- // particle 2 emits, particle 1 spectates
- else {
- iemit=1;
- pT_ = pT[1];
- y=yTemp[1];
- if(probTemp[1][0]>UseRandom::rnd()*(probTemp[1][0]+probTemp[1][1])) {
- x1 = x1Solution[1][0];
- x2 = x2Solution[1][0];
- }
- else {
- x1 = x1Solution[1][1];
- x2 = x2Solution[1][1];
- }
- }
- // find spectator
- unsigned int ispect = iemit == 0 ? 1 : 0;
- // Find the boost from the lab to the c.o.m with the spectator
- // along the -z axis, and then invert it.
- LorentzRotation eventFrame( ( quark_[0] + quark_[1] ).findBoostToCM() );
- Lorentz5Momentum spectator = eventFrame*quark_[ispect];
- eventFrame.rotateZ( -spectator.phi() );
- eventFrame.rotateY( -spectator.theta() - Constants::pi );
- eventFrame.invert();
- // spectator
- quark_[ispect].setT( 0.5*x2*mZ_ );
- quark_[ispect].setX( ZERO );
- quark_[ispect].setY( ZERO );
- quark_[ispect].setZ( -sqrt(0.25*mz2_*x2*x2-mz2_*mu2_) );
- // gluon
- gauge_.setT( pT_*cosh(y) );
- gauge_.setX( pT_*cos(phi) );
- gauge_.setY( pT_*sin(phi) );
- gauge_.setZ( pT_*sinh(y) );
- gauge_.setMass(ZERO);
- // emitter reconstructed from gluon & spectator
- quark_[iemit] = - gauge_ - quark_[ispect];
- quark_[iemit].setT( 0.5*mZ_*x1 );
- // boost constructed vectors into the event frame
- quark_[0] = eventFrame * quark_[0];
- quark_[1] = eventFrame * quark_[1];
- gauge_ = eventFrame * gauge_;
- // need to reset masses because for whatever reason the boost
- // touches the mass component of the five-vector and can make
- // zero mass objects acquire a floating point negative mass(!).
- gauge_.setMass( ZERO );
- quark_[iemit] .setMass(partons_[iemit ]->mass());
- quark_[ispect].setMass(partons_[ispect]->mass());
-
- return true;
-}
-
-InvEnergy SMZDecayer::calculateJacobian(double x1, double x2, Energy pT) const{
- double xPerp = abs(2.*pT/mZ_);
- Energy jac = mZ_*fabs((x1*x2-2.*mu2_*(x1+x2)+sqr(x2)-x2)/xPerp/pow(sqr(x2)-4.*mu2_,1.5));
- return 1./jac; //jacobian as defined is dptdy=jac*dx1dx2, therefore we have to divide by it
-}
-
-bool SMZDecayer::checkZMomenta(double x1, double x2, double x3, double y, Energy pT) const {
- double xPerp2 = 4.*pT*pT/mZ_/mZ_;
- static double tolerance = 1e-6;
- bool isMomentaReconstructed = false;
-
- if(pT*sinh(y)>ZERO) {
- if(abs(-sqrt(sqr(x2)-4.*mu2_)+sqrt(sqr(x3)-xPerp2) + sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance ||
- abs(-sqrt(sqr(x2)-4.*mu2_)+sqrt(sqr(x3)-xPerp2) - sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance) isMomentaReconstructed=true;
- }
- else if(pT*sinh(y) < ZERO){
- if(abs(-sqrt(sqr(x2)-4.*mu2_)-sqrt(sqr(x3)-xPerp2) + sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance ||
- abs(-sqrt(sqr(x2)-4.*mu2_)-sqrt(sqr(x3)-xPerp2) - sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance) isMomentaReconstructed=true;
- }
- else
- if(abs(-sqrt(sqr(x2)-4.*mu2_)+ sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance) isMomentaReconstructed=true;
-
- return isMomentaReconstructed;
-}
-
double SMZDecayer::calculateRealEmission(double x1, double x2,
- vector<PPtr> hardProcess,
- double phi,
- bool subtract) const {
+ vector<PPtr> hardProcess,
+ double phi,
+ bool subtract) const {
// make partons data object for meRatio
vector<cPDPtr> partons (3);
for(int ix=0; ix<3; ++ix)
partons[ix] = hardProcess[ix]->dataPtr();
partons.push_back(gluon_);
// calculate x3
double x3 = 2.-x1-x2;
double xT = sqrt(max(0.,sqr(x3) -0.25*sqr(sqr(x2)+sqr(x3)-sqr(x1))/(sqr(x2)-4.*mu2_)));
// calculate the momenta
Energy M = mZ_;
Lorentz5Momentum pspect(ZERO,ZERO,-0.5*M*sqrt(max(sqr(x2)-4.*mu2_,0.)),0.5*M*x2,M*mu_);
Lorentz5Momentum pemit (-0.5*M*xT*cos(phi),-0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x1)-sqr(xT)-4.*mu2_,0.)),0.5*M*x1,M*mu_);
Lorentz5Momentum pgluon(0.5*M*xT*cos(phi), 0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x3)-sqr(xT),0.)),0.5*M*x3,ZERO);
if(abs(pspect.z()+pemit.z()-pgluon.z())/M<1e-6)
pgluon.setZ(-pgluon.z());
else if(abs(pspect.z()-pemit.z()+pgluon.z())/M<1e-6)
pemit .setZ(- pemit.z());
// loop over the possible emitting partons
double realwgt(0.);
for(unsigned int iemit=0;iemit<2;++iemit) {
// boost and rotate momenta
LorentzRotation eventFrame( ( hardProcess[1]->momentum() +
hardProcess[2]->momentum() ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*hardProcess[iemit+1]->momentum();
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() );
eventFrame.invert();
vector<Lorentz5Momentum> momenta(3);
momenta[0] = hardProcess[0]->momentum();
if(iemit==0) {
momenta[2] = eventFrame*pspect;
momenta[1] = eventFrame*pemit ;
}
else {
momenta[1] = eventFrame*pspect;
momenta[2] = eventFrame*pemit ;
}
momenta.push_back(eventFrame*pgluon);
// calculate the weight
if(1.-x1>1e-5 && 1.-x2>1e-5)
realwgt += meRatio(partons,momenta,iemit,subtract);
}
// total real emission contribution
return realwgt;
}
double SMZDecayer::calculateRealEmission(double x1, double x2,
- vector<PPtr> hardProcess,
- double phi,
- bool subtract,
- int emitter) const {
+ vector<PPtr> hardProcess,
+ double phi,
+ bool subtract,
+ int emitter) const {
// make partons data object for meRatio
vector<cPDPtr> partons (3);
for(int ix=0; ix<3; ++ix)
partons[ix] = hardProcess[ix]->dataPtr();
partons.push_back(gluon_);
// calculate x3
double x3 = 2.-x1-x2;
double xT = sqrt(max(0.,sqr(x3) -0.25*sqr(sqr(x2)+sqr(x3)-sqr(x1))/(sqr(x2)-4.*mu2_)));
// calculate the momenta
Energy M = mZ_;
Lorentz5Momentum pspect(ZERO,ZERO,-0.5*M*sqrt(max(sqr(x2)-4.*mu2_,0.)),0.5*M*x2,M*mu_);
Lorentz5Momentum pemit (-0.5*M*xT*cos(phi),-0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x1)-sqr(xT)-4.*mu2_,0.)),0.5*M*x1,M*mu_);
Lorentz5Momentum pgluon( 0.5*M*xT*cos(phi), 0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x3)-sqr(xT),0.)),0.5*M*x3,ZERO);
if(abs(pspect.z()+pemit.z()-pgluon.z())/M<1e-6)
pgluon.setZ(-pgluon.z());
else if(abs(pspect.z()-pemit.z()+pgluon.z())/M<1e-6)
pemit .setZ(- pemit.z());
- // loop over the possible emitting partons
- double realwgt(0.);
-
// boost and rotate momenta
LorentzRotation eventFrame( ( hardProcess[1]->momentum() +
hardProcess[2]->momentum() ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*hardProcess[emitter+1]->momentum();
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() );
eventFrame.invert();
vector<Lorentz5Momentum> momenta(3);
momenta[0] = hardProcess[0]->momentum();
if(emitter==0) {
momenta[2] = eventFrame*pspect;
momenta[1] = eventFrame*pemit ;
}
else {
momenta[1] = eventFrame*pspect;
momenta[2] = eventFrame*pemit ;
}
momenta.push_back(eventFrame*pgluon);
// calculate the weight
+ double realwgt(0.);
if(1.-x1>1e-5 && 1.-x2>1e-5)
- realwgt += meRatio(partons,momenta,emitter,subtract);
+ realwgt = meRatio(partons,momenta,emitter,subtract);
// total real emission contribution
return realwgt;
}
diff --git a/Decay/Perturbative/SMZDecayer.h b/Decay/Perturbative/SMZDecayer.h
--- a/Decay/Perturbative/SMZDecayer.h
+++ b/Decay/Perturbative/SMZDecayer.h
@@ -1,562 +1,545 @@
// -*- C++ -*-
//
// SMZDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SMZDecayer_H
#define HERWIG_SMZDecayer_H
//
// This is the declaration of the SMZDecayer class.
//
#include "Herwig/Decay/PerturbativeDecayer.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "Herwig/Decay/DecayPhaseSpaceMode.h"
-#include "Herwig/Shower/Core/Couplings/ShowerAlpha.fh"
namespace Herwig {
using namespace ThePEG;
using namespace ThePEG::Helicity;
/** \ingroup Decay
*
* The <code>SMZDecayer</code> is designed to perform the decay of the
* Z boson to the Standard Model fermions. In principle it can also
* be used for these decays in any model.
*
* @see PerturbativeDecayer
*
*/
class SMZDecayer: public PerturbativeDecayer {
public:
/**
* Default constructor.
*/
SMZDecayer();
/**
* Virtual members to be overridden by inheriting classes
* which implement hard corrections
*/
//@{
/**
* Has an old fashioned ME correction
*/
virtual bool hasMECorrection() {return true;}
/**
* Initialize the ME correction
*/
virtual void initializeMECorrection(RealEmissionProcessPtr , double & ,
double & );
/**
* Apply the hard matrix element correction to a given hard process or decay
*/
virtual RealEmissionProcessPtr applyHardMatrixElementCorrection(RealEmissionProcessPtr);
/**
* Apply the soft matrix element correction
* @param initial The particle from the hard process which started the
* shower
* @param parent The initial particle in the current branching
* @param br The branching struct
* @return If true the emission should be vetoed
*/
virtual bool softMatrixElementVeto(ShowerProgenitorPtr initial,
ShowerParticlePtr parent,Branching br);
/**
* Has a POWHEG style correction
*/
- virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
-
- /**
- * Apply the POWHEG style correction
- */
- virtual RealEmissionProcessPtr generateHardest(RealEmissionProcessPtr);
+ virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
//@}
public:
/**
* Which of the possible decays is required
* @param cc Is this mode the charge conjugate
* @param parent The decaying particle
* @param children The decay products
*/
virtual int modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const;
/**
* Return the matrix element squared for a given mode and phase-space channel.
* @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
* @param decay The particles produced in the decay.
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
virtual double me2(const int ichan, const Particle & part,
const ParticleVector & decay,MEOption meopt) const;
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
*/
virtual void dataBaseOutput(ofstream & os,bool header) const;
/**
* Members for the generation of QED radiation in the decays
*/
//@{
/**
* The one-loop virtual correction.
* @param imode The mode required.
* @param part The decaying particle.
* @param products The decay products including the radiated photon.
* @return Whether the correction is implemented
*/
virtual double oneLoopVirtualME(unsigned int imode,
const Particle & part,
const ParticleVector & products);
/**
* The real emission matrix element
* @param imode The mode required
* @param part The decaying particle
* @param products The decay products including the radiated photon
* @param iemitter The particle which emitted the photon
* @param ctheta The cosine of the polar angle between the photon and the
* emitter
* @param stheta The sine of the polar angle between the photon and the
* emitter
* @param rot1 Rotation from rest frame to frame for real emission
* @param rot2 Rotation to place emitting particle along z
*/
virtual InvEnergy2 realEmissionME(unsigned int imode,
const Particle & part,
ParticleVector & products,
unsigned int iemitter,
double ctheta, double stheta,
const LorentzRotation & rot1,
const LorentzRotation & rot2);
//@}
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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
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 {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving and
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/**
* Apply the hard matrix element
*/
vector<Lorentz5Momentum> applyHard(const ParticleVector &p);
/**
* Get the weight for hard emission
*/
double getHard(double &, double &);
/**
* Set the \f$\rho\f$ parameter
*/
void setRho(double);
/**
* Set the \f$\tilde{\kappa}\f$ parameters symmetrically
*/
void setKtildeSymm();
/**
* Set second \f$\tilde{\kappa}\f$, given the first.
*/
void setKtilde2();
/**
* Translate the variables from \f$x_q,x_{\bar{q}}\f$ to \f$\tilde{\kappa},z\f$
*/
//@{
/**
* Calculate \f$z\f$.
*/
double getZfromX(double, double);
/**
* Calculate \f$\tilde{\kappa}\f$.
*/
double getKfromX(double, double);
//@}
/**
* Calculate \f$x_{q},x_{\bar{q}}\f$ from \f$\tilde{\kappa},z\f$.
* @param kt \f$\tilde{\kappa}\f$
* @param z \f$z\f$
* @param x \f$x_{q}\f$
* @param xbar \f$x_{\bar{q}}\f$
*/
void getXXbar(double kt, double z, double & x, double & xbar);
/**
* Soft weight
*/
//@{
/**
* Soft quark weight calculated from \f$x_{q},x_{\bar{q}}\f$
* @param x \f$x_{q}\f$
* @param xbar \f$x_{\bar{q}}\f$
*/
double qWeight(double x, double xbar);
/**
* Soft antiquark weight calculated from \f$x_{q},x_{\bar{q}}\f$
* @param x \f$x_{q}\f$
* @param xbar \f$x_{\bar{q}}\f$
*/
double qbarWeight(double x, double xbar);
/**
* Soft quark weight calculated from \f$\tilde{q},z\f$
* @param qtilde \f$\tilde{q}\f$
* @param z \f$z\f$
*/
double qWeightX(Energy qtilde, double z);
/**
* Soft antiquark weight calculated from \f$\tilde{q},z\f$
* @param qtilde \f$\tilde{q}\f$
* @param z \f$z\f$
*/
double qbarWeightX(Energy qtilde, double z);
//@}
/**
* ????
*/
double u(double);
/**
* Vector and axial vector parts of the matrix element
*/
//@{
/**
* Vector part of the matrix element
*/
double MEV(double, double);
/**
* Axial vector part of the matrix element
*/
double MEA(double, double);
/**
* The matrix element, given \f$x_1\f$, \f$x_2\f$.
* @param x1 \f$x_1\f$
* @param x2 \f$x_2\f$
*/
double PS(double x1, double x2);
-
- /**
- * Access to the strong coupling
- */
- ShowerAlphaPtr alphaS() const {return alpha_;}
//@}
protected:
/**
* Real emission term, for use in generating the hardest emission
*/
double calculateRealEmission(double x1, double x2,
vector<PPtr> hardProcess,
double phi,
bool subtract,
int emitter) const;
/**
* Real emission term, for use in generating the hardest emission
*/
double calculateRealEmission(double x1, double x2,
vector<PPtr> hardProcess,
double phi,
bool subtract) const;
/**
* Check the sign of the momentum in the \f$z\f$-direction is correct.
*/
bool checkZMomenta(double x1, double x2, double x3, double y, Energy pT) const;
/**
* Calculate the Jacobian
*/
InvEnergy calculateJacobian(double x1, double x2, Energy pT) const;
/**
* Calculate the ratio between NLO & LO ME
*/
double meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter,bool subtract) const;
+
+ /**
+ * Calculate matrix element ratio R/B
+ */
+ virtual double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
+ const ParticleVector & decay3, MEOption meopt);
+
/**
* Calculate the LO ME
*/
double loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const;
/**
* Calculate the NLO real emission piece of ME
*/
InvEnergy2 realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const;
/**
* Generate a real emission event
*/
bool getEvent(vector<PPtr> hardProcess);
private:
/**
* Private and non-existent assignment operator.
*/
SMZDecayer & operator=(const SMZDecayer &);
private:
/**
* Pointer to the fermion-antifermion Z vertex
*/
FFVVertexPtr FFZVertex_;
/**
* Pointer to the photon vertex
*/
AbstractFFVVertexPtr FFPVertex_;
/**
* Pointer to the fermion-antifermion G vertex
*/
AbstractFFVVertexPtr FFGVertex_;
/**
* maximum weights for the different integrations
*/
//@{
/**
* Weights for the Z to quarks decays.
*/
vector<double> quarkWeight_;
/**
* Weights for the Z to leptons decays.
*/
vector<double> leptonWeight_;
//@}
/**
* Spin density matrix for the decay
*/
mutable RhoDMatrix _rho;
/**
* Polarization vectors for the decay
*/
mutable vector<VectorWaveFunction> _vectors;
/**
* Spinors for the decay
*/
mutable vector<SpinorWaveFunction> _wave;
/**
* Barred spinors for the decay
*/
mutable vector<SpinorBarWaveFunction> _wavebar;
private:
/**
* CM energy
*/
Energy d_Q_;
/**
* Quark mass
*/
Energy d_m_;
/**
* The rho parameter
*/
double d_rho_;
/**
* The v parameter
*/
double d_v_;
/**
* The initial kappa-tilde values for radiation from the quark
*/
double d_kt1_;
/**
* The initial kappa-tilde values for radiation from the antiquark
*/
double d_kt2_;
/**
* Cut-off parameter
*/
static const double EPS_;
- /**
- * Pointer to the coupling
- */
- ShowerAlphaPtr alpha_;
-
private:
/**
* The colour factor
*/
double CF_;
/**
* The Z mass
*/
mutable Energy mZ_;
/**
* The reduced mass
*/
mutable double mu_;
/**
* The square of the reduced mass
*/
mutable double mu2_;
/**
* The strong coupling
*/
mutable double aS_;
/**
* The scale
*/
mutable Energy2 scale_;
/**
* Stuff for the POWHEG correction
*/
//@{
/**
* ParticleData object for the gluon
*/
tcPDPtr gluon_;
/**
- * The cut off on pt, assuming massless quarks.
- */
- Energy pTmin_;
-
- // radiative variables (pt,y)
- Energy pT_;
-
- /**
* The ParticleData objects for the fermions
*/
vector<tcPDPtr> partons_;
/**
* The fermion momenta
*/
vector<Lorentz5Momentum> quark_;
/**
* The momentum of the radiated gauge boson
*/
Lorentz5Momentum gauge_;
/**
* The Z boson
*/
PPtr zboson_;
/**
* Higgs mass squared
*/
Energy2 mz2_;
//@}
/**
* Whether or not to give an LO or NLO normalisation
*/
bool NLO_;
};
}
#endif /* HERWIG_SMZDecayer_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:09 PM (21 h, 16 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023392
Default Alt Text
(113 KB)

Event Timeline