Rivet analyses

title: H1_2002_I561885

Measurement of $D^{*\pm}$ meson production in deep inelastic scattering at HERA

Experiment: H1 (HERA)

Inspire ID: 561885

Status: VALIDATED

Authors: - Madhav Chithirasreemadam - Hannes Jung

References: none listed

Beams: p+ e-, e- p+

Beam energies: (820.0, 27.5); (27.5, 820.0)GeV

Run details: none listed

The inclusive production of $D^{\pm}$ (2010) mesons in deep-inelastic scattering is studied with the H1 detector at HERA. In the kinematic region $1 < Q^2 < 100 $ GeV$^2$ and $0.05 < y < 0.7$ in $e^+ p$ cross section for inclusive $D^{\pm}$ meson production is measured in the visible range $p_{T}^{D} > 1.5$ GeV and $|\eta^{D}| < 1.5 $.

Source code:H1_2002_I561885.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/DISKinematics.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief Measurement of D*+- meson production in deep inelastic scattering at HERA (H1) class H1_2002_I561885 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(H1_2002_I561885);


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

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


  declare(DISKinematics(), "Kinematics");
  declare(UnstableParticles(), "Dstars");
  //Cuts::abspid == PID::DSTARPLUS

  // Initialise and register projections


  Histo1DPtr dummy; //Introducing


  // Book histograms
  book(_h["W_GeV"], 2, 1, 1);
  book(_h["p_tD*"], 3, 1, 1);
  book(_h["logx"], 4, 1, 1);
  book(_h["etaD*"], 5, 1, 1);
  book(_h_Q2, 6, 1, 1);
  book(_h["Z_D*"], 7, 1, 1);

  book(_h_Q2eta, {-1.5, -0.5, 0.5, 1.5}, {"d08-x01-y01", "d08-x01-y02", "d08-x01-y03"});
  book(_h_Q2pt, {1.5, 4., 10.}, {"d09-x01-y01", "d09-x01-y02"});

  book(_h_eta1, {0., 0.25, 0.5, 1.}, {"d10-x01-y01", "d10-x01-y02", "d10-x01-y03"});
  book(_h_eta2, {1.5, 2.5, 4., 10.}, {"d11-x01-y01", "d11-x01-y02", "d11-x01-y03"});
  book(_h_zD1, {1.5, 2.5, 4., 10.}, {"d12-x01-y01", "d12-x01-y02", "d12-x01-y03"});
  _axes["Q2"] = YODA::Axis<double>{1., 2.371, 4.218, 10.0, 17.78, 31.62, 100.0};
  _axes["Q2eta"] = YODA::Axis<double>{1.0, 9.0, 20.0, 100.0};
  _axes["Q2pt"] = YODA::Axis<double>{1.0, 3.0, 8.0, 20.0, 100.0};
}


/// Perform the per-event analysis
void analyze(const Event& event) {
  if (_edges.empty()) {
    _edges["Q2"] = _h_Q2->xEdges();
    _edges["Q2eta"] = _h_Q2eta->bins()[1]->xEdges();
    _edges["Q2pt"] = _h_Q2pt->bins()[1]->xEdges();
  }
  const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");

  // Q2 and inelasticity cuts
  if (!inRange(kin.Q2(), 1.0 * GeV2, 100 * GeV2)) vetoEvent;
  if (!inRange(kin.y(), 0.05, 0.7)) vetoEvent;

  // D* reconstruction
  // const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5*GeV && Cuts::abseta < 1.5);
  const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5 * GeV);
  const Particles dstars = select(unstables,
                                  [](const Particle& p) { return p.abspid() == PID::DSTARPLUS; });
  if (dstars.empty()) vetoEvent;
  // MSG_DEBUG("#D* = " << dstars.size());
  //const Particle& dstar = dstars.front();
  for (const Particle& dstar : dstars) {
    const double zD = (dstar.E() - dstar.pz()) / (2 * kin.beamLepton().E() * kin.y());


    // Single-differential histograms
    _h["p_tD*"]->fill(dstar.pT() / GeV);
    _h["etaD*"]->fill(dstar.eta());
    _h["Z_D*"]->fill(zD / GeV);

    size_t idx = _axes["Q2"].index(kin.Q2() / GeV2);
    string edge = "OTHER";
    if (idx && idx <= _edges["Q2"].size()) edge = _edges["Q2"][idx - 1];
    _h_Q2->fill(edge);

    _h["W_GeV"]->fill(sqrt(kin.W2()) / GeV);
    _h["logx"]->fill(log10(kin.x()));

    // Double-differential (y,Q2) histograms

    idx = _axes["Q2eta"].index(kin.Q2() / GeV2);
    edge = "OTHER";
    if (idx && idx <= _edges["Q2eta"].size()) edge = _edges["Q2eta"][idx - 1];
    _h_Q2eta->fill(dstar.eta(), edge);
    idx = _axes["Q2pt"].index(kin.Q2() / GeV2);
    edge = "OTHER";
    if (idx && idx <= _edges["Q2pt"].size()) edge = _edges["Q2pt"][idx - 1];
    _h_Q2pt->fill(dstar.pT(), edge);
    _h_eta1->fill(zD, dstar.eta());
    _h_eta2->fill(dstar.pT(), dstar.eta());
    _h_zD1->fill(dstar.pT(), zD / GeV);
  }
}


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


  const double sf = crossSection() / nanobarn / sumOfWeights();
  scale(_h["p_tD*"], sf);
  scale(_h["etaD*"], sf);
  scale(_h["Z_D*"], sf);
  scale(_h_Q2, sf);
  for (auto& b : _h_Q2->bins()) {
    const size_t idx = b.index();
    b.scaleW(1. / _axes["Q2"].width(idx));
  }
  scale(_h["W_GeV"], sf);
  scale(_h["logx"], sf);

  scale(_h_Q2eta, sf);
  for (auto& histo : _h_Q2eta->bins()) {
    for (auto& b : histo->bins()) {
      const size_t idx = b.index();
      b.scaleW(1. / _axes["Q2eta"].width(idx));
    }
  }
  scale(_h_Q2pt, sf);
  for (auto& histo : _h_Q2pt->bins()) {
    for (auto& b : histo->bins()) {
      const size_t idx = b.index();
      b.scaleW(1. / _axes["Q2pt"].width(idx));
    }
  }
  scale(_h_eta1, sf);
  scale(_h_eta2, sf);
  scale(_h_zD1, sf);
  // distributionsa are really double differential divide by group width
  _h_Q2eta->divByGroupWidth();
  _h_Q2pt->divByGroupWidth();
  _h_eta1->divByGroupWidth();
  _h_eta2->divByGroupWidth();
  _h_zD1->divByGroupWidth();
}

///@}


/// @name Histograms
///@{
map<string, Histo1DPtr> _h;
map<string, Profile1DPtr> _p;
map<string, CounterPtr> _c;
HistoGroupPtr<double, string> _h_Q2eta, _h_Q2pt;
Histo1DGroupPtr _h_eta1, _h_eta2, _h_zD1;
BinnedHistoPtr<string> _h_Q2;
map<string, YODA::Axis<double>> _axes;
map<string, vector<string>> _edges;
///@}

};

RIVET_DECLARE_PLUGIN(H1_2002_I561885);

} ```