source: fact/tools/pyscripts/sandbox/dneise/rdietlic/ftcalsim.py@ 16858

Last change on this file since 16858 was 13368, checked in by neise, 13 years ago
initial commit of my sandbox
File size: 11.0 KB
Line 
1# simulation of the real data for timecalibration
2#
3# Remo Dietlicher
4# ETH Zürich / Semesterarbeit
5#
6#
7
8import pyfact
9from myhisto import *
10from hist import *
11import numpy as np
12import numpy.random as rnd
13from ROOT import *
14from time import time
15from optparse import OptionParser
16
17NROI = 1024 # length of the DRS pipeline
18NEvents = 200 # Number of evaluated Event
19fsampling = 1024./512.1 # sampling frequency
20freq = 250. # testfrequency
21P_nom = 1000./freq # nominal Period due to testfrequency
22DAMPING = 0.02 # Damping applied to correction
23NChip = 1 # Number of Chips
24
25t_s = time()
26
27parser = OptionParser()
28parser.add_option('-e', '--extended', action = 'store_true', dest='extend', default = False )
29(options, args) = parser.parse_args()
30
31variante = "dat"
32alg = "Remo"
33save = "no"
34
35if(options.extend):
36
37 variante = raw_input("Welche Daten? (sim/dat): ")
38 alg = raw_input("Welcher Algorithmus? (Remo/Oliver): ")
39 NChip = int(raw_input("Wieviele Chips? ([1, ... ,160]): "))
40 NEvents = int(raw_input("Wieviele Events? ([1, ... , 1000]) "))
41 save = raw_input("Soll CellTime gespeichert werden? (yes/no) ")
42 tag = raw_input("Soll ein spezieller tag an das File angehängt werden? ")
43
44NPix = 9*NChip
45
46h = hist_list(tcalsimHistograms, NChip, "Chip")
47MeanXingHist = np.zeros(NROI)
48
49def Datainit(variante, NEvents):
50
51 if(variante == "sim"):
52 y_off = 3.4
53
54 # Start is the first bin being written on.
55 Startp = rnd.random(NEvents)*(NROI - 2) # project [0, 1] onto int([0, lpipe-1])
56 Start = range(NEvents)
57 for i, val in enumerate(Startp):
58 Start[i] = int(round(val)) # getting integer values needed in roll later
59
60 ID = np.zeros(NROI)
61 x = np.linspace(0., float(NROI), NROI+1)
62 ID = np.sin(x*2*np.pi/NROI+np.pi)*5
63
64 dt_true = np.zeros(NROI)
65 for i in range(NROI):
66 dt_true[i] = ID[i+1]-ID[i]
67
68 versch = np.sum(dt_true)/1024
69 dt_true = dt_true - versch
70 print np.sum(dt_true)
71
72 #create time-array, time in nanoseconds
73 t_nom = np.ones(NROI)/fsampling
74 t_true = t_nom + dt_true
75
76 # simulation of the real pipeline
77 rCellTime = np.zeros([NEvents, NROI])
78 for Event in range(NEvents):
79 for i in range(NROI-1):
80 rCellTime[Event][i+1] = rCellTime[Event][i]+t_true[i]
81 rCellTime[Event] = np.roll(rCellTime[Event],-Start[Event])
82 rCellTime[-Start[Event]:] += NROI/fsampling
83
84 print "t_true tot = ", np.sum(t_true)
85
86 # calculation of the nominal pipeline
87 def CellTimef(t_nom):
88 CellTime = np.zeros(NROI+1)
89 for i in range(NROI):
90 CellTime[i+1] = CellTime[i] + t_nom[i]
91
92 return CellTime
93
94 CellTime = CellTimef(t_nom)
95
96 # Simulate Data
97 Data = np.zeros([NEvents, NROI])
98 for Event in range(NEvents):
99 Data[Event] = np.sin(rCellTime[Event]*2.*np.pi/P_nom)+y_off
100
101 if(variante == "dat"):
102 t_true = 0
103 Data = np.load("d1000.npy")
104 Start = np.load("start_cells_1.npy")
105
106
107 return Data, t_true, Start
108
109Data, t_true, Start = Datainit(variante, NEvents)
110
111
112# Functions needed for calibration
113# calculation of the nominal pipeline
114def CellTimef(t_nom):
115 CellTime = np.zeros(NROI+1)
116 for i in range(NROI):
117 CellTime[i+1] = CellTime[i] + t_nom[i]
118
119 return CellTime
120
121# getting t_nom from CellTime
122def t_nomf(CellTime):
123 t_nom = np.zeros(NROI)
124 for i in range(NROI):
125 t_nom[i] = CellTime[i+1]-CellTime[i]
126
127 return t_nom
128
129
130
131# find negative-going crossings of mean value
132def Crossing(Data, Mean, t_nom, Start, Chip):
133 TimeXing = np.zeros(NROI)
134 MeanXing = np.zeros(NROI)
135 NumXing = 0
136 CellTime = CellTimef(np.roll(t_nom, -Start))
137
138 for i in range(NROI-1):
139 if ((Data[i] > Mean) & (Data[i+1] < Mean)):
140 MeanXing[NumXing] = i
141 h.list[Chip].dict["lomeanxing"].Fill(np.mod(i+Start, NROI))
142
143 FirstCell = CellTime[i]
144 SecondCell = CellTime[i+1]
145
146 TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
147 NumXing += 1
148
149 # calculate the frequency of the testwave
150 freqp = fsampling*NumXing*1000./1024
151 h.list[Chip].dict["meanfreq"].Fill(freqp)
152 if(np.abs(freqp-freq) > 20):
153 print "This is not a timing calibration run!!"
154 raise SystemExit
155
156 return MeanXing, TimeXing, NumXing
157
158
159
160# find positive-going crossings of mean value
161def CrossingPos(Data, Mean, t_nom, Start, Chip):
162 TimeXing = np.zeros(NROI)
163 MeanXing = np.zeros(NROI)
164 NumXing = 0
165 CellTime = CellTimef(np.roll(t_nom, -Start))
166
167 for i in range(NROI-1):
168 if ((Data[i] < Mean) & (Data[i+1] > Mean)):
169 MeanXing[NumXing] = i
170 h.list[Chip].dict["lomeanxing"].Fill(np.mod(i+Start, NROI))
171
172 FirstCell = CellTime[i]
173 SecondCell = CellTime[i+1]
174
175 TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
176 NumXing += 1
177
178 # calculate the frequency of the testwave
179 freqp = fsampling*NumXing*1000./1024
180 h.list[Chip].dict["meanfreq"].Fill(freqp)
181 if(np.abs(freqp-freq) > 20):
182 print "This is not a timing calibration run!!"
183 raise SystemExit
184
185 return MeanXing, TimeXing, NumXing
186
187
188
189def CrossingOliver(Data, Mean, t_nom, Start, Chip):
190
191 TimeXing = np.zeros(NROI)
192 MeanXing = np.zeros(NROI)
193 NumXing = 0
194 CellTime = CellTimef(t_nom)
195
196
197 for i in range(NROI-1):
198 if ((Data[i] > Mean) & (Data[i+1] < Mean)):
199 MeanXing[NumXing] = np.mod(i+Start, NROI)
200
201 FirstCell = CellTime[np.mod(i+Start,NROI)]
202 SecondCell = CellTime[np.mod(i+1+Start, NROI)]
203
204 if(SecondCell < FirstCell): SecondCell =+ NROI/fsampl
205
206 TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
207 NumXing += 1
208
209 # calculate the frequency of the testwave
210 freqp = fsampling*NumXing*1000./1024
211 h.list[Chip].dict["meanfreq"].Fill(freqp)
212 if(np.abs(freqp-freq) > 20):
213 print "This is not a timing calibration run!!"
214 raise SystemExit
215
216 return MeanXing, TimeXing, NumXing
217
218# Determine time between crossings and apply correction
219def Corr(MeanXing, NumXing, TimeXing, t_nom, Start, Chip, Event):
220
221 Ctot = np.zeros(NumXing)
222
223 # loop over mean crossings
224 for i in range(int(NumXing-1)):
225 I1 = MeanXing[i]
226 I2 = MeanXing[i+1]
227 Period = TimeXing[i+1] - TimeXing[i]
228 if(Event > 300):
229 h.list[Chip].dict["period"].Fill(Period)
230
231
232 if(np.abs(1-Period/P_nom) > 0.2):
233 #print "Skipping zero crossing due to too large deviation of period."
234 h.list[Chip].dict["loleftout"].Fill(MeanXing[i])
235 h.list[Chip].dict["leftout"].Fill(np.mod(MeanXing[i]+Start, NROI))
236 continue
237
238 h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i])
239 h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[i]+Start, NROI))
240 MeanXingHist[np.mod(MeanXing[i]+Start, NROI)] += 1
241
242 # calculate the frequency of the testwave
243 freqp = 1000./Period
244 h.list[Chip].dict["freq"].Fill(freqp)
245
246 C = (P_nom-Period)*DAMPING
247 Ctot[i] = P_nom-Period
248
249 numnom = Period*fsampling
250 Correction = np.zeros(NROI)
251 Correction[I1+1:I2+1] += C/numnom
252 Correction[:I1+1] -= (I2-I1)*C/(numnom*(NROI-(I2-I1)))
253 Correction[I2+1:] -= (I2-I1)*C/(numnom*(NROI-(I2-I1)))
254 """
255 Correction = np.zeros(NROI)
256 Correction[I1+1:I2+1] += C/(I2-I1)
257 Correction[:I1+1] -= C/(NROI-(I2-I1))
258 Correction[I2+1:] -= C/(NROI-(I2-I1))
259
260 numnom = Period*fsampling
261 Correction = np.zeros(NROI)
262 Correction[I1+1:I2+1] += C/numnom
263 """
264
265 Correction = np.roll(Correction, Start)
266 t_nom += Correction
267 #t_nom = t_nom/np.sum(t_nom)*1024/fsampling
268
269 if(Event > 200):
270 h.list[Chip].dict["perioddevproj"].Fill(np.average(Ctot))
271 h.list[Chip].dict["perioddev"].SetBinContent(Event+1, np.abs(np.average(Ctot)))
272 h.list[Chip].dict["lomeanxing"].Fill(MeanXing[NumXing])
273 h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[NumXing]+Start, NROI))
274 MeanXingHist[np.mod(MeanXing[NumXing]+Start, NROI)] += 1
275
276
277 return t_nom
278
279def CorrOliver(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
280
281 CellTime = CellTimef(t_nom)
282 for i in range(int(NumXing-1)):
283 I1 = MeanXing[i]
284 I2 = MeanXing[i+1]
285 Period = TimeXing[i+1] - TimeXing[i]
286 if(Period<0): Period += NROI/fsampling
287 h.list[Chip].dict["period"].Fill(Period)
288
289 if(np.abs(1-Period/P_nom) > 0.2):
290 print "Skipping zero crossing due to too large deviation of period."
291 h.list[Chip].dict["loleftout"].Fill(MeanXing[i])
292 h.list[Chip].dict["leftout"].Fill(np.mod(MeanXing[i]+Start, NROI))
293 continue
294
295 # calculate the frequency of the testwave
296 freqp = 1000./Period
297 h.list[Chip].dict["freq"].Fill(freqp)
298
299
300 h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i])
301 h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[i]+Start, NROI))
302 MeanXingHist[np.mod(MeanXing[i]+Start, NROI)] += 1
303
304 Correction = (P_nom-Period)*DAMPING
305
306 # Positive correction from MeanXing[i] to MeanXing[i+1]
307 # and negative to all others
308 Time = 0
309 for j in range(NROI):
310 diff = CellTime[j+1] - CellTime[j]
311 if ((I2>I1 and j>I1 and j<=I2) or (I2<I1 and (j<=I2 or j>I1))):
312 diff += Correction/np.mod(I2-I1+NROI, NROI)
313 else:
314 diff -= Correction/(NROI-np.mod(I2-I1+NROI, NROI))
315 CellTime[j] = Time
316 Time += diff
317 CellTime[NROI] = Time
318
319 t_nom = t_nomf(CellTime)
320
321 h.list[Chip].dict["lomeanxing"].Fill(MeanXing[NumXing-1])
322 h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[NumXing-1]+Start, NROI))
323
324 return t_nom
325
326
327# create array to put the calibrated times for each cell in
328t_nom = np.ones(NROI)/fsampling
329# loop over events
330for Event in range( NEvents ):
331 print "Event ", Event
332 Chip = 0
333
334 # loop over every 9th channel in Chip
335 for Pix in range(NPix)[8::9]:
336 print "Pix ", Pix
337
338 # calculate t_nom for each pixel
339 if(alg == "Remo"):
340 MeanXing, TimeXing, NumXing = Crossing(Data[Event], np.average(Data[Event]), t_nom, int(Start[Event]), Chip)
341 h.list[Chip].dict["numxing"].Fill(NumXing)
342 t_nom = Corr(MeanXing, NumXing, TimeXing, t_nom, int(Start[Event]), Chip, Event)
343 h.list[Chip].dict["conv"].SetBinContent(Event+1, np.average(np.abs(t_nom-t_true)))
344
345
346 if(alg == "Oliver"):
347 MeanXing, TimeXing, NumXing = CrossingOliver(Data[Event], np.average(Data[Event]), t_nom, int(Start[Event]), Chip)
348 h.list[Chip].dict["numxing"].Fill(NumXing)
349 t_nom = CorrOliver(MeanXing, NumXing, TimeXing, t_nom, int(Start[Event]), Chip)
350 h.list[Chip].dict["conv"].SetBinContent(Event+1, np.average(np.abs(t_nom-t_true)))
351
352 Chip += 1
353
354CellTime = CellTimef(t_nom)
355
356print "t_tot = ", np.sum(t_nom)
357
358if(save == "yes"):
359 np.save(str(alg)+"_"+str(variante)+"_"+str(NEvents)+"x"+str(NChip)+str(tag), CellTime)
360 np.savetxt(str(alg)+"_"+str(variante)+"_"+str(NEvents)+"x"+str(NChip)+str(tag), CellTime[None], fmt="%7.4f")
361
362
363for Chip in range(NChip):
364 for i in range(NROI):
365 h.list[Chip].dict["int_dev"].SetBinContent(i+1, CellTime[i]-i/fsampling)
366 h.list[Chip].dict["diff_dev"].SetBinContent(i+1, t_nom[i]-1/fsampling)
367 h.list[Chip].dict["data"].SetBinContent(i+1, Data[-1][i])
368
369pyfact.SaveHistograms(h.list, "tcal_"+str(alg)+"_"+str(variante)+str(tag)+".root", "RECREATE")
370
371
372t_f = time()
373print "vergangene Zeit = ", (t_f - t_s)/60., " min."
374
375print "Variante = ", variante
376print "Algorithmus = ", alg
377print "Speichern = ", save
378print "NChip = ", NChip
379print "NEvents = ", NEvents
380print "Histo saved as = ", "tcal_"+str(alg)+"_"+str(variante)+str(tag)+".root"
381if(save == "yes"):
382 print "CellTime saved as = ", str(alg)+"_"+str(variante)+"_"+str(NEvents)+"x"+str(NChip)+str(tag)+".npy"
Note: See TracBrowser for help on using the repository browser.