Index: /fact/tools/pyscripts/sandbox/dneise/timecal.py =================================================================== --- /fact/tools/pyscripts/sandbox/dneise/timecal.py (revision 13403) +++ /fact/tools/pyscripts/sandbox/dneise/timecal.py (revision 13404) @@ -8,15 +8,5 @@ # D. Neise (TU Dortmund) April 2012 # -import pyfact -from myhisto import * -from hist import * import numpy as np -import numpy.random as rnd -from scipy import interpolate as ip -from ROOT import * -from time import time -from optparse import OptionParser - - from pyfact import RawData from drs_spikes import DRSSpikes @@ -36,15 +26,20 @@ self.calib_path = cpath self.run = RawData(dpath, cpath, return_dict = True) + + if self.run.nroi != 1024: + raise TypeError('Time Calibration Data must have ROI=1024, ROI of '+dpath+' is: '+str(self.run.nroi)) + self.despike = DRSSpikes() - self.zX = ZeroXing() + self.zero_crossing_finder = ZeroXing(slope = -1) # -1 means falling edge self.fsampling = 2. # sampling frequency self.freq = 250. # testfrequency self.P_nom = 1000./freq # nominal Period due to testfrequency self.DAMPING = 0.02 # Damping applied to correction - self.NChip = 1 # Number of Chips - self.NEvents = 300 # Number of Events to Process - self.t_nom = np.ones([run.npix, run.nroi])/fsampling - + # this ndarray will contain the best estimate of the true width of each slice after a __call__ + # but it is used to determine the current best estimate inside the __call__ as well, + # since the algorithm is iteratively implemented + #self.time_calib = np.linspace( 0. , 512, num=1024, endpoint=False).repeat(run.npix/9).reshape( (1024,run.npix/9) ).transpose() + self.time_calib = np.ones( (self.run.npix/9, 1024) )/self.fsampling def __call__(self): """ the class-call method performs the actual calculation based on the input data @@ -54,283 +49,108 @@ # loop over events for event in self.run: - data = despike( event['acal_data'] ) - start_cells = event['start_cells'] - MeanXing, TimeXing, NumXing = Calculate_Crossing(data, event['start_cells']) - - - def Calculate_Crossing(self, data, start_cells): - TimeXing = np.zeros(self.run.nroi) - MeanXing = np.zeros(self.run.nroi) - NumXing = 0 - CellTime = CellTimef(np.roll(t_nom, -Start)) - - for i in range(rd.NROI-1): - if (Data[i] > Mean) & (Data[i+1] < Mean): - MeanXing[NumXing] = i + # in the next() method, which is called during looping, + # the data is automatically amplitude calibrated. + # we just need every 9th pixel, so this is a place + # for optimization + data = event['acal_data'] + + # slice out every 9th channel + niners_ids = range(8,self.run.npix,9) # [8, 17, 26, ..., 1430, 1439] + data = data[niners_ids,:] + start_cell = event['start_cells'][niners_ids,:] + + data = despike( data ) + + # shift data down by the mean # possible in one line + # calculate mean + mean = data.mean(1) + # make mean in the same shape as data, for elementwise subtraction + mean = mean.repeat(data.shape[1]).reshape(data.shape) + # subtraction + data = data - mean + + # find zero crossings with falling edge + # the output is a list, of lists, of slices. + # each slice, shows a zero crossing. in fact, is is the slice where + # the data is still positive ... in the next slice it is negative, (or just zero). + # This was called 'MeanXing' before + all_crossings = self.zero_crossing_finder(data) + + # now we have to calculate the exact time of the crossings + # but in units of nanoseconds therefor we use self.time_calib + # rolled appropriately, and summed up in order to have the integral timing calibration + # + #the list, 'time_of_all_crossings' will contain sublists, + # which in turn will contain the times, in ns, for each zero-crossing + time_of_all_crossings = [] + # the list, 'float_slice_of_all_crossings' + # is similar to the list above, but i just contains the position of the crossings + # in the pipeline, but more precisely than 'all_crossings', the position is linear interpolated + float_slice_of_all_crossings = [] + + for chip_id,crossings in enumerate(all_crossings): + time_of_all_crossings.append([]) + float_slice_of_all_crossings.append([]) + time = np.roll(self.time_calib[chip_id], -start_cell[chip_id]).cumsum() + for crossing in crossings: + # performing linear interpolation of time of zero crossing + # what fraction of a slide is the interpolated crossing away from the integral slice id + fraction_of_slice = data[chip_id,crossing] / (data[chip_id,crossing]-data[chip_id,crossing+1]) + float_slice_of_all_crossings[-1].append(fraction_of_slice + crossing) + slope = time[crossing+1] - time[crossing] + time_of_crossing = time[crossing] + slope * fraction_of_slice + time_of_all_crossings[-1].append(time_of_crossing) + + # Now one can use those two lists, which were just created + # to calculate the number(float) of slices between two consecutive crossings + # as well as the *measured* time between two crossings. + # on the other hand, the period of the test signal is known to be self.P_nom + # so in case the measured time is longer or shorter than the known Period + # the associated time of all involved slices + # is assumed to be a bit longer or short than 1./self.fsampling + + # first we make the sublists of the lists to be numpy.arrays + for i in len(float_slice_of_all_crossings): #both list have equal length + float_slice_of_all_crossings[i] = np.array(float_slice_of_all_crossings[i]) + time_of_all_crossings[i] = np.array(time_of_all_crossings[i]) + + # now we can easily calculate the measured periods using np.diff() + all_measured_periods = [] + for chip_times in time_of_all_crossings: + all_measured_periods.append(np.diff(chip_times)) + + for chid,chip_periods in enumerate(all_measured_periods): - FirstCell = CellTime[i] - SecondCell = CellTime[i+1] - - TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/Data[i]) - NumXing += 1 - - # calculate the frequency of the testwave - freqp = fsampling*NumXing*1000./1024 - if(np.abs(freqp-freq) > 20): - print "This is not a timecalibration run!!" - raise SystemExit - - return MeanXing, TimeXing, NumXing - - - for Event in range( NEvents ): - print "Event ", Event - Chip = 0 - - # loop over every 9th channel in Chip - for Pix in range(NPIX)[8::9]: - print "Pix ", Pix - - # calculate t_nom for each pixel - MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip) - h.list[Chip].dict["numxing"].Fill(NumXing) - #h.list[Chip].dict["start"].Fill(rd.startCells[Pix]) - t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip) - - Chip += 1 - - rd.next() - - t_nom = t_nom[8::9] - - if(save == "yes"): - np.save(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom) - np.savetxt(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom) - - for Chip in range(NChip): - for i in range(rd.NROI): - h.list[Chip].dict["int_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][:i])-i*1/fsampling) - h.list[Chip].dict["diff_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][i] - 1/fsampling)) - - pyfact.SaveHistograms(h.list, str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip)+".root", "RECREATE") - - - t_f = time() - - - - print "time ellapsed = ", (t_f - t_s)/60., " min." - print "Speichern = ", save - print "NChip = ", NChip - print "NEvents = ", NEventsp - print "Histo saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".root" - if(save == "yes"): - print "CellTime saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".npy" - - - -parser = OptionParser() -parser.add_option('-e', '--extended', action = 'store_true', dest='extend', default = False ) -(options, args) = parser.parse_args() - -save = "yes" - -if(options.extend): - - NChip = int(raw_input("Wieviele Chips? ([1, ... , "+str(rd.NPIX/9)+"]): ")) - NEventsp = int(raw_input("Wieviele Events? ([1, ... , "+str(rd.NEvents)+"]) ")) - save = raw_input("Soll CellTime gespeichert werden? (yes/no) ") - - -h = hist_list(tcalHistograms, NChip, "Chip") - - -# Functions needed for calibration - -# find negative-going crossings of mean value -def Crossing(Data, Mean, t_nom, Start, Chip): - TimeXing = np.zeros(rd.NROI) - MeanXing = np.zeros(rd.NROI) - NumXing = 0 - CellTime = CellTimef(np.roll(t_nom, -Start)) - - for i in range(rd.NROI-1): - if ((Data[i] > Mean) & (Data[i+1] < Mean)): - MeanXing[NumXing] = i - - FirstCell = CellTime[i] - SecondCell = CellTime[i+1] - - TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i])) - NumXing += 1 - - # calculate the frequency of the testwave - freqp = fsampling*NumXing*1000./1024 - if(np.abs(freqp-freq) > 20): - print "This is not a timecalibration run!!" - raise SystemExit - - return MeanXing, TimeXing, NumXing - -def CrossingOliver(Data, Mean, t_nom, Start, Chip): - TimeXing = np.zeros(lpipe) - MeanXing = np.zeros(lpipe) - NumXing = 0 - CellTime = CellTimef(t_nom) - - - for i in range(lpipe-1): - if ((Data[i] > Mean) & (Data[i+1] < Mean)): - MeanXing[NumXing] = np.mod(i+Start, lpipe) - - FirstCell = CellTime[np.mod(i+Start,lpipe)] - SecondCell = CellTime[np.mod(i+1+Start, lpipe)] - - if(SecondCell < FirstCell): SecondCell =+ lpipe/fsampl - - TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i])) - NumXing += 1 - - - return MeanXing, TimeXing, NumXing - -# Determine time between crossings and apply correction -def Corr(MeanXing, NumXing, TimeXing, t_nom, Start, Chip): - - # loop over mean crossings - for i in range(int(NumXing-1)): - I1 = MeanXing[i] - I2 = MeanXing[i+1] - Period = TimeXing[i+1] - TimeXing[i] - - - if(np.abs(1-Period/P_nom) > 0.2): - #print "Skipping zero crossing due to too large deviation of period." - h.list[Chip].dict["loleftout"].Fill(MeanXing[i]) - #h.list[Chip].dict["leftout"].Fill(np.mod((MeanXing[i]+Start), rd.NROI)) - continue - - #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i]) - h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[i]+Start), rd.NROI)) - - # calculate the frequency of the testwave - freqp = 1000./Period - h.list[Chip].dict["freq"].Fill(freqp) - - C = (P_nom-Period)*DAMPING - numnom = Period*fsampling - Correction = np.zeros(rd.NROI) - Correction[I1+1:I2+1] += C/numnom - Correction[:I1+1] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1))) - Correction[I2+1:] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1))) - - Correction = np.roll(Correction, Start) - t_nom += Correction - #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[-1]) - h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[-1]+Start), rd.NROI)) - - return t_nom - -def CorrOliver(MeanXing, NumXing, TimeXing, t_nom, Start, Chip): - - CellTime = CellTimef(t_nom) - for i in range(int(NumXing-1)): - I1 = MeanXing[i] - I2 = MeanXing[i+1] - Period = TimeXing[i+1] - TimeXing[i] - if(Period<0): Period += rd.NROI/fsampl - - if(np.abs(1-Period/P_nom) > 0.2): - #print "Skipping zero crossing due to too large deviation of period." - continue - - Correction = (P_nom-Period)*DAMPING - - # Positive correction from MeanXing[i] to MeanXing[i+1] - # and negative to all others - Time = 0 - for j in range(rd.NROI): - diff = CellTime[j+1] - CellTime[j] - if ((I2>I1 and j>I1 and j<=I2) or (I2I1))): - diff += Correction/np.mod(I2-I1+rd.NROI, rd.NROI) - else: - diff -= Correction/(rd.NROI-np.mod(I2-I1+rd.NROI, rd.NROI)) - CellTime[j] = Time - Time += diff - CellTime[rd.NROI] = Time - - print CellTime[-1] - t_nom = t_nomf(CellTime) - - return t_nom - - - -# Summing up the nominal times for each cell to the total cell time. -def CellTimef(t_nom): - CellTime = np.zeros(rd.NROI+1) - for i in range(rd.NROI): - CellTime[i+1] = CellTime[i] + t_nom[i] - - return CellTime - -# getting t_nom from CellTime -def t_nomf(CellTime): - t_nom = np.zeros(rd.NROI) - for i in range(rd.NROI): - t_nom[i] = CellTime[i+1]-CellTime[i] - - return t_nom - - -rd.next() -NPIX = NChip*9 -NEvents = NEventsp - -# create array to put the calibrated times for each cell in -t_nom = np.ones([NPIX, rd.NROI])/fsampling -# loop over events -for Event in range( NEvents ): - print "Event ", Event - Chip = 0 - - # loop over every 9th channel in Chip - for Pix in range(NPIX)[8::9]: - print "Pix ", Pix - - # calculate t_nom for each pixel - MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip) - h.list[Chip].dict["numxing"].Fill(NumXing) - #h.list[Chip].dict["start"].Fill(rd.startCells[Pix]) - t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip) - - Chip += 1 - - rd.next() - -t_nom = t_nom[8::9] - -if(save == "yes"): - np.save(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom) - np.savetxt(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom) - -for Chip in range(NChip): - for i in range(rd.NROI): - h.list[Chip].dict["int_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][:i])-i*1/fsampling) - h.list[Chip].dict["diff_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][i] - 1/fsampling)) - -pyfact.SaveHistograms(h.list, str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip)+".root", "RECREATE") - - -t_f = time() - - - -print "time ellapsed = ", (t_f - t_s)/60., " min." -print "Speichern = ", save -print "NChip = ", NChip -print "NEvents = ", NEventsp -print "Histo saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".root" -if(save == "yes"): - print "CellTime saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".npy" - + for i,period in enumerate(chip_periods): + #shortcuts to the involved slices + start = all_crossings[chip][i] + end = all_crossings[chip][i+1] + interval = end-start + rest = 1024-interval + + delta_period = (self.P_nom - period) * self.DAMPING + # nominal_number_of_slices + nom_num_slices = period * self.fsampling + + # now the delta_period is distributed to all slices + # I do not exactly understand the +1 here ... error prone!!! + # the following explanation assumes delta_period is positive + # the slices between start and end get a positive correction value + # This correction value is given by delta_period / nom_num_slices + pos_fraction = delta_period / nom_num_slices + # the total positive correction sums up to + total_positive_correction = interval * pos_fraction + # the total negative correction must sum up to the sampe value + total_negative_correction = rest * (interval * (pos_fraction/rest)) + # so we can call the product following 'rest' the 'neg_fraction + neg_fraction = -1 * interval/rest * important + assert total_positive_correction - total_negative_correction == 0 + + correction = np.zeros(1024) + correction[start+1:end+1] = pos_fraction + correction[:start+1] = neg_fraction + correction[end+1:] = neg_fraction + assert correction.sum() == 0 + + # now we have to add these correction values to self.time_calib + self.time_calib[chip,:] += np.roll(correction, start_cell[chip])