CMS 3D CMS Logo

create_public_peakpu_plots.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 ######################################################################
4 ## File: create_public_lumi_plots.py
5 ######################################################################
6 
7 from __future__ import print_function
8 import sys
9 import csv
10 import os
11 import commands
12 import time
13 import datetime
14 import calendar
15 import copy
16 import math
17 import optparse
18 import ConfigParser
19 
20 import numpy as np
21 
22 import six
23 import matplotlib
24 matplotlib.use('Agg')
25 from matplotlib import pyplot as plt
26 # FIX FIX FIX
27 # This fixes a well-know bug with stepfilled logarithmic histograms in
28 # Matplotlib.
29 from RecoLuminosity.LumiDB.mpl_axes_hist_fix import hist
30 if matplotlib.__version__ != '1.0.1':
31  print("ERROR The %s script contains a hard-coded bug-fix " \
32  "for Matplotlib 1.0.1. The Matplotlib version loaded " \
33  "is %s" % (__file__, matplotlib.__version__), file=sys.stderr)
34  sys.exit(1)
35 matplotlib.axes.Axes.hist = hist
36 # FIX FIX FIX end
37 
38 from RecoLuminosity.LumiDB.public_plots_tools import ColorScheme
39 from RecoLuminosity.LumiDB.public_plots_tools import LatexifyUnits
40 from RecoLuminosity.LumiDB.public_plots_tools import AddLogo
41 from RecoLuminosity.LumiDB.public_plots_tools import InitMatplotlib
42 from RecoLuminosity.LumiDB.public_plots_tools import SavePlot
43 from RecoLuminosity.LumiDB.public_plots_tools import FONT_PROPS_SUPTITLE
44 from RecoLuminosity.LumiDB.public_plots_tools import FONT_PROPS_TITLE
45 from RecoLuminosity.LumiDB.public_plots_tools import FONT_PROPS_AX_TITLE
46 from RecoLuminosity.LumiDB.public_plots_tools import FONT_PROPS_TICK_LABEL
47 
48 try:
49  import debug_hook
50  import pdb
51 except ImportError:
52  pass
53 
54 
55 ######################################################################
56 
57 # Some global constants. Not nice, but okay.
58 DATE_FMT_STR_LUMICALC = "%m/%d/%y %H:%M:%S"
59 DATE_FMT_STR_LUMICALC_DAY = "%m/%d/%y"
60 DATE_FMT_STR_OUT = "%Y-%m-%d %H:%M"
61 DATE_FMT_STR_AXES = "%-d %b"
62 DATE_FMT_STR_CFG = "%Y-%m-%d"
63 NUM_SEC_IN_LS = 2**18 / 11246.
64 
65 KNOWN_ACCEL_MODES = ["PROTPHYS", "IONPHYS", "PAPHYS",
66  "2013_amode_bug_workaround"]
67 LEAD_SCALE_FACTOR = 82. / 208.
68 
69 def processdata(years):
70  results={}#{year:[[str,str],[pu,pu]]}
71  for y in years:
72  fills=[]#[fillnum]
73  maxputimelist=[]
74  maxpulist=[]
75  results[y]=[maxputimelist,maxpulist]
76  f=None
77  try:
78  lumifilename=str(y)+'lumibyls.csv'
79  f=open(lumifilename,'rb')
80  except IOError:
81  print('failed to open file ',lumifilename)
82  return result
83  freader=csv.reader(f,delimiter=',')
84  idx=0
85  for row in freader:
86  if idx==0:
87  idx=1
88  continue
89  if row[0].find('#')==1:
90  continue
91  [runnum,fillnum]=map(lambda i:int(i),row[0].split(':'))
92  avgpu=float(row[7])
93  if avgpu>50: continue
94  putime=row[2]
95  max_avgpu=avgpu
96  max_putime=putime
97  if fillnum not in fills:#new fill
98  fills.append(fillnum)
99  results[y][0].append(max_putime)
100  results[y][1].append(max_avgpu)
101  if avgpu>max_avgpu:
102  results[y][0][-1]=max_putime
103  results[y][1][-1]=max_avgpu
104  print(results)
105  return results
106 
107 
108 ######################################################################
109 
110 def CacheFilePath(cache_file_dir, day=None):
111  cache_file_path = os.path.abspath(cache_file_dir)
112  if day:
113  cache_file_name = "lumicalc_cache_%s.csv" % day.isoformat()
114  cache_file_path = os.path.join(cache_file_path, cache_file_name)
115  return cache_file_path
116 
117 
118 def GetXLocator(ax):
119  """Pick a DateLocator based on the range of the x-axis."""
120  (x_lo, x_hi) = ax.get_xlim()
121  num_days = x_hi - x_lo
122  min_num_ticks = min(num_days, 5)
123  locator = matplotlib.dates.AutoDateLocator(minticks=min_num_ticks,
124  maxticks=None)
125  # End of GetLocator().
126  return locator
127 
128 ######################################################################
129 
130 def TweakPlot(fig, ax, time_range,
131  add_extra_head_room=False):
132 
133  # Fiddle with axes ranges etc.
134  (time_begin, time_end) = time_range
135  ax.relim()
136  ax.autoscale_view(False, True, True)
137  for label in ax.get_xticklabels():
138  label.set_ha("right")
139  label.set_rotation(30.)
140 
141  # Bit of magic here: increase vertical scale by one tick to make
142  # room for the legend.
143  if add_extra_head_room:
144  y_ticks = ax.get_yticks()
145  (y_min, y_max) = ax.get_ylim()
146  is_log = (ax.get_yscale() == "log")
147  y_max_new = y_max
148  if is_log:
149  tmp = y_ticks[-1] / y_ticks[-2]
150  y_max_new = y_max * math.pow(tmp, add_extra_head_room)
151  else:
152  tmp = y_ticks[-1] - y_ticks[-2]
153  y_max_new = y_max + add_extra_head_room * tmp
154  ax.set_ylim(y_min, y_max_new)
155 
156  # Add a second vertical axis on the right-hand side.
157  ax_sec = ax.twinx()
158  ax_sec.set_ylim(ax.get_ylim())
159  ax_sec.set_yscale(ax.get_yscale())
160 
161  for ax_tmp in fig.axes:
162  for sub_ax in [ax_tmp.xaxis, ax_tmp.yaxis]:
163  for label in sub_ax.get_ticklabels():
164  label.set_font_properties(FONT_PROPS_TICK_LABEL)
165 
166  ax.set_xlim(time_begin, time_end)
167 
168  locator = GetXLocator(ax)
169  minorXlocator=matplotlib.ticker.AutoMinorLocator()
170  #ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
171  #ax.yaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(5))
172  ax.xaxis.set_major_locator(locator)
173  ax.xaxis.set_minor_locator(minorXlocator)
174 
175  formatter = matplotlib.dates.DateFormatter(DATE_FMT_STR_AXES)
176  ax.xaxis.set_major_formatter(formatter)
177 
178  fig.subplots_adjust(top=.85, bottom=.14, left=.13, right=.91)
179  # End of TweakPlot().
180 
181 ######################################################################
182 
183 if __name__ == "__main__":
184 
185  desc_str = "test pu" \
186 
187  arg_parser = optparse.OptionParser(description=desc_str)
188  arg_parser.add_option("--ignore-cache", action="store_true",
189  help="Ignore all cached lumiCalc results " \
190  "and re-query lumiCalc. " \
191  "(Rebuilds the cache as well.)")
192  (options, args) = arg_parser.parse_args()
193  if len(args) != 1:
194  print("ERROR Need exactly one argument: a config file name", file=sys.stderr)
195  sys.exit(1)
196  config_file_name = args[0]
197  ignore_cache = options.ignore_cache
198 
199  cfg_defaults = {
200  "lumicalc_flags" : "",
201  "date_end" : None,
202  "color_schemes" : "Joe, Greg",
203  "beam_energy" : None,
204  "beam_fluctuation" : None,
205  "verbose" : False,
206  "oracle_connection" : None
207  }
208  cfg_parser = ConfigParser.SafeConfigParser(cfg_defaults)
209  if not os.path.exists(config_file_name):
210  print("ERROR Config file '%s' does not exist" % config_file_name, file=sys.stderr)
211  sys.exit(1)
212  cfg_parser.read(config_file_name)
213 
214  # Which color scheme to use for drawing the plots.
215  color_scheme_names_tmp = cfg_parser.get("general", "color_schemes")
216  color_scheme_names = [i.strip() for i in color_scheme_names_tmp.split(",")]
217  # Where to store cache files containing the lumiCalc output.
218  cache_file_dir = cfg_parser.get("general", "cache_dir")
219  # Flag to turn on verbose output.
220  verbose = cfg_parser.getboolean("general", "verbose")
221 
222  # Some details on how to invoke lumiCalc.
223  lumicalc_script = cfg_parser.get("general", "lumicalc_script")
224  lumicalc_flags_from_cfg = cfg_parser.get("general", "lumicalc_flags")
225  accel_mode = cfg_parser.get("general", "accel_mode")
226  # Check if we know about this accelerator mode.
227  if not accel_mode in KNOWN_ACCEL_MODES:
228  print("ERROR Unknown accelerator mode '%s'" % \
229  accel_mode, file=sys.stderr)
230 
231  # WORKAROUND WORKAROUND WORKAROUND
232  amodetag_bug_workaround = False
233  if accel_mode == "2013_amode_bug_workaround":
234  amodetag_bug_workaround = True
235  accel_mode = "PAPHYS"
236  # WORKAROUND WORKAROUND WORKAROUND end
237 
238  beam_energy_tmp = cfg_parser.get("general", "beam_energy")
239  # If no beam energy specified, use the default(s) for this
240  # accelerator mode.
241  beam_energy = None
242  beam_energy_from_cfg = None
243  if not beam_energy_tmp:
244  print("No beam energy specified --> using defaults for '%s'" % \
245  accel_mode)
246  beam_energy_from_cfg = False
247  else:
248  beam_energy_from_cfg = True
249  beam_energy = float(beam_energy_tmp)
250 
251  beam_fluctuation_tmp = cfg_parser.get("general", "beam_fluctuation")
252  # If no beam energy fluctuation specified, use the default for
253  # this accelerator mode.
254  beam_fluctuation = None
255  beam_fluctuation_from_cfg = None
256  if not beam_fluctuation_tmp:
257  print("No beam energy fluctuation specified --> using the defaults to '%s'" % \
258  accel_mode)
259  beam_fluctuation_from_cfg = False
260  else:
261  beam_fluctuation_from_cfg = True
262  beam_fluctuation = float(beam_fluctuation_tmp)
263 
264  # Overall begin and end dates of all data to include.
265  tmp = cfg_parser.get("general", "date_begin")
266  date_begin = datetime.datetime.strptime(tmp, DATE_FMT_STR_CFG).date()
267  tmp = cfg_parser.get("general", "date_end")
268  date_end = None
269  if tmp:
270  date_end = datetime.datetime.strptime(tmp, DATE_FMT_STR_CFG).date()
271  # If no end date is given, use today.
272  today = datetime.datetime.utcnow().date()
273  if not date_end:
274  print("No end date given --> using today")
275  date_end = today
276  # If end date lies in the future, truncate at today.
277  if date_end > today:
278  print("End date lies in the future --> using today instead")
279  date_end = today
280  # If end date is before start date, give up.
281  if date_end < date_begin:
282  print("ERROR End date before begin date (%s < %s)" % \
283  (date_end.isoformat(), date_begin.isoformat()), file=sys.stderr)
284  sys.exit(1)
285 
286  # If an Oracle connection string is specified, use direct Oracle
287  # access. Otherwise access passes through the Frontier
288  # cache. (Fine, but much slower to receive the data.)
289  oracle_connection_string = cfg_parser.get("general", "oracle_connection")
290  use_oracle = (len(oracle_connection_string) != 0)
291 
292  ##########
293 
294  # Map accelerator modes (as fed to lumiCalc) to particle type
295  # strings to be used in plot titles etc.
296  particle_type_strings = {
297  "PROTPHYS" : "pp",
298  "IONPHYS" : "PbPb",
299  "PAPHYS" : "pPb"
300  }
301  particle_type_str = particle_type_strings[accel_mode]
302 
303  beam_energy_defaults = {
304  "PROTPHYS" : {2010 : 3500.,
305  2011 : 3500.,
306  2012 : 4000.,
307  2013 : 1380.1},
308  "IONPHYS" : {2010 : 3500.,
309  2011 : 3500.},
310  "PAPHYS" : {2013 : 4000.}
311  }
312  beam_fluctuation_defaults = {
313  "PROTPHYS" : {2010 : .15,
314  2011 : .15,
315  2012 : .15,
316  2013 : .15},
317  "IONPHYS" : {2010 : .15,
318  2011 : .15},
319  "PAPHYS" : {2013 : .15}
320  }
321 
322  ##########
323 
324  # Environment parameter for access to the Oracle DB.
325  if use_oracle:
326  os.putenv("TNS_ADMIN", "/afs/cern.ch/cms/lumi/DB")
327 
328  ##########
329 
330  # Tell the user what's going to happen.
331  print("Using configuration from file '%s'" % config_file_name)
332  if ignore_cache:
333  print("Ignoring all cached lumiCalc results (and rebuilding the cache)")
334  else:
335  print("Using cached lumiCalc results from %s" % \
336  CacheFilePath(cache_file_dir))
337  print("Using color schemes '%s'" % ", ".join(color_scheme_names))
338  print("Using lumiCalc script '%s'" % lumicalc_script)
339  print("Using additional lumiCalc flags from configuration: '%s'" % \
340  lumicalc_flags_from_cfg)
341  print("Selecting data for accelerator mode '%s'" % accel_mode)
342  if beam_energy_from_cfg:
343  print("Selecting data for beam energy %.0f GeV" % beam_energy)
344  else:
345  print("Selecting data for default beam energy for '%s' from:" % accel_mode)
346  for (key, val) in six.iteritems(beam_energy_defaults[accel_mode]):
347  print(" %d : %.1f GeV" % (key, val))
348  if beam_fluctuation_from_cfg:
349  print("Using beam energy fluctuation of +/- %.0f%%" % \
350  (100. * beam_fluctuation))
351  else:
352  print("Using default beam energy fluctuation for '%s' from:" % accel_mode)
353  for (key, val) in six.iteritems(beam_fluctuation_defaults[accel_mode]):
354  print(" %d : +/- %.0f%%" % (key, 100. * val))
355  if use_oracle:
356  print("Using direct access to the Oracle luminosity database")
357  else:
358  print("Using access to the luminosity database through the Frontier cache")
359 
360  ##########
361 
362  # See if the cache file dir exists, otherwise try to create it.
363  path_name = CacheFilePath(cache_file_dir)
364  if not os.path.exists(path_name):
365  if verbose:
366  print("Cache file path does not exist: creating it")
367  try:
368  os.makedirs(path_name)
369  except Exception as err:
370  print("ERROR Could not create cache dir: %s" % path_name, file=sys.stderr)
371  sys.exit(1)
372 
373  ##########
374 
375 
376  years=[2010,2011,2012]
377 
378  lumi_data_by_fill_per_year={} #{year:[[timestamp,timestamp],[max_pu,max_pu]]}
379  lumi_data_by_fill_per_year=processdata(years)
380 
381  #lumi_data_by_fill_per_year[2010]=[['05/05/10 05:59:58','06/02/10 10:47:25','07/02/10 12:47:25','08/02/10 11:47:25'],[10.0,2.5,11.3,4.5]]
382  #lumi_data_by_fill_per_year[2011]=[['05/05/11 05:59:58','06/02/11 10:47:25','07/02/11 12:47:25','08/02/11 11:47:25','09/05/11 05:59:58','10/02/11 10:47:25','11/02/11 12:47:25'],[20.0,27.4,30.5,40.,22.,15.,45.]]
383  #lumi_data_by_fill_per_year[2012]=[['05/05/12 05:59:58','06/02/12 10:47:25','07/02/12 12:47:25','08/02/12 11:47:25','09/05/12 05:59:58','10/02/12 10:47:25','11/02/12 12:47:25','11/02/12 12:47:25'],[10.0,17.4,20.5,30.,32.,25.,33.,42.]]
384 
385 
387 
388  ##########
389 
390  year_begin = date_begin.isocalendar()[0]
391  year_end = date_end.isocalendar()[0]
392  # DEBUG DEBUG DEBUG
393  assert year_end >= year_begin
394  ##########
395 
396  # And this is where the plotting starts.
397  print("Drawing things...")
398  ColorScheme.InitColors()
399 
400  #----------
401 
402  if len(years) > 1:
403  print(" peak interactions for %s together" % ", ".join([str(i) for i in years]))
404 
405  def PlotPeakPUAllYears(lumi_data_by_fill_per_year):
406  """Mode 1: years side-by-side"""
407 
408  # Loop over all color schemes and plot.
409  for color_scheme_name in color_scheme_names:
410  print(" color scheme '%s'" % color_scheme_name)
411  color_scheme = ColorScheme(color_scheme_name)
412  color_by_year = color_scheme.color_by_year
413  logo_name = color_scheme.logo_name
414  file_suffix = color_scheme.file_suffix
415 
416  for type in ["lin", "log"]:
417  is_log = (type == "log")
418  aspect_ratio = matplotlib.figure.figaspect(1. / 2.5)
419  fig = plt.figure(figsize=aspect_ratio)
420  ax = fig.add_subplot(111)
421 
422  time_begin_ultimate = datetime.datetime.strptime(lumi_data_by_fill_per_year[years[0]][0][0],DATE_FMT_STR_LUMICALC).date()
423  str_begin_ultimate = time_begin_ultimate.strftime(DATE_FMT_STR_OUT)
424  for (year_index, year) in enumerate(years):
425 
426  lumi_data = lumi_data_by_fill_per_year[year]
427  times_tmp = [datetime.datetime.strptime(tmp,DATE_FMT_STR_LUMICALC).date() for tmp in lumi_data[0]]
428  times = [matplotlib.dates.date2num(i) for i in times_tmp] #x_axis
429  maxpus = lumi_data[1] # y_axis
430 
431  # NOTE: Special case for 2010.
432  label = None
433  if year == 2010 or year == 2011 :
434  label = r"%d, %s" % \
435  (year,r'7TeV $\sigma$=71.5mb')
436  else:
437  label = r"%d, %s" % \
438  (year,r'8TeV $\sigma$=73mb')
439  ax.plot(times,maxpus,
440  color=color_by_year[year],
441 # marker="none", linestyle="solid",
442  marker="o", linestyle='none',
443  linewidth=4,
444  label=label)
445  if is_log:
446  ax.set_yscale("log")
447 
448  time_begin = datetime.datetime(years[0], 1, 1, 0, 0, 0)
449  time_end = datetime.datetime(years[-1], 12, 16, 20, 50,9)
450  str_begin = time_begin.strftime(DATE_FMT_STR_OUT)
451  str_end = time_end.strftime(DATE_FMT_STR_OUT)
452 
453  num_cols = None
454  num_cols = len(years)
455  tmp_x = 0.095
456  tmp_y = .95
457 
458  leg = ax.legend(loc="upper left", bbox_to_anchor=(tmp_x, 0., 1., tmp_y),
459  frameon=False, ncol=num_cols)
460  for t in leg.get_texts():
461  t.set_font_properties(FONT_PROPS_TICK_LABEL)
462 
463  # Set titles and labels.
464  fig.suptitle(r"CMS peak interactions per crossing, %s" % particle_type_str,
465  fontproperties=FONT_PROPS_SUPTITLE)
466  ax.set_title("Data included from %s to %s UTC \n" % \
467  (str_begin_ultimate, str_end),
468  fontproperties=FONT_PROPS_TITLE)
469  ax.set_xlabel(r"Date (UTC)", fontproperties=FONT_PROPS_AX_TITLE)
470  ax.set_ylabel(r"Peak interactions per crossing",\
471  fontproperties=FONT_PROPS_AX_TITLE)
472 
473  # Add the logo.
474  #zoom = 1.7
475  zoom = .95
476  AddLogo(logo_name, ax, zoom=zoom)
477  extra_head_room = 0
478  if is_log:
479  #if mode == 1:
480  extra_head_room = 1
481  #elif mode == 2:
482  # extra_head_room = 2
483  TweakPlot(fig, ax, (time_begin, time_end),
484 # TweakPlot(fig, ax, (time_begin_ultimate, time_end),
485  add_extra_head_room=True)
486 
487  log_suffix = ""
488  if is_log:
489  log_suffix = "_log"
490  SavePlot(fig, "peak_pu_%s_%s%s" % \
491  (particle_type_str.lower(),
492  log_suffix, file_suffix))
493  PlotPeakPUAllYears(lumi_data_by_fill_per_year)
494 
495  plt.close()
496 
497  print("Done")
498 
499 ######################################################################
def TweakPlot(fig, ax, time_range, add_extra_head_room=False)
S & print(S &os, JobReport::InputFile const &f)
Definition: JobReport.cc:66
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
T min(T a, T b)
Definition: MathUtil.h:58
static std::string join(char **cmd)
Definition: RemoteFile.cc:18
def AddLogo(logo_name, ax, zoom=1.2)
#define str(s)
def PlotPeakPUAllYears(lumi_data_by_fill_per_year)
def SavePlot(fig, file_name_base)
double split
Definition: MVATrainer.cc:139
def CacheFilePath(cache_file_dir, day=None)