rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

OPAL_1997_I447188

Polarization of $\Lambda^0$ baryons at LEP 1
Experiment: OPAL (LEP)
Inspire ID: 447188
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Eur.Phys.J. C2 (1998) 49-59
Beams: e- e+
Beam energies: ANY
Run details:
  • e+e- -> hadrons at 91.2 GeV

Measurement of the polarization of $\Lambda^0$ baryons at LEP 1. The $\cos\theta$ and $\cos\phi$ distributions are extracted and then the polarization obtained by fitting the resulting distributions.

Source code: OPAL_1997_I447188.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6#include "Rivet/Tools/BinnedHistogram.hh"
  7#include "Rivet/Projections/ChargedFinalState.hh"
  8#include <sstream>
  9namespace Rivet {
 10
 11
 12  /// @brief Lambda polarization at LEP1
 13  class OPAL_1997_I447188 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1997_I447188);
 18
 19
 20    /// @name Analysis methods
 21    //@{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Initialise and register projections
 27      declare(Beam(), "Beams");
 28      const ChargedFinalState cfs;
 29      const Thrust thrust(cfs);
 30      declare(thrust, "Thrust");
 31      declare(UnstableParticles(), "UFS");
 32
 33      // Book the histograms
 34      {Histo1DPtr temp; _h_ctheta.add(0.027, 0.05 ,book(temp, "/TMP/ctheta_0",20,-1.,1.));}
 35      {Histo1DPtr temp; _h_ctheta.add(0.05 , 0.08 ,book(temp, 4,1,1));}
 36      {Histo1DPtr temp; _h_ctheta.add(0.08 , 0.09 ,book(temp, "/TMP/ctheta_2",20,-1.,1.));}
 37      {Histo1DPtr temp; _h_ctheta.add(0.09 , 0.1  ,book(temp, "/TMP/ctheta_3",20,-1.,1.));}
 38      {Histo1DPtr temp; _h_ctheta.add(0.1  , 0.15 ,book(temp, 4,1,2));}
 39      {Histo1DPtr temp; _h_ctheta.add(0.15 , 0.2  ,book(temp, "/TMP/ctheta_5",20,-1.,1.));}
 40      {Histo1DPtr temp; _h_ctheta.add(0.2  , 0.3  ,book(temp, 4,1,3));}
 41      {Histo1DPtr temp; _h_ctheta.add(0.3  , 0.4  ,book(temp, "/TMP/ctheta_7",20,-1.,1.));}
 42      {Histo1DPtr temp; _h_ctheta.add(0.4  , 1.0  ,book(temp, "/TMP/ctheta_8",20,-1.,1.));}
 43      book(_h_ctheta_large, 4,1,4);
 44      
 45      {Histo1DPtr temp; _h_cphi .add(0.3,0.6,book(temp, "/TMP/cphiP_0",10,0.,1.));}
 46      {Histo1DPtr temp; _h_cphi .add(0.3,0.6,book(temp, 5,1,1));}
 47      {Histo1DPtr temp; _h_cphi .add(0.6,0.9,book(temp, 5,1,2));}
 48      {Histo1DPtr temp; _h_cphi .add(0.9,1.2,book(temp, "/TMP/cphiP_3",10,0.,1.));}
 49      {Histo1DPtr temp; _h_cphi .add(1.2,1.5,book(temp, "/TMP/cphiP_4",10,0.,1.));}
 50      
 51      book(_h_cphi_low , 5,1,4);
 52      book(_h_cphi_mid , "/TMP/cphiP_mid" ,10,0.,1.);
 53      book(_h_cphi_high, 5,1,3);
 54
 55      {Histo1DPtr temp; _h_plus_lam.add(0.027,0.05,book(temp, "/TMP/lamP_0",20,-1.,1.));}
 56      {Histo1DPtr temp; _h_plus_lam.add(0.05 ,0.08,book(temp, "/TMP/lamP_1",20,-1.,1.));}
 57      {Histo1DPtr temp; _h_plus_lam.add(0.08 ,0.09,book(temp, "/TMP/lamP_2",20,-1.,1.));}
 58      {Histo1DPtr temp; _h_plus_lam.add(0.09 ,0.1 ,book(temp, "/TMP/lamP_3",20,-1.,1.));}
 59      {Histo1DPtr temp; _h_plus_lam.add(0.1  ,0.15,book(temp, "/TMP/lamP_4",20,-1.,1.));}
 60      {Histo1DPtr temp; _h_plus_lam.add(0.15 ,0.2 ,book(temp, "/TMP/lamP_5",20,-1.,1.));}
 61      {Histo1DPtr temp; _h_plus_lam.add(0.2  ,0.3 ,book(temp, "/TMP/lamP_6",20,-1.,1.));}
 62      {Histo1DPtr temp; _h_plus_lam.add(0.3  ,0.4 ,book(temp, "/TMP/lamP_7",20,-1.,1.));}
 63      {Histo1DPtr temp; _h_plus_lam.add(0.4  ,1.0 ,book(temp, "/TMP/lamP_8",20,-1.,1.));}
 64      
 65      {Histo1DPtr temp; _h_minus_lam.add(0.027,0.05,book(temp, "/TMP/lamM_0",20,-1.,1.));}
 66      {Histo1DPtr temp; _h_minus_lam.add(0.05 ,0.08,book(temp, "/TMP/lamM_1",20,-1.,1.));}
 67      {Histo1DPtr temp; _h_minus_lam.add(0.08 ,0.09,book(temp, "/TMP/lamM_2",20,-1.,1.));}
 68      {Histo1DPtr temp; _h_minus_lam.add(0.09 ,0.1 ,book(temp, "/TMP/lamM_3",20,-1.,1.));}
 69      {Histo1DPtr temp; _h_minus_lam.add(0.1  ,0.15,book(temp, "/TMP/lamM_4",20,-1.,1.));}
 70      {Histo1DPtr temp; _h_minus_lam.add(0.15 ,0.2 ,book(temp, "/TMP/lamM_5",20,-1.,1.));}
 71      {Histo1DPtr temp; _h_minus_lam.add(0.2  ,0.3 ,book(temp, "/TMP/lamM_6",20,-1.,1.));}
 72      {Histo1DPtr temp; _h_minus_lam.add(0.3  ,0.4 ,book(temp, "/TMP/lamM_7",20,-1.,1.));}
 73      {Histo1DPtr temp; _h_minus_lam.add(0.4  ,1.0 ,book(temp, "/TMP/lamM_8",20,-1.,1.));}
 74
 75      {Histo1DPtr temp; _h_plus_lam_large1  = book(temp, "/TMP/lamP_large_1",20,-1.,1.);}
 76      {Histo1DPtr temp; _h_plus_lam_large2  = book(temp, "/TMP/lamP_large_2",20,-1.,1.);}
 77      {Histo1DPtr temp; _h_minus_lam_large1 = book(temp, "/TMP/lamM_large_1",20,-1.,1.);}
 78      {Histo1DPtr temp; _h_minus_lam_large2 = book(temp, "/TMP/lamM_large_2",20,-1.,1.);}
 79    }
 80
 81
 82    /// Perform the per-event analysis
 83    void analyze(const Event& event) {
 84
 85      // Get beams and average beam momentum
 86      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 87      const double meanBeamMom = ( beams.first.p3().mod() +
 88                                   beams.second.p3().mod() ) / 2.0;
 89      Vector3 beamAxis;
 90      if(beams.first.pid()==-11) {
 91	beamAxis = beams.first .momentum().p3().unit();
 92      }
 93      else {
 94	beamAxis = beams.second.momentum().p3().unit();
 95      }
 96	
 97      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 98      // thrust, to define an axis
 99      const Thrust& thrust = apply<Thrust>(event, "Thrust");
100      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
101
102      for(const Particle & lambda : ufs.particles(Cuts::abspid==3122)) {
103	double xE = lambda.momentum().t()/meanBeamMom;
104	int sign = lambda.pid()/3122;
105	Vector3 axis1 = lambda.momentum().p3().unit();
106	// assymetry
107	double cLam = axis1.dot(beamAxis);
108	if(sign>0)
109	  _h_plus_lam .fill(xE,cLam);
110	else
111	  _h_minus_lam.fill(xE,cLam);
112	if(xE>0.15) {
113	  if(sign>0)
114	    _h_plus_lam_large1 ->fill(cLam);
115	  else
116	    _h_minus_lam_large1->fill(cLam);
117	}
118	if(xE>0.3) {
119	  if(sign>0)
120	    _h_plus_lam_large2 ->fill(cLam);
121	  else
122	    _h_minus_lam_large2->fill(cLam);
123	}
124	if(lambda.children().size()!=2) continue;
125	// look at the decay products
126	Particle proton,pion;
127	if(lambda.children()[0].pid()==sign*2212 && 
128	   lambda.children()[1].pid()==-sign*211) {
129	  proton = lambda.children()[0];
130	  pion   = lambda.children()[1];
131	}
132	else if(lambda.children()[1].pid()==sign*2212 && 
133		lambda.children()[0].pid()==-sign*211) {
134	  proton = lambda.children()[1];
135	  pion   = lambda.children()[0];
136	}
137	else
138	  continue;
139	// boost to the lambda rest frame
140	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(lambda.momentum().betaVec());
141	FourMomentum pproton = boost.transform(proton.momentum());
142	// longitudinal polarization
143	double ctheta = axis1.dot(pproton.p3().unit());
144	_h_ctheta.fill(xE,ctheta);
145	if(xE>=0.3)  _h_ctheta_large->fill(ctheta);
146	// transverse polarization
147	Vector3 axis2;
148	if(lambda.momentum().p3().dot(thrust.thrustAxis())>=0.) {
149	  axis2 = thrust.thrustAxis();
150	}
151	else {
152	  axis2 =-thrust.thrustAxis();
153	}
154	Vector3 axis3 = axis2.cross(axis1).unit();	
155	double pT = sqrt(sqr(thrust.thrustMajorAxis().dot(lambda.momentum().p3()))+
156			 sqr(thrust.thrustMinorAxis().dot(lambda.momentum().p3())));
157	double cPhi = axis3.dot(pproton.p3().unit());
158	_h_cphi.fill(pT,cPhi);
159	if(pT>0.3) _h_cphi_low ->fill(cPhi);
160	if(pT>0.6) _h_cphi_mid ->fill(cPhi);
161	if(pT>1.5) _h_cphi_high->fill(cPhi);
162      }
163
164    }
165
166    pair<double,double> calcAlpha(Histo1DPtr hist) {
167      if(hist->numEntries()==0.) return make_pair(0.,0.);
168      double sum1(0.),sum2(0.);
169      for (auto bin : hist->bins() ) {
170	double Oi = bin.area();
171	if(Oi==0.) continue;
172	double ai = 0.5*(bin.xMax()-bin.xMin());
173	double bi = 0.5*ai*(bin.xMax()+bin.xMin());
174	double Ei = bin.areaErr();
175	sum1 += sqr(bi/Ei);
176	sum2 += bi/sqr(Ei)*(Oi-ai);
177      }
178      return make_pair(sum2/sum1,sqrt(1./sum1));
179    }
180    
181    pair<double,double> calcAsymmetry(Scatter2DPtr hist,unsigned int mode) {
182      double sum1(0.),sum2(0.);
183      for (auto bin : hist->points() ) {
184	double Oi = bin.y();
185	if(Oi==0.) continue;
186	double bi;
187	if(mode==0)
188	  bi = 0.25*(bin.xMax()-bin.xMin())*(bin.xMax()+bin.xMin());
189	else
190	  bi = 4.*(bin.xMax()+bin.xMin())/(3.+sqr(bin.xMax())+bin.xMax()*bin.xMin()+sqr(bin.xMin()));
191	double Ei = bin.yErrAvg();
192	sum1 += sqr(bi/Ei);
193	sum2 += bi/sqr(Ei)*Oi;
194      }
195      return make_pair(sum2/sum1,sqrt(1./sum1));
196    }
197
198
199    /// Normalise histograms etc., after the run
200    void finalize() {
201      // longitudinal polarization
202      vector<double> x_val = {3.850000e-02,6.500000e-02,8.500000e-02,9.500000e-02,1.250000e-01,
203			      1.750000e-01,2.500000e-01,3.500000e-01,7.000000e-01};
204      vector<double> x_err = {1.150000e-02,1.500000e-02,5.000000e-03,5.000000e-03,2.500000e-02,
205			      2.500000e-02,5.000000e-02,5.000000e-02,3.000000e-01};
206      unsigned int ipoint=0;
207      double aLam = 0.642;
208      Scatter2DPtr h_long;
209      book(h_long, 1,1,1);
210      for(Histo1DPtr & hist : _h_ctheta.histos()) {
211	normalize(hist);
212	pair<double,double> alpha = calcAlpha(hist);
213	alpha.first  /=aLam;
214	alpha.second /=aLam;
215	h_long->addPoint(x_val[ipoint], 100.*alpha.first, make_pair(x_err[ipoint],x_err[ipoint]),
216			   make_pair(100.*alpha.second,100.*alpha.second) );	
217	++ipoint;
218      }
219      normalize(_h_ctheta_large);
220      pair<double,double> alpha = calcAlpha(_h_ctheta_large);
221      alpha.first  /=aLam;
222      alpha.second /=aLam;
223      Scatter2DPtr h_long_l;
224      book(h_long_l,1,1,2);
225      h_long_l->addPoint(0.65, 100.*alpha.first, make_pair(0.35,0.35), make_pair(100.*alpha.second,100.*alpha.second) );
226      // transverse polarization
227      vector<double> pT_val = {0.15,0.45,0.75,1.05,1.35};
228      vector<double> pT_err = {0.15,0.15,0.15,0.15,0.15};
229      Scatter2DPtr h_trans;
230      book(h_trans,2,1,1);
231      ipoint=0;
232      for(Histo1DPtr & hist : _h_cphi.histos()) {
233	normalize(hist);
234	pair<double,double> alpha = calcAlpha(hist);
235	alpha.first  /=aLam;
236	alpha.second /=aLam;
237	h_trans->addPoint(pT_val[ipoint], 100.*alpha.first, make_pair(pT_err[ipoint],pT_err[ipoint]),
238			   make_pair(100.*alpha.second,100.*alpha.second) );	
239	++ipoint;
240      }
241      normalize(_h_cphi_low);
242      alpha = calcAlpha(_h_cphi_low);
243      alpha.first  /=aLam;
244      alpha.second /=aLam;
245      Scatter2DPtr h_trans_low;
246      book(h_trans_low,2,1,2);
247      h_trans_low->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(100.*alpha.second,100.*alpha.second) );
248      normalize(_h_cphi_mid);
249      alpha = calcAlpha(_h_cphi_mid);
250      alpha.first  /=aLam;
251      alpha.second /=aLam;
252      Scatter2DPtr h_trans_mid;
253      book(h_trans_mid,2,1,3);
254      h_trans_mid->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(100.*alpha.second,100.*alpha.second) );
255      normalize(_h_cphi_high);
256      alpha = calcAlpha(_h_cphi_high);
257      alpha.first  /=aLam;
258      alpha.second /=aLam;
259      Scatter2DPtr h_trans_high;
260      book(h_trans_high,2,1,4);
261      h_trans_high->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(100.*alpha.second,100.*alpha.second) );
262      // asyymetry
263      Scatter2DPtr h_asym;
264      book(h_asym,3,1,1);
265      for(unsigned int ix=0;ix<9;++ix) {
266       	normalize(_h_plus_lam.histos()[ix] );
267       	normalize(_h_minus_lam.histos()[ix]);
268	std::ostringstream title;
269       	title << "/TMP/a_lam_" << ix;
270       	Scatter2DPtr sTemp;
271	book(sTemp,title.str());
272       	asymm(_h_plus_lam.histos()[ix],_h_minus_lam.histos()[ix],sTemp);
273       	pair<double,double> alpha = calcAsymmetry(sTemp,1);
274       	h_asym->addPoint(x_val[ix],-alpha.first, make_pair(x_err[ix],x_err[ix]),
275			 make_pair(alpha.second,alpha.second) );
276      }
277      normalize(_h_plus_lam_large1 );
278      normalize(_h_minus_lam_large1);
279      Scatter2DPtr sTemp;
280      book(sTemp,"/TMP/a_lam_large1");
281      asymm(_h_plus_lam_large1,_h_minus_lam_large1,sTemp);
282      alpha = calcAsymmetry(sTemp,1);
283      book(h_asym,3,1,2);
284      h_asym->addPoint(0.575,-alpha.first, make_pair(0.425,0.425),
285		       make_pair(alpha.second,alpha.second) );
286      normalize(_h_plus_lam_large2 );
287      normalize(_h_minus_lam_large2);
288      book(sTemp,"/TMP/a_lam_large2");
289      asymm(_h_plus_lam_large2,_h_minus_lam_large2,sTemp);
290      alpha = calcAsymmetry(sTemp,1);
291      book(h_asym,3,1,3);
292      h_asym->addPoint(0.65,-alpha.first, make_pair(0.35,0.35),
293		       make_pair(alpha.second,alpha.second) );
294    }
295
296    //@}
297
298
299    /// @name Histograms
300    //@{
301    BinnedHistogram _h_ctheta,_h_cphi,_h_plus_lam,_h_minus_lam;
302    Histo1DPtr _h_ctheta_large;
303    Histo1DPtr _h_cphi_low, _h_cphi_mid, _h_cphi_high;
304    Histo1DPtr _h_plus_lam_large1,_h_plus_lam_large2;
305    Histo1DPtr _h_minus_lam_large1,_h_minus_lam_large2;
306    //@}
307
308  };
309
310
311  // The hook for the plugin system
312  RIVET_DECLARE_PLUGIN(OPAL_1997_I447188);
313
314
315}