1 files changed, 51 insertions, 46 deletions

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