[tutorials] Add Python tutorial for RDF h1 analysis

Author: GoodingJamie

tutorials/CMakeLists.txt

@@ -181,6 +181,7 @@ endif()
 
 if(NOT xrootd)
   set(xrootd_veto analysis/dataframe/df101_h1Analysis.C
+                  analysis/dataframe/df101_h1Analysis.py
                   analysis/dataframe/df102_NanoAODDimuonAnalysis.C
                   analysis/dataframe/df103_NanoAODHiggsAnalysis.C
                   analysis/dataframe/df106_HiggsToFourLeptons.C

tutorials/analysis/dataframe/df101_h1Analysis.py

@@ -0,0 +1,122 @@
+import ROOT
+import math
+
+ROOT.EnableImplicitMT(4)
+
+
+# Declare filters on RVec objects and JIT with Numba
+@ROOT.Numba.Declare(["int", "int", "RVecI"], "bool")
+def ik_ipi_nhitrp_cut(ik, ipi, nhitrp):
+    return nhitrp[ik - 1] * nhitrp[ipi - 1] > 1
+
+
+@ROOT.Numba.Declare(["int", "RVecF", "RVecF"], "bool")
+def ik_rstart_rend_cut(ik, rstart, rend):
+    return rend[ik - 1] - rstart[ik - 1] > 22
+
+
+@ROOT.Numba.Declare(["int", "RVecF", "RVecF"], "bool")
+def ipi_rstart_rend_cut(ipi, rstart, rend):
+    return rend[ipi - 1] - rstart[ipi - 1] > 22
+
+
+@ROOT.Numba.Declare(["int", "RVecF"], "bool")
+def ik_nlhk_cut(ik, nlhk):
+    return nlhk[ik - 1] > 0.1
+
+
+@ROOT.Numba.Declare(["int", "RVecF"], "bool")
+def ipi_nlhpi_cut(ipi, nlhpi):
+    return nlhpi[ipi - 1] > 0.1
+
+
+@ROOT.Numba.Declare(["int", "RVecF"], "bool")
+def ipis_nlhpi_cut(ipis, nlhpi):
+    return nlhpi[ipis - 1] > 0.1
+
+
+def select(rdf):
+    return (
+        rdf.Filter("TMath::Abs(md0_d - 1.8646) < 0.04")
+        .Filter("ptds_d > 2.5")
+        .Filter("TMath::Abs(etads_d) < 1.5")
+        .Filter("Numba::ik_ipi_nhitrp_cut(ik, ipi, nhitrp)")
+        .Filter("Numba::ik_rstart_rend_cut(ik, rstart, rend)")
+        .Filter("Numba::ipi_rstart_rend_cut(ipi, rstart, rend)")
+        .Filter("Numba::ik_nlhk_cut(ik, nlhk)")
+        .Filter("Numba::ipi_nlhpi_cut(ipi, nlhpi)")
+        .Filter("Numba::ipis_nlhpi_cut(ipis, nlhpi)")
+    )
+
+
+dxbin = (0.17 - 0.13) / 40
+condition = "x > 0.13957"
+xp3 = "(x - [3]) * (x - [3])"
+
+
+def FitAndPlotHdmd(hdmd: ROOT.TH1):
+    ROOT.gStyle.SetOptFit()
+    c1 = ROOT.TCanvas("c1", "h1analysis analysis", 10, 10, 800, 600)
+
+    hdmd.GetXaxis().SetTitleOffset(1.4)
+
+    hdraw = hdmd.DrawClone()
+
+    # Fit histogram hdmd with function f5 using the loglikelihood option
+    formula = f"{dxbin} * ([0] * pow(x - 0.13957, [1]) + [2] / 2.5066 / [4] * exp(-{xp3} / 2 / [4] / [4]))"
+    f5 = ROOT.TF1("f5", f"{condition} ? {formula} : 0", 0.139, 0.17, 5)
+    f5.SetParameters(1000000, 0.25, 2000, 0.1454, 0.001)
+    hdraw.Fit(f5, "lr")
+
+    return
+
+
+def FitAndPlotH2(h2: ROOT.TH2):
+    # Create the canvas for tau d0
+    c2 = ROOT.TCanvas("c2", "tauD0", 100, 100, 800, 600)
+
+    c2.SetGrid()
+    c2.SetBottomMargin(0.15)
+
+    # Project slices of 2-d histogram h2 along X , then fit each slice
+    # with function f2 and make a histogram for each fit parameter
+    # Note that the generated histograms are added to the list of objects
+    # in the current directory.
+    sigma = 0.0012
+    formula = f"{dxbin} * ([0] * pow(x - 0.13957, 0.25) + [1] / 2.5066 / {sigma} * exp(-{xp3} / 2 / {sigma} / {sigma}))"
+    print(f"TWO: {condition} ? {formula} : 0")
+
+    f2 = ROOT.TF1("f2", f"{condition} ? {formula} : 0", 0.139, 0.17, 2)
+    f2.SetParameters(10000, 10)
+    h2.FitSlicesX(f2, 0, -1, 1, "qln")
+
+    # See TH2::FitSlicesX documentation why h2_1 name is used
+    h2_1 = ROOT.gDirectory.Get("h2_1")
+    h2_1.SetDirectory(ROOT.nullptr)
+    h2_1.GetXaxis().SetTitle("#tau [ps]")
+    h2_1.SetMarkerStyle(21)
+    h2_1.Draw()
+
+    c2.Update()
+
+    line = ROOT.TLine(0, 0, 0, c2.GetUymax())
+    line.Draw()
+
+    return
+
+
+chain = ROOT.TChain("h42")
+chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarmb.root")
+chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarp1a.root")
+chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarp1b.root")
+chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarp2.root")
+
+df = ROOT.RDataFrame(chain)
+selected = select(df)
+# Note: The title syntax is "<Title>;<Label x axis>;<Label y axis>"
+hdmdARP = selected.Histo1D(("hdmd", "Dm_d;m_{K#pi#pi} - m_{K#pi}[GeV/c^{2}]", 40, 0.13, 0.17), "dm_d")
+selected_added_branch = selected.Define("h2_y", "rpd0_t / 0.029979f * 1.8646f / ptd0_d")
+h2ARP = selected_added_branch.Histo2D(("h2", "ptD0 vs Dm_d", 30, 0.135, 0.165, 30, -3, 6), "dm_d", "h2_y")
+
+FitAndPlotHdmd(hdmdARP)
+FitAndPlotH2(h2ARP)

tutorials/analysis/dataframe/index.md

@@ -129,7 +129,7 @@ With RDataFrame advanced analyses can be executed on large amounts of data. Thes
 | **Tutorial** || **Description** |
 |------|--------|-----------------|
 | df017_vecOpsHEP.C | df017_vecOpsHEP.py | Use RVecs to plot the transverse momentum of selected particles. |
-| df101_h1Analysis.C | | Express ROOT's standard H1 analysis. |
+| df101_h1Analysis.C | df101_h1Analysis.py | Express ROOT's standard H1 analysis. |
 | df102_NanoAODDimuonAnalysis.C | df102_NanoAODDimuonAnalysis.py | Process  NanoAOD files. |
 | df103_NanoAODHiggsAnalysis.C | df103_NanoAODHiggsAnalysis.py | An example of complex analysis: reconstructing the Higgs boson. |
 | | df104_HiggsToTwoPhotons.py | The Higgs to two photons analysis from the ATLAS Open Data 2020 release. |
@@ -138,4 +138,4 @@ With RDataFrame advanced analyses can be executed on large amounts of data. Thes
 | | df107_SingleTopAnalysis.py | A single top analysis using the ATLAS Open Data release of 2020. |
 
 \anchor df_alltutorials
-@}
+@}