### Exemplary CIAO script to process data of galaxy NGC4649 (M60)
### I tested all of this but there can always be bugs, so be careful and let me know if you find any

# # I suggest you comment out all the code now and then slowly uncomment step by step (note: for loops need to be finished by "done")
# # Also its probably better to start with just one OBSID and then add more into the IDS list

# # command that finds all obsids of a specific object and saves them into a file
find_chandra_obsid NGC4649 > obsids.txt

# # read the file and choose OBSIDs you want to analyze
# # in the first run, start with just one OBSID and look at the directories and files
IDS=(785)


# # Now we will loop over all OBSIDs and download and reprocess the data
for i in ${IDS[@]}; do
echo $i
# # donwload specific Chandra observation
download_chandra_obsid $i

# # # reprocess the data
chandra_repro $i $i/repro
done


# # reproject and merge all observations together, to include all observation use regular expression "*/repro/*evt2.fits"
# # this is update compared to the PDF - merge_obs is a script combining reproject_obs and flux_obs scripts
merge_obs */repro/*evt2.fits reproj/ bin=1 bands=broad psfecf=0.9 clobber=yes

# # detect point sources in the reprojected, merged image using also the exposure map and PSF map
# # many of the weird file are just output files, we care about the wav.reg file which contains regions with detected sources
# # warning: wavdetect might give errors on CIAO-4.16 -> use CIAO-4.15 either just for wavdetect or for the whole script
wavdetect infile="reproj/broad_thresh.img" expfile="reproj/broad_thresh.expmap" regfile=reproj/wav.reg outfile=reproj/wav.fits scellfile=reproj/wav_scell.fits imagefile=reproj/wav_img.fits defnbkgfile=reproj/wav_bkg.fits scales="2 4 8" psffile="reproj/broad_thresh.psfmap" clobber=yes
# # now look at the point sources in DS9: $ ds9 reproj/broad_thresh.img -region wav.reg
# # correct wrongly detected or add non-detected point sources (use circles or ellipses)
# # then save the sources into a different file e.g. sources.reg in the CIAO > WCS > FK5 format


for i in ${IDS[@]}; do
echo $i
# # extract background lightcurve excluding the point sources
dmextract "reproj/${i}_reproj_evt.fits[exclude sky=region(sources.reg)][bin time=::137.12]" outfile="${i}/repro/${i}.lc" opt=ltc1 clobber=yes

# # run deflaring algorithm to remove bins higher than 3 sigma -> produces Good Time Intervals (GTI) file
deflare "${i}/repro/${i}.lc" "${i}/repro/${i}.gti" plot=no method=sigma stddev=3.0 save="${i}/repro/${i}_deflared"

# # filter reprojected event files using GTI file and save it as _clean.fits
dmcopy "reproj/${i}_reproj_evt.fits[@${i}/repro/${i}.gti]" "${i}/repro/${i}_clean.fits" clobber=yes
done


# # merge all cleaned event files
merge_obs */repro/*_clean.fits merged/ bin=1 bands=broad clobber=yes

# # now you may want to fill the point sources in the images 
# # -> before doing that, exclude the central PS (or any other source you want) from sources.reg since you usually dont want to fill that one
# # then run the create_fill_backgrounds.py script to create backgrounds (you can play with a,b constants)
python3 create_fill_backgrounds.py
# # after you get the sources_bkg.reg file, run the following command
dmfilth "merged/broad_thresh.img" "broad_thresh_filled.fits" srclist="@sources.reg" bkglist="@sources_bkg.reg" method="POISSON"
# # now check the broad_thresh_filled.fits image and probably iterate (update sources.reg) because some PSs will still remain
# # smoothing the image in DS9 will also help you visually identify sources


for i in ${IDS[@]}; do
echo $i
# # create a background event file for spectral extraction (based on CalDB files)
blanksky evtfile="${i}/repro/${i}_clean.fits" outfile="${i}/repro/${i}_bkg.fits" clobber=yes

# # for spectral extraction, use DS9 to create a region file for the source (src.reg) and background (bkg.reg), save in CIAO > WCS > FK5 format
specextract "${i}/repro/${i}_clean.fits[sky=region(src.reg)]" spectra/spec bkgfile="${i}/repro/${i}_clean.fits[sky=region(bkg.reg)]" clobber=yes
# # if source is extedned, you can use blanksky background event file and extract background from the same region as the source (src.reg)
specextract "${i}/repro/${i}_clean.fits[sky=region(src.reg)]" spectra_blanksky/spec bkgfile="${i}/repro/${i}_bkg.fits[sky=region(src.reg)]" bkgresp=no clobber=yes
# # if the background is extracted from different part of the chip (bkg.reg), ARF and RMF files are also generated for background because it's different part of chip
# # if the background spectrum is extracted from blanksky bkg file, we use the same ARF and RMF as for the source -> bkgresp=no
done
