Rivet analyses referenceOPAL_1997_I447188Polarization of $\Lambda^0$ baryons at LEP 1Experiment: OPAL (LEP) Inspire ID: 447188 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Measurement of the polarization of $\Lambda^0$ baryons at LEP 1. The $\cos\theta$ and $\cos\phi$ distributions are extracted and then the polarization obtained by fitting the resulting distributions. Source code: OPAL_1997_I447188.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5#include "Rivet/Projections/Thrust.hh"
6#include "Rivet/Tools/BinnedHistogram.hh"
7#include "Rivet/Projections/ChargedFinalState.hh"
8#include <sstream>
9namespace Rivet {
10
11
12 /// @brief Lambda polarization at LEP1
13 class OPAL_1997_I447188 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1997_I447188);
18
19
20 /// @name Analysis methods
21 //@{
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25
26 // Initialise and register projections
27 declare(Beam(), "Beams");
28 const ChargedFinalState cfs;
29 const Thrust thrust(cfs);
30 declare(thrust, "Thrust");
31 declare(UnstableParticles(), "UFS");
32
33 // Book the histograms
34 {Histo1DPtr temp; _h_ctheta.add(0.027, 0.05 ,book(temp, "/TMP/ctheta_0",20,-1.,1.));}
35 {Histo1DPtr temp; _h_ctheta.add(0.05 , 0.08 ,book(temp, 4,1,1));}
36 {Histo1DPtr temp; _h_ctheta.add(0.08 , 0.09 ,book(temp, "/TMP/ctheta_2",20,-1.,1.));}
37 {Histo1DPtr temp; _h_ctheta.add(0.09 , 0.1 ,book(temp, "/TMP/ctheta_3",20,-1.,1.));}
38 {Histo1DPtr temp; _h_ctheta.add(0.1 , 0.15 ,book(temp, 4,1,2));}
39 {Histo1DPtr temp; _h_ctheta.add(0.15 , 0.2 ,book(temp, "/TMP/ctheta_5",20,-1.,1.));}
40 {Histo1DPtr temp; _h_ctheta.add(0.2 , 0.3 ,book(temp, 4,1,3));}
41 {Histo1DPtr temp; _h_ctheta.add(0.3 , 0.4 ,book(temp, "/TMP/ctheta_7",20,-1.,1.));}
42 {Histo1DPtr temp; _h_ctheta.add(0.4 , 1.0 ,book(temp, "/TMP/ctheta_8",20,-1.,1.));}
43 book(_h_ctheta_large, 4,1,4);
44
45 {Histo1DPtr temp; _h_cphi .add(0.3,0.6,book(temp, "/TMP/cphiP_0",10,0.,1.));}
46 {Histo1DPtr temp; _h_cphi .add(0.3,0.6,book(temp, 5,1,1));}
47 {Histo1DPtr temp; _h_cphi .add(0.6,0.9,book(temp, 5,1,2));}
48 {Histo1DPtr temp; _h_cphi .add(0.9,1.2,book(temp, "/TMP/cphiP_3",10,0.,1.));}
49 {Histo1DPtr temp; _h_cphi .add(1.2,1.5,book(temp, "/TMP/cphiP_4",10,0.,1.));}
50
51 book(_h_cphi_low , 5,1,4);
52 book(_h_cphi_mid , "/TMP/cphiP_mid" ,10,0.,1.);
53 book(_h_cphi_high, 5,1,3);
54
55 {Histo1DPtr temp; _h_plus_lam.add(0.027,0.05,book(temp, "/TMP/lamP_0",20,-1.,1.));}
56 {Histo1DPtr temp; _h_plus_lam.add(0.05 ,0.08,book(temp, "/TMP/lamP_1",20,-1.,1.));}
57 {Histo1DPtr temp; _h_plus_lam.add(0.08 ,0.09,book(temp, "/TMP/lamP_2",20,-1.,1.));}
58 {Histo1DPtr temp; _h_plus_lam.add(0.09 ,0.1 ,book(temp, "/TMP/lamP_3",20,-1.,1.));}
59 {Histo1DPtr temp; _h_plus_lam.add(0.1 ,0.15,book(temp, "/TMP/lamP_4",20,-1.,1.));}
60 {Histo1DPtr temp; _h_plus_lam.add(0.15 ,0.2 ,book(temp, "/TMP/lamP_5",20,-1.,1.));}
61 {Histo1DPtr temp; _h_plus_lam.add(0.2 ,0.3 ,book(temp, "/TMP/lamP_6",20,-1.,1.));}
62 {Histo1DPtr temp; _h_plus_lam.add(0.3 ,0.4 ,book(temp, "/TMP/lamP_7",20,-1.,1.));}
63 {Histo1DPtr temp; _h_plus_lam.add(0.4 ,1.0 ,book(temp, "/TMP/lamP_8",20,-1.,1.));}
64
65 {Histo1DPtr temp; _h_minus_lam.add(0.027,0.05,book(temp, "/TMP/lamM_0",20,-1.,1.));}
66 {Histo1DPtr temp; _h_minus_lam.add(0.05 ,0.08,book(temp, "/TMP/lamM_1",20,-1.,1.));}
67 {Histo1DPtr temp; _h_minus_lam.add(0.08 ,0.09,book(temp, "/TMP/lamM_2",20,-1.,1.));}
68 {Histo1DPtr temp; _h_minus_lam.add(0.09 ,0.1 ,book(temp, "/TMP/lamM_3",20,-1.,1.));}
69 {Histo1DPtr temp; _h_minus_lam.add(0.1 ,0.15,book(temp, "/TMP/lamM_4",20,-1.,1.));}
70 {Histo1DPtr temp; _h_minus_lam.add(0.15 ,0.2 ,book(temp, "/TMP/lamM_5",20,-1.,1.));}
71 {Histo1DPtr temp; _h_minus_lam.add(0.2 ,0.3 ,book(temp, "/TMP/lamM_6",20,-1.,1.));}
72 {Histo1DPtr temp; _h_minus_lam.add(0.3 ,0.4 ,book(temp, "/TMP/lamM_7",20,-1.,1.));}
73 {Histo1DPtr temp; _h_minus_lam.add(0.4 ,1.0 ,book(temp, "/TMP/lamM_8",20,-1.,1.));}
74
75 {Histo1DPtr temp; _h_plus_lam_large1 = book(temp, "/TMP/lamP_large_1",20,-1.,1.);}
76 {Histo1DPtr temp; _h_plus_lam_large2 = book(temp, "/TMP/lamP_large_2",20,-1.,1.);}
77 {Histo1DPtr temp; _h_minus_lam_large1 = book(temp, "/TMP/lamM_large_1",20,-1.,1.);}
78 {Histo1DPtr temp; _h_minus_lam_large2 = book(temp, "/TMP/lamM_large_2",20,-1.,1.);}
79 }
80
81
82 /// Perform the per-event analysis
83 void analyze(const Event& event) {
84
85 // Get beams and average beam momentum
86 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
87 const double meanBeamMom = ( beams.first.p3().mod() +
88 beams.second.p3().mod() ) / 2.0;
89 Vector3 beamAxis;
90 if(beams.first.pid()==-11) {
91 beamAxis = beams.first .momentum().p3().unit();
92 }
93 else {
94 beamAxis = beams.second.momentum().p3().unit();
95 }
96
97 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
98 // thrust, to define an axis
99 const Thrust& thrust = apply<Thrust>(event, "Thrust");
100 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
101
102 for(const Particle & lambda : ufs.particles(Cuts::abspid==3122)) {
103 double xE = lambda.momentum().t()/meanBeamMom;
104 int sign = lambda.pid()/3122;
105 Vector3 axis1 = lambda.momentum().p3().unit();
106 // assymetry
107 double cLam = axis1.dot(beamAxis);
108 if(sign>0)
109 _h_plus_lam .fill(xE,cLam);
110 else
111 _h_minus_lam.fill(xE,cLam);
112 if(xE>0.15) {
113 if(sign>0)
114 _h_plus_lam_large1 ->fill(cLam);
115 else
116 _h_minus_lam_large1->fill(cLam);
117 }
118 if(xE>0.3) {
119 if(sign>0)
120 _h_plus_lam_large2 ->fill(cLam);
121 else
122 _h_minus_lam_large2->fill(cLam);
123 }
124 if(lambda.children().size()!=2) continue;
125 // look at the decay products
126 Particle proton,pion;
127 if(lambda.children()[0].pid()==sign*2212 &&
128 lambda.children()[1].pid()==-sign*211) {
129 proton = lambda.children()[0];
130 pion = lambda.children()[1];
131 }
132 else if(lambda.children()[1].pid()==sign*2212 &&
133 lambda.children()[0].pid()==-sign*211) {
134 proton = lambda.children()[1];
135 pion = lambda.children()[0];
136 }
137 else
138 continue;
139 // boost to the lambda rest frame
140 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(lambda.momentum().betaVec());
141 FourMomentum pproton = boost.transform(proton.momentum());
142 // longitudinal polarization
143 double ctheta = axis1.dot(pproton.p3().unit());
144 _h_ctheta.fill(xE,ctheta);
145 if(xE>=0.3) _h_ctheta_large->fill(ctheta);
146 // transverse polarization
147 Vector3 axis2;
148 if(lambda.momentum().p3().dot(thrust.thrustAxis())>=0.) {
149 axis2 = thrust.thrustAxis();
150 }
151 else {
152 axis2 =-thrust.thrustAxis();
153 }
154 Vector3 axis3 = axis2.cross(axis1).unit();
155 double pT = sqrt(sqr(thrust.thrustMajorAxis().dot(lambda.momentum().p3()))+
156 sqr(thrust.thrustMinorAxis().dot(lambda.momentum().p3())));
157 double cPhi = axis3.dot(pproton.p3().unit());
158 _h_cphi.fill(pT,cPhi);
159 if(pT>0.3) _h_cphi_low ->fill(cPhi);
160 if(pT>0.6) _h_cphi_mid ->fill(cPhi);
161 if(pT>1.5) _h_cphi_high->fill(cPhi);
162 }
163
164 }
165
166 pair<double,double> calcAlpha(Histo1DPtr hist) {
167 if(hist->numEntries()==0.) return make_pair(0.,0.);
168 double sum1(0.),sum2(0.);
169 for (auto bin : hist->bins() ) {
170 double Oi = bin.area();
171 if(Oi==0.) continue;
172 double ai = 0.5*(bin.xMax()-bin.xMin());
173 double bi = 0.5*ai*(bin.xMax()+bin.xMin());
174 double Ei = bin.areaErr();
175 sum1 += sqr(bi/Ei);
176 sum2 += bi/sqr(Ei)*(Oi-ai);
177 }
178 return make_pair(sum2/sum1,sqrt(1./sum1));
179 }
180
181 pair<double,double> calcAsymmetry(Scatter2DPtr hist,unsigned int mode) {
182 double sum1(0.),sum2(0.);
183 for (auto bin : hist->points() ) {
184 double Oi = bin.y();
185 if(Oi==0.) continue;
186 double bi;
187 if(mode==0)
188 bi = 0.25*(bin.xMax()-bin.xMin())*(bin.xMax()+bin.xMin());
189 else
190 bi = 4.*(bin.xMax()+bin.xMin())/(3.+sqr(bin.xMax())+bin.xMax()*bin.xMin()+sqr(bin.xMin()));
191 double Ei = bin.yErrAvg();
192 sum1 += sqr(bi/Ei);
193 sum2 += bi/sqr(Ei)*Oi;
194 }
195 return make_pair(sum2/sum1,sqrt(1./sum1));
196 }
197
198
199 /// Normalise histograms etc., after the run
200 void finalize() {
201 // longitudinal polarization
202 vector<double> x_val = {3.850000e-02,6.500000e-02,8.500000e-02,9.500000e-02,1.250000e-01,
203 1.750000e-01,2.500000e-01,3.500000e-01,7.000000e-01};
204 vector<double> x_err = {1.150000e-02,1.500000e-02,5.000000e-03,5.000000e-03,2.500000e-02,
205 2.500000e-02,5.000000e-02,5.000000e-02,3.000000e-01};
206 unsigned int ipoint=0;
207 double aLam = 0.642;
208 Scatter2DPtr h_long;
209 book(h_long, 1,1,1);
210 for(Histo1DPtr & hist : _h_ctheta.histos()) {
211 normalize(hist);
212 pair<double,double> alpha = calcAlpha(hist);
213 alpha.first /=aLam;
214 alpha.second /=aLam;
215 h_long->addPoint(x_val[ipoint], 100.*alpha.first, make_pair(x_err[ipoint],x_err[ipoint]),
216 make_pair(100.*alpha.second,100.*alpha.second) );
217 ++ipoint;
218 }
219 normalize(_h_ctheta_large);
220 pair<double,double> alpha = calcAlpha(_h_ctheta_large);
221 alpha.first /=aLam;
222 alpha.second /=aLam;
223 Scatter2DPtr h_long_l;
224 book(h_long_l,1,1,2);
225 h_long_l->addPoint(0.65, 100.*alpha.first, make_pair(0.35,0.35), make_pair(100.*alpha.second,100.*alpha.second) );
226 // transverse polarization
227 vector<double> pT_val = {0.15,0.45,0.75,1.05,1.35};
228 vector<double> pT_err = {0.15,0.15,0.15,0.15,0.15};
229 Scatter2DPtr h_trans;
230 book(h_trans,2,1,1);
231 ipoint=0;
232 for(Histo1DPtr & hist : _h_cphi.histos()) {
233 normalize(hist);
234 pair<double,double> alpha = calcAlpha(hist);
235 alpha.first /=aLam;
236 alpha.second /=aLam;
237 h_trans->addPoint(pT_val[ipoint], 100.*alpha.first, make_pair(pT_err[ipoint],pT_err[ipoint]),
238 make_pair(100.*alpha.second,100.*alpha.second) );
239 ++ipoint;
240 }
241 normalize(_h_cphi_low);
242 alpha = calcAlpha(_h_cphi_low);
243 alpha.first /=aLam;
244 alpha.second /=aLam;
245 Scatter2DPtr h_trans_low;
246 book(h_trans_low,2,1,2);
247 h_trans_low->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(100.*alpha.second,100.*alpha.second) );
248 normalize(_h_cphi_mid);
249 alpha = calcAlpha(_h_cphi_mid);
250 alpha.first /=aLam;
251 alpha.second /=aLam;
252 Scatter2DPtr h_trans_mid;
253 book(h_trans_mid,2,1,3);
254 h_trans_mid->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(100.*alpha.second,100.*alpha.second) );
255 normalize(_h_cphi_high);
256 alpha = calcAlpha(_h_cphi_high);
257 alpha.first /=aLam;
258 alpha.second /=aLam;
259 Scatter2DPtr h_trans_high;
260 book(h_trans_high,2,1,4);
261 h_trans_high->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(100.*alpha.second,100.*alpha.second) );
262 // asyymetry
263 Scatter2DPtr h_asym;
264 book(h_asym,3,1,1);
265 for(unsigned int ix=0;ix<9;++ix) {
266 normalize(_h_plus_lam.histos()[ix] );
267 normalize(_h_minus_lam.histos()[ix]);
268 std::ostringstream title;
269 title << "/TMP/a_lam_" << ix;
270 Scatter2DPtr sTemp;
271 book(sTemp,title.str());
272 asymm(_h_plus_lam.histos()[ix],_h_minus_lam.histos()[ix],sTemp);
273 pair<double,double> alpha = calcAsymmetry(sTemp,1);
274 h_asym->addPoint(x_val[ix],-alpha.first, make_pair(x_err[ix],x_err[ix]),
275 make_pair(alpha.second,alpha.second) );
276 }
277 normalize(_h_plus_lam_large1 );
278 normalize(_h_minus_lam_large1);
279 Scatter2DPtr sTemp;
280 book(sTemp,"/TMP/a_lam_large1");
281 asymm(_h_plus_lam_large1,_h_minus_lam_large1,sTemp);
282 alpha = calcAsymmetry(sTemp,1);
283 book(h_asym,3,1,2);
284 h_asym->addPoint(0.575,-alpha.first, make_pair(0.425,0.425),
285 make_pair(alpha.second,alpha.second) );
286 normalize(_h_plus_lam_large2 );
287 normalize(_h_minus_lam_large2);
288 book(sTemp,"/TMP/a_lam_large2");
289 asymm(_h_plus_lam_large2,_h_minus_lam_large2,sTemp);
290 alpha = calcAsymmetry(sTemp,1);
291 book(h_asym,3,1,3);
292 h_asym->addPoint(0.65,-alpha.first, make_pair(0.35,0.35),
293 make_pair(alpha.second,alpha.second) );
294 }
295
296 //@}
297
298
299 /// @name Histograms
300 //@{
301 BinnedHistogram _h_ctheta,_h_cphi,_h_plus_lam,_h_minus_lam;
302 Histo1DPtr _h_ctheta_large;
303 Histo1DPtr _h_cphi_low, _h_cphi_mid, _h_cphi_high;
304 Histo1DPtr _h_plus_lam_large1,_h_plus_lam_large2;
305 Histo1DPtr _h_minus_lam_large1,_h_minus_lam_large2;
306 //@}
307
308 };
309
310
311 // The hook for the plugin system
312 RIVET_DECLARE_PLUGIN(OPAL_1997_I447188);
313
314
315}
|