Rivet analyses referenceATLAS_2015_I1376945Colour flow in hadronic top decay at 8 TeVExperiment: ATLAS (LHC) Inspire ID: 1376945 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
The distribution and orientation of energy inside jets is predicted to be an experimental handle on colour connections between the hard-scatter quarks and gluons initiating the jets. This is a measurement of the distribution of one such variable, the jet pull angle. The pull angle is measured for jets produced in $t\bar{t}$ events with one $W$ boson decaying leptonically and the other decaying to jets using 20.3\,$\text{fb}^{-1}$ of data recorded with the ATLAS detector at a centre-of-mass energy of $\sqrt{s} = 8$\,TeV at the LHC. The jet pull angle distribution is corrected for detector resolution and acceptance effects. Source code: ATLAS_2015_I1376945.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/PromptFinalState.hh"
5#include "Rivet/Projections/IdentifiedFinalState.hh"
6#include "Rivet/Projections/LeptonFinder.hh"
7#include "Rivet/Projections/VetoedFinalState.hh"
8#include "Rivet/Projections/FastJets.hh"
9
10namespace Rivet {
11
12
13 /// @brief Colour flow in hadronic top decay at 8 TeV
14 class ATLAS_2015_I1376945 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1376945);
19
20
21 /// @name Analysis methods
22 /// @{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 const FinalState fs;
28
29 PromptFinalState promptFs(fs);
30 promptFs.acceptTauDecays(true);
31 promptFs.acceptMuonDecays(false);
32
33 IdentifiedFinalState neutrino_fs(promptFs);
34 neutrino_fs.acceptNeutrinos();
35 declare(neutrino_fs, "NEUTRINO_FS");
36
37 IdentifiedFinalState Photon(fs);
38 Photon.acceptIdPair(PID::PHOTON);
39
40 IdentifiedFinalState bare_muons_fs(promptFs);
41 bare_muons_fs.acceptIdPair(PID::MUON);
42
43 IdentifiedFinalState bare_elecs_fs(promptFs);
44 bare_elecs_fs.acceptIdPair(PID::ELECTRON);
45
46 Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 1*MeV);
47 LeptonFinder muons(bare_muons_fs, Photon, 0.1, lep_cuts);
48 declare(muons, "MUONS");
49
50 LeptonFinder elecs(bare_elecs_fs, Photon, 0.1, lep_cuts);
51 declare(elecs, "ELECS");
52
53 VetoedFinalState vfs;
54 vfs.addVetoOnThisFinalState(muons);
55 vfs.addVetoOnThisFinalState(elecs);
56 vfs.addVetoOnThisFinalState(neutrino_fs);
57
58 FastJets fjets(vfs, JetAlg::ANTIKT, 0.4);
59 fjets.useInvisibles();
60 declare(fjets, "jets");
61
62 book(h_pull_all ,4,1,1);
63 book(h_pull_charged ,5,1,1);
64 }
65
66
67 /// Perform the per-event analysis
68 void analyze(const Event& event) {
69
70 /**************
71 * JETS *
72 **************/
73 const Jets& allJets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25.0*GeV && Cuts::absrap < 2.5);
74 const DressedLeptons& all_elecs = apply<LeptonFinder>(event, "ELECS").dressedLeptons();
75 const DressedLeptons& all_muons = apply<LeptonFinder>(event, "MUONS").dressedLeptons();
76 Jets goodJets;
77 for (const Jet & j : allJets) {
78 bool keep = true;
79 for (const DressedLepton & el : all_elecs) keep &= deltaR(j, el) >= 0.2;
80 if (keep) goodJets += j;
81 }
82 if ( goodJets.size() < 4 ) vetoEvent;
83
84 /****************
85 * LEPTONS *
86 ****************/
87 DressedLeptons muons, vetoMuons;
88 for (const DressedLepton & mu : all_muons) {
89 bool keep = true;
90 for (const Jet & j : goodJets) keep &= deltaR(j, mu) >= 0.4;
91 if (keep && mu.pt() > 15*GeV) {
92 vetoMuons.push_back(mu);
93 if (mu.pt() > 25*GeV) muons.push_back(mu);
94 }
95 }
96
97 DressedLeptons elecs, vetoElecs;
98 for (const DressedLepton & el : all_elecs) {
99 bool keep = true;
100 for (const Jet & j : goodJets) keep &= deltaR(j, el) >= 0.4;
101 if (keep && el.pt() > 15*GeV) {
102 vetoElecs.push_back(el);
103 if (el.pt() > 25*GeV) elecs.push_back(el);
104 }
105 }
106
107 if (muons.empty() && elecs.empty()) vetoEvent;
108
109 bool muCandidate = !( muons.size() < 1 || vetoMuons.size() > 1 || vetoElecs.size() > 0 );
110 bool elCandidate = !( elecs.size() < 1 || vetoElecs.size() > 1 || vetoMuons.size() > 0 );
111
112 if (!elCandidate && !muCandidate) vetoEvent;
113
114 /******************************
115 * ELECTRON-MUON OVERLAP *
116 ******************************/
117 for (const DressedLepton & electron : elecs) {
118 for (const DressedLepton & muon : muons) {
119 double d_theta = fabs(muon.theta() - electron.theta());
120 double d_phi = deltaPhi(muon.phi(), electron.phi());
121 if (d_theta < 0.005 && d_phi < 0.005) vetoEvent;
122 }
123 }
124
125 /****************
126 * NEUTRINOS *
127 ****************/
128 const Particles& neutrinos = apply<IdentifiedFinalState>(event, "NEUTRINO_FS").particlesByPt();
129 FourMomentum metVector = FourMomentum(0.,0.,0.,0.);
130 for (const Particle& n : neutrinos) {
131 metVector += n.momentum();
132 }
133 double met = metVector.pt();
134 if (met <= 20*GeV) vetoEvent;
135
136 if ( (_mT(muCandidate? muons[0] : elecs[0], metVector) + met) <= 60. ) vetoEvent;
137
138 /****************
139 * B-JETS *
140 ****************/
141 Jets bJets, wJets;
142 for(Jet & j : goodJets) {
143 bool b_tagged = false;
144 Particles bTags = j.bTags();
145 for ( Particle & b : bTags ) {
146 b_tagged |= b.pT() > 5*GeV;
147 }
148 if (b_tagged) bJets += j;
149 if (!b_tagged && j.abseta() < 2.1) wJets += j;
150 }
151
152 if ( bJets.size() < 2 || wJets.size() < 2 ) vetoEvent;
153
154 double pull_angle = fabs(CalculatePullAngle(wJets[0], wJets[1], 0));
155 h_pull_all->fill(pull_angle / Rivet::PI);
156
157 double pull_angle_charged = fabs(CalculatePullAngle(wJets[0], wJets[1], 1));
158 h_pull_charged->fill(pull_angle_charged / Rivet::PI);
159
160 }
161
162 Vector3 CalculatePull(Jet& jet, bool &isCharged) {
163 Vector3 pull(0.0, 0.0, 0.0);
164 double PT = jet.pT();
165 Particles& constituents = jet.particles();
166 Particles charged_constituents;
167 if (isCharged) {
168 for (Particle & p : constituents) {
169 if (p.charge3() != 0) charged_constituents += p;
170 }
171 constituents = charged_constituents;
172 }
173 // calculate axis
174 FourMomentum axis;
175 for (Particle& p : constituents) axis += p.momentum();
176 Vector3 J(axis.rap(), axis.phi(MINUSPI_PLUSPI), 0.0);
177 // calculate pull
178 for (Particle & p : constituents) {
179 Vector3 ri = Vector3(p.rap(), p.phi(MINUSPI_PLUSPI), 0.0) - J;
180 while (ri.y() > Rivet::PI) ri.setY(ri.y() - Rivet::TWOPI);
181 while (ri.y() < -Rivet::PI) ri.setY(ri.y() + Rivet::TWOPI);
182 pull.setX(pull.x() + (ri.mod() * ri.x() * p.pT()) / PT);
183 pull.setY(pull.y() + (ri.mod() * ri.y() * p.pT()) / PT);
184 }
185 return pull;
186 }
187
188 double CalculatePullAngle(Jet& jet1, Jet& axisjet, bool isCharged) {
189 Vector3 pull_vector = CalculatePull(jet1, isCharged);
190 pull_vector = Vector3(1000.*pull_vector.x(), 1000.*pull_vector.y(), 0.);
191 double drap = axisjet.rap() - jet1.rap();
192 double dphi = axisjet.phi(MINUSPI_PLUSPI) - jet1.phi(MINUSPI_PLUSPI);
193 Vector3 j2_vector(drap, dphi, 0.0);
194 return mapAngleMPiToPi(deltaPhi(pull_vector, j2_vector));
195 }
196
197 double _mT(const FourMomentum &l, FourMomentum &nu) const {
198 return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) );
199 }
200
201 /// Normalise histograms etc., after the run
202 void finalize() {
203 normalize(h_pull_all);
204 normalize(h_pull_charged);
205 }
206
207 /// @}
208
209
210 private:
211
212 Histo1DPtr h_pull_all;
213 Histo1DPtr h_pull_charged;
214
215 };
216
217
218 RIVET_DECLARE_PLUGIN(ATLAS_2015_I1376945);
219
220}
|