1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/InvisibleFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/LeptonFinder.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief ATLAS pTmiss+gamma measurement at 13 TeV
 13  class ATLAS_2018_I1698006 : public Analysis {
 14  public:
 15
 16    /// Default constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1698006);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22    void init() {
 23
 24      // Get options
 25      // Default (OFF) uses the extended phase space of the paper.
 26      // LVETO ON will add an additonal lepton veto to reject non-Zgamma events
 27      // (for conservative BSM limits, for example).
 28      _mode = 0;
 29      if ( getOption("LVETO") == "ON" ) _mode = 1;
 30
 31      //prompt photons
 32      const Cut photoncut = Cuts::abspid == PID::PHOTON && Cuts::Et > 150*GeV && Cuts::abseta < 2.37;
 33      const PromptFinalState photon_fs(photoncut);
 34      declare(photon_fs, "Photons");
 35
 36      //missing energy (prompt neutrinos)
 37      declare(InvisibleFinalState(OnlyPrompt::YES), "MET");
 38
 39      if (_mode==1) {
 40        FinalState allLeps(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
 41        FinalState photons(Cuts::abspid == PID::PHOTON);
 42        PromptFinalState promptLeps(allLeps);
 43        Cut dressedLep_cuts = (Cuts::abseta < 2.7) && (Cuts::pT > 7*GeV);
 44        LeptonFinder dressedLeps(promptLeps, photons, 0.1, dressedLep_cuts);
 45        declare(dressedLeps, "dressedLeptons");
 46      }
 47
 48      //jets. run the jet finder on a final state without the prompt photons, and without neutrinos or muons
 49      VetoedFinalState jet_fs(Cuts::abseta > 4.5);
 50      jet_fs.addVetoOnThisFinalState(photon_fs);
 51      FastJets fastjets(jet_fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 52      declare(fastjets, "Jets");
 53
 54
 55      //books histograms
 56      //fig.4a
 57      book(_h["Et_inc"],2,1,1);
 58      //fig.4b
 59      book(_h["Et_exc"],3,1,1);
 60      //fig.5a
 61      book(_h["pT_inc"],4,1,1);
 62      //fig.5b
 63      book(_h["pT_exc"],5,1,1);
 64      //fig.6
 65      book(_h["Njets"],6,1,1);
 66
 67    }
 68
 69
 70    void analyze(const Event& event) {
 71
 72      const Particles& photons = apply<PromptFinalState>(event,"Photons").particlesByPt();
 73      const Jets& jets = apply<FastJets>(event,"Jets").jetsByPt(Cuts::pT > 50*GeV);
 74      const FinalState& metfs = apply<InvisibleFinalState>(event,"MET");
 75      Vector3 met_vec;
 76      for (const Particle& p : metfs.particles()) met_vec += p.mom().perpVec();
 77
 78
 79      if (_mode==1) {
 80	const DressedLeptons &dressedLeptons = apply<LeptonFinder>(event, "dressedLeptons").dressedLeptons();
 81	if (dressedLeptons.size() > 0) vetoEvent;
 82      }
 83
 84      //Nγ==1 and Emiss > 150 GeV
 85      if (met_vec.mod() > 150*GeV && photons.size()==1){
 86
 87	  //inclusive case (Njet>=0)
 88	  //if (jets.size()>=0){ //< size is always >= 0
 89          {
 90	    bool dR_veto = any(jets, DeltaRLess(photons[0], 0.3));
 91	    if (not dR_veto) {
 92	      double Et_photon = photons[0].Et()/GeV;
 93	      _h["Et_inc"]->fill(Et_photon);
 94	      //fill in missing energy (neutrino) pT histogram (inclusive)
 95	      _h["pT_inc"]->fill(met_vec.mod()/GeV);
 96	    }
 97	  }
 98
 99	  //exclusive case (Njet==0)
100	  if (jets.size() == 0){
101	    double Et_photon = photons[0].Et()/GeV;
102	    _h["Et_exc"]->fill(Et_photon);
103	     //fill in missing energy (neutrino) pT histogram (exclusive)
104	    _h["pT_exc"]->fill(met_vec.mod()/GeV);
105	  }
106
107	  _h["Njets"]->fill(jets.size());
108
109      }
110
111    }
112
113
114    void finalize() {
115      const double sf = crossSection()/femtobarn/sumOfWeights();
116      scale(_h, sf);
117    }
118
119    /// @}
120
121
122  private:
123
124    map<string, Histo1DPtr> _h;
125
126    size_t _mode;
127
128  };
129
130
131  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1698006);
132
133}