Rivet analyses

title: ATLAS_2012_I1119557

Jet shapes and jet masses

Experiment: ATLAS (LHC 7TeV)

Inspire ID: 1119557

Status: VALIDATED

Authors: - Lily Asquith - Roman Lysak

References: - Expt page: ATLAS-STDM-2011-38 - Phys.Rev. D86 (2012) 072006 - DOI: 10.1103/PhysRevD.86.072006 - arXiv: 1206.5369

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - QCD events at 7 TeV, leading-pT jets with $\pT > 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

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

namespace Rivet {

namespace { // unnamed namespace // Forward declarations of calculator functions: implementations at bottom of file double getWidth(const Jet& jet); /// @todo Re-enable eccentricity calculation // double getEcc(const Jet& jet); double getPFlow(const Jet& jet); double getAngularity(const Jet& jet); }

class ATLAS_2012_I1119557 : public Analysis { public:

ATLAS_2012_I1119557()
    : Analysis("ATLAS_2012_I1119557") { }


void init() {
  const FinalState fs;
  declare(fs, "FinalState");

  FastJets fj06(fs, JetAlg::ANTIKT, 0.6);
  declare(fj06, "AntiKT06");
  FastJets fj10(fs, JetAlg::ANTIKT, 1.0);
  declare(fj10, "AntiKT10");

  for (size_t alg = 0; alg < 2; ++alg) {
    book(_hs_mass[alg], 1, alg + 1, 1);
    book(_hs_width[alg], 2, alg + 1, 1);
    /// @todo Commented eccentricity out for now: reinstate
    //book( _hs_eccentricity[alg] ,3, alg+1, 1);
  }
  book(_h_planarFlow, 4, 2, 1);
  book(_h_angularity, 5, 1, 1);
}


/// Perform the per-event analysis
void analyze(const Event& event) {
  Jets jetAr[2];
  jetAr[0] = apply<FastJets>(event, "AntiKT06").jetsByPt(Cuts::pT > 300 * GeV && Cuts::abseta < 2.0);
  jetAr[1] = apply<FastJets>(event, "AntiKT10").jetsByPt(Cuts::pT > 300 * GeV && Cuts::abseta < 2.0);

  for (size_t alg = 0; alg < 2; ++alg) {
    // Require at least one jet
    if (jetAr[alg].size() < 1) continue;

    // The leading jet
    const Jet& jet = jetAr[alg][0];
    const double m = jet.mass();
    const double eta = jet.eta();

    _hs_mass[alg]->fill(m / GeV);
    _hs_width[alg]->fill(getWidth(jet));
    /// @todo Commented eccentricity out for now: reinstate
    // if (fabs(eta) < 0.7 && m > 100*GeV) _hs_eccentricity[alg]->fill(getEcc(jet));

    if (fabs(eta) < 0.7) {
      if (alg == 0 && inRange(m / GeV, 100., 130.)) _h_angularity->fill(getAngularity(jet));
      if (alg == 1 && inRange(m / GeV, 130., 210.)) _h_planarFlow->fill(getPFlow(jet));
    }
  }
}


/// Normalise histograms etc., after the run
void finalize() {
  for (size_t alg = 0; alg < 2; ++alg) {
    normalize(_hs_mass[alg]);
    normalize(_hs_width[alg]);
    /// @todo Commented eccentricity out for now: reinstate
    // normalize(_hs_eccentricity[alg]);
  }
  normalize(_h_planarFlow);
  normalize(_h_angularity);
}

private:

Histo1DPtr _hs_mass[2];
Histo1DPtr _hs_width[2];
/// @todo Commented eccentricity out for now: reinstate
// Histo1DPtr _hs_eccentricity[2];
Histo1DPtr _h_planarFlow;
Histo1DPtr _h_angularity;

};

RIVET_DECLARE_PLUGIN(ATLAS_2012_I1119557);

namespace { // unnamed namespace

// Adapted code from Lily
/// @todo Convert to use the Rivet rotation matrix code (should be simpler)
FourMomentum RotateAxes(const Rivet::FourMomentum& p, double M[3][3]) {
  double px_rot = M[0][0] * (p.px()) + M[0][1] * (p.py()) + M[0][2] * (p.pz());
  double py_rot = M[1][0] * (p.px()) + M[1][1] * (p.py()) + M[1][2] * (p.pz());
  double pz_rot = M[2][0] * (p.px()) + M[2][1] * (p.py()) + M[2][2] * (p.pz());
  return FourMomentum(p.E(), px_rot, py_rot, pz_rot);
}

// Untouched code from Lily
/// @todo Convert to use the Rivet rotation matrix code (should be simpler)
void CalcRotationMatrix(double nvec[3], double rot_mat[3][3]) {
  // clear momentum tensor
  for (size_t i = 0; i < 3; i++) {
    for (size_t j = 0; j < 3; j++) {
      rot_mat[i][j] = 0.;
    }
  }
  double mag3 = sqrt(nvec[0] * nvec[0] + nvec[1] * nvec[1] + nvec[2] * nvec[2]);
  double mag2 = sqrt(nvec[0] * nvec[0] + nvec[1] * nvec[1]);
  /// @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
  assert(mag3 > 0);
  if (mag3 <= 0) {
    cout << "rotation axis is null" << '\n';
    return;
  }

  double ctheta0 = nvec[2] / mag3;
  double stheta0 = mag2 / mag3;
  double cphi0 = (mag2 > 0.) ? nvec[0] / mag2 : 0;
  double sphi0 = (mag2 > 0.) ? nvec[1] / mag2 : 0;

  rot_mat[0][0] = -ctheta0 * cphi0;
  rot_mat[0][1] = -ctheta0 * sphi0;
  rot_mat[0][2] = stheta0;
  rot_mat[1][0] = sphi0;
  rot_mat[1][1] = -1. * cphi0;
  rot_mat[1][2] = 0.;
  rot_mat[2][0] = stheta0 * cphi0;
  rot_mat[2][1] = stheta0 * sphi0;
  rot_mat[2][2] = ctheta0;
}


/// Jet width calculation
double getWidth(const Jet& jet) {
  const double phi_jet = jet.phi();
  const double eta_jet = jet.eta();
  double width(0), pTsum(0);
  for (const Particle& p : jet.particles()) {
    double pT = p.pT();
    double eta = p.eta();
    double phi = p.phi();
    width += sqrt(pow(phi_jet - phi, 2) + pow(eta_jet - eta, 2)) * pT;
    pTsum += pT;
  }
  const double rtn = (pTsum != 0.0) ? width / pTsum : -1;
  return rtn;
}


/// Eccentricity calculation, copied and adapted from Lily's code
/// @todo Re-enable eccentricity calculation
// double getEcc(const Jet& jet) {
//   vector<double> phis;
//   vector<double> etas;
//   vector<double> energies;

//   double etaSum(0), phiSum(0), eTot(0);
//   for (const Particle& p : jet.particles()) {
//     const double E = p.E();
//     const double eta = p.eta();

//     energies.push_back(E);
//     etas.push_back(jet.eta() - eta);

//     eTot   += E;
//     etaSum += eta * E;

//     /// @todo Replace with the Rivet deltaPhi function (or the mapAngleTo* functions)
//     double dPhi = jet.phi() - p.phi();
//     //if DPhi does not lie within 0 < DPhi < PI take 2*PI off DPhi
//     //this covers cases where DPhi is greater than PI
//     if( fabs( dPhi - TWOPI ) < fabs(dPhi) ) dPhi -= TWOPI;
//     //if DPhi does not lie within -PI < DPhi < 0 add 2*PI to DPhi
//     //this covers cases where DPhi is less than -PI
//     else if( fabs(dPhi + TWOPI) < fabs(dPhi) ) dPhi += TWOPI;
//     phis.push_back(dPhi);

//     phiSum += dPhi * E;
//   }

//   //these are the "pull" away from the jet axis
//   etaSum = etaSum/eTot;
//   phiSum = phiSum/eTot;

//   // now for every cluster we alter its position by moving it:
//   // away from the new axis if it is in the direction of -ve pull
//   // closer to the new axis if it is in the direction of +ve pull
//   // the effect of this will be that the new energy weighted center will be on the old jet axis.
//   double little_x(0), little_y(0);
//   for (size_t k = 0; k < jet.particles().size(); ++k) {
//     little_x+= etas[k]-etaSum;
//     little_y+= phis[k]-phiSum;
//     etas[k] = etas[k]-etaSum;
//     phis[k] = phis[k]-phiSum;
//   }

//   double x1(0), x2(0);
//   for (size_t i = 0; i < jet.particles().size(); ++i) {
//     x1 += 2. * energies[i]* etas[i] * phis[i]; // this is 2*X*Y
//     x2 += energies[i]*(phis[i] * phis[i] - etas[i] * etas[i] ); // this is X^2 - Y^2
//   }

//   // Variance calculations
//   double theta = .5*atan2(x1, x2);
//   double sinTheta =sin(theta);
//   double cosTheta = cos(theta);
//   double theta2 = theta + 0.5*PI;
//   double sinThetaPrime = sin(theta2);
//   double cosThetaPrime = cos(theta2);

//   double varX(0), varY(0);
//   for (size_t i = 0; i < jet.particles().size(); i++) {
//     const double x = sinTheta*etas[i] + cosTheta*phis[i];
//     const double y = sinThetaPrime*etas[i] + cosThetaPrime*phis[i];
//     varX += energies[i]* sqr(x);
//     varY += energies[i]* sqr(y);
//   }
//   const double varianceMax = max(varX, varY);
//   const double varianceMin = min(varX, varY);
//   const double ecc = (varianceMax != 0.0) ? 1 - varianceMin/varianceMax : -1;
//   return ecc;
// }


/// Planar flow calculation, copied and adapted from Lily's code
double getPFlow(const Jet& jet) {
  const double phi0 = jet.phi();
  const double eta0 = jet.eta();

  double nref[3]; ///< @todo 3-vector to rotate x to? Use Rivet vector classes
  nref[0] = cos(phi0) / cosh(eta0);
  nref[1] = sin(phi0) / cosh(eta0);
  nref[2] = tanh(eta0);

  // Rotation matrix to align with nref
  double rotationMatrix[3][3];
  CalcRotationMatrix(nref, rotationMatrix);

  double iw00(0.), iw01(0.), iw11(0.), iw10(0.);
  for (const Particle& p : jet.particles()) {
    double a = 1. / (p.E() * jet.mass());
    FourMomentum rotclus = RotateAxes(p.momentum(), rotationMatrix);
    iw00 += a * pow(rotclus.px(), 2);
    iw01 += a * rotclus.px() * rotclus.py();
    iw10 += a * rotclus.py() * rotclus.px();
    iw11 += a * pow(rotclus.py(), 2);
  }

  const double det = iw00 * iw11 - iw01 * iw10;
  const double trace = iw00 + iw11;

  const double pf = (trace != 0.0) ? (4.0 * det) / sqr(trace) : -1;
  return pf;
}


/// Angularity calculation, copied and adapted from Lily's code
double getAngularity(const Jet& jet) {
  double sum_a = 0.;
  // a can take any value < 2 (e.g. 1,0,-0.5 etc) for infrared safety
  const double a = -2.;
  for (const Particle& p : jet.particles()) {
    double e_i = p.E();
    double theta_i = jet.momentum().angle(p.momentum());
    double e_theta_i = e_i * pow(sin(theta_i), a) * pow(1 - cos(theta_i), 1 - a);
    sum_a += e_theta_i;
  }
  const double rtn = (jet.mass() != 0.0 && !std::isnan(sum_a)) ? sum_a / jet.mass() : -1;
  return rtn;
}

}

} ```