1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/DISKinematics.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Charged particle multiplicities in deep inelastic scattering at HERA (H1)
 12  class H1_1996_I422230 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1996_I422230);
 17
 18
 19    /// @name Analysis methods
 20    ///@{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Initialise and register projections
 26      const DISLepton disl;
 27      declare(disl, "Lepton");
 28      declare(DISKinematics(), "Kinematics");
 29
 30      // The basic final-state projection:
 31      // all final-state particles within
 32      // the given eta acceptance
 33      const FinalState fs;
 34      declare(fs,"FS");
 35      const ChargedFinalState cfs(disl.remainingFinalState());
 36      declare(cfs,"CFS");
 37
 38      //binned histograms to count for the multiplicity
 39      for (size_t ix = 0; ix < 4; ++ix) {
 40        book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts"+ to_string(ix));
 41      }
 42
 43      const vector<double> WEdges{80., 115, 150., 185., 220.};
 44      book(_g["mult1"], WEdges);
 45      book(_g["mult2"], WEdges);
 46      book(_g["mult3"], WEdges);
 47      book(_g["mult4"], WEdges);
 48      book(_g["mult_all"], WEdges);
 49      book(_g["mult10_all"], WEdges);
 50      book(_g["mult11_all"], WEdges);
 51      book(_g["mult12_all"], WEdges);
 52      for (size_t iW=0; iW < _g["mult1"]->numBins(); ++iW) {
 53        book(_g["mult1"]->bin(iW+1), iW+1, 1, 1);
 54        book(_g["mult2"]->bin(iW+1), iW+1, 1, 2);
 55        book(_g["mult3"]->bin(iW+1), iW+1, 1, 3);
 56        book(_g["mult4"]->bin(iW+1), iW+1, 1, 4);
 57        const auto& ref = refData<YODA::BinnedEstimate<int>>(4,1,4);
 58        book(_g["mult_all"]->bin(iW+1),   "TMP/dummy" + to_string(iW), ref);
 59        book(_g["mult10_all"]->bin(iW+1), "TMP/dummy1"+ to_string(iW), ref);
 60        book(_g["mult11_all"]->bin(iW+1), "TMP/dummy2"+ to_string(iW), ref);
 61        book(_g["mult12_all"]->bin(iW+1), "TMP/dummy3"+ to_string(iW), ref);
 62      }
 63
 64      //histograms for the statistical moments and the mean
 65      book(_e["mean0"],5,1,1);
 66      book(_e["D2_0"],5,1,2);
 67      book(_e["D3_0"],5,1,3);
 68      book(_e["D4_0"],5,1,4);
 69      book(_e["C2_0"],5,1,5);
 70      book(_e["C3_0"],5,1,6);
 71      book(_e["C4_0"],5,1,7);
 72      book(_e["R2_0"],5,1,8);
 73      book(_e["R3_0"],5,1,9);
 74
 75      book(_e["mean12"],6,1,1);
 76      book(_e["D2_12"],6,1,2);
 77      book(_e["D3_12"],6,1,3);
 78      book(_e["D4_12"],6,1,4);
 79      book(_e["C2_12"],6,1,5);
 80      book(_e["C3_12"],6,1,6);
 81      book(_e["C4_12"],6,1,7);
 82      book(_e["R2_12"],6,1,8);
 83      book(_e["R3_12"],6,1,9);
 84      book(_e["K3_12"],6,1,10);
 85
 86      book(_e["mean13"],7,1,1);
 87      book(_e["D2_13"],7,1,2);
 88      book(_e["D3_13"],7,1,3);
 89      book(_e["D4_13"],7,1,4);
 90      book(_e["C2_13"],7,1,5);
 91      book(_e["C3_13"],7,1,6);
 92      book(_e["C4_13"],7,1,7);
 93      book(_e["R2_13"],7,1,8);
 94      book(_e["R3_13"],7,1,9);
 95      book(_e["K3_13"],7,1,10);
 96
 97      book(_e["mean14"],8,1,1);
 98      book(_e["D2_14"],8,1,2);
 99      book(_e["D3_14"],8,1,3);
100      book(_e["D4_14"],8,1,4);
101      book(_e["C2_14"],8,1,5);
102      book(_e["C3_14"],8,1,6);
103      book(_e["C4_14"],8,1,7);
104      book(_e["R2_14"],8,1,8);
105      book(_e["R3_14"],8,1,9);
106
107      book(_e["mean15"],9,1,1);
108      book(_e["D2_15"],9,1,2);
109      book(_e["D3_15"],9,1,3);
110      book(_e["D4_15"],9,1,4);
111      book(_e["C2_15"],9,1,5);
112      book(_e["C3_15"],9,1,6);
113      book(_e["C4_15"],9,1,7);
114      book(_e["R2_15"],9,1,8);
115      book(_e["R3_15"],9,1,9);
116
117      book(_e["mean23"],10,1,1);
118      book(_e["D2_23"],10,1,2);
119      book(_e["D3_23"],10,1,3);
120      book(_e["D4_23"],10,1,4);
121      book(_e["C2_23"],10,1,5);
122      book(_e["C3_23"],10,1,6);
123      book(_e["C4_23"],10,1,7);
124      book(_e["R2_23"],10,1,8);
125      book(_e["R3_23"],10,1,9);
126      book(_e["K3_23"],10,1,10);
127
128      book(_e["mean34"],11,1,1);
129      book(_e["D2_34"],11,1,2);
130      book(_e["D3_34"],11,1,3);
131      book(_e["D4_34"],11,1,4);
132      book(_e["C2_34"],11,1,5);
133      book(_e["C3_34"],11,1,6);
134      book(_e["C4_34"],11,1,7);
135      book(_e["R2_34"],11,1,8);
136      book(_e["R3_34"],11,1,9);
137      book(_e["K3_34"],11,1,10);
138
139      book(_e["mean45"],12,1,1);
140      book(_e["D2_45"],12,1,2);
141      book(_e["D3_45"],12,1,3);
142      book(_e["D4_45"],12,1,4);
143      book(_e["C2_45"],12,1,5);
144      book(_e["C3_45"],12,1,6);
145      book(_e["C4_45"],12,1,7);
146      book(_e["R2_45"],12,1,8);
147      book(_e["R3_45"],12,1,9);
148      book(_e["K3_45"],12,1,10);
149    }
150
151
152    /// Perform the per-event analysis
153    void analyze(const Event& event) {
154
155      //apply the final state of all particles and the one for only charged particles
156      const FinalState& fs = apply<FinalState>(event, "FS");
157      const size_t numParticles = fs.particles().size();
158      const ChargedFinalState& cfs = apply<ChargedFinalState>(event,"CFS");
159      const Particles& particles = cfs.particles();
160
161      //because it does not make sense to have a collision with the numparticles is less than two,we use the vetoEvent so that if there is an event like this it does not run the analysis and goes to the next one
162
163      if (numParticles<2) {
164        MSG_DEBUG("Failed leptonic cut");
165        vetoEvent;
166      }
167
168      //apply DIS kinematics in the event
169      const DISKinematics& dk = apply<DISKinematics>(event,"Kinematics");
170      //const DISLepton& dl = apply<DISLepton>(event,"Lepton");
171
172      //get the DIS Kinematics but not in a loop because they are variables that descrube the type of event not the particles
173      double Q2 = dk.Q2();
174      double W2 = dk.W2();
175      double W = std::sqrt(W2);
176      bool cut1 = W<220. && W>80.;
177      if (!cut1) vetoEvent;
178      bool cut = Q2<1000. && Q2>10. && W>80. && W<220.;
179      if (!cut) vetoEvent;
180
181      double Efwd = 0. ;
182      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
183        const Particle& p = particles[ip1];
184        double theta = p.theta()/degree ;
185        if ( inRange(theta,4.4,15.) ) {
186          Efwd = Efwd + p.E() ;
187        }
188      }
189
190      bool cut_fwd = Efwd > 0.5 && dk.beamLepton().E() > 12. ;
191      if (!cut_fwd) vetoEvent ;
192
193      const size_t idx = _g["mult1"]->binAt(W).index();
194      if (0 < idx && idx < 5) {
195        _Nevt_after_cuts[idx-1]->fill();
196      }
197
198
199      //boost to the hadronic centre of mass frame
200      const LorentzTransform hcmboost = dk.boostHCM();
201      int kall = 0.;
202      int k1 = 0.;
203      int k2 = 0.;
204      int k3 = 0.;
205      int k4 = 0.;
206      int k10 = 0.;
207      int k11 = 0.;
208      int k12 = 0.;
209      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
210        const Particle& p = particles[ip1];
211        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
212        const double etahcm_char = hcmMom.eta();
213        if (etahcm_char>0.)  ++kall;
214        if (etahcm_char>1. && etahcm_char<2.) {
215          ++k1;
216        }
217        if (etahcm_char>1. && etahcm_char<3.) {
218          ++k2;
219        }
220        if (etahcm_char>1. && etahcm_char<4.) {
221          ++k3;
222        }
223        if (etahcm_char>1. && etahcm_char<5.) {
224          ++k4;
225        }
226        if (etahcm_char>2. && etahcm_char<3.) {
227          ++k10;
228        }
229        if (etahcm_char>3. && etahcm_char<4.) {
230          ++k11;
231        }
232        if (etahcm_char>4. && etahcm_char<5.) {
233          ++k12;
234        }
235
236      }
237      _g["mult_all"]->fill(W,kall);
238      _g["mult10_all"]->fill(W,k10);
239      _g["mult11_all"]->fill(W,k11);
240      _g["mult12_all"]->fill(W,k12);
241      _g["mult1"]->fill(W,k1);
242      _g["mult2"]->fill(W,k2);
243      _g["mult3"]->fill(W,k3);
244      _g["mult4"]->fill(W,k4);
245
246    }
247
248    /// Normalise histograms etc., after the run
249    void finalize() {
250      double iq;
251      double mean, dispersion, cq, R2, R3, K3;
252      for (auto& histo : _g["mult_all"]->bins()) {
253        iq = 2;
254        size_t iW = histo.index();
255        _histo_to_moments(histo, iq , mean, dispersion, cq, R2, R3, K3);
256
257        // just to have some values, needs to be corrected
258        _e["mean0"]->bin(iW).set(mean, 0.5*dispersion);
259        _e["D2_0"]->bin(iW).set(dispersion, dispersion/mean);
260        _e["C2_0"]->bin(iW).set(cq, cq/mean);
261        _e["R2_0"]->bin(iW).set(R2, R2/mean);
262        _e["R3_0"]->bin(iW).set(R3, R3/mean);
263        iq = 3 ;
264        dispersion = 0 ;
265        cq = 0;
266        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
267        _e["D3_0"]->bin(iW).set(dispersion, dispersion/mean);
268        _e["C3_0"]->bin(iW).set(cq, cq/mean);
269        iq = 4 ;
270        dispersion = 0 ;
271        cq = 0;
272        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
273        _e["D4_0"]->bin(iW).set(dispersion, dispersion/mean);
274        _e["C4_0"]->bin(iW).set(cq, cq/mean);
275      }
276      mean = 0.;
277      dispersion = 0.;
278      cq = 0.;
279      R2 = 0;
280      R3 = 0.;
281      K3 = 0.;
282      for (auto& histo : _g["mult1"]->bins()) {
283        // use scale to get Prob in %, and use counter to count events after cuts
284        size_t iW = histo.index();
285        scale(histo, 100.0/ *_Nevt_after_cuts[iW-1]);
286        //
287        iq = 2 ;
288        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
289
290        _e["mean12"]->bin(iW).set(mean, 0.5*dispersion);
291        _e["D2_12"]->bin(iW).set(dispersion, dispersion/mean);
292        _e["C2_12"]->bin(iW).set(cq, cq/mean);
293        _e["R2_12"]->bin(iW).set(R2, R2/mean);
294        _e["R3_12"]->bin(iW).set(R3, R3/mean);
295        _e["K3_12"]->bin(iW).set(K3, K3/mean);
296        iq = 3 ;
297        dispersion = 0 ;
298        cq = 0;
299        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
300        _e["D3_12"]->bin(iW).set(dispersion, dispersion/mean);
301        _e["C3_12"]->bin(iW).set(cq, cq/mean);
302        iq = 4 ;
303        dispersion = 0 ;
304        cq = 0;
305        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
306        _e["D4_12"]->bin(iW).set(dispersion, dispersion/mean);
307        _e["C4_12"]->bin(iW).set(cq, cq/mean);
308      }
309      mean = 0.;
310      dispersion = 0.;
311      cq = 0.;
312      R2 = 0;
313      R3 = 0.;
314      K3 = 0.;
315
316      for (auto& histo : _g["mult2"]->bins()) {
317        // use scale to get Prob in %, and use counter to count events after cuts
318        size_t iW = histo.index();
319        scale(histo, 100.0/ *_Nevt_after_cuts[iW-1]);
320        //
321        iq = 2 ;
322        _histo_to_moments(histo, iq, mean, dispersion,cq, R2, R3, K3);
323
324        _e["mean13"]->bin(iW).set(mean, 0.5*dispersion);
325        _e["D2_13"]->bin(iW).set(dispersion, dispersion/mean);
326        _e["C2_13"]->bin(iW).set(cq, cq/mean);
327        _e["R2_13"]->bin(iW).set(R2, R2/mean);
328        _e["R3_13"]->bin(iW).set(R3, R3/mean);
329        _e["K3_13"]->bin(iW).set(K3, K3/mean);
330        iq = 3 ;
331        dispersion = 0 ;
332        cq = 0;
333        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
334        _e["D3_13"]->bin(iW).set(dispersion, dispersion/mean);
335        _e["C3_13"]->bin(iW).set(cq, cq/mean);
336        iq = 4 ;
337        dispersion = 0 ;
338        cq = 0;
339        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
340        _e["D4_13"]->bin(iW).set(dispersion, dispersion/mean);
341        _e["C4_13"]->bin(iW).set(cq, cq/mean);
342      }
343      mean = 0.;
344      dispersion = 0.;
345      cq = 0.;
346      R2 = 0;
347      R3 = 0.;
348      K3 = 0.;
349
350      for (auto& histo : _g["mult3"]->bins()) {
351        // use scale to get Prob in %, and use counter to count events after cuts
352        size_t iW = histo.index();
353        scale(histo, 100.0/ *_Nevt_after_cuts[iW-1]);
354        //
355        iq = 2 ;
356        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
357        _e["mean14"]->bin(iW).set(mean, 0.5*dispersion);
358        _e["D2_14"]->bin(iW).set(dispersion, dispersion/mean);
359        _e["C2_14"]->bin(iW).set(cq, cq/mean);
360        _e["R2_14"]->bin(iW).set(R2, R2/mean);
361        _e["R3_14"]->bin(iW).set(R3, R3/mean);
362        iq = 3 ;
363        dispersion = 0 ;
364        cq = 0;
365        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
366        _e["D3_14"]->bin(iW).set(dispersion, dispersion/mean);
367        _e["C3_14"]->bin(iW).set(cq, cq/mean);
368        iq = 4 ;
369        dispersion = 0 ;
370        cq = 0;
371        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3) ;
372        _e["D4_14"]->bin(iW).set(dispersion, dispersion/mean);
373        _e["C4_14"]->bin(iW).set(cq, cq/mean);
374      }
375      mean = 0.;
376      dispersion = 0.;
377      cq = 0.;
378      R2 = 0;
379      R3 = 0.;
380      K3 = 0.;
381
382      for (auto& histo : _g["mult4"]->bins()) {
383        // use scale to get Prob in %, and use counter to count events after cuts
384        size_t iW = histo.index();
385        scale(histo, 100.0/ *_Nevt_after_cuts[iW-1]);
386        //
387        iq = 2 ;
388        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
389        _e["mean15"]->bin(iW).set(mean, 0.5*dispersion);
390        _e["D2_15"]->bin(iW).set(dispersion, dispersion/mean);
391        _e["C2_15"]->bin(iW).set(cq, cq/mean);
392        _e["R2_15"]->bin(iW).set(R2, R2/mean);
393        _e["R3_15"]->bin(iW).set(R3, R3/mean);
394        iq = 3 ;
395        dispersion = 0 ;
396        cq = 0;
397        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
398        _e["D3_15"]->bin(iW).set(dispersion, dispersion/mean);
399        _e["C3_15"]->bin(iW).set(cq, cq/mean);
400        iq = 4 ;
401        dispersion = 0 ;
402        cq = 0;
403        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
404        _e["D4_15"]->bin(iW).set(dispersion, dispersion/mean);
405        _e["C4_15"]->bin(iW).set(cq, cq/mean);
406      }
407      mean = 0.;
408      dispersion = 0.;
409      cq = 0.;
410      R2 = 0;
411      R3 = 0.;
412      K3 = 0.;
413
414      for (auto& histo : _g["mult10_all"]->bins()) {
415        // use scale to get Prob in %, and use counter to count events after cuts
416        size_t iW = histo.index();
417        scale(histo, 100.0/ *_Nevt_after_cuts[iW-1]);
418        //
419        iq = 2 ;
420        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
421        _e["mean23"]->bin(iW).set(mean, 0.5*dispersion);
422        _e["D2_23"]->bin(iW).set(dispersion, dispersion/mean);
423        _e["C2_23"]->bin(iW).set(cq, cq/mean);
424        _e["R2_23"]->bin(iW).set(R2, R2/mean);
425        _e["R3_23"]->bin(iW).set(R3, R3/mean);
426        _e["K3_23"]->bin(iW).set(K3, K3/mean);
427        iq = 3 ;
428        dispersion = 0 ;
429        cq = 0;
430        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
431        _e["D3_23"]->bin(iW).set(dispersion, dispersion/mean);
432        _e["C3_23"]->bin(iW).set(cq, cq/mean);
433        iq = 4 ;
434        dispersion = 0 ;
435        cq = 0;
436        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
437        _e["D4_23"]->bin(iW).set(dispersion, dispersion/mean);
438        _e["C4_23"]->bin(iW).set(cq, cq/mean);
439      }
440
441      mean = 0.;
442      dispersion = 0.;
443      cq = 0.;
444      R2 = 0;
445      R3 = 0.;
446      K3 = 0.;
447
448      for (auto& histo : _g["mult11_all"]->bins()) {
449        // use scale to get Prob in %, and use counter to count events after cuts
450        size_t iW = histo.index();
451        scale(histo, 100.0/ *_Nevt_after_cuts[iW-1]);
452        //
453        iq = 2 ;
454        _histo_to_moments(histo, iq , mean, dispersion, cq, R2, R3, K3);
455        _e["mean34"]->bin(iW).set(mean, 0.5*dispersion);
456        _e["D2_34"]->bin(iW).set(dispersion, dispersion/mean);
457        _e["C2_34"]->bin(iW).set(cq, cq/mean);
458        _e["R2_34"]->bin(iW).set(R2, R2/mean);
459        _e["R3_34"]->bin(iW).set(R3, R3/mean);
460        _e["K3_34"]->bin(iW).set(K3, K3/mean);
461        iq = 3 ;
462        dispersion = 0 ;
463        cq = 0;
464        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
465        _e["D3_34"]->bin(iW).set(dispersion, dispersion/mean);
466        _e["C3_34"]->bin(iW).set(cq, cq/mean);
467        iq = 4 ;
468        dispersion = 0 ;
469        cq = 0;
470        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
471        _e["D4_34"]->bin(iW).set(dispersion, dispersion/mean);
472        _e["C4_34"]->bin(iW).set(cq, cq/mean);
473      }
474      mean = 0.;
475      dispersion = 0.;
476      cq = 0.;
477      R2 = 0;
478      R3 = 0.;
479      K3 = 0.;
480
481      for (auto& histo : _g["mult12_all"]->bins()) {
482        // use scale to get Prob in %, and use counter to count events after cuts
483        size_t iW = histo.index();
484        scale(histo, 100.0/ *_Nevt_after_cuts[iW-1]);
485        //
486        iq = 2 ;
487        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
488        _e["mean45"]->bin(iW).set(mean, 0.5*dispersion);
489        _e["D2_45"]->bin(iW).set(dispersion, dispersion/mean);
490        _e["C2_45"]->bin(iW).set(cq, cq/mean);
491        _e["R2_45"]->bin(iW).set(R2, R2/mean);
492        _e["R3_45"]->bin(iW).set(R3, R3/mean);
493        _e["K3_45"]->bin(iW).set(K3, K3/mean);
494        iq = 3 ;
495        dispersion = 0 ;
496        cq = 0;
497        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
498        _e["D3_45"]->bin(iW).set(dispersion, dispersion/mean);
499        _e["C3_45"]->bin(iW).set(cq, cq/mean);
500        iq = 4 ;
501        dispersion = 0 ;
502        cq = 0;
503        _histo_to_moments(histo, iq, mean, dispersion, cq, R2, R3, K3);
504        _e["D4_45"]->bin(iW).set(dispersion, dispersion/mean);
505        _e["C4_45"]->bin(iW).set(cq, cq/mean);
506      }
507
508    }
509
510
511    inline void _histo_to_moments(const BinnedHistoPtr<int>& histo_input,
512                                  double iq, double& mean, double& dispersion, double& cq,
513                                  double& R2, double& R3, double& K3) {
514
515      double mysumWX = 0.;
516      double mysumW = 0.;
517
518      if (histo_input->effNumEntries() == 0 || histo_input->sumW() == 0) {
519        mean = dispersion = cq = R2 = R3 = K3 = std::numeric_limits<double>::quiet_NaN();
520        MSG_WARNING("Requested mean of a distribution with no net fill weights");
521      }
522      else {
523        // loop to calculate mean
524        for (const auto& b : histo_input->bins()) { // loop over points
525          mysumWX  += b.sumW()*b.xEdge();
526          mysumW   += b.sumW();
527        }
528        mean = mysumWX/mysumW;
529
530        // loop to calculate dispersion (variance)
531        double var = 0., c = 0., r2 = 0., r3 = 0.;
532        for (const auto& b : histo_input->bins()) { // loop over points
533          const double xval = b.xEdge();
534          const double weight = b.sumW();
535          var = var + weight * pow((xval - mean),iq) ;
536          c = c+pow(xval,iq)*weight;
537          r2 = r2+xval*(xval-1)*weight;
538          r3 = r3+xval*(xval-1)*(xval-2)*weight;
539        }
540        var = var/mysumW;
541        dispersion = pow(var,1./iq);
542        cq = c/(mysumW*pow(mean,iq));
543        R2 = r2/(mysumW*pow(mean,2));
544        R3 = r3/(mysumW*pow(mean,3));
545        K3 = R3 - 3*R2 + 2;
546      }
547    }
548
549    ///@}
550
551
552  private:
553
554    /// @name Histograms
555    ///@{
556
557    map<string,HistoGroupPtr<double,int>> _g;
558    map<string,Estimate1DPtr> _e;
559    CounterPtr _Nevt_after_cuts[4];
560
561    ///@}
562
563  };
564
565  RIVET_DECLARE_PLUGIN(H1_1996_I422230);
566
567}