rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_1994_S2919893

H1 energy flow and charged particle spectra in DIS
Experiment: H1 (HERA)
Inspire ID: 372350
Status: VALIDATED
Authors:
  • Peter Richardson
References: Beams: p+ e-, p+ e+
Beam energies: (820.0, 26.7) GeV
Run details:
  • $e^- p$ / $e^+ p$ deep inelastic scattering, 820 GeV protons colliding with 26.7 GeV electrons

Global properties of the hadronic final state in deep inelastic scattering events at HERA are investigated. The data are corrected for detector effects. Energy flows in both the laboratory frame and the hadronic centre of mass system, and energy-energy correlations in the laboratory frame are presented. Historically, the Ariadne colour dipole model provided the only satisfactory description of this data, hence making it a useful 'target' analysis for MC shower models.

Source code: H1_1994_S2919893.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Math/Constants.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief H1 energy flow and charged particle spectra
 11  ///
 12  /// @author Peter Richardson
 13  /// Based on the equivalent HZTool analysis
 14  class H1_1994_S2919893 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1994_S2919893);
 19
 20
 21    /// @name Analysis methods
 22    //@{
 23
 24    /// Initialise projections and histograms
 25    void init() {
 26      // Projections
 27      declare(DISLepton(), "Lepton");
 28      declare(DISKinematics(), "Kinematics");
 29      declare(FinalState(), "FS");
 30
 31      // Histos
 32      book(_histEnergyFlowLowX ,1, 1, 1);
 33      book(_histEnergyFlowHighX ,1, 1, 2);
 34
 35      book(_histEECLowX ,2, 1, 1);
 36      book(_histEECHighX ,2, 1, 2);
 37
 38      book(_histSpectraW77 ,3, 1, 1);
 39      book(_histSpectraW122 ,3, 1, 2);
 40      book(_histSpectraW169 ,3, 1, 3);
 41      book(_histSpectraW117 ,3, 1, 4);
 42
 43      book(_histPT2 ,4, 1, 1);
 44
 45      book(_w77 .first, "TMP/w77_1");
 46      book(_w122.first, "TMP/w122_1");
 47      book(_w169.first, "TMP/w169_1");
 48      book(_w117.first, "TMP/w117_1");
 49      book(_wEnergy.first, "TMP/wEnergy_1");
 50
 51      book(_w77 .second, "TMP/w77_2");
 52      book(_w122.second, "TMP/w122_2");
 53      book(_w169.second, "TMP/w169_2");
 54      book(_w117.second, "TMP/w117_2");
 55      book(_wEnergy.second, "TMP/wEnergy_2");
 56    }
 57
 58
 59    /// Analyse each event
 60    void analyze(const Event& event) {
 61
 62      // Get the DIS kinematics
 63      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 64      if ( dk.failed() ) vetoEvent;
 65      const double x  = dk.x();
 66      const double w2 = dk.W2();
 67      const double w = sqrt(w2);
 68
 69      // Momentum of the scattered lepton
 70      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 71      if ( dl.failed() ) return;
 72      const FourMomentum leptonMom = dl.out();
 73      const double ptel = leptonMom.pT();
 74      const double enel = leptonMom.E();
 75      const double thel = leptonMom.angle(dk.beamHadron().mom())/degree;
 76
 77      // Extract the particles other than the lepton
 78      const FinalState& fs = apply<FinalState>(event, "FS");
 79      Particles particles;
 80      particles.reserve(fs.particles().size());
 81      ConstGenParticlePtr dislepGP = dl.out().genParticle();
 82      for(const Particle& p: fs.particles()) {
 83        ConstGenParticlePtr loopGP = p.genParticle();
 84        if (loopGP == dislepGP) continue;
 85        particles.push_back(p);
 86      }
 87
 88      // Cut on the forward energy
 89      double efwd = 0.0;
 90      for (const Particle& p : particles) {
 91        const double th = p.angle(dk.beamHadron())/degree;
 92        if (inRange(th, 4.4, 15)) efwd += p.E();
 93      }
 94
 95      // Apply the cuts
 96      // Lepton energy and angle, w2 and forward energy
 97      MSG_DEBUG("enel/GeV = " << enel/GeV << ", thel = " << thel
 98                << ", w2 = " << w2 << ", efwd/GeV = " << efwd/GeV);
 99      bool cut = enel/GeV > 14. && thel > 157. && thel < 172.5 && w2 >= 3000. && efwd/GeV > 0.5;
100      if (!cut) vetoEvent;
101
102      // Weight of the event
103      (x < 1e-3 ? _wEnergy.first : _wEnergy.second)->fill();
104
105      // Boost to hadronic CM
106      const LorentzTransform hcmboost = dk.boostHCM();
107      // Loop over the particles
108      long ncharged(0);
109      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
110        const Particle& p = particles[ip1];
111
112        const double th = p.angle(dk.beamHadron().momentum()) / degree;
113        // Boost momentum to lab
114        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
115        // Angular cut
116        if (th <= 4.4) continue;
117
118        // Energy flow histogram
119        const double et = fabs(hcmMom.Et());
120        const double eta = hcmMom.eta();
121        (x < 1e-3 ? _histEnergyFlowLowX : _histEnergyFlowHighX)->fill(eta, et);
122        if (PID::charge3(p.pid()) != 0) {
123          /// @todo Use units in w comparisons... what are the units?
124          if (w > 50. && w <= 200.) {
125            double xf= 2 * hcmMom.z() / w;
126            double pt2 = hcmMom.pT2();
127            if (w > 50. && w <= 100.) {
128              _histSpectraW77 ->fill(xf);
129            } else if (w > 100. && w <= 150.) {
130              _histSpectraW122->fill(xf);
131            } else if (w > 150. && w <= 200.) {
132              _histSpectraW169->fill(xf);
133            }
134            _histSpectraW117->fill(xf);
135            /// @todo Is this profile meant to be filled with 2 weight factors?
136            _histPT2->fill(xf, pt2/GeV2);
137            ++ncharged;
138          }
139        }
140
141
142        // Energy-energy correlation
143        if (th <= 8.) continue;
144        double phi1 = p.phi(ZERO_2PI);
145        double eta1 = p.eta();
146        double et1 = fabs(p.momentum().Et());
147        for (size_t ip2 = ip1+1; ip2 < particles.size(); ++ip2) {
148          const Particle& p2 = particles[ip2];
149
150          //double th2 = beamAngle(p2.momentum(), order);
151          double th2 = p2.angle(dk.beamHadron().momentum()) / degree;
152          if (th2 <= 8.) continue;
153          double phi2 = p2.phi(ZERO_2PI);
154
155          /// @todo Use angle function
156          double deltaphi = phi1 - phi2;
157          if (fabs(deltaphi) > PI) deltaphi = fabs(fabs(deltaphi) - TWOPI);
158          double eta2 = p2.eta();
159          double omega = sqrt(sqr(eta1-eta2) + sqr(deltaphi));
160          double et2 = fabs(p2.momentum().Et());
161          double wt = et1*et2 / sqr(ptel);
162          (x < 1e-3 ? _histEECLowX : _histEECHighX)->fill(omega, wt);
163        }
164      }
165
166      // Factors for normalization
167      if (w > 50. && w <= 200.) {
168        if (w <= 100.) {
169          _w77.first ->fill(ncharged);
170          _w77.second->fill();
171        } else if (w <= 150.) {
172          _w122.first ->fill(ncharged);
173          _w122.second->fill();
174        } else {
175          _w169.first ->fill(ncharged);
176          _w169.second->fill();
177        }
178        _w117.first ->fill(ncharged);
179        _w117.second->fill();
180      }
181    }
182
183
184    // Normalize inclusive single particle distributions to the average number of charged particles per event.
185    void finalize() {
186      normalize(_histSpectraW77,  *_w77.first/ *_w77.second);
187      normalize(_histSpectraW122,  *_w122.first/ *_w122.second);
188      normalize(_histSpectraW169,  *_w169.first/ *_w169.second);
189      normalize(_histSpectraW117,  *_w117.first/ *_w117.second);
190
191      scale(_histEnergyFlowLowX , 1./ *_wEnergy.first );
192      scale(_histEnergyFlowHighX, 1./ *_wEnergy.second);
193
194      scale(_histEECLowX , 1./ *_wEnergy.first );
195      scale(_histEECHighX, 1./ *_wEnergy.second);
196    }
197
198    //@}
199
200
201  private:
202
203    /// Polar angle with right direction of the beam
204    inline double beamAngle(const FourVector& v, bool order) {
205      double thel = v.polarAngle()/degree;
206      if (thel < 0) thel += 180.;
207      if (!order) thel = 180 - thel;
208      return thel;
209    }
210
211    /// @name Histograms
212    /// @{
213    Histo1DPtr _histEnergyFlowLowX, _histEnergyFlowHighX;
214    Histo1DPtr _histEECLowX, _histEECHighX;
215    Histo1DPtr _histSpectraW77, _histSpectraW122, _histSpectraW169, _histSpectraW117;
216    Profile1DPtr _histPT2;
217    /// @}
218
219    /// @name Storage of weights to calculate averages for normalisation
220    /// @{
221    pair<CounterPtr,CounterPtr> _w77, _w122, _w169, _w117, _wEnergy;
222    /// @}
223
224  };
225
226
227
228  RIVET_DECLARE_ALIASED_PLUGIN(H1_1994_S2919893, H1_1994_I372350);
229
230}