Rivet analyses referenceATLAS_2016_I1449082Charge asymmetry in top quark pair production in dilepton channelExperiment: ATLAS (LHC) Inspire ID: 1449082 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
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
256 int nBins = hist->numBins();
257 if (nBins % 2 != 0) {
258 asym = -999; err = -999.;
259 return;
260 }
261
262 double Nneg = 0.;
263 double Npos = 0.;
264 double dNneg = 0.;
265 double dNpos = 0.;
266 for (int iB = 0; iB < nBins; ++iB) {
267 if (iB < nBins/2) {
268 Nneg += hist->bin(iB).sumW();
269 dNneg += hist->bin(iB).sumW2();
270 } else {
271 Npos += hist->bin(iB).sumW();
272 dNpos += hist->bin(iB).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}
|