1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Calo-based underlying event at 900 GeV and 7 TeV in ATLAS
  9  ///
 10  /// @author Jinlong Zhang
 11  class ATLAS_2011_S8994773 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_S8994773);
 15
 16
 17    void init() {
 18      const FinalState fs500((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  500*MeV));
 19      declare(fs500, "FS500");
 20      const FinalState fslead((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  1.0*GeV));
 21      declare(fslead, "FSlead");
 22
 23      // Get an index for the beam energy
 24      isqrts = -1;
 25      if (isCompatibleWithSqrtS(900)) isqrts = 0;
 26      else if (isCompatibleWithSqrtS( 7000)) isqrts = 1;
 27      assert(isqrts >= 0);
 28
 29      // N profiles, 500 MeV pT cut
 30      book(_hist_N_transverse_500 ,1+isqrts, 1, 1);
 31      // pTsum profiles, 500 MeV pT cut
 32      book(_hist_ptsum_transverse_500 ,3+isqrts, 1, 1);
 33      // N vs. Delta(phi) profiles, 500 MeV pT cut
 34      book(_hist_N_vs_dPhi_1_500 ,13+isqrts, 1, 1);
 35      book(_hist_N_vs_dPhi_2_500 ,13+isqrts, 1, 2);
 36      book(_hist_N_vs_dPhi_3_500 ,13+isqrts, 1, 3);
 37    }
 38
 39
 40    void analyze(const Event& event) {
 41      // Require at least one cluster in the event with pT >= 1 GeV
 42      const FinalState& fslead = apply<FinalState>(event, "FSlead");
 43      if (fslead.size() < 1) {
 44        vetoEvent;
 45      }
 46
 47      // These are the particles  with pT > 500 MeV
 48      const FinalState& chargedNeutral500 = apply<FinalState>(event, "FS500");
 49
 50      // Identify leading object and its phi and pT
 51      Particles particles500 = chargedNeutral500.particlesByPt();
 52      Particle p_lead = particles500[0];
 53      const double philead = p_lead.phi();
 54      const double etalead = p_lead.eta();
 55      const double pTlead  = p_lead.pT();
 56      MSG_DEBUG("Leading object: pT = " << pTlead << ", eta = " << etalead << ", phi = " << philead);
 57
 58      // Iterate over all > 500 MeV particles and count particles and scalar pTsum in the three regions
 59      vector<double> num500(3, 0), ptSum500(3, 0.0);
 60      // Temporary histos that bin N in dPhi.
 61      // NB. Only one of each needed since binnings are the same for the energies and pT cuts
 62      Histo1D hist_num_dphi_500(refData(13+isqrts,1,1));
 63      for (const Particle& p : particles500) {
 64        const double pT = p.pT();
 65        const double dPhi = deltaPhi(philead, p.phi());
 66        const int ir = region_index(dPhi);
 67        num500[ir] += 1;
 68        ptSum500[ir] += pT;
 69
 70        // Fill temp histos to bin N in dPhi
 71        if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
 72          hist_num_dphi_500.fill(dPhi, 1);
 73        }
 74      }
 75
 76
 77      // Now fill underlying event histograms
 78      // The densities are calculated by dividing the UE properties by dEta*dPhi
 79      // -- each region has a dPhi of 2*PI/3 and dEta is two times 2.5
 80      const double dEtadPhi = (2*2.5 * 2*PI/3.0);
 81      _hist_N_transverse_500->fill(pTlead/GeV,  num500[1]/dEtadPhi);
 82      _hist_ptsum_transverse_500->fill(pTlead/GeV, ptSum500[1]/GeV/dEtadPhi);
 83
 84      // Update the "proper" dphi profile histograms
 85      // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
 86      // The values tabulated in the note are for an (undefined) signed Delta(phi) rather than
 87      // |Delta(phi)| and so differ by a factor of 2: we have to actually norm for angular range = 2pi
 88      const size_t nbins = refData(13+isqrts,1,1).numPoints();
 89      for (size_t i = 0; i < nbins; ++i) {
 90        double mean = hist_num_dphi_500.bin(i).xMid();
 91        double value = 0.;
 92        if (hist_num_dphi_500.bin(i).numEntries() > 0) {
 93          mean = hist_num_dphi_500.bin(i).xMean();
 94          value = hist_num_dphi_500.bin(i).area()/hist_num_dphi_500.bin(i).xWidth()/10.0;
 95        }
 96        if (pTlead/GeV >= 1.0) _hist_N_vs_dPhi_1_500->fill(mean, value);
 97        if (pTlead/GeV >= 2.0) _hist_N_vs_dPhi_2_500->fill(mean, value);
 98        if (pTlead/GeV >= 3.0) _hist_N_vs_dPhi_3_500->fill(mean, value);
 99      }
100
101    }
102
103
104  private:
105
106    // Little helper function to identify Delta(phi) regions
107    inline int region_index(double dphi) {
108      assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
109      if (dphi < PI/3.0) return 0;
110      if (dphi < 2*PI/3.0) return 1;
111      return 2;
112    }
113
114
115    int isqrts;
116
117    Profile1DPtr _hist_N_transverse_500;
118
119    Profile1DPtr _hist_ptsum_transverse_500;
120
121    Profile1DPtr _hist_N_vs_dPhi_1_500;
122    Profile1DPtr _hist_N_vs_dPhi_2_500;
123    Profile1DPtr _hist_N_vs_dPhi_3_500;
124
125  };
126
127
128
129  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_S8994773, ATLAS_2011_I891834);
130
131}