Page MenuHomeHEPForge
Files F7879446

No OneTemporary
Actions

View Options

diff --git a/Decay/DecayPhaseSpaceChannel.cc b/Decay/DecayPhaseSpaceChannel.cc
--- a/Decay/DecayPhaseSpaceChannel.cc
+++ b/Decay/DecayPhaseSpaceChannel.cc
@@ -1,587 +1,587 @@
// -*- C++ -*-
//
// DecayPhaseSpaceChannel.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 DecayPhaseSpaceChannel class.
//
// Author: Peter Richardson
//
#include "DecayPhaseSpaceChannel.h"
#include "DecayPhaseSpaceMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include <ThePEG/Repository/CurrentGenerator.h>
using namespace Herwig;
DecayPhaseSpaceChannel::DecayPhaseSpaceChannel(tcDecayPhaseSpaceModePtr in)
: _mode(in) {}
void DecayPhaseSpaceChannel::persistentOutput(PersistentOStream & os) const {
os << _intpart << _jactype << ounit(_intmass,GeV) << ounit(_intwidth,GeV)
<< ounit(_intmass2,GeV2) << ounit(_intmwidth,GeV2) << _intpower
<< _intdau1 << _intdau2 << _intext << _mode;
}
void DecayPhaseSpaceChannel::persistentInput(PersistentIStream & is, int) {
is >> _intpart >> _jactype >> iunit(_intmass,GeV) >> iunit(_intwidth,GeV)
>> iunit(_intmass2,GeV2) >> iunit(_intmwidth,GeV2) >> _intpower
>> _intdau1 >> _intdau2 >> _intext >> _mode;
}
ClassDescription<DecayPhaseSpaceChannel>
DecayPhaseSpaceChannel::initDecayPhaseSpaceChannel;
// Definition of the static class description member.
void DecayPhaseSpaceChannel::Init() {
static RefVector<DecayPhaseSpaceChannel,ParticleData> interfaceIntermediateParticles
("IntermediateParticles",
"The intermediate particles in the decay chain.",
&DecayPhaseSpaceChannel::_intpart, 0, false, false, true, false);
static ParVector<DecayPhaseSpaceChannel,int> interfacejactype
("Jacobian",
"The type of Jacobian to use for the intermediate particle",
&DecayPhaseSpaceChannel::_jactype,
0, 0, 0, 0, 1, false, false, true);
static ParVector<DecayPhaseSpaceChannel,double> interfaceIntermediatePower
("IntermediatePower",
"The power to use in the Jacobian",
&DecayPhaseSpaceChannel::_intpower,
0, 0, 0, -10, 10, false, false, true);
static ParVector<DecayPhaseSpaceChannel,int> interfaceIntermediateDau1
("IntermediateDaughter1",
"First Daughter of the intermediate",
&DecayPhaseSpaceChannel::_intdau1,
0, 0, 0, -10, 10, false, false, true);
static ParVector<DecayPhaseSpaceChannel,int> interfaceIntermediateDau2
("IntermediateDaughter2",
"Second Daughter of the intermediate",
&DecayPhaseSpaceChannel::_intdau2,
0, 0, 0, -10, 10, false, false, true);
static ClassDocumentation<DecayPhaseSpaceChannel> documentation
("The DecayPhaseSpaceChannel class defines a channel"
" for the multichannel integration of the phase space for a decay.");
}
// generate the momenta of the external particles
vector<Lorentz5Momentum>
DecayPhaseSpaceChannel::generateMomenta(const Lorentz5Momentum & pin,
const vector<Energy> & massext) {
// integers for loops
unsigned int ix,iy,idau[2],iz;
// storage of the momenta of the external particles
vector<Lorentz5Momentum> pexternal;
// and the internal particles
vector<Lorentz5Momentum> pinter;
// copy the momentum of the incoming particle
pexternal.push_back(pin); pinter.push_back(pin);
pexternal.resize(_mode->numberofParticles());
pinter.resize(_intpart.size());
// masses of the intermediate particles
vector<Energy> massint;
massint.resize(_intpart.size());
massint[0]=pin.mass();
// generate all the decays in the chain
Energy lower,upper,lowerb[2];
for(ix=0;ix<_intpart.size();++ix) {
idau[0] = abs(_intdau1[ix]);
idau[1] = abs(_intdau2[ix]);
// if both decay products off-shell
if(_intdau1[ix]<0&&_intdau2[ix]<0) {
// lower limits on the masses of the two resonances
for(iy=0;iy<2;++iy) {
lowerb[iy]=ZERO;
for(iz=0;iz<_intext[idau[iy]].size();++iz) {
lowerb[iy]+=massext[_intext[idau[iy]][iz]];
}
}
// randomize the order
if(UseRandom::rnd()<0.5) {
// mass of the first resonance
upper = massint[ix]-lowerb[1];
lower = lowerb[0];
massint[idau[0]]=generateMass(idau[0],lower,upper);
// mass of the second resonance
upper = massint[ix]-massint[idau[0]];
lower = lowerb[1];
massint[idau[1]]=generateMass(idau[1],lower,upper);
}
else {
// mass of the second resonance
upper = massint[ix]-lowerb[0];
lower = lowerb[1];
massint[idau[1]]=generateMass(idau[1],lower,upper);
// mass of the first resonance
upper = massint[ix]-massint[idau[1]];
lower = lowerb[0];
massint[idau[0]]=generateMass(idau[0],lower,upper);
}
// generate the momenta of the decay products
twoBodyDecay(pinter[ix],massint[idau[0]],massint[idau[1]],
pinter[idau[0]],pinter[idau[1]]);
}
// only first off-shell
else if(_intdau1[ix]<0) {
// compute the limits of integration
upper = massint[ix]-massext[idau[1]];
lower = ZERO;
for(iy=0;iy<_intext[idau[0]].size();++iy) {
lower+=massext[_intext[idau[0]][iy]];
}
massint[idau[0]]=generateMass(idau[0],lower,upper);
// generate the momenta of the decay products
twoBodyDecay(pinter[ix],massint[idau[0]],massext[idau[1]],
pinter[idau[0]],pexternal[idau[1]]);
}
// only second off-shell
else if(_intdau2[ix]<0) {
// compute the limits of integration
upper = massint[ix]-massext[idau[0]];
lower = ZERO;
for(iy=0;iy<_intext[idau[1]].size();++iy) {
lower+=massext[_intext[idau[1]][iy]];
}
massint[idau[1]]=generateMass(idau[1],lower,upper);
// generate the momenta of the decay products
twoBodyDecay(pinter[ix],massext[idau[0]],massint[idau[1]],
pexternal[idau[0]],pinter[idau[1]]);
}
// both on-shell
else {
// generate the momenta of the decay products
twoBodyDecay(pinter[ix],massext[idau[0]],massext[idau[1]],
pexternal[idau[0]],pexternal[idau[1]]);
}
}
// return the external momenta
return pexternal;
}
// generate the weight for this channel given a phase space configuration
double DecayPhaseSpaceChannel::generateWeight(const vector<Lorentz5Momentum> & output) {
using Constants::pi;
// integers for loops
unsigned int ix,iy,idau[2],iz;
// include the prefactor due to the weight of the channel
double wgt=1.;
// work out the masses of the intermediate particles
vector<Energy2> intmass2(_intpart.size());
vector<Energy> intmass(_intpart.size());
Lorentz5Momentum pinter;
for(ix=0;ix<_intpart.size();++ix) {
pinter=output[_intext[ix][0]];
for(iz=1;iz<_intext[ix].size();++iz) pinter+=output[_intext[ix][iz]];
pinter.rescaleMass();
intmass[ix] = pinter.mass();
intmass2[ix] = sqr(intmass[ix]);
}
Energy2 scale(intmass2[0]);
// calculate the terms for each of the decays
Energy lower,upper,lowerb[2];
for(ix=0;ix<_intpart.size();++ix) {
idau[0] = abs(_intdau1[ix]);
idau[1] = abs(_intdau2[ix]);
// if both decay products off-shell
Energy pcm;
if(_intdau1[ix]<0&&_intdau2[ix]<0) {
// lower limits on the masses of the two resonances
for(iy=0;iy<2;++iy) {
lowerb[iy]=ZERO;
for(iz=0;iz<_intext[idau[iy]].size();++iz)
lowerb[iy]+=output[_intext[idau[iy]][iz]].mass();
}
// undo effect of randomising
// weight for the first order
// contribution of first resonance
upper = intmass[ix]-lowerb[1];
lower = lowerb[0];
InvEnergy2 wgta=massWeight(idau[0],intmass[idau[0]],lower,upper);
// contribution of second resonance
upper = intmass[ix]-intmass[idau[0]];
lower = lowerb[1];
InvEnergy4 wgta2 = wgta*massWeight(idau[1],intmass[idau[1]],lower,upper);
// weight for the second order
upper = intmass[ix]-lowerb[0];
lower = lowerb[1];
InvEnergy2 wgtb=massWeight(idau[1],intmass[idau[1]],lower,upper);
upper = intmass[ix]-intmass[idau[1]];
lower = lowerb[0];
InvEnergy4 wgtb2=wgtb*massWeight(idau[0],intmass[idau[0]],lower,upper);
// weight factor
wgt *=0.5*sqr(scale)*(wgta2+wgtb2);
// factor for the kinematics
pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],intmass[idau[0]],
intmass[idau[1]]);
wgt *= intmass[ix]*8.*pi*pi/pcm;
}
// only first off-shell
else if(_intdau1[ix]<0) {
// compute the limits of integration
upper = intmass[ix]-output[idau[1]].mass();
lower = ZERO;
for(iy=0;iy<_intext[idau[0]].size();++iy)
lower+=output[_intext[idau[0]][iy]].mass();
wgt *=scale*massWeight(idau[0],intmass[idau[0]],lower,upper);
pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],intmass[idau[0]],
output[idau[1]].mass());
wgt *= intmass[ix]*8.*pi*pi/pcm;
}
// only second off-shell
else if(_intdau2[ix]<0) {
// compute the limits of integration
upper = intmass[ix]-output[idau[0]].mass();
lower = ZERO;
for(iy=0;iy<_intext[idau[1]].size();++iy)
lower+=output[_intext[idau[1]][iy]].mass();
wgt *=scale*massWeight(idau[1],intmass[idau[1]],lower,upper);
pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],intmass[idau[1]],
output[idau[0]].mass());
wgt *=intmass[ix]*8.*pi*pi/pcm;
}
// both on-shell
else {
pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],output[idau[1]].mass(),
output[idau[0]].mass());
wgt *=intmass[ix]*8.*pi*pi/pcm;
}
}
// finally the overall factor
wgt /= pi;
// return the answer
return wgt;
}
// output the information to a stream
ostream & Herwig::operator<<(ostream & os, const DecayPhaseSpaceChannel & channel) {
// output of the external particles
os << "Channel for the decay of " << channel._mode->externalParticles(0)->PDGName()
<< " -> ";
for(unsigned int ix=1;ix<channel._mode->numberofParticles();++ix)
os << channel._mode->externalParticles(ix)->PDGName() << " ";
os << endl;
os << "Decay proceeds in following steps ";
for(unsigned int ix=0;ix<channel._intpart.size();++ix) {
os << channel._intpart[ix]->PDGName() << " -> ";
if(channel._intdau1[ix]>0) {
os << channel._mode->externalParticles(channel._intdau1[ix])->PDGName()
<< "(" << channel._intdau1[ix]<< ") ";
}
else {
os << channel._intpart[-channel._intdau1[ix]]->PDGName()
<< "(" << channel._intdau1[ix]<< ") ";
}
if(channel._intdau2[ix]>0) {
os << channel._mode->externalParticles(channel._intdau2[ix])->PDGName()
<< "(" <<channel._intdau2[ix] << ") ";
}
else{
os << channel._intpart[-channel._intdau2[ix]]->PDGName()
<< "(" <<channel._intdau2[ix] << ") ";
}
os << endl;
}
return os;
}
// doinit method
void DecayPhaseSpaceChannel::doinit() {
Interfaced::doinit();
// check if the mode pointer exists
if(!_mode){throw InitException() << "DecayPhaseSpaceChannel::doinit() the "
<< "channel must have a pointer to a decay mode "
<< Exception::abortnow;}
// masses and widths of the intermediate particles
for(unsigned int ix=0;ix<_intpart.size();++ix) {
_intmass.push_back(_intpart[ix]->mass());
_intwidth.push_back(_intpart[ix]->width());
_intmass2.push_back(_intpart[ix]->mass()*_intpart[ix]->mass());
_intmwidth.push_back(_intpart[ix]->mass()*_intpart[ix]->width());
}
// external particles for each intermediate particle
vector<int> temp;
_intext.resize(_intpart.size());
// loop over the intermediate particles
for(int ix=_intpart.size()-1;ix>=0;--ix) {
temp.clear();
// add the first daughter
if(_intdau1[ix]>=0) {
temp.push_back(_intdau1[ix]);
}
else {
int iy = -_intdau1[ix];
vector<int>::iterator istart=_intext[iy].begin();
vector<int>::iterator iend=_intext[iy].end();
for(;istart!=iend;++istart) temp.push_back(*istart);
}
// add the second daughter
if(_intdau2[ix]>=0) {
temp.push_back(_intdau2[ix]);
}
else {
int iy = -_intdau2[ix];
vector<int>::iterator istart=_intext[iy].begin();
vector<int>::iterator iend=_intext[iy].end();
for(;istart!=iend;++istart) temp.push_back(*istart);
}
_intext[ix]=temp;
}
// ensure intermediates either have the width set, or
// can't possibly be on-shell
Energy massmax;
if(_mode->testOnShell()) {
massmax = _mode->externalParticles(0)->mass();
for(unsigned int ix=1;ix<_mode->numberofParticles();++ix)
massmax -= _mode->externalParticles(ix)->mass();
}
else {
massmax = _mode->externalParticles(0)->massMax();
for(unsigned int ix=1;ix<_mode->numberofParticles();++ix)
massmax -= _mode->externalParticles(ix)->massMin();
}
for(unsigned int ix=0;ix<_intpart.size();++ix) {
if(_intwidth.back()==ZERO && ix>0 && _jactype[ix]==0 ) {
Energy massmin(ZERO);
for(unsigned int iy=0;iy<_intext[ix].size();++iy)
massmin += _mode->testOnShell() ?
_mode->externalParticles(_intext[ix][iy])->mass() :
_mode->externalParticles(_intext[ix][iy])->massMin();
// check if can be on-shell
if(_intmass[ix]>=massmin&&_intmass[ix]<=massmax+massmin) {
string modeout;
for(unsigned int iy=0;iy<_mode->numberofParticles();++iy) {
modeout += _mode->externalParticles(iy)->PDGName() + " ";
}
throw InitException() << "Width zero for " << _intpart[ix]->PDGName()
<< " in DecayPhaseSpaceChannel::doinit() "
<< modeout
<< Exception::runerror;
}
}
}
}
void DecayPhaseSpaceChannel::doinitrun() {
Interfaced::doinitrun();
if(!_mode->testOnShell()) return;
_intmass.clear();
_intwidth.clear();
_intmass2.clear();
_intmwidth.clear();
// masses and widths of the intermediate particles
for(unsigned int ix=0;ix<_intpart.size();++ix) {
_intmass.push_back(_intpart[ix]->mass());
_intwidth.push_back(_intpart[ix]->width());
_intmass2.push_back(_intpart[ix]->mass()*_intpart[ix]->mass());
_intmwidth.push_back(_intpart[ix]->mass()*_intpart[ix]->width());
}
// ensure intermediates either have the width set, or
// can't possibly be on-shell
Energy massmax = _mode->externalParticles(0)->massMax();
for(unsigned int ix=1;ix<_mode->numberofParticles();++ix)
massmax -= _mode->externalParticles(ix)->massMin();
for(unsigned int ix=0;ix<_intpart.size();++ix) {
if(_intwidth.back()==0.*MeV && ix>0 && _jactype[ix]==0 ) {
Energy massmin(0.*GeV);
for(unsigned int iy=0;iy<_intext[ix].size();++iy)
massmin += _mode->externalParticles(_intext[ix][iy])->massMin();
// check if can be on-shell
if(_intmass[ix]>=massmin&&_intmass[ix]<=massmax+massmin) {
string modeout;
for(unsigned int iy=0;iy<_mode->numberofParticles();++iy) {
modeout += _mode->externalParticles(iy)->PDGName() + " ";
}
throw Exception() << "Width zero for " << _intpart[ix]->PDGName()
<< " in DecayPhaseSpaceChannel::doinitrun() "
<< modeout
<< Exception::runerror;
}
}
}
}
// generate the final-state particles including the intermediate resonances
void DecayPhaseSpaceChannel::generateIntermediates(bool cc, const Particle & in,
ParticleVector & out) {
// integers for the loops
unsigned int ix,iz;
// create the particles
// incoming particle
ParticleVector external;
external.push_back(const_ptr_cast<tPPtr>(&in));
// outgoing
for(ix=0;ix<out.size();++ix) external.push_back(out[ix]);
out.clear();
// now create the intermediates
ParticleVector resonance; resonance.push_back(external[0]);
PPtr respart;
tcPDPtr parttemp;
Lorentz5Momentum pinter;
for(ix=1;ix<_intpart.size();++ix) {
pinter=external[_intext[ix][0]]->momentum();
for(iz=1;iz<_intext[ix].size();++iz)
pinter+=external[_intext[ix][iz]]->momentum();
pinter.rescaleMass();
respart = (cc&&_intpart[ix]->CC()) ?
_intpart[ix]->CC()->produceParticle(pinter) :
_intpart[ix] ->produceParticle(pinter);
resonance.push_back(respart);
}
// set up the mother daughter relations
for(ix=1;ix<_intpart.size();++ix) {
resonance[ix]->addChild( _intdau1[ix]<0 ?
resonance[-_intdau1[ix]] : external[_intdau1[ix]]);
resonance[ix]->addChild( _intdau2[ix]<0 ?
resonance[-_intdau2[ix]] : external[_intdau2[ix]]);
}
// construct the output with the particles in the first step
out.push_back( _intdau1[0]>0 ? external[_intdau1[0]] : resonance[-_intdau1[0]]);
out.push_back( _intdau2[0]>0 ? external[_intdau2[0]] : resonance[-_intdau2[0]]);
}
double DecayPhaseSpaceChannel::atanhelper_(int ires, Energy limit) {
return atan2( limit*limit-_intmass2[ires], _intmwidth[ires] );
}
// return the weight for a given resonance
InvEnergy2 DecayPhaseSpaceChannel::massWeight(int ires, Energy moff,
Energy lower,Energy upper) {
InvEnergy2 wgt = ZERO;
if(lower>upper) {
throw DecayPhaseSpaceError() << "DecayPhaseSpaceChannel::massWeight not allowed"
<< ires << " " << _intpart[ires]->id() << " "
<< moff/GeV << Exception::eventerror;
}
// use a Breit-Wigner
if ( _jactype[ires] == 0 ) {
double rhomin = atanhelper_(ires,lower);
double rhomax = atanhelper_(ires,upper) - rhomin;
if ( rhomax != 0.0 ) {
Energy2 moff2=moff*moff-_intmass2[ires];
wgt = _intmwidth[ires]/rhomax/(moff2*moff2+_intmwidth[ires]*_intmwidth[ires]);
}
else {
wgt = 1./((sqr(upper)-sqr(lower))*sqr(sqr(moff)-_intmass2[ires])/
(sqr(lower)-_intmass2[ires])/(sqr(upper)-_intmass2[ires]));
}
}
// power law
else if(_jactype[ires]==1) {
double rhomin = pow(sqr(lower/MeV),_intpower[ires]+1.);
double rhomax = pow(sqr(upper/MeV),_intpower[ires]+1.)-rhomin;
wgt = (_intpower[ires]+1.)/rhomax*pow(sqr(moff/MeV),_intpower[ires])
/MeV/MeV;
}
else if(_jactype[ires]==2) {
wgt = 1./Constants::pi/_intmwidth[ires];
}
else {
throw DecayPhaseSpaceError() << "Unknown type of Jacobian in "
<< "DecayPhaseSpaceChannel::massWeight"
<< Exception::eventerror;
}
return wgt;
}
Energy DecayPhaseSpaceChannel::generateMass(int ires,Energy lower,Energy upper) {
static const Energy eps=1e-9*MeV;
if(lower<eps) lower=eps;
Energy mass=ZERO;
if(lower>upper) throw DecayPhaseSpaceError() << "DecayPhaseSpaceChannel::generateMass"
<< " not allowed"
<< Exception::eventerror;
- if(abs(lower-upper)/(lower+upper)>1e-10) {
+ if(abs(lower-upper)/(lower+upper)>2e-10) {
lower +=1e-10*(lower+upper);
upper -=1e-10*(lower+upper);
}
else
return 0.5*(lower+upper);
// use a Breit-Wigner
if(_jactype[ires]==0) {
if(_intmwidth[ires]!=ZERO) {
Energy2 lower2 = sqr(lower);
Energy2 upper2 = sqr(upper);
double rhomin = atan2((lower2 - _intmass2[ires]),_intmwidth[ires]);
double rhomax = atan2((upper2 - _intmass2[ires]),_intmwidth[ires])-rhomin;
double rho = rhomin+rhomax*UseRandom::rnd();
Energy2 mass2 = max(lower2,min(upper2,_intmass2[ires]+_intmwidth[ires]*tan(rho)));
if(mass2<ZERO) mass2 = ZERO;
mass = sqrt(mass2);
}
else {
mass = sqrt(_intmass2[ires]+
(sqr(lower)-_intmass2[ires])*(sqr(upper)-_intmass2[ires])/
(sqr(lower)-_intmass2[ires]-UseRandom::rnd()*(sqr(lower)-sqr(upper))));
}
}
// use a power-law
else if(_jactype[ires]==1) {
double rhomin = pow(sqr(lower/MeV),_intpower[ires]+1.);
double rhomax = pow(sqr(upper/MeV),_intpower[ires]+1.)-rhomin;
double rho = rhomin+rhomax*UseRandom::rnd();
mass = pow(rho,0.5/(_intpower[ires]+1.))*MeV;
}
else if(_jactype[ires]==2) {
mass = _intmass[ires];
}
else {
throw DecayPhaseSpaceError() << "Unknown type of Jacobian in "
<< "DecayPhaseSpaceChannel::generateMass"
<< Exception::eventerror;
}
if(mass<lower) mass=lower+1e-10*(lower+upper);
else if(mass>upper) mass=upper-1e-10*(lower+upper);
return mass;
}
void DecayPhaseSpaceChannel::twoBodyDecay(const Lorentz5Momentum & p,
const Energy m1, const Energy m2,
Lorentz5Momentum & p1,
Lorentz5Momentum & p2 ) {
static const double eps=1e-6;
double ctheta,phi;
Kinematics::generateAngles(ctheta,phi);
Axis unitDir1=Kinematics::unitDirection(ctheta,phi);
Momentum3 pstarVector;
Energy min=p.m();
if ( min >= m1 + m2 && m1 >= ZERO && m2 >= ZERO ) {
pstarVector = unitDir1 * Kinematics::pstarTwoBodyDecay(min,m1,m2);
}
else if( m1 >= ZERO && m2 >= ZERO && (m1+m2-min)/(min+m1+m2)<eps) {
pstarVector = Momentum3();
}
else {
throw DecayPhaseSpaceError() << "Two body decay cannot proceed "
<< "p = " << p / GeV
<< " p.m() = " << min / GeV
<< " -> " << m1/GeV
<< ' ' << m2/GeV << Exception::eventerror;
}
p1 = Lorentz5Momentum(m1, pstarVector);
p2 = Lorentz5Momentum(m2,-pstarVector);
// boost from CM to LAB
Boost bv = p.vect() * (1./p.t());
p1.boost( bv );
p2.boost( bv );
}
diff --git a/MatrixElement/DIS/DISBase.cc b/MatrixElement/DIS/DISBase.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/DIS/DISBase.cc
@@ -0,0 +1,1371 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the DISBase class.
+//
+
+#include "DISBase.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Reference.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
+#include "Herwig++/Utilities/Maths.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "ThePEG/PDT/StandardMatchers.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Repository/CurrentGenerator.h"
+#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
+#include "Herwig++/PDT/StandardMatchers.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include <numeric>
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+namespace {
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+void debuggingMatrixElement(bool BGF,const Lorentz5Momentum & pin,
+ const Lorentz5Momentum & p1,
+ const Lorentz5Momentum & p2,
+ tcPDPtr gluon,
+ const Lorentz5Momentum & pl1,
+ const Lorentz5Momentum & pl2,
+ const Lorentz5Momentum & pq1,
+ const Lorentz5Momentum & pq2,
+ tcPDPtr lepton1,tcPDPtr lepton2,
+ tcPDPtr quark1 ,tcPDPtr quark2,
+ Energy2 Q2,double phi, double x2, double x3,
+ double xperp, double zp, double xp,
+ const vector<double> & azicoeff,
+ bool normalize) {
+ tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>
+ (CurrentGenerator::current().standardModel());
+ assert(hwsm);
+ vector<AbstractFFVVertexPtr> weakVertex;
+ vector<PDPtr> bosons;
+ AbstractFFVVertexPtr strongVertex = hwsm->vertexFFG();
+ if(lepton1->id()==lepton2->id()) {
+ weakVertex.push_back(hwsm->vertexFFZ());
+ bosons.push_back(hwsm->getParticleData(ParticleID::Z0));
+ weakVertex.push_back(hwsm->vertexFFP());
+ bosons.push_back(hwsm->getParticleData(ParticleID::gamma));
+ }
+ else {
+ weakVertex.push_back(hwsm->vertexFFW());
+ bosons.push_back(hwsm->getParticleData(ParticleID::Wplus));
+ }
+ if(!BGF) {
+ SpinorWaveFunction l1,q1,qp1;
+ SpinorBarWaveFunction l2,q2,qp2;
+ VectorWaveFunction gl(p2,gluon,outgoing);
+ if(lepton1->id()>0) {
+ l1 = SpinorWaveFunction (pl1,lepton1,incoming);
+ l2 = SpinorBarWaveFunction(pl2,lepton2,outgoing);
+ }
+ else {
+ l1 = SpinorWaveFunction (pl2,lepton2,outgoing);
+ l2 = SpinorBarWaveFunction(pl1,lepton1,incoming);
+ }
+ if(quark1->id()>0) {
+ q1 = SpinorWaveFunction (pq1,quark1,incoming);
+ q2 = SpinorBarWaveFunction(pq2,quark2,outgoing);
+ qp1 = SpinorWaveFunction (pin,quark1,incoming);
+ qp2 = SpinorBarWaveFunction(p1 ,quark2,outgoing);
+ }
+ else {
+ q1 = SpinorWaveFunction (pq2,quark2,outgoing);
+ q2 = SpinorBarWaveFunction(pq1,quark1,incoming);
+ qp1 = SpinorWaveFunction (p1 ,quark2,outgoing);
+ qp2 = SpinorBarWaveFunction(pin,quark1,incoming);
+ }
+ double lome(0.),realme(0.);
+ for(unsigned int lhel1=0;lhel1<2;++lhel1) {
+ l1.reset(lhel1);
+ for(unsigned int lhel2=0;lhel2<2;++lhel2) {
+ l2.reset(lhel2);
+ for(unsigned int qhel1=0;qhel1<2;++qhel1) {
+ q1.reset(qhel1);
+ qp1.reset(qhel1);
+ for(unsigned int qhel2=0;qhel2<2;++qhel2) {
+ q2.reset(qhel2);
+ qp2.reset(qhel2);
+ // leading order matrix element
+ Complex diagLO(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ diagLO += weakVertex[ix]->evaluate(Q2,q1,q2,inter);
+ }
+ lome += norm(diagLO);
+ // real emission matrix element
+ for(unsigned int ghel=0;ghel<2;++ghel) {
+ gl.reset(2*ghel);
+ Complex diagReal(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ SpinorWaveFunction off1 =
+ strongVertex->evaluate(Q2,5,qp1.particle(),qp1,gl);
+ Complex diag1 = weakVertex[ix]->evaluate(Q2,off1,qp2,inter);
+ SpinorBarWaveFunction off2 =
+ strongVertex->evaluate(Q2,5,qp2.particle(),qp2,gl);
+ Complex diag2 = weakVertex[ix]->evaluate(Q2,qp1,off2,inter);
+ diagReal += diag1+diag2;
+ }
+ realme += norm(diagReal);
+ }
+ }
+ }
+ }
+ }
+ double test1 = realme/lome/hwsm->alphaS(Q2)*Q2*UnitRemoval::InvE2;
+ double cphi(cos(phi));
+ double test2;
+ if(normalize) {
+ test2 = 8.*Constants::pi/(1.-xp)/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi))*
+ (1.+sqr(xp)*(sqr(x2)+1.5*sqr(xperp)));
+ }
+ else {
+ test2 = 8.*Constants::pi/(1.-xp)/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi));
+ }
+ cerr << "testing RATIO A " << test1/test2 << "\n";
+ }
+ else {
+ SpinorWaveFunction l1,q1,qp1;
+ SpinorBarWaveFunction l2,q2,qp2;
+ VectorWaveFunction gl(pin,gluon,incoming);
+ if(lepton1->id()>0) {
+ l1 = SpinorWaveFunction (pl1,lepton1,incoming);
+ l2 = SpinorBarWaveFunction(pl2,lepton2,outgoing);
+ }
+ else {
+ l1 = SpinorWaveFunction (pl2,lepton2,outgoing);
+ l2 = SpinorBarWaveFunction(pl1,lepton1,incoming);
+ }
+ if(quark1->id()>0) {
+ q1 = SpinorWaveFunction (pq1,quark1 ,incoming);
+ q2 = SpinorBarWaveFunction(pq2,quark2 ,outgoing);
+ qp2 = SpinorBarWaveFunction(p1 ,quark2 ,outgoing);
+ qp1 = SpinorWaveFunction (p2 ,quark1->CC(),outgoing);
+ }
+ else {
+ q1 = SpinorWaveFunction (pq2,quark2 ,outgoing);
+ q2 = SpinorBarWaveFunction(pq1,quark1 ,incoming);
+ qp2 = SpinorBarWaveFunction(p2 ,quark1->CC(),outgoing);
+ qp1 = SpinorWaveFunction (p1 ,quark2 ,outgoing);
+ }
+ double lome(0.),realme(0.);
+ for(unsigned int lhel1=0;lhel1<2;++lhel1) {
+ l1.reset(lhel1);
+ for(unsigned int lhel2=0;lhel2<2;++lhel2) {
+ l2.reset(lhel2);
+ for(unsigned int qhel1=0;qhel1<2;++qhel1) {
+ q1.reset(qhel1);
+ qp1.reset(qhel1);
+ for(unsigned int qhel2=0;qhel2<2;++qhel2) {
+ q2.reset(qhel2);
+ qp2.reset(qhel2);
+ // leading order matrix element
+ Complex diagLO(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ diagLO += weakVertex[ix]->evaluate(Q2,q1,q2,inter);
+ }
+ lome += norm(diagLO);
+ // real emission matrix element
+ for(unsigned int ghel=0;ghel<2;++ghel) {
+ gl.reset(2*ghel);
+ Complex diagReal(0.);
+ for(unsigned int ix=0;ix<weakVertex.size();++ix) {
+ VectorWaveFunction inter =
+ weakVertex[ix]->evaluate(Q2,3,bosons[ix],l1,l2);
+ SpinorWaveFunction off1 =
+ strongVertex->evaluate(Q2,5,qp1.particle(),qp1,gl);
+ Complex diag1 = weakVertex[ix]->evaluate(Q2,off1,qp2,inter);
+ SpinorBarWaveFunction off2 =
+ strongVertex->evaluate(Q2,5,qp2.particle(),qp2,gl);
+ Complex diag2 = weakVertex[ix]->evaluate(Q2,qp1,off2,inter);
+ diagReal += diag1+diag2;
+ }
+ realme += norm(diagReal);
+ }
+ }
+ }
+ }
+ }
+ double test1 = realme/lome/hwsm->alphaS(Q2)*Q2*UnitRemoval::InvE2;
+ double cphi(cos(phi));
+ double test2;
+ if(normalize) {
+ test2 = 8.*Constants::pi/zp/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi))*
+ sqr(xp)*(sqr(x3)+sqr(x2)+3.*sqr(xperp));
+ }
+ else {
+ test2 = 8.*Constants::pi/zp/(1.-zp)*
+ (azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi));
+ }
+ cerr << "testing RATIO B " << test1/test2 << "\n";
+ }
+}
+
+}
+
+DISBase::DISBase() : initial_(6.), final_(3.),
+ procProb_(0.35),
+ comptonInt_(0.), bgfInt_(0.),
+ comptonWeight_(50.), BGFWeight_(150.),
+ pTmin_(0.1*GeV),
+ scaleOpt_(1), muF_(100.*GeV), scaleFact_(1.),
+ contrib_(0), power_(0.1)
+{}
+
+DISBase::~DISBase() {}
+
+void DISBase::persistentOutput(PersistentOStream & os) const {
+ os << comptonInt_ << bgfInt_ << procProb_ << initial_ << final_ << alpha_
+ << ounit(pTmin_,GeV) << comptonWeight_ << BGFWeight_ << gluon_
+ << ounit(muF_,GeV) << scaleFact_ << scaleOpt_ << contrib_<< power_;
+}
+
+void DISBase::persistentInput(PersistentIStream & is, int) {
+ is >> comptonInt_ >> bgfInt_ >> procProb_ >> initial_ >> final_ >> alpha_
+ >> iunit(pTmin_,GeV) >> comptonWeight_ >> BGFWeight_ >> gluon_
+ >> iunit(muF_,GeV) >> scaleFact_ >> scaleOpt_ >> contrib_ >> power_;
+}
+
+AbstractClassDescription<DISBase> DISBase::initDISBase;
+// Definition of the static class description member.
+
+void DISBase::Init() {
+
+ static ClassDocumentation<DISBase> documentation
+ ("The DISBase class provides the base class for the "
+ "implementation of DIS type processes including the "
+ "hard corrections in either the old-fashioned matrix "
+ "element correction of POWHEG approaches");
+
+ static Parameter<DISBase,double> interfaceProcessProbability
+ ("ProcessProbability",
+ "The probabilty of the QCD compton process for the process selection",
+ &DISBase::procProb_, 0.3, 0.0, 1.,
+ false, false, Interface::limited);
+
+ static Reference<DISBase,ShowerAlpha> interfaceCoupling
+ ("Coupling",
+ "Pointer to the object to calculate the coupling for the correction",
+ &DISBase::alpha_, false, false, true, false, false);
+
+ static Parameter<DISBase,Energy> interfacepTMin
+ ("pTMin",
+ "The minimum pT",
+ &DISBase::pTmin_, GeV, 1.*GeV, 0.0*GeV, 10.0*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceComptonWeight
+ ("ComptonWeight",
+ "Weight for the overestimate ofthe compton channel",
+ &DISBase::comptonWeight_, 50.0, 0.0, 100.0,
+ false, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceBGFWeight
+ ("BGFWeight",
+ "Weight for the overestimate of the BGF channel",
+ &DISBase::BGFWeight_, 100.0, 0.0, 1000.0,
+ false, false, Interface::limited);
+
+ static Switch<DISBase,unsigned int> interfaceContribution
+ ("Contribution",
+ "Which contributions to the cross section to include",
+ &DISBase::contrib_, 0, false, false);
+ static SwitchOption interfaceContributionLeadingOrder
+ (interfaceContribution,
+ "LeadingOrder",
+ "Just generate the leading order cross section",
+ 0);
+ static SwitchOption interfaceContributionPositiveNLO
+ (interfaceContribution,
+ "PositiveNLO",
+ "Generate the positive contribution to the full NLO cross section",
+ 1);
+ static SwitchOption interfaceContributionNegativeNLO
+ (interfaceContribution,
+ "NegativeNLO",
+ "Generate the negative contribution to the full NLO cross section",
+ 2);
+
+ static Switch<DISBase,unsigned int> interfaceScaleOption
+ ("ScaleOption",
+ "Option for the choice of factorization (and renormalization) scale",
+ &DISBase::scaleOpt_, 1, false, false);
+ static SwitchOption interfaceDynamic
+ (interfaceScaleOption,
+ "Dynamic",
+ "Dynamic factorization scale equal to the current sqrt(sHat())",
+ 1);
+ static SwitchOption interfaceFixed
+ (interfaceScaleOption,
+ "Fixed",
+ "Use a fixed factorization scale set with FactorizationScaleValue",
+ 2);
+
+ static Parameter<DISBase,Energy> interfaceFactorizationScale
+ ("FactorizationScale",
+ "Value to use in the event of a fixed factorization scale",
+ &DISBase::muF_, GeV, 100.0*GeV, 1.0*GeV, 500.0*GeV,
+ true, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceScaleFactor
+ ("ScaleFactor",
+ "The factor used before Q2 if using a running scale",
+ &DISBase::scaleFact_, 1.0, 0.0, 10.0,
+ false, false, Interface::limited);
+
+ static Parameter<DISBase,double> interfaceSamplingPower
+ ("SamplingPower",
+ "Power for the sampling of xp",
+ &DISBase::power_, 0.6, 0.0, 1.,
+ false, false, Interface::limited);
+}
+
+void DISBase::doinit() {
+ HwMEBase::doinit();
+ // integrals of me over phase space
+ double r5=sqrt(5.),darg((r5-1.)/(r5+1.)),ath(0.5*log((1.+1./r5)/(1.-1./r5)));
+ comptonInt_ = 2.*(-21./20.-6./(5.*r5)*ath+sqr(Constants::pi)/3.
+ -2.*Math::ReLi2(1.-darg)-2.*Math::ReLi2(1.-1./darg));
+ bgfInt_ = 121./9.-56./r5*ath;
+ // extract the gluon ParticleData objects
+ gluon_ = getParticleData(ParticleID::g);
+}
+
+void DISBase::initializeMECorrection(ShowerTreePtr tree, double & initial,
+ double & final) {
+ initial = initial_;
+ final = final_;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ partons_[0] = cit->first->progenitor()->dataPtr();
+ pq_[0] = cit->first->progenitor()->momentum();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ leptons_[0] = cit->first->progenitor()->dataPtr();
+ pl_[0] = cit->first->progenitor()->momentum();
+ }
+ }
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ partons_[1] = cit->first->progenitor()->dataPtr();
+ pq_[1] = cit->first->progenitor()->momentum();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ leptons_[1] = cit->first->progenitor()->dataPtr();
+ pl_[1] = cit->first->progenitor()->momentum();
+ }
+ }
+ // extract the born variables
+ q_ =pl_[0]-pl_[1];
+ q2_ = -q_.m2();
+ double yB = (q_*pq_[0])/(pl_[0]*pq_[0]);
+ l_ = 2./yB-1.;
+ // calculate the A coefficient for the correlations
+ acoeff_ = A(leptons_[0],leptons_[1],
+ partons_[0],partons_[1],q2_);
+}
+
+void DISBase::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
+ static const double eps=1e-6;
+ // find the incoming and outgoing quarks and leptons
+ ShowerParticlePtr quark[2],lepton[2];
+ PPtr hadron;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ hadron = cit->first->original()->parents()[0];
+ quark [0] = cit->first->progenitor();
+ beam_ = cit->first->beam();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[0] = cit->first->progenitor();
+ }
+ }
+ pdf_ = beam_->pdf();
+ assert(beam_&&pdf_&&quark[0]&&lepton[0]);
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ quark [1] = cit->first->progenitor();
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[1] = cit->first->progenitor();
+ }
+ }
+ // momentum fraction
+ assert(quark[1]&&lepton[1]);
+ xB_ = quark[0]->x();
+ // calculate the matrix element
+ vector<double> azicoeff;
+ // select the type of process
+ bool BGF = UseRandom::rnd()>procProb_;
+ double xp,zp,wgt,x1,x2,x3,xperp;
+ // generate a QCD compton process
+ if(!BGF) {
+ wgt = generateComptonPoint(xp,zp);
+ if(xp<eps) return;
+ // common pieces
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
+ wgt *= 2./3./Constants::pi*alpha_->value(scale)/procProb_;
+ // PDF piece
+ wgt *= pdf_->xfx(beam_,quark[0]->dataPtr(),scale,xB_/xp)/
+ pdf_->xfx(beam_,quark[0]->dataPtr(),q2_ ,xB_);
+ // other bits
+ xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ x1 = -1./xp;
+ x2 = 1.-(1.-zp)/xp;
+ x3 = 2.+x1-x2;
+ // matrix element pieces
+ azicoeff = ComptonME(xp,x2,xperp,true);
+ }
+ // generate a BGF process
+ else {
+ wgt = generateBGFPoint(xp,zp);
+ if(xp<eps) return;
+ // common pieces
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1);
+ wgt *= 0.25/Constants::pi*alpha_->value(scale)/(1.-procProb_);
+ // PDF piece
+ wgt *= pdf_->xfx(beam_,gluon_ ,scale,xB_/xp)/
+ pdf_->xfx(beam_,quark[0]->dataPtr(),q2_ ,xB_);
+ // other bits
+ xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ x1 = -1./xp;
+ x2 = 1.-(1.-zp)/xp;
+ x3 = 2.+x1-x2;
+ // matrix element pieces
+ azicoeff = BGFME(xp,x2,x3,xperp,true);
+ }
+ // compute the azimuthal average of the weight
+ wgt *= (azicoeff[0]+0.5*azicoeff[2]);
+ // decide whether or not to accept the weight
+ if(UseRandom::rnd()>wgt) return;
+ // if generate generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi));
+ phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
+ << "::applyHardMatrixElementCorrection() to"
+ << " generate phi" << Exception::eventerror;
+ // construct lorentz transform from lab to breit frame
+ Lorentz5Momentum phadron = hadron->momentum();
+ phadron.setMass(0.*GeV);
+ phadron.rescaleEnergy();
+ Lorentz5Momentum pcmf = phadron+0.5/xB_*q_;
+ pcmf.rescaleMass();
+ LorentzRotation rot(-pcmf.boostVector());
+ Lorentz5Momentum pbeam = rot*phadron;
+ Axis axis(pbeam.vect().unit());
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ rot.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ Lorentz5Momentum pl = rot*pl_[0];
+ rot.rotateZ(-atan2(pl.y(),pl.x()));
+ pl_[0] *= rot;
+ pl_[1] *= rot;
+ pq_[0] *= rot;
+ pq_[1] *= rot;
+ // compute the new incoming and outgoing momenta
+ Energy Q(sqrt(q2_));
+ Lorentz5Momentum p1 = Lorentz5Momentum( 0.5*Q*xperp*cos(phi), 0.5*Q*xperp*sin(phi),
+ -0.5*Q*x2,0.*GeV,0.*GeV);
+ p1.rescaleEnergy();
+ Lorentz5Momentum p2 = Lorentz5Momentum(-0.5*Q*xperp*cos(phi),-0.5*Q*xperp*sin(phi),
+ -0.5*Q*x3,0.*GeV,0.*GeV);
+ p2.rescaleEnergy();
+ Lorentz5Momentum pin(0.*GeV,0.*GeV,-0.5*x1*Q,-0.5*x1*Q,0.*GeV);
+// debuggingMatrixElement(BGF,pin,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
+// lepton[0]->dataPtr(),lepton[1]->dataPtr(),
+// quark [0]->dataPtr(),quark [1]->dataPtr(),
+// q2_,phi,x2,x3,xperp,zp,xp,azicoeff,true);
+ // we need the Lorentz transform back to the lab
+ rot.invert();
+ // transform the momenta to lab frame
+ pin *= rot;
+ p1 *= rot;
+ p2 *= rot;
+ // test to ensure outgoing particles can be put on-shell
+ if(!BGF) {
+ if(p1.e()<quark[1]->dataPtr()->constituentMass()) return;
+ if(p2.e()<gluon_ ->constituentMass()) return;
+ }
+ else {
+ if(p1.e()<quark[1]->dataPtr() ->constituentMass()) return;
+ if(p2.e()<quark[0]->dataPtr()->CC()->constituentMass()) return;
+ }
+ // create the new particles and add to ShowerTree
+ bool isquark = quark[0]->colourLine();
+ if(!BGF) {
+ PPtr newin = new_ptr(Particle(*quark[0]));
+ newin->set5Momentum(pin);
+ PPtr newg = gluon_ ->produceParticle(p2 );
+ PPtr newout = quark[1]->dataPtr()->produceParticle(p1 );
+ ColinePtr col=isquark ?
+ quark[0]->colourLine() : quark[0]->antiColourLine();
+ ColinePtr newline=new_ptr(ColourLine());
+ // final-state emission
+ if(xp>zp) {
+ col->removeColoured(newout,!isquark);
+ col->addColoured(newin,!isquark);
+ col->addColoured(newg,!isquark);
+ newline->addColoured(newg,isquark);
+ newline->addColoured(newout,!isquark);
+ }
+ // initial-state emission
+ else {
+ col->removeColoured(newin ,!isquark);
+ col->addColoured(newout,!isquark);
+ col->addColoured(newg,isquark);
+ newline->addColoured(newg,!isquark);
+ newline->addColoured(newin,!isquark);
+ }
+ PPtr orig;
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[0]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newin);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
+ cit->first->progenitor(sp);
+ tree->incomingLines()[cit->first]=sp;
+ sp->x(xB_/xp);
+ cit->first->perturbative(xp>zp);
+ if(xp<=zp) orig=cit->first->original();
+ }
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[1]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newout);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
+ cit->first->progenitor(sp);
+ tree->outgoingLines()[cit->first]=sp;
+ cit->first->perturbative(xp<=zp);
+ if(xp>zp) orig=cit->first->original();
+ }
+ assert(orig);
+ // add the gluon
+ ShowerParticlePtr sg=new_ptr(ShowerParticle(*newg,1,true));
+ ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
+ gluon->perturbative(false);
+ tree->outgoingLines().insert(make_pair(gluon,sg));
+ tree->hardMatrixElementCorrection(true);
+ }
+ else {
+ PPtr newin = gluon_ ->produceParticle(pin);
+ PPtr newqbar = quark[0]->dataPtr()->CC()->produceParticle(p2 );
+ PPtr newout = quark[1]->dataPtr() ->produceParticle(p1 );
+ ColinePtr col=isquark ? quark[0]->colourLine() : quark[0]->antiColourLine();
+ ColinePtr newline=new_ptr(ColourLine());
+ col ->addColoured(newin ,!isquark);
+ newline->addColoured(newin , isquark);
+ col ->addColoured(newout ,!isquark);
+ newline->addColoured(newqbar, isquark);
+ PPtr orig;
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[0]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newin);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
+ cit->first->progenitor(sp);
+ tree->incomingLines()[cit->first]=sp;
+ sp->x(xB_/xp);
+ cit->first->perturbative(false);
+ orig=cit->first->original();
+ }
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(cit->first->progenitor()!=quark[1]) continue;
+ // remove old particles from colour line
+ col->removeColoured(cit->first->copy(),!isquark);
+ col->removeColoured(cit->first->progenitor(),!isquark);
+ // insert new particles
+ cit->first->copy(newout);
+ ShowerParticlePtr sp(new_ptr(ShowerParticle(*newout,1,true)));
+ cit->first->progenitor(sp);
+ tree->outgoingLines()[cit->first]=sp;
+ cit->first->perturbative(true);
+ }
+ assert(orig);
+ // add the (anti)quark
+ ShowerParticlePtr sqbar=new_ptr(ShowerParticle(*newqbar,1,true));
+ ShowerProgenitorPtr qbar=new_ptr(ShowerProgenitor(orig,newqbar,sqbar));
+ qbar->perturbative(false);
+ tree->outgoingLines().insert(make_pair(qbar,sqbar));
+ tree->hardMatrixElementCorrection(true);
+ }
+}
+
+bool DISBase::softMatrixElementVeto(ShowerProgenitorPtr initial,
+ ShowerParticlePtr parent, Branching br) {
+ bool veto = !UseRandom::rndbool(parent->isFinalState() ? 1./final_ : 1./initial_);
+ // check if me correction should be applied
+ long id[2]={initial->id(),parent->id()};
+ if(id[0]!=id[1]||id[1]==ParticleID::g) return veto;
+ // get the pT
+ Energy pT=br.kinematics->pT();
+ // check if hardest so far
+ if(pT<initial->highestpT()) return veto;
+ double kappa(sqr(br.kinematics->scale())/q2_),z(br.kinematics->z());
+ double zk((1.-z)*kappa);
+ // final-state
+ double wgt(0.);
+ if(parent->isFinalState()) {
+ double zp=z,xp=1./(1.+z*zk);
+ double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ double x2 = 1.-(1.-zp)/xp;
+ vector<double> azicoeff = ComptonME(xp,x2,xperp,false);
+ wgt = (azicoeff[0]+0.5*azicoeff[2])*xp/(1.+sqr(z))/final_;
+ if(wgt<.0||wgt>1.) {
+ ostringstream wstring;
+ wstring << "Soft ME correction weight too large or "
+ << "negative for FSR in DISBase::"
+ << "softMatrixElementVeto() soft weight "
+ << " xp = " << xp << " zp = " << zp
+ << " weight = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ else {
+ double xp = 2.*z/(1.+zk+sqrt(sqr(1.+zk)-4.*z*zk));
+ double zp = 0.5* (1.-zk+sqrt(sqr(1.+zk)-4.*z*zk));
+ double xperp = sqrt(4.*(1.-xp)*(1.-zp)*zp/xp);
+ double x1 = -1./xp, x2 = 1.-(1.-zp)/xp, x3 = 2.+x1-x2;
+ // compton
+ if(br.ids[0]!=ParticleID::g) {
+ vector<double> azicoeff = ComptonME(xp,x2,xperp,false);
+ wgt = (azicoeff[0]+0.5*azicoeff[2])*xp*(1.-z)/(1.-xp)/(1.+sqr(z))/
+ (1.-zp+xp-2.*xp*(1.-zp));
+ }
+ // BGF
+ else {
+ vector<double> azicoeff = BGFME(xp,x2,x3,xperp,true);
+ wgt = (azicoeff[0]+0.5*azicoeff[2])*xp/(1.-zp+xp-2.*xp*(1.-zp))/(sqr(z)+sqr(1.-z));
+ }
+ wgt /=initial_;
+ if(wgt<.0||wgt>1.) {
+ ostringstream wstring;
+ wstring << "Soft ME correction weight too large or "
+ << "negative for ISR in DISBase::"
+ << "softMatrixElementVeto() soft weight "
+ << " xp = " << xp << " zp = " << zp
+ << " weight = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ // if not vetoed
+ if(UseRandom::rndbool(wgt)) {
+ initial->highestpT(pT);
+ return false;
+ }
+ // otherwise
+ parent->setEvolutionScale(br.kinematics->scale());
+ return true;
+}
+
+double DISBase::generateComptonPoint(double &xp, double & zp) {
+ static const double maxwgt = 1.;
+ double wgt;
+ do {
+ xp = UseRandom::rnd();
+ double zpmin = xp, zpmax = 1./(1.+xp*(1.-xp));
+ zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+ wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+ if(UseRandom::rndbool()) swap(xp,zp);
+ double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp,x2=1.-(1.-zp)/xp;
+ wgt *= 2.*(1.+sqr(xp)*(sqr(x2)+1.5*xperp2))/(1.-xp)/(1.-zp);
+ if(wgt>maxwgt) {
+ ostringstream wstring;
+ wstring << "DISBase::generateComptonPoint "
+ << "Weight greater than maximum "
+ << "wgt = " << wgt << " maxwgt = 1\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(wgt<UseRandom::rnd()*maxwgt);
+ return comptonInt_;
+}
+
+double DISBase::generateBGFPoint(double &xp, double & zp) {
+ static const double maxwgt = 25.;
+ double wgt;
+ do {
+ xp = UseRandom::rnd();
+ double zpmax = 1./(1.+xp*(1.-xp)), zpmin = 1.-zpmax;
+ zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+ wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
+ wgt *= sqr(xp)/(1.-zp)*(sqr(x3)+sqr(x2)+3.*xperp2);
+ if(wgt>maxwgt) {
+ ostringstream wstring;
+ wstring << "DISBase::generateBGFPoint "
+ << "Weight greater than maximum "
+ << "wgt = " << wgt << " maxwgt = 1\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(wgt<UseRandom::rnd()*maxwgt);
+ return bgfInt_;
+// static const double maxwgt = 2.,npow=0.34,ac=1.0;
+// double wgt;
+// do {
+// double rho = UseRandom::rnd();
+// xp = 1.-pow(rho,1./(1.-npow));
+// wgt = (sqr(xp)+ac+sqr(1.-xp));
+// if(wgt>1.+ac) cerr << "testing violates BGF maxA " << wgt << "\n";
+// }
+// while(wgt<UseRandom::rnd()*(1.+ac));
+// double xpwgt = -((6.-5.*npow+sqr(npow))*ac-3.*npow+sqr(npow)+4)
+// /(sqr(npow)*(npow-6.)+11.*npow-6.);
+// xpwgt *= pow(1.-xp,npow)/wgt;
+// double xp2(sqr(xp)),lxp(log(xp)),xp4(sqr(xp2)),lxp1(log(1.-xp));
+// double zpwgt = (2.*xp4*(lxp+lxp1-3.)+4.*xp*xp2*(3.-lxp-lxp1)
+// +xp2*(-13.+lxp+lxp1)+xp*(+7.+lxp+lxp1)-lxp-lxp1-1.)/(1.+xp-xp2);
+// do {
+// double zpmax = 1./(1.+xp*(1.-xp)), zpmin = 1.-zpmax;
+// zp = 1.-pow((1.-zpmin)/(1.-zpmax),UseRandom::rnd())*(1.-zpmax);
+// wgt = log((1.-zpmin)/(1.-zpmax))*(1.-zp);
+// double x1 = -1./xp;
+// double x2 = 1.-(1.-zp)/xp;
+// double x3 = 2.+x1-x2;
+// double xperp2 = 4.*(1.-xp)*(1.-zp)*zp/xp;
+// wgt *= sqr(xp)/(1.-zp)*(sqr(x3)+sqr(x2)+3.*xperp2);
+// if(wgt>maxwgt*zpwgt) cerr << "testing violates BGF maxB " << wgt/xpwgt << "\n";
+// }
+// while(wgt<UseRandom::rnd()*maxwgt);
+// return zpwgt*xpwgt;
+}
+
+vector<double> DISBase::ComptonME(double xp, double x2, double xperp,
+ bool norm) {
+ vector<double> output(3,0.);
+ double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
+ double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
+ double root = sqrt(sqr(l_)-1.);
+ output[0] = sqr(cos2)+acoeff_*cos2*l_+sqr(l_);
+ output[1] = -acoeff_*cos2*root*sin2-2.*l_*root*sin2;
+ output[2] = sqr(root)*sqr(sin2);
+ double lo(1+acoeff_*l_+sqr(l_));
+ double denom = norm ? 1.+sqr(xp)*(sqr(x2)+1.5*sqr(xperp)) : 1.;
+ double fact = sqr(xp)*(sqr(x2)+sqr(xperp))/lo;
+ for(unsigned int ix=0;ix<output.size();++ix)
+ output[ix] = ((ix==0 ? 1. : 0.) +fact*output[ix])/denom;
+ return output;
+}
+
+vector<double> DISBase::BGFME(double xp, double x2, double x3,
+ double xperp, bool norm) {
+ vector<double> output(3,0.);
+ double cos2 = x2 /sqrt(sqr(x2)+sqr(xperp));
+ double sin2 = xperp/sqrt(sqr(x2)+sqr(xperp));
+ double fact2 = sqr(xp)*(sqr(x2)+sqr(xperp));
+ double cos3 = x3 /sqrt(sqr(x3)+sqr(xperp));
+ double sin3 = xperp/sqrt(sqr(x3)+sqr(xperp));
+ double fact3 = sqr(xp)*(sqr(x3)+sqr(xperp));
+ double root = sqrt(sqr(l_)-1.);
+ output[0] = fact2*(sqr(cos2)+acoeff_*cos2*l_+sqr(l_)) +
+ fact3*(sqr(cos3)-acoeff_*cos3*l_+sqr(l_));
+ output[1] = - fact2*(acoeff_*cos2*root*sin2+2.*l_*root*sin2)
+ - fact3*(acoeff_*cos3*root*sin3-2.*l_*root*sin3);
+ output[2] = fact2*(sqr(root)*sqr(sin2)) +
+ fact3*(sqr(root)*sqr(sin3));
+ double lo(1+acoeff_*l_+sqr(l_));
+ double denom = norm ? sqr(xp)*(sqr(x3)+sqr(x2)+3.*sqr(xperp))*lo : lo;
+ for(unsigned int ix=0;ix<output.size();++ix) output[ix] /= denom;
+ return output;
+}
+
+HardTreePtr DISBase::generateHardest(ShowerTreePtr tree) {
+ ShowerParticlePtr quark[2],lepton[2];
+ PPtr hadron;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data())) {
+ hadron = cit->first->original()->parents()[0];
+ quark [0] = cit->first->progenitor();
+ beam_ = cit->first->beam();
+ }
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[0] = cit->first->progenitor();
+ leptons_[0] = lepton[0]->dataPtr();
+ }
+ }
+ pdf_=beam_->pdf();
+ assert(beam_&&pdf_&&quark[0]&&lepton[0]);
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ quark [1] = cit->first->progenitor();
+ else if(LeptonMatcher::Check(cit->first->progenitor()->data())) {
+ lepton[1] = cit->first->progenitor();
+ leptons_[1] = lepton[1]->dataPtr();
+ }
+ }
+ assert(quark[1]&&lepton[1]);
+ // Particle data objects
+ for(unsigned int ix=0;ix<2;++ix) partons_[ix] = quark[ix]->dataPtr();
+ // extract the born variables
+ q_ =lepton[0]->momentum()-lepton[1]->momentum();
+ q2_ = -q_.m2();
+ xB_ = quark[0]->x();
+ double yB =
+ ( q_*quark[0]->momentum())/
+ (lepton[0]->momentum()*quark[0]->momentum());
+ l_ = 2./yB-1.;
+ // construct lorentz transform from lab to breit frame
+ Lorentz5Momentum phadron = hadron->momentum();
+ phadron.setMass(0.*GeV);
+ phadron.rescaleRho();
+ Lorentz5Momentum pb = quark[0]->momentum();
+ Lorentz5Momentum pbasis = phadron;
+ Axis axis(q_.vect().unit());
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ LorentzRotation rot_ = LorentzRotation();
+ if(axis.perp2()>1e-20) {
+ rot_.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ rot_.rotateX(Constants::pi);
+ }
+ if(abs(1.-q_.e()/q_.vect().mag())>1e-6) rot_.boostZ( q_.e()/q_.vect().mag());
+ pb *= rot_;
+ if(pb.perp2()/GeV2>1e-20) {
+ Boost trans = -1./pb.e()*pb.vect();
+ trans.setZ(0.);
+ rot_.boost(trans);
+ }
+ Lorentz5Momentum pl = rot_*lepton[0]->momentum();
+ rot_.rotateZ(-atan2(pl.y(),pl.x()));
+ // momenta of the particles
+ pl_[0]=rot_*lepton[0]->momentum();
+ pl_[1]=rot_*lepton[1]->momentum();
+ pq_[0]=rot_* quark[0]->momentum();
+ pq_[1]=rot_* quark[1]->momentum();
+ q_ *= rot_;
+ // coefficient for the matrix elements
+ acoeff_ = A(lepton[0]->dataPtr(),lepton[1]->dataPtr(),
+ quark [0]->dataPtr(),quark [1]->dataPtr(),q2_);
+ // generate a compton point
+ generateCompton();
+ generateBGF();
+ // no valid emission, return
+ if(pTCompton_<ZERO&&pTBGF_<ZERO) return HardTreePtr();
+ // type of emission, pick highest pT
+ bool isCompton=pTCompton_>pTBGF_;
+// // find the sudakov for the branching
+// SudakovPtr sudakov;
+// // ISR
+// if(ComptonISFS_||!isCompton) {
+// BranchingList branchings=evolver()->splittingGenerator()->initialStateBranchings();
+// long index = abs(partons_[0]->id());
+// IdList br(3);
+// if(isCompton) {
+// br[0] = index;
+// br[1] = index;
+// br[2] = ParticleID::g;
+// }
+// else {
+// br[0] = ParticleID::g;
+// br[1] = abs(partons_[0]->id());
+// br[2] = -abs(partons_[0]->id());
+// }
+// for(BranchingList::const_iterator cit = branchings.lower_bound(index);
+// cit != branchings.upper_bound(index); ++cit ) {
+// IdList ids = cit->second.second;
+// if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
+// sudakov=cit->second.first;
+// break;
+// }
+// }
+// }
+// // FSR
+// else {
+// BranchingList branchings =
+// evolver()->splittingGenerator()->finalStateBranchings();
+// long index=abs(partons_[1]->id());
+// for(BranchingList::const_iterator cit = branchings.lower_bound(index);
+// cit != branchings.upper_bound(index); ++cit ) {
+// IdList ids = cit->second.second;
+// if(ids[0]==index&&ids[1]==index&&ids[2]==ParticleID::g) {
+// sudakov = cit->second.first;
+// break;
+// }
+// }
+// }
+// if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
+// << "DISBase::generateHardest()"
+// << Exception::runerror;
+ // add the leptons
+ vector<HardBranchingPtr> spaceBranchings,allBranchings;
+ spaceBranchings.push_back(new_ptr(HardBranching(lepton[0],SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ allBranchings.push_back(spaceBranchings.back());
+ allBranchings.push_back(new_ptr(HardBranching(lepton[1],SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ // compton hardest
+ if(isCompton) {
+ rot_.invert();
+ for(unsigned int ix=0;ix<ComptonMomenta_.size();++ix) {
+ ComptonMomenta_[ix].transform(rot_);
+ }
+ ShowerParticlePtr newqout (new_ptr(ShowerParticle(partons_[1],true)));
+ newqout->set5Momentum(ComptonMomenta_[1]);
+ ShowerParticlePtr newg(new_ptr(ShowerParticle(gluon_,true)));
+ newg->set5Momentum(ComptonMomenta_[2]);
+ ShowerParticlePtr newqin (new_ptr(ShowerParticle(partons_[0],false )));
+ newqin->set5Momentum(ComptonMomenta_[0]);
+ if(ComptonISFS_) {
+ ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[0],false)));
+ newspace->set5Momentum(ComptonMomenta_[0]-ComptonMomenta_[2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),
+ spaceBranch,
+ HardBranching::Incoming)));
+ spaceBranch->addChild(offBranch);
+ HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),spaceBranch,
+ HardBranching::Outgoing)));
+ spaceBranch->addChild(g);
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ allBranchings.push_back(outBranch);
+ ColinePtr newin(new_ptr(ColourLine())),newout(new_ptr(ColourLine()));
+ newin->addColoured(newqin,partons_[0]->id()<0);
+ newin->addColoured(newg ,partons_[0]->id()<0);
+ newout->addColoured(newspace,partons_[0]->id()<0);
+ newout->addColoured(newqout,partons_[1]->id()<0);
+ newout->addColoured(newg ,partons_[1]->id()>0);
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newqout ,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ else {
+ ShowerParticlePtr newtime(new_ptr(ShowerParticle(partons_[1],true)));
+ newtime->set5Momentum(ComptonMomenta_[1]+ComptonMomenta_[2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newqin,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newtime,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ HardBranchingPtr g(new_ptr(HardBranching(newg,SudakovPtr(),offBranch,
+ HardBranching::Outgoing)));
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newqout,SudakovPtr(),offBranch,
+ HardBranching::Outgoing)));
+ offBranch->addChild(outBranch);
+ offBranch->addChild(g);
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newqin ,newqin->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newtime,newqin->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ }
+ // BGF hardest
+ else {
+ rot_.invert();
+ for(unsigned int ix=0;ix<BGFMomenta_.size();++ix) {
+ BGFMomenta_[ix].transform(rot_);
+ }
+ ShowerParticlePtr newq (new_ptr(ShowerParticle(partons_[1],true)));
+ newq->set5Momentum(BGFMomenta_[1]);
+ ShowerParticlePtr newqbar(new_ptr(ShowerParticle(partons_[0]->CC(),true)));
+ newqbar->set5Momentum(BGFMomenta_[2]);
+ ShowerParticlePtr newg (new_ptr(ShowerParticle(gluon_,false)));
+ newg->set5Momentum(BGFMomenta_[0]);
+ ShowerParticlePtr newspace(new_ptr(ShowerParticle(partons_[0],false)));
+ newspace->set5Momentum(BGFMomenta_[0]-BGFMomenta_[2]);
+ HardBranchingPtr spaceBranch(new_ptr(HardBranching(newg,SudakovPtr(),HardBranchingPtr(),
+ HardBranching::Incoming)));
+ HardBranchingPtr offBranch(new_ptr(HardBranching(newspace,SudakovPtr(),spaceBranch,
+ HardBranching::Incoming)));
+ HardBranchingPtr qbar(new_ptr(HardBranching(newqbar,SudakovPtr(),spaceBranch,
+ HardBranching::Outgoing)));
+ spaceBranch->addChild(offBranch);
+ spaceBranch->addChild(qbar);
+ HardBranchingPtr outBranch(new_ptr(HardBranching(newq,SudakovPtr(),
+ HardBranchingPtr(),
+ HardBranching::Outgoing)));
+ spaceBranchings.push_back(spaceBranch);
+ allBranchings.push_back(offBranch);
+ allBranchings.push_back(outBranch);
+ ColinePtr newline(new_ptr(ColourLine()));
+ newline->addColoured(newspace,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ newline->addColoured(newq ,newspace->dataPtr()->iColour()!=PDT::Colour3);
+ }
+ HardTreePtr newTree(new_ptr(HardTree(allBranchings,spaceBranchings,
+ ShowerInteraction::QCD)));
+ // Set the maximum pt for all other emissions and connect hard and shower tree
+ Energy pT = isCompton ? pTCompton_ : pTBGF_;
+ // incoming particles
+ for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
+ cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
+ // set maximum pT
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pT);
+ for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+ cjt!=newTree->branchings().end();++cjt) {
+ if(!(*cjt)->branchingParticle()->isFinalState()&&
+ (*cjt)->branchingParticle()->id()==cit->first->progenitor()->id()) {
+ newTree->connect(cit->first->progenitor(),*cjt);
+ tPPtr beam =cit->first->original();
+ if(!beam->parents().empty()) beam=beam->parents()[0];
+ (*cjt)->beam(beam);
+ HardBranchingPtr parent=(*cjt)->parent();
+ while(parent) {
+ parent->beam(beam);
+ parent=parent->parent();
+ };
+ }
+ }
+ }
+ // outgoing particles
+ for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
+ cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
+ // set maximum pT
+ if(QuarkMatcher::Check(cit->first->progenitor()->data()))
+ cit->first->maximumpT(pT);
+ for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+ cjt!=newTree->branchings().end();++cjt) {
+ if((*cjt)->branchingParticle()->isFinalState()&&
+ (*cjt)->branchingParticle()->id()==cit->first->progenitor()->id()) {
+ newTree->connect(cit->first->progenitor(),*cjt);
+ }
+ }
+ }
+ // set the evolution partners and scales
+ ShowerParticleVector particles;
+ for(set<HardBranchingPtr>::iterator cit=newTree->branchings().begin();
+ cit!=newTree->branchings().end();++cit) {
+ particles.push_back((*cit)->branchingParticle());
+ }
+// evolver()->showerModel()->partnerFinder()->
+// setInitialEvolutionScales(particles,true,ShowerInteraction::QCD,true);
+// // Calculate the shower variables:
+// evolver()->showerModel()->kinematicsReconstructor()->
+// deconstructHardJets(newTree,evolver(),ShowerInteraction::QCD);
+// for(set<HardBranchingPtr>::iterator cjt=newTree->branchings().begin();
+// cjt!=newTree->branchings().end();++cjt) {
+// if(cjt==newTree->branchings().begin()) {
+// (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
+// ++cjt;
+// (**cjt).showerMomentum((**cjt).branchingParticle()->momentum());
+// ++cjt;
+// }
+// // incoming
+// if((**cjt).status()==HardBranching::Incoming) {
+// quark[0]->set5Momentum((**cjt).showerMomentum());
+// }
+// // outgoing
+// else {
+// quark[1]->set5Momentum((**cjt).showerMomentum());
+// }
+// }
+ return newTree;
+}
+
+void DISBase::generateCompton() {
+ // maximum value of the xT
+ double xT = sqrt((1.-xB_)/xB_);
+ double xTMin = 2.*pTmin_/sqrt(q2_);
+ double zp;
+ // prefactor
+ double a = alpha_->overestimateValue()*comptonWeight_/Constants::twopi;
+ // loop to generate kinematics
+ double wgt(0.),xp(0.);
+ vector<double> azicoeff;
+ do {
+ wgt = 0.;
+ // intergration variables dxT/xT^3
+ xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
+ // zp
+ zp = UseRandom::rnd();
+ xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
+ // check allowed
+ if(xp<xB_||xp>1.) continue;
+ // phase-space piece of the weight
+ wgt = 8.*(1.-xp)*zp/comptonWeight_;
+ // PDF piece of the weight
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
+ wgt *= pdf_->xfx(beam_,partons_[0],scale,xB_/xp)/
+ pdf_->xfx(beam_,partons_[0],q2_ ,xB_);
+ // me piece of the weight
+ double x2 = 1.-(1.-zp)/xp;
+ azicoeff = ComptonME(xp,x2,xT,false);
+ wgt *= 4./3.*alpha_->ratio(0.25*q2_*sqr(xT))*(azicoeff[0]+0.5*azicoeff[2]);
+ if(wgt>1.||wgt<0.) {
+ ostringstream wstring;
+ wstring << "DISBase::generateCompton() "
+ << "Weight greater than one or less than zero"
+ << "wgt = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(xT>xTMin&&UseRandom::rnd()>wgt);
+ if(xT<=xTMin) {
+ pTCompton_=-GeV;
+ return;
+ }
+ // generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi));
+ phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
+ << "::generateCompton() to"
+ << " generate phi" << Exception::eventerror;
+ // momenta for the configuration
+ Energy Q(sqrt(q2_));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ pTCompton_ = 0.5*Q*xT;
+ ComptonMomenta_.resize(3);
+ ComptonMomenta_[0] = p0;
+ ComptonMomenta_[1] = p1;
+ ComptonMomenta_[2] = p2;
+ ComptonISFS_ = zp>xp;
+// debuggingMatrixElement(false,p0,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
+// leptons_[0],leptons_[1],
+// partons_[0],partons_[1],
+// q2_,phi,x2,x3,xT,zp,xp,azicoeff,false);
+}
+
+void DISBase::generateBGF() {
+ // maximum value of the xT
+ double xT = (1.-xB_)/xB_;
+ double xTMin = 2.*max(pTmin_,pTCompton_)/sqrt(q2_);
+ double zp;
+ // prefactor
+ double a = alpha_->overestimateValue()*BGFWeight_/Constants::twopi;
+ // loop to generate kinematics
+ double wgt(0.),xp(0.);
+ vector<double> azicoeff;
+ do {
+ wgt = 0.;
+ // intergration variables dxT/xT^3
+ xT *= 1./sqrt(1.-2.*log(UseRandom::rnd())/a*sqr(xT));
+ // zp
+ zp = UseRandom::rnd();
+ xp = 1./(1.+0.25*sqr(xT)/zp/(1.-zp));
+ // check allowed
+ if(xp<xB_||xp>1.) continue;
+ // phase-space piece of the weight
+ wgt = 8.*sqr(1.-xp)*zp/BGFWeight_;
+ // PDF piece of the weight
+ Energy2 scale = q2_*((1.-xp)*(1-zp)*zp/xp+1.);
+ wgt *= pdf_->xfx(beam_,gluon_ ,scale,xB_/xp)/
+ pdf_->xfx(beam_,partons_[0],q2_ ,xB_);
+ // me piece of the weight
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ azicoeff = BGFME(xp,x2,x3,xT,false);
+ wgt *= 0.5*alpha_->ratio(0.25*q2_*sqr(xT))*
+ (azicoeff[0]+0.5*azicoeff[2]);
+ if(wgt>1.||wgt<0.) {
+ ostringstream wstring;
+ wstring << "DISBase::generateBGF() "
+ << "Weight greater than one or less than zero"
+ << "wgt = " << wgt << "\n";
+ generator()->logWarning( Exception(wstring.str(),
+ Exception::warning) );
+ }
+ }
+ while(xT>xTMin&&UseRandom::rnd()>wgt);
+ if(xT<=xTMin) {
+ pTBGF_=-GeV;
+ return;
+ }
+ // generate phi
+ unsigned int itry(0);
+ double phimax = std::accumulate(azicoeff.begin(),azicoeff.end(),0.);
+ double phiwgt,phi;
+ do {
+ phi = UseRandom::rnd()*Constants::twopi;
+ double cphi(cos(phi));
+ phiwgt = azicoeff[0]+azicoeff[1]*cphi+azicoeff[2]*sqr(cphi);
+ ++itry;
+ }
+ while (phimax*UseRandom::rnd() > phiwgt && itry<200);
+ if(itry==200) throw Exception() << "Too many tries in DISMECorrection"
+ << "::generateBGF() to"
+ << " generate phi" << Exception::eventerror;
+ // momenta for the configuration
+ Energy Q(sqrt(q2_));
+ double x1 = -1./xp;
+ double x2 = 1.-(1.-zp)/xp;
+ double x3 = 2.+x1-x2;
+ Lorentz5Momentum p1( 0.5*Q*xT*cos(phi), 0.5*Q*xT*sin(phi),
+ -0.5*Q*x2, 0.5*Q*sqrt(sqr(xT)+sqr(x2)));
+ Lorentz5Momentum p2(-0.5*Q*xT*cos(phi), -0.5*Q*xT*sin(phi),
+ -0.5*Q*x3, 0.5*Q*sqrt(sqr(xT)+sqr(x3)));
+ Lorentz5Momentum p0(ZERO,ZERO,-0.5*Q*x1,-0.5*Q*x1);
+ pTBGF_=0.5*Q*xT;
+ BGFMomenta_.resize(3);
+ BGFMomenta_[0]=p0;
+ BGFMomenta_[1]=p1;
+ BGFMomenta_[2]=p2;
+// debuggingMatrixElement(true,p0,p1,p2,gluon_,pl_[0],pl_[1],pq_[0],pq_[1],
+// leptons_[0],leptons_[1],
+// partons_[0],partons_[1],
+// q2_,phi,x2,x3,xT,zp,xp,azicoeff,false);
+}
+
+int DISBase::nDim() const {
+ return HwMEBase::nDim() + (contrib_>0 ? 1 : 0 );
+}
+
+bool DISBase::generateKinematics(const double * r) {
+ // Born kinematics
+ if(!HwMEBase::generateKinematics(r)) return false;
+ if(contrib_!=0) {
+ // hadron and momentum fraction
+ if(HadronMatcher::Check(*lastParticles().first->dataPtr())) {
+ hadron_ = dynamic_ptr_cast<tcBeamPtr>(lastParticles().first->dataPtr());
+ xB_ = lastX1();
+ }
+ else {
+ hadron_ = dynamic_ptr_cast<tcBeamPtr>(lastParticles().second->dataPtr());
+ xB_ = lastX2();
+ }
+ // Q2
+ q2_ = -(meMomenta()[0]-meMomenta()[2]).m2();
+ // xp
+ int ndim=nDim();
+ double rhomin = pow(1.-xB_,1.-power_);
+ double rho = r[ndim-1]*rhomin;
+ xp_ = 1.-pow(rho,1./(1.-power_));
+ jac_ = rhomin/(1.-power_)*pow(1.-xp_,power_);
+ jacobian(jacobian()*jac_);
+ }
+ return true;
+}
+
+Energy2 DISBase::scale() const {
+ return scaleOpt_ == 1 ?
+ -sqr(scaleFact_)*tHat() : sqr(scaleFact_*muF_);
+}
+
+CrossSection DISBase::dSigHatDR() const {
+ return NLOWeight()*HwMEBase::dSigHatDR();
+}
+
+double DISBase::NLOWeight() const {
+ // If only leading order is required return 1:
+ if(contrib_==0) return 1.;
+ // scale and prefactors
+ Energy2 mu2(scale());
+ double aS = SM().alphaS(mu2);
+ double CFfact = 4./3.*aS/Constants::twopi;
+ double TRfact = 1./2.*aS/Constants::twopi;
+ // LO + dipole subtracted virtual + collinear quark bit with LO pdf
+ double virt = 1.+CFfact*(-4.5-1./3.*sqr(Constants::pi)+1.5*log(q2_/mu2/(1.-xB_))
+ +2.*log(1.-xB_)*log(q2_/mu2)+sqr(log(1.-xB_)));
+ virt /= jac_;
+ // PDF from leading-order
+ double loPDF = hadron_->pdf()->xfx(hadron_,mePartonData()[1],mu2,xB_)/xB_;
+ // NLO gluon PDF
+ tcPDPtr gluon = getParticleData(ParticleID::g);
+ double gPDF = hadron_->pdf()->xfx(hadron_,gluon,mu2,xB_/xp_)*xp_/xB_;
+ // NLO quark PDF
+ double qPDF = hadron_->pdf()->xfx(hadron_,mePartonData()[1],mu2,xB_/xp_)*xp_/xB_;
+ // collinear counterterms
+ // gluon
+ double collg =
+ TRfact/xp_*gPDF*(2.*xp_*(1.-xp_)+(sqr(xp_)+sqr(1.-xp_))*log((1.-xp_)*q2_/xp_/mu2));
+ // quark
+ double collq =
+ CFfact/xp_*qPDF*(1-xp_-2./(1.-xp_)*log(xp_)-(1.+xp_)*log((1.-xp_)/xp_*q2_/mu2))+
+ CFfact/xp_*(qPDF-xp_*loPDF)*(2./(1.-xp_)*log(q2_*(1.-xp_)/mu2)-1.5/(1.-xp_));
+ // calculate the A coefficient for the real pieces
+ double a(A(mePartonData()[0],mePartonData()[2],
+ mePartonData()[1],mePartonData()[3],q2_));
+ // cacluate lepton kinematic variables
+ Lorentz5Momentum q = meMomenta()[0]-meMomenta()[2];
+ double yB = (q*meMomenta()[1])/(meMomenta()[0]*meMomenta()[1]);
+ double l = 2./yB-1.;
+ // q -> qg term
+ double realq = CFfact/xp_/(1.+a*l+sqr(l))*qPDF/loPDF*
+ (2.+2.*sqr(l)-xp_+3.*xp_*sqr(l)+a*l*(2.*xp_+1.));
+ // g -> q qbar term
+ double realg =-TRfact/xp_/(1.+a*l+sqr(l))*gPDF/loPDF*
+ ((1.+sqr(l)+2.*(1.-3.*sqr(l))*xp_*(1.-xp_))
+ +2.*a*l*(1.-2.*xp_*(1.-xp_)));
+ // return the full result
+ double wgt = virt+((collq+collg)/loPDF+realq+realg);
+ // double f2g = gPDF/xp_*TRfact*((sqr(1-xp_)+sqr(xp_))*log((1-xp_)/xp_)+
+ // 8*xp_*(1.-xp_)-1.);
+ // double f2q =
+ // loPDF/jac_*(1.+CFfact*(-1.5*log(1.-xB_)+sqr(log(1.-xB_))
+ // -sqr(Constants::pi)/3.-4.5))
+ // +qPDF *CFfact/xp_*(3.+2.*xp_-(1.+xp_)*log(1.-xp_)
+ // -(1.+sqr(xp_))/(1.-xp_)*log(xp_))
+ // +(qPDF-xp_*loPDF)*CFfact/xp_*(2.*log(1.-xp_)/(1.-xp_)-1.5/(1.-xp_));
+ // double wgt = (f2g+f2q)/loPDF;
+ return contrib_ == 1 ? max(0.,wgt) : max(0.,-wgt);
+}
diff --git a/MatrixElement/DIS/DISBase.h b/MatrixElement/DIS/DISBase.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/DIS/DISBase.h
@@ -0,0 +1,467 @@
+// -*- C++ -*-
+#ifndef HERWIG_DISBase_H
+#define HERWIG_DISBase_H
+//
+// This is the declaration of the DISBase class.
+//
+
+#include "Herwig++/MatrixElement/HwMEBase.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The DISBase class is the base class for the implementation
+ * of DIS type processes including corrections in both the old
+ * fashioned matrix element and POWHEG approaches
+ *
+ * @see \ref DISBaseInterfaces "The interfaces"
+ * defined for DISBase.
+ */
+class DISBase: public HwMEBase {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ DISBase();
+
+ /**
+ * The default constructor.
+ */
+ virtual ~DISBase();
+
+ /**
+ * Members for the old-fashioned matrix element correction
+ */
+ //@{
+ /**
+ * Has an old fashioned ME correction
+ */
+ virtual bool hasMECorrection() {return true;}
+
+ /**
+ * Initialize the ME correction
+ */
+ virtual void initializeMECorrection(ShowerTreePtr, double &,
+ double & );
+
+ /**
+ * Apply the hard matrix element correction to a given hard process or decay
+ */
+ virtual void applyHardMatrixElementCorrection(ShowerTreePtr);
+
+ /**
+ * Apply the soft matrix element correction
+ * @param initial The particle from the hard process which started the
+ * shower
+ * @param parent The initial particle in the current branching
+ * @param br The branching struct
+ * @return If true the emission should be vetoed
+ */
+ virtual bool softMatrixElementVeto(ShowerProgenitorPtr,
+ ShowerParticlePtr,Branching);
+ //@}
+
+ /**
+ * Members for the POWHEG stype correction
+ */
+ //@{
+ /**
+ * Has a POWHEG style correction
+ */
+ virtual bool hasPOWHEGCorrection() {return true;}
+
+ /**
+ * Apply the POWHEG style correction
+ */
+ virtual HardTreePtr generateHardest(ShowerTreePtr);
+ //@}
+
+public:
+
+ /** @name Virtual functions required by the MEBase class. */
+ //@{
+ /**
+ * Return the scale associated with the last set phase space point.
+ */
+ virtual Energy2 scale() const;
+
+ /**
+ * The number of internal degrees of freedom used in the matrix
+ * element.
+ */
+ virtual int nDim() const;
+
+ /**
+ * Generate internal degrees of freedom given nDim() uniform
+ * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
+ * generator, the dSigHatDR should be a smooth function of these
+ * numbers, although this is not strictly necessary.
+ * @param r a pointer to the first of nDim() consecutive random numbers.
+ * @return true if the generation succeeded, otherwise false.
+ */
+ virtual bool generateKinematics(const double * r);
+
+ /**
+ * Return the matrix element squared differential in the variables
+ * given by the last call to generateKinematics().
+ */
+ virtual CrossSection dSigHatDR() const;
+ //@}
+
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ * @throws InitException if object could not be initialized properly.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is an abstract class with persistent data.
+ */
+ static AbstractClassDescription<DISBase> initDISBase;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ DISBase & operator=(const DISBase &);
+
+protected:
+
+ /**
+ * The NLO weight
+ */
+ double NLOWeight() const;
+
+ /**
+ * Calculate the coefficient A for the correlations
+ */
+ virtual double A(tcPDPtr lin, tcPDPtr lout, tcPDPtr qin, tcPDPtr qout,
+ Energy2 scale) const =0;
+
+ /**
+ * Members for the matrix element correction
+ */
+ //@{
+ /**
+ * Generate the values of \f$x_p\f$ and \f$z_p\f$
+ * @param xp The value of xp, output
+ * @param zp The value of zp, output
+ */
+ double generateComptonPoint(double &xp, double & zp);
+
+ /**
+ * Generate the values of \f$x_p\f$ and \f$z_p\f$
+ * @param xp The value of xp, output
+ * @param zp The value of zp, output
+ */
+ double generateBGFPoint(double &xp, double & zp);
+
+ /**
+ * Return the coefficients for the matrix element piece for
+ * the QCD compton case. The output is the \f$a_i\f$ coefficients to
+ * give the function as
+ * \f$a_0+a_1\cos\phi+a_2\sin\phi+a_3\cos^2\phi+a_4\sin^2\phi\f$
+ * @param xp \f$x_p\f$
+ * @param x2 \f$x_2\f$
+ * @param xperp \f$x_\perp\f$
+ * @param norm Normalise to the large $l$ value of the ME
+ */
+ vector<double> ComptonME(double xp, double x2, double xperp,
+ bool norm);
+
+ /**
+ * Return the coefficients for the matrix element piece for
+ * the QCD compton case. The output is the \f$a_i\f$ coefficients to
+ * give the function as
+ * \f$a_0+a_1\cos\phi+a_2\sin\phi+a_3\cos^2\phi+a_4\sin^2\phi\f$
+ * @param xp \f$x_p\f$
+ * @param x2 \f$x_3\f$
+ * @param x3 \f$x_2\f$
+ * @param xperp \f$x_\perp\f$
+ * @param norm Normalise to the large $l$ value of the ME
+ */
+ vector<double> BGFME(double xp, double x2, double x3, double xperp,
+ bool norm);
+ //@}
+
+ /**
+ * Members for the POWHEG correction
+ */
+ //@{
+ /**
+ * Generate a Compton process
+ */
+ void generateCompton();
+
+ /**
+ * Generate a BGF process
+ */
+ void generateBGF();
+ //@}
+
+private:
+
+ /**
+ * Parameters for the matrix element correction
+ */
+ //@{
+ /**
+ * Enchancement factor for ISR
+ */
+ double initial_;
+
+ /**
+ * Enchancement factor for FSR
+ */
+ double final_;
+
+ /**
+ * Relative fraction of compton and BGF processes to generate
+ */
+ double procProb_;
+
+ /**
+ * Integral for compton process
+ */
+ double comptonInt_;
+
+ /**
+ * Integral for BGF process
+ */
+ double bgfInt_;
+ //@}
+
+ /**
+ * Parameters for the POWHEG correction
+ */
+ //@{
+ /**
+ * Weight for the compton channel
+ */
+ double comptonWeight_;
+
+ /**
+ * Weight for the BGF channel
+ */
+ double BGFWeight_;
+
+ /**
+ * Minimum value of \f$p_T\f$
+ */
+ Energy pTmin_;
+ //@}
+
+ /**
+ * Parameters for the point being generated
+ */
+ //@{
+ /**
+ * \f$Q^2\f$
+ */
+ Energy2 q2_;
+
+ /**
+ *
+ */
+ double l_;
+
+ /**
+ * Borm momentum fraction
+ */
+ double xB_;
+
+ /**
+ * Beam particle
+ */
+ tcBeamPtr beam_;
+
+ /**
+ * Partons
+ */
+ tcPDPtr partons_[2];
+
+ /**
+ * Leptons
+ */
+ tcPDPtr leptons_[2];
+
+ /**
+ * PDF object
+ */
+ tcPDFPtr pdf_;
+ /**
+ * Rotation to the Breit frame
+ */
+ LorentzRotation rot_;
+
+ /**
+ * Lepton momenta
+ */
+ Lorentz5Momentum pl_[2];
+
+ /**
+ * Quark momenta
+ */
+ Lorentz5Momentum pq_[2];
+
+ /**
+ * q
+ */
+ Lorentz5Momentum q_;
+
+ /**
+ * Compton parameters
+ */
+ Energy pTCompton_;
+ bool ComptonISFS_;
+ vector<Lorentz5Momentum> ComptonMomenta_;
+
+ /**
+ * BGF parameters
+ */
+ Energy pTBGF_;
+ vector<Lorentz5Momentum> BGFMomenta_;
+ //@}
+
+ /**
+ * The coefficient for the correlations
+ */
+ double acoeff_;
+
+ /**
+ * Coupling
+ */
+ ShowerAlphaPtr alpha_;
+
+ /**
+ * Gluon particle data object
+ */
+ PDPtr gluon_;
+
+private:
+
+ /**
+ * The radiative variables
+ */
+ //@{
+ /**
+ * The \f$x_p\f$ or \f$z\f$ real integration variable
+ */
+ double xp_;
+ //@}
+
+ /**
+ * The hadron
+ */
+ tcBeamPtr hadron_;
+
+ /**
+ * Selects a dynamic or fixed factorization scale
+ */
+ unsigned int scaleOpt_;
+
+ /**
+ * The factorization scale
+ */
+ Energy muF_;
+
+ /**
+ * Prefactor if variable scale used
+ */
+ double scaleFact_;
+
+ /**
+ * Whether to generate the positive, negative or leading order contribution
+ */
+ unsigned int contrib_;
+
+ /**
+ * Power for sampling \f$x_p\f$
+ */
+ double power_;
+
+ /**
+ * Jacobian for \f$x_p\f$ integral
+ */
+ double jac_;
+
+};
+
+}
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of DISBase. */
+template <>
+struct BaseClassTrait<Herwig::DISBase,1> {
+ /** Typedef of the first base class of DISBase. */
+ typedef Herwig::HwMEBase NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the DISBase class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::DISBase>
+ : public ClassTraitsBase<Herwig::DISBase> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::DISBase"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * MENeutralCurrentDIS is implemented. It may also include several, space-separated,
+ * libraries if the class MENeutralCurrentDIS depends on other classes (base classes
+ * excepted). In this case the listed libraries will be dynamically
+ * linked in the order they are specified.
+ */
+ static string library() { return "HwMEDIS.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_DISBase_H */
diff --git a/MatrixElement/DIS/MEChargedCurrentDIS.cc b/MatrixElement/DIS/MEChargedCurrentDIS.cc
--- a/MatrixElement/DIS/MEChargedCurrentDIS.cc
+++ b/MatrixElement/DIS/MEChargedCurrentDIS.cc
@@ -1,297 +1,301 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEChargedCurrentDIS class.
//
#include "MEChargedCurrentDIS.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "Herwig++/MatrixElement/HardVertex.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
using namespace Herwig;
MEChargedCurrentDIS::MEChargedCurrentDIS()
: _maxflavour(5), _massopt(0) {
vector<unsigned int> mopt(2,1);
mopt[1] = _massopt;
massOption(mopt);
}
void MEChargedCurrentDIS::doinit() {
- HwMEBase::doinit();
+ DISBase::doinit();
_wp = getParticleData(ThePEG::ParticleID::Wplus );
_wm = getParticleData(ThePEG::ParticleID::Wminus);
// cast the SM pointer to the Herwig SM pointer
tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException()
<< "Must be the Herwig++ StandardModel class in "
<< "MEChargedCurrentDIS::doinit" << Exception::abortnow;
// vertices
_theFFWVertex = hwsm->vertexFFW();
}
void MEChargedCurrentDIS::getDiagrams() const {
// possible quarks
typedef std::vector<pair<long,long> > Pairvector;
Pairvector quarkpair;
quarkpair.reserve(6);
// don't even think of putting 'break' in here!
switch(_maxflavour) {
case 6:
quarkpair.push_back(make_pair(ParticleID::s, ParticleID::t));
quarkpair.push_back(make_pair(ParticleID::d, ParticleID::t));
quarkpair.push_back(make_pair(ParticleID::b, ParticleID::t));
case 5:
quarkpair.push_back(make_pair(ParticleID::b, ParticleID::c));
quarkpair.push_back(make_pair(ParticleID::b, ParticleID::u));
case 4:
quarkpair.push_back(make_pair(ParticleID::s, ParticleID::c));
quarkpair.push_back(make_pair(ParticleID::d, ParticleID::c));
case 3:
quarkpair.push_back(make_pair(ParticleID::s, ParticleID::u));
case 2:
quarkpair.push_back(make_pair(ParticleID::d, ParticleID::u));
default:
;
}
// create the diagrams
for(int il1=11;il1<=14;++il1) {
int il2 = il1%2==0 ? il1-1 : il1+1;
for(unsigned int iz=0;iz<2;++iz) {
tcPDPtr lepin = iz==1 ? getParticleData(il1) : getParticleData(-il1);
tcPDPtr lepout = iz==1 ? getParticleData(il2) : getParticleData(-il2);
tcPDPtr inter = lepin->iCharge()-lepout->iCharge()==3 ? _wp : _wm;
for(unsigned int iq=0;iq<quarkpair.size();++iq) {
tcPDPtr first = getParticleData(quarkpair[iq].first );
tcPDPtr second = getParticleData(quarkpair[iq].second);
if(inter==_wp) {
add(new_ptr((Tree2toNDiagram(3), lepin, inter, first ,
1, lepout, 2, second , -1)));
add(new_ptr((Tree2toNDiagram(3), lepin, inter, second->CC(),
1, lepout, 2, first->CC(), -2)));
}
else {
add(new_ptr((Tree2toNDiagram(3), lepin, inter, second ,
1, lepout, 2, first , -1)));
add(new_ptr((Tree2toNDiagram(3), lepin, inter, first->CC(),
1, lepout, 2, second->CC(), -2)));
}
}
}
}
}
-Energy2 MEChargedCurrentDIS::scale() const {
- return -tHat();
-}
-
unsigned int MEChargedCurrentDIS::orderInAlphaS() const {
return 0;
}
unsigned int MEChargedCurrentDIS::orderInAlphaEW() const {
return 2;
}
Selector<const ColourLines *>
MEChargedCurrentDIS::colourGeometries(tcDiagPtr diag) const {
static ColourLines c1("3 5");
static ColourLines c2("-3 -5");
Selector<const ColourLines *> sel;
if ( diag->id() == -1 )
sel.insert(1.0, &c1);
else
sel.insert(1.0, &c2);
return sel;
}
void MEChargedCurrentDIS::persistentOutput(PersistentOStream & os) const {
os << _theFFWVertex << _maxflavour << _wp << _wm << _massopt;
}
void MEChargedCurrentDIS::persistentInput(PersistentIStream & is, int) {
is >> _theFFWVertex >> _maxflavour >> _wp >> _wm >> _massopt;
}
ClassDescription<MEChargedCurrentDIS> MEChargedCurrentDIS::initMEChargedCurrentDIS;
// Definition of the static class description member.
void MEChargedCurrentDIS::Init() {
static ClassDocumentation<MEChargedCurrentDIS> documentation
("The MEChargedCurrentDIS class implements the matrix elements "
"for leading-order charged current deep inelastic scattering");
static Parameter<MEChargedCurrentDIS,unsigned int> interfaceMaxFlavour
( "MaxFlavour",
"The heaviest incoming quark flavour this matrix element is allowed to handle "
"(if applicable).",
&MEChargedCurrentDIS::_maxflavour, 5, 2, 6, false, false, true);
static Switch<MEChargedCurrentDIS,unsigned int> interfaceMassOption
("MassOption",
"Option for the treatment of the mass of the outgoing quarks",
&MEChargedCurrentDIS::_massopt, 0, false, false);
static SwitchOption interfaceMassOptionMassless
(interfaceMassOption,
"Massless",
"Treat the outgoing quarks as massless",
0);
static SwitchOption interfaceMassOptionMassive
(interfaceMassOption,
"Massive",
"Treat the outgoing quarks as massive",
1);
}
Selector<MEBase::DiagramIndex>
MEChargedCurrentDIS::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) sel.insert(1., i);
return sel;
}
double MEChargedCurrentDIS::helicityME(vector<SpinorWaveFunction> & f1,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder, bool calc) const {
// scale
Energy2 mb2(scale());
// matrix element to be stored
ProductionMatrixElement menew(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
// pick a W boson
tcPDPtr ipart = (mePartonData()[0]->iCharge()-mePartonData()[1]->iCharge())==3 ?
_wp : _wm;
// declare the variables we need
VectorWaveFunction inter;
double me(0.);
Complex diag;
// sum over helicities to get the matrix element
unsigned int hel[4];
unsigned int lhel1,lhel2,qhel1,qhel2;
for(lhel1=0;lhel1<2;++lhel1) {
for(lhel2=0;lhel2<2;++lhel2) {
// intermediate W
inter = _theFFWVertex->evaluate(mb2,3,ipart,f1[lhel1],a1[lhel2]);
for(qhel1=0;qhel1<2;++qhel1) {
for(qhel2=0;qhel2<2;++qhel2) {
hel[0] = lhel1;
hel[1] = qhel1;
hel[2] = lhel2;
hel[3] = qhel2;
if(!lorder) swap(hel[0],hel[2]);
if(!qorder) swap(hel[1],hel[3]);
diag = _theFFWVertex->evaluate(mb2,f2[qhel1],a2[qhel2],inter);
me += norm(diag);
menew(hel[0],hel[1],hel[2],hel[3]) = diag;
}
}
}
}
// spin and colour factor
me *= 0.25;
tcPolarizedBeamPDPtr beam[2] =
{dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[0]),
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[1])};
if( beam[0] || beam[1] ) {
RhoDMatrix rho[2] = {beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()),
beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())};
me = menew.average(rho[0],rho[1]);
}
if(calc) _me.reset(menew);
return me;
}
double MEChargedCurrentDIS::me2() const {
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
bool lorder,qorder;
SpinorWaveFunction l1,q1;
SpinorBarWaveFunction l2,q2;
// lepton wave functions
if(mePartonData()[0]->id()>0) {
lorder=true;
l1 = SpinorWaveFunction (meMomenta()[0],mePartonData()[0],incoming);
l2 = SpinorBarWaveFunction(meMomenta()[2],mePartonData()[2],outgoing);
}
else {
lorder=false;
l1 = SpinorWaveFunction (meMomenta()[2],mePartonData()[2],outgoing);
l2 = SpinorBarWaveFunction(meMomenta()[0],mePartonData()[0],incoming);
}
// quark wave functions
if(mePartonData()[1]->id()>0) {
qorder = true;
q1 = SpinorWaveFunction (meMomenta()[1],mePartonData()[1],incoming);
q2 = SpinorBarWaveFunction(meMomenta()[3],mePartonData()[3],outgoing);
}
else {
qorder = false;
q1 = SpinorWaveFunction (meMomenta()[3],mePartonData()[3],outgoing);
q2 = SpinorBarWaveFunction(meMomenta()[1],mePartonData()[1],incoming);
}
// wavefunctions for various helicities
for(unsigned int ix=0;ix<2;++ix) {
l1.reset(ix); f1.push_back(l1);
l2.reset(ix); a1.push_back(l2);
q1.reset(ix); f2.push_back(q1);
q2.reset(ix); a2.push_back(q2);
}
return helicityME(f1,f2,a1,a2,lorder,qorder,false);
}
void MEChargedCurrentDIS::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);
hard.push_back(sub->outgoing()[1]);
unsigned int order[4]={0,1,2,3};
bool lorder(true),qorder(true);
if(hard[0]->id()<0) {
swap(order[0],order[2]);
lorder = false;
}
if(hard[1]->id()<0) {
swap(order[1],order[3]);
qorder = false;
}
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
SpinorWaveFunction (f1,hard[order[0]],incoming,!lorder,true);
SpinorWaveFunction (f2,hard[order[1]],incoming,!qorder,true);
SpinorBarWaveFunction(a1,hard[order[2]],outgoing, lorder,true);
SpinorBarWaveFunction(a2,hard[order[3]],outgoing, qorder,true);
helicityME(f1,f2,a1,a2,lorder,qorder,false);
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix) {
tSpinPtr spin = hard[ix]->spinInfo();
if(ix<2) {
tcPolarizedBeamPDPtr beam =
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(hard[ix]->dataPtr());
if(beam) spin->rhoMatrix() = beam->rhoMatrix();
}
spin->productionVertex(hardvertex);
}
}
+
+double MEChargedCurrentDIS::A(tcPDPtr lin, tcPDPtr,
+ tcPDPtr qin, tcPDPtr, Energy2) const {
+ double output = 2.;
+ if(qin->id()<0) output *= -1.;
+ if(lin->id()<0) output *= -1;
+ return output;
+}
diff --git a/MatrixElement/DIS/MEChargedCurrentDIS.h b/MatrixElement/DIS/MEChargedCurrentDIS.h
--- a/MatrixElement/DIS/MEChargedCurrentDIS.h
+++ b/MatrixElement/DIS/MEChargedCurrentDIS.h
@@ -1,261 +1,263 @@
// -*- C++ -*-
#ifndef HERWIG_MEChargedCurrentDIS_H
#define HERWIG_MEChargedCurrentDIS_H
//
// This is the declaration of the MEChargedCurrentDIS class.
//
-#include "Herwig++/MatrixElement/HwMEBase.h"
+#include "DISBase.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MEChargedCurrentDIS class provides the matrix elements for
* charged current DIS.
*
* By default both the incoming and outgong quarks are assumed to be massless
* although the mass of the outgoing quark can be included if required. This
* option should be used if top production is included.
*
* @see \ref MEChargedCurrentDISInterfaces "The interfaces"
* defined for MEChargedCurrentDIS.
*/
-class MEChargedCurrentDIS: public HwMEBase {
+class MEChargedCurrentDIS: public DISBase {
public:
/**
* The default constructor.
*/
MEChargedCurrentDIS();
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* The matrix element for the kinematical configuration
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
* @return the matrix element scaled with sHat() to give a
* dimensionless number.
*/
virtual double me2() const;
/**
- * Return the scale associated with the last set phase space point.
- */
- virtual Energy2 scale() const;
-
- /**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Get diagram selector. With the information previously supplied with the
* setKinematics method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
* @param dv the diagrams to be weighted.
* @return a Selector relating the given diagrams to their weights.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
* @param diag the diagram chosen.
* @return the possible colour geometries weighted by their
* relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Construct the vertex of spin correlations.
*/
virtual void constructVertex(tSubProPtr);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Matrix element for \f$\ell q\to \gamma/Z \to \ell q\f$.
* @param f1 Fermion on lepton line
* @param a1 Anti-fermion on lepton line
* @param f2 Fermion on quark line
* @param a2 Anti-fermion on quark line
* @param lorder The order of particles on the lepton line
* @param qorder The order of particles on the quark line
* @param me Whether or not to calculate the matrix element for spin correlations
*/
double helicityME(vector<SpinorWaveFunction> & f1 ,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1 ,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder,
bool me) const;
+ /**
+ * Calculate the coefficient A for the correlations in the hard
+ * radiation
+ */
+ virtual double A(tcPDPtr lin, tcPDPtr lout, tcPDPtr qin, tcPDPtr qout,
+ Energy2 scale) const;
+
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MEChargedCurrentDIS> initMEChargedCurrentDIS;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEChargedCurrentDIS & operator=(const MEChargedCurrentDIS &);
private:
/**
* Pointer to the vertex for the helicity calculations
*/
AbstractFFVVertexPtr _theFFWVertex;
/**
* The allowed flavours of the incoming quarks
*/
unsigned int _maxflavour;
/**
* Option for the mass of the outgoing quarks
*/
unsigned int _massopt;
/**
* Matrix element for spin correlations
*/
ProductionMatrixElement _me;
/**
* Pointers to the intermediates resonances
*/
//@{
/**
* Pointer to the \f$W^+\f$
*/
tcPDPtr _wp;
/**
* Pointer to the \f$W^-\f$
*/
tcPDPtr _wm;
//@}
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MEChargedCurrentDIS. */
template <>
struct BaseClassTrait<Herwig::MEChargedCurrentDIS,1> {
/** Typedef of the first base class of MEChargedCurrentDIS. */
- typedef Herwig::HwMEBase NthBase;
+ typedef Herwig::DISBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MEChargedCurrentDIS class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MEChargedCurrentDIS>
: public ClassTraitsBase<Herwig::MEChargedCurrentDIS> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MEChargedCurrentDIS"; }
/**
* The name of a file containing the dynamic library where the class
* MEChargedCurrentDIS is implemented. It may also include several, space-separated,
* libraries if the class MEChargedCurrentDIS depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwMEDIS.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MEChargedCurrentDIS_H */
diff --git a/MatrixElement/DIS/MENeutralCurrentDIS.cc b/MatrixElement/DIS/MENeutralCurrentDIS.cc
--- a/MatrixElement/DIS/MENeutralCurrentDIS.cc
+++ b/MatrixElement/DIS/MENeutralCurrentDIS.cc
@@ -1,323 +1,357 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MENeutralCurrentDIS class.
//
#include "MENeutralCurrentDIS.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig++/MatrixElement/HardVertex.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
using namespace Herwig;
MENeutralCurrentDIS::MENeutralCurrentDIS()
- : _minflavour(1), _maxflavour(5), _gammaZ(0) {
+ : _minflavour(1), _maxflavour(5), _gammaZ(0),
+ _sinW(0.), _cosW(0.), _mz2(ZERO) {
vector<unsigned int> mopt(2,0);
mopt[0] = 1;
massOption(mopt);
}
void MENeutralCurrentDIS::doinit() {
- HwMEBase::doinit();
+ DISBase::doinit();
_z0 = getParticleData(ThePEG::ParticleID::Z0);
_gamma = getParticleData(ThePEG::ParticleID::gamma);
// cast the SM pointer to the Herwig SM pointer
tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException()
<< "Must be the Herwig++ StandardModel class in "
<< "MENeutralCurrentDIS::doinit" << Exception::abortnow;
// vertices
_theFFZVertex = hwsm->vertexFFZ();
_theFFPVertex = hwsm->vertexFFP();
+ // electroweak parameters
+ _sinW = generator()->standardModel()->sin2ThetaW();
+ _cosW = sqrt(1.-_sinW);
+ _sinW = sqrt(_sinW);
+ _mz2 = sqr(_z0->mass());
}
void MENeutralCurrentDIS::getDiagrams() const {
// which intermediates to include
bool gamma = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// create the diagrams
for(int ix=11;ix<=14;++ix) {
for(unsigned int iz=0;iz<2;++iz) {
tPDPtr lep = getParticleData(ix);
if(iz==1) lep = lep->CC();
for(int iy=_minflavour;iy<=_maxflavour;++iy) {
tPDPtr quark = getParticleData(iy);
// lepton quark scattering via gamma and Z
if(gamma) add(new_ptr((Tree2toNDiagram(3), lep, _gamma, quark,
1, lep, 2, quark, -1)));
if(Z0) add(new_ptr((Tree2toNDiagram(3), lep, _z0 , quark,
1, lep, 2, quark, -2)));
// lepton antiquark scattering via gamma and Z
quark = quark->CC();
if(gamma) add(new_ptr((Tree2toNDiagram(3), lep, _gamma, quark,
1, lep, 2, quark, -3)));
if(Z0) add(new_ptr((Tree2toNDiagram(3), lep, _z0 , quark,
1, lep, 2, quark, -4)));
}
}
}
}
-Energy2 MENeutralCurrentDIS::scale() const {
- return -tHat();
-}
-
unsigned int MENeutralCurrentDIS::orderInAlphaS() const {
return 0;
}
unsigned int MENeutralCurrentDIS::orderInAlphaEW() const {
return 2;
}
Selector<const ColourLines *>
MENeutralCurrentDIS::colourGeometries(tcDiagPtr diag) const {
static ColourLines c1("3 5");
static ColourLines c2("-3 -5");
Selector<const ColourLines *> sel;
if ( diag->id() == -1 || diag->id() == -2 )
sel.insert(1.0, &c1);
else
sel.insert(1.0, &c2);
return sel;
}
void MENeutralCurrentDIS::persistentOutput(PersistentOStream & os) const {
os << _minflavour << _maxflavour << _gammaZ << _theFFZVertex << _theFFPVertex
- << _gamma << _z0;
+ << _gamma << _z0 << _sinW << _cosW << ounit(_mz2,GeV2);
}
void MENeutralCurrentDIS::persistentInput(PersistentIStream & is, int) {
is >> _minflavour >> _maxflavour >> _gammaZ >> _theFFZVertex >> _theFFPVertex
- >> _gamma >> _z0;
+ >> _gamma >> _z0 >> _sinW >> _cosW >> iunit(_mz2,GeV2) ;
}
ClassDescription<MENeutralCurrentDIS> MENeutralCurrentDIS::initMENeutralCurrentDIS;
// Definition of the static class description member.
void MENeutralCurrentDIS::Init() {
static ClassDocumentation<MENeutralCurrentDIS> documentation
("The MENeutralCurrentDIS class implements the matrix elements for leading-order "
"neutral current deep inelastic scattering.");
static Parameter<MENeutralCurrentDIS,int> interfaceMaxFlavour
("MaxFlavour",
"The highest incoming quark flavour this matrix element is allowed to handle",
&MENeutralCurrentDIS::_maxflavour, 5, 1, 5,
false, false, Interface::limited);
static Parameter<MENeutralCurrentDIS,int> interfaceMinFlavour
("MinFlavour",
"The lightest incoming quark flavour this matrix element is allowed to handle",
&MENeutralCurrentDIS::_minflavour, 1, 1, 5,
false, false, Interface::limited);
static Switch<MENeutralCurrentDIS,unsigned int> interfaceGammaZ
("GammaZ",
"Which terms to include",
&MENeutralCurrentDIS::_gammaZ, 0, false, false);
static SwitchOption interfaceGammaZAll
(interfaceGammaZ,
"All",
"Include both gamma and Z terms",
0);
static SwitchOption interfaceGammaZGamma
(interfaceGammaZ,
"Gamma",
"Only include the photon",
1);
static SwitchOption interfaceGammaZZ
(interfaceGammaZ,
"Z",
"Only include the Z",
2);
}
Selector<MEBase::DiagramIndex>
MENeutralCurrentDIS::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if ( diags[i]->id() == -1 || diags[i]->id() == -3 ) sel.insert(meInfo()[0], i);
else if ( diags[i]->id() == -2 || diags[i]->id() == -4 ) sel.insert(meInfo()[1], i);
}
return sel;
}
double MENeutralCurrentDIS::helicityME(vector<SpinorWaveFunction> & f1,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder,
bool calc) const {
// scale
Energy2 mb2(scale());
// matrix element to be stored
ProductionMatrixElement menew (PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
ProductionMatrixElement gamma (PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
ProductionMatrixElement Zboson(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half);
// which intermediates to include
bool gam = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// declare the variables we need
VectorWaveFunction inter[2];
double me[3]={0.,0.,0.};
Complex diag1,diag2;
// sum over helicities to get the matrix element
unsigned int hel[4];
unsigned int lhel1,lhel2,qhel1,qhel2;
for(lhel1=0;lhel1<2;++lhel1) {
for(lhel2=0;lhel2<2;++lhel2) {
// intermediate for photon
if(gam) inter[0]=_theFFPVertex->evaluate(mb2,1,_gamma,f1[lhel1],a1[lhel2]);
// intermediate for Z
if(Z0) inter[1]=_theFFZVertex->evaluate(mb2,1,_z0 ,f1[lhel1],a1[lhel2]);
for(qhel1=0;qhel1<2;++qhel1) {
for(qhel2=0;qhel2<2;++qhel2) {
hel[0] = lhel1;
hel[1] = qhel1;
hel[2] = lhel2;
hel[3] = qhel2;
if(!lorder) swap(hel[0],hel[2]);
if(!qorder) swap(hel[1],hel[3]);
// first the photon exchange diagram
diag1 = gam ?
_theFFPVertex->evaluate(mb2,f2[qhel1],a2[qhel2],inter[0]) : 0.;
// then the Z exchange diagram
diag2 = Z0 ?
_theFFZVertex->evaluate(mb2,f2[qhel1],a2[qhel2],inter[1]) : 0.;
// add up squares of individual terms
me[1] += norm(diag1);
gamma (hel[0],hel[1],hel[2],hel[3]) = diag1;
me[2] += norm(diag2);
Zboson(hel[0],hel[1],hel[2],hel[3]) = diag2;
// the full thing including interference
diag1 += diag2;
me[0] += norm(diag1);
menew(hel[0],hel[1],hel[2],hel[3]) = diag1;
}
}
}
}
// spin and colour factor
double colspin = 0.25;
// results
for(int ix=0;ix<3;++ix) me[ix] *= colspin;
tcPolarizedBeamPDPtr beam[2] =
{dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[0]),
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[1])};
if( beam[0] || beam[1] ) {
RhoDMatrix rho[2] = {beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()),
beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())};
me[0] = menew .average(rho[0],rho[1]);
me[1] = gamma .average(rho[0],rho[1]);
me[2] = Zboson.average(rho[0],rho[1]);
}
DVector save;
save.push_back(me[1]);
save.push_back(me[2]);
meInfo(save);
if(calc) _me.reset(menew);
// analytic expression for testing
// double test = 8.*sqr(4.*Constants::pi*generator()->standardModel()->alphaEM(mb2))*
// sqr(double(mePartonData()[1]->iCharge())/3.)/sqr(tHat())
// *(sqr(sHat())+sqr(uHat())+4.*sqr(mePartonData()[0]->mass())*tHat())/4.;
// cerr << "testing me " << me[0]/test << "\n";
return me[0];
}
double MENeutralCurrentDIS::me2() const {
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
bool lorder,qorder;
SpinorWaveFunction l1,q1;
SpinorBarWaveFunction l2,q2;
// lepton wave functions
if(mePartonData()[0]->id()>0) {
lorder=true;
l1 = SpinorWaveFunction (meMomenta()[0],mePartonData()[0],incoming);
l2 = SpinorBarWaveFunction(meMomenta()[2],mePartonData()[2],outgoing);
}
else {
lorder=false;
l1 = SpinorWaveFunction (meMomenta()[2],mePartonData()[2],outgoing);
l2 = SpinorBarWaveFunction(meMomenta()[0],mePartonData()[0],incoming);
}
// quark wave functions
if(mePartonData()[1]->id()>0) {
qorder = true;
q1 = SpinorWaveFunction (meMomenta()[1],mePartonData()[1],incoming);
q2 = SpinorBarWaveFunction(meMomenta()[3],mePartonData()[3],outgoing);
}
else {
qorder = false;
q1 = SpinorWaveFunction (meMomenta()[3],mePartonData()[3],outgoing);
q2 = SpinorBarWaveFunction(meMomenta()[1],mePartonData()[1],incoming);
}
// wavefunctions for various helicities
for(unsigned int ix=0;ix<2;++ix) {
l1.reset(ix); f1.push_back(l1);
l2.reset(ix); a1.push_back(l2);
q1.reset(ix); f2.push_back(q1);
q2.reset(ix); a2.push_back(q2);
}
return helicityME(f1,f2,a1,a2,lorder,qorder,false);
}
void MENeutralCurrentDIS::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);
hard.push_back(sub->outgoing()[1]);
unsigned int order[4]={0,1,2,3};
bool lorder(true),qorder(true);
if(hard[0]->id()<0) {
swap(order[0],order[2]);
lorder = false;
}
if(hard[1]->id()<0) {
swap(order[1],order[3]);
qorder = false;
}
vector<SpinorWaveFunction> f1,f2;
vector<SpinorBarWaveFunction> a1,a2;
SpinorWaveFunction (f1,hard[order[0]],incoming,!lorder,true);
SpinorWaveFunction (f2,hard[order[1]],incoming,!qorder,true);
SpinorBarWaveFunction(a1,hard[order[2]],outgoing, lorder,true);
SpinorBarWaveFunction(a2,hard[order[3]],outgoing, qorder,true);
helicityME(f1,f2,a1,a2,lorder,qorder,false);
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix) {
tSpinPtr spin = hard[ix]->spinInfo();
if(ix<2) {
tcPolarizedBeamPDPtr beam =
dynamic_ptr_cast<tcPolarizedBeamPDPtr>(hard[ix]->dataPtr());
if(beam) spin->rhoMatrix() = beam->rhoMatrix();
}
spin->productionVertex(hardvertex);
}
}
-
-
+double MENeutralCurrentDIS::A(tcPDPtr lin, tcPDPtr,
+ tcPDPtr qin, tcPDPtr, Energy2 q2) const {
+ // photon only
+ if(_gammaZ==1) return 0.;
+ // everything else
+ double r = _gammaZ==0 || _gammaZ==2 ? double(q2/(q2+_mz2)) : 0.;
+ double eL,eQ,cvL,caL,cvQ,caQ;;
+ if(abs(lin->id())%2==0) {
+ eL = _gammaZ==0 ? generator()->standardModel()->enu() : 0.;
+ cvL = 0.25*generator()->standardModel()->vnu();
+ caL = 0.25*generator()->standardModel()->anu();
+ }
+ else {
+ eL = _gammaZ==0 ? generator()->standardModel()->ee() : 0.;
+ cvL = 0.25*generator()->standardModel()->ve();
+ caL = 0.25*generator()->standardModel()->ae();
+ }
+ if(abs(qin->id())%2==0) {
+ eQ = _gammaZ==0 ? generator()->standardModel()->eu() : 0.;
+ cvQ = 0.25*generator()->standardModel()->vu();
+ caQ = 0.25*generator()->standardModel()->au();
+ }
+ else {
+ eQ = _gammaZ==0 ? generator()->standardModel()->ed() : 0.;
+ cvQ = 0.25*generator()->standardModel()->vd();
+ caQ = 0.25*generator()->standardModel()->ad();
+ }
+ double output = 4.*r*caL*caQ*(eL*eQ+2.*r*cvL*cvQ/sqr(_sinW*_cosW))
+ /sqr(_sinW*_cosW)/(sqr(eL*eQ)+2.*eL*eQ*r/sqr(_cosW*_sinW)*cvL*cvQ
+ +sqr(r/sqr(_cosW*_sinW))*(sqr(cvL)+sqr(caL))*(sqr(cvQ)+sqr(caQ)));
+ if(qin->id()<0) output *= -1.;
+ if(lin->id()<0) output *= -1.;
+ return output;
+}
diff --git a/MatrixElement/DIS/MENeutralCurrentDIS.h b/MatrixElement/DIS/MENeutralCurrentDIS.h
--- a/MatrixElement/DIS/MENeutralCurrentDIS.h
+++ b/MatrixElement/DIS/MENeutralCurrentDIS.h
@@ -1,280 +1,307 @@
// -*- C++ -*-
#ifndef HERWIG_MENeutralCurrentDIS_H
#define HERWIG_MENeutralCurrentDIS_H
//
// This is the declaration of the MENeutralCurrentDIS class.
//
-#include "Herwig++/MatrixElement/HwMEBase.h"
+#include "DISBase.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MENeutralCurrentDIS class provides the matrix elements for
* neutral current DIS.
*
* For consistency both the incoming and outgoing quarks are assumed to be massless.
*
* @see \ref MENeutralCurrentDISInterfaces "The interfaces"
* defined for MENeutralCurrentDIS.
*/
-class MENeutralCurrentDIS: public HwMEBase {
+class MENeutralCurrentDIS: public DISBase {
public:
/**
* The default constructor.
*/
MENeutralCurrentDIS();
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* The matrix element for the kinematical configuration
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
* @return the matrix element scaled with sHat() to give a
* dimensionless number.
*/
virtual double me2() const;
/**
- * Return the scale associated with the last set phase space point.
- */
- virtual Energy2 scale() const;
-
- /**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Get diagram selector. With the information previously supplied with the
* setKinematics method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
* @param dv the diagrams to be weighted.
* @return a Selector relating the given diagrams to their weights.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
* @param diag the diagram chosen.
* @return the possible colour geometries weighted by their
* relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Construct the vertex of spin correlations.
*/
virtual void constructVertex(tSubProPtr);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Matrix element for \f$\ell q\to \gamma/Z \to \ell q\f$.
* @param f1 Fermion on lepton line
* @param a1 Anti-fermion on lepton line
* @param f2 Fermion on quark line
* @param a2 Anti-fermion on quark line
* @param lorder The order of particles on the lepton line
* @param qorder The order of particles on the quark line
* @param me Whether or not to calculate the matrix element for spin correlations
*/
double helicityME(vector<SpinorWaveFunction> & f1 ,
vector<SpinorWaveFunction> & f2,
vector<SpinorBarWaveFunction> & a1 ,
vector<SpinorBarWaveFunction> & a2,
bool lorder, bool qorder,
bool me) const;
+ /**
+ * Option for treatment of $\gamma/Z\f$ terms
+ */
+ inline unsigned int gammaZOption() const {return _gammaZ;}
+
+ /**
+ * Calculate the coefficient A for the correlations in the hard
+ * radiation
+ */
+ virtual double A(tcPDPtr lin, tcPDPtr lout, tcPDPtr qin, tcPDPtr qout,
+ Energy2 scale) const;
+
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MENeutralCurrentDIS> initMENeutralCurrentDIS;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MENeutralCurrentDIS & operator=(const MENeutralCurrentDIS &);
private:
/**
* Pointer to the vertices for the helicity calculations
*/
//@{
/**
* Pointer to the Z vertex
*/
AbstractFFVVertexPtr _theFFZVertex;
/**
* Pointer to the photon vertex
*/
AbstractFFVVertexPtr _theFFPVertex;
//@}
/**
* Pointers to the intermediate resonances
*/
//@{
/**
* Pointer to the Z ParticleData object
*/
tcPDPtr _z0;
/**
* Pointer to the photon ParticleData object
*/
tcPDPtr _gamma;
//@}
/**
* Switches to control the particles in the hard process
*/
//@{
/**
* Minimumflavour of the incoming quarks
*/
int _minflavour;
/**
* Maximum flavour of the incoming quarks
*/
int _maxflavour;
/**
* Whether to include both \f$Z^0\f$ and \f$\gamma\f$ or only one
*/
unsigned int _gammaZ;
//@}
/**
* Matrix element for spin correlations
*/
ProductionMatrixElement _me;
+ /**
+ * Electroweak parameters
+ */
+ //@{
+ /**
+ * \f$\sin\theta_W\f$
+ */
+ double _sinW;
+
+ /**
+ * \f$\cos\theta_W\f$
+ */
+ double _cosW;
+
+ /**
+ * The square of the Z mass
+ */
+ Energy2 _mz2;
+ //@}
+
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MENeutralCurrentDIS. */
template <>
struct BaseClassTrait<Herwig::MENeutralCurrentDIS,1> {
/** Typedef of the first base class of MENeutralCurrentDIS. */
- typedef Herwig::HwMEBase NthBase;
+ typedef Herwig::DISBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MENeutralCurrentDIS class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MENeutralCurrentDIS>
: public ClassTraitsBase<Herwig::MENeutralCurrentDIS> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MENeutralCurrentDIS"; }
/**
* The name of a file containing the dynamic library where the class
* MENeutralCurrentDIS is implemented. It may also include several, space-separated,
* libraries if the class MENeutralCurrentDIS depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwMEDIS.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MENeutralCurrentDIS_H */
diff --git a/MatrixElement/DIS/Makefile.am b/MatrixElement/DIS/Makefile.am
--- a/MatrixElement/DIS/Makefile.am
+++ b/MatrixElement/DIS/Makefile.am
@@ -1,5 +1,6 @@
pkglib_LTLIBRARIES = HwMEDIS.la
HwMEDIS_la_SOURCES = \
+DISBase.h DISBase.cc \
MENeutralCurrentDIS.cc MENeutralCurrentDIS.h \
MEChargedCurrentDIS.cc MEChargedCurrentDIS.h
HwMEDIS_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 3:0:0
diff --git a/Models/General/ModelGenerator.cc b/Models/General/ModelGenerator.cc
--- a/Models/General/ModelGenerator.cc
+++ b/Models/General/ModelGenerator.cc
@@ -1,468 +1,468 @@
// -*- C++ -*-
//
// ModelGenerator.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 ModelGenerator class.
//
#include "ModelGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "BSMWidthGenerator.h"
#include "Herwig++/PDT/GenericMassGenerator.h"
#include "Herwig++/Decay/DecayIntegrator.h"
#include "ThePEG/Repository/BaseRepository.h"
using namespace Herwig;
IBPtr ModelGenerator::clone() const {
return new_ptr(*this);
}
IBPtr ModelGenerator::fullclone() const {
return new_ptr(*this);
}
void ModelGenerator::persistentOutput(PersistentOStream & os) const {
os << hardProcessConstructors_ << _theDecayConstructor << particles_
<< offshell_ << Offsel_ << BRnorm_
<< Npoints_ << Iorder_ << BWshape_ << brMin_ << decayOutput_;
}
void ModelGenerator::persistentInput(PersistentIStream & is, int) {
is >> hardProcessConstructors_ >> _theDecayConstructor >> particles_
>> offshell_ >> Offsel_ >> BRnorm_
>> Npoints_ >> Iorder_ >> BWshape_ >> brMin_ >> decayOutput_;
}
ClassDescription<ModelGenerator> ModelGenerator::initModelGenerator;
// Definition of the static class description member.
void ModelGenerator::Init() {
static ClassDocumentation<ModelGenerator> documentation
("This class controls the the use of BSM physics.",
"BSM physics was produced using the algorithm of "
"\\cite{Gigg:2007cr,Gigg:2008yc}",
"\\bibitem{Gigg:2007cr} M.~Gigg and P.~Richardson, \n"
"Eur.\\ Phys.\\ J.\\ C {\\bf 51} (2007) 989.\n"
"%%CITATION = EPHJA,C51,989;%%\n"
" %\\cite{Gigg:2008yc}\n"
"\\bibitem{Gigg:2008yc}\n"
" M.~A.~Gigg and P.~Richardson,\n"
" %``Simulation of Finite Width Effects in Physics Beyond the Standard Model,''\n"
" arXiv:0805.3037 [hep-ph].\n"
" %%CITATION = ARXIV:0805.3037;%%\n"
);
static RefVector<ModelGenerator,HardProcessConstructor>
interfaceHardProcessConstructors
("HardProcessConstructors",
"The objects to construct hard processes",
&ModelGenerator::hardProcessConstructors_, -1,
false, false, true, false, false);
static Reference<ModelGenerator,Herwig::DecayConstructor>
interfaceDecayConstructor
("DecayConstructor",
"Pointer to DecayConstructor helper class",
&ModelGenerator::_theDecayConstructor, false, false, true, false);
static RefVector<ModelGenerator,ThePEG::ParticleData> interfaceModelParticles
("DecayParticles",
"ParticleData pointers to the particles requiring spin correlation "
"decayers. If decay modes do not exist they will also be created.",
&ModelGenerator::particles_, -1, false, false, true, false);
static RefVector<ModelGenerator,ParticleData> interfaceOffshell
("Offshell",
"The particles to treat as off-shell",
&ModelGenerator::offshell_, -1, false, false, true, false);
static Switch<ModelGenerator,int> interfaceWhichOffshell
("WhichOffshell",
"A switch to determine which particles to create mass and width "
"generators for.",
&ModelGenerator::Offsel_, 0, false, false);
static SwitchOption interfaceWhichOffshellSelected
(interfaceWhichOffshell,
"Selected",
"Only create mass and width generators for the particles specified",
0);
static SwitchOption interfaceWhichOffshellAll
(interfaceWhichOffshell,
"All",
"Treat all particles specified in the DecayParticles "
"list as off-shell",
1);
static Switch<ModelGenerator,bool> interfaceBRNormalize
("BRNormalize",
"Whether to normalize the partial widths to BR*total width for an "
"on-shell particle",
&ModelGenerator::BRnorm_, true, false, false);
static SwitchOption interfaceBRNormalizeNormalize
(interfaceBRNormalize,
"Yes",
"Normalize the partial widths",
true);
static SwitchOption interfaceBRNormalizeNoNormalize
(interfaceBRNormalize,
"No",
"Do not normalize the partial widths",
false);
static Parameter<ModelGenerator,int> interfacePoints
("InterpolationPoints",
"Number of points to use for interpolation tables when needed",
&ModelGenerator::Npoints_, 50, 5, 1000,
false, false, true);
static Parameter<ModelGenerator,unsigned int>
interfaceInterpolationOrder
("InterpolationOrder", "The interpolation order for the tables",
&ModelGenerator::Iorder_, 1, 1, 5,
false, false, Interface::limited);
static Switch<ModelGenerator,int> interfaceBreitWignerShape
("BreitWignerShape",
"Controls the shape of the mass distribution generated",
&ModelGenerator::BWshape_, 0, false, false);
static SwitchOption interfaceBreitWignerShapeDefault
(interfaceBreitWignerShape,
"Default",
"Running width with q in numerator and denominator width factor",
0);
static SwitchOption interfaceBreitWignerShapeFixedWidth
(interfaceBreitWignerShape,
"FixedWidth",
"Use a fixed width",
1);
static SwitchOption interfaceBreitWignerShapeNoq
(interfaceBreitWignerShape,
"Noq",
"Use M rather than q in the numerator and denominator width factor",
2);
static SwitchOption interfaceBreitWignerShapeNoNumerator
(interfaceBreitWignerShape,
"NoNumerator",
"Neglect the numerator factors",
3);
static Parameter<ModelGenerator,double> interfaceMinimumBR
("MinimumBR",
"The minimum branching fraction to include",
&ModelGenerator::brMin_, 1e-6, 0.0, 1.0,
false, false, Interface::limited);
static Switch<ModelGenerator,unsigned int> interfaceDecayOutput
("DecayOutput",
"Option to control the output of the decay mode information",
&ModelGenerator::decayOutput_, 1, false, false);
static SwitchOption interfaceDecayOutputNone
(interfaceDecayOutput,
"None",
"No output",
0);
static SwitchOption interfaceDecayOutputPlain
(interfaceDecayOutput,
"Plain",
"Default plain text output",
1);
static SwitchOption interfaceDecayOutputSLHA
(interfaceDecayOutput,
"SLHA",
"Output in the Susy Les Houches Accord format",
2);
}
namespace {
/// Helper function for sorting by mass
inline bool massIsLess(tcPDPtr a, tcPDPtr b) {
return a->mass() < b->mass();
}
}
void ModelGenerator::doinit() {
useMe();
Interfaced::doinit();
// make sure the model is initialized
Ptr<Herwig::StandardModel>::pointer model
= dynamic_ptr_cast<Ptr<Herwig::StandardModel>::pointer>(generator()->standardModel());
model->init();
// and the vertices
for(size_t iv = 0; iv < model->numberOfVertices(); ++iv)
model->vertex(iv)->init();
// sort DecayParticles list by mass
sort(particles_.begin(),particles_.end(),
massIsLess);
//create mass and width generators for the requested particles
PDVector::iterator pit, pend;
if( Offsel_ == 0 ) {
pit = offshell_.begin();
pend = offshell_.end();
}
else {
pit = particles_.begin();
pend = particles_.end();
}
for(; pit != pend; ++pit)
createWidthGenerator(*pit);
//create decayers and decaymodes (if necessary)
if( _theDecayConstructor ) {
_theDecayConstructor->init();
_theDecayConstructor->createDecayers(particles_,brMin_);
}
// write out decays with spin correlations
ostream & os = CurrentGenerator::current().misc();
ofstream ofs;
- if ( decayOutput_ >1 ) {
+ if ( decayOutput_ > 1 ) {
string filename
= CurrentGenerator::current().filename() + "-BR.spc";
ofs.open(filename.c_str());
}
if(decayOutput_!=0) {
if(decayOutput_==1) {
os << "# The decay modes listed below will have spin\n"
<< "# correlations included when they are generated.\n#\n#";
}
else {
ofs << "# Herwig++ decay tables in SUSY Les Houches accord format\n";
ofs << "Block DCINFO # Program information\n";
ofs << "1 Herwig++ # Decay Calculator\n";
ofs << "2 " << generator()->strategy()->versionstring()
<< " # Version number\n";
}
}
pit = particles_.begin();
pend = particles_.end();
for( ; pit != pend; ++pit) {
tPDPtr parent = *pit;
// Check decays for ones where quarks cannot be put on constituent
// mass-shell
checkDecays(parent);
parent->reset();
parent->update();
if( parent->CC() ) parent->CC()->synchronize();
if( parent->decaySelector().empty() ) {
parent->stable(true);
parent->width(ZERO);
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
else {
if ( decayOutput_ == 2 )
writeDecayModes(ofs, parent);
else
writeDecayModes(os, parent);
}
if( parent->massGenerator() ) {
parent->widthCut(5.*parent->width());
parent->massGenerator()->reset();
if(decayOutput_==1)
os << "# " <<parent->PDGName() << " will be considered off-shell.\n#\n";
}
if( parent->widthGenerator() ) parent->widthGenerator()->reset();
}
//Now construct hard processes given that we know which
//objects have running widths
for(unsigned int ix=0;ix<hardProcessConstructors_.size();++ix) {
hardProcessConstructors_[ix]->init();
hardProcessConstructors_[ix]->constructDiagrams();
}
}
void ModelGenerator::checkDecays(PDPtr parent) {
if( parent->stable() ) {
if(parent->coloured())
cerr << "Warning: No decays for coloured particle " << parent->PDGName() << "\n\n"
<< "have been calcluated in BSM model.\n"
<< "This may cause problems in the hadronization phase.\n"
<< "You may have forgotten to switch on the decay mode calculation using\n"
<< " set TwoBodyDC:CreateDecayModes Yes\n"
<< " set ThreeBodyDC:CreateDecayModes Yes\n"
<< " set WeakDecayConstructor:CreateDecayModes Yes\n"
<< "or the decays of this particle are missing from your\n"
<< "input spectrum and decay file in the SLHA format.\n\n";
return;
}
DecaySet::iterator dit = parent->decayModes().begin();
DecaySet::iterator dend = parent->decayModes().end();
Energy oldwidth(parent->width()), newwidth(ZERO);
bool rescalebrat(false);
double brsum(0.);
for(; dit != dend; ++dit ) {
if( !(**dit).on() ) continue;
Energy release((**dit).parent()->mass());
tPDVector::const_iterator pit = (**dit).orderedProducts().begin();
tPDVector::const_iterator pend =(**dit).orderedProducts().end();
for( ; pit != pend; ++pit ) {
release -= (**pit).constituentMass();
}
if( (**dit).brat() < brMin_ || release < ZERO ) {
if( release < ZERO )
cerr << "Warning: The shower cannot be generated using this decay "
<< (**dit).tag() << " because it is too close to threshold. It "
<< "will be switched off and the branching fractions of the "
<< "remaining modes rescaled.\n";
rescalebrat = true;
generator()->preinitInterface(*dit, "OnOff", "set", "Off");
generator()->preinitInterface(*dit, "BranchingRatio",
"set", "0.0");
DecayIntegratorPtr decayer = dynamic_ptr_cast<DecayIntegratorPtr>((**dit).decayer());
if(decayer) {
generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0");
}
}
else {
brsum += (**dit).brat();
newwidth += (**dit).brat()*oldwidth;
}
}
if( rescalebrat || (abs(brsum - 1.) > 1e-12) ) {
dit = parent->decayModes().begin();
dend = parent->decayModes().end();
double factor = oldwidth/newwidth;
brsum = 0.;
for( ; dit != dend; ++dit ) {
if( !(**dit).on() ) continue;
double newbrat = ((**dit).brat())*factor;
brsum += newbrat;
ostringstream brf;
brf << setprecision(13) << newbrat;
generator()->preinitInterface(*dit, "BranchingRatio",
"set", brf.str());
}
parent->width(newwidth);
if( newwidth > ZERO ) parent->cTau(hbarc/newwidth);
}
}
namespace {
struct DecayModeOrdering {
bool operator()(tcDMPtr m1, tcDMPtr m2) {
if(m1->brat()!=m2->brat()) {
return m1->brat()>m2->brat();
}
else {
if(m1->products().size()==m2->products().size()) {
ParticleMSet::const_iterator it1=m1->products().begin();
ParticleMSet::const_iterator it2=m2->products().begin();
do {
if((**it1).id()!=(**it2).id()) {
return (**it1).id()>(**it2).id();
}
++it1;
++it2;
}
while(it1!=m1->products().end()&&
it2!=m2->products().end());
assert(false);
}
else
return m1->products().size()<m2->products().size();
}
return false;
}
};
}
void ModelGenerator::writeDecayModes(ostream & os, tcPDPtr parent) const {
if(decayOutput_==0) return;
set<tcDMPtr,DecayModeOrdering> modes;
for(Selector<tDMPtr>::const_iterator dit = parent->decaySelector().begin();
dit != parent->decaySelector().end(); ++dit) {
modes.insert((*dit).second);
}
if(decayOutput_==1) {
os << " Parent: " << parent->PDGName() << " Mass (GeV): "
<< parent->mass()/GeV << " Total Width (GeV): "
<< parent->width()/GeV << endl;
os << std::left << std::setw(40) << '#'
<< std::left << std::setw(20) << "Partial Width/GeV"
<< "BR\n";
for(set<tcDMPtr,DecayModeOrdering>::iterator dit=modes.begin();
dit!=modes.end();++dit)
os << std::left << std::setw(40) << (**dit).tag()
<< std::left << std::setw(20) << (**dit).brat()*parent->width()/GeV
<< (**dit).brat() << '\n';
os << "#\n#";
}
else if(decayOutput_==2) {
os << "# \t PDG \t Width\n";
os << "DECAY\t" << parent->id() << "\t" << parent->width()/GeV << "\t # " << parent->PDGName() << "\n";
for(set<tcDMPtr,DecayModeOrdering>::iterator dit=modes.begin();
dit!=modes.end();++dit) {
os << "\t" << std::left << std::setw(10)
<< (**dit).brat() << "\t" << (**dit).orderedProducts().size()
<< "\t";
for(unsigned int ix=0;ix<(**dit).orderedProducts().size();++ix)
os << std::right << std::setw(10)
<< (**dit).orderedProducts()[ix]->id() ;
for(unsigned int ix=(**dit).orderedProducts().size();ix<4;++ix)
os << "\t";
os << "# " << (**dit).tag() << "\n";
}
}
}
void ModelGenerator::createWidthGenerator(tPDPtr p) {
string wn = p->fullName() + string("-WGen");
string mn = p->fullName() + string("-MGen");
GenericMassGeneratorPtr mgen = dynamic_ptr_cast<GenericMassGeneratorPtr>
(generator()->preinitCreate("Herwig::GenericMassGenerator", mn));
BSMWidthGeneratorPtr wgen = dynamic_ptr_cast<BSMWidthGeneratorPtr>
(generator()->preinitCreate("Herwig::BSMWidthGenerator", wn));
//set the particle interface
generator()->preinitInterface(mgen, "Particle", "set", p->fullName());
generator()->preinitInterface(wgen, "Particle", "set", p->fullName());
//set the generator interfaces in the ParticleData object
generator()->preinitInterface(p, "Mass_generator","set", mn);
generator()->preinitInterface(p, "Width_generator","set", wn);
//allow the branching fraction of this particle type to vary
p->variableRatio(true);
if( p->CC() ) p->CC()->variableRatio(true);
//initialize the generators
generator()->preinitInterface(mgen, "Initialize", "set", "Yes");
generator()->preinitInterface(wgen, "Initialize", "set", "Yes");
string norm = BRnorm_ ? "Yes" : "No";
generator()->preinitInterface(wgen, "BRNormalize", "set", norm);
ostringstream os;
os << Npoints_;
generator()->preinitInterface(wgen, "InterpolationPoints", "set",
os.str());
os.str("");
os << Iorder_;
generator()->preinitInterface(wgen, "InterpolationOrder", "set",
os.str());
os.str("");
os << BWshape_;
generator()->preinitInterface(mgen, "BreitWignerShape", "set",
os.str());
}
diff --git a/Models/General/TwoBodyDecay.h b/Models/General/TwoBodyDecay.h
--- a/Models/General/TwoBodyDecay.h
+++ b/Models/General/TwoBodyDecay.h
@@ -1,66 +1,80 @@
// -*- C++ -*-
//
// TwoBodyDecay.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_TwoBodyDecay_H
#define HERWIG_TwoBodyDecay_H
//
// This is the declaration of the TwoBodyDecay struct.
//
namespace Herwig {
using namespace ThePEG;
using Helicity::tVertexBasePtr;
/**
* Struct for the prototype of a two-body decay mode
*/
struct TwoBodyDecay {
public:
/**
* Constructor
* @param pa Decaying particle
* @param pb First decay product
* @param pc Second decay product
*/
TwoBodyDecay(tPDPtr pa, tPDPtr pb, tPDPtr pc,
tVertexBasePtr vertex) : parent_(pa), vertex_(vertex) {
ParticleOrdering order;
if( order(pb, pc) ) {
children_.first = pb;
children_.second = pc;
}
else {
children_.first = pc;
children_.second = pb;
}
}
/**
* The parent
*/
tPDPtr parent_;
/**
* The children
*/
tPDPair children_;
/**
* Vertex
*/
tVertexBasePtr vertex_;
private:
TwoBodyDecay();
};
+
+/**
+ * Test whether one TwoBodyDecay is less than another
+ */
+inline bool operator<(const TwoBodyDecay & x, const TwoBodyDecay & y) {
+ if(x.parent_->id()!=y.parent_->id())
+ return x.parent_->id()<y.parent_->id();
+ if(x.children_.first->id()!=y.children_.first->id())
+ return x.children_.first->id() < y.children_.first->id();
+ if(x.children_.second->id()!=y.children_.second->id())
+ return x.children_.second->id() < y.children_.second->id();
+ return x.vertex_<y.vertex_;
+}
+
}
#endif /* HERWIG_TwoBodyDecay_H */
diff --git a/Models/General/TwoBodyDecayConstructor.cc b/Models/General/TwoBodyDecayConstructor.cc
--- a/Models/General/TwoBodyDecayConstructor.cc
+++ b/Models/General/TwoBodyDecayConstructor.cc
@@ -1,284 +1,284 @@
// -*- C++ -*-
//
// TwoBodyDecayConstructor.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 TwoBodyDecayConstructor class.
//
#include "TwoBodyDecayConstructor.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "Herwig++/Decay/General/GeneralTwoBodyDecayer.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "DecayConstructor.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractSSTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractRFSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractRFVVertex.fh"
using namespace Herwig;
using ThePEG::Helicity::VertexBasePtr;
IBPtr TwoBodyDecayConstructor::clone() const {
return new_ptr(*this);
}
IBPtr TwoBodyDecayConstructor::fullclone() const {
return new_ptr(*this);
}
NoPIOClassDescription<TwoBodyDecayConstructor>
TwoBodyDecayConstructor::initTwoBodyDecayConstructor;
// Definition of the static class description member.
void TwoBodyDecayConstructor::Init() {
static ClassDocumentation<TwoBodyDecayConstructor> documentation
("The TwoBodyDecayConstructor implements to creation of 2 body decaymodes "
"and decayers that do not already exist for the given set of vertices.");
}
void TwoBodyDecayConstructor::DecayList(const set<PDPtr> & particles) {
if( particles.empty() ) return;
tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
unsigned int nv(model->numberOfVertices());
for(set<PDPtr>::const_iterator ip=particles.begin();
ip!=particles.end();++ip) {
tPDPtr parent = *ip;
for(unsigned int iv = 0; iv < nv; ++iv) {
if(excluded(model->vertex(iv)) ||
model->vertex(iv)->getNpoint()>3) continue;
for(unsigned int il = 0; il < 3; ++il) {
- vector<TwoBodyDecay> decays =
+ set<TwoBodyDecay> decays =
createModes(parent, model->vertex(iv), il);
if( !decays.empty() ) createDecayMode(decays);
}
}
}
}
-vector<TwoBodyDecay> TwoBodyDecayConstructor::
+set<TwoBodyDecay> TwoBodyDecayConstructor::
createModes(tPDPtr inpart, VertexBasePtr vertex,
unsigned int list) {
int id = inpart->id();
if( id < 0 || !vertex->isIncoming(inpart) || vertex->getNpoint() != 3 )
- return vector<TwoBodyDecay>();
+ return set<TwoBodyDecay>();
Energy m1(inpart->mass());
tPDVector decaylist = vertex->search(list, inpart);
- vector<TwoBodyDecay> decays;
+ set<TwoBodyDecay> decays;
tPDVector::size_type nd = decaylist.size();
for( tPDVector::size_type i = 0; i < nd; i += 3 ) {
tPDPtr pa(decaylist[i]), pb(decaylist[i + 1]), pc(decaylist[i + 2]);
if( pb->id() == id ) swap(pa, pb);
if( pc->id() == id ) swap(pa, pc);
//allowed on-shell decay?
if( m1 <= pb->mass() + pc->mass() ) continue;
//vertices are defined with all particles incoming
if( pb->CC() ) pb = pb->CC();
if( pc->CC() ) pc = pc->CC();
- decays.push_back( TwoBodyDecay(inpart,pb, pc, vertex) );
+ decays.insert( TwoBodyDecay(inpart,pb, pc, vertex) );
}
return decays;
}
-GeneralTwoBodyDecayerPtr TwoBodyDecayConstructor::createDecayer(TwoBodyDecay & decay) {
+GeneralTwoBodyDecayerPtr TwoBodyDecayConstructor::createDecayer(TwoBodyDecay decay) {
string name;
using namespace Helicity::VertexType;
PDT::Spin in = decay.parent_->iSpin();
// PDT::Spin out1 = decay.children_.first ->iSpin();
PDT::Spin out2 = decay.children_.second->iSpin();
switch(decay.vertex_->getName()) {
case FFV :
if(in == PDT::Spin1Half) {
name = "FFVDecayer";
if(out2==PDT::Spin1Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "VFFDecayer";
}
break;
case FFS :
if(in == PDT::Spin1Half) {
name = "FFSDecayer";
if(out2==PDT::Spin1Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "SFFDecayer";
}
break;
case VVS :
if(in == PDT::Spin1) {
name = "VVSDecayer";
if(out2==PDT::Spin1)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "SVVDecayer";
}
break;
case VSS :
if(in == PDT::Spin1) {
name = "VSSDecayer";
}
else {
name = "SSVDecayer";
if(out2==PDT::Spin0)
swap(decay.children_.first,decay.children_.second);
}
break;
case VVT :
name = in==PDT::Spin2 ? "TVVDecayer" : "Unknown";
break;
case FFT :
name = in==PDT::Spin2 ? "TFFDecayer" : "Unknown";
break;
case SST :
name = in==PDT::Spin2 ? "TSSDecayer" : "Unknown";
break;
case SSS :
name = "SSSDecayer";
break;
case VVV :
name = "VVVDecayer";
break;
case RFS :
if(in==PDT::Spin1Half) {
name = "FRSDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else if(in==PDT::Spin0) {
name = "SRFDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "Unknown";
}
break;
case RFV :
if(in==PDT::Spin1Half) {
name = "FRVDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else
name = "Unknown";
break;
default : Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName();
}
if(name=="Unknown")
Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName();
ostringstream fullname;
fullname << "/Herwig/Decays/" << name << "_" << decay.parent_->PDGName()
<< "_" << decay.children_.first ->PDGName()
<< "_" << decay.children_.second->PDGName();
string classname = "Herwig::" + name;
GeneralTwoBodyDecayerPtr decayer;
decayer = dynamic_ptr_cast<GeneralTwoBodyDecayerPtr>
(generator()->preinitCreate(classname,fullname.str()));
if(!decayer)
Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName();
decayer->setDecayInfo(decay.parent_,decay.children_,decay.vertex_);
decayer->init();
setDecayerInterfaces(fullname.str());
return decayer;
}
void TwoBodyDecayConstructor::
-createDecayMode(vector<TwoBodyDecay> & decays) {
- tPDPtr inpart = decays[0].parent_;
+createDecayMode(set<TwoBodyDecay> & decays) {
+ tPDPtr inpart = decays.begin()->parent_;
tEGPtr eg = generator();
- vector<TwoBodyDecay>::iterator dend = decays.end();
- for( vector<TwoBodyDecay>::iterator dit = decays.begin();
+ set<TwoBodyDecay>::iterator dend = decays.end();
+ for( set<TwoBodyDecay>::iterator dit = decays.begin();
dit != dend; ++dit ) {
tPDPtr pb((*dit).children_.first), pc((*dit).children_.second);
string tag = inpart->name() + "->" + pb->name() + "," +
pc->name() + ";";
// Does it exist already ?
tDMPtr dm = eg->findDecayMode(tag);
// Check if tag is one that should be disabled
if( decayConstructor()->disableDecayMode(tag) ) {
// If mode alread exists, ie has been read from file,
// disable it
if( dm ) {
eg->preinitInterface(dm, "BranchingRatio", "set", "0.0");
eg->preinitInterface(dm, "OnOff", "set", "Off");
}
continue;
}
// now create DecayMode objects that do not already exist
if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) {
tDMPtr ndm = eg->preinitCreateDecayMode(tag);
if(ndm) {
inpart->stable(false);
GeneralTwoBodyDecayerPtr decayer=createDecayer(*dit);
if(!decayer) continue;
eg->preinitInterface(ndm, "Decayer", "set",
decayer->fullName());
eg->preinitInterface(ndm, "OnOff", "set", "On");
Energy width =
decayer->partialWidth(make_pair(inpart,inpart->mass()),
make_pair(pb,pb->mass()) ,
make_pair(pc,pc->mass()));
setBranchingRatio(ndm, width);
if(ndm->brat()<decayConstructor()->minimumBR()) {
generator()->preinitInterface(decayer->fullName(),
"Initialize", "set","0");
}
}
else
throw NBodyDecayConstructorError()
<< "TwoBodyDecayConstructor::createDecayMode - Needed to create "
<< "new decaymode but one could not be created for the tag "
<< tag << Exception::warning;
}
else if( dm ) {
if(dm->brat()<decayConstructor()->minimumBR()) {
continue;
}
if((dm->decayer()->fullName()).find("Mambo") != string::npos) {
inpart->stable(false);
GeneralTwoBodyDecayerPtr decayer=createDecayer(*dit);
if(!decayer) continue;
eg->preinitInterface(dm, "Decayer", "set",
decayer->fullName());
}
}
}
// update CC mode if it exists
if( inpart->CC() ) inpart->CC()->synchronize();
}
diff --git a/Models/General/TwoBodyDecayConstructor.h b/Models/General/TwoBodyDecayConstructor.h
--- a/Models/General/TwoBodyDecayConstructor.h
+++ b/Models/General/TwoBodyDecayConstructor.h
@@ -1,158 +1,158 @@
// -*- C++ -*-
//
// TwoBodyDecayConstructor.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_TwoBodyDecayConstructor_H
#define HERWIG_TwoBodyDecayConstructor_H
//
// This is the declaration of the TwoBodyDecayConstructor class.
//
#include "NBodyDecayConstructorBase.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "Herwig++/Decay/General/GeneralTwoBodyDecayer.fh"
#include "TwoBodyDecay.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
using Helicity::tVertexBasePtr;
/**
* The TwoBodyDecayConstructor class inherits from the dummy base class
* NBodyDecayConstructorBase and implements the necessary functions in
* order to create the 2 body decay modes for a given set of vertices
* stored in a Model class.
*
* @see \ref TwoBodyDecayConstructorInterfaces "The interfaces"
* defined for TwoBodyDecayConstructor.
* @see NBodyDecayConstructor
**/
class TwoBodyDecayConstructor: public NBodyDecayConstructorBase {
public:
/**
* The default constructor.
*/
TwoBodyDecayConstructor() {}
/**
* Function used to determine allowed decaymodes
*@param part vector of ParticleData pointers containing particles in model
*/
virtual void DecayList(const set<PDPtr> & part);
/**
* Number of outgoing lines. Required for correct ordering.
*/
virtual unsigned int numBodies() const { return 2; }
public:
/**
* 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:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static NoPIOClassDescription<TwoBodyDecayConstructor> initTwoBodyDecayConstructor;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
TwoBodyDecayConstructor & operator=(const TwoBodyDecayConstructor &);
private:
/** @name Functions to create decayers and decaymodes. */
//@{
/**
* Function to create decays
* @param inpart Incoming particle
* @param vert The vertex to create decays for
* @param ilist Which list to search
* @param iv Row number in _theExistingDecayers member
* @return A vector a decay modes
*/
- vector<TwoBodyDecay> createModes(tPDPtr inpart, VertexBasePtr vert,
- unsigned int ilist);
+ set<TwoBodyDecay> createModes(tPDPtr inpart, VertexBasePtr vert,
+ unsigned int ilist);
/**
* Function to create decayer for specific vertex
* @param decay decay mode for this decay
* member variable
*/
- GeneralTwoBodyDecayerPtr createDecayer(TwoBodyDecay & decay);
+ GeneralTwoBodyDecayerPtr createDecayer(TwoBodyDecay decay);
/**
* Create decay mode(s) from given part and decay modes
* @param decays The vector of decay modes
* @param decayer The decayer responsible for this decay
*/
- void createDecayMode(vector<TwoBodyDecay> & decays);
+ void createDecayMode(set<TwoBodyDecay> & decays);
//@}
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of TwoBodyDecayConstructor. */
template <>
struct BaseClassTrait<Herwig::TwoBodyDecayConstructor,1> {
/** Typedef of the first base class of TwoBodyDecayConstructor. */
typedef Herwig::NBodyDecayConstructorBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the TwoBodyDecayConstructor class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::TwoBodyDecayConstructor>
: public ClassTraitsBase<Herwig::TwoBodyDecayConstructor> {
/** Return a platform-independent class name */
static string className() { return "Herwig::TwoBodyDecayConstructor"; }
};
/** @endcond */
}
#endif /* HERWIG_TwoBodyDecayConstructor_H */
diff --git a/Models/Susy/SSGOGOHVertex.cc b/Models/Susy/SSGOGOHVertex.cc
--- a/Models/Susy/SSGOGOHVertex.cc
+++ b/Models/Susy/SSGOGOHVertex.cc
@@ -1,234 +1,235 @@
// -*- C++ -*-
//
// SSGOGOHVertex.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 SSGOGOHVertex class.
//
#include "SSGOGOHVertex.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include <cassert>
using namespace ThePEG::Helicity;
using namespace Herwig;
SSGOGOHVertex::SSGOGOHVertex() : theMw(), theSij(2, vector<Complex>(2,0.0)),
theQij(2, vector<Complex>(2,0.0)),
theQijLp(4, vector<Complex>(2,0.0)),
theQijRp(4, vector<Complex>(2,0.0)),
theSijdp(4, vector<Complex>(4,0.0)),
theQijdp(4, vector<Complex>(4,0.0)),
theSa(0.0), theSb(0.0),
theCa(0.0), theCb(0.0), theCoupLast(0.0),
theLLast(0.0), theRLast(0.0), theHLast(0),
theID1Last(0), theID2Last(0), theq2last() {
orderInGem(1);
orderInGs(0);
}
void SSGOGOHVertex::doinit() {
long neu[4] = {1000022, 1000023, 1000025, 1000035};
long chg[2] = {1000024, 1000037};
long higgs[3] = {25, 35, 36};
for(unsigned int i = 0; i < 4; ++i) {
//neutralinos
for(unsigned int j = 0; j < 4; ++j) {
for(unsigned int k = 0; k < 4; ++k) {
if( i < 3 ) {
- addToList(neu[j], neu[k], higgs[i]);
+ if(k<=j)
+ addToList(neu[j], neu[k], higgs[i]);
//charginos
if( j < 2 && k < 2 ) {
addToList(-chg[j], chg[k], higgs[i]);
}
}
else {
if( k == 2 ) break;
addToList(-chg[k], neu[j], 37);
addToList( chg[k], neu[j],-37);
}
}
}
}
FFSVertex::doinit();
tMSSMPtr theMSSM = dynamic_ptr_cast<tMSSMPtr>(generator()->standardModel());
if( !theMSSM )
throw InitException()
<< "SSGOGOHVertex::doinit() - The pointer to the MSSM object is null!"
<< Exception::abortnow;
theMw = getParticleData(ParticleID::Wplus)->mass();
double theSw = sqrt(sin2ThetaW());
double tw = theSw/sqrt(1. - theSw*theSw);
double tanb = theMSSM->tanBeta();
theSb = tanb/sqrt(1. + sqr(tanb));
theCb = sqrt( 1. - sqr(theSb) );
theSa = sin(theMSSM->higgsMixingAngle());
theCa = sqrt(1. - sqr(theSa));
MixingMatrix nmix = *theMSSM->neutralinoMix();
MixingMatrix umix = *theMSSM->charginoUMix();
MixingMatrix vmix = *theMSSM->charginoVMix();
for(unsigned int i = 0; i < 4; ++i) {
for(unsigned int j = 0; j < 4; ++j) {
if( i < 2 && j < 2 ) {
theQij[i][j] = vmix(i,0)*umix(j,1)/sqrt(2);
theSij[i][j] = vmix(i,1)*umix(j,0)/sqrt(2);
}
if( j < 2 ) {
theQijLp[i][j] = conj(nmix(i, 3)*vmix(j,0)
+ (nmix(i,1) + nmix(i,0)*tw)*vmix(j,1)/sqrt(2));
theQijRp[i][j] = nmix(i, 2)*umix(j,0)
- (nmix(i,1) + nmix(i,0)*tw)*umix(j,1)/sqrt(2);
}
theQijdp[i][j] = 0.5*( nmix(i,2)*( nmix(j,1) - tw*nmix(j,0) )
+ nmix(j,2)*( nmix(i,1) - tw*nmix(i,0) ) );
theSijdp[i][j] = 0.5*( nmix(i,3)*( nmix(j,1) - tw*nmix(j,0) )
+ nmix(j,3)*( nmix(i,1) - tw*nmix(i,0) ) );
}
}
}
void SSGOGOHVertex::persistentOutput(PersistentOStream & os) const {
os << theSij << theQij << theQijLp << theQijRp << theSijdp
<< theQijdp << ounit(theMw,GeV) << theSa << theSb << theCa
<< theCb;
}
void SSGOGOHVertex::persistentInput(PersistentIStream & is, int) {
is >> theSij >> theQij >> theQijLp >> theQijRp >> theSijdp
>> theQijdp >> iunit(theMw,GeV) >> theSa >> theSb >> theCa
>> theCb;
}
ClassDescription<SSGOGOHVertex> SSGOGOHVertex::initSSGOGOHVertex;
// Definition of the static class description member.
void SSGOGOHVertex::Init() {
static ClassDocumentation<SSGOGOHVertex> documentation
("The coupling of the higgs bosons to SM fermions in the MSSM");
}
/// \todo fixme
void SSGOGOHVertex::setCoupling(Energy2 q2, tcPDPtr particle1,
tcPDPtr particle2,tcPDPtr particle3) {
long f1ID(particle1->id()), f2ID(particle2->id()), higgsID(particle3->id());
assert(higgsID == ParticleID::h0 || higgsID == ParticleID::H0 ||
higgsID == ParticleID::A0 || abs(higgsID) == ParticleID::Hplus);
if( f1ID < 0 ) swap(f1ID, f2ID);
if( q2 != theq2last || theCoupLast == 0. ) {
theCoupLast = weakCoupling(q2);
theq2last = q2;
}
if( higgsID == theHLast && f1ID == theID1Last && f2ID == theID2Last) {
norm(theCoupLast);
left(theLLast);
right(theRLast);
return;
}
theHLast = higgsID;
theID1Last = f1ID;
theID2Last = f2ID;
if( higgsID == ParticleID::h0 ) {
//charginos
if( abs(f2ID) == ParticleID::SUSY_chi_1plus ||
abs(f2ID) == ParticleID::SUSY_chi_2plus ) {
unsigned int ei = (abs(f1ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
unsigned int ej = (abs(f2ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
theLLast = conj(theQij[ej][ei])*theSa - conj(theSij[ej][ei])*theCa;
theRLast = theQij[ei][ej]*theSa - theSij[ei][ej]*theCa;
}
//neutralinos
else {
unsigned int ei(f1ID - ParticleID::SUSY_chi_10),
ej(f2ID - ParticleID::SUSY_chi_10);
if( ei > 1 )
ei = ( ei == 13 ) ? 3 : 2;
if( ej > 1 )
ej = ( ej == 13 ) ? 3 : 2;
theLLast = conj(theQijdp[ej][ei])*theSa + conj(theSijdp[ej][ei])*theCa;
theRLast = theQijdp[ei][ej]*theSa + theSijdp[ei][ej]*theCa ;
}
}
else if( higgsID == ParticleID::H0 ) {
//charginos
if( abs(f2ID) == ParticleID::SUSY_chi_1plus ||
abs(f2ID) == ParticleID::SUSY_chi_2plus ) {
unsigned int ei = (abs(f1ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
unsigned int ej = (abs(f2ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
theLLast = -conj(theQij[ej][ei])*theCa - conj(theSij[ej][ei])*theSa;
theRLast = -theQij[ei][ej]*theCa - theSij[ei][ej]*theSa;
}
//neutralinos
else {
unsigned int ei(f1ID - ParticleID::SUSY_chi_10),
ej(f2ID - ParticleID::SUSY_chi_10);
if( ei > 1 )
ei = ( ei == 13 ) ? 3 : 2;
if( ej > 1 )
ej = ( ej == 13 ) ? 3 : 2;
theLLast = -conj(theQijdp[ej][ei])*theCa + conj(theSijdp[ej][ei])*theSa;
theRLast = -theQijdp[ei][ej]*theCa + theSijdp[ei][ej]*theSa;
}
}
else if( higgsID == ParticleID::A0 ) {
if( abs(f2ID) == ParticleID::SUSY_chi_1plus ||
abs(f2ID) == ParticleID::SUSY_chi_2plus ) {
unsigned int ei = (abs(f1ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
unsigned int ej = (abs(f2ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
theLLast = Complex(0.,1.)*( conj(theQij[ej][ei])*theSb
+ conj(theSij[ej][ei])*theCb );
theRLast = -Complex(0.,1.)*(theQij[ei][ej]*theSb + theSij[ei][ej]*theCb);
}
//neutralinos
else {
unsigned int ei(f1ID - ParticleID::SUSY_chi_10),
ej(f2ID - ParticleID::SUSY_chi_10);
if( ei > 1 )
ei = ( ei == 13 ) ? 3 : 2;
if( ej > 1 )
ej = ( ej == 13 ) ? 3 : 2;
theLLast = Complex(0.,1.)*( conj(theQijdp[ej][ei])*theSb
- conj(theSijdp[ej][ei])*theCb );
theRLast = -Complex(0.,1.)*(theQijdp[ei][ej]*theSb - theSijdp[ei][ej]*theCb);
}
}
//charged higgs
else {
unsigned int ei(0),ej(0);
long chg(f2ID), neu(f1ID);
if( abs(neu) == ParticleID::SUSY_chi_1plus ||
abs(neu) == ParticleID::SUSY_chi_2plus ) swap(chg, neu);
ej = ( abs(chg) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
ei = neu - ParticleID::SUSY_chi_10;
if( ei > 1 ) ei = ( ei == 13 ) ? 3 : 2;
theLLast = -theQijLp[ei][ej]*theCb;
theRLast = -theQijRp[ei][ej]*theSb;
if( higgsID < 0 ) {
Complex tmp = theLLast;
theLLast = conj(theRLast);
theRLast = conj(tmp);
}
}
norm(theCoupLast);
left(theLLast);
right(theRLast);
}
diff --git a/NEWS b/NEWS
--- a/NEWS
+++ b/NEWS
@@ -1,1007 +1,1018 @@
Herwig++ News (emacs -*- outline -*- format)
+* Herwig++-2.5.2 release: 2011-mm-dd (tagged at SVN r....)
+
+** 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.
+
+** Ticket #353: Improved numerical stability in chargino decays
+
+
* 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/Tests/Makefile.am b/Tests/Makefile.am
--- a/Tests/Makefile.am
+++ b/Tests/Makefile.am
@@ -1,241 +1,253 @@
AUTOMAKE_OPTIONS = -Wno-portability
AM_LDFLAGS += -module -avoid-version -rpath /dummy/path/not/used
EXTRA_DIST = Inputs python Rivet
dist-hook:
rm -rf $(distdir)/Inputs/.svn
rm -rf $(distdir)/python/.svn
rm -rf $(distdir)/Rivet/.svn
EXTRA_LTLIBRARIES = LeptonTest.la GammaTest.la HadronTest.la DISTest.la
if WANT_LIBFASTJET
EXTRA_LTLIBRARIES += HadronJetTest.la LeptonJetTest.la
HadronJetTest_la_SOURCES = \
Hadron/VHTest.h Hadron/VHTest.cc\
Hadron/VTest.h Hadron/VTest.cc\
Hadron/HTest.h Hadron/HTest.cc
HadronJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
HadronJetTest_la_LIBADD = $(FASTJETLIBS)
LeptonJetTest_la_SOURCES = \
Lepton/TopDecay.h Lepton/TopDecay.cc
LeptonJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
LeptonJetTest_la_LIBADD = $(FASTJETLIBS)
endif
LeptonTest_la_SOURCES = \
Lepton/VVTest.h Lepton/VVTest.cc \
Lepton/VBFTest.h Lepton/VBFTest.cc \
Lepton/VHTest.h Lepton/VHTest.cc \
Lepton/FermionTest.h Lepton/FermionTest.cc
GammaTest_la_SOURCES = \
Gamma/GammaMETest.h Gamma/GammaMETest.cc \
Gamma/GammaPMETest.h Gamma/GammaPMETest.cc
DISTest_la_SOURCES = \
DIS/DISTest.h DIS/DISTest.cc
HadronTest_la_SOURCES = \
Hadron/HadronVVTest.h Hadron/HadronVVTest.cc\
Hadron/HadronVBFTest.h Hadron/HadronVBFTest.cc\
Hadron/WHTest.h Hadron/WHTest.cc\
Hadron/ZHTest.h Hadron/ZHTest.cc\
Hadron/VGammaTest.h Hadron/VGammaTest.cc\
Hadron/ZJetTest.h Hadron/ZJetTest.cc\
Hadron/WJetTest.h Hadron/WJetTest.cc\
Hadron/QQHTest.h Hadron/QQHTest.cc
REPO = $(top_builddir)/src/HerwigDefaults.rpo
HERWIG = $(top_builddir)/src/Herwig++
HWREAD = $(HERWIG) read -r $(REPO) -L $(builddir)/.libs
HWRUN = $(HERWIG) run
tests : tests-LEP tests-DIS tests-LHC tests-Gamma
if WANT_LIBFASTJET
tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons \
test-LEP-default test-LEP-Powheg test-LEP-TopDecay
else
tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons
endif
tests-DIS : test-DIS-Charged test-DIS-Neutral
if WANT_LIBFASTJET
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top test-LHC-Bottom \
test-LHC-WHJet test-LHC-ZHJet test-LHC-HJet test-LHC-ZShower test-LHC-WShower\
test-LHC-WHJet-Powheg test-LHC-ZHJet-Powheg test-LHC-HJet-Powheg \
test-LHC-ZShower-Powheg test-LHC-WShower-Powheg
else
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top
endif
tests-Anomalous : test-LHC-ZZ-Anomalous
tests-Gamma : test-Gamma-FF test-Gamma-WW test-Gamma-P
if WANT_LIBFASTJET
test-LEP-% : Inputs/LEP-%.in LeptonTest.la LeptonJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LEP-% : Inputs/LEP-%.in LeptonTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
Rivet-LEP-% : Rivet/LEP-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-TVT-% : Rivet/TVT-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-DIS-% : Rivet/DIS-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-LHC-% : Rivet/LHC-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-Star-% : Rivet/Star-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-LEP: Rivet-LEP-22 Rivet-LEP-35 Rivet-LEP-44 Rivet-LEP-91 \
Rivet-LEP-133 Rivet-LEP-161 Rivet-LEP-172 Rivet-LEP-177 \
Rivet-LEP-183 Rivet-LEP-189 Rivet-LEP-196 Rivet-LEP-197 \
Rivet-LEP-200 Rivet-LEP-206 Rivet-LEP-14 \
Rivet-LEP-Powheg-14 Rivet-LEP-Powheg-22 \
Rivet-LEP-Powheg-35 Rivet-LEP-Powheg-44 \
Rivet-LEP-Powheg-91 Rivet-LEP-Powheg-133 \
Rivet-LEP-Powheg-161 Rivet-LEP-Powheg-172 \
Rivet-LEP-Powheg-177 Rivet-LEP-Powheg-183 \
Rivet-LEP-Powheg-189 Rivet-LEP-Powheg-196 \
Rivet-LEP-Powheg-197 Rivet-LEP-Powheg-200 \
Rivet-LEP-Powheg-206
for i in LEP-*.aida; do rivet-rmgaps $$i; done;
rivet-rmgaps LEP-91.aida;
rivet-rmgaps LEP-Powheg-91.aida
rm -rf Rivet-LEP
python/merge-LEP LEP
python/merge-LEP LEP-Powheg
rivet-mkhtml -o Rivet-LEP LEP.aida:Hw++ LEP-Powheg.aida:Hw++-Powheg
Rivet-DIS: Rivet-DIS-e--LowQ2 \
- Rivet-DIS-e+-LowQ2 Rivet-DIS-e+-HighQ2
+ Rivet-DIS-e+-LowQ2 Rivet-DIS-e+-HighQ2\
+ Rivet-DIS-Powheg-e--LowQ2 \
+ Rivet-DIS-Powheg-e+-LowQ2 Rivet-DIS-Powheg-e+-HighQ2\
+ Rivet-DIS-NoME-e--LowQ2 \
+ Rivet-DIS-NoME-e+-LowQ2 Rivet-DIS-NoME-e+-HighQ2
rivet-rmgaps DIS-e+-LowQ2.aida
rivet-rmgaps DIS-e--LowQ2.aida
- rivet-rmgaps DIS-e--HighQ2.aida
- python/merge-DIS DIS
- rivet-mkhtml -o Rivet-DIS DIS.aida
+ rivet-rmgaps DIS-e+-HighQ2.aida
+ rivet-rmgaps DIS-Powheg-e+-LowQ2.aida
+ rivet-rmgaps DIS-Powheg-e--LowQ2.aida
+ rivet-rmgaps DIS-Powheg-e+-HighQ2.aida
+ rivet-rmgaps DIS-NoME-e+-LowQ2.aida
+ rivet-rmgaps DIS-NoME-e--LowQ2.aida
+ rivet-rmgaps DIS-NoME-e+-HighQ2.aida
+ python/merge-DIS DIS
+ python/merge-DIS DIS-Powheg
+ python/merge-DIS DIS-NoME
+ rivet-mkhtml -o Rivet-DIS DIS.aida:Hw++ DIS-Powheg.aida:Hw++-Powheg DIS-NoME.aida:Hw++-NoME
Rivet-TVT-WZ: Rivet-TVT-Run-I-Z Rivet-TVT-Powheg-Run-I-Z \
Rivet-TVT-Run-I-W Rivet-TVT-Powheg-Run-I-W \
Rivet-TVT-Run-I-WZ Rivet-TVT-Powheg-Run-I-WZ\
Rivet-TVT-Run-II-Z-e Rivet-TVT-Powheg-Run-II-Z-e \
Rivet-TVT-Run-II-Z-mu Rivet-TVT-Powheg-Run-II-Z-mu \
Rivet-TVT-Run-II-W Rivet-TVT-Powheg-Run-II-W
rivet-rmgaps TVT-Run-II-Z-e.aida;
rivet-rmgaps TVT-Powheg-Run-II-Z-e.aida;
rm -rf Rivet-TVT-WZ
python/merge-TVT-EW TVT-Run-II-W.aida TVT-Run-II-Z-{e,mu}.aida\
TVT-Run-I-{W,Z,WZ}.aida -o TVT-WZ.aida
python/merge-TVT-EW TVT-Powheg-Run-II-W.aida TVT-Powheg-Run-II-Z-{e,mu}.aida\
TVT-Powheg-Run-I-{W,Z,WZ}.aida -o TVT-Powheg-WZ.aida
rivet-mkhtml -o Rivet-TVT-WZ TVT-WZ.aida:Hw++ TVT-Powheg-WZ.aida:Hw++-Powheg
Rivet-TVT-Photon: Rivet-TVT-Run-II-DiPhoton Rivet-TVT-Run-II-PromptPhoton
# Rivet-TVT-Run-I-PromptPhoton
rm -rf Rivet-TVT-Photon
python/merge-aida TVT-Run-II-DiPhoton.aida TVT-Run-II-PromptPhoton.aida\
-o TVT-Photon.aida
rivet-mkhtml -o Rivet-TVT-Photon TVT-Photon.aida:Hw++
Rivet-TVT-Jets: Rivet-TVT-Run-II-Jets-1 Rivet-TVT-Run-II-Jets-2 \
Rivet-TVT-Run-II-Jets-3 Rivet-TVT-Run-II-Jets-4 \
Rivet-TVT-Run-II-Jets-5 Rivet-TVT-Run-II-Jets-6 \
Rivet-TVT-Run-II-Jets-7 Rivet-TVT-Run-II-Jets-8 \
Rivet-TVT-Run-II-Jets-9 Rivet-TVT-Run-II-Jets-10\
Rivet-TVT-Run-II-Jets-11 Rivet-TVT-Run-II-UE \
Rivet-TVT-Run-I-Jets-1 Rivet-TVT-Run-I-Jets-2 \
Rivet-TVT-Run-I-Jets-3 Rivet-TVT-Run-I-Jets-4 \
Rivet-TVT-Run-I-Jets-5 Rivet-TVT-Run-I-Jets-6 \
Rivet-TVT-Run-I-Jets-7 Rivet-TVT-Run-I-Jets-8\
Rivet-TVT-Run-I-UE\
Rivet-TVT-630-UE Rivet-TVT-630-Jets-1 \
Rivet-TVT-630-Jets-2 Rivet-TVT-630-Jets-3
rivet-rmgaps TVT-Run-I-Jets-4.aida
rm -rf Rivet-TVT-Jets
python/merge-TVT-Jets TVT
rivet-mkhtml -o Rivet-TVT-Jets TVT-Jets.aida:Hw++
Rivet-LHC-Jets: Rivet-LHC-7-Jets-1 Rivet-LHC-7-Jets-2 \
Rivet-LHC-7-Jets-3 Rivet-LHC-7-Jets-4 \
Rivet-LHC-7-Jets-5 Rivet-LHC-7-Jets-6 \
Rivet-LHC-7-Jets-7 Rivet-LHC-7-Jets-8 \
Rivet-LHC-7-Jets-9 Rivet-LHC-7-Jets-10 \
Rivet-LHC-7-Jets-11 Rivet-LHC-7-Jets-12 \
Rivet-LHC-7-Jets-13 Rivet-LHC-7-UE \
Rivet-LHC-2360-UE Rivet-LHC-900-UE
rm -rf Rivet-LHC-Jets
python/merge-LHC-Jets
rivet-mkhtml -o Rivet-LHC-Jets LHC-Jets.aida:Hw++
Rivet-Star: Rivet-Star-UE Rivet-Star-Jets-1 \
Rivet-Star-Jets-2 Rivet-Star-Jets-3 \
Rivet-Star-Jets-4
rm -rf Rivet-Star
rivet-rmgaps Star-UE.aida
python/merge-Star Star
rivet-mkhtml -o Rivet-Star Star.aida
Rivet-LHC-EW: Rivet-LHC-W-e Rivet-LHC-Powheg-W-e \
Rivet-LHC-W-mu Rivet-LHC-Powheg-W-mu \
Rivet-LHC-Z Rivet-LHC-Powheg-Z \
Rivet-LHC-WW Rivet-LHC-Powheg-WW\
Rivet-LHC-ZZ Rivet-LHC-Powheg-ZZ
rm -rf Rivet-LHC-EW;
python/merge-LHC-EW LHC-{W-e,W-mu,Z,WW,ZZ}.aida -o LHC-EW.aida;
python/merge-LHC-EW LHC-Powheg-{W-e,W-mu,Z,WW,ZZ}.aida -o LHC-Powheg-EW.aida;
rivet-mkhtml -o Rivet-LHC-EW LHC-EW.aida:Hw++ LHC-Powheg-EW.aida:Hw++-Powheg;
Rivet-LHC-Photon: Rivet-LHC-7-PromptPhoton Rivet-LHC-7-DiPhoton
rm -rf Rivet-LHC-Photon
python/merge-aida LHC-7-PromptPhoton.aida LHC-7-DiPhoton.aida -o LHC-Photon.aida
rivet-mkhtml -o Rivet-LHC-Photon LHC-Photon.aida:Hw++
test-Gamma-% : Inputs/Gamma-%.in GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
test-DIS-% : Inputs/DIS-%.in DISTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
if WANT_LIBFASTJET
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la HadronJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
clean-local:
rm -f *.out *.log *.tex *.top *.run *.dump *.mult *.Bmult *.aida
diff --git a/Tests/Rivet/DIS-NoME-e+-HighQ2.in b/Tests/Rivet/DIS-NoME-e+-HighQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-NoME-e+-HighQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-NoME.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 40.
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-NoME-e+-HighQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-NoME-e+-LowQ2.in b/Tests/Rivet/DIS-NoME-e+-LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-NoME-e+-LowQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-NoME.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-NoME-e+-LowQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-NoME-e--LowQ2.in b/Tests/Rivet/DIS-NoME-e--LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-NoME-e--LowQ2.in
@@ -0,0 +1,26 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-NoME.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e-
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 26.7*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1994_S2919893
+#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1995_S3167097
+
+
+cd /Herwig/Generators
+saverun DIS-NoME-e--LowQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-Powheg-e+-HighQ2.in b/Tests/Rivet/DIS-Powheg-e+-HighQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-Powheg-e+-HighQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-Powheg.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 40.
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-Powheg-e+-HighQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-Powheg-e+-LowQ2.in b/Tests/Rivet/DIS-Powheg-e+-LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-Powheg-e+-LowQ2.in
@@ -0,0 +1,24 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-Powheg.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 27.5*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_2000_S4129130
+
+cd /Herwig/Generators
+saverun DIS-Powheg-e+-LowQ2 DISGenerator
diff --git a/Tests/Rivet/DIS-Powheg-e--LowQ2.in b/Tests/Rivet/DIS-Powheg-e--LowQ2.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DIS-Powheg-e--LowQ2.in
@@ -0,0 +1,26 @@
+##################################################
+# Rivet analyses at HERA
+##################################################
+read DISBase-Powheg.in
+
+##################################################
+# DIS parameters
+##################################################
+cd /Herwig/Generators
+set /Herwig/EventHandlers/DISHandler:BeamB /Herwig/Particles/p+
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxB 820.*GeV
+set /Herwig/EventHandlers/DISHandler:BeamA /Herwig/Particles/e-
+set /Herwig/EventHandlers/DISLuminosity:BeamEMaxA 26.7*GeV
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 1000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxW2 1000000.
+set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 60.
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 2.5
+set /Herwig/Cuts/DISCuts:MHatMin 0.
+
+# H1 energy flow analysis
+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1994_S2919893
+#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 H1_1995_S3167097
+
+
+cd /Herwig/Generators
+saverun DIS-Powheg-e--LowQ2 DISGenerator
diff --git a/Tests/Rivet/DISBase-NoME.in b/Tests/Rivet/DISBase-NoME.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DISBase-NoME.in
@@ -0,0 +1,24 @@
+##################################################
+# Example generator based on DIS parameters
+# usage: Herwig++ read DIS.in
+##################################################
+
+set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
+set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
+
+cd /Herwig/MatrixElements/
+# Neutral current DIS
+insert SimpleDIS:MatrixElements[0] MEDISNC
+
+cd /Herwig/Generators
+set DISGenerator:NumberOfEvents 10000000
+set DISGenerator:RandomNumberGenerator:Seed 31122001
+set DISGenerator:DebugLevel 0
+set DISGenerator:PrintEvent 10
+set DISGenerator:MaxErrors 1000000
+set /Herwig/Shower/ShowerHandler:MPIHandler NULL
+set /Herwig/Shower/Evolver:MECorrMode 0
+cd /Herwig/Generators
+
+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
+insert DISGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
diff --git a/Tests/Rivet/DISBase-Powheg.in b/Tests/Rivet/DISBase-Powheg.in
new file mode 100644
--- /dev/null
+++ b/Tests/Rivet/DISBase-Powheg.in
@@ -0,0 +1,43 @@
+##################################################
+# Example generator based on DIS parameters
+# usage: Herwig++ read DIS.in
+##################################################
+
+set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
+set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
+
+##################################################
+# Need to use an NLO PDF
+##################################################
+set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+
+##################################################
+# Setup the POWHEG shower
+##################################################
+cd /Herwig/Shower
+set Evolver:IntrinsicPtGaussian 1.9*GeV
+set Evolver:HardEmissionMode POWHEG
+
+##################################################
+# and strong coupling
+##################################################
+create Herwig::O2AlphaS O2AlphaS
+set /Herwig/Model:QCD/RunningAlphaS O2AlphaS
+
+cd /Herwig/MatrixElements/
+# Neutral current DIS
+
+insert SimpleDIS:MatrixElements[0] /Herwig/MatrixElements/PowhegMEDISNC
+
+cd /Herwig/Generators
+set DISGenerator:NumberOfEvents 10000000
+set DISGenerator:RandomNumberGenerator:Seed 31122001
+set DISGenerator:DebugLevel 0
+set DISGenerator:PrintEvent 10
+set DISGenerator:MaxErrors 1000000
+set /Herwig/Shower/ShowerHandler:MPIHandler NULL
+cd /Herwig/Generators
+
+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
+insert DISGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,170 +1,170 @@
dnl Process this file with autoconf to produce a configure script.
AC_PREREQ([2.59])
AC_INIT([Herwig++],[SVN],[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
case "${host}" in
*-darwin[[0156]].*)
AC_MSG_ERROR([Herwig++ requires OS X 10.3 or later])
;;
*-darwin7.*)
if test "x$MACOSX_DEPLOYMENT_TARGET" != "x10.3"; then
AC_MSG_ERROR(
[Please add 'MACOSX_DEPLOYMENT_TARGET=10.3' to the configure line.])
fi
;;
esac
dnl === disable debug symbols by default =====
if test "x$CXXFLAGS" = "x"; then
CXXFLAGS=-O3
fi
if test "x$CFLAGS" = "x"; then
CFLAGS=-O3
fi
dnl Looptools manual requires optimization off
if test "x$FCFLAGS" = "x"; then
FCFLAGS=-O0
fi
dnl ==========================================
AC_LANG([C++])
AM_INIT_AUTOMAKE([1.9 gnu dist-bzip2 -Wall])
m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
dnl Checks for programs.
AC_PROG_CXX([g++])
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])
LT_INIT([disable-static dlopen pic-only])
dnl ####################################
dnl ####################################
dnl for Doc/fixinterfaces.pl
AC_PATH_PROG(PERL, perl)
HERWIG_CHECK_GSL
HERWIG_CHECK_THEPEG
HERWIG_COMPILERFLAGS
HERWIG_LOOPTOOLS
HERWIG_PDF_PATH
-HERWIG_CHECK_FASTJET
+FASTJET_CHECK_FASTJET
HERWIG_VERSIONSTRING
HERWIG_CHECK_ABS_BUG
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/Susy/MatrixElements/Makefile
Models/UED/Makefile
Models/LH/Makefile
Models/LHTP/Makefile
Models/HiddenValley/Makefile
Models/Transplanckian/Makefile
Models/Leptoquarks/Makefile
Models/AnomalousCouplings/Makefile
Models/Zprime/Makefile
Models/Feynrules/Makefile
Models/Feynrules/Makefile-FR
Models/ADD/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
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
Shower/SplittingFunctions/Makefile
Shower/Default/Makefile
Shower/Base/Makefile
Shower/Makefile
Utilities/Makefile
Hadronization/Makefile
lib/Makefile
include/Makefile
src/Makefile
src/defaults/Makefile
src/herwig-config
Doc/Makefile
Doc/HerwigDefaults.in
Looptools/Makefile
Analysis/Makefile
src/Makefile-UserModules
src/defaults/Analysis.in
Contrib/Makefile
Contrib/make_makefiles.sh
Tests/Makefile
Makefile])
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/m4/fastjet.m4 b/m4/fastjet.m4
new file mode 100644
--- /dev/null
+++ b/m4/fastjet.m4
@@ -0,0 +1,80 @@
+dnl CHECK FASTJET BEGIN
+dnl
+dnl This script can be used in configure scripts to check for the
+dnl usability of the FastJet librarty.
+dnl
+dnl By defaults, it searches the FastJet library in standard system
+dnl locations but an alternative path can be specified using the
+dnl --with-fastjet=... configure option
+dnl
+dnl If FastJet is found and functional, the variables FASTJET_CXXFLAGS
+dnl and FASTJET_LIBS are set
+dnl
+dnl modified for Herwig++ 2011-10-04 D.Grellscheid
+dnl
+AC_DEFUN([FASTJET_CHECK_FASTJET],
+[
+dnl ckeck if a directory is specified for FastJet
+AC_ARG_WITH(fastjet,
+ [AC_HELP_STRING([--with-fastjet=dir],
+ [Assume the given directory for FastJet])])
+
+dnl search for the fastjet-config script
+if test "$with_fastjet" = ""; then
+ AC_PATH_PROG(fjconfig, fastjet-config, no)
+else
+ AC_PATH_PROG(fjconfig, fastjet-config, no, ${with_fastjet}/bin)
+fi
+
+LOAD_FASTJET=""
+CREATE_FASTJET="#create"
+
+if test "${fjconfig}" = "no"; then
+ AC_MSG_CHECKING(FastJet)
+ AC_MSG_RESULT(no);
+ $2
+else
+
+ dnl now see if FastJet is functional
+ save_CXXFLAGS="$CXXFLAGS"
+ save_LIBS="$LIBS"
+
+ CXXFLAGS="${CXXFLAGS} `${fjconfig} --cxxflags`"
+ LIBS="${LIBS} `${fjconfig} --libs --plugins`"
+
+ AC_MSG_CHECKING([if FastJet is functional])
+ AC_LANG_PUSH(C++)
+ AC_COMPILE_IFELSE(AC_LANG_PROGRAM([[
+#include <fastjet/ClusterSequence.hh>
+ ]], [[
+fastjet::PseudoJet pj=fastjet::PtYPhiM(10.0,0.5,1.0,0.0);
+ ]]), [fjok='yes'], [fjok='no'])
+ AC_MSG_RESULT([$fjok])
+ AC_LANG_POP()
+ CXXFLAGS="$save_CXXFLAGS"
+ LIBS="$save_LIBS"
+
+ AC_MSG_CHECKING(FastJet)
+ if test "${fjok}" = "yes"; then
+ FASTJET_CXXFLAGS="`${fjconfig} --cxxflags`"
+ FASTJET_LIBS="`${fjconfig} --libs --plugins`"
+ LOAD_FASTJET="library HwLEPJetAnalysis.so"
+ CREATE_FASTJET="create"
+ AC_MSG_RESULT(yes)
+ $1
+ else
+ AC_MSG_RESULT(no)
+ $2
+ fi
+fi
+
+
+AC_SUBST(FASTJETINCLUDE,[$FASTJET_CXXFLAGS])
+AC_SUBST(CREATE_FASTJET)
+AC_SUBST(LOAD_FASTJET)
+AC_SUBST(FASTJETLIBS,[${FASTJET_LIBS/-Wl,-rpath,/-R}])
+
+AM_CONDITIONAL(WANT_LIBFASTJET,[test "x$CREATE_FASTJET" = "xcreate"])
+])
+
+dnl CHECK FASTJET END
diff --git a/m4/herwig.m4 b/m4/herwig.m4
--- a/m4/herwig.m4
+++ b/m4/herwig.m4
@@ -1,468 +1,399 @@
# check for gcc bug http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34130
AC_DEFUN([HERWIG_CHECK_ABS_BUG],
[
AC_REQUIRE([HERWIG_COMPILERFLAGS])
if test "$GCC" = "yes"; then
AC_MSG_CHECKING([for gcc abs bug])
AC_RUN_IFELSE([
AC_LANG_PROGRAM(
[[ int foo (int i) { return -2 * __builtin_abs(i - 2); } ]],
[[ if ( foo(1) != -2 || foo(3) != -2 ) return 1; ]]
)],
[ AC_MSG_RESULT([not found. Compiler is ok.]) ],
[
AC_MSG_RESULT([found. Builtin abs() is buggy.])
AC_MSG_CHECKING([if -fno-builtin-abs works])
oldcxxflags=$CXXFLAGS
CXXFLAGS="$CXXFLAGS -fno-builtin-abs"
AC_RUN_IFELSE([
AC_LANG_PROGRAM(
[[
#include <cstdlib>
int foo (int i) { return -2 * std::abs(i - 2); }
]],
[[
if (foo(1) != -2 || foo(3) != -2) return 1;
]]
)],
[
AC_MSG_RESULT([yes. Setting -fno-builtin-abs.])
AM_CXXFLAGS="$AM_CXXFLAGS -fno-builtin-abs"
AM_CFLAGS="$AM_CFLAGS -fno-builtin-abs"
],
[
AC_MSG_RESULT([no. Setting -fno-builtin.])
AC_MSG_WARN([
*****************************************************************************
For this version of gcc, -fno-builtin-abs alone did not work to avoid the
gcc abs() bug. Instead, all gcc builtin functions are now disabled.
Update gcc if possible.
*****************************************************************************])
AM_CXXFLAGS="$AM_CXXFLAGS -fno-builtin"
AM_CFLAGS="$AM_CFLAGS -fno-builtin"
]
)
CXXFLAGS=$oldcxxflags
]
)
fi
])
dnl ##### THEPEG #####
AC_DEFUN([HERWIG_CHECK_THEPEG],
[
defaultlocation="${prefix}"
test "x$defaultlocation" = xNONE && defaultlocation="${ac_default_prefix}"
AC_MSG_CHECKING([for libThePEG in])
AC_ARG_WITH(thepeg,
AC_HELP_STRING([--with-thepeg=DIR],[location of ThePEG installation]),
[],
[with_thepeg="${defaultlocation}"])
AC_MSG_RESULT([$with_thepeg])
if test "x$with_thepeg" = "xno"; then
AC_MSG_ERROR([Cannot build Herwig++ without ThePEG. Please set --with-thepeg.])
fi
THEPEGLDFLAGS="-L${with_thepeg}/lib/ThePEG"
THEPEGPATH="${with_thepeg}"
oldldflags="$LDFLAGS"
oldlibs="$LIBS"
LDFLAGS="$LDFLAGS $THEPEGLDFLAGS"
AC_CHECK_LIB([ThePEG],[debugThePEG],[],
[AC_MSG_ERROR([No ThePEG libraries in $THEPEGLDFLAGS. Please set --with-thepeg.])])
AC_SUBST([THEPEGLIB],[-lThePEG])
AC_SUBST(THEPEGLDFLAGS)
AC_SUBST(THEPEGPATH)
LIBS="$oldlibs"
LDFLAGS="$oldldflags"
AC_MSG_CHECKING([for ThePEG headers in])
AC_ARG_WITH([thepeg-headers],
AC_HELP_STRING([--with-thepeg-headers=DIR],[location of ThePEG include directory]),
[],
[with_thepeg_headers="${with_thepeg}/include"])
AC_MSG_RESULT([$with_thepeg_headers])
if test "x$with_thepeg_headers" = "xno"; then
AC_MSG_ERROR([Cannot build Herwig++ without ThePEG headers. Please set --with-thepeg-headers.])
fi
THEPEGINCLUDE="-I$with_thepeg_headers"
oldcppflags="$CPPFLAGS"
CPPFLAGS="$CPPFLAGS $THEPEGINCLUDE"
AC_CHECK_HEADER([ThePEG/Config/ThePEG.h],[],
[AC_MSG_ERROR([No ThePEG headers in $with_thepeg_headers. Please set --with-thepeg-headers.])])
CPPFLAGS="$oldcppflags"
AC_SUBST(THEPEGINCLUDE)
AC_MSG_CHECKING([for HepMCAnalysis.so in ThePEG])
if test -x "$THEPEGPATH/lib/ThePEG/HepMCAnalysis.so" ; then
CREATE_HEPMC="create"
AC_MSG_RESULT([found])
else
CREATE_HEPMC="# create"
AC_MSG_RESULT([not found])
fi
AC_SUBST([CREATE_HEPMC])
])
-dnl ##### FastJet #####
-AC_DEFUN([HERWIG_CHECK_FASTJET],[
-
-FASTJETPATH=""
-FASTJETLIBS=""
-FASTJETINCLUDE=""
-LOAD_FASTJET=""
-CREATE_FASTJET="#create"
-
-AC_MSG_CHECKING([for FastJet location])
-
-AC_ARG_WITH(fastjet,
- AC_HELP_STRING([--with-fastjet=DIR],[Location of FastJet installation @<:@default=system libs@:>@]),
- [],
- [with_fastjet=system])
-
-if test "x$with_fastjet" = "xno"; then
- AC_MSG_RESULT([FastJet support disabled.])
-elif test "x$with_fastjet" = "xsystem"; then
- AC_MSG_RESULT([in system libraries])
- oldlibs="$LIBS"
- AC_CHECK_LIB(fastjet,main,
- [],
- [with_fastjet=no
- AC_MSG_WARN([
-FastJet not found in system libraries])
- ])
- FASTJETLIBS="$LIBS"
- LIBS=$oldlibs
-else
- AC_MSG_RESULT([$with_fastjet])
- FASTJETPATH="$with_fastjet"
- FASTJETINCLUDE="-I$with_fastjet/include"
- FASTJETLIBS="-L$with_fastjet/lib -R$with_fastjet/lib -lfastjet"
-fi
-
-if test "x$with_fastjet" != "xno"; then
- # Now lets see if the libraries work properly
- oldLIBS="$LIBS"
- oldLDFLAGS="$LDFLAGS"
- oldCPPFLAGS="$CPPFLAGS"
- LIBS="$LIBS `echo $FASTJETLIBS | sed -e 's!-R.* ! !'`"
- LDFLAGS="$LDFLAGS"
- CPPFLAGS="$CPPFLAGS $FASTJETINCLUDE"
-
- AC_MSG_CHECKING([that FastJet works])
- AC_LINK_IFELSE([AC_LANG_PROGRAM([[#include <fastjet/ClusterSequence.hh>
-]],[[fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm, fastjet::E_scheme, fastjet::Best); ]])],
- [AC_MSG_RESULT([yes])],[AC_MSG_RESULT([no])
- AC_MSG_ERROR([Use '--with-fastjet=' to set a path or use '--without-fastjet'.])
- ])
-
- AC_CHECK_HEADERS([fastjet/ClusterSequence.hh])
-
- LIBS="$oldLIBS"
- LDFLAGS="$oldLDFLAGS"
- CPPFLAGS="$oldCPPFLAGS"
-
- CREATE_FASTJET="create"
- LOAD_FASTJET="library HwLEPJetAnalysis.so"
-fi
-AC_SUBST(FASTJETINCLUDE)
-AC_SUBST(CREATE_FASTJET)
-AC_SUBST(LOAD_FASTJET)
-AC_SUBST(FASTJETLIBS)
-
-AM_CONDITIONAL(WANT_LIBFASTJET,[test "x$CREATE_FASTJET" = "xcreate"])
-])
-
dnl ##### LOOPTOOLS #####
AC_DEFUN([HERWIG_LOOPTOOLS],
[
AC_REQUIRE([AC_PROG_FC])
AC_REQUIRE([AC_FC_LIBRARY_LDFLAGS])
AC_REQUIRE([AC_PROG_CC])
AC_REQUIRE([HERWIG_COMPILERFLAGS])
AC_MSG_CHECKING([if Looptools build works])
enable_looptools=yes
if test "x$GCC" = "xyes"; then
case "${host}" in
x86_64-*|*-darwin1*)
AM_FCFLAGS="$AM_FCFLAGS -fdefault-integer-8"
;;
esac
AC_LANG_PUSH([Fortran])
oldFCFLAGS="$FCFLAGS"
FCFLAGS="$AM_FCFLAGS"
AC_COMPILE_IFELSE(
AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
[],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([needs gfortran on 64bit machines])]
)
FCFLAGS="$oldFCFLAGS"
AC_LANG_POP([Fortran])
fi
AC_MSG_RESULT([$enable_looptools])
AC_SUBST([F77],[$FC])
AC_SUBST([FFLAGS],[$FCFLAGS])
AC_SUBST([AM_FFLAGS],[$AM_FCFLAGS])
AC_SUBST([FLIBS],[$FCLIBS])
])
dnl ##### PDF PATH #####
AC_DEFUN([HERWIG_PDF_PATH],
[
AC_MSG_CHECKING([which Herwig++ PDF data to use])
AC_ARG_WITH(pdf,
AC_HELP_STRING([--with-pdf=DIR],[installation path of Herwig++PDF data tarball]),
[],
[with_pdf=${prefix}]
)
HERWIG_PDF_PREFIX=${with_pdf}/share/Herwig++PDF
if test -f "${HERWIG_PDF_PREFIX}/mrst/2008/mrstMCal.dat"; then
AC_MSG_RESULT([$with_pdf])
localPDFneeded=false
else
AC_MSG_RESULT([Using built-in PDF data set. For other data sets, set --with-pdf.])
HERWIG_PDF_PREFIX=PDF
localPDFneeded=true
fi
HERWIG_PDF_DEFAULT=${HERWIG_PDF_PREFIX}/mrst/2008/mrstMCal.dat
HERWIG_PDF_NLO=${HERWIG_PDF_PREFIX}/mrst/2002/mrst2002nlo.dat
HERWIG_PDF_POMERON=${HERWIG_PDF_PREFIX}/diffraction/
AM_CONDITIONAL(WANT_LOCAL_PDF,[test "x$localPDFneeded" = "xtrue"])
AC_SUBST(HERWIG_PDF_DEFAULT)
AC_SUBST(HERWIG_PDF_NLO)
AC_SUBST(HERWIG_PDF_POMERON)
])
dnl ###### GSL ######
AC_DEFUN([HERWIG_CHECK_GSL],
[
AC_MSG_CHECKING([for gsl location])
GSLINCLUDE=""
GSLLIBS=""
AC_ARG_WITH(gsl,
AC_HELP_STRING([--with-gsl=DIR],[location of gsl installation @<:@default=system libs@:>@]),
[],
[with_gsl=system])
if test "x$with_gsl" = "xno"; then
AC_MSG_ERROR([libgsl is required. Please install the GNU scientific library and header files.])
fi
if test "x$with_gsl" = "xsystem"; then
AC_MSG_RESULT([in system libraries])
oldlibs="$LIBS"
AC_CHECK_LIB(m,main)
AC_CHECK_LIB(gslcblas,main)
AC_CHECK_LIB(gsl,main,[],
[
AC_MSG_ERROR([Cannot find libgsl. Please install the GNU scientific library and header files or use --with-gsl=.])
]
)
GSLLIBS="$LIBS"
LIBS=$oldlibs
else
if test "`uname -m`" = "x86_64" -a -e "$with_gsl/lib64/libgsl.a" -a -d "$with_gsl/include/gsl"; then
AC_MSG_RESULT([found in $with_gsl])
GSLLIBS="-L$with_gsl/lib64 -R$with_gsl/lib64 -lgslcblas -lgsl"
GSLINCLUDE="-I$with_gsl/include"
elif test -e "$with_gsl/lib/libgsl.a" -a -d "$with_gsl/include/gsl"; then
AC_MSG_RESULT([found in $with_gsl])
GSLLIBS="-L$with_gsl/lib -R$with_gsl/lib -lgslcblas -lgsl"
GSLINCLUDE="-I$with_gsl/include"
else
AC_MSG_RESULT([not found])
AC_MSG_ERROR([Can't find $with_gsl/lib/libgsl.a or the headers in $with_gsl/include])
fi
fi
AC_SUBST(GSLINCLUDE)
AC_SUBST(GSLLIBS)
])
AC_DEFUN([HERWIG_VERSIONSTRING],
[
if test -d $srcdir/.svn; then
AC_CHECK_PROG(have_svnversion,[svnversion],[yes],[no])
fi
AM_CONDITIONAL(USE_SVNVERSION,[test "x$have_svnversion" = "xyes"])
])
dnl ##### COMPILERFLAGS #####
AC_DEFUN([HERWIG_COMPILERFLAGS],
[
AC_REQUIRE([HERWIG_CHECK_THEPEG])
AM_CPPFLAGS="-I\$(top_builddir)/include $THEPEGINCLUDE \$(GSLINCLUDE)"
AC_MSG_CHECKING([for debugging mode])
AC_ARG_ENABLE(debug,
AC_HELP_STRING([--enable-debug],[debug mode, use --enable-debug=slow for additional options that slow down the run.]),
[],
[enable_debug=no]
)
AC_MSG_RESULT([$enable_debug])
if test "x$enable_debug" = "xno"; then
AM_CPPFLAGS="$AM_CPPFLAGS -DNDEBUG"
else
debugflags="-g"
fi
dnl -Wfloat-equal -fvisibility-inlines-hidden -Wctor-dtor-privacy -Weffc++
if test -n "$GCC"; then
warnflags="-ansi -pedantic -Wall -W"
if test "x$enable_debug" = "xslow"; then
debugflags="$debugflags -fno-inline"
AM_CPPFLAGS="$AM_CPPFLAGS -D_GLIBCXX_DEBUG"
fi
fi
dnl do an actual capability check on ld instead of this workaround
case "${host}" in
*-darwin*)
;;
*)
AM_LDFLAGS="-Wl,--enable-new-dtags"
;;
esac
AC_SUBST(AM_CPPFLAGS)
AC_SUBST(AM_CFLAGS, ["$warnflags $debugflags"])
AC_SUBST(AM_CXXFLAGS,["$warnflags $debugflags"])
AC_SUBST(AM_FCFLAGS, ["$debugflags"])
AC_SUBST(AM_LDFLAGS)
])
AC_DEFUN([HERWIG_ENABLE_MODELS],
[
AC_MSG_CHECKING([for BSM models to include])
LOAD_RS=""
LOAD_SUSY=""
LOAD_NMSSM=""
LOAD_TRP=""
LOAD_UED=""
LOAD_ADD=""
AC_ARG_ENABLE(models,
AC_HELP_STRING([--enable-models=LIST],[Comma-separated list of BSM models to enable. Options are (mssm nmssm ued rs add lh lhtp rpv trp anomalous leptoquarks) or --disable-models to turn them all off.]),
[],
[enable_models=all]
)
if test "x$enable_models" = "xyes" -o "x$enable_models" = "xall"; then
all=yes
fi
AC_MSG_RESULT([$enable_models])
if test ! "$all"; then
oldIFS="$IFS"
IFS=","
for i in $enable_models; do
declare $i=yes
done
IFS="$oldIFS"
fi
if test "$nmssm" -o "$all"; then
LOAD_NMSSM="library HwNMSSM.so"
mssm=yes
fi
AC_SUBST(LOAD_NMSSM)
if test "$rpv"; then
mssm=yes
fi
if test "$rs" -o "$all" ; then
LOAD_RS="library HwRSModel.so"
fi
AC_SUBST(LOAD_RS)
if test "$mssm" -o "$all"; then
LOAD_SUSY="library HwSusy.so"
fi
AC_SUBST(LOAD_SUSY)
if test "$trp" -o "$all"; then
LOAD_TRP="library HwTransplanck.so"
fi
AC_SUBST(LOAD_TRP)
if test "$ued" -o "$all"; then
LOAD_UED="library HwUED.so"
fi
AC_SUBST(LOAD_UED)
if test "$add" -o "$all"; then
LOAD_ADD="library HwADDModel.so"
fi
AC_SUBST(LOAD_ADD)
if test "$leptoquarks" -o "$all"; then
LOAD_LEPTOQUARKS="library HwLeptoquarkModel.so"
fi
AC_SUBST(LOAD_LEPTOQUARKS)
AM_CONDITIONAL(WANT_MSSM,[test "$mssm" -o "$all"])
AM_CONDITIONAL(WANT_NMSSM,[test "$nmssm" -o "$all"])
AM_CONDITIONAL(WANT_RPV,[test "$rpv" -o "$all"])
AM_CONDITIONAL(WANT_UED,[test "$ued" -o "$all"])
AM_CONDITIONAL(WANT_RS,[test "$rs" -o "$all"])
AM_CONDITIONAL(WANT_LH,[test "$lh" -o "$all"])
AM_CONDITIONAL(WANT_LHTP,[test "$lhtp" -o "$all"])
AM_CONDITIONAL(WANT_Anomalous,[test "$anomalous" -o "$all"])
AM_CONDITIONAL(WANT_Leptoquark,[test "$leptoquarks" -o "$all"])
AM_CONDITIONAL(WANT_TRP,[test "$trp" -o "$all"])
AM_CONDITIONAL(WANT_ADD,[test "$add" -o "$all"])
])
AC_DEFUN([HERWIG_OVERVIEW],
[
FCSTRING=`$FC --version | head -1`
CXXSTRING=`$CXX --version | head -1`
CCSTRING=`$CC --version | head -1`
cat << _HW_EOF_ > config.herwig
*****************************************************
*** $PACKAGE_STRING configuration summary
*** Please include this information in bug reports!
***--------------------------------------------------
*** Prefix: $prefix
***
*** BSM models: $enable_models
*** Herwig debug mode: $enable_debug
***
*** GSL: $with_gsl
***
*** ThePEG: $with_thepeg
*** ThePEG headers: $with_thepeg_headers
***
-*** Fastjet: $with_fastjet
+*** Fastjet: ${fjconfig}
***
*** Host: $host
*** CXX: $CXXSTRING
*** FC: $FCSTRING
*** CC: $CCSTRING
*****************************************************
_HW_EOF_
])
diff --git a/src/DIS.in b/src/DIS.in
--- a/src/DIS.in
+++ b/src/DIS.in
@@ -1,64 +1,100 @@
##################################################
# Example generator based on DIS parameters
# usage: Herwig++ read DIS.in
##################################################
##################################################
# Technical parameters for this run
##################################################
cd /Herwig/Generators
set DISGenerator:NumberOfEvents 10000000
set DISGenerator:RandomNumberGenerator:Seed 31122001
set DISGenerator:PrintEvent 10
set DISGenerator:MaxErrors 10000
set /Herwig/Shower/ShowerHandler:MPIHandler NULL
##################################################
# DIS physics parameters (override defaults here)
##################################################
##################################################
# Matrix Elements for lepton-hadron collisions
# (by default only neutral-current switched on)
##################################################
cd /Herwig/MatrixElements/
# Neutral current DIS
insert SimpleDIS:MatrixElements[0] MEDISNC
+# Charged current matrix element
+insert SimpleDIS:MatrixElements[0] MEDISCC
+
+##################################################
+# NLO IN POWHEG SCHEME
+##################################################
+
+###################################################
+## Need to use an NLO PDF
+###################################################
+#set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+#set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+#
+###################################################
+## Setup the POWHEG shower
+###################################################
+#cd /Herwig/Shower
+#set Evolver:HardEmissionMode POWHEG
+#
+##################################################
+# and strong coupling
+##################################################
+#create Herwig::O2AlphaS O2AlphaS
+#set /Herwig/Model:QCD/RunningAlphaS O2AlphaS
+#
+###################################################
+## NLO Matrix Elements for lepton-hadron collisions
+## in the POWHEG approach
+###################################################
+#
+#cd /Herwig/MatrixElements/
+#
+## Neutral current DIS
+#insert SimpleDIS:MatrixElements[0] PowhegMEDISNC
+## Charged current matrix element
+#insert SimpleDIS:MatrixElements[0] PowhegMEDISCC
cd /Herwig/Generators
##################################################
# Useful analysis handlers for lepton-hadron physics
##################################################
##################################################
# Useful analysis handlers for HepMC related output
##################################################
# Schematic overview of an event (requires --with-hepmc to be set at configure time
# and the graphviz program 'dot' to produce a plot)
# insert DISGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot
# A HepMC dump file (requires --with-hepmc to be set at configure time)
# insert DISGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
# set /Herwig/Analysis/HepMCFile:PrintEvent 100
# set /Herwig/Analysis/HepMCFile:Format GenEvent
# set /Herwig/Analysis/HepMCFile:Units GeV_mm
set /Herwig/Shower/KinematicsReconstructor:ReconstructionOption 1
##################################################
# Save run for later usage with 'Herwig++ run'
##################################################
saverun DIS DISGenerator
##################################################
# uncomment this section for an example batch run
# of two repeats with different parameters
#
# Note that a separate call of 'Herwig run'
# is not required in this case
##################################################
# set DISGenerator:NumberOfEvents 10
# run DIS-full DISGenerator
#
# run DIS-initial DISGenerator
diff --git a/src/defaults/MatrixElements.in b/src/defaults/MatrixElements.in
--- a/src/defaults/MatrixElements.in
+++ b/src/defaults/MatrixElements.in
@@ -1,215 +1,228 @@
##############################################################################
# Setup of default matrix elements.
#
# Only one ME is activated by default, but this file lists
# some alternatives. All available MEs can be found in the
# 'include/Herwig++/MatrixElements' subdirectory of your Herwig++
# installation.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Instead of editing this file directly, you should reset
# the matrix elements in your own input files:
#
# - create your custom SubProcessHandler
# - insert the MEs you need
# - set your SubProcessHandler instead of the default (see HerwigDefaults.in)
##############################################################################
mkdir /Herwig/MatrixElements
cd /Herwig/MatrixElements
library HwMELepton.so
library HwMEHadron.so
library HwMEDIS.so
############################################################
# e+e- matrix elements
############################################################
# e+e- > q qbar
create Herwig::MEee2gZ2qq MEee2gZ2qq
newdef MEee2gZ2qq:MinimumFlavour 1
newdef MEee2gZ2qq:MaximumFlavour 5
newdef MEee2gZ2qq:Coupling /Herwig/Shower/AlphaQCD
# e+e- -> l+l-
create Herwig::MEee2gZ2ll MEee2gZ2ll
newdef MEee2gZ2ll:Allowed Charged
# e+e- -> W+W- ZZ
create Herwig::MEee2VV MEee2VV
# e+e- -> ZH
create Herwig::MEee2ZH MEee2ZH
newdef MEee2ZH:Coupling /Herwig/Shower/AlphaQCD
# e+e- -> e+e-H/nu_enu_ebarH
create Herwig::MEee2HiggsVBF MEee2HiggsVBF
############################################################
# NLO (POWHEG e+e- matrix elements
############################################################
library HwPowhegMELepton.so
create Herwig::MEee2gZ2qqPowheg PowhegMEee2gZ2qq
newdef PowhegMEee2gZ2qq:MinimumFlavour 1
newdef PowhegMEee2gZ2qq:MaximumFlavour 5
newdef PowhegMEee2gZ2qq:Coupling /Herwig/Shower/AlphaQCD
############################################################
# hadron-hadron matrix elements
############################################################
###################################
# Electroweak processes
###################################
# q qbar -> gamma/Z -> l+l-
create Herwig::MEqq2gZ2ff MEqq2gZ2ff
newdef MEqq2gZ2ff:Process 3
newdef MEqq2gZ2ff:Coupling /Herwig/Shower/AlphaQCD
# q qbar to W -> l nu
create Herwig::MEqq2W2ff MEqq2W2ff
newdef MEqq2W2ff:Process 2
newdef MEqq2W2ff:Coupling /Herwig/Shower/AlphaQCD
# W+jet
create Herwig::MEPP2WJet MEWJet
newdef MEWJet:WDecay Leptons
# Z+jet
create Herwig::MEPP2ZJet MEZJet
newdef MEZJet:ZDecay ChargedLeptons
# PP->WW/WZ/ZZ
create Herwig::MEPP2VV MEPP2VV
# PP->WZ gamma
create Herwig::MEPP2VGamma MEPP2VGamma
###################################
# Photon and jet processes
###################################
# qqbar/gg -> gamma gamma
create Herwig::MEPP2GammaGamma MEGammaGamma
# hadron-hadron to gamma+jet
create Herwig::MEPP2GammaJet MEGammaJet
# QCD 2-to-2
create Herwig::MEQCD2to2 MEQCD2to2
# MinBias
create Herwig::MEMinBias MEMinBias
# qqbar/gg -> t tbar
create Herwig::MEPP2QQ MEHeavyQuark
###################################
# Higgs processes
###################################
# hadron-hadron to higgs
create Herwig::MEPP2Higgs MEHiggs
newdef MEHiggs:ShapeScheme MassGenerator
newdef MEHiggs:Process gg
newdef MEHiggs:Coupling /Herwig/Shower/AlphaQCD
# hadron-hadron to higgs+jet
create Herwig::MEPP2HiggsJet MEHiggsJet
# PP->ZH
create Herwig::MEPP2ZH MEPP2ZH
newdef MEPP2ZH:Coupling /Herwig/Shower/AlphaQCD
# PP->WH
create Herwig::MEPP2WH MEPP2WH
newdef MEPP2WH:Coupling /Herwig/Shower/AlphaQCD
# PP -> Higgs via VBF
create Herwig::MEPP2HiggsVBF MEPP2HiggsVBF
# PP -> t tbar Higgs
create Herwig::MEPP2QQHiggs MEPP2ttbarH
newdef MEPP2ttbarH:QuarkType Top
# PP -> b bbar Higgs
create Herwig::MEPP2QQHiggs MEPP2bbbarH
newdef MEPP2bbbarH:QuarkType Bottom
##########################################################
# Hadron-Hadron NLO matrix elements in the Powheg scheme
##########################################################
library HwPowhegMEHadron.so
# q qbar -> gamma/Z -> l+l-
create Herwig::MEqq2gZ2ffPowheg PowhegMEqq2gZ2ff
newdef PowhegMEqq2gZ2ff:Process 3
newdef PowhegMEqq2gZ2ff:Coupling /Herwig/Shower/AlphaQCD
# q qbar to W -> l nu
create Herwig::MEqq2W2ffPowheg PowhegMEqq2W2ff
newdef PowhegMEqq2W2ff:Process 2
newdef PowhegMEqq2W2ff:Coupling /Herwig/Shower/AlphaQCD
# PP->ZH
create Herwig::MEPP2ZHPowheg PowhegMEPP2ZH
newdef PowhegMEPP2ZH:Coupling /Herwig/Shower/AlphaQCD
# PP->WH
create Herwig::MEPP2WHPowheg PowhegMEPP2WH
newdef PowhegMEPP2WH:Coupling /Herwig/Shower/AlphaQCD
# hadron-hadron to higgs
create Herwig::MEPP2HiggsPowheg PowhegMEHiggs
newdef PowhegMEHiggs:ShapeScheme MassGenerator
newdef PowhegMEHiggs:Process gg
newdef PowhegMEHiggs:Coupling /Herwig/Shower/AlphaQCD
# PP->VV
create Herwig::MEPP2VVPowheg PowhegMEPP2VV
newdef PowhegMEPP2VV:Coupling /Herwig/Shower/AlphaQCD
##########################################################
# DIS matrix elements
##########################################################
# neutral current
create Herwig::MENeutralCurrentDIS MEDISNC
+set MEDISNC:Coupling /Herwig/Shower/AlphaQCD
+set MEDISNC:Contribution 0
# charged current
create Herwig::MEChargedCurrentDIS MEDISCC
+set MEDISCC:Coupling /Herwig/Shower/AlphaQCD
+set MEDISCC:Contribution 0
+
+# neutral current (POWHEG)
+create Herwig::MENeutralCurrentDIS PowhegMEDISNC
+set PowhegMEDISNC:Coupling /Herwig/Shower/AlphaQCD
+set PowhegMEDISNC:Contribution 1
+# charged current (POWHEG)
+create Herwig::MEChargedCurrentDIS PowhegMEDISCC
+set PowhegMEDISCC:Coupling /Herwig/Shower/AlphaQCD
+set PowhegMEDISCC:Contribution 1
##########################################################
# Gamma-Gamma matrix elements
##########################################################
# fermion-antiferimon
create Herwig::MEGammaGamma2ff MEgg2ff HwMEGammaGamma.so
# W+ W-
create Herwig::MEGammaGamma2WW MEgg2WW HwMEGammaGamma.so
##########################################################
# Gamma-Hadron matrix elements
##########################################################
# gamma parton -> 2 jets
create Herwig::MEGammaP2Jets MEGammaP2Jets HwMEGammaHadron.so
##########################################################
# Set up the Subprocesses
#
# For e+e-
##########################################################
create ThePEG::SubProcessHandler SimpleEE
newdef SimpleEE:PartonExtractor /Herwig/Partons/EEExtractor
##########################################################
# For hadron-hadron
##########################################################
create ThePEG::SubProcessHandler SimpleQCD
newdef SimpleQCD:PartonExtractor /Herwig/Partons/QCDExtractor
##########################################################
# For DIS
##########################################################
create ThePEG::SubProcessHandler SimpleDIS
newdef SimpleDIS:PartonExtractor /Herwig/Partons/DISExtractor

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:25 PM (1 d, 5 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806025
Default Alt Text
(248 KB)

Event Timeline