1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6#include "Rivet/Projections/DISFinalState.hh"
  7#include "Rivet/Projections/Thrust.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Event-shape variables in deep-inelastic scattering at HERA
 13  class H1_2006_I699835 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2006_I699835);
 18
 19
 20    /// @name Analysis methods
 21    ///@{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Initialise and register projections
 27
 28      // The basic final-state projection:
 29      // all final-state particles within
 30      // the given eta acceptance
 31      const FinalState fs(Cuts::abseta < 4.9);
 32      declare(fs, "FS");
 33
 34      const DISFinalState DISfs(DISFrame::BREIT);
 35
 36      const FinalState DISfsCut(DISfs, Cuts::eta < 0);
 37
 38      declare(Thrust(DISfsCut), "ThrustCut");
 39
 40      declare(DISKinematics(), "Kinematics");
 41
 42      //Book histograms for different Q ranges:
 43      book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");
 44
 45      const vector<double> QEdges{14., 16., 20., 30., 50., 70., 100., 200.};
 46      book(_h_tauc, QEdges);
 47      book(_h_tau, QEdges);
 48      book(_h_B, QEdges);
 49      book(_h_rho, QEdges);
 50      for (size_t iQ=0; iQ < QEdges.size()-1; ++iQ) {
 51        book(_h_tauc->bin(iQ+1), 1+iQ,1,1);
 52        book(_h_tau->bin(iQ+1),  8+iQ,1,1);
 53        book(_h_B->bin(iQ+1),   15+iQ,1,1);
 54        book(_h_rho->bin(iQ+1), 22+iQ,1,1);
 55        book(_Nevt_after_cuts_Q[iQ], "TMP/Nevt_after_cuts_Q"+ to_string(iQ));
 56      }
 57
 58    }
 59
 60
 61    /// Perform the per-event analysis
 62    void analyze(const Event& event) {
 63
 64      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 65
 66      // The kinematic region covered by the analysis is defined by ranges of Q² and y
 67      if (dk.Q2() < 196 or dk.Q2() > 40000 or dk.y() < 0.1 or dk.y() > 0.7) {
 68      	vetoEvent;
 69      }
 70
 71      const double Q = sqrt(dk.Q2());
 72      _Nevt_after_cuts->fill();
 73      const size_t iQ = _h_tauc->binAt(Q).index();
 74      if (0 < iQ && iQ < 8)  _Nevt_after_cuts_Q[iQ-1]->fill();
 75
 76
 77      // Boost to hadronic Breit frame
 78      const LorentzTransform breitboost = dk.boostBreit();
 79
 80      const FinalState& fs = apply<FinalState>(event, "FS");
 81
 82      /*Calculate event shape variables:
 83      thrust_num is \sum |pz_h|
 84      thrust_den is \sum |p_h|
 85      b_num is \sum |pt_h|
 86      sumE is the sum of energies
 87      (All sums run through the particles in the current hemisphere)
 88      */
 89      double thrust_num, thrust_den, b_num, sumE;
 90      thrust_num = thrust_den = b_num = sumE = 0;
 91      Vector3 sumMom;
 92
 93      for (const Particle& p : fs.particles()) {
 94      	// Boost to Breit frame
 95        const FourMomentum breitMom = breitboost.transform(p.momentum());
 96      	if (breitMom.eta() < 0) {
 97          thrust_num += abs(breitMom.pz());
 98          thrust_den += breitMom.p();
 99          b_num += abs(breitMom.pt());
100          sumMom.operator+=(breitMom.p3());
101          sumE += breitMom.E();
102      	}
103      }
104
105      //The energy in the current hemisphere must exceed a certain value.
106      if (sumE <= Q/10.0)  vetoEvent;
107
108      /* Comment from A. Galvan:
109      Thrust here is with respect to the z axis (tau in the paper), while the
110      one from Rivet projection is with respect to the maximum thrust
111      axis (1 - tau_c in the paper)
112      */
113
114      double thrust_mine = 1.0 - ((double)thrust_num)/((double)thrust_den);
115      double b_mine = ((double)b_num)/(2.0*(double)thrust_den);
116      double rho_num = thrust_den*thrust_den - sumMom.dot(sumMom);
117      double rho_mine = ((double)rho_num)/(4.0*(double)thrust_den*(double)thrust_den);
118
119      const Thrust& thrCut = apply<Thrust>(event, "ThrustCut");
120
121      //Fill histograms:
122      _h_tauc->fill(Q, 1.0 - thrCut.thrust());
123      _h_tau->fill(Q, thrust_mine);
124      _h_B->fill(Q, b_mine);
125      _h_rho->fill(Q, rho_mine);
126
127      /* Comment from A. Galvan:
128       As for the C-parameter, my results did not fit the reference data at
129       all. The formula given in the paper is a bit ambiguous, because there is
130       a sum that runs through all pairs of particles h,h' and I was not sure
131       whether each pair should be counted twice (h,h' and h',h) or not, or if
132       the pair of a particle with itself (hh) should be considered.
133       That is why I tried all of the possibilities, but none of them worked.
134      */
135    }
136
137
138    /// Normalise histograms etc., after the run
139    void finalize() {
140
141      // need to multiply by bin widths
142      for (size_t iQ=0; iQ < _h_tauc->numBins(); ++iQ) {
143        const double Nev = dbl(*_Nevt_after_cuts_Q[iQ]);
144        if (Nev != 0) {
145          scale(_h_tauc->bin(iQ+1), 1./Nev);
146          for(auto & bin : _h_tauc->bin(iQ+1)->bins())
147            bin.scaleW(bin.xWidth());
148          scale(_h_tau->bin(iQ+1), 1./Nev);
149          for(auto & bin : _h_tau->bin(iQ+1)->bins())
150            bin.scaleW(bin.xWidth());
151          scale(_h_B->bin(iQ+1), 1./Nev);
152          for(auto & bin : _h_B->bin(iQ+1)->bins())
153            bin.scaleW(bin.xWidth());
154          scale(_h_rho->bin(iQ+1), 1./Nev);
155          for(auto & bin : _h_rho->bin(iQ+1)->bins())
156            bin.scaleW(bin.xWidth());
157        }
158      }
159    }
160
161    ///@}
162
163
164    /// @name Histograms
165    ///@{
166    Histo1DGroupPtr _h_tauc, _h_tau, _h_B, _h_rho;
167    CounterPtr _Nevt_after_cuts;
168    CounterPtr _Nevt_after_cuts_Q[7];
169    ///@}
170
171  };
172
173
174  RIVET_DECLARE_PLUGIN(H1_2006_I699835);
175
176}