Rivet analyses

title: H1_2006_I699835

Measurement of Event Shape Variables in Deep-Inelastic Scattering at HERA

Experiment: H1 (HERA)

Inspire ID: 699835

Status: VALIDATED

Authors: - Alejandro Basilio Galvan - Hannes Jung

References: - Eur.Phys.J.C 46 (2006) 343-356 - DOI:10.1140/epjc/s2006-02493-x - arXiv: hep-ex/0512014v1

Beams: e+ p+

Beam energies: (27.6, 820.0); (27.6, 920.0)GeV

Run details: none listed

Deep-inelastic ep scattering data taken with the H1 detector at HERA and corresponding to an integrated luminosity of $106 pb^{-1}$ are used to study the differential distributions of event shape variables. These include thrust, jet broadening, jet mass and the C-parameter. The four-momentum transfer Q is taken to be the relevant energy scale and ranges between 14 GeV and 200 GeV.

Source code:H1_2006_I699835.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/DISFinalState.hh"

include "Rivet/Projections/DISKinematics.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/PromptFinalState.hh"

include "Rivet/Projections/Thrust.hh"

namespace Rivet {

/// @brief Event-shape variables in deep-inelastic scattering at HERA class H1_2006_I699835 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(H1_2006_I699835);


/// @name Analysis methods
///@{

/// Book histograms and initialise projections before the run
void init() {

  // Initialise and register projections

  // The basic final-state projection:
  // all final-state particles within
  // the given eta acceptance
  const FinalState fs(Cuts::abseta < 4.9);
  declare(fs, "FS");

  const DISFinalState DISfs(DISFrame::BREIT);

  const FinalState DISfsCut(DISfs, Cuts::eta < 0);

  declare(Thrust(DISfsCut), "ThrustCut");

  declare(DISKinematics(), "Kinematics");

  //Book histograms for different Q ranges:
  book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");

  const vector<double> QEdges{14., 16., 20., 30., 50., 70., 100., 200.};
  book(_h_tauc, QEdges);
  book(_h_tau, QEdges);
  book(_h_B, QEdges);
  book(_h_rho, QEdges);
  for (size_t iQ = 0; iQ < QEdges.size() - 1; ++iQ) {
    book(_h_tauc->bin(iQ + 1), 1 + iQ, 1, 1);
    book(_h_tau->bin(iQ + 1), 8 + iQ, 1, 1);
    book(_h_B->bin(iQ + 1), 15 + iQ, 1, 1);
    book(_h_rho->bin(iQ + 1), 22 + iQ, 1, 1);
    book(_Nevt_after_cuts_Q[iQ], "TMP/Nevt_after_cuts_Q" + to_string(iQ));
  }
}


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

  const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");

  // The kinematic region covered by the analysis is defined by ranges of Q² and y
  if (dk.Q2() < 196 or dk.Q2() > 40000 or dk.y() < 0.1 or dk.y() > 0.7) {
    vetoEvent;
  }

  const double Q = sqrt(dk.Q2());
  _Nevt_after_cuts->fill();
  const size_t iQ = _h_tauc->binAt(Q).index();
  if (0 < iQ && iQ < 8) _Nevt_after_cuts_Q[iQ - 1]->fill();


  // Boost to hadronic Breit frame
  const LorentzTransform breitboost = dk.boostBreit();

  const FinalState& fs = apply<FinalState>(event, "FS");

  /*Calculate event shape variables:
  thrust_num is \sum |pz_h|
  thrust_den is \sum |p_h|
  b_num is \sum |pt_h|
  sumE is the sum of energies
  (All sums run through the particles in the current hemisphere)
  */
  double thrust_num, thrust_den, b_num, sumE;
  thrust_num = thrust_den = b_num = sumE = 0;
  Vector3 sumMom;

  for (const Particle& p : fs.particles()) {
    // Boost to Breit frame
    const FourMomentum breitMom = breitboost.transform(p.momentum());
    if (breitMom.eta() < 0) {
      thrust_num += abs(breitMom.pz());
      thrust_den += breitMom.p();
      b_num += abs(breitMom.pt());
      sumMom.operator+=(breitMom.p3());
      sumE += breitMom.E();
    }
  }

  //The energy in the current hemisphere must exceed a certain value.
  if (sumE <= Q / 10.0) vetoEvent;

  /* Comment from A. Galvan:
  Thrust here is with respect to the z axis (tau in the paper), while the
  one from Rivet projection is with respect to the maximum thrust
  axis (1 - tau_c in the paper)
  */

  double thrust_mine = 1.0 - ((double)thrust_num) / ((double)thrust_den);
  double b_mine = ((double)b_num) / (2.0 * (double)thrust_den);
  double rho_num = thrust_den * thrust_den - sumMom.dot(sumMom);
  double rho_mine = ((double)rho_num) / (4.0 * (double)thrust_den * (double)thrust_den);

  const Thrust& thrCut = apply<Thrust>(event, "ThrustCut");

  //Fill histograms:
  _h_tauc->fill(Q, 1.0 - thrCut.thrust());
  _h_tau->fill(Q, thrust_mine);
  _h_B->fill(Q, b_mine);
  _h_rho->fill(Q, rho_mine);

  /* Comment from A. Galvan:
   As for the C-parameter, my results did not fit the reference data at
   all. The formula given in the paper is a bit ambiguous, because there is
   a sum that runs through all pairs of particles h,h' and I was not sure
   whether each pair should be counted twice (h,h' and h',h) or not, or if
   the pair of a particle with itself (hh) should be considered.
   That is why I tried all of the possibilities, but none of them worked.
  */
}


/// Normalise histograms etc., after the run
void finalize() {

  // need to multiply by bin widths
  for (size_t iQ = 0; iQ < _h_tauc->numBins(); ++iQ) {
    const double Nev = dbl(*_Nevt_after_cuts_Q[iQ]);
    if (Nev != 0) {
      scale(_h_tauc->bin(iQ + 1), 1. / Nev);
      for (auto& bin : _h_tauc->bin(iQ + 1)->bins()) bin.scaleW(bin.xWidth());
      scale(_h_tau->bin(iQ + 1), 1. / Nev);
      for (auto& bin : _h_tau->bin(iQ + 1)->bins()) bin.scaleW(bin.xWidth());
      scale(_h_B->bin(iQ + 1), 1. / Nev);
      for (auto& bin : _h_B->bin(iQ + 1)->bins()) bin.scaleW(bin.xWidth());
      scale(_h_rho->bin(iQ + 1), 1. / Nev);
      for (auto& bin : _h_rho->bin(iQ + 1)->bins()) bin.scaleW(bin.xWidth());
    }
  }
}

///@}


/// @name Histograms
///@{
Histo1DGroupPtr _h_tauc, _h_tau, _h_B, _h_rho;
CounterPtr _Nevt_after_cuts;
CounterPtr _Nevt_after_cuts_Q[7];
///@}

};

RIVET_DECLARE_PLUGIN(H1_2006_I699835);

} ```