source: fact/tools/pyscripts/sandbox/dneise/cleaning/begin.py

Last change on this file was 13147, checked in by neise, 13 years ago
enabled tab checking and excahnged python_rename with python
  • Property svn:executable set to *
File size: 4.2 KB
Line 
1#!/usr/bin/python -itt
2#
3# Dominik Neise
4#
5# cleaning a small step towards the truth
6from pyfact import RawData
7import os.path
8import matplotlib.pyplot as plt
9import numpy as np
10from fir_filter import *
11from extractor import *
12from drs_spikes import *
13from plotters import *
14import time as t
15from coor import Coordinator
16confirm_next_step = False# this is for user interaction
17
18data_file_name = '/media/daten_platte/FACT/data/20120229_144.fits.gz'
19calib_file_name = '/media/daten_platte/FACT/data/20120229_132.drs.fits.gz'
20if not os.path.isfile(data_file_name):
21 print 'not able to find file:', data_file_name
22 sys.exit(-1)
23if not os.path.isfile(calib_file_name ):
24 print 'not able to find file:', calib_file_name
25 sys.exit(-1)
26
27run = RawData(data_file_name, calib_file_name)
28despike = DRSSpikes()
29smooth = SlidingAverage(8)
30extract = GlobalMaxFinder(40,200)
31
32#plotA = CamPlotter('amplitudes')
33#plotT = CamPlotter('times')
34#plotCA = CamPlotter('cleaned amplitudes')
35
36#plotArea = HistPlotter('area', 1440, (0,1440) )
37#plotSize = HistPlotter('size', 1000, (0,10000) )
38
39#funnyPlot = Plotter('area vs log(size')
40
41
42coreTHR = 45 # copied from F. Temme
43edgeTHR = 18 # copied from F. Temme
44
45# get dictionary of next neighbors
46nn = Coordinator().nn
47
48print 'Core threshold=', coreTHR
49print 'Edge threshold=', edgeTHR
50print '*'*70
51
52areas = []
53sizes = []
54
55for data,startcell,tt in run:
56 if tt==4:
57 data = despike(data)
58 data = smooth(data)
59 amplitude, time_of_max = extract(data)
60
61 #plotA.name='amplitudes EvtID:' + str(run.event_id.value) + ' TT:' + str(tt)
62 #plotA(amplitude)
63 #plotT(time_of_max)
64
65 # Here we start the cleaning ... just like that...
66 #print 'cleaning away all pixel <' , coreTHR
67 core_chid_candidates = np.where( amplitude > coreTHR)[0]
68 core_chids = [] # coor chids, which survive Gauks step 1
69 survivors = [] # coor & edge pixel
70 #print 'number of core candidates:', len(core_chid_candidates)
71
72 #print 'throwing away all pixel w/o any neighbor'
73 # get rid of single coor pics
74 for cand in core_chid_candidates:
75 neighbor_found = False
76 # loop over all neigbors of cand'idate
77 for n in nn[cand]:
78 if n in core_chid_candidates:
79 neighbor_found = True
80 if neighbor_found:
81 core_chids.append(cand)
82
83 #print 'after deletion of single coor pixels'
84
85
86 #add edge pixel to the edge of the coors
87 #print 'resurrecting edge pixels ... i.e. all pixel >', edgeTHR
88 survivors = core_chids[:]
89 for coor in core_chids:
90 for n in nn[coor]:
91 # if neighbor is a core pixel, then do nothing
92 if n in core_chids:
93 pass
94 elif amplitude[n] > edgeTHR:
95 survivors.append(n)
96 survivors = np.array(survivors, dtype=int)
97 print '#survivors', len(survivors), 'evtID', run.event_id.value
98
99 #plotCA(data=amplitude, mask=survivors)
100
101 size = 0
102 for pixel in survivors:
103 size += amplitude[pixel]
104
105 if len(survivors) > 0:
106 areas.append( len(survivors) )
107 sizes.append( size )
108
109
110 #plotArea(areas, 'areas of ' + str(run.event_id.value) + 'events')
111 #plotSize(sizes, 'sizes of ' + str(run.event_id.value) + 'events')
112
113 if confirm_next_step:
114 user_input = raw_input("'q'-quit, 'r'-run, anything else goes one step")
115 number=None
116 try:
117 number=int(user_input)
118 except:
119 number=None
120 if user_input.find('q') != -1:
121 sys.exit(0)
122 elif user_input.find('r') != -1:
123 confirm_next_step = False
124 elif number!=None:
125 run += number
126
127 total_event_number = run.event_id.value
128
129plt.ion()
130myfig = plt.figure()
131myn = myfig.number
132logsize = np.log10(np.array(sizes))
133areas = np.array(areas)
134
135plt.figure(myn)
136plt.title('area vs. log10(size) of '+ str(total_event_number) + 'events')
137plt.xlabel('log10(size/1mV)')
138plt.ylabel('area [#pixel]')
139plt.plot( logsize,areas, '.')
Note: See TracBrowser for help on using the repository browser.