Page MenuHomeHEPForge
Files F7878364

No OneTemporary
Actions

View Options

diff --git a/Decay/Perturbative/SMHiggsGGHiggsPPDecayer.cc b/Decay/Perturbative/SMHiggsGGHiggsPPDecayer.cc
--- a/Decay/Perturbative/SMHiggsGGHiggsPPDecayer.cc
+++ b/Decay/Perturbative/SMHiggsGGHiggsPPDecayer.cc
@@ -1,423 +1,424 @@
// -*- C++ -*-
//
// SMHiggsGGHiggsPPDecayer.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 SMHiggsGGHiggsPPDecayer class.
//
#include "SMHiggsGGHiggsPPDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include "Herwig/Utilities/HiggsLoopFunctions.h"
using namespace Herwig;
using namespace Herwig::HiggsLoopFunctions;
using namespace ThePEG::Helicity;
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SMHiggsGGHiggsPPDecayer,PerturbativeDecayer>
describeHerwigSMHiggsGGHiggsPPDecayer("Herwig::SMHiggsGGHiggsPPDecayer",
"HwPerturbativeHiggsDecay.so");
bool SMHiggsGGHiggsPPDecayer::accept(tcPDPtr parent,
const tPDVector & children) const {
int idp = parent->id();
int id0 = children[0]->id();
int id1 = children[1]->id();
if((idp == ParticleID::h0 &&
id0 == ParticleID::g && id1 == ParticleID::g) ||
(idp == ParticleID::h0 &&
id0 == ParticleID::gamma && id1 == ParticleID::gamma)||
(idp == ParticleID::h0 &&
id0 == ParticleID::Z0 && id1 == ParticleID::gamma)||
(idp == ParticleID::h0 &&
id0 == ParticleID::gamma && id1 == ParticleID::Z0))
return true;
else
return false;
}
ParticleVector SMHiggsGGHiggsPPDecayer::decay(const Particle & parent,
const tPDVector & children) const {
int imode(2);
if(children[0]->id() == ParticleID::gamma &&
children[1]->id() == ParticleID::gamma)
imode = 1;
else if(children[0]->id() ==ParticleID::g)
imode = 0;
ParticleVector out(generate(true,false,imode,parent));
//colour flow
if(children[0]->id() == ParticleID::g &&
children[1]->id() == ParticleID::g) {
out[0]->colourNeighbour(out[1]);
out[0]->antiColourNeighbour(out[1]);
}
return out;
}
void SMHiggsGGHiggsPPDecayer::persistentOutput(PersistentOStream & os) const {
os << _hggvertex << _hppvertex << _hzpvertex << _h0wgt
<< _minloop << _maxloop << _massopt;
}
void SMHiggsGGHiggsPPDecayer::persistentInput(PersistentIStream & is, int) {
is >> _hggvertex >> _hppvertex >> _hzpvertex >> _h0wgt
>> _minloop >> _maxloop >> _massopt;
}
void SMHiggsGGHiggsPPDecayer::Init() {
static ClassDocumentation<SMHiggsGGHiggsPPDecayer> documentation
("This is an implentation of h0->gg or h0->gamma,gamma "
"decayer using the SMHGGVertex.");
static Reference<SMHiggsGGHiggsPPDecayer,AbstractVVSVertex>
interfaceSMHGGVertex
("SMHGGVertex",
"Pointer to SMHGGVertex",
&SMHiggsGGHiggsPPDecayer::_hggvertex, false, false, true,
false, false);
static Reference<SMHiggsGGHiggsPPDecayer,AbstractVVSVertex>
interfaceSMHPPVertex
("SMHPPVertex",
"Pointer to SMHPPVertex",
&SMHiggsGGHiggsPPDecayer::_hppvertex, false, false, true,
false, false);
static Reference<SMHiggsGGHiggsPPDecayer,AbstractVVSVertex>
interfaceSMHZPVertex
("SMHZPVertex",
"Pointer to SMHZPVertex",
&SMHiggsGGHiggsPPDecayer::_hzpvertex, false, false, true,
false, false);
static ParVector<SMHiggsGGHiggsPPDecayer,double> interfaceMaxWeights
("MaxWeights",
"Maximum weights for the various decays",
&SMHiggsGGHiggsPPDecayer::_h0wgt, 3, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<SMHiggsGGHiggsPPDecayer,int> interfaceMinimumInLoop
("MinimumInLoop",
"The minimum flavour of the quarks to include in the loops",
&SMHiggsGGHiggsPPDecayer::_minloop, 6, 4, 6,
false, false, Interface::limited);
static Parameter<SMHiggsGGHiggsPPDecayer,int> interfaceMaximumInLoop
("MaximumInLoop",
"The maximum flavour of the quarks to include in the loops",
&SMHiggsGGHiggsPPDecayer::_maxloop, 6, 4, 6,
false, false, Interface::limited);
static Switch<SMHiggsGGHiggsPPDecayer,unsigned int> interfaceMassOption
("MassOption",
"Option for the treatment of the masses in the loop diagrams",
&SMHiggsGGHiggsPPDecayer::_massopt, 0, false, false);
static SwitchOption interfaceMassOptionFull
(interfaceMassOption,
"Full",
"Include the full mass dependence",
0);
static SwitchOption interfaceMassOptionLarge
(interfaceMassOption,
"Large",
"Use the heavy mass limit",
1);
}
double SMHiggsGGHiggsPPDecayer::me2(const int,
const Particle & part,
const ParticleVector & decay,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1,PDT::Spin1)));
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);
for(unsigned int ix=0;ix<2;++ix)
VectorWaveFunction::constructSpinInfo(_vwave[ix],decay[ix],
outgoing,true,
decay[ix]->id()!=ParticleID::Z0);
return 0.;
}
for(unsigned int ix=0;ix<2;++ix)
VectorWaveFunction::
calculateWaveFunctions(_vwave[ix],decay[ix],outgoing,decay[ix]->id()!=ParticleID::Z0);
//Set up decay matrix
Energy2 scale(sqr(part.mass()));
unsigned int v1hel,v2hel;
AbstractVVSVertexPtr vertex;
unsigned int vstep1(2),vstep2(2);
double sym(1.);
if(decay[0]->id() == ParticleID::g &&
decay[1]->id() == ParticleID::g) {
vertex = _hggvertex;
sym = 2.;
}
else if(decay[0]->id() == ParticleID::gamma &&
decay[1]->id() == ParticleID::gamma) {
vertex = _hppvertex;
sym = 2.;
}
else if(decay[0]->id() == ParticleID::Z0 &&
decay[1]->id() == ParticleID::gamma) {
vertex = _hzpvertex;
vstep1 = 1;
}
else if(decay[1]->id() == ParticleID::Z0 &&
decay[0]->id() == ParticleID::gamma) {
vertex = _hzpvertex;
vstep2 = 1;
}
else
assert(false);
// loop over the helicities of the outgoing bosons
for(v1hel = 0;v1hel < 3;v1hel+=vstep1) {
for(v2hel = 0;v2hel < 3;v2hel+=vstep2) {
(*ME())(0,v1hel,v2hel) = vertex->evaluate(scale,_vwave[0][v1hel],
_vwave[1][v2hel],_swave);
}
}
//store matrix element
double output = ME()->contract(_rho).real()*UnitRemoval::E2/scale;
//colour factor (N^2 - 1)/4
if(decay[0]->id() == ParticleID::g) output *= 8.;
//symmetric final states
output /= sym;
// return the answer
return output;
}
void SMHiggsGGHiggsPPDecayer::doinit() {
PerturbativeDecayer::doinit();
// get the width generator for the higgs
tPDPtr higgs = getParticleData(ParticleID::h0);
if(_hggvertex) {
_hggvertex->init();
}
else {
throw InitException() << "SMHiggsGGHiggsPPDecayer::doinit() - "
<< "_hggvertex is null";
}
if(_hppvertex) {
_hppvertex->init();
}
else {
throw InitException() << "SMHiggsGGHiggsPPDecayer::doinit() - "
<< "_hppvertex is null";
}
if(_hzpvertex) {
_hzpvertex->init();
}
else {
throw InitException() << "SMHiggsGGHiggsZPDecayer::doinit() - "
<< "_hzpvertex is null";
}
//set up decay modes
DecayPhaseSpaceModePtr mode;
tPDVector extpart(3);
vector<double> wgt(0);
// glu,glu mode
extpart[0] = getParticleData(ParticleID::h0);
extpart[1] = getParticleData(ParticleID::g);
extpart[2] = getParticleData(ParticleID::g);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
addMode(mode,_h0wgt[0],wgt);
// gamma,gamma mode
extpart[1] = getParticleData(ParticleID::gamma);
extpart[2] = getParticleData(ParticleID::gamma);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
addMode(mode,_h0wgt[1],wgt);
// Z0,gamma mode
extpart[1] = getParticleData(ParticleID::Z0);
extpart[2] = getParticleData(ParticleID::gamma);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
addMode(mode,_h0wgt[2],wgt);
}
void SMHiggsGGHiggsPPDecayer::doinitrun() {
_hggvertex->initrun();
_hppvertex->initrun();
_hzpvertex->initrun();
PerturbativeDecayer::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<numberModes();++ix) {
_h0wgt[ix] = mode(ix)->maxWeight();
}
}
}
void SMHiggsGGHiggsPPDecayer::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<_h0wgt.size();++ix) {
os << "newdef " << name() << ":MaxWeights " << ix << " "
<< _h0wgt[ix] << "\n";
}
os << "newdef " << name() << ":SMHGGVertex " << _hggvertex->fullName() << "\n";
os << "newdef " << name() << ":SMHPPVertex " << _hppvertex->fullName() << "\n";
os << "newdef " << name() << ":SMHZPVertex " << _hzpvertex->fullName() << "\n";
PerturbativeDecayer::dataBaseOutput(os,false);
if(header) os << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
double SMHiggsGGHiggsPPDecayer::matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
const ParticleVector & decay3, MEOption,
ShowerInteraction inter) {
assert(inter==ShowerInteraction::QCD);
// 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());
}
Energy2 scale = sqr(inpart.mass());
Energy2 lome = loME(inpart.mass());
double reme = realME(partons,realmom);
double ratio = reme/lome*scale;
// // analytic value for mt -> infinity
// double x1 = 2.*decay3[0]->momentum().t()/inpart.mass();
// double x2 = 2.*decay3[1]->momentum().t()/inpart.mass();
// double x3 = 2.*decay3[2]->momentum().t()/inpart.mass();
// double test = 8.*Constants::pi*3.*(1.+pow(1-x1,4)+pow(1-x2,4)+pow(1-x3,4))
// /(1.-x1)/(1.-x2)/(1.-x3);
// generator()->log() << "TESTING RATIO " << test << " " << ratio << " " << ratio/test << "\n";
- return ratio;
+ // remember the symmetry factor
+ return ratio/3.;
}
double SMHiggsGGHiggsPPDecayer::realME(//const vector<cPDPtr> & partons,
const vector<cPDPtr> &,
const vector<Lorentz5Momentum> & momenta) const {
// using std::norm;
// ScalarWaveFunction hout(momenta[0],partons[0],outgoing);
// LorentzPolarizationVector g[3][2];
// // calculate the polarization vectors for the gluons
// for(unsigned int iw=0;iw<3;++iw) {
// VectorWaveFunction gwave(momenta[iw+1],partons[iw+1],outgoing);
// for(unsigned int ix=0;ix<2;++ix) {
// //if(iw==2) gwave.reset(10);
// //else
// gwave.reset(2*ix);
// g[iw][ix] = gwave.wave();
// }
// }
Energy2 mh2 = momenta[0].mass2();
Energy2 s = (momenta[1]+momenta[2]).m2();
Energy2 t = (momenta[1]+momenta[3]).m2();
Energy2 u = (momenta[2]+momenta[3]).m2();
// calculate the loop functions
Complex A4stu(0.),A2stu(0.),A2tsu(0.),A2ust(0.);
for(int ix=_minloop;ix<=_maxloop;++ix) {
// loop functions
if(_massopt==0) {
Energy2 mf2=sqr(getParticleData(ix)->mass());
A4stu+=A4(s,t,u,mf2);
A2stu+=A2(s,t,u,mf2);
A2tsu+=A2(t,s,u,mf2);
A2ust+=A2(u,s,t,mf2);
}
else {
A4stu=-1./3.;
A2stu=-sqr(s/mh2)/3.;
A2tsu=-sqr(t/mh2)/3.;
A2ust=-sqr(u/mh2)/3.;
}
}
// Complex A3stu=0.5*(A2stu+A2ust+A2tsu-A4stu);
// // compute the dot products for the matrix element
// // and polarization vector * momenta
// Energy2 pdot[3][3];
// complex<InvEnergy> eps[3][3][2];
// for(unsigned int ig=0;ig<3;++ig) {
// for(unsigned int ip=0;ip<3;++ip) {
// pdot[ig][ip]=momenta[ig+1]*momenta[ip+1];
// for(unsigned int ih=0;ih<2;++ih) {
// if(ig!=ip)
// eps[ig][ip][ih]=g[ig][ih].dot(momenta[ip+1])/pdot[ig][ip];
// else
// eps[ig][ip][ih]=ZERO;
// }
// }
// }
// prefactors
Energy mw(getParticleData(ParticleID::Wplus)->mass());
// Energy3 pre=sqr(mh2)/mw;
// // compute the matrix element
// double output(0.);
// complex<InvEnergy2> wdot[3][3];
// for(unsigned int ghel1=0;ghel1<2;++ghel1) {
// for(unsigned int ghel2=0;ghel2<2;++ghel2) {
// for(unsigned int ghel3=0;ghel3<2;++ghel3) {
// wdot[0][1]=g[0][ghel1].dot(g[1][ghel2])/pdot[0][1];
// wdot[0][2]=g[0][ghel1].dot(g[2][ghel3])/pdot[0][2];
// wdot[1][0]=wdot[0][1];
// wdot[1][2]=g[1][ghel2].dot(g[2][ghel3])/pdot[1][2];
// wdot[2][0]=wdot[0][2];
// wdot[2][1]=wdot[1][2];
// // last piece
// Complex diag=pre*A3stu*(eps[0][2][ghel1]*eps[1][0][ghel2]*eps[2][1][ghel3]-
// eps[0][1][ghel1]*eps[1][2][ghel2]*eps[2][0][ghel3]+
// (eps[2][0][ghel3]-eps[2][1][ghel3])*wdot[0][1]+
// (eps[1][2][ghel2]-eps[1][0][ghel2])*wdot[0][2]+
// (eps[0][1][ghel1]-eps[0][2][ghel1])*wdot[1][2]);
// // first piece
// diag+=pre*(+A2stu*(eps[0][1][ghel1]*eps[1][0][ghel2]-wdot[0][1])*
// (eps[2][0][ghel3]-eps[2][1][ghel3])
// +A2tsu*(eps[0][2][ghel1]*eps[2][0][ghel3]-wdot[0][2])*
// (eps[1][2][ghel2]-eps[1][0][ghel2])
// +A2ust*(eps[1][2][ghel2]*eps[2][1][ghel3]-wdot[1][2])*
// (eps[0][1][ghel1]-eps[0][2][ghel1]));
// output+=norm(diag);
// }
// }
// }
// // colour factor and what's left of the prefactor
// output *= 6.;
double me=4.*24./s/t/u*pow<4,1>(mh2)/sqr(mw)*
(norm(A2stu)+norm(A2ust)+norm(A2tsu)+norm(A4stu));
return me;
}
Energy2 SMHiggsGGHiggsPPDecayer::loME(Energy mh) const {
Complex loop(0.);
Energy2 mh2(sqr(mh));
Energy mw(getParticleData(ParticleID::Wplus)->mass());
for(int ix=_minloop;ix<=_maxloop;++ix) {
// loop functions
if(_massopt==0) {
Energy2 mf2=sqr(getParticleData(ix)->mass());
loop += A1(mh2,mf2);
}
else {
loop += 2./3.;
}
}
return 1./Constants::pi*sqr(mh2)/sqr(mw)*norm(loop);
}
diff --git a/Decay/PerturbativeDecayer.cc b/Decay/PerturbativeDecayer.cc
--- a/Decay/PerturbativeDecayer.cc
+++ b/Decay/PerturbativeDecayer.cc
@@ -1,1068 +1,1089 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PerturbativeDecayer class.
//
#include "PerturbativeDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/EnumIO.h"
using namespace Herwig;
void PerturbativeDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(pTmin_,GeV) << oenum(inter_) << alphaS_ << alphaEM_
<< useMEforT2_ << C_ << ymax_ << phaseOpt_;
}
void PerturbativeDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(pTmin_,GeV) >> ienum(inter_) >> alphaS_ >> alphaEM_
>> useMEforT2_ >> C_ >> ymax_ >> phaseOpt_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeAbstractClass<PerturbativeDecayer,DecayIntegrator>
describeHerwigPerturbativeDecayer("Herwig::PerturbativeDecayer",
"Herwig.so HwPerturbativeDecay.so");
void PerturbativeDecayer::Init() {
static ClassDocumentation<PerturbativeDecayer> documentation
("The PerturbativeDecayer class is the mase class for "
"perturbative decays in Herwig");
static Parameter<PerturbativeDecayer,Energy> interfacepTmin
("pTmin",
"Minimum transverse momentum from gluon radiation",
&PerturbativeDecayer::pTmin_, GeV, 1.0*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Switch<PerturbativeDecayer,ShowerInteraction> interfaceInteractions
("Interactions",
"which interactions to include for the hard corrections",
&PerturbativeDecayer::inter_, ShowerInteraction::QCD, false, false);
static SwitchOption interfaceInteractionsQCD
(interfaceInteractions,
"QCD",
"QCD Only",
ShowerInteraction::QCD);
static SwitchOption interfaceInteractionsQED
(interfaceInteractions,
"QED",
"QED only",
ShowerInteraction::QED);
static SwitchOption interfaceInteractionsQCDandQED
(interfaceInteractions,
"QCDandQED",
"Both QCD and QED",
ShowerInteraction::Both);
static Reference<PerturbativeDecayer,ShowerAlpha> interfaceAlphaS
("AlphaS",
"Object for the coupling in the generation of hard QCD radiation",
&PerturbativeDecayer::alphaS_, false, false, true, true, false);
static Reference<PerturbativeDecayer,ShowerAlpha> interfaceAlphaEM
("AlphaEM",
"Object for the coupling in the generation of hard QED radiation",
&PerturbativeDecayer::alphaEM_, false, false, true, true, false);
static Switch<PerturbativeDecayer,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",
&PerturbativeDecayer::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);
static Parameter<PerturbativeDecayer,double> interfacePrefactor
("Prefactor",
"The prefactor for the sampling of the powheg Sudakov",
&PerturbativeDecayer::C_, 6.3, 0.0, 1e10,
false, false, Interface::limited);
static Parameter<PerturbativeDecayer,double> interfaceYMax
("YMax",
"The maximum value for the rapidity",
&PerturbativeDecayer::ymax_, 10., 0.0, 100.,
false, false, Interface::limited);
static Switch<PerturbativeDecayer,unsigned int> interfacePhaseSpaceOption
("PhaseSpaceOption",
"Option for the phase-space sampling",
&PerturbativeDecayer::phaseOpt_, 0, false, false);
static SwitchOption interfacePhaseSpaceOptionFixedYLimits
(interfacePhaseSpaceOption,
"FixedYLimits",
"Use a fixed limit for the rapidity",
0);
static SwitchOption interfacePhaseSpaceOptionVariableYLimits
(interfacePhaseSpaceOption,
"VariableYLimits",
"Change limit for the rapidity with pT",
1);
}
double PerturbativeDecayer::matrixElementRatio(const Particle & ,
const ParticleVector & ,
const ParticleVector & ,
MEOption ,
ShowerInteraction ) {
throw Exception() << "Base class PerturbativeDecayer::matrixElementRatio() "
<< "called, should have an implementation in the inheriting class"
<< Exception::runerror;
return 0.;
}
RealEmissionProcessPtr PerturbativeDecayer::generateHardest(RealEmissionProcessPtr born) {
return getHardEvent(born,false,inter_);
}
RealEmissionProcessPtr PerturbativeDecayer::applyHardMatrixElementCorrection(RealEmissionProcessPtr born) {
return getHardEvent(born,true,ShowerInteraction::QCD);
}
RealEmissionProcessPtr PerturbativeDecayer::getHardEvent(RealEmissionProcessPtr born,
bool inDeadZone,
ShowerInteraction inter) {
// check one incoming
assert(born->bornIncoming().size()==1);
// check exactly two outgoing particles
assert(born->bornOutgoing().size()==2); // search for coloured particles
bool colouredParticles=born->bornIncoming()[0]->dataPtr()->coloured();
bool chargedParticles=born->bornIncoming()[0]->dataPtr()->charged();
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
if(born->bornOutgoing()[ix]->dataPtr()->coloured())
colouredParticles=true;
if(born->bornOutgoing()[ix]->dataPtr()->charged())
chargedParticles=true;
}
// if no coloured/charged particles return
if ( !colouredParticles && !chargedParticles ) return RealEmissionProcessPtr();
if ( !colouredParticles && inter==ShowerInteraction::QCD ) return RealEmissionProcessPtr();
if ( ! chargedParticles && inter==ShowerInteraction::QED ) return RealEmissionProcessPtr();
// for decay b -> a c
// set progenitors
PPtr cProgenitor = born->bornOutgoing()[0];
PPtr aProgenitor = born->bornOutgoing()[1];
// get the decaying particle
PPtr bProgenitor = born->bornIncoming()[0];
// identify which dipoles are required
vector<DipoleType> dipoles;
if(!identifyDipoles(dipoles,aProgenitor,bProgenitor,cProgenitor,inter)) {
return RealEmissionProcessPtr();
}
Energy trialpT = pTmin_;
LorentzRotation eventFrame;
vector<Lorentz5Momentum> momenta;
vector<Lorentz5Momentum> trialMomenta(4);
PPtr finalEmitter, finalSpectator;
PPtr trialEmitter, trialSpectator;
DipoleType finalType(FFa,ShowerInteraction::QCD);
for (int i=0; i<int(dipoles.size()); ++i) {
+ if(dipoles[i].type==FFg) continue;
// assign emitter and spectator based on current dipole
if (dipoles[i].type==FFc || dipoles[i].type==IFc || dipoles[i].type==IFbc) {
trialEmitter = cProgenitor;
trialSpectator = aProgenitor;
}
else if (dipoles[i].type==FFa || dipoles[i].type==IFa || dipoles[i].type==IFba) {
trialEmitter = aProgenitor;
trialSpectator = cProgenitor;
}
// find rotation from lab to frame with the spectator along -z
LorentzRotation trialEventFrame(bProgenitor->momentum().findBoostToCM());
Lorentz5Momentum pspectator = (trialEventFrame*trialSpectator->momentum());
trialEventFrame.rotateZ( -pspectator.phi() );
trialEventFrame.rotateY( -pspectator.theta() - Constants::pi );
// invert it
trialEventFrame.invert();
// try to generate an emission
pT_ = pTmin_;
vector<Lorentz5Momentum> trialMomenta
= hardMomenta(bProgenitor, trialEmitter, trialSpectator,
dipoles, i, inDeadZone);
// select dipole which gives highest pT emission
if(pT_>trialpT) {
trialpT = pT_;
momenta = trialMomenta;
eventFrame = trialEventFrame;
finalEmitter = trialEmitter;
finalSpectator = trialSpectator;
finalType = dipoles[i];
if (dipoles[i].type==FFc || dipoles[i].type==FFa ) {
if((momenta[3]+momenta[1]).m2()-momenta[1].m2()>
(momenta[3]+momenta[2]).m2()-momenta[2].m2()) {
swap(finalEmitter,finalSpectator);
swap(momenta[1],momenta[2]);
}
}
}
}
pT_ = trialpT;
// if no emission return
if(momenta.empty()) {
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QCD)
born->pT()[ShowerInteraction::QCD] = pTmin_;
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QED)
born->pT()[ShowerInteraction::QED] = pTmin_;
return born;
}
// rotate momenta back to the lab
for(unsigned int ix=0;ix<momenta.size();++ix) {
momenta[ix] *= eventFrame;
}
// set maximum pT for subsequent branchings
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QCD)
born->pT()[ShowerInteraction::QCD] = pT_;
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QED)
born->pT()[ShowerInteraction::QED] = pT_;
// get ParticleData objects
tcPDPtr b = bProgenitor ->dataPtr();
tcPDPtr e = finalEmitter ->dataPtr();
tcPDPtr s = finalSpectator->dataPtr();
tcPDPtr boson = getParticleData(finalType.interaction==ShowerInteraction::QCD ?
ParticleID::g : ParticleID::gamma);
// create new ShowerParticles
PPtr emitter = e ->produceParticle(momenta[1]);
PPtr spectator = s ->produceParticle(momenta[2]);
PPtr gauge = boson->produceParticle(momenta[3]);
PPtr incoming = b ->produceParticle(bProgenitor->momentum());
// insert the particles
born->incoming().push_back(incoming);
unsigned int iemit(0),ispect(0);
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
if(born->bornOutgoing()[ix]==finalEmitter) {
born->outgoing().push_back(emitter);
iemit = born->outgoing().size();
}
else if(born->bornOutgoing()[ix]==finalSpectator) {
born->outgoing().push_back(spectator);
ispect = born->outgoing().size();
}
}
born->outgoing().push_back(gauge);
if(!spectator->dataPtr()->coloured() ||
(finalType.type != FFa && finalType.type!=FFc) ) ispect = 0;
born->emitter(iemit);
born->spectator(ispect);
born->emitted(3);
// boost if being use as ME correction
if(inDeadZone) {
if(finalType.type==IFa || finalType.type==IFba) {
LorentzRotation trans(cProgenitor->momentum().findBoostToCM());
trans.boost(spectator->momentum().boostVector());
born->transformation(trans);
}
else if(finalType.type==IFc || finalType.type==IFbc) {
LorentzRotation trans(bProgenitor->momentum().findBoostToCM());
trans.boost(spectator->momentum().boostVector());
born->transformation(trans);
}
}
// set the interaction
born->interaction(finalType.interaction);
// set up colour lines
getColourLines(born);
// return the tree
return born;
}
bool PerturbativeDecayer::identifyDipoles(vector<DipoleType> & dipoles,
PPtr & aProgenitor,
PPtr & bProgenitor,
PPtr & cProgenitor,
ShowerInteraction inter) const {
+ enhance_ = 1.;
// identify any QCD dipoles
if(inter==ShowerInteraction::QCD ||
inter==ShowerInteraction::Both) {
PDT::Colour bColour = bProgenitor->dataPtr()->iColour();
PDT::Colour cColour = cProgenitor->dataPtr()->iColour();
PDT::Colour aColour = aProgenitor->dataPtr()->iColour();
// decaying colour singlet
if (bColour==PDT::Colour0 ) {
if ((cColour==PDT::Colour3 && aColour==PDT::Colour3bar) ||
(cColour==PDT::Colour3bar && aColour==PDT::Colour3) ||
(cColour==PDT::Colour8 && aColour==PDT::Colour8)){
dipoles.push_back(DipoleType(FFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc,ShowerInteraction::QCD));
+ if(aProgenitor->id()==ParticleID::g &&
+ cProgenitor->id()==ParticleID::g ) {
+ enhance_ = 1.5;
+ dipoles.push_back(DipoleType(FFg,ShowerInteraction::QCD));
+ }
}
}
// decaying colour triplet
else if (bColour==PDT::Colour3 ) {
if (cColour==PDT::Colour3 && aColour==PDT::Colour0){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour0 && aColour==PDT::Colour3){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour3){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour3 && aColour==PDT::Colour8){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
}
// decaying colour anti-triplet
else if (bColour==PDT::Colour3bar) {
if ((cColour==PDT::Colour3bar && aColour==PDT::Colour0)){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
}
else if ((cColour==PDT::Colour0 && aColour==PDT::Colour3bar)){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour3bar){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour3bar && aColour==PDT::Colour8){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
}
// decaying colour octet
else if (bColour==PDT::Colour8){
if ((cColour==PDT::Colour3 && aColour==PDT::Colour3bar) ||
(cColour==PDT::Colour3bar && aColour==PDT::Colour3)){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour0){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour0 && aColour==PDT::Colour8){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
}
}
}
// QED dipoles
if(inter==ShowerInteraction::Both ||
inter==ShowerInteraction::QED) {
const bool & bCharged = bProgenitor->dataPtr()->charged();
const bool & cCharged = cProgenitor->dataPtr()->charged();
const bool & aCharged = aProgenitor->dataPtr()->charged();
// initial-final
if(bCharged && aCharged) {
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QED));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QED));
}
if(bCharged && cCharged) {
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QED));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QED));
}
// final-state
if(aCharged && cCharged) {
dipoles.push_back(DipoleType(FFa,ShowerInteraction::QED));
dipoles.push_back(DipoleType(FFc,ShowerInteraction::QED));
}
}
// check colour structure is allowed
return !dipoles.empty();
}
vector<Lorentz5Momentum> PerturbativeDecayer::hardMomenta(PPtr in, PPtr emitter,
PPtr spectator,
const vector<DipoleType> &dipoles,
int i, bool inDeadZone) {
// get masses of the particles
mb_ = in ->momentum().mass();
e_ = emitter ->momentum().mass()/mb_;
s_ = spectator->momentum().mass()/mb_;
e2_ = sqr(e_);
s2_ = sqr(s_);
vector<Lorentz5Momentum> particleMomenta (4);
Energy2 lambda = sqr(mb_)*sqrt(1.+sqr(s2_)+sqr(e2_)-2.*s2_-2.*e2_-2.*s2_*e2_);
// calculate A
double A = 2.*C_/Constants::twopi;
// factor due sampling choice
if(phaseOpt_==0) A *=ymax_;
// coupling factor
if(dipoles[i].interaction==ShowerInteraction::QCD)
A *= alphaS() ->overestimateValue();
else
A *= alphaEM()->overestimateValue();
Energy pTmax = 0.5*mb_*(1.-sqr(s_+e_));
// if no possible branching return
if ( pTmax < pTmin_ ) {
particleMomenta.clear();
return particleMomenta;
}
Energy pT=pTmax;
while (pT >= pTmin_) {
double ymax;
// overestimate with flat y limit
if(phaseOpt_==0) {
pT *= pow(UseRandom::rnd(),(1./A));
ymax=ymax_;
}
// pT sampling including tighter pT dependent y limit
else {
pT = 2.*pTmax*exp(-sqrt(-2.*log(UseRandom::rnd())/A+sqr(log(2.*pTmax/pT))));
// choice of limit overestimate ln(2*pTmax/pT) (true limit acosh(pTmax/pT))
ymax = log(2.*pTmax/pT);
}
if (pT < pTmin_) {
particleMomenta.clear();
break;
}
double phi = UseRandom::rnd()*Constants::twopi;
double y = ymax*(2.*UseRandom::rnd()-1.);
double weight[2] = {0.,0.};
double xs[2], xe[2], xe_z[2], xg;
for (unsigned int j=0; j<2; j++) {
// check if the momenta are physical
if (!calcMomenta(j, pT, y, phi, xg, xs[j], xe[j], xe_z[j], particleMomenta))
continue;
// check if point lies within phase space
if (!psCheck(xg, xs[j]))
continue;
// check if point lies within the dead-zone (if required)
if(inDeadZone) {
if(!inTotalDeadZone(xg,xs[j],dipoles,i)) continue;
}
// decay products for 3 body decay
PPtr inpart = in ->dataPtr()->produceParticle(particleMomenta[0]);
ParticleVector decay3;
decay3.push_back(emitter ->dataPtr()->produceParticle(particleMomenta[1]));
decay3.push_back(spectator->dataPtr()->produceParticle(particleMomenta[2]));
if(dipoles[i].interaction==ShowerInteraction::QCD)
decay3.push_back(getParticleData(ParticleID::g )->produceParticle(particleMomenta[3]));
else
decay3.push_back(getParticleData(ParticleID::gamma)->produceParticle(particleMomenta[3]));
// decay products for 2 body decay
Lorentz5Momentum p1(ZERO,ZERO, lambda/2./mb_,(mb_/2.)*(1.+e2_-s2_),mb_*e_);
Lorentz5Momentum p2(ZERO,ZERO,-lambda/2./mb_,(mb_/2.)*(1.+s2_-e2_),mb_*s_);
ParticleVector decay2;
decay2.push_back(emitter ->dataPtr()->produceParticle(p1));
decay2.push_back(spectator->dataPtr()->produceParticle(p2));
if (dipoles[i].type==FFa || dipoles[i].type==IFa || dipoles[i].type==IFba) {
swap(decay2[0],decay2[1]);
swap(decay3[0],decay3[1]);
}
// calculate matrix element ratio R/B
double meRatio = matrixElementRatio(*inpart,decay2,decay3,Initialize,dipoles[i].interaction);
// calculate dipole factor
double dipoleSum(0.),numerator(0.);
for (int k=0; k<int(dipoles.size()); ++k) {
// skip dipoles which are not of the interaction being considered
if(dipoles[k].interaction!=dipoles[i].interaction) continue;
- double dipole = abs(calculateDipole(dipoles[k],*inpart,decay3));
- dipoleSum += dipole;
- if (k==i) numerator = dipole;
+ pair<double,double> dipole = calculateDipole(dipoles[k],*inpart,decay3);
+ dipoleSum += abs(dipole.first);
+ if (k==i) numerator = abs(dipole.second);
}
meRatio *= numerator/dipoleSum;
// calculate jacobian
Energy2 denom = (mb_-particleMomenta[3].e())*particleMomenta[2].vect().mag() -
particleMomenta[2].e()*particleMomenta[3].z();
InvEnergy2 J = (particleMomenta[2].vect().mag2())/(lambda*denom);
// calculate weight
- weight[j] = meRatio*fabs(sqr(pT)*J)/C_/Constants::twopi;
+ weight[j] = enhance_*meRatio*fabs(sqr(pT)*J)/C_/Constants::twopi;
if(dipoles[i].interaction==ShowerInteraction::QCD)
weight[j] *= alphaS() ->ratio(pT*pT);
else
weight[j] *= alphaEM()->ratio(pT*pT);
}
// accept point if weight > R
if (weight[0] + weight[1] > UseRandom::rnd()) {
if (weight[0] > (weight[0] + weight[1])*UseRandom::rnd()) {
particleMomenta[1].setE( (mb_/2.)*xe [0]);
particleMomenta[1].setZ( (mb_/2.)*xe_z[0]);
particleMomenta[2].setE( (mb_/2.)*xs [0]);
particleMomenta[2].setZ(-(mb_/2.)*sqrt(sqr(xs[0])-4.*s2_));
}
else {
particleMomenta[1].setE( (mb_/2.)*xe [1]);
particleMomenta[1].setZ( (mb_/2.)*xe_z[1]);
particleMomenta[2].setE( (mb_/2.)*xs [1]);
particleMomenta[2].setZ(-(mb_/2.)*sqrt(sqr(xs[1])-4.*s2_));
}
pT_ = pT;
break;
}
}
return particleMomenta;
}
bool PerturbativeDecayer::calcMomenta(int j, Energy pT, double y, double phi,
double& xg, double& xs, double& xe, double& xe_z,
vector<Lorentz5Momentum>& particleMomenta){
// calculate xg
xg = 2.*pT*cosh(y) / mb_;
if (xg>(1. - sqr(e_ + s_)) || xg<0.) return false;
// calculate the two values of xs
double xT = 2.*pT / mb_;
double A = 4.-4.*xg+sqr(xT);
double B = 4.*(3.*xg-2.+2.*e2_-2.*s2_-sqr(xg)-xg*e2_+xg*s2_);
double L = 1.+sqr(s2_)+sqr(e2_)-2.*s2_-2.*e2_-2.*s2_*e2_;
double det = 16.*( -L*sqr(xT)+pow(xT,4)*s2_+2.*xg*sqr(xT)*(1.-s2_-e2_)+
L*sqr(xg)-sqr(xg*xT)*(1.+s2_)+pow(xg,4)+
2.*pow(xg,3)*(-1.+s2_+e2_) );
if (det<0.) return false;
if (j==0) xs = (-B+sqrt(det))/(2.*A);
if (j==1) xs = (-B-sqrt(det))/(2.*A);
// check value of xs is physical
if (xs>(1.+s2_-e2_) || xs<2.*s_) return false;
// calculate xe
xe = 2.-xs-xg;
// check value of xe is physical
if (xe>(1.+e2_-s2_) || xe<2.*e_) return false;
// calculate xe_z
double root1 = sqrt(max(0.,sqr(xs)-4.*s2_)), root2 = sqrt(max(0.,sqr(xe)-4.*e2_-sqr(xT)));
double epsilon_p = -root1+xT*sinh(y)+root2;
double epsilon_m = -root1+xT*sinh(y)-root2;
// find direction of emitter
if (fabs(epsilon_p) < 1.e-10) xe_z = sqrt(sqr(xe)-4.*e2_-sqr(xT));
else if (fabs(epsilon_m) < 1.e-10) xe_z = -sqrt(sqr(xe)-4.*e2_-sqr(xT));
else return false;
// check the emitter is on shell
if (fabs((sqr(xe)-sqr(xT)-sqr(xe_z)-4.*e2_))>1.e-10) return false;
// calculate 4 momenta
particleMomenta[0].setE ( mb_);
particleMomenta[0].setX ( ZERO);
particleMomenta[0].setY ( ZERO);
particleMomenta[0].setZ ( ZERO);
particleMomenta[0].setMass( mb_);
particleMomenta[1].setE ( mb_*xe/2.);
particleMomenta[1].setX (-pT*cos(phi));
particleMomenta[1].setY (-pT*sin(phi));
particleMomenta[1].setZ ( mb_*xe_z/2.);
particleMomenta[1].setMass( mb_*e_);
particleMomenta[2].setE ( mb_*xs/2.);
particleMomenta[2].setX ( ZERO);
particleMomenta[2].setY ( ZERO);
particleMomenta[2].setZ (-mb_*sqrt(sqr(xs)-4.*s2_)/2.);
particleMomenta[2].setMass( mb_*s_);
particleMomenta[3].setE ( pT*cosh(y));
particleMomenta[3].setX ( pT*cos(phi));
particleMomenta[3].setY ( pT*sin(phi));
particleMomenta[3].setZ ( pT*sinh(y));
particleMomenta[3].setMass( ZERO);
return true;
}
bool PerturbativeDecayer::psCheck(const double xg, const double xs) {
// check is point is in allowed region of phase space
double xe_star = (1.-s2_+e2_-xg)/sqrt(1.-xg);
double xg_star = xg/sqrt(1.-xg);
if ((sqr(xe_star)-4.*e2_) < 1e-10) return false;
double xs_max = (4.+4.*s2_-sqr(xe_star+xg_star)+
sqr(sqrt(sqr(xe_star)-4.*e2_)+xg_star))/ 4.;
double xs_min = (4.+4.*s2_-sqr(xe_star+xg_star)+
sqr(sqrt(sqr(xe_star)-4.*e2_)-xg_star))/ 4.;
if (xs < xs_min || xs > xs_max) return false;
return true;
}
-double PerturbativeDecayer::calculateDipole(const DipoleType & dipoleId,
- const Particle & inpart,
- const ParticleVector & decay3) {
+pair<double,double> PerturbativeDecayer::calculateDipole(const DipoleType & dipoleId,
+ const Particle & inpart,
+ const ParticleVector & decay3) {
// calculate dipole for decay b->ac
- double dipole(0.);
+ pair<double,double> dipole = make_pair(0.,0.);
double x1 = 2.*decay3[0]->momentum().e()/mb_;
double x2 = 2.*decay3[1]->momentum().e()/mb_;
double xg = 2.*decay3[2]->momentum().e()/mb_;
double mu12 = sqr(decay3[0]->mass()/mb_);
double mu22 = sqr(decay3[1]->mass()/mb_);
tcPDPtr part[3] = {inpart.dataPtr(),decay3[0]->dataPtr(),decay3[1]->dataPtr()};
if(dipoleId.type==FFa || dipoleId.type == IFa || dipoleId.type == IFba) {
swap(part[1],part[2]);
swap(x1,x2);
swap(mu12,mu22);
}
// radiation from b with initial-final connection
if (dipoleId.type==IFba || dipoleId.type==IFbc) {
- dipole = -2./sqr(xg);
- dipole *= colourCoeff(part[0],part[1],part[2],dipoleId);
+ dipole.first = -2./sqr(xg);
+ dipole.first *= colourCoeff(part[0],part[1],part[2],dipoleId);
}
// radiation from a/c with initial-final connection
else if (dipoleId.type==IFa || dipoleId.type==IFc) {
double z = 1. - xg/(1.-mu22+mu12);
- dipole = (-2.*mu12/sqr(1.-x2+mu22-mu12) + (1./(1.-x2+mu22-mu12))*
+ dipole.first = (-2.*mu12/sqr(1.-x2+mu22-mu12) + (1./(1.-x2+mu22-mu12))*
(2./(1.-z)-dipoleSpinFactor(part[1],z)));
- dipole *= colourCoeff(part[1],part[0],part[2],dipoleId);
+ dipole.first *= colourCoeff(part[1],part[0],part[2],dipoleId);
}
// radiation from a/c with final-final connection
else if (dipoleId.type==FFa || dipoleId.type==FFc) {
double z = 1. + ((x1-1.+mu22-mu12)/(x2-2.*mu22));
double y = (1.-x2-mu12+mu22)/(1.-mu12-mu22);
double vt = sqrt((1.-sqr(e_+s_))*(1.-sqr(e_-s_)))/(1.-mu12-mu22);
double v = sqrt(sqr(2.*mu22+(1.-mu12-mu22)*(1.-y))-4.*mu22)
/(1.-y)/(1.-mu12-mu22);
if(part[1]->iSpin()!=PDT::Spin1) {
- dipole = (1./(1.-x2+mu22-mu12))*
+ dipole.first = (1./(1.-x2+mu22-mu12))*
((2./(1.-z*(1.-y)))-vt/v*(dipoleSpinFactor(part[1],z)+(2.*mu12/(1.+mu22-mu12-x2))));
}
else {
- dipole = (1./(1.-x2+mu22-mu12))*
+ dipole.first = (1./(1.-x2+mu22-mu12))*
(1./(1.-z*(1.-y))+1./(1.-(1.-z)*(1.-y))+(z*(1.-z)-2.)/v-vt/v*(2.*mu12/(1.+mu22-mu12-x2)));
+ dipole.second = (1./(1.-x2+mu22-mu12))*
+ (2./(1.-z*(1.-y))+(z*(1.-z)-2.)/v-vt/v*(2.*mu12/(1.+mu22-mu12-x2)));
+ dipole.second *= colourCoeff(part[1],part[2],part[0],dipoleId);
}
- dipole *= colourCoeff(part[1],part[2],part[0],dipoleId);
+ dipole.first *= colourCoeff(part[1],part[2],part[0],dipoleId);
+ }
+ // special for the case that all particles are gluons
+ else if(dipoleId.type==FFg) {
+ double z = (1.-x2)/xg;
+ double y = 1.-xg;
+ dipole.first = 1./(1.-xg)*(1./(1.-z*(1.-y))+1./(1.-(1.-z)*(1.-y))+(z*(1.-z)-2.));
+ dipole.first *= colourCoeff(part[1],part[2],part[0],dipoleId);
}
else
assert(false);
// coupling prefactors
- dipole *= 8.*Constants::pi;
+ if(dipole.second==0.) dipole.second=dipole.first;
+ dipole.first *= 8.*Constants::pi;
+ dipole.second *= 8.*Constants::pi;
// return the answer
return dipole;
}
double PerturbativeDecayer::dipoleSpinFactor(tcPDPtr part, double z){
// calculate the spin dependent component of the dipole
if (part->iSpin()==PDT::Spin0)
return 2.;
else if (part->iSpin()==PDT::Spin1Half)
return (1. + z);
else if (part->iSpin()==PDT::Spin1)
return -(z*(1.-z) - 1./(1.-z) + 1./z -2.);
return 0.;
}
double PerturbativeDecayer::colourCoeff(tcPDPtr emitter,
tcPDPtr spectator,
tcPDPtr other,
DipoleType dipole) {
if(dipole.interaction==ShowerInteraction::QCD) {
// calculate the colour factor of the dipole
double numerator=1.;
double denominator=1.;
if (emitter->iColour()!=PDT::Colour0 &&
spectator->iColour()!=PDT::Colour0 &&
other->iColour()!=PDT::Colour0) {
if (emitter->iColour() ==PDT::Colour3 ||
emitter->iColour() ==PDT::Colour3bar) numerator=-4./3;
else if (emitter->iColour() ==PDT::Colour8) numerator=-3. ;
denominator=-1.*numerator;
if (spectator->iColour()==PDT::Colour3 ||
spectator->iColour()==PDT::Colour3bar) numerator-=4./3;
else if (spectator->iColour()==PDT::Colour8) numerator-=3. ;
if (other->iColour() ==PDT::Colour3 ||
other->iColour() ==PDT::Colour3bar) numerator+=4./3;
else if (other->iColour() ==PDT::Colour8) numerator+=3. ;
numerator*=(-1./2.);
}
if (emitter->iColour()==PDT::Colour3 ||
emitter->iColour()== PDT::Colour3bar) numerator*=4./3.;
else if (emitter->iColour()==PDT::Colour8 &&
spectator->iColour()!=PDT::Colour8) numerator*=3.;
else if (emitter->iColour()==PDT::Colour8 &&
spectator->iColour()==PDT::Colour8) numerator*=6.;
return (numerator/denominator);
}
else {
double val = double(emitter->iCharge()*spectator->iCharge())/9.;
// FF dipoles
if(dipole.type==FFa || dipole.type == FFc) {
return val;
}
else {
return -val;
}
}
}
void PerturbativeDecayer::getColourLines(RealEmissionProcessPtr real) {
// extract the particles
vector<tPPtr> branchingPart;
branchingPart.push_back(real->incoming()[0]);
for(unsigned int ix=0;ix<real->outgoing().size();++ix) {
branchingPart.push_back(real->outgoing()[ix]);
}
vector<unsigned int> sing,trip,atrip,oct,sex,asex;
for (size_t ib=0;ib<branchingPart.size()-1;++ib) {
if (branchingPart[ib]->dataPtr()->iColour()==PDT::Colour0 ) sing. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour3 ) trip. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour3bar) atrip.push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour8 ) oct. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour6 ) sex. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour6bar) asex. push_back(ib);
}
// decaying colour singlet
if (branchingPart[0]->dataPtr()->iColour()==PDT::Colour0) {
// 0 -> 3 3bar
if (trip.size()==1 && atrip.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[atrip[0]]->colourConnect(branchingPart[ 3 ]);
branchingPart[ 3 ]->colourConnect(branchingPart[trip[0]]);
}
else {
branchingPart[atrip[0]]->colourConnect(branchingPart[trip[0]]);
}
}
// 0 -> 8 8
else if (oct.size()==2 ) {
if(real->interaction()==ShowerInteraction::QCD) {
bool col = UseRandom::rndbool();
branchingPart[oct[0]]->colourConnect(branchingPart[ 3 ],col);
branchingPart[ 3 ]->colourConnect(branchingPart[oct[1]],col);
branchingPart[oct[1]]->colourConnect(branchingPart[oct[0]],col);
}
else {
branchingPart[oct[0]]->colourConnect(branchingPart[oct[1]]);
branchingPart[oct[1]]->colourConnect(branchingPart[oct[0]]);
}
}
else
assert(real->interaction()==ShowerInteraction::QED);
}
// decaying colour triplet
else if (branchingPart[0]->dataPtr()->iColour()==PDT::Colour3 ){
// 3 -> 3 0
if (trip.size()==2 && sing.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[3]->incomingColour(branchingPart[trip[0]]);
branchingPart[3]-> colourConnect(branchingPart[trip[1]]);
}
else {
branchingPart[trip[1]]->incomingColour(branchingPart[trip[0]]);
}
}
// 3 -> 3 8
else if (trip.size()==2 && oct.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
// 8 emit incoming partner
if(real->emitter()==oct[0]&&real->spectator()==0) {
branchingPart[ 3 ]->incomingColour(branchingPart[trip[0]]);
branchingPart[ 3 ]-> colourConnect(branchingPart[oct[0] ]);
branchingPart[oct[0]]-> colourConnect(branchingPart[trip[1]]);
}
// 8 emit final spectator or vice veras
else {
branchingPart[oct[0]]->incomingColour(branchingPart[trip[0]]);
branchingPart[oct[0]]-> colourConnect(branchingPart[ 3 ]);
branchingPart[ 3 ]-> colourConnect(branchingPart[trip[1]]);
}
}
else {
branchingPart[oct[0]]->incomingColour(branchingPart[trip[0]]);
branchingPart[oct[0]]-> colourConnect(branchingPart[trip[1]]);
}
}
// 3 -> 3bar 3bar
else if(trip.size() ==1 && atrip.size()==2) {
if(real->interaction()==ShowerInteraction::QCD)
assert(false);
else {
tColinePtr col[3] = {ColourLine::create(branchingPart[ trip[0]],false),
ColourLine::create(branchingPart[atrip[0]],true ),
ColourLine::create(branchingPart[atrip[1]],true)};
col[0]->setSinkNeighbours(col[1],col[2]);
}
}
else
assert(false);
}
// decaying colour anti-triplet
else if (branchingPart[0]->dataPtr()->iColour()==PDT::Colour3bar) {
// 3bar -> 3bar 0
if (atrip.size()==2 && sing.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[3]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[3]-> colourConnect(branchingPart[atrip[1]],true);
}
else {
branchingPart[atrip[1]]->incomingColour(branchingPart[atrip[0]],true);
}
}
// 3 -> 3 8
else if (atrip.size()==2 && oct.size()==1){
if(real->interaction()==ShowerInteraction::QCD) {
// 8 emit incoming partner
if(real->emitter()==oct[0]&&real->spectator()==0) {
branchingPart[ 3 ]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[ 3 ]-> colourConnect(branchingPart[oct[0] ],true);
branchingPart[oct[0]]-> colourConnect(branchingPart[atrip[1]],true);
}
// 8 emit final spectator or vice veras
else {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[oct[0]]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[oct[0]]-> colourConnect(branchingPart[ 3 ],true);
branchingPart[3]-> colourConnect(branchingPart[atrip[1]] ,true);
}
}
}
else {
branchingPart[oct[0]]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[oct[0]]-> colourConnect(branchingPart[atrip[1]],true);
}
}
// 3bar -> bar bar
else if(atrip.size() ==1 && trip.size()==2) {
if(real->interaction()==ShowerInteraction::QCD)
assert(false);
else {
tColinePtr col[3] = {ColourLine::create(branchingPart[atrip[0]],true ),
ColourLine::create(branchingPart[ trip[0]],false),
ColourLine::create(branchingPart[ trip[1]],false)};
col[0]->setSourceNeighbours(col[1],col[2]);
}
}
else
assert(false);
}
// decaying colour octet
else if(branchingPart[0]->dataPtr()->iColour()==PDT::Colour8 ) {
// 8 -> 3 3bar
if (trip.size()==1 && atrip.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
// 3 emits
if(trip[0]==real->emitter()) {
branchingPart[3] ->incomingColour(branchingPart[oct[0]] );
branchingPart[3] -> colourConnect(branchingPart[trip[0]]);
branchingPart[atrip[0]]->incomingColour(branchingPart[oct[0]],true);
}
// 3bar emits
else {
branchingPart[3] ->incomingColour(branchingPart[oct[0]] ,true);
branchingPart[3] -> colourConnect(branchingPart[atrip[0]],true);
branchingPart[trip[0]]->incomingColour(branchingPart[oct[0]] );
}
}
else {
branchingPart[trip[0]]->incomingColour(branchingPart[oct[0]] );
branchingPart[atrip[0]]->incomingColour(branchingPart[oct[0]],true);
}
}
// 8 -> 8 0
else if (sing.size()==1 && oct.size()==2) {
if(real->interaction()==ShowerInteraction::QCD) {
bool col = UseRandom::rndbool();
branchingPart[ 3 ]->colourConnect (branchingPart[oct[1]], col);
branchingPart[ 3 ]->incomingColour(branchingPart[oct[0]], col);
branchingPart[oct[1]]->incomingColour(branchingPart[oct[0]],!col);
}
else {
branchingPart[oct[1]]->incomingColour(branchingPart[oct[0]]);
branchingPart[oct[1]]->incomingColour(branchingPart[oct[0]],true);
}
}
else
assert(false);
}
// sextet
else if(branchingPart[0]->dataPtr()->iColour() == PDT::Colour6) {
if(trip.size()==2) {
assert(real->interaction()==ShowerInteraction::QED);
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(branchingPart[0]->colourInfo());
for(unsigned int ix=0;ix<2;++ix) {
ColinePtr cline = new_ptr(ColourLine());
parentColour->colourLine(cline);
cline->addColoured(branchingPart[trip[ix]]);
}
}
else
assert(false);
}
// antisextet
else if(branchingPart[0]->dataPtr()->iColour() == PDT::Colour6bar) {
if(atrip.size()==2) {
assert(real->interaction()==ShowerInteraction::QED);
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(branchingPart[0]->colourInfo());
for(unsigned int ix=0;ix<2;++ix) {
ColinePtr cline = new_ptr(ColourLine());
parentColour->antiColourLine(cline);
cline->addColoured(branchingPart[atrip[ix]],true);
}
}
else
assert(false);
}
else
assert(false);
}
PerturbativeDecayer::phaseSpaceRegion
PerturbativeDecayer::inInitialFinalDeadZone(double xg, double xa,
double a, double c) const {
double lam = sqrt(1.+a*a+c*c-2.*a-2.*c-2.*a*c);
double kappab = 1.+0.5*(1.-a+c+lam);
double kappac = kappab-1.+c;
double kappa(0.);
// check whether or not in the region for emission from c
double r = 0.5;
if(c!=0.) r += 0.5*c/(1.+a-xa);
double pa = sqrt(sqr(xa)-4.*a);
double z = ((2.-xa)*(1.-r)+r*pa-xg)/pa;
if(z<1. && z>0.) {
kappa = (1.+a-c-xa)/(z*(1.-z));
if(kappa<kappac)
return emissionFromC;
}
// check in region for emission from b (T1)
double cq = sqr(1.+a-c)-4*a;
double bq = -2.*kappab*(1.-a-c);
double aq = sqr(kappab)-4.*a*(kappab-1);
double dis = sqr(bq)-4.*aq*cq;
z=1.-(-bq-sqrt(dis))/2./aq;
double w = 1.-(1.-z)*(kappab-1.);
double xgmax = (1.-z)*kappab;
// possibly in T1 region
if(xg<xgmax) {
z = 1.-xg/kappab;
kappa=kappab;
}
// possibly in T2 region
else {
aq = 4.*a;
bq = -4.*a*(2.-xg);
cq = sqr(1.+a-c-xg);
dis = sqr(bq)-4.*aq*cq;
z = (-bq-sqrt(dis))/2./aq;
kappa = xg/(1.-z);
}
// compute limit on xa
double u = 1.+a-c-(1.-z)*kappa;
w = 1.-(1.-z)*(kappa-1.);
double v = sqr(u)-4.*z*a*w;
if(v<0. && v>-1e-10) v= 0.;
v = sqrt(v);
if(xa<0.5*((u+v)/w+(u-v)/z)) {
if(xg<xgmax)
return emissionFromA1;
else if(useMEforT2_)
return deadZone;
else
return emissionFromA2;
}
else
return deadZone;
}
PerturbativeDecayer::phaseSpaceRegion
PerturbativeDecayer::inFinalFinalDeadZone(double xb, double xc,
double b, double c) const {
// basic kinematics
double lam = sqrt(1.+b*b+c*c-2.*b-2.*c-2.*b*c);
// check whether or not in the region for emission from b
double r = 0.5;
if(b!=0.) r+=0.5*b/(1.+c-xc);
double pc = sqrt(sqr(xc)-4.*c);
double z = -((2.-xc)*r-r*pc-xb)/pc;
if(z<1. and z>0.) {
if((1.-b+c-xc)/(z*(1.-z))<0.5*(1.+b-c+lam)) return emissionFromB;
}
// check whether or not in the region for emission from c
r = 0.5;
if(c!=0.) r+=0.5*c/(1.+b-xb);
double pb = sqrt(sqr(xb)-4.*b);
z = -((2.-xb)*r-r*pb-xc)/pb;
if(z<1. and z>0.) {
if((1.-c+b-xb)/(z*(1.-z))<0.5*(1.-b+c+lam)) return emissionFromC;
}
return deadZone;
}
bool PerturbativeDecayer::inTotalDeadZone(double xg, double xs,
const vector<DipoleType> & dipoles,
int i) {
double xb,xc,b,c;
if(dipoles[i].type==FFa || dipoles[i].type == IFa || dipoles[i].type == IFba) {
xc = xs;
xb = 2.-xg-xs;
b = e2_;
c = s2_;
}
else {
xb = xs;
xc = 2.-xg-xs;
b = s2_;
c = e2_;
}
for(unsigned int ix=0;ix<dipoles.size();++ix) {
if(dipoles[ix].interaction!=dipoles[i].interaction)
continue;
// should also remove negative QED dipoles but shouldn't be an issue unless we
// support QED ME corrections
switch (dipoles[ix].type) {
case FFa :
if(inFinalFinalDeadZone(xb,xc,b,c)!=deadZone) return false;
break;
case FFc :
if(inFinalFinalDeadZone(xc,xb,c,b)!=deadZone) return false;
break;
case IFa : case IFba:
if(inInitialFinalDeadZone(xg,xc,c,b)!=deadZone) return false;
break;
case IFc : case IFbc:
if(inInitialFinalDeadZone(xg,xb,b,c)!=deadZone) return false;
break;
+ case FFg:
+ break;
}
}
return true;
}
diff --git a/Decay/PerturbativeDecayer.h b/Decay/PerturbativeDecayer.h
--- a/Decay/PerturbativeDecayer.h
+++ b/Decay/PerturbativeDecayer.h
@@ -1,339 +1,345 @@
// -*- 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"
#include "Herwig/Shower/Core/ShowerInteraction.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};
+ enum dipoleType {FFa, FFc, IFa, IFc, IFba, IFbc, FFg};
/**
* Phase-space region for an emission (assumes \f$a\to b,c\f$
*/
enum phaseSpaceRegion {emissionFromB,emissionFromC,emissionFromA1,emissionFromA2,deadZone};
/**
* Type of dipole
*/
struct DipoleType {
DipoleType() {}
DipoleType(dipoleType a, ShowerInteraction b)
: type(a), interaction(b)
{}
dipoleType type;
ShowerInteraction interaction;
};
public:
/**
* The default constructor.
*/
PerturbativeDecayer() : inter_(ShowerInteraction::QCD),
pTmin_(GeV), useMEforT2_(true),
C_(6.3), ymax_(10.), phaseOpt_(0),
pT_(ZERO),mb_(ZERO), e_(0.),
- s_(0.), e2_(0.), s2_(0.)
+ s_(0.), e2_(0.), s2_(0.), enhance_(1.)
{}
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {return No;}
/**
* Member to generate the hardest emission in the POWHEG scheme
*/
virtual RealEmissionProcessPtr generateHardest(RealEmissionProcessPtr);
/**
* Apply the hard matrix element correction to a given hard process or decay
*/
virtual RealEmissionProcessPtr applyHardMatrixElementCorrection(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:
/**
* 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,
ShowerInteraction inter);
/**
* Work out the type of process
*/
bool identifyDipoles(vector<DipoleType> & dipoles,
PPtr & aProgenitor,
PPtr & bProgenitor,
PPtr & cProgenitor,
ShowerInteraction inter) const;
/**
* Coupling for the generation of hard QCD radiation
*/
ShowerAlphaPtr alphaS() {return alphaS_;}
/**
* Coupling for the generation of hard QED radiation
*/
ShowerAlphaPtr alphaEM() {return alphaEM_;}
/**
* Return the momenta including the hard emission
*/
vector<Lorentz5Momentum> hardMomenta(PPtr in, PPtr emitter,
PPtr spectator,
const vector<DipoleType> & dipoles,
int i, bool inDeadZone);
/**
* 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
*/
- double calculateDipole(const DipoleType & dipoleId, const Particle & inpart,
- const ParticleVector & decay3);
+ pair<double,double> calculateDipole(const DipoleType & dipoleId,
+ const Particle & inpart,
+ const ParticleVector & decay3);
/**
* Return contribution to dipole that depends on the spin of the emitter
*/
double dipoleSpinFactor(tcPDPtr emitter, double z);
/**
* Return the colour coefficient of the dipole
*/
double colourCoeff(tcPDPtr emitter, tcPDPtr spectator,
tcPDPtr other, DipoleType dipole);
/**
* Set up the colour lines
*/
void getColourLines(RealEmissionProcessPtr real);
/**
* Generate a hard emission
*/
RealEmissionProcessPtr getHardEvent(RealEmissionProcessPtr born,
bool inDeadZone,
ShowerInteraction inter);
/**
* Is the \f$x_g,x_s\f$ point in the dead-zone for all the dipoles
*/
bool inTotalDeadZone(double xg, double xs,
const vector<DipoleType> & dipoles,
int i);
/**
* Is the \f$x_g,x_a\f$ point in the dead-zone for an initial-final colour connection
*/
phaseSpaceRegion inInitialFinalDeadZone(double xg, double xa, double a, double c) const;
/**
* Is the \f$x_b,x_c\f$ point in the dead-zone for a final-final colour connection
*/
phaseSpaceRegion inFinalFinalDeadZone(double xb, double xc, double b, double c) const;
/**
* For me corrections use the shower or me for the T2 region
*/
bool useMEforT2() const {return useMEforT2_;}
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
*/
//@{
/**
* Which types of radiation to generate
*/
ShowerInteraction inter_;
/**
* Coupling for the generation of hard QCD radiation
*/
ShowerAlphaPtr alphaS_;
/**
* Coupling for the generation of hard QED radiation
*/
ShowerAlphaPtr alphaEM_;
/**
* Minimum \f$p_T\f$
*/
Energy pTmin_;
/**
* This flag determines whether the T2 region in the decay shower
* (JHEP12(2003)_045) is populated by the ME correction (true) or
* the shower from the decaying particle.
*/
bool useMEforT2_;
/**
* Prefactor for the sampling
*/
double C_;
/**
* Maximum value for y
*/
double ymax_;
/**
* Option for phase-space sampling
*/
unsigned int phaseOpt_;
//@}
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_;
+
+ /**
+ * Enhancement prefactor for special cases
+ */
+ mutable double enhance_;
//@}
};
}
#endif /* Herwig_PerturbativeDecayer_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:00 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805476
Default Alt Text
(64 KB)

Event Timeline