Skip to content

Commit 3310e52

Browse files
ocerriocerri
committed
Tracks Consistency and SiPMs positioning added
1 parent b27d792 commit 3310e52

File tree

5 files changed

+249
-310
lines changed

5 files changed

+249
-310
lines changed

‎DQM_SiPM.py‎

Lines changed: 225 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,13 @@
11
import numpy as np
22
import os, re, shutil
33
import argparse
4+
import itertools
45
import ROOT as rt
56
from root_numpy import tree2array
67
from lib.histo_utilities import create_TH1D, create_TH2D
7-
from lib.cebefo_style import cebefo_style
8+
from lib.cebefo_style import cebefo_style, Set_2D_colz_graphics
9+
10+
donotdelete = []
811

912

1013

@@ -55,6 +58,18 @@ def __init__(self, infile):
5558
print 'Cleaning tracks with:', cuts[:-4]
5659
else:
5760
self.TracksCleaning['cuts'] = ''
61+
elif '-->TracksConsistency' in l[0]:
62+
self.TracksConsistency = {}
63+
self.TracksConsistency['Bad'] = 0
64+
self.TracksConsistency['Checked'] = 0
65+
for aux in l[1:]:
66+
if aux[:5] == 'RefCh':
67+
self.TracksConsistency['RefCh'] = int(aux[5:])
68+
print 'Tracks consistency: ch', int(aux[5:]), 'will be checked'
69+
70+
if aux == 'List':
71+
self.TracksConsistency['List'] = 1
72+
5873
elif len(self.labels) == 0:
5974
self.labels = l[1:]
6075
else:
@@ -107,6 +122,36 @@ def define_range_around_peak(h, perc = [0.4, 0.4], Range=[0.0, 99999.0]):
107122
x_up = h.GetBinLowEdge(n_up) + h.GetBinWidth(n_up)
108123
return x_low, x_up, n_pk
109124

125+
def rootTH2_to_np(h, cut = None, Norm = False):
126+
nx = h.GetNbinsX()
127+
ny = h.GetNbinsY()
128+
129+
arr = np.zeros((ny, nx))
130+
pos = np.zeros((ny, nx, 2))
131+
132+
for ix in range(nx):
133+
for iy in range(ny):
134+
x = h.GetXaxis().GetBinCenter( ix+1 );
135+
y = h.GetYaxis().GetBinCenter( iy+1 );
136+
z = h.GetBinContent(h.GetBin(ix+1, iy+1))
137+
138+
if cut == None:
139+
arr[iy, ix] = z
140+
else:
141+
arr[iy, ix] = z if z > cut else 0
142+
pos[iy, ix] = [x,y]
143+
return arr, pos
144+
145+
def circle_filter(h_arr, cx, cy, rx, ry):
146+
p = 0
147+
for i,j in itertools.product(np.arange(h_arr.shape[0]), np.arange(h_arr.shape[1])):
148+
if (float(cx-j)/rx)**2 + (float(cy-i)/ry)**2 < 1:
149+
p += h_arr[i-1, j-1]
150+
151+
return p/(np.pi*rx*ry)
152+
153+
154+
110155
if __name__ == '__main__':
111156
cebefo_style()
112157

@@ -152,7 +197,7 @@ def define_range_around_peak(h, perc = [0.4, 0.4], Range=[0.0, 99999.0]):
152197
if os.path.isdir(save_loc):
153198
if save_loc[-1] != '/':
154199
save_loc += '/'
155-
out_dir = save_loc + 'Run' + str(flag) + '_plots'
200+
out_dir = save_loc + 'Run' + str(flag) + '_DQM'
156201
if os.path.exists(out_dir):
157202
shutil.rmtree(out_dir)
158203
os.mkdir(out_dir)
@@ -187,6 +232,7 @@ def define_range_around_peak(h, perc = [0.4, 0.4], Range=[0.0, 99999.0]):
187232
canvas['risetime'] = {}
188233
canvas['wave'] = {}
189234
canvas['pos'] = {}
235+
canvas['pos_sel'] = {}
190236
canvas['w_pos'] = {}
191237
canvas['t_res_raw'] = {}
192238
canvas['space_corr'] = {}
@@ -391,40 +437,180 @@ def define_range_around_peak(h, perc = [0.4, 0.4], Range=[0.0, 99999.0]):
391437
dy = configurations.xy_center[1]
392438
dx = configurations.xy_center[0]
393439
width = configurations.xy_center[2]
440+
N_bins = [100, 100]
394441
if 'PosRaw' in configurations.plots:
395442
name = 'h_pos_'+str(k)
396443
title = 'Track position at z_DUT[' + str(conf['idx_dut']) + '], for channel {} selected event'.format(k)
397444

398-
h = rt.TH2D(name, title, 100, -width+dx, width+dx, 100, -width+dy, width+dy)
445+
h = rt.TH2D(name, title, N_bins[0], -width+dx, width+dx, N_bins[1], -width+dy, width+dy)
399446
h.SetXTitle('x [mm]')
400447
h.SetYTitle('y [mm]')
401448

402449
chain.Project(name, var, selection)
403450

404451
canvas['pos'][k] = rt.TCanvas('c_pos_'+str(k), 'c_pos_'+str(k), 800, 600)
405452
h.DrawCopy('colz')
406-
407453
canvas['pos'][k].Update()
408454
canvas['pos'][k].SaveAs(out_dir + '/PositionXY_raw_ch{:02d}.png'.format(k))
409455

410-
if 'PosWeight' in configurations.plots:
411-
name = 'h_weight_pos_'+str(k)
412-
title = 'Track position at z_DUT[' + str(conf['idx_dut']) + '] weighted with channel {} selected event'.format(k)
413-
h_w = rt.TH2D(name, title, 100, -width+dx, width+dx, 100, -width+dy, width+dy)
414-
h_w.SetXTitle('x [mm]')
415-
h_w.SetYTitle('y [mm]')
416-
h_w.SetZTitle('Average Integral [pC]')
456+
if ('PosSel' in configurations.plots) or ('PosSel+' in configurations.plots):
457+
canvas['pos_sel'][k] = rt.TCanvas('c_pos_sel_'+str(k), 'c_pos_sel_'+str(k), 800, 600)
458+
459+
name = 'h_pos_sel'
460+
title = 'Events average selection efficiency'
461+
h = rt.TH3D(name, title, N_bins[0], -width+dx, width+dx, N_bins[1], -width+dy, width+dy, 2, -0.5, 1.5)
417462

418-
weights = '('+ selection +') * amp[' + str(k) + ']'
419-
chain.Project(name, var, weights)
463+
chain.Project(name, selection + ':' + var)
464+
h = h.Project3DProfile('yx')
465+
h.SetYTitle('y [mm]')
466+
h.SetXTitle('x [mm]')
467+
h.DrawCopy('colz')
420468

421-
if 'PosRaw' in configurations.plots:
422-
h_w.Divide(h)
469+
# if hasattr(configurations, 'TracksConsistency'):
470+
# if configurations.TracksConsistency['RefCh'] == k:
471+
# x_mean = h.GetMean(1)
472+
# bx = h.GetXaxis().FindBin(x_mean)
473+
# x_sx = h.GetRMS(1)
474+
# bx_sx = h.GetXaxis().GetBinCenter(1) - h.GetXaxis().GetBinWidth(1)*0.51 + x_sx
475+
# bx_sx = h.GetXaxis().FindBin(bx_sx)
476+
#
477+
# y_mean = h.GetMean(2)
478+
# by = h.GetYaxis().FindBin(y_mean)
479+
# y_sx = h.GetRMS(2)
480+
# by_sx = h.GetYaxis().GetBinCenter(1) - h.GetYaxis().GetBinWidth(1)*0.51 + y_sx
481+
# by_sx = h.GetYaxis().FindBin(by_sx)
482+
#
483+
# arr_h, pos_h = rootTH2_to_np(h, cut=0.2)
484+
#
485+
#
486+
# best_scale = 0.5
487+
# p_max = 0
488+
# for s in np.arange(0.5, 3, 0.05):
489+
# p = circle_filter(arr_h, bx, by, s*bx_sx, s*by_sx)
490+
# print p
491+
#
492+
# if p > p_max:
493+
# p_max = p
494+
# best_scale = s
495+
#
496+
# print best_scale
497+
#
498+
# ell = rt.TEllipse(x_mean, y_mean, best_scale*x_sx, best_scale*y_sx)
499+
# ell.SetLineColor(6)
500+
# ell.SetLineStyle(9)
501+
# ell.SetLineWidth(2)
502+
# ell.SetFillStyle(0)
503+
# ell.Draw()
504+
# donotdelete.append(ell)
505+
506+
507+
508+
if ('PosSel+' in configurations.plots) and ('shape' in conf.keys()):
509+
if not conf['shape'] == 'None':
510+
size = conf['shape'].split('x')
511+
512+
x = h.GetXaxis().GetBinCenter(1) - h.GetXaxis().GetBinWidth(1)*0.51 + float(size[0])
513+
nbx = h.GetXaxis().FindBin(x)
514+
515+
y = h.GetYaxis().GetBinCenter(1) - h.GetYaxis().GetBinWidth(1)*0.51 + float(size[1])
516+
nby = h.GetYaxis().FindBin(y)
517+
# print h.GetYaxis().GetBinCenter(1), y, nby
518+
519+
# print nby, nbx
520+
# x_mean = h.GetMean(1)
521+
# x_border = [x_mean - 0.5*float(size[0]), x_mean + 0.5*float(size[0])]
522+
# bx_border = [h.GetXaxis().FindBin(x_border[0]), h.GetXaxis().FindBin(x_border[1])]
523+
#
524+
# y_mean = h.GetMean(1)
525+
# y_border = [y_mean - 0.5*float(size[0]), y_mean + 0.5*float(size[0])]
526+
# by_border = [h.GetXaxis().FindBin(y_border[0]), h.GetXaxis().FindBin(y_border[1])]
527+
528+
arr_h, pos_h = rootTH2_to_np(h, cut=0.2)
529+
# print arr_h
530+
531+
p_max = 0
532+
idx_max = [0, 0]
533+
for iy in range(0, arr_h.shape[0]-nby):
534+
for ix in range(0, arr_h.shape[1]-nbx):
535+
p = np.sum(arr_h[iy:iy+nby, ix:ix+nbx])
536+
537+
if p > p_max:
538+
p_max = p
539+
idx_max = [iy, ix]
540+
541+
avg_prob = p_max/(nbx*nby)
542+
print avg_prob
543+
544+
if hasattr(configurations, 'TracksConsistency'):
545+
configurations.TracksConsistency['Checked'] += 1
546+
x_mean = h.GetMean(1)
547+
bx = h.GetXaxis().FindBin(x_mean)
548+
x_sx = h.GetRMS(1)
549+
bx_sx = h.GetXaxis().GetBinCenter(1) - h.GetXaxis().GetBinWidth(1)*0.51 + x_sx
550+
bx_sx = h.GetXaxis().FindBin(bx_sx)
551+
552+
y_mean = h.GetMean(2)
553+
by = h.GetYaxis().FindBin(y_mean)
554+
y_sx = h.GetRMS(2)
555+
by_sx = h.GetYaxis().GetBinCenter(1) - h.GetYaxis().GetBinWidth(1)*0.51 + y_sx
556+
by_sx = h.GetYaxis().FindBin(by_sx)
557+
558+
p_circle = -1
559+
if (by_sx != 0) and (bx_sx != 0):
560+
p_circle = circle_filter(arr_h, bx, by, 2*bx_sx, 2*by_sx)
561+
562+
print p_circle
563+
if p_circle > avg_prob*0.75 or p_circle == -1:
564+
print 'Possilble random scatter shape'
565+
configurations.TracksConsistency['Bad'] += 1
566+
567+
568+
569+
iy, ix = idx_max
570+
# arr_filter = np.zeros_like(arr_h)
571+
# arr_filter[iy:iy+nby, ix:ix+nbx] = 1.
572+
# print arr_filter
573+
574+
x_start = h.GetXaxis().GetBinCenter(ix+1) - 0.5*h.GetXaxis().GetBinWidth(ix+1)
575+
x_stop = h.GetXaxis().GetBinCenter(ix+nbx) + 0.5*h.GetXaxis().GetBinWidth(ix+nbx)
576+
y_start = h.GetYaxis().GetBinCenter(iy+1) - 0.5*h.GetYaxis().GetBinWidth(iy+1)
577+
y_stop = h.GetYaxis().GetBinCenter(iy+nby) + 0.5*h.GetYaxis().GetBinWidth(iy+nby)
578+
579+
line.SetLineStyle(9)
580+
line.SetLineColor(6)
581+
line.SetLineWidth(2)
582+
line.DrawLine(x_start, y_start, x_stop, y_start)
583+
line.DrawLine(x_start, y_stop, x_stop, y_stop)
584+
line.DrawLine(x_start, y_start, x_start, y_stop)
585+
line.DrawLine(x_stop, y_start, x_stop, y_stop)
586+
587+
Set_2D_colz_graphics()
588+
canvas['pos_sel'][k].Update()
589+
canvas['pos_sel'][k].SaveAs(out_dir + '/PositionXY_sel_ch{:02d}.png'.format(k))
590+
591+
592+
593+
594+
if 'PosWeight' in configurations.plots:
595+
name = 'h_weight_pos_'+str(k)
596+
title = 'Track position at z_DUT[' + str(conf['idx_dut']) + '] weighted with channel {} selected event'.format(k)
597+
h_w = rt.TH2D(name, title, N_bins[0], -width+dx, width+dx, N_bins[1], -width+dy, width+dy)
598+
h_w.SetXTitle('x [mm]')
599+
h_w.SetYTitle('y [mm]')
600+
h_w.SetZTitle('Average Integral [pC]')
423601

424-
canvas['w_pos'][k] = rt.TCanvas('c_w_pos_'+str(k), 'c_w_pos_'+str(k), 800, 600)
425-
h_w.DrawCopy('colz')
426-
canvas['w_pos'][k].Update()
427-
canvas['w_pos'][k].SaveAs(out_dir + '/PositionXY_amp_weight_ch{:02d}.png'.format(k))
602+
h = rt.TH2D('h_amp_aux'+str(k), title, N_bins[0], -width+dx, width+dx, N_bins[1], -width+dy, width+dy)
603+
chain.Project('h_amp_aux'+str(k), var, selection)
604+
605+
weights = '('+ selection +') * amp[' + str(k) + ']'
606+
chain.Project(name, var, weights)
607+
608+
h_w.Divide(h)
609+
610+
canvas['w_pos'][k] = rt.TCanvas('c_w_pos_'+str(k), 'c_w_pos_'+str(k), 800, 600)
611+
h_w.DrawCopy('colz')
612+
canvas['w_pos'][k].Update()
613+
canvas['w_pos'][k].SaveAs(out_dir + '/PositionXY_amp_weight_ch{:02d}.png'.format(k))
428614

429615
'''=========================== End Selections ==========================='''
430616
conf['sel'] = selection
@@ -638,3 +824,23 @@ def create_regression_input(x, y, amp):
638824

639825
canvas['dt_corr'][k].Update()
640826
canvas['dt_corr'][k].SaveAs(out_dir + '/TimeResolution_OneShot_ch{:02d}.png'.format(k))
827+
828+
'''End of channels loop'''
829+
if hasattr(configurations, 'TracksConsistency'):
830+
aux = float(configurations.TracksConsistency['Bad'])/configurations.TracksConsistency['Checked']
831+
if aux > 0.3:
832+
print '\n\n============ Run to be discarted!!!!! ===============\n\n'
833+
834+
if 'List' in configurations.TracksConsistency.keys():
835+
f_aux = flag[1:]
836+
file = args.save_loc
837+
if aux > 0.3:
838+
file += 'TracksConsistency_Bad.txt'
839+
else:
840+
file += 'TracksConsistency_Good.txt'
841+
if not os.path.exists(file):
842+
os.system('touch '+file)
843+
if not (int(f_aux) in np.loadtxt(file).astype(np.int)):
844+
cmd = 'echo ' + f_aux + ' >> ' + file
845+
print cmd
846+
os.system(cmd)

0 commit comments

Comments
 (0)