1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/LeptonFinder.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9#include "Rivet/Projections/InvisibleFinalState.hh"
 10#include "Rivet/Tools/RivetMT2.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// @brief WW production in pp at 13 TeV in electroweak SUSY-inspired phase space
 16  class ATLAS_2022_I2103950 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2022_I2103950);
 21
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    /// Book histograms and initialise projections before the run
 27    void init() {
 28
 29      const FinalState fs(Cuts::abseta < 5);
 30
 31      // Project photons for dressing
 32      FinalState photons(Cuts::abspid == PID::PHOTON);
 33
 34      // Cuts for leptons
 35      Cut lepton_cuts = ((Cuts::abseta < 2.6 && Cuts::abspid == PID::MUON ) ||
 36                         (Cuts::abseta < 2.47 && Cuts::abspid == PID::ELECTRON ) );
 37
 38      // Muons
 39      PromptFinalState bare_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 40      LeptonFinder all_dressed_mu(bare_mu, photons, 0.1);
 41
 42      // Electrons
 43      PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 44      LeptonFinder all_dressed_el(bare_el, photons, 0.1);
 45
 46      //Jet forming
 47      VetoedFinalState vfs(FinalState(Cuts::abseta < 5));
 48
 49      InvisibleFinalState prompt_invis(OnlyPrompt::YES, TauDecaysAs::PROMPT);
 50      vfs.addVetoOnThisFinalState(prompt_invis);
 51      vfs.addVetoOnThisFinalState(all_dressed_el);
 52      vfs.addVetoOnThisFinalState(all_dressed_mu);
 53
 54      FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 55      declare(jets, "jets");
 56
 57
 58      // Project dressed leptons (e/mu not from tau) with pT > 25 GeV
 59      PromptFinalState lep_bare(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 60      declare(lep_bare,"lep_bare");
 61      PromptFinalState prompt_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 62      PromptFinalState prompt_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 63
 64      LeptonFinder lep_dressed(lep_bare, photons, 0.1, lepton_cuts);
 65      declare(lep_dressed,"lep_dressed");
 66      LeptonFinder elecs(prompt_el, photons, 0.1, lepton_cuts);
 67      declare(elecs, "elecs");
 68      LeptonFinder muons(prompt_mu, photons, 0.1, lepton_cuts);
 69      declare(muons, "muons");
 70
 71      // Get MET from generic invisibles
 72      VetoedFinalState ivfs(fs);
 73      ivfs.addVetoOnThisFinalState(VisibleFinalState(fs));
 74      declare(ivfs, "InvisibleFS");
 75
 76      // fiducial differential cross sections
 77      book(_h["ptlead"], 13, 1, 1);
 78      book(_h["mll"], 15, 1, 1);
 79      book(_h["ptll"], 17, 1, 1);
 80      book(_h["yll"], 7, 1, 1);
 81      book(_h["dphill"], 9, 1, 1);
 82      book(_h["costhetastarll"], 11, 1, 1);
 83
 84    }
 85
 86
 87    /// Perform the per-event analysis
 88    void analyze(const Event& event) {
 89
 90
 91      const FinalState& ifs = apply<VetoedFinalState>(event, "InvisibleFS");
 92      const FourMomentum metinvisible = sum(ifs.particles(), FourMomentum());
 93
 94      // Get met and find leptons
 95      //const DressedLeptons &leptons = apply<LeptonFinder>(event, "lep_dressed").dressedLeptons();
 96      const DressedLeptons& all_elecs = apply<LeptonFinder>(event, "elecs").dressedLeptons();
 97      const DressedLeptons& all_muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
 98      //Particles bare_leps  = apply<IdentifiedFinalState>(event, "bare_leptons").particles();
 99
100      // Find jets and jets for simplified phase space (for the latter slightly different leptons are excluded from clustering)
101      Jets alljets = apply<FastJets>(event, "jets").jetsByPt(Cuts::abseta < 4.5 && Cuts::pT > 20*GeV);
102
103      DressedLeptons elecs_muonOR;
104      for (const DressedLepton& el : all_elecs) {
105        bool overlaps = false;
106        for (const DressedLepton& mu : all_muons) {
107          if (deltaR(el, mu) < 0.01) {
108            overlaps = true;
109            break;
110          }
111        }
112        if (overlaps) continue;
113        elecs_muonOR.push_back(el);
114      }
115
116      // jet electron overlap removal
117      for (const DressedLepton& e : elecs_muonOR) {
118        idiscard(alljets, deltaRLess(e, 0.2, RAPIDITY));
119      }
120
121      // muon jet overlap removal
122      DressedLeptons muons;
123      for (const DressedLepton& mu : all_muons) {
124        float dRcut=0.4;
125        if ((0.04+10.0/mu.pT()) < 0.4){
126          dRcut = 0.04+10.0/mu.pT();
127        }
128        bool overlaps = false;
129        for (const Jet& jet : alljets) {
130          if (deltaR(mu, jet) < dRcut) {
131            overlaps = true;
132            break;
133          }
134        }
135        if (overlaps) continue;
136        muons.push_back(mu);
137      }
138      // electron jet overlap removal
139      DressedLeptons elecs;
140      for (const DressedLepton& el : elecs_muonOR) {
141        float dRcut=0.4;
142        if ((0.04+10/el.pT()) < 0.4){
143          dRcut = 0.04+10/el.pT();
144        }
145        bool overlaps = false;
146        for (const Jet& jet : alljets) {
147          if (deltaR(el, jet) < dRcut) {
148            overlaps = true;
149            break;
150          }
151        }
152        if (overlaps) continue;
153        elecs.push_back(el);
154      }
155      Jets jets20;
156      for (const Jet& jet : alljets) {
157      	if (jet.abseta()<2.4 ) jets20 += jet;
158      }
159      // Remove events that do not contain 2 good leptons (either muons or electrons)
160      if ((elecs.size()+muons.size()) !=2) vetoEvent;
161      // only select electron-muon events
162      if (elecs.size() != 1) vetoEvent;
163      DressedLeptons leptons;
164      if (elecs[0].pT()>muons[0].pT()) {
165         leptons.push_back(elecs[0]);
166         leptons.push_back(muons[0]);
167      }
168      else {
169         leptons.push_back(muons[0]);
170         leptons.push_back(elecs[0]);
171      }
172
173
174      // Define observables
175      const FourMomentum dilep  = leptons.size()>1 ? leptons[0].momentum() + leptons[1].momentum() : FourMomentum(0,0,0,0);
176      const double ptll         = leptons.size()>1 ? dilep.pT()/GeV : -1;
177      const double Mll          = leptons.size()>1 ? dilep.mass()/GeV : -1;
178      const double Yll          = leptons.size()>1 ? dilep.absrap() : -5;
179      const double DPhill       = leptons.size()>1 ? fabs(deltaPhi(leptons[0], leptons[1])) : -1.;
180      const double costhetastar = leptons.size()>1 ? fabs(tanh((leptons[0].eta() - leptons[1].eta()) / 2)) : -0.2;
181
182      const double mt2 = leptons.size()>1 ? mT2(leptons[0],leptons[1], metinvisible,0) : 0;
183
184      // Event selection for proper fiducial phase space
185      if ( leptons.size() != 2)  vetoEvent;
186      if ( leptons[0].pT()< 25*GeV || leptons[1].pT()< 25*GeV )  vetoEvent;
187
188      // Veto same-flavour events
189      if ( leptons[0].abspid() == leptons[1].abspid())  vetoEvent;
190
191      // Veto same-charge events
192      if ( leptons[0].pid()*leptons[1].pid()>0)  vetoEvent;
193
194      // jet veto
195      if (!(jets20.size()==0))  vetoEvent;
196
197      if (metinvisible.pT() <= 60*GeV || metinvisible.pT() > 80*GeV )  vetoEvent;
198
199      // m_ll cut
200      if (dilep.mass() <= 100*GeV)  vetoEvent;
201
202      //mt2 cut
203      if( (mt2<=60*GeV) || (mt2>80*GeV)) vetoEvent;
204
205      // Jet veto at 35 GeV is the default
206      if (!jets20.empty())  vetoEvent;
207
208      // fill histograms
209
210      _h["ptlead"]->fill(leptons[0].pT()/GeV);
211      _h["ptll"]->fill(ptll);
212      _h["mll"]->fill(Mll);
213      _h["yll"]->fill(Yll);
214      _h["dphill"]->fill(DPhill);
215      _h["costhetastarll"]->fill(costhetastar);
216
217    }
218
219
220    /// Normalise histograms etc., after the run
221    void finalize() {
222      const double sf(crossSection()/femtobarn/sumOfWeights());
223      // scale to cross section
224      scale(_h, sf);
225    }
226
227    /// @}
228
229  private:
230
231    /// Histograms
232    map<string, Histo1DPtr> _h;
233
234  };
235
236
237  RIVET_DECLARE_PLUGIN(ATLAS_2022_I2103950);
238
239}