1
|
#!/usr/bin/env python
|
2
|
|
3
|
import libMausCpp
|
4
|
import ROOT
|
5
|
from ROOT import gROOT
|
6
|
import os
|
7
|
from os import listdir
|
8
|
import numpy as np
|
9
|
from math import floor
|
10
|
import _ctypes
|
11
|
|
12
|
class Analysis:
|
13
|
def __init__ (self):
|
14
|
|
15
|
# defines where to cut on bad channel noise
|
16
|
self.bc_noise_cut = 2.0 # twice the expected value
|
17
|
self.bc_quiet_cut = 0.5 # half the expected value
|
18
|
self.event_cut = -1 # number of events to process, -1 for all events
|
19
|
|
20
|
# Designed mostly for testing purposes, when turned on allows much quicker
|
21
|
# loading times on data that has already been analyized once.
|
22
|
self.load_previous = False
|
23
|
self.save_record = True
|
24
|
self.previous_filename = "bc_history.root"
|
25
|
|
26
|
# Maus data location
|
27
|
self.run_number = "7417"
|
28
|
self.data_directory = "/vols/fets2/heidt/offline/data/7417/"
|
29
|
self.scifi_map_file = "scifi_mapping_2015-07-28.txt"
|
30
|
self.scifi_map_directory = os.environ["MAUS_ROOT_DIR"]+"/files/cabling/"
|
31
|
|
32
|
# Analysis
|
33
|
self.process()
|
34
|
# Output
|
35
|
self.Output()
|
36
|
|
37
|
##################################################################################################
|
38
|
# Formats the data and calls Bad_Channels() which does the analysis
|
39
|
def process(self):
|
40
|
self.Load_mapp()
|
41
|
|
42
|
self.unit_hist = ROOT.TH1D("unit","unit_hist_longr",212,0,212)
|
43
|
for i in range(212):
|
44
|
self.unit_hist.SetBinContent(i,1)
|
45
|
|
46
|
if self.load_previous == True:
|
47
|
self.prev_file = ROOT.TFile(self.previous_filename, "READ")
|
48
|
self.Recreate_ROOT()
|
49
|
|
50
|
if self.load_previous == False:
|
51
|
# Creates the container ROOT files and loads up the MAUS processed file.
|
52
|
self.Make_ROOT()
|
53
|
file_in = self.Load_file()
|
54
|
print "Reading MAUS processed file: ",file_in
|
55
|
root_file = ROOT.TFile(file_in, "READ")
|
56
|
# Checks spill/event is good data
|
57
|
data = ROOT.MAUS.Data()
|
58
|
tree = root_file.Get("Spill")
|
59
|
if not tree:
|
60
|
return
|
61
|
tree.SetBranchAddress("data", data)
|
62
|
print "Filling Histogram"
|
63
|
peat_count = 0
|
64
|
if self.event_cut > tree.GetEntries() or self.event_cut <= 0:
|
65
|
self.event_cut = tree.GetEntries()
|
66
|
for i in range(self.event_cut):
|
67
|
peat_count += 1
|
68
|
if (peat_count % 100 == 0):
|
69
|
print "Filling event: ",peat_count, "/", self.event_cut
|
70
|
tree.GetEntry(i)
|
71
|
self.spill = data.GetSpill()
|
72
|
if not self.spill.GetDaqEventType() == "physics_event":
|
73
|
continue
|
74
|
if self.spill.GetReconEvents().size() == 0:
|
75
|
print "No Recon Events"
|
76
|
continue
|
77
|
|
78
|
# Fills the ROOT containers with data from the MAUS files
|
79
|
self.Fill_Hists()
|
80
|
|
81
|
check_out = open("check_out.txt", "w")
|
82
|
for tra in range(0,2):
|
83
|
for sta in range(1,6):
|
84
|
for pla in range(0,3):
|
85
|
for bin in range(len(self.dig_cont[tra][sta][pla])):
|
86
|
check_out.write("Tracker: ")
|
87
|
check_out.write(str(tra))
|
88
|
check_out.write(" Station: ")
|
89
|
check_out.write(str(sta))
|
90
|
check_out.write(" Plane: ")
|
91
|
check_out.write(str(pla))
|
92
|
check_out.write(" Bin: ")
|
93
|
check_out.write(str(bin))
|
94
|
check_out.write(" = ")
|
95
|
check_out.write(str(self.dig_cont[tra][sta][pla].GetBinContent(bin)))
|
96
|
check_out.write("\n")
|
97
|
check_out.close()
|
98
|
|
99
|
print "Running bad channel analysis"
|
100
|
self.Bad_Channels()
|
101
|
|
102
|
##################################################################################################
|
103
|
# Pulls out tracker digits and sorts them into histos sorted by plane/station/tracker
|
104
|
def Fill_Hists(self):
|
105
|
for rc in range(self.spill.GetReconEvents().size()):
|
106
|
digit = self.spill.GetReconEvents()[rc].GetSciFiEvent().digits()
|
107
|
for di in range(len(digit)):
|
108
|
trac = stat = chan = plan = nupe = -1
|
109
|
trac = digit[di].get_tracker()
|
110
|
stat = digit[di].get_station()
|
111
|
chan = digit[di].get_channel()
|
112
|
plan = digit[di].get_plane()
|
113
|
nupe = digit[di].get_npe()
|
114
|
cont = self.dig_cont[trac][stat][plan].GetBinContent(chan)
|
115
|
self.dig_cont[trac][stat][plan].SetBinContent(chan, cont+1)
|
116
|
self.Record[trac][stat][plan].SetBinContent(chan, cont+1)
|
117
|
|
118
|
##################################################################################################
|
119
|
# loads tracker mapping, used to call back to bank/board for bad channel list
|
120
|
# used by Map_to_VLSB()
|
121
|
def Load_mapp(self):
|
122
|
print "Loading Map"
|
123
|
mapfi = open(self.scifi_map_directory + self.scifi_map_file,'r')
|
124
|
self.map=[]
|
125
|
for line in mapfi:
|
126
|
line = line.strip('\n')
|
127
|
self.map.append(line.split(" "))
|
128
|
|
129
|
##################################################################################################
|
130
|
# Takes the zero channels from all previous functions and writes them out
|
131
|
# to a bad channels list.
|
132
|
def Dead_Channels(self):
|
133
|
self.dig_emp_list = []
|
134
|
for tra in range(0,2):
|
135
|
for sta in range(1,6):
|
136
|
for pla in range(0,3):
|
137
|
for bin in range(len(self.dig_cont[tra][sta][pla])):
|
138
|
if self.dig_cont[tra][sta][pla].GetBinContent(bin) == 0:
|
139
|
BB = self.Map_to_VLSB(tra,sta,pla,bin)
|
140
|
if BB != [-1,-1]:
|
141
|
self.dig_emp_list.append(BB)
|
142
|
|
143
|
##################################################################################################
|
144
|
# Uses map to determine board and bank info
|
145
|
def Map_to_VLSB(self,tracker,station,plane,channel):
|
146
|
board = channel_ro = -1
|
147
|
for mp in range(len(self.map)):
|
148
|
if int(self.map[mp][3])==tracker and int(self.map[mp][4])==station and \
|
149
|
int(self.map[mp][5])==plane and int(self.map[mp][6])==channel:
|
150
|
board=int(self.map[mp][0])*4+int(self.map[mp][1])
|
151
|
channel_ro=int(self.map[mp][2])
|
152
|
continue
|
153
|
if board == -1 or channel_ro == -1:
|
154
|
print "NOT FOUND! ",tracker," ",station," ",plane," ",channel
|
155
|
VLSB=[board,channel_ro]
|
156
|
return VLSB
|
157
|
|
158
|
##################################################################################################
|
159
|
# Calls the recursive functions Hot_Channels() and Shh_Channels() which look
|
160
|
# for hot and quite channels respectively.
|
161
|
def Bad_Channels(self):
|
162
|
for tra in self.dig_cont:
|
163
|
for sta in self.dig_cont[tra]:
|
164
|
for pla in self.dig_cont[tra][sta]:
|
165
|
# print "Bad channels analysis in: tracker ", tra," station ",sta," plane ",pla
|
166
|
fit_hist = self.dig_cont[tra][sta][pla]
|
167
|
self.Hot_Channels(fit_hist)
|
168
|
self.Shh_Channels(fit_hist)
|
169
|
self.dig_cont[tra][sta][pla] = fit_hist
|
170
|
self.Dead_Channels()
|
171
|
|
172
|
###################################################################################################
|
173
|
# Attempts to fit various functions to digit counts
|
174
|
def Fit_Hist(self, fit_hist):
|
175
|
status=fit_hist.Fit("gaus","SQ")
|
176
|
fit = "gaus"
|
177
|
if status.CovMatrixStatus() < 3:
|
178
|
# print "Problem with fitting, attempting a linear fit"
|
179
|
fit_hist.GetListOfFunctions().Clear()
|
180
|
status = fit_hist.Fit("pol1","SQ")
|
181
|
# fit_hist.Draw()
|
182
|
# raw_input("Press Enter to Exit")
|
183
|
fit = "pol1"
|
184
|
return fit
|
185
|
|
186
|
##################################################################################################
|
187
|
# Looks at channels/fit ratio and cuts above bc_noise_cut criteria then sets bin content
|
188
|
# equal to zero. Bad bins are read out in Dead_Channels()
|
189
|
def Hot_Channels(self, fit_hist):
|
190
|
fit = self.Fit_Hist(fit_hist)
|
191
|
hdiff=fit_hist.Clone("hdiff")
|
192
|
hdiff.Divide(hdiff.GetFunction(fit))
|
193
|
hdiff.GetListOfFunctions().Clear()
|
194
|
bin = hdiff.GetMaximumBin()
|
195
|
max = hdiff.GetBinContent(int(bin))
|
196
|
if (max > self.bc_noise_cut):
|
197
|
# The lines below exists for testing purposes only.
|
198
|
# fit_hist.Draw()
|
199
|
# raw_input("Press Enter to Exit")
|
200
|
fit_hist.SetBinContent(bin,0)
|
201
|
self.Hot_Channels(fit_hist)
|
202
|
|
203
|
##################################################################################################
|
204
|
# Looks at channel/fit ratio and cuts below bc_quite_cut criteria then sets bin content
|
205
|
# equal to zero. Bad bins are read out in Dead_Channels()
|
206
|
def Shh_Channels(self, fit_hist):
|
207
|
fit = self.Fit_Hist(fit_hist)
|
208
|
hdiff=fit_hist.Clone("hdiff")
|
209
|
hdiff.Divide(hdiff.GetFunction(fit))
|
210
|
hdiff.GetListOfFunctions().Clear()
|
211
|
hdiff.Divide(self.unit_hist,hdiff)
|
212
|
bin = hdiff.GetMaximumBin()
|
213
|
min = hdiff.GetBinContent(int(bin))
|
214
|
if (min > 1.0/self.bc_quiet_cut):
|
215
|
# The lines below exists for testing purposes only.
|
216
|
# fit_hist.Draw()
|
217
|
# raw_input("Press Enter to Exit")
|
218
|
fit_hist.SetBinContent(bin,0)
|
219
|
self.Shh_Channels(fit_hist)
|
220
|
|
221
|
##################################################################################################
|
222
|
# Creates the container root file that will be passed around script.
|
223
|
def Make_ROOT(self):
|
224
|
print "Creating empty ROOT file"
|
225
|
self.dig_cont={}
|
226
|
self.Record={}
|
227
|
for tra in range(0,2):
|
228
|
self.dig_cont[tra]={}
|
229
|
self.Record[tra]={}
|
230
|
if tra == 0:
|
231
|
trs = "U"
|
232
|
else:
|
233
|
trs = "D"
|
234
|
for sta in range(1,6):
|
235
|
self.dig_cont[tra][sta]={}
|
236
|
self.Record[tra][sta]={}
|
237
|
for pla in range(0,3):
|
238
|
if pla == 0 or pla == 1 or pla == 2:
|
239
|
dig_name="DC_Tk%s%d%d" %(trs,sta,pla)
|
240
|
dig_titl="Digit Scan Tk%s S%d P%d" %(trs,sta,pla)
|
241
|
self.dig_cont[tra][sta][pla]=ROOT.TH1D(dig_name,dig_titl,212,0,212)
|
242
|
|
243
|
dig_name="Rec_Tk%s%d%d" %(trs,sta,pla)
|
244
|
dig_titl="Digit Scan Tk%s S%d P%d" %(trs,sta,pla)
|
245
|
self.Record[tra][sta][pla]=ROOT.TH1D(dig_name,dig_titl,212,0,212)
|
246
|
|
247
|
##################################################################################################
|
248
|
# Loads a previously generated root file.
|
249
|
def Recreate_ROOT(self):
|
250
|
print "Loading previously used ROOT file"
|
251
|
self.dig_cont={}
|
252
|
for tra in range(0,2):
|
253
|
self.dig_cont[tra]={}
|
254
|
if tra == 0:
|
255
|
trs = "U"
|
256
|
else:
|
257
|
trs = "D"
|
258
|
for sta in range(1,6):
|
259
|
self.dig_cont[tra][sta]={}
|
260
|
for pla in range(0,3):
|
261
|
dig_name="Rec_Tk%s%d%d" %(trs,sta,pla)
|
262
|
self.dig_cont[tra][sta][pla]=self.prev_file.Get(dig_name)
|
263
|
|
264
|
##################################################################################################
|
265
|
# Outputter
|
266
|
def Output(self):
|
267
|
out_root = ROOT.TFile("bc_out_test.root",'RECREATE')
|
268
|
for tr in range(0,2):
|
269
|
for st in range(1,6):
|
270
|
for pl in range(0,3):
|
271
|
self.dig_cont[tr][st][pl].Write()
|
272
|
out_root.Close()
|
273
|
|
274
|
if self.save_record == True:
|
275
|
out_record = ROOT.TFile(self.previous_filename,'RECREATE')
|
276
|
for tr in range(0,2):
|
277
|
for st in range(1,6):
|
278
|
for pl in range(0,3):
|
279
|
self.Record[tr][st][pl].Write()
|
280
|
out_record.Close()
|
281
|
|
282
|
dead_list = open("dead_channel_list.txt","w")
|
283
|
for dl in range(len(self.dig_emp_list)):
|
284
|
dead_list.write(str(self.dig_emp_list[dl][0]))
|
285
|
dead_list.write(" ")
|
286
|
dead_list.write(str(self.dig_emp_list[dl][1]))
|
287
|
dead_list.write("\n")
|
288
|
raw_input("Press Enter to Exit")
|
289
|
|
290
|
##################################################################################################
|
291
|
# Searches set directory and finds a processed MAUS file.
|
292
|
def Load_file(self):
|
293
|
for input_file in listdir(self.data_directory):
|
294
|
if self.run_number in input_file and ".root" in input_file:
|
295
|
file_in = self.data_directory + input_file
|
296
|
return file_in
|
297
|
|
298
|
##################################################################################################
|
299
|
|
300
|
if __name__ == "__main__":
|
301
|
Analysis()
|
302
|
|