LoadSassenaParams.cpp

/*WIKI*

This algorithm loads a Sassena output file into a group workspace.
It will create a workspace for each scattering intensity and one workspace for the Q-values

*WIKI*/

#include<fstream>
#include <utility>
#include <bits/stl_pair.h>

#include "MantidDataHandling/LoadSassenaParams.h"
#include "MantidAPI/LoadAlgorithmFactory.h"
#include "MantidAPI/FileProperty.h"
#include "MantidKernel/Exception.h"

// special library headers
#include <boost/filesystem.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/regex.hpp>
#include <boost/program_options.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_on_sphere.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/thread.hpp>

namespace Mantid
{
namespace DataHandling
{

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(LoadSassenaParams)

// Initialize the loggers
Kernel::Logger& XMLInterface::g_log = Kernel::Logger::get("XMLInterface");
Kernel::Logger& ScatteringAverageOrientationVectorsParameters::g_log = Kernel::Logger::get("ScatteringAverageOrientationVectorsParameters");
Kernel::Logger& ScatteringAverageOrientationMultipoleMomentsParameters::g_log = Kernel::Logger::get("ScatteringAverageOrientationMultipoleMomentsParameters");
Kernel::Logger& ScatteringVectorsParameters::g_log = Kernel::Logger::get("ScatteringVectorsParameters");
Kernel::Logger& Params::g_log = Kernel::Logger::get("Params");

/**Constructor
 * @param ptr xml node, possibly containing children
 */
XMLElement::XMLElement(xmlNodePtr ptr) {
  node_ptr = ptr;
  if (node_ptr->children!=NULL) {
    xmlNodePtr cur = node_ptr->children;
    while(cur!=NULL) {
      children.push_back(cur);
      cur=cur->next;
    }
  }
}

void XMLElement::print(const std::string prepend, bool showchildren)
{
  std::cout << prepend << "type   = " << this->type() << std::endl;
  std::cout << prepend << "name   = " << this->name() << std::endl;
  std::cout << prepend << "content= " << this->content() << std::endl;
  if (!showchildren) return;
  if (this->children.size()==0) return;
  std::cout << prepend << "children: [" << this->children.size() << "]" << std::endl;
  for(size_t i = 0; i < this->children.size(); ++i){
    std::cout << prepend << "-----" << std::endl;
    this->children[i].print(std::string(prepend)+"  ", showchildren);
    std::cout << prepend << "-----" << std::endl;
  }
}

XMLInterface::XMLInterface(const std::string filename)
{
  xmlInitParser();
  try
  {
    p_doc = xmlParseFile(filename.c_str());
  }
  catch(xmlError const& p_excp)
  {
    throw Kernel::Exception::FileError(std::string(p_excp.message)+". Unable to parse File:", filename);
  }

  /* Create xpath evaluation context */
  p_xpathCtx = xmlXPathNewContext(p_doc);
  if(p_xpathCtx == NULL)
  {
    xmlFreeDoc(p_doc);
    throw Kernel::Exception::NullPointerException("p_xpathCtx","XMLInterface::XMLInterface");
  }

  maxrecursion = 20;
  size_t recursion = 1;
  while( xmlXIncludeProcess(p_doc)!=0 )
  {
    if (recursion >maxrecursion)
      throw Kernel::Exception::IndexError(recursion,maxrecursion,"XMLInterface::XMLInterface");
    recursion++;
  }
} // end of XMLInterface::XMLInterface

void XMLInterface::dump(std::vector<char>& c)
{
  int chars;
  xmlChar * data;
  xmlDocDumpMemory(p_doc,&data,&chars);
  for(int i = 0; i < chars; ++i) c.push_back(data[i]);
  xmlFree(data);
}

XMLInterface::~XMLInterface()
{
  xmlXPathFreeContext(p_xpathCtx);
  xmlFreeDoc(p_doc);
  xmlCleanupParser();
}

std::vector<XMLElement> XMLInterface::get(const char* xpathexp)
{
  return get(std::string(xpathexp));
}

std::vector<XMLElement> XMLInterface::get(std::string xpathexp)
{
  // PREPARE XPATH expressions first...
  const xmlChar* p_xpath_exp = xmlCharStrdup(xpathexp.c_str());
  xmlXPathObjectPtr p_xpathObj = xmlXPathEvalExpression(p_xpath_exp, p_xpathCtx);
  if (p_xpathObj==NULL)
  {
    g_log.error("XPath expression not found: "+xpathexp);
    throw;
  }
  xmlNodeSetPtr nodes = p_xpathObj->nodesetval; // nodes
  size_t size = (nodes) ? nodes->nodeNr : 0;
  std::vector<XMLElement> result;
  for(size_t i = 0; i < size; ++i)
  {
    result.push_back(nodes->nodeTab[i]);
  }
  xmlXPathFreeObject(p_xpathObj);
  return result;
}

bool XMLInterface::exists(const char* xpathexp)
{
  const xmlChar* p_xpath_exp = xmlCharStrdup(xpathexp);
  xmlXPathObjectPtr p_xpathObj = xmlXPathEvalExpression(p_xpath_exp, p_xpathCtx);
  if (p_xpathObj->nodesetval->nodeNr==0)
  {
  return false;
  }
  else
  {
    return true;
  }
}

template<> bool XMLInterface::get_value<bool>(const char* xpathexp)
{
  std::vector<XMLElement> elements = get(xpathexp);
  // if elements has more than one entry, then xpathexp is ambigious
  if (elements.size()==0)
  {
    return boost::lexical_cast<bool>("");
  }
  else if (elements.size()>1)
  {
    g_log.error("Xpathexp is ambigious, multiple fields matched: "+std::string(xpathexp));
    throw;
  }
  XMLElement& thisel = elements[0];
  if (thisel.type()!=XML_ELEMENT_NODE)
  {
    g_log.error("Xpathexp doesn't resolve to XML_ELEMENT_NODE: "+std::string(xpathexp));
  throw;
  }
  size_t textelements = 0;
  std::string result;
  for(size_t i = 0; i < thisel.children.size(); ++i)
  {
    if (thisel.children[0].type()==XML_TEXT_NODE)
    {
      textelements++;
      result = thisel.children[0].content();
    }
  }
  // textelements == 0 corresponds to an empty text field. That is legal!
  if (textelements>1)
  {
    g_log.error("Xpathexp resolves to more than one text field: "+std::string(xpathexp));
    throw;
  }
  // trim whitespaces!
  boost::trim(result);
  std::string result_upper = result;
  boost::to_upper(result_upper);
  // circumvent differences b/w lexical_cast and xml specification
  if (result_upper=="FALSE") result = "0";
  if (result_upper=="TRUE") result = "1";
  return boost::lexical_cast<bool>(result);
}

/// Provides type safe access to content element in XML files through XPATH.
template<class convT> convT XMLInterface::get_value(const char* xpathexp)
{
  std::vector<XMLElement> elements = get(xpathexp);
  // if elements has more than one entry, then xpathexp is ambigious
  if (elements.empty())
  {
    return boost::lexical_cast<convT>("");
  }
  else if (elements.size()>1)
  {
    g_log.error("Xpathexp is ambigious, multiple fields matched: "+std::string(xpathexp));
    throw;
  }
  XMLElement& thisel = elements[0];
  if (thisel.type()!=XML_ELEMENT_NODE)
  {
    g_log.error("Xpathexp doesn't resolve to XML_ELEMENT_NODE: "+std::string(xpathexp));
    throw;
  }
  size_t textelements = 0;
  std::string result;
  for(size_t i = 0; i < thisel.children.size(); ++i)
  {
    if (thisel.children[0].type()==XML_TEXT_NODE)
    {
      textelements++;
      result = thisel.children[0].content();
    }
  }
  // textelements == 0 corresponds to an empty text field. That is legal!
  if (textelements>1)
  {
    g_log.error("Xpathexp resolves to more than one text field: "+std::string(xpathexp));
    throw;
   }
  // trim whitespaces!
  boost::trim(result);
  // circumvent differences b/w lexical_cast and xml specification
  return boost::lexical_cast<convT>(result);
}

/**
 *  Overloading operator '<' to compare two CartesianCoor3D vectors by their components
 * @param that right hand side
 */
bool CartesianCoor3D::operator<(const CartesianCoor3D& that) const
{
  if (this->x < that.x)
    return true;

  if (this->x > that.x)
    return false;

  if (this->y < that.y)
    return true;

  if (this->y > that.y)
    return false;

  if (this->z < that.z)
    return true;

  if (this->z > that.z)
    return false;

  return false;
}

/**
 * Conversion constructor cylinder -> cartesian
 * @param cc a CylinderCoor3D vector
 */
CartesianCoor3D::CartesianCoor3D(CylinderCoor3D cc)
{
  x = cc.r*cos(cc.phi);
  y = cc.r*sin(cc.phi);
  z = cc.z;
}

/**
 * Conversion constructor spherical -> cartesian
 * @param cc a SphericalCoor3D vector
 */
CartesianCoor3D::CartesianCoor3D(SphericalCoor3D cc)
{
  x = cc.r*sin(cc.theta)*cos(cc.phi);
  y = cc.r*sin(cc.theta)*sin(cc.phi);
  z = cc.r*cos(cc.theta);
}

/// return the size of the vector
double CartesianCoor3D::length()
{
  return sqrt(pow(x,2)+pow(y,2)+pow(z,2));
}


std::ostream& operator<<(std::ostream& os, const CartesianCoor3D& cc)
{
  return os << "(x=" << cc.x << ",y=" << cc.y << ",z=" << cc.z << ")";
}


CartesianCoor3D& CartesianCoor3D::operator=(const CartesianCoor3D& that)
{
  if (this != &that) { x=that.x; y=that.y; z=that.z; }
  return *this;
}

CartesianCoor3D CartesianCoor3D::operator-(const CartesianCoor3D& that)
{
  return CartesianCoor3D(x-that.x,y-that.y,z-that.z);
}

CartesianCoor3D CartesianCoor3D::operator+(const CartesianCoor3D& that)
{
  return CartesianCoor3D(x+that.x,y+that.y,z+that.z);
}

double CartesianCoor3D::operator*(const CartesianCoor3D& that)
{
  return (x*that.x+y*that.y+z*that.z);
}

CartesianCoor3D CartesianCoor3D::cross_product(const CartesianCoor3D& that)
{
  return CartesianCoor3D(y*that.z-z*that.y,z*that.x-x*that.z,x*that.y-y*that.x);
}

CartesianCoor3D operator*(const double lambda, const CartesianCoor3D& that)
{
  return CartesianCoor3D(lambda*that.x,lambda*that.y,lambda*that.z);
}

CartesianCoor3D operator*(const CartesianCoor3D& that,const double lambda)
{
  return CartesianCoor3D(lambda*that.x,lambda*that.y,lambda*that.z);
}

CartesianCoor3D operator/(const CartesianCoor3D& that, const double lambda)
{
  return CartesianCoor3D(that.x/lambda,that.y/lambda,that.z/lambda);
}

/**
 * Conversion constructor cartesian -> cylinder
 */
CylinderCoor3D::CylinderCoor3D(double v1,double v2,double v3)
{
  if (v1<0){
    throw "Negative values not allowed for cylinder radius\n";
  }
  r   = v1;
  phi = v2;
  z   = v3;
}

/**
 * Conversion constructor cartesian -> cylinder
 */
CylinderCoor3D::CylinderCoor3D(CartesianCoor3D cc)
{
  float M_PIf = static_cast<float>(M_PI);
  float M_PI_2f = static_cast<float>(M_PI_2);
  float ccyf = static_cast<float>(cc.y);
  r = sqrt(pow(cc.x,2)+pow(cc.y,2));

  if (cc.x!=0.0)
  {
    phi = atan(cc.y/cc.x);
    if (cc.x<0.0)
    {
      phi = static_cast<double>( sign(M_PIf,ccyf) )+phi;
    }
  }
  else if (cc.y!=0.0)
  {
    phi = static_cast<double>( sign(M_PI_2f,ccyf) );
  }
  else
  {
    phi = 0.0;
  }

  phi = phi<0 ? 2*M_PI + phi : phi;

  if ( (phi<0) || (phi>2*M_PI))
  {
    std::cerr << "PHI OUT OF BOUNDS: " << r << ", " << phi << std::endl;
    std::cerr << cc << std::endl;
    throw;
  }
  z = cc.z;
}

/**
 * Conversion constructor spherical -> cylinder
 */
CylinderCoor3D::CylinderCoor3D(SphericalCoor3D cc)
{
  CartesianCoor3D c1(cc);
  CylinderCoor3D c2(c1);
  r = c2.r; phi = c2.phi; z = c2.z;
}

std::ostream& operator<<(std::ostream& os, const CylinderCoor3D& cc)
{
  return os << "(r=" << cc.r << ",phi=" << cc.phi << ",z=" << cc.z << ")";
}

CylinderCoor3D& CylinderCoor3D::operator=(const CylinderCoor3D& that)
{
  if (this != &that) { r=that.r; phi=that.phi; z=that.z; }
  return *this;
}

CylinderCoor3D CylinderCoor3D::operator-(const CylinderCoor3D& that)
{
  return CylinderCoor3D(r-that.r,phi-that.phi,z-that.z);
}

/**
 * Conversion constructor cartesian -> spherical
 */
SphericalCoor3D::SphericalCoor3D(CartesianCoor3D cc)
{
  r = sqrt(pow(cc.x,2)+pow(cc.y,2)+pow(cc.z,2));
  if (r==0)
  {
    theta = 0;
    phi=0;
    return;
  }
  theta = acos(cc.z/r);
  if (theta<0) throw;
  if (theta>M_PI) throw;
  if (cc.x!=0.0)
  {
    phi = atan(cc.y/cc.x);
    if (cc.x<0.0) phi += M_PI;
    else if (cc.y<0.0) phi += 2*M_PI;
  }
  else if (cc.y!=0.0)
  {
    if (cc.y>0) phi=M_PI_2;
    if (cc.y<0) phi=3*M_PI_2;
  }
  else
  {
    phi = 0.0;
  }
  if (phi<0) throw;
  if (phi>2*M_PI) throw;
}

/**
 * Conversion constructor cylinder -> spherical
 */
SphericalCoor3D::SphericalCoor3D(CylinderCoor3D cc)
{
  CartesianCoor3D c1(cc);
  SphericalCoor3D c2(c1);
  r = c2.r; phi = c2.phi; theta = c2.theta;
}

std::ostream& operator<<(std::ostream& os, const SphericalCoor3D& cc)
{
  return os << "(r=" << cc.r << ",phi=" << cc.phi << ",theta=" << cc.theta << ")";
}

SphericalCoor3D& SphericalCoor3D::operator=(const SphericalCoor3D& that)
{
  if (this != &that) { r=that.r; phi=that.phi; theta=that.theta;}
  return *this;
}

SphericalCoor3D SphericalCoor3D::operator-(const SphericalCoor3D& that)
{
  return SphericalCoor3D(r-that.r,phi-that.phi,theta-that.theta);
}

//transformations

CartesianCoor3D rotate(CartesianCoor3D c,std::string axis,double rad)
{
  CartesianCoor3D r;
  if (axis=="x") {
    double m[3][3];
    m[0][0] = 1; m[0][1] = 0; m[0][2] = 0;
    m[1][0] = 0; m[1][1] = cos(rad); m[1][2] = -sin(rad);
    m[2][0] = 0; m[2][1] = sin(rad); m[2][2] = cos(rad);
    r.x = m[0][0]*c.x + m[0][1]*c.y + m[0][2]*c.z;
    r.y = m[1][0]*c.x + m[1][1]*c.y + m[1][2]*c.z;
    r.z = m[2][0]*c.x + m[2][1]*c.y + m[2][2]*c.z;
  }
  if (axis=="y") {
    double m[3][3];
    m[0][0] = cos(rad); m[0][1] = 0; m[0][2] = sin(rad);
    m[1][0] = 0; m[1][1] = 1; m[1][2] = 0;
    m[2][0] = -sin(rad); m[2][1] = 0; m[2][2] = cos(rad);
    r.x = m[0][0]*c.x + m[0][1]*c.y + m[0][2]*c.z;
    r.y = m[1][0]*c.x + m[1][1]*c.y + m[1][2]*c.z;
    r.z = m[2][0]*c.x + m[2][1]*c.y + m[2][2]*c.z;
  }
  if (axis=="z") {
    double m[3][3];
    m[0][0] = cos(rad); m[0][1] = -sin(rad); m[0][2] = 0;
    m[1][0] = sin(rad); m[1][1] = cos(rad); m[1][2] = 0;
    m[2][0] = 0; m[2][1] = 0; m[2][2] = 1;
    r.x = m[0][0]*c.x + m[0][1]*c.y + m[0][2]*c.z;
    r.y = m[1][0]*c.x + m[1][1]*c.y + m[1][2]*c.z;
    r.z = m[2][0]*c.x + m[2][1]*c.y + m[2][2]*c.z;
  }
  return r;
}

CartesianVectorBase::CartesianVectorBase(CartesianCoor3D axis)
{
  CartesianCoor3D ek(0,0,1);
  CartesianCoor3D ej(0,1,0);
  CartesianCoor3D ez = axis/axis.length();
  CartesianCoor3D ekez = ek.cross_product(ez);
  if (ekez.length()==0)
  {
    ekez = ej.cross_product(ez);
  }
  CartesianCoor3D er = ekez/ekez.length();
  CartesianCoor3D ezer = ez.cross_product(er);
  CartesianCoor3D ephi = ezer/ezer.length();
  base_.push_back(er);
  base_.push_back(ephi);
  base_.push_back(ez);
}

CartesianCoor3D& CartesianVectorBase::operator[](size_t index)
{
  return base_.at(index);
}

/// Project the vector vec onto the base and return the 3 components along the three basis vectors
CartesianCoor3D CartesianVectorBase::project(CartesianCoor3D vec)
{
  return CartesianCoor3D(vec*base_[0],vec*base_[1],vec*base_[2]);
}

SampleParameters::~SampleParameters()
{
  for(std::map<std::string,SampleSelectionParameters*>::iterator i=selections.begin();i!=selections.end();++i)
  {
    delete i->second;
  }
}

void ScatteringAverageOrientationVectorsParameters::create()
{
  if (type=="file")
  {
    g_log.information("Reading orientations for orientational averaging from file: "+this->file);
    std::ifstream qqqfile( file.c_str() );
    double x,y,z;
    while (qqqfile >> x >> y >> z)
    {
      CartesianCoor3D q(x,y,z);
      double ql = q.length();
      if (ql!=0) q = (1.0/ql) * q;
      this->push_back(q);
    }
  }
  else if (type=="sphere")
  {
    if (algorithm=="boost_uniform_on_sphere")
    {
      boost::mt19937 rng; // that's my random number generator
      rng.seed( static_cast<unsigned int>(seed) );
      boost::uniform_on_sphere<double> s(3); // that's my distribution
      boost::variate_generator<boost::mt19937, boost::uniform_on_sphere<double> > mysphere(rng,s);
      g_log.information("Generating orientations for orientational averaging using sphere,boost_uniform_on_sphere");
      for(size_t i = 0; i < resolution; ++i)
      {
        vector<double> r = mysphere();
        this->push_back( CartesianCoor3D(r[0],r[1],r[2]) );
      }
    }
    else
    {
      g_log.error("Vectors algorithm not understood: "+this->algorithm);
      throw;
      }
    }
  else if (type=="cylinder")
  {
    if (algorithm=="boost_uniform_on_sphere")
    {
      boost::mt19937 rng; // that's my random number generator
      rng.seed( static_cast<unsigned int>(seed) );
      boost::uniform_on_sphere<double> s(2); // that's my distribution
      boost::variate_generator<boost::mt19937, boost::uniform_on_sphere<double> > mysphere(rng,s);
      g_log.information("Generating orientations for orientational averaging using cylinder,boost_uniform_on_sphere");
      for(size_t i = 0; i < resolution; ++i)
      {
        vector<double> r = mysphere();
        this->push_back( CartesianCoor3D(r[0],r[1],0) );
      }
    }
    else if (algorithm=="raster_linear")
    {
      const double M_2PI = 2*M_PI;
      const double radincr = (M_2PI) / ( 360.0*static_cast<double>(resolution) );
      g_log.information("Generating orientations for orientational averaging using cylinder,raster_linear");
      for (double phi=0;phi<M_2PI;phi+=radincr)
      {
        this->push_back( CartesianCoor3D(cos(phi),sin(phi),0) );
      }
    }
    else
      {
        g_log.error("Vectors algorithm not understood: "+this->algorithm);
        throw;
      }
    }
  else
  {
    g_log.error("Vectors orientation averaging type not understood: "+this->type);
    throw;
  }
  g_log.information("Initialized orientational averaging with "+boost::lexical_cast<std::string>(this->size())+" vectors.");
}// end of ScatteringAverageOrientationVectorsParameters::create()

void ScatteringAverageOrientationMultipoleMomentsParameters::create()
{
  if(type=="file")
  {
    g_log.information("Reading multipole moments for orientational averaging from file: "+this->file);
    std::ifstream mmfile(this->file.c_str());
    long major, minor;
    while (mmfile >> major >> minor)
    {
      this->push_back(std::make_pair(major,minor));
    }
  }
  else if(type=="resolution")
  {
    g_log.information("Generating multipole moments for orientational averaging using a maxium major="+boost::lexical_cast<std::string>(resolution));
    if (Params::Inst()->scattering.average.orientation.multipole.type=="sphere")
    {
      this->push_back(std::make_pair(0L,0L));
      for(long l = 1; l <= resolution; ++l)
      {
        for(long m = -l; m <= l; ++m)
        {
          this->push_back(std::make_pair(l,m));
        }
      }
    }
    else if(Params::Inst()->scattering.average.orientation.multipole.type=="cylinder")
    {
      this->push_back(std::make_pair(0L,0L));
      for(long l = 1; l <= resolution; ++l)
      {
        this->push_back(std::make_pair(l,0L));
        this->push_back(std::make_pair(l,1L));
        this->push_back(std::make_pair(l,2L));
        this->push_back(std::make_pair(l,3L));
      }
    }
    else
    {
      g_log.error("Type not understood: scattering.average.orientation.multipole.moments.type="+type);
      throw;
    }
  }
  else
  {
    g_log.error("Type not understood: scattering.average.orientation.multipole.type="+Params::Inst()->scattering.average.orientation.multipole.type);
        throw;
    }

  // check validity of moments
  g_log.information("Checking multipole moments.");
  if (Params::Inst()->scattering.average.orientation.multipole.type=="cylinder")
  {
    for(size_t i = 0; i < size(); ++i)
    {
      std::pair<long,long> mm = this->at(i);
      if (mm.first<0) {
        g_log.error("Major multipole moment must be >= 0!");
        g_log.error("major="+boost::lexical_cast<std::string>(mm.first)+", minor="+boost::lexical_cast<std::string>(mm.second));
        throw;
      }
      if ((mm.second<0) || (mm.second>3) )
      {
        g_log.error("Minor multipole moment must be between 0 and 3!");
        g_log.error("major="+boost::lexical_cast<std::string>(mm.first)+", minor="+boost::lexical_cast<std::string>(mm.second));
        throw;
      }
      if ((mm.first==0) && (mm.second!=0) )
      {
        g_log.error("Minor multipole moment must be 0 for Major 0!");
        g_log.error("major="+boost::lexical_cast<std::string>(mm.first)+", minor="+boost::lexical_cast<std::string>(mm.second));
        throw;
      }
    }
  }
  else if(Params::Inst()->scattering.average.orientation.multipole.type=="sphere")
  {
    for(size_t i = 0; i < size(); ++i)
    {
      std::pair<long,long> mm = this->at(i);
      if (mm.first<0)
      {
        g_log.error("Major multipole moment must be >= 0!");
        g_log.error("major="+boost::lexical_cast<std::string>(mm.first));
        throw;
      }
      if (labs(mm.second)>mm.first)
      {
        g_log.error("Minor multipole moment must be between -Major and +Major!");
        g_log.error("major="+boost::lexical_cast<std::string>(mm.first)+", minor="+boost::lexical_cast<std::string>(mm.second));
        throw;
      }
    }
  }
  g_log.information("Initialized orientational averaging with "+boost::lexical_cast<std::string>(size())+" moments.");
} // end of ScatteringAverageOrientationMultipoleMomentsParameters::create()

/// read out the scans and push the results onto the internal vector
void ScatteringVectorsParameters::create_from_scans()
{
  if (scans.size()>3)
  {
    g_log.error("More than 3 scan definitions are not supported.");
    throw;
  }
  // local vector unfolds each scan first, then we'll do element-wise vector addition
  vector< vector<CartesianCoor3D> > qvectors(scans.size());
  for(size_t i = 0; i < scans.size(); ++i)
  {
    if (scans[i].points==0) continue;
    if (scans[i].points==1)
    {
      double scal = (scans[i].from+scans[i].to)/2;
      qvectors[i].push_back(scal*scans[i].basevector);
      continue;
    }
    if (scans[i].points==2)
    {
      qvectors[i].push_back(scans[i].from*scans[i].basevector);
      qvectors[i].push_back(scans[i].to*scans[i].basevector);
      continue;
    }
    // else
    qvectors[i].push_back(scans[i].from*scans[i].basevector);
    for (size_t j=1;j<(scans[i].points-1);j++)
    {
      double base = static_cast<double>(j)*1.0/(static_cast<double>(scans[i].points)-1.0);
      double scal = scans[i].from + pow(base,scans[i].exponent)*(scans[i].to-scans[i].from);
      qvectors[i].push_back(scal*scans[i].basevector);
    }
    qvectors[i].push_back(scans[i].to*scans[i].basevector);
  }
  // trivial case: only one scan!
  if (scans.size()==1)
  {
    for(size_t i = 0; i < qvectors[0].size(); ++i)
    {
      this->push_back(qvectors[0][i]);
    }
    return;
  }
  // if scans.size()>1 , not trivial..
  if (scans.size()==2)
  {
    for(size_t i = 0; i < qvectors[0].size(); ++i)
    {
      for(size_t j = 0; j < qvectors[1].size(); ++j)
      {
        this->push_back(qvectors[0][i]+qvectors[1][j]);
      }
    }
  }
  if (scans.size()==3)
  {
    for(size_t i = 0; i < qvectors[0].size(); ++i)
    {
      for(size_t j = 0; j < qvectors[1].size(); ++j)
      {
        for(size_t k = 0; k < qvectors[2].size(); ++k)
        {
          this->push_back(qvectors[0][i]+qvectors[1][j]+qvectors[2][k]);
        }
      }
    }
  }
} // end of ScatteringVectorsParameters::create_from_scans()

std::string Params::get_filepath(std::string fname)
{
  using namespace boost::filesystem;
  path fpath(fname);
  std::string fdir;
  if(fpath.parent_path().is_complete())
    fdir = fpath.parent_path().string();
  else if (!config_rootpath.empty())
    fdir = (path(config_rootpath) / fpath.parent_path()).string();
  else
    fdir = (initial_path() / fpath.parent_path()).string();
  return (path(fdir) / fpath.filename()).string();
}

void Params::write_xml_to_file(std::string filename)
{
  std::ofstream conf(filename.c_str());
  conf << write_xml();
  conf.close();
}

/// store the configuration line by line into an internal buffer
void Params::read_xml(std::string filename)
{
  // this is for keeping history
  std::ifstream conffile(filename.c_str());
  char c;
  while (conffile.good())
  {
    conffile.get(c);
    if (conffile.good())
    {
      rawconfig.push_back(c);
    }
  }
  conffile.close();
  // START OF sample section //
  XMLInterface xmli(filename);
  g_log.information("Reading parameters from XML file: "+filename);
  xmli.dump(config);
  // now read the parameters
  sample.structure.file="sample.pdb";
  sample.structure.filepath = get_filepath(sample.structure.file);
  sample.structure.format="pdb";
  if (xmli.exists("//sample"))
  {
    if (xmli.exists("//sample/structure"))
    {
      if (xmli.exists("//sample/structure/file"))
      {
        sample.structure.file = get_filepath(xmli.get_value<std::string>("//sample/structure/file"));
        sample.structure.filepath = get_filepath(sample.structure.file);
      }
      if (xmli.exists("//sample/structure/format"))
      {
        sample.structure.format   = xmli.get_value<std::string>("//sample/structure/format");
      }
    }
    if (xmli.exists("//sample/selections"))
    {
      std::vector<XMLElement> selections = xmli.get("//sample/selections/selection");
      for(size_t i = 0; i < selections.size(); ++i)
      {
        xmli.set_current(selections[i]);
        std::string type = "index";
        if (xmli.exists("./type")) type = xmli.get_value<std::string>("./type");
        if (type=="index")
        {
          std::string sname = "_"+boost::lexical_cast<std::string>(i);
          std::vector<size_t> ids;
          if (xmli.exists("./name")) sname = xmli.get_value<std::string>("./name") ;
          std::vector<XMLElement> indexes = xmli.get("./index");
          for(size_t pi = 0; pi < indexes.size(); ++pi)
          {
            xmli.set_current(indexes[pi]);
            size_t index = xmli.get_value<size_t>(".");
            ids.push_back(index);
          }
          sample.selections[sname] = new SampleIndexSelectionParameters(ids);
          g_log.information("Creating Index Atomselection: "+sname+" , elements:"+boost::lexical_cast<std::string>(ids.size()));
        }
        else if (type=="range")
        {
          std::string sname = "_"+boost::lexical_cast<std::string>(i);
          size_t from = 0;
          size_t to = 0;
          if (xmli.exists("./name")) sname = xmli.get_value<std::string>("./name");
          if (xmli.exists("./from")) from = xmli.get_value<size_t>("./from");
          if (xmli.exists("./to")) to = xmli.get_value<size_t>("./to") ;
          sample.selections[sname] = new SampleRangeSelectionParameters(from,to);
          g_log.information("Creating Range Atomselection: " + sname +
            " , from:"+boost::lexical_cast<std::string>(from) +
            ", to: "+boost::lexical_cast<std::string>(to) );
        }
        else if (type=="lexical")
        {
          std::string sname = "_"+boost::lexical_cast<std::string>(i);
          std::string expression = "";
          if (xmli.exists("./name")) sname = xmli.get_value<std::string>("./name");
          if (xmli.exists("./expression")) expression = xmli.get_value<std::string>("./expression");
          sample.selections[sname] = new SampleLexicalSelectionParameters(expression);
          g_log.information("Creating Lexical Atomselection: "+sname+" , expression:"+expression);
        }
        else if (type=="file")
        {
          std::string sname = "_"+boost::lexical_cast<std::string>(i); // NOT used by ndx file format
          std::string file  = "selection.pdb";
          std::string filepath = get_filepath(file);
          std::string format = "pdb";
          std::string selector = "beta";
          std::string expression = "1|1\\.0|1\\.00";
          if (xmli.exists("./name")) sname = xmli.get_value<std::string>("./name");
          if (xmli.exists("./format")) format = xmli.get_value<std::string>("./format");
          // this is a convenience overwrite for the ndx file format
          if (format=="ndx")
          {
            filename = "index.ndx";
            selector = "name";
            expression = ".*";
          }
          if (xmli.exists("./file"))
          {
            filename = xmli.get_value<std::string>("./file");
            filepath = get_filepath(file);
          }
          if (xmli.exists("./selector")) selector = xmli.get_value<std::string>("./selector");
          if (xmli.exists("./expression")) expression = xmli.get_value<std::string>("./expression");
          sample.selections[sname] = new SampleFileSelectionParameters(filepath,format,selector,expression);
          g_log.information("Creating File Atomselection: "+sname+" , file: "+filepath+
            ", format: "+format+", selector: "+selector+", expression: "+expression);
        }
      }
    }
    if (xmli.exists("//sample/framesets"))
    {
      size_t def_first=0;
      size_t def_last=0;
      bool def_last_set=false;
      size_t def_stride = 1;
      size_t def_clones = 1;
      std::string def_format = "dcd";
      std::string def_file="sample.dcd";
      if (xmli.exists("//sample/framesets/first")) def_first = xmli.get_value<size_t>("//sample/framesets/first");
      if (xmli.exists("//sample/framesets/last"))
      {
        def_last = xmli.get_value<size_t>("//sample/framesets/last");
        def_last_set = true;
      }
      if (xmli.exists("//sample/framesets/stride")) def_stride = xmli.get_value<size_t>("//sample/framesets/stride");
      if (xmli.exists("//sample/framesets/clones")) def_clones = xmli.get_value<size_t>("//sample/framesets/clones");
      if (xmli.exists("//sample/framesets/format")) def_first = xmli.get_value<size_t>("//sample/framesets/format");
      // read in frame information
      std::vector<XMLElement> framesets = xmli.get("//sample/framesets/frameset");
      for(size_t i = 0; i < framesets.size(); ++i)
      {
        xmli.set_current(framesets[i]);
        SampleFramesetParameters fset;
        fset.first = def_first;
        fset.last = def_last;
        fset.last_set = def_last_set;
        fset.stride = def_stride;
        fset.format = def_format;
        fset.clones = def_clones;
        fset.file = def_file;
        fset.filepath = get_filepath(fset.file);
        if (xmli.exists("./file"))
        {
          fset.file = xmli.get_value<std::string>("./file");
          fset.filepath = get_filepath(fset.file);
        }
        boost::filesystem::path index_path = get_filepath(xmli.get_value<std::string>("./file"));
        fset.index = index_path.parent_path().string()+"/"+index_path.stem().string()+".tnx";
        fset.index_default = true;
        if (xmli.exists("./format")) fset.format = xmli.get_value<std::string>("./format");
        if (xmli.exists("./first")) fset.first = xmli.get_value<size_t>("./first");
        if (xmli.exists("./last"))
        {
          fset.last = xmli.get_value<size_t>("./last");
          fset.last_set = true;
        }
        if (xmli.exists("./stride"))  fset.stride = xmli.get_value<size_t>("./stride");
        if (xmli.exists("./clones"))  fset.clones = xmli.get_value<size_t>("./clones");
        if (xmli.exists("./index"))
        {
          fset.index = get_filepath(xmli.get_value<std::string>("./index"));
          fset.indexpath = get_filepath(fset.index);
          fset.index_default = false;
        }
        sample.framesets.push_back(fset);
        g_log.information("Added frames from "+fset.filepath+" using format: "+fset.format);
        g_log.information("Options: first="+boost::lexical_cast<std::string>(fset.first)+
          ", last="+boost::lexical_cast<std::string>(fset.last)+
          ", lastset="+boost::lexical_cast<std::string>(fset.last_set)+
          ", stride="+boost::lexical_cast<std::string>(fset.stride));
      }
    }
    if (xmli.exists("//sample/motions"))
    {
      std::vector<XMLElement> motions = xmli.get("//sample/motions/motion");
      for(size_t i = 0; i < motions.size(); ++i)
      {
        xmli.set_current(motions[i]);
        SampleMotionParameters motion;