3.2 Beyond the Basics: Tracking and Vertexing

Original authors: Emmanuel Francois Perez, Clement Helsens

Getting the FCCAnalyses¶

The FCCAnalyses framework is provided already compiled in the Key4hep stack, however once the stack is released this version of the FCCAnalyses can’t be updated anymore. If one wants to benefit from the changes in the upstream FCCAnalyses, they need to compile the latest version of it. The compilation process is quite streamlined and requires only few steps.

As a first step create a fork of the FCCAnalyses project on GitHub. More details about FCC Software development workflow can be found in Development workflow.

After a short while the forking should be done and you can download the FCCAnalyses to your machine. Go inside the area that you have setup for this tutorial (for example mkdir vtx-tutorial) and clone your FCCAnalyses repository:

cd vtx-tutorial
git clone --branch pre-edm4hep1 https://github.com/<your-github-handle>/FCCAnalyses.git

Go inside the FCCAnalyses directory and run the building of the FCCAnalyses with:

cd FCCAnalyses
source ./setup.sh
fccanalysis build

After the successful building of the FCCAnalyses you can return back to the main directory:

cd ..

To check if you are using the locally build version of FCCAnalyses use

which fccanalysis

It should lead to executable in your local install directory.

Adding analyzers to the FCCAnalyses¶

In order to extend the number of the available FCCAnalyses analyzers (C++ functions or functors used in the .Define() and .Filter methods of the ROOT RDataFrame) we can either provide them in the additional C++ header file(s) or edit the source code of the FCCAnalyses analyzers library.

For a quick addition of a few analyzers the first method is recommended, however since in this method the additional C++ header file(s) need to be JIT compiled every time the analysis is run the start-up time of the analysis can become very long.

To register the C++ header file with the additional analyzers in your analysis use includePaths attribute of the analysis script:

includePaths = ["analyzers.h"]

The second method involves directly editing the source code located in the analyzers/dataframe directory of FCCAnalyses and afterwards running fccanalysis build command.

Exercise¶

Let’s add a dummy analyzer which extracts $x$ component of the particle momentum from the collection of reconstructed particles.

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/dummy_analysis.py

vim analyzers.h

#ifndef FCCANA_ADDITIONAL_ANALYZERS_H
#define FCCANA_ADDITIONAL_ANALYZERS_H

#include "ROOT/RVec.hxx"
#include "edm4hep/ReconstructedParticleData.h"

namespace FCCAnalyses :: VtxAna {
  namespace rv = ROOT::VecOps;

  rv::RVec<float> dummy_analyzer(const rv::RVec<edm4hep::ReconstructedParticleData>& parts) {
    rv::RVec<float> output;
    for (size_t i = 0; i < parts.size(); ++i)
      output.emplace_back(parts.at(i).momentum.x);
    return output;
  }
}

#endif /* FCCANA_ADDITIONAL_ANALYZERS_H */

fccanalysis run dummy_analysis.py --test --nevents 10 --output dummy_result.root

Reconstruction of the primary vertex and of primary tracks¶

Let’s start by running primary vertex reconstruction on a few events of one test file:

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/analysis_primary_vertex.py
fccanalysis run analysis_primary_vertex.py --test --nevents 1000 --output primary_Zuds.root

Note: with the option --test, we process the file that is hard-coded in the attribute testFile inside of the analysis_primary_vertex.py script. In this case, it is a file of $Z \rightarrow q \bar{q}$ with $q=u,d,s$.

The resulting ntuple primary_Zuds.root contains the MC event vertex MC_PrimaryVertex, and the reconstructed primary vertex PrimaryVertex.

    # MC event primary vertex
    .Define("MC_PrimaryVertex",
            "FCCAnalyses::MCParticle::get_EventPrimaryVertex(21)(Particle)")

    # Number of tracks
    .Define("ntracks",
            "ReconstructedParticle2Track::getTK_n(EFlowTrack_1)")

    # Fit all tracks of the events to a common vertex --- here using a
    # beam-spot constraint:

    # VertexObject_allTracks is an object of type `VertexingUtils::FCCAnalysesVertex`
    # It contains in particular:
    #   - an `edm4hep::VertexData`:
    #     - std::int32_t primary{}; ///< boolean flag, if vertex is the primary vertex of the event
    #     - float chi2{}; ///< chi-squared of the vertex fit
    #     - edm4hep::Vector3f position{}; ///< [mm] position of the vertex.
    #     - std::array<float, 6> covMatrix{}; ///< covariance matrix of the position (stored as lower triangle matrix, i.e. cov(xx),cov(y,x),cov(z,x),cov(y,y),...)
    #   - `ROOT::VecOps::RVec<float> reco_chi2`: the contribution to
    #     the chi2 of all tracks used in the fit
    #   - `ROOT::VecOps::RVec<TVector3> updated_track_momentum_at_vertex`:
    #     the post-fit (px, py, pz ) of the tracks, at the vertex
    #     (and not at their d.c.a.)
    .Define("VertexObject_allTracks",
            "VertexFitterSimple::VertexFitter_Tk(1, EFlowTrack_1, true, 4.5, 20e-3, 300)")

    # `EFlowTrack_1` is the collection of all tracks (the fitting
    # method can of course be applied to a subset of tracks
    # (see later)).
    # "true" means that a beam-spot constraint is applied. Default is
    # no BSC. Following args are the BS size and position, in mum:
    #   - bool BeamSpotConstraint = false,
    #   - double sigmax=0., double sigmay=0., double sigmaz=0.,
    #   - double bsc_x=0., double bsc_y=0., double bsc_z=0. );

    # This returns the `edm4hep::VertexData`:
    .Define("Vertex_allTracks",
            "VertexingUtils::get_VertexData(VertexObject_allTracks)")  # primary vertex, in mm

    # This is not a good estimate of the primary vertex: even in a Z -> uds event, there are displaced tracks (e.g. Ks, Lambdas), which would bias the fit.
    # Below, we determine the "primary tracks" using an iterative algorithm - cf LCFI+.
    .Define("RecoedPrimaryTracks",
            "VertexFitterSimple::get_PrimaryTracks(EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0.)")

    # Now we run again the vertex fit, but only on the primary tracks :
    .Define("PrimaryVertexObject",
            "VertexFitterSimple::VertexFitter_Tk(1, RecoedPrimaryTracks, true, 4.5, 20e-3, 300)")
    .Define("PrimaryVertex",
            "VertexingUtils::get_VertexData(PrimaryVertexObject)")

    # It is often useful to retrieve the secondary (i.e. non-primary) tracks, for example to search for secondary vertices.
    # The method below simply "subtracts" the primary tracks from the full collection :
    .Define("SecondaryTracks",
            "VertexFitterSimple::get_NonPrimaryTracks(EFlowTrack_1, RecoedPrimaryTracks)")

To produce example plots, run the ROOT plotting macro plot_primary_vertex.C.

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/plot_primary_vertex.C

root -b -q 'plot_primary_vertex.C()'

This produces normalized $\chi^2$ of the primary vertex fit, the resolutions in $x$, $y$, $z$, and the pulls of the fitted vertex position.

Exercise 1¶

Add the number of primary and secondary tracks into the ntuple using the function ReconstructedParticle2Track::getTK_n(ROOT::VecOps::RVec<edm4hep::TrackState> x), see the definition here.

    # Number of primary and secondary tracks:
    .Define("n_RecoedPrimaryTracks",
            "ReconstructedParticle2Track::getTK_n(RecoedPrimaryTracks)")
    .Define("n_SecondaryTracks",
            "ReconstructedParticle2Track::getTK_n(SecondaryTracks)")
    # equivalent: (this is to show that a simple C++ statement can be
    # included in a ".Define")
    .Define("n_SecondaryTracks_v2",
            "return ntracks - n_RecoedPrimaryTracks;")

    branchList = [
        "MC_PrimaryVertex",
        "ntracks",
        "Vertex_allTracks",
        "PrimaryVertex",
        "n_RecoedPrimaryTracks",
        "n_SecondaryTracks",
        "n_SecondaryTracks_v2",
    ]

Exercise 2¶

Add the total $p_{T}$ that is carried by the primary tracks. This requires some simple analysis code to be written (added to our analyzers.h file). Then, the analysis script needs to be updated to include includePaths = ["analyzers.h"] attribute.

...
#include <cmath>
#include "FCCAnalyses/VertexingUtils.h"
...

  double sum_momentum_tracks(const VertexingUtils::FCCAnalysesVertex& vertex) {
    ROOT::VecOps::RVec<TVector3> momenta = vertex.updated_track_momentum_at_vertex;

    double sum = 0;
    for (const auto& p: momenta) {
       double px = p[0];
       double py = p[1];
       double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
       sum += pt;
    }

    return sum;
  }

includePaths = ["analyzers.h"]

            ...
            # Total pT carried by the primary tracks:
            .Define("sum_pt_primaries",
                    "VtxAna::sum_momentum_tracks( PrimaryVertexObject )")
            ...

        ...
        branchList = [
            ...
            "sum_pt_primaries",
        ]
        ...

Exercise 3¶

Compare these distributions in $Z \rightarrow uds$ events and in $Z \rightarrow b\bar{b}$ events.

Exercise 4¶

Exercise to go beyond

The reconstruction of all secondary vertices following the LCFI+ algorithm has been implemented in FCCAnalyses by Kunal Gautam and Armin Ilg in the pull request PR#206. It contains example analysis script examples/FCCee/vertex_lcfiplus/analysis_SV.py which:

• reconstructs the primary vertex and primary tracks as done above

• reconstructs jets using the Durham algorithm

• reconstructs secondary vertices within all jets, and determines some properties of these secondary vertices It is also possible to reconstruct all secondary vertices in an event, without reconstructing jets.

Reconstruction of displaced vertices in an exclusive decay chain¶

We consider here $Z \rightarrow b \bar{b}$ events. When a $B_s$ is produced, it is forced to decay into $J/\Psi \Phi$ with $J/\Psi \rightarrow \mu^+\mu^-$ and $\Phi \rightarrow K^+ K^-$. We want to reconstruct the $B_s$ decay vertex and determine the resolution on the position of this vertex. Here, we use the MC-matching information to figure out which are the reconstructed tracks that are matched to the $B_s$ decay products, and we fit these tracks to a common vertex. That means, we “seed” the vertex reconstruction using the MC-truth information. Let’s run the following:

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/analysis_Bs2JpsiPhi_MCseeded.py
fccanalysis run analysis_Bs2JpsiPhi_MCseeded.py --test --nevents 1000 --output Bs2JpsiPhi_MCseeded.root

The ntuple Bs2JpsiPhi_MCseeded.root contains the MC decay vertex of the $B_s$, and the reconstructed decay vertex.

    # MC indices of the decay
    # Bs (PDG = 531) -> mu+ (PDG = -13) mu- (PDG = 13) K+ (PDG = 321) K- (PDG = -321)
    # Retrieves a vector of integers which correspond to indices in the
    # Particle block
    # vector[0] = the mother, and then the daughters in the order
    # specified, i.e. here [1] = the mu+, [2] = the mu-, [3] = the K+,
    # [4] = the K-
    #
    # Boolean arguments:
    #   1st: `stableDaughters`, when set to true, the daughters specified
    #        in the list are looked for among the final, stable
    #        particles that come out from the mother, i.e. the decay
    #        tree is explored recursively if needed.
    #   2nd: `chargeConjugateMother`
    #   3rd: `chargeConjugateDaughters`
    #   4th: `inclusiveDecay`, when set to false, if a mother is found,
    #        that decays into the particles specified in the list plus
    #        other particle(s), this decay is not selected.
    # If the event contains more than one such decays, only the first
    # one is kept.
    .Define("Bs2MuMuKK_indices",
            "MCParticle::get_indices(531, {-13,13,321,-321}, true, true, true, false) (Particle, Particle1)")

    # select events for which the requested decay chain has been found:
    .Filter("Bs2MuMuKK_indices.size() > 0")

    # the mu+ (MCParticle) that comes from the Bs decay :
    .Define("MC_Muplus", "return Particle.at(Bs2MuMuKK_indices[1]);")
    # Decay vertex (an `edm4hep::Vector3d`) of the Bs (MC) = production
    # vertex of the muplus:
    .Define("BsMCDecayVertex", "return MC_Muplus.vertex;")

    # Returns the RecoParticles associated with the four Bs decay products.
    # The size of this collection is always 4 provided that
    # Bs2MuMuKK_indices is not empty, possibly including "dummy"
    # particles in case one of the legs did not make a RecoParticle
    # (e.g. because it is outside the tracker acceptance). This is done
    # on purpose, in order to maintain the mapping with the indices ---
    # i.e. the 1st particle in the list BsRecoParticles is the mu+,
    # then the mu-, etc.
    # (selRP_matched_to_list ignores the unstable MC particles that are
    # in the input list of indices hence the mother particle, which is
    # the [0] element of the Bs2MuMuKK_indices vector).
    #
    # The matching between RecoParticles and MCParticles requires 4
    # collections. For more detail, see
    # https://github.com/HEP-FCC/FCCAnalyses/tree/master/examples/basics
    .Define("BsRecoParticles",
            "ReconstructedParticle2MC::selRP_matched_to_list(Bs2MuMuKK_indices, MCRecoAssociations0, MCRecoAssociations1, ReconstructedParticles, Particle)")

    # the corresponding tracks --- here, dummy particles, if any, are
    # removed, i.e. one may have < 4 tracks, e.g. if one muon or kaon
    # was emitted outside of the acceptance
    .Define("BsTracks",
            "ReconstructedParticle2Track::getRP2TRK(BsRecoParticles, EFlowTrack_1)")

    # number of tracks in this BsTracks collection (= the #tracks
    # used to reconstruct the Bs vertex)
    .Define("n_BsTracks",
            "ReconstructedParticle2Track::getTK_n(BsTracks)")

    # Fit the tracks to a common vertex. That would be a secondary
    # vertex, hence we put a "2" as the first argument of
    # VertexFitter_Tk: First the full object, of type
    # Vertexing::FCCAnalysesVertex
    .Define("BsVertexObject", "VertexFitterSimple::VertexFitter_Tk(2, BsTracks)")
    # from which we extract the edm4hep::VertexData object, which
    # contains the vertex position in mm
    .Define("BsVertex",
            "VertexingUtils::get_VertexData(BsVertexObject)")

When you run the root macro plot_Bs2JsiPhi.C it will produce various plots showing the vertex $\chi^2$, the vertex resolutions and the pulls of the vertex fit

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/plots_Bs2JsiPhi.C
root -b -q "plot_Bs2JsiPhi.C()"

Analysis of $\tau \rightarrow 3 \mu$¶

The analysis showcased in Reconstruction of displaced vertices in an exclusive decay chain is used as a stepping stone for the following set of exercises.

Exercise 1¶

Start from the analysis_Bs2JpsiPhi_MCseeded.py analysis script and adapt it to the decay of $\tau \rightarrow 3 \mu$.

cp analysis_Bs2JpsiPhi_MCseeded.py analysis_Tau3Mu_MCseeded.py

# or

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/analysis_Bs2JpsiPhi_MCseeded.py -O analysis_Tau3Mu_MCseeded.py

change the testFile to:

testFile = "/eos/experiment/fcc/ee/generation/DelphesEvents/spring2021/IDEA/p8_noBES_ee_Ztautau_ecm91_EvtGen_TauMinus2MuMuMu/events_189205650.root"

Modify the call to MCParticle::get_indices to retrieve properly the indices of the decay of interest. Subsequently rename the Bs2MuMuKK_indices into the name you chose --- and, to have meaningful variable names, rename Bsxxx into Tauxxx.

    .Define("Tau3Mu_indices",
            "MCParticle::get_indices(15, {-13, 13, 13}, true, true, true, false) (Particle, Particle1)" )

Exercise 2¶

Add the reconstructed $\tau$ mass to the ntuple (you will need to write new code). Check that the mass resolution is improved when it is determined from the track momenta at the tau decay vertex, compared to a blunt 3-muon mass determined from the default track momenta (taken at the distance of closest approach).

Suggested implementation: analyzers with the following signatures need to be added to your analyzers header file analyzers_Tau3Mu.h (created by copying the analyzers.h from previous section):

double tau3mu_vertex_mass(const VertexingUtils::FCCAnalysesVertex& vertex);
double tau3mu_raw_mass(const ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>& legs);

[...]

#include "TLorentzVector.h"
#include "edm4hep/ReconstructedParticleData.h"
#include "FCCAnalyses/VertexingUtils.h"

[...]

  const double MUON_MASS = 0.1056;  // GeV

  double tau3mu_vertex_mass(const VertexingUtils::FCCAnalysesVertex& vertex) {
    TLorentzVector tau;
    ROOT::VecOps::RVec<TVector3> momenta = vertex.updated_track_momentum_at_vertex;
    int n = momenta.size();
    for (int ileg=0; ileg < n; ileg++) {
      TVector3 track_momentum = momenta[ileg];
      TLorentzVector leg;
      leg.SetXYZM(track_momentum[0], track_momentum[1], track_momentum[2], MUON_MASS);
      tau += leg;
    }

    return tau.M();
  }

  double tau3mu_raw_mass(const ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>& legs) {
    TLorentzVector tau;
    int n = legs.size();
    for (int ileg=0; ileg < n; ileg++) {
      TLorentzVector leg;
      leg.SetXYZM(legs[ileg].momentum.x, legs[ileg].momentum.y, legs[ileg].momentum.z, MUON_MASS);
      tau += leg;
    }

    return tau.M();
  }

    # The reco'ed tau mass --- from the post-VertxFit momenta, at the tau decay
    # vertex:
    .Define("TauMass", "VtxAna::tau3mu_vertex_mass(TauVertexObject)")
    # The "raw" mass - using the track  momenta at their dca:
    .Define("RawMass", "VtxAna::tau3mu_raw_mass(TauRecoParticles)")

includePaths = ["analyzers_Tau3Mu.h"]

fccanalysis run analysis_Tau3Mu_MCseeded.py --test --nevents 1000 --output Tau3Mu_MCseeded.root

root -b Tau3Mu_MCseeded.root

TH1F* h1 = new TH1F("h1", ";Tau Mass (GeV); a.u.", 20, 1.75, 1.8);
events->Draw("TauMass >> h1");
TH1F* h2 = new TH1F("h2", ";Raw Mass (GeV); a.u.", 20, 1.75, 1.8);
events->Draw("RawMass >> h2");
h1->SetLineColor(2);

TCanvas* c = new TCanvas("c", "c");
h1->Draw("hist");
h2->Draw("same, hist");
c->Print("tau_mass.png");

Exercise 3¶

So far, everything was done using “Monte-Carlo seeding”, which gives the resolutions that we expect, in the absence of possible combinatoric issues. The next step is to write a new analysis script which starts from th reconstructed muons.

cp analysis_Tau3Mu_MCseeded.py  analysis_Tau3Mu.py

and keep only the skeleton of the analysers() method such that you delete all of the defines and aliases, the result should look something like this:

class RDFanalysis():
    def analysers(df):
        df2 = (
            df
        )
        return df2

Clear out also the branchList, and insert new defines into the analysers() method:

    # Use the "AllMuons" collection, which contains also non-isolated muons (in
    # contrast to the "Muons" collection)
    #
    # Actually, "Muon" or ("AllMuon") just contain pointers (indices) to the
    # RecoParticle collections, hence one needs to first retrieve the
    # RecoParticles corresponding to these muons (for more detail about the
    # subset collections, see:
    # https://github.com/HEP-FCC/FCCAnalyses/tree/master/examples/basics)
    .Alias("Muon0", "AllMuon#0.index")
    .Define("muons", "ReconstructedParticle::get(Muon0, ReconstructedParticles)")
    .Define("n_muons", "ReconstructedParticle::get_n(muons)")

We now want to write a method that builds muon triplets --- actually, since the MC files produced for this tutorial only forced the decay of the $\tau^-$, we are interested in triplets with total charge equal to -1.

Sub-exercise 1: create a function in your analyzers_Tau3Mu.h that builds such triplet.

  ROOT::VecOps::RVec<ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>>
  build_triplets(const ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>& inParticles,
                 float total_charge) {
    ROOT::VecOps::RVec<ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>> result;

    int n = inParticles.size();
    if (n < 3) {
      return result;
    }

    for (int i = 0; i < n; ++i) {
      edm4hep::ReconstructedParticleData pi = inParticles[i];

      for (int j = i + 1; j < n; ++j) {
        edm4hep::ReconstructedParticleData pj = inParticles[j];

        for (int k=j+1; k < n; ++k) {
          edm4hep::ReconstructedParticleData pk = inParticles[k];
          float charge_tot = pi.charge + pj.charge + pk.charge;

          if (charge_tot == total_charge) {
            ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> a_triplet = {pi, pj, pk};
            result.push_back(a_triplet);
          }
        }  // end of the loop over k
      }  // end of the loop over j
    }  // end of the loop over i

    return result;
  }

    # Build triplets of muons.
    # We are interested in tau- -> mu- mu- mu+ (the MC files produced
    # for this tutorial only forced the decay of the tau- , not the
    # tau+). Hence we look for triples of total charge = -1:
    .Define("triplets_m", "VtxAna::build_triplets(muons, -1.)")
    .Define("n_triplets_m", "return triplets_m.size();")

# ----------------------------------------------------
#    # Considering only the 1st triplet:
#    .Define("the_muons_candidate_0",
#            "return triplets_m[0];")  # the_muons_candidates = a vector of 3 RecoParticles
#    # get the corresponding tracks:
#    .Define("the_muontracks_candidate_0",
#            "ReconstructedParticle2Track::getRP2TRK(the_muons_candidate_0, EFlowTrack_1)")
#    # and fit them to a common vertex:
#    .Define("TauVertexObject_candidate_0",
#            "VertexFitterSimple::VertexFitter_Tk(2, the_muontracks_candidate_0)")
#    # Now we can get the mass of this candidate, as before:
#    .Define("TauMass_candidate_0",
#            "myAnalysis::tau3mu_vertex_mass( TauVertexObject_candidate_0)")

Sub-exercise 2: create a function in your analyzers_Tau3Mu.h to retrieve all tau candidates, and their corresponding tau mass.

ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex>
build_AllTauVertexObject(const ROOT::VecOps::RVec<ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>>& triplets,
                         const ROOT::VecOps::RVec<edm4hep::TrackState>& allTracks) {
  ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex> results;
  int ntriplets = triplets.size();
  for (int i=0; i < ntriplets; i++) {
    ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> legs = triplets[i];

    ROOT::VecOps::RVec<edm4hep::TrackState> the_tracks = ReconstructedParticle2Track::getRP2TRK(legs, allTracks);
    VertexingUtils::FCCAnalysesVertex vertex = VertexFitterSimple::VertexFitter_Tk(2, the_tracks);
    results.push_back(vertex);
  }

  return results;
}


ROOT::VecOps::RVec<double>
build_AllTauMasses(const ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex>& vertices) {
  ROOT::VecOps::RVec<double> results;
  for (auto& v: vertices) {
    double mass = tau3mu_vertex_mass(v);
    results.push_back(mass);
  }

  return results;
}

#include "FCCAnalyses/ReconstructedParticle2Track.h"
#include "FCCAnalyses/VertexFitterSimple.h"

   # ----------------------------------------------------
   # Now consider all triplets :

   .Define("TauVertexObject_allCandidates",
           "VtxAna::build_AllTauVertexObject(triplets_m, EFlowTrack_1)")
   .Define("TauMass_allCandidates",
           "VtxAna::build_AllTauMasses(TauVertexObject_allCandidates)")

    branchList = [
                 "n_muons",
                 "n_triplets_m",
                 "TauMass_allCandidates"
                 ]

fccanalysis run analysis_Tau3Mu.py --test --nevents 1000 --output Tau3Mu.root

TCanvas* c = new TCanvas("c", "c");
// candidates at large mass pick up a muon from the "other" leg
events->Draw("TauMass_allCandidates")
// the genuine tau to 3mu candidates
events->Draw("TauMass_allCandidates", "TauMass_allCandidates < 2")
c->Print("tau_mass.png");

Exercise 4¶

We now want to look at the background. Create a new script file by copying your analysis_Tau3Mu.py into a new file:

cp analysis_Tau3Mu.py analysis_Tau3Mu_stage1.py

Note: The _stage1 suffix indicates that it is a first step in the chain of analysis scripts.

The main background is expected to come from $\tau \rightarrow 3 \pi \nu$ decays, when the charged pions are misidentified as muons. But there is no “fakes” in Delphes: all “Muon” objects that you have in the EDM4hep file do originate from genuine muons (which may, of course, come from a hadron decay). To alleviate this limitation, we first select the ReconstructedParticles that are matched to a stable, charged hadron. Edit your analysis_Tau3Mu_stage1.py script file and add the following lines right before the .Define("n_muons", ... ):

    # -----------------------------------------
    # Add fake muons from pi -> mu
    # This selects the charged hadrons:
    .Alias("MCRecoAssociations0", "MCRecoAssociations#0.index")
    .Alias("MCRecoAssociations1", "MCRecoAssociations#1.index")
    .Define("ChargedHadrons",
            "ReconstructedParticle2MC::selRP_ChargedHadrons(MCRecoAssociations0, MCRecoAssociations1, ReconstructedParticles, Particle)")

(As mentioned earlier, the matching between reconstructed particles and Monte Carlo particles requires 4 collections. See here for more detail).

and further select the ones that are above 2 GeV --- since only particles above 2 GeV will make it to the muon detector:

    # Only the ones with  p > 2 GeV could be selected as muons:
    .Define("ChargedHadrons_pgt2",
            "ReconstructedParticle::sel_p(2.)(ChargedHadrons)")

Now we want to apply a “flat” fake rate, i.e. accept a random fraction of the above particles as muons.

Exercise itself: create an analyzer (functor) in your analyzers_Tau3Mu.h that does that.

  selRP_Fakes::selRP_Fakes(float arg_fakeRate, float arg_mass): m_fakeRate(arg_fakeRate),
                                                                m_mass(arg_mass) {
    auto seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::default_random_engine generator(seed);
    m_generator = generator;
    std::uniform_real_distribution<float> flatdis(0., 1.);
    m_flat.param(flatdis.param());
  };

  ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>
  selRP_Fakes::operator() (const ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData>& inParticles) {
    ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> result;
    for (const auto& p: inParticles) {
      float arandom = m_flat(m_generator);
      if (arandom <= m_fakeRate) {
        edm4hep::ReconstructedParticleData reso = p;
        // overwrite the mass:
        reso.mass = m_mass;
        result.push_back(reso);
      }
    }

    return result;
  }

    # Build fake muons based on a flat fake rate (random selection) ---
    # HUGE fake rate used on purpose here:
    .Define("fakeMuons_5em2",
            ROOT.VtxAna.selRP_Fakes(5e-2, ROOT.VtxAna.MUON_MASS),
            ["ChargedHadrons_pgt2"])
    # Now we merge the collection of fake muons with the genuine muons:
    .Define("muons_with_fakes",
            "ReconstructedParticle::merge( muons, fakeMuons_5em2 )")
    # and we use this collection later on, instead of "muons":
    .Define("n_muons_with_fakes",
            "ReconstructedParticle::get_n(muons_with_fakes)")

    # returns a vector of triplets, i.e. of vectors of 3 RecoParticles
    .Define("triplets_m", "myAnalysis::build_triplets(muons_with_fakes, -1.)")

import ROOT

    # Total visible energy in the event:
    .Define("RecoPartEnergies",
            "ReconstructedParticle::get_e(ReconstructedParticles)")
    .Define("visible_energy", "Sum(RecoPartEnergies)")

    branchList = [
        ...
        "n_muons",
        "n_triplets_m",
        "TauMass_allCandidates",
        "visible_energy"
    ]

fccanalysis run analysis_Tau3Mu_stage1.py --test --nevents 1000 --output Tau3Mu.root

The files analysis_Tau3Mu_stage1.py and analyzers_Tau3Mu.h with all the changes discussed above can be downloaded here and here.

Exercise 5¶

We now have a simple analysis that can be used to process the signal and background samples, and plot the mass of the $\tau \rightarrow 3\mu$ candidates. We can now process the full statistics. In order for you to have access to the dataset metadata provided in the /afs/cern.ch/work/f/fccsw/public/FCCDicts/, you need to have CERN computing account and have the AFS Workspaces service activated.

All samples that have been centrally produced can be found on this web page. We use spring2021 samples (in campaign), and the files made with IDEA detector concept. If you enter TauMinus2MuMuMu and TauMinus2PiPiPinus in the search field, you will see the datasets produced for the signal and the $\tau \rightarrow 3\pi \nu$ background. The first column shows the dataset names, in this case p8_noBES_ee_Ztautau_ecm91_EvtGen_TauMinus2MuMuMu and p8_noBES_ee_Ztautau_ecm91_EvtGen_TauMinus2PiPiPinu.

To run fccanalysis over these datasets (and not anymore over one test file), the list of datasets to be processed should be inserted into your analysis analysis_Tau3Mu_stage1.py:

import ROOT

processList = {
    'p8_noBES_ee_Ztautau_ecm91_EvtGen_TauMinus2MuMuMu': {},
    'p8_noBES_ee_Ztautau_ecm91_EvtGen_TauMinus2PiPiPinu': {}
}

as well as information where to look for the dataset metadata information

# Mandatory: Production tag when running over centrally produced events, this
# points to the yaml files for getting the dataset statistics
prodTag = "FCCee/spring2021/IDEA/"

and we tell fccanalysis to put all files into the Tau3Mu directory:

# Optional: output directory, default is local running directory
outputDir = "Tau3Mu"

This produces flat ntuples:

fccanalysis run analysis_Tau3Mu_stage1.py

After a few minutes, you will see two ntuples, one for the signal, the other for the background, in the Tau3Mu directory.

Now that you have ntuples with the variables you need for your analysis, you may analyze them using the tools you prefer. As explained here, you can also use fccanalysis to produce histograms of some variables, with some sets of cuts, and to plot them. The normalisation is then taken care of by fccanalyses, thanks to a JSON file (FCCee_procDict_spring2021_IDEA.json) which contains the cross-section of the MC processes. To use these functionalities:

The analysis_Tau3Mu_final.py analysis script produces histograms of selected variables, with some selection:

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/analysis_Tau3Mu_final.py
fccanalysis final analysis_Tau3Mu_final.py

# Dictionnay of the list of cuts. The key is the name of the selection that
# will be added to the output file
cutList = {"nocut": "true",
           "sel0": "n_triplets_m > 0",
           "sel1": "Min(TauMass_allCandidates) < 3"
}


# Dictionary for the output variables/hitograms. The key is the name of the
# variable in the output files. "name" is the name of the variable in the
# input file, "title" is the x-axis label of the histogram, "bin" the number
# of bins of the histogram, "xmin" the minimum x-axis value and "xmax" the
# maximum x-axis value.
histoList = {
    "mTau": {
        "name": "TauMass_allCandidates",
        "title": "m_{tau} [GeV]",
        "bin": 100,
        "xmin": 0,
        "xmax": 3
    },
    "mTau_zoom": {
        "name": "TauMass_allCandidates",
        "title": "m_{tau} [GeV]",
        "bin": 100,
        "xmin": 1.7,
        "xmax": 1.9
    },
    "Evis": {
        "name": "visible_energy",
        "title": "E_{vis} [GeV]",
        "bin": 100,
        "xmin": 0,
        "xmax": 91.2
    },
}

This produces three histograms, for three sets of selections. The histogram files appear in the TauMu/final directory. There are now 6 files (2 processes, 3 sets of cuts). The histograms are normalised to a luminosity of one inverse pb.

Finally, the plotting script analysis_Tau3Mu_plots.py makes few plots:

wget https://fccsw.web.cern.ch/fccsw/tutorials/vtx-tutorial/analysis_Tau3Mu_plots.py
fccanalysis plots analysis_Tau3Mu_plots.py

which appear in the Tau3Mu/plots directory. Look for example at the plot mTaTau3Mu_nostack_log in Tau3Mu/plots/sel1. The candidates in the background sample, corresponding to $\tau \rightarrow 3 \pi \nu$ decays, are reconstructed at a mass that is usually below the $\tau$ mass due to the presence of the neutrino. Note that the background sample corresponds to a very low statistics (50M have been produced, we expect 14B) and we see fluctuations. The signal cross-section used here corresponds to $B( \tau \rightarrow 3 \mu) = 2 \times 10^{-8}$, which is roughly the current upper limit. The luminosity for these final plots is entered in analysis_Tau3Mu_plots.py.

Exercise 6¶

Exercise to go beyond

Interested in studying the FCC-ee sensitivity to $\tau \rightarrow 3 \mu$?
Please contact <alberto.lusiani@pi.infn.it> and <monteil@in2p3.fr>.

Note that the input files used in this tutorial are low statistics examples. Larger datasets that are more accurate (KKMC instead of Pythia) can be produced if someone is interested.