clean

2024-10-31 13:56:16 +03:00
parent 09ce08e69a
commit e02db833cc
313 changed files with 126 additions and 99566 deletions
--- a/scripts/01_bgdmodel.py
+++ b/scripts/01_bgdmodel.py
@@ -30,48 +30,27 @@ enkey = sys.argv[1]
 fn="detcnts.{}.fits".format(enkey)
 print("Reading {}".format(datadir+fn))

-#d = fits.getdata(datadir+fn)
-#df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
-
 with fits.open(datadir+fn) as data:
    df = pd.DataFrame(data[1].data)
    
-print(df.columns)
+# print(df.columns)

-#interp_a = np.interp(50, [0,100], [0, 1])
-#print(interp_a)
-#sys.exit()
-        
 plotme=False
-sigma=3
 ntotal=0
 nrev=0
 bgdmodel={}

+# will be populated
 ignored_scw=[]
-ignored_rev=[131,132,133,134,135,136,# flare
-             280,281,282,#flare
-             1057,1102,1150, # flare
-             1452,1459,1472, # flare
-             2303,2304, # flare
-             2549, # flare
-             334,1760] # ignore orbits with pecular slope over phase
-ignored_rev=[131,132,133,136,
-             280,281,282,
-             1057,1102,1150,1452,1459,
-             2303,2304,2305,
-             2483,2549]
+

 if not os.path.exists(proddir):
    os.makedirs(proddir)
    
 for rev in range(revmin,revmax):

-#    if not (rev==341):
-#        continue
-
    df0 = df.query('CLEAN > 0.0 & REV == {} & ( abs(LAT) > {} | abs(LON) > {}) & PHASE > {} & PHASE < {}'.format(rev,bmax,lmax,phmin,phmax))
-    #df0 = df.query('CLEAN > 0.0 & REV == {} & abs(LAT) > {} & abs(LON) > {} & PHASE > {} & PHASE < {}'.format(rev,bmax,lmax,phmin,phmax))
+
    nobs=len(df0)
    print("*** REV {} *** {} ScWs".format(rev,df0.shape[0]))
    if not(nobs):
@@ -89,7 +68,7 @@ for rev in range(revmin,revmax):
            c = 0
            """ run regression """
            print("*** Run regression for {}".format(rev))
-            #print(df0['CLEAN'].values)
+            
            x = x.reshape((-1, 1))
            # https://machinelearningmastery.com/robust-regression-for-machine-learning-in-python/
            #model = LinearRegression()
--- a/scripts/01_bgdmodel.sh
+++ b/scripts/01_bgdmodel.sh
@@ -1,49 +1,10 @@
-./01_bgdmodel.py E01
-#./01_bgdmodel.py E02
-#./01_bgdmodel.py E03
-#./01_bgdmodel.py E04
-#./01_bgdmodel.py E05
-#./01_bgdmodel.py E06
-#./01_bgdmodel.py E07
-#./01_bgdmodel.py E08
-#./01_bgdmodel.py E09
-#./01_bgdmodel.py E10
-#./01_bgdmodel.py E11
-#./01_bgdmodel.py E12
-./01_bgdmodel.py E13
-./01_bgdmodel.py E14
-#./01_bgdmodel.py E15
+for band in A01 A02 A03
+do
+./01_bgdmodel.py $band 
+done

-#./01_bgdmodel.py A01
-#./01_bgdmodel.py A02
-#./01_bgdmodel.py A03
-#./01_bgdmodel.py A04
-#./01_bgdmodel.py A05
-#./01_bgdmodel.py A06
-#./01_bgdmodel.py A07
-#./01_bgdmodel.py A08
-#./01_bgdmodel.py A09
-#./01_bgdmodel.py A10
-#./01_bgdmodel.py A11
+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
+./01_bgdmodel.py $band 
+done

-./01_bgdmodel.py B01
-./01_bgdmodel.py B02
-./01_bgdmodel.py B03
-./01_bgdmodel.py B04
-./01_bgdmodel.py B05
-./01_bgdmodel.py B06
-./01_bgdmodel.py B07
-./01_bgdmodel.py B08
-./01_bgdmodel.py B09
-./01_bgdmodel.py B10
-./01_bgdmodel.py B11
-./01_bgdmodel.py B12
-./01_bgdmodel.py B13
-./01_bgdmodel.py B14
-./01_bgdmodel.py B15
-./01_bgdmodel.py B16
-./01_bgdmodel.py B17
-./01_bgdmodel.py B18
-./01_bgdmodel.py B19
-./01_bgdmodel.py B20
-./01_bgdmodel.py B21
--- a/scripts/01_crabmodel.py
+++ b/scripts/01_crabmodel.py
@@ -41,18 +41,17 @@ enkey = sys.argv[1]
 sigma=3

 # for these bands, slope is taken from stacked profile
-fixed_slope = ['E04','E05','E06','E07','E08','E09','E10','E11','E12','E13','E14','E15',
-               'A02','A03','A04','A05','A06','A07','A08','A09','A10','A11',
+fixed_slope = ['A02','A03',
               'B02','B03','B04','B05','B06','B07','B08','B09','B10','B11','B12','B13','B14','B15','B16','B17','B18','B19','B20','B21']

 # for these bands, slope is free for each revolution
-free_slope = ['E01', 'E02', 'E03', 'A01','B01']
+free_slope = ['A01', 'B01']

 # for these bands, slope is fixed at constant (or if positive, which is not allowed)
-const_slope = ['E10','E11','E12','E13','E14','E15','A10','A11','B18','B19','B20','B21']
+const_slope = ['A02','A03','B18','B19','B20','B21']

 # for stacked profile, skip orbits>800 for energy channels <30 keV
-skip800 = ['E02','E03','A01','B01']
+skip800 = ['B01',]

 # some static revs/scws to be removed
 ignore_orbits=[352,834,912,1019,1021,1028,2275,2405,2493]
@@ -259,8 +258,6 @@ hdu.data=count_map
 hdu.writeto(proddir+fn.replace(".fits",".crab_count_map.fits"), overwrite=True)

 npoly=4
-if(enkey in ['E11','E12',]):
-    npoly=0
    
 z = np.polyfit(poly_x, poly_y, npoly)

@@ -282,17 +279,7 @@ indices = sorted(
 )

 coldefs = fits.ColDefs([
-    #fits.Column(name='OBSID', format='11A', array=[obs_id[index] for index in indices]), 
-    #fits.Column(name='RA', format='D', unit='deg',  array=[ra[index] for index in indices]),
-    #fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]),
-    #fits.Column(name='LON', format='D', unit='deg', array=[lon0[index] for index in indices]),
-    #fits.Column(name='LAT', format='D', unit='deg', array=[lat0[index] for index in indices]),
    fits.Column(name='REV', format='J', unit='', array=[rev0[index] for index in indices]),
-    #fits.Column(name='MJD', format='D', unit='', array=[mjd0[index] for index in indices]),
-    #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='RESID', format='D', unit='cts/s', array=[resid0[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='ERR', format='D', unit='', array=[err0[index] for index in indices]),
--- a/scripts/01_crabmodel.sh
+++ b/scripts/01_crabmodel.sh
@@ -1,49 +1,10 @@
-./01_crabmodel.py E01
-#./01_crabmodel.py E02
-#./01_crabmodel.py E03
-#./01_crabmodel.py E04
-#./01_crabmodel.py E05
-#./01_crabmodel.py E06
-#./01_crabmodel.py E07
-#./01_crabmodel.py E08
-#./01_crabmodel.py E09
-#./01_crabmodel.py E10
-#./01_crabmodel.py E11
-#./01_crabmodel.py E12
-./01_crabmodel.py E13
-./01_crabmodel.py E14
-#./01_crabmodel.py E15
+for band in A01 A02 A03
+do
+./01_crabmodel.py $band 
+done

-#./01_crabmodel.py A01
-#./01_crabmodel.py A02
-#./01_crabmodel.py A03
-#./01_crabmodel.py A04
-#./01_crabmodel.py A05
-#./01_crabmodel.py A06
-#./01_crabmodel.py A07
-#./01_crabmodel.py A08
-#./01_crabmodel.py A09
-#./01_crabmodel.py A10
-#./01_crabmodel.py A11
+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
+./01_crabmodel.py $band 
+done

-./01_crabmodel.py B01
-./01_crabmodel.py B02
-./01_crabmodel.py B03
-./01_crabmodel.py B04
-./01_crabmodel.py B05
-./01_crabmodel.py B06
-./01_crabmodel.py B07
-./01_crabmodel.py B08
-./01_crabmodel.py B09
-./01_crabmodel.py B10
-./01_crabmodel.py B11
-./01_crabmodel.py B12
-./01_crabmodel.py B13
-./01_crabmodel.py B14
-./01_crabmodel.py B15
-./01_crabmodel.py B16
-./01_crabmodel.py B17
-./01_crabmodel.py B18
-./01_crabmodel.py B19
-./01_crabmodel.py B20
-./01_crabmodel.py B21
--- a/scripts/02_grxe_resid.py
+++ b/scripts/02_grxe_resid.py
@@ -39,7 +39,7 @@ from ridge.config import *
 plotme=False
 enkey = sys.argv[1]
 outkey = sys.argv[2]
-#outkey = "BKG"
+

 fn="detcnts.{}.fits".format(enkey)

@@ -56,8 +56,6 @@ 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))
@@ -100,8 +98,6 @@ sx=int(hdulist[0].header['NAXIS1'])
 sy=int(hdulist[0].header['NAXIS2'])


-#d  = fits.getdata(datadir+fn)
-#df = pd.DataFrame(np.array(d).byteswap().newbyteorder())

 with fits.open(datadir+fn) as data:
    df = pd.DataFrame(data[1].data)
@@ -121,11 +117,7 @@ if(outkey=='ALL'):
    df = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & PHASE > {} & PHASE < {}'.format(revmin,revmax,phmin,phmax))
    print("N={}".format(df.shape[0]))

-"""
-GAL 70226 ScWs, 111.7 Ms
-BKG 61214 ScWs, 114.3 Ms
-Total 131440 ScWs, 225.9 Ms
-"""
+


 # fill AITOF map indexes
@@ -213,7 +205,7 @@ for i, row in df.iterrows():

 print("N={} ScWs, {:.1f} Ms".format(len(resid0),np.sum(texp0)/1e6))

-sigma=3
+
 rev_min=np.min(rev0)
 rev_max=np.max(rev0)
 orbits=np.array(rev0)
--- a/scripts/02_grxe_resid.sh
+++ b/scripts/02_grxe_resid.sh
@@ -1,79 +1,14 @@
-./02_grxe_resid.py E01 ALL
-./02_grxe_resid.py E01 BKG
-#./02_grxe_resid.py E02 ALL
-#./02_grxe_resid.py E03 ALL
-#./02_grxe_resid.py E04 ALL
-#./02_grxe_resid.py E05 ALL
-#./02_grxe_resid.py E06 ALL
-#./02_grxe_resid.py E07 ALL
-#./02_grxe_resid.py E08 ALL
-#./02_grxe_resid.py E09 ALL
-#./02_grxe_resid.py E10 ALL
-#./02_grxe_resid.py E11 ALL
-#./02_grxe_resid.py E12 ALL
+for band in A01 A02 A03
+do
+./02_grxe_resid.py $band ALL
+./02_grxe_resid.py $band BKG
+./02_grxe_resid.py $band GAL
+done

-./02_grxe_resid.py E13 ALL
-./02_grxe_resid.py E14 ALL
+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
+./02_grxe_resid.py $band ALL
+./02_grxe_resid.py $band BKG
+./02_grxe_resid.py $band GAL
+done

-./02_grxe_resid.py E13 BKG
-./02_grxe_resid.py E14 BKG
-
-#./02_grxe_resid.py E15 ALL
-
-#./02_grxe_resid.py A01 ALL
-#./02_grxe_resid.py A02 ALL
-#./02_grxe_resid.py A03 ALL
-#./02_grxe_resid.py A04 ALL
-#./02_grxe_resid.py A05 ALL
-#./02_grxe_resid.py A06 ALL
-#./02_grxe_resid.py A07 ALL
-#./02_grxe_resid.py A08 ALL
-#./02_grxe_resid.py A09 ALL
-#./02_grxe_resid.py A10 ALL
-#./02_grxe_resid.py A11 ALL
-
-#./02_grxe_resid.py B01 ALL
-#./02_grxe_resid.py B02 ALL
-#./02_grxe_resid.py B03 ALL
-#./02_grxe_resid.py B04 ALL
-#./02_grxe_resid.py B05 ALL
-#./02_grxe_resid.py B06 ALL
-#./02_grxe_resid.py B07 ALL
-#./02_grxe_resid.py B08 ALL
-#./02_grxe_resid.py B09 ALL
-#./02_grxe_resid.py B10 ALL
-#./02_grxe_resid.py B11 ALL
-#./02_grxe_resid.py B12 ALL
-#./02_grxe_resid.py B13 ALL
-#./02_grxe_resid.py B14 ALL
-#./02_grxe_resid.py B15 ALL
-#./02_grxe_resid.py B16 ALL
-#./02_grxe_resid.py B17 ALL
-#./02_grxe_resid.py B18 ALL
-#./02_grxe_resid.py B19 ALL
-#./02_grxe_resid.py B20 ALL
-#./02_grxe_resid.py B21 ALL
-
-./02_grxe_resid.py B01 ALL
-./02_grxe_resid.py B02 ALL
-./02_grxe_resid.py B03 ALL
-./02_grxe_resid.py B04 ALL
-./02_grxe_resid.py B05 ALL
-./02_grxe_resid.py B06 ALL
-./02_grxe_resid.py B07 ALL
-./02_grxe_resid.py B08 ALL
-./02_grxe_resid.py B09 ALL
-./02_grxe_resid.py B10 ALL
-./02_grxe_resid.py B11 ALL
-./02_grxe_resid.py B12 ALL
-./02_grxe_resid.py B13 ALL
-./02_grxe_resid.py B14 ALL
-./02_grxe_resid.py B15 ALL
-./02_grxe_resid.py B16 ALL
-./02_grxe_resid.py B17 ALL
-./02_grxe_resid.py B18 ALL
-./02_grxe_resid.py B19 ALL
-./02_grxe_resid.py B20 ALL
-./02_grxe_resid.py B21 ALL
-
-./02_grxe_resid.py B21 BKG
--- a/scripts/README.md
+++ b/scripts/README.md
@@ -1,32 +1,34 @@
-Подготовка рабочего окружения:
-```
-source <MY PATH>/venv/bin/activate.csh
-```

-Номер скрипта означает последовательность выполнения. Результаты работы скриптов помещаются в ```ridge/products/```.
+
+Don't forget to activate environment: `source <MY PATH>/venv/bin/activate`
+
+The script number indicates the sequence of execution. The results are placed in `ridge/products/`. The corresponding shell scripts run Python code with different parameters.

 ### 01_bgdmodel.py / 01_bgdmodel.sh

-Создает модель фона.
+Calibrates IBIS/ISGRI background model

 ### 01_crabmodel.py / 01_crabmodel.sh

-Создает модель скорости счета Краба.
+Calibrates Crab detector count rate model.

 ### 02_grxe_resid.py / 01_grxe_resid.sh

-Вычисляет остатки после вычитания модели из данных в единицах мКраб.
+Calculates difference between detector count rate and that predicted by background model in mCrab units for `GAL` (Galaxy), `BKG` (Extragalactic), or `ALL` (all-sky) regions.

 ### 02_grxe_map.py / 01_grxe_map.sh

-Создает карту остатков в единицах мКраб.
+Makes the map of the residuals in mCrab units (not covered in the paper).

 ### 03_grxe_spec.py

-Создает спектр выбранного региона неба, например:
+Extract spectrum of selected sky region, e.g.:

-```./03_grxe_spec.py GC``` или  ```./03_grxe_spec.py M81```
+```./03_grxe_spec.py GC``` or  ```./03_grxe_spec.py M81```

-Регионы неба задаются в ```config.py```. Спектры помещаются в директорию ```./products/spec```. Затем надо перейти в эту директорию и выполнить, например: ```perl ../../scripts/do_pha_igr_v3_mCrab.pl M81.spec```.
+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```.

-Если не задать регион неба, скрипт делает спектры для всех регионов.