Rivet analyses

title: ATLAS_2011_I945498

$Z$+jets in $pp$ at 7~TeV

Experiment: ATLAS (LHC)

Inspire ID: 945498

Status: VALIDATED

Authors: - Evelin Meoni - Holger Schulz - Roman Lysak

References: - Expt page: ATLAS-STDM-2011-27 - arXiv: 1111.2690 - CERN-PH-EP-2011-162

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - Z+jets, electronic and/or muonic Z-decays. Jets with transverse momentum $p_T > 30$ GeV and jet rapidity $|y| < 4.4$.

Production of jets in association with a $Z/\gamma^$ boson in proton--proton collisions at $\sqrt{s} = 7$ TeV with the ATLAS detector. The analysis includes the full 2010 data set, collected with a low rate of multiple proton--proton collisions in the accelerator, corresponding to an integrated luminosity of 36 pb$^{-1}$. Inclusive jet cross sections in $Z/\gamma^$ events, with $Z/\gamma^*$ decaying into electron or muon pairs, are measured for jets with transverse momentum $p_T > 30$ GeV and jet rapidity $|y| < 4.4$.

Source code:ATLAS_2011_I945498.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/IdentifiedFinalState.hh"

include "Rivet/Projections/LeadingParticlesFinalState.hh"

include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {

/// ATLAS Z+jets in pp at 7 TeV class ATLAS_2011_I945498 : public Analysis { public:

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


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

  // Variable initialisation
  _isZeeSample = false;
  _isZmmSample = false;
  for (size_t chn = 0; chn < 3; ++chn) {
    book(weights_nj0[chn], "weights_nj0_" + to_str(chn));
    book(weights_nj1[chn], "weights_nj1_" + to_str(chn));
    book(weights_nj2[chn], "weights_nj2_" + to_str(chn));
    book(weights_nj3[chn], "weights_nj3_" + to_str(chn));
    book(weights_nj4[chn], "weights_nj4_" + to_str(chn));
  }

  // Set up projections
  Cut cuts_mu = Cuts::abseta < 2.4 && Cuts::pT > 20 * GeV;
  DileptonFinder zfinder_mu(91.2 * GeV, 0.1, cuts_mu && Cuts::abspid == PID::MUON,
                            Cuts::massIn(66 * GeV, 116 * GeV));
  declare(zfinder_mu, "DileptonFinder_mu");

  Cut cuts_e = (Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47)) && Cuts::pT > 20 * GeV;
  DileptonFinder zfinder_el(91.2 * GeV, 0.1, cuts_e && Cuts::abspid == PID::ELECTRON,
                            Cuts::massIn(66 * GeV, 116 * GeV));
  declare(zfinder_el, "DileptonFinder_el");

  // For combined cross-sections (combined phase space + dressed level)
  Cut cuts25_20 = Cuts::abseta < 2.5 && Cuts::pT > 20 * GeV;
  DileptonFinder zfinder_comb_mu(91.2 * GeV, 0.1, cuts25_20 && Cuts::abspid == PID::MUON,
                                 Cuts::massIn(66.0 * GeV, 116.0 * GeV));
  declare(zfinder_comb_mu, "DileptonFinder_comb_mu");
  DileptonFinder zfinder_comb_el(91.2 * GeV, 0.1, cuts25_20 && Cuts::abspid == PID::ELECTRON,
                                 Cuts::massIn(66.0 * GeV, 116.0 * GeV));
  declare(zfinder_comb_el, "DileptonFinder_comb_el");

  // Define veto FS in order to prevent Z-decay products entering the jet algorithm
  VetoedFinalState remfs;
  remfs.addVetoOnThisFinalState(zfinder_el);
  remfs.addVetoOnThisFinalState(zfinder_mu);
  VetoedFinalState remfs_comb;
  remfs_comb.addVetoOnThisFinalState(zfinder_comb_el);
  remfs_comb.addVetoOnThisFinalState(zfinder_comb_mu);

  FastJets jets(remfs, JetAlg::ANTIKT, 0.4);
  jets.useInvisibles();
  declare(jets, "jets");
  FastJets jets_comb(remfs_comb, JetAlg::ANTIKT, 0.4);
  jets_comb.useInvisibles();
  declare(jets_comb, "jets_comb");

  // 0=el, 1=mu, 2=comb
  for (size_t chn = 0; chn < 3; ++chn) {
    book(_h_njet_incl[chn], 1, 1, chn + 1);
    book(_h_njet_ratio[chn], 2, 1, chn + 1);
    book(_h_ptjet[chn], 3, 1, chn + 1);
    book(_h_ptlead[chn], 4, 1, chn + 1);
    book(_h_ptseclead[chn], 5, 1, chn + 1);
    book(_h_yjet[chn], 6, 1, chn + 1);
    book(_h_ylead[chn], 7, 1, chn + 1);
    book(_h_yseclead[chn], 8, 1, chn + 1);
    book(_h_mass[chn], 9, 1, chn + 1);
    book(_h_deltay[chn], 10, 1, chn + 1);
    book(_h_deltaphi[chn], 11, 1, chn + 1);
    book(_h_deltaR[chn], 12, 1, chn + 1);
  }
}


// Jet selection criteria universal for electron and muon channel
/// @todo Replace with a Cut passed to jetsByPt
Jets selectJets(const DileptonFinder* zf, const FastJets* allJets) {
  const FourMomentum l1 = zf->constituents()[0].momentum();
  const FourMomentum l2 = zf->constituents()[1].momentum();
  Jets jets;
  for (const Jet& jet : allJets->jetsByPt(Cuts::pT > 30 * GeV)) {
    const FourMomentum jmom = jet.momentum();
    if (jmom.absrap() < 4.4 && deltaR(l1, jmom) > 0.5 && deltaR(l2, jmom) > 0.5) {
      jets.push_back(jet);
    }
  }
  return jets;
}


/// Perform the per-event analysis
void analyze(const Event& event) {
  vector<const DileptonFinder*> zfs;
  zfs.push_back(&(apply<DileptonFinder>(event, "DileptonFinder_el")));
  zfs.push_back(&(apply<DileptonFinder>(event, "DileptonFinder_mu")));
  zfs.push_back(&(apply<DileptonFinder>(event, "DileptonFinder_comb_el")));
  zfs.push_back(&(apply<DileptonFinder>(event, "DileptonFinder_comb_mu")));

  vector<const FastJets*> fjs;
  fjs.push_back(&(apply<FastJets>(event, "jets")));
  fjs.push_back(&(apply<FastJets>(event, "jets_comb")));

  // Determine what kind of MC sample this is
  const bool isZee = (zfs[0]->bosons().size() == 1) || (zfs[2]->bosons().size() == 1);
  const bool isZmm = (zfs[1]->bosons().size() == 1) || (zfs[3]->bosons().size() == 1);
  if (isZee) _isZeeSample = true;
  if (isZmm) _isZmmSample = true;

  // Require exactly one electronic or muonic Z-decay in the event
  bool isZeemm = ((zfs[0]->bosons().size() == 1 && zfs[1]->bosons().size() != 1)
                  || (zfs[1]->bosons().size() == 1 && zfs[0]->bosons().size() != 1));
  bool isZcomb = ((zfs[2]->bosons().size() == 1 && zfs[3]->bosons().size() != 1)
                  || (zfs[3]->bosons().size() == 1 && zfs[2]->bosons().size() != 1));
  if (!isZeemm && !isZcomb) vetoEvent;

  vector<int> zfIDs;
  vector<int> fjIDs;
  if (isZeemm) {
    int chn = zfs[0]->bosons().size() == 1 ? 0 : 1;
    zfIDs.push_back(chn);
    fjIDs.push_back(0);
  }
  if (isZcomb) {
    int chn = zfs[2]->bosons().size() == 1 ? 2 : 3;
    zfIDs.push_back(chn);
    fjIDs.push_back(1);
  }

  for (size_t izf = 0; izf < zfIDs.size(); ++izf) {
    int zfID = zfIDs[izf];
    int fjID = fjIDs[izf];

    int chn = zfID;
    if (zfID == 2 || zfID == 3) chn = 2;

    Jets jets = selectJets(zfs[zfID], fjs[fjID]);

    switch (jets.size()) {
      case 0: weights_nj0[chn]->fill(); break;
      case 1:
        weights_nj0[chn]->fill();
        weights_nj1[chn]->fill();
        break;
      case 2:
        weights_nj0[chn]->fill();
        weights_nj1[chn]->fill();
        weights_nj2[chn]->fill();
        break;
      case 3:
        weights_nj0[chn]->fill();
        weights_nj1[chn]->fill();
        weights_nj2[chn]->fill();
        weights_nj3[chn]->fill();
        break;
      default: // >= 4
        weights_nj0[chn]->fill();
        weights_nj1[chn]->fill();
        weights_nj2[chn]->fill();
        weights_nj3[chn]->fill();
        weights_nj4[chn]->fill();
    }

    // Require at least one jet
    if (jets.empty()) continue;

    // Fill jet multiplicities
    for (size_t ijet = 1; ijet <= jets.size(); ++ijet) {
      _h_njet_incl[chn]->fill(ijet);
    }

    // Loop over selected jets, fill inclusive jet distributions
    for (size_t ijet = 0; ijet < jets.size(); ++ijet) {
      _h_ptjet[chn]->fill(jets[ijet].pT() / GeV);
      _h_yjet[chn]->fill(fabs(jets[ijet].rapidity()));
    }

    // Leading jet histos
    const double ptlead = jets[0].pT() / GeV;
    const double yabslead = fabs(jets[0].rapidity());
    _h_ptlead[chn]->fill(ptlead);
    _h_ylead[chn]->fill(yabslead);

    if (jets.size() >= 2) {
      // Second jet histos
      const double pt2ndlead = jets[1].pT() / GeV;
      const double yabs2ndlead = fabs(jets[1].rapidity());
      _h_ptseclead[chn]->fill(pt2ndlead);
      _h_yseclead[chn]->fill(yabs2ndlead);

      // Dijet histos
      const double deltaphi = fabs(deltaPhi(jets[1], jets[0]));
      const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity());
      const double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY));
      const double mass = (jets[0].momentum() + jets[1].momentum()).mass();
      _h_mass[chn]->fill(mass / GeV);
      _h_deltay[chn]->fill(deltarap);
      _h_deltaphi[chn]->fill(deltaphi);
      _h_deltaR[chn]->fill(deltar);
    }
  }
}


/// @name Ratio calculator util functions
/// @{

/// Calculate the ratio, being careful about div-by-zero
double ratio(double a, double b) {
  return (b != 0) ? a / b : 0;
}

/// Calculate the ratio error, being careful about div-by-zero
double ratio_err(double a, double b) {
  return (b != 0) ? sqrt(a / sqr(b) + sqr(a) / (b * b * b)) : 0;
}

/// @}


void finalize() {
  // Fill ratio histograms
  for (size_t chn = 0; chn < 3; ++chn) {
    _h_njet_ratio[chn]->bin(1).set(ratio(weights_nj1[chn]->val(), weights_nj0[chn]->val()),
                                   ratio_err(weights_nj1[chn]->val(), weights_nj0[chn]->val()));
    _h_njet_ratio[chn]->bin(2).set(ratio(weights_nj2[chn]->val(), weights_nj1[chn]->val()),
                                   ratio_err(weights_nj2[chn]->val(), weights_nj1[chn]->val()));
    _h_njet_ratio[chn]->bin(3).set(ratio(weights_nj3[chn]->val(), weights_nj2[chn]->val()),
                                   ratio_err(weights_nj3[chn]->val(), weights_nj2[chn]->val()));
    _h_njet_ratio[chn]->bin(4).set(ratio(weights_nj4[chn]->val(), weights_nj3[chn]->val()),
                                   ratio_err(weights_nj4[chn]->val(), weights_nj3[chn]->val()));
  }

  // Scale other histos
  for (size_t chn = 0; chn < 3; ++chn) {
    // For ee and mumu channels: normalize to Njet inclusive cross-section
    double xs = crossSectionPerEvent() / picobarn;
    if (chn != 2 && weights_nj0[chn]->val() != 0.) xs = 1.0 / weights_nj0[chn]->val();
    // For inclusive MC sample(ee/mmu channels together) we want the single-lepton-flavor xsec
    if (_isZeeSample && _isZmmSample) xs *= 0.5;

    // Special case histogram: always not normalized
    scale(_h_njet_incl[chn], (chn < 2) ? crossSectionPerEvent() / picobarn : xs);

    scale(_h_ptjet[chn], xs);
    scale(_h_ptlead[chn], xs);
    scale(_h_ptseclead[chn], xs);
    scale(_h_yjet[chn], xs);
    scale(_h_ylead[chn], xs);
    scale(_h_yseclead[chn], xs);
    scale(_h_deltaphi[chn], xs);
    scale(_h_deltay[chn], xs);
    scale(_h_deltaR[chn], xs);
    scale(_h_mass[chn], xs);
  }
}

/// @}

private:

bool _isZeeSample;
bool _isZmmSample;

CounterPtr weights_nj0[3];
CounterPtr weights_nj1[3];
CounterPtr weights_nj2[3];
CounterPtr weights_nj3[3];
CounterPtr weights_nj4[3];

Estimate1DPtr _h_njet_ratio[3];
Histo1DPtr _h_njet_incl[3];
Histo1DPtr _h_ptjet[3];
Histo1DPtr _h_ptlead[3];
Histo1DPtr _h_ptseclead[3];
Histo1DPtr _h_yjet[3];
Histo1DPtr _h_ylead[3];
Histo1DPtr _h_yseclead[3];
Histo1DPtr _h_deltaphi[3];
Histo1DPtr _h_deltay[3];
Histo1DPtr _h_deltaR[3];
Histo1DPtr _h_mass[3];

};

RIVET_DECLARE_PLUGIN(ATLAS_2011_I945498);

} ```