Rivet analyses referenceCLEO_1985_I205668Identified Particle Spectra and rates in $\Upsilon(1S)$ decays and continuum at 10.49 GeVExperiment: CLEO (CESR) Inspire ID: 205668 Status: VALIDATED Authors:
Beam energies: (4.7, 4.7); (5.2, 5.2) GeV Run details:
Spectra and rates for $\pi^\pm$, $K^\pm$, $\pi^0$, $K^0$, $\Lambda$, $\Xi^-$, $\rho^0$, $K^{*\pm}$, $K^{*0}$ and $\phi$ production in $\Upsilon(1S)$ decays and continuum at 10.49 GeV. Source code: CLEO_1985_I205668.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief Spectra in Upsilon(1S) decay and nearby continuum
10 class CLEO_1985_I205668 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1985_I205668);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22 // projections
23 declare(FinalState(), "FS");
24 declare(UnstableParticles(), "UFS");
25 // histos
26 book(_weightSum_cont,"TMP/weightSumcont");
27 book(_weightSum_Ups1,"TMP/weightSumUps1");
28 // multiplcities
29 for (size_t ix=0; ix<2; ++ix) {
30 for (size_t iy=0; iy<12; ++iy) {
31 book(_mult[ix][iy],"/TMP/MULT_" +toString(ix) + "_" +toString(iy));
32 }
33 }
34 // cont spectra
35 book(_cont["pip"] , 1,1,1);
36 book(_cont["Kp"] , 2,1,1);
37 book(_cont["p"] , 3,1,1);
38 book(_cont["pi0"] , 4,1,1);
39 book(_cont["K0"] , 5,1,1);
40 book(_cont["lam"] , 6,1,1);
41 book(_cont["xi"] , 7,1,1);
42 book(_cont["rho"] , 8,1,1);
43 book(_cont["Kstarp"], 9,1,1);
44 book(_cont["Kstar0"],10,1,1);
45 book(_cont["phi"] ,11,1,1);
46 // ups spectra
47 book(_ups1["pip"] , 1,1,2);
48 book(_ups1["Kp"] , 2,1,2);
49 book(_ups1["p"] , 3,1,2);
50 book(_ups1["pi0"] , 4,1,2);
51 book(_ups1["K0"] , 5,1,2);
52 book(_ups1["lam"] , 6,1,2);
53 book(_ups1["xi"] , 7,1,2);
54 book(_ups1["rho"] , 8,1,2);
55 book(_ups1["Kstarp"], 9,1,2);
56 book(_ups1["Kstar0"],10,1,2);
57 book(_ups1["phi"] ,11,1,2);
58
59 _axes[0]["pip"] = YODA::Axis<double>({0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17,
60 0.19, 0.58, 0.68, 0.78, 0.98});
61 _axes[0]["Kp"] = YODA::Axis<double>({0.03, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19});
62 _axes[0]["p"] = YODA::Axis<double>({0.06, 0.14, 0.155, 0.185, 0.215, 0.245, 0.275});
63 _axes[0]["pi0"] = YODA::Axis<double>({0.1, 0.2, 0.3, 0.4, 0.5});
64 _axes[0]["K0"] = YODA::Axis<double>({0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
65 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9});
66 _axes[0]["lam"] = YODA::Axis<double>({0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4,
67 0.5, 0.65, 0.8, 0.95});
68 _axes[0]["xi"] = YODA::Axis<double>({0.2, 0.3, 0.4, 0.5, 0.6, 0.7});
69 _axes[0]["rho"] = YODA::Axis<double>({0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0});
70 _axes[0]["Kstarp"] = YODA::Axis<double>({0.06, 0.12, 0.24, 0.36, 0.48, 0.6});
71 _axes[0]["Kstar0"] = YODA::Axis<double>({0.0, 0.06, 0.12, 0.24, 0.36, 0.48});
72 _axes[0]["phi"] = YODA::Axis<double>({0.195, 0.385, 0.575, 0.945});
73
74 _axes[1]["pip"] = YODA::Axis<double>({0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19, 0.58, 0.68, 0.88});
75 _axes[1]["Kp"] = YODA::Axis<double>({0.02, 0.1, 0.11, 0.13, 0.15, 0.17, 0.19});
76 _axes[1]["p"] = _axes[0]["p"];
77 _axes[1]["pi0"] = _axes[0]["p0"];
78 _axes[1]["K0"] = _axes[0]["K0"];
79 _axes[1]["lam"] = _axes[0]["lam"];
80 _axes[1]["xi"] = _axes[0]["xi"];
81 _axes[1]["rho"] = YODA::Axis<double>({0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7});
82 _axes[1]["Kstarp"] = _axes[0]["Kstarp"];
83 _axes[1]["Kstar0"] = YODA::Axis<double>({0.06, 0.12, 0.24, 0.36, 0.48});
84 _axes[1]["phi"] = YODA::Axis<double>({0.28, 0.36, 0.7, 1.0});
85
86 }
87
88 /// Recursively walk the decay tree to find decay products of @a p
89 void findDecayProducts(Particle mother, Particles& unstable) {
90 for(const Particle & p: mother.children()) {
91 const int id = p.abspid();
92 if (id == PID::PIPLUS || id==PID::KPLUS || id==PID::PROTON ||
93 id==PID::PI0 || id==PID::K0S || id==PID::K0L ||
94 id==PID::LAMBDA || id==PID::XIMINUS || id==PID::RHO0 ||
95 id==323 || id==313 || id==225 || id==PID::PHI) {
96 unstable.push_back(p);
97 }
98 if(!p.children().empty())
99 findDecayProducts(p, unstable);
100 }
101 }
102
103 /// Perform the per-event analysis
104 void analyze(const Event& event) {
105 if (_edges[0].empty()) {
106 for (const auto& item : _cont) {
107 _edges[0][item.first] = item.second->xEdges();
108 _edges[1][item.first] = _ups1[item.first]->xEdges();
109 }
110 }
111 // Find the upsilons
112 // First in unstable final state
113 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
114 Particles upsilons = ufs.particles(Cuts::pid==553);
115 // continuum
116 if (upsilons.empty()) {
117 _weightSum_cont->fill();
118 const FinalState& fs = apply<FinalState>(event, "FS");
119 // FS particles
120 for (const Particle& p : fs.particles()) {
121 int id = p.abspid();
122 double xp = 2.*p.p3().mod()/sqrtS();
123 if(id==PID::PIPLUS) {
124 discfill("pip", xp, 0);
125 _mult[1][0]->fill();
126 }
127 else if(id==PID::KPLUS) {
128 discfill("Kp", xp, 0);
129 _mult[1][1]->fill();
130 }
131 else if(id==PID::PROTON) {
132 discfill("p", xp, 0);
133 _mult[1][2]->fill();
134 }
135 }
136 // Unstable particles
137 for (const Particle& p : ufs.particles()) {
138 int id = p.abspid();
139 double xp = 2.*p.p3().mod()/sqrtS();
140 if(id==PID::PI0) {
141 discfill("pi0", xp, 0);
142 _mult[1][3]->fill();
143 }
144 else if(id==PID::K0S || id==PID::K0L) {
145 discfill("K0", xp, 0);
146 _mult[1][4]->fill();
147 }
148 else if(id==PID::LAMBDA) {
149 discfill("lam", xp, 0);
150 _mult[1][5]->fill();
151 }
152 else if(id==PID::XIMINUS) {
153 discfill("xi", xp, 0);
154 _mult[1][6]->fill();
155 }
156 else if(id==PID::RHO0) {
157 discfill("rho", xp, 0);
158 _mult[1][7]->fill();
159 }
160 else if(id==323) {
161 discfill("Kstarp", xp, 0);
162 _mult[1][8]->fill();
163 }
164 else if(id==313) {
165 discfill("Kstar0", xp, 0);
166 _mult[1][9]->fill();
167 }
168 else if(id==PID::PHI) {
169 discfill("phi", xp, 0);
170 _mult[1][10]->fill();
171 }
172 else if(id==225) {
173 _mult[1][11]->fill();
174 }
175 }
176 }
177 else {
178 for (const Particle& ups : upsilons) {
179 _weightSum_Ups1->fill();
180 Particles unstable;
181 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
182 // Find the decay products we want
183 findDecayProducts(ups,unstable);
184 for (const Particle& p : unstable) {
185 int id = p.abspid();
186 double xp = 2.*boost.transform(p.momentum()).p3().mod()/ups.mass();
187 if(id==PID::PIPLUS) {
188 discfill("pip", xp, 1);
189 _mult[0][0]->fill();
190 }
191 else if(id==PID::KPLUS) {
192 discfill("Kp", xp, 1);
193 _mult[0][1]->fill();
194 }
195 else if(id==PID::PROTON) {
196 discfill("p", xp, 1);
197 _mult[0][2]->fill();
198 }
199 else if(id==PID::PI0) {
200 discfill("pi0", xp, 1);
201 _mult[0][3]->fill();
202 }
203 else if(id==PID::K0S || id==PID::K0L) {
204 discfill("K0", xp, 1);
205 _mult[0][4]->fill();
206 }
207 else if(id==PID::LAMBDA) {
208 discfill("lam", xp, 1);
209 _mult[0][5]->fill();
210 }
211 else if(id==PID::XIMINUS) {
212 discfill("xi", xp, 1);
213 _mult[0][6]->fill();
214 }
215 else if(id==PID::RHO0) {
216 discfill("rho", xp, 1);
217 _mult[0][7]->fill();
218 }
219 else if(id==323) {
220 discfill("Kstarp", xp, 1);
221 _mult[0][8]->fill();
222 }
223 else if(id==313) {
224 discfill("Kstar0", xp, 1);
225 _mult[0][9]->fill();
226 }
227 else if(id==PID::PHI) {
228 discfill("phi", xp, 1);
229 _mult[0][10]->fill();
230 }
231 else if(id==225) {
232 _mult[0][11]->fill();
233 }
234 }
235 }
236 }
237 }
238
239 void discfill(const string& name, const double value, const size_t k) {
240 string edge = "OTHER";
241 const size_t idx = _axes[k][name].index(value);
242 if (idx && idx <= _edges[k][name].size()) edge = _edges[k][name][idx-1];
243 (k? _ups1 : _cont)[name]->fill(edge, value);
244 }
245
246
247 /// Normalise histograms etc., after the run
248 void finalize() {
249 // multiplicities
250 const vector<CounterPtr> scales = {_weightSum_Ups1, _weightSum_cont};
251 for (size_t ix=0; ix<12; ++ix) {
252 BinnedEstimatePtr<string> est;
253 book(est, ix+12, 1, 1);
254 for (size_t iy=0; iy<2; ++iy) {
255 if (scales[iy]->val() > 0.) {
256 scale(_mult[iy][ix], 1./ *scales[iy]);
257 est->bin(iy+1).set(_mult[iy][ix]->val(), _mult[iy][ix]->err());
258 }
259 }
260 }
261 // spectra
262 if (_weightSum_cont->val() > 0.) {
263 scale(_cont, 1. / *_weightSum_cont);
264 }
265 if (_weightSum_Ups1->val() > 0.) {
266 scale(_ups1, 1. / *_weightSum_Ups1);
267 }
268 }
269
270 /// @}
271
272
273 /// @name Histograms
274 /// @{
275 map<string,BinnedHistoPtr<string>> _cont, _ups1;
276 map<string, YODA::Axis<double>> _axes[2];
277 map<string, vector<string>> _edges[2];
278 CounterPtr _weightSum_cont, _weightSum_Ups1;
279 CounterPtr _mult[2][12];
280 /// @}
281
282
283 };
284
285
286 RIVET_DECLARE_PLUGIN(CLEO_1985_I205668);
287
288}
|