1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/DileptonFinder.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// CMS Z rapidity measurement
 12  class CMS_2013_I1258128 : public Analysis {
 13  public:
 14
 15    // Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1258128);
 17
 18
 19    void init() {
 20
 21      // Z finders for electrons and muons
 22      Cut cuts = Cuts::abseta < 2.1 && Cuts::pT > 20*GeV;
 23      const DileptonFinder zfe(91.2*GeV, 0.1, cuts && Cuts::abspid == PID::ELECTRON, Cuts::massIn(76*GeV, 106*GeV));
 24      const DileptonFinder zfm(91.2*GeV, 0.1, cuts && Cuts::abspid == PID::MUON, Cuts::massIn(76*GeV, 106*GeV));
 25      declare(zfe, "ZFE");
 26      declare(zfm, "ZFM");
 27
 28      // Try to get the leading photon
 29      LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < 2.5 && Cuts::pT > 40*GeV));
 30      photonfs.addParticleId(PID::PHOTON);
 31      declare(photonfs, "LeadingPhoton");
 32
 33      // Jets
 34      const FastJets jets(FinalState(Cuts::abseta < 5), JetAlg::ANTIKT, 0.5);
 35      declare(jets, "JETS");
 36
 37      // Histograms
 38      book(_hist1YZ      ,1, 1, 1);
 39      book(_hist1YJet    ,2, 1, 1);
 40      book(_hist1YSum    ,3, 1, 1);
 41      book(_hist1YDif	,4, 1, 1);
 42      book(_hist2YPhoton ,5, 1, 1);
 43      book(_hist2YJet	,6, 1, 1);
 44      book(_hist2YSum	,7, 1, 1);
 45      book(_hist2YDif	,8, 1, 1);
 46    }
 47
 48
 49    void makeZCut(const Event& event) {
 50      // Apply the Z finders and veto if no Z found
 51      const DileptonFinder& zfe = apply<DileptonFinder>(event, "ZFE");
 52      const DileptonFinder& zfm = apply<DileptonFinder>(event, "ZFM");
 53      if (zfe.empty() && zfm.empty()) vetoEvent;
 54
 55      // Choose the Z candidate
 56      const Particles& z = (!zfm.empty()) ? zfm.bosons() : zfe.bosons();
 57      const Particles& clusteredConstituents = (!zfm.empty()) ? zfm.constituents() : zfe.constituents();
 58
 59      // Insist that the Z is in a high-pT (boosted) regime
 60      if (z[0].pT() < 40*GeV) return;
 61
 62      // Build the jets
 63      const FastJets& jetfs = apply<FastJets>(event, "JETS");
 64      Jets jets = jetfs.jetsByPt(Cuts::pT > 30*GeV && Cuts::abseta < 2.4);
 65      if (jets.empty()) return;
 66
 67      // Clean the jets against the lepton candidates with a DeltaR cut of 0.5
 68      vector<const Jet*> cleanedJets;
 69      for (const Jet& j : jets) {
 70        bool isolated = true;
 71        for (const Particle& p : clusteredConstituents) {
 72          if (deltaR(p, j) < 0.5) {
 73            isolated = false;
 74            break;
 75          }
 76        }
 77        if (isolated) cleanedJets.push_back(&j);
 78      }
 79      // Require exactly 1 isolated jet
 80      if (cleanedJets.size() != 1) return;
 81
 82      // Fill histos
 83      const double yz = z[0].rapidity();
 84      const double yjet = cleanedJets[0]->momentum().rapidity();
 85      _hist1YZ->fill(fabs(yz));
 86      _hist1YJet->fill(fabs(yjet));
 87      _hist1YSum->fill(0.5*fabs(yz + yjet));
 88      _hist1YDif->fill(0.5*fabs(yz - yjet));
 89    }
 90
 91
 92    void makePhotonCut(const Event& event) {
 93        // Get the photon
 94        const FinalState& photonfs = apply<FinalState>(event, "LeadingPhoton");
 95        if (photonfs.particles().size() < 1) return;
 96        const Particle& photon = photonfs.particles().front();
 97        if (photon.pT() < 40*GeV) return;
 98        if (fabs(photon.eta()) > 1.4442 ) return;
 99
100      // Build the jets
101      const FastJets& jetfs = apply<FastJets>(event, "JETS");
102      Jets jets = jetfs.jetsByPt(Cuts::pT > 30*GeV && Cuts::abseta < 2.4);
103      if (jets.empty()) return;
104
105      // Clean the jets against the photon candidate with a DeltaR cut of 0.5
106      vector<const Jet*> cleanedJets;
107      for (const Jet& j : jets)
108        if (deltaR(photon, j) > 0.5)
109          cleanedJets.push_back(&j);
110      // Require exactly 1 jet
111      if (cleanedJets.size() != 1) return;
112
113      // Fill histos
114      const double ypho = photon.rapidity();
115      const double yjet = cleanedJets[0]->momentum().rapidity();
116      _hist2YPhoton->fill(fabs(ypho));
117      _hist2YJet->fill(fabs(yjet));
118      _hist2YSum->fill(0.5*fabs(ypho + yjet));
119      _hist2YDif->fill(0.5*fabs(ypho - yjet));
120    }
121
122
123    void analyze(const Event& event) {
124      makeZCut(event);
125      makePhotonCut(event);
126    }
127
128
129    void finalize() {
130      normalizeByContents(_hist1YZ);
131      normalizeByContents(_hist1YJet);
132      normalizeByContents(_hist1YSum);
133      normalizeByContents(_hist1YDif);
134      normalizeByContents(_hist2YPhoton);
135      normalizeByContents(_hist2YJet);
136      normalizeByContents(_hist2YSum);
137      normalizeByContents(_hist2YDif);
138    }
139
140
141    // The CMS normalization in this analysis is that the sum over bin contents
142    // is equal to 1. This function normalizes to area = area*bin_width.  /
143    // @note This is a strange definition... why?
144    void normalizeByContents(Histo1DPtr h) {
145      normalize(h, h->bin(1).xWidth());
146    }
147
148
149  private:
150
151    Histo1DPtr _hist1YZ, _hist1YJet, _hist1YSum, _hist1YDif;
152    Histo1DPtr _hist2YPhoton, _hist2YJet, _hist2YSum, _hist2YDif;
153
154  };
155
156
157  RIVET_DECLARE_PLUGIN(CMS_2013_I1258128);
158
159}