File indexing completed on 2024-04-06 12:08:42
0001
0002
0003 from __future__ import print_function
0004 import sys
0005 import math
0006 from ROOT import gROOT, TClass, TFile
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()