diff --git a/data/xspec/GC06.xcm b/data/xspec/GC06.xcm new file mode 100644 index 0000000..4ab05cc --- /dev/null +++ b/data/xspec/GC06.xcm @@ -0,0 +1,12 @@ +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 diff --git a/data/xspec/GC10.xcm b/data/xspec/GC10.xcm new file mode 100644 index 0000000..b1c02c9 --- /dev/null +++ b/data/xspec/GC10.xcm @@ -0,0 +1,12 @@ +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 diff --git a/data/xspec/LON+20.xcm b/data/xspec/LON+20.xcm new file mode 100644 index 0000000..8706b9a --- /dev/null +++ b/data/xspec/LON+20.xcm @@ -0,0 +1,10 @@ +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 diff --git a/ridge/ridge/config.py b/ridge/ridge/config.py index 74d3081..654daae 100644 --- a/ridge/ridge/config.py +++ b/ridge/ridge/config.py @@ -45,6 +45,11 @@ SEM означает standatd error on mean (~RMS/sqrt(n)) """ sem_cut=98 +""" +Измерение GRXE ((data-model)/crab) со значением ошибки выше этого персентиля будет отброшено +""" +grxe_err_cut=90 + """ Координаты Краба """ crab_ra = 83.6287 crab_dec = 22.014 diff --git a/scripts/02_grxe_resid.py b/scripts/02_grxe_resid.py index 2588d41..0c50b84 100755 --- a/scripts/02_grxe_resid.py +++ b/scripts/02_grxe_resid.py @@ -43,8 +43,9 @@ with open(proddir+fn.replace(".fits",".ignored_rev.pkl"), 'rb') as fp: with open(proddir+fn.replace(".fits",".crabmodel.pkl"), 'rb') as fp: crabmodel, z = pickle.load(fp) p = np.poly1d(z) + crabmodel_keys = list(crabmodel.keys()) #print(crabmodel) - +#sys.exit() crab_rev_max = np.max(list(crabmodel.keys())) print("Crab is defined untill orbit {}".format(crab_rev_max)) @@ -66,11 +67,13 @@ clean0=[] model0=[] resid0=[] # residuals in cts/s grxe0=[] # mCrab +grxe_err0=[] # mCrab crab0=[] # Crab count rate mjd0=[] a0=[] b0=[] -err0=[] +model_err0=[] +crab_err0=[] lon0=[] lat0=[] base0=[] @@ -114,11 +117,22 @@ for i, row in df.iterrows(): a = bgdmodel[orbit]['a'] b = bgdmodel[orbit]['b'] c = bgdmodel[orbit]['c'] + + # Crab error is defined only for Crab resolutions, so we interpolate between + if(orbit in crabmodel_keys): + crab_err = crabmodel[orbit]['err'] + else: + left,right = find_nearest(crabmodel_keys, orbit) + crab_err = np.interp(orbit, [left,right], [crabmodel[left]['err'], crabmodel[right]['err']]) + #print() + #print(orbit, left, right) + #print(orbit, 'err', crabmodel[left]['err'], crab_err, crabmodel[right]['err']) + #sys.exit() err = bgdmodel[orbit]['err'] m = a*row['PHASE']+b r1 = bgdmodel[orbit]['r1'] # nearest left orbit used for calibration r2 = bgdmodel[orbit]['r2'] # nearest right orbit used for calibration - + c0.append(c) base0.append(abs(orbit - int(np.min([r1,r2])))) clean0.append(clean[i]) @@ -126,12 +140,14 @@ for i, row in df.iterrows(): model0.append(m) resid0.append(clean[i]-m) grxe0.append(1000*(clean[i]-m)/p(orbit)) + grxe_err0.append(1000*np.sqrt(err**2 + crab_err**2)/p(orbit)) crab0.append(p(orbit)) a0.append(a) b0.append(b) c0.append(c) - err0.append(err) + model_err0.append(err) + crab_err0.append(crab_err) phase0.append(row['PHASE']) rev0.append(orbit) lon0.append(row['LON']) @@ -155,13 +171,15 @@ coldefs = fits.ColDefs([ 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='MODEL_ERR', format='D', unit='', array=[model_err0[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='GRXE_ERR', format='D', unit='mCrab', array=[grxe_err0[index] for index in indices]), fits.Column(name='CRAB', format='D', unit='cts/s', array=[crab0[index] for index in indices]), + fits.Column(name='CRAB_ERR', format='D', unit='', array=[crab_err0[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='C', format='D', unit='', array=[c0[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)) diff --git a/scripts/03_grxe_spec.py b/scripts/03_grxe_spec.py index 6e8c16c..370d856 100755 --- a/scripts/03_grxe_spec.py +++ b/scripts/03_grxe_spec.py @@ -85,6 +85,13 @@ for skey in skeys: #plt.show() grxe = np.array(df['GRXE']) + grxe_err = np.array(df['GRXE_ERR']) + + perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False) + print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(skey,enkey,sem_cut,perc)) + idx=np.where(grxe_err < perc) + grxe=grxe[idx] + grxe_err=grxe_err[idx] filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=10, return_bounds=True) filtered_grxe = filtered_data[0] @@ -99,7 +106,7 @@ for skey in skeys: #sg_sem*=1.5 #if(sg_mean<0.0): - # sg_mean=1e-9 + # sg_mean=1e-6 # sg_sem*=2 ebands0[enkey]=[sg_mean,sg_sem] @@ -119,7 +126,7 @@ for skey in skeys: 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])) + fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(bands[enkey],ebands0[enkey][0],ebands0[enkey][1])) subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec]) try: diff --git a/scripts/03_plot_skyreg.py b/scripts/03_plot_skyreg.py index 1223493..3537a8a 100755 --- a/scripts/03_plot_skyreg.py +++ b/scripts/03_plot_skyreg.py @@ -71,6 +71,7 @@ plt.show() rev = np.array(df['REV']) grxe = np.array(df['GRXE']) +grxe_err = np.array(df['GRXE_ERR']) @@ -86,6 +87,35 @@ sstr = '%-14s mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f' n, (smin, smax), sm, sv, ss, sk = describe(grxe) print(sstr % ('sample:', sm, sv, ss, sk)) +# Calculate the percentiles across the x dimension +errmax=np.max(grxe_err) +perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False) +print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(skey,enkey,sem_cut,perc)) +idx=np.where(grxe_err < perc) + +plt.hist(grxe_err, bins=n_bins, range=[0,errmax], color="red") +rev=rev[idx] +grxe=grxe[idx] +grxe_err=grxe_err[idx] +plt.hist(grxe_err, bins=n_bins, range=[0,errmax], color="grey") + +plt.xlabel("GRXE ERROR, mCrab") +plt.title("Distribution of errors {}".format(skey)) +plt.show() + +grxe_sign=np.divide(grxe,grxe_err) +plt.hist(grxe_sign, bins=n_bins) +plt.xlabel("GRXE S/N") +plt.title("Distribution of significance {}".format(skey)) +plt.show() + +A=np.sum(grxe/grxe_err**2) +B=np.sum(1.0/grxe_err**2) +wgt_mean=A/B +wgt_mean_err=np.sqrt(1.0/B) + +print("Weighted mean: {:.2f}+/-{:.2f} normal mean: {:.2f}".format(wgt_mean,wgt_mean_err,np.mean(grxe))) + #filtered_grxe = sigma_clip(grxe, sigma=sigma, maxiters=10) filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=10, return_bounds=True) filtered_grxe = filtered_data[0]