Page MenuHomeHEPForge
Files F9501533

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Lepton/MEee2ff.cc b/MatrixElement/Lepton/MEee2ff.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Lepton/MEee2ff.cc
@@ -0,0 +1,776 @@
+// -*- C++ -*-
+//
+// MEee2ff.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
+// Copyright (C) 2002-2019 The Herwig Collaboration
+//
+// Herwig is licenced under version 3 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the MEee2ff class.
+//
+
+#include "MEee2ff.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/Helicity/WaveFunction/VectorWaveFunction.h"
+#include "ThePEG/Handlers/StandardXComb.h"
+#include "Herwig/MatrixElement/HardVertex.h"
+#include "ThePEG/PDF/PolarizedBeamParticleData.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include <numeric>
+#include "Herwig/Shower/RealEmissionProcess.h"
+
+using namespace Herwig;
+
+void MEee2ff::getDiagrams() const {
+ // specific the diagrams
+ tcPDPtr f, fbar;
+ if(MEee2ff::incoming_==0) {
+ f = getParticleData(ParticleID::eplus);
+ fbar = getParticleData(ParticleID::eminus);
+ }
+ else if(MEee2ff::incoming_==1) {
+ f = getParticleData(ParticleID::muplus);
+ fbar = getParticleData(ParticleID::muminus);
+ }
+ else
+ assert(false);
+ // setup the processes
+ for( int i =11;i<=16;++i) {
+ if(allowed_==0 || (allowed_==1 && i%2==1) || (allowed_==2&&i==11)
+ || (allowed_==3&&i==13) || (allowed_==4&&i==15)) {
+ tcPDPtr lm = getParticleData(i);
+ tcPDPtr lp = lm->CC();
+ add(new_ptr((Tree2toNDiagram(2), f, fbar, 1, gamma_, 3, lm, 3, lp, -1)));
+ add(new_ptr((Tree2toNDiagram(2), f, fbar, 1, Z0_, 3, lm, 3, lp, -2)));
+ }
+ }
+}
+
+Energy2 MEee2ff::scale() const {
+ return sHat();
+}
+
+unsigned int MEee2ff::orderInAlphaS() const {
+ return 0;
+}
+
+unsigned int MEee2ff::orderInAlphaEW() const {
+ return 2;
+}
+
+Selector<MEBase::DiagramIndex>
+MEee2ff::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 *>
+MEee2ff::colourGeometries(tcDiagPtr) const {
+ static ColourLines ctST(" ");
+ Selector<const ColourLines *> sel;
+ sel.insert(1.0, &ctST);
+ return sel;
+}
+
+void MEee2ff::persistentOutput(PersistentOStream & os) const {
+ os << allowed_ << FFZVertex_ << FFPVertex_ << gamma_ << Z0_
+ << alphaQED_ << ounit(pTmin_,GeV) << preFactor_
+ << helicity_ << incoming_;
+}
+
+void MEee2ff::persistentInput(PersistentIStream & is, int) {
+ is >> allowed_ >> FFZVertex_ >> FFPVertex_ >> gamma_ >> Z0_
+ >> alphaQED_ >> iunit(pTmin_,GeV) >> preFactor_
+ >> helicity_ >> incoming_;
+}
+
+// *** 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<MEee2ff,HwMEBase>
+describeMEee2ff("Herwig::MEee2ff", "HwMELepton.so");
+
+void MEee2ff::Init() {
+
+ static ClassDocumentation<MEee2ff> documentation
+ ("The MEee2ff class implements the matrix element for"
+ "e+e- to leptons via Z and photon exchange using helicity amplitude"
+ "techniques");
+
+ static Switch<MEee2ff,int> interfaceallowed
+ ("Outgoing",
+ "Allowed outgoing leptons",
+ &MEee2ff::allowed_, 0, false, false);
+ static SwitchOption interfaceallowedAll
+ (interfaceallowed,
+ "All",
+ "Allow all leptons as outgoing particles",
+ 0);
+ static SwitchOption interfaceallowedCharged
+ (interfaceallowed,
+ "Charged",
+ "Only charged leptons as outgoing particles",
+ 1);
+ static SwitchOption interfaceallowedElectron
+ (interfaceallowed,
+ "Electron",
+ "Only the electron and positron as outgoing leptons",
+ 2);
+ static SwitchOption interfaceallowedMuon
+ (interfaceallowed,
+ "Muon",
+ "Only muons as outgoing particles",
+ 3);
+ static SwitchOption interfaceallowedTau
+ (interfaceallowed,
+ "Tau",
+ "Only taus as outgoing particles",
+ 4);
+ static SwitchOption interfaceallowedElectronNeutrino
+ (interfaceallowed,
+ "ElectronNeutrino",
+ "Only electron neutrino as outgoing particles",
+ 5);
+ static SwitchOption interfaceallowedMuonNeutrino
+ (interfaceallowed,
+ "MuonNeutrino",
+ "Only muon neutrino as outgoing particles",
+ 6);
+ static SwitchOption interfaceallowedTauNeutrino
+ (interfaceallowed,
+ "TauNeutrino",
+ "Only tau neutrino as outgoing particles",
+ 7);
+
+ static Parameter<MEee2ff,Energy> interfacepTMin
+ ("pTMin",
+ "Minimum pT for hard radiation",
+ &MEee2ff::pTmin_, GeV, 1.0*GeV, 0.001*GeV, 10.0*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<MEee2ff,double> interfacePrefactor
+ ("Prefactor",
+ "Prefactor for the overestimate of the emission probability",
+ &MEee2ff::preFactor_, 6.0, 1.0, 100.0,
+ false, false, Interface::limited);
+
+ static Reference<MEee2ff,ShowerAlpha> interfaceEMCoupling
+ ("AlphaQED",
+ "Pointer to the object to calculate the EM coupling for the correction",
+ &MEee2ff::alphaQED_, false, false, true, false, false);
+
+ static Switch<MEee2ff,int> interfacehelicity
+ ("Helicity",
+ "Helicity of the outgoing fermions",
+ &MEee2ff::helicity_, 0, false, false);
+ static SwitchOption interfacehelicityall
+ (interfacehelicity,
+ "Any",
+ "Allow all helicity configurations for the outgoing particles",
+ 0);
+ static SwitchOption interfacehelicityll
+ (interfacehelicity,
+ "LL",
+ "Allow only left-left helicity configuration for the outgoing particles",
+ 1);
+ static SwitchOption interfacehelicityrr
+ (interfacehelicity,
+ "RR",
+ "Allow only right-right helicity configuration for the outgoing particles",
+ 2);
+
+ static Switch<MEee2ff,int> interfaceincoming
+ ("Incoming",
+ "Incoming fermionc beam",
+ &MEee2ff::incoming_, 0, false, false);
+ static SwitchOption interfaceincomingElectron
+ (interfaceincoming,
+ "Electron",
+ "Incoming elecron-positron beam",
+ 0);
+ static SwitchOption interfaceincomingMuon
+ (interfaceincoming,
+ "Muon",
+ "Incoming muon-antimuon beam",
+ 1);
+}
+
+double MEee2ff::me2() const {
+ return loME(mePartonData(),rescaledMomenta(),true);
+}
+
+double MEee2ff::loME(const vector<cPDPtr> & partons,
+ const vector<Lorentz5Momentum> & momenta,
+ bool first) const {
+ vector<SpinorWaveFunction> fin,aout;
+ vector<SpinorBarWaveFunction> ain,fout;
+ SpinorWaveFunction ein (momenta[0],partons[0],incoming);
+ SpinorBarWaveFunction pin (momenta[1],partons[1],incoming);
+ SpinorBarWaveFunction ilmout(momenta[2],partons[2],outgoing);
+ SpinorWaveFunction ilpout(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 );
+ ilmout.reset(ix);fout.push_back(ilmout);
+ ilpout.reset(ix);aout.push_back(ilpout);
+ }
+ // 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;
+}
+
+ProductionMatrixElement MEee2ff::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(sHat(),1,Z0_,fin[inhel1],ain[inhel2]);
+ // intermediate photon
+ interG = FFPVertex_->evaluate(sHat(),1,gamma_,fin[inhel1],ain[inhel2]);
+ for(outhel1=0;outhel1<2;++outhel1) {
+ for(outhel2=0;outhel2<2;++outhel2) {
+ // in case LL or RR helicity configurations are requested
+ if(MEee2ff::helicity_==1 && (outhel1==1 || outhel2==1)) //LL
+ continue;
+ if(MEee2ff::helicity_==2 && (outhel1==0 || outhel2==0)) //RR
+ continue;
+ // first the Z exchange diagram
+ diag1 = FFZVertex_->evaluate(sHat(),aout[outhel2],fout[outhel1],
+ interZ);
+ // then the photon exchange diagram
+ diag2 = FFPVertex_->evaluate(sHat(),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;
+ }
+ }
+ }
+ }
+ // results
+ 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]);
+ }
+ cont = total[2];
+ BW = total[1];
+ me = total[0];
+ return output;
+}
+
+void MEee2ff::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;
+ 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);
+ // 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 MEee2ff::doinit() {
+ HwMEBase::doinit();
+ // set the particle data objects
+ Z0_=getParticleData(ThePEG::ParticleID::Z0);
+ gamma_=getParticleData(ThePEG::ParticleID::gamma);
+ // 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 "
+ << "MEee2ff::doinit() the Herwig"
+ << " version must be used"
+ << Exception::runerror;
+ FFZVertex_ = hwsm->vertexFFZ();
+ FFPVertex_ = hwsm->vertexFFP();
+}
+
+void MEee2ff::rebind(const TranslationMap & trans) {
+ FFZVertex_ = trans.translate(FFZVertex_);
+ FFPVertex_ = trans.translate(FFPVertex_);
+ Z0_ = trans.translate(Z0_);
+ gamma_ = trans.translate(gamma_);
+ HwMEBase::rebind(trans);
+}
+
+IVector MEee2ff::getReferences() {
+ IVector ret = HwMEBase::getReferences();
+ ret.push_back(FFZVertex_);
+ ret.push_back(FFPVertex_);
+ ret.push_back(Z0_ );
+ ret.push_back(gamma_ );
+ return ret;
+}
+
+RealEmissionProcessPtr MEee2ff::generateHardest(RealEmissionProcessPtr born,
+ ShowerInteraction inter) {
+ // check if QED switched on
+ if(inter==ShowerInteraction::QCD) return RealEmissionProcessPtr();
+ // generate the momenta for the hard emission
+ vector<Lorentz5Momentum> emission;
+ unsigned int iemit,ispect;
+ Energy pTveto = generateHard(born,emission,iemit,ispect,false);
+ // check if charged
+ if(!partons_[2]->charged()) return RealEmissionProcessPtr();
+ // maximum pT of emission
+ if(pTveto<=ZERO) {
+ born->pT()[ShowerInteraction::QED] = pTmin_;
+ return born;
+ }
+ else {
+ born->pT()[ShowerInteraction::QED] = pTveto;
+ }
+ born->interaction(ShowerInteraction::QED);
+ // get the quark and antiquark
+ ParticleVector qq;
+ for(unsigned int ix=0;ix<2;++ix) qq.push_back(born->bornOutgoing()[ix]);
+ bool order = qq[0]->id()>0;
+ if(!order) swap(qq[0],qq[1]);
+ // create the new quark, antiquark and gluon
+ PPtr newq = qq[0]->dataPtr()->produceParticle(emission[2]);
+ PPtr newa = qq[1]->dataPtr()->produceParticle(emission[3]);
+ PPtr newg = gamma_->produceParticle(emission[4]);
+ // create the output real emission process
+ for(unsigned int ix=0;ix<born->bornIncoming().size();++ix) {
+ born->incoming().push_back(born->bornIncoming()[ix]->dataPtr()->
+ produceParticle(born->bornIncoming()[ix]->momentum()));
+ }
+ if(order) {
+ born->outgoing().push_back(newq);
+ born->outgoing().push_back(newa);
+ }
+ else {
+ born->outgoing().push_back(newa);
+ born->outgoing().push_back(newq);
+ swap(iemit,ispect);
+ }
+ born->outgoing().push_back(newg);
+ // set emitter and spectator
+ born->emitter (iemit);
+ born->spectator(ispect);
+ born->emitted(4);
+ return born;
+}
+
+double MEee2ff::meRatio(vector<cPDPtr> partons,
+ vector<Lorentz5Momentum> momenta,
+ unsigned int iemitter, 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);
+ }
+ InvEnergy2 ratio = realME(partons,momenta)
+ *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;
+}
+
+InvEnergy2 MEee2ff::realME(const vector<cPDPtr> & partons,
+ const vector<Lorentz5Momentum> & momenta) 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 photon(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);
+ photon.reset(2*ix);
+ gout.push_back(photon);
+ }
+ 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 =
+ FFPVertex_->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 =
+ FFPVertex_->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 and charge
+ total /= norm(FFPVertex_->norm())*
+ sqr(double(mePartonData()[2]->iCharge())/3.);
+ // return the total
+ return total*UnitRemoval::InvE2;
+}
+
+Energy MEee2ff::generateHard(RealEmissionProcessPtr born,
+ vector<Lorentz5Momentum> & emmision,
+ unsigned int & iemit, unsigned int & ispect,
+ bool applyVeto) {
+ // get the momenta of the incoming and outgoing particles
+ // incoming
+ tPPtr em = born->bornIncoming()[0];
+ tPPtr ep = born->bornIncoming()[1];
+ if(em->id()<0) swap(em,ep);
+ // outgoing
+ tPPtr qk = born->bornOutgoing()[0];
+ tPPtr qb = born->bornOutgoing()[1];
+ if(qk->id()<0) swap(qk,qb);
+ // extract the momenta
+ loMomenta_.clear();
+ loMomenta_.push_back(em->momentum());
+ loMomenta_.push_back(ep->momentum());
+ loMomenta_.push_back(qk->momentum());
+ loMomenta_.push_back(qb->momentum());
+ // and ParticleData objects
+ partons_.resize(5);
+ partons_[0]=em->dataPtr();
+ partons_[1]=ep->dataPtr();
+ partons_[2]=qk->dataPtr();
+ partons_[3]=qb->dataPtr();
+ partons_[4]=gamma_;
+ // 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
+ double ymax = acosh(pTmax/pTmin_);
+ double 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_) {
+ 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,false);
+ // coupling piece
+ 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_)
+ pT[ix] = -GeV;
+ }
+ // now pick the emmision with highest pT
+ Energy pTemit(ZERO);
+ // no emission
+ if(pT[0]<ZERO&&pT[1]<ZERO) return -GeV;
+ // which one emitted
+ if(pT[0]>pT[1]) {
+ iemit = 2;
+ ispect = 3;
+ pTemit = pT[0];
+ if(UseRandom::rnd()<contrib[0][0]/(contrib[0][0]+contrib[0][1]))
+ emmision = realMomenta[0][0];
+ else
+ emmision = realMomenta[0][1];
+ }
+ else {
+ iemit = 3;
+ ispect = 2;
+ pTemit = pT[1];
+ if(UseRandom::rnd()<contrib[1][0]/(contrib[1][0]+contrib[1][1]))
+ emmision = realMomenta[1][0];
+ else
+ emmision = realMomenta[1][1];
+ }
+ // return pT of emmision
+ return pTemit;
+}
diff --git a/MatrixElement/Lepton/MEee2ff.h b/MatrixElement/Lepton/MEee2ff.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Lepton/MEee2ff.h
@@ -0,0 +1,373 @@
+// -*- C++ -*-
+//
+// MEee2ff.h is a part of Herwig - A multi-purpose Monte Carlo event generator
+// Copyright (C) 2002-2019 The Herwig Collaboration
+//
+// Herwig is licenced under version 3 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+#ifndef HERWIG_MEee2ff_H
+#define HERWIG_MEee2ff_H
+//
+// This is the declaration of the MEee2ff class.
+//
+
+#include "Herwig/MatrixElement/HwMEBase.h"
+#include "Herwig/Models/StandardModel/StandardModel.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "Herwig/MatrixElement/ProductionMatrixElement.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
+#include "Herwig/Shower/ShowerAlpha.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The MEee2ff class provides the matrix element for
+ * \f$e^+e^-\to\ell^+\ell^-\f$. N.B. for the production of \f$e^+e^-\f$
+ * only the \f$s\f$-channel Z and photon diagrams are included.
+ *
+ * @see \ref MEee2ffInterfaces "The interfaces"
+ * defined for MEee2ff.
+ */
+class MEee2ff: public HwMEBase {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ MEee2ff() : allowed_(0), helicity_(0), incoming_(0), pTmin_(GeV),
+ preFactor_(6.) {
+ massOption(vector<unsigned int>(2,1));
+ }
+
+ /**
+ * Members for hard corrections to the emission of QCD radiation
+ */
+ //@{
+ /**
+ * Has a POWHEG style correction
+ */
+ virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
+
+ /**
+ * Has an old fashioned ME correction
+ */
+ virtual bool hasMECorrection() {return false;}
+
+ /**
+ * Apply the POWHEG style correction
+ */
+ virtual RealEmissionProcessPtr generateHardest(RealEmissionProcessPtr,
+ ShowerInteraction);
+ //@}
+
+public:
+
+ /** @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 \ell^+ \ell^-\f$.
+ * @param partons The incoming and outgoing particles
+ * @param momenta The momenta of the incoming and outgoing particles
+ */
+ 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
+ * \[\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)}}\]
+ * 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
+ */
+ double meRatio(vector<cPDPtr> partons,
+ vector<Lorentz5Momentum> momenta,
+ unsigned int iemittor,
+ bool subtract=false) const;
+
+ /**
+ * Calculate the matrix element for \f$e^-e^-\to q \bar q g$.
+ * @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) const;
+
+ /**
+ * Generate the momenta for a hard configuration
+ */
+ Energy generateHard(RealEmissionProcessPtr tree,
+ vector<Lorentz5Momentum> & emission,
+ unsigned int & iemit, unsigned int & ispect,
+ bool applyVeto);
+
+
+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_;}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ MEee2ff & operator=(const MEee2ff &) = delete;
+
+private:
+
+ /**
+ * 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 particle data object for the Z
+ */
+ PDPtr Z0_;
+
+ /**
+ * Pointer to the particle data object for the photon
+ */
+ PDPtr gamma_;
+
+ /**
+ * The allowed outgoing
+ */
+ int allowed_;
+
+ /**
+ * The allowed helicity configuration
+ */
+ int helicity_;
+
+ /**
+ * Set the incoming lepton beam, either electron or meuon beams
+ */
+ int incoming_;
+
+ /**
+ * Pointer to the EM coupling
+ */
+ ShowerAlphaPtr alphaQED_;
+
+ /**
+ * Variables for the POWHEG style corrections
+ */
+ //@{
+ /**
+ * The cut off on pt, assuming massless quarks.
+ */
+ Energy pTmin_;
+
+ /**
+ * 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_MEee2ff_H */
diff --git a/MatrixElement/Lepton/Makefile.am b/MatrixElement/Lepton/Makefile.am
--- a/MatrixElement/Lepton/Makefile.am
+++ b/MatrixElement/Lepton/Makefile.am
@@ -1,10 +1,11 @@
pkglib_LTLIBRARIES = HwMELepton.la
HwMELepton_la_SOURCES = \
MEee2gZ2qq.h MEee2gZ2qq.cc\
MEee2gZ2ll.h MEee2gZ2ll.cc\
MEee2ZH.h MEee2ZH.cc\
MEee2HiggsVBF.h MEee2HiggsVBF.cc \
MEee2VV.h MEee2VV.cc \
MEee2VectorMeson.h MEee2VectorMeson.cc \
-MEee2Higgs2SM.h MEee2Higgs2SM.cc
+MEee2Higgs2SM.h MEee2Higgs2SM.cc \
+MEee2ff.h MEee2ff.cc
HwMELepton_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 8:0:1
diff --git a/src/defaults/MatrixElements.in b/src/defaults/MatrixElements.in
--- a/src/defaults/MatrixElements.in
+++ b/src/defaults/MatrixElements.in
@@ -1,250 +1,254 @@
# -*- ThePEG-repository -*-
##############################################################################
# Setup of default matrix elements.
#
# Only one ME is activated by default, but this file lists
# some alternatives. All available MEs can be found in the
# 'include/Herwig/MatrixElements' subdirectory of your Herwig
# installation.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Instead of editing this file directly, you should reset
# the matrix elements in your own input files:
#
# - create your custom SubProcessHandler
# - insert the MEs you need
# - set your SubProcessHandler instead of the default (see HerwigDefaults.in)
##############################################################################
mkdir /Herwig/MatrixElements
cd /Herwig/MatrixElements
library HwMELepton.so
library HwMEHadron.so
library HwMEDIS.so
############################################################
# e+e- matrix elements
############################################################
# e+e- > q qbar
create Herwig::MEee2gZ2qq MEee2gZ2qq
newdef MEee2gZ2qq:MinimumFlavour 1
newdef MEee2gZ2qq:MaximumFlavour 5
newdef MEee2gZ2qq:AlphaQCD /Herwig/Shower/AlphaQCD
newdef MEee2gZ2qq:AlphaQED /Herwig/Shower/AlphaQED
# e+e- -> l+l-
create Herwig::MEee2gZ2ll MEee2gZ2ll
newdef MEee2gZ2ll:Allowed Charged
set MEee2gZ2ll:AlphaQED /Herwig/Shower/AlphaQED
+# e+e- -> l+l-
+create Herwig::MEee2ff MEee2ff
+set MEee2ff:AlphaQED /Herwig/Shower/AlphaQED
+
# e+e- -> W+W- ZZ
create Herwig::MEee2VV MEee2VV
# e+e- -> ZH
create Herwig::MEee2ZH MEee2ZH
newdef MEee2ZH:Coupling /Herwig/Shower/AlphaQCD
# e+e- -> e+e-H/nu_enu_ebarH
create Herwig::MEee2HiggsVBF MEee2HiggsVBF
############################################################
# NLO (POWHEG e+e- matrix elements
############################################################
library HwPowhegMELepton.so
create Herwig::MEee2gZ2qqPowheg PowhegMEee2gZ2qq
newdef PowhegMEee2gZ2qq:MinimumFlavour 1
newdef PowhegMEee2gZ2qq:MaximumFlavour 5
newdef PowhegMEee2gZ2qq:AlphaQCD /Herwig/Shower/AlphaQCD
newdef PowhegMEee2gZ2qq:AlphaQED /Herwig/Shower/AlphaQED
create Herwig::MEee2gZ2llPowheg PowhegMEee2gZ2ll
newdef PowhegMEee2gZ2ll:Allowed Charged
set PowhegMEee2gZ2ll:AlphaQED /Herwig/Shower/AlphaQED
############################################################
# hadron-hadron matrix elements
############################################################
###################################
# Electroweak processes
###################################
# q qbar -> gamma/Z -> l+l-
create Herwig::MEqq2gZ2ff MEqq2gZ2ff
newdef MEqq2gZ2ff:Process 3
newdef MEqq2gZ2ff:Coupling /Herwig/Shower/AlphaQCD
# q qbar to W -> l nu
create Herwig::MEqq2W2ff MEqq2W2ff
newdef MEqq2W2ff:Process 2
newdef MEqq2W2ff:Coupling /Herwig/Shower/AlphaQCD
# W+jet
create Herwig::MEPP2WJet MEWJet
newdef MEWJet:WDecay Leptons
# Z+jet
create Herwig::MEPP2ZJet MEZJet
newdef MEZJet:ZDecay ChargedLeptons
# PP->WW/WZ/ZZ
create Herwig::MEPP2VV MEPP2VV
# PP->WZ gamma
create Herwig::MEPP2VGamma MEPP2VGamma
###################################
# Photon and jet processes
###################################
# qqbar/gg -> gamma gamma
create Herwig::MEPP2GammaGamma MEGammaGamma
# hadron-hadron to gamma+jet
create Herwig::MEPP2GammaJet MEGammaJet
# QCD 2-to-2
create Herwig::MEQCD2to2 MEQCD2to2
# MinBias
create Herwig::MEMinBias MEMinBias
newdef MEMinBias:csNorm 0.01
newdef MEMinBias:Scale 2.0
###################################
# Heavy Quark
###################################
# qqbar/gg -> t tbar
create Herwig::MEPP2QQ MEHeavyQuark
create Herwig::MEPP2SingleTop MESingleTopTChannel
set MESingleTopTChannel:Process tChannel
create Herwig::MEPP2SingleTop MESingleTopSChannel
set MESingleTopSChannel:Process sChannel
create Herwig::MEPP2SingleTop MESingleTopTW
set MESingleTopTW:Process tW
###################################
# Higgs processes
###################################
# hadron-hadron to higgs
create Herwig::MEPP2Higgs MEHiggs
newdef MEHiggs:ShapeScheme MassGenerator
newdef MEHiggs:Process gg
newdef MEHiggs:Coupling /Herwig/Shower/AlphaQCD
# hadron-hadron to higgs+jet
create Herwig::MEPP2HiggsJet MEHiggsJet
# PP->ZH
create Herwig::MEPP2ZH MEPP2ZH
newdef MEPP2ZH:Coupling /Herwig/Shower/AlphaQCD
# PP->WH
create Herwig::MEPP2WH MEPP2WH
newdef MEPP2WH:Coupling /Herwig/Shower/AlphaQCD
# PP -> Higgs via VBF
create Herwig::MEPP2HiggsVBF MEPP2HiggsVBF
newdef MEPP2HiggsVBF:ShowerAlphaQCD /Herwig/Shower/AlphaQCD
# PP -> t tbar Higgs
create Herwig::MEPP2QQHiggs MEPP2ttbarH
newdef MEPP2ttbarH:QuarkType Top
# PP -> b bbar Higgs
create Herwig::MEPP2QQHiggs MEPP2bbbarH
newdef MEPP2bbbarH:QuarkType Bottom
##########################################################
# Hadron-Hadron NLO matrix elements in the Powheg scheme
##########################################################
library HwPowhegMEHadron.so
# q qbar -> gamma/Z -> l+l-
create Herwig::MEqq2gZ2ffPowheg PowhegMEqq2gZ2ff
newdef PowhegMEqq2gZ2ff:Process 3
newdef PowhegMEqq2gZ2ff:Coupling /Herwig/Shower/AlphaQCD
# q qbar to W -> l nu
create Herwig::MEqq2W2ffPowheg PowhegMEqq2W2ff
newdef PowhegMEqq2W2ff:Process 2
newdef PowhegMEqq2W2ff:Coupling /Herwig/Shower/AlphaQCD
# PP->ZH
create Herwig::MEPP2ZHPowheg PowhegMEPP2ZH
newdef PowhegMEPP2ZH:Coupling /Herwig/Shower/AlphaQCD
# PP->WH
create Herwig::MEPP2WHPowheg PowhegMEPP2WH
newdef PowhegMEPP2WH:Coupling /Herwig/Shower/AlphaQCD
# hadron-hadron to higgs
create Herwig::MEPP2HiggsPowheg PowhegMEHiggs
newdef PowhegMEHiggs:ShapeScheme MassGenerator
newdef PowhegMEHiggs:Process gg
newdef PowhegMEHiggs:Coupling /Herwig/Shower/AlphaQCD
# PP->VV
create Herwig::MEPP2VVPowheg PowhegMEPP2VV
newdef PowhegMEPP2VV:Coupling /Herwig/Shower/AlphaQCD
# PP -> Higgs via VBF
create Herwig::MEPP2HiggsVBFPowheg PowhegMEPP2HiggsVBF
newdef PowhegMEPP2HiggsVBF:ShowerAlphaQCD /Herwig/Shower/AlphaQCD
# PP -> diphoton NLO
create Herwig::MEPP2GammaGammaPowheg MEGammaGammaPowheg
set MEGammaGammaPowheg:Process 0
set MEGammaGammaPowheg:Contribution 1
set MEGammaGammaPowheg:ShowerAlphaQCD /Herwig/Shower/AlphaQCD
set MEGammaGammaPowheg:ShowerAlphaQED /Herwig/Shower/AlphaQED
##########################################################
# DIS matrix elements
##########################################################
# neutral current
create Herwig::MENeutralCurrentDIS MEDISNC
newdef MEDISNC:Coupling /Herwig/Shower/AlphaQCD
newdef MEDISNC:Contribution 0
# charged current
create Herwig::MEChargedCurrentDIS MEDISCC
newdef MEDISCC:Coupling /Herwig/Shower/AlphaQCD
newdef MEDISCC:Contribution 0
# neutral current (POWHEG)
create Herwig::MENeutralCurrentDIS PowhegMEDISNC
newdef PowhegMEDISNC:Coupling /Herwig/Shower/AlphaQCD
newdef PowhegMEDISNC:Contribution 1
# charged current (POWHEG)
create Herwig::MEChargedCurrentDIS PowhegMEDISCC
newdef PowhegMEDISCC:Coupling /Herwig/Shower/AlphaQCD
newdef PowhegMEDISCC:Contribution 1
##########################################################
# Gamma-Gamma matrix elements
##########################################################
# fermion-antiferimon
create Herwig::MEGammaGamma2ff MEgg2ff HwMEGammaGamma.so
# W+ W-
create Herwig::MEGammaGamma2WW MEgg2WW HwMEGammaGamma.so
##########################################################
# Gamma-Hadron matrix elements
##########################################################
# gamma parton -> 2 jets
create Herwig::MEGammaP2Jets MEGammaP2Jets HwMEGammaHadron.so
##########################################################
# Set up the Subprocesses
#
# Generic for all colliders
##########################################################
create ThePEG::SubProcessHandler SubProcess
newdef SubProcess:PartonExtractor /Herwig/Partons/PPExtractor

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:37 PM (16 h, 16 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486661
Default Alt Text
(46 KB)

Event Timeline