JetsPlusAnalysis.h

// -*- C++ -*-
#ifndef Herwig_JetsPlusAnalysis_H
#define Herwig_JetsPlusAnalysis_H
//
// This is the declaration of the JetsPlusAnalysis 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 JetsPlusAnalysis: public AnalysisHandler {
 
public:
 
  JetsPlusAnalysis();
 
public:
 
  virtual void analyze(tEventPtr event, long ieve, int loop, int state);
 
protected:
 
  void analyze(ParticleVector&, long, double);
 
  void clear() {
    theHardObjects.clear();
    theJets.clear();
  }
 
  virtual void reconstructHardObjects(ParticleVector&) {}
 
  LorentzMomentum& hardObjectMomentum(const string& id) {
    return theHardObjects[id];
  }
 
  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];
  }
 
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 pt_logx;
 
    Statistics::Histogram y;
 
    Statistics::Histogram phi;
 
    Statistics::Histogram mass;
 
    ObjectProperties() {}
 
    ObjectProperties(const string& name, Energy)
      : pt(name + "Pt",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false),
        pt_logx(name + "PtLogX",Statistics::Histogram::logBinEdges(0.1,1000,1000),true,false),
        y(name + "Y",Statistics::Histogram::regularBinEdges(-6,6,120),false,false),
        phi(name + "Phi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
            make_pair(-Constants::pi,Constants::pi)),
        mass(name + "Mass",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false) {}
 
    void count(const LorentzMomentum& p, double weight, unsigned int id) {
      pt.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),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,1.),id);
    }
 
    void count(Energy perp, double rapidity, 
               double xphi, Energy m,
               double weight, unsigned int id) {
      pt.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,1.),id);
    }
 
    void finalize(XML::Element& elem) {
      pt.finalize(); elem.append(pt.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());
    }
 
  };
 
  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());
    }
 
  };
 
private:
 
  bool theIsShowered;
 
  Ptr<JetFinder>::ptr theJetFinder;
 
  vector<Ptr<JetRegion>::ptr> theJetRegions;
 
  map<string,LorentzMomentum> theHardObjects;
 
  map<unsigned int,LorentzMomentum> theJets;
 
  map<string,ObjectProperties> theHardObjectProperties;
 
  map<unsigned int,ObjectProperties> theJetProperties;
 
  map<unsigned int,ObjectProperties> theExclusiveJetProperties;
 
  ObjectProperties theJetInclusiveProperties;
 
  ObjectProperties theJetSummedProperties;
 
  ObjectProperties theJetAverageProperties;
 
  Statistics::Histogram theNJetsInclusive;
 
  Statistics::Histogram theNJetsExclusive;
 
  map<pair<string,string>,PairProperties> theHardPairProperties;
 
  map<pair<unsigned int,unsigned int>,PairProperties> theJetPairProperties;
 
  map<pair<unsigned int,string>,PairProperties> theJetHardPairProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties> theThreeJetProperties;
 
  map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourJetProperties;
 
protected:
 
  ObjectProperties& hardObjectProperties(const string& id) {
    map<string,ObjectProperties>::iterator h = 
      theHardObjectProperties.find(id);
    if ( h != theHardObjectProperties.end() )
      return h->second;
    return 
      theHardObjectProperties[id] =
      ObjectProperties(id,generator()->maximumCMEnergy());
  }
 
  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);
  }
 
  PairProperties& hardPairProperties(const string& id, const string& jd) {
    map<pair<string,string>,PairProperties>::iterator h = 
      theHardPairProperties.find(make_pair(id,jd));
    if ( h != theHardPairProperties.end() )
      return h->second;
    return theHardPairProperties[make_pair(id,jd)] = 
      PairProperties(id+jd,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& jetHardPairProperties(const unsigned int id, const string& jd) {
    map<pair<unsigned int,string>,PairProperties>::iterator h = 
      theJetHardPairProperties.find(make_pair(id,jd));
    if ( h != theJetHardPairProperties.end() )
      return h->second;
    ostringstream ids; ids << "Jet" << id << jd;
    return theJetHardPairProperties[make_pair(id,jd)] = 
      PairProperties(ids.str(),generator()->maximumCMEnergy());
  }
 
  ObjectProperties& threeJetProperties(const unsigned int id1, const unsigned int id2,
                                       const unsigned int id3) {
    map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties>::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)] =
      ObjectProperties(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());
  }  
 
  virtual void analyzeSpecial(long, double) {}
 
  virtual void finalize(XML::Element&) {}
 
private:
 
  JetsPlusAnalysis & operator=(const JetsPlusAnalysis &) = delete;
 
};
 
}
 
#endif /* Herwig_JetsPlusAnalysis_H */