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(_e["Avg1"], 10,1,1);
 50      book(_e["Avg2"], 11,1,1);
 51
 52      book(_h["QT1"], "TMP/QT1",refData(10,1,1));
 53      book(_h["QT2"], "TMP/QT2",refData(11,1,1));
 54
 55      book(_h["AvgT1"], "TMP/AvgT1",refData(10,1,1));
 56      book(_h["AvgT2"], "TMP/AvgT2",refData(11,1,1));
 57
 58      for (size_t iE = 0; iE < 8; ++iE) {
 59        book(_h["E"+to_string(iE)], 12+iE, 1, 1);
 60        book(_Nevt_after_cuts_E[iE], "TMP/Nevt_after_cuts_E"+to_string(iE));
 61      }
 62
 63      for (size_t iN = 0; iN < 6; ++iN) {
 64        book(_h["N"+to_string(iN)], 20+iN, 1, 1);
 65        book(_Nevt_after_cuts_N[iN], "TMP/Nevt_after_cuts_N"+to_string(iN));
 66      }
 67
 68      book(_h["MeanTest1"], "TMP/Meantest1",20,5.,6.14);
 69      book(_h["MeanTest2"], "TMP/Meantest2",50,19.,8000.);
 70    }
 71
 72
 73    /// Perform the per-event analysis
 74    void analyze(const Event& event) {
 75
 76      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 77      Particles particles = cfs.particles();
 78      const size_t numParticles = particles.size();
 79
 80      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 81      double Q2 = dk.Q2()/GeV2;
 82      double y= dk.y();
 83      double x= dk.x();
 84      double W2= dk.W2();
 85      double Q = sqrt(Q2);
 86
 87      if (y < 0.05 or y > 0.6 ) vetoEvent;
 88      if (W2 < 4400) vetoEvent ;
 89
 90      if (numParticles < 2) {
 91        MSG_DEBUG("Failed leptonic event cut");
 92        vetoEvent;
 93      }
 94
 95      for (int iQ = 0; iQ < 11; ++iQ) {
 96        if (inRange(sqrt(Q2), QEdges[iQ],QEdges[iQ+1])) {
 97          _Nevt_after_cuts_Q[iQ] -> fill();
 98        }
 99      }
100
101      if (Q > 5 && Q < 6.14)  _Nevt_after_cuts_E[0]->fill();
102      if (Q >19 && Q < 8000)  _Nevt_after_cuts_E[1]->fill();
103      if (Q2 >12 && Q2 <  15) _Nevt_after_cuts_E[2]->fill();
104      if (Q2 >15 && Q2 <  25) _Nevt_after_cuts_E[3]->fill();
105      if (Q2 >20 && Q2 <  40) _Nevt_after_cuts_E[4]->fill();
106      if (Q2 >40 && Q2 <  60) _Nevt_after_cuts_E[5]->fill();
107      if (Q2 >60 && Q2 <  80) _Nevt_after_cuts_E[6]->fill();
108      if (Q2 >80 && Q2 < 100) _Nevt_after_cuts_E[7]->fill();
109
110      if (Q2 >  12 && Q2 <  30 && x>6e-4 && x<2e-3)  _Nevt_after_cuts_N[0]->fill();
111      if (Q2 >  12 && Q2 <  30 && x>2e-3 && x<1e-2)  _Nevt_after_cuts_N[1]->fill();
112      if (Q2 >  30 && Q2 <  80 && x>6e-4 && x<2e-3)  _Nevt_after_cuts_N[2]->fill();
113      if (Q2 >  30 && Q2 <  80 && x>2e-3 && x<1e-2)  _Nevt_after_cuts_N[3]->fill();
114      if (Q2 > 100 && Q2 < 500 && x>2e-3 && x<1e-2)  _Nevt_after_cuts_N[4]->fill();
115      if (Q2 > 100 && Q2 < 500 && x>1e-2 && x<2e-1)  _Nevt_after_cuts_N[5]->fill();
116
117      if ( Q2 >12 && Q2 < 100 )  {
118        _Nevt_after_cuts_Qlow->fill();
119      }
120
121      if ( Q2 >100 && Q2 < 8000 )  {
122        _Nevt_after_cuts_QHigh->fill();
123      }
124
125      double multi=0;
126      double multiperevent1 = 0.0;
127      double multiperevent2 = 0.0;
128      double multiperevent3 = 0.0;
129      double multiperevent4 = 0.0;
130      double multiperevent5 = 0.0;
131      double multiperevent6 = 0.0;
132
133
134      const LorentzTransform Breitboost = dk.boostBreit();
135
136
137      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
138        const Particle& p = particles[ip1];
139
140        const FourMomentum BreMom = Breitboost.transform(p.momentum());
141
142        if ( BreMom.pz() > 0. ) continue;
143        double pcal= std::sqrt(BreMom.px2() + BreMom.py2()+ BreMom.pz2());
144        double xp = 2*pcal/(sqrt(Q2));
145
146        double E = std::sqrt((.27*.27)+(pcal*pcal));
147        double dp = 4*M_PI*pcal*pcal;
148        double dE = pcal;
149        double factor = dE/dp;
150
151        if ( Q2 >12 && Q2 < 100 )  {
152          _h["xp"] -> fill(xp);
153          _h["xi"] -> fill(log(1/(xp)));
154        }
155
156        if ( Q2 >100 && Q2 < 8000 )  {
157          _h["xpQgt"] -> fill(xp);
158          _h["xiQgt"] -> fill(log(1/(xp)));
159        }
160
161        if (Q  >  5 && Q < 6.14)  _h["E0"]->fill(E, factor);
162        if (Q  > 19 && Q < 8000)  _h["E1"]->fill(E, factor);
163        if (Q2 > 12 && Q2 <  15)  _h["E2"]->fill(E, factor);
164        if (Q2 > 15 && Q2 <  25)  _h["E3"]->fill(E, factor);
165        if (Q2 > 20 && Q2 <  40)  _h["E4"]->fill(E, factor);
166        if (Q2 > 40 && Q2 <  60)  _h["E5"]->fill(E, factor);
167        if (Q2 > 60 && Q2 <  80)  _h["E6"]->fill(E, factor);
168        if (Q2 > 80 && Q2 < 100)  _h["E7"]->fill(E, factor);
169
170        _h_Q2_xp->fill(sqrt(Q2),xp);
171
172        multi = multi+1;
173
174        if (Q2 >  12 && Q2 <  30 && x>6e-4 && x<2e-3)   multiperevent1 = multiperevent1+1;
175        if (Q2 >  12 && Q2 <  30 && x>2e-3 && x<1e-2)   multiperevent2 = multiperevent2+1;
176        if (Q2 >  30 && Q2 <  80 && x>6e-4 && x<2e-3)   multiperevent3 = multiperevent3+1;
177        if (Q2 >  30 && Q2 <  80 && x>2e-3 && x<1e-2)   multiperevent4 = multiperevent4+1;
178        if (Q2 > 100 && Q2 < 500 && x>2e-3 && x<1e-2)   multiperevent5 = multiperevent5+1;
179        if (Q2 > 100 && Q2 < 500 && x>1e-2 && x<2e-1)   multiperevent6 = multiperevent6+1;
180
181      }
182
183      if (Q2 >  12 && Q2 <  30 && x>6e-4 && x<2e-3)   _h["N0"]->fill(multiperevent1);
184      if (Q2 >  12 && Q2 <  30 && x>2e-3 && x<1e-2)   _h["N1"]->fill(multiperevent2);
185      if (Q2 >  30 && Q2 <  80 && x>6e-4 && x<2e-3)   _h["N2"]->fill(multiperevent3);
186      if (Q2 >  30 && Q2 <  80 && x>2e-3 && x<1e-2)   _h["N3"]->fill(multiperevent4);
187      if (Q2 > 100 && Q2 < 500 && x>2e-3 && x<1e-2)   _h["N4"]->fill(multiperevent5);
188      if (Q2 > 100 && Q2 < 500 && x>1e-2 && x<2e-1)   _h["N5"]->fill(multiperevent6);
189
190      if (Q2 > 12) {
191       _h["AvgT2"] -> fill(sqrt(Q2), multi);
192       _h["QT2"] -> fill(sqrt(Q2));
193      }
194
195      _h["AvgT1"] -> fill(sqrt(Q2), multi);
196      _h["QT1"] -> fill(sqrt(Q2));
197
198      _h["MeanTest1"] -> fill(sqrt(Q2));
199      _h["MeanTest2"] -> fill(sqrt(Q2));
200
201    }
202
203
204    /// Normalise histograms etc., after the run
205    void finalize() {
206      MSG_DEBUG("Nevt Qlow " << dbl(*_Nevt_after_cuts_Qlow));
207      scale(_h["xp"], 1.0/ *_Nevt_after_cuts_Qlow);
208      scale(_h["xi"], 1.0/ *_Nevt_after_cuts_Qlow);
209      scale(_h["xpQgt"], 1.0/ *_Nevt_after_cuts_QHigh);
210      scale(_h["xiQgt"], 1.0/ *_Nevt_after_cuts_QHigh);
211
212      for(int iE=0 ; iE< 8 ; ++iE){
213        scale(_h["E"+to_string(iE)], 1.0/ *_Nevt_after_cuts_E[iE]);
214      }
215
216      for(int iN=0 ; iN< 8 ; ++iN){
217        scale(_h["N"+to_string(iN)], 1.0/ *_Nevt_after_cuts_N[iN]);
218      }
219
220      divide(_h["AvgT1"], _h["QT1"],_e["Avg1"] );
221      divide(_h["AvgT2"], _h["QT2"],_e["Avg2"] );
222
223
224
225      int iQ = 0;
226      double mean;
227      for (auto& histo :_h_Q2_xp->bins()) {
228        const double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
229        if (Nev != 0) scale(histo, 1./Nev);
230
231        for (size_t iP = 0; iP < iPmax; ++iP) {
232          mean = histo->bin(iP).sumW() ;
233          double mean_err = mean/100;
234          _e["Qxp"+to_string(iP)]->bin(iP+1).set(mean, mean_err);
235        }
236        ++iQ;
237      }
238      if(_h["MeanTest1"]->numEntries(false)>0 && _h["MeanTest1"]->effNumEntries(false)>0) {
239        const double x1 = _h["MeanTest1"]->xMean(false);
240        MSG_DEBUG("Mean of low Q = " << x1);
241      }
242      if(_h["MeanTest2"]->numEntries(false)>0 && _h["MeanTest2"]->effNumEntries(false)>0) {
243        const double x2 = _h["MeanTest2"]->xMean(false);
244        MSG_DEBUG("Mean of High Q = " << x2);
245      }
246    }
247
248    /// @}
249
250
251  private:
252
253    /// @name Histograms
254    /// @{
255
256    CounterPtr _Nevt_after_cuts_Qlow,_Nevt_after_cuts_QHigh,_Nevt_after_cuts_E[8],_Nevt_after_cuts_N[6];
257    CounterPtr _Nevt_after_cuts_Q[12];
258    map<string, Histo1DPtr> _h;
259    map<string, Estimate1DPtr> _e;
260    Estimate1DPtr _h_pt_06_ratio;
261
262    Histo1DGroupPtr _h_Q2_xp;
263
264    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};
265    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};
266    const vector<double> xp_range{0.02, 0.05, 0.10, 0.20, 0.3, 0.4, 0.5, 0.7};
267    const size_t iPmax = 7;
268
269    /// @}
270
271  };
272
273
274  RIVET_DECLARE_PLUGIN(H1_1997_I445116);
275
276}