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