rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_1996_I422230

Charged particle multiplicities in deep inelastic scattering at HERA
Experiment: H1 (HERA)
Inspire ID: 422230
Status: VALIDATED
Authors:
  • Ariadna Leon
  • Hannes Jung
References: Beams: e+ p+, p+ e+
Beam energies: (27.5, 820.0); (820.0, 27.5) GeV
Run details:
  • validation with 3000001 events generated RAPGAP33, 26.7x820 GeV, Q2>2, 0.01 < y <0.95, IPRO=12, pt2cut=9, mixing with O(alphas) process, Nf_QCDC=4, NFLA=5

Using the H1 detector at HERA, charged particle multiplicity distributions in deep inelastic ep scattering have been measured over a large kinematical region. The evolution with W and Q2 of the multiplicity distribution and of the multiplicity moments in pseudorapidity domains of varying size is studied in the current fragmentation region of the hadronic centre-of-mass frame. Data on multiplicity as well as the various moments of the distribution are given in a variety of psuedorapidity regions. The moments are the generalised dispersions (D), the factorial moments (R), the Muller moments (K) and the C factors (C).

Source code: H1_1996_I422230.cc
  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#include "Rivet/Tools/BinnedHistogram.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Charged particle multiplicities in deep inelastic scattering at HERA (H1)
 13  class H1_1996_I422230 : public Analysis {
 14  public:
 15
 16    const vector<double> WEdges {80., 115, 150., 185., 220.};
 17
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1996_I422230);
 21
 22
 23    /// @name Analysis methods
 24    ///@{
 25
 26    /// Book histograms and initialise projections before the run
 27    void init() {
 28
 29      // Initialise and register projections
 30      const DISLepton disl;
 31      declare(disl, "Lepton");
 32      declare(DISKinematics(), "Kinematics");
 33
 34      // The basic final-state projection:
 35      // all final-state particles within
 36      // the given eta acceptance
 37      const FinalState fs;
 38      declare(fs,"FS");
 39      const ChargedFinalState cfs(disl.remainingFinalState());
 40      declare(cfs,"CFS");
 41
 42
 43
 44      //binned histograms to count for the multiplicity
 45      for (size_t ix = 0; ix < 4; ++ix) {
 46        book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts"+ to_string(ix));
 47      }
 48      Histo1DPtr dummy;
 49      for (int iW = 0; iW < 4; ++iW) {
 50        // cout << " iW " << iW << " " << WEdges[iW+1] << " " << WEdges[iW] << endl;
 51        _h_mult1.add(WEdges[iW], WEdges[iW+1],book(dummy,iW+1,1,1));
 52        _h_mult2.add(WEdges[iW], WEdges[iW+1],book(dummy,iW+1,1,2));
 53        _h_mult3.add(WEdges[iW], WEdges[iW+1],book(dummy,iW+1,1,3));
 54        _h_mult4.add(WEdges[iW], WEdges[iW+1],book(dummy,iW+1,1,4));
 55        _h_mult_all.add(WEdges[iW], WEdges[iW+1],book(dummy,"TMP/dummy"+ to_string(iW),refData(4,1,4)));
 56        _h_mult10_all.add(WEdges[iW], WEdges[iW+1],book(dummy,"TMP/dummy1"+ to_string(iW),refData(4,1,4)));
 57        _h_mult11_all.add(WEdges[iW], WEdges[iW+1],book(dummy,"TMP/dummy2"+ to_string(iW),refData(4,1,4)));
 58        _h_mult12_all.add(WEdges[iW], WEdges[iW+1],book(dummy,"TMP/dummy3"+ to_string(iW),refData(4,1,4)));
 59      }
 60      //histograms for the statistical moments and the mean
 61
 62      book(_h_mean0,5,1,1);
 63      book(_h_D2_0,5,1,2);
 64      book(_h_D3_0,5,1,3);
 65      book(_h_D4_0,5,1,4);
 66      book(_h_C2_0,5,1,5);
 67      book(_h_C3_0,5,1,6);
 68      book(_h_C4_0,5,1,7);
 69      book(_h_R2_0,5,1,8);
 70      book(_h_R3_0,5,1,9);
 71
 72      book(_h_mean12,6,1,1);
 73      book(_h_D2_12,6,1,2);
 74      book(_h_D3_12,6,1,3);
 75      book(_h_D4_12,6,1,4);
 76      book(_h_C2_12,6,1,5);
 77      book(_h_C3_12,6,1,6);
 78      book(_h_C4_12,6,1,7);
 79      book(_h_R2_12,6,1,8);
 80      book(_h_R3_12,6,1,9);
 81      book(_h_K3_12,6,1,10);
 82
 83      book(_h_mean13,7,1,1);
 84      book(_h_D2_13,7,1,2);
 85      book(_h_D3_13,7,1,3);
 86      book(_h_D4_13,7,1,4);
 87      book(_h_C2_13,7,1,5);
 88      book(_h_C3_13,7,1,6);
 89      book(_h_C4_13,7,1,7);
 90      book(_h_R2_13,7,1,8);
 91      book(_h_R3_13,7,1,9);
 92      book(_h_K3_13,7,1,10);
 93
 94      book(_h_mean14,8,1,1);
 95      book(_h_D2_14,8,1,2);
 96      book(_h_D3_14,8,1,3);
 97      book(_h_D4_14,8,1,4);
 98      book(_h_C2_14,8,1,5);
 99      book(_h_C3_14,8,1,6);
100      book(_h_C4_14,8,1,7);
101      book(_h_R2_14,8,1,8);
102      book(_h_R3_14,8,1,9);
103
104      book(_h_mean15,9,1,1);
105      book(_h_D2_15,9,1,2);
106      book(_h_D3_15,9,1,3);
107      book(_h_D4_15,9,1,4);
108      book(_h_C2_15,9,1,5);
109      book(_h_C3_15,9,1,6);
110      book(_h_C4_15,9,1,7);
111      book(_h_R2_15,9,1,8);
112      book(_h_R3_15,9,1,9);
113
114      book(_h_mean23,10,1,1);
115      book(_h_D2_23,10,1,2);
116      book(_h_D3_23,10,1,3);
117      book(_h_D4_23,10,1,4);
118      book(_h_C2_23,10,1,5);
119      book(_h_C3_23,10,1,6);
120      book(_h_C4_23,10,1,7);
121      book(_h_R2_23,10,1,8);
122      book(_h_R3_23,10,1,9);
123      book(_h_K3_23,10,1,10);
124
125      book(_h_mean34,11,1,1);
126      book(_h_D2_34,11,1,2);
127      book(_h_D3_34,11,1,3);
128      book(_h_D4_34,11,1,4);
129      book(_h_C2_34,11,1,5);
130      book(_h_C3_34,11,1,6);
131      book(_h_C4_34,11,1,7);
132      book(_h_R2_34,11,1,8);
133      book(_h_R3_34,11,1,9);
134      book(_h_K3_34,11,1,10);
135
136      book(_h_mean45,12,1,1);
137      book(_h_D2_45,12,1,2);
138      book(_h_D3_45,12,1,3);
139      book(_h_D4_45,12,1,4);
140      book(_h_C2_45,12,1,5);
141      book(_h_C3_45,12,1,6);
142      book(_h_C4_45,12,1,7);
143      book(_h_R2_45,12,1,8);
144      book(_h_R3_45,12,1,9);
145      book(_h_K3_45,12,1,10);
146    }
147
148
149    /// Perform the per-event analysis
150    void analyze(const Event& event) {
151
152      //apply the final state of all particles and the one for only charged particles
153      const FinalState& fs = apply<FinalState>(event, "FS");
154      const size_t numParticles = fs.particles().size();
155      const ChargedFinalState& cfs = apply<ChargedFinalState>(event,"CFS");
156      const Particles& particles = cfs.particles();
157
158      //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
159
160      if (numParticles<2) {
161        MSG_DEBUG("Failed leptonic cut");
162        vetoEvent;
163      }
164
165      //apply DIS kinematics in the event
166      const DISKinematics& dk = apply<DISKinematics>(event,"Kinematics");
167      //const DISLepton& dl = apply<DISLepton>(event,"Lepton");
168
169      //get the DIS Kinematics but not in a loop because they are variables that descrube the type of event not the particles
170      double Q2 = dk.Q2();
171      double W2 = dk.W2();
172      double W = std::sqrt(W2);
173      bool cut1 = W<220. && W>80.;
174      if (!cut1) vetoEvent;
175      bool cut = Q2<1000. && Q2>10. && W>80. && W<220.;
176      if (!cut) vetoEvent;
177
178      double Efwd = 0. ;
179      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
180        const Particle& p = particles[ip1];
181        double theta = p.theta()/degree ;
182        if ( inRange(theta,4.4,15.) ) {
183          Efwd = Efwd + p.E() ;
184        }
185      }
186
187      bool cut_fwd = Efwd > 0.5 && dk.beamLepton().E() > 12. ;
188      if (!cut_fwd) vetoEvent ;
189
190      for (int iW = 0; iW < 4; ++iW) {
191        if (inRange(W, WEdges[iW],WEdges[iW+1])) {
192          _Nevt_after_cuts[iW] -> fill();
193        }
194      }
195
196
197      //boost to the hadronic centre of mass frame
198      const LorentzTransform hcmboost = dk.boostHCM();
199      double kall = 0.;
200      double k1 = 0.;
201      double k2 = 0.;
202      double k3 = 0.;
203      double k4 = 0.;
204      double k10 = 0.;
205      double k11 = 0.;
206      double k12 = 0.;
207      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
208        const Particle& p = particles[ip1];
209        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
210        const double etahcm_char = hcmMom.eta();
211        if (etahcm_char>0.) {
212          kall = kall+1;
213        }
214        if (etahcm_char>1. && etahcm_char<2.) {
215          k1 = k1+1;
216        }
217        if (etahcm_char>1. && etahcm_char<3.) {
218          k2 = k2+1;
219        }
220        if (etahcm_char>1. && etahcm_char<4.) {
221          k3 = k3+1;
222        }
223        if (etahcm_char>1. && etahcm_char<5.) {
224          k4 = k4+1;
225        }
226        if (etahcm_char>2. && etahcm_char<3.) {
227          k10 = k10+1;
228        }
229        if (etahcm_char>3. && etahcm_char<4.) {
230          k11 = k11+1;
231        }
232        if (etahcm_char>4. && etahcm_char<5.) {
233          k12 = k12+1;
234        }
235
236      }
237      // cout<<k1<<endl;
238      _h_mult_all.fill(W,kall);
239      _h_mult10_all.fill(W,k10);
240      _h_mult11_all.fill(W,k11);
241      _h_mult12_all.fill(W,k12);
242      _h_mult1.fill(W,k1);
243      _h_mult2.fill(W,k2);
244      _h_mult3.fill(W,k3);
245      _h_mult4.fill(W,k4);
246
247    }
248
249    /// Normalise histograms etc., after the run
250    void finalize() {
251      // cout << _h_mult1.histo() << endl;
252      int iW = 0 ;
253      double iq  ;
254      double mean, dispersion, cq ,R2,R3,K3 ;
255      for (Histo1DPtr histo : _h_mult_all.histos()) {
256        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
257        //cout << " iW " << iW << " Edges " << WEdges[iW+1] << " " << WEdges[iW] << " Werr = " << Werr << endl;
258        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
259        //cout << " new histo: W = "  << Wmid << " " << histo->name() << endl;
260        iq = 2 ;
261        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
262        //cout << " mean " << mean << " dispersion " << dispersion << " for iq = " << iq << endl;
263
264        // just to have some values, needs to be corrected
265
266        _h_mean0->addPoint(Wmid, mean, Werr, dispersion/2);
267        _h_D2_0->addPoint(Wmid, dispersion, Werr, dispersion/mean);
268        _h_C2_0->addPoint(Wmid, cq, Werr, cq/mean);
269        _h_R2_0->addPoint(Wmid, R2, Werr, R2/mean);
270        _h_R3_0->addPoint(Wmid,R3,Werr,R3/mean);
271        iq = 3 ;
272        dispersion = 0 ;
273        cq = 0;
274        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
275        _h_D3_0->addPoint(Wmid, dispersion, Werr, dispersion/mean);
276        _h_C3_0->addPoint(Wmid, cq, Werr, cq/mean);
277        iq = 4 ;
278        dispersion = 0 ;
279        cq = 0;
280        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
281        _h_D4_0->addPoint(Wmid, dispersion, Werr, dispersion/mean);
282        _h_C4_0->addPoint(Wmid, cq, Werr, cq/mean);
283
284        ++iW;
285      }
286      iW = 0 ;
287      mean = 0.;
288      dispersion = 0.;
289      cq = 0.;
290      R2 = 0;
291      R3 = 0.;
292      K3 = 0.;
293      for (Histo1DPtr histo : _h_mult1.histos()) {
294        // use scale to get Prob in %, and use counter to count events after cuts
295        //cout << " Nevt " << dbl(*_Nevt_after_cuts[iW]) << endl;
296        scale(histo, 100.0/ *_Nevt_after_cuts[iW]);
297        //
298        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
299        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
300        iq = 2 ;
301        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
302
303        _h_mean12->addPoint(Wmid, mean, Werr, dispersion/2);
304        _h_D2_12->addPoint(Wmid, dispersion, Werr, dispersion/mean);
305        _h_C2_12->addPoint(Wmid, cq, Werr, cq/mean);
306        _h_R2_12->addPoint(Wmid, R2, Werr, R2/mean);
307        _h_R3_12->addPoint(Wmid,R3,Werr,R3/mean);
308        _h_K3_12->addPoint(Wmid,K3,Werr,K3/mean);
309        iq = 3 ;
310        dispersion = 0 ;
311        cq = 0;
312        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
313        _h_D3_12->addPoint(Wmid, dispersion, Werr, dispersion/mean);
314        _h_C3_12->addPoint(Wmid, cq, Werr, cq/mean);
315        iq = 4 ;
316        dispersion = 0 ;
317        cq = 0;
318        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
319        _h_D4_12->addPoint(Wmid, dispersion, Werr, dispersion/mean);
320        _h_C4_12->addPoint(Wmid, cq, Werr, cq/mean);
321        ++iW;
322      }
323      iW = 0 ;
324      mean = 0.;
325      dispersion = 0.;
326      cq = 0.;
327      R2 = 0;
328      R3 = 0.;
329      K3 = 0.;
330
331      for (Histo1DPtr histo : _h_mult2.histos()) {
332        // use scale to get Prob in %, and use counter to count events after cuts
333        //cout << " Nevt " << dbl(*_Nevt_after_cuts[0]) << endl;
334        scale(histo, 100.0/ *_Nevt_after_cuts[iW]);
335        //
336        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
337        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
338        iq = 2 ;
339        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
340
341        _h_mean13->addPoint(Wmid, mean, Werr, dispersion/2);
342        _h_D2_13->addPoint(Wmid, dispersion, Werr, dispersion/mean);
343        _h_C2_13->addPoint(Wmid, cq, Werr, cq/mean);
344        _h_R2_13->addPoint(Wmid, R2, Werr, R2/mean);
345        _h_R3_13->addPoint(Wmid,R3,Werr,R3/mean);
346        _h_K3_13->addPoint(Wmid,K3,Werr,K3/mean);
347        iq = 3 ;
348        dispersion = 0 ;
349        cq = 0;
350        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
351        _h_D3_13->addPoint(Wmid, dispersion, Werr, dispersion/mean);
352        _h_C3_13->addPoint(Wmid, cq, Werr, cq/mean);
353        iq = 4 ;
354        dispersion = 0 ;
355        cq = 0;
356        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
357        _h_D4_13->addPoint(Wmid, dispersion, Werr, dispersion/mean);
358        _h_C4_13->addPoint(Wmid, cq, Werr, cq/mean);
359        ++iW;
360      }
361      iW = 0 ;
362      mean = 0.;
363      dispersion = 0.;
364      cq = 0.;
365      R2 = 0;
366      R3 = 0.;
367      K3 = 0.;
368
369      for (Histo1DPtr histo : _h_mult3.histos()) {
370        // use scale to get Prob in %, and use counter to count events after cuts
371        //cout << " Nevt " << dbl(*_Nevt_after_cuts[0]) << endl;
372        scale(histo, 100.0/ *_Nevt_after_cuts[iW]);
373        //
374        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
375        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
376        iq = 2 ;
377        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
378        _h_mean14->addPoint(Wmid, mean, Werr, dispersion/2);
379        _h_D2_14->addPoint(Wmid, dispersion, Werr, dispersion/mean);
380        _h_C2_14->addPoint(Wmid, cq, Werr, cq/mean);
381        _h_R2_14->addPoint(Wmid, R2, Werr, R2/mean);
382        _h_R3_14->addPoint(Wmid,R3,Werr,R3/mean);
383        iq = 3 ;
384        dispersion = 0 ;
385        cq = 0;
386        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
387        _h_D3_14->addPoint(Wmid, dispersion, Werr, dispersion/mean);
388        _h_C3_14->addPoint(Wmid, cq, Werr, cq/mean);
389        iq = 4 ;
390        dispersion = 0 ;
391        cq = 0;
392        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
393        _h_D4_14->addPoint(Wmid, dispersion, Werr, dispersion/mean);
394        _h_C4_14->addPoint(Wmid, cq, Werr, cq/mean);
395        ++iW;
396      }
397      iW = 0 ;
398      mean = 0.;
399      dispersion = 0.;
400      cq = 0.;
401      R2 = 0;
402      R3 = 0.;
403      K3 = 0.;
404
405      for (Histo1DPtr histo : _h_mult4.histos()) {
406        // use scale to get Prob in %, and use counter to count events after cuts
407        //cout << " Nevt " << dbl(*_Nevt_after_cuts[0]) << endl;
408        scale(histo, 100.0/ *_Nevt_after_cuts[iW]);
409        //
410        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
411        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
412        iq = 2 ;
413        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
414        _h_mean15->addPoint(Wmid, mean, Werr, dispersion/2);
415        _h_D2_15->addPoint(Wmid, dispersion, Werr, dispersion/mean);
416        _h_C2_15->addPoint(Wmid, cq, Werr, cq/mean);
417        _h_R2_15->addPoint(Wmid, R2, Werr, R2/mean);
418        _h_R3_15->addPoint(Wmid,R3,Werr,R3/mean);
419        iq = 3 ;
420        dispersion = 0 ;
421        cq = 0;
422        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
423        _h_D3_15->addPoint(Wmid, dispersion, Werr, dispersion/mean);
424        _h_C3_15->addPoint(Wmid, cq, Werr, cq/mean);
425        iq = 4 ;
426        dispersion = 0 ;
427        cq = 0;
428        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
429        _h_D4_15->addPoint(Wmid, dispersion, Werr, dispersion/mean);
430        _h_C4_15->addPoint(Wmid, cq, Werr, cq/mean);
431        ++iW;
432      }
433      iW = 0 ;
434      mean = 0.;
435      dispersion = 0.;
436      cq = 0.;
437      R2 = 0;
438      R3 = 0.;
439      K3 = 0.;
440
441      for (Histo1DPtr histo : _h_mult10_all.histos()) {
442        // use scale to get Prob in %, and use counter to count events after cuts
443        //cout << " Nevt " << dbl(*_Nevt_after_cuts[0]) << endl;
444        scale(histo, 100.0/ *_Nevt_after_cuts[iW]);
445        //
446        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
447        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
448        iq = 2 ;
449        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
450        _h_mean23->addPoint(Wmid, mean, Werr, dispersion/2);
451        _h_D2_23->addPoint(Wmid, dispersion, Werr, dispersion/mean);
452        _h_C2_23->addPoint(Wmid, cq, Werr, cq/mean);
453        _h_R2_23->addPoint(Wmid, R2, Werr,R2/mean);
454        _h_R3_23->addPoint(Wmid,R3,Werr,R3/mean);
455        _h_K3_23->addPoint(Wmid,K3,Werr,K3/mean);
456        iq = 3 ;
457        dispersion = 0 ;
458        cq = 0;
459        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
460        _h_D3_23->addPoint(Wmid, dispersion, Werr, dispersion/mean);
461        _h_C3_23->addPoint(Wmid, cq, Werr, cq/mean);
462        iq = 4 ;
463        dispersion = 0 ;
464        cq = 0;
465        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
466        _h_D4_23->addPoint(Wmid, dispersion, Werr, dispersion/mean);
467        _h_C4_23->addPoint(Wmid, cq, Werr, cq/mean);
468        ++iW;
469      }
470
471      iW = 0 ;
472      mean = 0.;
473      dispersion = 0.;
474      cq = 0.;
475      R2 = 0;
476      R3 = 0.;
477      K3 = 0.;
478
479      for (Histo1DPtr histo : _h_mult11_all.histos()) {
480        // use scale to get Prob in %, and use counter to count events after cuts
481        //cout << " Nevt " << dbl(*_Nevt_after_cuts[0]) << endl;
482        scale(histo, 100.0/ *_Nevt_after_cuts[iW]);
483        //
484        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
485        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
486        iq = 2 ;
487        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
488        _h_mean34->addPoint(Wmid, mean, Werr, dispersion/2);
489        _h_D2_34->addPoint(Wmid, dispersion, Werr, dispersion/mean);
490        _h_C2_34->addPoint(Wmid, cq, Werr, cq/mean);
491        _h_R2_34->addPoint(Wmid, R2, Werr, R2/mean);
492        _h_R3_34->addPoint(Wmid,R3,Werr,R3/mean);
493        _h_K3_34->addPoint(Wmid,K3,Werr,K3/mean);
494        iq = 3 ;
495        dispersion = 0 ;
496        cq = 0;
497        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
498        _h_D3_34->addPoint(Wmid, dispersion, Werr, dispersion/mean);
499        _h_C3_34->addPoint(Wmid, cq, Werr, cq/mean);
500        iq = 4 ;
501        dispersion = 0 ;
502        cq = 0;
503        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
504        _h_D4_34->addPoint(Wmid, dispersion, Werr, dispersion/mean);
505        _h_C4_34->addPoint(Wmid, cq, Werr, cq/mean);
506        ++iW;
507      }
508      iW = 0 ;
509      mean = 0.;
510      dispersion = 0.;
511      cq = 0.;
512      R2 = 0;
513      R3 = 0.;
514      K3 = 0.;
515
516      for (Histo1DPtr histo : _h_mult12_all.histos()) {
517        // use scale to get Prob in %, and use counter to count events after cuts
518        //cout << " Nevt " << dbl(*_Nevt_after_cuts[0]) << endl;
519        scale(histo, 100.0/ *_Nevt_after_cuts[iW]);
520        //
521        double Werr = (WEdges[iW+1] - WEdges[iW])/2. ;
522        double Wmid = (WEdges[iW+1] + WEdges[iW])/2. ;
523        iq = 2 ;
524        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
525        _h_mean45->addPoint(Wmid, mean, Werr, dispersion/2);
526        _h_D2_45->addPoint(Wmid, dispersion, Werr, dispersion/mean);
527        _h_C2_45->addPoint(Wmid, cq, Werr, cq/mean);
528        _h_R2_45->addPoint(Wmid, R2, Werr, R2/mean);
529        _h_R3_45->addPoint(Wmid,R3,Werr,R3/mean);
530        // cout<<"R2= "<<R2<<"R3= "<<R3<<"K3= "<<K3<<endl;
531        _h_K3_45->addPoint(Wmid,K3,Werr,K3/mean);
532        iq = 3 ;
533        dispersion = 0 ;
534        cq = 0;
535        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
536        _h_D3_45->addPoint(Wmid, dispersion, Werr, dispersion/mean);
537        _h_C3_45->addPoint(Wmid, cq, Werr, cq/mean);
538        iq = 4 ;
539        dispersion = 0 ;
540        cq = 0;
541        _histo_to_moments(histo, iq , mean, dispersion,cq,R2,R3,K3) ;
542        _h_D4_45->addPoint(Wmid, dispersion, Werr, dispersion/mean);
543        _h_C4_45->addPoint(Wmid, cq, Werr, cq/mean);
544        ++iW;
545      }
546
547    }
548
549
550    inline void _histo_to_moments(Histo1DPtr histo_input, double iq, double & mean, double & dispersion, double & cq,double &R2,double &R3, double &K3) {
551
552      //cout << " histo mean = " << histo_input->xMean() << " variance " << histo_input->xVariance() << endl;
553      double mysumWX = 0. ;
554      // double mysumW2X = 0. ;
555      // double mysumWX2 = 0. ;
556      // double mysumW2 = 0. ;
557      double mysumW = 0. ;
558
559      // cout << histo_input->name() << endl;
560      if (histo_input->effNumEntries() == 0 || histo_input->sumW() == 0) {
561        MSG_WARNING("Requested mean of a distribution with no net fill weights");
562      } else {
563        // loop to calcualte mean
564        for (size_t b = 0; b < histo_input->numBins(); ++b) { // loop over points
565          mysumWX  += histo_input->bin(b).height()      *  histo_input->bin(b).xMid() ;
566          //mysumW2X += sqr(histo_input->bin(b).height()) *  histo_input->bin(b).xMid() ;
567          //mysumWX2 += histo_input->bin(b).height()      *  sqr(histo_input->bin(b).xMid()) ;
568          //mysumW2  += sqr(histo_input->bin(b).height()) ;
569          mysumW   += histo_input->bin(b).height() ;
570        }
571        mean = mysumWX/mysumW ;
572
573        // loop to calculate dispersion (variance)
574        double var = 0.;
575        for (size_t b = 0; b < histo_input->numBins(); ++b) { // loop over points
576          double xval = histo_input->bin(b).xMid() ;
577          double weight = histo_input->bin(b).height() ;
578          var = var + weight * pow((xval - mean),iq) ;
579        }
580
581        var = var/mysumW ;
582
583        dispersion = pow(var,1./iq) ;
584        double c = 0.;
585        for (size_t b = 0; b < histo_input->numBins(); ++b) { // loop over points
586          double xval = histo_input->bin(b).xMid() ;
587          double weight = histo_input->bin(b).height() ;
588          c = c+pow(xval,iq)*weight;
589        }
590        cq = c/(mysumW*pow(mean,iq));
591        double r2 = 0.;
592        double r3 = 0.;
593        for (size_t b = 0; b < histo_input->numBins(); ++b) { // loop over points
594          double xval = histo_input->bin(b).xMid() ;
595          double weight = histo_input->bin(b).height() ;
596          r2 = r2+xval*(xval-1)*weight;
597          r3 = r3+xval*(xval-1)*(xval-2)*weight;
598        }
599        R2 = r2/(mysumW*pow(mean,2));
600        R3 = r3/(mysumW*pow(mean,3));
601        K3 = R3 - 3*R2 + 2;
602      }
603    }
604
605    ///@}
606
607
608  private:
609
610    /// @name Histograms
611    ///@{
612
613    BinnedHistogram _h_mult1;
614    BinnedHistogram _h_mult2;
615    BinnedHistogram _h_mult3;
616    BinnedHistogram _h_mult4;
617    BinnedHistogram _h_mult_all;
618    BinnedHistogram _h_mult10_all;
619    BinnedHistogram _h_mult11_all;
620    BinnedHistogram _h_mult12_all;
621
622    Histo1DPtr _h2_W5;
623    Histo1DPtr _h2_W6;
624    Histo1DPtr _h2_W7;
625    Histo1DPtr _h2_W8;
626    Histo1DPtr _h2_W9;
627    Histo1DPtr _h2_W10;
628    Histo1DPtr _h2_W11;
629    Histo1DPtr _h2_W12;
630    //Histo1DPtr _h2_W5_R2;
631    Scatter2DPtr _h_mean0;
632    Scatter2DPtr _h_D2_0;
633    Scatter2DPtr _h_D3_0;
634    Scatter2DPtr _h_D4_0;
635    Scatter2DPtr _h_C2_0;
636    Scatter2DPtr _h_C3_0;
637    Scatter2DPtr _h_C4_0;
638    Scatter2DPtr _h_R2_0;
639    Scatter2DPtr _h_R3_0;
640    Scatter2DPtr _h_mean12;
641    Scatter2DPtr _h_D2_12;
642    Scatter2DPtr _h_D3_12;
643    Scatter2DPtr _h_D4_12;
644    Scatter2DPtr _h_C2_12;
645    Scatter2DPtr _h_C3_12;
646    Scatter2DPtr _h_C4_12;
647    Scatter2DPtr _h_R2_12;
648    Scatter2DPtr _h_R3_12;
649    Scatter2DPtr _h_K3_12;
650    Scatter2DPtr _h_mean13;
651    Scatter2DPtr _h_D2_13;
652    Scatter2DPtr _h_D3_13;
653    Scatter2DPtr _h_D4_13;
654    Scatter2DPtr _h_C2_13;
655    Scatter2DPtr _h_C3_13;
656    Scatter2DPtr _h_C4_13;
657    Scatter2DPtr _h_R2_13;
658    Scatter2DPtr _h_R3_13;
659    Scatter2DPtr _h_K3_13;
660    Scatter2DPtr _h_mean14;
661    Scatter2DPtr _h_D2_14;
662    Scatter2DPtr _h_D3_14;
663    Scatter2DPtr _h_D4_14;
664    Scatter2DPtr _h_C2_14;
665    Scatter2DPtr _h_C3_14;
666    Scatter2DPtr _h_C4_14;
667    Scatter2DPtr _h_R2_14;
668    Scatter2DPtr _h_R3_14;
669    Scatter2DPtr _h_mean15;
670    Scatter2DPtr _h_D2_15;
671    Scatter2DPtr _h_D3_15;
672    Scatter2DPtr _h_D4_15;
673    Scatter2DPtr _h_C2_15;
674    Scatter2DPtr _h_C3_15;
675    Scatter2DPtr _h_C4_15;
676    Scatter2DPtr _h_R2_15;
677    Scatter2DPtr _h_R3_15;
678    Scatter2DPtr _h_mean23;
679    Scatter2DPtr _h_D2_23;
680    Scatter2DPtr _h_D3_23;
681    Scatter2DPtr _h_D4_23;
682    Scatter2DPtr _h_C2_23;
683    Scatter2DPtr _h_C3_23;
684    Scatter2DPtr _h_C4_23;
685    Scatter2DPtr _h_R2_23;
686    Scatter2DPtr _h_R3_23;
687    Scatter2DPtr _h_K3_23;
688    Scatter2DPtr _h_mean34;
689    Scatter2DPtr _h_D2_34;
690    Scatter2DPtr _h_D3_34;
691    Scatter2DPtr _h_D4_34;
692    Scatter2DPtr _h_C2_34;
693    Scatter2DPtr _h_C3_34;
694    Scatter2DPtr _h_C4_34;
695    Scatter2DPtr _h_R2_34;
696    Scatter2DPtr _h_R3_34;
697    Scatter2DPtr _h_K3_34;
698    Scatter2DPtr _h_mean45;
699    Scatter2DPtr _h_D2_45;
700    Scatter2DPtr _h_D3_45;
701    Scatter2DPtr _h_D4_45;
702    Scatter2DPtr _h_C2_45;
703    Scatter2DPtr _h_C3_45;
704    Scatter2DPtr _h_C4_45;
705    Scatter2DPtr _h_R2_45;
706    Scatter2DPtr _h_R3_45;
707    Scatter2DPtr _h_K3_45;
708
709    CounterPtr _Nevt_after_cuts[4];
710
711    ///@}
712
713  };
714
715
716  RIVET_DECLARE_PLUGIN(H1_1996_I422230);
717
718}