CMS 3D CMS Logo
overlapValidationPlot.py
Go to the documentation of this file.
1 from __future__ import absolute_import
2 from __future__ import print_function
3 
4 from builtins import range
5 import math
6 
7 import ROOT
8 
9 from .TkAlStyle import TkAlStyle
10 
11 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
12 detNameList = ("BPIX", "FPIX", "TIB", "TID", "TOB", "TEC")
13 def subdetCondition(subdet,module,overlap):
14  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)):
15  return True
16  else:
17  return False
18 def hist(tree_file_name, hist_name,subdet_ids,module_directions,overlap_directions,profile_directions):
19  f = ROOT.TFile(tree_file_name)
20  t = f.Get("analysis/Overlaps")
21  h = []
22  for subdet in range(6):#Creates a 4-D list of empty histograms for permutations of subdetector, overlap direction, module direction and profile direction.
23  h.append([])
24  for module in range(3):
25  h[subdet].append([])
26  for overlap in range(3):
27  h[subdet][module].append([])
28  for profile in range(4):
29  if subdetConditon(subdet,module,overlap):
30  h[subdet][module][overlap].append(0)
31  continue
32  name = hist_name + "{0}_{1}_{2}".format(dirNameList[module],dirNameList[overlap],detNameList[subdet])
33  if not ((subdet_ids[subdet]) and (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]): #puts a 0 in any unwanted plots
34  h[subdet][module][overlap].append(0)
35  continue
36  if (profile>0):
37  if profile == 3: bins, xmin, xmax = 10, -math.pi, math.pi
38  if subdet+1 == 1 and profile == 1: bins, xmin, xmax = 10, -30, -30
39  if subdet+1 == 2 and profile == 1: bins, xmin, xmax = 40, -60, 60
40  if subdet+1 == 3 and profile == 1: bins, xmin, xmax = 10, -70, 70
41  if subdet+1 == 4 and profile == 1: bins, xmin, xmax = 40, -110, 110
42  if subdet+1 == 5 and profile == 1: bins, xmin, xmax = 10, -110, 110
43  if subdet+1 == 6 and profile == 1: bins, xmin, xmax = 20, -280, 280
44  if subdet+1 == 1 and profile == 2: bins, xmin, xmax = 10, 0, 20
45  if subdet+1 == 2 and profile == 2: bins, xmin, xmax = 10, 0, 20
46  if subdet+1 == 3 and profile == 2: bins, xmin, xmax = 10, 20, 55
47  if subdet+1 == 4 and profile == 2: bins, xmin, xmax = 10, 20, 55
48  if subdet+1 == 5 and profile == 2: bins, xmin, xmax = 10, 55, 110
49  if subdet+1 == 6 and profile == 2: bins, xmin, xmax = 10, 20, 110
50  h[subdet][module][overlap].append(ROOT.TProfile(hist_name, hist_name, bins, xmin, xmax))
51  elif subdet+1==4 or subdet+1==6:
52  h[subdet][module][overlap].append(ROOT.TH1F(name, name, 100, -5000, 5000))
53  else:
54  h[subdet][module][overlap].append(ROOT.TH1F(name, name, 100, -300, 300))
55  h[subdet][module][overlap][profile].SetDirectory(0)
56 
57  nentries = t.GetEntries()
58 
59  for i, entry in enumerate(t, start=1):#loops through the tree, filling in relevant histograms for each event
60  if i % 10000 == 0 or i == nentries:
61  print(i, "/", nentries)
62 
63  subdet_id = t.subdetID
64  if subdet_ids[subdet_id-1]==False:
65  continue
66 
67  modulePhi0 = math.atan2(t.moduleY[0], t.moduleX[0])
68  modulePhi1 = math.atan2(t.moduleY[1], t.moduleX[1])
69  phidiff = min(abs(modulePhi0-modulePhi1), abs(math.pi - abs(modulePhi0-modulePhi1)))
70  moduleR0 = math.sqrt(t.moduleY[0]**2+t.moduleX[0]**2)
71  moduleR1 = math.sqrt(t.moduleY[1]**2+t.moduleX[1]**2)
72 
73  #determines the direction the modules are in respect to each other for each event
74  if subdet_id%2 == 1 and ((abs(t.moduleZ[0] - t.moduleZ[1]) > 1)):
75  module_direction = 0 #0 refers to Z, 1 to R and 2 to phi
76  elif subdet_id%2 == 0 and (abs(moduleR0 - moduleR1)>1):
77  module_direction = 1
78  elif ((moduleR0*phidiff)>1):
79  module_direction = 2
80  else: #prints if this method of selecting module_direction misses anything
81  print(str(i)+" skipped")
82  continue
83  if module_directions[module_direction]==False:
84  continue
85  avgList=[(t.moduleZ[0]+t.moduleZ[1])/2,(moduleR0+moduleR1)/2,(modulePhi0+modulePhi1)/2]
86  if overlap_directions[2]:
87  overlap_direction = 2
88  if modulePhi0 > modulePhi1:
89  hitXA = t.hitX[1]
90  hitXB = t.hitX[0]
91  predXA = t.predX[1]
92  predXB = t.predX[0]
93  overlapSignA = t.localxdotglobalphi[1]
94  overlapSignB = t.localxdotglobalphi[0]
95  else:
96  hitXA = t.hitX[0]
97  hitXB = t.hitX[1]
98  predXA = t.predX[0]
99  predXB = t.predX[1]
100  overlapSignA = t.localxdotglobalphi[0]
101  overlapSignB = t.localxdotglobalphi[1]
102  residualA = hitXA - predXA
103  residualB = hitXB - predXB
104  if overlapSignA < 0:
105  residualA *= -1
106  if overlapSignB < 0:
107  residualB *= -1
108  A = 10000*(residualA - residualB)
109  h[subdet_id-1][module_direction][overlap_direction][0].Fill(A)
110  for profile in range(3):
111  if profile_directions[profile+1]:
112  h[subdet_id-1][module_direction][overlap_direction][profile+1].Fill(avgList[profile],A)
113 
114 
115  if subdet_id==1 and module_direction != 1 and overlap_directions[0]:
116  overlap_direction = 0
117  if t.moduleZ[0] > t.moduleZ[1]:
118  hitXA = t.hitY[1]
119  hitXB = t.hitY[0]
120  predXA = t.predY[1]
121  predXB = t.predY[0]
122  overlapSignA = t.localydotglobalz[1]
123  overlapSignB = t.localydotglobalz[0]
124  else:
125  hitXA = t.hitY[0]
126  hitXB = t.hitY[1]
127  predXA = t.predY[0]
128  predXB = t.predY[1]
129  overlapSignA = t.localydotglobalz[0]
130  overlapSignB = t.localydotglobalz[1]
131  residualA = hitXA - predXA
132  residualB = hitXB - predXB
133  if overlapSignA < 0:
134  residualA *= -1
135  if overlapSignB < 0:
136  residualB *= -1
137  A = 10000*(residualA - residualB)
138  h[subdet_id-1][module_direction][overlap_direction][0].Fill(A)
139  for profile in range(3):
140  if profile_directions[profile+1]:
141  h[subdet_id-1][module_direction][overlap_direction][profile+1].Fill(avgList[profile],A)
142 
143  if subdet_id==2 and module_direction !=0 and overlap_directions[1]:
144  overlap_direction = 1
145  if moduleR0 > moduleR1:
146  hitXA = t.hitY[1]
147  hitXB = t.hitY[0]
148  predXA = t.predY[1]
149  predXB = t.predY[0]
150  overlapSignA = t.localydotglobalr[1]
151  overlapSignB = t.localydotglobalr[0]
152  else:
153  hitXA = t.hitY[0]
154  hitXB = t.hitY[1]
155  predXA = t.predY[0]
156  predXB = t.predY[1]
157  overlapSignA = t.localydotglobalr[0]
158  overlapSignB = t.localydotglobalr[1]
159 
160  residualA = hitXA - predXA
161  residualB = hitXB - predXB
162  if overlapSignA < 0:
163  residualA *= -1
164  if overlapSignB < 0:
165  residualB *= -1
166  A = 10000*(residualA - residualB)
167  h[subdet_id-1][module_direction][overlap_direction][0].Fill(A)
168  for profile in range(3):
169  if profile_directions[profile+1]:
170  h[subdet_id-1][module_direction][overlap_direction][profile+1].Fill(avgList[profile],A)
171  return h
172 
173 def plot(file_name,subdet_ids,module_directions,overlap_directions,profile_directions,*filesTitlesColorsStyles):
174 
175  legend=[]
176  hstack=[]
177 
178  for subdet in range(6):#creates lists of empty THStacks and legends to be filled later
179  hstack.append([])
180  legend.append([])
181  for module in range(3):
182  hstack[subdet].append([])
183  legend[subdet].append([])
184  for overlap in range(3):
185  hstack[subdet][module].append([])
186  legend[subdet][module].append([])
187  for profile in range(4):
188  if subdetCondition(subdet,module,overlap):
189  hstack[subdet][module][overlap].append(0)
190  legend[subdet][module][overlap].append(0)
191  continue
192  if not ((subdet_ids[subdet]) and (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]):
193  legend[subdet][module][overlap].append(0)
194  hstack[subdet][module][overlap].append(0)
195  continue
196  else:
197  hstack[subdet][module][overlap].append(ROOT.THStack("hstack",""))
198  legend[subdet][module][overlap].append(TkAlStyle.legend(len(filesTitlesColorsStyles), 1))
199  legend[subdet][module][overlap][profile].SetBorderSize(0)
200  legend[subdet][module][overlap][profile].SetFillStyle(0)
201  for files, title, color, style in filesTitlesColorsStyles:
202  h = hist(files,files.replace("/",""),subdet_ids,module_directions,overlap_directions,profile_directions)
203  for subdet in range(6):
204  for module in range(3):
205  for overlap in range(3):
206  if subdetCondition(subdet,module,overlap):
207  continue
208  for profile in range(4):
209  if not((subdet_ids[subdet]) and (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]):
210  continue
211  g = h[subdet][module][overlap][profile]
212  g.SetLineColor(color)
213  g.SetLineStyle(style)
214  hMean = g.GetMean(1)
215  hMeanError = g.GetMeanError(1)
216  if (profile==0):
217  legend[subdet][module][overlap][profile].AddEntry(g, title + ", mean = {0}#pm{1}#mum ".format(round(hMean,3),round(hMeanError,3)), "l")
218  g.Scale(1 / g.Integral())
219  else:
220  legend[subdet][module][overlap][profile].AddEntry(g, title, "l")
221  hstack[subdet][module][overlap][profile].Add(g)
222  for subdet in range(6):
223  for module in range(3):
224  for overlap in range(3):
225  if subdetCondition(subdet,module,overlap):
226  continue
227  for profile in range(4):
228  if not((subdet_ids[subdet]) and (module_directions[module]) and (profile_directions[profile]) and overlap_directions[overlap]):
229  continue
230  currLegend = legend[subdet][module][overlap][profile]
231  currhstack = hstack[subdet][module][overlap][profile]
232  currhstack.SetMaximum(currhstack.GetMaximum("nostack") * 1.2)
233  c = ROOT.TCanvas()
234  currhstack.Draw("hist nostack" if profile == 0 else "nostack")
235  currLegend.Draw()
236  xTitle = "hit_{A} - pred_{A} - (hit_{B} - pred_{B}) (#mum)"
237  yTitle="fraction of events"
238  save_as_file_name = file_name + "{0}_{1}_{2}".format(dirNameList[module],dirNameList[overlap],detNameList[subdet])
239  if profile>0:
240  save_as_file_name = file_name +"Profiles/profile_{0}_{1}_{2}_{3}".format(dirNameList[module],dirNameList[overlap],detNameList[subdet],dirNameList[profile-1])
241  yTitle= xTitle
242  xTitle= dirNameList[profile-1]
243  currhstack.GetXaxis().SetTitle(xTitle)
244  currhstack.GetYaxis().SetTitle(yTitle)
245  if profile==0:
246  currhstack.GetXaxis().SetNdivisions(404)
247 
248  TkAlStyle.drawStandardTitle()
249 
250  for ext in "png", "eps", "root", "pdf":
251  c.SaveAs(save_as_file_name+"." +ext)
252 
FastTimerService_cff.range
range
Definition: FastTimerService_cff.py:34
min
T min(T a, T b)
Definition: MathUtil.h:58
overlapValidationPlot.plot
def plot(file_name, subdet_ids, module_directions, overlap_directions, profile_directions, *filesTitlesColorsStyles)
Definition: overlapValidationPlot.py:173
str
#define str(s)
Definition: TestProcessor.cc:51
TkAlStyle::legend
static TLegend * legend(const int nEntries, const double relWidth=0.5)
Definition: TkAlStyle.cc:109
print
void print(TMatrixD &m, const char *label=nullptr, bool mathematicaFormat=false)
Definition: Utilities.cc:46
mps_setup.append
append
Definition: mps_setup.py:85
overlapValidationPlot.subdetCondition
def subdetCondition(subdet, module, overlap)
Definition: overlapValidationPlot.py:13
overlapValidationPlot.hist
def hist(tree_file_name, hist_name, subdet_ids, module_directions, overlap_directions, profile_directions)
Definition: overlapValidationPlot.py:18
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:1053
TkAlStyle::drawStandardTitle
static void drawStandardTitle()
Definition: TkAlStyle.cc:75
format
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Generated for CMSSW Reference Manual by doxygen 1.8.16