<dBR/dq2>(B0->K*ee) near threshold

[fdesse]

Dear all,
I am witnessing some strange behaviour of the observable <dBR/dq2>(B0->K*ee) near the threshold.
I understand that binned differential branching ratio in Flavio are normalized to the bin width and thus have units of GeV^(-2).

Very naively, I thought that if I want the "unormalized" differential branching fraction, I just have to multiply the result by the bin width. This seems to work fine, as long as I stay away from the threshold:

BR_0004_1 = flavio.sm_prediction("<dBR/dq2>(B0->K*ee)", .0004, 1.)*(1. - .0004)
BR_1_3 = flavio.sm_prediction("<dBR/dq2>(B0->K*ee)", 1., 3.)*(3. - 1.)
BR_0004_3 = flavio.sm_prediction("<dBR/dq2>(B0->K*ee)", .0004, 3.)*(3. - .0004)

gives:

In [81]: BR_0004_3
Out[81]: 3.790997798081069e-07

In [82]: BR_0004_1 + BR_1_3
Out[82]: 3.7912788564241227e-07

That seems reasonably close.

But near the threshold:

BR_th_1 = flavio.sm_prediction("<dBR/dq2>(B0->K*ee)", (2*.511/1000.)**2, 1.)*(1. - (2*.511/1000.)**2)
BR_1_3 = flavio.sm_prediction("<dBR/dq2>(B0->K*ee)", 1., 3.)*(3. - 1.)
BR_th_3 = flavio.sm_prediction("<dBR/dq2>(B0->K*ee)", (2*.511/1000.)**2, 3.)*(3. - (2*.511/1000.)**2)

gives:

In [86]: BR_th_3
Out[86]: 3.642236984062088e-07

In [87]: BR_th_1 + BR_1_3
Out[87]: 4.093330290405161e-07

I am probably doing something very silly, but I don't understand what. Thanks in advance !

[DavidMStraub]

Right, this should have worked ...

I looked into it and it seems the scipy.integrate.quadrature integrator fails by a large margin inspite of the fact that the relative tolerance is set to 5e-3. That's very suprising/annoying. I checked that using scipy.integrate.quad one gets the correct result, but I don't want to use that everywhere since it's often too slow. I'll have to think about it. In the meantime, as a workaround you can just use an integrator of your choice and directly integrate the function

lambda q2: flavio.sm_prediction('dBR/dq2(B+->Kee)', q2)

However note that you are really looking at a pathological case here: you are integrating numerically near a pole. If one looks at realistic bins used in experiments, everything is fine. Take for instance http://arxiv.org/pdf/1501.03038.pdf which has q2min=0.002: in this case your check works to excellent accuracy. Same goes for the RK* measurement with q2min=0.045.

[fdesse]

Thanks for your complete answer David.

Yes indeed, experimentally, we are not going that low in q2. But current analyses at LHCb aim at going down to q2min=0.0001, which, if you take into account the leakage (even if you cut at reconstructed q2min, you might have some events with lower true q2) we dangerously approach the limit.

In the meantime, I will use the workaround you are proposing :)

[fdesse]

Hi !
I am still a bit confused. I tried to use scipy.integrate.quad as you suggested.

dBR = lambda q2: flavio.sm_prediction('dBR/dq2(B0->K*ee)', q2) 
def BR(q2min, q2max):
     return scipy.integrate.quad(dBR, q2min, q2max)[0]*(q2max - q2min)

BRquad_th_1 = BR((2*.511/1000.)**2, 1.)
BRquad_1_3  = BR(1., 3.)
BRquad_th_3 = BR((2*.511/1000.)**2, 3.)

I get:

In [65]: BRquad_th_3
Out[65]: 1.7111898035244934e-06

In [66]: BRquad_th_1 + BRquad_1_3
Out[66]: 6.646987892392002e-07

Seems to me that its even worse. As a workaround, I now use .000025 as a lower limit. This seems to work more or less fine.

[DavidMStraub]

When you're integrating the differential branching ratio (units 1/GeV²) over q² (units GeV²), of course you should not divide by the bin width anymore. Remove *(q2max - q2min) and you get the right result.

[fdesse]

Ah but of course !! Silly me. Sorry for the noise.