Rivet analyses referenceCLEO_1998_I467595$D^{*+}$ polarization in $e^+e^-$ at 10.5 GeVExperiment: CLEO (CESR) Inspire ID: 467595 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Measurement of the polarization of $D^{*+}$ mesons in $e^+e^-$ collisions at 10.5 GeV by CLEO Source code: CLEO_1998_I467595.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* helicity in e+e- at 10.5 GeV
11 class CLEO_1998_I467595 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1998_I467595);
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 temp; _h_ctheta.add(0.25,0.45,book(temp,5,1,1));}
30 {Histo1DPtr temp; _h_ctheta.add(0.45,0.55,book(temp,5,1,2));}
31 {Histo1DPtr temp; _h_ctheta.add(0.55,0.65,book(temp,5,1,3));}
32 {Histo1DPtr temp; _h_ctheta.add(0.65,0.75,book(temp,5,1,4));}
33 {Histo1DPtr temp; _h_ctheta.add(0.75,0.85,book(temp,5,1,5));}
34 {Histo1DPtr temp; _h_ctheta.add(0.85,1. ,book(temp,5,1,6));}
35
36 }
37
38
39 /// Perform the per-event analysis
40 void analyze(const Event& event) {
41 // Get beams and average beam momentum
42 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
43 const double Emax = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
44 const double Pmax = sqrt(sqr(Emax)-sqr(2.01026));
45
46 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
47 for (const Particle& p : ufs.particles(Cuts::abspid==413)) {
48 if(p.children().size()!=2) continue;
49 int sign = p.pid()/413;
50 Particle D0;
51 if(p.children()[0].pid()==sign*421 && p.children()[1].pid()==sign*211) {
52 D0 = p.children()[0];
53 }
54 else if(p.children()[1].pid()==sign*421 && p.children()[0].pid()==sign*211) {
55 D0 = p.children()[1];
56 }
57 else
58 continue;
59 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
60 double xp = (p.momentum().p3().mod()+p.momentum().t())/(Pmax+Emax);
61 Vector3 e1z = p.momentum().p3().unit();
62 Vector3 axis1 = boost.transform(D0.momentum()).p3().unit();
63 double ctheta = e1z.dot(axis1);
64 _h_ctheta.fill(xp,ctheta);
65 }
66 }
67
68 pair<double,double> calcRho(Histo1DPtr hist) {
69 if(hist->numEntries()==0.) return make_pair(0.,0.);
70 double sum1(0.),sum2(0.);
71 for (auto bin : hist->bins() ) {
72 double Oi = bin.area();
73 if(Oi==0.) continue;
74 double ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
75 double bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
76 double Ei = bin.areaErr();
77 sum1 += sqr(bi/Ei);
78 sum2 += bi/sqr(Ei)*(Oi-ai);
79 }
80 return make_pair(sum2/sum1,sqrt(1./sum1));
81 }
82
83 pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
84 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
85 double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
86 double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
87 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
88 for (auto bin : hist->bins() ) {
89 double Oi = bin.area();
90 if(Oi==0.) continue;
91 double a = d*(bin.xMax() - bin.xMin());
92 double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
93 double Ei = bin.areaErr();
94 sum1 += a*Oi/sqr(Ei);
95 sum2 += b*Oi/sqr(Ei);
96 sum3 += sqr(a)/sqr(Ei);
97 sum4 += sqr(b)/sqr(Ei);
98 sum5 += a*b/sqr(Ei);
99 }
100 // calculate alpha
101 double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
102 // and error
103 double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
104 double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
105 double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
106 double dis = sqr(bb)-4.*aa*cc;
107 if(dis>0.) {
108 dis = sqrt(dis);
109 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
110 }
111 else {
112 return make_pair(alpha,make_pair(0.,0.));
113 }
114 }
115
116 /// Normalise histograms etc., after the run
117 void finalize() {
118 vector<double> x = {0.25,0.45,0.55,0.65,0.75,0.85,1.};
119 Scatter2DPtr h_alpha;
120 book(h_alpha,3,1,1);
121 Scatter2DPtr h_rho;
122 book(h_rho ,4,1,1);
123 for(unsigned int ix=0;ix<6;++ix) {
124 normalize(_h_ctheta.histos()[ix]);
125 pair<double,double> rho00 = calcRho(_h_ctheta.histos()[ix]);
126 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])),
127 make_pair(rho00.second,rho00.second) );
128 pair<double,pair<double,double> > alpha = calcAlpha(_h_ctheta.histos()[ix]);
129 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])),
130 alpha.second);
131 }
132 }
133
134 //@}
135
136
137 /// @name Histograms
138 //@{
139 BinnedHistogram _h_ctheta;
140 //@}
141
142
143 };
144
145
146 // The hook for the plugin system
147 RIVET_DECLARE_PLUGIN(CLEO_1998_I467595);
148
149
150}
|