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], "/hist_num_dphi");
111 Histo1D hist_pt_dphi(*_hist_pT_vs_dPhi[0], "/hist_pt_dphi");
112 hist_num_dphi.reset();
113 hist_pt_dphi .reset();
114
115 int tmpnch[2] = {0,0};
116 double tmpptsum[2] = {0,0};
117 for (const Particle& p : particles) {
118 const double pT = p.pT()/GeV;
119 const double dPhi = deltaPhi(philead, p.phi()); // in range (0,pi)
120 const int ir = region_index(dPhi); // gives just toward/away/trans
121
122 // Toward/away/trans region: just count
123 num [ir] += 1;
124 ptSum[ir] += pT;
125
126 // Determine which transverse side
127 if (ir == kTrans) {
128 const size_t iside = (mapAngleMPiToPi(p.phi() - philead) > 0) ? 0 : 1;
129 tmpnch [iside] += 1;
130 tmpptsum[iside] += p.pT();
131 }
132
133 // Fill temp histos to bin Nch and pT in dPhi
134 if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
135 hist_num_dphi.fill(dPhi/M_PI*180);
136 hist_pt_dphi .fill(dPhi/M_PI*180, pT/GeV);
137 }
138 }
139
140 // Construct max/min/diff regions
141 num[kTransMax ] = std::max(tmpnch[0], tmpnch[1]);
142 num[kTransMin ] = std::min(tmpnch[0], tmpnch[1]);
143 num[kTransDiff] = num[kTransMax ] - num[kTransMin ];
144 ptSum[kTransMax ] = std::max(tmpptsum[0], tmpptsum[1]);
145 ptSum[kTransMin ] = std::min(tmpptsum[0], tmpptsum[1]);
146 ptSum[kTransDiff] = ptSum[kTransMax ] - ptSum[kTransMin ];
147 avgpt[kToward] = (num[kToward] > 0 ) ? ptSum[kToward] / num[kToward] : 0. ;
148 avgpt[kAway] = (num[kAway ] > 0 ) ? ptSum[kAway] / num[kAway] : 0. ;
149 avgpt[kTrans] = (num[kTrans ] > 0 ) ? ptSum[kTrans] / num[kTrans] : 0. ;
150 // Avg pt max/min regions determined according sumpt max/min
151 int sumptMaxRegID = (tmpptsum[0] > tmpptsum[1]) ? 0 : 1 ;
152 int sumptMinRegID = (sumptMaxRegID == 0) ? 1 : 0;
153 avgpt[kTransMax ] = (tmpnch[sumptMaxRegID] > 0) ? tmpptsum[sumptMaxRegID] / tmpnch[sumptMaxRegID] : 0.;
154 avgpt[kTransMin ] = (tmpnch[sumptMinRegID] > 0) ? tmpptsum[sumptMinRegID] / tmpnch[sumptMinRegID] : 0.;
155 avgpt[kTransDiff] = ((tmpnch[sumptMaxRegID] > 0) && (tmpnch[sumptMinRegID] > 0)) ? avgpt[kTransMax ] - avgpt[kTransMin ] : 0.;
156
157
158 // Now fill underlying event histograms
159
160 // The densities are calculated by dividing the UE properties by dEta*dPhi
161 // -- each basic region has a dPhi of 2*PI/3 and dEta is two times 2.5
162 // min/max/diff regions are only half of that
163 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,
164 2*2.5 * PI/3.0, 2*2.5 * PI/3.0, 2*2.5 * PI/3.0 };
165 for (size_t iR = 0; iR < NREGIONS; ++iR) {
166
167 _hist_nch [iR]->fill(pTlead/GeV, num[iR] /dEtadPhi[iR] );
168 _hist_ptsum[iR]->fill(pTlead/GeV, ptSum[iR] /GeV/dEtadPhi[iR] );
169
170 // <pT> profiles vs. pT_lead (first 3 are the same!)
171 switch (iR) {
172 case kToward :
173 case kAway :
174 case kTrans :
175 if (num[iR] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
176 break;
177 case kTransMax :
178 if (tmpnch[sumptMaxRegID] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
179 break;
180 case kTransMin :
181 if (tmpnch[sumptMinRegID] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
182 break;
183 case kTransDiff :
184 break;
185 default: //should not get here!!!
186 MSG_WARNING("Unknown region in <pT> profiles vs.pt lead switch!!! : " << iR);
187 }
188
189 // <pT> profiles vs. Nch (first 3 are the same!)
190 switch (iR) {
191 case kToward :
192 case kAway :
193 case kTrans :
194 if (num[iR] > 0) _hist_dn_dpt[iR]->fill(num[iR] , avgpt[iR]/GeV);
195 break;
196 case kTransMax :
197 if (tmpnch[sumptMaxRegID] > 0) {
198 _hist_dn_dpt [iR]->fill(num[kTrans] , avgpt[iR]/GeV);
199 _hist_dn_dpt2[iR]->fill(tmpnch[sumptMaxRegID], avgpt[iR]/GeV);
200 }
201 break;
202 case kTransMin :
203 if (tmpnch[sumptMinRegID] > 0) {
204 _hist_dn_dpt [iR]->fill(num[kTrans] , avgpt[iR]/GeV);
205 _hist_dn_dpt2[iR]->fill(tmpnch[sumptMinRegID], avgpt[iR]/GeV);
206 }
207 break;
208 case kTransDiff :
209 break;
210 default : //should not get here!!!
211 MSG_INFO("unknown region in <pT> profiles vs. nch switch!!! : " << iR);
212 }
213
214 }
215
216
217 // Update the "proper" dphi profile histograms
218 // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
219 const double dEtadPhi2 = (2*2.5 * 2) * M_PI/180.;
220 for (size_t i = 0; i < hist_num_dphi.numBins(); ++i) {
221
222 // First Nch
223 double mean = hist_num_dphi.bin(i).xMid() ;
224 double value = 0.;
225 if (hist_num_dphi.bin(i).numEntries() > 0) {
226 mean = hist_num_dphi.bin(i).xMean() ;
227 value = hist_num_dphi.bin(i).area()/hist_num_dphi.bin(i).xWidth()/dEtadPhi2;
228 }
229 for (size_t iC = 0; iC < NCUTS; ++iC) {
230 if (pTlead >= PTCUTS[iC]*GeV) _hist_N_vs_dPhi[iC] ->fill(mean, value);
231 }
232
233 // Then pT
234 mean = hist_pt_dphi.bin(i).xMid();
235 value = 0.;
236 if (hist_pt_dphi.bin(i).numEntries() > 0) {
237 mean = hist_pt_dphi.bin(i).xMean() ;
238 value = hist_pt_dphi.bin(i).area()/hist_pt_dphi.bin(i).xWidth()/dEtadPhi2;
239 }
240 for (size_t iC = 0; iC < NCUTS; ++iC) {
241 if (pTlead >= PTCUTS[iC]*GeV) _hist_pT_vs_dPhi[iC] ->fill(mean, value);
242 }
243 }
244
245 }
246
247
248 void finalize() {
249 for (size_t iC = 0; iC < NCUTS; ++iC) {
250 if (iC == 0 || iC == 1) scale(_hist_ptLead[iC], 1.0/_counters[iC]->sumW());
251 }
252 }
253
254
255 private:
256
257 enum regionID {
258 kToward = 0,
259 kAway,
260 kTrans,
261 kTransMax,
262 kTransMin,
263 kTransDiff,
264 NREGIONS
265 };
266
267 // Little helper function to identify basic Delta(phi) regions: toward/away/trans
268 int region_index(double dphi) {
269 assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
270 if (dphi < PI/3.0) return kToward;
271 if (dphi < 2*PI/3.0) return kTrans;
272 return kAway;
273 }
274
275 const static size_t NCUTS = 3;
276 const vector<double> PTCUTS = { 1., 5., 10. };
277
278
279 /// @name Histograms
280 //@{
281
282 // Nch, sumpT, avgpT profiles vs. pTlead
283 Profile1DPtr _hist_nch [NREGIONS]; //for regions: all 6 regions
284 Profile1DPtr _hist_ptsum [NREGIONS]; //for regions: all 6 regions
285 Profile1DPtr _hist_ptavg [NREGIONS]; //for regions: trans towards/away/all/min/max
286 // Nch, sumpT, avgpT profiles vs. Nch
287 Profile1DPtr _hist_dn_dpt [NREGIONS]; //regions: towards/away/ vs nch(region) & trans all/min/max vs nch(trans)
288 Profile1DPtr _hist_dn_dpt2[NREGIONS]; //regions: trans min/max vs. nch(region)
289
290 Profile1DPtr _hist_N_vs_dPhi [NCUTS];
291 Profile1DPtr _hist_pT_vs_dPhi[NCUTS];
292 Histo1DPtr _hist_ptLead[NCUTS]; //for 1,5 GeV cuts only
293 CounterPtr _counters[NCUTS];
294
295 //@}
296
297 };
298
299
300 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1509919);
301
302}
|