1 files changed, 48 insertions, 47 deletions

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());