diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index ad21caa47d9a81b6c5fd73e44796fcb97ff38c0d..1ff4bdddd72281e964bafb0b1490ca6aca1b2bdf 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -90,7 +90,7 @@ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const s
 
     // we assume variable bin size and prepare [nbins+1] array of left edges of each bin plus right edge of the last bin
     for(size_t i_axis=0; i_axis<output.getNdimensions(); ++i_axis) {
-        const AxisDouble *axis = output.getAxis(i_axis);
+        const IAxis *axis = output.getAxis(i_axis);
         if( !axis ) throw("IsGISAXSTools::getOutputDataTH123D() -> Error! Can't cast axis");
         double dx(0);
         haxises[i_axis].nbins = axis->getSize();
@@ -154,7 +154,7 @@ TH1 *IsGISAXSTools::getOutputDataTH123D(const OutputData<double>& output, const
 
     // we assume variable bin size and prepare [nbins+1] array of left edges of each bin plus right edge of the last bin
     for(size_t i_axis=0; i_axis<output.getNdimensions(); ++i_axis) {
-        const AxisDouble *axis = output.getAxis(i_axis);
+        const IAxis *axis = output.getAxis(i_axis);
         if( !axis ) throw("IsGISAXSTools::getOutputDataTH123D() -> Error! Can't cast axis");
         double dx(0);
         haxises[i_axis].nbins = axis->getSize();
@@ -465,8 +465,8 @@ void IsGISAXSTools::exportOutputDataInVectors2D(const OutputData<double> &output
 {
     if (output_data.getRank() != 2) return;
 
-    const AxisDouble *p_axis0 = output_data.getAxis(0);
-    const AxisDouble *p_axis1 = output_data.getAxis(1);
+    const IAxis *p_axis0 = output_data.getAxis(0);
+    const IAxis *p_axis1 = output_data.getAxis(1);
     std::string axis0_name = p_axis0->getName();
     std::string axis1_name = p_axis1->getName();
     size_t axis0_size = p_axis0->getSize();
diff --git a/App/src/TestFormFactor.cpp b/App/src/TestFormFactor.cpp
index 45df3456288963c44fe0d7b24c97814b9e64126f..acd7f3844d91fa8e9a9445f4a45762c1ee8412b3 100644
--- a/App/src/TestFormFactor.cpp
+++ b/App/src/TestFormFactor.cpp
@@ -27,8 +27,8 @@ TestFormFactor::~TestFormFactor()
 void TestFormFactor::execute()
 {
     OutputData<double>::iterator it = mp_intensity_output->begin();
-    const AxisDouble *p_y_axis = mp_intensity_output->getAxis("detector y-axis");
-    const AxisDouble *p_z_axis = mp_intensity_output->getAxis("detector z-axis");
+    const IAxis *p_y_axis = mp_intensity_output->getAxis("detector y-axis");
+    const IAxis *p_z_axis = mp_intensity_output->getAxis("detector z-axis");
     double lambda = 1.0;
     double alpha_i = 0.2*M_PI/180.0;
     cvector_t k_i;
@@ -53,8 +53,8 @@ void TestFormFactor::draw()
     TCanvas *c1 = new TCanvas("c1_test_formfactor", "Cylinder Formfactor", 0, 0, 1024, 768);
     (void)c1;
 
-    const AxisDouble *p_y_axis = mp_intensity_output->getAxis("detector y-axis");
-    const AxisDouble *p_z_axis = mp_intensity_output->getAxis("detector z-axis");
+    const IAxis *p_y_axis = mp_intensity_output->getAxis("detector y-axis");
+    const IAxis *p_z_axis = mp_intensity_output->getAxis("detector z-axis");
     size_t y_size = p_y_axis->getSize();
     size_t z_size = p_z_axis->getSize();
     double y_start = (*p_y_axis)[0];
diff --git a/App/src/TestFresnelCoeff.cpp b/App/src/TestFresnelCoeff.cpp
index 2476675cb306a0f28bf5c2f08a026f1b9508d3a9..671108d6276cfd03db2e05aa99d5b8a77b0f0d4f 100644
--- a/App/src/TestFresnelCoeff.cpp
+++ b/App/src/TestFresnelCoeff.cpp
@@ -268,7 +268,7 @@ void TestFresnelCoeff::draw_roughness_set()
     size_t ncoeffs = 2;
     enum key_coeffs { kCoeffR, kCoeffT};
 
-    const AxisDouble *p_rough = mp_coeffs->getAxis("roughness");
+    const IAxis *p_rough = mp_coeffs->getAxis("roughness");
     size_t nroughness = p_rough->getSize();
 
 
diff --git a/App/src/TestIsGISAXS12.cpp b/App/src/TestIsGISAXS12.cpp
index fdac6ae83ecba87b9a8637f80a29e229ae20dc3b..524cab456813cf3ca5f586103565f60910b2e6cc 100644
--- a/App/src/TestIsGISAXS12.cpp
+++ b/App/src/TestIsGISAXS12.cpp
@@ -235,7 +235,7 @@ void TestIsGISAXS12::run_isgisaxs_fit()
         OutputData<double > *real = obj->getChiSquaredModule()->getRealData()->clone();
         OutputData<double > *simul = obj->getChiSquaredModule()->getSimulationData()->clone();
 
-        c2->cd(i_set+3);
+        c2->cd((int)(i_set+3));
         *simul /= *real;
         TH1D *hratio = IsGISAXSTools::getOutputDataScanHist(*simul,"gisasfw_real_simul_ratio");
         hratio->DrawCopy();
@@ -766,7 +766,7 @@ void TestIsGISAXS12::print_axes(DataSet &data)
     for(size_t i_set=0; i_set<data.size(); ++i_set) {
         std::cout << "scan #" << i_set << "  ";
         for(size_t i_axis=0; i_axis<data[i_set]->getNdimensions(); ++i_axis) {
-            const AxisDouble *axis = data[i_set]->getAxis(i_axis);
+            const IAxis *axis = data[i_set]->getAxis(i_axis);
             std::cout << "( " << axis->getName() << ", " << axis->getSize() << ", " << axis->getMin() << ", " << axis->getMax() << " )   ";
         }
         std::cout << std::endl;
diff --git a/App/src/TestMesoCrystal2.cpp b/App/src/TestMesoCrystal2.cpp
index 083b28231f99e8edd37f36b2417d12d08c1197c6..c6e2a5e20d63cec5633d430ab81214c398da09e4 100644
--- a/App/src/TestMesoCrystal2.cpp
+++ b/App/src/TestMesoCrystal2.cpp
@@ -226,8 +226,8 @@ void TestMesoCrystal2::initializeExperiment(const OutputData<double> *output_dat
         mp_experiment->setDetectorParameters(200, 0.3*Units::degree, 0.073, 200, -0.4*Units::degree, 0.066);
     } else {
         // if there is output_data as input parameter, build detector using output_data axises
-        const AxisDouble *axis0 = output_data->getAxis(0);
-        const AxisDouble *axis1 = output_data->getAxis(1);
+        const IAxis *axis0 = output_data->getAxis(0);
+        const IAxis *axis1 = output_data->getAxis(1);
         //std::cout << axis0->getSize() << " " << axis0->getMin() << " " << axis0->getMax() << " " << axis1->getSize() << " " << axis1->getMin() << " " << axis1->getMax() << std::endl;
         mp_experiment->setDetectorParameters(axis0->getSize(), axis0->getMin(), axis0->getMax(), axis1->getSize(), axis1->getMin(), axis1->getMax());
     }
diff --git a/App/src/TestMiscellaneous.cpp b/App/src/TestMiscellaneous.cpp
index c2ddc317ccc0612eeaf31f0d2c843858c1afd3bc..1add96973a7b1933b19e0cb9d0046ac9b077eef7 100644
--- a/App/src/TestMiscellaneous.cpp
+++ b/App/src/TestMiscellaneous.cpp
@@ -300,7 +300,7 @@ void TestMiscellaneous::test_DoubleToComplexInterpolatingFunction()
 
     OpticalFresnel fresnelCalculator;
 
-    const AxisDouble *p_alpha_axis = data_alpha->getAxis("alpha_f");
+    const IAxis *p_alpha_axis = data_alpha->getAxis("alpha_f");
     std::map<double, OpticalFresnel::MultiLayerCoeff_t> fresnel_coeff_map;
     for (size_t i=0; i<p_alpha_axis->getSize(); ++i) {
         double angle = (*p_alpha_axis)[i];
diff --git a/App/src/TestRootTree.cpp b/App/src/TestRootTree.cpp
index a2e50f9aea6a01e5b730d536a383673d9a061f62..ffbd26dfa4e67d8361b0b9cac58aefbf752f740c 100644
--- a/App/src/TestRootTree.cpp
+++ b/App/src/TestRootTree.cpp
@@ -259,8 +259,8 @@ void TestRootTree::simple_write()
 
         mp_data = mp_experiment->getOutputDataClone();
         // accessing to scattering data
-        const AxisDouble *axis0 = mp_data->getAxis(0);
-        const AxisDouble *axis1 = mp_data->getAxis(1);
+        const IAxis *axis0 = mp_data->getAxis(0);
+        const IAxis *axis1 = mp_data->getAxis(1);
         std::string axis0_name = axis0->getName();
         std::string axis1_name = axis1->getName();
 
diff --git a/Core/Algorithms/inc/Detector.h b/Core/Algorithms/inc/Detector.h
index 26ad22825b35bf4064d4dd8d6d5395d95149be79..8b6e371d956fc16743a47cfb2445ef17d691f4bf 100644
--- a/Core/Algorithms/inc/Detector.h
+++ b/Core/Algorithms/inc/Detector.h
@@ -16,6 +16,7 @@
 
 #include "IDetectorResolution.h"
 #include "IParameterized.h"
+#include "SafePointerVector.h"
 
 
 #include <vector>
@@ -33,8 +34,8 @@ public:
 
 	virtual ~Detector();
 
-	void addAxis(const AxisDouble &axis);
-	AxisDouble getAxis(size_t index) const;
+	void addAxis(const IAxis &axis);
+	const IAxis &getAxis(size_t index) const;
 	size_t getDimension() const { return m_axes.size(); }
 	void clear();
     void setDetectorResolution(IDetectorResolution *p_detector_resolution) { delete mp_detector_resolution; mp_detector_resolution = p_detector_resolution; }
@@ -52,7 +53,7 @@ private:
     //! swap function
     void swapContent(Detector &other);
 
-	std::vector<AxisDouble > m_axes;
+    SafePointerVector<IAxis> m_axes;
 	IDetectorResolution *mp_detector_resolution;
 };
 
diff --git a/Core/Algorithms/inc/GISASExperiment.h b/Core/Algorithms/inc/GISASExperiment.h
index ee2364bc3afe354eb75f8246b90c50e962ec75d8..03458059961657a51736da27f141d031b2b87f9b 100644
--- a/Core/Algorithms/inc/GISASExperiment.h
+++ b/Core/Algorithms/inc/GISASExperiment.h
@@ -58,7 +58,7 @@ private:
 	double deltaPhi(double alpha, double phi, double zeta);
 	void createZetaAndProbVectors(std::vector<double> &zetas, std::vector<double> &probs, size_t nbr_zetas, double zeta_sigma);
 	void addToIntensityMap(double alpha, double phi, double value);
-	int findClosestIndex(const AxisDouble *p_axis, double value);
+//	int findClosestIndex(const AxisDouble *p_axis, double value);
 };
 
 #endif /* GISASEXPERIMENT_H_ */
diff --git a/Core/Algorithms/src/ConvolutionDetectorResolution.cpp b/Core/Algorithms/src/ConvolutionDetectorResolution.cpp
index a77ee5a546764083e17759b78477ce1bf2cb7655..ecaac74237405286f944caca8e2cce51423c1659 100644
--- a/Core/Algorithms/src/ConvolutionDetectorResolution.cpp
+++ b/Core/Algorithms/src/ConvolutionDetectorResolution.cpp
@@ -98,7 +98,7 @@ void ConvolutionDetectorResolution::apply1dConvolution(OutputData<double>* p_int
     if (p_intensity_map->getRank() != 1) {
         throw LogicErrorException("ConvolutionDetectorResolution::apply1dConvolution() -> Error! Number of axes for intensity map does not correspond to the dimension of the map.");
     }
-    const AxisDouble *p_axis = p_intensity_map->getAxis(0);
+    const IAxis *p_axis = p_intensity_map->getAxis(0);
     // Construct source vector from original intensity map
     std::vector<double> source_vector = p_intensity_map->getRawDataVector();
     size_t data_size = source_vector.size();
@@ -131,8 +131,8 @@ void ConvolutionDetectorResolution::apply2dConvolution(OutputData<double>* p_int
     if (p_intensity_map->getRank() != 2) {
         throw LogicErrorException("ConvolutionDetectorResolution::apply2dConvolution() -> Error! Number of axes for intensity map does not correspond to the dimension of the map.");
     }
-    const AxisDouble *p_axis_1 = p_intensity_map->getAxis(0);
-    const AxisDouble *p_axis_2 = p_intensity_map->getAxis(1);
+    const IAxis *p_axis_1 = p_intensity_map->getAxis(0);
+    const IAxis *p_axis_2 = p_intensity_map->getAxis(1);
     size_t axis_size_1 = p_axis_1->getSize();
     size_t axis_size_2 = p_axis_2->getSize();
     if (axis_size_1 < 2 || axis_size_2 < 2) {
diff --git a/Core/Algorithms/src/Detector.cpp b/Core/Algorithms/src/Detector.cpp
index 133d45c58a56bf4a4e1c836127757e0c2d45f003..61233d92711e2da1473bc036d8d54b5d41ec97bd 100644
--- a/Core/Algorithms/src/Detector.cpp
+++ b/Core/Algorithms/src/Detector.cpp
@@ -47,15 +47,15 @@ void Detector::swapContent(Detector &other)
 /* ************************************************************************* */
 // other methods
 /* ************************************************************************* */
-void Detector::addAxis(const AxisDouble &axis)
+void Detector::addAxis(const IAxis &axis)
 {
-	m_axes.push_back(axis);
+	m_axes.push_back(axis.clone());
 }
 
-AxisDouble Detector::getAxis(size_t index) const
+const IAxis &Detector::getAxis(size_t index) const
 {
 	if (isCorrectAxisIndex(index)) {
-		return m_axes[index];
+		return *m_axes[index];
 	}
 	throw OutOfBoundsException("Not so many axes in this detector.");
 }
diff --git a/Core/Algorithms/src/Experiment.cpp b/Core/Algorithms/src/Experiment.cpp
index 9c3e22f8cb75c75e53424a32f406eb21177cf37b..5aeeebb8aa1b36c69a0f6f3762ff702256bd9ea7 100644
--- a/Core/Algorithms/src/Experiment.cpp
+++ b/Core/Algorithms/src/Experiment.cpp
@@ -184,7 +184,7 @@ void Experiment::setDetectorParameters(const OutputData<double > &output_data)
     //std::cout << "Experiment::setDetectorParameters() -> Info. Adjusting detector to have shape as in given output data" << std::endl;
     m_detector.clear();
     for(size_t i_axis=0; i_axis<output_data.getNdimensions(); ++i_axis) {
-        const AxisDouble *axis = output_data.getAxis(i_axis);
+        const IAxis *axis = output_data.getAxis(i_axis);
         m_detector.addAxis(*axis);
     }
     updateIntensityMapAxes();
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index f400c35ad38ce6267668f00a8071e6fc8f36a906..2b20a564dc73ca7b009c978295feeb4f006e2138 100644
--- a/Core/Algorithms/src/GISASExperiment.cpp
+++ b/Core/Algorithms/src/GISASExperiment.cpp
@@ -194,8 +194,8 @@ double GISASExperiment::getSolidAngle(size_t index) const
     const std::string s_alpha_f("alpha_f");
     const std::string s_phi_f("phi_f");
 
-    const AxisDouble *p_alpha_axis = m_intensity_map.getAxis(s_alpha_f);
-    const AxisDouble *p_phi_axis = m_intensity_map.getAxis(s_phi_f);
+    const IAxis *p_alpha_axis = m_intensity_map.getAxis(s_alpha_f);
+    const IAxis *p_phi_axis = m_intensity_map.getAxis(s_phi_f);
     size_t alpha_index = m_intensity_map.getIndexOfAxis(s_alpha_f, index);
     size_t alpha_size = p_alpha_axis->getSize();
     size_t phi_index = m_intensity_map.getIndexOfAxis(s_phi_f, index);
@@ -283,23 +283,23 @@ void GISASExperiment::createZetaAndProbVectors(std::vector<double>& zetas,
 
 void GISASExperiment::addToIntensityMap(double alpha, double phi, double value)
 {
-    const AxisDouble *p_alpha_axis = m_intensity_map.getAxis("alpha_f");
-    const AxisDouble *p_phi_axis = m_intensity_map.getAxis("phi_f");
+    const IAxis *p_alpha_axis = m_intensity_map.getAxis("alpha_f");
+    const IAxis *p_phi_axis = m_intensity_map.getAxis("phi_f");
     std::vector<int> coordinates;
-    coordinates.push_back(findClosestIndex(p_alpha_axis, alpha));
-    coordinates.push_back(findClosestIndex(p_phi_axis, phi));
+    coordinates.push_back((int)p_alpha_axis->findClosestIndex(alpha));
+    coordinates.push_back((int)p_phi_axis->findClosestIndex(phi));
     m_intensity_map[m_intensity_map.toIndex(coordinates)] += value;
 }
 
-int GISASExperiment::findClosestIndex(const AxisDouble *p_axis, double value)
-{
-    int result = 0;
-    double smallest_diff = std::abs(value-(*p_axis)[0]);
-    for (size_t i=1; i<p_axis->getSize(); ++i) {
-        double new_diff = std::abs(value-(*p_axis)[i]);
-        if (new_diff > smallest_diff) break;
-        result = (int)i;
-        smallest_diff = new_diff;
-    }
-    return result;
-}
+//int GISASExperiment::findClosestIndex(const AxisDouble *p_axis, double value)
+//{
+//    int result = 0;
+//    double smallest_diff = std::abs(value-(*p_axis)[0]);
+//    for (size_t i=1; i<p_axis->getSize(); ++i) {
+//        double new_diff = std::abs(value-(*p_axis)[i]);
+//        if (new_diff > smallest_diff) break;
+//        result = (int)i;
+//        smallest_diff = new_diff;
+//    }
+//    return result;
+//}
diff --git a/Core/Algorithms/src/MultiLayerDWBASimulation.cpp b/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
index a2de088391ad2588436d74c6ad0c321188afddad..5614ba974b293335c0b8219db4783392ffa4016b 100644
--- a/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
+++ b/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
@@ -60,7 +60,7 @@ void MultiLayerDWBASimulation::run()
     kvector_t m_ki_real(m_ki.x().real(), m_ki.y().real(), m_ki.z().real());
 
     m_dwba_intensity.setAllTo(0.0);
-    const AxisDouble *p_alpha_axis = getDWBAIntensity().getAxis("alpha_f");
+    const IAxis *p_alpha_axis = getDWBAIntensity().getAxis("alpha_f");
     double lambda = 2.0*M_PI/m_ki_real.mag();
 
     typedef std::pair<double, OpticalFresnel::MultiLayerCoeff_t > doublefresnel_t;
diff --git a/Core/Core.pro b/Core/Core.pro
index 90c504a8c0c605e0a4b5fa462f44aa955f67ccdc..9a5f7b5a2bf95698b1aa8286b6900b647babae63 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -88,8 +88,8 @@ SOURCES += \
     Samples/src/ParticleDecoration.cpp \
     Samples/src/ParticleInfo.cpp \
     \
-    Tools/src/Axisdouble.cpp \
-    Tools/src/BinAxis.cpp \
+    Tools/src/AxisBin.cpp \
+    Tools/src/AxisDouble.cpp \
     Tools/src/Convolve.cpp \
     Tools/src/CoreOptionsDescription.cpp \
     Tools/src/DoubleToComplexInterpolatingFunction.cpp \
@@ -224,9 +224,10 @@ HEADERS += \
     Samples/inc/ParticleDecoration.h \
     Samples/inc/ParticleInfo.h \
     \
+    Tools/inc/AxisBin.h \
+    Tools/inc/AxisDouble.h \
     Tools/inc/AttLimits.h \
-    Tools/inc/Axisdouble.h \
-    Tools/inc/BinAxis.h \
+    Tools/inc/Bin.h \
     Tools/inc/Convolve.h \
     Tools/inc/Coordinate3D.h \
     Tools/inc/DoubleToComplexInterpolatingFunction.h \
@@ -239,6 +240,7 @@ HEADERS += \
     Tools/inc/FitSuiteObjects.h \
     Tools/inc/FitSuiteParameters.h \
     Tools/inc/FitSuiteStrategy.h \
+    Tools/inc/IAxis.h \
     Tools/inc/IDoubleToComplexFunction.h \
     Tools/inc/IFactory.h \
     Tools/inc/IMinimizer.h \
diff --git a/Core/Tools/inc/AxisBin.h b/Core/Tools/inc/AxisBin.h
new file mode 100644
index 0000000000000000000000000000000000000000..8e2f0e2614a78bdb0bd4afd8d2a018b95b036fed
--- /dev/null
+++ b/Core/Tools/inc/AxisBin.h
@@ -0,0 +1,69 @@
+#ifndef AXISBIN_H_
+#define AXISBIN_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   BinAxis.h
+//! @brief  Definition of BinAxis class
+//! @author Scientific Computing Group at FRM II
+//! @date   Dec 6, 2012
+
+#include "IAxis.h"
+
+#include <string>
+#include <vector>
+
+//- -------------------------------------------------------------------
+//! @class BinAxis
+//! @brief Definition of BinAxis class that stores the bins of an axis
+//- -------------------------------------------------------------------
+class AxisBin : public IAxis
+{
+public:
+    //! constructors
+    AxisBin(std::string name);
+    AxisBin(std::string name, size_t nbr_bins, double start, double end);
+
+    //! clone function
+    virtual AxisBin *clone() const;
+
+    //! create a new axis with half the number of bins
+    virtual AxisBin *createDoubleBinSize() const;
+
+    //! destructor
+    virtual ~AxisBin() {}
+
+    //! add new bin limit to the end
+    void push_back(double limit) { m_value_vector.push_back(limit); }
+
+    virtual size_t getSize() const;
+
+   //! indexed accessor retrieves midpoint of given bin
+    virtual double operator[](size_t index) const;
+
+    virtual Bin1D getBin(size_t index) const;
+
+    virtual double getMin() const { return m_value_vector.front(); }
+
+    virtual double getMax() const { return m_value_vector.back(); }
+
+    //! initialize axis bins
+    void initBins(size_t nbr_bins, double start, double end);
+
+    virtual size_t findClosestIndex(double value) const;
+
+//    //! find the bin that contains the given value
+//    Bin1D findMatchingBin(double value) const;
+protected:
+    virtual bool equals(const IAxis &other) const;
+private:
+    std::vector<double> m_value_vector;  //!< vector containing the bin limits
+};
+
+#endif /* AXISBIN_H_ */
diff --git a/Core/Tools/inc/AxisDouble.h b/Core/Tools/inc/AxisDouble.h
index def10467db7e5c68f5067458cc83dad8f8639afb..6f49aadfe9b2f7e4b9c56980840dfda3644608e2 100644
--- a/Core/Tools/inc/AxisDouble.h
+++ b/Core/Tools/inc/AxisDouble.h
@@ -14,17 +14,19 @@
 //! @author Scientific Computing Group at FRM II
 //! @date   Dec 4, 2012
 
+#include "IAxis.h"
+
 #include <string>
 #include <vector>
 
 // Forward declaration of BinAxis class, as needed for conversion constructor
-class BinAxis;
+class AxisBin;
 
 //- -------------------------------------------------------------------
 //! @class AxisDouble
 //! @brief Definition of AxisDouble class that stores the points of an axis
 //- -------------------------------------------------------------------
-class AxisDouble
+class AxisDouble : public IAxis
 {
 public:
     //! constructors
@@ -33,46 +35,35 @@ public:
 
     //! explicit conversion from BinAxis
     //TODO: make explicit
-    AxisDouble(const BinAxis &source);
+    AxisDouble(const AxisBin &source);
 
-    //! clone function
-    AxisDouble* clone() const;
+    virtual AxisDouble *clone() const;
 
-    //! create a new axis with half the number of bins
-    AxisDouble createDoubleBinSize() const;
+    virtual AxisDouble *createDoubleBinSize() const;
 
     //! destructor
-    ~AxisDouble() {}
-
-    //! retrieve the number of bins
-    size_t getSize() const { return m_value_vector.size(); }
-
-    //! retrieve the label of the axis
-    std::string getName() const { return m_name; }
+    virtual ~AxisDouble() {}
 
-    //! set the axis label
-    void setName(std::string name) { m_name = name; }
+    virtual size_t getSize() const { return m_sample_vector.size(); }
 
     //! add point to the end of the axis
-    void push_back(double element) { m_value_vector.push_back(element); }
+    void push_back(double element) { m_sample_vector.push_back(element); }
 
-    //! indexed accessor
-    double& operator[](size_t index) { return m_value_vector.at(index); }
+    virtual double operator[](size_t index) const { return m_sample_vector[index]; }
 
-    //! const indexed accessor
-    const double& operator[](size_t index) const { return m_value_vector.at(index); }
+    virtual Bin1D getBin(size_t index) const;
 
     //! get value of first point of axis
-    double getMin() const { return m_value_vector.front(); }
+    virtual double getMin() const { return m_sample_vector.front(); }
 
     //! get value of last point of axis
-    double getMax() const { return m_value_vector.back(); }
+    virtual double getMax() const { return m_sample_vector.back(); }
 
     //! initialize axis points
     void initElements(size_t size, double start, double end);
 
     //! find number of bin which is closest to given value
-    size_t findClosestIndex(double value) const;
+    virtual size_t findClosestIndex(double value) const;
 
     //! find the index that corresponds to the given lower bound (index is inclusive)
     size_t getLowerBoundIndex(double value) const;
@@ -80,13 +71,9 @@ public:
     //! find the index that corresponds to the given upper bound (index is inclusive)
     size_t getUpperBoundIndex(double value) const;
 private:
-    std::string m_name;  //!< axis label
-    std::vector<double> m_value_vector;  //!< vector containing the axis points
+    virtual bool equals(const IAxis &other) const;
+    std::vector<double> m_sample_vector;  //!< vector containing the axis points
+    double m_bin_size;
 };
 
-//! global helper function for comparison of axes
-bool HaveSameNameAndShape(const AxisDouble &left, const AxisDouble &right);
-
-
-
 #endif /* AXISDOUBLE_H_ */
diff --git a/Core/Tools/inc/Bin.h b/Core/Tools/inc/Bin.h
new file mode 100644
index 0000000000000000000000000000000000000000..655d6c4bd1a1e39bb98b269306d592cc23cc07c2
--- /dev/null
+++ b/Core/Tools/inc/Bin.h
@@ -0,0 +1,46 @@
+#ifndef BIN_H_
+#define BIN_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   Bin.h
+//! @brief  Definition of Bin classes
+//! @author Scientific Computing Group at FRM II
+//! @date   Dec 10, 2012
+
+#include "Numeric.h"
+
+#include <cmath>
+
+//- -------------------------------------------------------------------
+//! @class Bin1D
+//! @brief Definition of Bin1D class that stores the bounds of a
+//! one-dimensional bin
+//- -------------------------------------------------------------------
+struct Bin1D
+{
+    double m_lower;  //!< lower bound of the bin
+    double m_upper;  //!< upper bound of the bin
+    double getMidPoint() const { return (m_lower + m_upper)/2.0; }
+};
+
+//! equality operator for bins
+inline bool operator==(const Bin1D &left, const Bin1D &right)
+{
+    if (std::abs(left.m_lower - right.m_lower) > Numeric::double_epsilon) return false;
+    if (std::abs(left.m_upper - right.m_upper) > Numeric::double_epsilon) return false;
+    return true;
+}
+
+//! inequality operator for bins
+inline bool operator!=(const Bin1D &left, const Bin1D &right) {
+    return !(left==right);
+}
+
+#endif /* BIN_H_ */
diff --git a/Core/Tools/inc/BinAxis.h b/Core/Tools/inc/BinAxis.h
deleted file mode 100644
index f9479160bc12452a3864bee56b990e4ec5b458ff..0000000000000000000000000000000000000000
--- a/Core/Tools/inc/BinAxis.h
+++ /dev/null
@@ -1,101 +0,0 @@
-#ifndef BINAXIS_H_
-#define BINAXIS_H_
-// ********************************************************************
-// * The BornAgain project                                            *
-// * Simulation of neutron and x-ray scattering at grazing incidence  *
-// *                                                                  *
-// * LICENSE AND DISCLAIMER                                           *
-// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
-// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
-// * mollis quis. Mauris commodo rhoncus porttitor.                   *
-// ********************************************************************
-//! @file   BinAxis.h
-//! @brief  Definition of BinAxis class
-//! @author Scientific Computing Group at FRM II
-//! @date   Dec 6, 2012
-
-#include <string>
-#include <vector>
-
-//- -------------------------------------------------------------------
-//! @class Bin
-//! @brief Definition of Bin class that stores the bounds of a bin
-//- -------------------------------------------------------------------
-struct Bin1D
-{
-    double m_lower;  //!< lower bound of the bin
-    double m_upper;  //!< upper bound of the bin
-    double getMidPoint() { return (m_lower + m_upper)/2.0; }
-};
-
-//! equality operator for bins
-bool operator==(const Bin1D &left, const Bin1D &right);
-
-//! inequality operator for bins
-bool operator!=(const Bin1D &left, const Bin1D &right) {
-    return !(left==right);
-}
-
-//- -------------------------------------------------------------------
-//! @class BinAxis
-//! @brief Definition of BinAxis class that stores the bins of an axis
-//- -------------------------------------------------------------------
-class BinAxis
-{
-public:
-    //! constructors
-    BinAxis(std::string name);
-    BinAxis(std::string name, size_t nbr_bins, double start, double end);
-
-    //! clone function
-    BinAxis* clone() const;
-
-    //! create a new axis with half the number of bins
-    BinAxis createDoubleBinSize() const;
-
-    //! destructor
-    ~BinAxis() {}
-
-    //! retrieve the number of bins
-    size_t getSize() const;
-
-    //! retrieve the label of the axis
-    std::string getName() const { return m_name; }
-
-    //! set the axis label
-    void setName(std::string name) { m_name = name; }
-
-    //! add new bin limit to the end
-    void push_back(double limit) { m_value_vector.push_back(limit); }
-
-    //! indexed accessor retrieves midpoint of given bin
-    Bin1D operator[](size_t index);
-
-    //! const indexed accessor retrieves midpoint of given bin
-    const Bin1D operator[](size_t index) const;
-
-    //! get minimum value of axis
-    double getMin() const { return m_value_vector.front(); }
-
-    //! get maximum value of axis
-    double getMax() const { return m_value_vector.back(); }
-
-    //! initialize axis bins
-    void initBins(size_t nbr_bins, double start, double end);
-
-    //! find number of bin for the given value
-    size_t findMatchingBinIndex(double value) const;
-
-    //! find the bin that contains the given value
-    Bin1D findMatchingBin(double value) const;
-
-private:
-    std::string m_name;  //!< axis label
-    std::vector<double> m_value_vector;  //!< vector containing the bin limits
-};
-
-//! global helper function for comparison of axes
-bool HaveSameNameAndShape(const BinAxis &left, const BinAxis &right);
-
-
-#endif /* BINAXIS_H_ */
diff --git a/Core/Tools/inc/IAxis.h b/Core/Tools/inc/IAxis.h
new file mode 100644
index 0000000000000000000000000000000000000000..ced92e9c13e7c5a32852a0c5bbd2e89ffa6ab37e
--- /dev/null
+++ b/Core/Tools/inc/IAxis.h
@@ -0,0 +1,96 @@
+#ifndef IAXIS_H_
+#define IAXIS_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   IAxis.h
+//! @brief  Definition of IAxis interface
+//! @author Scientific Computing Group at FRM II
+//! @date   Dec 10, 2012
+
+#include "Bin.h"
+
+#include <string>
+
+//- -------------------------------------------------------------------
+//! @class IAxis
+//! @brief Definition of IAxis interface for 1d axes
+//- -------------------------------------------------------------------
+
+class IAxis
+{
+public:
+    //! constructors
+    IAxis(std::string name) : m_name(name) {}
+
+    //! clone function
+    virtual IAxis *clone() const=0;
+
+    //! create a new axis with half the number of bins
+    virtual IAxis *createDoubleBinSize() const=0;
+
+    //! destructor
+    virtual ~IAxis() {}
+
+    //! retrieve the number of bins
+    virtual size_t getSize() const=0;
+
+    //! retrieve the label of the axis
+    std::string getName() const { return m_name; }
+
+    //! set the axis label
+    void setName(std::string name) { m_name = name; }
+
+    //! indexed accessor retrieves a sample
+    virtual double operator[](size_t index) const=0;
+
+    //! retrieve a 1d bin for the given index
+    virtual Bin1D getBin(size_t index) const=0;
+
+    //! get value of first point of axis
+    virtual double getMin() const=0;
+
+    //! get value of last point of axis
+    virtual double getMax() const=0;
+
+    //! find bin index which is best match for given value
+    virtual size_t findClosestIndex(double value) const=0;
+
+    //! test for equality
+    friend bool operator==(const IAxis &left, const IAxis &right) {
+        return left.equals(right);
+    }
+
+//    //! find the index that corresponds to the given lower bound (index is inclusive)
+//    virtual size_t getLowerBoundIndex(double value) const=0;
+//
+//    //! find the index that corresponds to the given upper bound (index is inclusive)
+//    virtual size_t getUpperBoundIndex(double value) const=0;
+protected:
+    virtual bool equals(const IAxis &other) const;
+    std::string m_name;  //!< axis label
+};
+
+inline bool IAxis::equals(const IAxis& other) const
+{
+    return getName()==other.getName();
+}
+
+//! test for inequality
+inline bool operator!=(const IAxis &left, const IAxis &right) {
+    return !(left ==  right);
+}
+
+//! global helper function for comparison of axes
+inline bool HaveSameNameAndShape(const IAxis &left, const IAxis &right)
+{
+    return left == right;
+}
+
+#endif /* IAXIS_H_ */
diff --git a/Core/Tools/inc/IMinimizer.h b/Core/Tools/inc/IMinimizer.h
index bdef292947e7f47b3d263b33b12f5d3ffcdcc727..886b5880745667412e31672998b662f00cf3e81e 100644
--- a/Core/Tools/inc/IMinimizer.h
+++ b/Core/Tools/inc/IMinimizer.h
@@ -104,7 +104,7 @@ public:
     //! return pointer to the parameters values at the minimum
     virtual double getValueOfVariableAtMinimum(size_t i) const
     {
-        std::map<int, double >::const_iterator pos = m_values.find(i);
+        std::map<int, double >::const_iterator pos = m_values.find((int)i);
         if(pos != m_values.end()){
             return pos->second;
         } else {
diff --git a/Core/Tools/inc/OutputData.h b/Core/Tools/inc/OutputData.h
index a2f4e236d8a25de531352612b95d59f5a7b2704c..f20129b7a623afabe627e972aac9a2f22481634a 100644
--- a/Core/Tools/inc/OutputData.h
+++ b/Core/Tools/inc/OutputData.h
@@ -19,6 +19,7 @@
 #include "Types.h"
 #include "LLData.h"
 #include "OutputDataIterator.h"
+#include "SafePointerVector.h"
 
 #include <string>
 #include <sstream>
@@ -38,11 +39,11 @@ public:
 
     void copyFrom(const OutputData<T> &x);
 
-    void addAxis(const AxisDouble &new_axis);
+    void addAxis(const IAxis &new_axis);
     void addAxis(std::string name, size_t size, double start, double end);
 
-    const AxisDouble *getAxis(size_t index) const;
-    const AxisDouble *getAxis(std::string label) const;
+    const IAxis *getAxis(size_t index) const;
+    const IAxis *getAxis(std::string label) const;
     size_t getAxisIndex(const std::string &label) const;
 
     // ---------------------------------
@@ -196,7 +197,7 @@ private:
     //! memory allocation for current dimensions configuration
     void allocate();
 
-    std::vector<AxisDouble> m_value_axes;
+    SafePointerVector<IAxis> m_value_axes;
     LLData<T> *mp_ll_data;
     Mask *mp_mask;
 };
@@ -247,12 +248,12 @@ template <class T> void OutputData<T>::copyFrom(const OutputData<T> &other)
 
 
 
-template <class T> void OutputData<T>::addAxis(const AxisDouble &new_axis)
+template <class T> void OutputData<T>::addAxis(const IAxis &new_axis)
 {
     if( getAxis(new_axis.getName()) ) throw LogicErrorException("OutputData<T>::addAxis(AxisDouble *) -> Error! Attempt to add axis with already existing name '"+ new_axis.getName()+std::string("'"));
     if (new_axis.getSize()>0)
     {
-        m_value_axes.push_back(new_axis);
+        m_value_axes.push_back(new_axis.clone());
         allocate();
     }
 }
@@ -265,18 +266,18 @@ void OutputData<T>::addAxis(std::string name, size_t size, double start, double
     addAxis(new_axis);
 }
 
-template <class T> const AxisDouble *OutputData<T>::getAxis(size_t index) const
+template <class T> const IAxis *OutputData<T>::getAxis(size_t index) const
 {
-    return &m_value_axes.at(index);
+    return m_value_axes[index];
 }
 
-template <class T> const AxisDouble *OutputData<T>::getAxis(std::string label) const
+template <class T> const IAxis *OutputData<T>::getAxis(std::string label) const
 {
-    for (std::vector<AxisDouble>::const_iterator it = m_value_axes.begin(); it != m_value_axes.end(); ++it)
+    for (size_t i = 0; i < m_value_axes.size(); ++i)
     {
-        if (it->getName() == label)
+        if (m_value_axes[i]->getName() == label)
         {
-            return &(*it);
+            return m_value_axes[i];
         }
     }
     return 0;
@@ -285,10 +286,9 @@ template <class T> const AxisDouble *OutputData<T>::getAxis(std::string label) c
 // return index of axis
 template <class T> size_t OutputData<T>::getAxisIndex(const std::string &label) const
 {
-    size_t index(0);
-    for (std::vector<AxisDouble>::const_iterator it = m_value_axes.begin(); it != m_value_axes.end(); ++it, ++index)
+    for (size_t i = 0; i < m_value_axes.size(); ++i)
     {
-        if (it->getName() == label) return index;
+        if (m_value_axes[i]->getName() == label) return i;
     }
     throw LogicErrorException("OutputData<T>::getIndexOfAxis() -> Error! Axis with given name not found '"+label+std::string("'"));
 }
@@ -374,8 +374,8 @@ template<class T> std::vector<int> OutputData<T>::toCoordinates(size_t index) co
     result.resize(mp_ll_data->getRank());
     for (size_t i=0; i<mp_ll_data->getRank(); ++i)
     {
-        result[mp_ll_data->getRank()-1-i] = (int)(remainder % m_value_axes[mp_ll_data->getRank()-1-i].getSize());
-        remainder /= m_value_axes[mp_ll_data->getRank()-1-i].getSize();
+        result[mp_ll_data->getRank()-1-i] = (int)(remainder % m_value_axes[mp_ll_data->getRank()-1-i]->getSize());
+        remainder /= m_value_axes[mp_ll_data->getRank()-1-i]->getSize();
     }
     return result;
 }
@@ -390,7 +390,7 @@ template <class T> size_t OutputData<T>::toIndex(std::vector<int> coordinates) c
     for (size_t i=mp_ll_data->getRank(); i>0; --i)
     {
         result += coordinates[i-1]*step_size;
-        step_size *= m_value_axes[i-1].getSize();
+        step_size *= m_value_axes[i-1]->getSize();
     }
     return result;
 }
@@ -399,7 +399,7 @@ template <class T> size_t OutputData<T>::getIndexOfAxis(std::string axis_name, s
 {
     std::vector<int> coordinates = toCoordinates(total_index);
     for (size_t i=0; i<m_value_axes.size(); ++i) {
-        if (m_value_axes[i].getName() == axis_name) {
+        if (m_value_axes[i]->getName() == axis_name) {
             return coordinates[i];
         }
     }
@@ -412,8 +412,8 @@ double OutputData<T>::getValueOfAxis(std::string axis_name, size_t index) const
 {
     std::vector<int> coordinates = toCoordinates(index);
     for (size_t i=0; i<m_value_axes.size(); ++i) {
-        if (m_value_axes[i].getName() == axis_name) {
-            return m_value_axes[i][coordinates[i]];
+        if (m_value_axes[i]->getName() == axis_name) {
+            return (*m_value_axes[i])[coordinates[i]];
         }
     }
     throw LogicErrorException("OutputData<T>::getValueOfAxis() -> Error! Axis with given name not found '"+axis_name+std::string("'"));
@@ -528,10 +528,10 @@ bool OutputData<T>::hasSameShape(const OutputData<T> &right) const
         throw LogicErrorException("OutputData<T>::hasSameShape() -> Panic! Inconsistent dimensions in LLData and axes");
     }
     for (size_t i=0; i<m_value_axes.size(); ++i) {
-        const AxisDouble axis_left = m_value_axes[i];
-        const AxisDouble axis_right = right.m_value_axes[i];
+        const IAxis *p_axis_left = m_value_axes[i];
+        const IAxis *p_axis_right = right.m_value_axes[i];
 
-        if( !HaveSameNameAndShape(axis_left, axis_right)) return false;
+        if( !HaveSameNameAndShape(*p_axis_left, *p_axis_right)) return false;
     }
     return true;
 }
diff --git a/Core/Tools/inc/SafePointerVector.h b/Core/Tools/inc/SafePointerVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..bfdc59626dc1d7d86fca410128100ac4e9e2e7a8
--- /dev/null
+++ b/Core/Tools/inc/SafePointerVector.h
@@ -0,0 +1,111 @@
+#ifndef SAFEPOINTERVECTOR_H_
+#define SAFEPOINTERVECTOR_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   SafePointerVector.h
+//! @brief  Definition of SafePointerVector template
+//! @author Scientific Computing Group at FRM II
+//! @date   Dec 10, 2012
+
+#include <vector>
+
+//- -------------------------------------------------------------------
+//! @class SafePointerVector
+//! @brief Definition of SafePointerVector<T> template for safe handling
+//! of vectors of pointers that are owned by the vector
+//- -------------------------------------------------------------------
+template <class T> class SafePointerVector
+{
+public:
+    typedef typename std::vector<T *>::iterator iterator;
+    typedef typename std::vector<T *>::const_iterator const_iterator;
+    SafePointerVector();
+    SafePointerVector(const SafePointerVector &other);
+    ~SafePointerVector();
+
+    SafePointerVector &operator=(const SafePointerVector &right);
+    size_t size() const;
+    void push_back(T *pointer);
+    T *operator[](size_t index);
+    const T *operator[](size_t index) const;
+    iterator begin() { return m_pointers.begin(); }
+    const_iterator begin() const { return m_pointers.begin(); }
+    iterator end() { return m_pointers.end(); }
+    const_iterator end() const { return m_pointers.end(); }
+
+    T *back() { return m_pointers.back(); }
+    T *const back() const { return m_pointers.back(); }
+    void clear();
+private:
+    std::vector<T *> m_pointers;
+};
+
+template<class T> SafePointerVector<T>::SafePointerVector()
+{
+}
+
+template<class T> SafePointerVector<T>::SafePointerVector(const SafePointerVector<T> &other)
+{
+    for (const_iterator it = other.begin();
+            it != other.end(); ++it) {
+        m_pointers.push_back((*it)->clone());
+    }
+}
+
+template<class T> SafePointerVector<T>::~SafePointerVector()
+{
+    clear();
+}
+
+template<class T> SafePointerVector<T> &SafePointerVector<T>::operator=(
+        const SafePointerVector<T> &right)
+{
+    if (this == &right) return *this;
+    clear();
+    for (const_iterator it = right.begin();
+            it != right.end(); ++it) {
+        m_pointers.push_back((*it)->clone());
+    }
+    return *this;
+}
+
+template<class T>
+inline size_t SafePointerVector<T>::size() const
+{
+    return m_pointers.size();
+}
+
+template<class T>
+inline void SafePointerVector<T>::push_back(T* pointer)
+{
+    m_pointers.push_back(pointer);
+}
+
+template<class T>
+inline T* SafePointerVector<T>::operator[](size_t index)
+{
+    return m_pointers[index];
+}
+
+template<class T>
+inline const T* SafePointerVector<T>::operator[](size_t index) const
+{
+    return m_pointers[index];
+}
+
+template<class T> void SafePointerVector<T>::clear()
+{
+    for (iterator it = begin(); it != end(); ++it) {
+        delete (*it);
+    }
+    m_pointers.clear();
+}
+
+#endif /* SAFEPOINTERVECTOR_H_ */
diff --git a/Core/Tools/src/AxisBin.cpp b/Core/Tools/src/AxisBin.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e9fb04e868ae434ac6cab0aec7ef93ec6ab471e2
--- /dev/null
+++ b/Core/Tools/src/AxisBin.cpp
@@ -0,0 +1,97 @@
+#include "AxisBin.h"
+#include "Numeric.h"
+#include "Exceptions.h"
+
+AxisBin::AxisBin(std::string name)
+: IAxis(name)
+{
+}
+
+AxisBin::AxisBin(std::string name, size_t nbr_bins, double start, double end)
+: IAxis(name)
+{
+    initBins(nbr_bins, start, end);
+}
+
+AxisBin *AxisBin::clone() const
+{
+    AxisBin *p_clone = new AxisBin(getName());
+    p_clone->m_value_vector = this->m_value_vector;
+    return p_clone;
+}
+
+AxisBin *AxisBin::createDoubleBinSize() const
+{
+    if (getSize() < 2) {
+        return clone();
+    }
+    AxisBin *p_result = new AxisBin(getName());
+    for (size_t source_index=0; source_index<getSize(); source_index+=2)
+    {
+        p_result->push_back(m_value_vector[source_index]);
+    }
+    p_result->push_back(m_value_vector[getSize()]);
+    return p_result;
+}
+
+size_t AxisBin::getSize() const
+{
+    if (m_value_vector.size()<2) {
+        return 0;
+    }
+    return m_value_vector.size() - 1;
+}
+
+double AxisBin::operator[](size_t index) const
+{
+    return getBin(index).getMidPoint();
+}
+
+Bin1D AxisBin::getBin(size_t index) const
+{
+    Bin1D result = { m_value_vector[index], m_value_vector[index+1] };
+    return result;
+}
+
+void AxisBin::initBins(size_t nbr_bins, double start, double end)
+{
+    double step = (end - start)/(nbr_bins);
+    for (size_t i=0; i<nbr_bins+1; ++i)
+    {
+        push_back(start + step*(int)i);
+    }
+}
+
+size_t AxisBin::findClosestIndex(double value) const
+{
+    if(m_value_vector.size()<2) {
+        throw ClassInitializationException("AxisBin not (yet) correctly initialized");
+    }
+    if (value < getMin() || value > getMax()) {
+        throw OutOfBoundsException("Given value not in any bin");
+    }
+    std::vector<double>::const_iterator top_limit = std::lower_bound(m_value_vector.begin(), m_value_vector.end(), value);
+    if(top_limit != m_value_vector.begin() ) --top_limit;
+    size_t nbin = top_limit - m_value_vector.begin();
+    return nbin;
+}
+
+//Bin1D AxisBin::findMatchingBin(double value) const
+//{
+//    return (*this)[findMatchingBinIndex(value)];
+//}
+//
+bool AxisBin::equals(const IAxis& other) const
+{
+    if (!IAxis::equals(other)) return false;
+    if (const AxisBin *p_other_cast = dynamic_cast<const AxisBin *>(&other)) {
+        if (getSize() != p_other_cast->getSize()) return false;
+        for(size_t i=0; i<m_value_vector.size(); ++i) {
+            if( std::abs(m_value_vector[i] - p_other_cast->m_value_vector[i]) > Numeric::double_epsilon ) {
+                return false;
+            }
+        }
+        return true;
+    }
+    return false;
+}
diff --git a/Core/Tools/src/AxisDouble.cpp b/Core/Tools/src/AxisDouble.cpp
index 5e55ca963550f0ce813924202c7d1fbc2d7f24ad..729fa2963c5aa668a84bbb49f0ba2bc2dece6cf1 100644
--- a/Core/Tools/src/AxisDouble.cpp
+++ b/Core/Tools/src/AxisDouble.cpp
@@ -1,53 +1,62 @@
 #include "AxisDouble.h"
-#include "BinAxis.h"
+#include "AxisBin.h"
 #include "Numeric.h"
 #include "Exceptions.h"
 
-#include <cmath>
-
 AxisDouble::AxisDouble(std::string name)
-: m_name(name)
+: IAxis(name)
+, m_bin_size(0)
 {
 }
 
 AxisDouble::AxisDouble(std::string name, size_t size, double start, double end)
-: m_name(name)
+: IAxis(name)
+, m_bin_size(0)
 {
     initElements(size, start, end);
 }
 
-AxisDouble::AxisDouble(const BinAxis& source)
-: m_name(source.getName())
+AxisDouble::AxisDouble(const AxisBin& source)
+: IAxis(source.getName())
+, m_bin_size(0)
 {
     for (size_t i=0; i<source.getSize(); ++i) {
-        m_value_vector.push_back(source[i].getMidPoint());
+        m_sample_vector.push_back(source[i]);
     }
 }
 
-AxisDouble* AxisDouble::clone() const
+AxisDouble *AxisDouble::clone() const
 {
     AxisDouble *p_clone = new AxisDouble(getName());
-    p_clone->m_value_vector = this->m_value_vector;
+    p_clone->m_sample_vector = this->m_sample_vector;
+    p_clone->m_bin_size = this->m_bin_size;
     return p_clone;
 }
 
-AxisDouble AxisDouble::createDoubleBinSize() const
+AxisDouble *AxisDouble::createDoubleBinSize() const
 {
     if (getSize() < 2) {
-        return *this;
+        return clone();
     }
-    AxisDouble result(getName());
+    AxisDouble *p_result = new AxisDouble(getName());
     for (size_t source_index=0; source_index<getSize(); source_index+=2)
     {
         double value;
         if (source_index==getSize()-1) {
-            value = (3.0*m_value_vector.at(source_index) - m_value_vector.at(source_index-1))/2.0;
+            value = (3.0*m_sample_vector.at(source_index) - m_sample_vector.at(source_index-1))/2.0;
         }
         else {
-            value =  (m_value_vector.at(source_index) + m_value_vector.at(source_index+1))/2.0;
+            value =  (m_sample_vector.at(source_index) + m_sample_vector.at(source_index+1))/2.0;
         }
-        result.push_back(value);
+        p_result->push_back(value);
     }
+    p_result->m_bin_size = m_bin_size*2.0;
+    return p_result;
+}
+
+Bin1D AxisDouble::getBin(size_t index) const
+{
+    Bin1D result = { m_sample_vector[index] - m_bin_size/2.0, m_sample_vector[index] + m_bin_size/2.0 };
     return result;
 }
 
@@ -62,38 +71,45 @@ void AxisDouble::initElements(size_t size, double start, double end)
 
 size_t AxisDouble::findClosestIndex(double value) const
 {
-    if(m_value_vector.size()<2) return 0;
-    std::vector<double>::const_iterator before = std::lower_bound(m_value_vector.begin(), m_value_vector.end(), value);
-    if(before == m_value_vector.end() ) --before;
-    if(before == m_value_vector.begin() ) ++before;
+    if(m_sample_vector.size()<2) return 0;
+    std::vector<double>::const_iterator before = std::lower_bound(m_sample_vector.begin(), m_sample_vector.end(), value);
+    if(before == m_sample_vector.end() ) --before;
+    if(before == m_sample_vector.begin() ) ++before;
     std::vector<double>::const_iterator after = before;
     --before;
     size_t nbin(0);
-    ( *after-value) < (value - *before) ? nbin = std::distance(m_value_vector.begin(), after) : nbin = std::distance(m_value_vector.begin(), before);
+    ( *after-value) < (value - *before) ? nbin = std::distance(m_sample_vector.begin(), after) : nbin = std::distance(m_sample_vector.begin(), before);
     return nbin;
 }
 
 size_t AxisDouble::getLowerBoundIndex(double value) const
 {
-    if(m_value_vector.size()<2) return 0;
-    std::vector<double>::const_iterator lbound_it = std::lower_bound(m_value_vector.begin(), m_value_vector.end(), value);
-    if(lbound_it == m_value_vector.end() ) {
+    if(m_sample_vector.size()<2) return 0;
+    std::vector<double>::const_iterator lbound_it = std::lower_bound(m_sample_vector.begin(), m_sample_vector.end(), value);
+    if(lbound_it == m_sample_vector.end() ) {
         throw RuntimeErrorException("Given lower bound higher than highest present value");
     }
-    return lbound_it - m_value_vector.begin();
+    return lbound_it - m_sample_vector.begin();
 }
 
 size_t AxisDouble::getUpperBoundIndex(double value) const
 {
-    if(m_value_vector.size()<2) return 0;
-    std::vector<double>::const_iterator lbound_it = std::upper_bound(m_value_vector.begin(), m_value_vector.end(), value);
-    return (lbound_it - m_value_vector.begin()) - 1;
+    if(m_sample_vector.size()<2) return 0;
+    std::vector<double>::const_iterator lbound_it = std::upper_bound(m_sample_vector.begin(), m_sample_vector.end(), value);
+    return (lbound_it - m_sample_vector.begin()) - 1;
 }
 
-bool HaveSameNameAndShape(const AxisDouble& left, const AxisDouble& right)
+bool AxisDouble::equals(const IAxis& other) const
 {
-    if(left.getSize() != right.getSize()) return false;
-    if(left.getName() != right.getName()) return false;
-    for(size_t i=0; i<left.getSize(); ++i) if( std::abs(left[i] - right[i]) > Numeric::double_epsilon)  return false;
-    return true;
+    if (!IAxis::equals(other)) return false;
+    if (const AxisDouble *p_other_cast = dynamic_cast<const AxisDouble *>(&other)) {
+        if (getSize() != p_other_cast->getSize()) return false;
+        for(size_t i=0; i<m_sample_vector.size(); ++i) {
+            if( std::abs(m_sample_vector[i] - p_other_cast->m_sample_vector[i]) > Numeric::double_epsilon ) {
+                return false;
+            }
+        }
+        return true;
+    }
+    return false;
 }
diff --git a/Core/Tools/src/BinAxis.cpp b/Core/Tools/src/BinAxis.cpp
deleted file mode 100644
index 547ff03f93277e8d21fa6a3de155121ab5b5d204..0000000000000000000000000000000000000000
--- a/Core/Tools/src/BinAxis.cpp
+++ /dev/null
@@ -1,100 +0,0 @@
-#include "BinAxis.h"
-#include "Numeric.h"
-#include "Exceptions.h"
-
-#include <cmath>
-
-bool operator==(const Bin1D &left, const Bin1D &right)
-{
-    if (std::abs(left.m_lower - right.m_lower) > Numeric::double_epsilon) return false;
-    if (std::abs(left.m_upper - right.m_upper) > Numeric::double_epsilon) return false;
-    return true;
-}
-
-BinAxis::BinAxis(std::string name)
-: m_name(name)
-{
-}
-
-BinAxis::BinAxis(std::string name, size_t nbr_bins, double start, double end)
-: m_name(name)
-{
-    initBins(nbr_bins, start, end);
-}
-
-BinAxis* BinAxis::clone() const
-{
-    BinAxis *p_clone = new BinAxis(getName());
-    p_clone->m_value_vector = this->m_value_vector;
-    return p_clone;
-}
-
-BinAxis BinAxis::createDoubleBinSize() const
-{
-    if (getSize() < 2) {
-        return *this;
-    }
-    BinAxis result(getName());
-    for (size_t source_index=0; source_index<getSize(); source_index+=2)
-    {
-        result.push_back(m_value_vector[source_index]);
-    }
-    result.push_back(m_value_vector[getSize()]);
-    return result;
-}
-
-size_t BinAxis::getSize() const
-{
-    if (m_value_vector.size()<2) {
-        return 0;
-    }
-    return m_value_vector.size()-1;
-}
-
-Bin1D BinAxis::operator[](size_t index)
-{
-    Bin1D result = { m_value_vector.at(index), m_value_vector.at(index+1) };
-    return result;
-}
-
-const Bin1D BinAxis::operator[](size_t index) const
-{
-    Bin1D result = { m_value_vector.at(index), m_value_vector.at(index+1) };
-    return result;
-}
-
-void BinAxis::initBins(size_t nbr_bins, double start, double end)
-{
-    double step = (end - start)/(nbr_bins);
-    for (size_t i=0; i<nbr_bins+1; ++i)
-    {
-        push_back(start + step*(int)i);
-    }
-}
-
-size_t BinAxis::findMatchingBinIndex(double value) const
-{
-    if(m_value_vector.size()<2) {
-        throw ClassInitializationException("BinAxis not (yet) correctly initialized");
-    }
-    if (value < getMin() || value > getMax()) {
-        throw OutOfBoundsException("Given value not in any bin");
-    }
-    std::vector<double>::const_iterator top_limit = std::lower_bound(m_value_vector.begin(), m_value_vector.end(), value);
-    if(top_limit != m_value_vector.begin() ) --top_limit;
-    size_t nbin = top_limit - m_value_vector.begin();
-    return nbin;
-}
-
-Bin1D BinAxis::findMatchingBin(double value) const
-{
-    return (*this)[findMatchingBinIndex(value)];
-}
-
-bool HaveSameNameAndShape(const BinAxis& left, const BinAxis& right)
-{
-    if(left.getSize() != right.getSize()) return false;
-    if(left.getName() != right.getName()) return false;
-    for(size_t i=0; i<left.getSize(); ++i) if( left[i] != right[i])  return false;
-    return true;
-}
diff --git a/Core/Tools/src/FitSuite.cpp b/Core/Tools/src/FitSuite.cpp
index f508af25b2ab0705f7ce78d98a50f1e7ea170e49..fc55b6e096b93ff56e0013ac3c176d028f9d8649 100644
--- a/Core/Tools/src/FitSuite.cpp
+++ b/Core/Tools/src/FitSuite.cpp
@@ -163,7 +163,7 @@ double FitSuite::functionToMinimize(const double *pars_current_values)
     m_fit_objects.runSimulation();
 
     // caclulate chi2 value
-    double chi_squared = m_fit_objects.getChiSquaredValue(m_fit_parameters.getNfreeParameters());
+    double chi_squared = m_fit_objects.getChiSquaredValue((int)m_fit_parameters.getNfreeParameters());
 
     notifyObservers();
     m_n_call++;
diff --git a/Core/Tools/src/FitSuiteObjects.cpp b/Core/Tools/src/FitSuiteObjects.cpp
index 2e1493edd87cd0bdba43b1c773cb5092e9d5893f..fb1ca886e5b85fb194233c0bbf6cc2202bff699b 100644
--- a/Core/Tools/src/FitSuiteObjects.cpp
+++ b/Core/Tools/src/FitSuiteObjects.cpp
@@ -51,7 +51,7 @@ double FitSuiteObjects::getChiSquaredValue(int n_free_fit_parameters)
     for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
         IChiSquaredModule *chi = (*it)->getChiSquaredModule();
 
-        chi->setNdegreeOfFreedom( m_fit_objects.size() * (*it)->getRealData()->getAllocatedSize() - n_free_fit_parameters);
+        chi->setNdegreeOfFreedom( (int)(m_fit_objects.size() * (*it)->getRealData()->getAllocatedSize() - n_free_fit_parameters) );
         // normalizing datasets to the maximum intensity over all fit objects defined
         OutputDataNormalizerScaleAndShift *data_normalizer =  dynamic_cast<OutputDataNormalizerScaleAndShift *>(chi->getOutputDataNormalizer());
         if( data_normalizer) data_normalizer->setMaximumIntensity( getSimulationMaxIntensity() );
diff --git a/Core/Tools/src/OutputDataFunctions.cpp b/Core/Tools/src/OutputDataFunctions.cpp
index 0971cc09d341daa9737a19c3d4e45265b3756023..84ea350454a68f6281b15b28dad83c17bab1ed0a 100644
--- a/Core/Tools/src/OutputDataFunctions.cpp
+++ b/Core/Tools/src/OutputDataFunctions.cpp
@@ -23,8 +23,10 @@ OutputData<double> *OutputDataFunctions::doubleBinSize(const OutputData<double>
     // create new axes
     for (size_t i=0; i<dimension; ++i) {
         needs_resizing.push_back(source_sizes[i] > 1);
-        const AxisDouble *source_axis = source.getAxis(i);
-        p_result->addAxis(source_axis->createDoubleBinSize());
+        const IAxis *source_axis = source.getAxis(i);
+        IAxis *p_new_axis = source_axis->createDoubleBinSize();
+        p_result->addAxis(*p_new_axis);
+        delete p_new_axis;
     }
     // calculate new data content
     OutputData<double>::const_iterator it_source = source.begin();
@@ -208,10 +210,10 @@ OutputData<double> *OutputDataFunctions::sliceAccrossOneAxis(const OutputData<do
 
     OutputData<double > *sliced_data = new OutputData<double >;
 
-    const AxisDouble *fixed_axis(0);
+    const IAxis *fixed_axis(0);
     int fixed_axis_index(-1);
     for(size_t i_axis=0; i_axis<data.getNdimensions(); i_axis++) {
-        const AxisDouble *axis = data.getAxis(i_axis);
+        const IAxis *axis = data.getAxis(i_axis);
         if( axis->getName() != fixed_axis_name ) {
             sliced_data->addAxis(*axis);
         } else {
@@ -250,7 +252,7 @@ OutputData<double> *OutputDataFunctions::selectRangeOnOneAxis(const OutputData<d
         throw LogicErrorException("OutputDataFunctions::selectRangeOnOneAxis() -> Error! It was checked only with number of dimensions equal 2.");
     }
 
-    const AxisDouble *selected_axis = data.getAxis(selected_axis_name);
+    const IAxis *selected_axis = data.getAxis(selected_axis_name);
     if( !selected_axis ) {
         throw LogicErrorException("OutputDataFunctions::selectRangeOnOneAxis() -> Error! No axis with name "+selected_axis_name);
     }
@@ -268,7 +270,7 @@ OutputData<double> *OutputDataFunctions::selectRangeOnOneAxis(const OutputData<d
     // preparing new data with modified axes
     OutputData<double > *new_data = new OutputData<double >;
     for(size_t i_axis=0; i_axis<data.getNdimensions(); i_axis++) {
-        const AxisDouble *axis = data.getAxis(i_axis);
+        const IAxis *axis = data.getAxis(i_axis);
         if( axis->getName() != selected_axis_name ) {
             new_data->addAxis(*axis);
         } else {
@@ -341,9 +343,9 @@ Mask* OutputDataFunctions::CreateRectangularMask(const OutputData<double>& data,
     int *maxima_i = new int[rank];
     int *dims_i = new int[rank];
     for (size_t i=0; i<rank; ++i) {
-        const AxisDouble *p_axis = data.getAxis(i);
-        minima_i[i] = (int)p_axis->getLowerBoundIndex(minima[i]);
-        maxima_i[i] = (int)p_axis->getUpperBoundIndex(maxima[i]);
+        const IAxis *p_axis = data.getAxis(i);
+        minima_i[i] = (int)p_axis->findClosestIndex(minima[i]);
+        maxima_i[i] = (int)p_axis->findClosestIndex(maxima[i]);
         dims_i[i] = (int)p_axis->getSize();
     }
     MaskCoordinateRectangleFunction *p_rectangle_function = new MaskCoordinateRectangleFunction(rank, minima_i, maxima_i);
@@ -364,9 +366,9 @@ Mask* OutputDataFunctions::CreateEllipticMask(const OutputData<double>& data,
     int *radii_i = new int[rank];
     int *dims_i = new int[rank];
     for (size_t i=0; i<rank; ++i) {
-        const AxisDouble *p_axis = data.getAxis(i);
-        center_i[i] = (int)p_axis->getLowerBoundIndex(center[i]);
-        int lower_index = (int)p_axis->getLowerBoundIndex((*p_axis)[center_i[i]] - radii[i]);
+        const IAxis *p_axis = data.getAxis(i);
+        center_i[i] = (int)p_axis->findClosestIndex(center[i]);
+        int lower_index = (int)p_axis->findClosestIndex((*p_axis)[center_i[i]] - radii[i]);
         radii_i[i] = center_i[i] - lower_index;
         dims_i[i] = (int)p_axis->getSize();
     }
diff --git a/UnitTests/TestCore/DetectorTest.h b/UnitTests/TestCore/DetectorTest.h
index 89be1d35b359d6b96280d46b18c7ea77742b40a8..b6d5def59ce449be9d248e596289b66b28f925b3 100644
--- a/UnitTests/TestCore/DetectorTest.h
+++ b/UnitTests/TestCore/DetectorTest.h
@@ -57,9 +57,10 @@ TEST_F(DetectorTest, DetectorConstruction)
     AxisDouble axis1("axis1", 20, 0.0, 20.0);
     constructedDetector.addAxis(axis0);
     constructedDetector.addAxis(axis1);
+
     EXPECT_EQ((size_t)2, constructedDetector.getDimension());
-    AxisDouble axis0copy = constructedDetector.getAxis(0);
-    AxisDouble axis1copy = constructedDetector.getAxis(1);
+    const IAxis &axis0copy = constructedDetector.getAxis(0);
+    const IAxis &axis1copy = constructedDetector.getAxis(1);
     ASSERT_TRUE(axis0.getMin() == axis0copy.getMin() );
     ASSERT_TRUE(axis0.getMax() == axis0copy.getMax() );
     ASSERT_TRUE(axis0.getName() == axis0copy.getName() );