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