File indexing completed on 2023-03-17 10:56:40
0001
0002
0003 from __future__ import print_function
0004 import sys
0005 import math
0006 from ROOT import *
0007 from copy import deepcopy
0008 from scipy import signal
0009
0010 gROOT.SetBatch()
0011
0012 class InefficientDeadROCs:
0013
0014
0015 def __TraverseDirTree(self, dir):
0016
0017 for obj in dir.GetListOfKeys():
0018 if not obj.IsFolder():
0019 if obj.ReadObjectAny(TClass.GetClass("TH2")):
0020 th1 = deepcopy(obj.ReadObj())
0021 name = th1.GetName()
0022 if name.startswith(self.lookForStr):
0023
0024 newName = name.split(self.lookForStr)[1]
0025 th1.SetName(newName)
0026
0027
0028 layer = 0
0029
0030 if newName.startswith("B"):
0031 layer = "B" + ((newName.split("_LYR"))[1])[0]
0032 else:
0033 layer = ((newName.split("_D"))[1])[0]
0034 if newName.startswith("FPix_Bm"):
0035 layer = "-" + layer
0036 layer = "F" + layer
0037
0038 if layer in self.dicOfModuleHistograms:
0039 self.dicOfModuleHistograms[layer].append(th1)
0040 else:
0041 self.dicOfModuleHistograms.update({layer : [th1]})
0042 else:
0043 self.__TraverseDirTree(obj.ReadObj())
0044
0045 def __init__(self, inputDQMName, outputFileName, noiseOutputFileName, dirs):
0046
0047 self.inputFileName = inputDQMName
0048 self.outputFileName = outputFileName
0049 self.noiseOutputFileName = noiseOutputFileName
0050 self.dirs = dirs
0051
0052 self.lookForStr = "digi_occupancy_per_col_per_row_"
0053
0054 self.rocMaxCol = 52
0055 self.rocMaxRow = 80
0056 self.rocsInRow = 8
0057 self.rocsInCol = 2
0058
0059 self.inputFile = TFile(self.inputFileName)
0060 self.dicOfModuleHistograms = {}
0061
0062
0063 self.pixelNoisynessTh = 6
0064 self.rocOccupancyTh = 200
0065
0066 self.barrelNoisyColumnTh = 1.35
0067 self.barrelNoisyColumnTh2 = 4.5
0068 self.endcapNoisyColumnTh = 1.5
0069
0070 self.barrelInefficientDColTh = 8
0071 self.endcapInefficientDColTh = 30
0072
0073
0074
0075 if self.inputFile.IsOpen():
0076 print("%s opened successfully!" % (self.inputFileName))
0077
0078 for dir in self.dirs:
0079 self.__TraverseDirTree(self.inputFile.Get(dir))
0080 print("Histograms to read: %d" % (len(self.dicOfModuleHistograms)))
0081
0082 self.detDict = {}
0083
0084 else:
0085 print("Unable to open file %s" % (self.inputFileName))
0086
0087 def __lmsExp(self, data, xMin, xMax):
0088 meanOfX = (xMax + xMin) * 0.5
0089 meanOfY = sum( [math.log(data[i]) for i in range(len(data))] ) / len(data)
0090
0091 D = 0
0092 for i in range(xMin, xMax + 1):
0093 D = D + (i - meanOfX)**2
0094
0095
0096 a = 0
0097 for i in range(len(data)):
0098 a = a + math.log(data[i]) * (xMin + i - meanOfX)
0099 a = a/D
0100
0101 lnb = meanOfY - a * meanOfX
0102
0103 return a, math.exp(lnb)
0104
0105 def __lmsLin(self, data, xMin, xMax):
0106 meanOfX = (xMax + xMin) * 0.5
0107 meanOfY = sum(data) / len(data)
0108
0109 D = 0
0110 for i in range(xMin, xMax + 1):
0111 D = D + (i - meanOfX)**2
0112
0113 a = 0
0114 for i in range(len(data)):
0115 a = a + data[i] * (xMin + i - meanOfX)
0116 a = a/D
0117
0118 b = meanOfY - a * meanOfX
0119
0120 return a, b, D
0121
0122 def __customMedianFilter(self, array, radius = 2):
0123
0124 filtered = [0 for i in range(len(array))]
0125 currArray = []
0126 for i in range(len(array)):
0127 if i - radius < 0:
0128 currArray = array[0 : i + radius + 1]
0129 elif i + radius + 1 >= len(array):
0130 currArray = array[i - radius : ]
0131
0132 currArray.sort()
0133 filtered[i] = currArray[len(currArray) // 2]
0134
0135 return filtered
0136
0137 def __getROCData(self, hist, startPixel, endPixel, row, repeatFilter = 3, filterKernelSize = 5):
0138 pixelArr = []
0139 columnsWithSuspiciouslyNoisyPixels = []
0140
0141 for x in range(startPixel, endPixel):
0142
0143 columnPixels = [hist.GetBinContent(x, y + 1) for y in range(row * self.rocMaxRow, (row + 1) * self.rocMaxRow)]
0144 columnSum = sum(columnPixels)
0145
0146 columnMean = columnSum / len(columnPixels)
0147 for i in range(len(columnPixels)):
0148 if columnPixels[i] > self.pixelNoisynessTh * columnMean:
0149
0150
0151 columnsWithSuspiciouslyNoisyPixels.append(x)
0152
0153
0154 break
0155
0156 pixelArr.append(columnSum)
0157
0158 if len(pixelArr) == 0:
0159 return None, None, None
0160
0161 medFiltRes, sciPyMedFiltRes = deepcopy(pixelArr), deepcopy(pixelArr)
0162 for i in range(repeatFilter):
0163 sciPyMedFiltRes = signal.medfilt(sciPyMedFiltRes, filterKernelSize)
0164 medFiltRes = self.__customMedianFilter(medFiltRes, filterKernelSize // 2)
0165
0166 return pixelArr, medFiltRes, columnsWithSuspiciouslyNoisyPixels, sciPyMedFiltRes
0167
0168 def __getPixelArrWithRemovedDrops(self, pixelArr, medFiltRes):
0169 return [ (pixelArr[i] if pixelArr[i] > medFiltRes[i] else medFiltRes[i]) if 0 < i < len(pixelArr) - 1 else min(medFiltRes) for i in range(len(pixelArr))]
0170
0171 def __normalizeArray(self, pixelArr):
0172 c_min, c_max = min(pixelArr), max(pixelArr)
0173 if c_min != c_max:
0174 c_diff_inv = 1.0 / (c_max - c_min)
0175 return [ (pixelArr[i] - c_min) * c_diff_inv for i in range(len(pixelArr))]
0176 return [0 for i in range(len(pixelArr))]
0177
0178 def __setNormalizedArrayZeroInThePoint(self, pixelArr, pt):
0179 c_diff_inv = 1.0 / (1.0 - pt)
0180
0181 return [ (pixelArr[i] - pt) * c_diff_inv for i in range(len(pixelArr))]
0182
0183 def __determineBarrelNoise(self, noiseFile, columnsWithSuspiciouslyNoisyPixels, histName, meanOfPixels, maxMed, val, pos, rocCol, rocRow):
0184 noisyROC = False;
0185 if meanOfPixels < self.rocOccupancyTh:
0186
0187 noisyROC = True
0188 else:
0189 th = self.barrelNoisyColumnTh * maxMed
0190 if val > th:
0191 if pos not in columnsWithSuspiciouslyNoisyPixels:
0192 rocNum, xCoordInROC = self.__convertCoordinatesFromHistToROCSpace(histName, pos, rocRow)
0193 noiseFile.write("%s\t(x, row)->[rocNum, xRoc]\t(%d, %d)->[%d, %d];\t{VAL, TH}\t{%f, %f}\n" % (histName, pos, rocRow+1, rocNum, xCoordInROC, val, th))
0194
0195 return 1, noisyROC
0196
0197
0198
0199 return 0, noisyROC
0200
0201 def __determineBarrelNoise2(self, noiseFile, columnsWithSuspiciouslyNoisyPixels, histName, meanOfPixels, normMeanOfPixels, normVal, pos, rocCol, rocRow):
0202 noisyROC = False;
0203 if meanOfPixels < self.rocOccupancyTh:
0204
0205 noisyROC = True
0206 else:
0207 th = self.barrelNoisyColumnTh2 * normMeanOfPixels
0208 if normVal > th:
0209 if pos not in columnsWithSuspiciouslyNoisyPixels:
0210 rocNum, xCoordInROC = self.__convertCoordinatesFromHistToROCSpace(histName, pos, rocRow)
0211 noiseFile.write("%s\t(x, row)->[rocNum, xRoc]\t(%d, %d)->[%d, %d];\t{NORMVAL, TH}\t{%f, %f}\n" % (histName, pos, rocRow+1, rocNum, xCoordInROC, normVal, th))
0212
0213 return 1, noisyROC
0214
0215
0216
0217 return 0, noisyROC
0218
0219 def __determineEndcapNoise(self, noiseFile, columnsWithSuspiciouslyNoisyPixels, histName, meanOfPixels, linVal, val, pos, rocCol, rocRow):
0220 noisyROC = False;
0221 if meanOfPixels < self.rocOccupancyTh:
0222
0223 noisyROC = True
0224
0225 else:
0226 th = self.endcapNoisyColumnTh * linVal
0227 if val > th:
0228 if pos not in columnsWithSuspiciouslyNoisyPixels:
0229 rocNum, xCoordInROC = self.__convertCoordinatesFromHistToROCSpace(histName, pos, rocRow)
0230 noiseFile.write("%s\t(x, row)->[rocNum, xRoc]\t(%d, %d)->[%d, %d];\t{VAL, TH}\t{%f, %f}\n" % (histName, pos, rocRow+1, rocNum, xCoordInROC, val, th))
0231
0232 return 1, noisyROC
0233
0234
0235
0236 return 0, noisyROC
0237
0238 def __convertCoordinatesFromHistToROCSpace(self, histName, histXpos, histRow):
0239 tempXROC = (histXpos / self.rocMaxCol)
0240 tempYROC = histRow
0241
0242 tempXCoordInROC = histXpos % self.rocMaxCol
0243
0244 realXROC, realYROC = tempXROC, tempYROC
0245 xCoordInROC = tempXCoordInROC
0246
0247 rocNum = 0
0248
0249 if histName.find("BPix_Bp") != -1:
0250 realYROC = 1 - tempYROC
0251 if realYROC == 1:
0252 rocNum = 15 - realXROC
0253 xCoordInROC = self.rocMaxCol - 1 - xCoordInROC
0254 else:
0255 rocNum = realXROC
0256 else:
0257 realXROC = 7 - tempXROC
0258 if realYROC == 1:
0259 rocNum = 15 - realXROC
0260 else:
0261 rocNum = realXROC
0262 xCoordInROC = self.rocMaxCol - 1 - xCoordInROC
0263
0264 return rocNum, xCoordInROC
0265
0266
0267 def __determineBarrelDColInefficiencyAndNoise(self, medFiltRes, histName, pixelArr, pixelArrWithoutDrops, startPixel, rocCol, rocRow, outputFile, columnsWithSuspiciouslyNoisyPixels, noiseFile):
0268 meanOfPixels = sum(medFiltRes) / len(medFiltRes)
0269 maxMed = max(medFiltRes)
0270 minMed = min(medFiltRes)
0271
0272 normMeanOfPixels = sum(pixelArrWithoutDrops) / len(pixelArrWithoutDrops)
0273
0274
0275 doubleDeadCols = 0
0276 noisyColsNum = 0
0277 noisyROC = 0
0278
0279
0280 for i in range(1, len(pixelArr) - 2):
0281
0282 bin1valDiff = minMed - pixelArr[i + 0]
0283 bin2valDiff = minMed - pixelArr[i + 1]
0284
0285 bin0valDiff = minMed - pixelArr[i - 1]
0286 bin3valDiff = minMed - pixelArr[i + 2]
0287
0288
0289 currentDoubleBinThreshold = math.sqrt(meanOfPixels) * self.barrelInefficientDColTh
0290
0291 if bin1valDiff > currentDoubleBinThreshold and bin2valDiff > currentDoubleBinThreshold and not bin3valDiff > currentDoubleBinThreshold and not bin0valDiff > currentDoubleBinThreshold:
0292
0293 doubleColInRoc = ((i + startPixel) % (self.rocMaxCol)) // 2 + 1
0294 doubleDeadCols = doubleDeadCols + 1
0295
0296
0297 rocNum, xCoordInROC = self.__convertCoordinatesFromHistToROCSpace(histName, startPixel + i, rocRow)
0298 outputFile.write("%s\t(x, row)->[rocNum, doubleXPixelColInROC]\t(%d, %d)->[%d, %d];\t{MIN - VAL, TH}\t{%f, %f}\n" % (histName, startPixel + i, rocRow + 1, rocNum, xCoordInROC / 2, bin1valDiff, currentDoubleBinThreshold))
0299
0300
0301 if noisyROC == True:
0302 continue
0303
0304
0305
0306
0307
0308
0309
0310
0311 res = self.__determineBarrelNoise2(noiseFile, columnsWithSuspiciouslyNoisyPixels, histName, meanOfPixels, normMeanOfPixels, pixelArrWithoutDrops[i], startPixel + i, rocCol, rocRow)
0312 noisyColsNum, noisyROC = noisyColsNum + res[0], res[1]
0313 if i == len(pixelArr) - 3:
0314 res = self.__determineBarrelNoise2(noiseFile, columnsWithSuspiciouslyNoisyPixels, histName, meanOfPixels, normMeanOfPixels, pixelArrWithoutDrops[i + 1], startPixel + i + 1, rocCol, rocRow)
0315 noisyColsNum, noisyROC = noisyColsNum + res[0], res[1]
0316
0317
0318
0319 return doubleDeadCols, noisyColsNum
0320
0321 def __determineEndcapDColInefficiencyAndNoise(self, medFiltRes, histName, pixelArr, startPixel, rocCol, rocRow, outputFile, columnsWithSuspiciouslyNoisyPixels, noiseFile):
0322 doubleDeadCols = 0
0323 noisyColsNum = 0
0324 noisyROC = 0
0325
0326 useLin = True
0327
0328 a, b, D = self.__lmsLin(medFiltRes, startPixel, len(medFiltRes) + startPixel)
0329
0330 meanOfPixels = sum(medFiltRes) / len(medFiltRes)
0331
0332
0333 for i in range(1, len(pixelArr) - 2):
0334
0335 if useLin == True:
0336 linVal1 = a * (i + startPixel + 0) + b
0337 linVal2 = a * (i + startPixel + 1) + b
0338
0339 linVal0 = a * (i + startPixel - 1) + b
0340 linVal3 = a * (i + startPixel + 2) + b
0341 else:
0342 linVal1 = b * math.exp(a * (i + startPixel + 0))
0343 linVal2 = b * math.exp(a * (i + startPixel + 1))
0344
0345 linVal0 = b * math.exp(a * (i + startPixel - 1))
0346 linVal3 = b * math.exp(a * (i + startPixel + 2))
0347
0348 bin1valDiff = linVal1 - pixelArr[i + 0]
0349 bin2valDiff = linVal2 - pixelArr[i + 1]
0350
0351 bin0valDiff = linVal0 - pixelArr[i - 1]
0352 bin3valDiff = linVal3 - pixelArr[i + 2]
0353
0354 try:
0355 currentDoubleBinThreshold = math.sqrt((linVal1 + linVal2) * 0.5) * self.endcapInefficientDColTh
0356 except:
0357
0358 continue
0359
0360 if bin1valDiff > currentDoubleBinThreshold and bin2valDiff > currentDoubleBinThreshold and not bin3valDiff > currentDoubleBinThreshold and not bin0valDiff > currentDoubleBinThreshold:
0361
0362 doubleColInRoc = ((i + startPixel) % (self.rocMaxCol)) // 2 + 1
0363 doubleDeadCols = doubleDeadCols + 1
0364
0365
0366 rocNum, xCoordInROC = self.__convertCoordinatesFromHistToROCSpace(histName, startPixel + i, rocRow)
0367 outputFile.write("%s\t(x, row)->[rocNum, doubleXPixelColInROC]\t(%d, %d)->[%d, %d];\t{LIN(x) - VAL, TH}\t{%f, %f}\n" % (histName, startPixel + i, rocRow + 1, rocNum, xCoordInROC / 2, bin1valDiff, currentDoubleBinThreshold))
0368
0369
0370
0371 if noisyROC == True:
0372 continue
0373
0374 res = self.__determineEndcapNoise(noiseFile, columnsWithSuspiciouslyNoisyPixels, histName, meanOfPixels, linVal1, pixelArr[i], i + startPixel, rocCol, rocRow)
0375 noisyColsNum, noisyROC = noisyColsNum + res[0], res[1]
0376 if i == len(pixelArr) - 3:
0377 res = self.__determineEndcapNoise(noiseFile, columnsWithSuspiciouslyNoisyPixels, histName, meanOfPixels, linVal2, pixelArr[i + 1], i + 1 + startPixel, rocCol, rocRow)
0378 noisyColsNum, noisyROC = noisyColsNum + res[0], res[1]
0379
0380 return doubleDeadCols, noisyColsNum
0381
0382 def ReadHistograms(self):
0383 doubleDeadCols, noisyColsNum = 0, 0
0384
0385 with open(self.noiseOutputFileName, "w") as noiseFile:
0386 with open(self.outputFileName, "w") as outputFile:
0387 for layer in self.dicOfModuleHistograms:
0388
0389 doubleDeadColsInLayer, noisyColsNumInLayer = 0, 0
0390
0391 outputFile.write("-> " + layer + "\n\n")
0392 noiseFile.write("-> " + layer + "\n\n")
0393
0394 for hist in self.dicOfModuleHistograms[layer]:
0395 for row in range(2):
0396 for rocNum in range(self.rocsInRow):
0397 startPixel = rocNum * self.rocMaxCol + 1
0398 endPixel = (rocNum + 1) * self.rocMaxCol + 1
0399
0400 rocCol = rocNum + 1
0401
0402 pixelArr, medFiltRes, columnsWithSuspiciouslyNoisyPixels, sciPyMedFiltRes = self.__getROCData(hist, startPixel, endPixel, row, 3, 5)
0403
0404 if pixelArr == None:
0405 continue
0406
0407
0408
0409
0410
0411
0412 if "F" not in layer:
0413
0414 pixelArrWithoutDrops = self.__getPixelArrWithRemovedDrops(pixelArr, sciPyMedFiltRes)
0415 pixelArrWithoutDropsNormalized = self.__normalizeArray(pixelArrWithoutDrops)
0416
0417
0418
0419
0420 result = self.__determineBarrelDColInefficiencyAndNoise(medFiltRes, hist.GetName(), pixelArr, pixelArrWithoutDropsNormalized, startPixel, rocCol, row, outputFile, columnsWithSuspiciouslyNoisyPixels, noiseFile)
0421 else:
0422 result = self.__determineEndcapDColInefficiencyAndNoise(medFiltRes, hist.GetName(), pixelArr, startPixel, rocCol, row, outputFile, columnsWithSuspiciouslyNoisyPixels, noiseFile)
0423
0424 doubleDeadCols, doubleDeadColsInLayer = doubleDeadCols + result[0], doubleDeadColsInLayer + result[0]
0425 noisyColsNum, noisyColsNumInLayer = noisyColsNum + result[1], noisyColsNumInLayer + result[1]
0426
0427 outputFile.write("\n\tTOTAL IN LAYER/DISK: %d\n\n" % (doubleDeadColsInLayer))
0428 noiseFile.write("\n\tTOTAL IN LAYER/DISK: %d\n\n" % (noisyColsNumInLayer))
0429
0430 print("Number of inefficient double columns: %d"%(doubleDeadCols))
0431 print("Number of noisy cols: %d"%(noisyColsNum))
0432
0433
0434
0435 for i in range(1, len(sys.argv), 1):
0436 if i == 1:
0437 inputFileName = sys.argv[i]
0438
0439 runNum = ((inputFileName.split("/")[-1].split("."))[0].split("_R000"))[1]
0440 print("Run number: %s"%(runNum))
0441 baseRootDir = ["DQMData/Run " + runNum + "/PixelPhase1/Run summary/Phase1_MechanicalView"]
0442 print(baseRootDir[0])
0443 outputFileName = "inefficientDPixelColumns.txt"
0444 noiseOutputFileName = "noisyPixelColumns.txt"
0445
0446 readerObj = InefficientDeadROCs(inputFileName, outputFileName, noiseOutputFileName, baseRootDir)
0447 readerObj.ReadHistograms()