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