Rivet analyses referenceDELPHI_2000_I511443$\tau$ polarization at LEP1Experiment: DELPHI (LEP) Inspire ID: 511443 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 DELPHI experiment at LEP1. Source code: DELPHI_2000_I511443.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 DELPHI_2000_I511443 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2000_I511443);
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> bins{-0.94,-0.732,-0.488,-0.2440,0.,0.244,0.488,0.732,0.94};
29 book(_h_e, bins);
30 book(_h_mu, bins);
31 book(_h_pi, bins);
32 book(_h_rho, bins);
33 for (size_t ix = 0; ix < _h_e->numBins(); ++ix) {
34 const string suff = 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 }
41
42 void findTau(const Particle & p, unsigned int & nprod, Particles& piP,
43 Particles& pi0, Particles& ell, Particles& nu_ell, Particles& nu_tau) {
44 for (const Particle& child : p.children()) {
45 if (child.pid()==PID::ELECTRON || child.pid()==PID::MUON) {
46 ++nprod;
47 ell.push_back(child);
48 }
49 else if (child.pid()==PID::NU_EBAR || child.pid()==PID::NU_MUBAR) {
50 ++nprod;
51 nu_ell.push_back(child);
52 }
53 else if (child.pid()==PID::PIMINUS) {
54 ++nprod;
55 piP.push_back(child);
56 }
57 else if (child.pid()==PID::PI0) {
58 ++nprod;
59 pi0.push_back(child);
60 }
61 else if (child.pid()==PID::NU_TAU) {
62 ++nprod;
63 nu_tau.push_back(child);
64 }
65 else if (child.pid()==PID::GAMMA) continue;
66 else if (child.children().empty() || child.pid()==221 || child.pid()==331) {
67 ++nprod;
68 }
69 else {
70 findTau(child,nprod,piP,pi0,ell,nu_ell,nu_tau);
71 }
72 }
73 }
74
75 /// Perform the per-event analysis
76 void analyze(const Event& event) {
77 // require 2 chanrged particles to veto hadronic events
78 if (apply<ChargedFinalState>(event, "FS").particles().size()!=2) vetoEvent;
79 // Get beams and average beam momentum
80 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
81 Vector3 axis;
82 if (beams.first.pid()>0)
83 axis = beams.first .momentum().p3().unit();
84 else
85 axis = beams.second.momentum().p3().unit();
86 // loop over tau leptons
87 for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==15)) {
88 unsigned int nprod(0);
89 Particles piP, pi0, ell, nu_ell, nu_tau;
90 findTau(p,nprod,piP, pi0, ell, nu_ell, nu_tau);
91 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
92 double cBeam = axis.dot(p.momentum().p3().unit());
93 if (nprod==2 && nu_tau.size()==1 && piP.size()==1) {
94 FourMomentum pPi = boost1.transform(piP[0].momentum());
95 double cTheta = pPi.p3().unit().dot(p.momentum().p3().unit());
96 _h_pi->fill(cBeam,cTheta);
97 }
98 else if (nprod==3 && nu_tau.size()==1 && ell.size()==1 && nu_ell.size()==1) {
99 if (ell[0].pid()==PID::ELECTRON) {
100 _h_e->fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
101 }
102 else {
103 _h_mu->fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
104 }
105 }
106 else if (nprod==3 && nu_tau.size()==1 && piP.size()==1&& pi0.size()==1) {
107 FourMomentum pRho = boost1.transform(piP[0].momentum()+pi0[0].momentum());
108 double cTheta = pRho.p3().unit().dot(p.momentum().p3().unit());
109 _h_rho->fill(cBeam,cTheta);
110 }
111 }
112 }
113
114 pair<double,double> calcP(Histo1DPtr hist,unsigned int imode) {
115 if(hist->numEntries()==0.) return make_pair(0.,0.);
116 double sum1(0.),sum2(0.);
117 for (const auto& bin : hist->bins() ) {
118 double Oi = bin.sumW();
119 if(Oi==0.) continue;
120 double ai(0.),bi(0.);
121 // tau -> pi/rho nu
122 if(imode==0) {
123 ai = 0.5*(bin.xMax()-bin.xMin());
124 bi = 0.5*ai*(bin.xMax()+bin.xMin());
125 }
126 // lepton mode
127 else {
128 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.;
129 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.;
130 }
131 double Ei = bin.errW();
132 sum1 += sqr(bi/Ei);
133 sum2 += bi/sqr(Ei)*(Oi-ai);
134 }
135 return make_pair(sum2/sum1,sqrt(1./sum1));
136 }
137
138 /// Normalise histograms etc., after the run
139 void finalize() {
140 Estimate1DPtr _h_P;
141 book(_h_P,1,1,1);
142 for (size_t ix=0; ix < _h_e->numBins(); ++ix) {
143 normalize(_h_e->bin(ix+1));
144 normalize(_h_mu->bin(ix+1));
145 normalize(_h_pi->bin(ix+1));
146 normalize(_h_rho->bin(ix+1));
147 pair<double,double> P_e = calcP(_h_e->bin(ix+1), 1);
148 double s1 = P_e.first/sqr(P_e.second);
149 double s2 = 1./sqr(P_e.second);
150 pair<double,double> P_mu = calcP(_h_mu->bin(ix+1), 1);
151 s1 += P_mu.first/sqr(P_mu.second);
152 s2 += 1./sqr(P_mu.second);
153 pair<double,double> P_pi = calcP(_h_pi->bin(ix+1), 0);
154 s1 += P_pi.first/sqr(P_pi.second);
155 s2 += 1./sqr(P_pi.second);
156 pair<double,double> P_rho = calcP(_h_rho->bin(ix+1), 0);
157 P_rho.first /=0.46;
158 P_rho.second /=0.46;
159 s1 += P_rho.first/sqr(P_rho.second);
160 s2 += 1./sqr(P_rho.second);
161 // average
162 _h_P->bin(ix+1).set(s1/s2, sqrt(1./s2));
163 }
164 }
165
166 /// @}
167
168
169 /// @name Histograms
170 /// @{
171 Histo1DGroupPtr _h_e,_h_mu,_h_pi,_h_rho;
172 /// @}
173
174
175 };
176
177
178 RIVET_DECLARE_PLUGIN(DELPHI_2000_I511443);
179
180}
|
