1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/DISLepton.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Evolution of e p fragmentation and multiplicity distributions in the Breit frame (H1)
 12  class H1_1997_I445116 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1997_I445116);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      declare(DISKinematics(), "Kinematics");
 26      const DISLepton dl;
 27      declare(ChargedFinalState(dl.remainingFinalState()), "CFS");
 28
 29      book(_Nevt_after_cuts_Qlow, "TMP/Nevt_after_cuts_Qlow");
 30      book(_Nevt_after_cuts_QHigh, "TMP/Nevt_after_cuts_QHigh");
 31
 32      book(_h["xp"], 1, 1, 1);
 33      book(_h["xpQgt"], 1, 1, 2);
 34      book(_h["xi"], 2, 1, 1);
 35      book(_h["xiQgt"], 2, 1, 2);
 36
 37      book(_h_Q2_xp, QEdges);
 38      for (auto& b : _h_Q2_xp->bins()) {
 39        const size_t iQ = b.index()-1;
 40        const string suff = to_string(iQ);
 41        book(b, "TMP/xpQ"+suff, xp_range);
 42        book(_Nevt_after_cuts_Q[iQ], "TMP/Nevt_after_cuts_Q"+suff);
 43      }
 44
 45      for (size_t iP = 0 ; iP < iPmax; ++iP) {
 46        book(_e["Qxp"+to_string(iP)], 3+iP, 1, 1);
 47      }
 48
 49      book(_p["Avg1"], 10,1,1);
 50      book(_p["Avg2"], 11,1,1);
 51
 52      for (size_t iE = 0; iE < 8; ++iE) {
 53        book(_h["E"+to_string(iE)], 12+iE, 1, 1);
 54        book(_Nevt_after_cuts_E[iE], "TMP/Nevt_after_cuts_E"+to_string(iE));
 55      }
 56
 57      for (size_t iN = 0; iN < 6; ++iN) {
 58        book(_h["N"+to_string(iN)], 20+iN, 1, 1);
 59        book(_Nevt_after_cuts_N[iN], "TMP/Nevt_after_cuts_N"+to_string(iN));
 60      }
 61
 62      book(_h["MeanTest1"], "TMP/Meantest1",20,5.,6.14);
 63      book(_h["MeanTest2"], "TMP/Meantest2",50,19.,8000.);
 64    }
 65
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69
 70      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 71      Particles particles = cfs.particles();
 72      const size_t numParticles = particles.size();
 73
 74      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 75      double Q2 = dk.Q2()/GeV2;
 76      double y= dk.y();
 77      double x= dk.x();
 78      double W2= dk.W2();
 79      double Q = sqrt(Q2);
 80
 81      if (y < 0.05 or y > 0.6 ) vetoEvent;
 82      if (W2 < 4400) vetoEvent ;
 83
 84      if (numParticles < 2) {
 85        MSG_DEBUG("Failed leptonic event cut");
 86        vetoEvent;
 87      }
 88
 89      for (int iQ = 0; iQ < 11; ++iQ) {
 90        if (inRange(sqrt(Q2), QEdges[iQ],QEdges[iQ+1])) {
 91          _Nevt_after_cuts_Q[iQ] -> fill();
 92        }
 93      }
 94
 95      if (Q > 5 && Q < 6.14)  _Nevt_after_cuts_E[0]->fill();
 96      if (Q >19 && Q < 8000)  _Nevt_after_cuts_E[1]->fill();
 97      if (Q2 >12 && Q2 <  15) _Nevt_after_cuts_E[2]->fill();
 98      if (Q2 >15 && Q2 <  25) _Nevt_after_cuts_E[3]->fill();
 99      if (Q2 >20 && Q2 <  40) _Nevt_after_cuts_E[4]->fill();
100      if (Q2 >40 && Q2 <  60) _Nevt_after_cuts_E[5]->fill();
101      if (Q2 >60 && Q2 <  80) _Nevt_after_cuts_E[6]->fill();
102      if (Q2 >80 && Q2 < 100) _Nevt_after_cuts_E[7]->fill();
103
104      if (Q2 >  12 && Q2 <  30 && x>6e-4 && x<2e-3)  _Nevt_after_cuts_N[0]->fill();
105      if (Q2 >  12 && Q2 <  30 && x>2e-3 && x<1e-2)  _Nevt_after_cuts_N[1]->fill();
106      if (Q2 >  30 && Q2 <  80 && x>6e-4 && x<2e-3)  _Nevt_after_cuts_N[2]->fill();
107      if (Q2 >  30 && Q2 <  80 && x>2e-3 && x<1e-2)  _Nevt_after_cuts_N[3]->fill();
108      if (Q2 > 100 && Q2 < 500 && x>2e-3 && x<1e-2)  _Nevt_after_cuts_N[4]->fill();
109      if (Q2 > 100 && Q2 < 500 && x>1e-2 && x<2e-1)  _Nevt_after_cuts_N[5]->fill();
110
111      if ( Q2 >12 && Q2 < 100 )  {
112        _Nevt_after_cuts_Qlow->fill();
113      }
114
115      if ( Q2 >100 && Q2 < 8000 )  {
116        _Nevt_after_cuts_QHigh->fill();
117      }
118
119      double multi=0;
120      double multiperevent1 = 0.0;
121      double multiperevent2 = 0.0;
122      double multiperevent3 = 0.0;
123      double multiperevent4 = 0.0;
124      double multiperevent5 = 0.0;
125      double multiperevent6 = 0.0;
126
127
128      const LorentzTransform Breitboost = dk.boostBreit();
129
130
131      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
132        const Particle& p = particles[ip1];
133
134        const FourMomentum BreMom = Breitboost.transform(p.momentum());
135
136        if ( BreMom.pz() > 0. ) continue;
137        double pcal= std::sqrt(BreMom.px2() + BreMom.py2()+ BreMom.pz2());
138        double xp = 2*pcal/(sqrt(Q2));
139
140        double E = std::sqrt((.27*.27)+(pcal*pcal));
141        double dp = 4*M_PI*pcal*pcal;
142        double dE = pcal;
143        double factor = dE/dp;
144
145        if ( Q2 >12 && Q2 < 100 )  {
146          _h["xp"] -> fill(xp);
147          _h["xi"] -> fill(log(1/(xp)));
148        }
149
150        if ( Q2 >100 && Q2 < 8000 )  {
151          _h["xpQgt"] -> fill(xp);
152          _h["xiQgt"] -> fill(log(1/(xp)));
153        }
154
155        if (Q  >  5 && Q < 6.14)  _h["E0"]->fill(E, factor);
156        if (Q  > 19 && Q < 8000)  _h["E1"]->fill(E, factor);
157        if (Q2 > 12 && Q2 <  15)  _h["E2"]->fill(E, factor);
158        if (Q2 > 15 && Q2 <  25)  _h["E3"]->fill(E, factor);
159        if (Q2 > 20 && Q2 <  40)  _h["E4"]->fill(E, factor);
160        if (Q2 > 40 && Q2 <  60)  _h["E5"]->fill(E, factor);
161        if (Q2 > 60 && Q2 <  80)  _h["E6"]->fill(E, factor);
162        if (Q2 > 80 && Q2 < 100)  _h["E7"]->fill(E, factor);
163
164        _h_Q2_xp->fill(sqrt(Q2),xp);
165
166        multi = multi+1;
167
168        if (Q2 >  12 && Q2 <  30 && x>6e-4 && x<2e-3)   multiperevent1 = multiperevent1+1;
169        if (Q2 >  12 && Q2 <  30 && x>2e-3 && x<1e-2)   multiperevent2 = multiperevent2+1;
170        if (Q2 >  30 && Q2 <  80 && x>6e-4 && x<2e-3)   multiperevent3 = multiperevent3+1;
171        if (Q2 >  30 && Q2 <  80 && x>2e-3 && x<1e-2)   multiperevent4 = multiperevent4+1;
172        if (Q2 > 100 && Q2 < 500 && x>2e-3 && x<1e-2)   multiperevent5 = multiperevent5+1;
173        if (Q2 > 100 && Q2 < 500 && x>1e-2 && x<2e-1)   multiperevent6 = multiperevent6+1;
174
175      }
176
177      if (Q2 >  12 && Q2 <  30 && x>6e-4 && x<2e-3)   _h["N0"]->fill(multiperevent1);
178      if (Q2 >  12 && Q2 <  30 && x>2e-3 && x<1e-2)   _h["N1"]->fill(multiperevent2);
179      if (Q2 >  30 && Q2 <  80 && x>6e-4 && x<2e-3)   _h["N2"]->fill(multiperevent3);
180      if (Q2 >  30 && Q2 <  80 && x>2e-3 && x<1e-2)   _h["N3"]->fill(multiperevent4);
181      if (Q2 > 100 && Q2 < 500 && x>2e-3 && x<1e-2)   _h["N4"]->fill(multiperevent5);
182      if (Q2 > 100 && Q2 < 500 && x>1e-2 && x<2e-1)   _h["N5"]->fill(multiperevent6);
183
184      if (Q2 > 12) _p["Avg2"] -> fill(sqrt(Q2), multi);
185
186      _p["Avg1"] -> fill(sqrt(Q2), multi);
187
188      _h["MeanTest1"] -> fill(sqrt(Q2));
189      _h["MeanTest2"] -> fill(sqrt(Q2));
190
191    }
192
193
194    /// Normalise histograms etc., after the run
195    void finalize() {
196      MSG_DEBUG("Nevt Qlow " << dbl(*_Nevt_after_cuts_Qlow));
197      scale(_h["xp"], 1.0/ *_Nevt_after_cuts_Qlow);
198      scale(_h["xi"], 1.0/ *_Nevt_after_cuts_Qlow);
199      scale(_h["xpQgt"], 1.0/ *_Nevt_after_cuts_QHigh);
200      scale(_h["xiQgt"], 1.0/ *_Nevt_after_cuts_QHigh);
201
202      for(int iE=0 ; iE< 8 ; ++iE){
203        scale(_h["E"+to_string(iE)], 1.0/ *_Nevt_after_cuts_E[iE]);
204      }
205
206      for(int iN=0 ; iN< 6 ; ++iN){
207        scale(_h["N"+to_string(iN)], 1.0/ *_Nevt_after_cuts_N[iN]);
208      }
209
210
211
212
213      int iQ = 0;
214      for (auto& histo :_h_Q2_xp->bins()) {
215        const double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
216        if (Nev != 0) scale(histo, 1./Nev);
217
218        for (size_t iP = 0; iP < iPmax; ++iP) {
219          double mean      = histo->bin(iP+1).sumW()       /histo->bin(iP+1).xWidth();
220          double mean_err = sqrt( histo->bin(iP+1).sumW2())/histo->bin(iP+1).xWidth();
221          _e["Qxp"+to_string(iP)]->bin(iQ+1).set(mean, mean_err);
222        }
223        ++iQ;
224      }
225      if(_h["MeanTest1"]->numEntries(false)>0 && _h["MeanTest1"]->effNumEntries(false)>0) {
226        const double x1 = _h["MeanTest1"]->xMean(false);
227        MSG_DEBUG("Mean of low Q = " << x1);
228      }
229      if(_h["MeanTest2"]->numEntries(false)>0 && _h["MeanTest2"]->effNumEntries(false)>0) {
230        const double x2 = _h["MeanTest2"]->xMean(false);
231        MSG_DEBUG("Mean of High Q = " << x2);
232      }
233    }
234
235    /// @}
236
237
238  private:
239
240    /// @name Histograms
241    /// @{
242
243    CounterPtr _Nevt_after_cuts_Qlow,_Nevt_after_cuts_QHigh,_Nevt_after_cuts_E[8],_Nevt_after_cuts_N[6];
244    CounterPtr _Nevt_after_cuts_Q[12];
245    map<string, Histo1DPtr> _h;
246    map<string, Estimate1DPtr> _e;
247    map<string, Profile1DPtr> _p;
248    Estimate1DPtr _h_pt_06_ratio;
249
250    Histo1DGroupPtr _h_Q2_xp;
251
252    const vector<double> QEdges {3.17, 3.915, 4.72, 6.13, 7.635, 8.85, 10.81, 13.415, 16.365, 21.745, 33.835, 50.745};
253    const vector<double> AvgEdges {3.13, 3.875, 4.675, 6.065, 7.55, 8.76, 9.785, 14.675, 16.25, 21.45, 33.06, 49.26};
254    const vector<double> xp_range{0.02, 0.05, 0.10, 0.20, 0.3, 0.4, 0.5, 0.7};
255    const size_t iPmax = 7;
256
257    /// @}
258
259  };
260
261
262  RIVET_DECLARE_PLUGIN(H1_1997_I445116);
263
264}