Rivet analyses referenceBELLE_2015_I1369998$\bar{B}^0\to D^{*+}\omega\pi^-$ decaysExperiment: BELLE (KEKB) Inspire ID: 1369998 Status: VALIDATED NOHEPDATA Authors:
Beam energies: ANY Run details:
Mass and aangular distributions in $\bar{B}^0\to D^{*+}\omega\pi^-$ decays. Data read from plots with the backgrounds given subtracted. Source code: BELLE_2015_I1369998.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include "Rivet/Projections/DecayedParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief B -> D* omega pi
10 class BELLE_2015_I1369998 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2015_I1369998);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22 UnstableParticles ufs = UnstableParticles(Cuts::abspid==511);
23 declare(ufs, "UFS");
24 DecayedParticles B0(ufs);
25 B0.addStable( 413);
26 B0.addStable(-413);
27 B0.addStable( 223);
28 declare(B0, "B0");
29 for(unsigned int ix=0;ix<4;++ix)
30 for(unsigned int iy=0;iy<6;++iy)
31 book(_h[ix][iy],1+ix,1,1+iy);
32 }
33
34 void findChildren(const Particle & p, Particles & pim, Particles & pip,
35 Particles & pi0, unsigned int &ncount) {
36 for( const Particle &child : p.children()) {
37 if(child.pid()==PID::PIPLUS) {
38 pip.push_back(child);
39 ncount+=1;
40 }
41 else if(child.pid()==PID::PIMINUS) {
42 pim.push_back(child);
43 ncount+=1;
44 }
45 else if(child.pid()==PID::PI0) {
46 pi0.push_back(child);
47 ncount+=1;
48 }
49 else if(child.children().empty()) {
50 ncount+=1;
51 }
52 else
53 findChildren(child,pim,pip,pi0,ncount);
54 }
55 }
56
57 /// Perform the per-event analysis
58 void analyze(const Event& event) {
59 static const map<PdgId,unsigned int> & mode = { { 413,1},{ 223,1}, {-211,1}};
60 static const map<PdgId,unsigned int> & modeCC = { {-413,1},{ 223,1}, { 211,1}};
61 DecayedParticles B0 = apply<DecayedParticles>(event, "B0");
62 // loop over particles
63 for(unsigned int ix=0;ix<B0.decaying().size();++ix) {
64 int sign = 1;
65 if(B0.decaying()[ix].pid()<0 && B0.modeMatches(ix,3,mode))
66 sign = 1;
67 else if(B0.decaying()[ix].pid()>0 && B0.modeMatches(ix,3,modeCC))
68 sign = -1;
69 else
70 continue;
71 const Particle & Dstar = B0.decayProducts()[ix].at( sign*413)[0];
72 const Particle & omega = B0.decayProducts()[ix].at( 223)[0];
73 const Particle & pim1 = B0.decayProducts()[ix].at(-sign*211)[0];
74 // mass hists, no cuts
75 double mOmegaPi2 = (omega.momentum()+pim1.momentum()).mass2();
76 _h[0][0]->fill(mOmegaPi2);
77 double mDstarpi2 = (Dstar.momentum()+pim1.momentum()).mass2();
78 _h[1][0]->fill(mDstarpi2);
79 // check the no of decay products
80 if(Dstar.children().size()!=2 || omega.children().size()!=3)
81 continue;
82 // find the children of the D* meson
83 Particle D0,pip1;
84 if(Dstar.children()[0].pid()==sign*211 &&
85 Dstar.children()[1].pid()==sign*421) {
86 pip1 = Dstar.children()[0];
87 D0 = Dstar.children()[1];
88 }
89 else if(Dstar.children()[1].pid()==sign*211 &&
90 Dstar.children()[0].pid()==sign*421) {
91 pip1 = Dstar.children()[1];
92 D0 = Dstar.children()[0];
93 }
94 else
95 continue;
96 // children of the omega
97 unsigned int ncount=0;
98 Particles pip,pim,pi0;
99 findChildren(omega,pim,pip,pi0,ncount);
100 if( ncount!=3 || !(pim.size()==1 && pip.size()==1 && pi0.size()==1)) continue;
101 // first bottom to the B frame
102 LorentzTransform boostB = LorentzTransform::mkFrameTransformFromBeta(B0.decaying()[ix].momentum().betaVec());
103 FourMomentum pOmega = boostB.transform(omega.momentum());
104 FourMomentum pDstar = boostB.transform(Dstar.momentum());
105 FourMomentum pD = boostB.transform(D0 .momentum());
106 FourMomentum ppim1 = boostB.transform(pim1 .momentum());
107 FourMomentum ppim2 = boostB.transform(pim[0].momentum());
108 FourMomentum ppip1 = boostB.transform(pip1 .momentum());
109 FourMomentum ppip2 = boostB.transform(pip[0].momentum());
110 // ---------------------- First set of angles --------------------------------------
111 // first the angles for D* (pi omega)
112 LorentzTransform boostD = LorentzTransform::mkFrameTransformFromBeta(pDstar.betaVec());
113 Vector3 axisD = boostD.transform(pD ).p3().unit();
114 Vector3 axispip1 = boostD.transform(ppip1).p3().unit();
115 Vector3 axisDstar = (pOmega+ppim1).p3().unit();
116 double cBeta1 = axisDstar.dot(axisD);
117 _h[0][3]->fill(cBeta1);
118 LorentzTransform boostWpi = LorentzTransform::mkFrameTransformFromBeta((pOmega+ppim1).betaVec());
119 FourMomentum pOmega2 = boostWpi.transform(pOmega);
120 Vector3 axisW = pOmega2.p3().unit();
121 Vector3 axisWpi = (pOmega+ppim1).p3().unit();
122 double cXi1 = axisWpi.dot(axisW);
123 _h[0][1]->fill(cXi1);
124 // now angle between the two planes
125 Vector3 transW = axisW-cXi1*axisWpi;
126 Vector3 transD = axisD-cBeta1*axisDstar;
127 double psi1 = atan2(transW.cross(transD).dot(axisDstar), transW.dot(transD));
128 _h[0][5]->fill(psi1);
129 // normal to omega decay plane
130 LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(pOmega2.betaVec());
131 FourMomentum ppim3 = boostW.transform(boostWpi.transform(ppim2));
132 FourMomentum ppip3 = boostW.transform(boostWpi.transform(ppip2));
133 Vector3 nW = ppim3.p3().cross(ppip3.p3()).unit();
134 // boost B decay products to omega rest frame
135 FourMomentum pOmegaPi = boostW.transform(boostWpi.transform(pOmega+ppim1));
136 FourMomentum pDstar2 = boostW.transform(boostWpi.transform(pDstar));
137 Vector3 axisWpi2 = pOmegaPi.p3().unit();
138 double cTheta1 = axisWpi2.dot(nW);
139 transW = nW-cTheta1*axisWpi2;
140 transD = pDstar2.p3().unit()-pDstar2.p3().unit().dot(axisWpi2)*axisWpi2;
141 double phi1 = atan2(transW.cross(transD).dot(axisWpi2), transW.dot(transD));
142 _h[0][2]->fill(cTheta1);
143 _h[0][4]->fill(phi1);
144 // ---------------------- Second set of angles --------------------------------------
145 // boost to D* pi frame
146 LorentzTransform boostDpi = LorentzTransform::mkFrameTransformFromBeta((pDstar+ppim1).betaVec());
147 pDstar2 = boostDpi.transform(pDstar);
148 pOmega2 = boostDpi.transform(pOmega);
149 axisW = pOmega2.p3().unit();
150 axisDstar = pDstar2.p3().unit();
151 double cXi2 = axisW.dot(axisDstar);
152 _h[1][1]->fill(cXi2);
153 // boost to D* rest frame
154 LorentzTransform boostDstar = LorentzTransform::mkFrameTransformFromBeta(pDstar2.betaVec());
155 axisW = boostDstar.transform(pOmega2).p3().unit();
156 Vector3 axisDSpi= boostDstar.transform(boostDpi.transform(pDstar+ppim1)).p3().unit();
157 axisD= boostDstar.transform(boostDpi.transform(pD)).p3().unit();
158 double cBeta2 = axisD.dot(axisDSpi);
159 _h[1][3]->fill(cBeta2);
160 transW = axisW-axisW.dot(axisDSpi)*axisDSpi;
161 transD = axisD-cBeta2*axisDSpi;
162 double psi2 = atan2(transW.cross(transD).dot(axisDSpi), transW.dot(transD));
163 _h[1][5]->fill(psi2);
164 // boost to omega frame
165 boostW = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
166 ppim3 = boostW.transform(ppim2);
167 ppip3 = boostW.transform(ppip2);
168 nW = ppim3.p3().cross(ppip3.p3()).unit();
169 axisDSpi = boostW.transform(pDstar+ppim1).p3().unit();
170 axisDstar = boostW.transform(pDstar).p3().unit();
171 double cTheta2 = axisDSpi.dot(nW);
172 _h[1][2]->fill(cTheta2);
173 transW = nW-cTheta2*axisDSpi;
174 transD = axisDstar.unit()-axisDstar.dot(axisDSpi)*axisDSpi;
175 double phi2 = atan2(transW.cross(transD).dot(axisDSpi), transW.dot(transD));
176 _h[1][4]->fill(psi2);
177 // restricted plots
178 if(abs(cTheta1)>.5) {
179 _h[2][0]->fill(mOmegaPi2);
180 }
181 else {
182 _h[2][1]->fill(mOmegaPi2);
183 _h[2][3]->fill(cBeta1);
184 _h[2][5]->fill(psi1);
185 _h[3][1]->fill(mDstarpi2);
186 _h[3][3]->fill(cTheta2);
187 _h[3][5]->fill(phi2);
188 }
189 if(cXi2>-.4) {
190 _h[2][2]->fill(cBeta1);
191 _h[2][4]->fill(psi1);
192 _h[3][0]->fill(mDstarpi2);
193 _h[3][2]->fill(cTheta2);
194 _h[3][4]->fill(phi2);
195 }
196 }
197 }
198
199
200 /// Normalise histograms etc., after the run
201 void finalize() {
202 for(unsigned int ix=0;ix<4;++ix)
203 for(unsigned int iy=0;iy<6;++iy)
204 normalize(_h[ix][iy],1.,false);
205 }
206
207 /// @}
208
209
210 /// @name Histograms
211 /// @{
212 Histo1DPtr _h[4][6];
213 /// @}
214
215
216 };
217
218
219 RIVET_DECLARE_PLUGIN(BELLE_2015_I1369998);
220
221}
|