Page MenuHomeHEPForge
Files F7878330

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/DIS/DISBase.cc b/MatrixElement/DIS/DISBase.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/DIS/DISBase.cc
@@ -0,0 +1,1378 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the DISBase class.
+//
+
+#include "DISBase.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Reference.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
+#include "Herwig++/Utilities/Maths.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "ThePEG/PDT/StandardMatchers.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Repository/CurrentGenerator.h"
+#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
+#include "Herwig++/PDT/StandardMatchers.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include <numeric>
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+namespace {
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+void debuggingMatrixElement(bool BGF,const Lorentz5Momentum & pin,
+ const Lorentz5Momentum & p1,
+ const Lorentz5Momentum & p2,
+ tcPDPtr gluon,
+ const Lorentz5Momentum & pl1,
+ const Lorentz5Momentum & pl2,
+ const Lorentz5Momentum & pq1,
+ const Lorentz5Momentum & pq2,
+ tcPDPtr lepton1,tcPDPtr lepton2,
+ tcPDPtr quark1 ,tcPDPtr quark2,
+ Energy2 Q2,double phi, double x2, double x3,
+ double xperp, double zp, double xp,
+ const vector<double> & azicoeff,
+ bool normalize) {
+ tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>
+ (CurrentGenerator::current().standardModel());
+ assert(hwsm);
+ vector<AbstractFFVVertexPtr> weakVertex;
+ vector<PDPtr> bosons;
+ AbstractFFVVertexPtr strongVertex = hwsm->vertexFFG();
+ if(lepton1->id()==lepton2->id()) {
+ weakVertex.push_back(hwsm->vertexFFZ());
+ bosons.push_back(hwsm->getParticleData(ParticleID::Z0));
+ weakVertex.push_back(hwsm->vertexFFP());
+ bosons.push_back(hwsm->getParticleData(ParticleID::gamma));
+ }
+ else {
+ weakVertex.push_back(hwsm->vertexFFW());
+ bosons.push_back(hwsm->getParticleData(ParticleID::Wplus));
+ }
+ if(!BGF) {
+ SpinorWaveFunction l1,q1,qp1;
+ SpinorBarWaveFunction l2,q2,qp2;
+ VectorWaveFunction gl(p2,gluon,outgoing);
+ if(lepton1->id()>0) {
+ l1 = SpinorWaveFunction (pl1,lepton1,incoming);
+ l2 = SpinorBarWaveFunction(pl2,lepton2,outgoing);
+ }
+ else {
+ l1 = SpinorWaveFunction (pl2,lepton2,outgoing);
+ l2 = SpinorBarWaveFunction(pl1,lepton1,incoming);
+ }
+ if(quark1->id()>0) {
+ q1 = SpinorWaveFunction (pq1,quark1,incoming);
+ q2 = SpinorBarWaveFunction(pq2,quark2,outgoing);
+ qp1 = SpinorWaveFunction (pin,quark1,incoming);
+ qp2 = SpinorBarWaveFunction(p1 ,quark2,outgoing);
+ }
+ else {
+ q1 = SpinorWaveFunction (pq2,quark2,outgoing);
+ q2 = SpinorBarWaveFunction(pq1,quark1,incoming);
+ qp1 = SpinorWaveFunction (p1 ,quark2,outgoing);
+ qp2 = SpinorBarWaveFunction(pin,quark1,incoming);
+ }
+ double lome(0.),realme(0.);
+ for(unsigned int lhel1=0;lhel1<2;++lhel1) {
+ l1.reset(lhel1);
+ for(unsigned int lhel2=0;lhel2<2;++lhel2) {
+ l2.reset(lhel2);
+ for(unsigned int qhel1=0;qhel1<2;++qhel1) {
+ q1.reset(qhel1);
+ qp1.reset(qhel1);
+ for(unsigned int qhel2=0;qhel2<2;++qhel2) {
+ q2.reset(qhel2);
+ qp2.reset(qhel2);
+ // leading order matrix element
+ Complex diagLO(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ diagLO += weakVertex[ix]->evaluate(Q2,q1,q2,inter);
+ }
+ lome += norm(diagLO);
+ // real emission matrix element
+ for(unsigned int ghel=0;ghel<2;++ghel) {
+ gl.reset(2*ghel);
+ Complex diagReal(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ SpinorWaveFunction off1 =
+ strongVertex->evaluate(Q2,5,qp1.particle(),qp1,gl);
+ Complex diag1 = weakVertex[ix]->evaluate(Q2,off1,qp2,inter);
+ SpinorBarWaveFunction off2 =
+ strongVertex->evaluate(Q2,5,qp2.particle(),qp2,gl);
+ Complex diag2 = weakVertex[ix]->evaluate(Q2,qp1,off2,inter);
+ diagReal += diag1+diag2;
+ }
+ realme += norm(diagReal);
+ }
+ }
+ }
+ }
+ }
+ double test1 = realme/lome/hwsm->alphaS(Q2)*Q2*UnitRemoval::InvE2;
+ double cphi(cos(phi));
+ double test2;
+ if(normalize) {
+ test2 = 8.*Constants::pi/(1.-xp)/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi))*
+ (1.+sqr(xp)*(sqr(x2)+1.5*sqr(xperp)));
+ }
+ else {
+ test2 = 8.*Constants::pi/(1.-xp)/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi));
+ }
+ cerr << "testing RATIO A " << test1/test2 << "\n";
+ }
+ else {
+ SpinorWaveFunction l1,q1,qp1;
+ SpinorBarWaveFunction l2,q2,qp2;
+ VectorWaveFunction gl(pin,gluon,incoming);
+ if(lepton1->id()>0) {
+ l1 = SpinorWaveFunction (pl1,lepton1,incoming);
+ l2 = SpinorBarWaveFunction(pl2,lepton2,outgoing);
+ }
+ else {
+ l1 = SpinorWaveFunction (pl2,lepton2,outgoing);
+ l2 = SpinorBarWaveFunction(pl1,lepton1,incoming);
+ }
+ if(quark1->id()>0) {
+ q1 = SpinorWaveFunction (pq1,quark1 ,incoming);
+ q2 = SpinorBarWaveFunction(pq2,quark2 ,outgoing);
+ qp2 = SpinorBarWaveFunction(p1 ,quark2 ,outgoing);
+ qp1 = SpinorWaveFunction (p2 ,quark1->CC(),outgoing);
+ }
+ else {
+ q1 = SpinorWaveFunction (pq2,quark2 ,outgoing);
+ q2 = SpinorBarWaveFunction(pq1,quark1 ,incoming);
+ qp2 = SpinorBarWaveFunction(p2 ,quark1->CC(),outgoing);
+ qp1 = SpinorWaveFunction (p1 ,quark2 ,outgoing);
+ }
+ double lome(0.),realme(0.);
+ for(unsigned int lhel1=0;lhel1<2;++lhel1) {
+ l1.reset(lhel1);
+ for(unsigned int lhel2=0;lhel2<2;++lhel2) {
+ l2.reset(lhel2);
+ for(unsigned int qhel1=0;qhel1<2;++qhel1) {
+ q1.reset(qhel1);
+ qp1.reset(qhel1);
+ for(unsigned int qhel2=0;qhel2<2;++qhel2) {
+ q2.reset(qhel2);
+ qp2.reset(qhel2);
+ // leading order matrix element
+ Complex diagLO(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ diagLO += weakVertex[ix]->evaluate(Q2,q1,q2,inter);
+ }
+ lome += norm(diagLO);
+ // real emission matrix element
+ for(unsigned int ghel=0;ghel<2;++ghel) {
+ gl.reset(2*ghel);
+ Complex diagReal(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ SpinorWaveFunction off1 =
+ strongVertex->evaluate(Q2,5,qp1.particle(),qp1,gl);
+ Complex diag1 = weakVertex[ix]->evaluate(Q2,off1,qp2,inter);
+ SpinorBarWaveFunction off2 =
+ strongVertex->evaluate(Q2,5,qp2.particle(),qp2,gl);
+ Complex diag2 = weakVertex[ix]->evaluate(Q2,qp1,off2,inter);
+ diagReal += diag1+diag2;
+ }
+ realme += norm(diagReal);
+ }
+ }
+ }
+ }
+ }
+ double test1 = realme/lome/hwsm->alphaS(Q2)*Q2*UnitRemoval::InvE2;
+ double cphi(cos(phi));
+ double test2;
+ if(normalize) {
+ test2 = 8.*Constants::pi/zp/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi))*
+ sqr(xp)*(sqr(x3)+sqr(x2)+3.*sqr(xperp));
+ }
+ else {
+ test2 = 8.*Constants::pi/zp/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi));
+ }
+ cerr << "testing RATIO B " << test1/test2 << "\n";
+ }
+}
+
+}
+
+DISBase::DISBase() : initial_(6.), final_(3.),
+ procProb_(0.35),
+ comptonInt_(0.), bgfInt_(0.),
+ comptonWeight_(50.), BGFWeight_(150.),
+ pTmin_(0.1*GeV),
+ scaleOpt_(1), muF_(100.*GeV), scaleFact_(1.),
+ contrib_(0), power_(0.1)
+{}
+
+DISBase::~DISBase() {}
+
+void DISBase::persistentOutput(PersistentOStream & os) const {
+ os << comptonInt_ << bgfInt_ << procProb_ << initial_ << final_ << alpha_
+ << ounit(pTmin_,GeV) << comptonWeight_ << BGFWeight_ << gluon_
+ << ounit(muF_,GeV) << scaleFact_ << scaleOpt_ << contrib_<< power_;
+}
+
+void DISBase::persistentInput(PersistentIStream & is, int) {
+ is >> comptonInt_ >> bgfInt_ >> procProb_ >> initial_ >> final_ >> alpha_
+ >> iunit(pTmin_,GeV) >> comptonWeight_ >> BGFWeight_ >> gluon_
+ >> iunit(muF_,GeV) >> scaleFact_ >> scaleOpt_ >> contrib_ >> power_;
+}
+
+AbstractClassDescription<DISBase> DISBase::initDISBase;
+// Definition of the static class description member.
+
+void DISBase::Init() {
+
+ static ClassDocumentation<DISBase> documentation
+ ("The DISBase class provides the base class for the "
+ "implementation of DIS type processes including the "
+ "hard corrections in either the old-fashioned matrix "
+ "element correction of POWHEG approaches");
+
+ static Parameter<DISBase,double> interfaceProcessProbability
+ ("ProcessProbability",
+ "The probabilty of the QCD compton process for the process selection",
+ &DISBase::procProb_, 0.3, 0.0, 1.,
+ false, false, Interface::limited);
+
+ static Reference<DISBase,ShowerAlpha> interfaceCoupling
+ ("Coupling",
+ "Pointer to the object to calculate the coupling for the correction",
+ &DISBase::alpha_, false, false, true, false, false);
+
+ static Parameter<DISBase,Energy> interfacepTMin
+ ("pTMin",
+ "The minimum pT",
+ &DISBase::pTmin_, GeV, 1.*GeV, 0.0*GeV, 10.0*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceComptonWeight
+ ("ComptonWeight",
+ "Weight for the overestimate ofthe compton channel",
+ &DISBase::comptonWeight_, 50.0, 0.0, 100.0,
+ false, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceBGFWeight
+ ("BGFWeight",
+ "Weight for the overestimate of the BGF channel",
+ &DISBase::BGFWeight_, 100.0, 0.0, 1000.0,
+ false, false, Interface::limited);
+
+ static Switch<DISBase,unsigned int> interfaceContribution
+ ("Contribution",
+ "Which contributions to the cross section to include",
+ &DISBase::contrib_, 0, false, false);
+ static SwitchOption interfaceContributionLeadingOrder
+ (interfaceContribution,
+ "LeadingOrder",
+ "Just generate the leading order cross section",
+ 0);
+ static SwitchOption interfaceContributionPositiveNLO
+ (interfaceContribution,
+ "PositiveNLO",
+ "Generate the positive contribution to the full NLO cross section",
+ 1);
+ static SwitchOption interfaceContributionNegativeNLO
+ (interfaceContribution,
+ "NegativeNLO",
+ "Generate the negative contribution to the full NLO cross section",
+ 2);
+
+ static Switch<DISBase,unsigned int> interfaceScaleOption
+ ("ScaleOption",
+ "Option for the choice of factorization (and renormalization) scale",
+ &DISBase::scaleOpt_, 1, false, false);
+ static SwitchOption interfaceDynamic
+ (interfaceScaleOption,
+ "Dynamic",
+ "Dynamic factorization scale equal to the current sqrt(sHat())",
+ 1);
+ static SwitchOption interfaceFixed
+ (interfaceScaleOption,
+ "Fixed",
+ "Use a fixed factorization scale set with FactorizationScaleValue",
+ 2);
+
+ static Parameter<DISBase,Energy> interfaceFactorizationScale
+ ("FactorizationScale",
+ "Value to use in the event of a fixed factorization scale",
+ &DISBase::muF_, GeV, 100.0*GeV, 1.0*GeV, 500.0*GeV,
+ true, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceScaleFactor
+ ("ScaleFactor",
+ "The factor used before Q2 if using a running scale",
+ &DISBase::scaleFact_, 1.0, 0.0, 10.0,
+ false, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceSamplingPower
+ ("SamplingPower",
+ "Power for the sampling of xp",
+ &DISBase::power_, 0.6, 0.0, 1.,
+ false, false, Interface::limited);
+}
+
+void DISBase::doinit() {
+ HwMEBase::doinit();
+ // integrals of me over phase space
+ double r5=sqrt(5.),darg((r5-1.)/(r5+1.)),ath(0.5*log((1.+1./r5)/(1.-1./r5)));
+ comptonInt_ = 2.*(-21./20.-6./(5.*r5)*ath+sqr(Constants::pi)/3.
+ -2.*Math::ReLi2(1.-darg)-2.*Math::ReLi2(1.-1./darg));
+ bgfInt_ = 121./9.-56./r5*ath;
+ // extract the gluon ParticleData objects
+ gluon_ = getParticleData(ParticleID::g);
+}
+
+void DISBase::initializeMECorrection(ShowerTreePtr tree, double & initial,
+ double & final) {
+ initial = initial_;
+ final = final_;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ partons_[0] = cit->first->progenitor()->dataPtr();
+ pq_[0] = cit->first->progenitor()->momentum();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ leptons_[0] = cit->first->progenitor()->dataPtr();
+ pl_[0] = cit->first->progenitor()->momentum();
+ }
+ }
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ partons_[1] = cit->first->progenitor()->dataPtr();
+ pq_[1] = cit->first->progenitor()->momentum();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ leptons_[1] = cit->first->progenitor()->dataPtr();
+ pl_[1] = cit->first->progenitor()->momentum();
+ }
+ }
+ // extract the born variables
+ q_ =pl_[0]-pl_[1];
+ q2_ = -q_.m2();
+ double yB = (q_*pq_[0])/(pl_[0]*pq_[0]);
+ l_ = 2./yB-1.;
+ // calculate the A coefficient for the correlations
+ acoeff_ = A(leptons_[0],leptons_[1],
+ partons_[0],partons_[1],q2_);
+}
+
+void DISBase::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
+ static const double eps=1e-6;
+ // find the incoming and outgoing quarks and leptons
+ ShowerParticlePtr quark[2],lepton[2];
+ PPtr hadron;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ hadron = cit->first->original()->parents()[0];
+ quark [0] = cit->first->progenitor();
+ beam_ = cit->first->beam();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[0] = cit->first->progenitor();
+ }
+ }
+ pdf_ = beam_->pdf();
+ assert(beam_&&pdf_&&quark[0]&&lepton[0]);
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ quark [1] = cit->first->progenitor();
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[1] = cit->first->progenitor();
+ }
+ }
+ // momentum fraction
+ assert(quark[1]&&lepton[1]);
+ xB_ = quark[0]->x();
+ // calculate the matrix element
+ vector<double> azicoeff;
+ // select the type of process
+ bool BGF = UseRandom::rnd()>procProb_;
+ double xp,zp,wgt,x1,x2,x3,xperp;
+ // generate a QCD compton process
+ if(!BGF) {
+ wgt = generateComptonPoint(xp,zp);
+ if(xp<eps) return;
+ // common pieces
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
+ wgt *= 2./3./Constants::pi*alpha_->value(scale)/procProb_;
+ // PDF piece
+ wgt *= pdf_->xfx(beam_,quark[0]->dataPtr(),scale,xB_/xp)/
+ pdf_->xfx(beam_,quark[0]->dataPtr(),q2_ ,xB_);
+ // other bits
+ xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ x1 = -1./xp;
+ x2 = 1.-(1.-zp)/xp;
+ x3 = 2.+x1-x2;
+ // matrix element pieces
+ azicoeff = ComptonME(xp,x2,xperp,true);
+ }
+ // generate a BGF process
+ else {
+ wgt = generateBGFPoint(xp,zp);
+ if(xp<eps) return;
+ // common pieces
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1);
+ wgt *= 0.25/Constants::pi*alpha_->value(scale)/(1.-procProb_);
+ // PDF piece
+ wgt *= pdf_->xfx(beam_,gluon_ ,scale,xB_/xp)/
+ pdf_->xfx(beam_,quark[0]->dataPtr(),q2_ ,xB_);
+ // other bits
+ xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ x1 = -1./xp;
+ x2 = 1.-(1.-zp)/xp;
+ x3 = 2.+x1-x2;
+ // matrix element pieces
+ azicoeff = BGFME(xp,x2,x3,xperp,true);
+ }
+ // compute the azimuthal average of the weight
+ wgt *= (azicoeff[0]+0.5*azicoeff[2]);
+ // decide whether or not to accept the weight
+ if(UseRandom::rnd()>wgt) return;
+ // if generate generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi));
+ phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
+ << "::applyHardMatrixElementCorrection() to"
+ << " generate phi" << Exception::eventerror;
+ // construct lorentz transform from lab to breit frame
+ Lorentz5Momentum phadron = hadron->momentum();
+ phadron.setMass(0.*GeV);
+ phadron.rescaleEnergy();
+ Lorentz5Momentum pcmf = phadron+0.5/xB_*q_;
+ pcmf.rescaleMass();
+ LorentzRotation rot(-pcmf.boostVector());
+ Lorentz5Momentum pbeam = rot*phadron;
+ Axis axis(pbeam.vect().unit());
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ rot.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ Lorentz5Momentum pl = rot*pl_[0];
+ rot.rotateZ(-atan2(pl.y(),pl.x()));
+ pl_[0] *= rot;
+ pl_[1] *= rot;
+ pq_[0] *= rot;
+ pq_[1] *= rot;
+ // compute the new incoming and outgoing momenta
+ Energy Q(sqrt(q2_));
+ Lorentz5Momentum p1 = Lorentz5Momentum( 0.5*Q*xperp*cos(phi), 0.5*Q*xperp*sin(phi),
+ -0.5*Q*x2,0.*GeV,0.*GeV);
+ p1.rescaleEnergy();
+ Lorentz5Momentum p2 = Lorentz5Momentum(-0.5*Q*xperp*cos(phi),-0.5*Q*xperp*sin(phi),
+ -0.5*Q*x3,0.*GeV,0.*GeV);
+ p2.rescaleEnergy();
+ Lorentz5Momentum pin(0.*GeV,0.*GeV,-0.5*x1*Q,-0.5*x1*Q,0.*GeV);
+// debuggingMatrixElement(BGF,pin,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
+// lepton[0]->dataPtr(),lepton[1]->dataPtr(),
+// quark [0]->dataPtr(),quark [1]->dataPtr(),
+// q2_,phi,x2,x3,xperp,zp,xp,azicoeff,true);
+ // we need the Lorentz transform back to the lab
+ rot.invert();
+ // transform the momenta to lab frame
+ pin *= rot;
+ p1 *= rot;
+ p2 *= rot;
+ // test to ensure outgoing particles can be put on-shell
+ if(!BGF) {
+ if(p1.e()<quark[1]->dataPtr()->constituentMass()) return;
+ if(p2.e()<gluon_ ->constituentMass()) return;
+ }
+ else {
+ if(p1.e()<quark[1]->dataPtr() ->constituentMass()) return;
+ if(p2.e()<quark[0]->dataPtr()->CC()->constituentMass()) return;
+ }
+ // create the new particles and add to ShowerTree
+ bool isquark = quark[0]->colourLine();
+ if(!BGF) {
+ PPtr newin = new_ptr(Particle(*quark[0]));
+ newin->set5Momentum(pin);
+ PPtr newg = gluon_ ->produceParticle(p2 );
+ PPtr newout = quark[1]->dataPtr()->produceParticle(p1 );
+ ColinePtr col=isquark ?
+ quark[0]->colourLine() : quark[0]->antiColourLine();
+ ColinePtr newline=new_ptr(ColourLine());
+ // final-state emission
+ if(xp>zp) {
+ col->removeColoured(newout,!isquark);
+ col->addColoured(newin,!isquark);
+ col->addColoured(newg,!isquark);
+ newline->addColoured(newg,isquark);
+ newline->addColoured(newout,!isquark);
+ }
+ // initial-state emission
+ else {
+ col->removeColoured(newin ,!isquark);
+ col->addColoured(newout,!isquark);
+ col->addColoured(newg,isquark);
+ newline->addColoured(newg,!isquark);
+ newline->addColoured(newin,!isquark);
+ }
+ PPtr orig;
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[0]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newin);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
+ cit->first->progenitor(sp);
+ tree->incomingLines()[cit->first]=sp;
+ sp->x(xB_/xp);
+ cit->first->perturbative(xp>zp);
+ if(xp<=zp) orig=cit->first->original();
+ }
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[1]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newout);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
+ cit->first->progenitor(sp);
+ tree->outgoingLines()[cit->first]=sp;
+ cit->first->perturbative(xp<=zp);
+ if(xp>zp) orig=cit->first->original();
+ }
+ assert(orig);
+ // add the gluon
+ ShowerParticlePtr sg=new_ptr(ShowerParticle(*newg,1,true));
+ ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
+ gluon->perturbative(false);
+ tree->outgoingLines().insert(make_pair(gluon,sg));
+ tree->hardMatrixElementCorrection(true);
+ }
+ else {
+ PPtr newin = gluon_ ->produceParticle(pin);
+ PPtr newqbar = quark[0]->dataPtr()->CC()->produceParticle(p2 );
+ PPtr newout = quark[1]->dataPtr() ->produceParticle(p1 );
+ ColinePtr col=isquark ? quark[0]->colourLine() : quark[0]->antiColourLine();
+ ColinePtr newline=new_ptr(ColourLine());
+ col ->addColoured(newin ,!isquark);
+ newline->addColoured(newin , isquark);
+ col ->addColoured(newout ,!isquark);
+ newline->addColoured(newqbar, isquark);
+ PPtr orig;
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[0]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newin);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
+ cit->first->progenitor(sp);
+ tree->incomingLines()[cit->first]=sp;
+ sp->x(xB_/xp);
+ cit->first->perturbative(false);
+ orig=cit->first->original();
+ }
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[1]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newout);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
+ cit->first->progenitor(sp);
+ tree->outgoingLines()[cit->first]=sp;
+ cit->first->perturbative(true);
+ }
+ assert(orig);
+ // add the (anti)quark
+ ShowerParticlePtr sqbar=new_ptr(ShowerParticle(*newqbar,1,true));
+ ShowerProgenitorPtr qbar=new_ptr(ShowerProgenitor(orig,newqbar,sqbar));
+ qbar->perturbative(false);
+ tree->outgoingLines().insert(make_pair(qbar,sqbar));
+ tree->hardMatrixElementCorrection(true);
+ }
+}
+
+bool DISBase::softMatrixElementVeto(ShowerProgenitorPtr initial,
+ ShowerParticlePtr parent, Branching br) {
+ bool veto = !UseRandom::rndbool(parent->isFinalState() ? 1./final_ : 1./initial_);
+ // check if me correction should be applied
+ long id[2]={initial->id(),parent->id()};
+ if(id[0]!=id[1]||id[1]==ParticleID::g) return veto;
+ // get the pT
+ Energy pT=br.kinematics->pT();
+ // check if hardest so far
+ if(pT<initial->highestpT()) return veto;
+ double kappa(sqr(br.kinematics->scale())/q2_),z(br.kinematics->z());
+ double zk((1.-z)*kappa);
+ // final-state
+ double wgt(0.);
+ if(parent->isFinalState()) {
+ double zp=z,xp=1./(1.+z*zk);
+ double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ double x2 = 1.-(1.-zp)/xp;
+ vector<double> azicoeff = ComptonME(xp,x2,xperp,false);
+ wgt = (azicoeff[0]+0.5*azicoeff[2])*xp/(1.+sqr(z))/final_;
+ if(wgt<.0||wgt>1.) {
+ ostringstream wstring;
+ wstring << "Soft ME correction weight too large or "
+ << "negative for FSR in DISBase::"
+ << "softMatrixElementVeto() soft weight "
+ << " xp = " << xp << " zp = " << zp
+ << " weight = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ else {
+ double xp = 2.*z/(1.+zk+sqrt(sqr(1.+zk)-4.*z*zk));
+ double zp = 0.5* (1.-zk+sqrt(sqr(1.+zk)-4.*z*zk));
+ double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ double x1 = -1./xp, x2 = 1.-(1.-zp)/xp, x3 = 2.+x1-x2;
+ // compton
+ if(br.ids[0]!=ParticleID::g) {
+ vector<double> azicoeff = ComptonME(xp,x2,xperp,false);
+ wgt = (azicoeff[0]+0.5*azicoeff[2])*xp*(1.-z)/(1.-xp)/(1.+sqr(z))/
+ (1.-zp+xp-2.*xp*(1.-zp));
+ }
+ // BGF
+ else {
+ vector<double> azicoeff = BGFME(xp,x2,x3,xperp,true);
+ wgt = (azicoeff[0]+0.5*azicoeff[2])*xp/(1.-zp+xp-2.*xp*(1.-zp))/(sqr(z)+sqr(1.-z));
+ }
+ wgt /=initial_;
+ if(wgt<.0||wgt>1.) {
+ ostringstream wstring;
+ wstring << "Soft ME correction weight too large or "
+ << "negative for ISR in DISBase::"
+ << "softMatrixElementVeto() soft weight "
+ << " xp = " << xp << " zp = " << zp
+ << " weight = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ // if not vetoed
+ if(UseRandom::rndbool(wgt)) {
+ initial->highestpT(pT);
+ return false;
+ }
+ // otherwise
+ parent->setEvolutionScale(br.kinematics->scale());
+ return true;
+}
+
+double DISBase::generateComptonPoint(double &xp, double & zp) {
+ static const double maxwgt = 1.;
+ double wgt;
+ do {
+ xp = UseRandom::rnd();
+ double zpmin = xp, zpmax = 1./(1.+xp*(1.-xp));
+ zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+ wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+ if(UseRandom::rndbool()) swap(xp,zp);
+ double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp,x2=1.-(1.-zp)/xp;
+ wgt *= 2.*(1.+sqr(xp)*(sqr(x2)+1.5*xperp2))/(1.-xp)/(1.-zp);
+ if(wgt>maxwgt) {
+ ostringstream wstring;
+ wstring << "DISBase::generateComptonPoint "
+ << "Weight greater than maximum "
+ << "wgt = " << wgt << " maxwgt = 1\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(wgt<UseRandom::rnd()*maxwgt);
+ return comptonInt_;
+}
+
+double DISBase::generateBGFPoint(double &xp, double & zp) {
+ static const double maxwgt = 25.;
+ double wgt;
+ do {
+ xp = UseRandom::rnd();
+ double zpmax = 1./(1.+xp*(1.-xp)), zpmin = 1.-zpmax;
+ zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+ wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
+ wgt *= sqr(xp)/(1.-zp)*(sqr(x3)+sqr(x2)+3.*xperp2);
+ if(wgt>maxwgt) {
+ ostringstream wstring;
+ wstring << "DISBase::generateBGFPoint "
+ << "Weight greater than maximum "
+ << "wgt = " << wgt << " maxwgt = 1\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(wgt<UseRandom::rnd()*maxwgt);
+ return bgfInt_;
+// static const double maxwgt = 2.,npow=0.34,ac=1.0;
+// double wgt;
+// do {
+// double rho = UseRandom::rnd();
+// xp = 1.-pow(rho,1./(1.-npow));
+// wgt = (sqr(xp)+ac+sqr(1.-xp));
+// if(wgt>1.+ac) cerr << "testing violates BGF maxA " << wgt << "\n";
+// }
+// while(wgt<UseRandom::rnd()*(1.+ac));
+// double xpwgt = -((6.-5.*npow+sqr(npow))*ac-3.*npow+sqr(npow)+4)
+// /(sqr(npow)*(npow-6.)+11.*npow-6.);
+// xpwgt *= pow(1.-xp,npow)/wgt;
+// double xp2(sqr(xp)),lxp(log(xp)),xp4(sqr(xp2)),lxp1(log(1.-xp));
+// double zpwgt = (2.*xp4*(lxp+lxp1-3.)+4.*xp*xp2*(3.-lxp-lxp1)
+// +xp2*(-13.+lxp+lxp1)+xp*(+7.+lxp+lxp1)-lxp-lxp1-1.)/(1.+xp-xp2);
+// do {
+// double zpmax = 1./(1.+xp*(1.-xp)), zpmin = 1.-zpmax;
+// zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+// wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+// double x1 = -1./xp;
+// double x2 = 1.-(1.-zp)/xp;
+// double x3 = 2.+x1-x2;
+// double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
+// wgt *= sqr(xp)/(1.-zp)*(sqr(x3)+sqr(x2)+3.*xperp2);
+// if(wgt>maxwgt*zpwgt) cerr << "testing violates BGF maxB " << wgt/xpwgt << "\n";
+// }
+// while(wgt<UseRandom::rnd()*maxwgt);
+// return zpwgt*xpwgt;
+}
+
+vector<double> DISBase::ComptonME(double xp, double x2, double xperp,
+ bool norm) {
+ vector<double> output(3,0.);
+ double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
+ double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
+ double root = sqrt(sqr(l_)-1.);
+ output[0] = sqr(cos2)+acoeff_*cos2*l_+sqr(l_);
+ output[1] = -acoeff_*cos2*root*sin2-2.*l_*root*sin2;
+ output[2] = sqr(root)*sqr(sin2);
+ double lo(1+acoeff_*l_+sqr(l_));
+ double denom = norm ? 1.+sqr(xp)*(sqr(x2)+1.5*sqr(xperp)) : 1.;
+ double fact = sqr(xp)*(sqr(x2)+sqr(xperp))/lo;
+ for(unsigned int ix=0;ix<output.size();++ix)
+ output[ix] = ((ix==0 ? 1. : 0.) +fact*output[ix])/denom;
+ return output;
+}
+
+vector<double> DISBase::BGFME(double xp, double x2, double x3,
+ double xperp, bool norm) {
+ vector<double> output(3,0.);
+ double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
+ double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
+ double fact2 = sqr(xp)*(sqr(x2)+sqr(xperp));
+ double cos3 = x3 /sqrt(sqr(x3)+sqr(xperp));
+ double sin3 = xperp/sqrt(sqr(x3)+sqr(xperp));
+ double fact3 = sqr(xp)*(sqr(x3)+sqr(xperp));
+ double root = sqrt(sqr(l_)-1.);
+ output[0] = fact2*(sqr(cos2)+acoeff_*cos2*l_+sqr(l_)) +
+ fact3*(sqr(cos3)-acoeff_*cos3*l_+sqr(l_));
+ output[1] = - fact2*(acoeff_*cos2*root*sin2+2.*l_*root*sin2)
+ - fact3*(acoeff_*cos3*root*sin3-2.*l_*root*sin3);
+ output[2] = fact2*(sqr(root)*sqr(sin2)) +
+ fact3*(sqr(root)*sqr(sin3));
+ double lo(1+acoeff_*l_+sqr(l_));
+ double denom = norm ? sqr(xp)*(sqr(x3)+sqr(x2)+3.*sqr(xperp))*lo : lo;
+ for(unsigned int ix=0;ix<output.size();++ix) output[ix] /= denom;
+ return output;
+}
+
+HardTreePtr DISBase::generateHardest(ShowerTreePtr tree,
+ vector<ShowerInteraction::Type> inter) {
+ bool found = false;
+ // check if generating QCD radiation
+ for(unsigned int ix=0;ix<inter.size();++ix) {
+ found |= inter[ix]==ShowerInteraction::QCD;
+ }
+ if(!found) return HardTreePtr();
+ ShowerParticlePtr quark[2],lepton[2];
+ PPtr hadron;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ hadron = cit->first->original()->parents()[0];
+ quark [0] = cit->first->progenitor();
+ beam_ = cit->first->beam();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[0] = cit->first->progenitor();
+ leptons_[0] = lepton[0]->dataPtr();
+ }
+ }
+ pdf_=beam_->pdf();
+ assert(beam_&&pdf_&&quark[0]&&lepton[0]);
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ quark [1] = cit->first->progenitor();
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[1] = cit->first->progenitor();
+ leptons_[1] = lepton[1]->dataPtr();
+ }
+ }
+ assert(quark[1]&&lepton[1]);
+ // Particle data objects
+ for(unsigned int ix=0;ix<2;++ix) partons_[ix] = quark[ix]->dataPtr();
+ // extract the born variables
+ q_ =lepton[0]->momentum()-lepton[1]->momentum();
+ q2_ = -q_.m2();
+ xB_ = quark[0]->x();
+ double yB =
+ ( q_*quark[0]->momentum())/
+ (lepton[0]->momentum()*quark[0]->momentum());
+ l_ = 2./yB-1.;
+ // construct lorentz transform from lab to breit frame
+ Lorentz5Momentum phadron = hadron->momentum();
+ phadron.setMass(0.*GeV);
+ phadron.rescaleRho();
+ Lorentz5Momentum pb = quark[0]->momentum();
+ Lorentz5Momentum pbasis = phadron;
+ Axis axis(q_.vect().unit());
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ LorentzRotation rot_ = LorentzRotation();
+ if(axis.perp2()>1e-20) {
+ rot_.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ rot_.rotateX(Constants::pi);
+ }
+ if(abs(1.-q_.e()/q_.vect().mag())>1e-6) rot_.boostZ( q_.e()/q_.vect().mag());
+ pb *= rot_;
+ if(pb.perp2()/GeV2>1e-20) {
+ Boost trans = -1./pb.e()*pb.vect();
+ trans.setZ(0.);
+ rot_.boost(trans);
+ }
+ Lorentz5Momentum pl = rot_*lepton[0]->momentum();
+ rot_.rotateZ(-atan2(pl.y(),pl.x()));
+ // momenta of the particles
+ pl_[0]=rot_*lepton[0]->momentum();
+ pl_[1]=rot_*lepton[1]->momentum();
+ pq_[0]=rot_* quark[0]->momentum();
+ pq_[1]=rot_* quark[1]->momentum();
+ q_ *= rot_;
+ // coefficient for the matrix elements
+ acoeff_ = A(lepton[0]->dataPtr(),lepton[1]->dataPtr(),
+ quark [0]->dataPtr(),quark [1]->dataPtr(),q2_);
+ // generate a compton point
+ generateCompton();
+ generateBGF();
+ // no valid emission, return
+ if(pTCompton_<ZERO&&pTBGF_<ZERO) return HardTreePtr();
+ // type of emission, pick highest pT
+ bool isCompton=pTCompton_>pTBGF_;
+// // find the sudakov for the branching
+// SudakovPtr sudakov;
+// // ISR
+// if(ComptonISFS_||!isCompton) {
+// BranchingList branchings=evolver()->splittingGenerator()->initialStateBranchings();
+// long index = abs(partons_[0]->id());
+// IdList br(3);
+// if(isCompton) {
+// br[0] = index;
+// br[1] = index;
+// br[2] = ParticleID::g;
+// }
+// else {
+// br[0] = ParticleID::g;
+// br[1] = abs(partons_[0]->id());
+// br[2] = -abs(partons_[0]->id());
+// }
+// for(BranchingList::const_iterator cit = branchings.lower_bound(index);
+// cit != branchings.upper_bound(index); ++cit ) {
+// IdList ids = cit->second.second;
+// if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
+// sudakov=cit->second.first;
+// break;
+// }
+// }
+// }
+// // FSR
+// else {
+// BranchingList branchings =
+// evolver()->splittingGenerator()->finalStateBranchings();
+// long index=abs(partons_[1]->id());
+// for(BranchingList::const_iterator cit = branchings.lower_bound(index);
+// cit != branchings.upper_bound(index); ++cit ) {
+// IdList ids = cit->second.second;
+// if(ids[0]==index&&ids[1]==index&&ids[2]==ParticleID::g) {
+// sudakov = cit->second.first;
+// break;
+// }
+// }
+// }
+// if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
+// << "DISBase::generateHardest()"
+// << Exception::runerror;
+ // add the leptons
+ vector<HardBranchingPtr> spaceBranchings,allBranchings;
+ spaceBranchings.push_back(new_ptr(HardBranching(lepton[0],SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ allBranchings.push_back(spaceBranchings.back());
+ allBranchings.push_back(new_ptr(HardBranching(lepton[1],SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ // compton hardest
+ if(isCompton) {
+ rot_.invert();
+ for(unsigned int ix=0;ix<ComptonMomenta_.size();++ix) {
+ ComptonMomenta_[ix].transform(rot_);
+ }
+ ShowerParticlePtr newqout (new_ptr(ShowerParticle(partons_[1],true)));
+ newqout->set5Momentum(ComptonMomenta_[1]);
+ ShowerParticlePtr newg(new_ptr(ShowerParticle(gluon_,true)));
+ newg->set5Momentum(ComptonMomenta_[2]);
+ ShowerParticlePtr newqin (new_ptr(ShowerParticle(partons_[0],false )));
+ newqin->set5Momentum(ComptonMomenta_[0]);
+ if(ComptonISFS_) {
+ ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[0],false)));
+ newspace->set5Momentum(ComptonMomenta_[0]-ComptonMomenta_[2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),
+ spaceBranch,
+ HardBranching::Incoming)));
+ spaceBranch->addChild(offBranch);
+ HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),spaceBranch,
+ HardBranching::Outgoing)));
+ spaceBranch->addChild(g);
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ allBranchings.push_back(outBranch);
+ ColinePtr newin(new_ptr(ColourLine())),newout(new_ptr(ColourLine()));
+ newin->addColoured(newqin,partons_[0]->id()<0);
+ newin->addColoured(newg ,partons_[0]->id()<0);
+ newout->addColoured(newspace,partons_[0]->id()<0);
+ newout->addColoured(newqout,partons_[1]->id()<0);
+ newout->addColoured(newg ,partons_[1]->id()>0);
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newqout ,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ else {
+ ShowerParticlePtr newtime(new_ptr(ShowerParticle(partons_[1],true)));
+ newtime->set5Momentum(ComptonMomenta_[1]+ComptonMomenta_[2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newtime,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),offBranch,
+ HardBranching::Outgoing)));
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),offBranch,
+ HardBranching::Outgoing)));
+ offBranch->addChild(outBranch);
+ offBranch->addChild(g);
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newqin ,newqin->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newtime,newqin->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ }
+ // BGF hardest
+ else {
+ rot_.invert();
+ for(unsigned int ix=0;ix<BGFMomenta_.size();++ix) {
+ BGFMomenta_[ix].transform(rot_);
+ }
+ ShowerParticlePtr newq (new_ptr(ShowerParticle(partons_[1],true)));
+ newq->set5Momentum(BGFMomenta_[1]);
+ ShowerParticlePtr newqbar(new_ptr(ShowerParticle(partons_[0]->CC(),true)));
+ newqbar->set5Momentum(BGFMomenta_[2]);
+ ShowerParticlePtr newg (new_ptr(ShowerParticle(gluon_,false)));
+ newg->set5Momentum(BGFMomenta_[0]);
+ ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[0],false)));
+ newspace->set5Momentum(BGFMomenta_[0]-BGFMomenta_[2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newg,SudakovPtr(),HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),spaceBranch,
+ HardBranching::Incoming)));
+ HardBranchingPtr qbar(new_ptr(HardBranching(newqbar,SudakovPtr(),spaceBranch,
+ HardBranching::Outgoing)));
+ spaceBranch->addChild(offBranch);
+ spaceBranch->addChild(qbar);
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newq,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ allBranchings.push_back(outBranch);
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newq ,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ HardTreePtr newTree(new_ptr(HardTree(allBranchings,spaceBranchings,
+ ShowerInteraction::QCD)));
+ // Set the maximum pt for all other emissions and connect hard and shower tree
+ Energy pT = isCompton ? pTCompton_ : pTBGF_;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ // set maximum pT
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pT);
+ for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+ cjt!=newTree->branchings().end();++cjt) {
+ if(!(*cjt)->branchingParticle()->isFinalState()&&
+ (*cjt)->branchingParticle()->id()==cit->first->progenitor()->id()) {
+ newTree->connect(cit->first->progenitor(),*cjt);
+ tPPtr beam =cit->first->original();
+ if(!beam->parents().empty()) beam=beam->parents()[0];
+ (*cjt)->beam(beam);
+ HardBranchingPtr parent=(*cjt)->parent();
+ while(parent) {
+ parent->beam(beam);
+ parent=parent->parent();
+ };
+ }
+ }
+ }
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ // set maximum pT
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pT);
+ for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+ cjt!=newTree->branchings().end();++cjt) {
+ if((*cjt)->branchingParticle()->isFinalState()&&
+ (*cjt)->branchingParticle()->id()==cit->first->progenitor()->id()) {
+ newTree->connect(cit->first->progenitor(),*cjt);
+ }
+ }
+ }
+ // set the evolution partners and scales
+ ShowerParticleVector particles;
+ for(set<HardBranchingPtr>::iterator cit=newTree->branchings().begin();
+ cit!=newTree->branchings().end();++cit) {
+ particles.push_back((*cit)->branchingParticle());
+ }
+// evolver()->showerModel()->partnerFinder()->
+// setInitialEvolutionScales(particles,true,ShowerInteraction::QCD,true);
+// // Calculate the shower variables:
+// evolver()->showerModel()->kinematicsReconstructor()->
+// deconstructHardJets(newTree,evolver(),ShowerInteraction::QCD);
+// for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+// cjt!=newTree->branchings().end();++cjt) {
+// if(cjt==newTree->branchings().begin()) {
+// (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
+// ++cjt;
+// (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
+// ++cjt;
+// }
+// // incoming
+// if((**cjt).status()==HardBranching::Incoming) {
+// quark[0]->set5Momentum((**cjt).showerMomentum());
+// }
+// // outgoing
+// else {
+// quark[1]->set5Momentum((**cjt).showerMomentum());
+// }
+// }
+ return newTree;
+}
+
+void DISBase::generateCompton() {
+ // maximum value of the xT
+ double xT = sqrt((1.-xB_)/xB_);
+ double xTMin = 2.*pTmin_/sqrt(q2_);
+ double zp;
+ // prefactor
+ double a = alpha_->overestimateValue()*comptonWeight_/Constants::twopi;
+ // loop to generate kinematics
+ double wgt(0.),xp(0.);
+ vector<double> azicoeff;
+ do {
+ wgt = 0.;
+ // intergration variables dxT/xT^3
+ xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
+ // zp
+ zp = UseRandom::rnd();
+ xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
+ // check allowed
+ if(xp<xB_||xp>1.) continue;
+ // phase-space piece of the weight
+ wgt = 8.*(1.-xp)*zp/comptonWeight_;
+ // PDF piece of the weight
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
+ wgt *= pdf_->xfx(beam_,partons_[0],scale,xB_/xp)/
+ pdf_->xfx(beam_,partons_[0],q2_ ,xB_);
+ // me piece of the weight
+ double x2 = 1.-(1.-zp)/xp;
+ azicoeff = ComptonME(xp,x2,xT,false);
+ wgt *= 4./3.*alpha_->ratio(0.25*q2_*sqr(xT))*(azicoeff[0]+0.5*azicoeff[2]);
+ if(wgt>1.||wgt<0.) {
+ ostringstream wstring;
+ wstring << "DISBase::generateCompton() "
+ << "Weight greater than one or less than zero"
+ << "wgt = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(xT>xTMin&&UseRandom::rnd()>wgt);
+ if(xT<=xTMin) {
+ pTCompton_=-GeV;
+ return;
+ }
+ // generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi));
+ phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
+ << "::generateCompton() to"
+ << " generate phi" << Exception::eventerror;
+ // momenta for the configuration
+ Energy Q(sqrt(q2_));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ pTCompton_ = 0.5*Q*xT;
+ ComptonMomenta_.resize(3);
+ ComptonMomenta_[0] = p0;
+ ComptonMomenta_[1] = p1;
+ ComptonMomenta_[2] = p2;
+ ComptonISFS_ = zp>xp;
+// debuggingMatrixElement(false,p0,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
+// leptons_[0],leptons_[1],
+// partons_[0],partons_[1],
+// q2_,phi,x2,x3,xT,zp,xp,azicoeff,false);
+}
+
+void DISBase::generateBGF() {
+ // maximum value of the xT
+ double xT = (1.-xB_)/xB_;
+ double xTMin = 2.*max(pTmin_,pTCompton_)/sqrt(q2_);
+ double zp;
+ // prefactor
+ double a = alpha_->overestimateValue()*BGFWeight_/Constants::twopi;
+ // loop to generate kinematics
+ double wgt(0.),xp(0.);
+ vector<double> azicoeff;
+ do {
+ wgt = 0.;
+ // intergration variables dxT/xT^3
+ xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
+ // zp
+ zp = UseRandom::rnd();
+ xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
+ // check allowed
+ if(xp<xB_||xp>1.) continue;
+ // phase-space piece of the weight
+ wgt = 8.*sqr(1.-xp)*zp/BGFWeight_;
+ // PDF piece of the weight
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
+ wgt *= pdf_->xfx(beam_,gluon_ ,scale,xB_/xp)/
+ pdf_->xfx(beam_,partons_[0],q2_ ,xB_);
+ // me piece of the weight
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ azicoeff = BGFME(xp,x2,x3,xT,false);
+ wgt *= 0.5*alpha_->ratio(0.25*q2_*sqr(xT))*
+ (azicoeff[0]+0.5*azicoeff[2]);
+ if(wgt>1.||wgt<0.) {
+ ostringstream wstring;
+ wstring << "DISBase::generateBGF() "
+ << "Weight greater than one or less than zero"
+ << "wgt = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(xT>xTMin&&UseRandom::rnd()>wgt);
+ if(xT<=xTMin) {
+ pTBGF_=-GeV;
+ return;
+ }
+ // generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi));
+ phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
+ << "::generateBGF() to"
+ << " generate phi" << Exception::eventerror;
+ // momenta for the configuration
+ Energy Q(sqrt(q2_));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ pTBGF_=0.5*Q*xT;
+ BGFMomenta_.resize(3);
+ BGFMomenta_[0]=p0;
+ BGFMomenta_[1]=p1;
+ BGFMomenta_[2]=p2;
+// debuggingMatrixElement(true,p0,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
+// leptons_[0],leptons_[1],
+// partons_[0],partons_[1],
+// q2_,phi,x2,x3,xT,zp,xp,azicoeff,false);
+}
+
+int DISBase::nDim() const {
+ return HwMEBase::nDim() + (contrib_>0 ? 1 : 0 );
+}
+
+bool DISBase::generateKinematics(const double * r) {
+ // Born kinematics
+ if(!HwMEBase::generateKinematics(r)) return false;
+ if(contrib_!=0) {
+ // hadron and momentum fraction
+ if(HadronMatcher::Check(*lastParticles().first->dataPtr())) {
+ hadron_ = dynamic_ptr_cast<tcBeamPtr>(lastParticles().first->dataPtr());
+ xB_ = lastX1();
+ }
+ else {
+ hadron_ = dynamic_ptr_cast<tcBeamPtr>(lastParticles().second->dataPtr());
+ xB_ = lastX2();
+ }
+ // Q2
+ q2_ = -(meMomenta()[0]-meMomenta()[2]).m2();
+ // xp
+ int ndim=nDim();
+ double rhomin = pow(1.-xB_,1.-power_);
+ double rho = r[ndim-1]*rhomin;
+ xp_ = 1.-pow(rho,1./(1.-power_));
+ jac_ = rhomin/(1.-power_)*pow(1.-xp_,power_);
+ jacobian(jacobian()*jac_);
+ }
+ return true;
+}
+
+Energy2 DISBase::scale() const {
+ return scaleOpt_ == 1 ?
+ -sqr(scaleFact_)*tHat() : sqr(scaleFact_*muF_);
+}
+
+CrossSection DISBase::dSigHatDR() const {
+ return NLOWeight()*HwMEBase::dSigHatDR();
+}
+
+double DISBase::NLOWeight() const {
+ // If only leading order is required return 1:
+ if(contrib_==0) return 1.;
+ // scale and prefactors
+ Energy2 mu2(scale());
+ double aS = SM().alphaS(mu2);
+ double CFfact = 4./3.*aS/Constants::twopi;
+ double TRfact = 1./2.*aS/Constants::twopi;
+ // LO + dipole subtracted virtual + collinear quark bit with LO pdf
+ double virt = 1.+CFfact*(-4.5-1./3.*sqr(Constants::pi)+1.5*log(q2_/mu2/(1.-xB_))
+ +2.*log(1.-xB_)*log(q2_/mu2)+sqr(log(1.-xB_)));
+ virt /= jac_;
+ // PDF from leading-order
+ double loPDF = hadron_->pdf()->xfx(hadron_,mePartonData()[1],mu2,xB_)/xB_;
+ // NLO gluon PDF
+ tcPDPtr gluon = getParticleData(ParticleID::g);
+ double gPDF = hadron_->pdf()->xfx(hadron_,gluon,mu2,xB_/xp_)*xp_/xB_;
+ // NLO quark PDF
+ double qPDF = hadron_->pdf()->xfx(hadron_,mePartonData()[1],mu2,xB_/xp_)*xp_/xB_;
+ // collinear counterterms
+ // gluon
+ double collg =
+ TRfact/xp_*gPDF*(2.*xp_*(1.-xp_)+(sqr(xp_)+sqr(1.-xp_))*log((1.-xp_)*q2_/xp_/mu2));
+ // quark
+ double collq =
+ CFfact/xp_*qPDF*(1-xp_-2./(1.-xp_)*log(xp_)-(1.+xp_)*log((1.-xp_)/xp_*q2_/mu2))+
+ CFfact/xp_*(qPDF-xp_*loPDF)*(2./(1.-xp_)*log(q2_*(1.-xp_)/mu2)-1.5/(1.-xp_));
+ // calculate the A coefficient for the real pieces
+ double a(A(mePartonData()[0],mePartonData()[2],
+ mePartonData()[1],mePartonData()[3],q2_));
+ // cacluate lepton kinematic variables
+ Lorentz5Momentum q = meMomenta()[0]-meMomenta()[2];
+ double yB = (q*meMomenta()[1])/(meMomenta()[0]*meMomenta()[1]);
+ double l = 2./yB-1.;
+ // q -> qg term
+ double realq = CFfact/xp_/(1.+a*l+sqr(l))*qPDF/loPDF*
+ (2.+2.*sqr(l)-xp_+3.*xp_*sqr(l)+a*l*(2.*xp_+1.));
+ // g -> q qbar term
+ double realg =-TRfact/xp_/(1.+a*l+sqr(l))*gPDF/loPDF*
+ ((1.+sqr(l)+2.*(1.-3.*sqr(l))*xp_*(1.-xp_))
+ +2.*a*l*(1.-2.*xp_*(1.-xp_)));
+ // return the full result
+ double wgt = virt+((collq+collg)/loPDF+realq+realg);
+ // double f2g = gPDF/xp_*TRfact*((sqr(1-xp_)+sqr(xp_))*log((1-xp_)/xp_)+
+ // 8*xp_*(1.-xp_)-1.);
+ // double f2q =
+ // loPDF/jac_*(1.+CFfact*(-1.5*log(1.-xB_)+sqr(log(1.-xB_))
+ // -sqr(Constants::pi)/3.-4.5))
+ // +qPDF *CFfact/xp_*(3.+2.*xp_-(1.+xp_)*log(1.-xp_)
+ // -(1.+sqr(xp_))/(1.-xp_)*log(xp_))
+ // +(qPDF-xp_*loPDF)*CFfact/xp_*(2.*log(1.-xp_)/(1.-xp_)-1.5/(1.-xp_));
+ // double wgt = (f2g+f2q)/loPDF;
+ return contrib_ == 1 ? max(0.,wgt) : max(0.,-wgt);
+}
diff --git a/MatrixElement/DIS/DISBase.h b/MatrixElement/DIS/DISBase.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/DIS/DISBase.h
@@ -0,0 +1,467 @@
+// -*- C++ -*-
+#ifndef HERWIG_DISBase_H
+#define HERWIG_DISBase_H
+//
+// This is the declaration of the DISBase class.
+//
+
+#include "Herwig++/MatrixElement/HwMEBase.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The DISBase class is the base class for the implementation
+ * of DIS type processes including corrections in both the old
+ * fashioned matrix element and POWHEG approaches
+ *
+ * @see \ref DISBaseInterfaces "The interfaces"
+ * defined for DISBase.
+ */
+class DISBase: public HwMEBase {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ DISBase();
+
+ /**
+ * The default constructor.
+ */
+ virtual ~DISBase();
+
+ /**
+ * Members for the old-fashioned matrix element correction
+ */
+ //@{
+ /**
+ * 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);
+ //@}
+
+ /**
+ * Members for the POWHEG stype correction
+ */
+ //@{
+ /**
+ * Has a POWHEG style correction
+ */
+ virtual bool hasPOWHEGCorrection() {return true;}
+
+ /**
+ * Apply the POWHEG style correction
+ */
+ virtual HardTreePtr generateHardest(ShowerTreePtr,vector<ShowerInteraction::Type>);
+ //@}
+
+public:
+
+ /** @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:
+
+ /** @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();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is an abstract class with persistent data.
+ */
+ static AbstractClassDescription<DISBase> initDISBase;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ DISBase & operator=(const DISBase &);
+
+protected:
+
+ /**
+ * The NLO weight
+ */
+ double NLOWeight() const;
+
+ /**
+ * Calculate the coefficient A for the correlations
+ */
+ virtual double A(tcPDPtr lin, tcPDPtr lout, tcPDPtr qin, tcPDPtr qout,
+ Energy2 scale) const =0;
+
+ /**
+ * Members for the matrix element correction
+ */
+ //@{
+ /**
+ * 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 norm Normalise to the large $l$ value of the ME
+ */
+ vector<double> ComptonME(double xp, double x2, double xperp,
+ bool norm);
+
+ /**
+ * 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 norm Normalise to the large $l$ value of the ME
+ */
+ vector<double> BGFME(double xp, double x2, double x3, double xperp,
+ bool norm);
+ //@}
+
+ /**
+ * Members for the POWHEG correction
+ */
+ //@{
+ /**
+ * Generate a Compton process
+ */
+ void generateCompton();
+
+ /**
+ * Generate a BGF process
+ */
+ void generateBGF();
+ //@}
+
+private:
+
+ /**
+ * 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_;
+ //@}
+
+ /**
+ * Parameters for the POWHEG correction
+ */
+ //@{
+ /**
+ * Weight for the compton channel
+ */
+ double comptonWeight_;
+
+ /**
+ * Weight for the BGF channel
+ */
+ double BGFWeight_;
+
+ /**
+ * Minimum value of \f$p_T\f$
+ */
+ Energy pTmin_;
+ //@}
+
+ /**
+ * Parameters for the point being generated
+ */
+ //@{
+ /**
+ * \f$Q^2\f$
+ */
+ Energy2 q2_;
+
+ /**
+ *
+ */
+ double l_;
+
+ /**
+ * Borm momentum fraction
+ */
+ double xB_;
+
+ /**
+ * Beam particle
+ */
+ tcBeamPtr beam_;
+
+ /**
+ * Partons
+ */
+ tcPDPtr partons_[2];
+
+ /**
+ * Leptons
+ */
+ tcPDPtr leptons_[2];
+
+ /**
+ * PDF object
+ */
+ tcPDFPtr pdf_;
+ /**
+ * Rotation to the Breit frame
+ */
+ LorentzRotation rot_;
+
+ /**
+ * Lepton momenta
+ */
+ Lorentz5Momentum pl_[2];
+
+ /**
+ * Quark momenta
+ */
+ Lorentz5Momentum pq_[2];
+
+ /**
+ * q
+ */
+ Lorentz5Momentum q_;
+
+ /**
+ * Compton parameters
+ */
+ Energy pTCompton_;
+ bool ComptonISFS_;
+ vector<Lorentz5Momentum> ComptonMomenta_;
+
+ /**
+ * BGF parameters
+ */
+ Energy pTBGF_;
+ vector<Lorentz5Momentum> BGFMomenta_;
+ //@}
+
+ /**
+ * The coefficient for the correlations
+ */
+ double acoeff_;
+
+ /**
+ * Coupling
+ */
+ ShowerAlphaPtr alpha_;
+
+ /**
+ * Gluon particle data object
+ */
+ PDPtr gluon_;
+
+private:
+
+ /**
+ * The radiative variables
+ */
+ //@{
+ /**
+ * The \f$x_p\f$ or \f$z\f$ real integration variable
+ */
+ double xp_;
+ //@}
+
+ /**
+ * The 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 DISBase. */
+template <>
+struct BaseClassTrait<Herwig::DISBase,1> {
+ /** Typedef of the first base class of DISBase. */
+ typedef Herwig::HwMEBase NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the DISBase class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::DISBase>
+ : public ClassTraitsBase<Herwig::DISBase> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::DISBase"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * MENeutralCurrentDIS is implemented. It may also include several, space-separated,
+ * libraries if the class MENeutralCurrentDIS 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 "HwMEDIS.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_DISBase_H */
diff --git a/MatrixElement/DIS/MEChargedCurrentDIS.cc b/MatrixElement/DIS/MEChargedCurrentDIS.cc
--- a/MatrixElement/DIS/MEChargedCurrentDIS.cc
+++ b/MatrixElement/DIS/MEChargedCurrentDIS.cc
@@ -1,297 +1,301 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEChargedCurrentDIS class.
//
#include "MEChargedCurrentDIS.h"
#include "ThePEG/Utilities/SimplePhaseSpace.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/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "Herwig++/MatrixElement/HardVertex.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
using namespace Herwig;
MEChargedCurrentDIS::MEChargedCurrentDIS()
: _maxflavour(5), _massopt(0) {
vector<unsigned int> mopt(2,1);
mopt[1] = _massopt;
massOption(mopt);
}
void MEChargedCurrentDIS::doinit() {
- HwMEBase::doinit();
+ DISBase::doinit();
_wp = getParticleData(ThePEG::ParticleID::Wplus );
_wm = getParticleData(ThePEG::ParticleID::Wminus);
// cast the SM pointer to the Herwig SM pointer
tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException()
<< "Must be the Herwig++ StandardModel class in "
<< "MEChargedCurrentDIS::doinit" << Exception::abortnow;
// vertices
_theFFWVertex = hwsm->vertexFFW();
}
void MEChargedCurrentDIS::getDiagrams() const {
// possible quarks
typedef std::vector<pair<long,long> > Pairvector;
Pairvector quarkpair;
quarkpair.reserve(6);
// don't even think of putting 'break' in here!
switch(_maxflavour) {
case 6:
quarkpair.push_back(make_pair(ParticleID::s, ParticleID::t));
quarkpair.push_back(make_pair(ParticleID::d, ParticleID::t));
quarkpair.push_back(make_pair(ParticleID::b, ParticleID::t));
case 5:
quarkpair.push_back(make_pair(ParticleID::b, ParticleID::c));
quarkpair.push_back(make_pair(ParticleID::b, ParticleID::u));
case 4:
quarkpair.push_back(make_pair(ParticleID::s, ParticleID::c));
quarkpair.push_back(make_pair(ParticleID::d, ParticleID::c));
case 3:
quarkpair.push_back(make_pair(ParticleID::s, ParticleID::u));
case 2:
quarkpair.push_back(make_pair(ParticleID::d, ParticleID::u));
default:
;
}
// create the diagrams
for(int il1=11;il1<=14;++il1) {
int il2 = il1%2==0 ? il1-1 : il1+1;
for(unsigned int iz=0;iz<2;++iz) {
tcPDPtr lepin = iz==1 ? getParticleData(il1) : getParticleData(-il1);
tcPDPtr lepout = iz==1 ? getParticleData(il2) : getParticleData(-il2);
tcPDPtr inter = lepin->iCharge()-lepout->iCharge()==3 ? _wp : _wm;
for(unsigned int iq=0;iq<quarkpair.size();++iq) {
tcPDPtr first = getParticleData(quarkpair[iq].first );
tcPDPtr second = getParticleData(quarkpair[iq].second);
if(inter==_wp) {
add(new_ptr((Tree2toNDiagram(3), lepin, inter, first ,
1, lepout, 2, second , -1)));
add(new_ptr((Tree2toNDiagram(3), lepin, inter, second->CC(),
1, lepout, 2, first->CC(), -2)));
}
else {
add(new_ptr((Tree2toNDiagram(3), lepin, inter, second ,
1, lepout, 2, first , -1)));
add(new_ptr((Tree2toNDiagram(3), lepin, inter, first->CC(),
1, lepout, 2, second->CC(), -2)));
}
}
}
}
}
-Energy2 MEChargedCurrentDIS::scale() const {
- return -tHat();
-}
-
unsigned int MEChargedCurrentDIS::orderInAlphaS() const {
return 0;
}
unsigned int MEChargedCurrentDIS::orderInAlphaEW() const {
return 2;
}
Selector<const ColourLines *>
MEChargedCurrentDIS::colourGeometries(tcDiagPtr diag) const {
static ColourLines c1("3 5");
static ColourLines c2("-3 -5");
Selector<const ColourLines *> sel;
if ( diag->id() == -1 )
sel.insert(1.0, &c1);
else
sel.insert(1.0, &c2);
return sel;
}
void MEChargedCurrentDIS::persistentOutput(PersistentOStream & os) const {
os << _theFFWVertex << _maxflavour << _wp << _wm << _massopt;
}
void MEChargedCurrentDIS::persistentInput(PersistentIStream & is, int) {
is >> _theFFWVertex >> _maxflavour >> _wp >> _wm >> _massopt;
}
ClassDescription<MEChargedCurrentDIS> MEChargedCurrentDIS::initMEChargedCurrentDIS;
// Definition of the static class description member.
void MEChargedCurrentDIS::Init() {
static ClassDocumentation<MEChargedCurrentDIS> documentation
("The MEChargedCurrentDIS class implements the matrix elements "
"for leading-order charged current deep inelastic scattering");
static Parameter<MEChargedCurrentDIS,unsigned int> interfaceMaxFlavour
( "MaxFlavour",
"The heaviest incoming quark flavour this matrix element is allowed to handle "
"(if applicable).",
&MEChargedCurrentDIS::_maxflavour, 5, 2, 6, false, false, true);
static Switch<MEChargedCurrentDIS,unsigned int> interfaceMassOption
("MassOption",
"Option for the treatment of the mass of the outgoing quarks",
&MEChargedCurrentDIS::_massopt, 0, false, false);
static SwitchOption interfaceMassOptionMassless
(interfaceMassOption,
"Massless",
"Treat the outgoing quarks as massless",
0);
static SwitchOption interfaceMassOptionMassive
(interfaceMassOption,
"Massive",
"Treat the outgoing quarks as massive",
1);
}
Selector<MEBase::DiagramIndex>
MEChargedCurrentDIS::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) sel.insert(1., i);
return sel;
}
double MEChargedCurrentDIS::helicityME(vector<SpinorWaveFunction> & f1,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder, bool calc) const {
// scale
Energy2 mb2(scale());
// matrix element to be stored
ProductionMatrixElement menew(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
// pick a W boson
tcPDPtr ipart = (mePartonData()[0]->iCharge()-mePartonData()[1]->iCharge())==3 ?
_wp : _wm;
// declare the variables we need
VectorWaveFunction inter;
double me(0.);
Complex diag;
// sum over helicities to get the matrix element
unsigned int hel[4];
unsigned int lhel1,lhel2,qhel1,qhel2;
for(lhel1=0;lhel1<2;++lhel1) {
for(lhel2=0;lhel2<2;++lhel2) {
// intermediate W
inter = _theFFWVertex->evaluate(mb2,3,ipart,f1[lhel1],a1[lhel2]);
for(qhel1=0;qhel1<2;++qhel1) {
for(qhel2=0;qhel2<2;++qhel2) {
hel[0] = lhel1;
hel[1] = qhel1;
hel[2] = lhel2;
hel[3] = qhel2;
if(!lorder) swap(hel[0],hel[2]);
if(!qorder) swap(hel[1],hel[3]);
diag = _theFFWVertex->evaluate(mb2,f2[qhel1],a2[qhel2],inter);
me += norm(diag);
menew(hel[0],hel[1],hel[2],hel[3]) = diag;
}
}
}
}
// spin and colour factor
me *= 0.25;
tcPolarizedBeamPDPtr beam[2] =
{dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[0]),
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[1])};
if( beam[0] || beam[1] ) {
RhoDMatrix rho[2] = {beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()),
beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())};
me = menew.average(rho[0],rho[1]);
}
if(calc) _me.reset(menew);
return me;
}
double MEChargedCurrentDIS::me2() const {
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
bool lorder,qorder;
SpinorWaveFunction l1,q1;
SpinorBarWaveFunction l2,q2;
// lepton wave functions
if(mePartonData()[0]->id()>0) {
lorder=true;
l1 = SpinorWaveFunction (meMomenta()[0],mePartonData()[0],incoming);
l2 = SpinorBarWaveFunction(meMomenta()[2],mePartonData()[2],outgoing);
}
else {
lorder=false;
l1 = SpinorWaveFunction (meMomenta()[2],mePartonData()[2],outgoing);
l2 = SpinorBarWaveFunction(meMomenta()[0],mePartonData()[0],incoming);
}
// quark wave functions
if(mePartonData()[1]->id()>0) {
qorder = true;
q1 = SpinorWaveFunction (meMomenta()[1],mePartonData()[1],incoming);
q2 = SpinorBarWaveFunction(meMomenta()[3],mePartonData()[3],outgoing);
}
else {
qorder = false;
q1 = SpinorWaveFunction (meMomenta()[3],mePartonData()[3],outgoing);
q2 = SpinorBarWaveFunction(meMomenta()[1],mePartonData()[1],incoming);
}
// wavefunctions for various helicities
for(unsigned int ix=0;ix<2;++ix) {
l1.reset(ix); f1.push_back(l1);
l2.reset(ix); a1.push_back(l2);
q1.reset(ix); f2.push_back(q1);
q2.reset(ix); a2.push_back(q2);
}
return helicityME(f1,f2,a1,a2,lorder,qorder,false);
}
void MEChargedCurrentDIS::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);
hard.push_back(sub->outgoing()[1]);
unsigned int order[4]={0,1,2,3};
bool lorder(true),qorder(true);
if(hard[0]->id()<0) {
swap(order[0],order[2]);
lorder = false;
}
if(hard[1]->id()<0) {
swap(order[1],order[3]);
qorder = false;
}
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
SpinorWaveFunction (f1,hard[order[0]],incoming,!lorder,true);
SpinorWaveFunction (f2,hard[order[1]],incoming,!qorder,true);
SpinorBarWaveFunction(a1,hard[order[2]],outgoing, lorder,true);
SpinorBarWaveFunction(a2,hard[order[3]],outgoing, qorder,true);
helicityME(f1,f2,a1,a2,lorder,qorder,false);
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix) {
tSpinPtr spin = hard[ix]->spinInfo();
if(ix<2) {
tcPolarizedBeamPDPtr beam =
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(hard[ix]->dataPtr());
if(beam) spin->rhoMatrix() = beam->rhoMatrix();
}
spin->productionVertex(hardvertex);
}
}
+
+double MEChargedCurrentDIS::A(tcPDPtr lin, tcPDPtr,
+ tcPDPtr qin, tcPDPtr, Energy2) const {
+ double output = 2.;
+ if(qin->id()<0) output *= -1.;
+ if(lin->id()<0) output *= -1;
+ return output;
+}
diff --git a/MatrixElement/DIS/MEChargedCurrentDIS.h b/MatrixElement/DIS/MEChargedCurrentDIS.h
--- a/MatrixElement/DIS/MEChargedCurrentDIS.h
+++ b/MatrixElement/DIS/MEChargedCurrentDIS.h
@@ -1,261 +1,263 @@
// -*- C++ -*-
#ifndef HERWIG_MEChargedCurrentDIS_H
#define HERWIG_MEChargedCurrentDIS_H
//
// This is the declaration of the MEChargedCurrentDIS class.
//
-#include "Herwig++/MatrixElement/HwMEBase.h"
+#include "DISBase.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MEChargedCurrentDIS class provides the matrix elements for
* charged current DIS.
*
* By default both the incoming and outgong quarks are assumed to be massless
* although the mass of the outgoing quark can be included if required. This
* option should be used if top production is included.
*
* @see \ref MEChargedCurrentDISInterfaces "The interfaces"
* defined for MEChargedCurrentDIS.
*/
-class MEChargedCurrentDIS: public HwMEBase {
+class MEChargedCurrentDIS: public DISBase {
public:
/**
* The default constructor.
*/
MEChargedCurrentDIS();
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* The matrix element for the kinematical configuration
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
* @return the matrix element scaled with sHat() to give a
* dimensionless number.
*/
virtual double me2() const;
/**
- * Return the scale associated with the last set phase space point.
- */
- virtual Energy2 scale() const;
-
- /**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Get diagram selector. With the information previously supplied with the
* setKinematics method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
* @param dv the diagrams to be weighted.
* @return a Selector relating the given diagrams to their weights.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
* @param diag the diagram chosen.
* @return the possible colour geometries weighted by their
* relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Construct the vertex of spin correlations.
*/
virtual void constructVertex(tSubProPtr);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Matrix element for \f$\ell q\to \gamma/Z \to \ell q\f$.
* @param f1 Fermion on lepton line
* @param a1 Anti-fermion on lepton line
* @param f2 Fermion on quark line
* @param a2 Anti-fermion on quark line
* @param lorder The order of particles on the lepton line
* @param qorder The order of particles on the quark line
* @param me Whether or not to calculate the matrix element for spin correlations
*/
double helicityME(vector<SpinorWaveFunction> & f1 ,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1 ,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder,
bool me) const;
+ /**
+ * Calculate the coefficient A for the correlations in the hard
+ * radiation
+ */
+ virtual double A(tcPDPtr lin, tcPDPtr lout, tcPDPtr qin, tcPDPtr qout,
+ Energy2 scale) const;
+
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MEChargedCurrentDIS> initMEChargedCurrentDIS;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEChargedCurrentDIS & operator=(const MEChargedCurrentDIS &);
private:
/**
* Pointer to the vertex for the helicity calculations
*/
AbstractFFVVertexPtr _theFFWVertex;
/**
* The allowed flavours of the incoming quarks
*/
unsigned int _maxflavour;
/**
* Option for the mass of the outgoing quarks
*/
unsigned int _massopt;
/**
* Matrix element for spin correlations
*/
ProductionMatrixElement _me;
/**
* Pointers to the intermediates resonances
*/
//@{
/**
* Pointer to the \f$W^+\f$
*/
tcPDPtr _wp;
/**
* Pointer to the \f$W^-\f$
*/
tcPDPtr _wm;
//@}
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MEChargedCurrentDIS. */
template <>
struct BaseClassTrait<Herwig::MEChargedCurrentDIS,1> {
/** Typedef of the first base class of MEChargedCurrentDIS. */
- typedef Herwig::HwMEBase NthBase;
+ typedef Herwig::DISBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MEChargedCurrentDIS class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MEChargedCurrentDIS>
: public ClassTraitsBase<Herwig::MEChargedCurrentDIS> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MEChargedCurrentDIS"; }
/**
* The name of a file containing the dynamic library where the class
* MEChargedCurrentDIS is implemented. It may also include several, space-separated,
* libraries if the class MEChargedCurrentDIS 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 "HwMEDIS.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MEChargedCurrentDIS_H */
diff --git a/MatrixElement/DIS/MENeutralCurrentDIS.cc b/MatrixElement/DIS/MENeutralCurrentDIS.cc
--- a/MatrixElement/DIS/MENeutralCurrentDIS.cc
+++ b/MatrixElement/DIS/MENeutralCurrentDIS.cc
@@ -1,323 +1,357 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MENeutralCurrentDIS class.
//
#include "MENeutralCurrentDIS.h"
#include "ThePEG/Utilities/SimplePhaseSpace.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++/MatrixElement/HardVertex.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
using namespace Herwig;
MENeutralCurrentDIS::MENeutralCurrentDIS()
- : _minflavour(1), _maxflavour(5), _gammaZ(0) {
+ : _minflavour(1), _maxflavour(5), _gammaZ(0),
+ _sinW(0.), _cosW(0.), _mz2(ZERO) {
vector<unsigned int> mopt(2,0);
mopt[0] = 1;
massOption(mopt);
}
void MENeutralCurrentDIS::doinit() {
- HwMEBase::doinit();
+ DISBase::doinit();
_z0 = getParticleData(ThePEG::ParticleID::Z0);
_gamma = getParticleData(ThePEG::ParticleID::gamma);
// cast the SM pointer to the Herwig SM pointer
tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException()
<< "Must be the Herwig++ StandardModel class in "
<< "MENeutralCurrentDIS::doinit" << Exception::abortnow;
// vertices
_theFFZVertex = hwsm->vertexFFZ();
_theFFPVertex = hwsm->vertexFFP();
+ // electroweak parameters
+ _sinW = generator()->standardModel()->sin2ThetaW();
+ _cosW = sqrt(1.-_sinW);
+ _sinW = sqrt(_sinW);
+ _mz2 = sqr(_z0->mass());
}
void MENeutralCurrentDIS::getDiagrams() const {
// which intermediates to include
bool gamma = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// create the diagrams
for(int ix=11;ix<=14;++ix) {
for(unsigned int iz=0;iz<2;++iz) {
tPDPtr lep = getParticleData(ix);
if(iz==1) lep = lep->CC();
for(int iy=_minflavour;iy<=_maxflavour;++iy) {
tPDPtr quark = getParticleData(iy);
// lepton quark scattering via gamma and Z
if(gamma) add(new_ptr((Tree2toNDiagram(3), lep, _gamma, quark,
1, lep, 2, quark, -1)));
if(Z0) add(new_ptr((Tree2toNDiagram(3), lep, _z0 , quark,
1, lep, 2, quark, -2)));
// lepton antiquark scattering via gamma and Z
quark = quark->CC();
if(gamma) add(new_ptr((Tree2toNDiagram(3), lep, _gamma, quark,
1, lep, 2, quark, -3)));
if(Z0) add(new_ptr((Tree2toNDiagram(3), lep, _z0 , quark,
1, lep, 2, quark, -4)));
}
}
}
}
-Energy2 MENeutralCurrentDIS::scale() const {
- return -tHat();
-}
-
unsigned int MENeutralCurrentDIS::orderInAlphaS() const {
return 0;
}
unsigned int MENeutralCurrentDIS::orderInAlphaEW() const {
return 2;
}
Selector<const ColourLines *>
MENeutralCurrentDIS::colourGeometries(tcDiagPtr diag) const {
static ColourLines c1("3 5");
static ColourLines c2("-3 -5");
Selector<const ColourLines *> sel;
if ( diag->id() == -1 || diag->id() == -2 )
sel.insert(1.0, &c1);
else
sel.insert(1.0, &c2);
return sel;
}
void MENeutralCurrentDIS::persistentOutput(PersistentOStream & os) const {
os << _minflavour << _maxflavour << _gammaZ << _theFFZVertex << _theFFPVertex
- << _gamma << _z0;
+ << _gamma << _z0 << _sinW << _cosW << ounit(_mz2,GeV2);
}
void MENeutralCurrentDIS::persistentInput(PersistentIStream & is, int) {
is >> _minflavour >> _maxflavour >> _gammaZ >> _theFFZVertex >> _theFFPVertex
- >> _gamma >> _z0;
+ >> _gamma >> _z0 >> _sinW >> _cosW >> iunit(_mz2,GeV2) ;
}
ClassDescription<MENeutralCurrentDIS> MENeutralCurrentDIS::initMENeutralCurrentDIS;
// Definition of the static class description member.
void MENeutralCurrentDIS::Init() {
static ClassDocumentation<MENeutralCurrentDIS> documentation
("The MENeutralCurrentDIS class implements the matrix elements for leading-order "
"neutral current deep inelastic scattering.");
static Parameter<MENeutralCurrentDIS,int> interfaceMaxFlavour
("MaxFlavour",
"The highest incoming quark flavour this matrix element is allowed to handle",
&MENeutralCurrentDIS::_maxflavour, 5, 1, 5,
false, false, Interface::limited);
static Parameter<MENeutralCurrentDIS,int> interfaceMinFlavour
("MinFlavour",
"The lightest incoming quark flavour this matrix element is allowed to handle",
&MENeutralCurrentDIS::_minflavour, 1, 1, 5,
false, false, Interface::limited);
static Switch<MENeutralCurrentDIS,unsigned int> interfaceGammaZ
("GammaZ",
"Which terms to include",
&MENeutralCurrentDIS::_gammaZ, 0, false, false);
static SwitchOption interfaceGammaZAll
(interfaceGammaZ,
"All",
"Include both gamma and Z terms",
0);
static SwitchOption interfaceGammaZGamma
(interfaceGammaZ,
"Gamma",
"Only include the photon",
1);
static SwitchOption interfaceGammaZZ
(interfaceGammaZ,
"Z",
"Only include the Z",
2);
}
Selector<MEBase::DiagramIndex>
MENeutralCurrentDIS::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if ( diags[i]->id() == -1 || diags[i]->id() == -3 ) sel.insert(meInfo()[0], i);
else if ( diags[i]->id() == -2 || diags[i]->id() == -4 ) sel.insert(meInfo()[1], i);
}
return sel;
}
double MENeutralCurrentDIS::helicityME(vector<SpinorWaveFunction> & f1,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder,
bool calc) const {
// scale
Energy2 mb2(scale());
// matrix element to be stored
ProductionMatrixElement menew (PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
ProductionMatrixElement gamma (PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
ProductionMatrixElement Zboson(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
// which intermediates to include
bool gam = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// declare the variables we need
VectorWaveFunction inter[2];
double me[3]={0.,0.,0.};
Complex diag1,diag2;
// sum over helicities to get the matrix element
unsigned int hel[4];
unsigned int lhel1,lhel2,qhel1,qhel2;
for(lhel1=0;lhel1<2;++lhel1) {
for(lhel2=0;lhel2<2;++lhel2) {
// intermediate for photon
if(gam) inter[0]=_theFFPVertex->evaluate(mb2,1,_gamma,f1[lhel1],a1[lhel2]);
// intermediate for Z
if(Z0) inter[1]=_theFFZVertex->evaluate(mb2,1,_z0 ,f1[lhel1],a1[lhel2]);
for(qhel1=0;qhel1<2;++qhel1) {
for(qhel2=0;qhel2<2;++qhel2) {
hel[0] = lhel1;
hel[1] = qhel1;
hel[2] = lhel2;
hel[3] = qhel2;
if(!lorder) swap(hel[0],hel[2]);
if(!qorder) swap(hel[1],hel[3]);
// first the photon exchange diagram
diag1 = gam ?
_theFFPVertex->evaluate(mb2,f2[qhel1],a2[qhel2],inter[0]) : 0.;
// then the Z exchange diagram
diag2 = Z0 ?
_theFFZVertex->evaluate(mb2,f2[qhel1],a2[qhel2],inter[1]) : 0.;
// add up squares of individual terms
me[1] += norm(diag1);
gamma (hel[0],hel[1],hel[2],hel[3]) = diag1;
me[2] += norm(diag2);
Zboson(hel[0],hel[1],hel[2],hel[3]) = diag2;
// the full thing including interference
diag1 += diag2;
me[0] += norm(diag1);
menew(hel[0],hel[1],hel[2],hel[3]) = diag1;
}
}
}
}
// spin and colour factor
double colspin = 0.25;
// results
for(int ix=0;ix<3;++ix) me[ix] *= colspin;
tcPolarizedBeamPDPtr beam[2] =
{dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[0]),
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[1])};
if( beam[0] || beam[1] ) {
RhoDMatrix rho[2] = {beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()),
beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())};
me[0] = menew .average(rho[0],rho[1]);
me[1] = gamma .average(rho[0],rho[1]);
me[2] = Zboson.average(rho[0],rho[1]);
}
DVector save;
save.push_back(me[1]);
save.push_back(me[2]);
meInfo(save);
if(calc) _me.reset(menew);
// analytic expression for testing
// double test = 8.*sqr(4.*Constants::pi*generator()->standardModel()->alphaEM(mb2))*
// sqr(double(mePartonData()[1]->iCharge())/3.)/sqr(tHat())
// *(sqr(sHat())+sqr(uHat())+4.*sqr(mePartonData()[0]->mass())*tHat())/4.;
// cerr << "testing me " << me[0]/test << "\n";
return me[0];
}
double MENeutralCurrentDIS::me2() const {
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
bool lorder,qorder;
SpinorWaveFunction l1,q1;
SpinorBarWaveFunction l2,q2;
// lepton wave functions
if(mePartonData()[0]->id()>0) {
lorder=true;
l1 = SpinorWaveFunction (meMomenta()[0],mePartonData()[0],incoming);
l2 = SpinorBarWaveFunction(meMomenta()[2],mePartonData()[2],outgoing);
}
else {
lorder=false;
l1 = SpinorWaveFunction (meMomenta()[2],mePartonData()[2],outgoing);
l2 = SpinorBarWaveFunction(meMomenta()[0],mePartonData()[0],incoming);
}
// quark wave functions
if(mePartonData()[1]->id()>0) {
qorder = true;
q1 = SpinorWaveFunction (meMomenta()[1],mePartonData()[1],incoming);
q2 = SpinorBarWaveFunction(meMomenta()[3],mePartonData()[3],outgoing);
}
else {
qorder = false;
q1 = SpinorWaveFunction (meMomenta()[3],mePartonData()[3],outgoing);
q2 = SpinorBarWaveFunction(meMomenta()[1],mePartonData()[1],incoming);
}
// wavefunctions for various helicities
for(unsigned int ix=0;ix<2;++ix) {
l1.reset(ix); f1.push_back(l1);
l2.reset(ix); a1.push_back(l2);
q1.reset(ix); f2.push_back(q1);
q2.reset(ix); a2.push_back(q2);
}
return helicityME(f1,f2,a1,a2,lorder,qorder,false);
}
void MENeutralCurrentDIS::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);
hard.push_back(sub->outgoing()[1]);
unsigned int order[4]={0,1,2,3};
bool lorder(true),qorder(true);
if(hard[0]->id()<0) {
swap(order[0],order[2]);
lorder = false;
}
if(hard[1]->id()<0) {
swap(order[1],order[3]);
qorder = false;
}
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
SpinorWaveFunction (f1,hard[order[0]],incoming,!lorder,true);
SpinorWaveFunction (f2,hard[order[1]],incoming,!qorder,true);
SpinorBarWaveFunction(a1,hard[order[2]],outgoing, lorder,true);
SpinorBarWaveFunction(a2,hard[order[3]],outgoing, qorder,true);
helicityME(f1,f2,a1,a2,lorder,qorder,false);
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix) {
tSpinPtr spin = hard[ix]->spinInfo();
if(ix<2) {
tcPolarizedBeamPDPtr beam =
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(hard[ix]->dataPtr());
if(beam) spin->rhoMatrix() = beam->rhoMatrix();
}
spin->productionVertex(hardvertex);
}
}
-
-
+double MENeutralCurrentDIS::A(tcPDPtr lin, tcPDPtr,
+ tcPDPtr qin, tcPDPtr, Energy2 q2) const {
+ // photon only
+ if(_gammaZ==1) return 0.;
+ // everything else
+ double r = _gammaZ==0 || _gammaZ==2 ? double(q2/(q2+_mz2)) : 0.;
+ double eL,eQ,cvL,caL,cvQ,caQ;;
+ if(abs(lin->id())%2==0) {
+ eL = _gammaZ==0 ? generator()->standardModel()->enu() : 0.;
+ cvL = 0.25*generator()->standardModel()->vnu();
+ caL = 0.25*generator()->standardModel()->anu();
+ }
+ else {
+ eL = _gammaZ==0 ? generator()->standardModel()->ee() : 0.;
+ cvL = 0.25*generator()->standardModel()->ve();
+ caL = 0.25*generator()->standardModel()->ae();
+ }
+ if(abs(qin->id())%2==0) {
+ eQ = _gammaZ==0 ? generator()->standardModel()->eu() : 0.;
+ cvQ = 0.25*generator()->standardModel()->vu();
+ caQ = 0.25*generator()->standardModel()->au();
+ }
+ else {
+ eQ = _gammaZ==0 ? generator()->standardModel()->ed() : 0.;
+ cvQ = 0.25*generator()->standardModel()->vd();
+ caQ = 0.25*generator()->standardModel()->ad();
+ }
+ double output = 4.*r*caL*caQ*(eL*eQ+2.*r*cvL*cvQ/sqr(_sinW*_cosW))
+ /sqr(_sinW*_cosW)/(sqr(eL*eQ)+2.*eL*eQ*r/sqr(_cosW*_sinW)*cvL*cvQ
+ +sqr(r/sqr(_cosW*_sinW))*(sqr(cvL)+sqr(caL))*(sqr(cvQ)+sqr(caQ)));
+ if(qin->id()<0) output *= -1.;
+ if(lin->id()<0) output *= -1.;
+ return output;
+}
diff --git a/MatrixElement/DIS/MENeutralCurrentDIS.h b/MatrixElement/DIS/MENeutralCurrentDIS.h
--- a/MatrixElement/DIS/MENeutralCurrentDIS.h
+++ b/MatrixElement/DIS/MENeutralCurrentDIS.h
@@ -1,280 +1,307 @@
// -*- C++ -*-
#ifndef HERWIG_MENeutralCurrentDIS_H
#define HERWIG_MENeutralCurrentDIS_H
//
// This is the declaration of the MENeutralCurrentDIS class.
//
-#include "Herwig++/MatrixElement/HwMEBase.h"
+#include "DISBase.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MENeutralCurrentDIS class provides the matrix elements for
* neutral current DIS.
*
* For consistency both the incoming and outgoing quarks are assumed to be massless.
*
* @see \ref MENeutralCurrentDISInterfaces "The interfaces"
* defined for MENeutralCurrentDIS.
*/
-class MENeutralCurrentDIS: public HwMEBase {
+class MENeutralCurrentDIS: public DISBase {
public:
/**
* The default constructor.
*/
MENeutralCurrentDIS();
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* The matrix element for the kinematical configuration
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
* @return the matrix element scaled with sHat() to give a
* dimensionless number.
*/
virtual double me2() const;
/**
- * Return the scale associated with the last set phase space point.
- */
- virtual Energy2 scale() const;
-
- /**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Get diagram selector. With the information previously supplied with the
* setKinematics method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
* @param dv the diagrams to be weighted.
* @return a Selector relating the given diagrams to their weights.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
* @param diag the diagram chosen.
* @return the possible colour geometries weighted by their
* relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Construct the vertex of spin correlations.
*/
virtual void constructVertex(tSubProPtr);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Matrix element for \f$\ell q\to \gamma/Z \to \ell q\f$.
* @param f1 Fermion on lepton line
* @param a1 Anti-fermion on lepton line
* @param f2 Fermion on quark line
* @param a2 Anti-fermion on quark line
* @param lorder The order of particles on the lepton line
* @param qorder The order of particles on the quark line
* @param me Whether or not to calculate the matrix element for spin correlations
*/
double helicityME(vector<SpinorWaveFunction> & f1 ,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1 ,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder,
bool me) const;
+ /**
+ * Option for treatment of $\gamma/Z\f$ terms
+ */
+ inline unsigned int gammaZOption() const {return _gammaZ;}
+
+ /**
+ * Calculate the coefficient A for the correlations in the hard
+ * radiation
+ */
+ virtual double A(tcPDPtr lin, tcPDPtr lout, tcPDPtr qin, tcPDPtr qout,
+ Energy2 scale) 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 {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<MENeutralCurrentDIS> initMENeutralCurrentDIS;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MENeutralCurrentDIS & operator=(const MENeutralCurrentDIS &);
private:
/**
* Pointer to the vertices for the helicity calculations
*/
//@{
/**
* Pointer to the Z vertex
*/
AbstractFFVVertexPtr _theFFZVertex;
/**
* Pointer to the photon vertex
*/
AbstractFFVVertexPtr _theFFPVertex;
//@}
/**
* Pointers to the intermediate resonances
*/
//@{
/**
* Pointer to the Z ParticleData object
*/
tcPDPtr _z0;
/**
* Pointer to the photon ParticleData object
*/
tcPDPtr _gamma;
//@}
/**
* Switches to control the particles in the hard process
*/
//@{
/**
* Minimumflavour of the incoming quarks
*/
int _minflavour;
/**
* Maximum flavour of the incoming quarks
*/
int _maxflavour;
/**
* Whether to include both \f$Z^0\f$ and \f$\gamma\f$ or only one
*/
unsigned int _gammaZ;
//@}
/**
* Matrix element for spin correlations
*/
ProductionMatrixElement _me;
+ /**
+ * Electroweak parameters
+ */
+ //@{
+ /**
+ * \f$\sin\theta_W\f$
+ */
+ double _sinW;
+
+ /**
+ * \f$\cos\theta_W\f$
+ */
+ double _cosW;
+
+ /**
+ * The square of the Z mass
+ */
+ Energy2 _mz2;
+ //@}
+
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MENeutralCurrentDIS. */
template <>
struct BaseClassTrait<Herwig::MENeutralCurrentDIS,1> {
/** Typedef of the first base class of MENeutralCurrentDIS. */
- typedef Herwig::HwMEBase NthBase;
+ typedef Herwig::DISBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MENeutralCurrentDIS class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MENeutralCurrentDIS>
: public ClassTraitsBase<Herwig::MENeutralCurrentDIS> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MENeutralCurrentDIS"; }
/**
* The name of a file containing the dynamic library where the class
* MENeutralCurrentDIS is implemented. It may also include several, space-separated,
* libraries if the class MENeutralCurrentDIS 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 "HwMEDIS.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MENeutralCurrentDIS_H */
diff --git a/MatrixElement/DIS/Makefile.am b/MatrixElement/DIS/Makefile.am
--- a/MatrixElement/DIS/Makefile.am
+++ b/MatrixElement/DIS/Makefile.am
@@ -1,5 +1,6 @@
pkglib_LTLIBRARIES = HwMEDIS.la
HwMEDIS_la_SOURCES = \
+DISBase.h DISBase.cc \
MENeutralCurrentDIS.cc MENeutralCurrentDIS.h \
MEChargedCurrentDIS.cc MEChargedCurrentDIS.h
HwMEDIS_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 3:0:0
diff --git a/Tests/Makefile.am b/Tests/Makefile.am
--- a/Tests/Makefile.am
+++ b/Tests/Makefile.am
@@ -1,239 +1,251 @@
AUTOMAKE_OPTIONS = -Wno-portability
AM_LDFLAGS += -module -avoid-version -rpath /dummy/path/not/used
EXTRA_DIST = Inputs python Rivet
dist-hook:
rm -rf $(distdir)/Inputs/.svn
rm -rf $(distdir)/python/.svn
rm -rf $(distdir)/Rivet/.svn
EXTRA_LTLIBRARIES = LeptonTest.la GammaTest.la HadronTest.la DISTest.la
if WANT_LIBFASTJET
EXTRA_LTLIBRARIES += HadronJetTest.la LeptonJetTest.la
HadronJetTest_la_SOURCES = \
Hadron/VHTest.h Hadron/VHTest.cc\
Hadron/VTest.h Hadron/VTest.cc\
Hadron/HTest.h Hadron/HTest.cc
HadronJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
HadronJetTest_la_LIBADD = $(FASTJETLIBS)
LeptonJetTest_la_SOURCES = \
Lepton/TopDecay.h Lepton/TopDecay.cc
LeptonJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
LeptonJetTest_la_LIBADD = $(FASTJETLIBS)
endif
LeptonTest_la_SOURCES = \
Lepton/VVTest.h Lepton/VVTest.cc \
Lepton/VBFTest.h Lepton/VBFTest.cc \
Lepton/VHTest.h Lepton/VHTest.cc \
Lepton/FermionTest.h Lepton/FermionTest.cc
GammaTest_la_SOURCES = \
Gamma/GammaMETest.h Gamma/GammaMETest.cc \
Gamma/GammaPMETest.h Gamma/GammaPMETest.cc
DISTest_la_SOURCES = \
DIS/DISTest.h DIS/DISTest.cc
HadronTest_la_SOURCES = \
Hadron/HadronVVTest.h Hadron/HadronVVTest.cc\
Hadron/HadronVBFTest.h Hadron/HadronVBFTest.cc\
Hadron/WHTest.h Hadron/WHTest.cc\
Hadron/ZHTest.h Hadron/ZHTest.cc\
Hadron/VGammaTest.h Hadron/VGammaTest.cc\
Hadron/ZJetTest.h Hadron/ZJetTest.cc\
Hadron/WJetTest.h Hadron/WJetTest.cc\
Hadron/QQHTest.h Hadron/QQHTest.cc
REPO = $(top_builddir)/src/HerwigDefaults.rpo
HERWIG = $(top_builddir)/src/Herwig++
HWREAD = $(HERWIG) read -r $(REPO) -L $(builddir)/.libs
HWRUN = $(HERWIG) run
tests : tests-LEP tests-DIS tests-LHC tests-Gamma
if WANT_LIBFASTJET
tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons \
test-LEP-default test-LEP-Powheg test-LEP-TopDecay
else
tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons
endif
tests-DIS : test-DIS-Charged test-DIS-Neutral
if WANT_LIBFASTJET
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top test-LHC-Bottom \
test-LHC-WHJet test-LHC-ZHJet test-LHC-HJet test-LHC-ZShower test-LHC-WShower\
test-LHC-WHJet-Powheg test-LHC-ZHJet-Powheg test-LHC-HJet-Powheg \
test-LHC-ZShower-Powheg test-LHC-WShower-Powheg
else
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top
endif
tests-Gamma : test-Gamma-FF test-Gamma-WW test-Gamma-P
if WANT_LIBFASTJET
test-LEP-% : Inputs/LEP-%.in LeptonTest.la LeptonJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LEP-% : Inputs/LEP-%.in LeptonTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
Rivet-LEP-% : Rivet/LEP-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-TVT-% : Rivet/TVT-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-DIS-% : Rivet/DIS-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-LHC-% : Rivet/LHC-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-Star-% : Rivet/Star-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-LEP: Rivet-LEP-22 Rivet-LEP-35 Rivet-LEP-44 Rivet-LEP-91 \
Rivet-LEP-133 Rivet-LEP-161 Rivet-LEP-172 Rivet-LEP-177 \
Rivet-LEP-183 Rivet-LEP-189 Rivet-LEP-196 Rivet-LEP-197 \
Rivet-LEP-200 Rivet-LEP-206 Rivet-LEP-14 \
Rivet-LEP-Powheg-14 Rivet-LEP-Powheg-22 \
Rivet-LEP-Powheg-35 Rivet-LEP-Powheg-44 \
Rivet-LEP-Powheg-91 Rivet-LEP-Powheg-133 \
Rivet-LEP-Powheg-161 Rivet-LEP-Powheg-172 \
Rivet-LEP-Powheg-177 Rivet-LEP-Powheg-183 \
Rivet-LEP-Powheg-189 Rivet-LEP-Powheg-196 \
Rivet-LEP-Powheg-197 Rivet-LEP-Powheg-200 \
Rivet-LEP-Powheg-206
for i in LEP-*.aida; do rivet-rmgaps $$i; done;
rivet-rmgaps LEP-91.aida;
rivet-rmgaps LEP-Powheg-91.aida
rm -rf Rivet-LEP
python/merge-LEP LEP
python/merge-LEP LEP-Powheg
rivet-mkhtml -o Rivet-LEP LEP.aida:Hw++ LEP-Powheg.aida:Hw++-Powheg
Rivet-DIS: Rivet-DIS-e--LowQ2 \
- Rivet-DIS-e+-LowQ2 Rivet-DIS-e+-HighQ2
+ Rivet-DIS-e+-LowQ2 Rivet-DIS-e+-HighQ2\
+ Rivet-DIS-Powheg-e--LowQ2 \
+ Rivet-DIS-Powheg-e+-LowQ2 Rivet-DIS-Powheg-e+-HighQ2\
+ Rivet-DIS-NoME-e--LowQ2 \
+ Rivet-DIS-NoME-e+-LowQ2 Rivet-DIS-NoME-e+-HighQ2
rivet-rmgaps DIS-e+-LowQ2.aida
rivet-rmgaps DIS-e--LowQ2.aida
- rivet-rmgaps DIS-e--HighQ2.aida
- python/merge-DIS DIS
- rivet-mkhtml -o Rivet-DIS DIS.aida
+ rivet-rmgaps DIS-e+-HighQ2.aida
+ rivet-rmgaps DIS-Powheg-e+-LowQ2.aida
+ rivet-rmgaps DIS-Powheg-e--LowQ2.aida
+ rivet-rmgaps DIS-Powheg-e+-HighQ2.aida
+ rivet-rmgaps DIS-NoME-e+-LowQ2.aida
+ rivet-rmgaps DIS-NoME-e--LowQ2.aida
+ rivet-rmgaps DIS-NoME-e+-HighQ2.aida
+ python/merge-DIS DIS
+ python/merge-DIS DIS-Powheg
+ python/merge-DIS DIS-NoME
+ rivet-mkhtml -o Rivet-DIS DIS.aida:Hw++ DIS-Powheg.aida:Hw++-Powheg DIS-NoME.aida:Hw++-NoME
Rivet-TVT-WZ: Rivet-TVT-Run-I-Z Rivet-TVT-Powheg-Run-I-Z \
Rivet-TVT-Run-I-W Rivet-TVT-Powheg-Run-I-W \
Rivet-TVT-Run-I-WZ Rivet-TVT-Powheg-Run-I-WZ\
Rivet-TVT-Run-II-Z-e Rivet-TVT-Powheg-Run-II-Z-e \
Rivet-TVT-Run-II-Z-mu Rivet-TVT-Powheg-Run-II-Z-mu \
Rivet-TVT-Run-II-W Rivet-TVT-Powheg-Run-II-W
rivet-rmgaps TVT-Run-II-Z-e.aida;
rivet-rmgaps TVT-Powheg-Run-II-Z-e.aida;
rm -rf Rivet-TVT-WZ
python/merge-TVT-EW TVT-Run-II-W.aida TVT-Run-II-Z-{e,mu}.aida\
TVT-Run-I-{W,Z,WZ}.aida -o TVT-WZ.aida
python/merge-TVT-EW TVT-Powheg-Run-II-W.aida TVT-Powheg-Run-II-Z-{e,mu}.aida\
TVT-Powheg-Run-I-{W,Z,WZ}.aida -o TVT-Powheg-WZ.aida
rivet-mkhtml -o Rivet-TVT-WZ TVT-WZ.aida:Hw++ TVT-Powheg-WZ.aida:Hw++-Powheg
Rivet-TVT-Photon: Rivet-TVT-Run-II-DiPhoton Rivet-TVT-Run-II-PromptPhoton
# Rivet-TVT-Run-I-PromptPhoton
rm -rf Rivet-TVT-Photon
python/merge-aida TVT-Run-II-DiPhoton.aida TVT-Run-II-PromptPhoton.aida\
-o TVT-Photon.aida
rivet-mkhtml -o Rivet-TVT-Photon TVT-Photon.aida:Hw++
Rivet-TVT-Jets: Rivet-TVT-Run-II-Jets-1 Rivet-TVT-Run-II-Jets-2 \
Rivet-TVT-Run-II-Jets-3 Rivet-TVT-Run-II-Jets-4 \
Rivet-TVT-Run-II-Jets-5 Rivet-TVT-Run-II-Jets-6 \
Rivet-TVT-Run-II-Jets-7 Rivet-TVT-Run-II-Jets-8 \
Rivet-TVT-Run-II-Jets-9 Rivet-TVT-Run-II-Jets-10\
Rivet-TVT-Run-II-Jets-11 Rivet-TVT-Run-II-UE \
Rivet-TVT-Run-I-Jets-1 Rivet-TVT-Run-I-Jets-2 \
Rivet-TVT-Run-I-Jets-3 Rivet-TVT-Run-I-Jets-4 \
Rivet-TVT-Run-I-Jets-5 Rivet-TVT-Run-I-Jets-6 \
Rivet-TVT-Run-I-Jets-7 Rivet-TVT-Run-I-Jets-8\
Rivet-TVT-Run-I-UE\
Rivet-TVT-630-UE Rivet-TVT-630-Jets-1 \
Rivet-TVT-630-Jets-2 Rivet-TVT-630-Jets-3
rivet-rmgaps TVT-Run-I-Jets-4.aida
rm -rf Rivet-TVT-Jets
python/merge-TVT-Jets TVT
rivet-mkhtml -o Rivet-TVT-Jets TVT-Jets.aida:Hw++
Rivet-LHC-Jets: Rivet-LHC-7-Jets-1 Rivet-LHC-7-Jets-2 \
Rivet-LHC-7-Jets-3 Rivet-LHC-7-Jets-4 \
Rivet-LHC-7-Jets-5 Rivet-LHC-7-Jets-6 \
Rivet-LHC-7-Jets-7 Rivet-LHC-7-Jets-8 \
Rivet-LHC-7-Jets-9 Rivet-LHC-7-Jets-10 \
Rivet-LHC-7-Jets-11 Rivet-LHC-7-Jets-12 \
Rivet-LHC-7-Jets-13 Rivet-LHC-7-UE \
Rivet-LHC-2360-UE Rivet-LHC-900-UE
rm -rf Rivet-LHC-Jets
python/merge-LHC-Jets
rivet-mkhtml -o Rivet-LHC-Jets LHC-Jets.aida:Hw++
Rivet-Star: Rivet-Star-UE Rivet-Star-Jets-1 \
Rivet-Star-Jets-2 Rivet-Star-Jets-3 \
Rivet-Star-Jets-4
rm -rf Rivet-Star
rivet-rmgaps Star-UE.aida
python/merge-Star Star
rivet-mkhtml -o Rivet-Star Star.aida
Rivet-LHC-EW: Rivet-LHC-W-e Rivet-LHC-Powheg-W-e \
Rivet-LHC-W-mu Rivet-LHC-Powheg-W-mu \
Rivet-LHC-Z Rivet-LHC-Powheg-Z \
Rivet-LHC-WW Rivet-LHC-Powheg-WW\
Rivet-LHC-ZZ Rivet-LHC-Powheg-ZZ
rm -rf Rivet-LHC-EW;
python/merge-LHC-EW LHC-{W-e,W-mu,Z,WW,ZZ}.aida -o LHC-EW.aida;
python/merge-LHC-EW LHC-Powheg-{W-e,W-mu,Z,WW,ZZ}.aida -o LHC-Powheg-EW.aida;
rivet-mkhtml -o Rivet-LHC-EW LHC-EW.aida:Hw++ LHC-Powheg-EW.aida:Hw++-Powheg;
Rivet-LHC-Photon: Rivet-LHC-7-PromptPhoton Rivet-LHC-7-DiPhoton
rm -rf Rivet-LHC-Photon
python/merge-aida LHC-7-PromptPhoton.aida LHC-7-DiPhoton.aida -o LHC-Photon.aida
rivet-mkhtml -o Rivet-LHC-Photon LHC-Photon.aida:Hw++
test-Gamma-% : Inputs/Gamma-%.in GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
test-DIS-% : Inputs/DIS-%.in DISTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
if WANT_LIBFASTJET
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la HadronJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
clean-local:
rm -f *.out *.log *.tex *.top *.run *.dump *.mult *.Bmult *.aida
diff --git a/Tests/Rivet/DIS-NoME-e+-HighQ2.in b/Tests/Rivet/DIS-NoME-e+-HighQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-NoME-e+-HighQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-NoME.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 40.
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-NoME-e+-HighQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-NoME-e+-LowQ2.in b/Tests/Rivet/DIS-NoME-e+-LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-NoME-e+-LowQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-NoME.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-NoME-e+-LowQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-NoME-e--LowQ2.in b/Tests/Rivet/DIS-NoME-e--LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-NoME-e--LowQ2.in
@@ -0,0 +1,26 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-NoME.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e-
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 26.7*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1994_S2919893
+#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1995_S3167097
+
+
+cd /Herwig/Generators
+saverun DIS-NoME-e--LowQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-Powheg-e+-HighQ2.in b/Tests/Rivet/DIS-Powheg-e+-HighQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-Powheg-e+-HighQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-Powheg.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 40.
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-Powheg-e+-HighQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-Powheg-e+-LowQ2.in b/Tests/Rivet/DIS-Powheg-e+-LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-Powheg-e+-LowQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-Powheg.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-Powheg-e+-LowQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-Powheg-e--LowQ2.in b/Tests/Rivet/DIS-Powheg-e--LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-Powheg-e--LowQ2.in
@@ -0,0 +1,26 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-Powheg.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e-
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 26.7*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1994_S2919893
+#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1995_S3167097
+
+
+cd /Herwig/Generators
+saverun DIS-Powheg-e--LowQ2 DISGenerator
diff --git a/Tests/Rivet/DISBase-NoME.in b/Tests/Rivet/DISBase-NoME.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DISBase-NoME.in
@@ -0,0 +1,24 @@
+##################################################
+# Example generator based on DIS parameters
+# usage: Herwig++ read DIS.in
+##################################################
+
+set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
+set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
+
+cd /Herwig/MatrixElements/
+# Neutral current DIS
+insert SimpleDIS:MatrixElements[0] MEDISNC
+
+cd /Herwig/Generators
+set DISGenerator:NumberOfEvents 10000000
+set DISGenerator:RandomNumberGenerator:Seed 31122001
+set DISGenerator:DebugLevel 0
+set DISGenerator:PrintEvent 10
+set DISGenerator:MaxErrors 1000000
+set /Herwig/Shower/ShowerHandler:MPIHandler NULL
+set /Herwig/Shower/Evolver:MECorrMode 0
+cd /Herwig/Generators
+
+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
+insert DISGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
diff --git a/Tests/Rivet/DISBase-Powheg.in b/Tests/Rivet/DISBase-Powheg.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DISBase-Powheg.in
@@ -0,0 +1,43 @@
+##################################################
+# Example generator based on DIS parameters
+# usage: Herwig++ read DIS.in
+##################################################
+
+set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
+set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
+
+##################################################
+# Need to use an NLO PDF
+##################################################
+set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+
+##################################################
+# Setup the POWHEG shower
+##################################################
+cd /Herwig/Shower
+set Evolver:IntrinsicPtGaussian 1.9*GeV
+set Evolver:HardEmissionMode POWHEG
+
+##################################################
+# and strong coupling
+##################################################
+create Herwig::O2AlphaS O2AlphaS
+set /Herwig/Model:QCD/RunningAlphaS O2AlphaS
+
+cd /Herwig/MatrixElements/
+# Neutral current DIS
+
+insert SimpleDIS:MatrixElements[0] /Herwig/MatrixElements/PowhegMEDISNC
+
+cd /Herwig/Generators
+set DISGenerator:NumberOfEvents 10000000
+set DISGenerator:RandomNumberGenerator:Seed 31122001
+set DISGenerator:DebugLevel 0
+set DISGenerator:PrintEvent 10
+set DISGenerator:MaxErrors 1000000
+set /Herwig/Shower/ShowerHandler:MPIHandler NULL
+cd /Herwig/Generators
+
+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
+insert DISGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
diff --git a/src/defaults/MatrixElements.in b/src/defaults/MatrixElements.in
--- a/src/defaults/MatrixElements.in
+++ b/src/defaults/MatrixElements.in
@@ -1,215 +1,228 @@
##############################################################################
# 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
# 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
##########################################################
# DIS matrix elements
##########################################################
# neutral current
create Herwig::MENeutralCurrentDIS MEDISNC
+set MEDISNC:Coupling /Herwig/Shower/AlphaQCD
+set MEDISNC:Contribution 0
# charged current
create Herwig::MEChargedCurrentDIS MEDISCC
+set MEDISCC:Coupling /Herwig/Shower/AlphaQCD
+set MEDISCC:Contribution 0
+
+# neutral current (POWHEG)
+create Herwig::MENeutralCurrentDIS PowhegMEDISNC
+set PowhegMEDISNC:Coupling /Herwig/Shower/AlphaQCD
+set PowhegMEDISNC:Contribution 1
+# charged current (POWHEG)
+create Herwig::MEChargedCurrentDIS PowhegMEDISCC
+set PowhegMEDISCC:Coupling /Herwig/Shower/AlphaQCD
+set 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, Nov 19, 5:50 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3799086
Default Alt Text
(129 KB)

Event Timeline