rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

OPAL_2001_I554583

$\tau$ polarization at LEP1
Experiment: OPAL (LEP)
Inspire ID: 554583
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Eur.Phys.J.C 21 (2001) 1-21
Beams: e+ e-
Beam energies: (45.6, 45.6) GeV
Run details:
  • e+ e- > tau+ tau-

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}