Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DTF via LoKi functors #3

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 67 additions & 82 deletions 22-decay-tree-fitter.md
Original file line number Diff line number Diff line change
@@ -1,118 +1,103 @@
---
layout: page
title: HowTo DecayTreeFitter
subtitle: How do I use DecayTreeFitter?
subtitle: DecayTreeFitter via LoKi functors
minutes: 10
---

> ## Learning Objectives {.objectives}
>
> * Add a kinematic fitter to a branch in the decay tree
> * Apply a mass constraint
> * Inspect the refitted decay tree

Once you have made a hypothesis on the chain of decays that lead to your final state, you then can incorporate the additional knowledge that comes with this hypothesis to get a new best estimate for the particle parameters -- in particular their momenta. The additional knowledge is represented as constraints, which your decay tree has to fulfill.
Let us consider the ```B0 -> J/psi(1S) K*(892)0``` decay, with ```J/psi(1S) -> e+ e-``` and ```K*(892)0 -> K+ pi-```.

For example, for the decay
```python
'[D*(2010)+ -> (D0 -> K- pi+) pi+]CC'
```
you can make the assumption that the (K- pi+) combine to form a D0 with a specific invariant mass. This results in a so called *mass-constraint*. In addition the kaon and the pion should originate from exactly the same point in space. If you know that your data only contains prompt D* candidates, you can constrain them to do come from the primary vertex. Boundary conditions like those are called *vertex-constraints*.
We can make the assumption that the `B0` originates from the primary vertex (vertex constraint) and that the `e+` and `e-` combine to form a `J/psi(1S)` with a specific invariant mass (mass constraint). We then want to add the corresponding new best estimates of some physical variables to our ntuple, in order to be able to use it in the subsequent analysis.

Applying such kinematic constraints leads to new best estimates for the track parameters of the final state particles. The process of calculating those is called a *kinematic refit* and the `TupleToolDecayTreeFitter` is the algorithm that performs this task for us.
First of all, we need to specify the corresponding decay descriptor and branches:

> ## The physics and mathematics behind DecayTreeFitter {.callout}
> For details of the method see the paper on [Decay chain fitting with a Kalman filter](http://arxiv.org/abs/physics/0503191).

So how do we use a `TupleToolDecayTreeFitter` to our DaVinci script? Let's create a branch to add the tool to. We'll just name it `'Dstar'`:
```python
dtt.addBranches({
'Dstar': '[D*(2010)+ -> (D0 -> K- pi+) pi+]CC',
})
```
To this branch we can now apply the `TupleToolDecayTreeFitter`.
```python
dtt.Dstar.addTupleTool('TupleToolDecayTreeFitter/ConsD')
```
Now we can proceed with the configuration of the fitter. We are going to constrain the decay tree to the primary vertex of origin. We want all the output available, so we set the `verbose` option. Finally we want to apply the mass constraint on the D0:
```python
dtt.Dstar.ConsD.constrainToOriginVertex = True
dtt.Dstar.ConsD.Verbose = True
dtt.Dstar.ConsD.daughtersToConstrain = ['D0']
myDecay = '[[B0]cc -> ^(J/psi(1S) -> ^e+ ^e-) ^(K*(892)0 -> ^K+ ^pi-)]CC'
myBranches = {
"K" : "[ [B0]cc -> (J/psi(1S) -> e+ e-) (K*(892)0 -> ^K+ pi-)]CC",
"Pi" : "[ [B0]cc -> (J/psi(1S) -> e+ e-) (K*(892)0 -> K+ ^pi-)]CC",
"Lplus" : "[ [B0]cc -> (J/psi(1S) -> ^e+ e-) (K*(892)0 -> K+ pi-)]CC",
"Lminus" : "[ [B0]cc -> (J/psi(1S) -> e+ ^e-) (K*(892)0 -> K+ pi-)]CC",
"Jpsi" : "[ [B0]cc -> ^(J/psi(1S) -> e+ e-) (K*(892)0 -> K+ pi-)]CC",
"Kstar" : "[ [B0]cc -> (J/psi(1S) -> e+ e-) ^(K*(892)0 -> K+ pi-)]CC",
"B" : "[[B0]cc -> (J/psi(1S) -> e+ e-) (K*(892)0 -> K+ pi-)]CC"
}
```
Note that you can constrain more than one intermediate state at once if that fits your decay.

> ## Which constraints to apply {.callout}
> It is important to be aware what assumptions you bake into your ntuple. For example, after you require the vertex constraint you can't use the `IP_CHI2_OWNPV` anymore, since the particle you are looking at is *forced* to point to the PV. Which constraints make most sense for you depends on the questions you want to ask in your analysis.
As seen in the [first-analysis-steps](https://lhcb.github.io/first-analysis-steps/), mass constraints and kinematic fitters can be applied to the decay by using the TupleToolDecayTreeFitter. In this case, some standard information concerning the mother particle and the daughters is added to the desired branch in the decay tree. However, if the daughters are not stable particles and decay further, the daughters of the daughters have no new variables associated to them. In some cases it might be useful to make this information available too. This can done by using the `DecayTreeFitter` via `LoKi functors`, as explained below. This new way of applying constraints has the advantage of flexibility, since the variables one wants to add to the ntuple can be individually specified.

First of all, an instance of the `LoKi::Hybrid::TupleTool` must be created, here called `LoKi_DTF`.

Once you have produced your ntuple you can have a look at the refitted variables.
```shell
root -l DVntuple.root
TupleDstToD0pi_D0ToKpi->cd()
DecayTree->StartViewer()
```
Plotting the raw mass of the D* (without the fit) `Dstar_MM` you should see a broad signal around 2 GeV:

<img src="./img/DstarRaw.png" alt="Dstar raw" style="width: 500px;"/>

Now let us look at the refitted mass of the D*, with the D0 constrained to its nominal mass. It is stored in the variable `Dstar_ConsD_M`. If you plot this you will note that some values are unphysical. So, let's restrict the range we look at to something that makes sense.

On the root prompt use the `arrow-up` key to get the last draw command and modify it to pipe the output into a histogram:
```shell
tv__tree->Draw("Dstar_ConsD_M>>h(200,2000,2030)","","");
LoKi_DTF = LoKi__Hybrid__TupleTool('LoKi_DTF')
```

<img src="./img/DstarRefit.png" alt="Dstar refitted" style="width: 500px;"/>
At this point, the desired variables can be specified in the form of a dictionary. The keys correspond to the names of the new branches that will be added to the ntuple, while the values correspond to the physical variables identified by the keys. Each value has the following syntax:

Note that this plot has 356 entries, although we only have 128 candidates in the raw mass spectrum. The reason for this is, that we typically have several primary vertices per event. When you use the vertex contraint, the fitter is run for each of the possible vertex hypothesis available in the event. So all the `Dstar_ConsD-xxx` variables are in fact arrays, where the first value corresponds to the *best PV* hypothesis. We can plot only those by doing
```shell
tv__tree->Draw("Dstar_ConsD_M[0]>>h(200,2000,2030)","","");
```
and we get the final kinematically refitted Dstar mass:

<img src="./img/DstarRefitBest.png" alt="Dstar refitted best PV" style="width: 500px;"/>

Finally, let's check how the D0 mass constraint has played out.

```shell
tv__tree->Draw("Dstar_ConsD_D0_M[0]>>h(100,1800,1900)","","", 128, 0);
"DTF_FUN(var,bool,particle)"
```

<img src="./img/D0Refit.png" alt="constrained D0" style="width: 500px;"/>
with
* `var` describing the physical variable, according to the same rules that apply to `LoKi functors`;
* `bool` being `True` or `False`, depending on whether the particle must be constrained to originate from the PV;
* `particle` describing the mass constraint.

As expected, the D0 candidates are forced onto their PDG mass value.


> ## Explore {.challenge}
> * Look at the `status` variable to check if the fits converged.
> * Look at the chi2 distribution of the fit
The `CHILD` `LoKi functor` can be used to access the information of the daughters and of the daughters of the daughters. The numbering scheme starts from 1 and corresponds to what specified in the decay descriptor.

For example, in our case:

`DecayTreeFitter` can be told to change some of the hypotheses in the decay tree. This is very useful if you want to slightly change which decays you want to look at. As an example let's say we want to examine the Cabibbo-suppressed decay of the D0 into pi- pi+ instead of K- pi+. For this we add a second fitter, giving it a new name `ConsDpipi`:
```python
dtt.Dstar.addTupleTool('TupleToolDecayTreeFitter/ConsDpipi')
dtt.Dstar.ConsDpipi.constrainToOriginVertex = True
dtt.Dstar.ConsDpipi.Verbose = True
dtt.Dstar.ConsDpipi.daughtersToConstrain = ['D0']
```
We now can tell the fitter to substitute the kaon in the D0 decay by a pion.
```python
dtt.Dstar.ConsDpipi.Substitutions = {
'Charm -> (D0 -> ^K- pi+) Meson': 'pi-',
'Charm -> (D~0 -> ^K+ pi-) Meson': 'pi+'
LoKi_DTF.Variables = {
# B variables
"DTF_B_M" : "DTF_FUN(M,True,'J/psi(1S)')",
"DTF_B_PT" : "DTF_FUN(PT,True,'J/psi(1S)')",
"DTF_B_THETA" : "DTF_FUN(atan(PT/PZ),True,'J/psi(1S)')",
# Jpsi variables
"DTF_Jpsi_M" : "DTF_FUN(CHILD(1,M),True,'J/psi(1S)')",
"DTF_Jpsi_PT" : "DTF_FUN(CHILD(1,PT),True,'J/psi(1S)')",
"DTF_Jpsi_THETA" : "DTF_FUN(CHILD(1,atan(PT/PZ)),True,'J/psi(1S)')",
# Kstar variables
"DTF_Kstar_M" : "DTF_FUN(CHILD(2,M),True,'J/psi(1S)')",
"DTF_Kstar_PT" : "DTF_FUN(CHILD(2,PT),True,'J/psi(1S)')",
"DTF_Kstar_THETA": "DTF_FUN(CHILD(2,atan(PT/PZ)),True,'J/psi(1S)')",
# Lplus variables
"DTF_Lplus_M" : "DTF_FUN(CHILD(1,CHILD(1,M)),True,'J/psi(1S)')",
"DTF_Lplus_PT" : "DTF_FUN(CHILD(1,CHILD(1,PT)),True,'J/psi(1S)')",
"DTF_Lplus_THETA" : "DTF_FUN(CHILD(1,CHILD(1,atan(PT/PZ))),True,'J/psi(1S)')",
# Lminus variables
"DTF_Lminus_M" : "DTF_FUN(CHILD(1,CHILD(2,M)),True,'J/psi(1S)')",
"DTF_Lminus_PT" : "DTF_FUN(CHILD(1,CHILD(2,PT)),True,'J/psi(1S)')",
"DTF_Lminus_THETA" : "DTF_FUN(CHILD(1,CHILD(2,atan(PT/PZ))),True,'J/psi(1S)')",
# K variables
"DTF_K_M" : "DTF_FUN(CHILD(2,CHILD(1,M)),True,'J/psi(1S)')",
"DTF_K_PT" : "DTF_FUN(CHILD(2,CHILD(1,PT)),True,'J/psi(1S)')",
"DTF_K_THETA" : "DTF_FUN(CHILD(2,CHILD(1,atan(PT/PZ))),True,'J/psi(1S)')",
# Pi variables
"DTF_Pi_M" : "DTF_FUN(CHILD(2,CHILD(2,M)),True,'J/psi(1S)')",
"DTF_Pi_PT" : "DTF_FUN(CHILD(2,CHILD(2,PT)),True,'J/psi(1S)')",
"DTF_Pi_THETA" : "DTF_FUN(CHILD(2,CHILD(2,atan(PT/PZ))),True,'J/psi(1S)')",
}
```
In the dictionary that is passed to the `Substitutions` property of the fitter, the keys are decay descriptors, where the respective particle to be substituted is marked with a `^`. The values are the respective new particle hypotheses. The substitution will only work if you start from a decay descriptor that actually matches your candidates. However, you are allowed to generalise parts of the decay. Here we replaced `D*(2010)` with the more general `Charm` and the bachelor `pi-` is just represented by a `Meson`.

Note that the substitution mechanism does not understand the `CC` symbol. Both charge states have to be specified explicitely.

Running the ntuple script again with these additions gives you fit results for the re-interpreted decay.
where we asked to save the constrained mass, transverse momentum, and polar angle of all the particles involved in the decay.

> ## Challenge {.challenge}
> * Compare the outcome of the two fits with the different mass hypothesis
> * Compare the fit quality between the correct and the the wrong hypothesis
> ## Refit {.callout}
>
> Be careful, the refit is performed for each and every variable you have in the list. This will easily make your script rather slow, so think about the variables you really need!

After this, the `LoKi::Hybrid::TupleTool` must be added to one of the branches in our decay tree (generally the mother particle):

```
tuple.ToolList += ['LoKi::Hybrid::TupleTool/LoKi_DTF']
tuple.addTool(LoKi_DTF)
```

The solution to this exercise `ntuple_DTF1.py`, is [available
here](./code/22-decay-tree-fitter/ntuple_DTF1.py).
We can use a similar syntax to constrain the `e+` and `e-` to combine to form a `Psi(2S)` or another particle, if that makes sense for the analysis we are working at.

To conclude, the `LoKi::Hybrid::TupleTool` provides an alternative way to apply constraints with respect to the `TupleToolDecayTreeFitter`. Which of the two approaches is the best depends on the particular situation you have to face. The former is more flexible, but has the disadvantage that it can be very slow if more than a few variables are added due to the fact that the kinematical refit is performed again for each functor.