Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/JetShape.hh"
#include "Rivet/Projections/FastJets.hh"

namespace Rivet {


  /// @brief D0 Run-1 jet shapes measurement
  class D0_1995_I398175 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(D0_1995_I398175);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      const FinalState fs((Cuts::etaIn(-4.0, 4.0)));
      declare(fs, "FS");
      //      FastJets jets(fs, JetAlg::ANTIKT, 0.6);
      FastJets jets(fs, JetAlg::D0ILCONE, 1.0);
      jets.useInvisibles();
      declare(jets, "Jets");


      // Specify jets pT bins
      _ptedges = {{ 45.0, 70.0, 105.0, 140.0, 1800.0}};

      // Book histograms
      for (size_t ptbin = 0; ptbin < 4; ++ptbin) {
        _jsnames_pT[ptbin] = "JetShape" + to_str(ptbin) ;
        const JetShape jsp(jets, 0.0, 1.0, 10, _ptedges[ptbin], _ptedges[ptbin+1], 0.0, 0.2, PSEUDORAPIDITY);
        declare(jsp, _jsnames_pT[ptbin]);
        book(   _h_Rho_pT_central[ptbin] ,ptbin+1, 1, 1);
      }

      const JetShape jspfwd0(jets, 0.0, 1.0, 10, 45, 70, 2.5, 3.5, PSEUDORAPIDITY);
      declare(jspfwd0, "JetShapeFwd0");
      const JetShape jspfwd1(jets, 0.0, 1.0, 10, 70, 105, 2.5, 3.5, PSEUDORAPIDITY);
      declare(jspfwd1, "JetShapeFwd1");
      book( _h_Rho_pT_forward[0] ,5, 1, 1);
      book( _h_Rho_pT_forward[1] ,6, 1, 1);

    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if(_edges.empty()) _edges=_h_Rho_pT_central[0]->xEdges();
      // Get jets and require at least one to pass pT and y cuts
      const Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::ptIn(_ptedges.front()*GeV, _ptedges.back()*GeV) );
      MSG_DEBUG("Selecting jets with pT> "<<_ptedges.front());
      MSG_DEBUG("Jet multiplicity before cuts = " << jets.size());
      if (jets.size() == 0){
    MSG_DEBUG("No jets found in required pT and rapidity range");
    vetoEvent;
      }

      // Calculate and histogram jet shapes
      for (size_t ipt = 0; ipt < 4; ++ipt) {
    const JetShape& jsipt = apply<JetShape>(event, _jsnames_pT[ipt]);
    for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) {
      for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
        const double r_rho = jsipt.rBinMid(rbin);
        MSG_DEBUG(ipt << " " << rbin << " (" << r_rho << ") " << jsipt.diffJetShape(ijet, rbin));
        /// @note Bin width Jacobian factor of 0.7/0.1 = 7 in the differential shapes plot
        //      _profhistRho_pT[ipt]->fill(r_rho/0.7, (0.7/0.1)*jsipt.diffJetShape(ijet, rbin));
        MSG_DEBUG(ipt << " " << rbin << " (" << r_rho << ") " << jsipt.intJetShape(ijet, rbin));
        _h_Rho_pT_central[ipt]->fill(_edges[rbin], jsipt.intJetShape(ijet, rbin));
      }
    }
      }


      const JetShape& jsiptfwd0 = apply<JetShape>(event, "JetShapeFwd0");
      for (size_t ijet = 0; ijet < jsiptfwd0.numJets(); ++ijet) {
    for (size_t rbin = 0; rbin < jsiptfwd0.numBins(); ++rbin) {
      _h_Rho_pT_forward[0]->fill(_edges[rbin], jsiptfwd0.intJetShape(ijet, rbin));
    }
      }

      const JetShape& jsiptfwd1 = apply<JetShape>(event, "JetShapeFwd1");
      for (size_t ijet = 0; ijet < jsiptfwd1.numJets(); ++ijet) {
        for (size_t rbin = 0; rbin < jsiptfwd1.numBins(); ++rbin) {
          _h_Rho_pT_forward[1]->fill(_edges[rbin], jsiptfwd1.intJetShape(ijet, rbin));
    }
      }





    }

    /// Normalise histograms etc., after the run
    void finalize() { }

    /// @}


  private:


    vector<double> _ptedges;
    string _jsnames_pT[4];
    vector<string> _edges;
    /// @name Histograms
    /// @{
    BinnedProfilePtr<string> _h_Rho_pT_central[4];
    BinnedProfilePtr<string> _h_Rho_pT_forward[2];
    /// @}


  };



  RIVET_DECLARE_PLUGIN(D0_1995_I398175);


}