Rivet analyses referenceOPAL_1997_I440103Polarization of $\phi$, $D^{*+}$ and $B^*$ mesons at LEP1Experiment: OPAL (LEP) Inspire ID: 440103 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
Measurement of the polarization of $\phi$, $D^{*+}$ and $B^*$ mesons at LEP1. Source code: OPAL_1997_I440103.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/ChargedFinalState.hh"
6#include "Rivet/Projections/UnstableParticles.hh"
7#include "Rivet/Projections/Thrust.hh"
8
9#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
10#include "Rivet/Projections/InitialQuarks.hh"
11
12namespace Rivet {
13
14
15 /// @brief phi D* and B* polarization
16 class OPAL_1997_I440103 : public Analysis {
17 public:
18
19 /// Constructor
20 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1997_I440103);
21
22
23 /// @name Analysis methods
24 /// @{
25
26 /// Book histograms and initialise projections before the run
27 void init() {
28 // Initialise and register projections
29 declare(Beam(), "Beams");
30 declare(Thrust(FinalState()), "Thrust");
31 declare(ChargedFinalState(), "FS");
32 declare(InitialQuarks(), "IQF");
33 declare(UnstableParticles(), "UFS" );
34
35 // Book histograms
36 // B*
37 book(_h_B , 8,1,1);
38 book(_h_B2, "/TMP/c_theta_B", 20, -1.,1.);
39 // phi
40 book(_h_phi_ctheta , 5,1,1);
41 book(_h_phi_ctheta2, "/TMP/c_theta_phi2", 20, -1.,1. );
42 book(_h_phi_ctheta3, "/TMP/c_theta_phi3", 20, -1.,1. );
43 book(_h_phi_ctheta4, "/TMP/c_theta_phi4", 20, -1.,1. );
44 book(_h_phi_alpha , 5,1,2);
45 book(_h_phi_alpha2 , "/TMP/alpha_phi2", 20, 0.,0.5*M_PI);
46 book(_h_phi_alpha3 , "/TMP/alpha_phi3", 20, 0.,0.5*M_PI);
47 book(_h_phi_alpha4 , "/TMP/alpha_phi4", 20, 0.,0.5*M_PI);
48 book(_h_phi_beta , 5,1,3);
49 book(_h_phi_beta2 , "/TMP/beta_phi2", 20, 0.,0.5*M_PI );
50 book(_h_phi_beta3 , "/TMP/beta_phi3", 20, 0.,0.5*M_PI );
51 book(_h_phi_beta4 , "/TMP/beta_phi4", 20, 0.,0.5*M_PI );
52 book(_c_phi_cos_plus , "/TMP/c_phi_cos_plus1");
53 book(_c_phi_cos_neg , "/TMP/c_phi_cos_neg1" );
54 book(_c_phi_sin_plus , "/TMP/c_phi_sin_plus1");
55 book(_c_phi_sin_neg , "/TMP/c_phi_sin_neg1" );
56 book(_c_phi_cos_plus2, "/TMP/c_phi_cos_plus2");
57 book(_c_phi_cos_neg2 , "/TMP/c_phi_cos_neg2" );
58 book(_c_phi_sin_plus2, "/TMP/c_phi_sin_plus2");
59 book(_c_phi_sin_neg2 , "/TMP/c_phi_sin_neg2" );
60 book(_c_phi_cos_plus3, "/TMP/c_phi_cos_plus3");
61 book(_c_phi_cos_neg3 , "/TMP/c_phi_cos_neg3" );
62 book(_c_phi_sin_plus3, "/TMP/c_phi_sin_plus3");
63 book(_c_phi_sin_neg3 , "/TMP/c_phi_sin_neg3" );
64 // D*
65 book(_h_DS_ctheta , 6,1,1);
66 book(_h_DS_ctheta2, "/TMP/c_theta_DS2", 20, -1.,1. );
67 book(_h_DS_alpha , 7,1,1);
68 book(_h_DS_alpha2 , "/TMP/alpha_DS2", 20, 0.,0.5*M_PI);
69 }
70
71
72 /// Perform the per-event analysis
73 void analyze(const Event& event) {
74 // First, veto on leptonic events by requiring at least 4 charged FS particles
75 const FinalState& fs = apply<FinalState>(event, "FS");
76 const size_t numParticles = fs.particles().size();
77
78 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
79 if (numParticles < 2) {
80 MSG_DEBUG("Failed leptonic event cut");
81 vetoEvent;
82 }
83 MSG_DEBUG("Passed leptonic event cut");
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 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
90 Vector3 axis;
91 if(beams.first.pid()>0)
92 axis = beams.first.momentum().p3().unit();
93 else
94 axis = beams.second.momentum().p3().unit();
95 // thrust, to define an axis
96 const Thrust& thrust = apply<Thrust>(event, "Thrust");
97
98 int flavour = 0;
99 const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
100
101 // If we only have two quarks (qqbar), just take the flavour.
102 // If we have more than two quarks, look for the highest energetic q-qbar pair.
103 /// @todo Yuck... does this *really* have to be quark-based?!?
104 if (iqf.particles().size() == 2) {
105 flavour = iqf.particles().front().abspid();
106 } else {
107 map<int, double> quarkmap;
108 for (const Particle& p : iqf.particles()) {
109 if (quarkmap[p.pid()] < p.E()) {
110 quarkmap[p.pid()] = p.E();
111 }
112 }
113 double maxenergy = 0.;
114 for (int i = 1; i <= 5; ++i) {
115 if (quarkmap[i]+quarkmap[-i] > maxenergy) {
116 flavour = i;
117 }
118 }
119 }
120
121 // loop over the particles
122 for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==513 or Cuts::abspid==523 or
123 Cuts::pid==333 or Cuts::abspid==413)) {
124 int sign = p.pid()/p.abspid();
125 Particle decay;
126 if(p.children().size()!=2) continue;
127 // B*
128 if(p.abspid()==513 or p.abspid()==523) {
129 int mid = p.abspid()-2;
130 if(p.children()[0].pid()==sign*mid &&
131 p.children()[1].pid()==22) {
132 decay = p.children()[1];
133 }
134 else if(p.children()[1].pid()==sign*mid &&
135 p.children()[0].pid()==22) {
136 decay = p.children()[0];
137 }
138 else {
139 continue;
140 }
141 }
142 // phi
143 else if(p.pid()==333) {
144 // cut x_E > 0.7
145 double xE = p.momentum().E()/meanBeamMom;
146 if(xE<0.7) continue;
147 if(p.children()[0].pid()== 321 &&
148 p.children()[1].pid()==-321) {
149 decay = p.children()[0];
150 }
151 else if(p.children()[1].pid()== 321 &&
152 p.children()[0].pid()==-321) {
153 decay = p.children()[1];
154 }
155 else {
156 continue;
157 }
158 }
159 // D*
160 else if(p.abspid()==413) {
161 double xE = p.momentum().E()/meanBeamMom;
162 if(xE<0.5 || flavour!=4) continue;
163 if(p.children()[0].pid()==sign*421 &&
164 p.children()[1].pid()==sign*211) {
165 decay = p.children()[1];
166 }
167 else if(p.children()[1].pid()==sign*421 &&
168 p.children()[0].pid()==sign*211) {
169 decay = p.children()[0];
170 }
171 else {
172 continue;
173 }
174 }
175 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
176 Vector3 e1z = p.p3().unit();
177 FourMomentum pp = boost.transform(decay.momentum());
178 Vector3 axis1 = boost.transform(decay.momentum()).p3().unit();
179 double ctheta = e1z.dot(axis1);
180 if(p.abspid()==513 or p.abspid()==523) {
181 _h_B ->fill(ctheta);
182 _h_B2->fill(ctheta);
183 }
184 // D*
185 else if(p.abspid()==413) {
186 // y and z axis
187 Vector3 e1y = e1z.cross(axis).unit();
188 Vector3 e1x = e1y.cross(e1z).unit();
189 // helicity beam axis, all phis
190 // cos theta_H
191 _h_DS_ctheta ->fill(ctheta);
192 _h_DS_ctheta2->fill(ctheta);
193 // alpha
194 double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
195 double alpha = abs(abs(phi)-0.5*M_PI);
196 _h_DS_alpha ->fill(alpha);
197 _h_DS_alpha2 ->fill(alpha);
198 }
199 else if(p.pid()==333) {
200 // y and z axis
201 Vector3 e1y = e1z.cross(axis).unit();
202 Vector3 e1x = e1y.cross(e1z).unit();
203 // helicity beam axis, all phis
204 // cos theta_H
205 _h_phi_ctheta->fill(abs(ctheta));
206 _h_phi_ctheta2->fill(ctheta);
207 // alpha and beta
208 double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
209 double alpha = abs(abs(phi)-0.5*M_PI);
210 double beta = abs(abs(phi+0.25*M_PI)-0.5*M_PI);
211 _h_phi_alpha ->fill(alpha);
212 _h_phi_alpha2 ->fill(alpha);
213 _h_phi_beta ->fill( beta);
214 _h_phi_beta2 ->fill( beta);
215 /// counters for asymmetries
216 double sin2H = 2.*ctheta*sqrt(1.-sqr(ctheta));
217 if(sin2H*cos(phi)>0.)
218 _c_phi_cos_plus->fill();
219 else
220 _c_phi_cos_neg->fill();
221 if(sin2H*sin(phi)>0.)
222 _c_phi_sin_plus->fill();
223 else
224 _c_phi_sin_neg->fill();
225 // whether or not is a primary hadron
226 Particle parent = p.parents()[0];
227 if(parent.children().size()==1 && parent.abspid()==p.abspid())
228 parent = parent.parents()[0];
229 bool primary = !PID::isHadron(parent.pid());
230 if (primary) {
231 // cos theta_H
232 _h_phi_ctheta3->fill(ctheta);
233 // alpha and beta
234 _h_phi_alpha3 ->fill(alpha);
235 _h_phi_beta3 ->fill( beta);
236 /// counters for asymmetries
237 if(sin2H*cos(phi)>0.)
238 _c_phi_cos_plus2->fill();
239 else
240 _c_phi_cos_neg2->fill();
241 if(sin2H*sin(phi)>0.)
242 _c_phi_sin_plus2->fill();
243 else
244 _c_phi_sin_neg2->fill();
245 }
246 // pT w.r.t thrust axis
247 double pT = sqrt(sqr(thrust.thrustMajorAxis().dot(p.momentum().p3()))+ sqr(thrust.thrustMinorAxis().dot(p.momentum().p3())));
248 // helicity-quark frame
249 if(pT>1.2) {
250 // cos theta H
251 _h_phi_ctheta4->fill(ctheta);
252 Vector3 axis2;
253 if(p.momentum().p3().dot(thrust.thrustAxis())>=0.) {
254 axis2 = thrust.thrustAxis();
255 }
256 else {
257 axis2 =-thrust.thrustAxis();
258 }
259 Vector3 e2y = e1z.cross(axis2).unit();
260 Vector3 e2x = e2y.cross(e1z).unit();
261 // alpha and beta
262 double phi = atan2(e2y.dot(axis1),e2x.dot(axis1));
263 double alpha = abs(abs(phi)-0.5*M_PI);
264 double beta = abs(abs(phi+0.25*M_PI)-0.5*M_PI);
265 _h_phi_alpha4 ->fill(alpha);
266 _h_phi_beta4 ->fill( beta);
267 /// counters for asymmetries
268 double sin2H = 2.*ctheta*sqrt(1.-sqr(ctheta));
269 if(sin2H*cos(phi)>0.)
270 _c_phi_cos_plus3->fill();
271 else
272 _c_phi_cos_neg3->fill();
273 if(sin2H*sin(phi)>0.)
274 _c_phi_sin_plus3->fill();
275 else
276 _c_phi_sin_neg3->fill();
277
278 }
279 }
280 }
281 }
282
283 pair<double,double> calcRho(Histo1DPtr hist,unsigned int mode) {
284 if(hist->numEntries()==0.) return make_pair(0.,0.);
285 double sum1(0.),sum2(0.);
286 for (const auto& bin : hist->bins() ) {
287 double Oi = bin.sumW();
288 if (Oi==0.) continue;
289 double ai(0.),bi(0.);
290 if(mode==0) {
291 ai = 0.25*( -bin.xMin()*(3.-sqr(bin.xMin())) + bin.xMax()*(3.-sqr(bin.xMax())));
292 bi =-0.75*( -bin.xMin()*(1.-sqr(bin.xMin())) + bin.xMax()*(1.-sqr(bin.xMax())));
293 }
294 else if(mode==1) {
295 ai = 0.125*( -bin.xMin()*(3.+sqr(bin.xMin())) + bin.xMax()*(3.+sqr(bin.xMax())));
296 bi = 0.375*( -bin.xMin()*(1.-sqr(bin.xMin())) + bin.xMax()*(1.-sqr(bin.xMax())));
297 }
298 else if(mode==2) {
299 ai = -2.*(bin.xMin()-bin.xMax())/M_PI;
300 bi = -2.*(sin(2.*bin.xMin())-sin(2.*bin.xMax()))/M_PI;
301 }
302 double Ei = bin.errW();
303 sum1 += sqr(bi/Ei);
304 sum2 += bi/sqr(Ei)*(Oi-ai);
305 }
306 return make_pair(sum2/sum1,sqrt(1./sum1));
307 }
308
309 /// Normalise histograms etc., after the run
310 void finalize() {
311 // B*
312 normalize(_h_B,1.,false);
313 normalize(_h_B2);
314 pair<double,double> rho = calcRho(_h_B2,1);
315 BinnedEstimatePtr<string> h_rhoB;
316 book(h_rhoB,4,1,1);
317 h_rhoB->bin(1).set(rho.first, rho.second);
318 // D*
319 normalize(_h_DS_ctheta );
320 normalize(_h_DS_ctheta2);
321 rho = calcRho(_h_DS_ctheta2,1);
322 BinnedEstimatePtr<string> h_rhoD;
323 book(h_rhoD,3,1,1);
324 h_rhoD->bin(1).set(rho.first, rho.second);
325 normalize(_h_DS_alpha );
326 normalize(_h_DS_alpha2);
327 BinnedEstimatePtr<string> h_reRho_D;
328 book(h_reRho_D,3,1,2);
329 rho = calcRho(_h_DS_alpha2,2);
330 h_reRho_D->bin(1).set(rho.first, rho.second);
331 // phi
332 // rho00
333 normalize(_h_phi_ctheta );
334 normalize(_h_phi_ctheta2);
335 normalize(_h_phi_ctheta3);
336 normalize(_h_phi_ctheta4);
337 BinnedEstimatePtr<string> hrho_phi;
338 book(hrho_phi,1,1,1);
339 rho = calcRho(_h_phi_ctheta2,0);
340 hrho_phi->bin(1).set(rho.first, rho.second);
341 rho = calcRho(_h_phi_ctheta3,0);
342 hrho_phi->bin(2).set(rho.first, rho.second);
343 rho = calcRho(_h_phi_ctheta4,0);
344 hrho_phi->bin(3).set(rho.first, rho.second);
345 // Re rho
346 normalize(_h_phi_alpha );
347 normalize(_h_phi_alpha2);
348 normalize(_h_phi_alpha3);
349 normalize(_h_phi_alpha4);
350 BinnedEstimatePtr<string> hreRho_phi;
351 book(hreRho_phi,1,1,2);
352 rho = calcRho(_h_phi_alpha2,2);
353 hreRho_phi->bin(1).set(rho.first, rho.second);
354 rho = calcRho(_h_phi_alpha3,2);
355 hreRho_phi->bin(2).set(rho.first, rho.second);
356 rho = calcRho(_h_phi_alpha4,2);
357 hreRho_phi->bin(3).set(rho.first, rho.second);
358 // Im rho
359 normalize(_h_phi_beta );
360 normalize(_h_phi_beta2);
361 normalize(_h_phi_beta3);
362 normalize(_h_phi_beta4);
363 BinnedEstimatePtr<string> himRho_phi;
364 book(himRho_phi,1,1,3);
365 rho = calcRho(_h_phi_beta2,2);
366 himRho_phi->bin(1).set(rho.first, rho.second);
367 rho = calcRho(_h_phi_beta3,2);
368 himRho_phi->bin(2).set(rho.first, rho.second);
369 rho = calcRho(_h_phi_beta4,2);
370 himRho_phi->bin(3).set(rho.first, rho.second);
371 // real diff
372 Estimate0D temp = ((*_c_phi_cos_plus-*_c_phi_cos_neg)/(*_c_phi_cos_plus+*_c_phi_cos_neg));
373 Estimate0D temp2 = ((*_c_phi_cos_plus2-*_c_phi_cos_neg2)/(*_c_phi_cos_plus2+*_c_phi_cos_neg2));
374 Estimate0D temp3 = ((*_c_phi_cos_plus3-*_c_phi_cos_neg3)/(*_c_phi_cos_plus3+*_c_phi_cos_neg3));
375 BinnedEstimatePtr<string> hreDiff_phi;
376 book(hreDiff_phi,1,1,4);
377 hreDiff_phi->bin(1) = temp;
378 hreDiff_phi->bin(2) = temp2;
379 hreDiff_phi->bin(3) = temp3;
380 // im diff
381 temp = ((*_c_phi_sin_plus-*_c_phi_sin_neg)/(*_c_phi_sin_plus+*_c_phi_sin_neg));
382 temp2 = ((*_c_phi_sin_plus2-*_c_phi_sin_neg2)/(*_c_phi_sin_plus2+*_c_phi_sin_neg2));
383 temp3 = ((*_c_phi_sin_plus3-*_c_phi_sin_neg3)/(*_c_phi_sin_plus3+*_c_phi_sin_neg3));
384 BinnedEstimatePtr<string> himDiff_phi;
385 book(himDiff_phi,1,1,5);
386 himDiff_phi->bin(1) = temp;
387 himDiff_phi->bin(2) = temp2;
388 himDiff_phi->bin(3) = temp3;
389 }
390
391 /// @}
392
393
394 /// @name Histograms
395 /// @{
396 Histo1DPtr _h_B,_h_B2;
397 Histo1DPtr _h_phi_ctheta, _h_phi_ctheta2, _h_phi_ctheta3, _h_phi_ctheta4;
398 Histo1DPtr _h_phi_alpha , _h_phi_alpha2 , _h_phi_alpha3 , _h_phi_alpha4 ;
399 Histo1DPtr _h_phi_beta , _h_phi_beta2 , _h_phi_beta3 , _h_phi_beta4 ;
400 CounterPtr _c_phi_cos_plus, _c_phi_cos_neg, _c_phi_cos_plus2, _c_phi_cos_neg2, _c_phi_cos_plus3, _c_phi_cos_neg3;
401 CounterPtr _c_phi_sin_plus, _c_phi_sin_neg, _c_phi_sin_plus2, _c_phi_sin_neg2, _c_phi_sin_plus3, _c_phi_sin_neg3;
402 Histo1DPtr _h_DS_ctheta, _h_DS_ctheta2;
403 Histo1DPtr _h_DS_alpha , _h_DS_alpha2 ;
404 /// @}
405
406
407 };
408
409
410 RIVET_DECLARE_PLUGIN(OPAL_1997_I440103);
411
412
413}
|
