Rivet analyses referenceH1_1996_I422230Charged particle multiplicities in deep inelastic scattering at HERAExperiment: H1 (HERA) Inspire ID: 422230 Status: VALIDATED Authors:
Beam energies: (27.5, 820.0); (820.0, 27.5) GeV Run details:
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
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}
|