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