1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12
 13  class ATLAS_2013_I1217867 : public Analysis {
 14  public:
 15
 16    /// @name Constructors etc.
 17    /// @{
 18
 19    /// Constructor
 20    ATLAS_2013_I1217867()
 21      : Analysis("ATLAS_2013_I1217867")
 22    {
 23      m_njet = 4;
 24      _h_dI.resize(2, std::vector<Histo1DPtr>(m_njet));
 25      _h_dI_ratio.resize(2, std::vector<Histo1DPtr>(m_njet-1));
 26    }
 27
 28    /// @}
 29
 30
 31  public:
 32
 33    /// @name Analysis methods
 34    /// @{
 35
 36    /// Book histograms and initialise projections before the run
 37    void init() {
 38      // Initialise projections
 39
 40      FinalState fs(Cuts::abseta < 5.0);
 41
 42      Cut cuts = (Cuts::absetaIn(0, 1.37) || Cuts::absetaIn(1.52, 2.47)) && Cuts::pT > 20*GeV;
 43      LeptonFinder electronClusters(0.1, Cuts::abspid == PID::ELECTRON && cuts);
 44      declare(electronClusters, "electronClusters");
 45
 46      Cut mucuts = Cuts::abseta < 2.4 && Cuts::pT > 20*GeV;
 47      LeptonFinder muonClusters(0.1, Cuts::abspid == PID::MUON && mucuts);
 48      declare(muonClusters, "muonClusters");
 49
 50      IdentifiedFinalState neutrinos(Cuts::pT > 25*GeV);
 51      neutrinos.acceptNeutrinos();
 52      declare(neutrinos, "neutrinos");
 53
 54      VetoedFinalState jetFS(fs);
 55      jetFS.addVetoOnThisFinalState(electronClusters);
 56      jetFS.addVetoOnThisFinalState(muonClusters);
 57      jetFS.addVetoOnThisFinalState(neutrinos);
 58      FastJets jetpro(jetFS, JetAlg::KT, 0.6, JetMuons::ALL, JetInvisibles::DECAY);
 59      declare(jetpro, "jets");
 60
 61      // Book histograms
 62      for (size_t flav = 0; flav < 2; ++flav) {
 63        for (size_t i = 0; i < m_njet; ++i)   book(_h_dI[flav][i],         i+1, 1, flav+1);
 64        for (size_t i = 0; i < m_njet-1; ++i) book(_h_dI_ratio[flav][i], 4+i+1, 1, flav+1);
 65      }
 66    }
 67
 68
 69    /// Perform the per-event analysis
 70    void analyze(const Event& e) {
 71      const LeptonFinder& electronClusters = apply<LeptonFinder>(e, "electronClusters");
 72      const LeptonFinder& muonClusters = apply<LeptonFinder>(e, "muonClusters");
 73      int ne = electronClusters.dressedLeptons().size();
 74      int nmu = muonClusters.dressedLeptons().size();
 75
 76      FourMomentum lepton;
 77      size_t flav = 2;
 78      if (ne==1) {
 79        lepton=electronClusters.dressedLeptons()[0].momentum();
 80        flav = 0;
 81        if (nmu > 0) vetoEvent;
 82      }
 83      else if (nmu == 1) {
 84        lepton=muonClusters.dressedLeptons()[0].momentum();
 85        flav = 1;
 86        if (ne > 0) vetoEvent;
 87      }
 88      else {
 89        vetoEvent;
 90      }
 91
 92      const Particles& neutrinos = apply<FinalState>(e, "neutrinos").particlesByPt();
 93      if (neutrinos.size() < 1) vetoEvent;
 94      FourMomentum neutrino = neutrinos[0].momentum();
 95
 96      double mtW=sqrt(2.0*lepton.pT()*neutrino.pT()*(1-cos(lepton.phi()-neutrino.phi())));
 97      if (mtW<40.0*GeV) vetoEvent;
 98
 99      const shared_ptr<fastjet::ClusterSequence> seq = apply<FastJets>(e, "jets").clusterSeq();
100      if (seq) {
101        for (size_t i = 0; i < min(m_njet,(size_t)seq->n_particles()); ++i) {
102          double d_ij = sqrt(seq->exclusive_dmerge_max(i));
103          _h_dI[flav][i]->fill(d_ij);
104
105          if (i<m_njet-1) {
106            if (d_ij>20.0*GeV) {
107              double d_ijplus1 = sqrt(seq->exclusive_dmerge_max(i+1));
108              _h_dI_ratio[flav][i]->fill(d_ijplus1/d_ij);
109            }
110          }
111        }
112      }
113
114    }
115
116
117    /// Normalise histograms etc., after the run
118    void finalize() {
119      for (size_t flav = 0; flav < 2; ++flav) {
120        for (size_t i = 0; i < m_njet; ++i) {
121          normalize(_h_dI[flav][i], 1.0, false);
122          if (i < m_njet-1) normalize(_h_dI_ratio[flav][i], 1.0, false);
123        }
124      }
125    }
126
127    /// @}
128
129
130  private:
131
132    /// @name Histograms
133    /// @{
134    vector< vector<Histo1DPtr> > _h_dI;
135    vector< vector<Histo1DPtr> > _h_dI_ratio;
136    /// @}
137
138    size_t m_njet;
139  };
140
141
142  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1217867);
143
144
145}