1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/MissingMomentum.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/WFinder.hh"
  6
  7namespace Rivet {
  8
  9
 10  ///@brief: Electroweak Wjj production at 8 TeV
 11  class ATLAS_2017_I1517194 : public Analysis {
 12  public:
 13
 14    /// @name Constructors etc.
 15    //@{
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1517194);
 19    //@}
 20
 21    /// @name Analysis methods
 22    //@{
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Get options from the new option system
 27      _mode = 0;
 28      if ( getOption("LMODE") == "EL" ) _mode = 0;
 29      if ( getOption("LMODE") == "MU" ) _mode = 1;
 30
 31      const FinalState fs;
 32      // W Selection
 33      WFinder wfinder(fs, Cuts::rap < 2.5 && Cuts::pT >= 25*GeV, _mode? PID::MUON : PID::ELECTRON, 0*GeV, 13*TeV, 0*GeV, 0.1);
 34      declare(wfinder, "WFinder");
 35
 36      FastJets jets( wfinder.remainingFinalState(), FastJets::ANTIKT, 0.4, JetAlg::Muons::DECAY, JetAlg::Invisibles::ALL);
 37      declare(jets, "Jets_w");
 38
 39      MissingMomentum missmom(FinalState(Cuts::eta < 5.0));
 40      declare(missmom, "mm");
 41
 42      const vector<string> phase_spaces = { "highmass15", "antiLC", "signal10",
 43                                            "highmass10", "inclusive", "highmass20",
 44                                            "antiLCantiJC", "antiJC", "signal",  };
 45      const vector<string> variables = { "dijetmass", "dijetpt", "dphi12", "dy12", "j1pt", "JC", "LC", "ngapjets" };
 46      size_t hepdataid = 10;
 47      for (size_t group = 0; group < 4; ++group) {
 48        for ( size_t ps=0; ps < phase_spaces.size(); ++ps ) {
 49          for ( size_t var=0; var < variables.size(); ++var ) {
 50            if (group < 2) {
 51              if ((ps == 0 || ps == 2 || ps == 3 || ps == 5) && var == 0)  continue;
 52              if ((ps == 1 || ps == 2 || ps  > 5) && var  > 4)  continue;
 53              if (group == 1 && ps == 7 && var == 3)  continue;
 54            }
 55            else {
 56              if (ps == 1 || ps == 4 || ps > 5)  continue;
 57              if ((ps == 0 || ps == 5) && var  < 2)  continue;
 58              if (ps == 2 && var  > 4)  continue;
 59              if ((ps == 0 || ps == 3 || ps == 5) && var == 5)  continue;
 60              if (group == 2) {
 61                if ((ps == 0 || ps == 5) && var == 3)  continue;
 62                if (ps == 3 && var == 1)  continue;
 63              }
 64              else {
 65                if ((ps == 0 || ps == 3 || ps == 5) && var == 6)  continue;
 66                if ((ps == 2 || ps == 3) && var  < 2)  continue;
 67              }
 68            }
 69
 70            ++hepdataid;
 71
 72            string label = variables[var]+"_"+phase_spaces[ps];
 73            //string pre = group > 1? "ew_" : "";
 74            //std::cout << "rivet -- " << pre << label << suff << " " << hepdataid << std::endl;
 75            string suff = group % 2? "" : "_norm";
 76            if (group > 1)  book(_hists["ew_" + label + suff], hepdataid, 1, 1);
 77            else            book(_hists[label + suff], hepdataid, 1, 1);
 78          }
 79        }
 80      }
 81    }
 82
 83
 84    /// Perform the per-event analysis
 85    void analyze(const Event& event) {
 86
 87      FourMomentum boson, lepton, neutrino;
 88      const WFinder& wfinder = apply<WFinder>(event, "WFinder");
 89      const FastJets& jetpro = apply<FastJets>(event, "Jets_w");
 90      const MissingMomentum& missmom = apply<MissingMomentum>(event, "mm");
 91
 92      if ( wfinder.bosons().size() != 1 ) { vetoEvent; }
 93
 94      boson    = wfinder.bosons().front().momentum();
 95      lepton   = wfinder.constituentLeptons().front().momentum();
 96      neutrino = wfinder.constituentNeutrinos().front().momentum();
 97
 98      vector<FourMomentum> jets;
 99      for (const Jet& jet : jetpro.jetsByPt(30*GeV)) {
100        if ( fabs(jet.momentum().rapidity()) > 4.4 ) continue;
101        if ( fabs(deltaR(jet, lepton)) < 0.3 ) continue;
102        jets.push_back(jet.momentum());
103      }
104
105      if (jets.size() < 2)  vetoEvent;
106
107      double mT = sqrt( 2.0*lepton.pT()*neutrino.Et()*( 1.0-cos( lepton.phi()-neutrino.phi() ) ) );
108      double DeltaRap = fabs( jets[0].rapidity()-jets[1].rapidity() );
109      double dijet_mass = FourMomentum(jets[0]+jets[1]).mass();
110      size_t nojets = jets.size();
111
112      if (missmom.vectorEt().mod() < 25*GeV)  vetoEvent;
113      if (jets[0].pT() < 80*GeV)  vetoEvent;
114      if (jets[1].pT() < 60*GeV)  vetoEvent;
115      if (dijet_mass < 500*GeV)  vetoEvent;
116      if (DeltaRap < 2.0)  vetoEvent;
117      if (mT < 40*GeV)  vetoEvent;
118
119      // By now, the event has passed all VBF cuts
120      double DijetPt = FourMomentum(jets[0]+jets[1]).pT();
121      double dphi = fabs(jets[0].phi() - jets[1].phi());
122      double DeltaPhi = ( dphi<=pi ) ? dphi/pi : (2.*pi-dphi)/pi;
123      double jet_1_rap = jets[0].rapidity();
124      double jet_2_rap = jets[1].rapidity();
125
126      // Njets in Gap Control Regions info
127      int njetsingap = 0;
128      bool firstIsForward = jet_1_rap > jet_2_rap;
129      int jF = (firstIsForward) ? 0 : 1; // sets most positive jet to be forward index
130      int jB = (firstIsForward) ? 1 : 0; // sets most negative jet to be backward index
131      for (size_t j = 2; j < nojets; ++j) {
132        if ( (jets[j].rapidity()<jets[jF].rapidity()) && (jets[j].rapidity()>jets[jB].rapidity()) ) njetsingap++;
133      }
134
135      // Third+ Jet Centrality Cut (Vetos any event that has a jet between raps of the two leading jets)
136      bool passJC = false;
137      std::vector<double> JCvals;
138      JCvals.clear();
139      if ( nojets < 3) passJC = true;
140      else if ( nojets > 2 ) {
141        passJC = true;
142        for (size_t j = 2; j < nojets; ++j) {
143          double jet_3_rap = jets[j].rapidity();
144          double jet3gap = fabs(( jet_3_rap - ((jet_1_rap + jet_2_rap)/2.0))/(jet_1_rap - jet_2_rap));
145          JCvals.push_back(jet3gap);
146          if ( jet3gap < 0.4 ) { passJC = false; }
147        }
148      }
149
150      double lepton_rap = lepton.rapidity();
151      double lep_cent = fabs((lepton_rap - ((jet_1_rap + jet_2_rap)/2.0) )/(jet_1_rap - jet_2_rap));
152      bool passLC = (lep_cent < 0.4);
153
154      map<string, bool> phaseSpaces;
155      phaseSpaces["inclusive"] = true;
156      phaseSpaces["highmass10"] = (dijet_mass>1000.0*GeV);
157      phaseSpaces["highmass15"] = (dijet_mass>1500.0*GeV);
158      phaseSpaces["highmass20"] = (dijet_mass>2000.0*GeV);
159      phaseSpaces["antiLC"] = ( !passLC && passJC );
160      phaseSpaces["antiJC"] = ( passLC && !passJC );
161      phaseSpaces["antiLCantiJC"] = ( !passLC && !passJC );
162      phaseSpaces["signal"] = ( passLC && passJC );
163      phaseSpaces["signal10"] = ( (dijet_mass>1000.0*GeV) && passLC && passJC );
164
165      for (const auto& ps : phaseSpaces) {
166        if (!ps.second)  continue;
167        fillHisto("dijetmass_"+ps.first, dijet_mass);
168        fillHisto("dijetpt_"+ps.first, DijetPt);
169        fillHisto("dy12_"+ps.first, fabs(jet_1_rap-jet_2_rap));
170        fillHisto("dphi12_"+ps.first, DeltaPhi);
171        fillHisto("j1pt_"+ps.first, jets[0].pT());
172        if (ps.first == "inclusive" || ps.first.find("highmass") != string::npos) {
173          fillHisto("LC_"+ps.first, lep_cent);
174          fillHisto("ngapjets_"+ps.first, njetsingap);
175          for (auto& jc : JCvals) {
176            fillHisto("JC_"+ps.first, jc);
177          }
178        }
179      }
180
181    }
182
183
184    /// Normalise histograms etc., after the run
185    void finalize() {
186      double factor = crossSection()/sumOfWeights()/femtobarn;  // refData is in fb
187      for (const auto& key_hist : _hists) {
188        scale(key_hist.second, factor);
189        if (key_hist.first.find("_norm") != string::npos) normalize(key_hist.second);
190      }
191    }
192
193    //@}
194
195
196    void fillHisto(const string& label, const double value) {
197      if (_hists.find(label) != _hists.end()) { // QCD+EW, absolute
198        _hists[label]->fill(value);
199      }
200      if (_hists.find(label + "_norm") != _hists.end()) { // QCD+EW, normalised
201        _hists[label + "_norm"]->fill(value);
202      }
203      if (_hists.find("ew_" + label) != _hists.end()) { // EW-only, absolute
204        _hists["ew_" + label]->fill(value);
205      }
206      if (_hists.find("ew_" + label + "_norm") != _hists.end()) { // EW-only, normalised
207        _hists["ew_" + label + "_norm"]->fill(value);
208      }
209    }
210
211
212  protected:
213
214    size_t _mode;
215
216  private:
217
218    map<string, Histo1DPtr> _hists;
219
220  };
221
222
223  // The hook for the plugin system
224  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1517194);
225
226
227}