Rivet analyses

title: ATLAS_2019_I1725190

Dilepton mass spectrum in 13 TeV pp collisions with 139/fb Run 2 dataset

Experiment: ATLAS (LHC)

Inspire ID: 1725190

Status: VALIDATED SINGLEWEIGHT NOREENTRY

Authors: - Andy Buckley

References: - Expt page: ATLAS-EXOT-2018-08 - Phys.Lett. B796 (2019) 68-87 - DOI:10.1016/j.physletb.2019.07.016 - arXiv: 1903.06248

Beams: p+ p+

Beam energies: (6500.0, 6500.0)GeV

Run details: - Dilepton events, either full process, BSM signal-only, or signal + SM interference.

A search for high-mass dielectron and dimuon resonances in the mass range of 250 GeV to 6 TeV. Functional forms have been fitted to the components of Crystal Ball + Gaussian dilepton invariant-mass distribution smearing functions for electrons and muons separately, which are encoded in this analysis to allow particle-level MC to be compared to the reco-level experimental data.

Source code:ATLAS_2019_I1725190.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/PromptFinalState.hh"

include "Rivet/Projections/Smearing.hh"

include "Rivet/Tools/Random.hh"

namespace Rivet {

/// @brief Dilepton high-mass resonance search /// /// @todo Use the proper smearing system rather than hand-rolled sampling class ATLAS_2019_I1725190 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1725190);


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

/// Book histograms and initialise projections before the run
void init() {

  PromptFinalState electrons(Cuts::abspid == PID::ELECTRON && Cuts::Et > 30 * GeV
                             && (Cuts::abseta < 2.47 && !Cuts::absetaIn(1.37, 1.52)));
  SmearedParticles recoelectrons(electrons, PARTICLE_EFF_ONE); //ELECTRON_EFF_ATLAS_RUN2_MEDIUM);
  declare(recoelectrons, "Elecs");

  PromptFinalState muons(Cuts::abspid == PID::MUON && Cuts::pT > 30 * GeV && Cuts::abseta < 2.5);
  SmearedParticles recomuons(muons, PARTICLE_EFF_ONE);
  // [](const Particle& m) -> double {
  //   if (m.pT() < 1*TeV) return 0.69;
  //   if (m.pT() > 2.5*TeV) return 0.57;
  //   return 0.69 - 0.12*(m.pT() - 1*TeV)/(2.5*TeV - 1*TeV);
  // });
  declare(recomuons, "Muons");

  // Book histograms
  book(_h_mee, 1, 1, 1);
  book(_h_mmm, 2, 1, 1);
}


/// Perform the per-event analysis
void analyze(const Event& event) {

  // Get leptons
  const Particles elecs = apply<ParticleFinder>(event, "Elecs").particlesByPt();
  const Particles muons = apply<ParticleFinder>(event, "Muons").particlesByPt();

  // MSG_INFO("Num e, mu = " << elecs.size() << ", " << muons.size());
  // for (const Particle& e : elecs) MSG_INFO(e.mom());
  // for (const Particle& m : muons) MSG_INFO(m.mom());

  // Isolation
  /// @todo Can't be done from provided information?
  // Particles isoelecs, isomuons;
  // for (const Particle& e : elecs) {
  // }

  // Choose the highest-pT lepton pair, preferring electrons
  if (elecs.size() < 2 && muons.size() < 2) vetoEvent;
  const Particle l1 = (elecs.size() >= 2) ? elecs[0] : muons[0];
  const Particle l2 = (elecs.size() >= 2) ? elecs[1] : muons[1];

  // Require opposite sign for muons only
  const bool mumu = (l1.abspid() == PID::MUON);
  if (mumu && l1.pid() * l2.pid() > 0) vetoEvent;

  // Get the true dilepton pair
  const FourMomentum pll = l1.mom() + l2.mom();
  const double mll = pll.mass();

  // Make sure we're in a region where the smearing and efficiencies are well-behaved
  if (mll < 200 * GeV) vetoEvent;

  // Apply dilepton efficiency curves (as function of mX ~ mll)
  const double eff = mumu
      ? (0.54 - (mll - 200 * GeV) / (6000 * GeV - 200 * GeV) * (0.54 - 0.38))
      : (0.74 - 0.04 * exp(-(mll - 200 * GeV) / 100 * GeV) - 0.08 * exp(-(mll - 200 * GeV) / 1000 * GeV));
  if (rand01() > eff) vetoEvent;

  // Smear the dilepton mass with a CB + Gauss function
  double muCB, sigCB, alpCB, nCB, muG, sigG, kappa;
  if (!mumu) {
    const double lnm = log(mll);
    static const vector<double> pmuCB = {0.13287, -0.410663, -0.0126743, 2.9547e-6};
    muCB = pmuCB[0] + pmuCB[1] / lnm + pmuCB[2] * lnm + pmuCB[3] * pow(lnm, 4);
    static const vector<double> psigCB = {0.0136624, 0.230678, 1.73254};
    sigCB = sqrt(pow(psigCB[0], 2) + pow(psigCB[1], 2) / mll + pow(psigCB[2] / mll, 2));
    alpCB = 1.59112;
    static const vector<double> pnCB = {1.13055, 0.76705, 0.00298312};
    nCB = pnCB[0] + pnCB[1] * exp(-pnCB[2] * mll);
    static const vector<double> pmuG = {-0.00402708, 0.814172, -3.94281e-7,
                                        7.97076e-6,  -87.6397, -1.64806e-11};
    muG = pmuG[0] + pmuG[1] / mll + pmuG[2] * mll + pmuG[3] * pow(lnm, 3) + pmuG[4] / sqr(mll)
        + pmuG[5] * sqr(mll);
    static const vector<double> psigG = {0.00553858, 0.140909, 0.644418};
    sigG = sqrt(sqr(psigG[0]) + sqr(psigG[1]) / mll + sqr(psigG[2] / mll));
    static const vector<double> pkappa = {0.347003, 0.135768, 0.00372497, -2.2822e-5, 5.06351e-13};
    kappa = pkappa[0] + pkappa[1] * exp(-pkappa[2] * mll) + pkappa[3] * mll + pkappa[4] * pow(mll, 3);
  }
  else {
    static const vector<double> pmuCB = {-0.0891397, 10.6169,     -951.712,    74775.3,
                                         5.60192e-5, -1.58827e-9, -3.81706e-13};
    muCB = pmuCB[0] + pmuCB[1] / mll + pmuCB[2] / sqr(mll) + pmuCB[3] / pow(mll, 3) + pmuCB[4] * mll
        + pmuCB[5] * sqr(mll) + pmuCB[6] * pow(mll, 3);
    static const vector<double> psigCB = {0.0836349, -8.98476, 491.19, 5.18068e-5, -3.45042e-10};
    sigCB = psigCB[0] + psigCB[1] / mll + psigCB[2] / sqr(mll) + psigCB[3] * mll + psigCB[4] * sqr(mll);
    static const vector<double> palpCB = {0.512577, 252.922, -79337.4, 7.31863e6, 0.000237883};
    alpCB = palpCB[0] + palpCB[1] / mll + palpCB[2] / sqr(mll) + palpCB[3] / pow(mll, 3)
        + palpCB[4] * mll;
    static const vector<double> pnCB = {1.13055, 0.76705, 0.00298312};
    nCB = 6.08818;
    static const vector<double> pmuG = {-0.00410659, 2.82352e-6, -6.264e-9, 1.25547e-12, -9.94431e-17};
    muG = pmuG[0] + pmuG[1] * mll + pmuG[2] * sqr(mll) + pmuG[3] * pow(mll, 3) + pmuG[4] * pow(mll, 4);
    static const vector<double> psigG = {0.0214264, -0.795058, 15.4726, 3.38205e-5, -1.64984e-9};
    sigG = psigG[0] + psigG[1] / mll + psigG[2] / sqr(mll) + psigG[3] * mll + psigG[4] * sqr(mll);
    static const vector<double> pkappa = {0.235477, -31.2227, 3447.34, 4.54408e-5, -3.25374e-9};
    kappa = pkappa[0] + pkappa[1] / mll + pkappa[2] / sqr(mll) + pkappa[3] * mll + pkappa[4] * sqr(mll);
  }

  // Do the smearing
  double mll_scale = -1;
  do {
    mll_scale = (rand01() > kappa) ? randnorm(muG, sigG) : randcrystalball(alpCB, nCB, muCB, sigCB);
  } while (fabs(mll_scale) > 0.5);
  const double mll_reco = (1 + mll_scale) * mll;

  // Require high-Mll to avoid the Z peak
  if (mll_reco < 225 * GeV) vetoEvent;

  // Fill Mll histograms
  (mumu ? _h_mmm : _h_mee)->fill(mll_reco / GeV);
}


/// Normalise histograms etc., after the run.
//  The plot is expressed as number of events in a 10 GeV bin.
//  Multiply by luminosity*cross section (in same units!) to get number
//  of events a 1 GeV bin, then by 10 to get number of events in a 10 GeV.
void finalize() {
  scale(_h_mee, 10. * crossSection() / picobarn * luminosity() / sumOfWeights());
  scale(_h_mmm, 10. * crossSection() / picobarn * luminosity() / sumOfWeights());
}

/// @}


/// @name Histograms
/// @{
Histo1DPtr _h_mee, _h_mmm;
/// @}

};

RIVET_DECLARE_PLUGIN(ATLAS_2019_I1725190);

} ```