Rivet analyses referenceATLAS_2010_I879407Track-based underlying event at 900 GeV and 7 TeV in ATLASExperiment: ATLAS (LHC) Inspire ID: 879407 Status: VALIDATED Authors:
Beam energies: (450.0, 450.0); (3500.0, 3500.0) GeV Run details:
The underlying event measurements with the ATLAS detector at the LHC at the center of mass energies of 900 GeV and 7 TeV. The observables sensitive to the underlying event, i.e the charged particle density and charged pT sum, as well as their standard deviations and the average pT, are measured as functions of the leading track. A track pT cut of 500 MeV is applied for most observables, but the main profile plots are also shown for a lower track cut of 100 MeV, which includes much more of the soft cross-section. The angular distribution of the charged tracks with respect to the leading track is also studied, as are the correlation between mean transverse momentum and charged particle multiplicity, and the `plateau' height as a function of the leading track $|\eta|$. Beam energy must be specified (in GeV) as analysis option "ENERGY" when rivet-merging samples.' Source code: ATLAS_2010_I879407.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4
5namespace Rivet {
6
7
8 class ATLAS_2010_I879407 : public Analysis {
9 public:
10
11 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2010_I879407);
12
13
14 void init() {
15 const ChargedFinalState cfs100((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >= 100*MeV));
16 declare(cfs100, "CFS100");
17 const ChargedFinalState cfs500((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >= 500*MeV));
18 declare(cfs500, "CFS500");
19 const ChargedFinalState cfslead((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >= 1.0*GeV));
20 declare(cfslead, "CFSlead");
21
22 // Get an index for the beam energy
23 int isqrts = -1;
24 if (isCompatibleWithSqrtS(900*GeV)) isqrts = 0;
25 else if (isCompatibleWithSqrtS(7000*GeV)) isqrts = 1;
26 assert(isqrts >= 0);
27
28 // Nch profiles, 500 MeV track pT cut
29 book(_hist_nch_transverse_500[0] ,1+isqrts, 1, 1);
30 book(_hist_nch_toward_500 ,1+isqrts, 1, 2);
31 book(_hist_nch_away_500 ,1+isqrts, 1, 3);
32
33 // pTsum profiles, 500 MeV track pT cut
34 book(_hist_ptsum_transverse_500[0] ,3+isqrts, 1, 1);
35 book(_hist_ptsum_toward_500 ,3+isqrts, 1, 2);
36 book(_hist_ptsum_away_500 ,3+isqrts, 1, 3);
37
38 // Standard deviation profiles
39 // First the higher moments of main profiles to calculate variance and error on variance...
40 for (size_t i = 1; i < 4; ++i) {
41 book(_hist_nch_transverse_500[i], "TMP/nch"+to_str(i), refData(1+isqrts, 1, 1));
42 book(_hist_ptsum_transverse_500[i], "TMP/ptsum"+to_str(i), refData(3+isqrts, 1, 1));
43 }
44 // Then the data point sets into which the results will be inserted
45 book(_dps_sdnch_transverse_500 , 5+isqrts, 1, 1);
46 book(_dps_sdptsum_transverse_500, 7+isqrts, 1, 1);
47
48 // <pT> profiles, 500 MeV track pT cut
49 book(_hist_ptavg_transverse_500 ,9+isqrts, 1, 1);
50 book(_hist_ptavg_toward_500 ,9+isqrts, 1, 2);
51 book(_hist_ptavg_away_500 ,9+isqrts, 1, 3);
52
53 // <pT> vs. Nch profiles, 500 MeV track pT cut
54 book(_hist_dn_dpt_transverse_500 ,11+isqrts, 1, 1);
55 book(_hist_dn_dpt_toward_500 ,11+isqrts, 1, 2);
56 book(_hist_dn_dpt_away_500 ,11+isqrts, 1, 3);
57
58 // Nch vs. Delta(phi) profiles, 500 MeV track pT cut
59 book(_hist_N_vs_dPhi_1_500 ,13+isqrts, 1, 1);
60 book(_hist_N_vs_dPhi_2_500 ,13+isqrts, 1, 2);
61 book(_hist_N_vs_dPhi_3_500 ,13+isqrts, 1, 3);
62 book(_hist_N_vs_dPhi_5_500 ,13+isqrts, 1, 4);
63 // pT vs. Delta(phi) profiles, 500 MeV track pT cut
64 book(_hist_pT_vs_dPhi_1_500 ,15+isqrts, 1, 1);
65 book(_hist_pT_vs_dPhi_2_500 ,15+isqrts, 1, 2);
66 book(_hist_pT_vs_dPhi_3_500 ,15+isqrts, 1, 3);
67 book(_hist_pT_vs_dPhi_5_500 ,15+isqrts, 1, 4);
68
69 // Nch and pTsum profiles, 100 MeV track pT cut
70 book(_hist_nch_transverse_100 ,17+isqrts, 1, 1);
71 book(_hist_nch_toward_100 ,17+isqrts, 1, 2);
72 book(_hist_nch_away_100 ,17+isqrts, 1, 3);
73 book(_hist_ptsum_transverse_100 ,19+isqrts, 1, 1);
74 book(_hist_ptsum_toward_100 ,19+isqrts, 1, 2);
75 book(_hist_ptsum_away_100 ,19+isqrts, 1, 3);
76
77 // Profiles vs. eta (7 TeV only)
78 if (isqrts == 1) {
79 book(_hist_nch_vs_eta_transverse_100 ,21, 1, 1);
80 book(_hist_ptsum_vs_eta_transverse_100 ,22, 1, 1);
81 }
82
83 }
84
85
86 // Little helper function to identify Delta(phi) regions
87 inline int region_index(double dphi) {
88 assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
89 if (dphi < PI/3.0) return 0;
90 if (dphi < 2*PI/3.0) return 1;
91 return 2;
92 }
93
94
95 void analyze(const Event& event) {
96 // Require at least one track in the event with pT >= 1 GeV
97 const ChargedFinalState& cfslead = apply<ChargedFinalState>(event, "CFSlead");
98 if (cfslead.size() < 1) {
99 vetoEvent;
100 }
101
102 // These are the charged particles (tracks) with pT > 500 MeV
103 const ChargedFinalState& charged500 = apply<ChargedFinalState>(event, "CFS500");
104
105 // Identify leading track and its phi and pT (this is the same for both the 100 MeV and 500 MeV track cuts)
106 Particles particles500 = charged500.particlesByPt();
107 Particle p_lead = particles500[0];
108 const double philead = p_lead.phi();
109 const double etalead = p_lead.eta();
110 const double pTlead = p_lead.perp();
111 MSG_DEBUG("Leading track: pT = " << pTlead << ", eta = " << etalead << ", phi = " << philead);
112
113 // Iterate over all > 500 MeV particles and count particles and scalar pTsum in the three regions
114 vector<double> num500(3, 0), ptSum500(3, 0.0);
115 // Temporary histos that bin Nch and pT in dPhi.
116 // NB. Only one of each needed since binnings are the same for the energies and pT cuts
117 Histo1D hist_num_dphi_500(refData(13,1,1).binning());
118 Histo1D hist_pt_dphi_500(refData(15,1,1).binning());
119 for (const Particle& p : particles500) {
120 const double pT = p.pT();
121 const double dPhi = deltaPhi(philead, p.phi());
122 const int ir = region_index(dPhi);
123 num500[ir] += 1;
124 ptSum500[ir] += pT;
125
126 // Fill temp histos to bin Nch and pT in dPhi
127 if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
128 hist_num_dphi_500.fill(dPhi, 1);
129 hist_pt_dphi_500.fill(dPhi, pT);
130 }
131 }
132
133
134 // Iterate over charged particles again for profiles against Nch
135 // This is necessary since the Nch are region-specific and so are only known after the first loop
136 for (const Particle& p : particles500) {
137 const double pT = p.pT();
138 const double dPhi = deltaPhi(philead, p.phi());
139 const int ir = region_index(dPhi);
140 switch (ir) {
141 case 0:
142 _hist_dn_dpt_toward_500->fill(num500[0], pT);
143 break;
144 case 1:
145 _hist_dn_dpt_transverse_500->fill(num500[1], pT);
146 break;
147 case 2:
148 _hist_dn_dpt_away_500->fill(num500[2], pT);
149 break;
150 default:
151 assert(false && "How did we get here?");
152 }
153 }
154
155
156 // Now fill underlying event histograms
157 // The densities are calculated by dividing the UE properties by dEta*dPhi
158 // -- each region has a dPhi of 2*PI/3 and dEta is two times 2.5
159 const double dEtadPhi = (2*2.5 * 2*PI/3.0);
160 // Transverse profiles need 4 orders of moments for stddev with errors
161 for (int i = 0; i < 4; ++i) {
162 _hist_nch_transverse_500[i]->fill(pTlead/GeV, intpow(num500[1]/dEtadPhi, i+1));
163 _hist_ptsum_transverse_500[i]->fill(pTlead/GeV, intpow(ptSum500[1]/GeV/dEtadPhi, i+1));
164 }
165 // Toward and away profiles only need the first moment (the mean)
166 _hist_nch_toward_500->fill(pTlead/GeV, num500[0]/dEtadPhi);
167 _hist_nch_away_500->fill(pTlead/GeV, num500[2]/dEtadPhi);
168 _hist_ptsum_toward_500->fill(pTlead/GeV, ptSum500[0]/GeV/dEtadPhi);
169 _hist_ptsum_away_500->fill(pTlead/GeV, ptSum500[2]/GeV/dEtadPhi);
170 // <pT> profiles
171 //MSG_INFO("Trans pT1, pTsum, Nch, <pT>" << pTlead/GeV << ", " << ptSum500[1]/GeV << ", " << num500[1] << ", " << ptSum500[1]/GeV/num500[1]);
172 if (num500[1] > 0) _hist_ptavg_transverse_500->fill(pTlead/GeV, ptSum500[1]/GeV/num500[1]);
173 if (num500[0] > 0) _hist_ptavg_toward_500->fill(pTlead/GeV, ptSum500[0]/GeV/num500[0]);
174 if (num500[2] > 0) _hist_ptavg_away_500->fill(pTlead/GeV, ptSum500[2]/GeV/num500[2]);
175
176
177 // Update the "proper" dphi profile histograms
178 // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
179 // The values tabulated in the note are for an (undefined) signed Delta(phi) rather than
180 // |Delta(phi)| and so differ by a factor of 2: we have to actually norm for angular range = 2pi
181 const size_t nbins = refData(13,1,1).numBins();
182 std::vector<double> ptcut;
183 if (isCompatibleWithSqrtS(900*GeV)) {
184 ptcut += 1.0; ptcut += 1.5; ptcut += 2.0; ptcut += 2.5;
185 }
186 else if (isCompatibleWithSqrtS(7000*GeV)) {
187 ptcut += 1.0; ptcut += 2.0; ptcut += 3.0; ptcut += 5.0;
188 }
189 assert(ptcut.size() == 4);
190 for (size_t i = 0; i < nbins; ++i) {
191 // First Nch
192 double mean = hist_num_dphi_500.bin(i).xMid();
193 double value = 0.;
194 if (hist_num_dphi_500.bin(i).numEntries() > 0) {
195 mean = hist_num_dphi_500.bin(i).xMean();
196 value = hist_num_dphi_500.bin(i).sumW()/hist_num_dphi_500.bin(i).xWidth()/10.0;
197 }
198 if (pTlead/GeV >= ptcut[0]) _hist_N_vs_dPhi_1_500->fill(mean, value);
199 if (pTlead/GeV >= ptcut[1]) _hist_N_vs_dPhi_2_500->fill(mean, value);
200 if (pTlead/GeV >= ptcut[2]) _hist_N_vs_dPhi_3_500->fill(mean, value);
201 if (pTlead/GeV >= ptcut[3]) _hist_N_vs_dPhi_5_500->fill(mean, value);
202
203 // Then pT
204 mean = hist_pt_dphi_500.bin(i).xMid();
205 value = 0.;
206 if (hist_pt_dphi_500.bin(i).numEntries() > 0) {
207 mean = hist_pt_dphi_500.bin(i).xMean();
208 value = hist_pt_dphi_500.bin(i).sumW()/hist_pt_dphi_500.bin(i).xWidth()/10.0;
209 }
210 if (pTlead/GeV >= ptcut[0]) _hist_pT_vs_dPhi_1_500->fill(mean, value);
211 if (pTlead/GeV >= ptcut[1]) _hist_pT_vs_dPhi_2_500->fill(mean, value);
212 if (pTlead/GeV >= ptcut[2]) _hist_pT_vs_dPhi_3_500->fill(mean, value);
213 if (pTlead/GeV >= ptcut[3]) _hist_pT_vs_dPhi_5_500->fill(mean, value);
214 }
215
216
217 //////////////////////
218
219
220 // These are the charged particles (tracks) with pT > 100 MeV
221 const ChargedFinalState& charged100 = apply<ChargedFinalState>(event, "CFS100");
222
223 // Iterate over all > 100 MeV particles and count particles and scalar pTsum in the three regions
224 vector<double> num100(3, 0), ptSum100(3, 0.0);
225 for (const Particle& p : charged100.particles()) {
226 const double pT = p.pT();
227 const double dPhi = deltaPhi(philead, p.phi());
228 const int ir = region_index(dPhi);
229 num100[ir] += 1;
230 ptSum100[ir] += pT;
231 }
232
233 // Now fill the two sets of 100 MeV underlying event histograms
234 _hist_nch_transverse_100->fill(pTlead/GeV, num100[1]/dEtadPhi);
235 _hist_nch_toward_100->fill(pTlead/GeV, num100[0]/dEtadPhi);
236 _hist_nch_away_100->fill(pTlead/GeV, num100[2]/dEtadPhi);
237 _hist_ptsum_transverse_100->fill(pTlead/GeV, ptSum100[1]/GeV/dEtadPhi);
238 _hist_ptsum_toward_100->fill(pTlead/GeV, ptSum100[0]/GeV/dEtadPhi);
239 _hist_ptsum_away_100->fill(pTlead/GeV, ptSum100[2]/GeV/dEtadPhi);
240
241 // And finally the Nch and pT vs eta_lead profiles (again from > 100 MeV tracks, and only at 7 TeV)
242 if (isCompatibleWithSqrtS(7000*GeV) && pTlead > 5*GeV) {
243 _hist_nch_vs_eta_transverse_100->fill(etalead, num100[1]/dEtadPhi);
244 _hist_ptsum_vs_eta_transverse_100->fill(etalead, ptSum100[1]/GeV/dEtadPhi);
245 }
246
247 }
248
249
250 void finalize() {
251 // Convert the various moments of the 500 MeV trans pT and Nch distributions to std devs with correct error
252 _moments_to_stddev(_hist_nch_transverse_500, _dps_sdnch_transverse_500);
253 _moments_to_stddev(_hist_ptsum_transverse_500, _dps_sdptsum_transverse_500);
254 }
255
256
257 private:
258
259
260 inline void _moments_to_stddev(Profile1DPtr moment_profiles[], Estimate1DPtr target_dps) {
261 for (size_t b = 1; b < moment_profiles[0]->numBins()+1; ++b) { // loop over points
262 /// @todo Assuming unit weights here! Should use N_effective = sumW**2/sumW2?
263 const double numentries = moment_profiles[0]->bin(b).numEntries();
264 double var = 0., sd = 0.;
265 if (numentries > 0) {
266 var = moment_profiles[1]->bin(b).yMean() - intpow(moment_profiles[0]->bin(b).yMean(), 2);
267 sd = fuzzyLessEquals(var,0.) ? 0 : sqrt(var); ///< Numerical safety check
268 }
269 if (sd == 0 || numentries < 3) {
270 MSG_WARNING("Need at least 3 bin entries and a non-zero central value to calculate "
271 << "an error on standard deviation profiles (bin " << b << ")");
272 target_dps->bin(b).set(sd, 0);
273 continue;
274 }
275 // c2(y) = m4(x) - 4 m3(x) m1(x) - m2(x)^2 + 8 m2(x) m1(x)^2 - 4 m1(x)^4
276 const double var_on_var = moment_profiles[3]->bin(b).yMean()
277 - 4 * moment_profiles[2]->bin(b).yMean() * moment_profiles[0]->bin(b).yMean()
278 - intpow(moment_profiles[1]->bin(b).yMean(), 2)
279 + 8 * moment_profiles[1]->bin(b).yMean() * intpow(moment_profiles[0]->bin(b).yMean(), 2)
280 - 4 * intpow(moment_profiles[0]->bin(b).yMean(), 4);
281 const double stderr_on_var = sqrt(var_on_var/(numentries-2.0));
282 const double stderr_on_sd = stderr_on_var / (2.0*sd);
283 target_dps->bin(b).set(sd, stderr_on_sd);
284 }
285 }
286
287
288 private:
289
290 Profile1DPtr _hist_nch_transverse_500[4];
291 Profile1DPtr _hist_nch_toward_500;
292 Profile1DPtr _hist_nch_away_500;
293
294 Profile1DPtr _hist_ptsum_transverse_500[4];
295 Profile1DPtr _hist_ptsum_toward_500;
296 Profile1DPtr _hist_ptsum_away_500;
297
298 Estimate1DPtr _dps_sdnch_transverse_500;
299 Estimate1DPtr _dps_sdptsum_transverse_500;
300
301 Profile1DPtr _hist_ptavg_transverse_500;
302 Profile1DPtr _hist_ptavg_toward_500;
303 Profile1DPtr _hist_ptavg_away_500;
304
305 Profile1DPtr _hist_dn_dpt_transverse_500;
306 Profile1DPtr _hist_dn_dpt_toward_500;
307 Profile1DPtr _hist_dn_dpt_away_500;
308
309 Profile1DPtr _hist_N_vs_dPhi_1_500;
310 Profile1DPtr _hist_N_vs_dPhi_2_500;
311 Profile1DPtr _hist_N_vs_dPhi_3_500;
312 Profile1DPtr _hist_N_vs_dPhi_5_500;
313
314 Profile1DPtr _hist_pT_vs_dPhi_1_500;
315 Profile1DPtr _hist_pT_vs_dPhi_2_500;
316 Profile1DPtr _hist_pT_vs_dPhi_3_500;
317 Profile1DPtr _hist_pT_vs_dPhi_5_500;
318
319 Profile1DPtr _hist_nch_transverse_100;
320 Profile1DPtr _hist_nch_toward_100;
321 Profile1DPtr _hist_nch_away_100;
322 Profile1DPtr _hist_ptsum_transverse_100;
323 Profile1DPtr _hist_ptsum_toward_100;
324 Profile1DPtr _hist_ptsum_away_100;
325
326 Profile1DPtr _hist_nch_vs_eta_transverse_100;
327 Profile1DPtr _hist_ptsum_vs_eta_transverse_100;
328
329 };
330
331
332
333 RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2010_I879407, ATLAS_2010_S8894728);
334
335}
|