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