Rivet analyses

title: MC_SUSY

Validate generic SUSY events, including various lepton invariant mass

Experiment: ()

Status: VALIDATED

Authors: - Andy Buckley

References: none listed

Beams: * *

Beam energies: ANY

Run details: - SUSY events at any energy. $\pT$ 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 $\pT$ 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

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/ChargedFinalState.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/IdentifiedFinalState.hh"

include "Rivet/Projections/LeadingParticlesFinalState.hh"

include "Rivet/Projections/MissingMomentum.hh"

namespace Rivet {

/// @brief MC validation analysis for SUSY events /// @author Andy Buckley class MC_SUSY : public Analysis { public:

/// Constructor
MC_SUSY()
    : Analysis("MC_SUSY") { }


/// @name Analysis methods
/// @{

// Book histograms
void init() {
  // Basic final state
  const FinalState fs((Cuts::etaIn(-4.0, 4.0) && Cuts::pT >= 10 * GeV));

  // Tracks and jets
  declare(ChargedFinalState(fs), "Tracks");
  declare(FastJets(fs, JetAlg::ANTIKT, 0.7), "Jets");

  IdentifiedFinalState photonfs(fs);
  photonfs.acceptId(PID::PHOTON);
  declare(photonfs, "AllPhotons");

  IdentifiedFinalState efs(fs);
  efs.acceptIdPair(PID::ELECTRON);
  declare(efs, "Electrons");

  IdentifiedFinalState mufs(fs);
  mufs.acceptIdPair(PID::MUON);
  declare(mufs, "Muons");

  MissingMomentum missing(fs);
  declare(missing, "MET");

  LeadingParticlesFinalState lpfs(fs);
  lpfs.addParticleIdPair(PID::ELECTRON);
  lpfs.addParticleIdPair(PID::MUON);
  declare(lpfs, "LeadingParticles");

  book(_hist_n_trk, "n-trk", 50, 0.5, 300.5);
  book(_hist_phi_trk, "phi-trk", 50, -PI, PI);
  book(_hist_eta_trk, "eta-trk", 50, -4, 4);
  book(_hist_pt_trk, "pt-trk", 100, 0.0, 1500);

  book(_hist_n_jet, "n-jet", 21, -0.5, 20.5);
  book(_hist_phi_jet, "phi-jet", 50, -PI, PI);
  book(_hist_eta_jet, "eta-jet", 50, -4, 4);
  book(_hist_pt_jet, "pt-jet", 100, 0.0, 1500);

  book(_hist_n_e, "n-e", 11, -0.5, 10.5);
  book(_hist_phi_e, "phi-e", 50, -PI, PI);
  book(_hist_eta_e, "eta-e", 50, -4, 4);
  book(_hist_pt_e, "pt-e", 100, 0.0, 500);

  book(_hist_n_mu, "n-mu", 11, -0.5, 10.5);
  book(_hist_phi_mu, "phi-mu", 50, -PI, PI);
  book(_hist_eta_mu, "eta-mu", 50, -4, 4);
  book(_hist_pt_mu, "pt-mu", 100, 0.0, 500);

  book(_hist_n_gamma, "n-gamma", 11, -0.5, 10.5);
  book(_hist_phi_gamma, "phi-gamma", 50, -PI, PI);
  book(_hist_eta_gamma, "eta-gamma", 50, -4, 4);
  book(_hist_pt_gamma, "pt-gamma", 100, 0.0, 500);

  book(_hist_n_gammaiso, "n-gamma-iso", 11, -0.5, 10.5);
  book(_hist_phi_gammaiso, "phi-gamma-iso", 50, -PI, PI);
  book(_hist_eta_gammaiso, "eta-gamma-iso", 50, -4, 4);
  book(_hist_pt_gammaiso, "pt-gamma-iso", 100, 0.0, 500);

  book(_hist_met, "Etmiss", 100, 0.0, 1500);

  book(_hist_mll_ossf_ee, "mll-ossf-ee", 50, 0.0, 500);
  book(_hist_mll_ossf_mumu, "mll-ossf-mumu", 50, 0.0, 500);
  book(_hist_mll_osof_emu, "mll-osof-emu", 50, 0.0, 500);

  book(_hist_mll_all_ossf_ee, "mll-all-ossf-ee", 50, 0.0, 500);
  book(_hist_mll_all_ossf_mumu, "mll-all-ossf-mumu", 50, 0.0, 500);
  book(_hist_mll_all_osof_emu, "mll-all-osof-emu", 50, 0.0, 500);

  book(_hist_mll_2_ossf_ee, "mll-2-ossf-ee", 50, 0.0, 500);
  book(_hist_mll_2_ossf_mumu, "mll-2-ossf-mumu", 50, 0.0, 500);
  book(_hist_mll_2_osof_emu, "mll-2-osof-emu", 50, 0.0, 500);

  /// @todo LSP eta, pT, phi, mass: no reliable cross-scenario LSP PID but
  /// maybe plot for all of chi^0_1, gravitino, sneutrino, gluino, ... or
  /// identify the LSP as any PID::isSUSY (?) particle with status = 1?
}


// Do the analysis
void analyze(const Event& evt) {
  const FinalState& tracks = apply<FinalState>(evt, "Tracks");
  if (tracks.particles().empty()) {
    MSG_DEBUG("Failed multiplicity cut");
    vetoEvent;
  }

  // Fill track histos
  _hist_n_trk->fill(tracks.size());
  for (const Particle& t : tracks.particles()) {
    const FourMomentum& p = t.momentum();
    _hist_phi_trk->fill(mapAngleMPiToPi(p.phi()));
    _hist_eta_trk->fill(p.eta());
    _hist_pt_trk->fill(p.pT() / GeV);
  }

  // Get jets and fill jet histos
  const FastJets& jetpro = apply<FastJets>(evt, "Jets");
  const Jets jets = jetpro.jetsByPt();
  MSG_DEBUG("Jet multiplicity = " << jets.size());
  _hist_n_jet->fill(jets.size());
  for (const Jet& j : jets) {
    const FourMomentum& pj = j.momentum();
    _hist_phi_jet->fill(mapAngleMPiToPi(pj.phi()));
    _hist_eta_jet->fill(pj.eta());
    _hist_pt_jet->fill(pj.pT() / GeV);
  }

  /// @todo Resum photons around electrons

  // Fill final state electron/positron histos
  const FinalState& efs = apply<FinalState>(evt, "Electrons");
  _hist_n_e->fill(efs.size());
  vector<FourMomentum> epluses, eminuses;
  for (const Particle& e : efs.particles()) {
    const FourMomentum& p = e.momentum();
    _hist_phi_e->fill(mapAngleMPiToPi(p.phi()));
    _hist_eta_e->fill(p.eta());
    _hist_pt_e->fill(p.pT() / GeV);
    // Add sufficiently hard leptons to collections for m_ll histo
    if (p.pT() / GeV > 20) {
      if (PID::charge3(e.pid()) > 0)
        epluses += p;
      else
        eminuses += p;
    }
  }

  /// @todo Resum photons around muons

  // Fill final state muon/antimuon histos
  const FinalState& mufs = apply<FinalState>(evt, "Muons");
  _hist_n_mu->fill(mufs.size());
  vector<FourMomentum> mupluses, muminuses;
  for (const Particle& mu : mufs.particles()) {
    const FourMomentum& p = mu.momentum();
    _hist_phi_mu->fill(mapAngleMPiToPi(p.phi()));
    _hist_eta_mu->fill(p.eta());
    _hist_pt_mu->fill(p.pT() / GeV);
    // Add sufficiently hard leptons to collections for m_ll histo
    if (p.pT() / GeV > 20) {
      if (PID::charge3(mu.pid()) > 0)
        mupluses += p;
      else
        muminuses += p;
    }
  }

  // Fill final state non-isolated photon histos
  const FinalState& allphotonfs = apply<FinalState>(evt, "AllPhotons");
  _hist_n_gamma->fill(allphotonfs.size());
  Particles isolatedphotons;
  for (const Particle& ph : allphotonfs.particles()) {
    const FourMomentum& p = ph.momentum();
    _hist_phi_gamma->fill(mapAngleMPiToPi(p.phi()));
    _hist_eta_gamma->fill(p.eta());
    _hist_pt_gamma->fill(p.pT() / GeV);
    // Select isolated photons
    bool isolated = true;
    for (const Jet& j : jets) {
      if (deltaR(j.momentum(), p) < 0.2) {
        isolated = false;
        break;
      }
    }
    if (isolated) isolatedphotons += ph;
  }


  // Fill final state isolated photon histos
  _hist_n_gammaiso->fill(isolatedphotons.size());
  for (const Particle& ph_iso : isolatedphotons) {
    const FourMomentum& p = ph_iso.momentum();
    _hist_phi_gammaiso->fill(mapAngleMPiToPi(p.phi()));
    _hist_eta_gammaiso->fill(p.eta());
    _hist_pt_gammaiso->fill(p.pT() / GeV);
  }

  // Calculate and fill missing Et histos
  const MissingMomentum& met = apply<MissingMomentum>(evt, "MET");
  _hist_met->fill(met.vectorEt().mod() / GeV);

  // Choose highest-pT leptons of each sign and flavour for dilepton mass edges
  const FinalState& lpfs = apply<FinalState>(evt, "LeadingParticles");
  bool eplus_ok(false), eminus_ok(false), muplus_ok(false), muminus_ok(false);
  FourMomentum peplus, peminus, pmuplus, pmuminus;
  for (const Particle& p : lpfs.particles()) {
    // Only use leptons above 20 GeV
    if (p.pT() / GeV < 20) continue;
    // Identify the PID
    const PdgId pid = p.pid();
    if (pid == PID::ELECTRON) {
      eminus_ok = true;
      peminus = p.momentum();
    }
    else if (pid == PID::POSITRON) {
      eplus_ok = true;
      peplus = p.momentum();
    }
    else if (pid == PID::MUON) {
      muminus_ok = true;
      pmuminus = p.momentum();
    }
    else if (pid == PID::ANTIMUON) {
      muplus_ok = true;
      pmuplus = p.momentum();
    }
    else {
      throw Error("Unexpected particle type in leading particles FS!");
    }
  }
  // m_ee
  if (eminus_ok && eplus_ok) {
    const double m_ee = FourMomentum(peplus + peminus).mass();
    _hist_mll_ossf_ee->fill(m_ee / GeV);
    if (epluses.size() == 1 && eminuses.size() == 1) _hist_mll_2_ossf_ee->fill(m_ee / GeV);
  }
  // m_mumu
  if (muminus_ok && muplus_ok) {
    const double m_mumu = FourMomentum(pmuplus + pmuminus).mass();
    _hist_mll_ossf_mumu->fill(m_mumu / GeV);
    if (mupluses.size() == 1 && muminuses.size() == 1) _hist_mll_2_ossf_mumu->fill(m_mumu / GeV);
  }
  // m_emu (both configurations)
  if (eminus_ok && muplus_ok) {
    const double m_emu = FourMomentum(pmuplus + peminus).mass();
    _hist_mll_osof_emu->fill(m_emu / GeV);
    if (mupluses.size() == 1 && eminuses.size() == 1) _hist_mll_2_osof_emu->fill(m_emu / GeV);
  }
  if (muminus_ok && eplus_ok) {
    const double m_mue = FourMomentum(peplus + pmuminus).mass();
    _hist_mll_osof_emu->fill(m_mue / GeV);
    if (epluses.size() == 1 && muminuses.size() == 1) _hist_mll_2_osof_emu->fill(m_mue / GeV);
  }


  // m_ll plots using *all* electrons, positrons, muons and antimuons
  // m_ee
  for (const FourMomentum& peplus : epluses) {
    for (const FourMomentum& peminus : eminuses) {
      const double m_ee = FourMomentum(peplus + peminus).mass();
      _hist_mll_all_ossf_ee->fill(m_ee / GeV);
    }
  }
  // m_mumu
  for (const FourMomentum& pmuplus : mupluses) {
    for (const FourMomentum& pmuminus : muminuses) {
      const double m_mumu = FourMomentum(pmuplus + pmuminus).mass();
      _hist_mll_all_ossf_mumu->fill(m_mumu / GeV);
    }
  }
  // m_emu (both configurations)
  for (const FourMomentum& pmuplus : mupluses) {
    for (const FourMomentum& peminus : eminuses) {
      const double m_emu = FourMomentum(pmuplus + peminus).mass();
      _hist_mll_all_osof_emu->fill(m_emu / GeV);
    }
  }
  for (const FourMomentum& peplus : epluses) {
    for (const FourMomentum& pmuminus : muminuses) {
      const double m_mue = FourMomentum(peplus + pmuminus).mass();
      _hist_mll_all_osof_emu->fill(m_mue / GeV);
    }
  }
}


void finalize() {
  /// @todo Normalisations
}

/// @}

private:

Histo1DPtr _hist_n_trk, _hist_phi_trk, _hist_eta_trk, _hist_pt_trk;
Histo1DPtr _hist_n_jet, _hist_phi_jet, _hist_eta_jet, _hist_pt_jet;
Histo1DPtr _hist_n_e, _hist_phi_e, _hist_eta_e, _hist_pt_e;
Histo1DPtr _hist_n_mu, _hist_phi_mu, _hist_eta_mu, _hist_pt_mu;
Histo1DPtr _hist_n_gamma, _hist_phi_gamma, _hist_eta_gamma, _hist_pt_gamma;
Histo1DPtr _hist_n_gammaiso, _hist_phi_gammaiso, _hist_eta_gammaiso, _hist_pt_gammaiso;
Histo1DPtr _hist_met;
Histo1DPtr _hist_mll_2_ossf_ee, _hist_mll_2_ossf_mumu, _hist_mll_2_osof_emu;
Histo1DPtr _hist_mll_ossf_ee, _hist_mll_ossf_mumu, _hist_mll_osof_emu;
Histo1DPtr _hist_mll_all_ossf_ee, _hist_mll_all_ossf_mumu, _hist_mll_all_osof_emu;

};

RIVET_DECLARE_PLUGIN(MC_SUSY);

} ```