1 | #!/usr/bin/python
|
---|
2 | #
|
---|
3 | # Werner Lustermann
|
---|
4 | # ETH Zurich
|
---|
5 | #
|
---|
6 | import pyfact
|
---|
7 | from myhisto import *
|
---|
8 | from hist import *
|
---|
9 | import numpy as np
|
---|
10 | import numpy.random as rnd
|
---|
11 | from scipy import interpolate as ip
|
---|
12 | from ROOT import *
|
---|
13 | from time import time
|
---|
14 | from optparse import OptionParser
|
---|
15 |
|
---|
16 | t_s = time()
|
---|
17 |
|
---|
18 | import os.path
|
---|
19 | import sys
|
---|
20 | #dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits'
|
---|
21 | #calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits'
|
---|
22 |
|
---|
23 |
|
---|
24 | #dfname = '/data00/fact-construction/raw/2011/11/24/20111124_113.fits.gz'
|
---|
25 | #calfname = '/data00/fact-construction/raw/2011/11/24/20111124_111.drs.fits.gz'
|
---|
26 |
|
---|
27 | dfname = '/data00/fact-construction/raw/2012/02/23/20120223_205.fits.gz'
|
---|
28 | calfname = '/data00/fact-construction/raw/2012/02/23/20120223_202.drs.fits.gz'
|
---|
29 |
|
---|
30 | if not os.path.isfile(dfname):
|
---|
31 | print dfname, 'is not a file!!!!'
|
---|
32 | sys.exit(1)
|
---|
33 | if not os.path.isfile(calfname):
|
---|
34 | print calfname, 'is not a file'
|
---|
35 | sys.exit(2)
|
---|
36 |
|
---|
37 | rd = pyfact.rawdata( dfname, calfname )
|
---|
38 |
|
---|
39 | fsampling = 2. # sampling frequency
|
---|
40 | freq = 250. # testfrequency
|
---|
41 | P_nom = 1000./freq # nominal Period due to testfrequency
|
---|
42 | DAMPING = 0.02 # Damping applied to correction
|
---|
43 | NChip = 1 # Number of Chips
|
---|
44 | NEvents = 300
|
---|
45 |
|
---|
46 |
|
---|
47 | parser = OptionParser()
|
---|
48 | parser.add_option('-e', '--extended', action = 'store_true', dest='extend', default = False )
|
---|
49 | (options, args) = parser.parse_args()
|
---|
50 |
|
---|
51 | save = "yes"
|
---|
52 |
|
---|
53 | if(options.extend):
|
---|
54 |
|
---|
55 | NChip = int(raw_input("Wieviele Chips? ([1, ... , "+str(rd.NPIX/9)+"]): "))
|
---|
56 | NEventsp = int(raw_input("Wieviele Events? ([1, ... , "+str(rd.NEvents)+"]) "))
|
---|
57 | save = raw_input("Soll CellTime gespeichert werden? (yes/no) ")
|
---|
58 |
|
---|
59 |
|
---|
60 | h = hist_list(tcalHistograms, NChip, "Chip")
|
---|
61 |
|
---|
62 |
|
---|
63 | # Functions needed for calibration
|
---|
64 |
|
---|
65 | # find negative-going crossings of mean value
|
---|
66 | def Crossing(Data, Mean, t_nom, Start, Chip):
|
---|
67 | TimeXing = np.zeros(rd.NROI)
|
---|
68 | MeanXing = np.zeros(rd.NROI)
|
---|
69 | NumXing = 0
|
---|
70 | CellTime = CellTimef(np.roll(t_nom, -Start))
|
---|
71 |
|
---|
72 | for i in range(rd.NROI-1):
|
---|
73 | if ((Data[i] > Mean) & (Data[i+1] < Mean)):
|
---|
74 | MeanXing[NumXing] = i
|
---|
75 |
|
---|
76 | FirstCell = CellTime[i]
|
---|
77 | SecondCell = CellTime[i+1]
|
---|
78 |
|
---|
79 | TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
|
---|
80 | NumXing += 1
|
---|
81 |
|
---|
82 | # calculate the frequency of the testwave
|
---|
83 | freqp = fsampling*NumXing*1000./1024
|
---|
84 | if(np.abs(freqp-freq) > 20):
|
---|
85 | print "This is not a timecalibration run!!"
|
---|
86 | raise SystemExit
|
---|
87 |
|
---|
88 | return MeanXing, TimeXing, NumXing
|
---|
89 |
|
---|
90 | def CrossingOliver(Data, Mean, t_nom, Start, Chip):
|
---|
91 | TimeXing = np.zeros(lpipe)
|
---|
92 | MeanXing = np.zeros(lpipe)
|
---|
93 | NumXing = 0
|
---|
94 | CellTime = CellTimef(t_nom)
|
---|
95 |
|
---|
96 |
|
---|
97 | for i in range(lpipe-1):
|
---|
98 | if ((Data[i] > Mean) & (Data[i+1] < Mean)):
|
---|
99 | MeanXing[NumXing] = np.mod(i+Start, lpipe)
|
---|
100 |
|
---|
101 | FirstCell = CellTime[np.mod(i+Start,lpipe)]
|
---|
102 | SecondCell = CellTime[np.mod(i+1+Start, lpipe)]
|
---|
103 |
|
---|
104 | if(SecondCell < FirstCell): SecondCell =+ lpipe/fsampl
|
---|
105 |
|
---|
106 | TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
|
---|
107 | NumXing += 1
|
---|
108 |
|
---|
109 |
|
---|
110 | return MeanXing, TimeXing, NumXing
|
---|
111 |
|
---|
112 | # Determine time between crossings and apply correction
|
---|
113 | def Corr(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
|
---|
114 |
|
---|
115 | # loop over mean crossings
|
---|
116 | for i in range(int(NumXing-1)):
|
---|
117 | I1 = MeanXing[i]
|
---|
118 | I2 = MeanXing[i+1]
|
---|
119 | Period = TimeXing[i+1] - TimeXing[i]
|
---|
120 |
|
---|
121 |
|
---|
122 | if(np.abs(1-Period/P_nom) > 0.2):
|
---|
123 | #print "Skipping zero crossing due to too large deviation of period."
|
---|
124 | h.list[Chip].dict["loleftout"].Fill(MeanXing[i])
|
---|
125 | #h.list[Chip].dict["leftout"].Fill(np.mod((MeanXing[i]+Start), rd.NROI))
|
---|
126 | continue
|
---|
127 |
|
---|
128 | #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i])
|
---|
129 | h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[i]+Start), rd.NROI))
|
---|
130 |
|
---|
131 | # calculate the frequency of the testwave
|
---|
132 | freqp = 1000./Period
|
---|
133 | h.list[Chip].dict["freq"].Fill(freqp)
|
---|
134 |
|
---|
135 | C = (P_nom-Period)*DAMPING
|
---|
136 | numnom = Period*fsampling
|
---|
137 | Correction = np.zeros(rd.NROI)
|
---|
138 | Correction[I1+1:I2+1] += C/numnom
|
---|
139 | Correction[:I1+1] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
|
---|
140 | Correction[I2+1:] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
|
---|
141 |
|
---|
142 | Correction = np.roll(Correction, Start)
|
---|
143 | t_nom += Correction
|
---|
144 | #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[-1])
|
---|
145 | h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[-1]+Start), rd.NROI))
|
---|
146 |
|
---|
147 | return t_nom
|
---|
148 |
|
---|
149 | def CorrOliver(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
|
---|
150 |
|
---|
151 | CellTime = CellTimef(t_nom)
|
---|
152 | for i in range(int(NumXing-1)):
|
---|
153 | I1 = MeanXing[i]
|
---|
154 | I2 = MeanXing[i+1]
|
---|
155 | Period = TimeXing[i+1] - TimeXing[i]
|
---|
156 | if(Period<0): Period += rd.NROI/fsampl
|
---|
157 |
|
---|
158 | if(np.abs(1-Period/P_nom) > 0.2):
|
---|
159 | #print "Skipping zero crossing due to too large deviation of period."
|
---|
160 | continue
|
---|
161 |
|
---|
162 | Correction = (P_nom-Period)*DAMPING
|
---|
163 |
|
---|
164 | # Positive correction from MeanXing[i] to MeanXing[i+1]
|
---|
165 | # and negative to all others
|
---|
166 | Time = 0
|
---|
167 | for j in range(rd.NROI):
|
---|
168 | diff = CellTime[j+1] - CellTime[j]
|
---|
169 | if ((I2>I1 and j>I1 and j<=I2) or (I2<I1 and (j<=I2 or j>I1))):
|
---|
170 | diff += Correction/np.mod(I2-I1+rd.NROI, rd.NROI)
|
---|
171 | else:
|
---|
172 | diff -= Correction/(rd.NROI-np.mod(I2-I1+rd.NROI, rd.NROI))
|
---|
173 | CellTime[j] = Time
|
---|
174 | Time += diff
|
---|
175 | CellTime[rd.NROI] = Time
|
---|
176 |
|
---|
177 | print CellTime[-1]
|
---|
178 | t_nom = t_nomf(CellTime)
|
---|
179 |
|
---|
180 | return t_nom
|
---|
181 |
|
---|
182 |
|
---|
183 |
|
---|
184 | # Summing up the nominal times for each cell to the total cell time.
|
---|
185 | def CellTimef(t_nom):
|
---|
186 | CellTime = np.zeros(rd.NROI+1)
|
---|
187 | for i in range(rd.NROI):
|
---|
188 | CellTime[i+1] = CellTime[i] + t_nom[i]
|
---|
189 |
|
---|
190 | return CellTime
|
---|
191 |
|
---|
192 | # getting t_nom from CellTime
|
---|
193 | def t_nomf(CellTime):
|
---|
194 | t_nom = np.zeros(rd.NROI)
|
---|
195 | for i in range(rd.NROI):
|
---|
196 | t_nom[i] = CellTime[i+1]-CellTime[i]
|
---|
197 |
|
---|
198 | return t_nom
|
---|
199 |
|
---|
200 |
|
---|
201 | rd.next()
|
---|
202 | NPIX = NChip*9
|
---|
203 | NEvents = NEventsp
|
---|
204 |
|
---|
205 | # create array to put the calibrated times for each cell in
|
---|
206 | t_nom = np.ones([NPIX, rd.NROI])/fsampling
|
---|
207 | # loop over events
|
---|
208 | for Event in range( NEvents ):
|
---|
209 | print "Event ", Event
|
---|
210 | Chip = 0
|
---|
211 |
|
---|
212 | # loop over every 9th channel in Chip
|
---|
213 | for Pix in range(NPIX)[8::9]:
|
---|
214 | print "Pix ", Pix
|
---|
215 |
|
---|
216 | # calculate t_nom for each pixel
|
---|
217 | MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip)
|
---|
218 | h.list[Chip].dict["numxing"].Fill(NumXing)
|
---|
219 | #h.list[Chip].dict["start"].Fill(rd.startCells[Pix])
|
---|
220 | t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip)
|
---|
221 |
|
---|
222 | Chip += 1
|
---|
223 |
|
---|
224 | rd.next()
|
---|
225 |
|
---|
226 | t_nom = t_nom[8::9]
|
---|
227 |
|
---|
228 | if(save == "yes"):
|
---|
229 | np.save(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
|
---|
230 | np.savetxt(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
|
---|
231 |
|
---|
232 | for Chip in range(NChip):
|
---|
233 | for i in range(rd.NROI):
|
---|
234 | h.list[Chip].dict["int_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][:i])-i*1/fsampling)
|
---|
235 | h.list[Chip].dict["diff_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][i] - 1/fsampling))
|
---|
236 |
|
---|
237 | pyfact.SaveHistograms(h.list, str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip)+".root", "RECREATE")
|
---|
238 |
|
---|
239 |
|
---|
240 | t_f = time()
|
---|
241 |
|
---|
242 |
|
---|
243 |
|
---|
244 | print "time ellapsed = ", (t_f - t_s)/60., " min."
|
---|
245 | print "Speichern = ", save
|
---|
246 | print "NChip = ", NChip
|
---|
247 | print "NEvents = ", NEventsp
|
---|
248 | print "Histo saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".root"
|
---|
249 | if(save == "yes"):
|
---|
250 | print "CellTime saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".npy"
|
---|
251 |
|
---|