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