diff options
author | Aaron LI <aaronly.me@gmail.com> | 2016-05-27 22:47:24 +0800 |
---|---|---|
committer | Aaron LI <aaronly.me@gmail.com> | 2016-05-27 22:47:24 +0800 |
commit | ffd178e0bd72562a3c2cff9747b6e656edc881dc (patch) | |
tree | 8800b7b5b2e8bc3df1a6760df5cd54eaaa686702 /mass_profile/fit_mt_bpl.cpp | |
parent | 5c35fad9240fb42c1371c721e0b2af7379bd9ea0 (diff) | |
download | chandra-acis-analysis-ffd178e0bd72562a3c2cff9747b6e656edc881dc.tar.bz2 |
Add mass_profile tools
* These tools are mainly use to calculate the total gravitational mass
profile, as well as the intermediate products (e.g., surface
brightness profile fitting, gas density profile, NFW fitting, etc.)
* There are additional tools for calculating the luminosity and flux.
* These tools mainly developed by Junhua GU, and contributed by
Weitian (Aaron) LI, and Zhenghao ZHU.
Diffstat (limited to 'mass_profile/fit_mt_bpl.cpp')
-rw-r--r-- | mass_profile/fit_mt_bpl.cpp | 350 |
1 files changed, 350 insertions, 0 deletions
diff --git a/mass_profile/fit_mt_bpl.cpp b/mass_profile/fit_mt_bpl.cpp new file mode 100644 index 0000000..7cf22b2 --- /dev/null +++ b/mass_profile/fit_mt_bpl.cpp @@ -0,0 +1,350 @@ +/* + Perform a double-beta density model fitting to the surface brightness data + Author: Junhua Gu + Last modified: 2011.01.01 + This code is distributed with no warrant +*/ + +//#define HAVE_X_ERROR +#include <iomanip> +#include <iostream> +#include <sstream> +#include <fstream> +#include <models/bpl1d.hpp> +#include <models/lin1d.hpp> +#include "statistics/chisq.hpp" +#include "statistics/logchisq.hpp" +#include "statistics/leastsq.hpp" +#include <data_sets/default_data_set.hpp> +#include <methods/powell/powell_method.hpp> +#include <methods/gsl_simplex/gsl_simplex.hpp> +#include <core/freeze_param.hpp> + +using namespace std; +using namespace opt_utilities; +//double s=5.63136645E20; +const double kpc=3.086E21;//kpc in cm +const double Mpc=kpc*1000; +const double pi=4*atan(1); +double std_norm_rand() +{ + double u=0; + double v=0; + while(u<=0||v<=0) + { + u=rand()/(double)RAND_MAX; + rand(); + v=rand()/(double)RAND_MAX; + } + double x=std::sqrt(-log(u))*cos(2*pi*v); + return x; +} + +double shuffle_data(double xc,double xl,double xu) +{ + double lxc=log(xc); + double lxl=log(xc-xl)-log(xc); + double lxu=log(xc+xu)-log(xc); + + if(std_norm_rand()>0) + { + double result=std::exp(lxc-std::abs(std_norm_rand()*lxl)); + return result; + } + else + { + double result=std::exp(lxc+std::abs(std_norm_rand()*lxu)); + return result; + } +} + +int main(int argc,char* argv[]) +{ + srand(time(0)); + if(argc!=4) + { + cerr<<"Usage:"<<argv[0]<<" <a 6 column file with T -Terr +Terr M -Merr +Merr> <initial broken temperature> <T lower limit>"<<endl; + return -1; + } + double T_lower_limit(atof(argv[3])); + double Tb=atof(argv[2]); + assert(Tb>0); + ifstream ifs_data(argv[1]); + default_data_set<double,double> ds; + ofstream ofs_result("m-t_bpl_result.qdp"); + ofs_result<<"read terr 1 2"<<endl; + ofs_result<<"skip single"<<endl; + ofs_result<<"log"<<endl; + //ofs_result<<"li on 2"<<endl; + ofs_result<<"time off"<<endl; + ofs_result<<"la f"<<endl; + ofs_result<<"la x temperature (keV)"<<endl; + ofs_result<<"la y mass (M\\d\\(2281)\\u)"<<endl; + double sxxl=0; + double s1l=0; + double sxl=0; + double syl=0; + double sxyl=0; + + double sxxu=0; + double s1u=0; + double sxu=0; + double syu=0; + double sxyu=0; + + bool is_first_nonono=true; + + for(;;) + { + double T,Tl,Tu; + double M,Ml,Mu; + std::string line; + getline(ifs_data,line); + //ifs_data>>T>>Tl>>Tu>>M>>Ml>>Mu; + if(!ifs_data.good()) + { + break; + } + line+=" "; + istringstream iss(line); + + if(line[0]=='#') + { + if(!is_first_nonono) + { + ofs_result<<"no no no"<<endl; + } + else + { + is_first_nonono=false; + } + continue; + } + iss>>T>>Tl>>Tu>>M>>Ml>>Mu; + //std::cerr<<L<<"\t"<<Lerr<<endl; + if(!iss.good()) + { + continue; + } + + if(T<T_lower_limit||M<0) + { + continue; + } + + if(std::abs(Mu)<M*.1||std::abs(Ml)<M*.1) + { + cerr<<"mass error less than 10%, skipped"<<endl; + cerr<<line<<endl; + continue; + } +#if 1 + if(std::abs(Tu)<.1||std::abs(Tl)<.1) + { + cerr<<"T error less than 10%, skipped"<<endl; + cerr<<line<<endl; + continue; + } +#endif + + if(std::abs(Mu)+std::abs(Ml)<M*.1) + { + double k=M*.1/(std::abs(Mu)+std::abs(Ml)); + Mu*=k; + Ml*=k; + } + if(std::abs(Ml)>std::abs(M)) + { + continue; + } + Tl=std::abs(Tl); + Tu=std::abs(Tu); + Ml=std::abs(Ml); + Mu=std::abs(Mu); + ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<M<<"\t"<<-std::abs(Ml)<<"\t"<<+std::abs(Mu)<<endl; + double x=(T); + double y=(M); + double xu=Tu; + double xl=Tl; + + double yu=Mu; + double yl=Ml; + if(T>Tb) + { + sxxl+=log(x)*log(x); + sxl+=log(x); + syl+=log(y); + sxyl+=log(y)*log(x); + s1l+=1; + } + else + { + sxxu+=log(x)*log(x); + sxu+=log(x); + syu+=log(y); + sxyu+=log(y)*log(x); + s1u+=1; + } + data<double,double> d(x,y,std::abs(yl),std::abs(yu), + std::abs(xl),std::abs(xu)); + ds.add_data(d); + } + + double Ml=sxxl*s1l-sxl*sxl; + double Mal=sxyl*s1l-syl*sxl; + double Mbl=sxxl*syl-sxl*sxyl; + double k0l=Mal/Ml; + double b0l=Mbl/Ml; + + double Mu=sxxu*s1u-sxu*sxu; + double Mau=sxyu*s1u-syu*sxu; + double Mbu=sxxu*syu-sxu*sxyu; + double k0u=Mau/Mu; + double b0u=Mbu/Mu; + + double gamma0l=k0l; + double gamma0u=k0u; + + double ampl0l=exp(b0l)*pow(Tb,gamma0l); + double ampl0u=exp(b0u)*pow(Tb,gamma0u);; + + + + ofs_result<<"no no no"<<endl; + fitter<double,double,vector<double>,double,std::string> fit; + fit.set_opt_method(powell_method<double,vector<double> >()); + + fit.set_statistic(logchisq<double,double,vector<double>,double,std::string>()); + //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>()); + fit.set_model(bpl1d<double>()); + fit.load_data(ds); + + cerr<<"k0l="<<k0l<<endl; + cerr<<"k0u="<<k0u<<endl; + cerr<<"Ampl0="<<(ampl0l+ampl0u)/2<<endl; + + fit.set_param_value("bpx",Tb); + fit.set_param_value("bpy",(ampl0l+ampl0u)/2); + fit.set_param_value("gamma1",gamma0l); + fit.set_param_value("gamma2",gamma0u); + + fit.fit(); + //fit.set_opt_method(gsl_simplex<double,vector<double> >()); + fit.fit(); + std::vector<double> p=fit.fit(); + Tb=fit.get_param_value("bpx"); + //std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl; + for(double i=.5;i<12;i*=1.01) + { + ofs_result<<i<<"\t0\t0\t"<<fit.eval_model_raw(i,p)<<"\t0\t0\n"; + } + + + std::vector<double> mean_p(p.size()); + std::vector<double> mean_p2(p.size()); + int cnt=0; + for(int n=0;n<100;++n) + { + ++cnt; + cerr<<"."; + double sxxl=0; + double s1l=0; + double sxl=0; + double syl=0; + double sxyl=0; + + double sxxu=0; + double s1u=0; + double sxu=0; + double syu=0; + double sxyu=0; + + opt_utilities::default_data_set<double,double> ds1; + for(int i=0;i<ds.size();++i) + { + double x=ds.get_data(i).get_x(); + double y=ds.get_data(i).get_y(); + double xl=ds.get_data(i).get_x_lower_err(); + double xu=ds.get_data(i).get_x_upper_err(); + double yl=ds.get_data(i).get_y_lower_err(); + double yu=ds.get_data(i).get_y_upper_err(); + + double new_x=shuffle_data(x, + xl, + xu); + double new_y=shuffle_data(y, + yl, + yu); + + ds1.add_data(data<double,double>(new_x,new_y, + yl/y*new_y, + yu/y*new_y, + xl/x*new_x, + xu/x*new_x)); + + x=new_x; + y=new_y; + if(x>Tb) + { + sxxl+=log(x)*log(x); + sxl+=log(x); + syl+=log(y); + sxyl+=log(y)*log(x); + s1l+=1; + } + else + { + sxxu+=log(x)*log(x); + sxu+=log(x); + syu+=log(y); + sxyu+=log(y)*log(x); + s1u+=1; + } + } + double Ml=sxxl*s1l-sxl*sxl; + double Mal=sxyl*s1l-syl*sxl; + double Mbl=sxxl*syl-sxl*sxyl; + double k0l=Mal/Ml; + double b0l=Mbl/Ml; + + double Mu=sxxu*s1u-sxu*sxu; + double Mau=sxyu*s1u-syu*sxu; + double Mbu=sxxu*syu-sxu*sxyu; + double k0u=Mau/Mu; + double b0u=Mbu/Mu; + + double gamma0l=k0l; + double gamma0u=k0u; + + double ampl0l=exp(b0l)*pow(Tb,gamma0l); + double ampl0u=exp(b0u)*pow(Tb,gamma0u);; + + fit.set_param_value("bpx",Tb); + fit.set_param_value("bpy",(ampl0l+ampl0u)/2); + fit.set_param_value("gamma1",gamma0l); + fit.set_param_value("gamma2",gamma0u); + + + fit.load_data(ds1); + + fit.fit(); + vector<double> p=fit.fit(); + for(int i=0;i<p.size();++i) + { + mean_p[i]+=p[i]; + mean_p2[i]+=p[i]*p[i]; + } + //cerr<<fit.get_param_value("gamma1")<<"\t"<<fit.get_param_value("gamma2")<<endl; + + } + vector<double> std_p(p.size()); + cerr<<endl; + for(int i=0;i<mean_p.size();++i) + { + mean_p[i]/=cnt; + mean_p2[i]/=cnt; + std_p[i]=std::sqrt(mean_p2[i]-mean_p[i]*mean_p[i]); + cout<<fit.get_param_info(i).get_name()<<"= "<<p[i]<<" +/- "<<std_p[i]<<endl; + } + std::cout<<"Num of sources:"<<ds.size()<<endl; +} |