Page MenuHomeHEPForge
Files F7879088

No OneTemporary
Actions

View Options

diff --git a/.hgtags b/.hgtags
new file mode 100644
--- /dev/null
+++ b/.hgtags
@@ -0,0 +1,1 @@
+6c1be2daaa215d44a15d33c701f2e17750c3a942 release 1.1
diff --git a/Amplitudes/AmplitudeBase.cc.in b/Amplitudes/AmplitudeBase.cc.in
--- a/Amplitudes/AmplitudeBase.cc.in
+++ b/Amplitudes/AmplitudeBase.cc.in
@@ -1,814 +1,814 @@
// -*- C++ -*-
//
// AmplitudeBase.cc is a part of HJets
// Copyright (C) 2011-2012
// Ken Arnold, Francisco Campanario, Terrance Figy, Simon Platzer and Malin Sjodahl
//
// HJets 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 AmplitudeBase class.
//
#include "AmplitudeBase.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace HJets;
AmplitudeBase::AmplitudeBase()
: MatchboxAmplitude(), complexMassScheme(true),
mu2Factor(1.0),
theKappaZ(1.0), theKappaW(1.0),
theStableEpsilon(1e-9) {}
AmplitudeBase::~AmplitudeBase() {}
bool AmplitudeBase::canHandle(const PDVector& proc) const {
int n = proc.size();
if ( n != nPartons() + 1 ) {
return false;
}
return HJetsProcessInfo::canHandle(proc,SM());
}
map<cPDVector,vector<AmplitudeInfo> >& AmplitudeBase::amplitudeInfos() {
static map<cPDVector,vector<AmplitudeInfo> > theAmplitudeInfos;
return theAmplitudeInfos;
}
map<cPDVector,map<int, pair<int, double> > >& AmplitudeBase::virtualInfos() {
static map<cPDVector,map<int, pair<int, double> > > theVirtualInfos;
return theVirtualInfos;
}
struct cfless {
inline bool operator() (const vector<size_t>& a,
const vector<size_t>& b) const {
if ( a.size() > b.size() )
return true;
return less<vector<size_t> >()(a,b);
}
};
size_t toColorFull(const set<vector<size_t>,cfless>& cs) {
ostringstream cfids;
cfids << "[";
for ( set<vector<size_t>,cfless>::const_iterator l = cs.begin();
l != cs.end(); ++l ) {
cfids << "{";
for ( vector<size_t>::const_iterator i = l->begin(); i != l->end(); ++i ) {
cfids << *i;
if ( i != l->end()-1 )
cfids << ",";
}
cfids << "}";
}
cfids << "]";
string cfid = cfids.str();
if ( cfid == "[{1,2}{3,4}]" ) return 0;
if ( cfid == "[{1,4}{3,2}]" ) return 1;
if ( cfid == "[{1,5,2}{3,4}]" ) return 0;
if ( cfid == "[{1,5,4}{3,2}]" ) return 1;
if ( cfid == "[{3,5,2}{1,4}]" ) return 2;
if ( cfid == "[{3,5,4}{1,2}]" ) return 3;
if ( cfid == "[{1,5,6,2}{3,4}]" ) return 0;
if ( cfid == "[{1,5,6,4}{3,2}]" ) return 1;
if ( cfid == "[{1,6,5,2}{3,4}]" ) return 2;
if ( cfid == "[{1,6,5,4}{3,2}]" ) return 3;
if ( cfid == "[{3,5,6,2}{1,4}]" ) return 4;
if ( cfid == "[{3,5,6,4}{1,2}]" ) return 5;
if ( cfid == "[{3,6,5,2}{1,4}]" ) return 6;
if ( cfid == "[{3,6,5,4}{1,2}]" ) return 7;
if ( cfid == "[{1,5,2}{3,6,4}]" ) return 8;
if ( cfid == "[{1,5,4}{3,6,2}]" ) return 9;
if ( cfid == "[{1,6,2}{3,5,4}]" ) return 10;
if ( cfid == "[{1,6,4}{3,5,2}]" ) return 11;
if ( cfid == "[{1,2}{3,4}{5,6}]" ) return 0;
if ( cfid == "[{1,2}{3,6}{5,4}]" ) return 1;
if ( cfid == "[{1,4}{3,2}{5,6}]" ) return 2;
if ( cfid == "[{1,4}{3,6}{5,2}]" ) return 3;
if ( cfid == "[{1,6}{3,2}{5,4}]" ) return 4;
if ( cfid == "[{1,6}{3,4}{5,2}]" ) return 5;
assert(false);
return 0;
}
vector<AmplitudeInfo>& AmplitudeBase::amplitudeInfo() {
const cPDVector& proc = mePartonData();
map<cPDVector,vector<AmplitudeInfo> >::iterator i =
amplitudeInfos().find(proc);
if ( i != amplitudeInfos().end() )
return i->second;
#ifdef AMPVERBOSE
// --------------------------------------------------------------------------------
generator()->log() << "--------------------------------------------------------------------------------\n";
generator()->log() << "\n" << name() << "\n";
generator()->log() << "getting amplitudeInfo:\nfor subprocess: ";
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p ) {
generator()->log() << (**p).PDGName() << " ";
}
generator()->log() << "\ncrossed to: ";
size_t ampcount = 0;
for ( cPDVector::const_iterator p = amplitudePartonData().begin();
p != amplitudePartonData().end(); ++p, ++ampcount ) {
generator()->log() << ampcount << " : " << (**p).PDGName()
<< " [" << crossingMap()[ampcount] << "] ";
}
generator()->log() << "\n" << flush;
// generator()->log() << "using crossing sign = "
// << crossingSign() << "\n" << flush;
generator()->log() << "amplitude ids map to colour basis ids as:\n" << flush;
for ( map<size_t,size_t>::const_iterator c = amplitudeToColourMap().begin();
c != amplitudeToColourMap().end(); ++c ) {
generator()->log() << "[" << c->first << " -> " << (c->second + 1)
<< "] ";
}
generator()->log() << "\n" << flush;
// --------------------------------------------------------------------------------
#endif
amplitudeInfos()[proc] =
HJetsProcessInfo::getConfigurations(mePartonData(),crossingMap(),SM(),complexMassScheme,
kappaZ(),kappaW());
i = amplitudeInfos().find(proc);
map<size_t,size_t>& cmap = amplitudeToColourMap();
#ifdef AMPVERBOSE
// --------------------------------------------------------------------------------
generator()->log() << "generated amplitude info:\n\n" << flush;
// --------------------------------------------------------------------------------
#endif
for ( vector<AmplitudeInfo>::iterator k = i->second.begin();
k != i->second.end(); ++k ) {
#ifdef AMPVERBOSE
// --------------------------------------------------------------------------------
k->print(generator()->log());
// --------------------------------------------------------------------------------
#endif
map<set<vector<size_t>,cfless>,double> colourStructures;
if ( !k->qqbarEmitted ) {
vector<size_t> ijLine;
ijLine.push_back(cmap[k->ijLine.first]+1);
if ( k->ijLineEmissions.first > 0 )
ijLine.push_back(cmap[k->ijLineEmissions.first]+1);
if ( k->ijLineEmissions.second > 0 )
ijLine.push_back(cmap[k->ijLineEmissions.second]+1);
ijLine.push_back(cmap[k->ijLine.second]+1);
vector<size_t> klLine;
klLine.push_back(cmap[k->klLine.first]+1);
if ( k->klLineEmissions.first > 0 )
klLine.push_back(cmap[k->klLineEmissions.first]+1);
if ( k->klLineEmissions.second > 0 )
klLine.push_back(cmap[k->klLineEmissions.second]+1);
klLine.push_back(cmap[k->klLine.second]+1);
set<vector<size_t>,cfless> cStructure;
cStructure.insert(ijLine);
cStructure.insert(klLine);
colourStructures.insert(make_pair(cStructure,1.));
} else {
if ( k->ijLineEmissions.first > 0 ) {
vector<size_t> klLine;
klLine.push_back(cmap[k->klLine.first]+1);
klLine.push_back(cmap[k->klLine.second]+1);
vector<size_t> ijLineLeading1;
ijLineLeading1.push_back(cmap[k->ijLine.first]+1);
ijLineLeading1.push_back(cmap[k->ijLineEmissions.second]+1);
vector<size_t> ijLineLeading2;
ijLineLeading2.push_back(cmap[k->ijLineEmissions.first]+1);
ijLineLeading2.push_back(cmap[k->ijLine.second]+1);
vector<size_t> ijLineSubleading1;
ijLineSubleading1.push_back(cmap[k->ijLine.first]+1);
ijLineSubleading1.push_back(cmap[k->ijLine.second]+1);
vector<size_t> ijLineSubleading2;
ijLineSubleading2.push_back(cmap[k->ijLineEmissions.first]+1);
ijLineSubleading2.push_back(cmap[k->ijLineEmissions.second]+1);
set<vector<size_t>,cfless> cStructureLeading;
cStructureLeading.insert(klLine);
cStructureLeading.insert(ijLineLeading1);
cStructureLeading.insert(ijLineLeading2);
set<vector<size_t>,cfless> cStructureSubleading;
cStructureSubleading.insert(klLine);
cStructureSubleading.insert(ijLineSubleading1);
cStructureSubleading.insert(ijLineSubleading2);
colourStructures.insert(make_pair(cStructureLeading,0.5));
colourStructures.insert(make_pair(cStructureSubleading,-0.5/3.));
}
if ( k->klLineEmissions.first > 0 ) {
vector<size_t> ijLine;
ijLine.push_back(cmap[k->ijLine.first]+1);
ijLine.push_back(cmap[k->ijLine.second]+1);
vector<size_t> klLineLeading1;
klLineLeading1.push_back(cmap[k->klLine.first]+1);
klLineLeading1.push_back(cmap[k->klLineEmissions.second]+1);
vector<size_t> klLineLeading2;
klLineLeading2.push_back(cmap[k->klLineEmissions.first]+1);
klLineLeading2.push_back(cmap[k->klLine.second]+1);
vector<size_t> klLineSubleading1;
klLineSubleading1.push_back(cmap[k->klLine.first]+1);
klLineSubleading1.push_back(cmap[k->klLine.second]+1);
vector<size_t> klLineSubleading2;
klLineSubleading2.push_back(cmap[k->klLineEmissions.first]+1);
klLineSubleading2.push_back(cmap[k->klLineEmissions.second]+1);
set<vector<size_t>,cfless> cStructureLeading;
cStructureLeading.insert(ijLine);
cStructureLeading.insert(klLineLeading1);
cStructureLeading.insert(klLineLeading2);
set<vector<size_t>,cfless> cStructureSubleading;
cStructureSubleading.insert(ijLine);
cStructureSubleading.insert(klLineSubleading1);
cStructureSubleading.insert(klLineSubleading2);
colourStructures.insert(make_pair(cStructureLeading,0.5));
colourStructures.insert(make_pair(cStructureSubleading,-0.5/3.));
}
}
#ifdef AMPVERBOSE
// --------------------------------------------------------------------------------
generator()->log() << "with colour structures\n";
for ( map<set<vector<size_t>,cfless>,double>::const_iterator c = colourStructures.begin();
c != colourStructures.end(); ++c ) {
generator()->log() << "[";
for ( set<vector<size_t>,cfless>::const_iterator l = c->first.begin();
l != c->first.end(); ++l ) {
generator()->log() << "{";
for ( vector<size_t>::const_iterator i = l->begin(); i != l->end(); ++i ) {
generator()->log() << *i;
if ( i != l->end()-1 )
generator()->log() << ",";
}
generator()->log() << "}";
}
generator()->log() << "]";
generator()->log() << " x " << c->second << "\n";
}
generator()->log() << "\n" << flush;
// --------------------------------------------------------------------------------
#endif
for ( map<set<vector<size_t>,cfless>,double>::const_iterator c = colourStructures.begin();
c != colourStructures.end(); ++c ) {
k->colourTensors[toColorFull(c->first)] = c->second;
}
#ifdef AMPVERBOSE
// --------------------------------------------------------------------------------
generator()->log() << "mapped to colourfull indices as\n";
for ( map<size_t,double>::const_iterator c = k->colourTensors.begin();
c != k->colourTensors.end(); ++c ) {
generator()->log() << c->first << " x " << c->second
<< " / ";
}
generator()->log() << "\n" << flush;
// --------------------------------------------------------------------------------
#endif
}
#ifdef AMPVERBOSE
// --------------------------------------------------------------------------------
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
// --------------------------------------------------------------------------------
#endif
return i->second;
}
Complex AmplitudeBase::bosonFactor(const AmplitudeInfo& c) const {
Complex ijPropagator(-sqr(c.bosonMass/amplitudeScale()),
(c.bosonMass*c.bosonWidth)/sqr(amplitudeScale()));
Complex klPropagator = ijPropagator;
if ( c.ijLineEmissions.first < 0 && c.ijLineEmissions.second < 0 ) {
ijPropagator +=
invariant(c.ijLine.first,c.ijLine.second);
} else if ( c.ijLineEmissions.first > 0 && c.ijLineEmissions.second < 0 ) {
ijPropagator +=
invariant(c.ijLine.first,c.ijLine.second) +
invariant(c.ijLine.first,c.ijLineEmissions.first) +
invariant(c.ijLine.second,c.ijLineEmissions.first);
} else if ( c.ijLineEmissions.first > 0 && c.ijLineEmissions.second > 0 ) {
ijPropagator +=
invariant(c.ijLine.first,c.ijLine.second) +
invariant(c.ijLine.first,c.ijLineEmissions.first) +
invariant(c.ijLine.second,c.ijLineEmissions.first) +
invariant(c.ijLine.first,c.ijLineEmissions.second) +
invariant(c.ijLine.second,c.ijLineEmissions.second) +
invariant(c.ijLineEmissions.first,c.ijLineEmissions.second);
}
if ( c.klLineEmissions.first < 0 && c.klLineEmissions.second < 0 ) {
klPropagator +=
invariant(c.klLine.first,c.klLine.second);
} else if ( c.klLineEmissions.first > 0 && c.klLineEmissions.second < 0 ) {
klPropagator +=
invariant(c.klLine.first,c.klLine.second) +
invariant(c.klLine.first,c.klLineEmissions.first) +
invariant(c.klLine.second,c.klLineEmissions.first);
} else if ( c.klLineEmissions.first > 0 && c.klLineEmissions.second > 0 ) {
klPropagator +=
invariant(c.klLine.first,c.klLine.second) +
invariant(c.klLine.first,c.klLineEmissions.first) +
invariant(c.klLine.second,c.klLineEmissions.first) +
invariant(c.klLine.first,c.klLineEmissions.second) +
invariant(c.klLine.second,c.klLineEmissions.second) +
invariant(c.klLineEmissions.first,c.klLineEmissions.second);
}
if (c.boson == getParticleData(ParticleID::Z0)) {
return (c.higgsCoupling/amplitudeScale())/(ijPropagator*klPropagator);
}
return (c.higgsCoupling/amplitudeScale())/(ijPropagator*klPropagator);
}
inline bool isQuark(cPDPtr p) {
if (p->id() > 0 && p->id() < 7) return true;
return false;
}
inline bool isAntiQuark(cPDPtr p) {
if (p->id() < 0 && p->id() > -7) return true;
return false;
}
inline bool isGluon(cPDPtr p) {
if (p->id() == 21 ) return true;
return false;
}
inline bool isHiggs(cPDPtr p) {
if (p->id() == 25 ) return true;
return false;
}
inline int sign(int a) {
return a < 0 ? -1 : 1;
}
const map<int, pair<int, double> >& AmplitudeBase::virtualInfo() const {
const cPDVector& proc = mePartonData();
map<cPDVector,map<int, pair<int, double> > >::const_iterator i =
virtualInfos().find(proc);
if ( i != virtualInfos().end() )
return i->second;
vector<int> in;
vector<int> out;
int gluon = 0;
int higgs = 0;
// find incoming and outgoing particles in the master topology
if (isQuark(proc[0])) in.push_back(1);
else if (isAntiQuark(proc[0])) out.push_back(-1);
else if (isGluon(proc[0])) gluon = -1;
if (isQuark(proc[1])) in.push_back(2);
else if (isAntiQuark(proc[1])) out.push_back(-2);
else if (isGluon(proc[1])) gluon = -2;
for (int i = 2; i != proc.size(); i++) {
if (isQuark(proc[i])) out.push_back(i+1);
if (isAntiQuark(proc[i])) in.push_back(-i-1);
if (isGluon(proc[i])) gluon = i+1;
if (isHiggs(proc[i])) higgs = i+1;
}
assert(in.size() == 2);
assert(out.size() == 2);
assert(higgs != 0);
//construct a map for the indices
map<int, pair<int, double> > XMap;
for (int i = 0; i < 2; i++) {
XMap[2*i+1] = make_pair (abs(in[i])-1, sign(in[i]));
}
for (int i = 0; i < 2; i++) {
XMap[2*i+2] = make_pair(abs(out[i])-1, sign(out[i]));
}
if ( gluon != 0 ){
XMap[5] = make_pair(abs(gluon)-1, sign(gluon));
XMap[6] = make_pair(higgs-1, 1.0);
}
else
XMap[5] = make_pair(higgs-1, 1.0);
//swap entries 2 and 4 if they cannot be connected to the incoming partons
//by a t-channel diagram
// int one = 1; int two = 2; int three = 3; int four = 4;
// if (proc[XMap.find(one)->second.first]->id() % 2 * XMap.find(one)->second.second
// != proc[XMap.find(two)->second.first]->id() % 2 * XMap.find(two)->second.second){
// // || (abs(proc[XMap.find(one)->second.first]->id())-1) / 2
// // != (abs(proc[XMap.find(two)->second.first]->id())-1) / 2) {
// assert(proc[XMap.find(two)->second.first]->id() % 2 * XMap.find(two)->second.second
// == proc[XMap.find(three)->second.first]->id() % 2 * XMap.find(three)->second.second);
// assert(proc[XMap.find(one)->second.first]->id() % 2 * XMap.find(one)->second.second
// == proc[XMap.find(four)->second.first]->id() % 2 * XMap.find(four)->second.second);
// // assert((abs(proc[XMap.find(two)->second.first]->id())-1) / 2
// // == (abs(proc[XMap.find(three)->second.first]->id())-1) / 2 );
// // assert((abs(proc[XMap.find(one)->second.first]->id())-1) / 2
// // == (abs(proc[XMap.find(four)->second.first]->id())-1) / 2 );
// pair<int,double> buffer = XMap.find(four)->second;
// XMap[four] = XMap.find(two)->second;
// XMap[two] = buffer;
// }
virtualInfos()[proc] = XMap;
i = virtualInfos().find(proc);
return i->second;
}
bool AmplitudeBase::topologyOneIsNC() const{
const cPDVector& proc = mePartonData();
map<int, pair<int, double> > XMap = virtualInfos().find(proc)->second;
if (proc[XMap.find(1)->second.first]->iCharge() * XMap.find(1)->second.second
- proc[XMap.find(2)->second.first]->iCharge() * XMap.find(2)->second.second == ZERO
&& (abs(proc[XMap.find(1)->second.first]->id())-1)/2
== (abs(proc[XMap.find(2)->second.first]->id())-1)/2){
assert(proc[XMap.find(3)->second.first]->iCharge() * XMap.find(3)->second.second
- proc[XMap.find(4)->second.first]->iCharge() * XMap.find(4)->second.second == ZERO);
return true;
}
return false;
}
bool AmplitudeBase::topologyOneIsCC() const{
const cPDVector& proc = mePartonData();
map<int, pair<int, double> > XMap = virtualInfos().find(proc)->second;
if (abs(proc[XMap.find(1)->second.first]->iCharge() * XMap.find(1)->second.second
- proc[XMap.find(2)->second.first]->iCharge() * XMap.find(2)->second.second) == 3
&& (abs(proc[XMap.find(1)->second.first]->id())-1)/2
== (abs(proc[XMap.find(2)->second.first]->id())-1)/2
&& (abs(proc[XMap.find(3)->second.first]->id())-1)/2
== (abs(proc[XMap.find(4)->second.first]->id())-1)/2){
assert(abs(proc[XMap.find(3)->second.first]->iCharge() * XMap.find(3)->second.second
- proc[XMap.find(4)->second.first]->iCharge() * XMap.find(4)->second.second) == 3);
return true;
}
return false;
}
bool AmplitudeBase::topologyTwoIsNC() const{
const cPDVector& proc = mePartonData();
map<int, pair<int, double> > XMap = virtualInfos().find(proc)->second;
if (proc[XMap.find(1)->second.first]->iCharge() * XMap.find(1)->second.second
- proc[XMap.find(4)->second.first]->iCharge() * XMap.find(4)->second.second == ZERO
&& (abs(proc[XMap.find(1)->second.first]->id())-1)/2
== (abs(proc[XMap.find(4)->second.first]->id())-1)/2){
assert(proc[XMap.find(3)->second.first]->iCharge() * XMap.find(3)->second.second
- proc[XMap.find(2)->second.first]->iCharge() * XMap.find(2)->second.second == ZERO);
return true;
}
return false;
}
bool AmplitudeBase::topologyTwoIsCC() const{
const cPDVector& proc = mePartonData();
map<int, pair<int, double> > XMap = virtualInfos().find(proc)->second;
if (abs(proc[XMap.find(1)->second.first]->iCharge() * XMap.find(1)->second.second
- proc[XMap.find(4)->second.first]->iCharge() * XMap.find(4)->second.second) == 3
&& (abs(proc[XMap.find(1)->second.first]->id())-1)/2
== (abs(proc[XMap.find(4)->second.first]->id())-1)/2
&& (abs(proc[XMap.find(3)->second.first]->id())-1)/2
== (abs(proc[XMap.find(2)->second.first]->id())-1)/2){
assert(abs(proc[XMap.find(3)->second.first]->iCharge() * XMap.find(3)->second.second
- proc[XMap.find(2)->second.first]->iCharge() * XMap.find(2)->second.second) == 3);
return true;
}
return false;
}
void AmplitudeBase::getCouplings(double* Recl1, double* Recr1, double* Recl3, double* Recr3,
double* Imcl1, double* Imcr1, double* Imcl3, double* Imcr3,
double* Reghvv, double* Imghvv,
double* mass, double* width, int* isVaW) const{
const cPDVector& proc = mePartonData();
map<int, pair<int, double> > XMap = virtualInfos().find(proc)->second;
for (int i = 0; i < 2; i++) {
Recl1[i]=0.0;
Recr1[i]=0.0;
Recl3[i]=0.0;
Recr3[i]=0.0;
Reghvv[i]=0.0;
Imcl1[i]=0.0;
Imcr1[i]=0.0;
Imcl3[i]=0.0;
Imcr3[i]=0.0;
Imghvv[i]=0.0;
mass[i]=10.0;
width[i]=10.0;
isVaW[i]=0;
}
// fix width
//double gem = sqrt(4.*Constants::pi*SM().alphaEMMZ());
//double s2tw = SM().sin2ThetaW();
//double stw = sqrt(s2tw);
//double ctw = sqrt(1.-s2tw);
- Energy MW = getParticleData(ParticleID::Wplus)->mass();
- Energy MZ = getParticleData(ParticleID::Z0)->mass();
+ Energy MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
+ Energy MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
- Energy WWidth = getParticleData(ParticleID::Wplus)->width();
- Energy ZWidth = getParticleData(ParticleID::Z0)->width();
+ Energy WWidth = getParticleData(ParticleID::Wplus)->hardProcessWidth();
+ Energy ZWidth = getParticleData(ParticleID::Z0)->hardProcessWidth();
//cout << "Wwidth="<< WWidth/GeV << "\n";
//cout << "Zwidth="<< ZWidth/GeV << "\n";
complex<double> alphaQED = SM().alphaEMMZ();
complex<double> c2tw = complex<double>(1.,0.) - SM().sin2ThetaW();
complex<double> s2tw = SM().sin2ThetaW();
complex<double> MW2c = double(sqr(MW/GeV));
complex<double> MWc = sqrt(MW2c);
complex<double> MZ2c = double(sqr(MZ/GeV));
complex<double> MZc = sqrt(MZ2c);
if( complexMassScheme ){
MW2c = complex<double>(sqr(MW/GeV),-MW*WWidth/GeV2);
MWc = sqrt(MW2c);
MZ2c = complex<double>(sqr(MZ/GeV),-MZ*ZWidth/GeV2);
MZc = sqrt(MZ2c);
c2tw = MW2c/MZ2c;
s2tw = complex<double>(1.,0.) - c2tw;
alphaQED = sqrt(2.)*SM().fermiConstant()*MW2c*s2tw/Constants::pi*GeV2;
}
/*cout << "s2tw= "<< s2tw << "\n";
cout << "MZc= " << MZc << "\n";
cout << "MWc= " << MWc << "\n";*/
complex<double> stw = sqrt(s2tw);
complex<double> ctw = sqrt(c2tw);
complex <double> gem = sqrt(4.*Constants::pi*alphaQED);
//get couplings to Z0
double Q = proc[XMap.find(1)->second.first]->iCharge()/3. * XMap.find(1)->second.second;
double I3 = proc[XMap.find(1)->second.first]->id() % 2 == 0 ? 0.5 : -0.5;
complex<double> CUpperR = -gem*Q*stw/ctw;
complex<double> CUpperL = gem*(I3-s2tw*Q)/(stw*ctw);
Q = proc[XMap.find(3)->second.first]->iCharge()/3. * XMap.find(3)->second.second;
I3 = proc[XMap.find(3)->second.first]->id() % 2 == 0 ? 0.5 : -0.5;
complex<double> CLowerR = -gem*Q*stw/ctw;
complex<double> CLowerL = gem*(I3-s2tw*Q)/(stw*ctw);
//find NC topology and fill in couplings
complex<double> ghzz = kappaZ()*gem*MZc/(stw*ctw);
if (topologyOneIsNC()) {
Recl1[0]=CUpperL.real();
Recr1[0]=CUpperR.real();
Recl3[0]=CLowerL.real();
Recr3[0]=CLowerR.real();
Reghvv[0]=ghzz.real();
if (complexMassScheme ){
Imcl1[0]=CUpperL.imag();
Imcr1[0]=CUpperR.imag();
Imcl3[0]=CLowerL.imag();
Imcr3[0]=CLowerR.imag();
Imghvv[0]=ghzz.imag();
// throw "go here";
}
mass[0]=MZ/GeV;
width[0]=ZWidth/GeV;
isVaW[0]=0;
}
if (topologyTwoIsNC()) {
Recl1[1]=CUpperL.real();
Recr1[1]=CUpperR.real();
Recl3[1]=CLowerL.real();
Recr3[1]=CLowerR.real();
Reghvv[1]=ghzz.real();
if (complexMassScheme){
Imcl1[1]=CUpperL.imag();
Imcr1[1]=CUpperR.imag();
Imcl3[1]=CLowerL.imag();
Imcr3[1]=CLowerR.imag();
Imghvv[1]=ghzz.imag();
}
mass[1]=MZ/GeV;
width[1]=ZWidth/GeV;
isVaW[1]=0;
}
//get couplings to W
CUpperL = gem/(sqrt(2.)*stw);
CLowerL = gem/(sqrt(2.)*stw);
complex<double> ghww = kappaW()*gem*MWc/stw;
assert(!(topologyOneIsCC() && topologyTwoIsCC()));
//find CC topology and fill in couplings
if (topologyOneIsCC()) {
Recl1[0]=CUpperL.real();
Recl3[0]=CLowerL.real();
Reghvv[0]=ghww.real();
if (complexMassScheme){
Imcl1[0]=CUpperL.imag();
Imcl3[0]=CLowerL.imag();
Imghvv[0]=ghww.imag();
}
mass[0]=MW/GeV;
width[0]=WWidth/GeV;
isVaW[0]=1;
}
else if (topologyTwoIsCC()) {
Recl1[1]=CUpperL.real();
Recl3[1]=CLowerL.real();
Reghvv[1]=ghww.real();
if (complexMassScheme){
Imcl1[1]=CUpperL.imag();
Imcl3[1]=CLowerL.imag();
Imghvv[1]=ghww.imag();
}
mass[1]=MW/GeV;
width[1]=WWidth/GeV;
isVaW[1]=1;
}
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void AmplitudeBase::persistentOutput(PersistentOStream & os) const {
os << complexMassScheme << mu2Factor
<< theKappaZ << theKappaW << theStableEpsilon;
}
void AmplitudeBase::persistentInput(PersistentIStream & is, int) {
is >> complexMassScheme >> mu2Factor
>> theKappaZ >> theKappaW >> theStableEpsilon;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<AmplitudeBase,Herwig::MatchboxAmplitude>
describeHJetsAmplitudeBase("HJets::AmplitudeBase", "HwMatchboxBuiltin.so HJets.so");
void AmplitudeBase::Init() {
static ClassDocumentation<AmplitudeBase> documentation
("Helicity amplitudes for 0 -> h + jets");
static Switch<AmplitudeBase,bool> interfaceComplexMassScheme
("ComplexMassScheme",
"Switch on or off the complex mass scheme.",
&AmplitudeBase::complexMassScheme, true, false, false);
static SwitchOption interfaceComplexMassSchemeOn
(interfaceComplexMassScheme,
"On",
"Switch on the complex mass scheme.",
true);
static SwitchOption interfaceComplexMassSchemeOff
(interfaceComplexMassScheme,
"Off",
"Switch off the complex mass scheme.",
false);
static Parameter<AmplitudeBase,double> interfaceMu2Factor
("Mu2Factor",
"The t'Hooft scale factor to be changed for validation.",
&AmplitudeBase::mu2Factor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<AmplitudeBase,double> interfaceKappaZ
("KappaZ",
"The factor to rescale the HZZ coupling.",
&AmplitudeBase::theKappaZ, 1.0, 0, 0,
false, false, Interface::nolimits);
static Parameter<AmplitudeBase,double> interfaceKappaW
("KappaW",
"The factor to rescale the HWW coupling.",
&AmplitudeBase::theKappaW, 1.0, 0, 0,
false, false, Interface::nolimits);
static Parameter<AmplitudeBase,double> interfaceStableEpsilon
("StableEpsilon",
"The threshold for unstable point classification.",
&AmplitudeBase::theStableEpsilon, 1e-9, 0.0, 0,
false, false, Interface::lowerlim);
}
namespace HJets {
struct Banner {
Banner() {
cout << "################################################################################\n"
<< "# #\n"
- << "# HJets++ 1.0 #\n"
+ << "# HJets++ 1.1 #\n"
<< "# ========================================================================== #\n"
<< "# A Matchbox plugin for electroweak Higgs plus #\n"
<< "# two and three jet production at NLO QCD #\n"
<< "# #\n"
<< "# #\n"
<< "# (C) 2012-2014 The HJets++ collaboration #\n"
<< "# ------------------------------------------------------------ #\n"
<< "# Francisco Campanario <francisco.campanario@ific.uv.es> #\n"
<< "# Terrance M. Figy <terrance.figy@hep.manchester.ac.uk> #\n"
<< "# Simon Platzer <simon.plaetzer@desy.de> #\n"
<< "# Malin Sjodahl <malin.sjodahl@thep.lu.se> #\n"
<< "# ------------------------------------------------------------ #\n"
<< "# #\n"
<< "# #\n"
<< "# ========================================================================== #\n"
<< "# #\n"
<< "# HJets++ is licenced under version 2 of the GPL, see COPYING for details. #\n"
<< "# Please respect the MCnet academic guidelines, see GUIDELINES for details. #\n"
<< "# #\n"
<< "# When using HJets++ please cite: #\n"
<< "# #\n"
<< "# F. Campanario, T.M. Figy, S. Platzer and M. Sjodahl: #\n"
<< "# 'Electroweak Higgs plus Three Jet Production at NLO QCD', #\n"
<< "# Phys.Rev.Lett. 111 (2013) 211802, arXiv:1308.2932 [hep-ph] #\n"
<< "# #\n"
<< "################################################################################\n"
<< flush;
Repository::appendReadDir("@prefix@/share/HJets");
}
};
}
HJets::Banner hjetsbanner;
diff --git a/Amplitudes/AmplitudeBase.h b/Amplitudes/AmplitudeBase.h
--- a/Amplitudes/AmplitudeBase.h
+++ b/Amplitudes/AmplitudeBase.h
@@ -1,239 +1,239 @@
// -*- C++ -*-
//
// AmplitudeBase.h is a part of HJets
// Copyright (C) 2011-2012
// Ken Arnold, Francisco Campanario, Terrance Figy, Simon Platzer and Malin Sjodahl
//
// HJets is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HJets_AmplitudeBase_H
#define HJets_AmplitudeBase_H
//
// This is the declaration of the AmplitudeBase class.
//
-#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
-#include "Herwig++/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxCurrents.h"
+#include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
+#include "Herwig/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxCurrents.h"
#include "HJetsProcessInfo.h"
namespace HJets {
using namespace ThePEG;
using namespace Herwig;
//#define AMPVERBOSE
/**
* \brief Tree level helicity amplitudes for 0 -> h + jets
*/
class AmplitudeBase:
public MatchboxAmplitude, public MatchboxCurrents {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
AmplitudeBase();
/**
* The destructor.
*/
virtual ~AmplitudeBase();
//@}
public:
/**
* Return true, if this amplitude can handle the given process.
*/
virtual bool canHandle(const PDVector&) const;
/**
* Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGem() const { return 3; }
/**
* Return the (tree-level) order in \f$g_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGs() const { return nPartons()-4; }
/**
* Flush all cashes.
*/
virtual void flushCaches() {
MatchboxCurrents::reset();
MatchboxAmplitude::flushCaches();
}
/**
* Return the number of partons attached to this amplitude.
*/
virtual int nPartons() const = 0;
/**
* Return true, if one-loop contributions will be evaluated at amplitude level.
*/
virtual bool oneLoopAmplitudes() const { return false; }
/**
* Return true, if one loop corrections are given in the convenctions
* of everything expanded.
*/
virtual bool isExpanded() const { return true; }
//virtual bool isDR() const{return false;}
/**
* Return the virtual amplitude info
*/
const map<int, pair<int, double> >& virtualInfo() const;
/**
* The virtuals info maps.
*/
static map<cPDVector,map<int, pair<int, double> > >& virtualInfos();
/**
* Return the value of the dimensional regularization
* parameter. Note that renormalization scale dependence is fully
* restored in DipoleIOperator.
*/
virtual Energy2 mu2() const { return mu2Factor*lastSHat(); }
protected:
/**
* Return the amplitude configurations to be considered.
*/
vector<AmplitudeInfo>& amplitudeInfo();
/**
* Get the boson propagators and Higgs coupling factor
*/
Complex bosonFactor(const AmplitudeInfo&) const;
/**
* Check if t-channel diagram is a neutral current diagram.
*/
bool topologyOneIsNC() const;
/**
* Check if u-channel diagram is a neutral current diagram.
*/
bool topologyTwoIsNC() const;
/**
* Check if t-channel diagram is a neutral current diagram.
*/
bool topologyOneIsCC() const;
/**
* Check if u-channel diagram is a neutral current diagram.
*/
bool topologyTwoIsCC() const;
/**
* Work out the couplings
*/
void getCouplings(double*, double*, double*, double*,
double*, double*, double*, double*, double*, double*, double*, double*, int*) 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();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/**
* Return the factor to rescale the HZZ coupling
*/
double kappaZ() const { return theKappaZ; }
/**
* Return the factor to rescale the HWW coupling
*/
double kappaW() const { return theKappaW; }
protected:
/**
* Return the threshold for classifying points as unstable
*/
double stableEpsilon() const { return theStableEpsilon; }
private:
/**
* The amplitude info map.
*/
static map<cPDVector,vector<AmplitudeInfo> >& amplitudeInfos();
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
AmplitudeBase & operator=(const AmplitudeBase &);
/**
* Wether or not the complex mass scheme should be used
*/
bool complexMassScheme;
/**
* A rescaling factor for the t'Hooft mass
*/
double mu2Factor;
/**
* A factor to rescale the HZZ coupling
*/
double theKappaZ;
/**
* A factor to rescale the HWW coupling
*/
double theKappaW;
/**
* The threshold for classifying points as unstable
*/
double theStableEpsilon;
};
}
#endif /* HJets_AmplitudeBase_H */
diff --git a/Amplitudes/Amplitudehqqbarkkbar.cc b/Amplitudes/Amplitudehqqbarkkbar.cc
--- a/Amplitudes/Amplitudehqqbarkkbar.cc
+++ b/Amplitudes/Amplitudehqqbarkkbar.cc
@@ -1,486 +1,489 @@
// -*- C++ -*-
//
// Amplitudehqqbarkkbar.cc is a part of HJets
// Copyright (C) 2011-2012
// Ken Arnold, Francisco Campanario, Terrance Figy, Simon Platzer and Malin Sjodahl
//
// HJets 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 Amplitudehqqbarkkbar class.
//
#include "Amplitudehqqbarkkbar.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
extern "C"
{
void qqhqqvbf_(double PBAR1[4], double PBAR2[4], double PBAR3[4],
double PBAR4[4], double PBAR5[4],
int SIGN[5], double ReCL1[2], double ReCR1[2],
double ReCL3[2], double ReCR3[2],
double ImCL1[2], double ImCR1[2],
double ImCL3[2], double ImCR3[2],
double Mass[2], double Width[2], int IsVaW[2],
int& divMax,
double ReGHVV[2], double ImGHVV[2], double& scale,
double& born, double virt[3]);
}
inline void F77_qqhqqvbf(double PBAR1[4], double PBAR2[4], double PBAR3[4],
double PBAR4[4], double PBAR5[4],
int SIGN[5], double ReCL1[2], double ReCR1[2],
double ReCL3[2], double ReCR3[2],
double ImCL1[2], double ImCR1[2],
double ImCL3[2], double ImCR3[2],
double Mass[2], double Width[2], int IsVaW[2], int& divMax,
double ReGHVV[2], double ImGHVV[2], double& scale,
double& born, double virt[3]){
qqhqqvbf_(PBAR1,PBAR2,PBAR3,PBAR4,PBAR5,SIGN,ReCL1,ReCR1,ReCL3,ReCR3,
ImCL1,ImCR1,ImCL3,ImCR3,Mass,Width,IsVaW,divMax,ReGHVV,ImGHVV,scale,born,virt);
return;
}
using namespace HJets;
Amplitudehqqbarkkbar::Amplitudehqqbarkkbar() {}
Amplitudehqqbarkkbar::~Amplitudehqqbarkkbar() {}
IBPtr Amplitudehqqbarkkbar::clone() const {
return new_ptr(*this);
}
IBPtr Amplitudehqqbarkkbar::fullclone() const {
return new_ptr(*this);
}
void Amplitudehqqbarkkbar::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
// inform the amplitude cache how many legs we are handling
nPoints(5);
// everything is expressed in dimensionless quantities, referring to
// sqrt(shat) of the collission as the mass unit
amplitudeScale(sqrt(lastSHat()));
// the higgs momentum
Lorentz5Momentum ph = amplitudeMomentum(0);
momentum(0,ph,true,ph.mass());
// the parton momenta
for ( size_t k = 1; k < 5; ++k )
momentum(k,amplitudeMomentum(k),true);
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex Amplitudehqqbarkkbar::evaluate(size_t colourTensor,
const vector<int>& helicities,
Complex& largeN) {
// get configurations we need to consider
const vector<AmplitudeInfo>& confs = amplitudeInfo();
// sum of all contributions
Complex result = 0;
largeN = 0.;
for ( vector<AmplitudeInfo>::const_iterator c = confs.begin();
c != confs.end(); ++c ) {
map<size_t,double>::const_iterator cx = c->colourTensors.find(colourTensor);
if ( cx == c->colourTensors.end() )
continue;
Complex diagram = 0.;
Complex ijL = c->ijLeftCoupling;
Complex ijR = c->ijRightCoupling;
Complex klL = c->klLeftCoupling;
Complex klR = c->klRightCoupling;
int i = c->ijLine.first;
int j = c->ijLine.second;
int k = c->klLine.first;
int l = c->klLine.second;
if ( c->ijLineEmissions.first < 0 && c->ijLineEmissions.second < 0 &&
c->klLineEmissions.first < 0 && c->klLineEmissions.second < 0 ) {
if ( ijL != 0. && klL != 0. ) {
diagram +=
ijL*klL*qqbarLeftCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbarLeftCurrent(k,helicities[k],
l,helicities[l]));
}
if ( ijL != 0. && klR != 0. ) {
diagram +=
ijL*klR*qqbarLeftCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbarRightCurrent(k,helicities[k],
l,helicities[l]));
}
if ( ijR != 0. && klL != 0. ) {
diagram +=
ijR*klL*qqbarRightCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbarLeftCurrent(k,helicities[k],
l,helicities[l]));
}
if ( ijR != 0. && klR != 0. ) {
diagram +=
ijR*klR*qqbarRightCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbarRightCurrent(k,helicities[k],
l,helicities[l]));
}
} else assert(false);
diagram *= bosonFactor(*c)*cx->second * c->fermionSign;
result += diagram;
if ( cx->second > 0. )
largeN += diagram;
}
return result;
}
inline void convert(const Lorentz5Momentum& p, double* res) {
res[0] = p.t()/GeV;
res[1] = p.x()/GeV;
res[2] = p.y()/GeV;
res[3] = p.z()/GeV;
}
double Amplitudehqqbarkkbar::oneLoopInterference() const {
+ if ( !calculateOneLoopInterference() )
+ return lastOneLoopInterference();
const map<int, pair<int, double> >& transmap = virtualInfo();
double pbar1[4];
double pbar2[4];
double pbar3[4];
double pbar4[4];
double pbar5[4];
convert(meMomenta()[transmap.find(1)->second.first],pbar1);
convert(meMomenta()[transmap.find(2)->second.first],pbar2);
convert(meMomenta()[transmap.find(3)->second.first],pbar3);
convert(meMomenta()[transmap.find(4)->second.first],pbar4);
convert(meMomenta()[transmap.find(5)->second.first],pbar5);
int sign[5];
for ( size_t k = 0; k < 5; ++k )
sign[k] = transmap.find(k+1)->second.second;
double Recl1[2];
double Recr1[2];
double Recl3[2];
double Recr3[2];
double Reghvv[2];
double Imcl1[2];
double Imcr1[2];
double Imcl3[2];
double Imcr3[2];
double Imghvv[2];
double mass[2];
double width[2];
int isVaW[2];
getCouplings(Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,Reghvv,Imghvv,mass,width,isVaW);
//int NF = nLight();
int divMax = 2;
double born;
double virt[3];
double scale = mu2()/GeV2;
F77_qqhqqvbf (pbar1,pbar2,pbar3,pbar4,pbar5,sign,Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,
mass,width,isVaW,divMax,Reghvv,Imghvv,scale,born,virt);
double res = (SM().alphaS()/(2.*Constants::pi))*virt[0]*(lastSHat()/GeV2);
#ifdef AMPVERBOSE
const cPDVector& proc = mePartonData();
const cPDVector& proc1 = amplitudePartonData();
map<int, pair<int, double> > XMap = virtualInfos().find(proc)->second;
generator()->log() << "c2="<< virt[2]/born << "\n";
generator()->log() << "c1="<< virt[1]/born << "\n";
generator()->log() << "c0="<< virt[0]/born << "\n";
born*=lastSHat()/GeV2;
generator()->log() << "--------------------- \n";
- generator()->log() << "rb=" << born/me2()/LastMatchboxXCombInfo::crossingSign() << " Base(TF) Diagram:"
+ generator()->log() << "rb=" << born/me2() << " Base(TF) Diagram:"
<< XMap.find(1)->second.first << (XMap.find(1)->second.second > 0 ? "+" : "-") << proc[XMap.find(1)->second.first]->PDGName() << " "
<< XMap.find(3)->second.first << (XMap.find(3)->second.second > 0 ? "+" : "-") << proc[XMap.find(3)->second.first]->PDGName() << " -> "
<< XMap.find(2)->second.first << (XMap.find(2)->second.second > 0 ? "+" : "-") << proc[XMap.find(2)->second.first]->PDGName() << " "
<< XMap.find(4)->second.first << (XMap.find(4)->second.second > 0 ? "+" : "-") << proc[XMap.find(4)->second.first]->PDGName() << " "
<< XMap.find(5)->second.first << (XMap.find(5)->second.second > 0 ? "+" : "-") << proc[XMap.find(5)->second.first]->PDGName() << "\n";
- generator()->log() << "rb=" << born/me2()/LastMatchboxXCombInfo::crossingSign() << " Base(SP) Diagram:"
+ generator()->log() << "rb=" << born/me2() << " Base(SP) Diagram:"
<< crossingMap()[0] << " " << proc1[0]->PDGName() << " "
<< crossingMap()[1] << " " << proc1[1]->PDGName() << " "
<< crossingMap()[2] << " " << proc1[2]->PDGName() << " "
<< crossingMap()[3] << " " << proc1[3]->PDGName() << " "
<< crossingMap()[4] << " " << proc1[4]->PDGName() << "\n";
generator()->log() << "OneisNC= " << topologyOneIsNC()
<< " OneisCC= " << topologyOneIsCC()
<< " TwoisNC= " << topologyTwoIsNC()
<< " TwoisCC= " << topologyTwoIsCC()
<< "\n";
generator()->log() << "sign[0]="<< sign[0]<< "\n";
generator()->log() << "sign[1]="<< sign[1]<< "\n";
generator()->log() << "sign[2]="<< sign[2]<< "\n";
generator()->log() << "sign[3]="<< sign[3]<< "\n";
generator()->log() << "sign[4]="<< sign[4]<< "\n";
generator()->log() << "Recl1[0]="<< Recl1[0] << " Recl3[0]=" << Recl3[0] << "\n";
generator()->log() << "Recr1[0]="<< Recr1[0] << " Recr3[0]=" << Recr3[0] << "\n";
generator()->log() << "Imcl1[0]="<< Imcl1[0] << " Imcl3[0]=" << Imcl3[0] << "\n";
generator()->log() << "Imcr1[0]="<< Imcr1[0] << " Imcr3[0]=" << Imcr3[0] << "\n";
generator()->log() << "Mass[0]="<< mass[0] << "\n";
generator()->log() << "Recl1[1]="<< Recl1[1] << " Recl3[0]=" << Recl3[1] << "\n";
generator()->log() << "Recr1[1]="<< Recr1[1] << " Recr3[0]=" << Recr3[1] << "\n";
generator()->log() << "Imcl1[1]="<< Imcl1[1] << " Imcl3[0]=" << Imcl3[1] << "\n";
generator()->log() << "Imcr1[1]="<< Imcr1[1] << " Imcr3[0]=" << Imcr3[1] << "\n";
generator()->log() << "Mass[1]="<< mass[1] << "\n";
generator()->log() << "born=" << born << "\n";
- generator()->log() << "me2=" << me2()*LastMatchboxXCombInfo::crossingSign() << "\n";
+ generator()->log() << "me2=" << me2() << "\n";
#endif
/*
- double mm = LastMatchboxXCombInfo::crossingSign()*me2();
+ double mm = me2();
born *= lastSHat()/GeV2;
double delta = abs((mm-born)/(mm+born));
if ( mm < 0.0 || born < 0.0 ) {
cerr << "in ";
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p )
cerr << (**p).PDGName() << " ";
cerr << setprecision(20);
cerr << "\nfatal : mm = " << mm << " born = " << born << "\n" << flush;
}
if ( delta > 1.e-12 ) {
cerr << "in ";
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p )
cerr << (**p).PDGName() << " ";
cerr << setprecision(20);
cerr << "\nproblem : mm = " << mm << " born = " << born
<< " -> delta = " << delta << "\n" << flush;
}
*/
- return LastMatchboxXCombInfo::crossingSign()*res;
+ lastOneLoopInterference(res);
+ return lastOneLoopInterference();
}
double Amplitudehqqbarkkbar::oneLoopDoublePole() const {
const map<int, pair<int, double> >& transmap = virtualInfo();
double pbar1[4];
double pbar2[4];
double pbar3[4];
double pbar4[4];
double pbar5[4];
convert(meMomenta()[transmap.find(1)->second.first],pbar1);
convert(meMomenta()[transmap.find(2)->second.first],pbar2);
convert(meMomenta()[transmap.find(3)->second.first],pbar3);
convert(meMomenta()[transmap.find(4)->second.first],pbar4);
convert(meMomenta()[transmap.find(5)->second.first],pbar5);
int sign[5];
for ( size_t k = 0; k < 5; ++k )
sign[k] = transmap.find(k+1)->second.second;
double Recl1[2];
double Recr1[2];
double Recl3[2];
double Recr3[2];
double Reghvv[2];
double Imcl1[2];
double Imcr1[2];
double Imcl3[2];
double Imcr3[2];
double Imghvv[2];
double mass[2];
double width[2];
int IsVaW[2];
getCouplings(Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,Reghvv,Imghvv,mass,width,IsVaW);
//int NF = nLight();
int divMax = 2;
double born;
double virt[3];
double scale = mu2()/GeV2;
F77_qqhqqvbf (pbar1,pbar2,pbar3,pbar4,pbar5,sign,Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,mass,width,IsVaW,divMax,Reghvv,Imghvv,scale,born,virt);
double res = (SM().alphaS()/(2.*Constants::pi))*virt[2]*(lastSHat()/GeV2);
#ifdef AMPVERBOSE
generator()->log() << "double pole="<< res/4./3./3. << "\n";
#endif
- return LastMatchboxXCombInfo::crossingSign()*res;
+ return res;
}
double Amplitudehqqbarkkbar::oneLoopSinglePole() const {
const map<int, pair<int, double> >& transmap = virtualInfo();
double pbar1[4];
double pbar2[4];
double pbar3[4];
double pbar4[4];
double pbar5[4];
convert(meMomenta()[transmap.find(1)->second.first],pbar1);
convert(meMomenta()[transmap.find(2)->second.first],pbar2);
convert(meMomenta()[transmap.find(3)->second.first],pbar3);
convert(meMomenta()[transmap.find(4)->second.first],pbar4);
convert(meMomenta()[transmap.find(5)->second.first],pbar5);
int sign[5];
for ( size_t k = 0; k < 5; ++k )
sign[k] = transmap.find(k+1)->second.second;
double Recl1[2];
double Recr1[2];
double Recl3[2];
double Recr3[2];
double Reghvv[2];
double Imcl1[2];
double Imcr1[2];
double Imcl3[2];
double Imcr3[2];
double Imghvv[2];
double mass[2];
double width[2];
int IsVaW[2];
getCouplings(Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,Reghvv,Imghvv,mass,width,IsVaW);
//int NF = nLight();
int divMax = 2;
double born;
double virt[3];
double scale = mu2()/GeV2;
F77_qqhqqvbf (pbar1,pbar2,pbar3,pbar4,pbar5,sign,Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,
mass,width,IsVaW,divMax,Reghvv,Imghvv,scale,born,virt);
double res = (SM().alphaS()/(2.*Constants::pi))*virt[1]*(lastSHat()/GeV2);
#ifdef AMPVERBOSE
generator()->log() << "single pole="<< res/4./3./3. << "\n";
#endif
- return LastMatchboxXCombInfo::crossingSign()*res;
+ return res;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Amplitudehqqbarkkbar::persistentOutput(PersistentOStream &) const {}
void Amplitudehqqbarkkbar::persistentInput(PersistentIStream &, int) {}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<Amplitudehqqbarkkbar,AmplitudeBase>
describeHJetsAmplitudehqqbarkkbar("HJets::Amplitudehqqbarkkbar", "HwMatchboxBuiltin.so HJets.so");
void Amplitudehqqbarkkbar::Init() {
static ClassDocumentation<Amplitudehqqbarkkbar> documentation
("Helicity amplitudes for 0 -> h q qbar k kbar");
}
diff --git a/Amplitudes/Amplitudehqqbarkkbarg.cc b/Amplitudes/Amplitudehqqbarkkbarg.cc
--- a/Amplitudes/Amplitudehqqbarkkbarg.cc
+++ b/Amplitudes/Amplitudehqqbarkkbarg.cc
@@ -1,691 +1,696 @@
// -*- C++ -*-
//
// Amplitudehqqbarkkbarg.cc is a part of HJets
// Copyright (C) 2011-2012
// Ken Arnold, Francisco Campanario, Terrance Figy, Simon Platzer and Malin Sjodahl
//
// HJets 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 Amplitudehqqbarkkbarg class.
//
#include "Amplitudehqqbarkkbarg.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
extern "C"
{
void qqhqqgvbf_(double PBAR1[4], double PBAR2[4], double PBAR3[4],
double PBAR4[4], double PBAR5[4], double PBAR6[4],
int SIGN[6], double ReCL1[2], double ReCR1[2],
double ReCL3[2], double ReCR3[2],
double ImCL1[2], double ImCR1[2],
double ImCL3[2], double ImCR3[2],
double Mass[2], double Width[2], int IsVaW[2],
int& divMax, int gaugegood[3], int gaugebad[3], double& gaugethres,
double& gaugeacc,
double ReGHVV[2], double ImGHVV[2], double& scale, double& NF,
double& born, double virt[3]);
}
inline void F77_qqhqqgvbf(double PBAR1[4], double PBAR2[4], double PBAR3[4],
double PBAR4[4], double PBAR5[4], double PBAR6[4],
int SIGN[6], double ReCL1[2], double ReCR1[2],
double ReCL3[2], double ReCR3[2],
double ImCL1[2], double ImCR1[2],
double ImCL3[2], double ImCR3[2],
double Mass[2], double Width[2], int IsVaW[2],
int& divMax, int gaugegood[3], int gaugebad[3], double& gaugethres,
double& gaugeacc,
double ReGHVV[2], double ImGHVV[2], double& scale, double& NF,
double& born, double virt[3]){
qqhqqgvbf_(PBAR1,PBAR2,PBAR3,PBAR4,PBAR5,PBAR6,SIGN,ReCL1,ReCR1,ReCL3,ReCR3,
ImCL1,ImCR1,ImCL3,ImCR3,Mass,Width,IsVaW,divMax,
gaugegood,gaugebad,gaugethres,gaugeacc,ReGHVV,ImGHVV,scale,NF,born,virt);
return;
}
using namespace HJets;
struct GaugeStatistics {
unsigned long gaugePassed[3];
unsigned long gaugeFailed[3];
GaugeStatistics();
~GaugeStatistics();
};
GaugeStatistics::GaugeStatistics() {
for ( unsigned int i = 0; i < 3; ++i ) {
gaugePassed[i] = 0;
gaugeFailed[i] = 0;
}
}
GaugeStatistics::~GaugeStatistics() {
if ( (double)gaugePassed[0] > 0 ||
(double)gaugePassed[1] > 0 ||
(double)gaugePassed[2] > 0 ) {
cerr << "\ngauge check summary\n"
<< "--------------------------------------------------------------------------------\n";
if ( gaugePassed[0] > 0 )
cerr << 100*((double)gaugeFailed[0]/(double)gaugePassed[0]) << "% points failing box gauge checks\n";
if ( gaugePassed[1] > 0 )
cerr << 100*((double)gaugeFailed[1]/(double)gaugePassed[1]) << "% points failing abelian hexagon gauge checks\n";
if ( gaugePassed[2] > 0 )
cerr << 100*((double)gaugeFailed[2]/(double)gaugePassed[2]) << "% points failing non-abelian hexagon gauge checks\n";
cerr << "--------------------------------------------------------------------------------\n"
<< flush;
} else if ( (double)gaugeFailed[0] > 0 ||
(double)gaugeFailed[1] > 0 ||
(double)gaugeFailed[2] > 0 ) {
cerr << "\ngauge check fatal: no points passed, some failed!\n" << flush;
}
}
static GaugeStatistics theGaugeStatistics;
Amplitudehqqbarkkbarg::Amplitudehqqbarkkbarg()
: theGaugeThreshold(0.1), theAccuracyCheck(false) {}
Amplitudehqqbarkkbarg::~Amplitudehqqbarkkbarg() {}
IBPtr Amplitudehqqbarkkbarg::clone() const {
return new_ptr(*this);
}
IBPtr Amplitudehqqbarkkbarg::fullclone() const {
return new_ptr(*this);
}
void Amplitudehqqbarkkbarg::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
// inform the amplitude cache how many legs we are handling
nPoints(6);
// everything is expressed in dimensionless quantities, referring to
// sqrt(shat) of the collission as the mass unit
amplitudeScale(sqrt(lastSHat()));
// the higgs momentum
Lorentz5Momentum ph = amplitudeMomentum(0);
momentum(0,ph,true,ph.mass());
// the parton momenta
for ( size_t k = 1; k < 6; ++k )
momentum(k,amplitudeMomentum(k),true);
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex Amplitudehqqbarkkbarg::evaluate(size_t colourTensor,
const vector<int>& helicities,
Complex& largeN) {
// get configurations we need to consider
const vector<AmplitudeInfo>& confs = amplitudeInfo();
// sum of all contributions
Complex result = 0;
largeN = 0.;
for ( vector<AmplitudeInfo>::const_iterator c = confs.begin();
c != confs.end(); ++c ) {
map<size_t,double>::const_iterator cx = c->colourTensors.find(colourTensor);
if ( cx == c->colourTensors.end() )
continue;
Complex diagram = 0.;
Complex ijL = c->ijLeftCoupling;
Complex ijR = c->ijRightCoupling;
Complex klL = c->klLeftCoupling;
Complex klR = c->klRightCoupling;
int i = c->ijLine.first;
int j = c->ijLine.second;
int k = c->klLine.first;
int l = c->klLine.second;
if ( c->ijLineEmissions.first > 0 && c->ijLineEmissions.second < 0 &&
c->klLineEmissions.first < 0 && c->klLineEmissions.second < 0 ) {
int g = c->ijLineEmissions.first;
if ( ijL != 0. && klL != 0. ) {
diagram +=
ijL*klL*qqbargLeftCurrent(i,helicities[i],
j,helicities[j],
g,helicities[g]).
dot(qqbarLeftCurrent(k,helicities[k],
l,helicities[l]));
}
if ( ijL != 0. && klR != 0. ) {
diagram +=
ijL*klR*qqbargLeftCurrent(i,helicities[i],
j,helicities[j],
g,helicities[g]).
dot(qqbarRightCurrent(k,helicities[k],
l,helicities[l]));
}
if ( ijR != 0. && klL != 0. ) {
diagram +=
ijR*klL*qqbargRightCurrent(i,helicities[i],
j,helicities[j],
g,helicities[g]).
dot(qqbarLeftCurrent(k,helicities[k],
l,helicities[l]));
}
if ( ijR != 0. && klR != 0. ) {
diagram +=
ijR*klR*qqbargRightCurrent(i,helicities[i],
j,helicities[j],
g,helicities[g]).
dot(qqbarRightCurrent(k,helicities[k],
l,helicities[l]));
}
} else if ( c->ijLineEmissions.first < 0 && c->ijLineEmissions.second < 0 &&
c->klLineEmissions.first > 0 && c->klLineEmissions.second < 0 ) {
int g = c->klLineEmissions.first;
if ( ijL != 0. && klL != 0. ) {
diagram +=
ijL*klL*qqbarLeftCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbargLeftCurrent(k,helicities[k],
l,helicities[l],
g,helicities[g]));
}
if ( ijL != 0. && klR != 0. ) {
diagram +=
ijL*klR*qqbarLeftCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbargRightCurrent(k,helicities[k],
l,helicities[l],
g,helicities[g]));
}
if ( ijR != 0. && klL != 0. ) {
diagram +=
ijR*klL*qqbarRightCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbargLeftCurrent(k,helicities[k],
l,helicities[l],
g,helicities[g]));
}
if ( ijR != 0. && klR != 0. ) {
diagram +=
ijR*klR*qqbarRightCurrent(i,helicities[i],
j,helicities[j]).
dot(qqbargRightCurrent(k,helicities[k],
l,helicities[l],
g,helicities[g]));
}
} else assert(false);
diagram *= bosonFactor(*c)*cx->second * c->fermionSign;
result += diagram;
if ( cx->second > 0. )
largeN += diagram;
}
largeN *= sqrt(4.*Constants::pi*SM().alphaS());
return sqrt(4.*Constants::pi*SM().alphaS())*result;
}
inline void convert(const Lorentz5Momentum& p, double* res) {
res[0] = p.t()/GeV;
res[1] = p.x()/GeV;
res[2] = p.y()/GeV;
res[3] = p.z()/GeV;
}
double Amplitudehqqbarkkbarg::oneLoopInterference() const {
+ if ( !calculateOneLoopInterference() )
+ return lastOneLoopInterference();
+
const map<int, pair<int, double> >& transmap = virtualInfo();
double pbar1[4];
double pbar2[4];
double pbar3[4];
double pbar4[4];
double pbar5[4];
double pbar6[4];
convert(meMomenta()[transmap.find(1)->second.first],pbar1);
convert(meMomenta()[transmap.find(3)->second.first],pbar3);
convert(meMomenta()[transmap.find(2)->second.first],pbar2);
convert(meMomenta()[transmap.find(4)->second.first],pbar4);
convert(meMomenta()[transmap.find(5)->second.first],pbar5);
convert(meMomenta()[transmap.find(6)->second.first],pbar6);
int sign[6];
for ( size_t k = 0; k < 6; ++k )
sign[k] = (int) transmap.find(k+1)->second.second;
double Recl1[2];
double Recr1[2];
double Recl3[2];
double Recr3[2];
double Reghvv[2];
double Imcl1[2];
double Imcr1[2];
double Imcl3[2];
double Imcr3[2];
double Imghvv[2];
double mass[2];
double width[2];
int IsVaW[2];
getCouplings(Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,Reghvv,Imghvv,
mass,width,IsVaW);
double NF = nLight();
int gaugegood[3];
int gaugebad[3];
int divMax = 2;
double born;
double virt[3];
double scale = mu2()/GeV2;
double gt = theGaugeThreshold;
double gacc(0.);
#ifdef AMPVERBOSE
generator()->log()<<"scale="<<scale << "\n";
#endif
F77_qqhqqgvbf (pbar1,pbar2,pbar3,pbar4,pbar5,pbar6,sign,Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,mass,width,IsVaW,divMax,gaugegood,gaugebad,gt,gacc,
Reghvv,Imghvv,scale,NF,born,virt);
for (int i = 0; i < 3; i++) {
theGaugeStatistics.gaugePassed[i] += gaugegood[i];
theGaugeStatistics.gaugeFailed[i] += gaugebad[i];
}
double res = 2.*sqr(SM().alphaS())*virt[0]*pow(lastSHat()/GeV2,2);
if (gaugebad[0]+gaugebad[1]+gaugebad[2]) {
throw Veto();
}
if( theAccuracyCheck ){
generator()->log() << "checkacc "<< log10 (abs(gacc))<<"\n";
}
#ifdef AMPVERBOSE
const cPDVector& proc = mePartonData();
map<int, pair<int, double> > XMap = virtualInfos().find(proc)->second;
generator()->log() << setprecision(20);
generator()->log() << "--------------------- \n";
generator()->log() << "c2="<< virt[2]/born << "\n";
generator()->log() << "c1="<< virt[1]/born << "\n";
generator()->log() << "c0="<< virt[0]/born << "\n";
double gs2(SM().alphaS()*4.*Constants::pi);
born*=pow(lastSHat()/GeV2,2)*gs2;
- generator()->log() << "rb=" << born/me2()/LastMatchboxXCombInfo::crossingSign() << " Base Diagram:"
+ generator()->log() << "rb=" << born/me2() << " Base Diagram:"
<< (XMap.find(1)->second.second > 0 ? "+" : "-") << proc[XMap.find(1)->second.first]->PDGName() << " "
<< (XMap.find(3)->second.second > 0 ? "+" : "-") << proc[XMap.find(3)->second.first]->PDGName() << " -> "
<< (XMap.find(2)->second.second > 0 ? "+" : "-") << proc[XMap.find(2)->second.first]->PDGName() << " "
<< (XMap.find(4)->second.second > 0 ? "+" : "-") << proc[XMap.find(4)->second.first]->PDGName() << " "
<< (XMap.find(5)->second.second > 0 ? "+" : "-") << proc[XMap.find(5)->second.first]->PDGName() << " "
<< (XMap.find(6)->second.second > 0 ? "+" : "-") << proc[XMap.find(6)->second.first]->PDGName() << "\n";
generator()->log() << "OneisNC= " << topologyOneIsNC()
<< " OneisCC= " << topologyOneIsCC()
<< " TwoisNC= " << topologyTwoIsNC()
<< " TwoisCC= " << topologyTwoIsCC()
<< "\n";
generator()->log() << "sign[0]="<< sign[0]<< " xmap(1)="<< XMap.find(1)->second.second << "\n";
generator()->log() << "sign[2]="<< sign[2]<< " xmap(3)="<< XMap.find(3)->second.second << "\n";
generator()->log() << "sign[1]="<< sign[1]<< "\n";
generator()->log() << "sign[3]="<< sign[3]<< "\n";
generator()->log() << "sign[4]="<< sign[4]<< "\n";
generator()->log() << "sign[5]="<< sign[5]<< "\n";
generator()->log() << "Recl1[0]="<< Recl1[0] << " Recl3[0]=" << Recl3[0] << "\n";
generator()->log() << "Recr1[0]="<< Recr1[0] << " Recr3[0]=" << Recr3[0] << "\n";
generator()->log() << "Imcl1[0]="<< Imcl1[0] << " Imcl3[0]=" << Imcl3[0] << "\n";
generator()->log() << "Imcr1[0]="<< Imcr1[0] << " Imcr3[0]=" << Imcr3[0] << "\n";
generator()->log() << "Reghvv[0]="<< Reghvv[0]<< " Imghvv[0]=" << Imghvv[0] << "\n";
generator()->log() << "Mass[0]="<< mass[0] << "\n";
generator()->log() << "Recl1[1]="<< Recl1[1] << " Recl3[0]=" << Recl3[1] << "\n";
generator()->log() << "Recr1[1]="<< Recr1[1] << " Recr3[0]=" << Recr3[1] << "\n";
generator()->log() << "Imcl1[1]="<< Imcl1[1] << " Imcl3[0]=" << Imcl3[1] << "\n";
generator()->log() << "Imcr1[1]="<< Imcr1[1] << " Imcr3[0]=" << Imcr3[1] << "\n";
generator()->log() << "Reghvv[1]="<< Reghvv[1] << " Imghvv[1]="<<Imghvv[1] << "\n";
generator()->log() << "Mass[1]="<< mass[1] << "\n";
generator()->log() << "born=" << born << "\n";
- generator()->log() << "me2=" << me2()*LastMatchboxXCombInfo::crossingSign() << "\n";
+ generator()->log() << "me2=" << me2() << "\n";
generator()->log() << "--------------------- \n";
#endif
- return LastMatchboxXCombInfo::crossingSign()*res;
+ lastOneLoopInterference(res);
+
+ return lastOneLoopInterference();
}
double Amplitudehqqbarkkbarg::oneLoopDoublePole() const {
const map<int, pair<int, double> >& transmap = virtualInfo();
double pbar1[4];
double pbar2[4];
double pbar3[4];
double pbar4[4];
double pbar5[4];
double pbar6[4];
convert(meMomenta()[transmap.find(1)->second.first],pbar1);
convert(meMomenta()[transmap.find(3)->second.first],pbar3);
convert(meMomenta()[transmap.find(2)->second.first],pbar2);
convert(meMomenta()[transmap.find(4)->second.first],pbar4);
convert(meMomenta()[transmap.find(5)->second.first],pbar5);
convert(meMomenta()[transmap.find(6)->second.first],pbar6);
int sign[6];
for ( size_t k = 0; k < 6; ++k )
sign[k] = (int) transmap.find(k+1)->second.second;
double Recl1[2];
double Recr1[2];
double Recl3[2];
double Recr3[2];
double Reghvv[2];
double Imcl1[2];
double Imcr1[2];
double Imcl3[2];
double Imcr3[2];
double Imghvv[2];
double mass[2];
double width[2];
int IsVaW[2];
getCouplings(Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,Reghvv,Imghvv,
mass,width,IsVaW);
double NF = nLight();
int gaugegood[3];
int gaugebad[3];
int divMax = 2;
double born;
double virt[3];
double scale = mu2()/GeV2;
double gt = theGaugeThreshold;
double gacc(0.);
#ifdef AMPVERBOSE
generator()->log()<<"scale="<<scale << "\n";
#endif
F77_qqhqqgvbf (pbar1,pbar2,pbar3,pbar4,pbar5,pbar6,sign,Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,mass,width,IsVaW,divMax,gaugegood,gaugebad,
gt,gacc,Reghvv,Imghvv,
scale,NF,born,virt);
double res = 2.*sqr(SM().alphaS())*virt[2]*pow(lastSHat()/GeV2,2);
if (gaugebad[0]+gaugebad[1]+gaugebad[2]) {
throw Veto();
}
/*generator()->log() << "c2="<< virt[2]/born << "\n";
generator()->log() << "c1="<< virt[1]/born << "\n";
generator()->log() << "c0="<< virt[0]/born << "\n";
*/
- return LastMatchboxXCombInfo::crossingSign()*res;
+ return res;
}
double Amplitudehqqbarkkbarg::oneLoopSinglePole() const {
const map<int, pair<int, double> >& transmap = virtualInfo();
double pbar1[4];
double pbar2[4];
double pbar3[4];
double pbar4[4];
double pbar5[4];
double pbar6[4];
convert(meMomenta()[transmap.find(1)->second.first],pbar1);
convert(meMomenta()[transmap.find(3)->second.first],pbar3);
convert(meMomenta()[transmap.find(2)->second.first],pbar2);
convert(meMomenta()[transmap.find(4)->second.first],pbar4);
convert(meMomenta()[transmap.find(5)->second.first],pbar5);
convert(meMomenta()[transmap.find(6)->second.first],pbar6);
int sign[6];
for ( size_t k = 0; k < 6; ++k )
sign[k] = (int) transmap.find(k+1)->second.second;
double Recl1[2];
double Recr1[2];
double Recl3[2];
double Recr3[2];
double Reghvv[2];
double Imcl1[2];
double Imcr1[2];
double Imcl3[2];
double Imcr3[2];
double Imghvv[2];
double mass[2];
double width[2];
int IsVaW[2];
getCouplings(Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,Reghvv,Imghvv,
mass,width,IsVaW);
double NF = nLight();
int gaugegood[3];
int gaugebad[3];
int divMax = 2;
double born;
double virt[3];
double scale = mu2()/GeV2;
double gt = theGaugeThreshold;
double gacc(0.);
#ifdef AMPVERBOSE
generator()->log()<<"scale="<<scale << "\n";
#endif
F77_qqhqqgvbf (pbar1,pbar2,pbar3,pbar4,pbar5,pbar6,sign,Recl1,Recr1,Recl3,Recr3,
Imcl1,Imcr1,Imcl3,Imcr3,
mass,width,IsVaW,divMax,gaugegood,gaugebad,gt,gacc,Reghvv,Imghvv,
scale,NF,born,virt);
double res = 2.*sqr(SM().alphaS())*virt[1]*pow(lastSHat()/GeV2,2);
if (gaugebad[0]+gaugebad[1]+gaugebad[2]) {
throw Veto();
}
/*
generator()->log() << "NLight= "<<NF<< "\n";
generator()->log() << "c2="<< virt[2]/born << "\n";
generator()->log() << "c1="<< virt[1]/born << "\n";
generator()->log() << "c0="<< virt[0]/born << "\n";
const cPDVector& proc = mePartonData();
double CA1 = 3.;
double CF1 = 4./3.;
double msq = me2();
double two(2.);
double cTaTb = colourCorrelatedME2(pair<int,int>(0,1))*
(proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/abs((meMomenta()[0] + meMomenta()[1]).m2()));
double cT1T2 = colourCorrelatedME2(pair<int,int>(2,3))*
(proc[2]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/abs((meMomenta()[2] + meMomenta()[3]).m2()));
double cT1T3 = colourCorrelatedME2(pair<int,int>(2,4))*
(proc[2]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[2] + meMomenta()[4]).m2())));
double cT2T3 = colourCorrelatedME2(pair<int,int>(3,4))*
(proc[3]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[3] + meMomenta()[4]).m2())));
double cTaT1 = colourCorrelatedME2(pair<int,int>(0,2))*
(proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[0] + meMomenta()[2]).m2())));
double cTaT2 = colourCorrelatedME2(pair<int,int>(0,3))*
(proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[0] + meMomenta()[3]).m2())));
double cTaT3 = colourCorrelatedME2(pair<int,int>(0,4))*
(proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[0] + meMomenta()[4]).m2())));
double cTbT1 = colourCorrelatedME2(pair<int,int>(1,2))*
(proc[1]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[1] + meMomenta()[2]).m2())));
double cTbT2 = colourCorrelatedME2(pair<int,int>(1,3))*
(proc[1]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[1] + meMomenta()[3]).m2())));
double cTbT3 = colourCorrelatedME2(pair<int,int>(1,4))*
(proc[1]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[1] + meMomenta()[4]).m2())));
generator()->log() <<"CA1"<<CA1<< "\n";
double gammaQuark = (3./2.)*CF1;
double gammaGluon = (11./6.)*CA1 - (1./3.)*NF;
double cfac(4.*gammaQuark+gammaGluon);
double res1 =-cfac*msq+2.*(cT1T2 + cT1T3 + cTaT1 + cTbT1 + cT2T3 + cTaT2+ cTbT2 + cTaT3 + cTbT3 + cTaTb);
double gs2(SM().alphaS()*4.*Constants::pi);
born*=pow(lastSHat()/GeV2,2)*gs2;
generator()->log()<<"res1/born="<<res1/me2()<<"\n";
*/
/*
- double mm = LastMatchboxXCombInfo::crossingSign()*me2();
+ double mm = me2();
born *= pow(lastSHat()/GeV2,2);
double delta = abs((mm-born)/(mm+born));
if ( mm < 0.0 || born < 0.0 ) {
cerr << "in ";
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p )
cerr << (**p).PDGName() << " ";
cerr << setprecision(20);
cerr << "\nfatal : mm = " << mm << " born = " << born << "\n" << flush;
}
if ( delta > 1.e-12 ) {
cerr << "in ";
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p )
cerr << (**p).PDGName() << " ";
cerr << setprecision(20);
cerr << "\nproblem : mm = " << mm << " born = " << born
<< " -> delta = " << delta << "\n" << flush;
}
*/
- return LastMatchboxXCombInfo::crossingSign()*res;
+ return res;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Amplitudehqqbarkkbarg::persistentOutput(PersistentOStream & os) const {
os << theGaugeThreshold;
}
void Amplitudehqqbarkkbarg::persistentInput(PersistentIStream & is, int) {
is >> theGaugeThreshold;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<Amplitudehqqbarkkbarg,AmplitudeBase>
describeHJetsAmplitudehqqbarkkbarg("HJets::Amplitudehqqbarkkbarg", "HwMatchboxBuiltin.so HJets.so");
void Amplitudehqqbarkkbarg::Init() {
static ClassDocumentation<Amplitudehqqbarkkbarg> documentation
("Helicity amplitudes for 0 -> h q qbar k kbar g");
static Parameter<Amplitudehqqbarkkbarg,double> interfaceGaugeThreshold
("GaugeThreshold",
"The threshold for gauge check failure.",
&Amplitudehqqbarkkbarg::theGaugeThreshold, 0.1, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<Amplitudehqqbarkkbarg,bool> interfaceAccuracyCheck
("AccuracyCheck",
"Turn on check of accuracy of box,pentagon, and hexagons.",
&Amplitudehqqbarkkbarg::theAccuracyCheck, false, false, false);
static SwitchOption interfaceAccuracyCheckOn
(interfaceAccuracyCheck,
"On",
"On",
true);
static SwitchOption interfaceAccuracyCheckOff
(interfaceAccuracyCheck,
"Off",
"Off",
false);
}
diff --git a/INSTALL b/INSTALL
deleted file mode 100644
--- a/INSTALL
+++ /dev/null
@@ -1,365 +0,0 @@
-Installation Instructions
-*************************
-
-Copyright (C) 1994, 1995, 1996, 1999, 2000, 2001, 2002, 2004, 2005,
-2006, 2007, 2008, 2009 Free Software Foundation, Inc.
-
- Copying and distribution of this file, with or without modification,
-are permitted in any medium without royalty provided the copyright
-notice and this notice are preserved. This file is offered as-is,
-without warranty of any kind.
-
-Basic Installation
-==================
-
- Briefly, the shell commands `./configure; make; make install' should
-configure, build, and install this package. The following
-more-detailed instructions are generic; see the `README' file for
-instructions specific to this package. Some packages provide this
-`INSTALL' file but do not implement all of the features documented
-below. The lack of an optional feature in a given package is not
-necessarily a bug. More recommendations for GNU packages can be found
-in *note Makefile Conventions: (standards)Makefile Conventions.
-
- The `configure' shell script attempts to guess correct values for
-various system-dependent variables used during compilation. It uses
-those values to create a `Makefile' in each directory of the package.
-It may also create one or more `.h' files containing system-dependent
-definitions. Finally, it creates a shell script `config.status' that
-you can run in the future to recreate the current configuration, and a
-file `config.log' containing compiler output (useful mainly for
-debugging `configure').
-
- It can also use an optional file (typically called `config.cache'
-and enabled with `--cache-file=config.cache' or simply `-C') that saves
-the results of its tests to speed up reconfiguring. Caching is
-disabled by default to prevent problems with accidental use of stale
-cache files.
-
- If you need to do unusual things to compile the package, please try
-to figure out how `configure' could check whether to do them, and mail
-diffs or instructions to the address given in the `README' so they can
-be considered for the next release. If you are using the cache, and at
-some point `config.cache' contains results you don't want to keep, you
-may remove or edit it.
-
- The file `configure.ac' (or `configure.in') is used to create
-`configure' by a program called `autoconf'. You need `configure.ac' if
-you want to change it or regenerate `configure' using a newer version
-of `autoconf'.
-
- The simplest way to compile this package is:
-
- 1. `cd' to the directory containing the package's source code and type
- `./configure' to configure the package for your system.
-
- Running `configure' might take a while. While running, it prints
- some messages telling which features it is checking for.
-
- 2. Type `make' to compile the package.
-
- 3. Optionally, type `make check' to run any self-tests that come with
- the package, generally using the just-built uninstalled binaries.
-
- 4. Type `make install' to install the programs and any data files and
- documentation. When installing into a prefix owned by root, it is
- recommended that the package be configured and built as a regular
- user, and only the `make install' phase executed with root
- privileges.
-
- 5. Optionally, type `make installcheck' to repeat any self-tests, but
- this time using the binaries in their final installed location.
- This target does not install anything. Running this target as a
- regular user, particularly if the prior `make install' required
- root privileges, verifies that the installation completed
- correctly.
-
- 6. You can remove the program binaries and object files from the
- source code directory by typing `make clean'. To also remove the
- files that `configure' created (so you can compile the package for
- a different kind of computer), type `make distclean'. There is
- also a `make maintainer-clean' target, but that is intended mainly
- for the package's developers. If you use it, you may have to get
- all sorts of other programs in order to regenerate files that came
- with the distribution.
-
- 7. Often, you can also type `make uninstall' to remove the installed
- files again. In practice, not all packages have tested that
- uninstallation works correctly, even though it is required by the
- GNU Coding Standards.
-
- 8. Some packages, particularly those that use Automake, provide `make
- distcheck', which can by used by developers to test that all other
- targets like `make install' and `make uninstall' work correctly.
- This target is generally not run by end users.
-
-Compilers and Options
-=====================
-
- Some systems require unusual options for compilation or linking that
-the `configure' script does not know about. Run `./configure --help'
-for details on some of the pertinent environment variables.
-
- You can give `configure' initial values for configuration parameters
-by setting variables in the command line or in the environment. Here
-is an example:
-
- ./configure CC=c99 CFLAGS=-g LIBS=-lposix
-
- *Note Defining Variables::, for more details.
-
-Compiling For Multiple Architectures
-====================================
-
- You can compile the package for more than one kind of computer at the
-same time, by placing the object files for each architecture in their
-own directory. To do this, you can use GNU `make'. `cd' to the
-directory where you want the object files and executables to go and run
-the `configure' script. `configure' automatically checks for the
-source code in the directory that `configure' is in and in `..'. This
-is known as a "VPATH" build.
-
- With a non-GNU `make', it is safer to compile the package for one
-architecture at a time in the source code directory. After you have
-installed the package for one architecture, use `make distclean' before
-reconfiguring for another architecture.
-
- On MacOS X 10.5 and later systems, you can create libraries and
-executables that work on multiple system types--known as "fat" or
-"universal" binaries--by specifying multiple `-arch' options to the
-compiler but only a single `-arch' option to the preprocessor. Like
-this:
-
- ./configure CC="gcc -arch i386 -arch x86_64 -arch ppc -arch ppc64" \
- CXX="g++ -arch i386 -arch x86_64 -arch ppc -arch ppc64" \
- CPP="gcc -E" CXXCPP="g++ -E"
-
- This is not guaranteed to produce working output in all cases, you
-may have to build one architecture at a time and combine the results
-using the `lipo' tool if you have problems.
-
-Installation Names
-==================
-
- By default, `make install' installs the package's commands under
-`/usr/local/bin', include files under `/usr/local/include', etc. You
-can specify an installation prefix other than `/usr/local' by giving
-`configure' the option `--prefix=PREFIX', where PREFIX must be an
-absolute file name.
-
- You can specify separate installation prefixes for
-architecture-specific files and architecture-independent files. If you
-pass the option `--exec-prefix=PREFIX' to `configure', the package uses
-PREFIX as the prefix for installing programs and libraries.
-Documentation and other data files still use the regular prefix.
-
- In addition, if you use an unusual directory layout you can give
-options like `--bindir=DIR' to specify different values for particular
-kinds of files. Run `configure --help' for a list of the directories
-you can set and what kinds of files go in them. In general, the
-default for these options is expressed in terms of `${prefix}', so that
-specifying just `--prefix' will affect all of the other directory
-specifications that were not explicitly provided.
-
- The most portable way to affect installation locations is to pass the
-correct locations to `configure'; however, many packages provide one or
-both of the following shortcuts of passing variable assignments to the
-`make install' command line to change installation locations without
-having to reconfigure or recompile.
-
- The first method involves providing an override variable for each
-affected directory. For example, `make install
-prefix=/alternate/directory' will choose an alternate location for all
-directory configuration variables that were expressed in terms of
-`${prefix}'. Any directories that were specified during `configure',
-but not in terms of `${prefix}', must each be overridden at install
-time for the entire installation to be relocated. The approach of
-makefile variable overrides for each directory variable is required by
-the GNU Coding Standards, and ideally causes no recompilation.
-However, some platforms have known limitations with the semantics of
-shared libraries that end up requiring recompilation when using this
-method, particularly noticeable in packages that use GNU Libtool.
-
- The second method involves providing the `DESTDIR' variable. For
-example, `make install DESTDIR=/alternate/directory' will prepend
-`/alternate/directory' before all installation names. The approach of
-`DESTDIR' overrides is not required by the GNU Coding Standards, and
-does not work on platforms that have drive letters. On the other hand,
-it does better at avoiding recompilation issues, and works well even
-when some directory options were not specified in terms of `${prefix}'
-at `configure' time.
-
-Optional Features
-=================
-
- If the package supports it, you can cause programs to be installed
-with an extra prefix or suffix on their names by giving `configure' the
-option `--program-prefix=PREFIX' or `--program-suffix=SUFFIX'.
-
- Some packages pay attention to `--enable-FEATURE' options to
-`configure', where FEATURE indicates an optional part of the package.
-They may also pay attention to `--with-PACKAGE' options, where PACKAGE
-is something like `gnu-as' or `x' (for the X Window System). The
-`README' should mention any `--enable-' and `--with-' options that the
-package recognizes.
-
- For packages that use the X Window System, `configure' can usually
-find the X include and library files automatically, but if it doesn't,
-you can use the `configure' options `--x-includes=DIR' and
-`--x-libraries=DIR' to specify their locations.
-
- Some packages offer the ability to configure how verbose the
-execution of `make' will be. For these packages, running `./configure
---enable-silent-rules' sets the default to minimal output, which can be
-overridden with `make V=1'; while running `./configure
---disable-silent-rules' sets the default to verbose, which can be
-overridden with `make V=0'.
-
-Particular systems
-==================
-
- On HP-UX, the default C compiler is not ANSI C compatible. If GNU
-CC is not installed, it is recommended to use the following options in
-order to use an ANSI C compiler:
-
- ./configure CC="cc -Ae -D_XOPEN_SOURCE=500"
-
-and if that doesn't work, install pre-built binaries of GCC for HP-UX.
-
- On OSF/1 a.k.a. Tru64, some versions of the default C compiler cannot
-parse its `<wchar.h>' header file. The option `-nodtk' can be used as
-a workaround. If GNU CC is not installed, it is therefore recommended
-to try
-
- ./configure CC="cc"
-
-and if that doesn't work, try
-
- ./configure CC="cc -nodtk"
-
- On Solaris, don't put `/usr/ucb' early in your `PATH'. This
-directory contains several dysfunctional programs; working variants of
-these programs are available in `/usr/bin'. So, if you need `/usr/ucb'
-in your `PATH', put it _after_ `/usr/bin'.
-
- On Haiku, software installed for all users goes in `/boot/common',
-not `/usr/local'. It is recommended to use the following options:
-
- ./configure --prefix=/boot/common
-
-Specifying the System Type
-==========================
-
- There may be some features `configure' cannot figure out
-automatically, but needs to determine by the type of machine the package
-will run on. Usually, assuming the package is built to be run on the
-_same_ architectures, `configure' can figure that out, but if it prints
-a message saying it cannot guess the machine type, give it the
-`--build=TYPE' option. TYPE can either be a short name for the system
-type, such as `sun4', or a canonical name which has the form:
-
- CPU-COMPANY-SYSTEM
-
-where SYSTEM can have one of these forms:
-
- OS
- KERNEL-OS
-
- See the file `config.sub' for the possible values of each field. If
-`config.sub' isn't included in this package, then this package doesn't
-need to know the machine type.
-
- If you are _building_ compiler tools for cross-compiling, you should
-use the option `--target=TYPE' to select the type of system they will
-produce code for.
-
- If you want to _use_ a cross compiler, that generates code for a
-platform different from the build platform, you should specify the
-"host" platform (i.e., that on which the generated programs will
-eventually be run) with `--host=TYPE'.
-
-Sharing Defaults
-================
-
- If you want to set default values for `configure' scripts to share,
-you can create a site shell script called `config.site' that gives
-default values for variables like `CC', `cache_file', and `prefix'.
-`configure' looks for `PREFIX/share/config.site' if it exists, then
-`PREFIX/etc/config.site' if it exists. Or, you can set the
-`CONFIG_SITE' environment variable to the location of the site script.
-A warning: not all `configure' scripts look for a site script.
-
-Defining Variables
-==================
-
- Variables not defined in a site shell script can be set in the
-environment passed to `configure'. However, some packages may run
-configure again during the build, and the customized values of these
-variables may be lost. In order to avoid this problem, you should set
-them in the `configure' command line, using `VAR=value'. For example:
-
- ./configure CC=/usr/local2/bin/gcc
-
-causes the specified `gcc' to be used as the C compiler (unless it is
-overridden in the site shell script).
-
-Unfortunately, this technique does not work for `CONFIG_SHELL' due to
-an Autoconf bug. Until the bug is fixed you can use this workaround:
-
- CONFIG_SHELL=/bin/bash /bin/bash ./configure CONFIG_SHELL=/bin/bash
-
-`configure' Invocation
-======================
-
- `configure' recognizes the following options to control how it
-operates.
-
-`--help'
-`-h'
- Print a summary of all of the options to `configure', and exit.
-
-`--help=short'
-`--help=recursive'
- Print a summary of the options unique to this package's
- `configure', and exit. The `short' variant lists options used
- only in the top level, while the `recursive' variant lists options
- also present in any nested packages.
-
-`--version'
-`-V'
- Print the version of Autoconf used to generate the `configure'
- script, and exit.
-
-`--cache-file=FILE'
- Enable the cache: use and save the results of the tests in FILE,
- traditionally `config.cache'. FILE defaults to `/dev/null' to
- disable caching.
-
-`--config-cache'
-`-C'
- Alias for `--cache-file=config.cache'.
-
-`--quiet'
-`--silent'
-`-q'
- Do not print messages saying which checks are being made. To
- suppress all normal output, redirect it to `/dev/null' (any error
- messages will still be shown).
-
-`--srcdir=DIR'
- Look for the package's source code in directory DIR. Usually
- `configure' can determine that directory automatically.
-
-`--prefix=DIR'
- Use DIR as the installation prefix. *note Installation Names::
- for more details, including other options available for fine-tuning
- the installation locations.
-
-`--no-create'
-`-n'
- Run the configure checks, but stop before creating any output
- files.
-
-`configure' also accepts some other, not widely useful, options. Run
-`configure --help' for more details.
-
diff --git a/NEWS b/NEWS
--- a/NEWS
+++ b/NEWS
@@ -0,0 +1,8 @@
+HJets++ News -*- outline -*-
+================================================================================
+
+* HJets 1.1 released
+** a major bug has been fixed as the 3-jet virtual contributions
+ used a convention which has not been supported anymore by Matchbox
+
+
diff --git a/README b/README
--- a/README
+++ b/README
@@ -1,43 +1,44 @@
################################################################################
# #
# HJets++ 1.0 #
# ========================================================================== #
# A Matchbox plugin for electroweak Higgs plus #
# two and three jet production at NLO QCD #
# #
# #
# (C) 2012-2014 The HJets++ collaboration #
# ------------------------------------------------------------ #
# Francisco Campanario <francisco.campanario@ific.uv.es> #
# Terrance M. Figy <terrance.figy@hep.manchester.ac.uk> #
# Simon Platzer <simon.plaetzer@desy.de> #
# Malin Sjodahl <malin.sjodahl@thep.lu.se> #
# ------------------------------------------------------------ #
# #
# #
# ========================================================================== #
# #
# HJets++ is licenced under version 2 of the GPL, see COPYING for details. #
# Please respect the MCnet academic guidelines, see GUIDELINES for details. #
# #
# When using HJets++ please cite: #
# #
# F. Campanario, T.M. Figy, S. Platzer and M. Sjodahl: #
# 'Electroweak Higgs plus Three Jet Production at NLO QCD', #
# Phys.Rev.Lett. 111 (2013) 211802, arXiv:1308.2932 [hep-ph] #
# #
################################################################################
Installation
============
-A Herwig++ installation with versions >= 3.0 is required to use HJets++
+A Herwig++ installation with versions >= 3.0 (i.e. Herwig 7.0) is required to
+use HJets++
--with-herwig=path
-informs HJets++'s configure where to look for the Herwig++
+informs HJets++'s configure where to look for the Herwig 7
installation. Installation then proceeds as usual with make and make
-install. Example usage of HJets++ is provided along with the Herwig++
-releases from version 3.0 onwards.
+install. Example usage of HJets++ is provided along with the Herwig 7
+releases.
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,48 +1,51 @@
dnl Process this file with autoconf to produce a configure script.
AC_PREREQ([2.59])
-AC_INIT([HJets],[SVN],[simon.plaetzer@desy.de],[HJets])
+AC_INIT([HJets],[1.1],[simon.plaetzer@desy.de],[HJets])
AC_CONFIG_AUX_DIR([config])
AC_CONFIG_MACRO_DIR([m4])
-AC_CONFIG_SRCDIR(["m4/$PACKAGE_TARNAME.m4"])
+AC_CONFIG_SRCDIR([Amplitudes/AmplitudeBase.h])
AC_CANONICAL_HOST
if test -n "$F77" ; then
AC_MSG_ERROR([F77 is set internally. Please use FC to specify a non-default Fortran compiler.])
fi
AC_LANG([C++])
AM_INIT_AUTOMAKE([1.9 gnu subdir-objects dist-bzip2 -Wall])
m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
dnl Checks for programs.
+
+AC_PROG_CC
+AC_PROG_CXX
+
AC_PROG_FC([gfortran])
-AC_PROG_F77([gfortran])
-AC_PROG_CXX
+AC_SUBST([F77],[$FC])
+
AC_PROG_INSTALL
AC_PROG_MAKE_SET
AC_PROG_LN_S
-AM_PROG_AR
+m4_ifdef([AM_PROG_AR], [AM_PROG_AR])
-AC_DISABLE_STATIC
-AC_LIBTOOL_DLOPEN
-AC_PROG_LIBTOOL
+LT_PREREQ([2.2.6])
+LT_INIT([disable-static dlopen pic-only])
HJETS_CHECK_HERWIG
HJETS_FORTRAN_FLAGS
AC_CONFIG_FILES([Amplitudes/Makefile
Amplitudes/AmplitudeBase.cc
Virtuals/Amplitudes/Makefile
Virtuals/Integrals/Makefile
Virtuals/MatrixElement/Makefile
Virtuals/TenRed/Makefile
Virtuals/Utilities/Makefile
Virtuals/OneLOop/Makefile
Virtuals/Makefile
include/Makefile
src/Makefile
Makefile])
AC_OUTPUT
diff --git a/m4/HJets.m4 b/m4/HJets.m4
--- a/m4/HJets.m4
+++ b/m4/HJets.m4
@@ -1,48 +1,48 @@
-dnl --- check for Herwig++ --
+dnl --- check for Herwig 7 --
AC_DEFUN([HJETS_CHECK_HERWIG],
[
defaultlocation="${prefix}"
test "x$defaultlocation" = xNONE && defaultlocation="${ac_default_prefix}"
-AC_MSG_CHECKING([for Herwig++/Matchbox in])
+AC_MSG_CHECKING([for Herwig 7 in])
AC_ARG_WITH(herwig,
- AC_HELP_STRING([--with-herwig=DIR],[location of Herwig++ installation]),
+ AC_HELP_STRING([--with-herwig=DIR],[location of Herwig 7 installation]),
[],
[with_herwig="${defaultlocation}"])
AC_MSG_RESULT([$with_herwig])
AS_IF([test "x$with_herwig" != "xno"],
[AC_CHECK_FILES(
${with_herwig}/bin/herwig-config,
[have_herwig=yes], [have_herwig=no])],
[have_herwig=no])
AS_IF([test "x$have_herwig" = "xyes"],
[HERWIGCPPFLAGS=`${with_herwig}/bin/herwig-config --cppflags`
HERWIGLDFLAGS=`${with_herwig}/bin/herwig-config --ldflags`
HERWIGLDLIBS=`${with_herwig}/bin/herwig-config --ldlibs`
],
[AS_IF([test "x$with_herwig" != "xno"],
[AC_MSG_ERROR([Cannot build HJets without Herwig. Please set --with-herwig.])
])
])
AM_CPPFLAGS="-I\$(top_builddir)/include $HERWIGCPPFLAGS"
AM_LDFLAGS="$HERWIGLDFLAGS $HERWIGLDLIBS"
AC_SUBST(AM_CPPFLAGS)
AC_SUBST(AM_LDFLAGS)
])
dnl --- set fortran flags ---
AC_DEFUN([HJETS_FORTRAN_FLAGS],
[
AM_FCFLAGS="$AM_FCFLAGS -fno-automatic -ffixed-line-length-none"
AM_FFLAGS="$AM_FFLAGS -fno-automatic -ffixed-line-length-none"
AC_SUBST(AM_FCFLAGS)
AC_SUBST(AM_FFLAGS)
])
diff --git a/src/HJetsProcesses.in b/src/HJetsProcesses.in
--- a/src/HJetsProcesses.in
+++ b/src/HJetsProcesses.in
@@ -1,25 +1,35 @@
+# -*- ThePEG-repository -*-
+
cd /Herwig/MatrixElements/Matchbox/Amplitudes
mkdir HJets
cd HJets
create HJets::Amplitudehqqbarkkbar Amplitudehqqbarkkbar
set Amplitudehqqbarkkbar:ColourBasis /Herwig/MatrixElements/Matchbox/Amplitudes/TraceBasis
create HJets::Amplitudehqqbarkkbarg Amplitudehqqbarkkbarg
set Amplitudehqqbarkkbarg:ColourBasis /Herwig/MatrixElements/Matchbox/Amplitudes/TraceBasis
create HJets::Amplitudehqqbarkkbargg Amplitudehqqbarkkbargg
set Amplitudehqqbarkkbargg:ColourBasis /Herwig/MatrixElements/Matchbox/Amplitudes/TraceBasis
create HJets::Amplitudehqqbarkkbarrrbar Amplitudehqqbarkkbarrrbar
set Amplitudehqqbarkkbarrrbar:ColourBasis /Herwig/MatrixElements/Matchbox/Amplitudes/TraceBasis
-insert /Herwig/MatrixElements/Matchbox/PPFactory:Amplitudes 0 Amplitudehqqbarkkbar
-insert /Herwig/MatrixElements/Matchbox/PPFactory:Amplitudes 0 Amplitudehqqbarkkbarg
-insert /Herwig/MatrixElements/Matchbox/PPFactory:Amplitudes 0 Amplitudehqqbarkkbargg
-insert /Herwig/MatrixElements/Matchbox/PPFactory:Amplitudes 0 Amplitudehqqbarkkbarrrbar
+insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehqqbarkkbar
+insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehqqbarkkbarg
+insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehqqbarkkbargg
+insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehqqbarkkbarrrbar
insert /Herwig/MatrixElements/Matchbox/Utility/DiagramGenerator:ExcludeVertices 0 /Herwig/Vertices/FFPVertex
insert /Herwig/MatrixElements/Matchbox/Utility/DiagramGenerator:ExcludeVertices 0 /Herwig/Vertices/FFHVertex
-insert /Herwig/MatrixElements/Matchbox/Utility/DiagramGenerator:ExcludeVertices 0 /Herwig/Vertices/HGGVertex
\ No newline at end of file
+insert /Herwig/MatrixElements/Matchbox/Utility/DiagramGenerator:ExcludeVertices 0 /Herwig/Vertices/HGGVertex
+
+cd /Herwig/Particles
+
+set h0:HardProcessWidth 0*GeV
+do W+:UnsetHardProcessWidth
+do W-:UnsetHardProcessWidth
+do Z0:UnsetHardProcessWidth
+

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:31 PM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805839
Default Alt Text
(97 KB)

Event Timeline