Rivet analyses referenceVENUS_1998_I453613Charged particle multiplicities in heavy and light quark initiated events at 58 GeVExperiment: VENUS (Tristan) Inspire ID: 453613 Status: VALIDATED No authors listed References:
Beam energies: (29.0, 29.0) GeV Run details:
Measurement of the mean charged multiplicities separately for $b\bar b$, and light quark ($uds$) initiated events in $e^+e^-$ interactions at 58 GeV. Source code: VENUS_1998_I453613.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5
6#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
7#include "Rivet/Projections/InitialQuarks.hh"
8
9namespace Rivet {
10
11
12 /// @brief charged multiplicity at 58 GeV
13 class VENUS_1998_I453613 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(VENUS_1998_I453613);
18
19
20 /// @name Analysis methods
21 //@{
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25 declare(Beam(), "Beams");
26 declare(ChargedFinalState(), "CFS");
27 declare(InitialQuarks(), "IQF");
28 book(_cLight , "/TMP/CLIGHT" );
29 book(_cBottom, "/TMP/CBOTTOM");
30 book(_weightLight , "/TMP/WLIGHT" );
31 book(_weightBottom, "/TMP/WBOTTOM");
32 }
33
34
35 /// Perform the per-event analysis
36 void analyze(const Event& event) {
37 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
38 const FinalState& cfs = apply<FinalState>(event, "CFS");
39 if (cfs.size() < 2) vetoEvent;
40
41
42 int flavour = 0;
43 const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
44
45 // If we only have two quarks (qqbar), just take the flavour.
46 // If we have more than two quarks, look for the highest energetic q-qbar pair.
47 if (iqf.particles().size() == 2) {
48 flavour = iqf.particles().front().abspid();
49 }
50 else {
51 map<int, double> quarkmap;
52 for ( const Particle& p : iqf.particles()) {
53 if (quarkmap[p.pid()] < p.E()) {
54 quarkmap[p.pid()] = p.E();
55 }
56 }
57 double maxenergy = 0.;
58 for (int i = 1; i <= 5; ++i) {
59 if (quarkmap[i]+quarkmap[-i] > maxenergy) {
60 flavour = i;
61 }
62 }
63 }
64 const size_t numParticles = cfs.particles().size();
65 switch (flavour) {
66 case 1: case 2: case 3:
67 _weightLight->fill(); ;
68 _cLight->fill(numParticles);
69 break;
70 case 5:
71 _weightBottom->fill();
72 _cBottom->fill(numParticles);
73 break;
74 }
75
76 }
77
78
79 /// Normalise histograms etc., after the run
80 void finalize() {
81 // calculate the averages and diffs
82 if(_weightLight ->effNumEntries()!=0) scale( _cLight, 1./ *_weightLight);
83 if(_weightBottom->effNumEntries()!=0) scale(_cBottom, 1./ *_weightBottom);
84 Counter _cDiff = *_cBottom - *_cLight;
85 // fill the histograms
86 for(unsigned int ix=1;ix<4;++ix) {
87 double val(0.), err(0.0);
88 if(ix==1) {
89 val = _cBottom->val();
90 err = _cBottom->err();
91 }
92 else if(ix==2) {
93 val = _cLight->val();
94 err = _cLight->err();
95 }
96 else if(ix==3) {
97 val = _cDiff.val();
98 err = _cDiff.err();
99 }
100 Scatter2D temphisto(refData(1, 1, ix));
101 Scatter2DPtr mult;
102 book(mult,1, 1, ix);
103 for (size_t b = 0; b < temphisto.numPoints(); b++) {
104 const double x = temphisto.point(b).x();
105 pair<double,double> ex = temphisto.point(b).xErrs();
106 pair<double,double> ex2 = ex;
107 if(ex2.first ==0.) ex2. first=0.0001;
108 if(ex2.second==0.) ex2.second=0.0001;
109 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
110 mult->addPoint(x, val, ex, make_pair(err,err));
111 }
112 else {
113 mult->addPoint(x, 0., ex, make_pair(0.,.0));
114 }
115 }
116 }
117
118 }
119
120 //@}
121
122 /// @name Multiplicities
123 //@{
124 CounterPtr _cLight;
125 CounterPtr _cCharm;
126 CounterPtr _cBottom;
127 //@}
128
129 /// @name Weights
130 //@{
131 CounterPtr _weightLight;
132 CounterPtr _weightBottom;
133 //@}
134
135 };
136
137
138 // The hook for the plugin system
139 RIVET_DECLARE_PLUGIN(VENUS_1998_I453613);
140
141
142}
|