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 +#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 +MEee2ff::diagrams(const DiagramVector & diags) const { + double lastCont(0.5),lastBW(0.5); + if ( lastXCombPtr() ) { + lastCont = meInfo()[0]; + lastBW = meInfo()[1]; + } + Selector 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 +MEee2ff::colourGeometries(tcDiagPtr) const { + static ColourLines ctST(" "); + Selector 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 +describeMEee2ff("Herwig::MEee2ff", "HwMELepton.so"); + +void MEee2ff::Init() { + + static ClassDocumentation documentation + ("The MEee2ff class implements the matrix element for" + "e+e- to leptons via Z and photon exchange using helicity amplitude" + "techniques"); + + static Switch 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 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 interfacePrefactor + ("Prefactor", + "Prefactor for the overestimate of the emission probability", + &MEee2ff::preFactor_, 6.0, 1.0, 100.0, + false, false, Interface::limited); + + static Reference interfaceEMCoupling + ("AlphaQED", + "Pointer to the object to calculate the EM coupling for the correction", + &MEee2ff::alphaQED_, false, false, true, false, false); + + static Switch 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 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 & partons, + const vector & momenta, + bool first) const { + vector fin,aout; + vector 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 & fin, + vector & ain, + vector & fout, + vector & 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(mePartonData()[0]), + dynamic_ptr_cast(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()id()) swap(hard[0],hard[1]); + if(hard[2]->id()id()) swap(hard[2],hard[3]); + vector fin,aout; + vector 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(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(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 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;ixbornIncoming().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 partons, + vector 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 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 & partons, + const vector & momenta) const { + // compute the spinors + vector fin,aout; + vector ain,fout; + vector 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 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(partons[0]), + dynamic_ptr_cast(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 & 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 realMomenta[2][2]= + {{vector(5),vector(5)}, + {vector(5),vector(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] 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]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]pT[1]) { + iemit = 2; + ispect = 3; + pTemit = pT[0]; + if(UseRandom::rnd()(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 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 + 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 & partons, + const vector & 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 & fin, + vector & ain, + vector & fout, + vector & 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 partons, + vector 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 & partons, + const vector & momenta) const; + + /** + * Generate the momenta for a hard configuration + */ + Energy generateHard(RealEmissionProcessPtr tree, + vector & 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 partons_; + + /** + * Momenta of the leading-order partons + */ + vector 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