1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/JetShape.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief D0 Run-1 jet shapes measurement
 11  class D0_1995_I398175 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(D0_1995_I398175);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      const FinalState fs((Cuts::etaIn(-4.0, 4.0)));
 25      declare(fs, "FS");
 26      //      FastJets jets(fs, JetAlg::ANTIKT, 0.6);
 27      FastJets jets(fs, JetAlg::D0ILCONE, 1.0);
 28      jets.useInvisibles();
 29      declare(jets, "Jets");
 30
 31
 32      // Specify jets pT bins
 33      _ptedges = {{ 45.0, 70.0, 105.0, 140.0, 1800.0}};
 34
 35      // Book histograms
 36      for (size_t ptbin = 0; ptbin < 4; ++ptbin) {
 37        _jsnames_pT[ptbin] = "JetShape" + to_str(ptbin) ;
 38        const JetShape jsp(jets, 0.0, 1.0, 10, _ptedges[ptbin], _ptedges[ptbin+1], 0.0, 0.2, PSEUDORAPIDITY);
 39        declare(jsp, _jsnames_pT[ptbin]);
 40        book(	_h_Rho_pT_central[ptbin] ,ptbin+1, 1, 1);
 41      }
 42
 43      const JetShape jspfwd0(jets, 0.0, 1.0, 10, 45, 70, 2.5, 3.5, PSEUDORAPIDITY);
 44      declare(jspfwd0, "JetShapeFwd0");
 45      const JetShape jspfwd1(jets, 0.0, 1.0, 10, 70, 105, 2.5, 3.5, PSEUDORAPIDITY);
 46      declare(jspfwd1, "JetShapeFwd1");
 47      book(	_h_Rho_pT_forward[0] ,5, 1, 1);
 48      book(	_h_Rho_pT_forward[1] ,6, 1, 1);
 49
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      // Get jets and require at least one to pass pT and y cuts
 57      const Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::ptIn(_ptedges.front()*GeV, _ptedges.back()*GeV) );
 58      MSG_DEBUG("Selecting jets with pT> "<<_ptedges.front());
 59      MSG_DEBUG("Jet multiplicity before cuts = " << jets.size());
 60      if (jets.size() == 0){
 61	MSG_DEBUG("No jets found in required pT and rapidity range");
 62	vetoEvent;
 63      }
 64
 65      // Calculate and histogram jet shapes
 66      for (size_t ipt = 0; ipt < 4; ++ipt) {
 67	const JetShape& jsipt = apply<JetShape>(event, _jsnames_pT[ipt]);
 68	for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) {
 69	  for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
 70	    const double r_rho = jsipt.rBinMid(rbin);
 71	    MSG_DEBUG(ipt << " " << rbin << " (" << r_rho << ") " << jsipt.diffJetShape(ijet, rbin));
 72	    /// @note Bin width Jacobian factor of 0.7/0.1 = 7 in the differential shapes plot
 73	    //	    _profhistRho_pT[ipt]->fill(r_rho/0.7, (0.7/0.1)*jsipt.diffJetShape(ijet, rbin));
 74	    const double r_Psi = jsipt.rBinMax(rbin);
 75	    MSG_DEBUG(ipt << " " << rbin << " (" << r_rho << ") " << jsipt.intJetShape(ijet, rbin));
 76	    _h_Rho_pT_central[ipt]->fill(r_Psi/1.0, jsipt.intJetShape(ijet, rbin));
 77	  }
 78	}
 79      }
 80
 81
 82      const JetShape& jsiptfwd0 = apply<JetShape>(event, "JetShapeFwd0");
 83      for (size_t ijet = 0; ijet < jsiptfwd0.numJets(); ++ijet) {
 84	for (size_t rbin = 0; rbin < jsiptfwd0.numBins(); ++rbin) {
 85	  const double r_Psi = jsiptfwd0.rBinMax(rbin);
 86	  _h_Rho_pT_forward[0]->fill(r_Psi/1.0, jsiptfwd0.intJetShape(ijet, rbin));
 87	}
 88      }
 89
 90      const JetShape& jsiptfwd1 = apply<JetShape>(event, "JetShapeFwd1");
 91      for (size_t ijet = 0; ijet < jsiptfwd1.numJets(); ++ijet) {
 92        for (size_t rbin = 0; rbin < jsiptfwd1.numBins(); ++rbin) {
 93	  const double r_Psi = jsiptfwd1.rBinMax(rbin);
 94          _h_Rho_pT_forward[1]->fill(r_Psi/1.0, jsiptfwd1.intJetShape(ijet, rbin));
 95	}
 96      }
 97
 98
 99
100
101
102    }
103
104    /// Normalise histograms etc., after the run
105    void finalize() {
106
107      // scale(_h_YYYY, crossSection()/picobarn/sumOfWeights()); // norm to cross section
108      // normalize(_h_YYYY); // normalize to unity
109
110    }
111
112    /// @}
113
114
115  private:
116
117
118    vector<double> _ptedges;
119    string _jsnames_pT[4];
120    /// @name Histograms
121    /// @{
122    Profile1DPtr _h_Rho_pT_central[4];
123    Profile1DPtr _h_Rho_pT_forward[2];
124
125    /// @}
126
127
128  };
129
130
131
132  RIVET_DECLARE_PLUGIN(D0_1995_I398175);
133
134
135}