Rivet analyses referenceLHCB_2024_I2824757Decay parameters in $\Lambda^0_b\to\Lambda_c^+(\pi^-,K^-)$ with $\Lambda_c^+\to \Lambda^0(\pi^+,K^+)$ or $\Lambda_c^+\to p K^0_S$Experiment: LHCB (LHC) Inspire ID: 2824757 Status: VALIDATED NOHEPDATA Authors:
Beam energies: ANY Run details:
Measurement of the decay parameters in in $\Lambda^0_b\to\Lambda_c^+(\pi^-,K^-)$ with $\Lambda_c^+\to \Lambda^0(\pi^+,K^+)$ or $\Lambda_c^+\to p K^0_S$, with $\Lambda^0\to p\pi^-$. The data was read from Tables 1 and 2 in the paper. Source code: LHCB_2024_I2824757.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4
5namespace Rivet {
6
7
8 /// @brief Lambda_b0 -> Lambda_c+ pi-,K-
9 class LHCB_2024_I2824757 : public Analysis {
10 public:
11
12 /// Constructor
13 RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2024_I2824757);
14
15
16 /// @name Analysis methods
17 /// @{
18
19 /// Book histograms and initialise projections before the run
20 void init() {
21 // Initialise and register projections
22 declare(UnstableParticles(Cuts::abspid==5122), "UFS" );
23 for(unsigned int ip=0;ip<2;++ip) {
24 for(unsigned int ix=0;ix<2;++ix) {
25 for(unsigned int iy=0;iy<3;++iy) {
26 string base = toString(ip)+"_"+toString(ix)+"_"+toString(iy);
27 book(_h_cos1 [ip][ix][iy],"h_cos_1_" +base,20,-1,1);
28 if(iy==2) continue;
29 book(_h_cos2 [ip][ix][iy],"h_cos_2_" +base,20,-1,1);
30 book(_p_cos12[ip][ix][iy],"p_cos_12_"+base,1,-1,1);
31 book(_h_phi2 [ip][ix][iy],"h_phi_2_" +base,20,-M_PI,M_PI);
32 }
33 }
34 }
35 }
36
37
38 /// Perform the per-event analysis
39 void analyze(const Event& event) {
40 // loop over Lambda_b0 baryons
41 for (const Particle& Lamb : apply<UnstableParticles>(event, "UFS").particles()) {
42 int sign = Lamb.pid()/5122;
43 if(Lamb.children().size()!=2) continue;
44 Particle baryon1,meson1;
45 if ( Lamb.children()[0].pid()==sign*4122 &&
46 (Lamb.children()[1].pid()==-sign*211 ||
47 Lamb.children()[1].pid()==-sign*321)) {
48 baryon1 = Lamb.children()[0];
49 meson1 = Lamb.children()[1];
50 }
51 else if ( Lamb.children()[1].pid()==sign*4122 &&
52 (Lamb.children()[0].pid()==-sign*211 ||
53 Lamb.children()[0].pid()==-sign*321)) {
54 baryon1 = Lamb.children()[1];
55 meson1 = Lamb.children()[0];
56 }
57 else continue;
58 if(baryon1.children().size()!=2) continue;
59 unsigned int ib = meson1.abspid()==211 ? 0 : 1;
60 Particle baryon2,meson2;
61 unsigned int ic=0;
62 if (baryon1.children()[0].pid()== sign*3122 && baryon1.children()[1].pid()== sign*211) {
63 baryon2 = baryon1.children()[0];
64 meson2 = baryon1.children()[1];
65 ic=0;
66 }
67 else if (baryon1.children()[1].pid()== sign*3122 && baryon1.children()[0].pid()== sign*211) {
68 baryon2 = baryon1.children()[1];
69 meson2 = baryon1.children()[0];
70 ic=0;
71 }
72 else if (baryon1.children()[0].pid()== sign*3122 && baryon1.children()[1].pid()== sign*321) {
73 baryon2 = baryon1.children()[0];
74 meson2 = baryon1.children()[1];
75 ic=1;
76 }
77 else if (baryon1.children()[1].pid()== sign*3122 && baryon1.children()[0].pid()== sign*321) {
78 baryon2 = baryon1.children()[1];
79 meson2 = baryon1.children()[0];
80 ic=1;
81 }
82 else if( baryon1.children()[0].pid()== sign*2212 &&
83 (baryon1.children()[1].pid()==-sign*311 ||
84 baryon1.children()[1].pid()== 310 )) {
85 baryon2 = baryon1.children()[0];
86 meson2 = baryon1.children()[1];
87 ic=2;
88 }
89 else if( baryon1.children()[1].pid()== sign*2212 &&
90 (baryon1.children()[1].pid()==-sign*311 ||
91 baryon1.children()[1].pid()== 310 )) {
92 baryon2 = baryon1.children()[1];
93 meson2 = baryon1.children()[0];
94 ic=2;
95 }
96 else continue;
97 // deal with Kbar0 in Lambda_c+ -> Kbar0 p
98 if (ic==2 && meson2.pid()!=310) {
99 if (meson2.children().size()==1) {
100 meson2=meson2.children()[0];
101 }
102 if (meson2.pid()!=310) continue;
103 }
104 // first boost to the Lambdab0 rest frame
105 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Lamb.mom().betaVec());
106 FourMomentum pbaryon1 = boost1.transform(baryon1.mom());
107 FourMomentum pbaryon2 = boost1.transform(baryon2.mom());
108 FourMomentum pmeson2 = boost1.transform(meson2.mom());
109 // boost to lambda_c rest frame
110 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pbaryon1.betaVec());
111 Vector3 axis1 = pbaryon1.p3().unit();
112 FourMomentum pp = boost2.transform(pbaryon2);
113 double cTheta1 = pp.p3().unit().dot(axis1);
114 _h_cos1[(1-sign)/2][ib][ic]->fill(cTheta1);
115 // that's it for Lambda_c+ -> p+ K_S0
116 if (ic==2) continue;
117 if (baryon2.children().size()!=2) continue;
118 Particle baryon3,meson3;
119 if (baryon2.children()[0].pid()== sign*2212 &&
120 baryon2.children()[1].pid()==-sign*211) {
121 baryon3 = baryon2.children()[0];
122 meson3 = baryon2.children()[1];
123 }
124 else if (baryon2.children()[1].pid()== sign*2212 &&
125 baryon2.children()[0].pid()==-sign*211) {
126 baryon3 = baryon2.children()[1];
127 meson3 = baryon2.children()[0];
128 }
129 else continue;
130 Vector3 axis2=pp.p3().unit();
131 LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pp.betaVec());
132 FourMomentum pbaryon3 = boost3.transform(boost2.transform(boost1.transform(baryon3.mom())));
133 double cTheta2 = pbaryon3.p3().unit().dot(axis2);
134 _h_cos2 [(1-sign)/2][ib][ic]->fill(cTheta2);
135 _p_cos12[(1-sign)/2][ib][ic]->fill(cTheta1*cTheta2,cTheta1*cTheta2);
136 Vector3 trans1 = axis1 - axis1.dot(axis2)*axis2;
137 Vector3 trans2 = pbaryon3.p3()-pbaryon3.p3().dot(axis2)*pbaryon3.p3();
138 double phi2 = atan2(trans1.cross(trans2).dot(axis2), trans1.dot(trans2));
139 _h_phi2[(1-sign)/2][ib][ic]->fill(phi2);
140 }
141 }
142
143 pair<double,double> calcAlpha(Histo1DPtr hist) const {
144 if (hist->numEntries()==0.) return make_pair(0.,0.);
145 double sum1(0.),sum2(0.);
146 for (const auto& bin : hist->bins()) {
147 double Oi = bin.sumW();
148 if (Oi==0.) continue;
149 double ai = 0.5*(bin.xMax()-bin.xMin());
150 double bi = 0.5*ai*(bin.xMax()+bin.xMin());
151 double Ei = bin.errW();
152 sum1 += sqr(bi/Ei);
153 sum2 += bi/sqr(Ei)*(Oi-ai);
154 }
155 return make_pair(sum2/sum1,sqrt(1./sum1));
156 }
157
158 vector<pair<double,double>> calcBetaGamma(Histo1DPtr hist) const {
159 if(hist->numEntries()==0.) return vector<pair<double,double>>(2,make_pair(0.,0.));;
160 double sum[7]={0.,0.,0.,0.,0.,0.,0.,};
161 for (const auto& bin : hist->bins()) {
162 double Oi = bin.sumW();
163 if(Oi==0.) continue;
164 double ai = (bin.xMax()-bin.xMin())/2./M_PI;
165 double bi = M_PI/32.*(cos(bin.xMin())-cos(bin.xMax()));
166 double ci =-M_PI/32.*(sin(bin.xMin())-sin(bin.xMax()));
167 double Ei = bin.errW();
168 sum[0] += bi*Oi/sqr(Ei);
169 sum[1] += ci*Oi/sqr(Ei);
170 sum[2] += ai*bi/sqr(Ei);
171 sum[3] += ai*ci/sqr(Ei);
172 sum[4] += sqr(bi)/sqr(Ei);
173 sum[5] += sqr(ci)/sqr(Ei);
174 sum[6] += bi*ci/sqr(Ei);
175 }
176 double gamma = (sum[0]/sum[4]-sum[2]/sum[4]-sum[1]/sum[6]+sum[3]/sum[6])/(sum[6]/sum[4]-sum[5]/sum[6]);
177 double beta = (sum[0]/sum[6]-sum[2]/sum[6]-sum[1]/sum[5]+sum[3]/sum[5])/(sum[4]/sum[6]-sum[6]/sum[5]);
178 return {make_pair(beta ,1./sqrt(sum[4])),make_pair(gamma,1./sqrt(sum[5]))};
179 }
180
181
182 /// Normalise histograms etc., after the run
183 void finalize() {
184 normalize(_h_cos1);
185 normalize(_h_cos2);
186 normalize(_h_phi2);
187 // compute values for specific modes
188 pair<double,double> a_b[2][2][2], a_c[2][2][2], a_s[2][2][2],beta[2][2][2],gamma[2][2][2];
189 for (unsigned int ip=0;ip<2;++ip) {
190 for (unsigned int ib=0;ib<2;++ib) {
191 for (unsigned int ic=0;ic<2;++ic) {
192 pair<double,double> p1 = calcAlpha(_h_cos1[ip][ib][ic]);
193 pair<double,double> p2 = calcAlpha(_h_cos2[ip][ib][ic]);
194 pair<double,double> p3 = make_pair(9.*_p_cos12[ip][ib][ic]->bin(1).mean(2),
195 9.*_p_cos12[ip][ib][ic]->bin(1).stdErr(2));
196 double ferr = sqrt(sqr(p1.second/p1.first)+sqr(p2.second/p2.first)+sqr(p3.second/p3.first));
197 // lamda decay first as sign well known
198 double l3 = p2.first*p3.first/p1.first;
199 double e3 = l3*ferr;
200 l3 = sqrt(max(l3,0.));
201 e3 = 0.5*e3;
202 // sign ambiguity so fix alpha for lambda0 -> p pi to known sign
203 if(ip==1) l3*=-1.;
204 a_s[ip][ib][ic] = make_pair(l3,e3);
205 // lambda_b decay
206 double l1 = p1.first*p3.first/p2.first;
207 double e1 = l1*ferr;
208 l1 = sqrt(max(l1,0.));
209 e1 = 0.5*e1;
210 if(p3.first*l3<0.) l1 *= -1.;
211 a_b[ip][ib][ic] = make_pair(l1,e1);
212 // lambda_c decay
213 double l2 = p1.first*p2.first/p3.first;
214 double e2 = l2*ferr;
215 l2 = sqrt(max(l2,0.));
216 e2 = 0.5*e2;
217 if(p2.first*l3<0.) l2 *= -1.;
218 a_c[ip][ib][ic] = make_pair(l2,e2);
219 vector<pair<double,double>> temp = calcBetaGamma(_h_phi2[ip][ib][ic]);
220 for(auto & val : temp) {
221 double error = sqrt(sqr(val.second/val.first)+sqr(p3.second/p3.first));
222 val.first /=p3.first;
223 val.second = val.first*error;
224 }
225 beta [ip][ib][ic] = temp[0];
226 gamma[ip][ib][ic] = temp[1];
227 }
228 }
229 }
230 // perform averages
231 // Lambda^0 -> p pi
232 pair<double,double> alam[2];
233 for(unsigned int ip=0;ip<2;++ip) {
234 double sum1(0.),sum2(0.);
235 for(unsigned int ib=0;ib<2;++ib) {
236 for(unsigned int ic=0;ic<2;++ic) {
237 sum1 += a_s[ip][ib][ic].first/sqr(a_s[ip][ib][ic].second);
238 sum2 += 1./sqr(a_s[ip][ib][ic].second);
239 }
240 }
241 alam[ip] = make_pair(sum1/sum2,sqrt(1./sum2));
242 Estimate0DPtr tmp;
243 book(tmp,1,6,1+ip);
244 tmp->set(alam[ip].first,alam[ip].second);
245 }
246 Estimate0DPtr tmp;
247 book(tmp,1,6,3);
248 tmp->set(0.5*(alam[0].first-alam[1].first),
249 0.5*sqrt(sqr(alam[0].second)+sqr(alam[1].second)));
250 book(tmp,1,6,4);
251 tmp->set((alam[0].first+alam[1].first)/(alam[0].first-alam[1].first),
252 4.*(sqr(alam[0].first*alam[1].second)+sqr(alam[1].first*alam[0].second))/pow(alam[0].first-alam[1].first,4));
253 // Lambda_b decays
254 pair<double,double> ab[2][2];
255 for(unsigned int ib=0;ib<2;++ib) {
256 for(unsigned int ip=0;ip<2;++ip) {
257 double sum1(0.),sum2(0.);
258 for(unsigned int ic=0;ic<2;++ic) {
259 sum1 += a_b[ip][ib][ic].first/sqr(a_b[ip][ib][ic].second);
260 sum2 += 1./sqr(a_b[ip][ib][ic].second);
261 }
262 ab[ip][ib] = make_pair(sum1/sum2,sqrt(1./sum2));
263 Estimate0DPtr tmp;
264 book(tmp,1,1+ib,1+ip);
265 tmp->set(ab[ip][ib].first,ab[ip][ib].second);
266 }
267 Estimate0DPtr tmp;
268 book(tmp,1,1+ib,3);
269 tmp->set(0.5*(ab[0][ib].first-ab[1][ib].first),
270 0.5*sqrt(sqr(ab[0][ib].second)+sqr(ab[1][ib].second)));
271 book(tmp,1,1+ib,4);
272 tmp->set((ab[0][ib].first+ab[1][ib].first)/(ab[0][ib].first-ab[1][ib].first),
273 4.*(sqr(ab[0][ib].first*ab[1][ib].second)+sqr(ab[1][ib].first*ab[0][ib].second))/pow(ab[0][ib].first-ab[1][ib].first,4));
274 }
275 // Lambda_c decays
276 for(unsigned int ic=0;ic<2;++ic) {
277 pair<double,double> ac[2],bc[2],gc[2];
278 for(unsigned int ip=0;ip<2;++ip) {
279 double sum1(0.),sum2(0.);
280 for(unsigned int ib=0;ib<2;++ib) {
281 sum1 += a_c[ip][ib][ic].first/sqr(a_c[ip][ib][ic].second);
282 sum2 += 1./sqr(a_c[ip][ib][ic].second);
283 }
284 ac[ip] = make_pair(sum1/sum2,sqrt(1./sum2));
285 Estimate0DPtr tmp;
286 book(tmp,1,3+ic,1+ip);
287 tmp->set(ac[ip].first,ac[ip].second);
288 sum1=sum2=0.;
289 for(unsigned int ib=0;ib<2;++ib) {
290 sum1 += beta[ip][ib][ic].first/sqr(beta[ip][ib][ic].second);
291 sum2 += 1./sqr(beta[ip][ib][ic].second);
292 }
293 bc[ip] = make_pair(sum1/sum2,sqrt(1./sum2));
294 book(tmp,2,1+ip,1+ic);
295 tmp->set(bc[ip].first,bc[ip].second);
296 sum1=sum2=0.;
297 for(unsigned int ib=0;ib<2;++ib) {
298 sum1 += gamma[ip][ib][ic].first/sqr(gamma[ip][ib][ic].second);
299 sum2 += 1./sqr(gamma[ip][ib][ic].second);
300 }
301 gc[ip] = make_pair(sum1/sum2,sqrt(1./sum2));
302 book(tmp,2,3+ip,1+ic);
303 tmp->set(gc[ip].first,gc[ip].second);
304
305
306
307
308
309 }
310 Estimate0DPtr tmp;
311 book(tmp,1,3+ic,3);
312 tmp->set(0.5*(ac[0].first-ac[1].first),
313 0.5*sqrt(sqr(ac[0].second)+sqr(ac[1].second)));
314 book(tmp,1,3+ic,4);
315 tmp->set((ac[0].first+ac[1].first)/(ac[0].first-ac[1].first),
316 4.*(sqr(ac[0].first*ac[1].second)+sqr(ac[1].first*ac[0].second))/pow(ac[0].first-ac[1].first,4));
317 }
318 // finally lambda_c -> p KS0
319 pair<double,double> ap[2];
320 for(unsigned int ip=0;ip<2;++ip) {
321 double sum1(0.),sum2(0.);
322 for(unsigned int ib=0;ib<2;++ib) {
323 pair<double,double> p1 = calcAlpha(_h_cos1[ip][ib][2]);
324 double val = p1.first/ab[ip][ib].first;
325 double err = sqr(val)*(sqr(p1.second/p1.first) + sqr(ab[ip][ib].second/ab[ip][ib].first));
326 sum1 += val/err;
327 sum2 += 1./err;
328 }
329 ap[ip] = make_pair(sum1/sum2,sqrt(1./sum2));
330 Estimate0DPtr tmp;
331 book(tmp,1,5,1+ip);
332 tmp->set(ap[ip].first,ap[ip].second);
333 }
334 book(tmp,1,5,3);
335 tmp->set(0.5*(ap[0].first-ap[1].first),
336 0.5*sqrt(sqr(ap[0].second)+sqr(ap[1].second)));
337 book(tmp,1,5,4);
338 tmp->set((ap[0].first+ap[1].first)/(ap[0].first-ap[1].first),
339 4.*(sqr(ap[0].first*ap[1].second)+sqr(ap[1].first*ap[0].second))/pow(ap[0].first-ap[1].first,4));
340 }
341 /// @}
342
343
344 /// @name Histograms
345 /// @{
346 Histo1DPtr _h_cos1[2][2][3],_h_cos2[2][2][2],_h_phi2[2][2][2];
347 Profile1DPtr _p_cos12[2][2][2];
348 /// @}
349
350
351 };
352
353
354 RIVET_DECLARE_PLUGIN(LHCB_2024_I2824757);
355
356}
|
