makeBrightnessMap.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 import ROOT
3 import numpy as np
4 import argparse
5 import bisect
6 import math
7 
8 def findClosestVal(val, quantList):
9  n = len( quantList )
10  i = bisect.bisect( quantList, val )
11  if i == 0:
12  return 0
13  if i == n:
14  return n-1
15  if math.fabs( val - quantList[i-1] ) < math.fabs( val - quantList[i] ):
16  return i-1
17  else:
18  return i
19 
20 if __name__=='__main__':
21 
22  parser = argparse.ArgumentParser(description='Process fiber brightness files')
23  parser.add_argument('-fd', '--far-detector', action='store_true')
24  parser.add_argument('-nd', '--near-detector', action='store_true')
25  args = parser.parse_args()
26 
27  if args.far_detector and args.near_detector:
28  print "Can only set --fd or --nd!"
29  exit(1)
30 
31  if not args.far_detector and not args.near_detector:
32  print "Must set either --fd or --nd!"
33  exit(1)
34 
35  if args.far_detector:
36  fIn = open("FDModuleData.txt", "r")
37  fOutName = "fdBrightness.root"
38  if args.near_detector:
39  fIn = open("NDModuleData.txt", "r")
40  fOutName = "ndBrightness.root"
41 
42  planeDict = {}
43  plane_moduleDict = {}
44  brightnessLst = []
45  grBright = ROOT.TGraph()
46 
47  maxModule = 0
48  maxCell = 0
49  maxPlane = 0
50 
51  maxBrightness = -9999
52  minBrightness = 9999
53  for line in fIn:
54  splitLine = line.split()
55  plane = int(splitLine[0])
56  cell = int(splitLine[1])
57  brightness = float(splitLine[2])
58  module = cell/32
59 
60  if brightness > 20 and brightness < 180:
61  brightnessLst.append(brightness/100.)
62 
63  if plane > maxPlane:
64  maxPlane = plane
65  if cell > maxCell:
66  maxCell = cell
67  if module > maxModule:
68  maxModule = module
69  if brightness < minBrightness:
70  minBrightness = brightness
71  if brightness > maxBrightness:
72  maxBrightness = brightness
73 
74  print "min brightness:", minBrightness, "max brightness:", maxBrightness
75 
76  brightness1D = ROOT.TH1D("brightness1D", "Brightness Level;Brightness;Cells", 1200, -1, 5)
77  brightnessValue = ROOT.TH1D("brightnessValue", "Brightness Level Bin Values;Brightness Bin Number;Values", 9, 0, 9)
78  brightness2D = ROOT.TH2D("brightness2D", ";Plane;Cell", maxPlane+1, 0, maxPlane+1, maxCell+1, 0, maxCell+1)
79  brightness2DX = ROOT.TH2D("brightness2DX", ";Plane;Cell", maxPlane/2+1, 0, maxPlane/2+1, maxCell+1, 0, maxCell+1)
80  brightness2DY = ROOT.TH2D("brightness2DY", ";Plane;Cell", maxPlane/2+1, 0, maxPlane/2+1, maxCell+1, 0, maxCell+1)
81 
82  brightness2D_module = ROOT.TH2D("brightness2D_module", ";Plane;Module", maxPlane+1, 0, maxPlane+1, maxModule+1, 0, maxModule+1)
83  brightness2D_bin = ROOT.TH2C("brightness2D_bin", "Brightness Level Bin Map;Plane;Cell", maxPlane+1, 0, maxPlane+1, maxCell+1, 0, maxCell+1)
84  brightness_byplane = ROOT.TH1D("brightness_plane", ";Plane", maxPlane+1, 0, maxPlane+1);
85 
86  print "maxplane: {}".format(maxPlane)
87  print "maxcell: {}".format(maxCell)
88 
89  brightnessLst.sort()
90  brightnessArr = np.array(brightnessLst)
91 
92  #quantList = []
93  #for i in [-2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5]:
94  # quantList.append( 100*0.5*(1.0 + ROOT.TMath.Erf(i/ROOT.TMath.Sqrt(2.0))) )
95 
96  quantList = [10, 20, 30, 40, 50, 60, 70, 80, 90]
97  quantVals = np.percentile(brightnessArr, quantList)
98 
99  brightnessValue.SetBinContent(1, 1)
100  for i in range(len(quantVals)):
101  brightnessValue.SetBinContent(i+1, quantVals[i])
102 
103  for i in xrange(len(quantList)):
104  print str(quantList[i])+"%:", quantVals[i]
105 
106  for iX in xrange(brightness2D.GetXaxis().GetNbins()):
107  for iY in xrange(brightness2D.GetYaxis().GetNbins()):
108  brightness2D.SetBinContent(iX+1, iY+1, 1.0)
109 
110  for iX in xrange(brightness2D_module.GetXaxis().GetNbins()):
111  for iY in xrange(brightness2D_module.GetYaxis().GetNbins()):
112  brightness2D_module.SetBinContent(iX+1, iY+1, 1.0)
113 
114  for iX in xrange(brightness2D_bin.GetXaxis().GetNbins()):
115  for iY in xrange(brightness2D_bin.GetYaxis().GetNbins()):
116  brightness2D_bin.SetBinContent(iX+1, iY+1, 1)
117 
118  position = fIn.seek(0, 0);
119  for line in fIn:
120  splitLine = line.split()
121  plane = int(splitLine[0])
122  cell = int(splitLine[1])
123  brightness = float(splitLine[2])
124  module = cell/32
125 
126  # Fill brightness 1D
127  brightness1D.Fill(brightness/100.0)
128 
129  # Fill brightness 2D, 2DX, 2DY
130  if brightness > 20 and brightness < 180:
131  iquantBrightness = findClosestVal(brightness/100., quantVals)
132  quantBrightness = quantVals[iquantBrightness]
133  brightness2D.SetBinContent(plane+1, cell+1, quantBrightness)
134 
135  if plane%2 == 0:
136  brightness2DX.SetBinContent((plane/2)+1, cell+1, quantBrightness)
137  else:
138  brightness2DY.SetBinContent((plane/2)+1, cell+1, quantBrightness)
139  if plane in planeDict:
140  planeDict[plane][0] += 1
141  planeDict[plane][1] += brightness
142  else:
143  planeDict[plane] = [1, brightness]
144  if (plane, module) in plane_moduleDict:
145  plane_moduleDict[(plane,module)][0] += 1
146  plane_moduleDict[(plane,module)][1] += brightness
147  else:
148  plane_moduleDict[(plane,module)] = [1, brightness]
149 
150  # Fill brightness by plane
151  for key in planeDict:
152  grBright.SetPoint(key, key, planeDict[key][1]/planeDict[key][0])
153  brightness_byplane.SetBinContent(key+1, planeDict[key][1]/(100.0*planeDict[key][0]));
154 
155  # Fill brightness by module
156  for key in plane_moduleDict:
157  plane = key[0]
158  module = key[1]
159  count = float(plane_moduleDict[key][0])
160  totBrightness = float(plane_moduleDict[key][1])/100.
161  brightness2D_module.SetBinContent(plane+1, module+1, totBrightness/count)
162 
163  # Fill brightness by bin, set cell brightness level to module average if it's not in [20, 180]
164  position = fIn.seek(0, 0);
165  for line in fIn:
166  splitLine = line.split()
167  plane = int(splitLine[0])
168  cell = int(splitLine[1])
169  brightness = float(splitLine[2])
170  module = cell/32
171 
172  if brightness > 20 and brightness < 180:
173  iquantBrightness = findClosestVal(brightness/100., quantVals)
174  brightness2D_bin.SetBinContent(plane+1, cell+1, iquantBrightness+1)
175  else:
176  iquantBrightness = findClosestVal(brightness_byplane.GetBinContent(plane+1, module+1), quantVals)
177  brightness2D_bin.SetBinContent(plane+1, cell+1, iquantBrightness+1)
178 
179  fOut = ROOT.TFile(fOutName, "recreate")
180  fOut.WriteTObject(brightness2D, "BrightnessByCell")
181  fOut.WriteTObject(brightnessValue, "BrightnessValue")
182  fOut.WriteTObject(brightness2D_module, "BrightnessByModule")
183  fOut.WriteTObject(brightness_byplane, "BrightnessByPlane")
184  fOut.WriteTObject(brightness2D_bin, "BrightnessByBin")
def findClosestVal(val, quantList)
std::string format(const int32_t &value, const int &ndigits=8)
Definition: HexUtils.cpp:14
procfile open("FD_BRL_v0.txt")
exit(0)