1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/FastJets.hh"
 5#include "Rivet/Projections/ChargedFinalState.hh"
 6
 7namespace Rivet {
 8
 9
10  /// Per-event yield of the highest transverse momentum charged particle and charged-particle jet
11  class CMS_2015_I1380605 : public Analysis {
12  public:
13
14    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2015_I1380605);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22      const ChargedFinalState cfs((Cuts::etaIn(-7., 7.)));
23      declare(cfs, "CFS");
24      declare(FastJets(cfs, JetAlg::ANTIKT, 0.5), "Jets");
25
26      book(_h_tracks, 1, 1, 1);
27      book(_h_jets  , 2, 1, 1);
28      book(_ntracks, "ntracks");
29    }
30
31
32    /// Perform the per-event analysis
33    void analyze(const Event& event) {
34
35      // Veto events without forward activity on both sides
36      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
37      const size_t count_plus  = cfs.particles(Cuts::eta > 5.3  && Cuts::eta < 6.5).size();
38      const size_t count_minus = cfs.particles(Cuts::eta < -5.3 && Cuts::eta > -6.5).size();
39      if (!(count_plus > 0 || count_minus > 0)) vetoEvent; //< @todo Should this really also veto the jet analysis?
40      /// @warning Needs migration to an AO Counter
41      /// @note This isn't the number of tracks, it's the sum of event weights passing the veto
42      _ntracks->fill();
43
44      // Do track analysis here
45      // Find pttrackmax
46      double track_ptmax = 0;
47      for (const Particle& p : cfs.particles(Cuts::abseta < 2.4)) track_ptmax = max(track_ptmax, p.pT());
48      // Fill track analysis histograms
49      for (size_t i = 1; i <= _h_tracks->numBins(); ++i) {
50        const double binlimitlow_t = _h_tracks->bin(i).xMin();
51        const double weightbw_t = _h_tracks->bin(i).xWidth();
52        const double xbin_t = _h_tracks->bin(i).xMid();
53        if (track_ptmax > binlimitlow_t) _h_tracks -> fill(xbin_t, weightbw_t);
54      }
55
56      // Do jet analysis here
57      const Jets jetsdeta = apply<FastJets>(event,"Jets").jets(Cuts::pT > 1*GeV && Cuts::pT < 60*GeV && Cuts::abseta < 1.9);
58      // Find ptjetmax
59      double jet_ptmax = 0;
60      for (const Jet& j : jetsdeta) jet_ptmax = max(jet_ptmax, j.pT());
61      // Fill jet analysis histograms
62      for (size_t i = 1; i <= _h_jets->numBins(); ++i) {
63        const double binlimitlow_j = _h_jets->bin(i).xMin();
64        const double weightbw_j = _h_jets->bin(i).xWidth();
65        const double xbin_j = _h_jets->bin(i).xMid();
66        if (jet_ptmax > binlimitlow_j) _h_jets -> fill(xbin_j, weightbw_j);
67      }
68    }
69
70
71    /// Normalise histograms etc., after the run
72    void finalize() {
73      const double norm_t0 = _h_tracks->bin(8).sumW()/2.056170e-03/_h_tracks->bin(8).xWidth();
74      //const double norm_t1 = _h_tracks->bin(7).sumW()/2.056170e-03;
75      const double norm_j0 = _h_jets->bin(14).sumW()/3.575290e-03/_h_jets->bin(14).xWidth();
76      //const double norm_j1 = _h_jets->bin(13).sumW()/3.575290e-03;
77      if (norm_t0 > 0 ) scale(_h_tracks, 1./ norm_t0);
78      if (norm_j0 > 0 ) scale(_h_jets, 1./ norm_j0);
79    }
80
81    /// @}
82
83
84  private:
85
86    /// @name Histograms
87    /// @{
88    Histo1DPtr _h_tracks, _h_jets;
89    CounterPtr _ntracks;
90    /// @}
91
92  };
93
94
95  RIVET_DECLARE_PLUGIN(CMS_2015_I1380605);
96
97}