diff --git a/.gitignore b/.gitignore index 608b3f9..e46bb43 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,7 @@ products*/ data/archive data/processed data/processed-oot +data/chandra wheels/ share/python-wheels/ *.egg-info/ diff --git a/arches/arches/config.py b/arches/arches/config.py index b62f448..0257f6d 100644 --- a/arches/arches/config.py +++ b/arches/arches/config.py @@ -1,6 +1,26 @@ from pathlib import Path -good=['0862470801','0862470701','0862470601'] +# Arches 2020, Maica Clavel: +# good=['0862470801','0862470701','0862470601'] + +# after 2020: MOS1 is OFF, PN at border +good=[ + #'0862470801', # 2020 Maica Clavel + #'0862470701', # 2020 Maica Clavel + #'0862470601', # 2020 Maica Clavel + '0860620401', # 2021-03-30 M1 off, M2 OK + '0860620201', # 2021-03-28 M1 off, M2 OK, PN border + '0860620301', # 2021-03-29 M1 off, M2 good, PN border + '0893811101', # 2022-03-18 M1 off, M2 OK, PN border Transient nearby + '0893811301', # 2022-03-26 M1 off, M2 OK, PN border, transient nearby + '0944580501', # 2024-03-29 M1 off, M2 good, PN border, note: "U" exposure + '0944580801', # 2024-03-28 M1 off, M2 OK + ##'0944580701', low exposure + ##'0860620501', # 2021-03-31 M1 off, M2 OK + ##'0932392301', # 2024-04-03 M1 off, M2 Bad, no exposure + ##'0944580601', # 2024-03-26 M1 off, M2 OK, but low exposure, PN border +] + """ The angular size of an XMM DETXY physical coordinate system is 0.05 arcseconds """ det_pix_as=0.05 @@ -26,11 +46,12 @@ skip=['0862471401','0862471501','0862470501'] # energy bands taken from Tatischeff 2012 emin=[ - '2000', # continuum, 2-4 keV - '4170', # continuum, 4.17-5.86 keV (T12) - '6300', # 6.4 keV line, 6.3-6.48 keV (T12) - '6564', # 6.7 keV line, 6.564-6.753 keV (T12) - '8000', # high-energy continuum, 8-12 keV + '2000', # [0] continuum, 2-4 keV + '4170', # [1] continuum, 4.17-5.86 keV (T12) + '6300', # [2] 6.4 keV line, 6.3-6.48 keV (T12) + '6564', # [3] 6.7 keV line, 6.564-6.753 keV (T12) + '8000', # [4] high-energy continuum, 8-12 keV + '8000', # [5] high-energy continuum, 8-10 keV ] emax=[ '4000', @@ -38,4 +59,17 @@ emax=[ '6480', '6753', '12000', + '10000', ] + +""" +nh=9.5, Gamma=3 + Model Flux 1.5475e-05 photons (1.2204e-13 ergs/cm^2/s) range (4.1700 - 5.8600 keV) + Model Flux 9.3046e-07 photons (9.5245e-15 ergs/cm^2/s) range (6.3000 - 6.4800 keV) +0.0601267 + Model Flux 8.7952e-07 photons (9.3813e-15 ergs/cm^2/s) range (6.5640 - 6.7530 keV) +0.0568349 + Model Flux 6.195e-06 photons (9.5642e-14 ergs/cm^2/s) range (8.0000 - 12.000 keV) +0.400323 +""" + diff --git a/arches/arches/utils.py b/arches/arches/utils.py index c85fd86..d0ec503 100644 --- a/arches/arches/utils.py +++ b/arches/arches/utils.py @@ -65,20 +65,30 @@ def remove_file(filename): def move_files(pattern, destination): for file in glob.glob(pattern): shutil.move(file, os.path.join(destination,file)) -def test_file(fn): + +def test_file(fn, alive=False): if not (os.path.isfile(fn)==True): print(f'requested filename {fn} is not found') cwd = os.getcwd() print(f"Current directory: {cwd}") - sys.exit() + if(alive): + return False + else: + print("### Abort ###") + sys.exit() + return True -def get_first_file(pattern): +def get_first_file(pattern, alive=False): files = glob.glob(pattern) if not (files): print(f"Files not found from pattern {pattern}") - sys.exit() - filename = files[0] - return filename + if(alive): + return None + else: + print("### Abort ###") + sys.exit() + return files[0] + # Function to plot Lightcurve @@ -628,12 +638,12 @@ def run_evqpb(key, work_dir=None, evtfile=None, attfile=None, pimin=[200,], pima result = subprocess.run(cmd, capture_output=True, text=True, check=True) """ -def run_mosaic(outfile_cts=None,outfile_qpb=None,outfile_exp=None, - outfile_sub=None,outfile_flx=None,outfile_pix=None, - cts=None,qpb=None,exp=None,nn=3,devmax=5.0,cutbox=100): +def run_mosaic(out_cts=None,out_qpb=None,out_exp=None, + out_sub=None,out_flx=None,out_pix=None, + cts=None,qpb=None,exp=None,nn=2,devmax=5.0,cutbox=100,fracexp=0.25): ref_crd = SkyCoord(arches_ra, arches_dec, frame=FK5(), unit="deg") nn2=nn*nn - print(outfile_cts) + print(out_cts) print(cts[0]) print(exp[0]) # take first image as reference @@ -684,19 +694,19 @@ def run_mosaic(outfile_cts=None,outfile_qpb=None,outfile_exp=None, map_cts[col_index][row_index]=2.0 ref_hdul[0].data=map_cts - ref_hdul.writeto(outfile_cts,overwrite=True) + ref_hdul.writeto(out_cts,overwrite=True) sys.exit() """ for idx,in_cts in enumerate(cts): - #if not (idx==2): - # continue + hdul0 = fits.open(in_cts) c_data = hdul0[0].data c_header = hdul0[0].header c_wcs = WCS(c_header) hdul0.close() + print(f"{in_cts} filter {c_header['FILTER']}") hdul0 = fits.open(qpb[idx]) @@ -710,7 +720,9 @@ def run_mosaic(outfile_cts=None,outfile_qpb=None,outfile_exp=None, e_header = hdul0[0].header e_wcs = WCS(e_header) hdul0.close() - print(exp[idx]) + cutexp=np.max(e_data)*fracexp + print(f'{exp[idx]}, fracexp: {fracexp} cutexp: {cutexp:.1f} sec') + xx, yy = c_wcs.world_to_pixel(ref_crd) @@ -740,7 +752,7 @@ def run_mosaic(outfile_cts=None,outfile_qpb=None,outfile_exp=None, exp_xx, exp_yy = e_wcs.world_to_pixel(c_sky) exp_x = int(np.rint(np.float64(exp_xx))) exp_y = int(np.rint(np.float64(exp_yy))) - if not (e_data[exp_y, exp_x]>0): + if not (e_data[exp_y, exp_x]>cutexp): continue sep_arcmin = c_sky.separation(ref_crd).arcmin @@ -775,39 +787,37 @@ def run_mosaic(outfile_cts=None,outfile_qpb=None,outfile_exp=None, # Convert summed sky back to normal units - for x in range(nx): - for y in range(ny): - if(map_exp[y][x]>0): - map_flx[y][x]=(map_cts[y][x]-map_qpb[y][x])/map_exp[y][x] - map_cts[y][x]=map_cts[world_to_pixely][x]/map_pix[y][x] - map_qpb[y][x]=map_qpb[y][x]/map_pix[y][x] - map_exp[y][x]=map_exp[y][x]/map_pix[y][x] - else: - map_flx[y][x]=0.0 - - if(outfile_pix): + + index = map_exp > 0 + map_flx[index]=(map_cts[index]-map_qpb[index])/map_exp[index] + #map_cts[index]=map_cts[index]/map_pix[index] + #map_qpb[index]=map_qpb[index]/map_pix[index] + #map_exp[index]=map_exp[index]/map_pix[index] + + + if(out_pix): ref_hdul[0].data=map_pix - ref_hdul.writeto(outfile_pix,overwrite=True) + ref_hdul.writeto(out_pix,overwrite=True) ref_hdul[0].data=map_cts - ref_hdul.writeto(outfile_cts,overwrite=True) + ref_hdul.writeto(out_cts,overwrite=True) ref_hdul[0].data=map_exp - ref_hdul.writeto(outfile_exp,overwrite=True) + ref_hdul.writeto(out_exp,overwrite=True) ref_hdul[0].data=map_qpb - ref_hdul.writeto(outfile_qpb,overwrite=True) + ref_hdul.writeto(out_qpb,overwrite=True) ref_hdul[0].data=map_flx - ref_hdul.writeto(outfile_flx,overwrite=True) + ref_hdul.writeto(out_flx,overwrite=True) #sys.exit() """ - cmd = ['farith', f'{outfile_cts}', f'{outfile_qpb}', f'{outfile_sub}', 'SUB', + cmd = ['farith', f'{out_cts}', f'{out_qpb}', f'{out_sub}', 'SUB', 'copyprime=yes', 'clobber=yes'] result = subprocess.run(cmd, capture_output=True, text=True, check=True) - cmd = ['farith', f'{outfile_sub}', f'{outfile_exp}', f'{outfile_rat}', 'DIV', + cmd = ['farith', f'{out_sub}', f'{out_exp}', f'{out_rat}', 'DIV', 'copyprime=yes', 'clobber=yes'] result = subprocess.run(cmd, capture_output=True, text=True, check=True) """ diff --git a/data/ds9reg/arches.grd b/data/ds9reg/arches.grd new file mode 100644 index 0000000..c138ba6 --- /dev/null +++ b/data/ds9reg/arches.grd @@ -0,0 +1,2 @@ +global grid +array set grid { grid,color white grid,style 1 view 1 axes,origin lll numlab,type interior border 1 tick,width 1 numlab,color white numlab,weight normal title,text {} axes,color white numlab,slant roman axes,style 0 textlab,color black numlab 1 skyformat degrees textlab,gap1 {} border,color cyan textlab,slant roman textlab,gap2 {} textlab,size 10 grid,gapunit1 degrees border,style 0 grid,gapunit2 degrees grid,gapunit3 pixels title,gap {} textlab,font helvetica grid,width 1 format1 .2f title 1 format2 .2f sky galactic textlab 1 title,color black axes,width 1 title,slant roman border,width 1 system wcs numlab,vertical 0 tick,color cyan textlab,def1 1 tick,style 0 textlab,def2 1 axes 0 type analysis frame Frame4 grid,gap1 0.05 tick 1 grid,gap2 0.05 numlab,gap1 {} grid,gap3 {} numlab,gap2 {} grid 1 numlab,size 12 numlab,gap3 {} axes,type interior textlab,weight normal title,size 12 numlab,font helvetica title,def 1 title,font helvetica textlab,text1 {} title,weight normal textlab,text2 {} } diff --git a/data/xspec/load_merged.xcm b/data/xspec/load_merged.xcm new file mode 100644 index 0000000..7eea333 --- /dev/null +++ b/data/xspec/load_merged.xcm @@ -0,0 +1,7 @@ +data EPIC_merged_spectrum_grp.fits +ign **-4. +ign 10.-** +setpl en +cpd /xs +@../../../data/xspec/pow.xcm +pl lda diff --git a/data/xspec/pow.xcm b/data/xspec/pow.xcm new file mode 100644 index 0000000..5eea46c --- /dev/null +++ b/data/xspec/pow.xcm @@ -0,0 +1,19 @@ +method leven 10 0.01 +abund wilm +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model wabs*pegpwrlw + gaussian + gaussian + 7 -1 0 0 100000 1e+06 + 2 -1 -3 -2 9 10 + 4.17 -0.01 -100 -100 1e+10 1e+10 + 5.86 -0.01 -100 -100 1e+10 1e+10 + 0.155061 0.01 0 0 1e+20 1e+24 + 6.4 -1 0 0 1e+06 1e+06 + 0.01 -1 0 0 10 20 + 4.20779e-07 0.01 0 0 1e+20 1e+24 + 6.68402 0.05 0 0 1e+06 1e+06 + 0.01 -1 0 0 10 20 + 6.01672e-06 0.01 0 0 1e+20 1e+24 +bayes off diff --git a/scripts/01_init_events-oot.py b/scripts/01_init_events-oot.py index af956e9..6eeed8d 100755 --- a/scripts/01_init_events-oot.py +++ b/scripts/01_init_events-oot.py @@ -23,11 +23,14 @@ inargs = ['--version'] t = w('sasver', inargs) t.run() -files = glob.glob(archive_dir+'/0862*') +files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - + + if not (obsid in after2020): + continue + local_ccf=f'{events_dir}/{obsid}/ccf.cif' work_dir=f'{events_dir}/{obsid}' diff --git a/scripts/01_init_events.py b/scripts/01_init_events.py index ac615d5..b4366e4 100755 --- a/scripts/01_init_events.py +++ b/scripts/01_init_events.py @@ -23,18 +23,26 @@ inargs = ['--version'] t = w('sasver', inargs) t.run() -files = glob.glob(archive_dir+'/0862*') +# Какие наклоны разрешены для какой температуры +# W1 жесткость от центра? Муно и др. +# есть ли двойные системы? Кларк и др. + +files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) + + if not (obsid in good): + continue local_ccf=f'{events_dir}/{obsid}/ccf.cif' work_dir=f'{events_dir}/{obsid}' + # downloads data #inargs = [f'odfid={obsid}',f'workdir={work_dir}'] #w('startsas', inargs).run() - os.environ['SAS_ODF'] = f'{archive_dir}/{obsid}' + os.environ['SAS_ODF'] = f'{archive_dir}/{obsid}/odf' # ??? create_folder(work_dir) diff --git a/scripts/02_sas_flares_mos.py b/scripts/02_sas_flares_mos.py index 7d10139..3d26e30 100755 --- a/scripts/02_sas_flares_mos.py +++ b/scripts/02_sas_flares_mos.py @@ -29,7 +29,7 @@ events_dir=root_path+'/data/processed' products_dir=root_path+'/products/sas' ds9reg_dir=root_path+'/data/ds9reg' -imos=1 +imos=2 create_folder(products_dir) @@ -42,14 +42,17 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + + if not (obsid in good): continue + + work_dir = init_work_dir(obsid, products_dir=products_dir) os.chdir(work_dir) - search_str = f'{events_dir}/{obsid}/????_{obsid}_EMOS{imos}_S???_ImagingEvts.ds' + search_str = f'{events_dir}/{obsid}/????_{obsid}_EMOS{imos}_????_ImagingEvts.ds' print(search_str) emfiles = glob.glob(search_str) if not (emfiles): @@ -312,33 +315,34 @@ for obsid in files: # Create Clean Events image - xmax=np.amax(gti_evt_data['X']) - xmin=np.amin(gti_evt_data['X']) - xmid=(xmax-xmin)/2.+xmin - ymax=np.amax(gti_evt_data['Y']) - ymin=np.amin(gti_evt_data['Y']) - xbin_size=80 - ybin_size=80 - NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size)) + if(len(gti_evt_data['X'])>10): + xmax=np.amax(gti_evt_data['X']) + xmin=np.amin(gti_evt_data['X']) + xmid=(xmax-xmin)/2.+xmin + ymax=np.amax(gti_evt_data['Y']) + ymin=np.amin(gti_evt_data['Y']) + xbin_size=80 + ybin_size=80 + NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size)) - plt.subplot(1, 2, pl) + plt.subplot(1, 2, pl) - img_zero_mpl = plt.hist2d(gti_evt_data['X'], gti_evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm()) + img_zero_mpl = plt.hist2d(gti_evt_data['X'], gti_evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm()) - cbar = plt.colorbar(ticks=[10.,100.,1000.]) - cbar.ax.set_yticklabels(['10','100','1000']) + cbar = plt.colorbar(ticks=[10.,100.,1000.]) + cbar.ax.set_yticklabels(['10','100','1000']) - plt.title(out_clean_evtFile) - plt.xlabel('x') - plt.ylabel('y') + plt.title(out_clean_evtFile) + plt.xlabel('x') + plt.ylabel('y') - plt.text(xmid, ymin-0.1*(ymax-ymin), expression, ha='center') + plt.text(xmid, ymin-0.1*(ymax-ymin), expression, ha='center') - pl=pl+1 + pl=pl+1 - gti_hdu_list.close() - hdu_list.close() - plt.show() + gti_hdu_list.close() + hdu_list.close() + plt.show() ############################################################################### # Define some parameters to produce the image and the name of the output file # diff --git a/scripts/02_sas_flares_pn.py b/scripts/02_sas_flares_pn.py index 071f9c6..f248607 100755 --- a/scripts/02_sas_flares_pn.py +++ b/scripts/02_sas_flares_pn.py @@ -35,23 +35,25 @@ inargs = ['--version'] t = w('sasver', inargs) t.run() -files = glob.glob(archive_dir+'/0862*') +files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) os.chdir(work_dir) - search_str = f'{events_dir}/{obsid}/????_{obsid}_EPN_S???_ImagingEvts.ds' + search_str = f'{events_dir}/{obsid}/????_{obsid}_EPN_????_ImagingEvts.ds' print(search_str) epfiles = glob.glob(search_str) if not (epfiles): print("*** run 01_init_events.py ***") - #w('epproc', []).run() + #w('epproc', []).run() + continue + eventfile = epfiles[0] print("Checking for EPIC-pn Event Files ..... \n") diff --git a/scripts/03_sas_image_mos.py b/scripts/03_sas_image_mos.py index 1cc960a..80aa41a 100755 --- a/scripts/03_sas_image_mos.py +++ b/scripts/03_sas_image_mos.py @@ -41,7 +41,7 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) diff --git a/scripts/03_sas_image_pn.py b/scripts/03_sas_image_pn.py index 3b627e1..96ae8b5 100755 --- a/scripts/03_sas_image_pn.py +++ b/scripts/03_sas_image_pn.py @@ -39,7 +39,7 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) diff --git a/scripts/04_sas_spectrum_mos.py b/scripts/04_sas_spectrum_mos.py index 0d656da..878b052 100755 --- a/scripts/04_sas_spectrum_mos.py +++ b/scripts/04_sas_spectrum_mos.py @@ -47,7 +47,7 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) diff --git a/scripts/04_sas_spectrum_pn.py b/scripts/04_sas_spectrum_pn.py index 3b3e934..98636c4 100755 --- a/scripts/04_sas_spectrum_pn.py +++ b/scripts/04_sas_spectrum_pn.py @@ -40,16 +40,16 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) os.chdir(work_dir) - eventfile=f'{products_dir}/{obsid}/EPIC_pn_gtiFilteredEvts.ds' + eventfile=f'{products_dir}/{obsid}/EPIC_PN_gtiFilteredEvts.ds' if(os.path.isfile(eventfile)==False): - print("*** run 02_filter_flares_pn.py ***") + print("*** run 02_sas_flares_pn.py ***") sys.exit() # EPIC_pn_gtiFilteredEvts.ds @@ -64,7 +64,7 @@ for obsid in files: # Visualize the image with ds9 src_fn=f'{ds9reg_dir}/{obsid}/src.reg' - bkg_fn=f'{ds9reg_dir}/{obsid}/bkg.reg' + bkg_fn=f'{ds9reg_dir}/{obsid}/ann.reg' x_source,y_source,r_source,x_bkg,y_bkg,r_bkg,r2_bkg = get_ds9_regions(out_IMFile, src_fn, bkg_fn) ######################################### diff --git a/scripts/05_sas_spectrum_merge.py b/scripts/05_sas_spectrum_merge.py index 565e92f..05e6bc5 100755 --- a/scripts/05_sas_spectrum_merge.py +++ b/scripts/05_sas_spectrum_merge.py @@ -41,7 +41,7 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) @@ -55,7 +55,7 @@ for obsid in files: bkg=[] resp=[] arf=[] - for imos in [1,2]: + for imos in [2]: in_SPSRCFile = work_dir+f'/EPIC_MOS{imos}_source_spectrum.fits' # Name of the output source spectrum in_SPBKGFile = work_dir+f'/EPIC_MOS{imos}_background_spectrum.fits' # Name of the output background spectrum in_RESPFile = work_dir+f'/EPIC_MOS{imos}.rmf' # Name of the output redistribution diff --git a/scripts/06_sas_QPB.py b/scripts/06_sas_QPB.py index 0db8f6e..7c1b751 100755 --- a/scripts/06_sas_QPB.py +++ b/scripts/06_sas_QPB.py @@ -41,7 +41,7 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) @@ -55,7 +55,8 @@ for obsid in files: attfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_AttHk.ds') #for imos in [1,2]: - for key in ['MOS1','MOS2','PN']: + #for key in ['MOS1','MOS2','PN']: + for key in ['MOS2']: evtfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_E{key}*_ImagingEvts.ds') # first element contains all energies diff --git a/scripts/07_sas_eimageget.py b/scripts/07_sas_eimageget.py index 3b73579..2142393 100755 --- a/scripts/07_sas_eimageget.py +++ b/scripts/07_sas_eimageget.py @@ -43,7 +43,7 @@ files = glob.glob(archive_dir+'/*') for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = init_work_dir(obsid, products_dir=products_dir) @@ -55,8 +55,11 @@ for obsid in files: attfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_AttHk.ds') - for key in ['MOS1','MOS2','PN']: + #for key in ['MOS1','MOS2','PN']: + for key in ['MOS2']: evtfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_E{key}*_ImagingEvts.ds') + if not (os.path.isfile(evtfile)==True): + continue fwcfile = get_first_file(fwc_dir+f'/{key}_closed_FF_2025_v1.fits') pimin=" ".join(emin) pimax=" ".join(emax) diff --git a/scripts/08_sas_mosaic.py b/scripts/08_sas_mosaic.py index 6f0368f..563fb8b 100755 --- a/scripts/08_sas_mosaic.py +++ b/scripts/08_sas_mosaic.py @@ -29,8 +29,9 @@ archive_dir=root_path+'/data/archive' events_dir=root_path+'/data/processed' products_dir=root_path+'/products/sas' create_folder(products_dir) -mosaic_dir=root_path+'/products/mosaic' -create_folder(mosaic_dir) + +mdir=root_path+'/products/mosaic_2021-2024' +create_folder(mdir) ds9reg_dir=root_path+'/data/ds9reg' @@ -44,6 +45,9 @@ files = glob.glob(archive_dir+'/*') # run over energies for idx,e in enumerate(emin): + #if not (idx==1): + # continue + label=f"_{idx}" e1=emin[idx] e2=emax[idx] @@ -75,24 +79,28 @@ for idx,e in enumerate(emin): for obsid in files: obsid = os.path.basename(obsid) - if(obsid in skip): + if not (obsid in good): continue work_dir = f'{products_dir}/{obsid}' for ikey,key in enumerate(['MOS1','MOS2','PN']): shortkey = key.replace("OS","") - expkey=f'S00{ikey+1}' + #expkey=f'00{ikey+1}' + #expkey=f'00{ikey+1}' # sci and qpb image are produced by 06_sas_QPB.py sci_image=work_dir+f'/EPIC_{key}_sci_image{label}.fits' - test_file(sci_image) + if not (test_file(sci_image, alive=True)): + continue qpb_image=work_dir+f'/EPIC_{key}_qpb_image{label}.fits' - test_file(qpb_image) + if not (test_file(qpb_image, alive=True)): + continue # exposure image is produced by eimageget (07_sas_eimageget.py) # example: P0862470801PNS003_ima_4exp.fits - exp_image=work_dir+f'/P{obsid}{shortkey}{expkey}_ima{label}exp.fits' - test_file(exp_image) - + exp_image = get_first_file(work_dir+f'/P{obsid}{shortkey}?00?_ima{label}exp.fits', alive=True) + if not (exp_image): + continue + ima.append(sci_image) qpb.append(qpb_image) exp.append(exp_image) @@ -114,36 +122,26 @@ for idx,e in enumerate(emin): mos_qpb.append(qpb_image) mos_exp.append(exp_image) + print("MOS1 ",len(m1_ima)) + print(m1_ima) + print("MOS2 ",len(m2_ima)) + print(m2_ima) + print("MOS1+2 ",len(mos_ima)) + print(mos_ima) + print("PN ",len(pn_ima)) + print(pn_ima) + print("All ",len(ima)) + print(ima) + + run_mosaic(out_cts=mdir+f'/m1_cts{label}.fits', out_exp=mdir+f'/m1_exp{label}.fits', out_qpb=mdir+f'/m1_qpb{label}.fits', out_sub=mdir+f'/m1_sub{label}.fits', out_flx=mdir+f'/m1_flx{label}.fits', out_pix=mdir+f'/m1_pix{label}.fits', cts=m1_ima,qpb=m1_qpb,exp=m1_exp) + + run_mosaic(out_cts=mdir+f'/m2_cts{label}.fits', out_exp=mdir+f'/m2_exp{label}.fits', out_qpb=mdir+f'/m2_qpb{label}.fits', out_sub=mdir+f'/m2_sub{label}.fits', out_flx=mdir+f'/m2_flx{label}.fits', cts=m2_ima,qpb=m2_qpb,exp=m2_exp) + + run_mosaic(out_cts=mdir+f'/mos_cts{label}.fits', out_exp=mdir+f'/mos_exp{label}.fits', out_qpb=mdir+f'/mos_qpb{label}.fits', out_sub=mdir+f'/mos_sub{label}.fits', out_flx=mdir+f'/mos_flx{label}.fits', cts=mos_ima,qpb=mos_qpb,exp=mos_exp) + + run_mosaic(out_cts=mdir+f'/pn_cts{label}.fits', out_exp=mdir+f'/pn_exp{label}.fits', out_qpb=mdir+f'/pn_qpb{label}.fits', out_sub=mdir+f'/pn_sub{label}.fits', out_flx=mdir+f'/pn_flx{label}.fits', cts=pn_ima,qpb=pn_qpb,exp=pn_exp) + + + run_mosaic(out_cts=mdir+f'/cts{label}.fits', out_exp=mdir+f'/exp{label}.fits', out_qpb=mdir+f'/qpb{label}.fits', out_sub=mdir+f'/sub{label}.fits', out_flx=mdir+f'/flx{label}.fits', cts=ima,qpb=qpb,exp=exp) - run_mosaic(outfile_cts=mosaic_dir+f'/m1_cts{label}.fits', - outfile_exp=mosaic_dir+f'/m1_exp{label}.fits', - outfile_qpb=mosaic_dir+f'/m1_qpb{label}.fits', - outfile_sub=mosaic_dir+f'/m1_sub{label}.fits', - outfile_flx=mosaic_dir+f'/m1_flx{label}.fits', - outfile_pix=mosaic_dir+f'/m1_pix{label}.fits', - cts=m1_ima,qpb=m1_qpb,exp=m1_exp) - run_mosaic(outfile_cts=mosaic_dir+f'/m2_cts{label}.fits', - outfile_exp=mosaic_dir+f'/m2_exp{label}.fits', - outfile_qpb=mosaic_dir+f'/m2_qpb{label}.fits', - outfile_sub=mosaic_dir+f'/m2_sub{label}.fits', - outfile_flx=mosaic_dir+f'/m2_flx{label}.fits', - cts=m2_ima,qpb=m2_qpb,exp=m2_exp) - run_mosaic(outfile_cts=mosaic_dir+f'/mos_cts{label}.fits', - outfile_exp=mosaic_dir+f'/mos_exp{label}.fits', - outfile_qpb=mosaic_dir+f'/mos_qpb{label}.fits', - outfile_sub=mosaic_dir+f'/mos_sub{label}.fits', - outfile_flx=mosaic_dir+f'/mos_flx{label}.fits', - cts=mos_ima,qpb=mos_qpb,exp=mos_exp) - run_mosaic(outfile_cts=mosaic_dir+f'/pn_cts{label}.fits', - outfile_exp=mosaic_dir+f'/pn_exp{label}.fits', - outfile_qpb=mosaic_dir+f'/pn_qpb{label}.fits', - outfile_sub=mosaic_dir+f'/pn_sub{label}.fits', - outfile_flx=mosaic_dir+f'/pn_flx{label}.fits', - cts=pn_ima,qpb=pn_qpb,exp=pn_exp) - run_mosaic(outfile_cts=mosaic_dir+f'/cts{label}.fits', - outfile_exp=mosaic_dir+f'/exp{label}.fits', - outfile_qpb=mosaic_dir+f'/qpb{label}.fits', - outfile_sub=mosaic_dir+f'/sub{label}.fits', - outfile_flx=mosaic_dir+f'/flx{label}.fits', - cts=ima,qpb=qpb,exp=exp) diff --git a/scripts/09_sas_convol.py b/scripts/09_sas_convol.py new file mode 100755 index 0000000..6527494 --- /dev/null +++ b/scripts/09_sas_convol.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python + +from pysas.wrapper import Wrapper as w +import os, sys +from os.path import dirname +import inspect +import glob + +import os.path +from os import path +import subprocess +import numpy as np +import matplotlib.pyplot as plt +from astropy.io import fits +from astropy.table import Table +from matplotlib.colors import LogNorm + +from astropy.convolution import interpolate_replace_nans, convolve, convolve_fft +from astropy.convolution import Gaussian2DKernel + +import re +import pyds9 + +import arches +from arches.utils import * +from arches.config import * + +root_path=dirname(dirname(dirname(inspect.getfile(arches)))) +print("Arches root path: {}".format(root_path)) + +archive_dir=root_path+'/data/archive' +events_dir=root_path+'/data/processed' +products_dir=root_path+'/products/sas' +create_folder(products_dir) +mdir=root_path+'/products/mosaic' +create_folder(mdir) + +ds9reg_dir=root_path+'/data/ds9reg' + + +inargs = ['--version'] +t = w('sasver', inargs) +t.run() + +# run over energies +for idx,e in enumerate(emin): + label=f"_{idx}" + e1=emin[idx] + e2=emax[idx] + + for ikey,key in enumerate(['','m1_','m2_','mos_','pn_']): + in_cts=f'{mdir}/{key}cts{label}.fits' + if not (test_file(in_cts,alive=True)): + continue + + hdul0 = fits.open(in_cts) + c_data = hdul0[0].data + c_header = hdul0[0].header + c_wcs = WCS(c_header) + hdul0.close() + + in_exp=f'{mdir}/{key}exp{label}.fits' + if not (test_file(in_exp,alive=True)): + continue + + hdul0 = fits.open(in_exp) + e_data = hdul0[0].data + e_header = hdul0[0].header + e_wcs = WCS(e_header) + hdul0.close() + + in_qpb=f'{mdir}/{key}qpb{label}.fits' + if not (test_file(in_qpb,alive=True)): + continue + + hdul0 = fits.open(in_qpb) + q_data = hdul0[0].data + q_header = hdul0[0].header + q_wcs = WCS(q_header) + hdul0.close() + + gauss2d = Gaussian2DKernel(x_stddev=0.530785) + + convol_cts = convolve(c_data, gauss2d) + out_cts = f'{mdir}/{key}cts{label}_gauss.fits' + hdul0[0].data=convol_cts + hdul0.writeto(out_cts,overwrite=True) + + convol_qpb = convolve(q_data, gauss2d) + out_qpb = f'{mdir}/{key}qpb{label}_gauss.fits' + hdul0[0].data=convol_qpb + hdul0.writeto(out_qpb,overwrite=True) + + (nx,ny) = c_data.shape + map_flx = np.zeros(c_data.shape,dtype=np.float64) + index = e_data > 0 + map_flx[index]=(convol_cts[index]-convol_qpb[index])/e_data[index] + out_flx = f'{mdir}/{key}flx{label}_gauss.fits' + hdul0[0].data=map_flx + hdul0.writeto(out_flx,overwrite=True) +