Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-11-25 02:29:06 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!