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