Rivet analyses referenceATLAS_2012_I1119557Jet shapes and jet massesExperiment: ATLAS (LHC 7TeV) Inspire ID: 1119557 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
Measurements are presented of the properties of high transverse momentum jets, produced in proton-proton collisions at a center-of-mass energy of $\sqrt{s} = 7$ TeV. Jet mass, width, eccentricity, planar flow and angularity are measured for jets reconstructed using the anti-$k_t$ algorithm with distance parameters $R = 0.6$ and 1.0, with transverse momentum $pT > 300$ GeV and pseudorapidity $|\eta| < 2$. Source code: ATLAS_2012_I1119557.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FastJets.hh"
4#include "Rivet/Projections/FinalState.hh"
5
6namespace Rivet {
7
8
9 namespace { // unnamed namespace
10 // Forward declarations of calculator functions: implementations at bottom of file
11 double getWidth(const Jet& jet);
12 /// @todo Re-enable eccentricity calculation
13 // double getEcc(const Jet& jet);
14 double getPFlow(const Jet& jet);
15 double getAngularity(const Jet& jet);
16 }
17
18
19 class ATLAS_2012_I1119557 : public Analysis {
20 public:
21
22 ATLAS_2012_I1119557()
23 : Analysis("ATLAS_2012_I1119557")
24 { }
25
26
27 void init() {
28 const FinalState fs;
29 declare(fs, "FinalState");
30
31 FastJets fj06(fs, JetAlg::ANTIKT, 0.6);
32 declare(fj06, "AntiKT06");
33 FastJets fj10(fs, JetAlg::ANTIKT, 1.0);
34 declare(fj10, "AntiKT10");
35
36 for (size_t alg = 0; alg < 2; ++alg) {
37 book(_hs_mass[alg] ,1, alg+1, 1);
38 book(_hs_width[alg] ,2, alg+1, 1);
39 /// @todo Commented eccentricity out for now: reinstate
40 //book( _hs_eccentricity[alg] ,3, alg+1, 1);
41 }
42 book(_h_planarFlow ,4, 2, 1);
43 book(_h_angularity ,5, 1, 1);
44
45 }
46
47
48 /// Perform the per-event analysis
49 void analyze(const Event& event) {
50 Jets jetAr[2];
51 jetAr[0] = apply<FastJets>(event, "AntiKT06").jetsByPt(Cuts::pT > 300*GeV && Cuts::abseta < 2.0);
52 jetAr[1] = apply<FastJets>(event, "AntiKT10").jetsByPt(Cuts::pT > 300*GeV && Cuts::abseta < 2.0);
53
54 for (size_t alg = 0; alg < 2; ++alg) {
55 // Require at least one jet
56 if (jetAr[alg].size() < 1) continue;
57
58 // The leading jet
59 const Jet& jet = jetAr[alg][0];
60 const double m = jet.mass();
61 const double eta = jet.eta();
62
63 _hs_mass[alg]->fill(m/GeV);
64 _hs_width[alg]->fill(getWidth(jet));
65 /// @todo Commented eccentricity out for now: reinstate
66 // if (fabs(eta) < 0.7 && m > 100*GeV) _hs_eccentricity[alg]->fill(getEcc(jet));
67
68 if (fabs(eta) < 0.7) {
69 if (alg == 0 && inRange(m/GeV, 100., 130.)) _h_angularity->fill(getAngularity(jet));
70 if (alg == 1 && inRange(m/GeV, 130., 210.)) _h_planarFlow->fill(getPFlow(jet));
71 }
72 }
73 }
74
75
76 /// Normalise histograms etc., after the run
77 void finalize() {
78 for (size_t alg = 0; alg < 2; ++alg) {
79 normalize(_hs_mass[alg]);
80 normalize(_hs_width[alg]);
81 /// @todo Commented eccentricity out for now: reinstate
82 // normalize(_hs_eccentricity[alg]);
83 }
84 normalize(_h_planarFlow);
85 normalize(_h_angularity);
86 }
87
88
89 private:
90
91 Histo1DPtr _hs_mass[2];
92 Histo1DPtr _hs_width[2];
93 /// @todo Commented eccentricity out for now: reinstate
94 // Histo1DPtr _hs_eccentricity[2];
95 Histo1DPtr _h_planarFlow;
96 Histo1DPtr _h_angularity;
97 };
98
99
100 RIVET_DECLARE_PLUGIN(ATLAS_2012_I1119557);
101
102
103
104 namespace { // unnamed namespace
105
106 // Adapted code from Lily
107 /// @todo Convert to use the Rivet rotation matrix code (should be simpler)
108 FourMomentum RotateAxes(const Rivet::FourMomentum& p, double M[3][3]){
109 double px_rot = M[0][0]*(p.px()) + M[0][1]*(p.py()) + M[0][2]*(p.pz());
110 double py_rot = M[1][0]*(p.px()) + M[1][1]*(p.py()) + M[1][2]*(p.pz());
111 double pz_rot = M[2][0]*(p.px()) + M[2][1]*(p.py()) + M[2][2]*(p.pz());
112 return FourMomentum(p.E(), px_rot, py_rot, pz_rot);
113 }
114
115 // Untouched code from Lily
116 /// @todo Convert to use the Rivet rotation matrix code (should be simpler)
117 void CalcRotationMatrix(double nvec[3],double rot_mat[3][3]){
118 // clear momentum tensor
119 for (size_t i = 0; i < 3; i++) {
120 for (size_t j = 0; j < 3; j++) {
121 rot_mat[i][j] = 0.;
122 }
123 }
124 double mag3 = sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1]+ nvec[2]*nvec[2]);
125 double mag2 = sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1]);
126 /// @todo cout is not a valid response to a numerical error! Is the error condition reported?!? Assert added by AB for Rivet 1.8.2
127 assert(mag3 > 0);
128 if (mag3 <= 0) {
129 cout << "rotation axis is null" << '\n';
130 return;
131 }
132
133 double ctheta0 = nvec[2]/mag3;
134 double stheta0 = mag2/mag3;
135 double cphi0 = (mag2 > 0.) ? nvec[0]/mag2 : 0;
136 double sphi0 = (mag2 > 0.) ? nvec[1]/mag2 : 0;
137
138 rot_mat[0][0] = -ctheta0*cphi0;
139 rot_mat[0][1] = -ctheta0*sphi0;
140 rot_mat[0][2] = stheta0;
141 rot_mat[1][0] = sphi0;
142 rot_mat[1][1] = -1.*cphi0;
143 rot_mat[1][2] = 0.;
144 rot_mat[2][0] = stheta0*cphi0;
145 rot_mat[2][1] = stheta0*sphi0;
146 rot_mat[2][2] = ctheta0;
147 }
148
149
150 /// Jet width calculation
151 double getWidth(const Jet& jet) {
152 const double phi_jet = jet.phi();
153 const double eta_jet = jet.eta();
154 double width(0), pTsum(0);
155 for (const Particle& p : jet.particles()) {
156 double pT = p.pT();
157 double eta = p.eta();
158 double phi = p.phi();
159 width += sqrt(pow(phi_jet - phi,2) + pow(eta_jet - eta, 2)) * pT;
160 pTsum += pT;
161 }
162 const double rtn = (pTsum != 0.0) ? width/pTsum : -1;
163 return rtn;
164 }
165
166
167 /// Eccentricity calculation, copied and adapted from Lily's code
168 /// @todo Re-enable eccentricity calculation
169 // double getEcc(const Jet& jet) {
170 // vector<double> phis;
171 // vector<double> etas;
172 // vector<double> energies;
173
174 // double etaSum(0), phiSum(0), eTot(0);
175 // for (const Particle& p : jet.particles()) {
176 // const double E = p.E();
177 // const double eta = p.eta();
178
179 // energies.push_back(E);
180 // etas.push_back(jet.eta() - eta);
181
182 // eTot += E;
183 // etaSum += eta * E;
184
185 // /// @todo Replace with the Rivet deltaPhi function (or the mapAngleTo* functions)
186 // double dPhi = jet.phi() - p.phi();
187 // //if DPhi does not lie within 0 < DPhi < PI take 2*PI off DPhi
188 // //this covers cases where DPhi is greater than PI
189 // if( fabs( dPhi - TWOPI ) < fabs(dPhi) ) dPhi -= TWOPI;
190 // //if DPhi does not lie within -PI < DPhi < 0 add 2*PI to DPhi
191 // //this covers cases where DPhi is less than -PI
192 // else if( fabs(dPhi + TWOPI) < fabs(dPhi) ) dPhi += TWOPI;
193 // phis.push_back(dPhi);
194
195 // phiSum += dPhi * E;
196 // }
197
198 // //these are the "pull" away from the jet axis
199 // etaSum = etaSum/eTot;
200 // phiSum = phiSum/eTot;
201
202 // // now for every cluster we alter its position by moving it:
203 // // away from the new axis if it is in the direction of -ve pull
204 // // closer to the new axis if it is in the direction of +ve pull
205 // // the effect of this will be that the new energy weighted center will be on the old jet axis.
206 // double little_x(0), little_y(0);
207 // for (size_t k = 0; k < jet.particles().size(); ++k) {
208 // little_x+= etas[k]-etaSum;
209 // little_y+= phis[k]-phiSum;
210 // etas[k] = etas[k]-etaSum;
211 // phis[k] = phis[k]-phiSum;
212 // }
213
214 // double x1(0), x2(0);
215 // for (size_t i = 0; i < jet.particles().size(); ++i) {
216 // x1 += 2. * energies[i]* etas[i] * phis[i]; // this is 2*X*Y
217 // x2 += energies[i]*(phis[i] * phis[i] - etas[i] * etas[i] ); // this is X^2 - Y^2
218 // }
219
220 // // Variance calculations
221 // double theta = .5*atan2(x1, x2);
222 // double sinTheta =sin(theta);
223 // double cosTheta = cos(theta);
224 // double theta2 = theta + 0.5*PI;
225 // double sinThetaPrime = sin(theta2);
226 // double cosThetaPrime = cos(theta2);
227
228 // double varX(0), varY(0);
229 // for (size_t i = 0; i < jet.particles().size(); i++) {
230 // const double x = sinTheta*etas[i] + cosTheta*phis[i];
231 // const double y = sinThetaPrime*etas[i] + cosThetaPrime*phis[i];
232 // varX += energies[i]* sqr(x);
233 // varY += energies[i]* sqr(y);
234 // }
235 // const double varianceMax = max(varX, varY);
236 // const double varianceMin = min(varX, varY);
237 // const double ecc = (varianceMax != 0.0) ? 1 - varianceMin/varianceMax : -1;
238 // return ecc;
239 // }
240
241
242 /// Planar flow calculation, copied and adapted from Lily's code
243 double getPFlow(const Jet& jet) {
244 const double phi0 = jet.phi();
245 const double eta0 = jet.eta();
246
247 double nref[3]; ///< @todo 3-vector to rotate x to? Use Rivet vector classes
248 nref[0] = cos(phi0)/cosh(eta0);
249 nref[1] = sin(phi0)/cosh(eta0);
250 nref[2] = tanh(eta0);
251
252 // Rotation matrix to align with nref
253 double rotationMatrix[3][3];
254 CalcRotationMatrix(nref, rotationMatrix);
255
256 double iw00(0.), iw01(0.), iw11(0.), iw10(0.);
257 for (const Particle& p : jet.particles()) {
258 double a = 1./(p.E()*jet.mass());
259 FourMomentum rotclus = RotateAxes(p.momentum(), rotationMatrix);
260 iw00 += a*pow(rotclus.px(), 2);
261 iw01 += a*rotclus.px()*rotclus.py();
262 iw10 += a*rotclus.py()*rotclus.px();
263 iw11 += a*pow(rotclus.py(), 2);
264 }
265
266 const double det = iw00*iw11 - iw01*iw10;
267 const double trace = iw00 + iw11;
268
269 const double pf = (trace != 0.0) ? (4.0*det)/sqr(trace) : -1;
270 return pf;
271 }
272
273
274 /// Angularity calculation, copied and adapted from Lily's code
275 double getAngularity(const Jet& jet) {
276 double sum_a = 0.;
277 // a can take any value < 2 (e.g. 1,0,-0.5 etc) for infrared safety
278 const double a = -2.;
279 for (const Particle& p : jet.particles()) {
280 double e_i = p.E();
281 double theta_i = jet.momentum().angle(p.momentum());
282 double e_theta_i = e_i * pow(sin(theta_i), a) * pow(1-cos(theta_i), 1-a);
283 sum_a += e_theta_i;
284 }
285 const double rtn = (jet.mass() != 0.0 && !std::isnan(sum_a)) ? sum_a/jet.mass() : -1;
286 return rtn;
287 }
288
289 }
290
291
292}
|