Rivet analyses

title: ATLAS_2020_I1808726

Hadronic event shapes in multijet final states

Experiment: ATLAS (LHC)

Inspire ID: 1808726

Status: VALIDATED

Authors: - Javier Llorente - Deepak Kar

References: - Expt page: ATLAS-STDM-2019-02 - arXiv: 2007.12600

Beams: p+ p+

Beam energies: (6500.0, 6500.0)GeV

Run details: - Jet production at 13000 GeV. Requires $H_{T2} \equiv \sum p_T$ of two lead jets $> 1$ TeV. Observables include up to 6th jet bin with $H_{T2} > 2$ TeV.

A measurement of event-shape variables in proton-proton collisions at large momentum transfer is presented using data collected at $\sqrt{s}$=13 TeV with the ATLAS detector at the Large Hadron Collider. Six event-shape variables calculated using hadronic jets are studied in inclusive multijet events using data corresponding to an integrated luminosity of 139 fb$^{-1}$. Measurements are performed in bins of jet multiplicity and in different ranges of the scalar sum of the transverse momenta of the two leading jets, reaching scales beyond 2 TeV. These measurements are compared with predictions from Monte Carlo event generators containing leading-order or next-to-leading order matrix elements matched to parton showers simulated to leading-logarithm accuracy. At low jet multiplicities, shape discrepancies between the measurements and the Monte Carlo predictions are observed. At high jet multiplicities, the shapes are better described but discrepancies in the normalisation are observed.

Source code:ATLAS_2020_I1808726.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/Sphericity.hh"

include "Rivet/Projections/Thrust.hh"

namespace Rivet {

/// @brief Multijet event shapes at 13 TeV class ATLAS_2020_I1808726 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2020_I1808726);

double xs1 = 0.0;
double xs2 = 0.0;
double xs3 = 0.0;

/// Initialization, called once before running
void init() {

  //Jet collection (excluding muons and neutrinos)
  const FinalState fs(Cuts::abseta < 4.5);

  FastJets jets(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
  declare(jets, "Jets");

  // Book histograms
  //Jet multiplicity
  book(_h["njet_h1"], 73, 1, 1);
  book(_h["njet_h2"], 74, 1, 1);
  book(_h["njet_h3"], 75, 1, 1);

  //Transverse thrust
  book(_h["transThrust_j3_h1"], 1, 1, 1);
  book(_h["transThrust_j3_h2"], 5, 1, 1);
  book(_h["transThrust_j3_h3"], 9, 1, 1);
  book(_h["transThrust_j4_h1"], 2, 1, 1);
  book(_h["transThrust_j4_h2"], 6, 1, 1);
  book(_h["transThrust_j4_h3"], 10, 1, 1);
  book(_h["transThrust_j5_h1"], 3, 1, 1);
  book(_h["transThrust_j5_h2"], 7, 1, 1);
  book(_h["transThrust_j5_h3"], 11, 1, 1);
  book(_h["transThrust_j6_h1"], 4, 1, 1);
  book(_h["transThrust_j6_h2"], 8, 1, 1);
  book(_h["transThrust_j6_h3"], 12, 1, 1);

  //Thrust minor
  book(_h["transMinor_j3_h1"], 13, 1, 1);
  book(_h["transMinor_j3_h2"], 17, 1, 1);
  book(_h["transMinor_j3_h3"], 21, 1, 1);
  book(_h["transMinor_j4_h1"], 14, 1, 1);
  book(_h["transMinor_j4_h2"], 18, 1, 1);
  book(_h["transMinor_j4_h3"], 22, 1, 1);
  book(_h["transMinor_j5_h1"], 15, 1, 1);
  book(_h["transMinor_j5_h2"], 19, 1, 1);
  book(_h["transMinor_j5_h3"], 23, 1, 1);
  book(_h["transMinor_j6_h1"], 16, 1, 1);
  book(_h["transMinor_j6_h2"], 20, 1, 1);
  book(_h["transMinor_j6_h3"], 24, 1, 1);

  //Transverse sphericity
  book(_h["transSphericity_j3_h1"], 25, 1, 1);
  book(_h["transSphericity_j3_h2"], 29, 1, 1);
  book(_h["transSphericity_j3_h3"], 33, 1, 1);
  book(_h["transSphericity_j4_h1"], 26, 1, 1);
  book(_h["transSphericity_j4_h2"], 30, 1, 1);
  book(_h["transSphericity_j4_h3"], 34, 1, 1);
  book(_h["transSphericity_j5_h1"], 27, 1, 1);
  book(_h["transSphericity_j5_h2"], 31, 1, 1);
  book(_h["transSphericity_j5_h3"], 35, 1, 1);
  book(_h["transSphericity_j6_h1"], 28, 1, 1);
  book(_h["transSphericity_j6_h2"], 32, 1, 1);
  book(_h["transSphericity_j6_h3"], 36, 1, 1);

  //Aplanarity
  book(_h["aplanarity_j3_h1"], 37, 1, 1);
  book(_h["aplanarity_j3_h2"], 41, 1, 1);
  book(_h["aplanarity_j3_h3"], 45, 1, 1);
  book(_h["aplanarity_j4_h1"], 38, 1, 1);
  book(_h["aplanarity_j4_h2"], 42, 1, 1);
  book(_h["aplanarity_j4_h3"], 46, 1, 1);
  book(_h["aplanarity_j5_h1"], 39, 1, 1);
  book(_h["aplanarity_j5_h2"], 43, 1, 1);
  book(_h["aplanarity_j5_h3"], 47, 1, 1);
  book(_h["aplanarity_j6_h1"], 40, 1, 1);
  book(_h["aplanarity_j6_h2"], 44, 1, 1);
  book(_h["aplanarity_j6_h3"], 48, 1, 1);

  //C
  book(_h["C_j3_h1"], 49, 1, 1);
  book(_h["C_j3_h2"], 53, 1, 1);
  book(_h["C_j3_h3"], 57, 1, 1);
  book(_h["C_j4_h1"], 50, 1, 1);
  book(_h["C_j4_h2"], 54, 1, 1);
  book(_h["C_j4_h3"], 58, 1, 1);
  book(_h["C_j5_h1"], 51, 1, 1);
  book(_h["C_j5_h2"], 55, 1, 1);
  book(_h["C_j5_h3"], 59, 1, 1);
  book(_h["C_j6_h1"], 52, 1, 1);
  book(_h["C_j6_h2"], 56, 1, 1);
  book(_h["C_j6_h3"], 60, 1, 1);

  //D
  book(_h["D_j3_h1"], 61, 1, 1);
  book(_h["D_j3_h2"], 65, 1, 1);
  book(_h["D_j3_h3"], 69, 1, 1);
  book(_h["D_j4_h1"], 62, 1, 1);
  book(_h["D_j4_h2"], 66, 1, 1);
  book(_h["D_j4_h3"], 70, 1, 1);
  book(_h["D_j5_h1"], 63, 1, 1);
  book(_h["D_j5_h2"], 67, 1, 1);
  book(_h["D_j5_h3"], 71, 1, 1);
  book(_h["D_j6_h1"], 64, 1, 1);
  book(_h["D_j6_h2"], 68, 1, 1);
  book(_h["D_j6_h3"], 72, 1, 1);
}

void analyze(const Event& event) {

  const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 7.0 * GeV);

  //Select jets passing kinematic cuts
  std::vector<const Jet*> goodJets;
  goodJets.clear();
  std::vector<Vector3> momenta2;
  momenta2.clear();
  std::vector<Vector3> momenta3;
  momenta3.clear();

  //foreach (const Jet& j, jets) {
  for (const Jet& j : jets) {
    if (j.abseta() < 2.4 && j.pt() > 100.0 * GeV) {
      goodJets.push_back(&j);

      Vector3 jet2 = j.p3();
      jet2.setZ(0.0);
      momenta2.push_back(jet2);

      Vector3 jet3 = j.p3();
      momenta3.push_back(jet3);
    }
  }

  //Dijet event selection
  if (goodJets.size() < 2) vetoEvent;
  double ht2 = goodJets[0]->pt() + goodJets[1]->pt();
  if (ht2 <= 1000.0 * GeV) vetoEvent;

  //Jet multiplicity
  if (ht2 > 1000.0 * GeV && ht2 < 1500.0 * GeV) _h["njet_h1"]->fill(goodJets.size());
  if (ht2 > 1500.0 * GeV && ht2 < 2000.0 * GeV) _h["njet_h2"]->fill(goodJets.size());
  if (ht2 > 2000.0 * GeV) _h["njet_h3"]->fill(goodJets.size());

  //Thrust calculation
  Thrust thrust;
  thrust.calc(momenta2);
  const double transThrust = 1.0 - thrust.thrust();
  const double transMinor = thrust.thrustMajor();

  //Linearized sphericity calculation (2D)
  double a11 = 0.0;
  double a22 = 0.0;
  double a12 = 0.0;
  double modSum2 = 0.0;

  for (size_t k = 0; k < momenta2.size(); ++k) {
    modSum2 += momenta2[k].mod();
    a11 += momenta2[k].x() * momenta2[k].x() / momenta2[k].mod();
    a22 += momenta2[k].y() * momenta2[k].y() / momenta2[k].mod();
    a12 += momenta2[k].x() * momenta2[k].y() / momenta2[k].mod();
  }

  double trc2 = (a11 + a22) / modSum2;
  double det2 = (a11 * a22 - a12 * a12) / pow(modSum2, 2);

  double eigen21 = (trc2 + sqrt(pow(trc2, 2) - 4 * det2)) / 2;
  double eigen22 = (trc2 - sqrt(pow(trc2, 2) - 4 * det2)) / 2;
  double transSphericity = 2 * eigen22 / (eigen21 + eigen22);

  //Linearized sphericity calculation (3D)
  double b11 = 0.0;
  double b12 = 0.0;
  double b13 = 0.0;
  double b22 = 0.0;
  double b23 = 0.0;
  double b33 = 0.0;
  double modSum3 = 0.0;

  for (size_t k = 0; k < momenta3.size(); ++k) {
    modSum3 += momenta3[k].mod();
    b11 += momenta3[k].x() * momenta3[k].x() / momenta3[k].mod();
    b22 += momenta3[k].y() * momenta3[k].y() / momenta3[k].mod();
    b33 += momenta3[k].z() * momenta3[k].z() / momenta3[k].mod();
    b12 += momenta3[k].x() * momenta3[k].y() / momenta3[k].mod();
    b13 += momenta3[k].x() * momenta3[k].z() / momenta3[k].mod();
    b23 += momenta3[k].y() * momenta3[k].z() / momenta3[k].mod();
  }

  Matrix3 sph3;
  sph3.set(0, 0, b11 / modSum3);
  sph3.set(0, 1, b12 / modSum3);
  sph3.set(0, 2, b13 / modSum3);
  sph3.set(1, 0, b12 / modSum3);
  sph3.set(1, 1, b22 / modSum3);
  sph3.set(1, 2, b23 / modSum3);
  sph3.set(2, 0, b13 / modSum3);
  sph3.set(2, 1, b23 / modSum3);
  sph3.set(2, 2, b33 / modSum3);

  double q = sph3.trace() / 3.;
  double p1 = sph3.get(0, 1) * sph3.get(0, 1) + sph3.get(0, 2) * sph3.get(0, 2)
      + sph3.get(1, 2) * sph3.get(1, 2);
  double p2 = (sph3.get(0, 0) - q) * (sph3.get(0, 0) - q) + (sph3.get(1, 1) - q) * (sph3.get(1, 1) - q)
      + (sph3.get(2, 2) - q) * (sph3.get(2, 2) - q) + 2 * p1;
  double p = sqrt(p2 / 6.);

  Matrix3 I3 = Matrix3::mkIdentity();
  double r = (1. / p * (sph3 - q * I3)).det() / 2.;

  double phi(0);
  if (r <= -1)
    phi = M_PI / 3.;
  else if (r >= 1)
    phi = 0;
  else
    phi = acos(r) / 3.;

  double eigen31 = q + 2 * p * cos(phi);
  double eigen33 = q + 2 * p * cos(phi + (2 * M_PI / 3.));
  double eigen32 = 3 * q - eigen31 - eigen33;

  double aplanarity = (3. / 2) * eigen33;
  double C = 3 * (eigen31 * eigen32 + eigen31 * eigen33 + eigen32 * eigen33);
  double D = 27 * eigen31 * eigen32 * eigen33;

  //Fill event-shape histograms
  if (ht2 > 1000.0 * GeV && ht2 < 1500.0 * GeV) {

    if (goodJets.size() == 3) {
      _h["transThrust_j3_h1"]->fill(transThrust);
      _h["transMinor_j3_h1"]->fill(transMinor);
      _h["transSphericity_j3_h1"]->fill(transSphericity);
      _h["aplanarity_j3_h1"]->fill(aplanarity);
      _h["C_j3_h1"]->fill(C);
      _h["D_j3_h1"]->fill(D);
    }

    if (goodJets.size() == 4) {
      _h["transThrust_j4_h1"]->fill(transThrust);
      _h["transMinor_j4_h1"]->fill(transMinor);
      _h["transSphericity_j4_h1"]->fill(transSphericity);
      _h["aplanarity_j4_h1"]->fill(aplanarity);
      _h["C_j4_h1"]->fill(C);
      _h["D_j4_h1"]->fill(D);
    }

    if (goodJets.size() == 5) {
      _h["transThrust_j5_h1"]->fill(transThrust);
      _h["transMinor_j5_h1"]->fill(transMinor);
      _h["transSphericity_j5_h1"]->fill(transSphericity);
      _h["aplanarity_j5_h1"]->fill(aplanarity);
      _h["C_j5_h1"]->fill(C);
      _h["D_j5_h1"]->fill(D);
    }

    if (goodJets.size() >= 6) {
      _h["transThrust_j6_h1"]->fill(transThrust);
      _h["transMinor_j6_h1"]->fill(transMinor);
      _h["transSphericity_j6_h1"]->fill(transSphericity);
      _h["aplanarity_j6_h1"]->fill(aplanarity);
      _h["C_j6_h1"]->fill(C);
      _h["D_j6_h1"]->fill(D);
    }
  }


  if (ht2 > 1500.0 * GeV && ht2 < 2000.0 * GeV) {

    if (goodJets.size() == 3) {
      _h["transThrust_j3_h2"]->fill(transThrust);
      _h["transMinor_j3_h2"]->fill(transMinor);
      _h["transSphericity_j3_h2"]->fill(transSphericity);
      _h["aplanarity_j3_h2"]->fill(aplanarity);
      _h["C_j3_h2"]->fill(C);
      _h["D_j3_h2"]->fill(D);
    }

    if (goodJets.size() == 4) {
      _h["transThrust_j4_h2"]->fill(transThrust);
      _h["transMinor_j4_h2"]->fill(transMinor);
      _h["transSphericity_j4_h2"]->fill(transSphericity);
      _h["aplanarity_j4_h2"]->fill(aplanarity);
      _h["C_j4_h2"]->fill(C);
      _h["D_j4_h2"]->fill(D);
    }

    if (goodJets.size() == 5) {
      _h["transThrust_j5_h2"]->fill(transThrust);
      _h["transMinor_j5_h2"]->fill(transMinor);
      _h["transSphericity_j5_h2"]->fill(transSphericity);
      _h["aplanarity_j5_h2"]->fill(aplanarity);
      _h["C_j5_h2"]->fill(C);
      _h["D_j5_h2"]->fill(D);
    }

    if (goodJets.size() >= 6) {
      _h["transThrust_j6_h2"]->fill(transThrust);
      _h["transMinor_j6_h2"]->fill(transMinor);
      _h["transSphericity_j6_h2"]->fill(transSphericity);
      _h["aplanarity_j6_h2"]->fill(aplanarity);
      _h["C_j6_h2"]->fill(C);
      _h["D_j6_h2"]->fill(D);
    }
  }

  if (ht2 > 2000.0 * GeV) {

    if (goodJets.size() == 3) {
      _h["transThrust_j3_h3"]->fill(transThrust);
      _h["transMinor_j3_h3"]->fill(transMinor);
      _h["transSphericity_j3_h3"]->fill(transSphericity);
      _h["aplanarity_j3_h3"]->fill(aplanarity);
      _h["C_j3_h3"]->fill(C);
      _h["D_j3_h3"]->fill(D);
    }

    if (goodJets.size() == 4) {
      _h["transThrust_j4_h3"]->fill(transThrust);
      _h["transMinor_j4_h3"]->fill(transMinor);
      _h["transSphericity_j4_h3"]->fill(transSphericity);
      _h["aplanarity_j4_h3"]->fill(aplanarity);
      _h["C_j4_h3"]->fill(C);
      _h["D_j4_h3"]->fill(D);
    }

    if (goodJets.size() == 5) {
      _h["transThrust_j5_h3"]->fill(transThrust);
      _h["transMinor_j5_h3"]->fill(transMinor);
      _h["transSphericity_j5_h3"]->fill(transSphericity);
      _h["aplanarity_j5_h3"]->fill(aplanarity);
      _h["C_j5_h3"]->fill(C);
      _h["D_j5_h3"]->fill(D);
    }

    if (goodJets.size() >= 6) {
      _h["transThrust_j6_h3"]->fill(transThrust);
      _h["transMinor_j6_h3"]->fill(transMinor);
      _h["transSphericity_j6_h3"]->fill(transSphericity);
      _h["aplanarity_j6_h3"]->fill(aplanarity);
      _h["C_j6_h3"]->fill(C);
      _h["D_j6_h3"]->fill(D);
    }
  }
}


void finalize() {

  const double xs1 = _h["njet_h1"]->sumW();
  const double xs2 = _h["njet_h2"]->sumW();
  const double xs3 = _h["njet_h3"]->sumW();

  for (auto& hist : _h) {
    if (hist.first.find("njet_") != string::npos)
      scale(hist.second, crossSectionPerEvent() / picobarn);
    else if (hist.first.find("_h1") != string::npos)
      scale(hist.second, 1.0 / xs1);
    else if (hist.first.find("_h2") != string::npos)
      scale(hist.second, 1.0 / xs2);
    else if (hist.first.find("_h3") != string::npos)
      scale(hist.second, 1.0 / xs3);
  }
}

private:

//Jet multiplicity
map<string, Histo1DPtr> _h;

};

RIVET_DECLARE_PLUGIN(ATLAS_2020_I1808726); } ```