Herwig/Analysis/LeptonsJetsAnalysis.h

0001 // -*- C++ -*-
0002 #ifndef Herwig_LeptonsJetsAnalysis_H
0003 #define Herwig_LeptonsJetsAnalysis_H
0004 //
0005 // This is the declaration of the LeptonsJetsAnalysis class.
0006 //
0007
0008 #include "ThePEG/Handlers/AnalysisHandler.h"
0009 #include "ThePEG/Cuts/JetFinder.h"
0010 #include "ThePEG/Cuts/JetRegion.h"
0011 #include "ThePEG/Repository/EventGenerator.h"
0012 #include "Herwig/Utilities/Statistics/Histogram.h"
0013
0014 namespace Herwig {
0015
0016 using namespace ThePEG;
0017
0018 /**
0019  * Here is the documentation of the LeptonsJetsAnalysis class.
0020  *
0021  * @see \ref LeptonsJetsAnalysisInterfaces "The interfaces"
0022  * defined for LeptonsJetsAnalysis.
0023  */
0024 class LeptonsJetsAnalysis: public AnalysisHandler {
0025
0026 public:
0027
0028   /**
0029    * The default constructor.
0030    */
0031   LeptonsJetsAnalysis();
0032
0033 public:
0034
0035   /** @name Virtual functions required by the AnalysisHandler class. */
0036   //@{
0037   /**
0038    * Analyze a given Event. Note that a fully generated event
0039    * may be presented several times, if it has been manipulated in
0040    * between. The default version of this function will call transform
0041    * to make a lorentz transformation of the whole event, then extract
0042    * all final state particles and call analyze(tPVector) of this
0043    * analysis object and those of all associated analysis objects. The
0044    * default version will not, however, do anything on events which
0045    * have not been fully generated, or have been manipulated in any
0046    * way.
0047    * @param event pointer to the Event to be analyzed.
0048    * @param ieve the event number.
0049    * @param loop the number of times this event has been presented.
0050    * If negative the event is now fully generated.
0051    * @param state a number different from zero if the event has been
0052    * manipulated in some way since it was last presented.
0053    */
0054   virtual void analyze(tEventPtr event, long ieve, int loop, int state);
0055   //@}
0056
0057 protected:
0058
0059   /**
0060    * Analyze one subprocess, given the event number it belongs to
0061    */
0062   void analyze(ParticleVector&, long, double);
0063
0064   /**
0065    * Clear the electroweak objects and jets for the next event
0066    */
0067   void clear() {
0068     theJets.clear();
0069     theEWIDs.clear();
0070     theChargedLeptons.clear();
0071     theNeutrinos.clear();
0072     theHiggs.clear();
0073   }
0074
0075   /**
0076    * Reconstruct the jets and fill the respective momenta.
0077    */
0078   virtual void reconstructJets(const ParticleVector&);
0079
0080   /**
0081    * The jet finder to use
0082    */
0083   Ptr<JetFinder>::tptr jetFinder() const {
0084     return theJetFinder;
0085   }
0086
0087   /**
0088    * The jet regions to match.
0089    */
0090   const vector<Ptr<JetRegion>::ptr>& jetRegions() const { return theJetRegions; }
0091
0092   /**
0093    * Return the number of matched jets
0094    */
0095   unsigned int nJets() const { return theJets.size(); }
0096
0097   /**
0098    * Set the momentum of the indicated jet.
0099    */
0100   LorentzMomentum& jetMomentum(const unsigned int id) {
0101     return theJets[id];
0102   }
0103
0104   /**
0105    * Reconstruct all the variables for EW particles and fill the respective momenta.
0106    */
0107   virtual void reconstructEWParticles(ParticleVector&);
0108
0109   /**
0110    * Set the momentum of the indicated electroweak particle.
0111    */
0112   LorentzMomentum& eWIDMomentum(const unsigned int id) {
0113     return theEWIDs[id];
0114   }
0115
0116   /**
0117    * Set the momentum of the indicated charged lepton.
0118    */
0119   LorentzMomentum& chargedLeptonMomentum(const unsigned int id) {
0120     return theChargedLeptons[id];
0121   }
0122
0123   /**
0124    * Set the momentum of the indicated neutrino.
0125    */
0126   LorentzMomentum& neutrinoMomentum(const unsigned int id) {
0127     return theNeutrinos[id];
0128   }
0129
0130   /**
0131    * Set the missing pT momentum.
0132    */
0133   LorentzMomentum& pTmissMomentum() {
0134     return thePTmiss;
0135   }
0136
0137   /**
0138    * Set the momentum of the indicated Higgs.
0139    */
0140   LorentzMomentum& higgsMomentum(const unsigned int id) {
0141     return theHiggs[id];
0142   }
0143
0144 protected:
0145
0146   /**
0147    * Finalize this object. Called in the run phase just after a
0148    * run has ended. Used eg. to write out statistics.
0149    */
0150   virtual void dofinish();
0151
0152 public:
0153
0154   /** @name Functions used by the persistent I/O system. */
0155   //@{
0156   /**
0157    * Function used to write out object persistently.
0158    * @param os the persistent output stream written to.
0159    */
0160   void persistentOutput(PersistentOStream & os) const;
0161
0162   /**
0163    * Function used to read in object persistently.
0164    * @param is the persistent input stream read from.
0165    * @param version the version number of the object when written.
0166    */
0167   void persistentInput(PersistentIStream & is, int version);
0168   //@}
0169
0170   /**
0171    * The standard Init function used to initialize the interfaces.
0172    * Called exactly once for each class by the class description system
0173    * before the main function starts or
0174    * when this class is dynamically loaded.
0175    */
0176   static void Init();
0177
0178 protected:
0179
0180   /** @name Clone Methods. */
0181   //@{
0182   /**
0183    * Make a simple clone of this object.
0184    * @return a pointer to the new object.
0185    */
0186   virtual IBPtr clone() const;
0187
0188   /** Make a clone of this object, possibly modifying the cloned object
0189    * to make it sane.
0190    * @return a pointer to the new object.
0191    */
0192   virtual IBPtr fullclone() const;
0193   //@}
0194
0195
0196 // If needed, insert declarations of virtual function defined in the
0197 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0198
0199 protected:
0200
0201   /**
0202    * Collection of object histograms; ranges are adjusted to the
0203    * maximum, so range constraints and rebinning can be applied later.
0204    */
0205   struct ObjectProperties {
0206
0207     /**
0208      * Transverse momentum
0209      */
0210     Statistics::Histogram pt;
0211     Statistics::Histogram ptlow;
0212     Statistics::Histogram pt_logx;
0213
0214     /**
0215      * Rapidity
0216      */
0217     Statistics::Histogram y;
0218
0219     /**
0220      * Azimuth
0221      */
0222     Statistics::Histogram phi;
0223
0224     /**
0225      * Mass
0226      */
0227     Statistics::Histogram mass;
0228     Statistics::Histogram masslow;
0229
0230     /**
0231      * Default constructor
0232      */
0233     ObjectProperties() {}
0234
0235     /**
0236      * Construct given Ecm
0237      */
0238     ObjectProperties(const string& name, Energy)
0239       : pt(name + "Pt",Statistics::Histogram::regularBinEdges(0,1000,200),true,false),
0240         ptlow(name + "Ptlow",Statistics::Histogram::regularBinEdges(0,200,100),true,false),
0241         pt_logx(name + "PtLogX",Statistics::Histogram::logBinEdges(0.1,1000,100),true,false),
0242     y(name + "Y",Statistics::Histogram::regularBinEdges(-6,6,120),false,false),
0243     phi(name + "Phi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,62),
0244         make_pair(-Constants::pi,Constants::pi)),
0245     mass(name + "Mass",Statistics::Histogram::regularBinEdges(0,5000,500),true,false),
0246     masslow(name + "Masslow",Statistics::Histogram::regularBinEdges(0,250,100),true,false) {}
0247
0248     /**
0249      * Count given momentum, weight and id
0250      */
0251     void count(const LorentzMomentum& p, double weight, unsigned int id) {
0252       pt.count(Statistics::EventContribution(p.perp()/GeV,weight,5.),id);
0253       ptlow.count(Statistics::EventContribution(p.perp()/GeV,weight,2.),id);
0254       pt_logx.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
0255       y.count(Statistics::EventContribution(p.rapidity(),weight,0.1),id);
0256       phi.count(Statistics::EventContribution(p.phi(),weight,0.1),id);
0257       mass.count(Statistics::EventContribution(p.m()/GeV,weight,10.),id);
0258       masslow.count(Statistics::EventContribution(p.m()/GeV,weight,2.5),id);
0259     }
0260
0261     /**
0262      * Count given momentum components, weight and id
0263      */
0264     void count(Energy perp, double rapidity,
0265            double xphi, Energy m,
0266            double weight, unsigned int id) {
0267       pt.count(Statistics::EventContribution(perp/GeV,weight,5.),id);
0268       ptlow.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
0269       pt_logx.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
0270       y.count(Statistics::EventContribution(rapidity,weight,0.1),id);
0271       phi.count(Statistics::EventContribution(xphi,weight,0.1),id);
0272       mass.count(Statistics::EventContribution(m/GeV,weight,5.),id);
0273       masslow.count(Statistics::EventContribution(m/GeV,weight,1.25),id);
0274     }
0275
0276     /**
0277      * Convert to XML
0278      */
0279     void finalize(XML::Element& elem) {
0280       pt.finalize(); elem.append(pt.toXML());
0281       ptlow.finalize(); elem.append(ptlow.toXML());
0282       pt_logx.finalize(); elem.append(pt_logx.toXML());
0283       y.finalize(); elem.append(y.toXML());
0284       phi.finalize(); elem.append(phi.toXML());
0285       mass.finalize(); elem.append(mass.toXML());
0286       masslow.finalize(); elem.append(masslow.toXML());
0287     }
0288
0289   };
0290
0291   /**
0292    * Collection of pair histograms; ranges are adjusted to the
0293    * maximum, so range constraints and rebinning can be applied later.
0294    */
0295   struct PairProperties
0296     : public ObjectProperties {
0297
0298     /**
0299      * Calculate deltaPhi
0300      */
0301     static double dPhi(const LorentzMomentum& a,
0302                const LorentzMomentum& b){
0303       double phi1 = a.phi();
0304       double phi2 = b.phi();
0305       double diff=phi1-phi2;
0306       if(diff<-Constants::pi){
0307     diff+=(2.0*Constants::pi);
0308       }
0309       else if (diff>Constants::pi){
0310     diff-=(2.0*Constants::pi);
0311       }
0312       return diff;
0313     }
0314
0315     /**
0316      * Calculate deltaY
0317      */
0318     static double dY(const LorentzMomentum& a,
0319              const LorentzMomentum& b){
0320       return abs(a.rapidity()-b.rapidity());
0321     }
0322
0323     /**
0324      * Calculate deltaR
0325      */
0326     static double dR(const LorentzMomentum& a,
0327              const LorentzMomentum& b){
0328       return sqrt(sqr(dPhi(a,b))+sqr(dY(a,b)));
0329     }
0330
0331     /**
0332      * Calculate ydoty
0333      */
0334     static double yy(const LorentzMomentum& a,
0335              const LorentzMomentum& b){
0336       double ya = a.rapidity();
0337       double yb = b.rapidity();
0338       double yres = sqrt(abs(ya*yb));
0339       return ya*yb < 0. ? -yres : yres;
0340     }
0341
0342     /**
0343      * Delta y
0344      */
0345     Statistics::Histogram deltaY;
0346
0347     /**
0348      * Delta phi
0349      */
0350     Statistics::Histogram deltaPhi;
0351
0352     /**
0353      * Delta phi
0354      */
0355     Statistics::Histogram deltaR;
0356
0357     /**
0358      * Product of the rapidities
0359      */
0360     Statistics::Histogram yDotY;
0361
0362     /**
0363      * Default constructor
0364      */
0365     PairProperties()
0366       : ObjectProperties() {}
0367
0368     /**
0369      * Construct given Ecm
0370      */
0371     PairProperties(const string& name, Energy ecm)
0372       : ObjectProperties(name,ecm),
0373     deltaY(name + "DeltaY",Statistics::Histogram::regularBinEdges(0,6,60),true,false),
0374     deltaPhi(name + "DeltaPhi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
0375          make_pair(-Constants::pi,Constants::pi)),
0376     deltaR(name + "DeltaR",Statistics::Histogram::regularBinEdges(0,10,100),true,false),
0377     yDotY(name + "YDotY",Statistics::Histogram::regularBinEdges(-6,6,120),false,false) {}
0378
0379     /**
0380      * Count given momentum, weight and id
0381      */
0382     void count(const LorentzMomentum& p, const LorentzMomentum& q, double weight, unsigned int id) {
0383       ObjectProperties::count(p+q,weight,id);
0384       deltaY.count(Statistics::EventContribution(dY(p,q),weight,0.1),id);
0385       deltaPhi.count(Statistics::EventContribution(dPhi(p,q),weight,0.1),id);
0386       deltaR.count(Statistics::EventContribution(dR(p,q),weight,0.1),id);
0387       yDotY.count(Statistics::EventContribution(yy(p,q),weight,0.1),id);
0388     }
0389
0390     /**
0391      * Convert to XML
0392      */
0393     void finalize(XML::Element& elem) {
0394       ObjectProperties::finalize(elem);
0395       deltaY.finalize(); elem.append(deltaY.toXML());
0396       deltaPhi.finalize(); elem.append(deltaPhi.toXML());
0397       deltaR.finalize(); elem.append(deltaR.toXML());
0398       yDotY.finalize(); elem.append(yDotY.toXML());
0399     }
0400
0401   };
0402
0403   /**
0404    * Collection of triple histograms; ranges are adjusted to the
0405    * maximum, so range constraints and rebinning can be applied later.
0406    */
0407   struct TripleProperties
0408     : public ObjectProperties {
0409
0410     /**
0411      * Calculate deltaY^*
0412      */
0413     static double dYstar(const LorentzMomentum& a,
0414                  const LorentzMomentum& b,
0415                          const LorentzMomentum& c){
0416       return c.rapidity()-(a.rapidity()+b.rapidity())/2.;
0417     }
0418
0419     /**
0420      * Calculate deltaZ^* -- normalized deltaY^*
0421      */
0422     static double dZstar(const LorentzMomentum& a,
0423                  const LorentzMomentum& b,
0424                          const LorentzMomentum& c){
0425       return dYstar(a,b,c)*2./abs(a.rapidity()-b.rapidity());
0426     }
0427
0428     /**
0429      * Delta y^*
0430      */
0431     Statistics::Histogram deltaYstar;
0432
0433     /**
0434      * Delta z^*
0435      */
0436     Statistics::Histogram deltaZstar;
0437
0438     /**
0439      * Default constructor
0440      */
0441     TripleProperties()
0442       : ObjectProperties() {}
0443
0444     /**
0445      * Construct given Ecm
0446      */
0447     TripleProperties(const string& name, Energy ecm)
0448       : ObjectProperties(name,ecm),
0449     deltaYstar(name + "DeltaYstar",Statistics::Histogram::regularBinEdges(-6,6,120),true,false),
0450     deltaZstar(name + "DeltaZstar",Statistics::Histogram::regularBinEdges(-3,3,120),true,false) {}
0451
0452     /**
0453      * Count given momentum, weight and id
0454      */
0455     void count(const LorentzMomentum& p, const LorentzMomentum& q, const LorentzMomentum& r,
0456                double weight, unsigned int id) {
0457       ObjectProperties::count(p+q+r,weight,id);
0458       deltaYstar.count(Statistics::EventContribution(dYstar(p,q,r),weight,0.1),id);
0459       deltaZstar.count(Statistics::EventContribution(dZstar(p,q,r),weight,0.05),id);
0460     }
0461
0462     /**
0463      * Convert to XML
0464      */
0465     void finalize(XML::Element& elem) {
0466       ObjectProperties::finalize(elem);
0467       deltaYstar.finalize(); elem.append(deltaYstar.toXML());
0468       deltaZstar.finalize(); elem.append(deltaZstar.toXML());
0469     }
0470
0471   };
0472
0473 private:
0474
0475   /**
0476    * Switch between fixed order and showered
0477    */
0478   bool theIsShowered;
0479
0480   /**
0481    * Switch whether to apply extra analysis cuts
0482    */
0483   bool theApplyCuts;
0484
0485   /**
0486    * The jet finder to use
0487    */
0488   Ptr<JetFinder>::ptr theJetFinder;
0489
0490   /**
0491    * The jet regions to match.
0492    */
0493   vector<Ptr<JetRegion>::ptr> theJetRegions;
0494
0495   /**
0496    * The reconstructed jets
0497    */
0498   map<unsigned int,LorentzMomentum> theJets;
0499
0500   /**
0501    * The reconstructed electroweak particles
0502    */
0503   map<unsigned int,LorentzMomentum> theEWIDs;
0504
0505   /**
0506    * The reconstructed charged leptons
0507    */
0508   map<unsigned int,LorentzMomentum> theChargedLeptons;
0509
0510   /**
0511    * The reconstructed neutrinos
0512    */
0513   map<unsigned int,LorentzMomentum> theNeutrinos;
0514
0515   /**
0516    * The reconstructed missing pT
0517    */
0518   LorentzMomentum thePTmiss;
0519
0520   /**
0521    * The reconstructed Higgs
0522    */
0523   map<unsigned int,LorentzMomentum> theHiggs;
0524
0525   /**
0526    * Jet properties
0527    */
0528   map<unsigned int,ObjectProperties> theJetProperties;
0529
0530   /**
0531    * Exclusive jet properties
0532    */
0533   map<unsigned int,ObjectProperties> theExclusiveJetProperties;
0534
0535   /**
0536    * Jet-inclusive properties
0537    */
0538   ObjectProperties theJetInclusiveProperties;
0539
0540   /**
0541    * Jet-summed properties
0542    */
0543   ObjectProperties theJetSummedProperties;
0544
0545   /**
0546    * Jet-average properties
0547    */
0548   ObjectProperties theJetAverageProperties;
0549
0550   /**
0551    * Inclusive jet multiplicities
0552    */
0553   Statistics::Histogram theNJetsInclusive;
0554
0555   /**
0556    * Exclusive jet multiplicities
0557    */
0558   Statistics::Histogram theNJetsExclusive;
0559
0560   /**
0561    * Electroweak properties
0562    */
0563   map<unsigned int,ObjectProperties> theEWIDProperties;
0564
0565   /**
0566    * Charged lepton properties
0567    */
0568   map<unsigned int,ObjectProperties> theChargedLeptonProperties;
0569
0570   /**
0571    * Neutrino properties
0572    */
0573   map<unsigned int,ObjectProperties> theNeutrinoProperties;
0574
0575   /**
0576    * missing pT properties
0577    */
0578   ObjectProperties thePTmissProperties;
0579
0580   /**
0581    * Higgs properties
0582    */
0583   map<unsigned int,ObjectProperties> theHiggsProperties;
0584
0585   /**
0586    * Jet pair properties
0587    */
0588   map<pair<unsigned int,unsigned int>,PairProperties> theJetPairProperties;
0589
0590   /**
0591    * Jet/electroweak pair properties
0592    */
0593   map<pair<unsigned int,unsigned int>,PairProperties> theJetEWIDPairProperties;
0594
0595   /**
0596    * Jet/charged lepton pair properties
0597    */
0598   map<pair<unsigned int,unsigned int>,PairProperties> theJetChargedLeptonPairProperties;
0599
0600   /**
0601    * Jet/neutrino pair properties
0602    */
0603   map<pair<unsigned int,unsigned int>,PairProperties> theJetNeutrinoPairProperties;
0604
0605   /**
0606    * Jet/missing pT pair properties
0607    */
0608   map<unsigned int,PairProperties> theJetPTmissPairProperties;
0609
0610   /**
0611    * Jet/Higgs pair properties
0612    */
0613   map<pair<unsigned int,unsigned int>,PairProperties> theJetHiggsPairProperties;
0614
0615   /**
0616    * Electroweak pair properties
0617    */
0618   map<pair<unsigned int,unsigned int>,PairProperties> theEWIDPairProperties;
0619
0620   /**
0621    * Charged lepton pair properties
0622    */
0623   map<pair<unsigned int,unsigned int>,PairProperties> theChargedLeptonPairProperties;
0624
0625   /**
0626    * Trijet properties
0627    */
0628   map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeJetProperties;
0629
0630   /**
0631    * Jet-pair/electroweak triple properties
0632    */
0633   map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairEWIDTripleProperties;
0634
0635   /**
0636    * Jet-pair/charged lepton triple properties
0637    */
0638   map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairChargedLeptonTripleProperties;
0639
0640   /**
0641    * Jet-pair/neutrino triple properties
0642    */
0643   map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairNeutrinoTripleProperties;
0644
0645   /**
0646    * Jet-pair/missing pT triple properties
0647    */
0648   map<pair<unsigned int,unsigned int>,TripleProperties> theJetPairPTmissTripleProperties;
0649
0650   /**
0651    * Jet-pair/Higgs triple properties
0652    */
0653   map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theJetPairHiggsTripleProperties;
0654
0655   /**
0656    * Triple electroweak properties
0657    */
0658   map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeEWIDProperties;
0659
0660   /**
0661    * Triple charged lepton properties
0662    */
0663   map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties> theThreeChargedLeptonProperties;
0664
0665   /**
0666    * Fourjet properties
0667    */
0668   map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourJetProperties;
0669
0670   /**
0671    * Four electroweak properties
0672    */
0673   map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourEWIDProperties;
0674
0675   /**
0676    * Four charged lepton properties
0677    */
0678   map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourChargedLeptonProperties;
0679
0680 protected:
0681
0682   /**
0683    * Jet properties
0684    */
0685   ObjectProperties& jetProperties(const unsigned int id) {
0686     map<unsigned int,ObjectProperties>::iterator h =
0687       theJetProperties.find(id);
0688     if ( h != theJetProperties.end() )
0689       return h->second;
0690     ostringstream ids; ids << "Jet" << id;
0691     return
0692       theJetProperties[id] =
0693       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0694   }
0695
0696   /**
0697    * Exclusive jet properties
0698    */
0699   ObjectProperties& exclusiveJetProperties(const unsigned int id) {
0700     map<unsigned int,ObjectProperties>::iterator h =
0701       theExclusiveJetProperties.find(id);
0702     if ( h != theExclusiveJetProperties.end() )
0703       return h->second;
0704     ostringstream ids; ids << "ExclusiveJet" << id;
0705     return
0706       theExclusiveJetProperties[id] =
0707       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0708   }
0709
0710   /**
0711    * Jet-inclusive properties
0712    */
0713   ObjectProperties& jetInclusiveProperties() {
0714     if ( !theJetInclusiveProperties.pt.bins().empty() )
0715       return theJetInclusiveProperties;
0716     return
0717       theJetInclusiveProperties =
0718       ObjectProperties("JetInclusive",generator()->maximumCMEnergy());
0719   }
0720
0721   /**
0722    * Jet-summed properties
0723    */
0724   ObjectProperties& jetSummedProperties() {
0725     if ( !theJetSummedProperties.pt.bins().empty() )
0726       return theJetSummedProperties;
0727     return
0728       theJetSummedProperties =
0729       ObjectProperties("JetSummed",generator()->maximumCMEnergy());
0730   }
0731
0732   /**
0733    * Jet-average properties
0734    */
0735   ObjectProperties& jetAverageProperties() {
0736     if ( !theJetAverageProperties.pt.bins().empty() )
0737       return theJetAverageProperties;
0738     return
0739       theJetAverageProperties =
0740       ObjectProperties("JetAverage",generator()->maximumCMEnergy());
0741   }
0742
0743
0744   /**
0745    * Inclusive jet multiplicities
0746    */
0747   Statistics::Histogram& nJetsInclusive() {
0748     if ( !theNJetsInclusive.bins().empty() )
0749       return theNJetsInclusive;
0750     return
0751       theNJetsInclusive =
0752       Statistics::Histogram("NJetsInclusive",
0753                 Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
0754                                    theJetRegions.size()+1),
0755                 true,true);
0756   }
0757
0758   /**
0759    * Exclusive jet multiplicities
0760    */
0761   Statistics::Histogram& nJetsExclusive() {
0762     if ( !theNJetsExclusive.bins().empty() )
0763       return theNJetsExclusive;
0764     return
0765       theNJetsExclusive =
0766       Statistics::Histogram("NJetsExclusive",
0767                 Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
0768                                    theJetRegions.size()+1),
0769                 true,true);
0770   }
0771
0772   /**
0773    * Lepton properties -- all sorted by ID
0774    */
0775   ObjectProperties& eWIDProperties(const unsigned int id) {
0776     map<unsigned int,ObjectProperties>::iterator h =
0777       theEWIDProperties.find(id);
0778     if ( h != theEWIDProperties.end() )
0779       return h->second;
0780     ostringstream ids; ids << "EWID" << id;
0781     return
0782       theEWIDProperties[id] =
0783       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0784   }
0785
0786   /**
0787    * Charged lepton properties
0788    */
0789   ObjectProperties& chargedLeptonProperties(const unsigned int id) {
0790     map<unsigned int,ObjectProperties>::iterator h =
0791       theChargedLeptonProperties.find(id);
0792     if ( h != theChargedLeptonProperties.end() )
0793       return h->second;
0794     ostringstream ids; ids << "ChargedLepton" << id;
0795     return
0796       theChargedLeptonProperties[id] =
0797       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0798   }
0799
0800   /**
0801    * Neutrino properties
0802    */
0803   ObjectProperties& neutrinoProperties(const unsigned int id) {
0804     map<unsigned int,ObjectProperties>::iterator h =
0805       theNeutrinoProperties.find(id);
0806     if ( h != theNeutrinoProperties.end() )
0807       return h->second;
0808     ostringstream ids; ids << "Neutrino" << id;
0809     return
0810       theNeutrinoProperties[id] =
0811       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0812   }
0813
0814   /**
0815    * Missing pT properties
0816    */
0817   ObjectProperties& pTmissProperties() {
0818     if ( !thePTmissProperties.pt.bins().empty() )
0819       return thePTmissProperties;
0820     return
0821       thePTmissProperties =
0822       ObjectProperties("PTmiss",generator()->maximumCMEnergy());
0823   }
0824
0825   /**
0826    * Higgs properties
0827    */
0828   ObjectProperties& higgsProperties(const unsigned int id) {
0829     map<unsigned int,ObjectProperties>::iterator h =
0830       theHiggsProperties.find(id);
0831     if ( h != theHiggsProperties.end() )
0832       return h->second;
0833     ostringstream ids; ids << "Higgs" << id;
0834     return
0835       theHiggsProperties[id] =
0836       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0837   }
0838
0839   /**
0840    * Jet pair properties
0841    */
0842   PairProperties& jetPairProperties(const unsigned int id, const unsigned int jd) {
0843     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0844       theJetPairProperties.find(make_pair(id,jd));
0845     if ( h != theJetPairProperties.end() )
0846       return h->second;
0847     ostringstream ids; ids << "Jet" << id << jd;
0848     return theJetPairProperties[make_pair(id,jd)] =
0849       PairProperties(ids.str(),generator()->maximumCMEnergy());
0850   }
0851
0852   /**
0853    * Jet/lepton(all sorted by ID) pair properties
0854    */
0855   PairProperties& jetEWIDPairProperties(const unsigned int id, const unsigned int jd) {
0856     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0857       theJetEWIDPairProperties.find(make_pair(id,jd));
0858     if ( h != theJetEWIDPairProperties.end() )
0859       return h->second;
0860     ostringstream ids; ids << "Jet" << id << "EWID" << jd;
0861     return theJetEWIDPairProperties[make_pair(id,jd)] =
0862       PairProperties(ids.str(),generator()->maximumCMEnergy());
0863   }
0864
0865   /**
0866    * Jet/charged lepton pair properties
0867    */
0868   PairProperties& jetChargedLeptonPairProperties(const unsigned int id, const unsigned int jd) {
0869     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0870       theJetChargedLeptonPairProperties.find(make_pair(id,jd));
0871     if ( h != theJetChargedLeptonPairProperties.end() )
0872       return h->second;
0873     ostringstream ids; ids << "Jet" << id << "ChargedLepton" << jd;
0874     return theJetChargedLeptonPairProperties[make_pair(id,jd)] =
0875       PairProperties(ids.str(),generator()->maximumCMEnergy());
0876   }
0877
0878   /**
0879    * Jet/neutrino pair properties
0880    */
0881   PairProperties& jetNeutrinoPairProperties(const unsigned int id, const unsigned int jd) {
0882     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0883       theJetNeutrinoPairProperties.find(make_pair(id,jd));
0884     if ( h != theJetNeutrinoPairProperties.end() )
0885       return h->second;
0886     ostringstream ids; ids << "Jet" << id << "Neutrino" << jd;
0887     return theJetNeutrinoPairProperties[make_pair(id,jd)] =
0888       PairProperties(ids.str(),generator()->maximumCMEnergy());
0889   }
0890
0891   /**
0892    * Jet/missing pT pair properties
0893    */
0894   PairProperties& jetPTmissPairProperties(const unsigned int id) {
0895     map<unsigned int,PairProperties>::iterator h =
0896       theJetPTmissPairProperties.find(id);
0897     if ( h != theJetPTmissPairProperties.end() )
0898       return h->second;
0899     ostringstream ids; ids << "Jet" << id << "PTmiss";
0900     return theJetPTmissPairProperties[id] =
0901       PairProperties(ids.str(),generator()->maximumCMEnergy());
0902   }
0903
0904   /**
0905    * Jet/Higgs pair properties
0906    */
0907   PairProperties& jetHiggsPairProperties(const unsigned int id, const unsigned int jd) {
0908     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0909       theJetHiggsPairProperties.find(make_pair(id,jd));
0910     if ( h != theJetHiggsPairProperties.end() )
0911       return h->second;
0912     ostringstream ids; ids << "Jet" << id << "Higgs" << jd;
0913     return theJetHiggsPairProperties[make_pair(id,jd)] =
0914       PairProperties(ids.str(),generator()->maximumCMEnergy());
0915   }
0916
0917   /**
0918    * Electroweak pair properties
0919    */
0920   PairProperties& eWIDPairProperties(const unsigned int id, const unsigned int jd) {
0921     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0922       theEWIDPairProperties.find(make_pair(id,jd));
0923     if ( h != theEWIDPairProperties.end() )
0924       return h->second;
0925     ostringstream ids; ids << "EWID" << id << jd;
0926     return theEWIDPairProperties[make_pair(id,jd)] =
0927       PairProperties(ids.str(),generator()->maximumCMEnergy());
0928   }
0929
0930   /**
0931    * Charged lepton pair properties
0932    */
0933   PairProperties& chargedLeptonPairProperties(const unsigned int id, const unsigned int jd) {
0934     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0935       theChargedLeptonPairProperties.find(make_pair(id,jd));
0936     if ( h != theChargedLeptonPairProperties.end() )
0937       return h->second;
0938     ostringstream ids; ids << "ChargedLepton" << id << jd;
0939     return theChargedLeptonPairProperties[make_pair(id,jd)] =
0940       PairProperties(ids.str(),generator()->maximumCMEnergy());
0941   }
0942
0943   /**
0944    * Trijet properties
0945    */
0946   TripleProperties& threeJetProperties(const unsigned int id1, const unsigned int id2,
0947                        const unsigned int id3) {
0948     map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
0949       theThreeJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
0950     if ( it != theThreeJetProperties.end() )
0951       return it->second;
0952     ostringstream ids;
0953     ids << "Jet" << id1 << id2 << id3;
0954     return theThreeJetProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
0955       TripleProperties(ids.str(),generator()->maximumCMEnergy());
0956   }
0957
0958   /**
0959    * Jet-pair/electroweak triple properties
0960    */
0961   TripleProperties& jetPairEWIDTripleProperties(const unsigned int id1, const unsigned int id2,
0962                               const unsigned int id3) {
0963     map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
0964       theJetPairEWIDTripleProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
0965     if ( it != theJetPairEWIDTripleProperties.end() )
0966       return it->second;
0967     ostringstream ids;
0968     ids << "Jet" << id1 << id2 << "EWID" << id3;
0969     return theJetPairEWIDTripleProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
0970       TripleProperties(ids.str(),generator()->maximumCMEnergy());
0971   }
0972
0973   /**
0974    * Jet-pair/charged lepton triple properties
0975    */
0976   TripleProperties& jetPairChargedLeptonTripleProperties(const unsigned int id1, const unsigned int id2,
0977                               const unsigned int id3) {
0978     map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
0979       theJetPairChargedLeptonTripleProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
0980     if ( it != theJetPairChargedLeptonTripleProperties.end() )
0981       return it->second;
0982     ostringstream ids;
0983     ids << "Jet" << id1 << id2 << "ChargedLepton" << id3;
0984     return theJetPairChargedLeptonTripleProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
0985       TripleProperties(ids.str(),generator()->maximumCMEnergy());
0986   }
0987
0988   /**
0989    * Jet-pair/neutrino triple properties
0990    */
0991   TripleProperties& jetPairNeutrinoTripleProperties(const unsigned int id1, const unsigned int id2,
0992                               const unsigned int id3) {
0993     map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
0994       theJetPairNeutrinoTripleProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
0995     if ( it != theJetPairNeutrinoTripleProperties.end() )
0996       return it->second;
0997     ostringstream ids;
0998     ids << "Jet" << id1 << id2 << "Neutrino" << id3;
0999     return theJetPairNeutrinoTripleProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
1000       TripleProperties(ids.str(),generator()->maximumCMEnergy());
1001   }
1002
1003   /**
1004    * Jet-pair/missing pT triple properties
1005    */
1006   TripleProperties& jetPairPTmissTripleProperties(const unsigned int id1, const unsigned int id2) {
1007     map<pair<unsigned int,unsigned int>,TripleProperties>::iterator it =
1008       theJetPairPTmissTripleProperties.find(pair<unsigned int,unsigned int>(id1,id2));
1009     if ( it != theJetPairPTmissTripleProperties.end() )
1010       return it->second;
1011     ostringstream ids;
1012     ids << "Jet" << id1 << id2 << "PTmiss";
1013     return theJetPairPTmissTripleProperties[pair<unsigned int,unsigned int>(id1,id2)] =
1014       TripleProperties(ids.str(),generator()->maximumCMEnergy());
1015   }
1016
1017   /**
1018    * Jet-pair/Higgs triple properties
1019    */
1020   TripleProperties& jetPairHiggsTripleProperties(const unsigned int id1, const unsigned int id2,
1021                               const unsigned int id3) {
1022     map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
1023       theJetPairHiggsTripleProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
1024     if ( it != theJetPairHiggsTripleProperties.end() )
1025       return it->second;
1026     ostringstream ids;
1027     ids << "Jet" << id1 << id2 << "Higgs" << id3;
1028     return theJetPairHiggsTripleProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
1029       TripleProperties(ids.str(),generator()->maximumCMEnergy());
1030   }
1031
1032   /**
1033    * Triple electroweak properties -- all sorted by ID
1034    */
1035   TripleProperties& threeEWIDProperties(const unsigned int id1, const unsigned int id2,
1036                         const unsigned int id3) {
1037     map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
1038       theThreeEWIDProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
1039     if ( it != theThreeEWIDProperties.end() )
1040       return it->second;
1041     ostringstream ids;
1042     ids << "EWID" << id1 << id2 << id3;
1043     return theThreeEWIDProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
1044       TripleProperties(ids.str(),generator()->maximumCMEnergy());
1045   }
1046
1047   /**
1048    * Triple charged lepton properties
1049    */
1050   TripleProperties& threeChargedLeptonProperties(const unsigned int id1, const unsigned int id2,
1051                        const unsigned int id3) {
1052     map<std::tuple<unsigned int,unsigned int,unsigned int>,TripleProperties>::iterator it =
1053       theThreeChargedLeptonProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
1054     if ( it != theThreeChargedLeptonProperties.end() )
1055       return it->second;
1056     ostringstream ids;
1057     ids << "ChargedLepton" << id1 << id2 << id3;
1058     return theThreeChargedLeptonProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
1059       TripleProperties(ids.str(),generator()->maximumCMEnergy());
1060   }
1061
1062   /**
1063    * Fourjet properties
1064    */
1065   ObjectProperties& fourJetProperties(const unsigned int id1, const unsigned int id2,
1066                       const unsigned int id3, const unsigned int id4) {
1067     map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
1068       theFourJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
1069     if ( it != theFourJetProperties.end() )
1070       return it->second;
1071     ostringstream ids;
1072     ids << "Jet" << id1 << id2 << id3 << id4;
1073     return theFourJetProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
1074       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
1075   }
1076
1077   /**
1078    * Four electroweak properties
1079    */
1080   ObjectProperties& fourEWIDProperties(const unsigned int id1, const unsigned int id2,
1081                       const unsigned int id3, const unsigned int id4) {
1082     map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
1083       theFourEWIDProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
1084     if ( it != theFourEWIDProperties.end() )
1085       return it->second;
1086     ostringstream ids;
1087     ids << "EWID" << id1 << id2 << id3 << id4;
1088     return theFourEWIDProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
1089       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
1090   }
1091
1092   /**
1093    * Four charged lepton properties
1094    */
1095   ObjectProperties& fourChargedLeptonProperties(const unsigned int id1, const unsigned int id2,
1096                       const unsigned int id3, const unsigned int id4) {
1097     map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
1098       theFourChargedLeptonProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
1099     if ( it != theFourChargedLeptonProperties.end() )
1100       return it->second;
1101     ostringstream ids;
1102     ids << "ChargedLepton" << id1 << id2 << id3 << id4;
1103     return theFourChargedLeptonProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
1104       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
1105   }
1106
1107   /**
1108    * Perform any additional analysis required
1109    */
1110   virtual void analyzeSpecial(long, double) {}
1111
1112   /**
1113    * Append any additional histograms to the given histogram element
1114    */
1115   virtual void finalize(XML::Element&) {}
1116
1117 private:
1118
1119   /**
1120    * The assignment operator is private and must never be called.
1121    * In fact, it should not even be implemented.
1122    */
1123   LeptonsJetsAnalysis & operator=(const LeptonsJetsAnalysis &) = delete;
1124
1125 };
1126
1127 }
1128
1129 #endif /* Herwig_LeptonsJetsAnalysis_H */