Changeset 13404
- Timestamp:
- 04/21/12 09:25:20 (13 years ago)
- File:
-
- 1 edited
-
fact/tools/pyscripts/sandbox/dneise/timecal.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/pyscripts/sandbox/dneise/timecal.py
r13368 r13404 8 8 # D. Neise (TU Dortmund) April 2012 9 9 # 10 import pyfact11 from myhisto import *12 from hist import *13 10 import numpy as np 14 import numpy.random as rnd15 from scipy import interpolate as ip16 from ROOT import *17 from time import time18 from optparse import OptionParser19 20 21 11 from pyfact import RawData 22 12 from drs_spikes import DRSSpikes … … 36 26 self.calib_path = cpath 37 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 38 32 self.despike = DRSSpikes() 39 self.z X = ZeroXing()33 self.zero_crossing_finder = ZeroXing(slope = -1) # -1 means falling edge 40 34 self.fsampling = 2. # sampling frequency 41 35 self.freq = 250. # testfrequency 42 36 self.P_nom = 1000./freq # nominal Period due to testfrequency 43 37 self.DAMPING = 0.02 # Damping applied to correction 44 self.NChip = 1 # Number of Chips45 self.NEvents = 300 # Number of Events to Process46 38 47 self.t_nom = np.ones([run.npix, run.nroi])/fsampling 48 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 49 44 def __call__(self): 50 45 """ the class-call method performs the actual calculation based on the input data … … 54 49 # loop over events 55 50 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 51 # in the next() method, which is called during looping, 52 # the data is automatically amplitude calibrated. 53 # we just need every 9th pixel, so this is a place 54 # for optimization 55 data = event['acal_data'] 56 57 # slice out every 9th channel 58 niners_ids = range(8,self.run.npix,9) # [8, 17, 26, ..., 1430, 1439] 59 data = data[niners_ids,:] 60 start_cell = event['start_cells'][niners_ids,:] 61 62 data = despike( data ) 63 64 # shift data down by the mean # possible in one line 65 # calculate mean 66 mean = data.mean(1) 67 # make mean in the same shape as data, for elementwise subtraction 68 mean = mean.repeat(data.shape[1]).reshape(data.shape) 69 # subtraction 70 data = data - mean 71 72 # find zero crossings with falling edge 73 # the output is a list, of lists, of slices. 74 # each slice, shows a zero crossing. in fact, is is the slice where 75 # the data is still positive ... in the next slice it is negative, (or just zero). 76 # This was called 'MeanXing' before 77 all_crossings = self.zero_crossing_finder(data) 78 79 # now we have to calculate the exact time of the crossings 80 # but in units of nanoseconds therefor we use self.time_calib 81 # rolled appropriately, and summed up in order to have the integral timing calibration 82 # 83 #the list, 'time_of_all_crossings' will contain sublists, 84 # which in turn will contain the times, in ns, for each zero-crossing 85 time_of_all_crossings = [] 86 # the list, 'float_slice_of_all_crossings' 87 # is similar to the list above, but i just contains the position of the crossings 88 # in the pipeline, but more precisely than 'all_crossings', the position is linear interpolated 89 float_slice_of_all_crossings = [] 90 91 for chip_id,crossings in enumerate(all_crossings): 92 time_of_all_crossings.append([]) 93 float_slice_of_all_crossings.append([]) 94 time = np.roll(self.time_calib[chip_id], -start_cell[chip_id]).cumsum() 95 for crossing in crossings: 96 # performing linear interpolation of time of zero crossing 97 # what fraction of a slide is the interpolated crossing away from the integral slice id 98 fraction_of_slice = data[chip_id,crossing] / (data[chip_id,crossing]-data[chip_id,crossing+1]) 99 float_slice_of_all_crossings[-1].append(fraction_of_slice + crossing) 100 slope = time[crossing+1] - time[crossing] 101 time_of_crossing = time[crossing] + slope * fraction_of_slice 102 time_of_all_crossings[-1].append(time_of_crossing) 103 104 # Now one can use those two lists, which were just created 105 # to calculate the number(float) of slices between two consecutive crossings 106 # as well as the *measured* time between two crossings. 107 # on the other hand, the period of the test signal is known to be self.P_nom 108 # so in case the measured time is longer or shorter than the known Period 109 # the associated time of all involved slices 110 # is assumed to be a bit longer or short than 1./self.fsampling 111 112 # first we make the sublists of the lists to be numpy.arrays 113 for i in len(float_slice_of_all_crossings): #both list have equal length 114 float_slice_of_all_crossings[i] = np.array(float_slice_of_all_crossings[i]) 115 time_of_all_crossings[i] = np.array(time_of_all_crossings[i]) 116 117 # now we can easily calculate the measured periods using np.diff() 118 all_measured_periods = [] 119 for chip_times in time_of_all_crossings: 120 all_measured_periods.append(np.diff(chip_times)) 121 122 for chid,chip_periods in enumerate(all_measured_periods): 70 123 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 132 parser = OptionParser() 133 parser.add_option('-e', '--extended', action = 'store_true', dest='extend', default = False ) 134 (options, args) = parser.parse_args() 135 136 save = "yes" 137 138 if(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 145 h = hist_list(tcalHistograms, NChip, "Chip") 146 147 148 # Functions needed for calibration 149 150 # find negative-going crossings of mean value 151 def 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 175 def 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 198 def 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 234 def 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. 270 def 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 278 def 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 286 rd.next() 287 NPIX = NChip*9 288 NEvents = NEventsp 289 290 # create array to put the calibrated times for each cell in 291 t_nom = np.ones([NPIX, rd.NROI])/fsampling 292 # loop over events 293 for 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 311 t_nom = t_nom[8::9] 312 313 if(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 317 for 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 322 pyfact.SaveHistograms(h.list, str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip)+".root", "RECREATE") 323 324 325 t_f = time() 326 327 328 329 print "time ellapsed = ", (t_f - t_s)/60., " min." 330 print "Speichern = ", save 331 print "NChip = ", NChip 332 print "NEvents = ", NEventsp 333 print "Histo saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".root" 334 if(save == "yes"): 335 print "CellTime saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".npy" 336 124 for i,period in enumerate(chip_periods): 125 #shortcuts to the involved slices 126 start = all_crossings[chip][i] 127 end = all_crossings[chip][i+1] 128 interval = end-start 129 rest = 1024-interval 130 131 delta_period = (self.P_nom - period) * self.DAMPING 132 # nominal_number_of_slices 133 nom_num_slices = period * self.fsampling 134 135 # now the delta_period is distributed to all slices 136 # I do not exactly understand the +1 here ... error prone!!! 137 # the following explanation assumes delta_period is positive 138 # the slices between start and end get a positive correction value 139 # This correction value is given by delta_period / nom_num_slices 140 pos_fraction = delta_period / nom_num_slices 141 # the total positive correction sums up to 142 total_positive_correction = interval * pos_fraction 143 # the total negative correction must sum up to the sampe value 144 total_negative_correction = rest * (interval * (pos_fraction/rest)) 145 # so we can call the product following 'rest' the 'neg_fraction 146 neg_fraction = -1 * interval/rest * important 147 assert total_positive_correction - total_negative_correction == 0 148 149 correction = np.zeros(1024) 150 correction[start+1:end+1] = pos_fraction 151 correction[:start+1] = neg_fraction 152 correction[end+1:] = neg_fraction 153 assert correction.sum() == 0 154 155 # now we have to add these correction values to self.time_calib 156 self.time_calib[chip,:] += np.roll(correction, start_cell[chip])
Note:
See TracChangeset
for help on using the changeset viewer.
