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