paulray / nicersoft Goto Github PK

Analysis software for the NICER mission

License: MIT License

Python 34.53% Jupyter Notebook 65.47%

astronomy python x-ray pulsars

nicersoft's Introduction

NICERsoft

User-contributed analysis tools for the NICER mission. This is separate from the NICER-team-developed tools that are part of the HEASoft distribution. Nothing here should be construed as formally recommended by the NICER instrument team.

Please contribute tools you develop!

Disclaimer

All of this code is user contributed! Use at your own risk! No warranty is expressed or implied. There may be bugs! Do not use any of this code without understanding what it is doing.

Documentation

See the Wiki for some documentation: https://github.com/paulray/NICERsoft/wiki

Usage

To use these scripts add <basedir>/NICERsoft/scripts to your PATH and add <basedir>/NICERsoft to your PYTHONPATH, where <basedir> is wherever you cloned NICERsoft.

NOTES:

Executable scripts go in the scripts subdirectory
Importable python modules go in the nicer subdirectory
Scripts in any language (Python, PERL, bash, tcsh, etc...) are OK to contribute!

PREREQUISITES:

pip-installable packages:
- numpy, matplotlib, scipy, astropy, pint-pulsar, loguru, tqdm
Heasoft (with NICER Data Analysis Software)
Scripts using PyXspec require Heasoft to be compiled from source.

nicersoft's People

Contributors

Stargazers

Watchers

Forkers

sguillot peterbult kerrm scottransom jcolen juliadeneva astrojuan zaynabghazi masonng-astro ttgsgmafe amrutajaodand champagne-cmd ktzhao krizmo1992 astrobuff xstarkit lmauviard liangbingkaer harej10

nicersoft's Issues

nitemplate.py is not python3 compatible

This is due to the several python2 style print lines.

Make BKG_PLOT faster

@paulray following your suggestion in #33
Let me know how I can help.
Looking at bkg_plots.py, most plots already come from the MKF file, except the overshootrate and undershootrate plots., right?

CFITSIO error in nifpmsel

When running psrpipe.py the nifpmsel step fails.

Task nifpmsel 1.8 terminating with status -1
Traceback (most recent call last):
  File "/scratch/software/nicer_20230814/NICERsoft/scripts/psrpipe.py", line 653, in <module>
    runcmd(cmd)
  File "/scratch/software/nicer_20230814/NICERsoft/scripts/psrpipe.py", line 256, in runcmd
    check_call(cmd, env=os.environ)
  File "/usr/lib/python3.11/subprocess.py", line 413, in check_call
    raise CalledProcessError(retcode, cmd)
subprocess.CalledProcessError: Command '['nifpmsel', 'infile=2070020493_pipe/intermediate.evt', 'outfile=2070020493_pipe/cleanfilt.evt', 'detlist=launch', 'mkfile=2070020493/auxil/ni2070020493.mkf', 'clobber=yes', 'history=yes']' returned non-zero exit status 255.

Rerunning this command by hand shows the following error:

nifpmsel infile=2070020493_pipe/intermediate.evt outfile=2070020493_pipe/cleanfilt.evt detlist=launch mkfile=2070020493/auxil/ni2070020493.mkf clobber=yes history=yes
nifpmsel 1.8
  FPM List: launch
  Recalculating GTIs
ERROR: (internal error)  GTI is missing extension name at /scratch/software/nicer_20230814/heasoft-6.32/x86_64-pc-linux-gnu-libc2.36/bin/nifpmsel line 966.
Dumping CFITSIO error stack:
--------------------------------------------------
Extension doesn't start with SIMPLE or XTENSION keyword. (ffrhdu)
END
Error while closing HDU number 0 (ffchdu).
--------------------------------------------------
CFITSIO error stack dump complete.
Task nifpmsel 1.8 terminating with status -1

--ultraclean in psrpipe.py

The --ultraclean option mentioned in the description of psrpipe.py is actually undefined.
I guess, it is supposed to say --dark instead of --ultraclean

Figtext for many GTIs in sci_plots (minor issue)

When the event file has a lot of GTIs, the display command
plt.figtext(str(gtitable['START'][:]))
of the astropy.Table automatically truncates the number of rows in the table, but there are still too many, and the list of GTIs overlaps with the title and goes out of the page.

That's a minor issue, but I can't find a way to truncate the display even more. There doesn't seem to be any way to tell 'print' or 'str()' to shorten the output of the astropy.Table

Any ideas?

Make NICERsoft a pip-installable package

Should make this so that it can be installed by pip and thus used in virtualenvs and conda without setting PYTHONPATH

Improving ENG plot

Some suggestions to improve the ENG plots:

rename "Reset rate" to "Undershoot" for clarity

Python3 compatibility

We should make sure this script runs under Python3 as well as Python2.7
Here are the best practices for doing that:
https://github.com/nanograv/PINT/wiki/Python-Compatibility

Updates the output figures (background, etc...)

Considering what has been learned about the mission in the past 5 years, the output diagnostics plots could be updated with more relevant parameters. We should think about:

which plots can be removed ?
which plots could be added (e.g., 3C50 parameters).

Add option to perform barycentric correction

Following en email discussion, I am adding this suggestion to the list of issues:
Having an option to barycentre the photons during the psrpipe processing would be useful.

This requires being able to:

specify the precise position to use
ephemeris to use

Improvements to cr_cut.py (C. Markwardt's suggestion)

The steps 3 to 6 of the script cr_cut.py could be replaced by a single command maketime:

maketime merged_detid14.lc merged_detid14.gti "RATE < 1.0" name=NAME value=VALUE compact=NO prefr=0.5 postfr=0.5 premax=4.0 postmax=4.0 time=TIME

Interpolation problem between orbit files from merge.py

Following an email discussion with @paulray:

When calling pint.photonphase from merge.py, for a merged event file with a list of orbit files, there seems to be an interpolation problem with the orbit files. This gives different results than when the PULSE_PHASES are calculated on each individual ObsID event files.
This only happens when using merge.py independently on events files that do not have the PULSE_PHASE column. This does not happen when the options '--par parfile.par' and '--merge' are used in psrpipe.py.

Replace basemap

matplotlib.basemap has been deprecated. We should switch to cartopy or some other package for making the maps.

Improving the BKG plots

Some discussed suggestions to improve the BKG plots:

move the COR limit (red dashed line) to 1.5 ?
Show COR plot as y-log scale for clearer plot
make COR markers smaller to see the values better?
show minor ticks grid lines ?

nicerfits2presto.py error when --dec DEC is a negative value

When I run nicerfits2presto.py for a target with a negative declination, I get an error:
> nicerfits2presto.py --dt 0.125 --ra 17:32:03.3 --dec -34:45:18 merged_cut_barycorr.evt usage: nicerfits2presto.py [-h] [--dt DT] [--ra RA] [--dec DEC] [--observer OBSERVER] evfile nicerfits2presto.py: error: argument --dec: expected one argument

This does not happen if the --dec is positive

psrpipe.py wants unzipped orb and att files

By defaults, files from the HEASARC archives are gzipped. Running nicerl2 creates unzipped *evt and *mkf files, but the *orb and *att files remained gzipped. This throws an list index out of range error at line 373 of psrpipe.py, which may be cryptic to new users.

  File "/Users/sebastien/Softwares/NICERsoft/scripts/psrpipe.py", line 373, in <module>
    orbfile = glob(path.join(obsdir, "auxil/ni*.orb"))[0]
IndexError: list index out of range

I can think of a few solutions:

Catching the error with a warning that *orb and *att files must be unzipped.
Replacing glob(path.join(obsdir, "auxil/ni*.orb"))[0] by glob(path.join(obsdir, "auxil/ni*.orb*"))[0], which would work with the downloaded files by default (with only the gzipped) but which may become more complicated if both .orb and .orb.gz files are present.
checking if glob(path.join(obsdir, "auxil/ni*.orb")) has length zero, and if it's the case using with glob(path.join(obsdir, "auxil/ni*.orb.gz"))

There are probably other ways to solve this.
Any preferences or suggestions ?

Improving SCI plots

Some suggestions to improve the science plots:

give the total clock time between start and end.
give amount of data thrown out by filtering (between prefilt and cleantfilt)
rename "Good Time Fraction" (can be confusing)
Fix table of GTI when list is too long (either reduce font, or show first and last few GTI of the list). See issue #14

Bug in --maxovershoot

I got this report from Mason Ng:

I did encounter an error with using “—maxovershoot”:

"UnicodeDecodeError: 'ascii' codec can't decode byte 0xa9 in position 1: ordinal not in range(128)”,

which if I read correctly, originated from line 324 of the script on reading the number of rows of the GTI file. It’s an odd error. I wrote in a different line [len(fits.open(gtiname3)[1].data)] just for a temporary fix, but I don’t think the flag actually made the right cuts (wanted overonly_range=0-2.0)? Anyways, on my end (locally), I just added a new “—nooverfilt” argument similar in structure to “—nounderfilt”.

Bkg plots from filtered (clean) events/gti files

I think it could also be useful to plot the BKG plots once all the filtering has been performed. This would be a useful diagnostic of what could still cause count rates variations.

But this means giving the housekeeping files as well as the clean event file...I haven't really figure out how to do this. Any suggestions?

Bug with option --usez in script ni_Htest_sortGTI.py

Weird bug with a data set happened when adding the option --usez; and only with some values of --emin and --emax.
For example, --emin=0.4 --emax=2.0 works fine, but --emin=0.3 --emax=2.0 doesn't.

Traceback (most recent call last):
  File "/Users/sebastien/Softwares/NICERsoft/scripts/ni_Htest_sortgti.py", line 305, in <module>
    hsig = sig2sigma(chi2.sf(hs,4))
  File "/Users/sebastien/Softwares/anaconda3/lib/python3.6/site-packages/pint-0.5.4+360.ga4b33ea-py3.6.egg/pint/eventstats.py", line 76, in sig2sigma
    results[isig] = fsolve(inverfc, [8], (sig,))
  File "/Users/sebastien/Softwares/anaconda3/lib/python3.6/site-packages/scipy/optimize/minpack.py", line 147, in fsolve
    res = _root_hybr(func, x0, args, jac=fprime, **options)
  File "/Users/sebastien/Softwares/anaconda3/lib/python3.6/site-packages/scipy/optimize/minpack.py", line 225, in _root_hybr
    ml, mu, epsfcn, factor, diag)
ValueError: The array returned by a function changed size between calls

Bug in new version of ni_Htest_sortGTI

There is weird behaviour in make_sn() which causes the code to crash.

In some cases, the length of the arrays gti_cts and gti_len have different sizes, with gti_cts having a size larger than gti_len by 1 extra value.

I think I have managed to track back the difference to the array gti_idx, which is sometimes different between the old and new versions of ni_Htest_sortGTI. So the TIMEZERO bug correction is causing searchsorted(t1s,times) to produce different arrays (I guess that's expected), but this somehow produces a gti_cts (from gti_cts = np.bincounts(gti_idx))with a different sizes between the two versions…

It’s happening with most ObsIDs that I have tried, but not all. I can’t figure out why…

bug in ni_Htest_sortGTI: discrepancy between grid search results and fixed energy range

I noticed another unexpected behaviour from ni_Htest_sortGTI.py.
When using the results from a grid search (i.e., optimal energy range) into a single run, the output significance, exposure time and Nb of events are different than what the grid search gives.

For example, with the following execution:
ni_Htest_sortgti.py @TEST.list newtest --gridsearch --emin=0.27 --maxemin=0.3 --minemax=1.22 --emax=1.26
gives as a result:

Maximum significance: 6.847 sigma
   obtained in 414.76 (out of 414.82 ksec)
   between 0.29 and 1.24 keV
   for 149230 events

Then, if I use this optimal energy range (0.29-1.24 keV), with the following command:
ni_Htest_sortgti.py @TEST.list newtest2 --emin=0.29 --emax=1.24
the result is:

Maximum significance: 6.628 sigma
   obtained in 406.95 ksec (out of 414.82 ksec)
   for 144175 events

But the significance, exposure time and Nb of events should be the same at the grid search run above. Note that this happens with different data sets (different pulsars).

I have found yet why this is different, and I'll investigate in the next few days.
@kerrm, if you spot something obvious, please let me know.

MaxOverShoot GTI filtering from MKF file instead of creating BKF file

@paulray
Following up on your suggestion (in #35) to enable Maximum Overshoot filtering from the MKF file instead of creating a BKF file:

I think this can easily be done in psrpipe.py with something like:

       gtiname3 = path.join(pipedir,'bkf.gti')
       mktable = Table.read(mkfile,hdu=1)
       mkf_expr = 'FPM_OVERONLY_COUNT.lt.{0}'.format(args.maxovershoot/52)
       cmd = ["maketime", mkfile, gtiname3, 'expr={0}'.format(mkf_expr),
              "compact=no", "time=TIME", "prefr=0", "postfr=0", "clobber=yes"]
       runcmd(cmd)

This appears to be working fine. See the 2 attached figures, the first one without MaxOverShoot filtering, the second with MaxOverShoot<20.0. The time binning of the MKF file is currently 1 sec.

Let me know what you think

Add --polycos option for phase calculation

The new option --polycos in the PINT photonphase script significantly improves the calculation time for large event files (nanograv/PINT#1199).

We should consider adding this option to psrpipe.py and merge.py.

Hard-coded count rate for nopulsetest in ni_Htest_sortGTI

Hi @kerrm,
I noticed that for the nopulsetest that uses only the S/N, the pred_rate was hard-coded, with a value that give very small exposure time for some data sets (e.g. for PSR J1012).
We should probably want to avoid that. Thoughts?

Improvements to GTI sort scripts

Some suggestions by Mike Wolff:

dump to text file the count rate of the selected (optimal) GTIs
show count rate light curve of the (optimal) GTIs

Other improvement:

reorganised the code a bit if possible (@kerrm)

Exception in sci_plots.py

I'm getting exceptions sometimes in sci_plots when it tries to format the stringtable.
What I see is that at this line:

NICERsoft/nicer/sci_plots.py

Line 114 in c9d6841

stringtable_duration = str(gtitable["DURATION"][:]).split('\n')

it grabs a table of the durations. Later it converts all the durations to float, but this fails because I'm getting '...' for one of the durations. I'm not sure why. Perhaps something is contracting down a longer table?

Speeding up ni_Htest_sortGTI

@kerrm, let me know how it goes, and how you are planning to do it.

I was thinking about implementing it in a multiprocessing.Pool, but if you have adjustments to the algorithms to do, I'll let you sort that you sort that out first.

Switch to logguru

PINT changed logging methods. We should change NICERsoft too:

import pint.logging                                                                                                                                                      
from loguru import logger as log                                                                                                                                         
                                                                                                                                                                         
log.remove()                                                                                                                                                             
log.add(                                                                                                                                                                 
    sys.stderr,                                                                                                                                                          
    level="WARNING",                                                                                                                                                     
    colorize=True,                                                                                                                                                       
    format=pint.logging.format,                                                                                                                                          
    filter=pint.logging.LogFilter(),                                                                                                                                     
)