Belle II Software  release-05-01-25
analyzeGainEff.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3 
4 # ---------------------------------------------------------------------------------------
5 # TOP MCP-PMT channel-by-channel gain/efficiency anslysis using laser data
6 # Usage: basf2 analyzeGainEff.py <input_file.sroot> --arg --slot (slotNum) --arg --PMT (PMTNum)
7 # ---------------------------------------------------------------------------------------
8 
9 from basf2 import *
10 import sys
11 import os
12 import argparse
13 import re
14 
15 parser = argparse.ArgumentParser(description="analyze gain and efficiency for laser run data")
16 parser.add_argument("inputFile", nargs='*', default=["NoInputFile"],
17  help="Input sroot files. Mulitple files can be given.")
18 parser.add_argument("--interimRootFile", default="NoInterimRootFile",
19  help="Interim root file name to store timing-height 2D histograms.")
20 parser.add_argument("--outputRootFile", default="NoOutputRootFile",
21  help="Output root file name to save TTree containing fit results.")
22 parser.add_argument("--outputPDFFile", default="NoOutputPDFFile",
23  help="Output PDF file name to save fitting results for each channel.")
24 parser.add_argument("--interimFW", action="store_true", default=False,
25  help="use when analyzing interimFW data")
26 parser.add_argument("--slotID", type=int, default=0,
27  help="slot number [1-16]")
28 parser.add_argument("--PMTID", type=int, default=0,
29  help="slot number [1-32]")
30 parser.add_argument("--calChannel", type=int, default=0,
31  help="asic channel number where calibration signals are injected [0-7]")
32 parser.add_argument("--threshold", type=float, default=100.,
33  help="threshold in pulse height fitting")
34 parser.add_argument("--globalDAQ", action="store_true", default=False,
35  help="Force to assume global DAQ data.")
36 parser.add_argument("--pocketDAQ", action="store_true", default=False,
37  help="Force to assume global DAQ data.")
38 parser.add_argument("--noOfflineFE", action="store_true", default=False,
39  help="Use only online FE hits, and do not use waveform data.")
40 parser.add_argument("--useSingleCalPulse", action="store_true", default=False,
41  help="Do not require double calibration pulses, but require only the first one.")
42 parser.add_argument("--windowSelection", type=int, default=0,
43  help="Select window number (All=0,Odd=2,Even=1)")
44 parser.add_argument("--includePrimaryChargeShare", action="store_true", default=False,
45  help="don't want to use primary charge share cut,type true.")
46 parser.add_argument("--includeAllChargeShare", action="store_true", default=False,
47  help="don't want to use all charge share cut,type true.")
48 parser.add_argument("--sumChargeShare", action="store_true", default=False,
49  help="When you want to use summation of charge share for gain distribution,type true.")
50 parser.add_argument("--fitoption", action="store", default='L',
51  help="fitting with chisquare=R, loglikelihood=L")
52 parser.add_argument("--timeCut", type=float, default=1,
53  help="cut of hittiming for chargeshare")
54 parser.add_argument("--HVHigherSetting", type=int, default=0,
55  help="HV positive difference from nominal value")
56 parser.add_argument("--HVLowerSetting", type=int, default=0,
57  help="HV negative difference from nominal value")
58 args = parser.parse_args()
59 
60 if (args.inputFile[0] == "NoInputFile") and (args.interimRootFile == "NoInterimRootFile"):
61  print("Steering file to study gain/efficiency analysis from laser data.")
62  print("In the first step, 2D histograms (hit time vs pulse height) were created for all the available channels, "
63  "and saved in an interim root file.")
64  print("Then, in the second process, 1D pulse height distribution is extracted "
65  "from the 2D histogram and fitted to evaluate gain and efficiency for each channel.")
66  print("The second process is done only for one give PMT as it takes time.")
67  print("usage:")
68  print("basf2 analyzeGainEff.py [input_filename1.sroot, input_filename2.sroot, ...]")
69  print(" [--arg --interimRootFile interim_histo_output.root]")
70  print(" [--arg --outputRootFile summary_tree_output.root]")
71  print(" [--arg --outputPDFFile summary_plot.pdf]")
72  print(" [--arg --slotID slotNum(1-16)] [--arg --PMTID PMTID(1-32)]")
73  print(" [--arg --calChannel asicCh(0-7)] [--arg --threshold threshold]")
74  print(" [--arg --globalDAQ] [--arg --pocketDAQ]")
75  print(" [--arg --noOfflineFE] [--arg --useSingleCalPulse]")
76  print(" [--arg --timeCut] [--arg --fitoption]")
77  print(" [--arg --windowSelection] [--arg --sumChargeShare]")
78  print(" [--arg --includePrimaryChargeShare] [--arg --includeAllChargeShare]")
79  print(" [--arg --HVHigherSetting] [--arg --HVLowerSetting]")
80  print("*When input sroot files are missing but interim root file is given,"
81  " skip the first process to create 2D histogram and start fitting.")
82  print("*Both the slot and PMT numbers are mandatory to proceed to the second processes.")
83  print("*default calibration asic channel is " + str(args.calChannel))
84  print("*default threshold is " + str(args.threshold))
85  print()
86  sys.exit()
87 
88 inputFiles = args.inputFile
89 interimRoot = args.interimRootFile
90 outputRoot = args.outputRootFile
91 outputPDF = args.outputPDFFile
92 isInterimFW = args.interimFW
93 slotId = args.slotID
94 pmtId = args.PMTID
95 calChannel = args.calChannel
96 threshold = args.threshold
97 isGlobalDAQ = False
98 isGlobalDAQForced = args.globalDAQ
99 isPocketDAQForced = args.pocketDAQ
100 isOfflineFEDisabled = args.noOfflineFE
101 windowSelect = args.windowSelection
102 includePrimaryChargeShare = args.includePrimaryChargeShare
103 includeAllChargeShare = args.includeAllChargeShare
104 sumChargeShare = args.sumChargeShare
105 HVHigherSetting = args.HVHigherSetting
106 HVLowerSetting = args.HVLowerSetting
107 fitOption = args.fitoption
108 timeCut = args.timeCut
109 
110 useSingleCalPulse = (True if isOfflineFEDisabled else args.useSingleCalPulse)
111 skipFirst = ((inputFiles[0] == "NoInputFile") and (interimRoot != "NoInterimRootFile"))
112 skipSecond = ((slotId < 1) or (slotId > 16) or (pmtId < 1) or (pmtId > 32))
113 pmtStr = "s" + ('%02d' % slotId) + "_PMT" + ('%02d' % pmtId)
114 if isGlobalDAQForced and isPocketDAQForced:
115  print("ERROR : both of --GlobalDAQ or --PocketDAQ can not be given.")
116  sys.exit()
117 if skipFirst and skipSecond:
118  print("ERROR : unable to run either of the first and second process:"
119  "no input sroot file and no valid slot ID or PMT ID")
120  sys.exit()
121 if (calChannel < 0) or (calChannel > 7):
122  print("ERROR : invalid asic channel with calibration pulses : " + str(calChannel))
123  sys.exit()
124 
125 if includePrimaryChargeShare and includeAllChargeShare:
126  print("ERROR : both of --includePrimaryChargeShare and --includeAllChargeShare can not be given.")
127  sys.exit()
128 
129 if sumChargeShare and includeAllChargeShare:
130  print("ERROR : both of --sumChargeShare and --includeAllChargeShare can not be given."
131  "While both of --sumChargeShare and --includePrimaryChargeShare can be given.")
132  sys.exit()
133 
134 if HVHigherSetting and HVLowerSetting:
135  print("ERROR : both of --HVHigherSetting and --HVLowerSetting can not be given.")
136  sys.exit()
137 
138 # data base
139 reset_database()
140 path_to_db = "/group/belle2/group/detector/TOP/calibration/combined/Combined_TBCrun417x_LocaT0run4855_AfterRelease01/localDB/"
141 use_local_database(path_to_db + '/localDB.txt', path_to_db)
142 
143 inputBase = inputFiles[0] if not skipFirst else interimRoot
144 dotPos = inputBase.rfind('.')
145 outputBase = inputBase[0:dotPos] if (dotPos > 0) else inputBase
146 
147 if re.search(r"run[0-9]+_slot[0-1][0-9]", inputFiles[0]):
148  outputBase = re.search(r"run[0-9]+_slot[0-1][0-9]", inputFiles[0]).group()
149 elif re.search(r"(top|cosmic|cdc|ecl|klm|test|debug|beam)\.[0-9]+\.[0-9]+", inputFiles[0]):
150  isGlobalDAQ = True
151  outputBase = re.search(r"(top|cosmic|cdc|ecl|klm|test|debug|beam)\.[0-9]+\.[0-9]+", inputFiles[0]).group()
152 
153 if interimRoot is "NoInterimRootFile":
154  interimRoot = outputBase + "_gain_histo.root"
155 if outputRoot is "NoOutputRootFile":
156  outputRoot = outputBase + "_gain_" + pmtStr + ".root"
157 if outputPDF is "NoOutputPDFFile":
158  outputPDF = outputBase + "_" + pmtStr
159 
160 if isGlobalDAQForced and (not isPocketDAQForced):
161  isGlobalDAQ = True
162 elif (not isGlobalDAQForced) and isPocketDAQForced:
163  isGlobalDAQ = False
164 
165 if not skipFirst:
166  print("*first process : " + str(inputFiles) + " --> " + interimRoot)
167 else:
168  print("*first process is skipped...")
169 if not skipSecond:
170  print("*second process : " + interimRoot + " --> " + outputRoot + ", " + outputPDF)
171 else:
172  print("*second process is skipped")
173 print("*Is global DAQ? : " + str(isGlobalDAQ))
174 print("*Offline FE : " + ("enabled" if not isOfflineFEDisabled else "disabled"))
175 print("*use double pulse : " + str(not useSingleCalPulse))
176 print("*cal. asic channel: " + str(calChannel))
177 print("*threshold : " + str(threshold))
178 print()
179 print("start process...")
180 
181 
182 if not skipFirst:
183  # Define a global tag
184  use_central_database('data_reprocessing_proc8')
185 
186  # Create path
187  first = create_path()
188 
189  srootInput = register_module('SeqRootInput')
190  srootInput.param('inputFileNames', inputFiles)
191  first.add_module(srootInput)
192 
193  # HistoManager
194  histoman = register_module('HistoManager')
195  histoman.param('histoFileName', interimRoot)
196  first.add_module(histoman)
197 
198  # conversion from RawCOPPER or RawDataBlock to RawDetector objects
199  if not isGlobalDAQ:
200  converter = register_module('Convert2RawDet')
201  first.add_module(converter)
202 
203  # geometry parameters
204  gearbox = register_module('Gearbox')
205  first.add_module(gearbox)
206 
207  # Geometry (only TOP needed)
208  geometry = register_module('Geometry')
209  first.add_module(geometry)
210 
211  # Unpacking
212  unpack = register_module('TOPUnpacker')
213  if isInterimFW:
214  unpack.param('swapBytes', True)
215  unpack.param('dataFormat', 0x0301)
216  first.add_module(unpack)
217 
218  # Add multiple hits by running feature extraction offline
219  if not isOfflineFEDisabled:
220  featureExtractor = register_module('TOPWaveformFeatureExtractor')
221  first.add_module(featureExtractor)
222 
223  # Convert to TOPDigits
224  converter = register_module('TOPRawDigitConverter')
225  converter.param('useSampleTimeCalibration', False)
226  converter.param('useChannelT0Calibration', False)
227  converter.param('useModuleT0Calibration', False)
228  converter.param('useCommonT0Calibration', False)
229  converter.param('calibrationChannel', calChannel) # if set, cal pulses will be flagged
230  converter.param('calpulseHeightMin', 100) # in [ADC counts]
231  converter.param('calpulseHeightMax', 650) # in [ADC counts]
232  converter.param('calpulseWidthMin', 0.8) # in [ns]
233  converter.param('calpulseWidthMax', 2.4) # in [ns]
234  first.add_module(converter)
235 
236  flagSetter = register_module('TOPXTalkChargeShareSetter')
237  flagSetter.param('sumChargeShare', sumChargeShare)
238  flagSetter.param('timeCut', timeCut) # in [nsec]
239  first.add_module(flagSetter)
240 
241  laserHitSelector = register_module('TOPLaserHitSelector')
242  laserHitSelector.param('useDoublePulse', (not useSingleCalPulse))
243  laserHitSelector.param('minHeightFirstCalPulse', 300) # in [ADC counts]
244  laserHitSelector.param('minHeightSecondCalPulse', 100) # in [ADC counts]
245  laserHitSelector.param('nominalDeltaT', 25.5) # in [ns]
246  laserHitSelector.param('nominalDeltaTRange', 2) # in [ns]
247  laserHitSelector.param('windowSelect', windowSelect)
248  laserHitSelector.param('includePrimaryChargeShare', includePrimaryChargeShare)
249  laserHitSelector.param('includeAllChargeShare', includeAllChargeShare)
250  # laserHitSelector.param('timeHistogramBinning', [100,-150,-50]) # number of bins, lower limit, upper limit
251  first.add_module(laserHitSelector)
252 
253  # Print progress
254  progress = register_module('Progress')
255  first.add_module(progress)
256 
257  # Process events
258  process(first)
259 
260 if not skipSecond:
261  # second process
262  print("start the second process...")
263  second = create_path()
264 
265  eventinfosetter = register_module('EventInfoSetter')
266  eventinfosetter.param('evtNumList', [1])
267  second.add_module(eventinfosetter)
268 
269  # HistoManager
270  histoman2 = register_module('HistoManager')
271  histoman2.param('histoFileName', outputRoot)
272  second.add_module(histoman2)
273 
274  analysis = register_module('TOPGainEfficiencyCalculator')
275  analysis.param('inputFile', interimRoot)
276  analysis.param('outputPDFFile', outputPDF)
277  analysis.param('targetSlotId', slotId)
278  analysis.param('targetPmtId', pmtId)
279  analysis.param('hvDiff', HVHigherSetting - HVLowerSetting)
280  analysis.param('threshold', threshold)
281  analysis.param('fitoption', fitOption)
282  second.add_module(analysis)
283 
284  process(second)
285 
286 # Print statistics
287 print(statistics)