Rivet analyses referenceBESIII_2016_I1422780Analysis of $J/\psi$ and $\psi(2S)$ decays to $\Xi^-\bar\Xi^+$ and $\Sigma^{*\mp}\bar\Sigma^{*\pm}$Experiment: BESIII (BEPC) Inspire ID: 1422780 Status: VALIDATED Authors:
Beam energies: (1.6, 1.6); (1.8, 1.8) GeV Run details:
Analysis of the angular distribution of the baryons produced in $e^+e^-\to J/\psi,\psi(2S) \to \Xi^-\bar\Xi^+$ and $\Sigma^{*\mp}\bar\Sigma^{*\pm}$. Gives information about the decay and is useful for testing correlations in hadron decays. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: BESIII_2016_I1422780.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/UnstableParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief Jpsi/psi2S baryon decay analysis
11 class BESIII_2016_I1422780 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2016_I1422780);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // Initialise and register projections
25 declare(Beam(), "Beams");
26 declare(UnstableParticles(), "UFS");
27 declare(FinalState(), "FS");
28
29 // Book histograms
30 if(isCompatibleWithSqrtS(3.1*GeV,1e-1)) {
31 book(_h_xi , 2, 1, 1);
32 book(_h_sigm, 2, 1, 2);
33 book(_h_sigp, 2, 1, 3);
34 }
35 else if (isCompatibleWithSqrtS(3.686*GeV, 1E-1)) {
36 book(_h_xi , 2, 1, 4);
37 book(_h_sigm, 2, 1, 5);
38 book(_h_sigp, 2, 1, 6);
39 }
40 }
41
42 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
43 for( const Particle &child : p.children()) {
44 if(child.children().empty()) {
45 nRes[child.pid()]-=1;
46 --ncount;
47 }
48 else
49 findChildren(child,nRes,ncount);
50 }
51 }
52
53 /// Perform the per-event analysis
54 void analyze(const Event& event) {
55 // get the axis, direction of incoming electron
56 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
57 Vector3 axis;
58 if(beams.first.pid()>0)
59 axis = beams.first .momentum().p3().unit();
60 else
61 axis = beams.second.momentum().p3().unit();
62 // types of final state particles
63 const FinalState& fs = apply<FinalState>(event, "FS");
64 map<long,int> nCount;
65 int ntotal(0);
66 for (const Particle& p : fs.particles()) {
67 nCount[p.pid()] += 1;
68 ++ntotal;
69 }
70
71 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
72 for (const Particle& p : ufs.particles(Cuts::abspid==3312 or Cuts::abspid==3224 or Cuts::abspid==3114)) {
73 if(p.children().empty()) continue;
74 map<long,int> nRes=nCount;
75 int ncount = ntotal;
76 findChildren(p,nRes,ncount);
77 bool matched=false;
78 // check for antiparticle
79 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
80 if(p2.children().empty()) continue;
81 map<long,int> nRes2=nRes;
82 int ncount2 = ncount;
83 findChildren(p2,nRes2,ncount2);
84 if(ncount2==0) {
85 matched = true;
86 for(auto const & val : nRes2) {
87 if(val.second!=0) {
88 matched = false;
89 break;
90 }
91 }
92 // fond baryon and antibaryon
93 if(matched) {
94 // calc cosine
95 double ctheta;
96 if(p.pid()>0)
97 ctheta = p .momentum().p3().unit().dot(axis);
98 else
99 ctheta = p2.momentum().p3().unit().dot(axis);
100 if(abs(p.pid())==3312)
101 _h_xi ->fill(ctheta);
102 else if(abs(p.pid())==3114)
103 _h_sigm->fill(ctheta);
104 else if(abs(p.pid())==3224)
105 _h_sigp->fill(ctheta);
106 break;
107 }
108 }
109 }
110 if(matched) break;
111 }
112 }
113
114 pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
115 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
116 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
117 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
118 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
119 for (const auto& bin : hist->bins() ) {
120 double Oi = bin.sumW();
121 if(Oi==0.) continue;
122 double a = d*(bin.xMax() - bin.xMin());
123 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
124 double Ei = bin.errW();
125 sum1 += a*Oi/sqr(Ei);
126 sum2 += b*Oi/sqr(Ei);
127 sum3 += sqr(a)/sqr(Ei);
128 sum4 += sqr(b)/sqr(Ei);
129 sum5 += a*b/sqr(Ei);
130 }
131 // calculate alpha
132 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
133 // and error
134 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
135 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
136 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
137 double dis = sqr(bb)-4.*aa*cc;
138 if(dis>0.) {
139 dis = sqrt(dis);
140 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
141 }
142 else {
143 return make_pair(alpha,make_pair(0.,0.));
144 }
145 }
146
147 /// Normalise histograms etc., after the run
148 void finalize() {
149 // find energy
150 //int ioff=-1;
151 //if(isCompatibleWithSqrtS(3.1*GeV,1e-1)) ioff=0;
152 //else if (isCompatibleWithSqrtS(3.686*GeV, 1E-1)) ioff=1;
153 vector< pair<double,pair<double,double> > > alpha;
154 normalize(_h_xi,1.,false);
155 alpha.push_back(calcAlpha(_h_xi));
156 normalize(_h_sigm,1.,false);
157 alpha.push_back(calcAlpha(_h_sigm));
158 normalize(_h_sigp,1.,false);
159 alpha.push_back(calcAlpha(_h_sigp));
160 Estimate1DPtr _h_alpha;
161 book(_h_alpha,1,1,3);
162 for (unsigned int ix=0;ix<3;++ix) {
163 _h_alpha->bin(ix+1).set(alpha[ix].first,
164 make_pair(alpha[ix].second.first,alpha[ix].second.second));
165 }
166 }
167 /// @}
168
169
170 /// @name Histograms
171 /// @{
172 Histo1DPtr _h_xi,_h_sigm,_h_sigp;
173 /// @}
174
175
176 };
177
178
179 RIVET_DECLARE_PLUGIN(BESIII_2016_I1422780);
180
181
182}
|