Rivet analyses referenceBESIII_2023_I2636760Decay assymetry in $\Sigma^+\to p\gamma$Experiment: BESIII () Inspire ID: 2636760 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.6, 1.6) GeV Run details:
Analysis of the angular distribution of the baryons, and decay products, produced in $e^+e^-\to J/\psi \to \Sigma^+\bar\Sigma^-$. Gives information about the decay and is useful for testing correlations in hadron decays. Source code: BESIII_2023_I2636760.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 Sigma + -> p gamma
11 class BESIII_2023_I2636760 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2636760);
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 for (unsigned int ix=0; ix<2; ++ix) {
31 book(_n[ix],"TMP/n_" + toString(ix+1));
32 book(_t[ix],"TMP/t_" + toString(ix+1));
33 }
34 book(_n[2],"TMP/n_3");
35 book(_t[2],"TMP/t_3");
36 }
37
38
39 void findChildren(const Particle& p,map<long,int>& nRes, int& ncount) {
40 for (const Particle& child : p.children()) {
41 if (child.children().empty()) {
42 nRes[child.pid()]-=1;
43 --ncount;
44 }
45 else {
46 findChildren(child,nRes,ncount);
47 }
48 }
49 }
50
51 /// Perform the per-event analysis
52 void analyze(const Event& event) {
53 // get the axis, direction of incoming electron
54 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
55 Vector3 axis;
56 if (beams.first.pid()>0) {
57 axis = beams.first .mom().p3().unit();
58 }
59 else {
60 axis = beams.second.mom().p3().unit();
61 }
62 // types of final state particles
63 const FinalState& fs = apply<FinalState>(event, "FS");
64 map<long,int> nCount;
65 int ntotal(0);
66 for (const Particle& p : fs.particles()) {
67 nCount[p.pid()] += 1;
68 ++ntotal;
69 }
70 // loop over Sigma+ baryons
71 const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
72 Particle Sigma, SigBar;
73 bool matched(false);
74 for (const Particle& p : ufs.particles(Cuts::abspid==3222)) {
75 if (p.children().empty()) continue;
76 map<long,int> nRes=nCount;
77 int ncount = ntotal;
78 findChildren(p, nRes, ncount);
79 matched=false;
80 // check for antiparticle
81 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
82 if(p2.children().empty()) continue;
83 map<long,int> nRes2=nRes;
84 int ncount2 = ncount;
85 findChildren(p2, nRes2, ncount2);
86 if (ncount2==0) {
87 matched = true;
88 for (const auto& val : nRes2) {
89 if (val.second!=0) {
90 matched = false;
91 break;
92 }
93 }
94 // fond baryon and antibaryon
95 if (matched) {
96 if (p.pid()>0) {
97 Sigma = p;
98 SigBar = p2;
99 }
100 else {
101 Sigma = p2;
102 SigBar = p;
103 }
104 break;
105 }
106 }
107 }
108 if (matched) break;
109 }
110 if (!matched) vetoEvent;
111 if (Sigma .children().size()!=2) vetoEvent;
112 if (SigBar.children().size()!=2) vetoEvent;
113 // find proton
114 bool radiative[2]={false,false};
115 // identify Sigma decay
116 Particle baryon1;
117 if (Sigma.children()[0].pid()==PID::PROTON && Sigma.children()[1].pid()==PID::PI0) {
118 radiative[0]=false;
119 baryon1 = Sigma.children()[0];
120 }
121 else if (Sigma.children()[1].pid()==PID::PROTON && Sigma.children()[0].pid()==PID::PI0) {
122 radiative[0]=false;
123 baryon1 = Sigma.children()[1];
124 }
125 else if (Sigma.children()[0].pid()==PID::PROTON && Sigma.children()[1].pid()==PID::PHOTON) {
126 radiative[0]=true;
127 baryon1 = Sigma.children()[0];
128 }
129 else if (Sigma.children()[1].pid()==PID::PROTON && Sigma.children()[0].pid()==PID::PHOTON ) {
130 radiative[0]=true;
131 baryon1 = Sigma.children()[1];
132 }
133 else {
134 vetoEvent;
135 }
136 // antisigma decay
137 Particle baryon2;
138 if (SigBar.children()[0].pid()==PID::ANTIPROTON && SigBar.children()[1].pid()==PID::PI0 ) {
139 radiative[1]=false;
140 baryon2 = SigBar.children()[0];
141 }
142 else if (SigBar.children()[1].pid()==PID::ANTIPROTON && SigBar.children()[0].pid()==PID::PI0 ) {
143 radiative[1]=false;
144 baryon2 = SigBar.children()[1];
145 }
146 else if (SigBar.children()[0].pid()==PID::ANTIPROTON && SigBar.children()[1].pid()==PID::PHOTON ) {
147 radiative[1]=true;
148 baryon2 = SigBar.children()[0];
149 }
150 else if (SigBar.children()[1].pid()==PID::ANTIPROTON && SigBar.children()[0].pid()==PID::PHOTON ) {
151 radiative[1]=true;
152 baryon2 = SigBar.children()[1];
153 }
154 else {
155 vetoEvent;
156 }
157 if (radiative[0] == radiative[1]) vetoEvent;
158 // boost to the Sigma rest frame
159 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Sigma.mom().betaVec());
160 Vector3 e1z = Sigma.mom().p3().unit();
161 Vector3 e1y = e1z.cross(axis).unit();
162 Vector3 e1x = e1y.cross(e1z).unit();
163 Vector3 axis1 = boost1.transform(baryon1.mom()).p3().unit();
164 double n1x(e1x.dot(axis1)), n1z(e1z.dot(axis1));
165 // boost to the Sigma bar
166 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(SigBar.mom().betaVec());
167 Vector3 axis2 = boost2.transform(baryon2.mom()).p3().unit();
168 double n2x(e1x.dot(axis2)),n2z(e1z.dot(axis2));
169 double cosL = axis.dot(Sigma.mom().p3().unit());
170 double sinL = sqrt(1.-sqr(cosL));
171 double T1 = sqr(sinL)*n1x*n2x+sqr(cosL)*n1z*n2z;
172 // sigma+ -> p gamma
173 if (radiative[0]) {
174 _n[0]->fill();
175 _n[2]->fill();
176 _t[0]->fill(T1);
177 _t[2]->fill(T1);
178 }
179 // sigmabar- -> pbar gamma
180 else {
181 _n[1]->fill();
182 _n[2]->fill();
183 _t[1]->fill(T1);
184 _t[2]->fill(T1);
185 }
186 }
187
188 /// Normalise histograms etc., after the run
189 void finalize() {
190 // values of constants
191 const double aPsi =-0.508;
192 const double aPlus =-0.998;
193 const double factor = 45.*(3. +aPsi)/(11. + 5.*aPsi)/aPlus;
194 // alpha from the moments
195 for (unsigned int ix=0;ix<3;++ix) {
196 double value = _t[ix]->val()/_n[ix]->val();
197 double error = _t[ix]->err()/_n[ix]->val();
198 value *= factor;
199 error *= abs(factor);
200 if (ix==1) value *=-1.;
201 Estimate0DPtr alpha;
202 book(alpha, 1, 1, 1+ix);
203 alpha->set(value, error);
204 }
205
206 }
207
208 /// @}
209
210
211 /// @name Histograms
212 /// @{
213 CounterPtr _n[3];
214 Histo1DPtr _h_ctheta[3];
215 CounterPtr _t[3];
216 /// @}
217
218
219 };
220
221
222 RIVET_DECLARE_PLUGIN(BESIII_2023_I2636760);
223
224}
|