From 9841cf5ded9aa4e4b46d2ff59358531b01f13f7b Mon Sep 17 00:00:00 2001 From: Roman Krivonos Date: Wed, 17 Apr 2024 15:48:08 +0300 Subject: [PATCH] spectra --- ridge/ridge/config.py | 19 ++++- scripts/03_grxe_flux.py | 7 +- scripts/03_grxe_spec.py | 124 ++++++++++++++++++--------------- scripts/README.md | 2 + scripts/do_pha_igr_v3_mCrab.pl | 7 +- 5 files changed, 99 insertions(+), 60 deletions(-) diff --git a/ridge/ridge/config.py b/ridge/ridge/config.py index 0bbbc88..74d3081 100644 --- a/ridge/ridge/config.py +++ b/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}, } diff --git a/scripts/03_grxe_flux.py b/scripts/03_grxe_flux.py index 335a856..881712e 100755 --- a/scripts/03_grxe_flux.py +++ b/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 diff --git a/scripts/03_grxe_spec.py b/scripts/03_grxe_spec.py index cc18745..6e8c16c 100755 --- a/scripts/03_grxe_spec.py +++ b/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 diff --git a/scripts/README.md b/scripts/README.md index d23d458..1ee16ab 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -28,3 +28,5 @@ source /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```. + +Если не задать регион неба, скрипт делает спектры для всех регионов. \ No newline at end of file diff --git a/scripts/do_pha_igr_v3_mCrab.pl b/scripts/do_pha_igr_v3_mCrab.pl index 2961a1e..871d9f6 100755 --- a/scripts/do_pha_igr_v3_mCrab.pl +++ b/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 <