Rivet analyses
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 pT and η 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 pT and angle of two leptons coming from the leading Z is also presented for both signal processes.
Source
code:ATLAS_2012_I1203852.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/LeptonFinder.hh"
#include "Rivet/Projections/MergedFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/InvMassFinalState.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);
}