From 18e50108246fefb6062538df452b8fbf1d235731 Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Tue, 7 Feb 2017 13:57:52 +0800
Subject: Update tools to use the new style sbp config file

---
 mass_profile/beta_cfg.cpp      |  18 +++----
 mass_profile/beta_cfg.hpp      |   7 ++-
 mass_profile/calc_lx_beta.cpp  |  95 ++++++++++++++++++-------------------
 mass_profile/calc_lx_dbeta.cpp |  97 ++++++++++++++++++++------------------
 mass_profile/fit_beta_sbp.cpp  | 103 +++++++++++++++++++++--------------------
 mass_profile/fit_dbeta_sbp.cpp | 100 ++++++++++++++++++++-------------------
 6 files changed, 213 insertions(+), 207 deletions(-)

(limited to 'mass_profile')

diff --git a/mass_profile/beta_cfg.cpp b/mass_profile/beta_cfg.cpp
index 408970f..c373688 100644
--- a/mass_profile/beta_cfg.cpp
+++ b/mass_profile/beta_cfg.cpp
@@ -20,29 +20,23 @@ cfg_map parse_cfg_file(std::istream& is)
       string key;
       istringstream iss(line);
       iss>>key;
-      if(key=="radius_file")
+      if(key=="sbp_data")
 	{
 	  string value;
 	  iss>>value;
-	  result.radius_file=value;
+	  result.sbp_data=value;
 	}
-      else if(key=="sbp_file")
+      else if(key=="cfunc_profile")
 	{
 	  string value;
 	  iss>>value;
-	  result.sbp_file=value;
+	  result.cfunc_profile=value;
 	}
-      else if(key=="cfunc_file")
+      else if(key=="tprofile")
 	{
 	  string value;
 	  iss>>value;
-	  result.cfunc_file=value;
-	}
-      else if(key=="T_file")
-	{
-	  string value;
-	  iss>>value;
-	  result.T_file=value;
+	  result.tprofile=value;
 	}
       else if(key=="z")
 	{
diff --git a/mass_profile/beta_cfg.hpp b/mass_profile/beta_cfg.hpp
index 5dc9b2d..27148df 100644
--- a/mass_profile/beta_cfg.hpp
+++ b/mass_profile/beta_cfg.hpp
@@ -8,10 +8,9 @@
 
 struct cfg_map
 {
-  std::string radius_file;
-  std::string sbp_file;
-  std::string cfunc_file;
-  std::string T_file;
+  std::string sbp_data;
+  std::string cfunc_profile;
+  std::string tprofile;
   double z;
   double cm_per_pixel;
   double rmin_kpc;
diff --git a/mass_profile/calc_lx_beta.cpp b/mass_profile/calc_lx_beta.cpp
index aae8b2b..59ea51c 100644
--- a/mass_profile/calc_lx_beta.cpp
+++ b/mass_profile/calc_lx_beta.cpp
@@ -5,12 +5,11 @@
   This code is distributed with no warrant
 */
 
-#include <unistd.h>
 #include <iostream>
 #include <fstream>
+#include <sstream>
 #include <list>
 #include <algorithm>
-using namespace std;
 #include "beta_cfg.hpp"
 #include "dump_fit_qdp.hpp"
 #include "report_error.hpp"
@@ -23,18 +22,18 @@ using namespace std;
 #include <core/freeze_param.hpp>
 #include <error_estimator/error_estimator.hpp>
 #include "spline.h"
+
+using namespace std;
 using namespace opt_utilities;
 //double s=5.63136645E20;
 const double kpc=3.086E21;//kpc in cm
 const double Mpc=kpc*1000;
-double beta_func(double r,
-		 double n0,double rc,double beta)
-{
 
-  return abs(n0)*pow(1+r*r/rc/rc,-3./2.*abs(beta));
+double beta_func(double r, double n0, double rc, double beta)
+{
+  return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta));
 }
 
-
 //A class enclosing the spline interpolation method of cooling function
 //check spline.h for more detailed information
 //this class is a thin wrapper for the spline class defined in spline.h
@@ -82,8 +81,6 @@ int main(int argc,char* argv[])
   assert(cfg_file.is_open());
   cfg_map cfg=parse_cfg_file(cfg_file);
 
-  //check the existence of following parameters
-
   const double z=cfg.z;
 
   //initialize the radius list, sbp list and sbp error list
@@ -93,44 +90,49 @@ int main(int argc,char* argv[])
   std::vector<double> radii_all;
   std::vector<double> sbps_all;
   std::vector<double> sbpe_all;
+  //prepend the zero point to radius list
+  radii.push_back(0.0);
+  radii_all.push_back(0.0);
   //read sbp and sbp error data
-  for(ifstream ifs(cfg.sbp_file.c_str());;)
+  ifstream ifs(cfg.sbp_data.c_str());
+  std::string line;
+  if (ifs.is_open())
     {
-      assert(ifs.is_open());
-      double x,xe;
-      ifs>>x>>xe;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      if(x/xe<2)
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<xe<<endl;
-      sbps.push_back(x);
-      sbpe.push_back(xe);
-      sbps_all.push_back(x);
-      sbpe_all.push_back(xe);
+      while(std::getline(ifs, line))
+        {
+          if (line.empty())
+            continue;
+
+          std::istringstream iss(line);
+          double x, xe, y, ye;
+          if ((iss >> x >> xe >> y >> ye))
+            {
+              std::cerr << "sbprofile data: "
+                        << x << ", " << xe << ", " << y << ", " << ye
+                        << std::endl;
+              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
+              radii_all.push_back(x+xe);
+              sbps.push_back(y);
+              sbps_all.push_back(y);
+              sbpe.push_back(ye);
+              sbpe_all.push_back(ye);
+            }
+          else
+            {
+              std::cerr << "skipped line: " << line << std::endl;
+            }
+        }
     }
-
-  //read radius data
-  for(ifstream ifs(cfg.radius_file.c_str());;)
+  else
     {
-      assert(ifs.is_open());
-      double x;
-      ifs>>x;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<endl;
-      radii.push_back(x);
-      radii_all.push_back(x);
+      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
+                << std::endl;
+      return 1;
     }
+
   //initialize the cm/pixel value
   double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin=5*kpc/cm_per_pixel;
+  double rmin;
   if(cfg.rmin_pixel>0)
     {
       rmin=cfg.rmin_pixel;
@@ -140,7 +142,7 @@ int main(int argc,char* argv[])
       rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
     }
 
-  cerr<<"rmin="<<rmin<<endl;
+  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
   std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
   radii_tmp.resize(radii.size());
   sbps_tmp.resize(sbps.size());
@@ -169,7 +171,7 @@ int main(int argc,char* argv[])
 
   //read cooling function data
   spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_file.c_str());;)
+  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
     {
       assert(ifs.is_open());
       double x,y;
@@ -183,15 +185,14 @@ int main(int argc,char* argv[])
 	{
 	  break;
 	}
-      //cf.add_point(x,y*2.1249719395939022e-68);//change with source
-      cf.add_point(x,y);//change with source
+      cf.add_point(x,y);
     }
   cf.gen_spline();
 
   //read temperature profile data
   spline_func_obj Tprof;
   int tcnt=0;
-  for(ifstream ifs1(cfg.T_file.c_str());;++tcnt)
+  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
     {
       assert(ifs1.is_open());
       double x,y;
@@ -286,8 +287,8 @@ int main(int argc,char* argv[])
       cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
       param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
     }
-  cerr<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
 
   std::vector<double> mv=f.eval_model_raw(radii_all,p);
   int sbps_inner_cut_size=int(sbps_all.size()-sbps.size());
diff --git a/mass_profile/calc_lx_dbeta.cpp b/mass_profile/calc_lx_dbeta.cpp
index f8f2c2a..5fa2cb6 100644
--- a/mass_profile/calc_lx_dbeta.cpp
+++ b/mass_profile/calc_lx_dbeta.cpp
@@ -8,8 +8,8 @@
 
 #include <iostream>
 #include <fstream>
+#include <sstream>
 #include <list>
-using namespace std;
 #include "vchisq.hpp"
 #include "dbeta.hpp"
 #include "beta_cfg.hpp"
@@ -18,16 +18,19 @@ using namespace std;
 #include <core/freeze_param.hpp>
 #include <error_estimator/error_estimator.hpp>
 #include "spline.h"
+
+using namespace std;
 using namespace opt_utilities;
 //double s=5.63136645E20;
 const double kpc=3.086E21;//kpc in cm
 const double Mpc=kpc*1000;
-double dbeta_func(double r,
-		  double n01,double rc1,double beta1,
-		  double n02,double rc2,double beta2)
-{
 
-  return abs(n01)*pow(1+r*r/rc1/rc1,-3./2.*abs(beta1))+abs(n02)*pow(1+r*r/rc2/rc2,-3./2.*abs(beta2));
+double dbeta_func(double r, double n01, double rc1, double beta1,
+                  double n02, double rc2, double beta2)
+{
+  double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1));
+  double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2));
+  return v1 + v2;
 }
 
 
@@ -78,8 +81,6 @@ int main(int argc,char* argv[])
   assert(cfg_file.is_open());
   cfg_map cfg=parse_cfg_file(cfg_file);
 
-  //check the existence of following parameters
-
   const double z=cfg.z;
 
   //initialize the radius list, sbp list and sbp error list
@@ -89,44 +90,49 @@ int main(int argc,char* argv[])
   std::vector<double> radii_all;
   std::vector<double> sbps_all;
   std::vector<double> sbpe_all;
+  //prepend the zero point to radius list
+  radii.push_back(0.0);
+  radii_all.push_back(0.0);
   //read sbp and sbp error data
-  for(ifstream ifs(cfg.sbp_file.c_str());;)
+  ifstream ifs(cfg.sbp_data.c_str());
+  std::string line;
+  if (ifs.is_open())
     {
-      assert(ifs.is_open());
-      double x,xe;
-      ifs>>x>>xe;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      if(x/xe<2)
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<xe<<endl;
-      sbps.push_back(x);
-      sbpe.push_back(xe);
-      sbps_all.push_back(x);
-      sbpe_all.push_back(xe);
+      while(std::getline(ifs, line))
+        {
+          if (line.empty())
+            continue;
+
+          std::istringstream iss(line);
+          double x, xe, y, ye;
+          if ((iss >> x >> xe >> y >> ye))
+            {
+              std::cerr << "sbprofile data: "
+                        << x << ", " << xe << ", " << y << ", " << ye
+                        << std::endl;
+              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
+              radii_all.push_back(x+xe);
+              sbps.push_back(y);
+              sbps_all.push_back(y);
+              sbpe.push_back(ye);
+              sbpe_all.push_back(ye);
+            }
+          else
+            {
+              std::cerr << "skipped line: " << line << std::endl;
+            }
+        }
     }
-
-  //read radius data
-  for(ifstream ifs(cfg.radius_file.c_str());;)
+  else
     {
-      assert(ifs.is_open());
-      double x;
-      ifs>>x;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<endl;
-      radii.push_back(x);
-      radii_all.push_back(x);
+      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
+                << std::endl;
+      return 1;
     }
+
   //initialize the cm/pixel value
   double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin=5*kpc/cm_per_pixel;
+  double rmin;
   if(cfg.rmin_pixel>0)
     {
       rmin=cfg.rmin_pixel;
@@ -136,7 +142,7 @@ int main(int argc,char* argv[])
       rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
     }
 
-  cerr<<"rmin="<<rmin<<endl;
+  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
   std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
   radii_tmp.resize(radii.size());
   sbps_tmp.resize(sbps.size());
@@ -165,7 +171,7 @@ int main(int argc,char* argv[])
 
   //read cooling function data
   spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_file.c_str());;)
+  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
     {
       assert(ifs.is_open());
       double x,y;
@@ -179,15 +185,14 @@ int main(int argc,char* argv[])
 	{
 	  break;
 	}
-      //cf.add_point(x,y*2.1249719395939022e-68);//change with source
-      cf.add_point(x,y);//change with source
+      cf.add_point(x,y);
     }
   cf.gen_spline();
 
   //read temperature profile data
   spline_func_obj Tprof;
   int tcnt=0;
-  for(ifstream ifs1(cfg.T_file.c_str());;++tcnt)
+  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
     {
       assert(ifs1.is_open());
       double x,y;
@@ -372,8 +377,8 @@ int main(int argc,char* argv[])
       cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
       param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
     }
-  cerr<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
 
   //c.verbose(false);
   //f.set_statistic(c);
diff --git a/mass_profile/fit_beta_sbp.cpp b/mass_profile/fit_beta_sbp.cpp
index 9db48cc..c16c092 100644
--- a/mass_profile/fit_beta_sbp.cpp
+++ b/mass_profile/fit_beta_sbp.cpp
@@ -5,12 +5,11 @@
   This code is distributed with no warrant
 */
 
-#include <unistd.h>
 #include <iostream>
 #include <fstream>
+#include <sstream>
 #include <list>
 #include <algorithm>
-using namespace std;
 #include "beta_cfg.hpp"
 #include "dump_fit_qdp.hpp"
 #include "report_error.hpp"
@@ -23,19 +22,19 @@ using namespace std;
 #include <core/freeze_param.hpp>
 #include <error_estimator/error_estimator.hpp>
 #include "spline.h"
+
+using namespace std;
 using namespace opt_utilities;
 //double s=5.63136645E20;
 const double kpc=3.086E21;//kpc in cm
 const double Mpc=kpc*1000;
-double beta_func(double r,
-		 double n0,double rc,double beta)
-{
 
-  return abs(n0)*pow(1+r*r/rc/rc,-3./2.*abs(beta));
+double beta_func(double r, double n0, double rc, double beta)
+{
+  return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta));
 }
 
-
-  //calculate critical density from z, under following cosmological constants
+//calculate critical density from z, under following cosmological constants
 static double calc_critical_density(double z,
 				    const double H0=2.3E-18,
 				      const double Omega_m=.27)
@@ -94,8 +93,6 @@ int main(int argc,char* argv[])
   assert(cfg_file.is_open());
   cfg_map cfg=parse_cfg_file(cfg_file);
 
-  //check the existence of following parameters
-
   const double z=cfg.z;
 
   //initialize the radius list, sbp list and sbp error list
@@ -105,44 +102,50 @@ int main(int argc,char* argv[])
   std::vector<double> radii_all;
   std::vector<double> sbps_all;
   std::vector<double> sbpe_all;
+  //prepend the zero point to radius list
+  radii.push_back(0.0);
+  radii_all.push_back(0.0);
   //read sbp and sbp error data
-  for(ifstream ifs(cfg.sbp_file.c_str());;)
+  cerr << "Read surface brightness profile data ..." << endl;
+  ifstream ifs(cfg.sbp_data.c_str());
+  std::string line;
+  if (ifs.is_open())
     {
-      assert(ifs.is_open());
-      double x,xe;
-      ifs>>x>>xe;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      if(x/xe<2)
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<xe<<endl;
-      sbps.push_back(x);
-      sbpe.push_back(xe);
-      sbps_all.push_back(x);
-      sbpe_all.push_back(xe);
+      while(std::getline(ifs, line))
+        {
+          if (line.empty())
+            continue;
+
+          std::istringstream iss(line);
+          double x, xe, y, ye;
+          if ((iss >> x >> xe >> y >> ye))
+            {
+              std::cerr << "sbprofile data: "
+                        << x << ", " << xe << ", " << y << ", " << ye
+                        << std::endl;
+              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
+              radii_all.push_back(x+xe);
+              sbps.push_back(y);
+              sbps_all.push_back(y);
+              sbpe.push_back(ye);
+              sbpe_all.push_back(ye);
+            }
+          else
+            {
+              std::cerr << "skipped line: " << line << std::endl;
+            }
+        }
     }
-
-  //read radius data
-  for(ifstream ifs(cfg.radius_file.c_str());;)
+  else
     {
-      assert(ifs.is_open());
-      double x;
-      ifs>>x;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<endl;
-      radii.push_back(x);
-      radii_all.push_back(x);
+      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
+                << std::endl;
+      return 1;
     }
+
   //initialize the cm/pixel value
   double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin=5*kpc/cm_per_pixel;
+  double rmin;
   if(cfg.rmin_pixel>0)
     {
       rmin=cfg.rmin_pixel;
@@ -152,7 +155,7 @@ int main(int argc,char* argv[])
       rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
     }
 
-  cerr<<"rmin="<<rmin<<endl;
+  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
   std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
   radii_tmp.resize(radii.size());
   sbps_tmp.resize(sbps.size());
@@ -180,8 +183,9 @@ int main(int argc,char* argv[])
   copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
 
   //read cooling function data
+  cerr << "Read cooling function profile data ..." << endl;
   spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_file.c_str());;)
+  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
     {
       assert(ifs.is_open());
       double x,y;
@@ -193,17 +197,18 @@ int main(int argc,char* argv[])
       cerr<<x<<"\t"<<y<<endl;
       if(x>radii.back())
 	{
+      cerr << "radius_max: " << radii.back() << endl;
 	  break;
 	}
-      //cf.add_point(x,y*2.1249719395939022e-68);//change with source
-      cf.add_point(x,y);//change with source
+      cf.add_point(x,y);
     }
   cf.gen_spline();
 
   //read temperature profile data
+  cerr << "Read temperature profile data ..." << endl;
   spline_func_obj Tprof;
   int tcnt=0;
-  for(ifstream ifs1(cfg.T_file.c_str());;++tcnt)
+  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
     {
       assert(ifs1.is_open());
       double x,y;
@@ -221,8 +226,6 @@ int main(int argc,char* argv[])
 #endif
       Tprof.add_point(x,y);
     }
-
-
   Tprof.gen_spline();
 
   default_data_set<std::vector<double>,std::vector<double> > ds;
@@ -295,8 +298,8 @@ int main(int argc,char* argv[])
       cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
       param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
     }
-  cerr<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
 
   std::vector<double> mv=f.eval_model_raw(radii_all,p);
   int sbps_inner_cut_size=int(sbps_all.size()-sbps.size());
diff --git a/mass_profile/fit_dbeta_sbp.cpp b/mass_profile/fit_dbeta_sbp.cpp
index 7fb1c7d..30c4c35 100644
--- a/mass_profile/fit_dbeta_sbp.cpp
+++ b/mass_profile/fit_dbeta_sbp.cpp
@@ -8,8 +8,8 @@
 
 #include <iostream>
 #include <fstream>
+#include <sstream>
 #include <list>
-using namespace std;
 #include "vchisq.hpp"
 #include "dbeta.hpp"
 #include "beta_cfg.hpp"
@@ -18,20 +18,22 @@ using namespace std;
 #include <core/freeze_param.hpp>
 #include <error_estimator/error_estimator.hpp>
 #include "spline.h"
+
+using namespace std;
 using namespace opt_utilities;
 //double s=5.63136645E20;
 const double kpc=3.086E21;//kpc in cm
 const double Mpc=kpc*1000;
-double dbeta_func(double r,
-		  double n01,double rc1,double beta1,
-		  double n02,double rc2,double beta2)
-{
 
-  return abs(n01)*pow(1+r*r/rc1/rc1,-3./2.*abs(beta1))+abs(n02)*pow(1+r*r/rc2/rc2,-3./2.*abs(beta2));
+double dbeta_func(double r, double n01, double rc1, double beta1,
+                  double n02, double rc2, double beta2)
+{
+  double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1));
+  double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2));
+  return v1 + v2;
 }
 
-
-  //calculate critical density from z, under following cosmological constants
+//calculate critical density from z, under following cosmological constants
 static double calc_critical_density(double z,
 				    const double H0=2.3E-18,
 				      const double Omega_m=.27)
@@ -90,8 +92,6 @@ int main(int argc,char* argv[])
   assert(cfg_file.is_open());
   cfg_map cfg=parse_cfg_file(cfg_file);
 
-  //check the existence of following parameters
-
   const double z=cfg.z;
 
   //initialize the radius list, sbp list and sbp error list
@@ -101,44 +101,49 @@ int main(int argc,char* argv[])
   std::vector<double> radii_all;
   std::vector<double> sbps_all;
   std::vector<double> sbpe_all;
+  //prepend the zero point to radius list
+  radii.push_back(0.0);
+  radii_all.push_back(0.0);
   //read sbp and sbp error data
-  for(ifstream ifs(cfg.sbp_file.c_str());;)
+  ifstream ifs(cfg.sbp_data.c_str());
+  std::string line;
+  if (ifs.is_open())
     {
-      assert(ifs.is_open());
-      double x,xe;
-      ifs>>x>>xe;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      if(x/xe<2)
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<xe<<endl;
-      sbps.push_back(x);
-      sbpe.push_back(xe);
-      sbps_all.push_back(x);
-      sbpe_all.push_back(xe);
+      while(std::getline(ifs, line))
+        {
+          if (line.empty())
+            continue;
+
+          std::istringstream iss(line);
+          double x, xe, y, ye;
+          if ((iss >> x >> xe >> y >> ye))
+            {
+              std::cerr << "sbprofile data: "
+                        << x << ", " << xe << ", " << y << ", " << ye
+                        << std::endl;
+              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
+              radii_all.push_back(x+xe);
+              sbps.push_back(y);
+              sbps_all.push_back(y);
+              sbpe.push_back(ye);
+              sbpe_all.push_back(ye);
+            }
+          else
+            {
+              std::cerr << "skipped line: " << line << std::endl;
+            }
+        }
     }
-
-  //read radius data
-  for(ifstream ifs(cfg.radius_file.c_str());;)
+  else
     {
-      assert(ifs.is_open());
-      double x;
-      ifs>>x;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<endl;
-      radii.push_back(x);
-      radii_all.push_back(x);
+      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
+                << std::endl;
+      return 1;
     }
+
   //initialize the cm/pixel value
   double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin=5*kpc/cm_per_pixel;
+  double rmin;
   if(cfg.rmin_pixel>0)
     {
       rmin=cfg.rmin_pixel;
@@ -148,7 +153,7 @@ int main(int argc,char* argv[])
       rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
     }
 
-  cerr<<"rmin="<<rmin<<endl;
+  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
   std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
   radii_tmp.resize(radii.size());
   sbps_tmp.resize(sbps.size());
@@ -177,7 +182,7 @@ int main(int argc,char* argv[])
 
   //read cooling function data
   spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_file.c_str());;)
+  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
     {
       assert(ifs.is_open());
       double x,y;
@@ -191,15 +196,14 @@ int main(int argc,char* argv[])
 	{
 	  break;
 	}
-      //cf.add_point(x,y*2.1249719395939022e-68);//change with source
-      cf.add_point(x,y);//change with source
+      cf.add_point(x,y);
     }
   cf.gen_spline();
 
   //read temperature profile data
   spline_func_obj Tprof;
   int tcnt=0;
-  for(ifstream ifs1(cfg.T_file.c_str());;++tcnt)
+  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
     {
       assert(ifs1.is_open());
       double x,y;
@@ -377,8 +381,8 @@ int main(int argc,char* argv[])
       cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
       param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
     }
-  cerr<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
+  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
 
   //c.verbose(false);
   //f.set_statistic(c);
-- 
cgit v1.2.2