feat: draw Cherenkov angle vs. p curves

benchmarks/rich/draw_benchmark.py

@@ -3,11 +3,26 @@
 # Subject to the terms in the LICENSE file found in the top-level directory.
 
 import ROOT as r
-import sys, getopt, pathlib
+import sys, getopt, pathlib, math
 
 # suppress graphics
 r.gROOT.SetBatch(True)
 
+# GLOBAL SETTINGS
+################################################################
+RESID_MAX = 10 # cherenkov angle |residual| maximum
+RADIATORS = {
+        'Aerogel': { 'p_max': 30, 'p_rebin': 12, 'rindex_ref': 1.0190},
+        'Gas':     { 'p_max': 70, 'p_rebin': 30, 'rindex_ref': 1.00076},
+        'Merged':  { 'p_max': 70, 'p_rebin': 30, },
+        }
+PARTICLE_MASS = {
+        'e':  0.00051,
+        'pi': 0.13957,
+        'K':  0.49368,
+        'p':  0.93827,
+        }
+
 # ARGUMENTS
 ################################################################
 
@@ -62,18 +77,50 @@ canv_dict = {
 }
 
 # helper functions
+### minimum momentum for Cherenkov radiation
+def calculate_mom_min(mass,n):
+    return mass / math.sqrt(n**2-1)
+### Cherenkov angle at saturation
+def calculate_theta_max(n):
+    return 1000 * math.acos(1/n)
+### execute `op` for every pad of `canv`
 def loop_pads(canv, op):
     for pad in canv.GetListOfPrimitives():
         if pad.InheritsFrom(r.TVirtualPad.Class()):
             op(pad)
-def draw_profile(hist, style="COLZ"):
+### draw a profile on 2D histogram `hist`
+def draw_profile(hist, rad, style="BOX"):
     prof = hist.ProfileX("_pfx", 1, -1, "i")
-    prof.SetLineColor(r.kBlack)
-    prof.SetLineWidth(3)
-    hist.Draw(style)
+    prof.SetMarkerStyle(r.kFullCircle)
+    prof.SetMarkerColor(r.kRed)
+    prof.SetMarkerSize(2)
+    hist_rebinned = hist.Clone(f'{hist.GetName()}_rebin')
+    hist_rebinned.RebinX(RADIATORS[rad]['p_rebin'])
+    hist_rebinned.SetLineColor(r.kBlue)
+    hist_rebinned.SetFillColor(r.kBlue)
+    for p in [hist_rebinned, prof]:
+        p.GetXaxis().SetRangeUser(0, RADIATORS[rad]['p_max'])
+        if 'rindex_ref' in RADIATORS[rad] and 'theta_vs_p' in hist.GetName():
+            p.GetYaxis().SetRangeUser(0, 1.5*calculate_theta_max(RADIATORS[rad]['rindex_ref']))
+    hist_rebinned.Draw(style)
     prof.Draw("SAME")
 
-################################################################
+# Cherenkov angle vs. p functions
+for rad, radHash in RADIATORS.items():
+    if rad == "Merged":
+        continue
+    radHash['cherenkov_curves'] = []
+    n = radHash['rindex_ref']
+    for part, mass in PARTICLE_MASS.items():
+        ftn = r.TF1(
+                f'ftn_theta_{part}_{rad}',
+                f'1000*TMath::ACos(TMath::Sqrt(x^2+{mass}^2)/({n}*x))',
+                calculate_mom_min(mass, n),
+                radHash['p_max']
+                )
+        ftn.SetLineColor(r.kBlack)
+        ftn.SetLineWidth(1)
+        radHash['cherenkov_curves'].append(ftn)
 
 # draw photon spectra
 canv = canv_dict["photon_spectra"]
@@ -122,18 +169,28 @@ for rad in ["Aerogel", "Gas", "Merged"]:
     hist.DrawCopy()
     canv.cd(4)
     hist.SetTitle(hist.GetTitle() + " - ZOOM")
-    RESID_MAX = 10
     hist.GetXaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
     hist.Draw()
-    hist.Fit("gaus", "", "", -RESID_MAX, RESID_MAX)
+    resid_mode = hist.GetBinCenter(hist.GetMaximumBin())
+    resid_dev  = hist.GetStdDev()
+    nsigma = 3
+    hist.Fit("gaus", "", "", resid_mode - nsigma*resid_dev, resid_mode + nsigma*resid_dev)
+    hist.GetFunction("gaus").SetLineWidth(5)
     canv.cd(5)
     ana_file.Get(f'{pid_name}/highestWeight_dist_{rad}').Draw()
     canv.cd(6)
-    draw_profile(ana_file.Get(f'{pid_name}/npe_vs_p_{rad}'))
+    hist = ana_file.Get(f'{pid_name}/npe_vs_p_{rad}')
+    draw_profile(hist, rad)
     canv.cd(7)
-    draw_profile(ana_file.Get(f'{pid_name}/theta_vs_p_{rad}'))
+    hist = ana_file.Get(f'{pid_name}/theta_vs_p_{rad}')
+    draw_profile(hist, rad)
+    if 'cherenkov_curves' in RADIATORS[rad]:
+        for ftn in RADIATORS[rad]['cherenkov_curves']:
+            ftn.Draw("SAME")
     canv.cd(8)
-    draw_profile(ana_file.Get(f'{pid_name}/thetaResid_vs_p_{rad}'))
+    hist = ana_file.Get(f'{pid_name}/thetaResid_vs_p_{rad}')
+    hist.GetYaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
+    draw_profile(hist, rad)
     canv.cd(9)
     ana_file.Get(f'{pid_name}/photonTheta_vs_photonPhi_{rad}').Draw("COLZ")
     canv.cd(10)