rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_2000_S4129130

H1 energy flow in DIS
Experiment: H1 (HERA)
Inspire ID: 503947
Status: VALIDATED
Authors:
  • Peter Richardson
References: Beams: p+ e+
Beam energies: (820.0, 27.5) GeV
Run details:
  • $e^+ p$ deep inelastic scattering with $p$ at 820 GeV, $e^+$ at 27.5 GeV \to $\sqrt{s} = 300 \text{GeV}$

Measurements of transverse energy flow for neutral current deep- inelastic scattering events produced in positron-proton collisions at HERA. The kinematic range covers squared momentum transfers $Q^2$ from 3.2 to 2200 GeV$^2$; the Bjorken scaling variable $x$ from $8 \times 10^{-5}$ to 0.11 and the hadronic mass $W$ from 66 to 233 GeV. The transverse energy flow is measured in the hadronic centre of mass frame and is studied as a function of $Q^2$, $x$, $W$ and pseudorapidity. The behaviour of the mean transverse energy in the central pseudorapidity region and an interval corresponding to the photon fragmentation region are analysed as a function of $Q^2$ and $W$. This analysis is useful for exploring the effect of photon PDFs and for tuning models of parton evolution and treatment of fragmentation and the proton remnant in DIS.

Source code: H1_2000_S4129130.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  ///
 14  /// Based on the HZTOOL analysis HZ99091
 15  class H1_2000_S4129130 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2000_S4129130);
 20
 21
 22    /// @name Analysis methods
 23    //@{
 24
 25    /// Initialise projections and histograms
 26    void init() {
 27      // Projections
 28      declare(DISLepton(), "Lepton");
 29      declare(DISKinematics(), "Kinematics");
 30      declare(FinalState(), "FS");
 31
 32      // Histograms and weight vectors for low Q^2 a
 33      _histETLowQa.resize(17);
 34      for (size_t ix = 0; ix < 17; ++ix) {
 35        book(_histETLowQa[ix], ix+1, 1, 1);
 36        book(_weightETLowQa[ix], "TMP/ETLowQa" + to_string(ix));
 37      }
 38
 39      // Histograms and weight vectors for high Q^2 a
 40      _histETHighQa.resize(7);
 41      for (size_t ix = 0; ix < 7; ++ix) {
 42        book(_histETHighQa[ix], ix+18, 1, 1);
 43        book(_weightETHighQa[ix], "TMP/ETHighQa" + to_string(ix));
 44      }
 45
 46      // Histograms and weight vectors for low Q^2 b
 47      _histETLowQb.resize(5);
 48      for (size_t ix = 0; ix < 5; ++ix) {
 49        book(_histETLowQb[ix], ix+25, 1, 1);
 50        book(_weightETLowQb[ix], "TMP/ETLowQb" + to_string(ix));
 51      }
 52
 53      // Histograms and weight vectors for high Q^2 b
 54      _histETHighQb.resize(5);
 55      for (size_t ix = 0; ix < 3; ++ix) {
 56        book(_histETHighQb[ix], 30+ix, 1, 1);
 57        book(_weightETHighQb[ix], "TMP/ETHighQb" + to_string(ix));
 58      }
 59
 60      // Histograms for the averages
 61      book(_histAverETCentral ,33,  1, 1);
 62      book(_histAverETFrag ,34,  1, 1);
 63    }
 64
 65
 66    /// Analyze each event
 67    void analyze(const Event& event) {
 68
 69      // DIS kinematics
 70      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 71      if ( dk.failed() ) vetoEvent;
 72      double q2  = dk.Q2();
 73      double x   = dk.x();
 74      double y   = dk.y();
 75      double w2  = dk.W2();
 76
 77      // Kinematics of the scattered lepton
 78      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 79      if ( dl.failed() ) vetoEvent;
 80      const FourMomentum leptonMom = dl.out();
 81      const double enel = leptonMom.E();
 82      const double thel = 180 - leptonMom.angle(dl.in().mom())/degree;
 83
 84      // Extract the particles other than the lepton
 85      const FinalState& fs = apply<FinalState>(event, "FS");
 86      Particles particles; particles.reserve(fs.size());
 87      ConstGenParticlePtr dislepGP = dl.out().genParticle(); ///< @todo Is the GenParticle stuff necessary? (Not included in Particle::==?)
 88      for(const Particle& p: fs.particles()) {
 89        ConstGenParticlePtr loopGP = p.genParticle();
 90        if (loopGP == dislepGP) continue;
 91        particles.push_back(p);
 92      }
 93
 94      // Cut on the forward energy
 95      double efwd = 0.;
 96      for (const Particle& p : particles) {
 97        const double th = 180 - p.angle(dl.in())/degree;
 98        if (inRange(th, 4.4, 15.0)) efwd += p.E();
 99      }
100      // There are four possible selections for events
101      bool evcut[4];
102      // Low  Q2 selection a
103      evcut[0] = enel/GeV > 12. && w2 >= 4400.*GeV2 && efwd/GeV > 0.5 && inRange(thel,157.,176.);
104      // Low  Q2 selection b
105      evcut[1] = enel/GeV > 12. && inRange(y,0.3,0.5);
106      // High Q2 selection a
107      evcut[2] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && w2 >= 4400.*GeV2 && efwd > 0.5;
108      // High Q2 selection b
109      evcut[3] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && inRange(w2,27110.*GeV2,45182.*GeV2);
110
111      // Veto if fails all cuts
112      /// @todo Can we use all()?
113      if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3]) ) vetoEvent;
114
115      // Find the bins
116      int bin[4] = {-1,-1,-1,-1};
117      // For the low Q2 selection a)
118      if (q2 > 2.5*GeV2 && q2 <= 5.*GeV2) {
119        if (x > 0.00005 && x <= 0.0001 ) bin[0] = 0;
120        if (x > 0.0001  && x <= 0.0002 ) bin[0] = 1;
121        if (x > 0.0002  && x <= 0.00035) bin[0] = 2;
122        if (x > 0.00035 && x <= 0.0010 ) bin[0] = 3;
123      }
124      else if (q2 > 5.*GeV2 && q2 <= 10.*GeV2) {
125        if (x > 0.0001  && x <= 0.0002 ) bin[0] = 4;
126        if (x > 0.0002  && x <= 0.00035) bin[0] = 5;
127        if (x > 0.00035 && x <= 0.0007 ) bin[0] = 6;
128        if (x > 0.0007  && x <= 0.0020 ) bin[0] = 7;
129      }
130      else if (q2 > 10.*GeV2 && q2 <= 20.*GeV2) {
131        if (x > 0.0002 && x <= 0.0005) bin[0] = 8;
132        if (x > 0.0005 && x <= 0.0008) bin[0] = 9;
133        if (x > 0.0008 && x <= 0.0015) bin[0] = 10;
134        if (x > 0.0015 && x <= 0.040 ) bin[0] = 11;
135      }
136      else if (q2 > 20.*GeV2 && q2 <= 50.*GeV2) {
137        if (x > 0.0005 && x <= 0.0014) bin[0] = 12;
138        if (x > 0.0014 && x <= 0.0030) bin[0] = 13;
139        if (x > 0.0030 && x <= 0.0100) bin[0] = 14;
140      }
141      else if (q2 > 50.*GeV2 && q2 <= 100.*GeV2) {
142        if (x >0.0008 && x <= 0.0030) bin[0] = 15;
143        if (x >0.0030 && x <= 0.0200) bin[0] = 16;
144      }
145      // check in one of the bins
146      evcut[0] &= bin[0] >= 0;
147      // For the low Q2 selection b)
148      if (q2 > 2.5*GeV2 && q2 <= 5.  *GeV2) bin[1] = 0;
149      if (q2 > 5. *GeV2 && q2 <= 10. *GeV2) bin[1] = 1;
150      if (q2 > 10.*GeV2 && q2 <= 20. *GeV2) bin[1] = 2;
151      if (q2 > 20.*GeV2 && q2 <= 50. *GeV2) bin[1] = 3;
152      if (q2 > 50.*GeV2 && q2 <= 100.*GeV2) bin[1] = 4;
153      // check in one of the bins
154      evcut[1] &= bin[1] >= 0;
155      // for the high Q2 selection a)
156      if (q2 > 100.*GeV2 && q2 <= 400.*GeV2) {
157        if (x > 0.00251 && x <= 0.00631) bin[2] = 0;
158        if (x > 0.00631 && x <= 0.0158 ) bin[2] = 1;
159        if (x > 0.0158  && x <= 0.0398 ) bin[2] = 2;
160      }
161      else if (q2 > 400.*GeV2 && q2 <= 1100.*GeV2) {
162        if (x > 0.00631 && x <= 0.0158 ) bin[2] = 3;
163        if (x > 0.0158  && x <= 0.0398 ) bin[2] = 4;
164        if (x > 0.0398  && x <= 1.     ) bin[2] = 5;
165      }
166      else if (q2 > 1100.*GeV2 && q2 <= 100000.*GeV2) {
167        if (x > 0. && x <= 1.) bin[2] = 6;
168      }
169      // check in one of the bins
170      evcut[2] &= bin[2] >= 0;
171      // for the high Q2 selection b)
172      if      (q2 > 100.*GeV2 && q2 <= 220.*GeV2) bin[3] = 0;
173      else if (q2 > 220.*GeV2 && q2 <= 400.*GeV2) bin[3] = 1;
174      else if (q2 > 400.              ) bin[3] = 2;
175      // check in one of*GeV the bins
176      evcut[3] &= bin[3] >= 0;
177
178      // Veto if fails all cuts after bin selection
179      /// @todo Can we use all()?
180      if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3])) vetoEvent;
181
182      // Increment the count for normalisation
183      if (evcut[0]) _weightETLowQa [bin[0]]->fill();
184      if (evcut[1]) _weightETLowQb [bin[1]]->fill();
185      if (evcut[2]) _weightETHighQa[bin[2]]->fill();
186      if (evcut[3]) _weightETHighQb[bin[3]]->fill();
187
188      // Boost to hadronic CoM
189      const LorentzTransform hcmboost = dk.boostHCM();
190
191      // Loop over the particles
192      double etcent = 0;
193      double etfrag = 0;
194      for (const Particle& p : particles) {
195        // Boost momentum to CMS
196        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
197        double et = fabs(hcmMom.Et());
198        double eta = hcmMom.eta();
199        // Averages in central and forward region
200        if (fabs(eta) < .5 ) etcent += et;
201        if (eta > 2 && eta <= 3.) etfrag += et;
202        // Histograms of Et flow
203        if (evcut[0]) _histETLowQa [bin[0]]->fill(eta, et);
204        if (evcut[1]) _histETLowQb [bin[1]]->fill(eta, et);
205        if (evcut[2]) _histETHighQa[bin[2]]->fill(eta, et);
206        if (evcut[3]) _histETHighQb[bin[3]]->fill(eta, et);
207      }
208      // Fill histograms for the average quantities
209      if (evcut[1] || evcut[3]) {
210        _histAverETCentral->fill(q2, etcent);
211        _histAverETFrag   ->fill(q2, etfrag);
212      }
213    }
214
215
216    // Finalize
217    void finalize() {
218      // Normalization of the Et distributions
219      /// @todo Simplify by using normalize() instead? Are all these being normalized to area=1?
220      for (size_t ix = 0; ix < 17; ++ix) if (_weightETLowQa[ix]->val()  != 0) scale(_histETLowQa[ix],  1/ *_weightETLowQa[ix]);
221      for (size_t ix = 0; ix <  7; ++ix) if (_weightETHighQa[ix]->val() != 0) scale(_histETHighQa[ix], 1/ *_weightETHighQa[ix]);
222      for (size_t ix = 0; ix <  5; ++ix) if (_weightETLowQb[ix]->val()  != 0) scale(_histETLowQb[ix],  1/ *_weightETLowQb[ix]);
223      for (size_t ix = 0; ix <  3; ++ix) if (_weightETHighQb[ix]->val() != 0) scale(_histETHighQb[ix], 1/ *_weightETHighQb[ix]);
224    }
225
226    //@}
227
228
229  private:
230
231    /// @name Histograms
232    /// @{
233    vector<Histo1DPtr> _histETLowQa;
234    vector<Histo1DPtr> _histETHighQa;
235    vector<Histo1DPtr> _histETLowQb;
236    vector<Histo1DPtr> _histETHighQb;
237    Profile1DPtr _histAverETCentral;
238    Profile1DPtr _histAverETFrag;
239    /// @}
240
241    /// @name Storage of weights for normalisation
242    /// @{
243    array<CounterPtr,17> _weightETLowQa;
244    array<CounterPtr, 7> _weightETHighQa;
245    array<CounterPtr, 5> _weightETLowQb;
246    array<CounterPtr, 3> _weightETHighQb;
247    /// @}
248
249  };
250
251
252
253  RIVET_DECLARE_ALIASED_PLUGIN(H1_2000_S4129130, H1_2000_I503947);
254
255}