Raw File
plot_accelcands.py
#!/usr/bin/env python
import glob
import os.path

import matplotlib.pyplot as plt
import matplotlib
import numpy as np

import presto
import infodata

FUDGEFACTOR = 1

class FreqInterval:
    def __init__(self, fcent, ferr, numel=1):
	self.ferr = ferr
	self.fcent = fcent
	self.flo = self.fcent - self.ferr
	self.fhi = self.fcent + self.ferr
	self.width = (self.fhi-self.flo)*FUDGEFACTOR
	self.numelements = numel

    def __contains__(self, other):
	if type(self) != type(other):
	   raise ValueError("Contains test must be made between two FreqInterval objects.")
	return ((self.flo) < other.flo < (self.fhi)) or \
		((self.flo) < other.fhi < (self.fhi)) or \
		((other.flo) < self.flo < (other.fhi)) or \
		((other.flo) < self.fhi < (other.fhi))

    def __add__(self, other):
	if type(self) != type(other):
	    raise ValueError("Addition must be between two FreqInterval objects.")
	flo = min([self.fcent-self.ferr, other.fcent-other.ferr])
	fhi = max([self.fcent+self.ferr, other.fcent+other.ferr])
	return FreqInterval((flo+fhi)/2.0, (fhi-flo)/2.0, \
			    numel=self.numelements+other.numelements)

    def __str__(self):
	return "<FreqInterval: flo=%g, fhi=%g, numelements=%d>" % \
		(self.flo, self.fhi, self.numelements)

    def zaplist_string(self):
	return "\t%f\t%f" % (self.fcent, self.width)

def main():
    freqs = []
    freqerrs = []
    filenums = []
    intervals = []

    filenum = 0
    inffiles = glob.glob("*.inf")
    for inffile in inffiles:
        accelfile = inffile[:-4]+"_ACCEL_0.cand"
   	if not os.path.exists(accelfile):
	    continue
	filenum += 1
    	rzws = presto.read_rzwcands(accelfile)
    	inf = infodata.infodata(inffile)
    	T = inf.dt*inf.N
    	for rzw in rzws:
	    freq = rzw.r/T
	    freqerr = rzw.rerr/T
	    freqs.append(freq)
	    freqerrs.append(freqerr)
	    filenums.append(filenum)

	    fint = FreqInterval(freq, freqerr)
	    # Traverse list of intervals backwards
	    for ii in range(len(intervals))[::-1]:
		if fint in intervals[ii]:
		    fint = fint+intervals[ii]
		    matchfound = True
		    del intervals[ii]
	    intervals.append(fint)

    freqs = np.array(freqs)
    freqerrs = np.array(freqerrs)
    filenums = np.array(filenums)

    plt.figure(figsize=(11,8.5))
    ebax = plt.axes((0.1, 0.1, 0.7, 0.7))
    plt.errorbar(freqs, filenums, xerr=freqerrs, fmt=None, zorder=1, ecolor='k')
    # Plot intervals worth masking
    for i in intervals:
	if i.numelements > 7:
	    r = matplotlib.patches.Rectangle((i.fcent-i.width/2.0,0), i.width, max(filenums), \
						fill=True, fc='r', ec='none', \
						alpha=0.25, zorder=-1)
	    plt.gca().add_patch(r)
	    print i.zaplist_string()
    plt.xlabel("Spin Frequency (Hz)")
    plt.ylabel("File number (index)")
    hax = plt.axes((0.8, 0.1, 0.15, 0.7), sharey=ebax)
    plt.hist(filenums, bins=max(filenums), range=(0, max(filenums)), orientation='horizontal', fc='none')
    plt.savefig("accelcands.ps", orientation="landscape", papertype="letter")
    plt.show()

if __name__=='__main__':
    main()
back to top