update

data/README.md

@@ -34,5 +34,5 @@ This directory contains preprocessed ISGRI detector count rate cleaned from the
 
 ### cobe_ibis_resp_lon.dat
 
-This file comes from [Krivonos+ 2007](https://ui.adsabs.harvard.edu/abs/2007A%26A...463..957K/abstract) paper
+This file comes from [Krivonos+ 2007](https://ui.adsabs.harvard.edu/abs/2007A%26A...463..957K/abstract) paper, and is used in Figure 1 of the current paper.
 

data/xspec/GC06.xcm

@@ -1,12 +0,0 @@
-method leven 10 0.01
-abund angr
-xsect vern
-cosmo 70 0 0.73
-xset delta 0.01
-systematic 0
-model  atable{../../data/polarmodel.fits} + powerlaw
-        0.73486          1          0          0          3          3
-    2.27115e-09       0.01          0          0      1e+20      1e+24
-        1.23053       0.01         -3         -2          9         10
-      0.0103601       0.01          0          0      1e+20      1e+24
-bayes off

data/xspec/GC10.xcm

@@ -1,12 +0,0 @@
-method leven 10 0.01
-abund angr
-xsect vern
-cosmo 70 0 0.73
-xset delta 0.01
-systematic 0
-model  atable{../../data/polarmodel.fits} + powerlaw
-       0.713934          1          0          0          3          3
-    2.27522e-09       0.01          0          0      1e+20      1e+24
-        1.80029       0.01         -3         -2          9         10
-       0.138919       0.01          0          0      1e+20      1e+24
-bayes off

data/xspec/LON+20.xcm

@@ -1,10 +0,0 @@
-method leven 10 0.01
-abund angr
-xsect vern
-cosmo 70 0 0.73
-xset delta 0.01
-systematic 0
-model  atable{../../data/polarmodel.fits}
-       0.724678          1          0          0          3          3
-    1.03624e-09       0.01          0          0      1e+20      1e+24
-bayes off

data/xspec/LON+20_bremms.xcm

@@ -1,18 +0,0 @@
-method leven 10 0.01
-abund angr
-xsect vern
-cosmo 70 0 0.73
-xset delta 0.01
-systematic 0
-model  cflux*bremss + cflux*powerlaw
-             20       -0.1          0          0      1e+06      1e+06
-             60       -0.1          0          0      1e+06      1e+06
-        -9.0642       0.01       -100       -100        100        100
-        20.3536          1     0.0001     0.0001        100        200
-          1e-05         -1          0          0      1e+20      1e+24
-             80       -0.1          0          0      1e+06      1e+06
-            200       -0.1          0          0      1e+06      1e+06
-       -11.2811       0.01       -100       -100        100        100
-           -1.5         -1         -3         -2          9         10
-          1e-05         -1          0          0      1e+20      1e+24
-bayes off

models/xspec/cutoffpl_gb.xcm

@@ -0,0 +1,19 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0
+model  cflux*cutoffpl + cflux*powerlaw
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.01266       0.01       -100       -100        100        100
+              0         -1         -3         -2          9         10
+        11.1866       0.01       0.01          1        500        500
+              1         -1          0          0      1e+20      1e+24
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.31022       0.01       -100       -100        100        100
+           1.55         -1         -3         -2          9         10
+              1         -1          0          0      1e+20      1e+24
+bayes off

models/xspec/cutoffpl_l+20.xcm

@@ -0,0 +1,19 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0
+model  cflux*cutoffpl + cflux*powerlaw
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.36476       0.01       -100       -100        100        100
+              0         -1         -3         -2          9         10
+        11.1662       0.01       0.01          1        500        500
+              1         -1          0          0      1e+20      1e+24
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.75981       0.01       -100       -100        100        100
+           1.55         -1         -3         -2          9         10
+              1         -1          0          0      1e+20      1e+24
+bayes off

models/xspec/cutoffpl_l-20.xcm

@@ -0,0 +1,19 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0
+model  cflux*cutoffpl + cflux*powerlaw
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.36741       0.01       -100       -100        100        100
+              0         -1         -3         -2          9         10
+        11.2463       0.01       0.01          1        500        500
+              1         -1          0          0      1e+20      1e+24
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.65316       0.01       -100       -100        100        100
+           1.55         -1         -3         -2          9         10
+              1         -1          0          0      1e+20      1e+24
+bayes off

models/xspec/load_gb_cutoffpl.xcm

@@ -0,0 +1,9 @@
+data GB.spec.pha
+setpl en
+cpd /xs
+setpl reb 2
+systematic 0
+setpl reb 3 3
+@cutoffpl_gb
+pl eeuf delchi
+

models/xspec/load_gb_polar.xcm

@@ -0,0 +1,9 @@
+data GB.spec.pha
+setpl en
+cpd /xs
+setpl reb 2
+systematic 0
+setpl reb 3 3
+@polar_gb_cflux
+pl eeuf delchi
+

models/xspec/load_l+20_cutoffpl.xcm

@@ -0,0 +1,8 @@
+data  L+20.spec.pha
+setpl en
+cpd /xs
+@cutoffpl_l+20
+setpl reb 3 3
+pl eeuf delchi
+#steppar 9  -10 -12 100
+#new 9 -10.34

models/xspec/load_l+20_polar.xcm

@@ -0,0 +1,8 @@
+data  L+20.spec.pha
+setpl en
+cpd /xs
+@polar_l+20_cflux
+setpl reb 3 3
+pl eeuf delchi
+#steppar 9  -10 -12 100
+#new 9 -10.34

models/xspec/load_l-20_cutoffpl.xcm

@@ -0,0 +1,8 @@
+data  L-20.spec.pha
+setpl en
+cpd /xs
+@cutoffpl_l-20
+setpl reb 3 3
+pl eeuf delchi
+#steppar 9  -10 -12 100
+#new 9 -10.34

models/xspec/load_l-20_polar.xcm

@@ -0,0 +1,8 @@
+data  L-20.spec.pha
+setpl en
+cpd /xs
+@polar_l-20_cflux
+setpl reb 2 2
+pl eeuf delchi
+#steppar 9  -10 -12 100
+#new 9 -10.34

models/xspec/polar_gb_cflux.xcm

@@ -0,0 +1,18 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0.05
+model  cflux*atable{polarmodel.fits} + cflux*powerlaw
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.01521       0.01       -100       -100        100        100
+       0.657219       0.01          0          0          3          3
+              1         -1          0          0      1e+20      1e+24
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+        -9.1965       0.01       -100       -100        100        100
+           1.55         -1         -3         -2          9         10
+              1         -1          0          0      1e+20      1e+24
+bayes off

models/xspec/polar_l+20_cflux.xcm

@@ -0,0 +1,18 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0.05
+model  cflux*atable{polarmodel.fits} + cflux*powerlaw
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.28956       0.01       -100       -100        100        100
+       0.734431       0.01          0          0          3          3
+              1         -1          0          0      1e+20      1e+24
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.95672       0.01       -100       -100        100        100
+           1.55         -1         -3         -2          9         10
+              1         -1          0          0      1e+20      1e+24
+bayes off

models/xspec/polar_l-20_cflux.xcm

@@ -0,0 +1,18 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0.05
+model  cflux*atable{polarmodel.fits} + cflux*powerlaw
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.25148       0.01       -100       -100        100        100
+        0.79495       0.01          0          0          3          3
+              1         -1          0          0      1e+20      1e+24
+             30       -0.1          0          0      1e+06      1e+06
+             80       -0.1          0          0      1e+06      1e+06
+       -9.96087       0.01       -100       -100        100        100
+           1.55         -1         -3         -2          9         10
+              1         -1          0          0      1e+20      1e+24
+bayes off

ridge/ridge/config.py

@@ -3,6 +3,7 @@ All scripts must be ru in ridge/scripts directory
 """
 
 # All directory paths below must be ended by /
+modelsdir="../models/"
 datadir="../data/"
 proddir="../products/"
 mapsdir="../products/Maps/"
@@ -124,18 +125,10 @@ skyreg={
     #'Geminga':{'lon':-164.866,'lat':4.265811,'wlon':10.0,'wlat':10.0},
 }
 
-""" Bootstrap settings """
-
-# Number of simulations
-nsim=100
-
-# Fraction of data, randomly selected during each simulation
-simfrac=10 
-
 sigma=3
 
 # template for map
-modelrxte="grxe_model/modelrxte_ait_3to20keV_flux_2deg.fits"
+modelrxte="grxe/modelrxte_ait_3to20keV_flux_2deg.fits"
 
 # Galactic longitude profile
 lon_max=173.0

ridge/ridge/utils.py

@@ -303,8 +303,7 @@ def get_spec(df, grxe_err_cut=None, skey=None, enkey=None, sigma=3, n_bins = 60,
     resid = np.array(df['RESID'])
     if(plotme):
         print("get_spec: Initial size: {}".format(len(grxe)))
-    #print(grxe)
+
-    #print(grxe_err)
     
     perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False)
     #print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(skey,enkey,grxe_err_cut,perc))
@@ -417,6 +416,7 @@ def get_spec(df, grxe_err_cut=None, skey=None, enkey=None, sigma=3, n_bins = 60,
         ydata = n
         (mu, sg) = norm.fit(grxe)
         
+        if(plotme==True):
             print("Initial Gaiss fit: mu={:.2f} sigma={:.2f}".format(mu,sg))
         
         y_peak = norm.pdf(mu, mu, sg)
@@ -443,8 +443,8 @@ def get_spec(df, grxe_err_cut=None, skey=None, enkey=None, sigma=3, n_bins = 60,
             
         #pfit, perr = fit_curvefit(params, xdata, ydata, const_gaussian_ff)
         pfit, perr = fit_leastsq(params, xdata, ydata, const_gaussian_ff)
-        print(pfit)
+        #print(pfit)
-        print(perr)
+        #print(perr)
 
         #pfit, perr = fit_leastsq(params, xdata, ydata, double_gaussian2_ff)
         #pfit, perr = fit_curvefit(params, xdata, ydata, double_gaussian_ff)

scripts/00_stats.py

@@ -27,7 +27,7 @@ from ridge.utils import *
 from ridge.config import *
 
 #enkey = sys.argv[1]
-enkey="E01"
+enkey="A01"
 inkey="ALL"
 
 fn="detcnts.{}.fits".format(enkey)
@@ -63,7 +63,7 @@ print("Total {} ScWs, {:.1f} Ms".format(df_tot.shape[0], np.sum(df_tot['TEXP'])/
 
 
 # read ignored orbits
-with open(proddir+'detcnts.B21.ignored_rev.resid.pkl', 'rb') as fp:
+with open(ignored_rev_file, 'rb') as fp:
     ignored_rev = pickle.load(fp)
     print("{} orbits ignored".format(len(ignored_rev)))
     ign=ignored_rev.tolist()

scripts/03_grxe_flux.py

@@ -39,7 +39,6 @@ inkey="ALL"
 
 
 
-sigma=3
 plotme=False
 
 ebands0={
@@ -51,7 +50,6 @@ ebands0={
 
 if len(sys.argv) > 1:
     skeys = [sys.argv[1]]
-    #simfrac = int(sys.argv[2])
 else:
     skeys = list(skyreg.keys())
 
@@ -67,24 +65,13 @@ with open(ignored_rev_file, 'rb') as fp:
     print("{} orbits ignored".format(len(ignored_rev)))
     ign=ignored_rev.tolist()
 
-
-
-
-
 for skey in skeys:
     if not skey in skyreg.keys():
         print("{} not found in {}".format(skey,list(skyreg.keys())))
         sys.exit()
 
-    # generate array for bootstrap
-    ebands_sim={}
-    for enkey in ebands0.keys():
-        ebands_sim[enkey] = []
 
     for enkey in ebands0.keys():
-        #bkg_fn="detcnts.{}.BKG.resid.fits".format(enkey,inkey)
-        #syserr, bkg_sem = get_syserror(proddir+bkg_fn)
-        
         fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
         print("Reading {}".format(proddir+fn))
         dat = Table.read(proddir+fn, unit_parse_strict='silent')
@@ -113,13 +100,6 @@ for skey in skeys:
         ebands0[enkey]=[sg_mean,sg_sem]
 
 
-        
-        nsel = int(df.shape[0]*simfrac/100)
-        for n in range(nsim):
-            df0=df.sample(nsel)
-            sg_mean,sg_sem,skew_val,skew_err = get_spec(df0, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey)
-            ebands_sim[enkey].append(sg_mean)           
-    ###
     fspec="{}{}.dat".format(fluxdir,skey)
     with open(fspec, 'w') as fp:
         for enkey in ebands0.keys():
@@ -129,33 +109,6 @@ for skey in skeys:
             fp.write("{} {} {} {:.6f} {:.6f}\n".format(skey,enkey,bands[enkey],flux,err))    
     
 
-    ###
+
-    fspec="{}{}.sim{:02d}.dat".format(fluxdir,skey,simfrac)
-    with open(fspec, 'w') as fp:
-        for enkey in ebands_sim.keys():
-            data=ebands_sim[enkey]
-
-            (mu, sg) = norm.fit(data)
-            fp.write("{:02d} {} {} {} {:.6f} {:.6f}\n".format(simfrac,skey,enkey,bands[enkey],mu,sg))    
-            print("[BOOT] {}: {} {:.6f} {:.6f}".format(skey,enkey,mu,sg))
-
-            if(plotme):
-                n, bins, patches = plt.hist(data, 60, density=True, facecolor='green', alpha=0.75)
-                # add a 'best fit' line
-                y = norm.pdf(bins, mu, sg)
-                l = plt.plot(bins, y, 'r--', linewidth=2)
-                #plot
-                plt.axvline(mu, color="black")
-                plt.axvline(ebands0[enkey][0], color="black", linestyle="dashed")
-                #plt.axvline(mu+sg_sem, color="black", linestyle="dashed")
-                #plt.axvline(mu-sg_sem, color="black", linestyle="dashed")
-                plt.axvline(mu+sg, color="blue", linestyle="dashed")
-                plt.axvline(mu-sg, color="blue", linestyle="dashed")
-
-                plt.xlabel('Flux, mCrab')
-                plt.ylabel('Probability')
-                plt.title("[BOOT] {}: {:.2f} {:.2f}".format(enkey, mu, sg))
-                plt.grid(True)
-                plt.show()
             
 

scripts/03_grxe_galprof.py

@@ -113,14 +113,6 @@ for i in range(lon_nbin):
     nsel = int(df0.shape[0]*simfrac/100)
     print("lon {:.2f} ".format(glon[i]),"nsel=",nsel,df0.shape[0])
 
-    """
-    for n in range(nsim):
-        df1=df0.sample(nsel)
-        sg_mean1,sg_sem1,skew_val,skew_err = get_spec(df1, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
-        mean_sim[dkey].append(sg_mean1)   
-    """
-    
-        
     #print('sg_sem',sg_sem)
     mean_map[i] = sg_mean
     sem_map[i]  = sg_sem
@@ -130,7 +122,7 @@ for i in range(lon_nbin):
         lon_arr.append(row['LON'])
         lat_arr.append(row['LAT'])
 
-#sys.exit()
+
 
         
 if not os.path.exists(profdir):

scripts/03_grxe_galprof.sh

@@ -1,4 +1,3 @@
-./03_grxe_galprof.py E01 21
+./03_grxe_galprof.py A01 21
-./03_grxe_galprof.py E13 21
+./03_grxe_galprof.py A02 21
-./03_grxe_galprof.py E14 21
+./03_grxe_galprof.py A03 21
-

scripts/03_grxe_galprof_plot.py

@@ -6,7 +6,7 @@ __copyright__   = "Space Research Institute (IKI)"
 from astropy.table import Table, Column
 from matplotlib import ticker
 import matplotlib.pyplot as plt
-import sys
+import sys, os
 
 from ridge.utils import *
 from ridge.config import *
@@ -14,15 +14,15 @@ from ridge.config import *
 
 scale = 100.0
 
-fn="detcnts.E01.ALL.resid.galprof.fits"
+fn="detcnts.A01.ALL.resid.galprof.fits"
 dat = Table.read(profdir+fn, unit_parse_strict='silent')
 df1 = dat.to_pandas().sort_values(by=['LON1'])
 
-fn="detcnts.E14.ALL.resid.galprof.fits"
+fn="detcnts.A02.ALL.resid.galprof.fits"
 dat = Table.read(profdir+fn, unit_parse_strict='silent')
 df2 = dat.to_pandas().sort_values(by=['LON1'])
 
-fn="detcnts.E13.ALL.resid.galprof.fits"
+fn="detcnts.A03.ALL.resid.galprof.fits"
 dat = Table.read(profdir+fn, unit_parse_strict='silent')
 df3 = dat.to_pandas().sort_values(by=['LON1'])
 
@@ -96,5 +96,10 @@ ax2.set_ylabel('GRXE flux, x100 mCrab',fontsize=14, fontweight='normal')
 #plt.xscale('linear')
 #plt.yscale('linear')
 
-plt.savefig(profdir+'galprof.png', bbox_inches='tight')
+if not os.path.exists(figdir):
+    os.makedirs(figdir)
+
+filename=figdir+'galprof.png'
+plt.savefig(filename, bbox_inches='tight')
 plt.close(fig)
+print("Result is saved as {}".format(filename))

scripts/03_grxe_map.py

@@ -64,7 +64,7 @@ with open(ignored_rev_file, 'rb') as fp:
 ign=ignored_rev.tolist()
 
 
-hdulist = fits.open(datadir+modelrxte)
+hdulist = fits.open(modelsdir+modelrxte)
 w = wcs.WCS(hdulist[0].header)
 smap =hdulist[0].data
 
@@ -72,7 +72,7 @@ sx=int(hdulist[0].header['NAXIS1'])
 sy=int(hdulist[0].header['NAXIS2'])
 
 # fill AITOF map indexes
-# Already done in 02_grxe_resid.py
+# -- Already done in 02_grxe_resid.py
 """
 ds9x=[]
 ds9y=[]
@@ -136,23 +136,8 @@ for i in range(sx):
         if  (df0.shape[0] < nscw_min):
             continue
 
-        #print("*** *** REV *** ***")
-        #print(df0["REV"])
 
-        # check coordinates
+        sg_mean,sg_sem,skew_val,skew_err = get_spec(df0, sigma=4,  grxe_err_cut=grxe_err_cut, enkey=enkey, plotme=False, bootstrap=False, gaussfit=True)
-        #print("***",i+1,j+1,lon,lat,smap[j][i])
-        #for i0,row0 in df0.iterrows():
-        #    print(row0['LON'],row0['LAT'],row0['GRXE'])
-        
-
-        sg_mean,sg_sem = get_spec(df0, sigma=sigma,  grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
-        
-        nsel = int(df0.shape[0]*simfrac/100)
-        #print("nsel=",nsel,df0.shape[0])
-        for n in range(nsim):
-            df1=df0.sample(nsel)
-            sg_mean1,sg_sem1 = get_spec(df1, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
-            mean_sim[dkey].append(sg_mean1)   
 
         
         #print('sg_sem',sg_sem)
@@ -161,19 +146,9 @@ for i in range(sx):
         sign_map[j][i] = sg_mean/sg_sem
         cnt_map[j][i]  = df0.shape[0]
 
-        """
-        obsid_map[dkey]  = obsid
-        grxe_map[dkey]  = grxe
-        grxe_err_map[dkey]  = grxe_err
-        """
-
-
-
-        
 """ Filter by error map """
 # 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(key,enkey,sem_cut,perc))
 
@@ -184,10 +159,6 @@ sem_map[idx]=0.0
 cnt_map[idx]=0
 sign_map[idx]=0.0
 
-#mean_sim_map[idx]=0.0
-#error_sim_map[idx]=0.0
-
-
 if not os.path.exists(mapsdir):
     os.makedirs(mapsdir)
 
@@ -203,76 +174,8 @@ 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)
 
-print("saving simulations")
-for i in range(sx):
-    for j in range(sy):
-        dkey="{:04d}{:04d}".format(j,i)
-        data=mean_sim[dkey]
-        #print("{} size {}".format(dkey,len(data)))
-        if(len(data)>10):
-            (mu, sg) = norm.fit(data)
-            mean_sim_map[j][i] = mu
-            error_sim_map[j][i] = sg
-    
-perc = np.percentile(error_sim_map, sem_cut, axis=(0, 1), keepdims=False)
-print("{} {}: {}% cut of SEM map: {:.2f} mCrab".format(key,enkey,sem_cut,perc))
-
-idx=np.where(error_sim_map > perc)
-print("index size {}".format(len(idx)))
-mean_sim_map[idx]=0.0
-error_sim_map[idx]=0.0
-
-hdulist[0].data=mean_sim_map
-hdulist.writeto(mapsdir+fn.replace(".fits",".sim.mean.fits"),overwrite=True)
-hdulist[0].data=error_sim_map
-hdulist.writeto(mapsdir+fn.replace(".fits",".sim.error.fits"),overwrite=True)
-
-sys.exit()
-
-print("Prepare data for fine map")
-
-obsid_list=[]
-grxe_list=[]
-grxe_err_list=[]
-for i in range(sx):
-    for j in range(sy):
-        """ Use mask """
-        if not (cnt_map[j][i]>0):
-            continue
-        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
-
-        dkey="{:04d}{:04d}".format(j,i)
-        for scw in obsid_map[dkey]:
-            obsid_list.append(scw.decode("UTF-8"))
-        for grxe in grxe_map[dkey]:
-            grxe_list.append(grxe)
-        for grxe in grxe_err_map[dkey]:
-            grxe_err_list.append(grxe)
 
 
-coldefs = fits.ColDefs([
-    fits.Column(name='OBSID', format='11A', array=obsid_list), 
-])
 
-fout = fn.replace(".fits",".grxe.fits")
-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() 

scripts/03_grxe_map.sh

@@ -1,37 +1,10 @@
-./03_grxe_map.py E01
+for band in A01 A02 A03
-#./03_grxe_map.py E02
+do
-#./03_grxe_map.py E03
+./03_grxe_map.py $band 
-#./03_grxe_map.py E04
+done
-#./03_grxe_map.py E05
-#./03_grxe_map.py E06
-#./03_grxe_map.py E07
-#./03_grxe_map.py E08
-#./03_grxe_map.py E09
-#./03_grxe_map.py E10
-./03_grxe_map.py E13
-./03_grxe_map.py E14
 
-
+# if somebody needs maps in narrow nands:
-
+#for band in B01 B02 B03 B04 B05 B06 B07 B08 B09 B10 B11 B12 B13 B14 B15 B16 B17 B18 B19 B20 B21
-
+#do
-./03_grxe_map.py B01
+#./03_grxe_map.py $band 
-./03_grxe_map.py B02
+#done
-./03_grxe_map.py B03
-./03_grxe_map.py B04
-./03_grxe_map.py B05
-./03_grxe_map.py B06
-./03_grxe_map.py B07
-./03_grxe_map.py B08
-./03_grxe_map.py B09
-./03_grxe_map.py B10
-./03_grxe_map.py B11
-./03_grxe_map.py B12
-./03_grxe_map.py B13
-./03_grxe_map.py B14
-./03_grxe_map.py B15
-./03_grxe_map.py B16
-./03_grxe_map.py B17
-./03_grxe_map.py B18
-./03_grxe_map.py B19
-./03_grxe_map.py B20
-./03_grxe_map.py B21

scripts/03_grxe_spec.py

@@ -38,7 +38,7 @@ from ridge.config import *
 inkey="ALL"
 
 
-sigma=3
+
 plotme=False
 
 """
@@ -117,11 +117,15 @@ else:
 if not os.path.exists(specdir):
     os.makedirs(specdir)
 
+fitsdir = "{}fits/".format(specdir)
+if not os.path.exists(fitsdir):
+    os.makedirs(fitsdir)
+
 with open(ignored_rev_file, 'rb') as fp:
     ignored_rev = pickle.load(fp)
     print(ignored_rev)
     print("{} orbits ignored".format(len(ignored_rev)))
-#sys.exit()
+
 ign=ignored_rev.tolist()
 
 
@@ -131,14 +135,7 @@ for skey in skeys:
         print("{} not found in {}".format(skey,list(skyreg.keys())))
         sys.exit()
 
-    # generate array for bootstrap
-    ebands_sim={}
     for enkey in ebands0.keys():
-        ebands_sim[enkey] = []
-
-    for enkey in ebands0.keys():
-        #bkg_fn="detcnts.{}.BKG.resid.fits".format(enkey,inkey)
-        #syserr, bkg_sem = get_syserror(proddir+bkg_fn)
         
         fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
         #d1 = fits.getdata(proddir+fn)
@@ -163,7 +160,7 @@ for skey in skeys:
             skyreg[skey]['lat'] - skyreg[skey]['wlat']/2,
             skyreg[skey]['lat'] + skyreg[skey]['wlat']/2)
 
-        #sys.exit()
+    
         df = df.query(query)
         print("{}, {}: {} N={}".format(skey, enkey, query, df.shape[0]))
 
@@ -178,22 +175,12 @@ for skey in skeys:
         if not (df.shape[0]>0):
             continue
         
-        #plt.scatter(df['LON'],df['LAT'])
+        sg_mean,sg_sem,skew_val,skew_err = get_spec(df, sigma=3,  grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=False, bootstrap=False, gaussfit=True)
-        #plt.show()
-        #print("*** {} {} Data Frame size {} ***".format(skey, enkey, df.size))
-        sg_mean,sg_sem,skew_val,skew_err = get_spec(df, sigma=3,  grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, bootstrap=False, gaussfit=True)
         ebands0[enkey]=[sg_mean,sg_sem]   
         skew0[enkey]=[skew_val,skew_err]   
 
-        nsel = int(df.shape[0]*simfrac/100)
-        for n in range(nsim):
-            df0=df.sample(nsel)
-            #sg_mean,sg_sem = get_spec(df0, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey)
-            #ebands_sim[enkey].append(sg_mean)   
-            #ebands_sim[enkey][1].append(sg_sem)   
         
 
-    ###
     fspec="{}{}.spec".format(specdir,skey)
     with open(fspec, 'w') as fp:
         for enkey in ebands0.keys():
@@ -203,8 +190,6 @@ for skey in skeys:
             fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(bands[enkey],flux,err))    
     subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec])
 
-    ###
-
     fspec="{}{}.skew".format(specdir,skey)
     with open(fspec, 'w') as fp:
         fp.write("read serr 4\n")
@@ -213,37 +198,6 @@ for skey in skeys:
             fp.write("{} {} {:.6f} {:.6f}\n".format(count,bands[enkey],skew0[enkey][0],skew0[enkey][1]))    
             count+=1
     
-    fspec="{}{}.sim.spec".format(specdir,skey)
-    with open(fspec, 'w') as fp:
-        for enkey in ebands_sim.keys():
-            data=ebands_sim[enkey]
-
-            (mu, sg) = norm.fit(data)
-            fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(bands[enkey],mu,sg))    
-            print("[BOOT] {}: {} {:.6f} {:.6f}".format(skey,enkey,mu,sg))
-
-            if(plotme):
-                n, bins, patches = plt.hist(data, 60, density=True, facecolor='green', alpha=0.75)
-                # add a 'best fit' line
-                y = norm.pdf(bins, mu, sg)
-                l = plt.plot(bins, y, 'r--', linewidth=2)
-                #plot
-                plt.axvline(mu, color="black")
-                plt.axvline(ebands0[enkey][0], color="black", linestyle="dashed")
-                #plt.axvline(mu+sg_sem, color="black", linestyle="dashed")
-                #plt.axvline(mu-sg_sem, color="black", linestyle="dashed")
-                plt.axvline(mu+sg, color="blue", linestyle="dashed")
-                plt.axvline(mu-sg, color="blue", linestyle="dashed")
-
-                plt.xlabel('Flux, mCrab')
-                plt.ylabel('Probability')
-                plt.title("[BOOT] {}: {:.2f} {:.2f}".format(enkey, mu, sg))
-                plt.grid(True)
-                plt.show()
-            
-    subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec])
-
-
     try:
         for remfile in ["cols","cols1","cols2","header",]:
             os.remove(remfile)

scripts/README.md

@@ -20,7 +20,7 @@ Plot long-term Crab detector count rate approximated by cubic polynomial functio
 
 ### 01_crabmodel_plot_sys.py
 
-Plots the normalized distribution of the residuals between the IBIS/ISGRI Crab count rate and the corresponding polynomial fit (Fig. B8 in paper). The distribution is approximated with a Gaussian function. 
+Plots the normalized distribution of the residuals between the IBIS/ISGRI Crab count rate and the corresponding polynomial fit (Fig. B8 in the paper). The distribution is approximated with a Gaussian function. 
 
 ### 02_grxe_resid.py
 
@@ -28,21 +28,46 @@ Calculates difference between detector count rate and that predicted by backgrou
 
 ### 02_grxe_resid_plot.py
 
-Plots normalized distribution of the relative residuals of the background model obtained in three energy bands for BKG region (Figs. A5 and A6 in paper).
+Plots normalized distribution of the relative residuals of the background model obtained in three energy bands for BKG region (Figs. A5 and A6 in the paper).
 
 ### 02_grxe_map.py
 
-Makes the map of the residuals in mCrab units (not covered in the paper).
+Makes the map of the residuals in mCrab units (not shown in the paper).
+
+### 03_grxe_flux.py 
+
+Estimates flux in regions defined in config.py, and places result in products/Flux, e.g.:
+
+```
+cat ../products/Flux/GB.dat
+GB A01  25.00   60.00 154.886019 3.087652
+GB A02  60.00   80.00 96.472798 9.218161
+GB A03  80.00  200.00 102.207302 13.291950
+```
+wgere GB means Galactic Bulge region (see Table 1 in the paper).
+
+### 03_grxe_galprof.py/03_grxe_galprof.sh and 03_grxe_galprof_plot.py
+
+Makes Galactic longitude profile and plots it.
 
 ### 03_grxe_spec.py
 
 Extract spectrum of selected sky region, e.g.:
 
-```./03_grxe_spec.py GC``` or  ```./03_grxe_spec.py M81```
+```./03_grxe_spec.py GB``` or  ```./03_grxe_spec.py L+20```
 
 If you do not specify the sky region, the script makes spectra for all regions.
 
 Sky regions are defined in ```config.py```.
 
-Spectra are stored in `./products/Spec`. After spectra files are ready, go to directory `./products/Spec` and run ```perl ../../scripts/do_pha_igr_v3_mCrab.pl M81.spec```.
+Spectra are stored in `../products/Spectra`.
 
+To load the spectrum, e.g. Galactic Bulge, do the following steps:
+
+```
+cd ../products/Spectra
+cp ../../models/xspec/* .
+xspec
+@load_gb_cutoffpl.xcm
+fit
+```