Rivet analyses

title: ATLAS_2014_I1312627

Ratios of $V$+jets observables between $W$ and $Z$ events

Experiment: ATLAS (LHC)

Inspire ID: 1312627

Status: VALIDATED

Authors: - Christian Gutschow

References: - Expt page: ATLAS-STDM-2012-25 - Eur.Phys.J. C74 (2014) 3168 - DOI: 10.1140/epjc/s10052-014-3168-9 - arXiv: 1408.6510

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - Inclusive $W$ and inclusive $Z$ in both electron and muon channels

Measurements of the ratio of the production cross sections for $W$ and $Z$ bosons in association with jets in proton-proton collisions at $\sqrt{s} = 7$~TeV with the ATLAS experiment at the Large Hadron Collider. The measurement is based on the entire 2011 dataset, corresponding to an integrated luminosity of 4.6$\text{fb}^{-1}$. Inclusive and differential cross-section ratios for massive vector bosons decaying to electrons and muons are measured in association with jets with transverse momentum $p_\text{T} > 30$~GeV and jet rapidity $|y| < 4.4$. The default routine will pick up the electron decay channel of the heavy bosons and compare it to the combined (muon and electron channel) data. Individual channels (for data) are available as well, use ATLAS_2014_I1312627_EL and ATLAS_2014_I1312627_MU to specify the decay channel directly. NB #1: The "x01" Scatter2D objects are constructed from the ratio of "x02" to "x03" Histo1D objects. If several output yoda files are merged with yodamerge, the merged "x01" objects will become meaningless. New "x01" Scatter2Ds can easil be constructed in a postprocessing step from the merged "x02" (nominator) and "x03" (denominator) objects. NB #2: Special care ought to be taken when evaluating theoretical uncertainties due to potential cancellations/correlations between numerator and denominator.

Source code:ATLAS_2014_I1312627.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/DileptonFinder.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/LeptonFinder.hh"

include "Rivet/Projections/MissingMomentum.hh"

include "Rivet/Projections/PromptFinalState.hh"

include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {

/// Measurement of V+jets distributions, taken as ratios between W and Z channels class ATLAS_2014_I1312627 : public Analysis { public:

/// @name Plotting helper types
/// @{

struct Plots {
  string ref;
  Histo1DPtr comp[2];  // (de)nominator components
  Estimate1DPtr ratio; // Rjets plot
};

/// @}


/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2014_I1312627);


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

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

  // get option
  _mode = 0;
  if (getOption("LMODE") == "EL") _mode = 1;
  if (getOption("LMODE") == "MU") _mode = 2;

  // Set up cuts
  Cut cuts;
  if (_mode == 2) { // muon channel
    cuts = Cuts::pT > 25 * GeV && Cuts::abseta < 2.4;
    ;
  }
  else if (_mode) { // electron channel
    cuts = Cuts::pT > 25 * GeV && (Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47));
  }
  else { // combined data extrapolated to common phase space
    cuts = Cuts::pT > 25 * GeV && Cuts::abseta < 2.5;
  }

  // Boson finders
  declare("MET", MissingMomentum());
  LeptonFinder lf(0.1, cuts && Cuts::abspid == (_mode > 1 ? PID::MUON : PID::ELECTRON));
  declare(lf, "Leptons");
  DileptonFinder zf(91.2 * GeV, 0.1, cuts && Cuts::abspid == (_mode > 1 ? PID::MUON : PID::ELECTRON),
                    Cuts::massIn(66 * GeV, 116 * GeV));
  declare(zf, "ZF");

  // Jets
  VetoedFinalState jet_fs;
  jet_fs.addVetoOnThisFinalState(lf);
  jet_fs.addVetoOnThisFinalState(zf);
  FastJets jets(jet_fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL); //< !!
  declare(jets, "Jets");


  // Book Rjets plots
  _suff = string("-y0") + to_str(_mode + 1);
  hInit("Njets_incl", "d01");  // inclusive number of jets
  hInit("Njets_excl", "d04");  // exclusive number of jets
  hInit("Pt1_N1incl", "d05");  // leading jet pT, at least 1 jet
  hInit("Pt1_N1excl", "d06");  // leading jet pT, exactly 1 jet
  hInit("Pt1_N2incl", "d07");  // leading jet pT, at least 2 jets
  hInit("Pt1_N3incl", "d08");  // leading jet pT, at least 3 jets
  hInit("Pt2_N2incl", "d09");  // subleading jet pT, at least 2 jets
  hInit("Pt3_N3incl", "d10");  // sub-subleading jet pT, at least 3 jets
  hInit("ST_N2incl", "d11");   // scalar jet pT sum, at least 2 jets
  hInit("ST_N2excl", "d12");   // scalar jet pT sum, exactly 2 jets
  hInit("ST_N3incl", "d13");   // scalar jet pT sum, at least 3 jets
  hInit("ST_N3excl", "d14");   // scalar jet pT sum, exactly 3 jets
  hInit("DR_N2incl", "d15");   // deltaR(j1, j2), at least 2 jets
  hInit("DPhi_N2incl", "d16"); // deltaPhi(j1, j2), at least 2 jets
  hInit("Mjj_N2incl", "d17");  // mjj, at least 2 jets
  hInit("Rap1_N1incl", "d18"); // leading jet rapidity, at least 1 jet
  hInit("Rap2_N2incl", "d19"); // subleading jet rapidity, at least 2 jets
  hInit("Rap3_N3incl", "d20"); // sub-subleading jet rapidity, at least 3 jets

  // Also book numerator and denominator for Rjets plots
  for (auto& item : _plots) {
    book(item.second.comp[0], item.second.ref + "2" + _suff, refData(item.second.ref + "1" + _suff));
    book(item.second.comp[1], item.second.ref + "3" + _suff, refData(item.second.ref + "1" + _suff));
  }
}


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

  // W reco, starting with MET
  const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
  const Particles& ls = apply<LeptonFinder>(event, "Leptons").particles();
  const Particles ls_mtfilt = select(ls, [&](const Particle& l) { return mT(l, pmiss) > 40 * GeV; });
  const int ifound = closestMatchIndex(ls_mtfilt, pmiss, Kin::mass, 80.4 * GeV);

  // Z reco
  const DileptonFinder& zf = apply<DileptonFinder>(event, "ZF");

  // Exit if no bosons found (note: overlapping W and Z reco is allowed)
  if (ifound < 0 && zf.empty()) vetoEvent;

  // Retrieve jets
  const JetFinder& jetfs = apply<JetFinder>(event, "Jets");
  Jets jets = jetfs.jetsByPt(Cuts::pT > 30 * GeV && Cuts::absrap < 4.4);

  // Apply boson cuts and fill histograms
  if (!zf.empty()) {
    const Particles& leptons = zf.constituents();
    if (oppSign(leptons[0], leptons[1]) && deltaR(leptons[0], leptons[1]) > 0.2)
      fillPlots(leptons, jets, 1);
  }
  if (ifound >= 0) {
    const Particle& lepton = ls_mtfilt[ifound];
    if (pmiss.pT() > 25 * GeV) fillPlots(Particles{lepton}, jets, 0);
  }
}


/// Normalise histograms etc., after the run
void finalize() {
  ///  Normalise, scale and otherwise manipulate histograms here
  const double sf(crossSection() / picobarn / sumOfWeights());
  for (auto& item : _plots) {
    scale(item.second.comp[0], sf);
    scale(item.second.comp[1], sf);
    divide(item.second.comp[0], item.second.comp[1], item.second.ratio);
  }
}

/// @}


/// Analysis helper functions
/// @{

/// Histogram filling function, to avoid duplication
void fillPlots(const Particles& leptons, Jets& jets, int isZ) {
  // Do jet-lepton overlap removal
  idiscardIfAnyDeltaRLess(jets, leptons, 0.5);

  // Calculate jet ST
  const size_t njets = jets.size();
  const double ST = sum(jets, pT, 0.0) / GeV;

  // Fill jet histos
  _plots["Njets_excl"].comp[isZ]->fill(njets + 0.5);
  for (size_t i = 0; i <= njets; ++i) _plots["Njets_incl"].comp[isZ]->fill(i + 0.5);

  if (njets > 0) {
    const double pT1 = jets[0].pT() / GeV;
    const double rap1 = jets[0].absrap();
    _plots["Pt1_N1incl"].comp[isZ]->fill(pT1);
    _plots["Rap1_N1incl"].comp[isZ]->fill(rap1);

    if (njets == 1) {
      _plots["Pt1_N1excl"].comp[isZ]->fill(pT1);
    }
    else if (njets > 1) {
      const double pT2 = jets[1].pT() / GeV;
      const double rap2 = jets[1].absrap();
      const double dR = deltaR(jets[0], jets[1]);
      const double dPhi = deltaPhi(jets[0], jets[1]);
      const double mjj = (jets[0].momentum() + jets[1].momentum()).mass() / GeV;
      _plots["Pt1_N2incl"].comp[isZ]->fill(pT1);
      _plots["Pt2_N2incl"].comp[isZ]->fill(pT2);
      _plots["Rap2_N2incl"].comp[isZ]->fill(rap2);
      _plots["DR_N2incl"].comp[isZ]->fill(dR);
      _plots["DPhi_N2incl"].comp[isZ]->fill(dPhi);
      _plots["Mjj_N2incl"].comp[isZ]->fill(mjj);
      _plots["ST_N2incl"].comp[isZ]->fill(ST);

      if (njets == 2) {
        _plots["ST_N2excl"].comp[isZ]->fill(ST);
      }
      else if (njets > 2) {
        const double pT3 = jets[2].pT() / GeV;
        const double rap3 = jets[2].absrap();
        _plots["Pt1_N3incl"].comp[isZ]->fill(pT1);
        _plots["Pt3_N3incl"].comp[isZ]->fill(pT3);
        _plots["Rap3_N3incl"].comp[isZ]->fill(rap3);
        _plots["ST_N3incl"].comp[isZ]->fill(ST);

        if (njets == 3) _plots["ST_N3excl"].comp[isZ]->fill(ST);
      }
    }
  }
}


/// Helper for histogram initialisation
void hInit(string label, string ident) {
  string pre = ident + "-x0";
  _plots[label].ref = pre;
  book(_plots[label].ratio, pre + "1" + _suff);
}

/// @}

protected:

// Data members
size_t _mode;
string _suff;

private:

/// @name Map of histograms
map<string, Plots> _plots;

};

RIVET_DECLARE_PLUGIN(ATLAS_2014_I1312627);

} ```