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/DressedLeptons.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      IdentifiedFinalState bareElectrons(fs);
 43      bareElectrons.acceptIdPair(PID::ELECTRON);
 44
 45      Cut cuts = (Cuts::absetaIn(0, 1.37) || Cuts::absetaIn(1.52, 2.47)) && Cuts::pT > 20*GeV;
 46
 47      DressedLeptons electronClusters(fs, bareElectrons, 0.1, cuts);
 48      declare(electronClusters, "electronClusters");
 49
 50      IdentifiedFinalState bareMuons(fs);
 51      bareMuons.acceptIdPair(PID::MUON);
 52      Cut mucuts = Cuts::abseta < 2.4 && Cuts::pT > 20*GeV;
 53      DressedLeptons muonClusters(fs, bareMuons, 0.1, mucuts);
 54      declare(muonClusters, "muonClusters");
 55
 56      IdentifiedFinalState neutrinos(Cuts::pT > 25*GeV);
 57      neutrinos.acceptNeutrinos();
 58      declare(neutrinos, "neutrinos");
 59
 60      VetoedFinalState jetFS(fs);
 61      jetFS.addVetoOnThisFinalState(electronClusters);
 62      jetFS.addVetoOnThisFinalState(muonClusters);
 63      jetFS.addVetoOnThisFinalState(neutrinos);
 64      FastJets jetpro(jetFS, FastJets::KT, 0.6, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
 65      declare(jetpro, "jets");
 66
 67      // Book histograms
 68      for (size_t flav = 0; flav < 2; ++flav) {
 69        for (size_t i = 0; i < m_njet; ++i)   book(_h_dI[flav][i],         i+1, 1, flav+1);
 70        for (size_t i = 0; i < m_njet-1; ++i) book(_h_dI_ratio[flav][i], 4+i+1, 1, flav+1);
 71      }
 72    }
 73
 74
 75    /// Perform the per-event analysis
 76    void analyze(const Event& e) {
 77      const DressedLeptons& electronClusters = apply<DressedLeptons>(e, "electronClusters");
 78      const DressedLeptons& muonClusters = apply<DressedLeptons>(e, "muonClusters");
 79      int ne = electronClusters.dressedLeptons().size();
 80      int nmu = muonClusters.dressedLeptons().size();
 81
 82      FourMomentum lepton;
 83      size_t flav = 2;
 84      if (ne==1) {
 85        lepton=electronClusters.dressedLeptons()[0].momentum();
 86        flav = 0;
 87        if (nmu > 0) vetoEvent;
 88      }
 89      else if (nmu == 1) {
 90        lepton=muonClusters.dressedLeptons()[0].momentum();
 91        flav = 1;
 92        if (ne > 0) vetoEvent;
 93      }
 94      else {
 95        vetoEvent;
 96      }
 97
 98      const Particles& neutrinos = apply<FinalState>(e, "neutrinos").particlesByPt();
 99      if (neutrinos.size() < 1) vetoEvent;
100      FourMomentum neutrino = neutrinos[0].momentum();
101
102      double mtW=sqrt(2.0*lepton.pT()*neutrino.pT()*(1-cos(lepton.phi()-neutrino.phi())));
103      if (mtW<40.0*GeV) vetoEvent;
104
105      const shared_ptr<fastjet::ClusterSequence> seq = apply<FastJets>(e, "jets").clusterSeq();
106      if (seq) {
107        for (size_t i = 0; i < min(m_njet,(size_t)seq->n_particles()); ++i) {
108          double d_ij = sqrt(seq->exclusive_dmerge_max(i));
109          _h_dI[flav][i]->fill(d_ij);
110
111          if (i<m_njet-1) {
112            if (d_ij>20.0*GeV) {
113              double d_ijplus1 = sqrt(seq->exclusive_dmerge_max(i+1));
114              _h_dI_ratio[flav][i]->fill(d_ijplus1/d_ij);
115            }
116          }
117        }
118      }
119
120    }
121
122
123    /// Normalise histograms etc., after the run
124    void finalize() {
125      for (size_t flav = 0; flav < 2; ++flav) {
126        for (size_t i = 0; i < m_njet; ++i) {
127          normalize(_h_dI[flav][i], 1.0, false);
128          if (i < m_njet-1) normalize(_h_dI_ratio[flav][i], 1.0, false);
129        }
130      }
131    }
132
133    //@}
134
135
136  private:
137
138    /// @name Histograms
139    //@{
140    vector< vector<Histo1DPtr> > _h_dI;
141    vector< vector<Histo1DPtr> > _h_dI_ratio;
142    //@}
143
144    size_t m_njet;
145  };
146
147
148  // The hook for the plugin system
149  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1217867);
150
151
152}