aboutsummaryrefslogtreecommitdiffstats
path: root/scripts/ciao_calc_ct.sh
blob: 2a277e999aca307e40ca2e03a0c72a6f8f70fc7d (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
#!/bin/sh
#
## Calculate the cooling time within the (0.-0.048 r500) region.
## Based on `ciao_r500avgt.sh'
##
## Junhua GU
## Created: 2012-08-22
## Weitian LI
## Updated: 2016-06-08
##
## Change logs:
## 2016-06-08, Weitian LI
##   * Drop 'calc_distance' in favor of 'cosmo_calc'
## v2.0, 2015/06/03, Weitian LI
##   * Copy needed pfiles to current working directory, and
##     set environment variable $PFILES to use these first.
##   * Replace 'grep' with '\grep', 'ls' with '\ls'
##   * replaced 'grppha' with 'dmgroup' to group spectra
##     (dmgroup will add history to fits file, while grppha NOT)
##   * ds9 colormap changed from 'sls' to 'he'
## v1.2, 2015/05/27, Weitian LI
##   update 'DFT_ARF' & 'DFT_RMF' to find '*.arf' & '*.rmf' files
##   (specextract only use .arf & .rmf extensions since revision 2014-12)
## v1.1, 2014/06/18
##   'cooling_time2.sh' -> 'ciao_calc_ct.sh'
##

## error code {{{
ERR_USG=1
ERR_DIR=11
ERR_EVT=12
ERR_BKG=13
ERR_REG=14
ERR_ASOL=21
ERR_BPIX=22
ERR_PBK=23
ERR_MSK=24
ERR_BKGTY=31
ERR_SPEC=32
ERR_ARF=51
ERR_RMF=52
ERR_UNI=61
## error code }}}

## usage, help {{{
case "$1" in
    -[hH]*|--[hH]*)
        printf "usage:\n"
        printf "    `basename $0` basedir=<repro_dir> evt=<evt2_clean> r500=<r500_kpc> regin=<input_reg> regout=<output_reg> bkgd=<blank_evt | lbkg_reg | bkg_spec> nh=<nH> z=<redshift> arf=<arf_file> rmf=<rmf_file> [ grouptype=<NUM_CTS|BIN> grouptypeval=<number> binspec=<binspec> log=<log_file> ]\n"
        exit ${ERR_USG}
        ;;
esac
## usage, help }}}

## default parameters {{{
# default `event file' which used to match `blanksky' files
#DFT_EVT="_NOT_EXIST_"
DFT_EVT="`\ls evt2*_clean.fits 2> /dev/null`"
# default `bkgd', use `bkgcorr_blanksky*' corrected bkg spectrum
DFT_BKGD="`\ls bkgcorr_blanksky_*.pi 2> /dev/null`"
# default basedir
DFT_BASEDIR="../.."
# default info_json pattern
DFT_JSON_PAT="*_INFO.json"
# default `radial region file'
#DFT_REG_IN="_NOT_EXIST_"
DFT_REG_IN="rspec.reg"
# default output region file (0.1-0.5 r500 region)
DFT_REG_OUT="cooling.reg"
# default ARF/RMF, the one of the outmost region
DFT_ARF="`\ls r1_*.warf r1_*.arf 2> /dev/null`"
DFT_RMF="`\ls r1_*.wrmf r1_*.rmf 2> /dev/null`"

# default parameters for 'dmgroup'
DFT_GROUPTYPE="NUM_CTS"
DFT_GROUPTYPEVAL="20"
#DFT_GROUPTYPE="BIN"
DFT_BINSPEC="1:128:2,129:256:4,257:512:8,513:1024:16"

INNER=0.0
OUTER=0.048
CT_RES="cooling_results.txt"
# default output xspec scripts
XSPEC_SCRIPT="xspec_cooling.xcm"
# deproj xspec script, generated by `deproj_spectra'
# from which get `nh' and `redshift'
XSPEC_DEPROJ="xspec_deproj.xcm"

# default `log file'
DFT_LOGFILE="cooling_`date '+%Y%m%d'`.log"
## default parameters }}}

## functions {{{
# process commandline arguments
# cmdline arg format: `KEY=VALUE'
getopt_keyval() {
    until [ -z "$1" ]
    do
        key=${1%%=*}                    # extract key
        val=${1#*=}                     # extract value
        keyval="${key}=\"${val}\""
        echo "## getopt: eval '${keyval}'"
        eval ${keyval}
        shift                           # shift, process next one
    done
}

## background renormalization (BACKSCAL) {{{
# renorm background according to particle background
# energy range: 9.5-12.0 keV (channel: 651-822)
CH_LOW=651
CH_HI=822
pb_flux() {
    punlearn dmstat
    COUNTS=`dmstat "$1[channel=${CH_LOW}:${CH_HI}][cols COUNTS]" | \grep -i 'sum:' | awk '{ print $2 }'`
    punlearn dmkeypar
    EXPTIME=`dmkeypar $1 EXPOSURE echo=yes`
    BACK=`dmkeypar $1 BACKSCAL echo=yes`
    # fix `scientific notation' bug for `bc'
    EXPTIME_B=`echo ${EXPTIME} | sed 's/[eE]/\*10\^/' | sed 's/+//'`
    BACK_B=`echo "( ${BACK} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'`
    PB_FLUX=`echo "scale = 16; ${COUNTS} / ${EXPTIME_B} / ${BACK_B}" | bc -l`
    echo ${PB_FLUX}
}

bkg_renorm() {
    # $1: src spectrum, $2: back spectrum
    PBFLUX_SRC=`pb_flux $1`
    PBFLUX_BKG=`pb_flux $2`
    BACK_OLD=`dmkeypar $2 BACKSCAL echo=yes`
    BACK_OLD_B=`echo "( ${BACK_OLD} )" | sed 's/[eE]/\*10\^/' | sed 's/+//'`
    BACK_NEW=`echo "scale = 16; ${BACK_OLD_B} * ${PBFLUX_BKG} / ${PBFLUX_SRC}" | bc -l`
    printf "\`$2': BACKSCAL:\n"
    printf "    ${BACK_OLD} --> ${BACK_NEW}\n"
    punlearn dmhedit
    dmhedit infile=$2 filelist=none operation=add \
        key=BACKSCAL value=${BACK_NEW} comment="old value: ${BACK_OLD}"
}
## bkg renorm }}}
## functions end }}}

## parameters {{{
# process cmdline args using `getopt_keyval'
getopt_keyval "$@"

## check log parameters {{{
if [ ! -z "${log}" ]; then
    LOGFILE="${log}"
else
    LOGFILE=${DFT_LOGFILE}
fi
printf "## use logfile: \`${LOGFILE}'\n"
[ -e "${LOGFILE}" ] && mv -fv ${LOGFILE} ${LOGFILE}_bak
TOLOG="tee -a ${LOGFILE}"
echo "process script: `basename $0`" >> ${LOGFILE}
echo "process date: `date`" >> ${LOGFILE}
## log }}}

# check given parameters
# process `nh' and `redshift' {{{
if [ ! -r "${XSPEC_DEPROJ}" ]; then
    printf "## the script \`${XSPEC_DEPROJ}' generated by \`deproj_spectra' NOT found\n" | ${TOLOG}
    printf "## please supply the value of \`nh' and \`redshift' manully\n" | ${TOLOG}
    read -p "> value of nh: " N_H
    read -p "> value of redshift: " REDSHIFT
else
    # get `nh' and `redshift' from xspec script
    LN=`\grep -n 'projct\*wabs\*apec' ${XSPEC_DEPROJ} | tail -n 1 | cut -d':' -f1`
    # calc the line number of which contains `nh'
    LN_NH=`expr ${LN} + 4`
    NH_XCM=`head -n ${LN_NH} ${XSPEC_DEPROJ} | tail -n 1 | awk '{ print $1 }'`
    # calc the line number of `redshift'
    LN_RS=`expr ${LN} + 7`
    RS_XCM=`head -n ${LN_RS} ${XSPEC_DEPROJ} | tail -n 1 | awk '{ print $1 }'`
    printf "## get value of nh: \`${NH_XCM}' (from \`${XSPEC_DEPROJ}')\n" | ${TOLOG}
    printf "## get value of redshift: \`${RS_XCM}' (from \`${XSPEC_DEPROJ}')\n" | ${TOLOG}

    ## if `nh' and `redshift' supplied in cmdline, then use them
    if [ ! -z "${nh}" ]; then
        N_H=${nh}
    else
        N_H=${NH_XCM}
    fi
    # redshift
    if  [ ! -z "${z}" ]; then
        REDSHIFT=${z}
    else
        REDSHIFT=${RS_XCM}
    fi
fi
printf "## use nH: ${N_H}\n" | ${TOLOG}
printf "## use redshift: ${REDSHIFT}\n" | ${TOLOG}
# nh & redshift }}}

# check basedir & json file
if [ -d "${basedir}" ]; then
    BASEDIR=${basedir}
else
    BASEDIR=${DFT_BASEDIR}
fi
# json file
if [ ! -z "${json}" ] && [ -r "${BASEDIR}/${json}" ]; then
    JSON_FILE="${BASEDIR}/${json}"
elif [ `\ls ${BASEDIR}/${DFT_JSON_PAT} 2> /dev/null | wc -l` -eq 1 ]; then
    JSON_FILE=`\ls ${BASEDIR}/${DFT_JSON_PAT}`
else
    read -p "> JSON_file: " JSON_FILE
    if [ ! -r "${JSON_FILE}" ]; then
        printf "ERROR: cannot access given \`${JSON_FILE}'\n"
        exit ${ERR_JSON}
    fi
fi
printf "## use json_file: \`${JSON_FILE}'\n" | ${TOLOG}

# process `r500' {{{
R500_RAW=`\grep '"R500.*kpc' ${JSON_FILE} | sed 's/.*"R500.*":\ //' | sed 's/\ *,$//'`
if [ ! -z "${r500}" ]; then
    R500_RAW=${r500}
fi
if [ -z "${R500_RAW}" ]; then
    printf "## input R500 followed with unit, e.g.: 800kpc, 400pix\n"
    read -p "> value of \`R500' (in pixel/kpc): " R500_RAW
fi
R500_VAL=`echo "${R500_RAW}" | tr -d 'a-zA-Z, '`
R500_UNI=`echo "${R500_RAW}" | tr -d '0-9, '`
printf "## get \`R500': ${R500_VAL} in unit \`${R500_UNI}'\n" | ${TOLOG}

# if in kpc, convert to pix
case "${R500_UNI}" in
    [pP]*)
        printf "## units in \`pixel', conversion not needed\n" | ${TOLOG}
        R500_PIX_B=`echo ${R500_VAL} | sed 's/[eE]/\*10\^/' | sed 's/+//'`
        ;;
    *)
        printf "## units in \`kpc', convert to \`Chandra pixel'\n" | ${TOLOG}
        KPC_PER_PIX=`cosmo_calc ${REDSHIFT} | \grep 'kpc.*pix' | tr -d 'a-zA-Z_#=(),:/ '`
        # convert scientific notation for `bc'
        KPC_PER_PIX_B=`echo ${KPC_PER_PIX} | sed 's/[eE]/\*10\^/' | sed 's/+//'`
        printf "## calculated \`kpc/pixel': ${KPC_PER_PIX_B}\n"
        R500_VAL_B=`echo ${R500_VAL} | sed 's/[eE]/\*10\^/' | sed 's/+//'`
        R500_PIX_B=`echo "scale = 4; ${R500_VAL_B} / ( ${KPC_PER_PIX_B} )" | bc -l`
        ;;
esac
# calc (inner-outer R500)
R_IN=`echo "scale = 4; ${INNER} * ${R500_PIX_B}" | bc -l`
R_OUT=`echo "scale = 4; ${OUTER} * ${R500_PIX_B}" | bc -l`
printf "## R500 in units pixel: ${R500_PIX_B}\n" | ${TOLOG}
printf "## (${INNER}-${OUTER} R500) range in pixel: ${R_IN} - ${R_OUT}\n" | ${TOLOG}
# r500 }}}

# check evt file
if [ -r "${evt}" ]; then
    EVT=${evt}
elif [ -r "${DFT_EVT}" ]; then
    EVT=${DFT_EVT}
else
    read -p "> clean evt2 file: " EVT
    if [ ! -r "${EVT}" ]; then
        printf "ERROR: cannot access given \`${EVT}' evt file\n"
        exit ${ERR_EVT}
    fi
fi
printf "## use evt file: \`${EVT}'\n" | ${TOLOG}

# input and output region files {{{
if [ -r "${regin}" ]; then
    REG_IN="${regin}"
elif [ -r "${DFT_REG_IN}" ]; then
    REG_IN=${DFT_REG_IN}
else
    read -p "> previous used radial spec regfile: " REG_IN
    if [ ! -r "${REG_IN}" ]; then
        printf "ERROR: cannot access given \`${REG_IN}' region file\n"
        exit ${ERR_REG}
    fi
fi
printf "## use previous regfile: \`${REG_IN}'\n" | ${TOLOG}
if [ ! -z "${regout}" ]; then
    REG_OUT="${regout}"
else
    REG_OUT=${DFT_REG_OUT}
fi
[ -e "${REG_OUT}" ] && mv -fv ${REG_OUT} ${REG_OUT}_bak
printf "## set output regfile (0.1-0.5 r500 region): \`${REG_OUT}'\n" | ${TOLOG}

# get center position from `regin'
# only consider `pie' or `annulus'-shaped region
TMP_REG=`\grep -iE '(pie|annulus)' ${REG_IN} | head -n 1`
XC=`echo ${TMP_REG} | tr -d ' ' | awk -F'[(),]' '{ print $2 }'`
YC=`echo ${TMP_REG} | tr -d ' ' | awk -F'[(),]' '{ print $3 }'`
printf "## get center coord: (${XC},${YC})\n" | ${TOLOG}
# region files }}}

# check given bkgd, determine background {{{
if [ -r "${bkgd}" ]; then
    BKGD=${bkgd}
elif [ -r "${DFT_BKGD}" ]; then
    BKGD=${DFT_BKGD}
else
    read -p "> background (blanksky_evt | lbkg_reg | bkg_spec): " BKGD
    if [ ! -r "${BKGD}" ]; then
        printf "ERROR: cannot access given \`${BKGD}'\n"
        exit ${ERR_BKG}
    fi
fi
printf "## use bkgd: \`${BKGD}'\n" | ${TOLOG}
# determine bkg type: blanksky, lbkg_reg, bkg_spec ?
# according to file type first: text / FITS
# if FITS, then get values of `HDUCLAS1' and `OBJECT'
if file -bL ${BKGD} | \grep -qi 'text'; then
    printf "## given \`${BKGD}' is a \`text file'\n"
    printf "##   use it as local bkg region file\n"
    printf "##   use *LOCAL BKG SPEC*\n" | ${TOLOG}
    # just set flags, extract spectrum later
    USE_LBKG_REG=YES
    USE_BLANKSKY=NO
    USE_BKG_SPEC=NO
elif file -bL ${BKGD} | \grep -qi 'FITS'; then
    printf "## given \`${BKGD}' is a \`FITS file'\n"
    # get FITS header keyword
    punlearn dmkeypar
    HDUCLAS1=`dmkeypar ${BKGD} HDUCLAS1 echo=yes`
    if [ "${HDUCLAS1}" = "EVENTS" ]; then
        # event file
        printf "##   given file is \`event'\n"
        # check if `blanksky' or `stowed bkg'
        BKG_OBJ=`dmkeypar ${BKGD} OBJECT echo=yes`
        if [ "${BKG_OBJ}" = "BACKGROUND DATASET" ] || [ "${BKG_OBJ}" = "ACIS STOWED" ]; then
            # valid bkg evt file
            printf "##   given FITS file is a valid bkgrnd file\n"
            printf "##   use *BLANKSKY*\n" | ${TOLOG}
            USE_BLANKSKY=YES
            USE_LBKG_REG=NO
            USE_BKG_SPEC=NO
            # specify `BLANKSKY'
            BLANKSKY=${BKGD}
        else
            # invalid bkg evt file
            printf "ERROR: invalid bkg evt file given\n"
            exit ${ERR_BKGTY}
        fi
    elif [ "${HDUCLAS1}" = "SPECTRUM" ]; then
        # spectrum file
        printf "##   given file is \`spectrum'\n"
        printf "##   use *BKG SPECTRUM*\n" | ${TOLOG}
        USE_BKG_SPEC=YES
        USE_BLANKSKY=NO
        USE_LBKG_REG=NO
        # specify `BKG_SPEC'
        BKG_SPEC=${BKGD}
    else
        # other type
        printf "ERROR: other type FITS given\n"
        exit ${ERR_BKGTY}
    fi
else
    printf "ERROR: given \`${BKGD}' type UNKNOWN\n"
    exit ${ERR_BKGTY}
fi
# bkgd }}}

# check `arf' and `rmf' {{{
if [ -r "${arf}" ]; then
    ARF=${arf}
elif [ -r "${DFT_ARF}" ]; then
    ARF=${DFT_ARF}
else
    read -p "> provide the ARF to use: " ARF
    if [ ! -r "${ARF}" ]; then
        printf "ERROR: cannot access given \`${ARF}'\n"
        exit ${ERR_ARF}
    fi
fi
printf "## use ARF: \`${ARF}'\n" | ${TOLOG}
# rmf
if [ -r "${rmf}" ]; then
    RMF=${rmf}
elif [ -r "${DFT_RMF}" ]; then
    RMF=${DFT_RMF}
else
    read -p "> provide the RMF to use: " RMF
    if [ ! -r "${RMF}" ]; then
        printf "ERROR: cannot access given \`${RMF}'\n"
        exit ${ERR_RMF}
    fi
fi
printf "## use RMF: \`${RMF}'\n" | ${TOLOG}
# arf & rmf }}}

# check given dmgroup parameters: grouptype, grouptypeval, binspec
if [ -z "${grouptype}" ]; then
    GROUPTYPE="${DFT_GROUPTYPE}"
elif [ "x${grouptype}" = "xNUM_CTS" ] || [ "x${grouptype}" = "xBIN" ]; then
    GROUPTYPE="${grouptype}"
else
    printf "ERROR: given grouptype \`${grouptype}' invalid.\n"
    exit ${ERR_GRPTYPE}
fi
printf "## use grouptype: \`${GROUPTYPE}'\n" | ${TOLOG}
if [ ! -z "${grouptypeval}" ]; then
    GROUPTYPEVAL="${grouptypeval}"
else
    GROUPTYPEVAL="${DFT_GROUPTYPEVAL}"
fi
printf "## use grouptypeval: \`${GROUPTYPEVAL}'\n" | ${TOLOG}
if [ ! -z "${binspec}" ]; then
    BINSPEC="${binspec}"
else
    BINSPEC="${DFT_BINSPEC}"
fi
printf "## use binspec: \`${BINSPEC}'\n" | ${TOLOG}
## parameters }}}

## prepare parameter files (pfiles) {{{
CIAO_TOOLS="dmstat dmkeypar dmhedit dmextract dmlist dmgroup"

# Copy necessary pfiles for localized usage
for tool in ${CIAO_TOOLS}; do
    pfile=`paccess ${tool}`
    [ -n "${pfile}" ] && punlearn ${tool} && cp -Lvf ${pfile} .
done

# Modify environment variable 'PFILES' to use local pfiles first
export PFILES="./:${PFILES}"
## pfiles }}}

### main ###

## D_A {{{
D_A_CM=`cosmo_calc ${REDSHIFT} | \grep -i 'd_a.*cm' | awk -F'=' '{ print $2 }' | awk '{ print $1 }'`
printf "D_A_CM(${REDSHIFT})= ${D_A_CM}\n"
## D_A }}}

## region related {{{
## generate the needed region file
printf "generate the output region file ...\n"
cat > ${REG_OUT} << _EOF_
# Region file format: CIAO version 1.0
pie(${XC},${YC},${R_IN},${R_OUT},0,360)
_EOF_

## open the evt file to verify or modify
printf "## check the generated pie region ...\n"
printf "## if modified, save with the same name \`${REG_OUT}' (overwrite)\n"
ds9 ${EVT} -regions ${REG_OUT} -cmap he -bin factor 4

## check the (modified) region (pie region end angle)
printf "check the above region (for pie region end angle) ...\n"
INVALID=`\grep -i 'pie' ${REG_OUT} | awk -F'[,()]' '$7 > 360'`
if [ "x${INVALID}" != "x" ]; then
    printf "*** WARNING: there are pie regions' END_ANGLE > 360\n" | ${TOLOG}
    printf "*** will to fix ...\n"
    mv -fv ${REG_OUT} ${REG_OUT}_tmp
    # using `awk' to fix
    awk -F'[,()]' '{
        if ($7 > 360) {
            printf "%s(%.2f,%.2f,%.2f,%.2f,%.2f,%.2f)\n", $1,$2,$3,$4,$5,$6,($7-360)
        }
        else {
            print $0
        }
    }' ${REG_OUT}_tmp > ${REG_OUT}
    rm -f ${REG_OUT}_tmp
fi
## region related }}}

## generate spectrum {{{
# check counts
punlearn dmlist
CNT_RC=`dmlist infile="${EVT}[sky=region(${REG_OUT})][energy=700:7000]" opt=block | \grep 'EVENTS' | awk '{ print $8 }'`
if [ ${CNT_RC} -lt 500 ]; then
    F_WC=true
    WC="LOW_COUNTS"
    printf "*** WARNING: counts_in_0.048R500=${CNT_RC} < 500 ***\n"
fi

# object
printf "extract object spectrum \`${AVGT_SPEC}' ...\n"
AVGT_SPEC="${REG_OUT%.reg}.pi"
AVGT_SPEC_GRP="${AVGT_SPEC%.pi}_grp.pi"
punlearn dmextract
dmextract infile="${EVT}[sky=region(${REG_OUT})][bin PI]" \
    outfile="${AVGT_SPEC}" wmap="[bin det=8]" clobber=yes
# group spectrum
printf "group object spectrum ...\n"
punlearn dmgroup
dmgroup infile="${AVGT_SPEC}" outfile="${AVGT_SPEC_GRP}" \
    grouptype="${GROUPTYPE}" grouptypeval=${GROUPTYPEVAL} \
    binspec="${BINSPEC}" xcolumn="CHANNEL" ycolumn="COUNTS" \
    clobber=yes

# background
printf "generate the background spectrum ...\n"
AVGT_BKG="${AVGT_SPEC%.pi}_bkg.pi"
if [ "${USE_BLANKSKY}" = "YES" ]; then
    # use blanksky as background file
    printf "extract spectrum from blanksky ...\n"
    punlearn dmextract
    dmextract infile="${BLANKSKY}[sky=region(${REG_OUT})][bin PI]" \
        outfile=${AVGT_BKG} wmap="[bin det=8]" clobber=yes
elif [ "${USE_LBKG_REG}" = "YES" ]; then
    printf "extract local background ...\n"
    punlearn dmextract
    dmextract infile="${EVT}[sky=region(${BKGD})][bin PI]" \
        outfile=${AVGT_BKG} wmap="[bin det=8]" clobber=yes
elif [ "${USE_BKG_SPEC}" = "YES" ]; then
    printf "copy specified background spectrum ...\n"
    cp -fv ${BKG_SPEC} ${AVGT_BKG}
fi

printf "renormalize the background ...\n"
bkg_renorm ${AVGT_SPEC} ${AVGT_BKG}
## spectrum }}}

## generate XSPEC script {{{
printf "generate a XSPEC script ...\n"
[ -e "${XSPEC_SCRIPT}" ] && mv -fv ${XSPEC_SCRIPT} ${XSPEC_SCRIPT}_bak
cat > ${XSPEC_SCRIPT} << _EOF_
## XSPEC script
## spectrum analysis to get the average temperatue with (0-0.048 R500)
##
## generated by: \``basename $0`'
## date: \``date`'
##

# xspec settings
statistic chi
abund grsa
query yes

# data
data ${AVGT_SPEC_GRP}
response ${RMF}
arf ${ARF}
backgrnd ${AVGT_BKG}

# fitting range
ignore bad
ignore 0.0-0.7,7.0-**

# plot related
setplot energy

method leven 10 0.01
xsect bcmc
cosmo 70 0 0.73
xset delta 0.01
systematic 0

# model
model wabs*apec
         ${N_H}     -0.001          0          0     100000      1e+06
            1.0       0.01      0.008      0.008         64         64
            0.4      0.001          0          0          5          5
    ${REDSHIFT}      -0.01     -0.999     -0.999         10         10
            0.0       0.01          0          0      1e+24      1e+24

## xspec script end

proc calc_cooling_time {} {
    set rout ${R_OUT}
    set d_a_cm ${D_A_CM}
    fit 1000
    tclout param 4
    set z  [ lindex \$xspec_tclout 0 ]
    tclout param 2
    set T [ lindex \$xspec_tclout 0 ]
    tclout param 5
    set norm [ lindex \$xspec_tclout 0 ]
    newpar 1 0
    dummyrsp .001 100
    flux .001 100

    tclout flux
    set flux [ lindex \$xspec_tclout 0 ]
    puts "flux(0.01-100kev): \$flux"
    set rout_cm [ expr \$rout*.492/3600/180*3.14159*\$d_a_cm ]
    set V [ expr 4./3.*3.14159*\$rout_cm*\$rout_cm*\$rout_cm ]
    set nenh [ expr \$norm*1E14*4*3.14159*\$d_a_cm*\$d_a_cm*(1+\$z)*(1+\$z)*(1+\$z)*(1+\$z)/\$V ]
    set d_l_cm [ expr \$d_a_cm*(1+\$z)*(1+\$z) ]
    set ne_np_ratio 1.2
    set ne [ expr sqrt(\$nenh*\$ne_np_ratio) ]
    set lx [ expr \$flux*4*3.14159*\$d_l_cm*\$d_l_cm ]
    set kb 1.602E-9
    set ct [ expr 3./2.*(\$ne+\$ne/\$ne_np_ratio)*\$kb*\$T*\$V/\$lx ]
    set ct_gyr [ expr \$ct/(3600*24*365.25*1E9) ]
    puts "Cooling_time= \$ct_gyr Gyr"
}

fit 1000
calc_cooling_time

tclexit
_EOF_
## xspec script }}}

## invoke xspec to calc
if [ "x${F_WC}" = "xtrue" ]; then
    printf "\n*** WC: LOW_COUNTS ***\n"
    printf "*** WARNING: counts_in_0.048R500=${CNT_RC} < 500 ***\n"
else
    [ -e "${CT_RES}" ] && mv -f ${CT_RES} ${CT_RES}_bak
    printf "invoking XSPEC to calculate cooling time ...\n"
    xspec - ${XSPEC_SCRIPT} | tee ${CT_RES}

    OBS_ID=`\grep '"Obs.*ID' ${JSON_FILE} | awk -F':' '{ print $2 }' | tr -d ' ,'`
    OBJ_NAME=`\grep '"Source\ Name' ${JSON_FILE} | awk -F':' '{ print $2 }' | sed -e 's/\ *"//' -e 's/"\ *,$//'`
    CT=`\grep -i '^Cooling_time' ${CT_RES} | awk '{ print $2 }'`

    printf "\n" | tee -a ${CT_RES}
    printf "# OBS_ID,OBJ_NAME,CT_gyr\n" | tee -a ${CT_RES}
    printf "# $OBS_ID,$OBJ_NAME,$CT\n" | tee -a ${CT_RES}
fi

exit 0