1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/InvisibleFinalState.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9
 10namespace Rivet {
 11
 12  /// @brief VBS Zy production at 13 TeV
 13  class ATLAS_2023_I2663725 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2023_I2663725);
 18
 19    /// @name Analysis methods
 20    //@{
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      FinalState fs;
 25      // Prompt photons
 26      const PromptFinalState photon_fs(Cuts::abspid == PID::PHOTON && Cuts::pT > 25 * GeV && Cuts::abseta < 2.37);
 27      declare(photon_fs, "Photons");
 28
 29      // Prompt leptons
 30      const PromptFinalState bareelectron_fs = Cuts::abspid == PID::ELECTRON;
 31      const PromptFinalState baremuon_fs = Cuts::abspid == PID::MUON;
 32
 33      // Dressed leptons
 34      const FinalState allphoton_fs(Cuts::abspid == PID::PHOTON);
 35      const Cut leptoncut = Cuts::pT > 20 * GeV && Cuts::abseta < 2.47;
 36      // use *all* photons for dressing
 37      const LeptonFinder dressedelectron_fs(bareelectron_fs, allphoton_fs, 0.1, leptoncut);
 38      const LeptonFinder dressedmuon_fs(baremuon_fs, allphoton_fs, 0.1, leptoncut);
 39
 40      declare(dressedelectron_fs, "Electrons");
 41      declare(dressedmuon_fs, "Muons");
 42
 43      // FS excluding the leading photon
 44      VetoedFinalState vfs;
 45      vfs.addVetoOnThisFinalState(photon_fs);
 46      vfs.addVetoOnThisFinalState(dressedmuon_fs);
 47      vfs.addVetoOnThisFinalState(InvisibleFinalState());
 48      declare(vfs, "isolatedFS");
 49
 50      VetoedFinalState hadrons(FinalState(Cuts::abseta < 4.4));
 51      hadrons.addVetoOnThisFinalState(dressedelectron_fs);
 52      hadrons.addVetoOnThisFinalState(dressedmuon_fs);
 53
 54      // Jets
 55      FastJets jets(hadrons, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 56      declare(jets, "jets");
 57
 58      // Histograms
 59      size_t idx = 0;
 60      const vector<string> dists{"pTlep1", "pTy", "pTjet1", "pTll", "pTlly", "mjj", "dRap", "dPhi", "cent"};
 61      for (const string& reg : vector<string>{"", "_ext"}) {
 62        for (const string& obs : dists) {
 63          if (reg == "" && (obs == "pTll" || obs == "cent"))  continue;
 64          book(_h[obs+reg], ++idx, 1, 1);
 65        }
 66      }
 67    }
 68
 69    /// Perform the per-event analysis
 70    void analyze(const Event &event) {
 71
 72      // Sort the dressed leptons by pt
 73      Particles electrons = apply<LeptonFinder>(event, "Electrons").particlesByPt();
 74      Particles muons = apply<LeptonFinder>(event, "Muons").particlesByPt();
 75
 76      if (electrons.size() < 2 && muons.size() < 2) vetoEvent;
 77
 78      Particles lep;
 79      if (electrons.size() >= 2) {
 80        lep.push_back(electrons[0]);
 81        lep.push_back(electrons[1]);
 82      }
 83      else if (muons.size() >= 2) {
 84        lep.push_back(muons[0]);
 85        lep.push_back(muons[1]);
 86      }
 87      if (lep.size() != 2)  vetoEvent;
 88      if (lep[0].pT() < 30*GeV) vetoEvent;
 89
 90      const double mll = (lep[0].momentum() + lep[1].momentum()).mass();
 91
 92      Particles photons = apply<PromptFinalState>(event, "Photons").particlesByPt();
 93      if (photons.empty())  vetoEvent;
 94
 95      idiscardIfAnyDeltaRLess(photons, electrons, 0.4);
 96      idiscardIfAnyDeltaRLess(photons, muons, 0.4);
 97
 98      Particles selectedPh;
 99      Particles fs = apply<VetoedFinalState>(event, "isolatedFS").particles();
100      for (const Particle& ph : photons) {
101        // check photon isolation
102        double coneEnergy(0.0);
103        for (const Particle& p : fs) {
104          if (deltaR(ph, p) < 0.2)  coneEnergy += p.Et(); // etcone20
105        }
106        if (coneEnergy / ph.pT() > 0.07)  continue;
107        selectedPh += ph;
108      }
109
110      if (selectedPh.size() < 1)  vetoEvent;
111
112      // Get jets
113      const Cut jetscut = (Cuts::pT > 25 * GeV && Cuts::absrap < 4.4);
114      Jets jetsMix = apply<FastJets>(event, "jets").jetsByPt(jetscut);
115
116      idiscardIfAnyDeltaRLess(jetsMix, photons, 0.4);
117      idiscardIfAnyDeltaRLess(jetsMix, electrons, 0.3);
118      idiscardIfAnyDeltaRLess(jetsMix, muons, 0.3);
119
120      // Jet multiplicity
121      const size_t njets = jetsMix.size();
122      if (njets < 2)  vetoEvent;
123
124      if (jetsMix[1].pT() < 50*GeV) vetoEvent; // sorted by pT, must all > 50GeV
125
126      if (mll < 40*GeV)   vetoEvent;
127
128      const FourMomentum lly_system = lep[0].momentum() + lep[1].momentum() + selectedPh[0].momentum();
129      const double mlly = lly_system.mass();
130      if (mll + mlly <= 182 * GeV) vetoEvent;
131
132      // other complex cuts
133      const FourMomentum jj_system = jetsMix[0].momentum() + jetsMix[1].momentum();
134      const double mjj = jj_system.mass();
135      if (mjj < 150*GeV)  vetoEvent;
136
137      const double Dyjj = fabs(deltaRap(jetsMix[0].momentum(), jetsMix[1].momentum()));
138      if (Dyjj < 1)  vetoEvent;
139
140      const double j1rap = jetsMix[0].rap();
141      const double j2rap = jetsMix[1].rap();
142
143      const double Centrality = fabs((lly_system.rap() - 0.5*(j1rap + j2rap)) / (j1rap - j2rap));
144      if (Centrality > 5) vetoEvent;
145
146      size_t njGap = 0;
147      for (size_t i=2; i < njets; ++i) {
148        const double jrap = jetsMix[i].rap();
149        if ((jrap < j1rap && jrap > j2rap) || (jrap < j2rap && jrap > j1rap))  ++njGap;
150      }
151      if (njGap > 0)  vetoEvent;
152
153      // now start to log hist
154      // first we select the extended SR
155      if (Centrality > 0.4) vetoEvent;
156
157      _h["pTjet1_ext"]->fill(jetsMix[0].pT() / GeV);
158      _h["pTy_ext"]->fill(selectedPh[0].pT() / GeV);
159      _h["pTlep1_ext"]->fill(lep[0].pT() / GeV);
160      _h["pTll_ext"]->fill((lep[0].momentum() + lep[1].momentum()).pT() / GeV);
161      _h["pTlly_ext"]->fill(lly_system.pT() / GeV);
162
163      _h["cent_ext"]->fill(Centrality);
164      _h["dRap_ext"]->fill(Dyjj);
165      _h["mjj_ext"]->fill(mjj);
166      _h["dPhi_ext"]->fill(fabs(deltaPhi(lly_system, jj_system)));
167
168      // deal with EWK SR500
169      if (mjj < 450 * GeV)  vetoEvent; // special since mjj in SR includes the underflow
170
171      _h["mjj"]->fill(mjj);
172
173      if (mjj < 500*GeV)  vetoEvent;
174
175      _h["pTjet1"]->fill(jetsMix[0].pT() / GeV);
176      _h["pTy"]->fill(selectedPh[0].pT() / GeV);
177      _h["pTlep1"]->fill(lep[0].pT() / GeV);
178      _h["pTlly"]->fill(lly_system.pT() / GeV);
179      _h["dRap"]->fill(Dyjj);
180      _h["dPhi"]->fill(fabs(deltaPhi(lly_system, jj_system)));
181    }
182
183    /// Normalise histograms etc., after the run
184    void finalize() {
185      scale(_h, crossSection() / femtobarn / sumOfWeights());
186    }
187
188  private:
189
190    /// Histograms
191    map<string, Histo1DPtr> _h;
192
193  };
194
195  RIVET_DECLARE_PLUGIN(ATLAS_2023_I2663725);
196
197}