Rivet analyses referenceNA48_2010_I868871Decay asymmetries in $\Xi^0\to\Lambda^0\gamma$, $\Lambda^0\pi^0$ and $\Sigma^0\gamma$Experiment: NA48 () Inspire ID: 868871 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Measurement of the decay asymmetries in $\Xi^0\to\Lambda^0\gamma$, $\Lambda^0\pi^0$ and $\Sigma^0\gamma$ by the NA48 experiment. The asymmetry parameter is extracted by fitting to normalised angular distribution. This analysis is useful for testing spin correlations in hadron decays. Source code: NA48_2010_I868871.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include <sstream>
5namespace Rivet {
6
7
8 /// @brief asymmetrics in Xi0 decays
9 class NA48_2010_I868871 : public Analysis {
10 public:
11
12 /// Constructor
13 RIVET_DEFAULT_ANALYSIS_CTOR(NA48_2010_I868871);
14
15
16 /// @name Analysis methods
17 //@{
18
19 /// Book histograms and initialise projections before the run
20 void init() {
21
22 // Initialise and register projections
23 declare(UnstableParticles(), "UFS" );
24 // Book histograms
25 book(_h_ctheta_pi0 , "ctheta_pi0" , 20,-1,1);
26 book(_h_ctheta_gamma, "ctheta_gamma", 20,-1,1);
27 double step=0.1;
28 double xmin=-1.;
29 for(unsigned int ix=0;ix<20;++ix) {
30 Histo1DPtr temp;
31 std::ostringstream title;
32 title << "ctheta_Sigma_" << ix;
33 book(temp,title.str(), 20,-1,1);
34 _h_ctheta_Sigma.add(xmin, xmin+step, temp);
35 xmin+=step;
36 }
37 _nSigma=0.;
38 }
39
40
41 /// Perform the per-event analysis
42 void analyze(const Event& event) {
43 // loop over Omega baryons
44 for(const Particle& Xi : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==3322)) {
45 int sign = Xi.pid()/3322;
46 if(Xi.children().size()!=2) continue;
47 Particle baryon1,meson1;
48 unsigned int mode(0);
49 if(Xi.children()[0].pid()==sign*3122 &&
50 Xi.children()[1].pid()==111) {
51 baryon1 = Xi.children()[0];
52 meson1 = Xi.children()[1];
53 mode=1;
54 }
55 else if(Xi.children()[1].pid()==sign*3122 &&
56 Xi.children()[0].pid()==111) {
57 baryon1 = Xi.children()[1];
58 meson1 = Xi.children()[0];
59 mode=1;
60 }
61 else if(Xi.children()[0].pid()==sign*3122 &&
62 Xi.children()[1].pid()==22) {
63 baryon1 = Xi.children()[0];
64 meson1 = Xi.children()[1];
65 mode=2;
66 }
67 else if(Xi.children()[1].pid()==sign*3122 &&
68 Xi.children()[0].pid()==22) {
69 baryon1 = Xi.children()[1];
70 meson1 = Xi.children()[0];
71 mode=2;
72 }
73 else if(Xi.children()[0].pid()==sign*3212 &&
74 Xi.children()[1].pid()==22) {
75 baryon1 = Xi.children()[0];
76 meson1 = Xi.children()[1];
77 mode=3;
78 }
79 else if(Xi.children()[1].pid()==sign*3212 &&
80 Xi.children()[0].pid()==22) {
81 baryon1 = Xi.children()[1];
82 meson1 = Xi.children()[0];
83 mode=3;
84 }
85 else
86 continue;
87 if(baryon1.children().size()!=2) continue;
88 Particle baryon2,meson2,baryon3,meson3;
89 if(mode==1 || mode ==2) {
90 if(baryon1.children()[0].pid()== sign*2212 &&
91 baryon1.children()[1].pid()==-sign*211) {
92 baryon2 = baryon1.children()[0];
93 meson2 = baryon1.children()[1];
94 }
95 else if(baryon1.children()[1].pid()== sign*2212 &&
96 baryon1.children()[0].pid()==-sign*211) {
97 baryon2 = baryon1.children()[1];
98 meson2 = baryon1.children()[0];
99 }
100 else
101 continue;
102 }
103 else if(mode==3) {
104 if(baryon1.children()[0].pid()== sign*3122 &&
105 baryon1.children()[1].pid()== 22) {
106 baryon2 = baryon1.children()[0];
107 meson2 = baryon1.children()[1];
108 }
109 else if(baryon1.children()[1].pid()== sign*3122 &&
110 baryon1.children()[0].pid()== 22) {
111 baryon2 = baryon1.children()[1];
112 meson2 = baryon1.children()[0];
113 }
114 else
115 continue;
116 if(baryon2.children()[0].pid()== sign*2212 &&
117 baryon2.children()[1].pid()==-sign*211) {
118 baryon3 = baryon2.children()[0];
119 meson3 = baryon2.children()[1];
120 }
121 else if(baryon2.children()[1].pid()== sign*2212 &&
122 baryon2.children()[0].pid()==-sign*211) {
123 baryon3 = baryon2.children()[1];
124 meson3 = baryon2.children()[0];
125 }
126 else
127 continue;
128 }
129 // first boost to the Xi rest frame
130 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Xi.momentum().betaVec());
131 FourMomentum pbaryon1 = boost1.transform(baryon1.momentum());
132 FourMomentum pbaryon2 = boost1.transform(baryon2.momentum());
133 // to lambda rest frame
134 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pbaryon1.betaVec());
135 Vector3 axis = pbaryon1.p3().unit();
136 FourMomentum pp = boost2.transform(pbaryon2);
137 // calculate angle
138 double cTheta = pp.p3().unit().dot(axis);
139 if(mode==1) {
140 _h_ctheta_pi0->fill(cTheta,1.);
141 }
142 else if(mode==2) {
143 _h_ctheta_gamma->fill(cTheta,1.);
144 }
145 else if(mode==3) {
146 FourMomentum pbaryon3 = boost1.transform(baryon3.momentum());
147 FourMomentum pp2 = boost2.transform(pbaryon3);
148 Vector3 axis2 = pp.p3().unit();
149 double cTheta2 = pp2.p3().unit().dot(axis2);
150 _h_ctheta_Sigma.fill(cTheta,cTheta2,1.);
151 _nSigma += 1.;
152 }
153 }
154 }
155
156 pair<double,double> calcAlpha(Histo1DPtr hist) {
157 if(hist->numEntries()==0.) return make_pair(0.,0.);
158 double sum1(0.),sum2(0.);
159 for (auto bin : hist->bins() ) {
160 double Oi = bin.area();
161 if(Oi==0.) continue;
162 double ai = 0.5*(bin.xMax()-bin.xMin());
163 double bi = 0.5*ai*(bin.xMax()+bin.xMin());
164 double Ei = bin.areaErr();
165 sum1 += sqr(bi/Ei);
166 sum2 += bi/sqr(Ei)*(Oi-ai);
167 }
168 return make_pair(sum2/sum1,sqrt(1./sum1));
169 }
170
171 pair<double,double> calcAlpha(BinnedHistogram hist) {
172 double sum1(0.),sum2(0.);
173 double step=0.1;
174 double xmin=-1.;
175 for(unsigned int ix=0;ix<20;++ix) {
176 double xsum=2.*xmin+step;
177 Histo1DPtr h2 = hist.histo(xmin+0.5*step);
178 for (auto bin : h2->bins() ) {
179 double Oi = bin.area();
180 if(Oi==0.) continue;
181 double ai = 0.25*(bin.xMax()-bin.xMin())*step;
182 double bi = 0.25*ai*(bin.xMax()+bin.xMin())*xsum;
183 double Ei = bin.areaErr();
184 sum1 += sqr(bi/Ei);
185 sum2 += bi/sqr(Ei)*(Oi-ai);
186 }
187 xmin+=step;
188 }
189 return make_pair(sum2/sum1,sqrt(1./sum1));
190 }
191
192 /// Normalise histograms etc., after the run
193 void finalize() {
194 // Xi0 -> Lambda0 pi0
195 normalize(_h_ctheta_pi0);
196 Scatter2DPtr _h_alpha_pi0;
197 book(_h_alpha_pi0,1,1,1);
198 pair<double,double> alpha = calcAlpha(_h_ctheta_pi0);
199 _h_alpha_pi0->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(alpha.second,alpha.second) );
200 // Xi0 -> Lambda gamma (N.B. sign due defns)
201 normalize(_h_ctheta_gamma);
202 Scatter2DPtr _h_alpha_gamma;
203 book(_h_alpha_gamma,2,1,1);
204 alpha = calcAlpha(_h_ctheta_gamma);
205 _h_alpha_gamma->addPoint(0.5,-alpha.first, make_pair(0.5,0.5), make_pair(alpha.second,alpha.second) );
206 // Xi0 -> Sigma gamma
207 _h_ctheta_Sigma.scale(1./_nSigma,this);
208 Scatter2DPtr _h_alpha_Sigma;
209 book(_h_alpha_Sigma,3,1,1);
210 alpha = calcAlpha(_h_ctheta_Sigma);
211 _h_alpha_Sigma->addPoint(0.5,alpha.first, make_pair(0.5,0.5), make_pair(alpha.second,alpha.second) );
212 }
213
214 //@}
215
216
217 /// @name Histograms
218 //@{
219 Histo1DPtr _h_ctheta_pi0,_h_ctheta_gamma;
220 BinnedHistogram _h_ctheta_Sigma;
221 double _nSigma;
222 //@}
223
224 };
225
226
227 // The hook for the plugin system
228 RIVET_DECLARE_PLUGIN(NA48_2010_I868871);
229
230
231}
|