Page MenuHomeHEPForge
Files F8725604

No OneTemporary
Actions

View Options

diff --git a/Analysis/LeptonsJetsAnalysis.cc b/Analysis/LeptonsJetsAnalysis.cc
new file mode 100644
--- /dev/null
+++ b/Analysis/LeptonsJetsAnalysis.cc
@@ -0,0 +1,594 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the LeptonsJetsAnalysis class.
+//
+
+#include "LeptonsJetsAnalysis.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+
+#include "ThePEG/Interface/Reference.h"
+#include "ThePEG/Interface/RefVector.h"
+#include "ThePEG/Interface/Switch.h"
+
+#include "ThePEG/EventRecord/SubProcess.h"
+#include "ThePEG/EventRecord/SubProcessGroup.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 Herwig;
+
+LeptonsJetsAnalysis::LeptonsJetsAnalysis()
+ : theIsShowered(false), theApplyCuts(false) {}
+
+LeptonsJetsAnalysis::~LeptonsJetsAnalysis() {}
+
+
+
+#ifndef LWH_AIAnalysisFactory_H
+#ifndef LWH
+#define LWH ThePEGLWH
+#endif
+#include "ThePEG/Analysis/LWH/AnalysisFactory.h"
+#endif
+
+struct SortPt {
+
+ inline bool operator()(const LorentzMomentum& a,
+ const LorentzMomentum& b) const {
+ return a.perp() > b.perp();
+ }
+
+};
+
+struct SortId {
+
+ inline bool operator()(const pair<PID,LorentzMomentum>& a,
+ const pair<PID,LorentzMomentum>& b) const {
+ // sort by abs(pid); if equal, first particle, then antiparticle
+ // this puts pairs forming bosons next to each other
+ long p1 = a.first;
+ long p2 = b.first;
+ if (abs(p1)==abs(p2)) {
+ return p1 > p2;
+ } else {
+ return abs(p1) < abs(p2);
+ }
+ }
+
+};
+
+void LeptonsJetsAnalysis::reconstructJets(const ParticleVector& parts) {
+
+ tcPDVector outType;
+ vector<LorentzMomentum> outMomenta;
+
+ for ( ParticleVector::const_iterator p = parts.begin();
+ p != parts.end(); ++p ) {
+ outType.push_back((**p).dataPtr());
+ outMomenta.push_back((**p).momentum());
+ }
+
+ jetFinder()->cluster(outType, outMomenta,
+ tcCutsPtr(), tcPDPtr(),tcPDPtr());
+
+ sort(outMomenta.begin(),outMomenta.end(),SortPt());
+
+ for ( vector<Ptr<JetRegion>::ptr>::iterator j =
+ theJetRegions.begin(); j != theJetRegions.end(); ++j )
+ (**j).reset();
+
+ for ( size_t k = 0; k < outMomenta.size(); ++k ) {
+ for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = jetRegions().begin();
+ r != jetRegions().end(); ++r ) {
+ if ( (**r).matches(tCutsPtr(),k+1,outMomenta[k])) {
+ jetMomentum(k+1) = outMomenta[k];
+ break;
+ }
+ }
+ }
+
+}
+
+void LeptonsJetsAnalysis::reconstructEWParticles(ParticleVector& parts) {
+
+ vector< pair<PID,LorentzMomentum> > partall;
+ vector<LorentzMomentum> partl, partnu, parth;
+ LorentzMomentum ptmiss = LorentzMomentum(ZERO,ZERO,ZERO,ZERO);;
+
+ ParticleVector::iterator p = parts.begin();
+ while (p != parts.end()) {
+ PID pid = (**p).id();
+ if ( ( static_cast<long>(pid) == ParticleID::eminus ) ||
+ ( static_cast<long>(pid) == ParticleID::eplus ) ||
+ ( static_cast<long>(pid) == ParticleID::muminus ) ||
+ ( static_cast<long>(pid) == ParticleID::muplus ) ||
+ ( static_cast<long>(pid) == ParticleID::tauminus ) ||
+ ( static_cast<long>(pid) == ParticleID::tauplus ) ) {
+ partall.push_back(pair<PID,LorentzMomentum>(pid,(**p).momentum()));
+ partl.push_back((**p).momentum());
+ p = parts.erase(p);
+ } else
+ if ( ( static_cast<long>(pid) == ParticleID::nu_e ) ||
+ ( static_cast<long>(pid) == ParticleID::nu_ebar ) ||
+ ( static_cast<long>(pid) == ParticleID::nu_mu ) ||
+ ( static_cast<long>(pid) == ParticleID::nu_mubar ) ||
+ ( static_cast<long>(pid) == ParticleID::nu_tau ) ||
+ ( static_cast<long>(pid) == ParticleID::nu_taubar ) ) {
+ partall.push_back(pair<PID,LorentzMomentum>(pid,(**p).momentum()));
+ partnu.push_back((**p).momentum());
+ ptmiss += (**p).momentum();
+ p = parts.erase(p);
+ } else
+ if ( static_cast<long>(pid) == ParticleID::h0 ) {
+ partall.push_back(pair<PID,LorentzMomentum>(pid,(**p).momentum()));
+ parth.push_back((**p).momentum());
+ p = parts.erase(p);
+ } else
+ p++;
+
+ }
+
+ sort(partall.begin(),partall.end(),SortId());
+ sort(partl.begin(),partl.end(),SortPt());
+ sort(partnu.begin(),partnu.end(),SortPt());
+ sort(parth.begin(),parth.end(),SortPt());
+
+ // make missing transverse momentum transverse and add as last entry in leptonID and neutrino
+ ptmiss.setE(ptmiss.perp());
+ ptmiss.setZ(0*GeV);
+ partall.push_back(pair<PID,LorentzMomentum>(ParticleID::nu_e,ptmiss));
+ partnu.push_back(ptmiss);
+
+ for ( size_t k = 0; k < partall.size(); ++k )
+ leptonIDMomentum(k+1) = partall[k].second;
+ for ( size_t k = 0; k < partl.size(); ++k )
+ leptonPTMomentum(k+1) = partl[k];
+ for ( size_t k = 0; k < partnu.size(); ++k )
+ neutrinoMomentum(k+1) = partnu[k];
+ for ( size_t k = 0; k < parth.size(); ++k )
+ higgsMomentum(k+1) = parth[k];
+
+}
+
+void LeptonsJetsAnalysis::analyze(ParticleVector& parts, long id, double weight) {
+
+ clear();
+ reconstructEWParticles(parts);
+ reconstructJets(parts);
+
+ if ( theApplyCuts ) {
+// VBF cuts
+ if ( nJets()<2 ) return;
+ if ( (jetMomentum(1)+jetMomentum(2)).m() < 600*GeV ) return;
+ if ( abs(jetMomentum(1).rapidity()-jetMomentum(2).rapidity()) < 3.6 ) return;
+ if ( jetMomentum(1).rapidity()*jetMomentum(2).rapidity() > 0 ) return;
+ for ( map<unsigned int,LorentzMomentum>::const_iterator h = theLeptonPTs.begin();
+ h != theLeptonPTs.end(); ++h ) {
+ if ( h->second.perp() < 20*GeV ) return;
+ if ( abs(h->second.rapidity()) > 2.5 ) return;
+ }
+ }
+
+ unsigned int njets = 0;
+ Energy jetSummedPerp = ZERO;
+ double jetSummedRapidity = 0.0;
+ double jetSummedPhi = 0.0;
+ Energy jetSummedM = ZERO;
+
+ nJetsInclusive().count(Statistics::EventContribution(njets,weight,0.0),id);
+ if ( njets == theJets.size() )
+ nJetsExclusive().count(Statistics::EventContribution(njets,weight,0.0),id);
+
+ for ( map<unsigned int,LorentzMomentum>::const_iterator h = theJets.begin();
+ h != theJets.end(); ++h ) {
+ njets += 1;
+ jetProperties(h->first).count(h->second,weight,id);
+ jetInclusiveProperties().count(h->second,weight,id);
+ nJetsInclusive().count(Statistics::EventContribution(njets,weight,0.0),id);
+ if ( njets == theJets.size() ) {
+ exclusiveJetProperties(h->first).count(h->second,weight,id);
+ nJetsExclusive().count(Statistics::EventContribution(njets,weight,0.0),id);
+ }
+ jetSummedPerp += h->second.perp();
+ jetSummedRapidity += h->second.rapidity();
+ jetSummedPhi += h->second.phi();
+ jetSummedM += h->second.m();
+ map<unsigned int,LorentzMomentum>::const_iterator g = h; ++g;
+ for ( ; g != theJets.end(); ++g ) {
+ jetPairProperties(h->first,g->first).count(h->second,g->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g1 = g; ++g1;
+ for ( ; g1 != theJets.end(); ++g1 ) {
+ threeJetProperties(h->first,g->first,g1->first).count(h->second,g->second,g1->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g2 = g1; ++g2;
+ for ( ; g2 != theJets.end(); ++g2 ) {
+ LorentzMomentum p1234 =
+ h->second + g->second + g1->second + g2->second;
+ fourJetProperties(h->first,g->first,g1->first,g2->first).count(p1234,weight,id);
+ }
+ }
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g1 = theLeptonIDs.begin();
+ g1 != theLeptonIDs.end(); ++g1 )
+ jetPairLeptonIDTripleProperties(h->first,g->first,g1->first).count(h->second,g->second,g1->second,weight,id);
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g1 = theLeptonPTs.begin();
+ g1 != theLeptonPTs.end(); ++g1 )
+ jetPairLeptonPTTripleProperties(h->first,g->first,g1->first).count(h->second,g->second,g1->second,weight,id);
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g1 = theNeutrinos.begin();
+ g1 != theNeutrinos.end(); ++g1 )
+ jetPairNeutrinoTripleProperties(h->first,g->first,g1->first).count(h->second,g->second,g1->second,weight,id);
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g1 = theHiggs.begin();
+ g1 != theHiggs.end(); ++g1 )
+ jetPairHiggsTripleProperties(h->first,g->first,g1->first).count(h->second,g->second,g1->second,weight,id);
+ }
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g = theLeptonIDs.begin();
+ g != theLeptonIDs.end(); ++g )
+ jetLeptonIDPairProperties(h->first,g->first).count(h->second,g->second,weight,id);
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g = theLeptonPTs.begin();
+ g != theLeptonPTs.end(); ++g )
+ jetLeptonPTPairProperties(h->first,g->first).count(h->second,g->second,weight,id);
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g = theNeutrinos.begin();
+ g != theNeutrinos.end(); ++g )
+ jetNeutrinoPairProperties(h->first,g->first).count(h->second,g->second,weight,id);
+ for ( map<unsigned int,LorentzMomentum>::const_iterator g = theHiggs.begin();
+ g != theHiggs.end(); ++g )
+ jetHiggsPairProperties(h->first,g->first).count(h->second,g->second,weight,id);
+ }
+
+ if ( njets > 0 )
+ jetSummedProperties().count(jetSummedPerp,jetSummedRapidity,
+ jetSummedPhi,jetSummedM,
+ weight,id);
+
+ if ( njets > 0 )
+ jetAverageProperties().count(jetSummedPerp/njets,jetSummedRapidity/njets,
+ jetSummedPhi/njets,jetSummedM/njets,
+ weight,id);
+
+ for ( map<unsigned int,LorentzMomentum>::const_iterator h = theLeptonIDs.begin();
+ h != theLeptonIDs.end(); ++h ) {
+ leptonIDProperties(h->first).count(h->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g = h; ++g;
+ for ( ; g != theLeptonIDs.end(); ++g ) {
+ leptonIDPairProperties(h->first,g->first).count(h->second,g->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g1 = g; ++g1;
+ for ( ; g1 != theLeptonIDs.end(); ++g1 ) {
+ threeLeptonIDProperties(h->first,g->first,g1->first).count(h->second,g->second,g1->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g2 = g1; ++g2;
+ for ( ; g2 != theLeptonIDs.end(); ++g2 ) {
+ LorentzMomentum p1234 =
+ h->second + g->second + g1->second + g2->second;
+ fourLeptonIDProperties(h->first,g->first,g1->first,g2->first).count(p1234,weight,id);
+ }
+ }
+ }
+ }
+
+ for ( map<unsigned int,LorentzMomentum>::const_iterator h = theLeptonPTs.begin();
+ h != theLeptonPTs.end(); ++h ) {
+ leptonPTProperties(h->first).count(h->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g = h; ++g;
+ for ( ; g != theLeptonPTs.end(); ++g ) {
+ leptonPTPairProperties(h->first,g->first).count(h->second,g->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g1 = g; ++g1;
+ for ( ; g1 != theLeptonPTs.end(); ++g1 ) {
+ threeLeptonPTProperties(h->first,g->first,g1->first).count(h->second,g->second,g1->second,weight,id);
+ map<unsigned int,LorentzMomentum>::const_iterator g2 = g1; ++g2;
+ for ( ; g2 != theLeptonPTs.end(); ++g2 ) {
+ LorentzMomentum p1234 =
+ h->second + g->second + g1->second + g2->second;
+ fourLeptonPTProperties(h->first,g->first,g1->first,g2->first).count(p1234,weight,id);
+ }
+ }
+ }
+ }
+
+ for ( map<unsigned int,LorentzMomentum>::const_iterator h = theNeutrinos.begin();
+ h != theNeutrinos.end(); ++h ) {
+ neutrinoProperties(h->first).count(h->second,weight,id);
+ }
+
+ for ( map<unsigned int,LorentzMomentum>::const_iterator h = theHiggs.begin();
+ h != theHiggs.end(); ++h ) {
+ higgsProperties(h->first).count(h->second,weight,id);
+ }
+
+ analyzeSpecial(id,weight);
+
+}
+
+void LeptonsJetsAnalysis::analyze(tEventPtr event, long ieve, int loop, int state) {
+ AnalysisHandler::analyze(event, ieve, loop, state);
+
+ Ptr<StandardEventHandler>::tptr seh =
+ dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
+ Ptr<GeneralSampler>::tptr sampler =
+ dynamic_ptr_cast<Ptr<GeneralSampler>::tptr>(seh->sampler());
+
+ double norm = sampler->maxXSec()/picobarn;
+
+ if ( !theIsShowered ) {
+
+ tSubProPtr sub = event->primarySubProcess();
+ Ptr<SubProcessGroup>::tptr grp =
+ dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(sub);
+
+ ParticleVector hfs = sub->outgoing();
+
+ analyze(hfs,ieve,norm*event->weight()*sub->groupWeight());
+
+ if ( grp ) {
+
+ for ( SubProcessVector::const_iterator s = grp->dependent().begin();
+ s != grp->dependent().end(); ++s ) {
+ ParticleVector fs = (**s).outgoing();
+ analyze(fs,ieve,norm*event->weight()*(**s).groupWeight());
+ }
+
+ }
+
+ } else {
+
+ ParticleVector fs;
+ event->getFinalState(fs);
+ analyze(fs,ieve,norm*event->weight());
+
+ }
+
+}
+
+void LeptonsJetsAnalysis::dofinish() {
+ AnalysisHandler::dofinish();
+
+ Ptr<StandardEventHandler>::tptr seh =
+ dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
+ Ptr<GeneralSampler>::tptr sampler =
+ dynamic_ptr_cast<Ptr<GeneralSampler>::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));
+
+ XML::Element xhistos(XML::ElementTypes::Element,"Histograms");
+
+ for ( map<unsigned int,ObjectProperties>::iterator h = theJetProperties.begin();
+ h != theJetProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<unsigned int,ObjectProperties>::iterator h = theExclusiveJetProperties.begin();
+ h != theExclusiveJetProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ if ( !theJetInclusiveProperties.pt.bins().empty() ) {
+ theJetInclusiveProperties.finalize(xhistos);
+ }
+
+ if ( !theJetSummedProperties.pt.bins().empty() ) {
+ theJetSummedProperties.finalize(xhistos);
+ }
+
+ if ( !theJetAverageProperties.pt.bins().empty() ) {
+ theJetAverageProperties.finalize(xhistos);
+ }
+
+ if ( !theNJetsInclusive.bins().empty() ) {
+ theNJetsInclusive.finalize();
+ xhistos.append(theNJetsInclusive.toXML());
+ }
+
+ if ( !theNJetsExclusive.bins().empty() ) {
+ theNJetsExclusive.finalize();
+ xhistos.append(theNJetsExclusive.toXML());
+ }
+
+ for ( map<unsigned int,ObjectProperties>::iterator h = theLeptonIDProperties.begin();
+ h != theLeptonIDProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<unsigned int,ObjectProperties>::iterator h = theLeptonPTProperties.begin();
+ h != theLeptonPTProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<unsigned int,ObjectProperties>::iterator h = theNeutrinoProperties.begin();
+ h != theNeutrinoProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<unsigned int,ObjectProperties>::iterator h = theHiggsProperties.begin();
+ h != theHiggsProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = theJetPairProperties.begin();
+ h != theJetPairProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = theJetLeptonIDPairProperties.begin();
+ h != theJetLeptonIDPairProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = theJetLeptonPTPairProperties.begin();
+ h != theJetLeptonPTPairProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = theJetNeutrinoPairProperties.begin();
+ h != theJetNeutrinoPairProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = theJetHiggsPairProperties.begin();
+ h != theJetHiggsPairProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = theLeptonIDPairProperties.begin();
+ h != theLeptonIDPairProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = theLeptonPTPairProperties.begin();
+ h != theLeptonPTPairProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator h =
+ theThreeJetProperties.begin(); h != theThreeJetProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator h =
+ theJetPairLeptonIDTripleProperties.begin(); h != theJetPairLeptonIDTripleProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator h =
+ theJetPairLeptonPTTripleProperties.begin(); h != theJetPairLeptonPTTripleProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator h =
+ theJetPairNeutrinoTripleProperties.begin(); h != theJetPairNeutrinoTripleProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator h =
+ theJetPairHiggsTripleProperties.begin(); h != theJetPairHiggsTripleProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator h =
+ theThreeLeptonIDProperties.begin(); h != theThreeLeptonIDProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator h =
+ theThreeLeptonPTProperties.begin(); h != theThreeLeptonPTProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator h =
+ theFourJetProperties.begin(); h != theFourJetProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator h =
+ theFourLeptonIDProperties.begin(); h != theFourLeptonIDProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ for ( map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator h =
+ theFourLeptonPTProperties.begin(); h != theFourLeptonPTProperties.end(); ++h ) {
+ h->second.finalize(xhistos);
+ }
+
+ finalize(xhistos);
+
+ elem.append(xhistos);
+
+ string fname = name() + ".xml";
+ ofstream runXML(fname.c_str());
+ XML::ElementIO::put(elem,runXML);
+
+}
+
+IBPtr LeptonsJetsAnalysis::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr LeptonsJetsAnalysis::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 LeptonsJetsAnalysis::persistentOutput(PersistentOStream & os) const {
+ os << theIsShowered << theApplyCuts << theJetFinder << theJetRegions;
+}
+
+void LeptonsJetsAnalysis::persistentInput(PersistentIStream & is, int) {
+ is >> theIsShowered >> theApplyCuts >> theJetFinder >> theJetRegions;
+}
+
+
+// *** 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<LeptonsJetsAnalysis,AnalysisHandler>
+ describeHerwigLeptonsJetsAnalysis("Herwig::LeptonsJetsAnalysis", "JetCuts.so HwJetsAnalysis.so");
+
+void LeptonsJetsAnalysis::Init() {
+
+ static ClassDocumentation<LeptonsJetsAnalysis> documentation
+ ("There is no documentation for the LeptonsJetsAnalysis class");
+
+ static Reference<LeptonsJetsAnalysis,JetFinder> interfaceJetFinder
+ ("JetFinder",
+ "",
+ &LeptonsJetsAnalysis::theJetFinder, false, false, true, false, false);
+
+ static RefVector<LeptonsJetsAnalysis,JetRegion> interfaceJetRegions
+ ("JetRegions",
+ "",
+ &LeptonsJetsAnalysis::theJetRegions, -1, false, false, true, false, false);
+
+ static Switch<LeptonsJetsAnalysis,bool> interfaceIsShowered
+ ("IsShowered",
+ "",
+ &LeptonsJetsAnalysis::theIsShowered, false, false, false);
+ static SwitchOption interfaceIsShoweredYes
+ (interfaceIsShowered,
+ "Yes",
+ "",
+ true);
+ static SwitchOption interfaceIsShoweredNo
+ (interfaceIsShowered,
+ "No",
+ "",
+ false);
+
+ static Switch<LeptonsJetsAnalysis,bool> interfaceApplyCuts
+ ("ApplyCuts",
+ "",
+ &LeptonsJetsAnalysis::theApplyCuts, false, false, false);
+ static SwitchOption interfaceApplyCutsYes
+ (interfaceApplyCuts,
+ "Yes",
+ "",
+ true);
+ static SwitchOption interfaceApplyCutsNo
+ (interfaceApplyCuts,
+ "No",
+ "",
+ false);
+
+}
+
diff --git a/Analysis/LeptonsJetsAnalysis.h b/Analysis/LeptonsJetsAnalysis.h
new file mode 100644
--- /dev/null
+++ b/Analysis/LeptonsJetsAnalysis.h
@@ -0,0 +1,1058 @@
+// -*- C++ -*-
+#ifndef Herwig_LeptonsJetsAnalysis_H
+#define Herwig_LeptonsJetsAnalysis_H
+//
+// This is the declaration of the LeptonsJetsAnalysis class.
+//
+
+#include "ThePEG/Handlers/AnalysisHandler.h"
+#include "ThePEG/Cuts/JetFinder.h"
+#include "ThePEG/Cuts/JetRegion.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "Herwig++/Utilities/Statistics/Histogram.h"
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/tuple/tuple_comparison.hpp>
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the LeptonsJetsAnalysis class.
+ *
+ * @see \ref LeptonsJetsAnalysisInterfaces "The interfaces"
+ * defined for LeptonsJetsAnalysis.
+ */
+class LeptonsJetsAnalysis: public AnalysisHandler {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ LeptonsJetsAnalysis();
+
+ /**
+ * The destructor.
+ */
+ virtual ~LeptonsJetsAnalysis();
+ //@}
+
+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);
+ //@}
+
+protected:
+
+ /**
+ * Analyze one subprocess, given the event number it belongs to
+ */
+ void analyze(ParticleVector&, long, double);
+
+ /**
+ * Clear the electroweak objects and jets for the next event
+ */
+ void clear() {
+ theJets.clear();
+ theLeptonIDs.clear();
+ theLeptonPTs.clear();
+ theNeutrinos.clear();
+ theHiggs.clear();
+ }
+
+ /**
+ * Reconstruct the jets and fill the respective momenta.
+ */
+ virtual void reconstructJets(const ParticleVector&);
+
+ /**
+ * The jet finder to use
+ */
+ Ptr<JetFinder>::tptr jetFinder() const {
+ return theJetFinder;
+ }
+
+ /**
+ * The jet regions to match.
+ */
+ const vector<Ptr<JetRegion>::ptr>& jetRegions() const { return theJetRegions; }
+
+ /**
+ * Return the number of matched jets
+ */
+ unsigned int nJets() const { return theJets.size(); }
+
+ /**
+ * Set the momentum of the indicated jet.
+ */
+ LorentzMomentum& jetMomentum(const unsigned int id) {
+ return theJets[id];
+ }
+
+ /**
+ * Reconstruct all the variables for EW particles and fill the respective momenta.
+ */
+ virtual void reconstructEWParticles(ParticleVector&);
+
+ /**
+ * Set the momentum of the indicated leptonID.
+ */
+ LorentzMomentum& leptonIDMomentum(const unsigned int id) {
+ return theLeptonIDs[id];
+ }
+
+ /**
+ * Set the momentum of the indicated leptonPT.
+ */
+ LorentzMomentum& leptonPTMomentum(const unsigned int id) {
+ return theLeptonPTs[id];
+ }
+
+ /**
+ * Set the momentum of the indicated neutrino.
+ */
+ LorentzMomentum& neutrinoMomentum(const unsigned int id) {
+ return theNeutrinos[id];
+ }
+
+ /**
+ * Set the momentum of the indicated Higgs.
+ */
+ LorentzMomentum& higgsMomentum(const unsigned int id) {
+ return theHiggs[id];
+ }
+
+protected:
+
+ /**
+ * 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:
+
+ /**
+ * Collection of object histograms; ranges are adjusted to the
+ * maximum, so range constraints and rebinning can be applied later.
+ */
+ struct ObjectProperties {
+
+ /**
+ * Transverse momentum
+ */
+ Statistics::Histogram pt;
+ Statistics::Histogram ptlow;
+ Statistics::Histogram pt_logx;
+
+ /**
+ * Rapidity
+ */
+ Statistics::Histogram y;
+
+ /**
+ * Azimuth
+ */
+ Statistics::Histogram phi;
+
+ /**
+ * Mass
+ */
+ Statistics::Histogram mass;
+ Statistics::Histogram masslow;
+
+ /**
+ * Default constructor
+ */
+ ObjectProperties() {}
+
+ /**
+ * Construct given Ecm
+ */
+ ObjectProperties(const string& name, Energy)
+ : pt(name + "Pt",Statistics::Histogram::regularBinEdges(0,1000,200),true,false),
+ ptlow(name + "Ptlow",Statistics::Histogram::regularBinEdges(0,200,100),true,false),
+ pt_logx(name + "PtLogX",Statistics::Histogram::logBinEdges(0.1,1000,100),true,false),
+ y(name + "Y",Statistics::Histogram::regularBinEdges(-6,6,120),false,false),
+ phi(name + "Phi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,62),
+ make_pair(-Constants::pi,Constants::pi)),
+ mass(name + "Mass",Statistics::Histogram::regularBinEdges(0,5000,500),true,false),
+ masslow(name + "Masslow",Statistics::Histogram::regularBinEdges(0,250,100),true,false) {}
+
+ /**
+ * Count given momentum, weight and id
+ */
+ void count(const LorentzMomentum& p, double weight, unsigned int id) {
+ pt.count(Statistics::EventContribution(p.perp()/GeV,weight,5.),id);
+ ptlow.count(Statistics::EventContribution(p.perp()/GeV,weight,2.),id);
+ pt_logx.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
+ y.count(Statistics::EventContribution(p.rapidity(),weight,0.1),id);
+ phi.count(Statistics::EventContribution(p.phi(),weight,0.1),id);
+ mass.count(Statistics::EventContribution(p.m()/GeV,weight,10.),id);
+ masslow.count(Statistics::EventContribution(p.m()/GeV,weight,2.5),id);
+ }
+
+ /**
+ * Count given momentum components, weight and id
+ */
+ void count(Energy perp, double rapidity,
+ double xphi, Energy m,
+ double weight, unsigned int id) {
+ pt.count(Statistics::EventContribution(perp/GeV,weight,5.),id);
+ ptlow.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
+ pt_logx.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
+ y.count(Statistics::EventContribution(rapidity,weight,0.1),id);
+ phi.count(Statistics::EventContribution(xphi,weight,0.1),id);
+ mass.count(Statistics::EventContribution(m/GeV,weight,5.),id);
+ masslow.count(Statistics::EventContribution(m/GeV,weight,1.25),id);
+ }
+
+ /**
+ * Convert to XML
+ */
+ void finalize(XML::Element& elem) {
+ pt.finalize(); elem.append(pt.toXML());
+ ptlow.finalize(); elem.append(ptlow.toXML());
+ pt_logx.finalize(); elem.append(pt_logx.toXML());
+ y.finalize(); elem.append(y.toXML());
+ phi.finalize(); elem.append(phi.toXML());
+ mass.finalize(); elem.append(mass.toXML());
+ masslow.finalize(); elem.append(masslow.toXML());
+ }
+
+ };
+
+ /**
+ * Collection of pair histograms; ranges are adjusted to the
+ * maximum, so range constraints and rebinning can be applied later.
+ */
+ struct PairProperties
+ : public ObjectProperties {
+
+ /**
+ * Calculate deltaPhi
+ */
+ static double dPhi(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;
+ }
+
+ /**
+ * Calculate deltaY
+ */
+ static double dY(const LorentzMomentum& a,
+ const LorentzMomentum& b){
+ return abs(a.rapidity()-b.rapidity());
+ }
+
+ /**
+ * Calculate deltaR
+ */
+ static double dR(const LorentzMomentum& a,
+ const LorentzMomentum& b){
+ return sqrt(sqr(dPhi(a,b))+sqr(dY(a,b)));
+ }
+
+ /**
+ * Calculate ydoty
+ */
+ static double yy(const LorentzMomentum& a,
+ const LorentzMomentum& b){
+ double ya = a.rapidity();
+ double yb = b.rapidity();
+ double yres = sqrt(abs(ya*yb));
+ return ya*yb < 0. ? -yres : yres;
+ }
+
+ /**
+ * Delta y
+ */
+ Statistics::Histogram deltaY;
+
+ /**
+ * Delta phi
+ */
+ Statistics::Histogram deltaPhi;
+
+ /**
+ * Delta phi
+ */
+ Statistics::Histogram deltaR;
+
+ /**
+ * Product of the rapidities
+ */
+ Statistics::Histogram yDotY;
+
+ /**
+ * Default constructor
+ */
+ PairProperties()
+ : ObjectProperties() {}
+
+ /**
+ * Construct given Ecm
+ */
+ PairProperties(const string& name, Energy ecm)
+ : ObjectProperties(name,ecm),
+ deltaY(name + "DeltaY",Statistics::Histogram::regularBinEdges(0,6,60),true,false),
+ deltaPhi(name + "DeltaPhi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
+ make_pair(-Constants::pi,Constants::pi)),
+ deltaR(name + "DeltaR",Statistics::Histogram::regularBinEdges(0,10,100),true,false),
+ yDotY(name + "YDotY",Statistics::Histogram::regularBinEdges(-6,6,120),false,false) {}
+
+ /**
+ * Count given momentum, weight and id
+ */
+ void count(const LorentzMomentum& p, const LorentzMomentum& q, double weight, unsigned int id) {
+ ObjectProperties::count(p+q,weight,id);
+ deltaY.count(Statistics::EventContribution(dY(p,q),weight,0.1),id);
+ deltaPhi.count(Statistics::EventContribution(dPhi(p,q),weight,0.1),id);
+ deltaR.count(Statistics::EventContribution(dR(p,q),weight,0.1),id);
+ yDotY.count(Statistics::EventContribution(yy(p,q),weight,0.1),id);
+ }
+
+ /**
+ * Convert to XML
+ */
+ void finalize(XML::Element& elem) {
+ ObjectProperties::finalize(elem);
+ deltaY.finalize(); elem.append(deltaY.toXML());
+ deltaPhi.finalize(); elem.append(deltaPhi.toXML());
+ deltaR.finalize(); elem.append(deltaR.toXML());
+ yDotY.finalize(); elem.append(yDotY.toXML());
+ }
+
+ };
+
+ /**
+ * Collection of triple histograms; ranges are adjusted to the
+ * maximum, so range constraints and rebinning can be applied later.
+ */
+ struct TripleProperties
+ : public ObjectProperties {
+
+ /**
+ * Calculate deltaY^*
+ */
+ static double dYstar(const LorentzMomentum& a,
+ const LorentzMomentum& b,
+ const LorentzMomentum& c){
+ return c.rapidity()-(a.rapidity()+b.rapidity())/2.;
+ }
+
+ /**
+ * Delta y^*
+ */
+ Statistics::Histogram deltaYstar;
+
+ /**
+ * Default constructor
+ */
+ TripleProperties()
+ : ObjectProperties() {}
+
+ /**
+ * Construct given Ecm
+ */
+ TripleProperties(const string& name, Energy ecm)
+ : ObjectProperties(name,ecm),
+ deltaYstar(name + "DeltaYstar",Statistics::Histogram::regularBinEdges(-6,6,120),true,false) {}
+
+ /**
+ * Count given momentum, weight and id
+ */
+ void count(const LorentzMomentum& p, const LorentzMomentum& q, const LorentzMomentum& r,
+ double weight, unsigned int id) {
+ ObjectProperties::count(p+q+r,weight,id);
+ deltaYstar.count(Statistics::EventContribution(dYstar(p,q,r),weight,0.1),id);
+ }
+
+ /**
+ * Convert to XML
+ */
+ void finalize(XML::Element& elem) {
+ ObjectProperties::finalize(elem);
+ deltaYstar.finalize(); elem.append(deltaYstar.toXML());
+ }
+
+ };
+
+private:
+
+ /**
+ * Switch between fixed order and showered
+ */
+ bool theIsShowered;
+
+ /**
+ * Switch whether to apply extra analysis cuts
+ */
+ bool theApplyCuts;
+
+ /**
+ * The jet finder to use
+ */
+ Ptr<JetFinder>::ptr theJetFinder;
+
+ /**
+ * The jet regions to match.
+ */
+ vector<Ptr<JetRegion>::ptr> theJetRegions;
+
+ /**
+ * The reconstructed jets
+ */
+ map<unsigned int,LorentzMomentum> theJets;
+
+ /**
+ * The reconstructed leptons -- all ordered by ID
+ */
+ map<unsigned int,LorentzMomentum> theLeptonIDs;
+
+ /**
+ * The reconstructed leptons --charged ordered by PT
+ */
+ map<unsigned int,LorentzMomentum> theLeptonPTs;
+
+ /**
+ * The reconstructed neutrinos
+ */
+ map<unsigned int,LorentzMomentum> theNeutrinos;
+
+ /**
+ * The reconstructed Higgs
+ */
+ map<unsigned int,LorentzMomentum> theHiggs;
+
+ /**
+ * Jet properties
+ */
+ map<unsigned int,ObjectProperties> theJetProperties;
+
+ /**
+ * Exclusive jet properties
+ */
+ map<unsigned int,ObjectProperties> theExclusiveJetProperties;
+
+ /**
+ * Jet-inclusive properties
+ */
+ ObjectProperties theJetInclusiveProperties;
+
+ /**
+ * Jet-summed properties
+ */
+ ObjectProperties theJetSummedProperties;
+
+ /**
+ * Jet-average properties
+ */
+ ObjectProperties theJetAverageProperties;
+
+ /**
+ * Inclusive jet multiplicities
+ */
+ Statistics::Histogram theNJetsInclusive;
+
+ /**
+ * Exclusive jet multiplicities
+ */
+ Statistics::Histogram theNJetsExclusive;
+
+ /**
+ * Lepton properties -- all sorted by ID
+ */
+ map<unsigned int,ObjectProperties> theLeptonIDProperties;
+
+ /**
+ * Lepton properties -- charged sorted by PT
+ */
+ map<unsigned int,ObjectProperties> theLeptonPTProperties;
+
+ /**
+ * Neutrino properties
+ */
+ map<unsigned int,ObjectProperties> theNeutrinoProperties;
+
+ /**
+ * Higgs properties
+ */
+ map<unsigned int,ObjectProperties> theHiggsProperties;
+
+ /**
+ * Jet pair properties
+ */
+ map<pair<unsigned int,unsigned int>,PairProperties> theJetPairProperties;
+
+ /**
+ * Jet/lepton(all sorted by ID) pair properties
+ */
+ map<pair<unsigned int,unsigned int>,PairProperties> theJetLeptonIDPairProperties;
+
+ /**
+ * Jet/lepton(charged sorted by PT) pair properties
+ */
+ map<pair<unsigned int,unsigned int>,PairProperties> theJetLeptonPTPairProperties;
+
+ /**
+ * Jet/neutrino pair properties
+ */
+ map<pair<unsigned int,unsigned int>,PairProperties> theJetNeutrinoPairProperties;
+
+ /**
+ * Jet/Higgs pair properties
+ */
+ map<pair<unsigned int,unsigned int>,PairProperties> theJetHiggsPairProperties;
+
+ /**
+ * Lepton pair properties -- all sorted by ID
+ */
+ map<pair<unsigned int,unsigned int>,PairProperties> theLeptonIDPairProperties;
+
+ /**
+ * Lepton pair properties -- charged sorted by PT
+ */
+ map<pair<unsigned int,unsigned int>,PairProperties> theLeptonPTPairProperties;
+
+ /**
+ * Trijet properties
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeJetProperties;
+
+ /**
+ * Jet-pair/lepton(all sorted by ID) triple properties
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairLeptonIDTripleProperties;
+
+ /**
+ * Jet-pair/lepton(charged sorted by PT) triple properties
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairLeptonPTTripleProperties;
+
+ /**
+ * Jet-pair/neutrino triple properties
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairNeutrinoTripleProperties;
+
+ /**
+ * Jet-pair/Higgs triple properties
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairHiggsTripleProperties;
+
+ /**
+ * Trilepton properties -- all sorted by ID
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeLeptonIDProperties;
+
+ /**
+ * Trilepton properties -- charged sorted by PT
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeLeptonPTProperties;
+
+ /**
+ * Fourjet properties
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourJetProperties;
+
+ /**
+ * Fourlepton properties -- all sorted by ID
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourLeptonIDProperties;
+
+ /**
+ * Fourlepton properties -- charged sorted by PT
+ */
+ map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourLeptonPTProperties;
+
+protected:
+
+ /**
+ * Jet properties
+ */
+ ObjectProperties& jetProperties(const unsigned int id) {
+ map<unsigned int,ObjectProperties>::iterator h =
+ theJetProperties.find(id);
+ if ( h != theJetProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Jet" << id;
+ return
+ theJetProperties[id] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Exclusive jet properties
+ */
+ ObjectProperties& exclusiveJetProperties(const unsigned int id) {
+ map<unsigned int,ObjectProperties>::iterator h =
+ theExclusiveJetProperties.find(id);
+ if ( h != theExclusiveJetProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "ExclusiveJet" << id;
+ return
+ theExclusiveJetProperties[id] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet-inclusive properties
+ */
+ ObjectProperties& jetInclusiveProperties() {
+ if ( !theJetInclusiveProperties.pt.bins().empty() )
+ return theJetInclusiveProperties;
+ return
+ theJetInclusiveProperties =
+ ObjectProperties("JetInclusive",generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet-summed properties
+ */
+ ObjectProperties& jetSummedProperties() {
+ if ( !theJetSummedProperties.pt.bins().empty() )
+ return theJetSummedProperties;
+ return
+ theJetSummedProperties =
+ ObjectProperties("JetSummed",generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet-average properties
+ */
+ ObjectProperties& jetAverageProperties() {
+ if ( !theJetAverageProperties.pt.bins().empty() )
+ return theJetAverageProperties;
+ return
+ theJetAverageProperties =
+ ObjectProperties("JetAverage",generator()->maximumCMEnergy());
+ }
+
+
+ /**
+ * Inclusive jet multiplicities
+ */
+ Statistics::Histogram& nJetsInclusive() {
+ if ( !theNJetsInclusive.bins().empty() )
+ return theNJetsInclusive;
+ return
+ theNJetsInclusive =
+ Statistics::Histogram("NJetsInclusive",
+ Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
+ theJetRegions.size()+1),
+ true,true);
+ }
+
+ /**
+ * Exclusive jet multiplicities
+ */
+ Statistics::Histogram& nJetsExclusive() {
+ if ( !theNJetsExclusive.bins().empty() )
+ return theNJetsExclusive;
+ return
+ theNJetsExclusive =
+ Statistics::Histogram("NJetsExclusive",
+ Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
+ theJetRegions.size()+1),
+ true,true);
+ }
+
+ /**
+ * Lepton properties -- all sorted by ID
+ */
+ ObjectProperties& leptonIDProperties(const unsigned int id) {
+ map<unsigned int,ObjectProperties>::iterator h =
+ theLeptonIDProperties.find(id);
+ if ( h != theLeptonIDProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "LeptonID" << id;
+ return
+ theLeptonIDProperties[id] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Lepton properties -- charged sorted by PT
+ */
+ ObjectProperties& leptonPTProperties(const unsigned int id) {
+ map<unsigned int,ObjectProperties>::iterator h =
+ theLeptonPTProperties.find(id);
+ if ( h != theLeptonPTProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "LeptonPT" << id;
+ return
+ theLeptonPTProperties[id] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Neutrino properties
+ */
+ ObjectProperties& neutrinoProperties(const unsigned int id) {
+ map<unsigned int,ObjectProperties>::iterator h =
+ theNeutrinoProperties.find(id);
+ if ( h != theNeutrinoProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Neutrino" << id;
+ return
+ theNeutrinoProperties[id] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Higgs properties
+ */
+ ObjectProperties& higgsProperties(const unsigned int id) {
+ map<unsigned int,ObjectProperties>::iterator h =
+ theHiggsProperties.find(id);
+ if ( h != theHiggsProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Higgs" << id;
+ return
+ theHiggsProperties[id] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet pair properties
+ */
+ PairProperties& jetPairProperties(const unsigned int id, const unsigned int jd) {
+ map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
+ theJetPairProperties.find(make_pair(id,jd));
+ if ( h != theJetPairProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Jet" << id << jd;
+ return theJetPairProperties[make_pair(id,jd)] =
+ PairProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet/lepton(all sorted by ID) pair properties
+ */
+ PairProperties& jetLeptonIDPairProperties(const unsigned int id, const unsigned int jd) {
+ map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
+ theJetLeptonIDPairProperties.find(make_pair(id,jd));
+ if ( h != theJetLeptonIDPairProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Jet" << id << "LeptonID" << jd;
+ return theJetLeptonIDPairProperties[make_pair(id,jd)] =
+ PairProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet/lepton(charged sorted by PT) pair properties
+ */
+ PairProperties& jetLeptonPTPairProperties(const unsigned int id, const unsigned int jd) {
+ map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
+ theJetLeptonPTPairProperties.find(make_pair(id,jd));
+ if ( h != theJetLeptonPTPairProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Jet" << id << "LeptonPT" << jd;
+ return theJetLeptonPTPairProperties[make_pair(id,jd)] =
+ PairProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet/neutrino pair properties
+ */
+ PairProperties& jetNeutrinoPairProperties(const unsigned int id, const unsigned int jd) {
+ map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
+ theJetNeutrinoPairProperties.find(make_pair(id,jd));
+ if ( h != theJetNeutrinoPairProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Jet" << id << "Neutrino" << jd;
+ return theJetNeutrinoPairProperties[make_pair(id,jd)] =
+ PairProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet/Higgs pair properties
+ */
+ PairProperties& jetHiggsPairProperties(const unsigned int id, const unsigned int jd) {
+ map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
+ theJetHiggsPairProperties.find(make_pair(id,jd));
+ if ( h != theJetHiggsPairProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "Jet" << id << "Higgs" << jd;
+ return theJetHiggsPairProperties[make_pair(id,jd)] =
+ PairProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Lepton pair properties -- all sorted by ID
+ */
+ PairProperties& leptonIDPairProperties(const unsigned int id, const unsigned int jd) {
+ map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
+ theLeptonIDPairProperties.find(make_pair(id,jd));
+ if ( h != theLeptonIDPairProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "LeptonID" << id << jd;
+ return theLeptonIDPairProperties[make_pair(id,jd)] =
+ PairProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Lepton pair properties -- charged sorted by PT
+ */
+ PairProperties& leptonPTPairProperties(const unsigned int id, const unsigned int jd) {
+ map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
+ theLeptonPTPairProperties.find(make_pair(id,jd));
+ if ( h != theLeptonPTPairProperties.end() )
+ return h->second;
+ ostringstream ids; ids << "LeptonPT" << id << jd;
+ return theLeptonPTPairProperties[make_pair(id,jd)] =
+ PairProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Trijet properties
+ */
+ TripleProperties& threeJetProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
+ theThreeJetProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
+ if ( it != theThreeJetProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "Jet" << id1 << id2 << id3;
+ return theThreeJetProperties[boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
+ TripleProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet-pair/lepton(all sorted by ID) triple properties
+ */
+ TripleProperties& jetPairLeptonIDTripleProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
+ theJetPairLeptonIDTripleProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
+ if ( it != theJetPairLeptonIDTripleProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "Jet" << id1 << id2 << "LeptonID" << id3;
+ return theJetPairLeptonIDTripleProperties[boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
+ TripleProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet-pair/lepton(charged sorted by PT) triple properties
+ */
+ TripleProperties& jetPairLeptonPTTripleProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
+ theJetPairLeptonPTTripleProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
+ if ( it != theJetPairLeptonPTTripleProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "Jet" << id1 << id2 << "LeptonPT" << id3;
+ return theJetPairLeptonPTTripleProperties[boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
+ TripleProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet-pair/neutrino triple properties
+ */
+ TripleProperties& jetPairNeutrinoTripleProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
+ theJetPairNeutrinoTripleProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
+ if ( it != theJetPairNeutrinoTripleProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "Jet" << id1 << id2 << "Neutrino" << id3;
+ return theJetPairNeutrinoTripleProperties[boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
+ TripleProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Jet-pair/Higgs triple properties
+ */
+ TripleProperties& jetPairHiggsTripleProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
+ theJetPairHiggsTripleProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
+ if ( it != theJetPairHiggsTripleProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "Jet" << id1 << id2 << "Higgs" << id3;
+ return theJetPairHiggsTripleProperties[boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
+ TripleProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Trilepton properties -- all sorted by ID
+ */
+ TripleProperties& threeLeptonIDProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
+ theThreeLeptonIDProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
+ if ( it != theThreeLeptonIDProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "LeptonID" << id1 << id2 << id3;
+ return theThreeLeptonIDProperties[boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
+ TripleProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Trilepton properties -- charged sorted by PT
+ */
+ TripleProperties& threeLeptonPTProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
+ theThreeLeptonPTProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
+ if ( it != theThreeLeptonPTProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "LeptonPT" << id1 << id2 << id3;
+ return theThreeLeptonPTProperties[boost::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
+ TripleProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Fourjet properties
+ */
+ ObjectProperties& fourJetProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3, const unsigned int id4) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
+ theFourJetProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
+ if ( it != theFourJetProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "Jet" << id1 << id2 << id3 << id4;
+ return theFourJetProperties[boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Fourlepton properties -- all sorted by ID
+ */
+ ObjectProperties& fourLeptonIDProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3, const unsigned int id4) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
+ theFourLeptonIDProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
+ if ( it != theFourLeptonIDProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "LeptonID" << id1 << id2 << id3 << id4;
+ return theFourLeptonIDProperties[boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Fourlepton properties -- charged sorted by PT
+ */
+ ObjectProperties& fourLeptonPTProperties(const unsigned int id1, const unsigned int id2,
+ const unsigned int id3, const unsigned int id4) {
+ map<boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
+ theFourLeptonPTProperties.find(boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
+ if ( it != theFourLeptonPTProperties.end() )
+ return it->second;
+ ostringstream ids;
+ ids << "LeptonPT" << id1 << id2 << id3 << id4;
+ return theFourLeptonPTProperties[boost::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
+ ObjectProperties(ids.str(),generator()->maximumCMEnergy());
+ }
+
+ /**
+ * Perform any additional analysis required
+ */
+ virtual void analyzeSpecial(long, double) {}
+
+ /**
+ * Append any additional histograms to the given histogram element
+ */
+ virtual void finalize(XML::Element&) {}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ LeptonsJetsAnalysis & operator=(const LeptonsJetsAnalysis &);
+
+};
+
+}
+
+#endif /* Herwig_LeptonsJetsAnalysis_H */
diff --git a/Analysis/Makefile.am b/Analysis/Makefile.am
--- a/Analysis/Makefile.am
+++ b/Analysis/Makefile.am
@@ -1,58 +1,59 @@
pkglib_LTLIBRARIES = HwAnalysis.la
HwAnalysis_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 9:3:0
HwAnalysis_la_SOURCES = \
EventShapes.cc EventShapes.fh EventShapes.h \
EventShapesMasterAnalysis.cc EventShapesMasterAnalysis.h \
BasicConsistency.cc BasicConsistency.h \
LEPMultiplicityCount.cc LEPMultiplicityCount.h \
MultiplicityInfo.h \
LEPBMultiplicity.cc\
LEPBMultiplicity.h \
SimpleLHCAnalysis.h SimpleLHCAnalysis.cc\
TTbarAnalysis.h TTbarAnalysis.cc\
LPairAnalysis.h LPairAnalysis.cc\
GammaGammaAnalysis.h GammaGammaAnalysis.cc\
GammaJetAnalysis.h GammaJetAnalysis.cc\
HiggsJetAnalysis.h HiggsJetAnalysis.cc \
ParallelRunAnalysis.h ParallelRunAnalysis.cc \
DrellYanPT.h DrellYanPT.cc
pkglib_LTLIBRARIES += HwJetsAnalysis.la
HwJetsAnalysis_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
HwJetsAnalysis_la_SOURCES = \
+LeptonsJetsAnalysis.h LeptonsJetsAnalysis.cc \
JetsPlusAnalysis.h JetsPlusAnalysis.cc \
HJetsAnalysis.h HJetsAnalysis.cc \
ZJetsAnalysis.h ZJetsAnalysis.cc \
TTJetsAnalysis.h TTJetsAnalysis.cc \
CrossSectionAnalysis.h CrossSectionAnalysis.cc
pkglib_LTLIBRARIES += HwLEPAnalysis.la
HwLEPAnalysis_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 6:0:0
HwLEPAnalysis_la_SOURCES = \
BFragmentationAnalysisHandler.cc BFragmentationAnalysisHandler.h\
SingleParticleAnalysis.cc SingleParticleAnalysis.h\
LEPEventShapes.cc LEPEventShapes.h\
IdentifiedParticleAnalysis.cc IdentifiedParticleAnalysis.h\
BELLECharmAnalysis.h BELLECharmAnalysis.cc\
CLEOCharmAnalysis.h CLEOCharmAnalysis.cc
# analysis code which depends on fastjet
if WANT_LIBFASTJET
pkglib_LTLIBRARIES += HwLEPJetAnalysis.la
HwLEPJetAnalysis_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 8:0:0
HwLEPJetAnalysis_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
HwLEPJetAnalysis_la_LIBADD = $(FASTJETLIBS)
HwLEPJetAnalysis_la_SOURCES = \
LEPJetAnalysis.cc LEPJetAnalysis.h\
LEPFourJetsAnalysis.cc LEPFourJetsAnalysis.h
endif
pkglib_LTLIBRARIES += HwTevatronAnalysis.la
HwTevatronAnalysis_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
HwTevatronAnalysis_la_CPPFLAGS = $(AM_CPPFLAGS)
HwTevatronAnalysis_la_SOURCES = \
ZpTRun2.cc ZpTRun2.h \
ZpTRun1.cc ZpTRun1.h \
Wpt.cc Wpt.h \
Zrapidity.cc Zrapidity.h

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:05 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4212337
Default Alt Text
(57 KB)

Event Timeline