Rivet analyses referenceCLEO_2001_I554175Mass and angular distributions in $B\to D^{(*)} \pi^+\pi^-\pi^-\pi^0$ decaysExperiment: CLEO (CESR) Inspire ID: 554175 Status: VALIDATED NOHEPDATA Authors:
Beam energies: ANY Run details:
Measurement of mass and angular distributions in $B\to D^{(*)} \pi^+\pi^-\pi^-\pi^0$ decays, primarily $B\to D^{(*)} \omega\pi^-$. The background subtracted data were read from the figures in the paper. Source code: CLEO_2001_I554175.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(*) pi+pi-pi-pi0
10 class CLEO_2001_I554175 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2001_I554175);
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 || Cuts::abspid==521);
23 declare(ufs, "UFS");
24 DecayedParticles BB(ufs);
25 BB.addStable(PID::PI0);
26 BB.addStable( 413);
27 BB.addStable(-413);
28 BB.addStable( 423);
29 BB.addStable(-423);
30 BB.addStable( 411);
31 BB.addStable(-411);
32 BB.addStable( 421);
33 BB.addStable(-421);
34 declare(BB, "BB");
35 // // histos
36 for (unsigned int ix=0; ix<6; ++ix) {
37 book(_h[ix],1+ix,1,1);
38 }
39 for (unsigned int ix=0; ix<3; ++ix) {
40 book(_h_angle[ix],7,1,1+ix);
41 }
42 book(_h_sum,8,1,1);
43 for (unsigned int ix=0;ix<2;++ix) {
44 book(_c[ix],"TMP/nB_"+toString(ix+1));
45 }
46 }
47
48
49 /// Perform the per-event analysis
50 void analyze(const Event& event) {
51 // loop over particles
52 DecayedParticles BB = apply<DecayedParticles>(event, "BB");
53 int imode = -1;
54 for (unsigned int ix=0; ix<BB.decaying().size(); ++ix) {
55 int sign = BB.decaying()[ix].pid()/BB.decaying()[ix].abspid();
56 if (BB.decaying()[ix].abspid()==511) _c[0]->fill();
57 else _c[1]->fill();
58 if ((sign== 1 && BB.modeMatches(ix,5,mode1) ) ||
59 (sign==-1 && BB.modeMatches(ix,5,mode1CC) )) {
60 imode=0;
61 }
62 else if ((sign== 1 && BB.modeMatches(ix,5,mode2) ) ||
63 (sign==-1 && BB.modeMatches(ix,5,mode2CC) )) {
64 imode=1;
65 }
66 else if ((sign== 1 && BB.modeMatches(ix,5,mode3) ) ||
67 (sign==-1 && BB.modeMatches(ix,5,mode3CC)))
68 imode=2;
69 else if ((sign== 1 && BB.modeMatches(ix,5,mode4) ) ||
70 (sign==-1 && BB.modeMatches(ix,5,mode4CC) )) {
71 imode=3;
72 }
73 else {
74 continue;
75 }
76 const Particles& pip = BB.decayProducts()[ix].at( sign*211);
77 const Particle & pim = BB.decayProducts()[ix].at(-sign*211)[0];
78 const Particle & pi0 = BB.decayProducts()[ix].at( 111)[0];
79 FourMomentum pOmegaPi = pip[0].mom()+pip[1].mom()+pim.mom()+pi0.mom();
80 const double mHad = pOmegaPi.mass();
81 if (imode==0) {
82 _h[0]->fill(mHad);
83 }
84 // find the children of the omega
85 Particles omegaDec;
86 for (const Particle& p : {pip[0],pip[1],pim,pi0} ) {
87 Particle parent = p;
88 while (parent.parents()[0].pid()!=BB.decaying()[ix].pid()) {
89 parent = parent.parents()[0];
90 if (parent.pid()==223) {
91 omegaDec.push_back(p);
92 break;
93 }
94 }
95 }
96 if (omegaDec.size()!=3) continue;
97 if (imode<2) _h[imode+1]->fill(mHad);
98 else _h[5 ]->fill(mHad);
99 if (imode==1) continue;
100 _h_sum->fill(mHad);
101 // boost to B rest frame
102 LorentzTransform boostB = LorentzTransform::mkFrameTransformFromBeta(BB.decaying()[ix].mom().betaVec());
103 /// D star angles
104 if (imode==0) {
105 const Particle& Dstar = BB.decayProducts()[ix].at(-sign*413)[0];
106 if (Dstar.children().size()==2) {
107 Particle D0;
108 if (Dstar.children()[0].pid() ==-sign*211 &&
109 Dstar.children()[1].abspid() ==-sign*421) {
110 D0 = Dstar.children()[1];
111 }
112 else if (Dstar.children()[1].pid() ==-sign*211 &&
113 Dstar.children()[0].abspid() ==-sign*421) {
114 D0 = Dstar.children()[0];
115 }
116 // if right decay mode
117 if (D0.abspid()==421) {
118 FourMomentum pDstar = boostB.transform(Dstar.mom());
119 FourMomentum pD0 = boostB.transform(D0 .mom());
120 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pDstar.betaVec());
121 pD0 = boost2.transform(pD0);
122 double c1 = pD0.p3().unit().dot(pDstar.p3().unit());
123 _h[3]->fill(c1);
124 }
125 }
126 }
127 // omega momenta in B rest frame
128 FourMomentum pOmega;
129 for (const Particle & p : omegaDec) pOmega+=p.mom();
130 pOmega = boostB.transform(pOmega);
131 // boost to A rest frame
132 pOmegaPi = boostB.transform(pOmegaPi);
133 FourMomentum pDstar = boostB.transform(BB.decaying()[ix].mom()) - pOmegaPi;
134 LorentzTransform boostWpi = LorentzTransform::mkFrameTransformFromBeta(pOmegaPi.betaVec());
135 pOmega = boostWpi.transform(pOmega);
136 Vector3 axisW = pOmega.p3().unit();
137 Vector3 axisWpi = pOmegaPi.p3().unit();
138 double cA = axisW.dot(axisWpi);
139 // omega angles
140 LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
141 pOmegaPi = boostW.transform(boostWpi.transform(pOmegaPi));
142 pDstar = boostW.transform(boostWpi.transform(pDstar));
143 axisWpi = pOmegaPi.p3().unit();
144 FourMomentum ppip = boostW.transform(boostWpi.transform(boostB.transform(omegaDec[0])));
145 FourMomentum ppim = boostW.transform(boostWpi.transform(boostB.transform(omegaDec[1])));
146 Vector3 nW = ppip.p3().cross(ppim.p3()).unit();
147 double cTheta1 = axisWpi.dot(nW);
148 if (imode==0) {
149 _h[4]->fill(cTheta1);
150 continue;
151 }
152 _h_angle[0]->fill(cA);
153 _h_angle[1]->fill(cTheta1);
154 Vector3 transW = nW-cTheta1*axisWpi;
155 Vector3 transD = pDstar.p3().unit()-pDstar.p3().unit().dot(axisWpi)*axisWpi;
156 double phi = abs(atan(transW.cross(transD).dot(axisWpi)/ transW.dot(transD)));
157 _h_angle[2]->fill(phi);
158 }
159 }
160
161
162 /// Normalise histograms etc., after the run
163 void finalize() {
164 // first hist is differential BR
165 scale(_h[0], 1./ *_c[0]);
166 // rest are unit normalized
167 for(unsigned int ix=1;ix<6;++ix)
168 normalize(_h[ix], 1.0, false);
169 normalize(_h_angle, 1.0, false);
170 normalize(_h_sum, 1.0, false);
171 }
172
173 /// @}
174
175
176 /// @name Histograms
177 /// @{
178 Histo1DPtr _h[6],_h_angle[3],_h_sum;
179 CounterPtr _c[2];
180 const map<PdgId,unsigned int> mode1 = { { -413,1}, { 211,2}, {-211,1}, {111,1}};
181 const map<PdgId,unsigned int> mode1CC = { { 413,1}, {-211,2}, { 211,1}, {111,1}};
182 const map<PdgId,unsigned int> mode2 = { { -423,1}, { 211,2}, {-211,1}, {111,1}};
183 const map<PdgId,unsigned int> mode2CC = { { 423,1}, {-211,2}, { 211,1}, {111,1}};
184 const map<PdgId,unsigned int> mode3 = { { -411,1}, { 211,2}, {-211,1}, {111,1}};
185 const map<PdgId,unsigned int> mode3CC = { { 411,1}, {-211,2}, { 211,1}, {111,1}};
186 const map<PdgId,unsigned int> mode4 = { { -421,1}, { 211,2}, {-211,1}, {111,1}};
187 const map<PdgId,unsigned int> mode4CC = { { 421,1}, {-211,2}, { 211,1}, {111,1}};
188 /// @}
189
190
191 };
192
193
194 RIVET_DECLARE_PLUGIN(CLEO_2001_I554175);
195
196}
|