source: fact/tools/pyscripts/sandbox/dneise/timecal.py

Last change on this file was 13634, checked in by neise, 12 years ago
added comments
File size: 8.9 KB
Line 
1#!/usr/bin/python -tti
2#
3# calculation of DRS Time Calibration constants
4# initial implementation by Remo Dietlicher and Werner Lustermann (ETH Zurich)
5# based on C++ classes by Oliver Grimm (ETH Zurich)
6#
7# this python implementation was coded by
8# D. Neise (TU Dortmund) April 2012
9#
10import numpy as np
11from pyfact import RawData
12from drs_spikes import DRSSpikes
13from extractor import ZeroXing
14class CalculateDrsTimeCalibrationConstants( object ):
15 """ basic analysis class for the calculation of DRS time aperture jitter calibration constants
16 """
17
18 def __init__(self, dpath, cpath):
19 """
20 *dpath*: data file path, file should contain a periodic signal in eahc 9th channel
21 *cpath*: std DRS amplitude calibration path
22 """
23 #dfname = '/data00/fact-construction/raw/2011/11/24/20111124_113.fits.gz'
24 #calfname = '/data00/fact-construction/raw/2011/11/24/20111124_111.drs.fits.gz'
25 self.data_path = dpath
26 self.calib_path = cpath
27 self.run = RawData(dpath, cpath, return_dict = True)
28
29 if self.run.nroi != 1024:
30 raise TypeError('Time Calibration Data must have ROI=1024, ROI of '+dpath+' is: '+str(self.run.nroi))
31
32 self.despike = DRSSpikes()
33 self.zero_crossing_finder = ZeroXing(slope = -1) # -1 means falling edge
34 self.fsampling = 2. # sampling frequency
35 self.freq = 250. # testfrequency
36 self.P_nom = 1000./self.freq # nominal Period due to testfrequency
37 self.DAMPING = 0.02 # Damping applied to correction
38
39 # this ndarray will contain the best estimate of the true width of each slice after a __call__
40 # but it is used to determine the current best estimate inside the __call__ as well,
41 # since the algorithm is iteratively implemented
42 #self.time_calib = np.linspace( 0. , 512, num=1024, endpoint=False).repeat(run.npix/9).reshape( (1024,run.npix/9) ).transpose()
43 self.time_calib = np.ones( (self.run.npix/9, 1024) )/self.fsampling
44 def __call__(self):
45 """ the class-call method performs the actual calculation based on the input data
46 returns: tuple of length 160 containing tuples of length 1024 containing the length of each cell in ns
47 the class contains self.result after the call.
48 """
49 # loop over events
50 counter = 0
51 for event in self.run:
52 print 'counter', counter
53 counter +=1
54 # in the next() method, which is called during looping,
55 # the data is automatically amplitude calibrated.
56 # we just need every 9th pixel, so this is a place
57 # for optimization
58 data = event['acal_data']
59
60 # slice out every 9th channel
61 niners_ids = range(8,self.run.npix,9) # [8, 17, 26, ..., 1430, 1439]
62 data = data[niners_ids,:]
63 start_cell = event['start_cells'][niners_ids,:]
64
65 data = self.despike( data )
66
67 # shift data down by the mean # possible in one line
68 # calculate mean
69 mean = data.mean(1)
70 # make mean in the same shape as data, for elementwise subtraction
71 mean = mean.repeat(data.shape[1]).reshape(data.shape)
72 # subtraction
73 data = data - mean
74
75 # find zero crossings with falling edge
76 # the output is a list, of lists, of slices.
77 # each slice, shows a zero crossing. in fact, is is the slice where
78 # the data is still positive ... in the next slice it is negative, (or just zero).
79 # This was called 'MeanXing' before
80 all_crossings = self.zero_crossing_finder(data)
81
82 # now we want to calculate the exact time of the crossings
83 # but in units of nanoseconds, not in units of slices.
84 # In order to transform slices into nanoseconds, we already use self.time_calib
85 # rolled appropriately, and summed up in order to have the integral timing calibration
86 #
87 # the list, 'time_of_all_crossings' will contain sublists, for each chip
88 # which in turn will contain the times, in ns, for each zero-crossing
89 time_of_all_crossings = []
90 # the list, 'float_slice_of_all_crossings'
91 # is similar to the list above, but i just contains the position of the crossings
92 # in the pipeline, but more precisely than 'all_crossings', the position is linear interpolated
93 float_slice_of_all_crossings = []
94
95 for chip_id,crossings in enumerate(all_crossings):
96 time_of_all_crossings.append([])
97 float_slice_of_all_crossings.append([])
98 time = np.roll(self.time_calib[chip_id], -start_cell[chip_id]).cumsum()
99 for crossing in crossings:
100 # performing linear interpolation of time of zero crossing
101 # what fraction of a slide is the interpolated crossing away from the integral slice id
102 fraction_of_slice = data[chip_id,crossing] / (data[chip_id,crossing]-data[chip_id,crossing+1])
103 float_slice_of_all_crossings[-1].append(fraction_of_slice + crossing)
104 slope = time[crossing+1] - time[crossing]
105 time_of_crossing = time[crossing] + slope * fraction_of_slice
106 time_of_all_crossings[-1].append(time_of_crossing)
107
108 # Now one can use those two lists, which were just created
109 # to calculate the number(float) of slices between two consecutive crossings
110 # as well as the *measured* time between two crossings.
111 # on the other hand, the period of the test signal is known to be self.P_nom
112 # so in case the measured time is longer or shorter than the known Period
113 # the associated time of all involved slices
114 # is assumed to be a bit longer or short than 1./self.fsampling
115
116 # first we make the sublists of the lists to be numpy.arrays
117 for i in range(len(float_slice_of_all_crossings)): #both list have equal length
118 float_slice_of_all_crossings[i] = np.array(float_slice_of_all_crossings[i])
119 time_of_all_crossings[i] = np.array(time_of_all_crossings[i])
120
121 # now we can easily calculate the measured periods using np.diff()
122 all_measured_periods = []
123 for chip_times in time_of_all_crossings:
124 all_measured_periods.append(np.diff(chip_times))
125
126 for chip,chip_periods in enumerate(all_measured_periods):
127
128 for i,period in enumerate(chip_periods):
129 #shortcuts to the involved slices
130 start = all_crossings[chip][i]
131 end = all_crossings[chip][i+1]
132 interval = end-start
133 rest = 1024-interval
134
135 delta_period = (self.P_nom - period) * self.DAMPING
136 # nominal_number_of_slices
137 nom_num_slices = period * self.fsampling
138
139 # now the delta_period is distributed to all slices
140 # I do not exactly understand the +1 here ... error prone!!!
141 # the following explanation assumes delta_period is positive
142 # the slices between start and end get a positive correction value
143 # This correction value is given by delta_period / nom_num_slices
144 pos_fraction = delta_period / nom_num_slices
145 # the total positive correction sums up to
146 total_positive_correction = interval * pos_fraction
147 # the total negative correction must sum up to the sampe value
148 total_negative_correction = rest * (interval * (pos_fraction/rest))
149 # so we can call the product following 'rest' the 'neg_fraction
150 neg_fraction = -1 * interval/rest * pos_fraction
151 #print '1. should be zero', total_positive_correction - total_negative_correction
152 #assert total_positive_correction - total_negative_correction == 0
153
154 correction = np.zeros(1024)
155 correction[start+1:end+1] = pos_fraction
156 correction[:start+1] = neg_fraction
157 correction[end+1:] = neg_fraction
158 #print '2. should be zero', correction.sum()
159 #assert correction.sum() == 0
160
161 # now we have to add these correction values to self.time_calib
162 self.time_calib[chip,:] += np.roll(correction, start_cell[chip])
Note: See TracBrowser for help on using the repository browser.