Rivet analyses

title: CMS_2014_I1305624

Study of hadronic event-shape variables in multijet final states in $pp$ collisions at $\sqrt{s} = 7$ TeV

Experiment: CMS (LHC)

Inspire ID: 1305624

Status: VALIDATED

Authors: - Sunanda Banerjee - Sandeep Bhowmik - Monoranjan Guchait - Gobinda Majumder - Manas Maity < - Debarati Roy

References: - 10.1007/JHEP10(2014)087 - arXiv: 1407.2856 - Expt page: CMS-PAS-SMP-12-022 - CERN-PH-EP-2014-146

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - $pp$ QCD interactions at $\sqrt{s} = 7$ TeV. Data collected by CMS during the year 2011.

Event-shape variables, which are sensitive to perturbative and nonperturbative aspects of quantum chromodynamic (QCD) interactions, are studied in multijet events recorded in proton-proton collisions at $\sqrt{s}=7$ TeV. Events are selected with at least one jet with transverse momentum $p_T > 110$ GeV and pseudorapidity $\mid\eta\mid < 2.4$, in a data sample corresponding to integrated luminosities of up to $5 fb^{-1}$. The distributions of five event-shape variables in various leading jet $p_T$ ranges are compared to predictions from different QCD Monte Carlo event generators. Five event-shape variables are analyzed in this paper: the transverse thrust $\tau_{\perp}$, the total jet broadening $B_{tot}$, the total jet mass $\rho_{tot}$, the total transverse jet mass $\rho^T_{tot}$ and the third-jet resolution parameter $Y_{23}$. In the formulae below, $p_{T,i}$, $\eta_i$, and $\phi_i$ represent the transverse momentum, pseudorapidity, and azimuthal angle of the $i$th jet, and $\hat{n}{T}$ is the unit vector that maximizes the sum of the projections of $\vec{p}{T,i}$. The transverse thrust axis $\hat{n}{T}$ and the beam form the so-called event plane. Based on the direction of $\hat{n}{T}$, the transverse region is separated into an upper side $\cal{C}U$, consisting of all jets with $\vec{p}{T}\cdot\hat{n}{T}$ $>$ 0, and a lower side $\cal{C}_L$, with $\vec{p}{T}\cdot\hat{n}{T}$ $<$ 0. The jet broadening and third-jet resolution variables require at least three jets, whereas the calculation of other variables requires at least two jets. The $\hat{n}{T}$ vector is defined only up to a global sign - choosing one sign or the other has no consequence since it simply exchanges the upper and lower events regions. Transverse Thrust : The event thrust observable in the transverse plane is defined as \begin{eqnarray} \tau_{\perp} & \equiv & 1 - \max_{\hat{n}{T}} \frac{\sum_i |\vec{p}{T,i} \cdot \hat{n}{T} |}{\sum_i p{T,i}} \end{eqnarray} This variable probes the hadronization process and is sensitive to the modeling of two-jet and multijet topologies. In this paper multijet' refers tomore-than-two-jet'. In the limit of a perfectly balanced two-jet event, $\tau_{\perp}$ is zero, while in isotropic multijet events it amounts to $(1-2/\pi)$. Jet Broadenings : The pseudorapidities and the azimuthal angles of the axes for the upper and lower event regions are defined by \begin{eqnarray} \eta_X & \equiv & \frac{\sum_{i\in{\cal{C}}X} p{T,i}\eta_i}{\sum_{i\in{\cal{C}}X} p{T,i}} ~, \ \phi_X & \equiv & \frac{\sum_{i\in{\cal{C}}X} p{T,i}\phi_i}{\sum_{i\in{\cal{C}}X} p{T,i}} ~, \end{eqnarray} where $X$ refers to upper ($U$) or lower ($L$) side. From these, the jet broadening variable in each region is defined as \begin{eqnarray} B_{X} & \equiv & \frac{1}{2P_{T}} \sum_{i\in{\cal{C}}X}p{T,i}\sqrt{(\eta_i -\eta_X)^2 + (\phi_i - \phi_X)^2} ~, \end{eqnarray} where $P_{T}$ is the scalar sum of the transverse momenta of all the jets. The total jet broadening is then defined as \begin{eqnarray} B_{tot} & \equiv & B_{U} + B_{L}. \end{eqnarray} Jet Masses : The normalized squared invariant mass of the jets in the upper and lower regions of the event is defined by \begin{eqnarray} \rho_X & \equiv & \frac{M^2_X}{P^2}, \end{eqnarray} where $M_X$ is the invariant mass of the constituents of the jets in the region $X$, and $P$ is the scalar sum of the momenta of all constituents in both sides. The jet mass variable is defined as the sum of the masses in the upper and lower regions, \begin{eqnarray} \rho_{tot} & \equiv & \rho_U + \rho_L ~. \end{eqnarray} The corresponding jet mass in the transverse plane, $\rho^T_{tot}$, is also similarly calculated in transverse plane. Third-jet resolution parameter : The third-jet resolution parameter is defined as \begin{eqnarray} Y_{23} \equiv \frac{\mathrm{min}(p_{T,3}^2,[\mathrm{min}(p_{T,i}, p_{T,j})^2 \times (\Delta R_{ij})^2/R^2])}{P_{12}^2} ~, \end{eqnarray} where i, j run over all three jets, $(\Delta R_{ij})^2 = (\eta_i - \eta_j)^2 + (\phi_i - \phi_j)^2$, and $p_{T,3}$ is the transverse momentum of the third jet in the event. If there are more than three jets in the event, they are iteratively merged using the $k_T$ algorithm with a distance parameter $R = 0.6$. To compute $P_{12}$, three jets are merged into two using the procedure described above and $P_{12}$ is then defined as the scalar sum of the transverse momenta of the two remaining jets. The $Y_{23}$ variable estimates the relative strength of the $p_T$ of the third jet with respect to the other two jets. It vanishes for two-jet events, and a nonzero value of $Y_{23}$ indicates the presence of hard parton emission, which tests the parton showering model of QCD event generators. A test like this is less sensitive to the details of the underlying event (UE) and parton hadronization models than the other event-shape variables.

Source code:`CMS_2014_I1305624.cc`

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

namespace Rivet {

namespace {

/// Number of event shape variables
/// @todo Move into the EventShape class
const int NEVTVAR = 5;

/// Number of leading jet pT thresholds
/// @todo Move into the analysis class
const int NJETPTMN = 5;
/// Leading jet pT thresholds
/// @todo Move into the analysis class
const double LEADINGPTTHRESHOLD[NJETPTMN] = {110.0, 170.0, 250.0, 320.0, 390.0};


// Helpers for event shape calculations in hidden namespace; implementation at bottom of file
/// @todo Why a class? Improve/remove this junk
class EventShape {
public:

  /// Constructor from vectors of four-vectors as input objects in the event to calculate the event shapes
  EventShape(const vector<double>& px_vector,
             const vector<double>& py_vector,
             const vector<double>& pz_vector,
             const vector<double>& e_vector,
             double eta_central,
             int irap,
             int nmn)
      : _object_px(px_vector),
        _object_py(py_vector),
        _object_pz(pz_vector),
        _object_e(e_vector),
        _eta_c(eta_central),
        _irap(irap),
        _nmnjet(nmn) { }

  /// @brief Returns the values of the five event shapes
  ///
  /// Event shape indices:
  /// 0. central transverse thrust
  /// 1. central total jet broadening
  /// 2. central total jet mass
  /// 3. central total transverse jet mass
  /// 4. central three-jet resolution threshold
  vector<double> getEventShapes() {
    _calculate(); ///< @todo There should be some test for success/failure!!
    return _event_shapes;
  }

  /// Returns the global thrust axis Nx, Ny, Nz=0
  vector<double> getThrustAxis() {
    _calculate(); ///< @todo There should be some test for success/failure!!
    return _thrust_axis;
  }

  /// Returns the central thrust axis Nx, Ny, Nz=0
  vector<double> getThrustAxisC() {
    _calculate(); ///< @todo There should be some test for success/failure!!
    return _thrust_axis_c;
  }

  // /// @brief Choice of the central region
  // void setEtaC(double eta_central) { _eta_c = eta_central; }

  // // Whether to use the rapidity y (rap==1)  or the pseudorapidity eta (rap==0)
  // void setRapType(int irap) { _irap = irap; }


private:

  /// Calculate everything
  int _calculate();

  /// Returns the difference in phi between two vectors
  double _delta_phi(double, double);

  /// The Lorentz scalar product
  double _lorentz_sp(const vector<double>&, const vector<double>&);

  // Calculates the three-jet resolutions
  double _three_jet_res(const vector<double>&,
                        const vector<double>&,
                        const vector<double>&,
                        const vector<double>&,
                        int);

  // Calculates the thrust axis and the tau values
  vector<double> _thrust(const vector<double>&, const vector<double>&);


  vector<double> _object_px, _object_py, _object_pz, _object_p;
  vector<double> _object_pt, _object_e, _object_phi, _object_eta;
  vector<double> _event_shapes;
  vector<double> _thrust_axis, _thrust_axis_c;

  double _eta_c;
  int _irap;
  size_t _nmnjet;
};

}

class CMS_2014_I1305624 : public Analysis { public:

/// Constructor
CMS_2014_I1305624()
    : Analysis("CMS_2014_I1305624") { }


/// @name Analysis methods

/// Book histograms and initialise projections before the run
void init() {
  const FastJets jets(FinalState(Cuts::abseta < 2.6), JetAlg::ANTIKT, 0.5);
  declare(jets, "Jets");

  for (int ij = 0; ij < NJETPTMN; ij++) {
    book(_h_thrustc[ij], 1, 1, ij + 1);
    book(_h_broadt[ij], 1, 2, ij + 1);
    book(_h_tot3dmass[ij], 1, 3, ij + 1);
    book(_h_tottrnsmass[ij], 1, 4, ij + 1);
    book(_h_y23c[ij], 1, 5, ij + 1);
    //
  }
  _needBinInit = true;
}


/// Perform the per-event analysis
void analyze(const Event& event) {

  if (_needBinInit) {
    for (int ij = 0; ij < NJETPTMN; ij++) {
      _alow1[ij] = _h_thrustc[ij]->xMin();
      _alow2[ij] = _h_broadt[ij]->xMin();
      _alow3[ij] = _h_tot3dmass[ij]->xMin();
      _alow4[ij] = _h_tottrnsmass[ij]->xMin();
      _alow5[ij] = _h_y23c[ij]->xMin();
      //
      _ahgh1[ij] = _h_thrustc[ij]->xMax();
      _ahgh2[ij] = _h_broadt[ij]->xMax();
      _ahgh3[ij] = _h_tot3dmass[ij]->xMax();
      _ahgh4[ij] = _h_tottrnsmass[ij]->xMax();
      _ahgh5[ij] = _h_y23c[ij]->xMax();

      _needBinInit = false;
    }
  }

  const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 30 * GeV);
  if (jets.size() < 2) vetoEvent;
  if (jets[0].abseta() > 2.4 || jets[1].abseta() > 2.4) vetoEvent;

  const double leadingpt = jets[0].pT();
  if (leadingpt < 110 * GeV) vetoEvent;

  vector<double> jtpx, jtpy, jtpz, jten;
  for (const Jet& j : jets) {
    if (j.abseta() < 2.4) {
      jtpx.push_back(j.px());
      jtpy.push_back(j.py());
      jtpz.push_back(j.pz());
      jten.push_back(j.E());
    }
  }

  EventShape eventshape(jtpx, jtpy, jtpz, jten, 2.4, 0, 2);
  const vector<double> eventvar = eventshape.getEventShapes();
  if (eventvar[NEVTVAR] < 0) vetoEvent; // Jets are not only one hemisphere

  for (int ij = NJETPTMN - 1; ij >= 0; --ij) {
    if (leadingpt / GeV > LEADINGPTTHRESHOLD[ij]) {
      if (inRange(eventvar[0], _alow1[ij], _ahgh1[ij])) _h_thrustc[ij]->fill(eventvar[0]);
      if (inRange(eventvar[2], _alow3[ij], _ahgh3[ij])) _h_tot3dmass[ij]->fill(eventvar[2]);
      if (inRange(eventvar[3], _alow4[ij], _ahgh4[ij])) _h_tottrnsmass[ij]->fill(eventvar[3]);
      if (eventvar[NEVTVAR] >= 3) {
        if (inRange(eventvar[1], _alow2[ij], _ahgh2[ij])) _h_broadt[ij]->fill(eventvar[1]);
        if (inRange(eventvar[4], _alow5[ij], _ahgh5[ij])) _h_y23c[ij]->fill(eventvar[4]);
      }
      break;
    }
  }
}


/// Normalise histograms etc., after the run
void finalize() {
  for (int ij = 0; ij < NJETPTMN; ij++) {
    normalize(_h_thrustc[ij]);
    normalize(_h_broadt[ij]);
    normalize(_h_tot3dmass[ij]);
    normalize(_h_tottrnsmass[ij]);
    normalize(_h_y23c[ij]);
  }
}

private:

/// @name Histograms
/// @{
Histo1DPtr _h_thrustc[NJETPTMN];
Histo1DPtr _h_broadt[NJETPTMN];
Histo1DPtr _h_tot3dmass[NJETPTMN];
Histo1DPtr _h_tottrnsmass[NJETPTMN];
Histo1DPtr _h_y23c[NJETPTMN];
/// @}

// Data members
bool _needBinInit;
double _alow1[NJETPTMN], _alow2[NJETPTMN], _alow3[NJETPTMN], _alow4[NJETPTMN], _alow5[NJETPTMN];
double _ahgh1[NJETPTMN], _ahgh2[NJETPTMN], _ahgh3[NJETPTMN], _ahgh4[NJETPTMN], _ahgh5[NJETPTMN];

};

RIVET_DECLARE_PLUGIN(CMS_2014_I1305624);

/////////////////////