rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2012_I1119557

Jet shapes and jet masses
Experiment: ATLAS (LHC 7TeV)
Inspire ID: 1119557
Status: VALIDATED
Authors:
  • Lily Asquith
  • Roman Lysak
References: Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • QCD events at 7 TeV, leading-pT jets with $p_\perp > 300 GeV$.

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, FastJets::ANTIKT, 0.6);
 32      declare(fj06, "AntiKT06");
 33      FastJets fj10(fs, FastJets::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}