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