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