1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/DISFinalState.hh"
  4#include "Rivet/Projections/CentralEtHCM.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// H1 energy flow in DIS
 10  ///
 11  /// @todo Make histograms match those in HepData and use autobooking
 12  ///
 13  /// @author Leif Lonnblad
 14  /// @author Andy Buckley
 15  class H1_1995_I396365 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1995_I396365);
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    void init() {
 26      // Projections
 27      declare(DISKinematics(), "Kinematics");
 28      const DISFinalState& fshcm = declare(DISFinalState(DISFrame::HCM), "FS");
 29      declare(CentralEtHCM(fshcm), "Y1HCM");
 30
 31      // Histograms
 32      /// @todo Convert to use autobooking and correspond to HepData data tables
 33
 34      _hEtFlow.resize(9);
 35      for (size_t i = 0; i < 9; ++i) {
 36        book(_sumw[i], "sumW_" + to_str(i));
 37        book(_hEtFlow[i] ,to_str(i), 24, -6, 6);
 38      }
 39      book(_tmphAvEt, "TMP/hAvEt", 9, 1.0, 10.0);
 40      book(_tmphAvX , "TMP/hAvX", 9, 1.0, 10.0);
 41      book(_tmphAvQ2, "TMP/hAvQ2", 9, 1.0, 10.0);
 42      book(_tmphN   , "TMP/hN", 9, 1.0, 10.0);
 43    }
 44
 45
 46    /// Calculate the bin number from the DISKinematics projection
 47    /// @todo Convert to use a HEPUtils Binning1D
 48    size_t _getbin(const DISKinematics& dk) {
 49      if (inRange(dk.Q2()/GeV2, 5.0, 10.0)) {
 50        if (inRange(dk.x(), 1e-4, 2e-4)) return 0;
 51        if (inRange(dk.x(), 2e-4, 5e-4) && dk.Q2() > 6.0*GeV2) return 1;
 52      } else if (inRange(dk.Q2()/GeV2, 10.0, 20.0)) {
 53        if (inRange(dk.x(), 2e-4, 5e-4)) return 2;
 54        if (inRange(dk.x(), 5e-4, 8e-4)) return 3;
 55        if (inRange(dk.x(), 8e-4, 1.5e-3)) return 4;
 56        if (inRange(dk.x(), 1.5e-3, 4e-3)) return 5;
 57      } else if (inRange(dk.Q2()/GeV2, 20.0, 50.0)) {
 58        if (inRange(dk.x(), 5e-4, 1.4e-3)) return 6;
 59        if (inRange(dk.x(), 1.4e-3, 3e-3)) return 7;
 60        if (inRange(dk.x(), 3e-3, 1e-2)) return 8;
 61      }
 62      return -1;
 63    }
 64
 65
 66    void analyze(const Event& event) {
 67      const FinalState& fs = apply<FinalState>(event, "FS");
 68      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 69      if ( dk.failed() ) vetoEvent;
 70      const CentralEtHCM& y1 = apply<CentralEtHCM>(event, "Y1HCM");
 71      if ( y1.failed() ) vetoEvent;
 72
 73      const int ibin = _getbin(dk);
 74      if (ibin < 0) vetoEvent;
 75
 76      _sumw[ibin]->fill();
 77
 78      for (size_t i = 0, N = fs.particles().size(); i < N; ++i) {
 79        const double rap = fs.particles()[i].rapidity();
 80        const double et = fs.particles()[i].Et();
 81        _hEtFlow[ibin]->fill(rap, et/GeV);
 82      }
 83
 84      _tmphAvEt->fill(ibin + 1.5, y1.sumEt()/GeV);
 85      _tmphAvX->fill(ibin + 1.5, dk.x());
 86      _tmphAvQ2->fill(ibin + 1.5, dk.Q2()/GeV2);
 87      _tmphN->fill(ibin + 1.5);
 88    }
 89
 90
 91    void finalize() {
 92      for (size_t ibin = 0; ibin < 9; ++ibin)
 93        scale(_hEtFlow[ibin], 0.5/ *_sumw[ibin]);
 94      /// @todo Improve this!
 95      Estimate1DPtr s21,s22,s23;
 96      book(s21, "21");
 97      divide(_tmphAvEt,_tmphN,s21);
 98      book(s22, "22");
 99      divide(_tmphAvX,_tmphN,s22);
100      book(s23, "23");
101      divide(_tmphAvQ2,_tmphN,s23);
102    }
103
104    /// @}
105
106
107  private:
108
109    /// Histograms for the \f$ E_T \f$ flow
110    vector<Histo1DPtr> _hEtFlow;
111
112    /// Temporary histograms for averages in different kinematical bins.
113    Histo1DPtr _tmphAvEt, _tmphAvX, _tmphAvQ2, _tmphN;
114
115    /// Weights counters for each kinematic bin
116    array<CounterPtr, 9> _sumw;
117
118  };
119
120
121
122  RIVET_DECLARE_ALIASED_PLUGIN(H1_1995_I396365, H1_1995_S3167097);
123
124}