Rivet analyses

title: H1_2007_I736052

Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1)

Experiment: H1 (HERA)

Inspire ID: 736052

Status: VALIDATED

Authors: - Maksim Davydov - Hannes Jung

References: - Eur.Phys.J.C 51 (2007) 271 - DOI: 10.1140/epjc/s10052-007-0296-5 - arXiv: hep-ex/0701023 - DESY-06-240

Beams: e+ p+, p+ e+

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

Run details: none listed

Inclusive $D^$ production is measured in deep-inelastic $ep$ scattering at HERA with the H1 detector. In addition, the production of dijets in events with a $D^$ meson is investigated. The analysis covers values of photon virtuality $2< Q^2 \leq 100$ GeV$^2$ and of inelasticity $0.05 \leq y \leq 0.7$. Differential cross sections are measured as a function of $Q^2$ and $x$ and of various $D^*$ meson and jet observables.

Source code:H1_2007_I736052.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/DISFinalState.hh"

include "Rivet/Projections/DISKinematics.hh"

include "Rivet/Projections/DISLepton.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1) class H1_2007_I736052 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(H1_2007_I736052);


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

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

  declare(DISLepton(), "Lepton");
  declare(DISKinematics(), "Kinematics");

  const FinalState fs;
  declare(fs, "FS");
  const UnstableParticles ufs;
  declare(ufs, "UFS");

  double jet_radius = 1.0;
  const DISFinalState DISfs(DISFrame::BREIT);
  declare(FastJets(DISfs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::Et_scheme,
                   jet_radius),
          "jets_fs");

  Histo1DPtr tmp;
  book(_h["411"], 4, 1, 1);
  book(_h["511"], 5, 1, 1);
  book(_h["611"], 6, 1, 1);
  book(_h["711"], 7, 1, 1);
  book(_h["811"], 8, 1, 1);
  book(_h["911"], 9, 1, 1);
  book(_h_binned["Q2xbj"], {2., 4.22, 10., 17.8, 31.6, 100.},
       {"d10-x01-y01", "d11-x01-y01", "d12-x01-y01", "d13-x01-y01", "d14-x01-y01"});
  book(_h["1511"], 15, 1, 1);
  book(_h["1611"], 16, 1, 1);
  book(_h["1711"], 17, 1, 1);
  book(_h["1811"], 18, 1, 1);
  book(_h["1911"], 19, 1, 1);
  book(_h["2011"], 20, 1, 1);
  book(_h["2111"], 21, 1, 1);
  book(_h["2211"], 22, 1, 1);

  book(_h_binned["Q2phi"], {2., 10., 100.}, {"d23-x01-y01", "d24-x01-y01"});
  book(_h["2511"], 25, 1, 1);
  book(_h["2611"], 26, 1, 1);
  book(_h["2711"], 27, 1, 1);
  book(_h["2811"], 28, 1, 1);
  book(_h_binned["Q2xgam"], {2., 5., 10., 100.}, {"d29-x01-y01", "d29-x01-y02", "d29-x01-y03"});
  book(_h["3011"], 30, 1, 1);
  book(_h_binned["Q2xglue"], {2., 5., 10., 100.}, {"d31-x01-y01", "d31-x01-y02", "d31-x01-y03"});
}


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

  const FinalState& fs = apply<FinalState>(event, "FS");
  const size_t numParticles = fs.particles().size();

  Jets jets_fs = apply<JetFinder>(event, "jets_fs").jetsByPt(); // Jets with cut on eta
  double jet_radius = 1.0;

  const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");

  if (numParticles < 2) {
    MSG_DEBUG("Failed leptonic event cut");
    vetoEvent;
  }

  Particles Dstar;
  for (const Particle& p : select(ufs.particles(),
                                  Cuts::pT > 1.5 * GeV and Cuts::pT < 15 * GeV and Cuts::abseta < 1.5
                                      and Cuts::abspid == 413)) {
    Dstar.push_back(p);
  }
  if (Dstar.size() == 0) { // Cut on Dstar
    MSG_DEBUG("Failed Dstar cut");
    vetoEvent;
  }

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


  double Q2 = dk.Q2();
  double y = dk.y();

  if (y < 0.05 or y > 0.7 or Q2 < 2 or Q2 > 100) { // Cut on event kinematics
    MSG_DEBUG("Failed kinematics cut");
    vetoEvent;
  }

  // Extract the particles other than the lepton
  Particles particles;
  particles.reserve(fs.particles().size());
  ConstGenParticlePtr dislepGP = dl.out().genParticle();
  for (const Particle& p : fs.particles()) {
    ConstGenParticlePtr loopGP = p.genParticle();
    if (loopGP == dislepGP) continue;
    particles.push_back(p);
  }


  const LorentzTransform hcmboost = dk.boostHCM();        // Hadron cm system
  const LorentzTransform breitboost = dk.boostBreit();    //Breit system
  const LorentzTransform labboost = breitboost.inverse(); //Labsystem from Breit

  double xbj = dk.x();
  double W2 = dk.W2();
  double pT_cm(0);

  _h["411"]->fill(Q2);
  _h["511"]->fill(xbj);
  _h["611"]->fill(std::sqrt(W2));

  _h_binned["Q2xbj"]->fill(Q2, xbj);

  for (const Particle& p : Dstar) {
    _h["711"]->fill(p.pT());
    _h["811"]->fill(p.eta());
    _h["911"]->fill((p.E() - p.pz()) / (2 * y * dk.beamLepton().E()));

    const FourMomentum hcmMom = hcmboost.transform(p.momentum());
    if (pT_cm < hcmMom.pT()) pT_cm = hcmMom.pT();

    if (hcmMom.pT() > 2) { //Cut for 1711 and 1811 histograms
      _h["1711"]->fill(p.pT());
      _h["1811"]->fill(p.eta());
    }
  }

  if (pT_cm > 2) {
    _h["1511"]->fill(Q2);
    _h["1611"]->fill(xbj);
  }

  bool two_jets_with_cut(false);

  Jets jet_cut;

  for (const Jet& j : jets_fs) {
    double etalab = labboost.transform(j.momentum()).eta();
    if (j.momentum().Et() > 3 and etalab > -1 and etalab < 2.5) { // Cut on jet energy in Breit sys.
      jet_cut.push_back(j);
    }
  }

  if (jet_cut.size() >= 2 && jet_cut[0].momentum().Et() > 4) two_jets_with_cut = true;

  double delta_phi;
  FourMomentum jet1;
  FourMomentum jet2;

  Jets jet_DJ, jet_OJ;
  bool found_DJ(false), found_OJ(false);
  //cout << " check D and other jets " << endl;
  if (two_jets_with_cut) {
    jet1 = jet_cut[0].momentum(); // momentum of jet #1 in Breit sys.
    jet2 = jet_cut[1].momentum(); // momentum of jet #2 in Breit sys.

    delta_phi = deltaPhi(jet1, jet2);
    _h["1911"]->fill(Q2);
    _h["2011"]->fill(xbj);
    _h["2111"]->fill(jet_cut[0].momentum().Et() / GeV);
    _h["2211"]->fill(FourMomentum(jet1 + jet2).mass() / GeV); // Jets invariant mass in Breit sys.
    _h_binned["Q2phi"]->fill(Q2, delta_phi);
    for (const Jet& jet : jet_cut) {
      for (const Particle& p : Dstar) {
        if (deltaR(breitboost.transform(p.momentum()), jet.momentum()) < jet_radius) {
          jet_DJ.push_back(jet);
          found_DJ = true;
        }
        else {
          jet_OJ.push_back(jet);
          found_OJ = true;
        }
      }

      if (found_DJ && found_OJ) break;
    }
    if (jet_DJ.size() > 0 && jet_OJ.size() > 0) {
      double eta_DJ_breit = jet_DJ[0].momentum().eta();
      double eta_OJ_breit = jet_OJ[0].momentum().eta();

      _h["2511"]->fill(eta_DJ_breit);
      _h["2611"]->fill(eta_OJ_breit);
      _h["2711"]->fill(abs(eta_DJ_breit - eta_OJ_breit));

      double x_gamma, x_gluon;
      double E_star_p_z_star(0); // E() - pz() sum for x_gamma in \gamma p cm
      double E_star_p_z_had(0);

      // for jets in had cms: boost first back to lab and then to hcm
      Jet jet_DJ_hcm, jet_OJ_hcm;
      jet_DJ_hcm = hcmboost.transform(labboost(jet_DJ[0].momentum()));
      jet_OJ_hcm = hcmboost.transform(labboost(jet_OJ[0].momentum()));

      // Need to change sign: by default hcmboost has gamma* in +z dir,
      // and p in -z dir, but here we need: gamma* -in -z and proton in +z.
      // so - -> +: watch out for E-pz -> E+pz  and exp(eta_OJ_hcm) -> exp(-eta_OJ_hcm)
      // this was already noted in  H1_2007_I746380

      double Et_DJ_hcm = jet_DJ_hcm.Et();
      double eta_DJ_hcm = jet_DJ_hcm.eta();
      double Et_OJ_hcm = jet_OJ_hcm.Et();
      double eta_OJ_hcm = jet_OJ_hcm.eta();

      //observed fraction of the photon momentum carried by the parton involved in the hard subprocess

      E_star_p_z_star = jet_DJ_hcm.momentum().E() + jet_DJ_hcm.momentum().pz() + jet_OJ_hcm.momentum().E()
          + jet_OJ_hcm.momentum().pz();

      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
        const Particle& p = particles[ip1];
        E_star_p_z_had += (hcmboost.transform(p.momentum())).E()
            + (hcmboost.transform(p.momentum())).pz();
      }
      x_gamma = E_star_p_z_star / E_star_p_z_had;
      x_gluon = (Et_OJ_hcm * exp(-eta_OJ_hcm) + Et_DJ_hcm * exp(-eta_DJ_hcm))
          / (2. * hcmboost.transform(dk.beamHadron().momentum()).E());
      //observed fraction of the proton momentum carried by the gluon

      _h["2811"]->fill(x_gamma);
      _h_binned["Q2xgam"]->fill(Q2, x_gamma);

      _h["3011"]->fill(log10(x_gluon));
      _h_binned["Q2xglue"]->fill(Q2, log10(x_gluon));
    }
  }
}


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

  const double norm = crossSection() / nanobarn / sumW();
  scale(_h, norm);
  scale(_h_binned, norm);
  divByGroupWidth(_h_binned);

  // new scaling needed, since x bins are in log10(x)
  for (auto& b : _h["3011"]->bins()) {
    const double scale_new = b.xWidth() / pow(10, b.xWidth());
    // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
    const double factor = pow(10, b.xMid()) / 2.3026;
    b.scaleW(scale_new / factor);
  }

  // new scaling needed, since x bins are in log10(x)
  for (auto& histo : _h_binned["Q2xglue"]->bins()) {
    for (auto& b : histo->bins()) {
      const double scale_new = b.xWidth() / pow(10, b.xWidth());
      // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
      const double factor = pow(10, b.xMid()) / 2.3026;
      b.scaleW(scale_new / factor);
    }
  }
}

///@}


/// @name Histograms
///@{
map<string, Histo1DPtr> _h;
map<string, Histo1DGroupPtr> _h_binned;
///@}

};

RIVET_DECLARE_PLUGIN(H1_2007_I736052);

} ```