File indexing completed on 2024-11-25 02:29:58
0001
0002 import argparse
0003 import RecoLuminosity.LumiDB.LumiConstants as LumiConstants
0004 import re
0005 from math import sqrt
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 if __name__ == '__main__':
0024 parameters = LumiConstants.ParametersObject()
0025
0026 parser = argparse.ArgumentParser(description="Script to estimate average and RMS pileup using the per-bunch luminosity information provided by brilcalc. The output is a JSON file containing a dictionary by runs with one entry per lumi section.")
0027 parser.add_argument('inputFile', help='CSV input file as produced from brilcalc')
0028 parser.add_argument('outputFile', help='Name of JSON output file')
0029 parser.add_argument('-b', '--selBX', metavar='BXLIST', help='Comma-separated list of BXs to use (will use all by default)')
0030 parser.add_argument('-n', '--no-threshold', action='store_true', help='By default, to avoid including spurious luminosity from afterglow in the pileup calculation, bunches with luminosity below a given threshold are excluded. This threshold is 8.0/ub/LS for HFOC at nominal energy, 2.0/ub/LS for HFOC at low energy, and 1.2/ub/LS for other luminometers. If the input data has already been preprocessed (e.g. by using the --xingTr argument in brilcalc) to exclude these bunches, or if you are running on special data with low overall luminosity, then use this flag to disable the application of the threshold.')
0031 args = parser.parse_args()
0032
0033 output = args.outputFile
0034
0035 sel_bx = set()
0036 if args.selBX:
0037 for ibx in args.selBX.split(","):
0038 try:
0039 bx = int(ibx)
0040 sel_bx.add(bx)
0041 except:
0042 print(ibx,"is not an int")
0043 print("Processing",args.inputFile,"with selected BXs:", sorted(sel_bx))
0044 else:
0045 print("Processing",args.inputFile,"with all BX")
0046
0047
0048
0049 sepRE = re.compile(r'[\]\[\s,;:]+')
0050 csv_input = open(args.inputFile, 'r')
0051 last_run = -1
0052
0053 last_valid_lumi = []
0054
0055 output_line = '{'
0056 for line in csv_input:
0057 if line[0] == '#':
0058 continue
0059
0060 pieces = sepRE.split(line.strip())
0061
0062 if len(pieces) < 15:
0063
0064
0065 raise RuntimeError("Not enough fields in input line; maybe you forgot to include --xing in your brilcalc command?\n"+line)
0066 try:
0067 run = int(pieces[0])
0068 lumi_section = int(pieces[2])
0069 beam_energy = float(pieces[10])
0070
0071
0072 luminometer = pieces[14]
0073
0074 if luminometer == "HFOC":
0075 if beam_energy > 3000:
0076
0077 threshold = 8.0
0078 else:
0079
0080 threshold = 2.0
0081 else:
0082 threshold = 1.2
0083
0084 xing_lumi_array = []
0085 for bxid, bunch_del_lumi, bunch_rec_lumi in zip(pieces[15::3], pieces[16::3], pieces[17::3]):
0086 if sel_bx and int(bxid) not in sel_bx:
0087 continue
0088 if args.no_threshold or float(bunch_del_lumi) > threshold:
0089 xing_lumi_array.append([int(bxid), float(bunch_del_lumi), float(bunch_rec_lumi)])
0090 except:
0091 print("Failed to parse line: check if the input format has changed")
0092 print(pieces[0],pieces[1],pieces[2],pieces[3],pieces[4],pieces[5],pieces[6],pieces[7],pieces[8],pieces[9])
0093 continue
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103 if run != last_run:
0104
0105 if last_run > 0:
0106 output_line = output_line[:-1] + '], '
0107 last_run = run
0108 output_line += ('\n"%d":' % run )
0109 output_line += ' ['
0110
0111
0112 total_lumi = 0
0113 total_int = 0
0114 total_int2 = 0
0115 total_weight2 = 0
0116
0117
0118 for bxid, bunch_del_lumi, bunch_rec_lumi in xing_lumi_array:
0119 total_lumi += bunch_rec_lumi
0120
0121 total_int += bunch_del_lumi*bunch_rec_lumi
0122
0123
0124
0125 total_int *= parameters.orbitLength / parameters.lumiSectionLength
0126 if total_lumi > 0:
0127 mean_int = total_int/total_lumi
0128 else:
0129 mean_int = 0
0130
0131
0132 for bxid, bunch_del_lumi, bunch_rec_lumi in xing_lumi_array:
0133 mean_pileup = bunch_del_lumi * parameters.orbitLength / parameters.lumiSectionLength
0134 if mean_pileup > 100:
0135 print("mean number of pileup events > 100 for run %d, lum %d : m %f l %f" % \
0136 (run, lumi_section, mean_pileup, bunch_del_lumi))
0137
0138
0139 total_int2 += bunch_rec_lumi*(mean_pileup-mean_int)*(mean_pileup-mean_int)
0140 total_weight2 += bunch_rec_lumi*bunch_rec_lumi
0141
0142
0143
0144 bunch_rms_lumi = 0
0145 denom = total_lumi*total_lumi-total_weight2
0146 if total_lumi > 0 and denom > 0:
0147 bunch_rms_lumi = sqrt(total_lumi*total_int2/denom)
0148
0149 output_line += "[%d,%2.4e,%2.4e,%2.4e]," % (lumi_section, total_lumi, bunch_rms_lumi, mean_int)
0150 last_valid_lumi = [lumi_section, total_lumi, bunch_rms_lumi, mean_int]
0151
0152 output_line = output_line[:-1] + ']}'
0153 csv_input.close()
0154
0155 outputfile = open(output,'w')
0156 if not outputfile:
0157 raise RuntimeError("Could not open '%s' as an output JSON file" % output)
0158
0159 outputfile.write(output_line)
0160 outputfile.close()
0161 print("Output written to", output)