diff --git a/AUTHORS b/AUTHORS --- a/AUTHORS +++ b/AUTHORS @@ -0,0 +1,4 @@ +Francisco Campanario +Terrance M. Figy +Simon Platzer +Malin Sjodahl diff --git a/Amplitudes/HJetsProcessInfo.cc b/Amplitudes/HJetsProcessInfo.cc --- a/Amplitudes/HJetsProcessInfo.cc +++ b/Amplitudes/HJetsProcessInfo.cc @@ -1,626 +1,626 @@ // -*- C++ -*- // // HJetsProcessInfo.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. // #include "HJetsProcessInfo.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/PDT/EnumParticles.h" using namespace HJets; void AmplitudeInfo::print(ostream& os) const { if ( ijLineEmissions.first >= 0 && ijLineEmissions.second >= 0 ) { os << " " << ijLineEmissions.first << " " << ijLineEmissions.second << " \n" << " \\ / \n" << " \\ / \n" << " \\ / \n"; } if ( ijLineEmissions.first >= 0 && ijLineEmissions.second < 0 ) { os << " " << ijLineEmissions.first << " \n" << " | \n" << " | \n" << " | \n"; } if ( ijLineEmissions.first >= 0 || ijLineEmissions.second >= 0 ) { os << ijLine.first << "---<----0----<---" << ijLine.second << "\n"; } else { os << ijLine.first << "---<---------<---" << ijLine.second << "\n"; } os << " \\ " << boson->PDGName() << "\n" << " / \n" << " \\........0\n" << " / \n" << " \\ \n"; if ( klLineEmissions.first >= 0 || klLineEmissions.second >= 0 ) { os << klLine.first << "---<----0----<---" << klLine.second << "\n"; } else { os << klLine.first << "---<---------<---" << klLine.second << "\n"; } if ( klLineEmissions.first >= 0 && klLineEmissions.second >= 0 ) { os << " / \\ \n" << " / \\ \n" << " / \\ \n" << " " << klLineEmissions.first << " " << klLineEmissions.second << " \n"; } if ( klLineEmissions.first >= 0 && klLineEmissions.second < 0 ) { os << " | \n" << " | \n" << " | \n" << " " << klLineEmissions.first << " \n"; } os << "x " << fermionSign << "\n\n"; } bool isZ(long a, long b) { return a + b == 0 && abs(a) < 6; } bool isZPair(long a, long b, long c, long d) { return isZ(a,b) && isZ(c,d); } bool isWplus(long a, long b) { if ( abs(a)%2 != 0 ) swap(a,b); if ( a < 0 || b > 0 ) return false; return a == ParticleID::u && b == ParticleID::dbar || a == ParticleID::c && b == ParticleID::sbar; } bool isWminus(long a, long b) { return isWplus(-a,-b); } bool isWPair(long a, long b, long c, long d) { return (isWplus(a,b) && isWminus(c,d)) || (isWminus(a,b) && isWplus(c,d)); } vector HJetsProcessInfo::configurations(const cPDVector& proc, const StandardModelBase& sm, bool complexMassScheme, double kappaZ, double kappaW) { /* cerr << "check if can handle "; for ( cPDVector::const_iterator p = proc.begin(); p != proc.end(); ++p ) cerr << (**p).PDGName() << " "; cerr << "\n" << flush; */ vector res; cPDVector xproc = proc; if ( xproc[0]->CC() ) xproc[0] = xproc[0]->CC(); if ( xproc[1]->CC() ) xproc[1] = xproc[1]->CC(); cPDVector::iterator h = xproc.begin(); for ( ; h != xproc.end(); ++h ) if ( (**h).id() == ParticleID::h0 ) break; if ( h == xproc.end() ) return res; map >, vector > > coreAndEmissions; vector > proto; for ( size_t k = 0; k < xproc.size(); ++k ) { if ( xproc[k]->id() == ParticleID::h0 ) continue; proto.push_back(pair(k,xproc[k])); } if ( proto.size() == 4 ) { coreAndEmissions.insert(make_pair(proto,vector >())); } if ( proto.size() > 4 ) { bool gotGlue = false; for ( vector >::const_iterator p = proto.begin(); p != proto.end(); ++p ) if ( p->second->id() == ParticleID::g ) { gotGlue = true; break; } if ( gotGlue ) { vector > core; vector > emissions; for ( vector >::const_iterator p = proto.begin(); p != proto.end(); ++p ) if ( p->second->id() == ParticleID::g ) { emissions.push_back(*p); } else { core.push_back(*p); } if ( core.size() != 4 || emissions.size() > 2 ) return res; coreAndEmissions.insert(make_pair(core,emissions)); } else { if ( proto.size() != 6 ) return res; for ( vector >::const_iterator i = proto.begin(); i != proto.end(); ++i ) { vector >::const_iterator j = i; ++j; for ( ; j != proto.end(); ++j ) { if ( abs(i->second->id()) < 6 && abs(j->second->id()) < 6 && i->second->id() + j->second->id() == 0 ) { vector > core; vector > emissions; for ( vector >::const_iterator k = proto.begin(); k != proto.end(); ++k ) if ( k != i && k != j ) core.push_back(*k); emissions.push_back(*i); emissions.push_back(*j); coreAndEmissions.insert(make_pair(core,emissions)); } } } } } - Energy MZ = sm.getParticleData(ParticleID::Z0)->mass(); - Energy GZ = sm.getParticleData(ParticleID::Z0)->width(); - Energy MW = sm.getParticleData(ParticleID::Wplus)->mass(); - Energy GW = sm.getParticleData(ParticleID::Wplus)->width(); + Energy MZ = sm.getParticleData(ParticleID::Z0)->hardProcessMass(); + Energy GZ = sm.getParticleData(ParticleID::Z0)->hardProcessWidth(); + Energy MW = sm.getParticleData(ParticleID::Wplus)->hardProcessMass(); + Energy GW = sm.getParticleData(ParticleID::Wplus)->hardProcessWidth(); complex alphaQED = sm.alphaEMMZ(); complex c2tw = complex(1.,0.) - sm.sin2ThetaW(); complex s2tw = sm.sin2ThetaW(); if ( complexMassScheme ) { c2tw = complex(sqr(MW),-MW*GW)/complex(sqr(MZ),-MZ*GZ); s2tw = complex(1.,0.) - c2tw; alphaQED = sqrt(2.)*sm.fermiConstant()*complex(sqr(MW),-MW*GW)*s2tw/Constants::pi; } complex stw = sqrt(s2tw); complex ctw = sqrt(c2tw); complex gem = sqrt(4.*Constants::pi*alphaQED); map >, vector > > coresSorted; for ( map >, vector > >::const_iterator c = coreAndEmissions.begin(); c != coreAndEmissions.end(); ++c ) { bool valid = true; vector > coreSorted(4,pair(-1,cPDPtr())); for ( size_t k = 0; k < 4; ++k ) { if ( c->first[k].second->id() > 0 ) if ( coreSorted[0].first < 0 ) { coreSorted[0].first = c->first[k].first; coreSorted[0].second = c->first[k].second; } else if ( coreSorted[2].first < 0 ) { coreSorted[2].first = c->first[k].first; coreSorted[2].second = c->first[k].second; } else { valid = false; } if ( c->first[k].second->id() < 0 ) if ( coreSorted[1].first < 0 ) { coreSorted[1].first = c->first[k].first; coreSorted[1].second = c->first[k].second; } else if ( coreSorted[3].first < 0 ) { coreSorted[3].first = c->first[k].first; coreSorted[3].second = c->first[k].second; } else { valid = false; } } if ( !valid ) continue; if ( coreSorted[0].first > coreSorted[2].first ) { swap(coreSorted[0],coreSorted[2]); swap(coreSorted[1],coreSorted[3]); } vector > emissionsSorted = c->second; if ( emissionsSorted.size() == 2 ) { if ( emissionsSorted[0].second->id() < emissionsSorted[1].second->id() ) swap(emissionsSorted[0],emissionsSorted[1]); if ( emissionsSorted[0].second->id() == emissionsSorted[1].second->id() ) { assert(emissionsSorted[0].second->id() == ParticleID::g); if ( emissionsSorted[0].first > emissionsSorted[1].first ) swap(emissionsSorted[0],emissionsSorted[1]); } } coresSorted.insert(make_pair(coreSorted,emissionsSorted)); } coreAndEmissions = coresSorted; /* cerr << "cores sorted are:\n"; for ( map >, vector > >::const_iterator c = coreAndEmissions.begin(); c != coreAndEmissions.end(); ++c ) { for ( vector >::const_iterator p = c->first.begin(); p != c->first.end(); ++p ) { cerr << p->first << " " << p->second->PDGName() << " "; } cerr << " -- "; for ( vector >::const_iterator p = c->second.begin(); p != c->second.end(); ++p ) { cerr << p->first << " " << p->second->PDGName() << " "; } cerr << "\n"; } */ for ( map >, vector > >::const_iterator c = coreAndEmissions.begin(); c != coreAndEmissions.end(); ++c ) { vector cores; AmplitudeInfo core; core.ijLineEmissions.first = -1; core.ijLineEmissions.second = -1; core.klLineEmissions.first = -1; core.klLineEmissions.second = -1; core.qqbarEmitted = false; core.fermionSign = 1.; if ( c->second.size() == 2 ) if ( abs(c->second[0].second->id()) < 6 ) core.qqbarEmitted = true; if ( isZPair(c->first[0].second->id(),c->first[1].second->id(), c->first[2].second->id(),c->first[3].second->id()) || isZPair(c->first[0].second->id(),c->first[3].second->id(), c->first[1].second->id(),c->first[2].second->id()) ) { double Q; double I3; core.boson = sm.getParticleData(ParticleID::Z0); core.bosonMass = MZ; core.bosonWidth = GZ; core.higgsCoupling = kappaZ*gem*MZ/(stw*ctw); if( complexMassScheme){ complex MZ2c = complex(sqr(MZ),-MZ*GZ); complex MZc = GeV*sqrt(MZ2c/GeV2); core.higgsCoupling = kappaZ*gem*MZc/(stw*ctw); } Q = c->first[0].second->iCharge()/3.; I3 = c->first[0].second->id() % 2 == 0 ? 0.5 : -0.5; core.ijRightCoupling = -gem*Q*stw/ctw; core.ijLeftCoupling = gem*(I3-s2tw*Q)/(stw*ctw); Q = c->first[2].second->iCharge()/3.; I3 = c->first[2].second->id() % 2 == 0 ? 0.5 : -0.5; core.klRightCoupling = -gem*Q*stw/ctw; core.klLeftCoupling = gem*(I3-s2tw*Q)/(stw*ctw); if ( isZPair(c->first[0].second->id(),c->first[1].second->id(), c->first[2].second->id(),c->first[3].second->id()) ) { core.ijLine.first = c->first[0].first; core.ijLine.second = c->first[1].first; core.klLine.first = c->first[2].first; core.klLine.second = c->first[3].first; cores.push_back(core); } if ( isZPair(c->first[0].second->id(),c->first[3].second->id(), c->first[1].second->id(),c->first[2].second->id()) ) { core.ijLine.first = c->first[0].first; core.ijLine.second = c->first[3].first; core.klLine.first = c->first[2].first; core.klLine.second = c->first[1].first; cores.push_back(core); } } if ( isWPair(c->first[0].second->id(),c->first[1].second->id(), c->first[2].second->id(),c->first[3].second->id()) || isWPair(c->first[0].second->id(),c->first[3].second->id(), c->first[1].second->id(),c->first[2].second->id()) ) { core.bosonMass = MW; core.bosonWidth = GW; core.higgsCoupling = kappaW*gem*MW/stw; if( complexMassScheme){ complex MW2c = complex(sqr(MW),-MW*GW); complex MWc = GeV*sqrt(MW2c/GeV2); core.higgsCoupling = kappaW*gem*MWc/stw; } core.ijLeftCoupling = gem/(sqrt(2.)*stw); core.ijRightCoupling = 0.; core.klLeftCoupling = gem/(sqrt(2.)*stw); core.klRightCoupling = 0.; if ( isWPair(c->first[0].second->id(),c->first[1].second->id(), c->first[2].second->id(),c->first[3].second->id()) ) { if ( isWplus(c->first[0].second->id(),c->first[1].second->id()) ) { core.boson = sm.getParticleData(ParticleID::Wplus); } else { core.boson = sm.getParticleData(ParticleID::Wminus); } core.ijLine.first = c->first[0].first; core.ijLine.second = c->first[1].first; core.klLine.first = c->first[2].first; core.klLine.second = c->first[3].first; cores.push_back(core); } if ( isWPair(c->first[0].second->id(),c->first[3].second->id(), c->first[1].second->id(),c->first[2].second->id()) ) { if ( isWplus(c->first[0].second->id(),c->first[3].second->id()) ) { core.boson = sm.getParticleData(ParticleID::Wplus); } else { core.boson = sm.getParticleData(ParticleID::Wminus); } core.ijLine.first = c->first[0].first; core.ijLine.second = c->first[3].first; core.klLine.first = c->first[2].first; core.klLine.second = c->first[1].first; cores.push_back(core); } } for ( vector::const_iterator conf = cores.begin(); conf != cores.end(); ++conf ) { AmplitudeInfo full = *conf; if ( c->second.size() == 0 ) { res.push_back(full); } if ( c->second.size() == 1 ) if ( c->second[0].second->id() == ParticleID::g ) { full.ijLineEmissions.first = c->second[0].first; res.push_back(full); full.ijLineEmissions.first = -1; full.klLineEmissions.first = c->second[0].first; res.push_back(full); full.klLineEmissions.first = -1; } if ( c->second.size() == 2 ) if ( c->second[0].second->id() == ParticleID::g && c->second[1].second->id() == ParticleID::g ) { full.ijLineEmissions.first = c->second[0].first; full.ijLineEmissions.second = c->second[1].first; res.push_back(full); swap(full.ijLineEmissions.first,full.ijLineEmissions.second); res.push_back(full); full.ijLineEmissions.first = -1; full.ijLineEmissions.second = -1; full.klLineEmissions.first = c->second[0].first; full.klLineEmissions.second = c->second[1].first; res.push_back(full); swap(full.klLineEmissions.first,full.klLineEmissions.second); res.push_back(full); full.klLineEmissions.first = -1; full.klLineEmissions.second = -1; full.ijLineEmissions.first = c->second[0].first; full.klLineEmissions.first = c->second[1].first; res.push_back(full); swap(full.ijLineEmissions.first,full.klLineEmissions.first); res.push_back(full); full.ijLineEmissions.first = -1; full.klLineEmissions.first = -1; } if ( c->second.size() == 2 ) if ( abs(c->second[0].second->id()) < 6 && c->second[0].second->id() + c->second[1].second->id() == 0 ) { full.ijLineEmissions.first = c->second[0].first; full.ijLineEmissions.second = c->second[1].first; res.push_back(full); full.ijLineEmissions.first = -1; full.ijLineEmissions.second = -1; full.klLineEmissions.first = c->second[0].first; full.klLineEmissions.second = c->second[1].first; res.push_back(full); full.klLineEmissions.first = -1; full.klLineEmissions.second = -1; } } } return res; } bool HJetsProcessInfo::canHandle(const PDVector& proc, const StandardModelBase& sm) { cPDVector cproc; copy(proc.begin(), proc.end(), back_inserter(cproc)); vector confs = configurations(cproc,sm); return !confs.empty(); } vector HJetsProcessInfo::getConfigurations(const cPDVector& proc, const vector& crossed2proc, const StandardModelBase& sm, bool complexMassScheme, double kappaZ, double kappaW) { vector procConfigs = configurations(proc,sm,complexMassScheme,kappaZ,kappaW); /* cerr << "configurations before crossing\n\n"; for ( vector::const_iterator c = procConfigs.begin(); c != procConfigs.end(); ++c ) { c->print(cerr); cerr << "\n"; } */ vector ampConfigs; map proc2crossed; for ( size_t crossedx = 0; crossedx < crossed2proc.size(); ++crossedx ) proc2crossed[crossed2proc[crossedx]] = crossedx; proc2crossed[-1] = -1; for ( vector::const_iterator c = procConfigs.begin(); c != procConfigs.end(); ++c ) { AmplitudeInfo crossed = *c; crossed.ijLine.first = proc2crossed[crossed.ijLine.first]; crossed.ijLine.second = proc2crossed[crossed.ijLine.second]; crossed.klLine.first = proc2crossed[crossed.klLine.first]; crossed.klLine.second = proc2crossed[crossed.klLine.second]; crossed.ijLineEmissions.first = proc2crossed[crossed.ijLineEmissions.first]; crossed.ijLineEmissions.second = proc2crossed[crossed.ijLineEmissions.second]; crossed.klLineEmissions.first = proc2crossed[crossed.klLineEmissions.first]; crossed.klLineEmissions.second = proc2crossed[crossed.klLineEmissions.second]; // work out the fermion signs set > ordering; ordering.insert(crossed.ijLine); ordering.insert(crossed.klLine); if ( crossed.qqbarEmitted ) { if ( crossed.ijLineEmissions.first >= 0 && crossed.ijLineEmissions.second >= 0 ) ordering.insert(crossed.ijLineEmissions); else if ( crossed.klLineEmissions.first >= 0 && crossed.klLineEmissions.second >= 0 ) ordering.insert(crossed.klLineEmissions); else assert(false); } set >::iterator first = ordering.begin(); set >::iterator second = first; ++second; set >::iterator third = ordering.end(); if ( crossed.qqbarEmitted ) { third = second; ++third; } if ( !crossed.qqbarEmitted ) { if ( *first == pair(1,2) && *second == pair(3,4) ) { crossed.fermionSign = 1.; } else if ( *first == pair(1,4) && *second == pair(3,2) ) { crossed.fermionSign = -1.; } else assert(false); } else { if ( *first == pair(1,2) && *second == pair(3,4) && *third == pair(5,6) ) { crossed.fermionSign = 1.; } else if ( *first == pair(1,4) && *second == pair(3,2) && *third == pair(5,6) ) { crossed.fermionSign = -1.; } else if ( *first == pair(1,6) && *second == pair(3,4) && *third == pair(5,2) ) { crossed.fermionSign = -1.; } else if ( *first == pair(1,2) && *second == pair(3,6) && *third == pair(5,4) ) { crossed.fermionSign = -1.; } else if ( *first == pair(1,4) && *second == pair(3,6) && *third == pair(5,2) ) { crossed.fermionSign = 1.; } else if ( *first == pair(1,6) && *second == pair(3,2) && *third == pair(5,4) ) { crossed.fermionSign = 1.; } else assert(false); } ampConfigs.push_back(crossed); } return ampConfigs; } diff --git a/Amplitudes/HJetsScaleChoice.cc b/Amplitudes/HJetsScaleChoice.cc deleted file mode 100644 --- a/Amplitudes/HJetsScaleChoice.cc +++ /dev/null @@ -1,193 +0,0 @@ -// -*- C++ -*- -// -// HJetsScaleChoice.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. -// -// -// This is the implementation of the non-inlined, non-templated member -// functions of the HJetsScaleChoice class. -// - -#include "HJetsScaleChoice.h" -#include "ThePEG/Interface/ClassDocumentation.h" -#include "ThePEG/Interface/Switch.h" -#include "ThePEG/Interface/Parameter.h" -#include "ThePEG/Interface/Reference.h" -#include "ThePEG/EventRecord/Particle.h" -#include "ThePEG/Repository/UseRandom.h" -#include "ThePEG/Repository/EventGenerator.h" -#include "ThePEG/Utilities/DescribeClass.h" - - -#include "ThePEG/Persistency/PersistentOStream.h" -#include "ThePEG/Persistency/PersistentIStream.h" - -using namespace HJets; - -HJetsScaleChoice::HJetsScaleChoice() - : theHtFactor(1.0), theJetsOnly(true), - theIncludeHiggs(false), theDoAverage(false) {} - -HJetsScaleChoice::~HJetsScaleChoice() {} - -IBPtr HJetsScaleChoice::clone() const { - return new_ptr(*this); -} - -IBPtr HJetsScaleChoice::fullclone() const { - return new_ptr(*this); -} - -Energy2 HJetsScaleChoice::renormalizationScale() const { - - tcPDVector pd (mePartonData().begin() + 2, mePartonData().end()); - vector p (meMomenta().begin() + 2, meMomenta().end()); - tcPDPtr t1 = mePartonData()[0]; - tcPDPtr t2 = mePartonData()[1]; - tcCutsPtr cuts = lastCutsPtr(); - - - theJetFinder->cluster(pd, p, cuts, t1, t2); - - - Energy htjets = ZERO; - Energy mthiggs = ZERO; - double yboost = 0.0; - int found = 0; - tcPDVector::const_iterator itpd = pd.begin(); - - for (vector::const_iterator itp = p.begin() ; - itp != p.end(); ++itp, ++itpd ) { - if ( theJetsOnly ) { - if ( (**itpd).coloured()){ - yboost += (*itp).rapidity(); - found++; - } - } - } - //cout << "found " << found << flush; - //cout << "yboost " << yboost << flush; - // IS THERE A BETTER WAY TO DO THIS? - tcPDVector::const_iterator itpd1 = pd.begin(); - - for (vector::const_iterator itp1 = p.begin() ; - itp1 != p.end(); ++itp1, ++itpd1 ) { - - if ( theJetsOnly ) { - if ( (**itpd1).coloured() ){ - htjets += (*itp1).perp()*exp(abs((*itp1).rapidity()-yboost/found)); - } - } - else { - htjets += (*itp1).perp(); - found++; - } - // transverse mass of onshell higgs boson - if (theIncludeHiggs ){ - if( (**itpd1).id() == 25 ){ - mthiggs += (*itp1).mt(); - // cout << "mass of higgs=" << (*itp).m()/GeV << "\n" << flush; - } - } - - } - if (theDoAverage && found !=0) return sqr((mthiggs + htjets/found)/theHtFactor); - return sqr((mthiggs+htjets)/theHtFactor); - -} - -Energy2 HJetsScaleChoice::factorizationScale() const { - return renormalizationScale(); -} - -// If needed, insert default implementations of virtual function defined -// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). - - -void HJetsScaleChoice::persistentOutput(PersistentOStream & os) const { - os << theJetFinder << theHtFactor << theJetsOnly << theDoAverage; -} - -void HJetsScaleChoice::persistentInput(PersistentIStream & is, int) { - is >> theJetFinder >> theHtFactor >> theJetsOnly << theDoAverage; -} - - -// *** 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 - describeHJetsScaleChoice("HJets::HJetsScaleChoice", "HJets.so"); - -void HJetsScaleChoice::Init() { - - static ClassDocumentation documentation - ("HJetsScaleChoice implements scale choices."); - - static Parameter interfaceHtFactor - ("HtFactor", - "The factor dividing HT", - &HJetsScaleChoice::theHtFactor, 1.0, 0.01, 10, true, false, true ); - - - static Reference interfaceJetFinder - ("JetFinder", - "A reference to the jet finder.", - &HJetsScaleChoice::theJetFinder, false, false, true, false, false); - - static Switch interfaceJetsOnly - ("JetsOnly", - "The mode to use.", - &HJetsScaleChoice::theJetsOnly, true, false, false); - static SwitchOption interfaceJetsOnlyTrue - (interfaceJetsOnly, - "True", - "Only include jets.", - true); - static SwitchOption interfaceJetsOnlyFalse - (interfaceJetsOnly, - "False", - "Include all outgoing particles in the matrix element.", - false); - - - static Switch interfaceIncludeHiggs - ("IncludeHiggs", - "Include Higgs transverse mass", - &HJetsScaleChoice::theIncludeHiggs, false, false, false); - - static SwitchOption interfaceIncludeHiggsTrue - (interfaceIncludeHiggs, - "True", - "Include Higgs transverse mass", - true); - static SwitchOption interfaceIncludeHiggsFalse - (interfaceIncludeHiggs, - "False", - "do not include higgs travsverse mass", - false); - - static Switch interfaceDoAverage - ("DoAverage", - "The mode to use.", - &HJetsScaleChoice::theDoAverage, false, false, false); - static SwitchOption interfaceDoAverageTrue - (interfaceDoAverage, - "True", - "Average over number of jets/particles", - true); - static SwitchOption interfaceDoAverageFalse - (interfaceDoAverage, - "False", - "Do not average.", - false); - - -} - diff --git a/Amplitudes/HJetsScaleChoice.h b/Amplitudes/HJetsScaleChoice.h deleted file mode 100644 --- a/Amplitudes/HJetsScaleChoice.h +++ /dev/null @@ -1,148 +0,0 @@ -// -*- C++ -*- -// -// HJetsScaleChoice.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_HJetsScaleChoice_H -#define HJets_HJetsScaleChoice_H -// -// This is the declaration of the HJetsScaleChoice class. -// - -#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h" -#include "ThePEG/Cuts/JetFinder.h" - -namespace HJets { - -using namespace ThePEG; -using namespace Herwig; - -/** - * \ingroup Matchbox - * \author Simon Platzer - * - * \brief HJetsScaleChoice implements scale choices related to transverse momenta. - * - */ -class HJetsScaleChoice: public MatchboxScaleChoice { - -public: - - /** @name Standard constructors and destructors. */ - //@{ - /** - * The default constructor. - */ - HJetsScaleChoice(); - - /** - * The destructor. - */ - virtual ~HJetsScaleChoice(); - //@} - -public: - - /** - * Return the renormalization scale. This default version returns - * shat. - */ - virtual Energy2 renormalizationScale() const; - - /** - * Return the factorization scale. This default version returns - * shat. - */ - virtual Energy2 factorizationScale() const; - -public: - - /** @name Functions used by the persistent I/O system. */ - //@{ - /** - * Function used to write out object persistently. - * @param os the persistent output stream written to. - */ - void persistentOutput(PersistentOStream & os) const; - - /** - * Function used to read in object persistently. - * @param is the persistent input stream read from. - * @param version the version number of the object when written. - */ - void persistentInput(PersistentIStream & is, int version); - //@} - - /** - * The standard Init function used to initialize the interfaces. - * Called exactly once for each class by the class description system - * before the main function starts or - * when this class is dynamically loaded. - */ - static void Init(); - -protected: - - /** @name Clone Methods. */ - //@{ - /** - * Make a simple clone of this object. - * @return a pointer to the new object. - */ - virtual IBPtr clone() const; - - /** Make a clone of this object, possibly modifying the cloned object - * to make it sane. - * @return a pointer to the new object. - */ - virtual IBPtr fullclone() const; - //@} - - -// If needed, insert declarations of virtual function defined in the -// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). - - -private: - - /** - * Reference to the jet finder - */ - Ptr::ptr theJetFinder; - - /** - * Factor dividing HT (i.e number of jets) - */ - double theHtFactor; - - /** - * Choose to use only jets or to include all other particles, too. - */ - bool theJetsOnly; - - /** - * Choose to use only jets or to include all other particles, too. - */ - bool theIncludeHiggs; - - /** - * Set to true to return the average scalar transverse momentum squared - */ - bool theDoAverage; - - - /** - * The assignment operator is private and must never be called. - * In fact, it should not even be implemented. - */ - HJetsScaleChoice & operator=(const HJetsScaleChoice &); - -}; - -} - -#endif /* HJets_HJetsScaleChoice_H */ diff --git a/Amplitudes/Makefile.am b/Amplitudes/Makefile.am --- a/Amplitudes/Makefile.am +++ b/Amplitudes/Makefile.am @@ -1,20 +1,18 @@ pkglib_LTLIBRARIES = HJets.la HJets_la_LDFLAGS = -module -version-info 11:1:0 HJets_la_LIBADD = $(top_builddir)/Virtuals/libHJetsVirtuals.la HJets_la_SOURCES = \ AmplitudeBase.h \ Amplitudehqqbarkkbarrrbar.h \ Amplitudehqqbarkkbargg.h \ Amplitudehqqbarkkbarg.h \ Amplitudehqqbarkkbar.h \ HJetsProcessInfo.h \ AmplitudeBase.cc \ Amplitudehqqbarkkbar.cc \ Amplitudehqqbarkkbarg.cc \ Amplitudehqqbarkkbargg.cc \ Amplitudehqqbarkkbarrrbar.cc \ -HJetsProcessInfo.cc \ -HJetsScaleChoice.h \ -HJetsScaleChoice.cc +HJetsProcessInfo.cc diff --git a/Analysis/HJetsAnalysis2.cc b/Analysis/HJetsAnalysis2.cc deleted file mode 100644 --- a/Analysis/HJetsAnalysis2.cc +++ /dev/null @@ -1,782 +0,0 @@ -// -*- C++ -*- -// -// This is the implementation of the non-inlined, non-templated member -// functions of the HJetsAnalysis2 class. -// - -#include "HJetsAnalysis2.h" -#include "ThePEG/Interface/ClassDocumentation.h" -#include "ThePEG/Interface/Reference.h" -#include "ThePEG/Interface/RefVector.h" -#include "ThePEG/Interface/Switch.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" - -#include "ThePEG/Handlers/EventHandler.h" -#include "ThePEG/Handlers/StandardEventHandler.h" -#include "Herwig++/Sampling/GeneralSampler.h" - -#include "Herwig++/Utilities/XML/ElementIO.h" - -using namespace HJets; -using namespace Herwig; -using namespace Statistics; - -HJetsAnalysis2::HJetsAnalysis2() - : theFullEvent(false), theBinFactor(1.0), theFuzzy(true), theEnergyCutWidth(1.0*GeV), theRapidityCutWidth(0.1),theRapidityGap(0.0), - theopposite_dir(false), thenjets_min(2), thehiggs_ingap(false), themjj_min(0.0*GeV),thehiggsPt_min(0.0*GeV), thehiggsY_max(1000.0) {} - -HJetsAnalysis2::~HJetsAnalysis2() {} - -#ifndef LWH_AIAnalysisFactory_H -#ifndef LWH -#define LWH ThePEGLWH -#endif -#include "ThePEG/Analysis/LWH/AnalysisFactory.h" -#endif - -inline double deltaPhi(const LorentzMomentum& a, - const LorentzMomentum& b){ - double phi1 = a.phi(); - double phi2 = b.phi(); - double diff=phi1-phi2; - if(diff<-Constants::pi){ - diff+=(2.0*Constants::pi); - } - else if (diff>Constants::pi){ - diff-=(2.0*Constants::pi); - } - return diff; -} - -inline double deltaY(const LorentzMomentum& a, - const LorentzMomentum& b){ - return abs(a.rapidity()-b.rapidity()); -} - -inline double deltaR(const LorentzMomentum& a, - const LorentzMomentum& b){ - return sqrt(sqr(deltaPhi(a,b))+sqr(deltaY(a,b))); -} - -inline double mass(const LorentzMomentum& a, - const LorentzMomentum& b){ - return (a+b).m()/GeV; -} - -inline double mass(const LorentzMomentum& a, - const LorentzMomentum& b, const LorentzMomentum& c){ - return (a+b+c).m()/GeV; -} - -void HJetsAnalysis2::analyze(tEventPtr event, long ieve, int loop, int state) { - AnalysisHandler::analyze(event, ieve, loop, state); - - - if ( !theFullEvent ) { - - tSubProPtr sub = event->primarySubProcess(); - Ptr::tptr grp = - dynamic_ptr_cast::tptr>(sub); - - analyze(sub->outgoing(),ieve,event->weight()*sub->groupWeight()); - - if ( grp ) { - - for ( SubProcessVector::const_iterator s = grp->dependent().begin(); - s != grp->dependent().end(); ++s ) { - analyze((**s).outgoing(),ieve,event->weight()*(**s).groupWeight()); - } - - } - - } else { - - ParticleVector fs; - event->getFinalState(fs); - analyze(fs,ieve,event->weight()); - - } - -} - -struct SortPt { - - inline bool operator()(const LorentzMomentum& a, - const LorentzMomentum& b) const { - return a.perp() > b.perp(); - } - -}; - - -double step(double r) { - if ( r < -0.5 ) - return 0.0; - if ( r > 0.5 ) - return 1.0; - return r+0.5; -} - -bool HJetsAnalysis2::lessThanEnergy(Energy a, Energy b, double& weight) const { - if ( !fuzzy() ) { - if ( a < b ) { - weight = 1.0; - return true; - } - weight = 0.0; - return false; - } - double w = step((b-a)/theEnergyCutWidth); - if ( w == 0.0 ) { - weight = 0.0; - return false; - } - weight *= w; - return true; -} - -bool HJetsAnalysis2::lessThanRapidity(double a, double b, double& weight) const { - if ( !fuzzy() ) { - if ( a < b ) { - weight = 1.0; - return true; - } - weight = 0.0; - return false; - } - double w = step((b-a)/theRapidityCutWidth); - if ( w == 0.0 ) { - weight = 0.0; - return false; - } - weight *= w; - return true; -} - - -bool HJetsAnalysis2::selectionCuts(const LorentzMomentum& h, - const vector& jets, double& weight){ - - - double theCutWeight(1.0); - - // higgs pt min and y_max - //std::cout << thehiggsPt_min/GeV << endl; - if ( !(lessThanRapidity(abs(h.rapidity()),thehiggsY_max,theCutWeight))){ - theCutWeight = 0.0; - return false; - } - - //std::cout << thehiggsPt_min/GeV << endl; - if ( !(lessThanEnergy(thehiggsPt_min,h.perp(),theCutWeight))){ - theCutWeight = 0.0; - return false; - } - - if(jets.size() > 1){ - // opposite detectors - double py = jets[0].rapidity()*jets[1].rapidity(); - - if (theopposite_dir){ - if( !(lessThanRapidity(py,0.0,theCutWeight)) ){ - theCutWeight=0.0; - return false; - } - } - // rapidity gap - if( !lessThanRapidity(theRapidityGap,deltaY(jets[0],jets[1]),theCutWeight)){ - theCutWeight = 0.0; - return false; - } - // invariant mass cut - - Energy massj1j2; - massj1j2 = (jets[0]+jets[1]).m(); - //std::cout << massj1j2/GeV << endl; - if( !lessThanEnergy(themjj_min,massj1j2,theCutWeight)){ - theCutWeight = 0.0; - return false; - } - - // higgs in the gap - - if(thehiggs_ingap){ - double minjrap = min(jets[0].rapidity(),jets[1].rapidity()); - double maxjrap = max(jets[0].rapidity(),jets[1].rapidity()); - if(!(lessThanRapidity(h.rapidity(),maxjrap,theCutWeight) && lessThanRapidity(minjrap,h.rapidity(),theCutWeight))){ - theCutWeight = 0.0; - return false; - } - } - - } - // should be a rescaling - weight=theCutWeight*weight; - return true; -} - -void HJetsAnalysis2::analyze(const ParticleVector& out, long ieve, double weight) { - - if ( weight == 0.0 ) - return; - - tcPDVector outType; - vector outMomenta; - - for ( ParticleVector::const_iterator p = out.begin(); - p != out.end(); ++p ) { - outType.push_back((**p).dataPtr()); - outMomenta.push_back((**p).momentum()); - } - - theJetFinder->cluster(outType, outMomenta, - tcCutsPtr(), tcPDPtr(),tcPDPtr()); - - for ( vector::ptr>::iterator j = - theJetRegions.begin(); j != theJetRegions.end(); ++j ) - (**j).reset(); - - LorentzMomentum higgs; - vector allJets; - - bool gotHiggs = false; - - tcPDVector::const_iterator pd = outType.begin(); - vector::const_iterator p = outMomenta.begin(); - for ( ; pd != outType.end(); ++pd, ++p ) { - if ( (**pd).id() == 25 ) { - gotHiggs = true; - higgs = *p; - } else if ( (**pd).coloured() ) { - allJets.push_back(*p); - } - } - - assert(gotHiggs); - - sort(allJets.begin(),allJets.end(),SortPt()); - // @TODO - // perhaps set a PT min threshold and max abs rapidity ABSY_MAX - // - vector selectedJets20; - for ( size_t k = 0; k < allJets.size(); ++k ) { - if (allJets[k].perp()/GeV > 20.0 ){ - selectedJets20.push_back(allJets[k]); - } - } - - sort(selectedJets20.begin(),selectedJets20.end(),SortPt()); - // - //std::cout << "number of jets from selected jets 20 " << selectedJets20.size() << endl; - // Should we only use this for plain NLO runs - size_t matchedJets = 0; - vector selectedJets; - - for ( size_t k = 0; k < allJets.size(); ++k ) { - for ( vector::ptr>::const_iterator r = theJetRegions.begin(); - r != theJetRegions.end(); ++r ) { - if ( (**r).matches(tCutsPtr(),k+1,allJets[k])) { - ++matchedJets; - selectedJets.push_back(allJets[k]); - break; - } - } - - } - // apply jets cuts according to the jet regions - - double cutWeight = 1.0; - bool pass = true; - bool pass_jetcuts = true; - - for ( vector::ptr>::const_iterator r = theJetRegions.begin(); - r != theJetRegions.end(); ++r ) { - pass &= (**r).didMatch(); - cutWeight *= (**r).cutWeight(); - if ( !pass ) { - pass_jetcuts = false; - cutWeight = 0.0; - break; - } - } - weight=cutWeight*weight; - - if ( matchedJets < theJetRegions.size() ) { - if ( !theFullEvent ) - generator()->log() << "ATTENTION less than number of expected jets!\n"; - return; - } - - sort(selectedJets.begin(),selectedJets.end(),SortPt()); - // - // for plain NLO, we do not want these cuts - // for NLO + PS, we do want these cuts - // For plain NLO, we do not use these cuts for the moment. - // @TODO: Looking for a better suggestion. - // - - bool cutsOK = false; - // must pass jet cuts - if (pass_jetcuts){ - cutsOK = selectionCuts(higgs,selectedJets,weight); - } - if(cutsOK){ - analyze(higgs,allJets,ieve,weight); - - } - - // -} - -void HJetsAnalysis2::analyze(const LorentzMomentum& h, - const vector& jets, - long xid, double weight) { - - const double pi = 3.14159265; - - const double massWidth = 4000.0/64.0*theBinFactor; - const double rapWidth1 = 60.0/64.0*theBinFactor; - const double rapWidth = 10.0/32.0*theBinFactor; - const double etaWidth = 10.0/32.0*theBinFactor; - const double ptWidth = 400.0/40.0*theBinFactor; - const double ptHiggsWidth = 400.0/40.0*theBinFactor; - const double phiWidth = 2.*pi/32.0*theBinFactor; - const double ysWidth = 4.0/32.0*theBinFactor; - const double jetMassWidth = 100.0/20.0*theBinFactor; - const double HTWidth = 1000.0/32.0*theBinFactor; - const double higgsDijetDeltaPhiWidth = pi/32.0*theBinFactor; - const double higgsDijetpTWidth=100/32.0*theBinFactor; - const double threeJetMassWidth=2000/32.0*theBinFactor; - const double jetSumDeltaPhiWidth=2.*pi/32.0*theBinFactor; - - // const double massWidth = 0.0; - // const double rapWidth = 0.0; - // const double etaWidth = 0.0; - // const double ptWidth = 0.0; - // const double ptHiggsWidth = 0.0; - // const double phiWidth = 0.0; - // const double ysWidth = 0.0; - // const double jetMassWidth = 0.0; - // const double HTWidth = 0.0; - // const double higgsDijetDeltaPhiWidth = 0.0; - // const double higgsDijetpTWidth=0.0; - // const double threeJetMassWidth=0.0; - // const double jetSumDeltaPhiWidth=0.0; - - - Ptr::tptr seh = - dynamic_ptr_cast::tptr>(generator()->eventHandler()); - Ptr::tptr sampler = - dynamic_ptr_cast::tptr>(seh->sampler()); - - CrossSection maxXSection = sampler->maxXSec(); - - weight = weight*maxXSection/nanobarn; - //std::cout << "weight " << weight << endl; - //std::cout << "maxXSection " << maxXSection/nanobarn << endl; - size_t njets = jets.size(); - size_t id = xid; - //std::cout << "number of jets from jets.size()" << njets << endl; - totalxs.count(EventContribution(0.0,weight,0.0),id); - size_t njets20 = 0; - for ( size_t k = 0; k< njets; ++k ) { - if (jets[k].perp()/GeV > 20.0 ) ++njets20; - } - //std::cout << "number of jets" << njets20 << endl; - njetsex_xs.count(EventContribution(njets20,weight,0.0),id); - - for (size_t i = 0; i < njets20; ++i) { - if (njets20 >= i) { - njetsinc_xs.count(EventContribution(i,weight,0.0),id); - } - } - - - higgsPt.count(EventContribution(h.perp()/GeV,weight,ptHiggsWidth),id); - higgsY.count(EventContribution(h.rapidity(),weight,rapWidth),id); - higgsEta.count(EventContribution(h.eta(),weight,etaWidth),id); - - // For containing HT - double HTval=0; - size_t n_max=5; // max number of jets to histogram - size_t njet_max = min(njets,n_max); - - for ( size_t k = 0; k < njets; ++k ) { - // Loop over jet momenta to calculate scalar HT - HTval+=jets[k].perp()/GeV; - } - for ( size_t k = 0; k < njet_max; ++k ) { - jetPt[k].count(EventContribution(jets[k].perp()/GeV,weight,ptWidth),id); - jetY[k].count(EventContribution(jets[k].rapidity(),weight,rapWidth),id); - jetHiggsDeltaPhi[k].count(EventContribution(deltaPhi(jets[k],h),weight,phiWidth),id); - jetHiggsDeltaY[k].count(EventContribution(deltaY(jets[k],h),weight,rapWidth),id); - jetHiggsDeltaR[k].count(EventContribution(deltaR(jets[k],h),weight,rapWidth),id); - jetHiggsMass[k].count(EventContribution(mass(jets[k],h),weight,massWidth),id); - jetMass[k].count(EventContribution(jets[k].m()/GeV,weight,jetMassWidth),id); - - - for ( size_t l = k + 1; l < njet_max; ++l ) { - jetDeltaPhi[k][l-k-1].count(EventContribution(deltaPhi(jets[k],jets[l]),weight,phiWidth),id); - jetDeltaY[k][l-k-1].count(EventContribution(deltaY(jets[k],jets[l]),weight,rapWidth),id); - jetYdotY[k][l-k-1].count(EventContribution(jets[k].rapidity()*jets[l].rapidity(),weight,rapWidth1),id); - jetDeltaR[k][l-k-1].count(EventContribution(deltaR(jets[k],jets[l]),weight,rapWidth),id); - jetPairMass[k][l-k-1].count(EventContribution(mass(jets[k],jets[l]),weight,massWidth),id); - } - } - HT.count(EventContribution(HTval,weight,HTWidth),id); - - - if ( njets > 2 ) { - double ystar = jets[2].rapidity() - 0.5 * (jets[0].rapidity()+jets[1].rapidity()); - ystar /= jets[0].rapidity()-jets[1].rapidity(); - yStarNormed.count(EventContribution(ystar,weight,ysWidth),id); - - } - - if( njets > 1 ){ - LorentzMomentum dijetSum=jets[0]+jets[1]; - double phiDiff=deltaPhi( dijetSum, h); - higgsDijetDeltaPhi.count(EventContribution( abs( phiDiff ), weight, higgsDijetDeltaPhiWidth),id); - - double higgsDijetpTSum=(dijetSum+h).perp()/GeV; - higgsDijetpT.count(EventContribution( higgsDijetpTSum, weight, higgsDijetpTWidth),id ); - } - - if ( njets > 2 ) { - threeJetMass.count( EventContribution(mass(jets[0],jets[1],jets[2]),weight,threeJetMassWidth),id); - } - - // jeppe DeltaPhi, 1001.3822 eq. 9. - if( njets > 1 ){ // need at lest one jet with higher y and one with lower y than the Higgs - - // to be changed (probably) - double smallest_rap=jets[0].rapidity(); - double largest_rap=jets[0].rapidity(); - - // Find largest and smallest delta rap w.r.t. Higgs - for ( size_t k = 0; k < njets; ++k ) { - if( smallest_rap > jets[k].rapidity() ) smallest_rap=jets[k].rapidity(); - if( largest_rap < jets[k].rapidity() ) largest_rap=jets[k].rapidity(); - } - - // If one jet has higher rapidity and one lower than Higgs - if( smallest_rap< h.rapidity() and largest_rap> h.rapidity() ){ - - // Collect jets into those having lower rapidity - // and those having higher rapidity than the Higgs - LorentzMomentum p_a; - LorentzMomentum p_b; - for ( size_t k = 0; k < njets; ++k ) { - if( jets[k].rapidity() < h.rapidity() ) p_a+=jets[k]; - else if( jets[k].rapidity() > h.rapidity() ) p_b+=jets[k]; - } - double phiDiff=deltaPhi( p_a, p_b ); - jetSumDeltaPhi.count( EventContribution(phiDiff,weight,jetSumDeltaPhiWidth),id ); - } - - } - -} - -void HJetsAnalysis2::finalizeHistogram(Histogram& h, XML::Element& elem) { - h.finalize(); - elem.append(h.toXML()); -} - -void HJetsAnalysis2::dofinish() { - AnalysisHandler::dofinish(); - - Ptr::tptr seh = - dynamic_ptr_cast::tptr>(generator()->eventHandler()); - Ptr::tptr sampler = - dynamic_ptr_cast::tptr>(seh->sampler()); - - unsigned long attemptedPoints = sampler->attempts(); - double sumOfWeights = sampler->sumWeights(); - double sumOfSquaredWeights = sampler->sumWeights2(); - CrossSection maxXSection = sampler->maxXSec(); - - XML::Element elem(XML::ElementTypes::Element,"Run"); - - elem.appendAttribute("name",generator()->runName()); - elem.appendAttribute("attemptedPoints",attemptedPoints); - elem.appendAttribute("sumOfWeights",sumOfWeights*maxXSection/nanobarn); - elem.appendAttribute("sumOfSquaredWeights",sumOfSquaredWeights*sqr(maxXSection/nanobarn)); - //elem.appendAttribute("sumOfWeights",sumOfWeights); - //elem.appendAttribute("sumOfSquaredWeights",sumOfSquaredWeights); - - XML::Element xhistos(XML::ElementTypes::Element,"Histograms"); - - finalizeHistogram(totalxs,xhistos); - finalizeHistogram(njetsex_xs,xhistos); - finalizeHistogram(njetsinc_xs,xhistos); - finalizeHistogram(higgsPt,xhistos); - finalizeHistogram(higgsY,xhistos); - finalizeHistogram(higgsEta,xhistos); - finalizeHistogram(HT,xhistos); - finalizeHistogram(higgsDijetpT,xhistos); - finalizeHistogram(threeJetMass,xhistos); - finalizeHistogram(jetSumDeltaPhi,xhistos); - - size_t n_min = 5; - - for ( size_t k = 0; k < n_min; ++k ) { - finalizeHistogram(jetPt[k],xhistos); - finalizeHistogram(jetY[k],xhistos); - finalizeHistogram(jetHiggsDeltaPhi[k],xhistos); - finalizeHistogram(jetHiggsDeltaY[k],xhistos); - finalizeHistogram(jetHiggsDeltaR[k],xhistos); - finalizeHistogram(jetHiggsMass[k],xhistos); - finalizeHistogram(jetMass[k],xhistos); - - for ( size_t l = k + 1; l < n_min; ++l ) { - finalizeHistogram(jetDeltaPhi[k][l-k-1],xhistos); - finalizeHistogram(jetDeltaY[k][l-k-1],xhistos); - finalizeHistogram(jetYdotY[k][l-k-1],xhistos); - finalizeHistogram(jetDeltaR[k][l-k-1],xhistos); - finalizeHistogram(jetPairMass[k][l-k-1],xhistos); - } - } - if ( theJetRegions.size() > 1 ) { - finalizeHistogram(higgsDijetDeltaPhi,xhistos); - } - if ( theJetRegions.size() > 2 ) { - finalizeHistogram(yStarNormed,xhistos); - } - - elem.append(xhistos); - - string fname = name() + ".xml"; - ofstream runXML(fname.c_str()); - XML::ElementIO::put(elem,runXML); - -} - -void HJetsAnalysis2::doinitrun() { - AnalysisHandler::doinitrun(); - - size_t njets = 5; //theJetRegions.size(); - - const double pi = 3.14159265; - - totalxs = Histogram("sigma_tot",Histogram::regularBinEdges(-0.5,0.5,1)); - njetsex_xs = Histogram("sigma_njetex",Histogram::regularBinEdges(-0.5,8.5,9)); - njetsinc_xs = Histogram("sigma_njetinc",Histogram::regularBinEdges(-0.5,8.5,9)); - higgsPt = Histogram("higgsPt",Histogram::regularBinEdges(0,400,40),true); - higgsY = Histogram("higgsY",Histogram::regularBinEdges(-5,5,32)); - higgsEta = Histogram("higgsEta",Histogram::regularBinEdges(-5,5,32)); - HT = Histogram("HT",Histogram::regularBinEdges(0,1000,32),true); - higgsDijetDeltaPhi = Histogram("higgsDijetDeltaPhi",Histogram::regularBinEdges(0,pi,32),true,true); - higgsDijetpT = Histogram("higgsDijetpT",Histogram::regularBinEdges(0,100,20),true); - threeJetMass = Histogram("threeJetMass",Histogram::regularBinEdges(0,2000,32),true,false); - jetSumDeltaPhi = Histogram("jetSumDeltaPhi",Histogram::regularBinEdges(-pi,pi,32),true,true); - - for ( size_t k = 0; k < njets; ++k ) { - ostringstream name; - name << "jet" << (k+1); - jetPt.push_back(Histogram(name.str() + "Pt",Histogram::regularBinEdges(0,400,40),true)); - jetY.push_back(Histogram(name.str() + "Y",Histogram::regularBinEdges(-5,5,32))); - jetHiggsDeltaPhi.push_back(Histogram(name.str() + "HiggsDeltaPhi",Histogram::regularBinEdges(-pi,pi,32),true,true)); - jetHiggsDeltaY.push_back(Histogram(name.str() + "HiggsDeltaY",Histogram::regularBinEdges(0,10,32),true)); - jetHiggsDeltaR.push_back(Histogram(name.str() + "HiggsDeltaR",Histogram::regularBinEdges(0,10,32),true)); - jetHiggsMass.push_back(Histogram(name.str() + "HiggsMass",Histogram::regularBinEdges(0,4000,64),true)); - jetMass.push_back(Histogram(name.str() + "JetMass",Histogram::regularBinEdges(0,100,32),true, false)); - - - vector xjetDeltaPhi; - vector xjetDeltaY; - vector xjetYdotY; - vector xjetDeltaR; - vector xjetPairMass; - for ( size_t l = k + 1; l < njets; ++l ) { - ostringstream xname; - xname << "jet" << (k+1) << (l+1); - - xjetDeltaPhi.push_back(Histogram(xname.str() + "DeltaPhi",Histogram::regularBinEdges(-pi,pi,32),true,true)); - xjetDeltaY.push_back(Histogram(xname.str() + "DeltaY",Histogram::regularBinEdges(0,10,32),true)); - xjetYdotY.push_back(Histogram(xname.str() + "YdotY",Histogram::regularBinEdges(-30,30,64),true)); - xjetDeltaR.push_back(Histogram(xname.str() + "DeltaR",Histogram::regularBinEdges(0,10,32),true)); - xjetPairMass.push_back(Histogram(xname.str() + "Mass",Histogram::regularBinEdges(0,4000,64),true)); - } - jetDeltaPhi.push_back(xjetDeltaPhi); - jetDeltaY.push_back(xjetDeltaY); - jetYdotY.push_back(xjetYdotY); - jetDeltaR.push_back(xjetDeltaR); - jetPairMass.push_back(xjetPairMass); - } - - yStarNormed = Histogram("yStarNormed",Histogram::regularBinEdges(-2,2,32)); - -} - - -IBPtr HJetsAnalysis2::clone() const { - return new_ptr(*this); -} - -IBPtr HJetsAnalysis2::fullclone() const { - return new_ptr(*this); -} - - -// If needed, insert default implementations of virtual function defined -// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). - - -void HJetsAnalysis2::persistentOutput(PersistentOStream & os) const { - os << theJetFinder << theJetRegions << theFullEvent - << theBinFactor << theRapidityGap << ounit(themjj_min,GeV) - << ounit(thehiggsPt_min,GeV) << thehiggsY_max << theopposite_dir << thenjets_min - << thehiggs_ingap << theFuzzy << ounit(theEnergyCutWidth,GeV) << theRapidityCutWidth; -} - -void HJetsAnalysis2::persistentInput(PersistentIStream & is, int) { - is >> theJetFinder >> theJetRegions >> theFullEvent - >> theBinFactor >> theRapidityGap >> iunit(themjj_min,GeV) - >> iunit(thehiggsPt_min,GeV) >> thehiggsY_max >> theopposite_dir >> thenjets_min - >> thehiggs_ingap >> theFuzzy >> iunit(theEnergyCutWidth,GeV) >> theRapidityCutWidth; -} - - -// *** 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 - describeHJetsHJetsAnalysis2("HJets::HJetsAnalysis2", "HJetsAnalysis2.so"); - -void HJetsAnalysis2::Init() { - - static ClassDocumentation documentation - ("There is no documentation for the HJetsAnalysis2 class"); - - - static Reference interfaceJetFinder - ("JetFinder", - "", - &HJetsAnalysis2::theJetFinder, false, false, true, false, false); - - static RefVector interfaceJetRegions - ("JetRegions", - "", - &HJetsAnalysis2::theJetRegions, -1, false, false, true, false, false); - - - static Switch interfaceFullEvent - ("FullEvent", - "Analyse a fully simulated event", - &HJetsAnalysis2::theFullEvent, false, false, false); - static SwitchOption interfaceFullEventYes - (interfaceFullEvent, - "Yes", - "", - true); - static SwitchOption interfaceFullEventNo - (interfaceFullEvent, - "No", - "", - false); - - static Parameter interfaceBinFactor - ("BinFactor", - "The factor multiplying the bin width", - &HJetsAnalysis2::theBinFactor, 1.0, 0.0, 1.0, true, false, true ); - - static Parameter interfaceRapidityGap - ("RapidityGap", - "The distance in pseudorapidity between the two tagging jets", - &HJetsAnalysis2::theRapidityGap, 0.0, 0.0, 20.0, true, false, true); - - - static Parameter interfacemjj_min - ("mjj_min", - "The minimum dijet invariant mass of the two taggin jets.", - &HJetsAnalysis2::themjj_min, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, - false, false, Interface::lowerlim); - - - static Parameter interfacehiggsPt_min - ("higgsPt_min", - "Minimun transverse momentum of higgs.", - &HJetsAnalysis2::thehiggsPt_min, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, - false, false, Interface::lowerlim); - - static Parameter interfacehiggsY_max - ("higgsY_max", - "maximum absolute rapidity of higgs", - &HJetsAnalysis2::thehiggsY_max, 0, 0, 14000, true, false, true); - - - - static Switch interfaceopposite_dir - ("opposite_dir", - "Should the tagging jets reside in opposite directions (Yes/No)?", - &HJetsAnalysis2::theopposite_dir, false, false, false); - static SwitchOption interfaceopposite_dirYes - (interfaceopposite_dir, - "Yes", - "", - true); - static SwitchOption interfaceopposite_dirNo - (interfaceopposite_dir, - "No", - "", - false); - - static Parameter interfacenjets_min - ("njets_min", - "The minimum number of jets", - &HJetsAnalysis2::thenjets_min, 2, 1, 10, true, false, true); - - - static Switch interfacehiggs_ingap - ("higgs_ingap", - "Should the higgs reside in the gap (Yes/No)?", - &HJetsAnalysis2::thehiggs_ingap, false, false, false); - static SwitchOption interfacehiggs_ingapYes - (interfacehiggs_ingap, - "Yes", - "", - true); - static SwitchOption interfacehiggs_ingapNo - (interfacehiggs_ingap, - "No", - "", - false); - - static Switch interfaceFuzzy - ("Fuzzy", - "Make this jet region a fuzzy cut", - &HJetsAnalysis2::theFuzzy, false, false, false); - static SwitchOption interfaceFuzzyOn - (interfaceFuzzy, - "Yes", - "", - true); - static SwitchOption interfaceFuzzyOff - (interfaceFuzzy, - "No", - "", - false); - - static Parameter interfaceEnergyCutWidth - ("EnergyCutWidth", - "The pt cut smearing.", - &HJetsAnalysis2::theEnergyCutWidth, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, - false, false, Interface::lowerlim); - - static Parameter interfaceRapidityCutWidth - ("RapidityCutWidth", - "The rapidity cut smearing.", - &HJetsAnalysis2::theRapidityCutWidth, 0.1, 0.0, 0, - false, false, Interface::lowerlim); - -} - diff --git a/Analysis/HJetsAnalysis2.h b/Analysis/HJetsAnalysis2.h deleted file mode 100644 --- a/Analysis/HJetsAnalysis2.h +++ /dev/null @@ -1,394 +0,0 @@ -// -*- C++ -*- -#ifndef HJets_HJetsAnalysis2_H -#define HJets_HJetsAnalysis2_H -// -// This is the declaration of the HJetsAnalysis2 class. -// - -#include "ThePEG/Handlers/AnalysisHandler.h" -#include "ThePEG/Cuts/JetFinder.h" -#include "ThePEG/Cuts/JetRegion.h" -#include "ThePEG/Interface/Parameter.h" -#include "ThePEG/EventRecord/SubProcess.h" -#include "ThePEG/EventRecord/SubProcessGroup.h" - -#include "Herwig++/Utilities/Statistics/Histogram.h" - -namespace HJets { - -using namespace ThePEG; - -/** - * Here is the documentation of the HJetsAnalysis2 class. - * - * @see \ref HJetsAnalysis2Interfaces "The interfaces" - * defined for HJetsAnalysis2. - */ -class HJetsAnalysis2: public AnalysisHandler { - -public: - - /** @name Standard constructors and destructors. */ - //@{ - /** - * The default constructor. - */ - HJetsAnalysis2(); - - /** - * The destructor. - */ - virtual ~HJetsAnalysis2(); - //@} - -public: - - /** @name Virtual functions required by the AnalysisHandler class. */ - //@{ - /** - * Analyze a given Event. Note that a fully generated event - * may be presented several times, if it has been manipulated in - * between. The default version of this function will call transform - * to make a lorentz transformation of the whole event, then extract - * all final state particles and call analyze(tPVector) of this - * analysis object and those of all associated analysis objects. The - * default version will not, however, do anything on events which - * have not been fully generated, or have been manipulated in any - * way. - * @param event pointer to the Event to be analyzed. - * @param ieve the event number. - * @param loop the number of times this event has been presented. - * If negative the event is now fully generated. - * @param state a number different from zero if the event has been - * manipulated in some way since it was last presented. - */ - virtual void analyze(tEventPtr event, long ieve, int loop, int state); - //@} - - /** - * Analyze one subprocess, given the event number it belongs to - */ - void analyze(const ParticleVector&, long, double); - - /** - * Analyze higgs and jet momenta, given the event number it belongs to - */ - void analyze(const LorentzMomentum& h, - const vector& jets, - long, double); - - -protected: - - /** - * Initialize this object. Called in the run phase just before - * a run begins. - */ - virtual void doinitrun(); - - /** - * Finalize this object. Called in the run phase just after a - * run has ended. Used eg. to write out statistics. - */ - virtual void dofinish(); - - - -public: - - /** @name Functions used by the persistent I/O system. */ - //@{ - /** - * Function used to write out object persistently. - * @param os the persistent output stream written to. - */ - void persistentOutput(PersistentOStream & os) const; - - /** - * Function used to read in object persistently. - * @param is the persistent input stream read from. - * @param version the version number of the object when written. - */ - void persistentInput(PersistentIStream & is, int version); - //@} - - /** - * The standard Init function used to initialize the interfaces. - * Called exactly once for each class by the class description system - * before the main function starts or - * when this class is dynamically loaded. - */ - static void Init(); - - - - -protected: - - /** @name Clone Methods. */ - //@{ - /** - * Make a simple clone of this object. - * @return a pointer to the new object. - */ - virtual IBPtr clone() const; - - /** Make a clone of this object, possibly modifying the cloned object - * to make it sane. - * @return a pointer to the new object. - */ - virtual IBPtr fullclone() const; - //@} - - -// If needed, insert declarations of virtual function defined in the -// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). - -protected: - - /** - * Return true, if this jet region is fuzzy - */ - bool fuzzy() const { return theFuzzy; } - - - -private: - - - /** - * Check for cuts - */ - bool selectionCuts(const LorentzMomentum& h, const vector& jets, double& weight); - - - /** - * Perform a (potentially) fuzzy check on energy-type quantities - */ - bool lessThanEnergy(Energy a, Energy b, double& weight) const; - - /** - * Perform a (potentially) fuzzy check on angular-type quantities - */ - bool lessThanRapidity(double a, double b, double& weight) const; - - - /** - * The jet finder to use - */ - Ptr::ptr theJetFinder; - - /** - * The jet regions to match. - */ - vector::ptr> theJetRegions; - - - /** - * True if full evnet is taken - */ - bool theFullEvent; - - - /** - * Factor multiplying bin width (bin smearing) - */ - double theBinFactor; - - /** - * Rapidity gap - */ - double theRapidityGap; - - /** - * di-jet mass cut - */ - - Energy themjj_min; - - /** - * Higgs Pt min - */ - Energy thehiggsPt_min; - - /** - * Higgs max rapidity - */ - - double thehiggsY_max; - - /** - * True if leading jets are opposite in rapidity - */ - bool theopposite_dir; - /** - * Mininum number of analyzed jets - */ - unsigned int thenjets_min; - - /** - * True if higgs is rapidity gap - */ - bool thehiggs_ingap; - - /** - * True if Fuzzy cuts - */ - bool theFuzzy; - - - /** - * The smearing width for the pt or mass cuts, if fuzzy - */ - Energy theEnergyCutWidth; - - /** - * The smearing width for the rapidity cut, if fuzzy - */ - double theRapidityCutWidth; - - private: - - /** - * The assignment operator is private and must never be called. - * In fact, it should not even be implemented. - */ - HJetsAnalysis2 & operator=(const HJetsAnalysis2 &); - - /** - * Histogram indices for total xs - */ - Statistics::Histogram totalxs; - - /** - * Histogram indices for exclusive xs - */ - Statistics::Histogram njetsex_xs; - - /** - * Histogram indices for inclusive njet xs - */ - Statistics::Histogram njetsinc_xs; - - - /** - * Histogram indices for higgs pts - */ - Statistics::Histogram higgsPt; - - /** - * Statistics::Histogram indices for higgs rapidities - */ - Statistics::Histogram higgsY; - - /** - * Statistics::Histogram indices for higgs pseudo rapidities - */ - Statistics::Histogram higgsEta; - - /** - * Statistics::Scalar sum of transverse momenta - */ - Statistics::Histogram HT; - - /** - * Statistics::Histogram of the transverse momentum - * of higgs plus the 2 hardest jets - */ - Statistics::Histogram higgsDijetpT; - - /** - * Statistics::Histogram of delta phi between higgs - * and the two hardest jets - */ - Statistics::Histogram higgsDijetDeltaPhi; - - - /** - * Statistics::Histogram of the invariant mass of - * the three hardest jets - */ - Statistics::Histogram threeJetMass; - - /** - * Statistics::Histogram of the azimuthal angle as defined - * in 1001.3822 eq. 9. - */ - Statistics::Histogram jetSumDeltaPhi; - - - /** - * Statistics::Histogram indices for jet pts - */ - vector jetPt; - - /** - * Statistics::Histogram indices for jet rapidities - */ - vector jetY; - - /** - * Statistics::Histogram indices for jet delta phis - */ - vector > jetDeltaPhi; - - /** - * Statistics::Histogram indices for jet y dot y - */ - vector > jetYdotY; - - /** - * Statistics::Histogram indices for jet delta y - */ - vector > jetDeltaY; - - /** - * Statistics::Histogram indices for jet delta R - */ - vector > jetDeltaR; - - /** - * Statistics::Histogram indices for jet pair masses - */ - vector > jetPairMass; - - /** - * Statistics::Histogram indices for jet higgs delta phis - */ - vector jetHiggsDeltaPhi; - - /** - * Statistics::Histogram indices for jet higgs delta y - */ - vector jetHiggsDeltaY; - - /** - * Statistics::Histogram indices for jet higgs delta R - */ - vector jetHiggsDeltaR; - - /** - * Statistics::Histogram indices for jet higgs masses - */ - vector jetHiggsMass; - - /** - * Statistics::Histogram for invariant mass of jets - */ - vector jetMass; - - /** - * y* - */ - Statistics::Histogram yStarNormed; - - /** - * Finalize histogram and add to XML element - */ - void finalizeHistogram(Statistics::Histogram&,XML::Element&); - -}; - -} - -#endif /* HJets_HJetsAnalysis2_H */ diff --git a/Analysis/Makefile.am b/Analysis/Makefile.am deleted file mode 100644 --- a/Analysis/Makefile.am +++ /dev/null @@ -1,6 +0,0 @@ -pkglib_LTLIBRARIES = HJetsAnalysis2.la -HJetsAnalysis2_la_LDFLAGS = -module -version-info 1:1:0 - -HJetsAnalysis2_la_SOURCES = \ -HJetsAnalysis2.cc \ -HJetsAnalysis2.h diff --git a/COPYING b/COPYING --- a/COPYING +++ b/COPYING @@ -0,0 +1,340 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License. diff --git a/Makefile.am b/Makefile.am --- a/Makefile.am +++ b/Makefile.am @@ -1,8 +1,7 @@ SUBDIRS = \ include \ Virtuals \ Amplitudes \ -Analysis \ src ACLOCAL_AMFLAGS = -I m4 diff --git a/README b/README --- a/README +++ b/README @@ -0,0 +1,43 @@ +################################################################################ +# # +# 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++ + +--with-herwig=path + +informs HJets++'s configure where to look for the Herwig++ +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. + diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -1,49 +1,48 @@ 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_CONFIG_AUX_DIR([config]) AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR(["m4/$PACKAGE_TARNAME.m4"]) 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_FC([gfortran]) AC_PROG_F77([gfortran]) AC_PROG_CXX AC_PROG_INSTALL AC_PROG_MAKE_SET AC_PROG_LN_S AM_PROG_AR AC_DISABLE_STATIC AC_LIBTOOL_DLOPEN AC_PROG_LIBTOOL -HJETS_CHECK_MATCHBOX +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 - Analysis/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++ -- -AC_DEFUN([HJETS_CHECK_MATCHBOX], +AC_DEFUN([HJETS_CHECK_HERWIG], [ defaultlocation="${prefix}" test "x$defaultlocation" = xNONE && defaultlocation="${ac_default_prefix}" -AC_MSG_CHECKING([for Matchbox/Matchbox in]) -AC_ARG_WITH(matchbox, - AC_HELP_STRING([--with-matchbox=DIR],[location of Matchbox installation]), +AC_MSG_CHECKING([for Herwig++/Matchbox in]) +AC_ARG_WITH(herwig, + AC_HELP_STRING([--with-herwig=DIR],[location of Herwig++ installation]), [], - [with_matchbox="${defaultlocation}"]) -AC_MSG_RESULT([$with_matchbox]) + [with_herwig="${defaultlocation}"]) +AC_MSG_RESULT([$with_herwig]) -AS_IF([test "x$with_matchbox" != "xno"], +AS_IF([test "x$with_herwig" != "xno"], [AC_CHECK_FILES( - ${with_matchbox}/bin/herwig-config, - [have_matchbox=yes], [have_matchbox=no])], - [have_matchbox=no]) + ${with_herwig}/bin/herwig-config, + [have_herwig=yes], [have_herwig=no])], + [have_herwig=no]) -AS_IF([test "x$have_matchbox" = "xyes"], - [MATCHBOXCPPFLAGS=`${with_matchbox}/bin/herwig-config --cppflags` - MATCHBOXLDFLAGS=`${with_matchbox}/bin/herwig-config --ldflags` - MATCHBOXLDLIBS=`${with_matchbox}/bin/herwig-config --ldlibs` +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_matchbox" != "xno"], - [AC_MSG_ERROR([Cannot build HJets without Matchbox. Please set --with-matchbox.]) + [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 $MATCHBOXCPPFLAGS" -AM_LDFLAGS="$MATCHBOXLDFLAGS $MATCHBOXLDLIBS" +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) ])