diff --git a/.hgtags b/.hgtags
--- a/.hgtags
+++ b/.hgtags
@@ -1,24 +1,25 @@
168ae2110e964d62fbc1331a1c2e095952a67748 release-2-5-2
3abb4fa42e20e332796c2572334c2d77204cd0e0 release-2-4-2
4796ca080aafd5daa3b7349b015cb1df944428a2 release-2-5-0
76da042f056eb153981b4d005d5474ffb90a5e88 release-2-4-1
81a684a558413c69df314365eabf09893ffd43d8 release-2-6-0
bd75cd00d99f4bdbaed992daf98f0a73c0f91e9b release-2-4-0
ff6ecc8d49ce10299303b050394bd5cb5837f1c3 release-2-5-1
d0389f5453b2c210923e1adc7b872b18269de668 release-2-6-1
f8998033021185942533b824607285feb3fbd2dc release-2-6-1a
cead23e428b9aacaf2d709e722624e54f844498b release-2-6-1b
191db4655439045f912cb21bd905e729d59ec7bc release-2-6-2
edb538156e9c3d64bb842934b4cebf0126aeb9ea release-2-6-3
eb4a104591859ecac18746b1ad54d6aa0c2a5d1a release-2-7-0
568971ac5b3c1d044c9259f2280a8304fc5a62e9 trunk-before-QED
6e3edb6cfeb4ee48687eb4eb3d016026fc59d602 trunk-after-QED
633abb80b571aa23088957df60e9b0000bbb8a22 release-2-7-1
1bdde095d2346c15ee548e5406a96f0fc6d6e0f1 beforeHQ
a0f9fb821396092bdbeee532bcb0bd624f58335b before_MB_merge
270c1e6b34aa7f758f1d9868c4d3e1ec4bf4e709 herwig-7-0-0
6e0f198c1c2603ecd1a0b6cfe40105cda4bd58c5 herwig-7-0-1
566c1de845a8070559cda45b1bdb40afa18cb2cc herwig-7-0-2
f5c4aa956880f2def763ebd57de7b5bfa55cb1db herwig-7-0-3
65282dedfc2e4bec184e68678dbf4c553c968f38 herwig-7-0-4
541e7790b65ed423c86780bf66ec30e6b99b5a18 herwig-7-1-0
+dd35a1c12d57c047169e8c5fb18644972d49c6ac herwig-7-1-1
diff --git a/AUTHORS b/AUTHORS
--- a/AUTHORS
+++ b/AUTHORS
@@ -1,56 +1,56 @@
================================================================================
-Herwig 7.1.0
+Herwig 7.1.1
================================================================================
Please contact <herwig@projects.hepforge.org> for any queries.
--------------------------------------------------------------------------------
Herwig is actively developed by:
--------------------------------------------------------------------------------
Johannes Bellm
Stefan Gieseke
David Grellscheid
Patrick Kirchgaeßer
Frashër Loshaj
Graeme Nail
Andreas Papaefstathiou
Simon Plätzer
Radek Podskubka
Michael Rauch
Christian Reuschle
Peter Richardson
Peter Schichtel
Mike Seymour
Andrzej Siödmok
Stephen Webster
--------------------------------------------------------------------------------
Former authors are:
--------------------------------------------------------------------------------
Ken Arnold
Manuel Bähr
Luca d'Errico
Martyn Gigg
Nadine Fischer
Keith Hamilton
Marco A. Harrendorf
Seyi Latunde-Dada
Daniel Rauch
Alberto Ribon
Christian Röhr
Pavel Růžička
Alex Schofield
Thomas Schuh
Alexander Sherstnev
Philip Stephens
Martin Stoll
Louise Suter
Jon Tully
Bryan Webber
Alix Wilcock
David Winn
Benedikt Zimmermann
diff --git a/MatrixElement/Hadron/MEDiffraction.cc b/MatrixElement/Hadron/MEDiffraction.cc
--- a/MatrixElement/Hadron/MEDiffraction.cc
+++ b/MatrixElement/Hadron/MEDiffraction.cc
@@ -1,876 +1,877 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEDiffraction class.
//
#include "MEDiffraction.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "Herwig/Utilities/Kinematics.h"
MEDiffraction::MEDiffraction()
: HwMEBase(),
deltaOnly(false),
isInRunPhase(false),
- theProtonMass(0.93827203*GeV)
+ theProtonMass(-MeV) // to be set in doinit
{}
void MEDiffraction::getDiagrams() const {
//incoming particles
cPDPair incomingHardons = generator()->eventHandler()->incoming();
tcPDPtr pom = getParticleData(990);
//get incoming particles
tcPDPtr prt11 = getParticleData(incomingHardons.first->id());
tcPDPtr prt12 = getParticleData(incomingHardons.second->id());
//get sign of id
int sign1=0, sign2=0;
sign1 = (incomingHardons.first->id() > 0) ? 1 : -1;
sign2 = (incomingHardons.second->id() > 0) ? 1 : -1;
tcPDPtr prt21 = getParticleData(sign1*2214);//Delta+
tcPDPtr prt22 = getParticleData(sign2*2214);//Delta+
//for the left side
tcPDPtr q11 = getParticleData(sign1*2); //u
tcPDPtr q21 = getParticleData(sign1*1); //d
//for the right side
tcPDPtr q12 = getParticleData(sign2*2); //u
tcPDPtr q22 = getParticleData(sign2*1); //d
//for the left side
tcPDPtr dq11 = getParticleData(sign1*2101); //ud_0
tcPDPtr dq111 = getParticleData(sign1*2103); //ud_1
tcPDPtr dq21 = getParticleData(sign1*2203); //uu_1
//for the right side
tcPDPtr dq12 = getParticleData(sign2*2101); //ud_0
tcPDPtr dq112 = getParticleData(sign2*2103); //ud_1
tcPDPtr dq22 = getParticleData(sign2*2203); //uu_1
tcPDPtr gl = getParticleData(21);//gluon
//switch between dissociation decays to different
//number of clusters or dissociation into delta only
//(Maybe can be automated???)
//(Should be generalized to ppbar, for example!!!)
switch(dissociationDecay){
case 0: //one cluster or only delta in the final state
if(deltaOnly) //only delta in the final state
{
switch (diffDirection){
case 0:
add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt21, 2, prt12, -1)));
break;
case 1:
add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt11, 2, prt22, -1)));
break;
case 2:
add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt21, 2, prt22, -1)));
break;
}
}else
{
//switch between direction of dissociated proton for single diffraction or
//double diffraction
switch (diffDirection){
case 0: //left
//u -- ud_0
add(new_ptr((Tree2toNDiagram(4), prt11, q11, pom, prt12, 3, prt12, 1, dq11, 2, q11, -1)));
//d -- uu_1
add(new_ptr((Tree2toNDiagram(4), prt11, q21, pom, prt12, 3, prt12, 1, dq21, 2, q21, -2)));
break;
case 1: //right
//u -- ud_0
add(new_ptr((Tree2toNDiagram(4), prt11, pom, q12, prt12, 1, prt11, 3, dq12, 2, q12, -1)));
//d -- uu_1
add(new_ptr((Tree2toNDiagram(4), prt11, pom, q22, prt12, 1, prt11, 3, dq22, 2, q22, -2)));
break;
case 2: //double
//u -- ud_0 left u -- ud_0 right
add(new_ptr((Tree2toNDiagram(5), prt11, q11, pom, q12, prt12, 1, dq11, 2, q11, 3, q12, 4, dq12, -1)));
//u -- ud_0 left d -- uu_1 right
add(new_ptr((Tree2toNDiagram(5), prt11, q11, pom, q22, prt12, 1, dq11, 2, q11, 3, q22, 4, dq22, -2)));
//d -- uu_1 left u -- ud_0 right
add(new_ptr((Tree2toNDiagram(5), prt11, q21, pom, q12, prt12, 1,dq21, 2, q21, 3, q12, 4, dq12, -3)));
//d -- uu_1 left d -- uu_1 right
add(new_ptr((Tree2toNDiagram(5), prt11, q21, pom, q22, prt12, 1, dq21, 2, q21, 3, q22, 4, dq22, -4)));
break;
}
}
break;
case 1: //two clusters (cases with ud_1 not included)
switch (diffDirection){
case 0: //left
//u -- ud_0
add(new_ptr((Tree2toNDiagram(5), prt11, q11, gl, pom, prt12, 1, dq11, 2, q11, 3, gl, 4, prt12, -1)));
//d -- uu_1
add(new_ptr((Tree2toNDiagram(5), prt11, q21, gl, pom, prt12, 1, dq21, 2, q21, 3, gl, 4, prt12, -2)));
break;
case 1: //right
//u -- ud_0
add(new_ptr((Tree2toNDiagram(5), prt11, pom, gl, q12, prt12, 1, prt11, 2, gl, 3, q12, 4, dq12, -1)));
//d -- ud_1
add(new_ptr((Tree2toNDiagram(5), prt11, pom, gl, q22, prt12, 1, prt11, 2, gl, 3, q22, 4, dq22, -2)));
break;
case 2: //double
//u -- ud_0 left u -- ud_0 right
add(new_ptr((Tree2toNDiagram(7), prt11, q11, gl, pom, gl, q12, prt12, 1, dq11, 2, q11, 3, gl, 4,
gl, 5, q12, 6, dq12, -1)));
//u -- ud_0 left d -- uu_1 right
add(new_ptr((Tree2toNDiagram(7), prt11, q11, gl, pom, gl, q22, prt12, 1, dq11, 2, q11, 3, gl, 4,
gl, 5, q22, 6, dq22, -2)));
//d -- uu_1 left u -- ud_0 right
add(new_ptr((Tree2toNDiagram(7), prt11, q21, gl, pom, gl, q12, prt12, 1, dq21, 2, q21, 3, gl, 4,
gl, 5, q12, 6, dq12, -3)));
//d -- uu_1 left d -- uu_1 right
add(new_ptr((Tree2toNDiagram(7), prt11, q21, gl, pom, gl, q22, prt12, 1, dq21, 2, q21, 3, gl, 4,
gl, 5, q22, 6, dq22, -4)));
break;
}
break;
}
}
Energy2 MEDiffraction::scale() const {
return sqr(10*GeV);
}
int MEDiffraction::nDim() const {
return 0;
}
void MEDiffraction::setKinematics() {
HwMEBase::setKinematics(); // Always call the base class method first
}
bool MEDiffraction::generateKinematics(const double * ) {
// generate the masses of the particles
for (size_t i = 2; i < meMomenta().size(); ++i)
meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass());
/* sample M12, M22 and t, characterizing the diffractive final state */
const pair<pair<Energy2,Energy2>,Energy2> point = diffractiveMassAndMomentumTransfer();
const Energy2 M12 (point.first.first);
const Energy2 M22 (point.first.second);
const Energy2 t(point.second);
/* construct the hadronic momenta in the lab frame */
const double phi = UseRandom::rnd() * Constants::twopi;
const Energy cmEnergy = generator()->maximumCMEnergy();
const Energy2 s = sqr(cmEnergy);
//proton mass
const Energy2 m2 = sqr( theProtonMass );
const Energy E3 = (s - M22 + M12) / (2.*cmEnergy);
const Energy E4 = (s + M22 - M12) / (2.*cmEnergy);
//Momentum of outgoing proton and dissociated proton
const Energy pprime = sqrt(kallen(s, M12, M22)) / (2.*cmEnergy);
//costheta of scattering angle
double costheta = s*(s + 2*t - 2*m2 - M12 - M22)
/ sqrt( kallen(s, M12, M22)*kallen(s, m2, m2) );
assert(abs(costheta)<=1.);
const Energy pzprime = pprime*costheta;
const Energy pperp = pprime*sqrt(1 - sqr(costheta));
/* momenta in the lab frame */
const Lorentz5Momentum p3 = Lorentz5Momentum(pperp*cos(phi), pperp*sin(phi), pzprime, E3);
const Lorentz5Momentum p4 = Lorentz5Momentum(-pperp*cos(phi), -pperp*sin(phi), -pzprime, E4);
/* decay dissociated proton into quark-diquark */
//squares of constituent masses of quark and diquark
const Energy2 mq2(sqr(mq()));
Energy2 Mx2;
switch(diffDirection){
case 0:
Mx2=M12;
break;
case 1:
Mx2=M22;
break;
}
/* Select between left/right single diffraction and double diffraction */
//check if we want only delta for the excited state
//pair of momenta for double decay for a two cluster case
pair<Lorentz5Momentum,Lorentz5Momentum> momPair, momPair1;
//fraction of momenta
double frac = UseRandom::rnd();
switch(dissociationDecay){
case 0:
if(!deltaOnly)
{
pair<Lorentz5Momentum,Lorentz5Momentum> decayMomenta;
pair<Lorentz5Momentum,Lorentz5Momentum> decayMomentaTwo;
//absolute collinear dissociation of the hadron
const double phiprime = phi;
//aligned with outgoing dissociated proton
const double costhetaprime = costheta;
const double sinthetaprime=sqrt(1-sqr(costhetaprime));
//axis along which diquark from associated proton is aligned
Axis dir = Axis(sinthetaprime*cos(phiprime), sinthetaprime*sin(phiprime), costhetaprime);
switch (diffDirection){
case 0://Left single diffraction
meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z())));
////////////////////////////////////////////////////
do{}
while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second));
///////////
meMomenta()[2].setVect(p4.vect());
meMomenta()[2].setT(p4.t());
meMomenta()[3].setVect(decayMomenta.first.vect());
meMomenta()[3].setT(decayMomenta.first.t());
meMomenta()[4].setVect(decayMomenta.second.vect());
meMomenta()[4].setT(decayMomenta.second.t());
meMomenta()[2].rescaleEnergy();
meMomenta()[3].rescaleEnergy();
meMomenta()[4].rescaleEnergy();
break;
case 1://Right single diffraction
meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z())));
////////////////////////////////////////////////////
do{}
while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomenta.first,decayMomenta.second));
meMomenta()[2].setVect(p3.vect());
meMomenta()[2].setT(p3.t());
meMomenta()[3].setVect(decayMomenta.first.vect());
meMomenta()[3].setT(decayMomenta.first.t());
meMomenta()[4].setVect(decayMomenta.second.vect());
meMomenta()[4].setT(decayMomenta.second.t());
meMomenta()[2].rescaleEnergy();
meMomenta()[3].rescaleEnergy();
meMomenta()[4].rescaleEnergy();
break;
case 2://double diffraction
do{}
while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second));
do{}
while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomentaTwo.first,decayMomentaTwo.second));
meMomenta()[2].setVect(decayMomenta.first.vect());
meMomenta()[2].setT(decayMomenta.first.t());
meMomenta()[3].setVect(decayMomenta.second.vect());
meMomenta()[3].setT(decayMomenta.second.t());
meMomenta()[4].setVect(decayMomentaTwo.second.vect());
meMomenta()[4].setT(decayMomentaTwo.second.t());
meMomenta()[5].setVect(decayMomentaTwo.first.vect());
meMomenta()[5].setT(decayMomentaTwo.first.t());
meMomenta()[2].rescaleEnergy();
meMomenta()[3].rescaleEnergy();
meMomenta()[4].rescaleEnergy();
meMomenta()[5].rescaleEnergy();
break;
}
}else
{
- meMomenta()[2+diffDirection].setVect(p3.vect());
- meMomenta()[2+diffDirection].setT(p3.t());
- meMomenta()[3-diffDirection].setVect(p4.vect());
- meMomenta()[3-diffDirection].setT(p4.t());
+ const auto tmp=diffDirection==1?1:0;
+ meMomenta()[2+tmp].setVect(p3.vect());
+ meMomenta()[2+tmp].setT(p3.t());
+ meMomenta()[3-tmp].setVect(p4.vect());
+ meMomenta()[3-tmp].setT(p4.t());
meMomenta()[2].rescaleEnergy();
meMomenta()[3].rescaleEnergy();
}
break;
case 1:
switch(diffDirection){
case 0:
//quark and diquark masses
meMomenta()[2].setMass(mqq());
meMomenta()[3].setMass(mq());
//gluon constituent mass
meMomenta()[4].setMass(getParticleData(21)->constituentMass());
//outgoing proton
meMomenta()[5].setVect(p4.vect());
meMomenta()[5].setT(p4.t());
//two body decay of the outgoing dissociation proton
do{}
while(!Kinematics::twoBodyDecay(p3,mqq()+mq(),getParticleData(21)->constituentMass(),
p3.vect().unit(),momPair.first,momPair.second));
//put gluon back-to-back with quark-diquark
//set momenta of quark and diquark
frac = mqq()/(mqq()+mq());
meMomenta()[2].setVect(frac*momPair.first.vect());
meMomenta()[2].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq())));
meMomenta()[3].setVect((1-frac)*momPair.first.vect());
meMomenta()[3].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq())));
//set momentum of gluon
meMomenta()[4].setVect(momPair.second.vect());
meMomenta()[4].setT(momPair.second.t());
break;
case 1:
//quark and diquark masses
meMomenta()[5].setMass(mqq());
meMomenta()[4].setMass(mq());
//gluon constituent mass
meMomenta()[3].setMass(getParticleData(21)->constituentMass());
//outgoing proton
meMomenta()[2].setVect(p3.vect());
meMomenta()[2].setT(p3.t());
//two body decay of the outgoing dissociation proton
do{}
while(!Kinematics::twoBodyDecay(p4,mqq()+mq(),getParticleData(21)->constituentMass(),
p4.vect().unit(),momPair.first,momPair.second));
//put gluon back-to-back with quark-diquark
//set momenta of quark and diquark
frac = mqq()/(mqq()+mq());
meMomenta()[5].setVect(frac*momPair.first.vect());
meMomenta()[5].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq())));
meMomenta()[4].setVect((1-frac)*momPair.first.vect());
meMomenta()[4].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq())));
//set momentum of gluon
meMomenta()[3].setVect(momPair.second.vect());
meMomenta()[3].setT(momPair.second.t());
break;
case 2:
//first dissociated proton constituents
meMomenta()[2].setMass(mqq());
meMomenta()[3].setMass(mq());
meMomenta()[4].setMass(getParticleData(21)->constituentMass());
//second dissociated proton constituents
meMomenta()[5].setMass(getParticleData(21)->constituentMass());
meMomenta()[6].setMass(mq());
meMomenta()[7].setMass(mqq());
//two body decay of the outgoing dissociation proton
do{}
while(!Kinematics::twoBodyDecay(p3,mqq()+mq(),getParticleData(21)->constituentMass(),
p3.vect().unit(),momPair.first,momPair.second));
do{}
while(!Kinematics::twoBodyDecay(p4,mqq()+mq(),getParticleData(21)->constituentMass(),
p4.vect().unit(),momPair1.first,momPair1.second));
//put gluon back-to-back with quark-diquark
frac = mqq()/(mqq()+mq());
//first dissociated proton
//set momenta of quark and diquark
meMomenta()[2].setVect(frac*momPair.first.vect());
meMomenta()[2].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq())));
meMomenta()[3].setVect((1-frac)*momPair.first.vect());
meMomenta()[3].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq())));
//set momentum of gluon
meMomenta()[4].setVect(momPair.second.vect());
meMomenta()[4].setT(momPair.second.t());
//first dissociated proton
//set momenta of quark and diquark
meMomenta()[7].setVect(frac*momPair1.first.vect());
meMomenta()[7].setT(sqrt(sqr(frac)*momPair1.first.vect().mag2()+sqr(mqq())));
meMomenta()[6].setVect((1-frac)*momPair1.first.vect());
meMomenta()[6].setT(sqrt(sqr(1-frac)*momPair1.first.vect().mag2()+sqr(mq())));
//set momentum of gluon
meMomenta()[5].setVect(momPair1.second.vect());
meMomenta()[5].setT(momPair1.second.t());
break;
}
meMomenta()[2].rescaleEnergy();
meMomenta()[3].rescaleEnergy();
meMomenta()[4].rescaleEnergy();
meMomenta()[5].rescaleEnergy();
if(diffDirection==2){
meMomenta()[6].rescaleEnergy();
meMomenta()[7].rescaleEnergy();
}
break;
}
jacobian(sqr(cmEnergy)/GeV2);
return true;
}
//Generate masses of dissociated protons and momentum transfer from probability f(M2,t)
//(for single diffraction). Sample according to f(M2,t)=f(M2)f(t|M2).
pair<pair<Energy2,Energy2>,Energy2> MEDiffraction::diffractiveMassAndMomentumTransfer() const {
Energy2 theM12(ZERO),theM22(ZERO), thet(ZERO);
int count = 0;
//proton mass squared
const Energy2 m2 = sqr(theProtonMass);
//delta mass squared
const Energy2 md2 = sqr(getParticleData(2214)->mass());
Energy2 M2;
bool condition = true;
do {
//check if we want only delta
if(deltaOnly) {
switch(diffDirection){
case 0:
theM12 = md2;
theM22 = m2;
M2 = md2;
thet = randomt(md2);
break;
case 1:
theM22 = md2;
theM12 = m2;
M2 = md2;
thet = randomt(md2);
break;
case 2:
theM12 = md2;
theM22 = md2;
M2 = md2;
thet = doublediffrandomt(theM12,theM22);
break;
}
}
else {
switch (diffDirection){
case 0:
M2=randomM2();
thet = randomt(M2);
theM12=M2;
theM22=m2;
break;
case 1:
theM12=m2;
M2=randomM2();
thet = randomt(M2);
theM22=M2;
break;
case 2:
theM12=randomM2();
theM22=randomM2();
M2=(theM12>theM22) ? theM12: theM22;
thet = doublediffrandomt(theM12,theM22);
break;
}
}
count++;
const Energy cmEnergy = generator()->maximumCMEnergy();
const Energy2 s = sqr(cmEnergy);
if(generator()->maximumCMEnergy()<sqrt(theM12)+sqrt(theM22)) {
condition = true;
}
else {
InvEnergy2 slope;
if(diffDirection==2){
slope = 2*softPomeronSlope()*log(.1+(sqr(cmEnergy)/softPomeronSlope())/(theM12*theM22));
}else{
slope = protonPomeronSlope()
+ 2*softPomeronSlope()*log(sqr(cmEnergy)/M2);
}
const double expmax = exp(slope*tmaxfun(s,m2,M2));
const double expmin = exp(slope*tminfun(s,m2,M2));
//without (1-M2/s) constraint
condition = (UseRandom::rnd()>(protonPomeronSlope()*GeV2)*(expmax-expmin)/(slope*GeV2))
||((theM12/GeV2)*(theM22/GeV2)>=(sqr(cmEnergy)/GeV2)/(softPomeronSlope()*GeV2));
}
}
while(condition);
return make_pair (make_pair(theM12,theM22),thet);
}
//Decay of the excited proton to quark-diquark
pair<Lorentz5Momentum,Lorentz5Momentum> MEDiffraction::twoBodyDecayMomenta(Lorentz5Momentum pp) const{
//Decay of the excited proton
const Energy2 Mx2(sqr(pp.mass())),mq2(sqr(mq())),mqq2(sqr(mqq()));
const Energy2 psq = ((Mx2-sqr(mq()+mqq()))*(Mx2-sqr(mq()-mqq())))/(4*Mx2);
assert(psq/GeV2>0);
const Energy p(sqrt(psq));
const double phi = UseRandom::rnd() * Constants::twopi;
const double costheta =1-2*UseRandom::rnd();
const double sintheta = sqrt(1-sqr(costheta));
Lorentz5Momentum k1=Lorentz5Momentum(p*sintheta*cos(phi), p*sintheta*sin(phi), p*costheta, sqrt(mq2+psq));
Lorentz5Momentum k2=Lorentz5Momentum(-p*sintheta*cos(phi), -p*sintheta*sin(phi), -p*costheta,sqrt(mqq2+psq));
//find boost to pp center of mass
const Boost betap3 = (pp).findBoostToCM();
//k1 and k2 calculated at p3 center of mass, so boost back
k1.boost(-betap3);
k2.boost(-betap3);
//first is quark, second diquark
return make_pair(k1,k2);
}
Energy2 MEDiffraction::randomt(Energy2 M2) const {
assert(protonPomeronSlope()*GeV2 > 0);
//proton mass
const Energy2 m2 = sqr( theProtonMass );
const Energy cmEnergy = generator()->maximumCMEnergy();
const Energy2 ttmin = tminfun(sqr(cmEnergy),m2,M2);
const Energy2 ttmax = tmaxfun(sqr(cmEnergy),m2,M2);
const InvEnergy2 slope = protonPomeronSlope()
+ 2*softPomeronSlope()*log(sqr(cmEnergy)/M2);
return log( exp(slope*ttmin) +
UseRandom::rnd()*(exp(slope*ttmax) - exp(slope*ttmin)) ) / slope;
}
Energy2 MEDiffraction::doublediffrandomt(Energy2 M12, Energy2 M22) const {
const Energy cmEnergy = generator()->maximumCMEnergy();
const double shift = 0.1;
const InvEnergy2 slope = 2*softPomeronSlope()*log(shift+(sqr(cmEnergy)/softPomeronSlope())/(M12*M22));
const Energy2 ttmin = tminfun(sqr(cmEnergy),M12,M22);
const Energy2 ttmax = tmaxfun(sqr(cmEnergy),M12,M22);
double r = UseRandom::rnd();
Energy2 newVal;
if(slope*ttmax>slope*ttmin) {
newVal = ttmax + log( r + (1.-r)*exp(slope*(ttmin-ttmax)) ) / slope;
}
else {
newVal = ttmin + log( 1. - r + r*exp(slope*(ttmax-ttmin))) / slope;
}
return newVal;
}
Energy2 MEDiffraction::randomM2() const {
const double tmp = 1 - softPomeronIntercept();
const Energy cmEnergy = generator()->maximumCMEnergy();
return sqr(cmEnergy) * pow( pow(M2min()/sqr(cmEnergy),tmp) +
UseRandom::rnd() * (pow(M2max()/sqr(cmEnergy),tmp) - pow(M2min()/sqr(cmEnergy),tmp)),
1.0/tmp );
}
Energy2 MEDiffraction::tminfun(Energy2 s, Energy2 M12, Energy2 M22) const {
const Energy2 m2 = sqr( theProtonMass );
return 0.5/s*(-sqrt(kallen(s, m2, m2)*kallen(s, M12, M22))-sqr(s)+2*s*m2+s*M12+s*M22);
}
Energy2 MEDiffraction::tmaxfun(Energy2 s, Energy2 M12, Energy2 M22) const {
const Energy2 m2 = sqr( theProtonMass );
return 0.5/s*(sqrt(kallen(s, m2, m2)*kallen(s, M12, M22))-sqr(s)+2*s*m2+s*M12+s*M22);
}
double MEDiffraction::me2() const{
return theme2;
}
CrossSection MEDiffraction::dSigHatDR() const {
return me2()*jacobian()/sHat()*sqr(hbarc);
}
unsigned int MEDiffraction::orderInAlphaS() const {
return 0;
}
unsigned int MEDiffraction::orderInAlphaEW() const {
return 0;
}
Selector<MEBase::DiagramIndex>
MEDiffraction::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
if(!deltaOnly){
if(diffDirection<2){
for(unsigned int i = 0; i < diags.size(); i++){
if(diags[0]->id()==-1)
sel.insert(2./3.,i);
else
sel.insert(1./3.,i);
}
}else{
for(unsigned int i = 0; i < diags.size(); i++){
if(diags[0]->id()==-1)
sel.insert(4./9.,i);
else if(diags[0]->id()==-2)
sel.insert(2./9.,i);
else if(diags[0]->id()==-3)
sel.insert(2./9.,i);
else
sel.insert(1./9.,i);
}
}
}else{
sel.insert(1.0,0);
}
return sel;
}
Selector<const ColourLines *>
MEDiffraction::colourGeometries(tcDiagPtr ) const {
Selector<const ColourLines *> sel;
int sign1=0, sign2=0;
sign1 = (generator()->eventHandler()->incoming().first->id() > 0) ? 1 : -1;
sign2 = (generator()->eventHandler()->incoming().second->id() > 0) ? 1 : -1;
switch(dissociationDecay){
case 0:
if(!deltaOnly)
{
if(diffDirection!=2){
if (diffDirection == 0){
if(sign1>0){
static ColourLines dqq0=ColourLines("-6 2 7");
sel.insert(1.0,&dqq0);
}else{
static ColourLines dqq0=ColourLines("6 -2 -7");
sel.insert(1.0,&dqq0);
}
}
else{
if(sign2>0){
static ColourLines dqq1=ColourLines("-6 3 7");
sel.insert(1.0,&dqq1);
}else{
static ColourLines dqq1=ColourLines("6 -3 -7");
sel.insert(1.0,&dqq1);
}
}
}else{
if(sign1>0 && sign2>0){
static ColourLines ddqq0=ColourLines("-6 2 7, -9 4 8");
sel.insert(1.0,&ddqq0);
}else if(sign1<0 && sign2>0){
static ColourLines ddqq0=ColourLines("6 -2 -7, -9 4 8");
sel.insert(1.0,&ddqq0);
}else if(sign1>0&& sign2<0){
static ColourLines ddqq0=ColourLines("-6 2 7, 9 -4 -8");
sel.insert(1.0,&ddqq0);
}else{
static ColourLines ddqq0=ColourLines("6 -2 -7, 9 -4 -8");
sel.insert(1.0,&ddqq0);
}
}
}else
{
static ColourLines cl("");
sel.insert(1.0, &cl);
}
break;
case 1:
switch(diffDirection){
case 0:
static ColourLines clleft("-6 2 3 8, -8 -3 7");
sel.insert(1.0, &clleft);
break;
case 1:
static ColourLines clright("-9 4 3 7, -7 -3 8");
sel.insert(1.0, &clright);
break;
case 2:
static ColourLines cldouble("-8 2 3 10, -10 -3 9, -13 6 5 11, -11 -5 12");
sel.insert(1.0, &cldouble);
break;
}
break;
}
return sel;
}
void MEDiffraction::doinit() {
HwMEBase::doinit();
theProtonMass = getParticleData(2212)->mass();
}
void MEDiffraction::doinitrun() {
HwMEBase::doinitrun();
isInRunPhase = true;
}
IBPtr MEDiffraction::clone() const {
return new_ptr(*this);
}
IBPtr MEDiffraction::fullclone() const {
return new_ptr(*this);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEDiffraction,HwMEBase>
describeHerwigMEDiffraction("Herwig::MEDiffraction", "HwMEHadron.so");
void MEDiffraction::persistentOutput(PersistentOStream & os) const {
os << theme2 << deltaOnly << diffDirection << theprotonPomeronSlope
<< thesoftPomeronIntercept << thesoftPomeronSlope << dissociationDecay
<< ounit(theProtonMass,GeV);
}
void MEDiffraction::persistentInput(PersistentIStream & is, int) {
is >> theme2 >> deltaOnly >> diffDirection >> theprotonPomeronSlope
>> thesoftPomeronIntercept >> thesoftPomeronSlope >> dissociationDecay
>> iunit(theProtonMass,GeV);
}
InvEnergy2 MEDiffraction::protonPomeronSlope() const{
return theprotonPomeronSlope/GeV2;
}
double MEDiffraction::softPomeronIntercept() const {
return thesoftPomeronIntercept;
}
InvEnergy2 MEDiffraction::softPomeronSlope() const {
return thesoftPomeronSlope/GeV2;
}
void MEDiffraction::Init() {
static ClassDocumentation<MEDiffraction> documentation
("There is no documentation for the MEDiffraction class");
static Parameter<MEDiffraction,double> interfaceme2
("DiffractionAmplitude",
"The square of the diffraction amplitude used to determine the "
"cross section.",
&MEDiffraction::theme2, 1.0, 0.00001, 100.0,
false, false, Interface::limited);
static Parameter<MEDiffraction,double> interfaceprotonPomeronSlope
("ProtonPomeronSlope",
"The proton-pomeron slope parameter.",
&MEDiffraction::theprotonPomeronSlope, 10.1, 0.00001, 100.0,
false, false, Interface::limited);
static Parameter<MEDiffraction,double> interfacesoftPomeronIntercept
("SoftPomeronIntercept",
"The soft pomeron intercept.",
&MEDiffraction::thesoftPomeronIntercept, 1.08, 0.00001, 100.0,
false, false, Interface::limited);
static Parameter<MEDiffraction,double> interfacesoftPomeronSlope
("SoftPomeronSlope",
"The soft pomeron slope parameter.",
&MEDiffraction::thesoftPomeronSlope, 0.25, 0.00001, 100.0,
false, false, Interface::limited);
static Switch<MEDiffraction, bool> interfaceDeltaOnly
("DeltaOnly",
"proton-proton to proton-delta only",
&MEDiffraction::deltaOnly, 0, false, false);
static SwitchOption interfaceDeltaOnly0
(interfaceDeltaOnly,"No","Final state with Delta only is OFF", 0);
static SwitchOption interfaceDeltaOnly1
(interfaceDeltaOnly,"Yes","Final state with Delta only is ON", 1);
//Select if the left, right or both protons are excited
static Switch<MEDiffraction, unsigned int> interfaceDiffDirection
("DiffDirection",
"Direction of the excited proton",
&MEDiffraction::diffDirection, 0, false, false);
static SwitchOption left
(interfaceDiffDirection,"Left","Proton moving in the positive z direction", 0);
static SwitchOption right
(interfaceDiffDirection,"Right","Proton moving in the negative z direction", 1);
static SwitchOption both
(interfaceDiffDirection,"Both","Both protons", 2);
//Select if two or three body decay
static Switch<MEDiffraction, unsigned int> interfaceDissociationDecay
("DissociationDecay",
"Number of clusters the dissociated proton decays",
&MEDiffraction::dissociationDecay, 0, false, false);
static SwitchOption one
(interfaceDissociationDecay,"One","Dissociated proton decays into one cluster", 0);
static SwitchOption two
(interfaceDissociationDecay,"Two","Dissociated proton decays into two clusters", 1);
}
diff --git a/MatrixElement/Hadron/Makefile.am b/MatrixElement/Hadron/Makefile.am
--- a/MatrixElement/Hadron/Makefile.am
+++ b/MatrixElement/Hadron/Makefile.am
@@ -1,28 +1,28 @@
pkglib_LTLIBRARIES = HwMEHadron.la
HwMEHadron_la_SOURCES = \
MEqq2gZ2ff.cc MEqq2gZ2ff.h \
MEqq2W2ff.cc MEqq2W2ff.h \
MEPP2GammaJet.h MEPP2GammaJet.cc\
MEQCD2to2.h MEQCD2to2.cc\
MEPP2HiggsJet.h MEPP2HiggsJet.cc\
MEPP2GammaGamma.h MEPP2GammaGamma.cc \
MEPP2QQ.h MEPP2QQ.cc \
MEPP2QQHiggs.h MEPP2QQHiggs.cc \
MEPP2Higgs.h MEPP2Higgs.cc\
MEPP2WH.h MEPP2WH.cc \
MEPP2ZH.h MEPP2ZH.cc \
MEPP2WJet.cc MEPP2WJet.h \
MEPP2ZJet.cc MEPP2ZJet.h \
MEPP2VV.cc MEPP2VV.h \
MEPP2VGamma.cc MEPP2VGamma.h \
MEPP2HiggsVBF.cc MEPP2HiggsVBF.h \
MEPP2SingleTop.cc MEPP2SingleTop.h \
MEMinBias.h MEMinBias.cc \
MEDiffraction.h MEDiffraction.cc
-HwMEHadron_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 8:0:0
+HwMEHadron_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 8:1:0
pkglib_LTLIBRARIES += HwMEHadronFast.la
HwMEHadronFast_la_SOURCES = \
MEQCD2to2Fast.h MEQCD2to2Fast.cc
HwMEHadronFast_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 6:0:0
diff --git a/MatrixElement/Matchbox/Scales/Makefile.am b/MatrixElement/Matchbox/Scales/Makefile.am
--- a/MatrixElement/Matchbox/Scales/Makefile.am
+++ b/MatrixElement/Matchbox/Scales/Makefile.am
@@ -1,32 +1,32 @@
pkglib_LTLIBRARIES = HwMatchboxScales.la
-HwMatchboxScales_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 3:0:0
+HwMatchboxScales_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 4:0:0
HwMatchboxScales_la_SOURCES = \
MatchboxHtScale.h \
MatchboxLeptonMassScale.h \
MatchboxLeptonPtScale.h \
MatchboxParticlePtScale.h \
MatchboxPtScale.h \
MatchboxHtScale.cc \
MatchboxLeptonMassScale.cc \
MatchboxLeptonPtScale.cc \
MatchboxParticlePtScale.cc \
MatchboxPtScale.cc \
MatchboxSHatScale.h \
MatchboxSHatScale.cc \
MatchboxTopMassScale.h \
MatchboxTopMassScale.cc \
MatchboxTopMTScale.h \
MatchboxTopMTScale.cc \
MatchboxTopSumMTScale.h \
MatchboxTopSumMTScale.cc \
MatchboxTopMinMTScale.h \
MatchboxTopMinMTScale.cc \
MatchboxTriVecScales.h \
MatchboxTriVecScales.cc \
MatchboxTopLinearSumMTScale.h \
MatchboxTopLinearSumMTScale.cc \
MatchboxTopIndividualMTScale.h \
MatchboxTopIndividualMTScale.cc
diff --git a/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.cc b/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.cc
--- a/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.cc
+++ b/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.cc
@@ -1,192 +1,192 @@
// -*- C++ -*-
//
// MatchboxTriVecScales.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxTriVecScales class.
//
#include "MatchboxTriVecScales.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
MatchboxTriVecScales::MatchboxTriVecScales()
: theTriVecScaleChoice(1) {}
MatchboxTriVecScales::~MatchboxTriVecScales() {}
IBPtr MatchboxTriVecScales::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxTriVecScales::fullclone() const {
return new_ptr(*this);
}
-Energy2 MatchboxTriVecScales::HtPrimeScale() const {
+Energy2 MatchboxTriVecScales::renormalizationScale() const {
+ // calculate scales
tcPDVector pd (mePartonData().begin() + 2, mePartonData().end());
vector<LorentzMomentum> p (meMomenta().begin() + 2, meMomenta().end());
tcPDPtr t1 = mePartonData()[0];
tcPDPtr t2 = mePartonData()[1];
tcCutsPtr cuts = lastCutsPtr();
theJetFinder->cluster(pd, p, cuts, t1, t2);
Energy sumpartpt = ZERO;
- int foundpart = 0;
- Energy vec1et = ZERO;
- Energy vec2et = ZERO;
- Energy vec3et = ZERO;
+ Energy vecet[3] = {ZERO,ZERO,ZERO};
Energy sumvecet = ZERO;
int foundlept[3] = {0,0,0}; // First entry for e-like, second entry for mu-like, third entry for tau-like
- int ivec1 = 0;
- int ivec2 = 0;
- int ivec3 = 0;
+ int ivec[3] = {0,0,0};
tcPDVector::const_iterator itpd = pd.begin();
int ip = 2;
+
+ double avgy12 = 0;
+ if ( theTriVecScaleChoice==2 ) {
+ LorentzMomentum p1 = LorentzMomentum();
+ LorentzMomentum p2 = LorentzMomentum();
+ int njets = 0;
+ for (vector<LorentzMomentum>::const_iterator itp = p.begin() ;
+ itp != p.end(); ++itp, ++itpd, ++ip )
+ if ( (**itpd).coloured() ) {
+ ++njets;
+ if ( itp->perp() > p1.perp() ) {
+ p1 = *itp;
+ p2 = p1;
+ } else if ( itp->perp() > p2.perp() )
+ p2 = *itp;
+ }
+ if ( njets < 2 )
+ throw Exception() << "MatchboxTriVecScales: Not enough jets in event for HtPrimeModScale!"
+ << Exception::runerror;
+ avgy12 = (p1.rapidity()+p2.rapidity())/2.;
+ }
+
for (vector<LorentzMomentum>::const_iterator itp = p.begin() ;
itp != p.end(); ++itp, ++itpd, ++ip ) {
if ( (**itpd).coloured() ) {
- sumpartpt += (*itp).perp();
- foundpart++;
+ if ( theTriVecScaleChoice==2 )
+ sumpartpt += (*itp).perp()*exp(abs((*itp).rapidity()-avgy12));
+ else
+ sumpartpt += (*itp).perp();
}
if ( !(**itpd).coloured() ) {
if ( abs((**itpd).id())==ParticleID::nu_e || abs((**itpd).id())==ParticleID::eminus ) {
- if ( foundlept[0]==0 ) {
- foundlept[0] +=1;
- ivec1 = ip;
- }
- if ( foundlept[0]==1 ) {
- vec1et = sqrt( (p[ivec1]+p[ip]).m2() + (p[ivec1]+p[ip]).perp2() );
- }
+ if ( ++foundlept[0] == 1 )
+ ivec[0] = ip;
+ if ( foundlept[0] == 2 )
+ vecet[0] = sqrt( (p[ivec[0]]+p[ip]).m2() + (p[ivec[0]]+p[ip]).perp2() );
}
if ( abs((**itpd).id())==ParticleID::nu_mu || abs((**itpd).id())==ParticleID::muminus ) {
- if ( foundlept[1]==0 ) {
- foundlept[1] +=1;
- ivec2 = ip;
- }
- if ( foundlept[1]==1 ) {
- vec2et = sqrt( (p[ivec2]+p[ip]).m2() + (p[ivec2]+p[ip]).perp2() );
- }
+ if ( ++foundlept[1] == 1 )
+ ivec[1] = ip;
+ if ( foundlept[1] == 2 )
+ vecet[1] = sqrt( (p[ivec[1]]+p[ip]).m2() + (p[ivec[1]]+p[ip]).perp2() );
}
if ( abs((**itpd).id())==ParticleID::nu_tau || abs((**itpd).id())==ParticleID::tauminus ) {
- if ( foundlept[2]==0 ) {
- foundlept[2] +=1;
- ivec3 = ip;
- }
- if ( foundlept[2]==1 ) {
- vec3et = sqrt( (p[ivec3]+p[ip]).m2() + (p[ivec3]+p[ip]).perp2() );
- }
+ if ( ++foundlept[2] == 1 )
+ ivec[2] = ip;
+ if ( foundlept[2] == 2 )
+ vecet[2] = sqrt( (p[ivec[2]]+p[ip]).m2() + (p[ivec[2]]+p[ip]).perp2() );
}
}
}
- sumvecet = vec1et+vec2et+vec3et;
// Check for consistency in number of lepton pairs and members therein
for (int i = 0; i<3; ++i ) {
- if (foundlept[i]>2 || foundlept[i]%2!=0)
- throw Exception() << "MatchboxTriVecScales::HtPrimeScale(): Inconsistency in number of lepton pairs and members therein!"
+ if ( foundlept[i] != 2 )
+ throw Exception() << "MatchboxTriVecScales: Inconsistency in number of lepton pairs and members therein!"
<< Exception::runerror;
}
- return sqr(sumpartpt+sumvecet);
+ sumvecet = vecet[0]+vecet[1]+vecet[2];
-}
-
-Energy2 MatchboxTriVecScales::HtPrimeModScale() const {
- throw Exception() << "MatchboxTriVecScales::HtPrimeModScale(): Not yet implemented!" << Exception::runerror;
-}
-
-Energy2 MatchboxTriVecScales::EtScale() const {
- throw Exception() << "MatchboxTriVecScales::EtScale(): Not yet implemented!" << Exception::runerror;
-}
-
-Energy2 MatchboxTriVecScales::renormalizationScale() const {
-// if ( theTriVecScaleChoice==1 ) return 0.5*HtPrimeScale();
-// else if ( theTriVecScaleChoice==2 ) return 0.5*HtPrimeModScale();
-// else if ( theTriVecScaleChoice==3 ) return 0.5*EtScale();
- if ( theTriVecScaleChoice==2 ) return 0.5*HtPrimeModScale();
- if ( theTriVecScaleChoice==3 ) return 0.5*EtScale();
- return 0.5*HtPrimeScale(); // The default is theTriVecScaleChoice==1
+ if ( theTriVecScaleChoice==1 || theTriVecScaleChoice==2 )
+ return sqr(0.5*(sumpartpt+sumvecet));
+ else if ( theTriVecScaleChoice==3 )
+ return sqr(0.5*sumvecet);
+ else
+ throw Exception() << "MatchboxTriVecScales: Scale choice out of range: "
+ << theTriVecScaleChoice
+ << Exception::runerror;
}
Energy2 MatchboxTriVecScales::factorizationScale() const {
return renormalizationScale();
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void MatchboxTriVecScales::persistentOutput(PersistentOStream & os) const {
os << theJetFinder << theTriVecScaleChoice;
}
void MatchboxTriVecScales::persistentInput(PersistentIStream & is, int) {
is >> theJetFinder >> theTriVecScaleChoice;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<MatchboxTriVecScales,MatchboxScaleChoice>
describeHerwigMatchboxTriVecScales("Herwig::MatchboxTriVecScales", "HwMatchboxScales.so");
void MatchboxTriVecScales::Init() {
static ClassDocumentation<MatchboxTriVecScales> documentation
("MatchboxTriVecScales implements scale choices related to transverse momenta and transverse energies for"
"events with up to three vector bosons in the final state plus additional jets, where the vector bosons"
"are associated to lepton pairs of distinct families.");
static Reference<MatchboxTriVecScales,JetFinder> interfaceJetFinder
("JetFinder",
"A reference to the jet finder.",
&MatchboxTriVecScales::theJetFinder, false, false, true, false, false);
static Switch<MatchboxTriVecScales,unsigned int> interfaceTriVecScaleChoice
("TriVecScaleChoice",
"The scale choice to use.",
&MatchboxTriVecScales::theTriVecScaleChoice, 1, false, false);
static SwitchOption interfaceTriVecScaleChoice1
(interfaceTriVecScaleChoice,
"HtPrimeScale",
"Sum of the transverse energies of the lepton pairs and the transverse momenta of the jets.",
1);
static SwitchOption interfaceTriVecScaleChoice2
(interfaceTriVecScaleChoice,
"HtPrimeModScale",
"Sum of the transverse energies of the lepton pairs and the transverse momenta of the jets."
"Each transverse jet momentum is thereby suppressed by an exponential of the rapidity difference to the average rapidity of the two hardest jets."
"This scale choice is defined only for two or more jets in the event.",
2);
static SwitchOption interfaceTriVecScaleChoice3
(interfaceTriVecScaleChoice,
"EtScale",
- "Sum of the transverse energies of the lepton pairs and the transverse energies of the jets."
- "This scale choice is defined only for two or more jets in the event.",
+ "Sum of the transverse energies of the lepton pairs.",
3);
}
diff --git a/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.h b/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.h
--- a/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.h
+++ b/MatrixElement/Matchbox/Scales/MatchboxTriVecScales.h
@@ -1,150 +1,135 @@
// -*- C++ -*-
//
// MatchboxTriVecScales.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_MatchboxTriVecScales_H
#define Herwig_MatchboxTriVecScales_H
//
// This is the declaration of the MatchboxTriVecScales class.
//
#include "Herwig/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
#include "ThePEG/PDT/MatcherBase.h"
#include "ThePEG/Cuts/JetFinder.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Christian Reuschle
*
* \brief MatchboxTriVecScales implements scale choices related to transverse momenta and transverse energies
* in events with up to three vector bosons in the final state, plus additional jets, where the vector bosons
* are associated to lepton pairs of distinct families.
*
*/
class MatchboxTriVecScales: public MatchboxScaleChoice {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxTriVecScales();
/**
* The destructor.
*/
virtual ~MatchboxTriVecScales();
//@}
public:
/**
- * Return the HtPrime scale.
- */
- virtual Energy2 HtPrimeScale() const;
-
- /**
- * Return the rapidity modified HtPrime scale.
- */
- virtual Energy2 HtPrimeModScale() const;
-
- /**
- * Return the Et scale.
- */
- virtual Energy2 EtScale() const;
-
- /**
* Return the renormalization scale. This default version returns
* shat.
*/
virtual Energy2 renormalizationScale() const;
/**
* Return the factorization scale. This default version returns
* shat.
*/
virtual Energy2 factorizationScale() 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 Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* Reference to the jet finder
*/
Ptr<JetFinder>::ptr theJetFinder;
/**
* The scale choice to use.
*/
unsigned int theTriVecScaleChoice;
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxTriVecScales & operator=(const MatchboxTriVecScales &);
};
}
#endif /* Herwig_MatchboxTriVecScales_H */
diff --git a/MatrixElement/Matchbox/Utility/MatchboxXCombData.h b/MatrixElement/Matchbox/Utility/MatchboxXCombData.h
--- a/MatrixElement/Matchbox/Utility/MatchboxXCombData.h
+++ b/MatrixElement/Matchbox/Utility/MatchboxXCombData.h
@@ -1,1106 +1,1112 @@
// -*- C++ -*-
//
// MatchboxXCombData.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_MatchboxXCombData_H
#define Herwig_MatchboxXCombData_H
//
// This is the declaration of the MatchboxXCombData class.
//
+// work around a Boost 1.64 bug where ublas headers would fail otherwise
+#include <boost/version.hpp>
+#if (BOOST_VERSION / 100 >= 1064)
+#include <boost/serialization/array_wrapper.hpp>
+#endif
+
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include "ThePEG/MatrixElement/MEBase.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/Persistency/PersistentOStream.fh"
#include "ThePEG/Persistency/PersistentIStream.fh"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define complex vector from boost::uBLAS
*/
typedef boost::numeric::ublas::vector<Complex> CVector;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define how amplitudes are stored
*/
typedef map<vector<int>,CVector> AmplitudeMap;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define amplitude iterators
*/
typedef map<vector<int>,CVector>::iterator AmplitudeIterator;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define amplitude const iterators
*/
typedef map<vector<int>,CVector>::const_iterator AmplitudeConstIterator;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Matchbox extensions to StandardXComb
*/
class MatchboxXCombData {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Standard constructor.
*/
explicit MatchboxXCombData(tMEPtr newME);
/**
* Default constructor.
*/
MatchboxXCombData();
/**
* Destructor.
*/
virtual ~MatchboxXCombData();
//@}
public:
/**
* Reset all cache flags
*/
void flushCaches();
public:
/**
* Get the factory
*/
Ptr<MatchboxFactory>::tcptr factory() const;
/**
* Get the matrix element; may return null
*/
Ptr<MatchboxMEBase>::tptr matchboxME() const;
/**
* Get the dipole; may return null
*/
Ptr<SubtractionDipole>::tptr subtractionDipole() const;
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
const vector<int>& crossingMap() const { return theCrossingMap; }
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
vector<int>& crossingMap() { return theCrossingMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
const map<size_t,size_t>& amplitudeToColourMap() const { return theAmplitudeToColourMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t>& amplitudeToColourMap() { return theAmplitudeToColourMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
const map<size_t,size_t>& colourToAmplitudeMap() const { return theColourToAmplitudeMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t>& colourToAmplitudeMap() { return theColourToAmplitudeMap; }
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
double crossingSign() const { return theCrossingSign; }
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
void crossingSign(double c) { theCrossingSign = c; }
/**
* The last renormalization scale
*/
Energy2 lastRenormalizationScale() const { return theLastRenormalizationScale; }
/**
* The last renormalization scale
*/
void lastRenormalizationScale(Energy2 lrs) { theLastRenormalizationScale = lrs; }
/**
* The amplitude parton data.
*/
const cPDVector& amplitudePartonData() const { return theAmplitudePartonData; }
/**
* The amplitude parton data.
*/
cPDVector& amplitudePartonData() { return theAmplitudePartonData; }
/**
* The crossed momenta
*/
const vector<Lorentz5Momentum>& amplitudeMomenta() const { return theAmplitudeMomenta; }
/**
* The crossed momenta
*/
vector<Lorentz5Momentum>& amplitudeMomenta() { return theAmplitudeMomenta; }
/**
* True, if the the tree level amplitudes need to be calculated
*/
bool calculateTreeAmplitudes() const { return theCalculateTreeAmplitudes; }
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastAmplitudes() const { return theLastAmplitudes; }
/**
* True, if the the tree level amplitudes need to be calculated
*/
void haveTreeAmplitudes(bool f = true) { theCalculateTreeAmplitudes = !f; }
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastAmplitudes() { return theLastAmplitudes; }
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastLargeNAmplitudes() const { return theLastLargeNAmplitudes; }
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastLargeNAmplitudes() { return theLastLargeNAmplitudes; }
/**
* True, if the the one-loop amplitudes need to be calculated
*/
bool calculateOneLoopAmplitudes() const { return theCalculateOneLoopAmplitudes; }
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastOneLoopAmplitudes() const { return theLastOneLoopAmplitudes; }
/**
* True, if the the one-loop amplitudes need to be calculated
*/
void haveOneLoopAmplitudes(bool f = true) { theCalculateOneLoopAmplitudes = !f; }
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastOneLoopAmplitudes() { return theLastOneLoopAmplitudes; }
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool calculateTreeME2() const { return theCalculateTreeME2; }
/**
* The last tree-level matrix element squared
*/
double lastTreeME2() const { return theLastTreeME2; }
/**
* The last tree-level matrix element squared
*/
void lastTreeME2(double v) {
theLastTreeME2 = v; theCalculateTreeME2 = false;
}
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool calculateLargeNME2() const { return theCalculateLargeNME2; }
/**
* The last tree-level matrix element squared
*/
double lastLargeNME2() const { return theLastLargeNME2; }
/**
* The last tree-level matrix element squared
*/
void lastLargeNME2(double v) {
theLastLargeNME2 = v; theCalculateLargeNME2 = false;
}
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool calculateOneLoopInterference() const { return theCalculateOneLoopInterference; }
/**
* The last one-loop/tree-level interference.
*/
double lastOneLoopInterference() const { return theLastOneLoopInterference; }
/**
* The last one-loop/tree-level interference.
*/
void lastOneLoopInterference(double v) {
theLastOneLoopInterference = v; theCalculateOneLoopInterference = false;
}
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool calculateOneLoopPoles() const { return theCalculateOneLoopPoles; }
/**
* The last one-loop/tree-level interference.
*/
pair<double,double> lastOneLoopPoles() const { return theLastOneLoopPoles; }
/**
* The last one-loop/tree-level interference.
*/
void lastOneLoopPoles(pair<double,double> v) {
theLastOneLoopPoles = v; theCalculateOneLoopPoles = false;
}
/**
* True, if the indexed colour correlated matrix element needs to be
* calculated.
*/
bool calculateColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,bool>::const_iterator f =
theCalculateColourCorrelators.find(ij);
if ( f == theCalculateColourCorrelators.end() )
return true;
return f->second;
}
/**
* The colour correlated matrix element.
*/
double lastColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,double>::const_iterator v =
theColourCorrelators.find(ij);
if ( v == theColourCorrelators.end() )
return 0.;
return v->second;
}
/**
* The colour correlated matrix element.
*/
void lastColourCorrelator(pair<int,int> ij, double v) {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
theColourCorrelators[ij] = v;
theCalculateColourCorrelators[ij] = false;
}
/**
* True, if the indexed large-N colour correlated matrix element needs to be
* calculated.
*/
bool calculateLargeNColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,bool>::const_iterator f =
theCalculateLargeNColourCorrelators.find(ij);
if ( f == theCalculateLargeNColourCorrelators.end() )
return true;
return f->second;
}
/**
* The large-N colour correlated matrix element.
*/
double lastLargeNColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,double>::const_iterator v =
theLargeNColourCorrelators.find(ij);
if ( v == theLargeNColourCorrelators.end() )
return 0.;
return v->second;
}
/**
* The large-N colour correlated matrix element.
*/
void lastLargeNColourCorrelator(pair<int,int> ij, double v) {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
theLargeNColourCorrelators[ij] = v;
theCalculateLargeNColourCorrelators[ij] = false;
}
/**
* True, if the indexed colour/spin correlated matrix element needs to be
* calculated.
*/
bool calculateColourSpinCorrelator(const pair<int,int>& ij) const {
map<pair<int,int>,bool>::const_iterator f =
theCalculateColourSpinCorrelators.find(ij);
if ( f == theCalculateColourSpinCorrelators.end() )
return true;
return f->second;
}
/**
* The colour/spin correlated matrix element.
*/
Complex lastColourSpinCorrelator(const pair<int,int>& ij) const {
map<pair<int,int>,Complex>::const_iterator v =
theColourSpinCorrelators.find(ij);
if ( v == theColourSpinCorrelators.end() )
return 0.;
return v->second;
}
/**
* The colour/spin correlated matrix element.
*/
void lastColourSpinCorrelator(const pair<int,int>& ij, Complex v) {
theColourSpinCorrelators[ij] = v;
theCalculateColourSpinCorrelators[ij] = false;
}
/**
* True, if the indexed spin correlated matrix element needs to be
* calculated.
*/
bool calculateSpinCorrelator(const pair<int,int>& ij) const {
map<pair<int,int>,bool>::const_iterator f =
theCalculateSpinCorrelators.find(ij);
if ( f == theCalculateSpinCorrelators.end() )
return true;
return f->second;
}
/**
* The spin correlated matrix element.
*/
Complex lastSpinCorrelator(const pair<int,int>& ij) const {
map<pair<int,int>,Complex>::const_iterator v =
theSpinCorrelators.find(ij);
if ( v == theSpinCorrelators.end() )
return 0.;
return v->second;
}
/**
* The spin correlated matrix element.
*/
void lastSpinCorrelator(const pair<int,int>& ij, Complex v) {
theSpinCorrelators[ij] = v;
theCalculateSpinCorrelators[ij] = false;
}
/**
* Return the number of light flavours to be considered for this process.
*/
unsigned int nLight() const { return theNLight; }
/**
* Set the number of light flavours to be considered for this process.
*/
void nLight(unsigned int n) { theNLight = n; }
/**
* Return the vector that contains the PDG ids of
* the light flavours, which are contained in the
* jet particle group.
*/
vector<long> nLightJetVec() const { return theNLightJetVec; }
/**
* Set the elements of the vector that contains the PDG
* ids of the light flavours, which are contained in the
* jet particle group.
*/
void nLightJetVec(long n) { theNLightJetVec.push_back(n); }
/**
* Return the vector that contains the PDG ids of
* the heavy flavours, which are contained in the
* jet particle group.
*/
vector<long> nHeavyJetVec() const { return theNHeavyJetVec; }
/**
* Set the elements of the vector that contains the PDG
* ids of the heavy flavours, which are contained in the
* jet particle group.
*/
void nHeavyJetVec(long n) { theNHeavyJetVec.push_back(n); }
/**
* Return the vector that contains the PDG ids of
* the light flavours, which are contained in the
* proton particle group.
*/
vector<long> nLightProtonVec() const { return theNLightProtonVec; }
/**
* Set the elements of the vector that contains the PDG
* ids of the light flavours, which are contained in the
* proton particle group.
*/
void nLightProtonVec(long n) { theNLightProtonVec.push_back(n); }
/**
* Get the dimensionality of the colour basis for this process.
*/
size_t colourBasisDim() const { return theColourBasisDim; }
/**
* Set the dimensionality of the colour basis for this process.
*/
void colourBasisDim(size_t d) { theColourBasisDim = d; }
/**
* Return the number of degrees of freedom required by the phase space generator
*/
int nDimPhasespace() const { return theNDimPhasespace; }
/**
* Set the number of degrees of freedom required by the phase space generator
*/
void nDimPhasespace(int d) { theNDimPhasespace = d; }
/**
* Return the number of degrees of freedom required by the amplitude
*/
int nDimAmplitude() const { return theNDimAmplitude; }
/**
* Set the number of degrees of freedom required by the amplitude
*/
void nDimAmplitude(int d) { theNDimAmplitude = d; }
/**
* Return the number of degrees of freedom required by the insertion operators
*/
int nDimInsertions() const { return theNDimInsertions; }
/**
* Set the number of degrees of freedom required by the insertion operators
*/
void nDimInsertions(int d) { theNDimInsertions = d; }
/**
* Get the additional random numbers required by the amplitude
*/
const vector<double>& amplitudeRandomNumbers() const { return theAmplitudeRandomNumbers; }
/**
* Access the additional random numbers required by the amplitude
*/
vector<double>& amplitudeRandomNumbers() { return theAmplitudeRandomNumbers; }
/**
* Get the additional random numbers required by the insertion operator
*/
const vector<double>& insertionRandomNumbers() const { return theInsertionRandomNumbers; }
/**
* Access the additional random numbers required by the insertion operator
*/
vector<double>& insertionRandomNumbers() { return theInsertionRandomNumbers; }
/**
* Return the diagram weights indexed by diagram id.
*/
const map<int,double>& diagramWeights() const { return theDiagramWeights; }
/**
* Access the diagram weights indexed by diagram id.
*/
map<int,double>& diagramWeights() { return theDiagramWeights; }
/**
* Return the singular limits
*/
const set<pair<size_t,size_t> >& singularLimits() const { return theSingularLimits; }
/**
* Access the singular limits
*/
set<pair<size_t,size_t> >& singularLimits() { return theSingularLimits; }
/**
* Return the last matched singular limit.
*/
const set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() const { return theLastSingularLimit; }
/**
* Access the last matched singular limit.
*/
set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() { return theLastSingularLimit; }
/**
* Set the Herwig StandardModel object
*/
void hwStandardModel(Ptr<StandardModel>::tcptr sm) { theStandardModel = sm; }
/**
* Get the Herwig StandardModel object
*/
Ptr<StandardModel>::tcptr hwStandardModel() const { return theStandardModel; }
/**
* Return the symmetry factor
*/
double symmetryFactor() const { return theSymmetryFactor; }
/**
* Set the symmetry factor
*/
void symmetryFactor(double f) { theSymmetryFactor = f; }
/**
* Return the OLP process ids
*/
const vector<int>& olpId() const { return theOLPId; }
/**
* Set the OLP process ids
*/
void olpId(int pType, int id) {
if ( theOLPId.empty() )
theOLPId.resize(4,0);
theOLPId[pType] = id;
}
/**
* Set the OLP process ids
*/
void olpId(const vector<int>& id) {
theOLPId = id;
}
/**
* Return the olp momentum vector
*/
double* olpMomenta() { return theOLPMomenta; }
/**
* Fill the olp momentum vector
*/
void fillOLPMomenta(const vector<Lorentz5Momentum>& mm,
const cPDVector& mePartonData = cPDVector(),
const map<long,Energy>& reshuffleMap = map<long,Energy>());
/**
* Helper struct to define the reshuffling equation
*/
struct ReshuffleEquation {
ReshuffleEquation(Energy xq,
cPDVector::const_iterator xdBegin,
cPDVector::const_iterator xdEnd,
vector<Lorentz5Momentum>::const_iterator xmBegin,
const map<long,Energy>* xreshuffleMap)
: q(xq), dBegin(xdBegin), dEnd(xdEnd), mBegin(xmBegin),
reshuffleMap(xreshuffleMap) {}
typedef double ArgType;
typedef double ValType;
static double aUnit() { return 1.; }
static double vUnit() { return 1.; }
double operator() (double xi) const;
Energy q;
cPDVector::const_iterator dBegin;
cPDVector::const_iterator dEnd;
vector<Lorentz5Momentum>::const_iterator mBegin;
const map<long,Energy>* reshuffleMap;
};
/**
* Perform a reshuffling to the mass shells contained in the map
*/
void reshuffle(vector<Lorentz5Momentum>& momenta,
const cPDVector& mePartonData,
const map<long,Energy>& reshuffleMap) const;
/**
* Return a generic process id to communicate with external codes
*/
int externalId() const { return theExternalId; }
/**
* Set a generic process id to communicate with external codes
*/
void externalId(int id) { theExternalId = id; }
/**
* True if the process has been initialized
*/
bool initialized() const { return theInitialized; }
/**
* True if the process has been initialized
*/
void isInitialized(bool is = true) { theInitialized = is; }
/**
* Return the external momentum vector
*/
const vector<double*>& externalMomenta() const { return theExternalMomenta; }
/**
* Fill the external momentum vector
*/
void fillExternalMomenta(const vector<Lorentz5Momentum>&);
/**
* Caching for the external madgraph colour structures
*/
const map<vector<int>,vector < complex<double> > >& heljamp() const { return theHelJamp; }
/**
* Caching for the external madgraph colour structures (large N)
*/
const map<vector<int>,vector < complex<double> > >& helLNjamp() const { return theLNHelJamp; }
/**
* pushback the madgraph colour structures
*/
void pushheljamp(const vector<int>& hel, const complex<double>& jamp) { theHelJamp[hel].push_back(jamp); }
/**
* clear the madgraph colour structures
*/
void clearheljamp() { theHelJamp.clear(); }
/**
* pushback the madgraph colour structures (large N)
*/
void pushhelLNjamp(const vector<int>& hel, const complex<double>& jamp) { theLNHelJamp[hel].push_back(jamp); }
/**
* clear the madgraph colour structures (large N)
*/
void clearhelLNjamp() { theLNHelJamp.clear(); }
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);
/**
* Put a CVector to the persistent ostream
*/
static void putCVector(PersistentOStream&, const CVector&);
/**
* Get a CVector from the persistent istream
*/
static void getCVector(PersistentIStream&, CVector&);
/**
* Put an amplitude map to the persistent ostream
*/
static void putAmplitudeMap(PersistentOStream&, const map<vector<int>,CVector>&);
/**
* Get an amplitude map from the persistent istream
*/
static void getAmplitudeMap(PersistentIStream&, map<vector<int>,CVector>&);
//@}
/**
* 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();
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxXCombData & operator=(const MatchboxXCombData &);
private:
/**
* The factory
*/
Ptr<MatchboxFactory>::tcptr theFactory;
/**
* The matrix element
*/
Ptr<MatchboxMEBase>::tptr theMatchboxME;
/**
* The dipole
*/
Ptr<SubtractionDipole>::tptr theSubtractionDipole;
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
vector<int> theCrossingMap;
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t> theAmplitudeToColourMap;
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t> theColourToAmplitudeMap;
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
double theCrossingSign;
/**
* The last renormalization scale
*/
Energy2 theLastRenormalizationScale;
/**
* The amplitude parton data.
*/
cPDVector theAmplitudePartonData;
/**
* The ccrossed momenta
*/
vector<Lorentz5Momentum> theAmplitudeMomenta;
/**
* True, if the the tree level amplitudes need to be calculated
*/
bool theCalculateTreeAmplitudes;
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastAmplitudes;
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastLargeNAmplitudes;
/**
* True, if the the one-loop amplitudes need to be calculated
*/
bool theCalculateOneLoopAmplitudes;
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastOneLoopAmplitudes;
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool theCalculateTreeME2;
/**
* The last tree-level matrix element squared
*/
double theLastTreeME2;
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool theCalculateLargeNME2;
/**
* The last tree-level matrix element squared
*/
double theLastLargeNME2;
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool theCalculateOneLoopInterference;
/**
* The last one-loop/tree-level interference.
*/
double theLastOneLoopInterference;
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool theCalculateOneLoopPoles;
/**
* The last one-loop/tree-level interference.
*/
pair<double,double> theLastOneLoopPoles;
/**
* True, if the indexed colour correlated matrix element needs to be
* calculated.
*/
map<pair<int,int>,bool> theCalculateColourCorrelators;
/**
* The colour correlated matrix element.
*/
map<pair<int,int>,double> theColourCorrelators;
/**
* True, if the indexed large-N colour correlated matrix element needs to be
* calculated.
*/
map<pair<int,int>,bool> theCalculateLargeNColourCorrelators;
/**
* The large-N colour correlated matrix element.
*/
map<pair<int,int>,double> theLargeNColourCorrelators;
/**
* True, if the indexed colour/spin correlated matrix element needs to be
* calculated.
*/
map<pair<int,int>,bool> theCalculateColourSpinCorrelators;
/**
* The colour/spin correlated matrix element.
*/
map<pair<int,int>,Complex> theColourSpinCorrelators;
/**
* True, if the indexed spin correlated matrix element needs to be
* calculated.
*/
map<pair<int,int>,bool> theCalculateSpinCorrelators;
/**
* The spin correlated matrix element.
*/
map<pair<int,int>,Complex> theSpinCorrelators;
/**
* The number of light flavours to be considered for this process.
*/
static unsigned int theNLight;
/**
* Vector with the PDG ids of the light quark flavours,
* which are contained in the jet particle group.
*/
static vector<long> theNLightJetVec;
/**
* Vector with the PDG ids of the heavy quark flavours,
* which are contained in the jet particle group.
*/
static vector<long> theNHeavyJetVec;
/**
* Vector with the PDG ids of the light quark flavours,
* which are contained in the proton particle group.
*/
static vector<long> theNLightProtonVec;
/**
* The dimensionality of the colour basis for this process.
*/
size_t theColourBasisDim;
/**
* The number of degrees of freedom required by the phase space generator
*/
int theNDimPhasespace;
/**
* The number of degrees of freedom required by the amplitude
*/
int theNDimAmplitude;
/**
* The number of degrees of freedom required by the insertion operators
*/
int theNDimInsertions;
/**
* Additional random numbers required by the amplitude
*/
vector<double> theAmplitudeRandomNumbers;
/**
* Additional random numbers required by the insertion operator
*/
vector<double> theInsertionRandomNumbers;
/**
* The diagram weights indexed by diagram id.
*/
map<int,double> theDiagramWeights;
/**
* If not empty, the entries here serve to limit phasespace
* generation to the corresponding collinear limits, or soft limits
* if both pair entries are equal.
*/
set<pair<size_t,size_t> > theSingularLimits;
/**
* The last matched singular limit.
*/
set<pair<size_t,size_t> >::const_iterator theLastSingularLimit;
/**
* The Herwig StandardModel object
*/
Ptr<StandardModel>::tcptr theStandardModel;
/**
* The symmetry factor
*/
double theSymmetryFactor;
/**
* The OLP process id
*/
vector<int> theOLPId;
/**
* Return the olp momentum vector
*/
double* theOLPMomenta;
/**
* True, if olp momenta have been filled
*/
bool filledOLPMomenta;
/**
* A generic process id to communicate with external codes
*/
int theExternalId;
/**
* True if the process has been initialized
*/
bool theInitialized;
/**
* The external momenta
*/
vector<double*> theExternalMomenta;
/**
* True, if external momenta have been filled
*/
bool filledExternalMomenta;
/**
* caching of different colour structures (MadGraph-Interface)
*/
map<vector<int>,vector < complex<double> > > theHelJamp;
/**
* caching of different colour structures (MadGraph-Interface)
*/
map<vector<int>,vector < complex<double> > > theLNHelJamp;
};
}
#endif /* Herwig_MatchboxXCombData_H */
diff --git a/Models/Makefile.am b/Models/Makefile.am
--- a/Models/Makefile.am
+++ b/Models/Makefile.am
@@ -1,181 +1,181 @@
SUBDIRS = RSModel StandardModel General Susy UED Zprime \
Transplanckian ADD Leptoquarks Sextet TTbAsymm \
LH LHTP Feynrules
noinst_LTLIBRARIES = libHwStandardModel.la
nodist_libHwStandardModel_la_SOURCES = \
StandardModel/SM__all.cc
libHwStandardModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/StandardModel
noinst_LTLIBRARIES += libHwModelGenerator.la
nodist_libHwModelGenerator_la_SOURCES = \
General/ModelGenerator__all.cc
libHwModelGenerator_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/General
if WANT_BSM
pkglib_LTLIBRARIES = \
HwRSModel.la \
HwUED.la \
HwSusy.la \
HwNMSSM.la \
HwRPV.la \
HwZprimeModel.la \
HwTransplanck.la \
HwADDModel.la \
HwLeptoquarkModel.la \
HwSextetModel.la \
HwTTbAModel.la \
HwLHModel.la \
HwLHTPModel.la
endif
#############
HwRSModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwRSModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/RSModel
nodist_HwRSModel_la_SOURCES = \
RSModel/RS__all.cc
#############
HwUED_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwUED_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/UED
nodist_HwUED_la_SOURCES = \
UED/UED__all.cc
#############
HwSusy_la_LDFLAGS = \
-$(AM_LDFLAGS) -module -version-info 12:0:0
+$(AM_LDFLAGS) -module -version-info 13:0:0
HwSusy_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Susy
nodist_HwSusy_la_SOURCES = \
Susy/Susy__all.cc
#############
HwNMSSM_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 5:0:0
HwNMSSM_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Susy/NMSSM
nodist_HwNMSSM_la_SOURCES = \
Susy/NMSSM/NMSSM__all.cc
#############
HwRPV_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 3:0:0
HwRPV_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Susy/RPV
nodist_HwRPV_la_SOURCES = \
Susy/RPV/RPV__all.cc
#############
HwZprimeModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 3:0:0
HwZprimeModel_la_SOURCES = \
Zprime/ZprimeModel.cc Zprime/ZprimeModelZPQQVertex.cc
#############
HwTransplanck_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 4:0:0
HwTransplanck_la_SOURCES = \
Transplanckian/METRP2to2.cc
#############
HwADDModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 4:0:0
HwADDModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/ADD
nodist_HwADDModel_la_SOURCES = \
ADD/ADD__all.cc
#############
HwLeptoquarkModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 5:0:0
HwLeptoquarkModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Leptoquarks
nodist_HwLeptoquarkModel_la_SOURCES = \
Leptoquarks/Leptoquark__all.cc
#############
HwSextetModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 3:0:0
HwSextetModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Sextet
nodist_HwSextetModel_la_SOURCES = \
Sextet/Sextet__all.cc
#############
HwTTbAModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 3:0:0
HwTTbAModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/TTbAsymm
nodist_HwTTbAModel_la_SOURCES = \
TTbAsymm/TTbA__all.cc
#############
HwLHModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 5:0:0
HwLHModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/LH
nodist_HwLHModel_la_SOURCES = \
LH/LH__all.cc
#############
HwLHTPModel_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 5:0:0
HwLHTPModel_la_LIBADD = \
$(GSLLIBS)
HwLHTPModel_la_CPPFLAGS = \
$(AM_CPPFLAGS) $(GSLINCLUDE) -I$(srcdir)/LHTP
nodist_HwLHTPModel_la_SOURCES = \
LHTP/LHTP__all.cc
#############
diff --git a/NEWS b/NEWS
--- a/NEWS
+++ b/NEWS
@@ -1,1623 +1,1644 @@
Herwig News -*- outline -*-
================================================================================
+* Herwig 7.1.1 release: 2017-07-14
+
+** Snippets are now all installed
+
+** Fixed broken ufo2herwig output and LHC-MB.in
+
+** UFO improvements
+*** More robust SLHA file handling
+*** option of creating diagonal mixing matrices, needed for ATLAS simplfied models
+*** Improved warnings about resetting standard model particles
+*** Fixed certain cases where the wrong lorentz structure was picked in VVS vertices
+
+** Improved error message for unhandled beam particles
+
+** Fix for Dipole Shower chain selection
+
+** Fixed crash in double diffractive delta resonances
+
+
+
* Herwig 7.1.0 release: 2017-05-19
** Major new release
- For a more detailed overview and further references please see the release note arXiv:1705.06919
+ For a more detailed overview and further references please see the
+ release note arXiv:1705.06919
** NLO multijet merging with the dipole shower
** A new soft model
** An interface to EvtGen
** Improved calculation of mass effects in the dipole shower
** Top decays in the dipole shower, and NLO corrections to the decay
** An implementation of the KrkNLO method for simple processes
** Major restructuring and cleanup of default input files
** C++11 is now mandatory and heavily used in the code
** Many smaller bugfixes and improvements
* Herwig 7.0.4 release: 2016-10-24
** API
The high level API is now properly available as a library, providing an
alternative to the Herwig main program.
** Dipole shower
Added nloops() function to the AlphaS classes to return the number of loops
in the running coupling
** Matchbox
Improved error handling and clearer messages
** BSM models
Initialize W mass correctly from SLHA file. Improved reading in of decay
modes if they already exist.
** Sampling
Introduced option to reduce reference weight in AlmostUnweighted mode by
Kappa factor. Useful for processes where full unweighting is infeasible.
** Qtilde shower
Better control of scale in splitting functions using the
SplittingFunction:ScaleChoice interface.
** Tests
New NLO fixed-order input files for testing Matchbox in Tests/ExternNLO.
** Input files
Set diagonal CKM options consistently. Added TauTauHVertex and MuMuHVertex
to MatchboxDefaults.
* Herwig 7.0.3 release: 2016-07-21
** subprocess generation filters
A number of filters has been implemented to speed up process and diagram
generation for the case of Standard Model like processes; this is
particularly helpful for processes involving a large number of electroweak
combinatorics. See the documentation for more details.
** fix to directory hierarchy
A bug in the Herwig-scratch directory hierarchy when integrating with
different seeds has been fixed.
** fix for relocating MadGraph generated amplitudes
The directory and symbolic link structure for MadGraph-built amplitudes has
been changed to allow for relocation of a run directory to a different file
system hierarchy.
** z* variable added to LeptonsJetsAnalysis
The builtin LeptonsJetsAnalysis has been extented to include the normalized
relative rapidity z*
** detuning parameter for shower reweighting
An efficiency tweak when using the shower reweighting facilities has been
introduced for the veto algorithms in both the QTilde and the Dipole
shower. See the documentation for more details.
** BaryonThreeQuarkModelFormFactor
A scaling bug in the form factor series expansion was fixed.
* Herwig 7.0.2 release: 2016-04-29
** Event reweighting
New option of calculating the weights for the effect of varying the scale
used in both the default and dipole showers. The method is described in
arXiv:1605.08256
** mergegrids mode
A main mode for the Herwig executable has been added which merges the
integration grids from parallel runs.
** NLO Diphoton production
The NLO matrix elements in the POWHEG approach for diphoton production as
described in JHEP 1202 (2012) 130, arXiv:1106.3939 have been included as
MEPP2GammaGammaPowheg.
** BSM Hard process constructor
A missing Feynman diagram with a t-channel vector boson has been added to
the matrix element for vv2ss
** BSM Decay Modes calculation
The behaviour of the option to disable decay modes in BSM Models has been
changed so that the partial width for the ignored modes is now calculated
and included in the total width, but the branching ratio is set to
zero. This is more physical than the previous option where the mode was
t otally ignored and hence not included in the calculation of the width.
** Mass Generation
The behaviour of the GenericMassGenerator has been changed so that modes
which have been disabled are only included in the calculation of the total
width and not in the partial width used in the numerator of the weight used
to select the off-shell mass.
** Boost detection
Boost could not detect the compiler version for gcc-5.3 and gcc-6.1
* Herwig 7.0.1 release: 2016-01-17
** Version number written to log file
The Herwig version number is now included in addition to ThePEG's version.
** Tau lifetimes
A bug with lifetimes for hard process Taus is fixed. Thanks to ATLAS
for the report!
** Shower FSR retries
Rare events could take a long time due to an extremely large number
of FSR retries. These are now capped at a configurable number.
** Dipole shower
Reweighting for non-radiating events; fixes for shower profile
handling; option to downgrade large-Q expansion for alphaS
** Matchbox builtins
Added massive currents for q-qbar
** ShowerAlphaQCD
can now use user-defined thresholds
** Input snippets
W/Z/H on-shell now split into three files; 4-flavour scheme added
** UFO converter
The converter now has experimental support for writing out param
cards of its current settings.
** LEPJetAnalysis loading fixed
** Contrib
HJets++ has moved to a stand-alone project, FxFx has been added
* Herwig 7.0.0 (Herwig++ 3.0.0) release: 2015-12-04
** Major new release
A major new release of the Monte Carlo event generator Herwig++ (version
3.0) is now available. This release marks the end of distinguishing
Herwig++ and HERWIG development and therefore constitutes the first major
release of version 7 of the Herwig event generator family. The new version
features a number of significant improvements to the event simulation,
including: built-in NLO hard process calculation for all Standard Model
processes, with matching to both angular ordered and dipole shower modules
via variants of both subtractive (MC@NLO-type) and multiplicative
(Powheg-type) algorithms; QED radiation and spin correlations in the
angular ordered shower; a consistent treatment of perturbative
uncertainties within the hard process and parton showering, as well as a
vastly improved documentation. This version includes (for a more detailed
overview and further references please see the release note
arXiv:1512.01178):
** A long list of improvements and fixes for the Matchbox module
*** includes MC@NLO and Powheg matching to both showers with truncated showering
** A long list of improvements and fixes for both of the shower modules
*** includes improvements of numerics issues relevant to 100 TeV pp collisions
** NLO event simulation and Matchbox development
*** Interfaces to a number of external libraries
*** A new workflow for event generation
*** Electroweak corrections to VV production
** Parton shower development
*** QED radiation in the angular ordered shower
*** Spin correlations in the angular ordered shower
*** New scale choices in gluon branchings
** Improvements to event generation workflow
*** Re-organized and streamlined input files for the new NLO development
*** A unified treatment of shower and matching uncertainties
*** New integrator modules featuring parallel integration
** New default tunes for both shower modules
** New contrib modules
*** Electroweak Higgs plus jets production
*** FxFx merging support
*** Higgs boson pair production
* Herwig++-2.7.1 release: 2014-07-07
** New shower switches to select schemes for momentum reconstruction
*** QTildeReconstructor:FinalStateReconOption
has the following options:
*** Default
All momenta are rescaled in the rest frame.
*** MostOffShell
Put all particles on the new-mass shell and the most off-shell and
recoiling system are rescaled to ensure 4-momentum is conserved.
*** Recursive
Recursively put the most off-shell particle which hasn't yet been
rescaled on-shell by rescaling the particles and the recoiling
system.
*** RestMostOffShell
The most off-shell is put on shell by rescaling it and the
recoiling system, the recoiling system is then put on-shell in its
rest frame.
*** RestRecursive
As above, but recursively treat the currently most-off shell
(only makes a difference for more than 3 partons)
** Ticket #378: Hadronization of baryon number violating clusters involving diquarks
Fixed by only considering non-diquarks to be combined in the
ClusterFinder.
** UFO converter can now parse SLHA files for parameter settings
The UFO converter code can now use SLHA files for modifying
parameters. The first pass "ufo2herwig" produces the model to be
compiled. For each parameter card, run "slha2herwig" to get the
matching input file.
** Fix for systems using lib64
The repository is now initialized correctly on systems using lib64
as the library location.
** Efficiency optimization
Better allocation of internal vector variables for a noticeable
speed increase of 10-20% with LHC events.
* Herwig++-2.7.0 release: 2013-10-28
** UFO interface to Feynman rules generators
Herwig++ now includes "ufo2herwig", a tool that automatically
creates all required files to run a BSM model from a UFO
directory. The conversion has been extensively tested against
Feynrules models MSSM, NMSSM, RS, Technicolor,
and less extensively with most of the
other models in the Feynrules model database.
We expect that following this release there will be no further
hard-coded new physics models added to Herwig++ and that future
models will be included using the UFO interface.
** Shower uncertainties
A first set of scaling parameters to estimate shower uncertainties
is provided for both the angular ordered as well as the dipole
shower; they are Evolver:HardScaleFactor and ShowerAlphaQCD:
RenormalizationScaleFactor.
** Rewrite of Matchbox NLO matching
The NLO matching implementation has been rewritten and is now more
flexible and consistent. Profile scales are provided for the
hardest emission both for the dipole shower and matrix element
correction matching.
** BLHA2 Interface and new processes
Matchbox now features a generic BLHA2 interface to one-loop
amplitude codes and now also includes W and W+jet production as
well as Higss production in gluon fusion as builtin processes.
** Impoved dipole shower kinematics parametrization
The kinematics parametrization for emissions in the dipole shower
has been made more efficient.
** W and Z Powheg decays
Decays of W and Z bosons now use the Powheg decayers by default.
** Improved treatment of beam remnants
The handling of beam remnants has been improved in multiple
contexts, leading to a much lower error rate at p+/p- collisions.
An additional value "VeryHard" for ClusterFissioner:RemnantOption
can be used to disable any special treatment of beam remnant
clusters.
** New underlying event tune
Herwig++ now uses tune UE-EE-5-MRST by default. Other related tunes
can be obtained from the Herwig++ tunes page
** Improvements in BSM code
The UFO development identified many sign fixes in rarely used BSM
vertices; many improvements were made to general decayers, allowing
four-body decays in BSM for the first time; Powheg is enabled in
General two-body decayers; and the handling of colour
sextets has been improved.
** A new HiggsPair matrix element in Contrib.
** A new matrix element for single top production.
** The Higgs mass is now set to 125.9 GeV (from PDG 2013 update).
** C++-11 testing
To help with the coming transition to C++-11, we provide the new
--enable-stdcxx11 configure flag. Please try to test builds with
this flag enabled and let us know any problems, but do not use this
in production code yet. In future releases, this flag will be on by
default.
** Other changes
*** Many new Rivet analyses have been included in the Tests directory.
*** Cleaned Shower header structure; grouped shower parameters into one struct.
*** The boolean Powheg flag in HwMEBase changed to an enumeration.
* Herwig++-2.6.3 release: 2013-02-22
** Decay vertex positioning in HepMC output
Pseudo-vertices that Herwig++ inserts for technical reasons will
now not contribute to the Lorentz positions of downstream vertices.
Thanks to ATLAS for the bug report!
** Updated Rivet tests
Herwig's library of Rivet test runs has been brought up-to-date
with new analyses that were recently published by the Rivet
collaboration.
* Herwig++-2.6.2 release: 2013-01-30
** Fixes for PDF and scale choices in POWHEG events
Scale settings for MPI and the regular shower are now correct in
POWHEG events. This should fix reported anomalies in POWHEG jet rates.
NLO PDFs are now also set consistently in the example input files.
** Ticket #373: Branching ratio factors in cross-section
If any decay modes are selectively disabled, setting the following
post-handler will cause all reported cross-sections to include the
branching ratio factor(s) from the previous stages correctly:
create Herwig::BranchingRatioReweighter BRreweight
insert LHCGenerator:EventHandler:PostDecayHandlers 0 BRreweight
** Anomalous vertices now possible in MEfftoVH
** Interactive shell does not quit on error
** Better warning messages for events from inconsistent LHEF files
** Possible division by zero error fixed in BSM branching ratio calculations
** Decayer and ME changes to improve checkpointing
The checkpointing changes in ThePEG 1.8.2 are implemented here, too. Regular
dump files are consistent now.
* Herwig++-2.6.1 release: 2012-10-16
** Configure switches
The various switches to turn off compilation of BSM models have
been unified into a single '--disable-models'. A new flag
'--disable-dipole' can be used to turn off the compilation of the
Dipole and Matchbox codes.
** Ticket #348: Search path for repository 'read' command
The search path for the 'read' command is configurable on the
command line with the -i and -I switches. By default, the
installation location is now included in the search path, so that
'Herwig++ read LEP.in' will work in an empty directory. The current
working directory will always be searched first.
The rarely used "Herwig++ init" command has been made consistent
with 'read' and 'run' and should now be used without the '-i' flag.
** Width treatment in BSM
The width treatment in BSM decay chains has been greatly improved
and is now switched on by default in the .model files. To get the
old behaviour, use
set /Herwig/NewPhysics/NewModel:WhichOffshell Selected
** New BSM models
Little Higgs models with and without T-parity are now available.
** Resonance photon lifetime
A lifetime bug affecting decays of pi0 to e+e-X was fixed. The
virtual photon is not part of the event record anymore.
** Ticket #371: Hard diffraction FPE
Herwig++ 2.6.0 introduced a bug into the diffraction code which
would abort any runs. This is now fixed.
** O2AlphaS
Support for setting quark masses different from the particle data
objects as introduced in ThePEG 1.8.1 has been enabled.
** Matchbox
Several improvements and bug fixes are included for
Matchbox. Amplitudes relevant to pp -> Z+jet and crossed processes
at NLO are now available, and various scale choices have been added
in a more flexible way. All subtraction dipoles for massive quarks
are now included.
** Dipole shower
Parameters to perform scale variations in the shower have been
added to estimate uncertainties. A bug in showering off gg -> h has
been fixed.
** Minor fixes
*** Two broken colour structures in GeneralHardME
*** Susy Higgs mixing matrix
*** BaryonFactorizedDecayer out-of-bounds access
*** Mass values in SimpleLHCAnalysis
* Herwig++-2.6.0 release: 2012-05-21 (tagged at SVN r7407)
** New NLO framework
Matchbox, a flexible and very general framework for performing NLO
calculations at fixed order or matched to parton showers is
provided with this release.
** Dipole shower algorithm
A first implementation of the coherent dipole shower algorithm by
Plätzer and Gieseke (arXiv:0909.5593 and arXiv:1109.6256) is
available.
** Alternative samplers and the ExSample library
The ExSample library by Plätzer (arXiv:1108.6182) is shipped along
with Herwig++ in an extended version. The extended version provides
SamplerBase objects which can be used alternatively to the default
ACDCSampler.
** New BSM models
*** New colour sextet diquark model
A colour sextet diquark model has been included, as described in
Richardson and Winn (arXiv:1108.6154).
*** Models reproducing the CDF t-tbar asymmetry
Four models that can reproduce the reported t-tbar asymmetry have
been included.
*** Zprime
A simple standard model extension by one additional heavy neutral
vector boson.
** Interface to AlpGen, with MLM merging
The Contrib directory contains a new interface to the AlpGen matrix
element generator. AlpGen events must be preprocessed with the
provided AlpGenToLH.exe tool before they can be used. More
information can be found in the Herwig++ 2.6 release note.
** HiggsVBF Powheg
Higgs boson production by vector boson fusion is available at NLO
in the POWHEG scheme, as described in d'Errico, Richardson
(arXiv:1106.2983). The Powheg DIS processes were available in
Herwig++-2.5.2 already.
** Statistical colour reconnection
Alternative mechanisms to minimize the colour length Sum(m_clu)
before the hadronization stage, based on Metropolis and annealing
algorithms.
** Energy extrapolation of underlying-event tunes
To describe underlying-event data at different c.m. energies, the
energy-dependent parameter pT_min will now be adjusted
automatically, following a power-law. The new tune parameters are
the value at 7000 GeV "MPIHandler:pTmin0", and MPIHandler:Power.
** Ticket #239: Reporting of minimum-bias cross-section
When simulating minimum-bias events using the MEMinBias matrix
element, the correct unitarized cross section can now be reported
via the standard facilities; it is no longer necessary to extract
it from the .log file of the run. The corresponding functionality
is enabled by inserting a MPIXSecReweighter object as a
post-subprocess handler:
create Herwig::MPIXSecReweighter MPIXSecReweighter
insert LHCHandler:PostSubProcessHandlers 0 MPIXSecReweighter
** Dependency on 'boost'
Herwig++ now requires the boost headers to build; if not detected
in standard locations, specify with the --with-boost configure
option.
** Tests directory
The Tests directory now contains input cards for almost all Rivet
analyses. A full comparison run can be initiated with 'make tests'.
** Minor changes
*** Default LHC energy now 8 TeV
All LHC-based defaults have now been updated to use 8 TeV as the
center-of-mass energy.
*** Herwig::ExtraParticleID -> ThePEG::ParticleID
The namespace for additional particles has been unified into
ThePEG::ParticleID
*** MEee2VectorMeson
The e+e- -> vector meson matrix element has moved from Contrib into
HwMELepton.so
*** SUSY numerics fixes
Better handling of rare numerical instabilities.
*** YODA output for Rivet
The built-in histogramming handler can now output data in the YODA
format used by Rivet.
*** Consistency checks in SLHA file reader
Better warnings for inconsistent SusyLHA files
*** better colour flow checking for development
** Bug fixes
*** Extremely offshell W from top decay
Numerical improvements for very off-shell W bosons coming from top
decays.
*** Ticket #367: problems in using SUSY + LHE
Susy events from Les Houches event files are now handled better.
*** Infinite loop in remnant decayer
The remnant decayer will now abort after 100 tries.
*** Fix to HiggsVBF LO diagrams
The diagram structure of HiggsVBF LO matrix elements has been fixed.
*** LEP thrust fix
The calculation of the transverse momentum of a branching from the
evolution variable in final-state radiation can now be
changed. While formally a sub-leading choice this enables a better
description of the thrust distribution in e+e- collisions at
small values of the thrust. Currently the default behaviour, where
the cut-off masses are used in the calculation, remains the same
as previous versions.
* Herwig++-2.5.2 release: 2011-11-01 (tagged at SVN r6928)
** Optional new jet vetoing model
The jet vetoing model by Schofield and Seymour (arXiv:1103.4811) is
available via Evolver:ColourEvolutionMethod,
PartnerFinder:PartnerMethod and SplittingFunction:SplittingColourMethod.
The default behaviour is unchanged.
** MPI tune
Version 3 of the MPI tunes is now the default. Please note that the
pT parameter is energy-dependent and needs to be modified when an
LHC run is not at 7 TeV.
The latest tunes are always available at
http://projects.hepforge.org/herwig/trac/wiki/MB_UE_tunes
** MPI PDFs
MPI PDFs can now be controlled independently.
** Initialization time speedup
A new BSMModel base class was introduced between StandardModel and
the BSM model classes. Together with a restructured decay mode
initialization, this offers significantly faster startup times for
BSM runs. ThreeBodyDecays can now always be switched on without a
large speed penalty.
** Decay mode file
Decay mode files in the SLHA format can now be read separately in
any BSM model with 'set Model:DecayFileName filename'
** Powheg DIS
Charged- and neutral-current DIS processes implementing the POWHEG
method are now available.
** Diffraction models
Xi cut implemented in PomeronFlux
** Ticket #352: Colour reconnection fixed in DIS
** Ticket #353: Improved numerical stability in chargino decays
** Ticket #358: Infinite loop in top events with pT cut in shower
** Ticket #361: Problem with duplicate 2-body modes in BSM
** Tickets #362 / #363: Crashes with baryon number violating models
Particle decays in SUSY models with RPV now work correctly in the
colour 8 -> 3,3,3 case. Colour reshuffling now works for RPV
clusters.
** Improved Fastjet detection
The configure step uses fastjet-config to make sure all header file
paths are seen.
** Darwin 11 / OS X Lion
A configure bug was fixed which prevented 'make check' from
succeeding on OS X Lion.
** Vertex classes
The specification of QED / QCD orders has been moved to the vertex
constructors, to allow ThePEG consistency checks. WWHH vertices in
MSSM and NMSSM were fixed. Some Leptoquark and UED vertices fixed.
** Hadronization
Cleanup of obsolete code.
* Herwig++-2.5.1 release: 2011-06-24 (tagged at SVN r6609)
** Example input files at 7 TeV
All our example input files for LHC now have their beam energy set
to 7 TeV instead of 14 TeV.
** Colour reconnection on by default
The colour reconnection tunes are now the default setup. Version 2
of the tunes replaces the *-1 tunes, which had a problem with LEP
event shapes.
** Run name tags
Aded possibility to add a tag to the run name when running with the
'-t' option. One run file can thus be run with different seeds and
results stored in different output files.
** Floating point exceptions
The new command line option -D enables floating point error checking.
** General improvements to WeakCurrent decays
** Remnant decayer
Hardwired gluon mass was removed.
** WeakCurrentDecayConstructor
Instead of specifying separate Particle1...Particle5 vectors for
the decay modes, the new interface DecayModes can be filled with
decay tags in the standard syntax.
** BSM: improvements to handling of vertex and model initialisation
** Powheg Higgs
Option to use pT or mT as the scale in alphaS and for the
factorization scale in the PDFs
** Ticket #337: Tau polarization wrong in charged Higgs decay
** Ticket #339: Colour flows in GeneralThreeBody Decayers for 3bar -> 8 3bar 1
** Ticket #340: Crash for resonant zero-width particles
** Ticket #341: Varying scale for BSM processes
The scale used is now ResonantProcessConstructor:ScaleFactor or
TwoToTwoProcessConstructor:ScaleFactor multiplied by sHat.
** Ticket #346: Chargino decays
Chargino decayers now automatically switch between the mesonic
decays for mass differences less than 2 GeV and the normal partonic
decays above 2 GeV.
** Ticket #349: Stop by default on input file errors
The '--exitonerror' flag is now the default behaviour for the
Herwig++ binary. To switch back to the old behaviour,
'--noexitonerror' is required.
** Ticket #351: Four-body stop decays
** Tested with gcc-4.6
* Herwig++-2.5.0 release: 2011-02-08 (tagged at SVN r6274)
** Uses ThePEG-1.7.0
Herwig++ 2.5.0 requires ThePEG 1.7.0 to benefit from various
improvements, particularly: handling of diffractive processes;
respecting LD_LIBRARY_PATH when loading dynamic libraries,
including LHAPDF; improvements to repository commands for decay
modes. See ThePEG's NEWS file for more details.
** POWHEG improvements
*** New POWHEG processes
Simulation at NLO accuracy using the POWHEG method is now
available for hadronic diboson production (pp to WW,WZ,ZZ), Higgs
decays to heavy quarks, and e+e- to two jets or ttbar, including
full mass dependence.
*** Input file changes
The input files for setting up POWHEG process simulation have been
simplified. See the example files LHC-Powheg.in and TVT-Powheg.in
for the improved command list.
*** Structural changes
The POWHEG backend in the shower code has been restructured to
make future additions easier: PowhegEvolver has merged with
Evolver; both the matrix element corrections and real corrections
in the POWHEG scheme are implemented directly in the ME or Decayer
classes.
** New processes at leading order
*** Photon initiated processes
We have added a matrix element for dijet production in gamma
hadron collisions.
*** Bottom and charm in heavy quark ME
The option of bottom and charm quarks is now supported for heavy
quark production in MEHeavyQuark.
** Colour reconnection
The cluster hadronization model has been extended by an option to
reconnect coloured constituents between clusters with a given
probability. This new model is different from the colour
reconnection model used in FORTRAN HERWIG, and improves the
description of minimum bias and underlying event data.
** Diffractive Processes
Both single and double diffractive processes are now supported in
Herwig++. The Pomeron PDF is implemented using a fit to HERA data,
and a pion PDF can be used to model reggeon flux.
** BSM physics
*** New models
We have added new BSM models, particularly ADD-type extra
dimension models and the next-to-minimal supersymmetric standard
model (NMSSM). Effects of leptoquarks can as well be simulated.
*** Vertex additions
We have added flavour changing stop interactions (stop -
neutralino - charm) and gravitino interactions with particular
emphasis on numerical stability for very light gravitinos.
Tri-linear Higgs and Higgs-Higgs/Vector-Vector four-vertices are
available as well.
*** Input file changes
The SUSY model can now also extract the SLHA information from the
header of a Les Houches event file: replace the SLHA file name
in the example input files with the LH file name.
*** Structure
The backend structure of the HardProcessConstructor has changed,
to allow easier inclusion of new process constructors. Some 2->3
BSM scattering processes involving neutral higgs bosons are now
included. The spin handling has been improved in the background.
** Shower splitting code reorganized
The selection of spin structures has been decoupled from the choice
of colour structure. This gives more flexibility in implementing
new splittings. Selected splittings can be disabled in the input
files.
** B mixing
B mixing, and indirect CP violation in the B meson system are
included now.
** Looptools
The Looptools directory has been updated to reflect T.Hahn's
Looptools 2.6.
** Contrib changes
The ROOT interface has been removed as deprecated. The MCPWNLO code
has temporarily been removed from the Contrib directory as a major
review of this code is required. Additionally, there are various
fixes to all other codes shipped in Contrib.
** DIS improvements
The momentum reshuffling in DIS events has been improved.
** mu and nu beams
mu, nu_e and nu_mu and their antiparticles are now available as
beam particles. They are all supported in the DIS matrix
elements. mu+ mu- collisions are supported in the general
matrix element code for BSM models, but not yet in the hard-coded
matrix elements for lepton-lepton scattering.
** Structural changes
*** Inline code
Inline code has been merged into the header files, .icc files were
removed.
*** Silent build
By default, Herwig++ now builds with silent build rules. To get
the old behaviour, run 'make V=1'.
*** Debug level
The debug level on the command line will now always have priority.
*** Event counter
The event counter has been simplified.
*** Interpolator persistency
Interpolators can now be written persistently.
** Ticket #307: Momentum violation check in BasicConsistency
Added parameters AbsoluteMomentumTolerance and
RelativeMomentumTolerance
** Example POWHEG input files
The example input files for Powheg processes now set the NLO
alpha_S correctly, and are run as part of 'make check'.
** Truncated shower
A problem which lead to the truncated shower not being applied in
some cases has been fixed.
** Fixes to numerical problems
Minor problems with values close to zero were fixed in several
locations.
** Remove duplicated calculation of event shapes
An accidental duplication in the calculation of event shapes was
removed, they are now only calculated once per event. Several other
minor issues in the event shape calculations have also been fixed.
** MRST PDFs fixed
An initialization problem in the internal MRST PDFs was fixed.
** Vertex scale choice
The scale in the Vertex classes can now be zero where
possible.
** Treatment of -N flag
The Herwig++ main program now correctly treats the -N flag
as optional.
** Numerical stability improved
The numerical stability in the 'RunningMass' and
'QTildeReconstructor' classes has been improved. The
stability of the boosts in the SOPTHY code for the
simulation of QED radiation has been improved.
The accuracy of boosts in the z-direction has been improved to
fix problems with extremely high p_T partons.
** Bugfix in initial state splittings
A bug in the implementation of the PDF weight in initial-state
qbar -> qbar g splittings has been fixed.
** Bugfix in chargino neutralino vertices
A bug in the 'chi+- chi0 W-+' and charged
Higgs-sfermions vertices has been fixed.
** Remove uninitialized variables written to repository
A number of uninitialised variables which were written to the
repository have been initialised to zero to avoid problems on some
systems.
** Fix to QED radiation in hadronic collisions
The longitudinal boost of the centre-of-mass frame in hadronic
collisions is correctly accounted for now in the generation of QED
radiation.
** Fix to numerical problems in two-body decays
Numerical problems have been fixed, which appeared in the rare case
that the three-momenta of the decay products in two-body decays are
zero in the rest frame of the decay particle.
** A problem with forced splittings in the Remnant was fixed.
** ME correction for W+- decays applied properly
The matrix element correction for QCD radiation in W+- decays
which was not being applied is now correctly used.
** Top quark decays from SLHA file
The presence of top quark decay modes in SLHA files is now handled
correctly.
** Exceptional shower reconstruction kinematics
Additional protection against problems due to the shower
reconstruction leading to partons with x>1 has been added.
** Ordering of particles in BSM processes
Changes have been made to allow arbitrary ordering of the outgoing
particles in BSM processes.
** Bugfixes in tau decays
Two bugs involving tau decays have been fixed. The wrong masses
were used in the 'KPiCurrent' class for the scalar form factors
and a mistake in the selection of decay products lead to
tau- --> pi0 K- being generated instead of tau- --> eta K-.
** Avoid crashes in baryon number violating processes.
To avoid crashes, better protection has been introduced for the
case where diquarks cannot be formed from the quarks in a
baryon-number violating process. In addition, the parents of the
baryon-number violating clusters have been changed to avoid
problems with the conversion of the events to HepMC.
** QED radiation in W- decays
A bug in the 'QEDRadiationHandler' class which resulted
in no QED radiation being generated in W- decays has been fixed.
** A number of minor fixes to the SUSY models have been made.
** Partial width calculations in BSM models
A fix for the direction of the incoming particle in the calculation
of two-body partial widths in BSM models has been made.
** LoopTools improvements
The LoopTools cache is now cleared more frequently to
reduce the amount of memory used by the particle.
** Negative gluino masses are now correctly handled.
** A problem with mixing matrices which are not square has been fixed.
** Removed duplicate diagram
The 'MEee2gZ2ll' class has been fixed to only include the
photon exchange diagram once rather than twice as previously.
** Fix for duplicate particles in DecayConstructor
A problem has been fixed which occurred if the same particle was
included in the list of DecayConstructor:DecayParticles.
** Fixes for UED model vertices
A number of minor problems in the vertices for the UED model have
been fixed.
** Include missing symmetry factor
The missing identical-particle symmetry factor in
'MEPP2GammaGamma' has been included.
** Fix floating point problem in top decays
A floating point problem in the matrix element correction for top
decays has been fixed.
* Herwig++-2.4.2 release: 2009-12-11 (tagged at SVN r5022)
** Ticket #292: Tau decay numerical instability
The momentum assignment in tau decays contained numerical
instabilities which have been fixed by postponing the tau decay
until after the parton shower. A new interface setting
DecayHandler:Excluded is available to prevent decays in the shower
step. This is enabled by default for tau only.
** Ticket #290: Missing MSSM colour structure
The missing colour structure for gluino -> gluon neutralino was added.
** Ticket #294: Zero momentum in some decays
Some rare phase space points lead to zero momentum in two-body
decays. This has been fixed.
** Ticket #295: Stability of QED radiation for lepton collider processes
The numerical stability of QED radiation momenta was improved
further.
** Ticket #296: K0 oscillation vertex was wrong
The oscillation from K0 to K0_L/S now takes place at the production
vertex of K0.
** Ticket #289: Undefined variables in repository
On some system configurations, undefined variables were written to
the repository. These have been fixed.
** Fixed QED radiation for hadron processes
The longitudinal boost of the centre-of-mass frame in hadronic
collisions is correctly accounted for now.
** Numerical stability fixes
Small fixes in RunningMass and QTildeReconstructor.
** Powheg example input files
The example input files for Powheg processes now set the NLO
alpha_S correctly, and are run as part of 'make check'.
** OS X builds for Snow Leopard
Snow Leopard machines will now be recognized as a 64bit
architecture.
* Herwig++-2.4.1 release: 2009-11-19 (tagged at SVN r4932)
** Uses ThePEG-1.6.0
Herwig++ now requires ThePEG-1.6.0 to benefit from the improved
helicity code there. If you have self-written vertex classes, see
ThePEG's NEWS file for conversion instructions.
** Vertex improvements
ThePEG's new helicity code allowed major simplification of the vertex
implementations for all Standard Model and BSM physics models.
** New Transplanckian scattering model
An example configuration is in LHC-TRP.in
** BSM ModelGenerator as branching ratio calculator
The BSM ModelGenerator has a new switch to output the branching
ratios for a given SLHA file in SLHA format, which can then be used
elsewhere.
** BSM debugging: HardProcessConstructor
New interface 'Excluded' to exclude certain particles from
intermediate lines.
** Chargino-Neutralino-W vertex fixed
** Spin correlations
are now switched on by default for all perturbative decays.
** Ticket #276: Scale choice in BSM models' HardProcessConstructor
New interface 'ScaleChoice' to choose process scale between
- sHat (default for colour neutral intermediates) and
- transverse mass (default for all other processes).
** Ticket #287: Powheg process scale choice
The default choice is now the mass of the colour-singlet system.
** Ticket #278: QED radiation for BSM
Soft QED radiation is now enabled in BSM decays and all
perturbative decays by default.
** Ticket #279: Full 1-loop QED radiation for Z decays
Soft QED radiation in Z decays is now fully 1-loop by default.
** Ticket #280: Redirect all files to stdout
This is now implemented globally. The files previously ending in
-UE.out and -BSMinfo.out are now appended to the log file. They now
also obey the EventGenerator:UseStdout flag.
** Ticket #270: LaTeX output updated
After each run, a LaTeX file is produced that contains the full
list of citations. Please include the relevant ones in publications.
** Ticket #256: Mac OS X problems
An initialization problem that affected only some configurations has
been identified and fixed.
** Tests directory added
This contains many .in files, to exercise most matrix
elements.
** Minor fixes
*** Prevent rare x>1 partons in shower reconstruction.
*** SUSY-LHA parameter EXTPAR can be used to set tan beta
*** Improved Fastjet detection at configure time
* Herwig++-2.4.0 release: 2009-09-01 (tagged at SVN r4616)
** New matrix elements
We have added a built-in implementation of several new matrix elements:
PP --> WW / WZ / ZZ
PP --> W gamma / Z gamma
PP --> VBF Higgs
PP --> Higgs tt / Higgs bb
e+e- --> WW / ZZ
gamma gamma --> ff / WW
** Base code improvements
*** Ticket #257: Remnant handling
A problem with forced splittings in the Remnant was fixed.
*** Ticket #264: Soft matrix element correction
A problem with emissions form antiquarks was fixed.
** PDF sets
*** New default set
MRST LO** is the new default PDF set. LO* is also available built-in.
*** Shower PDFs can be set separately from the hard process
Use the 'ShowerHandler:PDF' interface.
** Parameter tunes
Shower, hadronization and underlying event parameters were retuned
against LEP and Tevatron data respectively.
** BSM module improvements
*** Ticket #259: read error for some UED models
Arbitrary ordering of outgoing lines in the process description is now
possible.
*** Ticket #266: branching ratio sums
The warning threshold for branching ratios not summing to 1 has
been relaxed. It is now a user interface parameter.
*** Ticket #267: Top decay modes
Top decay modes listed in SLHA files are now handled correctly.
** QED radiation
*** Ticket #241: Soft QED radiation is now enabled by default
*** Ticket #265: Radiation off W+ and W- is now handled correctly
** Interfaces
*** Ticket #243: Fastjet
Fastjet is now the only supported jet finder code. All example
analyses have been converted to use Fastjet.
*** KtJet and CLHEP interfaces have been removed.
*** New interfaces to AcerDet and PGS available in Contrib
*** MCPWnlo distributed in Contrib
*** HepMC and Rivet interfaces moved to ThePEG
** Ticket #239: Inelastic cross-section for MinBias
This information is now available in the ...-UE.out files.
** Technical changes
*** Ticket #186
Configure now looks for ThePEG in the --prefix location first.
*** Configure information
Important configuration information is listed at the end of the
'configure' run and in the file 'config.thepeg'. Please provide
this file in any bug reports.
*** New ZERO object
The ZERO object can be used to set any dimensionful quantity to
zero. This avoids explicit constructs like 0.0*GeV.
*** Exception specifiers removed
Client code changes are needed in doinit() etc., simply remove the
exception specifier after the function name.
*** Ticket #263: Tau polarizations can be forced in TauDecayer
* Herwig++-2.3.2 release: 2009-05-08 (tagged at SVN r4249)
** SUSY enhancements
*** Ticket #245: Select inclusive / exclusive production
Using the new 'HardProcessConstructor:Processes' switch options
'SingleParticleInclusive', 'TwoParticleInclusive' or 'Exclusive'
*** Improved three-body decay generation
Several problems were fixed, incl. tickets #249 #250 #251
Thanks to J.Tattersall and K.Rolbiecki for the stress-testing!
*** Looptools fix
Release 2.3.1 had broken the Looptools initialization.
*** Improved warning message texts
** Ticket #237:
Values of q2last can now be zero where possible.
** Ticket #240:
The Herwig++ main program now correctly treats the -N flag as optional.
** Ticket #246:
Extreme pT partons fixed by improving accuracy of z boosts.
** DIS
Improved parton shower momentum reshuffling.
** Minimum Bias events
The zero-momentum interacting particle used for
bookkeeping is now labelled as a pomeron.
** User Makefile
Makefile-UserModules does not enable -pedantic anymore. User's ROOT
code will not compile otherwise.
** Build system
Small fixes in the build system.
* Herwig++-2.3.1 release: 2009-03-31 (tagged at SVN r4140)
** Initial state showers
The PDF veto was wrongly applied to qbar->qbar g splittings.
** User interaction
The Makefile-UserModules now includes the Herwig version number.
The -N flag to 'Herwig++ run' is optional now, as was always intended.
** Contrib
The contrib directory is now included in the tarball. The omission
was accidental.
** Numerical accuracy
Minor problems with values close to zero were fixed in several
locations.
** LEP event shapes
An accidental duplication was removed, they are now only calculated
once per event.
** MRST PDF code
Initialization problem fixed.
** Mac OS X
The configure script was improved to detect libraries better.
** Libtool
Updated to version 2.2.6
* Herwig++-2.3.0 release: 2008-12-02 (tagged at SVN r3939)
** Major release, with many new features and bug fixes
** Extension to lepton-hadron collisions
** Inclusion of several processes accurate at next-to-leading order
in the POsitive Weight Hardest Emission Generator (POWHEG) scheme
** Inclusion of three-body decays and finite-width effects
in BSM processes
** New procedure for reconstructing kinematics of the parton shower
based on the colour structure of the hard scattering process
** New model for baryon decays including excited baryon multiplets
** Addition of a soft component to the multiple scattering model
of the underlying event and the option to choose more than one hard
scattering explicitly
** New matrix elements for DIS and e+e- processes
** New /Contrib directory added
containing external modules that will hopefully be of use to some
users but are not expected to be needed by most users and are not
supported at the same level as the main Herwig++ code
** Minor changes to improve the physics simulation:
*** IncomingPhotonEvolver added
to allow the simulation of partonic processes with incoming photons
in hadron collisions
*** KTRapidityCut added
to allow cuts on the p_T and rapidity, rather than just the p_T and
pseudorapidity used in SimpleKTCut. This is now used by default for
cuts on massive particles such as the $W^\pm$, $Z^0$ and Higgs
bosons and the top quark
*** Several changes to the decayers of B mesons
both to resolve problems with the modelling of partonic decays and
improve agreement with $\Upsilon(4s)$ data
*** Changes to allow values other than transverse mass of final-state particles as maximum transverse momentum for radiation in parton shower
either SCALUP for Les Houches events or the scale of the hard
process for internally generated hard processes
*** Changed defaults for intrinsic transverse momentum in hadron collisions
to 1.9GeV, 2.1GeV and 2.2GeV for the Tevatron and LHC at 10 TeV and
14 TeV, respectively
*** Pdfinfo object is now created in the HepMC interface
However in order to support all versions of HepMC containing this
feature the PDF set is not specified as not all versions contain
this information
*** New option of only decaying particles with lifetimes below user specified value
*** New options for the cut-off in the shower
and some obsolete parameters removed
*** Added option of switching off certain decay modes in BSM models
*** Added a Matcher for Higgs boson
to allow cuts to be placed on it
*** Diffractive particles deleted from default input files
they were not previously used
** Technical changes:
*** Some AnalysisHandler classes comparing to LEP data have been renamed
e.g. MultiplicityCount becomes LEPMultiplicityCount to avoid
confusion with those supplied in /Contrib for observables at the
Upsilon(4s) resonance
*** Reorganisation to remove the majority of the .icc files
by moving inlined functions to headers in an effort to improve
compile time
*** Restructured the decay libraries to reduce the amount of memory allocation
and de-allocation which improves run-time performance
*** The switch to turn off LoopTools has been removed
because LoopTools is now used by several core modules. As LoopTools
does not work on 64-bit platforms with g77 this build option is not
supported
*** Removed support for obsolete version of HepMC supplied with CLHEP
and improved the support for different units options with HepMC
*** EvtGen interface has been removed until it is more stable
*** Support for ROOT has been removed
it was not previously used
*** CKKW infrastructure has been removed from the release
until a concrete implementation is available
*** Default optimisation has been increased from -O2 to -O3
*** Handling of the fortran compiler has been improved
mainly due to improvements in the autotools
*** Use of FixedAllocator for Particle objects in ThePEG has been removed
as it had no performance benefits
** Bugs fixed:
*** Problems with the mother/daughter relations in the hard process
and diagram selection in W+- and Z0 production in association with a
hard jet
*** In general matrix element code for fermion-vector to fermion-scalar
where the outgoing fermion is coloured and the scalar neutral
*** In the selection of diagrams in some associated squark gaugino processes
*** h0->mu+mu- was being generated when h0->tau+tau-
*** Normalisation in the Histogram class for non unit-weight events
*** Protection against negative PDF values has been improved
these can occur when using NLO PDF sets
*** Lifetime for BSM particles is now automatically calculated
at the same time as the width
*** Hadrons containing a top quark now treated like hadrons containing BSM particles
in order to support this possibility
*** Several ambiguous uses of unsigned int
*** Several variables that may have been used undefined
*** Several memory leaks at initialisation
*** The configuration now aborts if no fortran compiler is found
as this is required to compile Looptools
*** Several minor floating point errors that did not affect results
* Herwig++-2.2.1 release: 2008-07-09 (tagged at SVN r3434)
** Ticket #181: BSM shower with a decay close to threshold
Now fixed.
** Ticket #191: Split SUSY crash
Improved error message.
** Ticket #192: using SCALUP as the pT veto in the shower
Now implemented.
** Ticket #194: production processes of ~chi_1(2)-
Fixed bug in the diagram creation.
** Removed unused particles
DiffractiveParticles.in was removed, they were never produced.
** Hadronization
Top quark clusters now possible, handled as 'exotic' clusters.
** Improved handling of decay modes
See ThePEG-1.3.0. 'defaultparticle' command is now obsolete.
** Multi-Parton interactions
Increased phase space sampling to have less than 1% uncertainty on
average multiplicity.
** New libtool version
gfortran is now used as default if it is available. Set FC=g77 to
override this.
** Fixed several memory leaks
** Memory allocation
Now using plain 'new' and 'delete'.
* Herwig++-2.2.0 release: 2008-04-18 (tagged at SVN r3195)
** Major release: now as stand-alone library
Herwig++ is now a stand-alone dlopen() plugin to ThePEG.
No compile-time linking to Herwig code is required. The Herwig++
binary is a simple executable steering ThePEG, which can
be replaced by other frontends (such as setupThePEG / runThePEG).
** New matrix elements
p p -> W + jet / Z + jet / W + higgs / Z + higgs
e+ e- -> Z + higgs
** Looptools
Updated to version 2.2.
** Ticket #141: segfault from using 'run' command
Fixed by using default allocators in Herwig++, and the
Repository::cleanup() method in ThePEG 1.2.0.
** Ticket #157: broken gsl library path on some 64bit systems
Paths with lib64 are correctly identified now.
** Ticket #159: p_t spectrum of ttbar pair
Fixed identical particle multiplier in Sudakov form factor.
** Ticket #161: glibc segfault
Rare segfault in MPI handler fixed.
** Ticket #165: rare infinite loop in four-body decay
All 4-body decays without dedicated decayers now use the Mambo algorithm.
A loop guard has been introduced to 3-body decays to avoid infinite retries.
** Ticket #166: rare infinite loop in top ME correction
These very rare events (O(1) in 10^7) close to mass threshold
now are discarded.
** Higgs width fixes
** SatPDF
Optionally, the PDF extrapolation behaviour outside a given range
can now be specified.
** gcc 4.3
Herwig++-2.2 compiles cleanly with the new gcc 4.3 series.
* Herwig++-2.1.4 release: 2008-03-03 (tagged at SVN r3024)
** Ticket #152: Vertex positions
All vertex positions of unphysical particles are set to zero until
a fix for the previous nonsensical values can be implemented.
* Herwig++-2.1.3 release: 2008-02-25 (tagged at SVN r2957)
** Ticket #129: Baryon decays
Fix for baryon decay modes.
** Ticket #131: HepMC
Check if IO_GenEvent exists
** Ticket #134: Hadronization
Smearing of hadron directions in cluster decay fixed.
** Ticket #137: HepMC
HepMC conversion allows specification of energy and length units to
be used.
** Ticket #139: Neutral kaons
Ratio K_L / K_S corrected.
** Ticket #140 / #141: Crash on shutdown
Event generation from the 'read' stage or an interface now shuts
down cleanly. Fixes a crash bug introduced in 2.1.1 which affected
external APIs to ThePEG / Herwig.
** Ticket #146: BSM models can be disabled
To save build time, some or all of the BSM models can be disabled
using the '--enable-models' configure switch.
** Reorganised .model files
The .model files now include the model-specific particles, too.
** Re-tune
Re-tuned hadronization parameters to LEP data.
** Other fixes in
QSPAC implementation in Shower; Multi-parton interaction tuning;
MRST initialization
* Herwig++-2.1.2 release: 2008-01-05 (tagged at SVN r2694)
** Ticket #127
Thanks to a patch submitted by Fred Stober, HepMCFile now can
output event files in all supported formats.
** Ticket #128
Fixed incorrect value of pi in histogram limits.
** Other fixes in
CKKW Qtilde clusterers, BSM width cut, SUSY mixing matrices.
* Herwig++-2.1.1 release: 2007-12-08 (tagged at SVN r2589)
** Bug #123
Fixed a bug with particle lifetimes which resulted in nan for some
vertex positions.
** Secondary scatters
Fixed bug which gave intrinsic pT to secondary scatters.
** gcc abs bug detection
configure now checks for and works around
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34130
** CKKW reweighting
Fixed wrong check for top quarks.
** MPIHandler
Fixed call order ambiguity.
* Herwig++-2.1.0 release: 2007-11-20 (tagged at SVN r2542)
** Major new release
Herwig++-2.1 includes significant improvements, including
multi-parton interactions, BSM physics and a new hadronic decay
model, tuned to LEP data.
For an overview of the changes, please see the release note
arXiv:0711.3137
* Herwig++-2.0.3 release: 2007-08-21 (tagged at SVN r2101)
** Bug #90
nan in top decay ME corrections fixed.
** unlisted
Colour flow fix in LightClusterDecayer
** unlisted
Updated version of MultiplicityCount analysis handler.
* Herwig++-2.0.2 release: 2007-07-06 (tagged at SVN r1716)
** Bug #80
Separation of HepMC from CLHEP is handled properly now.
** Bug #83
Workaround for OS X header problem
** unlisted
Veto on very hard emissions from Shower.
** unlisted
Detailed documentation in .in files
* Herwig++-2.0.1 release: 2006-12-05 (tagged at SVN r1195)
** Bug #54
ClusterFissioner vertex calculation fixed.
** Bug #57
Crash when showering W+jet events supplied by Les Houches interface.
** Bug #59
Fix for #57 applied to LHC events.
** Bug #60
Segfault when PDF is set to NoPDF.
** Bug #61
Missing weight factor for I=0 mesons
** Bug #62
Spinor vertex calculations broken when spinor rep is not default rep.
** Bug #63
Top decay never produces tau.
** Bug #69
TTbar and HiggsJet analysis handlers fixed.
** unlisted
Reorganization of Hadronization module gives 30% speedup.
Thanks to Vincenzo Innocente at CMS for his profiling work!
** unlisted
cleaner automake files in include/ and src/
** unlisted
Hw64 hadron selection algorithm 'abortnow' fixed.
** unlisted
Top/LeptonDalitzAnalysis removed (only worked with modified code).
** unlisted
removed f'_0 from particle list, decays were not handled
* Herwig++-2.0.0 release: 2006-09-28 (tagged at SVN r1066)
** Full simulation of hadron collisions
diff --git a/PDT/GenericWidthGenerator.cc b/PDT/GenericWidthGenerator.cc
--- a/PDT/GenericWidthGenerator.cc
+++ b/PDT/GenericWidthGenerator.cc
@@ -1,889 +1,889 @@
// -*- C++ -*-
//
// GenericWidthGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GenericWidthGenerator class.
//
#include "GenericWidthGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "TwoBodyAllOnCalculator.h"
#include "OneOffShellCalculator.h"
#include "TwoOffShellCalculator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <ctime>
using namespace Herwig;
DescribeClass<GenericWidthGenerator,WidthGenerator>
describeHerwigGenericWidthGenerator("Herwig::GenericWidthGenerator","");
HERWIG_INTERPOLATOR_CLASSDESC(GenericWidthGenerator,Energy,Energy)
void GenericWidthGenerator::persistentOutput(PersistentOStream & os) const {
os << particle_ << ounit(mass_,GeV) << prefactor_ << MEtype_ << MEcode_
<< ounit(MEmass1_,GeV) << ounit(MEmass2_,GeV) << MEcoupling_ << modeOn_
<< ounit(interMasses_,GeV) << ounit(interWidths_,GeV)
<< noOfEntries_ << initialize_ << output_ << BRnorm_ << twoBodyOnly_
<< npoints_ << decayModes_ << decayTags_ << ounit(minMass_,GeV)
<< BRminimum_ << intOrder_ << interpolators_;
}
void GenericWidthGenerator::persistentInput(PersistentIStream & is, int) {
is >> particle_ >> iunit(mass_,GeV) >> prefactor_ >> MEtype_ >> MEcode_
>> iunit(MEmass1_,GeV) >> iunit(MEmass2_,GeV) >> MEcoupling_ >>modeOn_
>> iunit(interMasses_,GeV) >> iunit(interWidths_,GeV)
>> noOfEntries_ >> initialize_ >> output_ >> BRnorm_ >> twoBodyOnly_
>> npoints_ >> decayModes_ >> decayTags_ >> iunit(minMass_,GeV)
>> BRminimum_ >> intOrder_ >> interpolators_;
}
void GenericWidthGenerator::setParticle(string p) {
if ( (particle_ = Repository::GetPtr<tPDPtr>(p)) ) return;
particle_ = Repository::findParticle(StringUtils::basename(p));
if ( ! particle_ )
Throw<InterfaceException>()
<< "Could not set Particle interface "
<< "for the object \"" << name()
<< "\". Particle \"" << StringUtils::basename(p) << "\" not found."
<< Exception::runerror;
}
string GenericWidthGenerator::getParticle() const {
return particle_ ? particle_->fullName() : "";
}
void GenericWidthGenerator::Init() {
static ClassDocumentation<GenericWidthGenerator> documentation
("The GenericWidthGenerator class is the base class for running widths");
static Parameter<GenericWidthGenerator,string> interfaceParticle
("Particle",
"The particle for which this is the width generator",
0, "", true, false,
&GenericWidthGenerator::setParticle,
&GenericWidthGenerator::getParticle);
static Switch<GenericWidthGenerator,bool> interfaceInitialize
("Initialize",
"Initialize the width using the particle data object",
&GenericWidthGenerator::initialize_, false, false, false);
static SwitchOption interfaceInitializeInitialization
(interfaceInitialize,
"Yes",
"Do the initialization",
true);
static SwitchOption interfaceInitializeNoInitialization
(interfaceInitialize,
"No",
"Don't do the initalization",
false);
static Switch<GenericWidthGenerator,bool> interfaceOutput
("Output",
"Output the setup",
&GenericWidthGenerator::output_, false, false, false);
static SwitchOption interfaceOutputYes
(interfaceOutput,
"Yes",
"Output the data",
true);
static SwitchOption interfaceOutputNo
(interfaceOutput,
"No",
"Don't output the data",
false);
static ParVector<GenericWidthGenerator,int> interfacemetype
("MEtype",
"The type of matrix element either 2-body from this class or higher from"
" class inheriting from this",
&GenericWidthGenerator::MEtype_,
0, 0, 0, 0, 3, false, false, true);
static ParVector<GenericWidthGenerator,int> interfacemecode
("MEcode",
"The code of matrix element either 2-body from this class or higher from"
" class inheriting from this",
&GenericWidthGenerator::MEcode_,
0, 0, 0, -1, 200, false, false, true);
static ParVector<GenericWidthGenerator,Energy> interfaceMinimumMasses
("MinimumMasses",
"The minimum mass of the decay products",
&GenericWidthGenerator::minMass_,
GeV, 0, ZERO, ZERO, 1.E12*GeV, false, false, true);
static ParVector<GenericWidthGenerator,double> interfaceMEcoupling
("MEcoupling",
"The coupling for a given ME",
&GenericWidthGenerator::MEcoupling_,
0, 0, 0, 0, 1.E12, false, false, true);
static ParVector<GenericWidthGenerator,bool> interfaceModeOn
("ModeOn",
"Is this mode included in the total width calculation",
&GenericWidthGenerator::modeOn_,
0, 0, 0, 0, 1, false, false, true);
static ParVector<GenericWidthGenerator,Energy> interfaceMEmass1
("MEmass1",
"The mass for first particle in a two body mode",
&GenericWidthGenerator::MEmass1_,
GeV, 0, ZERO, ZERO, 1.E12*GeV, false, false, true);
static ParVector<GenericWidthGenerator,Energy> interfaceMEmass2
("MEmass2",
"The mass for second particle in a two body mode",
&GenericWidthGenerator::MEmass2_,
GeV, 0, ZERO, ZERO, 1.E12*GeV, false, false, true);
static ParVector<GenericWidthGenerator,Energy> interfaceInterpolationMasses
("InterpolationMasses",
"The masses for interpolation table",
&GenericWidthGenerator::interMasses_,
GeV, 0, ZERO, ZERO, 1E12*GeV, false, false, true);
static ParVector<GenericWidthGenerator,Energy> interfaceInterpolationWidths
("InterpolationWidths",
"The widths for interpolation table",
&GenericWidthGenerator::interWidths_,
GeV, 0, ZERO, ZERO, 1E12*GeV, false, false, true);
static ParVector<GenericWidthGenerator,int> interfacenoofenteries
("NumberofEntries",
"The number of entries in the table after this mode",
&GenericWidthGenerator::noOfEntries_,
0, 0, 0, 0, 100000000, false, false, true);
static Switch<GenericWidthGenerator,bool> interfaceBRNormalize
("BRNormalize",
"Normalize the partial widths so that they have the value BR*Total Width"
" for an on-shell particle",
&GenericWidthGenerator::BRnorm_, false, false, false);
static SwitchOption interfaceBRNormalizeNormalize
(interfaceBRNormalize,
"Yes",
"Perform the normalization",
true);
static SwitchOption interfaceBRNormalizeNoNormalisation
(interfaceBRNormalize,
"No",
"Do not perform the normalization",
false);
static Parameter<GenericWidthGenerator,double> interfaceBRMinimum
("BRMinimum",
"Minimum branching ratio for inclusion in the running width calculation.",
&GenericWidthGenerator::BRminimum_, 0.01, 0.0, 1.0,
false, false, true);
static Parameter<GenericWidthGenerator,double> interfacePrefactor
("Prefactor",
"The prefactor to get the correct on-shell width",
&GenericWidthGenerator::prefactor_, 1.0, 0., 1000.,
false, false, false);
static Parameter<GenericWidthGenerator,int> interfacePoints
("Points",
"Number of points to use for interpolation tables when needed",
&GenericWidthGenerator::npoints_, 50, 5, 1000,
false, false, true);
static ParVector<GenericWidthGenerator,string> interfaceDecayModes
("DecayModes",
"The tags for the decay modes used in the width generator",
&GenericWidthGenerator::decayTags_, -1, "", "", "",
false, false, Interface::nolimits);
static Parameter<GenericWidthGenerator,unsigned int> interfaceInterpolationOrder
("InterpolationOrder",
"The interpolation order for the tables",
&GenericWidthGenerator::intOrder_, 1, 1, 5,
false, false, Interface::limited);
static Switch<GenericWidthGenerator,bool> interfaceTwoBodyOnly
("TwoBodyOnly",
"Only Use two-body modes for the calculation of the running "
"width, higher multiplicity modes fixed partial width",
&GenericWidthGenerator::twoBodyOnly_, false, false, false);
static SwitchOption interfaceTwoBodyOnlyYes
(interfaceTwoBodyOnly,
"Yes",
"Only include two-body modes",
true);
static SwitchOption interfaceTwoBodyOnlyNo
(interfaceTwoBodyOnly,
"No",
"Include all modes",
false);
}
Energy GenericWidthGenerator::width(const ParticleData &, Energy m) const {
Energy gamma= ZERO;
for(unsigned int ix =0;ix<MEcoupling_.size();++ix) {
if(modeOn_[ix]) gamma += partialWidth(ix,m);
}
return gamma*prefactor_;
}
void GenericWidthGenerator::doinit() {
WidthGenerator::doinit();
if(particle()->widthGenerator()!=this) return;
// make sure the particle data object was initialized
particle_->init();
tDecayIntegratorPtr decayer;
// mass of the decaying particle
mass_ = particle_->mass();
if(initialize_) {
// the initial prefactor
prefactor_=1.;
// resize all the storage vectors
MEtype_.clear();
MEcode_.clear();
MEmass1_.clear();
MEmass2_.clear();
MEcoupling_.clear();
modeOn_.clear();
minMass_.clear();
interMasses_.clear();
interWidths_.clear();
noOfEntries_.clear();
decayTags_.clear();
// integrators that we may need
WidthCalculatorBasePtr widthptr;
// get the list of decay modes as a decay selector
DecayMap modes=particle_->decaySelector();
- if ( Debug::level > 0 )
+ if ( Debug::level > 1 )
Repository::cout() << "Width generator for "
<< particle_->PDGName() << endl;
DecayMap::const_iterator start=modes.begin();
DecayMap::const_iterator end=modes.end();
tPDPtr part1,part2;
tGenericMassGeneratorPtr massgen1,massgen2;
// loop over the decay modes to get the partial widths
for(;start!=end;++start) {
// the decay mode
tcDMPtr mode=(*start).second;
clock_t time = std::clock();
if ( Debug::level > 1 ) {
Repository::cout() << "Partial width "
<< left
<< std::setw(40)
<< mode->tag() << flush;
}
decayModes_.push_back(const_ptr_cast<DMPtr>(mode));
decayTags_.push_back(decayModes_.back()->tag());
ParticleMSet::const_iterator pit(mode->products().begin());
// minimum mass for the decaymode
Energy minmass = ZERO;
for(;pit!=mode->products().end();++pit) {
(**pit).init();
minmass+=(**pit).massMin();
}
minMass_.push_back(minmass);
pit=mode->products().begin();
// its decayer
decayer=dynamic_ptr_cast<tDecayIntegratorPtr>(mode->decayer());
if(decayer) decayer->init();
// if there's no decayer then set the partial width to the br times the
// on-shell value
if(!decayer) {
MEtype_.push_back(0);
MEcode_.push_back(0);
MEcoupling_.push_back(mode->brat());
MEmass1_.push_back(ZERO);
MEmass2_.push_back(ZERO);
noOfEntries_.push_back(interMasses_.size());
modeOn_.push_back(mode->brat()>BRminimum_);
setupMode(mode,decayer,MEtype_.size()-1);
}
else if(mode->products().size()==2) {
// the outgoing particles
ParticleMSet::const_iterator pit = mode->products().begin();
part1=*pit;++pit;
part2=*pit;
// mass generators
if( part1->massGenerator() )
massgen1=dynamic_ptr_cast<tGenericMassGeneratorPtr>(part1->massGenerator());
else
massgen1=tGenericMassGeneratorPtr();
if( part2->massGenerator() )
massgen2=dynamic_ptr_cast<tGenericMassGeneratorPtr>(part2->massGenerator());
else
massgen2=tGenericMassGeneratorPtr();
if(massgen1) massgen1->init();
if(massgen2) massgen2->init();
double coupling(0.);
int mecode(-1);
bool order(decayer->twoBodyMEcode(*mode,mecode,coupling));
MEcode_.push_back(mecode);
MEcoupling_.push_back(coupling);
modeOn_.push_back(mode->brat()>BRminimum_);
if(order) {
MEmass1_.push_back(part1->mass());
MEmass2_.push_back(part2->mass());
}
else {
MEmass1_.push_back(part2->mass());
MEmass2_.push_back(part1->mass());
}
// perform setup in the inheriting class
setupMode(mode,decayer,MEcode_.size()-1);
// both particles on shell
if(!massgen1&&!massgen2) {
MEtype_.push_back(1);
noOfEntries_.push_back(interMasses_.size());
if(BRnorm_) {
if(mass_>MEmass1_[MEtype_.size()-1]+MEmass2_[MEtype_.size()-1]) {
Energy gamma(partial2BodyWidth(MEtype_.size()-1,mass_));
if(gamma==ZERO) {
cerr << "Partial width for " << mode->tag()
<< " is zero in GenericWidthGenerator::doinit().\n"
<< "If doing BSM physics this is probably a problem with your input "
<< "parameters.\n"
<< "Zeroing mode\n";
MEcoupling_.back() = 0.;
}
else {
double ratio(mode->brat()*mode->parent()->width()/gamma);
ratio=sqrt(ratio);
MEcoupling_.back() *=ratio;
}
}
}
}
else {
// one off-shell particle
if(!massgen1||!massgen2) {
// create the width calculator
tGenericWidthGeneratorPtr
ttthis(const_ptr_cast<tGenericWidthGeneratorPtr>(this));
WidthCalculatorBasePtr twobody
(new_ptr(TwoBodyAllOnCalculator(ttthis,MEcode_.size()-1,
MEmass1_[MEcode_.size()-1],
MEmass2_[MEcode_.size()-1])));
int ioff = ((part1->massGenerator()&&!order)||
(part2->massGenerator()&&order)) ? 2 : 1;
if(massgen1)
widthptr=new_ptr(OneOffShellCalculator(ioff,twobody,massgen1,ZERO));
else
widthptr=new_ptr(OneOffShellCalculator(ioff,twobody,massgen2,ZERO));
}
else {
int ioff = order ? 1 : 2;
int iother = order ? 2 : 1;
// create the width calculator
tGenericWidthGeneratorPtr
ttthis(const_ptr_cast<tGenericWidthGeneratorPtr>(this));
// this is the both on-shell case
WidthCalculatorBasePtr twobody
(new_ptr(TwoBodyAllOnCalculator(ttthis,MEcode_.size()-1,
MEmass1_[MEcode_.size()-1],
MEmass2_[MEcode_.size()-1])));
// this is the first off-shell
WidthCalculatorBasePtr widthptr2=
new_ptr(OneOffShellCalculator(ioff,twobody,massgen1,ZERO));
widthptr=new_ptr(TwoOffShellCalculator(iother,widthptr2,massgen2,
ZERO,massgen1->lowerLimit()));
}
// set up the interpolation table
Energy test(part1->massMin()+part2->massMin());
Energy min(max(particle_->massMin(),test)),upp(particle_->massMax());
Energy step((upp-min)/(npoints_-1));
Energy moff(min);
Energy2 moff2;
// additional points to improve the interpolation
if(min==test) {
interMasses_.push_back(moff-2.*step);interWidths_.push_back(ZERO);
interMasses_.push_back(moff- step);interWidths_.push_back(ZERO);
interMasses_.push_back(moff );interWidths_.push_back(ZERO);
double fact(exp(0.1*log(1.+step/moff)));
for(unsigned int ix=0;ix<10;++ix) {
moff*=fact;
moff2=sqr(moff);
interMasses_.push_back(moff);
interWidths_.push_back(widthptr->partialWidth(moff2));
}
moff+=step;
}
else if(test>min-2.*step) {
interMasses_.push_back(moff-2.*step);interWidths_.push_back(ZERO);
interMasses_.push_back(test );interWidths_.push_back(ZERO);
}
else {
interMasses_.push_back(moff-2.*step);
interWidths_.push_back(widthptr->partialWidth((moff-2.*step)*
(moff-2.*step)));
interMasses_.push_back(moff- step);
interWidths_.push_back(widthptr->partialWidth((moff- step)*
(moff- step)));
}
for(; moff<upp+2.5*step;moff+=step) {
moff2=moff*moff;
interMasses_.push_back(moff);
interWidths_.push_back(widthptr->partialWidth(moff2));
}
if(BRnorm_) {
double ratio(1.);
if((massgen1&&massgen2&&
mass_>massgen1->lowerLimit()+massgen2->lowerLimit())||
(massgen1&&!massgen2&&
mass_>massgen1->lowerLimit()+part2->mass())||
(massgen2&&!massgen1&&
mass_>massgen2->lowerLimit()+part1->mass())||
(!massgen1&&!massgen2&&
mass_>part1->mass()+part2->mass())) {
Energy gamma(widthptr->partialWidth(mass_*mass_));
if(gamma==ZERO) {
cerr << "Partial width for " << mode->tag()
<< " is zero in GenericWidthGenerator::doinit()"
<< " if doing BSM physics this is probably a problem with your input "
<< "parameters.\n"
<< "Zeroing mode\n";
ratio = 0.;
}
else {
ratio=mode->brat()*mode->parent()->width()/gamma;
}
}
MEcoupling_.back()=ratio;
}
else MEcoupling_.back()=1.;
MEtype_.push_back(2);
MEcode_.back()=0;
noOfEntries_.push_back(interMasses_.size());
interpolators_.resize(MEtype_.size());
// get the vectors we will need
vector<Energy>::iterator istart= interMasses_.begin();
if(MEtype_.size()>1){istart+=noOfEntries_[MEtype_.size()-2];}
vector<Energy>::iterator iend=interMasses_.end();
vector<Energy> masses(istart,iend);
istart= interWidths_.begin();
if(MEtype_.size()>1){istart+=noOfEntries_[MEtype_.size()-2];}
iend=interWidths_.end();
vector<Energy> widths(istart,iend);
interpolators_.back() = make_InterpolatorPtr(widths,masses,intOrder_);
}
}
// higher multiplicities
else {
setupMode(mode,decayer,MEcode_.size());
widthptr = twoBodyOnly_ ? WidthCalculatorBasePtr() : decayer->threeBodyMEIntegrator(*mode);
if(!widthptr) {
MEtype_.push_back(0);
MEcode_.push_back(0);
MEcoupling_.push_back(mode->brat());
MEmass1_.push_back(ZERO);
MEmass2_.push_back(ZERO);
noOfEntries_.push_back(interMasses_.size());
modeOn_.push_back(mode->brat()>BRminimum_);
}
else {
Energy step((particle_->widthUpCut()+particle_->widthLoCut())/
(npoints_-1));
Energy moff(particle_->massMin()),upp(particle_->massMax());
for( ; moff<upp+0.5*step;moff+=step) {
Energy2 moff2=sqr(moff);
Energy wtemp=widthptr->partialWidth(moff2);
interMasses_.push_back(moff);
interWidths_.push_back(wtemp);
}
double coupling(1.);
if(BRnorm_) {
Energy gamma = widthptr->partialWidth(mass_*mass_);
if(gamma==ZERO) {
cerr << "Partial width for " << mode->tag()
<< " is zero in GenericWidthGenerator::doinit()"
<< " if doing BSM physics this is probably a problem with your input "
<< "parameters.\n"
<< "Zeroing mode\n";
coupling = 0.;
}
else {
coupling = mode->brat()*mode->parent()->width()/gamma;
}
}
MEtype_.push_back(2);
MEcode_.push_back(0);
MEcoupling_.push_back(coupling);
MEmass1_.push_back(ZERO);
MEmass2_.push_back(ZERO);
modeOn_.push_back(mode->brat()>BRminimum_);
noOfEntries_.push_back(interMasses_.size());
interpolators_.resize(MEtype_.size());
// get the vectors we will need
vector<Energy>::iterator istart( interMasses_.begin()),
iend(interMasses_.end());
if(MEtype_.size()>1){istart+=noOfEntries_[MEtype_.size()-2];}
vector<Energy> masses(istart,iend);
istart= interWidths_.begin();
if(MEtype_.size()>1){istart+=noOfEntries_[MEtype_.size()-2];}
iend=interWidths_.end();
vector<Energy> widths(istart,iend);
interpolators_.back() = make_InterpolatorPtr(widths,masses,intOrder_);
}
}
if ( Debug::level > 1 ) {
double diff = double(std::clock()-time)/CLOCKS_PER_SEC;
if ( diff > 0.2 )
Repository::cout() << ' ' << diff << " s";
Repository::cout() << endl;
}
}
// now check the overall normalisation of the running width
Energy gamma = width(*particle_,mass_);
if(gamma>ZERO) prefactor_ = particle_->width()/gamma;
// output the info so it can be read back in
}
else {
// get the decay modes from the tags
if(decayTags_.size()!=0) {
decayModes_.clear();
for(unsigned int ix=0;ix<decayTags_.size();++ix) {
decayModes_.push_back(CurrentGenerator::current().findDecayMode(decayTags_[ix]));
if(!decayModes_.back())
generator()->log() << "Error in GenericWidthGenerator::doinit(). "
<< "Failed to find DecayMode for tag"
<< decayTags_[ix] << "\n";
}
}
// otherwise just use the modes from the selector
else {
DecaySet modes(particle_->decayModes());
DecaySet::const_iterator start(modes.begin()),end(modes.end());
tcDMPtr mode;
for(;start!=end;++start) {
decayModes_.push_back(const_ptr_cast<DMPtr>(*start));
}
}
// set up the interpolators
interpolators_.resize(MEtype_.size());
vector<Energy>::iterator estart(interMasses_.begin()),eend;
vector<Energy>::iterator wstart(interWidths_.begin()),wend;
vector<Energy> masses,widths;
for(unsigned int ix=0;ix<MEtype_.size();++ix) {
eend=interMasses_.begin()+noOfEntries_[ix];
wend=interWidths_.begin()+noOfEntries_[ix];
if(MEtype_[ix]==2) {
masses.assign(estart,eend);
widths.assign(wstart,wend);
interpolators_[ix]= make_InterpolatorPtr(widths,masses,intOrder_);
}
estart=eend;
wstart=wend;
}
}
// setup the partial widths in the decayers for normalization
tDecayIntegratorPtr temp;
for(unsigned int ix=0;ix<decayModes_.size();++ix) {
if(!decayModes_[ix]) continue;
decayModes_[ix]->init();
decayer=dynamic_ptr_cast<tDecayIntegratorPtr>(decayModes_[ix]->decayer());
if(!decayer) continue;
decayer->init();
if(particle_->widthGenerator() &&
particle_->widthGenerator()==this ) decayer->setPartialWidth(*decayModes_[ix],ix);
}
if ( Debug::level > 29 ) {
// code to output plots
string fname = CurrentGenerator::current().filename() +
string("-") + name() + string(".top");
ofstream output(fname.c_str());
Energy step = (particle_->massMax()-particle_->massMin())/100.;
output << "SET FONT DUPLEX\n";
output << "TITLE TOP \"Width for " << particle_->name() << "\"\n";
output << "TITLE BOTTOM \"m/GeV\"\n";
output << "TITLE LEFT \"G/GeV\"\n";
output << "CASE \"F \"\n";
output << "SET LIMITS X "
<< (particle_->massMin()-10.*step)/GeV << " "
<< particle_->massMax()/GeV << "\n";
Energy upper(ZERO);
for(Energy etest=particle_->massMin();etest<particle_->massMax();etest+=step) {
Energy gamma=width(*particle_,etest);
upper = max(gamma,upper);
output << etest/GeV << "\t" << gamma/GeV << "\n";
}
output << "SET LIMITS Y 0. " << upper/GeV << "\n";
output << "JOIN\n";
output << (particle_->massMin()-9.*step)/GeV << "\t"
<< upper*(MEcode_.size()+1)/(MEcode_.size()+2)/GeV << "\n";
output << (particle_->massMin()-7.*step)/GeV << "\t"
<< upper*(MEcode_.size()+1)/(MEcode_.size()+2)/GeV << "\n";
output << "JOIN\n";
output << "TITLE DATA "
<< (particle_->massMin()-6.*step)/GeV << "\t"
<< upper*(MEcode_.size()+1)/(MEcode_.size()+2)/GeV
<< " \"total\"\n";
for(unsigned int ix=0;ix<MEcode_.size();++ix) {
for(Energy etest=particle_->massMin();etest<particle_->massMax();etest+=step) {
output << etest/GeV << "\t" << partialWidth(ix,etest)*prefactor_/GeV << "\n";
}
switch(ix) {
case 0: output << "join red\n" ; break;
case 1: output << "join blue\n" ; break;
case 2: output << "join green\n" ; break;
case 3: output << "join yellow\n" ; break;
case 4: output << "join magenta\n"; break;
case 5: output << "join cyan\n" ; break;
case 6: output << "join dashes\n" ; break;
case 7: output << "join dotted\n" ; break;
case 8: output << "join dotdash\n"; break;
default: output << "join daashes space\n"; break;
}
output << (particle_->massMin()-9.*step)/GeV << "\t"
<< upper*(MEcode_.size()-ix)/(MEcode_.size()+2)/GeV << "\n";
output << (particle_->massMin()-7.*step)/GeV << "\t"
<< upper*(MEcode_.size()-ix)/(MEcode_.size()+2)/GeV << "\n";
switch(ix) {
case 0: output << "join red\n" ; break;
case 1: output << "join blue\n" ; break;
case 2: output << "join green\n" ; break;
case 3: output << "join yellow\n" ; break;
case 4: output << "join magenta\n"; break;
case 5: output << "join cyan\n" ; break;
case 6: output << "join dashes\n" ; break;
case 7: output << "join dotted\n" ; break;
case 8: output << "join dotdash\n"; break;
default: output << "join daashes space\n"; break;
}
output << "TITLE DATA "
<< (particle_->massMin()-6.*step)/GeV << "\t"
<< upper*(MEcode_.size()-ix)/(MEcode_.size()+2)/GeV
<< " \"" << decayTags_[ix] << "\"\n";
}
}
}
void GenericWidthGenerator::dataBaseOutput(ofstream & output, bool header) {
if(header) output << "update Width_Generators set parameters=\"";
// prefactor and general switiches
output << "newdef " << name() << ":Prefactor " << prefactor_ << "\n";
output << "newdef " << name() << ":BRNormalize " << BRnorm_ << "\n";
output << "newdef " << name() << ":BRMinimum " << BRminimum_ << "\n";
output << "newdef " << name() << ":Points " << npoints_ << "\n";
output << "newdef " << name() << ":InterpolationOrder " << intOrder_ << "\n";
// the type of the matrix element
for(unsigned int ix=0;ix<MEtype_.size();++ix) {
output << "insert " << name() << ":MEtype " << ix << " "
<< MEtype_[ix] << "\n";
}
// the code for thew two body matrix elements
for(unsigned int ix=0;ix<MEcode_.size();++ix) {
output << "insert " << name() << ":MEcode "
<< ix << " " << MEcode_[ix] << "\n";
}
// the coupling for trhe two body matrix elements
for(unsigned int ix=0;ix<MEcoupling_.size();++ix) {
output << "insert " << name() << ":MEcoupling "
<< ix << " " << MEcoupling_[ix] << "\n";
}
// use this mode for the running width
for(unsigned int ix=0;ix<modeOn_.size();++ix) {
output << "insert " << name() << ":ModeOn "
<< ix << " " << modeOn_[ix] << "\n";
}
// first outgoing mass
for(unsigned int ix=0;ix<minMass_.size();++ix) {
output << "insert " << name() << ":MinimumMasses "
<< ix << " " << minMass_[ix]/GeV << "\n";
}
// first outgoing mass
for(unsigned int ix=0;ix<MEmass1_.size();++ix) {
output << "insert " << name() << ":MEmass1 "
<< ix << " " << MEmass1_[ix]/GeV << "\n";
}
// second outgoing mass
for(unsigned int ix=0;ix<MEmass2_.size();++ix) {
output << "insert " << name() << ":MEmass2 "
<< ix << " " << MEmass2_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<decayModes_.size();++ix) {
output << "insert " << name() << ":DecayModes "
<< ix << " " << decayTags_[ix] << " \n";
}
// data for the interpolation tables
std::streamsize curpre=output.precision();
output.precision(curpre+2);
for(unsigned int ix=0;ix<interMasses_.size();++ix) {
output << "insert " << name()
<< ":InterpolationMasses "
<< ix << " " << interMasses_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<interWidths_.size();++ix) {
output << "insert " << name()
<< ":InterpolationWidths "
<< ix << " " << interWidths_[ix]/GeV << "\n";
}
output.precision(curpre);
for(unsigned int ix=0;ix<noOfEntries_.size();++ix) {
output << "insert " << name()
<< ":NumberofEntries "
<< ix << " " << noOfEntries_[ix] << "\n";
}
if(header) output << "\n\" where BINARY ThePEGName=\"" << name()
<< "\";" << endl;
}
DecayMap GenericWidthGenerator::rate(const Particle & p) {
// return default if not using running widths
if(!particle_->variableRatio()) return p.data().decaySelector();
// use the running widths to generate the branching ratios
Energy scale(p.mass());
DecayMap dm;
Energy width = particle_->width();
for(unsigned int ix=0;ix<decayModes_.size();++ix) {
dm.insert(partialWidth(ix,scale)/width,
p.id()==particle_->id() ?
decayModes_[ix] : decayModes_[ix]->CC());
}
return dm;
}
void GenericWidthGenerator::setupMode(tcDMPtr, tDecayIntegratorPtr,
unsigned int)
{}
Energy GenericWidthGenerator::partialWidth(int imode,Energy q) const {
if(q<minMass_[imode]) return ZERO;
Energy gamma;
if(MEtype_[imode]==0) {
gamma=MEcoupling_[imode]*particle_->width();
}
else if(MEtype_[imode]==1) {
gamma=partial2BodyWidth(imode,q);
}
else if(MEtype_[imode]==2) {
gamma=MEcoupling_[imode]*(*interpolators_[imode])(q);
}
else {
throw Exception() << "Unknown type of mode " << MEtype_[imode]
<< "in GenericWidthGenerator::partialWidth()"
<< Exception::runerror;
}
return max(gamma,ZERO);
}
void GenericWidthGenerator::dofinish() {
if(output_) {
string fname = CurrentGenerator::current().filename() +
string("-") + name() + string(".output");
ofstream output(fname.c_str());
dataBaseOutput(output,true);
}
WidthGenerator::dofinish();
}
void GenericWidthGenerator::rebind(const TranslationMap & trans) {
particle_ = trans.translate(particle_);
WidthGenerator::rebind(trans);
}
IVector GenericWidthGenerator::getReferences() {
IVector ret = WidthGenerator::getReferences();
ret.push_back(particle_);
return ret;
}
Length GenericWidthGenerator::lifeTime(const ParticleData &, Energy m, Energy w) const {
if(m<particle_->massMin()) w = width(*particle_,particle_->massMin());
else if(m>particle_->massMax()) w = width(*particle_,particle_->massMax());
return UseRandom::rndExp(hbarc/w);
}
Energy GenericWidthGenerator::partial2BodyWidth(int imode, Energy q,Energy m1,
Energy m2) const {
using Constants::pi;
if(q<m1+m2) return ZERO;
// calcluate the decay momentum
Energy2 q2(q*q),m02(mass_*mass_),m12(m1*m1),m22(m2*m2),
pcm2(0.25*(q2*(q2-2.*m12-2.*m22)+(m12-m22)*(m12-m22))/q2);
if(MEcode_[imode]==-1) return q/mass_*sqr(MEcoupling_[imode])*particle_->width();
Energy pcm(sqrt(pcm2));
double gam(0.);
switch(MEcode_[imode]) {
// V -> P P
case 0: gam = pcm2/6./q2;
break;
// V -> P V
case 1: gam = pcm2/12./m02;
break;
// V -> f fbar
case 2: gam = 1./12.*(q2*(2.*q2-m12-m22+6.*m1*m2)
-(m12-m22)*(m12-m22))/q2/q2;
break;
// P -> VV
case 3: gam = 0.25*pcm2/m02;
break;
// A -> VP
case 4: gam = (2.*pcm2+3.*m12)/24./m02;
break;
// V -> VV
case 5: gam = pcm2/3./q2*(1.+m12/q2+m22/q2);
break;
// S -> SS
case 6: gam = 0.125/q2*m02;
break;
// T -> PP
case 7: gam = pcm2*pcm2/60./q2/m02;
break;
// T -> VP
case 8: gam = pcm2*pcm2/40./m02/m02;
break;
// T -> VV
case 9: gam = 1./30./q2/q2/m02*
(3.*q2*(8.*pcm2*pcm2+5.*(m12*m22+pcm2*(m12+m22)))
-5.*(m12-m22)*(m12-m22)*pcm2);
break;
// P -> PV
case 10: gam = 0.5*pcm2/m22;
break;
// P -> PT
case 11: gam = sqr(pcm2)/12.*q2/m12/m12/m02;
break;
// S -> VV
case 12: gam = 0.125*(2.*pcm2+3.*m12*m22/q2)/m02;
break;
// unknown
default:
throw Exception() << "Unknown type of mode " << MEcode_[imode]
<< " in GenericWidthGenerator::partial2BodyWidth() "
<< Exception::abortnow;
}
return gam*pcm*sqr(MEcoupling_[imode])/pi;
}
pair<Energy,Energy>
GenericWidthGenerator::width(Energy m, const ParticleData & ) const {
pair<Energy,Energy> gamma(make_pair(ZERO,ZERO));
for(unsigned int ix=0;ix<decayModes_.size();++ix) {
if(modeOn_[ix]) {
Energy partial = partialWidth(ix,m);
gamma.second += partial;
if(decayModes_[ix]->on())
gamma.first += partial;
}
}
gamma.first *= prefactor_;
gamma.second *= prefactor_;
return gamma;
}
diff --git a/README b/README
--- a/README
+++ b/README
@@ -1,14 +1,14 @@
========
Herwig 7
========
-This is the release of Herwig 7.1.0, a multi purpose event
+This is the release of Herwig 7.1.1, a multi purpose event
generator for high energy physics.
The Herwig++ distribution contains an adapted version of LoopTools 2.6
<http://www.feynarts.de/looptools/>.
BUILD AND INSTALL
=================
Please refer to https://herwig.hepforge.org/tutorials/ for detailed
instructions.
diff --git a/Sampling/MonacoSampler.h b/Sampling/MonacoSampler.h
--- a/Sampling/MonacoSampler.h
+++ b/Sampling/MonacoSampler.h
@@ -1,195 +1,201 @@
// -*- C++ -*-
//
// MonacoSampler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_MonacoSampler_H
#define Herwig_MonacoSampler_H
//
// This is the declaration of the MonacoSampler class.
//
#include "Herwig/Sampling/BinSampler.h"
#include "Herwig/Utilities/XML/Element.h"
+// work around a Boost 1.64 bug where ublas headers would fail otherwise
+#include <boost/version.hpp>
+#if (BOOST_VERSION / 100 >= 1064)
+#include <boost/serialization/array_wrapper.hpp>
+#endif
+
#include <boost/numeric/ublas/matrix.hpp>
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Michael Rauch
*
* \brief MonacoSampler samples XCombs bins using using Monaco
*
* @see \ref MonacoSamplerInterfaces "The interfaces"
* defined for MonacoSampler.
*/
class MonacoSampler: public Herwig::BinSampler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MonacoSampler();
/**
* The destructor.
*/
virtual ~MonacoSampler();
//@}
public:
/**
* Clone this object.
*/
Ptr<MonacoSampler>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MonacoSampler>::ptr>(clone());
}
public:
/**
* Generate the next point and return its weight; store the point in
* lastPoint().
*/
virtual double generate();
/**
* Adapt this sampler after an iteration has been run
*/
virtual void adapt();
/**
* Return true, if grid data exists for this sampler.
*/
virtual bool existsGrid() const;
/**
* Save grid data
*/
virtual void saveGrid() const;
/**
* Initialize this bin sampler. This default version calls runIteration.
*/
virtual void initialize(bool progress);
/**
* Finalize this sampler.
*/
virtual void finalize(bool);
/**
* Fill Monaco grid data from an XML element
*/
virtual void fromXML(const XML::Element&);
/**
* Return an XML element for the data of the Monaco grid
*/
virtual XML::Element toXML() 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 Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* Rate of grid modification (0 for no modification)
*/
double theAlpha;
/**
* Number of grid divisions per dimension
*/
size_t theGridDivisions;
/**
* Grid boundaries
* (first index: dimension of random numbers,
* second index: dimension of partitions per random number)
*/
boost::numeric::ublas::matrix<double> theGrid;
/**
* Collected value per grid bin
*/
boost::numeric::ublas::matrix<double> theGridData;
private:
/**
* Number of points collected in iteration so far
*/
size_t theIterationPoints;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MonacoSampler & operator=(const MonacoSampler &);
};
}
#endif /* Herwig_MonacoSampler_H */
diff --git a/Shower/Dipole/Base/DipoleEventRecord.cc b/Shower/Dipole/Base/DipoleEventRecord.cc
--- a/Shower/Dipole/Base/DipoleEventRecord.cc
+++ b/Shower/Dipole/Base/DipoleEventRecord.cc
@@ -1,1531 +1,1537 @@
// -*- C++ -*-
//
// DipoleEventRecord.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleEventRecord class.
//
#include "DipoleEventRecord.h"
#include "Herwig/Shower/Dipole/Utility/DipolePartonSplitter.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Decay/HwDecayerBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include <boost/utility.hpp>
#include <algorithm>
#include <iterator>
using namespace Herwig;
PList DipoleEventRecord::colourOrdered(PPair & in,
PList & out) {
PList colour_ordered;
size_t done_size = out.size();
if (in.first->coloured())
++done_size;
if (in.second && in.second->coloured())
++done_size;
while (colour_ordered.size() != done_size) {
PPtr current;
// start with singlets, as long as we have some
if (find(colour_ordered.begin(),colour_ordered.end(),in.first) ==
colour_ordered.end() && in.first->coloured()) {
if (!in.first->hasColour() || !in.first->hasAntiColour())
current = in.first;
}
if (!current) {
for (PList::iterator p = out.begin();
p != out.end(); ++p) {
if (find(colour_ordered.begin(),colour_ordered.end(),*p) ==
colour_ordered.end() && (**p).coloured()) {
if (!(**p).hasColour() || !(**p).hasAntiColour()) {
current = *p;
break;
}
}
}
}
if (!current) {
if (in.second && find(colour_ordered.begin(),colour_ordered.end(),in.second) ==
colour_ordered.end() && in.second->coloured()) {
if (!in.second->hasColour() || !in.second->hasAntiColour())
current = in.second;
}
}
// then go on with anything else
if (!current) {
if (find(colour_ordered.begin(),colour_ordered.end(),in.first) ==
colour_ordered.end() && in.first->coloured()) {
current = in.first;
}
}
if (!current) {
for (PList::iterator p = out.begin();
p != out.end(); ++p) {
if (find(colour_ordered.begin(),colour_ordered.end(),*p) ==
colour_ordered.end() && (**p).coloured()) {
current = *p;
break;
}
}
}
if (!current) {
if (in.second && find(colour_ordered.begin(),colour_ordered.end(),in.second) ==
colour_ordered.end() && in.second->coloured()) {
current = in.second;
}
}
assert(current);
PPtr next;
Ptr<ColourLine>::ptr walk_the_line;
while (true) {
if (!walk_the_line) {
if (current->hasColour()) {
walk_the_line = current->colourLine();
}
else if (current->hasAntiColour()) {
walk_the_line = current->antiColourLine();
}
}
if (!next)
for (tPVector::const_iterator p = walk_the_line->coloured().begin();
p != walk_the_line->coloured().end(); ++p) {
if (*p == current)
continue;
if (find(out.begin(),out.end(),*p) != out.end() ||
*p == in.first ||
(in.second && *p == in.second)) {
next = *p;
if (next->hasColour() && next->hasAntiColour()) {
walk_the_line = walk_the_line == next->colourLine() ? next->antiColourLine() : next->colourLine();
}
break;
}
}
if (!next)
for (tPVector::const_iterator p = walk_the_line->antiColoured().begin();
p != walk_the_line->antiColoured().end(); ++p) {
if (*p == current)
continue;
if (find(out.begin(),out.end(),*p) != out.end() ||
*p == in.first ||
(in.second && *p == in.second)) {
next = *p;
if (next->hasColour() && next->hasAntiColour()) {
walk_the_line = walk_the_line == next->colourLine() ? next->antiColourLine() : next->colourLine();
}
break;
}
}
assert(next);
colour_ordered.push_back(current);
current = next;
// done if next is not a gluon or next is already in colour_ordered
if ((current->hasColour() && !current->hasAntiColour()) ||
(!current->hasColour() && current->hasAntiColour())) {
colour_ordered.push_back(current);
break;
}
if (next->hasColour() && next->hasAntiColour()) {
if (find(colour_ordered.begin(),colour_ordered.end(),next) != colour_ordered.end())
break;
}
next = PPtr();
}
}
return colour_ordered;
}
void DipoleEventRecord::popChain() {
assert(!theChains.empty());
theDoneChains.push_back(DipoleChain());
theDoneChains.back().dipoles().splice(theDoneChains.back().dipoles().begin(),theChains.front().dipoles());
theChains.pop_front();
}
void DipoleEventRecord::popChain(list<DipoleChain>::iterator ch) {
assert(!theChains.empty());
theDoneChains.push_back(DipoleChain());
theDoneChains.back().dipoles().splice(theDoneChains.back().dipoles().begin(),ch->dipoles());
theChains.erase(ch);
}
void DipoleEventRecord::popChains(const list<list<DipoleChain>::iterator>& chs) {
assert(!theChains.empty());
for ( list<list<DipoleChain>::iterator>::const_iterator ch =
chs.begin(); ch != chs.end(); ++ch ) {
theDoneChains.push_back(DipoleChain());
theDoneChains.back().dipoles().splice(theDoneChains.back().dipoles().begin(),(*ch)->dipoles());
}
for ( list<list<DipoleChain>::iterator>::const_iterator ch =
chs.begin(); ch != chs.end(); ++ch )
theChains.erase(*ch);
}
DipoleIndex
DipoleEventRecord::mergeIndex(list<Dipole>::iterator firstDipole, const pair<bool,bool>& whichFirst,
list<Dipole>::iterator secondDipole, const pair<bool,bool>& whichSecond) const {
tcPDPtr emitterData =
whichFirst.first ? firstDipole->leftParticle()->dataPtr() : firstDipole->rightParticle()->dataPtr();
tcPDPtr spectatorData =
whichSecond.first ? secondDipole->leftParticle()->dataPtr() : secondDipole->rightParticle()->dataPtr();
const PDF& emitterPDF =
whichFirst.first ? firstDipole->leftPDF() : firstDipole->rightPDF();
const PDF& spectatorPDF =
whichSecond.first ? secondDipole->leftPDF() : secondDipole->rightPDF();
return DipoleIndex(emitterData,spectatorData,emitterPDF,spectatorPDF);
}
SubleadingSplittingInfo
DipoleEventRecord::mergeSplittingInfo(list<DipoleChain>::iterator firstChain, list<Dipole>::iterator firstDipole,
const pair<bool,bool>& whichFirst,
list<DipoleChain>::iterator secondChain, list<Dipole>::iterator secondDipole,
const pair<bool,bool>& whichSecond) const {
SubleadingSplittingInfo res;
res.index(mergeIndex(firstDipole,whichFirst,secondDipole,whichSecond));
res.emitter(whichFirst.first ? firstDipole->leftParticle() : firstDipole->rightParticle());
res.spectator(whichSecond.first ? secondDipole->leftParticle() : secondDipole->rightParticle());
res.emitterX(whichFirst.first ? firstDipole->leftFraction() : firstDipole->rightFraction());
res.spectatorX(whichSecond.first ? secondDipole->leftFraction() : secondDipole->rightFraction());
res.configuration(whichFirst);
res.spectatorConfiguration(whichSecond);
res.emitterChain(firstChain);
res.emitterDipole(firstDipole);
res.spectatorChain(secondChain);
res.spectatorDipole(secondDipole);
return res;
}
void DipoleEventRecord::getSubleadingSplittings(list<SubleadingSplittingInfo>& res) {
static pair<bool,bool> left(true,false);
static pair<bool,bool> right(false,true);
res.clear();
for ( list<DipoleChain>::iterator cit = theChains.begin();
cit != theChains.end(); ++cit ) {
for ( list<Dipole>::iterator dit = cit->dipoles().begin();
dit != cit->dipoles().end(); ++dit ) {
for ( list<Dipole>::iterator djt = dit;
djt != cit->dipoles().end(); ++djt ) {
res.push_back(mergeSplittingInfo(cit,dit,left,cit,djt,left));
res.push_back(mergeSplittingInfo(cit,dit,right,cit,djt,right));
if ( dit != djt ) {
res.push_back(mergeSplittingInfo(cit,dit,left,cit,djt,right));
res.push_back(mergeSplittingInfo(cit,dit,right,cit,djt,left));
}
}
}
list<DipoleChain>::iterator cjt = cit; ++cjt;
for ( ; cjt != theChains.end(); ++cjt ) {
for ( list<Dipole>::iterator dit = cit->dipoles().begin();
dit != cit->dipoles().end(); ++dit ) {
for ( list<Dipole>::iterator djt = cjt->dipoles().begin();
djt != cjt->dipoles().end(); ++djt ) {
res.push_back(mergeSplittingInfo(cit,dit,left,cjt,djt,left));
res.push_back(mergeSplittingInfo(cit,dit,right,cjt,djt,right));
res.push_back(mergeSplittingInfo(cit,dit,left,cjt,djt,right));
res.push_back(mergeSplittingInfo(cit,dit,right,cjt,djt,left));
}
}
}
}
}
void DipoleEventRecord::splitSubleading(SubleadingSplittingInfo& dsplit,
pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
DipoleChain*& firstChain, DipoleChain*& secondChain) {
if ( dsplit.emitterDipole() == dsplit.spectatorDipole() ) {
assert(dsplit.emitterChain() == dsplit.spectatorChain());
split(dsplit.emitterDipole(),dsplit.emitterChain(),dsplit,
childIterators,firstChain,secondChain,false);
} else {
// first need to recoil, then split
recoil(dsplit.spectatorDipole(),dsplit.spectatorChain(),dsplit);
split(dsplit.emitterDipole(),dsplit.emitterChain(),dsplit,
childIterators,firstChain,secondChain,true);
}
}
void DipoleEventRecord::findChains(const PList& ordered, const bool decay) {
// All uses of findChains should guarantee
// a non-empty list of particles
assert( !ordered.empty() );
theChains.clear();
theDoneChains.clear();
DipoleChain current_chain;
// this whole thing needs to have a more elegant implementation at some point
bool startIsTriplet =
(ordered.front()->hasColour() && !ordered.front()->hasAntiColour()) ||
(!ordered.front()->hasColour() && ordered.front()->hasAntiColour());
bool endIsTriplet =
(ordered.back()->hasColour() && !ordered.back()->hasAntiColour()) ||
(!ordered.back()->hasColour() && ordered.back()->hasAntiColour());
if (!( ordered.size() == 2 && startIsTriplet && endIsTriplet)) {
PList::const_iterator theStart = ordered.begin();
bool onceMore = false;
for (PList::const_iterator p = ordered.begin();
p != ordered.end(); ++p) {
PList::const_iterator next_it =
p != --ordered.end() ? std::next(p) : ordered.begin();
if (!DipolePartonSplitter::colourConnected(*p,*next_it)) {
// it may have happened that we need to close the chain due to another
// chain starting right now; see the above global comment for this fix
bool startIsOctet =
(**theStart).hasColour() && (**theStart).hasAntiColour();
bool endIsOctet =
(**p).hasColour() && (**p).hasAntiColour();
if ( DipolePartonSplitter::colourConnected(*p,*theStart) &&
startIsOctet && endIsOctet ) {
swap(next_it,theStart);
onceMore = true;
} else {
- theStart = next_it;
- current_chain.check();
- theChains.push_back(current_chain);
- current_chain.dipoles().clear();
- continue;
+ theStart = next_it;
+ current_chain.check();
+ // Randomize the chains agains biasing of directions.
+ if(UseRandom::rndbool()) theChains.push_back(current_chain);
+ else theChains.insert(theChains.begin(),current_chain);
+ current_chain.dipoles().clear();
+ continue;
}
}
pair<bool,bool> initial_state (false,false);
initial_state.first = (*p == incoming().first || *p == incoming().second);
initial_state.second = (*next_it == incoming().first || *next_it == incoming().second);
pair<int,int> which_in (-1,-1);
if (initial_state.first)
which_in.first = *p == incoming().first ? 0 : 1;
if (initial_state.second)
which_in.second = *next_it == incoming().first ? 0 : 1;
pair<double,double> xs (1.,1.);
if (initial_state.first)
xs.first = *p == incoming().first ? fractions().first : fractions().second;
if (initial_state.second)
xs.second = *next_it == incoming().first ? fractions().first : fractions().second;
pair<PDF,PDF> pdf;
if ( which_in.first == 0 )
pdf.first = pdfs().first;
else if ( which_in.first == 1 )
pdf.first = pdfs().second;
if ( which_in.second == 0 )
pdf.second = pdfs().first;
else if ( which_in.second == 1 )
pdf.second = pdfs().second;
// In the case of a decay process register which
// parton is incoming to the decay
pair<bool,bool> decayed_parton (false,false);
if (decay) {
decayed_parton.first = (*p == currentDecay()->incoming()[0].first);
decayed_parton.second = (*next_it == currentDecay()->incoming()[0].first);
}
current_chain.dipoles().push_back(Dipole({*p,*next_it},pdf,xs,decayed_parton));
if ( onceMore ) {
- next_it = theStart;
- current_chain.check();
- theChains.push_back(current_chain);
- current_chain.dipoles().clear();
- onceMore = false;
+ next_it = theStart;
+ current_chain.check();
+ // Randomize the chains agains biasing of directions.
+ if(UseRandom::rndbool()) theChains.push_back(current_chain);
+ else theChains.insert(theChains.begin(),current_chain);
+ current_chain.dipoles().clear();
+ onceMore = false;
}
}
} else {
// treat 2 -> singlet, singlet -> 2 and 1 + singlet -> 1 + singlet special
// to prevent duplicate dipole
assert(DipolePartonSplitter::colourConnected(ordered.front(),ordered.back()));
pair<bool,bool> initial_state (false,false);
initial_state.first = (ordered.front() == incoming().first || ordered.front() == incoming().second);
initial_state.second = (ordered.back() == incoming().first || ordered.back() == incoming().second);
pair<int,int> which_in (-1,-1);
if (initial_state.first)
which_in.first = ordered.front() == incoming().first ? 0 : 1;
if (initial_state.second)
which_in.second = ordered.back() == incoming().first ? 0 : 1;
pair<double,double> xs (1.,1.);
if (initial_state.first)
xs.first = ordered.front() == incoming().first ? fractions().first : fractions().second;
if (initial_state.second)
xs.second = ordered.back() == incoming().first ? fractions().first : fractions().second;
pair<PDF,PDF> pdf;
if ( which_in.first == 0 )
pdf.first = pdfs().first;
else if ( which_in.first == 1 )
pdf.first = pdfs().second;
if ( which_in.second == 0 )
pdf.second = pdfs().first;
else if ( which_in.second == 1 )
pdf.second = pdfs().second;
// In the case of a decay process register which
// parton is incoming to the decay
pair<bool,bool> decayed_parton (false,false);
if (decay) {
decayed_parton.first = (ordered.front() == currentDecay()->incoming()[0].first);
decayed_parton.second = (ordered.back() == currentDecay()->incoming()[0].first);
}
current_chain.dipoles().push_back(Dipole({ordered.front(),ordered.back()},pdf,xs,decayed_parton));
}
if (!current_chain.dipoles().empty()) {
current_chain.check();
- theChains.push_back(current_chain);
+ // Randomize the chains agains biasing of directions.
+ if(UseRandom::rndbool()) theChains.push_back(current_chain);
+ else theChains.insert(theChains.begin(),current_chain);
}
}
const map<PPtr,PPtr>&
DipoleEventRecord::prepare(tSubProPtr subpro,
tStdXCombPtr xc,
StepPtr step,
const pair<PDF,PDF>& pdf,tPPair beam,
bool firstInteraction, bool dipoles) {
// set the subprocess
subProcess(subpro);
// clear the event record
outgoing().clear();
theHard.clear();
theOriginals.clear();
theDecays.clear();
theCurrentDecay = PerturbativeProcessPtr();
// extract incoming particles
PPair in = subpro->incoming();
// get the incoming momentum fractions
// don't take these from the XComb as it may be null
pair<double,double> xs;
ThePEG::Direction<0> dir(true);
xs.first = in.first->momentum().dirPlus()/beam.first->momentum().dirPlus();
dir.reverse();
xs.second = in.second->momentum().dirPlus()/beam.second->momentum().dirPlus();
xcombPtr(xc);
pdfs() = pdf;
fractions() = xs;
// use ShowerHandler to split up the hard process
PerturbativeProcessPtr hard;
DecayProcessMap decay;
// Special handling for the first interaction:
// If a post subprocess handler (e.g. QED radiation)
// is applied, there may be particles in the step object not
// present in the subprocess object (other than any remnants).
// These need to be included in any transformations due to
// II splittings in ::update.
if ( firstInteraction ) {
// Initialise a PVector for the outgoing
tPVector hardProcOutgoing;
// Include all outgoing particles that are not remnants
for ( auto & part : step->particles() )
if ( part->id() != 82 ) {
hardProcOutgoing.push_back(part);
}
ShowerHandler::currentHandler()->splitHardProcess(hardProcOutgoing,
hard, decay);
}
// For secondary collisions we must use the
// subProcess object and not the step as the
// step stores all outgoing from the entire collision
else
ShowerHandler::currentHandler()->splitHardProcess(tPVector(subpro->outgoing().begin(),
subpro->outgoing().end()),
hard,decay);
// vectors for originals and copies of the particles
vector<PPtr> original;
vector<PPtr> copies;
// fill originals
for(unsigned int ix=0;ix<2;++ix)
original.push_back(hard->incoming()[ix].first);
for(unsigned int ix=0;ix<hard->outgoing().size();++ix)
original.push_back(hard->outgoing()[ix].first);
for(DecayProcessMap::const_iterator it=decay.begin();it!=decay.end();++it) {
fillFromDecays(it->second, original);
}
// and make copies
for ( vector<PPtr>::const_iterator p = original.begin();
p != original.end(); ++p ) {
PPtr copy = new_ptr(Particle(**p));
copies.push_back(copy);
theOriginals[*p] = copy;
}
// isolate the colour of the copies from the originals
colourIsolate(original,copies);
// set the incoming particles
incoming().first = copies[0];
ParticleVector children = incoming().first->children();
for ( ParticleVector::const_iterator c = children.begin();
c != children.end(); ++c )
incoming().first->abandonChild(*c);
incoming().second = copies[1];
children = incoming().second->children();
for ( ParticleVector::const_iterator c = children.begin();
c != children.end(); ++c )
incoming().second->abandonChild(*c);
// set the outgoing particles for the hard process
for(unsigned int ix=0;ix<hard->outgoing().size();++ix) {
if(hard->outgoing()[ix].first->coloured())
outgoing().push_back(theOriginals[hard->outgoing()[ix].first]);
else
theHard.push_back(theOriginals[hard->outgoing()[ix].first]);
}
if ( dipoles ) {
PList cordered = colourOrdered(incoming(),outgoing());
if ( !cordered.empty() )
findChains(cordered,false);
}
// sort out the decays
for(auto const & dec : decay) {
// If the decay particle is in original it needs
// to be added to the decays and the decay needs to be
// changed to the copied particles.
if ( theOriginals.find(dec.second->incoming()[0].first) != theOriginals.end() ) {
theDecays[theOriginals[dec.second->incoming()[0].first]] = dec.second;
PerturbativeProcessPtr decayProc = theDecays[theOriginals[dec.second->incoming()[0].first]];
separateDecay(decayProc);
}
else {
assert( find( copies.begin(), copies.end(), dec.second->incoming()[0].first ) != copies.end() );
theDecays[dec.second->incoming()[0].first] = dec.second;
}
}
PList::const_iterator XFirst, XLast;
if ( !theHard.empty() ) {
XFirst = theHard.begin();
XLast = theHard.end();
} else {
XFirst = outgoing().begin();
XLast = outgoing().end();
}
thePX = (**XFirst).momentum();
++XFirst;
for ( ; XFirst != XLast; ++XFirst )
thePX += (**XFirst).momentum();
identifyEventType();
return theOriginals;
}
void DipoleEventRecord::slimprepare(tSubProPtr subpro,
tStdXCombPtr xc,
const pair<PDF,PDF>& pdf,tPPair beam,
bool dipoles) {
// set the subprocess
subProcess(subpro);
// clear the event record
outgoing().clear();
theHard.clear();
theOriginals.clear();
theDecays.clear();
theCurrentDecay = PerturbativeProcessPtr();
// extract incoming particles
PPair in = subpro->incoming();
// get the beam
// get the incoming momentum fractions
// don't take these from the XComb as it may be null
pair<double,double> xs;
ThePEG::Direction<0> dir(true);
xs.first = in.first->momentum().dirPlus()/beam.first->momentum().dirPlus();
dir.reverse();
xs.second = in.second->momentum().dirPlus()/beam.second->momentum().dirPlus();
xcombPtr(xc);
pdfs() = pdf;
fractions() = xs;
incoming() = in;
for(unsigned int ix=0;ix<subpro->outgoing().size();++ix) {
if(subpro->outgoing()[ix]->coloured())
outgoing().push_back(subpro->outgoing()[ix]);
}
if ( dipoles ) {
PList cordered = colourOrdered(incoming(),outgoing());
if ( !cordered.empty() )
findChains(cordered,false);
}
}
void DipoleEventRecord::fillFromDecays(PerturbativeProcessPtr decayProc, vector<PPtr>& original) {
// Loop over the outgoing of the given perturbative process
for ( auto const & outIt : decayProc->outgoing() ) {
// Add the outgoing particle to the vector of original particles
original.push_back(outIt.first);
// Iterate through the outgoing
if ( outIt.second )
fillFromDecays( outIt.second, original);
}
}
void DipoleEventRecord::separateDecay(PerturbativeProcessPtr decayProc) {
// Iteratively replace all entries in the incoming
// with their copies.
for ( auto & inIt : decayProc->incoming() ) {
if ( theOriginals.find( inIt.first ) != theOriginals.end() )
inIt.first = theOriginals[inIt.first];
}
// Iteratively replace all entries in the outgoing
// with their copies.
for ( auto & outIt : decayProc->outgoing()) {
if ( theOriginals.count( outIt.first ) )
outIt.first = theOriginals[outIt.first];
if ( outIt.second )
separateDecay(outIt.second);
}
}
void DipoleEventRecord::clear() {
ShowerEventRecord::clear();
theDecays.clear();
theHard.clear();
theChains.clear();
theDoneChains.clear();
theOriginals.clear();
}
pair<PVector,PVector> DipoleEventRecord::tmpupdate(DipoleSplittingInfo& dsplit) {
PVector inc;
PVector out;
tcPPtr IF = incoming().first;
tcPPtr IS = incoming().second;
tcPPtr DE = dsplit.emitter();
tcPPtr DS = dsplit.spectator();
if ( IF != DE && IF != DS ) {
PPtr p = IF->data().produceParticle(IF->momentum());
inc.push_back(p);
}
else if ( IF == DE ) inc.push_back( dsplit.splitEmitter() );
else if ( IF == DS ) inc.push_back( dsplit.splitSpectator() );
if ( IS != DE && IS != DS ) {
PPtr p = IS->data().produceParticle(IS->momentum());
inc.push_back(p);
}
else if ( IS == DE ) inc.push_back( dsplit.splitEmitter() );
else if ( IS == DS ) inc.push_back( dsplit.splitSpectator() );
if ( IF != DE && IS != DE)
out.push_back( dsplit.splitEmitter());
if ( IF != DS && IS != DS)
out.push_back( dsplit.splitSpectator());
out.push_back( dsplit.emission());
for ( tcPPtr h : theHard ){
PPtr p = h->data().produceParticle(h->momentum());
if ( dsplit.splittingKinematics()->doesTransform() ) {
p->set5Momentum( dsplit.splittingKinematics()->transform(p->momentum()) );
}
out.push_back(p);
}
for ( tcPPtr p : outgoing() )
if ( p != DE &&
p != DS &&
p != dsplit.emission() ){
PPtr ou = p->data().produceParticle(p->momentum());;
if ( dsplit.splittingKinematics()->doesTransform() ){
ou->set5Momentum( dsplit.splittingKinematics()->transform(ou->momentum()) );
}
out.push_back(ou);
}
return {inc,out};
}
void DipoleEventRecord::update(DipoleSplittingInfo& dsplit) {
if ( incoming().first == dsplit.emitter() ) {
intermediates().push_back(dsplit.emitter());
incoming().first = dsplit.splitEmitter();
fractions().first /= dsplit.lastEmitterZ();
} else if ( incoming().first == dsplit.spectator() ) {
intermediates().push_back(dsplit.spectator());
incoming().first = dsplit.splitSpectator();
fractions().first /= dsplit.lastSpectatorZ();
}
if ( incoming().second == dsplit.emitter() ) {
intermediates().push_back(dsplit.emitter());
incoming().second = dsplit.splitEmitter();
fractions().second /= dsplit.lastEmitterZ();
} else if ( incoming().second == dsplit.spectator() ) {
intermediates().push_back(dsplit.spectator());
incoming().second = dsplit.splitSpectator();
fractions().second /= dsplit.lastSpectatorZ();
}
PList::iterator pos;
pos = find(outgoing().begin(), outgoing().end(), dsplit.emitter());
if (pos != outgoing().end()) {
intermediates().push_back(*pos);
*pos = dsplit.splitEmitter();
}
pos = find(outgoing().begin(), outgoing().end(), dsplit.spectator());
if (pos != outgoing().end()) {
intermediates().push_back(*pos);
*pos = dsplit.splitSpectator();
}
outgoing().push_back(dsplit.emission());
if (dsplit.splittingKinematics()->doesTransform()) {
for (PList::iterator p = intermediates().begin();
p != intermediates().end(); ++p) {
(**p).set5Momentum(dsplit.splittingKinematics()->transform((**p).momentum()));
}
for (PList::iterator h = theHard.begin();
h != theHard.end(); ++h) {
(**h).set5Momentum(dsplit.splittingKinematics()->transform((**h).momentum()));
}
for (PList::iterator p = outgoing().begin();
p != outgoing().end(); ++p) {
if ((*p) != dsplit.splitEmitter() &&
(*p) != dsplit.splitSpectator() &&
(*p) != dsplit.emission())
(**p).set5Momentum(dsplit.splittingKinematics()->transform((**p).momentum()));
}
}
// Handle updates related to decays
// Showering of decay processes
// Treat the evolution of the incoming
// decayed particle as in backward evolution
if ( dsplit.isDecayProc() ) {
// Create a pointer to the decay process
PerturbativeProcessPtr decayProc = currentDecay();
// Add the emission to the outgoing of the decay process
decayProc->outgoing().push_back( {dsplit.emission(), PerturbativeProcessPtr() });
// Bools to be used throughout
const bool decayedEmtr = dsplit.index().incomingDecayEmitter();
const bool decayedSpec = dsplit.index().incomingDecaySpectator();
/*
In the current implementation, **following the hard process**
all particles in theDecays evolve independently
e.g. if we have W -> XYZ where all X, Y and Z need to be
showered and decayed, we only identify them as needing decaying
(and hence put them in theDecays) AFTER showering the decay of W.
Hence, XYZ are not even in theDecays until W has been fully
showered and then they are decayed and showered completely independently
KEY POINT - Never need to update other entries of theDecays
Note: The PPtr in theDecays should remain unchanged and all changes
should be made to the relative PerturbativeProcess.
*/
// Splittings from dipoles in the decay process which
// do not have the decayed parton as emitter or spectator.
// Update the decay process in theDecays
if ( !decayedEmtr && !decayedSpec ) {
// Find and replace the old spectator and
// emitter in the outgoing of the decay process
bool decayProcEm = false;
bool decayProcSp = false;
for ( auto & outIt : decayProc->outgoing() ) {
if ( !decayProcEm && outIt.first == dsplit.emitter() ) {
outIt = {dsplit.splitEmitter(), PerturbativeProcessPtr()};
decayProcEm = true;
}
if ( !decayProcSp && outIt.first == dsplit.spectator() ) {
outIt = {dsplit.splitSpectator(), PerturbativeProcessPtr() };
decayProcSp = true;
}
if ( decayProcEm && decayProcSp )
break;
}
// Test that nothing strange is happening
assert( (decayProcEm && decayProcSp) );
return;
}
// The spectator is the decayed particle
else if ( decayedSpec ) {
// Update the dipole event record intermediates
intermediates().push_back(dsplit.splitSpectator());
// Update the the decayProcess incoming
decayProc->incoming().clear();
decayProc->incoming().push_back({dsplit.splitSpectator(),decayProc});
// Update the decay process outgoing
// Replace the old emitter with the new emitter
for ( auto & outEmtrIt : decayProc->outgoing() ) {
if ( outEmtrIt.first == dsplit.emitter() ){
outEmtrIt = {dsplit.splitEmitter(), PerturbativeProcessPtr() };
break;
}
}
// Perform the recoil transformation
// Find all particles in the recoil system
PList recoilSystem;
for ( auto const & outIt : decayProc->outgoing() ) {
if ( outIt.first != dsplit.splitEmitter() && outIt.first != dsplit.emission() ) {
recoilSystem.push_back(outIt.first);
}
}
dsplit.splittingKinematics()->decayRecoil( recoilSystem );
return;
}
// The emitter is the decayed particle
else {
throw Exception()
<< "DipoleEventRecord: The emitter as a decayed particle is currently not implemented."
<< Exception::runerror;
assert( currentDecay()->incoming()[0].first == dsplit.emitter() && decayedEmtr && !decayedSpec );
// Update the dipole event record intermediates
intermediates().push_back(dsplit.splitEmitter());
// Update the the decayProcess incoming
decayProc->incoming().clear();
decayProc->incoming().push_back({dsplit.splitEmitter(),decayProc});
// Update the decay process outgoing
// Replace the old spectator with the new spectator
for (auto & outSpecIt : decayProc->outgoing() ) {
if ( outSpecIt.first == dsplit.spectator() ){
outSpecIt = { dsplit.splitSpectator(), PerturbativeProcessPtr() };
break;
}
}
// Perform the recoil transformation
assert(dsplit.splittingKinematics()->isDecay());
// Find all particles in the recoil system
PList recoilSystem;
for ( auto const & outIt : decayProc->outgoing() ) {
if ( outIt.first != dsplit.splitSpectator() && outIt.first != dsplit.emission() ) {
recoilSystem.push_back(outIt.first);
}
}
dsplit.splittingKinematics()->decayRecoil( recoilSystem );
return;
}
}
}
void
DipoleEventRecord::split(list<Dipole>::iterator dip,
list<DipoleChain>::iterator ch,
DipoleSplittingInfo& dsplit,
pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
DipoleChain*& firstChain, DipoleChain*& secondChain,
bool colourSpectator) {
static DipoleChain empty;
pair<Dipole,Dipole> children = dip->split(dsplit,colourSpectator);
list<Dipole>::iterator breakup =
ch->insertSplitting(dip,children,childIterators);
if ( breakup == ch->dipoles().end() ) {
firstChain = &(*ch);
secondChain = ∅
} else {
DipoleChain other;
other.dipoles().splice(other.dipoles().end(),ch->dipoles(),breakup,ch->dipoles().end());
chains().push_back(other);
firstChain = &(*ch);
secondChain = &(chains().back());
// explicitly fix iterators in case the splice implementation
// at hand does invalidate iterators (the SGI docu says, it doesn't,
// but it seems that this behaviour is not part of the standard)
childIterators.first = --firstChain->dipoles().end();
childIterators.second = secondChain->dipoles().begin(); }
if ( !colourSpectator ) {
update(dsplit); // otherwise done by recoil(...)
}
}
pair<PVector,PVector> DipoleEventRecord::tmpsplit(list<Dipole>::iterator dip,
list<DipoleChain>::iterator ,
DipoleSplittingInfo& dsplit,
pair<list<Dipole>::iterator,list<Dipole>::iterator>& ,
DipoleChain*& , DipoleChain*& ,
bool colourSpectator) {
dip->tmpsplit(dsplit,colourSpectator);
return tmpupdate(dsplit); // otherwise done by recoil(...)
}
void DipoleEventRecord::recoil(list<Dipole>::iterator dip,
list<DipoleChain>::iterator ch,
DipoleSplittingInfo& dsplit) {
dip->recoil(dsplit);
ch->updateDipole(dip);
update(dsplit);
}
list<pair<list<Dipole>::iterator,list<DipoleChain>::iterator> >
DipoleEventRecord::inDipoles() {
list<pair<list<Dipole>::iterator,list<DipoleChain>::iterator> > res;
for ( list<DipoleChain>::iterator chit = theDoneChains.begin();
chit != theDoneChains.end(); ++chit ) {
bool haveOne = false;
for ( list<Dipole>::iterator dit = chit->dipoles().begin();
dit != chit->dipoles().end(); ++dit ) {
if ( dit->leftPDF().pdf() || dit->rightPDF().pdf() ) {
haveOne = true;
break;
}
}
if ( haveOne ) {
theChains.splice(theChains.begin(),theDoneChains,chit);
for ( list<Dipole>::iterator dit = theChains.front().dipoles().begin();
dit != theChains.front().dipoles().end(); ++dit ) {
if ( dit->leftPDF().pdf() || dit->rightPDF().pdf() ) {
res.push_back({dit,theChains.begin()});
}
}
}
}
return res;
}
void DipoleEventRecord::transform(const SpinOneLorentzRotation& rot) {
Lorentz5Momentum tmp;
for (PList::iterator p = intermediates().begin();
p != intermediates().end(); ++p) {
tmp = (**p).momentum(); tmp = rot * tmp;
(**p).set5Momentum(tmp);
}
for (PList::iterator h = theHard.begin();
h != theHard.end(); ++h) {
tmp = (**h).momentum(); tmp = rot * tmp;
(**h).set5Momentum(tmp);
}
for (PList::iterator p = outgoing().begin();
p != outgoing().end(); ++p) {
tmp = (**p).momentum(); tmp = rot * tmp;
(**p).set5Momentum(tmp);
}
}
tPPair DipoleEventRecord::fillEventRecord(StepPtr step, bool firstInteraction, bool) {
PPtr inSubPro = subProcess()->incoming().first;
PPtr inParticle;
if ( !(inSubPro->parents().empty()) )
inParticle = inSubPro->parents()[0];
else
inParticle = inSubPro;
PPtr inParton = theOriginals[inSubPro];
theOriginals.erase(inSubPro);
updateColour(incoming().first,true);
if ( inParticle != inSubPro )
inParticle->abandonChild(inSubPro);
inParton->addChild(inSubPro);
if ( inParticle != inSubPro )
inParticle->addChild(incoming().first);
intermediates().push_back(inSubPro);
intermediates().push_back(inParton);
// Repeat all the above for the second incoming particle
inSubPro = subProcess()->incoming().second;
if ( !(inSubPro->parents().empty()) )
inParticle = inSubPro->parents()[0];
else
inParticle = inSubPro;
inParton = theOriginals[inSubPro];
theOriginals.erase(inSubPro);
updateColour(incoming().second,true);
if ( inParticle != inSubPro )
inParticle->abandonChild(inSubPro);
inParton->addChild(inSubPro);
if ( inParticle != inSubPro )
inParticle->addChild(incoming().second);
intermediates().push_back(inSubPro);
intermediates().push_back(inParton);
// theOriginals is populated in ::prepare and contains all of the incoming and outgoing particles of the original hard process
// Here outgoing particles from theOriginals are added into the intermediates()
while ( !theOriginals.empty() ) {
PPtr outSubPro = theOriginals.begin()->first;
PPtr outParton = theOriginals.begin()->second;
// workaround for OS X Mavericks LLVM libc++
#ifdef _LIBCPP_VERSION
map<PPtr,PPtr>::const_iterator beg = theOriginals.begin();
#else
map<PPtr,PPtr>::iterator beg = theOriginals.begin();
#endif
theOriginals.erase(beg);
updateColour(outParton,true);
outSubPro->addChild(outParton);
intermediates().push_back(outSubPro);
}
// Update the intermediates of the step
step->addIntermediates(intermediates().begin(),intermediates().end());
for (auto const & p : outgoing())
step->addDecayProduct( p );
for (auto const & p : theHard)
step->addDecayProduct( p );
if ( firstInteraction &&
(incoming().first->coloured() ||
incoming().second->coloured() ) ) {
ShowerHandler::currentHandler()->lastExtractor()
->newRemnants(subProcess()->incoming(),incoming(),step);
}
step->addIntermediate(incoming().first);
step->addIntermediate(incoming().second);
return incoming();
}
bool DipoleEventRecord::prepareDecay( PerturbativeProcessPtr decayProc ) {
// Create objects containing the incoming and outgoing partons,
// required as inputs for colourOrdered.
PList out;
for( auto const & dec : decayProc->outgoing()) {
if(dec.first->coloured()) {
out.push_back(dec.first);
}
}
// Only need to shower if we have coloured outgoing particles
if ( out.empty() )
return false;
else {
// For the incoming, use a PPair containing the incoming and a null pointer
PPair in;
in.first = decayProc->incoming()[0].first;
// Create an ordered list of particles
PList cordered;
cordered = colourOrdered(in,out);
// Find the dipole chains for this decay
findChains(cordered,true);
return true;
}
}
Energy DipoleEventRecord::decay(PPtr incoming, bool& powhegEmission) {
// get the process
PerturbativeProcessPtr process = theDecays[incoming];
assert(process);
//tDMPtr decayMode = new_ptr(DecayMode());
tDMPtr decayMode = DMPtr();
// Do not decay particles that have already been decayed
// Note the herwig decayer deals with colour connections
if ( process->outgoing().empty() ) {
process->incoming()[0].first = incoming;
DecayProcessMap decay;
// Decay the particle, returning a pointer to the decay mode
decayMode = ShowerHandler::currentHandler()->decay(process,decay,true);
}
// Sort out the colour connections of particles already decayed
else {
// sort out the colour of the incoming
map<tColinePtr,tColinePtr> cmap;
if(incoming->colourLine())
cmap[process->incoming()[0].first->colourLine()] = incoming->colourLine();
if(incoming->antiColourLine())
cmap[process->incoming()[0].first->antiColourLine()] = incoming->antiColourLine();
// fix colours of outgoing
for(auto const & outg : process->outgoing()) {
map<tColinePtr,tColinePtr>::iterator it =
cmap.find(outg.first->colourLine());
if(it!=cmap.end()) {
ColinePtr c1=outg.first->colourLine();
c1->removeColoured(outg.first);
it->second->addColoured(outg.first);
}
it = cmap.find(outg.first->antiColourLine());
if(it!=cmap.end()) {
ColinePtr c1=outg.first->antiColourLine();
c1->removeAntiColoured(outg.first);
it->second->addAntiColoured(outg.first);
}
}
// swap the incoming
process->incoming()[0].first = incoming;
}
// Set the scale of all particles involved in the decay process to the
// mass of the decaying particle
// Initialise the scale for the evolution of
// the parton shower following the decay
Energy showerScale = ZERO;
// Set the scale for the evolution of the shower
showerScale = process->incoming()[0].first->momentum().m();
Energy2 decayScaleSqr = sqr( showerScale );
process->incoming()[0].first->scale( decayScaleSqr );
for(auto & outg : process->outgoing()) {
outg.first->scale( decayScaleSqr );
}
// Update the decaying particle in the process and the event
PList::iterator posOut = find(outgoing().begin(), outgoing().end(), incoming);
PList::iterator posHard = find(hard().begin(), hard().end(), incoming);
assert((posOut!=outgoing().end() && posHard==hard().end()) ||
(posOut==outgoing().end() && posHard!=hard().end()) );
if ( posOut!=outgoing().end() ) {
outgoing().erase(posOut);
}
else {
hard().erase(posHard);
}
intermediates().push_back(process->incoming()[0].first);
// Populate the children of the incoming
for(auto const & outg : process->outgoing()) {
PPtr outgoing = outg.first;
process->incoming()[0].first->addChild(outgoing);
}
// If a decayed particle is not decayed above,
// e.g. a W in a 3-body top decay, find its decaymode.
if ( powhegEmission && !decayMode ) {
string tag = incoming->dataPtr()->name() + "->";
// Must use OrderedParticles for a tag search
ShowerHandler::OrderedParticles decayOut;
for(auto const & outg : process->outgoing()) {
decayOut.insert(outg.first->dataPtr());
}
// Construct the tag
for(auto const & dec : decayOut) {
if( dec!=*decayOut.begin() ) tag += ",";
tag +=dec->name();
}
tag += ";";
// Find the decay mode
decayMode = ShowerHandler::currentHandler()->findDecayMode(tag);
}
// Perform the powheg emission
if ( powhegEmission ) {
if ( decayMode ) {
HwDecayerBasePtr decayer;
decayer = dynamic_ptr_cast<HwDecayerBasePtr>(decayMode->decayer());
if ( decayer->hasPOWHEGCorrection() ) {
// Construct a real emission process and populate its
// incoming and outcoming prior to any powheg emission
RealEmissionProcessPtr born = new_ptr( RealEmissionProcess() );
born->bornIncoming().push_back( incoming );
for(auto const & outg : process->outgoing()) {
born->bornOutgoing().push_back(outg.first);
}
// Generate any powheg emission, returning 'real'
RealEmissionProcessPtr real = decayer->generateHardest( born );
// If an emission has been attempted
// (Note if the emission fails, a null ptr is returned)
if ( real ) {
showerScale = real->pT()[ShowerInteraction::QCD];
// If an emission is generated sort out the particles
if ( !real->outgoing().empty() ) {
// Update the decay process
// Note: Do not use the new incoming particle
PPtr oldEmitter;
PPtr newEmitter;
// Use the name recoiler to avoid confusion with
// the spectator in the POWHEGDecayer
// i.e. the recoiler can be coloured or non-coloured
PPtr oldRecoiler;
PPtr newRecoiler;
if ( real->emitter() == 1 ) {
oldEmitter = real->bornOutgoing()[0];
oldRecoiler = real->bornOutgoing()[1];
newEmitter = real->outgoing()[0];
newRecoiler = real->outgoing()[1];
}
else if ( real->emitter() == 2) {
oldEmitter = real->bornOutgoing()[1];
oldRecoiler = real->bornOutgoing()[0];
newEmitter = real->outgoing()[1];
newRecoiler = real->outgoing()[0];
}
PPtr emitted = real->outgoing()[ real->emitted()-1];
// Update the scales
newRecoiler->scale(oldRecoiler->scale());
newEmitter->scale(sqr(showerScale));
emitted->scale(sqr(showerScale));
// Update the colour flow of the new outgoing particles
// Note the emitted and newEmitter are already colour
// connected by the powheg emission function
emitted->incomingColour(oldEmitter, oldEmitter->id()<0);
if ( newRecoiler->coloured() )
newRecoiler->incomingColour(oldRecoiler, oldRecoiler->id()<0);
// Update the children of the outgoing
oldRecoiler->addChild( newRecoiler );
oldEmitter->addChild( newEmitter );
oldEmitter->addChild( emitted );
// Note: The particles in the pert proc outgoing and both outgoing
// vectors of the real emission proc are in the same order
for(unsigned int ix=0;ix<real->bornOutgoing().size();++ix) {
// Update the decay process
assert(process->outgoing()[ix].first == real->bornOutgoing()[ix]);
process->outgoing()[ix].first = real->outgoing()[ix];
// Add the outgoing from the born
// decay to the event intermediates
intermediates().push_back(real->bornOutgoing()[ix]);
}
// Add the emitted to the outgoing of the decay process
process->outgoing().push_back( { emitted, PerturbativeProcessPtr() } );
}
// Else, if no emission above pTmin, set particle scales
else {
for(auto & outg : process->outgoing()) {
outg.first->scale( sqr(showerScale) );
}
powhegEmission = false;
}
}
// No powheg emission occurred:
else
powhegEmission = false;
}
// No powheg emission occurred:
else
powhegEmission = false;
}
// No powheg emission occurred:
else
powhegEmission = false;
}
// Copy the outgoing from the decay
// process to the event record
for(auto const & outg : process->outgoing()) {
if ( outg.first->coloured() )
outgoing().push_back(outg.first);
else
hard().push_back(outg.first);
}
return showerScale;
}
void DipoleEventRecord::updateDecayMom( PPtr decayParent, PerturbativeProcessPtr decayProc ) {
// Only particles that have already been decayed
// should be passed to this function
assert( !(decayProc->outgoing().empty()) );
// Create a list of the children to update their momenta
PList children;
for ( auto const & outg : decayProc->outgoing() ) {
children.push_back( outg.first );
}
// Boost the children
PList::iterator beginChildren = children.begin();
PList::iterator endChildren = children.end();
ThePEG::UtilityBase::setMomentum(beginChildren, endChildren, decayParent->momentum().vect() );
}
void DipoleEventRecord::updateDecayChainMom( PPtr decayParent, PerturbativeProcessPtr decayProc ) {
// Note - this updates the momenta of the
// outgoing of the given decay process
// Update the momenta of the outgoing from this decay
updateDecayMom( decayParent, decayProc );
// Iteratively update the momenta of the rest of the decay chain
for ( auto & outg : decayProc->outgoing() ) {
// If a child has a corresponding pert proc
// then it has decay products
if ( outg.second ) {
for ( auto & dec : theDecays ) {
if(dec.second==outg.second) {
dec.first->setMomentum(outg.first->momentum());
break;
}
}
// Iteratively update any decay products
if ( !outg.second->outgoing().empty() )
updateDecayChainMom( outg.first, outg.second );
}
}
}
void DipoleEventRecord::updateDecays(PerturbativeProcessPtr decayProc, bool iterate) {
// Note - This does not update the momenta of the outgoing
// of decayProc.
// i.e. it is for use following the (non-)showering
// of a decay when the daughter momentum are correct.
// With iterate = true, this updates the rest of the decay chain.
// Loop over the outgoing from this decay
for ( auto & outg : decayProc->outgoing() ) {
if ( outg.second && !outg.second->outgoing().empty() ) {
// Outgoing particles which have already been decayed
PPtr newDecayed = outg.first;
PerturbativeProcessPtr newDecayProc = outg.second;
// Update the outgoing momenta from this decay
updateDecayMom( newDecayed, newDecayProc);
// If this decay is already in theDecays then erase it
for ( auto const & dec : theDecays ) {
if(dec.second==newDecayProc) {
theDecays.erase(dec.first);
break;
}
}
// Add to theDecays
theDecays[newDecayed] = newDecayProc;
//assert(theDecays[newDecayed]->incoming()[0].second==decayProc);
// Iteratively update theDecays from the decay chain
if ( iterate )
updateDecays( newDecayProc );
}
// Deal with any outgoing which need to be decayed
else if ( ShowerHandler::currentHandler()->decaysInShower(outg.first->id()) ) {
PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess());
newDecay->incoming().push_back({ outg.first , decayProc } );
theDecays[outg.first] = newDecay;
}
}
}
void DipoleEventRecord::debugLastEvent(ostream& os) const {
bool first = ShowerHandler::currentHandler()->firstInteraction();
os << "--- DipoleEventRecord ----------------------------------------------------------\n";
os << " the " << (first ? "hard" : "secondary") << " subprocess is:\n"
<< (*subProcess());
os << " using PDF's " << pdfs().first.pdf() << " and "
<< pdfs().second.pdf() << "\n";
os << " chains showering currently:\n";
for ( list<DipoleChain>::const_iterator chit = theChains.begin();
chit != theChains.end(); ++chit )
os << (*chit);
os << " chains which finished showering:\n";
for ( list<DipoleChain>::const_iterator chit = theDoneChains.begin();
chit != theDoneChains.end(); ++chit )
os << (*chit);
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
diff --git a/Shower/Dipole/Makefile.am b/Shower/Dipole/Makefile.am
--- a/Shower/Dipole/Makefile.am
+++ b/Shower/Dipole/Makefile.am
@@ -1,22 +1,22 @@
SUBDIRS = Base Kernels Kinematics Utility AlphaS Merging
pkglib_LTLIBRARIES = HwDipoleShower.la
-HwDipoleShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 8:0:0
+HwDipoleShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 8:1:0
HwDipoleShower_la_LIBADD = \
Base/libHwDipoleShowerBase.la \
Kernels/libHwDipoleShowerKernels.la \
Kinematics/libHwDipoleShowerKinematics.la \
Utility/libHwDipoleShowerUtility.la \
Merging/libHwDipoleShowerMerging.la
HwDipoleShower_la_SOURCES = \
DipoleShowerHandler.h DipoleShowerHandler.fh DipoleShowerHandler.cc
pkglib_LTLIBRARIES += HwKrknloEventReweight.la
HwKrknloEventReweight_la_SOURCES = \
KrkNLO/KrknloEventReweight.h KrkNLO/KrknloEventReweight.cc
HwKrknloEventReweight_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am
--- a/Shower/QTilde/Makefile.am
+++ b/Shower/QTilde/Makefile.am
@@ -1,29 +1,29 @@
SUBDIRS = Matching
pkglib_LTLIBRARIES = HwShower.la
-HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 24:0:0
+HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 25:0:0
HwShower_la_SOURCES = \
Couplings/ShowerAlphaQCD.h Couplings/ShowerAlphaQCD.cc \
Couplings/ShowerAlphaQED.h Couplings/ShowerAlphaQED.cc\
QTildeShowerHandler.h QTildeShowerHandler.fh QTildeShowerHandler.cc \
SplittingFunctions/HalfHalfOneSplitFn.h SplittingFunctions/HalfHalfOneSplitFn.cc\
SplittingFunctions/OneOneOneSplitFn.h SplittingFunctions/OneOneOneSplitFn.cc\
SplittingFunctions/ZeroZeroOneSplitFn.h SplittingFunctions/ZeroZeroOneSplitFn.cc\
SplittingFunctions/OneHalfHalfSplitFn.h SplittingFunctions/OneHalfHalfSplitFn.cc\
SplittingFunctions/HalfOneHalfSplitFn.h SplittingFunctions/HalfOneHalfSplitFn.cc\
SplittingFunctions/CMWOneOneOneSplitFn.h SplittingFunctions/CMWOneOneOneSplitFn.cc\
SplittingFunctions/CMWHalfHalfOneSplitFn.h SplittingFunctions/CMWHalfHalfOneSplitFn.cc\
Default/QTildeSudakov.cc Default/QTildeSudakov.h\
Default/QTildeModel.cc Default/QTildeModel.h\
Default/Decay_QTildeShowerKinematics1to2.cc \
Default/Decay_QTildeShowerKinematics1to2.h \
Default/IS_QTildeShowerKinematics1to2.cc Default/IS_QTildeShowerKinematics1to2.h \
Default/FS_QTildeShowerKinematics1to2.cc Default/FS_QTildeShowerKinematics1to2.h \
Default/QTildeFinder.cc Default/QTildeFinder.h\
Default/QTildeReconstructor.cc Default/QTildeReconstructor.h Default/QTildeReconstructor.tcc \
Base/KinematicsReconstructor.cc \
Base/KinematicsReconstructor.h \
Base/KinematicsReconstructor.fh \
Base/ShowerModel.cc Base/ShowerModel.h Base/ShowerModel.fh \
Base/PartnerFinder.h Base/PartnerFinder.fh Base/PartnerFinder.cc \
Base/ShowerVeto.h Base/ShowerVeto.fh Base/ShowerVeto.cc \
Base/FullShowerVeto.h Base/FullShowerVeto.fh Base/FullShowerVeto.cc
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,240 +1,242 @@
dnl Process this file with autoconf to produce a configure script.
AC_PREREQ([2.63])
AC_INIT([Herwig],[devel],[herwig@projects.hepforge.org],[Herwig])
AC_CONFIG_SRCDIR([Utilities/HerwigStrategy.cc])
AC_CONFIG_AUX_DIR([Config])
AC_CONFIG_MACRO_DIR([m4])
AC_CONFIG_HEADERS([Config/config.h])
dnl AC_PRESERVE_HELP_ORDER
AC_CANONICAL_HOST
dnl === disable debug symbols by default =====
if test "x$CXXFLAGS" = "x"; then
CXXFLAGS=-O2
fi
if test "x$CFLAGS" = "x"; then
CFLAGS=-O2
fi
AC_LANG([C++])
AM_INIT_AUTOMAKE([1.11 subdir-objects gnu dist-bzip2 no-dist-gzip -Wall -Wno-portability])
m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
m4_ifdef([AM_PROG_AR], [AM_PROG_AR])
dnl Checks for C++ compiler. Handle C++11 flags.
AC_PROG_CXX
AX_CXX_COMPILE_STDCXX([11],[noext],[mandatory])
dnl check for POSIX
AC_CHECK_HEADER([unistd.h],[],
[AC_MSG_ERROR([Herwig needs "unistd.h". Non-POSIX systems are not supported.])])
dnl Checks for programs.
AC_PROG_INSTALL
AC_PROG_MAKE_SET
AC_PROG_LN_S
dnl modified search order
AC_PROG_FC([gfortran g95 g77])
dnl xlf95 f95 fort ifort ifc efc pgf95 lf95 ftn xlf90 f90 pgf90 pghpf epcf90 xlf f77 frt pgf77 cf77 fort77 fl32 af77])
AC_LANG_PUSH([Fortran])
AC_MSG_CHECKING([if the Fortran compiler ($FC) works])
AC_COMPILE_IFELSE(
AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
[AC_MSG_RESULT([yes])],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([A Fortran compiler is required to build Herwig.])
]
)
AC_LANG_POP([Fortran])
LT_PREREQ([2.2.6])
LT_INIT([disable-static dlopen pic-only])
dnl ####################################
dnl ####################################
dnl for Doc/fixinterfaces.pl
AC_PATH_PROG(PERL, perl)
dnl for Models/Feynrules
AM_PATH_PYTHON([2.6],, [:])
AM_CONDITIONAL([HAVE_PYTHON], [test "x$PYTHON" != "x:"])
HERWIG_CHECK_GSL
HERWIG_CHECK_THEPEG
BOOST_REQUIRE([1.41])
BOOST_FIND_HEADER([boost/numeric/ublas/io.hpp])
-BOOST_FIND_HEADER([boost/numeric/ublas/matrix.hpp])
-BOOST_FIND_HEADER([boost/numeric/ublas/matrix_proxy.hpp])
-BOOST_FIND_HEADER([boost/numeric/ublas/matrix_sparse.hpp])
-BOOST_FIND_HEADER([boost/numeric/ublas/symmetric.hpp])
-BOOST_FIND_HEADER([boost/numeric/ublas/vector.hpp])
+dnl Boost 1.64 is missing a required header to make these work
+dnl we just assume they're there if io.hpp has been found OK above
+dnl BOOST_FIND_HEADER([boost/numeric/ublas/matrix.hpp])
+dnl BOOST_FIND_HEADER([boost/numeric/ublas/matrix_proxy.hpp])
+dnl BOOST_FIND_HEADER([boost/numeric/ublas/matrix_sparse.hpp])
+dnl BOOST_FIND_HEADER([boost/numeric/ublas/symmetric.hpp])
+dnl BOOST_FIND_HEADER([boost/numeric/ublas/vector.hpp])
BOOST_FIND_HEADER([boost/operators.hpp])
BOOST_FIND_HEADER([boost/progress.hpp])
BOOST_FILESYSTEM([mt])
BOOST_TEST()
HERWIG_CHECK_VBFNLO
HERWIG_CHECK_NJET
HERWIG_CHECK_GOSAM
HERWIG_CHECK_GOSAM_CONTRIB
HERWIG_CHECK_OPENLOOPS
HERWIG_CHECK_MADGRAPH
HERWIG_CHECK_EVTGEN
HERWIG_CHECK_PYTHIA
HERWIG_COMPILERFLAGS
HERWIG_LOOPTOOLS
FASTJET_CHECK_FASTJET
HERWIG_ENABLE_MODELS
SHARED_FLAG=-shared
AM_CONDITIONAL(NEED_APPLE_FIXES,
[test "xx${host/darwin/foundit}xx" != "xx${host}xx"])
if test "xx${host/darwin/foundit}xx" != "xx${host}xx"; then
APPLE_DSO_FLAGS=-Wl,-undefined,dynamic_lookup
SHARED_FLAG=-bundle
fi
AC_SUBST([APPLE_DSO_FLAGS])
AC_SUBST([SHARED_FLAG])
AC_CONFIG_FILES([UnderlyingEvent/Makefile
Models/Makefile
Models/StandardModel/Makefile
Models/RSModel/Makefile
Models/General/Makefile
Models/Susy/Makefile
Models/Susy/NMSSM/Makefile
Models/Susy/RPV/Makefile
Models/UED/Makefile
Models/LH/Makefile
Models/LHTP/Makefile
Models/Transplanckian/Makefile
Models/Leptoquarks/Makefile
Models/Zprime/Makefile
Models/TTbAsymm/Makefile
Models/Feynrules/Makefile
Models/Feynrules/python/Makefile-FR
Models/ADD/Makefile
Models/Sextet/Makefile
Decay/Makefile
Decay/FormFactors/Makefile
Decay/Tau/Makefile
Decay/Baryon/Makefile
Decay/VectorMeson/Makefile
Decay/Perturbative/Makefile
Decay/ScalarMeson/Makefile
Decay/TensorMeson/Makefile
Decay/WeakCurrents/Makefile
Decay/Partonic/Makefile
Decay/General/Makefile
Decay/Radiation/Makefile
Decay/EvtGen/Makefile
Doc/refman.conf
Doc/refman.h
PDT/Makefile
PDF/Makefile
MatrixElement/Makefile
MatrixElement/General/Makefile
MatrixElement/Lepton/Makefile
MatrixElement/Hadron/Makefile
MatrixElement/DIS/Makefile
MatrixElement/Powheg/Makefile
MatrixElement/Gamma/Makefile
MatrixElement/Reweighters/Makefile
MatrixElement/Matchbox/Makefile
MatrixElement/Matchbox/Base/Makefile
MatrixElement/Matchbox/Utility/Makefile
MatrixElement/Matchbox/Phasespace/Makefile
MatrixElement/Matchbox/Dipoles/Makefile
MatrixElement/Matchbox/InsertionOperators/Makefile
MatrixElement/Matchbox/Matching/Makefile
MatrixElement/Matchbox/Cuts/Makefile
MatrixElement/Matchbox/Scales/Makefile
MatrixElement/Matchbox/ColorFull/Makefile
MatrixElement/Matchbox/CVolver/Makefile
MatrixElement/Matchbox/Builtin/Makefile
MatrixElement/Matchbox/Builtin/Amplitudes/Makefile
MatrixElement/Matchbox/Tests/Makefile
MatrixElement/Matchbox/External/Makefile
MatrixElement/Matchbox/External/BLHAGeneric/Makefile
MatrixElement/Matchbox/External/VBFNLO/Makefile
MatrixElement/Matchbox/External/NJet/Makefile
MatrixElement/Matchbox/External/GoSam/Makefile
MatrixElement/Matchbox/External/OpenLoops/Makefile
MatrixElement/Matchbox/External/MadGraph/Makefile
MatrixElement/Matchbox/External/MadGraph/mg2herwig.py
Sampling/Makefile
Sampling/CellGrids/Makefile
Shower/Makefile
Shower/QTilde/Makefile
Shower/QTilde/Matching/Makefile
Shower/Dipole/Makefile
Shower/Dipole/Base/Makefile
Shower/Dipole/Kernels/Makefile
Shower/Dipole/Kinematics/Makefile
Shower/Dipole/Utility/Makefile
Shower/Dipole/AlphaS/Makefile
Utilities/Makefile
Utilities/XML/Makefile
Utilities/Statistics/Makefile
Hadronization/Makefile
lib/Makefile
include/Makefile
src/Makefile
src/defaults/Makefile
src/snippets/Makefile
src/Matchbox/Makefile
src/herwig-config
Doc/Makefile
Doc/HerwigDefaults.in
Looptools/Makefile
Analysis/Makefile
API/Makefile
src/Makefile-UserModules
src/defaults/Analysis.in
src/defaults/MatchboxDefaults.in
src/defaults/Decays.in
src/defaults/decayers.in
src/defaults/setup.gosam.in
src/Matchbox/LO-DefaultShower.in
src/Matchbox/LO-DipoleShower.in
src/Matchbox/MCatLO-DefaultShower.in
src/Matchbox/MCatLO-DipoleShower.in
src/Matchbox/LO-NoShower.in
src/Matchbox/MCatNLO-DefaultShower.in
src/Matchbox/MCatNLO-DipoleShower.in
src/Matchbox/NLO-NoShower.in
src/Matchbox/Powheg-DefaultShower.in
src/Matchbox/Powheg-DipoleShower.in
src/Merging/Makefile
Shower/Dipole/Merging/Makefile
src/defaults/MatchboxMergingDefaults.in
Contrib/Makefile
Contrib/make_makefiles.sh
Tests/Makefile
Makefile])
AC_CONFIG_LINKS([Doc/BSMlibs.in:Doc/BSMlibs.in])
AC_CONFIG_FILES([Doc/fixinterfaces.pl],[chmod +x Doc/fixinterfaces.pl])
HERWIG_OVERVIEW
AC_CONFIG_COMMANDS([summary],[cat config.herwig])
AC_OUTPUT
diff --git a/lib/Makefile.am b/lib/Makefile.am
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -1,49 +1,49 @@
pkglib_LTLIBRARIES = Herwig.la
Herwig_la_SOURCES =
Herwig_la_LIBTOOLFLAGS = --tag=CXX
-Herwig_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 21:0:0
+Herwig_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 22:0:0
Herwig_la_LDFLAGS += $(THEPEGLDFLAGS) $(BOOST_SYSTEM_LDFLAGS) $(BOOST_FILESYSTEM_LDFLAGS) $(FCLIBS)
Herwig_la_LIBADD = \
$(top_builddir)/Hadronization/libHwHadronization.la \
$(top_builddir)/Models/libHwStandardModel.la \
$(top_builddir)/Decay/libHwDecay.la \
$(top_builddir)/Decay/libHwFormFactor.la \
$(top_builddir)/Decay/libHwDecRad.la \
$(top_builddir)/Utilities/libHwUtils.la \
$(top_builddir)/Models/libHwModelGenerator.la \
$(top_builddir)/Decay/General/libHwGeneralDecay.la \
$(top_builddir)/MatrixElement/General/libHwGeneralME.la \
$(top_builddir)/MatrixElement/libHwME.la \
$(top_builddir)/MatrixElement/Reweighters/libHwReweighters.la \
$(top_builddir)/MatrixElement/Matchbox/libHwMatchbox.la \
$(top_builddir)/Decay/libHwWeakCurrent.la \
$(top_builddir)/Looptools/libHwLooptools.la \
$(top_builddir)/Shower/libHwShower.la \
$(THEPEGLIB) $(BOOST_SYSTEM_LIBS) $(BOOST_FILESYSTEM_LIBS) -ldl
dist_noinst_SCRIPTS = fix-osx-path
POSTPROCESSING = done-all-links
if NEED_APPLE_FIXES
POSTPROCESSING += apple-fixes
endif
all-local: $(POSTPROCESSING)
done-all-links: Herwig.la
find $(top_builddir) \( -name '*.so.*' -or -name '*.so' \) \
-not -name 'lib*' -not -path '$(top_builddir)/lib/*' \
-not -path '$(top_builddir)/.hg/*' -exec $(LN_S) -f \{\} \;
$(LN_S) -f .libs/Herwig*so* .
echo "stamp" > $@
apple-fixes: fix-osx-path done-all-links
./$<
echo "stamp" > $@
clean-local:
rm -f *.so *.so.* done-all-links apple-fixes
diff --git a/src/DIS-Matchbox.in b/src/DIS-Matchbox.in
--- a/src/DIS-Matchbox.in
+++ b/src/DIS-Matchbox.in
@@ -1,143 +1,150 @@
# -*- ThePEG-repository -*-
##################################################
## Herwig/Matchbox example input file
##################################################
##################################################
## Collider type
##################################################
read snippets/Matchbox.in
read snippets/EPCollider.in
##################################################
## Beam energies
##################################################
cd /Herwig/EventHandlers
##################################################
## Process selection
##################################################
## Note that event generation may fail if no matching matrix element has
## been found. Coupling orders are with respect to the Born process,
## i.e. NLO QCD does not require an additional power of alphas.
## Model assumptions
read Matchbox/StandardModelLike.in
read Matchbox/DiagonalCKM.in
+## EP collider parameters
+cd /Herwig/EventHandlers
+set EventHandler:BeamA /Herwig/Particles/e-
+set Luminosity:BeamEMaxA 30.*GeV
+set EventHandler:BeamB /Herwig/Particles/p+
+set Luminosity:BeamEMaxB 920.*GeV
+
## Set the order of the couplings
cd /Herwig/MatrixElements/Matchbox
set Factory:OrderInAlphaS 0
set Factory:OrderInAlphaEW 2
## Select the process
## You may use identifiers such as p, pbar, j, l, mu+, h0 etc.
do Factory:Process e- p -> e- j
## Special settings required for on-shell production of unstable particles
## enable for on-shell top production
# read Matchbox/OnShellTopProduction.in
## enable for on-shell W, Z or h production
# read Matchbox/OnShellWProduction.in
# read Matchbox/OnShellZProduction.in
# read Matchbox/OnShellHProduction.in
##################################################
## Matrix element library selection
##################################################
## Select a generic tree/loop combination or a
## specialized NLO package
# read Matchbox/MadGraph-GoSam.in
# read Matchbox/MadGraph-MadGraph.in
# read Matchbox/MadGraph-NJet.in
# read Matchbox/MadGraph-OpenLoops.in
##################################################
## Cut selection
## See the documentation for more options
##################################################
## cuts on additional jets
# read Matchbox/DefaultEPJets.in
# insert JetCuts:JetRegions 0 FirstJet
# insert JetCuts:JetRegions 1 SecondJet
# insert JetCuts:JetRegions 2 ThirdJet
##################################################
## Scale choice
## See the documentation for more options
##################################################
cd /Herwig/MatrixElements/Matchbox
set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonQ2Scale
##################################################
## Matching and shower selection
## Please also see flavour scheme settings
## towards the end of the input file.
##################################################
read Matchbox/MCatNLO-DefaultShower.in
# read Matchbox/Powheg-DefaultShower.in
## use for strict LO/NLO comparisons
# read Matchbox/MCatLO-DefaultShower.in
## use for improved LO showering
# read Matchbox/LO-DefaultShower.in
# read Matchbox/MCatNLO-DipoleShower.in
# read Matchbox/Powheg-DipoleShower.in
## use for strict LO/NLO comparisons
# read Matchbox/MCatLO-DipoleShower.in
## use for improved LO showering
# read Matchbox/LO-DipoleShower.in
# read Matchbox/NLO-NoShower.in
# read Matchbox/LO-NoShower.in
##################################################
## Scale uncertainties
##################################################
# read Matchbox/MuDown.in
# read Matchbox/MuUp.in
##################################################
## Shower scale uncertainties
##################################################
# read Matchbox/MuQDown.in
# read Matchbox/MuQUp.in
##################################################
## PDF choice
##################################################
read Matchbox/FiveFlavourScheme.in
## required for dipole shower and fixed order in five flavour scheme
# read Matchbox/FiveFlavourNoBMassScheme.in
read Matchbox/MMHT2014.in
##################################################
## Analyses
##################################################
# cd /Herwig/Analysis
# insert Rivet:Analyses 0 XXX_2015_ABC123
# insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 Rivet
# insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 HepMC
##################################################
## Save the generator
##################################################
do /Herwig/MatrixElements/Matchbox/Factory:ProductionMode
cd /Herwig/Generators
saverun DIS-Matchbox EventGenerator
diff --git a/src/defaults/MatchboxDefaults.in.in b/src/defaults/MatchboxDefaults.in.in
--- a/src/defaults/MatchboxDefaults.in.in
+++ b/src/defaults/MatchboxDefaults.in.in
@@ -1,707 +1,709 @@
# -*- ThePEG-repository -*-
################################################################################
#
# Default setup for Matchbox matrix element generation.
# You do not need to make any change in here; processes of
# interest can be chosen in the standard input files.
#
################################################################################
################################################################################
# Load libraries
################################################################################
library JetCuts.so
library FastJetFinder.so
library HwMatchboxScales.so
library HwMatchboxCuts.so
library HwColorFull.so
library HwMatchboxBuiltin.so
################################################################################
# Setup the factory object
################################################################################
mkdir /Herwig/MatrixElements/Matchbox
cd /Herwig/MatrixElements/Matchbox
create Herwig::MatchboxFactory Factory
do Factory:StartParticleGroup p
insert Factory:ParticleGroup 0 /Herwig/Particles/b
insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
insert Factory:ParticleGroup 0 /Herwig/Particles/c
insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
insert Factory:ParticleGroup 0 /Herwig/Particles/s
insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
insert Factory:ParticleGroup 0 /Herwig/Particles/d
insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
insert Factory:ParticleGroup 0 /Herwig/Particles/u
insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
insert Factory:ParticleGroup 0 /Herwig/Particles/g
do Factory:EndParticleGroup
do Factory:StartParticleGroup pbar
insert Factory:ParticleGroup 0 /Herwig/Particles/b
insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
insert Factory:ParticleGroup 0 /Herwig/Particles/c
insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
insert Factory:ParticleGroup 0 /Herwig/Particles/s
insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
insert Factory:ParticleGroup 0 /Herwig/Particles/d
insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
insert Factory:ParticleGroup 0 /Herwig/Particles/u
insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
insert Factory:ParticleGroup 0 /Herwig/Particles/g
do Factory:EndParticleGroup
do Factory:StartParticleGroup j
insert Factory:ParticleGroup 0 /Herwig/Particles/b
insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
insert Factory:ParticleGroup 0 /Herwig/Particles/c
insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
insert Factory:ParticleGroup 0 /Herwig/Particles/s
insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
insert Factory:ParticleGroup 0 /Herwig/Particles/d
insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
insert Factory:ParticleGroup 0 /Herwig/Particles/u
insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
insert Factory:ParticleGroup 0 /Herwig/Particles/g
do Factory:EndParticleGroup
do Factory:StartParticleGroup u
insert Factory:ParticleGroup 0 /Herwig/Particles/u
do Factory:EndParticleGroup
do Factory:StartParticleGroup ubar
insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
do Factory:EndParticleGroup
do Factory:StartParticleGroup d
insert Factory:ParticleGroup 0 /Herwig/Particles/d
do Factory:EndParticleGroup
do Factory:StartParticleGroup dbar
insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
do Factory:EndParticleGroup
do Factory:StartParticleGroup s
insert Factory:ParticleGroup 0 /Herwig/Particles/s
do Factory:EndParticleGroup
do Factory:StartParticleGroup sbar
insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
do Factory:EndParticleGroup
do Factory:StartParticleGroup c
insert Factory:ParticleGroup 0 /Herwig/Particles/c
do Factory:EndParticleGroup
do Factory:StartParticleGroup cbar
insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
do Factory:EndParticleGroup
do Factory:StartParticleGroup b
insert Factory:ParticleGroup 0 /Herwig/Particles/b
do Factory:EndParticleGroup
do Factory:StartParticleGroup bbar
insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
do Factory:EndParticleGroup
do Factory:StartParticleGroup t
insert Factory:ParticleGroup 0 /Herwig/Particles/t
do Factory:EndParticleGroup
do Factory:StartParticleGroup tbar
insert Factory:ParticleGroup 0 /Herwig/Particles/tbar
do Factory:EndParticleGroup
do Factory:StartParticleGroup g
insert Factory:ParticleGroup 0 /Herwig/Particles/g
do Factory:EndParticleGroup
do Factory:StartParticleGroup gamma
insert Factory:ParticleGroup 0 /Herwig/Particles/gamma
do Factory:EndParticleGroup
do Factory:StartParticleGroup h0
insert Factory:ParticleGroup 0 /Herwig/Particles/h0
do Factory:EndParticleGroup
do Factory:StartParticleGroup W+
insert Factory:ParticleGroup 0 /Herwig/Particles/W+
do Factory:EndParticleGroup
do Factory:StartParticleGroup W-
insert Factory:ParticleGroup 0 /Herwig/Particles/W-
do Factory:EndParticleGroup
do Factory:StartParticleGroup Z0
insert Factory:ParticleGroup 0 /Herwig/Particles/Z0
do Factory:EndParticleGroup
do Factory:StartParticleGroup e+
insert Factory:ParticleGroup 0 /Herwig/Particles/e+
do Factory:EndParticleGroup
do Factory:StartParticleGroup e-
insert Factory:ParticleGroup 0 /Herwig/Particles/e-
do Factory:EndParticleGroup
do Factory:StartParticleGroup mu+
insert Factory:ParticleGroup 0 /Herwig/Particles/mu+
do Factory:EndParticleGroup
do Factory:StartParticleGroup mu-
insert Factory:ParticleGroup 0 /Herwig/Particles/mu-
do Factory:EndParticleGroup
do Factory:StartParticleGroup tau+
insert Factory:ParticleGroup 0 /Herwig/Particles/tau+
do Factory:EndParticleGroup
do Factory:StartParticleGroup tau-
insert Factory:ParticleGroup 0 /Herwig/Particles/tau-
do Factory:EndParticleGroup
do Factory:StartParticleGroup nu_e
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_e
do Factory:EndParticleGroup
do Factory:StartParticleGroup nu_mu
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mu
do Factory:EndParticleGroup
do Factory:StartParticleGroup nu_tau
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_tau
do Factory:EndParticleGroup
do Factory:StartParticleGroup nu_ebar
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_ebar
do Factory:EndParticleGroup
do Factory:StartParticleGroup nu_mubar
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mubar
do Factory:EndParticleGroup
do Factory:StartParticleGroup nu_taubar
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_taubar
do Factory:EndParticleGroup
do Factory:StartParticleGroup l
insert Factory:ParticleGroup 0 /Herwig/Particles/e+
insert Factory:ParticleGroup 0 /Herwig/Particles/mu+
insert Factory:ParticleGroup 0 /Herwig/Particles/e-
insert Factory:ParticleGroup 0 /Herwig/Particles/mu-
do Factory:EndParticleGroup
do Factory:StartParticleGroup nu
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_e
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mu
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_ebar
insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mubar
do Factory:EndParticleGroup
do Factory:StartParticleGroup l+
insert Factory:ParticleGroup 0 /Herwig/Particles/e+
insert Factory:ParticleGroup 0 /Herwig/Particles/mu+
do Factory:EndParticleGroup
do Factory:StartParticleGroup l-
insert Factory:ParticleGroup 0 /Herwig/Particles/e-
insert Factory:ParticleGroup 0 /Herwig/Particles/mu-
do Factory:EndParticleGroup
################################################################################
# Default settings for hard process widths
################################################################################
set /Herwig/Particles/mu+:HardProcessWidth 0*GeV
set /Herwig/Particles/mu-:HardProcessWidth 0*GeV
set /Herwig/Particles/tau+:HardProcessWidth 0*GeV
set /Herwig/Particles/tau-:HardProcessWidth 0*GeV
################################################################################
# Setup amplitudes
################################################################################
cd /Herwig/MatrixElements/Matchbox
mkdir Amplitudes
cd Amplitudes
create ColorFull::TraceBasis TraceBasis
create Herwig::MatchboxHybridAmplitude GenericProcesses
@LOAD_MADGRAPH@ HwMatchboxMadGraph.so
@CREATE_MADGRAPH@ Herwig::MadGraphAmplitude MadGraph
@SET_MADGRAPH@ MadGraph:ColourBasis TraceBasis
@LOAD_GOSAM@ HwMatchboxGoSam.so
@CREATE_GOSAM@ Herwig::GoSamAmplitude GoSam
@LOAD_NJET@ HwMatchboxNJet.so
@CREATE_NJET@ Herwig::NJetsAmplitude NJet
@DO_NJET@ NJet:Massless 5
@DO_NJET@ NJet:Massless -5
@LOAD_OPENLOOPS@ HwMatchboxOpenLoops.so
@CREATE_OPENLOOPS@ Herwig::OpenLoopsAmplitude OpenLoops
@LOAD_VBFNLO@ HwMatchboxVBFNLO.so
@CREATE_VBFNLO@ Herwig::VBFNLOAmplitude VBFNLO
mkdir Builtin
cd Builtin
create Herwig::SimpleColourBasis SimpleColourBasis
create Herwig::SimpleColourBasis2 SimpleColourBasis2
create Herwig::MatchboxAmplitudellbarqqbar Amplitudellbarqqbar
set Amplitudellbarqqbar:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudellbarqqbarg Amplitudellbarqqbarg
set Amplitudellbarqqbarg:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudellbarqqbargg Amplitudellbarqqbargg
set Amplitudellbarqqbargg:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudellbarqqbarqqbar Amplitudellbarqqbarqqbar
set Amplitudellbarqqbarqqbar:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudelnuqqbar Amplitudelnuqqbar
set Amplitudelnuqqbar:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudelnuqqbarg Amplitudelnuqqbarg
set Amplitudelnuqqbarg:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudelnuqqbargg Amplitudelnuqqbargg
set Amplitudelnuqqbargg:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudelnuqqbarqqbar Amplitudelnuqqbarqqbar
set Amplitudelnuqqbarqqbar:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudehgg Amplitudehgg
set Amplitudehgg:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudehggg Amplitudehggg
set Amplitudehggg:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudehqqbarg Amplitudehqqbarg
set Amplitudehqqbarg:ColourBasis SimpleColourBasis
create Herwig::MatchboxAmplitudeqqbarttbar Amplitudeqqbarttbar
set Amplitudeqqbarttbar:ColourBasis SimpleColourBasis2
create Herwig::MatchboxAmplitudeqqbarttbarg Amplitudeqqbarttbarg
set Amplitudeqqbarttbarg:ColourBasis SimpleColourBasis2
create Herwig::MatchboxAmplitudeggttbar Amplitudeggttbar
set Amplitudeggttbar:ColourBasis SimpleColourBasis2
create Herwig::MatchboxAmplitudeggttbarg Amplitudeggttbarg
set Amplitudeggttbarg:ColourBasis SimpleColourBasis2
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbar
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbarg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbargg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbarqqbar
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbar
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbarg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbargg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbarqqbar
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehgg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehggg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehqqbarg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeqqbarttbar
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeqqbarttbarg
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeggttbar
insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeggttbarg
################################################################################
# Setup phasespace generators
################################################################################
cd /Herwig/MatrixElements/Matchbox
mkdir Phasespace
cd Phasespace
create Herwig::PhasespaceCouplings PhasespaceCouplings
create Herwig::MatchboxRambo Rambo
set Rambo:CouplingData PhasespaceCouplings
create Herwig::FlatInvertiblePhasespace InvertiblePhasespace
set InvertiblePhasespace:CouplingData PhasespaceCouplings
create Herwig::FlatInvertibleLabframePhasespace InvertibleLabframePhasespace
set InvertibleLabframePhasespace:CouplingData PhasespaceCouplings
set InvertibleLabframePhasespace:LogSHat No
create Herwig::TreePhasespaceChannels TreePhasespaceChannels
create Herwig::TreePhasespace TreePhasespace
set TreePhasespace:ChannelMap TreePhasespaceChannels
set TreePhasespace:M0 0.0001*GeV
set TreePhasespace:MC 0.00005*GeV
set TreePhasespace:CouplingData PhasespaceCouplings
do TreePhasespace:SetPhysicalCoupling 21 -1 1 0.059
do TreePhasespace:SetPhysicalCoupling 21 -2 2 0.059
do TreePhasespace:SetPhysicalCoupling 21 -3 3 0.059
do TreePhasespace:SetPhysicalCoupling 21 -4 4 0.059
do TreePhasespace:SetPhysicalCoupling 21 -5 5 0.059
do TreePhasespace:SetPhysicalCoupling 21 -6 6 0.059
do TreePhasespace:SetPhysicalCoupling 21 1 -1 0.059
do TreePhasespace:SetPhysicalCoupling 21 2 -2 0.059
do TreePhasespace:SetPhysicalCoupling 21 3 -3 0.059
do TreePhasespace:SetPhysicalCoupling 21 4 -4 0.059
do TreePhasespace:SetPhysicalCoupling 21 5 -5 0.059
do TreePhasespace:SetPhysicalCoupling 21 6 -6 0.059
do TreePhasespace:SetPhysicalCoupling 1 21 1 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 2 21 2 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 3 21 3 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 4 21 4 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 5 21 5 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 6 21 6 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -1 21 -1 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -2 21 -2 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -3 21 -3 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -4 21 -4 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -5 21 -5 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -6 21 -6 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 1 1 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 2 2 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 3 3 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 4 4 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 5 5 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling 6 6 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -1 -1 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -2 -2 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -3 -3 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -4 -4 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -5 -5 21 0.15733333333333333333
do TreePhasespace:SetPhysicalCoupling -6 -6 21 0.15733333333333333333
do TreePhasespace:SetCoupling 25 -1 1 0
do TreePhasespace:SetCoupling 25 -2 2 0
do TreePhasespace:SetCoupling 25 -3 3 0.00000001184279069851
do TreePhasespace:SetCoupling 25 -4 4 0.00000205034465001885
do TreePhasespace:SetCoupling 25 -5 5 0.00002314757096085280
do TreePhasespace:SetCoupling 25 -6 6 0.03982017320025470767
do TreePhasespace:SetCoupling 25 -11 11 0.00000000000034264835
do TreePhasespace:SetCoupling 25 -12 12 0
do TreePhasespace:SetCoupling 25 -13 13 0.00000001464912263400
do TreePhasespace:SetCoupling 25 -14 14 0
do TreePhasespace:SetCoupling 25 -15 15 0.00000414359033108195
do TreePhasespace:SetCoupling 25 -16 16 0
do TreePhasespace:SetCoupling 22 -1 1 0.00083932358497608365
do TreePhasespace:SetCoupling 22 -2 2 0.00335729433990433461
do TreePhasespace:SetCoupling 22 -3 3 0.00083932358497608365
do TreePhasespace:SetCoupling 22 -4 4 0.00335729433990433461
do TreePhasespace:SetCoupling 22 -5 5 0.00083932358497608365
do TreePhasespace:SetCoupling 22 -6 6 0.00335729433990433461
do TreePhasespace:SetCoupling 22 -11 11 0.00755391226478475287
do TreePhasespace:SetCoupling 22 -13 13 0.00755391226478475287
do TreePhasespace:SetCoupling 22 -15 15 0.00755391226478475287
do TreePhasespace:SetCoupling 24 -2 1 0.01652748072644379386
do TreePhasespace:SetCoupling 24 -4 1 0.00382028458188709739
do TreePhasespace:SetCoupling 24 -6 1 0.00014707756360995175
do TreePhasespace:SetCoupling 24 -2 3 0.00382265953677814621
do TreePhasespace:SetCoupling 24 -4 3 0.01651340063673257587
do TreePhasespace:SetCoupling 24 -6 3 0.00068534412570265868
do TreePhasespace:SetCoupling 24 -2 5 0.00005954351191129535
do TreePhasespace:SetCoupling 24 -4 5 0.00069891529650865192
do TreePhasespace:SetCoupling 24 -6 5 0.01694947628265615369
do TreePhasespace:SetCoupling 24 -12 11 0.01696396350749155147
do TreePhasespace:SetCoupling 24 -14 13 0.01696396350749155147
do TreePhasespace:SetCoupling 24 -16 15 0.01696396350749155147
do TreePhasespace:SetCoupling -24 2 -1 0.01652748072644379386
do TreePhasespace:SetCoupling -24 4 -1 0.00382028458188709739
do TreePhasespace:SetCoupling -24 6 -1 0.00014707756360995175
do TreePhasespace:SetCoupling -24 2 -3 0.00382265953677814621
do TreePhasespace:SetCoupling -24 4 -3 0.01651340063673257587
do TreePhasespace:SetCoupling -24 6 -3 0.00068534412570265868
do TreePhasespace:SetCoupling -24 2 -5 0.00005954351191129535
do TreePhasespace:SetCoupling -24 4 -5 0.00069891529650865192
do TreePhasespace:SetCoupling -24 6 -5 0.01694947628265615369
do TreePhasespace:SetCoupling -24 12 -11 0.01696396350749155147
do TreePhasespace:SetCoupling -24 14 -13 0.01696396350749155147
do TreePhasespace:SetCoupling -24 16 -15 0.01696396350749155147
do TreePhasespace:SetCoupling 23 -1 1 0.00407649129960709158
do TreePhasespace:SetCoupling 23 -2 2 0.00317809816318353030
do TreePhasespace:SetCoupling 23 -3 3 0.00407649129960709158
do TreePhasespace:SetCoupling 23 -4 4 0.00317809816318353030
do TreePhasespace:SetCoupling 23 -5 5 0.00407649129960709158
do TreePhasespace:SetCoupling 23 -6 6 0.00317809816318353030
do TreePhasespace:SetCoupling 23 -11 11 0.00276049468148072129
do TreePhasespace:SetCoupling 23 -12 12 0.00545567409075140513
do TreePhasespace:SetCoupling 23 -13 13 0.00276049468148072129
do TreePhasespace:SetCoupling 23 -14 14 0.00545567409075140513
do TreePhasespace:SetCoupling 23 -15 15 0.00276049468148072129
do TreePhasespace:SetCoupling 23 -16 16 0.00545567409075140513
do TreePhasespace:SetCoupling 21 21 21 0.354
do TreePhasespace:SetCoupling 25 21 21 0.00000000016160437564
do TreePhasespace:SetCoupling 25 25 25 0.18719783125611995353
do TreePhasespace:SetCoupling 25 22 22 0.00000000006295673620
do TreePhasespace:SetCoupling 25 24 -24 219.30463760755686425818
do TreePhasespace:SetCoupling 25 23 23 362.91922658249853887524
do TreePhasespace:SetCoupling 22 24 -24 0.00755391226478475287
do TreePhasespace:SetCoupling 23 24 -24 0.02637401475019835008
@CREATE_VBFNLO@ Herwig::VBFNLOPhasespace VBFNLOPhasespace
@SET_VBFNLO@ VBFNLOPhasespace:CouplingData PhasespaceCouplings
set /Herwig/MatrixElements/Matchbox/Factory:Phasespace TreePhasespace
################################################################################
# Setup utilities for matching
################################################################################
cd /Herwig/MatrixElements/Matchbox
create Herwig::HardScaleProfile HardScaleProfile
create Herwig::MEMatching MEMatching
set MEMatching:RestrictPhasespace Yes
set MEMatching:HardScaleProfile /Herwig/MatrixElements/Matchbox/HardScaleProfile
set MEMatching:BornScaleInSubtraction BornScale
set MEMatching:RealEmissionScaleInSubtraction RealScale
set MEMatching:EmissionScaleInSubtraction RealScale
set MEMatching:BornScaleInSplitting ShowerScale
set MEMatching:RealEmissionScaleInSplitting ShowerScale
set MEMatching:EmissionScaleInSplitting ShowerScale
set MEMatching:TruncatedShower Yes
set MEMatching:MaxPtIsMuF Yes
set MEMatching:FFPtCut 1.0*GeV
set MEMatching:FIPtCut 1.0*GeV
set MEMatching:IIPtCut 1.0*GeV
set MEMatching:SafeCut 0.*GeV
create Herwig::ShowerApproximationGenerator MECorrectionHandler
set MECorrectionHandler:ShowerApproximation MEMatching
set MECorrectionHandler:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/InvertiblePhasespace
set MECorrectionHandler:PresamplingPoints 50000
set MECorrectionHandler:FreezeGrid 100000
create Herwig::DipoleMatching DipoleMatching HwDipoleMatching.so
# set in DipoleShowerDefaults.in as not available at this point
# set DipoleMatching:ShowerHandler /Herwig/DipoleShower/DipoleShowerHandler
set DipoleMatching:BornScaleInSubtraction BornScale
set DipoleMatching:RealEmissionScaleInSubtraction BornScale
set DipoleMatching:EmissionScaleInSubtraction BornScale
set DipoleMatching:FFPtCut 1.0*GeV
set DipoleMatching:FIPtCut 1.0*GeV
set DipoleMatching:IIPtCut 1.0*GeV
set DipoleMatching:SafeCut 4.*GeV
create Herwig::QTildeMatching QTildeMatching HwQTildeMatching.so
set QTildeMatching:ShowerHandler /Herwig/Shower/ShowerHandler
set QTildeMatching:BornScaleInSubtraction BornScale
set QTildeMatching:RealEmissionScaleInSubtraction BornScale
set QTildeMatching:EmissionScaleInSubtraction BornScale
set QTildeMatching:QTildeFinder /Herwig/Shower/PartnerFinder
set QTildeMatching:SafeCut 4.*GeV
# just a dummy, since SudakovCommonn can't be used
# it's only used to get the value of the kinCutoffScale
set QTildeMatching:QTildeSudakov /Herwig/Shower/QtoQGSudakov
################################################################################
# Setup utilities for process generation
################################################################################
cd /Herwig/MatrixElements/Matchbox
mkdir Utility
cd Utility
create Herwig::Tree2toNGenerator DiagramGenerator
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFGVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/GGGVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFPVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFZVertex
cp /Herwig/Vertices/FFWVertex /Herwig/Vertices/FFWMatchboxVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFWMatchboxVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/WWHVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/WWWVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/HGGVertex
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/HHHVertex
cp /Herwig/Vertices/FFHVertex /Herwig/Vertices/TTHVertex
set /Herwig/Vertices/TTHVertex:Fermion 6
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/TTHVertex
cp /Herwig/Vertices/FFHVertex /Herwig/Vertices/BBHVertex
set /Herwig/Vertices/BBHVertex:Fermion 5
cp /Herwig/Vertices/FFHVertex /Herwig/Vertices/TauTauHVertex
set /Herwig/Vertices/TauTauHVertex:Fermion 15
insert DiagramGenerator:Vertices 0 /Herwig/Vertices/TauTauHVertex
cp /Herwig/Vertices/FFHVertex /Herwig/Vertices/MuMuHVertex
set /Herwig/Vertices/MuMuHVertex:Fermion 13
create Herwig::ProcessData ProcessData
set /Herwig/MatrixElements/Matchbox/Factory:DiagramGenerator DiagramGenerator
set /Herwig/MatrixElements/Matchbox/Factory:ProcessData ProcessData
################################################################################
# Setup jet cuts
################################################################################
cd /Herwig/Cuts
create Herwig::MatchboxFactoryMatcher MatchboxJetMatcher
set MatchboxJetMatcher:Group j
create ThePEG::FastJetFinder JetFinder
set JetFinder:UnresolvedMatcher MatchboxJetMatcher
set JetFinder:Variant AntiKt
set JetFinder:RecombinationScheme E
set JetFinder:Mode Inclusive
set JetFinder:ConeRadius 0.7
create ThePEG::JetRegion FirstJet
set FirstJet:PtMin 20.*GeV
do FirstJet:YRange -5.0 5.0
set FirstJet:Fuzzy Yes
set FirstJet:EnergyCutWidth 4.0*GeV
set FirstJet:RapidityCutWidth 0.4
insert FirstJet:Accepts[0] 1
create ThePEG::JetRegion SecondJet
set SecondJet:PtMin 20.*GeV
do SecondJet:YRange -5.0 5.0
set SecondJet:Fuzzy Yes
set SecondJet:EnergyCutWidth 4.0*GeV
set SecondJet:RapidityCutWidth 0.4
insert SecondJet:Accepts[0] 2
create ThePEG::JetRegion ThirdJet
set ThirdJet:PtMin 20.*GeV
do ThirdJet:YRange -5.0 5.0
set ThirdJet:Fuzzy Yes
set ThirdJet:EnergyCutWidth 4.0*GeV
set ThirdJet:RapidityCutWidth 0.4
insert ThirdJet:Accepts[0] 3
create ThePEG::JetRegion FourthJet
set FourthJet:PtMin 20.*GeV
do FourthJet:YRange -5.0 5.0
set FourthJet:Fuzzy Yes
set FourthJet:EnergyCutWidth 4.0*GeV
set FourthJet:RapidityCutWidth 0.4
insert FourthJet:Accepts[0] 4
create ThePEG::FuzzyTheta FuzzyTheta
set FuzzyTheta:EnergyWidth 4.0*GeV
set FuzzyTheta:RapidityWidth 0.4
set FuzzyTheta:AngularWidth 0.4
create ThePEG::NJetsCut NJetsCut
set NJetsCut:UnresolvedMatcher MatchboxJetMatcher
set NJetsCut:NJetsMin 2
create ThePEG::JetCuts JetCuts
set JetCuts:UnresolvedMatcher MatchboxJetMatcher
set JetCuts:Ordering OrderPt
create Herwig::IdentifiedParticleCut IdentifiedParticleCut
cp IdentifiedParticleCut LeptonCut
set LeptonCut:Matcher /Herwig/Matchers/Lepton
cp IdentifiedParticleCut ChargedLeptonCut
set ChargedLeptonCut:Matcher /Herwig/Matchers/ChargedLepton
cp IdentifiedParticleCut BottomQuarkCut
set BottomQuarkCut:Matcher /Herwig/Matchers/Bottom
cp IdentifiedParticleCut TopQuarkCut
set TopQuarkCut:Matcher /Herwig/Matchers/Top
cp IdentifiedParticleCut WBosonCut
set WBosonCut:Matcher /Herwig/Matchers/WBoson
cp IdentifiedParticleCut ZBosonCut
set ZBosonCut:Matcher /Herwig/Matchers/ZBoson
cp IdentifiedParticleCut HiggsBosonCut
set HiggsBosonCut:Matcher /Herwig/Matchers/HiggsBoson
cp IdentifiedParticleCut PhotonCut
set PhotonCut:Matcher /Herwig/Matchers/Photon
create Herwig::FrixionePhotonSeparationCut PhotonIsolationCut
set PhotonIsolationCut:UnresolvedMatcher MatchboxJetMatcher
create Herwig::MatchboxDeltaRCut MatchboxDeltaRCut
cp MatchboxDeltaRCut LeptonDeltaRCut
set LeptonDeltaRCut:FirstMatcher /Herwig/Matchers/Lepton
set LeptonDeltaRCut:SecondMatcher /Herwig/Matchers/Lepton
cp MatchboxDeltaRCut ChargedLeptonDeltaRCut
set ChargedLeptonDeltaRCut:FirstMatcher /Herwig/Matchers/ChargedLepton
set ChargedLeptonDeltaRCut:SecondMatcher /Herwig/Matchers/ChargedLepton
create Herwig::InvariantMassCut InvariantMassCut
cp InvariantMassCut LeptonPairMassCut
set LeptonPairMassCut:FirstMatcher /Herwig/Matchers/Lepton
set LeptonPairMassCut:SecondMatcher /Herwig/Matchers/Lepton
cp InvariantMassCut ChargedLeptonPairMassCut
set ChargedLeptonPairMassCut:FirstMatcher /Herwig/Matchers/ChargedLepton
set ChargedLeptonPairMassCut:SecondMatcher /Herwig/Matchers/ChargedLepton
create Herwig::MissingPtCut MissingPtCut
set MissingPtCut:Matcher /Herwig/Matchers/Neutrino
################################################################################
# Setup scale choices
################################################################################
cd /Herwig/MatrixElements/Matchbox
mkdir Scales
cd Scales
create Herwig::MatchboxScaleChoice SHatScale
cp SHatScale FixedScale
set FixedScale:FixedScale 100.*GeV
create Herwig::MatchboxPtScale MaxJetPtScale
set MaxJetPtScale:JetFinder /Herwig/Cuts/JetFinder
create Herwig::MatchboxLeptonMassScale LeptonPairMassScale
create Herwig::MatchboxLeptonPtScale LeptonPairPtScale
create Herwig::MatchboxHtScale HTScale
create Herwig::MatchboxTopMassScale TopPairMassScale
create Herwig::MatchboxTopMTScale TopPairMTScale
create Herwig::MatchboxTopLinearSumMTScale TopPairLinearMTScale
create Herwig::MatchboxTopIndividualMTScale TopPairIndividualMTScale
+create Herwig::MatchboxTriVecScales TriVecScale
set HTScale:JetFinder /Herwig/Cuts/JetFinder
set HTScale:IncludeMT No
set HTScale:JetPtCut 15.*GeV
cp HTScale HTPrimeScale
set HTPrimeScale:IncludeMT Yes
set HTPrimeScale:JetPtCut 15.*GeV
+set TriVecScale:JetFinder /Herwig/Cuts/JetFinder
cp LeptonPairMassScale LeptonQ2Scale
set /Herwig/MatrixElements/Matchbox/Factory:ScaleChoice LeptonPairMassScale
cd /
diff --git a/src/snippets/DipoleShowerFourFlavours.in b/src/snippets/DipoleShowerFourFlavours.in
--- a/src/snippets/DipoleShowerFourFlavours.in
+++ b/src/snippets/DipoleShowerFourFlavours.in
@@ -1,23 +1,27 @@
# -*- ThePEG-repository -*-
cd /Herwig/Particles
set d:NominalMass 0*GeV
set dbar:NominalMass 0*GeV
set u:NominalMass 0*GeV
set ubar:NominalMass 0*GeV
set s:NominalMass 0*GeV
set sbar:NominalMass 0*GeV
set c:NominalMass 0*GeV
set cbar:NominalMass 0*GeV
cd /Herwig/DipoleShower/Kernels
set IFgx2bbbarxDipoleKernel:UseKernel No
set IFgx2bbarbxDipoleKernel:UseKernel No
set IFMgx2bbbarxDipoleKernel:UseKernel No
set IFMgx2bbarbxDipoleKernel:UseKernel No
set IIgx2bbbarxDipoleKernel:UseKernel No
-set IIgx2bbarbxDipoleKernel:UseKernel No
\ No newline at end of file
+set IIgx2bbarbxDipoleKernel:UseKernel No
+
+
+cd /Herwig/UnderlyingEvent
+set MEQCD2to2Fast:MaximumFlavour 4