Rename fit_nfwmass.sh to fit_mass; Update to the new style config

1 files changed, 289 insertions, 0 deletions

diff --git a/mass_profile/fit_mass.sh b/mass_profile/fit_mass.sh
new file mode 100755
index 0000000..2854b5a
--- /dev/null
+++ b/mass_profile/fit_mass.sh
@@ -0,0 +1,289 @@
+#!/bin/sh
+#
+# Front-end script used to calculate the mass and related values.
+#
+# Output:
+#   * final_result.txt / center_only_results.txt
+#   * beta_param.txt / dbeta_param_center.txt
+#   * gas_mass_int_center.qdp
+#   * mass_int_center.qdp
+#   * nfw_fit_center.qdp
+#   * nfw_param_center.txt
+#   * overdensity_center.qdp
+#   * rho_fit_center.dat
+#   * rho_fit_center.qdp
+#   * sbp_fit_center.qdp
+#   * entropy_center.qdp
+#   * wang2012_param_center.txt
+#   * tprofile_dump_center.qdp
+#   * tprofile_fit_center.qdp
+#   * summary_mass_profile.qdp
+#   * summary_overdensity.qdp
+#   * summary_gas_mass_profile.qdp
+#   * summary_entropy.qdp
+#
+# Junhua Gu
+# Weitian LI
+# 2016-06-07
+#
+# Change logs:
+# 2017-02-07, Weitian LI
+#   * Rename from `fit_nfwmass.sh` to `fit_mass.sh`
+#   * Update to use the new style configuration files
+#
+
+if [ $# -eq 1 ] || [ $# -eq 2 ]; then
+    :
+else
+    echo "usage:"
+    echo "    `basename $0` <mass.conf> [c]"
+    echo ""
+    echo "arguments:"
+    echo "    <mass.conf>: config file for mass profile calculation"
+    echo "    [c]: optional; if specified, do not calculate the errors"
+    exit 1
+fi
+
+if ! which xspec > /dev/null; then
+    printf "*** ERROR: please initialize HEASOFT first\n"
+    exit 2
+fi
+
+if [ -z "${HEADAS}" ]; then
+    printf "*** ERROR: variable \`HEADAS' not properly set\n"
+    exit 3
+fi
+
+export PATH="/usr/local/bin:/usr/bin:/bin:$PATH"
+
+base_path=$(dirname $(realpath $0))
+printf "## base_path: \`${base_path}'\n"
+mass_cfg="$1"
+printf "## use configuration file: \`${mass_cfg}'\n"
+case "$2" in
+    [cC])
+        F_C="YES"
+        ;;
+    *)
+        F_C="NO"
+        ;;
+esac
+
+nh=`grep            '^nh'            ${mass_cfg} | awk '{ print $2 }'`
+abund=`grep         '^abund'         ${mass_cfg} | awk '{ print $2 }'`
+nfw_rmin_kpc=`grep  '^nfw_rmin_kpc'  ${mass_cfg} | awk '{ print $2 }'`
+tprofile_data=`grep '^tprofile_data' ${mass_cfg} | awk '{ print $2 }'`
+tprofile_cfg=`grep  '^tprofile_cfg'  ${mass_cfg} | awk '{ print $2 }'`
+sbp_cfg=`grep       '^sbp_cfg'       ${mass_cfg} | awk '{ print $2 }'`
+
+# sbp config file
+sbp_data=`grep      '^sbp_data'      ${sbp_cfg} | awk '{ print $2 }'`
+tprofile=`grep      '^tprofile'      ${sbp_cfg} | awk '{ print $2 }'`
+cfunc_profile=`grep '^cfunc_profile' ${sbp_cfg} | awk '{ print $2 }'`
+z=`grep             '^z'             ${sbp_cfg} | awk '{ print $2 }'`
+cm_per_pixel=`cosmo_calc ${z} | grep 'cm/pixel' | awk -F':' '{ print $2 }'`
+sed -i'' "s/^cm_per_pixel.*$/cm_per_pixel   ${cm_per_pixel}/" ${sbp_cfg}
+
+da=`python -c "print($cm_per_pixel/(.492/3600/180*3.1415926))"`
+dl=`python -c "print($da*(1+$z)**2)"`
+printf "da= ${da}\n"
+printf "dl= ${dl}\n"
+
+if grep -q '^beta2' $sbp_cfg; then
+    MODEL="dbeta"
+    MODEL_NAME="double-beta"
+else
+    MODEL="beta"
+    MODEL_NAME="single-beta"
+fi
+
+PROG_TPROFILE="fit_wang2012_model"
+tprofile_dump="wang2012_dump.qdp"
+tprofile_param_center="wang2012_param_center.txt"
+tprofile_fit_center="tprofile_fit_center.qdp"
+tprofile_center="tprofile_dump_center.qdp"
+
+${base_path}/${PROG_TPROFILE} ${tprofile_data} ${tprofile_cfg} \
+            ${cm_per_pixel} 2> /dev/null | tee ${tprofile_param_center}
+cp -fv ${tprofile_dump} ${tprofile}
+mv -fv ${tprofile_dump} ${tprofile_center}
+mv -fv fit_result.qdp ${tprofile_fit_center}
+
+$base_path/coolfunc_calc.sh ${tprofile_center} \
+                            ${abund} ${nh} ${z} ${cfunc_profile} cfunc_bolo.dat
+cfunc_profile_center="coolfunc_profile_center.txt"
+cp -f ${cfunc_profile} ${cfunc_profile_center}
+mv -fv flux_cnt_ratio.txt flux_cnt_ratio_center.txt
+
+PROG_SBPFIT="fit_${MODEL}_sbp"
+RES_SBPFIT="${MODEL}_param.txt"
+RES_SBPFIT_CENTER="${MODEL}_param_center.txt"
+${base_path}/${PROG_SBPFIT} ${sbp_cfg} 2> /dev/null
+mv -fv ${RES_SBPFIT} ${RES_SBPFIT_CENTER}
+cat ${RES_SBPFIT_CENTER}
+mv -fv sbp_fit.qdp sbp_fit_center.qdp
+mv -fv rho_fit.qdp rho_fit_center.qdp
+mv -fv rho_fit.dat rho_fit_center.dat
+mv -fv entropy.qdp entropy_center.qdp
+${base_path}/fit_nfw_mass mass_int.dat ${z} ${nfw_rmin_kpc} 2> /dev/null
+mv -fv nfw_param.txt      nfw_param_center.txt
+mv -fv nfw_fit_result.qdp nfw_fit_center.qdp
+mv -fv nfw_dump.qdp       mass_int_center.qdp
+mv -fv overdensity.qdp    overdensity_center.qdp
+mv -fv gas_mass_int.qdp   gas_mass_int_center.qdp
+
+#exit 233
+
+## cooling time (-> use 'ciao_calc_ct.sh')
+${base_path}/cooling_time rho_fit_center.dat \
+            ${tprofile_center} cfunc_bolo.dat \
+            ${dl} ${cm_per_pixel} > cooling_time.dat
+## radius to calculate tcool, not the cooling time!
+rcool=`${base_path}/analyze_mass_profile.py 500 c | grep '^r500' | awk -F'=' '{ print .048*$2 }'`
+printf "rcool= ${rcool} (kpc)\n"
+
+## only calculate central value {{{
+if [ "${F_C}" = "YES" ]; then
+    RES_CENTER="center_only_results.txt"
+    [ -e "${RES_CENTER}" ] && mv -f ${RES_CENTER} ${RES_CENTER}_bak
+    ${base_path}/analyze_mass_profile.py  200 c | tee -a ${RES_CENTER}
+    ${base_path}/analyze_mass_profile.py  500 c | tee -a ${RES_CENTER}
+    ${base_path}/analyze_mass_profile.py 1500 c | tee -a ${RES_CENTER}
+    ${base_path}/analyze_mass_profile.py 2500 c | tee -a ${RES_CENTER}
+    ${base_path}/fg_2500_500.py c               | tee -a ${RES_CENTER}
+    ${base_path}/extract_tcool.py cooling_time.dat ${rcool} | tee -a ${RES_CENTER}
+    exit 0
+fi
+## central value }}}
+
+## ------------------------------------------------------------------
+
+# clean previous files
+rm -f summary_overdensity.qdp
+rm -f summary_mass_profile.qdp
+rm -f summary_gas_mass_profile.qdp
+rm -f summary_entropy.qdp
+
+# Estimate the errors of Lx and Fx by Monte Carlo simulation
+printf "\n+++++++++++++++++++ Monte Carlo +++++++++++++++++++++\n"
+MC_TIMES=100
+for i in `seq 1 ${MC_TIMES}`; do
+    ${base_path}/shuffle_profile.py ${tprofile_data} tmp_tprofile.txt
+    ${base_path}/shuffle_profile.py ${sbp_data} tmp_sbprofile.txt
+
+    # temperature profile
+    ${base_path}/${PROG_TPROFILE} tmp_tprofile.txt ${tprofile_cfg} \
+                ${cm_per_pixel} 2> /dev/null
+    mv -f ${tprofile_dump} ${tprofile}
+
+    TMP_SBP_CFG="tmp_sbp.cfg"
+    [ -e "${TMP_SBP_CFG}" ] && rm -f ${TMP_SBP_CFG}
+    cat ${sbp_cfg} | while read l; do
+        if echo "${l}" | grep -q '^sbp_data' >/dev/null; then
+            echo "sbp_data  tmp_sbprofile.txt" >> ${TMP_SBP_CFG}
+        elif echo "${l}" | grep -q '^tprofile' >/dev/null; then
+            echo "tprofile  ${tprofile}" >> ${TMP_SBP_CFG}
+        else
+            echo "${l}" >> ${TMP_SBP_CFG}
+        fi
+    done
+
+    printf "## ${i} / ${MC_TIMES} ##\n"
+    printf "## `pwd -P` ##\n"
+    ${base_path}/coolfunc_calc.sh ${tprofile} ${abund} ${nh} ${z} ${cfunc_profile}
+    ${base_path}/${PROG_SBPFIT} ${TMP_SBP_CFG} 2> /dev/null
+    cat ${RES_SBPFIT}
+    ${base_path}/fit_nfw_mass mass_int.dat ${z} ${nfw_rmin_kpc} 2> /dev/null
+    cat nfw_dump.qdp     >> summary_mass_profile.qdp
+    echo "no no no"      >> summary_mass_profile.qdp
+    cat overdensity.qdp  >> summary_overdensity.qdp
+    echo "no no no"      >> summary_overdensity.qdp
+    cat gas_mass_int.qdp >> summary_gas_mass_profile.qdp
+    echo "no no no"      >> summary_gas_mass_profile.qdp
+    cat entropy.qdp      >> summary_entropy.qdp
+    echo "no no no"      >> summary_entropy.qdp
+done  # end `while'
+
+# recover the files of original center values
+cp -f ${cfunc_profile_center} ${cfunc_profile}
+cp -f ${tprofile_center} ${tprofile}
+printf "\n+++++++++++++++++ MONTE CARLO END +++++++++++++++++++\n"
+
+## analyze results
+RES_TMP="_tmp_result_mrl.txt"
+RES_FINAL="final_result.txt"
+[ -e "${RES_TMP}" ]   && mv -fv ${RES_TMP}   ${RES_TMP}_bak
+[ -e "${RES_FINAL}" ] && mv -fv ${RES_FINAL} ${RES_FINAL}_bak
+
+${base_path}/analyze_mass_profile.py  200 | tee -a ${RES_TMP}
+${base_path}/analyze_mass_profile.py  500 | tee -a ${RES_TMP}
+${base_path}/analyze_mass_profile.py 1500 | tee -a ${RES_TMP}
+${base_path}/analyze_mass_profile.py 2500 | tee -a ${RES_TMP}
+
+R200_VAL=`grep  '^r200'  ${RES_TMP} | awk '{ print $2 }'`
+R500_VAL=`grep  '^r500'  ${RES_TMP} | awk '{ print $2 }'`
+R1500_VAL=`grep '^r1500' ${RES_TMP} | awk '{ print $2 }'`
+R2500_VAL=`grep '^r2500' ${RES_TMP} | awk '{ print $2 }'`
+printf "## R200:  ${R200_VAL}\n"
+printf "## R500:  ${R500_VAL}\n"
+printf "## R1500: ${R1500_VAL}\n"
+printf "## R2500: ${R2500_VAL}\n"
+L200E=`${base_path}/calc_lx  ${sbprofile} flux_cnt_ratio_center.txt $z $R200_VAL  ${tprofile_data} | grep '^Lx' | awk '{ print $2,$3,$4 }'`
+L500E=`${base_path}/calc_lx  ${sbprofile} flux_cnt_ratio_center.txt $z $R500_VAL  ${tprofile_data} | grep '^Lx' | awk '{ print $2,$3,$4 }'`
+L1500E=`${base_path}/calc_lx ${sbprofile} flux_cnt_ratio_center.txt $z $R1500_VAL ${tprofile_data} | grep '^Lx' | awk '{ print $2,$3,$4 }'`
+L2500E=`${base_path}/calc_lx ${sbprofile} flux_cnt_ratio_center.txt $z $R2500_VAL ${tprofile_data} | grep '^Lx' | awk '{ print $2,$3,$4 }'`
+R200E=`grep   '^r200'             ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+R500E=`grep   '^r500'             ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+R1500E=`grep  '^r1500'            ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+R2500E=`grep  '^r2500'            ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+M200E=`grep   '^m200'             ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+M500E=`grep   '^m500'             ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+M1500E=`grep  '^m1500'            ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+M2500E=`grep  '^m2500'            ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+MG200E=`grep  '^gas_m200'         ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+MG500E=`grep  '^gas_m500'         ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+MG1500E=`grep '^gas_m1500'        ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+MG2500E=`grep '^gas_m2500'        ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+FG200E=`grep  '^gas_fraction200'  ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+FG500E=`grep  '^gas_fraction500'  ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+FG1500E=`grep '^gas_fraction1500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+FG2500E=`grep '^gas_fraction2500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
+
+printf "\n+++++++++++++++ RESULTS (${MODEL_NAME}) +++++++++++++++\n"
+printf "model: ${MODEL_NAME}\n"                | tee -a ${RES_FINAL}
+printf "\n"                                    | tee -a ${RES_FINAL}
+cat ${RES_SBPFIT_CENTER}                       | tee -a ${RES_FINAL}
+printf "\n"                                    | tee -a ${RES_FINAL}
+printf "r200= ${R200E} kpc\n"                  | tee -a ${RES_FINAL}
+printf "m200= ${M200E} M_sun\n"                | tee -a ${RES_FINAL}
+printf "L200= ${L200E} erg/s\n"                | tee -a ${RES_FINAL}
+printf "gas_m200= ${MG200E} M_sun\n"           | tee -a ${RES_FINAL}
+printf "gas_fraction200= ${FG200E} x100%%\n"   | tee -a ${RES_FINAL}
+printf "r500= ${R500E} kpc\n"                  | tee -a ${RES_FINAL}
+printf "m500= ${M500E} M_sun\n"                | tee -a ${RES_FINAL}
+printf "L500= ${L500E} erg/s\n"                | tee -a ${RES_FINAL}
+printf "gas_m500= ${MG500E} M_sun\n"           | tee -a ${RES_FINAL}
+printf "gas_fraction500= ${FG500E} x100%%\n"   | tee -a ${RES_FINAL}
+printf "r1500= ${R1500E} kpc\n"                | tee -a ${RES_FINAL}
+printf "m1500= ${M1500E} M_sun\n"              | tee -a ${RES_FINAL}
+printf "L1500= ${L1500E} erg/s\n"              | tee -a ${RES_FINAL}
+printf "gas_m1500= ${MG1500E} M_sun\n"         | tee -a ${RES_FINAL}
+printf "gas_fraction1500= ${FG1500E} x100%%\n" | tee -a ${RES_FINAL}
+printf "r2500= ${R2500E} kpc\n"                | tee -a ${RES_FINAL}
+printf "m2500= ${M2500E} M_sun\n"              | tee -a ${RES_FINAL}
+printf "L2500= ${L2500E} erg/s\n"              | tee -a ${RES_FINAL}
+printf "gas_m2500= ${MG2500E} M_sun\n"         | tee -a ${RES_FINAL}
+printf "gas_fraction2500= ${FG2500E} x100%%\n" | tee -a ${RES_FINAL}
+printf "\n"                                    | tee -a ${RES_FINAL}
+printf "gas mass 200= ${MG200E} M_sun\n"       | tee -a ${RES_FINAL}
+printf "gas fractho 200= ${FG200E} x100%%\n"   | tee -a ${RES_FINAL}
+printf "gas mass 500= ${MG500E} M_sun\n"       | tee -a ${RES_FINAL}
+printf "gas fractho 500= ${FG500E} x100%%\n"   | tee -a ${RES_FINAL}
+printf "gas mass 1500= ${MG1500E} M_sun\n"     | tee -a ${RES_FINAL}
+printf "gas fractho 1500= ${FG1500E} x100%%\n" | tee -a ${RES_FINAL}
+printf "gas mass 2500= ${MG2500E} M_sun\n"     | tee -a ${RES_FINAL}
+printf "gas fractho 2500= ${FG2500E} x100%%\n" | tee -a ${RES_FINAL}
+printf "\n"                                    | tee -a ${RES_FINAL}
+${base_path}/fg_2500_500.py                    | tee -a ${RES_FINAL}
+${base_path}/extract_tcool.py cooling_time.dat ${rcool} | tee -a ${RES_FINAL}
+printf "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++\n"