Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:46:00

0001 from __future__ import absolute_import
0002 from __future__ import print_function
0003 
0004 from builtins import range
0005 import math
0006 
0007 import ROOT
0008 
0009 from .TkAlStyle import TkAlStyle
0010 
0011 dirNameList=["z","r","phi"]# in general directions are labeled z=0 r =1 phi =2 throughout this, I should probably think of something more elegant
0012 detNameList = ("BPIX", "FPIX", "TIB", "TID", "TOB", "TEC")
0013 def subdetCondition(subdet,module,overlap):
0014         if (subdet+1 == 1 and (module == 1 or overlap ==1))or (subdet+1== 2 and (module == 0 or overlap == 0))or ((subdet+1== 3 or subdet+1== 5) and (overlap != 2 or module == 1))or ((subdet+1== 4 or subdet+1== 6) and (overlap != 2 or module == 0)):
0015                 return True
0016         else:
0017                 return False
0018 def hist(tree_file_name, hist_name,subdet_ids,module_directions,overlap_directions,profile_directions):
0019     f = ROOT.TFile(tree_file_name)
0020     t = f.Get("analysis/Overlaps")
0021     h = []
0022     for subdet in range(6):#Creates a 4-D list of empty histograms for permutations of subdetector, overlap direction, module direction and profile direction.
0023                 h.append([])
0024                 for module in range(3):
0025                         h[subdet].append([])
0026                         for overlap in range(3):
0027                                 h[subdet][module].append([])
0028                                 for profile in range(4):
0029                                         if subdetConditon(subdet,module,overlap):
0030                                                 h[subdet][module][overlap].append(0)
0031                                                 continue
0032                                         name = hist_name + "{0}_{1}_{2}".format(dirNameList[module],dirNameList[overlap],detNameList[subdet])
0033                                         if not ((subdet_ids[subdet]) and  (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]): #puts a 0 in any unwanted plots
0034                                                 h[subdet][module][overlap].append(0)
0035                                                 continue
0036                                         if (profile>0):
0037                                                 if profile == 3: bins, xmin, xmax = 10, -math.pi, math.pi
0038                                                 if subdet+1 == 1 and profile == 1: bins, xmin, xmax = 10, -30, -30
0039                                                 if subdet+1 == 2 and profile == 1: bins, xmin, xmax = 40, -60, 60 
0040                                                 if subdet+1 == 3 and profile == 1: bins, xmin, xmax = 10, -70, 70 
0041                                                 if subdet+1 == 4 and profile == 1: bins, xmin, xmax = 40, -110, 110 
0042                                                 if subdet+1 == 5 and profile == 1: bins, xmin, xmax = 10, -110, 110
0043                                                 if subdet+1 == 6 and profile == 1: bins, xmin, xmax = 20, -280, 280 
0044                                                 if subdet+1 == 1 and profile == 2: bins, xmin, xmax = 10, 0, 20 
0045                                                 if subdet+1 == 2 and profile == 2: bins, xmin, xmax = 10, 0, 20
0046                                                 if subdet+1 == 3 and profile == 2: bins, xmin, xmax = 10, 20, 55 
0047                                                 if subdet+1 == 4 and profile == 2: bins, xmin, xmax = 10, 20, 55 
0048                                                 if subdet+1 == 5 and profile == 2: bins, xmin, xmax = 10, 55, 110 
0049                                                 if subdet+1 == 6 and profile == 2: bins, xmin, xmax = 10, 20, 110 
0050                                                 h[subdet][module][overlap].append(ROOT.TProfile(hist_name, hist_name, bins, xmin, xmax)) 
0051                                         elif subdet+1==4 or subdet+1==6:
0052                                             h[subdet][module][overlap].append(ROOT.TH1F(name, name, 100, -5000, 5000))
0053                                         else:
0054                                             h[subdet][module][overlap].append(ROOT.TH1F(name, name, 100, -300, 300))                                            
0055                                         h[subdet][module][overlap][profile].SetDirectory(0)
0056 
0057     nentries = t.GetEntries()
0058 
0059     for i, entry in enumerate(t, start=1):#loops through the tree, filling in relevant histograms for each event
0060         if i % 10000 == 0 or i == nentries:
0061             print(i, "/", nentries)
0062         
0063         subdet_id = t.subdetID
0064         if subdet_ids[subdet_id-1]==False:
0065                 continue
0066 
0067         modulePhi0 = math.atan2(t.moduleY[0], t.moduleX[0]) 
0068         modulePhi1 = math.atan2(t.moduleY[1], t.moduleX[1])
0069         phidiff = min(abs(modulePhi0-modulePhi1), abs(math.pi - abs(modulePhi0-modulePhi1)))
0070         moduleR0 = math.sqrt(t.moduleY[0]**2+t.moduleX[0]**2)
0071         moduleR1 = math.sqrt(t.moduleY[1]**2+t.moduleX[1]**2)
0072 
0073         #determines the direction the modules are in respect to each other for each event
0074         if subdet_id%2 == 1 and ((abs(t.moduleZ[0] - t.moduleZ[1]) > 1)): 
0075             module_direction = 0 #0 refers to Z, 1 to R and 2 to phi
0076         elif subdet_id%2 == 0 and (abs(moduleR0 - moduleR1)>1):
0077             module_direction = 1        
0078         elif ((moduleR0*phidiff)>1):
0079             module_direction = 2
0080         else: #prints if this method of selecting module_direction misses anything
0081             print(str(i)+" skipped")
0082             continue
0083         if module_directions[module_direction]==False:
0084                 continue
0085         avgList=[(t.moduleZ[0]+t.moduleZ[1])/2,(moduleR0+moduleR1)/2,(modulePhi0+modulePhi1)/2]
0086         if overlap_directions[2]:
0087              overlap_direction = 2
0088              if modulePhi0 > modulePhi1:
0089                  hitXA = t.hitX[1]
0090                  hitXB = t.hitX[0]
0091                  predXA = t.predX[1]
0092                  predXB = t.predX[0]
0093                  overlapSignA = t.localxdotglobalphi[1]
0094                  overlapSignB = t.localxdotglobalphi[0]
0095              else:
0096                  hitXA = t.hitX[0]
0097                  hitXB = t.hitX[1]
0098                  predXA = t.predX[0]
0099                  predXB = t.predX[1]
0100                  overlapSignA = t.localxdotglobalphi[0]
0101                  overlapSignB = t.localxdotglobalphi[1]
0102              residualA = hitXA - predXA
0103              residualB = hitXB - predXB
0104              if overlapSignA < 0:
0105                  residualA *= -1
0106              if overlapSignB < 0:
0107                  residualB *= -1
0108              A = 10000*(residualA - residualB)
0109              h[subdet_id-1][module_direction][overlap_direction][0].Fill(A)
0110              for profile in range(3):
0111                  if profile_directions[profile+1]:
0112                         h[subdet_id-1][module_direction][overlap_direction][profile+1].Fill(avgList[profile],A)
0113                  
0114                         
0115         if subdet_id==1 and module_direction != 1 and overlap_directions[0]:
0116              overlap_direction = 0
0117              if t.moduleZ[0] > t.moduleZ[1]:
0118                  hitXA = t.hitY[1]
0119                  hitXB = t.hitY[0]
0120                  predXA = t.predY[1]
0121                  predXB = t.predY[0]
0122                  overlapSignA = t.localydotglobalz[1]
0123                  overlapSignB = t.localydotglobalz[0]
0124              else:
0125                  hitXA = t.hitY[0]
0126                  hitXB = t.hitY[1]
0127                  predXA = t.predY[0]
0128                  predXB = t.predY[1]
0129                  overlapSignA = t.localydotglobalz[0]
0130                  overlapSignB = t.localydotglobalz[1]
0131              residualA = hitXA - predXA
0132              residualB = hitXB - predXB
0133              if overlapSignA < 0:
0134                 residualA *= -1
0135              if overlapSignB < 0:
0136                 residualB *= -1
0137              A = 10000*(residualA - residualB)
0138              h[subdet_id-1][module_direction][overlap_direction][0].Fill(A)
0139              for profile in range(3):
0140                  if profile_directions[profile+1]:
0141                         h[subdet_id-1][module_direction][overlap_direction][profile+1].Fill(avgList[profile],A)
0142 
0143         if subdet_id==2 and module_direction !=0 and overlap_directions[1]:
0144              overlap_direction = 1
0145              if moduleR0 > moduleR1:
0146                  hitXA = t.hitY[1]
0147                  hitXB = t.hitY[0]
0148                  predXA = t.predY[1]
0149                  predXB = t.predY[0]
0150                  overlapSignA = t.localydotglobalr[1]
0151                  overlapSignB = t.localydotglobalr[0]
0152              else:
0153                  hitXA = t.hitY[0]
0154                  hitXB = t.hitY[1]
0155                  predXA = t.predY[0]
0156                  predXB = t.predY[1]
0157                  overlapSignA = t.localydotglobalr[0]
0158                  overlapSignB = t.localydotglobalr[1]
0159              
0160              residualA = hitXA - predXA
0161              residualB = hitXB - predXB
0162              if overlapSignA < 0:
0163                  residualA *= -1
0164              if overlapSignB < 0:
0165                  residualB *= -1
0166              A = 10000*(residualA - residualB)
0167              h[subdet_id-1][module_direction][overlap_direction][0].Fill(A)
0168              for profile in range(3):
0169                  if profile_directions[profile+1]:
0170                         h[subdet_id-1][module_direction][overlap_direction][profile+1].Fill(avgList[profile],A)
0171     return h
0172 
0173 def plot(file_name,subdet_ids,module_directions,overlap_directions,profile_directions,*filesTitlesColorsStyles):
0174 
0175         legend=[]
0176         hstack=[]
0177         
0178         for subdet in range(6):#creates lists of empty THStacks and legends to be filled later
0179                 hstack.append([])
0180                 legend.append([])
0181                 for module in range(3):
0182                         hstack[subdet].append([])
0183                         legend[subdet].append([])
0184                         for overlap in range(3):
0185                                 hstack[subdet][module].append([])
0186                                 legend[subdet][module].append([])
0187                                 for profile in range(4):
0188                                         if subdetCondition(subdet,module,overlap):
0189                                                 hstack[subdet][module][overlap].append(0)
0190                                                 legend[subdet][module][overlap].append(0)
0191                                                 continue
0192                                         if not ((subdet_ids[subdet]) and  (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]):
0193                                                 legend[subdet][module][overlap].append(0)
0194                                                 hstack[subdet][module][overlap].append(0)
0195                                                 continue
0196                                         else:
0197                                                 hstack[subdet][module][overlap].append(ROOT.THStack("hstack",""))
0198                                                 legend[subdet][module][overlap].append(TkAlStyle.legend(len(filesTitlesColorsStyles), 1))
0199                                                 legend[subdet][module][overlap][profile].SetBorderSize(0)
0200                                                 legend[subdet][module][overlap][profile].SetFillStyle(0)    
0201         for files, title, color, style in filesTitlesColorsStyles:
0202                 h = hist(files,files.replace("/",""),subdet_ids,module_directions,overlap_directions,profile_directions)
0203                 for subdet in range(6):
0204                         for module in range(3):
0205                                 for overlap in range(3):
0206                                         if subdetCondition(subdet,module,overlap):
0207                                                         continue
0208                                         for profile in range(4):
0209                                                 if not((subdet_ids[subdet]) and  (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]):
0210                                                         continue
0211                                                 g = h[subdet][module][overlap][profile]
0212                                                 g.SetLineColor(color)
0213                                                 g.SetLineStyle(style)
0214                                                 hMean = g.GetMean(1)
0215                                                 hMeanError = g.GetMeanError(1)
0216                                                 if (profile==0):
0217                                                         legend[subdet][module][overlap][profile].AddEntry(g, title + ", mean = {0}#pm{1}#mum ".format(round(hMean,3),round(hMeanError,3)), "l")
0218                                                         g.Scale(1 / g.Integral())
0219                                                 else:
0220                                                         legend[subdet][module][overlap][profile].AddEntry(g, title, "l")
0221                                                 hstack[subdet][module][overlap][profile].Add(g)
0222         for subdet in range(6):
0223                 for module in range(3):
0224                         for overlap in range(3):
0225                                 if subdetCondition(subdet,module,overlap):
0226                                         continue
0227                                 for profile in range(4):
0228                                         if not((subdet_ids[subdet]) and  (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]):
0229                                                 continue
0230                                         currLegend = legend[subdet][module][overlap][profile]
0231                                         currhstack = hstack[subdet][module][overlap][profile]
0232                                         currhstack.SetMaximum(currhstack.GetMaximum("nostack") * 1.2)    
0233                                         c = ROOT.TCanvas()
0234                                         currhstack.Draw("hist nostack" if profile == 0 else "nostack")
0235                                         currLegend.Draw()
0236                                         xTitle = "hit_{A} - pred_{A} - (hit_{B} - pred_{B}) (#mum)"
0237                                         yTitle="fraction of events"
0238                                         save_as_file_name = file_name +  "{0}_{1}_{2}".format(dirNameList[module],dirNameList[overlap],detNameList[subdet])
0239                                         if profile>0:
0240                                                 save_as_file_name = file_name +"Profiles/profile_{0}_{1}_{2}_{3}".format(dirNameList[module],dirNameList[overlap],detNameList[subdet],dirNameList[profile-1])
0241                                                 yTitle= xTitle
0242                                                 xTitle= dirNameList[profile-1]
0243                                         currhstack.GetXaxis().SetTitle(xTitle)
0244                                         currhstack.GetYaxis().SetTitle(yTitle)
0245                                         if profile==0:
0246                                                 currhstack.GetXaxis().SetNdivisions(404)
0247                                                 
0248                                         TkAlStyle.drawStandardTitle()        
0249                                                         
0250                                         for ext in "png", "eps", "root", "pdf":
0251                                                 c.SaveAs(save_as_file_name+"." +ext)
0252