Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8724417
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
134 KB
Subscribers
None
View Options
diff --git a/MatrixElement/Lepton/MEee2gZ2qq.cc b/MatrixElement/Lepton/MEee2gZ2qq.cc
--- a/MatrixElement/Lepton/MEee2gZ2qq.cc
+++ b/MatrixElement/Lepton/MEee2gZ2qq.cc
@@ -1,1089 +1,1104 @@
// -*- C++ -*-
//
// MEee2gZ2qq.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 MEee2gZ2qq class.
//
#include "MEee2gZ2qq.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig++/MatrixElement/HardVertex.h"
#include "Herwig++/Shower/Base/Evolver.h"
#include "Herwig++/Shower/Base/KinematicsReconstructor.h"
#include "Herwig++/Shower/Base/PartnerFinder.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
#include <numeric>
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
const double MEee2gZ2qq::EPS_=0.00000001;
void MEee2gZ2qq::doinit() {
HwMEBase::doinit();
massOption(vector<unsigned int>(2,massopt_));
rescalingOption(3);
if(minflav_>maxflav_)
throw InitException() << "The minimum flavour " << minflav_
<< "must be lower the than maximum flavour " << maxflav_
<< " in MEee2gZ2qq::doinit() "
<< Exception::runerror;
// set the particle data objects
Z0_ = getParticleData(ParticleID::Z0);
gamma_ = getParticleData(ParticleID::gamma);
gluon_ = getParticleData(ParticleID::g);
// cast the SM pointer to the Herwig SM pointer
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(!hwsm) throw InitException()
<< "Wrong type of StandardModel object in "
<< "MEee2gZ2qq::doinit() the Herwig++ version must be used"
<< Exception::runerror;
FFZVertex_ = hwsm->vertexFFZ();
FFPVertex_ = hwsm->vertexFFP();
FFGVertex_ = hwsm->vertexFFG();
}
void MEee2gZ2qq::getDiagrams() const {
// specific the diagrams
tcPDPtr ep = getParticleData(ParticleID::eplus);
tcPDPtr em = getParticleData(ParticleID::eminus);
tcPDPtr gamma = getParticleData(ParticleID::gamma);
tcPDPtr Z0 = getParticleData(ParticleID::Z0);
// setup the processes
for ( int i =minflav_; i<=maxflav_; ++i ) {
tcPDPtr qk = getParticleData(i);
tcPDPtr qb = qk->CC();
add(new_ptr((Tree2toNDiagram(2), em, ep, 1, gamma, 3, qk, 3, qb, -1)));
add(new_ptr((Tree2toNDiagram(2), em, ep, 1, Z0 , 3, qk, 3, qb, -2)));
}
}
Energy2 MEee2gZ2qq::scale() const {
return sqr(getParticleData(ParticleID::Z0)->mass());
// return sHat();
}
unsigned int MEee2gZ2qq::orderInAlphaS() const {
return 0;
}
unsigned int MEee2gZ2qq::orderInAlphaEW() const {
return 2;
}
Selector<MEBase::DiagramIndex>
MEee2gZ2qq::diagrams(const DiagramVector & diags) const {
double lastCont(0.5),lastBW(0.5);
if ( lastXCombPtr() ) {
lastCont = meInfo()[0];
lastBW = meInfo()[1];
}
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if ( diags[i]->id() == -1 ) sel.insert(lastCont, i);
else if ( diags[i]->id() == -2 ) sel.insert(lastBW, i);
}
return sel;
}
Selector<const ColourLines *>
MEee2gZ2qq::colourGeometries(tcDiagPtr ) const {
static const ColourLines c("-5 4");
Selector<const ColourLines *> sel;
sel.insert(1.0, &c);
return sel;
}
void MEee2gZ2qq::persistentOutput(PersistentOStream & os) const {
os << FFZVertex_ << FFPVertex_ << FFGVertex_
<< Z0_ << gamma_ << gluon_ << minflav_
<< maxflav_ << massopt_ << alphaQCD_ << alphaQED_
- << ounit(pTmin_,GeV) << preFactor_;
+ << ounit(pTminQED_,GeV) << ounit(pTminQCD_,GeV) << preFactor_;
}
void MEee2gZ2qq::persistentInput(PersistentIStream & is, int) {
is >> FFZVertex_ >> FFPVertex_ >> FFGVertex_
>> Z0_ >> gamma_ >> gluon_ >> minflav_
>> maxflav_ >> massopt_ >> alphaQCD_ >> alphaQED_
- >> iunit(pTmin_,GeV) >> preFactor_;
+ >> iunit(pTminQED_,GeV) >> iunit(pTminQCD_,GeV) >> preFactor_;
}
// *** 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<MEee2gZ2qq,HwMEBase>
describeMEee2gZ2qq("Herwig::MEee2gZ2qq", "HwMELepton.so");
void MEee2gZ2qq::Init() {
static ClassDocumentation<MEee2gZ2qq> documentation
("The MEee2gZ2qq class implements the matrix element for e+e- -> q qbar");
static Parameter<MEee2gZ2qq,int> interfaceMinimumFlavour
("MinimumFlavour",
"The PDG code of the quark with the lowest PDG code to produce.",
&MEee2gZ2qq::minflav_, 1, 1, 6,
false, false, Interface::limited);
static Parameter<MEee2gZ2qq,int> interfaceMaximumFlavour
("MaximumFlavour",
"The PDG code of the quark with the highest PDG code to produce",
&MEee2gZ2qq::maxflav_, 5, 1, 6,
false, false, Interface::limited);
static Switch<MEee2gZ2qq,unsigned int> interfaceTopMassOption
("TopMassOption",
"Option for the treatment of the top quark mass",
&MEee2gZ2qq::massopt_, 1, false, false);
static SwitchOption interfaceTopMassOptionOnMassShell
(interfaceTopMassOption,
"OnMassShell",
"The top is produced on its mass shell",
1);
static SwitchOption interfaceTopMassOption2
(interfaceTopMassOption,
"OffShell",
"The top is generated off-shell using the mass and width generator.",
2);
- static Parameter<MEee2gZ2qq,Energy> interfacepTMin
- ("pTMin",
- "Minimum pT for hard radiation",
- &MEee2gZ2qq::pTmin_, GeV, 1.0*GeV, 0.001*GeV, 10.0*GeV,
+ static Parameter<MEee2gZ2qq,Energy> interfacepTMinQED
+ ("pTMinQED",
+ "Minimum pT for hard QED radiation",
+ &MEee2gZ2qq::pTminQED_, GeV, 1.0*GeV, 0.001*GeV, 10.0*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<MEee2gZ2qq,Energy> interfacepTMinQCD
+ ("pTMinQCD",
+ "Minimum pT for hard QCD radiation",
+ &MEee2gZ2qq::pTminQCD_, GeV, 1.0*GeV, 0.001*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<MEee2gZ2qq,double> interfacePrefactor
("Prefactor",
"Prefactor for the overestimate of the emission probability",
&MEee2gZ2qq::preFactor_, 6.0, 1.0, 100.0,
false, false, Interface::limited);
static Reference<MEee2gZ2qq,ShowerAlpha> interfaceQCDCoupling
("AlphaQCD",
"Pointer to the object to calculate the strong coupling for the correction",
&MEee2gZ2qq::alphaQCD_, false, false, true, false, false);
static Reference<MEee2gZ2qq,ShowerAlpha> interfaceEMCoupling
("AlphaQED",
"Pointer to the object to calculate the EM coupling for the correction",
&MEee2gZ2qq::alphaQED_, false, false, true, false, false);
}
double MEee2gZ2qq::me2() const {
return loME(mePartonData(),rescaledMomenta(),true);
}
ProductionMatrixElement MEee2gZ2qq::HelicityME(vector<SpinorWaveFunction> & fin,
vector<SpinorBarWaveFunction> & ain,
vector<SpinorBarWaveFunction> & fout,
vector<SpinorWaveFunction> & aout,
double & me,
double & cont,
double & BW ) const {
// the particles should be in the order
// for the incoming
// 0 incoming fermion (u spinor)
// 1 incoming antifermion (vbar spinor)
// for the outgoing
// 0 outgoing fermion (ubar spinor)
// 1 outgoing antifermion (v spinor)
// me to be returned
ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
ProductionMatrixElement gamma (PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
ProductionMatrixElement Zboson(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
// wavefunctions for the intermediate particles
VectorWaveFunction interZ,interG;
// temporary storage of the different diagrams
Complex diag1,diag2;
// sum over helicities to get the matrix element
unsigned int inhel1,inhel2,outhel1,outhel2;
double total[3]={0.,0.,0.};
for(inhel1=0;inhel1<2;++inhel1) {
for(inhel2=0;inhel2<2;++inhel2) {
// intermediate Z
interZ = FFZVertex_->evaluate(scale(),1,Z0_,fin[inhel1],ain[inhel2]);
// intermediate photon
interG = FFPVertex_->evaluate(scale(),1,gamma_,fin[inhel1],ain[inhel2]);
for(outhel1=0;outhel1<2;++outhel1) {
for(outhel2=0;outhel2<2;++outhel2) {
// first the Z exchange diagram
diag1 = FFZVertex_->evaluate(scale(),aout[outhel2],fout[outhel1],
interZ);
// then the photon exchange diagram
diag2 = FFPVertex_->evaluate(scale(),aout[outhel2],fout[outhel1],
interG);
// add up squares of individual terms
total[1] += norm(diag1);
Zboson(inhel1,inhel2,outhel1,outhel2) = diag1;
total[2] += norm(diag2);
gamma (inhel1,inhel2,outhel1,outhel2) = diag2;
// the full thing including interference
diag1 += diag2;
total[0] += norm(diag1);
output(inhel1,inhel2,outhel1,outhel2)=diag1;
}
}
}
}
for(int ix=0;ix<3;++ix) total[ix] *= 0.25;
tcPolarizedBeamPDPtr beam[2] =
{dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[0]),
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[1])};
if( beam[0] || beam[1] ) {
RhoDMatrix rho[2] =
{beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()),
beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())};
total[0] = output.average(rho[0],rho[1]);
total[1] = Zboson.average(rho[0],rho[1]);
total[2] = gamma .average(rho[0],rho[1]);
}
// results
for(int ix=0;ix<3;++ix) total[ix]*= 3.;
cont = total[2];
BW = total[1];
me = total[0];
return output;
}
void MEee2gZ2qq::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);
hard.push_back(sub->outgoing()[1]);
if(hard[0]->id()<hard[1]->id()) swap(hard[0],hard[1]);
if(hard[2]->id()<hard[3]->id()) swap(hard[2],hard[3]);
vector<SpinorWaveFunction> fin,aout;
vector<SpinorBarWaveFunction> ain,fout;
// get wave functions for off-shell momenta for later on
SpinorWaveFunction( fin ,hard[0],incoming,false,true);
SpinorBarWaveFunction(ain ,hard[1],incoming,false,true);
SpinorBarWaveFunction(fout,hard[2],outgoing,true ,true);
SpinorWaveFunction( aout,hard[3],outgoing,true ,true);
// now rescale the momenta and compute the matrix element with the
// rescaled momenta for correlations
vector<Lorentz5Momentum> momenta;
cPDVector data;
for(unsigned int ix=0;ix<4;++ix) {
momenta.push_back(hard[ix]->momentum());
data .push_back(hard[ix]->dataPtr());
}
rescaleMomenta(momenta,data);
SpinorWaveFunction ein (rescaledMomenta()[0],data[0],incoming);
SpinorBarWaveFunction pin (rescaledMomenta()[1],data[1],incoming);
SpinorBarWaveFunction qkout(rescaledMomenta()[2],data[2],outgoing);
SpinorWaveFunction qbout(rescaledMomenta()[3],data[3],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
ein.reset(ix) ; fin [ix] = ein ;
pin.reset(ix) ; ain [ix] = pin ;
qkout.reset(ix); fout[ix] = qkout;
qbout.reset(ix); aout[ix] = qbout;
}
// calculate the matrix element
double me,cont,BW;
ProductionMatrixElement prodme=HelicityME(fin,ain,fout,aout,me,cont,BW);
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(prodme);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix) {
tSpinPtr spin = hard[ix]->spinInfo();
if(ix<2) {
tcPolarizedBeamPDPtr beam =
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(hard[ix]->dataPtr());
if(beam) spin->rhoMatrix() = beam->rhoMatrix();
}
spin->productionVertex(hardvertex);
}
}
void MEee2gZ2qq::rebind(const TranslationMap & trans) {
FFZVertex_ = trans.translate(FFZVertex_);
FFPVertex_ = trans.translate(FFPVertex_);
FFGVertex_ = trans.translate(FFGVertex_);
Z0_ = trans.translate(Z0_);
gamma_ = trans.translate(gamma_);
gluon_ = trans.translate(gluon_);
HwMEBase::rebind(trans);
}
IVector MEee2gZ2qq::getReferences() {
IVector ret = HwMEBase::getReferences();
ret.push_back(FFZVertex_);
ret.push_back(FFPVertex_);
ret.push_back(FFGVertex_);
ret.push_back(Z0_ );
ret.push_back(gamma_ );
ret.push_back(gluon_ );
return ret;
}
void MEee2gZ2qq::initializeMECorrection(ShowerTreePtr , double & initial,
double & final) {
d_Q_ = sqrt(sHat());
d_m_ = 0.5*(meMomenta()[2].mass()+meMomenta()[3].mass());
// set the other parameters
d_rho_ = sqr(d_m_/d_Q_);
d_v_ = sqrt(1.-4.*d_rho_);
// maximum evolution scale
d_kt1_ = (1. + sqrt(1. - 4.*d_rho_))/2.;
double num = d_rho_ * d_kt1_ + 0.25 * d_v_ *(1.+d_v_)*(1.+d_v_);
double den = d_kt1_ - d_rho_;
d_kt2_ = num/den;
// maximums for reweighting
initial = 1.;
final = 1.;
}
void MEee2gZ2qq::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
vector<Lorentz5Momentum> emission;
unsigned int iemit,ispect;
- pair<Energy,ShowerInteraction::Type> output =
- generateHard(tree,emission,iemit,ispect,true,
- vector<ShowerInteraction::Type>(1,ShowerInteraction::QCD));
+ generateHard(tree,emission,iemit,ispect,true,
+ vector<ShowerInteraction::Type>(1,ShowerInteraction::QCD));
if(emission.empty()) return;
// get the quark and antiquark
ParticleVector qq;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit)
qq.push_back(cit->first->copy());
// ensure quark first
if(qq[0]->id()<0) swap(qq[0],qq[1]);
// perform final check to ensure energy greater than constituent mass
for (int i=0; i<2; i++) {
if (emission[i+2].e() < qq[i]->data().constituentMass()) return;
}
if (emission[4].e() < gluon_->constituentMass()) return;
// set masses
for (int i=0; i<2; i++) emission[i+2].setMass(qq[i]->mass());
emission[4].setMass(ZERO);
// decide which particle emits
bool firstEmits=
emission[4].vect().perp2(emission[2].vect())<
emission[4].vect().perp2(emission[3].vect());
// create the new quark, antiquark and gluon
PPtr newg = gluon_->produceParticle(emission[4]);
PPtr newq,newa;
if(firstEmits) {
newq = qq[0]->dataPtr()->produceParticle(emission[2]);
newa = new_ptr(Particle(*qq[1]));
qq[1]->antiColourLine()->removeAntiColoured(newa);
newa->set5Momentum(emission[3]);
}
else {
newq = new_ptr(Particle(*qq[0]));
qq[0]->colourLine()->removeColoured(newq);
newq->set5Momentum(emission[2]);
newa = qq[1]->dataPtr()->produceParticle(emission[3]);
}
// get the original colour line
ColinePtr col;
if(qq[0]->id()>0) col=qq[0]->colourLine();
else col=qq[0]->antiColourLine();
// set the colour lines
if(firstEmits) {
col->addColoured(newq);
col->addAntiColoured(newg);
newa->colourNeighbour(newg);
}
else {
col->addAntiColoured(newa);
col->addColoured(newg);
newq->antiColourNeighbour(newg);
}
// change the existing quark and antiquark
PPtr orig;
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(cit->first->progenitor()->id()==newq->id()) {
// remove old particles from colour line
col->removeColoured(cit->first->copy());
col->removeColoured(cit->first->progenitor());
// insert new particles
cit->first->copy(newq);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newq,1,true)));
cit->first->progenitor(sp);
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(!firstEmits);
if(firstEmits) orig=cit->first->original();
}
else {
// remove old particles from colour line
col->removeAntiColoured(cit->first->copy());
col->removeColoured(cit->first->progenitor());
// insert new particles
cit->first->copy(newa);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newa,1,true)));
cit->first->progenitor(sp);
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(firstEmits);
if(!firstEmits) orig=cit->first->original();
}
}
// 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);
}
bool MEee2gZ2qq::softMatrixElementVeto(ShowerProgenitorPtr initial,
ShowerParticlePtr parent,Branching br) {
// check we should be applying the veto
if(parent->id()!=initial->progenitor()->id()||
br.ids[0]!=br.ids[1]||
br.ids[2]!=ParticleID::g) return false;
// calculate pt
double d_z = br.kinematics->z();
Energy d_qt = br.kinematics->scale();
Energy2 d_m2 = parent->momentum().m2();
Energy pPerp = (1.-d_z)*sqrt( sqr(d_z*d_qt) - d_m2);
// if not hardest so far don't apply veto
if(pPerp<initial->highestpT()) return false;
// calculate x and xb
double kti = sqr(d_qt/d_Q_);
double w = sqr(d_v_) + kti*(-1. + d_z)*d_z*(2. + kti*(-1. + d_z)*d_z);
if (w < 0) {
initial->highestpT(pPerp);
return false;
}
double x = (1. + sqr(d_v_)*(-1. + d_z) + sqr(kti*(-1. + d_z))*d_z*d_z*d_z
+ d_z*sqrt(w)
- kti*(-1. + d_z)*d_z*(2. + d_z*(-2 + sqrt(w))))/
(1. - kti*(-1. + d_z)*d_z + sqrt(w));
double xb = 1. + kti*(-1. + d_z)*d_z;
// calculate the weight
if(parent->id()<0) swap(x,xb);
// if exceptionally out of phase space, leave this emission, as there
// is no good interpretation for the soft ME correction.
if( x<0 || xb<0) {
initial->highestpT(pPerp);
return false;
}
double xg = 2. - xb - x;
// always return one in the soft gluon region
if(xg < EPS_) {
initial->highestpT(pPerp);
return false;
}
// check it is in the phase space
if((1.-x)*(1.-xb)*(1.-xg) < d_rho_*xg*xg) {
parent->vetoEmission(parent->id()>0 ? ShowerPartnerType::QCDColourLine :
ShowerPartnerType::QCDColourLine,br.kinematics->scale());
return true;
}
double k1 = getKfromX(x, xb);
double k2 = getKfromX(xb, x);
double weight = 1.;
// quark emission
if(parent->id() > 0 && k1 < d_kt1_) {
weight = MEV(x, xb)/PS(x, xb);
// is it also in the anti-quark emission zone?
if(k2 < d_kt2_) weight *= 0.5;
}
// antiquark emission
if(parent->id() < 0 && k2 < d_kt2_) {
weight = MEV(x, xb)/PS(xb, x);
// is it also in the quark emission zone?
if(k1 < d_kt1_) weight *= 0.5;
}
// compute veto from weight
bool veto = !UseRandom::rndbool(weight);
// if not vetoed reset max
if(!veto) initial->highestpT(pPerp);
// if vetoing reset the scale
if(veto) {
parent->vetoEmission(parent->id()>0 ? ShowerPartnerType::QCDColourLine :
ShowerPartnerType::QCDColourLine,br.kinematics->scale());
}
// return the veto
return veto;
}
double MEee2gZ2qq::getKfromX(double x1, double x2) {
double uval = 0.5*(1. + d_rho_/(1.-x2+d_rho_));
double num = x1 - (2. - x2)*uval;
double den = sqrt(x2*x2 - 4.*d_rho_);
double zval = uval + num/den;
return (1.-x2)/(zval*(1.-zval));
}
double MEee2gZ2qq::MEV(double x1, double x2) {
// Vector part
double num = (x1+2.*d_rho_)*(x1+2.*d_rho_) + (x2+2.*d_rho_)*(x2+2.*d_rho_)
- 8.*d_rho_*(1.+2.*d_rho_);
double den = (1.+2.*d_rho_)*(1.-x1)*(1.-x2);
return (num/den - 2.*d_rho_/((1.-x1)*(1.-x1))
- 2*d_rho_/((1.-x2)*(1.-x2)))/d_v_;
}
double MEee2gZ2qq::PS(double x, double xbar) {
double u = 0.5*(1. + d_rho_ / (1.-xbar+d_rho_));
double z = u + (x - (2.-xbar)*u)/sqrt(xbar*xbar - 4.*d_rho_);
double brack = (1.+z*z)/(1.-z)- 2.*d_rho_/(1-xbar);
// interesting: the splitting function without the subtraction
// term. Actually gives a much worse approximation in the collinear
// limit. double brack = (1.+z*z)/(1.-z);
double den = (1.-xbar)*sqrt(xbar*xbar - 4.*d_rho_);
return brack/den;
}
pair<Energy,ShowerInteraction::Type>
MEee2gZ2qq::generateHard(ShowerTreePtr tree,
vector<Lorentz5Momentum> & emmision,
unsigned int & iemit, unsigned int & ispect,
bool applyVeto,
vector<ShowerInteraction::Type> inter) {
// get the momenta of the incoming and outgoing partons
// incoming
ShowerProgenitorPtr
emProgenitor = tree->incomingLines().begin() ->first,
epProgenitor = tree->incomingLines().rbegin()->first;
// outgoing
ShowerProgenitorPtr
qkProgenitor = tree->outgoingLines().begin() ->first,
qbProgenitor = tree->outgoingLines().rbegin()->first;
// get the order right
if(emProgenitor->id()<0) swap(emProgenitor,epProgenitor);
if(qkProgenitor->id()<0) swap(qkProgenitor,qbProgenitor);
loMomenta_.resize(0);
// extract the momenta
loMomenta_.push_back(emProgenitor->progenitor()->momentum());
loMomenta_.push_back(epProgenitor->progenitor()->momentum());
loMomenta_.push_back(qkProgenitor->progenitor()->momentum());
loMomenta_.push_back(qbProgenitor->progenitor()->momentum());
// and ParticleData objects
partons_.resize(5);
partons_[0]=emProgenitor->progenitor()->dataPtr();
partons_[1]=epProgenitor->progenitor()->dataPtr();
partons_[2]=qkProgenitor->progenitor()->dataPtr();
partons_[3]=qbProgenitor->progenitor()->dataPtr();
// boost from lab to CMS frame with outgoing particles
// along the z axis
LorentzRotation eventFrame( ( loMomenta_[2] + loMomenta_[3] ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*loMomenta_[2];
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() );
eventFrame.invert();
// mass of the final-state system
Energy2 M2 = (loMomenta_[2]+loMomenta_[3]).m2();
Energy M = sqrt(M2);
double mu1 = loMomenta_[2].mass()/M;
double mu2 = loMomenta_[3].mass()/M;
double mu12 = sqr(mu1), mu22 = sqr(mu2);
double lambda = sqrt(1.+sqr(mu12)+sqr(mu22)-2.*mu12-2.*mu22-2.*mu12*mu22);
// max pT
Energy pTmax = 0.5*sqrt(M2)*
(1.-sqr(loMomenta_[2].mass()+loMomenta_[3].mass())/M2);
// max y
- if ( pTmax < pTmin_ ) return make_pair(ZERO,ShowerInteraction::QCD);
- double ymax = acosh(pTmax/pTmin_);
+ if ( pTmax < pTminQED_ && pTmax < pTminQCD_ ) return make_pair(ZERO,ShowerInteraction::QCD);
vector<Energy> pTemit;
vector<vector<Lorentz5Momentum> > emittedMomenta;;
vector<unsigned int> iemitter,ispectater;
for(unsigned int iinter=0;iinter<inter.size();++iinter) {
- double a;
+ Energy pTmin(ZERO);
+ double a,ymax;
if(inter[iinter]==ShowerInteraction::QCD) {
+ pTmin = pTminQCD_;
+ ymax = acosh(pTmax/pTmin);
partons_[4] = gluon_;
// prefactor for the overestimate of the Sudakov
a = 4./3.*alphaQCD_->overestimateValue()/Constants::twopi*
2.*ymax*preFactor_;
}
else {
+ pTmin = pTminQED_;
+ ymax = acosh(pTmax/pTmin);
partons_[4] = gamma_;
a = alphaQED_->overestimateValue()/Constants::twopi*
2.*ymax*preFactor_*sqr(double(mePartonData()[2]->iCharge())/3.);
}
// variables for the emission
Energy pT[2];
double y[2],phi[2],x3[2],x1[2][2],x2[2][2];
double contrib[2][2];
// storage of the real emission momenta
vector<Lorentz5Momentum> realMomenta[2][2]=
{{vector<Lorentz5Momentum>(5),vector<Lorentz5Momentum>(5)},
{vector<Lorentz5Momentum>(5),vector<Lorentz5Momentum>(5)}};
for(unsigned int ix=0;ix<2;++ix)
for(unsigned int iy=0;iy<2;++iy)
for(unsigned int iz=0;iz<2;++iz)
realMomenta[ix][iy][iz] = loMomenta_[iz];
// generate the emission
for(unsigned int ix=0;ix<2;++ix) {
if(ix==1) {
swap(mu1 ,mu2 );
swap(mu12,mu22);
}
pT[ix] = pTmax;
y [ix] = 0.;
bool reject = true;
do {
// generate pT
pT[ix] *= pow(UseRandom::rnd(),1./a);
- if(pT[ix]<pTmin_) {
+ if(pT[ix]<pTmin) {
pT[ix] = -GeV;
break;
}
// generate y
y[ix] = -ymax+2.*UseRandom::rnd()*ymax;
// generate phi
phi[ix] = UseRandom::rnd()*Constants::twopi;
// calculate x3 and check in allowed region
x3[ix] = 2.*pT[ix]*cosh(y[ix])/M;
if(x3[ix] < 0. || x3[ix] > 1. -sqr( mu1 + mu2 ) ) continue;
// find the possible solutions for x1
double xT2 = sqr(2./M*pT[ix]);
double root = (-sqr(x3[ix])+xT2)*
(xT2*mu22+2.*x3[ix]-sqr(mu12)+2.*mu22+2.*mu12-sqr(x3[ix])-1.
+2.*mu12*mu22-sqr(mu22)-2.*mu22*x3[ix]-2.*mu12*x3[ix]);
double c1=2.*sqr(x3[ix])-4.*mu22-6.*x3[ix]+4.*mu12-xT2*x3[ix]
+2.*xT2-2.*mu12*x3[ix]+2.*mu22*x3[ix]+4.;
if(root<0.) continue;
x1[ix][0] = 1./(4.-4.*x3[ix]+xT2)*(c1-2.*sqrt(root));
x1[ix][1] = 1./(4.-4.*x3[ix]+xT2)*(c1+2.*sqrt(root));
// change sign of y if 2nd particle emits
if(ix==1) y[ix] *=-1.;
// loop over the solutions
for(unsigned int iy=0;iy<2;++iy) {
contrib[ix][iy]=0.;
// check x1 value allowed
if(x1[ix][iy]<2.*mu1||x1[ix][iy]>1.+mu12-mu22) continue;
// calculate x2 value and check allowed
x2[ix][iy] = 2.-x3[ix]-x1[ix][iy];
double root = max(0.,sqr(x1[ix][iy])-4.*mu12);
root = sqrt(root);
double x2min = 1.+mu22-mu12
-0.5*(1.-x1[ix][iy]+mu12-mu22)/(1.-x1[ix][iy]+mu12)*(x1[ix][iy]-2.*mu12+root);
double x2max = 1.+mu22-mu12
-0.5*(1.-x1[ix][iy]+mu12-mu22)/(1.-x1[ix][iy]+mu12)*(x1[ix][iy]-2.*mu12-root);
if(x2[ix][iy]<x2min||x2[ix][iy]>x2max) continue;
// check the z components
double z1 = sqrt(sqr(x1[ix][iy])-4.*mu12-xT2);
double z2 = -sqrt(sqr(x2[ix][iy])-4.*mu22);
double z3 = pT[ix]*sinh(y[ix])*2./M;
if(ix==1) z3 *=-1.;
if(abs(-z1+z2+z3)<1e-9) z1 *= -1.;
if(abs(z1+z2+z3)>1e-5) continue;
// if using as an ME correction the veto
if(applyVeto) {
double xb = x1[ix][iy], xc = x2[ix][iy];
double b = mu12, c = mu22;
double r = 0.5*(1.+b/(1.+c-xc));
double z1 = r + (xb-(2.-xc)*r)/sqrt(sqr(xc)-4.*c);
double kt1 = (1.-b+c-xc)/z1/(1.-z1);
r = 0.5*(1.+c/(1.+b-xb));
double z2 = r + (xc-(2.-xb)*r)/sqrt(sqr(xb)-4.*b);
double kt2 = (1.-c+b-xb)/z2/(1.-z2);
if(ix==1) {
swap(z1 ,z2);
swap(kt1,kt2);
}
// veto the shower region
if( kt1 < d_kt1_ || kt2 < d_kt2_ ) continue;
}
// construct the momenta
realMomenta[ix][iy][4] =
Lorentz5Momentum(pT[ix]*cos(phi[ix]),pT[ix]*sin(phi[ix]),
pT[ix]*sinh(y[ix]) ,pT[ix]*cosh(y[ix]),ZERO);
if(ix==0) {
realMomenta[ix][iy][2] =
Lorentz5Momentum(-pT[ix]*cos(phi[ix]),-pT[ix]*sin(phi[ix]),
z1*0.5*M,x1[ix][iy]*0.5*M,M*mu1);
realMomenta[ix][iy][3] =
Lorentz5Momentum(ZERO,ZERO, z2*0.5*M,x2[ix][iy]*0.5*M,M*mu2);
}
else {
realMomenta[ix][iy][2] =
Lorentz5Momentum(ZERO,ZERO,-z2*0.5*M,x2[ix][iy]*0.5*M,M*mu2);
realMomenta[ix][iy][3] =
Lorentz5Momentum(-pT[ix]*cos(phi[ix]),-pT[ix]*sin(phi[ix]),
-z1*0.5*M,x1[ix][iy]*0.5*M,M*mu1);
}
// boost the momenta back to the lab
for(unsigned int iz=2;iz<5;++iz)
realMomenta[ix][iy][iz] *= eventFrame;
// jacobian and prefactors for the weight
Energy J = M/sqrt(xT2)*abs(-x1[ix][iy]*x2[ix][iy]+2.*mu22*x1[ix][iy]
+x2[ix][iy]+x2[ix][iy]*mu12+mu22*x2[ix][iy]
-sqr(x2[ix][iy]))
/pow(sqr(x2[ix][iy])-4.*mu22,1.5);
// prefactors etc
contrib[ix][iy] = 0.5*pT[ix]/J/preFactor_/lambda;
// matrix element piece
contrib[ix][iy] *= meRatio(partons_,realMomenta[ix][iy],
ix,inter[iinter],false);
// coupling piece
if(inter[iinter]==ShowerInteraction::QCD)
contrib[ix][iy] *= alphaQCD_->ratio(sqr(pT[ix]));
else
contrib[ix][iy] *= alphaQED_->ratio(sqr(pT[ix]));
}
if(contrib[ix][0]+contrib[ix][1]>1.) {
ostringstream s;
s << "MEee2gZ2qq::generateHardest weight for channel " << ix
<< "is " << contrib[ix][0]+contrib[ix][1]
<< " which is greater than 1";
generator()->logWarning( Exception(s.str(), Exception::warning) );
}
reject = UseRandom::rnd() > contrib[ix][0] + contrib[ix][1];
}
while (reject);
- if(pT[ix]<pTmin_)
+ if(pT[ix]<pTmin)
pT[ix] = -GeV;
}
// pt of emission
if(pT[0]<ZERO && pT[1]<ZERO) {
pTemit.push_back(-GeV);
emittedMomenta.push_back(vector<Lorentz5Momentum>());
iemitter .push_back(0);
ispectater.push_back(0);
continue;
}
// now pick the emission with highest pT
vector<Lorentz5Momentum> emission;
if(pT[0]>pT[1]) {
iemitter .push_back(2);
ispectater.push_back(3);
pTemit.push_back(pT[0]);
if(UseRandom::rnd()<contrib[0][0]/(contrib[0][0]+contrib[0][1]))
emission = realMomenta[0][0];
else
emission = realMomenta[0][1];
}
else {
iemitter .push_back(3);
ispectater.push_back(2);
pTemit.push_back(pT[1]);
if(UseRandom::rnd()<contrib[1][0]/(contrib[1][0]+contrib[1][1]))
emission = realMomenta[1][0];
else
emission = realMomenta[1][1];
}
emittedMomenta.push_back(emission);
}
// select the type of emission
int iselect=-1;
pTmax = ZERO;
for(unsigned int ix=0;ix<inter.size();++ix) {
if(pTemit[ix]>pTmax) {
iselect = ix;
pTmax = pTemit[ix];
}
}
// no emission return
if(iselect<0) {
return make_pair(ZERO,ShowerInteraction::QCD);
}
partons_[4] = inter[iselect]==ShowerInteraction::QCD ? gluon_ : gamma_;
iemit = iemitter[iselect];
ispect = ispectater[iselect];
emmision = emittedMomenta[iselect];
// return pT of emission
return make_pair(pTmax,inter[iselect]);
}
HardTreePtr MEee2gZ2qq::generateHardest(ShowerTreePtr tree,
vector<ShowerInteraction::Type> inter) {
// generate the momenta for the hard emission
vector<Lorentz5Momentum> emmision;
unsigned int iemit,ispect;
pair<Energy,ShowerInteraction::Type> output
= generateHard(tree,emmision,iemit,ispect,false,inter);
Energy pTveto = output.first;
ShowerInteraction::Type force = output.second;
// incoming progenitors
ShowerProgenitorPtr
ePProgenitor = tree->incomingLines().begin() ->first,
eMProgenitor = tree->incomingLines().rbegin()->first;
if(eMProgenitor->id()<0) swap(eMProgenitor,ePProgenitor);
// outgoing progenitors
ShowerProgenitorPtr
qkProgenitor = tree->outgoingLines().begin() ->first,
qbProgenitor = tree->outgoingLines().rbegin()->first;
if(qkProgenitor->id()<0) swap(qkProgenitor,qbProgenitor);
// maximum pT of emission
if(emmision.empty()) {
for(unsigned int ix=0;ix<inter.size();++ix) {
- qkProgenitor->maximumpT(pTmin_,inter[ix]);
- qbProgenitor->maximumpT(pTmin_,inter[ix]);
+ if(inter[ix]==ShowerInteraction::QCD) {
+ qkProgenitor->maximumpT(pTminQCD_,inter[ix]);
+ qbProgenitor->maximumpT(pTminQCD_,inter[ix]);
+ }
+ else {
+ qkProgenitor->maximumpT(pTminQED_,inter[ix]);
+ qbProgenitor->maximumpT(pTminQED_,inter[ix]);
+ }
}
return HardTreePtr();
}
else {
for(unsigned int ix=0;ix<inter.size();++ix) {
qkProgenitor->maximumpT(pTveto,inter[ix]);
qbProgenitor->maximumpT(pTveto,inter[ix]);
}
}
// Make the particles for the hard tree
ShowerParticleVector hardParticles;
for(unsigned int ix=0;ix<partons_.size();++ix) {
hardParticles.push_back(new_ptr(ShowerParticle(partons_[ix],ix>=2)));
hardParticles.back()->set5Momentum(emmision[ix]);
}
ShowerParticlePtr parent(new_ptr(ShowerParticle(partons_[iemit],true)));
Lorentz5Momentum parentMomentum(emmision[iemit]+emmision[4]);
parentMomentum.setMass(partons_[iemit]->mass());
parent->set5Momentum(parentMomentum);
// Create the vectors of HardBranchings to create the HardTree:
vector<HardBranchingPtr> spaceBranchings,allBranchings;
// Incoming boson:
for(unsigned int ix=0;ix<2;++ix) {
spaceBranchings.push_back(new_ptr(HardBranching(hardParticles[ix],SudakovPtr(),
HardBranchingPtr(),
HardBranching::Incoming)));
allBranchings.push_back(spaceBranchings.back());
}
// Outgoing particles from hard emission:
HardBranchingPtr spectatorBranch(new_ptr(HardBranching(hardParticles[ispect],
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing)));
HardBranchingPtr emitterBranch(new_ptr(HardBranching(parent,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
if(force==ShowerInteraction::QED) {
emitterBranch->type(ShowerPartnerType::QED);
}
else {
emitterBranch->type(emitterBranch->branchingParticle()->id()>0 ?
ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine);
}
emitterBranch->addChild(new_ptr(HardBranching(hardParticles[iemit],
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing)));
emitterBranch->addChild(new_ptr(HardBranching(hardParticles[4],
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing)));
if(iemit==0) {
allBranchings.push_back(emitterBranch);
allBranchings.push_back(spectatorBranch);
}
else {
allBranchings.push_back( spectatorBranch );
allBranchings.push_back( emitterBranch );
}
emitterBranch ->branchingParticle()->partner(spectatorBranch->branchingParticle());
spectatorBranch->branchingParticle()->partner(emitterBranch ->branchingParticle());
if(force==ShowerInteraction::QED) {
spaceBranchings[0]->branchingParticle()->partner(spaceBranchings[1]->branchingParticle());
spaceBranchings[1]->branchingParticle()->partner(spaceBranchings[0]->branchingParticle());
}
// Make the HardTree from the HardBranching vectors.
HardTreePtr hardtree = new_ptr(HardTree(allBranchings,spaceBranchings,
force));
hardtree->partnersSet(true);
// Connect the particles with the branchings in the HardTree
hardtree->connect( eMProgenitor->progenitor(), allBranchings[0] );
tPPtr beam = eMProgenitor->original();
if(!beam->parents().empty()) beam = beam->parents()[0];
allBranchings[0]->beam(beam);
hardtree->connect( ePProgenitor->progenitor(), allBranchings[1] );
beam = ePProgenitor->original();
if(!beam->parents().empty()) beam = beam->parents()[0];
allBranchings[1]->beam(beam);
hardtree->connect( qkProgenitor->progenitor(), allBranchings[2] );
hardtree->connect( qbProgenitor->progenitor(), allBranchings[3] );
// colour flow
ColinePtr newline=new_ptr(ColourLine());
for(set<HardBranchingPtr>::const_iterator cit=hardtree->branchings().begin();
cit!=hardtree->branchings().end();++cit) {
if((**cit).branchingParticle()->dataPtr()->iColour()==PDT::Colour3)
newline->addColoured((**cit).branchingParticle());
else if((**cit).branchingParticle()->dataPtr()->iColour()==PDT::Colour3bar)
newline->addAntiColoured((**cit).branchingParticle());
}
// Return the HardTree
return hardtree;
}
double MEee2gZ2qq::meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter,
ShowerInteraction::Type inter,
bool subtract) const {
Lorentz5Momentum q = momenta[2]+momenta[3]+momenta[4];
Energy2 Q2=q.m2();
Energy2 lambda = sqrt((Q2-sqr(momenta[2].mass()+momenta[3].mass()))*
(Q2-sqr(momenta[2].mass()-momenta[3].mass())));
InvEnergy2 D[2];
double lome[2];
for(unsigned int iemit=0;iemit<2;++iemit) {
unsigned int ispect = iemit==0 ? 1 : 0;
Energy2 pipj = momenta[4 ] * momenta[2+iemit ];
Energy2 pipk = momenta[4 ] * momenta[2+ispect];
Energy2 pjpk = momenta[2+iemit] * momenta[2+ispect];
double y = pipj/(pipj+pipk+pjpk);
double z = pipk/( pipk+pjpk);
Energy mij = sqrt(2.*pipj+sqr(momenta[2+iemit].mass()));
Energy2 lamB = sqrt((Q2-sqr(mij+momenta[2+ispect].mass()))*
(Q2-sqr(mij-momenta[2+ispect].mass())));
Energy2 Qpk = q*momenta[2+ispect];
Lorentz5Momentum pkt =
lambda/lamB*(momenta[2+ispect]-Qpk/Q2*q)
+0.5/Q2*(Q2+sqr(momenta[2+ispect].mass())-sqr(momenta[2+ispect].mass()))*q;
Lorentz5Momentum pijt =
q-pkt;
double muj = momenta[2+iemit ].mass()/sqrt(Q2);
double muk = momenta[2+ispect].mass()/sqrt(Q2);
double vt = sqrt((1.-sqr(muj+muk))*(1.-sqr(muj-muk)))/(1.-sqr(muj)-sqr(muk));
double v = sqrt(sqr(2.*sqr(muk)+(1.-sqr(muj)-sqr(muk))*(1.-y))-4.*sqr(muk))
/(1.-y)/(1.-sqr(muj)-sqr(muk));
// dipole term
D[iemit] = 0.5/pipj*(2./(1.-(1.-z)*(1.-y))
-vt/v*(2.-z+sqr(momenta[2+iemit].mass())/pipj));
// matrix element
vector<Lorentz5Momentum> lomom(4);
lomom[0] = momenta[0];
lomom[1] = momenta[1];
if(iemit==0) {
lomom[2] = pijt;
lomom[3] = pkt ;
}
else {
lomom[3] = pijt;
lomom[2] = pkt ;
}
lome[iemit] = loME(partons,lomom,false)/3.;
}
InvEnergy2 ratio = realME(partons,momenta,inter)
*abs(D[iemitter])/(abs(D[0]*lome[0])+abs(D[1]*lome[1]));
double output = Q2*ratio;
if(subtract) output -= 2.*Q2*D[iemitter];
return output;
}
double MEee2gZ2qq::loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta,
bool first) const {
// compute the spinors
vector<SpinorWaveFunction> fin,aout;
vector<SpinorBarWaveFunction> ain,fout;
SpinorWaveFunction ein (momenta[0],partons[0],incoming);
SpinorBarWaveFunction pin (momenta[1],partons[1],incoming);
SpinorBarWaveFunction qkout(momenta[2],partons[2],outgoing);
SpinorWaveFunction qbout(momenta[3],partons[3],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
ein.reset(ix) ;
fin.push_back( ein );
pin.reset(ix) ;
ain.push_back( pin );
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
}
// compute the matrix element
double me,lastCont,lastBW;
HelicityME(fin,ain,fout,aout,me,lastCont,lastBW);
// save the components
if(first) {
DVector save;
save.push_back(lastCont);
save.push_back(lastBW);
meInfo(save);
}
// return the answer
return me;
}
InvEnergy2 MEee2gZ2qq::realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta,
ShowerInteraction::Type inter) const {
// compute the spinors
vector<SpinorWaveFunction> fin,aout;
vector<SpinorBarWaveFunction> ain,fout;
vector<VectorWaveFunction> gout;
SpinorWaveFunction ein (momenta[0],partons[0],incoming);
SpinorBarWaveFunction pin (momenta[1],partons[1],incoming);
SpinorBarWaveFunction qkout(momenta[2],partons[2],outgoing);
SpinorWaveFunction qbout(momenta[3],partons[3],outgoing);
VectorWaveFunction gluon(momenta[4],partons[4],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
ein.reset(ix) ;
fin.push_back( ein );
pin.reset(ix) ;
ain.push_back( pin );
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
gluon.reset(2*ix);
gout.push_back(gluon);
}
AbstractFFVVertexPtr vertex = inter == ShowerInteraction::QCD ?
FFGVertex_ : FFPVertex_;
vector<Complex> diag(4,0.);
ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1);
double total(0.);
for(unsigned int inhel1=0;inhel1<2;++inhel1) {
for(unsigned int inhel2=0;inhel2<2;++inhel2) {
// intermediate Z
VectorWaveFunction interZ =
FFZVertex_->evaluate(scale(),1,Z0_,fin[inhel1],ain[inhel2]);
// intermediate photon
VectorWaveFunction interG =
FFPVertex_->evaluate(scale(),1,gamma_,fin[inhel1],ain[inhel2]);
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
for(unsigned int outhel3=0;outhel3<2;++outhel3) {
SpinorBarWaveFunction off1 =
vertex->evaluate(scale(),3,partons[2],fout[outhel1],gout[outhel3]);
diag[0] = FFZVertex_->evaluate(scale(),aout[outhel2],off1,interZ);
diag[1] = FFPVertex_->evaluate(scale(),aout[outhel2],off1,interG);
SpinorWaveFunction off2 =
vertex->evaluate(scale(),3,partons[3],aout[outhel2],gout[outhel3]);
diag[2] = FFZVertex_->evaluate(scale(),off2,fout[outhel1],interZ);
diag[3] = FFPVertex_->evaluate(scale(),off2,fout[outhel1],interG);
// sum of diagrams
Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.));
// matrix element
output(inhel1,inhel2,outhel1,outhel2,outhel3)=sum;
// me2
total += norm(sum);
}
}
}
}
}
// spin average
total *= 0.25;
tcPolarizedBeamPDPtr beam[2] =
{dynamic_ptr_cast<tcPolarizedBeamPDPtr>(partons[0]),
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(partons[1])};
if( beam[0] || beam[1] ) {
RhoDMatrix rho[2] =
{beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()),
beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())};
total = output.average(rho[0],rho[1]);
}
// divide out the coupling
total /= norm(vertex->norm());
// and charge (if needed)
if(inter==ShowerInteraction::QED)
total /= sqr(double(mePartonData()[2]->iCharge())/3.);
// return the total
return total*UnitRemoval::InvE2;
}
diff --git a/MatrixElement/Lepton/MEee2gZ2qq.h b/MatrixElement/Lepton/MEee2gZ2qq.h
--- a/MatrixElement/Lepton/MEee2gZ2qq.h
+++ b/MatrixElement/Lepton/MEee2gZ2qq.h
@@ -1,491 +1,496 @@
// -*- C++ -*-
//
// MEee2gZ2qq.h 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.
//
#ifndef HERWIG_MEee2gZ2qq_H
#define HERWIG_MEee2gZ2qq_H
//
// This is the declaration of the MEee2gZ2qq class.
//
#include "Herwig++/MatrixElement/HwMEBase.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MEee2gZ2qq class implements the matrix element
* for \f$e^+e^-\to Z/\gamma \to q\bar{q}\f$ including spin correlations.
* The class includes greater control over the type of quark produced than is available
* in the corresponding matrix element from ThePEG, in addition to spin correlations.
*
* @see \ref MEee2gZ2qqInterfaces "The interfaces"
* defined for MEee2gZ2qq.
*/
class MEee2gZ2qq: public HwMEBase {
public:
/**
* The default constructor.
*/
- MEee2gZ2qq() : minflav_(1), maxflav_(5), massopt_(1), pTmin_(GeV),
+ MEee2gZ2qq() : minflav_(1), maxflav_(5), massopt_(1),
+ pTminQED_(GeV), pTminQCD_(GeV),
preFactor_(6.)
{}
/**
* Members for hard corrections to the emission of QCD radiation
*/
//@{
/**
* Has a POWHEG style correction
*/
virtual bool hasPOWHEGCorrection() {return true;}
/**
* Has an old fashioned ME correction
*/
virtual bool hasMECorrection() {return true;}
/**
* Initialize the ME correction
*/
virtual void initializeMECorrection(ShowerTreePtr, double &,
double & );
/**
* Apply the hard matrix element correction to a given hard process or decay
*/
virtual void applyHardMatrixElementCorrection(ShowerTreePtr);
/**
* Apply the soft matrix element correction
* @param initial The particle from the hard process which started the
* shower
* @param parent The initial particle in the current branching
* @param br The branching struct
* @return If true the emission should be vetoed
*/
virtual bool softMatrixElementVeto(ShowerProgenitorPtr initial,
ShowerParticlePtr parent,
Branching br);
/**
* Apply the POWHEG style correction
*/
virtual HardTreePtr generateHardest(ShowerTreePtr,
vector<ShowerInteraction::Type>);
//@}
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* The matrix element for the kinematical configuration
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
* @return the matrix element scaled with sHat() to give a
* dimensionless number.
*/
virtual double me2() const;
/**
* Return the scale associated with the last set phase space point.
*/
virtual Energy2 scale() const;
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Get diagram selector. With the information previously supplied with the
* setKinematics method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
* @param dv the diagrams to be weighted.
* @return a Selector relating the given diagrams to their weights.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
* @param diag the diagram chosen.
* @return the possible colour geometries weighted by their
* relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Construct the vertex of spin correlations.
*/
virtual void constructVertex(tSubProPtr);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans)
;
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
//@}
protected:
/**
* Calculate the matrix element for \f$e^-e^-\to q \bar q\f$.
* @param partons The incoming and outgoing particles
* @param momenta The momenta of the incoming and outgoing particles
* @param first Whether or not to calculate the spin correlations
*/
double loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta,
bool first) const;
/**
* Member to calculate the matrix element
* @param fin Spinors for incoming fermion
* @param ain Spinors for incoming antifermion
* @param fout Spinors for outgoing fermion
* @param aout Spinors for outgong antifermion
* @param me Spin summed Matrix element
* @param cont The continuum piece of the matrix element
* @param BW The Z piece of the matrix element
*/
ProductionMatrixElement HelicityME(vector<SpinorWaveFunction> & fin,
vector<SpinorBarWaveFunction> & ain,
vector<SpinorBarWaveFunction> & fout,
vector<SpinorWaveFunction> & aout,
double & me,
double & cont,
double & BW ) const;
/**
* The ratio of the matrix element for one additional jet over the
* leading order result. In practice
* \f[\frac{\hat{s}|\overline{\mathcal{M}}|^2_2|D_{\rm emit}|}{4\pi C_F\alpha_S|\overline{\mathcal{M}}|^2_3\left(|D_{\rm emit}|+|D_{\rm spect}|\right)}\f]
* is returned where \f$\|\overline{\mathcal{M}}|^2\f$ is
* the spin and colour summed/averaged matrix element.
* @param partons The incoming and outgoing particles
* @param momenta The momenta of the incoming and outgoing particles
* @param iemitter Whether the quark or antiquark is regardede as the emitter
* @param inter The type of interaction
* @param subtract Whether or not to subtract the relevant dipole term
*/
double meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter,
ShowerInteraction::Type inter,
bool subtract =false) const;
/**
* Calculate the matrix element for \f$e^-e^-\to q \bar q g\f$.
* @param partons The incoming and outgoing particles
* @param momenta The momenta of the incoming and outgoing particles
* @param inter The type of interaction
*/
InvEnergy2 realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta,
ShowerInteraction::Type inter) const;
private:
/**
* Generate the momenta for a hard configuration
*/
pair<Energy,ShowerInteraction::Type>
generateHard(ShowerTreePtr tree,
vector<Lorentz5Momentum> & emission,
unsigned int & iemit, unsigned int & ispect,
bool applyVeto,
vector<ShowerInteraction::Type>);
/**
* Calculate \f$\tilde{\kappa}\f$.
*/
double getKfromX(double, double);
/**
* Vector and axial vector parts of the matrix element
*/
//@{
/**
* Vector part of the matrix element
*/
double MEV(double, double);
/**
* The matrix element, given \f$x_1\f$, \f$x_2\f$.
* @param x1 \f$x_1\f$
* @param x2 \f$x_2\f$
*/
double PS(double x1, double x2);
//@}
protected:
/**
* Pointer to the fermion-antifermion Z vertex
*/
AbstractFFVVertexPtr FFZVertex() const {return FFZVertex_;}
/**
* Pointer to the fermion-antifermion photon vertex
*/
AbstractFFVVertexPtr FFPVertex() const {return FFPVertex_;}
/**
* Pointer to the particle data object for the Z
*/
PDPtr Z0() const {return Z0_;}
/**
* Pointer to the particle data object for the photon
*/
PDPtr gamma() const {return gamma_;}
/**
* Pointer to the particle data object for the gluon
*/
PDPtr gluon() const {return gluon_;}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEee2gZ2qq & operator=(const MEee2gZ2qq &);
private:
/**
* Parameters controlling the loead-order process
*/
//@{
/**
* The minimum PDG of the quarks to be produced
*/
int minflav_;
/**
* The maximum PDG of the quarks to be produced
*/
int maxflav_;
/**
* Option for the treatment of the top quark mass
*/
unsigned int massopt_;
//@}
/**
* Pointers to the vertices
*/
//@{
/**
* Pointer to the fermion-antifermion Z vertex
*/
AbstractFFVVertexPtr FFZVertex_;
/**
* Pointer to the fermion-antifermion photon vertex
*/
AbstractFFVVertexPtr FFPVertex_;
/**
* Pointer to the fermion-antifermion photon vertex
*/
AbstractFFVVertexPtr FFGVertex_;
//@}
/**
* Pointer to the ParticleData objects
*/
//@{
/**
* Pointer to the particle data object for the Z
*/
PDPtr Z0_;
/**
* Pointer to the particle data object for the photon
*/
PDPtr gamma_;
/**
* Pointer to the particle data object for the gluon
*/
PDPtr gluon_;
//@}
/**
* CM energy
*/
Energy d_Q_;
/**
* Quark mass
*/
Energy d_m_;
/**
* The rho parameter
*/
double d_rho_;
/**
* The v parameter
*/
double d_v_;
/**
* The initial kappa-tilde values for radiation from the quark
*/
double d_kt1_;
/**
* The initial kappa-tilde values for radiation from the antiquark
*/
double d_kt2_;
/**
* Cut-off parameter
*/
static const double EPS_;
/**
* Pointer to the strong coupling
*/
ShowerAlphaPtr alphaQCD_;
/**
* Pointer to the EM coupling
*/
ShowerAlphaPtr alphaQED_;
private:
/**
* Variables for the POWHEG style corrections
*/
//@{
/**
- * The cut off on pt, assuming massless quarks.
+ * The cut off on pt for QED, assuming massless quarks.
*/
- Energy pTmin_;
+ Energy pTminQED_;
+ /**
+ * The cut off on pt for QCD, assuming massless quarks.
+ */
+ Energy pTminQCD_;
/**
* Overestimate for the prefactor
*/
double preFactor_;
/**
* ParticleData objects for the partons
*/
vector<cPDPtr> partons_;
/**
* Momenta of the leading-order partons
*/
vector<Lorentz5Momentum> loMomenta_;
//@}
};
}
#endif /* HERWIG_MEee2gZ2qq_H */
diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
--- a/Shower/Base/Evolver.cc
+++ b/Shower/Base/Evolver.cc
@@ -1,2123 +1,2126 @@
// -*- C++ -*-
//
// Evolver.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 Evolver class.
//
#include "Evolver.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ShowerKinematics.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Utilities/Throw.h"
#include "ShowerTree.h"
#include "ShowerProgenitor.h"
#include "KinematicsReconstructor.h"
#include "PartnerFinder.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig++/Shower/ShowerHandler.h"
using namespace Herwig;
namespace {
void findChildren(tShowerParticlePtr parent,set<ShowerParticlePtr> & fs) {
for(unsigned int ix=0;ix<parent->children().size();++ix) {
tShowerParticlePtr child=
dynamic_ptr_cast<tShowerParticlePtr>(parent->children()[ix]);
if(child) findChildren(child,fs);
}
if(parent->children().empty()) {
if(parent->isFinalState()) fs.insert(parent);
}
}
}
IBPtr Evolver::clone() const {
return new_ptr(*this);
}
IBPtr Evolver::fullclone() const {
return new_ptr(*this);
}
void Evolver::persistentOutput(PersistentOStream & os) const {
os << _model << _splittingGenerator << _maxtry
<< _meCorrMode << _hardVetoMode << _hardVetoRead << _hardVetoReadOption
<< _limitEmissions
<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
<< _vetoes << _trunc_Mode << _hardEmissionMode
<< _colourEvolutionMethod << _reconOpt
<< interaction_<< interactions_.size();
for(unsigned int ix=0;ix<interactions_.size();++ix)
os << oenum(interactions_[ix]);
}
void Evolver::persistentInput(PersistentIStream & is, int) {
unsigned int isize;
is >> _model >> _splittingGenerator >> _maxtry
>> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _hardVetoReadOption
>> _limitEmissions
>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
>> _vetoes >> _trunc_Mode >> _hardEmissionMode
>> _colourEvolutionMethod >> _reconOpt
>> interaction_ >> isize;
interactions_.resize(isize);
for(unsigned int ix=0;ix<interactions_.size();++ix)
is >> ienum(interactions_[ix]);
}
void Evolver::doinit() {
Interfaced::doinit();
if(interaction_==0) {
interactions_.push_back(ShowerInteraction::QCD);
interactions_.push_back(ShowerInteraction::QED);
}
else if(interaction_==1) {
interactions_.push_back(ShowerInteraction::QCD);
}
else if(interaction_==2) {
interactions_.push_back(ShowerInteraction::QED);
interactions_.push_back(ShowerInteraction::QCD);
}
else if(interaction_==3) {
interactions_.push_back(ShowerInteraction::QED);
}
else if(interaction_==4) {
interactions_.push_back(ShowerInteraction::Both);
}
}
ClassDescription<Evolver> Evolver::initEvolver;
// Definition of the static class description member.
void Evolver::Init() {
static ClassDocumentation<Evolver> documentation
("This class is responsible for carrying out the showering,",
"including the kinematics reconstruction, in a given scale range,"
"including the option of the POWHEG approach to simulated next-to-leading order"
" radiation\\cite{Nason:2004rx}.",
"%\\cite{Nason:2004rx}\n"
"\\bibitem{Nason:2004rx}\n"
" P.~Nason,\n"
" ``A new method for combining NLO QCD with shower Monte Carlo algorithms,''\n"
" JHEP {\\bf 0411} (2004) 040\n"
" [arXiv:hep-ph/0409146].\n"
" %%CITATION = JHEPA,0411,040;%%\n");
static Reference<Evolver,SplittingGenerator>
interfaceSplitGen("SplittingGenerator",
"A reference to the SplittingGenerator object",
&Herwig::Evolver::_splittingGenerator,
false, false, true, false);
static Reference<Evolver,ShowerModel> interfaceShowerModel
("ShowerModel",
"The pointer to the object which defines the shower evolution model.",
&Evolver::_model, false, false, true, false, false);
static Parameter<Evolver,unsigned int> interfaceMaxTry
("MaxTry",
"The maximum number of attempts to generate the shower from a"
" particular ShowerTree",
&Evolver::_maxtry, 100, 1, 1000,
false, false, Interface::limited);
static Switch<Evolver, unsigned int> ifaceMECorrMode
("MECorrMode",
"Choice of the ME Correction Mode",
&Evolver::_meCorrMode, 1, false, false);
static SwitchOption off
(ifaceMECorrMode,"No","MECorrections off", 0);
static SwitchOption on
(ifaceMECorrMode,"Yes","hard+soft on", 1);
static SwitchOption hard
(ifaceMECorrMode,"Hard","only hard on", 2);
static SwitchOption soft
(ifaceMECorrMode,"Soft","only soft on", 3);
static Switch<Evolver, unsigned int> ifaceHardVetoMode
("HardVetoMode",
"Choice of the Hard Veto Mode",
&Evolver::_hardVetoMode, 1, false, false);
static SwitchOption HVoff
(ifaceHardVetoMode,"No","hard vetos off", 0);
static SwitchOption HVon
(ifaceHardVetoMode,"Yes","hard vetos on", 1);
static SwitchOption HVIS
(ifaceHardVetoMode,"Initial", "only IS emissions vetoed", 2);
static SwitchOption HVFS
(ifaceHardVetoMode,"Final","only FS emissions vetoed", 3);
static Switch<Evolver, unsigned int> ifaceHardVetoRead
("HardVetoScaleSource",
"If hard veto scale is to be read",
&Evolver::_hardVetoRead, 0, false, false);
static SwitchOption HVRcalc
(ifaceHardVetoRead,"Calculate","Calculate from hard process", 0);
static SwitchOption HVRread
(ifaceHardVetoRead,"Read","Read from XComb->lastScale", 1);
static Switch<Evolver, bool> ifaceHardVetoReadOption
("HardVetoReadOption",
"Apply read-in scale veto to all collisions or just the primary one?",
&Evolver::_hardVetoReadOption, false, false, false);
static SwitchOption AllCollisions
(ifaceHardVetoReadOption,
"AllCollisions",
"Read-in pT veto applied to primary and secondary collisions.",
false);
static SwitchOption PrimaryCollision
(ifaceHardVetoReadOption,
"PrimaryCollision",
"Read-in pT veto applied to primary but not secondary collisions.",
true);
static Parameter<Evolver, Energy> ifaceiptrms
("IntrinsicPtGaussian",
"RMS of intrinsic pT of Gaussian distribution:\n"
"2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
&Evolver::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
false, false, Interface::limited);
static Parameter<Evolver, double> ifacebeta
("IntrinsicPtBeta",
"Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
"(1-Beta) is the proportion of Gaussian distribution",
&Evolver::_beta, 0, 0, 1,
false, false, Interface::limited);
static Parameter<Evolver, Energy> ifacegamma
("IntrinsicPtGamma",
"Parameter for inverse quadratic:\n"
"2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))",
&Evolver::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static Parameter<Evolver, Energy> ifaceiptmax
("IntrinsicPtIptmax",
"Upper bound on intrinsic pT for inverse quadratic",
&Evolver::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static RefVector<Evolver,ShowerVeto> ifaceVetoes
("Vetoes",
"The vetoes to be checked during showering",
&Evolver::_vetoes, -1,
false,false,true,true,false);
static Switch<Evolver,unsigned int> interfaceLimitEmissions
("LimitEmissions",
"Limit the number and type of emissions for testing",
&Evolver::_limitEmissions, 0, false, false);
static SwitchOption interfaceLimitEmissionsNoLimit
(interfaceLimitEmissions,
"NoLimit",
"Allow an arbitrary number of emissions",
0);
static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
(interfaceLimitEmissions,
"OneInitialStateEmission",
"Allow one emission in the initial state and none in the final state",
1);
static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
(interfaceLimitEmissions,
"OneFinalStateEmission",
"Allow one emission in the final state and none in the initial state",
2);
static SwitchOption interfaceLimitEmissionsHardOnly
(interfaceLimitEmissions,
"HardOnly",
"Only allow radiation from the hard ME correction",
3);
static SwitchOption interfaceLimitEmissionsOneEmission
(interfaceLimitEmissions,
"OneEmission",
"Allow one emission in either the final state or initial state, but not both",
4);
static Switch<Evolver,bool> interfaceTruncMode
("TruncatedShower", "Include the truncated shower?",
&Evolver::_trunc_Mode, 1, false, false);
static SwitchOption interfaceTruncMode0
(interfaceTruncMode,"No","Truncated Shower is OFF", 0);
static SwitchOption interfaceTruncMode1
(interfaceTruncMode,"Yes","Truncated Shower is ON", 1);
static Switch<Evolver,unsigned int> interfaceHardEmissionMode
("HardEmissionMode",
"Whether to use ME corrections or POWHEG for the hardest emission",
&Evolver::_hardEmissionMode, 0, false, false);
static SwitchOption interfaceHardEmissionModeMECorrection
(interfaceHardEmissionMode,
"MECorrection",
"Old fashioned ME correction",
0);
static SwitchOption interfaceHardEmissionModePOWHEG
(interfaceHardEmissionMode,
"POWHEG",
"Powheg style hard emission",
1);
static Switch<Evolver,int> interfaceColourEvolutionMethod
("ColourEvolutionMethod",
"Choice of method for choosing the colour factor in gluon evolution",
&Evolver::_colourEvolutionMethod, 0, false, false);
static SwitchOption interfaceColourEvolutionMethodDefault
(interfaceColourEvolutionMethod,
"Default",
"Colour factor is CA for all scales",
0);
static SwitchOption interfaceColourEvolutionMethodHalfCA
(interfaceColourEvolutionMethod,
"HalfCA",
"Only use half the normal radiation until second scale is reached",
1);
static Switch<Evolver,unsigned int > interfaceInteractions
("Interactions",
"The interactions to be used in the shower",
&Evolver::interaction_, 1, false, false);
static SwitchOption interfaceInteractionsQCDFirst
(interfaceInteractions,
"QCDFirst",
"QCD first then QED",
0);
static SwitchOption interfaceInteractionsQCDOnly
(interfaceInteractions,
"QCDOnly",
"Only QCD",
1);
static SwitchOption interfaceInteractionsQEDFirst
(interfaceInteractions,
"QEDFirst",
"QED first then QCD",
2);
static SwitchOption interfaceInteractionsQEDOnly
(interfaceInteractions,
"QEDOnly",
"Only QED",
3);
static SwitchOption interfaceInteractionsBothAtOnce
(interfaceInteractions,
"BothAtOnce",
"Generate both at the same time",
4);
static Switch<Evolver,unsigned int> interfaceReconstructionOption
("ReconstructionOption",
"Treatment of the reconstruction of the transverse momentum of "
"a branching from the evolution scale.",
&Evolver::_reconOpt, 0, false, false);
static SwitchOption interfaceReconstructionOptionCutOff
(interfaceReconstructionOption,
"CutOff",
"Use the cut-off masses in the calculation",
0);
static SwitchOption interfaceReconstructionOptionOffShell
(interfaceReconstructionOption,
"OffShell",
"Use the off-shell masses in the calculation",
1);
}
void Evolver::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
_intrinsic.clear();
if ( !ipTon() || !isISRadiationON() ) return;
// don't do anything for the moment for secondary scatters
if( !ShowerHandler::currentHandler()->firstInteraction() ) return;
// generate intrinsic pT
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
// only consider initial-state particles
if(particlesToShower[ix]->progenitor()->isFinalState()) continue;
if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue;
Energy ipt;
if(UseRandom::rnd() > _beta) {
ipt=_iptrms*sqrt(-log(UseRandom::rnd()));
}
else {
ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.);
}
pair<Energy,double> pt = make_pair(ipt,UseRandom::rnd(Constants::twopi));
_intrinsic[particlesToShower[ix]] = pt;
}
}
void Evolver::setupMaximumScales(vector<ShowerProgenitorPtr> & p) {
// return if no vetos
if (!hardVetoOn()) return;
// find out if hard partonic subprocess.
bool isPartonic(false);
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = _currenttree->incomingLines().begin();
Lorentz5Momentum pcm;
for(; cit!=currentTree()->incomingLines().end(); ++cit) {
pcm += cit->first->progenitor()->momentum();
isPartonic |= cit->first->progenitor()->coloured();
}
// find maximum pt from hard process, the maximum pt from all outgoing
// coloured lines (this is simpler and more general than
// 2stu/(s^2+t^2+u^2)). Maximum scale for scattering processes will
// be transverse mass.
Energy ptmax = -1.0*GeV;
// general case calculate the scale
if (!hardVetoXComb()||
(hardVetoReadOption()&&
!ShowerHandler::currentHandler()->firstInteraction())) {
// scattering process
if(currentTree()->isHard()) {
// coloured incoming particles
if (isPartonic) {
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cjt = currentTree()->outgoingLines().begin();
for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) {
if (cjt->first->progenitor()->coloured())
ptmax = max(ptmax,cjt->first->progenitor()->momentum().mt());
}
}
if (ptmax < ZERO) ptmax = pcm.m();
if(hardVetoXComb()&&hardVetoReadOption()&&
!ShowerHandler::currentHandler()->firstInteraction()) {
ptmax=min(ptmax,sqrt(ShowerHandler::currentHandler()
->lastXCombPtr()->lastScale()));
}
}
// decay, incoming() is the decaying particle.
else {
ptmax = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
}
// hepeup.SCALUP is written into the lastXComb by the
// LesHouchesReader itself - use this by user's choice.
// Can be more general than this.
else {
ptmax = sqrt( ShowerHandler::currentHandler()
->lastXCombPtr()->lastScale() );
}
// set maxHardPt for all progenitors. For partonic processes this
// is now the max pt in the FS, for non-partonic processes or
// processes with no coloured FS the invariant mass of the IS
vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax);
}
void Evolver::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
_hardme = HwMEBasePtr();
// extract the matrix element
tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
if(lastXC) {
_hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
}
_decayme = HwDecayerBasePtr();
// set the current tree
currentTree(hard);
hardTree(HardTreePtr());
// number of attempts if more than one interaction switched on
unsigned int interactionTry=0;
do {
try {
// generate the showering
doShowering(true);
// if no vetos return
return;
}
catch (InteractionVeto) {
currentTree()->clear();
++interactionTry;
}
}
while(interactionTry<=5);
throw Exception() << "Too many tries for shower in "
<< "Evolver::showerHardProcess()"
<< Exception::eventerror;
}
void Evolver::hardMatrixElementCorrection(bool hard) {
// set the initial enhancement factors for the soft correction
_initialenhance = 1.;
_finalenhance = 1.;
// if hard matrix element switched off return
if(!MECOn()) return;
// see if we can get the correction from the matrix element
// or decayer
if(hard) {
if(_hardme&&_hardme->hasMECorrection()) {
_hardme->initializeMECorrection(_currenttree,
_initialenhance,_finalenhance);
if(hardMEC())
_hardme->applyHardMatrixElementCorrection(_currenttree);
}
}
else {
if(_decayme&&_decayme->hasMECorrection()) {
_decayme->initializeMECorrection(_currenttree,
_initialenhance,_finalenhance);
if(hardMEC())
_decayme->applyHardMatrixElementCorrection(_currenttree);
}
}
}
bool Evolver::timeLikeShower(tShowerParticlePtr particle,
ShowerInteraction::Type type,
bool first) {
// don't do anything if not needed
if(_limitEmissions == 1 || hardOnly() ||
( _limitEmissions == 2 && _nfs != 0) ||
( _limitEmissions == 4 && _nfs + _nis != 0) ) return false;
ShowerParticleVector theChildren;
int ntry=0;
do {
++ntry;
// generate the emission
Branching fb;
while (true) {
fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
// no emission return
if(!fb.kinematics) return false;
// if emission OK break
if(!timeLikeVetoed(fb,particle)) break;
// otherwise reset scale and continue - SO IS involved in veto algorithm
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
tcPDPtr pdata[2];
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
// update the children
particle->showerKinematics()->
updateChildren(particle, theChildren,fb.type);
// update number of emissions
++_nfs;
if(_limitEmissions!=0) return true;
// shower the first particle
timeLikeShower(theChildren[0],type,false);
// shower the second particle
timeLikeShower(theChildren[1],type,false);
// that's if for old approach
if(_reconOpt==0) break;
// branching has happened
particle->showerKinematics()->
updateParent(particle, theChildren,fb.type);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
}
}
while(particle->virtualMass()==ZERO&&ntry<50);
if(first)
particle->showerKinematics()->resetChildren(particle,theChildren);
return true;
}
bool
Evolver::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
ShowerInteraction::Type type) {
//using the pdf's associated with the ShowerHandler assures, that
//modified pdf's are used for the secondary interactions via
//CascadeHandler::resetPDFs(...)
tcPDFPtr pdf;
if(ShowerHandler::currentHandler()->firstPDF().particle() == _beam)
pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
if(ShowerHandler::currentHandler()->secondPDF().particle() == _beam)
pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
// don't do anything if not needed
if(_limitEmissions == 2 || hardOnly() ||
( _limitEmissions == 1 && _nis != 0 ) ||
( _limitEmissions == 4 && _nis + _nfs != 0 ) ) return false;
Branching bb;
// generate branching
while (true) {
bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
_initialenhance,
_beam,type,
pdf,freeze);
// return if no emission
if(!bb.kinematics) return false;
// if not vetoed break
if(!spaceLikeVetoed(bb,particle)) break;
// otherwise reset scale and continue
particle->vetoEmission(bb.type,bb.kinematics->scale());
}
// assign the splitting function and shower kinematics
particle->showerKinematics(bb.kinematics);
// For the time being we are considering only 1->2 branching
// particles as in Sudakov form factor
tcPDPtr part[2]={getParticleData(bb.ids[0]),
getParticleData(bb.ids[2])};
if(particle->id()!=bb.ids[1]) {
if(part[0]->CC()) part[0]=part[0]->CC();
if(part[1]->CC()) part[1]=part[1]->CC();
}
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent=new_ptr(ShowerParticle(part[0],false));
ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true));
ShowerParticleVector theChildren;
theChildren.push_back(particle);
theChildren.push_back(otherChild);
//this updates the evolution scale
particle->showerKinematics()->
updateParent(newParent, theChildren,bb.type);
// update the history if needed
_currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
_currenttree->addInitialStateBranching(particle,newParent,otherChild);
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
++_nis;
bool emitted = _limitEmissions==0 ?
spaceLikeShower(newParent,beam,type) : false;
// now reconstruct the momentum
if(!emitted) {
if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
bb.kinematics->updateLast(newParent,ZERO,ZERO);
}
else {
pair<Energy,double> kt=_intrinsic[_progenitor];
bb.kinematics->updateLast(newParent,
kt.first*cos(kt.second),
kt.first*sin(kt.second));
}
}
particle->showerKinematics()->
updateChildren(newParent, theChildren,bb.type);
if(_limitEmissions!=0) return true;
// perform the shower of the final-state particle
timeLikeShower(otherChild,type,true);
// return the emitted
return true;
}
void Evolver::showerDecay(ShowerTreePtr decay) {
_decayme = HwDecayerBasePtr();
_hardme = HwMEBasePtr();
// find the decayer
// try the normal way if possible
tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode();
if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode();
if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
// otherwise make a string and look it up
if(!dm) {
string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name()
+ "->";
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
if(it!=decay->outgoingLines().begin()) tag += ",";
tag += it->first->original()->dataPtr()->name();
}
tag += ";";
dm = generator()->findDecayMode(tag);
}
if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
// set the ShowerTree to be showered
currentTree(decay);
decay->applyTransforms();
hardTree(HardTreePtr());
unsigned int interactionTry=0;
do {
try {
// generate the showering
doShowering(false);
// if no vetos return
return;
}
catch (InteractionVeto) {
currentTree()->clear();
++interactionTry;
}
}
while(interactionTry<=5);
throw Exception() << "Too many tries for QED shower in Evolver::showerDecay()"
<< Exception::eventerror;
}
bool Evolver::spaceLikeDecayShower(tShowerParticlePtr particle,
const ShowerParticle::EvolutionScales & maxScales,
Energy minmass,ShowerInteraction::Type type) {
Branching fb;
while (true) {
fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
_initialenhance,type);
// return if no radiation
if(!fb.kinematics) return false;
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle)) break;
// otherwise reset scale and continue
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
tcPDPtr pdata[2];
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
// some code moved to updateChildren
particle->showerKinematics()->
updateChildren(particle, theChildren, fb.type);
// In the case of splittings which involves coloured particles,
// set properly the colour flow of the branching.
// update the history if needed
_currenttree->updateInitialStateShowerProduct(_progenitor,theChildren[0]);
_currenttree->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the first particle
spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
// shower the second particle
timeLikeShower(theChildren[1],type,true);
// branching has happened
return true;
}
vector<ShowerProgenitorPtr> Evolver::setupShower(bool hard) {
// generate POWHEG hard emission if needed
if(_hardEmissionMode==1) hardestEmission(hard);
+ ShowerInteraction::Type inter = interactions_[0];
+ if(_hardtree&&inter!=ShowerInteraction::Both) {
+ inter = _hardtree->interaction();
+ }
// set the initial colour partners
- setEvolutionPartners(hard,_hardtree ?
- _hardtree->interaction() : interactions_[0],false);
+ setEvolutionPartners(hard,inter,false);
// generate hard me if needed
if(_hardEmissionMode==0) hardMatrixElementCorrection(hard);
// get the particles to be showered
vector<ShowerProgenitorPtr> particlesToShower =
currentTree()->extractProgenitors();
// remake the colour partners if needed
if(_hardEmissionMode==0 && _currenttree->hardMatrixElementCorrection()) {
setEvolutionPartners(hard,interactions_[0],true);
_currenttree->resetShowerProducts();
}
// return the answer
return particlesToShower;
}
void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type type,
bool clear) {
// match the particles in the ShowerTree and hardTree
if(hardTree() && !hardTree()->connect(currentTree()))
throw Exception() << "Can't match trees in "
<< "Evolver::setEvolutionPartners()"
<< Exception::eventerror;
// extract the progenitors
vector<ShowerParticlePtr> particles =
currentTree()->extractProgenitorParticles();
// sort out the colour partners
if(hardTree()) {
// find the partner
for(unsigned int ix=0;ix<particles.size();++ix) {
tHardBranchingPtr partner =
hardTree()->particles()[particles[ix]]->colourPartner();
if(!partner) continue;
for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
it=hardTree()->particles().begin();
it!=hardTree()->particles().end();++it) {
if(it->second==partner) particles[ix]->partner(it->first);
}
if(!particles[ix]->partner())
throw Exception() << "Can't match partners in "
<< "Evolver::setEvolutionPartners()"
<< Exception::eventerror;
}
}
// Set the initial evolution scales
if(clear) {
for(unsigned int ix=0;ix<particles.size();++ix) {
particles[ix]->partner(ShowerParticlePtr());
particles[ix]->clearPartners();
}
}
showerModel()->partnerFinder()->
setInitialEvolutionScales(particles,!hard,type,!_hardtree);
}
void Evolver::updateHistory(tShowerParticlePtr particle) {
if(!particle->children().empty()) {
ShowerParticleVector theChildren;
for(unsigned int ix=0;ix<particle->children().size();++ix) {
ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
(particle->children()[ix]);
theChildren.push_back(part);
}
// update the history if needed
if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
_currenttree->updateFinalStateShowerProduct(_progenitor,
particle,theChildren);
_currenttree->addFinalStateBranching(particle,theChildren);
for(unsigned int ix=0;ix<theChildren.size();++ix)
updateHistory(theChildren[ix]);
}
}
bool Evolver::startTimeLikeShower(ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && !mit->second->children().empty() ) {
bool output=truncatedTimeLikeShower(progenitor()->progenitor(),
mit->second ,type);
if(output) updateHistory(progenitor()->progenitor());
return output;
}
}
bool output = hardOnly() ? false :
timeLikeShower(progenitor()->progenitor() ,type,true) ;
if(output) updateHistory(progenitor()->progenitor());
return output;
}
bool Evolver::startSpaceLikeShower(PPtr parent, ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit =hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && mit->second->parent() ) {
return truncatedSpaceLikeShower( progenitor()->progenitor(),
parent, mit->second->parent(), type );
}
}
return hardOnly() ? false :
spaceLikeShower(progenitor()->progenitor(),parent,type);
}
bool Evolver::
startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
Energy minimumMass,ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit =hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && mit->second->parent() ) {
HardBranchingPtr branch=mit->second;
while(branch->parent()) branch=branch->parent();
return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
minimumMass, branch ,type);
}
}
return hardOnly() ? false :
spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type);
}
bool Evolver::timeLikeVetoed(const Branching & fb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction::Type type = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// check whether emission was harder than largest pt of hard subprocess
if ( hardVetoFS() && fb.kinematics->pT() > _progenitor->maxHardPt() )
return true;
// soft matrix element correction veto
if( softMEC()) {
if(_hardme && _hardme->hasMECorrection()) {
if(_hardme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
else if(_decayme && _decayme->hasMECorrection()) {
if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
}
// veto on maximum pt
if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
// general vetos
if (fb.kinematics && !_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoTimeLike(_progenitor,particle,fb);
switch((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
}
if(vetoed) return true;
}
return false;
}
bool Evolver::spaceLikeVetoed(const Branching & bb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction::Type type = bb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// check whether emission was harder than largest pt of hard subprocess
if (hardVetoIS() && bb.kinematics->pT() > _progenitor->maxHardPt())
return true;
// apply the soft correction
if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
if(_hardme->softMatrixElementVeto(_progenitor,particle,bb))
return true;
}
// the more general vetos
// check vs max pt for the shower
if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,bb);
switch ((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
}
if (vetoed) return true;
}
return false;
}
bool Evolver::spaceLikeDecayVetoed( const Branching & fb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction::Type type = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// apply the soft correction
if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
// veto on hardest pt in the shower
if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
// general vetos
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,fb);
switch((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
if (vetoed) return true;
}
}
return false;
}
void Evolver::hardestEmission(bool hard) {
if( ( _hardme && _hardme->hasPOWHEGCorrection()) ||
(_decayme && _decayme->hasPOWHEGCorrection())) {
if(_hardme) {
assert(hard);
if(interaction_==4) {
vector<ShowerInteraction::Type> inter(2);
inter[0] = ShowerInteraction::QCD;
inter[1] = ShowerInteraction::QED;
_hardtree = _hardme->generateHardest( currentTree(),inter );
}
else {
_hardtree = _hardme->generateHardest( currentTree(),interactions_ );
}
}
else {
assert(!hard);
_hardtree = _decayme->generateHardest( currentTree() );
}
if(!_hardtree) return;
// join up the two trees
connectTrees(currentTree(),_hardtree,hard);
}
else {
_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());
}
}
bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
HardBranchingPtr branch,
ShowerInteraction::Type type) {
int ntry=0;
do {
++ntry;
Branching fb;
unsigned int iout=0;
tcPDPtr pdata[2];
while (true) {
// no truncated shower break
if(!isTruncatedShowerON()||hardOnly()) break;
// generate emission
fb=splittingGenerator()->chooseForwardBranching(*particle,1.,type);
// no emission break
if(!fb.kinematics) break;
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// get the particle data objects
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
// find the truncated line
iout=0;
if(pdata[0]->id()!=pdata[1]->id()) {
if(pdata[0]->id()==particle->id()) iout=1;
else if (pdata[1]->id()==particle->id()) iout=2;
}
else if(pdata[0]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
// only if same interaction for forced branching
ShowerInteraction::Type type2 = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// and evolution
if(type2==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
// should do base class vetos as well
if(timeLikeVetoed(fb,particle)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
break;
}
// if no branching force trunctaed emission
if(!fb.kinematics) {
// construct the kinematics for the hard emission
ShoKinPtr showerKin=
branch->sudakov()->createFinalStateBranching(branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT());
showerKin->initialize( *particle,PPtr() );
IdList idlist(3);
idlist[0] = particle->id();
idlist[1] = branch->children()[0]->branchingParticle()->id();
idlist[2] = branch->children()[1]->branchingParticle()->id();
fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics( fb.kinematics );
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
branchingParticle()->dataPtr(),true)));
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
branchingParticle()->dataPtr(),true)));
particle->showerKinematics()->
updateChildren(particle, theChildren,fb.type);
// shower the first particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() )
timeLikeShower(theChildren[0],type,false);
}
else {
truncatedTimeLikeShower( theChildren[0],branch->children()[0],type);
}
// shower the second particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() )
timeLikeShower( theChildren[1] , type,false);
}
else {
truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
}
// that's if for old approach
if(_reconOpt==0) return true;
// branching has happened
particle->showerKinematics()->updateParent(particle, theChildren,fb.type);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
}
else return true;
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back( new_ptr( ShowerParticle( pdata[0], true ) ) );
theChildren.push_back( new_ptr( ShowerParticle( pdata[1], true ) ) );
particle->showerKinematics()->
updateChildren( particle, theChildren , fb.type);
// shower the first particle
if( iout == 1 ) truncatedTimeLikeShower( theChildren[0], branch , type );
else timeLikeShower( theChildren[0] , type,false);
// shower the second particle
if( iout == 2 ) truncatedTimeLikeShower( theChildren[1], branch , type );
else timeLikeShower( theChildren[1] , type,false);
// that's if for old approach
if(_reconOpt==0) return true;
// branching has happened
particle->showerKinematics()->updateParent(particle, theChildren,fb.type);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
}
else return true;
}
while(ntry<50);
return false;
}
bool Evolver::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
HardBranchingPtr branch,
ShowerInteraction::Type type) {
tcPDFPtr pdf;
if(ShowerHandler::currentHandler()->firstPDF().particle() == beamParticle())
pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
if(ShowerHandler::currentHandler()->secondPDF().particle() == beamParticle())
pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
Branching bb;
// parameters of the force branching
double z(0.);
HardBranchingPtr timelike;
for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
timelike = branch->children()[ix];
}
if( branch->children()[ix]->status() ==HardBranching::Incoming )
z = branch->children()[ix]->z();
}
// generate truncated branching
tcPDPtr part[2];
if(z>=0.&&z<=1.) {
while (true) {
if( !isTruncatedShowerON() || hardOnly() ) break;
bb = splittingGenerator()->chooseBackwardBranching( *particle,
beam, 1., beamParticle(),
type , pdf,freeze);
if( !bb.kinematics || bb.kinematics->scale() < branch->scale() ) {
bb = Branching();
break;
}
// particles as in Sudakov form factor
part[0] = getParticleData( bb.ids[0] );
part[1] = getParticleData( bb.ids[2] );
//is emitter anti-particle
if( particle->id() != bb.ids[1]) {
if( part[0]->CC() ) part[0] = part[0]->CC();
if( part[1]->CC() ) part[1] = part[1]->CC();
}
double zsplit = bb.kinematics->z();
// apply the vetos for the truncated shower
// if doesn't carry most of momentum
ShowerInteraction::Type type2 = bb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
if(type2==branch->sudakov()->interactionType() &&
zsplit < 0.5) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// others
if( part[0]->id() != particle->id() || // if particle changes type
bb.kinematics->pT() > progenitor()->maximumpT(type2) || // pt veto
bb.kinematics->scale() < branch->scale()) { // angular ordering veto
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// and those from the base class
if(spaceLikeVetoed(bb,particle)) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
break;
}
}
if( !bb.kinematics ) {
//do the hard emission
ShoKinPtr kinematics =
branch->sudakov()->createInitialStateBranching( branch->scale(), z, branch->phi(),
branch->children()[0]->pT() );
kinematics->initialize( *particle, beam );
// assign the splitting function and shower kinematics
particle->showerKinematics( kinematics );
// For the time being we are considering only 1->2 branching
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent =
new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
ShowerParticlePtr otherChild =
new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
true, true ) );
ShowerParticleVector theChildren;
theChildren.push_back( particle );
theChildren.push_back( otherChild );
particle->showerKinematics()->
updateParent( newParent, theChildren, branch->type());
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted=false;
if(!hardOnly()) {
if( branch->parent() ) {
emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
}
else {
emitted = spaceLikeShower( newParent, beam , type);
}
}
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[progenitor()];
kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->
updateChildren( newParent, theChildren,bb.type);
if(hardOnly()) return true;
// perform the shower of the final-state particle
if( timelike->children().empty() ) {
timeLikeShower( otherChild , type,true);
}
else {
truncatedTimeLikeShower( otherChild, timelike , type);
}
// return the emitted
return true;
}
// assign the splitting function and shower kinematics
particle->showerKinematics( bb.kinematics );
// For the time being we are considering only 1->2 branching
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
ShowerParticleVector theChildren;
theChildren.push_back( particle );
theChildren.push_back( otherChild );
particle->showerKinematics()->
updateParent( newParent, theChildren, bb.type);
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
// now reconstruct the momentum
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
bb.kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
bb.kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->
updateChildren( newParent, theChildren, bb.type);
// perform the shower of the final-state particle
timeLikeShower( otherChild , type,true);
// return the emitted
return true;
}
bool Evolver::
truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
const ShowerParticle::EvolutionScales & maxScales,
Energy minmass, HardBranchingPtr branch,
ShowerInteraction::Type type) {
Branching fb;
unsigned int iout=0;
tcPDPtr pdata[2];
while (true) {
// no truncated shower break
if(!isTruncatedShowerON()||hardOnly()) break;
fb=splittingGenerator()->chooseDecayBranching(*particle,maxScales,minmass,1.,type);
// return if no radiation
if(!fb.kinematics) break;
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// get the particle data objects
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
// find the truncated line
iout=0;
if(pdata[0]->id()!=pdata[1]->id()) {
if(pdata[0]->id()==particle->id()) iout=1;
else if (pdata[1]->id()==particle->id()) iout=2;
}
else if(pdata[0]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
ShowerInteraction::Type type2 = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
if(type2==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
// should do base class vetos as well
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle)) break;
// otherwise reset scale and continue
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// if no branching insert hard emission and continue
if(!fb.kinematics) {
// construct the kinematics for the hard emission
ShoKinPtr showerKin=
branch->sudakov()->createDecayBranching(branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT());
showerKin->initialize( *particle,PPtr() );
IdList idlist(3);
idlist[0] = particle->id();
idlist[1] = branch->children()[0]->branchingParticle()->id();
idlist[2] = branch->children()[1]->branchingParticle()->id();
assert(false);
fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics( fb.kinematics );
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
branchingParticle()->dataPtr(),true)));
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
branchingParticle()->dataPtr(),true)));
particle->showerKinematics()->
updateChildren(particle, theChildren,fb.type);
if(theChildren[0]->id()==particle->id()) {
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the space-like particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() ) spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[0],maxScales,minmass,
branch->children()[0],type);
}
// shower the second particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() ) timeLikeShower( theChildren[1] , type,true);
}
else {
truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
}
}
else {
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[1]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the space-like particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() ) spaceLikeDecayShower(theChildren[1],maxScales,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[1],maxScales,minmass,
branch->children()[1],type);
}
// shower the second particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() ) timeLikeShower( theChildren[0] , type,true);
}
else {
truncatedTimeLikeShower( theChildren[0],branch->children()[0] ,type);
}
}
return true;
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
particle->showerKinematics()->updateChildren(particle, theChildren,fb.type);
// In the case of splittings which involves coloured particles,
// set properly the colour flow of the branching.
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the first particle
truncatedSpaceLikeDecayShower(theChildren[0],maxScales,minmass,branch,type);
// shower the second particle
timeLikeShower(theChildren[1],type,true);
// branching has happened
return true;
}
bool Evolver::constructDecayTree(vector<ShowerProgenitorPtr> & particlesToShower,
ShowerInteraction::Type inter) {
Energy ptmax(-GeV);
// get the maximum pt is all ready a hard tree
if(hardTree()) {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->maximumpT(inter)>ptmax&&
particlesToShower[ix]->progenitor()->isFinalState())
ptmax = particlesToShower[ix]->maximumpT(inter);
}
}
vector<HardBranchingPtr> spaceBranchings,allBranchings;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Outgoing));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing));
}
allBranchings.push_back(newBranch);
}
else {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Decay));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
HardBranchingPtr last=newBranch;
do {
for(unsigned int ix=0;ix<last->children().size();++ix) {
if(last->children()[ix]->branchingParticle()->id()==
particlesToShower[ix]->id()) {
last = last->children()[ix];
continue;
}
}
}
while(!last->children().empty());
last->status(HardBranching::Incoming);
spaceBranchings.push_back(newBranch);
allBranchings .push_back(last);
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Incoming));
spaceBranchings.push_back(newBranch);
allBranchings .push_back(newBranch);
}
}
}
HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,inter));
// set the charge partners
ShowerParticleVector particles;
particles.push_back(spaceBranchings.back()->branchingParticle());
for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
cit!=QCDTree->branchings().end();++cit) {
if((*cit)->status()==HardBranching::Outgoing)
particles.push_back((*cit)->branchingParticle());
}
// get the partners
showerModel()->partnerFinder()->setInitialEvolutionScales(particles,true,inter,true);
// do the inverse recon
if(!showerModel()->kinematicsReconstructor()->
deconstructDecayJets(QCDTree,this,inter)) {
return false;
}
// clear the old shower
currentTree()->clear();
// set the hard tree
hardTree(QCDTree);
// set the charge partners
setEvolutionPartners(false,inter,false);
// get the particles to be showered
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
particlesToShower.clear();
// incoming particles
for(cit=currentTree()->incomingLines().begin();
cit!=currentTree()->incomingLines().end();++cit)
particlesToShower.push_back(((*cit).first));
assert(particlesToShower.size()==1);
// outgoing particles
for(cjt=currentTree()->outgoingLines().begin();
cjt!=currentTree()->outgoingLines().end();++cjt) {
particlesToShower.push_back(((*cjt).first));
if(ptmax>ZERO) particlesToShower.back()->maximumpT(ptmax,inter);
}
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
if( mit != eit) {
if(mit->second->status()==HardBranching::Outgoing)
particlesToShower[ix]->progenitor()->set5Momentum(mit->second->pVector());
}
}
return true;
}
bool Evolver::constructHardTree(vector<ShowerProgenitorPtr> & particlesToShower,
ShowerInteraction::Type inter) {
bool noEmission = true;
vector<HardBranchingPtr> spaceBranchings,allBranchings;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
noEmission = false;
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Outgoing));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing));
}
allBranchings.push_back(newBranch);
}
else {
HardBranchingPtr first,last;
if(!particlesToShower[ix]->progenitor()->parents().empty()) {
noEmission = false;
constructSpaceLikeLine(particlesToShower[ix]->progenitor(),
first,last,SudakovPtr(),
particlesToShower[ix]->original()->parents()[0]);
}
else {
first = new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Incoming));
if(particlesToShower[ix]->original()->parents().empty())
first->beam(particlesToShower[ix]->original());
else
first->beam(particlesToShower[ix]->original()->parents()[0]);
last = first;
}
spaceBranchings.push_back(first);
allBranchings.push_back(last);
}
}
if(!noEmission) {
HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,
inter));
// set the charge partners
ShowerParticleVector particles;
for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
cit!=QCDTree->branchings().end();++cit) {
particles.push_back((*cit)->branchingParticle());
}
// get the partners
showerModel()->partnerFinder()->setInitialEvolutionScales(particles,false,
inter,true);
// do the inverse recon
if(!showerModel()->kinematicsReconstructor()->
deconstructHardJets(QCDTree,this,inter))
throw Exception() << "Can't to shower deconstruction for QED shower in"
<< "QEDEvolver::showerHard" << Exception::eventerror;
// set the hard tree
hardTree(QCDTree);
}
// clear the old shower
currentTree()->clear();
// set the charge partners
setEvolutionPartners(true,inter,false);
// get the particles to be showered
particlesToShower = currentTree()->extractProgenitors();
// reset momenta
if(hardTree()) {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
if( mit != eit) {
particlesToShower[ix]->progenitor()->set5Momentum(mit->second->showerMomentum());
}
}
}
return true;
}
void Evolver::constructTimeLikeLine(tHardBranchingPtr branch,
tShowerParticlePtr particle) {
for(unsigned int ix=0;ix<particle->children().size();++ix) {
HardBranching::Status status = branch->status();
tShowerParticlePtr child =
dynamic_ptr_cast<ShowerParticlePtr>(particle->children()[ix]);
if(child->children().empty()) {
HardBranchingPtr newBranch =
new_ptr(HardBranching(child,SudakovPtr(),branch,status));
branch->addChild(newBranch);
}
else {
HardBranchingPtr newBranch =
new_ptr(HardBranching(child,child->showerKinematics()->SudakovFormFactor(),
branch,status));
constructTimeLikeLine(newBranch,child);
branch->addChild(newBranch);
}
}
// sort out the type of interaction
if(!branch->children().empty()) {
if(branch->branchingParticle()->id()==ParticleID::gamma ||
branch->children()[0]->branchingParticle()->id()==ParticleID::gamma ||
branch->children()[1]->branchingParticle()->id()==ParticleID::gamma)
branch->type(ShowerPartnerType::QED);
else {
if(branch->branchingParticle()->id()==
branch->children()[0]->branchingParticle()->id()) {
if(branch->branchingParticle()->dataPtr()->iColour()==PDT::Colour8) {
tShowerParticlePtr emittor =
branch->branchingParticle()->showerKinematics()->z()>0.5 ?
branch->children()[0]->branchingParticle() :
branch->children()[1]->branchingParticle();
if(branch->branchingParticle()->colourLine()==emittor->colourLine())
branch->type(ShowerPartnerType::QCDAntiColourLine);
else if(branch->branchingParticle()->antiColourLine()==emittor->antiColourLine())
branch->type(ShowerPartnerType::QCDColourLine);
else
assert(false);
}
else if(branch->branchingParticle()->colourLine()) {
branch->type(ShowerPartnerType::QCDColourLine);
}
else if(branch->branchingParticle()->antiColourLine()) {
branch->type(ShowerPartnerType::QCDAntiColourLine);
}
else
assert(false);
}
else if(branch->branchingParticle()->id()==ParticleID::g &&
branch->children()[0]->branchingParticle()->id()==
-branch->children()[1]->branchingParticle()->id()) {
if(branch->branchingParticle()->showerKinematics()->z()>0.5)
branch->type(ShowerPartnerType::QCDAntiColourLine);
else
branch->type(ShowerPartnerType::QCDColourLine);
}
else
assert(false);
}
}
}
void Evolver::constructSpaceLikeLine(tShowerParticlePtr particle,
HardBranchingPtr & first,
HardBranchingPtr & last,
SudakovPtr sud,PPtr beam) {
if(!particle) return;
if(!particle->parents().empty()) {
tShowerParticlePtr parent =
dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]);
SudakovPtr newSud=particle->showerKinematics()->SudakovFormFactor();
constructSpaceLikeLine(parent,first,last,newSud,beam);
}
HardBranchingPtr newBranch =
new_ptr(HardBranching(particle,sud,last,HardBranching::Incoming));
newBranch->beam(beam);
if(!first) {
first=newBranch;
last =newBranch;
return;
}
last->addChild(newBranch);
tShowerParticlePtr timeChild =
dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]->children()[1]);
HardBranchingPtr timeBranch;
if(!timeChild->children().empty()) {
timeBranch =
new_ptr(HardBranching(timeChild,
timeChild->showerKinematics()->SudakovFormFactor(),
last,HardBranching::Outgoing));
constructTimeLikeLine(timeBranch,timeChild);
}
else {
timeBranch =
new_ptr(HardBranching(timeChild,SudakovPtr(),last,HardBranching::Outgoing));
}
last->addChild(timeBranch);
// sort out the type
if(last->branchingParticle() ->id() == ParticleID::gamma ||
newBranch->branchingParticle() ->id() == ParticleID::gamma ||
timeBranch->branchingParticle()->id() == ParticleID::gamma) {
last->type(ShowerPartnerType::QED);
}
else if(last->branchingParticle()->id()==newBranch->branchingParticle()->id()) {
if(last->branchingParticle()->id()==ParticleID::g) {
if(last->branchingParticle()->colourLine()==
newBranch->branchingParticle()->colourLine()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else {
last->type(ShowerPartnerType::QCDColourLine);
}
}
else if(last->branchingParticle()->hasColour()) {
last->type(ShowerPartnerType::QCDColourLine);
}
else if(last->branchingParticle()->hasAntiColour()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else
assert(false);
}
else if(newBranch->branchingParticle()->id()==ParticleID::g) {
if(last->branchingParticle()->hasColour()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else if(last->branchingParticle()->hasAntiColour()) {
last->type(ShowerPartnerType::QCDColourLine);
}
else
assert(false);
}
else if(newBranch->branchingParticle()->hasColour()) {
last->type(ShowerPartnerType::QCDColourLine);
}
else if(newBranch->branchingParticle()->hasAntiColour()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else {
assert(false);
}
last=newBranch;
}
void Evolver::connectTrees(ShowerTreePtr showerTree,
HardTreePtr hardTree, bool hard ) {
ShowerParticleVector particles;
// find the Sudakovs
for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
// Sudakovs for ISR
if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) {
++_nis;
IdList br(3);
br[0] = (**cit).parent()->branchingParticle()->id();
br[1] = (**cit). branchingParticle()->id();
br[2] = (**cit).parent()->children()[0]==*cit ?
(**cit).parent()->children()[1]->branchingParticle()->id() :
(**cit).parent()->children()[0]->branchingParticle()->id();
BranchingList branchings = splittingGenerator()->initialStateBranchings();
if(br[1]<0&&br[0]==br[1]) {
br[0] = abs(br[0]);
br[1] = abs(br[1]);
}
else if(br[1]<0) {
br[1] = -br[1];
br[2] = -br[2];
}
long index = abs(br[1]);
SudakovPtr sudakov;
for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
cjt != branchings.upper_bound(index); ++cjt ) {
IdList ids = cjt->second.second;
if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
sudakov=cjt->second.first;
break;
}
}
if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
<< "Evolver::connectTrees() for ISR"
<< Exception::runerror;
(**cit).parent()->sudakov(sudakov);
}
// Sudakovs for FSR
else if(!(**cit).children().empty()) {
++_nfs;
IdList br(3);
br[0] = (**cit) .branchingParticle()->id();
br[1] = (**cit).children()[0]->branchingParticle()->id();
br[2] = (**cit).children()[1]->branchingParticle()->id();
BranchingList branchings = splittingGenerator()->finalStateBranchings();
if(br[0]<0) {
br[0] = abs(br[0]);
br[1] = abs(br[1]);
br[2] = abs(br[2]);
}
long index = br[0];
SudakovPtr sudakov;
for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
cjt != branchings.upper_bound(index); ++cjt ) {
IdList ids = cjt->second.second;
if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
sudakov=cjt->second.first;
break;
}
}
if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
<< "Evolver::connectTrees()"
<< Exception::runerror;
(**cit).sudakov(sudakov);
}
}
// calculate the evolution scale
for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
particles.push_back((*cit)->branchingParticle());
}
showerModel()->partnerFinder()->
setInitialEvolutionScales(particles,!hard,hardTree->interaction(),
!hardTree->partnersSet());
hardTree->partnersSet(true);
// inverse reconstruction
if(hard)
showerModel()->kinematicsReconstructor()->
deconstructHardJets(hardTree,ShowerHandler::currentHandler()->evolver(),
hardTree->interaction());
else
showerModel()->kinematicsReconstructor()->
deconstructDecayJets(hardTree,ShowerHandler::currentHandler()->evolver(),
hardTree->interaction());
// now reset the momenta of the showering particles
vector<ShowerProgenitorPtr> particlesToShower;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=showerTree->incomingLines().begin();
cit!=showerTree->incomingLines().end();++cit )
particlesToShower.push_back(cit->first);
// extract the showering particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=showerTree->outgoingLines().begin();
cit!=showerTree->outgoingLines().end();++cit )
particlesToShower.push_back(cit->first);
// match them
vector<bool> matched(particlesToShower.size(),false);
for(set<HardBranchingPtr>::const_iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
Energy2 dmin( 1e30*GeV2 );
int iloc(-1);
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(matched[ix]) continue;
if( (**cit).branchingParticle()->id() != particlesToShower[ix]->progenitor()->id() ) continue;
if( (**cit).branchingParticle()->isFinalState() !=
particlesToShower[ix]->progenitor()->isFinalState() ) continue;
Energy2 dtest =
sqr( particlesToShower[ix]->progenitor()->momentum().x() - (**cit).showerMomentum().x() ) +
sqr( particlesToShower[ix]->progenitor()->momentum().y() - (**cit).showerMomentum().y() ) +
sqr( particlesToShower[ix]->progenitor()->momentum().z() - (**cit).showerMomentum().z() ) +
sqr( particlesToShower[ix]->progenitor()->momentum().t() - (**cit).showerMomentum().t() );
// add mass difference for identical particles (e.g. Z0 Z0 production)
dtest += 1e10*sqr(particlesToShower[ix]->progenitor()->momentum().m()-
(**cit).showerMomentum().m());
if( dtest < dmin ) {
iloc = ix;
dmin = dtest;
}
}
if(iloc<0) throw Exception() << "Failed to match shower and hard trees in Evolver::hardestEmission"
<< Exception::eventerror;
particlesToShower[iloc]->progenitor()->set5Momentum((**cit).showerMomentum());
matched[iloc] = true;
}
// correction boosts for daughter trees
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit = showerTree->treelinks().begin();
tit != showerTree->treelinks().end();++tit) {
ShowerTreePtr decayTree = tit->first;
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = decayTree->incomingLines().begin();
// reset the momentum of the decay particle
Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
Lorentz5Momentum newMomentum = tit->second.second->momentum();
LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
boost.boost (-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
decayTree->transform(boost,true);
}
}
void Evolver::doShowering(bool hard) {
// order of the interactions
bool showerOrder(true);
// zero number of emissions
_nis = _nfs = 0;
// extract particles to shower
vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
// setup the maximum scales for the shower
setupMaximumScales(particlesToShower);
// specifc stuff for hard processes and decays
Energy minmass(ZERO), mIn(ZERO);
// hard process generate the intrinsic p_T once and for all
if(hard) {
generateIntrinsicpT(particlesToShower);
}
// decay compute the minimum mass of the final-state
else {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
if(particlesToShower[ix]->progenitor()->dataPtr()->stable())
minmass += particlesToShower[ix]->progenitor()->dataPtr()->constituentMass();
else
minmass += particlesToShower[ix]->progenitor()->mass();
}
else {
mIn = particlesToShower[ix]->progenitor()->mass();
}
}
- // throw exception if decay can'ty happen
+ // throw exception if decay can't happen
if ( minmass > mIn ) {
throw Exception() << "Evolver.cc: Mass of decaying particle is "
<< "below constituent masses of decay products."
<< Exception::eventerror;
}
}
// check if interactions in right order
if(hardTree() && interaction_!=4 &&
hardTree()->interaction()!=interactions_[0]) {
assert(interactions_.size()==2);
showerOrder = false;
swap(interactions_[0],interactions_[1]);
}
// loop over possible interactions
for(unsigned int inter=0;inter<interactions_.size();++inter) {
// set up for second pass if required
if(inter!=0) {
// zero intrinsic pt so only added first time round
intrinsicpT().clear();
// construct the tree and throw veto if not possible
if(!(hard ?
constructHardTree (particlesToShower,interactions_[inter]) :
constructDecayTree(particlesToShower,interactions_[inter])))
throw InteractionVeto();
}
// main shower loop
unsigned int ntry(0);
bool reconstructed = false;
do {
// clear results of last attempt if needed
if(ntry!=0) {
currentTree()->clear();
setEvolutionPartners(hard,interactions_[inter],false);
_nis = _nfs = 0;
}
// generate the shower
// pick random starting point
unsigned int istart=UseRandom::irnd(particlesToShower.size());
unsigned int istop = particlesToShower.size();
// loop over particles with random starting point
for(unsigned int ix=istart;ix<=istop;++ix) {
if(ix==particlesToShower.size()) {
if(istart!=0) {
istop = istart-1;
ix=0;
}
else break;
}
// extract the progenitor
progenitor(particlesToShower[ix]);
// final-state radiation
if(progenitor()->progenitor()->isFinalState()) {
if(!isFSRadiationON()) continue;
// perform shower
progenitor()->hasEmitted(startTimeLikeShower(interactions_[inter]));
}
// initial-state radiation
else {
if(!isISRadiationON()) continue;
// hard process
if(hard) {
// get the PDF
setBeamParticle(_progenitor->beam());
assert(beamParticle());
// perform the shower
// set the beam particle
tPPtr beamparticle=progenitor()->original();
if(!beamparticle->parents().empty())
beamparticle=beamparticle->parents()[0];
// generate the shower
progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
interactions_[inter]));
}
// decay
else {
// skip colour and electrically neutral particles
if(!progenitor()->progenitor()->dataPtr()->coloured() &&
!progenitor()->progenitor()->dataPtr()->charged()) {
progenitor()->hasEmitted(false);
continue;
}
// perform shower
// set the scales correctly. The current scale is the maximum scale for
// emission not the starting scale
ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
if(progenitor()->progenitor()->dataPtr()->charged()) {
progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
}
if(progenitor()->progenitor()->hasColour()) {
progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass();
}
if(progenitor()->progenitor()->hasAntiColour()) {
progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
}
// perform the shower
progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
interactions_[inter]));
}
}
}
// do the kinematic reconstruction, checking if it worked
reconstructed = hard ?
showerModel()->kinematicsReconstructor()->
reconstructHardJets (currentTree(),intrinsicpT(),interactions_[inter]) :
showerModel()->kinematicsReconstructor()->
reconstructDecayJets(currentTree(),interactions_[inter]);
}
while(!reconstructed&&maximumTries()>++ntry);
// check if failed to generate the shower
if(ntry==maximumTries()) {
if(hard)
throw ShowerHandler::ShowerTriesVeto(ntry);
else
throw Exception() << "Failed to generate the shower after "
<< ntry << " attempts in Evolver::showerDecay()"
<< Exception::eventerror;
}
}
// tree has now showered
_currenttree->hasShowered(true);
if(!showerOrder) swap(interactions_[0],interactions_[1]);
hardTree(HardTreePtr());
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Mon, Jan 20, 11:27 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242847
Default Alt Text
(134 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment