1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// CMS 4-jet production at 7 TeV
 10  class CMS_2013_I1273574 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    CMS_2013_I1273574()
 15      : Analysis("CMS_2013_I1273574")
 16    {    }
 17
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      const FinalState cnfs((Cuts::etaIn(-4.7, 4.7)));
 22      declare(FastJets(cnfs, JetAlg::ANTIKT, 0.5), "Jets");
 23
 24      // Modified to match the HEPDATA record.
 25      // eta of highest pT jet
 26      //book(_h_jetetas[0]     ,1,1,1);
 27      book(_h_jetetas[0]     ,6,1,1);
 28
 29      // pt of the highest pT jet
 30      book(_h_jetpts[0]      ,2,1,1);
 31
 32      //book(_h_DeltaS         ,3,1,1);
 33      book(_h_DeltaS         ,12,1,1);
 34
 35      //book(_h_DeltaPhiSoft   ,4,1,1);
 36      book(_h_DeltaPhiSoft   ,10,1,1);
 37
 38      //book(_h_DeltaPtRelSoft ,5,1,1);
 39      book(_h_DeltaPtRelSoft ,11,1,1);
 40
 41      // eta and pT of 3rd highest pT jet
 42      //book(_h_jetetas[2]     ,6,1,1);
 43      //book(_h_jetpts[2]      ,7,1,1);
 44      book(_h_jetetas[2]     ,8,1,1);
 45      book(_h_jetpts[2]      ,4,1,1);
 46
 47      // eta and pT of 4th highest pT jet
 48      //book(_h_jetetas[3]     ,8,1,1);
 49      //book(_h_jetpts[3]      ,9,1,1);
 50      book(_h_jetetas[3]     ,9,1,1);
 51      book(_h_jetpts[3]      ,5,1,1);
 52
 53      // eta and pT of 2nd highest pT jet
 54      //book(_h_jetetas[1]     ,10,1,1);
 55      //book(_h_jetpts[1]      ,11,1,1);
 56      book(_h_jetetas[1]     ,7,1,1);
 57      book(_h_jetpts[1]      ,3,1,1);
 58
 59    }
 60
 61
 62    /// Perform the per-event analysis
 63    void analyze(const Event& event) {
 64      const Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.7);
 65      if (jets.size() < 4) vetoEvent;
 66
 67      // Ensure that there are exactly 4 jets > 20 GeV, with two above 50 GeV
 68      Jets hardjets, alljets;
 69      for (const Jet& j : jets) {
 70        if (j.pT() > 50*GeV) hardjets.push_back(j);
 71        alljets.push_back(j);
 72      }
 73      if (hardjets.size() < 2 || alljets.size() != 4) vetoEvent;
 74
 75      // Histogram pT and eta of all 4 jets
 76      for (size_t i = 0; i < 4; ++i) {
 77        _h_jetpts[i]->fill(alljets[i].pT()/GeV);
 78        _h_jetetas[i]->fill(alljets[i].eta());
 79      }
 80
 81      // Create vector sums of the hard and soft pairs of jets
 82      const FourMomentum p12 = alljets[0].momentum() + alljets[1].momentum();
 83      const FourMomentum p34 = alljets[2].momentum() + alljets[3].momentum();
 84
 85      // Fill the delta(phi) between the soft jets
 86      const double dphisoft = deltaPhi(alljets[2], alljets[3]);
 87      _h_DeltaPhiSoft->fill(dphisoft);
 88
 89      // Fill the pT balance between the soft jets
 90      const double ptbalanceSoft = p34.pT() / (alljets[2].pT() + alljets[3].pT());
 91      _h_DeltaPtRelSoft->fill(ptbalanceSoft);
 92
 93      // Fill the azimuthal angle difference between the two jet pairs
 94      const double p12p34_trans = p12.px()*p34.px() + p12.py()*p34.py();
 95      const double DeltaS = acos( p12p34_trans / p12.pT() / p34.pT() );
 96      _h_DeltaS->fill(DeltaS);
 97    }
 98
 99
100    /// Normalise histograms (mostly to cross-section)
101    void finalize() {
102      const double invlumi = crossSection()/picobarn/sumOfWeights();
103      for (size_t i = 0; i < 4; ++i) {
104        scale(_h_jetpts[i], invlumi);
105        scale(_h_jetetas[i], invlumi);
106      }
107      normalize(_h_DeltaPtRelSoft);
108      normalize(_h_DeltaPhiSoft);
109      normalize(_h_DeltaS);
110    }
111
112
113  private:
114
115    Histo1DPtr _h_jetpts[4], _h_jetetas[4];
116    Histo1DPtr _h_DeltaS, _h_DeltaPhiSoft, _h_DeltaPtRelSoft;
117
118  };
119
120
121  RIVET_DECLARE_PLUGIN(CMS_2013_I1273574);
122
123}