[tutorials] Add Python tutorial for RDF h1 analysis

Author: GoodingJamie

tutorials/analysis/dataframe/df101_h1Analysis.py

@@ -0,0 +1,124 @@
+
+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 # Bin-width
+
+def fdm5(xx, par):
+    x = xx[0]
+    if x <= 0.13957:
+        return 0
+    xp3 = (x - par[3]) ** 2
+    res = dxbin * (par[0] * (x - 0.13957) ** par[1]) + par[2] / 2.5066 / par[4] * math.exp(-xp3 / 2 / par[4] ** 2)
+    return res
+
+def fdm2(xx, par):
+    sigma = 0.0012
+    x = xx[0]
+    if x <= 0.13957:
+        return 0
+    xp3 = (x - 0.1454) ** 2
+    res = dxbin * (par[0] * (x - 0.13957) ** 0.25) + par[1] / 2.5066 / sigma * math.exp(-xp3 / 2 / sigma **2)
+    return res
+
+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
+    f5 = ROOT.TF1("f5", fdm5, 0.139, 0.17, 5)
+    f5.SetParameters(1000000, .25, 2000, .1454, .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.
+    
+    f2 = ROOT.TF1("f2", fdm2, 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
-@}
+@}