1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// CMS UE charged particles vs. leading jet at 2.76 TeV
 11  class CMS_2015_I1385107 : public Analysis {
 12  public:
 13    /// Constructor
 14    CMS_2015_I1385107() : Analysis("CMS_2015_I1385107"),
 15                          ETACUT(2.0),
 16                          AREATOT(2*ETACUT * 2*M_PI),
 17                          AREA3(AREATOT / 3.),
 18                          AREA6(AREATOT / 6.)
 19    {   }
 20
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      const ChargedFinalState cfs(Cuts::abseta < 2 && Cuts::pT > 500*MeV);
 26      declare(cfs, "CFS");
 27
 28      const ChargedFinalState cfsforjet(Cuts::abseta < 2.5 && Cuts::pT > 500*MeV);
 29      const FastJets jetpro(cfsforjet, FastJets::SISCONE, 0.5);
 30      declare(jetpro, "Jets");
 31
 32      book(_h_Nch_TransAVE_vs_pT ,1, 1, 1); // Nch vs. pT_max      (TransAVE)
 33      book(_h_Sum_TransAVE_vs_pT ,2, 1, 1); // sum(pT) vs. pT_max  (TransAVE)
 34      book(_h_Nch_TransMAX_vs_pT ,3, 1, 1); // Nch vs. pT_max      (TransMAX)
 35      book(_h_Sum_TransMAX_vs_pT ,4, 1, 1); // sum(pT) vs. pT_max  (TransMAX)
 36      book(_h_Nch_TransMIN_vs_pT ,5, 1, 1); // Nch vs. pT_max      (TransMIN)
 37      book(_h_Sum_TransMIN_vs_pT ,6, 1, 1); // sum(pT) vs. pT_max  (TransMIN)
 38      book(_h_Nch_TransDIF_vs_pT ,7, 1, 1); // Nch vs. pT_max      (TransDIF)
 39      book(_h_Sum_TransDIF_vs_pT ,8, 1, 1); // sum(pT) vs. pT_max  (TransDIF)
 40    }
 41
 42
 43    /// Local definition of a signed dphi, for use in differentating L and R trans regions
 44    double signedDeltaPhi(double jetphi, double partphi) {
 45      double delta = partphi - jetphi;
 46      while (delta <= -PI) delta += 2 * PI;
 47      while (delta > PI) delta -= 2 * PI;
 48      return delta;
 49    }
 50
 51    /// Perform the per-event analysis
 52    void analyze(const Event& event) {
 53
 54      // Find the lead jet, applying a restriction that the jets must be within |eta| < 2.
 55      FourMomentum p_lead;
 56      for (const Jet& j : apply<FastJets>(event, "Jets").jetsByPt(1*GeV)) {
 57        if (j.abseta() < 2.0) {
 58          p_lead = j.momentum();
 59          break;
 60        }
 61      }
 62      if (p_lead.isZero()) vetoEvent;
 63      const double phi_lead = p_lead.phi();
 64      const double pT_lead  = p_lead.pT();
 65
 66      // Loop on charged particles and separate Left and Right transverse regions
 67      Particles particles = apply<ChargedFinalState>(event, "CFS").particlesByPt();
 68      int nch_TransLeft = 0, nch_TransRight = 0;
 69      double ptSum_TransLeft = 0., ptSum_TransRight = 0.;
 70      for (const Particle& p : particles) {
 71        const double dphi = signedDeltaPhi(phi_lead, p.momentum().phi());
 72        if (!inRange(fabs(dphi), PI/3, 2*PI/3.)) continue; //< only fill trans regions
 73        if (dphi < 0) {  // Transverse Right region
 74          nch_TransRight += 1;
 75          ptSum_TransRight += p.pT() / GeV;
 76        } else if (dphi > 0) {  // Transverse Left region
 77          nch_TransLeft += 1;
 78          ptSum_TransLeft += p.pT() / GeV;
 79        }
 80      }
 81
 82      // Translate to min and max (+sum and diff) Transverse regions
 83      const int nch_TransMIN = std::min(nch_TransLeft, nch_TransRight);
 84      const int nch_TransMAX = std::max(nch_TransLeft, nch_TransRight);
 85      const int nch_TransSUM = nch_TransMAX + nch_TransMIN;
 86      const int nch_TransDIF = nch_TransMAX - nch_TransMIN;
 87      //
 88      const double ptSum_TransMIN = std::min(ptSum_TransLeft, ptSum_TransRight);
 89      const double ptSum_TransMAX = std::max(ptSum_TransLeft, ptSum_TransRight);
 90      const double ptSum_TransSUM = ptSum_TransMAX + ptSum_TransMIN;
 91      const double ptSum_TransDIF = ptSum_TransMAX - ptSum_TransMIN;
 92
 93      // Fill profiles
 94      _h_Nch_TransMIN_vs_pT->fill(pT_lead/GeV, 1/AREA6 * nch_TransMIN);
 95      _h_Sum_TransMIN_vs_pT->fill(pT_lead/GeV, 1/AREA6 * ptSum_TransMIN);
 96      //
 97      _h_Nch_TransMAX_vs_pT->fill(pT_lead/GeV, 1/AREA6 * nch_TransMAX);
 98      _h_Sum_TransMAX_vs_pT->fill(pT_lead/GeV, 1/AREA6 * ptSum_TransMAX);
 99      //
100      _h_Nch_TransAVE_vs_pT->fill(pT_lead/GeV, 1/AREA3 * nch_TransSUM);
101      _h_Sum_TransAVE_vs_pT->fill(pT_lead/GeV, 1/AREA3 * ptSum_TransSUM);
102      //
103      _h_Nch_TransDIF_vs_pT->fill(pT_lead/GeV, 1/AREA6 * nch_TransDIF);
104      _h_Sum_TransDIF_vs_pT->fill(pT_lead/GeV, 1/AREA6 * ptSum_TransDIF);
105    }
106
107
108  private:
109
110    // Data members like post-cuts event weight counters go here
111    const double ETACUT, AREATOT, AREA3, AREA6;
112
113    /// Histograms
114    Profile1DPtr _h_Nch_TransAVE_vs_pT, _h_Sum_TransAVE_vs_pT;
115    Profile1DPtr _h_Nch_TransDIF_vs_pT, _h_Sum_TransDIF_vs_pT;
116    Profile1DPtr _h_Nch_TransMIN_vs_pT, _h_Sum_TransMIN_vs_pT;
117    Profile1DPtr _h_Nch_TransMAX_vs_pT, _h_Sum_TransMAX_vs_pT;
118
119  };
120
121
122  // The hook for the plugin system
123  RIVET_DECLARE_PLUGIN(CMS_2015_I1385107);
124
125}