From e3923265d0d6949407a83726e9a9bd5d97079221 Mon Sep 17 00:00:00 2001 From: Weitian LI Date: Wed, 18 Jun 2014 22:20:59 +0800 Subject: Initial commit Added files: * mass_profile: developed by Junhua GU, modified by Weitian LI * opt_utilities: developed by Junhua GU * tools/cosmo_calc: originated from 'calc_distance', modified * scripts: scripts used to process Chandra ACIS data * files: useful files used in processing * HOWTO_chandra_acis_process.txt * README.md --- mod_ly/fit_dbeta_nfw_mass_profile.sh | 235 +++++++++++++++++++++++++++++++++++ 1 file changed, 235 insertions(+) create mode 100755 mod_ly/fit_dbeta_nfw_mass_profile.sh (limited to 'mod_ly/fit_dbeta_nfw_mass_profile.sh') diff --git a/mod_ly/fit_dbeta_nfw_mass_profile.sh b/mod_ly/fit_dbeta_nfw_mass_profile.sh new file mode 100755 index 0000000..be0404f --- /dev/null +++ b/mod_ly/fit_dbeta_nfw_mass_profile.sh @@ -0,0 +1,235 @@ +#!/bin/sh + +# modified by LIweitaNux, 2012/09/06 +export PATH="/usr/local/bin:/usr/bin:/bin:$PATH" + +if [ $# -eq 1 ]; then + cfg_file="$1" +elif [ $# -eq 2 ]; then + cfg_file="$1" + CENTER_VAL="YES" +else + echo "Usage: $0 [c]" + 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 PGPLOT_FONT="${HEADAS}/lib/grfont.dat" +printf "## PGPLOT_FONT: \`${PGPLOT_FONT}'\n" + +if [ "$0" = `basename $0` ]; then + script_path=`which $0` + base_path=`dirname ${script_path}` +else + base_path=`dirname $0` +fi +printf "## base_path: \`${base_path}'\n" +printf "## use configuration file: \`${cfg_file}'\n" + +#initialize profile type name +t_profile_type=`grep '^t_profile' $cfg_file | awk '{print $2}'` +printf "## t_profile_type: \`$t_profile_type'\n" +#initialize data file name +t_data_file=`grep '^t_data_file' $cfg_file | awk '{print $2}'` +t_param_file=`grep '^t_param_file' $cfg_file | awk '{print $2}'` +#initialize sbp config file +sbp_cfg=`grep '^sbp_cfg' $cfg_file | awk '{print $2}'` +#initialize the temperature profile file +T_file=`grep '^T_file' $sbp_cfg | awk '{print $2}'` +cfunc_file=`grep '^cfunc_file' ${sbp_cfg} |awk '{print $2}'` +abund=`grep '^abund' ${cfg_file} |awk '{print $2}'` +nh=`grep '^nh' ${cfg_file} |awk '{print $2}'` +## calc `cm_per_pixel' instead {{{ +# cm_per_pixel=`grep '^cm_per_pixel' $sbp_cfg|awk '{print $2}'` +z=`grep '^z' $sbp_cfg | awk '{ print $2 }'` +cm_per_pixel=`${base_path}/calc_distance ${z} | grep 'cm_per_pixel' | awk '{ print $2 }'` +sed -i'' "s/^cm_per_pixel.*$/cm_per_pixel ${cm_per_pixel}/" ${sbp_cfg} +printf "## redshift: ${z}, cm_per_pixel: ${cm_per_pixel}\n" +## cm_per_pixel }}} +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" +## sbp {{{ +sbp_data_file=`grep '^sbp_file' $sbp_cfg | awk '{print $2}'` +radius_sbp_file=`grep '^radius_sbp_file' ${cfg_file} | awk '{print $2}'` + +if [ "x$radius_sbp_file" = "x" ]; then + echo "ERROR, must have radius_sbp_file assigned, this file should be a 4-column file, which contains the radius, radius err, sbp, and sbp err" + exit 200 +fi + +TMP_RSBP="_tmp_rsbp.txt" +[ -e "${TMP_RSBP}" ] && rm -f ${TMP_RSBP} +cat ${radius_sbp_file} | sed 's/#.*$//' | grep -Ev '^\s*$' > ${TMP_RSBP} +# mv -f _tmp_rsbp.txt ${radius_sbp_file} +radius_sbp_file="${TMP_RSBP}" +## sbp }}} + +# determine which temperature profile to be used, and fit the T profile {{{ +if [ "$t_profile_type" = "zyy" ]; then + $base_path/fit_zyy_model $t_data_file $t_param_file $cm_per_pixel + # mv -f zyy_dump.qdp ${T_file} + mv -f zyy_dump.qdp zyy_dump_center.qdp +elif [ "$t_profile_type" = "m0603246" ]; then + $base_path/fit_m0603246 $t_data_file $cm_per_pixel + # mv -f m0603246_dump.qdp ${T_file} + mv -f m0603246_dump.qdp m0603246_dump_center.qdp +elif [ "$t_profile_type" = "wang2012" ]; then + $base_path/fit_wang2012_model $t_data_file $t_param_file $cm_per_pixel 2> /dev/null | tee wang2012_center_param.txt + # mv -f wang2012_dump.qdp ${T_file} + cp -fv wang2012_dump.qdp ${T_file} + T_file_center="wang2012_dump_center.qdp" + mv -fv wang2012_dump.qdp ${T_file_center} + mv -fv fit_result.qdp wang2012_fit_center.qdp +elif [ "$t_profile_type" = "allen" ]; then + $base_path/fit_allen_model $t_data_file $cm_per_pixel + # mv -f allen_dump.qdp ${T_file} + mv -f allen_dump.qdp allen_dump_center.qdp +elif [ "$t_profile_type" = "zzl" ]; then + $base_path/fit_zzl_model $t_data_file $t_param_file + # mv -f zzl_dump.qdp ${T_file} + mv -f zzl_dump.qdp zzl_dump_center.qdp +else + printf "ERROR: temperature profile name \`${t_profile_type}' invalid!\n" + exit 10 +fi +# temp profile }}} + +$base_path/coolfunc_calc_bolo.sh ${T_file_center} $abund $nh $z cfunc_bolo.dat +$base_path/coolfunc_calc.sh ${T_file_center} $abund $nh $z $cfunc_file +mv -fv flux_cnt_ratio.txt flux_cnt_ratio_center.txt +# fit sbp (double-beta) +$base_path/fit_dbeta_sbp $sbp_cfg 2> /dev/null +mv -fv dbeta_param.txt dbeta_param_center.txt +cat dbeta_param_center.txt +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 +$base_path/fit_nfw_mass mass_int.dat $z 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 +$base_path/cooling_time rho_fit_center.dat ${T_file_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}\n" + +## center value {{{ +if [ "${CENTER_VAL}" = "YES" ]; then + $base_path/analyze_mass_profile.py 200 c + $base_path/analyze_mass_profile.py 500 c + $base_path/extract_tcool.py $rcool | tee cooling_time_result.txt + exit 0 +fi +## center value }}} + +# clean previous files +rm -f summary_shuffle_mass_profile.qdp +rm -f summary_overdensity.qdp +rm -f summary_mass_profile.qdp +rm -f summary_gas_mass_profile.qdp + +## count +COUNT=1 + +#100 times of Monte-carlo simulation to determine error +#just repeat above steps +printf "\n+++++++++++++++++++ Monte Carlo +++++++++++++++++++++\n" +for i in `seq 1 100`; do + # echo $t_data_file + $base_path/shuffle_T.py $t_data_file temp_shuffled_t.dat + $base_path/shuffle_sbp.py $sbp_data_file temp_shuffled_sbp.dat + # t_data_file=temp_shuffled_t.dat +#exit + + if [ "$t_profile_type" = "zyy" ]; then + $base_path/fit_zyy_model temp_shuffled_t.dat $t_param_file $cm_per_pixel + mv -f zyy_dump.qdp ${T_file} + elif [ "$t_profile_type" = "m0603246" ]; then + $base_path/fit_m0603246 temp_shuffled_t.dat $cm_per_pixel + mv -f m0603246_dump.qdp ${T_file} + elif [ "$t_profile_type" = "wang2012" ]; then + $base_path/fit_wang2012_model temp_shuffled_t.dat $t_param_file $cm_per_pixel 2> /dev/null + mv -f wang2012_dump.qdp ${T_file} + elif [ "$t_profile_type" = "allen" ]; then + $base_path/fit_allen_model temp_shuffled_t.dat $cm_per_pixel + mv -f allen_dump.qdp ${T_file} + elif [ "$t_profile_type" = "zzl" ]; then + $base_path/fit_zzl_model temp_shuffled_t.dat $t_param_file + mv -f zzl_dump.qdp ${T_file} + else + printf "ERROR: temperature profile name \`${t_profile_type}' invalid!\n" + exit 10 + fi + +#exit + : > temp_sbp.cfg + + cat $sbp_cfg | while read l; do + if echo $l | grep -q 'sbp_file'; then + echo "sbp_file temp_shuffled_sbp.dat" >> temp_sbp.cfg + elif echo $l | grep -q 'T_file'; then + echo "T_file ${T_file}" >> temp_sbp.cfg + else + echo $l >> temp_sbp.cfg + fi + done + + ## count + printf "## ${COUNT} ##\n" + COUNT=`expr ${COUNT} + 1` + printf "## `pwd` ##\n" + $base_path/coolfunc_calc.sh ${T_file} ${abund} ${nh} ${z} ${cfunc_file} + + # double-beta model + $base_path/fit_dbeta_sbp temp_sbp.cfg 2> /dev/null + [ -r "dbeta_param.txt" ] && cat dbeta_param.txt + $base_path/fit_nfw_mass mass_int.dat $z 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 + +done # end `while' +printf "\n+++++++++++++++++++ Monte Carlo +++++++++++++++++++++\n" + +#analys the errors +printf "\n+++++++++++++++ RESULTS (double-beta) +++++++++++++++\n" +RESULT="results_dbeta_mrl.txt" +[ -e "${RESULT}" ] && mv -fv ${RESULT} ${RESULT}_bak + +$base_path/analyze_mass_profile.py 200 | tee -a ${RESULT} +$base_path/analyze_mass_profile.py 500 | tee -a ${RESULT} +$base_path/analyze_mass_profile.py 1500 | tee -a ${RESULT} +$base_path/analyze_mass_profile.py 2500 | tee -a ${RESULT} + +R200_VAL=`grep '^r200' ${RESULT} | awk '{ print $2 }'` +R500_VAL=`grep '^r500' ${RESULT} | awk '{ print $2 }'` +printf "\n## R200: ${R200_VAL}, R500: ${R500_VAL}\n" +L200=`$base_path/calc_lx $radius_sbp_file flux_cnt_ratio_center.txt $z ${R200_VAL} $t_data_file | grep '^Lx' | awk '{ print $2,$3,$4 }'` +L500=`$base_path/calc_lx $radius_sbp_file flux_cnt_ratio_center.txt $z ${R500_VAL} $t_data_file | grep '^Lx' | awk '{ print $2,$3,$4 }'` +printf "L200= ${L200} erg/s\n" | tee -a ${RESULT} +printf "L500= ${L500} erg/s\n" | tee -a ${RESULT} +$base_path/extract_tcool.py $rcool | tee cooling_time_result.txt +printf "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++\n" + -- cgit v1.2.2