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