def loadObs(self,stackLabel):
        timestampList = self.params[stackLabel+'Sequence']
        run = self.params['run']
        sunsetDate = self.params[stackLabel+'SunsetDate']
        utcDate = self.params[stackLabel+'UtcDate']
        intTime = self.params[stackLabel+'IntTime']
        wvlLowerCutoff = self.params[stackLabel+'WvlLowerCutoff']
        wvlUpperCutoff = self.params[stackLabel+'WvlUpperCutoff']

        calTimestamp = self.params[stackLabel+'WvlTimestamp']
        print stackLabel,calTimestamp
        wvlSolnFileName = FileName(run=run,date=sunsetDate,tstamp=calTimestamp).calSoln()
        wvlCalFileName = FileName(run=run,date=self.params[stackLabel+'WvlSunsetDate'],tstamp=calTimestamp).cal()
        flatSolnFileName = FileName(run=run,date=self.params[stackLabel+'FlatCalSunsetDate'],tstamp=self.params[stackLabel+'FlatCalTimestamp']).flatSoln()

        obsFileNames = [FileName(run=run,date=sunsetDate,tstamp=timestamp).obs() for timestamp in timestampList]
        obList = [ObsFile(obsFn) for obsFn in obsFileNames]
        for ob in obList:
            ob.loadWvlCalFile(wvlSolnFileName)
            ob.loadFlatCalFile(flatSolnFileName)

        self.stackObsFileLists[stackLabel] = obList
        
        cal = ObsFile(wvlCalFileName)
        cal.loadWvlCalFile(wvlSolnFileName)
        cal.loadFlatCalFile(flatSolnFileName)
        self.stackWvlCals[stackLabel] = cal

    def findv1(self):
        populationMax=2000
        ySum = np.zeros(populationMax)
        frameSum = 'none'
        seq5 = self.s['seq5'].split()
        for seq in seq5:
            print "seq=",seq
            outfileName = "cosmicTimeList-"+seq+".pkl"
            if not os.path.exists(outfileName):
                fn = FileName(self.s['run'], 
                              self.s['sundownDate'], 
                              self.s['obsDate']+"-"+str(seq))
                cosmic = Cosmic(fn, 
                                beginTime=self.s['beginTime'],
                                endTime=self.s['endTime'],
                                loggingLevel = logging.INFO)

                fc = cosmic.findCosmics(stride=int(self.s['stride']), 
                                        threshold=int(self.s['threshold']), 
                                        populationMax=populationMax,
                                        nSigma=float(self.s['nSigma']))
                outfile = open(outfileName, "wb")
                pickle.dump(fc['cosmicTimeList'],outfile)
                pickle.dump(fc['binContents'],outfile)
                outfile.close()
                cfn = "cosmicMask-%s.h5"%seq
                ObsFile.writeCosmicIntervalToFile(fc['interval'],1.0e6, cfn,
                                                  self.s['beginTime'],
                                                  self.s['endTime'],
                                                  int(self.s['stride']),
                                                  int(self.s['threshold']),
                                                  float(self.s['nSigma']),
                                                  populationMax)

                del cosmic

def testLoadBeammap():
    '''
    Test if a remapped beammap is actually remapped
    '''
    #open an obs file from PAL2012,the sky file for hr9087
    #we'll use the embedded beammap file, which has some misplaced pixels
    run = 'PAL2012'
    date = '20121210'
    obsTimestamp = '20121211-051650'
    obsFN = FileName(run=run,date=date,tstamp=obsTimestamp)
    obsFileName = obsFN.obs()
    obs = ObsFile(obsFileName)
    
    beammapFileName = obsFN.beammap()
    
    #load the PixelMap for PAL2012 to know which pixels should be remapped
    pixMap = remapPixels.PixelMap(obsFN.pixRemap())
    pixMapSourceList,pixMapDestList = pixMap.getRemappedPix()
    
    #load the corrected beammap into the obs
    obs.loadBeammapFile(beammapFileName)

    #check that each pixel that should be moved is moved
    #by comparing the embedded beammap and the loaded corrected one
    for source,dest in zip(pixMapSourceList,pixMapDestList):
        assert obs.beamImage[dest] == obs.file.root.beammap.beamimage[source]

    obs.file.close()

def main():
    np.set_printoptions(threshold=np.nan)
    testPixelRow = 24
    testPixelCol = 17
    #obs_20120919-131142.h5,obs_20120919-131346.h5
    #create a cal file from a twilight flat
    cal = FlatCal('../../params/flatCal.dict')
    #open another twilight flat as an observation and apply a wavelength cal and the new flat cal
#    run='LICK2012'
#    obsFileName = FileName(run=run,date='20120918',tstamp='20120919-131142').flat()
#    flatCalFileName = FileName(run=run,date='20120918',tstamp='20120919-131448').flatSoln()
#    wvlCalFileName = FileName(run=run,date='20120916',tstamp='20120917-072537').calSoln()

    run = 'PAL2012'
    obsFileName = FileName(run=run,date='20121211',tstamp='20121212-140003').obs()
    flatCalFileName = FileName(run=run,date='20121210',tstamp='').flatSoln()
    wvlCalFileName = FileName(run=run,date='20121210',tstamp='20121211-133056').calSoln()
    flatCalPath = os.path.dirname(flatCalFileName)

    ob = ObsFile(obsFileName)#('obs_20120919-131142.h5')
    ob.loadWvlCalFile(wvlCalFileName)#('calsol_20120917-072537.h5')
    ob.loadFlatCalFile(flatCalFileName)#('flatsol_20120919-131142.h5')

    #plot some uncalibrated and calibrated spectra for one pixel
    fig = plt.figure()
    ax = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    print ob.getPixelCount(testPixelRow,testPixelCol)

    #flatSpectrum,wvlBinEdges = ob.getPixelSpectrum(testPixelRow,testPixelCol,weighted=False)
    spectrum,wvlBinEdges = ob.getPixelSpectrum(testPixelRow,testPixelCol,wvlStart=cal.wvlStart,wvlStop=cal.wvlStop,wvlBinWidth=cal.wvlBinWidth,weighted=False,firstSec=0,integrationTime=-1)

    weightedSpectrum,wvlBinEdges = ob.getPixelSpectrum(testPixelRow,testPixelCol,weighted=True)
    #flatSpectrum,wvlBinEdges = cal.flatFile.getPixelSpectrum(testPixelRow,testPixelCol,wvlStart=cal.wvlStart,wvlStop=cal.wvlStop,wvlBinWidth=cal.wvlBinWidth,weighted=False,firstSec=0,integrationTime=-1)
    flatSpectrum = cal.spectra[testPixelRow,testPixelCol]
    x = wvlBinEdges[0:-1]
    ax.plot(x,cal.wvlMedians,label='median spectrum',alpha=.5)
    ax2.plot(x,cal.flatFactors[testPixelRow,testPixelCol,:],label='pixel weights',alpha=.5)
    ax2.set_title('flat weights for pixel %d,%d'%(testPixelRow,testPixelCol))
    ax.plot(x,spectrum+20,label='unweighted spectrum for pixel %d,%d'%(testPixelRow,testPixelCol),alpha=.5)
    ax.plot(x,weightedSpectrum+10,label='weighted %d,%d'%(testPixelRow,testPixelCol),alpha=.5)
    ax.plot(x,flatSpectrum+30,label='flatFile %d,%d'%(testPixelRow,testPixelCol),alpha=.5)

    ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.3),fancybox=True,ncol=3)
    plt.show()

    #display a time-flattened image of the twilight flat as it is and after using itself as it's flat cal
    #cal.flatFile.loadFlatCalFile(flatCalFileName)#('flatsol_20120919-131142.h5')
    #cal.flatFile.displaySec(weighted=True,integrationTime=-1)
    #ob.displaySec(integrationTime=-1)
    #ob.displaySec(weighted=True,integrationTime=-1)


    for idx in range(0,100,20):
        factors10 = cal.flatFactors[:,:,idx]
        plt.matshow(factors10,vmax=np.mean(factors10)+1.5*np.std(factors10))
        plt.title('Flat weights at %d'%cal.wvlBinEdges[idx])
        plt.colorbar()
        plt.savefig('plots/factors%d.png'%idx)
        plt.show()

def centroidObs(obsPath,centroidPath,centroidRa,centroidDec,haOffset,xGuess,yGuess,hotPath,flatPath):
    obs = ObsFile(obsPath)
    print obsPath,obs.getFromHeader('exptime'),obs
    if not os.path.exists(hotPath):
        hp.findHotPixels(obsFile=obs,outputFileName=hotPath)
    obs.loadHotPixCalFile(hotPath,switchOnMask=False)
    obs.loadBestWvlCalFile()
    obs.loadFlatCalFile(flatPath)
    obs.setWvlCutoffs(3000,8000)
    cc.centroidCalc(obs,centroidRa,centroidDec,guessTime=300,integrationTime=30,secondMaxCountsForDisplay=2000,HA_offset=haOffset,xyapprox=[xGuess,yGuess],outputFileName=centroidPath)
    print 'done centroid',centroidPath
    del obs

    def makemovie1(self):
        run = self.s['run']
        sundownDate = self.s['sundownDate']
        obsDate = self.s['obsDate']
        stride = int(self.s['stride'])
        seq5 = self.s['seq5'].split()
        for seq in seq5:
            inFile = open("cosmicTimeList-%s.pkl"%(seq),"rb")
            cosmicTimeList = pickle.load(inFile)
            binContents = pickle.load(inFile)
            cfn = "cosmicMask-%s.h5"%seq
            intervals = ObsFile.readCosmicIntervalFromFile(cfn)
            for interval in intervals:
                print "interval=",interval
                fn = FileName(run, sundownDate,obsDate+"-"+seq)
                obsFile = ObsFile(fn.obs())
                obsFile.loadTimeAdjustmentFile(fn.timeAdjustments())
                i0=interval[0]
                i1=interval[1]
                intervalTime = i1-i0
                dt = intervalTime/2
                beginTime = max(0,i0-0.000200)
                endTime = beginTime + 0.001
                integrationTime = endTime-beginTime
                nBins = int(np.round(obsFile.ticksPerSec*(endTime-beginTime)+1))
                timeHgValues = np.zeros(nBins, dtype=np.int64)
                ymax = sys.float_info.max/100.0
                for iRow in range(obsFile.nRow):
                    for iCol in range(obsFile.nCol):
                        gtpl = obsFile.getTimedPacketList(iRow,iCol,
                                                          beginTime,integrationTime)
                        ts = (gtpl['timestamps'] - beginTime)*obsFile.ticksPerSec
                        ts64 = np.round(ts).astype(np.uint64)
                        tsBinner.tsBinner(ts64, timeHgValues)
                plt.clf()
                plt.plot(timeHgValues, label="data")
                x0 = (i0-beginTime)*obsFile.ticksPerSec
                x1 = (i1-beginTime)*obsFile.ticksPerSec
                plt.fill_between((x0,x1),(0,0), (ymax,ymax), alpha=0.2, color='red')   
                plt.yscale("symlog",linthreshy=0.9)
                plt.xlim(0,1000)
                plt.ylim(-0.1,300)
                tick0 = int(np.round(i0*obsFile.ticksPerSec))
                plotfn = "cp-%05d-%s-%s-%s-%09d"%(timeHgValues.sum(),run,obsDate,seq,tick0)
                plt.title(plotfn)
                plt.legend()
                plt.savefig(plotfn+".png")
                plt.xlabel("nSigma=%d stride=%d threshold=%d"%(int(self.s['nSigma']),int(self.s['stride']),int(self.s['threshold'])))
                print "plotfn=",plotfn
                

        os.system("convert -delay 0 `ls -r cp*png` cp.gif")

def centroidObs(obsPath,centroidPath,centroidRa,centroidDec,haOffset,xGuess,yGuess,savePath,tstamp):
    obs = ObsFile(obsPath)
#    if not os.path.exists(hotPath):
#        hp.findHotPixels(obsFile=obs,outputFileName=hotPath)
    obs.loadAllCals()
#    obs.loadHotPixCalFile(hotPath,switchOnMask=True)
#    obs.loadBestWvlCalFile()
#    obs.loadFlatCalFile(flatPath)
    obs.setWvlCutoffs(3000,11000)
    obs.loadCentroidListFile(centroidPath)
    
    ctrdFile = obs.centroidListFile
    sliceTimes = ctrdFile.root.centroidlist.times.read()
    xPositions = ctrdFile.root.centroidlist.xPositions.read()
    yPositions = ctrdFile.root.centroidlist.yPositions.read()
    intTime = sliceTimes[1]-sliceTimes[0]
    
    for iTime,time in enumerate(sliceTimes):
        x = xPositions[iTime]
        y = yPositions[iTime]
        title='centroid_{}_{}s'.format(tstamp,time)
        imgDict = obs.getPixelCountImage(firstSec=time,integrationTime=intTime,weighted=True)
        imgPath=os.path.join(savePath,title+'.png')
        pop = PopUp(showMe=False)
        pop.plotArray(imgDict['image'],title=title)
        pop.axes.plot(x,y,color='g',marker='d')
        pop.fig.savefig(imgPath)
        print 'saved to',imgPath
        
    del obs

    def openImage(self):

        timestampList = [self.params['obsUtcDate']+'-'+ts for ts in self.params['obsSequence']] 
        run = self.params['run']
        sunsetDate = self.params['obsSunsetDate']
        utcDate = self.params['obsUtcDate']
        self.intTime = self.params['intTime']
        wvlLowerCutoff = self.params['wvlLowerCutoff']
        wvlUpperCutoff = self.params['wvlUpperCutoff']

        calTimestamp = self.params['wvlTimestamp']
        wfn = FileName(run=run,date=sunsetDate,tstamp=calTimestamp).calSoln()
        calfn = FileName(run=run,date=self.params['wvlSunsetDate'],tstamp=calTimestamp).cal()
        ffn = FileName(run=run,date=self.params['flatCalSunsetDate'],tstamp='').flatSoln()

        obsFns = [FileName(run=run,date=sunsetDate,tstamp=timestamp).obs() for timestamp in timestampList]
        self.obList = [ObsFile(obsFn) for obsFn in obsFns]
        for ob in self.obList:
            print 'Loading ',ob.fullFileName
            ob.loadWvlCalFile(wfn)
            ob.loadFlatCalFile(ffn)

        self.cal = ObsFile(calfn)
        self.cal.loadWvlCalFile(wfn)
        self.cal.loadFlatCalFile(ffn)
        self.loadSpectra()

    def __init__(self, fn, beginTime=0, endTime='exptime',
                 nBinsPerSec=10, flashMergeTime=1.0, 
                 applyCosmicMask = False,
                 loggingLevel=logging.CRITICAL,
                 loggingHandler=logging.StreamHandler()):
        
        """
        Opens fileName in MKID_RAW_PATH, sets roachList
        endTime is exclusive
        """
        self.logger = logging.getLogger("cosmic")
        self.logger.setLevel(loggingLevel)
        formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
        loggingHandler.setFormatter(formatter)
        self.logger.addHandler(loggingHandler)
        self.logger.info("Cosmic:  begin init for obsFile=%s"%fn.obs())
        self.fn = fn
        self.fileName = fn.obs();
        self.file = ObsFile(self.fileName)
        # apply Matt's time fix
        timeAdjustments = self.fn.timeAdjustments()
        if os.path.exists(timeAdjustments):
            self.file.loadTimeAdjustmentFile(timeAdjustments)
        # apply Julian's time masks
        timeMaskFile = self.fn.timeMask();
        if os.path.exists(timeMaskFile):
            self.file.loadHotPixCalFile(timeMaskFile,switchOnMask=True)

        # apply standard mask
        if applyCosmicMask:
            self.file.loadStandardCosmicMask()

        self._setRoachList()
        self._setAllSecs()
        self.exptime = self.file.getFromHeader('exptime')
        if endTime =='exptime':
            self.endTime = float(self.exptime)
        else:
            self.endTime = float(endTime)
        if ( (self.endTime > self.exptime) or (endTime < 0)):
            raise RuntimeError("bad endTime:  endTime=%s exptime=%s" % \
                                   (str(endTime),str(self.exptime)))

        self.beginTime = float(beginTime)
        self.timeHgs = "none"
        self.nBinsPerSec = nBinsPerSec
        self.flashMergeTime = flashMergeTime

        self.times = \
        np.arange(self.beginTime, self.endTime, 1.0/self.nBinsPerSec)

        # for measuring flashes, indexed by roach name
        self.rMean = {}   # mean from meanclip
        self.rSigma = {}  # sigma from meanclip
        self.rNSurvived = {} # number of survivors from meanclip
        self.rNormed = {} # (value-mean)/sigma
        self.flashInterval = {}
        self.logger.info("Cosmic:  end init:  beginTime=%s endTime=%s"%(str(self.beginTime),str(self.endTime)))

def timePL(tstamp,obsPath,centroidPath):
    obs = ObsFile(obsPath)
    obs.loadAllCals()
    obs.setWvlCutoffs(3000,11000)
    obs.loadCentroidListFile(centroidPath)
    writePhotonList(obs,photListDescription=PulsarPhotonList,checkForExisting=False)
    del obs

def main():

    objectName = "hz21"
    fileNum=0
    energyBinWidth = 0.1
    #make bins for 3000 to 13000
    wvlStart = 3000
    wvlStop = 13000
    wvlBinEdges = ObsFile.makeWvlBins(energyBinWidth,wvlStart,wvlStop)
    nWvlBins = len(wvlBinEdges)-1
    binWidths = np.empty(nWvlBins)
    print "Showing bin widths for %i bins"%(nWvlBins)
    for i in xrange(nWvlBins):
        binWidths[i] = wvlBinEdges[i+1]-wvlBinEdges[i]
    print binWidths

    nVirtPixX=250
    nVirtPixY=250
    cube = np.zeros((nVirtPixX,nVirtPixY,nWvlBins),dtype=float)

    for n in xrange(nWvlBins):
        print "Making image for wvls %i to %i"%(wvlBinEdges[n], wvlBinEdges[n+1])
        virtualImage, imageStack, medImage = makeImageStack(fileNames='photons_*.h5', dir=os.getenv('MKID_PROC_PATH', default="/Scratch")+'/photonLists/20131209',
                   detImage=False, saveFileName='stackedImage.pkl', wvlMin=wvlBinEdges[n],
                   wvlMax=wvlBinEdges[n+1], doWeighted=True, medCombine=False, vPlateScale=0.2,
                   nPixRA=nVirtPixX,nPixDec=nVirtPixY)
        print virtualImage
        print virtualImage.image
        print np.shape(virtualImage.image)
        cube[:,:,n] = virtualImage.image
    
    #calculate midpoints of wvl bins for plotting
    wvls = np.empty((nWvlBins),dtype=float)
    for n in xrange(nWvlBins):
        binsize=wvlBinEdges[n+1]-wvlBinEdges[n]
        wvls[n] = (wvlBinEdges[n]+(binsize/2.0))

    print "wvls ",wvls
    #reshape cube for makeMovie
    movieCube = np.zeros((nWvlBins,np.shape(cube)[0],np.shape(cube)[1]),dtype=float)
    for i in xrange(nWvlBins):
        movieCube[i,:,:] = cube[:,:,i]
        #show individual frames as they are made to debug
        #plt.matshow(movieCube[i],vmin = 0, vmax = 100)
        #plt.show()
    print "movieCube shape ", np.shape(movieCube)
    print "wvls shape ", np.shape(wvls)

    #print cube
    #print "--------------------------"
    #print movieCube

    np.savez('%s_raw_%s.npz'%(objectName,fileNum),stack=movieCube,wvls=wvls)
    utils.makeMovie(movieCube,frameTitles=wvls,cbar=True,outName='%s_pl_raw_%s.gif'%(objectName,fileNum), normMin=0, normMax=1000)

    '''

def centroidObs(obsPath,centroidPath,centroidRa,centroidDec,haOffset,xGuess,yGuess):
    obs = ObsFile(obsPath)
    print obsPath,obs.getFromHeader('exptime'),obs
    obs.loadAllCals()
#    obs.loadBestWvlCalFile()
#    obs.loadFlatCalFile(flatPath)
    obs.setWvlCutoffs(3000,11000)
#    if not os.path.exists(hotPath):
#        hp.findHotPixels(obsFile=obs,outputFileName=hotPath,display=True,fwhm=2.,boxSize=5, nSigmaHot=4.0,)
#    obs.loadHotPixCalFile(hotPath,switchOnMask=True)
    cc.centroidCalc(obs,centroidRa,centroidDec,guessTime=300,integrationTime=30,secondMaxCountsForDisplay=2000,HA_offset=haOffset,xyapprox=[xGuess,yGuess],outputFileName=centroidPath,usePsfFit=True,radiusOfSearch=8)
    print 'done centroid',centroidPath
    del obs

def main():
    """
    params = []
    paramfile = sys.argv[1]

    f = open(paramfile,'r')
    for line in f:
        params.append(line)
    f.close()
    datadir = params[0].split('=')[1].strip()
    flatdir = params[1].split('=')[1].strip()
    fluxdir = params[2].split('=')[1].strip()
    wvldir = params[3].split('=')[1].strip()
    obsfile = params[4].split('=')[1].strip()
    skyfile = params[5].split('=')[1].strip()
    flatfile = params[6].split('=')[1].strip()
    fluxfile = params[7].split('=')[1].strip()
    wvlfile = params[8].split('=')[1].strip()
    objectName = params[9].split('=')[1].strip()
    fluxCalObject = params[10].split('=')[1].strip()

    obsFileName = os.path.join(datadir, obsfile)
    skyFileName = os.path.join(datadir, skyfile)
    wvlCalFileName = os.path.join(wvldir, wvlfile)
    flatCalFileName = os.path.join(flatdir, flatfile)
    fluxCalFileName = os.path.join(fluxdir, fluxfile)
    """

    if len(sys.argv) > 2:
        fileNum = str(sys.argv[2])
    else:
        fileNum = "0"

    # science object parameter file
    params = []
    paramfile = sys.argv[1]
    f = open(paramfile, "r")
    for line in f:
        params.append(line)
    f.close()

    datadir = params[0].split("=")[1].strip()
    flatdir = params[1].split("=")[1].strip()
    wvldir = params[2].split("=")[1].strip()
    obsfile = params[3].split("=")[1].strip()
    skyfile = params[4].split("=")[1].strip()
    flatfile = params[5].split("=")[1].strip()
    wvlfile = params[6].split("=")[1].strip()
    objectName = params[9].split("=")[1].strip()

    # wvldir = "/Scratch/waveCalSolnFiles/oldbox_numbers/20121205"
    # objectName = "crabNight1"

    if len(params) > 10:
        xpix = int(params[10].split("=")[1].strip())
        ypix = int(params[11].split("=")[1].strip())
        apertureRadius = int(params[12].split("=")[1].strip())
        startTime = int(params[13].split("=")[1].strip())
        intTime = int(params[14].split("=")[1].strip())

    obsFileName = os.path.join(datadir, obsfile)
    skyFileName = os.path.join(datadir, skyfile)
    wvlCalFileName = os.path.join(wvldir, wvlfile)
    flatCalFileName = os.path.join(flatdir, flatfile)

    obs = ObsFile(obsFileName)
    # obs.loadWvlCalFile(wvlCalFileName)
    obs.loadBestWvlCalFile()
    obs.loadFlatCalFile(flatCalFileName)
    print "analyzing file %s" % (obsFileName)
    print "loaded data file and calibrations\n---------------------\n"

    nRow = obs.nRow
    nCol = obs.nCol
    obsTime = obs.getFromHeader("exptime")
    # wvlBinEdges,obsSpectra = loadSpectra(obs,nCol,nRow)
    # nWvlBins=len(wvlBinEdges)-1

    # print np.shape(obsSpectra)
    # print nRow
    # print nCol
    # print nWvlBins

    # Apply Hot pixel masking before getting dead time correction
    # HotPixFile = getTimeMaskFileName(obsFileName)
    HotPixFile = FileName(obsFile=obs).timeMask()
    print "making hot pixel file ", HotPixFile
    if not os.path.exists(HotPixFile):  # check if hot pix file already exists
        hp.findHotPixels(inputFileName=obsFileName, outputFileName=HotPixFile)
        print "Flux file pixel mask saved to %s" % (HotPixFile)
    obs.loadHotPixCalFile(HotPixFile)
    print "Hot pixel mask loaded %s" % (HotPixFile)

    # GET RAW PIXEL COUNT IMAGE TO CALCULATE CORRECTION FACTORS
    print "Making raw cube to get dead time correction"
    cubeDict = obs.getSpectralCube(firstSec=startTime, integrationTime=intTime, weighted=False)

    cube = np.array(cubeDict["cube"], dtype=np.double)
    wvlBinEdges = cubeDict["wvlBinEdges"]
    effIntTime = cubeDict["effIntTime"]
#.........这里部分代码省略.........

def main():

    run = 'PAL2014'
    year = '2014'
    initialPath = '/Scratch'
    packetMasterLogDir = '/LABTEST/PacketMasterLogs/'
    initialPath = os.path.join(initialPath,run)
    outPath = FileName(run=run,mkidDataDir='/Scratch/').timeAdjustments()
    outFile = tables.openFile(outPath,'w')
    timeAdjustGroup = outFile.createGroup('/','timeAdjust','Times to add to timestamps')
    firmwareDelayTable = outFile.createTable(timeAdjustGroup,'firmwareDelay',firmwareDelayDescription,'Times to add to all timestamps taken with a firmware bof')
    newFirmwareEntry = firmwareDelayTable.row
    newFirmwareEntry['firmwareName']='chan_svf_2014_Aug_06_1839.bof'
    newFirmwareEntry['firmwareDelay']=-41e-6 #s, subtract 41 us from timestamps
    newFirmwareEntry.append()
    firmwareDelayTable.flush()
    firmwareDelayTable.close()
    roachDelayTable = outFile.createTable(timeAdjustGroup,'roachDelays',roachDelaysDescription,'Times to add to each roach\'s timestamps')


    for sunsetDatePath in sorted(glob.glob(os.path.join(initialPath,year+'*'))):
        sunsetDate = os.path.basename(sunsetDatePath)
        for fullObsPath in sorted(glob.glob(os.path.join(sunsetDatePath,'obs*.h5'))):
            obsFileName = os.path.basename(fullObsPath)
            obsTStamp = obsFileName.split('.')[0].split('_')[1]
            print obsFileName
        

            obsFN = FileName(run=run,date=sunsetDate,tstamp=obsTStamp,packetMasterLogDir=packetMasterLogDir)
            pmLogFileName = obsFN.packetMasterLog()
            try:
                ob = ObsFile(fullObsPath)
            except:
                print 'can\'t open file'
                continue
                
            try:
                if os.path.getsize(pmLogFileName) <= 0:
                    continue
            except:
                print 'can\'t open Packet Master Log ',pmLogFileName
                continue

            try:
                f = open(pmLogFileName,'r')
            except:
                print 'can\'t open Packet Master Log ',pmLogFileName
                continue
            lastTstampLines = np.zeros(8)
            firstTstampLines = np.zeros(8)

            tstampLine = ''
            for line in f:
                if 'here' in line:
                    #skip over large log files with debug info
                    print 'skipping file with "here"'
                    continue
                if line.split(' ')[0] == 'bundle':
                    tstampLine = line
                    break

            if tstampLine == '':
                print 'skipping file without "bundle"'
                #didn't find lines with 'bundle' in them
                continue

            f.seek(0)
            for line in f:
                if line.split(' ')[0] == 'bundle':
                    try:
                        at = float(line.split('at')[1].split()[0])
                    except:
                        break
                    lastTstampLine = at
                    iRoach = int(line.split('roach')[1].split('took')[0].strip())
                    lastTstampLines[iRoach] = at
                    if firstTstampLines[iRoach] == 0:
                        firstTstampLines[iRoach] = at
                    
            packetReceivedUnixTimestamp = float((tstampLine.split('took')[1].split('total')[0].strip()))
            firstPacketDelay = packetReceivedUnixTimestamp-int(ob.getFromHeader('unixtime'))
                

            roachSecDelays =np.array(np.floor(lastTstampLines+firstPacketDelay-ob.getFromHeader('exptime')),dtype=np.int)
            print 'roach delays',roachSecDelays
            if np.all(roachSecDelays >= 0):
                newEntry = roachDelayTable.row
                newEntry['obsFileName'] = os.path.basename(fullObsPath)
                newEntry['roachDelays'] = roachSecDelays
                newEntry.append()
                roachDelayTable.flush()
            else:
                print 'obs was aborted midway, delays cannot be calculated'
    roachDelayTable.close()
    outFile.close()

inCounts = []
for line in powermeter:
    if '#' in line:
        continue
    inCounts.append(float(line.strip()))
inCounts = np.array(inCounts)

pixToUse = [[0,26],[2,25],[2,27],[2,31],[2,34],[2,39],[3,17],[3,26],[3,32],
            [3,36],[4,26],[4,32],[5,32],[6,22],[6,29],[6,30],[7,25],[7,27],[7,28],
            [7,32],[7,37],[8,30],[10,32],[10,42],[11,30],[12,29],[12,31],[13,29],[13,31],
            [14,27],[14,32],[15,25],[15,26],[15,30],[15,31],[15,32],[16,23],[16,30],
            [17,23],[17,25],[17,26],[17,39],[18,24],[18,28],[19,23],[19,26],[19,28],[19,30]]

cmap = matplotlib.cm.jet
obsFileName = '/Scratch/linearityTestData/obs_20130612-003423.h5'
obs = ObsFile(obsFileName)


row,col=pixToUse[0]
allResolutions = []
allModResolutions = []
allModResolutions2 = []
allCountRates = []
timeBinStarts = timeBinStarts[0:-3]
for (row,col) in pixToUse[0:3]:
    firstSec=timeBinStarts[0]
    intTime=timeBinWidth
    resolutions = []
    modResolutions = []
    modResolutions2 = []
    countRates = []

    def findCosmics(self, stride=10, threshold=100, 
                    populationMax=2000, nSigma=5, writeCosmicMask=False,
                    ppsStride=10000):
        """
        Find cosmics ray suspects.  Histogram the number of photons
        recorded at each timeStamp.  When the number of photons in a
        group of stride timeStamps is greater than threshold in second
        iSec, add (iSec,timeStamp) to cosmicTimeLists.  Also keep
        track of the histogram of the number of photons per stride
        timeStamps.

        return a dictionary of 'populationHg', 'cosmicTimeLists',
        'binContents', 'timeHgValues', 'interval', 'frameSum', and 'pps'
      
         populationHg is a histogram of the number of photons in each time bin.
        This is a poisson distribution with a long tail due to cosmic events

        
        cosmicTimeLists is a numpy array  of all the sequences that are
        suspects for cosmic rays


        binContents corresponds to cosmicTimeLists.  For each time in
        cosmicTimeLists, binContents is the number of photons detected
        at that time.
        
        timeHgValues is a histogram of the number of photons in each time
        interval
        
        frameSum is a two dimensional  numpy array of the number of photons
        detected by each pixel

        interval is the interval of data to be masked out

        pps is photons per second, calculated every ppsStride bins.

        """

        self.logger.info("findCosmics: begin stride=%d threshold=%d populationMax=%d nSigma=%d writeCosmicMask=%s"%(stride,threshold,populationMax,nSigma,writeCosmicMask))
        exptime = self.endTime-self.beginTime
        nBins = int(np.round(self.file.ticksPerSec*exptime+1))
        bins = np.arange(0, nBins, 1)
        timeHgValues,frameSum = self.getTimeHgAndFrameSum(self.beginTime,self.endTime)
        remainder = len(timeHgValues)%ppsStride
        if remainder > 0:
            temp = timeHgValues[:-remainder]
        else:
            temp = timeHgValues
        ppsTime = (ppsStride*self.file.tickDuration)
        pps = np.sum(temp.reshape(-1, ppsStride), axis=1)/ppsTime
        self.logger.info("findCosmics:  call populationFromTimeHgValues")
        pfthgv = Cosmic.populationFromTimeHgValues\
            (timeHgValues,populationMax,stride,threshold)
        #now build up all of the intervals in seconds
        self.logger.info("findCosmics:  build up intervals:  nCosmicTime=%d"%len(pfthgv['cosmicTimeList']))
        i = interval()
        iCount = 0
        secondsPerTick = self.file.tickDuration
        for cosmicTime in pfthgv['cosmicTimeList']:
            #t0 = max(0,self.beginTime+(cosmicTime-50)/1.e6)
            #t1 = min(self.endTime,self.beginTime+(cosmicTime+50)/1.e6)
            #intTime = t1-t0            
            t0 = self.beginTime+cosmicTime*secondsPerTick
            dt = stride*secondsPerTick
            t1 = t0+dt
            left = max(self.beginTime, t0-nSigma*dt)
            right = min(self.endTime, t1+2*nSigma*dt)
            i = i | interval[left,right]
            self.logger.debug("findCosmics:  iCount=%d t0=%f t1=%f left=%f right=%f"%(iCount,t0,t1,left,right))
            iCount+=1

        tMasked = Cosmic.countMaskedBins(i)
        ppmMasked = 1000000*tMasked/(self.endTime-self.beginTime)

        retval = {}
        retval['timeHgValues'] = timeHgValues
        retval['populationHg'] = pfthgv['populationHg']
        retval['cosmicTimeList'] = pfthgv['cosmicTimeList']
        retval['binContents'] = pfthgv['binContents']
        retval['frameSum'] = frameSum
        retval['interval'] = i
        retval['ppmMasked'] = ppmMasked
        retval['pps'] = pps
        retval['ppsTime'] = ppsTime
        if writeCosmicMask:
            cfn = self.fn.cosmicMask()
            self.logger.info("findCosmics:  write masks to =%s"%cfn)
            ObsFile.writeCosmicIntervalToFile(i, self.file.ticksPerSec, 
                                              cfn,self.beginTime, self.endTime, 
                                              stride, threshold, nSigma, populationMax)
        self.logger.info("findCosmics:  end with ppm masked=%d"%ppmMasked)
        return retval

NumFrames = 31
#nFiles = 13 #j0926
#nFiles=25 #crab
nFiles = 1 #try with only 1 file
curves = np.zeros((nFiles,NumFrames),dtype=float)

for k in xrange(nFiles):
    #FileName = '/home/srmeeker/scratch/standards/crabNight2_fit_%s.npz'%(fileNum)
    FileName = '/home/srmeeker/scratch/standards/crabNight1_fit_%s.npz'%(fileNum)
    print FileName
    t = np.load(FileName)

    energyBinWidth = 0.1
    wvlStart = 3000
    wvlStop = 13000
    wvlBinEdges = ObsFile.makeWvlBins(energyBinWidth,wvlStart,wvlStop)
    nWvlBins = len(wvlBinEdges)-1
    binWidths = np.empty(nWvlBins)
    for i in xrange(nWvlBins):
        binWidths[i] = wvlBinEdges[i+1]-wvlBinEdges[i]
    #print binWidths
    params = t['params']
    wvls = t['wvls']
    amps = params[:,1]
    widths = params[:,4]
    xpos = params[:,2]
    ypos = params[:,3]

    #print len(wvls)
    #print len(binWidths)

sunsetDate = '20121208'
utcDate = '20121209'
# Specify which wavelength calibration file to use.
calTimestamp = '20121209-131132'
# Create wavelength and flat cal file names
wvlCalFilename = FileName(run=run,date=sunsetDate,tstamp=calTimestamp).calSoln()
flatCalFilename = FileName(run=run,date=sunsetDate,tstamp=calTimestamp).flatSoln()
# No twilights taken on December 8, using December 7 twilights to do flat cal instead.
flatCalFilename = '/Scratch/flatCalSolnFiles/20121207/flatsol_20121207.h5'

# Run standard ObsFile functions.
# Create ObsFile instance.
tic = time()
print 'Loading obs file and performing calibrations ...'
obsFn = FileName(run=run,date=sunsetDate,tstamp=obsTimestamp).obs()
ob = ObsFile(obsFn)
# Load roach time delay corrections.
ob.loadTimeAdjustmentFile(FileName(run=run).timeAdjustments())
# Search for time mask for given observation file.  If the time mask does not exist, create it.
index1 = obsFn.find('_')
hotPixFn = '/Scratch/timeMasks/timeMask' + obsFn[index1:]
if not os.path.exists(hotPixFn):
    hp.findHotPixels(obsFn,hotPixFn)
    print "Flux file pixel mask saved to %s"%(hotPixFn)        
# Load time mask, wavelength calibration, and flat calibration and set wavelenth cutoffs.
ob.loadHotPixCalFile(hotPixFn,switchOnMask=True)
ob.loadWvlCalFile(wvlCalFilename)
ob.loadFlatCalFile(flatCalFilename)
ob.setWvlCutoffs(3000,5000)
print 'Total load time: ' + str(time()-tic) + 's'

    sunsetDate = sunsetDates[iSeq]

    for position in ObjPosFile:
        if iSeq + 1 == position[0]:
            print "finding aperture and sky masks"
            guessX = position[1]
            guessY = position[2]
            apertureMask = circleAperture(guessX, guessY, radius1)
            bigMask = circleAperture(guessX, guessY, radius2)
            skyMask = bigMask - apertureMask
    #            print skyMask[0:20][5:25]

    for i, ts in enumerate(timestampList):
        print "loading", ts
        obsFn = FileName(run=run, date=sunsetDate, tstamp=ts).obs()
        ob = ObsFile(obsFn)
        ob.loadTimeAdjustmentFile(FileName(run=run).timeAdjustments())
        index1 = obsFn.find("_")
        hotPixFn = "/Scratch/timeMasks/timeMask" + obsFn[index1:]
        if not os.path.exists(hotPixFn):
            hp.findHotPixels(obsFn, hotPixFn)
            print "Flux file pixel mask saved to %s" % (hotPixFn)
        ob.loadHotPixCalFile(hotPixFn, switchOnMask=True)
        ob.loadWvlCalFile(wfn)
        ob.loadFlatCalFile(ffn)
        ob.setWvlCutoffs(wvlLowerCutoff, wvlUpperCutoff)

        bad_solution_mask = np.zeros((46, 44))
        bad_count = 0
        for y in range(46):
            for x in range(44):

 #cosmic = Cosmic(fn, endTime='exptime')
 cosmic = Cosmic(fn, beginTime=0)
 fc = cosmic.findCosmics(stride=stride, 
                         threshold=threshold, 
                         populationMax=populationMax,
                         nSigma=nSigma)
 if frameSum == 'none':
     frameSum = fc['frameSum']
 else:
     frameSum += fc['frameSum']
 outfile = open("cosmicTimeList-"+seq+".pkl", "wb")
 pickle.dump(fc['cosmicTimeList'],outfile)
 pickle.dump(fc['binContents'],outfile)
 outfile.close()
 cfn = "cosm