rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_SUSY

Validate generic SUSY events, including various lepton invariant mass
Experiment: ()
Status: VALIDATED
Authors:
  • Andy Buckley
No references listed
Beams: * *
Beam energies: ANY
Run details:
  • SUSY events at any energy. $p_\perp$ cutoff at 10 GeV may be advised.

Analysis of generic SUSY events at the LHC, based on Atlas Herwig++ validation analysis contents. Plotted are $\eta$, $\phi$ and $p_\perp$ observables for charged tracks, photons, isolated photons, electrons, muons, and jets, as well as various dilepton mass `edge' plots for different event selection criteria.

Source code: MC_SUSY.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/IdentifiedFinalState.hh"
  7#include "Rivet/Projections/MissingMomentum.hh"
  8#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  9
 10namespace Rivet {
 11
 12  
 13
 14
 15  /// @brief MC validation analysis for SUSY events
 16  /// @author Andy Buckley
 17  class MC_SUSY : public Analysis {
 18  public:
 19
 20    /// Constructor
 21    MC_SUSY()
 22      : Analysis("MC_SUSY")
 23    {    }
 24
 25
 26    /// @name Analysis methods
 27    //@{
 28
 29    // Book histograms
 30    void init() {
 31      // Basic final state
 32      const FinalState fs((Cuts::etaIn(-4.0, 4.0) && Cuts::pT >=  10*GeV));
 33
 34      // Tracks and jets
 35      declare(ChargedFinalState(fs), "Tracks");
 36      declare(FastJets(fs, FastJets::ANTIKT, 0.7), "Jets");
 37
 38      IdentifiedFinalState photonfs(fs);
 39      photonfs.acceptId(PID::PHOTON);
 40      declare(photonfs, "AllPhotons");
 41
 42      IdentifiedFinalState efs(fs);
 43      efs.acceptIdPair(PID::ELECTRON);
 44      declare(efs, "Electrons");
 45
 46      IdentifiedFinalState mufs(fs);
 47      mufs.acceptIdPair(PID::MUON);
 48      declare(mufs, "Muons");
 49
 50      MissingMomentum missing(fs);
 51      declare(missing, "MET");
 52
 53      LeadingParticlesFinalState lpfs(fs);
 54      lpfs.addParticleIdPair(PID::ELECTRON);
 55      lpfs.addParticleIdPair(PID::MUON);
 56      declare(lpfs, "LeadingParticles");
 57
 58      book(_hist_n_trk   ,"n-trk", 50, 0.5, 300.5);
 59      book(_hist_phi_trk ,"phi-trk", 50, -PI, PI);
 60      book(_hist_eta_trk ,"eta-trk", 50, -4, 4);
 61      book(_hist_pt_trk  ,"pt-trk", 100, 0.0, 1500);
 62
 63      book(_hist_n_jet   ,"n-jet", 21, -0.5, 20.5);
 64      book(_hist_phi_jet ,"phi-jet", 50, -PI, PI);
 65      book(_hist_eta_jet ,"eta-jet", 50, -4, 4);
 66      book(_hist_pt_jet  ,"pt-jet", 100, 0.0, 1500);
 67
 68      book(_hist_n_e   ,"n-e", 11, -0.5, 10.5);
 69      book(_hist_phi_e ,"phi-e", 50, -PI, PI);
 70      book(_hist_eta_e ,"eta-e", 50, -4, 4);
 71      book(_hist_pt_e  ,"pt-e", 100, 0.0, 500);
 72
 73      book(_hist_n_mu   ,"n-mu", 11, -0.5, 10.5);
 74      book(_hist_phi_mu ,"phi-mu", 50, -PI, PI);
 75      book(_hist_eta_mu ,"eta-mu", 50, -4, 4);
 76      book(_hist_pt_mu  ,"pt-mu", 100, 0.0, 500);
 77
 78      book(_hist_n_gamma   ,"n-gamma", 11, -0.5, 10.5);
 79      book(_hist_phi_gamma ,"phi-gamma", 50, -PI, PI);
 80      book(_hist_eta_gamma ,"eta-gamma", 50, -4, 4);
 81      book(_hist_pt_gamma  ,"pt-gamma", 100, 0.0, 500);
 82
 83      book(_hist_n_gammaiso   ,"n-gamma-iso", 11, -0.5, 10.5);
 84      book(_hist_phi_gammaiso ,"phi-gamma-iso", 50, -PI, PI);
 85      book(_hist_eta_gammaiso ,"eta-gamma-iso", 50, -4, 4);
 86      book(_hist_pt_gammaiso  ,"pt-gamma-iso", 100, 0.0, 500);
 87
 88      book(_hist_met ,"Etmiss", 100, 0.0, 1500);
 89
 90      book(_hist_mll_ossf_ee   ,"mll-ossf-ee", 50, 0.0, 500);
 91      book(_hist_mll_ossf_mumu ,"mll-ossf-mumu", 50, 0.0, 500);
 92      book(_hist_mll_osof_emu  ,"mll-osof-emu", 50, 0.0, 500);
 93
 94      book(_hist_mll_all_ossf_ee   ,"mll-all-ossf-ee", 50, 0.0, 500);
 95      book(_hist_mll_all_ossf_mumu ,"mll-all-ossf-mumu", 50, 0.0, 500);
 96      book(_hist_mll_all_osof_emu  ,"mll-all-osof-emu", 50, 0.0, 500);
 97
 98      book(_hist_mll_2_ossf_ee   ,"mll-2-ossf-ee", 50, 0.0, 500);
 99      book(_hist_mll_2_ossf_mumu ,"mll-2-ossf-mumu", 50, 0.0, 500);
100      book(_hist_mll_2_osof_emu  ,"mll-2-osof-emu", 50, 0.0, 500);
101
102      /// @todo LSP eta, pT, phi, mass: no reliable cross-scenario LSP PID but
103      /// maybe plot for all of chi^0_1, gravitino, sneutrino, gluino, ... or
104      /// identify the LSP as any PID::isSUSY (?) particle with status = 1?
105    }
106
107
108    // Do the analysis
109    void analyze(const Event& evt) {
110      const FinalState& tracks = apply<FinalState>(evt, "Tracks");
111      if (tracks.particles().empty()) {
112        MSG_DEBUG("Failed multiplicity cut");
113        vetoEvent;
114      }
115
116      // Fill track histos
117      _hist_n_trk->fill(tracks.size());
118      for (const Particle& t : tracks.particles()) {
119        const FourMomentum& p = t.momentum();
120        _hist_phi_trk->fill(mapAngleMPiToPi(p.phi()));
121        _hist_eta_trk->fill(p.eta());
122        _hist_pt_trk->fill(p.pT()/GeV);
123      }
124
125      // Get jets and fill jet histos
126      const FastJets& jetpro = apply<FastJets>(evt, "Jets");
127      const Jets jets = jetpro.jetsByPt();
128      MSG_DEBUG("Jet multiplicity = " << jets.size());
129      _hist_n_jet->fill(jets.size());
130      for (const Jet& j : jets) {
131        const FourMomentum& pj = j.momentum();
132        _hist_phi_jet->fill(mapAngleMPiToPi(pj.phi()));
133        _hist_eta_jet->fill(pj.eta());
134        _hist_pt_jet->fill(pj.pT()/GeV);
135      }
136
137      /// @todo Resum photons around electrons
138
139      // Fill final state electron/positron histos
140      const FinalState& efs = apply<FinalState>(evt, "Electrons");
141      _hist_n_e->fill(efs.size());
142      vector<FourMomentum> epluses, eminuses;
143      for (const Particle& e : efs.particles()) {
144        const FourMomentum& p = e.momentum();
145        _hist_phi_e->fill(mapAngleMPiToPi(p.phi()));
146        _hist_eta_e->fill(p.eta());
147        _hist_pt_e->fill(p.pT()/GeV);
148        // Add sufficiently hard leptons to collections for m_ll histo
149        if (p.pT()/GeV > 20) {
150          if (PID::charge3(e.pid()) > 0) epluses += p; else eminuses += p;
151        }
152      }
153
154      /// @todo Resum photons around muons
155
156      // Fill final state muon/antimuon histos
157      const FinalState& mufs = apply<FinalState>(evt, "Muons");
158      _hist_n_mu->fill(mufs.size());
159      vector<FourMomentum> mupluses, muminuses;
160      for (const Particle& mu : mufs.particles()) {
161        const FourMomentum& p = mu.momentum();
162        _hist_phi_mu->fill(mapAngleMPiToPi(p.phi()));
163        _hist_eta_mu->fill(p.eta());
164        _hist_pt_mu->fill(p.pT()/GeV);
165        // Add sufficiently hard leptons to collections for m_ll histo
166        if (p.pT()/GeV > 20) {
167          if (PID::charge3(mu.pid()) > 0) mupluses += p; else muminuses += p;
168        }
169      }
170
171      // Fill final state non-isolated photon histos
172      const FinalState& allphotonfs = apply<FinalState>(evt, "AllPhotons");
173      _hist_n_gamma->fill(allphotonfs.size());
174      Particles isolatedphotons;
175      for (const Particle& ph : allphotonfs.particles()) {
176        const FourMomentum& p = ph.momentum();
177        _hist_phi_gamma->fill(mapAngleMPiToPi(p.phi()));
178        _hist_eta_gamma->fill(p.eta());
179        _hist_pt_gamma->fill(p.pT()/GeV);
180        // Select isolated photons
181        bool isolated = true;
182        for (const Jet& j : jets) {
183          if (deltaR(j.momentum(), p) < 0.2) {
184            isolated = false;
185            break;
186          }
187        }
188        if (isolated) isolatedphotons += ph;
189      }
190
191
192      // Fill final state isolated photon histos
193      _hist_n_gammaiso->fill(isolatedphotons.size());
194      for (const Particle& ph_iso : isolatedphotons) {
195        const FourMomentum& p = ph_iso.momentum();
196        _hist_phi_gammaiso->fill(mapAngleMPiToPi(p.phi()));
197        _hist_eta_gammaiso->fill(p.eta());
198        _hist_pt_gammaiso->fill(p.pT()/GeV);
199      }
200
201      // Calculate and fill missing Et histos
202      const MissingMomentum& met = apply<MissingMomentum>(evt, "MET");
203      _hist_met->fill(met.vectorEt().mod()/GeV);
204
205      // Choose highest-pT leptons of each sign and flavour for dilepton mass edges
206      const FinalState& lpfs = apply<FinalState>(evt, "LeadingParticles");
207      bool eplus_ok(false), eminus_ok(false), muplus_ok(false), muminus_ok(false);
208      FourMomentum peplus, peminus, pmuplus, pmuminus;
209      for (const Particle& p : lpfs.particles()) {
210        // Only use leptons above 20 GeV
211        if (p.pT()/GeV < 20) continue;
212        // Identify the PID
213        const PdgId pid = p.pid();
214        if (pid == PID::ELECTRON) {
215          eminus_ok = true;
216          peminus = p.momentum();
217        } else if (pid == PID::POSITRON) {
218          eplus_ok = true;
219          peplus = p.momentum();
220        } else if (pid == PID::MUON) {
221          muminus_ok = true;
222          pmuminus = p.momentum();
223        } else if (pid == PID::ANTIMUON) {
224          muplus_ok = true;
225          pmuplus = p.momentum();
226        } else {
227          throw Error("Unexpected particle type in leading particles FS!");
228        }
229      }
230      // m_ee
231      if (eminus_ok && eplus_ok) {
232        const double m_ee = FourMomentum(peplus + peminus).mass();
233        _hist_mll_ossf_ee->fill(m_ee/GeV);
234        if (epluses.size() == 1 && eminuses.size() == 1)
235          _hist_mll_2_ossf_ee->fill(m_ee/GeV);
236      }
237      // m_mumu
238      if (muminus_ok && muplus_ok) {
239        const double m_mumu = FourMomentum(pmuplus + pmuminus).mass();
240        _hist_mll_ossf_mumu->fill(m_mumu/GeV);
241        if (mupluses.size() == 1 && muminuses.size() == 1)
242          _hist_mll_2_ossf_mumu->fill(m_mumu/GeV);
243      }
244      // m_emu (both configurations)
245      if (eminus_ok && muplus_ok) {
246        const double m_emu = FourMomentum(pmuplus + peminus).mass();
247        _hist_mll_osof_emu->fill(m_emu/GeV);
248        if (mupluses.size() == 1 && eminuses.size() == 1)
249          _hist_mll_2_osof_emu->fill(m_emu/GeV);
250
251      }
252      if (muminus_ok && eplus_ok) {
253        const double m_mue = FourMomentum(peplus + pmuminus).mass();
254        _hist_mll_osof_emu->fill(m_mue/GeV);
255        if (epluses.size() == 1 && muminuses.size() == 1)
256          _hist_mll_2_osof_emu->fill(m_mue/GeV);
257      }
258
259
260      // m_ll plots using *all* electrons, positrons, muons and antimuons
261      // m_ee
262      for (const FourMomentum& peplus : epluses) {
263        for (const FourMomentum& peminus : eminuses) {
264          const double m_ee = FourMomentum(peplus + peminus).mass();
265          _hist_mll_all_ossf_ee->fill(m_ee/GeV);
266        }
267      }
268      // m_mumu
269      for (const FourMomentum& pmuplus : mupluses) {
270        for (const FourMomentum& pmuminus : muminuses) {
271          const double m_mumu = FourMomentum(pmuplus + pmuminus).mass();
272          _hist_mll_all_ossf_mumu->fill(m_mumu/GeV);
273        }
274      }
275      // m_emu (both configurations)
276      for (const FourMomentum& pmuplus : mupluses) {
277        for (const FourMomentum& peminus : eminuses) {
278          const double m_emu = FourMomentum(pmuplus + peminus).mass();
279          _hist_mll_all_osof_emu->fill(m_emu/GeV);
280        }
281      }
282      for (const FourMomentum& peplus : epluses) {
283        for (const FourMomentum& pmuminus : muminuses) {
284          const double m_mue = FourMomentum(peplus + pmuminus).mass();
285          _hist_mll_all_osof_emu->fill(m_mue/GeV);
286        }
287      }
288
289    }
290
291
292    void finalize() {
293      /// @todo Normalisations
294    }
295
296    //@}
297
298
299  private:
300
301    Histo1DPtr _hist_n_trk, _hist_phi_trk, _hist_eta_trk, _hist_pt_trk;
302    Histo1DPtr _hist_n_jet, _hist_phi_jet, _hist_eta_jet, _hist_pt_jet;
303    Histo1DPtr _hist_n_e, _hist_phi_e, _hist_eta_e, _hist_pt_e;
304    Histo1DPtr _hist_n_mu, _hist_phi_mu, _hist_eta_mu, _hist_pt_mu;
305    Histo1DPtr _hist_n_gamma, _hist_phi_gamma, _hist_eta_gamma, _hist_pt_gamma;
306    Histo1DPtr _hist_n_gammaiso, _hist_phi_gammaiso, _hist_eta_gammaiso, _hist_pt_gammaiso;
307    Histo1DPtr _hist_met;
308    Histo1DPtr _hist_mll_2_ossf_ee, _hist_mll_2_ossf_mumu, _hist_mll_2_osof_emu;
309    Histo1DPtr _hist_mll_ossf_ee, _hist_mll_ossf_mumu, _hist_mll_osof_emu;
310    Histo1DPtr _hist_mll_all_ossf_ee, _hist_mll_all_ossf_mumu, _hist_mll_all_osof_emu;
311  };
312
313
314
315  // The hook for the plugin system
316  RIVET_DECLARE_PLUGIN(MC_SUSY);
317
318}