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  /// Differential multijet cross-section measurement in different variables
 10  class ATLAS_2015_I1394679 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1394679);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections here
 24      const FinalState fs;
 25      declare(fs, "FinalState");
 26      FastJets fj04(fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 27      declare(fj04, "AntiKt4jets");
 28
 29      // Histograms
 30      book(_h["pt1"] ,1, 1, 1);
 31      book(_h["pt2"] ,2, 1, 1);
 32      book(_h["pt3"] ,3, 1, 1);
 33      book(_h["pt4"] ,4, 1, 1);
 34      book(_h["HT"]  ,5, 1, 1);
 35      book(_h["M4j"] ,6, 1, 1);
 36
 37      // Histograms with different pt/m4j cuts
 38      for (size_t i_hist = 0; i_hist < 4; ++i_hist) {
 39        book(_h["M2jratio_"+to_str(i_hist)] , 7 + i_hist, 1, 1);
 40        book(_h["dPhiMin2j_"+to_str(i_hist)] ,11 + i_hist, 1, 1);
 41        book(_h["dPhiMin3j_"+to_str(i_hist)] ,15 + i_hist, 1, 1);
 42        book(_h["dYMin2j_"+to_str(i_hist)] ,19 + i_hist, 1, 1);
 43        book(_h["dYMin3j_"+to_str(i_hist)] ,23 + i_hist, 1, 1);
 44        book(_h["dYMax2j_"+to_str(i_hist)] ,27 + i_hist, 1, 1);
 45        for (size_t ygap = 0; ygap < 4; ++ygap) {
 46          book(_h["sumPtCent_"+to_str(ygap)+to_str(i_hist)] ,31 + i_hist + ygap * 4, 1, 1);
 47        }
 48      }
 49
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      const Jets& alljetskt4 = apply<FastJets>(event, "AntiKt4jets").jetsByPt(Cuts::pT > 60*GeV && Cuts::absrap < 2.8);
 57
 58      // Require 4 jets with rap < 2.8 and passing pT cuts
 59      int nJets = alljetskt4.size();
 60      if (nJets < 4) vetoEvent;
 61      Jets leadingJetskt4 = alljetskt4; leadingJetskt4.resize(4);
 62      Jet jet1(leadingJetskt4[0]), jet2(leadingJetskt4[1]), jet3(leadingJetskt4[2]), jet4(leadingJetskt4[3]);
 63      if (jet1.pT() < 100*GeV) vetoEvent;
 64      if (jet4.pT() <  64*GeV) vetoEvent;
 65
 66      // dR cut
 67      const double dRcut = 0.65;
 68      double drmin = 9999;
 69      for (int ijet = 0; ijet < 4; ++ijet) {
 70        for (int jjet = ijet + 1; jjet < 4; ++jjet) {
 71          double myDR = deltaR(alljetskt4[ijet], alljetskt4[jjet], RAPIDITY);
 72          if (myDR < drmin) drmin = myDR;
 73        }
 74      }
 75      if (drmin < dRcut) vetoEvent;
 76
 77      // Variables for calculation in loops over jets
 78      FourMomentum sum_alljets;
 79      double HT = 0; // scalar sum of pt of 4 leading jets
 80      double Mjj = 99999; // minimum combined mass of 2 jets
 81      double minDphi_ij = 999, minDphi_ijk = 999; // minimum azimuthal distance btw 2 & 3 jets
 82      double maxDrap_ij = -999;  // maximum rapidity distance btw 2 jets
 83      double minDrap_ij = 999, minDrap_ijk = 999;  // minimum rapidity distance btw 2 & 3 jets
 84      size_t maxY_i = -1, maxY_j = -1;
 85
 86      // Loop over 4 leading jets
 87      for (size_t ij = 0; ij< 4; ++ij) {
 88        Jet& jeti = leadingJetskt4.at(ij);
 89        sum_alljets += jeti.mom();
 90        HT += jeti.pT();
 91
 92        for (size_t jj = 0; jj< 4; ++jj) {
 93          if ( ij == jj )  continue;
 94          Jet& jetj = leadingJetskt4.at(jj);
 95
 96          const double auxDphi = fabs(deltaPhi(jeti, jetj));
 97          minDphi_ij = std::min(auxDphi, minDphi_ij);
 98
 99          const double auxDrap = fabs(deltaRap(jeti, jetj));
100          minDrap_ij = std::min(auxDrap, minDrap_ij);
101          if (auxDrap > maxDrap_ij) {
102            maxDrap_ij = auxDrap;
103            maxY_i = ij;
104            maxY_j = jj;
105          }
106
107          FourMomentum sum_twojets = jeti.mom() + jetj.mom();
108          Mjj = std::min(Mjj, sum_twojets.mass());
109
110          for (size_t kj = 0; kj < 4; ++kj) {
111            if (kj == ij || kj == jj) continue;
112            Jet& jetk = leadingJetskt4.at(kj);
113
114            const double auxDphi2 = auxDphi + fabs(deltaPhi(jeti, jetk));
115            minDphi_ijk = std::min(auxDphi2, minDphi_ijk);
116
117            const double auxDrap2 = auxDrap + fabs(deltaRap(jeti, jetk));
118            minDrap_ijk = std::min(auxDrap2, minDrap_ijk);
119          }
120        }
121      } //end loop over 4 leading jets
122
123
124      // Combined mass of 4 leading jets
125      const double Mjjjj = sum_alljets.mass();
126
127      // Sum of central jets pT
128      double sumpt_twojets_cent = 0; // Scalar sum pt of central jets, with different rapidity gaps
129      for (size_t ijet = 0; ijet < 4; ++ijet) {
130        if (ijet == maxY_i || ijet == maxY_j) continue; // these are the forward jets
131        sumpt_twojets_cent += leadingJetskt4.at(ijet).pT();
132      }
133
134
135      // Fill histos
136      // Mass and pT cuts in which the analysis tables are split; values are in GeV and cuts are inclusive
137      const double m4jcuts[4]   = {500, 1000, 1500, 2000};
138      const double pt1cutA[4]   = {100,  400,  700, 1000};
139      const double pt1cutB[4]   = {100,  250,  400,  550};
140      const double rapGapCut[4] = {1, 2, 3, 4};
141
142      _h["pt1"]->fill(jet1.pt());
143      _h["pt2"]->fill(jet2.pt());
144      _h["pt3"]->fill(jet3.pt());
145      _h["pt4"]->fill(jet4.pt());
146      _h["HT"] ->fill(HT);
147      _h["M4j"]->fill(Mjjjj);
148
149      for (size_t i_cut = 0; i_cut < 4; ++i_cut) {
150        const string icutstr = to_str(i_cut);
151
152        if (Mjjjj > m4jcuts[i_cut])
153          _h["M2jratio_"+icutstr]->fill( Mjj/Mjjjj );
154
155        if (jet1.pT() > pt1cutA[i_cut]) {
156          _h["dPhiMin2j_"+icutstr]->fill(minDphi_ij );
157          _h["dPhiMin3j_"+icutstr]->fill(minDphi_ijk);
158          _h["dYMin2j_"+icutstr]->fill(minDrap_ij );
159          _h["dYMin3j_"+icutstr]->fill(minDrap_ijk);
160        }
161
162        if (jet1.pt() > pt1cutB[i_cut]) {
163          _h["dYMax2j_"+icutstr]->fill( maxDrap_ij );
164          for (size_t yy = 0; yy < 4; ++yy) {
165            if (maxDrap_ij > rapGapCut[yy])
166              _h["sumPtCent_"+to_str(yy)+icutstr]->fill(sumpt_twojets_cent);
167          }
168        }
169
170      } //end loop over pt/m4j cuts
171
172    }
173
174
175
176    /// Normalise histograms etc., after the run
177    void finalize() {
178      scale(_h, crossSection()/femtobarn / sumOfWeights());
179    }
180
181    /// @}
182
183
184  private:
185
186    /// @name Histograms
187    /// @{
188    map<string, Histo1DPtr> _h;
189    /// @}
190
191  };
192
193
194  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1394679);
195
196}