def add_circle(self,ra,dec,radius,pmax,type='radial'):
     logger.info('add a circle at ra=%f, dec=%f, radius=%f. PMAX=%f'%(ra,dec,rad,pmax))
     for i in range(self.nxpix):
         for j in range(self.nxpix):
             # radial:
             #posx = nxpix/2.-i
             #posy = -nxpix/2.+j
             # circular:
             #p_x   = i #-nxpix+j
             #p_y   = j #-nxpix+i
             world = self.w.wcs_pix2world([[i,j]], 0)
             w_ra  = world[0][0]          
             w_dec = world[0][1]
             dx    =  (w_ra-ra)*numpy.cos(numpy.deg2rad(dec)) # Effect of the projection
             dy    =  w_dec-dec
             p_x = -dy
             p_y = +dx
             dist = numpy.sqrt(dx*dx+dy*dy)
             #print w_ra, w_dec, ra, dec, dist, radius
             if (dist>radius):
                 p_x=0
                 p_y=0
                 pass
             self.pol_x[i,j]=p_x
             self.pol_y[i,j]=p_y            
             pass
         pass
     pol_deg=numpy.sqrt(self.pol_x*self.pol_x+self.pol_y*self.pol_y)
     self.pol_x*=pmax/pol_deg.max()
     self.pol_y*=pmax/pol_deg.max()
     pass

def analyze():
    """Analyze the data.
    """
    logger.info('Opening output file %s...' % ANALYSIS_FILE_PATH)
    analysis_file = open(ANALYSIS_FILE_PATH, 'w')
    for i, (_min, _max) in enumerate(zip(TIME_BINNING[:-1],
                                         TIME_BINNING[1:])):
        _mcube = xBinnedModulationCube(_mcube_file_path(i))
        _mcube.fit()
        _fit_results = _mcube.fit_results[0]
        _sel_file = xEventFile(_sel_file_path(i))
        _time = numpy.average(_sel_file.event_data['TIME'])
        _sel_file.close()
        _time_errp = _max - _time
        _time_errm = _time - _min
        _pol_deg = _fit_results.polarization_degree
        _pol_deg_err = _fit_results.polarization_degree_error
        _pol_angle = _fit_results.phase
        _pol_angle_err = _fit_results.phase_error
        _spec_fitter = PIPELINE.xpxspec(_pha1_file_path(i), plot=False)
        (_index, _index_err), (_norm, _norm_err) = _spec_fitter.fit_parameters()
        # The division by the phase interval is a workaround and we should
        # keep track of that in xpselect.
        _norm /= (_max - _min)
        _norm_err /= (_max - _min)
        _data = (_time, _time_errp, _time_errm, _pol_deg, _pol_deg_err,
                 _pol_angle, _pol_angle_err, _index, _index_err, _norm,
                 _norm_err)
        _fmt = ('%.4e   ' * len(_data)).strip()
        _fmt = '%s\n' % _fmt
        _line = _fmt % _data
        analysis_file.write(_line)
    analysis_file.close()

    def create(self):
        """Overloaded method.
        """
        self.xref  = self.get('xref')
        self.yref  = self.get('yref')
        self.nxpix = self.get('nxpix')
        self.binsz = self.get('binsz') # Dimension of the pimap (degrees)
        self.proj = self.get('proj')

        
        sidex = self.nxpix*self.binsz
        logger.info('Output image dimensions are %.1f x %.1f arcmin.' %\
                    (sidex*60, sidex*60))
                    
        logger.info('Center of the image is in R.A.=%.3f Dec.=%.3f' % (self.xref,self.yref))                    
        # Build the WCS object
        self.w = wcs.WCS(naxis=2)
        self.w.wcs.crpix = [(self.nxpix+1)/2,(self.nxpix+1)/2]
        self.w.wcs.cdelt = [-self.binsz, self.binsz]
        self.w.wcs.crval = [self.xref, self.yref]
        self.w.wcs.ctype = ['RA---%s' % self.proj, 'DEC--%s' % self.proj]
        self.w.wcs.equinox = 2000.0
        #w.wcs.radesys = 'ICRS'
        self.header = self.w.to_header()
        self.pol_x=numpy.zeros((self.nxpix,self.nxpix))
        self.pol_y=numpy.zeros((self.nxpix,self.nxpix))
        pass

 def __init__(self, file_path):
     """Constructor.
     """
     assert(file_path.endswith('.fits'))
     logger.info('Opening input binned file %s...' % file_path)
     self.hdu_list = fits.open(file_path)
     self.hdu_list.info()

 def bin_(self):
     """Overloaded method.
     """
     evt_header = self.event_file.hdu_list['PRIMARY'].header
     num_chans = evt_header['DETCHANS']
     total_time = self.event_file.total_good_time()
     binning = numpy.linspace(-0.5, num_chans -0.5, num_chans+1)
     n, bins = numpy.histogram(self.event_data['PHA'], bins=binning)
     primary_hdu = self.build_primary_hdu()
     data = [numpy.arange(num_chans),
             n/total_time,
             numpy.sqrt(n)/total_time
     ]
     spec_hdu = xBinTableHDUPHA1(data)
     spec_hdu.setup_header(self.event_file.primary_keywords())
     irf_name = evt_header['IRFNAME']
     keywords = [('EXPOSURE', total_time, 'exposure time'),
                 ('RESPFILE', irf_file_path(irf_name, 'rmf')),
                 ('ANCRFILE', irf_file_path(irf_name, 'arf'))]
     spec_hdu.setup_header(keywords)
     hdu_list = fits.HDUList([primary_hdu, spec_hdu])
     hdu_list.info()
     logger.info('Writing binned PHA1 data to %s...' % self.get('outfile'))
     hdu_list.writeto(self.get('outfile'), clobber=True)
     logger.info('Done.')

 def bin_(self):
     """Overloaded method.
     """
     evt_header = self.event_file.hdu_list['PRIMARY'].header
     if self.get('mc'):
         energy = self.event_data['MC_ENERGY']
     else:
         energy = self.event_data['ENERGY']
     phi = self.event_data['PE_ANGLE']
     counts, xedges, yedges = numpy.histogram2d(energy, phi,
                                                bins=self.make_binning())
     primary_hdu = self.build_primary_hdu()
     emin, emax = xedges[:-1], xedges[1:]
     emean = []
     for _emin, _emax in zip(emin, emax):
         emean.append(numpy.mean(energy[(energy > _emin)*(energy < _emax)]))
     data = [emin, emax, emean, counts]
     xBinTableHDUMCUBE.set_phi_spec(self.get('phibins'))
     mcube_hdu = xBinTableHDUMCUBE(data)
     mcube_hdu.setup_header(self.event_file.primary_keywords())
     gti_hdu = self.event_file.hdu_list['GTI']
     hdu_list = fits.HDUList([primary_hdu, mcube_hdu, gti_hdu])
     hdu_list.info()
     logger.info('Writing binned MCUBE data to %s...' % self.get('outfile'))
     hdu_list.writeto(self.get('outfile'), clobber=True)
     logger.info('Done.')

 def add_circle(self,ra,dec,radius,pmax,ptype='circular',angle=0.0):
     ''' This implement the case of a circular shape'''
     logger.info('add a circle at ra=%f, dec=%f, radius=%f. PMAX=%f'%(ra,dec,rad,pmax))
     for i in range(self.nxpix):
         for j in range(self.nxpix):
             world = self.w.wcs_pix2world([[i,j]], 0)
             w_ra  = world[0][0]          
             w_dec = world[0][1]
             dx    =  (w_ra-ra)*numpy.cos(numpy.deg2rad(dec)) # Effect of the projection
             dy    =  w_dec-dec
             if ptype is 'linear':
                 p_x = numpy.cos(numpy.deg2rad(angle))
                 p_y = numpy.sin(numpy.deg2rad(angle))
             elif ptype is 'circular':
                 p_x = -dy
                 p_y = +dx
             elif ptype is 'radial':
                 p_x = dx
                 p_y = dy
                 pass
             dist = numpy.sqrt(dx*dx+dy*dy)
             #print w_ra, w_dec, ra, dec, dist, radius
             if (dist>radius):
                 p_x=0
                 p_y=0
                 pass
             self.pol_x[i,j]=self.pol_x[i,j]+p_x
             self.pol_y[i,j]=self.pol_y[i,j]+p_y            
             pass
         pass
     pol_deg=numpy.sqrt(self.pol_x*self.pol_x+self.pol_y*self.pol_y)
     self.pol_x*=pmax/pol_deg.max()
     self.pol_y*=pmax/pol_deg.max()
     pass

def plot_swift_lc(grb_list,show=True):
    """Plots Swift GRB light curves.
    """
    plt.figure(figsize=(10, 8), dpi=80)
    plt.title('Swift XRT light curves')
    num_grb = 0
    for grb_name in grb_list:
        flux_outfile = download_swift_grb_lc_file(grb_name, min_obs_time=21600)
        if flux_outfile is not None:
            integral_flux_spline = parse_light_curve(flux_outfile)
            if integral_flux_spline is not None:
                if grb_name == 'GRB 130427A':
                    integral_flux_spline.plot(num_points=1000,logx=True,\
                                              logy=True,show=False,\
                                              color="red",linewidth=1.0)
                    num_grb += 1
                else:
                    c = random.uniform(0.4,0.8)
                    integral_flux_spline.plot(num_points=1000,logx=True,\
                                              logy=True,show=False,\
                                              color='%f'%c,linewidth=1.0)
                    num_grb += 1
        else:
            continue
    logger.info('%i GRBs included in the plot.'%num_grb)
    if show:
        plt.show()

    def rvs_event_list(self, aeff, psf, modf, edisp, **kwargs):
        """Extract an event list for the full ROI.

        Arguments
        ---------
        aeff : :py:class:`ximpol.irf.arf.xEffectiveArea` object.
            The effective area to be used.

        psf : :py:class:`ximpol.irf.psf.xPointSpreadFunction` object.
            The PSF to be used.

        modf : :py:class:`ximpol.irf.mrf.xModulationFactor` object.
            The modulation factor to the used.

        edisp : :py:class:`ximpol.irf.rmf.xEnergyDispersion` object.
            The energy dispersion to be used.

        sampling_time : array
            The array to sample the source light curve.

        Warning
        -------
        The sampling_time should not be the same for all sources, and each
        source should be able to decide its own in a sensible way.
        (See issue #44.)
        """
        event_list = xMonteCarloEventList()
        for source in self.values():
            logger.info('Generating event list for source "%s"...' %\
                        source.name)
            event_list += source.rvs_event_list(aeff, psf, modf, edisp,
                                                **kwargs)
            event_list.apply_vignetting(aeff, self.ra, self.dec)
        event_list.sort()
        return event_list

 def __init__(self, psf_file_path):
     """Constructor.
     """
     logger.info("Reading PSF data from %s..." % psf_file_path)
     self.hdu_list = fits.open(psf_file_path)
     self.hdu_list.info()
     _data = self.hdu_list["PSF"].data
     W = _data["W"]
     sigma = _data["SIGMA"]
     N = _data["N"]
     r_c = _data["R_C"]
     eta = _data["ETA"]
     self.__params = (W, sigma, N, r_c, eta)
     # Tabulate the actual PSF values.
     _r = numpy.linspace(0, self.MAX_RADIUS, 250)
     _y = gauss_king(_r, *self.__params)
     fmt = dict(xname="r", xunits="arcsec", yname="PSF", yunits="sr$^{-1}$")
     xInterpolatedUnivariateSpline.__init__(self, _r, _y, k=2, **fmt)
     # Include the solid angle for the actual underlying random generator.
     _y *= 2 * numpy.pi * _r
     fmt = dict(rvname="r", rvunits="arcsec", pdfname="$2 \\pi r \\times$ PSF", pdfunits="")
     self.generator = xUnivariateGenerator(_r, _y, k=1, **fmt)
     # Finally, calculate the
     self.eef, self.hew = self.build_eef()
     logger.info(self)

def parse_blazar_list(PRIORITY_ONLY):
    """
    """
    logger.info("Parsing input file %s..." % BLAZAR_LIST_PATH)
    src_list = []
    input_file = open(BLAZAR_LIST_PATH)
    for i in range(6):
        input_file.next()
    for line in input_file:
        if line.startswith("-"):
            if PRIORITY_ONLY is False:
                for i in range(3):
                    line = input_file.next()
                PRIORITY_ONLY = True
            else:
                return src_list
        line = line.replace("BL Lac", "BL_Lac")
        name, line, notes = line[:12].strip(), line[12:95], line[95:]
        ra, dec, opt_class, sed_class, flux_max, flux_min, p_opt_max, p_opt_min = line.split()
        flux_max = float(flux_max)
        flux_min = float(flux_min)
        p_opt_max = float(p_opt_max)
        p_opt_min = float(p_opt_min)
        if p_opt_min < 0.5:
            p_opt_min = 0.5
        if p_opt_max > 0.5:
            src = {
                "name": name,
                "flux_min": flux_min,
                "flux_max": flux_max,
                "p_opt_max": p_opt_max,
                "p_opt_min": p_opt_min,
            }
            src_list.append(src)

def analyze():
    """
    """
    if os.path.exists(ANALYSIS_FILE_PATH):
        logger.info('%s exists, delete it if you want to recreate it.' %\
                    ANALYSIS_FILE_PATH)
        return
    logger.info('Opening output file %s...' % ANALYSIS_FILE_PATH)
    analysis_file = open(ANALYSIS_FILE_PATH, 'w')
    for i, (_min, _max) in enumerate(PHASE_BINNING):
        _mcube = xBinnedModulationCube(_mcube_file_path(i))
        _mcube.fit()
        _fit_results = _mcube.fit_results[-1]
        print _fit_results
        _phase = 0.5*(_min + _max)
        _phase_err = 0.5*(_max - _min)
        _pol_deg = _fit_results.polarization_degree
        _pol_deg_err = _fit_results.polarization_degree_error
        _pol_angle = _fit_results.phase
        _pol_angle_err = _fit_results.phase_error
        _data = (_phase, _phase_err, _pol_deg, _pol_deg_err, _pol_angle,
                 _pol_angle_err)
        _fmt = ('%.4e   ' * len(_data)).strip()
        _fmt = '%s\n' % _fmt
        _line = _fmt % _data
        analysis_file.write(_line)
    analysis_file.close()

def build_grb_fits_file(data,outfile):
    primary_hdu = xPrimaryHDU()
    grb_info_hdu = xBinTableGRBmain(data)
    hdu_list = fits.HDUList([primary_hdu, grb_info_hdu])
    hdu_list.info()
    logger.info('Writing GRB main infos table to %s...' % outfile)
    hdu_list.writeto(outfile, clobber=True)
    logger.info('Done.')

def parse(file_path, emin=1., emax=15., flux_scale=4.3576):
    """Parse the input file with the complete spectral and polarization model.
    """
    logger.info('Parsing input file %s...' % file_path)
    energy, flux, angle, degree = numpy.loadtxt(file_path, unpack=True)
    flux *= flux_scale
    _mask = (energy >= emin)*(energy <= emax)
    return energy[_mask], flux[_mask], angle[_mask], degree[_mask]

def xpsrccoords(source_name):
    """
    """
    logger.info('Querying CDS name resolver for "%s"...' % source_name)
    coords = SkyCoord.from_name(source_name)
    print(coords.icrs)
    print(coords.galactic)
    logger.info('Done, bye!')

 def __init__(self, mrf_file_path):
     """Constructor.
     """
     logger.info('Reading energy dispersion data from %s...' % mrf_file_path)
     self.hdu_list = fits.open(mrf_file_path)
     self.hdu_list.info()
     self.matrix = xEnergyDispersionMatrix(self.hdu_list['MATRIX'])
     self.ebounds = xEnergyDispersionBounds(self.hdu_list['EBOUNDS'])

def parse(file_path, emin=1., emax=10.):
    """Parse the input file with the complete spectral and polarization model.
    """
    logger.info('Parsing input file %s...' % file_path)
    energy, flux, degree, angle = numpy.loadtxt(file_path, unpack=True,
                                                        usecols = (0,1,5,6))
    angle += 90.
    _mask = (energy >= emin)*(energy <= emax)
    return energy[_mask], flux[_mask], degree[_mask], angle[_mask]

 def __init__(self, file_path):
     """Constructor.
     """
     assert(file_path.endswith('.fits'))
     logger.info('Opening input event file %s...' % file_path)
     self.hdu_list = fits.open(file_path)
     self.hdu_list.info()
     self.event_data = self.hdu_list['EVENTS'].data
     self.roi_table = self.build_roi_table()

 def fit_bin(self, i):
     """Fit the azimuthal distribution for the i-th energy slice.
     """
     hist = (self.phi_y[i], self.phi_binning, None)
     _fit_results = xAzimuthalResponseGenerator.fit_histogram(hist)
     _fit_results.set_polarization(self.modf(self.emean[i]))
     logger.info(_fit_results)
     self.fit_results.append(_fit_results)
     return _fit_results

def run():
    PIPELINE.xpobssim(configfile=CFG_FILE, duration=SIM_DURATION,
                      outfile=EVT_FILE_PATH)
    pha1_file_path = PIPELINE.xpbin(EVT_FILE_PATH, algorithm='PHA1')
    spec_fitter = PIPELINE.xpxspec(pha1_file_path)
    (index, index_err), (norm, norm_err) = spec_fitter.fit_parameters()
    logger.info('Fitted PL norm = %.4f +- %4f (input = %.4f)' %\
                (norm, norm_err, PL_NORM))
    logger.info('Fitted PL index = %.4f +- %4f (input = %.4f)' %\
                (index, index_err, PL_INDEX))