Skip to content

Commit 382f46c

Browse files
ocerriocerri
committed
Space re-weighting and Tiles space dependence
1 parent d8a6acb commit 382f46c

File tree

5 files changed

+995
-84
lines changed

5 files changed

+995
-84
lines changed

‎Analysis_SiPM_Tile.py‎

Lines changed: 58 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66

77
import ROOT as rt
88
from root_numpy import tree2array, tree2rec
9-
from lib.histo_utilities import create_TH1D, create_TH2D, quantile
9+
from lib.histo_utilities import create_TH1D, create_TH2D, quantile, rootTH2_to_np
1010
from lib.cebefo_style import cebefo_style, Set_2D_colz_graphics
1111

1212
donotdelete = []
@@ -20,6 +20,7 @@ def parsing():
2020
parser.add_argument("-S", "--save_loc", type=str, default='./out_plots/', help="Saving location")
2121

2222
parser.add_argument("-B", "--batch", default=True, action='store_false', help="Root batch mode")
23+
parser.add_argument("-v", "--verbose", default=False, action='store_true', help="Verbose mode")
2324

2425
parser.add_argument("-N", "--runs_interval", default=None, help="Runs to run", nargs='+')
2526

@@ -130,34 +131,6 @@ def define_range_around_peak(h, perc = [0.4, 0.4], Range=[0.0, 99999.0]):
130131
x_up = h.GetBinLowEdge(n_up) + h.GetBinWidth(n_up)
131132
return x_low, x_up, n_pk
132133

133-
def rootTH2_to_np(h, cut = None, Norm = False):
134-
nx = h.GetNbinsX()
135-
ny = h.GetNbinsY()
136-
137-
arr = np.zeros((ny, nx))
138-
pos = np.zeros((ny, nx, 2))
139-
140-
for ix in range(nx):
141-
for iy in range(ny):
142-
x = h.GetXaxis().GetBinCenter( ix+1 );
143-
y = h.GetYaxis().GetBinCenter( iy+1 );
144-
z = h.GetBinContent(h.GetBin(ix+1, iy+1))
145-
146-
if cut == None:
147-
arr[iy, ix] = z
148-
else:
149-
arr[iy, ix] = z if z > cut else 0
150-
pos[iy, ix] = [x,y]
151-
return arr, pos
152-
153-
def circle_filter(h_arr, cx, cy, rx, ry):
154-
p = 0
155-
for i,j in itertools.product(np.arange(h_arr.shape[0]), np.arange(h_arr.shape[1])):
156-
if (float(cx-j)/rx)**2 + (float(cy-i)/ry)**2 < 1:
157-
p += h_arr[i-1, j-1]
158-
159-
return p/(np.pi*rx*ry)
160-
161134
def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canvas=True):
162135
q_up, e_up = quantile(1000*dt, 0.15)
163136
q_dwn, e_dwn = quantile(1000*dt, 0.85)
@@ -473,7 +446,7 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
473446
dy = configurations.xy_center[1]
474447
dx = configurations.xy_center[0]
475448
width = configurations.xy_center[2]
476-
N_bins = [100, 100]
449+
N_bins = [70, 70]
477450
if 'PosRaw' in configurations.plots:
478451
name = 'h_pos_'+str(k)
479452
title = 'Track position at z_DUT[' + str(conf['idx_dut']) + '], for channel {} selected event'.format(k)
@@ -486,6 +459,7 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
486459

487460
canvas['pos'][k] = rt.TCanvas('c_pos_'+str(k), 'c_pos_'+str(k), 800, 600)
488461
h.DrawCopy('colz')
462+
rt.gPad.SetRightMargin(0.15)
489463
canvas['pos'][k].Update()
490464
canvas['pos'][k].SaveAs(out_dir + '/PositionXY_raw_ch{:02d}.png'.format(k))
491465

@@ -600,9 +574,12 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
600574
h_w.GetZaxis().SetRangeUser(h_w.GetMinimum(10), h_w.GetMaximum())
601575
h_w.SetStats(0)
602576

603-
canvas['amp_w_pos'][k] = rt.TCanvas('c_w_pos_'+str(k), 'c_w_pos_'+str(k), 800, 600)
577+
can = rt.TCanvas('c_w_pos', 'c_w_pos', 1200, 500)
578+
if 'space_sel' in conf.keys():
579+
can.Divide(2,1)
580+
pad = can.cd(1)
604581
h_w.DrawCopy('colz')
605-
Set_2D_colz_graphics()
582+
rt.gPad.SetRightMargin(0.15)
606583
if not conf['SiPM'] == 'None':
607584
size = conf['SiPM'].split('x')
608585

@@ -613,8 +590,6 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
613590
nby = h_w.GetYaxis().FindBin(y)
614591

615592
arr_h, pos_h = rootTH2_to_np(h_w, cut=100)
616-
# print arr_h
617-
618593
p_max = 0
619594
idx_max = [0, 0]
620595
for iy in range(0, arr_h.shape[0]-nby):
@@ -627,8 +602,8 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
627602

628603
avg_prob = p_max/(nbx*nby)
629604
iy, ix = idx_max
630-
if k == 24:
631-
ix -= 2
605+
# if k == 24:
606+
# ix -= 2
632607
x_start = h.GetXaxis().GetBinCenter(ix+1) - 0.5*h.GetXaxis().GetBinWidth(ix+1)
633608
x_stop = h.GetXaxis().GetBinCenter(ix+nbx) + 0.5*h.GetXaxis().GetBinWidth(ix+nbx)
634609
y_start = h.GetYaxis().GetBinCenter(iy+1) - 0.5*h.GetYaxis().GetBinWidth(iy+1)
@@ -641,17 +616,50 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
641616
line.DrawLine(x_start, y_stop, x_stop, y_stop)
642617
line.DrawLine(x_start, y_start, x_start, y_stop)
643618
line.DrawLine(x_stop, y_start, x_stop, y_stop)
644-
l = rt.TLatex()
645-
l.SetTextColor(46)
646-
l.SetTextSize(0.02)
647-
l.DrawLatex(x_start, y_stop, 'Avg amp {:.0f} mV'.format(avg_prob))
648619

649620
conf['SiPM_sel'] = {}
650621
conf['SiPM_sel']['limits'] = [x_start, x_stop, y_start, y_stop]
651622
conf['SiPM_sel']['cut'] = '({y} < {h} && {y} > {l}'.format(y=conf['y_rot'], l=y_start, h=y_stop)
652623
conf['SiPM_sel']['cut'] += ' && {x} < {h} && {x} > {l})'.format(x=conf['x_rot'], l=x_start, h=x_stop)
653-
canvas['amp_w_pos'][k].Update()
654-
canvas['amp_w_pos'][k].SaveAs(out_dir + '/PositionXY_amp_weight_ch{:02d}.png'.format(k))
624+
625+
if 'space_sel' in conf.keys():
626+
x_start, x_stop, y_start, y_stop = conf['space_sel']['limits']
627+
y_sec, y_step = np.linspace(y_start, y_stop, 6, retstep=True)
628+
629+
leg = rt.TLegend(0.7,0.7,0.98,0.93)
630+
can.obj = [leg]
631+
632+
for i_s, (yd, yu) in enumerate(zip(y_sec[:-1], y_sec[1:])):
633+
pad = can.cd(1)
634+
line.SetLineColor(16)
635+
line.SetLineStyle(2)
636+
line.DrawLine(-width+dx, yu, width+dx, yu)
637+
if i_s == 0:
638+
line.DrawLine(-width+dx, yd, width+dx, yd)
639+
640+
h_aux = rt.TProfile('h_aux_strip_'+str(i_s), 'Average amplitude in selected events', N_bins[0], -width+dx, width+dx)
641+
h_aux.SetXTitle('x [mm]')
642+
h_aux.SetYTitle('Average amplitude [mV]')
643+
var_x = var[var.find(':')+1:]
644+
var_y = var[:var.find(':')]
645+
tmp_sel = selection + ' && {y}<{u} && {y}>{d}'.format(y=var_y, u=yu, d=yd)
646+
tmp_sel += ' && {x}<{u} && {x}>{d}'.format(x=var_x, u=x_stop, d=x_start)
647+
chain.Project('h_aux_strip_'+str(i_s), 'amp[' + str(k) + ']' + ':' + var_x, tmp_sel, 'prof')
648+
h_aux.SetStats(0)
649+
h_aux.GetYaxis().SetRangeUser(h_w.GetMinimum(10), 1.15*h_aux.GetMaximum())
650+
pad = can.cd(2)
651+
opt = 'PLC PMC'
652+
if i_s != 0:
653+
opt += ' SAME'
654+
h_aux.Draw(opt)
655+
leg.AddEntry(h_aux, ' y #in [{:.1f}, {:.1f}]'.format(yd, yu), 'l')
656+
can.obj.append(h_aux)
657+
658+
leg.Draw()
659+
660+
661+
can.Update()
662+
can.SaveAs(out_dir + '/PositionXY_amp_weight_ch{:02d}.png'.format(k))
655663

656664
''' -------------Avg integral-----------'''
657665
name = 'h_int_weight_pos_'+str(k)
@@ -765,9 +773,12 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
765773
s2 = '{y} < {y2}-{d} && {y} > {y1}+{d} && {x} < {x2}-{d} && {x} > {x1}+{d}'
766774
s = '!({} && !({}))'.format(s1, s2)
767775
# s = '!({})'.format(s1)
768-
selection += ' && ' + s.format(x=conf['x_rot'], y = conf['y_rot'], d=0.25, x1=x1, x2=x2, y1=y1, y2=y2)
776+
# selection += ' && ' + s.format(x=conf['x_rot'], y = conf['y_rot'], d=0.25, x1=x1, x2=x2, y1=y1, y2=y2)
769777

770778
# selection += ' && ' + configurations.channel[kk]['SiPM_sel']['cut']
779+
if args.verbose:
780+
print 'x_var:', conf['x_rot']
781+
print 'y_var:', conf['y_rot']
771782

772783
if chain.GetEntries(selection) < 5:
773784
print 'Not enought stat ({})'.format(chain.GetEntries(selection))
@@ -806,6 +817,10 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
806817
canvas['t_res_raw'][k].Update()
807818
canvas['t_res_raw'][k].SaveAs(out_dir + '/TimeResolution_raw_ch{:02d}.png'.format(k))
808819

820+
selection += '&& ({v} < {h} && {v} > {l})'.format(v=var_dT, h=median+2*width, l=median-2*width)
821+
if args.verbose:
822+
print 'selection:', selection
823+
809824
'''=========================== Time resolution 2D ==========================='''
810825
if 'TimeResRaw2D' in configurations.plots:
811826
if not 'shape' in conf.keys():

0 commit comments

Comments
 (0)