Rivet analyses referenceBELLE_2011_I878228$e^+e^-\to D^{(*)+}_sD^{(*)-}_s$ cross sections near thresholdExperiment: BELLE (KEKB) Inspire ID: 878228 Status: VALIDATED Authors:
Beam energies: ANY Run details:
$e^+e^-\to D^{(*)+}_sD^{(*)-}_s$ cross sections measured near threshold by BELLE using ISR. Cross sections for individual states, and $R$ for the sum Source code: BELLE_2011_I878228.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 e+e- > Ds(*)+ Ds*()-
10 class BELLE_2011_I878228 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2011_I878228);
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
27 // Book histograms
28 book( _c_DpDm, "/TMP/sigma_DpDm" );
29 book( _c_DpDmS, "/TMP/sigma_DpDmS" );
30 book(_c_DpSDmS, "/TMP/sigma_DpSDmS");
31 book( _c_All, "/TMP/sigma_All" );
32 book( _c_mu, "/TMP/sigma_mu" );
33 }
34
35 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
36 for(const Particle &child : p.children()) {
37 if(child.children().empty()) {
38 nRes[child.pid()]-=1;
39 --ncount;
40 }
41 else
42 findChildren(child,nRes,ncount);
43 }
44 }
45
46 /// Perform the per-event analysis
47 void analyze(const Event& event) {
48 const FinalState& fs = apply<FinalState>(event, "FS");
49 // total hadronic and muonic cross sections
50 map<long,int> nCount;
51 int ntotal(0);
52 for (const Particle& p : fs.particles()) {
53 nCount[p.pid()] += 1;
54 ++ntotal;
55 }
56 // mu+mu- + photons
57 if(nCount[-13]==1 and nCount[13]==1 &&
58 ntotal==2+nCount[22]) {
59 _c_mu->fill();
60 return;
61 }
62 // unstable charm analysis
63
64 Particles ds = apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==431 or Cuts::abspid==433);
65 for(unsigned int ix=0;ix<ds.size();++ix) {
66 const Particle& p1 = ds[ix];
67 int id1 = abs(p1.pid());
68 // check fs
69 bool fs = true;
70 for (const Particle & child : p1.children()) {
71 if(child.pid()==p1.pid()) {
72 fs = false;
73 break;
74 }
75 }
76 if(!fs) continue;
77 // find the children
78 map<long,int> nRes = nCount;
79 int ncount = ntotal;
80 findChildren(p1,nRes,ncount);
81 bool matched=false;
82 int sign = p1.pid()/id1;
83 // loop over the other fs particles
84 for(unsigned int iy=ix+1;iy<ds.size();++iy) {
85 const Particle& p2 = ds[iy];
86 fs = true;
87 for (const Particle & child : p2.children()) {
88 if(child.pid()==p2.pid()) {
89 fs = false;
90 break;
91 }
92 }
93 if(!fs) continue;
94 if(p2.pid()/abs(p2.pid())==sign) continue;
95 int id2 = abs(p2.pid());
96 if(!p2.parents().empty() && p2.parents()[0].pid()==p1.pid())
97 continue;
98 map<long,int> nRes2 = nRes;
99 int ncount2 = ncount;
100 findChildren(p2,nRes2,ncount2);
101 if(ncount2!=0) continue;
102 matched=true;
103 for(auto const & val : nRes2) {
104 if(val.second!=0) {
105 matched = false;
106 break;
107 }
108 }
109 if(matched) {
110 if(id1==431 && id2==431) {
111 _c_DpDm->fill();
112 }
113 else if(id1==433 && id2==433) {
114 _c_DpSDmS->fill();
115 }
116 else if((id1==431 && id2==433) ||
117 (id1==433 && id2==431)) {
118 _c_DpDmS->fill();
119 }
120 _c_All->fill();
121 break;
122 }
123 }
124 if(matched) break;
125 }
126 }
127
128
129 /// Normalise histograms etc., after the run
130 void finalize() {
131 double fact = crossSection()/ sumOfWeights()/nanobarn;
132 for(unsigned int iy=1;iy<4;++iy) {
133 double sigma = 0.0, error = 0.0;
134 if(iy==1) {
135 sigma = _c_DpDm->val()*fact;
136 error = _c_DpDm->err()*fact;
137 }
138 else if(iy==2) {
139 sigma = _c_DpDmS->val()*fact;
140 error = _c_DpDmS->err()*fact;
141 }
142 else if(iy==3) {
143 sigma = _c_DpSDmS->val()*fact;
144 error = _c_DpSDmS->err()*fact;
145 }
146 Scatter2D temphisto(refData(1, 1, iy));
147 Scatter2DPtr mult;
148 book(mult, 1, 1, iy);
149 for (size_t b = 0; b < temphisto.numPoints(); b++) {
150 const double x = temphisto.point(b).x();
151 pair<double,double> ex = temphisto.point(b).xErrs();
152 pair<double,double> ex2 = ex;
153 if(ex2.first ==0.) ex2. first=0.0001;
154 if(ex2.second==0.) ex2.second=0.0001;
155 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
156 mult ->addPoint(x, sigma, ex, make_pair(error,error));
157 }
158 else {
159 mult ->addPoint(x, 0., ex, make_pair(0.,.0));
160 }
161 }
162 }
163 // R
164 if(_c_mu->val()==0. || _c_All->val()==0.) return;
165 Scatter1D R = *_c_All/ *_c_mu;
166 double rval = R.point(0).x();
167 pair<double,double> rerr = R.point(0).xErrs();
168 Scatter2D temphisto(refData(2, 1, 1));
169 Scatter2DPtr mult;
170 book(mult, 2,1,1);
171 for (size_t b = 0; b < temphisto.numPoints(); b++) {
172 const double x = temphisto.point(b).x();
173 pair<double,double> ex = temphisto.point(b).xErrs();
174 pair<double,double> ex2 = ex;
175 if(ex2.first ==0.) ex2. first=0.0001;
176 if(ex2.second==0.) ex2.second=0.0001;
177 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
178 mult ->addPoint(x, rval, ex, rerr);
179 }
180 else {
181 mult ->addPoint(x, 0., ex, make_pair(0.,.0));
182 }
183 }
184 }
185
186 //@}
187
188
189 /// @name Histograms
190 //@{
191 CounterPtr _c_DpDm,_c_DpDmS,_c_DpSDmS,_c_All,_c_mu;
192 //@}
193
194
195 };
196
197
198 RIVET_DECLARE_PLUGIN(BELLE_2011_I878228);
199
200}
|