Project

General

Profile

Analysis Code » bad_channel.py

Heidt, Christopher, 19 January 2016 13:42

 
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

    
(3-3/3)