diff --git a/data/polarmodel.fits b/data/polarmodel.fits new file mode 100644 index 0000000..056f5d0 Binary files /dev/null and b/data/polarmodel.fits differ diff --git a/ridge/ridge/config.py b/ridge/ridge/config.py index b641462..641008b 100644 --- a/ridge/ridge/config.py +++ b/ridge/ridge/config.py @@ -43,3 +43,17 @@ sem_cut=98 """ Координаты Краба """ crab_ra = 83.6287 crab_dec = 22.014 + +bands={ + 'E02':'17.00 22.10', + 'E03':'22.10 28.73', + 'E04':'28.73 37.35', + 'E05':'37.35 48.55', + 'E06':'48.55 63.12', + 'E07':'63.12 82.06', + 'E08':'82.06 106.67', + 'E09':'106.67 138.67', + 'E10':'138.67 180.28', + 'E11':'180.28 234.36', + 'E12':'234.36 304.67', +} diff --git a/scripts/03_grxe_flux.py b/scripts/03_grxe_flux.py new file mode 100755 index 0000000..fcb0596 --- /dev/null +++ b/scripts/03_grxe_flux.py @@ -0,0 +1,95 @@ +#!/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 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 * + +key="GAL" +enkey = sys.argv[1] +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(3,3)) + +sz=3 +lon0=0.0 +lat0=0.0 +df = df.query('LON > {} & LON < {} & LAT > {} & LAT < {}'.format(lon0-sz,lon0+sz,lat0-sz,lat0+sz)) + +print("Number of ScWs: {}".format(df.shape[0])) + +n_bins = 80 +sigma=3 + +grxe = np.array(df['GRXE']) + +mean = np.mean(grxe) +std = np.std(grxe) +print("\n*** Unfiltered:") +print("mean {:.2f} std {:.2f}".format(mean,std)) +print("min {:.2f}".format(grxe.min())) +print("max {:.2f}".format(grxe.max())) +print("mean {:.2f}".format(grxe.mean())) +print("median {:.2f}".format(np.median(grxe))) +print("var {:.2f}".format(grxe.var())) +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)) + +print("\n***") + +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("Sigma clipping: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f}".format(sg_mean, sg_med, sg_std, sg_sem)) + + +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() diff --git a/scripts/03_grxe_map.py b/scripts/03_grxe_map.py index 837ec89..30ee53f 100755 --- a/scripts/03_grxe_map.py +++ b/scripts/03_grxe_map.py @@ -131,7 +131,7 @@ for i in range(sx): # 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)) +print("{} {}: {}% cut of SEM map: {:.2f} mCrab".format(key,enkey,sem_cut,perc)) idx=np.where(sem_map > perc) mean_map[idx]=0.0 diff --git a/scripts/03_grxe_map.sh b/scripts/03_grxe_map.sh index b3b7863..face240 100644 --- a/scripts/03_grxe_map.sh +++ b/scripts/03_grxe_map.sh @@ -9,3 +9,28 @@ ./03_grxe_map.py E10 ALL ./03_grxe_map.py E11 ALL ./03_grxe_map.py E12 ALL + +./03_grxe_map.py E02 GAL +./03_grxe_map.py E03 GAL +./03_grxe_map.py E04 GAL +./03_grxe_map.py E05 GAL +./03_grxe_map.py E06 GAL +./03_grxe_map.py E07 GAL +./03_grxe_map.py E08 GAL +./03_grxe_map.py E09 GAL +./03_grxe_map.py E10 GAL +./03_grxe_map.py E11 GAL +./03_grxe_map.py E12 GAL + +./03_grxe_map.py E02 BKG +./03_grxe_map.py E03 BKG +./03_grxe_map.py E04 BKG +./03_grxe_map.py E05 BKG +./03_grxe_map.py E06 BKG +./03_grxe_map.py E07 BKG +./03_grxe_map.py E08 BKG +./03_grxe_map.py E09 BKG +./03_grxe_map.py E10 BKG +./03_grxe_map.py E11 BKG +./03_grxe_map.py E12 BKG + diff --git a/scripts/03_grxe_spec.py b/scripts/03_grxe_spec.py new file mode 100755 index 0000000..e1e3706 --- /dev/null +++ b/scripts/03_grxe_spec.py @@ -0,0 +1,103 @@ +#!/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 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 * + +key="GAL" +n_bins = 80 +sigma=3 +plotme=False + +ebands0={#'E02':[0.0,0.0], + 'E03':[0.0,0.0], + 'E04':[0.0,0.0], + 'E05':[0.0,0.0], + 'E06':[0.0,0.0], + 'E07':[0.0,0.0], + 'E08':[0.0,0.0], + 'E09':[0.0,0.0], + 'E10':[0.0,0.0], + 'E11':[0.0,0.0], + 'E12':[0.0,0.0],} + +sz=3 +lon0=0.0 +lat0=0.0 + +for enkey in ebands0.keys(): + 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('LON > {} & LON < {} & LAT > {} & LAT < {}'.format(lon0-sz,lon0+sz,lat0-sz,lat0+sz)) + + #print("Number of ScWs: {}".format(df.shape[0])) + + 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}".format(enkey, sg_mean, sg_med, sg_std, sg_sem)) + 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() + + + +with open(proddir+"GRXE.spec", '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])) diff --git a/scripts/03_plot_bkg.py b/scripts/03_plot_bkg.py new file mode 100755 index 0000000..c10977a --- /dev/null +++ b/scripts/03_plot_bkg.py @@ -0,0 +1,78 @@ +#!/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 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 = "BKG" + + +fn='detcnts.{}.{}.resid.fits'.format(enkey,key) + +n_bins = 80 +minmax=[-1800,1800] +sigma=3 +with fits.open(proddir+fn) as hdul: + data=hdul[1].data + grxe = np.array(data['GRXE']) + mean = np.mean(grxe) + std = np.std(grxe) + print("mean {:.2f} std {:.2f}".format(mean,std)) + print("min {:.2f}".format(grxe.min())) + print("max {:.2f}".format(grxe.max())) + print("mean {:.2f}".format(grxe.mean())) + print("median {:.2f}".format(np.median(grxe))) + print("var {:.2f}".format(grxe.var())) + 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)) + + #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] + filtered_min = filtered_data[1] + filtered_max = filtered_data[2] + + sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=10) + 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)) + 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_std, color="black", linestyle="dashed") + plt.axvline(sg_mean-sg_std, color="black", linestyle="dashed") + plt.xlabel("mCrab") + plt.show() diff --git a/scripts/README.md b/scripts/README.md index 0bd67c5..3d867a6 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -7,15 +7,15 @@ source /venv/bin/activate.csh ### 01_bgdmodel.py / 01_bgdmodel.sh -Создает модель фона +Создает модель фона. ### 01_crabmodel.py / 01_crabmodel.sh -Создает модель скорости счета Краба +Создает модель скорости счета Краба. ### 02_grxe_resid.py / 01_grxe_resid.sh -Вычисляет остатки после вычитания модели из данных в единицах мКраб +Вычисляет остатки после вычитания модели из данных в единицах мКраб. ### 02_grxe_map.py / 01_grxe_map.sh diff --git a/scripts/do_pha_igr_v3_mCrab.pl b/scripts/do_pha_igr_v3_mCrab.pl new file mode 100755 index 0000000..2961a1e --- /dev/null +++ b/scripts/do_pha_igr_v3_mCrab.pl @@ -0,0 +1,186 @@ +#! /usr/bin/perl +#require "utils.pl"; +#require "getopts.pl"; + +#&Getopts('i:o:n:c:e'); + +############################################# +$alpha=2.1; +$norm=10; + + +$mean_rms=0.73; + + +$inp=$ARGV[0]; +$out=$inp; +open INPP,">$inp.dat"; +@spec=`cat $inp`; +chomp @spec; +$cnt1=0; +foreach $sp (@spec){ + ($null,$e1[$cnt1],$e2[$cnt1], $rate[$cnt1], $drate[$cnt1], $null)=split / +/,$sp; + $rate[$cnt1]=$rate[$cnt1]/1000; + $drate[$cnt1]=$drate[$cnt1]/1000; + $crrate[$cnt1]=1; + $dcrrate[$cnt1]=0.001; + print INPP "$cnt1 $rate[$cnt1] $drate[$cnt1]\n"; + $cnt1++; +} +close(INPP); + + + +$TFORM2=sprintf "%dE",$cnt1; +$rmff="$out.rmf"; +open COLS,">cols1"; +print COLS "ENERG_LO 1e keV +ENERG_HI 1e keV +N_GRP I +F_CHAN 1I +N_CHAN 1I +MATRIX $TFORM2 +"; +close COLS; +open COLS,">cols2"; +print COLS "CHANNEL J +E_MIN E keV +E_MAX E keV +"; +close COLS; + +open FITS, "|fcreate cols1 datafile=\"-\" outfile=\"$rmff\" tbltype=\"binary\" extname=\"SPECRESP MATRIX\" clobber=yes"; + + +$beta=1-$alpha; +for($i=0;$i<$cnt1;){ + $ener=($e1[$i]+$e2[$i])/2; + $crfl[$i]=$norm*(exp($beta*log($e2[$i]))-exp($beta*log($e1[$i])))/$beta; + $rmf[$i]=$crrate[$i]/$crfl[$i]; + print"$i $e1[$i] $e2[$i] $rmf[$i] $crfl[$i]\n"; + $num=$i+1; + print FITS "$e1[$i] $e2[$i] 1 0 $cnt1"; + for($j=0;$j<$cnt1;){ + $koef=0; + if($i==$j){$koef=$rmf[$i];} + print FITS " $koef"; + $j++; + } + + print FITS "\n"; + + $i++; + +} + +close FITS; + +open FITS, "|fcreate cols2 datafile=\"-\" outfile=\"temp.fits\" tbltype=\"binary\" extname=\"EBOUNDS\" clobber=yes"; + +for($i=0;$i<$cnt1;){ + $num=$i+1; + print FITS "$i $e1[$i] $e2[$i]\n"; + + $i++; + +} + +close FITS; + +system "fappend temp.fits[1] $rmff"; +system "rm -f temp.fits"; + +system "fparkey OGIP $rmff\[1] HDUCLASS add=yes"; +system "fparkey RESPONSE $rmff\[1] HDUCLAS1 add=yes"; +system "fparkey RSP_MATRIX $rmff\[1] HDUCLAS2 add=yes"; +system "fparkey $cnt1 $rmff\[1] DETCHANS add=yes"; +system "fparkey PHA $rmff\[1] CHANTYPE add=yes"; +system "fparkey 1.0 $rmff\[1] EFFAREA add=yes"; +system "fparkey FULL_RES $rmff\[1] BINNING add=yes"; +system "fparkey 1992a $rmff\[1] RMFVERSN add=yes"; +system "fparkey INTEGRAL $rmff\[1] TELESCOP add=yes"; +system "fparkey ISGRI $rmff\[1] INSTRUME add=yes"; +system "fparkey DAL_TABLE $rmff\[1] BASETYPE add=yes"; +system "fparkey 1 $rmff\[1] EXTVER add=yes"; +system "fparkey 1 $rmff\[1] GRPID1 add=yes"; +system "fparkey 8 $rmff\[1] TLMAX4 add=yes"; +system "fparkey 0 $rmff\[1] TLMIN4 add=yes"; + +system "fparkey OGIP $rmff\[2] HDUCLASS add=yes"; +system "fparkey RESPONSE $rmff\[2] HDUCLAS1 add=yes"; +system "fparkey EBOUNDS $rmff\[2] HDUCLAS2 add=yes"; +system "fparkey $cnt1 $rmff\[2] DETCHANS add=yes"; +system "fparkey PHA $rmff\[2] CHANTYPE add=yes"; +system "fparkey 1.0 $rmff\[2] EFFAREA add=yes"; +system "fparkey DAL_TABLE $rmff\[2] BASETYPE add=yes"; +system "fparkey 1992a $rmff\[2] RMFVERSN add=yes"; +system "fparkey 1.3.0 $rmff\[2] HDUVERS add=yes"; +system "fparkey 1.0.0 $rmff\[2] HDUVERS1 add=yes"; +system "fparkey 1.3.0 $rmff\[2] HDUVERS2 add=yes"; +system "fparkey 8 $rmff\[2] TLMAX1 add=yes"; +system "fparkey 0 $rmff\[2] TLMIN1 add=yes"; + +system "fparkey 8 $rmff\[0] BITPIX add=yes"; + +#system "cphead crab_1.15.rmf $rmff"; +#system "cphead crab_1.15.rmf\[1] $rmff\[1]"; +#system "cphead crab_1.15.rmf\[2] $rmff\[2]"; +#system "fparkey $cnt1 $rmff\[2] DETCHANS add=yes"; + + +###### pha file ######### + + +open(COLS,">cols"); +print COLS "CHANNEL 1I\n"; +print COLS "RATE 1D counts\n"; +print COLS "STAT_ERR 1D counts\n"; + +close(COLS); + +open(HEADER,">header"); +print HEADER <