1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/DileptonFinder.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Projections/VetoedFinalState.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// Measurement of V+jets distributions, taken as ratios between W and Z channels
 15  class ATLAS_2014_I1312627 : public Analysis {
 16  public:
 17
 18    /// @name Plotting helper types
 19    /// @{
 20
 21    struct Plots {
 22      string ref;
 23      Histo1DPtr comp[2]; // (de)nominator components
 24      Estimate1DPtr ratio; // Rjets plot
 25    };
 26
 27    /// @}
 28
 29
 30    /// Constructor
 31    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2014_I1312627);
 32
 33
 34    /// @name Analysis methods
 35    /// @{
 36
 37    /// Book histograms and initialise projections before the run
 38    void init() {
 39
 40      // get option
 41      _mode = 0;
 42      if ( getOption("LMODE") == "EL" )  _mode = 1;
 43      if ( getOption("LMODE") == "MU" )  _mode = 2;
 44
 45      // Set up cuts
 46      Cut cuts;
 47      if (_mode == 2) { // muon channel
 48        cuts = Cuts::pT > 25*GeV && Cuts::abseta < 2.4;;
 49      } else if (_mode) { // electron channel
 50        cuts = Cuts::pT > 25*GeV && ( Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47) );
 51      } else { // combined data extrapolated to common phase space
 52        cuts = Cuts::pT > 25*GeV && Cuts::abseta < 2.5;
 53      }
 54
 55      // Boson finders
 56      declare("MET", MissingMomentum());
 57      LeptonFinder lf(0.1, cuts && Cuts::abspid == (_mode > 1 ? PID::MUON : PID::ELECTRON));
 58      declare(lf, "Leptons");
 59      DileptonFinder zf(91.2*GeV, 0.1, cuts && Cuts::abspid == (_mode > 1 ? PID::MUON : PID::ELECTRON), Cuts::massIn(66*GeV, 116*GeV));
 60      declare(zf, "ZF");
 61
 62      // Jets
 63      VetoedFinalState jet_fs;
 64      jet_fs.addVetoOnThisFinalState(lf);
 65      jet_fs.addVetoOnThisFinalState(zf);
 66      FastJets jets(jet_fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL); //< !!
 67      declare(jets, "Jets");
 68
 69
 70      // Book Rjets plots
 71      _suff = string("-y0") + to_str(_mode + 1);
 72      hInit("Njets_incl",  "d01"); // inclusive number of jets
 73      hInit("Njets_excl",  "d04"); // exclusive number of jets
 74      hInit("Pt1_N1incl",  "d05"); // leading jet pT, at least 1 jet
 75      hInit("Pt1_N1excl",  "d06"); // leading jet pT, exactly 1 jet
 76      hInit("Pt1_N2incl",  "d07"); // leading jet pT, at least 2 jets
 77      hInit("Pt1_N3incl",  "d08"); // leading jet pT, at least 3 jets
 78      hInit("Pt2_N2incl",  "d09"); // subleading jet pT, at least 2 jets
 79      hInit("Pt3_N3incl",  "d10"); // sub-subleading jet pT, at least 3 jets
 80      hInit("ST_N2incl",   "d11"); // scalar jet pT sum, at least 2 jets
 81      hInit("ST_N2excl",   "d12"); // scalar jet pT sum, exactly 2 jets
 82      hInit("ST_N3incl",   "d13"); // scalar jet pT sum, at least 3 jets
 83      hInit("ST_N3excl",   "d14"); // scalar jet pT sum, exactly 3 jets
 84      hInit("DR_N2incl",   "d15"); // deltaR(j1, j2), at least 2 jets
 85      hInit("DPhi_N2incl", "d16"); // deltaPhi(j1, j2), at least 2 jets
 86      hInit("Mjj_N2incl",  "d17"); // mjj, at least 2 jets
 87      hInit("Rap1_N1incl", "d18"); // leading jet rapidity, at least 1 jet
 88      hInit("Rap2_N2incl", "d19"); // subleading jet rapidity, at least 2 jets
 89      hInit("Rap3_N3incl", "d20"); // sub-subleading jet rapidity, at least 3 jets
 90
 91      // Also book numerator and denominator for Rjets plots
 92      for (auto& item : _plots) {
 93        book(item.second.comp[0], item.second.ref + "2" + _suff, refData(item.second.ref + "1" + _suff));
 94        book(item.second.comp[1], item.second.ref + "3" + _suff, refData(item.second.ref + "1" + _suff));
 95      }
 96    }
 97
 98
 99    /// Perform the per-event analysis
100    void analyze(const Event& event) {
101
102      // W reco, starting with MET
103      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
104      const Particles& ls = apply<LeptonFinder>(event, "Leptons").particles();
105      const Particles ls_mtfilt = select(ls, [&](const Particle& l){ return mT(l, pmiss) > 40*GeV; });
106      const int ifound = closestMatchIndex(ls_mtfilt, pmiss, Kin::mass, 80.4*GeV);
107
108      // Z reco
109      const DileptonFinder& zf = apply<DileptonFinder>(event, "ZF");
110
111      // Exit if no bosons found (note: overlapping W and Z reco is allowed)
112      if (ifound < 0 && zf.empty()) vetoEvent;
113
114      // Retrieve jets
115      const JetFinder& jetfs = apply<JetFinder>(event, "Jets");
116      Jets jets = jetfs.jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4);
117
118      // Apply boson cuts and fill histograms
119      if (!zf.empty()) {
120        const Particles& leptons = zf.constituents();
121        if (oppSign(leptons[0], leptons[1]) && deltaR(leptons[0], leptons[1]) > 0.2)
122          fillPlots(leptons, jets, 1);
123      }
124      if (ifound >= 0) {
125      const Particle& lepton = ls_mtfilt[ifound];
126      if (pmiss.pT() > 25*GeV)
127          fillPlots(Particles{lepton}, jets, 0);
128      }
129    }
130
131
132    /// Normalise histograms etc., after the run
133    void finalize() {
134      ///  Normalise, scale and otherwise manipulate histograms here
135      const double sf( crossSection()/picobarn / sumOfWeights() );
136      for (auto& item : _plots) {
137        scale(item.second.comp[0], sf);
138        scale(item.second.comp[1], sf);
139        divide(item.second.comp[0], item.second.comp[1], item.second.ratio);
140      }
141    }
142
143    /// @}
144
145
146    /// Analysis helper functions
147    /// @{
148
149    /// Histogram filling function, to avoid duplication
150    void fillPlots(const Particles& leptons, Jets& jets, int isZ) {
151      // Do jet-lepton overlap removal
152      idiscardIfAnyDeltaRLess(jets, leptons, 0.5);
153
154      // Calculate jet ST
155      const size_t njets = jets.size();
156      const double ST = sum(jets, pT, 0.0)/GeV;
157
158      // Fill jet histos
159      _plots["Njets_excl"].comp[isZ]->fill(njets + 0.5);
160      for (size_t i = 0; i <= njets; ++i)
161        _plots["Njets_incl"].comp[isZ]->fill(i + 0.5);
162
163      if (njets > 0) {
164        const double pT1  = jets[0].pT()/GeV;
165        const double rap1 = jets[0].absrap();
166        _plots["Pt1_N1incl" ].comp[isZ]->fill(pT1);
167        _plots["Rap1_N1incl"].comp[isZ]->fill(rap1);
168
169        if (njets == 1) {
170          _plots["Pt1_N1excl"].comp[isZ]->fill(pT1);
171        } else if (njets > 1) {
172          const double pT2  = jets[1].pT()/GeV;
173          const double rap2 = jets[1].absrap();
174          const double dR   = deltaR(jets[0], jets[1]);
175          const double dPhi = deltaPhi(jets[0], jets[1]);
176          const double mjj  = (jets[0].momentum() + jets[1].momentum()).mass()/GeV;
177          _plots["Pt1_N2incl" ].comp[isZ]->fill(pT1);
178          _plots["Pt2_N2incl" ].comp[isZ]->fill(pT2);
179          _plots["Rap2_N2incl"].comp[isZ]->fill(rap2);
180          _plots["DR_N2incl"  ].comp[isZ]->fill(dR);
181          _plots["DPhi_N2incl"].comp[isZ]->fill(dPhi);
182          _plots["Mjj_N2incl" ].comp[isZ]->fill(mjj);
183          _plots["ST_N2incl"  ].comp[isZ]->fill(ST);
184
185          if (njets == 2) {
186            _plots["ST_N2excl"].comp[isZ]->fill(ST);
187          } else if (njets > 2) {
188            const double pT3  = jets[2].pT()/GeV;
189            const double rap3 = jets[2].absrap();
190            _plots["Pt1_N3incl" ].comp[isZ]->fill(pT1);
191            _plots["Pt3_N3incl" ].comp[isZ]->fill(pT3);
192            _plots["Rap3_N3incl"].comp[isZ]->fill(rap3);
193            _plots["ST_N3incl"  ].comp[isZ]->fill(ST);
194
195            if (njets == 3)
196              _plots["ST_N3excl"].comp[isZ]->fill(ST);
197          }
198        }
199      }
200    }
201
202
203    /// Helper for histogram initialisation
204    void hInit(string label, string ident) {
205      string pre = ident + "-x0";
206      _plots[label].ref = pre;
207      book(_plots[label].ratio, pre + "1" + _suff);
208    }
209
210    /// @}
211
212
213  protected:
214
215    // Data members
216    size_t _mode;
217    string _suff;
218
219
220  private:
221
222    /// @name Map of histograms
223    map<string, Plots> _plots;
224
225  };
226
227
228  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1312627);
229
230}