diff --git a/Decay/Perturbative/SMHiggsFermionsDecayer.cc b/Decay/Perturbative/SMHiggsFermionsDecayer.cc
--- a/Decay/Perturbative/SMHiggsFermionsDecayer.cc
+++ b/Decay/Perturbative/SMHiggsFermionsDecayer.cc
@@ -1,645 +1,405 @@
// -*- C++ -*-
//
// SMHiggsFermionsDecayer.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 SMHiggsFermionsDecayer class.
//
#include "SMHiggsFermionsDecayer.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/ScalarSpinInfo.h"
#include "ThePEG/Helicity/FermionSpinInfo.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include "Herwig/Utilities/Maths.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Shower/Core/Couplings/ShowerAlpha.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
SMHiggsFermionsDecayer::SMHiggsFermionsDecayer() :
CF_(4./3.), pTmin_(1.*GeV), NLO_(false) {
_maxwgt.resize(9);
_maxwgt[0]=0.;
_maxwgt[1]=0;
_maxwgt[2]=0;
_maxwgt[3]=0.0194397;
_maxwgt[4]=0.463542;
_maxwgt[5]=0.;
_maxwgt[6]=6.7048e-09;
_maxwgt[7]=0.00028665;
_maxwgt[8]=0.0809643;
}
void SMHiggsFermionsDecayer::doinit() {
PerturbativeDecayer::doinit();
// get the vertices from the Standard Model object
tcHwSMPtr hwsm=dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm)
throw InitException() << "SMHiggsFermionsDecayer needs the StandardModel class"
<< " to be either the Herwig one or a class inheriting"
<< " from it";
_hvertex = hwsm->vertexFFH();
// make sure they are initialized
_hvertex->init();
// get the width generator for the higgs
tPDPtr higgs = getParticleData(ParticleID::h0);
// set up the decay modes
vector<double> wgt(0);
unsigned int imode=0;
tPDVector extpart(3);
DecayPhaseSpaceModePtr mode;
int iy;
extpart[0]=higgs;
for(unsigned int istep=0;istep<11;istep+=10) {
for(unsigned ix=1;ix<7;++ix) {
if(istep<10||ix%2!=0) {
iy = ix+istep;
extpart[1]=getParticleData( iy);
extpart[2]=getParticleData(-iy);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
addMode(mode,_maxwgt[imode],wgt);
++imode;
}
}
}
gluon_ = getParticleData(ParticleID::g);
// Energy quarkMass = getParticleData(ParticleID::b )->mass();
// Energy higgsMass = getParticleData(ParticleID::h0)->mass();
// double mu = quarkMass/higgsMass;
// double beta = sqrt(1.-4.*sqr(mu));
// double beta2 = sqr(beta);
// double aS = SM().alphaS(sqr(higgsMass));
// double L = log((1.+beta)/(1.-beta));
// cerr << "testing " << beta << " " << mu << "\n";
// cerr << "testing " << aS << " " << L << "\n";
// double fact =
// 6.-0.75*(1.+beta2)/beta2+12.*log(mu)-8.*log(beta)
// +(5./beta-2.*beta+0.375*sqr(1.-beta2)/beta2/beta)*L
// +(1.+beta2)/beta*(4.*L*log(0.5*(1.+beta)/beta)
// -2.*log(0.5*(1.+beta))*log(0.5*(1.-beta))
// +8.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))
// -4.*Herwig::Math::ReLi2(0.5*(1.-beta)));
// cerr << "testing correction "
// << 1.+4./3.*aS/Constants::twopi*fact
// << "\n";
// double real = 4./3.*aS/Constants::twopi*
// (8.-0.75*(1.+beta2)/beta2+8.*log(mu)-8.*log(beta)
// +(3./beta+0.375*sqr(1.-beta2)/pow(beta,3))*L
// +(1.+beta2)/beta*(-0.5*sqr(L)+4.*L*log(0.5*(1.+beta))
// -2.*L*log(beta)-2.*log(0.5*(1.+beta))*log(0.5*(1.-beta))
// +6.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))
// -4.*Herwig::Math::ReLi2(0.5*(1.-beta))
// -2./3.*sqr(Constants::pi)));
// double virt = 4./3.*aS/Constants::twopi*
// (-2.+4.*log(mu)+(2./beta-2.*beta)*L
// +(1.+beta2)/beta*(0.5*sqr(L)-2.*L*log(beta)+2.*sqr(Constants::pi)/3.
// +2.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))));
// cerr << "testing real " << real << "\n";
// cerr << "testing virtual " << virt << "\n";
// cerr << "testing total no mb corr " << 1.+real+virt << "\n";
// cerr << "testing total mb corr " << 1.+real+virt +(8./3. - 2.*log(sqr(mu)))*aS/Constants::pi << "\n";
// InvEnergy2 Gf = 1.166371e-5/GeV2;
// Gf = sqrt(2.)*4*Constants::pi*SM().alphaEM(sqr(higgsMass))/8./SM().sin2ThetaW()/
// sqr(getParticleData(ParticleID::Wplus)->mass());
// cerr << "testing GF " << Gf*GeV2 << "\n";
// Energy LO = (3./8./Constants::pi)*sqrt(2)*sqr(quarkMass)*Gf*higgsMass*beta*beta*beta;
// cerr << "testing LO " << LO/GeV << "\n";
// cerr << "testing quark mass " << quarkMass/GeV << "\n";
// cerr << "testing gamma " << (1.+real+virt)*LO/MeV << "\n";
}
bool SMHiggsFermionsDecayer::accept(tcPDPtr parent, const tPDVector & children) const {
if(parent->id()!=ParticleID::h0||children.size()!=2) return false;
tPDVector::const_iterator pit = children.begin();
int id1=(**pit).id();
++pit;
int id2=(**pit).id();
if(id1==-id2&&(abs(id1)<=6||(abs(id1)>=11&&abs(id1)<=16)))
return true;
else
return false;
}
ParticleVector SMHiggsFermionsDecayer::decay(const Particle & parent,
const tPDVector & children) const {
// id's of the decaying particles
tPDVector::const_iterator pit(children.begin());
int id1((**pit).id());
int imode=-1;
if(abs(id1)<=6) imode = abs(id1)-1;
else if(abs(id1)>=11&&abs(id1)<=16) imode = (abs(id1)-11)/2+6;
ParticleVector output(generate(false,false,imode,parent));
// set up the colour flow
if(output[0]->hasColour()) output[0]->antiColourNeighbour(output[1]);
else if(output[1]->hasColour()) output[1]->antiColourNeighbour(output[0]);
return output;
}
void SMHiggsFermionsDecayer::persistentOutput(PersistentOStream & os) const {
os << _maxwgt << _hvertex << NLO_
<< alphaS_ << gluon_ << ounit( pTmin_, GeV );
}
void SMHiggsFermionsDecayer::persistentInput(PersistentIStream & is, int) {
is >> _maxwgt >> _hvertex >> NLO_
>> alphaS_ >> gluon_ >> iunit( pTmin_, GeV );
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SMHiggsFermionsDecayer,PerturbativeDecayer>
describeHerwigSMHiggsFermionsDecayer("Herwig::SMHiggsFermionsDecayer", "HwPerturbativeHiggsDecay.so");
void SMHiggsFermionsDecayer::Init() {
static ClassDocumentation<SMHiggsFermionsDecayer> documentation
("The SMHiggsFermionsDecayer class implements the decat of the Standard Model"
" Higgs boson to the Standard Model fermions");
static ParVector<SMHiggsFermionsDecayer,double> interfaceMaxWeights
("MaxWeights",
"Maximum weights for the various decays",
&SMHiggsFermionsDecayer::_maxwgt, 9, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static Reference<SMHiggsFermionsDecayer,ShowerAlpha> interfaceCoupling
("Coupling",
"The object calculating the strong coupling constant",
&SMHiggsFermionsDecayer::alphaS_, false, false, true, false, false);
static Parameter<SMHiggsFermionsDecayer, Energy> interfacePtMin
("minpT",
"The pt cut on hardest emision generation",
&SMHiggsFermionsDecayer::pTmin_, GeV, 1.*GeV, 0*GeV, 100000.0*GeV,
false, false, Interface::limited);
static Switch<SMHiggsFermionsDecayer,bool> interfaceNLO
("NLO",
"Whether to return the LO or NLO result",
&SMHiggsFermionsDecayer::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 SMHiggsFermionsDecayer::me2(const int, const Particle & part,
const ParticleVector & decay,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half)));
int iferm(1),ianti(0);
if(decay[0]->id()>0) swap(iferm,ianti);
if(meopt==Initialize) {
ScalarWaveFunction::
calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
_swave = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming);
}
if(meopt==Terminate) {
ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
incoming,true);
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);
Energy2 scale(sqr(part.mass()));
unsigned int ifm,ia;
for(ifm=0;ifm<2;++ifm) {
for(ia=0;ia<2;++ia) {
if(iferm>ianti)
(*ME())(0,ia,ifm)=_hvertex->evaluate(scale,_wave[ia],
_wavebar[ifm],_swave);
else
(*ME())(0,ifm,ia)=_hvertex->evaluate(scale,_wave[ia],
_wavebar[ifm],_swave);
}
}
int id = abs(decay[0]->id());
double output=(ME()->contract(_rho)).real()*UnitRemoval::E2/scale;
if(id <=6) output*=3.;
// test of the partial width
// Ptr<Herwig::StandardModel>::transient_const_pointer
// hwsm=dynamic_ptr_cast<Ptr<Herwig::StandardModel>::transient_const_pointer>(standardModel());
// double g2(hwsm->alphaEM(scale)*4.*Constants::pi/hwsm->sin2ThetaW());
// Energy mass(hwsm->mass(scale,decay[0]->dataPtr())),
// mw(getParticleData(ParticleID::Wplus)->mass());
// double beta(sqrt(1.-4.*decay[0]->mass()*decay[0]->mass()/scale));
// cerr << "testing alpha " << hwsm->alphaEM(scale) << "\n";
// Energy test(g2*mass*mass*beta*beta*beta*part.mass()/32./Constants::pi/mw/mw);
// if(abs(decay[0]->id())<=6){test *=3.;}
// cout << "testing the answer " << output << " "
// << test/GeV
// << endl;
// leading-order result
if(!NLO_) return output;
// fermion mass
Energy particleMass = decay[0]->dataPtr()->mass();
// check decay products coloured, otherwise return
if(!decay[0]->dataPtr()->coloured()||
particleMass==ZERO) return output;
// inital masses, couplings etc
// higgs mass
mHiggs_ = part.mass();
// strong coupling
aS_ = SM().alphaS(sqr(mHiggs_));
// reduced mass
mu_ = particleMass/mHiggs_;
mu2_ = sqr(mu_);
// generate y
double yminus = 0.;
double yplus = 1.-2.*mu_*(1.-mu_)/(1.-2*mu2_);
double y = yminus + UseRandom::rnd()*(yplus-yminus);
//generate z for D31,2
double v = sqrt(sqr(2.*mu2_+(1.-2.*mu2_)*(1.-y))-4.*mu2_)/(1.-2.*mu2_)/(1.-y);
double zplus = (1.+v)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
double zminus = (1.-v)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
double z = zminus + UseRandom::rnd()*(zplus-zminus);
// map y,z to x1,x2 for both possible emissions
double x2 = 1. - y*(1.-2.*mu2_);
double x1 = 1. - z*(x2-2.*mu2_);
//get the dipoles
InvEnergy2 D1 = dipoleSubtractionTerm( x1, x2);
InvEnergy2 D2 = dipoleSubtractionTerm( x2, x1);
InvEnergy2 dipoleSum = abs(D1) + abs(D2);
//jacobian
double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
//calculate real
Energy2 realPrefactor = 0.25*sqr(mHiggs_)*sqr(1.-2.*mu2_)
/sqrt(calculateLambda(1,mu2_,mu2_))/sqr(Constants::twopi);
- InvEnergy2 realEmission = calculateRealEmission( x1, x2);
+ InvEnergy2 realEmission = 4.*Constants::pi*aS_*calculateRealEmission( x1, x2);
// calculate the virtual
double virtualTerm = calculateVirtualTerm();
// running mass correction
virtualTerm += (8./3. - 2.*log(mu2_))*aS_/Constants::pi;
//answer = (born + virtual + real)/born * LO
output *= 1. + virtualTerm + 2.*jac*realPrefactor*(realEmission*abs(D1)/dipoleSum - D1);
// return the answer
return output;
}
void SMHiggsFermionsDecayer::dataBaseOutput(ofstream & os,bool header) const {
if(header) os << "update decayers set parameters=\"";
// parameters for the PerturbativeDecayer base class
for(unsigned int ix=0;ix<_maxwgt.size();++ix) {
os << "newdef " << name() << ":MaxWeights " << ix << " "
<< _maxwgt[ix] << "\n";
}
PerturbativeDecayer::dataBaseOutput(os,false);
if(header) os << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
void SMHiggsFermionsDecayer::doinitrun() {
PerturbativeDecayer::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<numberModes();++ix) {
_maxwgt[ix] = mode(ix)->maxWeight();
}
}
}
-RealEmissionProcessPtr SMHiggsFermionsDecayer::
-generateHardest(RealEmissionProcessPtr born) {
- // check coloured
- if(!born->bornOutgoing()[0]->dataPtr()->coloured()) return RealEmissionProcessPtr();
- assert(born->bornOutgoing().size()==2);
- // extract required info
- higgs_ = born->bornIncoming()[0];
- partons_.resize(2);
- quark_.resize(2);
- 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());
- }
- bool order = partons_[0]->id()<0;
- if(order) {
- swap(partons_[0] ,partons_[1] );
- swap(quark_[0] ,quark_[1] );
- }
- gauge_.setMass(0.*MeV);
- // Get the Higgs boson mass.
- mh2_ = (quark_[0] + quark_[1]).m2();
- mHiggs_ = sqrt(mh2_);
- aS_ = SM().alphaS(sqr(mHiggs_));
- Energy particleMass = partons_[0]->mass();
- mu_ = particleMass/mHiggs_;
- mu2_ = sqr(mu_);
- // Generate emission and set _quark[0,1] and _gauge to be the
- // momenta of q, qbar and g after the hardest emission:
- if(!getEvent()) {
- born->pT()[ShowerInteraction::QCD] = pTmin_;
- return born;
- }
- // 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 hboson = higgs_->dataPtr()->produceParticle(higgs_->momentum());
- born->incoming().push_back(hboson);
- 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;
-}
//calculate lambda
double SMHiggsFermionsDecayer::calculateLambda(double x, double y, double z) const{
return sqr(x)+sqr(y)+sqr(z)-2.*x*y-2.*x*z-2.*y*z;
}
//calculates the dipole subtraction term for x1, D31,2 (Dij,k),
// 2 is the spectator anti-fermion and 3 is the gluon
InvEnergy2 SMHiggsFermionsDecayer::
dipoleSubtractionTerm(double x1, double x2) const{
InvEnergy2 commonPrefactor = CF_*8.*Constants::pi*aS_/sqr(mHiggs_);
return commonPrefactor/(1.-x2)*
(2.*(1.-2.*mu2_)/(2.-x1-x2)-
sqrt((1.-4.*mu2_)/(sqr(x2)-4.*mu2_))*
(x2-2.*mu2_)*(2.+(x1-1.)/(x2-2.*mu2_)+2.*mu2_/(1.-x2))/(1.-2.*mu2_));
}
//return ME for real emission
InvEnergy2 SMHiggsFermionsDecayer::
calculateRealEmission(double x1, double x2) const {
- InvEnergy2 prefactor = CF_*8.*Constants::pi*aS_/sqr(mHiggs_)/(1.-4.*mu2_);
+ InvEnergy2 prefactor = 2.*CF_/sqr(mHiggs_)/(1.-4.*mu2_);
return prefactor*(2. + (1.-x1)/(1.-x2) + (1.-x2)/(1.-x1)
+ 2.*(1.-2.*mu2_)*(1.-4.*mu2_)/(1.-x1)/(1.-x2)
- 2.*(1.-4.*mu2_)*(1./(1.-x2)+1./(1.-x1))
- 2.*mu2_*(1.-4.*mu2_)*(1./sqr(1.-x2)+1./sqr(1.-x1)));
}
double SMHiggsFermionsDecayer::
calculateVirtualTerm() const {
// logs and prefactors
double beta = sqrt(1.-4.*mu2_);
double L = log((1.+beta)/(1.-beta));
double prefactor = CF_*aS_/Constants::twopi;
// non-singlet piece
double nonSingletTerm = calculateNonSingletTerm(beta, L);
double virtualTerm =
-2.+4.*log(mu_)+(2./beta - 2.*beta)*L
+ (2.-4.*mu2_)/beta*(0.5*sqr(L) - 2.*L*log(beta)
+ 2.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))
+ 2.*sqr(Constants::pi)/3.);
double iEpsilonTerm =
2.*(3.-sqr(Constants::pi)/2. + 0.5*log(mu2_) - 1.5*log(1.-2.*mu2_)
-(1.-2.*mu2_)/beta*(0.5*sqr(L)+sqr(Constants::pi)/6.
-2.*L*log(1.-2.*mu2_))
+ nonSingletTerm);
return prefactor*(virtualTerm+iEpsilonTerm);
}
//non-singlet piece of I(epsilon) insertion operator
double SMHiggsFermionsDecayer::
calculateNonSingletTerm(double beta, double L) const {
return 1.5*log(1.-2.*mu2_)
+ (1.-2.*mu2_)/beta*(- 2.*L*log(4.*(1.-2.*mu2_)/sqr(1.+beta))+
+ 2.*Herwig::Math::ReLi2(sqr((1.-beta)/(1.+beta)))
- 2.*Herwig::Math::ReLi2(2.*beta/(1.+beta))
- sqr(Constants::pi)/6.)
+ log(1.-mu_)
- 2.*log(1.-2.*mu_)
- 2.*mu2_/(1.-2.*mu2_)*log(mu_/(1.-mu_))
- mu_/(1.-mu_)
+ 2.*mu_*(2*mu_-1.)/(1.-2.*mu2_)
+ 0.5*sqr(Constants::pi);
}
-bool SMHiggsFermionsDecayer::getEvent() {
- // max pT
- Energy pTmax = 0.5*sqrt(mh2_);
- // 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()*4./3./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.};
- for(int i=0; i<2; i++) {
- do {
- // reject the emission
- reject = true;
- // generate pt
- pT[i] *= pow(UseRandom::rnd(),1./(prefactor*(ymax-ymin)));
- 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])/mHiggs_;
- // prefactor
- Energy2 weightPrefactor = mh2_/16./sqr(Constants::pi)/sqrt(1.-4.*mu2_);
- weightPrefactor /= prefactor;
- // calculate x1 & x2 solutions
- Energy4 discrim2 = (sqr(x3Solution[i]*mHiggs_) - 4.*pT2)*
- (mh2_*(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.*mh2_*x3Solution[i]-2.*mh2_+2.*pT2*x3Solution[i]-4.*pT2-mh2_*sqr(x3Solution[i]);
- Energy2 fact2 = 2.*mh2_*(x3Solution[i]-1.)-2.*pT2;
- // two solns for x1
- x1Solution[i][0] = (fact1 + discriminant)/fact2;
- x1Solution[i][1] = (fact1 - discriminant)/fact2;
- x2Solution[i][0] = 2.-x3Solution[i]-x1Solution[i][0];
- x2Solution[i][1] = 2.-x3Solution[i]-x1Solution[i][1];
- bool found = false;
- for(unsigned int j=0;j<2;++j) {
- if(x1Solution[i][j]>=x1Minus && x1Solution[i][j]<=x1Plus &&
- checkZMomenta(x1Solution[i][j], x2Solution[i][j], x3Solution[i], yTemp[i], pT[i])) {
- InvEnergy2 D1 = dipoleSubtractionTerm( x1Solution[i][j], x2Solution[i][j]);
- InvEnergy2 D2 = dipoleSubtractionTerm( x2Solution[i][j], x1Solution[i][j]);
- double dipoleFactor = abs(D1)/(abs(D1) + abs(D2));
- probTemp[i][j] = weightPrefactor*pT[i]*dipoleFactor*
- calculateJacobian(x1Solution[i][j], x2Solution[i][j], pT[i])*
- calculateRealEmission(x1Solution[i][j], x2Solution[i][j]);
-
- found = true;
- }
- else {
- probTemp[i][j] = 0.;
- }
- }
- if(!found) continue;
- // alpha S piece
- double wgt = (probTemp[i][0]+probTemp[i][1])*alphaS_->value(sqr(pT[i]))/aS_;
- // 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();
- //generation of phi
- double phi = UseRandom::rnd() * Constants::twopi;
- // spectator
- quark_[ispect].setT( 0.5*x2*mHiggs_ );
- quark_[ispect].setX( ZERO );
- quark_[ispect].setY( ZERO );
- quark_[ispect].setZ( -sqrt(0.25*mh2_*x2*x2-mh2_*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*mHiggs_*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;
+double SMHiggsFermionsDecayer::matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
+ const ParticleVector & decay3, MEOption) {
+ mHiggs_ = inpart.mass();
+ mu_ = decay2[0]->mass()/mHiggs_;
+ mu2_ = sqr(mu_);
+ double x1 = 2.*decay3[0]->momentum().t()/mHiggs_;
+ double x2 = 2.*decay3[1]->momentum().t()/mHiggs_;
+ return calculateRealEmission(x1,x2)*4.*Constants::pi*sqr(mHiggs_);
}
-
-InvEnergy SMHiggsFermionsDecayer::calculateJacobian(double x1, double x2, Energy pT) const{
- double xPerp = abs(2.*pT/mHiggs_);
- Energy jac = mHiggs_*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 SMHiggsFermionsDecayer::checkZMomenta(double x1, double x2, double x3, double y, Energy pT) const {
- double xPerp2 = 4.*pT*pT/mHiggs_/mHiggs_;
- 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;
-}
-
diff --git a/Decay/Perturbative/SMHiggsFermionsDecayer.h b/Decay/Perturbative/SMHiggsFermionsDecayer.h
--- a/Decay/Perturbative/SMHiggsFermionsDecayer.h
+++ b/Decay/Perturbative/SMHiggsFermionsDecayer.h
@@ -1,311 +1,297 @@
// -*- C++ -*-
//
// SMHiggsFermionsDecayer.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_SMHiggsFermionsDecayer_H
#define HERWIG_SMHiggsFermionsDecayer_H
//
// This is the declaration of the SMHiggsFermionsDecayer class.
//
#include "Herwig/Decay/PerturbativeDecayer.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "Herwig/Decay/DecayPhaseSpaceMode.h"
#include "Herwig/Shower/Core/Couplings/ShowerAlpha.fh"
namespace Herwig {
using namespace ThePEG;
/**
* The SMHiggsFermionsDecayer class is designed to decay the Standard Model Higgs
* to the Standard Model fermions.
*
* @see PerturbativeDecayer
*/
class SMHiggsFermionsDecayer: public PerturbativeDecayer {
public:
/**
* The default constructor.
*/
SMHiggsFermionsDecayer();
/**
* Which of the possible decays is required
*/
virtual int modeNumber(bool & , tcPDPtr , const tPDVector & ) const {return -1;}
/**
* Check if this decayer can perfom the decay for a particular mode.
* Uses the modeNumber member but can be overridden
* @param parent The decaying particle
* @param children The decay products
*/
virtual bool accept(tcPDPtr parent, const tPDVector & children) const;
/**
* For a given decay mode and a given particle instance, perform the
* decay and return the decay products. As this is the base class this
* is not implemented.
* @return The vector of particles produced in the decay.
*/
virtual ParticleVector decay(const Particle & 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;
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
/**
- * Apply the POWHEG style correction
+ * Calculate matrix element ratio R/B
*/
- virtual RealEmissionProcessPtr generateHardest(RealEmissionProcessPtr);
+ virtual double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
+ const ParticleVector & decay3, MEOption meopt);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {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:
/**
* Calcluate the Kallen function
*/
double calculateLambda(double x, double y, double z) const;
/**
* Dipole subtraction term
*/
InvEnergy2 dipoleSubtractionTerm(double x1, double x2) const;
/**
* Real emission term
*/
InvEnergy2 calculateRealEmission(double x1, double x2) const;
/**
* Virtual term
*/
double calculateVirtualTerm() const;
/**
* Non-singlet term
*/
double calculateNonSingletTerm(double beta, double L) 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;
-
- /**
- * Generate a real emission event
- */
- bool getEvent();
-
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SMHiggsFermionsDecayer & operator=(const SMHiggsFermionsDecayer &);
private:
/**
* Pointer to the Higgs vertex
*/
AbstractFFSVertexPtr _hvertex;
/**
* maximum weights for the different decay modes
*/
vector<double> _maxwgt;
/**
* Spin density matrix
*/
mutable RhoDMatrix _rho;
/**
* Scalar wavefunction
*/
mutable ScalarWaveFunction _swave;
/**
* Spinor wavefunction
*/
mutable vector<SpinorWaveFunction> _wave;
/**
* Barred spinor wavefunction
*/
mutable vector<SpinorBarWaveFunction> _wavebar;
private:
/**
* The colour factor
*/
double CF_;
/**
* The Higgs mass
*/
mutable Energy mHiggs_;
/**
* The reduced mass
*/
mutable double mu_;
/**
* The square of the reduced mass
*/
mutable double mu2_;
/**
* The strong coupling
*/
mutable double aS_;
/**
* Stuff for the POWHEG correction
*/
//@{
/**
* Pointer to the object calculating the strong coupling
*/
ShowerAlphaPtr alphaS_;
/**
* 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 Higgs boson
*/
PPtr higgs_;
/**
* Higgs mass squared
*/
Energy2 mh2_;
//@}
/**
* LO or NLO ?
*/
bool NLO_;
};
}
#endif /* HERWIG_SMHiggsFermionsDecayer_H */
diff --git a/Decay/Perturbative/SMTopDecayer.cc b/Decay/Perturbative/SMTopDecayer.cc
--- a/Decay/Perturbative/SMTopDecayer.cc
+++ b/Decay/Perturbative/SMTopDecayer.cc
@@ -1,1198 +1,1198 @@
// -*- C++ -*-
//
// SMTopDecayer.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 SMTopDecayer class.
//
#include "SMTopDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Decay/DecayVertex.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig/PDT/ThreeBodyAllOn1IntegralCalculator.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Shower/Core/Base/ShowerProgenitor.h"
#include "Herwig/Shower/Core/Base/ShowerParticle.h"
#include "Herwig/Shower/Core/Base/Branching.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
SMTopDecayer::SMTopDecayer()
: _wquarkwgt(6,0.),_wleptonwgt(3,0.), _xg_sampling(1.5),
_initialenhance(1.), _finalenhance(2.3), _useMEforT2(true) {
_wleptonwgt[0] = 0.302583;
_wleptonwgt[1] = 0.301024;
_wleptonwgt[2] = 0.299548;
_wquarkwgt[0] = 0.851719;
_wquarkwgt[1] = 0.0450162;
_wquarkwgt[2] = 0.0456962;
_wquarkwgt[3] = 0.859839;
_wquarkwgt[4] = 3.9704e-06;
_wquarkwgt[5] = 0.000489657;
generateIntermediates(true);
}
bool SMTopDecayer::accept(tcPDPtr parent, const tPDVector & children) const {
if(abs(parent->id()) != ParticleID::t) return false;
int id0(0),id1(0),id2(0);
for(tPDVector::const_iterator it = children.begin();
it != children.end();++it) {
int id=(**it).id(),absid(abs(id));
if(absid==ParticleID::b&&double(id)/double(parent->id())>0) {
id0=id;
}
else {
switch (absid) {
case ParticleID::nu_e:
case ParticleID::nu_mu:
case ParticleID::nu_tau:
id1 = id;
break;
case ParticleID::eminus:
case ParticleID::muminus:
case ParticleID::tauminus:
id2 = id;
break;
case ParticleID::b:
case ParticleID::d:
case ParticleID::s:
id1 = id;
break;
case ParticleID::u:
case ParticleID::c:
id2=id;
break;
default :
break;
}
}
}
if(id0==0||id1==0||id2==0) return false;
if(double(id1)/double(id2)>0) return false;
return true;
}
ParticleVector SMTopDecayer::decay(const Particle & parent,
const tPDVector & children) const {
int id1(0),id2(0);
for(tPDVector::const_iterator it = children.begin();
it != children.end();++it) {
int id=(**it).id(),absid=abs(id);
if(absid == ParticleID::b && double(id)/double(parent.id())>0) continue;
//leptons
if(absid > 10 && absid%2==0) id1=absid;
if(absid > 10 && absid%2==1) id2=absid;
//quarks
if(absid < 10 && absid%2==0) id2=absid;
if(absid < 10 && absid%2==1) id1=absid;
}
unsigned int imode(0);
if(id2 >=11 && id2<=16) imode = (id1-12)/2;
else imode = id1+1+id2/2;
bool cc = parent.id() == ParticleID::tbar;
ParticleVector out(generate(true,cc,imode,parent));
//arrange colour flow
PPtr pparent=const_ptr_cast<PPtr>(&parent);
out[1]->incomingColour(pparent,out[1]->id()<0);
ParticleVector products = out[0]->children();
if(products[0]->hasColour())
products[0]->colourNeighbour(products[1],true);
else if(products[0]->hasAntiColour())
products[0]->colourNeighbour(products[1],false);
return out;
}
void SMTopDecayer::persistentOutput(PersistentOStream & os) const {
os << _wvertex << _wquarkwgt << _wleptonwgt << _wplus
<< _initialenhance << _finalenhance << _xg_sampling << _useMEforT2;
}
void SMTopDecayer::persistentInput(PersistentIStream & is, int) {
is >> _wvertex >> _wquarkwgt >> _wleptonwgt >> _wplus
>> _initialenhance >> _finalenhance >> _xg_sampling >> _useMEforT2;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SMTopDecayer,PerturbativeDecayer>
describeHerwigSMTopDecayer("Herwig::SMTopDecayer", "HwPerturbativeDecay.so");
void SMTopDecayer::Init() {
static ClassDocumentation<SMTopDecayer> documentation
("This is the implementation of the SMTopDecayer which "
"decays top quarks into bottom quarks and either leptons "
"or quark-antiquark pairs including the matrix element for top decay",
"The matrix element correction for top decay \\cite{Hamilton:2006ms}.",
"%\\cite{Hamilton:2006ms}\n"
"\\bibitem{Hamilton:2006ms}\n"
" K.~Hamilton and P.~Richardson,\n"
" ``A simulation of QCD radiation in top quark decays,''\n"
" JHEP {\\bf 0702}, 069 (2007)\n"
" [arXiv:hep-ph/0612236].\n"
" %%CITATION = JHEPA,0702,069;%%\n");
static ParVector<SMTopDecayer,double> interfaceQuarkWeights
("QuarkWeights",
"Maximum weights for the hadronic decays",
&SMTopDecayer::_wquarkwgt, 6, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static ParVector<SMTopDecayer,double> interfaceLeptonWeights
("LeptonWeights",
"Maximum weights for the semi-leptonic decays",
&SMTopDecayer::_wleptonwgt, 3, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<SMTopDecayer,double> interfaceEnhancementFactor
("InitialEnhancementFactor",
"The enhancement factor for initial-state radiation in the shower to ensure"
" the weight for the matrix element correction is less than one.",
&SMTopDecayer::_initialenhance, 1.0, 1.0, 10000.0,
false, false, Interface::limited);
static Parameter<SMTopDecayer,double> interfaceFinalEnhancementFactor
("FinalEnhancementFactor",
"The enhancement factor for final-state radiation in the shower to ensure"
" the weight for the matrix element correction is less than one",
&SMTopDecayer::_finalenhance, 1.6, 1.0, 1000.0,
false, false, Interface::limited);
static Parameter<SMTopDecayer,double> interfaceSamplingTopHardMEC
("SamplingTopHardMEC",
"The importance sampling power for choosing an initial xg, "
"to sample xg according to xg^-_xg_sampling",
&SMTopDecayer::_xg_sampling, 1.5, 1.2, 2.0,
false, false, Interface::limited);
static Switch<SMTopDecayer,bool> interfaceUseMEForT2
("UseMEForT2",
"Use the matrix element correction, if available to fill the T2"
" region for the decay shower and don't fill using the shower",
&SMTopDecayer::_useMEforT2, true, false, false);
static SwitchOption interfaceUseMEForT2Shower
(interfaceUseMEForT2,
"Shower",
"Use the shower to fill the T2 region",
false);
static SwitchOption interfaceUseMEForT2ME
(interfaceUseMEForT2,
"ME",
"Use the Matrix element to fill the T2 region",
true);
}
double SMTopDecayer::me2(const int, const Particle & inpart,
const ParticleVector & decay,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half)));
// spinors etc for the decaying particle
if(meopt==Initialize) {
// spinors and rho
if(inpart.id()>0)
SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho,
const_ptr_cast<tPPtr>(&inpart),
incoming);
else
SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho,
const_ptr_cast<tPPtr>(&inpart),
incoming);
}
// setup spin info when needed
if(meopt==Terminate) {
// for the decaying particle
if(inpart.id()>0) {
SpinorWaveFunction::
constructSpinInfo(_inHalf,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true);
SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[1],outgoing,true);
SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[2],outgoing,true);
}
else {
SpinorBarWaveFunction::
constructSpinInfo(_inHalfBar,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true);
SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[1],outgoing,true);
SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[2],outgoing,true);
}
}
if ( ( decay[1]->momentum() + decay[2]->momentum() ).m()
< decay[1]->data().constituentMass() + decay[2]->data().constituentMass() )
return 0.0;
// spinors for the decay product
if(inpart.id()>0) {
SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar ,decay[0],outgoing);
SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[1],outgoing);
SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[2],outgoing);
}
else {
SpinorWaveFunction ::calculateWaveFunctions(_inHalf ,decay[0],outgoing);
SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[1],outgoing);
SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[2],outgoing);
}
Energy2 scale(sqr(inpart.mass()));
if(inpart.id() == ParticleID::t) {
//Define intermediate vector wave-function for Wplus
tcPDPtr Wplus(getParticleData(ParticleID::Wplus));
VectorWaveFunction inter;
unsigned int thel,bhel,fhel,afhel;
for(thel = 0;thel<2;++thel){
for(bhel = 0;bhel<2;++bhel){
inter = _wvertex->evaluate(scale,1,Wplus,_inHalf[thel],
_inHalfBar[bhel]);
for(afhel=0;afhel<2;++afhel){
for(fhel=0;fhel<2;++fhel){
(*ME())(thel,bhel,afhel,fhel) =
_wvertex->evaluate(scale,_outHalf[afhel],
_outHalfBar[fhel],inter);
}
}
}
}
}
else if(inpart.id() == ParticleID::tbar) {
VectorWaveFunction inter;
tcPDPtr Wminus(getParticleData(ParticleID::Wminus));
unsigned int tbhel,bbhel,afhel,fhel;
for(tbhel = 0;tbhel<2;++tbhel){
for(bbhel = 0;bbhel<2;++bbhel){
inter = _wvertex->
evaluate(scale,1,Wminus,_inHalf[bbhel],_inHalfBar[tbhel]);
for(afhel=0;afhel<2;++afhel){
for(fhel=0;fhel<2;++fhel){
(*ME())(tbhel,bbhel,fhel,afhel) =
_wvertex->evaluate(scale,_outHalf[afhel],
_outHalfBar[fhel],inter);
}
}
}
}
}
double output = (ME()->contract(_rho)).real();
if(abs(decay[1]->id())<=6) output *=3.;
return output;
}
void SMTopDecayer::doinit() {
PerturbativeDecayer::doinit();
//get vertices from SM object
tcHwSMPtr hwsm = dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException() << "Must have Herwig::StandardModel in "
<< "SMTopDecayer::doinit()";
_wvertex = hwsm->vertexFFW();
//initialise
_wvertex->init();
//set up decay modes
_wplus = getParticleData(ParticleID::Wplus);
DecayPhaseSpaceModePtr mode;
DecayPhaseSpaceChannelPtr Wchannel;
tPDVector extpart(4);
vector<double> wgt(1,1.0);
extpart[0] = getParticleData(ParticleID::t);
extpart[1] = getParticleData(ParticleID::b);
//lepton modes
for(int i=11; i<17;i+=2) {
extpart[2] = getParticleData(-i);
extpart[3] = getParticleData(i+1);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
Wchannel = new_ptr(DecayPhaseSpaceChannel(mode));
Wchannel->addIntermediate(extpart[0],0,0.0,-1,1);
Wchannel->addIntermediate(_wplus,0,0.0,2,3);
Wchannel->init();
mode->addChannel(Wchannel);
addMode(mode,_wleptonwgt[(i-11)/2],wgt);
}
//quark modes
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(_wvertex->allowed(-ix,iy,ParticleID::Wminus)) {
extpart[2] = getParticleData(-ix);
extpart[3] = getParticleData( iy);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
Wchannel = new_ptr(DecayPhaseSpaceChannel(mode));
Wchannel->addIntermediate(extpart[0],0,0.0,-1,1);
Wchannel->addIntermediate(_wplus,0,0.0,2,3);
Wchannel->init();
mode->addChannel(Wchannel);
addMode(mode,_wquarkwgt[iz],wgt);
++iz;
}
else {
throw InitException() << "SMTopDecayer::doinit() the W vertex"
<< "cannot handle all the quark modes"
<< Exception::abortnow;
}
}
}
}
void SMTopDecayer::dataBaseOutput(ofstream & os,bool header) const {
if(header) os << "update decayers set parameters=\"";
// parameters for the PerturbativeDecayer base class
for(unsigned int ix=0;ix<_wquarkwgt.size();++ix) {
os << "newdef " << name() << ":QuarkWeights " << ix << " "
<< _wquarkwgt[ix] << "\n";
}
for(unsigned int ix=0;ix<_wleptonwgt.size();++ix) {
os << "newdef " << name() << ":LeptonWeights " << ix << " "
<< _wleptonwgt[ix] << "\n";
}
PerturbativeDecayer::dataBaseOutput(os,false);
if(header) os << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
void SMTopDecayer::doinitrun() {
PerturbativeDecayer::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<numberModes();++ix) {
if(ix<3) _wleptonwgt[ix ] = mode(ix)->maxWeight();
else _wquarkwgt [ix-3] = mode(ix)->maxWeight();
}
}
}
WidthCalculatorBasePtr SMTopDecayer::threeBodyMEIntegrator(const DecayMode & dm) const {
// identify W decay products
int sign = dm.parent()->id() > 0 ? 1 : -1;
int iferm(0),ianti(0);
for(ParticleMSet::const_iterator pit=dm.products().begin();
pit!=dm.products().end();++pit) {
int id = (**pit).id();
if(id*sign != ParticleID::b) {
if (id*sign > 0 ) iferm = id*sign;
else ianti = id*sign;
}
}
assert(iferm!=0&&ianti!=0);
// work out which mode we are doing
int imode(-1);
for(unsigned int ix=0;ix<numberModes();++ix) {
if(mode(ix)->externalParticles(2)->id() == ianti &&
mode(ix)->externalParticles(3)->id() == iferm ) {
imode = ix;
break;
}
}
assert(imode>=0);
// get the masses we need
Energy m[3] = {mode(imode)->externalParticles(1)->mass(),
mode(imode)->externalParticles(3)->mass(),
mode(imode)->externalParticles(2)->mass()};
return
new_ptr(ThreeBodyAllOn1IntegralCalculator<SMTopDecayer>
(3,_wplus->mass(),_wplus->width(),0.0,*this,imode,m[0],m[1],m[2]));
}
InvEnergy SMTopDecayer::threeBodydGammads(const int imode, const Energy2 mt2,
const Energy2 mffb2, const Energy mb,
const Energy mf, const Energy mfb) const {
Energy mffb(sqrt(mffb2));
Energy mw(_wplus->mass());
Energy2 mw2(sqr(mw)),gw2(sqr(_wplus->width()));
Energy mt(sqrt(mt2));
Energy Eb = 0.5*(mt2-mffb2-sqr(mb))/mffb;
Energy Ef = 0.5*(mffb2-sqr(mfb)+sqr(mf))/mffb;
Energy Ebm = sqrt(sqr(Eb)-sqr(mb));
Energy Efm = sqrt(sqr(Ef)-sqr(mf));
Energy2 upp = sqr(Eb+Ef)-sqr(Ebm-Efm);
Energy2 low = sqr(Eb+Ef)-sqr(Ebm+Efm);
InvEnergy width=(dGammaIntegrand(mffb2,upp,mt,mb,mf,mfb,mw)-
dGammaIntegrand(mffb2,low,mt,mb,mf,mfb,mw))
/32./mt2/mt/8/pow(Constants::pi,3)/(sqr(mffb2-mw2)+mw2*gw2);
// couplings
width *= 0.25*sqr(4.*Constants::pi*generator()->standardModel()->alphaEM(mt2)/
generator()->standardModel()->sin2ThetaW());
width *= generator()->standardModel()->CKM(*mode(imode)->externalParticles(0),
*mode(imode)->externalParticles(1));
if(abs(mode(imode)->externalParticles(2)->id())<=6) {
width *=3.;
if(abs(mode(imode)->externalParticles(2)->id())%2==0)
width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(2),
*mode(imode)->externalParticles(3));
else
width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(3),
*mode(imode)->externalParticles(2));
}
// final spin average
assert(!std::isnan(double(width*MeV)));
return 0.5*width;
}
Energy6 SMTopDecayer::dGammaIntegrand(Energy2 mffb2, Energy2 mbf2, Energy mt,
Energy mb, Energy mf, Energy mfb, Energy mw) const {
Energy2 mt2(sqr(mt)) ,mb2(sqr(mb)) ,mf2(sqr(mf )),mfb2(sqr(mfb )),mw2(sqr(mw ));
Energy4 mt4(sqr(mt2)),mb4(sqr(mb2)),mf4(sqr(mf2)),mfb4(sqr(mfb2)),mw4(sqr(mw2));
return -mbf2 * ( + 6 * mb2 * mf2 * mfb2 * mffb2 + 6 * mb2 * mt2 * mfb2 * mffb2
+ 6 * mb2 * mt2 * mf2 * mffb2 + 12 * mb2 * mt2 * mf2 * mfb2
- 3 * mb2 * mfb4 * mffb2 + 3 * mb2 * mf2 * mffb2 * mffb2
- 3 * mb2 * mf4 * mffb2 - 6 * mb2 * mt2 * mfb4
- 6 * mb2 * mt2 * mf4 - 3 * mb4 * mfb2 * mffb2
- 3 * mb4 * mf2 * mffb2 - 6 * mb4 * mf2 * mfb2
+ 3 * mt4 * mf4 + 3 * mb4 * mfb4
+ 3 * mb4 * mf4 + 3 * mt4 * mfb4
+ 3 * mb2 * mfb2 * mffb2 * mffb2 + 3 * mt2 * mfb2 * mffb2 * mffb2
- 3 * mt2 * mfb4 * mffb2 + 3 * mt2 * mf2 * mffb2 * mffb2
- 3 * mt2 * mf4 * mffb2 - 3 * mt4 * mfb2 * mffb2
- 3 * mt4 * mf2 * mffb2 - 6 * mt4 * mf2 * mfb2
+ 6 * mt2 * mf2 * mfb2 * mffb2 + 12 * mt2 * mf2 * mw4
+ 12 * mb2 * mfb2 * mw4 + 12 * mb2 * mt2 * mw4
+ 6 * mw2 * mt2 * mfb2 * mbf2 - 12 * mw2 * mt2 * mf2 * mffb2
- 6 * mw2 * mt2 * mf2 * mbf2 - 12 * mw2 * mt2 * mf2 * mfb2
- 12 * mw2 * mb2 * mfb2 * mffb2 - 6 * mw2 * mb2 * mfb2 * mbf2
+ 6 * mw2 * mb2 * mf2 * mbf2 - 12 * mw2 * mb2 * mf2 * mfb2
- 12 * mw2 * mb2 * mt2 * mfb2 - 12 * mw2 * mb2 * mt2 * mf2
+ 12 * mf2 * mfb2 * mw4 + 4 * mbf2 * mbf2 * mw4
- 6 * mfb2 * mbf2 * mw4 - 6 * mf2 * mbf2 * mw4
- 6 * mt2 * mbf2 * mw4 - 6 * mb2 * mbf2 * mw4
+ 12 * mw2 * mt2 * mf4 + 12 * mw2 * mt4 * mf2
+ 12 * mw2 * mb2 * mfb4 + 12 * mw2 * mb4 * mfb2) /mw4 / 3.;
}
void SMTopDecayer::initializeMECorrection(RealEmissionProcessPtr born, double & initial,
double & final) {
// check the outgoing particles
PPtr part[2];
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
part[ix]= born->bornOutgoing()[ix];
}
// check the final-state particles and get the masses
if(abs(part[0]->id())==ParticleID::Wplus&&abs(part[1]->id())==ParticleID::b) {
_ma=part[0]->mass();
_mc=part[1]->mass();
}
else if(abs(part[1]->id())==ParticleID::Wplus&&abs(part[0]->id())==ParticleID::b) {
_ma=part[1]->mass();
_mc=part[0]->mass();
}
else {
return;
}
// set the top mass
_mt=born->bornIncoming()[0]->mass();
// set the gluon mass
_mg=getParticleData(ParticleID::g)->constituentMass();
// set the radiation enhancement factors
initial = _initialenhance;
final = _finalenhance;
// reduced mass parameters
_a=sqr(_ma/_mt);
_g=sqr(_mg/_mt);
_c=sqr(_mc/_mt);
double lambda = sqrt(1.+sqr(_a)+sqr(_c)-2.*_a-2.*_c-2.*_a*_c);
_ktb = 0.5*(3.-_a+_c+lambda);
_ktc = 0.5*(1.-_a+3.*_c+lambda);
useMe();
}
RealEmissionProcessPtr SMTopDecayer::applyHardMatrixElementCorrection(RealEmissionProcessPtr born) {
// Get b and a and put them in particle vector ba in that order...
ParticleVector ba;
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix)
ba.push_back(born->bornOutgoing()[ix]);
if(abs(ba[0]->id())!=5) swap(ba[0],ba[1]);
assert(born->bornIncoming().size()==1);
// Now decide if we get an emission into the dead region.
// If there is an emission newfs stores momenta for a,c,g
// according to NLO decay matrix element.
vector<Lorentz5Momentum> newfs = applyHard(ba,_ktb,_ktc);
// If there was no gluon emitted return.
if(newfs.size()!=3) return RealEmissionProcessPtr();
// Sanity checks to ensure energy greater than mass etc :)
bool check = true;
tcPDPtr gluondata=getParticleData(ParticleID::g);
if (newfs[0].e()<ba[0]->data().constituentMass()) check = false;
if (newfs[1].e()<ba[1]->mass()) check = false;
if (newfs[2].e()<gluondata->constituentMass()) check = false;
// Return if insane:
if (!check) return RealEmissionProcessPtr();
// // Set masses in 5-vectors:
newfs[0].setMass(ba[0]->mass());
newfs[1].setMass(ba[1]->mass());
newfs[2].setMass(ZERO);
// The next part of this routine sets the colour structure.
// To do this for decays we assume that the gluon comes from c!
// First create new particle objects for c, a and gluon:
PPtr newg = gluondata->produceParticle(newfs[2]);
PPtr newc = ba[0]->data().produceParticle(newfs[0]);
PPtr newa = ba[1]->data().produceParticle(newfs[1]);
born->spectator(0);
born->emitted(3);
// decaying particle
born->incoming().push_back(born->bornIncoming()[0]->dataPtr()->
produceParticle(born->bornIncoming()[0]->momentum()));
// colour flow
newg->incomingColour(born->incoming()[0],ba[0]->id()<0);
newg->colourConnect(newc ,ba[0]->id()<0);
if(born->bornOutgoing()[0]->id()==newc->id()) {
born->outgoing().push_back(newc);
born->outgoing().push_back(newa);
born->emitter(1);
}
else {
born->outgoing().push_back(newa);
born->outgoing().push_back(newc);
born->emitter(2);
}
born->outgoing().push_back(newg);
// boost for the W
LorentzRotation trans(ba[1]->momentum().findBoostToCM());
trans.boost(newfs[1].boostVector());
born->transformation(trans);
if(!inTheDeadRegion(_xg,_xa,_ktb,_ktc)) {
generator()->log()
<< "SMTopDecayer::applyHardMatrixElementCorrection()\n"
<< "Just found a point that escaped from the dead region!\n"
<< " _xg: " << _xg << " _xa: " << _xa
<< " newfs.size(): " << newfs.size() << endl;
}
born->interaction(ShowerInteraction::QCD);
return born;
}
vector<Lorentz5Momentum> SMTopDecayer::
applyHard(const ParticleVector &p,double ktb, double ktc) {
// ********************************* //
// First we see if we get a dead //
// region event: _xa,_xg //
// ********************************* //
vector<Lorentz5Momentum> fs;
// Return if there is no (NLO) gluon emission:
double weight = getHard(ktb,ktc);
if(weight>1.) {
generator()->log() << "Weight greater than 1 for hard emission in "
<< "SMTopDecayer::applyHard xg = " << _xg
<< " xa = " << _xa << "\n";
weight=1.;
}
// Accept/Reject
if (weight<UseRandom::rnd()||p.size()!= 2) return fs;
// Drop events if getHard returned a negative weight
// as in events that, somehow have escaped from the dead region
// or, worse, the allowed region.
if(weight<0.) return fs;
// Calculate xc by momentum conservation:
_xc = 2.-_xa-_xg;
// ************************************ //
// Now we get the boosts & rotations to //
// go from lab to top rest frame with //
// a in the +z direction. //
// ************************************ //
Lorentz5Momentum pa_lab,pb_lab,pc_lab,pg_lab;
// Calculate momentum of b:
pb_lab = p[0]->momentum() + p[1]->momentum();
// Define/assign momenta of c,a and the gluon:
if(abs(p[0]->id())==5) {
pc_lab = p[0]->momentum();
pa_lab = p[1]->momentum();
} else {
pc_lab = p[1]->momentum();
pa_lab = p[0]->momentum();
}
// Calculate the boost to the b rest frame:
SpinOneLorentzRotation rot0(pb_lab.findBoostToCM());
// Calculate the rotation matrix to position a along the +z direction
// in the rest frame of b and does a random rotation about z:
SpinOneLorentzRotation rot1 = rotateToZ(rot0*pa_lab);
// Calculate the boost from the b rest frame back to the lab:
// and the inverse of the random rotation about the z-axis and the
// rotation required to align a with +z:
SpinOneLorentzRotation invrot = rot0.inverse()*rot1.inverse();
// ************************************ //
// Now we construct the momenta in the //
// b rest frame using _xa,_xg. //
// First we construct b, then c and g, //
// finally we generate a by momentum //
// conservation. //
// ************************************ //
Lorentz5Momentum pa_brf, pb_brf(_mt), pc_brf, pg_brf;
// First we set the top quark to being on-shell and at rest.
// Second we set the energies of c and g,
pc_brf.setE(0.5*_mt*(2.-_xa-_xg));
pg_brf.setE(0.5*_mt*_xg);
// then their masses,
pc_brf.setMass(_mc);
pg_brf.setMass(ZERO);
// Now set the z-component of c and g. For pg we simply start from
// _xa and _xg, while for pc we assume it is equal to minus the sum
// of the z-components of a (assumed to point in the +z direction) and g.
double root=sqrt(_xa*_xa-4.*_a);
pg_brf.setZ(_mt*(1.-_xa-_xg+0.5*_xa*_xg-_c+_a)/root);
pc_brf.setZ(-1.*( pg_brf.z()+_mt*0.5*root));
// Now set the y-component of c and g's momenta
pc_brf.setY(ZERO);
pg_brf.setY(ZERO);
// Now set the x-component of c and g's momenta
pg_brf.setX(sqrt(sqr(pg_brf.t())-sqr(pg_brf.z())));
pc_brf.setX(-pg_brf.x());
// Momenta b,c,g are now set. Now we obtain a from momentum conservation,
pa_brf = pb_brf-pc_brf-pg_brf;
pa_brf.setMass(pa_brf.m());
pa_brf.rescaleEnergy();
// ************************************ //
// Now we orient the momenta and boost //
// them back to the original lab frame. //
// ************************************ //
// As in herwig6507 we assume that, in the rest frame
// of b, we have aligned the W boson momentum in the
// +Z direction by rot1*rot0*pa_lab, therefore
// we obtain the new pa_lab by applying:
// invrot*pa_brf.
pa_lab = invrot*pa_brf;
pb_lab = invrot*pb_brf;
pc_lab = invrot*pc_brf;
pg_lab = invrot*pg_brf;
fs.push_back(pc_lab);
fs.push_back(pa_lab);
fs.push_back(pg_lab);
return fs;
}
double SMTopDecayer::getHard(double ktb, double ktc) {
// zero the variables
_xg = 0.;
_xa = 0.;
_xc = 0.;
// Get a phase space point in the dead region:
double volume_factor = deadRegionxgxa(ktb,ktc);
// if outside region return -1
if(volume_factor<0) return volume_factor;
// Compute the weight for this phase space point:
double weight = volume_factor*me(_xa,_xg)*(1.+_a-_c-_xa);
// Alpha_S and colour factors - this hard wired Alpha_S needs removing.
weight *= (4./3.)/Constants::pi
*(coupling()->value(_mt*_mt*_xg*(1.-_xa+_a-_c)/(2.-_xg-_xa-_c)));
return weight;
}
bool SMTopDecayer::softMatrixElementVeto(ShowerProgenitorPtr initial,
ShowerParticlePtr parent,Branching br) {
// check if we need to apply the full correction
long id[2]={abs(initial->progenitor()->id()),abs(parent->id())};
// the initial-state correction
if(id[0]==ParticleID::t&&id[1]==ParticleID::t) {
Energy pt=br.kinematics->pT();
// check if hardest so far
// if not just need to remove effect of enhancement
bool veto(false);
// if not hardest so far
if(pt<initial->highestpT())
veto=!UseRandom::rndbool(1./_initialenhance);
// if hardest so far do calculation
else {
// values of kappa and z
double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt));
// parameters for the translation
double w(1.-(1.-z)*(kappa-1.)),u(1.+_a-_c-(1.-z)*kappa),v(sqr(u)-4.*_a*w*z);
// veto if outside phase space
if(v<0.)
veto=true;
// otherwise calculate the weight
else {
v = sqrt(v);
double xa((0.5*(u+v)/w+0.5*(u-v)/z)),xg((1.-z)*kappa);
double f(me(xa,xg)),
J(0.5*(u+v)/sqr(w)-0.5*(u-v)/sqr(z)+_a*sqr(w-z)/(v*w*z));
double wgt(f*J*2./kappa/(1.+sqr(z)-2.*z/kappa)/_initialenhance);
// This next `if' prevents the hardest emission from the
// top shower ever entering the so-called T2 region of the
// phase space if that region is to be populated by the hard MEC.
if(_useMEforT2&&xg>xgbcut(_ktb)) wgt = 0.;
if(wgt>1.) {
generator()->log() << "Violation of maximum for initial-state "
<< " soft veto in "
<< "SMTopDecayer::softMatrixElementVeto"
<< "xg = " << xg << " xa = " << xa
<< "weight = " << wgt << "\n";
wgt=1.;
}
// compute veto from weight
veto = !UseRandom::rndbool(wgt);
}
// if not vetoed reset max
if(!veto) initial->highestpT(pt);
}
// if vetoing reset the scale
if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
// return the veto
return veto;
}
// final-state correction
else if(id[0]==ParticleID::b&&id[1]==ParticleID::b) {
Energy pt=br.kinematics->pT();
// check if hardest so far
// if not just need to remove effect of enhancement
bool veto(false);
// if not hardest so far
if(pt<initial->highestpT()) return !UseRandom::rndbool(1./_finalenhance);
// if hardest so far do calculation
// values of kappa and z
double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt));
// momentum fractions
double xa(1.+_a-_c-z*(1.-z)*kappa),r(0.5*(1.+_c/(1.+_a-xa))),root(sqr(xa)-4.*_a);
if(root<0.) {
generator()->log() << "Imaginary root for final-state veto in "
<< "SMTopDecayer::softMatrixElementVeto"
<< "\nz = " << z << "\nkappa = " << kappa
<< "\nxa = " << xa
<< "\nroot^2= " << root;
parent->vetoEmission(br.type,br.kinematics->scale());
return true;
}
root=sqrt(root);
double xg((2.-xa)*(1.-r)-(z-r)*root);
// xfact (below) is supposed to equal xg/(1-z).
double xfact(z*kappa/2./(z*(1.-z)*kappa+_c)*(2.-xa-root)+root);
// calculate the full result
double f(me(xa,xg));
// jacobian
double J(z*root);
double wgt(f*J*2.*kappa/(1.+sqr(z)-2.*_c/kappa/z)/sqr(xfact)/_finalenhance);
if(wgt>1.) {
generator()->log() << "Violation of maximum for final-state soft veto in "
<< "SMTopDecayer::softMatrixElementVeto"
<< "xg = " << xg << " xa = " << xa
<< "weight = " << wgt << "\n";
wgt=1.;
}
// compute veto from weight
veto = !UseRandom::rndbool(wgt);
// if vetoing reset the scale
if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
// return the veto
return veto;
}
// otherwise don't veto
else return !UseRandom::rndbool(1./_finalenhance);
}
double SMTopDecayer::me(double xw,double xg) {
double prop(1.+_a-_c-xw),xg2(sqr(xg));
double lambda=sqrt(1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c);
double denom=(1.-2*_a*_a+_a+_c*_a+_c*_c-2.*_c);
double wgt=-_c*xg2/prop+(1.-_a+_c)*xg-(prop*(1 - xg)+xg2)
+(0.5*(1.+2.*_a+_c)*sqr(prop-xg)*xg+2.*_a*prop*xg2)/denom;
return wgt/(lambda*prop);
}
// This function is auxiliary to the xab function.
double SMTopDecayer::xgbr(int toggle) {
return 1.+toggle*sqrt(_a)-_c*(1.-toggle*sqrt(_a))/(1.-_a);
}
// This function is auxiliary to the xab function.
double SMTopDecayer::ktr(double xgb, int toggle) {
return 2.*xgb/
(xgb+toggle*sqrt((1.-1./_a)
*(xgb-xgbr( 1))
*(xgb-xgbr(-1))));
}
// Function xab determines xa (2*W energy fraction) for a given value
// of xg (2*gluon energy fraction) and kappa tilde (q tilde squared over
// m_top squared). Hence this function allows you to draw 1: the total
// phase space volume in the xa vs xg plane 2: for a given value of
// kappa tilde (i.e. starting evolution scale) the associated contour
// in the xa vs xg plane (and hence the regions that either shower can
// populate). This calculation is done assuming the emission came from
// the top quark i.e. kappa tilde here is the q tilde squared of the TOP
// quark divided by m_top squared.
double SMTopDecayer::xab(double xgb, double kt, int toggle) {
double xab;
if(toggle==2) {
// This applies for g==0.&&kt==ktr(a,c,0.,xgb,1).
xab = -2.*_a*(xgb-2.)/(1.+_a-_c-xgb);
} else if(toggle==1) {
// This applies for kt==1&&g==0.
double lambda = sqrt(sqr(xgb-1.+_a+_c)-4.*_a*_c);
xab = (0.5/(kt-xgb))*(kt*(1.+_a-_c-xgb)-lambda)
+ (0.5/(kt+xgb*(1.-kt)))*(kt*(1.+_a-_c-xgb)+lambda);
} else {
// This is the form of xab FOR _g=0.
double ktmktrpktmktrm = kt*kt - 4.*_a*(kt-1.)*xgb*xgb
/ (sqr(1.-_a-_c-xgb)-4.*_a*_c);
if(fabs(kt-(2.*xgb-2.*_g)/(xgb-sqrt(xgb*xgb-4.*_g)))/kt>1.e-6) {
double lambda = sqrt((sqr(1.-_a-_c-xgb)-4.*_a*_c)*ktmktrpktmktrm);
xab = (0.5/(kt-xgb))*(kt*(1.+_a-_c-xgb)-lambda)
+ (0.5/(kt+xgb*(1.-kt)))*(kt*(1.+_a-_c-xgb)+lambda);
}
else {
// This is the value of xa as a function of xb when kt->infinity.
// Where we take any kt > (2.*xgb-2.*_g)/(xgb-sqrt(xgb*xgb-4.*_g))
// as being effectively infinite. This kt value is actually the
// maximum allowed value kt can have if the phase space is calculated
// without the approximation of _g=0 (massless gluon). This formula
// for xab below is then valid for _g=0 AND kt=infinity only.
xab = ( 2.*_c+_a*(xgb-2.)
+ 3.*xgb
- xgb*(_c+xgb+sqrt(_a*_a-2.*(_c-xgb+1.)*_a+sqr(_c+xgb-1.)))
- 2.
)/2./(xgb-1.);
}
}
if(std::isnan(xab)) {
double ktmktrpktmktrm = ( sqr(xgb*kt-2.*(xgb-_g))
-kt*kt*(1.-1./_a)*(xgb-xgbr( 1)-_g/(1.+sqrt(_a)))
*(xgb-xgbr(-1)-_g/(1.-sqrt(_a)))
)/
(xgb*xgb-(1.-1./_a)*(xgb-xgbr( 1)-_g/(1.+sqrt(_a)))
*(xgb-xgbr(-1)-_g/(1.-sqrt(_a)))
);
double lambda = sqrt((xgb-1.+sqr(sqrt(_a)+sqrt(_c-_g)))
*(xgb-1.+sqr(sqrt(_a)-sqrt(_c-_g)))*
ktmktrpktmktrm);
xab = (0.5/(kt-xgb+_g))*(kt*(1.+_a-_c+_g-xgb)-lambda)
+ (0.5/(kt+xgb*(1.-kt)-_g))*(kt*(1.+_a-_c+_g-xgb)+lambda);
if(std::isnan(xab))
throw Exception() << "TopMECorrection::xab complex x_a value.\n"
<< " xgb = " << xgb << "\n"
<< " xab = " << xab << "\n"
<< " toggle = " << toggle << "\n"
<< " ktmktrpktmktrm = " << ktmktrpktmktrm
<< Exception::eventerror;
}
return xab;
}
// xgbcut is the point along the xg axis where the upper bound on the
// top quark (i.e. b) emission phase space goes back on itself in the
// xa vs xg plane i.e. roughly mid-way along the xg axis in
// the xa vs xg Dalitz plot.
double SMTopDecayer::xgbcut(double kt) {
double lambda2 = 1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c;
double num1 = kt*kt*(1.-_a-_c);
double num2 = 2.*kt*sqrt(_a*(kt*kt*_c+lambda2*(kt-1.)));
return (num1-num2)/(kt*kt-4.*_a*(kt-1.));
}
double SMTopDecayer::xaccut(double kt) {
return 1.+_a-_c-0.25*kt;
}
double SMTopDecayer::z(double xac, double kt,
int toggle1, int toggle2) {
double z = -1.0;
if(toggle2==0) {
z = (kt+toggle1*sqrt(kt*(kt-4.*(1.+_a-_c-xac))))/(2.*kt);
} else if(toggle2==1) {
z = ((1.+_a+_c-xac)+toggle1*(1.+_a-_c-xac))
/(2.*(1.+_a-xac));
} else if(toggle2==2) {
z = 0.5;
} else {
throw Exception() << "Cannot determine z in SMTopDecayer::z()"
<< Exception::eventerror;
}
return z;
}
double SMTopDecayer::xgc(double xac, double kt,
int toggle1, int toggle2) {
double tiny(1.e-6);
double xaToMinBoundary(xac*xac-4.*_a);
if(xaToMinBoundary<0) {
if(fabs(xaToMinBoundary/(1.-_a)/(1.-_a))<tiny)
xaToMinBoundary *= -1.;
else
throw Exception() << "SMTopDecayer::xgc xa not in phase space!"
<< Exception::eventerror;
}
return (2.-xac)*(1.-0.5*(1.+_c/(1.+_a-xac)))
-(z(xac,kt,toggle1,toggle2)-0.5*(1.+_c/(1.+_a-xac)))
*sqrt(xaToMinBoundary);
}
double SMTopDecayer::xginvc0(double xg , double kt) {
// The function xg(kappa_tilde_c,xa) surely, enough, draws a
// line of constant kappa_tilde_c in the xg, xa Dalitz plot.
// Such a function can therefore draw the upper and lower
// edges of the phase space for emission from c (the b-quark).
// However, to sample the soft part of the dead zone effectively
// we want to generate a value of xg first and THEN distribute
// xa in the associated allowed part of the dead zone. Hence, the
// function we want, to define the dead zone in xa for a given
// xg, is the inverse of xg(kappa_tilde_c,xa). The full expression
// for xg(kappa_tilde_c,xa) is complicated and, sure enough,
// does not invert. Therefore we try to overestimate the size
// of the dead zone initially, rejecting events which do not
// fall exactly inside it afterwards, with the immediate aim
// of getting an approximate version of xg(kappa_tilde_c,xa)
// that can be inverted. We do this by simply setting c=0 i.e.
// the b-quark mass to zero (and the gluon mass of course), in
// the full expression xg(...). The result of inverting this
// function is the output of this routine (a value of xa) hence
// the name xginvc0. xginvc0 is calculated to be,
// xginvc0 = (1./3.)*(1.+a+pow((U+sqrt(4.*V*V*V+U*U))/2.,1./3.)
// -V*pow(2./(U+sqrt(4.*V*V*V+U*U)),1./3.)
// )
// U = 2.*a*a*a - 66.*a*a + 9.*a*kt*xg + 18.*a*kt
// - 66.*a + 27.*kt*xg*xg - 45.*kt*xg +18.*kt +2. ;
// V = -1.-a*a-14.*a-3.kt*xg+3.*kt;
// This function, as with many functions in this ME correction,
// is plagued by cuts that have to handled carefully in numerical
// implementation. We have analysed the cuts and hence we implement
// it in the following way, with a series of 'if' statements.
//
// A useful -definition- to know in deriving the v<0 terms is
// that tanh^-1(z) = 0.5*(log(1.+z)-log(1.-z)).
double u,v,output;
u = 2.*_a*_a*_a-66.*_a*_a
+9.*xg*kt*_a+18.*kt*_a
-66.*_a+27.*xg*xg*kt
-45.*xg*kt+18.*kt+2.;
v = -_a*_a-14.*_a-3.*xg*kt+3.*kt-1.;
double u2=u*u,v3=v*v*v;
if(v<0.) {
if(u>0.&&(4.*v3+u2)<0.) output = cos( atan(sqrt(-4.*v3-u2)/u)/3.);
else if(u>0.&&(4.*v3+u2)>0.) output = cosh(atanh(sqrt( 4.*v3+u2)/u)/3.);
else output = cos(( atan(sqrt(-4.*v3-u2)/u)
+Constants::pi)/3.);
output *= 2.*sqrt(-v);
} else {
output = sinh(log((u+sqrt(4.*v3+u2))/(2.*sqrt(v3)))/3.);
output *= 2.*sqrt(v);
}
if(!isfinite(output)) {
throw Exception() << "TopMECorrection::xginvc0:\n"
<< "possible numerical instability detected.\n"
<< "\n v = " << v << " u = " << u << "\n4.*v3+u2 = " << 4.*v3+u2
<< "\n_a = " << _a << " ma = " << sqrt(_a*_mt*_mt/GeV2)
<< "\n_c = " << _c << " mc = " << sqrt(_c*_mt*_mt/GeV2)
<< "\n_g = " << _g << " mg = " << sqrt(_g*_mt*_mt/GeV2)
<< Exception::eventerror;
}
return ( 1.+_a +output)/3.;
}
double SMTopDecayer::approxDeadMaxxa(double xg,double ktb,double ktc) {
double maxxa(0.);
double x = min(xginvc0(xg,ktc),
xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0));
double y(-9999999999.);
if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) {
y = max(xab(xg,ktb,0),xab(xg,1.,1));
} else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1));
}
if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
if(x>=y) { maxxa = x ; }
else { maxxa = -9999999.; }
} else {
maxxa = -9999999.;
}
return maxxa;
}
double SMTopDecayer::approxDeadMinxa(double xg,double ktb,double ktc) {
double minxa(0.);
double x = min(xginvc0(xg,ktc),
xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0));
double y(-9999999999.);
if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) {
y = max(xab(xg,ktb,0),xab(xg,1.,1));
} else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
if(_useMEforT2) y = xab(xg,1.,1);
else y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1));
}
if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
if(x>=y) { minxa = y ; }
else { minxa = 9999999.; }
} else {
minxa = 9999999.;
}
return minxa;
}
// This function returns true if the phase space point (xg,xa) is in the
// kinematically allowed phase space.
bool SMTopDecayer::inTheAllowedRegion(double xg , double xa) {
bool output(true);
if(xg<2.*sqrt(_g)||xg>1.-sqr(sqrt(_a)+sqrt(_c))) output = false;
if(xa<xab(xg,1.,1)) output = false;
if(xa>xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0)) output = false;
return output;
}
// This function returns true if the phase space point (xg,xa) is in the
// approximate (overestimated) dead region.
bool SMTopDecayer::inTheApproxDeadRegion(double xg , double xa,
double ktb, double ktc) {
bool output(true);
if(!inTheAllowedRegion(xg,xa)) output = false;
if(xa<approxDeadMinxa(xg,ktb,ktc)) output = false;
if(xa>approxDeadMaxxa(xg,ktb,ktc)) output = false;
return output;
}
// This function returns true if the phase space point (xg,xa) is in the
// dead region.
bool SMTopDecayer::inTheDeadRegion(double xg , double xa,
double ktb, double ktc) {
bool output(true);
if(!inTheApproxDeadRegion(xg,xa,ktb,ktc)) output = false;
if(xa>xaccut(ktc)) {
if(xg<xgc(max(xaccut(ktc),2.*sqrt(_a)),ktc, 1,2)&&
xg>xgc(xa,ktc, 1,0)) { output = false; }
if(xg>xgc(max(xaccut(ktc),2.*sqrt(_a)),ktc,-1,2)&&
xg<xgc(xa,ktc,-1,0)) { output = false; }
}
return output;
}
// This function attempts to generate a phase space point in the dead
// region and returns the associated phase space volume factor needed for
// the associated event weight.
double SMTopDecayer::deadRegionxgxa(double ktb,double ktc) {
_xg=0.;
_xa=0.;
// Here we set limits on xg and generate a value inside the bounds.
double xgmin(2.*sqrt(_g)),xgmax(1.-sqr(sqrt(_a)+sqrt(_c)));
// Generate _xg.
if(_xg_sampling==2.) {
_xg=xgmin*xgmax/(xgmin+UseRandom::rnd()*(xgmax-xgmin));
} else {
_xg=xgmin*xgmax/pow(( pow(xgmin,_xg_sampling-1.)
+ UseRandom::rnd()*(pow(xgmax,_xg_sampling-1.)
-pow(xgmin,_xg_sampling-1.))
),1./(_xg_sampling-1.));
}
// Here we set the bounds on _xa for given _xg.
if(_xg<xgmin||xgmin>xgmax)
throw Exception() << "TopMECorrection::deadRegionxgxa:\n"
<< "upper xg bound is less than the lower xg bound.\n"
<< "\n_xg = " << _xg
<< "\n2.*sqrt(_g) = " << 2.*sqrt(_g)
<< "\n_a = " << _a << " ma = " << sqrt(_a*_mt*_mt/GeV2)
<< "\n_c = " << _c << " mc = " << sqrt(_c*_mt*_mt/GeV2)
<< "\n_g = " << _g << " mg = " << sqrt(_g*_mt*_mt/GeV2)
<< Exception::eventerror;
double xamin(approxDeadMinxa(_xg,ktb,ktc));
double xamax(approxDeadMaxxa(_xg,ktb,ktc));
// Are the bounds sensible? If not return.
if(xamax<=xamin) return -1.;
_xa=1.+_a-(1.+_a-xamax)*pow((1.+_a-xamin)/(1.+_a-xamax),UseRandom::rnd());
// If outside the allowed region return -1.
if(!inTheDeadRegion(_xg,_xa,ktb,ktc)) return -1.;
// The integration volume for the weight
double xg_vol,xa_vol;
if(_xg_sampling==2.) {
xg_vol = (xgmax-xgmin)
/ (xgmax*xgmin);
} else {
xg_vol = (pow(xgmax,_xg_sampling-1.)-pow(xgmin,_xg_sampling-1.))
/ ((_xg_sampling-1.)*pow(xgmax*xgmin,_xg_sampling-1.));
}
xa_vol = log((1.+_a-xamin)/(1.+_a-xamax));
// Here we return the integral volume factor multiplied by the part of the
// weight left over which is not included in the BRACES function, i.e.
// the part of _xg^-2 which is not absorbed in the integration measure.
return xg_vol*xa_vol*pow(_xg,_xg_sampling-2.);
}
LorentzRotation SMTopDecayer::rotateToZ(Lorentz5Momentum v) {
// compute the rotation matrix
LorentzRotation trans;
// rotate so in z-y plane
trans.rotateZ(-atan2(v.y(),v.x()));
// rotate so along Z
trans.rotateY(-acos(v.z()/v.vect().mag()));
// generate random rotation
double c,s,cs;
do
{
c = 2.*UseRandom::rnd()-1.;
s = 2.*UseRandom::rnd()-1.;
cs = c*c+s*s;
}
while(cs>1.||cs==0.);
double cost=(c*c-s*s)/cs,sint=2.*c*s/cs;
// apply random azimuthal rotation
trans.rotateZ(atan2(sint,cost));
return trans;
}
bool SMTopDecayer::deadZoneCheck(double xw, double xg){
//veto events not in the dead cone
double Lambda = sqrt(1. + sqr(s2()) + sqr(e2()) - 2.*s2() - 2.*e2() - 2.*s2()*e2());
double kappa = e2() + 0.5*(1. - s2() + e2() + Lambda);
//invert xw for z values
double A = 1.;
double B = -1.;
double C = (1.+s2()-e2()-xw)/kappa;
if((sqr(B) - 4.*A*C) >= 0.){
double z[2];
z[0] = (-B + sqrt(sqr(B) - 4.*A*C))/(2.*A);
z[1] = (-B - sqrt(sqr(B) - 4.*A*C))/(2.*A);
double r = 0.5*(1. + e2()/(1. + s2()- xw));
double xg_lims [2];
xg_lims[0] = (2. - xw)*(1.-r) - (z[0]-r)*sqrt(sqr(xw) - 4.*s2());
xg_lims[1] = (2. - xw)*(1.-r) - (z[1]-r)*sqrt(sqr(xw) - 4.*s2());
double xg_low_lim = min(xg_lims[0], xg_lims[1]);
double xg_upp_lim = max(xg_lims[0], xg_lims[1]);
if (xg>=xg_low_lim && xg<=xg_upp_lim) return false;
}
double kappa_t = 1. + 0.5*(1. - s2() + e2() + Lambda);
double z = 1. - xg/kappa_t;
double u = 1. + s2() - e2() - (1.-z)*kappa_t;
double y = 1. - (1.-z)*(kappa_t-1.);
if (sqr(u) - 4.*s2()*y*z >= 0.){
double v = sqrt(sqr(u) - 4.*s2()*y*z);
double xw_lim = (u + v) / (2.*y) + (u - v) / (2.*z);
if (xw <= xw_lim) return false;
}
else if (sqr(u) - 4.*s2()*y*z < 0.){
double xg_lim = (8.*s2() -2.*xw*(1-e2()+s2()))/(4.*s2()-2.*xw);
if (xg>=xg_lim) return false;
}
return true;
}
double SMTopDecayer::matrixElementRatio(const Particle & inpart,
const ParticleVector & ,
const ParticleVector & decay3,
MEOption ) {
double f = (1. + sqr(e2()) - 2.*sqr(s2()) + s2() + s2()*e2() - 2.*e2());
double Nc = standardModel()->Nc();
double Cf = (sqr(Nc) - 1.) / (2.*Nc);
double B = f/s2();
Energy2 PbPg = decay3[0]->momentum()*decay3[2]->momentum();
Energy2 PtPg = inpart.momentum()*decay3[2]->momentum();
Energy2 PtPb = inpart.momentum()*decay3[0]->momentum();
double R = Cf *((-4.*sqr(mb())*f/s2()) * ((sqr(mb())*e2()/sqr(PbPg)) +
(sqr(mb())/sqr(PtPg)) - 2.*(PtPb/(PtPg*PbPg))) +
(16. + 8./s2() + 8.*e2()/s2()) * ((PtPg/PbPg) + (PbPg/PtPg)) -
(16./s2()) * (1. + e2()));
- return R/B;
+ return R/B*Constants::pi;
}
diff --git a/Decay/Perturbative/SMWDecayer.cc b/Decay/Perturbative/SMWDecayer.cc
--- a/Decay/Perturbative/SMWDecayer.cc
+++ b/Decay/Perturbative/SMWDecayer.cc
@@ -1,924 +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.),
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_ << NLO_;
}
void SMWDecayer::persistentInput(PersistentIStream & is, int) {
is >> FFWVertex_ >> quarkWeight_ >> leptonWeight_
>> 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 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());
+ double ratio = CF_*reme/lome*sqr(inpart.mass())*4.*Constants::pi;
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/SMZDecayer.cc b/Decay/Perturbative/SMZDecayer.cc
--- a/Decay/Perturbative/SMZDecayer.cc
+++ b/Decay/Perturbative/SMZDecayer.cc
@@ -1,1331 +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.),
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_ << NLO_
<< gluon_;
}
void SMZDecayer::persistentInput(PersistentIStream & is, int) {
is >> FFZVertex_ >> FFPVertex_ >> FFGVertex_
>> 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 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);
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,
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*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;
}
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());
}
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());
+ double ratio = CF_*reme/lome*sqr(inpart.mass())*4.*Constants::pi;
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;
}
double SMZDecayer::calculateRealEmission(double x1, double x2,
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 {
// 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());
// 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);
// total real emission contribution
return realwgt;
}
diff --git a/Decay/PerturbativeDecayer.h b/Decay/PerturbativeDecayer.h
--- a/Decay/PerturbativeDecayer.h
+++ b/Decay/PerturbativeDecayer.h
@@ -1,245 +1,245 @@
// -*- C++ -*-
#ifndef Herwig_PerturbativeDecayer_H
#define Herwig_PerturbativeDecayer_H
//
// This is the declaration of the PerturbativeDecayer class.
//
#include "Herwig/Decay/DecayIntegrator.h"
#include "Herwig/Shower/Core/Couplings/ShowerAlpha.h"
namespace Herwig {
using namespace ThePEG;
/**
* The PerturbativeDecayer class is the base class for perturbative decays in
* Herwig and implements the functuality for the POWHEG corrections
*
* @see \ref PerturbativeDecayerInterfaces "The interfaces"
* defined for PerturbativeDecayer.
*/
class PerturbativeDecayer: public DecayIntegrator {
protected:
/**
* Type of dipole
*/
enum dipoleType {FFa, FFc, IFa, IFc, IFba, IFbc};
public:
/**
* The default constructor.
*/
PerturbativeDecayer() : pTmin_(GeV), pT_(ZERO), mb_(ZERO),
e_(0.), s_(0.), e2_(0.), s2_(0.)
{}
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {return No;}
/**
* Member to generate the hardest emission in the POWHEG scheme
*/
virtual RealEmissionProcessPtr generateHardest(RealEmissionProcessPtr);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Three-body matrix element including additional QCD radiation
*/
virtual double threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay, MEOption meopt);
/**
- * Calculate matrix element ratio R/B
+ * Calculate matrix element ratio \f$\frac{M^2}{\alpha_S}\frac{|\overline{\rm{ME}}_3|}{|\overline{\rm{ME}}_2|}\f$
*/
virtual double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
const ParticleVector & decay3, MEOption meopt);
/**
* Work out the type of process
*/
bool identifyDipoles(vector<dipoleType> & dipoles,
PPtr & aProgenitor,
PPtr & bProgenitor,
PPtr & cProgenitor) const;
/**
* Coupling for the generation of hard radiation
*/
ShowerAlphaPtr coupling() {return coupling_;}
/**
* Return the momenta including the hard emission
*/
vector<Lorentz5Momentum> hardMomenta(PPtr in, PPtr emitter,
PPtr spectator,
const vector<dipoleType> & dipoles, int i);
/**
* Calculate momenta of all the particles
*/
bool calcMomenta(int j, Energy pT, double y, double phi, double& xg,
double& xs, double& xe, double& xe_z,
vector<Lorentz5Momentum>& particleMomenta);
/**
* Check the calculated momenta are physical
*/
bool psCheck(const double xg, const double xs);
/**
* Return dipole corresponding to the dipoleType dipoleId
*/
InvEnergy2 calculateDipole(const dipoleType & dipoleId, const Particle & inpart,
const ParticleVector & decay3, const dipoleType & emittingDipole);
/**
* Return contribution to dipole that depends on the spin of the emitter
*/
double dipoleSpinFactor(const PPtr & emitter, double z);
/**
* Return the colour coefficient of the dipole
*/
double colourCoeff(const PDT::Colour emitter, const PDT::Colour spectator,
const PDT::Colour other);
/**
* Set up the colour lines
*/
void getColourLines(RealEmissionProcessPtr real);
protected:
/**
* Access to the kinematics for inheriting classes
*/
//@{
/**
* Transverse momentum of the emission
*/
const Energy & pT() const { return pT_;}
/**
* Mass of decaying particle
*/
const Energy & mb() const {return mb_;}
/**
* Reduced mass of emitter child particle
*/
const double & e() const {return e_;}
/**
* Reduced mass of spectator child particle
*/
const double & s() const {return s_;}
/**
* Reduced mass of emitter child particle squared
*/
const double & e2() const {return e2_;}
/**
* Reduced mass of spectator child particle squared
*/
const double & s2() const {return s2_;}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
PerturbativeDecayer & operator=(const PerturbativeDecayer &);
private:
/**
* Members for the generation of the hard radiation
*/
//@{
/**
* Coupling for the generation of hard radiation
*/
ShowerAlphaPtr coupling_;
/**
* Minimum \f$p_T\f$
*/
Energy pTmin_;
//@}
private:
/**
* Mmeber variables for the kinematics of the hard emission
*/
//@{
/**
* Transverse momentum of the emission
*/
Energy pT_;
/**
* Mass of decaying particle
*/
Energy mb_;
/**
* Reduced mass of emitter child particle
*/
double e_;
/**
* Reduced mass of spectator child particle
*/
double s_;
/**
* Reduced mass of emitter child particle squared
*/
double e2_;
/**
* Reduced mass of spectator child particle squared
*/
double s2_;
//@}
};
}
#endif /* Herwig_PerturbativeDecayer_H */