rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

TPC_1991_I316132

$D^{*+}$ polarization in $e^+e^-$ at 29 GeV
Experiment: TPC (PEP)
Inspire ID: 316132
Status: UNVALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev. D43 (1991) 29-33, 1991
Beams: e+ e-
Beam energies: ANY
Run details:
  • continuum e+e- -> hadrons at 29 GeV

Measurement of the polarization of $D^{*+}$ mesons in $e^+e^-$ collisions at 29 GeV by the TPC experiment

Source code: TPC_1991_I316132.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Tools/BinnedHistogram.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief D* poliarzation at 29 GeV
 11  class TPC_1991_I316132 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(TPC_1991_I316132);
 16
 17
 18    /// @name Analysis methods
 19    //@{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(Beam(), "Beams");
 26      declare(UnstableParticles(), "UFS");
 27
 28      // Book histograms
 29      {Histo1DPtr tmp; _h_ctheta.add(0.3,0.4,book(tmp, "ctheta_0",20,-1,1));}
 30      {Histo1DPtr tmp; _h_ctheta.add(0.4,0.5,book(tmp, "ctheta_1",20,-1,1));}
 31      {Histo1DPtr tmp; _h_ctheta.add(0.5,0.6,book(tmp, "ctheta_2",20,-1,1));}
 32      {Histo1DPtr tmp; _h_ctheta.add(0.6,0.7,book(tmp, "ctheta_3",20,-1,1));}
 33      {Histo1DPtr tmp; _h_ctheta.add(0.7,0.8,book(tmp, "ctheta_4",20,-1,1));}
 34      {Histo1DPtr tmp; _h_ctheta.add(0.8,1.0,book(tmp, "ctheta_5",20,-1,1));}
 35      book(_h_ctheta_all, "ctheta_all",20,-1,1);
 36
 37      {Histo1DPtr tmp; _h_phi.add(0.3,0.4,book(tmp, "phi_0",20,-M_PI,M_PI));}
 38      {Histo1DPtr tmp; _h_phi.add(0.4,0.5,book(tmp, "phi_1",20,-M_PI,M_PI));}
 39      {Histo1DPtr tmp; _h_phi.add(0.5,0.6,book(tmp, "phi_2",20,-M_PI,M_PI));}
 40      {Histo1DPtr tmp; _h_phi.add(0.6,0.7,book(tmp, "phi_3",20,-M_PI,M_PI));}
 41      {Histo1DPtr tmp; _h_phi.add(0.7,0.8,book(tmp, "phi_4",20,-M_PI,M_PI));}
 42      {Histo1DPtr tmp; _h_phi.add(0.8,1.0,book(tmp, "phi_5",20,-M_PI,M_PI));}
 43      book(_h_phi_all, "phi_all",20,-M_PI,M_PI);
 44
 45      {Histo1DPtr tmp; _h_01.add(0.3,0.4,book(tmp, "h_01_0",20,-1,1));}
 46      {Histo1DPtr tmp; _h_01.add(0.4,0.5,book(tmp, "h_01_1",20,-1,1));}
 47      {Histo1DPtr tmp; _h_01.add(0.5,0.6,book(tmp, "h_01_2",20,-1,1));}
 48      {Histo1DPtr tmp; _h_01.add(0.6,0.7,book(tmp, "h_01_3",20,-1,1));}
 49      {Histo1DPtr tmp; _h_01.add(0.7,0.8,book(tmp, "h_01_4",20,-1,1));}
 50      {Histo1DPtr tmp; _h_01.add(0.8,1.0,book(tmp, "h_01_5",20,-1,1));}
 51      book(_h_01_all, "h_01_all",20,-1,1);
 52    }
 53
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57      // Get beams and average beam momentum
 58      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 59      const double Emax = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
 60      Vector3 axis;
 61      if(beams.first.pid()>0)
 62	axis = beams.first .momentum().p3().unit();
 63      else
 64	axis = beams.second.momentum().p3().unit();
 65
 66      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 67      for  (const Particle& p : ufs.particles(Cuts::abspid==413)) {
 68	if(p.children().size()!=2) continue;
 69	int sign = p.pid()/413;
 70	Particle D0;
 71	if(p.children()[0].pid()==sign*421 && p.children()[1].pid()==sign*211) {
 72	  D0 = p.children()[0];
 73	}
 74	else if(p.children()[1].pid()==sign*421 && p.children()[0].pid()==sign*211) {
 75	  D0 = p.children()[1];
 76	}
 77	else
 78	  continue;
 79	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 80	double xE = p.momentum().t()/Emax;
 81	Vector3 e1z = p.momentum().p3().unit();
 82	Vector3 e1y = e1z.cross(axis).unit();
 83	Vector3 e1x = e1y.cross(e1z).unit();
 84	Vector3 axis1 = boost.transform(D0.momentum()).p3().unit();
 85	double ctheta = e1z.dot(axis1);
 86	_h_ctheta.fill(xE,ctheta);
 87	double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
 88	_h_phi.fill(xE,phi);
 89	_h_01.fill(xE,ctheta,cos(phi));
 90	
 91	if(xE>0.3) {
 92	  _h_ctheta_all->fill(ctheta);
 93	  _h_phi_all->fill(phi);
 94	  _h_01_all->fill(ctheta,cos(phi));
 95	}
 96      }
 97    }
 98
 99    pair<double,double> calcRho(Histo1DPtr hist,unsigned int mode) {
100      if(hist->numEntries()==0.) return make_pair(0.,0.);
101      double sum1(0.),sum2(0.);
102      for (auto bin : hist->bins() ) {
103	double Oi = bin.area();
104	if(Oi==0.) continue;
105	double ai,bi;
106	if(mode==0) {
107	  ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
108	  bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
109	}
110	else if(mode==1) {
111	  ai = 0.5/M_PI*(bin.xMax()-bin.xMin());
112	  bi = 0.5/M_PI*(sin(2.*bin.xMin())-sin(2.*bin.xMax()));
113	}
114	else  {
115	  ai=0.;
116	  bi= sqrt(0.5)*((1.-sqr(bin.xMax()))*sqrt(1.-sqr(bin.xMax()))-
117			 (1.-sqr(bin.xMin()))*sqrt(1.-sqr(bin.xMin())));
118	}
119	double Ei = bin.areaErr();
120	sum1 += sqr(bi/Ei);
121	sum2 += bi/sqr(Ei)*(Oi-ai);
122      }
123      return make_pair(sum2/sum1,sqrt(1./sum1));
124    }
125
126    pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
127      if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
128      double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
129      double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
130      double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
131      for (auto bin : hist->bins() ) {
132       	double Oi = bin.area();
133	if(Oi==0.) continue;
134	double a =  d*(bin.xMax() - bin.xMin());
135	double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
136       	double Ei = bin.areaErr();
137	sum1 +=   a*Oi/sqr(Ei);
138	sum2 +=   b*Oi/sqr(Ei);
139	sum3 += sqr(a)/sqr(Ei);
140	sum4 += sqr(b)/sqr(Ei);
141	sum5 +=    a*b/sqr(Ei);
142      }
143      // calculate alpha
144      double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
145      // and error
146      double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
147      double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
148      double aa =  sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);      
149      double dis = sqr(bb)-4.*aa*cc;
150      if(dis>0.) {
151	dis = sqrt(dis);
152	return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
153      }
154      else {
155	return make_pair(alpha,make_pair(0.,0.));
156      }
157    }
158    
159    /// Normalise histograms etc., after the run
160    void finalize() {
161      vector<double> x = {0.3,0.4,0.5,0.6,0.7,0.8,1.};
162      Scatter2DPtr h_alpha, h_rho00, h_rhooff, h_01;
163      book(h_alpha , 1,1,1);
164      book(h_rho00 , 2,1,1);
165      book(h_rhooff, 2,1,2);
166      book(h_01    , 2,1,3);
167      for(unsigned int ix=0;ix<6;++ix) {
168	double integral = _h_ctheta.histos()[ix]->integral();
169	scale( _h_ctheta.histos()[ix],1./integral);
170	scale( _h_01.histos()[ix],1./integral);
171	normalize(_h_phi.histos()[ix]);
172	pair<double,double> rho00 = calcRho(_h_ctheta.histos()[ix],0);
173	h_rho00->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])),
174			  make_pair(rho00.second,rho00.second) );
175	pair<double,pair<double,double> > alpha = calcAlpha(_h_ctheta.histos()[ix]);
176	h_alpha->addPoint(0.5*(x[ix]+x[ix+1]), alpha.first, make_pair(0.5*(x[ix+1]-x[ix]),0.5*(x[ix+1]-x[ix])),
177			  alpha.second);
178	pair<double,double> rhooff = calcRho(_h_phi.histos()[ix],1);
179	h_rhooff->addPoint(0.5*(x[ix]+x[ix+1]), rhooff.first, make_pair(0.5*(x[ix+1]-x[ix]),0.5*(x[ix+1]-x[ix])),
180			  make_pair(rhooff.second,rhooff.second) );
181	pair<double,double> rho01 = calcRho(_h_01.histos()[ix],2);
182	h_01->addPoint(0.5*(x[ix]+x[ix+1]), rho01.first, make_pair(0.5*(x[ix+1]-x[ix]),0.5*(x[ix+1]-x[ix])),
183		       make_pair(rho01.second,rho01.second) );
184      }
185      // integral over z
186      book(h_alpha , 1,2,1);
187      book(h_rho00 , 2,2,1);
188      book(h_rhooff, 2,2,2);
189      book(h_01    , 2,2,3);
190      double integral = _h_ctheta_all->integral();
191      scale( _h_ctheta_all,1./integral);
192      scale( _h_01_all,1./integral);
193      normalize(_h_phi_all);
194      pair<double,double> rho00 = calcRho(_h_ctheta_all,0);
195      h_rho00->addPoint(.65, rho00.first, make_pair(0.35,0.35),
196			make_pair(rho00.second,rho00.second) );
197      pair<double,pair<double,double> > alpha = calcAlpha(_h_ctheta_all);
198      h_alpha->addPoint(.65, alpha.first, make_pair(0.35,0.35),
199			alpha.second);
200      pair<double,double> rhooff = calcRho(_h_phi_all,1);
201      h_rhooff->addPoint(.65, rhooff.first, make_pair(0.35,0.35),
202			 make_pair(rhooff.second,rhooff.second) );
203      pair<double,double> rho01 = calcRho(_h_01_all,2);
204      h_01->addPoint(.65, rho01.first, make_pair(0.35,0.35),
205		     make_pair(rho01.second,rho01.second) );
206    }
207
208    //@}
209
210
211    /// @name Histograms
212    //@{
213    BinnedHistogram _h_ctheta,_h_phi,_h_01;
214    Histo1DPtr _h_ctheta_all,_h_phi_all,_h_01_all;
215    //@}
216
217
218  };
219
220
221  // The hook for the plugin system
222  RIVET_DECLARE_PLUGIN(TPC_1991_I316132);
223
224
225}