diff --git a/data/modelrxte_ait_3to20keV_flux_1deg.fits b/data/modelrxte_ait_3to20keV_flux_1deg.fits new file mode 100644 index 0000000..5614c10 Binary files /dev/null and b/data/modelrxte_ait_3to20keV_flux_1deg.fits differ diff --git a/data/modelrxte_ait_3to20keV_flux_2deg.fits b/data/modelrxte_ait_3to20keV_flux_2deg.fits new file mode 100644 index 0000000..1968102 Binary files /dev/null and b/data/modelrxte_ait_3to20keV_flux_2deg.fits differ diff --git a/ridge/ridge/config.py b/ridge/ridge/config.py index f0a3bb6..b641462 100644 --- a/ridge/ridge/config.py +++ b/ridge/ridge/config.py @@ -5,6 +5,7 @@ """ datadir="../data/" proddir="../products/" +mapsdir="../products/maps/" """ Определение областей "Галактика" и "Внегалактика" """ bmax=30.0 @@ -39,3 +40,6 @@ SEM означает standatd error on mean (~RMS/sqrt(n)) """ sem_cut=98 +""" Координаты Краба """ +crab_ra = 83.6287 +crab_dec = 22.014 diff --git a/ridge/ridge/utils.py b/ridge/ridge/utils.py index 2b9d7f9..3de9c5a 100644 --- a/ridge/ridge/utils.py +++ b/ridge/ridge/utils.py @@ -63,7 +63,7 @@ def plot_best_fit(X, y, model): y_pred = model.predict(X) resid = y - y_pred err = np.sqrt(np.sum((resid)**2))/len(resid) - print("ax+b: a={:.2e}, b={:.2e}, std={:.2f}".format(a,b,err)) + print("ax+b: a={:.2e}, b={:.2e}, std={:.2e}".format(a,b,err)) """ plt.scatter(X, y) xaxis = arange(X.min(), X.max(), 0.01) diff --git a/scripts/01_bgdmodel.py b/scripts/01_bgdmodel.py index ad9cbce..8f9f8e2 100755 --- a/scripts/01_bgdmodel.py +++ b/scripts/01_bgdmodel.py @@ -7,7 +7,7 @@ import numpy as np import pandas as pd from astropy.io import fits import matplotlib.pyplot as plt -import math, sys +import math, sys, os import pickle from sklearn.linear_model import LinearRegression @@ -42,6 +42,9 @@ nrev=0 bgdmodel={} ignored_scw=[] +if not os.path.exists(proddir): + os.makedirs(proddir) + for rev in range(revmin,revmax): # if not (rev==341): diff --git a/scripts/01_crabmodel.py b/scripts/01_crabmodel.py new file mode 100755 index 0000000..b6fa9c5 --- /dev/null +++ b/scripts/01_crabmodel.py @@ -0,0 +1,307 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import pandas as pd +from astropy.io import fits +from astropy import wcs +import matplotlib.pyplot as plt +import math, sys, os +import pickle +from numpy.polynomial import Polynomial + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + + +from astropy.coordinates import SkyCoord # High-level coordinates +from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames +from astropy.coordinates import Angle, Latitude, Longitude # Angles +import astropy.units as u + +#from statsmodels.robust.scale import huber +from astropy.stats import sigma_clip +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +if not os.path.exists(proddir): + os.makedirs(proddir) + +enkey = sys.argv[1] + +sigma=3 + +# some static revs/scws to be removed +ignore_orbits=[352,834,912,1019,1021,1028,2275,2405,2493] +ignore_scws=['066600420020','066600420030','132800350010','090200390010','269500190010'] + +fn="detcnts.{}.fits".format(enkey) + +with open(proddir+fn.replace(".fits",".ignored_scw.pkl"), 'rb') as fp: + ignored_scw = pickle.load(fp) + +d = fits.getdata(datadir+fn) +df=pd.DataFrame(np.array(d).byteswap().newbyteorder()) +#print(df.columns) + +crab_crd = SkyCoord(crab_ra, crab_dec, frame=FK5(), unit="deg") + +plotme=False + +npix = 50 +sw = 30.0 # deg +pix = sw/npix # pixel size in degrees + + +mean_map = [[0.0 for i in range(npix)] for j in range(npix)] +count_map = [[0 for i in range(npix)] for j in range(npix)] + + +crabmodel={} +rota_arr=[] +a0=[] +b0=[] +a_full0=[] +b_full0=[] +b_est0=[] +err0=[] +rev0=[] + +totx=[] +toty=[] + +for i,rec in df.iterrows(): + obsid = rec['OBSID'].decode("utf-8") + if(obsid in ignore_scws): + print("Remove {}".format(obsid)) + df.drop(index=i, inplace=True) + + if (obsid in ignored_scw): + print("Skip ScW",obsid) + continue + +# accumulate full data set +for rev in range(revmin,revmax): + if(rev in ignore_orbits): + continue + if(enkey == 'E02'): + if(rev > 800): + continue + if(enkey == 'E03'): + if(rev > 800): + continue + #if(rev > 1000): + # continue + df0 = df.query('SRC > 0.0 & REV == {} & PHASE > {} & PHASE < {} & CRAB_SEP < {}'.format(rev,phmin,phmax,crab_sep_max)) + nobs=len(df0) + if not (nobs> crab_nmax): + continue + for n in df0['CRAB_SEP'].values: + totx.append(n) + for n in df0['SRC'].values: + toty.append(n) + #if(np.min(df0['SRC'])<1e-3): + # print(rev) + # sys.exit() + + +x = np.array(totx) +y = np.array(toty) +x = x.reshape((-1, 1)) +model = LinearRegression() +#model = TheilSenRegressor() +results = evaluate_model(x, y, model) +a_full,b_full,err_full = plot_best_fit(x, y, model) +if(plotme): + plot_ab(x, y, a_full, b_full, err_full, title="REGRESSION") + +# go over orbits +poly_x=[] +poly_y=[] + +ntotal=0 +nrev=0 +for rev in range(revmin,revmax): + if(rev in ignore_orbits): + continue + df0 = df.query('SRC > 0.0 & REV == {} & PHASE > {} & PHASE < {} & CRAB_SEP < {}'.format(rev,phmin,phmax,crab_sep_max)) + nobs=len(df0) + if not (nobs> crab_nmax): + continue + cen_ra = np.array(df0['RA'].values) + cen_dec = np.array(df0['DEC'].values) + print("*** Orbit ",rev) + x = np.array(df0['CRAB_SEP'].values) + y = np.array(df0['SRC'].values) + rota_deg = np.array(df0['ROTA'].values) + rota = np.array(df0['ROTA'].values) * np.pi / 180. # in radians + detx = np.cos(rota)*x + dety = np.sin(rota)*x + + for i in range(rota_deg.shape[0]): + rota_arr.append(rota_deg[i]) + + w = wcs.WCS(naxis=2) + w.wcs.crpix = [npix/2, npix/2] + w.wcs.cdelt = np.array([-pix, pix]) + w.wcs.crota = [rota_deg[i], rota_deg[i]] + w.wcs.crval = [cen_ra[i], cen_dec[i]] + w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + #w.wcs.set_pv([(2, 1, 45.0)]) + + #header = w.to_header() + #hdu = fits.PrimaryHDU(header=header) + + + pixcrd = w.wcs_world2pix([(crab_ra,crab_dec)], 1) + ix=int(np.round(pixcrd[0][0])) + iy=int(np.round(pixcrd[0][1])) + #print(ix,iy) + mean_map[ix][iy] += y[i] + count_map[ix][iy] += 1 + + x = x.reshape((-1, 1)) + model = LinearRegression() + #model = TheilSenRegressor() + results = evaluate_model(x, y, model) + a,b,err = plot_best_fit(x, y, model) + + b_est = np.mean(y - a_full*x) + + if(enkey in ['E10','E11','E12',] or a > 0.0): + filtered_data = sigma_clip(y, sigma=sigma, maxiters=10, return_bounds=True) + filtered_y = filtered_data[0] + filtered_min = filtered_data[1] + filtered_max = filtered_data[2] + b = np.mean(filtered_y) + a = 0.0 + err = np.sqrt(np.sum((b-filtered_y)**2))/len(filtered_y) + + b_est = b + a_full = a + + a_full0.append(a_full) + b_full0.append(b_full) + b_est0.append(b_est) + + #if(b>2.0e-4): + # print(rev) + # sys.exit() + + a0.append(a) + b0.append(b) + err0.append(err) + rev0.append(rev) + + + poly_x.append(rev) + poly_y.append(b_est) + + if(enkey in ['E04','E05','E06','E07','E08','E09','E10','E11','E12']): + crabmodel[rev]={'a':a_full, 'b':b_est, 'err':err} + if(plotme): + plot_ab(x, y, a_full, b_est, err, title="REGRESSION rev {}".format(rev)) + + if(enkey in ['E02', 'E03']): + crabmodel[rev]={'a':a, 'b':b, 'err':err} + if(plotme): + plot_ab(x, y, a, b, err, title="REGRESSION rev {}".format(rev)) + + print("ax+b: a={:.2e}, b={:.2e}, std={:.2e}\n".format(a,b,err)) + ntotal=ntotal+nobs + nrev=nrev+1 + +print("Orbits: {}, obs: {}".format(nrev,ntotal)) + +for i in range(npix): + for j in range(npix): + if(count_map[i][j]>0): + mean_map[i][j]=mean_map[i][j]/count_map[i][j] + +w = wcs.WCS(naxis=2) +w.wcs.crpix = [npix/2, npix/2] +w.wcs.cdelt = np.array([-pix, pix]) +w.wcs.crval = [0.0, 0.0] +w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + +header = w.to_header() +hdu = fits.PrimaryHDU(header=header) + +hdu.data=mean_map +hdu.writeto(proddir+fn.replace(".fits",".crab_mean_map.fits"), overwrite=True) +hdu.data=count_map +hdu.writeto(proddir+fn.replace(".fits",".crab_count_map.fits"), overwrite=True) + +npoly=4 +if(enkey in ['E11','E12',]): + npoly=0 + +z = np.polyfit(poly_x, poly_y, npoly) + +p = np.poly1d(z) +poly_z=[] +for t in poly_x: + poly_z.append(p(t)) +plt.scatter(poly_x, poly_y) +plt.plot(poly_x, poly_z, color='r') +plt.title("Crab detector count rate evolution") +plt.ylabel("Crab count rate cts/s/pix") +plt.xlabel("INTEGRAL orbit") +#plt.show() +plt.savefig(proddir+fn.replace(".fits",".crab_rate.png")) + +indices = sorted( + range(len(a0)), + key=lambda index: a0[index] +) + +coldefs = fits.ColDefs([ + #fits.Column(name='OBSID', format='11A', array=[obs_id[index] for index in indices]), + #fits.Column(name='RA', format='D', unit='deg', array=[ra[index] for index in indices]), + #fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]), + #fits.Column(name='LON', format='D', unit='deg', array=[lon0[index] for index in indices]), + #fits.Column(name='LAT', format='D', unit='deg', array=[lat0[index] for index in indices]), + fits.Column(name='REV', format='J', unit='', array=[rev0[index] for index in indices]), + #fits.Column(name='MJD', format='D', unit='', array=[mjd0[index] for index in indices]), + #fits.Column(name='PHASE', format='D', unit='', array=[phase0[index] for index in indices]), + #fits.Column(name='CLEAN', format='D', unit='cts/s', array=[clean0[index] for index in indices]), + #fits.Column(name='MODEL', format='D', unit='cts/s', array=[model0[index] for index in indices]), + #fits.Column(name='RESID', format='D', unit='cts/s', array=[resid0[index] for index in indices]), + fits.Column(name='A', format='D', unit='', array=[a0[index] for index in indices]), + fits.Column(name='B', format='D', unit='', array=[b0[index] for index in indices]), + fits.Column(name='ERR', format='D', unit='', array=[err0[index] for index in indices]), + fits.Column(name='A_FULL', format='D', unit='', array=[a_full0[index] for index in indices]), + fits.Column(name='B_FULL', format='D', unit='', array=[b_full0[index] for index in indices]), + fits.Column(name='B_EST', format='D', unit='', array=[b_est0[index] for index in indices]), + fits.Column(name='B_POLY', format='D', unit='', array=[poly_z[index] for index in indices]), +]) + +fout = fn.replace(".fits",".crabmodel.fits") +hdu = fits.BinTableHDU.from_columns(coldefs, name='GRXE') +hdu.header['MISSION'] = ('INTEGRAL', '') +hdu.header['TELESCOP'] = ('IBIS', '') +hdu.header['INSTITUT'] = ('IKI', 'Affiliation') +hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') +hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') +#hdu.add_checksum() +print(hdu.columns) +hdu.writeto(proddir+fout, overwrite=True) + + +with open(proddir+fn.replace(".fits",".crabmodel.pkl"), 'wb') as fp: + pickle.dump([crabmodel, z], fp, protocol=pickle.HIGHEST_PROTOCOL) + + +with open(proddir+fn.replace(".fits",".rota.dat"), 'w') as fp: + for i in range(len(rota_arr)): + fp.write("{} {}\n".format(i,rota_arr[i])) + diff --git a/scripts/01_crabmodel.sh b/scripts/01_crabmodel.sh new file mode 100755 index 0000000..7fb25d2 --- /dev/null +++ b/scripts/01_crabmodel.sh @@ -0,0 +1,11 @@ +./01_crabmodel.py E02 +./01_crabmodel.py E03 +./01_crabmodel.py E04 +./01_crabmodel.py E05 +./01_crabmodel.py E06 +./01_crabmodel.py E07 +./01_crabmodel.py E08 +./01_crabmodel.py E09 +./01_crabmodel.py E10 +./01_crabmodel.py E11 +./01_crabmodel.py E12 diff --git a/scripts/02_grxe_resid.py b/scripts/02_grxe_resid.py new file mode 100755 index 0000000..3678652 --- /dev/null +++ b/scripts/02_grxe_resid.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import pandas as pd +from astropy.io import fits +from astropy.table import Table, Column +import matplotlib.pyplot as plt +import math, sys +import pickle + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + +#from statsmodels.robust.scale import huber +from astropy.stats import sigma_clip +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +enkey = sys.argv[1] +outkey = sys.argv[2] + +fn="detcnts.{}.fits".format(enkey) + +with open(proddir+fn.replace(".fits",".pkl"), 'rb') as fp: + bgdmodel = pickle.load(fp) + +with open(proddir+fn.replace(".fits",".ignored_scw.pkl"), 'rb') as fp: + ignored_scw = pickle.load(fp) + +with open(proddir+fn.replace(".fits",".crabmodel.pkl"), 'rb') as fp: + crabmodel, z = pickle.load(fp) + p = np.poly1d(z) + #print(crabmodel) + +crab_rev_max = np.max(list(crabmodel.keys())) +print("Crab is defined untill orbit {}".format(crab_rev_max)) + +with fits.open(datadir+fn) as hdul: + data=hdul[1].data + +#print(data.columns) + +rev = data.field('rev') +mjd = data.field('mjd') +clean = data.field('clean') +phase = data.field('phase') + + +rev0=[] +phase0=[] +clean0=[] +model0=[] +resid0=[] # residuals in cts/s +grxe0=[] # mCrab +crab0=[] # Crab count rate +mjd0=[] +a0=[] +b0=[] +err0=[] +lon0=[] +lat0=[] + +d = fits.getdata(datadir+fn) +df = pd.DataFrame(np.array(d).byteswap().newbyteorder()) + +# BKG +if(outkey == 'BKG'): + df = df.query('CLEAN > 0.0 & ( abs(LAT) > {} | abs(LON) > {}) & PHASE > {} & PHASE < {}'.format(bmax,lmax,phmin,phmax)) + +# GAL +if(outkey=='GAL'): + df = df.query('CLEAN > 0.0 & abs(LAT) < {} & abs(LON) < {} & PHASE > {} & PHASE < {}'.format(bmax,lmax,phmin,phmax)) + +# ALL +if(outkey=='ALL'): + df = df.query('CLEAN > 0.0 & PHASE > {} & PHASE < {}'.format(phmin,phmax)) + +for i, row in df.iterrows(): + orbit=row['REV'] + obsid=row['OBSID'].decode("UTF-8") + + if not (orbit > revmin and orbit < revmax): + print("Skip orbit",orbit,row['OBSID']) + continue + + if not (orbit < crab_rev_max): + print("Skip orbit",orbit,row['OBSID']) + continue + + if (obsid in ignored_scw): + print("Skip ScW",obsid) + continue + + a = bgdmodel[orbit]['a'] + b = bgdmodel[orbit]['b'] + err = bgdmodel[orbit]['err'] + m = a*row['PHASE']+b + + clean0.append(clean[i]) + mjd0.append(mjd[i]) + model0.append(m) + resid0.append(clean[i]-m) + grxe0.append(1000*(clean[i]-m)/p(orbit)) + crab0.append(p(orbit)) + + a0.append(a) + b0.append(b) + err0.append(err) + phase0.append(row['PHASE']) + rev0.append(orbit) + lon0.append(row['LON']) + lat0.append(row['LAT']) + + +indices = sorted( + range(len(mjd0)), + key=lambda index: mjd0[index] +) + +coldefs = fits.ColDefs([ + #fits.Column(name='OBSID', format='11A', array=[obs_id[index] for index in indices]), + #fits.Column(name='RA', format='D', unit='deg', array=[ra[index] for index in indices]), + #fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]), + fits.Column(name='LON', format='D', unit='deg', array=[lon0[index] for index in indices]), + fits.Column(name='LAT', format='D', unit='deg', array=[lat0[index] for index in indices]), + fits.Column(name='REV', format='J', unit='', array=[rev0[index] for index in indices]), + fits.Column(name='MJD', format='D', unit='', array=[mjd0[index] for index in indices]), + fits.Column(name='PHASE', format='D', unit='', array=[phase0[index] for index in indices]), + fits.Column(name='CLEAN', format='D', unit='cts/s', array=[clean0[index] for index in indices]), + fits.Column(name='MODEL', format='D', unit='cts/s', array=[model0[index] for index in indices]), + fits.Column(name='RESID', format='D', unit='cts/s', array=[resid0[index] for index in indices]), + fits.Column(name='GRXE', format='D', unit='mCrab', array=[grxe0[index] for index in indices]), + fits.Column(name='CRAB', format='D', unit='cts/s', array=[crab0[index] for index in indices]), + fits.Column(name='A', format='D', unit='', array=[a0[index] for index in indices]), + fits.Column(name='B', format='D', unit='', array=[b0[index] for index in indices]), + fits.Column(name='ERR', format='D', unit='', array=[err0[index] for index in indices]), +]) + +fout = fn.replace(".fits",".{}.resid.fits".format(outkey)) +hdu = fits.BinTableHDU.from_columns(coldefs, name='GRXE') +hdu.header['MISSION'] = ('INTEGRAL', '') +hdu.header['TELESCOP'] = (outkey, '') +hdu.header['INSTITUT'] = ('IKI', 'Affiliation') +hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') +hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') +#hdu.add_checksum() +print(hdu.columns) +hdu.writeto(proddir+fout, overwrite=True) + +with fits.open(proddir+fout, mode='update') as hdus: + hdus[1].add_checksum() + + + diff --git a/scripts/02_grxe_resid.sh b/scripts/02_grxe_resid.sh new file mode 100755 index 0000000..189d3a6 --- /dev/null +++ b/scripts/02_grxe_resid.sh @@ -0,0 +1,35 @@ +./02_grxe_resid.py E02 GAL +./02_grxe_resid.py E03 GAL +./02_grxe_resid.py E04 GAL +./02_grxe_resid.py E05 GAL +./02_grxe_resid.py E06 GAL +./02_grxe_resid.py E07 GAL +./02_grxe_resid.py E08 GAL +./02_grxe_resid.py E09 GAL +./02_grxe_resid.py E10 GAL +./02_grxe_resid.py E11 GAL +./02_grxe_resid.py E12 GAL + +./02_grxe_resid.py E02 BKG +./02_grxe_resid.py E03 BKG +./02_grxe_resid.py E04 BKG +./02_grxe_resid.py E05 BKG +./02_grxe_resid.py E06 BKG +./02_grxe_resid.py E07 BKG +./02_grxe_resid.py E08 BKG +./02_grxe_resid.py E09 BKG +./02_grxe_resid.py E10 BKG +./02_grxe_resid.py E11 BKG +./02_grxe_resid.py E12 BKG + +./02_grxe_resid.py E02 ALL +./02_grxe_resid.py E03 ALL +./02_grxe_resid.py E04 ALL +./02_grxe_resid.py E05 ALL +./02_grxe_resid.py E06 ALL +./02_grxe_resid.py E07 ALL +./02_grxe_resid.py E08 ALL +./02_grxe_resid.py E09 ALL +./02_grxe_resid.py E10 ALL +./02_grxe_resid.py E11 ALL +./02_grxe_resid.py E12 ALL diff --git a/scripts/03_grxe_map.py b/scripts/03_grxe_map.py new file mode 100755 index 0000000..837ec89 --- /dev/null +++ b/scripts/03_grxe_map.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import numpy.ma as ma +import pandas as pd +from astropy.wcs import WCS +from astropy import wcs +from astropy.io import fits +from astropy.table import Table, Column +import matplotlib.pyplot as plt +import math, sys, os +import pickle + +from astropy.coordinates import SkyCoord # High-level coordinates +from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames +from astropy.coordinates import Angle, Latitude, Longitude # Angles +import astropy.units as u + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + +from astropy.stats import sigma_clip +from astropy.stats import sigma_clipped_stats + +from scipy.stats import norm +from scipy.stats import describe +from scipy.stats import sem + + +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +enkey = sys.argv[1] +key = sys.argv[2] +#key="ALL" +fn='detcnts.{}.{}.resid.fits'.format(enkey,key) + +d = fits.getdata(proddir+fn) +df=pd.DataFrame(np.array(d).byteswap().newbyteorder()) +#print(df.columns) + +#df = df.query('abs(LAT) < {} & abs(LON) < {}'.format(5,5)) + +n_bins = 80 +minmax=[-300,800] +sigma=3 +maxiters=10 + +modelrxte="modelrxte_ait_3to20keV_flux_2deg.fits" +hdulist = fits.open(datadir+modelrxte) +w = wcs.WCS(hdulist[0].header) +smap =hdulist[0].data + +sx=int(hdulist[0].header['NAXIS1']) +sy=int(hdulist[0].header['NAXIS2']) + +# fill AITOF map indexes +ds9x=[] +ds9y=[] +for i,row in df.iterrows(): + lon=row['LON'] + lat=row['LAT'] + world = SkyCoord(lon,lat, frame=Galactic, unit="deg") + ra=world.fk5.ra.deg + dec=world.fk5.dec.deg + pixcrd = w.wcs_world2pix([(lon,lat)], 1) + y=int(pixcrd[0][0]) + x=int(pixcrd[0][1]) + ds9x.append(x) + ds9y.append(y) + #print(x,y,smap[y-1,x-1]) +df['DS9Y']=ds9x +df['DS9X']=ds9y + +mean_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) +sign_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) +sem_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) +cnt_map = np.array([[0 for i in range(sx)] for j in range(sy)]) + +for i in range(sx): + for j in range(sy): + world = w.wcs_pix2world([(i+1,j+1)], 1) + lon = world[0][0] + lat = world[0][1] + if(np.isnan(lon) or np.isnan(lat)): + continue + ds9i=i+1 + ds9j=j+1 + df0 = df.query('DS9X == {} & DS9Y == {}'.format(ds9i,ds9j)) + if (len(df0) <= nscw_min): + continue + + # check coordinates + #print("***",i+1,j+1,lon,lat,smap[j][i]) + #for i0,row0 in df0.iterrows(): + # print(row0['LON'],row0['LAT'],row0['GRXE']) + + grxe = np.array(df0['GRXE']) + + sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=maxiters) + filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=maxiters, return_bounds=True) + filtered_grxe = filtered_data[0] + filtered_min = filtered_data[1] + filtered_max = filtered_data[2] + + # final error on flux measurement ~RMS/sqrt(n) + sg_sem = sem(grxe[filtered_grxe.mask==False]) + + #print("Sigma clipping: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f}".format(sg_mean, sg_med, sg_std, sg_sem)) + + #plt.hist(grxe, bins=n_bins, range=minmax) + #plt.hist(grxe[filtered_grxe.mask], bins=n_bins, range=minmax) + #plt.show() + + mean_map[j][i] = sg_mean + sem_map[j][i] = sg_sem + sign_map[j][i] = sg_mean/sg_sem + cnt_map[j][i] = df0.shape[0] + + +# Calculate the percentiles across the x and y dimension +perc = np.percentile(sem_map, sem_cut, axis=(0, 1), keepdims=False) + +print("{}: {}% cut of SEM map: {:.2f} mCrab".format(enkey,sem_cut,perc)) + +idx=np.where(sem_map > perc) +mean_map[idx]=0.0 +sem_map[idx]=0.0 +cnt_map[idx]=0 +sign_map[idx]=0.0 + +if not os.path.exists(mapsdir): + os.makedirs(mapsdir) + +hdulist[0].data=mean_map +hdulist.writeto(mapsdir+fn.replace(".fits",".mean.fits"),overwrite=True) + +hdulist[0].data=sem_map +hdulist.writeto(mapsdir+fn.replace(".fits",".error.fits"),overwrite=True) + +hdulist[0].data=cnt_map +hdulist.writeto(mapsdir+fn.replace(".fits",".cnt.fits"),overwrite=True) + +hdulist[0].data=sign_map +hdulist.writeto(mapsdir+fn.replace(".fits",".sign.fits"),overwrite=True) diff --git a/scripts/03_grxe_map.sh b/scripts/03_grxe_map.sh new file mode 100644 index 0000000..b3b7863 --- /dev/null +++ b/scripts/03_grxe_map.sh @@ -0,0 +1,11 @@ +./03_grxe_map.py E02 ALL +./03_grxe_map.py E03 ALL +./03_grxe_map.py E04 ALL +./03_grxe_map.py E05 ALL +./03_grxe_map.py E06 ALL +./03_grxe_map.py E07 ALL +./03_grxe_map.py E08 ALL +./03_grxe_map.py E09 ALL +./03_grxe_map.py E10 ALL +./03_grxe_map.py E11 ALL +./03_grxe_map.py E12 ALL diff --git a/scripts/README.md b/scripts/README.md index 0968533..0bd67c5 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -3,41 +3,20 @@ source /venv/bin/activate.csh ``` +Номер скрипта означает последовательность выполнения. Результаты работы скриптов помещаются в ```ridge/products/```. -### 01_init_events.py +### 01_bgdmodel.py / 01_bgdmodel.sh -Создает начальные списки событий и помещает их в ```uds/data/processed``` -Оригинальные файлы со списками событий задаются в файлах ```uds/data/evtlists/*.txt``` +Создает модель фона -### 02_merge_events.py +### 01_crabmodel.py / 01_crabmodel.sh -Создает объедененный список событий и помещает его в ```uds/products```. Этот список событий нужен, в основном для извлечения спектров с помощью ```srctool```. +Создает модель скорости счета Краба -Попутно этот скрипт унифицирует оригинальные списки событий для последующей обработки. А именно, корректируются слова OBS_MODE=POINING/SURVEY в зависимости от типа наблюдения и производится центрирование на одни и те же координаты с помощью команды ```radec2xy```. +### 02_grxe_resid.py / 01_grxe_resid.sh -Для запуска адаптивного сглаживания ```do_adapt = True``` требуется запустить окружение ```ciao```, так как нужна команда ```dmimgadapt``` +Вычисляет остатки после вычитания модели из данных в единицах мКраб -``` -conda activate ciao-4.15 -source /eSASS4EDR/bin/esass-init.csh -``` +### 02_grxe_map.py / 01_grxe_map.sh -### 03_init_obs.py - -1) Подготавливает списки событий в разных энергетических диапазонах. -2) Запускает ```erbox``` в три этапа, чтобы получить рабочий список источников для ```ermldet```. -3) Запускает ```ermldet``` -4) Делает кросс-корреляцию с каталогом Gaia-unWISE ```do_cross_match=True```, которая создает три файла: ```.cross``` -- все пересечения, и ```.ref``` / ```.src``` -- входные каталоги для последующей команды ```wcs_match``` -5) Делает матрицу преобразования координат и корректирует списки событий. Для запуска команд```wcs_match/wcs_update``` требуется запустить окружение ```ciao``` (см. выше) - -### 04_mosaics.py - -Создает сборные изображения (мозайки) в разных энергетических диапазонах. - -### 05_scrtool.py - -Запускает scrtool для самого широкого канала 0.2-10 кэВ, чтобы спектры имели самое полное покрытие по энергиям. Список источников берется из 0.3-2.3 кэВ. - -Вычисляет ECF для всех диапазонов. - -Делает принудительную фотометрию в выбранных каналах (параметр```forced=True```). Внимание! ermldet из eSASS4EDR не делает ассимитричные ошибки на потоки. Мы запускаем более последнюю версию ermldet (v1.56/2.18 esass_200412 Jul 2 12:04:46 2022). Для этого используется параметр ```local_run=True```, который высвечивает какую команду надо запустить на другой машине и ждет ввода. \ No newline at end of file +Создает карту остатков в единицах мКраб. \ No newline at end of file