## for XCalibur meeting
## Hitomi Analysis script
## version 1 (2019.6.25)
## Misaki Mizumoto (misaki.mizumoto@durham.ac.uk)


#####################
## (0) Download
#####################
## from HEASARC (777MB)
wget -q -nH --no-check-certificate --cut-dirs=5 -r -l0 -c -N -np -R 'index*' -erobots=off --retr-symlinks https://heasarc.gsfc.nasa.gov/FTP/hitomi/data/obs/1//100040030/sxs/
wget -q -nH --no-check-certificate --cut-dirs=5 -r -l0 -c -N -np -R 'index*' -erobots=off --retr-symlinks https://heasarc.gsfc.nasa.gov/FTP/hitomi/data/obs/1//100040030/auxil/

## NXB files (218MB)
mkdir NXB
cd NXB
wget https://darts.jaxa.jp/pub/hitomi/data/nxb_20170510/ah_sxs_nxbafmar4_20140101v001.evt.gz 
wget https://darts.jaxa.jp/pub/hitomi/data/nxb_20170510/ah_gen_nxbehk_20140101v002.fits.gz 
cd ..

#####################
## (1) Screening
#####################
mkdir products_sxs 
cd products_sxs

ftselect infile='../100040030/sxs/event_cl/ah100040030sxs_p0px1010_cl.evt.gz[events]' \
outfile=ah100040030sxs_p0px1010_cl2.evt \
expression="(PI>=400)&&((RISE_TIME>=40&&RISE_TIME<=60&&ITYPE<4)||(ITYPE==4))&&STATUS[4]==b0"

#####################
## (2) Extract image
#####################
fkeyprint ah100040030sxs_p0px1010_cl2.evt+1 RA_NOM
## RA_NOM  =     49.9324156058081 / [deg] Nominal aspect point R.A.
fkeyprint ah100040030sxs_p0px1010_cl2.evt+1 DEC_NOM
## DEC_NOM =     41.5201310551332 / [deg] Nominal aspect point Dec.
fkeyprint ah100040030sxs_p0px1010_cl2.evt+1 PA_NOM
## PA_NOM  =      253.06395107642 / [deg] Nominal aspect point roll.
## 253.06395107642-360=-106.936048924

ra=49.9324156058081
dec=41.5201310551332
roll=-106.936048924

ahmkregion instrume=SXS ra=${ra} dec=${dec} roll=${roll}

## extract image
xselect
xsel
read events ah100040030sxs_p0px1010_cl2.evt
./
yes
extract image
save image ah100040030sxs_p0px1010_cl2_sky.img 
set xyname detx dety
extract image
save image ah100040030sxs_p0px1010_cl2_det.img
exit
no

## check the image
ds9 ah100040030sxs_p0px1010_cl2_det.img -region file SXS.DET.box.reg -region file SXS.DET.text.reg -zoom to fit -cmap bb &

## center point of the detector
coordpnt "3.5,3.5" outfile=none telescop=HITOMI instrume=SXS ra=${ra} dec=${dec} roll=${roll} startsys=det stopsys=radec
## coordpnt: OUTX OUTY=   49.93350633    41.52159242

ra_center=49.93350633
dec_center=41.52159242

cat << EOF >! perseus_offset_sky.reg
fk5
circle(${ra_center},${dec_center},150.0")
EOF

## check the image
ds9 ah100040030sxs_p0px1010_cl2_sky.img -region file perseus_offset_sky.reg -cmap bb &

#####################
## (3) Extract full-array data
#####################

## 35-pixel detector region and exposure map
sxsregext infile=ah100040030sxs_p0px1010_cl2.evt \
regmode=RADEC region=./perseus_offset_sky.reg \
resolist=0 outroot=ah100040030sxs_detreg \
outexp=ah100040030sxs.expo \
ehkfile=../100040030/auxil/ah100040030.ehk.gz \
delta=20 numphi=1 clobber=yes

## check image
ds9 ah100040030sxs_detreg.img -region file SXS.DET.box.reg -region file SXS.DET.text.reg -zoom to fit -cmap bb &

## check light curve
lcurve 
1
ah100040030sxs_detreg.lc
-
INDEF
INDEF
default
yes
/XW
quit

## extract spectra
xselect
xsel
read events ah100040030sxs_p0px1010_cl2.evt
./
yes
filter column "PIXEL=0:11,13:35" 
filter GRADE "0:0"
extract spectrum
save spectrum ah100040030sxs_p0px1010_cl2_HP.pi
exit
no


#####################
## (4) Extract NXB spectrum
#####################

sxsnxbgen infile=ah100040030sxs_p0px1010_cl2.evt \
ehkfile=../100040030/auxil/ah100040030.ehk.gz regfile=NONE \
innxbfile=../NXB/ah_sxs_nxbafmar4_20140101v001.evt.gz \
innxbehk=../NXB/ah_gen_nxbehk_20140101v002.fits.gz \
outpifile=ah100040030sxsnxb_cl2.pi pixels="-" \
cleanup=yes chatter=3 clobber=yes mode=hl \
logfile=ah100040030sxsnxb_cl2.log \
sortbin=0,4,5,6,7,8,9,10,11,12,13,99 \
expr="PI>=400&&RISE_TIME>=40&&RISE_TIME<=60&&ITYPE<4&&STATUS[4]==b0"

fthedit ah100040030sxsnxb_cl2.pi+1 BACKSCAL add 1.000000E+00


#####################
## (5) Generating response files
#####################

## (5.1) RMF file
sxsmkrmf infile=ah100040030sxs_p0px1010_cl2.evt \
outfile=ah100040030_sxs_cl2_HP_small.rmf resolist=0 \
regmode=det regionfile=ah100040030sxs_detreg.reg whichrmf=m

## (5.2) recalculate exposure map
## to assure that map includes only a single attitude bin

ahexpmap ehkfile=../100040030/auxil/ah100040030.ehk.gz \
gtifile=ah100040030sxs_p0px1010_cl2.evt instrume=SXS badimgfile=NONE \
pixgtifile=../100040030/sxs/event_uf/ah100040030sxs_px1010_exp.gti.gz \
outfile=ah100040030sxs_p0px1010.expo outmaptype=EXPOSURE \
delta=20.0 numphi=1 stopsys=SKY instmap=CALDB qefile=CALDB \
contamifile=CALDB vigfile=CALDB obffile=CALDB fwfile=CALDB \
gvfile=CALDB maskcalsrc=yes fwtype=DEFAULT specmode=MONO \
specfile=spec.fits specform=FITS energy=1.5 evperchan=DEFAULT \
abund=1 cols=0 covfac=1 clobber=yes chatter=1 \
logfile=make_expo_ah100040030sxs_p0px1010.log


## (5.3) ARF file
## You can choose the assumption about PSF; point source, extended source (beta model), and extended source (XMM image)
## Here the point source version is adopted.
## Calculation time: ~30 min (depending on your machine)

aharfgen xrtevtfile=raytrace_ah100040030sxs_p0px1010_ptsrc.fits \
source_ra=49.9507 source_dec=41.5128 telescop=HITOMI \
instrume=SXS emapfile=ah100040030sxs_p0px1010.expo \
regmode=DET regionfile=ah100040030sxs_detreg.reg \
sourcetype=POINT rmffile=ah100040030_sxs_cl2_HP_small.rmf \
erange="0.5 17.0 0 0" outfile=ah100040030sxs_p0px1010_ptsrc.arf \
numphoton=300000 minphoton=1 teldeffile=CALDB qefile=CALDB \
contamifile=CALDB obffile=CALDB fwfile=CALDB gatevalvefile=CALDB \
onaxisffile=CALDB onaxiscfile=CALDB mirrorfile=CALDB \
obstructfile=CALDB frontreffile=CALDB backreffile=CALDB \
pcolreffile=CALDB scatterfile=CALDB auxtransfile=NONE mode=h \
clobber=yes seed=7 logfile=make_arf_ah100040030sxs_p0px1010_ptsrc.log

#####################
## (6) Spectral fitting
#####################

## xspec script
xspec

statistic cstat
abund lodd
xsect bcmc

data 1:1 ah100040030sxs_p0px1010_cl2_HP.pi
resp 1 ah100040030_sxs_cl2_HP_small.rmf
arf 1 ah100040030sxs_p0px1010_ptsrc.arf
ig **-1.8,9.0-**
setplot rebin 1e6 4
pl ld

model tbabs(bapec+pegpwrlw)
0.138 -1
3.96
0.63
1.75e-2 0.01
185 0.01
0.162 0.01
1.8 -1
2 -1
10 -1
30
query yes
fit
pl ld ra

ipl ld ra
win all
log x off
r x 6.4 6.7
r y 0.09 4