diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt index 4c597d4964271..8cfa9f5016e40 100644 --- a/tutorials/CMakeLists.txt +++ b/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 diff --git a/tutorials/analysis/dataframe/df101_h1Analysis.py b/tutorials/analysis/dataframe/df101_h1Analysis.py new file mode 100644 index 0000000000000..96cb4ad8509a3 --- /dev/null +++ b/tutorials/analysis/dataframe/df101_h1Analysis.py @@ -0,0 +1,123 @@ +import ROOT + +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") + + c1.Update() + + 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 "