Rivet analyses
Decay kinematics of semileptonic B0 → D*− decays.
Experiment: BELLE (KEKB)
Inspire ID: 1693396
Status: VALIDATED NOHEPDATA
Authors: - Peter Richardson
References: - Phys.Rev.D 100 (2019) 5, 052007
Beams: * *
Beam energies: ANY
Run details: - Any process producing B0
Measurement of recoil w, helicity and decay plane angles of semileptonc B̄0 to D*+ decays. N.B. the data is not unfolded but the efficiencies and response matrices used in the paper are applied to the results of the simulation.
Source
code:BELLE_2019_I1693396.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/UnstableParticles.hh"
namespace Rivet {
/// @brief B0 -> D*- semileptonic
class BELLE_2019_I1693396 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2019_I1693396);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections
declare(UnstableParticles(Cuts::pid==511), "UFS");
// Book histograms
for(unsigned int ix=0;ix<2;++ix) {
book(_h[0][ix], "TMP/h_w_" +toString(ix+1), refData(1, 1, ix+1));
book(_h[1][ix], "TMP/h_costhl_"+toString(ix+1), refData(3, 1, ix+1));
book(_h[2][ix], "TMP/h_costhv_"+toString(ix+1), refData(2, 1, ix+1));
book(_h[3][ix], "TMP/h_chi_" +toString(ix+1), refData(4, 1, ix+1));
}
}
/// Perform the per-event analysis
bool analyzeDecay(Particle mother, vector<int> ids) {
// There is no point in looking for decays with less particles than to be analysed
if (mother.children().size() == ids.size()) {
bool decayfound = true;
for (int id : ids) {
if (!contains(mother, id)) decayfound = false;
}
return decayfound;
}
return false;
}
bool contains(Particle& mother, int id) {
return any(mother.children(), HasPID(id));
}
double recoilW(const Particle& mother) {
FourMomentum lepton, neutrino, meson, q;
for (const Particle& c : mother.children()) {
if (c.isNeutrino()) neutrino = c.mom();
else if (c.isChargedLepton()) lepton = c.mom();
else if (c.isHadron()) meson = c.mom();
}
q = lepton + neutrino; //no hadron before
double mb2= mother.mom()*mother.mom();
double mD2 = meson*meson;
return (mb2 + mD2 - q*q )/ (2. * sqrt(mb2) * sqrt(mD2) );
}
/// Perform the per-event analysis
void analyze(const Event& event) {
FourMomentum pl, pnu, pB, pD, pDs, ppi;
// Iterate of B0bar mesons
for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
pB = p.momentum();
// Find semileptonic decays
int iloc=-1;
if (analyzeDecay(p, {PID::DSTARMINUS,12,-11})) iloc = 0;
else if(analyzeDecay(p, {PID::DSTARMINUS,14,-13}) ) iloc = 1;
else continue;
_h[0][iloc]->fill(recoilW(p));
// Get the necessary momenta for the angles
bool foundDdecay=false;
for (const Particle & c : p.children()) {
if ( (c.pid() == PID::DSTARMINUS) && (analyzeDecay(c, {PID::PIMINUS, PID::D0BAR}) || analyzeDecay(c, {PID::PI0, PID::DMINUS})) ) {
foundDdecay=true;
pDs = c.momentum();
for (const Particle & dc : c.children()) {
if (dc.hasCharm()) pD = dc.momentum();
else ppi = dc.momentum();
}
}
if (c.pid() == -11 || c.pid() == -13) pl = c.momentum();
if (c.pid() == 12 || c.pid() == 14) pnu = c.momentum();
}
// This is the angle analysis
if (!foundDdecay) continue;
// First boost all relevant momenta into the B-rest frame
const LorentzTransform B_boost = LorentzTransform::mkFrameTransformFromBeta(pB.betaVec());
// Momenta in B rest frame:
FourMomentum lv_brest_Dstar = B_boost.transform(pDs);
FourMomentum lv_brest_w = B_boost.transform(pB - pDs);
FourMomentum lv_brest_D = B_boost.transform(pD);
FourMomentum lv_brest_lep = B_boost.transform(pl);
const LorentzTransform Ds_boost = LorentzTransform::mkFrameTransformFromBeta(lv_brest_Dstar.betaVec());
FourMomentum lv_Dstarrest_D = Ds_boost.transform(lv_brest_D);
const LorentzTransform W_boost = LorentzTransform::mkFrameTransformFromBeta(lv_brest_w.betaVec());
FourMomentum lv_wrest_lep = W_boost.transform(lv_brest_lep);
double cos_thetaV = cos(lv_brest_Dstar.p3().angle(lv_Dstarrest_D.p3()));
_h[2][iloc]->fill(cos_thetaV);
double cos_thetaL = cos(lv_brest_w.p3().angle(lv_wrest_lep.p3()));
_h[1][iloc]->fill(cos_thetaL);
Vector3 LTrans = lv_wrest_lep.p3() - cos_thetaL*lv_wrest_lep.p3().perp()*lv_brest_w.p3().unit();
Vector3 VTrans = lv_Dstarrest_D.p3() - cos_thetaV*lv_Dstarrest_D.p3().perp()*lv_brest_Dstar.p3().unit();
float chi = atan2(LTrans.cross(VTrans).dot(lv_brest_w.p3().unit()), LTrans.dot(VTrans));
_h[3][iloc]->fill(chi);
}
}
/// Normalise histograms etc., after the run
void finalize() {
// efficiencies
vector<double> eff[4][2]={{{2.72,5.72,7.7,9.1,10.03,10.61,10.74,10.67,10.23,9.1},
{2.68,5.66,7.66,9.05,9.91,10.43,10.6,10.52,10.04,9.14}},
{{3.12,3.97,5.73,7.96,9.31,9.85,10.23,10.59,11.06,11.21},
{3.16,3.52,5.19,7.59,9.1,9.78,10.27,10.43,11.0,11.36}},
{{11.72,11.52,11.35,10.88,10.2,9.34,8.29,7.16,6.05,4.82},
{11.54,11.43,11.14,10.74,10.09,9.29,8.25,7.1,5.97,4.72}},
{{8.6,8.74,8.96,9.3,9.81,9.82,9.33,9.0,8.77,8.59},
{8.51,8.67,8.82,9.15,9.7,9.73,9.2,8.83,8.62,8.54}}};
vector<double> efe[4][2]={{{0.02,0.02,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03},
{0.02,0.02,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03}},
{{0.03,0.02,0.02,0.03,0.03,0.03,0.03,0.03,0.03,0.03},
{0.03,0.02,0.02,0.03,0.03,0.03,0.03,0.03,0.03,0.03}},
{{0.03,0.03,0.03,0.04,0.04,0.04,0.03,0.03,0.02,0.02},
{0.03,0.03,0.03,0.04,0.04,0.04,0.03,0.03,0.02,0.02}},
{{0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03},
{0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03}}};
// response matricesdouble
double response[4][2][10][10] = {{{{0.803,0.053,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.197,0.778,0.098,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.168,0.717,0.126,0.002,0.0,0.0,0.0,0.0,0.0},
{0.0,0.001,0.182,0.667,0.149,0.006,0.0,0.0,0.0,0.0},
{0.0,0.0,0.004,0.199,0.626,0.167,0.011,0.0,0.0,0.0},
{0.0,0.0,0.0,0.009,0.207,0.592,0.177,0.015,0.0,0.0},
{0.0,0.0,0.0,0.0,0.016,0.215,0.575,0.183,0.018,0.0},
{0.0,0.0,0.0,0.0,0.0,0.021,0.213,0.567,0.186,0.017},
{0.0,0.0,0.0,0.0,0.0,0.0,0.024,0.214,0.598,0.186},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.022,0.198,0.797}},
{{0.961,0.024,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.038,0.952,0.027,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.021,0.948,0.041,0.001,0.0,0.0,0.0,0.0,0.0},
{0.0,0.001,0.023,0.918,0.067,0.003,0.001,0.001,0.001,0.0},
{0.0,0.001,0.001,0.04,0.871,0.097,0.005,0.001,0.001,0.0},
{0.0,0.0,0.0,0.001,0.06,0.817,0.129,0.006,0.001,0.0},
{0.0,0.0,0.0,0.0,0.001,0.082,0.758,0.164,0.007,0.001},
{0.0,0.0,0.0,0.0,0.0,0.001,0.106,0.698,0.196,0.008},
{0.0,0.0,0.0,0.0,0.0,0.0,0.001,0.128,0.657,0.212},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.002,0.137,0.777}}},
{{{0.918,0.077,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.082,0.806,0.095,0.001,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.115,0.761,0.101,0.002,0.0,0.0,0.0,0.0,0.0},
{0.0,0.001,0.141,0.735,0.105,0.002,0.0,0.0,0.0,0.0},
{0.0,0.0,0.002,0.16,0.719,0.1,0.001,0.0,0.0,0.0},
{0.0,0.0,0.0,0.003,0.17,0.722,0.093,0.001,0.0,0.0},
{0.0,0.0,0.0,0.0,0.003,0.173,0.738,0.08,0.001,0.0},
{0.0,0.0,0.0,0.0,0.0,0.002,0.166,0.771,0.072,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.001,0.147,0.819,0.064},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.001,0.108,0.936}},
{{0.659,0.129,0.011,0.003,0.002,0.002,0.002,0.004,0.013,0.144},
{0.151,0.691,0.132,0.012,0.004,0.002,0.002,0.002,0.004,0.016},
{0.015,0.141,0.697,0.147,0.016,0.005,0.002,0.002,0.002,0.005},
{0.005,0.012,0.134,0.671,0.162,0.018,0.005,0.002,0.002,0.002},
{0.002,0.004,0.013,0.14,0.634,0.155,0.016,0.004,0.002,0.002},
{0.002,0.002,0.004,0.015,0.155,0.633,0.141,0.013,0.004,0.003},
{0.002,0.002,0.002,0.004,0.018,0.163,0.67,0.136,0.012,0.004},
{0.005,0.002,0.002,0.002,0.005,0.015,0.147,0.695,0.14,0.015},
{0.016,0.004,0.002,0.002,0.002,0.004,0.013,0.132,0.691,0.15},
{0.142,0.013,0.003,0.002,0.002,0.002,0.003,0.012,0.13,0.659}}},
{{{0.812,0.051,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.188,0.784,0.096,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.164,0.728,0.126,0.002,0.0,0.0,0.0,0.0,0.0},
{0.0,0.001,0.172,0.676,0.149,0.006,0.0,0.0,0.0,0.0},
{0.0,0.0,0.004,0.19,0.631,0.165,0.01,0.0,0.0,0.0},
{0.0,0.0,0.0,0.008,0.203,0.6,0.181,0.016,0.0,0.0},
{0.0,0.0,0.0,0.0,0.014,0.209,0.578,0.187,0.019,0.0},
{0.0,0.0,0.0,0.0,0.0,0.02,0.209,0.573,0.195,0.017},
{0.0,0.0,0.0,0.0,0.0,0.0,0.022,0.205,0.6,0.195},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.019,0.186,0.788}},
{{0.959,0.022,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.039,0.955,0.012,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.021,0.96,0.022,0.001,0.0,0.0,0.0,0.0,0.0},
{0.001,0.001,0.026,0.931,0.043,0.001,0.0,0.0,0.0,0.0},
{0.0,0.0,0.001,0.047,0.889,0.07,0.002,0.0,0.0,0.0},
{0.0,0.001,0.0,0.0,0.067,0.837,0.103,0.002,0.001,0.0},
{0.0,0.0,0.0,0.0,0.0,0.091,0.778,0.138,0.003,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.117,0.715,0.174,0.004},
{0.0,0.0,0.0,0.0,0.0,0.0,0.001,0.142,0.672,0.193},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.002,0.151,0.803}}},
{{{0.918,0.077,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.082,0.805,0.091,0.001,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.117,0.763,0.101,0.002,0.0,0.0,0.0,0.0,0.0},
{0.0,0.001,0.142,0.735,0.103,0.002,0.0,0.0,0.0,0.0},
{0.0,0.0,0.003,0.159,0.723,0.098,0.001,0.0,0.0,0.0},
{0.0,0.0,0.0,0.004,0.169,0.726,0.091,0.001,0.0,0.0},
{0.0,0.0,0.0,0.0,0.004,0.172,0.745,0.082,0.001,0.0},
{0.0,0.0,0.0,0.0,0.0,0.002,0.161,0.771,0.074,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.001,0.145,0.817,0.066},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.107,0.934}},
{{0.653,0.129,0.012,0.004,0.003,0.002,0.002,0.004,0.014,0.144},
{0.152,0.686,0.13,0.013,0.004,0.003,0.002,0.002,0.005,0.017},
{0.016,0.143,0.693,0.147,0.016,0.006,0.003,0.002,0.003,0.005},
{0.005,0.013,0.138,0.667,0.16,0.018,0.005,0.002,0.002,0.003},
{0.003,0.004,0.013,0.142,0.63,0.156,0.015,0.004,0.002,0.002},
{0.002,0.002,0.004,0.015,0.158,0.629,0.142,0.013,0.004,0.003},
{0.003,0.002,0.002,0.005,0.018,0.164,0.667,0.138,0.013,0.005},
{0.005,0.003,0.002,0.003,0.006,0.016,0.148,0.692,0.141,0.016},
{0.017,0.004,0.002,0.002,0.003,0.005,0.013,0.131,0.686,0.152},
{0.144,0.014,0.004,0.002,0.002,0.002,0.004,0.012,0.129,0.654}}}};
// correct the values
for(unsigned int ix=0;ix<4;++ix) {
for(unsigned int iy=0;iy<2;++iy) {
Estimate1DPtr corrected;
book(corrected,ix+1,1,iy+1);
// first extract values and errors applying efficiency
Vector<10> val,err;
for(unsigned int ibin=0;ibin<_h[ix][iy]->bins().size();++ibin) {
val[ibin] = eff[ix][iy][ibin]/100. * _h[ix][iy]->bins()[ibin+1].sumW();
err[ibin] = sqr(eff[ix][iy][ibin]/100. * _h[ix][iy]->bins()[ibin+1].errW());
sqr(efe[ix][iy][ibin]/100. * _h[ix][iy]->bins()[ibin+1].sumW());
}
// put response into a matrix
Matrix<10> R,R2;
for(unsigned int i1=0;i1<10;++i1) {
for(unsigned int i2=0;i2<10;++i2) {
R .set(i1,i2, response[ix][iy][i1][i2]);
R2.set(i1,i2,sqr(response[ix][iy][i1][i2]));
}
}
// multiply to get value and error^2
val = multiply(R ,val);
err = multiply(R2,err);
// total for normalization
double total=0.;
for(unsigned int i1=0;i1<10;++i1) total+=val[i1];
// finally the output scatter
for(unsigned int ibin=0;ibin<_h[ix][iy]->bins().size();++ibin) {
double dx = 0.5*_h[ix][iy]->bins()[ibin+1].xWidth();
double dy = sqrt(err[ibin])/total/2./dx;
corrected->bin(ibin+1).set(val[ibin]/total/2./dx, dy);
}
}
}
}
/// @}
/// @name Histograms
/// @{
Histo1DPtr _h[4][2];
/// @}
};
RIVET_DECLARE_PLUGIN(BELLE_2019_I1693396);
}