1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/InvMassFinalState.hh"
  8#include "Rivet/Projections/MissingMomentum.hh"
  9#include "Rivet/Projections/WFinder.hh"
 10#include "Rivet/Projections/DressedLeptons.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// @brief Differential cross-section of W bosons + jets in pp collisions at sqrt(s)=7 TeV
 16  /// @author Darin Baumgartel (darinb@cern.ch)
 17  ///
 18  /// Based on Rivet analysis originally created by Anil Singh (anil@cern.ch), Lovedeep Saini (lovedeep@cern.ch)
 19  class CMS_2014_I1303894 : public Analysis {
 20  public:
 21
 22    /// Constructor
 23    CMS_2014_I1303894()
 24      : Analysis("CMS_2014_I1303894")
 25    {   }
 26
 27
 28    // Book histograms and initialise projections before the run
 29    void init() {
 30      // Prompt leptons only, no test on nu flavour.
 31      // Projections
 32      const FinalState fs;
 33      declare(fs, "FS");
 34
 35      MissingMomentum missing(fs);
 36      declare(missing, "MET");
 37
 38      PromptFinalState pfs(fs);
 39      IdentifiedFinalState bareMuons(pfs);
 40      bareMuons.acceptIdPair(PID::MUON);
 41      DressedLeptons muonClusters(fs, bareMuons, -1); //, Cuts::open(), false, false);
 42      declare(muonClusters, "muonClusters");
 43
 44      IdentifiedFinalState neutrinos(pfs);
 45      neutrinos.acceptIdPair(PID::NU_MU);
 46      declare(neutrinos, "neutrinos");
 47
 48      VetoedFinalState jetFS(fs);
 49      jetFS.addVetoOnThisFinalState(muonClusters);
 50      jetFS.addVetoOnThisFinalState(neutrinos);
 51      jetFS.vetoNeutrinos();
 52      FastJets jetprojection(jetFS, FastJets::ANTIKT, 0.5);
 53      declare(jetprojection, "Jets");
 54
 55      // Histograms
 56      book(_histDPhiMuJet1 ,1,1,1);
 57      book(_histDPhiMuJet2 ,2,1,1);
 58      book(_histDPhiMuJet3 ,3,1,1);
 59      book(_histDPhiMuJet4 ,4,1,1);
 60
 61      book(_histEtaJet1 ,5,1,1);
 62      book(_histEtaJet2 ,6,1,1);
 63      book(_histEtaJet3 ,7,1,1);
 64      book(_histEtaJet4 ,8,1,1);
 65
 66      book(_histHT1JetInc ,9,1,1);
 67      book(_histHT2JetInc ,10,1,1);
 68      book(_histHT3JetInc ,11,1,1);
 69      book(_histHT4JetInc ,12,1,1);
 70
 71      book(_histJet30MultExc  ,13,1,1);
 72      book(_histJet30MultInc  ,14,1,1);
 73
 74      book(_histPtJet1 ,15,1,1);
 75      book(_histPtJet2 ,16,1,1);
 76      book(_histPtJet3 ,17,1,1);
 77      book(_histPtJet4 ,18,1,1);
 78
 79      // Counters
 80      book(_n_1jet, "n_1jet");
 81      book(_n_2jet, "n_2jet");
 82      book(_n_3jet, "n_3jet");
 83      book(_n_4jet, "n_4jet");
 84      book(_n_inclusivebinsummation, "n_inclusivebinsummation");
 85    }
 86
 87
 88    void analyze(const Event& event) {
 89      // Get the dressed muon
 90      const DressedLeptons& muonClusters = apply<DressedLeptons>(event, "muonClusters");
 91      int nmu = muonClusters.dressedLeptons().size();
 92      if (nmu < 1) vetoEvent;
 93      DressedLepton dressedmuon = muonClusters.dressedLeptons()[0];
 94      if (dressedmuon.momentum().abseta() > 2.1) vetoEvent;
 95      if (dressedmuon.momentum().pT() < 25.0*GeV) vetoEvent;
 96
 97      // Get the muon neutrino
 98      //const Particles& neutrinos = apply<FinalState>(event, "neutrinos").particlesByPt();
 99
100      // Check that the muon and neutrino are not decay products of tau
101      if (dressedmuon.constituentLepton().hasAncestor( PID::TAU)) vetoEvent;
102      if (dressedmuon.constituentLepton().hasAncestor(-PID::TAU)) vetoEvent;
103
104      // Get the missing momentum
105      const MissingMomentum& met = apply<MissingMomentum>(event, "MET");
106      const double ptmet = met.visibleMomentum().pT();
107      const double phimet = (-met.visibleMomentum()).phi();
108
109      // Calculate MET and MT(mu,MET), and remove events with MT < 50 GeV
110      const double ptmuon = dressedmuon.pT();
111      const double phimuon = dressedmuon.phi();
112      const double mt_mumet = sqrt(2*ptmuon*ptmet*(1.0 - cos(phimet-phimuon)));
113
114      // Remove events in MT < 50 region
115      if (mt_mumet < 50*GeV) vetoEvent;
116
117      // Loop over jets and fill pt/eta/phi quantities in vectors
118      const Jets& jets_filtered = apply<FastJets>(event, "Jets").jetsByPt(0.0*GeV);
119      vector<float> finaljet_pT_list, finaljet_eta_list, finaljet_phi_list;
120      double htjets = 0.0;
121      for (size_t ii = 0; ii < jets_filtered.size(); ++ii) {
122        // Jet pT/eta/phi
123        double jet_pt = jets_filtered[ii].pT();
124        double jet_eta = jets_filtered[ii].eta();
125        double jet_phi = jets_filtered[ii].phi();
126
127        // Kinemetic cuts for jet acceptance
128        if (fabs(jet_eta) > 2.4) continue;
129        if (jet_pt < 30.0*GeV) continue;
130        if (deltaR(dressedmuon, jets_filtered[ii]) < 0.5) continue;
131
132        // Add jet to jet list and increases the HT variable
133        finaljet_pT_list.push_back(jet_pt);
134        finaljet_eta_list.push_back(jet_eta);
135        finaljet_phi_list.push_back(jet_phi);
136        htjets += fabs(jet_pt);
137      }
138
139
140      // Filling of histograms:
141      // Fill as many jets as there are into the exclusive jet multiplicity
142      if (!finaljet_pT_list.empty())
143        _histJet30MultExc->fill(finaljet_pT_list.size());
144
145      for (size_t ij = 0; ij < finaljet_pT_list.size(); ++ij) {
146        _histJet30MultInc->fill(ij+1);
147        _n_inclusivebinsummation->fill();
148      }
149
150      if (finaljet_pT_list.size() >= 1) {
151        _histPtJet1->fill(finaljet_pT_list[0]);
152        _histEtaJet1->fill(fabs(finaljet_eta_list[0]));
153        _histDPhiMuJet1->fill(deltaPhi(finaljet_phi_list[0], phimuon));
154        _histHT1JetInc->fill(htjets);
155        _n_1jet->fill();
156      }
157
158      if (finaljet_pT_list.size() >= 2) {
159        _histPtJet2->fill(finaljet_pT_list[1]);
160        _histEtaJet2->fill(fabs(finaljet_eta_list[1]));
161        _histDPhiMuJet2->fill(deltaPhi(finaljet_phi_list[1], phimuon));
162        _histHT2JetInc->fill(htjets);
163        _n_2jet->fill();
164      }
165
166      if (finaljet_pT_list.size() >= 3) {
167        _histPtJet3->fill(finaljet_pT_list[2]);
168        _histEtaJet3->fill(fabs(finaljet_eta_list[2]));
169        _histDPhiMuJet3->fill(deltaPhi(finaljet_phi_list[2], phimuon));
170        _histHT3JetInc->fill(htjets);
171        _n_3jet->fill();
172      }
173
174      if (finaljet_pT_list.size() >=4 ) {
175        _histPtJet4->fill(finaljet_pT_list[3]);
176        _histEtaJet4->fill(fabs(finaljet_eta_list[3]));
177        _histDPhiMuJet4->fill(deltaPhi(finaljet_phi_list[3], phimuon));
178        _histHT4JetInc-> fill(htjets);
179        _n_4jet->fill();
180      }
181
182    }
183
184
185    // Finalize the histograms.
186    void finalize() {
187
188      const double inclusive_cross_section = crossSection();
189      const double norm_1jet_histo = inclusive_cross_section*dbl(*_n_1jet)/sumOfWeights();
190      const double norm_2jet_histo = inclusive_cross_section*dbl(*_n_2jet)/sumOfWeights();
191      const double norm_3jet_histo = inclusive_cross_section*dbl(*_n_3jet)/sumOfWeights();
192      const double norm_4jet_histo = inclusive_cross_section*dbl(*_n_4jet)/sumOfWeights();
193      const double norm_incmultiplicity = inclusive_cross_section*dbl(*_n_inclusivebinsummation)/sumOfWeights();
194
195      normalize(_histJet30MultExc, norm_1jet_histo);
196      normalize(_histJet30MultInc, norm_incmultiplicity);
197
198      normalize(_histPtJet1, norm_1jet_histo);
199      normalize(_histHT1JetInc, norm_1jet_histo);
200      normalize(_histEtaJet1, norm_1jet_histo);
201      normalize(_histDPhiMuJet1, norm_1jet_histo);
202
203      normalize(_histPtJet2, norm_2jet_histo);
204      normalize(_histHT2JetInc, norm_2jet_histo);
205      normalize(_histEtaJet2, norm_2jet_histo);
206      normalize(_histDPhiMuJet2, norm_2jet_histo);
207
208      normalize(_histPtJet3, norm_3jet_histo);
209      normalize(_histHT3JetInc, norm_3jet_histo);
210      normalize(_histEtaJet3, norm_3jet_histo);
211      normalize(_histDPhiMuJet3, norm_3jet_histo);
212
213      normalize(_histPtJet4, norm_4jet_histo);
214      normalize(_histHT4JetInc, norm_4jet_histo);
215      normalize(_histEtaJet4, norm_4jet_histo);
216      normalize(_histDPhiMuJet4, norm_4jet_histo);
217    }
218
219
220  private:
221
222    Histo1DPtr  _histJet30MultExc, _histJet30MultInc;
223    Histo1DPtr  _histPtJet1, _histPtJet2, _histPtJet3, _histPtJet4;
224    Histo1DPtr  _histEtaJet1, _histEtaJet2, _histEtaJet3, _histEtaJet4;
225    Histo1DPtr  _histDPhiMuJet1, _histDPhiMuJet2, _histDPhiMuJet3, _histDPhiMuJet4;
226    Histo1DPtr  _histHT1JetInc, _histHT2JetInc, _histHT3JetInc, _histHT4JetInc;
227
228    CounterPtr _n_1jet, _n_2jet, _n_3jet, _n_4jet, _n_inclusivebinsummation;
229
230  };
231
232
233  RIVET_DECLARE_PLUGIN(CMS_2014_I1303894);
234
235}