1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// ATLAS charged particle jet underlying event and jet radius dependence
 10  class ATLAS_2012_I1125575 : public Analysis {
 11  public:
 12
 13    /// @name Constructors etc.
 14    /// @{
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1125575);
 18
 19    /// @}
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28      const ChargedFinalState jet_input((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  0.5*GeV));
 29      declare(jet_input, "JET_INPUT");
 30
 31      const ChargedFinalState track_input((Cuts::etaIn(-1.5, 1.5) && Cuts::pT >=  0.5*GeV));
 32      declare(track_input, "TRACK_INPUT");
 33
 34      const FastJets jets02(jet_input, JetAlg::ANTIKT, 0.2);
 35      declare(jets02, "JETS_02");
 36
 37      const FastJets jets04(jet_input, JetAlg::ANTIKT, 0.4);
 38      declare(jets04, "JETS_04");
 39
 40      const FastJets jets06(jet_input, JetAlg::ANTIKT, 0.6);
 41      declare(jets06, "JETS_06");
 42
 43      const FastJets jets08(jet_input, JetAlg::ANTIKT, 0.8);
 44      declare(jets08, "JETS_08");
 45
 46      const FastJets jets10(jet_input, JetAlg::ANTIKT, 1.0);
 47      declare(jets10, "JETS_10");
 48
 49      // Mean number of tracks
 50      initializeProfiles(_h_meanNch, 1);
 51
 52      // Mean of the average track pT in each region
 53      initializeProfiles(_h_meanPtAvg, 2);
 54
 55      // Mean of the scalar sum of track pT in each region
 56      initializeProfiles(_h_meanPtSum, 3);
 57
 58      // Distribution of Nch, in bins of leading track-jet pT
 59      initializeHistograms(_h_Nch, 4);
 60
 61      // Distribution of average track-jet pT, in bins of leading track-jet pT
 62      initializeHistograms(_h_PtAvg, 5);
 63
 64      // Distribution of sum of track-jet pT, in bins of leading track-jet pT
 65      initializeHistograms(_h_PtSum, 6);
 66
 67      for (int i = 0; i < 5; ++i)
 68        book(_nEvents[i], "nEvents_"+to_str(i));
 69    }
 70
 71
 72    void initializeProfiles(Profile1DPtr plots[5][2], int distribution) {
 73      for (int i = 0; i < 5; ++i) {
 74        for (int j = 0; j < 2; ++j) {
 75          book(plots[i][j], distribution, i+1, j+1);
 76        }
 77      }
 78    }
 79
 80
 81    void initializeHistograms(Histo1DGroupPtr plots[5][2], int distribution) {
 82      const Estimate1D& refest = refData(1, 1, 1);
 83      for (size_t i = 0; i < 5; ++i) {
 84        for (size_t y = 0; y < 2; ++y) {
 85          book(plots[i][y], refest.xEdges());
 86          for (auto& b : plots[i][y]->bins()) {
 87            size_t histogram_number = (b.index()*2)-((y+1)%2);
 88            book(b, distribution, i+1, histogram_number);
 89          }
 90        }
 91      }
 92    }
 93
 94
 95    /// Perform the per-event analysis
 96    void analyze(const Event& event) {
 97      vector<Jets*> all_jets;
 98      Jets jets_02 = apply<FastJets>(event, "JETS_02").jetsByPt(Cuts::pT > 4*GeV && Cuts::abseta < 1.5);
 99      all_jets.push_back(&jets_02);
100      Jets jets_04 = apply<FastJets>(event, "JETS_04").jetsByPt(Cuts::pT > 4*GeV && Cuts::abseta < 1.5);
101      all_jets.push_back(&jets_04);
102      Jets jets_06 = apply<FastJets>(event, "JETS_06").jetsByPt(Cuts::pT > 4*GeV && Cuts::abseta < 1.5);
103      all_jets.push_back(&jets_06);
104      Jets jets_08 = apply<FastJets>(event, "JETS_08").jetsByPt(Cuts::pT > 4*GeV && Cuts::abseta < 1.5);
105      all_jets.push_back(&jets_08);
106      Jets jets_10 = apply<FastJets>(event, "JETS_10").jetsByPt(Cuts::pT > 4*GeV && Cuts::abseta < 1.5);
107      all_jets.push_back(&jets_10);
108
109      // Count the number of tracks in the away and transverse regions, for each set of jets
110      double n_ch[5][2] = { {0,0}, {0,0}, {0,0}, {0,0}, {0,0} };
111      // Also add up the sum pT
112      double sumpt[5][2] = { {0,0}, {0,0}, {0,0}, {0,0}, {0,0} };
113      // ptmean = sumpt / n_ch
114      double ptavg[5][2] = { {0,0}, {0,0}, {0,0}, {0,0}, {0,0} };
115      // lead jet pT defines which bin we want to fill
116      double lead_jet_pts[5] = {0.0};
117
118      // Loop over each of the jet radii:
119      for (int i = 0; i < 5; ++i) {
120        if (all_jets[i]->size() < 1) continue;
121
122        // Find the lead jet pT
123        lead_jet_pts[i] = all_jets[i]->at(0).pT();
124
125        // Loop over each of the charged particles
126        const Particles& tracks = apply<ChargedFinalState>(event, "TRACK_INPUT").particlesByPt();
127        for(const Particle& t : tracks) {
128
129          // Get the delta-phi between the track and the leading jet
130          double dphi = deltaPhi(all_jets[i]->at(0), t);
131
132          // Find out which region this puts it in.
133          // 0 = away region, 1 = transverse region, 2 = toward region
134          int region = region_index(dphi);
135
136          // If the track is in the toward region, ignore it.
137          if (region == 2) continue;
138
139          // Otherwise, increment the relevant counters
140          ++n_ch[i][region];
141          sumpt[i][region] += t.pT();
142
143        }
144        // Calculate the pT_avg for the away and transverse regions.
145        // (And make sure we don't try to divide by zero.)
146        ptavg[i][0] = (n_ch[i][0] == 0 ? 0.0 : sumpt[i][0] / n_ch[i][0]);
147        ptavg[i][1] = (n_ch[i][1] == 0 ? 0.0 : sumpt[i][1] / n_ch[i][1]);
148
149        _nEvents[i]->fill();
150      }
151
152      fillProfiles(_h_meanNch,    n_ch, lead_jet_pts, 1.0 / (2*PI));
153      fillProfiles(_h_meanPtAvg, ptavg, lead_jet_pts, 1.0);
154      fillProfiles(_h_meanPtSum, sumpt, lead_jet_pts, 1.0 / (2*PI));
155
156      fillHistograms(_h_Nch,    n_ch, lead_jet_pts);
157      fillHistograms(_h_PtAvg, ptavg, lead_jet_pts);
158      fillHistograms(_h_PtSum, sumpt, lead_jet_pts);
159    }
160
161
162    void fillProfiles(Profile1DPtr plots[5][2], double var[5][2], double lead_pt[5], double scale) {
163      for (int i=0; i<5; ++i) {
164        double pt = lead_pt[i];
165        for (int j=0; j<2; ++j) {
166          double v = var[i][j];
167          plots[i][j]->fill(pt, v*scale);
168        }
169      }
170    }
171
172
173    void fillHistograms(Histo1DGroupPtr plots[5][2], double var[5][2], double lead_pt[5]) {
174      for (int i=0; i<5; ++i) {
175        double pt = lead_pt[i];
176        for (int j=0; j<2; ++j) {
177          double v = var[i][j];
178          plots[i][j]->fill(pt, v);
179        }
180      }
181    }
182
183
184    int region_index(double dphi) {
185      assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
186      if (dphi < PI/3.0) return 2;
187      if (dphi < 2*PI/3.0) return 1;
188      return 0;
189    }
190
191
192    /// Normalise histograms etc., after the run
193    void finalize() {
194      for (size_t i = 0; i < 5; ++i) {
195        scale(_h_Nch[i], 1.0/ *_nEvents[i]);
196        scale(_h_PtAvg[i], 1.0/ *_nEvents[i]);
197        scale(_h_PtSum[i], 1.0/ *_nEvents[i]);
198      }
199    }
200
201
202    /// @}
203
204
205  private:
206
207    // Data members like post-cuts event weight counters go here
208    CounterPtr _nEvents[5];
209
210    Profile1DPtr _h_meanNch[5][2];
211    Profile1DPtr _h_meanPtAvg[5][2];
212    Profile1DPtr _h_meanPtSum[5][2];
213
214    Histo1DGroupPtr _h_Nch[5][2];
215    Histo1DGroupPtr _h_PtAvg[5][2];
216    Histo1DGroupPtr _h_PtSum[5][2];
217
218  };
219
220
221
222  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1125575);
223
224}