Rivet analyses referenceATLAS_2017_I1509919Track-based underlying event at 13 TeV in ATLASExperiment: ATLAS (LHC) Inspire ID: 1509919 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
We present charged-particle distributions sensitive to the underlying event, measured by the ATLAS detector in proton-proton collisions at a centre-of-mass energy of 13 TeV, in low-luminosity Large Hadron Collider fills corresponding to an integrated luminosity of 1.6 nb${}^{-1}$. The distributions were constructed using charged particles with absolute pseudorapidity less than 2.5 and with transverse momentum greater than 500 MeV, in events with at least one such charged particle with transverse momentum above 1 GeV. These distributions characterise the angular distribution of energy and particle flows with respect to the charged particle with highest transverse momentum, as a function of both that momentum and of charged-particle multiplicity. The results have been corrected for detector effects and are compared to the predictions of various Monte Carlo event generators, experimentally establishing the level of underlying-event activity at LHC Run 2 energies and providing inputs for the development of event generator modelling. The current models in use for UE modelling typically describe this data to 5% accuracy, compared with data uncertainties of less than 1%. Source code: ATLAS_2017_I1509919.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4
5namespace Rivet {
6
7
8 /// Track-based underlying event at 13 TeV in ATLAS
9 class ATLAS_2017_I1509919 : public Analysis {
10 public:
11
12 // Constructor
13 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1509919);
14
15
16 // Pre-run histogram and projection booking
17 void init() {
18
19 declare(ChargedFinalState(Cuts::abseta < 2.5 && Cuts::pT > 500*MeV), "CFS500");
20
21 // Nch profiles vs. pT_lead
22 book(_hist_nch[0], 22, 1, 1);
23 book(_hist_nch[1], 23, 1, 1);
24 book(_hist_nch[2], 21, 1, 1);
25 book(_hist_nch[3], 3, 1, 1);
26 book(_hist_nch[4], 2, 1, 1);
27 book(_hist_nch[5], 4, 1, 1);
28
29 // pTsum profiles vs. pT_lead
30 book(_hist_ptsum[0], 25, 1, 1);
31 book(_hist_ptsum[1], 26, 1, 1);
32 book(_hist_ptsum[2], 24, 1, 1);
33 book(_hist_ptsum[3], 6, 1, 1);
34 book(_hist_ptsum[4], 5, 1, 1);
35 book(_hist_ptsum[5], 7, 1, 1);
36
37 // <pT> profiles vs pT_lead (not measured for trans diff)
38 book(_hist_ptavg[0], 29, 1, 1);
39 book(_hist_ptavg[1], 30, 1, 1);
40 book(_hist_ptavg[2], 11, 1, 1);
41 book(_hist_ptavg[3], 13, 1, 1);
42 book(_hist_ptavg[4], 12, 1, 1);
43
44 // <pT> profiles vs. Nch (not measured for trans diff)
45 book(_hist_dn_dpt[0], 27, 1, 1);
46 book(_hist_dn_dpt[1], 28, 1, 1);
47 book(_hist_dn_dpt[2], 8, 1, 1);
48 book(_hist_dn_dpt[3], 10, 1, 1);
49 book(_hist_dn_dpt[4], 9, 1, 1);
50
51 // Only measured for trans max/min
52 book(_hist_dn_dpt2[3], 32, 1, 1);
53 book(_hist_dn_dpt2[4], 31, 1, 1);
54
55 // Nch vs. Delta(phi) profiles
56 book(_hist_N_vs_dPhi[0], 15, 1, 1);
57 book(_hist_N_vs_dPhi[1], 16, 1, 1);
58 book(_hist_N_vs_dPhi[2], 17, 1, 1);
59
60 // pT vs. Delta(phi) profiles
61 book(_hist_pT_vs_dPhi[0], 18, 1, 1);
62 book(_hist_pT_vs_dPhi[1], 19, 1, 1);
63 book(_hist_pT_vs_dPhi[2], 20, 1, 1);
64
65 //ptLead histos only for 1 and 5 GeV cuts
66 book(_hist_ptLead[0], 1, 1, 1);
67 book(_hist_ptLead[1], 14, 1, 1);
68
69 for (size_t iC = 0; iC < NCUTS; ++iC) {
70 book(_counters[iC], "Ctr_cut_" + toString(iC));
71 }
72
73 }
74
75
76 void analyze(const Event& event) {
77
78 // Get charged particles (tracks) with pT > 500 MeV
79 const ChargedFinalState& charged500 = apply<ChargedFinalState>(event, "CFS500");
80 const Particles& particlesAll = charged500.particlesByPt();
81 MSG_DEBUG("Num tracks: " << particlesAll.size());
82
83 const Cut& pcut = ( (Cuts::abspid != PID::SIGMAMINUS) && (Cuts::abspid != PID::SIGMAPLUS) &&
84 (Cuts::abspid != PID::XIMINUS) && (Cuts::abspid != PID::OMEGAMINUS) );
85 const Particles& particles = charged500.particlesByPt(pcut);
86 MSG_DEBUG("Num tracks without strange baryons: " << particles.size());
87
88 // Require at least one track in the event for pTlead histograms
89 if (particles.empty()) vetoEvent;
90 for (size_t iC = 0; iC < 2; ++iC) {
91 if (particles[0].pT() < PTCUTS[iC]*GeV) continue;
92 _counters[iC]->fill();
93 _hist_ptLead[iC]->fill( particles[0].pT()/GeV);
94 }
95
96 // Require at least one track in the event with pT >= 1 GeV for the rest
97 if (particles[0].pT() < 1*GeV) vetoEvent;
98
99 // Identify leading track and its phi and pT
100 const Particle& p_lead = particles[0];
101 const double philead = p_lead.phi();
102 const double etalead = p_lead.eta();
103 const double pTlead = p_lead.perp();
104 MSG_DEBUG("Leading track: pT = " << pTlead << ", eta = " << etalead << ", phi = " << philead);
105
106 // Iterate over all particles and count particles and scalar pTsum in three basic regions
107 vector<double> num(NREGIONS, 0), ptSum(NREGIONS, 0.0), avgpt(NREGIONS, 0.0);
108
109 // Temporary histos that bin Nch and pT in dPhi.
110 Histo1D hist_num_dphi(_hist_N_vs_dPhi[0]->binning());
111 hist_num_dphi.setPath("/hist_num_dphi");
112 Histo1D hist_pt_dphi(_hist_pT_vs_dPhi[0]->binning());
113 hist_pt_dphi.setPath("/hist_pt_dphi");
114 hist_num_dphi.reset();
115 hist_pt_dphi.reset();
116
117 int tmpnch[2] = {0,0};
118 double tmpptsum[2] = {0,0};
119 for (const Particle& p : particles) {
120 const double pT = p.pT()/GeV;
121 const double dPhi = deltaPhi(philead, p.phi()); // in range (0,pi)
122 const int ir = region_index(dPhi); // gives just toward/away/trans
123
124 // Toward/away/trans region: just count
125 num [ir] += 1;
126 ptSum[ir] += pT;
127
128 // Determine which transverse side
129 if (ir == kTrans) {
130 const size_t iside = (mapAngleMPiToPi(p.phi() - philead) > 0) ? 0 : 1;
131 tmpnch [iside] += 1;
132 tmpptsum[iside] += p.pT();
133 }
134
135 // Fill temp histos to bin Nch and pT in dPhi
136 if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
137 hist_num_dphi.fill(dPhi/M_PI*180);
138 hist_pt_dphi .fill(dPhi/M_PI*180, pT/GeV);
139 }
140 }
141
142 // Construct max/min/diff regions
143 num[kTransMax ] = std::max(tmpnch[0], tmpnch[1]);
144 num[kTransMin ] = std::min(tmpnch[0], tmpnch[1]);
145 num[kTransDiff] = num[kTransMax ] - num[kTransMin ];
146 ptSum[kTransMax ] = std::max(tmpptsum[0], tmpptsum[1]);
147 ptSum[kTransMin ] = std::min(tmpptsum[0], tmpptsum[1]);
148 ptSum[kTransDiff] = ptSum[kTransMax ] - ptSum[kTransMin ];
149 avgpt[kToward] = (num[kToward] > 0 ) ? ptSum[kToward] / num[kToward] : 0. ;
150 avgpt[kAway] = (num[kAway ] > 0 ) ? ptSum[kAway] / num[kAway] : 0. ;
151 avgpt[kTrans] = (num[kTrans ] > 0 ) ? ptSum[kTrans] / num[kTrans] : 0. ;
152 // Avg pt max/min regions determined according sumpt max/min
153 int sumptMaxRegID = (tmpptsum[0] > tmpptsum[1]) ? 0 : 1 ;
154 int sumptMinRegID = (sumptMaxRegID == 0) ? 1 : 0;
155 avgpt[kTransMax ] = (tmpnch[sumptMaxRegID] > 0) ? tmpptsum[sumptMaxRegID] / tmpnch[sumptMaxRegID] : 0.;
156 avgpt[kTransMin ] = (tmpnch[sumptMinRegID] > 0) ? tmpptsum[sumptMinRegID] / tmpnch[sumptMinRegID] : 0.;
157 avgpt[kTransDiff] = ((tmpnch[sumptMaxRegID] > 0) && (tmpnch[sumptMinRegID] > 0)) ? avgpt[kTransMax ] - avgpt[kTransMin ] : 0.;
158
159
160 // Now fill underlying event histograms
161
162 // The densities are calculated by dividing the UE properties by dEta*dPhi
163 // -- each basic region has a dPhi of 2*PI/3 and dEta is two times 2.5
164 // min/max/diff regions are only half of that
165 const double dEtadPhi[NREGIONS] = { 2*2.5 * 2*PI/3.0, 2*2.5 * 2*PI/3.0, 2*2.5 * 2*PI/3.0,
166 2*2.5 * PI/3.0, 2*2.5 * PI/3.0, 2*2.5 * PI/3.0 };
167 for (size_t iR = 0; iR < NREGIONS; ++iR) {
168
169 _hist_nch [iR]->fill(pTlead/GeV, num[iR] /dEtadPhi[iR] );
170 _hist_ptsum[iR]->fill(pTlead/GeV, ptSum[iR] /GeV/dEtadPhi[iR] );
171
172 // <pT> profiles vs. pT_lead (first 3 are the same!)
173 switch (iR) {
174 case kToward :
175 case kAway :
176 case kTrans :
177 if (num[iR] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
178 break;
179 case kTransMax :
180 if (tmpnch[sumptMaxRegID] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
181 break;
182 case kTransMin :
183 if (tmpnch[sumptMinRegID] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
184 break;
185 case kTransDiff :
186 break;
187 default: //should not get here!!!
188 MSG_WARNING("Unknown region in <pT> profiles vs.pt lead switch!!! : " << iR);
189 }
190
191 // <pT> profiles vs. Nch (first 3 are the same!)
192 switch (iR) {
193 case kToward :
194 case kAway :
195 case kTrans :
196 if (num[iR] > 0) _hist_dn_dpt[iR]->fill(num[iR] , avgpt[iR]/GeV);
197 break;
198 case kTransMax :
199 if (tmpnch[sumptMaxRegID] > 0) {
200 _hist_dn_dpt [iR]->fill(num[kTrans] , avgpt[iR]/GeV);
201 _hist_dn_dpt2[iR]->fill(tmpnch[sumptMaxRegID], avgpt[iR]/GeV);
202 }
203 break;
204 case kTransMin :
205 if (tmpnch[sumptMinRegID] > 0) {
206 _hist_dn_dpt [iR]->fill(num[kTrans] , avgpt[iR]/GeV);
207 _hist_dn_dpt2[iR]->fill(tmpnch[sumptMinRegID], avgpt[iR]/GeV);
208 }
209 break;
210 case kTransDiff :
211 break;
212 default : //should not get here!!!
213 MSG_INFO("unknown region in <pT> profiles vs. nch switch!!! : " << iR);
214 }
215
216 }
217
218
219 // Update the "proper" dphi profile histograms
220 // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
221 const double dEtadPhi2 = (2*2.5 * 2) * M_PI/180.;
222 for (size_t i = 0; i < hist_num_dphi.numBins(); ++i) {
223
224 // First Nch
225 double mean = hist_num_dphi.bin(i).xMid() ;
226 double value = 0.;
227 if (hist_num_dphi.bin(i).numEntries() > 0) {
228 mean = hist_num_dphi.bin(i).xMean() ;
229 value = hist_num_dphi.bin(i).sumW()/hist_num_dphi.bin(i).xWidth()/dEtadPhi2;
230 }
231 for (size_t iC = 0; iC < NCUTS; ++iC) {
232 if (pTlead >= PTCUTS[iC]*GeV) _hist_N_vs_dPhi[iC] ->fill(mean, value);
233 }
234
235 // Then pT
236 mean = hist_pt_dphi.bin(i).xMid();
237 value = 0.;
238 if (hist_pt_dphi.bin(i).numEntries() > 0) {
239 mean = hist_pt_dphi.bin(i).xMean() ;
240 value = hist_pt_dphi.bin(i).sumW()/hist_pt_dphi.bin(i).xWidth()/dEtadPhi2;
241 }
242 for (size_t iC = 0; iC < NCUTS; ++iC) {
243 if (pTlead >= PTCUTS[iC]*GeV) _hist_pT_vs_dPhi[iC] ->fill(mean, value);
244 }
245 }
246
247 }
248
249
250 void finalize() {
251 for (size_t iC = 0; iC < NCUTS; ++iC) {
252 if (iC == 0 || iC == 1) scale(_hist_ptLead[iC], 1.0/_counters[iC]->sumW());
253 }
254 }
255
256
257 private:
258
259 enum regionID {
260 kToward = 0,
261 kAway,
262 kTrans,
263 kTransMax,
264 kTransMin,
265 kTransDiff,
266 NREGIONS
267 };
268
269 // Little helper function to identify basic Delta(phi) regions: toward/away/trans
270 int region_index(double dphi) {
271 assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
272 if (dphi < PI/3.0) return kToward;
273 if (dphi < 2*PI/3.0) return kTrans;
274 return kAway;
275 }
276
277 const static size_t NCUTS = 3;
278 const vector<double> PTCUTS = { 1., 5., 10. };
279
280
281 /// @name Histograms
282 /// @{
283
284 // Nch, sumpT, avgpT profiles vs. pTlead
285 Profile1DPtr _hist_nch [NREGIONS]; //for regions: all 6 regions
286 Profile1DPtr _hist_ptsum [NREGIONS]; //for regions: all 6 regions
287 Profile1DPtr _hist_ptavg [NREGIONS]; //for regions: trans towards/away/all/min/max
288 // Nch, sumpT, avgpT profiles vs. Nch
289 Profile1DPtr _hist_dn_dpt [NREGIONS]; //regions: towards/away/ vs nch(region) & trans all/min/max vs nch(trans)
290 Profile1DPtr _hist_dn_dpt2[NREGIONS]; //regions: trans min/max vs. nch(region)
291
292 Profile1DPtr _hist_N_vs_dPhi [NCUTS];
293 Profile1DPtr _hist_pT_vs_dPhi[NCUTS];
294 Histo1DPtr _hist_ptLead[NCUTS]; //for 1,5 GeV cuts only
295 CounterPtr _counters[NCUTS];
296
297 /// @}
298
299 };
300
301
302 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1509919);
303
304}
|