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.48, 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.48, 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]["pi0"];
78 _axes[1]["K0"] = YODA::Axis<double>({0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
79 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9});
80 _axes[1]["lam"] = _axes[0]["lam"];
81 _axes[1]["xi"] = _axes[0]["xi"];
82 _axes[1]["rho"] = YODA::Axis<double>({0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7});
83 _axes[1]["Kstarp"] = _axes[0]["Kstarp"];
84 _axes[1]["Kstar0"] = YODA::Axis<double>({0.06, 0.12, 0.24, 0.36, 0.48});
85 _axes[1]["phi"] = YODA::Axis<double>({0.28, 0.36, 0.7, 1.0});
86
87 }
88
89 /// Recursively walk the decay tree to find decay products of @a p
90 void findDecayProducts(Particle mother, Particles& unstable) {
91 for(const Particle & p: mother.children()) {
92 const int id = p.abspid();
93 if (id == PID::PIPLUS || id==PID::KPLUS || id==PID::PROTON ||
94 id==PID::PI0 || id==PID::K0S || id==PID::K0L ||
95 id==PID::LAMBDA || id==PID::XIMINUS || id==PID::RHO0 ||
96 id==323 || id==313 || id==225 || id==PID::PHI) {
97 unstable.push_back(p);
98 }
99 if(!p.children().empty())
100 findDecayProducts(p, unstable);
101 }
102 }
103
104 /// Perform the per-event analysis
105 void analyze(const Event& event) {
106 if (_edges[0].empty()) {
107 for (const auto& item : _cont) {
108 _edges[0][item.first] = item.second->xEdges();
109 _edges[1][item.first] = _ups1[item.first]->xEdges();
110 }
111 }
112 // Find the upsilons
113 // First in unstable final state
114 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
115 Particles upsilons = ufs.particles(Cuts::pid==553);
116 // continuum
117 if (upsilons.empty()) {
118 _weightSum_cont->fill();
119 const FinalState& fs = apply<FinalState>(event, "FS");
120 // FS particles
121 for (const Particle& p : fs.particles()) {
122 int id = p.abspid();
123 double xp = 2.*p.p3().mod()/sqrtS();
124 if(id==PID::PIPLUS) {
125 discfill("pip", xp, 0);
126 _mult[1][0]->fill();
127 }
128 else if(id==PID::KPLUS) {
129 discfill("Kp", xp, 0);
130 _mult[1][1]->fill();
131 }
132 else if(id==PID::PROTON) {
133 discfill("p", xp, 0);
134 _mult[1][2]->fill();
135 }
136 }
137 // Unstable particles
138 for (const Particle& p : ufs.particles()) {
139 int id = p.abspid();
140 double xp = 2.*p.p3().mod()/sqrtS();
141 if(id==PID::PI0) {
142 discfill("pi0", xp, 0);
143 _mult[1][3]->fill();
144 }
145 else if(id==PID::K0S || id==PID::K0L) {
146 discfill("K0", xp, 0);
147 _mult[1][4]->fill();
148 }
149 else if(id==PID::LAMBDA) {
150 discfill("lam", xp, 0);
151 _mult[1][5]->fill();
152 }
153 else if(id==PID::XIMINUS) {
154 discfill("xi", xp, 0);
155 _mult[1][6]->fill();
156 }
157 else if(id==PID::RHO0) {
158 discfill("rho", xp, 0);
159 _mult[1][7]->fill();
160 }
161 else if(id==323) {
162 discfill("Kstarp", xp, 0);
163 _mult[1][8]->fill();
164 }
165 else if(id==313) {
166 discfill("Kstar0", xp, 0);
167 _mult[1][9]->fill();
168 }
169 else if(id==PID::PHI) {
170 discfill("phi", xp, 0);
171 _mult[1][10]->fill();
172 }
173 else if(id==225) {
174 _mult[1][11]->fill();
175 }
176 }
177 }
178 else {
179 for (const Particle& ups : upsilons) {
180 _weightSum_Ups1->fill();
181 Particles unstable;
182 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
183 // Find the decay products we want
184 findDecayProducts(ups,unstable);
185 for (const Particle& p : unstable) {
186 int id = p.abspid();
187 double xp = 2.*boost.transform(p.momentum()).p3().mod()/ups.mass();
188 if(id==PID::PIPLUS) {
189 discfill("pip", xp, 1);
190 _mult[0][0]->fill();
191 }
192 else if(id==PID::KPLUS) {
193 discfill("Kp", xp, 1);
194 _mult[0][1]->fill();
195 }
196 else if(id==PID::PROTON) {
197 discfill("p", xp, 1);
198 _mult[0][2]->fill();
199 }
200 else if(id==PID::PI0) {
201 discfill("pi0", xp, 1);
202 _mult[0][3]->fill();
203 }
204 else if(id==PID::K0S || id==PID::K0L) {
205 discfill("K0", xp, 1);
206 _mult[0][4]->fill();
207 }
208 else if(id==PID::LAMBDA) {
209 discfill("lam", xp, 1);
210 _mult[0][5]->fill();
211 }
212 else if(id==PID::XIMINUS) {
213 discfill("xi", xp, 1);
214 _mult[0][6]->fill();
215 }
216 else if(id==PID::RHO0) {
217 discfill("rho", xp, 1);
218 _mult[0][7]->fill();
219 }
220 else if(id==323) {
221 discfill("Kstarp", xp, 1);
222 _mult[0][8]->fill();
223 }
224 else if(id==313) {
225 discfill("Kstar0", xp, 1);
226 _mult[0][9]->fill();
227 }
228 else if(id==PID::PHI) {
229 discfill("phi", xp, 1);
230 _mult[0][10]->fill();
231 }
232 else if(id==225) {
233 _mult[0][11]->fill();
234 }
235 }
236 }
237 }
238 }
239
240 void discfill(const string& name, const double value, const size_t k) {
241 string edge = "OTHER";
242 size_t idx = _axes[k][name].index(value);
243 if (name=="pip") {
244 if(idx==8) idx=0;
245 else if(idx>8) idx-=1;
246 }
247 if (idx && idx <= _edges[k][name].size()) edge = _edges[k][name][idx-1];
248 (k? _ups1 : _cont)[name]->fill(edge);
249 }
250
251
252 /// Normalise histograms etc., after the run
253 void finalize() {
254 // multiplicities
255 const vector<CounterPtr> scales = {_weightSum_Ups1, _weightSum_cont};
256 for (size_t ix=0; ix<12; ++ix) {
257 BinnedEstimatePtr<string> est;
258 book(est, ix+12, 1, 1);
259 for (size_t iy=0; iy<2; ++iy) {
260 if (scales[iy]->val() > 0.) {
261 unsigned int iz = iy==0 ? 2 : 1;
262 scale(_mult[iy][ix], 1./ *scales[iy]);
263 est->bin(iz).set(_mult[iy][ix]->val(), _mult[iy][ix]->err());
264 }
265 }
266 }
267 // spectra
268 if (_weightSum_cont->val() > 0.) {
269 scale(_cont, 1. / *_weightSum_cont);
270 for( auto & hist : _cont) {
271 for(auto & b: hist.second->bins()) {
272 size_t idx = b.index();
273 if(hist.first=="pip" && idx>=8) idx+=1;
274 b.scaleW(1./_axes[0][hist.first].width(idx));
275 }
276 }
277 }
278 if (_weightSum_Ups1->val() > 0.) {
279 scale(_ups1, 1. / *_weightSum_Ups1);
280 for( auto & hist : _ups1) {
281 for(auto & b: hist.second->bins()) {
282 size_t idx = b.index();
283 if(hist.first=="pip" && idx>=8) idx+=1;
284 b.scaleW(1./_axes[1][hist.first].width(idx));
285 }
286 }
287 }
288 }
289
290 /// @}
291
292
293 /// @name Histograms
294 /// @{
295 map<string,BinnedHistoPtr<string>> _cont, _ups1;
296 map<string, YODA::Axis<double>> _axes[2];
297 map<string, vector<string>> _edges[2];
298 CounterPtr _weightSum_cont, _weightSum_Ups1;
299 CounterPtr _mult[2][12];
300 /// @}
301
302
303 };
304
305
306 RIVET_DECLARE_PLUGIN(CLEO_1985_I205668);
307
308}
|