Rivet analyses referenceOPAL_2001_I554583$\tau$ polarization at LEP1Experiment: OPAL (LEP) Inspire ID: 554583 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
Measurement of the $\tau$ lepton polarization in $e^+e^-\to\tau^+\tau^-$ at the $Z^0$ pole by the OPAL experiment at LEP1. Source code: OPAL_2001_I554583.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/UnstableParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief e+e- > tau+ tau-
11 class OPAL_2001_I554583 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2001_I554583);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // Initialise and register projections
24 declare(Beam(), "Beams");
25 declare(ChargedFinalState(), "FS");
26 declare(UnstableParticles(), "UFS");
27 // book histos
28 const vector<double> edges{-0.9, -0.72, -0.54, -0.36, -0.18, 0., 0.18, 0.36, 0.54, 0.72, 0.9};
29 book(_h_e, edges);
30 book(_h_mu, edges);
31 book(_h_pi, edges);
32 book(_h_rho, edges);
33 for (size_t ix=0; _h_e->numBins(); ++ix) {
34 const string suff = std::to_string(ix);
35 book(_h_e->bin(ix+1), "_h_e_"+suff, 20, -1.0, 1.0);
36 book(_h_mu->bin(ix+1), "_h_mu_"+suff, 20, -1.0, 1.0);
37 book(_h_pi->bin(ix+1), "_h_pi_"+suff, 20, -1.0, 1.0);
38 book(_h_rho->bin(ix+1), "_h_rho_"+suff, 20, -1.0, 1.0);
39 }
40 book(_t_e ,"_t_e " , 20,-1,1);
41 book(_t_mu ,"_t_mu" , 20,-1,1);
42 book(_t_pi ,"_t_pi" , 20,-1,1);
43 book(_t_rho,"_t_rho", 20,-1,1);
44 }
45
46 void findTau(const Particle & p, unsigned int & nprod, Particles& piP,
47 Particles& pi0, Particles& ell, Particles& nu_ell, Particles& nu_tau) {
48 for (const Particle & child : p.children()) {
49 if(child.pid()==PID::ELECTRON || child.pid()==PID::MUON) {
50 ++nprod;
51 ell.push_back(child);
52 }
53 else if(child.pid()==PID::NU_EBAR || child.pid()==PID::NU_MUBAR) {
54 ++nprod;
55 nu_ell.push_back(child);
56 }
57 else if(child.pid()==PID::PIMINUS) {
58 ++nprod;
59 piP.push_back(child);
60 }
61 else if(child.pid()==PID::PI0) {
62 ++nprod;
63 pi0.push_back(child);
64 }
65 else if(child.pid()==PID::NU_TAU) {
66 ++nprod;
67 nu_tau.push_back(child);
68 }
69 else if(child.pid()==PID::GAMMA) {
70 continue;
71 }
72 else if(child.children().empty() || child.pid()==221 || child.pid()==331) {
73 ++nprod;
74 }
75 else {
76 findTau(child,nprod,piP,pi0,ell,nu_ell,nu_tau);
77 }
78 }
79 }
80
81 /// Perform the per-event analysis
82 void analyze(const Event& event) {
83 // require 2 chanrged particles to veto hadronic events
84 if(apply<ChargedFinalState>(event, "FS").particles().size()!=2) vetoEvent;
85 // Get beams and average beam momentum
86 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
87 Vector3 axis;
88 if (beams.first.pid()>0) {
89 axis = beams.first .momentum().p3().unit();
90 }
91 else {
92 axis = beams.second.momentum().p3().unit();
93 }
94 // loop over tau leptons
95 for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==15)) {
96 unsigned int nprod(0);
97 Particles piP, pi0, ell, nu_ell, nu_tau;
98 findTau(p,nprod,piP, pi0, ell, nu_ell, nu_tau);
99 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
100 double cBeam = axis.dot(p.momentum().p3().unit());
101 if (nprod==2 && nu_tau.size()==1 && piP.size()==1) {
102 FourMomentum pPi = boost1.transform(piP[0].momentum());
103 double cTheta = pPi.p3().unit().dot(p.momentum().p3().unit());
104 _h_pi->fill(cBeam, cTheta);
105 _t_pi->fill(cTheta);
106 }
107 else if (nprod==3 && nu_tau.size()==1 && ell.size()==1 && nu_ell.size()==1) {
108 if (ell[0].pid()==PID::ELECTRON) {
109 _h_e->fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
110 _t_e ->fill(2.*ell[0].momentum().t()/sqrtS());
111 }
112 else {
113 _h_mu->fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
114 _t_mu->fill(2.*ell[0].momentum().t()/sqrtS());
115 }
116 }
117 else if(nprod==3 && nu_tau.size()==1 && piP.size()==1&& pi0.size()==1) {
118 FourMomentum pRho = boost1.transform(piP[0].momentum()+pi0[0].momentum());
119 double cTheta = pRho.p3().unit().dot(p.momentum().p3().unit());
120 _h_rho->fill(cBeam, cTheta);
121 _t_rho->fill(cTheta);
122 }
123 }
124 }
125
126 pair<double,double> calcP(Histo1DPtr hist,unsigned int imode) {
127 if(hist->numEntries()==0.) return make_pair(0.,0.);
128 double sum1(0.),sum2(0.);
129 for (const auto& bin : hist->bins() ) {
130 double Oi = bin.sumW();
131 if(Oi==0.) continue;
132 double ai(0.),bi(0.);
133 // tau -> pi/rho nu
134 if(imode==0) {
135 ai = 0.5*(bin.xMax()-bin.xMin());
136 bi = 0.5*ai*(bin.xMax()+bin.xMin());
137 }
138 // lepton mode
139 else {
140 ai = (-5*bin.xMin() + 3*pow(bin.xMin(),3) - pow(bin.xMin(),4) + 5*bin.xMax() - 3*pow(bin.xMax(),3) + pow(bin.xMax(),4))/3.;
141 bi = ( -bin.xMin() + 3*pow(bin.xMin(),3) - 2*pow(bin.xMin(),4) + bin.xMax() - 3*pow(bin.xMax(),3) + 2*pow(bin.xMax(),4))/3.;
142 }
143 double Ei = bin.errW();
144 sum1 += sqr(bi/Ei);
145 sum2 += bi/sqr(Ei)*(Oi-ai);
146 }
147 return make_pair(sum2/sum1,sqrt(1./sum1));
148 }
149
150 /// Normalise histograms etc., after the run
151 void finalize() {
152 Estimate1DPtr _h_P;
153 book(_h_P,2,1,1);
154 Estimate1DPtr _t_P;
155 book(_t_P,1,1,5);
156 for (size_t ix=0; ix < _h_e->numBins()+1; ++ix) {
157 Histo1DPtr& he = ix<10 ? _h_e->bin(ix+1) : _t_e;
158 normalize(he);
159 pair<double,double> P_e = calcP(he, 1);
160 double s1 = P_e.first/sqr(P_e.second);
161 double s2 = 1./sqr(P_e.second);
162 Histo1DPtr& hmu = ix<10 ? _h_mu->bin(ix+1) : _t_mu;
163 normalize(hmu);
164 pair<double,double> P_mu = calcP(hmu,1);
165 s1 += P_mu.first/sqr(P_mu.second);
166 s2 += 1./sqr(P_mu.second);
167 Histo1DPtr& hpi = ix<10 ? _h_pi->bin(ix+1) : _t_pi;
168 normalize(hpi);
169 pair<double,double> P_pi = calcP(hpi,0);
170 s1 += P_pi.first/sqr(P_pi.second);
171 s2 += 1./sqr(P_pi.second);
172 Histo1DPtr& hrho = ix<10 ? _h_rho->bin(ix+1) : _t_rho;
173 normalize(hrho);
174 pair<double,double> P_rho = calcP(hrho,0);
175 s1 += P_rho.first/sqr(P_rho.second);
176 s2 += 1./sqr(P_rho.second);
177 // average
178 _h_P->bin(ix+1).set(s1/s2, sqrt(1./s2));
179 }
180 }
181
182 /// @}
183
184
185 /// @name Histograms
186 /// @{
187 Histo1DGroupPtr _h_e,_h_mu,_h_pi,_h_rho;
188 Histo1DPtr _t_e,_t_mu,_t_pi,_t_rho;
189 /// @}
190
191
192 };
193
194
195 RIVET_DECLARE_PLUGIN(OPAL_2001_I554583);
196
197}
|