diff --git a/Decay/Perturbative/SMWFermionsPOWHEGDecayer.cc b/Decay/Perturbative/SMWFermionsPOWHEGDecayer.cc
--- a/Decay/Perturbative/SMWFermionsPOWHEGDecayer.cc
+++ b/Decay/Perturbative/SMWFermionsPOWHEGDecayer.cc
@@ -1,636 +1,637 @@
//-*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SMWFermionsPOWHEGDecayer class.
//
#include "SMWFermionsPOWHEGDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <numeric>
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Shower/Core/Couplings/ShowerAlpha.h"
using namespace Herwig;
SMWFermionsPOWHEGDecayer::SMWFermionsPOWHEGDecayer()
: CF_(4./3.), pTmin_(1.*GeV)
{ }
IBPtr SMWFermionsPOWHEGDecayer::clone() const {
return new_ptr(*this);
}
IBPtr SMWFermionsPOWHEGDecayer::fullclone() const {
return new_ptr(*this);
}
void SMWFermionsPOWHEGDecayer::persistentOutput(PersistentOStream & os) const {
os << FFGVertex_ << FFWVertex_ << gluon_ << ounit( pTmin_, GeV );
}
void SMWFermionsPOWHEGDecayer::persistentInput(PersistentIStream & is, int) {
is >> FFGVertex_ >> FFWVertex_ >> gluon_ >> iunit( pTmin_, GeV );
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SMWFermionsPOWHEGDecayer,SMWDecayer>
describeHerwigSMWFermionsPOWHEGDecayer("Herwig::SMWFermionsPOWHEGDecayer", "HwPerturbativeDecay.so");
void SMWFermionsPOWHEGDecayer::Init() {
static ClassDocumentation<SMWFermionsPOWHEGDecayer> documentation
("There is no documentation for the SMWFermionsPOWHEGDecayer class");
static Parameter<SMWFermionsPOWHEGDecayer, Energy> interfacePtMin
("minpT",
"The pt cut on hardest emision generation",
&SMWFermionsPOWHEGDecayer::pTmin_, GeV, 1.*GeV, 0*GeV, 100000.0*GeV,
false, false, Interface::limited);
}
RealEmissionProcessPtr SMWFermionsPOWHEGDecayer::
generateHardest(RealEmissionProcessPtr born) {
assert(born->bornOutgoing().size()==2);
// check coloured
if(!born->bornOutgoing()[0]->dataPtr()->coloured()) return RealEmissionProcessPtr();
// extract required info
partons_.resize(2);
quark_.resize(2);
vector<PPtr> hardProcess;
wboson_ = born->bornIncoming()[0];
hardProcess.push_back(wboson_);
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
partons_[ix] = born->bornOutgoing()[ix]->dataPtr();
quark_[ix] = born->bornOutgoing()[ix]->momentum();
quark_[ix].setMass(partons_[ix]->mass());
hardProcess.push_back(born->bornOutgoing()[ix]);
}
bool order = partons_[0]->id()<0;
if(order) {
swap(partons_[0] ,partons_[1] );
swap(quark_[0] ,quark_[1] );
swap(hardProcess[1],hardProcess[2]);
}
gauge_.setMass(0.*MeV);
// Get the W boson mass.
mw2_ = (quark_[0] + quark_[1]).m2();
// Generate emission and set _quark[0,1] and _gauge to be the
// momenta of q, qbar and g after the hardest emission:
if(!getEvent(hardProcess)) {
born->pT()[ShowerInteraction::QCD] = pTmin_;
return born;
}
// Ensure the energies are greater than the constituent masses:
- for (int i=0; i<2; i++) {
- if (quark_[i].e() < partons_[i]->constituentMass()) return RealEmissionProcessPtr();
- }
- if (gauge_.e() < gluon_ ->constituentMass()) return RealEmissionProcessPtr();
+ // for (int i=0; i<2; i++) {
+ // if (quark_[i].e() < partons_[i]->constituentMass()) return RealEmissionProcessPtr();
+ // }
+ // if (gauge_.e() < gluon_ ->constituentMass()) return RealEmissionProcessPtr();
// set masses
quark_[0].setMass( partons_[0]->mass() );
quark_[1].setMass( partons_[1]->mass() );
gauge_ .setMass( ZERO );
// // assign the emitter based on evolution scales
unsigned int iemitter = quark_[0]*gauge_ > quark_[1]*gauge_ ? 2 : 1;
unsigned int ispectator = iemitter==1 ? 1 : 2;
// create new partices and insert
PPtr wboson = wboson_->dataPtr()->produceParticle(wboson_->momentum());
born->incoming().push_back(wboson);
PPtr newq = partons_[0]->produceParticle(quark_[0]);
PPtr newa = partons_[1]->produceParticle(quark_[1]);
PPtr newg = gluon_->produceParticle(gauge_);
// make colour connections
newg->colourNeighbour(newq);
newa->colourNeighbour(newg);
// insert in output structure
if(!order) {
born->outgoing().push_back(newq);
born->outgoing().push_back(newa);
}
else {
born->outgoing().push_back(newa);
born->outgoing().push_back(newq);
swap(iemitter,ispectator);
}
born->outgoing().push_back(newg);
born->emitter (iemitter );
born->spectator(ispectator);
born->emitted (3);
born->pT()[ShowerInteraction::QCD] = pT_;
// return process
born->interaction(ShowerInteraction::QCD);
return born;
}
double SMWFermionsPOWHEGDecayer::
me2(const int ichan, const Particle & part,
const ParticleVector & decay, MEOption meopt) const {
// leading-order result
double output = SMWDecayer::me2(ichan,part,decay,meopt);
// check decay products coloured, otherwise return
if(!decay[0]->dataPtr()->coloured()) return output;
// inital masses, couplings etc
// W mass
mW_ = part.mass();
// strong coupling
aS_ = SM().alphaS(sqr(mW_));
// reduced mass
double mu1_ = (decay[0]->dataPtr()->mass())/mW_;
double mu2_ = (decay[1]->dataPtr()->mass())/mW_;
// scale
scale_ = sqr(mW_);
// now for the nlo loop correction
double virt = CF_*aS_/Constants::pi;
// now for the real correction
double realFact=0.;
for(int iemit=0;iemit<2;++iemit) {
double phi = UseRandom::rnd()*Constants::twopi;
// set the emitter and the spectator
double muj = iemit==0 ? mu1_ : mu2_;
double muk = iemit==0 ? mu2_ : mu1_;
double muj2 = sqr(muj);
double muk2 = sqr(muk);
// calculate y
double yminus = 0.;
double yplus = 1.-2.*muk*(1.-muk)/(1.-muj2-muk2);
double y = yminus + UseRandom::rnd()*(yplus-yminus);
double v = sqrt(sqr(2.*muk2 + (1.-muj2-muk2)*(1.-y))-4.*muk2)
/(1.-muj2-muk2)/(1.-y);
double zplus = (1.+v)*(1.-muj2-muk2)*y/2./(muj2+(1.-muj2-muk2)*y);
double zminus = (1.-v)*(1.-muj2-muk2)*y/2./(muj2+(1.-muj2-muk2)*y);
double z = zminus + UseRandom::rnd()*(zplus-zminus);
double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
// calculate x1,x2,x3,xT
double x2 = 1.-y*(1.-muj2-muk2)-muj2+muk2;
double x1 = 1.+muj2-muk2-z*(x2-2.*muk2);
// copy the particle objects over for calculateRealEmission
vector<PPtr> hardProcess(3);
hardProcess[0] = const_ptr_cast<PPtr>(&part);
hardProcess[1] = decay[0];
hardProcess[2] = decay[1];
realFact = 0.25*jac*sqr(1.-muj2-muk2)/
sqrt((1.-sqr(muj-muk))*(1.-sqr(muj+muk)))/Constants::twopi
*2.*CF_*aS_*calculateRealEmission(x1, x2, hardProcess, phi,
muj, muk, iemit, true);
}
// the born + virtual + real
output = output*(1. + virt + realFact);
return output;
}
double SMWFermionsPOWHEGDecayer::meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter, bool subtract) const {
Lorentz5Momentum q = momenta[1]+momenta[2]+momenta[3];
Energy2 Q2=q.m2();
Energy2 lambda = sqrt((Q2-sqr(momenta[1].mass()+momenta[2].mass()))*
(Q2-sqr(momenta[1].mass()-momenta[2].mass())));
InvEnergy2 D[2];
double lome[2];
for(unsigned int iemit=0;iemit<2;++iemit) {
unsigned int ispect = iemit==0 ? 1 : 0;
Energy2 pipj = momenta[3 ] * momenta[1+iemit ];
Energy2 pipk = momenta[3 ] * momenta[1+ispect];
Energy2 pjpk = momenta[1+iemit] * momenta[1+ispect];
double y = pipj/(pipj+pipk+pjpk);
double z = pipk/( pipk+pjpk);
Energy mij = sqrt(2.*pipj+sqr(momenta[1+iemit].mass()));
Energy2 lamB = sqrt((Q2-sqr(mij+momenta[1+ispect].mass()))*
(Q2-sqr(mij-momenta[1+ispect].mass())));
Energy2 Qpk = q*momenta[1+ispect];
Lorentz5Momentum pkt =
lambda/lamB*(momenta[1+ispect]-Qpk/Q2*q)
+0.5/Q2*(Q2+sqr(momenta[1+ispect].mass())-sqr(momenta[1+ispect].mass()))*q;
Lorentz5Momentum pijt =
q-pkt;
double muj = momenta[1+iemit ].mass()/sqrt(Q2);
double muk = momenta[1+ispect].mass()/sqrt(Q2);
double vt = sqrt((1.-sqr(muj+muk))*(1.-sqr(muj-muk)))/(1.-sqr(muj)-sqr(muk));
double v = sqrt(sqr(2.*sqr(muk)+(1.-sqr(muj)-sqr(muk))*(1.-y))-4.*sqr(muk))
/(1.-y)/(1.-sqr(muj)-sqr(muk));
// dipole term
D[iemit] = 0.5/pipj*(2./(1.-(1.-z)*(1.-y))
- -vt/v*(2.-z+sqr(momenta[1+iemit].mass())/pipj));
+ -vt/v*(2.-z+sqr(momenta[1+iemit].mass())/pipj));
// matrix element
vector<Lorentz5Momentum> lomom(3);
lomom[0] = momenta[0];
if(iemit==0) {
lomom[1] = pijt;
lomom[2] = pkt ;
}
else {
lomom[2] = pijt;
lomom[1] = pkt ;
}
lome[iemit] = loME(partons,lomom);
}
InvEnergy2 ratio = realME(partons,momenta)*abs(D[iemitter])
/(abs(D[0]*lome[0])+abs(D[1]*lome[1]));
if(subtract)
return Q2*(ratio-2.*D[iemitter]);
else
return Q2*ratio;
}
double SMWFermionsPOWHEGDecayer::loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
VectorWaveFunction win (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
}
for(unsigned int ix=0;ix<3;++ix){
win.reset(ix);
vin.push_back(win);
}
// temporary storage of the different diagrams
// sum over helicities to get the matrix element
double total(0.);
for(unsigned int inhel=0;inhel<3;++inhel) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
Complex diag1 = FFWVertex()->evaluate(scale_,aout[outhel2],fout[outhel1],vin[inhel]);
total += norm(diag1);
}
}
}
// return the answer
return total;
}
InvEnergy2 SMWFermionsPOWHEGDecayer::realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
vector<VectorWaveFunction> gout;
VectorWaveFunction win (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
VectorWaveFunction gluon(momenta[3],partons[3],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
gluon.reset(2*ix);
gout.push_back(gluon);
}
for(unsigned int ix=0;ix<3;++ix){
win.reset(ix);
vin.push_back(win);
}
vector<Complex> diag(2,0.);
double total(0.);
for(unsigned int inhel1=0;inhel1<3;++inhel1) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
for(unsigned int outhel3=0;outhel3<2;++outhel3) {
SpinorBarWaveFunction off1 =
FFGVertex()->evaluate(scale_,3,partons[1],fout[outhel1],gout[outhel3]);
diag[0] = FFWVertex()->evaluate(scale_,aout[outhel2],off1,vin[inhel1]);
SpinorWaveFunction off2 =
FFGVertex()->evaluate(scale_,3,partons[2],aout[outhel2],gout[outhel3]);
diag[1] = FFWVertex()->evaluate(scale_,off2,fout[outhel1],vin[inhel1]);
// sum of diagrams
Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.));
// me2
total += norm(sum);
}
}
}
}
// divide out the coupling
total /= norm(FFGVertex()->norm());
// return the total
return total*UnitRemoval::InvE2;
}
void SMWFermionsPOWHEGDecayer::doinit() {
// cast the SM pointer to the Herwig SM pointer
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(!hwsm) throw InitException()
<< "Wrong type of StandardModel object in "
<< "SMWFermionsPOWHEGDecayer::doinit() "
<< "the Herwig version must be used."
<< Exception::runerror;
// cast the vertices
FFWVertex_ = hwsm->vertexFFW();
FFWVertex_->init();
FFGVertex_ = hwsm->vertexFFG();
FFGVertex_->init();
SMWDecayer::doinit();
gluon_ = getParticleData(ParticleID::g);
}
bool SMWFermionsPOWHEGDecayer::getEvent(vector<PPtr> hardProcess) {
vector<Energy> particleMass;
for(unsigned int ix=0;ix<hardProcess.size();++ix) {
if(abs(hardProcess[ix]->id())==ParticleID::Wplus) {
mW_ = hardProcess[ix]->mass();
}
else {
particleMass.push_back(hardProcess[ix]->mass());
}
}
if (particleMass.size()!=2) {
throw Exception()
<< "Number of outgoing particles is not equal to 2 in "
<< "SMWFermionPOWHEGDecayer::getEvent()"
<< Exception::runerror;
}
// reduced mass
mu1_ = particleMass[0]/mW_;
mu2_ = particleMass[1]/mW_;
// scale
scale_ = sqr(mW_);
// max pT
Energy pTmax = 0.5*sqrt(mw2_);
if(pTmax<pTmin_) return false;
// Define over valued y_max & y_min according to the associated pt_min cut.
- double ymax = acosh(pTmax/pTmin_);
+ //double ymax = acosh(pTmax/pTmin_);
+ double ymax = 10.;
double ymin = -ymax;
// pt of the emmission
pT_ = pTmax;
// prefactor
double overEst = 4.;
double prefactor = overEst*alphaS()->overestimateValue()*CF_*
(ymax-ymin)/Constants::twopi;
// loop to generate the pt and rapidity
bool reject;
//arrays to hold the temporary probabilities whilst the for loop progresses
double probTemp[2][2]={{0.,0.},{0.,0.}};
probTemp[0][0]=probTemp[0][1]=probTemp[1][0]=probTemp[1][1]=0.;
double x1Solution[2][2] = {{0.,0.},{0.,0.}};
double x2Solution[2][2] = {{0.,0.},{0.,0.}};
double x3Solution[2] = {0.,0.};
Energy pT[2] = {pTmax,pTmax};
double yTemp[2] = {0.,0.};
double phi = 0.;
// do the competition
for(int i=0; i<2; i++) {
// set the emitter and the spectator
double muj = i==0 ? mu1_ : mu2_;
double muk = i==0 ? mu2_ : mu1_;
double muj2 = sqr(muj);
double muk2 = sqr(muk);
do {
// generation of phi
phi = UseRandom::rnd() * Constants::twopi;
// reject the emission
reject = true;
// generate pt
pT[i] *= pow(UseRandom::rnd(),1./prefactor);
if(pT[i]<pTmin_) {
pT[i] = -GeV;
break;
}
// generate xT
double xT2 = sqr(2./mW_*pT[i]);
// generate y
yTemp[i] = ymin + UseRandom::rnd()*(ymax-ymin);
// generate x3 & x1 from pT & y
double x1Plus = 1-muk2+muj2;
double x1Minus = 2.*muj;
x3Solution[i] = 2.*pT[i]*cosh(yTemp[i])/mW_;
// prefactor
double weightPrefactor = 0.5/sqrt((1.-sqr(muj-muk))*(1.-sqr(muj+muk)))/overEst;
// calculate x1 & x2 solutions
double discrim2 = (-sqr(x3Solution[i])+xT2)*
(xT2*muk2+2.*x3Solution[i]-sqr(muj2)+2.*muk2+2.*muj2-sqr(x3Solution[i])-1.
+2.*muj2*muk2-sqr(muk2)-2.*muk2*x3Solution[i]-2.*muj2*x3Solution[i]);
// check discrim2 is > 0
if( discrim2 < ZERO) continue;
double fact1 =2.*sqr(x3Solution[i])-4.*muk2-6.*x3Solution[i]+4.*muj2-xT2*x3Solution[i]
+2.*xT2-2.*muj2*x3Solution[i]+2.*muk2*x3Solution[i]+4.;
double fact2 = (4.-4.*x3Solution[i]+xT2);
double discriminant = sqrt(discrim2);
// two solns for x1
x1Solution[i][0] = (fact1 + 2.*discriminant)/fact2;
x1Solution[i][1] = (fact1 - 2.*discriminant)/fact2;
bool found = false;
for(unsigned int j=0;j<2;++j) {
// calculate x2
x2Solution[i][j] = 2.-x3Solution[i]-x1Solution[i][j];
// set limits on x2
double root = max(0.,sqr(x1Solution[i][j])-4.*muj2);
root = sqrt(root);
double x2Plus = 1.+muk2-muj2
-0.5*(1.-x1Solution[i][j]+muj2-muk2)/(1.-x1Solution[i][j]+muj2)
*(x1Solution[i][j]-2.*muj2-root);
double x2Minus = 1.+muk2-muj2
-0.5*(1.-x1Solution[i][j]+muj2-muk2)/(1.-x1Solution[i][j]+muj2)
*(x1Solution[i][j]-2.*muj2+root);
-
+
if(x1Solution[i][j]>=x1Minus && x1Solution[i][j]<=x1Plus &&
x2Solution[i][j]>=x2Minus && x2Solution[i][j]<=x2Plus &&
checkZMomenta(x1Solution[i][j], x2Solution[i][j], x3Solution[i], yTemp[i], pT[i],
muj, muk)) {
- probTemp[i][j] = weightPrefactor*pT[i]*
- calculateJacobian(x1Solution[i][j], x2Solution[i][j], pT[i], muj, muk)*
- calculateRealEmission(x1Solution[i][j], x2Solution[i][j],
+ probTemp[i][j] = weightPrefactor*pT[i]*
+ abs(calculateJacobian(x1Solution[i][j], x2Solution[i][j], pT[i], muj, muk))*
+ calculateRealEmission(x1Solution[i][j], x2Solution[i][j],
hardProcess, phi, muj, muk, i, false);
found = true;
}
else {
probTemp[i][j] = 0.;
}
}
if(!found) continue;
// alpha S piece
double wgt = (probTemp[i][0]+probTemp[i][1])*alphaS()->ratio(sqr(pT[i]));
// matrix element weight
reject = UseRandom::rnd()>wgt;
}
while(reject);
} // end of emitter for loop
// no emission
if(pT[0]<ZERO&&pT[1]<ZERO) return false;
//pick the spectator and x1 x2 values
double x1,x2,y;
// particle 1 emits, particle 2 spectates
unsigned int iemit=0;
if(pT[0]>pT[1]){
pT_ = pT[0];
y=yTemp[0];
if(probTemp[0][0]>UseRandom::rnd()*(probTemp[0][0]+probTemp[0][1])) {
x1 = x1Solution[0][0];
x2 = x2Solution[0][0];
}
else {
x1 = x1Solution[0][1];
x2 = x2Solution[0][1];
}
}
// particle 2 emits, particle 1 spectates
else {
iemit=1;
pT_ = pT[1];
y=yTemp[1];
if(probTemp[1][0]>UseRandom::rnd()*(probTemp[1][0]+probTemp[1][1])) {
x1 = x1Solution[1][0];
x2 = x2Solution[1][0];
}
else {
x1 = x1Solution[1][1];
x2 = x2Solution[1][1];
}
}
// find spectator
unsigned int ispect = iemit == 0 ? 1 : 0;
double muk = iemit == 0 ? mu2_ : mu1_;
double muk2 = sqr(muk);
double muj = iemit == 0 ? mu1_ : mu2_;
double muj2 = sqr(muj);
double xT2 = sqr(2./mW_*pT_);
// Find the boost from the lab to the c.o.m with the spectator
// along the -z axis, and then invert it.
LorentzRotation eventFrame( ( quark_[0] + quark_[1] ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*quark_[ispect];
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() - Constants::pi );
eventFrame.invert();
// spectator
quark_[ispect].setT( 0.5*x2*mW_ );
quark_[ispect].setX( ZERO );
quark_[ispect].setY( ZERO );
quark_[ispect].setZ( -sqrt(0.25*mw2_*x2*x2-mw2_*muk2) );
// gluon
gauge_.setT( pT_*cosh(y) );
gauge_.setX( pT_*cos(phi) );
gauge_.setY( pT_*sin(phi) );
gauge_.setZ( pT_*sinh(y) );
gauge_.setMass(ZERO);
// emitter
quark_[iemit].setX( -pT_*cos(phi) );
quark_[iemit].setY( -pT_*sin(phi) );
quark_[iemit].setZ( 0.5*mW_*sqrt(sqr(x1)-xT2-4.*muj2) );
if(sqrt(0.25*mw2_*x2*x2-mw2_*muk2)-pT_*sinh(y)<ZERO)
quark_[iemit].setZ(-quark_[iemit].z());
quark_[iemit].setT( 0.5*mW_*x1 );
// boost constructed vectors into the event frame
quark_[0] = eventFrame * quark_[0];
quark_[1] = eventFrame * quark_[1];
gauge_ = eventFrame * gauge_;
// need to reset masses because for whatever reason the boost
// touches the mass component of the five-vector and can make
// zero mass objects acquire a floating point negative mass(!).
gauge_.setMass( ZERO );
quark_[iemit] .setMass(partons_[iemit ]->mass());
quark_[ispect].setMass(partons_[ispect]->mass());
return true;
}
InvEnergy SMWFermionsPOWHEGDecayer::calculateJacobian(double x1, double x2, Energy pT,
double muj, double muk) const{
double xPerp = abs(2.*pT/mW_);
Energy jac = mW_/xPerp*fabs((x2*sqr(muj)+2.*sqr(muk)*x1
+sqr(muk)*x2-x1*x2-sqr(x2)+x2)/pow((sqr(x2)-4.*sqr(muk)),1.5));
return 1./jac; //jacobian as defined is dptdy=jac*dx1dx2, therefore we have to divide by it
}
bool SMWFermionsPOWHEGDecayer::checkZMomenta(double x1, double x2, double x3,
double y, Energy pT, double muj,
double muk) const {
double xPerp2 = 4.*pT*pT/mW_/mW_;
double root1 = sqrt(max(0.,sqr(x2)-4.*sqr(muk)));
double root2 = sqrt(max(0.,sqr(x1)-xPerp2 - 4.*sqr(muj)));
static double tolerance = 1e-6;
bool isMomentaReconstructed = false;
if(pT*sinh(y) > ZERO) {
if(abs(-root1 + sqrt(sqr(x3)-xPerp2) + root2) <= tolerance ||
abs(-root1 + sqrt(sqr(x3)-xPerp2) - root2) <= tolerance)
isMomentaReconstructed=true;
}
else if(pT*sinh(y) < ZERO){
if(abs(-root1 - sqrt(sqr(x3)-xPerp2) + root2) <= tolerance ||
abs(-root1 - sqrt(sqr(x3)-xPerp2) - root2) <= tolerance)
isMomentaReconstructed=true;
}
else
if(abs(-root1+ sqrt(sqr(x1)-xPerp2 - 4.*(muj))) <= tolerance)
isMomentaReconstructed=true;
return isMomentaReconstructed;
}
double SMWFermionsPOWHEGDecayer::calculateRealEmission(double x1, double x2,
vector<PPtr> hardProcess,
double phi, double muj,
double muk, int iemit,
bool subtract) const {
// make partons data object for meRatio
vector<cPDPtr> partons (3);
for(int ix=0; ix<3; ++ix)
partons[ix] = hardProcess[ix]->dataPtr();
partons.push_back(gluon_);
// calculate x3
double x3 = 2.-x1-x2;
double xT = sqrt(max(0.,sqr(x3)-0.25*sqr(sqr(x2)+sqr(x3)-sqr(x1)-4.*sqr(muk)+4.*sqr(muj))
/(sqr(x2)-4.*sqr(muk))));
// calculate the momenta
Energy M = mW_;
Lorentz5Momentum pspect(ZERO,ZERO,-0.5*M*sqrt(max(sqr(x2)-4.*sqr(muk),0.)),
0.5*M*x2,M*muk);
Lorentz5Momentum pemit (-0.5*M*xT*cos(phi),-0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x1)-sqr(xT)-4.*sqr(muj),0.)),
0.5*M*x1,M*muj);
Lorentz5Momentum pgluon(0.5*M*xT*cos(phi), 0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x3)-sqr(xT),0.)),0.5*M*x3,ZERO);
if(abs(pspect.z()+pemit.z()-pgluon.z())/M<1e-6)
pgluon.setZ(-pgluon.z());
else if(abs(pspect.z()-pemit.z()+pgluon.z())/M<1e-6)
pemit .setZ(- pemit.z());
// loop over the possible emitting partons
double realwgt(0.);
// boost and rotate momenta
LorentzRotation eventFrame( ( hardProcess[1]->momentum() +
hardProcess[2]->momentum() ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*hardProcess[iemit+1]->momentum();
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() );
eventFrame.invert();
vector<Lorentz5Momentum> momenta(3);
momenta[0] = hardProcess[0]->momentum();
if(iemit==0) {
momenta[2] = eventFrame*pspect;
momenta[1] = eventFrame*pemit ;
}
else {
momenta[1] = eventFrame*pspect;
momenta[2] = eventFrame*pemit ;
}
momenta.push_back(eventFrame*pgluon);
// calculate the weight
if(1.-x1>1e-5 && 1.-x2>1e-5)
realwgt += meRatio(partons,momenta,iemit,subtract);
return realwgt;
}
diff --git a/Decay/Perturbative/SMZFermionsPOWHEGDecayer.cc b/Decay/Perturbative/SMZFermionsPOWHEGDecayer.cc
--- a/Decay/Perturbative/SMZFermionsPOWHEGDecayer.cc
+++ b/Decay/Perturbative/SMZFermionsPOWHEGDecayer.cc
@@ -1,742 +1,743 @@
//-*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SMZFermionsPOWHEGDecayer class.
//
#include "SMZFermionsPOWHEGDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <numeric>
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Shower/Core/Couplings/ShowerAlpha.h"
using namespace Herwig;
SMZFermionsPOWHEGDecayer::SMZFermionsPOWHEGDecayer()
: CF_(4./3.), pTmin_(1.*GeV)
{ }
IBPtr SMZFermionsPOWHEGDecayer::clone() const {
return new_ptr(*this);
}
IBPtr SMZFermionsPOWHEGDecayer::fullclone() const {
return new_ptr(*this);
}
void SMZFermionsPOWHEGDecayer::persistentOutput(PersistentOStream & os) const {
os << FFGVertex_ << FFZVertex_ << gluon_ << ounit( pTmin_, GeV );
}
void SMZFermionsPOWHEGDecayer::persistentInput(PersistentIStream & is, int) {
is >> FFGVertex_ >> FFZVertex_ >> gluon_ >> iunit( pTmin_, GeV );
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SMZFermionsPOWHEGDecayer,SMZDecayer>
describeHerwigSMZFermionsPOWHEGDecayer("Herwig::SMZFermionsPOWHEGDecayer", "HwPerturbativeDecay.so");
void SMZFermionsPOWHEGDecayer::Init() {
static ClassDocumentation<SMZFermionsPOWHEGDecayer> documentation
("There is no documentation for the SMZFermionsPOWHEGDecayer class");
static Parameter<SMZFermionsPOWHEGDecayer, Energy> interfacePtMin
("minpT",
"The pt cut on hardest emision generation",
&SMZFermionsPOWHEGDecayer::pTmin_, GeV, 1.*GeV, 0*GeV, 100000.0*GeV,
false, false, Interface::limited);
}
RealEmissionProcessPtr SMZFermionsPOWHEGDecayer::
generateHardest(RealEmissionProcessPtr born) {
assert(born->bornOutgoing().size()==2);
// check coloured
if(!born->bornOutgoing()[0]->dataPtr()->coloured()) return RealEmissionProcessPtr();
// extract required info
partons_.resize(2);
quark_.resize(2);
vector<PPtr> hardProcess;
zboson_ = born->bornIncoming()[0];
hardProcess.push_back(zboson_);
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
partons_[ix] = born->bornOutgoing()[ix]->dataPtr();
quark_[ix] = born->bornOutgoing()[ix]->momentum();
quark_[ix].setMass(partons_[ix]->mass());
hardProcess.push_back(born->bornOutgoing()[ix]);
}
bool order = partons_[0]->id()<0;
if(order) {
swap(partons_[0] ,partons_[1] );
swap(quark_[0] ,quark_[1] );
swap(hardProcess[1],hardProcess[2]);
}
gauge_.setMass(0.*MeV);
// Get the Z boson mass.
mz2_ = (quark_[0] + quark_[1]).m2();
// Generate emission and set _quark[0,1] and _gauge to be the
// momenta of q, qbar and g after the hardest emission:
if(!getEvent(hardProcess)) {
born->pT()[ShowerInteraction::QCD] = pTmin_;
return born;
}
// Ensure the energies are greater than the constituent masses:
- for (int i=0; i<2; i++) {
- if (quark_[i].e() < partons_[i]->constituentMass()) return RealEmissionProcessPtr();
- }
- if (gauge_.e() < gluon_ ->constituentMass()) return RealEmissionProcessPtr();
+ // for (int i=0; i<2; i++) {
+ // if (quark_[i].e() < partons_[i]->constituentMass()) return RealEmissionProcessPtr();
+ // }
+ // if (gauge_.e() < gluon_ ->constituentMass()) return RealEmissionProcessPtr();
// set masses
quark_[0].setMass( partons_[0]->mass() );
quark_[1].setMass( partons_[1]->mass() );
gauge_ .setMass( ZERO );
// assign the emitter based on evolution scales
unsigned int iemitter = quark_[0]*gauge_ > quark_[1]*gauge_ ? 2 : 1;
unsigned int ispectator = iemitter==1 ? 1 : 2;
// create new partices and insert
PPtr zboson = zboson_->dataPtr()->produceParticle(zboson_->momentum());
born->incoming().push_back(zboson);
PPtr newq = partons_[0]->produceParticle(quark_[0]);
PPtr newa = partons_[1]->produceParticle(quark_[1]);
PPtr newg = gluon_->produceParticle(gauge_);
// make colour connections
newg->colourNeighbour(newq);
newa->colourNeighbour(newg);
// insert in output structure
if(!order) {
born->outgoing().push_back(newq);
born->outgoing().push_back(newa);
}
else {
born->outgoing().push_back(newa);
born->outgoing().push_back(newq);
swap(iemitter,ispectator);
}
born->outgoing().push_back(newg);
born->emitter (iemitter );
born->spectator(ispectator);
born->emitted (3);
born->pT()[ShowerInteraction::QCD] = pT_;
// return process
born->interaction(ShowerInteraction::QCD);
return born;
}
double SMZFermionsPOWHEGDecayer::
me2(const int ichan, const Particle & part,
const ParticleVector & decay, MEOption meopt) const {
// leading-order result
double output = SMZDecayer::me2(ichan,part,decay,meopt);
// check decay products coloured, otherwise return
if(!decay[0]->dataPtr()->coloured()) return output;
// inital masses, couplings etc
// fermion mass
Energy particleMass = decay[0]->dataPtr()->mass();
// Z mass
mZ_ = part.mass();
// strong coupling
aS_ = SM().alphaS(sqr(mZ_));
// reduced mass
mu_ = particleMass/mZ_;
mu2_ = sqr(mu_);
// scale
scale_ = sqr(mZ_);
// cast the vertices
tcFFVVertexPtr Zvertex = dynamic_ptr_cast<tcFFVVertexPtr>(FFZVertex());
// compute the spinors
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
vector<VectorWaveFunction> vin;
SpinorBarWaveFunction qkout(decay[0]->momentum(),decay[0]->dataPtr(),outgoing);
SpinorWaveFunction qbout(decay[1]->momentum(),decay[1]->dataPtr(),outgoing);
VectorWaveFunction zin (part.momentum() ,part.dataPtr() ,incoming);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
}
for(unsigned int ix=0;ix<3;++ix){
zin.reset(ix);
vin.push_back(zin);
}
// temporary storage of the different diagrams
// sum over helicities to get the matrix element
double total=0.;
if(mu_!=0.) {
LorentzPolarizationVector momDiff =
(decay[0]->momentum()-decay[1]->momentum())/2./
(decay[0]->momentum().mass()+decay[1]->momentum().mass());
// scalars
Complex scalar1 = zin.wave().dot(momDiff);
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
for(unsigned int inhel=0;inhel<3;++inhel) {
// first the LO bit
Complex diag1 = FFZVertex()->evaluate(scale_,aout[outhel2],fout[outhel1],vin[inhel]);
// extra stuff for NLO
LorentzPolarizationVector left =
aout[outhel2].wave().leftCurrent(fout[outhel1].wave());
LorentzPolarizationVector right =
aout[outhel2].wave().rightCurrent(fout[outhel1].wave());
Complex scalar =
aout[outhel2].wave().scalar(fout[outhel1].wave());
// nlo specific pieces
Complex diag3 =
Complex(0.,1.)*Zvertex->norm()*
(Zvertex->right()*( left.dot(zin.wave())) +
Zvertex-> left()*(right.dot(zin.wave())) -
( Zvertex-> left()+Zvertex->right())*scalar1*scalar);
// nlo piece
total += real(diag1*conj(diag3) + diag3*conj(diag1));
}
}
}
// rescale
total *= UnitRemoval::E2/scale_;
}
else {
total = ZERO;
}
// now for the NLO bit
double mu4 = sqr(mu2_);
double lmu = mu_!=0. ? log(mu_) : 0.;
double v = sqrt(1.-4.*mu2_),v2(sqr(v));
double omv = 4.*mu2_/(1.+v);
double f1,f2,fNS,VNS;
double r = omv/(1.+v);
double lr = mu_!=0. ? log(r) : 0.;
// normal form
if(mu_>1e-4) {
f1 = CF_*aS_/Constants::pi*
( +1. + 3.*log(0.5*(1.+v)) - 1.5*log(0.5*(1.+v2)) + sqr(Constants::pi)/6.
- 0.5*sqr(lr) - (1.+v2)/v*(lr*log(1.+v2) + sqr(Constants::pi)/12.
-0.5*log(4.*mu2_)*lr + 0.25*sqr(lr)));
fNS = -0.5*(1.+2.*v2)*lr/v + 1.5*lr - 2./3.*sqr(Constants::pi) + 0.5*sqr(lr)
+ (1.+v2)/v*(Herwig::Math::ReLi2(r) + sqr(Constants::pi)/3. - 0.25*sqr(lr)
+ lr*log((2.*v/ (1.+v))));
VNS = 1.5*log(0.5*(1.+v2))
+ 0.5*(1.+v2)/v*( 2.*lr*log(2.*(1.+v2)/sqr(1.+v))
+ 2.*Herwig::Math::ReLi2(sqr(r))
- 2.*Herwig::Math::ReLi2(2.*v/(1.+v)) - sqr(Constants::pi)/6.)
+ log(1.-mu_) - 2.*log(1.-2.*mu_) - 4.*mu2_/(1.+v2)*log(mu_/(1.-mu_))
- mu_/(1.-mu_)
+ 4.*(2.*mu2_-mu_)/(1.+v2) + 0.5*sqr(Constants::pi);
f2 = CF_*aS_/Constants::pi*mu2_*lr/v;
}
// small mass limit
else {
f1 = -CF_*aS_/Constants::pi/6.*
( - 6. - 24.*lmu*mu2_ - 15.*mu4 - 12.*mu4*lmu - 24.*mu4*sqr(lmu)
+ 2.*mu4*sqr(Constants::pi) - 12.*mu2_*mu4 - 96.*mu2_*mu4*sqr(lmu)
+ 8.*mu2_*mu4*sqr(Constants::pi) - 80.*mu2_*mu4*lmu);
fNS = - mu2_/18.*( + 36.*lmu - 36. - 45.*mu2_ + 216.*lmu*mu2_ - 24.*mu2_*sqr(Constants::pi)
+ 72.*mu2_*sqr(lmu) - 22.*mu4 + 1032.*mu4 * lmu
- 96.*mu4*sqr(Constants::pi) + 288.*mu4*sqr(lmu));
VNS = - mu2_/1260.*(-6930. + 7560.*lmu + 2520.*mu_ - 16695.*mu2_
+ 1260.*mu2_*sqr(Constants::pi)
+ 12600.*lmu*mu2_ + 1344.*mu_*mu2_ - 52780.*mu4 + 36960.*mu4*lmu
+ 5040.*mu4*sqr(Constants::pi) - 12216.*mu_*mu4);
f2 = CF_*aS_*mu2_/Constants::pi*( 2.*lmu + 4.*mu2_*lmu + 2.*mu2_ + 12.*mu4*lmu + 7.*mu4);
}
// add up bits for f1
f1 += CF_*aS_/Constants::pi*(fNS+VNS);
// now for the real correction
double phi = UseRandom::rnd()*Constants::twopi;
// calculate y
double yminus = 0.;
double yplus = 1.-2.*mu_*(1.-mu_)/(1.-2*mu2_);
double y = yminus + UseRandom::rnd()*(yplus-yminus);
// calculate z
double v1 = sqrt(sqr(2.*mu2_+(1.-2.*mu2_)*(1.-y))-4.*mu2_)/(1.-2.*mu2_)/(1.-y);
double zplus = (1.+v1)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
double zminus = (1.-v1)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
double z = zminus + UseRandom::rnd()*(zplus-zminus);
double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
// calculate x1,x2
double x2 = 1. - y*(1.-2.*mu2_);
double x1 = 1. - z*(x2-2.*mu2_);
// copy the particle objects over for calculateRealEmission
vector<PPtr> hardProcess(3);
hardProcess[0] = const_ptr_cast<PPtr>(&part);
hardProcess[1] = decay[0];
hardProcess[2] = decay[1];
// total real emission contribution
double realFact = 0.25*jac*sqr(1.-2.*mu2_)/
sqrt(1.-4.*mu2_)/Constants::twopi
*2.*CF_*aS_*calculateRealEmission(x1, x2, hardProcess, phi, true);
// the born + virtual + real
output = output*(1. + f1 + realFact) + f2*total;
return output;
}
double SMZFermionsPOWHEGDecayer::meRatio(vector<cPDPtr> partons,
vector<Lorentz5Momentum> momenta,
unsigned int iemitter, bool subtract) const {
Lorentz5Momentum q = momenta[1]+momenta[2]+momenta[3];
Energy2 Q2=q.m2();
Energy2 lambda = sqrt((Q2-sqr(momenta[1].mass()+momenta[2].mass()))*
(Q2-sqr(momenta[1].mass()-momenta[2].mass())));
InvEnergy2 D[2];
double lome[2];
for(unsigned int iemit=0;iemit<2;++iemit) {
unsigned int ispect = iemit==0 ? 1 : 0;
Energy2 pipj = momenta[3 ] * momenta[1+iemit ];
Energy2 pipk = momenta[3 ] * momenta[1+ispect];
Energy2 pjpk = momenta[1+iemit] * momenta[1+ispect];
double y = pipj/(pipj+pipk+pjpk);
double z = pipk/( pipk+pjpk);
Energy mij = sqrt(2.*pipj+sqr(momenta[1+iemit].mass()));
Energy2 lamB = sqrt((Q2-sqr(mij+momenta[1+ispect].mass()))*
(Q2-sqr(mij-momenta[1+ispect].mass())));
Energy2 Qpk = q*momenta[1+ispect];
Lorentz5Momentum pkt =
lambda/lamB*(momenta[1+ispect]-Qpk/Q2*q)
+0.5/Q2*(Q2+sqr(momenta[1+ispect].mass())-sqr(momenta[1+ispect].mass()))*q;
Lorentz5Momentum pijt =
q-pkt;
double muj = momenta[1+iemit ].mass()/sqrt(Q2);
double muk = momenta[1+ispect].mass()/sqrt(Q2);
double vt = sqrt((1.-sqr(muj+muk))*(1.-sqr(muj-muk)))/(1.-sqr(muj)-sqr(muk));
double v = sqrt(sqr(2.*sqr(muk)+(1.-sqr(muj)-sqr(muk))*(1.-y))-4.*sqr(muk))
/(1.-y)/(1.-sqr(muj)-sqr(muk));
// dipole term
D[iemit] = 0.5/pipj*(2./(1.-(1.-z)*(1.-y))
-vt/v*(2.-z+sqr(momenta[1+iemit].mass())/pipj));
// matrix element
vector<Lorentz5Momentum> lomom(3);
lomom[0] = momenta[0];
if(iemit==0) {
lomom[1] = pijt;
lomom[2] = pkt ;
}
else {
lomom[2] = pijt;
lomom[1] = pkt ;
}
lome[iemit] = loME(partons,lomom);
}
InvEnergy2 ratio = realME(partons,momenta)*abs(D[iemitter])
/(abs(D[0]*lome[0])+abs(D[1]*lome[1]));
if(subtract)
return Q2*(ratio-2.*D[iemitter]);
else
return Q2*ratio;
}
double SMZFermionsPOWHEGDecayer::loME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
VectorWaveFunction zin (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
}
for(unsigned int ix=0;ix<3;++ix){
zin.reset(ix);
vin.push_back(zin);
}
// temporary storage of the different diagrams
// sum over helicities to get the matrix element
double total(0.);
for(unsigned int inhel=0;inhel<3;++inhel) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
Complex diag1 = FFZVertex()->evaluate(scale_,aout[outhel2],fout[outhel1],vin[inhel]);
total += norm(diag1);
}
}
}
// return the answer
return total;
}
InvEnergy2 SMZFermionsPOWHEGDecayer::realME(const vector<cPDPtr> & partons,
const vector<Lorentz5Momentum> & momenta) const {
// compute the spinors
vector<VectorWaveFunction> vin;
vector<SpinorWaveFunction> aout;
vector<SpinorBarWaveFunction> fout;
vector<VectorWaveFunction> gout;
VectorWaveFunction zin (momenta[0],partons[0],incoming);
SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
SpinorWaveFunction qbout(momenta[2],partons[2],outgoing);
VectorWaveFunction gluon(momenta[3],partons[3],outgoing);
for(unsigned int ix=0;ix<2;++ix){
qkout.reset(ix);
fout.push_back(qkout);
qbout.reset(ix);
aout.push_back(qbout);
gluon.reset(2*ix);
gout.push_back(gluon);
}
for(unsigned int ix=0;ix<3;++ix){
zin.reset(ix);
vin.push_back(zin);
}
vector<Complex> diag(2,0.);
double total(0.);
for(unsigned int inhel1=0;inhel1<3;++inhel1) {
for(unsigned int outhel1=0;outhel1<2;++outhel1) {
for(unsigned int outhel2=0;outhel2<2;++outhel2) {
for(unsigned int outhel3=0;outhel3<2;++outhel3) {
SpinorBarWaveFunction off1 =
FFGVertex()->evaluate(scale_,3,partons[1],fout[outhel1],gout[outhel3]);
diag[0] = FFZVertex()->evaluate(scale_,aout[outhel2],off1,vin[inhel1]);
SpinorWaveFunction off2 =
FFGVertex()->evaluate(scale_,3,partons[2],aout[outhel2],gout[outhel3]);
diag[1] = FFZVertex()->evaluate(scale_,off2,fout[outhel1],vin[inhel1]);
// sum of diagrams
Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.));
// me2
total += norm(sum);
}
}
}
}
// divide out the coupling
total /= norm(FFGVertex()->norm());
// return the total
return total*UnitRemoval::InvE2;
}
void SMZFermionsPOWHEGDecayer::doinit() {
// cast the SM pointer to the Herwig SM pointer
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(!hwsm) throw InitException()
<< "Wrong type of StandardModel object in "
<< "SMZFermionsPOWHEGDecayer::doinit() "
<< "the Herwig version must be used."
<< Exception::runerror;
// cast the vertices
FFZVertex_ = hwsm->vertexFFZ();
FFZVertex_->init();
FFGVertex_ = hwsm->vertexFFG();
FFGVertex_->init();
SMZDecayer::doinit();
gluon_ = getParticleData(ParticleID::g);
}
bool SMZFermionsPOWHEGDecayer::getEvent(vector<PPtr> hardProcess) {
Energy particleMass = ZERO;
for(unsigned int ix=0;ix<hardProcess.size();++ix) {
if(hardProcess[ix]->id()==ParticleID::Z0) {
mZ_ = hardProcess[ix]->mass();
}
else {
particleMass = hardProcess[ix]->mass();
}
}
// reduced mass
mu_ = particleMass/mZ_;
mu2_ = sqr(mu_);
// scale
scale_ = sqr(mZ_);
// max pT
Energy pTmax = 0.5*sqrt(mz2_);
if(pTmax<pTmin_) return false;
// Define over valued y_max & y_min according to the associated pt_min cut.
- double ymax = acosh(pTmax/pTmin_);
+ //double ymax = acosh(pTmax/pTmin_);
+ double ymax=10.;
double ymin = -ymax;
// pt of the emmission
pT_ = pTmax;
// prefactor
double overEst = 4.;
double prefactor = overEst*alphaS()->overestimateValue()*CF_*
(ymax-ymin)/Constants::twopi;
// loop to generate the pt and rapidity
bool reject;
//arrays to hold the temporary probabilities whilst the for loop progresses
double probTemp[2][2]={{0.,0.},{0.,0.}};
probTemp[0][0]=probTemp[0][1]=probTemp[1][0]=probTemp[1][1]=0.;
double x1Solution[2][2] = {{0.,0.},{0.,0.}};
double x2Solution[2][2] = {{0.,0.},{0.,0.}};
double x3Solution[2] = {0.,0.};
Energy pT[2] = {pTmax,pTmax};
double yTemp[2] = {0.,0.};
double phi = 0.;
// do the competition
for(int i=0; i<2; i++) {
do {
//generation of phi
phi = UseRandom::rnd() * Constants::twopi;
// reject the emission
reject = true;
// generate pt
pT[i] *= pow(UseRandom::rnd(),1./prefactor);
Energy2 pT2 = sqr(pT[i]);
if(pT[i]<pTmin_) {
pT[i] = -GeV;
break;
}
// generate y
yTemp[i] = ymin + UseRandom::rnd()*(ymax-ymin);
//generate x3 & x1 from pT & y
double x1Plus = 1.;
double x1Minus = 2.*mu_;
x3Solution[i] = 2.*pT[i]*cosh(yTemp[i])/mZ_;
// prefactor
double weightPrefactor = 0.5/sqrt(1.-4.*mu2_)/overEst;
// calculate x1 & x2 solutions
Energy4 discrim2 = (sqr(x3Solution[i]*mZ_) - 4.*pT2)*
(mz2_*(x3Solution[i]-1.)*(4.*mu2_+x3Solution[i]-1.)-4.*mu2_*pT2);
//check discriminant2 is > 0
if( discrim2 < ZERO) continue;
Energy2 discriminant = sqrt(discrim2);
Energy2 fact1 = 3.*mz2_*x3Solution[i]-2.*mz2_+2.*pT2*x3Solution[i]
-4.*pT2-mz2_*sqr(x3Solution[i]);
Energy2 fact2 = 2.*mz2_*(x3Solution[i]-1.)-2.*pT2;
// two solns for x1
x1Solution[i][0] = (fact1 + discriminant)/fact2;
x1Solution[i][1] = (fact1 - discriminant)/fact2;
bool found = false;
for(unsigned int j=0;j<2;++j) {
x2Solution[i][0] = 2.-x3Solution[i]-x1Solution[i][0];
x2Solution[i][1] = 2.-x3Solution[i]-x1Solution[i][1];
// set limits on x2
double root = max(0.,sqr(x1Solution[i][j])-4.*mu2_);
root = sqrt(root);
double x2Plus = 1.-0.5*(1.-x1Solution[i][j])/(1.-x1Solution[i][j]+mu2_)
*(x1Solution[i][j]-2.*mu2_-root);
double x2Minus = 1.-0.5*(1.-x1Solution[i][j])/(1.-x1Solution[i][j]+mu2_)
*(x1Solution[i][j]-2.*mu2_+root);
if(x1Solution[i][j]>=x1Minus && x1Solution[i][j]<=x1Plus &&
x2Solution[i][j]>=x2Minus && x2Solution[i][j]<=x2Plus &&
checkZMomenta(x1Solution[i][j], x2Solution[i][j], x3Solution[i], yTemp[i], pT[i])) {
probTemp[i][j] = weightPrefactor*pT[i]*
- calculateJacobian(x1Solution[i][j], x2Solution[i][j], pT[i])*
+ abs(calculateJacobian(x1Solution[i][j], x2Solution[i][j], pT[i]))*
calculateRealEmission(x1Solution[i][j], x2Solution[i][j],
hardProcess, phi, false, i);
found = true;
}
else {
probTemp[i][j] = 0.;
}
}
if(!found) continue;
// alpha S piece
double wgt = (probTemp[i][0]+probTemp[i][1])*alphaS()->ratio(sqr(pT[i]));
// matrix element weight
reject = UseRandom::rnd()>wgt;
}
while(reject);
} //end of emitter for loop
// no emission
if(pT[0]<ZERO&&pT[1]<ZERO) return false;
//pick the spectator and x1 x2 values
double x1,x2,y;
//particle 1 emits, particle 2 spectates
unsigned int iemit=0;
if(pT[0]>pT[1]){
pT_ = pT[0];
y=yTemp[0];
if(probTemp[0][0]>UseRandom::rnd()*(probTemp[0][0]+probTemp[0][1])) {
x1 = x1Solution[0][0];
x2 = x2Solution[0][0];
}
else {
x1 = x1Solution[0][1];
x2 = x2Solution[0][1];
}
}
// particle 2 emits, particle 1 spectates
else {
iemit=1;
pT_ = pT[1];
y=yTemp[1];
if(probTemp[1][0]>UseRandom::rnd()*(probTemp[1][0]+probTemp[1][1])) {
x1 = x1Solution[1][0];
x2 = x2Solution[1][0];
}
else {
x1 = x1Solution[1][1];
x2 = x2Solution[1][1];
}
}
// find spectator
unsigned int ispect = iemit == 0 ? 1 : 0;
// Find the boost from the lab to the c.o.m with the spectator
// along the -z axis, and then invert it.
LorentzRotation eventFrame( ( quark_[0] + quark_[1] ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*quark_[ispect];
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() - Constants::pi );
eventFrame.invert();
// spectator
quark_[ispect].setT( 0.5*x2*mZ_ );
quark_[ispect].setX( ZERO );
quark_[ispect].setY( ZERO );
quark_[ispect].setZ( -sqrt(0.25*mz2_*x2*x2-mz2_*mu2_) );
// gluon
gauge_.setT( pT_*cosh(y) );
gauge_.setX( pT_*cos(phi) );
gauge_.setY( pT_*sin(phi) );
gauge_.setZ( pT_*sinh(y) );
gauge_.setMass(ZERO);
// emitter reconstructed from gluon & spectator
quark_[iemit] = - gauge_ - quark_[ispect];
quark_[iemit].setT( 0.5*mZ_*x1 );
// boost constructed vectors into the event frame
quark_[0] = eventFrame * quark_[0];
quark_[1] = eventFrame * quark_[1];
gauge_ = eventFrame * gauge_;
// need to reset masses because for whatever reason the boost
// touches the mass component of the five-vector and can make
// zero mass objects acquire a floating point negative mass(!).
gauge_.setMass( ZERO );
quark_[iemit] .setMass(partons_[iemit ]->mass());
quark_[ispect].setMass(partons_[ispect]->mass());
return true;
}
InvEnergy SMZFermionsPOWHEGDecayer::calculateJacobian(double x1, double x2, Energy pT) const{
double xPerp = abs(2.*pT/mZ_);
Energy jac = mZ_*fabs((x1*x2-2.*mu2_*(x1+x2)+sqr(x2)-x2)/xPerp/pow(sqr(x2)-4.*mu2_,1.5));
return 1./jac; //jacobian as defined is dptdy=jac*dx1dx2, therefore we have to divide by it
}
bool SMZFermionsPOWHEGDecayer::checkZMomenta(double x1, double x2, double x3, double y, Energy pT) const {
double xPerp2 = 4.*pT*pT/mZ_/mZ_;
static double tolerance = 1e-6;
bool isMomentaReconstructed = false;
if(pT*sinh(y)>ZERO) {
if(abs(-sqrt(sqr(x2)-4.*mu2_)+sqrt(sqr(x3)-xPerp2) + sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance ||
abs(-sqrt(sqr(x2)-4.*mu2_)+sqrt(sqr(x3)-xPerp2) - sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance) isMomentaReconstructed=true;
}
else if(pT*sinh(y) < ZERO){
if(abs(-sqrt(sqr(x2)-4.*mu2_)-sqrt(sqr(x3)-xPerp2) + sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance ||
abs(-sqrt(sqr(x2)-4.*mu2_)-sqrt(sqr(x3)-xPerp2) - sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance) isMomentaReconstructed=true;
}
else
if(abs(-sqrt(sqr(x2)-4.*mu2_)+ sqrt(sqr(x1)-xPerp2 - 4.*mu2_)) <= tolerance) isMomentaReconstructed=true;
return isMomentaReconstructed;
}
double SMZFermionsPOWHEGDecayer::calculateRealEmission(double x1, double x2,
vector<PPtr> hardProcess,
double phi,
bool subtract) const {
// make partons data object for meRatio
vector<cPDPtr> partons (3);
for(int ix=0; ix<3; ++ix)
partons[ix] = hardProcess[ix]->dataPtr();
partons.push_back(gluon_);
// calculate x3
double x3 = 2.-x1-x2;
double xT = sqrt(max(0.,sqr(x3) -0.25*sqr(sqr(x2)+sqr(x3)-sqr(x1))/(sqr(x2)-4.*mu2_)));
// calculate the momenta
Energy M = mZ_;
Lorentz5Momentum pspect(ZERO,ZERO,-0.5*M*sqrt(max(sqr(x2)-4.*mu2_,0.)),0.5*M*x2,M*mu_);
Lorentz5Momentum pemit (-0.5*M*xT*cos(phi),-0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x1)-sqr(xT)-4.*mu2_,0.)),0.5*M*x1,M*mu_);
Lorentz5Momentum pgluon(0.5*M*xT*cos(phi), 0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x3)-sqr(xT),0.)),0.5*M*x3,ZERO);
if(abs(pspect.z()+pemit.z()-pgluon.z())/M<1e-6)
pgluon.setZ(-pgluon.z());
else if(abs(pspect.z()-pemit.z()+pgluon.z())/M<1e-6)
pemit .setZ(- pemit.z());
// loop over the possible emitting partons
double realwgt(0.);
for(unsigned int iemit=0;iemit<2;++iemit) {
// boost and rotate momenta
LorentzRotation eventFrame( ( hardProcess[1]->momentum() +
hardProcess[2]->momentum() ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*hardProcess[iemit+1]->momentum();
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() );
eventFrame.invert();
vector<Lorentz5Momentum> momenta(3);
momenta[0] = hardProcess[0]->momentum();
if(iemit==0) {
momenta[2] = eventFrame*pspect;
momenta[1] = eventFrame*pemit ;
}
else {
momenta[1] = eventFrame*pspect;
momenta[2] = eventFrame*pemit ;
}
momenta.push_back(eventFrame*pgluon);
// calculate the weight
if(1.-x1>1e-5 && 1.-x2>1e-5)
realwgt += meRatio(partons,momenta,iemit,subtract);
}
// total real emission contribution
return realwgt;
}
double SMZFermionsPOWHEGDecayer::calculateRealEmission(double x1, double x2,
vector<PPtr> hardProcess,
double phi,
bool subtract,
int emitter) const {
// make partons data object for meRatio
vector<cPDPtr> partons (3);
for(int ix=0; ix<3; ++ix)
partons[ix] = hardProcess[ix]->dataPtr();
partons.push_back(gluon_);
// calculate x3
double x3 = 2.-x1-x2;
double xT = sqrt(max(0.,sqr(x3) -0.25*sqr(sqr(x2)+sqr(x3)-sqr(x1))/(sqr(x2)-4.*mu2_)));
// calculate the momenta
Energy M = mZ_;
Lorentz5Momentum pspect(ZERO,ZERO,-0.5*M*sqrt(max(sqr(x2)-4.*mu2_,0.)),0.5*M*x2,M*mu_);
Lorentz5Momentum pemit (-0.5*M*xT*cos(phi),-0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x1)-sqr(xT)-4.*mu2_,0.)),0.5*M*x1,M*mu_);
Lorentz5Momentum pgluon( 0.5*M*xT*cos(phi), 0.5*M*xT*sin(phi),
0.5*M*sqrt(max(sqr(x3)-sqr(xT),0.)),0.5*M*x3,ZERO);
if(abs(pspect.z()+pemit.z()-pgluon.z())/M<1e-6)
pgluon.setZ(-pgluon.z());
else if(abs(pspect.z()-pemit.z()+pgluon.z())/M<1e-6)
pemit .setZ(- pemit.z());
// loop over the possible emitting partons
double realwgt(0.);
// boost and rotate momenta
LorentzRotation eventFrame( ( hardProcess[1]->momentum() +
hardProcess[2]->momentum() ).findBoostToCM() );
Lorentz5Momentum spectator = eventFrame*hardProcess[emitter+1]->momentum();
eventFrame.rotateZ( -spectator.phi() );
eventFrame.rotateY( -spectator.theta() );
eventFrame.invert();
vector<Lorentz5Momentum> momenta(3);
momenta[0] = hardProcess[0]->momentum();
if(emitter==0) {
momenta[2] = eventFrame*pspect;
momenta[1] = eventFrame*pemit ;
}
else {
momenta[1] = eventFrame*pspect;
momenta[2] = eventFrame*pemit ;
}
momenta.push_back(eventFrame*pgluon);
// calculate the weight
if(1.-x1>1e-5 && 1.-x2>1e-5)
realwgt += meRatio(partons,momenta,emitter,subtract);
// total real emission contribution
return realwgt;
}
diff --git a/Tests/Inputs/GammaP.common b/Tests/Inputs/GammaP.common
--- a/Tests/Inputs/GammaP.common
+++ b/Tests/Inputs/GammaP.common
@@ -1,39 +1,39 @@
##################################################
# Example generator based for gamma hadron collisions
##################################################
cd /Herwig
create Herwig::AlphaEM AlphaEM2
set Model:EW/RunningAlphaEM AlphaEM2
create Herwig::O2AlphaS AlphaS2
set Model:QCD/RunningAlphaS AlphaS2
# Create GammaHadronHandler
cd /Herwig/EventHandlers
set Luminosity:Energy 1000.
set EventHandler:BeamA /Herwig/Particles/gamma
set EventHandler:BeamB /Herwig/Particles/p+
set EventHandler:Sampler /Herwig/ACDCSampler
-set /Herwig/Partons/Extractor:FlatSHatY 1
+set /Herwig/Partons/PPExtractor:FlatSHatY 1
# the cuts
cd /Herwig/Cuts
set Cuts:X1Min 1.0e-5
set Cuts:X2Min 1.0e-5
set Cuts:MHatMin 0.*GeV
insert Cuts:OneCuts[0] JetKtCut
set JetKtCut:MinKT 20.
##################################################
# Technical parameters for this run
##################################################
cd /Herwig/Generators
set EventGenerator:NumberOfEvents 100000000
set EventGenerator:RandomNumberGenerator:Seed 31122001
set EventGenerator:DebugLevel 1
set EventGenerator:PrintEvent 0
set EventGenerator:MaxErrors 10000
set EventGenerator:EventHandler:StatLevel Full
set EventGenerator:EventHandler:CascadeHandler:MPIHandler NULL
##################################################
# LEP physics parameters (override defaults)
##################################################
set EventGenerator:EventHandler:LuminosityFunction:Energy 1000.
diff --git a/Tests/Inputs/LEP-BB.in b/Tests/Inputs/LEP-BB.in
deleted file mode 100644
--- a/Tests/Inputs/LEP-BB.in
+++ /dev/null
@@ -1,21 +0,0 @@
-# -*- ThePEG-repository -*-
-read snippets/EECollider.in
-cd /Herwig/MatrixElements
-set MEee2gZ2qq:MinimumFlavour 4
-set MEee2gZ2qq:MaximumFlavour 4
-insert SubProcess:MatrixElements 0 MEee2gZ2qq
-
-read LEP.common
-
-cd /Herwig/Generators
-
-set EventGenerator:EventHandler:LuminosityFunction:Energy 10.53
-set EventGenerator:EventHandler:CascadeHandler:MPIHandler NULL
-
-cd /Herwig/Generators
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/BELLECharm
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/CLEOCharm
-
-set /Herwig/Cuts/Cuts:MHatMin 10.5299
-
-saverun LEP-BB EventGenerator
diff --git a/Tests/Inputs/LEP-Leptons.in b/Tests/Inputs/LEP-Leptons.in
--- a/Tests/Inputs/LEP-Leptons.in
+++ b/Tests/Inputs/LEP-Leptons.in
@@ -1,39 +1,41 @@
+# -*- ThePEG-repository -*-
+read snippets/EECollider.in
cd /Herwig/MatrixElements
insert SubProcess:MatrixElements[0] MEee2gZ2ll
set /Herwig/Particles/tau+:Stable 1
set /Herwig/Particles/tau-:Stable 1
read LEP.common
cd /Herwig/Generators
set EventGenerator:EventHandler:CascadeHandler NULL
set EventGenerator:EventHandler:HadronizationHandler NULL
set EventGenerator:EventHandler:DecayHandler NULL
set EventGenerator:EventHandler:LuminosityFunction:Energy 500.0*GeV
set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
create Herwig::FermionTest LeptonsTest LeptonTest.so
insert EventGenerator:AnalysisHandlers 0 LeptonsTest
set /Herwig/Analysis/Basics:CheckQuark No
# parameters to make the same as fortran
cd /Herwig
create Herwig::O2AlphaS AlphaS2
set Model:QCD/RunningAlphaS AlphaS2
set Model:EW/CKM:theta_12 0.22274457
set Model:EW/CKM:theta_13 0.
set Model:EW/CKM:theta_23 0.
set Model:EW/CKM:delta 0.
set Model:EW/Sin2ThetaW 0.22254916
create Herwig::AlphaEM AlphaEM2
set Model:EW/RunningAlphaEM AlphaEM2
set /Herwig/ACDCSampler:Ntry 100000
cd /Herwig/Generators
saverun LEP-Leptons EventGenerator
diff --git a/Tests/Inputs/LEP-Powheg.in b/Tests/Inputs/LEP-Powheg.in
deleted file mode 100644
--- a/Tests/Inputs/LEP-Powheg.in
+++ /dev/null
@@ -1,26 +0,0 @@
-# -*- ThePEG-repository -*-
-read snippets/EECollider.in
-cd /Herwig/MatrixElements
-create Herwig::O2AlphaS O2AlphaS
-set /Herwig/Generators/EventGenerator:StandardModelParameters:QCD/RunningAlphaS O2AlphaS
-set /Herwig/Shower/ShowerHandler:HardEmission POWHEG
-set /Herwig/Shower/AlphaQCD:AlphaMZ 0.118
-insert SubProcess:MatrixElements 0 PowhegMEee2gZ2qq
-
-read LEP.common
-
-cd /Herwig/Generators
-
-set EventGenerator:EventHandler:CascadeHandler:MPIHandler NULL
-
-set EventGenerator:EventHandler:LuminosityFunction:Energy 91.2
-
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPMultiplicity
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/BMultiplicity
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/BFrag
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/Shapes
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPIdent
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPFourJet
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPJet
-
-saverun LEP-Powheg EventGenerator
diff --git a/Tests/Inputs/LEP-default.in b/Tests/Inputs/LEP-default.in
deleted file mode 100644
--- a/Tests/Inputs/LEP-default.in
+++ /dev/null
@@ -1,22 +0,0 @@
-# -*- ThePEG-repository -*-
-read snippets/EECollider.in
-cd /Herwig/MatrixElements
-insert SubProcess:MatrixElements 0 MEee2gZ2qq
-
-read LEP.common
-
-cd /Herwig/Generators
-
-set EventGenerator:EventHandler:CascadeHandler:MPIHandler NULL
-
-set EventGenerator:EventHandler:LuminosityFunction:Energy 91.2
-
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPMultiplicity
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/BMultiplicity
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/BFrag
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/Shapes
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPIdent
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPFourJet
-insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/LEPJet
-
-saverun LEP-default EventGenerator
diff --git a/Tests/Makefile.am b/Tests/Makefile.am
--- a/Tests/Makefile.am
+++ b/Tests/Makefile.am
@@ -1,365 +1,365 @@
AM_LDFLAGS += -module -avoid-version -rpath /dummy/path/not/used
EXTRA_DIST = Inputs python Rivet
EXTRA_LTLIBRARIES = LeptonTest.la GammaTest.la HadronTest.la DISTest.la
if WANT_LIBFASTJET
EXTRA_LTLIBRARIES += HadronJetTest.la LeptonJetTest.la
HadronJetTest_la_SOURCES = \
Hadron/VHTest.h Hadron/VHTest.cc\
Hadron/VTest.h Hadron/VTest.cc\
Hadron/HTest.h Hadron/HTest.cc
HadronJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
HadronJetTest_la_LIBADD = $(FASTJETLIBS)
LeptonJetTest_la_SOURCES = \
Lepton/TopDecay.h Lepton/TopDecay.cc
LeptonJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
LeptonJetTest_la_LIBADD = $(FASTJETLIBS)
endif
LeptonTest_la_SOURCES = \
Lepton/VVTest.h Lepton/VVTest.cc \
Lepton/VBFTest.h Lepton/VBFTest.cc \
Lepton/VHTest.h Lepton/VHTest.cc \
Lepton/FermionTest.h Lepton/FermionTest.cc
GammaTest_la_SOURCES = \
Gamma/GammaMETest.h Gamma/GammaMETest.cc \
Gamma/GammaPMETest.h Gamma/GammaPMETest.cc
DISTest_la_SOURCES = \
DIS/DISTest.h DIS/DISTest.cc
HadronTest_la_SOURCES = \
Hadron/HadronVVTest.h Hadron/HadronVVTest.cc\
Hadron/HadronVBFTest.h Hadron/HadronVBFTest.cc\
Hadron/WHTest.h Hadron/WHTest.cc\
Hadron/ZHTest.h Hadron/ZHTest.cc\
Hadron/VGammaTest.h Hadron/VGammaTest.cc\
Hadron/ZJetTest.h Hadron/ZJetTest.cc\
Hadron/WJetTest.h Hadron/WJetTest.cc\
Hadron/QQHTest.h Hadron/QQHTest.cc
REPO = $(top_builddir)/src/HerwigDefaults.rpo
HERWIG = $(top_builddir)/src/Herwig
HWREAD = $(HERWIG) read --repo $(REPO) -L $(builddir)/.libs -i $(top_builddir)/src
HWRUN = $(HERWIG) run
tests : tests-LEP tests-DIS tests-LHC tests-Gamma
if WANT_LIBFASTJET
-tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons \
- test-LEP-default test-LEP-Powheg test-LEP-TopDecay
+tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-Quarks test-LEP-Leptons \
+ test-LEP-TopDecay
else
tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons
endif
tests-DIS : test-DIS-Charged test-DIS-Neutral
if WANT_LIBFASTJET
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top test-LHC-Bottom \
test-LHC-WHJet test-LHC-ZHJet test-LHC-HJet test-LHC-ZShower test-LHC-WShower\
test-LHC-WHJet-Powheg test-LHC-ZHJet-Powheg test-LHC-HJet-Powheg \
test-LHC-ZShower-Powheg test-LHC-WShower-Powheg
else
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top
endif
tests-Gamma : test-Gamma-FF test-Gamma-WW test-Gamma-P
if WANT_LIBFASTJET
test-LEP-% : Inputs/LEP-%.in LeptonTest.la LeptonJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LEP-% : Inputs/LEP-%.in LeptonTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
Rivet-LHC-Matchbox-% : Rivet/LHC-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-TVT-Matchbox-% : Rivet/TVT-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-TVT-Dipole-% : Rivet/TVT-Dipole-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-LHC-Dipole-% : Rivet/LHC-Dipole-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet/LEP-%.in :
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/DIS-%.in :
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/BFactory-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/TVT-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/LHC-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/Star-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/SppS-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/ISR-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet-LEP-Matchbox-% : Rivet/LEP-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-LEP-Dipole-% : Rivet/LEP-Dipole-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-BFactory-Matchbox-% : Rivet/BFactory-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-LEP-% : Rivet/LEP-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-BFactory-% : Rivet/BFactory-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-TVT-% : Rivet/TVT-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-DIS-% : Rivet/DIS-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-LHC-% : Rivet/LHC-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-Star-% : Rivet/Star-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-SppS-% : Rivet/SppS-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-ISR-% : Rivet/ISR-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-inputfiles: $(shell echo Rivet/LEP{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{9.4,12,13,17,27.6,29,30.2,30.7,30.75,30,31.3,34.8,43.6,50,52,55,56,57,60.8,60,61.4,10,12.8,22,26.8,35,44,48.0,91,93.0,130,133,136,161,172,177,183,189,192,196,197,200,202,206,91-nopi}.in) \
$(shell echo Rivet/LEP{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Powheg,-Matchbox-Powheg}-14.in) \
$(shell echo Rivet/LEP{,-Dipole}-{10.5,11.96,12.8,13.96,16.86,21.84,26.8,28.48,35.44,48.0,97.0}-gg.in) \
$(shell echo Rivet/BFactory{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg}-{10.52,10.52-sym,10.54,10.45}.in) \
$(shell echo Rivet/BFactory{,-Dipole}-{Upsilon,Upsilon2,Upsilon4,Tau,10.58-res}.in) \
$(shell echo Rivet/DIS{,-NoME,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{e--LowQ2,e+-LowQ2,e+-HighQ2}.in) \
$(shell echo Rivet/TVT{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{Run-I-Z,Run-I-W,Run-I-WZ,Run-II-Z-e,Run-II-Z-{,LowMass-,HighMass-}mu,Run-II-W}.in) \
$(shell echo Rivet/TVT{,-Dipole}-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet,PromptPhoton}.in) \
$(shell echo Rivet/TVT-Powheg-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.in) \
$(shell echo Rivet/TVT{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{Run-II-Jets-{0..11},Run-I-Jets-{1..8}}.in ) \
$(shell echo Rivet/TVT{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{630-Jets-{1..3},300-Jets-1,900-Jets-1}.in ) \
$(shell echo Rivet/TVT{,-Dipole}-{Run-I,Run-II,300,630,900}-UE.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-DiJets-{1..7}-{A,B,C}.in ) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{7,8,13}-Jets-{0..10}.in ) \
$(shell echo Rivet/LHC{,-Dipole}-{900,2360,2760,7,8,13}-UE.in ) \
$(shell echo Rivet/LHC{,-Dipole}-{900,7,13}-UE-Long.in ) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Charm-{1..5}.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Bottom-{0..8}.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Top-{L,SL}.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{7,8,13}-Top-All.in) \
$(shell echo Rivet/Star{,-Dipole}-{UE,Jets-{1..4}}.in ) \
$(shell echo Rivet/SppS{,-Dipole}-{200,500,900,546}-UE.in ) \
$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{W-{e,mu},13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},Z-{e,mu,mu-SOPHTY},Z-LowMass-{e,mu},Z-MedMass-e,WZ,WW-{emu,ll},ZZ-{ll,lv},8-ZZ-lv,8-WW-ll,W-Z-{e,mu}}.in) \
$(shell echo Rivet/LHC{,-Dipole}-7-{W,Z}Gamma-{e,mu}.in) \
$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{7-W-Jet-{1..3}-e,7-Z-Jet-{0..3}-e,7-Z-Jet-0-mu}.in) \
$(shell echo Rivet/LHC{-Matchbox,-Matchbox-Powheg,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{Z-b,Z-bb,W-b,8-Z-jj}.in) \
$(shell echo Rivet/LHC{,-Dipole}-{7,8}-PromptPhoton-{1..4}.in) Rivet/LHC-GammaGamma-7.in \
$(shell echo Rivet/LHC{,-Powheg}-{7,8}-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.in) \
$(shell echo Rivet/LHC{,-Powheg,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{ggH,VBF,WH,ZH}.in) \
$(shell echo Rivet/LHC{,-Powheg,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-8-{{ggH,VBF,WH,ZH}{,-GammaGamma},ggH-WW}.in) \
$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-ggHJet.in)
# $(shell echo Rivet/ISR-{30,44,53,62}-UE.in ) $(shell echo Rivet/SppS-{53,63}-UE.in )
Rivet-LEP: $(shell echo Rivet-LEP{,-Powheg,-Matchbox,-Dipole}-{9.4,12,13,17,27.6,29,30.2,30.7,30.75,30,31.3,34.8,43.6,50,52,55,56,57,60.8,60,61.4,10,12.8,22,26.8,35,44,48.0,91,93.0,130,133,136,161,172,177,183,189,192,196,197,200,202,206,91-nopi}) \
$(shell echo Rivet-LEP-{10.5,11.96,12.8,13.96,16.86,21.84,26.8,28.48,35.44,48.0,97.0}-gg) \
$(shell echo Rivet-LEP{,-Powheg,-Matchbox-Powheg}-14.in)
rm -rf Rivet-LEP
python/merge-LEP --with-gg LEP
python/merge-LEP LEP-Powheg
python/merge-LEP LEP-Matchbox
python/merge-LEP LEP-Dipole
rivet-mkhtml -o Rivet-LEP LEP.yoda:Hw LEP-Powheg.yoda:Hw-Powheg LEP-Matchbox.yoda:Hw-Matchbox LEP-Dipole.yoda:Hw-Dipole
Rivet-BFactory: $(shell echo Rivet-BFactory{,-Powheg,-Matchbox,-Dipole}-{10.52,10.52-sym,10.54,10.45}) \
$(shell echo Rivet-BFactory-{Upsilon,Upsilon2,Upsilon4,Tau,10.58-res,10.58})
rm -rf Rivet-BFactory
python/merge-BFactory BFactory
python/merge-BFactory BFactory-Powheg
python/merge-BFactory BFactory-Matchbox
python/merge-BFactory BFactory-Dipole
rivet-mkhtml -o Rivet-BFactory BFactory.yoda:Hw BFactory-Powheg.yoda:Hw-Powheg BFactory-Matchbox.yoda:Hw-Matchbox BFactory-Dipole.yoda:Hw-Dipole
Rivet-DIS: $(shell echo Rivet-DIS{,-NoME,-Powheg,-Matchbox,-Dipole}-{e--LowQ2,e+-LowQ2,e+-HighQ2})
rm -rf Rivet-DIS
python/merge-DIS DIS
python/merge-DIS DIS-Powheg
python/merge-DIS DIS-NoME
python/merge-DIS DIS-Matchbox
python/merge-DIS DIS-Dipole
rivet-mkhtml -o Rivet-DIS DIS.yoda:Hw DIS-Powheg.yoda:Hw-Powheg DIS-NoME.yoda:Hw-NoME DIS-Matchbox.yoda:Hw-Matchbox DIS-Dipole.yoda:Hw-Dipole
Rivet-TVT-WZ: $(shell echo Rivet-TVT{,-Powheg,-Matchbox,-Dipole}-{Run-I-Z,Run-I-W,Run-I-WZ,Run-II-Z-{e,{,LowMass-,HighMass-}mu},Run-II-W})
rm -rf Rivet-TVT-WZ
python/merge-TVT-EW TVT-Run-II-W.yoda TVT-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Run-I-{W,Z,WZ}.yoda -o TVT-WZ.yoda
python/merge-TVT-EW TVT-Powheg-Run-II-W.yoda TVT-Powheg-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Powheg-Run-I-{W,Z,WZ}.yoda -o TVT-Powheg-WZ.yoda
python/merge-TVT-EW TVT-Matchbox-Run-II-W.yoda TVT-Matchbox-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Matchbox-Run-I-{W,Z,WZ}.yoda -o TVT-Matchbox-WZ.yoda
python/merge-TVT-EW TVT-Dipole-Run-II-W.yoda TVT-Dipole-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Dipole-Run-I-{W,Z,WZ}.yoda -o TVT-Dipole-WZ.yoda
rivet-mkhtml -o Rivet-TVT-WZ TVT-WZ.yoda:Hw TVT-Powheg-WZ.yoda:Hw-Powheg TVT-Matchbox-WZ.yoda:Hw-Matchbox TVT-Dipole-WZ.yoda:Hw-Dipole
Rivet-TVT-Photon: $(shell echo Rivet-TVT-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet,PromptPhoton}) \
$(shell echo Rivet-TVT-Powheg-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet})
rm -rf Rivet-TVT-Photon
python/merge-TVT-Photon TVT -o TVT-Photon.yoda
python/merge-TVT-Photon TVT-Powheg -o TVT-Powheg-Photon.yoda
rivet-mkhtml -o Rivet-TVT-Photon TVT-Photon.yoda:Hw TVT-Powheg-Photon.yoda:Hw-Powheg
Rivet-TVT-Jets: $(shell echo Rivet-TVT-{Run-II-Jets-{0..11},Run-I-Jets-{1..8}} ) \
$(shell echo Rivet-TVT-{630-Jets-{1..3},300-Jets-1,900-Jets-1} ) \
$(shell echo Rivet-TVT-{Run-I,Run-II,300,630,900}-UE)
rm -rf Rivet-TVT-Jets
python/merge-TVT-Jets TVT
rivet-mkhtml -o Rivet-TVT-Jets TVT-Jets.yoda:Hw
#python/merge-TVT-Energy TVT
#rivet-merge-CDF_2012_NOTE10874 TVT-300-Energy.yoda TVT-900-Energy.yoda TVT-1960-Energy.yoda
#flat2yoda RatioPlots.dat -o TVT-RatioPlots.yoda
Rivet-LHC-Jets: $(shell echo Rivet-LHC-7-DiJets-{1..7}-{A,B,C} ) \
$(shell echo Rivet-LHC-{7,8,13}-Jets-{0..10} ) \
$(shell echo Rivet-LHC-{900,2360,2760,7,8,13}-UE ) \
$(shell echo Rivet-LHC-{900,7,13}-UE-Long ) \
$(shell echo Rivet-LHC-7-Charm-{1..5}) \
$(shell echo Rivet-LHC-7-Bottom-{0..8}) \
$(shell echo Rivet-LHC-7-Top-{L,SL})\
$(shell echo Rivet-LHC-{7,8,13}-Top-All)
rm -rf Rivet-LHC-Jets
python/merge-LHC-Jets LHC
rivet-mkhtml -o Rivet-LHC-Jets LHC-Jets.yoda:Hw
Rivet-Star: $(shell echo Rivet-Star-{UE,Jets-{1..4}} )
rm -rf Rivet-Star
python/merge-Star Star
rivet-mkhtml -o Rivet-Star Star.yoda
Rivet-SppS: $(shell echo Rivet-ISR-{30,44,53,62}-UE ) \
$(shell echo Rivet-SppS-{53,63,200,500,900,546}-UE )
rm -rf Rivet-SppS
python/merge-SppS SppS
rivet-mkhtml -o Rivet-SppS SppS.yoda
Rivet-LHC-EW: $(shell echo Rivet-LHC{,-Matchbox,-Powheg,-Dipole}-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-SOPHTY},Z-LowMass-{e,mu},Z-MedMass-e,WZ,WW-{emu,ll},ZZ-{ll,lv},8-ZZ-lv,,8-WW-llW-Z-{e,mu}}) \
$(shell echo Rivet-LHC{,-Matchbox,-Dipole}-{7-W-Jet-{1..3}-e,7-Z-Jet-{0..3}-e,7-Z-Jet-0-mu}) \
$(shell echo Rivet-LHC{-Matchbox,-Dipole}-{Z-b,Z-bb,W-b,8-Z-jj}) \
$(shell echo Rivet-LHC-7-{W,Z}Gamma-{e,mu}) \
rm -rf Rivet-LHC-EW;
python/merge-LHC-EW LHC-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-Short},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda LHC-7-{W,Z}-Jet-{1,2,3}-e.yoda LHC-7-{W,Z}Gamma-{e,mu}.yoda -o LHC-EW.yoda;
python/merge-LHC-EW LHC-Matchbox-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-Short},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda LHC-Matchbox-7-{W,Z}-Jet-{1,2,3}-e.yoda -o LHC-Matchbox-EW.yoda;
python/merge-LHC-EW LHC-Dipole-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-Short},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda LHC-Dipole-7-{W,Z}-Jet-{1,2,3}-e.yoda -o LHC-Dipole-EW.yoda;
python/merge-LHC-EW LHC-Powheg-{W-{e,mu},Z-{e,mu},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda -o LHC-Powheg-EW.yoda;
rivet-mkhtml -o Rivet-LHC-EW LHC-EW.yoda:Hw LHC-Powheg-EW.yoda:Hw-Powheg LHC-Matchbox-EW.yoda:Hw-Matchbox LHC-Matchbox-Z-b.yoda:Hw-Matchbox-Zb \
LHC-Matchbox-Z-bb.yoda:Hw-Matchbox-Zbb LHC-Matchbox-W-b.yoda:Hw-Matchbox-W-bb LHC-Dipole-EW.yoda:Hw-Dipole \
LHC-Dipole-Z-b.yoda:Hw-Dipole-Zb LHC-Dipole-Z-bb.yoda:Hw-Dipole-Zbb LHC-Dipole-W-b.yoda:Hw-Dipole-W-bb \
LHC-Z-mu-SOPHTY.yoda:Hw LHC-Powheg-Z-mu-SOPHTY.yoda:Hw-Powheg LHC-Matchbox-Z-mu-SOPHTY.yoda:Hw-Matchbox
Rivet-LHC-Photon: $(shell echo Rivet-LHC-{7,8}-PromptPhoton-{1..4}) Rivet-LHC-GammaGamma-7 \
$(shell echo Rivet-LHC{,-Powheg}-{7,8}-{DiPhoton-GammaGamma,DiPhoton-GammaJet})
rm -rf Rivet-LHC-Photon
python/merge-LHC-Photon LHC -o LHC-Photon.yoda
python/merge-LHC-Photon LHC-Powheg -o LHC-Powheg-Photon.yoda
rivet-mkhtml -o Rivet-LHC-Photon LHC-Photon.yoda:Hw LHC-Powheg-Photon.yoda:Hw-Powheg
Rivet-LHC-Higgs: $(shell echo Rivet-LHC{,-Powheg}-{ggH,VBF,WH,ZH})\
$(shell echo Rivet-LHC{,-Powheg}-8-{{ggH,VBF,WH,ZH}{,-GammaGamma},ggH-WW}) Rivet-LHC-ggHJet
rm -rf Rivet-LHC-Higgs
rivet-mkhtml -o Rivet-LHC-Higgs LHC-Powheg-ggH.yoda:gg-Powheg LHC-ggH.yoda:gg LHC-ggHJet.yoda:HJet \
LHC-Powheg-VBF.yoda:VBF-Powheg LHC-VBF.yoda:VBF LHC-WH.yoda:WH LHC-ZH.yoda:ZH \
LHC-Powheg-WH.yoda:WH-Powheg LHC-Powheg-ZH.yoda:ZH-Powheg
tests-Rivet : Rivet-LEP Rivet-BFactory Rivet-DIS Rivet-TVT-WZ Rivet-TVT-Photon Rivet-TVT-Jets Rivet-LHC-Jets Rivet-Star Rivet-SppS Rivet-LHC-EW Rivet-LHC-Photon Rivet-LHC-Higgs
test-Gamma-% : Inputs/Gamma-%.in GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
test-DIS-% : Inputs/DIS-%.in DISTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
if WANT_LIBFASTJET
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la HadronJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
clean-local:
rm -f *.out *.log *.tex *.top *.run *.dump *.mult *.Bmult *.yoda