Rivet analyses referenceBESIII_2019_I1765606Analysis of $J/\psi$ decays to $\Xi^{*-}\bar\Xi^{+}+\text{c.c.}$Experiment: BESIII (BEPC) Inspire ID: 1765606 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 \Xi^{*-}\bar\Xi^{+}+\text{c.c.}$. Gives information about the decay and is useful for testing correlations in hadron decays. Source code: BESIII_2019_I1765606.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 BES psi2s -> Xi*- Xibar+ + c.c.
11 class BESIII_2019_I1765606 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2019_I1765606);
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_cTheta, 1, 1, 1);
31
32 }
33
34 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
35 for(const Particle &child : p.children()) {
36 if(child.children().empty()) {
37 nRes[child.pid()]-=1;
38 --ncount;
39 }
40 else
41 findChildren(child,nRes,ncount);
42 }
43 }
44
45 /// Perform the per-event analysis
46 void analyze(const Event& event) {
47 // get the axis, direction of incoming electron
48 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
49 Vector3 axis;
50 if(beams.first.pid()<0)
51 axis = beams.first .momentum().p3().unit();
52 else
53 axis = beams.second.momentum().p3().unit();
54 // types of final state particles
55 const FinalState& fs = apply<FinalState>(event, "FS");
56 map<long,int> nCount;
57 int ntotal(0);
58 for (const Particle& p : fs.particles()) {
59 nCount[p.pid()] += 1;
60 ++ntotal;
61 }
62 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
63 for (const Particle& p : ufs.particles(Cuts::abspid==3314)) {
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 int sign = p.pid()/3314;
70 for (const Particle& p2 : ufs.particles(Cuts::pid==-sign*3312)) {
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 cTheta= (sign>0) ? p .momentum().p3().unit().dot(axis) : p2.momentum().p3().unit().dot(axis);
87 _h_cTheta->fill(cTheta);
88 break;
89 }
90
91 }
92 }
93 }
94 }
95
96 pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
97 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
98 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
99 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
100 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
101 for (const auto& bin : hist->bins() ) {
102 double Oi = bin.sumW();
103 if(Oi==0.) continue;
104 double a = d*(bin.xMax() - bin.xMin());
105 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
106 double Ei = bin.errW();
107 sum1 += a*Oi/sqr(Ei);
108 sum2 += b*Oi/sqr(Ei);
109 sum3 += sqr(a)/sqr(Ei);
110 sum4 += sqr(b)/sqr(Ei);
111 sum5 += a*b/sqr(Ei);
112 }
113 // calculate alpha
114 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
115 // and error
116 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
117 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
118 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
119 double dis = sqr(bb)-4.*aa*cc;
120 if(dis>0.) {
121 dis = sqrt(dis);
122 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
123 }
124 else {
125 return make_pair(alpha,make_pair(0.,0.));
126 }
127 }
128
129 /// Normalise histograms etc., after the run
130 void finalize() {
131 normalize(_h_cTheta);
132 Estimate0DPtr h_alpha_xi;
133 book(h_alpha_xi,2,1,1);
134 pair<double,pair<double,double> > alpha = calcAlpha(_h_cTheta);
135 h_alpha_xi->set(alpha.first, alpha.second);
136 }
137
138 /// @}
139
140
141 /// @name Histograms
142 /// @{
143 Histo1DPtr _h_cTheta;
144 /// @}
145
146 };
147
148
149 RIVET_DECLARE_PLUGIN(BESIII_2019_I1765606);
150
151}
|