spectra

ridge/ridge/config.py

@@ -70,8 +70,23 @@ basemin=20
 skyreg={
     'M81':{'lon':142.091816,'lat':40.900054,'wlon':20.0,'wlat':20.0},
     'BKG':{'lon':-70.0,'lat':70.0,'wlon':40.0,'wlat':40.0},
-    'GC':{'lon':0.0,'lat':0.0,'wlon':10.0,'wlat':10.0},
+    'GC10':{'lon':0.0,'lat':0.0,'wlon':10.0,'wlat':10.0},
+    'GC06':{'lon':0.0,'lat':0.0,'wlon':6.0,'wlat':6.0},
+    'LAT+20':{'lon':-34.0,'lat':20.0,'wlon':10.0,'wlat':10.0},
+    'CenA':{'lon':-50,'lat':20.0,'wlon':15.0,'wlat':15.0},
+    'CygX1':{'lon':71.0,'lat':1.0,'wlon':15.0,'wlat':20.0},
+    'ScoX1':{'lon':0.0,'lat':23.0,'wlon':15.0,'wlat':15.0},
+    'LON-75':{'lon':-75.0,'lat':0.0,'wlon':20.0,'wlat':20.0},
+    'LON+115':{'lon':115.0,'lat':0.0,'wlon':20.0,'wlat':10.0},
+    'LON+40':{'lon':43.0,'lat':0.0,'wlon':20.0,'wlat':10.0},
     'LON+20':{'lon':20.0,'lat':0.0,'wlon':20.0,'wlat':10.0},
     'LON-20':{'lon':-20.0,'lat':0.0,'wlon':20.0,'wlat':10.0},
-    'Geminga':{'lon':195.133837,'lat':4.265811,'wlon':10.0,'wlat':10.0},
+    'GP':{'lon':0.0,'lat':0.0,'wlon':100.0,'wlat':10.0},
+    'HerX1':{'lon':58.1500665,'lat':37.5219674,'wlon':10.0,'wlat':10.0},
+    '3C273':{'lon':-71,'lat':69.0,'wlon':20.0,'wlat':20.0},
+    'SMC':{'lon':-59,'lat':-44,'wlon':20.0,'wlat':20.0},
+    'LMC':{'lon':-80,'lat':-32,'wlon':20.0,'wlat':20.0},
+    'LAT-30':{'lon':0.0,'lat':-30,'wlon':10.0,'wlat':10.0},
+    '3C120':{'lon':-159,'lat':-21,'wlon':25.0,'wlat':25.0},
+    'Geminga':{'lon':-164.866,'lat':4.265811,'wlon':10.0,'wlat':10.0},
 }

scripts/03_grxe_flux.py

@@ -42,7 +42,12 @@ d = fits.getdata(proddir+fn)
 df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
 #print(df.columns)
 
-key="GC"
+#key="GC"
+#sz=5
+#lon0=0.0
+#lat0=0.0
+
+key="BKG"
 sz=5
 lon0=0.0
 lat0=0.0

scripts/03_grxe_spec.py

@@ -26,6 +26,7 @@ from scipy.stats import norm
 from scipy.stats import describe
 from scipy.stats import sem
 
+import subprocess
 
 from numpy import absolute
 from numpy import arange
@@ -52,64 +53,77 @@ ebands0={#'E02':[0.0,0.0],
          'E12':[0.0,0.0],}
 
 #skey='Geminga'
-skey = sys.argv[1]
+if len(sys.argv) > 1:
+    skeys = [sys.argv[1]]
+else:
+    skeys = list(skyreg.keys())
 
-if not skey in skyreg.keys():
-    print("{} not found in {}".format(skey,list(skyreg.keys())))
-    sys.exit()
 
-for enkey in ebands0.keys():
-    fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
-    d = fits.getdata(proddir+fn)
-    df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
-    #print(df.columns)
-
-    df = df.query('LON > {} & LON < {} & LAT > {} & LAT < {}'.format(skyreg[skey]['lon'] - skyreg[skey]['wlon']/2,
-                                                                     skyreg[skey]['lon'] + skyreg[skey]['wlon']/2,
-                                                                     skyreg[skey]['lat'] - skyreg[skey]['wlat']/2,
-                                                                     skyreg[skey]['lat'] + skyreg[skey]['wlat']/2))
-
-    #print("Number of ScWs: {}".format(df.shape[0]))
-    if not (df.shape[0]>0):
-        continue
-    #plt.scatter(df['LON'],df['LAT'])
-    #plt.show()
-
-    grxe = np.array(df['GRXE'])
-
-    filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=10, return_bounds=True)
-    filtered_grxe = filtered_data[0]
-    filtered_min  = filtered_data[1]
-    filtered_max  = filtered_data[2]
-
-    #print("length orig: {} taken: {} filtered: {}".format(len(grxe),len(grxe[filtered_grxe.mask==False]),len(grxe[filtered_grxe.mask==True])))
-
-    sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=10)
-    sg_sem = sem(grxe[filtered_grxe.mask==False])
-    print("{}: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f} N={}".format(enkey, sg_mean, sg_med, sg_std, sg_sem, len(grxe[filtered_grxe.mask==False])))
-    sg_sem*=1.5
-    if(sg_mean<0.0):
-        sg_mean=1e-9
-        sg_sem*=2
-
-    ebands0[enkey]=[sg_mean,sg_sem]   
-
-    if(plotme):
-        k=1.2
-        plt.hist(grxe, bins=n_bins, range=[filtered_min*k, filtered_max*k])
-        plt.hist(grxe[filtered_grxe.mask], bins=n_bins, range=[filtered_min*k, filtered_max*k])
-        plt.axvline(sg_mean, color="black")
-        plt.axvline(sg_mean+sg_sem, color="black", linestyle="dashed")
-        plt.axvline(sg_mean-sg_sem, color="black", linestyle="dashed")
-        plt.axvline(sg_mean+sg_std, color="blue", linestyle="dashed")
-        plt.axvline(sg_mean-sg_std, color="blue", linestyle="dashed")
-        plt.xlabel("mCrab")
-        plt.show()
-
-        
 if not os.path.exists(specdir):
     os.makedirs(specdir)
 
-with open("{}{}.spec".format(specdir,skey), 'w') as fp:
+
+for skey in skeys:
+    if not skey in skyreg.keys():
+        print("{} not found in {}".format(skey,list(skyreg.keys())))
+        sys.exit()
     for enkey in ebands0.keys():
-        fp.write("0 {} {:.2f} {:.2f} 0.0\n".format(bands[enkey],ebands0[enkey][0],ebands0[enkey][1]))
+        fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
+        d = fits.getdata(proddir+fn)
+        df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
+        #print(df.columns)
+
+        df = df.query('LON > {} & LON < {} & LAT > {} & LAT < {}'.format(skyreg[skey]['lon'] - skyreg[skey]['wlon']/2,
+                                                                         skyreg[skey]['lon'] + skyreg[skey]['wlon']/2,
+                                                                         skyreg[skey]['lat'] - skyreg[skey]['wlat']/2,
+                                                                         skyreg[skey]['lat'] + skyreg[skey]['wlat']/2))
+
+        #print("Number of ScWs: {}".format(df.shape[0]))
+        if not (df.shape[0]>0):
+            continue
+        #plt.scatter(df['LON'],df['LAT'])
+        #plt.show()
+
+        grxe = np.array(df['GRXE'])
+
+        filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=10, return_bounds=True)
+        filtered_grxe = filtered_data[0]
+        filtered_min  = filtered_data[1]
+        filtered_max  = filtered_data[2]
+
+        #print("length orig: {} taken: {} filtered: {}".format(len(grxe),len(grxe[filtered_grxe.mask==False]),len(grxe[filtered_grxe.mask==True])))
+
+        sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=10)
+        sg_sem = sem(grxe[filtered_grxe.mask==False])
+        print("{}: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f} N={}".format(enkey, sg_mean, sg_med, sg_std, sg_sem, len(grxe[filtered_grxe.mask==False])))
+
+        #sg_sem*=1.5
+        #if(sg_mean<0.0):
+        #    sg_mean=1e-9
+        #    sg_sem*=2
+
+        ebands0[enkey]=[sg_mean,sg_sem]   
+
+        if(plotme):
+            k=1.2
+            plt.hist(grxe, bins=n_bins, range=[filtered_min*k, filtered_max*k])
+            plt.hist(grxe[filtered_grxe.mask], bins=n_bins, range=[filtered_min*k, filtered_max*k])
+            plt.axvline(sg_mean, color="black")
+            plt.axvline(sg_mean+sg_sem, color="black", linestyle="dashed")
+            plt.axvline(sg_mean-sg_sem, color="black", linestyle="dashed")
+            plt.axvline(sg_mean+sg_std, color="blue", linestyle="dashed")
+            plt.axvline(sg_mean-sg_std, color="blue", linestyle="dashed")
+            plt.xlabel("mCrab")
+            plt.show()
+
+    fspec="{}{}.spec".format(specdir,skey)
+    with open(fspec, 'w') as fp:
+        for enkey in ebands0.keys():
+            fp.write("0 {} {:.2f} {:.2f} 0.0\n".format(bands[enkey],ebands0[enkey][0],ebands0[enkey][1]))
+        
+    subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec])
+    try:
+        for remfile in ["cols","cols1","cols2","header"]:
+            os.remove(remfile)
+    except OSError:
+        pass

scripts/README.md

@@ -28,3 +28,5 @@ source <MY PATH>/venv/bin/activate.csh
 ```./03_grxe_spec.py GC``` или  ```./03_grxe_spec.py M81```
 
 Регионы неба задаются в ```config.py```. Спектры помещаются в директорию ```./products/spec```. Затем надо перейти в эту директорию и выполнить, например: ```perl ../../scripts/do_pha_igr_v3_mCrab.pl M81.spec```.
+
+Если не задать регион неба, скрипт делает спектры для всех регионов.

scripts/do_pha_igr_v3_mCrab.pl

@@ -1,7 +1,7 @@
 #! /usr/bin/perl
 #require "utils.pl";
 #require "getopts.pl";
-
+use File::Basename;
 #&Getopts('i:o:n:c:e');
 
 #############################################
@@ -138,6 +138,9 @@ print COLS "STAT_ERR                   1D                  counts\n";
 
 close(COLS);
 
+
+$base = basename($rmff);
+
 open(HEADER,">header");
 print HEADER <<EOHEAD;
 EXTNAME = 'SPECTRUM'           / name of this binary table extension
@@ -159,7 +162,7 @@ BACKSCAL=       1.00000000E+00 / Background scale factor
 CORRSCAL=       0.00000000E+00 / Correlation scale factor
 BACKFILE= 'NONE    '           / Background FITS file for this object
 CORRFILE= 'NONE    '           / Correlation FITS file for this object
-RESPFILE= '$rmff'         / Redistribution matrix file (RMF)
+RESPFILE= '$base'         / Redistribution matrix file (RMF)
 ANCRFILE= 'NONE    '           / Ancillary response file (ARF)
 XFLT0001= 'NONE    '           / XSPEC selection filter description
 CHANTYPE= 'PHA     '           / Channels assigned by detector electronics