Analyzing Pulsars with the Fermi-LAT¶
Similar to other LAT sources, pulsars display both temporal and spectral signatures in LAT data. However, the periodicity for gamma-ray pulsars is much faster than the frequency at which photons are detected by the LAT.
Indeed, the Vela pulsar, the brightest persistent source in the LAT sky, produces (on average) one LAT event every 2 minutes. During that time, the pulsar has rotated 1766 times...emitting a radio pulse each time. These radio pulses provide a signal that can be easily used to characterize the pulsar.
While a number of pulsars have been discovered in the LAT gamma-ray data, the analysis required for that discovery is complex and requires significant computing power. Here we will discuss what can be done with pulsars that have already been detected and characterized for timing.
Once the timing solution for a pulsar has been found, it opens the door to analysis in the gamma rays. In order to leverage this knowledge, we first need a timing model that is valid during at least part of the Fermi mission. The easiest place to look is the most recent Pulsar Catalog, which includes all the pulsar timing files that were used in the process of developing the catalog.
Overview of steps¶
1) Perform Binned Likelihood analysis to acquire a fitted XML model file. This will take a while, depending on how much data you are analyzing.
2) Apply probabilities to the events and remove those with low probabilities. You can also weight the events by the probabiity.
3) Apply the timing model to the data and look for pulsations
Products already in hand¶
You should already have the following products for PSR J0248+86021, which is our source of interest for this analysis.
- Pulsar ephemeris from the Second Pulsar Catalog
- Filtered and GTI-adjusted event file for a 15-degree ROI centered at 42.0, 60.4
- Spacecraft file covering the entire validity range
Starting your analysis¶
As a rule, the E-dot (spindown power) for a pulsar needs to be greater than about 10^33 for it to have a chance of being detected in the gamma rays. Below that threshold, it is speculated, there is not enough energy in the system to spawn the processes that generate gamma rays.
The first step of searching in gamma rays for a pulsar with a known radio ephemeris is to determine if that source is even a pulsed gamma-ray source. To determine that, you should extract the events associated with that pulsar...
Initial Spectral Analysis¶
Many energetic young radio pulsars are faint in gamma rays, likely due to geometry...the part of the magnetosphere that produces radio emission only partially overlaps with the part that produces gamma rays. The same is true for bright young gamma-ray pulsars...nearly all the new pulsars discovered in the gamma rays are completely radio quiet. Only 4 of them have measurable radio emission, and one is the faintest radio pulsar ever detected.
However, this is not the case for recycled pulsars. In pulsars with millisecond periods, the two emission zones appear to have much more overlap. To date, only one gamma-ray MSP has been discovered that is completely undetectable in the radio.
Because pulsars are often faint, they can require careful treatment for gamma-ray spectral analysis.
The spectral fit for a faint source will be significantly affected if nearby, brighter sources are not properly modeled. This means you will need to free more parameters in your model, making the spectral fit take longer.
Let's first perform a binned analysis to characterize the sources in our ROI.
import gt_apps
dir(gt_apps)
from gt_apps import *
Create a binned counts cube from your data¶
This will generate the data file that is used in the Binned Likelihood analysis.
%system punlearn gtbin
counts_map.pars()
counts_map['evfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime.fits'
counts_map['scfile'] = 'NONE'
counts_map['outfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_ccube.fits'
counts_map['algorithm'] = 'CCUBE'
counts_map['xref'] = 42.0
counts_map['yref'] = 60.4
counts_map['coordsys'] = 'CEL'
counts_map['nxpix'] = 200
counts_map['nypix'] = 200
counts_map['binsz'] = 0.1
counts_map['axisrot'] = 0.0
counts_map['proj'] = 'AIT'
counts_map['emin'] = 300
counts_map['emax'] = 50000
counts_map['enumbins'] = 30
counts_map['ebinalg'] = 'LOG'
counts_map['rafield'] = 'RA'
counts_map['decfield'] = 'DEC'
counts_map.run()
Create a livetime cube (aka exposure cube)¶
Remember, this is where we should apply the zmax cut, since it was not accounted for in the maketime step.
expCube.pars()
expCube['evfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime.fits'
expCube['scfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/spacecraft.fits'
expCube['outfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_ltcube.fits'
expCube['dcostheta'] = 0.025
expCube['binsz'] = 1.0
expCube['zmax'] = 90
expCube.run()
What makes up an IRF?¶
Calculating the exposure cube is the first time you will need to use the Instrument Response Functions (IRFs). You can find details of the full suite of IRFs at the LAT Performance Page
If you install the Fermi Science Tools, the IRFs are also available in the following directory:
$FERMI_DIR/refdata/fermi/caldb/CALDB/data/glast/lat/bcf
from IPython.display import Image,HTML
HTML("<iframe src='http://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm' width='1000' height='500'></iframe>")
from os import environ
environ['CALDB']
ls /usr/local/ScienceTools/refdata/fermi/caldb/CALDB/data/glast/lat/bcf
ls '/usr/local/ScienceTools/refdata/fermi/caldb/CALDB/data/glast/lat/bcf/ea'
import pyfits
aeff_P8R2_CLEAN_V6_FB_HDUs = pyfits.open('/usr/local/ScienceTools/refdata/fermi/caldb/CALDB/data/glast/lat/bcf/ea/aeff_P7CLEAN_V6_back.fits')
print aeff_P8R2_CLEAN_V6_FB_HDUs.info()
Calculate an exposure map cube¶
You need to understand the exposure even for areas outside your ROI, so you should increase the size of the exposure map well above the size of your counts cube. Often, it's easier to just calculate the exposure map for the entire sky.
You also want to be sure that the geometry for your exposure cube exactly matches that of your counts cube. Otherwise, the mismatch will cause an analysis error.
gtexpcube2.pars()
gtexpcube2['infile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_ltcube.fits'
gtexpcube2['cmap'] = 'none'
gtexpcube2['outfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_expmap.fits'
gtexpcube2['irfs'] = 'P8R2_SOURCE_V6'
gtexpcube2['xref'] = 42.0
gtexpcube2['yref'] = 60.4
gtexpcube2['coordsys'] = 'CEL'
gtexpcube2['nxpix'] = 500
gtexpcube2['nypix'] = 500
gtexpcube2['binsz'] = 0.1
gtexpcube2['axisrot'] = 0.0
gtexpcube2['proj'] = 'AIT'
gtexpcube2['emin'] = 300
gtexpcube2['emax'] = 50000
gtexpcube2['enumbins'] = 30
gtexpcube2.run()
Generate an XML file for your region
The next step requires an XML model. Your model should reflect all the sources in your region of interest, plus sources that lie in a 10-degree ring around the outside of your ROI. Because sources outside your ROI can contribute photons to your data, they must be considered, even though all their parameters should be fixed.
You will need to download the Galactic and Isotropic diffuse models for Pass 8 and be certain they are included in the model file for your analysis. You can find the template files on the FSSC Background Models page, or in the Fermi Science Tools installation under $FERMI_DIR/refdata/fermi/galdiffuse/
ls $FERMI_DIR/refdata/fermi/galdiffuse/
Once you have the templates downloaded, be sure that your XML model reflects the correct filenames. Here are my entries for the galactic and isotropic diffuse models:
Image(filename='images/P8_J0248_diffuse_xml.png')
Calculate source maps for each source in your XML file¶
This creates a map for each source which convolves the IRFs with the exposure and the binned geometry. The resulting source maps file contains a plane for each source which is a representation of the data. These are then fitted to your model file in the Likelihood fit.
srcMaps.pars()
srcMaps['scfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/spacecraft.fits'
srcMaps['expcube'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_ltcube.fits'
srcMaps['cmap'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_ccube.fits'
srcMaps['bexpmap'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_expmap.fits'
srcMaps['srcmdl'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_pass8.xml'
srcMaps['outfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_srcmaps.fits'
srcMaps['irfs'] = 'P8R2_SOURCE_V6'
srcMaps.run()
Check that you have all the input files needed for the Binned Analysis
You will need the following files:
- Counts cube
- Exposure (Livetime) Cube
- Binned Exposure Map
- XML model file
- Source Maps file
%system ls /Users/eferrara/FSSC/VM/shared/psr_data
Perform a Binned Likelihood fit¶
import BinnedAnalysis
dir(BinnedAnalysis)
from BinnedAnalysis import *
help(BinnedObs)
obs = BinnedObs(binnedExpMap='/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_expmap.fits',
expCube='/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_ltcube.fits',
srcMaps='/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_srcmaps.fits',
irfs='P8R2_SOURCE_V6')
analysis = BinnedAnalysis(obs,srcModel='/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_pass8.xml', optimizer='MINUIT')
likeObj = pyLike.Minuit(analysis.logLike)
analysis.tol = 0.01
analysis.tol
analysis.fit(covar=True,optObject=likeObj)
Check for convergence - Write XML¶
analysis.writeXml('/Users/eferrara/FSSC/VM/shared/psr_data/P8_PSRJ0248_output_model.xml')
print likeObj.getRetCode()
Fitter convergence gives a return code of zero. With any other value, you may need to try again, using modified inputs. If you don't get convergence easily, you can try:
- Change the optimizer. Sometimes Minuit works, other times NewMinuit. Try both.
- Reduce the tolerance. Increase it iteratively.
- Reduce the number of free parameters.
- Make sure no parameters in your model have moved to the limit values.
analysis.Ts('PSRJ0248+6021')
analysis.model['PSRJ0248+6021']
analysis.NpredValue('PSRJ0248+6021')
print "PSRJ0248+6021 Integral Flux = {} +- {}".format(analysis.flux('PSRJ0248+6021'), analysis.fluxError('PSRJ0248+6021'))
There is a scond pulsar in the ROI. How does it compare?
analysis.Ts('PSRJ0205+6449')
analysis.model['PSRJ0205+6449']
analysis.NpredValue('PSRJ0205+6449')
print "PSRJ0205+6449 Integral Flux = {} +- {}".format(analysis.flux('PSRJ0205+6449'), analysis.fluxError('PSRJ0205+6449'))
Iterate as needed¶
Here I loosened my tolerance to allow the process to run more quickly. However, that's not appropriate for a scientific result. The default tolerance is 0.001, but you may want to go tighter.
When you iterate, you should use the output model from this fit so that you are closer to the (hopefully) correct result. However, you will want to move any parameters that have hit their limits away from the limit value. Unfortunately, there is not an easy way to do this within the pyLikelihood framework. This is an excellent reason to use quickLike instead, which is described in the Analysis Threads.
Event Probabilities¶
Once the initial spectral analysis is complete, it is simple to run gtsrcprob. The tool generates a new file with additional columns, one for each source you specify in the sourcelist.txt file. It then determines, for each event, how likely it is for that event to have originated from that source, and records the fractional probability in each column. Each event should have probabilities that add up to one.
Next you can use ftselect to filter out events below a particular probability threshold of your choosing. (Here I used 20%.) While this reduces the statistics for your search, it greatly improves the S/N ratio and greatly improves the chance of detecting significant pulsations.
BE AWARE!! Changing the number of events using a method where the GTIs cannot be corrected, means that the resulting file is not appropriate for measuring the flux or spectral parameters of your source. We are searching for a temporal signature, so this is not a concern for us.
import gt_apps
from gt_apps import *
filter['infile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_filtered.fits'
filter['outfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_filtered_3deg.fits'
filter['ra'] = 42.0
filter['dec'] = 60.4
filter['rad'] = 3
filter['tmin'] = 'INDEF'
filter['tmax'] = 'INDEF'
filter['emin'] = 300
filter['emax'] = 50000
filter['zmax'] = 100
filter['evclass'] = 128
filter['evtype'] = 3
filter.run()
maketime['evfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_filtered_3deg.fits'
maketime['outfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime_3deg.fits'
maketime['scfile'] = '/Users/eferrara/FSSC/VM/shared/psr_data/spacecraft.fits'
maketime['filter'] = 'DATA_QUAL>0 && LAT_CONFIG==1'
maketime['apply_filter'] = 'yes'
maketime['roicut'] = 'no'
maketime.run()
import pyfits
original_data = pyfits.open('/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime_3deg.fits')
original_data.info()
%system cp /Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime_3deg.fits /Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime_3deg.fits.bak
%system cat /Users/eferrara/FSSC/VM/shared/psr_data/sourcelist.txt
Calculate diffuse responses¶
Since gtsrcprob treats each event separately, it uses unbinned likelihood. As a result, each event will need an associated diffuse response for every diffuse component in your model. (This includes extended sources.)
To calculate the diffuse responses, run gtdiffrsp. This will add columns to your events file, so make a copy of the file first, if you're concerned.
%system gtdiffrsp \
evfile=/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime_3deg.fits \
scfile=/Users/eferrara/FSSC/VM/shared/psr_data/spacecraft.fits \
srcmdl=/Users/eferrara/FSSC/VM/shared/psr_data/P8_PSRJ0248_output_model.xml \
irfs=P8R2_SOURCE_V6
%system gtsrcprob \
srclist=/Users/eferrara/FSSC/VM/shared/psr_data/sourcelist.txt \
evfile=/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_maketime_3deg.fits \
scfile=/Users/eferrara/FSSC/VM/shared/psr_data/spacecraft.fits \
srcmdl=/Users/eferrara/FSSC/VM/shared/psr_data/P8_PSRJ0248_output_model.xml \
irfs=P8R2_SOURCE_V6 \
outfile=/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_prob_3deg.fits
Here I have edited the output FITS file (J0248_3deg_prob.fits) because ftselect does not recognize column names containing a "+" symbol. Frustrating!
Fixing this would have required changing the XML model and re-running the analysis.
%system ftselect \
/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_prob_3deg.fits \
/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_downselected_3deg.fits \
'PSRJ0248 > 0.2'
downselected_data = pyfits.open('/Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_downselected_3deg.fits')
downselected_data.info()
Apply the Timing Model¶
Once you have purified your event list, you can apply the radio timing model to see if there are significant gamma-ray pulsations. This is best performed in tempo2, using the user-contributed tempo2 plugin.
The tempo2 plugin does not provide any pulse search capability. It simply applies the radio timing model to the gamma-ray data, and then plots the results, including the H-test result. The tool also assigns phase values to the events in the LAT data file.
Let's investigate what's available with tempo2:
%system tempo2 -h
%system tempo2 -gr fermi -h
The tempo2 command is:
%system tempo2 -gr fermi -ft1 /Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248_downselected_3deg.fits -ft2 /Users/eferrara/FSSC/VM/shared/psr_data/spacecraft.fits -f /Users/eferrara/FSSC/VM/shared/psr_data/PSRJ0248+6021_2PC.par -phase
Image(filename='images/J0248_tempo2_out.png')
As you can see, the H-test for this source is about 350. Anything over 25 is considered a detection.
Image(filename='images/PSRJ0248+6021_2PC_LC.png')
%system plist gtpphase