Page MenuHomeHEPForge
Files F8725711

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Hadron/MEPP2HiggsVBF.cc b/MatrixElement/Hadron/MEPP2HiggsVBF.cc
--- a/MatrixElement/Hadron/MEPP2HiggsVBF.cc
+++ b/MatrixElement/Hadron/MEPP2HiggsVBF.cc
@@ -1,144 +1,1617 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEPP2HiggsVBF class.
//
#include "MEPP2HiggsVBF.h"
#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
+#include "ThePEG/PDT/StandardMatchers.h"
+#include <numeric>
+#include "Herwig++/Utilities/Maths.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "ThePEG/Repository/CurrentGenerator.h"
using namespace Herwig;
-MEPP2HiggsVBF::MEPP2HiggsVBF()
+namespace {
+using namespace Herwig;
+using namespace ThePEG;
+using namespace ThePEG::Helicity;
+
+void debuggingMatrixElement(bool BGF,
+ tcPDPtr partons1, tcPDPtr partons2,
+ tcPDPtr partons3, tcPDPtr partons4,
+ const Lorentz5Momentum & psystem0,
+ const Lorentz5Momentum & psystem1,
+ const Lorentz5Momentum & pother0,
+ const Lorentz5Momentum & pother1,
+ const Lorentz5Momentum & p0,
+ const Lorentz5Momentum & p1,
+ const Lorentz5Momentum & p2,
+ const Lorentz5Momentum & phiggs,
+ Energy2 Q12, Energy2 scale,
+ double old) {
+ // get the vertex and the boson
+ tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>
+ (CurrentGenerator::current().standardModel());
+ assert(hwsm);
+ tcPDPtr boson;
+ AbstractFFVVertexPtr weakVertex;
+ AbstractFFVVertexPtr strongVertex = hwsm->vertexFFG();
+ AbstractVVSVertexPtr higgsVertex = hwsm->vertexWWH();
+ if(partons1->id()==partons2->id()) {
+ weakVertex = hwsm->vertexFFZ();
+ boson = hwsm->getParticleData(ParticleID::Z0);
+ }
+ else {
+ weakVertex = hwsm->vertexFFW();
+ boson = hwsm->getParticleData(ParticleID::Wplus);
+ }
+ tcPDPtr hdata = hwsm->getParticleData(ParticleID::h0);
+ tcPDPtr gluon = hwsm->getParticleData(ParticleID::g);
+ SpinorWaveFunction q1,q2;
+ SpinorBarWaveFunction qbar1,qbar2;
+ if(partons1->id()>0) {
+ q1 = SpinorWaveFunction (psystem0,partons1,incoming);
+ qbar1 = SpinorBarWaveFunction(psystem1,partons2,outgoing);
+ }
+ else {
+ q1 = SpinorWaveFunction (psystem1,partons2,outgoing);
+ qbar1 = SpinorBarWaveFunction(psystem0,partons1,incoming);
+ }
+ if(partons3->id()>0) {
+ q2 = SpinorWaveFunction (pother0,partons3,incoming);
+ qbar2 = SpinorBarWaveFunction(pother1,partons4,outgoing);
+ }
+ else {
+ q2 = SpinorWaveFunction (pother1,partons4,outgoing);
+ qbar2 = SpinorBarWaveFunction(pother0,partons3,incoming);
+ }
+ ScalarWaveFunction higgs(phiggs,hdata,outgoing);
+ if(!BGF) {
+ SpinorWaveFunction q1p;
+ SpinorBarWaveFunction qbar1p;
+ if(partons1->id()>0) {
+ q1p = SpinorWaveFunction (p0 ,partons1,incoming);
+ qbar1p = SpinorBarWaveFunction(p1 ,partons2,outgoing);
+ }
+ else {
+ q1p = SpinorWaveFunction (p1 ,partons2,outgoing);
+ qbar1p = SpinorBarWaveFunction(p0 ,partons1,incoming);
+ }
+ VectorWaveFunction gl(p2,gluon,outgoing);
+ double lome(0.),realme(0.);
+ for(unsigned int lhel1=0;lhel1<2;++lhel1) {
+ q2.reset(lhel1);
+ for(unsigned int lhel2=0;lhel2<2;++lhel2) {
+ qbar2.reset(lhel2);
+ VectorWaveFunction off1
+ = weakVertex->evaluate(scale,3,boson,q2,qbar2);
+ VectorWaveFunction off2
+ = higgsVertex->evaluate(scale,3,boson,off1,higgs);
+ for(unsigned int qhel1=0;qhel1<2;++qhel1) {
+ q1.reset(qhel1);
+ q1p.reset(qhel1);
+ for(unsigned int qhel2=0;qhel2<2;++qhel2) {
+ qbar1.reset(qhel2);
+ qbar1p.reset(qhel2);
+ Complex diag = weakVertex->evaluate(scale,q1,qbar1,off2);
+ lome += norm(diag);
+ for(unsigned int ghel=0;ghel<2;++ghel) {
+ gl.reset(2*ghel);
+ SpinorWaveFunction inter1 =
+ strongVertex->evaluate(Q12,5,q1p.particle(),q1p,gl);
+ Complex diag1 = weakVertex->evaluate(scale,inter1,qbar1p,off2);
+ SpinorBarWaveFunction inter2 =
+ strongVertex->evaluate(Q12,5,qbar1p.particle(),qbar1p,gl);
+ Complex diag2 = weakVertex->evaluate(scale,q1p,inter2,off2);
+ realme += norm(diag1+diag2);
+ }
+ }
+ }
+ }
+ }
+ double test1 = realme/lome/hwsm->alphaS(Q12)*Q12*UnitRemoval::InvE2;
+ cerr << "testing ratio A " << old/test1 << "\n";
+ }
+ else {
+ SpinorWaveFunction q1p;
+ SpinorBarWaveFunction qbar1p;
+ if(partons1->id()>0) {
+ q1p = SpinorWaveFunction (p2,partons1->CC(),outgoing);
+ qbar1p = SpinorBarWaveFunction(p1,partons2 ,outgoing);
+ }
+ else {
+ q1p = SpinorWaveFunction (p1,partons2 ,outgoing);
+ qbar1p = SpinorBarWaveFunction(p2,partons1->CC(),outgoing);
+ }
+ VectorWaveFunction gl(p0,gluon,incoming);
+ double lome(0.),realme(0.);
+ for(unsigned int lhel1=0;lhel1<2;++lhel1) {
+ q2.reset(lhel1);
+ for(unsigned int lhel2=0;lhel2<2;++lhel2) {
+ qbar2.reset(lhel2);
+ VectorWaveFunction off1
+ = weakVertex->evaluate(scale,3,boson,q2,qbar2);
+ VectorWaveFunction off2
+ = higgsVertex->evaluate(scale,3,boson,off1,higgs);
+ for(unsigned int qhel1=0;qhel1<2;++qhel1) {
+ q1.reset(qhel1);
+ q1p.reset(qhel1);
+ for(unsigned int qhel2=0;qhel2<2;++qhel2) {
+ qbar1.reset(qhel2);
+ qbar1p.reset(qhel2);
+ Complex diag = weakVertex->evaluate(scale,q1,qbar1,off2);
+ lome += norm(diag);
+ for(unsigned int ghel=0;ghel<2;++ghel) {
+ gl.reset(2*ghel);
+ SpinorWaveFunction inter1 =
+ strongVertex->evaluate(Q12,5,q1p.particle(),q1p,gl);
+ Complex diag1 = weakVertex->evaluate(scale,inter1,qbar1p,off2);
+ SpinorBarWaveFunction inter2 =
+ strongVertex->evaluate(Q12,5,qbar1p.particle(),qbar1p,gl);
+ Complex diag2 = weakVertex->evaluate(scale,q1p,inter2,off2);
+ realme += norm(diag1+diag2);
+ }
+ }
+ }
+ }
+ }
+ double test1 = realme/lome/hwsm->alphaS(Q12)*Q12*UnitRemoval::InvE2;
+ cerr << "testing ratio B " << old/test1 << "\n";
+ }
+}
+}
+
+MEPP2HiggsVBF::MEPP2HiggsVBF() : comptonWeight_(8.), BGFWeight_(30.),
+ pTmin_(1.*GeV),initial_(10.),final_(8.),
+ procProb_(0.5), comptonInt_(0.), bgfInt_(0.),
+ nover_(0),maxwgt_(make_pair(0.,0.))
{}
+void MEPP2HiggsVBF::doinit() {
+ gluon_ = getParticleData(ParticleID::g);
+ // integrals of me over phase space
+ double r5=sqrt(5.),darg((r5-1.)/(r5+1.)),ath(0.5*log((1.+1./r5)/(1.-1./r5)));
+ comptonInt_ = 2.*(-21./20.-6./(5.*r5)*ath+sqr(Constants::pi)/3.
+ -2.*Math::ReLi2(1.-darg)-2.*Math::ReLi2(1.-1./darg));
+ bgfInt_ = 121./9.-56./r5*ath;
+ // get the vertex pointers from the SM object
+ tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
+ if(!hwsm)
+ throw InitException() << "Wrong type of StandardModel object in "
+ << "MEPP2HiggsVBF::doinit() the Herwig++"
+ << " version must be used"
+ << Exception::runerror;
+ // set the vertex
+ setWWHVertex(hwsm->vertexWWH());
+ higgs(getParticleData(ParticleID::h0));
+ MEfftoffH::doinit();
+}
+
+void MEPP2HiggsVBF::dofinish() {
+ MEfftoffH::dofinish();
+ if(nover_==0) return;
+ generator()->log() << "VBFMECorrection when applying the hard correction "
+ << nover_ << " weights larger than one were generated of which"
+ << " the largest was " << maxwgt_.first << " for the QCD compton"
+ << " processes and " << maxwgt_.second << " for the BGF process\n";
+}
+
void MEPP2HiggsVBF::getDiagrams() const {
// get the quark particle data objects as we'll be using them
tcPDPtr q[6],qbar[6];
for ( int ix=0; ix<5; ++ix ) {
q [ix] = getParticleData(ix+1);
qbar[ix] = q[ix]->CC();
}
// WW processes
if(process()==0||process()==1) {
std::vector<pair<tcPDPtr,tcPDPtr> > parentpair;
parentpair.reserve(6);
// don't even think of putting 'break' in here!
switch(maxFlavour()) {
case 5:
if (minFlavour()<=4)
parentpair.push_back(make_pair(getParticleData(ParticleID::b),
getParticleData(ParticleID::c)));
if (minFlavour()<=2)
parentpair.push_back(make_pair(getParticleData(ParticleID::b),
getParticleData(ParticleID::u)));
case 4:
if (minFlavour()<=3)
parentpair.push_back(make_pair(getParticleData(ParticleID::s),
getParticleData(ParticleID::c)));
if (minFlavour()<=1)
parentpair.push_back(make_pair(getParticleData(ParticleID::d),
getParticleData(ParticleID::c)));
case 3:
if (minFlavour()<=2)
parentpair.push_back(make_pair(getParticleData(ParticleID::s),
getParticleData(ParticleID::u)));
case 2:
if (minFlavour()<=1)
parentpair.push_back(make_pair(getParticleData(ParticleID::d),
getParticleData(ParticleID::u)));
default:
;
}
for(unsigned int ix=0;ix<parentpair.size();++ix) {
for(unsigned int iy=0;iy<parentpair.size();++iy) {
// q1 q2 -> q1' q2' h
if(parentpair[ix].first->id()<parentpair[iy].second->id()) {
add(new_ptr((Tree2toNDiagram(4), parentpair[ix].first, WMinus(), WPlus(),
parentpair[iy].second, 1, parentpair[ix].second, 4,
parentpair[iy].first, 2, higgs(),-1)));
}
else {
add(new_ptr((Tree2toNDiagram(4), parentpair[iy].second, WPlus(), WMinus(),
parentpair[ix].first, 1, parentpair[iy].first, 4,
parentpair[ix].second, 2, higgs(),-1)));
}
// q1 qbar2 -> q1' qbar2' h
add(new_ptr((Tree2toNDiagram(4), parentpair[ix].first, WMinus(), WPlus(),
parentpair[iy].first->CC(), 1,
parentpair[ix].second, 4, parentpair[iy].second->CC(),
2, higgs(),-1)));
add(new_ptr((Tree2toNDiagram(4),parentpair[iy].second, WPlus(), WMinus(),
parentpair[ix].second->CC(), 1, parentpair[iy].first,
4, parentpair[ix].first->CC(),
2, higgs(),-1)));
// qbar1 qbar2 -> qbar1' qbar2' h
if(parentpair[ix].first->id()<parentpair[ix].second->id()) {
add(new_ptr((Tree2toNDiagram(4), parentpair[ix].first->CC(), WPlus(), WMinus(),
parentpair[iy].second->CC(), 1,
parentpair[ix].second->CC(), 4, parentpair[iy].first->CC(),
2, higgs(),-1)));
}
else {
add(new_ptr((Tree2toNDiagram(4), parentpair[iy].second->CC(), WMinus(), WPlus(),
parentpair[ix].first->CC(), 1,
parentpair[iy].first->CC(), 4, parentpair[ix].second->CC(),
2, higgs(),-1)));
}
}
}
}
// ZZ processes
if(process()==0||process()==2) {
for(unsigned int ix=minFlavour()-1;ix<maxFlavour();++ix) {
for(unsigned int iy=ix;iy<maxFlavour();++iy) {
// q q -> q q H
add(new_ptr((Tree2toNDiagram(4), q[ix], Z0(), Z0(), q[iy],
1, q[ix], 4, q[iy], 2, higgs(),-2)));
// qbar qbar -> qbar qbar H
add(new_ptr((Tree2toNDiagram(4), qbar[ix], Z0(), Z0(), qbar[iy],
1, qbar[ix], 4, qbar[iy], 2, higgs(),-2)));
}
// q qbar -> q qbar H
for(unsigned int iy=minFlavour()-1;iy<maxFlavour();++iy) {
add(new_ptr((Tree2toNDiagram(4), q[ix], Z0(), Z0(), qbar[iy],
1, q[ix], 4, qbar[iy], 2, higgs(),-2)));
}
}
}
}
-void MEPP2HiggsVBF::persistentOutput(PersistentOStream & ) const {
+void MEPP2HiggsVBF::persistentOutput(PersistentOStream & os) const {
+ os << initial_ << final_
+ << alpha_ << ounit(pTmin_,GeV) << comptonWeight_ << BGFWeight_ << gluon_
+ << comptonInt_ << bgfInt_ << procProb_;
}
-void MEPP2HiggsVBF::persistentInput(PersistentIStream & , int) {
+void MEPP2HiggsVBF::persistentInput(PersistentIStream & is, int) {
+ is >> initial_ >> final_
+ >> alpha_ >> iunit(pTmin_,GeV) >> comptonWeight_ >> BGFWeight_ >> gluon_
+ >> comptonInt_ >> bgfInt_ >> procProb_;
}
ClassDescription<MEPP2HiggsVBF> MEPP2HiggsVBF::initMEPP2HiggsVBF;
// Definition of the static class description member.
void MEPP2HiggsVBF::Init() {
static ClassDocumentation<MEPP2HiggsVBF> documentation
("The MEPP2HiggsVBF class implements Higgs production via vector-boson fusion");
+ static Reference<MEPP2HiggsVBF,ShowerAlpha> interfaceShowerAlphaQCD
+ ("ShowerAlphaQCD",
+ "The object calculating the strong coupling constant",
+ &MEPP2HiggsVBF::alpha_, false, false, true, false, false);
+
+ static Parameter<MEPP2HiggsVBF,Energy> interfacepTMin
+ ("pTMin",
+ "The minimum pT",
+ &MEPP2HiggsVBF::pTmin_, GeV, 1.*GeV, 0.0*GeV, 10.0*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<MEPP2HiggsVBF,double> interfaceComptonWeight
+ ("ComptonWeight",
+ "Weight for the overestimate ofthe compton channel",
+ &MEPP2HiggsVBF::comptonWeight_, 50.0, 0.0, 100.0,
+ false, false, Interface::limited);
+
+ static Parameter<MEPP2HiggsVBF,double> interfaceBGFWeight
+ ("BGFWeight",
+ "Weight for the overestimate of the BGF channel",
+ &MEPP2HiggsVBF::BGFWeight_, 100.0, 0.0, 1000.0,
+ false, false, Interface::limited);
+
+ static Parameter<MEPP2HiggsVBF,double> interfaceProcessProbability
+ ("ProcessProbability",
+ "The probabilty of the QCD compton process for the process selection",
+ &MEPP2HiggsVBF::procProb_, 0.3, 0.0, 1.,
+ false, false, Interface::limited);
+
}
-void MEPP2HiggsVBF::doinit() {
- // get the vedrtex pointers from the SM object
- tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
- if(!hwsm)
- throw InitException() << "Wrong type of StandardModel object in "
- << "MEPP2HiggsVBF::doinit() the Herwig++"
- << " version must be used"
- << Exception::runerror;
- // set the vertex
- setWWHVertex(hwsm->vertexWWH());
- higgs(getParticleData(ParticleID::h0));
- MEfftoffH::doinit();
+HardTreePtr MEPP2HiggsVBF::generateHardest(ShowerTreePtr tree) {
+ pair< tShowerParticlePtr, tShowerParticlePtr> first,second;
+ pair<tcBeamPtr,tcBeamPtr> beams;
+ pair<tPPtr,tPPtr> hadrons;
+ // get the incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(!first.first) {
+ first.first = cit->first->progenitor();
+ beams.first = cit->first->beam();
+ hadrons.first = cit->first->original()->parents()[0];
+ }
+ else {
+ second.first = cit->first->progenitor();
+ beams.second = cit->first->beam();
+ hadrons.second = cit->first->original()->parents()[0];
+ }
+ }
+ // and the outgoing
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cjt=tree->outgoingLines().begin();cjt!=tree->outgoingLines().end();++cjt) {
+ if(cjt->first->progenitor()->id()==ParticleID::h0) {
+ higgs_ = cjt->first->progenitor();
+ }
+ else {
+ if(abs(cjt->first->progenitor()->id())>5) continue;
+ if(cjt->first->progenitor()->colourLine()&&
+ cjt->first->progenitor()->colourLine()==first.first->colourLine()) {
+ first.second = cjt->first->progenitor();
+ continue;
+ }
+ if(cjt->first->progenitor()->antiColourLine()&&
+ cjt->first->progenitor()->antiColourLine()==first.first->antiColourLine()) {
+ first.second = cjt->first->progenitor();
+ continue;
+ }
+ if(cjt->first->progenitor()->colourLine()&&
+ cjt->first->progenitor()->colourLine()==second.first->colourLine()) {
+ second.second = cjt->first->progenitor();
+ continue;
+ }
+ if(cjt->first->progenitor()->antiColourLine()&&
+ cjt->first->progenitor()->antiColourLine()==second.first->antiColourLine()) {
+ second.second = cjt->first->progenitor();
+ continue;
+ }
+ }
+ }
+ // loop over the two possible emitting systems
+ q_ [0] = first .second->momentum()-first .first->momentum();
+ q2_[0] = -q_[0].m2();
+ q_ [1] = second.second->momentum()-second.first->momentum();
+ q2_[1] = -q_[1].m2();
+ for(unsigned int ix=0;ix<2;++ix) {
+ if(ix==1) {
+ swap(first,second);
+ swap(beams.first,beams.second);
+ }
+ // check beam, all particles
+ assert(beams.first && higgs_ &&
+ first .first && first.second &&
+ second.first && second.second);
+ // beam and pdf
+ beam_[ix] = beams.first;
+ pdf_ [ix] = beam_[ix]->pdf();
+ assert(beam_[ix] && pdf_[ix] );
+ // Particle data objects
+ partons_[ix][0] = first. first->dataPtr();
+ partons_[ix][1] = first.second->dataPtr();
+ partons_[ix][2] = second. first->dataPtr();
+ partons_[ix][3] = second.second->dataPtr();
+ // extract the born variables
+ xB_[ix] = first.first->x();
+ Lorentz5Momentum pb = first.first->momentum();
+ Axis axis(q_[ix].vect().unit());
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ rot_[ix] = LorentzRotation();
+ if(axis.perp2()>1e-20) {
+ rot_[ix].setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ rot_[ix].rotateX(Constants::pi);
+ }
+ if(abs(1.-q_[ix].e()/q_[ix].vect().mag())>1e-6)
+ rot_[ix].boostZ( q_[ix].e()/q_[ix].vect().mag());
+ pb *= rot_[ix];
+ if(pb.perp2()/GeV2>1e-20) {
+ Boost trans = -1./pb.e()*pb.vect();
+ trans.setZ(0.);
+ rot_[ix].boost(trans);
+ }
+ // momenta of the particles
+ phiggs_ [ix] = rot_[ix]*higgs_->momentum();
+ pother_ [ix][0] = rot_[ix]*second. first->momentum();
+ pother_ [ix][1] = rot_[ix]*second.second->momentum();
+ psystem_[ix][0] = rot_[ix]* first. first->momentum();
+ psystem_[ix][1] = rot_[ix]* first.second->momentum();
+ q_[ix] *= rot_[ix];
+ pTCompton_[ix] = pTBGF_[ix] = ZERO;
+ // generate a compton point
+ generateCompton(ix);
+ // generate a BGF point
+ generateBGF(ix);
+ }
+ // no valid emissions, return
+ if(pTCompton_[0]<ZERO && pTCompton_[1]<ZERO&&
+ pTBGF_ [0]<ZERO && pTBGF_ [1]<ZERO) {
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pTmin_);
+ }
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pTmin_);
+ }
+ return HardTreePtr();
+ }
+ // find the maximum pT emission
+ unsigned int system = 0;
+ bool isCompton = false;
+ Energy pTmax = -GeV;
+ for(unsigned int ix=0;ix<2;++ix) {
+ if(pTCompton_[ix]>pTmax) {
+ pTmax = pTCompton_[ix];
+ isCompton = true;
+ system = ix;
+ }
+ if(pTBGF_[ix]>pTmax) {
+ pTmax = pTBGF_[ix];
+ isCompton = false;
+ system = ix;
+ }
+ }
+ if(system==0) {
+ swap(first,second);
+ swap(beams.first,beams.second);
+ }
+ // add the non emitting particles
+ vector<HardBranchingPtr> spaceBranchings,allBranchings;
+ spaceBranchings.push_back(new_ptr(HardBranching(second.first,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ allBranchings.push_back(spaceBranchings.back());
+ allBranchings.push_back(new_ptr(HardBranching(second.second,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ allBranchings.push_back(new_ptr(HardBranching(higgs_,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ rot_[system].invert();
+ // compton hardest
+ if(isCompton) {
+ for(unsigned int ix=0;ix<ComptonMomenta_[system].size();++ix) {
+ ComptonMomenta_[system][ix].transform(rot_[system]);
+ }
+ ShowerParticlePtr newqout (new_ptr(ShowerParticle(partons_[system][1],true)));
+ newqout->set5Momentum(ComptonMomenta_[system][1]);
+ ShowerParticlePtr newg(new_ptr(ShowerParticle(gluon_,true)));
+ newg->set5Momentum(ComptonMomenta_[system][2]);
+ ShowerParticlePtr newqin (new_ptr(ShowerParticle(partons_[system][0],false )));
+ newqin->set5Momentum(ComptonMomenta_[system][0]);
+ if(ComptonISFS_[system]) {
+ ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[system][0],false)));
+ newspace->set5Momentum(ComptonMomenta_[system][0]-ComptonMomenta_[system][2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),
+ spaceBranch,
+ HardBranching::Incoming)));
+ spaceBranch->addChild(offBranch);
+ HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),spaceBranch,
+ HardBranching::Outgoing)));
+ spaceBranch->addChild(g);
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ allBranchings.push_back(outBranch);
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newqout ,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ else {
+ ShowerParticlePtr newtime(new_ptr(ShowerParticle(partons_[system][1],true)));
+ newtime->set5Momentum(ComptonMomenta_[system][1]+ComptonMomenta_[system][2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newtime,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),offBranch,
+ HardBranching::Outgoing)));
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),offBranch,
+ HardBranching::Outgoing)));
+ offBranch->addChild(outBranch);
+ offBranch->addChild(g);
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newqin ,newqin->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newtime,newqin->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ }
+ // BGF hardest
+ else {
+ for(unsigned int ix=0;ix<BGFMomenta_[system].size();++ix) {
+ BGFMomenta_[system][ix].transform(rot_[system]);
+ }
+ ShowerParticlePtr newq (new_ptr(ShowerParticle(partons_[system][1],true)));
+ newq->set5Momentum(BGFMomenta_[system][1]);
+ ShowerParticlePtr newqbar(new_ptr(ShowerParticle(partons_[system][0]->CC(),true)));
+ newqbar->set5Momentum(BGFMomenta_[system][2]);
+ ShowerParticlePtr newg (new_ptr(ShowerParticle(gluon_,false)));
+ newg->set5Momentum(BGFMomenta_[system][0]);
+ ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[system][0],false)));
+ newspace->set5Momentum(BGFMomenta_[system][0]-BGFMomenta_[system][2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newg,SudakovPtr(),HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),spaceBranch,
+ HardBranching::Incoming)));
+ HardBranchingPtr qbar(new_ptr(HardBranching(newqbar,SudakovPtr(),spaceBranch,
+ HardBranching::Outgoing)));
+ spaceBranch->addChild(offBranch);
+ spaceBranch->addChild(qbar);
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newq,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ allBranchings.push_back(outBranch);
+
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newq ,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ HardTreePtr newTree(new_ptr(HardTree(allBranchings,spaceBranchings,
+ ShowerInteraction::QCD)));
+ // Set the maximum pt for all other emissions and connect hard and shower tree
+ Energy pT = isCompton ? pTCompton_[system] : pTBGF_[system];
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ // set maximum pT
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pT);
+ set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+ if(cit->first->progenitor()==first.first) {
+ ++cjt;++cjt;++cjt;
+ }
+ newTree->connect(cit->first->progenitor(),*cjt);
+ tPPtr beam =cit->first->original();
+ if(!beam->parents().empty()) beam=beam->parents()[0];
+ (*cjt)->beam(beam);
+ HardBranchingPtr parent=(*cjt)->parent();
+ while(parent) {
+ parent->beam(beam);
+ parent=parent->parent();
+ };
+ }
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ // set maximum pT
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pT);
+ for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+ cjt!=newTree->branchings().end();++cjt) {
+ if((*cjt)->branchingParticle()->isFinalState()&&
+ (*cjt)->branchingParticle()->id()==cit->first->progenitor()->id()) {
+ newTree->connect(cit->first->progenitor(),*cjt);
+ }
+ }
+ }
+ // set the evolution partners and scales
+ ShowerParticleVector particles;
+ for(set<HardBranchingPtr>::iterator cit=newTree->branchings().begin();
+ cit!=newTree->branchings().end();++cit) {
+ particles.push_back((*cit)->branchingParticle());
+ }
+ for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+ cjt!=newTree->branchings().end();++cjt) {
+ if(cjt==newTree->branchings().begin()) {
+ (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
+ ++cjt;
+ (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
+ ++cjt;
+ (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
+ ++cjt;
+ }
+ }
+ return newTree;
}
+
+void MEPP2HiggsVBF::generateCompton(unsigned int system) {
+ // calculate the A coefficient for the correlations
+ acoeff_ = A(partons_[system][0],partons_[system][1],
+ partons_[system][2],partons_[system][3]);
+ // maximum value of the xT
+ double xT = sqrt((1.-xB_[system])/xB_[system]);
+ double xTMin = 2.*pTmin_/sqrt(q2_[system]);
+ double zp;
+ // prefactor
+ double a = alpha_->overestimateValue()*comptonWeight_/Constants::twopi;
+ // loop to generate kinematics
+ double wgt(0.),xp(0.);
+ l_ = 2.*pother_[system][0]/sqrt(q2_[system]);
+ m_ = 2.*pother_[system][1]/sqrt(q2_[system]);
+ vector<double> azicoeff;
+ do {
+ wgt = 0.;
+ // intergration variables dxT/xT^3
+ xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
+ // dz
+ zp = UseRandom::rnd();
+ xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
+ // check allowed
+ if(xp<xB_[system]||xp>1.) continue;
+ // phase-space piece of the weight
+ wgt = 8.*(1.-xp)*zp/comptonWeight_;
+ // PDF piece of the weight
+ Energy2 mu2 = q2_[system]*((1.-xp)*(1-zp)*zp/xp+1.);
+ double pdf = pdf_[system]->xfx(beam_[system],partons_[system][0],
+ mu2 ,xB_[system]/xp)/
+ pdf_[system]->xfx(beam_[system],partons_[system][0],
+ scale(),xB_[system] );
+ wgt *= max(pdf,0.);
+ // me piece of the weight
+ // double me = comptonME(system,xT,xp,zp,phi);
+ double x2 = 1.-(1.-zp)/xp;
+ azicoeff = ComptonME(xp,x2,xT,l_,m_);
+ double me = 4./3.*alpha_->ratio(0.25*q2_[system]*sqr(xT))*
+ (azicoeff[0]+0.5*azicoeff[2]+0.5*azicoeff[4]);
+ wgt *= me;
+ if(wgt>1.||wgt<0.) {
+ ostringstream wstring;
+ wstring << "MEPP2HiggsVBF::generateCompton() "
+ << "Weight greater than one or less than zero"
+ << "wgt = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(xT>xTMin&&UseRandom::rnd()>wgt);
+ if(xT<=xTMin) {
+ pTCompton_[system]=-GeV;
+ return;
+ }
+ // generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi)),sphi(sin(phi));
+ phiwgt = azicoeff[0]+azicoeff[5]*sphi*cphi
+ +azicoeff[1]*cphi+azicoeff[2]*sqr(cphi)
+ +azicoeff[3]*sphi+azicoeff[4]*sqr(sphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in MEPP2HiggsVBF"
+ << "::generateCompton() to"
+ << " generate phi" << Exception::eventerror;
+ // momenta for the configuration
+ Energy Q(sqrt(q2_[system]));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ pTCompton_[system] = 0.5*Q*xT;
+ ComptonMomenta_[system].resize(3);
+ ComptonMomenta_[system][0] = p0;
+ ComptonMomenta_[system][1] = p1;
+ ComptonMomenta_[system][2] = p2;
+ ComptonISFS_[system] = zp>xp;
+}
+
+double MEPP2HiggsVBF::comptonME(unsigned int system, double xT,
+ double xp, double zp, double phi) {
+ // scale and prefactors
+ double CFfact = 4./3.*alpha_->ratio(0.25*q2_[system]*sqr(xT));
+ Energy Q(sqrt(q2_[system]));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ //set NLO momenta
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ Lorentz5Momentum qnlo = p2+p1-p0;
+ // Breit frame variables
+ Lorentz5Momentum r1 = -p0/x1;
+ Lorentz5Momentum r2 = p1/x2;
+ // electroweak parameters
+ Energy2 mz2 = sqr(Z0()->mass());
+ Energy2 mw2 = sqr(WPlus()->mass());
+ double c0L,c1L,c0R,c1R;
+ Energy2 mb2;
+ // W
+ if(partons_[system][0]->id()!=partons_[system][1]->id()) {
+ mb2 = mw2;
+ c0L = sqrt(0.5);
+ c0R = 0;
+ c1L = sqrt(0.5);
+ c1R = 0;
+ }
+ // Z
+ else {
+ mb2 = mz2;
+ if(abs(partons_[system][0]->id())%2==0) {
+ c0L =
+ generator()->standardModel()->vu()+
+ generator()->standardModel()->au();
+ c0R =
+ generator()->standardModel()->vu()-
+ generator()->standardModel()->au();
+ }
+ else {
+ c0L =
+ generator()->standardModel()->vd()+
+ generator()->standardModel()->ad();
+ c0R =
+ generator()->standardModel()->vd()-
+ generator()->standardModel()->ad();
+ }
+ if(abs(partons_[system][2]->id())%2==0) {
+ c1L =
+ generator()->standardModel()->vu()+
+ generator()->standardModel()->au();
+ c1R =
+ generator()->standardModel()->vu()-
+ generator()->standardModel()->au();
+ }
+ else {
+ c1L =
+ generator()->standardModel()->vd()+
+ generator()->standardModel()->ad();
+ c1R =
+ generator()->standardModel()->vd()-
+ generator()->standardModel()->ad();
+ }
+ c0L *= 0.25;
+ c0R *= 0.25;
+ c1L *= 0.25;
+ c1R *= 0.25;
+ }
+ // Matrix element variables
+ double G1 = sqr(c0L*c1L)+sqr(c0R*c1R);
+ double G2 = sqr(c0L*c1R)+sqr(c0R*c1L);
+ Energy4 term1,term2,term3,loME;
+ if(partons_[system][0]->id()>0) {
+ if(partons_[system][2]->id()>0) {
+ term1 = loMatrixElement(r1 ,pother_[system][0],
+ qnlo+r1 ,pother_[system][1],G1,G2);
+ term2 = loMatrixElement(r2-qnlo ,pother_[system][0],
+ r2 ,pother_[system][1],G1,G2);
+ loME = loMatrixElement(psystem_[system][0],pother_[system][0],
+ psystem_[system][1],pother_[system][1],G1,G2);
+ }
+ else {
+ term1 = loMatrixElement(r1 ,pother_[system][1],
+ qnlo+r1 ,pother_[system][0],G1,G2);
+ term2 = loMatrixElement(r2-qnlo ,pother_[system][1],
+ r2 ,pother_[system][0],G1,G2);
+ loME = loMatrixElement(psystem_[system][0],pother_[system][1],
+ psystem_[system][1],pother_[system][0],G1,G2);
+ }
+ }
+ else {
+ if(partons_[system][2]->id()>0) {
+ term1 = loMatrixElement(qnlo+r1 ,pother_[system][0],
+ r1 ,pother_[system][1],G1,G2);
+ term2 = loMatrixElement(r2 ,pother_[system][0],
+ r2-qnlo ,pother_[system][1],G1,G2);
+ loME = loMatrixElement(psystem_[system][1],pother_[system][0],
+ psystem_[system][0],pother_[system][1],G1,G2);
+ }
+ else {
+ term1 = loMatrixElement(qnlo+r1,pother_[system][1],r1 ,
+ pother_[system][0],G1,G2);
+ term2 = loMatrixElement(r2 ,pother_[system][1],r2-qnlo,
+ pother_[system][0],G1,G2);
+ loME = loMatrixElement(psystem_[system][1],pother_[system][1],
+ psystem_[system][0],pother_[system][0],G1,G2);
+ }
+ }
+ double R1 = term1/loME;
+ double R2 = sqr(x2)/(sqr(x2)+sqr(xT))*(term2/loME);
+// debuggingMatrixElement(false,
+// partons_[system][0],partons_[system][1],
+// partons_[system][2],partons_[system][3],
+// psystem_[system][0],psystem_[system][1],
+// pother_ [system][0],pother_ [system][1],
+// p0,p1,p2,phiggs_[system],q2_[system],scale(),
+// 8.*Constants::pi/(1.-xp)/(1.-zp)*(R1+sqr(xp)*(sqr(x2)+sqr(xT))*R2));
+// cerr << "testing pieces A " << R1 << " " << sqr(xp)*(sqr(x2)+sqr(xT)) << " " << R2 << "\n";
+ return CFfact*(R1+sqr(xp)*(sqr(x2)+sqr(xT))*R2);
+}
+
+void MEPP2HiggsVBF::generateBGF(unsigned int system) {
+ // maximum value of the xT
+ double xT = (1.-xB_[system])/xB_[system];
+ double xTMin = 2.*pTmin_/sqrt(q2_[system]);
+ double zp;
+ // prefactor
+ double a = alpha_->overestimateValue()*BGFWeight_/Constants::twopi;
+ // loop to generate kinematics
+ double wgt(0.),xp(0.);
+ l_ = 2.*pother_[system][0]/sqrt(q2_[system]);
+ m_ = 2.*pother_[system][1]/sqrt(q2_[system]);
+ vector<double> azicoeff;
+ do {
+ wgt = 0.;
+ // intergration variables dxT/xT^3
+ xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
+ // dzp
+ zp = UseRandom::rnd();
+ xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
+ // check allowed
+ if(xp<xB_[system]||xp>1.) continue;
+ // phase-space piece of the weight
+ wgt = 8.*sqr(1.-xp)*zp/BGFWeight_;
+ // PDF piece of the weight
+ Energy2 mu2 = q2_[system]*((1.-xp)*(1-zp)*zp/xp+1.);
+ wgt *= pdf_[system]->xfx(beam_[system],gluon_ ,
+ mu2 ,xB_[system]/xp)/
+ pdf_[system]->xfx(beam_[system],partons_[system][0],
+ scale(),xB_[system]);
+ // me piece of the weight
+ //double me = BGFME(system,xT,xp,zp,phi);
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ azicoeff = BGFME(xp,x2,x3,xT,l_,m_);
+ double me = 0.5*alpha_->ratio(0.25*q2_[system]*sqr(xT))*
+ (azicoeff[0]+0.5*azicoeff[2]+0.5*azicoeff[4]);
+ wgt *= me;
+ if(wgt>1.||wgt<0.) {
+ ostringstream wstring;
+ wstring << "MEPP2HiggsVBF::generateBGF() "
+ << "Weight greater than one or less than zero"
+ << "wgt = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(xT>xTMin&&UseRandom::rnd()>wgt);
+ if(xT<=xTMin) {
+ pTBGF_[system] = -GeV;
+ return;
+ }
+ // generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi)),sphi(sin(phi));
+ phiwgt = azicoeff[0]+azicoeff[5]*sphi*cphi
+ +azicoeff[1]*cphi+azicoeff[2]*sqr(cphi)
+ +azicoeff[3]*sphi+azicoeff[4]*sqr(sphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in MEPP2HiggsVBF"
+ << "::generateBGF() to"
+ << " generate phi" << Exception::eventerror;
+ // momenta for the configuration
+ Energy Q(sqrt(q2_[system]));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ pTBGF_[system] = 0.5*Q*xT;
+ BGFMomenta_[system].resize(3);
+ BGFMomenta_[system][0] = p0;
+ BGFMomenta_[system][1] = p1;
+ BGFMomenta_[system][2] = p2;
+}
+
+double MEPP2HiggsVBF::BGFME(unsigned int system, double xT,
+ double xp, double zp, double phi) {
+ // scale and prefactors
+ double TRfact = 0.5*alpha_->ratio(0.25*q2_[system]*sqr(xT));
+ Energy Q(sqrt(q2_[system]));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ // Set NLO momenta
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ Lorentz5Momentum qnlo = p2+p1-p0;
+ // Breit frame variables
+ Lorentz5Momentum r2 = p1/x2;
+ Lorentz5Momentum r3 = -p2/x3;
+ // electroweak parameters
+ Energy2 mz2 = sqr(Z0()->mass());
+ Energy2 mw2 = sqr(WPlus()->mass());
+ double c0L,c1L,c0R,c1R;
+ Energy2 mb2;
+ // W
+ if(partons_[system][0]->id()!=partons_[system][1]->id()) {
+ mb2 = mw2;
+ c0L = sqrt(0.5);
+ c0R = 0;
+ c1L = sqrt(0.5);
+ c1R = 0;
+ }
+ // Z
+ else {
+ mb2 = mz2;
+ if(abs(partons_[system][0]->id())%2==0) {
+ c0L =
+ generator()->standardModel()->vu()+
+ generator()->standardModel()->au();
+ c0R =
+ generator()->standardModel()->vu()-
+ generator()->standardModel()->au();
+ }
+ else {
+ c0L =
+ generator()->standardModel()->vd()+
+ generator()->standardModel()->ad();
+ c0R =
+ generator()->standardModel()->vd()-
+ generator()->standardModel()->ad();
+ }
+ if(abs(partons_[system][2]->id())%2==0) {
+ c1L =
+ generator()->standardModel()->vu()+
+ generator()->standardModel()->au();
+ c1R =
+ generator()->standardModel()->vu()-
+ generator()->standardModel()->au();
+ }
+ else {
+ c1L =
+ generator()->standardModel()->vd()+
+ generator()->standardModel()->ad();
+ c1R =
+ generator()->standardModel()->vd()-
+ generator()->standardModel()->ad();
+ }
+ c0L *= 0.25;
+ c0R *= 0.25;
+ c1L *= 0.25;
+ c1R *= 0.25;
+ }
+ // Matrix element variables
+ double G1 = sqr(c0L*c1L)+sqr(c0R*c1R);
+ double G2 = sqr(c0L*c1R)+sqr(c0R*c1L);
+ Energy4 term1,term2,term3,loME;
+ if(partons_[system][0]->id()>0) {
+ if(partons_[system][2]->id()>0) {
+ term2 = loMatrixElement(r2-qnlo,pother_[system][0],
+ r2 ,pother_[system][1],G1,G2);
+ term3 = loMatrixElement(r3 ,pother_[system][0],
+ qnlo+r3,pother_[system][1],G1,G2);
+ loME = loMatrixElement(psystem_[system][0],pother_[system][0],
+ psystem_[system][1],pother_[system][1],G1,G2);
+ }
+ else {
+ term2 = loMatrixElement(r2-qnlo,pother_[system][1],
+ r2 ,pother_[system][0],G1,G2);
+ term3 = loMatrixElement(r3 ,pother_[system][1],
+ qnlo+r3,pother_[system][0],G1,G2);
+ loME = loMatrixElement(psystem_[system][0],pother_[system][1],
+ psystem_[system][1],pother_[system][0],G1,G2);
+ }
+ }
+ else {
+ if(partons_[system][2]->id()>0) {
+ term2 = loMatrixElement(r2 ,pother_[system][0],
+ r2-qnlo,pother_[system][1],G1,G2);
+ term3 = loMatrixElement(qnlo+r3,pother_[system][0],
+ r3 ,pother_[system][1],G1,G2);
+ loME = loMatrixElement(psystem_[system][1],pother_[system][0],
+ psystem_[system][0],pother_[system][1],G1,G2);
+ }
+ else {
+ term2 = loMatrixElement(r2 ,pother_[system][1],
+ r2-qnlo,pother_[system][0],G1,G2);
+ term3 = loMatrixElement(qnlo+r3,pother_[system][1],
+ r3 ,pother_[system][0],G1,G2);
+ loME = loMatrixElement(psystem_[system][1],pother_[system][1],
+ psystem_[system][0],pother_[system][0],G1,G2);
+ }
+ }
+ double R3 = sqr(x3)/(sqr(x3)+sqr(xT))*(term3/loME);
+ double R2 = sqr(x2)/(sqr(x2)+sqr(xT))*(term2/loME);
+// debuggingMatrixElement(true,
+// partons_[system][0],partons_[system][1],
+// partons_[system][2],partons_[system][3],
+// psystem_[system][0],psystem_[system][1],
+// pother_ [system][0],pother_ [system][1],
+// p0,p1,p2,phiggs_[system],q2_[system],
+// 8.*Constants::pi/zp/(1.-zp)*(sqr(xp)*(sqr(x3)+sqr(xT))*R3+
+// sqr(xp)*(sqr(x2)+sqr(xT))*R2));
+ return TRfact*
+ (sqr(xp)*(sqr(x3)+sqr(xT))*R3+
+ sqr(xp)*(sqr(x2)+sqr(xT))*R2);
+
+}
+
+Energy4 MEPP2HiggsVBF::loMatrixElement(const Lorentz5Momentum &p1,
+ const Lorentz5Momentum &p2,
+ const Lorentz5Momentum &q1,
+ const Lorentz5Momentum &q2,
+ double G1, double G2) const {
+ return G1*(p1*p2)*(q1*q2) + G2*(p1*q2)*(q1*p2);
+}
+
+void MEPP2HiggsVBF::initializeMECorrection(ShowerTreePtr tree, double & initial,
+ double & final) {
+ systems_.clear();
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ systems_.push_back(tChannelPair());
+ systems_.back().hadron = cit->first->original()->parents()[0];
+ systems_.back().beam = cit->first->beam();
+ systems_.back().incoming = cit->first->progenitor();
+ systems_.back().pdf = systems_.back().beam->pdf();
+ }
+ }
+ vector<ShowerParticlePtr> outgoing;
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cjt=tree->outgoingLines().begin();cjt!=tree->outgoingLines().end();++cjt) {
+ if(cjt->first->progenitor()->id()==ParticleID::h0)
+ higgs_ = cjt->first->progenitor();
+ else if(QuarkMatcher::Check(cjt->first->progenitor()->data()))
+ outgoing.push_back(cjt->first->progenitor());
+ }
+ assert(outgoing.size()==2&&higgs_);
+ // match up the quarks
+ for(unsigned int ix=0;ix<systems_.size();++ix) {
+ if(systems_[ix].incoming->colourLine()) {
+ for(unsigned int iy=0;iy<outgoing.size();++iy) {
+ if(outgoing[iy]->colourLine()==systems_[ix].incoming->colourLine()) {
+ systems_[ix].outgoing=outgoing[iy];
+ break;
+ }
+ }
+ }
+ else {
+ for(unsigned int iy=0;iy<outgoing.size();++iy) {
+ if(outgoing[iy]->antiColourLine()==systems_[ix].incoming->antiColourLine()) {
+ systems_[ix].outgoing=outgoing[iy];
+ break;
+ }
+ }
+ }
+ }
+ assert(systems_[0].outgoing&&systems_[1].outgoing);
+ assert(systems_.size()==2);
+ initial = initial_;
+ final = final_;
+}
+
+void MEPP2HiggsVBF::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
+ static const double eps = 1e-6;
+ // select emitting line
+ if(UseRandom::rndbool()) swap(systems_[0],systems_[1]);
+ // extract the born variables
+ q_[0] = systems_[0].outgoing->momentum()-systems_[0].incoming->momentum();
+ q2_[0] = -q_[0].m2();
+ Energy Q = sqrt(q2_[0]);
+ xB_[0] = systems_[0].incoming->x();
+ // construct lorentz transform from lab to breit frame
+ Lorentz5Momentum phadron = systems_[0].hadron->momentum();
+ phadron.setMass(0.*GeV);
+ phadron.rescaleEnergy();
+ Lorentz5Momentum pcmf = phadron+0.5/xB_[0]*q_[0];
+ pcmf.rescaleMass();
+ LorentzRotation rot(-pcmf.boostVector());
+ Lorentz5Momentum pbeam = rot*phadron;
+ Axis axis(pbeam.vect().unit());
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ rot.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ Lorentz5Momentum pout = rot*(systems_[1].outgoing->momentum()+higgs_->momentum());
+ rot.rotateZ(-atan2(pout.y(),pout.x()));
+ // calculate the A coefficient for the correlations
+ acoeff_ = A(systems_[0].incoming->dataPtr(),systems_[0].outgoing->dataPtr(),
+ systems_[1].incoming->dataPtr(),systems_[1].outgoing->dataPtr());
+ vector<double> azicoeff;
+ // select the type of process
+ bool BGF = UseRandom::rnd()>procProb_;
+ double wgt,xp,zp,x1,x2,x3,xperp;
+ l_ = 2.*(rot*systems_[1].incoming->momentum())/Q;
+ m_ = 2.*(rot*systems_[1].outgoing->momentum())/Q;
+ // compton process
+ if(!BGF) {
+ wgt = generateComptonPoint(xp,zp);
+ if(xp<eps) return;
+ // common pieces
+ Energy2 mu2 = q2_[0]*((1.-xp)*(1-zp)*zp/xp+1);
+ wgt *= 2./3./Constants::pi*alpha_->value(mu2)/procProb_;
+ // PDF piece
+ wgt *= systems_[0].pdf->xfx(systems_[0].beam,
+ systems_[0].incoming->dataPtr(),mu2 ,xB_[0]/xp)/
+ systems_[0].pdf->xfx(systems_[0].beam,
+ systems_[0].incoming->dataPtr(),scale(),xB_[0] );
+ // numerator factors
+ wgt /= (1.-xp)*(1.-zp);
+ // other bits
+ xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ x1 = -1./xp;
+ x2 = 1.-(1.-zp)/xp;
+ x3 = 2.+x1-x2;
+ // matrix element pieces
+ azicoeff = ComptonME(xp,x2,xperp,l_,m_);
+ }
+ else {
+ wgt = generateBGFPoint(xp,zp);
+ if(xp<1e-6) return;
+ // common pieces
+ Energy2 mu2 = q2_[0]*((1.-xp)*(1-zp)*zp/xp+1);
+ wgt *= 0.25/Constants::pi*alpha_->value(mu2)/(1.-procProb_);
+ // PDF piece
+ wgt *= systems_[0].pdf->xfx(systems_[0].beam,
+ gluon_ ,mu2 ,xB_[0]/xp)/
+ systems_[0].pdf->xfx(systems_[0].beam,
+ systems_[0].incoming->dataPtr(),scale(),xB_[0] );
+ // numerator factors
+ wgt /= (1.-zp);
+ // other bits
+ xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ x1 = -1./xp;
+ x2 = 1.-(1.-zp)/xp;
+ x3 = 2.+x1-x2;
+ // matrix element pieces
+ azicoeff = BGFME(xp,x2,x3,xperp,l_,m_);
+ }
+ // compute the azimuthal average of the weight
+ wgt *= azicoeff[0]+0.5*(azicoeff[2]+azicoeff[4]);
+ // finally factor as picked one line
+ wgt *= 2.;
+ // decide whether or not to accept the weight
+ if(UseRandom::rnd()>wgt) return;
+ // if accepted generate generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi)),sphi(sin(phi));
+ phiwgt = azicoeff[0]+azicoeff[5]*sphi*cphi
+ +azicoeff[1]*cphi+azicoeff[2]*sqr(cphi)
+ +azicoeff[3]*sphi+azicoeff[4]*sqr(sphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in VBFMECorrection"
+ << "::applyHardMatrixElementCorrection() to"
+ << " generate phi" << Exception::eventerror;
+ // compute the new incoming and outgoing momenta
+ Lorentz5Momentum p1 = Lorentz5Momentum( 0.5*Q*xperp*cos(phi), 0.5*Q*xperp*sin(phi),
+ -0.5*Q*x2,0.*GeV,0.*GeV);
+ p1.rescaleEnergy();
+ Lorentz5Momentum p2 = Lorentz5Momentum(-0.5*Q*xperp*cos(phi),-0.5*Q*xperp*sin(phi),
+ -0.5*Q*x3,0.*GeV,0.*GeV);
+ p2.rescaleEnergy();
+ Lorentz5Momentum pin(0.*GeV,0.*GeV,-0.5*x1*Q,-0.5*x1*Q,0.*GeV);
+ // debugging code to test vs helicity amplitude expression for matrix elements
+// double cphi(cos(phi)),sphi(sin(phi));
+// double old = (azicoeff[0]+azicoeff[5]*sphi*cphi
+// +azicoeff[1]*cphi+azicoeff[2]*sqr(cphi)
+// +azicoeff[3]*sphi+azicoeff[4]*sqr(sphi));
+// if(!BGF) {
+// old *= 8.*Constants::pi/(1.-xp)/(1.-zp);
+// }
+// else {
+// old *= 8.*Constants::pi/zp/(1.-zp);
+// }
+// debuggingMatrixElement(BGF,
+// systems_[0].incoming->dataPtr(),
+// systems_[0].outgoing->dataPtr(),
+// systems_[1].incoming->dataPtr(),
+// systems_[1].outgoing->dataPtr(),
+// rot*systems_[0].incoming->momentum(),
+// rot*systems_[0].outgoing->momentum(),
+// rot*systems_[1].incoming->momentum(),
+// rot*systems_[1].outgoing->momentum(),
+// pin,p1,p2,rot*higgs_->momentum(),
+// q2_[0],scale(),old);
+ // we need inverse of the rotation, i.e back to lab from breit
+ rot.invert();
+ // transform the momenta to lab frame
+ pin *= rot;
+ p1 *= rot;
+ p2 *= rot;
+ // test to ensure outgoing particles can be put on-shell
+ if(!BGF) {
+ if(p1.e()<systems_[0].outgoing->dataPtr()->constituentMass()) return;
+ if(p2.e()<gluon_ ->constituentMass()) return;
+ }
+ else {
+ if(p1.e()<systems_[0].outgoing->dataPtr() ->constituentMass()) return;
+ if(p2.e()<systems_[0].incoming->dataPtr()->CC()->constituentMass()) return;
+ }
+ // stats for weights > 1
+ if(wgt>1.) {
+ ++nover_;
+ if(!BGF) maxwgt_.first = max(maxwgt_.first ,wgt);
+ else maxwgt_.second = max(maxwgt_.second,wgt);
+ }
+ // create the new particles and add to ShowerTree
+ bool isquark = systems_[0].incoming->colourLine();
+ if(!BGF) {
+ PPtr newin = new_ptr(Particle(*systems_[0].incoming));
+ newin->set5Momentum(pin);
+ PPtr newg = gluon_ ->produceParticle(p2 );
+ PPtr newout = systems_[0].outgoing->dataPtr()->produceParticle(p1 );
+ ColinePtr col=isquark ?
+ systems_[0].incoming->colourLine() : systems_[0].incoming->antiColourLine();
+ ColinePtr newline=new_ptr(ColourLine());
+ // final-state emission
+ if(xp>zp) {
+ col->removeColoured(newout,!isquark);
+ col->addColoured(newin,!isquark);
+ col->addColoured(newg,!isquark);
+ newline->addColoured(newg,isquark);
+ newline->addColoured(newout,!isquark);
+ }
+ // initial-state emission
+ else {
+ col->removeColoured(newin ,!isquark);
+ col->addColoured(newout,!isquark);
+ col->addColoured(newg,isquark);
+ newline->addColoured(newg,!isquark);
+ newline->addColoured(newin,!isquark);
+ }
+ PPtr orig;
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(cit->first->progenitor()!=systems_[0].incoming) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newin);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
+ cit->first->progenitor(sp);
+ tree->incomingLines()[cit->first]=sp;
+ sp->x(xB_[0]/xp);
+ cit->first->perturbative(xp>zp);
+ if(xp<=zp) orig=cit->first->original();
+ }
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(cit->first->progenitor()!=systems_[0].outgoing) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newout);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
+ cit->first->progenitor(sp);
+ tree->outgoingLines()[cit->first]=sp;
+ cit->first->perturbative(xp<=zp);
+ if(xp>zp) orig=cit->first->original();
+ }
+ assert(orig);
+ // add the gluon
+ ShowerParticlePtr sg=new_ptr(ShowerParticle(*newg,1,true));
+ ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
+ gluon->perturbative(false);
+ tree->outgoingLines().insert(make_pair(gluon,sg));
+ tree->hardMatrixElementCorrection(true);
+ }
+ else {
+ PPtr newin = gluon_ ->produceParticle(pin);
+ PPtr newqbar = systems_[0].incoming->dataPtr()->CC()->produceParticle(p2 );
+ PPtr newout = systems_[0].outgoing->dataPtr() ->produceParticle(p1 );
+ ColinePtr col=isquark ? systems_[0].incoming->colourLine() : systems_[0].incoming->antiColourLine();
+ ColinePtr newline=new_ptr(ColourLine());
+ col ->addColoured(newin ,!isquark);
+ newline->addColoured(newin , isquark);
+ col ->addColoured(newout ,!isquark);
+ newline->addColoured(newqbar, isquark);
+ PPtr orig;
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(cit->first->progenitor()!=systems_[0].incoming) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newin);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
+ cit->first->progenitor(sp);
+ tree->incomingLines()[cit->first]=sp;
+ sp->x(xB_[0]/xp);
+ cit->first->perturbative(false);
+ orig=cit->first->original();
+ }
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(cit->first->progenitor()!=systems_[0].outgoing) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newout);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
+ cit->first->progenitor(sp);
+ tree->outgoingLines()[cit->first]=sp;
+ cit->first->perturbative(true);
+ }
+ assert(orig);
+ // add the (anti)quark
+ ShowerParticlePtr sqbar=new_ptr(ShowerParticle(*newqbar,1,true));
+ ShowerProgenitorPtr qbar=new_ptr(ShowerProgenitor(orig,newqbar,sqbar));
+ qbar->perturbative(false);
+ tree->outgoingLines().insert(make_pair(qbar,sqbar));
+ tree->hardMatrixElementCorrection(true);
+ }
+}
+
+double MEPP2HiggsVBF::A(tcPDPtr qin1, tcPDPtr qout1,
+ tcPDPtr qin2, tcPDPtr ) {
+ double output;
+ // charged current
+ if(qin1->id()!=qout1->id()) {
+ output = 2;
+ }
+ // neutral current
+ else {
+ double cvl,cal,cvq,caq;
+ if(abs(qin2->id())%2==0) {
+ cvl = generator()->standardModel()->vu();
+ cal = generator()->standardModel()->au();
+ }
+ else {
+ cvl = generator()->standardModel()->vd();
+ cal = generator()->standardModel()->ad();
+ }
+ if(abs(qin1->id())%2==0) {
+ cvq = generator()->standardModel()->vu();
+ caq = generator()->standardModel()->au();
+ }
+ else {
+ cvq = generator()->standardModel()->vd();
+ caq = generator()->standardModel()->ad();
+ }
+ output = 8.*cvl*cal*cvq*caq/(sqr(cvl)+sqr(cal))/(sqr(cvq)+sqr(caq));
+ }
+ if(qin1->id()<0) output *= -1.;
+ if(qin2->id()<0) output *= -1;
+ return output;
+}
+
+double MEPP2HiggsVBF::generateComptonPoint(double &xp, double & zp) {
+ static const double maxwgt = 50.;
+ double wgt,xperp2,x2;
+ do {
+ xp = UseRandom::rnd();
+ double zpmin = xp, zpmax = 1./(1.+xp*(1.-xp));
+ zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+ wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+ if(UseRandom::rndbool()) swap(xp,zp);
+ xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
+ x2 = 1.-(1.-zp)/xp;
+ wgt *= 2.*(1.+sqr(xp)*(sqr(x2)+1.5*xperp2))/(1.-xp)/(1.-zp);
+ if(wgt>maxwgt)
+ if(wgt>maxwgt) {
+ ostringstream wstring;
+ wstring << "MEPP2HiggsVBF::generateComptonPoint() "
+ << "Weight greater than maximum"
+ << "wgt = " << wgt << " maxwgt = " << maxwgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(wgt<UseRandom::rnd()*maxwgt);
+ return comptonInt_/((1.+sqr(xp)*(sqr(x2)+1.5*xperp2))/(1.-xp)/(1.-zp));
+}
+
+double MEPP2HiggsVBF::generateBGFPoint(double &xp, double & zp) {
+ static const double maxwgt = 25.;
+ double wgt;
+ double x2,x3,xperp2;
+ do {
+ xp = UseRandom::rnd();
+ double zpmax = 1./(1.+xp*(1.-xp)), zpmin = 1.-zpmax;
+ zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+ wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+ double x1 = -1./xp;
+ x2 = 1.-(1.-zp)/xp;
+ x3 = 2.+x1-x2;
+ xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
+ wgt *= sqr(xp)/(1.-zp)*(sqr(x3)+sqr(x2)+3.*xperp2);
+ if(wgt>maxwgt) {
+ ostringstream wstring;
+ wstring << "DISBase::generateBGFPoint "
+ << "Weight greater than maximum "
+ << "wgt = " << wgt << " maxwgt = 1\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(wgt<UseRandom::rnd()*maxwgt);
+ return bgfInt_/sqr(xp)*(1.-zp)/(sqr(x3)+sqr(x2)+3.*xperp2);
+}
+
+bool MEPP2HiggsVBF::softMatrixElementVeto(ShowerProgenitorPtr initial,
+ ShowerParticlePtr parent,Branching br) {
+ bool veto = !UseRandom::rndbool(parent->isFinalState() ? 1./final_ : 1./initial_);
+ // check if me correction should be applied
+ long id[2]={initial->id(),parent->id()};
+ if(id[0]!=id[1]||id[1]==ParticleID::g) return veto;
+ // if not from the right side
+ if(initial->progenitor()!=systems_[0].incoming &&
+ initial->progenitor()!=systems_[0].outgoing) return veto;
+ // get the pT
+ Energy pT=br.kinematics->pT();
+ // check if hardest so far
+ if(pT<initial->highestpT()) return veto;
+ double kappa(sqr(br.kinematics->scale())/q2_[0]),z(br.kinematics->z());
+ double zk((1.-z)*kappa);
+ // final-state
+ double wgt(0.);
+ if(parent->isFinalState()) {
+ double zp=z,xp=1./(1.+z*zk);
+ double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ double x2 = 1.-(1.-zp)/xp;
+ vector<double> azicoeff = ComptonME(xp,x2,xperp,l_,m_);
+ wgt = (azicoeff[0]+0.5*azicoeff[2]+0.5*azicoeff[4])*
+ xp/(1.+sqr(z))/final_;
+ if(wgt<.0||wgt>1.) {
+ ostringstream wstring;
+ wstring << "Soft ME correction weight too large or "
+ << "negative for FSR in MEPP2HiggsVBF::"
+ << "softMatrixElementVeto() soft weight "
+ << " xp = " << xp << " zp = " << zp
+ << " weight = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ else {
+ double xp = 2.*z/(1.+zk+sqrt(sqr(1.+zk)-4.*z*zk));
+ double zp = 0.5* (1.-zk+sqrt(sqr(1.+zk)-4.*z*zk));
+ double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ double x1 = -1./xp, x2 = 1.-(1.-zp)/xp, x3 = 2.+x1-x2;
+ // compton
+ if(br.ids[0]!=ParticleID::g) {
+ vector<double> azicoeff = ComptonME(xp,x2,xperp,l_,m_);
+ wgt = (azicoeff[0]+0.5*azicoeff[2]+0.5*azicoeff[4])*
+ xp*(1.-z)/(1.-xp)/(1.+sqr(z))/(1.-zp+xp-2.*xp*(1.-zp));
+ }
+ // BGF
+ else {
+ vector<double> azicoeff = BGFME(xp,x2,x3,xperp,l_,m_);
+ wgt = (azicoeff[0]+0.5*azicoeff[2]+0.5*azicoeff[4])*
+ xp/(1.-zp+xp-2.*xp*(1.-zp))/(sqr(z)+sqr(1.-z));
+ }
+ wgt /=initial_;
+ if(wgt<.0||wgt>1.) {
+ ostringstream wstring;
+ wstring << "Soft ME correction weight too large or "
+ << "negative for ISR in MEPP2HiggsVBF::"
+ << "softMatrixElementVeto() soft weight "
+ << " xp = " << xp << " zp = " << zp
+ << " weight = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ // if not vetoed
+ if(UseRandom::rndbool(wgt)) {
+ initial->highestpT(pT);
+ return false;
+ }
+ // otherwise
+ parent->setEvolutionScale(br.kinematics->scale());
+ return true;
+}
+
+vector<double> MEPP2HiggsVBF::ComptonME(double xp, double x2, double xperp,
+ LorentzVector<double> l,
+ LorentzVector<double> m) {
+ vector<double> output(6,0.);
+ double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
+ double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
+ // no phi dependence
+ output[0] = l.t()*m.t()-l.z()*m.z()*sqr(cos2)+0.5*acoeff_*cos2*(l.t()*m.z()-l.z()*m.t());
+ // cos(phi)
+ output[1] = sin2*(-l.x()*m.t()-l.t()*m.x()
+ + 0.5*acoeff_*cos2*(l.z()*m.x()-m.z()*l.x()));
+ // cos(phi)^2
+ output[2] = +sqr(sin2)*l.x()*m.x();
+ // sin(phi)
+ output[3] = sin2*(-l.t()*m.y()-l.y()*m.t()
+ + 0.5*acoeff_*cos2*(l.z()*m.y()-m.z()*l.y()));
+ // sin(phi)^2
+ output[4] = +sqr(sin2)*l.y()*m.y();
+ // sin(phi)cos(phi)
+ output[5] = +sqr(sin2)*(m.y()*l.x()+m.x()*l.y());
+ // additional factors
+ double denom = -l.z()*m.z()+l.t()*m.t()+0.5*acoeff_*(l.t()*m.z()-l.z()*m.t());
+ double fact = sqr(xp)*(sqr(x2)+sqr(xperp))/denom;
+ for(unsigned int ix=0;ix<output.size();++ix) output[ix] *=fact;
+ output[0] += 1.;
+ return output;
+}
+
+vector<double> MEPP2HiggsVBF::BGFME(double xp, double x2, double x3,
+ double xperp,
+ LorentzVector<double> l,
+ LorentzVector<double> m) {
+ vector<double> output(6,0.);
+ double denom = -l.z()*m.z()+l.t()*m.t()+0.5*acoeff_*(l.t()*m.z()-l.z()*m.t());
+ double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
+ double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
+ double fact2 = sqr(xp)*(sqr(x2)+sqr(xperp))/denom;
+ double cos3 = x3 /sqrt(sqr(x3)+sqr(xperp));
+ double sin3 = xperp/sqrt(sqr(x3)+sqr(xperp));
+ double fact3 = sqr(xp)*(sqr(x3)+sqr(xperp))/denom;
+ // no phi dependence
+ output[0] =
+ fact2*(l.t()*m.t()-l.z()*m.z()*sqr(cos2)
+ + 0.5*acoeff_*cos2*(l.t()*m.z()-l.z()*m.t())) +
+ fact3*(l.t()*m.t()-l.z()*m.z()*sqr(cos3)
+ - 0.5*acoeff_*cos3*(l.t()*m.z()-l.z()*m.t()));
+ // cos(phi)
+ output[1] =
+ fact2*sin2*( - l.x()*m.t()-l.t()*m.x()
+ + 0.5*acoeff_*cos2*(l.z()*m.x()-m.z()*l.x())) -
+ fact3*sin3*( - l.x()*m.t()-l.t()*m.x()
+ - 0.5*acoeff_*cos3*(l.z()*m.x()-m.z()*l.x())) ;
+ // cos(phi)^2
+ output[2] = (fact2*sqr(sin2)+fact3*sqr(sin3))*l.x()*m.x();
+ // sin(phi)
+ output[3] =
+ fact2*sin2*( - l.t()*m.y()-l.y()*m.t()
+ + 0.5*acoeff_*cos2*(l.z()*m.y()-m.z()*l.y())) -
+ fact3*sin3*( - l.t()*m.y()-l.y()*m.t()
+ - 0.5*acoeff_*cos3*(l.z()*m.y()-m.z()*l.y()));
+ // sin(phi)^2
+ output[4] = (fact2*sqr(sin2)+fact3*sqr(sin3))*l.y()*m.y();
+ // sin(phi)cos(phi)
+ output[5] = (fact2*sqr(sin2)+fact3*sqr(sin3))*(m.y()*l.x()+m.x()*l.y());
+ // return the answer
+ return output;
+}
diff --git a/MatrixElement/Hadron/MEPP2HiggsVBF.h b/MatrixElement/Hadron/MEPP2HiggsVBF.h
--- a/MatrixElement/Hadron/MEPP2HiggsVBF.h
+++ b/MatrixElement/Hadron/MEPP2HiggsVBF.h
@@ -1,148 +1,500 @@
// -*- C++ -*-
#ifndef HERWIG_MEPP2HiggsVBF_H
#define HERWIG_MEPP2HiggsVBF_H
//
// This is the declaration of the MEPP2HiggsVBF class.
//
#include "Herwig++/MatrixElement/MEfftoffH.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MEPP2HiggsVBF class provides the matrix elements for the
* production of the Higgs boson via the vector boson fusion mechanism
* in hadron collisions
*
* @see \ref MEPP2HiggsVBFInterfaces "The interfaces"
* defined for MEPP2HiggsVBF.
*/
class MEPP2HiggsVBF: public MEfftoffH {
+ /**
+ * Struct to contain the hadronic system
+ */
+ struct tChannelPair{
+
+ /**
+ * The hadron
+ */
+ PPtr hadron;
+
+ /**
+ * The beam particle data object
+ */
+ tcBeamPtr beam;
+
+ /**
+ * The incoming particle
+ */
+ ShowerParticlePtr incoming;
+
+ /**
+ * The outgoing particle
+ */
+ ShowerParticlePtr outgoing;
+
+ /**
+ * The PDF
+ */
+ tcPDFPtr pdf;
+ };
+
public:
/**
* The default constructor.
*/
MEPP2HiggsVBF();
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
//@}
+ /**
+ * Virtual members to be overridden by inheriting classes
+ * which implement hard corrections
+ */
+ //@{
+ /**
+ * 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,
+ ShowerParticlePtr,Branching);
+
+ /**
+ * Apply the POWHEG style correction
+ */
+ virtual HardTreePtr generateHardest(ShowerTreePtr);
+ //@}
+
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:
+ /**
+ * Generate the hardest emission in the POWHEG approach
+ */
+ //@{
+ /**
+ * Generate a Compton process
+ */
+ void generateCompton(unsigned int system);
+
+ /**
+ * Generate a BGF process
+ */
+ void generateBGF(unsigned int system);
+
+ /**
+ * Matrix element piece for the Compton process
+ */
+ double comptonME(unsigned int system,
+ double xT,double xp, double zp, double phi);
+
+ /**
+ * Matrix element piece for the Compton process
+ */
+ double BGFME(unsigned int system,
+ double xT,double xp, double zp, double phi);
+
+ /**
+ * Leading order matrix element
+ */
+ Energy4 loMatrixElement(const Lorentz5Momentum &p1,
+ const Lorentz5Momentum &p2,
+ const Lorentz5Momentum &q1,
+ const Lorentz5Momentum &q2,
+ double G1, double G2) const;
+ //@}
+
+ /**
+ * Generate the hard emission in the old-fashioned matrix element correction approach
+ */
+ //@{
+ /**
+ * Generate the values of \f$x_p\f$ and \f$z_p\f$
+ * @param xp The value of xp, output
+ * @param zp The value of zp, output
+ */
+ double generateComptonPoint(double &xp, double & zp);
+
+ /**
+ * Generate the values of \f$x_p\f$ and \f$z_p\f$
+ * @param xp The value of xp, output
+ * @param zp The value of zp, output
+ */
+ double generateBGFPoint(double &xp, double & zp);
+
+ /**
+ * Return the coefficients for the matrix element piece for
+ * the QCD compton case. The output is the \f$a_i\f$ coefficients to
+ * give the function as
+ * \f$a_0+a_1\cos\phi+a_2\sin\phi+a_3\cos^2\phi+a_4\sin^2\phi\f$
+ * @param xp \f$x_p\f$
+ * @param x2 \f$x_2\f$
+ * @param xperp \f$x_\perp\f$
+ * @param l Scaled momentum of incoming spectator
+ * @param m Scaled momentum of outgoing spectator
+ *
+ */
+ vector<double> ComptonME(double xp, double x2, double xperp,
+ LorentzVector<double> l,
+ LorentzVector<double> m);
+
+ /**
+ * Return the coefficients for the matrix element piece for
+ * the QCD compton case. The output is the \f$a_i\f$ coefficients to
+ * give the function as
+ * \f$a_0+a_1\cos\phi+a_2\sin\phi+a_3\cos^2\phi+a_4\sin^2\phi\f$
+ * @param xp \f$x_p\f$
+ * @param x2 \f$x_3\f$
+ * @param x3 \f$x_2\f$
+ * @param xperp \f$x_\perp\f$
+ * @param l Scaled momentum of incoming spectator
+ * @param m Scaled momentum of outgoing spectator
+ *
+ */
+ vector<double> BGFME(double xp, double x2, double x3, double xperp,
+ LorentzVector<double> l,
+ LorentzVector<double> m);
+
+ /**
+ * Calculate the coefficient A for the correlations
+ */
+ double A(tcPDPtr qin1, tcPDPtr qout1, tcPDPtr qin2, tcPDPtr qout2);
+ //@}
+
+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();
+
+ /**
+ * Finalize this object. Called in the run phase just after a
+ * run has ended. Used eg. to write out statistics.
+ */
+ virtual void dofinish();
//@}
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); }
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MEPP2HiggsVBF> initMEPP2HiggsVBF;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEPP2HiggsVBF & operator=(const MEPP2HiggsVBF &);
+private:
+
+ /**
+ * Parameters for the hard POWHEG emission
+ */
+ //@{
+ /**
+ * Pointer to the object calculating the strong coupling
+ */
+ ShowerAlphaPtr alpha_;
+
+ /**
+ * Weight for the compton channel
+ */
+ double comptonWeight_;
+
+ /**
+ * Weight for the BGF channel
+ */
+ double BGFWeight_;
+
+ /**
+ * Minimum value of \f$p_T\f$
+ */
+ Energy pTmin_;
+
+ /**
+ * Gluon particle data object
+ */
+ PDPtr gluon_;
+ //@}
+
+ /**
+ * Properties of the emission
+ */
+ //@{
+ /**
+ * Beam particle
+ */
+ tcBeamPtr beam_[2];
+
+ /**
+ * PDF object
+ */
+ tcPDFPtr pdf_[2];
+
+ /**
+ * Partons
+ */
+ tcPDPtr partons_[2][4];
+
+ /**
+ * q
+ */
+ Lorentz5Momentum q_[2];
+
+ /**
+ * \f$Q^2\f$
+ */
+ Energy2 q2_[2];
+
+ /**
+ * Coupling factor
+ */
+ double acoeff_;
+
+ /**
+ * Lorentz vectors for the matrix element
+ */
+ LorentzVector<double> l_;
+
+ /**
+ * Lorentz vectors for the matrix element
+ */
+ LorentzVector<double> m_;
+
+ /**
+ * Born momentum fraction
+ */
+ double xB_[2];
+
+ /**
+ * Rotation to the Breit frame
+ */
+ LorentzRotation rot_[2];
+
+ /**
+ * Quark momenta for spectator system
+ */
+ Lorentz5Momentum pother_[2][2];
+
+ /**
+ * Quark momenta for emitting system
+ */
+ Lorentz5Momentum psystem_[2][2];
+
+ /**
+ * Higgs momenta
+ */
+ Lorentz5Momentum phiggs_[2];
+
+ /**
+ * Transverse momenta for the compton emissions
+ */
+ Energy pTCompton_[2];
+
+ /**
+ * Transverse momenta for the BGF emissions
+ */
+ Energy pTBGF_[2];
+
+ /**
+ * Whether the Compton radiation is ISR or FSR
+ */
+ bool ComptonISFS_[2];
+
+ /**
+ * Momenta of the particles for a compton emission
+ */
+ vector<Lorentz5Momentum> ComptonMomenta_[2];
+
+ /**
+ * Momenta of the particles for a BGF emission
+ */
+ vector<Lorentz5Momentum> BGFMomenta_[2];
+
+ /**
+ * the systems
+ */
+ vector<tChannelPair> systems_;
+
+ /**
+ * Higgs boson
+ */
+ ShowerParticlePtr higgs_;
+ //@}
+
+ /**
+ * Parameters for the matrix element correction
+ */
+ //@{
+ /**
+ * Enchancement factor for ISR
+ */
+ double initial_;
+
+ /**
+ * Enchancement factor for FSR
+ */
+ double final_;
+
+ /**
+ * Relative fraction of compton and BGF processes to generate
+ */
+ double procProb_;
+
+ /**
+ * Integral for compton process
+ */
+ double comptonInt_;
+
+ /**
+ * Integral for BGF process
+ */
+ double bgfInt_;
+
+ /**
+ * Number of weights greater than 1
+ */
+ unsigned int nover_;
+
+ /**
+ * Maximum weight
+ */
+ pair<double,double> maxwgt_;
+ //@}
+
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MEPP2HiggsVBF. */
template <>
struct BaseClassTrait<Herwig::MEPP2HiggsVBF,1> {
/** Typedef of the first base class of MEPP2HiggsVBF. */
typedef Herwig::MEfftoffH NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MEPP2HiggsVBF class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MEPP2HiggsVBF>
: public ClassTraitsBase<Herwig::MEPP2HiggsVBF> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MEPP2HiggsVBF"; }
/**
* The name of a file containing the dynamic library where the class
* MEPP2HiggsVBF is implemented. It may also include several, space-separated,
* libraries if the class MEPP2HiggsVBF depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwMEHadron.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MEPP2HiggsVBF_H */
diff --git a/MatrixElement/Powheg/MEPP2HiggsVBFPowheg.cc b/MatrixElement/Powheg/MEPP2HiggsVBFPowheg.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Powheg/MEPP2HiggsVBFPowheg.cc
@@ -0,0 +1,721 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the MEPP2HiggsVBFPowheg class.
+//
+
+#include "MEPP2HiggsVBFPowheg.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
+#include "Herwig++/PDT/StandardMatchers.h"
+#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
+//ACDC test
+//#include "ThePEG/Repository/UseRandom.h"
+
+using namespace Herwig;
+
+MEPP2HiggsVBFPowheg::MEPP2HiggsVBFPowheg()
+ : scaleOpt_(1), muF_(100.*GeV), scaleFact_(1.), contrib_(1), power_(0.1)
+{}
+
+int MEPP2HiggsVBFPowheg::nDim() const {
+ //ACDC test
+ // return 0;
+ return MEPP2HiggsVBF::nDim()+3;
+}
+//ACDC test (delete r)
+bool MEPP2HiggsVBFPowheg::generateKinematics(const double * r) {
+ //bool MEPP2HiggsVBFPowheg::generateKinematics(const double * ) {
+ int ndim = nDim();
+ //ACDC test
+ double r3 = r[ndim-3];
+ //double r3 = UseRandom::rnd();
+
+ // Born kinematics
+ //ACDC test (comment next line out)
+ if(!MEPP2HiggsVBF::generateKinematics(r)) return false;
+
+ // hadron and momentum fraction
+ // set Q2 process momenta
+
+ if(r3 > 0.5) {
+ r3 = 2.*(r3-0.5);
+ if(lastPartons().first ->dataPtr()==mePartonData()[0]&&
+ lastPartons().second->dataPtr()==mePartonData()[1]) {
+ _hadron = dynamic_ptr_cast<tcBeamPtr>(lastParticles().first->dataPtr());
+ _xB = lastX1();
+ }
+ else {
+ _hadron = dynamic_ptr_cast<tcBeamPtr>(lastParticles().second->dataPtr());
+ _xB = lastX2();
+ }
+ _partons[0] = mePartonData()[0];
+ _partons[1] = mePartonData()[1];
+ _partons[4] = mePartonData()[4];
+ if(!swapOrder()) {
+ _pb = meMomenta()[0];
+ _pc = meMomenta()[2];
+ _pbother = meMomenta()[1];
+ _pcother = meMomenta()[3];
+ _partons[2] = mePartonData()[2];
+ _partons[3] = mePartonData()[3];
+ }
+ else {
+ _pb = meMomenta()[0];
+ _pc = meMomenta()[3];
+ _pbother = meMomenta()[1];
+ _pcother = meMomenta()[2];
+ _partons[2] = mePartonData()[3];
+ _partons[3] = mePartonData()[2];
+ }
+ }
+else {
+ r3 = 2.*r3;
+ if(lastPartons().first ->dataPtr()==mePartonData()[0]&&
+ lastPartons().second->dataPtr()==mePartonData()[1]) {
+ _hadron = dynamic_ptr_cast<tcBeamPtr>(lastParticles().second->dataPtr());
+ _xB = lastX2();
+ }
+ else {
+ _hadron = dynamic_ptr_cast<tcBeamPtr>(lastParticles().first->dataPtr());
+ _xB = lastX1();
+ }
+ _partons[0] = mePartonData()[1];
+ _partons[1] = mePartonData()[0];
+ _partons[4] = mePartonData()[4];
+ if(!swapOrder()) {
+ _pb = meMomenta()[1];
+ _pc = meMomenta()[3];
+ _pbother = meMomenta()[0];
+ _pcother = meMomenta()[2];
+ _partons[2] = mePartonData()[3];
+ _partons[3] = mePartonData()[2];
+ }
+ else {
+ _pb = meMomenta()[1];
+ _pc = meMomenta()[2];
+ _pbother = meMomenta()[0];
+ _pcother = meMomenta()[3];
+ _partons[2] = mePartonData()[2];
+ _partons[3] = mePartonData()[3];
+ }
+ }
+ // LO Momenta assignment
+ _loMomenta[0] = _pb;
+ _loMomenta[1] = _pbother;
+ _loMomenta[2] = _pc;
+ _loMomenta[3] = _pcother;
+ _pa = _pc-_pb;
+ // xp
+ double rhomin = pow(1.-_xB,1.-power_);
+
+ //ACDC test
+ double rho = r[ndim-1]*rhomin;
+ //double rho = UseRandom::rnd()*rhomin;
+
+ _xp = 1.-pow(rho,1./(1.-power_));
+ // zp
+ //ACDC test
+ _zp = r[ndim-2];
+ //_zp = UseRandom::rnd();
+
+ // phi
+ _phi = r3*Constants::twopi;
+ jac_ = rhomin/(1.-power_)*pow(1.-_xp,power_);
+ return true;
+}
+
+IBPtr MEPP2HiggsVBFPowheg::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr MEPP2HiggsVBFPowheg::fullclone() const {
+ return new_ptr(*this);
+}
+
+
+void MEPP2HiggsVBFPowheg::persistentOutput(PersistentOStream & os) const {
+ os << ounit(muF_,GeV) << scaleFact_ << scaleOpt_ << contrib_
+ << ounit(_mz2,GeV2) << ounit(_mw2,GeV2) << power_;
+}
+
+void MEPP2HiggsVBFPowheg::persistentInput(PersistentIStream & is, int) {
+ is >> iunit(muF_,GeV) >> scaleFact_ >> scaleOpt_ >> contrib_
+ >> iunit(_mz2,GeV2) >> iunit(_mw2,GeV2) >> power_;
+}
+
+ClassDescription<MEPP2HiggsVBFPowheg>
+ MEPP2HiggsVBFPowheg::initMEPP2HiggsVBFPowheg;
+// Definition of the static class description member.
+
+void MEPP2HiggsVBFPowheg::Init() {
+
+ static ClassDocumentation<MEPP2HiggsVBFPowheg> documentation
+ ("The MENeutralCurrentDISPowheg class implements the NLO matrix element"
+ " for neutral current DIS in the Powheg scheme.");
+
+ static Switch<MEPP2HiggsVBFPowheg,unsigned int> interfaceContribution
+ ("Contribution",
+ "Which contributions to the cross section to include",
+ &MEPP2HiggsVBFPowheg::contrib_, 1, false, false);
+ static SwitchOption interfaceContributionLeadingOrder
+ (interfaceContribution,
+ "LeadingOrder",
+ "Just generate the leading order cross section",
+ 0);
+ static SwitchOption interfaceContributionPositiveNLO
+ (interfaceContribution,
+ "PositiveNLO",
+ "Generate the positive contribution to the full NLO cross section",
+ 1);
+ static SwitchOption interfaceContributionNegativeNLO
+ (interfaceContribution,
+ "NegativeNLO",
+ "Generate the negative contribution to the full NLO cross section",
+ 2);
+
+ static Switch<MEPP2HiggsVBFPowheg,unsigned int> interfaceScaleOption
+ ("ScaleOption",
+ "Option for the choice of factorization (and renormalization) scale",
+ &MEPP2HiggsVBFPowheg::scaleOpt_, 1, false, false);
+ static SwitchOption interfaceDynamic
+ (interfaceScaleOption,
+ "Dynamic",
+ "Dynamic factorization scale equal to the current sqrt(sHat())",
+ 1);
+ static SwitchOption interfaceFixed
+ (interfaceScaleOption,
+ "Fixed",
+ "Use a fixed factorization scale set with FactorizationScaleValue",
+ 2);
+
+ static Parameter<MEPP2HiggsVBFPowheg,Energy> interfaceFactorizationScale
+ ("FactorizationScale",
+ "Value to use in the event of a fixed factorization scale",
+ &MEPP2HiggsVBFPowheg::muF_, GeV, 100.0*GeV, 1.0*GeV, 500.0*GeV,
+ true, false, Interface::limited);
+
+ static Parameter<MEPP2HiggsVBFPowheg,double> interfaceScaleFactor
+ ("ScaleFactor",
+ "The factor used before Q2 if using a running scale",
+ &MEPP2HiggsVBFPowheg::scaleFact_, 1.0, 0.0, 10.0,
+ false, false, Interface::limited);
+
+ static Parameter<MEPP2HiggsVBFPowheg,double> interfaceSamplingPower
+ ("SamplingPower",
+ "Power for the sampling of xp",
+ &MEPP2HiggsVBFPowheg::power_, 0.6, 0.0, 1.,
+ false, false, Interface::limited);
+
+}
+
+Energy2 MEPP2HiggsVBFPowheg::scale() const {
+ return scaleOpt_ == 1 ?
+ sqr(scaleFact_)*MEPP2HiggsVBF::scale() : sqr(scaleFact_*muF_);
+}
+
+CrossSection MEPP2HiggsVBFPowheg::dSigHatDR() const {
+ return NLOWeight()*MEPP2HiggsVBF::dSigHatDR();
+}
+
+double MEPP2HiggsVBFPowheg::NLOWeight() const {
+ // If only leading order is required return 1:
+ if(contrib_==0) return 1.;
+ // Boost
+ Axis axis(_pa.vect().unit());
+ LorentzRotation rot;
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ rot = LorentzRotation();
+ if(axis.perp2()>1e-20) {
+ rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ rot.rotateX(Constants::pi);
+ }
+ if(abs(1.-_pa.e()/_pa.vect().mag())>1e-6)
+ rot.boostZ(_pa.e()/_pa.vect().mag());
+ _pb *= rot;
+ if(_pb.perp2()/GeV2>1e-20) {
+ Boost trans = -1./_pb.e()*_pb.vect();
+ trans.setZ(0.);
+ rot.boost(trans);
+ }
+ // momenta of particles
+ Lorentz5Momentum p1, p2,p1other,p2other;
+ _pa *= rot;
+ _q2 = -_pa.m2();
+ p1 = rot*_loMomenta[0];
+ p2 = rot*_loMomenta[2];
+ p1other = rot*_loMomenta[1];
+ p2other = rot*_loMomenta[3];
+
+ // scale and prefactors
+ Energy2 mu2 = scale();
+ Energy Q = sqrt(_q2);
+ // double aS = 2.*SM().alphaS(mu2);
+ double aS = 2*0.113291076960184;
+ double CFfact = 4./3.*aS/Constants::twopi;
+ double TRfact = 1./2.*aS/Constants::twopi;
+ double wgt = 0.;
+ // Breit frame variables
+ double x1 = -1./_xp;
+ double x2 = 1.-(1.-_zp)/_xp;
+ double x3 = 2.+x1-x2;
+ double xT = 2*sqrt((1.-_xp)*(1.-_zp)*_zp/_xp);
+
+ vector<Lorentz5Momentum> nloMomenta;
+ nloMomenta.resize(3);
+
+ nloMomenta[0] = Lorentz5Momentum(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ nloMomenta[1] = Lorentz5Momentum( 0.5*Q*xT*cos(_phi), 0.5*Q*xT*sin(_phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ nloMomenta[2] = Lorentz5Momentum(-0.5*Q*xT*cos(_phi), -0.5*Q*xT*sin(_phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum qnlo = nloMomenta[2]+nloMomenta[1]-nloMomenta[0];
+ Lorentz5Momentum r1 = -nloMomenta[0]/x1;
+ Lorentz5Momentum r2 = nloMomenta[1]/x2;
+ Lorentz5Momentum r3 = -nloMomenta[2]/x3;
+
+ // LO + dipole subtracted virtual + collinear quark bit with LO pdf
+ double virt = 1.+CFfact*(-4.5-1./3.*sqr(Constants::pi)+
+ 1.5*log(_q2/mu2/(1.-_xB))+
+ 2.*log(1.-_xB)*log(_q2/mu2)+
+ sqr(log(1.-_xB)));
+ // PDF from leading-order
+ double loPDF =
+ _hadron->pdf()->xfx(_hadron,_partons[0],mu2,_xB)/_xB;
+ // NLO gluon PDF
+
+ tcPDPtr gluon = getParticleData(ParticleID::g);
+ double gPDF =
+ _hadron->pdf()->xfx(_hadron,gluon,mu2,_xB/_xp)*_xp/_xB;
+ // NLO quark PDF
+ double qPDF =
+ _hadron->pdf()->xfx(_hadron,_partons[0],mu2,_xB/_xp)*_xp/_xB;
+ // collinear counterterms
+ // gluon
+ double collg =
+ TRfact/_xp*gPDF*(2.*_xp*(1.-_xp)+(sqr(_xp)+sqr(1.-_xp))*
+ log((1.-_xp)*_q2/_xp/mu2));
+ // quark
+ double collq =
+ CFfact/_xp*qPDF*(1-_xp-2./(1.-_xp)*log(_xp)-(1.+_xp)*
+ log((1.-_xp)/_xp*_q2/mu2))+
+ CFfact/_xp*(qPDF-_xp*loPDF)*(2./(1.-_xp)*
+ log(_q2*(1.-_xp)/mu2)-1.5/(1.-_xp));
+ // Electroweak coefficients
+ double c0L,c1L,c0R,c1R;
+ Energy2 mb2;
+ // W
+ if(_partons[0]->id()!=_partons[2]->id()) {
+ mb2 = _mw2;
+ c0L = sqrt(0.5);
+ c0R = 0;
+ c1L = sqrt(0.5);
+ c1R = 0;
+ }
+ // Z
+ else {
+ mb2 = _mz2;
+ if(abs(_partons[0]->id())%2==0) {
+ c0L =
+ generator()->standardModel()->vu()+
+ generator()->standardModel()->au();
+ c0R =
+ generator()->standardModel()->vu()-
+ generator()->standardModel()->au();
+ }
+ else {
+ c0L =
+ generator()->standardModel()->vd()+
+ generator()->standardModel()->ad();
+ c0R =
+ generator()->standardModel()->vd()-
+ generator()->standardModel()->ad();
+ }
+ if(abs(_partons[1]->id())%2==0) {
+ c1L =
+ generator()->standardModel()->vu()+
+ generator()->standardModel()->au();
+ c1R =
+ generator()->standardModel()->vu()-
+ generator()->standardModel()->au();
+ }
+ else {
+ c1L =
+ generator()->standardModel()->vd()+
+ generator()->standardModel()->ad();
+ c1R =
+ generator()->standardModel()->vd()-
+ generator()->standardModel()->ad();
+ }
+ c0L *= 0.25;
+ c0R *= 0.25;
+ c1L *= 0.25;
+ c1R *= 0.25;
+ }
+ // Matrix element variables
+ double G1 = sqr(c0L*c1L)+sqr(c0R*c1R);
+ double G2 = sqr(c0L*c1R)+sqr(c0R*c1L);
+ Energy4 term1,term2,term3,loME;
+ if(_partons[0]->id()>0) {
+ if(_partons[1]->id()>0) {
+ term1 = loMatrixElement(r1 ,p1other,qnlo+r1,p2other,G1,G2);
+ term2 = loMatrixElement(r2-qnlo,p1other,r2 ,p2other,G1,G2);
+ term3 = loMatrixElement(r3 ,p1other,qnlo+r3,p2other,G1,G2);
+ loME = loMatrixElement(p1 ,p1other,p2 ,p2other,G1,G2);
+ }
+ else {
+ term1 = loMatrixElement(r1 ,p2other,qnlo+r1,p1other,G1,G2);
+ term2 = loMatrixElement(r2-qnlo,p2other,r2 ,p1other,G1,G2);
+ term3 = loMatrixElement(r3 ,p2other,qnlo+r3,p1other,G1,G2);
+ loME = loMatrixElement(p1 ,p2other,p2 ,p1other,G1,G2);
+ }
+ }
+ else {
+ if(_partons[1]->id()>0) {
+ term1 = loMatrixElement(qnlo+r1,p1other,r1 ,p2other,G1,G2);
+ term2 = loMatrixElement(r2 ,p1other,r2-qnlo,p2other,G1,G2);
+ term3 = loMatrixElement(qnlo+r3,p1other,r3 ,p2other,G1,G2);
+ loME = loMatrixElement(p2 ,p1other,p1 ,p2other,G1,G2);
+ }
+ else {
+ term1 = loMatrixElement(qnlo+r1,p2other,r1 ,p1other,G1,G2);
+ term2 = loMatrixElement(r2 ,p2other,r2-qnlo,p1other,G1,G2);
+ term3 = loMatrixElement(qnlo+r3,p2other,r3 ,p1other,G1,G2);
+ loME = loMatrixElement(p2 ,p2other,p1 ,p1other,G1,G2);
+ }
+ }
+ if(1-_xp > 1e-10 && 1.-_zp > 1e-10){
+ // q -> qg term
+ double real1 = (term1+sqr(_xp)*sqr(x2)*term2)/loME;
+ double dipole1 = (sqr(_xp)+sqr(_zp));
+ double realq =
+ CFfact*qPDF/loPDF/_xp/((1.-_xp)*(1.-_zp))*(real1-dipole1);
+
+ // g -> q qbar term
+ double real2 = sqr(_xp)/loME*(sqr(x2)*term2+sqr(x3)*term3);
+ double dipole2 = sqr(_xp)+sqr(1.-_xp);
+ double realg = TRfact*gPDF/loPDF/_xp/(1.-_zp)*(real2-dipole2);
+
+ // Return The Full Result
+ wgt = virt+jac_*((collq+collg)/loPDF+realq+realg);
+ }
+
+
+ return contrib_ == 1 ? max(0.,wgt) : max(0.,-wgt);
+}
+
+void MEPP2HiggsVBFPowheg::doinit() {
+ MEPP2HiggsVBF::doinit();
+ // electroweak parameters
+ _mz2 = sqr(getParticleData(ParticleID::Z0)->mass());
+ _mw2 = sqr(getParticleData(ParticleID::Wplus)->mass());
+ tcPDPtr gluon = getParticleData(ParticleID::g);
+}
+
+Energy4 MEPP2HiggsVBFPowheg::
+loMatrixElement(const Lorentz5Momentum &p1,
+ const Lorentz5Momentum &p2,
+ const Lorentz5Momentum &q1,
+ const Lorentz5Momentum &q2,
+ double G1, double G2) const {
+ return G1*(p1*p2)*(q1*q2) + G2*(p1*q2)*(q1*p2);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+// double g;
+// g = 1./sqrt(SM().sin2ThetaW());
+// g = 1./sqrt((1.-SM().sin2ThetaW())*SM().sin2ThetaW());
+
+
+// vector<SpinorWaveFunction> f1,f2;
+// vector<SpinorBarWaveFunction> a1,a2;
+// Lorentz5Momentum phiggs = rot*meMomenta()[4];
+// ScalarWaveFunction higgs(phiggs,mePartonData()[4],1.,outgoing);
+// SpinorWaveFunction fin1,fin2;
+// SpinorBarWaveFunction ain1,ain2;
+// if(_partons[0]->id()>0) {
+// fin1 = SpinorWaveFunction(nloMomenta[0],_partons[0],incoming);
+// ain1 = SpinorBarWaveFunction(nloMomenta[1],_partons[2],outgoing);
+// }
+// else {
+// fin1 = SpinorWaveFunction(nloMomenta[1],_partons[2],outgoing);
+// ain1 = SpinorBarWaveFunction(nloMomenta[0],_partons[0],incoming);
+// }
+// if(_partons[1]->id()>0) {
+// fin2 = SpinorWaveFunction(p1other,_partons[1],incoming);
+// ain2 = SpinorBarWaveFunction(p2other,_partons[3],outgoing);
+// }
+// else {
+// fin2 = SpinorWaveFunction(p2other,_partons[3],outgoing);
+// ain2 = SpinorBarWaveFunction(p1other,_partons[1],incoming);
+// }
+// VectorWaveFunction gwave(nloMomenta[2],gluon,outgoing);
+// vector<VectorWaveFunction> g1;
+// for(unsigned int ix=0;ix<2;++ix) {
+// fin1.reset(ix); f1.push_back(fin1);
+// fin2.reset(ix); f2.push_back(fin2);
+// ain1.reset(ix); a1.push_back(ain1);
+// ain2.reset(ix); a2.push_back(ain2);
+// gwave.reset(2*ix); g1.push_back(gwave);
+// }
+// AbstractFFVVertexPtr vertex[2];
+// tcPDPtr vec[2];
+// for(unsigned int ix=0;ix<2;++ix) {
+// int icharge;
+// icharge = _partons[ix]->iCharge()-_partons[ix+2]->iCharge();
+// if(icharge==0) vec[ix] = _z0;
+// else if(icharge>0) vec[ix] = _wplus;
+// else vec[ix] = _wminus;
+// vertex[ix] = vec[ix]==_z0 ? _vertexFFZ : _vertexFFW;
+// }
+// tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
+// AbstractFFVVertexPtr gvertex = hwsm->vertexFFG();
+// VectorWaveFunction inter[2];
+// Complex diag[2];
+// double me(0.);
+// nloMomenta[0].setMass(ZERO);
+// nloMomenta[1].setMass(ZERO);
+// nloMomenta[2].setMass(ZERO);
+// for(unsigned int i1=0;i1<2;++i1) {
+// unsigned int i2=i1;
+// // wavefunction for the 1st intermediate vector
+// inter[0] = vertex[1]->evaluate(scale(),1,vec[1],f2[i1],a2[i2]);
+// for(unsigned int i3=0;i3<2;++i3) {
+// unsigned int i4=i3;
+// for(unsigned int ig=0;ig<2;++ig) {
+// SpinorWaveFunction soff =
+// gvertex->evaluate(scale(),5,f1[i3].getParticle(),
+// f1[i3],g1[ig]);
+// SpinorBarWaveFunction aoff =
+// gvertex->evaluate(scale(),5,a1[i4].getParticle(),
+// a1[i4],g1[ig]);
+// inter[1] = vertex[0]->evaluate(scale(),1,vec[0],soff,a1[i4]);
+// diag[0] = _vertexWWH->evaluate(scale(),inter[0],inter[1],higgs);
+// inter[1] = vertex[0]->evaluate(scale(),1,vec[0],f1[i3],aoff);
+// diag[1] = _vertexWWH->evaluate(scale(),inter[0],inter[1],higgs);
+// // cerr << "testing helicity "
+// // << i1 << " " << i2 << " " << i3 << " " << i4 << " "
+// // << ig << " " << (diag[0]+diag[1])/(abs(diag[0])+abs(diag[1]))
+// // << "\n";
+
+// me += norm(diag[0]+diag[1]);
+// }
+// }
+// }
+// // spin factor
+// me *=0.25;
+
+
+
+
+
+
+// cerr << "testing the quark emission "
+// << fact*8.*Constants::pi*SM().alphaS(scale())/(-1.-x1)/(1.-x2)/_q2
+// *(sqr(x1)*term1+sqr(x2)*term2)*sqr(MeV2)/me << "\n";
+
+
+
+
+
+
+
+
+
+
+// vector<SpinorWaveFunction> f1,f2;
+// vector<SpinorBarWaveFunction> a1,a2;
+// Lorentz5Momentum phiggs = rot*meMomenta()[4];
+// ScalarWaveFunction higgs(phiggs,mePartonData()[4],1.,outgoing);
+// SpinorWaveFunction fin1,fin2;
+// SpinorBarWaveFunction ain1,ain2;
+// if(_partons[0]->id()>0) {
+// fin1 = SpinorWaveFunction(nloMomenta[2],_partons[0]->CC(),outgoing);
+// ain1 = SpinorBarWaveFunction(nloMomenta[1],_partons[2],outgoing);
+// }
+// else {
+// fin1 = SpinorWaveFunction(nloMomenta[1],_partons[2],outgoing);
+// ain1 = SpinorBarWaveFunction(nloMomenta[2],_partons[0]->CC(),outgoing);
+// }
+// if(_partons[1]->id()>0) {
+// fin2 = SpinorWaveFunction(p1other,_partons[1],incoming);
+// ain2 = SpinorBarWaveFunction(p2other,_partons[3],outgoing);
+// }
+// else {
+// fin2 = SpinorWaveFunction(p2other,_partons[3],outgoing);
+// ain2 = SpinorBarWaveFunction(p1other,_partons[1],incoming);
+// }
+// VectorWaveFunction gwave(nloMomenta[0],gluon,incoming);
+// vector<VectorWaveFunction> g1;
+// for(unsigned int ix=0;ix<2;++ix) {
+// fin1.reset(ix); f1.push_back(fin1);
+// fin2.reset(ix); f2.push_back(fin2);
+// ain1.reset(ix); a1.push_back(ain1);
+// ain2.reset(ix); a2.push_back(ain2);
+// gwave.reset(2*ix);
+// g1.push_back(gwave);
+// }
+// AbstractFFVVertexPtr vertex[2];
+// tcPDPtr vec[2];
+// for(unsigned int ix=0;ix<2;++ix) {
+// int icharge;
+// icharge = _partons[ix]->iCharge()-_partons[ix+2]->iCharge();
+// if(icharge==0) vec[ix] = _z0;
+// else if(icharge>0) vec[ix] = _wplus;
+// else vec[ix] = _wminus;
+// vertex[ix] = vec[ix]==_z0 ? _vertexFFZ : _vertexFFW;
+// }
+// tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
+// AbstractFFVVertexPtr gvertex = hwsm->vertexFFG();
+// VectorWaveFunction inter[2];
+// Complex diag[2];
+// double me(0.);
+// nloMomenta[0].setMass(ZERO);
+// nloMomenta[1].setMass(ZERO);
+// nloMomenta[2].setMass(ZERO);
+// unsigned int o[2]={1,0};
+// for(unsigned int i1=0;i1<2;++i1) {
+// unsigned int i2=i1;
+// // wavefunction for the 1st intermediate vector
+// inter[0] = vertex[1]->evaluate(scale(),1,vec[1],f2[i1],a2[i2]);
+// for(unsigned int i3=0;i3<2;++i3) {
+// unsigned int i4=o[i3];
+// for(unsigned int ig=0;ig<2;++ig) {
+// SpinorWaveFunction soff =
+// gvertex->evaluate(scale(),5,f1[i3].getParticle(),
+// f1[i3],g1[ig]);
+// SpinorBarWaveFunction aoff =
+// gvertex->evaluate(scale(),5,a1[i4].getParticle(),
+// a1[i4],g1[ig]);
+// inter[1] = vertex[0]->evaluate(scale(),1,vec[0],soff,a1[i4]);
+// diag[0] = _vertexWWH->evaluate(scale(),inter[0],inter[1],higgs);
+// inter[1] = vertex[0]->evaluate(scale(),1,vec[0],f1[i3],aoff);
+// diag[1] = _vertexWWH->evaluate(scale(),inter[0],inter[1],higgs);
+// // cerr << "testing helicity "
+// // << i1 << " " << i2 << " " << i3 << " " << i4 << " "
+// // << ig << " " << (diag[0]+diag[1])/(abs(diag[0])+abs(diag[1]))
+// // << "\n";
+
+// me += norm(diag[0]+diag[1]);
+// }
+// }
+// }
+// // spin factor
+// me *=0.25;
+// cerr << "testing the gluon emission "
+// << fact*8.*Constants::pi*SM().alphaS(scale())/(1.-x2)/(1.-x3)/_q2
+// *(sqr(x3)*term3+sqr(x2)*term2)*sqr(MeV2)/me << "\n";
+
+
+
+
+
+
+
+
+
+
+
+
+
+// Energy2 D1 = -_q2-mb2;
+// Energy2 D2 = (p1other-p2other).m2()-mb2;
+// double e = sqrt(4.*Constants::pi*SM().alphaEM(scale()));
+
+// InvEnergy6 fact = 4.*pow(e*g,6)*mb2/sqr(D1)/sqr(D2);
+
+// cerr << "testing LO ME in NLO code "
+// << fact*loME*MeV2/_mestore << "\n";
+
+
+
+// vector<SpinorWaveFunction> f1,f2;
+// vector<SpinorBarWaveFunction> a1,a2;
+// Lorentz5Momentum phiggs = rot*meMomenta()[4];
+// ScalarWaveFunction higgs(phiggs,mePartonData()[4],1.,outgoing);
+// SpinorWaveFunction fin1,fin2;
+// SpinorBarWaveFunction ain1,ain2;
+// if(_partons[0]->id()>0) {
+// fin1 = SpinorWaveFunction(p1,_partons[0],incoming);
+// ain1 = SpinorBarWaveFunction(p2,_partons[2],outgoing);
+// }
+// else {
+// fin1 = SpinorWaveFunction(p2,_partons[2],outgoing);
+// ain1 = SpinorBarWaveFunction(p1,_partons[0],incoming);
+// }
+// if(_partons[1]->id()>0) {
+// fin2 = SpinorWaveFunction(p1other,_partons[1],incoming);
+// ain2 = SpinorBarWaveFunction(p2other,_partons[3],outgoing);
+// }
+// else {
+// fin2 = SpinorWaveFunction(p2other,_partons[3],outgoing);
+// ain2 = SpinorBarWaveFunction(p1other,_partons[1],incoming);
+// }
+// for(unsigned int ix=0;ix<2;++ix) {
+// fin1.reset(ix); f1.push_back(fin1);
+// fin2.reset(ix); f2.push_back(fin2);
+// ain1.reset(ix); a1.push_back(ain1);
+// ain2.reset(ix); a2.push_back(ain2);
+// }
+// AbstractFFVVertexPtr vertex[2];
+// tcPDPtr vec[2];
+// for(unsigned int ix=0;ix<2;++ix) {
+// int icharge;
+// icharge = _partons[ix]->iCharge()-_partons[ix+2]->iCharge();
+// if(icharge==0) vec[ix] = _z0;
+// else if(icharge>0) vec[ix] = _wplus;
+// else vec[ix] = _wminus;
+// vertex[ix] = vec[ix]==_z0 ? _vertexFFZ : _vertexFFW;
+// }
+// VectorWaveFunction inter[2];
+// Complex diag;
+// double me(0.);
+// for(unsigned int i1=0;i1<2;++i1) {
+// for(unsigned int i2=0;i2<2;++i2) {
+// // wavefunction for the 1st intermediate vector
+// inter[0] = vertex[0]->evaluate(scale(),1,vec[1],f2[i1],a2[i2]);
+// for(unsigned int i3=0;i3<2;++i3) {
+// for(unsigned int i4=0;i4<2;++i4) {
+// // wavefunction for the 2nd intermediate vector
+// inter[1] = vertex[1]->evaluate(scale(),1,vec[0],f1[i3],a1[i4]);
+// // matrix element
+// diag = _vertexWWH->evaluate(scale(),inter[0],inter[1],higgs);
+// me += norm(diag);
+// }
+// }
+// }
+// }
+// // spin factor
+// me *=0.25;
+
+// cerr << "testing helicity computation " << me/_mestore << "\n";
diff --git a/MatrixElement/Powheg/MEPP2HiggsVBFPowheg.h b/MatrixElement/Powheg/MEPP2HiggsVBFPowheg.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Powheg/MEPP2HiggsVBFPowheg.h
@@ -0,0 +1,314 @@
+// -*- C++ -*-
+#ifndef HERWIG_MEPP2HiggsVBFPowheg_H
+#define HERWIG_MEPP2HiggsVBFPowheg_H
+//
+// This is the declaration of the MEPP2HiggsVBFPowheg class.
+//
+
+#include "Herwig++/MatrixElement/Hadron/MEPP2HiggsVBF.h"
+#include "ThePEG/PDF/BeamParticleData.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Typedef for BeamParticleData
+ */
+typedef Ptr<BeamParticleData>::transient_const_pointer tcBeamPtr;
+
+/**
+ * Here is the documentation of the MEPP2HiggsVBFPowheg class.
+ *
+ * @see \ref MEPP2HiggsVBFPowhegInterfaces "The interfaces"
+ * defined for MEPP2HiggsVBFPowheg.
+ */
+class MEPP2HiggsVBFPowheg: public Herwig::MEPP2HiggsVBF {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ MEPP2HiggsVBFPowheg();
+
+ /** @name Virtual functions required by the MEBase class. */
+ //@{
+
+ /**
+ * Return the scale associated with the last set phase space point.
+ */
+ virtual Energy2 scale() const;
+
+ /**
+ * The number of internal degrees of freedom used in the matrix
+ * element.
+ */
+ virtual int nDim() const;
+
+ /**
+ * Generate internal degrees of freedom given nDim() uniform
+ * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
+ * generator, the dSigHatDR should be a smooth function of these
+ * numbers, although this is not strictly necessary.
+ * @param r a pointer to the first of nDim() consecutive random numbers.
+ * @return true if the generation succeeded, otherwise false.
+ */
+ virtual bool generateKinematics(const double * r);
+
+ /**
+ * Return the matrix element squared differential in the variables
+ * given by the last call to generateKinematics().
+ */
+ virtual CrossSection dSigHatDR() const;
+ //@}
+
+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:
+
+ /**
+ * The NLO weight
+ */
+ double NLOWeight() const;
+
+ /**
+ * Leading order matrix element
+ */
+ Energy4 loMatrixElement(const Lorentz5Momentum &p1,
+ const Lorentz5Momentum &p2,
+ const Lorentz5Momentum &q1,
+ const Lorentz5Momentum &q2,
+ double G1, double G2) const;
+
+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();
+ //@}
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** 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;
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<MEPP2HiggsVBFPowheg> initMEPP2HiggsVBFPowheg;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ MEPP2HiggsVBFPowheg & operator=(const MEPP2HiggsVBFPowheg &);
+
+private:
+
+ /**
+ * The Born variables
+ */
+ //@{
+ /**
+ * \f$x_B\f$
+ */
+ double _xB;
+
+ /**
+ * Partons
+ */
+ tcPDPtr _partons[5];
+
+ mutable Energy2 _q2;
+ //@}
+
+ /**
+ * The radiative variables
+ */
+ //@{
+ /**
+ * The \f$x_p\f$ or \f$z\f$ real integration variable
+ */
+ double _xp;
+
+ /**
+ * The \f$z_p\f$ real integration variable
+ */
+ double _zp;
+
+ /**
+ * The \f$fi\f$ real integration variable
+ */
+ double _phi;
+ //@}
+
+ /**
+ * The variables to get the right boost
+ */
+ //@{
+ /**
+ * LO momenta
+ */
+ Lorentz5Momentum _loMomenta[4];
+ /**
+ * The transfered (virtual boson) momentum
+ */
+ mutable Lorentz5Momentum _pa;
+
+ /**
+ * The incoming quark momentum
+ */
+ mutable Lorentz5Momentum _pb;
+
+ /**
+ * The outgoing quark momentum
+ */
+ Lorentz5Momentum _pc;
+
+ /**
+ * The incoming quark momentum
+ */
+ Lorentz5Momentum _pbother;
+
+ /**
+ * The outgoing quark momentum
+ */
+ Lorentz5Momentum _pcother;
+ //@}
+
+ /**
+ * Electroweak parameters
+ */
+ //@{
+ /**
+ * The square of the Z mass
+ */
+ Energy2 _mz2;
+
+ /**
+ * The square of the W mass
+ */
+ Energy2 _mw2;
+ //@}
+
+ /**
+ * The first hadron
+ */
+ tcBeamPtr _hadron;
+
+ /**
+ * Selects a dynamic or fixed factorization scale
+ */
+ unsigned int scaleOpt_;
+
+ /**
+ * The factorization scale
+ */
+ Energy muF_;
+
+ /**
+ * Prefactor if variable scale used
+ */
+ double scaleFact_;
+
+ /**
+ * Whether to generate the positive, negative or leading order contribution
+ */
+ unsigned int contrib_;
+
+ /**
+ * Power for sampling \f$x_p\f$
+ */
+ double power_;
+
+ /**
+ * Jacobian for \f$x_p\f$ integral
+ */
+ double jac_;
+
+
+};
+
+}
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of MEPP2HiggsVBFPowheg. */
+template <>
+struct BaseClassTrait<Herwig::MEPP2HiggsVBFPowheg,1> {
+ /** Typedef of the first base class of MEPP2HiggsVBFPowheg. */
+ typedef Herwig::MEPP2HiggsVBF NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the MEPP2HiggsVBFPowheg class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::MEPP2HiggsVBFPowheg>
+ : public ClassTraitsBase<Herwig::MEPP2HiggsVBFPowheg> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::MEPP2HiggsVBFPowheg"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * MEPP2HiggsVBFPowheg is implemented. It may also include several, space-separated,
+ * libraries if the class MEPP2HiggsVBFPowheg depends on other classes (base classes
+ * excepted). In this case the listed libraries will be dynamically
+ * linked in the order they are specified.
+ */
+ static string library() { return "HwMEHadron.so HwPowhegMEHadron.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_MEPP2HiggsVBFPowheg_H */
diff --git a/MatrixElement/Powheg/Makefile.am b/MatrixElement/Powheg/Makefile.am
--- a/MatrixElement/Powheg/Makefile.am
+++ b/MatrixElement/Powheg/Makefile.am
@@ -1,14 +1,15 @@
pkglib_LTLIBRARIES = HwPowhegMEHadron.la HwPowhegMELepton.la
HwPowhegMEHadron_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 4:1:0
HwPowhegMEHadron_la_SOURCES = \
MEqq2gZ2ffPowheg.cc MEqq2gZ2ffPowheg.h \
MEqq2W2ffPowheg.cc MEqq2W2ffPowheg.h \
MEPP2HiggsPowheg.cc MEPP2HiggsPowheg.h \
MEPP2WHPowheg.cc MEPP2WHPowheg.h \
MEPP2ZHPowheg.cc MEPP2ZHPowheg.h \
MEPP2VVPowheg.cc MEPP2VVPowheg.h \
-VVKinematics.cc VVKinematics.h
+VVKinematics.cc VVKinematics.h \
+MEPP2HiggsVBFPowheg.h MEPP2HiggsVBFPowheg.cc
HwPowhegMELepton_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
HwPowhegMELepton_la_SOURCES = \
MEee2gZ2qqPowheg.cc MEee2gZ2qqPowheg.h
diff --git a/src/LHC-Powheg.in b/src/LHC-Powheg.in
--- a/src/LHC-Powheg.in
+++ b/src/LHC-Powheg.in
@@ -1,108 +1,110 @@
##################################################
# Example generator based on LHC parameters
# using NLO matrix elements and matching in
# the Powheg formalism
# usage: Herwig++ read LHC.in
##################################################
##################################################
# Technical parameters for this run
##################################################
cd /Herwig/Generators
set LHCGenerator:NumberOfEvents 10000000
set LHCGenerator:RandomNumberGenerator:Seed 31122001
set LHCGenerator:PrintEvent 10
set LHCGenerator:MaxErrors 10000
##################################################
# Need to use an NLO PDF
##################################################
set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
##################################################
# and strong coupling
##################################################
create Herwig::O2AlphaS O2AlphaS
set /Herwig/Generators/LHCGenerator:StandardModelParameters:QCD/RunningAlphaS O2AlphaS
##################################################
# Setup the POWHEG shower
##################################################
cd /Herwig/Shower
set Evolver:HardEmissionMode POWHEG
##################################################
# LHC physics parameters (override defaults here)
##################################################
cd /Herwig/Generators
set LHCGenerator:EventHandler:LuminosityFunction:Energy 7000.0
# Intrinsic pT tune extrapolated to LHC energy
set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV
##################################################
# Matrix Elements for hadron-hadron collisions
# (by default only gamma/Z switched on)
##################################################
cd /Herwig/MatrixElements/
# Drell-Yan Z/gamma
insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff
#
# Drell-Yan W
# insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff
# higgs + W (N.B. if considering all W decay modes useful to set )
# (jet pT cut to zero so no cut on W decay products )
# insert SimpleQCD:MatrixElements[0] PowhegMEPP2WH
# set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
#
# higgs + Z (N.B. if considering all Z decay modes useful to set )
# (jet pT cut to zero so no cut on Z decay products )
# insert SimpleQCD:MatrixElements[0] PowhegMEPP2ZH
# set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
#
# gg/qqbar -> Higgs
# insert SimpleQCD:MatrixElements[0] PowhegMEHiggs
#
# Weak boson pair production: WW / ZZ / WZ / W+Z [WpZ] / W-Z [WmZ]
# insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV
# set PowhegMEPP2VV:Process WpZ
#
+# Higgs production via VBF
+# insert SimpleQCD:MatrixElements[0] PowhegMEPP2HiggsVBF
cd /Herwig/Generators
##################################################
# Useful analysis handlers for hadron-hadron physics
##################################################
# analysis of W/Z events
# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/DrellYan
##################################################
# Useful analysis handlers for HepMC related output
##################################################
# Schematic overview of an event (requires --with-hepmc to be set at configure time
# and the graphviz program 'dot' to produce a plot)
# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot
# A HepMC dump file (requires --with-hepmc to be set at configure time)
# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
# set /Herwig/Analysis/HepMCFile:PrintEvent 100
# set /Herwig/Analysis/HepMCFile:Format GenEvent
# set /Herwig/Analysis/HepMCFile:Units GeV_mm
##################################################
# Save run for later usage with 'Herwig++ run'
##################################################
saverun LHC-Powheg LHCGenerator
##################################################
# uncomment this section for an example batch run
# of two repeats with different parameters
#
# Note that a separate call of 'Herwig run'
# is not required in this case
##################################################
# set LHCGenerator:NumberOfEvents 10
# run LHC-full LHCGenerator
#
# set LHCGenerator:EventHandler:LuminosityFunction:Energy 900.0
# run LHC-initial LHCGenerator
diff --git a/src/defaults/MatrixElements.in b/src/defaults/MatrixElements.in
--- a/src/defaults/MatrixElements.in
+++ b/src/defaults/MatrixElements.in
@@ -1,228 +1,233 @@
##############################################################################
# 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:Coupling /Herwig/Shower/AlphaQCD
# e+e- -> l+l-
create Herwig::MEee2gZ2ll MEee2gZ2ll
newdef MEee2gZ2ll:Allowed Charged
# 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:Coupling /Herwig/Shower/AlphaQCD
############################################################
# 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
# qqbar/gg -> t tbar
create Herwig::MEPP2QQ MEHeavyQuark
###################################
# 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
+
##########################################################
# DIS matrix elements
##########################################################
# neutral current
create Herwig::MENeutralCurrentDIS MEDISNC
-set MEDISNC:Coupling /Herwig/Shower/AlphaQCD
-set MEDISNC:Contribution 0
+newdef MEDISNC:Coupling /Herwig/Shower/AlphaQCD
+newdef MEDISNC:Contribution 0
# charged current
create Herwig::MEChargedCurrentDIS MEDISCC
-set MEDISCC:Coupling /Herwig/Shower/AlphaQCD
-set MEDISCC:Contribution 0
+newdef MEDISCC:Coupling /Herwig/Shower/AlphaQCD
+newdef MEDISCC:Contribution 0
# neutral current (POWHEG)
create Herwig::MENeutralCurrentDIS PowhegMEDISNC
-set PowhegMEDISNC:Coupling /Herwig/Shower/AlphaQCD
-set PowhegMEDISNC:Contribution 1
+newdef PowhegMEDISNC:Coupling /Herwig/Shower/AlphaQCD
+newdef PowhegMEDISNC:Contribution 1
# charged current (POWHEG)
create Herwig::MEChargedCurrentDIS PowhegMEDISCC
-set PowhegMEDISCC:Coupling /Herwig/Shower/AlphaQCD
-set PowhegMEDISCC:Contribution 1
+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
#
# For e+e-
##########################################################
create ThePEG::SubProcessHandler SimpleEE
newdef SimpleEE:PartonExtractor /Herwig/Partons/EEExtractor
##########################################################
# For hadron-hadron
##########################################################
create ThePEG::SubProcessHandler SimpleQCD
newdef SimpleQCD:PartonExtractor /Herwig/Partons/QCDExtractor
##########################################################
# For DIS
##########################################################
create ThePEG::SubProcessHandler SimpleDIS
newdef SimpleDIS:PartonExtractor /Herwig/Partons/DISExtractor

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:20 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243655
Default Alt Text
(115 KB)

Event Timeline