LeptonsJetsAnalysis.h

// -*- 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"
 
namespace Herwig {
 
using namespace ThePEG;
 
class LeptonsJetsAnalysis: public AnalysisHandler {
 
public:
 
  LeptonsJetsAnalysis();
 
public:
 
  virtual void analyze(tEventPtr event, long ieve, int loop, int state);
 
protected:
 
  void analyze(ParticleVector&, long, double);
 
  void clear() {
    theJets.clear();
    theEWIDs.clear();
    theChargedLeptons.clear();
    theNeutrinos.clear();
    theHiggs.clear();
  }
 
  virtual void reconstructJets(const ParticleVector&);
 
  Ptr<JetFinder>::tptr jetFinder() const {
    return theJetFinder;
  }
 
  const vector<Ptr<JetRegion>::ptr>& jetRegions() const { return theJetRegions; }
 
  unsigned int nJets() const { return theJets.size(); }
 
  LorentzMomentum& jetMomentum(const unsigned int id) {
    return theJets[id];
  }
 
  virtual void reconstructEWParticles(ParticleVector&);
 
  LorentzMomentum& eWIDMomentum(const unsigned int id) {
    return theEWIDs[id];
  }
 
  LorentzMomentum& chargedLeptonMomentum(const unsigned int id) {
    return theChargedLeptons[id];
  }
 
  LorentzMomentum& neutrinoMomentum(const unsigned int id) {
    return theNeutrinos[id];
  }
 
  LorentzMomentum& pTmissMomentum() {
    return thePTmiss;
  }
 
  LorentzMomentum& higgsMomentum(const unsigned int id) {
    return theHiggs[id];
  }
 
protected:
 
  virtual void dofinish();
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
 
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
protected:
 
  struct ObjectProperties {
 
    Statistics::Histogram pt;
    Statistics::Histogram ptlow;
    Statistics::Histogram pt_logx;
 
    Statistics::Histogram y;
 
    Statistics::Histogram phi;
 
    Statistics::Histogram mass;
    Statistics::Histogram masslow;
 
    ObjectProperties() {}
 
    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) {}
 
    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);
    }
 
    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);
    }
 
    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());
    }
 
  };
 
  struct PairProperties
    : public ObjectProperties {
 
    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;
    }
 
    static double dY(const LorentzMomentum& a,
                     const LorentzMomentum& b){
      return abs(a.rapidity()-b.rapidity());
    }
 
    static double dR(const LorentzMomentum& a,
                     const LorentzMomentum& b){
      return sqrt(sqr(dPhi(a,b))+sqr(dY(a,b)));
    }
 
    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;
    }
 
    Statistics::Histogram deltaY;
 
    Statistics::Histogram deltaPhi;
 
    Statistics::Histogram deltaR;
 
    Statistics::Histogram yDotY;
 
    PairProperties() 
      : ObjectProperties() {}
 
    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) {}
 
    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);
    }
 
    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());
    }
 
  };
 
  struct TripleProperties
    : public ObjectProperties {
 
    static double dYstar(const LorentzMomentum& a,
                         const LorentzMomentum& b,
                         const LorentzMomentum& c){
      return c.rapidity()-(a.rapidity()+b.rapidity())/2.;
    }
 
    static double dZstar(const LorentzMomentum& a,
                         const LorentzMomentum& b,
                         const LorentzMomentum& c){
      return dYstar(a,b,c)*2./abs(a.rapidity()-b.rapidity());
    }
 
    Statistics::Histogram deltaYstar;
 
    Statistics::Histogram deltaZstar;
 
    TripleProperties() 
      : ObjectProperties() {}
 
    TripleProperties(const string& name, Energy ecm)
      : ObjectProperties(name,ecm),
        deltaYstar(name + "DeltaYstar",Statistics::Histogram::regularBinEdges(-6,6,120),true,false),
        deltaZstar(name + "DeltaZstar",Statistics::Histogram::regularBinEdges(-3,3,120),true,false) {}
 
    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);
      deltaZstar.count(Statistics::EventContribution(dZstar(p,q,r),weight,0.05),id);
    }
 
    void finalize(XML::Element& elem) {
      ObjectProperties::finalize(elem);
      deltaYstar.finalize(); elem.append(deltaYstar.toXML());
      deltaZstar.finalize(); elem.append(deltaZstar.toXML());
    }
 
  };
 
private:
 
  bool theIsShowered;
 
  bool theApplyCuts;
 
  Ptr<JetFinder>::ptr theJetFinder;
 
  vector<Ptr<JetRegion>::ptr> theJetRegions;
 
  map<unsigned int,LorentzMomentum> theJets;
 
  map<unsigned int,LorentzMomentum> theEWIDs;
 
  map<unsigned int,LorentzMomentum> theChargedLeptons;
 
  map<unsigned int,LorentzMomentum> theNeutrinos;
 
  LorentzMomentum thePTmiss;
 
  map<unsigned int,LorentzMomentum> theHiggs;
 
  map<unsigned int,ObjectProperties> theJetProperties;
 
  map<unsigned int,ObjectProperties> theExclusiveJetProperties;
 
  ObjectProperties theJetInclusiveProperties;
 
  ObjectProperties theJetSummedProperties;
 
  ObjectProperties theJetAverageProperties;
 
  Statistics::Histogram theNJetsInclusive;
 
  Statistics::Histogram theNJetsExclusive;
 
  map<unsigned int,ObjectProperties> theEWIDProperties;
 
  map<unsigned int,ObjectProperties> theChargedLeptonProperties;
 
  map<unsigned int,ObjectProperties> theNeutrinoProperties;
 
  ObjectProperties thePTmissProperties;
 
  map<unsigned int,ObjectProperties> theHiggsProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theJetPairProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theJetEWIDPairProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theJetChargedLeptonPairProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theJetNeutrinoPairProperties;
 
  map<unsigned int,PairProperties> theJetPTmissPairProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theJetHiggsPairProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theEWIDPairProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theChargedLeptonPairProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeJetProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairEWIDTripleProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairChargedLeptonTripleProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairNeutrinoTripleProperties;
 
  map<pair<unsigned int,unsigned int>,TripleProperties> theJetPairPTmissTripleProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairHiggsTripleProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeEWIDProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeChargedLeptonProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourJetProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourEWIDProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourChargedLeptonProperties;
 
protected:
 
  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());
  }
 
  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());
  }
 
  ObjectProperties& jetInclusiveProperties() {
    if ( !theJetInclusiveProperties.pt.bins().empty() )
      return theJetInclusiveProperties;
    return
      theJetInclusiveProperties = 
      ObjectProperties("JetInclusive",generator()->maximumCMEnergy());
  }
 
  ObjectProperties& jetSummedProperties() {
    if ( !theJetSummedProperties.pt.bins().empty() )
      return theJetSummedProperties;
    return
      theJetSummedProperties = 
      ObjectProperties("JetSummed",generator()->maximumCMEnergy());
  }
 
  ObjectProperties& jetAverageProperties() {
    if ( !theJetAverageProperties.pt.bins().empty() )
      return theJetAverageProperties;
    return
      theJetAverageProperties = 
      ObjectProperties("JetAverage",generator()->maximumCMEnergy());
  }
 
 
  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);
  }
 
  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);
  }
 
  ObjectProperties& eWIDProperties(const unsigned int id) {
    map<unsigned int,ObjectProperties>::iterator h = 
      theEWIDProperties.find(id);
    if ( h != theEWIDProperties.end() )
      return h->second;
    ostringstream ids; ids << "EWID" << id;
    return 
      theEWIDProperties[id] = 
      ObjectProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  ObjectProperties& chargedLeptonProperties(const unsigned int id) {
    map<unsigned int,ObjectProperties>::iterator h = 
      theChargedLeptonProperties.find(id);
    if ( h != theChargedLeptonProperties.end() )
      return h->second;
    ostringstream ids; ids << "ChargedLepton" << id;
    return 
      theChargedLeptonProperties[id] = 
      ObjectProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  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());
  }
 
  ObjectProperties& pTmissProperties() {
    if ( !thePTmissProperties.pt.bins().empty() )
      return thePTmissProperties;
    return
      thePTmissProperties = 
      ObjectProperties("PTmiss",generator()->maximumCMEnergy());
  }
 
  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());
  }
 
  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());
  }
 
  PairProperties& jetEWIDPairProperties(const unsigned int id, const unsigned int jd) {
    map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = 
      theJetEWIDPairProperties.find(make_pair(id,jd));
    if ( h != theJetEWIDPairProperties.end() )
      return h->second;
    ostringstream ids; ids << "Jet" << id << "EWID" << jd;
    return theJetEWIDPairProperties[make_pair(id,jd)] = 
      PairProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  PairProperties& jetChargedLeptonPairProperties(const unsigned int id, const unsigned int jd) {
    map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = 
      theJetChargedLeptonPairProperties.find(make_pair(id,jd));
    if ( h != theJetChargedLeptonPairProperties.end() )
      return h->second;
    ostringstream ids; ids << "Jet" << id << "ChargedLepton" << jd;
    return theJetChargedLeptonPairProperties[make_pair(id,jd)] = 
      PairProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  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());
  }
 
  PairProperties& jetPTmissPairProperties(const unsigned int id) {
    map<unsigned int,PairProperties>::iterator h = 
      theJetPTmissPairProperties.find(id);
    if ( h != theJetPTmissPairProperties.end() )
      return h->second;
    ostringstream ids; ids << "Jet" << id << "PTmiss";
    return theJetPTmissPairProperties[id] = 
      PairProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  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());
  }
 
  PairProperties& eWIDPairProperties(const unsigned int id, const unsigned int jd) {
    map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = 
      theEWIDPairProperties.find(make_pair(id,jd));
    if ( h != theEWIDPairProperties.end() )
      return h->second;
    ostringstream ids; ids << "EWID" << id << jd;
    return theEWIDPairProperties[make_pair(id,jd)] = 
      PairProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  PairProperties& chargedLeptonPairProperties(const unsigned int id, const unsigned int jd) {
    map<pair<unsigned int,unsigned int>,PairProperties>::iterator h = 
      theChargedLeptonPairProperties.find(make_pair(id,jd));
    if ( h != theChargedLeptonPairProperties.end() )
      return h->second;
    ostringstream ids; ids << "ChargedLepton" << id << jd;
    return theChargedLeptonPairProperties[make_pair(id,jd)] = 
      PairProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& threeJetProperties(const unsigned int id1, const unsigned int id2,
                                       const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
      theThreeJetProperties.find(std::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[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& jetPairEWIDTripleProperties(const unsigned int id1, const unsigned int id2,
                                              const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
      theJetPairEWIDTripleProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
    if ( it != theJetPairEWIDTripleProperties.end() )
      return it->second;
    ostringstream ids; 
    ids << "Jet" << id1 << id2 << "EWID" << id3;
    return theJetPairEWIDTripleProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& jetPairChargedLeptonTripleProperties(const unsigned int id1, const unsigned int id2,
                                              const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
      theJetPairChargedLeptonTripleProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
    if ( it != theJetPairChargedLeptonTripleProperties.end() )
      return it->second;
    ostringstream ids; 
    ids << "Jet" << id1 << id2 << "ChargedLepton" << id3;
    return theJetPairChargedLeptonTripleProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& jetPairNeutrinoTripleProperties(const unsigned int id1, const unsigned int id2,
                                              const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
      theJetPairNeutrinoTripleProperties.find(std::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[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& jetPairPTmissTripleProperties(const unsigned int id1, const unsigned int id2) {
    map<pair<unsigned int,unsigned int>,TripleProperties>::iterator it =
      theJetPairPTmissTripleProperties.find(pair<unsigned int,unsigned int>(id1,id2));
    if ( it != theJetPairPTmissTripleProperties.end() )
      return it->second;
    ostringstream ids; 
    ids << "Jet" << id1 << id2 << "PTmiss";
    return theJetPairPTmissTripleProperties[pair<unsigned int,unsigned int>(id1,id2)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& jetPairHiggsTripleProperties(const unsigned int id1, const unsigned int id2,
                                              const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
      theJetPairHiggsTripleProperties.find(std::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[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& threeEWIDProperties(const unsigned int id1, const unsigned int id2,
                                        const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
      theThreeEWIDProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
    if ( it != theThreeEWIDProperties.end() )
      return it->second;
    ostringstream ids; 
    ids << "EWID" << id1 << id2 << id3;
    return theThreeEWIDProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  TripleProperties& threeChargedLeptonProperties(const unsigned int id1, const unsigned int id2,
                                       const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
      theThreeChargedLeptonProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
    if ( it != theThreeChargedLeptonProperties.end() )
      return it->second;
    ostringstream ids; 
    ids << "ChargedLepton" << id1 << id2 << id3;
    return theThreeChargedLeptonProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
      TripleProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  ObjectProperties& fourJetProperties(const unsigned int id1, const unsigned int id2,
                                      const unsigned int id3, const unsigned int id4) {
    map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
      theFourJetProperties.find(std::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[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
      ObjectProperties(ids.str(),generator()->maximumCMEnergy());
  }  
 
  ObjectProperties& fourEWIDProperties(const unsigned int id1, const unsigned int id2,
                                      const unsigned int id3, const unsigned int id4) {
    map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
      theFourEWIDProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
    if ( it != theFourEWIDProperties.end() )
      return it->second;
    ostringstream ids; 
    ids << "EWID" << id1 << id2 << id3 << id4;
    return theFourEWIDProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
      ObjectProperties(ids.str(),generator()->maximumCMEnergy());
  }  
 
  ObjectProperties& fourChargedLeptonProperties(const unsigned int id1, const unsigned int id2,
                                      const unsigned int id3, const unsigned int id4) {
    map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
      theFourChargedLeptonProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
    if ( it != theFourChargedLeptonProperties.end() )
      return it->second;
    ostringstream ids; 
    ids << "ChargedLepton" << id1 << id2 << id3 << id4;
    return theFourChargedLeptonProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
      ObjectProperties(ids.str(),generator()->maximumCMEnergy());
  }  
 
  virtual void analyzeSpecial(long, double) {}
 
  virtual void finalize(XML::Element&) {}
 
private:
 
  LeptonsJetsAnalysis & operator=(const LeptonsJetsAnalysis &) = delete;
 
};
 
}
 
#endif /* Herwig_LeptonsJetsAnalysis_H */