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

Last change on this file since 13385 was 13368, checked in by neise, 13 years ago
initial commit of my sandbox
File size: 10.7 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 pyfact
11from myhisto import *
12from hist import *
13import numpy as np
14import numpy.random as rnd
15from scipy import interpolate as ip
16from ROOT import *
17from time import time
18from optparse import OptionParser
19
20
21from pyfact import RawData
22from drs_spikes import DRSSpikes
23from extractors import ZeroXing
24class CalculateDrsTimeCalibrationConstants( object ):
25 """ basic analysis class for the calculation of DRS time aperture jitter calibration constants
26 """
27
28 def __init__(self, dpath, cpath):
29 """
30 *dpath*: data file path, file should contain a periodic signal in eahc 9th channel
31 *cpath*: std DRS amplitude calibration path
32 """
33 #dfname = '/data00/fact-construction/raw/2011/11/24/20111124_113.fits.gz'
34 #calfname = '/data00/fact-construction/raw/2011/11/24/20111124_111.drs.fits.gz'
35 self.data_path = dpath
36 self.calib_path = cpath
37 self.run = RawData(dpath, cpath, return_dict = True)
38 self.despike = DRSSpikes()
39 self.zX = ZeroXing()
40 self.fsampling = 2. # sampling frequency
41 self.freq = 250. # testfrequency
42 self.P_nom = 1000./freq # nominal Period due to testfrequency
43 self.DAMPING = 0.02 # Damping applied to correction
44 self.NChip = 1 # Number of Chips
45 self.NEvents = 300 # Number of Events to Process
46
47 self.t_nom = np.ones([run.npix, run.nroi])/fsampling
48
49 def __call__(self):
50 """ the class-call method performs the actual calculation based on the input data
51 returns: tuple of length 160 containing tuples of length 1024 containing the length of each cell in ns
52 the class contains self.result after the call.
53 """
54 # loop over events
55 for event in self.run:
56 data = despike( event['acal_data'] )
57 start_cells = event['start_cells']
58 MeanXing, TimeXing, NumXing = Calculate_Crossing(data, event['start_cells'])
59
60
61 def Calculate_Crossing(self, data, start_cells):
62 TimeXing = np.zeros(self.run.nroi)
63 MeanXing = np.zeros(self.run.nroi)
64 NumXing = 0
65 CellTime = CellTimef(np.roll(t_nom, -Start))
66
67 for i in range(rd.NROI-1):
68 if (Data[i] > Mean) & (Data[i+1] < Mean):
69 MeanXing[NumXing] = i
70
71 FirstCell = CellTime[i]
72 SecondCell = CellTime[i+1]
73
74 TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/Data[i])
75 NumXing += 1
76
77 # calculate the frequency of the testwave
78 freqp = fsampling*NumXing*1000./1024
79 if(np.abs(freqp-freq) > 20):
80 print "This is not a timecalibration run!!"
81 raise SystemExit
82
83 return MeanXing, TimeXing, NumXing
84
85
86 for Event in range( NEvents ):
87 print "Event ", Event
88 Chip = 0
89
90 # loop over every 9th channel in Chip
91 for Pix in range(NPIX)[8::9]:
92 print "Pix ", Pix
93
94 # calculate t_nom for each pixel
95 MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip)
96 h.list[Chip].dict["numxing"].Fill(NumXing)
97 #h.list[Chip].dict["start"].Fill(rd.startCells[Pix])
98 t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip)
99
100 Chip += 1
101
102 rd.next()
103
104 t_nom = t_nom[8::9]
105
106 if(save == "yes"):
107 np.save(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
108 np.savetxt(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
109
110 for Chip in range(NChip):
111 for i in range(rd.NROI):
112 h.list[Chip].dict["int_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][:i])-i*1/fsampling)
113 h.list[Chip].dict["diff_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][i] - 1/fsampling))
114
115 pyfact.SaveHistograms(h.list, str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip)+".root", "RECREATE")
116
117
118 t_f = time()
119
120
121
122 print "time ellapsed = ", (t_f - t_s)/60., " min."
123 print "Speichern = ", save
124 print "NChip = ", NChip
125 print "NEvents = ", NEventsp
126 print "Histo saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".root"
127 if(save == "yes"):
128 print "CellTime saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".npy"
129
130
131
132parser = OptionParser()
133parser.add_option('-e', '--extended', action = 'store_true', dest='extend', default = False )
134(options, args) = parser.parse_args()
135
136save = "yes"
137
138if(options.extend):
139
140 NChip = int(raw_input("Wieviele Chips? ([1, ... , "+str(rd.NPIX/9)+"]): "))
141 NEventsp = int(raw_input("Wieviele Events? ([1, ... , "+str(rd.NEvents)+"]) "))
142 save = raw_input("Soll CellTime gespeichert werden? (yes/no) ")
143
144
145h = hist_list(tcalHistograms, NChip, "Chip")
146
147
148# Functions needed for calibration
149
150# find negative-going crossings of mean value
151def Crossing(Data, Mean, t_nom, Start, Chip):
152 TimeXing = np.zeros(rd.NROI)
153 MeanXing = np.zeros(rd.NROI)
154 NumXing = 0
155 CellTime = CellTimef(np.roll(t_nom, -Start))
156
157 for i in range(rd.NROI-1):
158 if ((Data[i] > Mean) & (Data[i+1] < Mean)):
159 MeanXing[NumXing] = i
160
161 FirstCell = CellTime[i]
162 SecondCell = CellTime[i+1]
163
164 TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
165 NumXing += 1
166
167 # calculate the frequency of the testwave
168 freqp = fsampling*NumXing*1000./1024
169 if(np.abs(freqp-freq) > 20):
170 print "This is not a timecalibration run!!"
171 raise SystemExit
172
173 return MeanXing, TimeXing, NumXing
174
175def CrossingOliver(Data, Mean, t_nom, Start, Chip):
176 TimeXing = np.zeros(lpipe)
177 MeanXing = np.zeros(lpipe)
178 NumXing = 0
179 CellTime = CellTimef(t_nom)
180
181
182 for i in range(lpipe-1):
183 if ((Data[i] > Mean) & (Data[i+1] < Mean)):
184 MeanXing[NumXing] = np.mod(i+Start, lpipe)
185
186 FirstCell = CellTime[np.mod(i+Start,lpipe)]
187 SecondCell = CellTime[np.mod(i+1+Start, lpipe)]
188
189 if(SecondCell < FirstCell): SecondCell =+ lpipe/fsampl
190
191 TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
192 NumXing += 1
193
194
195 return MeanXing, TimeXing, NumXing
196
197# Determine time between crossings and apply correction
198def Corr(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
199
200 # loop over mean crossings
201 for i in range(int(NumXing-1)):
202 I1 = MeanXing[i]
203 I2 = MeanXing[i+1]
204 Period = TimeXing[i+1] - TimeXing[i]
205
206
207 if(np.abs(1-Period/P_nom) > 0.2):
208 #print "Skipping zero crossing due to too large deviation of period."
209 h.list[Chip].dict["loleftout"].Fill(MeanXing[i])
210 #h.list[Chip].dict["leftout"].Fill(np.mod((MeanXing[i]+Start), rd.NROI))
211 continue
212
213 #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i])
214 h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[i]+Start), rd.NROI))
215
216 # calculate the frequency of the testwave
217 freqp = 1000./Period
218 h.list[Chip].dict["freq"].Fill(freqp)
219
220 C = (P_nom-Period)*DAMPING
221 numnom = Period*fsampling
222 Correction = np.zeros(rd.NROI)
223 Correction[I1+1:I2+1] += C/numnom
224 Correction[:I1+1] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
225 Correction[I2+1:] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
226
227 Correction = np.roll(Correction, Start)
228 t_nom += Correction
229 #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[-1])
230 h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[-1]+Start), rd.NROI))
231
232 return t_nom
233
234def CorrOliver(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
235
236 CellTime = CellTimef(t_nom)
237 for i in range(int(NumXing-1)):
238 I1 = MeanXing[i]
239 I2 = MeanXing[i+1]
240 Period = TimeXing[i+1] - TimeXing[i]
241 if(Period<0): Period += rd.NROI/fsampl
242
243 if(np.abs(1-Period/P_nom) > 0.2):
244 #print "Skipping zero crossing due to too large deviation of period."
245 continue
246
247 Correction = (P_nom-Period)*DAMPING
248
249 # Positive correction from MeanXing[i] to MeanXing[i+1]
250 # and negative to all others
251 Time = 0
252 for j in range(rd.NROI):
253 diff = CellTime[j+1] - CellTime[j]
254 if ((I2>I1 and j>I1 and j<=I2) or (I2<I1 and (j<=I2 or j>I1))):
255 diff += Correction/np.mod(I2-I1+rd.NROI, rd.NROI)
256 else:
257 diff -= Correction/(rd.NROI-np.mod(I2-I1+rd.NROI, rd.NROI))
258 CellTime[j] = Time
259 Time += diff
260 CellTime[rd.NROI] = Time
261
262 print CellTime[-1]
263 t_nom = t_nomf(CellTime)
264
265 return t_nom
266
267
268
269# Summing up the nominal times for each cell to the total cell time.
270def CellTimef(t_nom):
271 CellTime = np.zeros(rd.NROI+1)
272 for i in range(rd.NROI):
273 CellTime[i+1] = CellTime[i] + t_nom[i]
274
275 return CellTime
276
277# getting t_nom from CellTime
278def t_nomf(CellTime):
279 t_nom = np.zeros(rd.NROI)
280 for i in range(rd.NROI):
281 t_nom[i] = CellTime[i+1]-CellTime[i]
282
283 return t_nom
284
285
286rd.next()
287NPIX = NChip*9
288NEvents = NEventsp
289
290# create array to put the calibrated times for each cell in
291t_nom = np.ones([NPIX, rd.NROI])/fsampling
292# loop over events
293for Event in range( NEvents ):
294 print "Event ", Event
295 Chip = 0
296
297 # loop over every 9th channel in Chip
298 for Pix in range(NPIX)[8::9]:
299 print "Pix ", Pix
300
301 # calculate t_nom for each pixel
302 MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip)
303 h.list[Chip].dict["numxing"].Fill(NumXing)
304 #h.list[Chip].dict["start"].Fill(rd.startCells[Pix])
305 t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip)
306
307 Chip += 1
308
309 rd.next()
310
311t_nom = t_nom[8::9]
312
313if(save == "yes"):
314 np.save(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
315 np.savetxt(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
316
317for Chip in range(NChip):
318 for i in range(rd.NROI):
319 h.list[Chip].dict["int_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][:i])-i*1/fsampling)
320 h.list[Chip].dict["diff_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][i] - 1/fsampling))
321
322pyfact.SaveHistograms(h.list, str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip)+".root", "RECREATE")
323
324
325t_f = time()
326
327
328
329print "time ellapsed = ", (t_f - t_s)/60., " min."
330print "Speichern = ", save
331print "NChip = ", NChip
332print "NEvents = ", NEventsp
333print "Histo saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".root"
334if(save == "yes"):
335 print "CellTime saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".npy"
336
Note: See TracBrowser for help on using the repository browser.