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_I891834 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I891834);
 15
 16
 17    void init() {
 18      const FinalState fs500(Cuts::abseta < 2.5 && Cuts::pT >=  500*MeV);
 19      declare(fs500, "FS500");
 20      const FinalState fslead(Cuts::abseta < 2.5 && Cuts::pT >=  1.0*GeV);
 21      declare(fslead, "FSlead");
 22
 23      for (double eVal : allowedEnergies()) {
 24
 25        const string en = toString(int(eVal));
 26        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
 27        bool is7TeV(en == "7000");
 28
 29        // N profiles, 500 MeV pT cut
 30        book(_h[en+"N_transverse"], 1+is7TeV, 1, 1);
 31        // pTsum profiles, 500 MeV pT cut
 32        book(_h[en+"ptsum_transverse"], 3+is7TeV, 1, 1);
 33        // N vs. Delta(phi) profiles, 500 MeV pT cut
 34        book(_h[en+"N_vs_dPhi_1"], 13+is7TeV, 1, 1);
 35        book(_h[en+"N_vs_dPhi_2"], 13+is7TeV, 1, 2);
 36        book(_h[en+"N_vs_dPhi_3"], 13+is7TeV, 1, 3);
 37      }
 38      if (_sqs == "" && !merging()) {
 39        throw BeamError("Invalid beam energy for " + name() + "\n");
 40      }
 41    }
 42
 43
 44    void analyze(const Event& event) {
 45      // Require at least one cluster in the event with pT >= 1 GeV
 46      const FinalState& fslead = apply<FinalState>(event, "FSlead");
 47      if (fslead.size() < 1)  vetoEvent;
 48
 49      // These are the particles  with pT > 500 MeV
 50      const FinalState& chargedNeutral500 = apply<FinalState>(event, "FS500");
 51
 52      // Identify leading object and its phi and pT
 53      Particles particles500 = chargedNeutral500.particlesByPt();
 54      Particle p_lead = particles500[0];
 55      const double philead = p_lead.phi();
 56      const double etalead = p_lead.eta();
 57      const double pTlead  = p_lead.pT();
 58      MSG_DEBUG("Leading object: pT = " << pTlead << ", eta = " << etalead << ", phi = " << philead);
 59
 60      // Iterate over all > 500 MeV particles and count particles and scalar pTsum in the three regions
 61      vector<double> num500(3, 0), ptSum500(3, 0.0);
 62      // Temporary histos that bin N in dPhi.
 63      // NB. Only one of each needed since binnings are the same for the energies and pT cuts
 64      YODA::Histo1D htemp(_h[_sqs+"N_vs_dPhi_1"]->xEdges(),"tmp");
 65      for (const Particle& p : particles500) {
 66        const double pT = p.pT();
 67        const double dPhi = deltaPhi(philead, p.phi());
 68        const int ir = region_index(dPhi);
 69        num500[ir] += 1;
 70        ptSum500[ir] += pT;
 71
 72        // Fill temp histos to bin N in dPhi
 73        if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
 74          htemp.fill(dPhi);
 75        }
 76      }
 77
 78
 79      // Now fill underlying event histograms
 80      // The densities are calculated by dividing the UE properties by dEta*dPhi
 81      // -- each region has a dPhi of 2*PI/3 and dEta is two times 2.5
 82      const double dEtadPhi = (2*2.5 * 2*PI/3.0);
 83      _h[_sqs+"N_transverse"]->fill(pTlead/GeV,  num500[1]/dEtadPhi);
 84      _h[_sqs+"ptsum_transverse"]->fill(pTlead/GeV, ptSum500[1]/GeV/dEtadPhi);
 85
 86      // Update the "proper" dphi profile histograms
 87      // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
 88      // The values tabulated in the note are for an (undefined) signed Delta(phi) rather than
 89      // |Delta(phi)| and so differ by a factor of 2: we have to actually norm for angular range = 2pi
 90      for (const auto& b : htemp.bins()) {
 91        double mean = b.xMid(), value = 0.;
 92        if (b.numEntries() > 0) {
 93          mean = b.xMean();
 94          value = b.sumW()/b.xWidth()/10.0;
 95        }
 96        if (pTlead/GeV >= 1.0) _h[_sqs+"N_vs_dPhi_1"]->fill(mean, value);
 97        if (pTlead/GeV >= 2.0) _h[_sqs+"N_vs_dPhi_2"]->fill(mean, value);
 98        if (pTlead/GeV >= 3.0) _h[_sqs+"N_vs_dPhi_3"]->fill(mean, value);
 99      }
100
101    }
102
103
104    // Little helper function to identify Delta(phi) regions
105    inline int region_index(double dphi) {
106      assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
107      if (dphi < PI/3.0) return 0;
108      if (dphi < 2*PI/3.0) return 1;
109      return 2;
110    }
111
112  private:
113
114    string _sqs = "";
115
116    map<string,Profile1DPtr> _h;
117
118  };
119
120
121
122  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I891834, ATLAS_2011_S8994773);
123
124}