1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/HeavyHadrons.hh"
  7
  8namespace Rivet {
  9
 10  /// Top pair production with central jet veto
 11  class ATLAS_2012_I1094568 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1094568);
 16
 17    struct Plots {
 18      // Track which veto region this is, to match the autobooked histograms
 19      int region_index;
 20
 21      // Lower rapidity boundary or veto region
 22      double y_low;
 23      // Upper rapidity boundary or veto region
 24      double y_high;
 25
 26      double veto_Q0;
 27      double veto_Qsum;
 28
 29      // Histograms to store the veto jet pT and sum(veto jet pT) histograms.
 30      Histo1DPtr h_veto_Q0;
 31      Histo1DPtr h_veto_Qsum;
 32
 33      // Estimate1Ds for the gap fractions
 34      Estimate1DPtr gapFrac_Q0;
 35      Estimate1DPtr gapFrac_Qsum;
 36    };
 37
 38    /// Book histograms and initialise projections before the run
 39    void init() {
 40
 41      const FinalState fs(Cuts::abseta < 4.5);
 42
 43      /// Get electrons from truth record
 44      FinalState elec_fs(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && Cuts::pT > 25*GeV);
 45      declare(elec_fs, "ELEC_FS");
 46
 47      /// Get muons which pass the initial kinematic cuts:
 48      FinalState muon_fs(Cuts::abspid == PID::MUON && Cuts::abseta < 2.5 && Cuts::pT > 20*GeV);
 49      declare(muon_fs, "MUON_FS");
 50
 51      /// Get all neutrinos. These will not be used to form jets.
 52      /// We'll use the highest 2 pT neutrinos to calculate the MET
 53      IdentifiedFinalState neutrino_fs(Cuts::abseta < 4.5);
 54      neutrino_fs.acceptNeutrinos();
 55      declare(neutrino_fs, "NEUTRINO_FS");
 56
 57      // Get the jets
 58      FastJets jets(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 59      declare(fs, "jet_input");
 60      declare(jets, "JETS");
 61
 62      // get b-hadrons
 63      declare(HeavyHadrons(Cuts::pT > 5*GeV), "BHadrons");
 64
 65      // Initialise weight counter
 66      book(m_total_weight, "_total_weight");
 67
 68      // Init histogramming for the various regions
 69      m_plots[0].region_index = 1;
 70      m_plots[0].y_low = 0.0;
 71      m_plots[0].y_high = 0.8;
 72      initializePlots(m_plots[0]);
 73      //
 74      m_plots[1].region_index = 2;
 75      m_plots[1].y_low = 0.8;
 76      m_plots[1].y_high = 1.5;
 77      initializePlots(m_plots[1]);
 78      //
 79      m_plots[2].region_index = 3;
 80      m_plots[2].y_low = 1.5;
 81      m_plots[2].y_high = 2.1;
 82      initializePlots(m_plots[2]);
 83      //
 84      m_plots[3].region_index = 4;
 85      m_plots[3].y_low = 0.0;
 86      m_plots[3].y_high = 2.1;
 87      initializePlots(m_plots[3]);
 88    }
 89
 90
 91    void initializePlots(Plots& plots) {
 92      plots.veto_Q0 = 0.0;
 93      const string veto_Q0_name = "TMP/vetoJetPt_Q0_" + to_str(plots.region_index);
 94      book(plots.h_veto_Q0, veto_Q0_name, 200, 0.0, 1000.0);
 95      book(plots.gapFrac_Q0, plots.region_index, 1, 1);
 96
 97      plots.veto_Qsum = 0.0;
 98      const string veto_Qsum_name = "TMP/vetoJetPt_Qsum_" + to_str(plots.region_index);
 99      book(plots.h_veto_Qsum, veto_Qsum_name, 200, 0.0, 1000.0);
100      book(plots.gapFrac_Qsum, plots.region_index, 2, 1);
101    }
102
103
104    /// Perform the per-event analysis
105    void analyze(const Event& event) {
106
107      /// Get the various sets of final state particles
108      const Particles& elecFS = apply<FinalState>(event, "ELEC_FS").particlesByPt();
109      const Particles& muonFS = apply<FinalState>(event, "MUON_FS").particlesByPt();
110      const Particles& neutrinoFS = apply<IdentifiedFinalState>(event, "NEUTRINO_FS").particlesByPt();
111
112      // Get all jets with pT > 25 GeV and |y| < 2.4
113      Jets jets = apply<FastJets>(event, "JETS").jetsByPt(Cuts::pT > 25*GeV && Cuts::absrap < 2.4);
114
115      // For each of the jets that pass the rapidity cut, only keep those that are not
116      // too close to any leptons
117      idiscardIfAnyDeltaRLess(jets, elecFS, 0.4);
118      idiscardIfAnyDeltaRLess(jets, muonFS, 0.4);
119
120      // Get b hadrons with pT > 5 GeV
121      const Particles& bHadrons = apply<HeavyHadrons>(event, "BHadrons").bHadrons();
122
123      // For each of the good jets, check whether any are b-jets (via dR matching)
124      size_t nMatches = 0;
125      Jets bJets, vetoJets;
126      for (const Jet& jet : jets) {
127        bool isBjet = any(bHadrons, DeltaRLess(jet, 0.3));
128        if (isBjet) { ++nMatches; bJets += jet; }
129        if (!isBjet || nMatches > 2)  vetoJets += jet;
130      }
131
132      // Get the MET by taking the vector sum of all neutrinos
133      /// @todo Use MissingMomentum instead?
134      double MET = 0;
135      FourMomentum p_MET;
136      for(const Particle& p: neutrinoFS) {
137        p_MET = p_MET + p.momentum();
138      }
139      MET = p_MET.pT();
140
141      // Now we have everything we need to start doing the event selections
142      bool passed_ee = false;
143
144      // We want exactly 2 electrons...
145      if (elecFS.size() == 2) {
146        // ... with opposite sign charges.
147        if (charge(elecFS[0]) != charge(elecFS[1])) {
148          // Check the MET
149          if (MET >= 40*GeV) {
150            // Do some dilepton mass cuts
151            const double dilepton_mass = (elecFS[0].momentum() + elecFS[1].momentum()).mass();
152            if (dilepton_mass >= 15*GeV) {
153              if (fabs(dilepton_mass - 91.0*GeV) >= 10.0*GeV) {
154                // We need at least 2 b-jets
155                passed_ee = bJets.size() > 1;
156              }
157            }
158          }
159        }
160      }
161
162      bool passed_mumu = false;
163      // Now do the same checks for the mumu channel
164      // So we now want 2 good muons...
165      if (muonFS.size() == 2) {
166        // ...with opposite sign charges.
167        if (charge(muonFS[0]) != charge(muonFS[1])) {
168          // Check the MET
169          if (MET >= 40*GeV) {
170            // and do some di-muon mass cuts
171            const double dilepton_mass = (muonFS.at(0).momentum() + muonFS.at(1).momentum()).mass();
172            if (dilepton_mass >= 15*GeV) {
173              if (fabs(dilepton_mass - 91.0*GeV) >= 10.0*GeV) {
174                // Need at least 2 b-jets
175                passed_mumu = bJets.size() > 1;
176              }
177            }
178          }
179        }
180      }
181
182      bool passed_emu = false;
183      // Finally, the same again with the emu channel
184      // We want exactly 1 electron and 1 muon
185      if (elecFS.size() == 1 && muonFS.size() == 1) {
186        // With opposite sign charges
187        if (charge(elecFS[0]) != charge(muonFS[0])) {
188          // Calculate HT: scalar sum of the pTs of the leptons and all good jets
189          double HT = sum(jets, pT, 0.);
190          HT += elecFS[0].pT();
191          HT += muonFS[0].pT();
192          // Keep events with HT > 130 GeV
193          if (HT > 130.0*GeV) {
194            // And again we want 2 or more b-jets
195            passed_emu = bJets.size() > 1;
196          }
197        }
198      }
199
200      if (passed_ee || passed_mumu || passed_emu) {
201        // If the event passes the selection, we use it for all gap fractions
202        m_total_weight->fill();
203
204        // Loop over each veto jet
205        for (const Jet& j : vetoJets) {
206          const double pt = j.pT();
207          const double rapidity = j.absrap();
208          // Loop over each region
209          for (size_t i = 0; i < 4; ++i) {
210            // If the jet falls into this region, get its pT and increment sum(pT)
211            if (inRange(rapidity, m_plots[i].y_low, m_plots[i].y_high)) {
212              m_plots[i].veto_Qsum += pt;
213              // If we've already got a veto jet, don't replace it
214              if (m_plots[i].veto_Q0 == 0.0)  m_plots[i].veto_Q0 = pt;
215            }
216          }
217        }
218        for (size_t i = 0; i < 4; ++i) {
219          m_plots[i].h_veto_Q0->fill(m_plots[i].veto_Q0);
220          m_plots[i].h_veto_Qsum->fill(m_plots[i].veto_Qsum);
221          m_plots[i].veto_Q0 = 0.0;
222          m_plots[i].veto_Qsum = 0.0;
223        }
224      }
225    }
226
227
228    /// Normalise histograms etc., after the run
229    void finalize() {
230      const double totalWeight = m_total_weight->val();
231      for (size_t i = 0; i < 4; ++i) {
232        finalizeGapFraction(totalWeight, m_plots[i].gapFrac_Q0,   m_plots[i].h_veto_Q0);
233        finalizeGapFraction(totalWeight, m_plots[i].gapFrac_Qsum, m_plots[i].h_veto_Qsum);
234      }
235    }
236
237
238    /// Convert temporary histos to cumulative efficiency scatters
239    /// @todo Should be possible to replace this with a couple of YODA one-lines for diff -> integral and "efficiency division"
240    void finalizeGapFraction(const double total_weight, Estimate1DPtr gapFrac, Histo1DPtr vetoPt) {
241      // Stores the cumulative frequency of the veto jet pT histogram
242      double vetoPtWeightSum = 0.0;
243
244      // Keep track of which gap fraction point we're currently populating (#final_points != #tmp_bins)
245      size_t fgap_point = 0;
246      for (size_t i = 0; i < vetoPt->numBins(); ++i) {
247        // If we've done the last "final" point, stop
248        if (fgap_point == gapFrac->numBins())  break;
249
250        // Increment the cumulative vetoPt counter for this temp histo bin
251        /// @todo Get rid of this and use vetoPt->integral(i+1) when points and bins line up?
252        vetoPtWeightSum += vetoPt->bin(i).sumW();
253
254        // If this temp histo bin's upper edge doesn't correspond to the reference point, don't finalise the scatter.
255        // Note that points are ON the bin edges and have no width: they represent the integral up to exactly that point.
256        if ( !fuzzyEquals(vetoPt->bin(i).xMax(), gapFrac->bin(fgap_point+1).xMid()) )  continue;
257
258        // Calculate the gap fraction and its uncertainty
259        const double frac = (total_weight != 0.0) ? vetoPtWeightSum/total_weight : 0;
260        const double fracErr = (total_weight != 0.0) ? sqrt(frac*(1-frac)/total_weight) : 0;
261        gapFrac->bin(fgap_point+1).set(frac, fracErr);
262
263        ++fgap_point;
264      }
265    }
266
267
268  private:
269
270    CounterPtr m_total_weight;
271    Plots m_plots[4];
272
273  };
274
275  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1094568);
276}