Rivet analyses referenceBESIII_2020_I1791570Analysis of $J/\psi$, $\psi(2S)$ decays to $\Sigma^+\bar\Sigma^-$Experiment: BESIII (BEPC) Inspire ID: 1791570 Status: VALIDATED Authors:
Beam energies: (1.6, 1.6); (1.8, 1.8) GeV Run details:
Analysis of the angular distribution of the baryons, and decay products, produced in $e^+e^-\to J/\psi, \psi(2S) \to \Sigma^+\bar\Sigma^-$. Gives information about the decay and is useful for testing correlations in hadron decays. N.B. The moment data is not corrected for efficiency/acceptance and should therefore only be used qualatively. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: BESIII_2020_I1791570.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 J/Psi, psi(2S) -> Sigma+ Sigmabar-
11 class BESIII_2020_I1791570 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2020_I1791570);
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_T1, "/TMP/T1",20,-1.,1.);
31 book(_h_T2, "/TMP/T2",20,-1.,1.);
32 book(_h_T3, "/TMP/T3",20,-1.,1.);
33 book(_h_T4, "/TMP/T4",20,-1.,1.);
34 book(_h_T5, "/TMP/T5",20,-1.,1.);
35
36 book(_h_cThetaL,"/TMP/cThetaL",20,-1.,1.);
37 if(isCompatibleWithSqrtS(3.1*GeV,1e-2))
38 book(_h_mu,1,1,1);
39 else if(isCompatibleWithSqrtS(3.686*GeV,1e-2))
40 book(_h_mu,1,1,2);
41 else
42 throw Error("Unexpected sqrtS ! Only 3.1 and 3.686 GeV atr supported");
43 book(_wsum,"/TMP/wsum");
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 /// Perform the per-event analysis
58 void analyze(const Event& event) {
59 // get the axis, direction of incoming electron
60 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
61 Vector3 axis;
62 if(beams.first.pid()>0)
63 axis = beams.first .momentum().p3().unit();
64 else
65 axis = beams.second.momentum().p3().unit();
66 // types of final state particles
67 const FinalState& fs = apply<FinalState>(event, "FS");
68 map<long,int> nCount;
69 int ntotal(0);
70 for (const Particle& p : fs.particles()) {
71 nCount[p.pid()] += 1;
72 ++ntotal;
73 }
74 // loop over Sigma+ baryons
75 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
76 Particle Sigma,SigBar;
77 bool matched(false);
78 for (const Particle& p : ufs.particles(Cuts::abspid==3222)) {
79 if(p.children().empty()) continue;
80 map<long,int> nRes=nCount;
81 int ncount = ntotal;
82 findChildren(p,nRes,ncount);
83 matched=false;
84 // check for antiparticle
85 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
86 if(p2.children().empty()) continue;
87 map<long,int> nRes2=nRes;
88 int ncount2 = ncount;
89 findChildren(p2,nRes2,ncount2);
90 if(ncount2==0) {
91 matched = true;
92 for(auto const & val : nRes2) {
93 if(val.second!=0) {
94 matched = false;
95 break;
96 }
97 }
98 // fond baryon and antibaryon
99 if(matched) {
100 if(p.pid()>0) {
101 Sigma = p;
102 SigBar = p2;
103 }
104 else {
105 Sigma = p2;
106 SigBar = p;
107 }
108 break;
109 }
110 }
111 }
112 if(matched) break;
113 }
114 if(!matched) vetoEvent;
115 // find proton
116 Particle proton;
117 matched = false;
118 for (const Particle & p : Sigma.children()) {
119 if(p.pid()==2212) {
120 matched=true;
121 proton=p;
122 }
123 else if(p.pid()!=111) {
124 matched = false;
125 break;
126 }
127 }
128 if(!matched) vetoEvent;
129 // find antiproton
130 Particle pbar;
131 matched = false;
132 for (const Particle & p : SigBar.children()) {
133 if(p.pid()==-2212) {
134 matched=true;
135 pbar=p;
136 }
137 else if(p.pid()!=111) {
138 matched = false;
139 break;
140 }
141 }
142 if(!matched) vetoEvent;
143 // boost to the Sigma rest frame
144 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Sigma.momentum().betaVec());
145 Vector3 e1z = Sigma.momentum().p3().unit();
146 Vector3 e1y = e1z.cross(axis).unit();
147 Vector3 e1x = e1y.cross(e1z).unit();
148 Vector3 axis1 = boost1.transform(proton.momentum()).p3().unit();
149 double n1x(e1x.dot(axis1)),n1y(e1y.dot(axis1)),n1z(e1z.dot(axis1));
150 // boost to the Sigma bar
151 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(SigBar.momentum().betaVec());
152 Vector3 axis2 = boost2.transform(pbar.momentum()).p3().unit();
153 double n2x(e1x.dot(axis2)),n2y(e1y.dot(axis2)),n2z(e1z.dot(axis2));
154 double cosL = axis.dot(Sigma.momentum().p3().unit());
155 double sinL = sqrt(1.-sqr(cosL));
156 double T1 = sqr(sinL)*n1x*n2x+sqr(cosL)*n1z*n2z;
157 double T2 = -sinL*cosL*(n1x*n2z+n1z*n2x);
158 double T3 = -sinL*cosL*n1y;
159 double T4 = -sinL*cosL*n2y;
160 double T5 = n1z*n2z-sqr(sinL)*n1y*n2y;
161 double mu = -(n1y-n2y);
162 _h_T1->fill(cosL,T1);
163 _h_T2->fill(cosL,T2);
164 _h_T3->fill(cosL,T3);
165 _h_T4->fill(cosL,T4);
166 _h_T5->fill(cosL,T5);
167 _h_mu->fill(cosL,mu);
168 _wsum->fill();
169 _h_cThetaL->fill(cosL);
170 }
171
172
173 pair<double,pair<double,double> > calcAlpha0(Histo1DPtr hist) {
174 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
175 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
176 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
177 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
178 for (const auto& bin : hist->bins() ) {
179 double Oi = bin.sumW();
180 if(Oi==0.) continue;
181 double a = d*(bin.xMax() - bin.xMin());
182 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
183 double Ei = bin.errW();
184 sum1 += a*Oi/sqr(Ei);
185 sum2 += b*Oi/sqr(Ei);
186 sum3 += sqr(a)/sqr(Ei);
187 sum4 += sqr(b)/sqr(Ei);
188 sum5 += a*b/sqr(Ei);
189 }
190 // calculate alpha
191 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
192 // and error
193 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
194 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
195 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
196 double dis = sqr(bb)-4.*aa*cc;
197 if(dis>0.) {
198 dis = sqrt(dis);
199 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
200 }
201 else {
202 return make_pair(alpha,make_pair(0.,0.));
203 }
204 }
205
206 pair<double,double> calcCoeff(unsigned int imode,Histo1DPtr hist) {
207 if(hist->numEntries()==0.) return make_pair(0.,0.);
208 double sum1(0.),sum2(0.);
209 for (const auto& bin : hist->bins() ) {
210 double Oi = bin.sumW();
211 if(Oi==0.) continue;
212 double ai(0.),bi(0.);
213 if(imode==0) {
214 bi = (pow(1.-sqr(bin.xMin()),1.5) - pow(1.-sqr(bin.xMax()),1.5))/3.;
215 }
216 else if(imode>=2 && imode<=4) {
217 bi = ( pow(bin.xMin(),3)*( -5. + 3.*sqr(bin.xMin())) +
218 pow(bin.xMax(),3)*( 5. - 3.*sqr(bin.xMax())))/15.;
219 }
220 else
221 assert(false);
222 double Ei = bin.errW();
223 sum1 += sqr(bi/Ei);
224 sum2 += bi/sqr(Ei)*(Oi-ai);
225 }
226 return make_pair(sum2/sum1,sqrt(1./sum1));
227 }
228
229 /// Normalise histograms etc., after the run
230 void finalize() {
231 normalize(_h_cThetaL);
232 scale(_h_T1, 1./ *_wsum);
233 scale(_h_T2, 1./ *_wsum);
234 scale(_h_T3, 1./ *_wsum);
235 scale(_h_T4, 1./ *_wsum);
236 scale(_h_T5, 1./ *_wsum);
237 scale(_h_mu, 2./ *_wsum);
238 // histos for J/psi or psi(2s)
239 int ih = isCompatibleWithSqrtS(3.1*GeV,1e-2) ? 1 : 2;
240 // calculate alpha0
241 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h_cThetaL);
242 Estimate1DPtr _h_alpha0;
243 book(_h_alpha0,4,1,ih);
244 _h_alpha0->bin(1).set(alpha0.first, alpha0.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 = calcCoeff(2,_h_T2);
250 pair<double,double> c_T3 = calcCoeff(3,_h_T3);
251 pair<double,double> c_T4 = calcCoeff(4,_h_T4);
252 double s4 = sqr(c_T2.first);
253 double s5 = sqr(c_T3.first);
254 double s6 = sqr(c_T4.first);
255 double disc = s1*s5*s6*(-9.*s2*s4 + 4.*s1*s5*s6);
256 if(disc>=0.) {
257 disc = sqrt(disc);
258 double aM = -sqrt(-1./s2/s6*(2.*s1*s5*s6+disc));
259 double aP = c_T4.first/c_T3.first*aM;
260 double aM_P = (2*(alpha0.first*c_T4.first*alpha0.second.first + c_T4.second*s2)*(disc + 2*s1*s5*s6)
261 - c_T4.first*s2*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.first +s3*c_T4.first*c_T3.second +s3*c_T3.first*c_T4.second) +
262 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
263 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.first - s2*s3*c_T4.first*c_T3.second
264 - s2*s3*c_T3.first*c_T4.second)* s4
265 + 8*(c_T3.first*c_T4.first*alpha0.second.first + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))
266 /(4*pow(3 + alpha0.first,3)*pow(c_T3.first,3)*pow(c_T4.first,3) -9*s2*s3*c_T3.first*c_T4.first*s4)))/
267 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
268 double aM_M = (2*(alpha0.first*c_T4.first*alpha0.second.second + c_T4.second*s2)*(disc + 2*s1*s5*s6)
269 - c_T4.first*s2*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.second +s3*c_T4.first*c_T3.second +s3*c_T3.first*c_T4.second) +
270 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
271 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.second - s2*s3*c_T4.first*c_T3.second
272 - s2*s3*c_T3.first*c_T4.second)* s4
273 + 8*(c_T3.first*c_T4.first*alpha0.second.second + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))
274 /(4*pow(3 + alpha0.first,3)*pow(c_T3.first,3)*pow(c_T4.first,3) -9*s2*s3*c_T3.first*c_T4.first*s4)))/
275 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
276 double aP_M = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
277 (-2*c_T3.second - (2*alpha0.first*c_T3.first*alpha0.second.first)/s2 + (c_T3.first*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.first + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)
278 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
279 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.first - s2*s3*c_T4.first*c_T3.second
280 - s2*s3*c_T3.first*c_T4.second)* s4 +
281 8*(c_T3.first*c_T4.first*alpha0.second.first + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))/
282 (4* pow(3 + alpha0.first,3)* pow(c_T3.first,3)* pow(c_T4.first,3) - 9*s2*s3*c_T3.first*c_T4.first*s4)))/
283 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
284 double aP_P = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
285 (-2*c_T3.second - (2*alpha0.first*c_T3.first*alpha0.second.second)/s2 + (c_T3.first*(4*s3*c_T3.first*c_T4.first*(c_T3.first*c_T4.first*alpha0.second.second + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)
286 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
287 + 9*((1 - alpha0.first*(3 + 2*alpha0.first))* c_T3.first*c_T4.first*alpha0.second.second - s2*s3*c_T4.first*c_T3.second
288 - s2*s3*c_T3.first*c_T4.second)* s4 +
289 8*(c_T3.first*c_T4.first*alpha0.second.second + s3*c_T4.first*c_T3.second + s3*c_T3.first*c_T4.second)* s1*s5*s6))/
290 (4* pow(3 + alpha0.first,3)* pow(c_T3.first,3)* pow(c_T4.first,3) - 9*s2*s3*c_T3.first*c_T4.first*s4)))/
291 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
292 Estimate1DPtr _h_alphaM;
293 book(_h_alphaM,2,1,1);
294 _h_alphaM->bin(1).set(aM, make_pair(-aM_M , -aM_P ));
295 Estimate1DPtr _h_alphaP;
296 book(_h_alphaP,2,1,2);
297 _h_alphaP->bin(1).set(aP, make_pair(-aP_M , -aP_P ));
298 Estimate1DPtr _h_alphabar;
299 book(_h_alphabar,2,1,3);
300 _h_alphabar->bin(1).set(0.5*(aM-aP),
301 make_pair(0.5*sqrt(sqr(aM_M)+sqr(aP_P)) ,
302 0.5*sqrt(sqr(aM_P)+sqr(aP_M))));
303 // now for Delta
304 double sDelta = (-2.*(3. + alpha0.first)*c_T3.first)/(aM*sqrt(1 - sqr(alpha0.first)));
305 double cDelta = (-3*(3 + alpha0.first)*c_T2.first)/(aM*aP*sqrt(1 - sqr(alpha0.first)));
306
307 double Delta = asin(sDelta);
308 if(cDelta<0.) Delta = M_PI-Delta;
309 double ds_P = (-9*c_T2.first*((-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T3.first*c_T4.first*c_T2.second + c_T2.first*c_T4.first*(c_T3.first*(alpha0.second.first + 3*alpha0.first*alpha0.second.first) -(-1 + alpha0.first)*(1 + alpha0.first)*(3 + alpha0.first)*c_T3.second)
310 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2.first*c_T3.first*c_T4.second)*disc)/
311 (pow(1 - pow(alpha0.first,2),1.5)*pow(c_T4.first,3)*pow(-((disc + 2*s1*s5*s6)/ (s2*s6)),1.5)*(-9*s2*s4 + 4*s1*s5*s6));
312 double ds_M = (-9*c_T2.first*((-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T3.first*c_T4.first*c_T2.second + c_T2.first*c_T4.first*(c_T3.first*(alpha0.second.second + 3*alpha0.first*alpha0.second.second) -(-1 + alpha0.first)*(1 + alpha0.first)*(3 + alpha0.first)*c_T3.second)
313 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2.first*c_T3.first*c_T4.second)*disc)/
314 (pow(1 - pow(alpha0.first,2),1.5)*pow(c_T4.first,3)*pow(-((disc + 2*s1*s5*s6)/ (s2*s6)),1.5)*(-9*s2*s4 + 4*s1*s5*s6));
315 ds_P /= sqrt(1.-sqr(sDelta));
316 ds_M /= sqrt(1.-sqr(sDelta));
317 Estimate1DPtr _h_sin;
318 book(_h_sin,3,1,ih);
319 _h_sin->bin(1).set(Delta/M_PI*180., make_pair( -ds_P/M_PI*180., -ds_M/M_PI*180. ));
320 }
321 }
322
323 /// @}
324
325
326 /// @name Histograms
327 /// @{
328 Histo1DPtr _h_T1,_h_T2,_h_T3,_h_T4,_h_T5;
329 Histo1DPtr _h_cThetaL;
330 Histo1DPtr _h_mu;
331 CounterPtr _wsum;
332 /// @}
333
334 };
335
336
337 RIVET_DECLARE_PLUGIN(BESIII_2020_I1791570);
338
339}
|