rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

DELPHI_2000_I511443

$\tau$ polarization at LEP1
Experiment: DELPHI (LEP)
Inspire ID: 511443
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Eur.Phys.J.C 14 (2000) 585-611
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 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}