1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  class ATLAS_2014_I1298811 : public Analysis {
 10  public:
 11
 12
 13    ATLAS_2014_I1298811()
 14      : Analysis("ATLAS_2014_I1298811") {    }
 15
 16
 17    void init() {
 18      // Configure projections
 19      const FinalState fs((Cuts::etaIn(-4.8, 4.8)));
 20      declare(fs, "FS");
 21      const FastJets jets(fs, JetAlg::ANTIKT, 0.4);
 22      declare(jets, "Jets");
 23
 24      // Book histograms
 25      for (size_t itopo = 0; itopo < 2; ++itopo) {
 26        // Profiles
 27        for (size_t iregion = 0; iregion < 3; ++iregion) {
 28          book(_p_ptsumch_vs_ptlead[itopo][iregion] ,1+iregion, 1, itopo+1);
 29          book(_p_nch_vs_ptlead[itopo][iregion] ,4+iregion, 1, itopo+1);
 30        }
 31        book(_p_etsum25_vs_ptlead_trans[itopo] ,7, 1, itopo+1);
 32        book(_p_etsum48_vs_ptlead_trans[itopo] ,8, 1, itopo+1);
 33        book(_p_chratio_vs_ptlead_trans[itopo] ,9, 1, itopo+1);
 34        book(_p_ptmeanch_vs_ptlead_trans[itopo] ,10, 1, itopo+1);
 35        // 1D histos
 36        for (size_t iregion = 0; iregion < 3; ++iregion) {
 37          for (size_t ipt = 0; ipt < 4; ++ipt) {
 38            book(_h_ptsumch[ipt][itopo][iregion] ,13+3*ipt+iregion, 1, itopo+1);
 39            book(_h_nch[ipt][itopo][iregion] ,25+3*ipt+iregion, 1, itopo+1);
 40          }
 41        }
 42      }
 43      book(_p_ptmeanch_vs_nch_trans[0], 11, 1, 1);
 44      book(_p_ptmeanch_vs_nch_trans[1], 12, 1, 1);
 45
 46    }
 47
 48
 49
 50    void analyze(const Event& event) {
 51      // Find the jets with pT > 20 GeV and *rapidity* within 2.8
 52      /// @todo Use Cuts instead rather than an eta cut in the proj and a y cut after
 53      const Jets alljets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV);
 54      Jets jets;
 55      for (const Jet& j : alljets)
 56        if (j.absrap() < 2.8) jets.push_back(j);
 57      // Require at least one jet in the event
 58      if (jets.empty()) vetoEvent;
 59
 60      // Identify the leading jet and its phi and pT
 61      const FourMomentum plead = jets[0].momentum();
 62      const double philead = plead.phi();
 63      const double etalead = plead.eta();
 64      const double ptlead  = plead.pT();
 65      MSG_DEBUG("Leading object: pT = " << ptlead << ", eta = " << etalead << ", phi = " << philead);
 66
 67      // Sum particle properties in the transverse regions
 68      int tmpnch[2] = {0,0};
 69      double tmpptsum[2] = {0,0};
 70      double tmpetsum48[2] = {0,0};
 71      double tmpetsum25[2] = {0,0};
 72      const Particles particles = apply<FinalState>(event, "FS").particles();
 73      for (const Particle& p : particles) {
 74        // Only consider the transverse region(s), not toward or away
 75        if (!inRange(deltaPhi(p.phi(), philead), PI/3.0, TWOPI/3.0)) continue;
 76        // Work out which transverse side this particle is on
 77        const size_t iside = (mapAngleMPiToPi(p.phi() - philead) > 0) ? 0 : 1;
 78        MSG_TRACE(p.phi() << " vs. " << philead << ": " << iside);
 79        // Charged or neutral particle?
 80        const bool charged = p.charge3() != 0;
 81        // Track observables
 82        if (charged && fabs(p.eta()) < 2.5 && p.pT() > 500*MeV) {
 83          tmpnch[iside] += 1;
 84          tmpptsum[iside] += p.pT();
 85        }
 86        // Cluster observables
 87        if ((charged && p.p3().mod() > 200*MeV) || (!charged && p.p3().mod() > 500*MeV)) {
 88          tmpetsum48[iside] += p.pT();
 89          if (fabs(p.eta()) < 2.5) tmpetsum25[iside] += p.pT();
 90        }
 91      }
 92
 93      // Construct tot/max/min counts (for trans/max/min, indexed by iregion)
 94      const int nch[3] = { tmpnch[0] + tmpnch[1],
 95                           std::max(tmpnch[0], tmpnch[1]),
 96                           std::min(tmpnch[0], tmpnch[1]) };
 97      const double ptsum[3] = { tmpptsum[0] + tmpptsum[1],
 98                                std::max(tmpptsum[0], tmpptsum[1]),
 99                                std::min(tmpptsum[0], tmpptsum[1]) };
100      const double etsum48[3] = { tmpetsum48[0] + tmpetsum48[1],
101                                  std::max(tmpetsum48[0], tmpetsum48[1]),
102                                  std::min(tmpetsum48[0], tmpetsum48[1]) };
103      const double etsum25[3] = { tmpetsum25[0] + tmpetsum25[1],
104                                  std::max(tmpetsum25[0], tmpetsum25[1]),
105                                  std::min(tmpetsum25[0], tmpetsum25[1]) };
106
107
108      //////////////////////////////////////////////////////////
109      // Now fill the histograms with the computed quantities
110
111      // phi sizes of each trans/max/min region (for indexing by iregion)
112      const double dphi[3] = { 2*PI/3.0, PI/3.0, PI/3.0 };
113
114      // Loop over inclusive jet and exclusive dijet configurations
115      for (size_t itopo = 0; itopo < 2; ++itopo) {
116
117        // Exit early if in the exclusive dijet iteration and the exclusive dijet cuts are not met
118        if (itopo == 1) {
119          if (jets.size() != 2) continue;
120          const FourMomentum psublead = jets[1].momentum();
121          // Delta(phi) cut
122          const double phisublead = psublead.phi();
123          if (deltaPhi(philead, phisublead) < 2.5) continue;
124          // pT fraction cut
125          const double ptsublead  = psublead.pT();
126          if (ptsublead < 0.5*ptlead) continue;
127          MSG_DEBUG("Exclusive dijet event");
128        }
129
130        // Plot profiles and distributions which have no max/min region definition
131        _p_etsum25_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum25[0]/5.0/dphi[0]/GeV);
132        _p_etsum48_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum48[0]/9.6/dphi[0]/GeV);
133        if (etsum25[0] > 0) {
134          _p_chratio_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptsum[0]/etsum25[0]);
135        }
136        const double ptmean = safediv(ptsum[0], nch[0], -1); ///< Return -1 if div by zero
137        if (ptmean >= 0) {
138          _p_ptmeanch_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptmean/GeV);
139          _p_ptmeanch_vs_nch_trans[itopo]->fill(nch[0], ptmean/GeV);
140        }
141
142        // Plot remaining profile and 1D observables, which are defined in all 3 tot/max/min regions
143        for (size_t iregion = 0; iregion < 3; ++iregion) {
144          _p_ptsumch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, ptsum[iregion]/5.0/dphi[iregion]/GeV);
145          _p_nch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, nch[iregion]/5.0/dphi[iregion]);
146          for (size_t ipt = 0; ipt < 4; ++ipt) {
147            if (ipt == 1 && !inRange(ptlead/GeV, 20, 60)) continue;
148            if (ipt == 2 && !inRange(ptlead/GeV, 60, 210)) continue;
149            if (ipt == 3 && ptlead/GeV < 210) continue;
150            _h_ptsumch[ipt][itopo][iregion]->fill(ptsum[iregion]/5.0/dphi[iregion]/GeV);
151            _h_nch[ipt][itopo][iregion]->fill(nch[iregion]/5.0/dphi[iregion]);
152          }
153        }
154      }
155
156    }
157
158
159
160    void finalize() {
161      for (size_t iregion = 0; iregion < 3; ++iregion) {
162        for (size_t itopo = 0; itopo < 2; ++itopo) {
163          for (size_t ipt = 0; ipt < 4; ++ipt) {
164            normalize(_h_ptsumch[ipt][itopo][iregion], 1.0);
165            normalize(_h_nch[ipt][itopo][iregion], 1.0);
166          }
167        }
168      }
169    }
170
171
172  private:
173
174    /// @name Histogram arrays
175    /// @{
176
177    Profile1DPtr _p_ptsumch_vs_ptlead[2][3];
178    Profile1DPtr _p_nch_vs_ptlead[2][3];
179    Profile1DPtr _p_ptmeanch_vs_ptlead_trans[2];
180    Profile1DPtr _p_etsum25_vs_ptlead_trans[2];
181    Profile1DPtr _p_etsum48_vs_ptlead_trans[2];
182    Profile1DPtr _p_chratio_vs_ptlead_trans[2];
183    Profile1DPtr _p_ptmeanch_vs_nch_trans[2];
184    Histo1DPtr _h_ptsumch[4][2][3];
185    Histo1DPtr _h_nch[4][2][3];
186
187    /// @}
188
189  };
190
191
192
193  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1298811);
194
195}