Rivet analyses
Mass and angular distributions in B → D(*)π+π−π−π0 decays
Experiment: CLEO (CESR)
Inspire ID: 554175
Status: VALIDATED NOHEPDATA
Authors: - Peter Richardson
References: - Phys.Rev.D 64 (2001) 092001
Beams: * *
Beam energies: ANY
Run details: - Any process producing B mesons, originally Upsilon(4S) decays
Measurement of mass and angular distributions in B → D(*)π+π−π−π0 decays, primarily B → D(*)ωπ−. The background subtracted data were read from the figures in the paper.
Source
code:CLEO_2001_I554175.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DecayedParticles.hh"
#include "Rivet/Projections/UnstableParticles.hh"
namespace Rivet {
/// @brief B -> D(*) pi+pi-pi-pi0
class CLEO_2001_I554175 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2001_I554175);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
UnstableParticles ufs = UnstableParticles(Cuts::abspid == 511 || Cuts::abspid == 521);
declare(ufs, "UFS");
DecayedParticles BB(ufs);
BB.addStable(PID::PI0);
BB.addStable(413);
BB.addStable(-413);
BB.addStable(423);
BB.addStable(-423);
BB.addStable(411);
BB.addStable(-411);
BB.addStable(421);
BB.addStable(-421);
declare(BB, "BB");
// // histos
for (unsigned int ix = 0; ix < 6; ++ix) {
book(_h[ix], 1 + ix, 1, 1);
}
for (unsigned int ix = 0; ix < 3; ++ix) {
book(_h_angle[ix], 7, 1, 1 + ix);
}
book(_h_sum, 8, 1, 1);
for (unsigned int ix = 0; ix < 2; ++ix) {
book(_c[ix], "TMP/nB_" + toString(ix + 1));
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// loop over particles
DecayedParticles BB = apply<DecayedParticles>(event, "BB");
int imode = -1;
for (unsigned int ix = 0; ix < BB.decaying().size(); ++ix) {
int sign = BB.decaying()[ix].pid() / BB.decaying()[ix].abspid();
if (BB.decaying()[ix].abspid() == 511)
_c[0]->fill();
else
_c[1]->fill();
if ((sign == 1 && BB.modeMatches(ix, 5, mode1)) || (sign == -1 && BB.modeMatches(ix, 5, mode1CC))) {
imode = 0;
}
else if ((sign == 1 && BB.modeMatches(ix, 5, mode2))
|| (sign == -1 && BB.modeMatches(ix, 5, mode2CC))) {
imode = 1;
}
else if ((sign == 1 && BB.modeMatches(ix, 5, mode3))
|| (sign == -1 && BB.modeMatches(ix, 5, mode3CC)))
imode = 2;
else if ((sign == 1 && BB.modeMatches(ix, 5, mode4))
|| (sign == -1 && BB.modeMatches(ix, 5, mode4CC))) {
imode = 3;
}
else {
continue;
}
const Particles& pip = BB.decayProducts()[ix].at(sign * 211);
const Particle& pim = BB.decayProducts()[ix].at(-sign * 211)[0];
const Particle& pi0 = BB.decayProducts()[ix].at(111)[0];
FourMomentum pOmegaPi = pip[0].mom() + pip[1].mom() + pim.mom() + pi0.mom();
const double mHad = pOmegaPi.mass();
if (imode == 0) {
_h[0]->fill(mHad);
}
// find the children of the omega
Particles omegaDec;
for (const Particle& p : {pip[0], pip[1], pim, pi0}) {
Particle parent = p;
while (parent.parents()[0].pid() != BB.decaying()[ix].pid()) {
parent = parent.parents()[0];
if (parent.pid() == 223) {
omegaDec.push_back(p);
break;
}
}
}
if (omegaDec.size() != 3) continue;
if (imode < 2)
_h[imode + 1]->fill(mHad);
else
_h[5]->fill(mHad);
if (imode == 1) continue;
_h_sum->fill(mHad);
// boost to B rest frame
LorentzTransform boostB = LorentzTransform::mkFrameTransformFromBeta(
BB.decaying()[ix].mom().betaVec());
/// D star angles
if (imode == 0) {
const Particle& Dstar = BB.decayProducts()[ix].at(-sign * 413)[0];
if (Dstar.children().size() == 2) {
Particle D0;
if (Dstar.children()[0].pid() == -sign * 211 && Dstar.children()[1].abspid() == -sign * 421) {
D0 = Dstar.children()[1];
}
else if (Dstar.children()[1].pid() == -sign * 211
&& Dstar.children()[0].abspid() == -sign * 421) {
D0 = Dstar.children()[0];
}
// if right decay mode
if (D0.abspid() == 421) {
FourMomentum pDstar = boostB.transform(Dstar.mom());
FourMomentum pD0 = boostB.transform(D0.mom());
LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pDstar.betaVec());
pD0 = boost2.transform(pD0);
double c1 = pD0.p3().unit().dot(pDstar.p3().unit());
_h[3]->fill(c1);
}
}
}
// omega momenta in B rest frame
FourMomentum pOmega;
for (const Particle& p : omegaDec) pOmega += p.mom();
pOmega = boostB.transform(pOmega);
// boost to A rest frame
pOmegaPi = boostB.transform(pOmegaPi);
FourMomentum pDstar = boostB.transform(BB.decaying()[ix].mom()) - pOmegaPi;
LorentzTransform boostWpi = LorentzTransform::mkFrameTransformFromBeta(pOmegaPi.betaVec());
pOmega = boostWpi.transform(pOmega);
Vector3 axisW = pOmega.p3().unit();
Vector3 axisWpi = pOmegaPi.p3().unit();
double cA = axisW.dot(axisWpi);
// omega angles
LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
pOmegaPi = boostW.transform(boostWpi.transform(pOmegaPi));
pDstar = boostW.transform(boostWpi.transform(pDstar));
axisWpi = pOmegaPi.p3().unit();
FourMomentum ppip = boostW.transform(boostWpi.transform(boostB.transform(omegaDec[0])));
FourMomentum ppim = boostW.transform(boostWpi.transform(boostB.transform(omegaDec[1])));
Vector3 nW = ppip.p3().cross(ppim.p3()).unit();
double cTheta1 = axisWpi.dot(nW);
if (imode == 0) {
_h[4]->fill(cTheta1);
continue;
}
_h_angle[0]->fill(cA);
_h_angle[1]->fill(cTheta1);
Vector3 transW = nW - cTheta1 * axisWpi;
Vector3 transD = pDstar.p3().unit() - pDstar.p3().unit().dot(axisWpi) * axisWpi;
double phi = abs(atan(transW.cross(transD).dot(axisWpi) / transW.dot(transD)));
_h_angle[2]->fill(phi);
}
}
/// Normalise histograms etc., after the run
void finalize() {
// first hist is differential BR
scale(_h[0], 1. / *_c[0]);
// rest are unit normalized
for (unsigned int ix = 1; ix < 6; ++ix) normalize(_h[ix], 1.0, false);
normalize(_h_angle, 1.0, false);
normalize(_h_sum, 1.0, false);
}
/// @}
/// @name Histograms
/// @{
Histo1DPtr _h[6], _h_angle[3], _h_sum;
CounterPtr _c[2];
const map<PdgId, unsigned int> mode1 = {{-413, 1}, {211, 2}, {-211, 1}, {111, 1}};
const map<PdgId, unsigned int> mode1CC = {{413, 1}, {-211, 2}, {211, 1}, {111, 1}};
const map<PdgId, unsigned int> mode2 = {{-423, 1}, {211, 2}, {-211, 1}, {111, 1}};
const map<PdgId, unsigned int> mode2CC = {{423, 1}, {-211, 2}, {211, 1}, {111, 1}};
const map<PdgId, unsigned int> mode3 = {{-411, 1}, {211, 2}, {-211, 1}, {111, 1}};
const map<PdgId, unsigned int> mode3CC = {{411, 1}, {-211, 2}, {211, 1}, {111, 1}};
const map<PdgId, unsigned int> mode4 = {{-421, 1}, {211, 2}, {-211, 1}, {111, 1}};
const map<PdgId, unsigned int> mode4CC = {{421, 1}, {-211, 2}, {211, 1}, {111, 1}};
/// @}
};
RIVET_DECLARE_PLUGIN(CLEO_2001_I554175);
}