from matplotlib.pyplot import *
import glob
import os

import heasoftpy as hsp
import xspec

from astropy.table import Table
from astropy.time import Time
from astropy.io import fits
import astropy.units as u

Set up the NICER obsid directory

We'll assume the user has downloaded the NICER observation directory for OBSID 4142010107 in a directory called /Users/mcorcora/tmp/nicer

That is the directory /Users/mcorcora/tmp/nicer/4142010107 contains the xti, auxil and hk directories containing the standard archived data for the observation

We'll put the output from nicerl2 in an output directory separate from the input directory

cwd = os.getcwd()
os.chdir(nicerdatadir)

nicerdatadir = os.path.join(os.environ['HOME'],'tmp/nicer')
nicerobsID = '4020180460'
obsdir = os.path.join(nicerdatadir, nicerobsID)

# place cleaned output in a separate directory
outdir =  os.path.join(os.environ['HOME'],'tmp/nicer/','nicerl2_output/'+nicerobsID+'_out')
# if outdir doesn't exist, create it
if not os.path.exists(outdir):
    os.makedirs(outdir)
    print(f'Created {outdir}')
# copy the mkf file from the input directory to the outdir

mkf = os.path.join(obsdir,'auxil',f'ni{nicerobsID}.mkf')
if os.path.exists(mkf):
    # see if mkf is gzipped    
    cmd = f'cp {mkf} {outdir}/.'
    stat=os.system(cmd)
    mkf = os.path.join(outdir, os.path.split(mkf)[1])
    print(f'Setting mkf file to {mkf}')
elif os.path.exists(mkf+'.gz'):
    #try to copy gzipped mkf     
    cmd = f'cp {mkf}.gz {outdir}/.'
    print(cmd)
    os.system(cmd)
    mkf = os.path.join(outdir, os.path.split(mkf)[1])
    print(f'Setting mkf file to {mkf}')
    #cmd = f'gunzip -f {mkf}.gz'
    #print(cmd)
    #stat=os.system(cmd)
cmd = f'chmod u+w {mkf}*'
print(cmd)
stat = os.system(cmd)

cp /Users/mcorcora/tmp/nicer/4020180460/auxil/ni4020180460.mkf.gz /Users/mcorcora/tmp/nicer/nicerl2_output/4020180460_out/.
Setting mkf file to /Users/mcorcora/tmp/nicer/nicerl2_output/4020180460_out/ni4020180460.mkf
chmod u+w /Users/mcorcora/tmp/nicer/nicerl2_output/4020180460_out/ni4020180460.mkf*

Create the nicerl2 task

tstart = Time.now()
print(f'Start at: {tstart.iso[:19]}')
nicerl2 = hsp.HSPTask('nicerl2')

nicerl2.clobber="yes"
#nicerl2.cldir=outdir
#nicerl2.mkffile=mkf
# add the KP values to the mkf file during nicerl2 processing
nicerl2.geomag_path="https://heasarc.gsfc.nasa.gov/FTP/caldb/data/gen/pcf/geomag/"
nicerl2.geomag_columns="kp_noaa.fits(KP)"

resl2 = nicerl2(indir=nicerobsID, noprompt=True, cldir=outdir, mkfile=mkf)

tend = Time.now()
print(f'End at: {tend.iso[:19]}')
print(f'nicerl2 took: {(tend.mjd-tstart.mjd)*86400} seconds')

if resl2.returncode != 0:
    print('\n')
    for o in resl2.output[:]:
        print(o)

Start at: 2021-12-18 18:08:30
End at: 2021-12-18 18:11:36
nicerl2 took: 185.77480078674853 seconds

Extract products from cleaned events file

clevt = f'{outdir}/ni{nicerobsID}_0mpu7_cl.evt'
phafile = f'{outdir}/ni{nicerobsID}_0mpu7_cl.pha'
lcfile = f'{outdir}/ni{nicerobsID}_0mpu7_cl.lc'
res = hsp.extractor(filename=clevt, phafile=phafile, clobber='yes', binlc=10.0,fitsbinlc=lcfile, 
                    eventsout='NONE', imgfile='NONE', regionfile='NONE', timefile='NONE', tcol='TIME',
                   ecol='PI', xcolf='RAWX', xcolh='RAWX',ycolf='RAWY', ycolh='RAWY', )

Analyzing NICER spectra

# get the on-axis rmf
res = hsp.quzcif(mission='nicer', instrument='xti',detector='-',
             filter='-', date='-', time='-',expr='-',codename='MATRIX')
rmf = [x.split()[0] for x in res.output if 'nixtiref'  in x][0]

# get the on-axis arf
res = hsp.quzcif(mission='nicer', instrument='xti',detector='-',
             filter='-', date='-', time='-',expr='-',codename='SPECRESP')
arf = [x.split()[0] for x in res.output if 'nixtiaveonaxis'  in x][0]

xspec.AllData.clear()
spec = xspec.Spectrum(phafile)
spec.response = rmf
spec.response.arf = arf
spec.ignore('0.0-0.3, 10.0-**')

1 spectrum  in use
 
Spectral Data File: /Users/mcorcora/tmp/nicer/nicerl2_output/4020180460_out/ni4020180460_0mpu7_cl.pha  Spectrum 1
Net count rate (cts/s) for Spectrum:1  7.015e+01 +/- 2.220e-01
 Assigned to Data Group 1 and Plot Group 1
  Noticed Channels:  1-1501
  Telescope: NICER Instrument: XTI  Channel Type: PI
  Exposure Time: 1423 sec
 Using fit statistic: chi
 No response loaded.

***Warning!  One or more spectra are missing responses,
               and are not suitable for fit.
Response successfully loaded.

Fit statistic  : Chi-Squared                99491.94     using 1501 bins.

***Warning: Chi-square may not be valid due to bins with zero variance
            in spectrum number: 1

Test statistic : Chi-Squared                99491.94     using 1501 bins.

***Warning: Chi-square may not be valid due to bins with zero variance
            in spectrum number(s): 1 

 Null hypothesis probability of 0.00e+00 with 1496 degrees of freedom
 Current data and model not fit yet.
Arf successfully loaded.

Fit statistic  : Chi-Squared                27610.54     using 1501 bins.

***Warning: Chi-square may not be valid due to bins with zero variance
            in spectrum number: 1

Test statistic : Chi-Squared                27610.54     using 1501 bins.

***Warning: Chi-square may not be valid due to bins with zero variance
            in spectrum number(s): 1 

 Null hypothesis probability of 0.00e+00 with 1496 degrees of freedom
 Current data and model not fit yet.
    30 channels (1-30) ignored in spectrum #     1
   502 channels (1000-1501) ignored in spectrum #     1

Fit statistic  : Chi-Squared                 3127.73     using 969 bins.

Test statistic : Chi-Squared                 3127.73     using 969 bins.
 Null hypothesis probability of 2.72e-226 with 964 degrees of freedom
 Current data and model not fit yet.

model = xspec.Model('wabs*bknpow')
xspec.Fit.perform()

========================================================================
Model wabs<1>*bknpower<2> Source No.: 1   Active/On
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   wabs       nH         10^22    1.00000      +/-  0.0          
   2    2   bknpower   PhoIndx1            1.00000      +/-  0.0          
   3    2   bknpower   BreakE     keV      5.00000      +/-  0.0          
   4    2   bknpower   PhoIndx2            2.00000      +/-  0.0          
   5    2   bknpower   norm                1.00000      +/-  0.0          
________________________________________________________________________


Fit statistic  : Chi-Squared              1.146821e+08     using 969 bins.

Test statistic : Chi-Squared              1.146821e+08     using 969 bins.
 Null hypothesis probability of 0.000000e+00 with 964 degrees of freedom
 Current data and model not fit yet.
                                   Parameters
Chi-Squared  |beta|/N    Lvl          1:nH    2:PhoIndx1      3:BreakE    4:PhoIndx2        5:norm
44865        1293.48      -1      0.161880       1.12883       4.91288     -0.625071     0.0104577
40493        253023       -1     0.0528104       1.00933       4.75666      -2.42040    0.00793820
36526.8      368439       -2     0.0118557      0.814400       5.18901      -2.39618    0.00758772
18922.4      379942       -3    0.00102210       1.23639       5.11146      -2.49539     0.0113204
3632.22      319860       -4   0.000353044       1.65880       4.69324      -2.49812     0.0169442
2744.36      84924.3      -5    0.00587835       1.77760       4.38823      -2.42043     0.0188889
2720.56      4508.26      -6    0.00855859       1.81690       4.21915      -2.26013     0.0190699
2665.99      6978.14      -7    0.00979611       1.83766       4.09490      -2.06221     0.0191276
2604.57      5457.51      -8     0.0106566       1.85162       3.99463      -1.82726     0.0191670
2584.8       1246.82      -9     0.0113103       1.86224       3.90286      -1.66137     0.0191955
Number of trials exceeded: continue fitting? 2573.02      159.566     -10     0.0119420       1.87245       3.82317      -1.55576     0.0192219
2563.26      77.9134     -11     0.0125737       1.88267       3.75182      -1.46534     0.0192475
2554.59      60.9        -12     0.0131747       1.89243       3.68552      -1.38275     0.0192712
2546.53      53.7017     -13     0.0137550       1.90186       3.62221      -1.30494     0.0192934
2539.33      49.4844     -14     0.0143106       1.91091       3.56206      -1.23186     0.0193141
2533.6       43.2747     -15     0.0148108       1.91908       3.50810      -1.16682     0.0193323
2529.44      34.4934     -16     0.0152374       1.92607       3.46215      -1.11173     0.0193473
2526.58      26.0789     -17     0.0155914       1.93188       3.42411      -1.06632     0.0193596
2524.67      19.0932     -18     0.0158803       1.93664       3.39316      -1.02951     0.0193693
2523.4       13.9864     -19     0.0161162       1.94052       3.36797     -0.999669     0.0193772
Number of trials exceeded: continue fitting? 2522.5       10.4501     -20     0.0163141       1.94378       3.34702     -0.974993     0.0193836
2521.84      8.2021      -21     0.0164859       1.94661       3.32906     -0.953978     0.0193892
2521.31      6.8267      -22     0.0166395       1.94914       3.31315     -0.935491     0.0193940
2520.86      5.93054     -23     0.0167813       1.95148       3.29863     -0.918728     0.0193985
2520.46      5.38035     -24     0.0169151       1.95368       3.28504     -0.903129     0.0194027
2520.09      5.02325     -25     0.0170436       1.95580       3.27208     -0.888316     0.0194067
2519.74      4.78649     -26     0.0171689       1.95786       3.25953     -0.874038     0.0194105
2519.41      4.64973     -27     0.0172916       1.95988       3.24727     -0.860140     0.0194142
2519.08      4.53807     -28     0.0174126       1.96188       3.23523     -0.846534     0.0194178
2518.76      4.41661     -29     0.0175321       1.96385       3.22339     -0.833184     0.0194214
Number of trials exceeded: continue fitting? 2518.45      4.31972     -30     0.0176499       1.96579       3.21175     -0.820096     0.0194249
2518.16      4.25059     -30     0.0177651       1.96770       3.20037     -0.807317     0.0194283
2517.88      4.10151     -30     0.0178771       1.96955       3.18932     -0.794926     0.0194315
2517.62      3.9538      -30     0.0179848       1.97133       3.17868     -0.783020     0.0194347
2517.38      3.75082     -30     0.0180872       1.97302       3.16858     -0.771710     0.0194376
2517.18      3.52919     -30     0.0181832       1.97461       3.15909     -0.761098     0.0194403
2517         3.23613     -30     0.0182723       1.97609       3.15030     -0.751261     0.0194428
2516.85      3.0112      -30     0.0183536       1.97744       3.14223     -0.742242     0.0194451
2516.73      2.67473     -30     0.0184275       1.97866       3.13491     -0.734059     0.0194472
2516.63      2.40285     -30     0.0184939       1.97976       3.12833     -0.726699     0.0194490
Number of trials exceeded: continue fitting? 2516.56      2.14897     -30     0.0185531       1.98075       3.12246     -0.720134     0.0194507
2516.5       1.8837      -30     0.0186054       1.98162       3.11726     -0.714321     0.0194521
2516.45      1.6569      -30     0.0186514       1.98238       3.11269     -0.709208     0.0194534
2516.41      1.4482      -30     0.0186917       1.98305       3.10869     -0.704734     0.0194545
2516.39      1.25523     -30     0.0187266       1.98363       3.10520     -0.700838     0.0194554
2516.37      1.08172     -30     0.0187570       1.98414       3.10218     -0.697461     0.0194562
2516.35      0.933294    -30     0.0187832       1.98457       3.09957     -0.694544     0.0194569
2516.34      0.808599    -30     0.0188057       1.98495       3.09733     -0.692033     0.0194575
2516.33      0.672596    -30     0.0188252       1.98527       3.09540     -0.689878     0.0194581
============================================================
 Variances and Principal Axes
                 1        2        3        4        5  
 4.0877E-09| -0.1069   0.0011  -0.0020  -0.0012   0.9943  
 2.7299E-07|  0.9902  -0.0890  -0.0143  -0.0088   0.1065  
 5.3004E-05|  0.0799   0.7832   0.5193   0.3323   0.0092  
 3.2857E-04|  0.0416   0.6122  -0.6049  -0.5075   0.0020  
 1.9478E-03|  0.0039   0.0624  -0.6035   0.7949   0.0001  
------------------------------------------------------------

=========================================================

%matplotlib notebook
xspec.Plot.device='/null'
xspec.Plot.xAxis='keV'
xspec.Plot('lda')
cr=xspec.Plot.y()
crerr = xspec.Plot.yErr()
en = xspec.Plot.x()
enwid = xspec.Plot.xErr()
mop = xspec.Plot.model()
target = fits.open(spec.fileName)[1].header['OBJECT']


fig = figure(figsize=[8,6])
ylabel('Cts/s/keV', fontsize=12)
xlabel('Energy (keV)', fontsize=12)
title('Target = '+target+' OBSID = '+nicerobsID+' wabs*bknpow', fontsize=12)
yscale('log')
xscale('log')
errorbar(en, cr, xerr=enwid, yerr=crerr, fmt='k.', alpha=0.2)
plot(en, mop,'r-')

[<matplotlib.lines.Line2D at 0x1762fe790>]

target = fits.open(spec.fileName)[1].header['OBJECT']
target

'PSR_B0833-45'

Plot the lightcurve

lctab = Table.read(lcfile,hdu='RATE')
gtitab = Table.read(lcfile, hdu='GTI')
gtitab['START']=gtitab['START']-lctab.meta['TSTART']
gtitab['STOP']=gtitab['STOP']-lctab.meta['TSTART']

gtitab

Table length=5

START	STOP
s	s
float64	float64
0.0	1.0
3.0	6.0
8.0	607.0
77995.0	77998.0
78000.0	78817.0

# remove rows with no values from gtitab

row2remove=[]
for j in enumerate(gtitab):
    i=j[0]
    tsel = (lctab['TIME']>=gtitab[i]['START']) & (lctab['TIME']<=gtitab[i]['STOP'])
    if len(lctab[tsel]) < 1:
        row2remove.append(i)
gtitab.remove_rows(row2remove)
numgtis = len(gtitab)

%matplotlib notebook
fig, ax =subplots(1,numgtis,figsize=[10,3])
for i in range(numgtis):
    tsel = (lctab['TIME']>gtitab[i]['START']) & (lctab['TIME']<gtitab[i]['STOP'])
    t = lctab[tsel]['TIME']
    print(i, len(lctab[tsel]))
    r = lctab[tsel]['RATE']
    re = lctab[tsel]['ERROR']
    ax[i].set_ylabel('Cts/s', fontsize=12)
    ax[i].set_xlabel('Time (s)', fontsize=12)
    ax[i].set_yscale('log')
    ax[i].set_ylim(40, 500)
    ax[i].errorbar(t, r,yerr=re, fmt='k.')
    tight_layout()

0 59
1 80