umd-lhcb / lhcb-ntuples-gen

ntuples generation with DaVinci and in-house offline components
BSD 2-Clause "Simplified" License
1 stars 0 forks source link

Rest frame approximation #10

Closed yipengsun closed 5 years ago

yipengsun commented 5 years ago

We should try to include RFA in our first-stage ntuple generation, using SemileptonicCommonTools.

The tool that is responsible to RFA is called TupleToolTauMuDiscrVars, and can be imported as follows:

from Configurables import TupleToolTauMuDiscrVars

A simple documentation can be found at twiki: 1.

yipengsun commented 5 years ago

According to @CoffeeIntoScience, we just need to add this TupleTool to the YCands.

yipengsun commented 5 years ago

After using the tuple tool above, in the generated ntuples, I see the following variables: Y_FlightDir_Zanlge, Y_Estimated_P, Y_Estimated_PX, Y_Estimated_PY, Y_Estimated_PZ, Y_Boost_Beta.

yipengsun commented 5 years ago

We need to compute the missing mass based on these additional variables (which are supposedly the reconstructed momentum of the B mesons) and check that these are done correctly.

Just to remind myself, in RFA, we take the missing neutrinos into account and try to reconstruct B mesons momentum.

CoffeeIntoScience commented 5 years ago

Ok. First dynamics then logistics.

I think I sent you this already, but to do the computation its best to not reinvent the wheel and to make judicious use of the TLorentzVector classes (or the HepLorentzVector in CLHEP libraries if you're allergic to ROOT).

You can do things like TLorentzVector p4_B; p4_B.SetPxPyPzM(Y_Estimated_PX, Y_Estimated_Py, Y_Estimated_Pz); TVector3 boost_B=p4_B.BoostVector();

which gets you the boost that gets you from the B frame to the lab frame (i.e. p4_whatever.Boost( -1* boost_B) gives you p4_whatever in the B frame)

Then you just need to remember what kinematics you need and can easily do things like q2 = (p4_B - p4_Dst).M2()

CoffeeIntoScience commented 5 years ago

Now logistics

For this you're going to need to step through the tuples row by row and compute the fit variables. As always there a bunch of ways to do this with the ROOT classes.

The fastest way to get started is to use TTree::MakeSelector(char) to build a skeleton class to load the tree and provide methods that are called before starting the loop, event by event, and when the loop terminates.

You also are going to want to save this output in some way. There are a couple ways to do this. I put them all in a new tree along with relevant info for final sorting into signal/control samples and cuts that were not final and so and and so on. This got very bloated quickly. I think there's a way to add them as new leaves to the original tree. This is dangerous enough you should have a clean backup copy of your tuples. A third interesting option is to do a new tree with these variables, but only also save the EventNumber and RunNumber. You can then make this new tree a friend of the old tree and cross reference info for cuts in the original tree with the computed fit variables in the new tree.

yipengsun commented 5 years ago

The script snippet you sent me contains both how to do RFA and how to compute the missing mass with these additional variables:

TLorentzVector p4B, p4Y, p4mu p4c;
p4Y.SetXYZT(Y_PX,Y_PY,Y_PZ,Y_PE); //D(*)mu combination 4-momentum
p4c.SetXYZT(D0_PX,D0_PY,D0_PZ,D0_PE); //charm 4 momentum. D or D* (or Lc or=
 whatever)
p4mu.SetXYZT(muplus_PX,muplus_PY,muplus_PZ,muplus_PE); //muon 4-momentum

TVector3 flight;
flight.SetXYZ(Y_ENDVERTEX_X-Y_OWNPV_X, Y_ENDVERTEX_Y-Y_OWNPV_Y, Y_ENDVERTEX=
_Z-Y_OWNPV_Z); //B flight dir

double tantheta =3D (flight.Unit()).Perp()/((flight.Unit()).Z());  //needed=
 below.

/*here is the actual approximation itself.
  it's easier to work with pt than longitudinal momentum
  estimate because of existing TLorentzVector methods.
  b_m is just a constant holding the B mass in MeV
*/
double pt_estimate =3D (b_m/Y_M)*tantheta*Y_PZ;

p4B.SetPtEtaPhiM(pt_estimate,flight.Eta(),flight.Phi(),b_m);

//missing mass
double mm2 =3D (p4B-p4Y).M2();

//E_{#mu}
TLorentzVector p4muCM;
p4muCM=3Dp4mu;
p4muCM.Boost(-1*p4B.BoostVector());
double Emu =3D p4muCM.E();

//q^{2}
double q2 =3D (p4B-p4c).M2();
CoffeeIntoScience commented 5 years ago

Going back to dynamics, for the fit variables recall that missing mass and q2 are lorentz-invariant quantities so if you're worrying about boosts you are making life harder than it ought to be. E_\ell is the only one you need to worry about which frame you're in.

But you'll want to start getting the hang of it now because the theoretical description of the decays involves things like "the angle in lepton-neutrino center of mass frame between the lepton and the opposite of the Dst momentum" where you have to boost first to the B frame then to this other frame. (boosts don't commute!)

CoffeeIntoScience commented 5 years ago

The script you sent me contains both how to do RFA and how to compute the missing mass with these additional variables:

It does. Be careful though. I sent you some code, not a complete script. It won't actually run in any meaningful way without a good amount of other code around it

yipengsun commented 5 years ago

I'm trying to compute the missing mass and q^2 with DaVinci. So far there might be two routes:

Personally I prefer the LoKi route, because according to the LHCb starter kit 1, we can also add Loki functors as tuple tools, and it seems way easier.

yipengsun commented 5 years ago

1 and 2 provide examples on using LoKi functors as TupleTool. But I'm no longer sure if LoKi functors can access variables created by the TupleToolTauMuDiscrVars. Hopefully we can use the VALUE functor to extract these values.

At the end of 3, there's a brief instruction on writing our own TupleTool.

yipengsun commented 5 years ago

Giving up on the LoKi route, due to the inability to find documentation on the usage of VALUE functor.

yipengsun commented 5 years ago

From the TupleToolTauMuDiscrVars.cpp:

    // Filling of the branches.
    test &= tuple->column( head+"_FlightDir_Zangle" , B_angle_Z );
    test &= tuple->column( head+"_Estimated_P" , B_P );
    test &= tuple->column( head+"_Estimated_PX" , B_PXYZE.X() );
    test &= tuple->column( head+"_Estimated_PY" , B_PXYZE.Y() );
    test &= tuple->column( head+"_Estimated_PZ" , B_PXYZE.Z() );
    test &= tuple->column( head+"_Boost_Beta" , B_beta );
    test &= tuple->column( "FitVar_El" , El );
    test &= tuple->column( "FitVar_Mmiss2" , Mmiss2 );
    test &= tuple->column( "FitVar_q2" , q2 );

The missing mass and q2 should already present in the ntuple.

yipengsun commented 5 years ago

I can confirm that FitVar_El, FitVar_Mmiss2, FitVar_q2 are already in the ntuple.

yipengsun commented 5 years ago

Here I convince myself why we need the B_M/observed_B_M:

    // Bd estimated momentum.
    const double B_M = 5279.61;
    const double observed_B_M = P->measuredMass();
    const double observed_B_PZ = P->momentum().Z();
    const double B_P = B_M/observed_B_M*observed_B_PZ*::sqrt(1.+::tan(B_angle_Z)*::tan(B_angle_Z));
    const LoKi::LorentzVector B_PXYZE (B_P*B_cosangle_X,B_P*B_cosangle_Y,B_P*B_cosangle_Z,::sqrt(B_P*B_P+B_M*B_M));

Comparing with @CoffeeIntoScience's snippet:

double pt_estimate =3D (b_m/Y_M)*tantheta*Y_PZ;

Note that In @CoffeeIntoScience's code, Y_M is the mass of the D*-lepton system, not the B.

I claim: In the official TupleTool, the observed_B_M is equivalent to Y_M. observed_B_M is probably a misnomer just a name.

@CoffeeIntoScience Do you agree with the claim above?

CoffeeIntoScience commented 5 years ago

eh, i'd avoid calling anything a misnomer. my choice to rename the D(st)mu combination to not call it a B is personal preference, and Y is anything but standard. But yes, oberved_B_M = m_visibile = Y_M