rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2016_I1449082

Charge asymmetry in top quark pair production in dilepton channel
Experiment: ATLAS (LHC)
Inspire ID: 1449082
Status: VALIDATED
Authors:
  • Roman Lysak
References: Beams: p+ p+
Beam energies: (4000.0, 4000.0) GeV
Run details:
  • pp --> top + antitop events at 8 TeV

Measurements of the top-antitop quark pair production charge asymmetry in the dilepton channel are presented using data corresponding to an integrated luminosity of $20.3 \text{fb}^{-1}$ from $pp$ collisions at a center-of-mass energy of $\sqrt{s} = 8$ TeV collected with the ATLAS detector at the Large Hadron Collider at CERN. Inclusive and differential measurements as a function of the invariant mass, transverse momentum, and longitudinal boost of the $t\bar{t}$ system are performed both in the full phase space and in a fiducial phase space closely matching the detector acceptance (at least 2 leptons and 2 jets with $p_\text{T} > 25$ GeV and $|\eta| < 2.5$). Two observables are studied: $A_C^{\ell\ell}$ based on the selected leptons and $A_C^{t\bar{t}}$ based on the reconstructed $t\bar{t}$ final state. The unfolded distributions of $\Delta|\eta| = |\eta|_{\ell^+} - |\eta|_{\ell^-}$ and $\Delta|y| = |y|_\text{top} - |y|_\text{antitop}$ are provided. USERS SHOULD NOTE THAT EXPLICIT RECONSTRUCTION OF INDIVIDUAL STANDARD MODEL NEUTRINOS IS USED IN THIS ANALYSIS ROUTINE TO MATCH THE MONTE-CARLO-BASED CORRECTION TO THE FIDUCIAL REGION APPLIED IN THE PAPER.

Source code: ATLAS_2016_I1449082.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/DressedLeptons.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// Charge asymmetry in top quark pair production in dilepton channel
 14  class ATLAS_2016_I1449082 : public Analysis {
 15  public:
 16
 17    const double MW = 80.300*GeV;
 18    const double MTOP = 172.5*GeV;
 19
 20    enum MeasureType { kInclMeas, kmttMeas, kbetaMeas, kptMeas, kNmeas };
 21    const size_t kNbins = 2;
 22    const double bins[kNmeas][3];
 23    const string measStr[kNmeas] = {"incl","mtt","beta","ptt"};
 24    const string rangeStr[kNmeas][2];
 25
 26
 27    /// Constructor
 28    //RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1449082);
 29    ATLAS_2016_I1449082() : Analysis("ATLAS_2016_I1449082"),
 30                            // inclusive (dummy), mtt [GeV], beta, pTtt
 31                            bins{ { 0., 1., 2. }, { 0., 500., 2000.}, { 0., 0.6 , 1.0}, { 0., 30. , 1000.} },
 32                            rangeStr{ { "0_1", "1_2"}, { "0_500", "500_2000"}, { "0_0.6", "0.6_1.0"}, { "0_30" , "30_1000"} }
 33    {  }
 34
 35    /// @name Analysis methods
 36    //@{
 37
 38    /// Book histograms and initialise projections before the run
 39    void init() {
 40
 41      // Cuts
 42      const Cut eta_full = Cuts::abseta < 5.0;
 43      const Cut lep_cuts = Cuts::abseta < 2.5 && Cuts::pT > 25*GeV;
 44      // All final state particles
 45      FinalState fs(eta_full);
 46      // Get photons to dress leptons
 47      IdentifiedFinalState photons(fs, PID::PHOTON);
 48
 49
 50
 51      // Electron projections
 52      // ---------------------
 53      // Electron/muons are defined from electron/muon and photons within a cone of DR = 0.1.
 54      // No isolation condition is imposed. The parent of the electron/muon is required not to be a hadron or quark.
 55      IdentifiedFinalState el_id(fs, {PID::ELECTRON,-PID::ELECTRON});
 56      PromptFinalState electrons(el_id);
 57      electrons.acceptTauDecays(true);
 58      // Electron dressing
 59      DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts, true);
 60      declare(dressedelectrons, "dressedelectrons");
 61      DressedLeptons dressedelectrons_full(photons, electrons, 0.1, eta_full, true);
 62
 63      // Muon projections
 64      // ---------------------
 65      IdentifiedFinalState mu_id(fs, {PID::MUON,-PID::MUON});
 66      PromptFinalState muons(mu_id);
 67      muons.acceptTauDecays(true);
 68      // Muon dressing
 69      DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts, true);
 70      declare(dressedmuons, "dressedmuons");
 71      DressedLeptons dressedmuons_full(photons, muons, 0.1, eta_full, true);
 72
 73      // Neutrino projections
 74      // ---------------------
 75      // Missing ET is calculated as the 4–vector sum of neutrinos from W/Z-boson decays. Tau decays are
 76      // included. A neutrino is treated as a detectable particle and is selected for consideration in the same
 77      // way as electrons or muons, i.e. the parent is required not to be a hadron or quark (u − b).
 78      IdentifiedFinalState nu_id;
 79      nu_id.acceptNeutrinos();
 80      PromptFinalState neutrinos(nu_id);
 81      neutrinos.acceptTauDecays(true);
 82      declare(neutrinos, "neutrinos");
 83
 84      // Jets projections
 85      // ---------------------
 86      // Jets are defined with the anti-kt algorithm, clustering all stable particles excluding the electrons,
 87      // muons, neutrinos, and photons used in the definition of the selected leptons.
 88      VetoedFinalState vfs(fs);
 89      vfs.addVetoOnThisFinalState(dressedelectrons_full);
 90      vfs.addVetoOnThisFinalState(dressedmuons_full);
 91      vfs.addVetoOnThisFinalState(neutrinos);
 92      declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "Jets");
 93
 94      // Book histograms
 95      book(_h_dEta     , 1, 1, 1);
 96      book(_h_dY       , 2, 1, 1);
 97      for (size_t iM = 0; iM < kNmeas; ++iM) {
 98        book(_h_Acll[iM], 3+iM, 1, 1);
 99        book(_h_Actt[iM], 7+iM, 1, 1);
100      }
101      for (size_t iM = 0; iM < kNmeas; ++iM) {
102        for (size_t iB = 0; iB < kNbins; ++iB) {
103          book(    _h_dEta_asym[iM][iB],  "_dEta_asym_" + measStr[iM] + "_bin" + rangeStr[iM][iB],  2,  -10., 10.);
104          book(    _h_dY_asym  [iM][iB],  "_dY_asym_"   + measStr[iM] + "_bin" + rangeStr[iM][iB],  2,  -10., 10.);
105        }
106      }
107
108    }
109
110
111    /// Perform the per-event analysis
112    void analyze(const Event& event) {
113
114      // Get the electrons and muons
115      const vector<DressedLepton> dressedelectrons = apply<DressedLeptons>(event, "dressedelectrons").dressedLeptons();
116      const vector<DressedLepton> dressedmuons     = apply<DressedLeptons>(event, "dressedmuons").dressedLeptons();
117      const vector<DressedLepton> leptons = dressedelectrons + dressedmuons;
118      // Require at least 2 leptons in the event
119      if (leptons.size() < 2) vetoEvent;
120
121      // Get the neutrinos
122      const Particles neutrinos = apply<PromptFinalState>(event, "neutrinos").particlesByPt();
123      // Require at least 2 neutrinos in the event (ick)
124      if (neutrinos.size() < 2)  vetoEvent;
125
126      // Get jets and apply selection
127      const Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
128      // Require at least 2 jets in the event
129      if (jets.size() < 2)  vetoEvent;
130
131      // Remaining selections
132      // Events where leptons and jets overlap, within dR = 0.4, are rejected.
133      for (const DressedLepton& lepton : leptons) {
134        if (any(jets, deltaRLess(lepton, 0.4))) vetoEvent;
135      }
136
137      // Construct pseudo-tops
138      // Exactly 2 opposite-sign leptons are required (e/mu)
139      if (leptons.size() != 2) vetoEvent;
140      if ( (leptons[0].charge() * leptons[1].charge()) > 0.) vetoEvent;
141      const FourMomentum lep_p = (leptons[0].charge3() > 0) ? leptons[0] : leptons[1];
142      const FourMomentum lep_n = (leptons[0].charge3() > 0) ? leptons[1] : leptons[0];
143
144      // Only the 2 leading pT selected neutrinos are considered
145      const FourMomentum nu1 = neutrinos[0].momentum();
146      const FourMomentum nu2 = neutrinos[1].momentum();
147
148      // Two jets correspond to the two leading jets in the event.
149      // If there is any b-tagged jet in the event, then the b-tagged jets
150      // are preferentially selected over the non-tagged jets without taking into account its pT.
151      // A jet is a b–jet if any B–hadron is included in the jet.
152      // Only B-hadrons with an initial pT > 5 GeV are considered.
153      Jets bjets, lightjets;
154      for (const Jet& jet : jets) {
155        (jet.bTagged(Cuts::pT > 5*GeV) ? bjets : lightjets) += jet;
156      }
157      // Already sorted by construction, since jets is sorted by decreasing pT
158      // std::sort(bjets.begin()    , bjets.end()    , cmpMomByPt);
159      // std::sort(lightjets.begin(), lightjets.end(), cmpMomByPt);
160
161      // Initially take 2 highest pT jets
162      FourMomentum bjet1 = jets[0];
163      FourMomentum bjet2 = jets[1];
164      if (!bjets.empty()) {
165        bjet1 = bjets[0];
166        bjet2 = (bjets.size() > 1) ? bjets[1] : lightjets[0]; //< We should have a light jet because >=2 jets requirement
167      } else {
168        // No btagged jets --> should have >= 2 light jets
169        bjet1 = lightjets[0];
170        bjet2 = lightjets[1];
171      }
172
173      // Construct pseudo-W bosons from lepton-neutrino combinations
174      // Minimize the difference between the mass computed from each lepton-neutrino combination and the W boson mass
175      const double massDiffW1 = fabs( (nu1 + lep_p).mass() - MW ) + fabs( (nu2 + lep_n).mass() - MW );
176      const double massDiffW2 = fabs( (nu1 + lep_n).mass() - MW ) + fabs( (nu2 + lep_p).mass() - MW );
177      const FourMomentum Wp = (massDiffW1 < massDiffW2) ? nu1+lep_p : nu2+lep_p;
178      const FourMomentum Wn = (massDiffW1 < massDiffW2) ? nu2+lep_n : nu1+lep_n;
179
180      // Construct pseudo-tops from jets and pseudo-W bosons
181      // Minimize the difference between the mass computed from each W-boson and b-jet combination and the top mass
182      const double massDiffT1 = fabs( (Wp+bjet1).mass()*GeV - MTOP ) + fabs( (Wn+bjet2).mass()*GeV - MTOP );
183      const double massDiffT2 = fabs( (Wp+bjet2).mass()*GeV - MTOP ) + fabs( (Wn+bjet1).mass()*GeV - MTOP );
184      const FourMomentum top_p = (massDiffT1 < massDiffT2) ? Wp+bjet1 : Wp+bjet2;
185      const FourMomentum top_n = (massDiffT1 < massDiffT2) ? Wn+bjet2 : Wn+bjet1;
186
187      // Calculate d|eta|, d|y|, etc.
188      double dEta = lep_p.abseta() - lep_n.abseta();
189      double dY   = top_p.absrapidity() - top_n.absrapidity();
190      double mtt  = (top_p + top_n).mass()*GeV;
191      double beta = fabs( (top_p + top_n).pz() ) / (top_p + top_n).E();
192      double pttt = (top_p + top_n).pt()*GeV;
193
194      // Fill histos, counters
195      _h_dEta->fill(dEta);
196      _h_dY  ->fill(dY  );
197      // Histos for inclusive and differential asymmetries
198      int mttBinID  = getBinID(kmttMeas , mtt);
199      int betaBinID = getBinID(kbetaMeas, beta);
200      int ptttBinID = getBinID(kptMeas  , pttt);
201      for (int iM = 0; iM < kNmeas; ++iM) {
202        int binID = -1;
203        switch (iM) {
204          case kInclMeas : binID = 0;         break;
205          case kmttMeas  : binID = mttBinID ; break;
206          case kbetaMeas : binID = betaBinID; break;
207          case kptMeas   : binID = ptttBinID; break;
208          default: binID = -1; break;
209        }
210        if (binID >= 0) {
211          _h_dY_asym  [iM][binID] ->fill(dY  );
212          _h_dEta_asym[iM][binID] ->fill(dEta);
213        }
214      }
215    }
216
217
218    /// Normalise histograms etc., after the run
219    void finalize() {
220
221      // Calculate charge asymmetries and fill them to hists
222      // Just for cross-check, calculate asymmetries from original dEta/dY histos
223      double asym = 0, err = 0;
224      calcAsymAndError(_h_dEta, asym, err);
225      MSG_INFO("Lepton inclusive asymmetry from histo:  = " << asym << " +- " << err );
226      calcAsymAndError(_h_dY, asym, err);
227      MSG_INFO("ttbar inclusive asymmetry from histo:  = "  << asym << " +- " << err );
228
229      // dEta/dY distributions: normalize to unity
230      normalize(_h_dEta);
231      normalize(_h_dY);
232
233      // Build asymm scatters
234      for (size_t iM = 0; iM < kNmeas; ++iM) {
235        for (size_t iB = 0; iB < kNbins; ++iB) {
236          // Only one bin for inclusive measurement
237          if ( (iM == kInclMeas) && (iB != 0)) continue;
238
239          calcAsymAndError(_h_dEta_asym[iM][iB], asym, err);
240          _h_Acll[iM]->addPoint( (bins[iM][iB+1] + bins[iM][iB])/2., asym, (bins[iM][iB+1] - bins[iM][iB])/2., err);
241
242          calcAsymAndError(_h_dY_asym[iM][iB], asym, err);
243          _h_Actt[iM]->addPoint( (bins[iM][iB+1] + bins[iM][iB])/2., asym, (bins[iM][iB+1] - bins[iM][iB])/2., err);
244        }
245      }
246    }
247
248    //@}
249
250
251  private:
252
253    void calcAsymAndError(Histo1DPtr hist, double& asym, double& err)  {
254
255      int nBins = hist->numBins();
256      if (nBins % 2 != 0) {
257      	asym = -999; err = -999.;
258        return;
259      }
260
261      double Nneg  = 0.;
262      double Npos  = 0.;
263      double dNneg = 0.;
264      double dNpos = 0.;
265      for (int iB = 0; iB < nBins; ++iB) {
266        if (iB < nBins/2) {
267          Nneg  += hist->bin(iB).sumW();
268          dNneg += hist->bin(iB).sumW2();
269        } else {
270          Npos  += hist->bin(iB).sumW();
271          dNpos += hist->bin(iB).sumW2();
272        }
273      }
274      dNneg = sqrt(dNneg);
275      dNpos = sqrt(dNpos);
276      asym = (Npos + Nneg) != 0.0 ? (Npos - Nneg) / (Npos + Nneg) : -999.;
277
278      const double Ntot = Npos + Nneg;
279      const double Ntot2 = Ntot * Ntot;
280      const double dNpos2 = dNpos * dNpos;
281      const double dNneg2 = dNneg * dNneg;
282      err = Ntot2 != 0. ? 2. * sqrt( (dNneg2 * Npos * Npos + dNpos2 * Nneg * Nneg) / (Ntot2 * Ntot2)) : -999.;
283    }
284
285
286    int getBinID(MeasureType type, double value)  {
287      /// @todo Use Rivet binIndex() function
288      for (size_t iBin = 0; iBin < kNbins; ++iBin) {
289        if (value <= bins[type][iBin+1]) return iBin;
290      }
291      return -1;
292    }
293
294
295    /// @name Histograms
296    //@{
297    Histo1DPtr _h_dEta;
298    Histo1DPtr _h_dY;
299    Scatter2DPtr _h_Actt[kNmeas];
300    Scatter2DPtr _h_Acll[kNmeas];
301    // Histograms to calculate the asymmetries from
302    /// @todo Use /TMP histos?
303    Histo1DPtr _h_dEta_asym[kNmeas][2];
304    Histo1DPtr _h_dY_asym  [kNmeas][2];
305    //@}
306
307    // Not-scaled histos
308    Histo1DPtr _h_dEta_notscaled, _h_dY_notscaled;
309
310  };
311
312
313  // The hook for the plugin system
314  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1449082);
315}