ATLAS Offline Software
yearwise_luminosity.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 # Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 
4 """
5 Plot comparisons of Zee/Zmumu and Z/ATLAS over entire data-periods.
6 This can be done as a function of time or pileup
7 """
8 
9 import numpy as np
10 import ROOT as R
11 import python_tools as pt
12 import math
13 from array import array
14 import argparse
15 import csv
16 
17 parser = argparse.ArgumentParser()
18 parser.add_argument('--year', type=str, help='15,16,17,18,22,23,24,25 or run3 or year1_year2_...')
19 parser.add_argument('--channel', type=str, help='Zee or Zmumu')
20 parser.add_argument('--comp', action='store_true', help='Compare Zee and Zmumu?')
21 parser.add_argument('--absolute', action='store_true', help='Compare absolute luminosity')
22 parser.add_argument('--indir', type=str, help='Input directory for CSV files')
23 parser.add_argument('--outdir', type=str, help='Output directory for plots')
24 parser.add_argument('--outcsv', action='store_true', help='Create short CSV with plot content')
25 
26 args = parser.parse_args()
27 year = args.year
28 channel = args.channel
29 absolute = args.absolute
30 indir = args.indir
31 outdir = args.outdir
32 outcsv = args.outcsv
33 
34 # Do all of the ugly plot stlying here
35 yval = 0.85
36 xval = 0.2
37 if args.absolute: ymin, ymax = 0.91, 1.09
38 else: ymin, ymax = 0.93, 1.07
39 
40 if year == "run3":
41  years = ["22", "23", "24", "25"]
42  out_tag = "run3"
43  time_format = "%m/%y"
44  xtitle = 'Month / Year'
45  date_tag = "Run 3,#kern[-0.5]{ }#sqrt{s} = 13.6 TeV"
46  labelsize = 44
47  norm_type = "Run3"
48  multiyear = True
49 elif len(year.split("_")) > 1:
50  years = year.split("_")
51  out_tag = "data"+year
52  time_format = "%m/%y"
53  xtitle = 'Month / Year'
54  date_tag = "Run 3,#kern[-0.5]{ }#sqrt{s} = 13.6 TeV"
55  labelsize = 44
56  norm_type = "Run3"
57  multiyear = True
58 else:
59  years = [year]
60  out_tag = "data"+year
61  time_format = "%d/%m"
62  xtitle = 'Date in 20' + year
63  date_tag = "Data 20" + year + ",#kern[-0.5]{ }#sqrt{s} = 13.6 TeV"
64  norm_type = "year"
65  labelsize = 22
66  multiyear = False
67 
68 def main():
69  if args.comp:
70  channel_comparison(years)
71  else:
72  zcounting_vs_atlas(channel, years)
73 
74 def channel_comparison(years):
75 
76  print("Lumi Channel Comparison vs Time for years: ", years)
77 
78  dict_zlumi = {}
79  for year in years:
80 
81  grl = pt.get_grl(year)
82 
83  for channel in ["Zee", "Zmumu"]:
84 
85  for run in grl:
86  livetime, zlumi, zerr, olumi, timestamp, dfz_small = pt.get_dfz(args.indir, year, run, channel)
87  # Cut out short runs
88  if livetime < pt.runlivetimecut:
89  if livetime >= 0.: print(f"Skip Run {run} because of live time {livetime/60:.1f} min")
90  continue
91 
92  dict_zlumi[channel, run] = (zlumi, zerr, timestamp)
93 
94  vec_times = array('d')
95  vec_ratio = array('d')
96  vec_ratio_err = array('d')
97  keys = [key[1] for key in dict_zlumi if "Zee" in key]
98 
99  # If plotting vs. date simply calculate integrated lumi per run and fill array
100  for key in sorted(keys):
101  try:
102  ratio = dict_zlumi["Zee", key][0]/dict_zlumi["Zmumu", key][0]
103  error = ratio * math.sqrt( pow(dict_zlumi["Zee", key][1]/dict_zlumi["Zee", key][0], 2) + pow(dict_zlumi["Zmumu", key][1]/dict_zlumi["Zmumu", key][0], 2) )
104  date = dict_zlumi["Zee", key][2]
105 
106  if ratio < ymin or ratio > ymax:
107  print("WARNING: Run", key, "has Zee/Zmumu ratio %.2f +- %.2f, outside of y-axis range" % (ratio, error))
108  else:
109  vec_times.append(date)
110  vec_ratio.append(ratio)
111  vec_ratio_err.append(error)
112  except KeyError:
113  print("Cannot do ratio for", key)
114 
115  tg = R.TGraphErrors(len(vec_times), vec_times, vec_ratio, R.nullptr, vec_ratio_err)
116  leg = R.TLegend(0.645, 0.72, 0.805, 0.91)
117  leg.SetFillStyle(0)
118 
119  # Depending if we're plotting over whole Run-3, change canvas size
120  if multiyear:
121  c1 = R.TCanvas("c1", "c1", 2000, 1000)
122  else:
123  c1 = R.TCanvas()
124 
125  tg.Draw('ap')
126  tg.GetYaxis().SetTitle(pt.Leemumuratiolabel)
127  tg.Fit('pol0', '0q')
128  tg.GetFunction('pol0').SetLineColor(R.kRed)
129 
130  mean = tg.GetFunction('pol0').GetParameter(0)
131 
132  # Plot 68% percentile band
133  stdev = np.percentile(abs(vec_ratio - np.median(vec_ratio)), 68)
134  line1 = pt.make_bands(vec_times, stdev, mean)
135  line1.Draw("same 3")
136  tg.GetFunction('pol0').Draw("same l")
137  tg.Draw('same e0p')
138 
139  print("Pol0 fit mean +- 68% percentile = ", round(mean,3), " +- ", round(stdev, 3))
140 
141  leg.SetBorderSize(0)
142  leg.SetTextSize(0.05)
143  leg.AddEntry(tg, pt.Leemumuratiolabel, "ep")
144  leg.AddEntry(tg.GetFunction("pol0"), "Mean = %.3f" % mean, "l")
145  leg.AddEntry(line1, "68%% band (#pm %.3f)" % stdev, "f")
146  leg.Draw()
147 
148  pt.drawAtlasLabel(xval, 0.88, "Internal")
149  pt.drawText(xval, 0.82, date_tag, size=labelsize)
150 
151  R.gPad.Update()
152 
153  tg.GetYaxis().SetRangeUser(ymin, ymax)
154  tg.GetXaxis().SetTitle(xtitle)
155  tg.GetXaxis().SetTimeDisplay(2)
156  tg.GetXaxis().SetNdivisions(9,R.kFALSE)
157  tg.GetXaxis().SetTimeFormat(time_format)
158  tg.GetXaxis().SetTimeOffset(0,"gmt")
159 
160  if years == ["22", "23", "24", "25"]:
161  plot_title = "Ratio of Electron and Muon channel Z-counting Luminosities across Run 3"
162  else:
163  plot_title = "Ratio of Electron and Muon channel Z-counting Luminosities across 20" + years[0]
164 
165  tg.SetTitle(plot_title)
166  c1.Update()
167  c1.SaveAs(outdir + "ZeeZmm_ratio_vs_time_"+out_tag+".pdf")
168 
169 
170 def zcounting_vs_atlas(channel, years):
171  """
172  Plot normalised comparison of Z-counting luminosity to ATLAS luminosity.
173  This can be done as a function of time and pileup.
174  """
175 
176  zstring = pt.plotlabel[channel]+" counting"
177  ytitle = "L_{"+pt.plotlabel[channel]+"}/L_{ATLAS}"
178  if args.absolute:
179  leg_entry = "L_{"+pt.plotlabel[channel]+"}"
180  else:
181  leg_entry = "L_{"+pt.plotlabel[channel]+"}^{"+norm_type+"-normalised}/L_{ATLAS}"
182 
183  print("Channel", channel, "comparison to ATLAS vs time for years: ", years)
184 
185  arr_date = []
186  arr_olumi = []
187  arr_zlumi = []
188  arr_zerr = []
189  run_num = []
190  fill_num = []
191 
192  for year in years:
193 
194  grl = pt.get_grl(year)
195 
196  for run in grl:
197  livetime, zlumi, zerr, olumi, timestamp, dfz_small = pt.get_dfz(args.indir, year, run, channel)
198 
199  # Cut out short runs
200  if livetime < pt.runlivetimecut:
201  if livetime >= 0.: print(f"Skip Run {run} because of live time {livetime/60:.1f} min")
202  continue
203 
204  prelratio = zlumi/olumi
205  if prelratio < ymin or prelratio > ymax:
206  print("WARNING: Run", run, "has", channel, "/ATLAS ratio %.2f +- %.2f, outside of y-axis range" % (prelratio, zerr/olumi) )
207 
208  # If plotting vs. date simply fill the arrays here
209  arr_date.append(timestamp)
210  arr_olumi.append(olumi)
211  arr_zlumi.append(zlumi)
212  arr_zerr.append(zerr)
213  run_num.append(run)
214  fill = dfz_small['FillNum'].median()
215  fill_num.append(int(fill))
216 
217 
218  # for ROOT plotting we need Python arrays
219  arr_date = array('d', arr_date)
220  # convert lists to numpy arrays
221  arr_olumi = np.array(arr_olumi)
222  arr_zlumi = np.array(arr_zlumi)
223  arr_zerr = np.array(arr_zerr)
224  total_lumi = arr_olumi.sum()/1000000
225  total_lumi_string = "Official DQ "
226  if year == "25": total_lumi_string = "Preliminary DQ "
227  total_lumi_string += str(round(total_lumi, 1)) + " fb^{-1}"
228 
229  #-----------Normalisation------------
230 
231  # Calculate and apply overall normalisation
232  if args.absolute:
233  normalisation = 1.0
234  else:
235  normalisation = np.sum(arr_zlumi) / np.sum(arr_olumi)
236 
237  # do normalisation to period integral
238  arr_zlumi /= normalisation
239  arr_zerr /= normalisation
240 
241  # calculate ratio to ATLAS preferred lumi
242  arr_zlumi_ratio = arr_zlumi/arr_olumi
243  arr_zerr_ratio = arr_zerr/arr_olumi
244 
245 #-----------Normalisation------------
246 
247  tg = R.TGraphErrors(len(arr_date), arr_date, array('d',arr_zlumi_ratio), R.nullptr, array('d',arr_zerr_ratio))
248 
249  # Depending if we're plotting over whole Run-3, change canvas size
250  if multiyear:
251  c1 = R.TCanvas("c1", "c1", 2000, 1000)
252  leg = R.TLegend(0.5, 0.18, 0.65, 0.4)
253  else:
254  c1 = R.TCanvas()
255  leg = R.TLegend(0.6, 0.18, 0.75, 0.4)
256 
257  tg.Draw('ap e0')
258  tg.GetYaxis().SetRangeUser(ymin, ymax)
259  if args.absolute:
260  plot_title = "Absolute L_{"+ zstring +"} to ATLAS across " + norm_type
261  else:
262  plot_title = "Normalised L_{"+ zstring +"} to ATLAS across " + norm_type
263 
264  # Plot 68% percentile band
265  stdev = np.percentile(abs(arr_zlumi_ratio - np.median(arr_zlumi_ratio)), 68)
266  print("68% band =", stdev)
267  tg.Fit('pol0', '0q')
268  mean = tg.GetFunction('pol0').GetParameter(0)
269  print("const of pol0 fit", mean)
270  print("median", np.median(arr_zlumi_ratio))
271  print("mean", np.mean(arr_zlumi_ratio))
272 
273  line1 = pt.make_bands(arr_date, stdev, np.median(arr_zlumi_ratio))
274  line1.Draw("same 3")
275  tg.Draw('same ep')
276 
277  leg.SetFillStyle(0)
278  leg.SetBorderSize(0)
279  leg.SetTextSize(0.05)
280  leg.AddEntry(tg, leg_entry, "ep")
281  leg.AddEntry(line1, "68%% band (#pm %.3f)" % stdev, "f")
282  leg.Draw()
283 
284  pt.drawAtlasLabel(xval, yval-0.47, "Internal")
285  pt.drawText(xval, yval-0.53, date_tag, size=labelsize)
286  pt.drawText(xval, yval-0.59, zstring, size=labelsize)
287  pt.drawText(xval, yval-0.65, "OflLumi-Run3-006", size=labelsize)
288  pt.drawText(xval, yval-0.04, total_lumi_string, size=labelsize)
289 
290  pt.drawText(xval, 0.88, plot_title, size=labelsize)
291 
292  tg.GetXaxis().SetTitle(xtitle)
293  tg.GetYaxis().SetTitle(ytitle)
294  tg.GetYaxis().SetRangeUser(ymin, ymax)
295  tg.GetXaxis().SetTimeDisplay(2)
296  tg.GetXaxis().SetNdivisions(9,R.kFALSE)
297  tg.GetXaxis().SetTimeFormat(time_format)
298  tg.GetXaxis().SetTimeOffset(0,"gmt")
299 
300  c1.Update()
301  c1.Modified()
302 
303  filename = outdir + channel + "ATLAS_ratio_vs_time_"+out_tag
304  if args.absolute:
305  c1.SaveAs(filename+"_abs.pdf")
306  else:
307  c1.SaveAs(filename+".pdf")
308 
309  if outcsv:
310  if args.absolute:
311  csvfile = open(filename+"_abs.csv", 'w')
312  else:
313  csvfile = open(filename+".csv", 'w')
314  csvwriter = csv.writer(csvfile, delimiter=',')
315  csvwriter.writerow(['FillNum','RunNum','Time','OffLumi','ZLumi','ZLumiErr','OffZlumi','OffZlumiErr'])
316  for i in range(len(run_num)):
317  csvwriter.writerow([fill_num[i], run_num[i], arr_date[i], arr_olumi[i], arr_zlumi[i], arr_zerr[i], arr_zlumi_ratio[i], arr_zerr_ratio[i]])
318  csvfile.close()
319 
320 
321 
322 if __name__ == "__main__":
323  pt.setAtlasStyle()
324  R.gROOT.SetBatch(R.kTRUE)
325  main()
DerivationFramework::TriggerMatchingUtils::sorted
std::vector< typename R::value_type > sorted(const R &r, PROJ proj={})
Helper function to create a sorted vector from an unsorted range.
MuonGM::round
float round(const float toRound, const unsigned int decimals)
Definition: Mdt.cxx:27
InDet::median
float median(std::vector< float > &Vec)
Definition: BTagVrtSec.cxx:35
plotting.yearwise_luminosity.zcounting_vs_atlas
def zcounting_vs_atlas(channel, years)
Definition: yearwise_luminosity.py:170
plotting.yearwise_luminosity.channel_comparison
def channel_comparison(years)
Definition: yearwise_luminosity.py:74
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
print
void print(char *figname, TCanvas *c1)
Definition: TRTCalib_StrawStatusPlots.cxx:25
array
Trk::open
@ open
Definition: BinningType.h:40
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
str
Definition: BTagTrackIpAccessor.cxx:11
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
plotting.yearwise_luminosity.main
def main()
Definition: yearwise_luminosity.py:68
Generated on Mon Sep 1 2025 21:24:21 for ATLAS Offline Software by   doxygen 1.8.18