Rivet analyses referenceBESIII_2023_I2637702Transverse $\Lambda$ polarization in $e^+e^-\to\Lambda^0\bar\Lambda^0$ for $\sqrt{s}=3.68\to3.71\,$GeVExperiment: BESIII (BEPC) Inspire ID: 2637702 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.8, 1.8); (1.8, 1.8); (1.8, 1.8); (1.8, 1.8); (1.8, 1.8); (1.8, 1.8); (1.9, 1.9) GeV Run details:
Transverse $\Lambda$ polarization in $e^+e^-\to\Lambda^0\bar\Lambda^0$ for $\sqrt{s}=3.68\to3.71\,$GeV. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: BESIII_2023_I2637702.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 e+ e- > Lambda0 Lambdabar0
11 class BESIII_2023_I2637702 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2637702);
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(Cuts::abspid==3122), "UFS");
27 declare(FinalState(), "FS");
28
29 // Book histograms
30 book(_h_T2, "TMP/T2",20,-1.,1.);
31 book(_h_T3, "TMP/T3",20,-1.,1.);
32 book(_h_cThetaL,"cThetaL",20,-1.,1.);
33 book(_wsum,"TMP/wsum");
34
35 vector<string> energies({"3.68", "3.683", "3.684", "3.685", "3.687", "3.691", "3.71"});
36 for(const string& en : energies) {
37 double end = std::stod(en)*GeV;
38 if(isCompatibleWithSqrtS(end)) {
39 _ecms = en;
40 break;
41 }
42 }
43 if(!inRange(sqrtS(),3.68,3.71) && _ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
44 }
45
46 void findChildren (const Particle& p,map<long,int> & nRes, int &ncount) {
47 for (const Particle& child : p.children()) {
48 if (child.children().empty()) {
49 nRes[child.pid()]-=1;
50 --ncount;
51 }
52 else {
53 findChildren(child,nRes,ncount);
54 }
55 }
56 }
57
58 /// Perform the per-event analysis
59 void analyze(const Event& event) {
60 // get the axis, direction of incoming electron
61 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
62 Vector3 axis;
63 if (beams.first.pid()>0) {
64 axis = beams.first .mom().p3().unit();
65 }
66 else {
67 axis = beams.second.mom().p3().unit();
68 }
69 // types of final state particles
70 const FinalState& fs = apply<FinalState>(event, "FS");
71 map<long,int> nCount;
72 int ntotal(0);
73 for (const Particle& p : fs.particles()) {
74 nCount[p.pid()] += 1;
75 ++ntotal;
76 }
77 // loop over lambda0 baryons
78 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
79 Particle Lambda,LamBar;
80 bool matched(false);
81 for (const Particle& p : ufs.particles(Cuts::abspid==3122)) {
82 if (p.children().empty()) continue;
83 map<long,int> nRes=nCount;
84 int ncount = ntotal;
85 findChildren(p,nRes,ncount);
86 matched=false;
87 // check for antiparticle
88 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
89 if (p2.children().empty()) continue;
90 map<long,int> nRes2=nRes;
91 int ncount2 = ncount;
92 findChildren(p2,nRes2,ncount2);
93 if (ncount2==0) {
94 matched = true;
95 for (const auto& val : nRes2) {
96 if (val.second!=0) {
97 matched = false;
98 break;
99 }
100 }
101 // found baryon and antibaryon
102 if (matched) {
103 if (p.pid()>0) {
104 Lambda = p;
105 LamBar = p2;
106 }
107 else {
108 Lambda = p2;
109 LamBar = p;
110 }
111 break;
112 }
113 }
114 }
115 if (matched) break;
116 }
117 if (!matched) vetoEvent;
118 Particle proton;
119 matched = false;
120 for (const Particle& p : Lambda.children()) {
121 if (p.pid()==2212) {
122 matched=true;
123 proton=p;
124 }
125 else if (p.pid()==PID::PHOTON) {
126 vetoEvent;
127 }
128 }
129 if (!matched) vetoEvent;
130 Particle baryon;
131 matched = false;
132 for (const Particle & p : LamBar.children()) {
133 if (p.pid()==-2212) {
134 baryon=p;
135 matched=true;
136 }
137 else if (p.pid()==PID::PHOTON) {
138 vetoEvent;
139 }
140 }
141 if (!matched) vetoEvent;
142 // boost to the Lambda rest frame
143 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Lambda.mom().betaVec());
144 Vector3 e1z = Lambda.mom().p3().unit();
145 Vector3 e1y = e1z.cross(axis).unit();
146 Vector3 e1x = e1y.cross(e1z).unit();
147 Vector3 axis1 = boost1.transform(proton.mom()).p3().unit();
148 double n1x(e1x.dot(axis1)),n1y(e1y.dot(axis1)),n1z(e1z.dot(axis1));
149 // boost to the Lambda bar
150 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(LamBar.mom().betaVec());
151 Vector3 axis2 = boost2.transform(baryon.mom()).p3().unit();
152 double n2x(e1x.dot(axis2)),n2z(e1z.dot(axis2));
153 double cosL = axis.dot(Lambda.mom().p3().unit());
154 double sinL = sqrt(1.-sqr(cosL));
155 double T2 = -sinL*cosL*(n1x*n2z+n1z*n2x);
156 double T3 = -sinL*cosL*n1y;
157 _h_T2->fill(cosL,T2);
158 _h_T3->fill(cosL,T3);
159 _wsum->fill();
160 _h_cThetaL->fill(cosL);
161 }
162
163 pair<double,pair<double,double> > calcAlpha0(Histo1DPtr hist) {
164 if (hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
165 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
166 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
167 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
168 for (const auto& bin : hist->bins() ) {
169 double Oi = bin.sumW();
170 if (Oi==0.) continue;
171 double a = d*(bin.xMax() - bin.xMin());
172 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
173 double Ei = bin.errW();
174 sum1 += a*Oi/sqr(Ei);
175 sum2 += b*Oi/sqr(Ei);
176 sum3 += sqr(a)/sqr(Ei);
177 sum4 += sqr(b)/sqr(Ei);
178 sum5 += a*b/sqr(Ei);
179 }
180 // calculate alpha
181 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
182 // and error
183 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
184 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
185 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
186 double dis = sqr(bb)-4.*aa*cc;
187 if (dis>0.) {
188 dis = sqrt(dis);
189 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
190 }
191 else {
192 return make_pair(alpha,make_pair(0.,0.));
193 }
194 }
195
196 pair<double,double> calcCoeff(unsigned int imode, Histo1DPtr hist) {
197 if (hist->numEntries()==0.) return make_pair(0.,0.);
198 double sum1(0.),sum2(0.);
199 for (const auto& bin : hist->bins()) {
200 double Oi = bin.sumW();
201 if (Oi==0.) continue;
202 double ai(0.),bi(0.);
203 if (imode==0) {
204 bi = (pow(1.-sqr(bin.xMin()),1.5) - pow(1.-sqr(bin.xMax()),1.5))/3.;
205 }
206 else if(imode>=2 && imode<=4) {
207 bi = (pow(bin.xMin(),3)*( -5. + 3.*sqr(bin.xMin())) + pow(bin.xMax(),3)*(5. - 3.*sqr(bin.xMax())))/15.;
208 }
209 else {
210 assert(false);
211 }
212 double Ei = bin.errW();
213 sum1 += sqr(bi/Ei);
214 sum2 += bi/sqr(Ei)*(Oi-ai);
215 }
216 return make_pair(sum2/sum1,sqrt(1./sum1));
217 }
218
219 /// Normalise histograms etc., after the run
220 void finalize() {
221 const double aLambda = 0.754;
222 normalize(_h_cThetaL);
223 scale(_h_T2,1./ *_wsum);
224 scale(_h_T3,1./ *_wsum);
225 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h_cThetaL);
226 pair<double,pair<double,double> > R;
227 double tau = sqr(sqrtS()/(2*1.115683));
228 R.first = sqrt(tau*(1-alpha0.first)/(1+alpha0.first));
229 R.second.first = R.second.second = R.first/(1.-sqr(alpha0.first));
230
231 pair<double,double> c_T2 = calcCoeff(2,_h_T2);
232 pair<double,double> c_T3 = calcCoeff(3,_h_T3);
233
234 double sDelta = (-2.*(3. + alpha0.first)*c_T3.first)/(-aLambda*sqrt(1 - sqr(alpha0.first)));
235 double cDelta = (-3*(3 + alpha0.first)*c_T2.first)/(-sqr(aLambda)*sqrt(1 - sqr(alpha0.first)));
236
237 pair<double,pair<double,double> > Delta;
238 Delta.first = asin(sDelta);
239 if(cDelta<0.) Delta.first = M_PI-Delta.first;
240 Delta.second.first = (-4*(sqr(c_T3.second)*sqr(1. + alpha0.first)*
241 sqr(1 + alpha0.first)*sqr(3.+alpha0.first) +
242 sqr(alpha0.second.first)*sqr(c_T3.first)*sqr(1 + 3*alpha0.first)))
243 /(sqr(1.-sqr(alpha0.first))*(4*sqr(c_T3.first)*sqr(3 + alpha0.first) + sqr(aLambda)*
244 (-1 + sqr(alpha0.first))));
245 Delta.second.second = (-4*(sqr(c_T3.second)*sqr(1. + alpha0.first)*
246 sqr(1 + alpha0.first)*sqr(3.+alpha0.first) +
247 sqr(alpha0.second.second)*sqr(c_T3.first)*sqr(1 + 3*alpha0.first)))
248 /(sqr(1.-sqr(alpha0.first))*(4*sqr(c_T3.first)*sqr(3 + alpha0.first) + sqr(aLambda)*
249 (-1 + sqr(alpha0.first))));
250 Delta.first *=180./M_PI;
251 Delta.second.first *=180./M_PI;
252 Delta.second.second *=180./M_PI;
253 for(unsigned int iy=0;iy<3;++iy) {
254 double val;
255 pair<double,double> err;
256 if(iy==0) {
257 val = alpha0.first;
258 err = alpha0.second;
259 }
260 else if(iy==1) {
261 val = Delta.first ;
262 err = Delta.second;
263 }
264 else {
265 val = R.first;
266 err = R.second;
267 }
268 Estimate1DPtr tmp1;
269 book(tmp1,1,1,1+iy);
270 tmp1->bin(1).set(val,err);
271 BinnedEstimatePtr<string> tmp2;
272 book(tmp2,2,1,1+iy);
273 tmp2->binAt(_ecms).set(val,err);
274 }
275 }
276
277 /// @}
278
279
280 /// @name Histograms
281 /// @{
282 Histo1DPtr _h_T2,_h_T3;
283 Histo1DPtr _h_cThetaL;
284 CounterPtr _wsum;
285 string _ecms;
286 /// @}
287
288
289 };
290
291
292 RIVET_DECLARE_PLUGIN(BESIII_2023_I2637702);
293
294}
|