diff options
| -rwxr-xr-x | mass_profile/analyze_entropy_profile.py | 54 | ||||
| -rwxr-xr-x | mass_profile/calc_all_entropy.sh | 39 | ||||
| -rwxr-xr-x | mass_profile/calc_entropy.py | 104 | ||||
| -rwxr-xr-x | mass_profile/fit_beta_entropy.sh | 150 | ||||
| -rwxr-xr-x | mass_profile/fit_dbeta_entropy.sh | 149 |
5 files changed, 104 insertions, 392 deletions
diff --git a/mass_profile/analyze_entropy_profile.py b/mass_profile/analyze_entropy_profile.py deleted file mode 100755 index c123681..0000000 --- a/mass_profile/analyze_entropy_profile.py +++ /dev/null @@ -1,54 +0,0 @@ -#!/usr/bin/env python - -import sys -import numpy -import scipy.interpolate - -center_entropy_file=open('entropy_center.qdp') -entropy_file=open('summary_entropy.qdp') -confidence_level=.68 -rout=float(sys.argv[1]) - -center_s=0 -for i in center_entropy_file: - r,s=i.split() - r=float(r) - s=float(s) - if r>rout: - center_s=s - break - -new_data=True - - -s_list=[] -for i in entropy_file: - if i[0]=='n': - new_data=True - continue - if new_data==False: - continue - r,s=i.split() - r=float(r) - s=float(s) - if r>rout: - new_data=False - s_list.append(s) - -s_idx=-1 - -s_list.sort() -for i in range(len(s_list)-1): - if (center_s-s_list[i])*(center_s-s_list[i+1])<=0: - m_idx=i - break - - -slidx=int(s_idx*(1-confidence_level)) -suidx=s_idx-1+int((len(s_list)-s_idx)*confidence_level) - - -serr1=s_list[slidx]-center_s -serr2=s_list[suidx]-center_s - -print("S=\t%e\t %e/+%e keV cm^2 (1 sigma)"%(center_s,serr1,serr2)) diff --git a/mass_profile/calc_all_entropy.sh b/mass_profile/calc_all_entropy.sh deleted file mode 100755 index 5506ab8..0000000 --- a/mass_profile/calc_all_entropy.sh +++ /dev/null @@ -1,39 +0,0 @@ -#!/usr/bin/env bash - -if [ $# -eq 1 ] -then - : -else - echo "Usage: <global.cfg list file>" - exit -fi - -bdir=`pwd` - -for i in `cat $1` -do - dname=`dirname $i` - #dname=`dirname $dname` - echo $dname - cd $dname ||continue -# mkdir -p profile_entropy -# cp -rf profile/* profile_entropy/ -# cd profile/ - - [ -e result ] && r200file=result - [ -e result ] || r200file=result_checked - - r200=`grep r200 $r200file|awk -F = '{print $2}' |awk -F '-' '{print $1}'` - rout=`python -c "print($r200*.1)"` - echo $rout - if grep n01 source.cfg >/dev/null - then - echo "dbeta" - fit_dbeta_entropy.sh global.cfg $rout - else - echo "beta" - fit_beta_entropy.sh global.cfg $rout - fi - - cd $bdir -done diff --git a/mass_profile/calc_entropy.py b/mass_profile/calc_entropy.py new file mode 100755 index 0000000..736fe7e --- /dev/null +++ b/mass_profile/calc_entropy.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python3 +# +# Calculate the entropy within the specified radius. +# +# Junhua GU +# Weitian LI +# 2016-06-07 +# + +import argparse +import re +from itertools import groupby +import numpy as np + + +def isplit(iterable, splitters): + """ + Credit: https://stackoverflow.com/a/4322780/4856091 + """ + return [list(g) for k, g in groupby(iterable, + lambda x:x in splitters) if not k] + + +def get_entropy(data, r): + """ + Get the entropy *at* the specified radius. + + XXX: whether to interpolate first? + """ + radius = data[:, 0] + entropy = data[:, 1] + s = np.min(entropy[radius > r]) + return s + + +def read_merged_qdp(infile): + """ + Read merged QDP with multiple group of data separated by "no no no". + """ + lines = map(lambda line: re.sub(r"^\s*no\s+no\s+no.*$", "X", + line.strip(), flags=re.I), + open(infile).readlines()) + lines = isplit(lines, ("X",)) + data_groups = [] + for block in lines: + data = [list(map(float, l.split())) for l in block] + data.append(np.row_stack(data)) + return data_groups + + +def calc_error(center_value, mc_values, ci=0.683): + """ + Calculate the uncertainties/errors. + """ + data = np.concatenate([[center_value], mc_values]) + median, q_lower, q_upper = np.percentile(data, q=(50, 50-50*ci, 50+50*ci)) + mean = np.mean(data) + std = np.std(data) + return { + "mean": mean, + "std": std, + "median": median, + "q_lower": q_lower, + "q_upper": q_upper, + } + + +def main(): + parser = argparse.ArgumentParser( + description="Calculate the entropy within the given radius") + parser.add_argument("-C", "--confidence-level", dest="ci", + type=float, default=0.683, + help="confidence level to estimate the errors") + parser.add_argument("center_data", + help="calculate central entropy profile " + + "(e.g., entropy_center.qdp)") + parser.add_argument("mc_data", + help="Merged QDP file of all the Monte Carlo " + + "simulated entropy profiles " + + "(e.g., summary_entropy.qdp)") + parser.add_argument("rout", type=float, help="outer radius (kpc)") + args = parser.parse_args() + + center_data = np.loadtxt(args.center_data) + center_s = get_entropy(center_data, r=args.rout) + + data_groups = read_merged_qdp(args.mc_data) + entropy_list = [] + for dg in data_groups: + s = get_entropy(dg, r=args.rout) + entropy_list.append(s) + results = calc_error(center_s, entropy_list, ci=args.ci) + s_err_lower = results["q_lower"] - center_s + s_err_upper = results["q_upper"] - center_s + + print("entropy= %e %+e/%+e keV cm^2 (ci=%.1f%%)" % + (center_s, s_err_lower, s_err_upper, args.ci * 100)) + print("entropy(mean)= %e" % results["mean"]) + print("entropy(median)= %e" % results["median"]) + print("entropy(std)= %e" % results["std"]) + + +if __name__ == "__main__": + main() diff --git a/mass_profile/fit_beta_entropy.sh b/mass_profile/fit_beta_entropy.sh deleted file mode 100755 index f52e4a1..0000000 --- a/mass_profile/fit_beta_entropy.sh +++ /dev/null @@ -1,150 +0,0 @@ -#!/bin/bash - -echo $# -if [ $# -eq 2 ] -then - : -else - echo "Usage:$0 <cfg file> <radius in kpc>" - exit -fi -export PGPLOT_FONT=`locate grfont.dat|head -1` - -mkdir -p entropy_files - -cfg_file=$1 -base_path=`dirname $0` -echo $base_path -#initialize profile type name -t_profile_type=`grep t_profile $cfg_file|awk '{print $2}'` -#initialize data file name -t_data_file=`grep t_data_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}'` -#echo $t_profile_type -cm_per_pixel=`grep '^cm_per_pixel' $sbp_cfg|awk '{print $2}'` - -#determine which temperature profile to be used, and fit the T profile -if [ $t_profile_type == zyy ] -then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_zyy_model $t_data_file $t_param_file $cm_per_pixel - mv -f zyy_dump.qdp ${T_file} -elif [ $t_profile_type == m0603246 ] -then - $base_path/fit_m0603246 $t_data_file $cm_per_pixel - mv -f m0603246_dump.qdp ${T_file} -elif [ $t_profile_type == wang2012 ] -then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_wang2012_model $t_data_file $t_param_file $cm_per_pixel - mv -f wang2012_dump.qdp ${T_file} -elif [ $t_profile_type == allen ] -then - $base_path/fit_allen_model $t_data_file $cm_per_pixel - mv -f allen_dump.qdp ${T_file} -elif [ $t_profile_type == zzl ] -then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_zzl_model $t_data_file $t_param_file - mv -f zzl_dump.qdp ${T_file} -else - echo temperature profile name invalid! - exit -fi - -cfunc_file=`grep '^cfunc_file' ${sbp_cfg} |awk '{print $2}'` -z=`grep '^z' ${sbp_cfg}|awk '{print $2}'` -abund=`grep '^abund' ${cfg_file} |awk '{print $2}'` -nh=`grep '^nh' ${cfg_file} |awk '{print $2}'` -$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file -mv flux_cnt_ratio.txt flux_cnt_ratio_center_entropy.txt -#fit sbp -$base_path/fit_beta_sbp $sbp_cfg -echo $cfunc_file -#exit - -#store central valu -mv sbp_fit.qdp sbp_fit_center_entropy.qdp -mv mass_int.qdp mass_int_center_entropy.qdp -mv overdensity.qdp overdensity_center_entropy.qdp -mv gas_mass_int.qdp gas_mass_int_center_entropy.qdp -mv entropy.qdp entropy_center.qdp -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 -fi - -cat ${radius_sbp_file} | sed 's/#.*$//' | grep -Ev '^\s*$' > .tmp.txt -mv .tmp.txt ${radius_sbp_file} - -rm -f summary_entropy.qdp -#100 times of Monte-carlo simulation to determine error -#just repeat above steps -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 - -#exit - if [ $t_profile_type == zyy ] - then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $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 - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_wang2012_model temp_shuffled_t.dat $t_param_file $cm_per_pixel - 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 - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_zzl_model temp_shuffled_t.dat $t_param_file - mv -f zzl_dump.qdp ${T_file} - else - echo temperature profile name invalid! - exit - fi - -#exit - echo >temp_sbp.cfg - - cat $sbp_cfg|while read l -do - if echo $l|grep sbp_file >/dev/null - then - echo sbp_file temp_shuffled_sbp.dat >>temp_sbp.cfg - elif echo $l|grep T_file >/dev/null - then - echo T_file ${T_file} >>temp_sbp.cfg - else - echo $l >>temp_sbp.cfg - fi - -done - -#$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file - -$base_path/fit_beta_sbp temp_sbp.cfg - -cat entropy.qdp >>summary_entropy.qdp -echo no no no >>summary_entropy.qdp -done - -$base_path/analyze_entropy_profile.py $2 |tee entropy_result.txt diff --git a/mass_profile/fit_dbeta_entropy.sh b/mass_profile/fit_dbeta_entropy.sh deleted file mode 100755 index 3f1bbda..0000000 --- a/mass_profile/fit_dbeta_entropy.sh +++ /dev/null @@ -1,149 +0,0 @@ -#!/bin/bash - -echo $# -if [ $# -eq 2 ] -then - : -else - echo "Usage:$0 <cfg file> <radius in kpc>" - exit -fi -export PGPLOT_FONT=`locate grfont.dat|head -1` - -cfg_file=$1 -base_path=`dirname $0` -echo $base_path -#initialize profile type name -t_profile_type=`grep t_profile $cfg_file|awk '{print $2}'` -#initialize data file name -t_data_file=`grep t_data_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}'` -#echo $t_profile_type -cm_per_pixel=`grep '^cm_per_pixel' $sbp_cfg|awk '{print $2}'` -#determine which temperature profile to be used, and fit the T profile -if [ $t_profile_type == zyy ] -then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_zyy_model $t_data_file $t_param_file $cm_per_pixel - mv -f zyy_dump.qdp ${T_file} -elif [ $t_profile_type == m0603246 ] -then - $base_path/fit_m0603246 $t_data_file $cm_per_pixel - mv -f m0603246_dump.qdp ${T_file} -elif [ $t_profile_type == wang2012 ] -then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_wang2012_model $t_data_file $t_param_file $cm_per_pixel - mv -f wang2012_dump.qdp ${T_file} -elif [ $t_profile_type == allen ] -then - $base_path/fit_allen_model $t_data_file $cm_per_pixel - mv -f allen_dump.qdp ${T_file} -elif [ $t_profile_type == zzl ] -then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_zzl_model $t_data_file $t_param_file - mv -f zzl_dump.qdp ${T_file} -else - echo temperature profile name invalid! - exit -fi - -cfunc_file=`grep '^cfunc_file' ${sbp_cfg} |awk '{print $2}'` -z=`grep '^z' ${sbp_cfg}|awk '{print $2}'` -abund=`grep '^abund' ${cfg_file} |awk '{print $2}'` -nh=`grep '^nh' ${cfg_file} |awk '{print $2}'` -$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file -mv flux_cnt_ratio.txt flux_cnt_ratio_center_entropy.txt -#fit sbp -$base_path/fit_dbeta_sbp $sbp_cfg -echo $cfunc_file -#exit - -#store central valu -mv sbp_fit.qdp sbp_fit_center_entropy.qdp -mv mass_int.qdp mass_int_center_entropy.qdp -mv overdensity.qdp overdensity_center_entropy.qdp -mv gas_mass_int.qdp gas_mass_int_center_entropy.qdp -mv entropy.qdp entropy_center.qdp -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 -fi - -cat ${radius_sbp_file} | sed 's/#.*$//' | grep -Ev '^\s*$' > .tmp.txt -mv .tmp.txt ${radius_sbp_file} - -rm -f summary_entropy.qdp - -#100 times of Monte-carlo simulation to determine error -#just repeat above steps -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 - -#exit - if [ $t_profile_type == zyy ] - then - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $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 - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_wang2012_model temp_shuffled_t.dat $t_param_file $cm_per_pixel - 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 - t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'` - $base_path/fit_zzl_model temp_shuffled_t.dat $t_param_file - mv -f zzl_dump.qdp ${T_file} - else - echo temperature profile name invalid! - exit - fi - -#exit - echo >temp_sbp.cfg - - cat $sbp_cfg|while read l -do - if echo $l|grep sbp_file >/dev/null - then - echo sbp_file temp_shuffled_sbp.dat >>temp_sbp.cfg - elif echo $l|grep T_file >/dev/null - then - echo T_file ${T_file} >>temp_sbp.cfg - else - echo $l >>temp_sbp.cfg - fi - -done - -#$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file - -$base_path/fit_dbeta_sbp temp_sbp.cfg - -cat entropy.qdp >>summary_entropy.qdp -echo no no no >>summary_entropy.qdp -done - -$base_path/analyze_entropy_profile.py $2 |tee entropy_result.txt - |