1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/JetShape.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief CDF Run II b-jet shape paper
 11  class CDF_2008_I787780 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2008_I787780);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    void init() {
 22      // Set up projections
 23      const FinalState fs((Cuts::etaIn(-3.6, 3.6)));
 24      declare(fs, "FS");
 25      FastJets jetproj(fs, JetAlg::CDFMIDPOINT, 0.7);
 26      jetproj.useInvisibles();
 27      declare(jetproj, "Jets");
 28
 29      // Book histograms and corresponding jet shape projections
 30      _ptedges = {{ 52, 80, 104, 142, 300 }};
 31      for (size_t i = 0; i < 4; ++i) {
 32        stringstream ss; ss << "JetShape" << i;
 33        const string pname = ss.str();
 34        _jsnames_pT[i] = pname;
 35        const JetShape jsp(jetproj, 0.0, 0.7, 7, _ptedges[i], _ptedges[i+1], 0.0, 0.7, RAPIDITY);
 36        declare(jsp, pname);
 37        book(_h_Psi_pT[i] ,i+1, 2, 1);
 38      }
 39      book(_h_OneMinusPsi_vs_pT, 5, 1, 1);
 40    }
 41
 42
 43    // Do the analysis
 44    void analyze(const Event& event) {
 45      const FastJets& fjs = apply<FastJets>(event, "Jets");
 46      const Jets& jets = fjs.jets(Cuts::ptIn(_ptedges.front()*GeV, _ptedges.back()*GeV) && Cuts::absrap < 0.7);
 47      if (jets.size() == 0) {
 48        MSG_DEBUG("No jets found in required pT range");
 49        vetoEvent;
 50      }
 51
 52      // Filter to just get a vector of b-jets
 53      Jets bjets;
 54      for (const Jet& j : jets) {
 55        if (j.bTagged()) bjets += j;
 56      }
 57      if (bjets.empty())  {
 58        MSG_DEBUG("No b-jet axes in acceptance");
 59        vetoEvent;
 60      }
 61
 62      // Bin b-jets in pT
 63      Jets bjets_ptbinned[4];
 64      for (const Jet& bj : bjets) {
 65        const FourMomentum pbj = bj.momentum();
 66        const int ipt = binIndex(pbj.pT(), _ptedges);
 67        if (ipt == -1) continue; ///< Out of pT range (somehow!)
 68        bjets_ptbinned[ipt] += bj;
 69      }
 70
 71      // Loop over jet pT bins and fill shape profiles
 72      for (size_t ipt = 0; ipt < 4; ++ipt) {
 73        if (bjets_ptbinned[ipt].empty()) continue;
 74        // Don't use the cached result: copy construct and calculate for provided b-jets only
 75        JetShape jsipt = apply<JetShape>(event, _jsnames_pT[ipt]);
 76        jsipt.calc(bjets_ptbinned[ipt]);
 77        for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) {
 78          for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
 79            const double r_Psi = jsipt.rBinMax(rbin);
 80            _h_Psi_pT[ipt]->fill(r_Psi/0.7, jsipt.intJetShape(ijet, rbin));
 81          }
 82        }
 83      }
 84
 85    }
 86
 87
 88    /// Finalize
 89    void finalize() {
 90
 91      // Construct final 1-Psi(0.3/0.7) profile from Psi profiles
 92      for (size_t i = 0; i < _ptedges.size()-1; ++i) {
 93        // Get entry for rad_Psi = 0.2 bin
 94        Profile1DPtr ph_i = _h_Psi_pT[i];
 95        double y  = 0; // This is to protect against exceptions
 96        double ey = 0; // thrown by YODA when calling mean and
 97        if (ph_i->bin(1).effNumEntries() > 1) { // stdErr at
 98          y = 1.0 - ph_i->bin(1).yMean();        // low stats
 99          ey= ph_i->bin(1).yStdErr();
100        }
101        _h_OneMinusPsi_vs_pT->bin(i+1).set(y, ey);
102      }
103
104    }
105
106    /// @}
107
108
109  private:
110
111    /// @name Analysis data
112    /// @{
113    /// Jet \f$ p_\perp\f$ bins.
114    vector<double> _ptedges; // This can't be a raw array if we want to initialise it non-painfully
115    /// JetShape projection name for each \f$p_\perp\f$ bin.
116    string _jsnames_pT[4];
117    /// @}
118
119
120    /// @name Histograms
121    /// @{
122    Profile1DPtr _h_Psi_pT[4];
123    Estimate1DPtr _h_OneMinusPsi_vs_pT;
124    /// @}
125
126  };
127
128
129
130  RIVET_DECLARE_ALIASED_PLUGIN(CDF_2008_I787780, CDF_2008_S7782535);
131
132}