1
# Author: Sebastien Viret <viret@in2p3.fr>
2
# Python Porter: Christopher Tunnell <tunnell@hep.uchicago.edu>
6
# This macro compares two LASER runs and plots a variation map
10
from src.GenericWorker import *
11
from src.region import *
16
class Stability(GenericWorker):
17
"This worker computes the stability of laser constants between two runs"
23
# make this runs eventually
24
def __init__(self, run_type, referenceRun, mainRun):
25
self.run_type = run_type
26
self.referenceRun = referenceRun
27
self.mainRun = mainRun
30
# Get the channel position (SV: how to get that from region.py???)
31
def GetNumber(self, hash):
32
split = hash.split('_')[1:]
35
if split[0] == 'LBA': number.append(1)
36
elif split[0] == 'LBC': number.append(2)
37
elif split[0] == 'EBA': number.append(3)
38
elif split[0] == 'EBC': number.append(4)
42
number.append(int(split[1][1:]))
44
number.append(int(split[2][1:]))
46
if split[3] == 'lowgain': number.append(0)
47
elif split[3] == 'highgain': number.append(1)
54
# Reorder the PMTs (SV: how to get that from region.py???)
56
def get_PMT_index(self,part,j):
58
chan2PMT_LB=[1,2,3,4,5,6,7,8,9,10,
59
11,12,13,14,15,16,17,18,19,20,
60
21,22,23,24,27,26,25,30,29,28,
61
-33,-32,31,36,35,34,39,38,37,42,41,
62
40,45,-44,43,48,47,46]
64
chan2PMT_EB=[1,2,3,4,5,6,7,8,9,10,
65
11,12,13,14,15,16,17,18,-19,-20,
66
21,22,23,24,-25,-26,-27,-28,-31,-32,
67
33,29,30,-35,-36,34,44,38,37,43,42,
68
41,-39,-40,-45,-46,-47,-48]
71
chan = chan2PMT_LB[j]-1
73
chan = chan2PMT_EB[j]-1
77
# Macro necessary for the fiber re-ordering
78
def get_index(self, part, module, channel):
81
return 2*module+part+2*(0.5-part)*(channel%2)
83
return 2*module+(channel%2)
85
# Also for fiber stuff
86
def get_part(self, part):
93
def ProcessRegion(self, region):
97
# This macro just works for TileCal
99
if 'TILECAL' == region.GetHash():
103
# First got all the LASER events
105
for event in region.RecursiveGetEvents():
106
if event.runType == self.run_type: # Just collect the LASER stuff
109
# Only process events we should, but we return main_events + other_events + reference
111
reference_events = set()
114
channel_s_brut = ROOT.TH2F("channels variation map","",128,0.,128.,96,0.,96.)
115
channel_s_brut1 = ROOT.TH1F("variation summary","",200,-2.,2.)
117
# Initialization of variables
119
n_channels_tot = 9856 # Total number of LASER channels (per gain)
120
n_evts = 1000 # Minimum number of events recorded for a channel
121
n_sig = 2 # How many sigmas do we tolerate for filter correction (for outlier rejection ONLY) ?
122
n_sig2 = 2 # How many sigmas do we tolerate for fiber correction (for outlier rejection ONLY) ?
124
isBadLimit = 20 # Variation which put the channel automatically bad, after filter correction tough
125
isProblem = 2 # Maximum accpeted variation (after all corrections)
129
rat_ref = [0 for x in range(12288)]
133
n_FIL = [0 for x in range(4)]
134
e_FIL = [0 for x in range(4)]
135
e_FIL_n = [0 for x in range(4)]
136
sig_FIL = [0 for x in range(4)]
137
sig_FIL_n = [100000 for x in range(4)]
139
# Patch panel fibers shifts
141
n_FIB = [0 for x in range(384)]
142
e_FIB = [0 for x in range(384)]
143
e_FIB_n = [0 for x in range(384)]
144
sig_FIB = [0 for x in range(384)]
145
sig_FIB_n = [100000 for x in range(384)]
147
m_gain = 0 # Low gain by default
150
# Start by getting the variation
162
if event.runNumber == self.mainRun:
163
[p, i, j, w] = self.GetNumber(event.data['region'])
164
if self.get_PMT_index(p-1,j) >= 0: # a good channel
165
main_events.add(event)
166
filter_main = event.data['wheelpos']
167
amp_main = event.data['requamp']
168
elif event.runNumber == self.referenceRun:
169
[p, i, j, w] = self.GetNumber(event.data['region'])
170
if self.get_PMT_index(p-1,j) >= 0: # a good channel
171
index = 3072*(p-1) + 48*(i-1) + j
172
rat_ref[index] = event.data['calib_const']
173
filter_ref = event.data['wheelpos']
174
amp_ref = event.data['requamp']
176
other_events.add(event) # las
183
if filter_main != filter_ref or amp_main != amp_ref:
184
print 'Run don\'t have the same characteristics: abort!'
187
# Then compute the variation w.r.t. the reference run
188
for event in main_events:
189
[p, i, j, w] = self.GetNumber(event.data['region'])
190
index = 3072*(p-1) + 48*(i-1) + j
191
if rat_ref[index] == 0: # No reference data, skip the event
195
event.data['variation'] = 200*(event.data['calib_const']-rat_ref[index])/(event.data['calib_const']+rat_ref[index])
196
n_channels = n_channels + 1
198
main_events = main_events - to_remove
200
#################################################
202
## LOOP 1 : run over all the channels just to
203
## get the overall variation due to the
204
## filter wheel (one per partition)
206
## We do a certain number of iterations here
208
#################################################
211
for iter in range(4): # Iteration
213
n_FIL = [0 for x in range(4)]
214
e_FIL = [0 for x in range(4)]
215
sig_FIL = [0 for x in range(4)]
217
for event in main_events:
219
[part_num, i, j, w] = self.GetNumber(event.data['region'])
222
var = event.data['variation'] - e_FIL_n[part_num]
225
if iter>0 and math.fabs(var)>isBadLimit:
228
if math.fabs(var) < n_sig*sig_FIL_n[part_num]:
230
e_FIL[part_num] += var
237
for event in main_events:
239
[part_num, i, j, w] = self.GetNumber(event.data['region'])
242
var = event.data['variation'] - e_FIL[part_num] - e_FIL_n[part_num]
244
if iter>0 and math.fabs(var)>isBadLimit:
247
if math.fabs(var) < n_sig*sig_FIL_n[part_num]:
248
sig_FIL[part_num] += var*var
252
#print n_channels, '/', n_channels_tot, ' channels are initially considered'
257
sig_FIL_n[i] = math.sqrt(sig_FIL[i]/(n_FIL[i]-1))
258
e_FIL_n[i] += e_FIL[i]
259
# print 'Filter wheel is introducing a global ', e_FIL_n[i], '+/-', sig_FIL_n[i], '% shift in partition ', i
262
#################################################
265
## LOOP 2 : now we will start to tackle the PP fibers
266
## discrepancies. We collect the mean calibration
267
## value for all the fibers, corrected from the global shift
268
## Once again we will need to do some outlier removal. So we compute the
269
## sigma for all the fibers, and perform some iterations
271
#################################################
274
for iter in range(4): # Iteration
276
n_FIB = [0 for x in range(384)]
277
e_FIB = [0 for x in range(384)]
278
sig_FIB = [0 for x in range(384)]
280
for event in main_events:
282
[part_num, i, j, w] = self.GetNumber(event.data['region'])
284
j = self.get_PMT_index(part_num,j)
286
index = self.get_index(part_num,i-1,j)
287
part = self.get_part(part_num)
288
indice = int(128*part+index)
290
var = event.data['variation'] - e_FIL_n[part_num]
292
if (math.fabs(var) > isBadLimit):
295
var = var - e_FIB_n[indice]
298
if math.fabs(var) < n_sig2*sig_FIB_n[indice]:
307
for event in main_events:
309
[part_num, i, j, w] = self.GetNumber(event.data['region'])
310
#[part_num, i, j, w] = event.region.GetNumber()
312
j = self.get_PMT_index(part_num,j)
314
index = self.get_index(part_num,i-1,j)
315
part = self.get_part(part_num)
316
indice = int(128*part+index)
318
var = event.data['variation'] - e_FIL_n[part_num]
320
if (math.fabs(var) > isBadLimit):
323
var = var - e_FIB_n[indice]
325
if math.fabs(var) < n_sig2*sig_FIB_n[indice]:
326
sig_FIB[indice] += var*var
330
sig_FIB_n[i] = math.sqrt(sig_FIB[i]/(n_FIB[i]-1))
331
e_FIB_n[i] += e_FIB[i]
333
# print 'Filter wheel is introducing a global ', e_FIL_n[0], '+/-', sig_FIL_n[0], '% shift '
334
# print 'Fiber is introducing a global ', e_FIB_n[0], '+/-', sig_FIB_n[0], '% shift '
336
#################################################
339
## LOOP 3 : get the final comparison, removing the
342
#################################################
344
for event in main_events:
346
[part_num, i, j, w] = self.GetNumber(event.data['region'])
348
j = self.get_PMT_index(part_num,j)
350
index = self.get_index(part_num,i-1,j)
351
part = self.get_part(part_num)
352
indice = int(128*part+index)
354
var = event.data['variation'] - e_FIL_n[part_num] - e_FIB_n[indice]
356
channel_s_brut.Fill(i-0.5+64*(part_num%2),j+0.5+48*((part_num/2)%2),var)
357
channel_s_brut1.Fill(var)
360
## Then do the cosmetics
362
##### STAT ###### #TODO! DO NOT ACCESS gStyle directly! Hurts everybody else. fixme
363
ROOT.gStyle.SetStatW(0.2)
364
ROOT.gStyle.SetStatH(0.2)
365
ROOT.gStyle.SetStatX(0.8)
366
ROOT.gStyle.SetStatY(0.8)
367
ROOT.gStyle.SetStatFont(42)
368
ROOT.gStyle.SetOptStat(0)
369
ROOT.gStyle.SetOptFit(0)
375
Red = [ 1.0, 1.0, 1.0, 1.0, 0.0 ]
376
Green = [ 1.0, 1.0, 0.6, 0.0, 0.0 ]
377
Blue = [ 1.0, 0.0, 0.1, 0.0, 0.0 ]
379
# Define the length of the (color)-interval between this points
380
Length = [ 0.00, 0.25, 0.5, 0.75, 1.00 ]
382
ROOT.TColor.CreateGradientColorTable(Number, array('d', Length), array('d', Red),
383
array('d', Green), array('d', Blue), 300);
385
c1 = ROOT.TCanvas("c1","Filter+fiber correction",200,200,1000,700)
391
channel_s_brut.GetXaxis().SetTitle("Module number")
392
channel_s_brut.GetYaxis().SetTitle("Channel number")
393
channel_s_brut.Draw("colz")
394
pt = ROOT.TPaveText(20.,95.,110.,105.,"br")
397
filename = "Diff. bet. runs %d and %d : all partitions" % (self.mainRun, self.referenceRun)
399
text = pt.AddText(filename);
404
filename = "Diff. bet. runs %d and %d : all partitions (in %%)" % (self.mainRun, self.referenceRun)
406
channel_s_brut1.GetXaxis().SetTitle(filename)
407
channel_s_brut1.Draw()
410
c1.Print('plots/laser_comparison1.eps')
added
removed
workers/Laser/Stability.py
