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 >& AmplitudeBase::amplitudeInfos() { static map > theAmplitudeInfos; return theAmplitudeInfos; } map > >& AmplitudeBase::virtualInfos() { static map > > theVirtualInfos; return theVirtualInfos; } struct cfless { inline bool operator() (const vector& a, const vector& b) const { if ( a.size() > b.size() ) return true; return less >()(a,b); } }; size_t toColorFull(const set,cfless>& cs) { ostringstream cfids; cfids << "["; for ( set,cfless>::const_iterator l = cs.begin(); l != cs.end(); ++l ) { cfids << "{"; for ( vector::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& AmplitudeBase::amplitudeInfo() { const cPDVector& proc = mePartonData(); map >::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::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& cmap = amplitudeToColourMap(); #ifdef AMPVERBOSE // -------------------------------------------------------------------------------- generator()->log() << "generated amplitude info:\n\n" << flush; // -------------------------------------------------------------------------------- #endif for ( vector::iterator k = i->second.begin(); k != i->second.end(); ++k ) { #ifdef AMPVERBOSE // -------------------------------------------------------------------------------- k->print(generator()->log()); // -------------------------------------------------------------------------------- #endif map,cfless>,double> colourStructures; if ( !k->qqbarEmitted ) { vector 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 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,cfless> cStructure; cStructure.insert(ijLine); cStructure.insert(klLine); colourStructures.insert(make_pair(cStructure,1.)); } else { if ( k->ijLineEmissions.first > 0 ) { vector klLine; klLine.push_back(cmap[k->klLine.first]+1); klLine.push_back(cmap[k->klLine.second]+1); vector ijLineLeading1; ijLineLeading1.push_back(cmap[k->ijLine.first]+1); ijLineLeading1.push_back(cmap[k->ijLineEmissions.second]+1); vector ijLineLeading2; ijLineLeading2.push_back(cmap[k->ijLineEmissions.first]+1); ijLineLeading2.push_back(cmap[k->ijLine.second]+1); vector ijLineSubleading1; ijLineSubleading1.push_back(cmap[k->ijLine.first]+1); ijLineSubleading1.push_back(cmap[k->ijLine.second]+1); vector ijLineSubleading2; ijLineSubleading2.push_back(cmap[k->ijLineEmissions.first]+1); ijLineSubleading2.push_back(cmap[k->ijLineEmissions.second]+1); set,cfless> cStructureLeading; cStructureLeading.insert(klLine); cStructureLeading.insert(ijLineLeading1); cStructureLeading.insert(ijLineLeading2); set,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 ijLine; ijLine.push_back(cmap[k->ijLine.first]+1); ijLine.push_back(cmap[k->ijLine.second]+1); vector klLineLeading1; klLineLeading1.push_back(cmap[k->klLine.first]+1); klLineLeading1.push_back(cmap[k->klLineEmissions.second]+1); vector klLineLeading2; klLineLeading2.push_back(cmap[k->klLineEmissions.first]+1); klLineLeading2.push_back(cmap[k->klLine.second]+1); vector klLineSubleading1; klLineSubleading1.push_back(cmap[k->klLine.first]+1); klLineSubleading1.push_back(cmap[k->klLine.second]+1); vector klLineSubleading2; klLineSubleading2.push_back(cmap[k->klLineEmissions.first]+1); klLineSubleading2.push_back(cmap[k->klLineEmissions.second]+1); set,cfless> cStructureLeading; cStructureLeading.insert(ijLine); cStructureLeading.insert(klLineLeading1); cStructureLeading.insert(klLineLeading2); set,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,cfless>,double>::const_iterator c = colourStructures.begin(); c != colourStructures.end(); ++c ) { generator()->log() << "["; for ( set,cfless>::const_iterator l = c->first.begin(); l != c->first.end(); ++l ) { generator()->log() << "{"; for ( vector::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,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::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 >& AmplitudeBase::virtualInfo() const { const cPDVector& proc = mePartonData(); map > >::const_iterator i = virtualInfos().find(proc); if ( i != virtualInfos().end() ) return i->second; vector in; vector 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 > 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 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 > 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 > 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 > 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 > 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 > 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 alphaQED = SM().alphaEMMZ(); complex c2tw = complex(1.,0.) - SM().sin2ThetaW(); complex s2tw = SM().sin2ThetaW(); complex MW2c = double(sqr(MW/GeV)); complex MWc = sqrt(MW2c); complex MZ2c = double(sqr(MZ/GeV)); complex MZc = sqrt(MZ2c); if( complexMassScheme ){ MW2c = complex(sqr(MW/GeV),-MW*WWidth/GeV2); MWc = sqrt(MW2c); MZ2c = complex(sqr(MZ/GeV),-MZ*ZWidth/GeV2); MZc = sqrt(MZ2c); c2tw = MW2c/MZ2c; s2tw = complex(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 stw = sqrt(s2tw); complex ctw = sqrt(c2tw); complex 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 CUpperR = -gem*Q*stw/ctw; complex 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 CLowerR = -gem*Q*stw/ctw; complex CLowerL = gem*(I3-s2tw*Q)/(stw*ctw); //find NC topology and fill in couplings complex 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 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 describeHJetsAmplitudeBase("HJets::AmplitudeBase", "HwMatchboxBuiltin.so HJets.so"); void AmplitudeBase::Init() { static ClassDocumentation documentation ("Helicity amplitudes for 0 -> h + jets"); static Switch 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 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 interfaceKappaZ ("KappaZ", "The factor to rescale the HZZ coupling.", &AmplitudeBase::theKappaZ, 1.0, 0, 0, false, false, Interface::nolimits); static Parameter interfaceKappaW ("KappaW", "The factor to rescale the HWW coupling.", &AmplitudeBase::theKappaW, 1.0, 0, 0, false, false, Interface::nolimits); static Parameter 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 #\n" << "# Terrance M. Figy #\n" << "# Simon Platzer #\n" << "# Malin Sjodahl #\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 >& virtualInfo() const; /** * The virtuals info maps. */ static map > >& 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(); /** * 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 >& 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::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& helicities, Complex& largeN) { // get configurations we need to consider const vector& confs = amplitudeInfo(); // sum of all contributions Complex result = 0; largeN = 0.; for ( vector::const_iterator c = confs.begin(); c != confs.end(); ++c ) { map::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 >& 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 > 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 >& 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 >& 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 describeHJetsAmplitudehqqbarkkbar("HJets::Amplitudehqqbarkkbar", "HwMatchboxBuiltin.so HJets.so"); void Amplitudehqqbarkkbar::Init() { static ClassDocumentation 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::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& helicities, Complex& largeN) { // get configurations we need to consider const vector& confs = amplitudeInfo(); // sum of all contributions Complex result = 0; largeN = 0.; for ( vector::const_iterator c = confs.begin(); c != confs.end(); ++c ) { map::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 >& 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="<log() << "checkacc "<< log10 (abs(gacc))<<"\n"; } #ifdef AMPVERBOSE const cPDVector& proc = mePartonData(); map > 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]="<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 >& 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="<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 >& 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="<log() << "NLight= "<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(0,1))* (proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/abs((meMomenta()[0] + meMomenta()[1]).m2())); double cT1T2 = colourCorrelatedME2(pair(2,3))* (proc[2]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/abs((meMomenta()[2] + meMomenta()[3]).m2())); double cT1T3 = colourCorrelatedME2(pair(2,4))* (proc[2]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[2] + meMomenta()[4]).m2()))); double cT2T3 = colourCorrelatedME2(pair(3,4))* (proc[3]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[3] + meMomenta()[4]).m2()))); double cTaT1 = colourCorrelatedME2(pair(0,2))* (proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[0] + meMomenta()[2]).m2()))); double cTaT2 = colourCorrelatedME2(pair(0,3))* (proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[0] + meMomenta()[3]).m2()))); double cTaT3 = colourCorrelatedME2(pair(0,4))* (proc[0]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[0] + meMomenta()[4]).m2()))); double cTbT1 = colourCorrelatedME2(pair(1,2))* (proc[1]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[1] + meMomenta()[2]).m2()))); double cTbT2 = colourCorrelatedME2(pair(1,3))* (proc[1]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[1] + meMomenta()[3]).m2()))); double cTbT3 = colourCorrelatedME2(pair(1,4))* (proc[1]->id() == ParticleID::g ? CA1 : CF1)*log(scale*GeV2/(abs((meMomenta()[1] + meMomenta()[4]).m2()))); generator()->log() <<"CA1"<log()<<"res1/born="< 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 describeHJetsAmplitudehqqbarkkbarg("HJets::Amplitudehqqbarkkbarg", "HwMatchboxBuiltin.so HJets.so"); void Amplitudehqqbarkkbarg::Init() { static ClassDocumentation documentation ("Helicity amplitudes for 0 -> h q qbar k kbar g"); static Parameter interfaceGaugeThreshold ("GaugeThreshold", "The threshold for gauge check failure.", &Amplitudehqqbarkkbarg::theGaugeThreshold, 0.1, 0.0, 0, false, false, Interface::lowerlim); static Switch 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 `' 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 # # Terrance M. Figy # # Simon Platzer # # 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. # # # # 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 +