From 777c92ef2e7dc9c4565ca1f44b9f360952537c31 Mon Sep 17 00:00:00 2001 From: Roman Krivonos Date: Tue, 28 Oct 2025 00:11:18 +0300 Subject: [PATCH] good night --- arches/arches/config.py | 16 +++ arches/arches/utils.py | 270 ++++++++++++++++++++++++++++++++++++ data/FWC/README.md | 3 + scripts/06_sas_QPB.py | 81 ++--------- scripts/07_sas_eimageget.py | 77 ++++++++++ scripts/08_sas_mosaic.py | 149 ++++++++++++++++++++ 6 files changed, 526 insertions(+), 70 deletions(-) create mode 100644 data/FWC/README.md create mode 100755 scripts/07_sas_eimageget.py create mode 100755 scripts/08_sas_mosaic.py diff --git a/arches/arches/config.py b/arches/arches/config.py index 762a4c3..b62f448 100644 --- a/arches/arches/config.py +++ b/arches/arches/config.py @@ -23,3 +23,19 @@ Radius=15" """ 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 +] +emax=[ + '4000', + '5860', + '6480', + '6753', + '12000', +] diff --git a/arches/arches/utils.py b/arches/arches/utils.py index 7d909af..c85fd86 100644 --- a/arches/arches/utils.py +++ b/arches/arches/utils.py @@ -15,6 +15,7 @@ from pathlib import Path import pandas import pickle import pyds9 +import subprocess from os.path import dirname import inspect @@ -27,6 +28,7 @@ from astropy.coordinates import SkyCoord # High-level coordinates from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames from astropy.coordinates import Angle, Latitude, Longitude # Angles from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch +from astropy.convolution import interpolate_replace_nans import statistics import shutil @@ -63,6 +65,12 @@ 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): + 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() def get_first_file(pattern): files = glob.glob(pattern) @@ -541,3 +549,265 @@ def get_ds9_regions(image, src_fn, bkg_fn): r2_bkg = regionbkg.split(",")[3].replace('\n','') print(" r2_bkg = ", r2_bkg, "(physical, annulus)") return x_source,y_source,r_source,x_bkg,y_bkg,r_bkg,r2_bkg + + +def run_evqpb(key, work_dir=None, evtfile=None, attfile=None, pimin=[200,], pimax=[12000,], label=['',], exposurefactor=10.0): + + test_file(evtfile) + test_file(attfile) + + + # Run the SAS task evqpb over each one of _original_ event files. + evqpb_outset=work_dir+f'/EPIC_{key}_QPB.fits' # Name of the output file + inargs = [f'table={evtfile}', f'attfile={attfile}', + f'outset={evqpb_outset}', + f'exposurefactor={exposurefactor}',] + w("evqpb", inargs).run() + + gti=f'EPIC_{key}_gti.fit' + test_file(gti) + + expression=f'gti({gti},TIME)' + # filter the EPIC event files to create cleaned and filtered for + # flaring particle background event files for your observation + sci_clean=work_dir+f'/EPIC_{key}_filtered.fits' # Name of the output file + inargs = [f'table={evtfile}', 'withfilteredset=yes', + f'filteredset={sci_clean}', + 'destruct=yes','keepfilteroutput=true', + f'expression={expression}'] + w("evselect", inargs).run() + + # Use the same expression to filter the FWC event files: + qpb_clean=work_dir+f'/EPIC_{key}_QPB_clean.fits' # Name of the output file + inargs = [f'table={evqpb_outset}', 'withfilteredset=yes', + f'filteredset={qpb_clean}', + 'destruct=yes','keepfilteroutput=true', + f'expression={expression}'] + w("evselect", inargs).run() + + for idx, lab in enumerate(label): + if(key=='PN'): + expression=f'(#XMMEA_EP)&&(PATTERN<=4)&&(PI>={pimin[idx]})&&(PI<={pimax[idx]})' + else: + expression=f'(#XMMEA_EM)&&(PATTERN<=12)&&(PI>={pimin[idx]})&&(PI<={pimax[idx]})' + print(lab,expression) + + sci_image=work_dir+f'/EPIC_{key}_sci_image{lab}.fits' # Name of the output file + # Extract an image for the science exposure + inargs = [f'table={sci_clean}', f'expression={expression}', + 'imagebinning=binSize', f'imageset={sci_image}', + 'withimageset=yes','xcolumn=X', 'ycolumn=Y', + 'ximagebinsize=80', 'yimagebinsize=80', + ] + w("evselect", inargs).run() + + qpb_image=work_dir+f'/EPIC_{key}_qpb_image{lab}.fits' # Name of the output file + # Extract an image for the FWC exposure + inargs = [f'table={qpb_clean}', f'expression={expression}', + 'imagebinning=binSize', f'imageset={qpb_image}', + 'withimageset=yes','xcolumn=X', 'ycolumn=Y', + 'ximagebinsize=80', 'yimagebinsize=80','zcolumn=EWEIGHT', + ] + w("evselect", inargs).run() + + #hdul = fits.open(qpb_image) + #primary_hdu = hdul[0] + #qpb_data = hdul[0].data + #qpb_header = hdul[0].header + #hdul.close() + + #result = interpolate_replace_nans(image, kernel) + + # -- Subtract background -- + # skip, due to low statistics + """ + cor_image=work_dir+f'/EPIC_{key}_cor_image{lab}.fits' # Name of the output file + cmd = ['farith', f'{sci_image}', f'{qpb_image}', f'{cor_image}', 'SUB', + 'copyprime=yes', + 'clobber=yes'] + 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): + ref_crd = SkyCoord(arches_ra, arches_dec, frame=FK5(), unit="deg") + nn2=nn*nn + print(outfile_cts) + print(cts[0]) + print(exp[0]) + # take first image as reference + ref_hdul = fits.open(cts[0]) + ref_data = ref_hdul[0].data + ref_header = ref_hdul[0].header + ref_hdul.close() + ref_wcs = WCS(ref_header) + + (nx,ny) = ref_data.shape + map_flx = np.zeros(ref_data.shape,dtype=np.float64) + map_cts = np.zeros(ref_data.shape,dtype=np.float64) + map_qpb = np.zeros(ref_data.shape,dtype=np.float64) + map_exp = np.zeros(ref_data.shape,dtype=np.float64) + map_pix = np.zeros(ref_data.shape,dtype=np.float64) # number of subpixels inside + + # take reference exposure + exp_hdul = fits.open(exp[0]) + exp_data = exp_hdul[0].data + exp_header = exp_hdul[0].header + exp_hdul.close() + exp_wcs = WCS(exp_header) + + """ + for row_index in range(len(ref_data)): + for col_index in range(len(ref_data[row_index])): + refx, refy = ref_wcs.world_to_pixel(ref_crd) + + if(np.absolute(refx-row_index)>cutbox): + continue + if(np.absolute(refy-col_index)>cutbox): + continue + + + sky = ref_wcs.pixel_to_world(row_index, col_index) + + exp_xx, exp_yy = exp_wcs.world_to_pixel(sky) + exp_x = int(np.round(exp_xx)) + exp_y = int(np.round(exp_yy)) + + if not (exp_data[exp_y, exp_x]>0): + continue + map_cts[col_index][row_index]=1.0 + + sep_arcmin = sky.separation(ref_crd).arcmin + if(sep_arcmin > devmax): + continue + map_cts[col_index][row_index]=2.0 + + ref_hdul[0].data=map_cts + ref_hdul.writeto(outfile_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]) + q_data = hdul0[0].data + q_header = hdul0[0].header + q_wcs = WCS(q_header) + hdul0.close() + + hdul0 = fits.open(exp[idx]) + e_data = hdul0[0].data + e_header = hdul0[0].header + e_wcs = WCS(e_header) + hdul0.close() + print(exp[idx]) + + + xx, yy = c_wcs.world_to_pixel(ref_crd) + ref_x = int(np.rint(xx)) + ref_y = int(np.rint(yy)) + + xmin = ref_x - cutbox + xmax = ref_x + cutbox + ymin = ref_y - cutbox + ymax = ref_y + cutbox + + if(xmin<0): + xmin=0 + if(xmax>=nx): + xmax=nx + if(ymin<0): + ymin=0 + if(ymax>=ny): + ymax=ny + + #for row_index in range(len(c_data)): + # for col_index in range(len(c_data[row_index])): + for row_index in range(xmin,xmax): + for col_index in range(ymin,ymax): + c_sky = c_wcs.pixel_to_world(row_index, col_index) + + 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): + continue + + sep_arcmin = c_sky.separation(ref_crd).arcmin + if(sep_arcmin > devmax): + continue + + #xx, yy = ref_wcs.world_to_pixel(c_sky) + #refx = int(np.round(xx)) + #refy = int(np.round(yy)) + + # do subpixeling + for ii in range(nn): + for jj in range(nn): + xx0=np.float64(row_index)-0.5+(np.float64(ii+1)-0.5)/np.float64(nn) + yy0=np.float64(col_index)-0.5+(np.float64(jj+1)-0.5)/np.float64(nn) + subpix_sky = c_wcs.pixel_to_world(xx0, yy0) + + sub_xx, sub_yy = ref_wcs.world_to_pixel(subpix_sky) + refx = np.int32(np.rint(np.float64(sub_xx))) + refy = np.int32(np.rint(np.float64(sub_yy))) + + exp_xx, exp_yy = e_wcs.world_to_pixel(subpix_sky) + exp_x = np.int32(np.rint(np.float64(exp_xx))) + exp_y = np.int32(np.rint(np.float64(exp_yy))) + + if (refx>=0 and refx=0 and refy0): + 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): + ref_hdul[0].data=map_pix + ref_hdul.writeto(outfile_pix,overwrite=True) + + ref_hdul[0].data=map_cts + ref_hdul.writeto(outfile_cts,overwrite=True) + + ref_hdul[0].data=map_exp + ref_hdul.writeto(outfile_exp,overwrite=True) + + ref_hdul[0].data=map_qpb + ref_hdul.writeto(outfile_qpb,overwrite=True) + + ref_hdul[0].data=map_flx + ref_hdul.writeto(outfile_flx,overwrite=True) + #sys.exit() + + """ + cmd = ['farith', f'{outfile_cts}', f'{outfile_qpb}', f'{outfile_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', + 'copyprime=yes', 'clobber=yes'] + result = subprocess.run(cmd, capture_output=True, text=True, check=True) + """ diff --git a/data/FWC/README.md b/data/FWC/README.md new file mode 100644 index 0000000..4b79e46 --- /dev/null +++ b/data/FWC/README.md @@ -0,0 +1,3 @@ +Files taken from here: + +https://www.cosmos.esa.int/web/xmm-newton/filter-closed \ No newline at end of file diff --git a/scripts/06_sas_QPB.py b/scripts/06_sas_QPB.py index d51002f..0db8f6e 100755 --- a/scripts/06_sas_QPB.py +++ b/scripts/06_sas_QPB.py @@ -51,83 +51,24 @@ for obsid in files: # Create QPB # https://www.cosmos.esa.int/web/xmm-newton/sas-thread-background # - - attfiles = glob.glob(events_dir+f'/{obsid}/*{obsid}_AttHk.ds') - if(attfiles): - print(f"*** {attfiles} ***") - else: - print("AttHk?") - sys.exit() - attfile = attfiles[0] - if not os.path.isfile(attfile): - print ("File "+attfile+" does not exist, please check. \n") - sys.exit() + attfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_AttHk.ds') - spec=[] - bkg=[] - resp=[] - arf=[] #for imos in [1,2]: for key in ['MOS1','MOS2','PN']: evtfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_E{key}*_ImagingEvts.ds') - # Run the SAS task evqpb over each one of these event files. - evqpb_outset=work_dir+f'/EPIC_{key}_QPB.fits' # Name of the output file - inargs = [f'table={evtfile}', f'attfile={attfile}', - f'outset={evqpb_outset}', - f'exposurefactor=2.0',] - w("evqpb", inargs).run() - - gti=f'EPIC_{key}_gti.fit' - expression=f'gti({gti},TIME)' - # filter the EPIC event files to create cleaned and filtered for - # flaring particle background event files for your observation - sci_clean=work_dir+f'/EPIC_{key}_filtered.fits' # Name of the output file - inargs = [f'table={evtfile}', 'withfilteredset=yes', - f'filteredset={sci_clean}', - 'destruct=yes','keepfilteroutput=true', - f'expression={expression}'] - w("evselect", inargs).run() - - # Use the same expression to filter the FWC event files: - qpb_clean=work_dir+f'/EPIC_{key}_QPB_clean.fits' # Name of the output file - inargs = [f'table={evqpb_outset}', 'withfilteredset=yes', - f'filteredset={qpb_clean}', - 'destruct=yes','keepfilteroutput=true', - f'expression={expression}'] - w("evselect", inargs).run() + # first element contains all energies + label=['',] + pimin=[200,] + pimax=[12000,] + for idx,e1 in enumerate(emin): + label.append(f"_{idx}") + pimin.append(emin[idx]) + pimax.append(emax[idx]) + + run_evqpb(key, work_dir=work_dir, evtfile=evtfile, attfile=attfile, label=label, pimin=pimin, pimax=pimax) - if(key=='PN'): - expression='(#XMMEA_EP)&&(PATTERN<=4)&&(PI>=200)&&(PI<=12000)' - else: - expression='(#XMMEA_EM)&&(PATTERN<=12)&&(PI>=200)&&(PI<=12000)' - - sci_image=work_dir+f'/EPIC_{key}_sci_image.fits' # Name of the output file - # Extract an image for the science exposure - inargs = [f'table={sci_clean}', f'expression={expression}', - 'imagebinning=binSize', f'imageset={sci_image}', - 'withimageset=yes','xcolumn=X', 'ycolumn=Y', - 'ximagebinsize=80', 'yimagebinsize=80', - ] - w("evselect", inargs).run() - - qpb_image=work_dir+f'/EPIC_{key}_qpb_image.fits' # Name of the output file - # Extract an image for the FWC exposure - inargs = [f'table={qpb_clean}', f'expression={expression}', - 'imagebinning=binSize', f'imageset={qpb_image}', - 'withimageset=yes','xcolumn=X', 'ycolumn=Y', - 'ximagebinsize=80', 'yimagebinsize=80','zcolumn=EWEIGHT', - ] - w("evselect", inargs).run() - - # Subtract background - cor_image=work_dir+f'/EPIC_{key}_cor_image.fits' # Name of the output file - cmd = ['farith', f'{sci_image}', f'{qpb_image}', f'{cor_image}', 'SUB', - 'copyprime=yes', - 'clobber=yes'] - result = subprocess.run(cmd, capture_output=True, text=True, check=True) - continue # diff --git a/scripts/07_sas_eimageget.py b/scripts/07_sas_eimageget.py new file mode 100755 index 0000000..3b73579 --- /dev/null +++ b/scripts/07_sas_eimageget.py @@ -0,0 +1,77 @@ +#!/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 + +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)) + +fwc_dir=root_path+'/data/FWC' +archive_dir=root_path+'/data/archive' +events_dir=root_path+'/data/processed' +events_dir_oot=root_path+'/data/processed-oot' +products_dir=root_path+'/products/sas' +ds9reg_dir=root_path+'/data/ds9reg' + + +create_folder(products_dir) + +inargs = ['--version'] +t = w('sasver', inargs) +t.run() + +files = glob.glob(archive_dir+'/*') + +for obsid in files: + obsid = os.path.basename(obsid) + if(obsid in skip): + continue + + work_dir = init_work_dir(obsid, products_dir=products_dir) + os.chdir(work_dir) + + # + # https://www.cosmos.esa.int/web/xmm-newton/sas-thread-images + # + + attfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_AttHk.ds') + + for key in ['MOS1','MOS2','PN']: + evtfile = get_first_file(events_dir+f'/{obsid}/*{obsid}_E{key}*_ImagingEvts.ds') + fwcfile = get_first_file(fwc_dir+f'/{key}_closed_FF_2025_v1.fits') + pimin=" ".join(emin) + pimax=" ".join(emax) + gti=f'EPIC_{key}_gti.fit' + sci_clean=work_dir+f'/EPIC_{key}_filtered.fits' # Name of the output file + inargs = [f'evtfile={evtfile}','withemtaglenoise=no', + f'fwcfile={fwcfile}','withfwcimages=yes', + f'gtifile={gti}', + f'attfile={attfile}', + f'pimin={pimin}', + f'pimax={pimax}', + ] + if(key=='PN'): + evtfile_oot = get_first_file(events_dir_oot+f'/{obsid}/*{obsid}_E{key}*_ImagingEvts.ds') + inargs.append(f'ootfile={evtfile_oot}') + w("eimageget", inargs).run() + + #sys.exit() diff --git a/scripts/08_sas_mosaic.py b/scripts/08_sas_mosaic.py new file mode 100755 index 0000000..6f0368f --- /dev/null +++ b/scripts/08_sas_mosaic.py @@ -0,0 +1,149 @@ +#!/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 + +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) +mosaic_dir=root_path+'/products/mosaic' +create_folder(mosaic_dir) + +ds9reg_dir=root_path+'/data/ds9reg' + + +inargs = ['--version'] +t = w('sasver', inargs) +t.run() + +files = glob.glob(archive_dir+'/*') + + +# run over energies +for idx,e in enumerate(emin): + label=f"_{idx}" + e1=emin[idx] + e2=emax[idx] + + # all MOS1 + m1_ima=[] + m1_qpb=[] + m1_exp=[] + + # all MOS2 + m2_ima=[] + m2_qpb=[] + m2_exp=[] + + # all MOS + mos_ima=[] + mos_qpb=[] + mos_exp=[] + + # all PNs + pn_ima=[] + pn_qpb=[] + pn_exp=[] + + # all MOS1,2,PN + ima=[] + qpb=[] + exp=[] + + for obsid in files: + obsid = os.path.basename(obsid) + if(obsid in skip): + continue + + work_dir = f'{products_dir}/{obsid}' + + for ikey,key in enumerate(['MOS1','MOS2','PN']): + shortkey = key.replace("OS","") + expkey=f'S00{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) + qpb_image=work_dir+f'/EPIC_{key}_qpb_image{label}.fits' + test_file(qpb_image) + # 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) + + ima.append(sci_image) + qpb.append(qpb_image) + exp.append(exp_image) + + if(key=='PN'): + pn_ima.append(sci_image) + pn_qpb.append(qpb_image) + pn_exp.append(exp_image) + if(key=='MOS1'): + m1_ima.append(sci_image) + m1_qpb.append(qpb_image) + m1_exp.append(exp_image) + if(key=='MOS2'): + m2_ima.append(sci_image) + m2_qpb.append(qpb_image) + m2_exp.append(exp_image) + if(key=='MOS1' or key=='MOS2'): + mos_ima.append(sci_image) + mos_qpb.append(qpb_image) + mos_exp.append(exp_image) + + + 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) +