#!/usr/bin/python -tti # # calculation of DRS Time Calibration constants # initial implementation by Remo Dietlicher and Werner Lustermann (ETH Zurich) # based on C++ classes by Oliver Grimm (ETH Zurich) # # this python implementation was coded by # 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 from extractors import ZeroXing class CalculateDrsTimeCalibrationConstants( object ): """ basic analysis class for the calculation of DRS time aperture jitter calibration constants """ def __init__(self, dpath, cpath): """ *dpath*: data file path, file should contain a periodic signal in eahc 9th channel *cpath*: std DRS amplitude calibration path """ #dfname = '/data00/fact-construction/raw/2011/11/24/20111124_113.fits.gz' #calfname = '/data00/fact-construction/raw/2011/11/24/20111124_111.drs.fits.gz' self.data_path = dpath self.calib_path = cpath self.run = RawData(dpath, cpath, return_dict = True) self.despike = DRSSpikes() self.zX = ZeroXing() 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 def __call__(self): """ the class-call method performs the actual calculation based on the input data returns: tuple of length 160 containing tuples of length 1024 containing the length of each cell in ns the class contains self.result after the call. """ # 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 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"