Rivet analyses referenceBESIII_2019_I1726357Measurement of $e^+e^-\to\Lambda^0\bar{\Lambda}^0$ at 2.396 GeVExperiment: BESIII (BEPC) Inspire ID: 1726357 Status: VALIDATED Authors:
Beam energies: (1.2, 1.2) GeV Run details:
Measurement of the cross section, angular distribution and polarization for $e^+e^-\to\Lambda^0\bar{\Lambda}^0$ at 2.396 GeV by BESIII. Source code: BESIII_2019_I1726357.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 Lambda Lambdabar cross section
11 class BESIII_2019_I1726357 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2019_I1726357);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // Initialise and register projections
24 declare(Beam(), "Beams");
25 declare(FinalState(), "FS");
26 declare(UnstableParticles(), "UFS");
27 // histograms
28 book(_h_sigma ,1,1,1);
29 book(_h_cTheta,2,1,1);
30 book(_h_pol, {-1, -0.8, -0.6, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8, 1.0});
31 for (auto& b : _h_pol->bins()) {
32 book(b, "/TMP/h_pol_"+to_string(b.index()+1), 20, -1.0, 1.0);
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 const FinalState& fs = apply<FinalState>(event, "FS");
57 // total hadronic and muonic cross sections
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 // find the Lambdas
65 bool matched = false;
66 const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
67 Particle Lambda;
68 for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
69 const Particle& p1 = ufs.particles()[ix];
70 if(abs(p1.pid())!=3122) continue;
71 // check fs
72 bool fs = true;
73 for (const Particle & child : p1.children()) {
74 if(child.pid()==p1.pid()) {
75 fs = false;
76 break;
77 }
78 }
79 if(!fs) continue;
80 // find the children
81 map<long,int> nRes = nCount;
82 int ncount = ntotal;
83 findChildren(p1,nRes,ncount);
84 for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
85 matched=false;
86 const Particle& p2 = ufs.particles()[iy];
87 if(abs(p2.pid())!=3122) continue;
88 // check fs
89 bool fs = true;
90 for (const Particle & child : p2.children()) {
91 if(child.pid()==p2.pid()) {
92 fs = false;
93 break;
94 }
95 }
96 if (!fs) continue;
97 map<long,int> nRes2 = nRes;
98 int ncount2 = ncount;
99 findChildren(p2, nRes2, ncount2);
100 if (ncount2!=0) continue;
101 matched=true;
102 for (const auto& val : nRes2) {
103 if(val.second!=0) {
104 matched = false;
105 break;
106 }
107 }
108 if (matched) {
109 _h_sigma->fill(2.396);
110 if(p1.pid()==PID::LAMBDA) {
111 Lambda=p1;
112 }
113 else {
114 Lambda=p2;
115 }
116 break;
117 }
118 }
119 if (matched) break;
120 }
121 // now for the polarization measurements
122 if (matched) {
123 double cTheta = Lambda.momentum().p3().unit().dot(axis);
124 _h_cTheta->fill(cTheta);
125 Particle proton;
126 if(Lambda.children().size()!=2) return;
127 if(Lambda.children()[0].pid()==PID::PROTON &&
128 Lambda.children()[1].pid()==PID::PIMINUS)
129 proton = Lambda.children()[0];
130 else if(Lambda.children()[1].pid()==PID::PROTON &&
131 Lambda.children()[0].pid()==PID::PIMINUS)
132 proton = Lambda.children()[1];
133 else return;
134 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Lambda.momentum().betaVec());
135 Vector3 axis1 = boost1.transform(proton.momentum()).p3().unit();
136 double cPhi = axis1.dot(Lambda.momentum().p3().unit());
137 _h_pol->fill(cTheta, cPhi);
138 }
139 }
140
141 pair<double,double> calcAlpha(Histo1DPtr hist) {
142 if(hist->numEntries()==0.) return make_pair(0.,0.);
143 double sum1(0.),sum2(0.);
144 for (const auto& bin : hist->bins()) {
145 double Oi = bin.sumW();
146 if(Oi==0.) continue;
147 double ai = 0.5*bin.xWidth();
148 double bi = 0.5*ai*(bin.xMax()+bin.xMin());
149 double Ei = bin.errW();
150 sum1 += sqr(bi/Ei);
151 sum2 += bi/sqr(Ei)*(Oi-ai);
152 }
153 return make_pair(sum2/sum1,sqrt(1./sum1));
154 }
155
156 /// Normalise histograms etc., after the run
157 void finalize() {
158 scale(_h_sigma, crossSection()/sumOfWeights()/picobarn);
159 normalize(_h_cTheta);
160 Estimate1DPtr _h_alpha;
161 book(_h_alpha, 3, 1, 1);
162 for (auto& b : _h_pol->bins()) {
163 normalize(b);
164 pair<double,double> alpha = calcAlpha(b);
165 _h_alpha->bin(b.index()).set(alpha.first, alpha.second);
166 }
167 }
168 /// @}
169
170
171 /// @name Histograms
172 /// @{
173 Histo1DPtr _h_sigma,_h_cTheta;
174 Histo1DGroupPtr _h_pol;
175 /// @}
176
177
178 };
179
180
181 RIVET_DECLARE_PLUGIN(BESIII_2019_I1726357);
182
183}
|