SiStripMonitorClient/scripts/TH2PolyOfflineMaps.py

0001 #!/usr/bin/env python3
0002
0003 from __future__ import print_function
0004 import sys
0005 import os
0006 from ROOT import *
0007 from copy import deepcopy
0008 from array import array
0009
0010 gROOT.SetBatch()        # don't pop up canvases
0011
0012 #Find Data files
0013
0014 def getFileInPath(rfile):
0015    import os
0016    for dir in os.environ['CMSSW_SEARCH_PATH'].split(":"):
0017      if os.path.exists(os.path.join(dir,rfile)): return os.path.join(dir,rfile)
0018    return None
0019
0020
0021 # Default values
0022 inputFileName = "DQM_V0013_R000292154__StreamExpressCosmics__Commissioning2017-Express-v1__DQMIO.root"
0023 limitsFileName = "limits.dat"
0024 outputDirectoryName = "OUT/"
0025 minMaxFileName = "minmax.out"
0026 #detIDsFileName = "DATA/detids.dat"
0027 detIDsFileName = getFileInPath('DQM/SiStripMonitorClient/data/detids.dat')
0028 #default one
0029 baseRootDirs = ["DQMData/Run 292154/PixelPhase1/Run summary/Phase1_MechanicalView"
0030                 ,"DQMData/Run 292154/PixelPhase1/Run summary/Tracks"
0031                 ]
0032
0033
0034 maxPxBarrel = 4
0035 maxPxForward = 3
0036 barrelLadderShift = [0, 14, 44, 90]
0037
0038 forwardDiskXShift = [25, 75, 125]
0039 forwardDiskYShift = 45; # to make +DISK on top in the 'strip-like' layout
0040
0041 plotWidth, plotHeight = 3000, 2000
0042 extremeBinsNum = 20
0043
0044 limits = ["num_digis 0.01 90 1 0",
0045           "num_clusters 0.01 25 1 0",
0046           "Trechitsize_y 0.01 10 0 0",
0047           "Trechitsize_x 0.01 10 0 0",
0048           "Tresidual_y 0.0000001 0.004 0 1",
0049           "Tresidual_x 0.0000001 0.004 0 1",
0050           "Tcharge 2000 80000 0 0",
0051           "Thitefficiency 0.95 1 0 0",
0052           #"Tmissing 0.01 500 0 0",
0053           "Tnum_clusters_ontrack 0.01 15 1 0",
0054           "Tsize 0.01 15 0 0",
0055           #"Tvalid 0.01 90 0 0",
0056           "adc 0.01 256 0 0",
0057           "charge 2000 80000 0 0",
0058           "size 0.01 15 0 0",]
0059
0060 class TH2PolyOfflineMaps:
0061
0062   ###
0063   # LOTS OF CODE BORROWED FROM: PYTHONBINREADER, PIXELTRACKERMAP
0064   ###
0065
0066   ############################################################################
0067
0068   def __TraverseDirTree(self, dir):
0069
0070     try:
0071       currPath = (dir.GetPath().split(":/"))[1]
0072     except:
0073       print("Exception raised: Path not found in the input file")
0074       return
0075
0076     for obj in dir.GetListOfKeys():
0077       if not obj.IsFolder():
0078         if obj.ReadObjectAny(TClass.GetClass("TH2")):
0079           th2 = deepcopy(obj.ReadObj())
0080           name = th2.GetName()
0081           if 6 < th2.GetNbinsX() < 10 and name.find("per") != -1 and name.find("Lumisection") == -1: #take only module lvl plots
0082             print(''.join([dir.GetPath(), '/', name]))
0083
0084             # fix when there are plots starting with the same strings in different directories
0085             prefix = ""
0086             for i in self.dirs:
0087               if currPath.startswith(i):
0088                 prefix = self.dirsAliases[i]
0089                 break
0090             # print(currPath, prefix)
0091             th2.SetName(prefix + th2.GetName())
0092             self.listOfNumHistograms.append(th2)
0093       else:
0094         self.__TraverseDirTree(obj.ReadObj())
0095
0096   def __GetPartStr(self, isXlowerThanZero, isYlowerThanZero):
0097     if isXlowerThanZero and isYlowerThanZero:
0098       return "mO"
0099     if isXlowerThanZero and isYlowerThanZero == False:
0100       return "mI"
0101     if isXlowerThanZero == False and isYlowerThanZero:
0102       return "pO"
0103     if isXlowerThanZero == False and isYlowerThanZero == False:
0104       return "pI"
0105
0106   def __GetBarrelSector(self, layer, signedLadder, signedModule): #adapted from PixelBarrelName
0107     theLadder = abs(signedLadder)
0108     theModule = abs(signedModule)
0109
0110     sector = 0
0111
0112     if layer == 1:
0113
0114       if theLadder == 1:
0115         if theModule >= 2:
0116           return 1
0117         else:
0118           return 2
0119       if theLadder == 2:
0120         if theModule >= 3:
0121           return 2
0122         else:
0123           return 3
0124       if theLadder == 3:
0125         if theModule >= 4:
0126           return 3
0127         else:
0128           return 4
0129       if theLadder == 4:
0130         if theModule >= 2:
0131           return 5
0132         else:
0133           return 6
0134       if theLadder == 5:
0135         if theModule >= 3:
0136           return 6
0137         else:
0138           return 7
0139       if theLadder == 6:
0140         if theModule >= 4:
0141           return 7
0142         else:
0143           return 8
0144     # here is used simplified form of assignment, see source file for reference
0145     elif layer == 2:
0146       i = theLadder // 5
0147       sector = i * 3
0148       shortLadder = theLadder - 5 * i
0149       for i in range(0, shortLadder, 2):
0150         sector = sector + 1
0151       return sector
0152     elif layer == 3:
0153       sector = 1
0154       for i in range(2, theLadder, 3):
0155         if (i + 1) % 3 == 0:
0156           sector = sector + 1
0157       return sector
0158     elif layer == 4:
0159       sector = (theLadder + 3) // 4
0160       return sector
0161
0162   def __BuildOnlineBarrelName(self, signedModule, signedLadder, layer): #in Phase1 it is assumed that there are only full modules
0163     thePart = self.__GetPartStr(signedModule < 0, signedLadder < 0)
0164     theSector = str(self.__GetBarrelSector(layer, signedLadder, signedModule))
0165     return "BPix_B" + thePart + "_SEC" + theSector + "_LYR" + str(layer) + "_LDR" + str(abs(signedLadder)) + "F_MOD" + str(abs(signedModule))
0166
0167   def __BuildOnlineDiskName(self, signedDisk, signedBlade, panel, ring):
0168     thePart = self.__GetPartStr(signedDisk < 0, signedBlade < 0)
0169     return "FPix_B" + thePart + "_D" + str(abs(signedDisk)) + "_BLD" + str(abs(signedBlade)) + "_PNL" + str(panel) + "_RNG" + str(ring)
0170
0171   def __GroupHistograms(self):
0172     currentGroupName = ""
0173     groupOfHists = []
0174     self.groupedHistograms = []
0175
0176     ##### GROUP ALL LAYERS/RINGS HAVING SIMILAR INFORMATION
0177     for obj in self.listOfNumHistograms:
0178       objName = obj.GetName()
0179       objNameSplit = objName.split("_")
0180       objNameCollected = ''.join(objNameSplit[0:-1])
0181       if objNameCollected != currentGroupName:
0182         if len(groupOfHists):
0183           self.groupedHistograms.append(groupOfHists)
0184           groupOfHists = []
0185
0186         currentGroupName = objNameCollected
0187       groupOfHists.append(obj)
0188     self.groupedHistograms.append(groupOfHists) #the last group
0189
0190   def __AddNamedBins(self, geoFile, tX, tY, sX, sY, applyModuleRotation = False):
0191
0192     for line in geoFile:
0193       lineSpl = line.strip().split("\"")
0194
0195       detId = lineSpl[0].split(" ")[0]
0196
0197       vertices = lineSpl[1]
0198       xy = vertices.split(" ")
0199       x, y = array('d'), array('d')
0200       verNum = 1
0201       for coord in xy:
0202         coordSpl = coord.split(",")
0203         if applyModuleRotation:
0204           x.append(-(float(coordSpl[0]) * sX + tX))
0205           y.append((float(coordSpl[1]) * sY + tY))
0206         else:
0207           x.append(float(coordSpl[0]) * sX + tX)
0208           y.append(float(coordSpl[1]) * sY + tY)
0209         verNum = verNum + 1
0210       #close polygon
0211       x.append(x[0])
0212       y.append(y[0])
0213
0214       # print(detId, vertices)
0215       # print(x)
0216       # print(y)
0217       if applyModuleRotation:
0218         bin = TGraph(verNum, y, x)
0219       else:
0220         bin = TGraph(verNum, x, y)
0221       # bin = TGraph(verNum, y, x) # rotation by 90 deg (so that it had the same layout as for the strips)
0222       bin.SetName(detId)
0223
0224       self.__BaseTrackerMap.AddBin(bin)
0225
0226   def __CreateTrackerBaseMap(self):
0227
0228     self.__BaseTrackerMap = TH2Poly("Summary", "", -10, 160, -70, 70)
0229     # self.__BaseTrackerMap = TH2Poly("Summary", "Tracker Map", 0, 0, 0, 0)
0230     self.__BaseTrackerMap.SetFloat(1)
0231     self.__BaseTrackerMap.GetXaxis().SetTitle("")
0232     self.__BaseTrackerMap.GetYaxis().SetTitle("")
0233     self.__BaseTrackerMap.SetOption("COLZ L")
0234     self.__BaseTrackerMap.SetStats(0)
0235
0236     # BARREL FIRST
0237     for i in range(maxPxBarrel):
0238       with open(self.geometryFilenames[i], "r") as geoFile:
0239         currBarrelTranslateX = 0
0240         currBarrelTranslateY = barrelLadderShift[i]
0241
0242         self.__AddNamedBins(geoFile, currBarrelTranslateX, currBarrelTranslateY, 1, 1, True)
0243
0244       # break # debug only 1st layer
0245
0246     # MINUS FORWARD
0247     for i in range(-maxPxForward, 0):
0248       with open(self.geometryFilenames[maxPxBarrel + maxPxForward + i], "r") as geoFile:
0249         currForwardTranslateX = forwardDiskXShift[-i - 1]
0250         currForwardTranslateY = -forwardDiskYShift
0251
0252         self.__AddNamedBins(geoFile, currForwardTranslateX, currForwardTranslateY, 1, 1)
0253
0254     # PLUS FORWARD
0255     for i in range(maxPxForward):
0256       with open(self.geometryFilenames[maxPxBarrel + maxPxForward + i], "r") as geoFile:
0257         currForwardTranslateX = forwardDiskXShift[i]
0258         currForwardTranslateY = forwardDiskYShift
0259
0260         self.__AddNamedBins(geoFile, currForwardTranslateX, currForwardTranslateY, 1, 1)
0261
0262     # self.__BaseTrackerMap.Fill("305139728", 2)
0263
0264     print("Base Tracker Map: constructed")
0265
0266   ############################################################################
0267   def __init__(self, inputDQMName, outputDirName, minMaxFileName, limits,  modDicName, runNumber, dirs, dirsAliases):
0268 #  def __init__(self, inputDQMName, outputDirName, minMaxFileName, limitsFileName, modDicName, runNumber, dirs, dirsAliases):
0269     self.inputFileName = inputDQMName
0270     self.outputDirName = outputDirName
0271     self.minMaxFileName = minMaxFileName
0272 #    self.limitsFileName = limitsFileName
0273     self.detIDsFileName = modDicName
0274     self.limits = limits
0275
0276     self.runNumber = runNumber
0277     self.dirs = dirs
0278     self.dirsAliases = dirsAliases
0279
0280     self.inputFile = TFile(self.inputFileName)
0281     self.listOfNumHistograms = []
0282     self.availableNames = []
0283
0284     self.maxLadderToLayer = {6:1, 14:2, 22:3, 32:4}
0285     self.maxBladeToRing = {11:1, 17:2}
0286
0287     self.geometryFilenames = []
0288     for i in range(maxPxBarrel):
0289        self.geometryFilenames.append(getFileInPath("DQM/SiStripMonitorClient/data/Geometry/vertices_barrel_" + str(i + 1)))
0290 #      self.geometryFilenames.append("DATA/Geometry/vertices_barrel_" + str(i + 1))
0291     for i in range(-maxPxForward, maxPxForward + 1):
0292       if i == 0:
0293         continue #there is no 0 disk
0294       self.geometryFilenames.append(getFileInPath("DQM/SiStripMonitorClient/data/Geometry/vertices_forward_" + str(i)))
0295 #      self.geometryFilenames.append("DATA/Geometry/vertices_forward_" + str(i))
0296
0297     self.internalData = {}
0298
0299     if self.inputFile.IsOpen():
0300       print("%s opened successfully!" % (self.inputFileName))
0301       #Get all neeeded histograms
0302       for dir in self.dirs:
0303         self.__TraverseDirTree(self.inputFile.Get(dir))
0304       # print("Histograms to read %d" % (len(self.listOfNumHistograms)))
0305
0306       self.detDict = {}
0307
0308       with open(self.detIDsFileName, "r") as detIDs:  # create dictionary online -> rawid
0309         for entry in detIDs:
0310           items = entry.replace("\n", " ").split(" ")
0311           self.detDict.update({items[1] : int(items[0])})
0312           # init internal data structure
0313           self.internalData.update({int(items[0]) : {}})
0314
0315       self.rawToOnlineDict = dict((v,k) for k,v in self.detDict.items())
0316
0317       self.__GroupHistograms()
0318
0319       self.__CreateTrackerBaseMap()
0320
0321     else:
0322       print("Unable to open file %s" % (self.inputFileName))
0323
0324     ### CREATE LIMITS DICTIONARY
0325
0326     self.limitsDic = {}
0327     for y in limits:
0328
0329       lineSpl = y.strip().split(" ")
0330
0331       if len(lineSpl) < 5:
0332         continue
0333
0334       currName = lineSpl[0]
0335       zMin = float(lineSpl[1])
0336       zMax = float(lineSpl[2])
0337       isLog = False if lineSpl[3] == "0" else True
0338       isAbs = False if lineSpl[4] == "0" else True
0339
0340       self.limitsDic.update({currName : {"zMin" : zMin, "zMax" : zMax, "isLog" : isLog, "isAbs" : isAbs}})
0341  #     print limitsDic
0342
0343   def ReadHistograms(self):
0344     if self.inputFile.IsOpen():
0345       for group in self.groupedHistograms:
0346         # name = ''.join(group[0].GetName().split("_")[0:-1])
0347         if len(group) == 0:
0348           return
0349         print(group[0].GetName())
0350         name = ''.join(group[0].GetName().split("_per_")[0])
0351         self.availableNames.append(name)
0352         # print(name)
0353         for obj in group:
0354           nbinsX = obj.GetNbinsX()
0355           nbinsY = obj.GetNbinsY()
0356
0357           if nbinsX == 9: # BARREL
0358             maxX = nbinsX // 2
0359             maxY = nbinsY // 2
0360
0361             for x in range(-maxX, maxX + 1):
0362               if x == 0:
0363                 continue
0364               for y in range(-maxY, maxY + 1, 1):
0365                 if y == 0:
0366                   continue
0367                 onlineName = self.__BuildOnlineBarrelName(x, y, self.maxLadderToLayer[maxY])
0368                 self.internalData[self.detDict[onlineName]].update({name : obj.GetBinContent(x + maxX + 1, y + maxY + 1)})
0369
0370           elif nbinsX == 7: # FORWARD
0371             maxX = nbinsX // 2
0372             maxY = nbinsY // 4
0373
0374             for x in range(-maxX, maxX + 1):
0375               if x == 0:
0376                 continue
0377               for y in range(-maxY, maxY + 1):
0378                 if int(y) == 0:
0379                   continue
0380                 for panel in range(1, 3):
0381                   onlineName = self.__BuildOnlineDiskName(x, y, panel, self.maxBladeToRing[maxY])
0382                   self.internalData[self.detDict[onlineName]].update({name : obj.GetBinContent(x + maxX + 1, (y + maxY) * 2 + (3-panel))})
0383           else:
0384             print("Unrecognized plot")
0385       else:
0386         print("Histograms saved to internal data structure")
0387
0388   def DumpData(self):
0389     for key in self.internalData:
0390       print("#"*20)
0391       print(key)
0392       module = self.internalData[key]
0393       for d in module:
0394         print((d, module[d]))
0395
0396     print(len(self.internalData))
0397
0398     for i in self.availableNames:
0399       print(i)
0400     print(len(self.availableNames))
0401
0402   def PrintTrackerMaps(self):
0403     monitoredValues = []
0404     gStyle.SetPalette(1)
0405     for key in self.internalData:
0406       monitoredValues = self.internalData[key].keys()
0407       # print(monitoredValues)
0408       break
0409
0410     if os.path.exists(self.outputDirName) == False: # check whether directory exists
0411       os.system("mkdir " + self.outputDirName)
0412
0413     with open(self.outputDirName + self.minMaxFileName, "w") as minMaxFile:
0414
0415       for mv in monitoredValues:
0416         currentHist = deepcopy(self.__BaseTrackerMap)
0417         # currentHist.SetTitle("Run " + self.runNumber + ": Tracker Map for " + mv) // to make it compatible between ROOT v.
0418         histoTitle = "Run " + self.runNumber + ": Tracker Map for " + mv
0419
0420         applyLogScale = False
0421         applyAbsValue = False
0422         if mv in self.limitsDic:
0423           limitsElem = self.limitsDic[mv]
0424
0425           print(mv + " found in limits dictionary - applying custom limits...")
0426
0427           currentHist.SetMinimum(limitsElem["zMin"])
0428           currentHist.SetMaximum(limitsElem["zMax"])
0429           applyLogScale = limitsElem["isLog"]
0430           applyAbsValue = limitsElem["isAbs"]
0431
0432         listOfVals = []
0433         onlineName = ""
0434         for detId in self.internalData:
0435           val = (self.internalData[detId])[mv]
0436           onlineName = self.rawToOnlineDict[detId]
0437           listOfVals.append([val, detId, onlineName])
0438           if applyAbsValue:
0439              currentHist.Fill(str(detId), abs(val))
0440           else:
0441              currentHist.Fill(str(detId), val)
0442
0443         listOfVals = sorted(listOfVals, key = lambda item: item[0])
0444
0445         minMaxFile.write("\n" + mv + "\n\n")
0446
0447         minMaxFile.write("MIN:\n")
0448         for i in range(extremeBinsNum):
0449           minMaxFile.write("\t" + str(listOfVals[i][1]) + " " + str(listOfVals[i][2]) + " " + str(listOfVals[i][0]) + "\n")
0450
0451         minMaxFile.write("MAX:\n")
0452         for i in range(extremeBinsNum):
0453           minMaxFile.write("\t" + str(listOfVals[-i - 1][1]) + " " + str(listOfVals[-i - 1][2]) + " " + str(listOfVals[-i - 1][0]) + "\n")
0454
0455         c1 = TCanvas(mv, mv, plotWidth , plotHeight)
0456
0457         if applyLogScale:
0458           c1.SetLogz()
0459
0460
0461
0462         currentHist.Draw("AC COLZ L")
0463
0464         gPad.Update()
0465         palette = currentHist.FindObject("palette");
0466         palette.SetX1NDC(0.89);
0467         palette.SetX2NDC(0.91);
0468         palette.SetLabelSize(0.05);
0469         gPad.Update()
0470
0471
0472
0473         ### IMPORTANT - REALTIVE POSITIONING IS MESSY IN CURRENT VERION OF PYROOT
0474         ### IT CAN CHANGE FROM VERSION TO VERSION, SO YOU HAVE TO ADJUST IT FOR YOUR NEEDS
0475         ### !!!!!!!!!!!!!
0476
0477         # draw axes (z, phi -> BARREL; x, y -> FORWARD)
0478         ###################################################
0479
0480         ### z arrow
0481         arrow = TArrow(0.05, 27.0, 0.05, -30.0, 0.02, "|>")
0482         arrow.SetLineWidth(4)
0483         arrow.Draw()
0484         ### phi arrow
0485         phiArrow = TArrow(0.0, 27.0, 30.0, 27.0, 0.02, "|>")
0486         phiArrow.SetLineWidth(4)
0487         phiArrow.Draw()
0488         ### x arror
0489         xArrow = TArrow(25.0, 44.5, 50.0, 44.5, 0.02, "|>")
0490         xArrow.SetLineWidth(4)
0491         xArrow.Draw()
0492         ### y arror
0493         yArrow = TArrow(25.0, 44.5, 25.0, 69.5, 0.02, "|>")
0494         yArrow.SetLineWidth(4)
0495         yArrow.Draw()
0496
0497         ###################################################
0498
0499         # add some captions
0500         txt = TLatex()
0501         txt.SetNDC()
0502         txt.SetTextFont(1)
0503         txt.SetTextColor(1)
0504         txt.SetTextAlign(22)
0505         txt.SetTextAngle(0)
0506
0507         # draw new-style title
0508         txt.SetTextSize(0.05)
0509         txt.DrawLatex(0.5, 0.95, histoTitle)
0510
0511         txt.SetTextSize(0.03)
0512
0513         txt.DrawLatex(0.5, 0.125, "-DISK")
0514         txt.DrawLatex(0.5, 0.075, "NUMBER ->")
0515         txt.DrawLatex(0.5, 0.875, "+DISK")
0516
0517         txt.DrawLatex(0.17, 0.35, "+z")
0518         txt.DrawLatexNDC(0.36, 0.685, "+phi") # WAY TO FORCE IT TO DRAW LATEX CORRECTLY NOT FOUND ('#' DOESN'T WORK)
0519         txt.DrawLatex(0.38, 0.73, "+x")
0520         txt.DrawLatex(0.26, 0.875, "+y")
0521
0522         txt.SetTextAngle(90)
0523         txt.DrawLatex(0.17, 0.5, "BARREL")
0524
0525         #save to the png
0526         c1.Print(self.outputDirName + mv + ".png")
0527
0528   def __del__(self):
0529     if self.inputFile :
0530       if self.inputFile.IsOpen():
0531         self.inputFile.Close()
0532
0533 #--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
0534 for i in range(1, len(sys.argv), 1):
0535   if i == 1:
0536     inputFileName = sys.argv[i]
0537   elif i == 2:
0538     plotWidth = int(sys.argv[i])
0539   elif i == 3:
0540     plotHeight = int(sys.argv[i])
0541 #  elif i == 4:
0542 #    limitsFileName = sys.argv[i]
0543 #  elif i == 5:
0544   elif i == 4:
0545     detIDsFileName = sys.argv[i]
0546
0547 deductedRunNumber = inputFileName.split("_R000")[1][0:6]
0548 print(deductedRunNumber)
0549
0550 baseRootDirs = ["DQMData/Run " + deductedRunNumber + "/PixelPhase1/Run summary/Phase1_MechanicalView"    #maybe read it from the input file???
0551                 ,"DQMData/Run " + deductedRunNumber + "/PixelPhase1/Run summary/Tracks"
0552                 ]
0553
0554 baseRootDirsAliases = {baseRootDirs[0]:""
0555                     , baseRootDirs[1]:"T"
0556                     }
0557
0558 readerObj = TH2PolyOfflineMaps(inputFileName, outputDirectoryName, minMaxFileName, limits, detIDsFileName, deductedRunNumber, baseRootDirs, baseRootDirsAliases)
0559 #readerObj = TH2PolyOfflineMaps(inputFileName, outputDirectoryName, minMaxFileName, limitsFileName, detIDsFileName, deductedRunNumber, baseRootDirs, baseRootDirsAliases)
0560 readerObj.ReadHistograms()
0561 # readerObj.DumpData()
0562 readerObj.PrintTrackerMaps()