-
Notifications
You must be signed in to change notification settings - Fork 2
/
create_public_lumi_plots.py
executable file
·2661 lines (2293 loc) · 115 KB
/
create_public_lumi_plots.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python
# coding: utf-8
######################################################################
## File: create_public_lumi_plots.py
######################################################################
from __future__ import print_function
import sys
import os
import commands
import time
import datetime
import calendar
import copy
import math
import optparse
import ConfigParser
# CS : This funky module is in CMSSW but we want to be independent of
# CMSSW also if this costs some performance. We go for the
# standard json decoder...
# Introduced plt.close(fig) calls to release memory that was getting
# a bit tight...
# Adjusted the minimum of the log y axis to the value calculated and
# and previously used in the matplotlib fix.
# import cjson
import json
import numpy as np
import six
import matplotlib
# CS : To be compliant with the code below:
cjson = json.JSONDecoder()
matplotlib.use('Agg')
from matplotlib import pyplot as plt
# CS : This fix was for the old matplotlib version in CMSSW. We now
# go for the brilconda version which today (26/1/2022) features
# matplotlib 1.5.3. The sources reveal that the bug has been fixed.
# FIX FIX FIX
#
# This fixes a well-know bug with stepfilled logarithmic histograms in
# Matplotlib.
# from mpl_axes_hist_fix import hist
# if matplotlib.__version__ != '1.0.1':
# print("ERROR The %s script contains a hard-coded bug-fix " \
# "for Matplotlib 1.0.1. The Matplotlib version loaded " \
# "is %s" % (__file__, matplotlib.__version__), file=sys.stderr)
# sys.exit(1)
# matplotlib.axes.Axes.hist = hist
#
# FIX FIX FIX end
try:
from brokenaxes import brokenaxes
has_brokenaxes = True
except ImportError:
print("Warning: brokenaxes package is not installed, cannot produce multi-year cumulative plots with cutouts")
print("Use pip install --user brokenaxes to install")
has_brokenaxes = False
from public_plots_tools import ColorScheme
from public_plots_tools import LatexifyUnits
from public_plots_tools import AddLogo
from public_plots_tools import InitMatplotlib
from public_plots_tools import SavePlot
from public_plots_tools import FONT_PROPS_SUPTITLE
from public_plots_tools import FONT_PROPS_TITLE
from public_plots_tools import FONT_PROPS_AX_TITLE
from public_plots_tools import FONT_PROPS_TICK_LABEL
from public_plots_tools import FONT_PROPS_PLOT_LABEL
try:
import debug_hook
import pdb
except ImportError:
pass
######################################################################
# Some global constants. Not nice, but okay.
DATE_FMT_STR_LUMICALC = "%m/%d/%y %H:%M:%S"
DATE_FMT_STR_LUMICALC_DAY = "%m/%d/%y"
DATE_FMT_STR_OUT = "%Y-%m-%d %H:%M"
DATE_FMT_STR_AXES = "%-d %b"
DATE_FMT_STR_CFG = "%Y-%m-%d"
NUM_SEC_IN_LS = 2**18 / 11246.
KNOWN_ACCEL_MODES = ["PROTPHYS", "IONPHYS", "PAPHYS", "ALLIONS",
"2013_amode_bug_workaround"]
LEAD_SCALE_FACTOR = 82. / 208.
######################################################################
class LumiDataPoint(object):
"""Holds info from one line of lumiCalc lumibyls output."""
def __init__(self, line, json_file_name=None, data_scale_factor=1):
# Decode the comma-separated line from lumiCalc.
line_split = line.split(",")
tmp = line_split[0].split(":")
self.run_number = int(tmp[0])
self.fill_number = int(tmp[1])
tmp = line_split[1].split(":")
self.ls = int(tmp[0])
tmp = line_split[2]
self.timestamp = datetime.datetime.strptime(tmp, DATE_FMT_STR_LUMICALC)
# The 1e6 is to convert from ub^{-1} to b^{-1}, and any other factor can be applied here as well.
scale_factor = 1.e6*data_scale_factor
self.lum_del = scale_factor * float(line_split[5])
self.lum_rec = scale_factor * float(line_split[6])
#IR add avgpu
self.avgpu = scale_factor * float(line_split[7])
# Adding lum_cert for the data certification information
if json_file_name:
addcertls = bool(checkCertification(self.run_number, self.ls))
if addcertls:
self.lum_cert = scale_factor * float(line_split[6])
self.avgpu_cert = scale_factor * float(line_split[7])
else:
self.lum_cert = 0.
self.avgpu_cert = 0.
else:
self.lum_cert = 0.
self.avgpu_cert = 0.
#IR add avgpu
# End of __init__().
# End of class LumiDataPoint.
######################################################################
class LumiDataBlock(object):
"""A supposedly coherent block of LumiDataPoints.
NOTE: No checks on duplicates, sorting, etc.
"""
scale_factors = {
"fb^{-1}" : 1.e-15,
"pb^{-1}" : 1.e-12,
"nb^{-1}" : 1.e-9,
"ub^{-1}" : 1.e-6,
"mb^{-1}" : 1.e-3,
"b^{-1}" : 1.,
"Hz/fb" : 1.e-15,
"Hz/pb" : 1.e-12,
"Hz/nb" : 1.e-9,
"Hz/ub" : 1.e-6,
"Hz/mb" : 1.e-3,
"Hz/b" : 1.
}
def __init__(self, data_point=None):
if not data_point:
self.data_points = []
else:
self.data_points = [data_point]
# End of __init__().
def __iadd__(self, other):
self.data_points.extend(other.data_points)
# End of __iadd__().
return self
def __lt__(self, other):
# End of __lt__().
return self.time_mid() < other.time_mid()
def add(self, new_point):
self.data_points.append(new_point)
# End of add().
def len(self):
# End of add().
return len(self.data_points)
def copy(self):
# End of copy().
return copy.deepcopy(self)
def is_empty(self):
# End of is_empty().
return not len(self.data_points)
def lum_del_tot(self, units="b^{-1}"):
res = sum([i.lum_del for i in self.data_points])
res *= LumiDataBlock.scale_factors[units]
# End of lum_del_tot().
return res
def lum_rec_tot(self, units="b^{-1}"):
res = sum([i.lum_rec for i in self.data_points])
res *= LumiDataBlock.scale_factors[units]
# End of lum_rec_tot().
return res
def lum_cert_tot(self, units="b^{-1}"):
res = sum([i.lum_cert for i in self.data_points])
res *= LumiDataBlock.scale_factors[units]
# End of lum_cert_tot().
return res
def max_inst_lum(self, units="Hz/b"):
res = 0.
if len(self.data_points):
res = max([i.lum_del for i in self.data_points])
res /= NUM_SEC_IN_LS
res *= LumiDataBlock.scale_factors[units]
# End of max_inst_lum().
return res
def straighten(self):
self.data_points.sort()
# End of straighten().
def time_begin(self):
res = min([i.timestamp for i in self.data_points])
# End of time_begin().
return res
def time_end(self):
res = max([i.timestamp for i in self.data_points])
# End of time_end().
return res
def time_mid(self):
delta = self.time_end() - self.time_begin()
delta_sec = delta.days * 24 * 60 * 60 + delta.seconds
res = self.time_begin() + datetime.timedelta(seconds=.5*delta_sec)
# End of time_mid().
return res
# End of class LumiDataBlock.
######################################################################
class LumiDataBlockCollection(object):
"""A collection of LumiDataBlocks."""
def __init__(self, data_block=None):
if not data_block:
self.data_blocks = []
else:
self.data_blocks = [data_block]
# End of __init__().
def __len__(self):
# End of __len__().
return len(self.data_blocks)
def add(self, new_block):
self.data_blocks.append(new_block)
# End of add().
def sort(self):
self.data_blocks.sort()
# End of sort().
def time_begin(self):
res = datetime.datetime.max
if len(self.data_blocks):
res = min([i.time_begin() for i in self.data_blocks])
# End of time_begin().
return res
def time_end(self):
res = datetime.datetime.min
if len(self.data_blocks):
res = max([i.time_end() for i in self.data_blocks])
# End of time_end().
return res
def times(self):
res = [i.time_mid() for i in self.data_blocks]
# End of times().
return res
def lum_del(self, units="b^{-1}"):
res = [i.lum_del_tot(units) for i in self.data_blocks]
# End of lum_del().
return res
def lum_rec(self, units="b^{-1}"):
res = [i.lum_rec_tot(units) for i in self.data_blocks]
# End of lum_rec().
return res
def lum_cert(self, units="b^{-1}"):
res = [i.lum_cert_tot(units) for i in self.data_blocks]
# End of lum_cert().
return res
def lum_del_tot(self, units="b^{-1}"):
# End of lum_del().
return sum(self.lum_del(units))
def lum_rec_tot(self, units="b^{-1}"):
# End of lum_rec().
return sum(self.lum_rec(units))
def lum_cert_tot(self, units="b^{-1}"):
# End of lum_cert().
return sum(self.lum_cert(units))
def lum_inst_max(self, units="Hz/b"):
res = [i.max_inst_lum(units) for i in self.data_blocks]
# End of lum_inst_max().
return res
# End of class LumiDataBlockCollection.
######################################################################
def CacheFilePath(cache_file_dir, day=None):
cache_file_path = os.path.abspath(cache_file_dir)
if day:
cache_file_name = "lumicalc_cache_%s.csv" % day.isoformat()
cache_file_path = os.path.join(cache_file_path, cache_file_name)
return cache_file_path
######################################################################
def AtMidnight(datetime_in):
res = datetime.datetime.combine(datetime_in.date(), datetime.time())
# End of AtMidnight().
return res
######################################################################
def AtMidWeek(datetime_in):
# NOTE: The middle of the week is on Thursday according to our
# definition
tmp = datetime_in.date()
date_tmp = tmp - \
datetime.timedelta(days=tmp.weekday()) + \
datetime.timedelta(days=3)
res = datetime.datetime.combine(date_tmp, datetime.time())
# End of AtMidWeek().
return res
######################################################################
def GetUnits(year, accel_mode, mode):
# First check to see if units were specified in the config file.
# If so, use them!
if (cfg_parser.get("general", "units")):
cfg_units = cjson.decode(cfg_parser.get("general", "units"))
if mode in cfg_units:
return cfg_units[mode]
units_spec = {
"PROTPHYS" : {
2010 : {
"cum_day" : "pb^{-1}",
"cum_week" : "pb^{-1}",
"cum_year" : "pb^{-1}",
"max_inst" : "Hz/ub",
},
2011 : {
"cum_day" : "pb^{-1}",
"cum_week" : "pb^{-1}",
"cum_year" : "fb^{-1}",
"max_inst" : "Hz/nb",
},
2012 : {
"cum_day" : "pb^{-1}",
"cum_week" : "fb^{-1}",
"cum_year" : "fb^{-1}",
"max_inst" : "Hz/nb",
},
2013 : {
"cum_day" : "pb^{-1}",
"cum_week" : "pb^{-1}",
"cum_year" : "pb^{-1}",
"max_inst" : "Hz/ub",
},
2015 : {
"cum_day" : "fb^{-1}",
"cum_week" : "fb^{-1}",
"cum_year" : "fb^{-1}",
"max_inst" : "Hz/nb",
},
2016 : {
"cum_day" : "pb^{-1}",
"cum_week" : "fb^{-1}",
"cum_year" : "fb^{-1}",
"max_inst" : "Hz/nb",
},
2017 : {
"cum_day" : "pb^{-1}",
"cum_week" : "fb^{-1}",
"cum_year" : "fb^{-1}",
"max_inst" : "Hz/nb",
},
2018 : {
"cum_day" : "pb^{-1}",
"cum_week" : "fb^{-1}",
"cum_year" : "fb^{-1}",
"max_inst" : "Hz/nb",
}
},
"IONPHYS" : {
2011 : {
"cum_day" : "ub^{-1}",
"cum_week" : "ub^{-1}",
"cum_year" : "ub^{-1}",
"max_inst" : "Hz/mb",
},
2015 : { #brilcalc is coming in giving barn instead of microbarn
"cum_day" : "ub^{-1}",
"cum_week" : "ub^{-1}",
"cum_year" : "ub^{-1}",
"max_inst" : "Hz/mb",
},
2018 : {
"cum_day" : "ub^{-1}",
"cum_week" : "ub^{-1}",
"cum_year" : "ub^{-1}",
"max_inst" : "Hz/mb",
}
},
"PAPHYS" : {
2013 : {
"cum_day" : "nb^{-1}",
"cum_week" : "nb^{-1}",
"cum_year" : "nb^{-1}",
"max_inst" : "Hz/mb",
},
2016 : {
"cum_day" : "nb^{-1}",
"cum_week" : "nb^{-1}",
"cum_year" : "nb^{-1}",
"max_inst" : "Hz/mb",
}
}
}
units = None
try:
units = units_spec[accel_mode][year][mode]
except KeyError:
if mode == "cum_day":
units = "pb^{-1}"
elif mode == "cum_week":
units = "pb^{-1}"
elif mode == "cum_year":
units = "fb^{-1}"
elif mode == "max_inst":
units = "Hz/ub"
# DEBUG DEBUG DEBUG
assert not units is None
# DEBUG DEBUG DEBUG end
# End of GetUnits().
return units
######################################################################
def GetEnergyPerNucleonScaleFactor(amodetag):
assert amodetag in ["IONPHYS", "PAPHYS"]
res = LEAD_SCALE_FACTOR
if amodetag == "PAPHYS":
res = math.sqrt(res)
# End of GetEnergyPerNucleonScaleFactor().
return res
######################################################################
# Given the beam energy in GeV, produce a string version of the CMS energy in TeV,
# including the pPb or PbPb scale factor if necessary.
def FormatCMSEnergy(beam_energy, accel_mode, year, include_units=True):
cms_energy = 2 * beam_energy
cms_energy_str = "???"
if accel_mode == "PROTPHYS":
width = 0
if year == 2013:
width = 2
if year >= 2021:
# We do not run on an integer enery anymore...
width = 1
cms_energy_str = "%.*f" % \
(width, 1.e-3 * cms_energy)
if include_units:
cms_energy_str += " TeV"
elif accel_mode in ["IONPHYS", "PAPHYS", "ALLIONS"]:
this_accel_mode = accel_mode
if accel_mode == "ALLIONS":
# Get the actual mode used for this year from the config file.
try:
accel_mode_by_year = cjson.decode(cfg_parser.get("general", "accel_mode_by_year"))
this_accel_mode = accel_mode_by_year[str(year)]
except:
print("Error: accelerator mode for year",year,"not specified in accel_mode_by_year parameter in config file")
return cms_energy_str
cms_energy_str = "%.2f" % \
(1.e-3 * GetEnergyPerNucleonScaleFactor(this_accel_mode) * cms_energy)
if include_units:
if accel_mode == "ALLIONS":
cms_energy_str = particle_type_strings[this_accel_mode] + " " + cms_energy_str
cms_energy_str += " TeV/nucleon"
return cms_energy_str
######################################################################
def NumDaysInYear(year):
"""Returns the number of days in the given year."""
date_lo = datetime.date(year, 1, 1)
date_hi = datetime.date(year + 1, 1, 1)
num_days = (date_hi - date_lo).days
# End of NumDaysInYear().
return num_days
######################################################################
def GetXLocator(ax):
"""Pick a DateLocator based on the range of the x-axis."""
(x_lo, x_hi) = ax.get_xlim()
num_days = x_hi - x_lo
min_num_ticks = min(num_days, 5)
max_num_ticks = min(num_days, 20)
locator = matplotlib.dates.AutoDateLocator(minticks=min_num_ticks,
maxticks=max_num_ticks)
# End of GetLocator().
return locator
######################################################################
def TweakPlot(fig, ax, time_range,
add_extra_head_room=False, headroom_factor=1.3):
# Fiddle with axes ranges etc.
(time_begin, time_end) = time_range
ax.relim()
ax.autoscale_view(False, True, True)
for label in ax.get_xticklabels():
label.set_ha("right")
label.set_rotation(30.)
# Bit of magic here: increase vertical scale by one tick to make
# room for the legend.
if add_extra_head_room:
y_ticks = ax.get_yticks()
(y_min, y_max) = ax.get_ylim()
is_log = (ax.get_yscale() == "log")
y_max_new = y_max
if is_log:
tmp = y_ticks[-1] / y_ticks[-2]
y_max_new = y_max * math.pow(tmp, (add_extra_head_room + 2))
else:
tmp = y_ticks[-1] - y_ticks[-2]
y_max_new = y_max + add_extra_head_room * (headroom_factor*tmp)
ax.set_ylim(y_min, y_max_new)
## Add a second vertical axis on the right-hand side.
# CS: Even though useful, not CMS standard therefore commented out.
#ax_sec = ax.twinx()
#ax_sec.set_ylim(ax.get_ylim())
#ax_sec.set_yscale(ax.get_yscale())
for ax_tmp in fig.axes:
for sub_ax in [ax_tmp.xaxis, ax_tmp.yaxis]:
for label in sub_ax.get_ticklabels():
label.set_font_properties(FONT_PROPS_TICK_LABEL)
ax.set_xlim(time_begin, time_end)
locator = GetXLocator(ax)
ax.xaxis.set_major_locator(locator)
formatter = matplotlib.dates.DateFormatter(DATE_FMT_STR_AXES)
ax.xaxis.set_major_formatter(formatter)
fig.subplots_adjust(top=.85, bottom=.14, left=.13, right=.95)
# End of TweakPlot().
######################################################################
def checkCertification(run_number, ls):
"""Check if this run and LS are certified as good and return a boolean parameter."""
try:
ls_ranges = certification_data[run_number]
for ranges in ls_ranges:
if (ls >= ranges[0]) and (ls <= ranges[1]):
return True
except KeyError:
return False
return False
######################################################################
def loadCertificationJSON(json_file_name):
full_file = open(json_file_name, "r")
#full_file_content = ["".join(l.rstrip()) for l in full_file.readlines()]
#full_object = cjson.decode(full_file_content[0])
full_file_content=full_file.read().replace("\n","")
#print str(full_file_content)
full_object = cjson.decode(full_file_content)
# Now turn this into a dictionary for easier handling.
tmp = full_object.keys()
tmp = [int(i) for i in tmp]
run_list = sorted(tmp)
certification_data = {}
for run in run_list:
ls_ranges = full_object.get(str(run), None)
certification_data[run] = ls_ranges
return certification_data
######################################################################
if __name__ == "__main__":
desc_str = "This script creates the official CMS luminosity plots " \
"based on the output from the lumiCalc family of scripts."
arg_parser = optparse.OptionParser(description=desc_str)
arg_parser.add_option("--ignore-cache", action="store_true",
help="Ignore all cached brilcalc results " \
"and re-query brilcalc. " \
"(Rebuilds the cache as well.)")
(options, args) = arg_parser.parse_args()
if len(args) != 1:
print("ERROR Need exactly one argument: a config file name", file=sys.stderr)
sys.exit(1)
config_file_name = args[0]
ignore_cache = options.ignore_cache
cfg_defaults = {
"lumicalc_flags": "",
"date_end": None,
"color_schemes": "Joe, Greg",
"beam_energy": None,
"verbose": False,
"oracle_connection": "",
"json_file": None,
"file_suffix": "",
"plot_label": " ",
"units": None,
"display_units": None,
"normtag_file": None,
"display_scale_factor": None,
"data_scale_factor": None,
"skip_years": None,
"plot_multiple_years": "False",
"filter_brilcalc_results": "True",
}
cfg_parser = ConfigParser.SafeConfigParser(cfg_defaults)
if not os.path.exists(config_file_name):
print("ERROR Config file '%s' does not exist" % config_file_name, file=sys.stderr)
sys.exit(1)
cfg_parser.read(config_file_name)
print("Using configuration from file '%s'" % config_file_name)
# See if we're running in single-year or multiple-year mode. See the README for more
# details on how this works.
plot_multiple_years = cfg_parser.getboolean("general", "plot_multiple_years")
print("plot multiple years mode:",plot_multiple_years)
filter_brilcalc_results = cfg_parser.getboolean("general", "filter_brilcalc_results")
# Which color scheme to use for drawing the plots.
color_scheme_names_tmp = cfg_parser.get("general", "color_schemes")
color_scheme_names = [i.strip() for i in color_scheme_names_tmp.split(",")]
# Where to store cache files containing the brilcalc output.
cache_file_dir = cfg_parser.get("general", "cache_dir")
# Flag to turn on verbose output.
verbose = cfg_parser.getboolean("general", "verbose")
# Suffix to append to all file names.
file_suffix2 = cfg_parser.get("general", "file_suffix")
# Some details on how to invoke brilcalc.
lumicalc_script = cfg_parser.get("general", "lumicalc_script")
# Don't let people try to use lumiCalc2.py or pixelLumiCalc.py or lcr2.py or
# really anything other than brilcalc -- sadness will ensue.
if not "brilcalc" in lumicalc_script:
print("ERROR: Lumi calculation scripts other than brilcalc are no longer supported.", file=sys.stderr)
print("Please update your config file appropriately.", file=sys.stderr)
sys.exit(1)
lumicalc_flags_from_cfg = cfg_parser.get("general", "lumicalc_flags")
accel_mode = cfg_parser.get("general", "accel_mode")
# Check if we know about this accelerator mode.
if not accel_mode in KNOWN_ACCEL_MODES:
print("ERROR Unknown accelerator mode '%s'" % accel_mode, file=sys.stderr)
if accel_mode == "ALLIONS" and not plot_multiple_years:
print("Accelerator mode",accel_mode,"is only meaningful for multiple-year plots, sorry!", file=sys.stderr)
sys.exit(1)
# WORKAROUND WORKAROUND WORKAROUND
amodetag_bug_workaround = False
if accel_mode == "2013_amode_bug_workaround":
amodetag_bug_workaround = True
accel_mode = "PAPHYS"
# WORKAROUND WORKAROUND WORKAROUND end
beam_energy_tmp = cfg_parser.get("general", "beam_energy")
# If no beam energy specified, use the default(s) for this
# accelerator mode.
beam_energy = None
beam_energy_from_cfg = None
if not beam_energy_tmp:
print("No beam energy specified --> using defaults for '%s'" % \
accel_mode)
beam_energy_from_cfg = False
else:
beam_energy_from_cfg = True
beam_energy = float(beam_energy_tmp)
# get the directory where to put the plots
plot_directory_tmp = cfg_parser.get("general", "plot_directory")
if not plot_directory_tmp:
plot_directory = "plots"
print("No plot directory specified --> using default value '%s'" % plot_directory)
else:
plot_directory = plot_directory_tmp
print("Plots will be stored in directory '%s'." % plot_directory)
# Overall begin and end dates of all data to include.
tmp = cfg_parser.get("general", "date_begin")
date_begin = datetime.datetime.strptime(tmp, DATE_FMT_STR_CFG).date()
tmp = cfg_parser.get("general", "date_end")
date_end = None
if tmp:
date_end = datetime.datetime.strptime(tmp, DATE_FMT_STR_CFG).date()
# If no end date is given, use today.
today = datetime.datetime.utcnow().date()
if not date_end:
print("No end date given --> using today")
date_end = today
# If end date lies in the future, truncate at today.
if date_end > today and date_end.isocalendar()[0] != 2015 :
print("End date lies in the future --> using today instead")
date_end = today
# If end date is before start date, give up.
if date_end < date_begin:
print("ERROR End date before begin date (%s < %s)" % \
(date_end.isoformat(), date_begin.isoformat()), file=sys.stderr)
sys.exit(1)
# Years to skip if making a multiyear plot.
skip_years = []
if (cfg_parser.get("general", "skip_years")):
skip_years = cjson.decode(cfg_parser.get("general", "skip_years"))
# Oracle connection strings are no longer supported in brilcalc.
# (If you really want to specify a specific service use the
# -c option directly in the flags.)
if cfg_parser.get("general", "oracle_connection"):
print("WARNING: You have specified an Oracle connection string but these are no longer supported by brilcalc.", file=sys.stderr)
print("If you want to specify a particular database service, add the -c option to lumicalc_flags.", file=sys.stderr)
# Check if display units are specified. This is to work around the fact that in the PbPb runs
# everything is scaled by 1e6, so using the units as is will give wrong values. If display_units is set,
# then we do all of our calculations in the regular units but display the display units so everything
# works out correctly (albeit annoyingly).
display_units = None
if (cfg_parser.get("general", "display_units")):
display_units = cjson.decode(cfg_parser.get("general", "display_units"))
# See if the data needs to be scaled (e.g. for the 2015 PbPb data). This can be specified as either a float
# (in which case it will apply to everything) or as a dictionary (in which case you can specify different factors
# for different years, as in the multi-year PbPb plot).
data_scale_factor = None
if (cfg_parser.get("general", "data_scale_factor")):
data_scale_factor = cjson.decode(cfg_parser.get("general", "data_scale_factor"))
# For multiple-year plots, check to see if we should scale the data for particular
# years by a particular factor.
display_scale_factor = {}
if (cfg_parser.get("general", "display_scale_factor")):
display_scale_factor = cjson.decode(cfg_parser.get("general", "display_scale_factor"))
# Check to see if this is well-formed.
for year in display_scale_factor:
if type(display_scale_factor[year]) is not dict:
print("Error in display_scale_factor for "+year+": expected dictionary as entry", file=sys.stderr)
sys.exit(1)
if "integrated" not in display_scale_factor[year] or "peak" not in display_scale_factor[year]:
print("Error in display_scale_factor for "+year+": dictionary does not contain integrated and peak keys", file=sys.stderr)
sys.exit(1)
# If a JSON file is specified, use the JSON file to add in the
# plot data certified as good for physics.
json_file_name = cfg_parser.get("general", "json_file")
if len(str(json_file_name)) < 1:
json_file_name = None
if json_file_name:
if not os.path.exists(json_file_name):
print("ERROR Requested JSON file '%s' is not available" % json_file_name, file=sys.stderr)
sys.exit(1)
normtag_file = cfg_parser.get("general", "normtag_file")
##########
certification_data = None
if json_file_name:
certification_data = loadCertificationJSON(json_file_name)
##########
# Map accelerator modes (as fed to brilcalc) to particle type
# strings to be used in plot titles etc.
particle_type_strings = {
"PROTPHYS" : "pp",
"IONPHYS" : "PbPb",
"PAPHYS" : "pPb",
"ALLIONS": "PbPb+pPb"
}
particle_type_str = particle_type_strings[accel_mode]
beam_energy_defaults = {
"PROTPHYS" : {2010 : 3500.,
2011 : 3500.,
2012 : 4000.,
2013 : 1380.1,
2015 : 6500.,
2016 : 6500.,
2017 : 6500.,
2018 : 6500.,
2022 : 6800.,
2023 : 6800.},
"IONPHYS" : {2010 : 3500.,
2011 : 3500.,
2015 : 6369.,
2018 : 6370.},
"PAPHYS" : {2013 : 4000.,
2016 : 2500},
"ALLIONS": {2015 : 6369.,
2016 : 6500,
2018 : 6370.},
}
##########
# Tell the user what's going to happen.
print("Accelerator mode is '%s'" % accel_mode)
if ignore_cache:
print("Ignoring all cached brilcalc results (and rebuilding the cache)")
else:
print("Using cached brilcalc results from %s" % \
CacheFilePath(cache_file_dir))
# We only use brilcalc for the single-year plots, so don't bother printing out info if we're making a
# multi-year plot
if not plot_multiple_years:
print("Using brilcalc script '%s'" % lumicalc_script)
print("Using additional brilcalc flags from configuration: '%s'" % \
lumicalc_flags_from_cfg)
if normtag_file:
print("Normtag file selected:", normtag_file)
else:
print("No normtag file selected, online luminosity will be used")
if json_file_name:
print("Using JSON file '%s' for certified data" % json_file_name)
else:
print("No certification JSON file will be applied.")
if beam_energy_from_cfg:
print("Beam energy is %.0f GeV" % beam_energy)
else:
print("Using default beam energy for '%s' from:" % accel_mode)
for (key, val) in sorted(six.iteritems(beam_energy_defaults[accel_mode])):
print(" %d : %.1f GeV" % (key, val))
print("Using color schemes '%s'" % ", ".join(color_scheme_names))
##########
# See if the cache file dir exists, otherwise try to create it.
path_name = CacheFilePath(cache_file_dir)
if not os.path.exists(path_name):
if verbose:
print("Cache file path does not exist: creating it")
try:
os.makedirs(path_name)
except Exception as err:
print("ERROR Could not create cache dir: %s" % path_name, file=sys.stderr)
sys.exit(1)
##########
InitMatplotlib()
##########
week_begin = date_begin.isocalendar()[1]
week_end = date_end.isocalendar()[1]
year_begin = date_begin.isocalendar()[0]
year_end = date_end.isocalendar()[0]
# DEBUG DEBUG DEBUG
assert year_end >= year_begin
# DEBUG DEBUG DEBUG end
print("Building a list of days to include in the plots")
print(" first day to consider: %s (%d, week %d)" % \
(date_begin.isoformat(), year_begin, week_begin))
print(" last day to consider: %s (%d, week %d)" % \
(date_end.isoformat(), year_end, week_end))
num_days = (date_end - date_begin).days + 1
days = [date_begin + datetime.timedelta(days=i) for i in range(num_days)]
years = list(range(year_begin, year_end + 1))
weeks = []
day_cur = date_begin
while day_cur <= date_end:
year = day_cur.isocalendar()[0]
week = day_cur.isocalendar()[1]
weeks.append((year, week))
day_cur += datetime.timedelta(days=7)
if num_days <= 7:
year = date_end.isocalendar()[0]
week = date_end.isocalendar()[1]
weeks.append((year, week))
weeks = list(set(weeks))
weeks.sort()
# Figure out the last day we want to read back from the cache.
# NOTE: The above checking ensures that date_end is <= today, so
# the below only assumes that we're never more than N days
# behind on our luminosity numbers. For online we can use a smaller N,
# but if we're using a normtag file use a larger N to account for cases
# where we may have fallen behind on updating the fill validation.
day_margin = 3
if normtag_file:
day_margin = 7
last_day_from_cache = min(today - datetime.timedelta(days=day_margin), date_end)
if verbose:
print("Last day for which the cache will be used: %s" % \
last_day_from_cache.isoformat())
# First run brilcalc to get the luminosity for these days. This is only applicable
# for single-year plots; it's assumed that if you're making the multi-year plots
# you already have the luminosity data on hand.
if not plot_multiple_years:
print("Running brilcalc for all requested days")
for day in days:
print(" %s" % day.isoformat())
use_cache = (not ignore_cache) and (day <= last_day_from_cache)
cache_file_path = CacheFilePath(cache_file_dir, day)
cache_file_tmp = cache_file_path.replace(".csv", "_tmp.csv")
if (not os.path.exists(cache_file_path)) or (not use_cache):
date_begin_str = day.strftime(DATE_FMT_STR_LUMICALC)
date_begin_day_str = day.strftime(DATE_FMT_STR_LUMICALC_DAY)
date_end_str = (day + datetime.timedelta(days=1)).strftime(DATE_FMT_STR_LUMICALC)
date_previous_str = (day - datetime.timedelta(days=0)).strftime(DATE_FMT_STR_LUMICALC)
year = day.isocalendar()[0]
if year != 2014:
if not beam_energy_from_cfg:
beam_energy = beam_energy_defaults[accel_mode][year]
# WORKAROUND WORKAROUND WORKAROUND
# Trying to work around the issue with the unfilled
# accelerator mode in the RunInfo database.
if amodetag_bug_workaround:
# Don't use the amodetag in this case. Scary, but
# works for the moment.
lumicalc_flags = "%s " \
"--beamenergy %.0f "% \
(lumicalc_flags_from_cfg,
beam_energy)
else:
lumicalc_flags = lumicalc_flags_from_cfg
if filter_brilcalc_results:
lumicalc_flags += " --beamenergy %.0f --amodetag %s" % (beam_energy, accel_mode)
if normtag_file:
lumicalc_flags += " --normtag "+normtag_file
lumicalc_flags = lumicalc_flags.strip()
lumicalc_cmd = "%s %s" % (lumicalc_script, lumicalc_flags)
cmd = "%s --begin '%s' --end '%s' -o %s" % \
(lumicalc_cmd, date_previous_str, date_end_str, cache_file_tmp)
if verbose:
print(" running lumicalc as '%s'" % cmd)
(status, output) = commands.getstatusoutput(cmd)
# BUG BUG BUG
# Trying to track down the bad-cache problem.
output_0 = copy.deepcopy(output)
# BUG BUG BUG end
print(status)