File indexing completed on 2023-03-17 10:40:25
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"]
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):
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 subdetCondition(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]):
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):
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
0074 if subdet_id%2 == 1 and ((abs(t.moduleZ[0] - t.moduleZ[1]) > 1)):
0075 module_direction = 0
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:
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):
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
