1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief CDF leading track underlying event at 300, 900 and 1960 GeV
  9  ///
 10  /// @author Orestes Tumbarell Aranda (Havana), Hannes Jung (DESY)
 11  class CDF_2015_I1388868 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2015_I1388868);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Energy selection
 25      int isqrts = -1;
 26      if (isCompatibleWithSqrtS(300*GeV)) {
 27        isqrts = 2;
 28      } else if (isCompatibleWithSqrtS(900*GeV)) {
 29        isqrts = 1;
 30      } else if (isCompatibleWithSqrtS(1960*GeV)) {
 31        isqrts = 0;
 32      } else {
 33        throw UserError("Unexpected sqrtS ! Only 300, 900, 1960 GeV is supported by CDF_2015_I1388868");
 34      }
 35      MSG_DEBUG("CDF Tevatron UE: running with " << sqrtS()/GeV);
 36
 37      // Book projection
 38      const ChargedFinalState cfs(Cuts::abseta < 0.8 && Cuts::pT > 0.5*GeV);
 39      declare(cfs, "Tracks");
 40
 41      // Book profile histos
 42      book(_NchgPDFden1,  8*isqrts+4,1,1);
 43      book(_NchgPMNden1,  8*isqrts+2,1,1);
 44      book(_NchgPMXden1,  8*isqrts+1,1,1);
 45      book(_NchgPden1,    8*isqrts+3,1,1);
 46      book(_PTsumPDFden1, 8*isqrts+8,1,1);
 47      book(_PTsumPMNden1, 8*isqrts+6,1,1);
 48      book(_PTsumPMXden1, 8*isqrts+5,1,1);
 49      book(_PTsumPden1,   8*isqrts+7,1,1);
 50
 51    }
 52
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56
 57      // Require at least one track in the event with pT >= 0.5 GeV
 58      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "Tracks");
 59      if (cfs.empty()) vetoEvent;
 60      const Particles trks = cfs.particlesByPt();
 61
 62      // Get lead track
 63      const Particle p_lead = trks[0];
 64      const double philead = p_lead.phi();
 65      const double ptlead  = p_lead.pT();
 66
 67      // Loop over tracks and compute variables
 68      double NchgP1 = 0, NchgP2 = 0, PTsumP1 = 0, PTsumP2 = 0;
 69      for (const Particle& p : trks) {
 70
 71        // Region definition -- if not in transverse region, ignore
 72        const double dphi = mapAngle0To2Pi(p.phi() - philead);
 73        if (!inRange(dphi, PI/3, 2*PI/3) && !inRange(dphi, 4*PI/3, 5*PI/3)) continue;
 74
 75        // Transverse region 1
 76        if (inRange(dphi, PI/3, 2*PI/3)) {
 77          NchgP1 += 1;
 78          PTsumP1 += p.pT();
 79        }
 80        // Transverse region 2
 81        else if (inRange(dphi, 4*PI/3, 5*PI/3)) {
 82          NchgP2 += 1;
 83          PTsumP2 += p.pT();
 84        }
 85      }
 86
 87      // Calculate total variables
 88      const double NchgPtot = (NchgP1 + NchgP2)/2;
 89      const double NchgPmax = max(NchgP1,NchgP2);
 90      const double NchgPmin = min(NchgP1,NchgP2);
 91      const double PTsumPtot = (PTsumP1 + PTsumP2)/2;
 92      const double PTsumPmax = max(PTsumP1,PTsumP2);
 93      const double PTsumPmin = min(PTsumP1,PTsumP2);
 94      //
 95      const double PTsumPMXden = PTsumPmax/AREA;
 96      const double PTsumPMNden = PTsumPmin/AREA;
 97      const double NchgPMXden = NchgPmax/AREA;
 98      const double NchgPMNden = NchgPmin/AREA;
 99      //
100      const double NchgPDFden = NchgPMXden - NchgPMNden;
101      const double PTsumPDFden = PTsumPMXden - PTsumPMNden;
102
103      // Fill histograms
104      _NchgPden1  ->fill(ptlead/GeV, NchgPtot/AREA);
105      _NchgPMXden1->fill(ptlead/GeV, NchgPmax/AREA);
106      _NchgPMNden1->fill(ptlead/GeV, NchgPmin/AREA);
107      _NchgPDFden1->fill(ptlead/GeV, NchgPDFden );
108      _PTsumPden1  ->fill(ptlead/GeV, PTsumPtot/AREA);
109      _PTsumPMXden1->fill(ptlead/GeV, PTsumPmax/AREA);
110      _PTsumPMNden1->fill(ptlead/GeV, PTsumPmin/AREA);
111      _PTsumPDFden1->fill(ptlead/GeV, PTsumPDFden );
112    }
113
114    /// @}
115
116
117    /// eta-phi area of the transverse region
118    constexpr static double AREA = 2*0.8 * M_PI/3;
119
120    /// Histograms
121    Profile1DPtr _NchgPden1, _NchgPMXden1,_NchgPMNden1,_NchgPDFden1,_PTsumPden1,_PTsumPMXden1,_PTsumPMNden1,_PTsumPDFden1;
122
123  };
124
125
126  RIVET_DECLARE_PLUGIN(CDF_2015_I1388868);
127
128}