Rivet analyses

title: ATLAS_2012_I1203852

Measurement of the $ZZ(*)$ production cross-section in $pp$ collisions at 7 TeV with ATLAS

Experiment: ATLAS (LHC)

Inspire ID: 1203852

Status: VALIDATED

Authors: - Oldrich Kepka - Katerina Moudra

References: - Expt page: ATLAS-STDM-2012-02 - arXiv: 1211.6096

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - Run with inclusive $Z$ events, with $Z$ decays to 4 leptons or 2 leptons + MET.

Measurement of the fiducial cross section for $ZZ(*)$ production in proton proton collisions at a centre-of mass energy of 7 TeV, is presented, using data corresponding to an integrated luminosity of 4.6/fb collected by the ATLAS experiment at the Large Hadron Collider. The cross-section is measured using processes with two $Z$ bosons decaying to electrons or muons or with one $Z$ boson decaying to electrons or muons and a second $Z$ boson decaying to neutrinos. The fiducial region contains dressed leptons in restricted $p_T$ and $\eta$ ranges. The selection has specific requirements for both production processes. A measurement of the normalized fiducial cross-section as a function of $ZZ$ invariant mass, leading $Z$ $p_T$ and angle of two leptons coming from the leading $Z$ is also presented for both signal processes.

Source code:ATLAS_2012_I1203852.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/IdentifiedFinalState.hh"

include "Rivet/Projections/InvMassFinalState.hh"

include "Rivet/Projections/LeptonFinder.hh"

include "Rivet/Projections/MergedFinalState.hh"

include "Rivet/Projections/MissingMomentum.hh"

include "Rivet/Projections/PromptFinalState.hh"

include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {

/// Generic Z candidate struct Zstate : public ParticlePair { Zstate() { } Zstate(ParticlePair _particlepair) : ParticlePair(_particlepair) { } FourMomentum mom() const { return first.momentum() + second.momentum(); } operator FourMomentum() const { return mom(); } static bool cmppT(const Zstate& lx, const Zstate& rx) { return lx.mom().pT() < rx.mom().pT(); } };

/// ZZ analysis class ATLAS_2012_I1203852 : public Analysis { public:

/// Default constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1203852);


void init() {

  // Get options
  // Default does everything
  _mode = 0;
  if (getOption("LMODE") == "LL") _mode = 1;
  if (getOption("LMODE") == "LNU") _mode = 2;

  // NB Missing ET is not required to be neutrinos
  FinalState fs(Cuts::abseta < 5.0);
  PromptFinalState pfs(fs);

  // Final states to form Z bosons
  vids.push_back(make_pair(PID::ELECTRON, PID::POSITRON));
  vids.push_back(make_pair(PID::MUON, PID::ANTIMUON));

  if (_mode != 2) {

    // Selection 1: ZZ-> llll selection
    Cut etaranges_lep = Cuts::abseta < 3.16 && Cuts::pT > 7 * GeV;

    LeptonFinder electron_sel4l(0.1, etaranges_lep && Cuts::abspid == PID::ELECTRON);
    declare(electron_sel4l, "ELECTRON_sel4l");
    LeptonFinder muon_sel4l(0.1, etaranges_lep && Cuts::abspid == PID::MUON);
    declare(muon_sel4l, "MUON_sel4l");

    // Both ZZ on-shell histos
    book(_h_ZZ_xsect, 1, 1, 1);
    book(_h_ZZ_ZpT, 3, 1, 1);
    book(_h_ZZ_phill, 5, 1, 1);
    book(_h_ZZ_mZZ, 7, 1, 1);

    // One Z off-shell (ZZstar) histos
    book(_h_ZZs_xsect, 1, 1, 2);
  }

  if (_mode != 1) {

    // Selection 2: ZZ-> llnunu selection
    Cut etaranges_lep2 = Cuts::abseta < 2.5 && Cuts::pT > 10 * GeV;

    LeptonFinder electron_sel2l2nu(0.1, etaranges_lep2 && Cuts::abspid == PID::ELECTRON);
    declare(electron_sel2l2nu, "ELECTRON_sel2l2nu");
    LeptonFinder muon_sel2l2nu(0.1, etaranges_lep2 && Cuts::abspid == PID::MUON);
    declare(muon_sel2l2nu, "MUON_sel2l2nu");

    /// Get all neutrinos. These will not be used to form jets.
    IdentifiedFinalState neutrino_fs(Cuts::abseta < 4.5);
    neutrino_fs.acceptNeutrinos();
    declare(neutrino_fs, "NEUTRINO_FS");

    // Calculate missing ET from the visible final state, not by requiring neutrinos
    declare(MissingMomentum(Cuts::abseta < 4.5), "MISSING");

    VetoedFinalState jetinput;
    jetinput.addVetoOnThisFinalState(neutrino_fs);

    FastJets jetpro(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE);
    declare(jetpro, "jet");

    // ZZ -> llnunu histos
    book(_h_ZZnunu_xsect, 1, 1, 3);
    book(_h_ZZnunu_ZpT, 4, 1, 1);
    book(_h_ZZnunu_phill, 6, 1, 1);
    book(_h_ZZnunu_mZZ, 8, 1, 1);
  }
}


/// Do the analysis
void analyze(const Event& e) {

  if (_mode != 2) {

    ////////////////////////////////////////////////////////////////////
    // preselection of leptons for ZZ-> llll final state
    ////////////////////////////////////////////////////////////////////

    Particles leptons_sel4l;

    const DressedLeptons& mu_sel4l = apply<LeptonFinder>(e, "MUON_sel4l").dressedLeptons();
    const DressedLeptons& el_sel4l = apply<LeptonFinder>(e, "ELECTRON_sel4l").dressedLeptons();

    DressedLeptons leptonsFS_sel4l;
    leptonsFS_sel4l.insert(leptonsFS_sel4l.end(), mu_sel4l.begin(), mu_sel4l.end());
    leptonsFS_sel4l.insert(leptonsFS_sel4l.end(), el_sel4l.begin(), el_sel4l.end());

    ////////////////////////////////////////////////////////////////////
    // OVERLAP removal dR(l,l)>0.2
    ////////////////////////////////////////////////////////////////////
    for (const DressedLepton& l1 : leptonsFS_sel4l) {
      bool isolated = true;
      for (DressedLepton& l2 : leptonsFS_sel4l) {
        const double dR = deltaR(l1, l2);
        if (dR < 0.2 && !isSame(l1, l2)) {
          isolated = false;
          break;
        }
      }
      if (isolated) leptons_sel4l.push_back(l1);
    }

    //////////////////////////////////////////////////////////////////
    // Exactly two opposite charged leptons
    //////////////////////////////////////////////////////////////////

    // calculate total 'flavour' charge
    double totalcharge = 0;
    for (const Particle& l : leptons_sel4l) totalcharge += l.pid();

    // Analyze 4 lepton events
    if (leptons_sel4l.size() == 4 && totalcharge == 0) {
      Zstate Z1, Z2;

      // Identifies Z states from 4 lepton pairs
      identifyZstates(Z1, Z2, leptons_sel4l);

      ////////////////////////////////////////////////////////////////////////////
      // Z MASS WINDOW
      //  -ZZ: for both Z: 66<mZ<116 GeV
      //  -ZZ*: one Z on-shell: 66<mZ<116 GeV, one Z off-shell: mZ>20 GeV
      ///////////////////////////////////////////////////////////////////////////

      Zstate leadPtZ = std::max(Z1, Z2, Zstate::cmppT);

      double mZ1 = Z1.mom().mass();
      double mZ2 = Z2.mom().mass();
      double ZpT = leadPtZ.mom().pT();
      double phill = fabs(deltaPhi(leadPtZ.first, leadPtZ.second));
      if (phill > M_PI) phill = 2 * M_PI - phill;
      double mZZ = (Z1.mom() + Z2.mom()).mass();

      if (mZ1 > 20 * GeV && mZ2 > 20 * GeV) {
        // ZZ* selection
        if (inRange(mZ1, 66 * GeV, 116 * GeV) || inRange(mZ2, 66 * GeV, 116 * GeV)) {
          _h_ZZs_xsect->fill(sqrtS() * GeV); ///< @todo xsec * GeV??
        }

        // ZZ selection
        if (inRange(mZ1, 66 * GeV, 116 * GeV) && inRange(mZ2, 66 * GeV, 116 * GeV)) {
          _h_ZZ_xsect->fill(sqrtS() / GeV); ///< @todo xsec * GeV??
          _h_ZZ_ZpT->fill(ZpT);
          _h_ZZ_phill->fill(phill);
          _h_ZZ_mZZ->fill(mZZ);
        }
      }
    }
  }

  if (_mode != 1) {

    ////////////////////////////////////////////////////////////////////
    /// preselection of leptons for ZZ-> llnunu final state
    ////////////////////////////////////////////////////////////////////

    Particles leptons_sel2l2nu; // output
    const DressedLeptons& mu_sel2l2nu = apply<LeptonFinder>(e, "MUON_sel2l2nu").dressedLeptons();
    const DressedLeptons& el_sel2l2nu = apply<LeptonFinder>(e, "ELECTRON_sel2l2nu").dressedLeptons();

    DressedLeptons leptonsFS_sel2l2nu;
    leptonsFS_sel2l2nu.insert(leptonsFS_sel2l2nu.end(), mu_sel2l2nu.begin(), mu_sel2l2nu.end());
    leptonsFS_sel2l2nu.insert(leptonsFS_sel2l2nu.end(), el_sel2l2nu.begin(), el_sel2l2nu.end());

    // Lepton preselection for ZZ-> llnunu
    if ((mu_sel2l2nu.empty() || el_sel2l2nu.empty()) // cannot have opposite flavour
        && (leptonsFS_sel2l2nu.size() == 2)          // exactly two leptons
        && (leptonsFS_sel2l2nu[0].charge() * leptonsFS_sel2l2nu[1].charge() < 1) // opposite charge
        && (deltaR(leptonsFS_sel2l2nu[0], leptonsFS_sel2l2nu[1]) > 0.3)          // overlap removal
        && (leptonsFS_sel2l2nu[0].pT() > 20 * GeV
            && leptonsFS_sel2l2nu[1].pT() > 20 * GeV)) { // trigger requirement
      leptons_sel2l2nu.insert(leptons_sel2l2nu.end(), leptonsFS_sel2l2nu.begin(),
                              leptonsFS_sel2l2nu.end());
    }
    if (leptons_sel2l2nu.empty()) vetoEvent; // no further analysis, fine to veto

    Particles leptons_sel2l2nu_jetveto;
    for (const DressedLepton& l : mu_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.bareLepton());
    for (const DressedLepton& l : el_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.bareLepton());
    double ptll = (leptons_sel2l2nu[0].momentum() + leptons_sel2l2nu[1].momentum()).pT();

    // Find Z1-> ll
    FinalState fs2((Cuts::etaIn(-3.2, 3.2)));
    InvMassFinalState imfs(fs2, vids, 20 * GeV, sqrtS());
    imfs.calc(leptons_sel2l2nu);
    if (imfs.particlePairs().size() != 1) vetoEvent;
    const ParticlePair& Z1constituents = imfs.particlePairs()[0];
    FourMomentum Z1 = Z1constituents.first.momentum() + Z1constituents.second.momentum();

    // Z to neutrinos candidate from missing ET
    const MissingMomentum& missmom = apply<MissingMomentum>(e, "MISSING");
    const FourMomentum Z2 = missmom.missingMomentum(ZMASS);
    double met_Znunu = missmom.missingEt(); //Z2.pT();

    // mTZZ
    const double mT2_1st_term = add_quad(ZMASS, ptll) + add_quad(ZMASS, met_Znunu);
    const double mT2_2nd_term = Z1.px() + Z2.px();
    const double mT2_3rd_term = Z1.py() + Z2.py();
    const double mTZZ = sqrt(sqr(mT2_1st_term) - sqr(mT2_2nd_term) - sqr(mT2_3rd_term));

    if (!inRange(Z2.mass(), 66 * GeV, 116 * GeV)) vetoEvent;
    if (!inRange(Z1.mass(), 76 * GeV, 106 * GeV)) vetoEvent;

    /////////////////////////////////////////////////////////////
    // AXIAL MET < 75 GeV
    ////////////////////////////////////////////////////////////

    double dPhiZ1Z2 = fabs(deltaPhi(Z1, Z2));
    if (dPhiZ1Z2 > M_PI) dPhiZ1Z2 = 2 * M_PI - dPhiZ1Z2;
    const double axialEtmiss = -Z2.pT() * cos(dPhiZ1Z2);
    if (axialEtmiss < 75 * GeV) vetoEvent;

    const double ZpT = Z1.pT();
    double phill = fabs(deltaPhi(Z1constituents.first, Z1constituents.second));
    if (phill > M_PI) phill = 2 * M_PI - phill;


    ////////////////////////////////////////////////////////////////////////////
    // JETS
    //    -"j": found by "jetpro" projection && pT() > 25 GeV && |eta| < 4.5
    //    -"goodjets": "j"  && dR(electron/muon,jet) > 0.3
    //
    // JETVETO: veto all events with at least one good jet
    ///////////////////////////////////////////////////////////////////////////
    vector<Jet> good_jets;
    for (const Jet& j : apply<FastJets>(e, "jet").jetsByPt(Cuts::pT > 25 * GeV && Cuts::abseta < 4.5)) {
      bool isLepton = 0;
      for (const Particle& l : leptons_sel2l2nu_jetveto) {
        const double dR = deltaR(l.momentum(), j.momentum());
        if (dR < 0.3) {
          isLepton = true;
          break;
        }
      }
      if (!isLepton) good_jets.push_back(j);
    }
    size_t n_sel_jets = good_jets.size();
    if (n_sel_jets != 0) vetoEvent;


    /////////////////////////////////////////////////////////////
    // Fractional MET and lepton pair difference: "RatioMet"< 0.4
    ////////////////////////////////////////////////////////////
    double ratioMet = fabs(Z2.pT() - Z1.pT()) / Z1.pT();
    if (ratioMet > 0.4) vetoEvent;


    // End of ZZllnunu selection: now fill histograms
    _h_ZZnunu_xsect->fill(sqrtS() / GeV); ///< @todo xsec / GeV??
    _h_ZZnunu_ZpT->fill(ZpT);
    _h_ZZnunu_phill->fill(phill);
    _h_ZZnunu_mZZ->fill(mTZZ);
  }
}


/// Finalize
void finalize() {
  const double norm = crossSection() / sumOfWeights() / femtobarn;

  if (_mode != 2) {
    scale(_h_ZZ_xsect, norm);
    normalize(_h_ZZ_ZpT);
    normalize(_h_ZZ_phill);
    normalize(_h_ZZ_mZZ);
    scale(_h_ZZs_xsect, norm);
  }

  if (_mode != 1) {
    scale(_h_ZZnunu_xsect, norm);
    normalize(_h_ZZnunu_ZpT);
    normalize(_h_ZZnunu_phill);
    normalize(_h_ZZnunu_mZZ);
  }
}

protected:

// Data members like post-cuts event weight counters go here
size_t _mode;

private:

void identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& leptons_sel4l);
Histo1DPtr _h_ZZ_xsect, _h_ZZ_ZpT, _h_ZZ_phill, _h_ZZ_mZZ;
Histo1DPtr _h_ZZs_xsect;
Histo1DPtr _h_ZZnunu_xsect, _h_ZZnunu_ZpT, _h_ZZnunu_phill, _h_ZZnunu_mZZ;
vector<pair<PdgId, PdgId>> vids;
const double ZMASS = 91.1876; // GeV

};

/// 4l to ZZ assignment -- algorithm void ATLAS_2012_I1203852::identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& leptons_sel4l) {

/////////////////////////////////////////////////////////////////////////////
/// ZZ->4l pairing
/// - Exactly two same flavour opposite charged leptons
/// - Ambiguities in pairing are resolved by choosing the combination
///     that results in the smaller value of the sum |mll - mZ| for the two pairs
/////////////////////////////////////////////////////////////////////////////

Particles part_pos_el, part_neg_el, part_pos_mu, part_neg_mu;
for (const Particle& l : leptons_sel4l) {
  if (l.abspid() == PID::ELECTRON) {
    if (l.pid() < 0) part_neg_el.push_back(l);
    if (l.pid() > 0) part_pos_el.push_back(l);
  }
  else if (l.abspid() == PID::MUON) {
    if (l.pid() < 0) part_neg_mu.push_back(l);
    if (l.pid() > 0) part_pos_mu.push_back(l);
  }
}

// ee/mm channel
if (part_neg_el.size() == 2 || part_neg_mu.size() == 2) {

  Zstate Zcand_1, Zcand_2, Zcand_3, Zcand_4;
  if (part_neg_el.size() == 2) { // ee
    Zcand_1 = Zstate(ParticlePair(part_neg_el[0], part_pos_el[0]));
    Zcand_2 = Zstate(ParticlePair(part_neg_el[0], part_pos_el[1]));
    Zcand_3 = Zstate(ParticlePair(part_neg_el[1], part_pos_el[0]));
    Zcand_4 = Zstate(ParticlePair(part_neg_el[1], part_pos_el[1]));
  }
  else { // mumu
    Zcand_1 = Zstate(ParticlePair(part_neg_mu[0], part_pos_mu[0]));
    Zcand_2 = Zstate(ParticlePair(part_neg_mu[0], part_pos_mu[1]));
    Zcand_3 = Zstate(ParticlePair(part_neg_mu[1], part_pos_mu[0]));
    Zcand_4 = Zstate(ParticlePair(part_neg_mu[1], part_pos_mu[1]));
  }

  // We can have the following pairs: (Z1 + Z4) || (Z2 + Z3)
  double minValue_1, minValue_2;
  minValue_1 = fabs(Zcand_1.mom().mass() - ZMASS) + fabs(Zcand_4.mom().mass() - ZMASS);
  minValue_2 = fabs(Zcand_2.mom().mass() - ZMASS) + fabs(Zcand_3.mom().mass() - ZMASS);
  if (minValue_1 < minValue_2) {
    Z1 = Zcand_1;
    Z2 = Zcand_4;
  }
  else {
    Z1 = Zcand_2;
    Z2 = Zcand_3;
  }

  // emu channel
}
else if (part_neg_mu.size() == 1 && part_neg_el.size() == 1) {
  Z1 = Zstate(ParticlePair(part_neg_mu[0], part_pos_mu[0]));
  Z2 = Zstate(ParticlePair(part_neg_el[0], part_pos_el[0]));
}

}

RIVET_DECLARE_PLUGIN(ATLAS_2012_I1203852);

} ```