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|$.' 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::abseta < 2.5 && Cuts::pT >= 100*MeV);
16 declare(cfs100, "CFS100");
17 const ChargedFinalState cfs500(Cuts::abseta < 2.5 && Cuts::pT >= 500*MeV);
18 declare(cfs500, "CFS500");
19 const ChargedFinalState cfslead(Cuts::abseta < 2.5 && Cuts::pT >= 1.0*GeV);
20 declare(cfslead, "CFSlead");
21
22 for (double eVal : allowedEnergies()) {
23
24 const string en = toString(int(eVal));
25 if (isCompatibleWithSqrtS(eVal)) _sqs = en;
26 bool ih(en == "7000"s);
27
28 // Nch profiles, 500 MeV track pT cut
29 book(_p[en+"nch_transverse_500[0]"], 1+ih, 1, 1);
30 book(_p[en+"nch_toward_500"], 1+ih, 1, 2);
31 book(_p[en+"nch_away_500"], 1+ih, 1, 3);
32
33 // pTsum profiles, 500 MeV track pT cut
34 book(_p[en+"ptsum_transverse_500[0]"], 3+ih, 1, 1);
35 book(_p[en+"ptsum_toward_500"], 3+ih, 1, 2);
36 book(_p[en+"ptsum_away_500"], 3+ih, 1, 3);
37
38 // Standard deviation profiles
39 // First the higher moments of main profiles
40 // to calculate variance and error on variance...
41 for (size_t i = 1; i < 4; ++i) {
42 const string mom = toString(i);
43 book(_p[en+"nch_transverse_500_"+mom], "TMP/nch"+mom, refData(1+ih, 1, 1));
44 book(_p[en+"ptsum_transverse_500_"+mom], "TMP/ptsum"+mom, refData(3+ih, 1, 1));
45 }
46 // Then the data point sets into which the results will be inserted
47 book(_e[en+"nch_transverse_500"] , 5+ih, 1, 1);
48 book(_e[en+"ptsum_transverse_500"], 7+ih, 1, 1);
49
50 // <pT> profiles, 500 MeV track pT cut
51 book(_p[en+"ptavg_transverse_500"], 9+ih, 1, 1);
52 book(_p[en+"ptavg_toward_500"], 9+ih, 1, 2);
53 book(_p[en+"ptavg_away_500"], 9+ih, 1, 3);
54
55 // <pT> vs. Nch profiles, 500 MeV track pT cut
56 book(_p[en+"dn_dpt_transverse_500"], 11+ih, 1, 1);
57 book(_p[en+"dn_dpt_toward_500"], 11+ih, 1, 2);
58 book(_p[en+"dn_dpt_away_500"], 11+ih, 1, 3);
59
60 // Nch vs. Delta(phi) profiles, 500 MeV track pT cut
61 book(_p[en+"N_vs_dPhi_1_500"], 13+ih, 1, 1);
62 book(_p[en+"N_vs_dPhi_2_500"], 13+ih, 1, 2);
63 book(_p[en+"N_vs_dPhi_3_500"], 13+ih, 1, 3);
64 book(_p[en+"N_vs_dPhi_5_500"], 13+ih, 1, 4);
65 // pT vs. Delta(phi) profiles, 500 MeV track pT cut
66 book(_p[en+"pT_vs_dPhi_1_500"], 15+ih, 1, 1);
67 book(_p[en+"pT_vs_dPhi_2_500"], 15+ih, 1, 2);
68 book(_p[en+"pT_vs_dPhi_3_500"], 15+ih, 1, 3);
69 book(_p[en+"pT_vs_dPhi_5_500"], 15+ih, 1, 4);
70
71 // Nch and pTsum profiles, 100 MeV track pT cut
72 book(_p[en+"nch_transverse_100"], 17+ih, 1, 1);
73 book(_p[en+"nch_toward_100"], 17+ih, 1, 2);
74 book(_p[en+"nch_away_100"], 17+ih, 1, 3);
75 book(_p[en+"ptsum_transverse_100"], 19+ih, 1, 1);
76 book(_p[en+"ptsum_toward_100"], 19+ih, 1, 2);
77 book(_p[en+"ptsum_away_100"], 19+ih, 1, 3);
78 }
79 if (_sqs == "" && !merging()) {
80 throw BeamError("Invalid beam energy for " + name() + "\n");
81 }
82
83 // Profiles vs. eta (7 TeV only)
84 book(_p["nch_vs_eta_transverse_100"], 21, 1, 1);
85 book(_p["ptsum_vs_eta_transverse_100"], 22, 1, 1);
86
87 }
88
89
90 // Little helper function to identify Delta(phi) regions
91 inline int region_index(double dphi) const {
92 assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
93 if (dphi < PI/3.0) return 0;
94 if (dphi < 2*PI/3.0) return 1;
95 return 2;
96 }
97
98
99 void analyze(const Event& event) {
100 // Require at least one track in the event with pT >= 1 GeV
101 const ChargedFinalState& cfslead = apply<ChargedFinalState>(event, "CFSlead");
102 if (cfslead.size() < 1) vetoEvent;
103
104 // These are the charged particles (tracks) with pT > 500 MeV
105 const ChargedFinalState& charged500 = apply<ChargedFinalState>(event, "CFS500");
106
107 // Identify leading track and its phi and pT (this is the same for both the 100 MeV and 500 MeV track cuts)
108 Particles particles500 = charged500.particlesByPt();
109 Particle p_lead = particles500[0];
110 const double philead = p_lead.phi();
111 const double etalead = p_lead.eta();
112 const double pTlead = p_lead.perp();
113 MSG_DEBUG("Leading track: pT = " << pTlead << ", eta = " << etalead << ", phi = " << philead);
114
115 // Iterate over all > 500 MeV particles and count particles and scalar pTsum in the three regions
116 vector<double> num500(3, 0), ptSum500(3, 0.0);
117 // Temporary histos that bin Nch and pT in dPhi.
118 // NB. Only one of each needed since binnings are the same for the energies and pT cuts
119 Histo1D htemp_num(_p[_sqs+"N_vs_dPhi_1_500"]->binning());
120 Histo1D htemp_pt(_p[_sqs+"pT_vs_dPhi_1_500"]->binning());
121 for (const Particle& p : particles500) {
122 const double pT = p.pT();
123 const double dPhi = deltaPhi(philead, p.phi());
124 const int ir = region_index(dPhi);
125 num500[ir] += 1;
126 ptSum500[ir] += pT;
127
128 // Fill temp histos to bin Nch and pT in dPhi
129 if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
130 htemp_num.fill(dPhi, 1);
131 htemp_pt.fill(dPhi, pT);
132 }
133 }
134
135
136 // Iterate over charged particles again for profiles against Nch
137 // This is necessary since the Nch are region-specific and so are only known after the first loop
138 for (const Particle& p : particles500) {
139 const double pT = p.pT();
140 const double dPhi = deltaPhi(philead, p.phi());
141 const int ir = region_index(dPhi);
142 switch (ir) {
143 case 0:
144 _p[_sqs+"dn_dpt_toward_500"]->fill(num500[0], pT);
145 break;
146 case 1:
147 _p[_sqs+"dn_dpt_transverse_500"]->fill(num500[1], pT);
148 break;
149 case 2:
150 _p[_sqs+"dn_dpt_away_500"]->fill(num500[2], pT);
151 break;
152 default:
153 assert(false && "How did we get here?");
154 }
155 }
156
157
158 // Now fill underlying event histograms
159 // The densities are calculated by dividing the UE properties by dEta*dPhi
160 // -- each region has a dPhi of 2*PI/3 and dEta is two times 2.5
161 const double dEtadPhi = (2*2.5 * 2*PI/3.0);
162 // Transverse profiles need 4 orders of moments for stddev with errors
163 for (int i = 0; i < 4; ++i) {
164 const string mom = toString(i);
165 _p[_sqs+"nch_transverse_500_"+mom]->fill(pTlead/GeV, intpow(num500[1]/dEtadPhi, i+1));
166 _p[_sqs+"ptsum_transverse_500_"+mom]->fill(pTlead/GeV, intpow(ptSum500[1]/GeV/dEtadPhi, i+1));
167 }
168 // Toward and away profiles only need the first moment (the mean)
169 _p[_sqs+"nch_toward_500"]->fill(pTlead/GeV, num500[0]/dEtadPhi);
170 _p[_sqs+"nch_away_500"]->fill(pTlead/GeV, num500[2]/dEtadPhi);
171 _p[_sqs+"ptsum_toward_500"]->fill(pTlead/GeV, ptSum500[0]/GeV/dEtadPhi);
172 _p[_sqs+"ptsum_away_500"]->fill(pTlead/GeV, ptSum500[2]/GeV/dEtadPhi);
173 // <pT> profiles
174 //MSG_INFO("Trans pT1, pTsum, Nch, <pT>" << pTlead/GeV << ", " << ptSum500[1]/GeV << ", " << num500[1] << ", " << ptSum500[1]/GeV/num500[1]);
175 if (num500[1] > 0) _p[_sqs+"ptavg_transverse_500"]->fill(pTlead/GeV, ptSum500[1]/GeV/num500[1]);
176 if (num500[0] > 0) _p[_sqs+"ptavg_toward_500"]->fill(pTlead/GeV, ptSum500[0]/GeV/num500[0]);
177 if (num500[2] > 0) _p[_sqs+"ptavg_away_500"]->fill(pTlead/GeV, ptSum500[2]/GeV/num500[2]);
178
179
180 // Update the "proper" dphi profile histograms
181 // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
182 // The values tabulated in the note are for an (undefined) signed Delta(phi) rather than
183 // |Delta(phi)| and so differ by a factor of 2: we have to actually norm for angular range = 2pi
184 std::vector<double> ptcut;
185 if (_sqs == "900"s) {
186 ptcut += 1.0; ptcut += 1.5; ptcut += 2.0; ptcut += 2.5;
187 }
188 else if (_sqs == "7000"s) {
189 ptcut += 1.0; ptcut += 2.0; ptcut += 3.0; ptcut += 5.0;
190 }
191 assert(ptcut.size() == 4);
192 for (size_t i = 1; i < htemp_num.numBins()+1; ++i) {
193 // First Nch
194 double mean = htemp_num.bin(i).xMid();
195 double value = 0.;
196 if (htemp_num.bin(i).numEntries() > 0) {
197 mean = htemp_num.bin(i).xMean();
198 value = htemp_num.bin(i).sumW()/htemp_num.bin(i).xWidth()/10.0;
199 }
200 if (pTlead/GeV >= ptcut[0]) _p[_sqs+"N_vs_dPhi_1_500"]->fill(mean, value);
201 if (pTlead/GeV >= ptcut[1]) _p[_sqs+"N_vs_dPhi_2_500"]->fill(mean, value);
202 if (pTlead/GeV >= ptcut[2]) _p[_sqs+"N_vs_dPhi_3_500"]->fill(mean, value);
203 if (pTlead/GeV >= ptcut[3]) _p[_sqs+"N_vs_dPhi_5_500"]->fill(mean, value);
204
205 // Then pT
206 mean = htemp_pt.bin(i).xMid();
207 value = 0.;
208 if (htemp_pt.bin(i).numEntries() > 0) {
209 mean = htemp_pt.bin(i).xMean();
210 value = htemp_pt.bin(i).sumW()/htemp_pt.bin(i).xWidth()/10.0;
211 }
212 if (pTlead/GeV >= ptcut[0]) _p[_sqs+"pT_vs_dPhi_1_500"]->fill(mean, value);
213 if (pTlead/GeV >= ptcut[1]) _p[_sqs+"pT_vs_dPhi_2_500"]->fill(mean, value);
214 if (pTlead/GeV >= ptcut[2]) _p[_sqs+"pT_vs_dPhi_3_500"]->fill(mean, value);
215 if (pTlead/GeV >= ptcut[3]) _p[_sqs+"pT_vs_dPhi_5_500"]->fill(mean, value);
216 }
217
218
219 //////////////////////
220
221
222 // These are the charged particles (tracks) with pT > 100 MeV
223 const ChargedFinalState& charged100 = apply<ChargedFinalState>(event, "CFS100");
224
225 // Iterate over all > 100 MeV particles and count particles and scalar pTsum in the three regions
226 vector<double> num100(3, 0), ptSum100(3, 0.0);
227 for (const Particle& p : charged100.particles()) {
228 const double pT = p.pT();
229 const double dPhi = deltaPhi(philead, p.phi());
230 const int ir = region_index(dPhi);
231 num100[ir] += 1;
232 ptSum100[ir] += pT;
233 }
234
235 // Now fill the two sets of 100 MeV underlying event histograms
236 _p[_sqs+"nch_transverse_100"]->fill(pTlead/GeV, num100[1]/dEtadPhi);
237 _p[_sqs+"nch_toward_100"]->fill(pTlead/GeV, num100[0]/dEtadPhi);
238 _p[_sqs+"nch_away_100"]->fill(pTlead/GeV, num100[2]/dEtadPhi);
239 _p[_sqs+"ptsum_transverse_100"]->fill(pTlead/GeV, ptSum100[1]/GeV/dEtadPhi);
240 _p[_sqs+"ptsum_toward_100"]->fill(pTlead/GeV, ptSum100[0]/GeV/dEtadPhi);
241 _p[_sqs+"ptsum_away_100"]->fill(pTlead/GeV, ptSum100[2]/GeV/dEtadPhi);
242
243 // And finally the Nch and pT vs eta_lead profiles (again from > 100 MeV tracks, and only at 7 TeV)
244 if (_sqs == "7000"s && pTlead > 5*GeV) {
245 _p["nch_vs_eta_transverse_100"]->fill(etalead, num100[1]/dEtadPhi);
246 _p["ptsum_vs_eta_transverse_100"]->fill(etalead, ptSum100[1]/GeV/dEtadPhi);
247 }
248
249 }
250
251
252 void finalize() {
253 // Convert the various moments of the 500 MeV trans pT and Nch distributions to std devs with correct error
254 for (double eVal : allowedEnergies()) {
255 const string en = toString(int(eVal));
256 moments_to_stddev(en+"nch_transverse_500", _e["nch_transverse_500"]);
257 moments_to_stddev(en+"ptsum_transverse_500", _e["ptsum_transverse_500"]);
258 }
259 }
260
261 void moments_to_stddev(const string& label, Estimate1DPtr target_dps) {
262 for (size_t b = 1; b < _p[label+"_0"]->numBins()+1; ++b) { // loop over points
263 /// @todo Assuming unit weights here! Should use N_effective = sumW**2/sumW2?
264 const double numentries = _p[label+"_0"]->bin(b).numEntries();
265 double var = 0., sd = 0.;
266 if (numentries > 0) {
267 var = _p[label+"_1"]->bin(b).yMean() - intpow(_p[label+"_0"]->bin(b).yMean(), 2);
268 sd = fuzzyLessEquals(var,0.) ? 0 : sqrt(var); ///< Numerical safety check
269 }
270 if (sd == 0 || numentries < 3) {
271 MSG_WARNING("Need at least 3 bin entries and a non-zero central value to calculate "
272 << "an error on standard deviation profiles (bin " << b << ")");
273 target_dps->bin(b).set(sd, 0);
274 continue;
275 }
276 // c2(y) = m4(x) - 4 m3(x) m1(x) - m2(x)^2 + 8 m2(x) m1(x)^2 - 4 m1(x)^4
277 const double var_on_var = _p[label+"_3"]->bin(b).yMean()
278 - 4 * _p[label+"_2"]->bin(b).yMean() * _p[label+"_0"]->bin(b).yMean()
279 - intpow(_p[label+"_1"]->bin(b).yMean(), 2)
280 + 8 * _p[label+"_1"]->bin(b).yMean() * intpow(_p[label+"_0"]->bin(b).yMean(), 2)
281 - 4 * intpow(_p[label+"_0"]->bin(b).yMean(), 4);
282 const double stderr_on_var = sqrt(var_on_var/(numentries-2.0));
283 const double stderr_on_sd = stderr_on_var / (2.0*sd);
284 target_dps->bin(b).set(sd, stderr_on_sd);
285 }
286 }
287
288
289 private:
290
291 map<string,Profile1DPtr> _p;
292 map<string,Estimate1DPtr> _e;
293
294 string _sqs = "";
295
296 };
297
298
299
300 RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2010_I879407, ATLAS_2010_S8894728);
301
302}
|