Rivet analyses

title: H1_2015_I1343110

Diffractive dijets in DIS and photoproduction

Experiment: H1 (HERA)

Inspire ID: 1343110

Status: UNVALIDATED

Authors: - Christine O. Rasmussen - Ilkka Helenius

References: - JHEP 05 (2015) 056 - arXiv: 1502.01683

Beams: p+ e+

Beam energies: (920.0, 27.5)GeV

Run details: - 920 GeV protons colliding with 27.5 GeV positrons; Diffractive photoproduction of dijets; Tagged protons.; Jet $pT > 5$ GeV; Note that CM energy is WRONG in HepData table.; Has been changed by hand in YODA refdata file to correct value.

H1 diffractive jets from proton--positron collisions at beam energies of 920~GeV and 27.5~GeV. The cross section of the diffractive process e+p -> e+X+p is measured at a centre-of-mass energy of 318 GeV, where the system X contains at least two jets and the leading final state proton p is detected in the H1 Very Forward Proton Spectrometer. The measurement is performed in photoproduction with photon virtualities Q^2 < 2 GeV^2 and in deep-inelastic scattering with 4 GeV^2 < Q2 < 80 GeV^2.

Source code:H1_2015_I1343110.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/DISDiffHadron.hh"

include "Rivet/Projections/DISFinalState.hh"

include "Rivet/Projections/DISKinematics.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

namespace Rivet {

/// @brief H1 diffractive dijets /// /// Diffractive dijets H1 with 920 GeV p and 27.5 GeV e /// Tagged protons & jets found in gamma*p rest frame. /// /// @author Christine O. Rasmussen class H1_2015_I1343110 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(H1_2015_I1343110);

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

// Book projections and histograms
void init() {

  declare(DISKinematics(), "Kinematics");
  declare(DISFinalState(DISFrame::LAB), "FSLab");
  declare(FastJets(DISFinalState(DISFrame::HCM), fastjet::JetAlgorithm::kt_algorithm,
                   fastjet::RecombinationScheme::pt_scheme, 1.0, JetMuons::ALL, JetInvisibles::NONE,
                   nullptr),
          "DISFSJets");
  declare(DISDiffHadron(), "Hadron");

  // Book histograms from REF data
  book(_h_PHO_sig_sqrts, 1, 1, 1);
  book(_h_DIS_sig_sqrts, 2, 1, 1);
  book(_h_PHODIS_sqrts, 3, 1, 1);

  book(_h_DIS_dsigdz, 4, 1, 1);
  book(_h_DIS_dsigdxPom, 5, 1, 1);
  book(_h_DIS_dsigdy, 6, 1, 1);
  book(_h_DIS_dsigdQ2, 7, 1, 1);
  book(_h_DIS_dsigdEtj1, 8, 1, 1);
  book(_h_DIS_dsigdMX, 9, 1, 1);
  book(_h_DIS_dsigdDeltaEta, 10, 1, 1);
  book(_h_DIS_dsigdAvgEta, 11, 1, 1);

  book(_h_PHO_dsigdz, 12, 1, 1);
  book(_h_PHO_dsigdxPom, 13, 1, 1);
  book(_h_PHO_dsigdy, 14, 1, 1);
  book(_h_PHO_dsigdxGam, 15, 1, 1);
  book(_h_PHO_dsigdEtj1, 16, 1, 1);
  book(_h_PHO_dsigdMX, 17, 1, 1);
  book(_h_PHO_dsigdDeltaEta, 18, 1, 1);
  book(_h_PHO_dsigdAvgEta, 19, 1, 1);

  book(_h_PHODIS_deltaEta, 20, 1, 1);
  book(_h_PHODIS_y, 21, 1, 1);
  book(_h_PHODIS_z, 22, 1, 1);
  book(_h_PHODIS_Etj1, 23, 1, 1);

  isPHO = false;
  nVeto0 = 0;
  nVeto1 = 0;
  nVeto2 = 0;
  nVeto3 = 0;
  nVeto4 = 0;
  nVeto5 = 0;
  nVeto6 = 0;
  nPHO = 0;
  nDIS = 0;
}

// Do the analysis
void analyze(const Event& event) {

  // Event weight
  isPHO = false;

  // Projections - special handling of events where no proton found:
  const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
  if (kin.failed()) vetoEvent;
  Particle hadronOut;
  Particle hadronIn;
  try {
    const DISDiffHadron& diffhadr = apply<DISDiffHadron>(event, "Hadron");
    hadronOut = diffhadr.out();
    hadronIn = diffhadr.in();
  }
  catch (const Error& e) {
    vetoEvent;
  }
  ++nVeto0;

  // Determine kinematics: H1 has +z = proton direction
  int dir = kin.orientation();
  double y = kin.y();
  double Q2 = kin.Q2();

  // Separate into DIS and PHO regimes else veto
  if (Q2 < 2. * GeV2 && inRange(y, 0.2, 0.70)) {
    isPHO = true;
    ++nPHO;
  }
  else if (inRange(Q2, 4.0 * GeV2, 80. * GeV2) && inRange(y, 0.2, 0.7)) {
    isPHO = false;
    ++nDIS;
  }
  else
    vetoEvent;
  ++nVeto1;

  // Find diffractive variables as defined in paper.
  // Note tagged protons in VFPS => smaller allowed xPom range
  // xPom = 1 - E'/E, t = (P-P')^2
  const double abst = abs((hadronIn.mom() - hadronOut.mom()).mass2());
  const double xPom = 1. - hadronOut.energy() / hadronIn.energy();
  // Veto if outside allowed region
  if (abst > 0.6 * GeV2) vetoEvent;
  ++nVeto2;
  if (!inRange(xPom, 0.010, 0.024)) vetoEvent;
  ++nVeto3;

  // Jet selection. Note jets are found in HCM frame, but
  // eta cut is applied in lab frame!
  Jets jets = apply<FastJets>(event, "DISFSJets").jets(Cuts::Et > 4. * GeV, cmpMomByEt);
  // Veto if not dijets and if Et_j1 < 5.5
  if (jets.size() < 2) vetoEvent;
  if (jets[0].Et() < 5.5 * GeV) vetoEvent;
  ++nVeto4;

  const double EtJet1 = jets[0].Et();

  // Transform from HCM to LAB and apply eta cut
  const LorentzTransform boost = kin.boostHCM();
  for (int i = 0; i < 2; ++i) {
    jets[i].transformBy(boost.inverse());
    if (!inRange(dir * jets[i].eta(), -1, 2.5)) vetoEvent;
  }
  ++nVeto5;

  double sumE(0.), sumPz(0.);
  FourMomentum tot_mom;
  for (const Particle& part : apply<FinalState>(event, "FSLab").particles()) {
    if (part.genParticle() != hadronOut.genParticle()) {
      sumE += part.energy();
      sumPz += part.pz();
      tot_mom += part.mom();
    }
  }
  double EplusPz = sumE + dir * sumPz;
  double EminusPz = sumE - dir * sumPz;
  double M2X = tot_mom.mass2();

  // Pseudorapidity distributions are examined in lab frame:
  double deltaEtaJets = abs(dir * jets[0].eta() - dir * jets[1].eta());
  double avgEtaJets = 0.5 * (dir * jets[0].eta() + dir * jets[1].eta());

  // Evaluate observables
  double zPomJets(0.), xGamJets(0.);
  FourMomentum pj = jets[0].momentum() + jets[1].momentum();
  if (isPHO) {
    zPomJets = (pj.E() + dir * pj.pz()) / EplusPz;
    xGamJets = (pj.E() - dir * pj.pz()) / EminusPz;
  }
  else {
    zPomJets = (Q2 + pj.mass2()) / (Q2 + M2X);
  }
  // Veto events with zPom > 0.8
  if (zPomJets > 0.8) vetoEvent;
  ++nVeto6;

  // Now fill histograms
  if (isPHO) {
    _h_PHO_sig_sqrts->fill(sqrtS() / GeV);
    _h_PHO_dsigdz->fill(zPomJets);
    _h_PHO_dsigdxPom->fill(xPom);
    _h_PHO_dsigdy->fill(y);
    _h_PHO_dsigdxGam->fill(xGamJets);
    _h_PHO_dsigdEtj1->fill(EtJet1);
    _h_PHO_dsigdMX->fill(sqrt(M2X) / GeV);
    _h_PHO_dsigdDeltaEta->fill(deltaEtaJets);
    _h_PHO_dsigdAvgEta->fill(avgEtaJets);
  }
  else {
    _h_DIS_sig_sqrts->fill(sqrtS() / GeV);
    _h_DIS_dsigdz->fill(zPomJets);
    _h_DIS_dsigdxPom->fill(xPom);
    _h_DIS_dsigdy->fill(y);
    _h_DIS_dsigdQ2->fill(Q2);
    _h_DIS_dsigdEtj1->fill(EtJet1);
    _h_DIS_dsigdMX->fill(sqrt(M2X) / GeV);
    _h_DIS_dsigdDeltaEta->fill(deltaEtaJets);
    _h_DIS_dsigdAvgEta->fill(avgEtaJets);
  }
}

// Finalize
void finalize() {
  // Normalise to cross section
  // Remember to manually scale the cross section afterwards with
  // the number of rejected events.
  const double norm = crossSection() / picobarn / sumOfWeights();

  scale(_h_PHO_sig_sqrts, norm);
  scale(_h_PHO_dsigdz, norm);
  scale(_h_PHO_dsigdxPom, norm);
  scale(_h_PHO_dsigdy, norm);
  scale(_h_PHO_dsigdxGam, norm);
  scale(_h_PHO_dsigdEtj1, norm);
  scale(_h_PHO_dsigdMX, norm);
  scale(_h_PHO_dsigdDeltaEta, norm);
  scale(_h_PHO_dsigdAvgEta, norm);

  scale(_h_DIS_sig_sqrts, norm);
  scale(_h_DIS_dsigdz, norm);
  scale(_h_DIS_dsigdxPom, norm);
  scale(_h_DIS_dsigdy, norm);
  scale(_h_DIS_dsigdQ2, norm);
  scale(_h_DIS_dsigdEtj1, norm);
  scale(_h_DIS_dsigdMX, norm);
  scale(_h_DIS_dsigdDeltaEta, norm);
  scale(_h_DIS_dsigdAvgEta, norm);

  if (_h_DIS_sig_sqrts->numEntries() != 0) divide(_h_PHO_sig_sqrts, _h_DIS_sig_sqrts, _h_PHODIS_sqrts);
  if (_h_DIS_dsigdDeltaEta->numEntries() != 0)
    divide(_h_PHO_dsigdDeltaEta, _h_DIS_dsigdDeltaEta, _h_PHODIS_deltaEta);
  if (_h_DIS_dsigdy->numEntries() != 0) divide(_h_PHO_dsigdy, _h_DIS_dsigdy, _h_PHODIS_y);
  if (_h_DIS_dsigdz->numEntries() != 0) divide(_h_PHO_dsigdz, _h_DIS_dsigdz, _h_PHODIS_z);
  if (_h_DIS_dsigdEtj1->numEntries() != 0) divide(_h_PHO_dsigdEtj1, _h_DIS_dsigdEtj1, _h_PHODIS_Etj1);

  const double dPHO = numEvents();
  MSG_INFO("H1_2015_I1343110");
  MSG_INFO("Cross section = " << crossSection() / picobarn << " pb");
  MSG_INFO("Number of events = " << numEvents() << ", sumW = " << sumOfWeights());
  MSG_INFO("Number of PHO = " << nPHO << ", number of DIS = " << nDIS);
  MSG_INFO("Events passing event setup     = " << nVeto0 << " (" << nVeto0 / dPHO * 100. << "%)");
  MSG_INFO("Events passing electron veto   = " << nVeto1 << " (" << nVeto1 / dPHO * 100. << "%)");
  MSG_INFO("Events passing t veto          = " << nVeto2 << " (" << nVeto2 / dPHO * 100. << "%)");
  MSG_INFO("Events passing xPom            = " << nVeto3 << " (" << nVeto3 / dPHO * 100. << "%)");
  MSG_INFO("Events passing jet Et   veto   = " << nVeto4 << " (" << nVeto4 / dPHO * 100. << "%)");
  MSG_INFO("Events passing jet eta veto    = " << nVeto5 << " (" << nVeto5 / dPHO * 100. << "%)");
  MSG_INFO("Events passing zPom veto       = " << nVeto6 << " (" << nVeto6 / dPHO * 100. << "%)");
}

/// @}

private:

/// @name Histograms
/// @{
// Book histograms from REF data
BinnedHistoPtr<int> _h_PHO_sig_sqrts;
BinnedHistoPtr<int> _h_DIS_sig_sqrts;
BinnedEstimatePtr<int> _h_PHODIS_sqrts;

Histo1DPtr _h_DIS_dsigdz;
Histo1DPtr _h_DIS_dsigdxPom;
Histo1DPtr _h_DIS_dsigdy;
Histo1DPtr _h_DIS_dsigdQ2;
Histo1DPtr _h_DIS_dsigdEtj1;
Histo1DPtr _h_DIS_dsigdMX;
Histo1DPtr _h_DIS_dsigdDeltaEta;
Histo1DPtr _h_DIS_dsigdAvgEta;

Histo1DPtr _h_PHO_dsigdz;
Histo1DPtr _h_PHO_dsigdxPom;
Histo1DPtr _h_PHO_dsigdy;
Histo1DPtr _h_PHO_dsigdxGam;
Histo1DPtr _h_PHO_dsigdEtj1;
Histo1DPtr _h_PHO_dsigdMX;
Histo1DPtr _h_PHO_dsigdDeltaEta;
Histo1DPtr _h_PHO_dsigdAvgEta;

Estimate1DPtr _h_PHODIS_deltaEta;
Estimate1DPtr _h_PHODIS_y;
Estimate1DPtr _h_PHODIS_z;
Estimate1DPtr _h_PHODIS_Etj1;
/// @}

bool isPHO;
int nVeto0, nVeto1, nVeto2, nVeto3, nVeto4, nVeto5, nVeto6;
int nPHO, nDIS;

};

RIVET_DECLARE_PLUGIN(H1_2015_I1343110);

} ```