1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/VisibleFinalState.hh"
  6#include "Rivet/Projections/JetShape.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief CDF Run II jet shape analysis
 12  ///
 13  /// @author Andy Buckley
 14  class CDF_2005_I682179 : public Analysis {
 15  public:
 16
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2005_I682179);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    void init() {
 24      // Set up projections
 25      const FinalState fs(Cuts::abseta < 2);
 26      declare(fs, "FS");
 27      FastJets fj(fs, JetAlg::CDFMIDPOINT, 0.7, JetMuons::ALL, JetInvisibles::ALL);
 28      declare(fj, "Jets");
 29
 30      // Register a jet shape projection and histogram for each pT bin
 31      for (size_t i = 0; i < 6; ++i) {
 32        for (size_t j = 0; j < 3; ++j) {
 33          const size_t k = i*3 + j;
 34          stringstream ss; ss << "JetShape" << k;
 35          const string pname = ss.str();
 36          _jsnames_pT[k] = pname;
 37          const JetShape jsp(fj, 0.0, 0.7, 7, PTEDGES[k], PTEDGES[k+1], 0.1, 0.7, RAPIDITY);
 38          declare(jsp, pname);
 39          book(_profhistRho_pT[k], i+1, 1, j+1);
 40          book(_profhistPsi_pT[k], 6+i+1, 1, j+1);
 41        }
 42      }
 43
 44      // Final histo
 45      book(_profhistPsi_vs_pT, 13, 1, 1);
 46    }
 47
 48
 49    /// Do the analysis
 50    void analyze(const Event& evt) {
 51
 52      // Get jets and require at least one to pass pT and y cuts
 53      const Jets jets = apply<FastJets>(evt, "Jets")
 54        .jetsByPt(Cuts::ptIn(PTEDGES.front()*GeV, PTEDGES.back()*GeV) && Cuts::absrap < 0.7);
 55      MSG_DEBUG("Jet multiplicity before cuts = " << jets.size());
 56      if (jets.size() == 0) {
 57        MSG_DEBUG("No jets found in required pT and rapidity range");
 58        vetoEvent;
 59      }
 60
 61      // Calculate and histogram jet shapes
 62      for (size_t ipt = 0; ipt < 18; ++ipt) {
 63        const JetShape& jsipt = apply<JetShape>(evt, _jsnames_pT[ipt]);
 64        for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) {
 65          for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
 66            const double r_rho = jsipt.rBinMid(rbin);
 67            MSG_DEBUG(ipt << " " << rbin << " (" << r_rho << ") " << jsipt.diffJetShape(ijet, rbin));
 68            /// @note Bin width Jacobian factor of 0.7/0.1 = 7 in the differential shapes plot
 69            _profhistRho_pT[ipt]->fill(r_rho/0.7, (0.7/0.1)*jsipt.diffJetShape(ijet, rbin));
 70            const double r_Psi = jsipt.rBinMax(rbin);
 71            _profhistPsi_pT[ipt]->fill(r_Psi/0.7, jsipt.intJetShape(ijet, rbin));
 72          }
 73        }
 74      }
 75
 76    }
 77
 78
 79    // Finalize
 80    void finalize() {
 81      // Construct final 1-Psi(0.3/0.7) profile from Psi profiles
 82      for (auto& e : _profhistPsi_vs_pT->bins()) {
 83        // Get entry for rad_Psi = 0.2 bin
 84        const auto& b = _profhistPsi_pT[e.index()]->bin(3);
 85        if (!b.effNumEntries()) continue;
 86        e.set(b.yMean(), b.yStdErr());
 87      }
 88    }
 89
 90    /// @}
 91
 92
 93  private:
 94
 95    /// Jet \f$ p_\perp\f$ bins.
 96    const array<double, 19> PTEDGES =
 97      {{ 37.0, 45.0, 55.0, 63.0, 73.0, 84.0, 97.0, 112.0, 128.0, 148.0,
 98         166.0, 186.0, 208.0, 229.0, 250.0, 277.0, 304.0, 340.0, 380.0 }};
 99
100    /// JetShape projection name for each \f$p_\perp\f$ bin.
101    array<string, 18> _jsnames_pT;
102
103
104    /// @name Histograms
105    /// @{
106    Profile1DPtr _profhistRho_pT[18];
107    Profile1DPtr _profhistPsi_pT[18];
108    Estimate1DPtr _profhistPsi_vs_pT;
109    /// @}
110
111  };
112
113
114
115  RIVET_DECLARE_ALIASED_PLUGIN(CDF_2005_I682179, CDF_2005_S6217184);
116
117}