Rivet analyses referenceALEPH_1996_I415745Polarization of $\Lambda^0$ baryons at LEP 1Experiment: ALEPH (LEP) Inspire ID: 415745 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: ALEPH_1996_I415745.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>
9
10namespace Rivet {
11
12
13 /// @brief Lambda polarization at LEP1
14 class ALEPH_1996_I415745 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_1996_I415745);
19
20
21 /// @name Analysis methods
22 //@{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 // Initialise and register projections
28 declare(Beam(), "Beams");
29 const ChargedFinalState cfs;
30 const Thrust thrust(cfs);
31 declare(thrust, "Thrust");
32 declare(UnstableParticles(), "UFS");
33
34 // Book histograms
35 {Histo1DPtr temp; _h_ctheta.add(0.1 ,0.15,book(temp, "/TMP/ctheta_0",20,-1.,1.));}
36 {Histo1DPtr temp; _h_ctheta.add(0.15,0.2 ,book(temp, "/TMP/ctheta_1",20,-1.,1.));}
37 {Histo1DPtr temp; _h_ctheta.add(0.2 ,0.3 ,book(temp, "/TMP/ctheta_2",20,-1.,1.));}
38 {Histo1DPtr temp; _h_ctheta.add(0.3 ,0.4 ,book(temp, "/TMP/ctheta_3",20,-1.,1.));}
39 {Histo1DPtr temp; _h_ctheta.add(0.4 ,1. ,book(temp, "/TMP/ctheta_4",20,-1.,1.));}
40 book(_h_ctheta_large,"/TMP/ctheta_large",20,-1.,1.);
41
42 {Histo1DPtr temp; _h_plus_cphi .add(0.3,0.6,book(temp, "/TMP/cphiP_0",10,0.,1.));}
43 {Histo1DPtr temp; _h_plus_cphi .add(0.6,0.9,book(temp, "/TMP/cphiP_1",10,0.,1.));}
44 {Histo1DPtr temp; _h_plus_cphi .add(0.9,1.2,book(temp, "/TMP/cphiP_2",10,0.,1.));}
45 {Histo1DPtr temp; _h_plus_cphi .add(1.2,1.5,book(temp, "/TMP/cphiP_3",10,0.,1.));}
46 {Histo1DPtr temp; _h_minus_cphi.add(0.3,0.6,book(temp, "/TMP/cphiM_0",10,0.,1.));}
47 {Histo1DPtr temp; _h_minus_cphi.add(0.6,0.9,book(temp, "/TMP/cphiM_1",10,0.,1.));}
48 {Histo1DPtr temp; _h_minus_cphi.add(0.9,1.2,book(temp, "/TMP/cphiM_2",10,0.,1.));}
49 {Histo1DPtr temp; _h_minus_cphi.add(1.2,1.5,book(temp, "/TMP/cphiM_3",10,0.,1.));}
50 book(_h_plus_cphi_low , "/TMP/cphiP_low" ,10,0.,1.);
51 book(_h_plus_cphi_mid , "/TMP/cphiP_mid" ,10,0.,1.);
52 book(_h_plus_cphi_high , "/TMP/cphiP_high",10,0.,1.);
53 book(_h_minus_cphi_low , "/TMP/cphiM_low" ,10,0.,1.);
54 book(_h_minus_cphi_mid , "/TMP/cphiM_mid" ,10,0.,1.);
55 book(_h_minus_cphi_high, "/TMP/cphiM_high",10,0.,1.);
56
57 {Histo1DPtr temp; _h_plus_lam.add(0.1 ,0.15,book(temp, "/TMP/lamP_0",20,-1.,1.));}
58 {Histo1DPtr temp; _h_plus_lam.add(0.15,0.2 ,book(temp, "/TMP/lamP_1",20,-1.,1.));}
59 {Histo1DPtr temp; _h_plus_lam.add(0.2 ,0.3 ,book(temp, "/TMP/lamP_2",20,-1.,1.));}
60 {Histo1DPtr temp; _h_plus_lam.add(0.3 ,0.4 ,book(temp, "/TMP/lamP_3",20,-1.,1.));}
61 {Histo1DPtr temp; _h_plus_lam.add(0.4 ,1. ,book(temp, "/TMP/lamP_4",20,-1.,1.));}
62
63 {Histo1DPtr temp; _h_minus_lam.add(0.1 ,0.15,book(temp, "/TMP/lamM_0",20,-1.,1.));}
64 {Histo1DPtr temp; _h_minus_lam.add(0.15,0.2 ,book(temp, "/TMP/lamM_1",20,-1.,1.));}
65 {Histo1DPtr temp; _h_minus_lam.add(0.2 ,0.3 ,book(temp, "/TMP/lamM_2",20,-1.,1.));}
66 {Histo1DPtr temp; _h_minus_lam.add(0.3 ,0.4 ,book(temp, "/TMP/lamM_3",20,-1.,1.));}
67 {Histo1DPtr temp; _h_minus_lam.add(0.4 ,1. ,book(temp, "/TMP/lamM_4",20,-1.,1.));}
68
69 }
70
71
72 /// Perform the per-event analysis
73 void analyze(const Event& event) {
74
75 // Get beams and average beam momentum
76 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
77 const double meanBeamMom = ( beams.first.p3().mod() +
78 beams.second.p3().mod() ) / 2.0;
79 Vector3 beamAxis;
80 if(beams.first.pid()==-11) {
81 beamAxis = beams.first .momentum().p3().unit();
82 }
83 else {
84 beamAxis = beams.second.momentum().p3().unit();
85 }
86
87 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
88 // thrust, to define an axis
89 const Thrust& thrust = apply<Thrust>(event, "Thrust");
90 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
91
92 for(const Particle & lambda : ufs.particles(Cuts::abspid==3122)) {
93 double z = lambda.momentum().p3().mod()/meanBeamMom;
94 int sign = lambda.pid()/3122;
95 Vector3 axis1 = lambda.momentum().p3().unit();
96 // assymetry
97 double cLam = axis1.dot(beamAxis);
98 if(sign>0)
99 _h_plus_lam .fill(z,cLam);
100 else
101 _h_minus_lam.fill(z,cLam);
102 if(lambda.children().size()!=2) continue;
103 // look at the decay products
104 Particle proton,pion;
105 if(lambda.children()[0].pid()==sign*2212 &&
106 lambda.children()[1].pid()==-sign*211) {
107 proton = lambda.children()[0];
108 pion = lambda.children()[1];
109 }
110 else if(lambda.children()[1].pid()==sign*2212 &&
111 lambda.children()[0].pid()==-sign*211) {
112 proton = lambda.children()[1];
113 pion = lambda.children()[0];
114 }
115 else
116 continue;
117 // boost to the lambda rest frame
118 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(lambda.momentum().betaVec());
119 FourMomentum pproton = boost.transform(proton.momentum());
120 // longitudinal polarization
121 double ctheta = axis1.dot(pproton.p3().unit());
122 _h_ctheta.fill(z,ctheta);
123 if(z>=0.3) _h_ctheta_large->fill(ctheta);
124
125 // transverse polarization
126 if(z>0.15) {
127 Vector3 axis2;
128 if(lambda.momentum().p3().dot(thrust.thrustAxis())>=0.) {
129 axis2 = thrust.thrustAxis();
130 }
131 else {
132 axis2 =-thrust.thrustAxis();
133 }
134 Vector3 axis3 = axis2.cross(axis1).unit();
135
136 double pT = sqrt(sqr(thrust.thrustMajorAxis().dot(lambda.momentum().p3()))+
137 sqr(thrust.thrustMinorAxis().dot(lambda.momentum().p3())));
138 double cPhi = axis3.dot(pproton.p3().unit());
139 if(cPhi>0.) {
140 _h_plus_cphi .fill(pT,cPhi);
141 if(pT>0.3)
142 _h_plus_cphi_low ->fill(cPhi);
143 if(pT>0.6)
144 _h_plus_cphi_mid ->fill(cPhi);
145 if(pT>1.5)
146 _h_plus_cphi_high->fill(cPhi);
147 }
148 else {
149 _h_minus_cphi.fill(pT,abs(cPhi));
150 if(pT>0.3)
151 _h_minus_cphi_low ->fill(abs(cPhi));
152 if(pT>0.6)
153 _h_minus_cphi_mid ->fill(abs(cPhi));
154 if(pT>1.5)
155 _h_minus_cphi_high->fill(abs(cPhi));
156 }
157 }
158 }
159 }
160
161 pair<double,double> calcAlpha(Histo1DPtr hist) {
162 if(hist->numEntries()==0.) return make_pair(0.,0.);
163 double sum1(0.),sum2(0.);
164 for (auto bin : hist->bins() ) {
165 double Oi = bin.area();
166 if(Oi==0.) continue;
167 double ai = 0.5*(bin.xMax()-bin.xMin());
168 double bi = 0.5*ai*(bin.xMax()+bin.xMin());
169 double Ei = bin.areaErr();
170 sum1 += sqr(bi/Ei);
171 sum2 += bi/sqr(Ei)*(Oi-ai);
172 }
173 return make_pair(sum2/sum1,sqrt(1./sum1));
174 }
175
176 pair<double,double> calcAsymmetry(Scatter2DPtr hist,unsigned int mode) {
177 double sum1(0.),sum2(0.);
178 for (auto bin : hist->points() ) {
179 double Oi = bin.y();
180 if(Oi==0.) continue;
181 double bi;
182 if(mode==0)
183 bi = 0.25*(bin.xMax()-bin.xMin())*(bin.xMax()+bin.xMin());
184 else
185 bi = 4.*(bin.xMax()+bin.xMin())/(3.+sqr(bin.xMax())+bin.xMax()*bin.xMin()+sqr(bin.xMin()));
186 double Ei = bin.yErrAvg();
187 sum1 += sqr(bi/Ei);
188 sum2 += bi/sqr(Ei)*Oi;
189 }
190 return make_pair(sum2/sum1,sqrt(1./sum1));
191 }
192
193 /// Normalise histograms etc., after the run
194 void finalize() {
195 // longitudinal polarization
196 vector<double> x_val = {0.125 , 0.175,0.25,0.35,0.70};
197 vector<double> x_err = {0.025 , 0.025,0.05,0.05,0.30};
198 unsigned int ipoint=0;
199 double aLam = 0.642;
200 Scatter2DPtr h_long;
201 book(h_long,1,1,1);
202 for(Histo1DPtr & hist : _h_ctheta.histos()) {
203 normalize(hist);
204 pair<double,double> alpha = calcAlpha(hist);
205 alpha.first /=aLam;
206 alpha.second /=aLam;
207 h_long->addPoint(x_val[ipoint], alpha.first, make_pair(x_err[ipoint],x_err[ipoint]),
208 make_pair(alpha.second,alpha.second) );
209 ++ipoint;
210 }
211 normalize(_h_ctheta_large);
212 pair<double,double> alpha = calcAlpha(_h_ctheta_large);
213 alpha.first /=aLam;
214 alpha.second /=aLam;
215 Scatter2DPtr h_long_l;
216 book(h_long_l,1,2,1);
217 h_long_l->addPoint(0.65, alpha.first, make_pair(0.35,0.35), make_pair(alpha.second,alpha.second) );
218 // transverse polarization
219 vector<double> pT_val = {0.45,0.75,1.05,1.35};
220 vector<double> pT_err = {0.15,0.15,0.15,0.15};
221 Scatter2DPtr h_trans;
222 book(h_trans,2,1,1);
223 for(unsigned int ix=0;ix<4;++ix) {
224 normalize(_h_plus_cphi.histos()[ix] );
225 normalize(_h_minus_cphi.histos()[ix]);
226 std::ostringstream title;
227 title << "/TMP/a_cphi_" << ix;
228 Scatter2DPtr sTemp;
229 book(sTemp,title.str());
230 asymm(_h_plus_cphi.histos()[ix],_h_minus_cphi.histos()[ix],sTemp);
231 pair<double,double> alpha = calcAsymmetry(sTemp,0);
232 alpha.first /=aLam;
233 alpha.second /=aLam;
234 h_trans->addPoint(pT_val[ix], alpha.first, make_pair(pT_err[ix],pT_err[ix]),
235 make_pair(alpha.second,alpha.second) );
236 }
237 Scatter2DPtr sLow;
238 book(sLow,"/TMP/a_cphi_low");
239 asymm(_h_plus_cphi_low,_h_minus_cphi_low,sLow);
240 alpha = calcAsymmetry(sLow,0);
241 alpha.first /=aLam;
242 alpha.second /=aLam;
243 Scatter2DPtr h_trans_low;
244 book(h_trans_low,2,3,1);
245 h_trans_low->addPoint(5.15, alpha.first, make_pair(4.85,4.85), make_pair(alpha.second,alpha.second) );
246
247 Scatter2DPtr sMid;
248 book(sMid,"/TMP/a_cphi_mid");
249 asymm(_h_plus_cphi_mid,_h_minus_cphi_mid,sMid);
250 alpha = calcAsymmetry(sMid,0);
251 alpha.first /=aLam;
252 alpha.second /=aLam;
253 Scatter2DPtr h_trans_mid;
254 book(h_trans_mid,2,4,1);
255 h_trans_mid->addPoint(5.3, alpha.first, make_pair(4.7,4.7), make_pair(alpha.second,alpha.second) );
256
257 Scatter2DPtr sHigh;
258 book(sHigh,"/TMP/a_cphi_high");
259 asymm(_h_plus_cphi_high,_h_minus_cphi_high,sHigh);
260 alpha = calcAsymmetry(sHigh,0);
261 alpha.first /=aLam;
262 alpha.second /=aLam;
263 Scatter2DPtr h_trans_high;
264 book(h_trans_high,2,2,1);
265 h_trans_high->addPoint(5.75, alpha.first, make_pair(4.25,4.25), make_pair(alpha.second,alpha.second) );
266
267 // asyymetry
268 vector<double> x_val2 = {1.220000e-01,1.730000e-01,2.410000e-01,3.420000e-01,4.950000e-01};
269 vector<double> x_err2 = {2.200000e-02,2.300000e-02,4.100000e-02,4.200000e-02,9.500000e-02};
270 vector<double> x_err3 = {2.800000e-02,2.700000e-02,5.900000e-02,5.800000e-02,5.050000e-01};
271 Scatter2DPtr h_asym;
272 book(h_asym,3,1,1);
273 for(unsigned int ix=0;ix<5;++ix) {
274 normalize(_h_plus_lam.histos()[ix] );
275 normalize(_h_minus_lam.histos()[ix]);
276 std::ostringstream title;
277 title << "/TMP/a_lam_" << ix;
278 Scatter2DPtr sTemp;
279 book(sTemp,title.str());
280 asymm(_h_plus_lam.histos()[ix],_h_minus_lam.histos()[ix],sTemp);
281 pair<double,double> alpha = calcAsymmetry(sTemp,1);
282 h_asym->addPoint(x_val2[ix], alpha.first, make_pair(x_err2[ix],x_err3[ix]),
283 make_pair(alpha.second,alpha.second) );
284 }
285 }
286
287 //@}
288
289
290 /// @name Histograms
291 //@{
292 BinnedHistogram _h_ctheta,_h_plus_cphi,_h_minus_cphi,_h_plus_lam,_h_minus_lam;
293 Histo1DPtr _h_ctheta_large;
294 Histo1DPtr _h_minus_cphi_low,_h_minus_cphi_mid,_h_minus_cphi_high;
295 Histo1DPtr _h_plus_cphi_low ,_h_plus_cphi_mid ,_h_plus_cphi_high ;
296 //@}
297
298
299 };
300
301
302 // The hook for the plugin system
303 RIVET_DECLARE_PLUGIN(ALEPH_1996_I415745);
304
305
306}
|