25.4.9. Full Event Interpretation

Overview

Teaching: 10 min

Exercises: 30 min (including processing)

Prerequisites:

• Lessons 1-2

Questions:

• What are tagged analyses?

• How can I reconstruct a generic B meson?

• How can I use information from both B mesons for my analysis?

Objectives:

• Learn how to use the FEI in your steering file

• Understand the output from the FEI

Introduction

In each ϒ(4S) decay at Belle II, two B mesons are produced. If one of the B mesons decays to a final state involving neutrinos, this B meson cannot be reconstructed completely as the neutrinos are invisible to the detector. The information can, however, be recovered from the second B meson. This is called tagging, the second B meson is referred to as B_tag. The Full Event Interpretation (FEI; pronounced like phi) is an algorithm which reconstructs this B_tag in roughly 10 000 channels. It uses a hierarchical approach to do so, first reconstructing low-level particles such as Kaons and Pions, then intermediate D mesons and finally B mesons. At each stage, the most-likely particle candidates are selected using multivariate classifiers and then used for the next stage. Usually, we do not have to train these classifiers ourselves as this requires large computing resources. Instead, the classifiers are pre-trained and can be loaded from a central database.

Fortunately, for you as a user this is all contained within a few lines of code, allowing you to easily include the FEI in your steering file. In the first part of this exercise, we will use the FEI to reconstruct a B⁺_tag and have a look at its properties. In the second part we will combine the B_tag with a B meson reconstructed in a decay channel involving neutrinos, referred to as B_sig.

Exercise

Find the documentation of the FEI and look up “hadronic tagging” and “semileptonic tagging”. What is the difference between both analysis strategies? What are advantages and disadvantages?

Hint

The documentation can be found in the Advanced Topics section of the analysis module. The definitions can be found in the section titled Hadronic, Semileptonic and Inclusive Tagging.

Solution

The documentation is here: Full event interpretation

Hadronic tagging uses only B meson decays without neutrinos and is therefore more precise (as the exact momentum of the B meson can be known).

Semileptonic tagging uses decays that include one neutrino, these have higher branching fractions but are less precise as the kinematic information is incomplete.

In this exercise we will use hadronic tagging, the procedure for semileptonic tagging is analogous but as it is less commonly used it is not discussed in this exercise.

Reconstructing a B⁺_tag

Prerequisites

As usual, we start our steering file by importing the necessary python packages, creating a basf2 path and loading input data. In addition to the usual python packages (basf2 and modularAnalysis) we also import the fei package.

Exercise

Import basf2, modularAnalysis and fei. In addition, we will be using variable aliases so be sure to import the variable manager. Then create a path and load input data from an mdst file.

You can again use mdst files from /group/belle2/users/tenchini/prerelease-05-00-00a/charged/ on KEKCC.

Solution

#!/usr/bin/env python3

import sys
import basf2 as b2
import fei
import modularAnalysis as ma
from variables import variables as vm

# get input file number from the command line
filenumber = sys.argv[1]

# create path
main = b2.Path()

# load input data from mdst/udst file
ma.inputMdst(
    environmentType="default",
    filename=b2.find_file(f"starterkit/2021/1111540100_eph3_BGx0_{filenumber}.root", "examples"),
    path=main,
)

Now we need the Global Tag in which the weight files for the FEI can be found. This can change once a new central training of the FEI is released, so please check the recommended versions. One way to do so is to use the b2conditionsdb-recommend tool with the mdst file as argument. But there is also a handy method in the modularAnalysis package. Can you find it?

The correct Global Tag must then be used in your steering file by prepending it to the conditions.globaltags list in the basf2 namespace. There is also a convenience function for that!

Exercise

Include the recommended Global Tag in your steering file. For this you need to get the recommended tag using a method found in modularAnalysis and then prepend it to the list using a function documented at Conditions Database Overview.

Hint

You can get the recommended tag using modularAnalysis.getAnalysisGlobaltag

Hint

The function to prepend to the list of global tags is conditions.prepend_globaltag. Now combine this with the last hint!

Solution

Include it in the steering file like this:

Configuring the FEI

Now that all the prerequisites are there, we have to configure the FEI for our purposes. To do this, we need to configure two objects: The fei.get_default_channels function and the fei.FeiConfiguration object.

The fei.get_default_channels function controls which channels the FEI reconstructs. Disabling channels speeds up the FEI so it makes sense to only select what you need by specifying the appropriate arguments. As only charged B mesons are reconstructed in this exercise, the chargedB argument has to be set to True (default) and the neutralB argument to False. The same applies to the hadronic and semileptonic arguments, set hadronic=True and semileptonic=False as you will only reconstruct B mesons in hadronic decay channels. You should also enable B meson decay channels involving baryons with the baryonic=True argument as this increases efficiency.

Exercise

Configure fei.get_default_channels as described above and assign it to a variable.

Solution

    chargedB=True,
    neutralB=False,
    hadronic=True,
    semileptonic=False,
    baryonic=True,
)

The fei.FeiConfiguration class controls the other configuration options of the FEI. Here, the FEI monitoring should be disabled with the appropriate argument (monitor=False) as we are not interested in the internal performance characteristics of the FEI stages. We also have to specify the FEI prefix argument here. This prefix allows distinguishing between different trainings in a single Global Tag and is prefix=FEIv4_2020_MC13_release_04_01_01 for the current central training.

Exercise

Configure fei.FeiConfiguration as described above and assign it to a variable.

Solution

# Set up FEI configuration specifying the FEI prefix
configuration = fei.FeiConfiguration(
    prefix="FEIv4_2021_MC14_release_05_01_12", monitor=False
)

The configuration created above must now be turned into a basf2 path which can be appended to the main path. This is done with the fei.get_path function which takes the channel configuration and the general FEI configuration as arguments and returns a FeiState object. The path attribute of this newly-created FeiState (e.g. feistate.path) is then appended to the main path with the basf2.Path.add_path method.

Exercise

Create the FEI path with fei.get_path() and use its basf2.Path.add_path method to add it to your main path.

Hint

The syntax is mainpath.add_path(feipath).

Solution

feistate = fei.get_path(particles, configuration)

# Add FEI path to the path to be processed
main.add_path(feistate.path)

You have now successfully added the FEI to the main path. The FEI will add a particle list to the datastore called B+:generic. In addition to the usual variables, the B meson candidates in the particle list will also have two extraInfo variables set: The extraInfo(decayModeID) and the extraInfo(SignalProbability) variables. The first one specifies the decay channel in which the B meson has been reconstructed, the second one is the output of the B meson classifier. This value can be used to select B meson candidates to which the FEI assigns a higher probability.

Final Steps

What remains now is adding Monte Carlo matching and creating variable aliases for the extraInfos created by the FEI. Then, you can write the B+:generic list to a .root file. You should already be familiar with these topics from the previous exercises.

Exercise

Add MC matching to the B+:generic particle list and create variable aliases for extraInfo(decayModeID) and extraInfo(SignalProbability).

Then, write out the B+:generic particle list to a root file. Interesting variables are Mbc, deltaE, extraInfo(decayModeID), extraInfo(SignalProbability) and isSignal (or their aliases).

Finally, start the event loop with a call to basf2.process.

Hint

Go back to First steering file to see the step-by-step instructions on how to add MC matching, write the ntuple and start the event loop.

Aliases are introduced in Various additions, check there if you are unsure how to use them.

Solution

# Rank B+ candidates by signal classifier output
vm.addAlias("decayModeID", "extraInfo(decayModeID)")

# Store tag-side variables of interest.
ma.variablesToNtuple(
    "B+:generic",
    [
        "Mbc",
        "deltaE",
        "mcErrors",
        "SigProb",
        "decayModeID",
        "FEIProbRank",
    ],
    filename="B_charged_hadronic.root",
    path=main,
)

The FEI returns not only one B meson candidate for each event but up to 20. Using the modularAnalysis.rankByHighest function, it is possible to rank the candidates by the B meson classifier output in the extraInfo(SignalProbability) variable. This is optional but often useful to select the best, i.e. most likely correct, candidate.

Exercise (optional)

Use rankByHighest to rank the B meson candidates in the B+:generic list by the extraInfo(SignalProbability) variable. Write the rank to a new variable called FEIProbabilityRank.

Don’t forget to create an alias for this variable (within an extraInfo metavariable) and write this to the nTuple.

Hint

You should already be familiar with Best Candidate Selection from the Various additions lesson. The documentation on rankByHighest can be found here: modularAnalysis.rankByHighest.

If the outputVariable argument of modularAnalysis.rankByHighest is called FEIProbabilityRank, the alias should be created for extraInfo(FEIProbabilityRank).

Solution

ma.rankByHighest(
    particleList="B+:generic",
    variable="extraInfo(SignalProbability)",
    outputVariable="FEIProbabilityRank",
    path=main,
)
vm.addAlias("FEIProbRank", "extraInfo(FEIProbabilityRank)")

vm.addAlias("SigProb", "extraInfo(SignalProbability)")
vm.addAlias("decayModeID", "extraInfo(decayModeID)")

# Store tag-side variables of interest.
ma.variablesToNtuple(
    "B+:generic",
    [
        "Mbc",
        "deltaE",
        "mcErrors",
        "SigProb",
        "decayModeID",
        "FEIProbRank",
        "isSignal",
    ],
    filename="B_charged_hadronic.root",
    path=main,
)

You can now execute your steering file which should look somewhat like this:

Final steering file

#!/usr/bin/env python3

import sys
import basf2 as b2
import fei
import modularAnalysis as ma
from variables import variables as vm

# get input file number from the command line
filenumber = sys.argv[1]

# create path
main = b2.Path()

# load input data from mdst/udst file
ma.inputMdst(
    environmentType="default",
    filename=b2.find_file(f"starterkit/2021/1111540100_eph3_BGx0_{filenumber}.root", "examples"),
    path=main,
)

# Add the database with the classifier weight files for the FEI
b2.conditions.prepend_globaltag(ma.getAnalysisGlobaltag())

# Get FEI default channels.
# Utilise the arguments to toggle on and off certain channels
particles = fei.get_default_channels(
    chargedB=True,
    neutralB=False,
    hadronic=True,
    semileptonic=False,
    baryonic=True,
)

# Set up FEI configuration specifying the FEI prefix
configuration = fei.FeiConfiguration(
    prefix="FEIv4_2021_MC14_release_05_01_12", monitor=False
)

# Get FEI path
feistate = fei.get_path(particles, configuration)

# Add FEI path to the path to be processed
main.add_path(feistate.path)

# Add MC matching when applying to MC.
# This is required for variables like isSignal and mcErrors below
ma.matchMCTruth(list_name="B+:generic", path=main)

# Rank B+ candidates by signal classifier output
ma.rankByHighest(
    particleList="B+:generic",
    variable="extraInfo(SignalProbability)",
    outputVariable="FEIProbabilityRank",
    path=main,
)
vm.addAlias("FEIProbRank", "extraInfo(FEIProbabilityRank)")

vm.addAlias("SigProb", "extraInfo(SignalProbability)")
vm.addAlias("decayModeID", "extraInfo(decayModeID)")

# Store tag-side variables of interest.
ma.variablesToNtuple(
    "B+:generic",
    [
        "Mbc",
        "deltaE",
        "mcErrors",
        "SigProb",
        "decayModeID",
        "FEIProbRank",
        "isSignal",
    ],
    filename="B_charged_hadronic.root",
    path=main,
)

# Process events
b2.process(main)
print(b2.statistics)

Offline analysis

Now that you have created your ntuple, we can have a look at the properties of the B mesons we have created.

You have already looked at the beam-constrained mass Mbc in First steering file. For correctly reconstructed B mesons, this variable should peak at the B meson mass. It is therefore a good indicator for the quality of the B mesons we have reconstructed.

Exercise

Load your ntuple file into python, either using root_pandas or uproot. Then, plot the distribution of Mbc from 5.15 – 5.3 GeV.

You should see broad peak with a sharp drop-off below 5.2 GeV. This drop-off is caused by a fixed pre-cut in the FEI. Candidates below this threshold are rejected before classification as they are almost always incorrect.

Solution

# Include this only if running in a Jupyter notebook
%matplotlib inline

import matplotlib.pyplot as plt
from root_pandas import read_root

df = read_root('B_charged_hadronic.root')

fig, ax = plt.subplots()
n, bins, patches = ax.hist(df['Mbc'], bins=30, range=(5.15, 5.3))
ax.set_xlabel(r'$\mathrm{M}_{\mathrm{bc}}$ in GeV/c^2')
ax.set_ylabel('Number of candidates')
fig.savefig('m_bc.pdf')

Question

The distribution of Mbc you have just plotted doesn’t peak at the B meson mass of 5.28 GeV. Can you explain this?

Solution

We haven’t really used the classifier output of the FEI yet. The up to 20 candidates in each event are selected by FEI Signal Probability but many still have low absolute classifier values and by definition almost all of them are misreconstructed.

Exercise

Increase the signal purity by requiring minimum values of extraInfo__boSignalProbability__bc. Try cuts at 0.001, 0.01 and 0.1!

You can also try selecting the candidate with the highest FEI Probability in each event by using the extraInfo__boFEIProbabilityRank__bc variable created with modularAnalysis.rankByHighest (if you have done this).

These selections should increase the purity of the B meson candidates under consideration and lead to a sharper peak at the B mass.

You can also have a look at the correctly reconstructed B mesons by requiring isSignal == 1. By comparing this to the cuts on the FEI classifier you can see how well the FEI identifies correctly reconstructed B mesons.

Solution

# Include this only if running in a Jupyter notebook
%matplotlib inline

import matplotlib.pyplot as plt
from root_pandas import read_root

df = read_root('B_charged_hadronic.root')

fig, ax = plt.subplots()

# If you didn't do the optional exercise, the 'FEIProbRank' column won't be there.
# Just remove this part of the query in this case.
sig_prob_cut_df = df.query('(SigProb > 0.01)
                            & (FEIProbRank == 1)')

n, bins, patches = ax.hist(sig_prob_cut_df['Mbc'], bins=30, range=(5.15, 5.3))
ax.set_xlabel(r'$\mathrm{M}_{\mathrm{bc}}$ in GeV/c^2')
ax.set_ylabel('Total number of candidates')
ax.set_title('SigProb > 0.01')
fig.savefig('m_bc_cut_0_01.pdf')

Congratulations, you have now discovered the B meson in Monte Carlo data! This concludes the first part of this lesson. The second part of this lesson will show you now how to use the reconstructed B_tag in your own analysis.

Reconstructing the full ϒ(4S) event

In this part of the lesson, we will reconstruct the full ϒ(4S) event using the FEI. The B⁰ decay process we will be looking at is the semileptonic decay to a charged pion, a muon and a muon neutrino.

The FEI skim

You might have noticed that applying the FEI takes some time, even for the small file we have just processed. For this reason and to save computing resources, datasets with pre-applied FEI tagging exist. These preprocessed datasets are called skims.

We will be using a FEI-skimmed file in this exercise in which the B0:generic particle list already exists. If you would like to know more about skimming, you can have a look into Skimming.

The input file we will be using, found at /home/belle2/mbauer/fei_tutorial/fei_skimmed_xulnu.udst.root, only contains decays of B0 mesons to a light lepton and a charged pi or rho meson. This way we don’t have to process as much data as we would have to for a file containing decays in all B decay channels.

Prerequisites (part 2)

Lets get started with the usual steps. Nothing here should be new to you.

Exercise

Start a new steering file. In this file, you won’t need the fei package so you can skip this import. Create a path and load the udst file /home/belle2/mbauer/fei_tutorial/fei_skimmed_xulnu.udst.root.

NOTE: You can still use modularAnalysis.inputMdst to do this, even though it’s a udst file.

Then, fill two particle lists with muons and charged pions. For the muons, you can require a muonID above 0.9, for the pions a pionID above 0.5. For both you should apply some requirements on the track, you can use dr < 0.5 and abs(dz) < 2 and nCDCHits > 20 and thetaInCDCAcceptance

Hint

You should already be familiar with this from First steering file.

Solution

#!/usr/bin/env python3

import basf2 as b2
import modularAnalysis as ma
from variables import variables as vm

main = b2.Path()

ma.inputMdst(
    "default",
    b2.find_file("starterkit/2021/fei_skimmed_xulnu.udst.root", "examples"),
    path=main,
)

good_track = (
    "dr < 0.5 and abs(dz) < 2 and nCDCHits > 20 and thetaInCDCAcceptance"
)

ma.fillParticleList("mu-", "muonID > 0.9 and " + good_track, path=main)
ma.fillParticleList("pi-", "pionID > 0.5 and " + good_track, path=main)

ϒ(4S) Reconstruction

Now, let’s get started with the reconstruction. We will first create the signal B meson, then combine it with the B meson provided by the FEI to get the ϒ(4S).

As the semileptonic decay we are analysing contains a neutrino, a few things have to be considered. The neutrino is not used the reconstruction of the B_sig as it can’t be seen in the Belle II detector. This leads to a discrepancy between the simulated decay and the reconstructed decay, causing the isSignal variable to always be zero.

To solve this you can tell the MC matching algorithm to accept missing neutrinos by using the ?nu flag, similiar to the ?addbrems flag introduced in Various additions. Just add this flag to the end of the decay string in modularAnalysis.reconstructDecay. 1

Exercise

Reconstruct a B⁰ meson particle list called B0:signal from a positive pion and a muon. Add the ?nu flag as described above.

Solution

ma.reconstructDecay("B0:signal -> pi- mu+ ?nu", cut="", path=main)

So far, we have not used the FEI. Now, we will use the B_tag in the udst file and combine it with the B_sig we have just created.

Exercise

Combine the B0:generic from the udst file with the anti-B0:signal to a list called Upsilon(4S):opposite_cp .

NOTE: The B0:generic should come first in the decay string, otherwise the missing mass squared variable we’re using later won’t know which of the B mesons is the tag and which is the signal.

Can you think of a reason for the opposite_cp identifier of the Upsilon(4S) list?

Have we forgotten something?

Hint

Think of special properties of B⁰ mesons compared to B⁺ mesons.

For the implementation: You will most likely need one more modularAnalysis.reconstructDecay to create the second particle list and the modularAnalysis.copyLists function to combine both lists.

Solution

To account for B⁰ meson mixing, you should also combine same-sign B⁰ mesons as the anti-B⁰ can oscillate to a B⁰.

ma.reconstructDecay(
    "Upsilon(4S):opposite_cp -> B0:generic anti-B0:signal", cut="", path=main
)
ma.reconstructDecay(
    decayString="Upsilon(4S):same_cp -> B0:generic B0:signal",
    cut="",
    path=main,
)
# Combine the two Upsilon(4S) lists to one. Note: Duplicates are removed.
ma.copyLists(
    outputListName="Upsilon(4S)",
    inputListNames=["Upsilon(4S):opposite_cp", "Upsilon(4S):same_cp"],
    path=main,
)

Now that we have reconstructed the full ϒ(4S), we will create a Rest of Event. You have already done this in The Rest of Event (ROE) for a B meson, here however we will create a Rest of Event for the ϒ(4S). This allows us to count the tracks left over after reconstructing the ϒ(4S), for two correctly reconstructed B mesons there should be no tracks left over.

For this to work we have to use a slightly different ROE mask than in the ROE chapter. To match the cuts used by the FEI to reconstruct B_tag candidates, we have to tighten the cut on dr to below 2 and the cut on the absolute value of dz (abs(dz)) to below 4. The two other cuts (on pt and thetaInCDCAcceptance) can be left as-is.

Exercise

Create a Rest of Event for the Upsilon(4S) list. Then, append an ROE mask using the cuts introduced in the chapter The Rest of Event (ROE) and the additional track cuts mentioned above.

Hint

The track cuts should be thetaInCDCAcceptance and pt > 0.075 and dr < 2 and abs(dz) < 4 and the ecl cuts thetaInCDCAcceptance and E > 0.05. To create the ROE, use modularAnalysis.buildRestOfEvent. To append the ROE mask, use modularAnalysis.appendROEMasks.

Solution

ma.buildRestOfEvent("Upsilon(4S)", path=main)
track_based_cuts = "thetaInCDCAcceptance and pt > 0.075 and dr < 2 and abs(dz) < 4"
ecl_based_cuts = "thetaInCDCAcceptance and E > 0.05"
roe_mask = ("my_mask", track_based_cuts, ecl_based_cuts)
ma.appendROEMasks("Upsilon(4S)", [roe_mask], path=main)

Writing out the nTuple

We are now ready to add MC matching and write out the nTuple.

While the MC matching is applied to the Upsilon(4S) list, you can also access the daughter’s truth variables using the daughter metavariable.

This is especially useful in FEI analyses:

As you might have noticed in the first part of the exercise, the number of perfectly (i.e. isSignal == 1.0) reconstructed B_tag mesons is not very large As we are only really interested in the B_sig, the isSignal variable of this B meson can be a better signal definition than the isSignal variable of the Upsilon(4S)

In addition to the properties of the B mesons, we can now also use information from the full event.

An example here is the missing mass squared in the variable m2RecoilSignalSide. This quantity should peak at zero for decays in which only one neutrino is missing and thus provides high separating power in (semi-)leptonic analyses. There are different implementations of the missing mass squared in basf2, this version uses the explicit B_tag momentum (here reconstructed by the FEI) and has therefore a high resolution.

A second variable with high separating power is the number of additional charged tracks introduced above. This variable is called nROE_Charged in basf2. It needs the ROE mask name defined above as an argument in brackets.

Exercise

Add MC matching for the Upsilon(4S) to your path.

Write the Upsilon(4S) particle list to an nTuple. Include the variables used in the first part of the lesson, the signal variables mentioned above (both for the B_sig and the ϒ(4S)), the missing mass squared and the number of additional tracks.

It is recommended to define aliases for the variables as otherwise the variables will be long and unwieldy in your offline analysis.

Finally, don’t forget to process the path!

Hint

If you have forgotten how to use aliases: This is introduced in Various additions.

If you have adapted the solution given here, the variable containing the number of additional tracks will be nROE_Charged(my_mask).

Solution

ma.matchMCTruth(list_name="Upsilon(4S)", path=main)

vm.addAlias("Btag_SigProb", "daughter(0, extraInfo(SignalProbability))")
vm.addAlias("Btag_decayModeID", "daughter(0, extraInfo(decayModeID))")
vm.addAlias("Btag_Mbc", "daughter(0, Mbc)")

vm.addAlias("Bsig_isSignal",
            "daughter(1, isSignal)")
vm.addAlias("nCharged", "nROE_Charged(my_mask)")

ma.variablesToNtuple(
    "Upsilon(4S)",
    variables=[
        "Btag_SigProb",
        "Btag_decayModeID",
        "Btag_Mbc",
        "Bsig_isSignal",
        "isSignal",
        "m2RecoilSignalSide",
        "nCharged"
    ],
    filename='Upsilon4S.root',
    path=main,
)
# Process events
b2.process(main)

You can now execute your steering file which should look somewhat like this:

Final steering file

#!/usr/bin/env python3

import basf2 as b2
import modularAnalysis as ma
from variables import variables as vm

main = b2.Path()

ma.inputMdst(
    "default",
    b2.find_file("starterkit/2021/fei_skimmed_xulnu.udst.root", "examples"),
    path=main,
)

good_track = (
    "dr < 0.5 and abs(dz) < 2 and nCDCHits > 20 and thetaInCDCAcceptance"
)

ma.fillParticleList("mu-", "muonID > 0.9 and " + good_track, path=main)
ma.fillParticleList("pi-", "pionID > 0.5 and " + good_track, path=main)

ma.reconstructDecay("B0:signal -> pi- mu+ ?nu", cut="", path=main)

ma.reconstructDecay(
    "Upsilon(4S):opposite_cp -> B0:generic anti-B0:signal", cut="", path=main
)
ma.reconstructDecay(
    decayString="Upsilon(4S):same_cp -> B0:generic B0:signal",
    cut="",
    path=main,
)
# Combine the two Upsilon(4S) lists to one. Note: Duplicates are removed.
ma.copyLists(
    outputListName="Upsilon(4S)",
    inputListNames=["Upsilon(4S):opposite_cp", "Upsilon(4S):same_cp"],
    path=main,
)

ma.buildRestOfEvent("Upsilon(4S)", path=main)
track_based_cuts = "thetaInCDCAcceptance and pt > 0.075 and dr < 2 and abs(dz) < 4"
ecl_based_cuts = "thetaInCDCAcceptance and E > 0.05"
roe_mask = ("my_mask", track_based_cuts, ecl_based_cuts)
ma.appendROEMasks("Upsilon(4S)", [roe_mask], path=main)

ma.matchMCTruth(list_name="Upsilon(4S)", path=main)

vm.addAlias("Btag_SigProb", "daughter(0, extraInfo(SignalProbability))")
vm.addAlias("Btag_decayModeID", "daughter(0, extraInfo(decayModeID))")
vm.addAlias("Btag_Mbc", "daughter(0, Mbc)")

vm.addAlias("Bsig_isSignal",
            "daughter(1, isSignal)")
vm.addAlias("nCharged", "nROE_Charged(my_mask)")

ma.variablesToNtuple(
    "Upsilon(4S)",
    variables=[
        "Btag_SigProb",
        "Btag_decayModeID",
        "Btag_Mbc",
        "Bsig_isSignal",
        "isSignal",
        "m2RecoilSignalSide",
        "nCharged"
    ],
    filename='Upsilon4S.root',
    path=main,
)
# Process events
b2.process(main)
print(b2.statistics)

Offline analysis (part 2)

Like in the first part of this lesson, you can now analyse your nTuple. As before, you can use the FEI signal probability (now under the alias Btag_SigProb if you have adapted the example) to select more pure ϒ(4S) candidates and plot m2RecoilSignalSide for different values of the classifier.

You can also see how cuts on the number of additional tracks change the histogram.

NOTE: A histogram of Mbc will look quite different in this part of the exercise, this is because in the last exercise we have used a generic MC sample and in this exercise we are using an MC sample with only four decay channels.

Exercise

Plot m2RecoilSignalSide for Bsig_isSignal == 1.0 to see what well-reconstructed B -> pi l nu decays should look like. Then, plot m2RecoilSignalSide for different cuts on the FEI classifier output. See also how the shape of m2RecoilSignalSide changes when requiring zero additional tracks.

Solution

# Include this only if running in a Jupyter notebook
%matplotlib inline

import matplotlib.pyplot as plt
from root_pandas import read_root

df = read_root('Upsilon4S.root')

fig, ax = plt.subplots()

# Requiring only the B_sig to be correctly reconstructed.
# Try what happens if you require the whole Y(4S) to be correct!
signal_df = df.query('(Bsig_isSignal  == 1.0)')

n, bins, patches = ax.hist(signal_df['m2RecoilSignalSide'], bins=30, range=(-2, 5))
ax.set_xlabel(r'$\mathrm{M}_{\mathrm{bc}}$ in GeV/c^2')
ax.set_ylabel('Total number of candidates')
ax.set_title('Bsig_isSignal == 1')
fig.savefig('m2RSS_Signal.pdf')

fig, ax = plt.subplots()

# This is just an example cut, you can try without the cut on nCharged
# and for different values on BtagSigProb
cut_df = df.query('(Btag_SigProb > 0.01) & (nCharged == 0.0)')

n, bins, patches = ax.hist(cut_df['m2RecoilSignalSide'], bins=30, range=(-2, 5))
ax.set_xlabel(r'$\mathrm{M}_{\mathrm{bc}}$ in GeV/c^2')
ax.set_ylabel('Total number of candidates')
ax.set_title('SigProb > 0.01')
fig.savefig('m2RSS_FEIcut_0_01_nCharged_0.pdf')

Congratulations, you now know how to run the FEI and how to use it in your analysis. If you would like to know more you can always read the extensive documentation of the FEI. Here you can also find instructions on how to train the FEI and explanations on the code structure.

Key points

• Get the weight files from the conditions database

• Add the FEI to your path with fei.get_default_channels and fei.FeiConfiguration.

• FEI Purity and efficiency are controlled by a cut on extraInfo(SignalProbability)

• The B_tag from the FEI can be used to construct a complete ϒ(4S) event.

Stuck? We can help!

If you get stuck or have any questions to the online book material, the #starterkit-workshop channel in our chat is full of nice people who will provide fast help.

Refer to Collaborative Tools. for other places to get help if you have specific or detailed questions about your own analysis.

Improving things!

If you know how to do it, we recommend you to report bugs and other requests with JIRA. Make sure to use the documentation-training component of the Belle II Software project.

If you just want to give very quick feedback, use the last box “Quick feedback”.

Please make sure to be as precise as possible to make it easier for us to fix things! So for example:

• typos (where?)

• missing bits of information (what?)

• bugs (what did you do? what goes wrong?)

• too hard exercises (which one?)

• etc.

If you are familiar with git and want to create your first pull request for the software, take a look at How to contribute. We’d be happy to have you on the team!

Quick feedback!

Author of this lesson

Moritz Bauer

Footnotes

1

Alternatively, you can also use the isSignalAcceptMissingNeutrino variable instead of isSignal but we will be going with the ?nu flag and isSignal in our example.