Rivet analyses referenceBESIII_2017_I1510563Analysis of $J/\psi$ and $\psi(2S)$ decays to $\Lambda^0\bar\Lambda^0$ and $\Sigma^0\bar\Sigma^0$Experiment: BESIII (BEPC) Inspire ID: 1510563 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 \Lambda^0\bar\Lambda^0, \Sigma^0\bar\Sigma^0$. Gives information about the decay and is useful for testing correlations in hadron decay. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: BESIII_2017_I1510563.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_2017_I1510563 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2017_I1510563);
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,1e-1)) {
31 book(_h_lam, 1, 1, 1);
32 book(_h_sig, 1, 1, 3);
33 }
34 else if (isCompatibleWithSqrtS(3.686, 1E-1)) {
35 book(_h_lam, 1, 1, 2);
36 book(_h_sig, 1, 1, 4);
37 }
38 }
39
40 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
41 for( const Particle &child : p.children()) {
42 if(child.children().empty()) {
43 nRes[child.pid()]-=1;
44 --ncount;
45 }
46 else
47 findChildren(child,nRes,ncount);
48 }
49 }
50
51 /// Perform the per-event analysis
52 void analyze(const Event& event) {
53 // get the axis, direction of incoming electron
54 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
55 Vector3 axis;
56 if(beams.first.pid()>0)
57 axis = beams.first .momentum().p3().unit();
58 else
59 axis = beams.second.momentum().p3().unit();
60 // types of final state particles
61 const FinalState& fs = apply<FinalState>(event, "FS");
62 map<long,int> nCount;
63 int ntotal(0);
64 for (const Particle& p : fs.particles()) {
65 nCount[p.pid()] += 1;
66 ++ntotal;
67 }
68 // loop over lambda0 and sigma0 baryons
69 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
70 for (const Particle& p : ufs.particles(Cuts::abspid==3122 or Cuts::abspid==3212)) {
71 if(p.children().empty()) continue;
72 map<long,int> nRes=nCount;
73 int ncount = ntotal;
74 findChildren(p,nRes,ncount);
75 bool matched=false;
76 // check for antiparticle
77 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
78 if(p2.children().empty()) continue;
79 map<long,int> nRes2=nRes;
80 int ncount2 = ncount;
81 findChildren(p2,nRes2,ncount2);
82 if(ncount2==0) {
83 matched = true;
84 for(auto const & val : nRes2) {
85 if(val.second!=0) {
86 matched = false;
87 break;
88 }
89 }
90 // fond baryon and antibaryon
91 if(matched) {
92 // calc cosine
93 double ctheta;
94 if(p.pid()>0)
95 ctheta = p .momentum().p3().unit().dot(axis);
96 else
97 ctheta = p2.momentum().p3().unit().dot(axis);
98 if(abs(p.pid())==3122)
99 _h_lam->fill(ctheta);
100 else
101 _h_sig->fill(ctheta);
102 break;
103 }
104 }
105 }
106 if(matched) break;
107 }
108 }
109
110 pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
111 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
112 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
113 for (auto bin : hist->bins() ) {
114 double Oi = bin.area();
115 if(Oi==0.) continue;
116 double a = 1.5*(bin.xMax() - bin.xMin());
117 double b = 0.5*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
118 double Ei = bin.areaErr();
119 sum1 += a*Oi/sqr(Ei);
120 sum2 += b*Oi/sqr(Ei);
121 sum3 += sqr(a)/sqr(Ei);
122 sum4 += sqr(b)/sqr(Ei);
123 sum5 += a*b/sqr(Ei);
124 }
125 // calculate alpha
126 double alpha = (-3*sum1 + 9*sum2 + sum3 - 3*sum5)/(sum1 - 3*sum2 + 3*sum4 - sum5);
127 // and error
128 double cc = -pow((sum3 + 9*sum4 - 6*sum5),3);
129 double bb = -2*sqr(sum3 + 9*sum4 - 6*sum5)*(sum1 - 3*sum2 + 3*sum4 - sum5);
130 double aa = sqr(sum1 - 3*sum2 + 3*sum4 - sum5)*(-sum3 - 9*sum4 + sqr(sum1 - 3*sum2 + 3*sum4 - sum5) + 6*sum5);
131 double dis = sqr(bb)-4.*aa*cc;
132 if(dis>0.) {
133 dis = sqrt(dis);
134 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
135 }
136 else {
137 return make_pair(alpha,make_pair(0.,0.));
138 }
139 }
140
141 /// Normalise histograms etc., after the run
142 void finalize() {
143 // find energy
144 int ioff=-1;
145 if (isCompatibleWithSqrtS(3.1,1e-1)) ioff=0;
146 else if (isCompatibleWithSqrtS(3.686, 1E-1)) ioff=1;
147 // normalize
148 normalize(_h_lam);
149 pair<double,pair<double,double> > alpha = calcAlpha(_h_lam);
150 Scatter2DPtr _h_alpha_lam;
151 book(_h_alpha_lam, 2,2*ioff+1,1);
152 _h_alpha_lam->addPoint(0.5, alpha.first, make_pair(0.5,0.5),
153 make_pair(alpha.second.first,alpha.second.second) );
154 normalize(_h_sig);
155 alpha = calcAlpha(_h_sig);
156 Scatter2DPtr _h_alpha_sig;
157 book(_h_alpha_sig, 2,2*ioff+2,1);
158 _h_alpha_sig->addPoint(0.5, alpha.first, make_pair(0.5,0.5),
159 make_pair(alpha.second.first,alpha.second.second) );
160 }
161
162 //@}
163
164
165 /// @name Histograms
166 //@{
167 Histo1DPtr _h_lam,_h_sig;
168 //@}
169
170
171 };
172
173 // The hook for the plugin system
174 RIVET_DECLARE_PLUGIN(BESIII_2017_I1510563);
175
176}
|