Page MenuHomeHEPForge
Files F7879704

No OneTemporary
Actions

View Options

diff --git a/Decay/Perturbative/SMTopDecayer.cc b/Decay/Perturbative/SMTopDecayer.cc
--- a/Decay/Perturbative/SMTopDecayer.cc
+++ b/Decay/Perturbative/SMTopDecayer.cc
@@ -1,1193 +1,1193 @@
// -*- C++ -*-
//
// SMTopDecayer.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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/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"
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 << _alpha
<< _initialenhance << _finalenhance << _xg_sampling << _useMEforT2;
}
void SMTopDecayer::persistentInput(PersistentIStream & is, int) {
is >> _wvertex >> _wquarkwgt >> _wleptonwgt >> _wplus >> _alpha
>> _initialenhance >> _finalenhance >> _xg_sampling >> _useMEforT2;
}
ClassDescription<SMTopDecayer> SMTopDecayer::initSMTopDecayer;
// Definition of the static class description member.
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);
static Reference<SMTopDecayer,ShowerAlpha> interfaceCoupling
("Coupling",
"Pointer to the object to calculate the coupling for the correction",
&SMTopDecayer::_alpha, false, false, true, false, false);
}
double SMTopDecayer::me2(const int, const Particle & inpart,
const ParticleVector & decay,
MEOption meopt) const {
// 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);
ME(DecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
}
// 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);
}
}
// 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() {
DecayIntegrator::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 DecayIntegrator 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";
}
DecayIntegrator::dataBaseOutput(os,false);
if(header) os << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
void SMTopDecayer::doinitrun() {
DecayIntegrator::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(!isnan(width*GeV));
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(ShowerTreePtr tree, double & initial,
double & final) {
// check the outgoing particles
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
ShowerParticlePtr part[2];
unsigned int ix(0);
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
part[ix]=cit->first->progenitor();
++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=tree->incomingLines().begin()->first->progenitor()->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);
useMe();
}
void SMTopDecayer::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
// Get b and a and put them in particle vector ba in that order...
ParticleVector ba;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit)
ba.push_back(cit->first->copy());
PPtr temp;
if(abs(ba[0]->id())!=5) swap(ba[0],ba[1]);
// Get the starting scales for the showers $\tilde{\kappa}_{b}$
// $\tilde{\kappa}_{c}$. These are needed in order to know the boundary
// of the dead zone.
double ktb(0.),ktc(0.);
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cjt;
for(cjt = tree->incomingLines().begin();
cjt!= tree->incomingLines().end();++cjt) {
if(abs(cjt->first->progenitor()->id())==6)
ktb=sqr(cjt->first->progenitor()->evolutionScale()/_mt);
}
for(cit = tree->outgoingLines().begin();
cit!= tree->outgoingLines().end();++cit) {
if(abs(cit->first->progenitor()->id())==5)
ktc=sqr(cit->first->progenitor()->evolutionScale()/_mt);
}
if (ktb<=0.||ktc<=0.) {
throw Exception()
<< "SMTopDecayer::applyHardMatrixElementCorrection()"
<< " did not set ktb,ktc"
<< Exception::abortnow;
}
_ktb = ktb;
_ktc = ktc;
// 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;
// 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;
// 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,newa;
newc = ba[0]->data().produceParticle(newfs[0]);
newa = ba[1];
newa->set5Momentum(newfs[1]);
// set the colour lines
ColinePtr col;
if(ba[0]->id()>0) {
col=ba[0]->colourLine();
col->addColoured(newg);
newg->colourNeighbour(newc);
}
else {
col=ba[0]->antiColourLine();
col->addAntiColoured(newg);
newg->antiColourNeighbour(newc);
}
// change the existing quark and antiquark
PPtr orig;
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(cit->first->progenitor()->id()==newc->id()) {
// remove old particles from colour line
if(newc->id()>0) {
col->removeColoured(cit->first->copy());
col->removeColoured(cit->first->progenitor());
}
else {
col->removeAntiColoured(cit->first->copy());
col->removeAntiColoured(cit->first->progenitor());
}
// insert new particles
cit->first->copy(newc);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newc,2,true)));
cit->first->progenitor(sp);
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(false);
orig=cit->first->original();
}
else {
cit->first->copy(newa);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newa,2,true)));
map<tShowerTreePtr,pair<tShowerProgenitorPtr,
tShowerParticlePtr> >::const_iterator tit;
for(tit = tree->treelinks().begin();
tit != tree->treelinks().end();++tit) {
if(tit->second.first && tit->second.second==cit->first->progenitor())
break;
}
cit->first->progenitor(sp);
if(tit!=tree->treelinks().end())
tree->updateLink(tit->first,make_pair(cit->first,sp));
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(true);
}
}
// Add the gluon to the shower:
ShowerParticlePtr sg =new_ptr(ShowerParticle(*newg,2,true));
ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
gluon->perturbative(false);
tree->outgoingLines().insert(make_pair(gluon,sg));
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;
}
tree->hardMatrixElementCorrection(true);
}
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
*(_alpha->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
- unsigned int id[2]={abs(initial->progenitor()->id()),abs(parent->id())};
+ 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->setEvolutionScale(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->setEvolutionScale(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 not vetoed reset max
if(!veto) initial->highestpT(pt);
// if vetoing reset the scale
if(veto) parent->setEvolutionScale(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(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(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(isnan(output)||isinf(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;
}
diff --git a/MatrixElement/DIS/DISBase.cc b/MatrixElement/DIS/DISBase.cc
--- a/MatrixElement/DIS/DISBase.cc
+++ b/MatrixElement/DIS/DISBase.cc
@@ -1,1371 +1,1370 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DISBase class.
//
#include "DISBase.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig++/Utilities/Maths.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "Herwig++/PDT/StandardMatchers.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include <numeric>
using namespace Herwig;
using namespace ThePEG::Helicity;
namespace {
using namespace Herwig;
using namespace ThePEG::Helicity;
void debuggingMatrixElement(bool BGF,const Lorentz5Momentum & pin,
const Lorentz5Momentum & p1,
const Lorentz5Momentum & p2,
tcPDPtr gluon,
const Lorentz5Momentum & pl1,
const Lorentz5Momentum & pl2,
const Lorentz5Momentum & pq1,
const Lorentz5Momentum & pq2,
tcPDPtr lepton1,tcPDPtr lepton2,
tcPDPtr quark1 ,tcPDPtr quark2,
Energy2 Q2,double phi, double x2, double x3,
double xperp, double zp, double xp,
const vector<double> & azicoeff,
bool normalize) {
tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>
(CurrentGenerator::current().standardModel());
assert(hwsm);
vector<AbstractFFVVertexPtr> weakVertex;
vector<PDPtr> bosons;
AbstractFFVVertexPtr strongVertex = hwsm->vertexFFG();
if(lepton1->id()==lepton2->id()) {
weakVertex.push_back(hwsm->vertexFFZ());
bosons.push_back(hwsm->getParticleData(ParticleID::Z0));
weakVertex.push_back(hwsm->vertexFFP());
bosons.push_back(hwsm->getParticleData(ParticleID::gamma));
}
else {
weakVertex.push_back(hwsm->vertexFFW());
bosons.push_back(hwsm->getParticleData(ParticleID::Wplus));
}
if(!BGF) {
SpinorWaveFunction l1,q1,qp1;
SpinorBarWaveFunction l2,q2,qp2;
VectorWaveFunction gl(p2,gluon,outgoing);
if(lepton1->id()>0) {
l1 = SpinorWaveFunction (pl1,lepton1,incoming);
l2 = SpinorBarWaveFunction(pl2,lepton2,outgoing);
}
else {
l1 = SpinorWaveFunction (pl2,lepton2,outgoing);
l2 = SpinorBarWaveFunction(pl1,lepton1,incoming);
}
if(quark1->id()>0) {
q1 = SpinorWaveFunction (pq1,quark1,incoming);
q2 = SpinorBarWaveFunction(pq2,quark2,outgoing);
qp1 = SpinorWaveFunction (pin,quark1,incoming);
qp2 = SpinorBarWaveFunction(p1 ,quark2,outgoing);
}
else {
q1 = SpinorWaveFunction (pq2,quark2,outgoing);
q2 = SpinorBarWaveFunction(pq1,quark1,incoming);
qp1 = SpinorWaveFunction (p1 ,quark2,outgoing);
qp2 = SpinorBarWaveFunction(pin,quark1,incoming);
}
double lome(0.),realme(0.);
for(unsigned int lhel1=0;lhel1<2;++lhel1) {
l1.reset(lhel1);
for(unsigned int lhel2=0;lhel2<2;++lhel2) {
l2.reset(lhel2);
for(unsigned int qhel1=0;qhel1<2;++qhel1) {
q1.reset(qhel1);
qp1.reset(qhel1);
for(unsigned int qhel2=0;qhel2<2;++qhel2) {
q2.reset(qhel2);
qp2.reset(qhel2);
// leading order matrix element
Complex diagLO(0.);
for(unsigned int ix=0;ix<weakVertex.size();++ix) {
VectorWaveFunction inter =
weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
diagLO += weakVertex[ix]->evaluate(Q2,q1,q2,inter);
}
lome += norm(diagLO);
// real emission matrix element
for(unsigned int ghel=0;ghel<2;++ghel) {
gl.reset(2*ghel);
Complex diagReal(0.);
for(unsigned int ix=0;ix<weakVertex.size();++ix) {
VectorWaveFunction inter =
weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
SpinorWaveFunction off1 =
strongVertex->evaluate(Q2,5,qp1.particle(),qp1,gl);
Complex diag1 = weakVertex[ix]->evaluate(Q2,off1,qp2,inter);
SpinorBarWaveFunction off2 =
strongVertex->evaluate(Q2,5,qp2.particle(),qp2,gl);
Complex diag2 = weakVertex[ix]->evaluate(Q2,qp1,off2,inter);
diagReal += diag1+diag2;
}
realme += norm(diagReal);
}
}
}
}
}
double test1 = realme/lome/hwsm->alphaS(Q2)*Q2*UnitRemoval::InvE2;
double cphi(cos(phi));
double test2;
if(normalize) {
test2 = 8.*Constants::pi/(1.-xp)/(1.-zp)*
(azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi))*
(1.+sqr(xp)*(sqr(x2)+1.5*sqr(xperp)));
}
else {
test2 = 8.*Constants::pi/(1.-xp)/(1.-zp)*
(azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi));
}
cerr << "testing RATIO A " << test1/test2 << "\n";
}
else {
SpinorWaveFunction l1,q1,qp1;
SpinorBarWaveFunction l2,q2,qp2;
VectorWaveFunction gl(pin,gluon,incoming);
if(lepton1->id()>0) {
l1 = SpinorWaveFunction (pl1,lepton1,incoming);
l2 = SpinorBarWaveFunction(pl2,lepton2,outgoing);
}
else {
l1 = SpinorWaveFunction (pl2,lepton2,outgoing);
l2 = SpinorBarWaveFunction(pl1,lepton1,incoming);
}
if(quark1->id()>0) {
q1 = SpinorWaveFunction (pq1,quark1 ,incoming);
q2 = SpinorBarWaveFunction(pq2,quark2 ,outgoing);
qp2 = SpinorBarWaveFunction(p1 ,quark2 ,outgoing);
qp1 = SpinorWaveFunction (p2 ,quark1->CC(),outgoing);
}
else {
q1 = SpinorWaveFunction (pq2,quark2 ,outgoing);
q2 = SpinorBarWaveFunction(pq1,quark1 ,incoming);
qp2 = SpinorBarWaveFunction(p2 ,quark1->CC(),outgoing);
qp1 = SpinorWaveFunction (p1 ,quark2 ,outgoing);
}
double lome(0.),realme(0.);
for(unsigned int lhel1=0;lhel1<2;++lhel1) {
l1.reset(lhel1);
for(unsigned int lhel2=0;lhel2<2;++lhel2) {
l2.reset(lhel2);
for(unsigned int qhel1=0;qhel1<2;++qhel1) {
q1.reset(qhel1);
qp1.reset(qhel1);
for(unsigned int qhel2=0;qhel2<2;++qhel2) {
q2.reset(qhel2);
qp2.reset(qhel2);
// leading order matrix element
Complex diagLO(0.);
for(unsigned int ix=0;ix<weakVertex.size();++ix) {
VectorWaveFunction inter =
weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
diagLO += weakVertex[ix]->evaluate(Q2,q1,q2,inter);
}
lome += norm(diagLO);
// real emission matrix element
for(unsigned int ghel=0;ghel<2;++ghel) {
gl.reset(2*ghel);
Complex diagReal(0.);
for(unsigned int ix=0;ix<weakVertex.size();++ix) {
VectorWaveFunction inter =
weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
SpinorWaveFunction off1 =
strongVertex->evaluate(Q2,5,qp1.particle(),qp1,gl);
Complex diag1 = weakVertex[ix]->evaluate(Q2,off1,qp2,inter);
SpinorBarWaveFunction off2 =
strongVertex->evaluate(Q2,5,qp2.particle(),qp2,gl);
Complex diag2 = weakVertex[ix]->evaluate(Q2,qp1,off2,inter);
diagReal += diag1+diag2;
}
realme += norm(diagReal);
}
}
}
}
}
double test1 = realme/lome/hwsm->alphaS(Q2)*Q2*UnitRemoval::InvE2;
double cphi(cos(phi));
double test2;
if(normalize) {
test2 = 8.*Constants::pi/zp/(1.-zp)*
(azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi))*
sqr(xp)*(sqr(x3)+sqr(x2)+3.*sqr(xperp));
}
else {
test2 = 8.*Constants::pi/zp/(1.-zp)*
(azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi));
}
cerr << "testing RATIO B " << test1/test2 << "\n";
}
}
}
DISBase::DISBase() : initial_(6.), final_(3.),
procProb_(0.35),
comptonInt_(0.), bgfInt_(0.),
comptonWeight_(50.), BGFWeight_(150.),
pTmin_(0.1*GeV),
scaleOpt_(1), muF_(100.*GeV), scaleFact_(1.),
contrib_(0), power_(0.1)
{}
DISBase::~DISBase() {}
void DISBase::persistentOutput(PersistentOStream & os) const {
os << comptonInt_ << bgfInt_ << procProb_ << initial_ << final_ << alpha_
<< ounit(pTmin_,GeV) << comptonWeight_ << BGFWeight_ << gluon_
<< ounit(muF_,GeV) << scaleFact_ << scaleOpt_ << contrib_<< power_;
}
void DISBase::persistentInput(PersistentIStream & is, int) {
is >> comptonInt_ >> bgfInt_ >> procProb_ >> initial_ >> final_ >> alpha_
>> iunit(pTmin_,GeV) >> comptonWeight_ >> BGFWeight_ >> gluon_
>> iunit(muF_,GeV) >> scaleFact_ >> scaleOpt_ >> contrib_ >> power_;
}
AbstractClassDescription<DISBase> DISBase::initDISBase;
// Definition of the static class description member.
void DISBase::Init() {
static ClassDocumentation<DISBase> documentation
("The DISBase class provides the base class for the "
"implementation of DIS type processes including the "
"hard corrections in either the old-fashioned matrix "
"element correction of POWHEG approaches");
static Parameter<DISBase,double> interfaceProcessProbability
("ProcessProbability",
"The probabilty of the QCD compton process for the process selection",
&DISBase::procProb_, 0.3, 0.0, 1.,
false, false, Interface::limited);
static Reference<DISBase,ShowerAlpha> interfaceCoupling
("Coupling",
"Pointer to the object to calculate the coupling for the correction",
&DISBase::alpha_, false, false, true, false, false);
static Parameter<DISBase,Energy> interfacepTMin
("pTMin",
"The minimum pT",
&DISBase::pTmin_, GeV, 1.*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<DISBase,double> interfaceComptonWeight
("ComptonWeight",
"Weight for the overestimate ofthe compton channel",
&DISBase::comptonWeight_, 50.0, 0.0, 100.0,
false, false, Interface::limited);
static Parameter<DISBase,double> interfaceBGFWeight
("BGFWeight",
"Weight for the overestimate of the BGF channel",
&DISBase::BGFWeight_, 100.0, 0.0, 1000.0,
false, false, Interface::limited);
static Switch<DISBase,unsigned int> interfaceContribution
("Contribution",
"Which contributions to the cross section to include",
&DISBase::contrib_, 0, false, false);
static SwitchOption interfaceContributionLeadingOrder
(interfaceContribution,
"LeadingOrder",
"Just generate the leading order cross section",
0);
static SwitchOption interfaceContributionPositiveNLO
(interfaceContribution,
"PositiveNLO",
"Generate the positive contribution to the full NLO cross section",
1);
static SwitchOption interfaceContributionNegativeNLO
(interfaceContribution,
"NegativeNLO",
"Generate the negative contribution to the full NLO cross section",
2);
static Switch<DISBase,unsigned int> interfaceScaleOption
("ScaleOption",
"Option for the choice of factorization (and renormalization) scale",
&DISBase::scaleOpt_, 1, false, false);
static SwitchOption interfaceDynamic
(interfaceScaleOption,
"Dynamic",
"Dynamic factorization scale equal to the current sqrt(sHat())",
1);
static SwitchOption interfaceFixed
(interfaceScaleOption,
"Fixed",
"Use a fixed factorization scale set with FactorizationScaleValue",
2);
static Parameter<DISBase,Energy> interfaceFactorizationScale
("FactorizationScale",
"Value to use in the event of a fixed factorization scale",
&DISBase::muF_, GeV, 100.0*GeV, 1.0*GeV, 500.0*GeV,
true, false, Interface::limited);
static Parameter<DISBase,double> interfaceScaleFactor
("ScaleFactor",
"The factor used before Q2 if using a running scale",
&DISBase::scaleFact_, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<DISBase,double> interfaceSamplingPower
("SamplingPower",
"Power for the sampling of xp",
&DISBase::power_, 0.6, 0.0, 1.,
false, false, Interface::limited);
}
void DISBase::doinit() {
HwMEBase::doinit();
// integrals of me over phase space
double r5=sqrt(5.),darg((r5-1.)/(r5+1.)),ath(0.5*log((1.+1./r5)/(1.-1./r5)));
comptonInt_ = 2.*(-21./20.-6./(5.*r5)*ath+sqr(Constants::pi)/3.
-2.*Math::ReLi2(1.-darg)-2.*Math::ReLi2(1.-1./darg));
bgfInt_ = 121./9.-56./r5*ath;
// extract the gluon ParticleData objects
gluon_ = getParticleData(ParticleID::g);
}
void DISBase::initializeMECorrection(ShowerTreePtr tree, double & initial,
double & final) {
initial = initial_;
final = final_;
// incoming particles
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
partons_[0] = cit->first->progenitor()->dataPtr();
pq_[0] = cit->first->progenitor()->momentum();
}
else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
leptons_[0] = cit->first->progenitor()->dataPtr();
pl_[0] = cit->first->progenitor()->momentum();
}
}
// outgoing particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
partons_[1] = cit->first->progenitor()->dataPtr();
pq_[1] = cit->first->progenitor()->momentum();
}
else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
leptons_[1] = cit->first->progenitor()->dataPtr();
pl_[1] = cit->first->progenitor()->momentum();
}
}
// extract the born variables
q_ =pl_[0]-pl_[1];
q2_ = -q_.m2();
double yB = (q_*pq_[0])/(pl_[0]*pq_[0]);
l_ = 2./yB-1.;
// calculate the A coefficient for the correlations
acoeff_ = A(leptons_[0],leptons_[1],
partons_[0],partons_[1],q2_);
}
void DISBase::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
static const double eps=1e-6;
// find the incoming and outgoing quarks and leptons
ShowerParticlePtr quark[2],lepton[2];
PPtr hadron;
// incoming particles
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
hadron = cit->first->original()->parents()[0];
quark [0] = cit->first->progenitor();
beam_ = cit->first->beam();
}
else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
lepton[0] = cit->first->progenitor();
}
}
pdf_ = beam_->pdf();
assert(beam_&&pdf_&&quark[0]&&lepton[0]);
// outgoing particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(QuarkMatcher::Check(cit->first->progenitor()->data()))
quark [1] = cit->first->progenitor();
else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
lepton[1] = cit->first->progenitor();
}
}
// momentum fraction
assert(quark[1]&&lepton[1]);
xB_ = quark[0]->x();
// calculate the matrix element
vector<double> azicoeff;
// select the type of process
bool BGF = UseRandom::rnd()>procProb_;
double xp,zp,wgt,x1,x2,x3,xperp;
// generate a QCD compton process
if(!BGF) {
wgt = generateComptonPoint(xp,zp);
if(xp<eps) return;
// common pieces
Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
wgt *= 2./3./Constants::pi*alpha_->value(scale)/procProb_;
// PDF piece
wgt *= pdf_->xfx(beam_,quark[0]->dataPtr(),scale,xB_/xp)/
pdf_->xfx(beam_,quark[0]->dataPtr(),q2_ ,xB_);
// other bits
xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
x1 = -1./xp;
x2 = 1.-(1.-zp)/xp;
x3 = 2.+x1-x2;
// matrix element pieces
azicoeff = ComptonME(xp,x2,xperp,true);
}
// generate a BGF process
else {
wgt = generateBGFPoint(xp,zp);
if(xp<eps) return;
// common pieces
Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1);
wgt *= 0.25/Constants::pi*alpha_->value(scale)/(1.-procProb_);
// PDF piece
wgt *= pdf_->xfx(beam_,gluon_ ,scale,xB_/xp)/
pdf_->xfx(beam_,quark[0]->dataPtr(),q2_ ,xB_);
// other bits
xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
x1 = -1./xp;
x2 = 1.-(1.-zp)/xp;
x3 = 2.+x1-x2;
// matrix element pieces
azicoeff = BGFME(xp,x2,x3,xperp,true);
}
// compute the azimuthal average of the weight
wgt *= (azicoeff[0]+0.5*azicoeff[2]);
// decide whether or not to accept the weight
if(UseRandom::rnd()>wgt) return;
// if generate generate phi
unsigned int itry(0);
double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
double phiwgt,phi;
do {
phi = UseRandom::rnd()*Constants::twopi;
double cphi(cos(phi));
phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
++itry;
}
while (phimax*UseRandom::rnd() > phiwgt && itry<200);
if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
<< "::applyHardMatrixElementCorrection() to"
<< " generate phi" << Exception::eventerror;
// construct lorentz transform from lab to breit frame
Lorentz5Momentum phadron = hadron->momentum();
phadron.setMass(0.*GeV);
phadron.rescaleEnergy();
Lorentz5Momentum pcmf = phadron+0.5/xB_*q_;
pcmf.rescaleMass();
LorentzRotation rot(-pcmf.boostVector());
Lorentz5Momentum pbeam = rot*phadron;
Axis axis(pbeam.vect().unit());
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
Lorentz5Momentum pl = rot*pl_[0];
rot.rotateZ(-atan2(pl.y(),pl.x()));
pl_[0] *= rot;
pl_[1] *= rot;
pq_[0] *= rot;
pq_[1] *= rot;
// compute the new incoming and outgoing momenta
Energy Q(sqrt(q2_));
Lorentz5Momentum p1 = Lorentz5Momentum( 0.5*Q*xperp*cos(phi), 0.5*Q*xperp*sin(phi),
-0.5*Q*x2,0.*GeV,0.*GeV);
p1.rescaleEnergy();
Lorentz5Momentum p2 = Lorentz5Momentum(-0.5*Q*xperp*cos(phi),-0.5*Q*xperp*sin(phi),
-0.5*Q*x3,0.*GeV,0.*GeV);
p2.rescaleEnergy();
Lorentz5Momentum pin(0.*GeV,0.*GeV,-0.5*x1*Q,-0.5*x1*Q,0.*GeV);
// debuggingMatrixElement(BGF,pin,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
// lepton[0]->dataPtr(),lepton[1]->dataPtr(),
// quark [0]->dataPtr(),quark [1]->dataPtr(),
// q2_,phi,x2,x3,xperp,zp,xp,azicoeff,true);
// we need the Lorentz transform back to the lab
rot.invert();
// transform the momenta to lab frame
pin *= rot;
p1 *= rot;
p2 *= rot;
// test to ensure outgoing particles can be put on-shell
if(!BGF) {
if(p1.e()<quark[1]->dataPtr()->constituentMass()) return;
if(p2.e()<gluon_ ->constituentMass()) return;
}
else {
if(p1.e()<quark[1]->dataPtr() ->constituentMass()) return;
if(p2.e()<quark[0]->dataPtr()->CC()->constituentMass()) return;
}
// create the new particles and add to ShowerTree
bool isquark = quark[0]->colourLine();
if(!BGF) {
PPtr newin = new_ptr(Particle(*quark[0]));
newin->set5Momentum(pin);
PPtr newg = gluon_ ->produceParticle(p2 );
PPtr newout = quark[1]->dataPtr()->produceParticle(p1 );
ColinePtr col=isquark ?
quark[0]->colourLine() : quark[0]->antiColourLine();
ColinePtr newline=new_ptr(ColourLine());
// final-state emission
if(xp>zp) {
col->removeColoured(newout,!isquark);
col->addColoured(newin,!isquark);
col->addColoured(newg,!isquark);
newline->addColoured(newg,isquark);
newline->addColoured(newout,!isquark);
}
// initial-state emission
else {
col->removeColoured(newin ,!isquark);
col->addColoured(newout,!isquark);
col->addColoured(newg,isquark);
newline->addColoured(newg,!isquark);
newline->addColoured(newin,!isquark);
}
PPtr orig;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
if(cit->first->progenitor()!=quark[0]) continue;
// remove old particles from colour line
col->removeColoured(cit->first->copy(),!isquark);
col->removeColoured(cit->first->progenitor(),!isquark);
// insert new particles
cit->first->copy(newin);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
cit->first->progenitor(sp);
tree->incomingLines()[cit->first]=sp;
sp->x(xB_/xp);
cit->first->perturbative(xp>zp);
if(xp<=zp) orig=cit->first->original();
}
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(cit->first->progenitor()!=quark[1]) continue;
// remove old particles from colour line
col->removeColoured(cit->first->copy(),!isquark);
col->removeColoured(cit->first->progenitor(),!isquark);
// insert new particles
cit->first->copy(newout);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
cit->first->progenitor(sp);
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(xp<=zp);
if(xp>zp) orig=cit->first->original();
}
assert(orig);
// add the gluon
ShowerParticlePtr sg=new_ptr(ShowerParticle(*newg,1,true));
ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
gluon->perturbative(false);
tree->outgoingLines().insert(make_pair(gluon,sg));
tree->hardMatrixElementCorrection(true);
}
else {
PPtr newin = gluon_ ->produceParticle(pin);
PPtr newqbar = quark[0]->dataPtr()->CC()->produceParticle(p2 );
PPtr newout = quark[1]->dataPtr() ->produceParticle(p1 );
ColinePtr col=isquark ? quark[0]->colourLine() : quark[0]->antiColourLine();
ColinePtr newline=new_ptr(ColourLine());
col ->addColoured(newin ,!isquark);
newline->addColoured(newin , isquark);
col ->addColoured(newout ,!isquark);
newline->addColoured(newqbar, isquark);
PPtr orig;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
if(cit->first->progenitor()!=quark[0]) continue;
// remove old particles from colour line
col->removeColoured(cit->first->copy(),!isquark);
col->removeColoured(cit->first->progenitor(),!isquark);
// insert new particles
cit->first->copy(newin);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
cit->first->progenitor(sp);
tree->incomingLines()[cit->first]=sp;
sp->x(xB_/xp);
cit->first->perturbative(false);
orig=cit->first->original();
}
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(cit->first->progenitor()!=quark[1]) continue;
// remove old particles from colour line
col->removeColoured(cit->first->copy(),!isquark);
col->removeColoured(cit->first->progenitor(),!isquark);
// insert new particles
cit->first->copy(newout);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
cit->first->progenitor(sp);
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(true);
}
assert(orig);
// add the (anti)quark
ShowerParticlePtr sqbar=new_ptr(ShowerParticle(*newqbar,1,true));
ShowerProgenitorPtr qbar=new_ptr(ShowerProgenitor(orig,newqbar,sqbar));
qbar->perturbative(false);
tree->outgoingLines().insert(make_pair(qbar,sqbar));
tree->hardMatrixElementCorrection(true);
}
}
bool DISBase::softMatrixElementVeto(ShowerProgenitorPtr initial,
ShowerParticlePtr parent, Branching br) {
bool veto = !UseRandom::rndbool(parent->isFinalState() ? 1./final_ : 1./initial_);
// check if me correction should be applied
long id[2]={initial->id(),parent->id()};
if(id[0]!=id[1]||id[1]==ParticleID::g) return veto;
// get the pT
Energy pT=br.kinematics->pT();
// check if hardest so far
if(pT<initial->highestpT()) return veto;
double kappa(sqr(br.kinematics->scale())/q2_),z(br.kinematics->z());
double zk((1.-z)*kappa);
// final-state
double wgt(0.);
if(parent->isFinalState()) {
double zp=z,xp=1./(1.+z*zk);
double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
double x2 = 1.-(1.-zp)/xp;
vector<double> azicoeff = ComptonME(xp,x2,xperp,false);
wgt = (azicoeff[0]+0.5*azicoeff[2])*xp/(1.+sqr(z))/final_;
if(wgt<.0||wgt>1.) {
ostringstream wstring;
wstring << "Soft ME correction weight too large or "
<< "negative for FSR in DISBase::"
<< "softMatrixElementVeto() soft weight "
<< " xp = " << xp << " zp = " << zp
<< " weight = " << wgt << "\n";
generator()->logWarning( Exception(wstring.str(),
Exception::warning) );
}
}
else {
double xp = 2.*z/(1.+zk+sqrt(sqr(1.+zk)-4.*z*zk));
double zp = 0.5* (1.-zk+sqrt(sqr(1.+zk)-4.*z*zk));
double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
double x1 = -1./xp, x2 = 1.-(1.-zp)/xp, x3 = 2.+x1-x2;
// compton
if(br.ids[0]!=ParticleID::g) {
vector<double> azicoeff = ComptonME(xp,x2,xperp,false);
wgt = (azicoeff[0]+0.5*azicoeff[2])*xp*(1.-z)/(1.-xp)/(1.+sqr(z))/
(1.-zp+xp-2.*xp*(1.-zp));
}
// BGF
else {
vector<double> azicoeff = BGFME(xp,x2,x3,xperp,true);
wgt = (azicoeff[0]+0.5*azicoeff[2])*xp/(1.-zp+xp-2.*xp*(1.-zp))/(sqr(z)+sqr(1.-z));
}
wgt /=initial_;
if(wgt<.0||wgt>1.) {
ostringstream wstring;
wstring << "Soft ME correction weight too large or "
<< "negative for ISR in DISBase::"
<< "softMatrixElementVeto() soft weight "
<< " xp = " << xp << " zp = " << zp
<< " weight = " << wgt << "\n";
generator()->logWarning( Exception(wstring.str(),
Exception::warning) );
}
}
// if not vetoed
if(UseRandom::rndbool(wgt)) {
initial->highestpT(pT);
return false;
}
// otherwise
parent->setEvolutionScale(br.kinematics->scale());
return true;
}
double DISBase::generateComptonPoint(double &xp, double & zp) {
static const double maxwgt = 1.;
double wgt;
do {
xp = UseRandom::rnd();
double zpmin = xp, zpmax = 1./(1.+xp*(1.-xp));
zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
if(UseRandom::rndbool()) swap(xp,zp);
double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp,x2=1.-(1.-zp)/xp;
wgt *= 2.*(1.+sqr(xp)*(sqr(x2)+1.5*xperp2))/(1.-xp)/(1.-zp);
if(wgt>maxwgt) {
ostringstream wstring;
wstring << "DISBase::generateComptonPoint "
<< "Weight greater than maximum "
<< "wgt = " << wgt << " maxwgt = 1\n";
generator()->logWarning( Exception(wstring.str(),
Exception::warning) );
}
}
while(wgt<UseRandom::rnd()*maxwgt);
return comptonInt_;
}
double DISBase::generateBGFPoint(double &xp, double & zp) {
static const double maxwgt = 25.;
double wgt;
do {
xp = UseRandom::rnd();
double zpmax = 1./(1.+xp*(1.-xp)), zpmin = 1.-zpmax;
zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
double x1 = -1./xp;
double x2 = 1.-(1.-zp)/xp;
double x3 = 2.+x1-x2;
double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
wgt *= sqr(xp)/(1.-zp)*(sqr(x3)+sqr(x2)+3.*xperp2);
if(wgt>maxwgt) {
ostringstream wstring;
wstring << "DISBase::generateBGFPoint "
<< "Weight greater than maximum "
<< "wgt = " << wgt << " maxwgt = 1\n";
generator()->logWarning( Exception(wstring.str(),
Exception::warning) );
}
}
while(wgt<UseRandom::rnd()*maxwgt);
return bgfInt_;
// static const double maxwgt = 2.,npow=0.34,ac=1.0;
// double wgt;
// do {
// double rho = UseRandom::rnd();
// xp = 1.-pow(rho,1./(1.-npow));
// wgt = (sqr(xp)+ac+sqr(1.-xp));
// if(wgt>1.+ac) cerr << "testing violates BGF maxA " << wgt << "\n";
// }
// while(wgt<UseRandom::rnd()*(1.+ac));
// double xpwgt = -((6.-5.*npow+sqr(npow))*ac-3.*npow+sqr(npow)+4)
// /(sqr(npow)*(npow-6.)+11.*npow-6.);
// xpwgt *= pow(1.-xp,npow)/wgt;
// double xp2(sqr(xp)),lxp(log(xp)),xp4(sqr(xp2)),lxp1(log(1.-xp));
// double zpwgt = (2.*xp4*(lxp+lxp1-3.)+4.*xp*xp2*(3.-lxp-lxp1)
// +xp2*(-13.+lxp+lxp1)+xp*(+7.+lxp+lxp1)-lxp-lxp1-1.)/(1.+xp-xp2);
// do {
// double zpmax = 1./(1.+xp*(1.-xp)), zpmin = 1.-zpmax;
// zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
// wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
// double x1 = -1./xp;
// double x2 = 1.-(1.-zp)/xp;
// double x3 = 2.+x1-x2;
// double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
// wgt *= sqr(xp)/(1.-zp)*(sqr(x3)+sqr(x2)+3.*xperp2);
// if(wgt>maxwgt*zpwgt) cerr << "testing violates BGF maxB " << wgt/xpwgt << "\n";
// }
// while(wgt<UseRandom::rnd()*maxwgt);
// return zpwgt*xpwgt;
}
vector<double> DISBase::ComptonME(double xp, double x2, double xperp,
bool norm) {
vector<double> output(3,0.);
double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
double root = sqrt(sqr(l_)-1.);
output[0] = sqr(cos2)+acoeff_*cos2*l_+sqr(l_);
output[1] = -acoeff_*cos2*root*sin2-2.*l_*root*sin2;
output[2] = sqr(root)*sqr(sin2);
double lo(1+acoeff_*l_+sqr(l_));
double denom = norm ? 1.+sqr(xp)*(sqr(x2)+1.5*sqr(xperp)) : 1.;
double fact = sqr(xp)*(sqr(x2)+sqr(xperp))/lo;
for(unsigned int ix=0;ix<output.size();++ix)
output[ix] = ((ix==0 ? 1. : 0.) +fact*output[ix])/denom;
return output;
}
vector<double> DISBase::BGFME(double xp, double x2, double x3,
double xperp, bool norm) {
vector<double> output(3,0.);
double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
double fact2 = sqr(xp)*(sqr(x2)+sqr(xperp));
double cos3 = x3 /sqrt(sqr(x3)+sqr(xperp));
double sin3 = xperp/sqrt(sqr(x3)+sqr(xperp));
double fact3 = sqr(xp)*(sqr(x3)+sqr(xperp));
double root = sqrt(sqr(l_)-1.);
output[0] = fact2*(sqr(cos2)+acoeff_*cos2*l_+sqr(l_)) +
fact3*(sqr(cos3)-acoeff_*cos3*l_+sqr(l_));
output[1] = - fact2*(acoeff_*cos2*root*sin2+2.*l_*root*sin2)
- fact3*(acoeff_*cos3*root*sin3-2.*l_*root*sin3);
output[2] = fact2*(sqr(root)*sqr(sin2)) +
fact3*(sqr(root)*sqr(sin3));
double lo(1+acoeff_*l_+sqr(l_));
double denom = norm ? sqr(xp)*(sqr(x3)+sqr(x2)+3.*sqr(xperp))*lo : lo;
for(unsigned int ix=0;ix<output.size();++ix) output[ix] /= denom;
return output;
}
HardTreePtr DISBase::generateHardest(ShowerTreePtr tree) {
ShowerParticlePtr quark[2],lepton[2];
PPtr hadron;
// incoming particles
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
hadron = cit->first->original()->parents()[0];
quark [0] = cit->first->progenitor();
beam_ = cit->first->beam();
}
else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
lepton[0] = cit->first->progenitor();
leptons_[0] = lepton[0]->dataPtr();
}
}
pdf_=beam_->pdf();
assert(beam_&&pdf_&&quark[0]&&lepton[0]);
// outgoing particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(QuarkMatcher::Check(cit->first->progenitor()->data()))
quark [1] = cit->first->progenitor();
else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
lepton[1] = cit->first->progenitor();
leptons_[1] = lepton[1]->dataPtr();
}
}
assert(quark[1]&&lepton[1]);
// Particle data objects
for(unsigned int ix=0;ix<2;++ix) partons_[ix] = quark[ix]->dataPtr();
// extract the born variables
q_ =lepton[0]->momentum()-lepton[1]->momentum();
q2_ = -q_.m2();
xB_ = quark[0]->x();
double yB =
( q_*quark[0]->momentum())/
(lepton[0]->momentum()*quark[0]->momentum());
l_ = 2./yB-1.;
// construct lorentz transform from lab to breit frame
Lorentz5Momentum phadron = hadron->momentum();
phadron.setMass(0.*GeV);
phadron.rescaleRho();
Lorentz5Momentum pb = quark[0]->momentum();
- Lorentz5Momentum pbasis = phadron;
Axis axis(q_.vect().unit());
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
LorentzRotation rot_ = LorentzRotation();
if(axis.perp2()>1e-20) {
rot_.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
rot_.rotateX(Constants::pi);
}
if(abs(1.-q_.e()/q_.vect().mag())>1e-6) rot_.boostZ( q_.e()/q_.vect().mag());
pb *= rot_;
if(pb.perp2()/GeV2>1e-20) {
Boost trans = -1./pb.e()*pb.vect();
trans.setZ(0.);
rot_.boost(trans);
}
Lorentz5Momentum pl = rot_*lepton[0]->momentum();
rot_.rotateZ(-atan2(pl.y(),pl.x()));
// momenta of the particles
pl_[0]=rot_*lepton[0]->momentum();
pl_[1]=rot_*lepton[1]->momentum();
pq_[0]=rot_* quark[0]->momentum();
pq_[1]=rot_* quark[1]->momentum();
q_ *= rot_;
// coefficient for the matrix elements
acoeff_ = A(lepton[0]->dataPtr(),lepton[1]->dataPtr(),
quark [0]->dataPtr(),quark [1]->dataPtr(),q2_);
// generate a compton point
generateCompton();
generateBGF();
// no valid emission, return
if(pTCompton_<ZERO&&pTBGF_<ZERO) return HardTreePtr();
// type of emission, pick highest pT
bool isCompton=pTCompton_>pTBGF_;
// // find the sudakov for the branching
// SudakovPtr sudakov;
// // ISR
// if(ComptonISFS_||!isCompton) {
// BranchingList branchings=evolver()->splittingGenerator()->initialStateBranchings();
// long index = abs(partons_[0]->id());
// IdList br(3);
// if(isCompton) {
// br[0] = index;
// br[1] = index;
// br[2] = ParticleID::g;
// }
// else {
// br[0] = ParticleID::g;
// br[1] = abs(partons_[0]->id());
// br[2] = -abs(partons_[0]->id());
// }
// for(BranchingList::const_iterator cit = branchings.lower_bound(index);
// cit != branchings.upper_bound(index); ++cit ) {
// IdList ids = cit->second.second;
// if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
// sudakov=cit->second.first;
// break;
// }
// }
// }
// // FSR
// else {
// BranchingList branchings =
// evolver()->splittingGenerator()->finalStateBranchings();
// long index=abs(partons_[1]->id());
// for(BranchingList::const_iterator cit = branchings.lower_bound(index);
// cit != branchings.upper_bound(index); ++cit ) {
// IdList ids = cit->second.second;
// if(ids[0]==index&&ids[1]==index&&ids[2]==ParticleID::g) {
// sudakov = cit->second.first;
// break;
// }
// }
// }
// if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
// << "DISBase::generateHardest()"
// << Exception::runerror;
// add the leptons
vector<HardBranchingPtr> spaceBranchings,allBranchings;
spaceBranchings.push_back(new_ptr(HardBranching(lepton[0],SudakovPtr(),
HardBranchingPtr(),
HardBranching::Incoming)));
allBranchings.push_back(spaceBranchings.back());
allBranchings.push_back(new_ptr(HardBranching(lepton[1],SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
// compton hardest
if(isCompton) {
rot_.invert();
for(unsigned int ix=0;ix<ComptonMomenta_.size();++ix) {
ComptonMomenta_[ix].transform(rot_);
}
ShowerParticlePtr newqout (new_ptr(ShowerParticle(partons_[1],true)));
newqout->set5Momentum(ComptonMomenta_[1]);
ShowerParticlePtr newg(new_ptr(ShowerParticle(gluon_,true)));
newg->set5Momentum(ComptonMomenta_[2]);
ShowerParticlePtr newqin (new_ptr(ShowerParticle(partons_[0],false )));
newqin->set5Momentum(ComptonMomenta_[0]);
if(ComptonISFS_) {
ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[0],false)));
newspace->set5Momentum(ComptonMomenta_[0]-ComptonMomenta_[2]);
HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Incoming)));
HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),
spaceBranch,
HardBranching::Incoming)));
spaceBranch->addChild(offBranch);
HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),spaceBranch,
HardBranching::Outgoing)));
spaceBranch->addChild(g);
HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
spaceBranchings.push_back(spaceBranch);
allBranchings.push_back(offBranch);
allBranchings.push_back(outBranch);
ColinePtr newin(new_ptr(ColourLine())),newout(new_ptr(ColourLine()));
newin->addColoured(newqin,partons_[0]->id()<0);
newin->addColoured(newg ,partons_[0]->id()<0);
newout->addColoured(newspace,partons_[0]->id()<0);
newout->addColoured(newqout,partons_[1]->id()<0);
newout->addColoured(newg ,partons_[1]->id()>0);
ColinePtr newline(new_ptr(ColourLine()));
newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
newline->addColoured(newqout ,newspace->dataPtr()->iColour()!=PDT::Colour3);
}
else {
ShowerParticlePtr newtime(new_ptr(ShowerParticle(partons_[1],true)));
newtime->set5Momentum(ComptonMomenta_[1]+ComptonMomenta_[2]);
HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Incoming)));
HardBranchingPtr offBranch(new_ptr(HardBranching(newtime,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),offBranch,
HardBranching::Outgoing)));
HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),offBranch,
HardBranching::Outgoing)));
offBranch->addChild(outBranch);
offBranch->addChild(g);
spaceBranchings.push_back(spaceBranch);
allBranchings.push_back(spaceBranch);
allBranchings.push_back(offBranch);
ColinePtr newline(new_ptr(ColourLine()));
newline->addColoured(newqin ,newqin->dataPtr()->iColour()!=PDT::Colour3);
newline->addColoured(newtime,newqin->dataPtr()->iColour()!=PDT::Colour3);
}
}
// BGF hardest
else {
rot_.invert();
for(unsigned int ix=0;ix<BGFMomenta_.size();++ix) {
BGFMomenta_[ix].transform(rot_);
}
ShowerParticlePtr newq (new_ptr(ShowerParticle(partons_[1],true)));
newq->set5Momentum(BGFMomenta_[1]);
ShowerParticlePtr newqbar(new_ptr(ShowerParticle(partons_[0]->CC(),true)));
newqbar->set5Momentum(BGFMomenta_[2]);
ShowerParticlePtr newg (new_ptr(ShowerParticle(gluon_,false)));
newg->set5Momentum(BGFMomenta_[0]);
ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[0],false)));
newspace->set5Momentum(BGFMomenta_[0]-BGFMomenta_[2]);
HardBranchingPtr spaceBranch(new_ptr(HardBranching(newg,SudakovPtr(),HardBranchingPtr(),
HardBranching::Incoming)));
HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),spaceBranch,
HardBranching::Incoming)));
HardBranchingPtr qbar(new_ptr(HardBranching(newqbar,SudakovPtr(),spaceBranch,
HardBranching::Outgoing)));
spaceBranch->addChild(offBranch);
spaceBranch->addChild(qbar);
HardBranchingPtr outBranch(new_ptr(HardBranching(newq,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
spaceBranchings.push_back(spaceBranch);
allBranchings.push_back(offBranch);
allBranchings.push_back(outBranch);
ColinePtr newline(new_ptr(ColourLine()));
newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
newline->addColoured(newq ,newspace->dataPtr()->iColour()!=PDT::Colour3);
}
HardTreePtr newTree(new_ptr(HardTree(allBranchings,spaceBranchings,
ShowerInteraction::QCD)));
// Set the maximum pt for all other emissions and connect hard and shower tree
Energy pT = isCompton ? pTCompton_ : pTBGF_;
// incoming particles
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
// set maximum pT
if(QuarkMatcher::Check(cit->first->progenitor()->data()))
cit->first->maximumpT(pT);
for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
cjt!=newTree->branchings().end();++cjt) {
if(!(*cjt)->branchingParticle()->isFinalState()&&
(*cjt)->branchingParticle()->id()==cit->first->progenitor()->id()) {
newTree->connect(cit->first->progenitor(),*cjt);
tPPtr beam =cit->first->original();
if(!beam->parents().empty()) beam=beam->parents()[0];
(*cjt)->beam(beam);
HardBranchingPtr parent=(*cjt)->parent();
while(parent) {
parent->beam(beam);
parent=parent->parent();
};
}
}
}
// outgoing particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
// set maximum pT
if(QuarkMatcher::Check(cit->first->progenitor()->data()))
cit->first->maximumpT(pT);
for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
cjt!=newTree->branchings().end();++cjt) {
if((*cjt)->branchingParticle()->isFinalState()&&
(*cjt)->branchingParticle()->id()==cit->first->progenitor()->id()) {
newTree->connect(cit->first->progenitor(),*cjt);
}
}
}
// set the evolution partners and scales
ShowerParticleVector particles;
for(set<HardBranchingPtr>::iterator cit=newTree->branchings().begin();
cit!=newTree->branchings().end();++cit) {
particles.push_back((*cit)->branchingParticle());
}
// evolver()->showerModel()->partnerFinder()->
// setInitialEvolutionScales(particles,true,ShowerInteraction::QCD,true);
// // Calculate the shower variables:
// evolver()->showerModel()->kinematicsReconstructor()->
// deconstructHardJets(newTree,evolver(),ShowerInteraction::QCD);
// for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
// cjt!=newTree->branchings().end();++cjt) {
// if(cjt==newTree->branchings().begin()) {
// (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
// ++cjt;
// (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
// ++cjt;
// }
// // incoming
// if((**cjt).status()==HardBranching::Incoming) {
// quark[0]->set5Momentum((**cjt).showerMomentum());
// }
// // outgoing
// else {
// quark[1]->set5Momentum((**cjt).showerMomentum());
// }
// }
return newTree;
}
void DISBase::generateCompton() {
// maximum value of the xT
double xT = sqrt((1.-xB_)/xB_);
double xTMin = 2.*pTmin_/sqrt(q2_);
double zp;
// prefactor
double a = alpha_->overestimateValue()*comptonWeight_/Constants::twopi;
// loop to generate kinematics
double wgt(0.),xp(0.);
vector<double> azicoeff;
do {
wgt = 0.;
// intergration variables dxT/xT^3
xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
// zp
zp = UseRandom::rnd();
xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
// check allowed
if(xp<xB_||xp>1.) continue;
// phase-space piece of the weight
wgt = 8.*(1.-xp)*zp/comptonWeight_;
// PDF piece of the weight
Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
wgt *= pdf_->xfx(beam_,partons_[0],scale,xB_/xp)/
pdf_->xfx(beam_,partons_[0],q2_ ,xB_);
// me piece of the weight
double x2 = 1.-(1.-zp)/xp;
azicoeff = ComptonME(xp,x2,xT,false);
wgt *= 4./3.*alpha_->ratio(0.25*q2_*sqr(xT))*(azicoeff[0]+0.5*azicoeff[2]);
if(wgt>1.||wgt<0.) {
ostringstream wstring;
wstring << "DISBase::generateCompton() "
<< "Weight greater than one or less than zero"
<< "wgt = " << wgt << "\n";
generator()->logWarning( Exception(wstring.str(),
Exception::warning) );
}
}
while(xT>xTMin&&UseRandom::rnd()>wgt);
if(xT<=xTMin) {
pTCompton_=-GeV;
return;
}
// generate phi
unsigned int itry(0);
double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
double phiwgt,phi;
do {
phi = UseRandom::rnd()*Constants::twopi;
double cphi(cos(phi));
phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
++itry;
}
while (phimax*UseRandom::rnd() > phiwgt && itry<200);
if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
<< "::generateCompton() to"
<< " generate phi" << Exception::eventerror;
// momenta for the configuration
Energy Q(sqrt(q2_));
double x1 = -1./xp;
double x2 = 1.-(1.-zp)/xp;
double x3 = 2.+x1-x2;
Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
-0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
-0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
pTCompton_ = 0.5*Q*xT;
ComptonMomenta_.resize(3);
ComptonMomenta_[0] = p0;
ComptonMomenta_[1] = p1;
ComptonMomenta_[2] = p2;
ComptonISFS_ = zp>xp;
// debuggingMatrixElement(false,p0,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
// leptons_[0],leptons_[1],
// partons_[0],partons_[1],
// q2_,phi,x2,x3,xT,zp,xp,azicoeff,false);
}
void DISBase::generateBGF() {
// maximum value of the xT
double xT = (1.-xB_)/xB_;
double xTMin = 2.*max(pTmin_,pTCompton_)/sqrt(q2_);
double zp;
// prefactor
double a = alpha_->overestimateValue()*BGFWeight_/Constants::twopi;
// loop to generate kinematics
double wgt(0.),xp(0.);
vector<double> azicoeff;
do {
wgt = 0.;
// intergration variables dxT/xT^3
xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
// zp
zp = UseRandom::rnd();
xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
// check allowed
if(xp<xB_||xp>1.) continue;
// phase-space piece of the weight
wgt = 8.*sqr(1.-xp)*zp/BGFWeight_;
// PDF piece of the weight
Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
wgt *= pdf_->xfx(beam_,gluon_ ,scale,xB_/xp)/
pdf_->xfx(beam_,partons_[0],q2_ ,xB_);
// me piece of the weight
double x1 = -1./xp;
double x2 = 1.-(1.-zp)/xp;
double x3 = 2.+x1-x2;
azicoeff = BGFME(xp,x2,x3,xT,false);
wgt *= 0.5*alpha_->ratio(0.25*q2_*sqr(xT))*
(azicoeff[0]+0.5*azicoeff[2]);
if(wgt>1.||wgt<0.) {
ostringstream wstring;
wstring << "DISBase::generateBGF() "
<< "Weight greater than one or less than zero"
<< "wgt = " << wgt << "\n";
generator()->logWarning( Exception(wstring.str(),
Exception::warning) );
}
}
while(xT>xTMin&&UseRandom::rnd()>wgt);
if(xT<=xTMin) {
pTBGF_=-GeV;
return;
}
// generate phi
unsigned int itry(0);
double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
double phiwgt,phi;
do {
phi = UseRandom::rnd()*Constants::twopi;
double cphi(cos(phi));
phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
++itry;
}
while (phimax*UseRandom::rnd() > phiwgt && itry<200);
if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
<< "::generateBGF() to"
<< " generate phi" << Exception::eventerror;
// momenta for the configuration
Energy Q(sqrt(q2_));
double x1 = -1./xp;
double x2 = 1.-(1.-zp)/xp;
double x3 = 2.+x1-x2;
Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
-0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
-0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
pTBGF_=0.5*Q*xT;
BGFMomenta_.resize(3);
BGFMomenta_[0]=p0;
BGFMomenta_[1]=p1;
BGFMomenta_[2]=p2;
// debuggingMatrixElement(true,p0,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
// leptons_[0],leptons_[1],
// partons_[0],partons_[1],
// q2_,phi,x2,x3,xT,zp,xp,azicoeff,false);
}
int DISBase::nDim() const {
return HwMEBase::nDim() + (contrib_>0 ? 1 : 0 );
}
bool DISBase::generateKinematics(const double * r) {
// Born kinematics
if(!HwMEBase::generateKinematics(r)) return false;
if(contrib_!=0) {
// hadron and momentum fraction
if(HadronMatcher::Check(*lastParticles().first->dataPtr())) {
hadron_ = dynamic_ptr_cast<tcBeamPtr>(lastParticles().first->dataPtr());
xB_ = lastX1();
}
else {
hadron_ = dynamic_ptr_cast<tcBeamPtr>(lastParticles().second->dataPtr());
xB_ = lastX2();
}
// Q2
q2_ = -(meMomenta()[0]-meMomenta()[2]).m2();
// xp
int ndim=nDim();
double rhomin = pow(1.-xB_,1.-power_);
double rho = r[ndim-1]*rhomin;
xp_ = 1.-pow(rho,1./(1.-power_));
jac_ = rhomin/(1.-power_)*pow(1.-xp_,power_);
jacobian(jacobian()*jac_);
}
return true;
}
Energy2 DISBase::scale() const {
return scaleOpt_ == 1 ?
-sqr(scaleFact_)*tHat() : sqr(scaleFact_*muF_);
}
CrossSection DISBase::dSigHatDR() const {
return NLOWeight()*HwMEBase::dSigHatDR();
}
double DISBase::NLOWeight() const {
// If only leading order is required return 1:
if(contrib_==0) return 1.;
// scale and prefactors
Energy2 mu2(scale());
double aS = SM().alphaS(mu2);
double CFfact = 4./3.*aS/Constants::twopi;
double TRfact = 1./2.*aS/Constants::twopi;
// LO + dipole subtracted virtual + collinear quark bit with LO pdf
double virt = 1.+CFfact*(-4.5-1./3.*sqr(Constants::pi)+1.5*log(q2_/mu2/(1.-xB_))
+2.*log(1.-xB_)*log(q2_/mu2)+sqr(log(1.-xB_)));
virt /= jac_;
// PDF from leading-order
double loPDF = hadron_->pdf()->xfx(hadron_,mePartonData()[1],mu2,xB_)/xB_;
// NLO gluon PDF
tcPDPtr gluon = getParticleData(ParticleID::g);
double gPDF = hadron_->pdf()->xfx(hadron_,gluon,mu2,xB_/xp_)*xp_/xB_;
// NLO quark PDF
double qPDF = hadron_->pdf()->xfx(hadron_,mePartonData()[1],mu2,xB_/xp_)*xp_/xB_;
// collinear counterterms
// gluon
double collg =
TRfact/xp_*gPDF*(2.*xp_*(1.-xp_)+(sqr(xp_)+sqr(1.-xp_))*log((1.-xp_)*q2_/xp_/mu2));
// quark
double collq =
CFfact/xp_*qPDF*(1-xp_-2./(1.-xp_)*log(xp_)-(1.+xp_)*log((1.-xp_)/xp_*q2_/mu2))+
CFfact/xp_*(qPDF-xp_*loPDF)*(2./(1.-xp_)*log(q2_*(1.-xp_)/mu2)-1.5/(1.-xp_));
// calculate the A coefficient for the real pieces
double a(A(mePartonData()[0],mePartonData()[2],
mePartonData()[1],mePartonData()[3],q2_));
// cacluate lepton kinematic variables
Lorentz5Momentum q = meMomenta()[0]-meMomenta()[2];
double yB = (q*meMomenta()[1])/(meMomenta()[0]*meMomenta()[1]);
double l = 2./yB-1.;
// q -> qg term
double realq = CFfact/xp_/(1.+a*l+sqr(l))*qPDF/loPDF*
(2.+2.*sqr(l)-xp_+3.*xp_*sqr(l)+a*l*(2.*xp_+1.));
// g -> q qbar term
double realg =-TRfact/xp_/(1.+a*l+sqr(l))*gPDF/loPDF*
((1.+sqr(l)+2.*(1.-3.*sqr(l))*xp_*(1.-xp_))
+2.*a*l*(1.-2.*xp_*(1.-xp_)));
// return the full result
double wgt = virt+((collq+collg)/loPDF+realq+realg);
// double f2g = gPDF/xp_*TRfact*((sqr(1-xp_)+sqr(xp_))*log((1-xp_)/xp_)+
// 8*xp_*(1.-xp_)-1.);
// double f2q =
// loPDF/jac_*(1.+CFfact*(-1.5*log(1.-xB_)+sqr(log(1.-xB_))
// -sqr(Constants::pi)/3.-4.5))
// +qPDF *CFfact/xp_*(3.+2.*xp_-(1.+xp_)*log(1.-xp_)
// -(1.+sqr(xp_))/(1.-xp_)*log(xp_))
// +(qPDF-xp_*loPDF)*CFfact/xp_*(2.*log(1.-xp_)/(1.-xp_)-1.5/(1.-xp_));
// double wgt = (f2g+f2q)/loPDF;
return contrib_ == 1 ? max(0.,wgt) : max(0.,-wgt);
}
diff --git a/Models/Sextet/SextetFFSVertex.cc b/Models/Sextet/SextetFFSVertex.cc
--- a/Models/Sextet/SextetFFSVertex.cc
+++ b/Models/Sextet/SextetFFSVertex.cc
@@ -1,216 +1,216 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SextetFFSVertex class.
//
#include "SextetFFSVertex.h"
#include "SextetModel.h"
#include "SextetParticles.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"
using namespace Herwig;
IBPtr SextetFFSVertex::clone() const {
return new_ptr(*this);
}
IBPtr SextetFFSVertex::fullclone() const {
return new_ptr(*this);
}
void SextetFFSVertex::persistentOutput(PersistentOStream & os) const {
os << g1L_ << g1R_ << g1pR_ << g1ppR_ << g3L_;
}
void SextetFFSVertex::persistentInput(PersistentIStream & is, int) {
is >> g1L_ >> g1R_ >> g1pR_ >> g1ppR_ >> g3L_;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<SextetFFSVertex,Helicity::FFSVertex>
describeSextetFFSVertex("Herwig::SextetFFSVertex", "HwSextetModel.so");
void SextetFFSVertex::Init() {
static ClassDocumentation<SextetFFSVertex> documentation
("The SextetFFSVertex class implements the coupling of two "
"fermions to a scalar sextet particle.");
}
void SextetFFSVertex::doinit() {
orderInGs (0);
orderInGem(1);
SextetModelPtr model =
dynamic_ptr_cast<SextetModelPtr>(generator()->standardModel());
if(!model) throw Exception() << "Must be using the SextetModel"
<< " in SextetGSSVertex::doinit()"
<< Exception::runerror;
// extract the couplings
g1L_ = model->g1L();
g1R_ = model->g1R();
g1pR_ = model->g1pR();
g1ppR_ = model->g1ppR();
g3L_ = model->g3L();
// add the enabled particles
if(model->ScalarSingletY43Enabled()) {
for(long ix=0;ix<3;++ix) {
long iu = 2*ix + 2;
if(g1ppR_[ix]!=0.) {
addToList( iu, iu, ParticleID::ScalarDQSingletY43bar);
addToList( -iu, -iu, ParticleID::ScalarDQSingletY43);
}
}
}
if(model->ScalarSingletY13Enabled()) {
for(long ix=0;ix<3;++ix) {
long iu = 2*ix + 2;
long id = 2*ix + 1;
if(g1L_[ix]!=0. || g1R_[ix]!=0.) {
addToList( id, iu, ParticleID::ScalarDQSingletY13bar);
addToList( -id, -iu, ParticleID::ScalarDQSingletY13);
}
}
}
if(model->ScalarSingletY23Enabled()) {
for(long ix=0;ix<3;++ix) {
long id = 2*ix + 1;
if(g1pR_[ix]!=0. ) {
addToList( id, id,ParticleID::ScalarDQSingletY23bar);
addToList(-id,-id,ParticleID::ScalarDQSingletY23);
}
}
}
if(model->ScalarTripletY13Enabled()) {
for(long ix=0;ix<3;++ix) {
long iu = 2*ix + 2;
long id = 2*ix + 1;
if(g3L_[ix]!=0. ) {
addToList( iu, iu, ParticleID::ScalarDQTripletPbar);
addToList( -iu, -iu, ParticleID::ScalarDQTripletP);
addToList( iu, id, ParticleID::ScalarDQTriplet0bar);
addToList( -iu, -id, ParticleID::ScalarDQTriplet0);
addToList( id, id, ParticleID::ScalarDQTripletMbar);
addToList( -id, -id, ParticleID::ScalarDQTripletM);
}
}
}
Helicity::FFSVertex::doinit();
}
-void SextetFFSVertex::setCoupling(Energy2 q2,tcPDPtr part1,
- tcPDPtr part2,tcPDPtr part3) {
+void SextetFFSVertex::setCoupling(Energy2,tcPDPtr part1,
+ tcPDPtr part2,tcPDPtr part3) {
long q1ID=(abs(part1->id())), q2ID=(abs(part2->id())),
sDQID=(abs(part3->id()));
//check scalar diquark
assert( sDQID == ParticleID::ScalarDQSingletY43 ||
sDQID == ParticleID::ScalarDQSingletY13 ||
sDQID == ParticleID::ScalarDQSingletY23 ||
sDQID == ParticleID::ScalarDQTripletP ||
sDQID == ParticleID::ScalarDQTriplet0 ||
sDQID == ParticleID::ScalarDQTripletM);
//check quarks
assert(!(q1ID>6) && !(q2ID>6));
bool part1Up = (q1ID==2 || q1ID==4 || q1ID==6) ? true : false;
bool part2Up = (q2ID==2 || q2ID==4 || q2ID==6) ? true : false;
Complex cRight, cLeft, prefactor(1.);
if(sDQID==ParticleID::ScalarDQSingletY43){
//should both be up type
assert(part1Up && part2Up);
if(q1ID==2)
cRight=Complex(g1ppR_[0]);
else if(q1ID==4)
cRight=Complex(g1ppR_[1]);
else
cRight=Complex(g1ppR_[2]);
cLeft=Complex(0.);
}
if(sDQID==ParticleID::ScalarDQSingletY13){
//should be one up one down type
assert((part1Up && !part2Up) || (!part1Up && part2Up));
long upType;
if(part1Up)
upType=q1ID;
else
upType=q2ID;
if(upType==2){
cRight=Complex(g1R_[0]);
cLeft=Complex(2.*g1L_[0]);
}
else if(upType==4){
cRight=Complex(g1R_[1]);
cLeft=Complex(2.*g1L_[1]);
}
else
cRight=Complex(g1R_[2]);{
cLeft=Complex(2.*g1L_[2]);
}
}
if(sDQID==ParticleID::ScalarDQSingletY23){
//should both be down type
assert(!part1Up && !part2Up);
if(q1ID==1)
cRight=Complex(g1pR_[0]);
else if(q1ID==3)
cRight=Complex(g1pR_[1]);
else
cRight=Complex(g1pR_[2]);
cLeft=Complex(0.);
}
if(sDQID==ParticleID::ScalarDQTripletP){
//should both be up type
assert(part1Up && part2Up);
if(q1ID==2)
cLeft=Complex(g3L_[0]);
else if(q1ID==4)
cLeft=Complex(g3L_[1]);
else
cLeft=Complex(g3L_[2]);
cRight=Complex(0.);
}
if(sDQID==ParticleID::ScalarDQTriplet0){
//should both one up and down type
assert((part1Up && !part2Up) || (!part1Up && part2Up));
//possibly doesn't couple
long upType;
if(part1Up)
upType=q1ID;
else
upType=q2ID;
if(upType==2)
cLeft=Complex(g3L_[0]);
else if(upType==4)
cLeft=Complex(g3L_[1]);
else
cLeft=Complex(g3L_[2]);
cRight=Complex(0.);
}
if(sDQID==ParticleID::ScalarDQTripletM){
//should one both be down type
assert(!part1Up && !part2Up);
if(q1ID==1)
cLeft=Complex(g3L_[0]);
else if(q1ID==3)
cLeft=Complex(g3L_[1]);
else
cLeft=Complex(g3L_[2]);
cRight=Complex(0.);
prefactor=Complex(-1.);
}
left(cLeft);
right(cRight);
norm(prefactor);
}
diff --git a/Models/Sextet/SextetFFVVertex.cc b/Models/Sextet/SextetFFVVertex.cc
--- a/Models/Sextet/SextetFFVVertex.cc
+++ b/Models/Sextet/SextetFFVVertex.cc
@@ -1,171 +1,171 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SextetFFVVertex class.
//
#include "SextetFFVVertex.h"
#include "SextetModel.h"
#include "SextetParticles.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"
using namespace Herwig;
IBPtr SextetFFVVertex::clone() const {
return new_ptr(*this);
}
IBPtr SextetFFVVertex::fullclone() const {
return new_ptr(*this);
}
void SextetFFVVertex::persistentOutput(PersistentOStream & os) const {
os << g2_ << g2p_ ;
}
void SextetFFVVertex::persistentInput(PersistentIStream & is, int) {
is >> g2_ >> g2p_ ;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<SextetFFVVertex,Helicity::FFVVertex>
describeSextetFFVVertex("Herwig::SextetFFVVertex", "HwSextetModel.so");
void SextetFFVVertex::Init() {
static ClassDocumentation<SextetFFVVertex> documentation
("The SextetFFVVertex class implements the coupling of two "
"fermions to a scalar sextet particle.");
}
void SextetFFVVertex::doinit() {
orderInGs (0);
orderInGem(1);
SextetModelPtr model =
dynamic_ptr_cast<SextetModelPtr>(generator()->standardModel());
if(!model) throw Exception() << "Must be using the SextetModel"
<< " in SextetGSSVertex::doinit()"
<< Exception::runerror;
// extract the couplings
g2_ = model->g2();
g2p_= model->g2p();
// add the enabled particles
if(model->VectorDoubletY16Enabled()) {
for(long ix=0;ix<3;++ix) {
long iu = 2*ix + 2;
long id = 2*ix + 1;
if(g2_[ix]!=0.) {
addToList( -id, -iu, ParticleID::VectorDQY16P);
addToList( id, iu, ParticleID::VectorDQY16Pbar);
addToList( -id, -id, ParticleID::VectorDQY16M);
addToList( id, id, ParticleID::VectorDQY16Mbar);
}
}
}
if(model->VectorDoubletY56Enabled()) {
for(long ix=0;ix<3;++ix) {
long iu = 2*ix + 2;
long id = 2*ix + 1;
if(g2p_[ix]!=0.) {
addToList( -iu, -iu, ParticleID::VectorDQY56P);
addToList( iu, iu, ParticleID::VectorDQY56Pbar);
addToList( -id, -iu, ParticleID::VectorDQY56M);
addToList( id, iu, ParticleID::VectorDQY56Mbar);
}
}
}
Helicity::FFVVertex::doinit();
}
-void SextetFFVVertex::setCoupling(Energy2 q2,tcPDPtr part1,
- tcPDPtr part2,tcPDPtr part3) {
+void SextetFFVVertex::setCoupling(Energy2, tcPDPtr part1,
+ tcPDPtr part2,tcPDPtr part3) {
long q1ID=(abs(part1->id())), q2ID=(abs(part2->id())),
vDQID=(abs(part3->id()));
//check scalar diquark
assert( vDQID == ParticleID::VectorDQY16P ||
vDQID == ParticleID::VectorDQY16M ||
vDQID == ParticleID::VectorDQY56P ||
vDQID == ParticleID::VectorDQY56M);
//check quarks
assert(!(q1ID>6) && !(q2ID>6));
bool part1Up = (q1ID==2 || q1ID==4 || q1ID==6) ? true : false;
bool part2Up = (q2ID==2 || q2ID==4 || q2ID==6) ? true : false;
Complex cRight(1.,1.), cLeft(1.,1.), prefactor(1.,0.);
if(vDQID==ParticleID::VectorDQY16P){
//should be one up and down type
assert((!part1Up && part2Up) || (part1Up && !part2Up));
long upType;
if(part1Up)
upType=q1ID;
else
upType=q2ID;
if(upType==2)
cRight=Complex(g2_[0]);
else if(upType==4)
cRight=Complex(g2_[1]);
else
cRight=Complex(g2_[2]);
cLeft=Complex(0.);
}
if(vDQID==ParticleID::VectorDQY16M) {
//should be both be down type
assert(!part1Up && !part2Up);
if(q1ID==1)
cRight=Complex(g2_[0]);
else if(q1ID==2)
cRight=Complex(g2_[1]);
else
cRight=Complex(g2_[2]);
cLeft=Complex(0.);
}
if(vDQID==ParticleID::VectorDQY56P) {
//should both be up type
assert(part1Up && part2Up);
if(q1ID==2)
cRight=Complex(g2p_[0]);
else if(q1ID==4)
cRight=Complex(g2p_[1]);
else
cRight=Complex(g2p_[2]);
cLeft=Complex(0.);
}
if(vDQID==ParticleID::VectorDQY56M){
//should be one up and down type
assert((!part1Up && part2Up) || (part1Up && !part2Up));
long upType;
if(part1Up)
upType=q1ID;
else
upType=q2ID;
if(upType==2)
cRight=Complex(g2p_[0]);
else if(upType==4)
cRight=Complex(g2p_[1]);
else
cRight=Complex(g2p_[2]);
cLeft=Complex(0.);
}
left(cLeft);
right(cRight);
norm(prefactor);
}
diff --git a/Models/StandardModel/SMWWWWVertex.cc b/Models/StandardModel/SMWWWWVertex.cc
--- a/Models/StandardModel/SMWWWWVertex.cc
+++ b/Models/StandardModel/SMWWWWVertex.cc
@@ -1,164 +1,164 @@
// -*- C++ -*-
//
// SMWWWWVertex.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 SMWWWWVertex class.
//
#include "SMWWWWVertex.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
using namespace ThePEG;
SMWWWWVertex::SMWWWWVertex()
: _couplast(0.0), _q2last(sqr(Constants::MaxEnergy)),
_vfact(4,0.0), _sw2(0.), _cw2(0.) {
orderInGem(2);
orderInGs(0);
}
void SMWWWWVertex::doinit() {
// particles
addToList(24, -24, 24, -24);
addToList(23, 24, 23, -24);
addToList(22, 24, 22, -24);
addToList(22, 24, 23, -24);
VVVVVertex::doinit();
// couplings
_sw2 = sin2ThetaW();
_cw2 = 1.-_sw2;
double sw = sqrt(_sw2);
double cw = sqrt(_cw2);
_vfact[0] = -1./_sw2;
_vfact[1] = _cw2/_sw2;
_vfact[2] = 1.;
_vfact[3] = cw/sw;
// pointer for intermediate particles
_gamma = getParticleData(ThePEG::ParticleID::gamma);
_Z0 = getParticleData(ThePEG::ParticleID::Z0);
_wplus = getParticleData(ThePEG::ParticleID::Wplus);
_wminus = getParticleData(ThePEG::ParticleID::Wminus);
}
void SMWWWWVertex::persistentOutput(PersistentOStream & os) const {
os << _gamma << _Z0 << _wplus << _wminus
<< _vfact << _sw2 << _cw2;
}
void SMWWWWVertex::persistentInput(PersistentIStream & is, int) {
is >> _gamma >> _Z0 >> _wplus >> _wminus
>> _vfact >> _sw2 >> _cw2;
}
ClassDescription<SMWWWWVertex>SMWWWWVertex::initSMWWWWVertex;
// Definition of the static class description member.
void SMWWWWVertex::Init() {
static ClassDocumentation<SMWWWWVertex> documentation
("The SMWWWWVertex class is the implementation of the"
" Standard Model quartic electroweka gauge boson coupling.");
}
// couplings for the WWWW vertex
void SMWWWWVertex::setCoupling(Energy2 q2,tcPDPtr a,tcPDPtr b,
tcPDPtr c,tcPDPtr d) {
// id's of the particles
- int id[4]={a->id(),b->id(),c->id(),d->id()};
+ long id[4]={a->id(),b->id(),c->id(),d->id()};
// order the particles
int ngamma(0),nz(0);
int iorder[4];
for(int ix=0;ix<4;++ix) {
if (id[ix]==22) ++ngamma;
else if (id[ix]==23) ++nz;
}
// if photons or Z's
if(ngamma!=0 || nz!=0) {
int iy=0;
// put the photons first
for(int ix=0;iy<ngamma&&ix<4;++ix) {
if(id[ix]==22) {
iorder[iy]=ix;
++iy;
}
}
// then the Z bosons
for(int ix=0;iy<ngamma+nz&&ix<4;++ix) {
if(id[ix]==23) {
iorder[iy]=ix;
++iy;
}
}
// then the W+
for(int ix=0;iy<3&&ix<4;++ix) {
if(id[ix]==24) {
iorder[iy]=ix;
++iy;
}
}
assert(iy==3);
// finally the W-
for(int ix=0;iy<4&&ix<4;++ix) {
if(id[ix]==-24) {
iorder[iy]=ix;
++iy;
}
}
assert(iy==4);
}
else {
int iy=0;
// first the W+
for(int ix=0;iy<3&&ix<4;++ix) {
if(id[ix]==24) {
iorder[iy]=ix;
++iy;
}
}
assert(iy==2);
// finally the W-
for(int ix=0;iy<4&&ix<4;++ix) {
if(id[ix]==-24) {
iorder[iy]=ix;
++iy;
}
}
assert(iy==4);
setIntermediate(_gamma,_Z0,_sw2,_cw2);
}
setOrder(iorder[0],iorder[1],iorder[2],iorder[3]);
setType(2);
// first the overall normalisation
if(q2!=_q2last||_couplast==0.) {
_couplast = sqr(electroMagneticCoupling(q2));
_q2last=q2;
}
// id's of the first two particles
int ida(0),idb(0);
if(iorder[0]==0) ida = abs(a->id());
else if(iorder[0]==1) ida = abs(b->id());
else if(iorder[0]==2) ida = abs(c->id());
else if(iorder[0]==3) ida = abs(d->id());
if(iorder[1]==0) idb = abs(a->id());
else if(iorder[1]==1) idb = abs(b->id());
else if(iorder[1]==2) idb = abs(c->id());
else if(iorder[1]==3) idb = abs(d->id());
// WWWW coupling
if(ida==24) norm(_vfact[0]*_couplast);
// ZZWW coupling
else if(ida==23&&idb==23) norm(_vfact[1]*_couplast);
// gamma gamma WW coupling
else if(ida==22&&idb==22) norm(_couplast);
// gamma Z WW coupling
else if(ida==22&&idb==23) norm(_vfact[3]*_couplast);
else assert(false);
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:42 PM (1 d, 8 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806136
Default Alt Text
(114 KB)

Event Timeline