Rivet analyses referenceCLEOC_2008_I777917Measurement of charm final-states, the total hadronic cross section and $R$ for energies between 3.92 and 4.26 GeVExperiment: CLEOC (CESR) Inspire ID: 777917 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Measurement of charm final-states, the total hadronic cross section and $R$ for energies between 3.92 and 4.26 GeV The muonic cross section is also outputted to the yoda file so that ratio $R$ can be recalcuated if runs are combined. Source code: CLEOC_2008_I777917.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 Add a short analysis description here
10 class CLEOC_2008_I777917 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(CLEOC_2008_I777917);
15
16
17 /// @name Analysis methods
18 //@{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22
23 // Initialise and register projections
24 declare(FinalState(), "FS");
25 declare(UnstableParticles(), "UFS");
26 book(_c_hadrons, "/TMP/sigma_hadrons");
27 book(_c_muons, "/TMP/sigma_muons");
28 book(_c_D0D0, "/TMP/sigma_D0D0");
29 book(_c_DpDm, "/TMP/sigma_DpDm");
30 book(_c_DsDs, "/TMP/sigma_DsDs");
31 book(_c_D0D0S, "/TMP/sigma_D0D0S");
32 book(_c_DpDmS, "/TMP/sigma_DpDmS");
33 book(_c_DsDsS, "/TMP/sigma_DsDsS");
34 book(_c_D0SD0S, "/TMP/sigma_D0SD0S");
35 book(_c_DpSDmS, "/TMP/sigma_DpSDmS");
36 book(_c_DsSDsS, "/TMP/sigma_DsSDsS");
37 book(_c_DD, "/TMP/sigma_DD");
38 book(_c_DDX, "/TMP/sigma_DDX");
39 book(_c_DSDpi, "/TMP/sigma_DSDpi");
40 book(_c_DSDSpi, "/TMP/sigma_DSDSpi");
41 }
42
43 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
44 for (const Particle &child : p.children()) {
45 if(child.children().empty()) {
46 nRes[child.pid()]-=1;
47 --ncount;
48 }
49 else
50 findChildren(child,nRes,ncount);
51 }
52 }
53
54
55 /// Perform the per-event analysis
56 void analyze(const Event& event) {
57 const FinalState& fs = apply<FinalState>(event, "FS");
58 // total hadronic and muonic cross sections
59 map<long,int> nCount;
60 int ntotal(0);
61 for (const Particle& p : fs.particles()) {
62 nCount[p.pid()] += 1;
63 ++ntotal;
64 }
65 // mu+mu- + photons
66 if(nCount[-13]==1 and nCount[13]==1 &&
67 ntotal==2+nCount[22])
68 _c_muons->fill();
69 // everything else
70 else
71 _c_hadrons->fill();
72 // identified final state with D mesons
73 const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
74 for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
75 const Particle& p1 = ufs.particles()[ix];
76 int id1 = abs(p1.pid());
77 if(id1 != 411 && id1 != 413 && id1 != 421 && id1 != 423 &&
78 id1 != 431 && id1 != 433)
79 continue;
80 // check fs
81 bool fs = true;
82 for (const Particle & child : p1.children()) {
83 if(child.pid()==p1.pid()) {
84 fs = false;
85 break;
86 }
87 }
88 if(!fs) continue;
89 // find the children
90 map<long,int> nRes = nCount;
91 int ncount = ntotal;
92 findChildren(p1,nRes,ncount);
93 bool matched=false;
94 int sign = p1.pid()/id1;
95 // loop over the other fs particles
96 for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
97 const Particle& p2 = ufs.particles()[iy];
98 fs = true;
99 for (const Particle & child : p2.children()) {
100 if(child.pid()==p2.pid()) {
101 fs = false;
102 break;
103 }
104 }
105 if(!fs) continue;
106 if(p2.pid()/abs(p2.pid())==sign) continue;
107 int id2 = abs(p2.pid());
108 if(id2 != 411 && id2 != 413 && id2 != 421 && id2 != 423 &&
109 id2 != 431 && id2 != 433)
110 continue;
111 if(!p2.parents().empty() && p2.parents()[0].pid()==p1.pid())
112 continue;
113 if((id1==411 || id1==421 || id1==431) && (id2==411 || id2==421 || id2==431 ))
114 _c_DDX->fill();
115 map<long,int> nRes2 = nRes;
116 int ncount2 = ncount;
117 findChildren(p2,nRes2,ncount2);
118 if(ncount2==0) {
119 matched=true;
120 for(auto const & val : nRes2) {
121 if(val.second!=0) {
122 matched = false;
123 break;
124 }
125 }
126 if(matched) {
127 if(id1==411 && id2==411) {
128 _c_DpDm->fill();
129 _c_DD ->fill();
130 }
131 else if(id1==421&& id2==421) {
132 _c_D0D0->fill();
133 _c_DD ->fill();
134 }
135 else if(id1==431&& id2==431) {
136 _c_DsDs->fill();
137 }
138 else if(id1==413 && id2==413) {
139 _c_DpSDmS->fill();
140 }
141 else if(id1==423&& id2==423) {
142 _c_D0SD0S->fill();
143 }
144 else if(id1==433&& id2==433) {
145 _c_DsSDsS->fill();
146 }
147 else if((id1==421 && id2==423) ||
148 (id1==423 && id2==421)) {
149 _c_D0D0S->fill();
150 }
151 else if((id1==411 && id2==413) ||
152 (id1==413 && id2==411)) {
153 _c_DpDmS->fill();
154 }
155 else if((id1==431 && id2==433) ||
156 (id1==433 && id2==431)) {
157 _c_DsDsS->fill();
158 }
159 }
160 }
161 else if(ncount2==1) {
162 int ipi=0;
163 if(nRes2[111]==1 && nRes2[211]==0 && nRes[-211]==0 )
164 ipi = 111;
165 else if(nRes2[111]==0 && nRes2[211]==1 && nRes[-211]==0 )
166 ipi = 211;
167 else if(nRes2[111]==0 && nRes2[211]==0 && nRes[-211]==1 )
168 ipi =-211;
169 if(ipi==0) continue;
170 matched=true;
171 for(auto const & val : nRes2) {
172 if(val.first==ipi)
173 continue;
174 else if(val.second!=0) {
175 matched = false;
176 break;
177 }
178 }
179 if(matched) {
180 bool Ddecay = false;
181 Particle mother = p1;
182 while (!mother.parents().empty()) {
183 mother = mother.parents()[0];
184 if(PID::isCharmMeson(mother.pid()) && mother.pid()!=p1.pid()) {
185 Ddecay = true;
186 break;
187 }
188 }
189 mother = p2;
190 while (!mother.parents().empty()) {
191 mother = mother.parents()[0];
192 if(PID::isCharmMeson(mother.pid()) && mother.pid()!=p1.pid()) {
193 Ddecay = true;
194 break;
195 }
196 }
197 if(Ddecay) continue;
198 if((id1==413 || id1==423 ) &&
199 (id2==413 || id2==423 )) {
200 _c_DSDSpi->fill();
201 }
202 else if((id1==411 || id1==421 ) &&
203 (id2==413 || id2==423 )) {
204 _c_DSDpi->fill();
205 }
206 else if((id1==413 || id1==423 ) &&
207 (id2==411 || id2==421 )) {
208 _c_DSDpi->fill();
209 }
210 }
211 }
212 }
213 }
214 }
215
216 /// Normalise histograms etc., after the run
217 void finalize() {
218 // R
219 Scatter1D R = *_c_hadrons/ *_c_muons;
220 double rval = R.point(0).x();
221 pair<double,double> rerr = R.point(0).xErrs();
222 double fact = crossSection()/ sumOfWeights() /nanobarn;
223 double sig_h = _c_hadrons->val()*fact;
224 double err_h = _c_hadrons->err()*fact;
225 double sig_c = _c_DDX->val()*fact;
226 double err_c = _c_DDX->err()*fact;
227 double sig_m = _c_muons ->val()*fact;
228 double err_m = _c_muons ->err()*fact;
229 Scatter2D temphisto(refData(6, 1, 2));
230 Scatter2DPtr charm;
231 book(charm, 6,1,1);
232 Scatter2DPtr hadrons;
233 book(hadrons, "sigma_hadrons");
234 Scatter2DPtr muons;
235 book(muons, "sigma_muons" );
236 Scatter2DPtr mult;
237 book(mult, 6,1,2);
238 for (size_t b = 0; b < temphisto.numPoints(); b++) {
239 const double x = temphisto.point(b).x();
240 pair<double,double> ex = temphisto.point(b).xErrs();
241 pair<double,double> ex2 = ex;
242 if(ex2.first ==0.) ex2. first=0.0001;
243 if(ex2.second==0.) ex2.second=0.0001;
244 if (inRange(sqrtS()/MeV, x-ex2.first, x+ex2.second)) {
245 mult ->addPoint(x, rval, ex, rerr);
246 hadrons->addPoint(x, sig_h, ex, make_pair(err_h,err_h));
247 charm ->addPoint(x, sig_c, ex, make_pair(err_c,err_c));
248 muons ->addPoint(x, sig_m, ex, make_pair(err_m,err_m));
249 }
250 else {
251 mult ->addPoint(x, 0., ex, make_pair(0.,.0));
252 hadrons->addPoint(x, 0., ex, make_pair(0.,.0));
253 charm ->addPoint(x, 0., ex, make_pair(0.,.0));
254 muons ->addPoint(x, 0., ex, make_pair(0.,.0));
255 }
256 }
257 for(unsigned int ix=1;ix<6;++ix) {
258 unsigned int imax(0);
259 if (ix<=3) imax = 4;
260 else imax = 3;
261 for(unsigned int iy=1;iy<imax;++iy) {
262 double sigma(0),error(0);
263 if(ix==1) {
264 if(iy==1) {
265 sigma = _c_D0D0->val()/picobarn;
266 error = _c_D0D0->err()/picobarn;
267 }
268 else if(iy==2) {
269 sigma = _c_D0D0S->val()/picobarn;
270 error = _c_D0D0S->err()/picobarn;
271 }
272 else if(iy==3) {
273 sigma = _c_D0SD0S->val()/picobarn;
274 error = _c_D0SD0S->err()/picobarn;
275 }
276 }
277 else if(ix==2) {
278 if(iy==1) {
279 sigma = _c_DpDm->val()/picobarn;
280 error = _c_DpDm->err()/picobarn;
281 }
282 else if(iy==2) {
283 sigma = _c_DpDmS->val()/picobarn;
284 error = _c_DpDmS->err()/picobarn;
285 }
286 else if(iy==3) {
287 sigma = _c_DpSDmS->val()/picobarn;
288 error = _c_DpSDmS->err()/picobarn;
289 }
290 }
291 else if(ix==3) {
292 if(iy==1) {
293 sigma = _c_DsDs->val()/picobarn;
294 error = _c_DsDs->err()/picobarn;
295 }
296 else if(iy==2) {
297 sigma = _c_DsDsS->val()/picobarn;
298 error = _c_DsDsS->err()/picobarn;
299 }
300 else if(iy==3) {
301 sigma = _c_DsSDsS->val()/picobarn;
302 error = _c_DsSDsS->err()/picobarn;
303 }
304 }
305 else if(ix==4) {
306 if(iy==1) {
307 sigma = _c_DSDpi->val()/picobarn;
308 error = _c_DSDpi->err()/picobarn;
309 }
310 else if(iy==2) {
311 sigma = _c_DSDSpi->val()/picobarn;
312 error = _c_DSDSpi->err()/picobarn;
313 }
314 }
315 else if(ix==5) {
316 if(iy==1) {
317 sigma = _c_DD->val()/nanobarn;
318 error = _c_DD->err()/nanobarn;
319 }
320 else if(iy==2) {
321 sigma = _c_DDX->val()/nanobarn;
322 error = _c_DDX->err()/nanobarn;
323 }
324 }
325 sigma *= crossSection()/ sumOfWeights();
326 error *= crossSection()/ sumOfWeights();
327 Scatter2D temphisto(refData(ix, 1, iy));
328 Scatter2DPtr mult;
329 book(mult, ix,1,iy);
330 for (size_t b = 0; b < temphisto.numPoints(); b++) {
331 const double x = temphisto.point(b).x();
332 pair<double,double> ex = temphisto.point(b).xErrs();
333 pair<double,double> ex2 = ex;
334 if(ex2.first ==0.) ex2. first=0.0001;
335 if(ex2.second==0.) ex2.second=0.0001;
336 if (inRange(sqrtS()/MeV, x-ex2.first, x+ex2.second)) {
337 mult ->addPoint(x, sigma, ex, make_pair(error,error));
338 }
339 else {
340 mult ->addPoint(x, 0., ex, make_pair(0.,.0));
341 }
342 }
343 }
344 }
345 }
346
347 //@}
348
349
350 /// @name Histograms
351 //@{
352 CounterPtr _c_D0D0, _c_DpDm,_c_DsDs;
353 CounterPtr _c_D0D0S, _c_DpDmS,_c_DsDsS;
354 CounterPtr _c_D0SD0S, _c_DpSDmS,_c_DsSDsS;
355 CounterPtr _c_DD, _c_DDX;
356 CounterPtr _c_DSDpi, _c_DSDSpi;
357 CounterPtr _c_hadrons, _c_muons;
358 //@}
359
360
361 };
362
363
364 // The hook for the plugin system
365 RIVET_DECLARE_PLUGIN(CLEOC_2008_I777917);
366
367
368}
|