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,1e-2))
38 book(_h_mu,1,1,1);
39 else if(isCompatibleWithSqrtS(3.686,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 (auto bin : hist->bins() ) {
179 double Oi = bin.area();
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.areaErr();
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 (auto bin : hist->bins() ) {
210 double Oi = bin.area();
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.areaErr();
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,1e-2) ? 1 : 2;
240 // calculate alpha0
241 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h_cThetaL);
242 Scatter2DPtr _h_alpha0;
243 book(_h_alpha0,4,1,ih);
244 _h_alpha0->addPoint(0.5, alpha0.first, make_pair(0.5,0.5),
245 make_pair(alpha0.second.first,alpha0.second.second) );
246 double s2 = -1. + sqr(alpha0.first);
247 double s3 = 3 + alpha0.first;
248 double s1 = sqr(s3);
249 // alpha- and alpha+ from proton data
250 pair<double,double> c_T2 = calcCoeff(2,_h_T2);
251 pair<double,double> c_T3 = calcCoeff(3,_h_T3);
252 pair<double,double> c_T4 = calcCoeff(4,_h_T4);
253 double s4 = sqr(c_T2.first);
254 double s5 = sqr(c_T3.first);
255 double s6 = sqr(c_T4.first);
256 double disc = s1*s5*s6*(-9.*s2*s4 + 4.*s1*s5*s6);
257 if(disc>=0.) {
258 disc = sqrt(disc);
259 double aM = -sqrt(-1./s2/s6*(2.*s1*s5*s6+disc));
260 double aP = c_T4.first/c_T3.first*aM;
261 double aM_P = (2*(alpha0.first*c_T4.first*alpha0.second.first + c_T4.second*s2)*(disc + 2*s1*s5*s6)
262 - 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) +
263 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
264 + 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
265 - s2*s3*c_T3.first*c_T4.second)* s4
266 + 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))
267 /(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)))/
268 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
269 double aM_M = (2*(alpha0.first*c_T4.first*alpha0.second.second + c_T4.second*s2)*(disc + 2*s1*s5*s6)
270 - 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) +
271 (disc*(- 9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
272 + 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
273 - s2*s3*c_T3.first*c_T4.second)* s4
274 + 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))
275 /(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)))/
276 (2.*pow(c_T4.first,3)*pow(s2,2)*sqrt(-((disc + 2*s1*s5*s6)/(s2*s6))));
277 double aP_M = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
278 (-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)
279 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
280 + 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
281 - s2*s3*c_T3.first*c_T4.second)* s4 +
282 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))/
283 (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)))/
284 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
285 double aP_P = (c_T4.first*sqrt(-((disc + 2*s1*s5*s6)/ (s2*s6)))*
286 (-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)
287 + (disc*(-9*s2*s3*c_T2.first*c_T3.first*c_T4.first* c_T2.second
288 + 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
289 - s2*s3*c_T3.first*c_T4.second)* s4 +
290 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))/
291 (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)))/
292 (disc + 2*s1*s5*s6)))/(2.*pow(c_T3.first,2));
293 Scatter2DPtr _h_alphaM;
294 book(_h_alphaM,2,1,1);
295 _h_alphaM->addPoint(0.5, aM, make_pair(0.5,0.5),
296 make_pair(-aM_M , -aM_P ) );
297 Scatter2DPtr _h_alphaP;
298 book(_h_alphaP,2,1,2);
299 _h_alphaP->addPoint(0.5, aP, make_pair(0.5,0.5),
300 make_pair(-aP_M , -aP_P ) );
301 Scatter2DPtr _h_alphabar;
302 book(_h_alphabar,2,1,3);
303 _h_alphabar->addPoint(0.5, 0.5*(aM-aP), make_pair(0.5,0.5),
304 make_pair(0.5*sqrt(sqr(aM_M)+sqr(aP_P)) ,
305 0.5*sqrt(sqr(aM_P)+sqr(aP_M))) );
306 // now for Delta
307 double sDelta = (-2.*(3. + alpha0.first)*c_T3.first)/(aM*sqrt(1 - sqr(alpha0.first)));
308 double cDelta = (-3*(3 + alpha0.first)*c_T2.first)/(aM*aP*sqrt(1 - sqr(alpha0.first)));
309
310 double Delta = asin(sDelta);
311 if(cDelta<0.) Delta = M_PI-Delta;
312 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)
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 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)
316 - (-1 + alpha0.first)*(1 + alpha0.first)* (3 + alpha0.first)*c_T2.first*c_T3.first*c_T4.second)*disc)/
317 (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));
318 ds_P /= sqrt(1.-sqr(sDelta));
319 ds_M /= sqrt(1.-sqr(sDelta));
320 Scatter2DPtr _h_sin;
321 book(_h_sin,3,1,ih);
322 _h_sin->addPoint(0.5, Delta/M_PI*180., make_pair(0.5,0.5), make_pair( -ds_P/M_PI*180., -ds_M/M_PI*180. ) );
323 }
324 }
325
326 ///@}
327
328
329 /// @name Histograms
330 ///@{
331 Histo1DPtr _h_T1,_h_T2,_h_T3,_h_T4,_h_T5;
332 Histo1DPtr _h_cThetaL;
333 Histo1DPtr _h_mu;
334 CounterPtr _wsum;
335 ///@}
336
337
338 };
339
340
341 RIVET_DECLARE_PLUGIN(BESIII_2020_I1791570);
342
343}
|