1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/LeptonFinder.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief ZZ production at 13.6 TeV
 13  class ATLAS_2023_I2723369 : public Analysis {
 14  public:
 15
 16    /// Default constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2023_I2723369);
 18
 19    void init() {
 20
 21      /// Dressed leptons
 22      Cut cut_lep = (Cuts::abspid == PID::MUON && Cuts::abseta < 2.5 && Cuts::pT > 5*GeV) ||
 23                    (Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && Cuts::pT > 7*GeV);
 24	    PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 25      PromptFinalState leptons(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
 26	    LeptonFinder dLeptons(leptons, photons, 0.1, cut_lep);
 27      declare(dLeptons, "AllLeptons");
 28
 29      /// Jet inputs
 30      FinalState fs_jet(Cuts::abseta < 5.0);
 31      VetoedFinalState jet_input(fs_jet);
 32
 33      // reject all leptons dressed with only prompt photons from jet input
 34      FinalState all_leptons(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
 35   		LeptonFinder reject_leptons(all_leptons, photons, 0.1);
 36      jet_input.addVetoOnThisFinalState(reject_leptons);
 37
 38      // reject prompt invisibles, including from tau decays
 39      VetoedFinalState invis_fs_jet(fs_jet);
 40      invis_fs_jet.addVetoOnThisFinalState(VisibleFinalState(fs_jet));
 41      PromptFinalState invis_pfs_jet(invis_fs_jet, TauDecaysAs::PROMPT);
 42      jet_input.addVetoOnThisFinalState(invis_pfs_jet);
 43
 44      // declare jets
 45      FastJets jets(jet_input, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::DECAY);
 46      declare(jets, "Jets");
 47
 48      // Book histograms
 49      book(_h["m4l"],  1, 1, 1);
 50      book(_h["pt4l"], 2, 1, 1);
 51
 52    }
 53
 54    /// Do the per-event analysis
 55    void analyze(const Event& e) {
 56
 57
 58      const DressedLeptons& all_leps = apply<LeptonFinder>(e, "AllLeptons").dressedLeptons();
 59      size_t n_parts = all_leps.size();
 60      size_t n_OSSF_pairs = 0;
 61
 62      // Form Z candidate (opposite-sign same-flavour) lepton pairs
 63      std::vector<Zstate> dileptons;
 64      for (size_t i = 0; i < n_parts; ++i) {
 65        for (size_t j = i + 1; j < n_parts; ++j) {
 66          if (isOSSF(all_leps[i], all_leps[j])){
 67            n_OSSF_pairs += 1;
 68            // Set positive charge first for later ME calculation
 69            if (all_leps[i].charge() > 0) {
 70              dileptons.emplace_back(make_pair(all_leps[i], all_leps[j]));
 71            } else {
 72              dileptons.emplace_back(make_pair(all_leps[j], all_leps[i]));
 73            }
 74          }
 75        }
 76      }
 77
 78      // At least two pairs required to select ZZ->llll final state
 79      if (n_OSSF_pairs < 2) vetoEvent;
 80
 81      // Form the quadruplet of two lepon pairs passing kinematics cuts
 82      std::vector<Quadruplet> quadruplets;
 83      for (size_t i = 0; i < dileptons.size(); ++i) {
 84        for (size_t j = i+1; j < dileptons.size(); ++j) {
 85          // Only use unique leptons
 86          if (isSame( dileptons[i].first , dileptons[j].first  )) continue;
 87          if (isSame( dileptons[i].first , dileptons[j].second )) continue;
 88          if (isSame( dileptons[i].second, dileptons[j].first  )) continue;
 89          if (isSame( dileptons[i].second, dileptons[j].second )) continue;
 90
 91          Particles leptons{ dileptons[i].first, dileptons[i].second,
 92                             dileptons[j].first, dileptons[j].second };
 93          isortByPt(leptons);
 94
 95          // Apply kinematic cuts
 96          if (leptons[0].pT() < 20*GeV)  continue;
 97          if (leptons[1].pT() < 10*GeV)  continue;
 98
 99          // Form the quad with pair closest to Z pole first
100          if (dileptons[i].Zdist() < dileptons[j].Zdist()) {
101            quadruplets.emplace_back(dileptons[i], dileptons[j]);
102          } else {
103            quadruplets.emplace_back(dileptons[j], dileptons[i]);
104          }
105        }
106      }
107
108      // Veto if no quad passes kinematic selection
109      size_t n_quads = quadruplets.size();
110      if (n_quads == 0) vetoEvent;
111
112      // To resolve ambiguities in lepton pairing order quads by channel priority first, then m12 - mz and m34 - mz
113      // The first in every channel is considered nominal
114      std::sort(quadruplets.begin(), quadruplets.end(),
115        [](const Quadruplet& q1, const Quadruplet& q2) {
116          if (q1.ch_priority == q2.ch_priority) {
117            // if rarely, Z1 the same distance from the Z pole, compare Z2
118            if (fabs( q1.Z1().Zdist() - q2.Z1().Zdist() ) < 1.e-5) {
119              return q1.Z2().Zdist() < q2.Z2().Zdist();
120            }
121            return q1.Z1().Zdist() < q2.Z1().Zdist();
122          }
123          return q1.ch_priority < q2.ch_priority;
124      });
125
126
127      // Select the best quad
128      Particles leptons_sel4l;
129      Quadruplet quadSel;
130      bool atleastonequadpassed = false;
131
132      for (size_t iquad = 0; iquad < n_quads; ++iquad) {
133
134        // Veto event if nominal quad was not selected in 4 lepton case
135        if (n_parts == 4 && iquad > 0) vetoEvent;
136
137        const Quadruplet& quad = quadruplets[iquad];
138
139        // Z invariant mass requirements
140        if (!(inRange(quad.Z1().mom().mass(), 66*GeV, 116*GeV))) continue;
141        if (!(inRange(quad.Z2().mom().mass(), 66*GeV, 116*GeV))) continue;
142
143        // Lepton separation and J/Psi veto
144        bool b_pass_leptonseparation = true;
145        bool b_pass_jpsi = true;
146        leptons_sel4l.clear();
147        leptons_sel4l.push_back(quad.Z1().first);
148        leptons_sel4l.push_back(quad.Z1().second);
149        leptons_sel4l.push_back(quad.Z2().first);
150        leptons_sel4l.push_back(quad.Z2().second);
151
152        for (size_t i = 0; i < 4; ++i) {
153          for (size_t j = i+1; j < 4; ++j) {
154            if ( deltaR( leptons_sel4l[i], leptons_sel4l[j]) < 0.1) b_pass_leptonseparation = false;
155            if ( isOSSF(leptons_sel4l[i], leptons_sel4l[j]) && \
156                 (leptons_sel4l[i].mom() + leptons_sel4l[j].mom()).mass() <= 5.*GeV) b_pass_jpsi = false;
157          }
158        }
159        if (b_pass_leptonseparation == false || b_pass_jpsi == false) continue;
160
161        isortByPt(leptons_sel4l);
162        atleastonequadpassed = true;
163        quadSel = quad;
164        break;
165
166      }
167
168      if (!atleastonequadpassed)  vetoEvent;
169
170      // Veto if quad not in Higgs mass window
171      const FourMomentum ZZ = quadSel.mom();
172      _h["pt4l"]->fill(ZZ.pT()/GeV);
173      _h["m4l"]->fill(ZZ.mass()/GeV);
174
175    }
176
177
178    void finalize() {
179      scale(_h, crossSection() / femtobarn  / sumOfWeights());
180    }
181
182  private:
183
184    map<string, Histo1DPtr> _h;
185
186    /// Generic Z candidate
187    struct Zstate : public ParticlePair {
188
189      Zstate() { }
190
191      Zstate(ParticlePair&& _particlepair) : ParticlePair(std::move(_particlepair)) { }
192
193      FourMomentum mom() const { return first.momentum() + second.momentum(); }
194
195      double Zdist() const { return fabs(mom().mass() -  91.1876*GeV); }
196
197      int flavour() const { return first.abspid(); }
198    };
199
200    /// Generic quadruplet
201    struct Quadruplet {
202
203      // find out which type it is: 4mu = 0, 4e = 1, 2mu2e = 2, 2e2mu = 3 (mm, ee, me, em)
204      // channel priority is 4m, 2e2m, 2m2e, 4e
205      enum class FlavCombi { mm=0, ee, me, em, undefined };
206
207      Zstate _z1, _z2;
208
209      FlavCombi _type;
210
211      int ch_priority;
212
213      Quadruplet() { }
214
215      Quadruplet(const Zstate& z1, const Zstate& z2) : _z1(z1), _z2(z2) {
216        if (     _z1.flavour() == 13 && _z2.flavour() == 13) { _type = FlavCombi::mm; ch_priority = 0; }
217        else if (_z1.flavour() == 11 && _z2.flavour() == 11) { _type = FlavCombi::ee; ch_priority = 3; }
218        else if (_z1.flavour() == 13 && _z2.flavour() == 11) { _type = FlavCombi::me; ch_priority = 2; }
219        else if (_z1.flavour() == 11 && _z2.flavour() == 13) { _type = FlavCombi::em; ch_priority = 1; }
220        else  { _type = FlavCombi::undefined; }
221      }
222
223      const Zstate& Z1() const { return _z1; }
224
225      const Zstate& Z2() const { return _z2; }
226
227      FourMomentum mom() const { return _z1.mom() + _z2.mom(); }
228
229      const FlavCombi& type() const { return _type; }
230    };
231
232
233  };
234
235  RIVET_DECLARE_PLUGIN(ATLAS_2023_I2723369);
236}