diff --git a/Base/Types/OwningVector.h b/Base/Types/OwningVector.h
index a13fdd556c1b9b23ad726aca97d7b9660f489bb7..8daf571bc204d35f159d8af86ad554f0723e0d03 100644
--- a/Base/Types/OwningVector.h
+++ b/Base/Types/OwningVector.h
@@ -35,13 +35,21 @@ template <class T>
 class OwningVector {
 public:
     OwningVector() = default;
-    ~OwningVector() { clear(); }
+    //! Constructor that takes over ownership of elements in given vector
+    OwningVector(const std::vector<T*>& v)
+    {
+        m_v.reserve(v.size());
+        for (T* e : v)
+            m_v.emplace_back(e);
+    }
+    //! Constructor that clones elements in given vector
     OwningVector(const OwningVector& other)
     {
         m_v.reserve(other.size());
         for (T* e : other)
             m_v.emplace_back(e->clone());
     }
+    ~OwningVector() { clear(); }
     OwningVector& operator=(const OwningVector& other)
     {
         if (this == &other)
diff --git a/Device/Data/Powerfield.cpp b/Device/Data/Powerfield.cpp
index 798184a8630534ea027fd678723586f659db7fe0..c01305eb8cb815aacc00768cf4652f7f8ab4abff 100644
--- a/Device/Data/Powerfield.cpp
+++ b/Device/Data/Powerfield.cpp
@@ -19,6 +19,118 @@
 
 #include "Base/Py/PyCore.h"
 
+IField::IField(const std::vector<IAxis*>& axes)
+    : m_axes(axes)
+{
+}
+
+const IAxis& IField::axis(size_t serial_number) const
+{
+    return *m_axes[serial_number];
+}
+
+double IField::getAxisValue(size_t global_index, size_t i_selected_axis) const
+{
+    auto axis_index = getAxisBinIndex(global_index, i_selected_axis);
+    return (*m_axes[i_selected_axis])[axis_index];
+}
+
+double IField::getAxisValue(size_t global_index, const std::string& axis_name) const
+{
+    return getAxisValue(global_index, getAxisIndex(axis_name));
+}
+
+std::vector<double> IField::getAxesValues(size_t global_index) const
+{
+    std::vector<int> indices = getAxesBinIndices(global_index);
+    std::vector<double> result;
+    for (size_t i_axis = 0; i_axis < indices.size(); ++i_axis)
+        result.push_back((*m_axes[i_axis])[indices[i_axis]]);
+    return result;
+}
+
+Bin1D IField::getAxisBin(size_t global_index, size_t i_selected_axis) const
+{
+    auto axis_index = getAxisBinIndex(global_index, i_selected_axis);
+    return m_axes[i_selected_axis]->bin(axis_index);
+}
+
+Bin1D IField::getAxisBin(size_t global_index, const std::string& axis_name) const
+{
+    return getAxisBin(global_index, getAxisIndex(axis_name));
+}
+
+size_t IField::getAxisIndex(const std::string& axis_name) const
+{
+    for (size_t i = 0; i < m_axes.size(); ++i)
+        if (m_axes[i]->axisName() == axis_name)
+            return i;
+    ASSERT(0);
+}
+
+size_t IField::getAxisBinIndex(size_t global_index, const std::string& axis_name) const
+{
+    return getAxisBinIndex(global_index, getAxisIndex(axis_name));
+}
+
+std::vector<int> IField::getAxesBinIndices(size_t global_index) const
+{
+    size_t remainder = global_index;
+    std::vector<int> result;
+    result.resize(rank());
+    for (size_t i = 0; i < rank(); ++i) {
+        result[rank() - 1 - i] = (int)(remainder % m_axes[rank() - 1 - i]->size());
+        remainder /= m_axes[rank() - 1 - i]->size();
+    }
+    return result;
+}
+
+size_t IField::getAxisBinIndex(size_t global_index, size_t i_selected_axis) const
+{
+    size_t remainder(global_index);
+    for (size_t i = 0; i < rank(); ++i) {
+        size_t i_axis = rank() - 1 - i;
+        size_t result = remainder % m_axes[i_axis]->size();
+        if (i_selected_axis == i_axis)
+            return result;
+        remainder /= m_axes[i_axis]->size();
+    }
+    ASSERT(0);
+}
+
+size_t IField::toGlobalIndex(const std::vector<unsigned>& axes_indices) const
+{
+    ASSERT(axes_indices.size() == rank());
+    size_t result = 0;
+    size_t step_size = 1;
+    for (size_t i = rank(); i > 0; --i) {
+        ASSERT(axes_indices[i - 1] < m_axes[i - 1]->size());
+        result += axes_indices[i - 1] * step_size;
+        step_size *= m_axes[i - 1]->size();
+    }
+    return result;
+}
+
+size_t IField::findGlobalIndex(const std::vector<double>& coordinates) const
+{
+    ASSERT(coordinates.size() == rank());
+    std::vector<unsigned> axes_indexes;
+    axes_indexes.resize(rank());
+    for (size_t i = 0; i < rank(); ++i)
+        axes_indexes[i] = static_cast<unsigned>(m_axes[i]->findClosestIndex(coordinates[i]));
+    return toGlobalIndex(axes_indexes);
+}
+
+
+bool IField::axisNameExists(const std::string& axis_name) const
+{
+    for (size_t i = 0; i < m_axes.size(); ++i)
+        if (m_axes[i]->axisName() == axis_name)
+            return true;
+    return false;
+}
+
+
 template <>
 PyObject* Powerfield<double>::getArray() const
 {
@@ -77,7 +189,7 @@ double Powerfield<CumulativeValue>::getValue(size_t index) const
 template <>
 const Powerfield<double>& Powerfield<double>::normalizedToMaximum()
 {
-    double maxval = max();
+    double maxval = m_ll_data->max();
     ASSERT(maxval > 0);
     for (size_t i = 0; i < m_ll_data->getTotalSize(); ++i)
         (*this)[i] /= maxval;
diff --git a/Device/Data/Powerfield.h b/Device/Data/Powerfield.h
index 41b897a9412b4c1b09afcb6bc473a3432153893d..f7da207a0e559d96252df4dc58fb54c87ad86186 100644
--- a/Device/Data/Powerfield.h
+++ b/Device/Data/Powerfield.h
@@ -27,20 +27,104 @@
 
 using std::size_t;
 
+//! Holds one or two axes. Base class for Powerfield.
+
+class IField {
+public:
+    IField() = default;
+    IField(const std::vector<IAxis*>& axes);
+    virtual ~IField() {}
+
+    //! Returns number of dimensions.
+    size_t rank() const { return m_axes.size(); }
+
+    //! Returns axis with given serial number
+    const IAxis& axis(size_t serial_number) const;
+
+    //! Returns the value of selected axis for given global_index.
+    //! @param global_index The global index of this data structure.
+    //! @param i_selected_axis Serial number of selected axis.
+    //! @return corresponding bin center of selected axis
+    double getAxisValue(size_t global_index, size_t i_selected_axis) const;
+
+    //! Returns the value of selected axis for given global_index.
+    //! @param global_index The global index of this data structure.
+    //! @param axis_name The name of selected axis.
+    //! @return corresponding bin center of selected axis
+    double getAxisValue(size_t global_index, const std::string& axis_name) const;
+
+    //! Returns values on all defined axes for given globalbin number
+    //! @param global_index The global index of this data structure.
+    //! @return Vector of corresponding bin centers
+    std::vector<double> getAxesValues(size_t global_index) const;
+
+    //! Returns bin of selected axis for given global_index.
+    //! @param global_index The global index of this data structure.
+    //! @param i_selected_axis Serial number of selected axis.
+    //! @return Corresponding Bin1D object
+    Bin1D getAxisBin(size_t global_index, size_t i_selected_axis) const;
+
+    //! Returns bin of selected axis for given global_index.
+    //! @param global_index The global index of this data structure.
+    //! @param axis_name The name of selected axis.
+    //! @return Corresponding Bin1D object
+    Bin1D getAxisBin(size_t global_index, const std::string& axis_name) const;
+
+    //! Returns axis bin index for given global index
+    //! @param global_index The global index of this data structure.
+    //! @param axis_name The name of selected axis.
+    //! @return Corresponding bin index for selected axis
+    size_t getAxisBinIndex(size_t global_index, const std::string& axis_name) const;
+
+    //! Returns vector of axes indices for given global index
+    //! @param global_index The global index of this data structure.
+    //! @return Vector of bin indices for all axes defined
+    std::vector<int> getAxesBinIndices(size_t global_index) const;
+
+    //! Returns axis bin index for given global index
+    //! @param global_index The global index of this data structure.
+    //! @param i_selected_axis Serial number of selected axis.
+    //! @return Corresponding bin index for selected axis
+    size_t getAxisBinIndex(size_t global_index, size_t i_selected_axis) const;
+
+
+    //! Returns global index for specified indices of axes
+    //! @param axes_indices Vector of axes indices for all specified axes in this dataset
+    //! @return Corresponding global index
+    size_t toGlobalIndex(const std::vector<unsigned>& axes_indices) const;
+
+    //! Returns global index for specified axes values
+    //! @param coordinates Vector of axes coordinates for all specified axes in this dataset
+    //! @return Closest global index
+    size_t findGlobalIndex(const std::vector<double>& coordinates) const;
+
+
+protected:
+    //! Returns serial number of axis with given name
+    size_t getAxisIndex(const std::string& axis_name) const;
+
+    //! Checks if given axis name exists
+    bool axisNameExists(const std::string& axis_name) const;
+
+    OwningVector<IAxis> m_axes;
+};
+
 //! Templated class to store data in multi-dimensional space.
 //! Used with type T = double, CumulativeValue, bool
 //! @ingroup tools
 
 template <class T>
-class Powerfield {
+class Powerfield : public IField {
 public:
     using value_type = T;
 
     Powerfield();
+    Powerfield(const IAxis& xAxis);
+    Powerfield(const IAxis& xAxis, const IAxis& yAxis);
     Powerfield(const Powerfield&) = delete;
     Powerfield(Powerfield&&);
     const Powerfield& operator=(const Powerfield&) = delete;
-    ~Powerfield();
+    ~Powerfield() override;
     Powerfield* clone() const;
 
     void copyFrom(const Powerfield<T>& other);
@@ -51,14 +135,8 @@ public:
     void addAxis(const IAxis& new_axis);
     void addAxis(const std::string& name, size_t size, double start, double end);
 
-    //! Returns axis with given serial number
-    const IAxis& axis(size_t serial_number) const;
-
     // retrieve basic info
 
-    //! Returns number of dimensions.
-    size_t rank() const { return m_axes.size(); }
-
     //! Returns total size of data buffer (product of bin number in every dimension).
     size_t getAllocatedSize() const
     {
@@ -76,9 +154,6 @@ public:
     //! Returns sum of all values in the data structure
     T totalSum() const;
 
-    //! Returns maximum value in the data structure
-    T max() const;
-
     // iterators
 
     friend class PowerfieldIterator<T, Powerfield<T>>;
@@ -104,62 +179,6 @@ public:
 
     // coordinate and index functions
 
-    //! Returns vector of axes indices for given global index
-    //! @param global_index The global index of this data structure.
-    //! @return Vector of bin indices for all axes defined
-    std::vector<int> getAxesBinIndices(size_t global_index) const;
-
-    //! Returns axis bin index for given global index
-    //! @param global_index The global index of this data structure.
-    //! @param i_selected_axis Serial number of selected axis.
-    //! @return Corresponding bin index for selected axis
-    size_t getAxisBinIndex(size_t global_index, size_t i_selected_axis) const;
-
-    //! Returns axis bin index for given global index
-    //! @param global_index The global index of this data structure.
-    //! @param axis_name The name of selected axis.
-    //! @return Corresponding bin index for selected axis
-    size_t getAxisBinIndex(size_t global_index, const std::string& axis_name) const;
-
-    //! Returns global index for specified indices of axes
-    //! @param axes_indices Vector of axes indices for all specified axes in this dataset
-    //! @return Corresponding global index
-    size_t toGlobalIndex(const std::vector<unsigned>& axes_indices) const;
-
-    //! Returns global index for specified axes values
-    //! @param coordinates Vector of axes coordinates for all specified axes in this dataset
-    //! @return Closest global index
-    size_t findGlobalIndex(const std::vector<double>& coordinates) const;
-
-    //! Returns the value of selected axis for given global_index.
-    //! @param global_index The global index of this data structure.
-    //! @param i_selected_axis Serial number of selected axis.
-    //! @return corresponding bin center of selected axis
-    double getAxisValue(size_t global_index, size_t i_selected_axis) const;
-
-    //! Returns the value of selected axis for given global_index.
-    //! @param global_index The global index of this data structure.
-    //! @param axis_name The name of selected axis.
-    //! @return corresponding bin center of selected axis
-    double getAxisValue(size_t global_index, const std::string& axis_name) const;
-
-    //! Returns values on all defined axes for given globalbin number
-    //! @param global_index The global index of this data structure.
-    //! @return Vector of corresponding bin centers
-    std::vector<double> getAxesValues(size_t global_index) const;
-
-    //! Returns bin of selected axis for given global_index.
-    //! @param global_index The global index of this data structure.
-    //! @param i_selected_axis Serial number of selected axis.
-    //! @return Corresponding Bin1D object
-    Bin1D getAxisBin(size_t global_index, size_t i_selected_axis) const;
-
-    //! Returns bin of selected axis for given global_index.
-    //! @param global_index The global index of this data structure.
-    //! @param axis_name The name of selected axis.
-    //! @return Corresponding Bin1D object
-    Bin1D getAxisBin(size_t global_index, const std::string& axis_name) const;
-
     // modifiers
 
     //! Sets object into initial state (no dimensions, data)
@@ -168,9 +187,6 @@ public:
     //! Sets content of output data to specific value
     void setAllTo(const T& value);
 
-    //! Adds 'rank' axes with indicated sizes
-    void setAxisSizes(size_t rank, int* n_dims);
-
     //! Sets new values to raw data vector
     void setRawDataVector(const std::vector<T>& data_vector);
 
@@ -178,16 +194,16 @@ public:
     void setRawDataArray(const T* source);
 
     //! addition-assignment operator for two output data
-    const Powerfield<T>& operator+=(const Powerfield<T>& right);
+    const Powerfield<T>& operator+=(const Powerfield<T>& other);
 
     //! substraction-assignment operator for two output data
-    const Powerfield<T>& operator-=(const Powerfield<T>& right);
+    const Powerfield<T>& operator-=(const Powerfield<T>& other);
 
     //! division-assignment operator for two output data
-    const Powerfield<T>& operator/=(const Powerfield<T>& right);
+    const Powerfield<T>& operator/=(const Powerfield<T>& other);
 
     //! multiplication-assignment operator for two output data
-    const Powerfield<T>& operator*=(const Powerfield<T>& right);
+    const Powerfield<T>& operator*=(const Powerfield<T>& other);
 
     //! Returns value or summed value, depending on T
     double getValue(size_t index) const;
@@ -210,11 +226,11 @@ public:
 
     //! Returns true if object have same dimensions and number of axes bins
     template <class U>
-    bool hasSameDimensions(const Powerfield<U>& right) const;
+    bool hasSameDimensions(const Powerfield<U>& other) const;
 
     //! Returns true if objects a) have same dimensions b) bin boundaries of axes coincide
     template <class U>
-    bool hasSameShape(const Powerfield<U>& right) const;
+    bool hasSameShape(const Powerfield<U>& other) const;
 
 #ifdef BORNAGAIN_PYTHON
     //! Returns data as Python numpy array
@@ -230,14 +246,7 @@ public:
     const Powerfield<double>& normalizedToMaximum();
 
 private:
-    //! Returns serial number of axis with given name
-    size_t getAxisIndex(const std::string& axis_name) const;
-
-    //! Checks if given axis name exists
-    bool axisNameExists(const std::string& axis_name) const;
-
-    OwningVector<IAxis> m_axes;
-    LLData<T>* m_ll_data;
+    LLData<T>* m_ll_data{nullptr};
 };
 
 #ifndef USER_API
@@ -247,22 +256,40 @@ private:
 // --------------------------------------------------------------------- //
 
 template <class T>
-Powerfield<T>::Powerfield(Powerfield<T>&& other)
+Powerfield<T>::Powerfield()
+    : m_ll_data(nullptr)
 {
-    m_axes = std::move(other.m_axes);
-    m_ll_data = other.m_ll_data;
+    allocate();
+}
 
-    other.m_ll_data = nullptr;
-    other.allocate();
+template <class T>
+Powerfield<T>::Powerfield(const IAxis& xAxis)
+    : IField({xAxis.clone()})
+{
+    ASSERT(xAxis.size() > 0);
+    allocate();
 }
 
 template <class T>
-Powerfield<T>::Powerfield()
-    : m_ll_data(nullptr)
+Powerfield<T>::Powerfield(const IAxis& xAxis, const IAxis& yAxis)
+    : IField({xAxis.clone(), yAxis.clone()})
 {
+    ASSERT(xAxis.size() > 0);
+    ASSERT(yAxis.size() > 0);
+    ASSERT(xAxis.axisName() != yAxis.axisName());
     allocate();
 }
 
+template <class T>
+Powerfield<T>::Powerfield(Powerfield<T>&& other)
+    : IField(other)
+{
+    m_ll_data = other.m_ll_data;
+
+    other.m_ll_data = nullptr;
+    other.allocate();
+}
+
 template <class T>
 Powerfield<T>::~Powerfield()
 {
@@ -315,10 +342,9 @@ template <class T>
 void Powerfield<T>::addAxis(const IAxis& new_axis)
 {
     ASSERT(!axisNameExists(new_axis.axisName()));
-    if (new_axis.size() > 0) {
-        m_axes.emplace_back(new_axis.clone());
-        allocate();
-    }
+    ASSERT(new_axis.size() > 0);
+    m_axes.emplace_back(new_axis.clone());
+    allocate();
 }
 
 template <class T>
@@ -329,12 +355,6 @@ void Powerfield<T>::addAxis(const std::string& name, size_t size, double start,
     addAxis(new_axis);
 }
 
-template <class T>
-const IAxis& Powerfield<T>::axis(size_t serial_number) const
-{
-    return *m_axes[serial_number];
-}
-
 template <class T>
 inline std::vector<size_t> Powerfield<T>::getAllSizes() const
 {
@@ -371,105 +391,6 @@ typename Powerfield<T>::const_iterator Powerfield<T>::begin() const
     return result;
 }
 
-template <class T>
-std::vector<int> Powerfield<T>::getAxesBinIndices(size_t global_index) const
-{
-    ASSERT(m_ll_data);
-    size_t remainder = global_index;
-    std::vector<int> result;
-    result.resize(m_ll_data->rank());
-    for (size_t i = 0; i < m_ll_data->rank(); ++i) {
-        result[m_ll_data->rank() - 1 - i] =
-            (int)(remainder % m_axes[m_ll_data->rank() - 1 - i]->size());
-        remainder /= m_axes[m_ll_data->rank() - 1 - i]->size();
-    }
-    return result;
-}
-
-template <class T>
-size_t Powerfield<T>::getAxisBinIndex(size_t global_index, size_t i_selected_axis) const
-{
-    ASSERT(m_ll_data);
-    size_t remainder(global_index);
-    for (size_t i = 0; i < m_ll_data->rank(); ++i) {
-        size_t i_axis = m_ll_data->rank() - 1 - i;
-        size_t result = remainder % m_axes[i_axis]->size();
-        if (i_selected_axis == i_axis)
-            return result;
-        remainder /= m_axes[i_axis]->size();
-    }
-    ASSERT(0);
-}
-
-template <class T>
-size_t Powerfield<T>::getAxisBinIndex(size_t global_index, const std::string& axis_name) const
-{
-    return getAxisBinIndex(global_index, getAxisIndex(axis_name));
-}
-
-template <class T>
-size_t Powerfield<T>::toGlobalIndex(const std::vector<unsigned>& axes_indices) const
-{
-    ASSERT(m_ll_data);
-    ASSERT(axes_indices.size() == m_ll_data->rank());
-    size_t result = 0;
-    size_t step_size = 1;
-    for (size_t i = m_ll_data->rank(); i > 0; --i) {
-        ASSERT(axes_indices[i - 1] < m_axes[i - 1]->size());
-        result += axes_indices[i - 1] * step_size;
-        step_size *= m_axes[i - 1]->size();
-    }
-    return result;
-}
-
-template <class T>
-size_t Powerfield<T>::findGlobalIndex(const std::vector<double>& coordinates) const
-{
-    ASSERT(m_ll_data);
-    ASSERT(coordinates.size() == m_ll_data->rank());
-    std::vector<unsigned> axes_indexes;
-    axes_indexes.resize(m_ll_data->rank());
-    for (size_t i = 0; i < m_ll_data->rank(); ++i)
-        axes_indexes[i] = static_cast<unsigned>(m_axes[i]->findClosestIndex(coordinates[i]));
-    return toGlobalIndex(axes_indexes);
-}
-
-template <class T>
-double Powerfield<T>::getAxisValue(size_t global_index, size_t i_selected_axis) const
-{
-    auto axis_index = getAxisBinIndex(global_index, i_selected_axis);
-    return (*m_axes[i_selected_axis])[axis_index];
-}
-
-template <class T>
-double Powerfield<T>::getAxisValue(size_t global_index, const std::string& axis_name) const
-{
-    return getAxisValue(global_index, getAxisIndex(axis_name));
-}
-
-template <class T>
-std::vector<double> Powerfield<T>::getAxesValues(size_t global_index) const
-{
-    std::vector<int> indices = getAxesBinIndices(global_index);
-    std::vector<double> result;
-    for (size_t i_axis = 0; i_axis < indices.size(); ++i_axis)
-        result.push_back((*m_axes[i_axis])[indices[i_axis]]);
-    return result;
-}
-
-template <class T>
-Bin1D Powerfield<T>::getAxisBin(size_t global_index, size_t i_selected_axis) const
-{
-    auto axis_index = getAxisBinIndex(global_index, i_selected_axis);
-    return m_axes[i_selected_axis]->bin(axis_index);
-}
-
-template <class T>
-Bin1D Powerfield<T>::getAxisBin(size_t global_index, const std::string& axis_name) const
-{
-    return getAxisBin(global_index, getAxisIndex(axis_name));
-}
-
 template <class T>
 inline T Powerfield<T>::totalSum() const
 {
@@ -477,13 +398,6 @@ inline T Powerfield<T>::totalSum() const
     return m_ll_data->getTotalSum();
 }
 
-template <class T>
-inline T Powerfield<T>::max() const
-{
-    ASSERT(m_ll_data);
-    return m_ll_data->max();
-}
-
 template <class T>
 void Powerfield<T>::clear()
 {
@@ -499,38 +413,26 @@ void Powerfield<T>::setAllTo(const T& value)
 }
 
 template <class T>
-void Powerfield<T>::setAxisSizes(size_t rank, int* n_dims)
-{
-    clear();
-    std::string basename("axis");
-    for (size_t i = 0; i < rank; ++i) {
-        std::ostringstream name;
-        name << basename << i;
-        addAxis(name.str(), n_dims[i], 0.0, (double)(n_dims[i] - 1));
-    }
-}
-
-template <class T>
-const Powerfield<T>& Powerfield<T>::operator+=(const Powerfield<T>& right)
+const Powerfield<T>& Powerfield<T>::operator+=(const Powerfield<T>& other)
 {
     ASSERT(m_ll_data);
-    *this->m_ll_data += *right.m_ll_data;
+    *this->m_ll_data += *other.m_ll_data;
     return *this;
 }
 
 template <class T>
-const Powerfield<T>& Powerfield<T>::operator-=(const Powerfield<T>& right)
+const Powerfield<T>& Powerfield<T>::operator-=(const Powerfield<T>& other)
 {
     ASSERT(m_ll_data);
-    *this->m_ll_data -= *right.m_ll_data;
+    *this->m_ll_data -= *other.m_ll_data;
     return *this;
 }
 
 template <class T>
-const Powerfield<T>& Powerfield<T>::operator*=(const Powerfield<T>& right)
+const Powerfield<T>& Powerfield<T>::operator*=(const Powerfield<T>& other)
 {
     ASSERT(m_ll_data);
-    *this->m_ll_data *= *right.m_ll_data;
+    *this->m_ll_data *= *other.m_ll_data;
     return *this;
 }
 
@@ -547,10 +449,10 @@ bool Powerfield<T>::isInitialized() const
 }
 
 template <class T>
-const Powerfield<T>& Powerfield<T>::operator/=(const Powerfield<T>& right)
+const Powerfield<T>& Powerfield<T>::operator/=(const Powerfield<T>& other)
 {
     ASSERT(m_ll_data);
-    *this->m_ll_data /= *right.m_ll_data;
+    *this->m_ll_data /= *other.m_ll_data;
     return *this;
 }
 
@@ -586,16 +488,16 @@ inline void Powerfield<T>::setRawDataArray(const T* source)
 //! Returns true if object have same dimensions
 template <class T>
 template <class U>
-inline bool Powerfield<T>::hasSameDimensions(const Powerfield<U>& right) const
+inline bool Powerfield<T>::hasSameDimensions(const Powerfield<U>& other) const
 {
     if (!isInitialized())
         return false;
-    if (!right.isInitialized())
+    if (!other.isInitialized())
         return false;
-    if (rank() != right.rank())
+    if (rank() != other.rank())
         return false;
     for (size_t i_axis = 0; i_axis < rank(); ++i_axis)
-        if (axis(i_axis).size() != right.axis(i_axis).size())
+        if (axis(i_axis).size() != other.axis(i_axis).size())
             return false;
     return true;
 }
@@ -603,13 +505,13 @@ inline bool Powerfield<T>::hasSameDimensions(const Powerfield<U>& right) const
 //! Returns true if object have same dimensions and shape of axis
 template <class T>
 template <class U>
-bool Powerfield<T>::hasSameShape(const Powerfield<U>& right) const
+bool Powerfield<T>::hasSameShape(const Powerfield<U>& other) const
 {
-    if (!hasSameDimensions(right))
+    if (!hasSameDimensions(other))
         return false;
 
     for (size_t i = 0; i < m_axes.size(); ++i)
-        if (axis(i) != right.axis(i))
+        if (axis(i) != other.axis(i))
             return false;
     return true;
 }
@@ -620,24 +522,5 @@ template <>
 PyObject* Powerfield<double>::getArray() const;
 #endif
 
-//! Returns index of axis
-template <class T>
-size_t Powerfield<T>::getAxisIndex(const std::string& axis_name) const
-{
-    for (size_t i = 0; i < m_axes.size(); ++i)
-        if (m_axes[i]->axisName() == axis_name)
-            return i;
-    ASSERT(0);
-}
-
-template <class T>
-bool Powerfield<T>::axisNameExists(const std::string& axis_name) const
-{
-    for (size_t i = 0; i < m_axes.size(); ++i)
-        if (m_axes[i]->axisName() == axis_name)
-            return true;
-    return false;
-}
-
 #endif // USER_API
 #endif // BORNAGAIN_DEVICE_DATA_POWERFIELD_H
diff --git a/Device/Detector/IDetector2D.cpp b/Device/Detector/IDetector2D.cpp
index 19e502ff91308c9fc64566117ad76c81b04cf60c..35b074e3a210f4263e5a474de6aff9ff922083e2 100644
--- a/Device/Detector/IDetector2D.cpp
+++ b/Device/Detector/IDetector2D.cpp
@@ -39,9 +39,7 @@ IDetector2D::IDetector2D()
 {
 }
 
-IDetector2D::IDetector2D(const IDetector2D& other)
-
-    = default;
+IDetector2D::IDetector2D(const IDetector2D& other) = default;
 
 IDetector2D::~IDetector2D() = default;
 
diff --git a/Device/Histo/Histogram1D.cpp b/Device/Histo/Histogram1D.cpp
index ef6dfe84d10643c5e9c88021d8f56c447dd1984a..e80f585e00c2e614896f7936d162c44f3c6b2a84 100644
--- a/Device/Histo/Histogram1D.cpp
+++ b/Device/Histo/Histogram1D.cpp
@@ -18,13 +18,13 @@
 #include <memory>
 
 Histogram1D::Histogram1D(int nbinsx, double xlow, double xup)
+    : IHistogram(FixedBinAxis("x-axis", nbinsx, xlow, xup))
 {
-    m_data.addAxis(FixedBinAxis("x-axis", nbinsx, xlow, xup));
 }
 
 Histogram1D::Histogram1D(int nbinsx, const std::vector<double>& xbins)
+    : IHistogram(VariableBinAxis("x-axis", nbinsx, xbins))
 {
-    m_data.addAxis(VariableBinAxis("x-axis", nbinsx, xbins));
 }
 
 Histogram1D::Histogram1D(const IAxis& axis)
diff --git a/Device/Histo/Histogram2D.cpp b/Device/Histo/Histogram2D.cpp
index b8beca54ffca08b193c6cd5b852bee53a59dffe4..343f638de82986d4605b145786b426df5e83b6f7 100644
--- a/Device/Histo/Histogram2D.cpp
+++ b/Device/Histo/Histogram2D.cpp
@@ -18,16 +18,15 @@
 #include <memory>
 
 Histogram2D::Histogram2D(int nbinsx, double xlow, double xup, int nbinsy, double ylow, double yup)
+    : IHistogram(FixedBinAxis("x-axis", nbinsx, xlow, xup),
+                 FixedBinAxis("y-axis", nbinsy, ylow, yup))
 {
-    m_data.addAxis(FixedBinAxis("x-axis", nbinsx, xlow, xup));
-    m_data.addAxis(FixedBinAxis("y-axis", nbinsy, ylow, yup));
 }
 
 Histogram2D::Histogram2D(int nbinsx, const std::vector<double>& xbins, int nbinsy,
                          const std::vector<double>& ybins)
+    : IHistogram(VariableBinAxis("x-axis", nbinsx, xbins), VariableBinAxis("y-axis", nbinsy, ybins))
 {
-    m_data.addAxis(VariableBinAxis("x-axis", nbinsx, xbins));
-    m_data.addAxis(VariableBinAxis("y-axis", nbinsy, ybins));
 }
 
 Histogram2D::Histogram2D(const IAxis& axis_x, const IAxis& axis_y)
@@ -40,27 +39,6 @@ Histogram2D::Histogram2D(const Powerfield<double>& data)
     init_from_data(data);
 }
 
-// IMPORTANT intentionally passed by copy to avoid problems on Python side
-Histogram2D::Histogram2D(std::vector<std::vector<double>> data)
-{
-    initFromShape(data);
-    this->setContent(data);
-}
-
-void Histogram2D::initFromShape(const std::vector<std::vector<double>>& data)
-{
-    auto shape = DataUtils::Array::getShape(data);
-    const size_t nrows = shape.first;
-    const size_t ncols = shape.second;
-
-    if (nrows == 0 || ncols == 0)
-        throw std::runtime_error("Histogram2D::Histogram2D() -> Error. "
-                                 "Not a two-dimensional numpy array");
-
-    m_data.addAxis(FixedBinAxis("x-axis", ncols, 0.0, static_cast<double>(ncols)));
-    m_data.addAxis(FixedBinAxis("y-axis", nrows, 0.0, static_cast<double>(nrows)));
-}
-
 Histogram2D* Histogram2D::clone() const
 {
     return new Histogram2D(*this);
diff --git a/Device/Histo/Histogram2D.h b/Device/Histo/Histogram2D.h
index cb6dade65ef2b0d421a2456db147d9f6183ae64f..aeea4fb8fe79e9893688cf57681b27fca2b38887 100644
--- a/Device/Histo/Histogram2D.h
+++ b/Device/Histo/Histogram2D.h
@@ -49,7 +49,7 @@ public:
     Histogram2D(const Powerfield<double>& data);
 
     //! Constructor for 2D histograms from numpy array (thanks to swig)
-    Histogram2D(std::vector<std::vector<double>> data);
+    //    Histogram2D(std::vector<std::vector<double>> data);
 
     //! Returns clone of other histogram
     Histogram2D* clone() const override;
@@ -100,8 +100,6 @@ public:
     void addContent(const std::vector<std::vector<double>>& data);
 
 protected:
-    void initFromShape(const std::vector<std::vector<double>>& data);
-
     //! Creates projection along X. The projections is made by collecting the data in the range
     //! between [ybinlow, ybinup].
     Histogram1D* create_projectionX(int ybinlow, int ybinup);
diff --git a/Device/Histo/IHistogram.cpp b/Device/Histo/IHistogram.cpp
index af2a1b2a735b018fbc5c61eaaedc1b3dba8ad812..dd8ab28a25426aa5780f16e74520db353d0d1de4 100644
--- a/Device/Histo/IHistogram.cpp
+++ b/Device/Histo/IHistogram.cpp
@@ -26,14 +26,13 @@ IHistogram::IHistogram(const IHistogram& other)
 }
 
 IHistogram::IHistogram(const IAxis& axis_x)
+    : m_data(axis_x)
 {
-    m_data.addAxis(axis_x);
 }
 
 IHistogram::IHistogram(const IAxis& axis_x, const IAxis& axis_y)
+    : m_data(axis_x, axis_y)
 {
-    m_data.addAxis(axis_x);
-    m_data.addAxis(axis_y);
 }
 
 size_t IHistogram::getTotalNumberOfBins() const
@@ -246,16 +245,6 @@ IHistogram* IHistogram::createHistogram(const Powerfield<double>& source)
     throw std::runtime_error(message.str());
 }
 
-IHistogram* IHistogram::createFrom(const std::string& filename)
-{
-    return IOFactory::readIntensityData(filename);
-}
-
-IHistogram* IHistogram::createFrom(const std::vector<std::vector<double>>& data)
-{
-    return new Histogram2D(data);
-}
-
 void IHistogram::check_x_axis() const
 {
     if (rank() < 1) {
diff --git a/Device/Histo/IHistogram.h b/Device/Histo/IHistogram.h
index d9da9fff8340023f65e2fd5a4b769283426f1f50..9d8fbd247b4ef8208168510d74df3ae904f75dd7 100644
--- a/Device/Histo/IHistogram.h
+++ b/Device/Histo/IHistogram.h
@@ -149,12 +149,6 @@ public:
 
     static IHistogram* createHistogram(const Powerfield<double>& source);
 
-    //! Creates new histogram from file content.
-    static IHistogram* createFrom(const std::string& filename);
-
-    //! Create new histogram from numpy array.
-    static IHistogram* createFrom(const std::vector<std::vector<double>>& data);
-
     //! Creates new Powerfield with histogram's shape and values corresponding to DataType.
     Powerfield<double>* createPowerfield(DataType dataType = DataType::INTEGRAL) const;
 
diff --git a/Device/Histo/SimulationResult.cpp b/Device/Histo/SimulationResult.cpp
index 882390673de8a01b78f0eed3b53269bf925858b9..e7dff16840ff21fc1693879941362f0850d7f549 100644
--- a/Device/Histo/SimulationResult.cpp
+++ b/Device/Histo/SimulationResult.cpp
@@ -106,11 +106,6 @@ size_t SimulationResult::size() const
     return m_data ? m_data->getAllocatedSize() : 0;
 }
 
-double SimulationResult::max() const
-{
-    return m_data->max();
-}
-
 SimulationResult SimulationResult::relativeToMaximum() const
 {
     std::unique_ptr<Powerfield<double>> data2{m_data->clone()};
diff --git a/Device/Histo/SimulationResult.h b/Device/Histo/SimulationResult.h
index 44e0d8e9e043940c5d1036724d2e36985e5a281f..9b5ac566f847337551ac8b640197354ca59b92b0 100644
--- a/Device/Histo/SimulationResult.h
+++ b/Device/Histo/SimulationResult.h
@@ -67,7 +67,6 @@ public:
     double& operator[](size_t i);
     const double& operator[](size_t i) const;
     size_t size() const;
-    double max() const;
     bool empty() const { return size() == 0; }
 
     //! Returns modified SimulationResult: all intensities dvided by maximum intensity
diff --git a/Device/Mask/DetectorMask.cpp b/Device/Mask/DetectorMask.cpp
index c90daaa88b789d7e0eb32f85fc55f1a1039500ef..de536bcfdd2f64225ca2a341850e4a1be036fb33 100644
--- a/Device/Mask/DetectorMask.cpp
+++ b/Device/Mask/DetectorMask.cpp
@@ -22,6 +22,17 @@ DetectorMask::DetectorMask()
 {
 }
 
+DetectorMask::DetectorMask(const IAxis& xAxis, const IAxis& yAxis)
+{
+    ASSERT(m_shapes.size() == m_mask_of_shape.size());
+    m_mask_data.clear();
+
+    m_mask_data.addAxis(xAxis);
+    m_mask_data.addAxis(yAxis);
+
+    process_masks();
+}
+
 DetectorMask::~DetectorMask() = default;
 
 DetectorMask::DetectorMask(const DetectorMask& other)
diff --git a/Device/Mask/DetectorMask.h b/Device/Mask/DetectorMask.h
index faf1ff80661c5814b7caf8cdb101ad1e54cbec0a..ea63d292a47c5f73116cf72bc9f14d16c2ccfcfa 100644
--- a/Device/Mask/DetectorMask.h
+++ b/Device/Mask/DetectorMask.h
@@ -33,6 +33,7 @@ class Powerfield;
 class DetectorMask {
 public:
     DetectorMask();
+    DetectorMask(const IAxis& xAxis, const IAxis& yAxis);
     ~DetectorMask();
     DetectorMask(const DetectorMask& other);
     DetectorMask& operator=(const DetectorMask& other);
diff --git a/GUI/Model/Data/DataProperties.cpp b/GUI/Model/Data/DataProperties.cpp
index 593c3608fe921dadbf47e163945e0771e85d56ec..bbbe44a4e4296c5112a28a3920846e49400266ff 100644
--- a/GUI/Model/Data/DataProperties.cpp
+++ b/GUI/Model/Data/DataProperties.cpp
@@ -38,8 +38,7 @@ const QMap<QString, QCPScatterStyle::ScatterShape> scatter_map = {
     {"Circle", QCPScatterStyle::ScatterShape::ssCircle},
     {"Cross", QCPScatterStyle::ScatterShape::ssCross},
     {"Diamond", QCPScatterStyle::ScatterShape::ssDiamond},
-    {"Star", QCPScatterStyle::ScatterShape::ssStar}
-};
+    {"Star", QCPScatterStyle::ScatterShape::ssStar}};
 
 // connection lines for representation of 1D graphs
 const QMap<QString, QCPGraph::LineStyle> line_map = {
@@ -48,8 +47,7 @@ const QMap<QString, QCPGraph::LineStyle> line_map = {
     {"StepLeft", QCPGraph::LineStyle::lsStepLeft},
     {"StepRight", QCPGraph::LineStyle::lsStepRight},
     {"StepCenter", QCPGraph::LineStyle::lsStepCenter},
-    {"Impulse", QCPGraph::LineStyle::lsImpulse}
-};
+    {"Impulse", QCPGraph::LineStyle::lsImpulse}};
 
 struct ColorNameComparator {
     ColorNameComparator(QString value_to_comp)
diff --git a/GUI/Model/Data/DataProperties.h b/GUI/Model/Data/DataProperties.h
index 9818b55d05c9db2efcffbde866a393b4063bf806..b64b9d0822dfe2fbbe9ffd942576da88d73d1871 100644
--- a/GUI/Model/Data/DataProperties.h
+++ b/GUI/Model/Data/DataProperties.h
@@ -15,8 +15,8 @@
 #ifndef BORNAGAIN_GUI_MODEL_DATA_DATAPROPERTIES_H
 #define BORNAGAIN_GUI_MODEL_DATA_DATAPROPERTIES_H
 
-#include "GUI/Model/BaseItem/SessionItem.h"
 #include "3rdparty/GUI/qcustomplot/qcustomplot.h"
+#include "GUI/Model/BaseItem/SessionItem.h"
 
 class DataItem;
 
diff --git a/GUI/Model/Data/DataPropertyContainer.cpp b/GUI/Model/Data/DataPropertyContainer.cpp
index ccf9b5784a82a91e8b0f6961187f1b0c3e0229f0..be4762481a8a947471379599131c2086fcfeaaac 100644
--- a/GUI/Model/Data/DataPropertyContainer.cpp
+++ b/GUI/Model/Data/DataPropertyContainer.cpp
@@ -14,8 +14,8 @@
 
 #include "GUI/Model/Data/DataPropertyContainer.h"
 #include "GUI/Model/Data/DataItem.h"
-#include "GUI/Model/Data/RealDataItem.h"
 #include "GUI/Model/Data/DataProperties.h"
+#include "GUI/Model/Data/RealDataItem.h"
 #include "GUI/Util/Error.h"
 
 DataPropertyContainer::DataPropertyContainer()
@@ -66,7 +66,7 @@ void DataPropertyContainer::addItem(DataItem* data_item)
     property_item->setDataItem(data_item);
     property_item->setColorProperty(Data1DProperties::nextColorName(previous_item));
 
-    if(dynamic_cast<RealDataItem*>(data_item->parent())) {
+    if (dynamic_cast<RealDataItem*>(data_item->parent())) {
         property_item->setScatterProperty("Disc");
         property_item->setLineProperty("None");
     } else {
diff --git a/GUI/View/Mask/MaskResultsPresenter.cpp b/GUI/View/Mask/MaskResultsPresenter.cpp
index 05582401b3924eafbd07c2bb2ea1b6efd12f82b9..0ac5c4cebdd0a4426a9e95af963dfbc9656905b4 100644
--- a/GUI/View/Mask/MaskResultsPresenter.cpp
+++ b/GUI/View/Mask/MaskResultsPresenter.cpp
@@ -88,7 +88,8 @@ Powerfield<double>* MaskResultsPresenter::createMaskPresentation() const
 {
     // Requesting mask information
     std::unique_ptr<IShape2D> roi;
-    DetectorMask detectorMask;
+    Powerfield<double>* result = m_intensityDataItem->getPowerfield()->clone();
+    DetectorMask detectorMask(result->axis(0), result->axis(1));
     const double scale = 1.0;
     for (int i_row = m_maskModel->rowCount(m_maskContainerIndex); i_row > 0; --i_row) {
         QModelIndex itemIndex = m_maskModel->index(i_row - 1, 0, m_maskContainerIndex);
@@ -109,9 +110,6 @@ Powerfield<double>* MaskResultsPresenter::createMaskPresentation() const
     if (!detectorMask.hasMasks())
         return nullptr;
 
-    Powerfield<double>* result = m_intensityDataItem->getPowerfield()->clone();
-    detectorMask.initMaskData(result->axis(0), result->axis(1));
-
     for (size_t i = 0; i < result->getAllocatedSize(); ++i)
         if (detectorMask.isMasked(i))
             (*result)[i] = 0.0;
diff --git a/GUI/View/PlotSpecular/SpecularPlot.cpp b/GUI/View/PlotSpecular/SpecularPlot.cpp
index d30152b44cdb769ecf5830226f70c152e9f4cda0..11f98cf5967803f97625b13f206ff1a641748b39 100644
--- a/GUI/View/PlotSpecular/SpecularPlot.cpp
+++ b/GUI/View/PlotSpecular/SpecularPlot.cpp
@@ -12,10 +12,10 @@
 //
 //  ************************************************************************************************
 
+#include "GUI/View/PlotSpecular/SpecularPlot.h"
+#include "GUI/Model/Data/RealDataItem.h"
 #include "GUI/Model/Data/SpecularDataItem.h"
 #include "GUI/Model/Device/AxesItems.h"
-#include "GUI/Model/Data/RealDataItem.h"
-#include "GUI/View/PlotSpecular/SpecularPlot.h"
 #include "GUI/View/PlotUtil/PlotConstants.h"
 #include "GUI/View/PlotUtil/PlotEventInfo.h"
 #include "GUI/View/PlotUtil/RangeUtils.h"
@@ -142,12 +142,12 @@ void SpecularPlot::initPlot()
     m_custom_plot->xAxis->setTickLabelFont(
         QFont(QFont().family(), GUI::Constants::plot_tick_label_size()));
     m_custom_plot->yAxis->setTickLabelFont(
-                QFont(QFont().family(), GUI::Constants::plot_tick_label_size()));
+        QFont(QFont().family(), GUI::Constants::plot_tick_label_size()));
 }
 
 void SpecularPlot::initScatter()
 {
-    if(dynamic_cast<RealDataItem*>(specularItem()->parent())) {
+    if (dynamic_cast<RealDataItem*>(specularItem()->parent())) {
         m_custom_plot->graph()->setScatterStyle(QCPScatterStyle::ssDisc);
     } else {
         m_custom_plot->graph()->setScatterStyle(QCPScatterStyle::ssNone);
diff --git a/Tests/PyUnit/histogram2d.py b/Tests/PyUnit/histogram2d.py
deleted file mode 100644
index ae61004816210060c1183538e6566538ea850bc8..0000000000000000000000000000000000000000
--- a/Tests/PyUnit/histogram2d.py
+++ /dev/null
@@ -1,126 +0,0 @@
-import unittest
-import numpy
-import bornagain as ba
-
-
-class Histogram2DTest(unittest.TestCase):
-
-    def test_constructFromNumpyInt(self):
-        """
-        Testing construction of 2D histogram from numpy array.
-        Automatic conversion under Python3 is not working, that's why it is float64 and not int64
-        """
-        arr = numpy.array(
-            [[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15]],
-            numpy.float64)
-        hist = ba.Histogram2D(arr)
-
-        self.assertEqual(hist.getNbinsX(), 5)
-        self.assertEqual(hist.xMin(), 0)
-        self.assertEqual(hist.xMax(), 5)
-
-        self.assertEqual(hist.getNbinsY(), 3)
-        self.assertEqual(hist.yMin(), 0)
-        self.assertEqual(hist.yMax(), 3)
-
-        self.assertEqual(hist.binContent(0, 0), 11)
-        self.assertEqual(hist.binContent(0, 1), 6)
-        self.assertEqual(hist.binContent(4, 2), 5)
-
-        arr_from_hist = hist.array()
-
-        for (x, y), element in numpy.ndenumerate(arr):
-            self.assertEqual(element, arr_from_hist[x][y])
-
-    def test_constructFromNumpyDouble(self):
-        """
-        Testing construction of 2D histogram from numpy array
-        """
-        arr = numpy.array([[1, 2, 3, 4, 5.0], [6, 7, 8, 9, 10.0],
-                           [11, 12, 13, 14, 15.0]],
-                          dtype=numpy.float64)
-        hist = ba.Histogram2D(arr)
-
-        self.assertEqual(hist.getNbinsX(), 5)
-        self.assertEqual(hist.xMin(), 0)
-        self.assertEqual(hist.xMax(), 5)
-
-        self.assertEqual(hist.getNbinsY(), 3)
-        self.assertEqual(hist.yMin(), 0)
-        self.assertEqual(hist.yMax(), 3)
-
-        self.assertEqual(hist.binContent(0, 0), 11)
-        self.assertEqual(hist.binContent(0, 1), 6)
-        self.assertEqual(hist.binContent(4, 2), 5)
-
-        arr_from_hist = hist.array()
-
-        for (x, y), element in numpy.ndenumerate(arr):
-            self.assertEqual(element, arr_from_hist[x][y])
-
-    def test_constructAndAddFromNumpyInt(self):
-        """
-        Adding to the histogram content from numpy array
-        """
-        arr = numpy.array(
-            [[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15]],
-            dtype=numpy.float64)
-        print(type(arr))
-        hist = ba.Histogram2D(arr)
-        # adding same content once again
-        hist.addContent(arr)
-        arr_from_hist = hist.array()
-
-        for (x, y), element in numpy.ndenumerate(arr):
-            self.assertEqual(element*2, arr_from_hist[x][y])
-
-    def test_constructAndAddFromNumpyDouble(self):
-        """
-        Adding to the histogram content from numpy array
-        """
-        arr = numpy.array([[1, 2, 3, 4, 5.0], [6, 7, 8, 9, 10.0],
-                           [11, 12, 13, 14, 15.0]],
-                          dtype=numpy.float64)
-        hist = ba.Histogram2D(arr)
-        # adding same content once again
-        hist.addContent(arr)
-        arr_from_hist = hist.array()
-
-        for (x, y), element in numpy.ndenumerate(arr):
-            self.assertEqual(element*2, arr_from_hist[x][y])
-
-    def create_histogram(self, arr):
-        """
-        Returns newly created object
-        """
-        return ba.IHistogram.createFrom(arr)
-
-    def test_createFromInt(self):
-        """
-        Testing newly create object
-        """
-        arr = numpy.array(
-            [[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15]],
-            dtype=numpy.float64)
-        hist = self.create_histogram(arr)
-        arr_from_hist = hist.array()
-
-        for (x, y), element in numpy.ndenumerate(arr):
-            self.assertEqual(element, arr_from_hist[x][y])
-
-    def test_createFromDouble(self):
-        """
-        Testing newly create object
-        """
-        arr = numpy.array([[1, 2, 3, 4, 5.0], [6, 7, 8, 9, 10.0],
-                           [11, 12, 13, 14, 15.0]],
-                          dtype=numpy.float64)
-        hist = self.create_histogram(arr)
-        arr_from_hist = hist.array()
-
-        for (x, y), element in numpy.ndenumerate(arr):
-            self.assertEqual(element, arr_from_hist[x][y])
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/Tests/Unit/Device/PowerfieldIteratorTest.cpp b/Tests/Unit/Device/PowerfieldIteratorTest.cpp
index 02e81f62e478b8d46bd496494737862a0bc71e20..777580e4b1e3bab602db8312e0262e75b04663cb 100644
--- a/Tests/Unit/Device/PowerfieldIteratorTest.cpp
+++ b/Tests/Unit/Device/PowerfieldIteratorTest.cpp
@@ -13,25 +13,19 @@ PowerfieldIteratorTest::PowerfieldIteratorTest()
     auto* dims = new int[2];
     dims[0] = 3;
     dims[1] = 5;
-    _data.setAxisSizes(2, dims);
-    double value = 0.0;
-    Powerfield<double>::iterator it = _data.begin();
-    while (it != _data.end()) {
-        *it = value;
-        value += 1.0;
-        ++it;
-    }
-    delete[] dims;
+    _data.addAxis("x", 3, 0, 1);
+    _data.addAxis("y", 5, 0, 1);
+    _data.setRawDataVector({1., 2., 3., 4., 5., 11., 12., 13., 14., 15., 21., 22., 23., 24., 25.});
 }
 
 TEST_F(PowerfieldIteratorTest, Iterate)
 {
     Powerfield<double>::iterator it = _data.begin();
-    EXPECT_EQ(0.0, *it);
+    EXPECT_EQ(1.0, *it);
     for (size_t i = 0; i < 14; ++i) {
         ++it;
     }
-    EXPECT_DOUBLE_EQ(14.0, *it);
+    EXPECT_DOUBLE_EQ(25.0, *it);
     ++it;
     EXPECT_EQ(it, _data.end());
     ++it;
@@ -41,11 +35,11 @@ TEST_F(PowerfieldIteratorTest, Iterate)
 TEST_F(PowerfieldIteratorTest, ConstIterate)
 {
     Powerfield<double>::const_iterator it = _data.begin();
-    EXPECT_EQ(0.0, *it);
+    EXPECT_EQ(1.0, *it);
     for (size_t i = 0; i < 14; ++i) {
         ++it;
     }
-    EXPECT_DOUBLE_EQ(14.0, *it);
+    EXPECT_DOUBLE_EQ(25.0, *it);
     ++it;
     EXPECT_EQ(it, _data.end());
     ++it;
diff --git a/Tests/Unit/Device/PowerfieldTest.cpp b/Tests/Unit/Device/PowerfieldTest.cpp
index c40c121e775b35a0eda2315634b51f743fc4c090..14f929eee55fa983f27bffe9af20da6578fef0f1 100644
--- a/Tests/Unit/Device/PowerfieldTest.cpp
+++ b/Tests/Unit/Device/PowerfieldTest.cpp
@@ -90,25 +90,6 @@ TEST_F(PowerfieldTest, DataCopying)
     EXPECT_TRUE(data2.hasSameShape(data2));
 }
 
-TEST_F(PowerfieldTest, MaxElement)
-{
-    Powerfield<double> data;
-    data.addAxis("axis1", 10, 0., 10.);
-    data.addAxis("axis2", 2, 0., 10.);
-
-    Powerfield<double>::iterator it = data.begin();
-    for (size_t i = 0; i < data.getAllocatedSize(); ++i)
-        if (i == 10)
-            (*it) = 10.0;
-
-    Powerfield<double>::const_iterator cit = std::max_element(data.begin(), data.end());
-    EXPECT_EQ(10., (*cit));
-
-    EXPECT_EQ(10., data.max());
-
-    EXPECT_EQ(1., data.normalizedToMaximum().max());
-}
-
 // y |
 // --------------------------------------------
 // 1 | 1   3   5   7   9   11  13  15  17  19 |
diff --git a/Tests/Unit/Device/SimulationResultTest.cpp b/Tests/Unit/Device/SimulationResultTest.cpp
deleted file mode 100644
index 3e13a47a95aee701ee4dbfef7ce2a07b3cad6052..0000000000000000000000000000000000000000
--- a/Tests/Unit/Device/SimulationResultTest.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#include "Device/Histo/SimulationResult.h"
-#include "Base/Const/Units.h"
-#include "Device/Beam/Beam.h"
-#include "Device/Coord/CoordSystem2D.h"
-#include "Device/Data/Powerfield.h"
-#include "Device/Detector/SphericalDetector.h"
-#include "Tests/GTestWrapper/google_test.h"
-
-class SimulationResultTest : public ::testing::Test {
-};
-
-TEST_F(SimulationResultTest, max)
-{
-    Powerfield<double> data;
-    data.addAxis("axis1", 10, 0., 10.);
-    data.addAxis("axis2", 2, 0., 10.);
-    Powerfield<double>::iterator it = data.begin();
-    for (size_t i = 0; i < data.getAllocatedSize(); ++i)
-        if (i == 10)
-            (*it) = 10.0;
-
-    SphericalDetector det(100, 0.0, 5.0 * Units::deg, 70, -2.0 * Units::deg, 1.5);
-    Beam beam(1.0, 1.0, {1 * Units::deg, 0});
-    SphericalCoords coords(det.axesClippedToRegionOfInterest(), beam.direction(),
-                           beam.wavelength());
-
-    SimulationResult sr1(data, coords.clone());
-
-    SimulationResult sr2 = sr1.relativeToMaximum();
-
-    EXPECT_EQ(sr2.max(), 1.);
-}
diff --git a/Wrap/Swig/libBornAgainDevice.i b/Wrap/Swig/libBornAgainDevice.i
index b297dd35c8b2782ccf988d3918e87ec1850d0bc6..710e4f1fa03d5c4847d74ab802d44f6a6d8b433a 100644
--- a/Wrap/Swig/libBornAgainDevice.i
+++ b/Wrap/Swig/libBornAgainDevice.i
@@ -62,8 +62,6 @@
 %newobject DetectorMask::createHistogram() const;
 
 %newobject DataUtils::Data::importArrayToPowerfield;
-%newobject IHistogram::createFrom(const std::string& filename);
-%newobject IHistogram::createFrom(const std::vector<std::vector<double>>& data);
 
 %include "Device/Data/Powerfield.h"
 %template(IntensityData) Powerfield<double>;
diff --git a/auto/Wrap/doxygenBase.i b/auto/Wrap/doxygenBase.i
index 11697d8f5e728e555d89f9519b9a9758e836b672..2cbc6532cfc0b66f1229045eff79b40d8da95fa3 100644
--- a/auto/Wrap/doxygenBase.i
+++ b/auto/Wrap/doxygenBase.i
@@ -570,10 +570,17 @@ C++ includes: OwningVector.h
 %feature("docstring")  OwningVector::OwningVector "OwningVector< T >::OwningVector()=default
 ";
 
-%feature("docstring")  OwningVector::~OwningVector "OwningVector< T >::~OwningVector()
+%feature("docstring")  OwningVector::OwningVector "OwningVector< T >::OwningVector(const std::vector< T * > &v)
+
+Constructor that takes over ownership of elements in given vector. 
 ";
 
 %feature("docstring")  OwningVector::OwningVector "OwningVector< T >::OwningVector(const OwningVector &other)
+
+Constructor that clones elements in given vector. 
+";
+
+%feature("docstring")  OwningVector::~OwningVector "OwningVector< T >::~OwningVector()
 ";
 
 %feature("docstring")  OwningVector::emplace_back "void OwningVector< T >::emplace_back(T *e)
diff --git a/auto/Wrap/doxygenDevice.i b/auto/Wrap/doxygenDevice.i
index e5af0eded26b530841e1d8c23f503d660028c387..1eb9c324bc3380094bb92510ea70a5332f536a6a 100644
--- a/auto/Wrap/doxygenDevice.i
+++ b/auto/Wrap/doxygenDevice.i
@@ -408,6 +408,9 @@ C++ includes: DetectorMask.h
 %feature("docstring")  DetectorMask::DetectorMask "DetectorMask::DetectorMask()
 ";
 
+%feature("docstring")  DetectorMask::DetectorMask "DetectorMask::DetectorMask(const IAxis &xAxis, const IAxis &yAxis)
+";
+
 %feature("docstring")  DetectorMask::~DetectorMask "DetectorMask::~DetectorMask()
 ";
 
@@ -739,14 +742,11 @@ Constructor for 2D histogram with custom axes.
 Constructor for 2D histograms from basic  Powerfield object. 
 ";
 
-%feature("docstring")  Histogram2D::Histogram2D "Histogram2D::Histogram2D(std::vector< std::vector< double >> data)
-
-Constructor for 2D histograms from numpy array (thanks to swig) 
-";
-
 %feature("docstring")  Histogram2D::clone "Histogram2D * Histogram2D::clone() const override
 
-Returns clone of other histogram. 
+Constructor for 2D histograms from numpy array (thanks to swig)
+
+Returns clone of other histogram 
 ";
 
 %feature("docstring")  Histogram2D::rank "size_t Histogram2D::rank() const override
@@ -1194,6 +1194,182 @@ Apply the resolution function to the intensity data.
 ";
 
 
+// File: classIField.xml
+%feature("docstring") IField "
+
+Holds one or two axes. Base class for  Powerfield.
+
+C++ includes: Powerfield.h
+";
+
+%feature("docstring")  IField::IField "IField::IField()=default
+";
+
+%feature("docstring")  IField::IField "IField::IField(const std::vector< IAxis * > &axes)
+";
+
+%feature("docstring")  IField::~IField "virtual IField::~IField()
+";
+
+%feature("docstring")  IField::rank "size_t IField::rank() const
+
+Returns number of dimensions. 
+";
+
+%feature("docstring")  IField::axis "const IAxis & IField::axis(size_t serial_number) const
+
+Returns axis with given serial number. 
+";
+
+%feature("docstring")  IField::getAxisValue "double IField::getAxisValue(size_t global_index, size_t i_selected_axis) const
+
+Returns the value of selected axis for given global_index.
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+i_selected_axis: 
+Serial number of selected axis.
+
+corresponding bin center of selected axis 
+";
+
+%feature("docstring")  IField::getAxisValue "double IField::getAxisValue(size_t global_index, const std::string &axis_name) const
+
+Returns the value of selected axis for given global_index.
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+axis_name: 
+The name of selected axis.
+
+corresponding bin center of selected axis 
+";
+
+%feature("docstring")  IField::getAxesValues "std::vector< double > IField::getAxesValues(size_t global_index) const
+
+Returns values on all defined axes for given globalbin number
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+Vector of corresponding bin centers 
+";
+
+%feature("docstring")  IField::getAxisBin "Bin1D IField::getAxisBin(size_t global_index, size_t i_selected_axis) const
+
+Returns bin of selected axis for given global_index.
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+i_selected_axis: 
+Serial number of selected axis.
+
+Corresponding Bin1D object 
+";
+
+%feature("docstring")  IField::getAxisBin "Bin1D IField::getAxisBin(size_t global_index, const std::string &axis_name) const
+
+Returns bin of selected axis for given global_index.
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+axis_name: 
+The name of selected axis.
+
+Corresponding Bin1D object 
+";
+
+%feature("docstring")  IField::getAxisBinIndex "size_t IField::getAxisBinIndex(size_t global_index, const std::string &axis_name) const
+
+Returns axis bin index for given global index
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+axis_name: 
+The name of selected axis.
+
+Corresponding bin index for selected axis 
+";
+
+%feature("docstring")  IField::getAxesBinIndices "std::vector< int > IField::getAxesBinIndices(size_t global_index) const
+
+Returns vector of axes indices for given global index
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+Vector of bin indices for all axes defined 
+";
+
+%feature("docstring")  IField::getAxisBinIndex "size_t IField::getAxisBinIndex(size_t global_index, size_t i_selected_axis) const
+
+Returns axis bin index for given global index
+
+Parameters:
+-----------
+
+global_index: 
+The global index of this data structure.
+
+i_selected_axis: 
+Serial number of selected axis.
+
+Corresponding bin index for selected axis 
+";
+
+%feature("docstring")  IField::toGlobalIndex "size_t IField::toGlobalIndex(const std::vector< unsigned > &axes_indices) const
+
+Returns global index for specified indices of axes
+
+Parameters:
+-----------
+
+axes_indices: 
+Vector of axes indices for all specified axes in this dataset
+
+Corresponding global index 
+";
+
+%feature("docstring")  IField::findGlobalIndex "size_t IField::findGlobalIndex(const std::vector< double > &coordinates) const
+
+Returns global index for specified axes values
+
+Parameters:
+-----------
+
+coordinates: 
+Vector of axes coordinates for all specified axes in this dataset
+
+Closest global index 
+";
+
+
 // File: classIFootprintFactor.xml
 %feature("docstring") IFootprintFactor "
 
@@ -1854,13 +2030,19 @@ C++ includes: Powerfield.h
 %feature("docstring")  Powerfield::Powerfield "Powerfield< T >::Powerfield()
 ";
 
+%feature("docstring")  Powerfield::Powerfield "Powerfield< T >::Powerfield(const IAxis &xAxis)
+";
+
+%feature("docstring")  Powerfield::Powerfield "Powerfield< T >::Powerfield(const IAxis &xAxis, const IAxis &yAxis)
+";
+
 %feature("docstring")  Powerfield::Powerfield "Powerfield< T >::Powerfield(const Powerfield &)=delete
 ";
 
 %feature("docstring")  Powerfield::Powerfield "Powerfield< T >::Powerfield(Powerfield &&)
 ";
 
-%feature("docstring")  Powerfield::~Powerfield "Powerfield< T >::~Powerfield()
+%feature("docstring")  Powerfield::~Powerfield "Powerfield< T >::~Powerfield() override
 ";
 
 %feature("docstring")  Powerfield::clone "Powerfield< T > * Powerfield< T >::clone() const
@@ -1881,16 +2063,6 @@ C++ includes: Powerfield.h
 %feature("docstring")  Powerfield::addAxis "void Powerfield< T >::addAxis(const std::string &name, size_t size, double start, double end)
 ";
 
-%feature("docstring")  Powerfield::axis "const IAxis & Powerfield< T >::axis(size_t serial_number) const
-
-Returns axis with given serial number. 
-";
-
-%feature("docstring")  Powerfield::rank "size_t Powerfield< T >::rank() const
-
-Returns number of dimensions. 
-";
-
 %feature("docstring")  Powerfield::getAllocatedSize "size_t Powerfield< T >::getAllocatedSize() const
 
 Returns total size of data buffer (product of bin number in every dimension). 
@@ -1911,11 +2083,6 @@ Returns copy of raw data vector.
 Returns sum of all values in the data structure. 
 ";
 
-%feature("docstring")  Powerfield::max "T Powerfield< T >::max() const
-
-Returns maximum value in the data structure. 
-";
-
 %feature("docstring")  Powerfield::begin "Powerfield< T >::iterator Powerfield< T >::begin()
 
 Returns read/write iterator that points to the first element. 
@@ -1936,154 +2103,6 @@ Returns read/write iterator that points to the one past last element.
 Returns read-only iterator that points to the one past last element. 
 ";
 
-%feature("docstring")  Powerfield::getAxesBinIndices "std::vector< int > Powerfield< T >::getAxesBinIndices(size_t global_index) const
-
-Returns vector of axes indices for given global index
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-Vector of bin indices for all axes defined 
-";
-
-%feature("docstring")  Powerfield::getAxisBinIndex "size_t Powerfield< T >::getAxisBinIndex(size_t global_index, size_t i_selected_axis) const
-
-Returns axis bin index for given global index
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-i_selected_axis: 
-Serial number of selected axis.
-
-Corresponding bin index for selected axis 
-";
-
-%feature("docstring")  Powerfield::getAxisBinIndex "size_t Powerfield< T >::getAxisBinIndex(size_t global_index, const std::string &axis_name) const
-
-Returns axis bin index for given global index
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-axis_name: 
-The name of selected axis.
-
-Corresponding bin index for selected axis 
-";
-
-%feature("docstring")  Powerfield::toGlobalIndex "size_t Powerfield< T >::toGlobalIndex(const std::vector< unsigned > &axes_indices) const
-
-Returns global index for specified indices of axes
-
-Parameters:
------------
-
-axes_indices: 
-Vector of axes indices for all specified axes in this dataset
-
-Corresponding global index 
-";
-
-%feature("docstring")  Powerfield::findGlobalIndex "size_t Powerfield< T >::findGlobalIndex(const std::vector< double > &coordinates) const
-
-Returns global index for specified axes values
-
-Parameters:
------------
-
-coordinates: 
-Vector of axes coordinates for all specified axes in this dataset
-
-Closest global index 
-";
-
-%feature("docstring")  Powerfield::getAxisValue "double Powerfield< T >::getAxisValue(size_t global_index, size_t i_selected_axis) const
-
-Returns the value of selected axis for given global_index.
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-i_selected_axis: 
-Serial number of selected axis.
-
-corresponding bin center of selected axis 
-";
-
-%feature("docstring")  Powerfield::getAxisValue "double Powerfield< T >::getAxisValue(size_t global_index, const std::string &axis_name) const
-
-Returns the value of selected axis for given global_index.
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-axis_name: 
-The name of selected axis.
-
-corresponding bin center of selected axis 
-";
-
-%feature("docstring")  Powerfield::getAxesValues "std::vector< double > Powerfield< T >::getAxesValues(size_t global_index) const
-
-Returns values on all defined axes for given globalbin number
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-Vector of corresponding bin centers 
-";
-
-%feature("docstring")  Powerfield::getAxisBin "Bin1D Powerfield< T >::getAxisBin(size_t global_index, size_t i_selected_axis) const
-
-Returns bin of selected axis for given global_index.
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-i_selected_axis: 
-Serial number of selected axis.
-
-Corresponding Bin1D object 
-";
-
-%feature("docstring")  Powerfield::getAxisBin "Bin1D Powerfield< T >::getAxisBin(size_t global_index, const std::string &axis_name) const
-
-Returns bin of selected axis for given global_index.
-
-Parameters:
------------
-
-global_index: 
-The global index of this data structure.
-
-axis_name: 
-The name of selected axis.
-
-Corresponding Bin1D object 
-";
-
 %feature("docstring")  Powerfield::clear "void Powerfield< T >::clear()
 
 Sets object into initial state (no dimensions, data) 
@@ -2094,11 +2113,6 @@ Sets object into initial state (no dimensions, data)
 Sets content of output data to specific value. 
 ";
 
-%feature("docstring")  Powerfield::setAxisSizes "void Powerfield< T >::setAxisSizes(size_t rank, int *n_dims)
-
-Adds 'rank' axes with indicated sizes. 
-";
-
 %feature("docstring")  Powerfield::setRawDataVector "void Powerfield< T >::setRawDataVector(const std::vector< T > &data_vector)
 
 Sets new values to raw data vector. 
@@ -2114,14 +2128,14 @@ Sets new values to raw data array.
 Returns value or summed value, depending on T. 
 ";
 
-%feature("docstring")  Powerfield::hasSameDimensions "bool Powerfield< T >::hasSameDimensions(const Powerfield< U > &right) const
+%feature("docstring")  Powerfield::hasSameDimensions "bool Powerfield< T >::hasSameDimensions(const Powerfield< U > &other) const
 
 Returns true if object have same dimensions and number of axes bins.
 
 Returns true if object have same dimensions. 
 ";
 
-%feature("docstring")  Powerfield::hasSameShape "bool Powerfield< T >::hasSameShape(const Powerfield< U > &right) const
+%feature("docstring")  Powerfield::hasSameShape "bool Powerfield< T >::hasSameShape(const Powerfield< U > &other) const
 
 Returns true if objects a) have same dimensions b) bin boundaries of axes coincide.
 
@@ -2516,9 +2530,6 @@ Returns underlying unit converter.
 %feature("docstring")  SimulationResult::size "size_t SimulationResult::size() const
 ";
 
-%feature("docstring")  SimulationResult::max "double SimulationResult::max() const
-";
-
 %feature("docstring")  SimulationResult::empty "bool SimulationResult::empty() const
 ";
 
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index d438f0f7c93fab151445bfd0d29e26fb1982bd36..707ea5b481273471d7fc975073e8ebd383328fe8 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2048,11 +2048,11 @@ class vector_R3(object):
 _libBornAgainDevice.vector_R3_swigregister(vector_R3)
 
 import libBornAgainParam
-class IntensityData(object):
+class IField(object):
     r"""
 
 
-    Templated class to store data in multi-dimensional space. Used with type T = double,  CumulativeValue, bool
+    Holds one or two axes. Base class for  Powerfield.
 
     C++ includes: Powerfield.h
 
@@ -2063,143 +2063,100 @@ class IntensityData(object):
 
     def __init__(self, *args):
         r"""
-        __init__(IntensityData self) -> IntensityData
-        __init__(IntensityData self, IntensityData arg2) -> IntensityData
-        Powerfield< T >::Powerfield(Powerfield &&)
-
-        """
-        _libBornAgainDevice.IntensityData_swiginit(self, _libBornAgainDevice.new_IntensityData(*args))
-    __swig_destroy__ = _libBornAgainDevice.delete_IntensityData
-
-    def clone(self):
-        r"""
-        clone(IntensityData self) -> IntensityData
-        Powerfield< T > * Powerfield< T >::clone() const
+        __init__(IField self) -> IField
+        __init__(IField self, std::vector< IAxis *,std::allocator< IAxis * > > const & axes) -> IField
+        IField::IField(const std::vector< IAxis * > &axes)
 
         """
-        return _libBornAgainDevice.IntensityData_clone(self)
+        _libBornAgainDevice.IField_swiginit(self, _libBornAgainDevice.new_IField(*args))
+    __swig_destroy__ = _libBornAgainDevice.delete_IField
 
-    def copyFrom(self, other):
-        r"""
-        copyFrom(IntensityData self, IntensityData other)
-        void Powerfield< T >::copyFrom(const Powerfield< T > &other)
-
-        """
-        return _libBornAgainDevice.IntensityData_copyFrom(self, other)
-
-    def meanValues(self):
+    def rank(self):
         r"""
-        meanValues(IntensityData self) -> IntensityData
-        Powerfield< double > * Powerfield< T >::meanValues() const
+        rank(IField self) -> size_t
+        size_t IField::rank() const
 
-        """
-        return _libBornAgainDevice.IntensityData_meanValues(self)
-
-    def addAxis(self, *args):
-        r"""
-        addAxis(IntensityData self, IAxis new_axis)
-        addAxis(IntensityData self, std::string const & name, size_t size, double start, double end)
-        void Powerfield< T >::addAxis(const std::string &name, size_t size, double start, double end)
+        Returns number of dimensions. 
 
         """
-        return _libBornAgainDevice.IntensityData_addAxis(self, *args)
+        return _libBornAgainDevice.IField_rank(self)
 
     def axis(self, serial_number):
         r"""
-        axis(IntensityData self, size_t serial_number) -> IAxis
-        const IAxis & Powerfield< T >::axis(size_t serial_number) const
+        axis(IField self, size_t serial_number) -> IAxis
+        const IAxis & IField::axis(size_t serial_number) const
 
         Returns axis with given serial number. 
 
         """
-        return _libBornAgainDevice.IntensityData_axis(self, serial_number)
+        return _libBornAgainDevice.IField_axis(self, serial_number)
 
-    def rank(self):
+    def getAxisValue(self, *args):
         r"""
-        rank(IntensityData self) -> size_t
-        size_t Powerfield< T >::rank() const
-
-        Returns number of dimensions. 
-
-        """
-        return _libBornAgainDevice.IntensityData_rank(self)
+        getAxisValue(IField self, size_t global_index, size_t i_selected_axis) -> double
+        getAxisValue(IField self, size_t global_index, std::string const & axis_name) -> double
+        double IField::getAxisValue(size_t global_index, const std::string &axis_name) const
 
-    def getAllocatedSize(self):
-        r"""
-        getAllocatedSize(IntensityData self) -> size_t
-        size_t Powerfield< T >::getAllocatedSize() const
+        Returns the value of selected axis for given global_index.
 
-        Returns total size of data buffer (product of bin number in every dimension). 
+        Parameters:
+        -----------
 
-        """
-        return _libBornAgainDevice.IntensityData_getAllocatedSize(self)
+        global_index: 
+        The global index of this data structure.
 
-    def getAllSizes(self):
-        r"""
-        getAllSizes(IntensityData self) -> std::vector< size_t,std::allocator< size_t > >
-        std::vector< size_t > Powerfield< T >::getAllSizes() const
+        axis_name: 
+        The name of selected axis.
 
-        Returns all sizes of its axes. 
+        corresponding bin center of selected axis 
 
         """
-        return _libBornAgainDevice.IntensityData_getAllSizes(self)
+        return _libBornAgainDevice.IField_getAxisValue(self, *args)
 
-    def getRawDataVector(self):
+    def getAxesValues(self, global_index):
         r"""
-        getRawDataVector(IntensityData self) -> vdouble1d_t
-        std::vector< T > Powerfield< T >::getRawDataVector() const
+        getAxesValues(IField self, size_t global_index) -> vdouble1d_t
+        std::vector< double > IField::getAxesValues(size_t global_index) const
 
-        Returns copy of raw data vector. 
+        Returns values on all defined axes for given globalbin number
 
-        """
-        return _libBornAgainDevice.IntensityData_getRawDataVector(self)
+        Parameters:
+        -----------
 
-    def totalSum(self):
-        r"""
-        totalSum(IntensityData self) -> double
-        T Powerfield< T >::totalSum() const
+        global_index: 
+        The global index of this data structure.
 
-        Returns sum of all values in the data structure. 
+        Vector of corresponding bin centers 
 
         """
-        return _libBornAgainDevice.IntensityData_totalSum(self)
+        return _libBornAgainDevice.IField_getAxesValues(self, global_index)
 
-    def max(self):
+    def getAxisBin(self, *args):
         r"""
-        max(IntensityData self) -> double
-        T Powerfield< T >::max() const
-
-        Returns maximum value in the data structure. 
-
-        """
-        return _libBornAgainDevice.IntensityData_max(self)
+        getAxisBin(IField self, size_t global_index, size_t i_selected_axis) -> Bin1D
+        getAxisBin(IField self, size_t global_index, std::string const & axis_name) -> Bin1D
+        Bin1D IField::getAxisBin(size_t global_index, const std::string &axis_name) const
 
-    def begin(self, *args):
-        r"""
-        begin(IntensityData self) -> Powerfield< double >::iterator
-        begin(IntensityData self) -> Powerfield< double >::const_iterator
-        Powerfield< T >::const_iterator Powerfield< T >::begin() const
+        Returns bin of selected axis for given global_index.
 
-        Returns read-only iterator that points to the first element. 
+        Parameters:
+        -----------
 
-        """
-        return _libBornAgainDevice.IntensityData_begin(self, *args)
+        global_index: 
+        The global index of this data structure.
 
-    def end(self, *args):
-        r"""
-        end(IntensityData self) -> Powerfield< double >::iterator
-        end(IntensityData self) -> Powerfield< double >::const_iterator
-        const_iterator Powerfield< T >::end() const
+        axis_name: 
+        The name of selected axis.
 
-        Returns read-only iterator that points to the one past last element. 
+        Corresponding Bin1D object 
 
         """
-        return _libBornAgainDevice.IntensityData_end(self, *args)
+        return _libBornAgainDevice.IField_getAxisBin(self, *args)
 
     def getAxesBinIndices(self, global_index):
         r"""
-        getAxesBinIndices(IntensityData self, size_t global_index) -> vector_integer_t
-        std::vector< int > Powerfield< T >::getAxesBinIndices(size_t global_index) const
+        getAxesBinIndices(IField self, size_t global_index) -> vector_integer_t
+        std::vector< int > IField::getAxesBinIndices(size_t global_index) const
 
         Returns vector of axes indices for given global index
 
@@ -2212,13 +2169,13 @@ class IntensityData(object):
         Vector of bin indices for all axes defined 
 
         """
-        return _libBornAgainDevice.IntensityData_getAxesBinIndices(self, global_index)
+        return _libBornAgainDevice.IField_getAxesBinIndices(self, global_index)
 
     def getAxisBinIndex(self, *args):
         r"""
-        getAxisBinIndex(IntensityData self, size_t global_index, size_t i_selected_axis) -> size_t
-        getAxisBinIndex(IntensityData self, size_t global_index, std::string const & axis_name) -> size_t
-        size_t Powerfield< T >::getAxisBinIndex(size_t global_index, const std::string &axis_name) const
+        getAxisBinIndex(IField self, size_t global_index, std::string const & axis_name) -> size_t
+        getAxisBinIndex(IField self, size_t global_index, size_t i_selected_axis) -> size_t
+        size_t IField::getAxisBinIndex(size_t global_index, size_t i_selected_axis) const
 
         Returns axis bin index for given global index
 
@@ -2228,18 +2185,18 @@ class IntensityData(object):
         global_index: 
         The global index of this data structure.
 
-        axis_name: 
-        The name of selected axis.
+        i_selected_axis: 
+        Serial number of selected axis.
 
         Corresponding bin index for selected axis 
 
         """
-        return _libBornAgainDevice.IntensityData_getAxisBinIndex(self, *args)
+        return _libBornAgainDevice.IField_getAxisBinIndex(self, *args)
 
     def toGlobalIndex(self, axes_indices):
         r"""
-        toGlobalIndex(IntensityData self, std::vector< unsigned int,std::allocator< unsigned int > > const & axes_indices) -> size_t
-        size_t Powerfield< T >::toGlobalIndex(const std::vector< unsigned > &axes_indices) const
+        toGlobalIndex(IField self, std::vector< unsigned int,std::allocator< unsigned int > > const & axes_indices) -> size_t
+        size_t IField::toGlobalIndex(const std::vector< unsigned > &axes_indices) const
 
         Returns global index for specified indices of axes
 
@@ -2252,12 +2209,12 @@ class IntensityData(object):
         Corresponding global index 
 
         """
-        return _libBornAgainDevice.IntensityData_toGlobalIndex(self, axes_indices)
+        return _libBornAgainDevice.IField_toGlobalIndex(self, axes_indices)
 
     def findGlobalIndex(self, coordinates):
         r"""
-        findGlobalIndex(IntensityData self, vdouble1d_t coordinates) -> size_t
-        size_t Powerfield< T >::findGlobalIndex(const std::vector< double > &coordinates) const
+        findGlobalIndex(IField self, vdouble1d_t coordinates) -> size_t
+        size_t IField::findGlobalIndex(const std::vector< double > &coordinates) const
 
         Returns global index for specified axes values
 
@@ -2270,69 +2227,130 @@ class IntensityData(object):
         Closest global index 
 
         """
-        return _libBornAgainDevice.IntensityData_findGlobalIndex(self, coordinates)
+        return _libBornAgainDevice.IField_findGlobalIndex(self, coordinates)
 
-    def getAxisValue(self, *args):
+# Register IField in _libBornAgainDevice:
+_libBornAgainDevice.IField_swigregister(IField)
+
+class IntensityData(IField):
+    r"""
+
+
+    Templated class to store data in multi-dimensional space. Used with type T = double,  CumulativeValue, bool
+
+    C++ includes: Powerfield.h
+
+    """
+
+    thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
+    __repr__ = _swig_repr
+
+    def __init__(self, *args):
         r"""
-        getAxisValue(IntensityData self, size_t global_index, size_t i_selected_axis) -> double
-        getAxisValue(IntensityData self, size_t global_index, std::string const & axis_name) -> double
-        double Powerfield< T >::getAxisValue(size_t global_index, const std::string &axis_name) const
+        __init__(IntensityData self) -> IntensityData
+        __init__(IntensityData self, IAxis xAxis) -> IntensityData
+        __init__(IntensityData self, IAxis xAxis, IAxis yAxis) -> IntensityData
+        __init__(IntensityData self, IntensityData arg2) -> IntensityData
+        Powerfield< T >::Powerfield(Powerfield &&)
 
-        Returns the value of selected axis for given global_index.
+        """
+        _libBornAgainDevice.IntensityData_swiginit(self, _libBornAgainDevice.new_IntensityData(*args))
+    __swig_destroy__ = _libBornAgainDevice.delete_IntensityData
 
-        Parameters:
-        -----------
+    def clone(self):
+        r"""
+        clone(IntensityData self) -> IntensityData
+        Powerfield< T > * Powerfield< T >::clone() const
 
-        global_index: 
-        The global index of this data structure.
+        """
+        return _libBornAgainDevice.IntensityData_clone(self)
 
-        axis_name: 
-        The name of selected axis.
+    def copyFrom(self, other):
+        r"""
+        copyFrom(IntensityData self, IntensityData other)
+        void Powerfield< T >::copyFrom(const Powerfield< T > &other)
 
-        corresponding bin center of selected axis 
+        """
+        return _libBornAgainDevice.IntensityData_copyFrom(self, other)
+
+    def meanValues(self):
+        r"""
+        meanValues(IntensityData self) -> IntensityData
+        Powerfield< double > * Powerfield< T >::meanValues() const
 
         """
-        return _libBornAgainDevice.IntensityData_getAxisValue(self, *args)
+        return _libBornAgainDevice.IntensityData_meanValues(self)
 
-    def getAxesValues(self, global_index):
+    def addAxis(self, *args):
         r"""
-        getAxesValues(IntensityData self, size_t global_index) -> vdouble1d_t
-        std::vector< double > Powerfield< T >::getAxesValues(size_t global_index) const
+        addAxis(IntensityData self, IAxis new_axis)
+        addAxis(IntensityData self, std::string const & name, size_t size, double start, double end)
+        void Powerfield< T >::addAxis(const std::string &name, size_t size, double start, double end)
 
-        Returns values on all defined axes for given globalbin number
+        """
+        return _libBornAgainDevice.IntensityData_addAxis(self, *args)
 
-        Parameters:
-        -----------
+    def getAllocatedSize(self):
+        r"""
+        getAllocatedSize(IntensityData self) -> size_t
+        size_t Powerfield< T >::getAllocatedSize() const
 
-        global_index: 
-        The global index of this data structure.
+        Returns total size of data buffer (product of bin number in every dimension). 
 
-        Vector of corresponding bin centers 
+        """
+        return _libBornAgainDevice.IntensityData_getAllocatedSize(self)
+
+    def getAllSizes(self):
+        r"""
+        getAllSizes(IntensityData self) -> std::vector< size_t,std::allocator< size_t > >
+        std::vector< size_t > Powerfield< T >::getAllSizes() const
+
+        Returns all sizes of its axes. 
 
         """
-        return _libBornAgainDevice.IntensityData_getAxesValues(self, global_index)
+        return _libBornAgainDevice.IntensityData_getAllSizes(self)
 
-    def getAxisBin(self, *args):
+    def getRawDataVector(self):
         r"""
-        getAxisBin(IntensityData self, size_t global_index, size_t i_selected_axis) -> Bin1D
-        getAxisBin(IntensityData self, size_t global_index, std::string const & axis_name) -> Bin1D
-        Bin1D Powerfield< T >::getAxisBin(size_t global_index, const std::string &axis_name) const
+        getRawDataVector(IntensityData self) -> vdouble1d_t
+        std::vector< T > Powerfield< T >::getRawDataVector() const
 
-        Returns bin of selected axis for given global_index.
+        Returns copy of raw data vector. 
 
-        Parameters:
-        -----------
+        """
+        return _libBornAgainDevice.IntensityData_getRawDataVector(self)
 
-        global_index: 
-        The global index of this data structure.
+    def totalSum(self):
+        r"""
+        totalSum(IntensityData self) -> double
+        T Powerfield< T >::totalSum() const
 
-        axis_name: 
-        The name of selected axis.
+        Returns sum of all values in the data structure. 
 
-        Corresponding Bin1D object 
+        """
+        return _libBornAgainDevice.IntensityData_totalSum(self)
+
+    def begin(self, *args):
+        r"""
+        begin(IntensityData self) -> Powerfield< double >::iterator
+        begin(IntensityData self) -> Powerfield< double >::const_iterator
+        Powerfield< T >::const_iterator Powerfield< T >::begin() const
+
+        Returns read-only iterator that points to the first element. 
 
         """
-        return _libBornAgainDevice.IntensityData_getAxisBin(self, *args)
+        return _libBornAgainDevice.IntensityData_begin(self, *args)
+
+    def end(self, *args):
+        r"""
+        end(IntensityData self) -> Powerfield< double >::iterator
+        end(IntensityData self) -> Powerfield< double >::const_iterator
+        const_iterator Powerfield< T >::end() const
+
+        Returns read-only iterator that points to the one past last element. 
+
+        """
+        return _libBornAgainDevice.IntensityData_end(self, *args)
 
     def clear(self):
         r"""
@@ -2354,16 +2372,6 @@ class IntensityData(object):
         """
         return _libBornAgainDevice.IntensityData_setAllTo(self, value)
 
-    def setAxisSizes(self, rank, n_dims):
-        r"""
-        setAxisSizes(IntensityData self, size_t rank, int * n_dims)
-        void Powerfield< T >::setAxisSizes(size_t rank, int *n_dims)
-
-        Adds 'rank' axes with indicated sizes. 
-
-        """
-        return _libBornAgainDevice.IntensityData_setAxisSizes(self, rank, n_dims)
-
     def setRawDataVector(self, data_vector):
         r"""
         setRawDataVector(IntensityData self, vdouble1d_t data_vector)
@@ -2384,13 +2392,13 @@ class IntensityData(object):
         """
         return _libBornAgainDevice.IntensityData_setRawDataArray(self, source)
 
-    def __iadd__(self, right):
-        r"""__iadd__(IntensityData self, IntensityData right) -> IntensityData"""
-        return _libBornAgainDevice.IntensityData___iadd__(self, right)
+    def __iadd__(self, other):
+        r"""__iadd__(IntensityData self, IntensityData other) -> IntensityData"""
+        return _libBornAgainDevice.IntensityData___iadd__(self, other)
 
-    def __isub__(self, right):
-        r"""__isub__(IntensityData self, IntensityData right) -> IntensityData"""
-        return _libBornAgainDevice.IntensityData___isub__(self, right)
+    def __isub__(self, other):
+        r"""__isub__(IntensityData self, IntensityData other) -> IntensityData"""
+        return _libBornAgainDevice.IntensityData___isub__(self, other)
 
     def __itruediv__(self, *args):
         return _libBornAgainDevice.IntensityData___itruediv__(self, *args)
@@ -2398,9 +2406,9 @@ class IntensityData(object):
 
 
 
-    def __imul__(self, right):
-        r"""__imul__(IntensityData self, IntensityData right) -> IntensityData"""
-        return _libBornAgainDevice.IntensityData___imul__(self, right)
+    def __imul__(self, other):
+        r"""__imul__(IntensityData self, IntensityData other) -> IntensityData"""
+        return _libBornAgainDevice.IntensityData___imul__(self, other)
 
     def getValue(self, index):
         r"""
@@ -3585,6 +3593,7 @@ class DetectorMask(object):
     def __init__(self, *args):
         r"""
         __init__(DetectorMask self) -> DetectorMask
+        __init__(DetectorMask self, IAxis xAxis, IAxis yAxis) -> DetectorMask
         __init__(DetectorMask self, DetectorMask other) -> DetectorMask
         DetectorMask::DetectorMask(const DetectorMask &other)
 
@@ -4902,14 +4911,6 @@ class IHistogram(object):
         r"""createHistogram(IntensityData source) -> IHistogram"""
         return _libBornAgainDevice.IHistogram_createHistogram(source)
 
-    @staticmethod
-    def createFrom(*args):
-        r"""
-        createFrom(std::string const & filename) -> IHistogram
-        createFrom(vdouble2d_t data) -> IHistogram
-        """
-        return _libBornAgainDevice.IHistogram_createFrom(*args)
-
     def createPowerfield(self, *args):
         r"""
         createPowerfield(IHistogram self, IHistogram::DataType dataType=DataType::INTEGRAL) -> IntensityData
@@ -4983,13 +4984,6 @@ def IHistogram_createHistogram(source):
     r"""IHistogram_createHistogram(IntensityData source) -> IHistogram"""
     return _libBornAgainDevice.IHistogram_createHistogram(source)
 
-def IHistogram_createFrom(*args):
-    r"""
-    IHistogram_createFrom(std::string const & filename) -> IHistogram
-    IHistogram_createFrom(vdouble2d_t data) -> IHistogram
-    """
-    return _libBornAgainDevice.IHistogram_createFrom(*args)
-
 class Histogram1D(IHistogram):
     r"""
 
@@ -5133,10 +5127,9 @@ class Histogram2D(IHistogram):
         __init__(Histogram2D self, int nbinsx, vdouble1d_t xbins, int nbinsy, vdouble1d_t ybins) -> Histogram2D
         __init__(Histogram2D self, IAxis axis_x, IAxis axis_y) -> Histogram2D
         __init__(Histogram2D self, IntensityData data) -> Histogram2D
-        __init__(Histogram2D self, vdouble2d_t data) -> Histogram2D
-        Histogram2D::Histogram2D(std::vector< std::vector< double >> data)
+        Histogram2D::Histogram2D(const Powerfield< double > &data)
 
-        Constructor for 2D histograms from numpy array (thanks to swig) 
+        Constructor for 2D histograms from basic  Powerfield object. 
 
         """
         _libBornAgainDevice.Histogram2D_swiginit(self, _libBornAgainDevice.new_Histogram2D(*args))
@@ -5146,7 +5139,9 @@ class Histogram2D(IHistogram):
         clone(Histogram2D self) -> Histogram2D
         Histogram2D * Histogram2D::clone() const override
 
-        Returns clone of other histogram. 
+        Constructor for 2D histograms from numpy array (thanks to swig)
+
+        Returns clone of other histogram 
 
         """
         return _libBornAgainDevice.Histogram2D_clone(self)
@@ -5436,14 +5431,6 @@ class SimulationResult(object):
         """
         return _libBornAgainDevice.SimulationResult_size(self)
 
-    def max(self):
-        r"""
-        max(SimulationResult self) -> double
-        double SimulationResult::max() const
-
-        """
-        return _libBornAgainDevice.SimulationResult_max(self)
-
     def empty(self):
         r"""
         empty(SimulationResult self) -> bool
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index d16f80f89d016830a201851d78d3d45f51c77c88..32df9be11570459cae4ce20493d8dccea26340fb 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -3116,95 +3116,97 @@ namespace Swig {
 #define SWIGTYPE_p_IDetector swig_types[16]
 #define SWIGTYPE_p_IDetector2D swig_types[17]
 #define SWIGTYPE_p_IDetectorResolution swig_types[18]
-#define SWIGTYPE_p_IFootprintFactor swig_types[19]
-#define SWIGTYPE_p_IHistogram swig_types[20]
-#define SWIGTYPE_p_INode swig_types[21]
-#define SWIGTYPE_p_IOFactory swig_types[22]
-#define SWIGTYPE_p_IPixel swig_types[23]
-#define SWIGTYPE_p_IResolutionFunction2D swig_types[24]
-#define SWIGTYPE_p_IShape2D swig_types[25]
-#define SWIGTYPE_p_Instrument swig_types[26]
-#define SWIGTYPE_p_Line swig_types[27]
-#define SWIGTYPE_p_OwningVectorT_IAxis_t swig_types[28]
-#define SWIGTYPE_p_PolMatrices swig_types[29]
-#define SWIGTYPE_p_Polygon swig_types[30]
-#define SWIGTYPE_p_PolygonPrivate swig_types[31]
-#define SWIGTYPE_p_PowerfieldIteratorT_double_PowerfieldT_double_t_t swig_types[32]
-#define SWIGTYPE_p_PowerfieldIteratorT_double_const_PowerfieldT_double_t_const_t swig_types[33]
-#define SWIGTYPE_p_PowerfieldT_CumulativeValue_t swig_types[34]
-#define SWIGTYPE_p_PowerfieldT_bool_t swig_types[35]
-#define SWIGTYPE_p_PowerfieldT_double_t swig_types[36]
-#define SWIGTYPE_p_RealLimits swig_types[37]
-#define SWIGTYPE_p_Rectangle swig_types[38]
-#define SWIGTYPE_p_RectangularDetector swig_types[39]
-#define SWIGTYPE_p_RectangularPixel swig_types[40]
-#define SWIGTYPE_p_ResolutionFunction2DGaussian swig_types[41]
-#define SWIGTYPE_p_SimulationResult swig_types[42]
-#define SWIGTYPE_p_SphericalDetector swig_types[43]
-#define SWIGTYPE_p_SpinMatrix swig_types[44]
-#define SWIGTYPE_p_Vec3T_double_t swig_types[45]
-#define SWIGTYPE_p_Vec3T_int_t swig_types[46]
-#define SWIGTYPE_p_Vec3T_std__complexT_double_t_t swig_types[47]
-#define SWIGTYPE_p_VerticalLine swig_types[48]
-#define SWIGTYPE_p_allocator_type swig_types[49]
-#define SWIGTYPE_p_bool swig_types[50]
-#define SWIGTYPE_p_char swig_types[51]
-#define SWIGTYPE_p_const_iterator swig_types[52]
-#define SWIGTYPE_p_corr_matrix_t swig_types[53]
-#define SWIGTYPE_p_difference_type swig_types[54]
-#define SWIGTYPE_p_double swig_types[55]
-#define SWIGTYPE_p_first_type swig_types[56]
-#define SWIGTYPE_p_int swig_types[57]
-#define SWIGTYPE_p_iterator swig_types[58]
-#define SWIGTYPE_p_key_type swig_types[59]
-#define SWIGTYPE_p_long_long swig_types[60]
-#define SWIGTYPE_p_mapped_type swig_types[61]
-#define SWIGTYPE_p_p_ICoordSystem swig_types[62]
-#define SWIGTYPE_p_p_PyObject swig_types[63]
-#define SWIGTYPE_p_parameters_t swig_types[64]
-#define SWIGTYPE_p_second_type swig_types[65]
-#define SWIGTYPE_p_short swig_types[66]
-#define SWIGTYPE_p_signed_char swig_types[67]
-#define SWIGTYPE_p_size_type swig_types[68]
-#define SWIGTYPE_p_std__allocatorT_Vec3T_double_t_t swig_types[69]
-#define SWIGTYPE_p_std__allocatorT_double_t swig_types[70]
-#define SWIGTYPE_p_std__allocatorT_int_t swig_types[71]
-#define SWIGTYPE_p_std__allocatorT_std__complexT_double_t_t swig_types[72]
-#define SWIGTYPE_p_std__allocatorT_std__pairT_double_double_t_t swig_types[73]
-#define SWIGTYPE_p_std__allocatorT_std__pairT_std__string_const_double_t_t swig_types[74]
-#define SWIGTYPE_p_std__allocatorT_std__string_t swig_types[75]
-#define SWIGTYPE_p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t swig_types[76]
-#define SWIGTYPE_p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t swig_types[77]
-#define SWIGTYPE_p_std__allocatorT_unsigned_long_t swig_types[78]
-#define SWIGTYPE_p_std__complexT_double_t swig_types[79]
-#define SWIGTYPE_p_std__functionT_void_fSimulationAreaIterator_const_RF_t swig_types[80]
-#define SWIGTYPE_p_std__invalid_argument swig_types[81]
-#define SWIGTYPE_p_std__lessT_std__string_t swig_types[82]
-#define SWIGTYPE_p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t swig_types[83]
-#define SWIGTYPE_p_std__pairT_double_double_t swig_types[84]
-#define SWIGTYPE_p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t swig_types[85]
-#define SWIGTYPE_p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t swig_types[86]
-#define SWIGTYPE_p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t swig_types[87]
-#define SWIGTYPE_p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t swig_types[88]
-#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[89]
-#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[90]
-#define SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t swig_types[91]
-#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[92]
-#define SWIGTYPE_p_std__vectorT_std__pairT_double_double_t_std__allocatorT_std__pairT_double_double_t_t_t swig_types[93]
-#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[94]
-#define SWIGTYPE_p_std__vectorT_std__unique_ptrT_DiffuseElement_t_std__allocatorT_std__unique_ptrT_DiffuseElement_t_t_t swig_types[95]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[96]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[97]
-#define SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t swig_types[98]
-#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[99]
-#define SWIGTYPE_p_swig__SwigPyIterator swig_types[100]
-#define SWIGTYPE_p_unsigned_char swig_types[101]
-#define SWIGTYPE_p_unsigned_int swig_types[102]
-#define SWIGTYPE_p_unsigned_long_long swig_types[103]
-#define SWIGTYPE_p_unsigned_short swig_types[104]
-#define SWIGTYPE_p_value_type swig_types[105]
-static swig_type_info *swig_types[107];
-static swig_module_info swig_module = {swig_types, 106, 0, 0, 0, 0};
+#define SWIGTYPE_p_IField swig_types[19]
+#define SWIGTYPE_p_IFootprintFactor swig_types[20]
+#define SWIGTYPE_p_IHistogram swig_types[21]
+#define SWIGTYPE_p_INode swig_types[22]
+#define SWIGTYPE_p_IOFactory swig_types[23]
+#define SWIGTYPE_p_IPixel swig_types[24]
+#define SWIGTYPE_p_IResolutionFunction2D swig_types[25]
+#define SWIGTYPE_p_IShape2D swig_types[26]
+#define SWIGTYPE_p_Instrument swig_types[27]
+#define SWIGTYPE_p_Line swig_types[28]
+#define SWIGTYPE_p_OwningVectorT_IAxis_t swig_types[29]
+#define SWIGTYPE_p_PolMatrices swig_types[30]
+#define SWIGTYPE_p_Polygon swig_types[31]
+#define SWIGTYPE_p_PolygonPrivate swig_types[32]
+#define SWIGTYPE_p_PowerfieldIteratorT_double_PowerfieldT_double_t_t swig_types[33]
+#define SWIGTYPE_p_PowerfieldIteratorT_double_const_PowerfieldT_double_t_const_t swig_types[34]
+#define SWIGTYPE_p_PowerfieldT_CumulativeValue_t swig_types[35]
+#define SWIGTYPE_p_PowerfieldT_bool_t swig_types[36]
+#define SWIGTYPE_p_PowerfieldT_double_t swig_types[37]
+#define SWIGTYPE_p_RealLimits swig_types[38]
+#define SWIGTYPE_p_Rectangle swig_types[39]
+#define SWIGTYPE_p_RectangularDetector swig_types[40]
+#define SWIGTYPE_p_RectangularPixel swig_types[41]
+#define SWIGTYPE_p_ResolutionFunction2DGaussian swig_types[42]
+#define SWIGTYPE_p_SimulationResult swig_types[43]
+#define SWIGTYPE_p_SphericalDetector swig_types[44]
+#define SWIGTYPE_p_SpinMatrix swig_types[45]
+#define SWIGTYPE_p_Vec3T_double_t swig_types[46]
+#define SWIGTYPE_p_Vec3T_int_t swig_types[47]
+#define SWIGTYPE_p_Vec3T_std__complexT_double_t_t swig_types[48]
+#define SWIGTYPE_p_VerticalLine swig_types[49]
+#define SWIGTYPE_p_allocator_type swig_types[50]
+#define SWIGTYPE_p_bool swig_types[51]
+#define SWIGTYPE_p_char swig_types[52]
+#define SWIGTYPE_p_const_iterator swig_types[53]
+#define SWIGTYPE_p_corr_matrix_t swig_types[54]
+#define SWIGTYPE_p_difference_type swig_types[55]
+#define SWIGTYPE_p_double swig_types[56]
+#define SWIGTYPE_p_first_type swig_types[57]
+#define SWIGTYPE_p_int swig_types[58]
+#define SWIGTYPE_p_iterator swig_types[59]
+#define SWIGTYPE_p_key_type swig_types[60]
+#define SWIGTYPE_p_long_long swig_types[61]
+#define SWIGTYPE_p_mapped_type swig_types[62]
+#define SWIGTYPE_p_p_ICoordSystem swig_types[63]
+#define SWIGTYPE_p_p_PyObject swig_types[64]
+#define SWIGTYPE_p_parameters_t swig_types[65]
+#define SWIGTYPE_p_second_type swig_types[66]
+#define SWIGTYPE_p_short swig_types[67]
+#define SWIGTYPE_p_signed_char swig_types[68]
+#define SWIGTYPE_p_size_type swig_types[69]
+#define SWIGTYPE_p_std__allocatorT_Vec3T_double_t_t swig_types[70]
+#define SWIGTYPE_p_std__allocatorT_double_t swig_types[71]
+#define SWIGTYPE_p_std__allocatorT_int_t swig_types[72]
+#define SWIGTYPE_p_std__allocatorT_std__complexT_double_t_t swig_types[73]
+#define SWIGTYPE_p_std__allocatorT_std__pairT_double_double_t_t swig_types[74]
+#define SWIGTYPE_p_std__allocatorT_std__pairT_std__string_const_double_t_t swig_types[75]
+#define SWIGTYPE_p_std__allocatorT_std__string_t swig_types[76]
+#define SWIGTYPE_p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t swig_types[77]
+#define SWIGTYPE_p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t swig_types[78]
+#define SWIGTYPE_p_std__allocatorT_unsigned_long_t swig_types[79]
+#define SWIGTYPE_p_std__complexT_double_t swig_types[80]
+#define SWIGTYPE_p_std__functionT_void_fSimulationAreaIterator_const_RF_t swig_types[81]
+#define SWIGTYPE_p_std__invalid_argument swig_types[82]
+#define SWIGTYPE_p_std__lessT_std__string_t swig_types[83]
+#define SWIGTYPE_p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t swig_types[84]
+#define SWIGTYPE_p_std__pairT_double_double_t swig_types[85]
+#define SWIGTYPE_p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t swig_types[86]
+#define SWIGTYPE_p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t swig_types[87]
+#define SWIGTYPE_p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t swig_types[88]
+#define SWIGTYPE_p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t swig_types[89]
+#define SWIGTYPE_p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t swig_types[90]
+#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[91]
+#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[92]
+#define SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t swig_types[93]
+#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[94]
+#define SWIGTYPE_p_std__vectorT_std__pairT_double_double_t_std__allocatorT_std__pairT_double_double_t_t_t swig_types[95]
+#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[96]
+#define SWIGTYPE_p_std__vectorT_std__unique_ptrT_DiffuseElement_t_std__allocatorT_std__unique_ptrT_DiffuseElement_t_t_t swig_types[97]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[98]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[99]
+#define SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t swig_types[100]
+#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[101]
+#define SWIGTYPE_p_swig__SwigPyIterator swig_types[102]
+#define SWIGTYPE_p_unsigned_char swig_types[103]
+#define SWIGTYPE_p_unsigned_int swig_types[104]
+#define SWIGTYPE_p_unsigned_long_long swig_types[105]
+#define SWIGTYPE_p_unsigned_short swig_types[106]
+#define SWIGTYPE_p_value_type swig_types[107]
+static swig_type_info *swig_types[109];
+static swig_module_info swig_module = {swig_types, 108, 0, 0, 0, 0};
 #define SWIG_TypeQuery(name) SWIG_TypeQueryModule(&swig_module, &swig_module, name)
 #define SWIG_MangledTypeQuery(name) SWIG_MangledTypeQueryModule(&swig_module, &swig_module, name)
 
@@ -26991,87 +26993,86 @@ SWIGINTERN PyObject *vector_R3_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_IntensityData__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+SWIGINTERN PyObject *_wrap_new_IField__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
   PyObject *resultobj = 0;
-  Powerfield< double > *result = 0 ;
+  IField *result = 0 ;
   
   if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Powerfield< double > *)new Powerfield< double >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NEW |  0 );
+  result = (IField *)new IField();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_IField, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_new_IntensityData__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_IField__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = 0 ;
+  std::vector< IAxis *,std::allocator< IAxis * > > *arg1 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  Powerfield< double > *result = 0 ;
+  IField *result = 0 ;
   
   if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_PowerfieldT_double_t,  0 );
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t,  0  | 0);
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_IntensityData" "', argument " "1"" of type '" "Powerfield< double > &&""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_IField" "', argument " "1"" of type '" "std::vector< IAxis *,std::allocator< IAxis * > > const &""'"); 
   }
   if (!argp1) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_IntensityData" "', argument " "1"" of type '" "Powerfield< double > &&""'"); 
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_IField" "', argument " "1"" of type '" "std::vector< IAxis *,std::allocator< IAxis * > > const &""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = (Powerfield< double > *)new Powerfield< double >((Powerfield< double > &&)*arg1);
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NEW |  0 );
+  arg1 = reinterpret_cast< std::vector< IAxis *,std::allocator< IAxis * > > * >(argp1);
+  result = (IField *)new IField((std::vector< IAxis *,std::allocator< IAxis * > > const &)*arg1);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_IField, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_new_IntensityData(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_new_IField(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[2] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "new_IntensityData", 0, 1, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "new_IField", 0, 1, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_IntensityData__SWIG_0(self, argc, argv);
+    return _wrap_new_IField__SWIG_0(self, argc, argv);
   }
   if (argc == 1) {
     int _v;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NO_NULL);
+    int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t, SWIG_POINTER_NO_NULL | 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_IntensityData__SWIG_1(self, argc, argv);
+      return _wrap_new_IField__SWIG_1(self, argc, argv);
     }
   }
   
 fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_IntensityData'.\n"
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_IField'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Powerfield< double >::Powerfield()\n"
-    "    Powerfield< double >::Powerfield(Powerfield< double > &&)\n");
+    "    IField::IField()\n"
+    "    IField::IField(std::vector< IAxis *,std::allocator< IAxis * > > const &)\n");
   return 0;
 }
 
 
-SWIGINTERN PyObject *_wrap_delete_IntensityData(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_delete_IField(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
   PyObject *swig_obj[1] ;
   
   if (!args) SWIG_fail;
   swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_DISOWN |  0 );
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, SWIG_POINTER_DISOWN |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "delete_IntensityData" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "delete_IField" "', argument " "1"" of type '" "IField *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
+  arg1 = reinterpret_cast< IField * >(argp1);
   delete arg1;
   resultobj = SWIG_Py_Void();
   return resultobj;
@@ -27080,305 +27081,471 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_clone(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_rank(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
   PyObject *swig_obj[1] ;
-  Powerfield< double > *result = 0 ;
+  size_t result;
   
   if (!args) SWIG_fail;
   swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_clone" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_rank" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = (Powerfield< double > *)((Powerfield< double > const *)arg1)->clone();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  arg1 = reinterpret_cast< IField * >(argp1);
+  result = ((IField const *)arg1)->rank();
+  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_copyFrom(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_axis(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  Powerfield< double > *arg2 = 0 ;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
+  size_t val2 ;
+  int ecode2 = 0 ;
   PyObject *swig_obj[2] ;
+  IAxis *result = 0 ;
   
-  if (!SWIG_Python_UnpackTuple(args, "IntensityData_copyFrom", 2, 2, swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!SWIG_Python_UnpackTuple(args, "IField_axis", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_copyFrom" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
-  }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_PowerfieldT_double_t,  0  | 0);
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_copyFrom" "', argument " "2"" of type '" "Powerfield< double > const &""'"); 
-  }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_copyFrom" "', argument " "2"" of type '" "Powerfield< double > const &""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_axis" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg2 = reinterpret_cast< Powerfield< double > * >(argp2);
-  (arg1)->copyFrom((Powerfield< double > const &)*arg2);
-  resultobj = SWIG_Py_Void();
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_axis" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  result = (IAxis *) &((IField const *)arg1)->axis(arg2);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_IAxis, 0 |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_meanValues(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_getAxisValue__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
+  size_t arg3 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  Powerfield< double > *result = 0 ;
+  size_t val2 ;
+  int ecode2 = 0 ;
+  size_t val3 ;
+  int ecode3 = 0 ;
+  double result;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_meanValues" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxisValue" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = (Powerfield< double > *)((Powerfield< double > const *)arg1)->meanValues();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxisValue" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
+  if (!SWIG_IsOK(ecode3)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IField_getAxisValue" "', argument " "3"" of type '" "size_t""'");
+  } 
+  arg3 = static_cast< size_t >(val3);
+  result = (double)((IField const *)arg1)->getAxisValue(arg2,arg3);
+  resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_addAxis__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IField_getAxisValue__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  IAxis *arg2 = 0 ;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
+  std::string *arg3 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
+  size_t val2 ;
+  int ecode2 = 0 ;
+  int res3 = SWIG_OLDOBJ ;
+  double result;
   
-  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_addAxis" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxisValue" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_IAxis,  0  | 0);
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "IAxis const &""'"); 
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxisValue" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  {
+    std::string *ptr = (std::string *)0;
+    res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
+    if (!SWIG_IsOK(res3)) {
+      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IField_getAxisValue" "', argument " "3"" of type '" "std::string const &""'"); 
+    }
+    if (!ptr) {
+      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IField_getAxisValue" "', argument " "3"" of type '" "std::string const &""'"); 
+    }
+    arg3 = ptr;
   }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "IAxis const &""'"); 
+  result = (double)((IField const *)arg1)->getAxisValue(arg2,(std::string const &)*arg3);
+  resultobj = SWIG_From_double(static_cast< double >(result));
+  if (SWIG_IsNewObj(res3)) delete arg3;
+  return resultobj;
+fail:
+  if (SWIG_IsNewObj(res3)) delete arg3;
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_IField_getAxisValue(PyObject *self, PyObject *args) {
+  Py_ssize_t argc;
+  PyObject *argv[4] = {
+    0
+  };
+  
+  if (!(argc = SWIG_Python_UnpackTuple(args, "IField_getAxisValue", 0, 3, argv))) SWIG_fail;
+  --argc;
+  if (argc == 3) {
+    int _v;
+    void *vptr = 0;
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IField, 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      {
+        int res = SWIG_AsVal_size_t(argv[1], NULL);
+        _v = SWIG_CheckState(res);
+      }
+      if (_v) {
+        {
+          int res = SWIG_AsVal_size_t(argv[2], NULL);
+          _v = SWIG_CheckState(res);
+        }
+        if (_v) {
+          return _wrap_IField_getAxisValue__SWIG_0(self, argc, argv);
+        }
+      }
+    }
   }
-  arg2 = reinterpret_cast< IAxis * >(argp2);
-  (arg1)->addAxis((IAxis const &)*arg2);
-  resultobj = SWIG_Py_Void();
+  if (argc == 3) {
+    int _v;
+    void *vptr = 0;
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IField, 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      {
+        int res = SWIG_AsVal_size_t(argv[1], NULL);
+        _v = SWIG_CheckState(res);
+      }
+      if (_v) {
+        int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
+        _v = SWIG_CheckState(res);
+        if (_v) {
+          return _wrap_IField_getAxisValue__SWIG_1(self, argc, argv);
+        }
+      }
+    }
+  }
+  
+fail:
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IField_getAxisValue'.\n"
+    "  Possible C/C++ prototypes are:\n"
+    "    IField::getAxisValue(size_t,size_t) const\n"
+    "    IField::getAxisValue(size_t,std::string const &) const\n");
+  return 0;
+}
+
+
+SWIGINTERN PyObject *_wrap_IField_getAxesValues(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  size_t val2 ;
+  int ecode2 = 0 ;
+  PyObject *swig_obj[2] ;
+  std::vector< double,std::allocator< double > > result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "IField_getAxesValues", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxesValues" "', argument " "1"" of type '" "IField const *""'"); 
+  }
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxesValues" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  result = ((IField const *)arg1)->getAxesValues(arg2);
+  resultobj = swig::from(static_cast< std::vector< double,std::allocator< double > > >(result));
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_addAxis__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IField_getAxisBin__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  std::string *arg2 = 0 ;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
   size_t arg3 ;
-  double arg4 ;
-  double arg5 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  int res2 = SWIG_OLDOBJ ;
+  size_t val2 ;
+  int ecode2 = 0 ;
   size_t val3 ;
   int ecode3 = 0 ;
-  double val4 ;
-  int ecode4 = 0 ;
-  double val5 ;
-  int ecode5 = 0 ;
+  Bin1D result;
   
-  if ((nobjs < 5) || (nobjs > 5)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_addAxis" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxisBin" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxisBin" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
+  if (!SWIG_IsOK(ecode3)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IField_getAxisBin" "', argument " "3"" of type '" "size_t""'");
+  } 
+  arg3 = static_cast< size_t >(val3);
+  result = ((IField const *)arg1)->getAxisBin(arg2,arg3);
+  resultobj = SWIG_NewPointerObj((new Bin1D(static_cast< const Bin1D& >(result))), SWIGTYPE_p_Bin1D, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_IField_getAxisBin__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
+  std::string *arg3 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  size_t val2 ;
+  int ecode2 = 0 ;
+  int res3 = SWIG_OLDOBJ ;
+  Bin1D result;
+  
+  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxisBin" "', argument " "1"" of type '" "IField const *""'"); 
+  }
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxisBin" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
   {
     std::string *ptr = (std::string *)0;
-    res2 = SWIG_AsPtr_std_string(swig_obj[1], &ptr);
-    if (!SWIG_IsOK(res2)) {
-      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "std::string const &""'"); 
+    res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
+    if (!SWIG_IsOK(res3)) {
+      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IField_getAxisBin" "', argument " "3"" of type '" "std::string const &""'"); 
     }
     if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "std::string const &""'"); 
+      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IField_getAxisBin" "', argument " "3"" of type '" "std::string const &""'"); 
     }
-    arg2 = ptr;
+    arg3 = ptr;
   }
-  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
-  if (!SWIG_IsOK(ecode3)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IntensityData_addAxis" "', argument " "3"" of type '" "size_t""'");
-  } 
-  arg3 = static_cast< size_t >(val3);
-  ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
-  if (!SWIG_IsOK(ecode4)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "IntensityData_addAxis" "', argument " "4"" of type '" "double""'");
-  } 
-  arg4 = static_cast< double >(val4);
-  ecode5 = SWIG_AsVal_double(swig_obj[4], &val5);
-  if (!SWIG_IsOK(ecode5)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "IntensityData_addAxis" "', argument " "5"" of type '" "double""'");
-  } 
-  arg5 = static_cast< double >(val5);
-  (arg1)->addAxis((std::string const &)*arg2,arg3,arg4,arg5);
-  resultobj = SWIG_Py_Void();
-  if (SWIG_IsNewObj(res2)) delete arg2;
+  result = ((IField const *)arg1)->getAxisBin(arg2,(std::string const &)*arg3);
+  resultobj = SWIG_NewPointerObj((new Bin1D(static_cast< const Bin1D& >(result))), SWIGTYPE_p_Bin1D, SWIG_POINTER_OWN |  0 );
+  if (SWIG_IsNewObj(res3)) delete arg3;
   return resultobj;
 fail:
-  if (SWIG_IsNewObj(res2)) delete arg2;
+  if (SWIG_IsNewObj(res3)) delete arg3;
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_addAxis(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_getAxisBin(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[6] = {
+  PyObject *argv[4] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_addAxis", 0, 5, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "IField_getAxisBin", 0, 3, argv))) SWIG_fail;
   --argc;
-  if (argc == 2) {
+  if (argc == 3) {
     int _v;
     void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IField, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      int res = SWIG_ConvertPtr(argv[1], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
-      _v = SWIG_CheckState(res);
+      {
+        int res = SWIG_AsVal_size_t(argv[1], NULL);
+        _v = SWIG_CheckState(res);
+      }
       if (_v) {
-        return _wrap_IntensityData_addAxis__SWIG_0(self, argc, argv);
+        {
+          int res = SWIG_AsVal_size_t(argv[2], NULL);
+          _v = SWIG_CheckState(res);
+        }
+        if (_v) {
+          return _wrap_IField_getAxisBin__SWIG_0(self, argc, argv);
+        }
       }
     }
   }
-  if (argc == 5) {
+  if (argc == 3) {
     int _v;
     void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IField, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      int res = SWIG_AsPtr_std_string(argv[1], (std::string**)(0));
-      _v = SWIG_CheckState(res);
+      {
+        int res = SWIG_AsVal_size_t(argv[1], NULL);
+        _v = SWIG_CheckState(res);
+      }
       if (_v) {
-        {
-          int res = SWIG_AsVal_size_t(argv[2], NULL);
-          _v = SWIG_CheckState(res);
-        }
+        int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
+        _v = SWIG_CheckState(res);
         if (_v) {
-          {
-            int res = SWIG_AsVal_double(argv[3], NULL);
-            _v = SWIG_CheckState(res);
-          }
-          if (_v) {
-            {
-              int res = SWIG_AsVal_double(argv[4], NULL);
-              _v = SWIG_CheckState(res);
-            }
-            if (_v) {
-              return _wrap_IntensityData_addAxis__SWIG_1(self, argc, argv);
-            }
-          }
+          return _wrap_IField_getAxisBin__SWIG_1(self, argc, argv);
         }
       }
     }
   }
   
 fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_addAxis'.\n"
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IField_getAxisBin'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Powerfield< double >::addAxis(IAxis const &)\n"
-    "    Powerfield< double >::addAxis(std::string const &,size_t,double,double)\n");
+    "    IField::getAxisBin(size_t,size_t) const\n"
+    "    IField::getAxisBin(size_t,std::string const &) const\n");
   return 0;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_axis(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_getAxisBinIndex__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
   size_t arg2 ;
+  std::string *arg3 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
   size_t val2 ;
   int ecode2 = 0 ;
-  PyObject *swig_obj[2] ;
-  IAxis *result = 0 ;
+  int res3 = SWIG_OLDOBJ ;
+  size_t result;
   
-  if (!SWIG_Python_UnpackTuple(args, "IntensityData_axis", 2, 2, swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_axis" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxisBinIndex" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
+  arg1 = reinterpret_cast< IField * >(argp1);
   ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
   if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_axis" "', argument " "2"" of type '" "size_t""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxisBinIndex" "', argument " "2"" of type '" "size_t""'");
   } 
   arg2 = static_cast< size_t >(val2);
-  result = (IAxis *) &((Powerfield< double > const *)arg1)->axis(arg2);
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_IAxis, 0 |  0 );
+  {
+    std::string *ptr = (std::string *)0;
+    res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
+    if (!SWIG_IsOK(res3)) {
+      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IField_getAxisBinIndex" "', argument " "3"" of type '" "std::string const &""'"); 
+    }
+    if (!ptr) {
+      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IField_getAxisBinIndex" "', argument " "3"" of type '" "std::string const &""'"); 
+    }
+    arg3 = ptr;
+  }
+  result = ((IField const *)arg1)->getAxisBinIndex(arg2,(std::string const &)*arg3);
+  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
+  if (SWIG_IsNewObj(res3)) delete arg3;
   return resultobj;
 fail:
+  if (SWIG_IsNewObj(res3)) delete arg3;
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_rank(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_getAxesBinIndices(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  size_t result;
+  size_t val2 ;
+  int ecode2 = 0 ;
+  PyObject *swig_obj[2] ;
+  std::vector< int,std::allocator< int > > result;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!SWIG_Python_UnpackTuple(args, "IField_getAxesBinIndices", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_rank" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxesBinIndices" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = ((Powerfield< double > const *)arg1)->rank();
-  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxesBinIndices" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  result = ((IField const *)arg1)->getAxesBinIndices(arg2);
+  resultobj = swig::from(static_cast< std::vector< int,std::allocator< int > > >(result));
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAllocatedSize(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_getAxisBinIndex__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
+  size_t arg2 ;
+  size_t arg3 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
+  size_t val2 ;
+  int ecode2 = 0 ;
+  size_t val3 ;
+  int ecode3 = 0 ;
   size_t result;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAllocatedSize" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_getAxisBinIndex" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = ((Powerfield< double > const *)arg1)->getAllocatedSize();
+  arg1 = reinterpret_cast< IField * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IField_getAxisBinIndex" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
+  if (!SWIG_IsOK(ecode3)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IField_getAxisBinIndex" "', argument " "3"" of type '" "size_t""'");
+  } 
+  arg3 = static_cast< size_t >(val3);
+  result = ((IField const *)arg1)->getAxisBinIndex(arg2,arg3);
   resultobj = SWIG_From_size_t(static_cast< size_t >(result));
   return resultobj;
 fail:
@@ -27386,449 +27553,561 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAllSizes(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_getAxisBinIndex(PyObject *self, PyObject *args) {
+  Py_ssize_t argc;
+  PyObject *argv[4] = {
+    0
+  };
+  
+  if (!(argc = SWIG_Python_UnpackTuple(args, "IField_getAxisBinIndex", 0, 3, argv))) SWIG_fail;
+  --argc;
+  if (argc == 3) {
+    int _v;
+    void *vptr = 0;
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IField, 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      {
+        int res = SWIG_AsVal_size_t(argv[1], NULL);
+        _v = SWIG_CheckState(res);
+      }
+      if (_v) {
+        {
+          int res = SWIG_AsVal_size_t(argv[2], NULL);
+          _v = SWIG_CheckState(res);
+        }
+        if (_v) {
+          return _wrap_IField_getAxisBinIndex__SWIG_1(self, argc, argv);
+        }
+      }
+    }
+  }
+  if (argc == 3) {
+    int _v;
+    void *vptr = 0;
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IField, 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      {
+        int res = SWIG_AsVal_size_t(argv[1], NULL);
+        _v = SWIG_CheckState(res);
+      }
+      if (_v) {
+        int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
+        _v = SWIG_CheckState(res);
+        if (_v) {
+          return _wrap_IField_getAxisBinIndex__SWIG_0(self, argc, argv);
+        }
+      }
+    }
+  }
+  
+fail:
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IField_getAxisBinIndex'.\n"
+    "  Possible C/C++ prototypes are:\n"
+    "    IField::getAxisBinIndex(size_t,std::string const &) const\n"
+    "    IField::getAxisBinIndex(size_t,size_t) const\n");
+  return 0;
+}
+
+
+SWIGINTERN PyObject *_wrap_IField_toGlobalIndex(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
+  std::vector< unsigned int,std::allocator< unsigned int > > *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  SwigValueWrapper< std::vector< size_t,std::allocator< size_t > > > result;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  size_t result;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!SWIG_Python_UnpackTuple(args, "IField_toGlobalIndex", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAllSizes" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_toGlobalIndex" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = ((Powerfield< double > const *)arg1)->getAllSizes();
-  resultobj = SWIG_NewPointerObj((new std::vector< size_t,std::allocator< size_t > >(static_cast< const std::vector< size_t,std::allocator< size_t > >& >(result))), SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t, SWIG_POINTER_OWN |  0 );
+  arg1 = reinterpret_cast< IField * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IField_toGlobalIndex" "', argument " "2"" of type '" "std::vector< unsigned int,std::allocator< unsigned int > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IField_toGlobalIndex" "', argument " "2"" of type '" "std::vector< unsigned int,std::allocator< unsigned int > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< std::vector< unsigned int,std::allocator< unsigned int > > * >(argp2);
+  result = ((IField const *)arg1)->toGlobalIndex((std::vector< unsigned int,std::allocator< unsigned int > > const &)*arg2);
+  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getRawDataVector(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IField_findGlobalIndex(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IField *arg1 = (IField *) 0 ;
+  std::vector< double,std::allocator< double > > *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  std::vector< double,std::allocator< double > > result;
+  int res2 = SWIG_OLDOBJ ;
+  PyObject *swig_obj[2] ;
+  size_t result;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!SWIG_Python_UnpackTuple(args, "IField_findGlobalIndex", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IField, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getRawDataVector" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IField_findGlobalIndex" "', argument " "1"" of type '" "IField const *""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = ((Powerfield< double > const *)arg1)->getRawDataVector();
-  resultobj = swig::from(static_cast< std::vector< double,std::allocator< double > > >(result));
+  arg1 = reinterpret_cast< IField * >(argp1);
+  {
+    std::vector< double,std::allocator< double > > *ptr = (std::vector< double,std::allocator< double > > *)0;
+    res2 = swig::asptr(swig_obj[1], &ptr);
+    if (!SWIG_IsOK(res2)) {
+      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IField_findGlobalIndex" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
+    }
+    if (!ptr) {
+      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IField_findGlobalIndex" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
+    }
+    arg2 = ptr;
+  }
+  result = ((IField const *)arg1)->findGlobalIndex((std::vector< double,std::allocator< double > > const &)*arg2);
+  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
+  if (SWIG_IsNewObj(res2)) delete arg2;
   return resultobj;
 fail:
+  if (SWIG_IsNewObj(res2)) delete arg2;
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_totalSum(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *IField_swigregister(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *obj;
+  if (!SWIG_Python_UnpackTuple(args, "swigregister", 1, 1, &obj)) return NULL;
+  SWIG_TypeNewClientData(SWIGTYPE_p_IField, SWIG_NewClientData(obj));
+  return SWIG_Py_Void();
+}
+
+SWIGINTERN PyObject *IField_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  return SWIG_Python_InitShadowInstance(args);
+}
+
+SWIGINTERN PyObject *_wrap_new_IntensityData__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
+  Powerfield< double > *result = 0 ;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_totalSum" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
-  }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = (double)((Powerfield< double > const *)arg1)->totalSum();
-  resultobj = SWIG_From_double(static_cast< double >(result));
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Powerfield< double > *)new Powerfield< double >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_max(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_new_IntensityData__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IAxis *arg1 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
+  Powerfield< double > *result = 0 ;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_IAxis,  0  | 0);
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_max" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_IntensityData" "', argument " "1"" of type '" "IAxis const &""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = (double)((Powerfield< double > const *)arg1)->max();
-  resultobj = SWIG_From_double(static_cast< double >(result));
+  if (!argp1) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_IntensityData" "', argument " "1"" of type '" "IAxis const &""'"); 
+  }
+  arg1 = reinterpret_cast< IAxis * >(argp1);
+  result = (Powerfield< double > *)new Powerfield< double >((IAxis const &)*arg1);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_begin__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_IntensityData__SWIG_2(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  IAxis *arg1 = 0 ;
+  IAxis *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SwigValueWrapper< PowerfieldIterator< double,Powerfield< double > > > result;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  Powerfield< double > *result = 0 ;
   
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_IAxis,  0  | 0);
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_begin" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_IntensityData" "', argument " "1"" of type '" "IAxis const &""'"); 
   }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = (arg1)->begin();
-  resultobj = SWIG_NewPointerObj((new Powerfield< double >::iterator(static_cast< const Powerfield< double >::iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_PowerfieldT_double_t_t, SWIG_POINTER_OWN |  0 );
+  if (!argp1) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_IntensityData" "', argument " "1"" of type '" "IAxis const &""'"); 
+  }
+  arg1 = reinterpret_cast< IAxis * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_IAxis,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "new_IntensityData" "', argument " "2"" of type '" "IAxis const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_IntensityData" "', argument " "2"" of type '" "IAxis const &""'"); 
+  }
+  arg2 = reinterpret_cast< IAxis * >(argp2);
+  result = (Powerfield< double > *)new Powerfield< double >((IAxis const &)*arg1,(IAxis const &)*arg2);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_begin__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_IntensityData__SWIG_3(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  Powerfield< double > *arg1 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SwigValueWrapper< PowerfieldIterator< double const,Powerfield< double > const > > result;
+  Powerfield< double > *result = 0 ;
   
   if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_PowerfieldT_double_t,  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_begin" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_IntensityData" "', argument " "1"" of type '" "Powerfield< double > &&""'"); 
+  }
+  if (!argp1) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_IntensityData" "', argument " "1"" of type '" "Powerfield< double > &&""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = ((Powerfield< double > const *)arg1)->begin();
-  resultobj = SWIG_NewPointerObj((new Powerfield< double >::const_iterator(static_cast< const Powerfield< double >::const_iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_const_PowerfieldT_double_t_const_t, SWIG_POINTER_OWN |  0 );
+  result = (Powerfield< double > *)new Powerfield< double >((Powerfield< double > &&)*arg1);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_begin(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_new_IntensityData(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[2] = {
+  PyObject *argv[3] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_begin", 0, 1, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "new_IntensityData", 0, 2, argv))) SWIG_fail;
   --argc;
+  if (argc == 0) {
+    return _wrap_new_IntensityData__SWIG_0(self, argc, argv);
+  }
   if (argc == 1) {
     int _v;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
+    int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_IntensityData_begin__SWIG_0(self, argc, argv);
+      return _wrap_new_IntensityData__SWIG_1(self, argc, argv);
     }
   }
   if (argc == 1) {
     int _v;
     void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_NO_NULL);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_IntensityData_begin__SWIG_1(self, argc, argv);
+      return _wrap_new_IntensityData__SWIG_3(self, argc, argv);
+    }
+  }
+  if (argc == 2) {
+    int _v;
+    int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      int res = SWIG_ConvertPtr(argv[1], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
+      _v = SWIG_CheckState(res);
+      if (_v) {
+        return _wrap_new_IntensityData__SWIG_2(self, argc, argv);
+      }
     }
   }
   
 fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_begin'.\n"
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_IntensityData'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Powerfield< double >::begin()\n"
-    "    Powerfield< double >::begin() const\n");
+    "    Powerfield< double >::Powerfield()\n"
+    "    Powerfield< double >::Powerfield(IAxis const &)\n"
+    "    Powerfield< double >::Powerfield(IAxis const &,IAxis const &)\n"
+    "    Powerfield< double >::Powerfield(Powerfield< double > &&)\n");
   return 0;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_end__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_delete_IntensityData(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SwigValueWrapper< PowerfieldIterator< double,Powerfield< double > > > result;
+  PyObject *swig_obj[1] ;
   
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, SWIG_POINTER_DISOWN |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_end" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "delete_IntensityData" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = (arg1)->end();
-  resultobj = SWIG_NewPointerObj((new Powerfield< double >::iterator(static_cast< const Powerfield< double >::iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_PowerfieldT_double_t_t, SWIG_POINTER_OWN |  0 );
+  delete arg1;
+  resultobj = SWIG_Py_Void();
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_end__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IntensityData_clone(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SwigValueWrapper< PowerfieldIterator< double const,Powerfield< double > const > > result;
+  PyObject *swig_obj[1] ;
+  Powerfield< double > *result = 0 ;
   
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_end" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_clone" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  result = ((Powerfield< double > const *)arg1)->end();
-  resultobj = SWIG_NewPointerObj((new Powerfield< double >::const_iterator(static_cast< const Powerfield< double >::const_iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_const_PowerfieldT_double_t_const_t, SWIG_POINTER_OWN |  0 );
+  result = (Powerfield< double > *)((Powerfield< double > const *)arg1)->clone();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_end(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[2] = {
-    0
-  };
+SWIGINTERN PyObject *_wrap_IntensityData_copyFrom(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  Powerfield< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_end", 0, 1, argv))) SWIG_fail;
-  --argc;
-  if (argc == 1) {
-    int _v;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      return _wrap_IntensityData_end__SWIG_0(self, argc, argv);
-    }
+  if (!SWIG_Python_UnpackTuple(args, "IntensityData_copyFrom", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_copyFrom" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
   }
-  if (argc == 1) {
-    int _v;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      return _wrap_IntensityData_end__SWIG_1(self, argc, argv);
-    }
+  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_PowerfieldT_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_copyFrom" "', argument " "2"" of type '" "Powerfield< double > const &""'"); 
   }
-  
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_copyFrom" "', argument " "2"" of type '" "Powerfield< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Powerfield< double > * >(argp2);
+  (arg1)->copyFrom((Powerfield< double > const &)*arg2);
+  resultobj = SWIG_Py_Void();
+  return resultobj;
 fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_end'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    Powerfield< double >::end()\n"
-    "    Powerfield< double >::end() const\n");
-  return 0;
+  return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxesBinIndices(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IntensityData_meanValues(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  PyObject *swig_obj[2] ;
-  std::vector< int,std::allocator< int > > result;
+  PyObject *swig_obj[1] ;
+  Powerfield< double > *result = 0 ;
   
-  if (!SWIG_Python_UnpackTuple(args, "IntensityData_getAxesBinIndices", 2, 2, swig_obj)) SWIG_fail;
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxesBinIndices" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_meanValues" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxesBinIndices" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  result = ((Powerfield< double > const *)arg1)->getAxesBinIndices(arg2);
-  resultobj = swig::from(static_cast< std::vector< int,std::allocator< int > > >(result));
+  result = (Powerfield< double > *)((Powerfield< double > const *)arg1)->meanValues();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisBinIndex__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IntensityData_addAxis__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  size_t arg3 ;
+  IAxis *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  size_t val3 ;
-  int ecode3 = 0 ;
-  size_t result;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
   
-  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxisBinIndex" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_addAxis" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxisBinIndex" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
-  if (!SWIG_IsOK(ecode3)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IntensityData_getAxisBinIndex" "', argument " "3"" of type '" "size_t""'");
-  } 
-  arg3 = static_cast< size_t >(val3);
-  result = ((Powerfield< double > const *)arg1)->getAxisBinIndex(arg2,arg3);
-  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_IAxis,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "IAxis const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "IAxis const &""'"); 
+  }
+  arg2 = reinterpret_cast< IAxis * >(argp2);
+  (arg1)->addAxis((IAxis const &)*arg2);
+  resultobj = SWIG_Py_Void();
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisBinIndex__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IntensityData_addAxis__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  std::string *arg3 = 0 ;
+  std::string *arg2 = 0 ;
+  size_t arg3 ;
+  double arg4 ;
+  double arg5 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  int res3 = SWIG_OLDOBJ ;
-  size_t result;
+  int res2 = SWIG_OLDOBJ ;
+  size_t val3 ;
+  int ecode3 = 0 ;
+  double val4 ;
+  int ecode4 = 0 ;
+  double val5 ;
+  int ecode5 = 0 ;
   
-  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  if ((nobjs < 5) || (nobjs > 5)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxisBinIndex" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_addAxis" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxisBinIndex" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
   {
     std::string *ptr = (std::string *)0;
-    res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
-    if (!SWIG_IsOK(res3)) {
-      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IntensityData_getAxisBinIndex" "', argument " "3"" of type '" "std::string const &""'"); 
+    res2 = SWIG_AsPtr_std_string(swig_obj[1], &ptr);
+    if (!SWIG_IsOK(res2)) {
+      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "std::string const &""'"); 
     }
     if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_getAxisBinIndex" "', argument " "3"" of type '" "std::string const &""'"); 
+      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_addAxis" "', argument " "2"" of type '" "std::string const &""'"); 
     }
-    arg3 = ptr;
+    arg2 = ptr;
   }
-  result = ((Powerfield< double > const *)arg1)->getAxisBinIndex(arg2,(std::string const &)*arg3);
-  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
-  if (SWIG_IsNewObj(res3)) delete arg3;
+  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
+  if (!SWIG_IsOK(ecode3)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IntensityData_addAxis" "', argument " "3"" of type '" "size_t""'");
+  } 
+  arg3 = static_cast< size_t >(val3);
+  ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
+  if (!SWIG_IsOK(ecode4)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "IntensityData_addAxis" "', argument " "4"" of type '" "double""'");
+  } 
+  arg4 = static_cast< double >(val4);
+  ecode5 = SWIG_AsVal_double(swig_obj[4], &val5);
+  if (!SWIG_IsOK(ecode5)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "IntensityData_addAxis" "', argument " "5"" of type '" "double""'");
+  } 
+  arg5 = static_cast< double >(val5);
+  (arg1)->addAxis((std::string const &)*arg2,arg3,arg4,arg5);
+  resultobj = SWIG_Py_Void();
+  if (SWIG_IsNewObj(res2)) delete arg2;
   return resultobj;
 fail:
-  if (SWIG_IsNewObj(res3)) delete arg3;
+  if (SWIG_IsNewObj(res2)) delete arg2;
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisBinIndex(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_IntensityData_addAxis(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[4] = {
+  PyObject *argv[6] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_getAxisBinIndex", 0, 3, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_addAxis", 0, 5, argv))) SWIG_fail;
   --argc;
-  if (argc == 3) {
+  if (argc == 2) {
     int _v;
     void *vptr = 0;
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      {
-        int res = SWIG_AsVal_size_t(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
+      int res = SWIG_ConvertPtr(argv[1], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
+      _v = SWIG_CheckState(res);
       if (_v) {
-        {
-          int res = SWIG_AsVal_size_t(argv[2], NULL);
-          _v = SWIG_CheckState(res);
-        }
-        if (_v) {
-          return _wrap_IntensityData_getAxisBinIndex__SWIG_0(self, argc, argv);
-        }
+        return _wrap_IntensityData_addAxis__SWIG_0(self, argc, argv);
       }
     }
   }
-  if (argc == 3) {
+  if (argc == 5) {
     int _v;
     void *vptr = 0;
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      {
-        int res = SWIG_AsVal_size_t(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
+      int res = SWIG_AsPtr_std_string(argv[1], (std::string**)(0));
+      _v = SWIG_CheckState(res);
       if (_v) {
-        int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
-        _v = SWIG_CheckState(res);
+        {
+          int res = SWIG_AsVal_size_t(argv[2], NULL);
+          _v = SWIG_CheckState(res);
+        }
         if (_v) {
-          return _wrap_IntensityData_getAxisBinIndex__SWIG_1(self, argc, argv);
+          {
+            int res = SWIG_AsVal_double(argv[3], NULL);
+            _v = SWIG_CheckState(res);
+          }
+          if (_v) {
+            {
+              int res = SWIG_AsVal_double(argv[4], NULL);
+              _v = SWIG_CheckState(res);
+            }
+            if (_v) {
+              return _wrap_IntensityData_addAxis__SWIG_1(self, argc, argv);
+            }
+          }
         }
       }
     }
   }
   
 fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_getAxisBinIndex'.\n"
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_addAxis'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Powerfield< double >::getAxisBinIndex(size_t,size_t) const\n"
-    "    Powerfield< double >::getAxisBinIndex(size_t,std::string const &) const\n");
+    "    Powerfield< double >::addAxis(IAxis const &)\n"
+    "    Powerfield< double >::addAxis(std::string const &,size_t,double,double)\n");
   return 0;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_toGlobalIndex(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IntensityData_getAllocatedSize(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  std::vector< unsigned int,std::allocator< unsigned int > > *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
-  PyObject *swig_obj[2] ;
+  PyObject *swig_obj[1] ;
   size_t result;
   
-  if (!SWIG_Python_UnpackTuple(args, "IntensityData_toGlobalIndex", 2, 2, swig_obj)) SWIG_fail;
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_toGlobalIndex" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAllocatedSize" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t,  0  | 0);
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_toGlobalIndex" "', argument " "2"" of type '" "std::vector< unsigned int,std::allocator< unsigned int > > const &""'"); 
-  }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_toGlobalIndex" "', argument " "2"" of type '" "std::vector< unsigned int,std::allocator< unsigned int > > const &""'"); 
-  }
-  arg2 = reinterpret_cast< std::vector< unsigned int,std::allocator< unsigned int > > * >(argp2);
-  result = ((Powerfield< double > const *)arg1)->toGlobalIndex((std::vector< unsigned int,std::allocator< unsigned int > > const &)*arg2);
+  result = ((Powerfield< double > const *)arg1)->getAllocatedSize();
   resultobj = SWIG_From_size_t(static_cast< size_t >(result));
   return resultobj;
 fail:
@@ -27836,73 +28115,68 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_findGlobalIndex(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IntensityData_getAllSizes(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  std::vector< double,std::allocator< double > > *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  int res2 = SWIG_OLDOBJ ;
-  PyObject *swig_obj[2] ;
-  size_t result;
+  PyObject *swig_obj[1] ;
+  SwigValueWrapper< std::vector< size_t,std::allocator< size_t > > > result;
   
-  if (!SWIG_Python_UnpackTuple(args, "IntensityData_findGlobalIndex", 2, 2, swig_obj)) SWIG_fail;
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_findGlobalIndex" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAllSizes" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  {
-    std::vector< double,std::allocator< double > > *ptr = (std::vector< double,std::allocator< double > > *)0;
-    res2 = swig::asptr(swig_obj[1], &ptr);
-    if (!SWIG_IsOK(res2)) {
-      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IntensityData_findGlobalIndex" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_findGlobalIndex" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    arg2 = ptr;
+  result = ((Powerfield< double > const *)arg1)->getAllSizes();
+  resultobj = SWIG_NewPointerObj((new std::vector< size_t,std::allocator< size_t > >(static_cast< const std::vector< size_t,std::allocator< size_t > >& >(result))), SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_IntensityData_getRawDataVector(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  std::vector< double,std::allocator< double > > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getRawDataVector" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
-  result = ((Powerfield< double > const *)arg1)->findGlobalIndex((std::vector< double,std::allocator< double > > const &)*arg2);
-  resultobj = SWIG_From_size_t(static_cast< size_t >(result));
-  if (SWIG_IsNewObj(res2)) delete arg2;
+  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
+  result = ((Powerfield< double > const *)arg1)->getRawDataVector();
+  resultobj = swig::from(static_cast< std::vector< double,std::allocator< double > > >(result));
   return resultobj;
 fail:
-  if (SWIG_IsNewObj(res2)) delete arg2;
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisValue__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IntensityData_totalSum(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  size_t arg3 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  size_t val3 ;
-  int ecode3 = 0 ;
+  PyObject *swig_obj[1] ;
   double result;
   
-  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxisValue" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_totalSum" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxisValue" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
-  if (!SWIG_IsOK(ecode3)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IntensityData_getAxisValue" "', argument " "3"" of type '" "size_t""'");
-  } 
-  arg3 = static_cast< size_t >(val3);
-  result = (double)((Powerfield< double > const *)arg1)->getAxisValue(arg2,arg3);
+  result = (double)((Powerfield< double > const *)arg1)->totalSum();
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -27910,273 +28184,158 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisValue__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IntensityData_begin__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  std::string *arg3 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  int res3 = SWIG_OLDOBJ ;
-  double result;
+  SwigValueWrapper< PowerfieldIterator< double,Powerfield< double > > > result;
   
-  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxisValue" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_begin" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxisValue" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  {
-    std::string *ptr = (std::string *)0;
-    res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
-    if (!SWIG_IsOK(res3)) {
-      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IntensityData_getAxisValue" "', argument " "3"" of type '" "std::string const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_getAxisValue" "', argument " "3"" of type '" "std::string const &""'"); 
-    }
-    arg3 = ptr;
+  result = (arg1)->begin();
+  resultobj = SWIG_NewPointerObj((new Powerfield< double >::iterator(static_cast< const Powerfield< double >::iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_PowerfieldT_double_t_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_IntensityData_begin__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  SwigValueWrapper< PowerfieldIterator< double const,Powerfield< double > const > > result;
+  
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_begin" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
-  result = (double)((Powerfield< double > const *)arg1)->getAxisValue(arg2,(std::string const &)*arg3);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  if (SWIG_IsNewObj(res3)) delete arg3;
+  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
+  result = ((Powerfield< double > const *)arg1)->begin();
+  resultobj = SWIG_NewPointerObj((new Powerfield< double >::const_iterator(static_cast< const Powerfield< double >::const_iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_const_PowerfieldT_double_t_const_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
-  if (SWIG_IsNewObj(res3)) delete arg3;
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisValue(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_IntensityData_begin(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[4] = {
+  PyObject *argv[2] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_getAxisValue", 0, 3, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_begin", 0, 1, argv))) SWIG_fail;
   --argc;
-  if (argc == 3) {
+  if (argc == 1) {
     int _v;
     void *vptr = 0;
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      {
-        int res = SWIG_AsVal_size_t(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        {
-          int res = SWIG_AsVal_size_t(argv[2], NULL);
-          _v = SWIG_CheckState(res);
-        }
-        if (_v) {
-          return _wrap_IntensityData_getAxisValue__SWIG_0(self, argc, argv);
-        }
-      }
+      return _wrap_IntensityData_begin__SWIG_0(self, argc, argv);
     }
   }
-  if (argc == 3) {
+  if (argc == 1) {
     int _v;
     void *vptr = 0;
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      {
-        int res = SWIG_AsVal_size_t(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
-        _v = SWIG_CheckState(res);
-        if (_v) {
-          return _wrap_IntensityData_getAxisValue__SWIG_1(self, argc, argv);
-        }
-      }
+      return _wrap_IntensityData_begin__SWIG_1(self, argc, argv);
     }
   }
   
 fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_getAxisValue'.\n"
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_begin'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Powerfield< double >::getAxisValue(size_t,size_t) const\n"
-    "    Powerfield< double >::getAxisValue(size_t,std::string const &) const\n");
+    "    Powerfield< double >::begin()\n"
+    "    Powerfield< double >::begin() const\n");
   return 0;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxesValues(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  PyObject *swig_obj[2] ;
-  std::vector< double,std::allocator< double > > result;
-  
-  if (!SWIG_Python_UnpackTuple(args, "IntensityData_getAxesValues", 2, 2, swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxesValues" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
-  }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxesValues" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  result = ((Powerfield< double > const *)arg1)->getAxesValues(arg2);
-  resultobj = swig::from(static_cast< std::vector< double,std::allocator< double > > >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisBin__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IntensityData_end__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  size_t arg3 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  size_t val3 ;
-  int ecode3 = 0 ;
-  Bin1D result;
+  SwigValueWrapper< PowerfieldIterator< double,Powerfield< double > > > result;
   
-  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxisBin" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_end" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxisBin" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
-  if (!SWIG_IsOK(ecode3)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IntensityData_getAxisBin" "', argument " "3"" of type '" "size_t""'");
-  } 
-  arg3 = static_cast< size_t >(val3);
-  result = ((Powerfield< double > const *)arg1)->getAxisBin(arg2,arg3);
-  resultobj = SWIG_NewPointerObj((new Bin1D(static_cast< const Bin1D& >(result))), SWIGTYPE_p_Bin1D, SWIG_POINTER_OWN |  0 );
+  result = (arg1)->end();
+  resultobj = SWIG_NewPointerObj((new Powerfield< double >::iterator(static_cast< const Powerfield< double >::iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_PowerfieldT_double_t_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisBin__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_IntensityData_end__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  std::string *arg3 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  int res3 = SWIG_OLDOBJ ;
-  Bin1D result;
+  SwigValueWrapper< PowerfieldIterator< double const,Powerfield< double > const > > result;
   
-  if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_getAxisBin" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_end" "', argument " "1"" of type '" "Powerfield< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_getAxisBin" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  {
-    std::string *ptr = (std::string *)0;
-    res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
-    if (!SWIG_IsOK(res3)) {
-      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IntensityData_getAxisBin" "', argument " "3"" of type '" "std::string const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IntensityData_getAxisBin" "', argument " "3"" of type '" "std::string const &""'"); 
-    }
-    arg3 = ptr;
-  }
-  result = ((Powerfield< double > const *)arg1)->getAxisBin(arg2,(std::string const &)*arg3);
-  resultobj = SWIG_NewPointerObj((new Bin1D(static_cast< const Bin1D& >(result))), SWIGTYPE_p_Bin1D, SWIG_POINTER_OWN |  0 );
-  if (SWIG_IsNewObj(res3)) delete arg3;
+  result = ((Powerfield< double > const *)arg1)->end();
+  resultobj = SWIG_NewPointerObj((new Powerfield< double >::const_iterator(static_cast< const Powerfield< double >::const_iterator& >(result))), SWIGTYPE_p_PowerfieldIteratorT_double_const_PowerfieldT_double_t_const_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
-  if (SWIG_IsNewObj(res3)) delete arg3;
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_getAxisBin(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_IntensityData_end(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[4] = {
+  PyObject *argv[2] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_getAxisBin", 0, 3, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "IntensityData_end", 0, 1, argv))) SWIG_fail;
   --argc;
-  if (argc == 3) {
+  if (argc == 1) {
     int _v;
     void *vptr = 0;
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      {
-        int res = SWIG_AsVal_size_t(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        {
-          int res = SWIG_AsVal_size_t(argv[2], NULL);
-          _v = SWIG_CheckState(res);
-        }
-        if (_v) {
-          return _wrap_IntensityData_getAxisBin__SWIG_0(self, argc, argv);
-        }
-      }
+      return _wrap_IntensityData_end__SWIG_0(self, argc, argv);
     }
   }
-  if (argc == 3) {
+  if (argc == 1) {
     int _v;
     void *vptr = 0;
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_PowerfieldT_double_t, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      {
-        int res = SWIG_AsVal_size_t(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
-        _v = SWIG_CheckState(res);
-        if (_v) {
-          return _wrap_IntensityData_getAxisBin__SWIG_1(self, argc, argv);
-        }
-      }
+      return _wrap_IntensityData_end__SWIG_1(self, argc, argv);
     }
   }
   
 fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_getAxisBin'.\n"
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IntensityData_end'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Powerfield< double >::getAxisBin(size_t,size_t) const\n"
-    "    Powerfield< double >::getAxisBin(size_t,std::string const &) const\n");
+    "    Powerfield< double >::end()\n"
+    "    Powerfield< double >::end() const\n");
   return 0;
 }
 
@@ -28234,43 +28393,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IntensityData_setAxisSizes(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
-  size_t arg2 ;
-  int *arg3 = (int *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  void *argp3 = 0 ;
-  int res3 = 0 ;
-  PyObject *swig_obj[3] ;
-  
-  if (!SWIG_Python_UnpackTuple(args, "IntensityData_setAxisSizes", 3, 3, swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_PowerfieldT_double_t, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_setAxisSizes" "', argument " "1"" of type '" "Powerfield< double > *""'"); 
-  }
-  arg1 = reinterpret_cast< Powerfield< double > * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IntensityData_setAxisSizes" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  res3 = SWIG_ConvertPtr(swig_obj[2], &argp3,SWIGTYPE_p_int, 0 |  0 );
-  if (!SWIG_IsOK(res3)) {
-    SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IntensityData_setAxisSizes" "', argument " "3"" of type '" "int *""'"); 
-  }
-  arg3 = reinterpret_cast< int * >(argp3);
-  (arg1)->setAxisSizes(arg2,arg3);
-  resultobj = SWIG_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_IntensityData_setRawDataVector(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Powerfield< double > *arg1 = (Powerfield< double > *) 0 ;
@@ -33065,6 +33187,41 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_DetectorMask__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  IAxis *arg1 = 0 ;
+  IAxis *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  DetectorMask *result = 0 ;
+  
+  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_IAxis,  0  | 0);
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_DetectorMask" "', argument " "1"" of type '" "IAxis const &""'"); 
+  }
+  if (!argp1) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_DetectorMask" "', argument " "1"" of type '" "IAxis const &""'"); 
+  }
+  arg1 = reinterpret_cast< IAxis * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_IAxis,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "new_DetectorMask" "', argument " "2"" of type '" "IAxis const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_DetectorMask" "', argument " "2"" of type '" "IAxis const &""'"); 
+  }
+  arg2 = reinterpret_cast< IAxis * >(argp2);
+  result = (DetectorMask *)new DetectorMask((IAxis const &)*arg1,(IAxis const &)*arg2);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_DetectorMask, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_delete_DetectorMask(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   DetectorMask *arg1 = (DetectorMask *) 0 ;
@@ -33087,7 +33244,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_DetectorMask__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_DetectorMask__SWIG_2(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   DetectorMask *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -33113,11 +33270,11 @@ fail:
 
 SWIGINTERN PyObject *_wrap_new_DetectorMask(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[2] = {
+  PyObject *argv[3] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "new_DetectorMask", 0, 1, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "new_DetectorMask", 0, 2, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
     return _wrap_new_DetectorMask__SWIG_0(self, argc, argv);
@@ -33127,7 +33284,19 @@ SWIGINTERN PyObject *_wrap_new_DetectorMask(PyObject *self, PyObject *args) {
     int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_DetectorMask, SWIG_POINTER_NO_NULL | 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_DetectorMask__SWIG_1(self, argc, argv);
+      return _wrap_new_DetectorMask__SWIG_2(self, argc, argv);
+    }
+  }
+  if (argc == 2) {
+    int _v;
+    int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      int res = SWIG_ConvertPtr(argv[1], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
+      _v = SWIG_CheckState(res);
+      if (_v) {
+        return _wrap_new_DetectorMask__SWIG_1(self, argc, argv);
+      }
     }
   }
   
@@ -33135,6 +33304,7 @@ fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_DetectorMask'.\n"
     "  Possible C/C++ prototypes are:\n"
     "    DetectorMask::DetectorMask()\n"
+    "    DetectorMask::DetectorMask(IAxis const &,IAxis const &)\n"
     "    DetectorMask::DetectorMask(DetectorMask const &)\n");
   return 0;
 }
@@ -38683,96 +38853,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IHistogram_createFrom__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  std::string *arg1 = 0 ;
-  int res1 = SWIG_OLDOBJ ;
-  IHistogram *result = 0 ;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  {
-    std::string *ptr = (std::string *)0;
-    res1 = SWIG_AsPtr_std_string(swig_obj[0], &ptr);
-    if (!SWIG_IsOK(res1)) {
-      SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IHistogram_createFrom" "', argument " "1"" of type '" "std::string const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IHistogram_createFrom" "', argument " "1"" of type '" "std::string const &""'"); 
-    }
-    arg1 = ptr;
-  }
-  result = (IHistogram *)IHistogram::createFrom((std::string const &)*arg1);
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_IHistogram, SWIG_POINTER_OWN |  0 );
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return resultobj;
-fail:
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_IHistogram_createFrom__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *arg1 = 0 ;
-  int res1 = SWIG_OLDOBJ ;
-  IHistogram *result = 0 ;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  {
-    std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *ptr = (std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *)0;
-    res1 = swig::asptr(swig_obj[0], &ptr);
-    if (!SWIG_IsOK(res1)) {
-      SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IHistogram_createFrom" "', argument " "1"" of type '" "std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IHistogram_createFrom" "', argument " "1"" of type '" "std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &""'"); 
-    }
-    arg1 = ptr;
-  }
-  result = (IHistogram *)IHistogram::createFrom((std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &)*arg1);
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_IHistogram, SWIG_POINTER_OWN |  0 );
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return resultobj;
-fail:
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_IHistogram_createFrom(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[2] = {
-    0
-  };
-  
-  if (!(argc = SWIG_Python_UnpackTuple(args, "IHistogram_createFrom", 0, 1, argv))) SWIG_fail;
-  --argc;
-  if (argc == 1) {
-    int _v;
-    int res = SWIG_AsPtr_std_string(argv[0], (std::string**)(0));
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      return _wrap_IHistogram_createFrom__SWIG_0(self, argc, argv);
-    }
-  }
-  if (argc == 1) {
-    int _v;
-    int res = swig::asptr(argv[0], (std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > >**)(0));
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      return _wrap_IHistogram_createFrom__SWIG_1(self, argc, argv);
-    }
-  }
-  
-fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'IHistogram_createFrom'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    IHistogram::createFrom(std::string const &)\n"
-    "    IHistogram::createFrom(std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &)\n");
-  return 0;
-}
-
-
 SWIGINTERN PyObject *_wrap_IHistogram_createPowerfield__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   IHistogram *arg1 = (IHistogram *) 0 ;
@@ -39823,29 +39903,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Histogram2D__SWIG_4(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > arg1 ;
-  Histogram2D *result = 0 ;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  {
-    std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *ptr = (std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *)0;
-    int res = swig::asptr(swig_obj[0], &ptr);
-    if (!SWIG_IsOK(res) || !ptr) {
-      SWIG_exception_fail(SWIG_ArgError((ptr ? res : SWIG_TypeError)), "in method '" "new_Histogram2D" "', argument " "1"" of type '" "std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > >""'"); 
-    }
-    arg1 = *ptr;
-    if (SWIG_IsNewObj(res)) delete ptr;
-  }
-  result = (Histogram2D *)new Histogram2D(arg1);
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Histogram2D, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_new_Histogram2D(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[7] = {
@@ -39862,14 +39919,6 @@ SWIGINTERN PyObject *_wrap_new_Histogram2D(PyObject *self, PyObject *args) {
       return _wrap_new_Histogram2D__SWIG_3(self, argc, argv);
     }
   }
-  if (argc == 1) {
-    int _v;
-    int res = swig::asptr(argv[0], (std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > >**)(0));
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      return _wrap_new_Histogram2D__SWIG_4(self, argc, argv);
-    }
-  }
   if (argc == 2) {
     int _v;
     int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_IAxis, SWIG_POINTER_NO_NULL | 0);
@@ -39953,8 +40002,7 @@ fail:
     "    Histogram2D::Histogram2D(int,double,double,int,double,double)\n"
     "    Histogram2D::Histogram2D(int,std::vector< double,std::allocator< double > > const &,int,std::vector< double,std::allocator< double > > const &)\n"
     "    Histogram2D::Histogram2D(IAxis const &,IAxis const &)\n"
-    "    Histogram2D::Histogram2D(Powerfield< double > const &)\n"
-    "    Histogram2D::Histogram2D(std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > >)\n");
+    "    Histogram2D::Histogram2D(Powerfield< double > const &)\n");
   return 0;
 }
 
@@ -41599,29 +41647,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_SimulationResult_max(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  SimulationResult *arg1 = (SimulationResult *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
-  
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_SimulationResult, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SimulationResult_max" "', argument " "1"" of type '" "SimulationResult const *""'"); 
-  }
-  arg1 = reinterpret_cast< SimulationResult * >(argp1);
-  result = (double)((SimulationResult const *)arg1)->max();
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_SimulationResult_empty(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   SimulationResult *arg1 = (SimulationResult *) 0 ;
@@ -42723,106 +42748,87 @@ static PyMethodDef SwigMethods[] = {
 	 { "delete_vector_R3", _wrap_delete_vector_R3, METH_O, "delete_vector_R3(vector_R3 self)"},
 	 { "vector_R3_swigregister", vector_R3_swigregister, METH_O, NULL},
 	 { "vector_R3_swiginit", vector_R3_swiginit, METH_VARARGS, NULL},
-	 { "new_IntensityData", _wrap_new_IntensityData, METH_VARARGS, "\n"
-		"IntensityData()\n"
-		"new_IntensityData(IntensityData arg1) -> IntensityData\n"
-		"Powerfield< T >::Powerfield(Powerfield &&)\n"
-		"\n"
-		""},
-	 { "delete_IntensityData", _wrap_delete_IntensityData, METH_O, "\n"
-		"delete_IntensityData(IntensityData self)\n"
-		"Powerfield< T >::~Powerfield()\n"
+	 { "new_IField", _wrap_new_IField, METH_VARARGS, "\n"
+		"IField()\n"
+		"new_IField(std::vector< IAxis *,std::allocator< IAxis * > > const & axes) -> IField\n"
+		"IField::IField(const std::vector< IAxis * > &axes)\n"
 		"\n"
 		""},
-	 { "IntensityData_clone", _wrap_IntensityData_clone, METH_O, "\n"
-		"IntensityData_clone(IntensityData self) -> IntensityData\n"
-		"Powerfield< T > * Powerfield< T >::clone() const\n"
-		"\n"
-		""},
-	 { "IntensityData_copyFrom", _wrap_IntensityData_copyFrom, METH_VARARGS, "\n"
-		"IntensityData_copyFrom(IntensityData self, IntensityData other)\n"
-		"void Powerfield< T >::copyFrom(const Powerfield< T > &other)\n"
+	 { "delete_IField", _wrap_delete_IField, METH_O, "\n"
+		"delete_IField(IField self)\n"
+		"virtual IField::~IField()\n"
 		"\n"
 		""},
-	 { "IntensityData_meanValues", _wrap_IntensityData_meanValues, METH_O, "\n"
-		"IntensityData_meanValues(IntensityData self) -> IntensityData\n"
-		"Powerfield< double > * Powerfield< T >::meanValues() const\n"
+	 { "IField_rank", _wrap_IField_rank, METH_O, "\n"
+		"IField_rank(IField self) -> size_t\n"
+		"size_t IField::rank() const\n"
 		"\n"
-		""},
-	 { "IntensityData_addAxis", _wrap_IntensityData_addAxis, METH_VARARGS, "\n"
-		"IntensityData_addAxis(IntensityData self, IAxis new_axis)\n"
-		"IntensityData_addAxis(IntensityData self, std::string const & name, size_t size, double start, double end)\n"
-		"void Powerfield< T >::addAxis(const std::string &name, size_t size, double start, double end)\n"
+		"Returns number of dimensions. \n"
 		"\n"
 		""},
-	 { "IntensityData_axis", _wrap_IntensityData_axis, METH_VARARGS, "\n"
-		"IntensityData_axis(IntensityData self, size_t serial_number) -> IAxis\n"
-		"const IAxis & Powerfield< T >::axis(size_t serial_number) const\n"
+	 { "IField_axis", _wrap_IField_axis, METH_VARARGS, "\n"
+		"IField_axis(IField self, size_t serial_number) -> IAxis\n"
+		"const IAxis & IField::axis(size_t serial_number) const\n"
 		"\n"
 		"Returns axis with given serial number. \n"
 		"\n"
 		""},
-	 { "IntensityData_rank", _wrap_IntensityData_rank, METH_O, "\n"
-		"IntensityData_rank(IntensityData self) -> size_t\n"
-		"size_t Powerfield< T >::rank() const\n"
+	 { "IField_getAxisValue", _wrap_IField_getAxisValue, METH_VARARGS, "\n"
+		"IField_getAxisValue(IField self, size_t global_index, size_t i_selected_axis) -> double\n"
+		"IField_getAxisValue(IField self, size_t global_index, std::string const & axis_name) -> double\n"
+		"double IField::getAxisValue(size_t global_index, const std::string &axis_name) const\n"
 		"\n"
-		"Returns number of dimensions. \n"
+		"Returns the value of selected axis for given global_index.\n"
 		"\n"
-		""},
-	 { "IntensityData_getAllocatedSize", _wrap_IntensityData_getAllocatedSize, METH_O, "\n"
-		"IntensityData_getAllocatedSize(IntensityData self) -> size_t\n"
-		"size_t Powerfield< T >::getAllocatedSize() const\n"
+		"Parameters:\n"
+		"-----------\n"
 		"\n"
-		"Returns total size of data buffer (product of bin number in every dimension). \n"
+		"global_index: \n"
+		"The global index of this data structure.\n"
 		"\n"
-		""},
-	 { "IntensityData_getAllSizes", _wrap_IntensityData_getAllSizes, METH_O, "\n"
-		"IntensityData_getAllSizes(IntensityData self) -> std::vector< size_t,std::allocator< size_t > >\n"
-		"std::vector< size_t > Powerfield< T >::getAllSizes() const\n"
+		"axis_name: \n"
+		"The name of selected axis.\n"
 		"\n"
-		"Returns all sizes of its axes. \n"
+		"corresponding bin center of selected axis \n"
 		"\n"
 		""},
-	 { "IntensityData_getRawDataVector", _wrap_IntensityData_getRawDataVector, METH_O, "\n"
-		"IntensityData_getRawDataVector(IntensityData self) -> vdouble1d_t\n"
-		"std::vector< T > Powerfield< T >::getRawDataVector() const\n"
+	 { "IField_getAxesValues", _wrap_IField_getAxesValues, METH_VARARGS, "\n"
+		"IField_getAxesValues(IField self, size_t global_index) -> vdouble1d_t\n"
+		"std::vector< double > IField::getAxesValues(size_t global_index) const\n"
 		"\n"
-		"Returns copy of raw data vector. \n"
+		"Returns values on all defined axes for given globalbin number\n"
 		"\n"
-		""},
-	 { "IntensityData_totalSum", _wrap_IntensityData_totalSum, METH_O, "\n"
-		"IntensityData_totalSum(IntensityData self) -> double\n"
-		"T Powerfield< T >::totalSum() const\n"
+		"Parameters:\n"
+		"-----------\n"
 		"\n"
-		"Returns sum of all values in the data structure. \n"
+		"global_index: \n"
+		"The global index of this data structure.\n"
+		"\n"
+		"Vector of corresponding bin centers \n"
 		"\n"
 		""},
-	 { "IntensityData_max", _wrap_IntensityData_max, METH_O, "\n"
-		"IntensityData_max(IntensityData self) -> double\n"
-		"T Powerfield< T >::max() const\n"
+	 { "IField_getAxisBin", _wrap_IField_getAxisBin, METH_VARARGS, "\n"
+		"IField_getAxisBin(IField self, size_t global_index, size_t i_selected_axis) -> Bin1D\n"
+		"IField_getAxisBin(IField self, size_t global_index, std::string const & axis_name) -> Bin1D\n"
+		"Bin1D IField::getAxisBin(size_t global_index, const std::string &axis_name) const\n"
 		"\n"
-		"Returns maximum value in the data structure. \n"
+		"Returns bin of selected axis for given global_index.\n"
 		"\n"
-		""},
-	 { "IntensityData_begin", _wrap_IntensityData_begin, METH_VARARGS, "\n"
-		"IntensityData_begin(IntensityData self) -> Powerfield< double >::iterator\n"
-		"IntensityData_begin(IntensityData self) -> Powerfield< double >::const_iterator\n"
-		"Powerfield< T >::const_iterator Powerfield< T >::begin() const\n"
+		"Parameters:\n"
+		"-----------\n"
 		"\n"
-		"Returns read-only iterator that points to the first element. \n"
+		"global_index: \n"
+		"The global index of this data structure.\n"
 		"\n"
-		""},
-	 { "IntensityData_end", _wrap_IntensityData_end, METH_VARARGS, "\n"
-		"IntensityData_end(IntensityData self) -> Powerfield< double >::iterator\n"
-		"IntensityData_end(IntensityData self) -> Powerfield< double >::const_iterator\n"
-		"const_iterator Powerfield< T >::end() const\n"
+		"axis_name: \n"
+		"The name of selected axis.\n"
 		"\n"
-		"Returns read-only iterator that points to the one past last element. \n"
+		"Corresponding Bin1D object \n"
 		"\n"
 		""},
-	 { "IntensityData_getAxesBinIndices", _wrap_IntensityData_getAxesBinIndices, METH_VARARGS, "\n"
-		"IntensityData_getAxesBinIndices(IntensityData self, size_t global_index) -> vector_integer_t\n"
-		"std::vector< int > Powerfield< T >::getAxesBinIndices(size_t global_index) const\n"
+	 { "IField_getAxesBinIndices", _wrap_IField_getAxesBinIndices, METH_VARARGS, "\n"
+		"IField_getAxesBinIndices(IField self, size_t global_index) -> vector_integer_t\n"
+		"std::vector< int > IField::getAxesBinIndices(size_t global_index) const\n"
 		"\n"
 		"Returns vector of axes indices for given global index\n"
 		"\n"
@@ -42835,10 +42841,10 @@ static PyMethodDef SwigMethods[] = {
 		"Vector of bin indices for all axes defined \n"
 		"\n"
 		""},
-	 { "IntensityData_getAxisBinIndex", _wrap_IntensityData_getAxisBinIndex, METH_VARARGS, "\n"
-		"IntensityData_getAxisBinIndex(IntensityData self, size_t global_index, size_t i_selected_axis) -> size_t\n"
-		"IntensityData_getAxisBinIndex(IntensityData self, size_t global_index, std::string const & axis_name) -> size_t\n"
-		"size_t Powerfield< T >::getAxisBinIndex(size_t global_index, const std::string &axis_name) const\n"
+	 { "IField_getAxisBinIndex", _wrap_IField_getAxisBinIndex, METH_VARARGS, "\n"
+		"IField_getAxisBinIndex(IField self, size_t global_index, std::string const & axis_name) -> size_t\n"
+		"IField_getAxisBinIndex(IField self, size_t global_index, size_t i_selected_axis) -> size_t\n"
+		"size_t IField::getAxisBinIndex(size_t global_index, size_t i_selected_axis) const\n"
 		"\n"
 		"Returns axis bin index for given global index\n"
 		"\n"
@@ -42848,15 +42854,15 @@ static PyMethodDef SwigMethods[] = {
 		"global_index: \n"
 		"The global index of this data structure.\n"
 		"\n"
-		"axis_name: \n"
-		"The name of selected axis.\n"
+		"i_selected_axis: \n"
+		"Serial number of selected axis.\n"
 		"\n"
 		"Corresponding bin index for selected axis \n"
 		"\n"
 		""},
-	 { "IntensityData_toGlobalIndex", _wrap_IntensityData_toGlobalIndex, METH_VARARGS, "\n"
-		"IntensityData_toGlobalIndex(IntensityData self, std::vector< unsigned int,std::allocator< unsigned int > > const & axes_indices) -> size_t\n"
-		"size_t Powerfield< T >::toGlobalIndex(const std::vector< unsigned > &axes_indices) const\n"
+	 { "IField_toGlobalIndex", _wrap_IField_toGlobalIndex, METH_VARARGS, "\n"
+		"IField_toGlobalIndex(IField self, std::vector< unsigned int,std::allocator< unsigned int > > const & axes_indices) -> size_t\n"
+		"size_t IField::toGlobalIndex(const std::vector< unsigned > &axes_indices) const\n"
 		"\n"
 		"Returns global index for specified indices of axes\n"
 		"\n"
@@ -42869,9 +42875,9 @@ static PyMethodDef SwigMethods[] = {
 		"Corresponding global index \n"
 		"\n"
 		""},
-	 { "IntensityData_findGlobalIndex", _wrap_IntensityData_findGlobalIndex, METH_VARARGS, "\n"
-		"IntensityData_findGlobalIndex(IntensityData self, vdouble1d_t coordinates) -> size_t\n"
-		"size_t Powerfield< T >::findGlobalIndex(const std::vector< double > &coordinates) const\n"
+	 { "IField_findGlobalIndex", _wrap_IField_findGlobalIndex, METH_VARARGS, "\n"
+		"IField_findGlobalIndex(IField self, vdouble1d_t coordinates) -> size_t\n"
+		"size_t IField::findGlobalIndex(const std::vector< double > &coordinates) const\n"
 		"\n"
 		"Returns global index for specified axes values\n"
 		"\n"
@@ -42884,57 +42890,84 @@ static PyMethodDef SwigMethods[] = {
 		"Closest global index \n"
 		"\n"
 		""},
-	 { "IntensityData_getAxisValue", _wrap_IntensityData_getAxisValue, METH_VARARGS, "\n"
-		"IntensityData_getAxisValue(IntensityData self, size_t global_index, size_t i_selected_axis) -> double\n"
-		"IntensityData_getAxisValue(IntensityData self, size_t global_index, std::string const & axis_name) -> double\n"
-		"double Powerfield< T >::getAxisValue(size_t global_index, const std::string &axis_name) const\n"
+	 { "IField_swigregister", IField_swigregister, METH_O, NULL},
+	 { "IField_swiginit", IField_swiginit, METH_VARARGS, NULL},
+	 { "new_IntensityData", _wrap_new_IntensityData, METH_VARARGS, "\n"
+		"IntensityData()\n"
+		"IntensityData(IAxis xAxis)\n"
+		"IntensityData(IAxis xAxis, IAxis yAxis)\n"
+		"new_IntensityData(IntensityData arg1) -> IntensityData\n"
+		"Powerfield< T >::Powerfield(Powerfield &&)\n"
 		"\n"
-		"Returns the value of selected axis for given global_index.\n"
+		""},
+	 { "delete_IntensityData", _wrap_delete_IntensityData, METH_O, "\n"
+		"delete_IntensityData(IntensityData self)\n"
+		"Powerfield< T >::~Powerfield() override\n"
 		"\n"
-		"Parameters:\n"
-		"-----------\n"
+		""},
+	 { "IntensityData_clone", _wrap_IntensityData_clone, METH_O, "\n"
+		"IntensityData_clone(IntensityData self) -> IntensityData\n"
+		"Powerfield< T > * Powerfield< T >::clone() const\n"
 		"\n"
-		"global_index: \n"
-		"The global index of this data structure.\n"
+		""},
+	 { "IntensityData_copyFrom", _wrap_IntensityData_copyFrom, METH_VARARGS, "\n"
+		"IntensityData_copyFrom(IntensityData self, IntensityData other)\n"
+		"void Powerfield< T >::copyFrom(const Powerfield< T > &other)\n"
 		"\n"
-		"axis_name: \n"
-		"The name of selected axis.\n"
+		""},
+	 { "IntensityData_meanValues", _wrap_IntensityData_meanValues, METH_O, "\n"
+		"IntensityData_meanValues(IntensityData self) -> IntensityData\n"
+		"Powerfield< double > * Powerfield< T >::meanValues() const\n"
 		"\n"
-		"corresponding bin center of selected axis \n"
+		""},
+	 { "IntensityData_addAxis", _wrap_IntensityData_addAxis, METH_VARARGS, "\n"
+		"IntensityData_addAxis(IntensityData self, IAxis new_axis)\n"
+		"IntensityData_addAxis(IntensityData self, std::string const & name, size_t size, double start, double end)\n"
+		"void Powerfield< T >::addAxis(const std::string &name, size_t size, double start, double end)\n"
 		"\n"
 		""},
-	 { "IntensityData_getAxesValues", _wrap_IntensityData_getAxesValues, METH_VARARGS, "\n"
-		"IntensityData_getAxesValues(IntensityData self, size_t global_index) -> vdouble1d_t\n"
-		"std::vector< double > Powerfield< T >::getAxesValues(size_t global_index) const\n"
+	 { "IntensityData_getAllocatedSize", _wrap_IntensityData_getAllocatedSize, METH_O, "\n"
+		"IntensityData_getAllocatedSize(IntensityData self) -> size_t\n"
+		"size_t Powerfield< T >::getAllocatedSize() const\n"
 		"\n"
-		"Returns values on all defined axes for given globalbin number\n"
+		"Returns total size of data buffer (product of bin number in every dimension). \n"
 		"\n"
-		"Parameters:\n"
-		"-----------\n"
+		""},
+	 { "IntensityData_getAllSizes", _wrap_IntensityData_getAllSizes, METH_O, "\n"
+		"IntensityData_getAllSizes(IntensityData self) -> std::vector< size_t,std::allocator< size_t > >\n"
+		"std::vector< size_t > Powerfield< T >::getAllSizes() const\n"
 		"\n"
-		"global_index: \n"
-		"The global index of this data structure.\n"
+		"Returns all sizes of its axes. \n"
 		"\n"
-		"Vector of corresponding bin centers \n"
+		""},
+	 { "IntensityData_getRawDataVector", _wrap_IntensityData_getRawDataVector, METH_O, "\n"
+		"IntensityData_getRawDataVector(IntensityData self) -> vdouble1d_t\n"
+		"std::vector< T > Powerfield< T >::getRawDataVector() const\n"
+		"\n"
+		"Returns copy of raw data vector. \n"
 		"\n"
 		""},
-	 { "IntensityData_getAxisBin", _wrap_IntensityData_getAxisBin, METH_VARARGS, "\n"
-		"IntensityData_getAxisBin(IntensityData self, size_t global_index, size_t i_selected_axis) -> Bin1D\n"
-		"IntensityData_getAxisBin(IntensityData self, size_t global_index, std::string const & axis_name) -> Bin1D\n"
-		"Bin1D Powerfield< T >::getAxisBin(size_t global_index, const std::string &axis_name) const\n"
+	 { "IntensityData_totalSum", _wrap_IntensityData_totalSum, METH_O, "\n"
+		"IntensityData_totalSum(IntensityData self) -> double\n"
+		"T Powerfield< T >::totalSum() const\n"
 		"\n"
-		"Returns bin of selected axis for given global_index.\n"
+		"Returns sum of all values in the data structure. \n"
 		"\n"
-		"Parameters:\n"
-		"-----------\n"
+		""},
+	 { "IntensityData_begin", _wrap_IntensityData_begin, METH_VARARGS, "\n"
+		"IntensityData_begin(IntensityData self) -> Powerfield< double >::iterator\n"
+		"IntensityData_begin(IntensityData self) -> Powerfield< double >::const_iterator\n"
+		"Powerfield< T >::const_iterator Powerfield< T >::begin() const\n"
 		"\n"
-		"global_index: \n"
-		"The global index of this data structure.\n"
+		"Returns read-only iterator that points to the first element. \n"
 		"\n"
-		"axis_name: \n"
-		"The name of selected axis.\n"
+		""},
+	 { "IntensityData_end", _wrap_IntensityData_end, METH_VARARGS, "\n"
+		"IntensityData_end(IntensityData self) -> Powerfield< double >::iterator\n"
+		"IntensityData_end(IntensityData self) -> Powerfield< double >::const_iterator\n"
+		"const_iterator Powerfield< T >::end() const\n"
 		"\n"
-		"Corresponding Bin1D object \n"
+		"Returns read-only iterator that points to the one past last element. \n"
 		"\n"
 		""},
 	 { "IntensityData_clear", _wrap_IntensityData_clear, METH_O, "\n"
@@ -42951,13 +42984,6 @@ static PyMethodDef SwigMethods[] = {
 		"Sets content of output data to specific value. \n"
 		"\n"
 		""},
-	 { "IntensityData_setAxisSizes", _wrap_IntensityData_setAxisSizes, METH_VARARGS, "\n"
-		"IntensityData_setAxisSizes(IntensityData self, size_t rank, int * n_dims)\n"
-		"void Powerfield< T >::setAxisSizes(size_t rank, int *n_dims)\n"
-		"\n"
-		"Adds 'rank' axes with indicated sizes. \n"
-		"\n"
-		""},
 	 { "IntensityData_setRawDataVector", _wrap_IntensityData_setRawDataVector, METH_VARARGS, "\n"
 		"IntensityData_setRawDataVector(IntensityData self, vdouble1d_t data_vector)\n"
 		"void Powerfield< T >::setRawDataVector(const std::vector< T > &data_vector)\n"
@@ -42972,10 +42998,10 @@ static PyMethodDef SwigMethods[] = {
 		"Sets new values to raw data array. \n"
 		"\n"
 		""},
-	 { "IntensityData___iadd__", _wrap_IntensityData___iadd__, METH_VARARGS, "IntensityData___iadd__(IntensityData self, IntensityData right) -> IntensityData"},
-	 { "IntensityData___isub__", _wrap_IntensityData___isub__, METH_VARARGS, "IntensityData___isub__(IntensityData self, IntensityData right) -> IntensityData"},
-	 { "IntensityData___itruediv__", _wrap_IntensityData___itruediv__, METH_VARARGS, "IntensityData___itruediv__(IntensityData self, IntensityData right) -> IntensityData"},
-	 { "IntensityData___imul__", _wrap_IntensityData___imul__, METH_VARARGS, "IntensityData___imul__(IntensityData self, IntensityData right) -> IntensityData"},
+	 { "IntensityData___iadd__", _wrap_IntensityData___iadd__, METH_VARARGS, "IntensityData___iadd__(IntensityData self, IntensityData other) -> IntensityData"},
+	 { "IntensityData___isub__", _wrap_IntensityData___isub__, METH_VARARGS, "IntensityData___isub__(IntensityData self, IntensityData other) -> IntensityData"},
+	 { "IntensityData___itruediv__", _wrap_IntensityData___itruediv__, METH_VARARGS, "IntensityData___itruediv__(IntensityData self, IntensityData other) -> IntensityData"},
+	 { "IntensityData___imul__", _wrap_IntensityData___imul__, METH_VARARGS, "IntensityData___imul__(IntensityData self, IntensityData other) -> IntensityData"},
 	 { "IntensityData_getValue", _wrap_IntensityData_getValue, METH_VARARGS, "\n"
 		"IntensityData_getValue(IntensityData self, size_t index) -> double\n"
 		"double Powerfield< double >::getValue(size_t index) const\n"
@@ -43639,6 +43665,7 @@ static PyMethodDef SwigMethods[] = {
 		""},
 	 { "new_DetectorMask", _wrap_new_DetectorMask, METH_VARARGS, "\n"
 		"DetectorMask()\n"
+		"DetectorMask(IAxis xAxis, IAxis yAxis)\n"
 		"new_DetectorMask(DetectorMask other) -> DetectorMask\n"
 		"DetectorMask::DetectorMask(const DetectorMask &other)\n"
 		"\n"
@@ -44486,10 +44513,6 @@ static PyMethodDef SwigMethods[] = {
 		"\n"
 		""},
 	 { "IHistogram_createHistogram", _wrap_IHistogram_createHistogram, METH_O, "IHistogram_createHistogram(IntensityData source) -> IHistogram"},
-	 { "IHistogram_createFrom", _wrap_IHistogram_createFrom, METH_VARARGS, "\n"
-		"IHistogram_createFrom(std::string const & filename) -> IHistogram\n"
-		"IHistogram_createFrom(vdouble2d_t data) -> IHistogram\n"
-		""},
 	 { "IHistogram_createPowerfield", _wrap_IHistogram_createPowerfield, METH_VARARGS, "\n"
 		"IHistogram_createPowerfield(IHistogram self, IHistogram::DataType dataType=DataType::INTEGRAL) -> IntensityData\n"
 		"Powerfield< double > * IHistogram::createPowerfield(DataType dataType=DataType::INTEGRAL) const\n"
@@ -44617,18 +44640,19 @@ static PyMethodDef SwigMethods[] = {
 		"Histogram2D(int nbinsx, double xlow, double xup, int nbinsy, double ylow, double yup)\n"
 		"Histogram2D(int nbinsx, vdouble1d_t xbins, int nbinsy, vdouble1d_t ybins)\n"
 		"Histogram2D(IAxis axis_x, IAxis axis_y)\n"
-		"Histogram2D(IntensityData data)\n"
-		"new_Histogram2D(vdouble2d_t data) -> Histogram2D\n"
-		"Histogram2D::Histogram2D(std::vector< std::vector< double >> data)\n"
+		"new_Histogram2D(IntensityData data) -> Histogram2D\n"
+		"Histogram2D::Histogram2D(const Powerfield< double > &data)\n"
 		"\n"
-		"Constructor for 2D histograms from numpy array (thanks to swig) \n"
+		"Constructor for 2D histograms from basic  Powerfield object. \n"
 		"\n"
 		""},
 	 { "Histogram2D_clone", _wrap_Histogram2D_clone, METH_O, "\n"
 		"Histogram2D_clone(Histogram2D self) -> Histogram2D\n"
 		"Histogram2D * Histogram2D::clone() const override\n"
 		"\n"
-		"Returns clone of other histogram. \n"
+		"Constructor for 2D histograms from numpy array (thanks to swig)\n"
+		"\n"
+		"Returns clone of other histogram \n"
 		"\n"
 		""},
 	 { "Histogram2D_rank", _wrap_Histogram2D_rank, METH_O, "\n"
@@ -44781,11 +44805,6 @@ static PyMethodDef SwigMethods[] = {
 		"size_t SimulationResult::size() const\n"
 		"\n"
 		""},
-	 { "SimulationResult_max", _wrap_SimulationResult_max, METH_O, "\n"
-		"SimulationResult_max(SimulationResult self) -> double\n"
-		"double SimulationResult::max() const\n"
-		"\n"
-		""},
 	 { "SimulationResult_empty", _wrap_SimulationResult_empty, METH_O, "\n"
 		"SimulationResult_empty(SimulationResult self) -> bool\n"
 		"bool SimulationResult::empty() const\n"
@@ -44850,6 +44869,9 @@ static void *_p_SphericalDetectorTo_p_IDetector(void *x, int *SWIGUNUSEDPARM(new
 static void *_p_IDetector2DTo_p_IDetector(void *x, int *SWIGUNUSEDPARM(newmemory)) {
     return (void *)((IDetector *)  ((IDetector2D *) x));
 }
+static void *_p_PowerfieldT_double_tTo_p_IField(void *x, int *SWIGUNUSEDPARM(newmemory)) {
+    return (void *)((IField *)  ((Powerfield< double > *) x));
+}
 static void *_p_PolygonTo_p_IShape2D(void *x, int *SWIGUNUSEDPARM(newmemory)) {
     return (void *)((IShape2D *)  ((Polygon *) x));
 }
@@ -44989,6 +45011,7 @@ static swig_type_info _swigt__p_ICoordSystem = {"_p_ICoordSystem", "ICoordSystem
 static swig_type_info _swigt__p_IDetector = {"_p_IDetector", "IDetector *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_IDetector2D = {"_p_IDetector2D", "IDetector2D *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_IDetectorResolution = {"_p_IDetectorResolution", "IDetectorResolution *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_IField = {"_p_IField", "IField *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_IFootprintFactor = {"_p_IFootprintFactor", "IFootprintFactor *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_IHistogram = {"_p_IHistogram", "IHistogram *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_INode = {"_p_INode", "INode *", 0, 0, (void*)0, 0};
@@ -45056,6 +45079,7 @@ static swig_type_info _swigt__p_std__lessT_std__string_t = {"_p_std__lessT_std__
 static swig_type_info _swigt__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t = {"_p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t", "std::map< std::string,double,std::less< std::string >,std::allocator< std::pair< std::string const,double > > > *|std::map< std::string,double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__pairT_double_double_t = {"_p_std__pairT_double_double_t", "std::pair< double,double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t = {"_p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t", "std::vector< AxisInfo,std::allocator< AxisInfo > > *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t = {"_p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t", "std::vector< IAxis *,std::allocator< IAxis * > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t = {"_p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t", "std::vector< INode const *,std::allocator< INode const * > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t = {"_p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t", "std::vector< ParaMeta,std::allocator< ParaMeta > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t = {"_p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t", "std::vector< Vec3< double > > *|std::vector< Vec3< double >,std::allocator< Vec3< double > > > *", 0, 0, (void*)0, 0};
@@ -45097,6 +45121,7 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_IDetector,
   &_swigt__p_IDetector2D,
   &_swigt__p_IDetectorResolution,
+  &_swigt__p_IField,
   &_swigt__p_IFootprintFactor,
   &_swigt__p_IHistogram,
   &_swigt__p_INode,
@@ -45164,6 +45189,7 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t,
   &_swigt__p_std__pairT_double_double_t,
   &_swigt__p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t,
+  &_swigt__p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t,
   &_swigt__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t,
   &_swigt__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t,
   &_swigt__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t,
@@ -45205,6 +45231,7 @@ static swig_cast_info _swigc__p_ICoordSystem[] = {  {&_swigt__p_ICoordSystem, 0,
 static swig_cast_info _swigc__p_IDetector[] = {  {&_swigt__p_IDetector, 0, 0, 0},  {&_swigt__p_RectangularDetector, _p_RectangularDetectorTo_p_IDetector, 0, 0},  {&_swigt__p_SphericalDetector, _p_SphericalDetectorTo_p_IDetector, 0, 0},  {&_swigt__p_IDetector2D, _p_IDetector2DTo_p_IDetector, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_IDetector2D[] = {  {&_swigt__p_RectangularDetector, _p_RectangularDetectorTo_p_IDetector2D, 0, 0},  {&_swigt__p_SphericalDetector, _p_SphericalDetectorTo_p_IDetector2D, 0, 0},  {&_swigt__p_IDetector2D, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_IDetectorResolution[] = {  {&_swigt__p_IDetectorResolution, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_IField[] = {  {&_swigt__p_PowerfieldT_double_t, _p_PowerfieldT_double_tTo_p_IField, 0, 0},  {&_swigt__p_IField, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_IFootprintFactor[] = {  {&_swigt__p_FootprintSquare, _p_FootprintSquareTo_p_IFootprintFactor, 0, 0},  {&_swigt__p_IFootprintFactor, 0, 0, 0},  {&_swigt__p_FootprintGauss, _p_FootprintGaussTo_p_IFootprintFactor, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_IHistogram[] = {  {&_swigt__p_IHistogram, 0, 0, 0},  {&_swigt__p_Histogram2D, _p_Histogram2DTo_p_IHistogram, 0, 0},  {&_swigt__p_Histogram1D, _p_Histogram1DTo_p_IHistogram, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_INode[] = {  {&_swigt__p_INode, 0, 0, 0},  {&_swigt__p_FootprintSquare, _p_FootprintSquareTo_p_INode, 0, 0},  {&_swigt__p_IDetectorResolution, _p_IDetectorResolutionTo_p_INode, 0, 0},  {&_swigt__p_Instrument, _p_InstrumentTo_p_INode, 0, 0},  {&_swigt__p_IFootprintFactor, _p_IFootprintFactorTo_p_INode, 0, 0},  {&_swigt__p_Beam, _p_BeamTo_p_INode, 0, 0},  {&_swigt__p_IResolutionFunction2D, _p_IResolutionFunction2DTo_p_INode, 0, 0},  {&_swigt__p_ResolutionFunction2DGaussian, _p_ResolutionFunction2DGaussianTo_p_INode, 0, 0},  {&_swigt__p_IDetector, _p_IDetectorTo_p_INode, 0, 0},  {&_swigt__p_RectangularDetector, _p_RectangularDetectorTo_p_INode, 0, 0},  {&_swigt__p_SphericalDetector, _p_SphericalDetectorTo_p_INode, 0, 0},  {&_swigt__p_FootprintGauss, _p_FootprintGaussTo_p_INode, 0, 0},  {&_swigt__p_IDetector2D, _p_IDetector2DTo_p_INode, 0, 0},{0, 0, 0, 0}};
@@ -45272,6 +45299,7 @@ static swig_cast_info _swigc__p_std__lessT_std__string_t[] = {  {&_swigt__p_std_
 static swig_cast_info _swigc__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t[] = {  {&_swigt__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__pairT_double_double_t[] = {  {&_swigt__p_std__pairT_double_double_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t[] = {  {&_swigt__p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t[] = {  {&_swigt__p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t[] = {  {&_swigt__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t[] = {  {&_swigt__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t[] = {  {&_swigt__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
@@ -45313,6 +45341,7 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_IDetector,
   _swigc__p_IDetector2D,
   _swigc__p_IDetectorResolution,
+  _swigc__p_IField,
   _swigc__p_IFootprintFactor,
   _swigc__p_IHistogram,
   _swigc__p_INode,
@@ -45380,6 +45409,7 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t,
   _swigc__p_std__pairT_double_double_t,
   _swigc__p_std__vectorT_AxisInfo_std__allocatorT_AxisInfo_t_t,
+  _swigc__p_std__vectorT_IAxis_p_std__allocatorT_IAxis_p_t_t,
   _swigc__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t,
   _swigc__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t,
   _swigc__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t,