rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

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: 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
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/VetoedFinalState.hh"
  8#include "Rivet/Projections/LeptonFinder.hh"
  9#include "Rivet/Projections/MergedFinalState.hh"
 10#include "Rivet/Projections/MissingMomentum.hh"
 11#include "Rivet/Projections/InvMassFinalState.hh"
 12
 13namespace Rivet {
 14
 15
 16  /// Generic Z candidate
 17  struct Zstate : public ParticlePair {
 18    Zstate() { }
 19    Zstate(ParticlePair _particlepair) : ParticlePair(_particlepair) { }
 20    FourMomentum mom() const { return first.momentum() + second.momentum(); }
 21    operator FourMomentum() const { return mom(); }
 22    static bool cmppT(const Zstate& lx, const Zstate& rx) { return lx.mom().pT() < rx.mom().pT(); }
 23  };
 24
 25
 26
 27  /// ZZ analysis
 28  class ATLAS_2012_I1203852 : public Analysis {
 29  public:
 30
 31    /// Default constructor
 32    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1203852);
 33
 34
 35    void init() {
 36
 37      // Get options
 38      // Default does everything
 39      _mode = 0;
 40      if ( getOption("LMODE") == "LL" )  _mode = 1;
 41      if ( getOption("LMODE") == "LNU" ) _mode = 2;
 42
 43      // NB Missing ET is not required to be neutrinos
 44      FinalState fs(Cuts::abseta < 5.0);
 45      PromptFinalState pfs(fs);
 46
 47      // Final states to form Z bosons
 48      vids.push_back(make_pair(PID::ELECTRON, PID::POSITRON));
 49      vids.push_back(make_pair(PID::MUON, PID::ANTIMUON));
 50
 51      if (_mode!=2) {
 52
 53        // Selection 1: ZZ-> llll selection
 54        Cut etaranges_lep = Cuts::abseta < 3.16 && Cuts::pT > 7*GeV;
 55
 56        LeptonFinder electron_sel4l(0.1, etaranges_lep && Cuts::abspid == PID::ELECTRON);
 57        declare(electron_sel4l, "ELECTRON_sel4l");
 58        LeptonFinder muon_sel4l(0.1, etaranges_lep && Cuts::abspid == PID::MUON);
 59        declare(muon_sel4l, "MUON_sel4l");
 60
 61        // Both ZZ on-shell histos
 62        book(_h_ZZ_xsect ,1, 1, 1);
 63        book(_h_ZZ_ZpT   ,3, 1, 1);
 64        book(_h_ZZ_phill ,5, 1, 1);
 65        book(_h_ZZ_mZZ   ,7, 1, 1);
 66
 67        // One Z off-shell (ZZstar) histos
 68        book(_h_ZZs_xsect ,1, 1, 2);
 69
 70      }
 71
 72      if (_mode!=1) {
 73
 74        // Selection 2: ZZ-> llnunu selection
 75        Cut etaranges_lep2 = Cuts::abseta < 2.5 && Cuts::pT > 10*GeV;
 76
 77        LeptonFinder electron_sel2l2nu(0.1, etaranges_lep2 && Cuts::abspid == PID::ELECTRON);
 78        declare(electron_sel2l2nu, "ELECTRON_sel2l2nu");
 79        LeptonFinder muon_sel2l2nu(0.1, etaranges_lep2 && Cuts::abspid == PID::MUON);
 80        declare(muon_sel2l2nu, "MUON_sel2l2nu");
 81
 82        /// Get all neutrinos. These will not be used to form jets.
 83        IdentifiedFinalState neutrino_fs(Cuts::abseta < 4.5);
 84        neutrino_fs.acceptNeutrinos();
 85        declare(neutrino_fs, "NEUTRINO_FS");
 86
 87        // Calculate missing ET from the visible final state, not by requiring neutrinos
 88        declare(MissingMomentum(Cuts::abseta < 4.5), "MISSING");
 89
 90        VetoedFinalState jetinput;
 91        jetinput.addVetoOnThisFinalState(neutrino_fs);
 92
 93        FastJets jetpro(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE);
 94        declare(jetpro, "jet");
 95
 96        // ZZ -> llnunu histos
 97        book(_h_ZZnunu_xsect ,1, 1, 3);
 98        book(_h_ZZnunu_ZpT   ,4, 1, 1);
 99        book(_h_ZZnunu_phill ,6, 1, 1);
100        book(_h_ZZnunu_mZZ   ,8, 1, 1);
101
102      }
103
104    }
105
106
107    /// Do the analysis
108    void analyze(const Event& e) {
109
110      if (_mode != 2) {
111
112        ////////////////////////////////////////////////////////////////////
113        // preselection of leptons for ZZ-> llll final state
114        ////////////////////////////////////////////////////////////////////
115
116        Particles leptons_sel4l;
117
118        const DressedLeptons& mu_sel4l = apply<LeptonFinder>(e, "MUON_sel4l").dressedLeptons();
119        const DressedLeptons& el_sel4l = apply<LeptonFinder>(e, "ELECTRON_sel4l").dressedLeptons();
120
121        DressedLeptons leptonsFS_sel4l;
122        leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), mu_sel4l.begin(), mu_sel4l.end() );
123        leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), el_sel4l.begin(), el_sel4l.end() );
124
125        ////////////////////////////////////////////////////////////////////
126        // OVERLAP removal dR(l,l)>0.2
127        ////////////////////////////////////////////////////////////////////
128        for (const DressedLepton& l1 : leptonsFS_sel4l) {
129          bool isolated = true;
130          for (DressedLepton& l2 : leptonsFS_sel4l) {
131            const double dR = deltaR(l1, l2);
132            if (dR < 0.2 && !isSame(l1, l2)) { isolated = false; break; }
133          }
134          if (isolated) leptons_sel4l.push_back(l1);
135        }
136
137        //////////////////////////////////////////////////////////////////
138        // Exactly two opposite charged leptons
139        //////////////////////////////////////////////////////////////////
140
141        // calculate total 'flavour' charge
142        double totalcharge = 0;
143        for (const Particle& l : leptons_sel4l) totalcharge += l.pid();
144
145        // Analyze 4 lepton events
146        if (leptons_sel4l.size() == 4 && totalcharge == 0  ) {
147          Zstate Z1, Z2;
148
149          // Identifies Z states from 4 lepton pairs
150          identifyZstates(Z1, Z2,leptons_sel4l);
151
152          ////////////////////////////////////////////////////////////////////////////
153          // Z MASS WINDOW
154          //  -ZZ: for both Z: 66<mZ<116 GeV
155          //  -ZZ*: one Z on-shell: 66<mZ<116 GeV, one Z off-shell: mZ>20 GeV
156          ///////////////////////////////////////////////////////////////////////////
157
158          Zstate leadPtZ = std::max(Z1, Z2, Zstate::cmppT);
159
160          double mZ1   = Z1.mom().mass();
161          double mZ2   = Z2.mom().mass();
162          double ZpT   = leadPtZ.mom().pT();
163          double phill = fabs(deltaPhi(leadPtZ.first, leadPtZ.second));
164          if (phill > M_PI) phill = 2*M_PI-phill;
165          double mZZ   = (Z1.mom() + Z2.mom()).mass();
166
167          if (mZ1 > 20*GeV && mZ2 > 20*GeV) {
168            // ZZ* selection
169            if (inRange(mZ1, 66*GeV, 116*GeV) || inRange(mZ2, 66*GeV, 116*GeV)) {
170              _h_ZZs_xsect  -> fill(sqrtS()*GeV); ///< @todo xsec * GeV??
171            }
172
173            // ZZ selection
174            if (inRange(mZ1, 66*GeV, 116*GeV) && inRange(mZ2, 66*GeV, 116*GeV)) {
175              _h_ZZ_xsect  -> fill(sqrtS()/GeV); ///< @todo xsec * GeV??
176              _h_ZZ_ZpT    -> fill(ZpT);
177              _h_ZZ_phill  -> fill(phill);
178              _h_ZZ_mZZ    -> fill(mZZ);
179            }
180          }
181        }
182      }
183
184      if (_mode != 1) {
185
186        ////////////////////////////////////////////////////////////////////
187        /// preselection of leptons for ZZ-> llnunu final state
188        ////////////////////////////////////////////////////////////////////
189
190        Particles leptons_sel2l2nu; // output
191        const DressedLeptons& mu_sel2l2nu = apply<LeptonFinder>(e, "MUON_sel2l2nu").dressedLeptons();
192        const DressedLeptons& el_sel2l2nu = apply<LeptonFinder>(e, "ELECTRON_sel2l2nu").dressedLeptons();
193
194        DressedLeptons leptonsFS_sel2l2nu;
195        leptonsFS_sel2l2nu.insert( leptonsFS_sel2l2nu.end(), mu_sel2l2nu.begin(), mu_sel2l2nu.end() );
196        leptonsFS_sel2l2nu.insert( leptonsFS_sel2l2nu.end(), el_sel2l2nu.begin(), el_sel2l2nu.end() );
197
198        // Lepton preselection for ZZ-> llnunu
199        if ((mu_sel2l2nu.empty() || el_sel2l2nu.empty()) // cannot have opposite flavour
200            && (leptonsFS_sel2l2nu.size() == 2) // exactly two leptons
201            && (leptonsFS_sel2l2nu[0].charge() * leptonsFS_sel2l2nu[1].charge() < 1 ) // opposite charge
202            && (deltaR(leptonsFS_sel2l2nu[0], leptonsFS_sel2l2nu[1]) > 0.3) // overlap removal
203            && (leptonsFS_sel2l2nu[0].pT() > 20*GeV && leptonsFS_sel2l2nu[1].pT() > 20*GeV)) { // trigger requirement
204          leptons_sel2l2nu.insert(leptons_sel2l2nu.end(), leptonsFS_sel2l2nu.begin(), leptonsFS_sel2l2nu.end());
205        }
206        if (leptons_sel2l2nu.empty()) vetoEvent; // no further analysis, fine to veto
207
208        Particles leptons_sel2l2nu_jetveto;
209        for (const DressedLepton& l : mu_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.bareLepton());
210        for (const DressedLepton& l : el_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.bareLepton());
211        double ptll = (leptons_sel2l2nu[0].momentum() + leptons_sel2l2nu[1].momentum()).pT();
212
213        // Find Z1-> ll
214        FinalState fs2((Cuts::etaIn(-3.2, 3.2)));
215        InvMassFinalState imfs(fs2, vids, 20*GeV, sqrtS());
216        imfs.calc(leptons_sel2l2nu);
217        if (imfs.particlePairs().size() != 1) vetoEvent;
218        const ParticlePair& Z1constituents = imfs.particlePairs()[0];
219        FourMomentum Z1 = Z1constituents.first.momentum() + Z1constituents.second.momentum();
220
221        // Z to neutrinos candidate from missing ET
222        const MissingMomentum & missmom = apply<MissingMomentum>(e, "MISSING");
223        const FourMomentum Z2 = missmom.missingMomentum(ZMASS);
224        double met_Znunu = missmom.missingEt(); //Z2.pT();
225
226        // mTZZ
227        const double mT2_1st_term = add_quad(ZMASS, ptll) + add_quad(ZMASS, met_Znunu);
228        const double mT2_2nd_term = Z1.px() + Z2.px();
229        const double mT2_3rd_term = Z1.py() + Z2.py();
230        const double mTZZ = sqrt(sqr(mT2_1st_term) - sqr(mT2_2nd_term) - sqr(mT2_3rd_term));
231
232        if (!inRange(Z2.mass(), 66*GeV, 116*GeV)) vetoEvent;
233        if (!inRange(Z1.mass(), 76*GeV, 106*GeV)) vetoEvent;
234
235        /////////////////////////////////////////////////////////////
236        // AXIAL MET < 75 GeV
237        ////////////////////////////////////////////////////////////
238
239        double dPhiZ1Z2 = fabs(deltaPhi(Z1, Z2));
240        if (dPhiZ1Z2 > M_PI) dPhiZ1Z2 = 2*M_PI - dPhiZ1Z2;
241        const double axialEtmiss = -Z2.pT()*cos(dPhiZ1Z2);
242        if (axialEtmiss < 75*GeV) vetoEvent;
243
244        const double ZpT   = Z1.pT();
245        double phill = fabs(deltaPhi(Z1constituents.first, Z1constituents.second));
246        if (phill > M_PI) phill = 2*M_PI - phill;
247
248
249        ////////////////////////////////////////////////////////////////////////////
250        // JETS
251        //    -"j": found by "jetpro" projection && pT() > 25 GeV && |eta| < 4.5
252        //    -"goodjets": "j"  && dR(electron/muon,jet) > 0.3
253        //
254        // JETVETO: veto all events with at least one good jet
255        ///////////////////////////////////////////////////////////////////////////
256        vector<Jet> good_jets;
257        for (const Jet& j : apply<FastJets>(e, "jet").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 4.5)) {
258          bool isLepton = 0;
259          for (const Particle& l : leptons_sel2l2nu_jetveto) {
260            const double dR = deltaR(l.momentum(), j.momentum());
261            if (dR < 0.3) { isLepton = true; break; }
262          }
263          if (!isLepton) good_jets.push_back(j);
264        }
265        size_t n_sel_jets = good_jets.size();
266        if (n_sel_jets != 0) vetoEvent;
267
268
269        /////////////////////////////////////////////////////////////
270        // Fractional MET and lepton pair difference: "RatioMet"< 0.4
271        ////////////////////////////////////////////////////////////
272        double ratioMet = fabs(Z2.pT() - Z1.pT()) / Z1.pT();
273        if (ratioMet  > 0.4 ) vetoEvent;
274
275
276        // End of ZZllnunu selection: now fill histograms
277        _h_ZZnunu_xsect->fill(sqrtS()/GeV); ///< @todo xsec / GeV??
278        _h_ZZnunu_ZpT  ->fill(ZpT);
279        _h_ZZnunu_phill->fill(phill);
280        _h_ZZnunu_mZZ  ->fill(mTZZ);
281
282      }
283    }
284
285
286    /// Finalize
287    void finalize() {
288      const double norm = crossSection()/sumOfWeights()/femtobarn;
289
290      if (_mode!=2) {
291        scale(_h_ZZ_xsect, norm);
292        normalize(_h_ZZ_ZpT);
293        normalize(_h_ZZ_phill);
294        normalize(_h_ZZ_mZZ);
295        scale(_h_ZZs_xsect, norm);
296      }
297
298      if (_mode!=1) {
299        scale(_h_ZZnunu_xsect, norm);
300        normalize(_h_ZZnunu_ZpT);
301        normalize(_h_ZZnunu_phill);
302        normalize(_h_ZZnunu_mZZ);
303      }
304    }
305
306
307  protected:
308    // Data members like post-cuts event weight counters go here
309    size_t _mode;
310
311
312  private:
313
314    void identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& leptons_sel4l);
315    Histo1DPtr _h_ZZ_xsect, _h_ZZ_ZpT, _h_ZZ_phill, _h_ZZ_mZZ;
316    Histo1DPtr _h_ZZs_xsect;
317    Histo1DPtr _h_ZZnunu_xsect, _h_ZZnunu_ZpT, _h_ZZnunu_phill, _h_ZZnunu_mZZ;
318    vector< pair<PdgId,PdgId> > vids;
319    const double ZMASS = 91.1876; // GeV
320
321
322  };
323
324
325  /// 4l to ZZ assignment -- algorithm
326  void ATLAS_2012_I1203852::identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& leptons_sel4l) {
327
328    /////////////////////////////////////////////////////////////////////////////
329    /// ZZ->4l pairing
330    /// - Exactly two same flavour opposite charged leptons
331    /// - Ambiguities in pairing are resolved by choosing the combination
332    ///     that results in the smaller value of the sum |mll - mZ| for the two pairs
333    /////////////////////////////////////////////////////////////////////////////
334
335    Particles part_pos_el, part_neg_el, part_pos_mu, part_neg_mu;
336    for (const Particle& l : leptons_sel4l) {
337      if (l.abspid() == PID::ELECTRON) {
338        if (l.pid() < 0) part_neg_el.push_back(l);
339        if (l.pid() > 0) part_pos_el.push_back(l);
340      }
341      else if (l.abspid() == PID::MUON) {
342        if (l.pid() < 0) part_neg_mu.push_back(l);
343        if (l.pid() > 0) part_pos_mu.push_back(l);
344      }
345    }
346
347    // ee/mm channel
348    if ( part_neg_el.size() == 2 || part_neg_mu.size() == 2) {
349
350      Zstate Zcand_1, Zcand_2, Zcand_3, Zcand_4;
351      if (part_neg_el.size() == 2) { // ee
352        Zcand_1 = Zstate( ParticlePair( part_neg_el[0],  part_pos_el[0] ) );
353        Zcand_2 = Zstate( ParticlePair( part_neg_el[0],  part_pos_el[1] ) );
354        Zcand_3 = Zstate( ParticlePair( part_neg_el[1],  part_pos_el[0] ) );
355        Zcand_4 = Zstate( ParticlePair( part_neg_el[1],  part_pos_el[1] ) );
356      } else { // mumu
357        Zcand_1 = Zstate( ParticlePair( part_neg_mu[0],  part_pos_mu[0] ) );
358        Zcand_2 = Zstate( ParticlePair( part_neg_mu[0],  part_pos_mu[1] ) );
359        Zcand_3 = Zstate( ParticlePair( part_neg_mu[1],  part_pos_mu[0] ) );
360        Zcand_4 = Zstate( ParticlePair( part_neg_mu[1],  part_pos_mu[1] ) );
361      }
362
363      // We can have the following pairs: (Z1 + Z4) || (Z2 + Z3)
364      double minValue_1, minValue_2;
365      minValue_1 = fabs( Zcand_1.mom().mass() - ZMASS ) + fabs( Zcand_4.mom().mass() - ZMASS);
366      minValue_2 = fabs( Zcand_2.mom().mass() - ZMASS ) + fabs( Zcand_3.mom().mass() - ZMASS);
367      if (minValue_1 < minValue_2 ) {
368        Z1 = Zcand_1;
369        Z2 = Zcand_4;
370      } else {
371        Z1 = Zcand_2;
372        Z2 = Zcand_3;
373      }
374
375    // emu channel
376    } else if (part_neg_mu.size() == 1 && part_neg_el.size() == 1) {
377      Z1 = Zstate ( ParticlePair (part_neg_mu[0],  part_pos_mu[0] ) );
378      Z2 = Zstate ( ParticlePair (part_neg_el[0],  part_pos_el[0] ) );
379    }
380
381  }
382
383
384  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1203852);
385
386}