rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_1999_I504022

Forward pi0 meson production at HERA
Experiment: H1 (HERA)
Inspire ID: 504022
Status: VALIDATED
Authors:
  • Keila Moral Figueroa
  • Hannes Jung
References:
  • Phys.Lett.B462:440-452,1999
  • DOI: 10.1016/S0370-2693(99)00906-5
  • arXiv: hep-ex/9907030
Beams: p+ e-, e- p+
Beam energies: (820.0, 27.5); (27.5, 820.0) GeV
Run details:
  • $\pi^0$ production in the forward region detected by H1 in DIS ep scattering events at HERA. Cuts:0.1 $< y <$ 0.6, 2 $< Q^2 <$ 70 $GeV^2$, $\pi^0$ transverse momentum in the hadronic CMS. $p^*_T$ $ > $ 2.5 GeV, 5$^o< \theta < 25^o$, where $\theta$ is the polar angle of the produced $\pi^0$,$x_{\pi^0} > 0.001$ where $ x_{\pi^0} = E_{\pi^0}/ E_{proton}$ and $E_{proton} = 820GeV$ is the proton beam energy. RAPGAP 26.7 configuration parameters: 10^6 events.

High transverse momentum $\pi^0$ mesons have been measured with the H1 detector at HERA in deep-inelastic ep scattering events at low Bjorken-x, down to $x \approx 4 \cdot 10^{-5}$ .The measurement is performed in a region of small angles with respect to the proton remnant in the laboratory frame of reference, namely the forward region, and corresponds to central rapidity in the centre of mass system of the virtual photon and proton. This region is expected to be particularly sensitive to QCD effects in hadronic final states. Differential cross-sections for inclusive $\pi^0$ meson production are presented as a function of Bjorken-x and the four-momentum transfer $Q^2$, and as a function of transverse momentum and pseudorapidity.

Source code: H1_1999_I504022.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Forward pi0 meson production at HERA (H1)
 11  class H1_1999_I504022 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1999_I504022);
 16
 17
 18    /// @name Analysis methods
 19    ///@{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25
 26      // The basic final-state projection:
 27
 28      declare(FinalState(Cuts::abseta < 7 ), "FS");
 29      declare(DISKinematics(), "Kinematics");
 30      declare(UnstableParticles(), "UFS");
 31
 32
 33      // Book histograms
 34      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
 35
 36      book(_h["p_T>2.5&x"], 1, 1, 1);
 37      book(_h["Q:2.0-4.5-eta"], 2, 1, 1);
 38      book(_h["Q:2.0-4.5-p_T"], 3, 1, 1);
 39      book(_h["Q:4.5-15.0-x"], 4, 1, 1);
 40      book(_h["Q:4.5-15.0-eta"], 5, 1, 1);
 41      book(_h["Q:4.5-15.0-p_T"], 6, 1, 1);
 42      book(_h["Q:15.0-70.0-x"], 7, 1, 1);
 43      book(_h["Q:15.0-70.0-eta"], 8, 1, 1);
 44      book(_h["Q:15.0-70.0-p_T"], 9, 1, 1);
 45      book(_h["p_T>2.5&Q"], 10, 1, 1);
 46      book(_h["p_T>3.5&x"], 11, 1, 1);
 47      book(_h["p_T>3.5&Q"], 12, 1, 1);
 48    }
 49
 50
 51    /// Perform the per-event analysis
 52    void analyze(const Event& event) {
 53
 54      /// @todo Do the event by event analysis here
 55
 56        const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 57
 58        // Get the DIS kinematics
 59        double xbj  = dk.x();
 60        double ybj = dk.y();
 61        double Q2 = dk.Q2()/GeV2;
 62
 63        // Q2 and inelasticity cuts
 64        //cout << " after xbj " << xbj << endl;
 65        if (!inRange(ybj, 0.1, 0.6)) vetoEvent;
 66        //cout << " after ybj " << ybj << endl;
 67        if (!inRange(Q2, 2.0*GeV2, 70.0*GeV2)) vetoEvent;
 68        //cout << " after Q2 " << Q2 << endl;
 69
 70      const FinalState& fs = apply<FinalState>(event, "FS");
 71      const size_t numParticles = fs.particles().size();
 72      //cout << " Num all     final state particles " << numParticles << endl;
 73      //proton beam energy: 820 GeV
 74      double e_proton = dk.beamHadron().E()/GeV;
 75      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 76      if (numParticles < 2) {
 77        MSG_DEBUG("Failed leptonic event cut");
 78        vetoEvent;
 79      }
 80
 81      // Extracting the pi0
 82      const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");
 83      //Get the hadronic CMS kinematics
 84      const LorentzTransform hcmboost = dk.boostHCM();
 85
 86      for (const Particle& p : ufs.particles(Cuts::pid==PID::PI0)) {
 87        //cout << " Pid = " << p.pid() << endl;
 88        //Get the LAB kinematics
 89        double theta = p.theta();
 90        double eta = p.pseudorapidity();
 91        // double pT = p.pT()/GeV;
 92        //Boost hcm
 93        const FourMomentum hcmMom = hcmboost.transform(p.mom());
 94        double pThcm = hcmMom.pT();
 95
 96        double e_pi0 = p.E()/GeV;
 97        double x_pi0_proton = e_pi0/e_proton;
 98
 99        // epi0/e_proton, theta and pThcm cuts
100        if(x_pi0_proton < 0.01) continue;
101        //cout << " theta " << theta << " in deg " << theta/degree << endl;
102        if (!inRange(theta/degree, 5, 25)) continue;
103        if(pThcm < 2.5*GeV) continue;
104
105        //Three cuts for Q2:
106        if (Q2 > 2.0*GeV2 && Q2 < 4.5*GeV2){
107          _h["Q:2.0-4.5-eta"]->fill(eta);
108          _h["Q:2.0-4.5-p_T"]->fill(pThcm);
109        }
110        if (Q2 > 4.5*GeV2 && Q2 < 15.0*GeV2){
111          _h["Q:4.5-15.0-x"]->fill(xbj);
112          _h["Q:4.5-15.0-eta"]->fill(eta);
113          _h["Q:4.5-15.0-p_T"]->fill(pThcm);
114        }
115        if (Q2 > 15.0*GeV2 && Q2 < 70.0*GeV2){
116          _h["Q:15.0-70.0-x"]->fill(xbj);
117          _h["Q:15.0-70.0-eta"]->fill(eta);
118          _h["Q:15.0-70.0-p_T"]->fill(pThcm);
119        }
120
121        //Two cuts for p_T:
122        if (pThcm > 2.5*GeV){
123          _h["p_T>2.5&x"]->fill(xbj);
124          _h["p_T>2.5&Q"]->fill(Q2);
125        }
126        if (pThcm > 3.5*GeV){
127          _h["p_T>3.5&x"]->fill(xbj);
128          _h["p_T>3.5&Q"]->fill(Q2);
129        }
130
131      }
132
133    }
134
135
136    /// Normalise histograms etc., after the run
137    void finalize() {
138      const double sn = crossSection()/nanobarn/sumW();
139      const double sp = crossSection()/picobarn/sumW();
140      scale(_h["p_T>2.5&x"], sn);
141      scale(_h["Q:2.0-4.5-eta"], sp);
142      scale(_h["Q:2.0-4.5-p_T"], sp);
143      scale(_h["Q:4.5-15.0-x"], sn);
144      scale(_h["Q:4.5-15.0-eta"], sp);
145      scale(_h["Q:4.5-15.0-p_T"], sp);
146      scale(_h["Q:15.0-70.0-x"], sn);
147      scale(_h["Q:15.0-70.0-eta"], sp);
148      scale(_h["Q:15.0-70.0-p_T"], sp);
149      scale(_h["p_T>2.5&Q"], sp);
150      scale(_h["p_T>3.5&x"], sn);
151      scale(_h["p_T>3.5&Q"], sp);
152    }
153
154    ///@}
155
156
157    /// @name Histograms
158    ///@{
159    map<string, Histo1DPtr> _h;
160    ///@}
161
162
163  };
164
165
166  RIVET_DECLARE_PLUGIN(H1_1999_I504022);
167
168}