Rivet analyses referenceALEPH_2001_I555653$\tau$ polarization at LEP1Experiment: ALEPH (LEP) Inspire ID: 555653 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 ALEPH experiment at LEP1. Source code: ALEPH_2001_I555653.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 ALEPH_2001_I555653 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_2001_I555653);
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 hists
29 double x = -0.9, step=0.2;
30 for(unsigned int ix=0;ix<9;++ix) {
31 Histo1DPtr temp;
32 std::ostringstream title1;
33 title1 << "_h_e_" << ix;
34 book(temp,title1.str(), 20,-1,1);
35 _h_e .add(x,x+step,temp);
36 std::ostringstream title2;
37 title2 << "_h_mu_" << ix;
38 book(temp,title2.str(), 20,-1,1);
39 _h_mu.add(x,x+step,temp);
40 std::ostringstream title3;
41 title3 << "_h_pi_" << ix;
42 book(temp,title3.str(), 20,-1,1);
43 _h_pi.add(x,x+step,temp);
44 std::ostringstream title4;
45 title4 << "_h_rho_" << ix;
46 book(temp,title4.str(), 20,-1,1);
47 _h_rho.add(x,x+step,temp);
48 x+=step;
49 }
50 }
51
52 void findTau(const Particle & p, unsigned int & nprod,
53 Particles & piP,Particles & pi0, Particles & ell, Particles & nu_ell,
54 Particles & nu_tau) {
55 for(const Particle & child : p.children()) {
56 if(child.pid()==PID::ELECTRON || child.pid()==PID::MUON) {
57 ++nprod;
58 ell.push_back(child);
59 }
60 else if(child.pid()==PID::NU_EBAR || child.pid()==PID::NU_MUBAR) {
61 ++nprod;
62 nu_ell.push_back(child);
63 }
64 else if(child.pid()==PID::PIMINUS) {
65 ++nprod;
66 piP.push_back(child);
67 }
68 else if(child.pid()==PID::PI0) {
69 ++nprod;
70 pi0.push_back(child);
71 }
72 else if(child.pid()==PID::NU_TAU) {
73 ++nprod;
74 nu_tau.push_back(child);
75 }
76 else if(child.pid()==PID::GAMMA)
77 continue;
78 else if(child.children().empty() || child.pid()==221 || child.pid()==331) {
79 ++nprod;
80 }
81 else {
82 findTau(child,nprod,piP,pi0,ell,nu_ell,nu_tau);
83 }
84 }
85 }
86
87 /// Perform the per-event analysis
88 void analyze(const Event& event) {
89 // require 2 chanrged particles to veto hadronic events
90 if(apply<ChargedFinalState>(event, "FS").particles().size()!=2) vetoEvent;
91 // Get beams and average beam momentum
92 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
93 Vector3 axis;
94 if(beams.first.pid()>0)
95 axis = beams.first .momentum().p3().unit();
96 else
97 axis = beams.second.momentum().p3().unit();
98 // loop over tau leptons
99 for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==15)) {
100 unsigned int nprod(0);
101 Particles piP, pi0, ell, nu_ell, nu_tau;
102 findTau(p,nprod,piP, pi0, ell, nu_ell, nu_tau);
103 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
104 double cBeam = axis.dot(p.momentum().p3().unit());
105 if(nprod==2 && nu_tau.size()==1 && piP.size()==1) {
106 FourMomentum pPi = boost1.transform(piP[0].momentum());
107 double cTheta = pPi.p3().unit().dot(p.momentum().p3().unit());
108 _h_pi. fill(cBeam,cTheta);
109 }
110 else if(nprod==3 && nu_tau.size()==1 && ell.size()==1 && nu_ell.size()==1) {
111 if(ell[0].pid()==PID::ELECTRON) {
112 _h_e . fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
113 }
114 else {
115 _h_mu. fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
116 }
117 }
118 else if(nprod==3 && nu_tau.size()==1 && piP.size()==1&& pi0.size()==1) {
119 FourMomentum pRho = boost1.transform(piP[0].momentum()+pi0[0].momentum());
120 double cTheta = pRho.p3().unit().dot(p.momentum().p3().unit());
121 _h_rho. fill(cBeam,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 (auto bin : hist->bins() ) {
130 double Oi = bin.area();
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.areaErr();
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 double x = -0.8, step=0.2;
153 Scatter2DPtr _h_P;
154 book(_h_P,2,1,1);
155 for(unsigned int ix=0;ix<9;++ix) {
156 normalize(_h_e.histos()[ix]);
157 pair<double,double> P_e = calcP(_h_e.histos()[ix],1);
158 double s1 = P_e.first/sqr(P_e.second);
159 double s2 = 1./sqr(P_e.second);
160 normalize(_h_mu.histos()[ix]);
161 pair<double,double> P_mu = calcP(_h_mu.histos()[ix],1);
162 s1 += P_mu.first/sqr(P_mu.second);
163 s2 += 1./sqr(P_mu.second);
164 normalize(_h_pi.histos()[ix]);
165 pair<double,double> P_pi = calcP(_h_pi.histos()[ix],0);
166 s1 += P_pi.first/sqr(P_pi.second);
167 s2 += 1./sqr(P_pi.second);
168 normalize(_h_rho.histos()[ix]);
169 pair<double,double> P_rho = calcP(_h_rho.histos()[ix],0);
170 s1 += P_rho.first/sqr(P_rho.second);
171 s2 += 1./sqr(P_rho.second);
172 P_rho.first /=0.46;
173 P_rho.second /=0.46;
174 // average
175 pair<double,double> P_aver = make_pair(s1/s2,sqrt(1./s2));
176 _h_P->addPoint(x,P_aver.first, make_pair(0.5*step,0.5*step),
177 make_pair(P_aver.second,P_aver.second));
178 x+=step;
179 }
180 }
181
182 ///@}
183
184
185 /// @name Histograms
186 ///@{
187 BinnedHistogram _h_e,_h_mu,_h_pi,_h_rho;
188 ///@}
189
190
191 };
192
193
194 RIVET_DECLARE_PLUGIN(ALEPH_2001_I555653);
195
196}
|