Run5

CODEOWNERS

@@ -0,0 +1,9 @@
+# Code owners
+#
+# Lines starting with '#' are comments.
+# Each line is a file pattern followed by one or more owners.
+
+# These owners will be the default owners for everything in the repo.
+* @ginnocen @vkucera
+
+# Order is important. The last matching pattern has the most precedence.

FirstAnalysis/README.md

@@ -1,19 +1,29 @@
-# O2 First analysis utilities
-
-## Table of contents
-
-* [Introduction](#introduction)
-* [Available samples ALICE3](#available-samples-alice3)
-
+# O2 first analysis utilities
 
 ## Introduction
 The main purpose of this folder is to provide analysis tools to perform analysis of O2 analysis task outputs.
 It includes for the moment a set of QA plotting macros to inspect `AnalysisResults.root` of MC and data samples.
-It also includes some preliminary functionalities for performing efficiency corrections of HF analyses. It is currently being used
-also for performing Run5 analysis.
+It also includes some preliminary functionalities for performing efficiency corrections of HF analyses.
+It is currently being used also for performing Run5 analysis.
 
+## Instructions for contributors
+
+* Write documentation comments to make the code easy to understand.
+* Every file should contain the following details at the top:
+  * what the code does,
+  * which input it needs (e.g. output of which O2 workflow),
+  * how to run it (in case simply executing the script is not enough),
+  * names of authors, contributors, maintainers (to contact in case of problems or questions).
+* Follow [PEP 8 naming conventions](https://www.python.org/dev/peps/pep-0008/#naming-conventions) in Python code.
+* Do not use hard-coded paths! A default path should point to a file directly produced by the framework (e.g. `../codeHF/AnalysisResults_O2.root`).
 
 ## Available samples ALICE3
 * D2H KrKr production (X events):
   * AOD O2 tables: `/data/Run5data/EMBEDDING_KrKr_CCBARLc_scenario3/TEST`
   * Analysis result: `/data/Run5data_samples/EMBEDDING_KrKr_CCBARLc_scenario3_20210224/AnalysisResults_O2.root`
+
+## Plotting
+The analyses validation is supposed to use the [HFPlot tool](https://github.com/benedikt-voelkel/HFPlot) for plotting purposes in the future.
+It is basically a *pythonic* wrapper around `ROOT` with some automated features.
+As a first use case it is implemented in [distrib_studies.py](distrib_studies.py).
+A more extensive README can be found [here](https://github.com/benedikt-voelkel/HFPlot) and there are also some examples provided.

FirstAnalysis/contamination.py

@@ -0,0 +1,149 @@
+#!/usr/bin/env python
+
+from ROOT import TCanvas, TFile, gROOT, gStyle
+
+path_file = "../codeHF/AnalysisResults_O2.root"
+
+gStyle.SetOptStat(0)
+gStyle.SetErrorX(0)
+gStyle.SetFrameLineWidth(1)
+gStyle.SetTitleSize(0.045, "x")
+gStyle.SetTitleSize(0.045, "y")
+gStyle.SetMarkerSize(1)
+gStyle.SetLabelOffset(0.015, "x")
+gStyle.SetLabelOffset(0.02, "y")
+gStyle.SetTickLength(-0.02, "x")
+gStyle.SetTickLength(-0.02, "y")
+gStyle.SetTitleOffset(1.1, "x")
+gStyle.SetTitleOffset(1.0, "y")
+
+gROOT.SetBatch(True)
+
+
+def saveCanvas(canvas, title):
+    format_list = [".pdf", ".png"]
+    for fileFormat in format_list:
+        canvas.SaveAs(title + fileFormat)
+
+
+def kinematic_plots(var, particle, detector, hp):
+    fileo2 = TFile(path_file)
+    cres = TCanvas("cres", "resolution distribution")
+    cres.SetCanvasSize(1600, 1000)
+    cres.cd()
+    num = fileo2.Get(
+        "qa-tracking-rejection-%s/tracking%ssel%s/%seta" % (particle, detector, hp, var)
+    )
+    den = fileo2.Get("qa-tracking-rejection-%s/tracking/%seta" % (particle, var))
+    # gPad.SetLogz()
+    num.Divide(den)
+    num.Draw("coltz")
+    num.GetYaxis().SetTitle("#eta")
+    num.GetYaxis().SetTitleOffset(1.0)
+    num.GetZaxis().SetRangeUser(0.01, 1)
+    num.SetTitle("Fraction of %s selected by %s as %s" % (particle, detector, hp))
+    saveCanvas(
+        cres, "contamination/%seta_%sSelfrom%sas%s" % (var, particle, detector, hp)
+    )
+
+
+def ratioparticle(
+    var="pt",
+    numname="Electron",
+    selnum="RICHSelHpElLoose",
+    denname="Pion",
+    selden="RICHSelHpElLoose",
+    label="e/#pi",
+):
+    fileo2 = TFile(path_file)
+    cres = TCanvas("cres", "resolution distribution")
+    cres.SetCanvasSize(1600, 1000)
+    cres.cd()
+    cres.SetLogy()
+
+    num2d = fileo2.Get("qa-rejection-general/h%s%s/%seta" % (numname, selnum, var))
+    den2d = fileo2.Get("qa-rejection-general/h%s%s/%seta" % (denname, selden, var))
+    num = num2d.ProjectionX("num", 1, num2d.GetXaxis().GetNbins())
+    den = den2d.ProjectionX("den", 1, den2d.GetXaxis().GetNbins())
+    num.Divide(den)
+    num.Draw("coltz")
+    num.GetYaxis().SetTitle(label)
+    num.GetXaxis().SetTitle("p_{T}")
+    # num.SetMinimum(0.001)
+    # num.SetMaximum(2.0)
+    num.GetYaxis().SetTitleOffset(1.0)
+    # num.GetZaxis().SetRangeUser(0.01, 1)
+    # nameresult = "Fraction of %s selected by %s over %s selected by %s" % (num,selnum,den,selden)
+    canvas = "Fractionof%s%sOver%s%s" % (numname, selnum, denname, selden)
+    # num.SetTitle(nameresult)
+    saveCanvas(cres, "rejection/%s" % (canvas))
+
+
+def is_e_not_pi_plots(particle):
+    fileo2 = TFile(path_file)
+    task = "qa-rejection-general"
+    folder_gm = "h%sRICHSelHpElTight" % particle
+    folder_alt = "h%sRICHSelHpElTightAlt" % particle
+    folder_diff = "h%sRICHSelHpElTightAltDiff" % particle
+    hist = "pteta"
+    hist_gm = fileo2.Get("%s/%s/%s" % (task, folder_gm, hist))
+    hist_gm.SetTitle("%s isRICHElTight" % particle)
+    hist_alt = fileo2.Get("%s/%s/%s" % (task, folder_alt, hist))
+    hist_alt.SetTitle("%s isElectronAndNotPion" % particle)
+    hist_diff = fileo2.Get("%s/%s/%s" % (task, folder_diff, hist))
+    hist_diff.SetTitle("%s isRICHElTight != isElectronAndNotPion" % particle)
+    cepi = TCanvas("cepi", "e not pi selection")
+    cepi.SetCanvasSize(1600, 1000)
+    cepi.Divide(2, 2)
+    cepi.cd(1)
+    hist_gm.Draw("colz")
+    cepi.cd(2)
+    hist_alt.Draw("colz")
+    cepi.cd(3)
+    hist_diff.Draw("colz")
+    # num.GetYaxis().SetTitleOffset(1.0)
+    # num.GetZaxis().SetRangeUser(0.01, 1)
+    saveCanvas(cepi, "contamination/is_e_not_pi_%s" % particle)
+
+
+# kinematic_plots("p", "pion", "MID", "Muon")
+# kinematic_plots("p", "mu", "MID", "Muon")
+# kinematic_plots("p", "pion", "TOF", "Electron")
+# kinematic_plots("p", "pion", "RICH", "Electron")
+# kinematic_plots("pt", "pion", "RICH", "Electron")
+# kinematic_plots("p", "kaon", "RICH", "Electron")
+# kinematic_plots("pt", "kaon", "RICH", "Electron")
+# kinematic_plots("p", "pion", "TOF", "Kaon")
+# kinematic_plots("p", "pion", "RICH", "Kaon")
+
+ratioparticle(
+    "pt", "Electron", "RICHSelHpElTight", "Electron", "NoSel", "e/e RICHSelHpElTight"
+)
+ratioparticle("pt", "Electron", "NoSel", "Pion", "NoSel", "e/#pi No cuts")
+ratioparticle(
+    "pt", "Electron", "RICHSelHpElTight", "Pion", "RICHSelHpElTight", "Tight e/#pi "
+)
+ratioparticle("pt", "Muon", "MID", "Pion", "MID", "MIDSel")
+ratioparticle("pt", "Pion", "RICHSelHpElTight", "Pion", "NoSel", "Contamination")
+ratioparticle("pt", "Pion", "MID", "Pion", "NoSel", "Contamination MID")
+
+ratioparticle(
+    "pt",
+    "Electron",
+    "RICHSelHpElTightAlt",
+    "Electron",
+    "RICHSelHpElTight",
+    "e isElectronAndNotPion/RICHSelHpElTight",
+)
+ratioparticle(
+    "pt",
+    "Electron",
+    "RICHSelHpElTightAlt",
+    "Pion",
+    "RICHSelHpElTightAlt",
+    "isElectronAndNotPion e/#pi",
+)
+
+
+for p in ("Electron", "Pion", "Kaon", "Muon"):
+    is_e_not_pi_plots(p)

FirstAnalysis/database.yaml

@@ -8,7 +8,7 @@ Xi_cc:
       normalized: 1
       varlist: ["mass", "d0Prong0", "d0Prong1", "CPA", "Y", "DecLength", "Ct", "Eta", "chi2PCA"]
       varlatex: ["inv. mass (GeV)", "d_{0} prong 0 (cm)", "d_{0} prong 1 (cm)", "cos(#theta_{p})", "rapidity", "dec. length (cm)", "c#tau (cm)", "#eta", "chi^{2} decay"]
-      histonamesig: ["hmassSig", "hd0Prong0Sig", "hd0Prong1Sig", "hCPASig", "hYSig", "hDecLengthSig", "hCtSig", "hEtaSig", "hChi2PCABg"]
+      histonamesig: ["hmassSig", "hd0Prong0Sig", "hd0Prong1Sig", "hCPASig", "hYSig", "hDecLengthSig", "hCtSig", "hEtaSig", "hChi2PCASig"]
       histonamebkg: ["hmassBg", "hd0Prong0Bg", "hd0Prong1Bg", "hCPABg", "hYBg", "hDecLengthBg", "hCtBg", "hEtaBg", "hChi2PCABg"]
       ymin: [0.0001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
       ymax: [0.4, 5., 5., 5., 5., 5., 10., 5., 5.]
@@ -18,3 +18,45 @@ Xi_cc:
       dolog: [0, 1, 1, 1, 1, 1, 1, 1, 1]
       dologx: [0, 1, 0, 0, 0, 1, 0, 0, 1]
       legtoleft: [1, 1, 1, 1, 1, 1, 1, 1, 1]
+
+X3872:
+  pp14p0:
+    absy1p44:
+      latexcand: "X(3872)"
+      inputBkg: /data/Run5/results/Xanalysis_0p7pt_RICHTOF/AnalysisResults_O2_MB.root
+      inputSig: /data/Run5/results/Xanalysis_0p7pt_RICHTOF/AnalysisResults_O2_enriched.root
+      dirname: hf-task-x-mc
+      normalized: 1
+      varlist: ["mass", "d0Prong0", "d0Prong1", "CPA", "Y", "DecLength", "Ct", "Eta", "chi2PCA"]
+      varlatex: ["inv. mass (GeV)", "d_{0} prong 0 (cm)", "d_{0} prong 1 (cm)", "cos(#theta_{p})", "rapidity", "dec. length (cm)", "c#tau (cm)", "#eta", "chi^{2} decay"]
+      histonamesig: ["hMassRecSig", "hd0Prong0RecSig", "hd0Prong1RecSig", "hCPARecSig", "hYSig", "hDeclengthRecSig", "hCtSig", "hEtaRecSig", "hChi2PCASig"]
+      histonamebkg: ["hMassRecBg", "hd0Prong0RecBg", "hd0Prong1RecBg", "hCPARecBg", "hYBg", "hDeclengthRecBg", "hCtBg", "hEtaRecBg", "hChi2PCABg"]
+      ymin: [0.0001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
+      ymax: [0.4, 5., 5., 5., 5., 5., 10., 5., 5.]
+      xmin: [3.2, -0.01, -0.01, 0.5, -3., 0.00001, 0., -3., 1.e-6]
+      xmax: [4., 0.04, 0.04, 1.1, 3., 3., 0.1, 3., 0.1]
+      rebin: [5, 2, 5, 5, 5, 5, 5, 5, 5]
+      dolog: [0, 1, 1, 1, 1, 1, 1, 1, 1]
+      dologx: [0, 1, 0, 0, 0, 1, 0, 0, 1]
+      legtoleft: [1, 1, 1, 1, 1, 1, 1, 1, 1]
+
+Jpsi:
+  pp14p0:
+    absy1p44:
+      latexcand: "Jpsi"
+      inputBkg: /data/Run5/results/Xanalysis_0p7pt_RICHTOF/AnalysisResults_O2_MB.root
+      inputSig: /data/Run5/results/Xanalysis_0p7pt_RICHTOF/AnalysisResults_O2_enriched.root
+      dirname: hf-task-jpsi-mc
+      normalized: 1
+      varlist: ["mass", "d0Prong0", "d0Prong1", "Y", "DecLength", "Ct", "d0xd0", "chi2PCA"]
+      varlatex: ["inv. mass (GeV)", "d_{0} prong 0 (cm)", "d_{0} prong 1 (cm)", "rapidity", "dec. length (cm)", "c#tau (cm)", "d_{0} x d_{1}", "chi^{2} decay"]
+      histonamesig: ["hmassSig", "hd0Prong0Sig", "hd0Prong1Sig", "hYSig", "hdeclengthSig", "hCtSig", "hd0d0Sig", "hChi2PCASig"]
+      histonamebkg: ["hmassBg", "hd0Prong0Bg", "hd0Prong1Bg", "hYBg", "hdeclengthBg", "hCtBg", "hd0d0Bg", "hChi2PCABg"]
+      ymin: [0.0001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001]
+      ymax: [0.4, 5., 5., 5., 5., 10., 5., 5.]
+      xmin: [3.2, -0.01, -0.01, -3., 0.00001, 0., -0.02, 1.e-6]
+      xmax: [4., 0.04, 0.04, 3., 3., 0.1, 0.02, 0.1]
+      rebin: [5, 5, 5, 5, 1, 5, 5, 1]
+      dolog: [0, 1, 1, 1, 1, 1, 1, 1]
+      dologx: [0, 1, 0, 0, 1, 0, 0, 1]
+      legtoleft: [1, 1, 1, 1, 1, 1, 1, 1]

FirstAnalysis/distrib_studies.py

@@ -3,10 +3,10 @@
 from math import ceil, sqrt
 
 import yaml
+from hfplot.plot_spec_root import ROOTFigure
+from hfplot.style import StyleObject1D
 from ROOT import TFile
 
-from hfplot.plot_spec_root import ROOTFigure
-from hfplot.style import ROOTStyle1D
 
 def makeSavePaths(title, *fileFormats, outputdir="outputPlots"):
     """
@@ -17,7 +17,7 @@ def makeSavePaths(title, *fileFormats, outputdir="outputPlots"):
     return [outputdir + "/" + title + fileFormat for fileFormat in fileFormats]
 
 
-def distr_studies(hadron="Xi_cc", collision="pp14p0", yrange="absy1p44"):
+def distr_studies(hadron="X3872", collision="pp14p0", yrange="absy1p44"):
     """
     Make distribution comparisons
     """
@@ -47,17 +47,16 @@ def distr_studies(hadron="Xi_cc", collision="pp14p0", yrange="absy1p44"):
     lptMax = []
 
     # Define some styles
-    style_sig = ROOTStyle1D()
+    style_sig = StyleObject1D()
     style_sig.markercolor = 2
     style_sig.markerstyle = 21
-    style_sig.markersize = 1
+    style_sig.markersize = 2
     style_sig.draw_options = "P"
-    style_bkg = ROOTStyle1D()
+    style_bkg = StyleObject1D()
     style_bkg.markerstyle = 23
-    style_bkg.markersize = 1
+    style_bkg.markersize = 2
     style_bkg.draw_options = "P"
 
-
     for index, var in enumerate(lvarlist):
 
         lhistosigvar = []
@@ -66,12 +65,15 @@ def distr_studies(hadron="Xi_cc", collision="pp14p0", yrange="absy1p44"):
         hbkg = fileBkg.Get(f"{dirname}/{lhistonamebkg[index]}")
 
         nPtBins = hsig.GetNbinsY()
+        print(var)
         for iptBin in range(nPtBins):
             # Collect the histogram projections in bins of pT for each variable
-            hsig_px = hsig.ProjectionX(f"hsig_px_var{var}_pt{iptBin}",
-                                       iptBin + 1, iptBin + 1)
-            hbkg_px = hbkg.ProjectionX(f"hbkg_px_var{var}_pt{iptBin}",
-                                       iptBin + 1, iptBin + 1)
+            hsig_px = hsig.ProjectionX(
+                f"hsig_px_var{var}_pt{iptBin}", iptBin + 1, iptBin + 1
+            )
+            hbkg_px = hbkg.ProjectionX(
+                f"hbkg_px_var{var}_pt{iptBin}", iptBin + 1, iptBin + 1
+            )
             lhistosigvar.append(hsig_px)
             lhistobkgvar.append(hbkg_px)
             if index == 0:
@@ -86,9 +88,15 @@ def distr_studies(hadron="Xi_cc", collision="pp14p0", yrange="absy1p44"):
     for index, var in enumerate(lvarlist):
         # sort out required number of columns and rows for a squared grid and create a figure
         n_cols_rows = ceil(sqrt(nPtBins))
-        figure = ROOTFigure(n_cols_rows, n_cols_rows, row_margin=0.05, column_margin=0.05, size=(1500, 900))
+        figure = ROOTFigure(
+            n_cols_rows,
+            n_cols_rows,
+            row_margin=0.04,
+            column_margin=0.04,
+            size=(1500, 1100),
+        )
         # can adjust some axis properties globally
-        figure.axes(label_size=0.02, title_size=0.02)
+        figure.axes(label_size=0.015, title_size=0.015)
         # here we use the feature to only apply to certain axes
         figure.axes("x", title=lvarlatex[index])
         figure.axes("y", title="Entries")
@@ -109,7 +117,9 @@ def distr_studies(hadron="Xi_cc", collision="pp14p0", yrange="absy1p44"):
             nBkgEntries = hist_bkg.Integral()
 
             if not nSigEntries or not nBkgEntries:
-                print(f"ERROR: Found empty signal or background distribution for variable={var} in pT bin={iptBin}")
+                print(
+                    f"ERROR: Found empty signal or background distribution for variable={var} in pT bin={iptBin}"
+                )
                 continue
 
             if normalized:
@@ -124,12 +134,27 @@ def distr_studies(hadron="Xi_cc", collision="pp14p0", yrange="absy1p44"):
                     hist_bkg.SetBinError(ibin + 1, 0.0)
 
             figure.define_plot(x_log=ldologx[index], y_log=ldolog[index])
-            figure.add_object(hist_sig, style=style_sig, label=f"Sig before norm ({int(nSigEntries)} entries)")
-            figure.add_object(hist_bkg, style=style_bkg, label=f"Bkg before norm ({int(nBkgEntries)} entries)")
-            figure.add_text(f"{lptMin[iptBin]:.1f} GeV < p_{{T}} ({latexcand}) < {lptMax[iptBin]:.1f} GeV", 0.1, 0.1)
+            figure.add_object(
+                hist_sig,
+                style=style_sig,
+                label=f"Sig before norm ({int(nSigEntries)} entries)",
+            )
+            figure.add_object(
+                hist_bkg,
+                style=style_bkg,
+                label=f"Bkg before norm ({int(nBkgEntries)} entries)",
+            )
+            figure.add_text(
+                f"{lptMin[iptBin]:.1f} GeV < p_{{T}} ({latexcand}) < {lptMax[iptBin]:.1f} GeV",
+                0.1,
+                0.85,
+            )
 
         figure.create()
-        for save_paths in makeSavePaths(f"distribution_{var}", *formats, outputdir=f"output_{hadron}"):
+        for save_paths in makeSavePaths(
+            f"distribution_{var}", *formats, outputdir=f"output_{hadron}"
+        ):
             figure.save(save_paths)
 
+
 distr_studies()

FirstAnalysis/efficiency_studies.py

@@ -11,13 +11,7 @@ def saveCanvas(canvas, title):
 def efficiencytracking(var):
     # plots the efficiency vs pT, eta and phi for all the species(it extracts the
     # Efficiency from qa - tracking - efficiency if you have ran with-- make - eff)
-    hadron_list = [
-        "pion",
-        "proton",
-        "kaon",
-        "electron",
-        "muon",
-    ]
+    hadron_list = ["pion", "proton", "kaon", "electron", "muon"]
     color_list = [1, 2, 4, 6, 8]
     marker_list = [20, 21, 22, 34, 45]
     fileo2 = TFile("../codeHF/AnalysisResults_O2.root")