#!/usr/bin/python # # Werner Lustermann # ETH Zurich # 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 t_s = time() import os.path import sys #dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits' #calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits' #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' dfname = '/data00/fact-construction/raw/2012/02/23/20120223_205.fits.gz' calfname = '/data00/fact-construction/raw/2012/02/23/20120223_202.drs.fits.gz' if not os.path.isfile(dfname): print dfname, 'is not a file!!!!' sys.exit(1) if not os.path.isfile(calfname): print calfname, 'is not a file' sys.exit(2) rd = pyfact.rawdata( dfname, calfname ) fsampling = 2. # sampling frequency freq = 250. # testfrequency P_nom = 1000./freq # nominal Period due to testfrequency DAMPING = 0.02 # Damping applied to correction NChip = 1 # Number of Chips NEvents = 300 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"