rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

OPAL_2000_I502750

Polarization of $\rho^\pm$ and $\omega^0$ mesons at LEP1
Experiment: OPAL (LEP)
Inspire ID: 502750
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Eur.Phys.J. C16 (2000) 61-70
Beams: e+ e-
Beam energies: (45.6, 45.6) GeV
Run details:
  • Hadronic Z decay events generated on the Z pole ($\sqrt{s} = 91.2$ GeV)

The measurement of the polarization of $\rho^\pm$ and $\omega^0$ mesons at LEP1 by the OPAL experiment.

Source code: OPAL_2000_I502750.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#include "Rivet/Tools/BinnedHistogram.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief rho+/- and omega polarization
 12  class OPAL_2000_I502750 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2000_I502750);
 17
 18
 19    /// @name Analysis methods
 20    //@{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Initialise and register projections
 26      declare(Beam(), "Beams");
 27      declare(ChargedFinalState(), "FS");
 28      declare(UnstableParticles(), "UFS");
 29
 30      // Book histograms
 31      {Histo1DPtr temp; _h_ctheta_rho  .add(0.025,0.05,book(temp, "ctheta_rho_0",20,-1.,1.));}
 32      {Histo1DPtr temp; _h_ctheta_rho  .add(0.05 ,0.1 ,book(temp, "ctheta_rho_1",20,-1.,1.));}
 33      {Histo1DPtr temp; _h_ctheta_rho  .add(0.1  ,0.15,book(temp, "ctheta_rho_2",20,-1.,1.));}
 34      {Histo1DPtr temp; _h_ctheta_rho  .add(0.15 ,0.3 ,book(temp, "ctheta_rho_3",20,-1.,1.));}
 35      {Histo1DPtr temp; _h_ctheta_rho  .add(0.3  ,0.6 ,book(temp, "ctheta_rho_4",20,-1.,1.));}
 36      {Histo1DPtr temp; _h_ctheta_omega.add(0.025,0.05,book(temp, "ctheta_omega_0",20,-1.,1.));}
 37      {Histo1DPtr temp; _h_ctheta_omega.add(0.05 ,0.1 ,book(temp, "ctheta_omega_1",20,-1.,1.));}
 38      {Histo1DPtr temp; _h_ctheta_omega.add(0.1  ,0.15,book(temp, "ctheta_omega_2",20,-1.,1.));}
 39      {Histo1DPtr temp; _h_ctheta_omega.add(0.15 ,0.3 ,book(temp, "ctheta_omega_3",20,-1.,1.));}
 40      {Histo1DPtr temp; _h_ctheta_omega.add(0.3  ,0.6 ,book(temp, "ctheta_omega_4",20,-1.,1.));}
 41      book(_h_ctheta_omega_all, "ctheta_omega_all",20,-1.,1.);
 42    }
 43  
 44    pair<double,double> calcRho(Histo1DPtr hist) {
 45      if(hist->numEntries()==0.) return make_pair(0.,0.);
 46      double sum1(0.),sum2(0.);
 47      for (auto bin : hist->bins() ) {
 48	double Oi = bin.area();
 49	if(Oi==0.) continue;
 50	double ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
 51	double bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
 52	double Ei = bin.areaErr();
 53	sum1 += sqr(bi/Ei);
 54	sum2 += bi/sqr(Ei)*(Oi-ai);
 55      }
 56      return make_pair(sum2/sum1,sqrt(1./sum1));
 57    }
 58
 59    bool findOmegaDecay(Particle omega,Particles & pi0, Particles & pip, Particles & pim) {
 60      for(const Particle & child : omega.children()) {
 61	if(child.pid()==211)
 62	  pip.push_back(child);
 63	else if(child.pid()==-211)
 64	  pim.push_back(child);
 65	else if(child.pid()==111)
 66	  pi0.push_back(child);
 67	else if(!child.children().empty()) {
 68	  if(!findOmegaDecay(child,pi0,pip,pim)) return false;
 69	}
 70	else
 71	  return false;
 72      }
 73      return true;
 74    }
 75    
 76    /// Perform the per-event analysis
 77    void analyze(const Event& event) {
 78      // First, veto on leptonic events by requiring at least 4 charged FS particles
 79      const FinalState& fs = apply<FinalState>(event, "FS");
 80      const size_t numParticles = fs.particles().size();
 81
 82      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 83      if (numParticles < 2) {
 84        MSG_DEBUG("Failed leptonic event cut");
 85        vetoEvent;
 86      }
 87      MSG_DEBUG("Passed leptonic event cut");
 88      // Get beams and average beam momentum
 89      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 90      const double meanBeamMom = ( beams.first.p3().mod() +
 91                                   beams.second.p3().mod() ) / 2.0;
 92      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 93      // loop over rho and omega mesons
 94      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 95      for (const Particle& p : ufs.particles(Cuts::abspid==213 || Cuts::abspid==223)) {
 96	double xE = p.momentum().t()/meanBeamMom;
 97	Vector3 e1z = p.momentum().p3().unit();
 98	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 99	if(p.abspid()==213) {
100	  if(p.children().size()!=2) continue;
101	  int sign = p.pid()/213;
102	  Particle pion;
103	  if(p.children()[0].pid()==sign*211 && p.children()[1].pid()==111) {
104	    pion = p.children()[0];
105	  }
106	  else if(p.children()[1].pid()==sign*211 && p.children()[0].pid()==111) {
107	    pion = p.children()[1];
108	  }
109	  else
110	    continue;
111	  Vector3 axis1 = boost.transform(pion.momentum()).p3().unit();
112	  double ctheta = e1z.dot(axis1);
113	  _h_ctheta_rho.fill(xE,ctheta);
114	}
115	else {
116	  Particles pi0,pip,pim;
117	  bool three_pi = findOmegaDecay(p,pi0,pip,pim);
118	  if(!three_pi || pi0.size()!=1 || pip.size()!=1 || pim.size()!=1)
119	    continue;
120	  Vector3 v1 = boost.transform(pi0[0].momentum()).p3().unit();
121	  Vector3 v2 = boost.transform(pip[0].momentum()).p3().unit();
122	  Vector3 norm = v1.cross(v2).unit();
123	  double ctheta = e1z.dot(norm);
124	  _h_ctheta_omega.fill(xE,ctheta);
125	  if(xE>0.025) _h_ctheta_omega_all->fill(ctheta);
126	}
127      }
128    }
129
130
131    /// Normalise histograms etc., after the run
132    void finalize() {
133      vector<double> x = {0.025,0.05,0.1,0.15,0.3,0.6};
134      Scatter2DPtr h_rho  ;
135      book(h_rho, 1,1,1);
136      Scatter2DPtr h_omega;
137      book(h_omega, 2,1,1);
138      for(unsigned int ix=0;ix<_h_ctheta_rho.histos().size();++ix) {
139	// rho
140	normalize(_h_ctheta_rho.histos()[ix]);
141	pair<double,double> rho00 = calcRho(_h_ctheta_rho.histos()[ix]);
142	h_rho->addPoint(0.5*(x[ix]+x[ix+1]), rho00.first, make_pair(0.5*(x[ix+1]-x[ix]),0.5*(x[ix+1]-x[ix])),
143			make_pair(rho00.second,rho00.second) );
144	// omega
145	normalize(_h_ctheta_omega.histos()[ix]);
146	rho00 = calcRho(_h_ctheta_omega.histos()[ix]);
147	h_omega->addPoint(0.5*(x[ix]+x[ix+1]), rho00.first, make_pair(0.5*(x[ix+1]-x[ix]),0.5*(x[ix+1]-x[ix])),
148			make_pair(rho00.second,rho00.second) );
149      }
150      // omega over whole range
151      Scatter2DPtr h_omega_all;
152      book(h_omega_all,2,2,1);
153      normalize(_h_ctheta_omega_all);
154      pair<double,double> rho00 = calcRho(_h_ctheta_omega_all);
155      h_omega_all->addPoint(0.5125, rho00.first, make_pair(0.4875,0.4875),
156			    make_pair(rho00.second,rho00.second) );
157    }
158
159    //@}
160
161
162    /// @name Histograms
163    //@{
164    BinnedHistogram _h_ctheta_rho,_h_ctheta_omega;
165    Histo1DPtr _h_ctheta_omega_all;
166    //@}
167
168
169  };
170
171
172  // The hook for the plugin system
173  RIVET_DECLARE_PLUGIN(OPAL_2000_I502750);
174
175
176}