1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief CMS underlying event in leading track events at 7 TeV
  9  /// @author Paolo Gunnellini (DESY)
 10  ///
 11  /// CMS measurement of the underlying event in "leading track" events.
 12  class CMS_2012_PAS_FSQ_12_020 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2012_PAS_FSQ_12_020);
 17
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      const ChargedFinalState cfs(Cuts::abseta < 0.8 && Cuts::pT > 0.5*GeV);
 23      declare(cfs, "Tracks");
 24
 25      book(_NchgPDFden1 ,7,1,1);
 26      book(_NchgPMNden1 ,6,1,1);
 27      book(_NchgPMXden1 ,5,1,1);
 28      book(_PTsumPDFden1,10,1,1);
 29      book(_PTsumPMNden1,9,1,1);
 30      book(_PTsumPMXden1,8,1,1);
 31    }
 32
 33
 34    /// Perform the per-event analysis
 35    void analyze(const Event& event) {
 36
 37
 38      // Require at least one track in the event with pT >= 0.5 GeV
 39      const FinalState& cfs = apply<ChargedFinalState>(event, "Tracks");
 40      if (cfs.empty()) vetoEvent;
 41      const Particles trks = cfs.particlesByPt();
 42
 43      // Identify leading track and its phi and pT
 44      const Particle p_lead = trks[0];
 45      const double philead = p_lead.momentum().phi();
 46      const double ptlead  = p_lead.momentum().pT();
 47
 48      // Loop over particles and build transverse side variables
 49      double NchgP1 = 0, NchgP2 = 0, PTsumP1 = 0, PTsumP2 = 0;
 50      for (const Particle& p : trks) {
 51
 52        // Region definition -- if not in transverse region, ignore
 53        const double dphi = mapAngle0To2Pi(p.phi() - philead);
 54        if (!inRange(dphi, PI/3, 2*PI/3) && !inRange(dphi, 4*PI/3, 5*PI/3)) continue;
 55
 56        // Transverse region 1
 57        if (inRange(dphi, PI/3, 2*PI/3)) {
 58          NchgP1 += 1;
 59          PTsumP1 += p.pT();
 60        }
 61        // Transverse region 2
 62        else if (inRange(dphi, 4*PI/3, 5*PI/3)) {
 63          NchgP2 += 1;
 64          PTsumP2 += p.pT();
 65        }
 66      }
 67
 68      // Calculate total variables
 69      // const double NchgPtot = (NchgP1 + NchgP2)/2;
 70      const double NchgPmax = max(NchgP1,NchgP2);
 71      const double NchgPmin = min(NchgP1,NchgP2);
 72      // const double PTsumPtot = (PTsumP1 + PTsumP2)/2;
 73      const double PTsumPmax = max(PTsumP1,PTsumP2);
 74      const double PTsumPmin = min(PTsumP1,PTsumP2);
 75      //
 76      const double PTsumPMXden = PTsumPmax/AREA;
 77      const double PTsumPMNden = PTsumPmin/AREA;
 78      const double NchgPMXden = NchgPmax/AREA;
 79      const double NchgPMNden = NchgPmin/AREA;
 80      //
 81      const double NchgPDFden = NchgPMXden - NchgPMNden;
 82      const double PTsumPDFden = PTsumPMXden - PTsumPMNden;
 83
 84      // Fill histograms
 85      _NchgPMXden1->fill(ptlead/GeV, NchgPmax/AREA);
 86      _NchgPMNden1->fill(ptlead/GeV, NchgPmin/AREA);
 87      _NchgPDFden1->fill(ptlead/GeV, NchgPDFden);
 88      _PTsumPMXden1->fill(ptlead/GeV, PTsumPmax/AREA);
 89      _PTsumPMNden1->fill(ptlead/GeV, PTsumPmin/AREA);
 90      _PTsumPDFden1->fill(ptlead/GeV, PTsumPDFden);
 91
 92    }
 93
 94
 95    /// eta-phi area of the transverse region
 96    constexpr static double AREA = 2*0.8 * M_PI/3;
 97
 98    /// Histograms
 99    Profile1DPtr _NchgPden1, _NchgPMXden1, _NchgPMNden1, _NchgPDFden1, _PTsumPden1, _PTsumPMXden1, _PTsumPMNden1, _PTsumPDFden1;
100
101  };
102
103
104  RIVET_DECLARE_PLUGIN(CMS_2012_PAS_FSQ_12_020);
105
106}