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
 13  /// @brief Measurement of Event Shape Variables in Deep-Inelastic Scattering at HERA (H1)
 14  class H1_2006_I699835 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2006_I699835);
 19
 20
 21    /// @name Analysis methods
 22    ///@{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27      //cerr << endl << endl << "Initializing --------------------------------" << endl << endl;
 28      
 29      // Initialise and register projections
 30
 31      // The basic final-state projection:
 32      // all final-state particles within
 33      // the given eta acceptance
 34      const FinalState fs(Cuts::abseta < 4.9);
 35      declare(fs, "FS");
 36
 37      const DISFinalState DISfs(DISFinalState::BoostFrame::BREIT);
 38
 39      const FinalState DISfsCut(DISfs, Cuts::eta < 0);
 40
 41      declare(Thrust(DISfsCut), "ThrustCut");
 42      
 43      declare(DISKinematics(), "Kinematics");
 44      
 45      //Book histograms for different Q ranges:
 46      book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");
 47      
 48      Histo1DPtr dummy;
 49      for(int iQ = 0; iQ < 7; ++iQ) {
 50         // cout << " iQ " << iQ << " " << QEdges[iQ] << " " << QEdges[iQ+1] << endl;
 51        _h_tauc.add(QEdges[iQ], QEdges[iQ+1], book(dummy,1+iQ,1,1));
 52        _h_tau.add(QEdges[iQ], QEdges[iQ+1], book(dummy,8+iQ,1,1));
 53        _h_B.add(QEdges[iQ], QEdges[iQ+1], book(dummy,15+iQ,1,1));
 54        _h_rho.add(QEdges[iQ], QEdges[iQ+1], book(dummy,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      
 65      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 66      
 67      // The kinematic region covered by the analysis is defined by ranges of Q² and y
 68      if (dk.Q2() < 196 or dk.Q2() > 40000 or dk.y() < 0.1 or dk.y() > 0.7) {
 69      	// cerr << "Event out of the kinematic region covered by the analysis" << endl;
 70      	vetoEvent;
 71      }
 72
 73      double Q=sqrt(dk.Q2());
 74      _Nevt_after_cuts -> fill();
 75       for(int iQ = 0; iQ < 7; ++iQ) {
 76          if (inRange(Q, QEdges[iQ],QEdges[iQ+1])) {
 77         //  cout << " Q " << Q << " " << iQ << endl;
 78          _Nevt_after_cuts_Q[iQ] -> fill(); }
 79       }
 80
 81      
 82      // Boost to hadronic Breit frame
 83      const LorentzTransform breitboost = dk.boostBreit();
 84      
 85      const FinalState& fs = apply<FinalState>(event, "FS");
 86      
 87      /*Calculate event shape variables:
 88      thrust_num is \sum |pz_h|
 89      thrust_den is \sum |p_h|
 90      b_num is \sum |pt_h|
 91      sumE is the sum of energies
 92      (All sums run through the particles in the current hemisphere)
 93      */
 94      double thrust_num, thrust_den, b_num, sumE;
 95      thrust_num = thrust_den = b_num = sumE = 0;
 96      Vector3 sumMom;
 97      
 98      for (const Particle& p : fs.particles()) {
 99      	// Boost to Breit frame
100        const FourMomentum BreitMom = breitboost.transform(p.momentum());
101      	if (BreitMom.eta() < 0) {
102      		thrust_num += abs(BreitMom.pz());
103      		thrust_den += BreitMom.p();
104      		b_num += abs(BreitMom.pt());
105      		sumMom.operator+=(BreitMom.p3());
106      		sumE += BreitMom.E();
107      	}
108      }
109      
110      //The energy in the current hemisphere must exceed a certain value.
111      if (sumE <= Q/10.0) {
112      	// cerr << "Energy in the current hemisphere too low" << endl;
113      	vetoEvent;
114      }
115      
116      /* Comment from A. Galvan:
117      Thrust here is with respect to the z axis (tau in the paper), while the
118      one from Rivet projection is with respect to the maximum thrust 
119      axis (1 - tau_c in the paper)
120      */
121      
122      double thrust_mine = 1.0 - ((double)thrust_num)/((double)thrust_den);
123      double b_mine = ((double)b_num)/(2.0*(double)thrust_den);
124      double rho_num = thrust_den*thrust_den - sumMom.dot(sumMom);
125      double rho_mine = ((double)rho_num)/(4.0*(double)thrust_den*(double)thrust_den);
126      
127      const Thrust& thrCut = applyProjection<Thrust>(event, "ThrustCut");
128
129      //Fill histograms:
130      _h_tauc.fill(Q, 1.0 - thrCut.thrust());
131      _h_tau.fill(Q, thrust_mine);
132      _h_B.fill(Q, b_mine);
133      _h_rho.fill(Q, rho_mine);
134      
135      /* Comment from A. Galvan:
136       As for the C-parameter, my results did not fit the reference data at
137       all. The formula given in the paper is a bit ambiguous, because there is
138       a sum that runs through all pairs of particles h,h' and I was not sure
139       whether each pair should be counted twice (h,h' and h',h) or not, or if
140       the pair of a particle with itself (hh) should be considered. 
141       That is why I tried all of the possibilities, but none of them worked.
142      */
143    }
144
145
146    /// Normalise histograms etc., after the run
147    void finalize() {
148
149      //cerr << endl << endl << "Finalizing --------------------------------" << endl << endl;
150      
151      // cerr << "Cross section: " << crossSection() << endl;
152      
153      // double Nev = dbl(*_Nevt_after_cuts) ;
154      // cout << " N total " << Nev << endl; 
155
156
157      int iQ=0;
158      for (Histo1DPtr histo : _h_tauc.histos()) { 
159          double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
160          // cout << " Nev " << Nev << " " << iQ << endl;
161          if (Nev != 0) scale(histo, 1./Nev);
162          vector<YODA::HistoBin1D>& bins = histo -> bins();
163          for (auto & b : bins) b.scaleW(b.xWidth());
164          ++iQ;
165      }
166      
167      iQ=0;
168      for (Histo1DPtr histo : _h_tau.histos()) { 
169          double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
170          //cout << " Nev " << Nev << " " << iQ << endl;
171          if (Nev != 0) scale(histo, 1./Nev);
172          vector<YODA::HistoBin1D>& bins = histo -> bins();
173          for (auto & b : bins) b.scaleW(b.xWidth());
174          ++iQ;
175      }
176      
177      iQ = 0;
178      for (Histo1DPtr histo : _h_B.histos()) { 
179          double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
180          // cout << " Nev " << Nev << " " << iQ << endl;
181          if (Nev != 0) scale(histo, 1./Nev);
182          vector<YODA::HistoBin1D>& bins = histo -> bins();
183          for (auto & b : bins) b.scaleW(b.xWidth());
184          ++iQ;
185      }
186      
187      iQ = 0;
188      for (Histo1DPtr histo : _h_rho.histos()) { 
189          double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
190          // cout << " Nev " << Nev << " " << iQ << endl;
191          if (Nev != 0) scale(histo, 1./Nev);
192          vector<YODA::HistoBin1D>& bins = histo -> bins();
193          for (auto & b : bins) b.scaleW(b.xWidth());
194          ++iQ;
195      }
196
197    }
198
199    ///@}
200
201
202    /// @name Histograms
203    ///@{
204    map<string, Histo1DPtr> _h;
205    map<string, Profile1DPtr> _p;
206    map<string, CounterPtr> _c;
207    BinnedHistogram _h_tauc, _h_tau, _h_B, _h_rho;
208    CounterPtr _Nevt_after_cuts;
209    CounterPtr _Nevt_after_cuts_Q[7];
210    ///@}
211    const vector<double> QEdges {14., 16., 20., 30., 50., 70., 100., 200.};
212
213
214  };
215
216
217  RIVET_DECLARE_PLUGIN(H1_2006_I699835);
218
219}