Rivet analyses referenceBESIII_2021_I1974025Measurement of e+e−→Λ0ˉΛ0 at 3.773 GeVExperiment: BESIII (BEPC) Inspire ID: 1974025 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.9, 1.9) GeV Run details:
Measurement of the angular distribution and polarization for e+e−→Λ0ˉΛ0 at 3.773 GeV by BESIII. Source code: BESIII_2021_I1974025.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- > Lambda, Lambdabar
11 class BESIII_2021_I1974025 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2021_I1974025);
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(_wsum,"TMP/wsum");
29 // for(unsigned int ix=0;ix<6;++ix)
30 // book(_h_F[ix],1,1,1+ix);
31 for(unsigned int ix=0;ix<6;++ix)
32 book(_h_F[ix],"TMP/F_"+toString(ix+1),20,-1.,1.);
33 book(_h_F[5],1,1,6);
34 book(_h_mu,2,1,1);
35 }
36
37 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
38 for (const Particle &child : p.children()) {
39 if(child.children().empty()) {
40 nRes[child.pid()]-=1;
41 --ncount;
42 }
43 else
44 findChildren(child,nRes,ncount);
45 }
46 }
47
48
49 /// Perform the per-event analysis
50 void analyze(const Event& event) {
51 // get the axis, direction of incoming electron
52 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
53 Vector3 axis;
54 if(beams.first.pid()>0)
55 axis = beams.first .momentum().p3().unit();
56 else
57 axis = beams.second.momentum().p3().unit();
58 const FinalState& fs = apply<FinalState>(event, "FS");
59 // total hadronic and muonic cross sections
60 map<long,int> nCount;
61 int ntotal(0);
62 for (const Particle& p : fs.particles()) {
63 nCount[p.pid()] += 1;
64 ++ntotal;
65 }
66 // find the Lambdas
67 bool matched = false;
68 const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
69 Particle Lambda,LamBar;
70 for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
71 const Particle& p1 = ufs.particles()[ix];
72 if(abs(p1.pid())!=3122) continue;
73 // check fs
74 bool fs = true;
75 for (const Particle & child : p1.children()) {
76 if(child.pid()==p1.pid()) {
77 fs = false;
78 break;
79 }
80 }
81 if(!fs) continue;
82 // find the children
83 map<long,int> nRes = nCount;
84 int ncount = ntotal;
85 findChildren(p1,nRes,ncount);
86 for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
87 matched=false;
88 const Particle& p2 = ufs.particles()[iy];
89 if(abs(p2.pid())!=3122) continue;
90 // check fs
91 bool fs = true;
92 for (const Particle & child : p2.children()) {
93 if(child.pid()==p2.pid()) {
94 fs = false;
95 break;
96 }
97 }
98 if(!fs) continue;
99 map<long,int> nRes2 = nRes;
100 int ncount2 = ncount;
101 findChildren(p2,nRes2,ncount2);
102 if(ncount2!=0) continue;
103 matched=true;
104 for(auto const & val : nRes2) {
105 if(val.second!=0) {
106 matched = false;
107 break;
108 }
109 }
110 if(matched) {
111 if(p1.pid()==PID::LAMBDA) {
112 Lambda=p1;
113 LamBar=p2;
114 }
115 else {
116 Lambda=p2;
117 LamBar=p1;
118 }
119 break;
120 }
121 }
122 if(matched) break;
123 }
124 // and the children
125 Particle proton;
126 matched = false;
127 for (const Particle & p : Lambda.children()) {
128 if(p.pid()==2212) {
129 matched=true;
130 proton=p;
131 }
132 else if(p.pid()==PID::PHOTON)
133 vetoEvent;
134 }
135 if(!matched) vetoEvent;
136 Particle baryon;
137 matched = false;
138 for (const Particle & p : LamBar.children()) {
139 if(p.pid()==-2212) {
140 baryon=p;
141 matched=true;
142 }
143 else if(p.pid()==PID::PHOTON)
144 vetoEvent;
145 }
146 if(!matched) vetoEvent;
147 // now for the polarization measurements
148 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Lambda.momentum().betaVec());
149 Vector3 e1z = Lambda.momentum().p3().unit();
150 Vector3 e1y = e1z.cross(axis).unit();
151 Vector3 e1x = e1y.cross(e1z).unit();
152 Vector3 axis1 = boost1.transform(proton.momentum()).p3().unit();
153 double n1x(e1x.dot(axis1)),n1y(e1y.dot(axis1)),n1z(e1z.dot(axis1));
154 // boost to the Lambda bar
155 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(LamBar.momentum().betaVec());
156 Vector3 axis2 = boost2.transform(baryon.momentum()).p3().unit();
157 double n2x(e1x.dot(axis2)),n2y(e1y.dot(axis2)),n2z(e1z.dot(axis2));
158 double cosL = -axis.dot(Lambda.momentum().p3().unit());
159 double sinL = sqrt(1.-sqr(cosL));
160 double T1 = sqr(sinL)*n1x*n2x+sqr(cosL)*n1z*n2z;
161 double T2 = sinL*cosL*(n1x*n2z+n1z*n2x);
162 double T3 = sinL*cosL*n1y;
163 double T4 = sinL*cosL*n2y;
164 double T5 = n1z*n2z-sqr(sinL)*n1y*n2y;
165 double mu = n1y-n2y;
166 _h_F[0]->fill(cosL,T1);
167 _h_F[1]->fill(cosL,T2);
168 _h_F[2]->fill(cosL,T3);
169 _h_F[3]->fill(cosL,T4);
170 _h_F[4]->fill(cosL,T5);
171 _h_F[5]->fill(cosL);
172 _h_mu->fill(cosL,mu);
173 _wsum->fill();
174 }
175
176 pair<double,pair<double,double> > calcAlpha0(Histo1DPtr hist) {
177 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
178 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
179 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
180 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
181 for (const auto& bin : hist->bins() ) {
182 double Oi = bin.sumW();
183 if (Oi==0.) continue;
184 double a = d*(bin.xMax() - bin.xMin());
185 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
186 double Ei = bin.errW();
187 sum1 += a*Oi/sqr(Ei);
188 sum2 += b*Oi/sqr(Ei);
189 sum3 += sqr(a)/sqr(Ei);
190 sum4 += sqr(b)/sqr(Ei);
191 sum5 += a*b/sqr(Ei);
192 }
193 // calculate alpha
194 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
195 // and error
196 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
197 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
198 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
199 double dis = sqr(bb)-4.*aa*cc;
200 if(dis>0.) {
201 dis = sqrt(dis);
202 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
203 }
204 else {
205 return make_pair(alpha,make_pair(0.,0.));
206 }
207 }
208
209 pair<double,double> calcCoeff(unsigned int imode,Histo1DPtr hist) {
210 if(hist->numEntries()==0.) return make_pair(0.,0.);
211 double sum1(0.),sum2(0.);
212 for (const auto& bin : hist->bins() ) {
213 double Oi = bin.sumW();
214 if(Oi==0.) continue;
215 double ai(0.),bi(0.);
216 if(imode==0) {
217 bi = (pow(1.-sqr(bin.xMin()),1.5) - pow(1.-sqr(bin.xMax()),1.5))/3.;
218 }
219 else if(imode>=2 && imode<=4) {
220 bi = ( pow(bin.xMin(),3)*( -5. + 3.*sqr(bin.xMin())) +
221 pow(bin.xMax(),3)*( 5. - 3.*sqr(bin.xMax())))/15.;
222 }
223 else
224 assert(false);
225 double Ei = bin.errW();
226 sum1 += sqr(bi/Ei);
227 sum2 += bi/sqr(Ei)*(Oi-ai);
228 }
229 return make_pair(sum2/sum1,sqrt(1./sum1));
230 }
231
232 /// Normalise histograms etc., after the run
233 void finalize() {
234 // normalize histograms
235 for(unsigned int ix=0;ix<6;++ix)
236 scale(_h_F[ix], 1./ *_wsum);
237 scale(_h_mu, 10./ *_wsum);
238 // value of aLambda assumed in paper
239 double aLambda = 0.754;
240 // calculate alpha0
241 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h_F[5]);
242 Estimate0DPtr _h_alpha0;
243 book(_h_alpha0,3,1,1);
244 _h_alpha0->set(alpha0.first, make_pair(-alpha0.second.first,alpha0.second.second));
245 double s2 = -1. + sqr(alpha0.first);
246 double s3 = 3 + alpha0.first;
247 double s1 = sqr(s3);
248 // alpha- and alpha+ from proton data
249 pair<double,double> c_T2_p = calcCoeff(2,_h_F[1]);
250 pair<double,double> c_T3_p = calcCoeff(3,_h_F[2]);
251 pair<double,double> c_T4_p = calcCoeff(4,_h_F[3]);
252 double s4 = sqr(c_T2_p.first);
253 double s5 = sqr(c_T3_p.first);
254 double s6 = sqr(c_T4_p.first);
255 double disc = s1*s5*s6*(-9.*s2*s4 + 4.*s1*s5*s6);
256 // now for Delta
257 if(disc>0) {
258 double sDelta = (-2.*(3. + alpha0.first)*c_T3_p.first)/(aLambda*sqrt(1 - sqr(alpha0.first)));
259 double cDelta = (-3*(3 + alpha0.first)*c_T2_p.first)/(-aLambda*aLambda*sqrt(1 - sqr(alpha0.first)));
260 double Delta = asin(sDelta);
261 if(cDelta<0.) Delta = M_PI-Delta;
262 double ds_P = (-9*c_T2_p.first*((-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T3_p.first*c_T4_p.first*c_T2_p.second + c_T2_p.first*c_T4_p.first*(c_T3_p.first*(alpha0.second.first + 3*alpha0.first*alpha0.second.first) -(-1 + alpha0.first)*(1 + alpha0.first)*(3 + alpha0.first)*c_T3_p.second)
263 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2_p.first*c_T3_p.first*c_T4_p.second)*disc)/
264 (pow(1 - pow(alpha0.first,2),1.5)*pow(c_T4_p.first,3)*pow(-((disc + 2*s1*s5*s6)/ (s2*s6)),1.5)*(-9*s2*s4 + 4*s1*s5*s6));
265 double ds_M = (-9*c_T2_p.first*((-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T3_p.first*c_T4_p.first*c_T2_p.second + c_T2_p.first*c_T4_p.first*(c_T3_p.first*(alpha0.second.second + 3*alpha0.first*alpha0.second.second) -(-1 + alpha0.first)*(1 + alpha0.first)*(3 + alpha0.first)*c_T3_p.second)
266 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2_p.first*c_T3_p.first*c_T4_p.second)*disc)/
267 (pow(1 - pow(alpha0.first,2),1.5)*pow(c_T4_p.first,3)*pow(-((disc + 2*s1*s5*s6)/ (s2*s6)),1.5)*(-9*s2*s4 + 4*s1*s5*s6));
268 ds_P /= sqrt(1.-sqr(sDelta));
269 ds_M /= sqrt(1.-sqr(sDelta));
270 Estimate0DPtr _h_sin;
271 book(_h_sin,3,1,2);
272 _h_sin->set(Delta/M_PI*180., make_pair( ds_M/M_PI*180., -ds_P/M_PI*180. ));
273 }
274 }
275
276 /// @}
277
278
279 /// @name Histograms
280 /// @{
281 Histo1DPtr _h_F[6],_h_mu;
282 CounterPtr _wsum;
283 /// @}
284
285
286 };
287
288
289 RIVET_DECLARE_PLUGIN(BESIII_2021_I1974025);
290
291}
|
