diff --git a/Core/PythonAPI/src/IntensityDataFunctions.pypp.cpp b/Core/PythonAPI/src/IntensityDataFunctions.pypp.cpp
index e494ecaf3b982016961cb9294d5a4d1422d21e6d..35df429022b2162e2780f6c519a0515cf57cf32e 100644
--- a/Core/PythonAPI/src/IntensityDataFunctions.pypp.cpp
+++ b/Core/PythonAPI/src/IntensityDataFunctions.pypp.cpp
@@ -20,6 +20,26 @@ void register_IntensityDataFunctions_class(){
         typedef bp::class_< IntensityDataFunctions > IntensityDataFunctions_exposer_t;
         IntensityDataFunctions_exposer_t IntensityDataFunctions_exposer = IntensityDataFunctions_exposer_t( "IntensityDataFunctions" );
         bp::scope IntensityDataFunctions_scope( IntensityDataFunctions_exposer );
+        { //::IntensityDataFunctions::addEllipticMask
+        
+            typedef void ( *addEllipticMask_function_type )( ::OutputData< double > &,double,double,double,double,bool );
+            
+            IntensityDataFunctions_exposer.def( 
+                "addEllipticMask"
+                , addEllipticMask_function_type( &::IntensityDataFunctions::addEllipticMask )
+                , ( bp::arg("data"), bp::arg("xc"), bp::arg("yc"), bp::arg("rx"), bp::arg("ry"), bp::arg("invert_flag")=(bool)(false) ) );
+        
+        }
+        { //::IntensityDataFunctions::addRectangularMask
+        
+            typedef void ( *addRectangularMask_function_type )( ::OutputData< double > &,double,double,double,double,bool );
+            
+            IntensityDataFunctions_exposer.def( 
+                "addRectangularMask"
+                , addRectangularMask_function_type( &::IntensityDataFunctions::addRectangularMask )
+                , ( bp::arg("data"), bp::arg("x1"), bp::arg("y1"), bp::arg("x2"), bp::arg("y2"), bp::arg("invert_flag")=(bool)(false) ) );
+        
+        }
         { //::IntensityDataFunctions::createClippedDataSet
         
             typedef ::OutputData< double > * ( *createClippedDataSet_function_type )( ::OutputData< double > const &,double,double,double,double );
@@ -43,24 +63,26 @@ void register_IntensityDataFunctions_class(){
         }
         { //::IntensityDataFunctions::setEllipticMask
         
-            typedef void ( *setEllipticMask_function_type )( ::OutputData< double > &,double,double,double,double );
+            typedef void ( *setEllipticMask_function_type )( ::OutputData< double > &,double,double,double,double,bool );
             
             IntensityDataFunctions_exposer.def( 
                 "setEllipticMask"
                 , setEllipticMask_function_type( &::IntensityDataFunctions::setEllipticMask )
-                , ( bp::arg("data"), bp::arg("xc"), bp::arg("yc"), bp::arg("rx"), bp::arg("ry") ) );
+                , ( bp::arg("data"), bp::arg("xc"), bp::arg("yc"), bp::arg("rx"), bp::arg("ry"), bp::arg("invert_flag")=(bool)(false) ) );
         
         }
         { //::IntensityDataFunctions::setRectangularMask
         
-            typedef void ( *setRectangularMask_function_type )( ::OutputData< double > &,double,double,double,double );
+            typedef void ( *setRectangularMask_function_type )( ::OutputData< double > &,double,double,double,double,bool );
             
             IntensityDataFunctions_exposer.def( 
                 "setRectangularMask"
                 , setRectangularMask_function_type( &::IntensityDataFunctions::setRectangularMask )
-                , ( bp::arg("data"), bp::arg("x1"), bp::arg("y1"), bp::arg("x2"), bp::arg("y2") ) );
+                , ( bp::arg("data"), bp::arg("x1"), bp::arg("y1"), bp::arg("x2"), bp::arg("y2"), bp::arg("invert_flag")=(bool)(false) ) );
         
         }
+        IntensityDataFunctions_exposer.staticmethod( "addEllipticMask" );
+        IntensityDataFunctions_exposer.staticmethod( "addRectangularMask" );
         IntensityDataFunctions_exposer.staticmethod( "createClippedDataSet" );
         IntensityDataFunctions_exposer.staticmethod( "getRelativeDifference" );
         IntensityDataFunctions_exposer.staticmethod( "setEllipticMask" );
diff --git a/Core/PythonAPI/src/PythonModule.cpp b/Core/PythonAPI/src/PythonModule.cpp
index f8da0464b555b8049ec71c3ac703561ca93ac0b6..0d80ba940cdeae9e0b9345cdfcc2e631533fe1a4 100644
--- a/Core/PythonAPI/src/PythonModule.cpp
+++ b/Core/PythonAPI/src/PythonModule.cpp
@@ -10,129 +10,129 @@ GCC_DIAG_ON(missing-field-initializers)
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #define PY_ARRAY_UNIQUE_SYMBOL BORNAGAIN_PYTHONAPI_ARRAY
 #include "numpy/arrayobject.h"
-#include "FormFactorFullSpheroid.pypp.h"
-#include "DistributionGate.pypp.h"
-#include "SimpleSelectionRule.pypp.h"
-#include "RealParameterWrapper.pypp.h"
-#include "vdouble1d_t.pypp.h"
+#include "FormFactorFullSphere.pypp.h"
+#include "Lattice.pypp.h"
+#include "FormFactorInfLongBox.pypp.h"
+#include "Layer.pypp.h"
+#include "FormFactorHemiEllipsoid.pypp.h"
+#include "FormFactorPrism3.pypp.h"
 #include "SimulationParameters.pypp.h"
-#include "Transform3D.pypp.h"
-#include "ThreadInfo.pypp.h"
-#include "InterferenceFunction2DLattice.pypp.h"
-#include "LayerInterface.pypp.h"
+#include "FormFactorInfLongRipple1.pypp.h"
+#include "FTDistribution1DVoigt.pypp.h"
+#include "FormFactorRipple1.pypp.h"
+#include "ParticleInfo.pypp.h"
+#include "ICompositeSample.pypp.h"
+#include "OffSpecSimulation.pypp.h"
+#include "IResolutionFunction2D.pypp.h"
+#include "vector_kvector_t.pypp.h"
+#include "FTDistribution1DGate.pypp.h"
+#include "LatticeBasis.pypp.h"
+#include "FormFactorRipple2.pypp.h"
+#include "FTDistribution2DVoigt.pypp.h"
+#include "IFormFactorBorn.pypp.h"
+#include "FormFactorPyramid.pypp.h"
+#include "ISample.pypp.h"
 #include "ILayout.pypp.h"
-#include "FormFactorCone6.pypp.h"
 #include "FormFactorTetrahedron.pypp.h"
-#include "FTDistribution1DCosine.pypp.h"
-#include "FTDistribution1DTriangle.pypp.h"
-#include "FormFactorWeighted.pypp.h"
-#include "InterferenceFunctionRadialParaCrystal.pypp.h"
-#include "DistributionGaussian.pypp.h"
-#include "IDetectorResolution.pypp.h"
-#include "FormFactorCylinder.pypp.h"
-#include "Crystal.pypp.h"
-#include "FTDistribution1DCauchy.pypp.h"
-#include "IFormFactorBorn.pypp.h"
+#include "Simulation.pypp.h"
+#include "ParticleLayout.pypp.h"
+#include "FTDistribution2DCone.pypp.h"
 #include "FormFactorEllipsoidalCylinder.pypp.h"
-#include "InterferenceFunctionNone.pypp.h"
-#include "FTDistribution2DGate.pypp.h"
-#include "vector_kvector_t.pypp.h"
-#include "FormFactorTruncatedSpheroid.pypp.h"
-#include "Particle.pypp.h"
-#include "FormFactorTrivial.pypp.h"
-#include "ConstKBinAxis.pypp.h"
-#include "FTDistribution2DCauchy.pypp.h"
-#include "FormFactorCrystal.pypp.h"
-#include "ParticleDistribution.pypp.h"
-#include "vector_longinteger_t.pypp.h"
-#include "ResolutionFunction2DGaussian.pypp.h"
+#include "vector_IFormFactorPtr_t.pypp.h"
 #include "FTDistribution1DGauss.pypp.h"
-#include "FTDistribution1DGate.pypp.h"
-#include "FormFactorAnisoPyramid.pypp.h"
-#include "FixedBinAxis.pypp.h"
-#include "MultiLayer.pypp.h"
-#include "IFormFactor.pypp.h"
-#include "kvector_t.pypp.h"
-#include "FormFactorSphereUniformRadius.pypp.h"
-#include "OffSpecSimulation.pypp.h"
-#include "FormFactorRipple1.pypp.h"
-#include "Simulation.pypp.h"
+#include "HomogeneousMaterial.pypp.h"
+#include "cvector_t.pypp.h"
+#include "FTDistribution2DGauss.pypp.h"
 #include "IObservable.pypp.h"
+#include "FormFactorBox.pypp.h"
+#include "DistributionLogNormal.pypp.h"
+#include "FormFactorCylinder.pypp.h"
+#include "Detector.pypp.h"
+#include "IObserver.pypp.h"
+#include "ResolutionFunction2DGaussian.pypp.h"
+#include "IParameterized.pypp.h"
+#include "FormFactorGauss.pypp.h"
+#include "FormFactorDecoratorDebyeWaller.pypp.h"
+#include "IDetectorResolution.pypp.h"
+#include "FormFactorWeighted.pypp.h"
+#include "InterferenceFunction2DLattice.pypp.h"
+#include "RealParameterWrapper.pypp.h"
+#include "LayerInterface.pypp.h"
+#include "ParticleCoreShell.pypp.h"
+#include "InterferenceFunction1DLattice.pypp.h"
+#include "FormFactorCrystal.pypp.h"
+#include "Instrument.pypp.h"
+#include "PythonInterface_global_variables.pypp.h"
+#include "ParticleDistribution.pypp.h"
+#include "IFTDistribution1D.pypp.h"
+#include "FormFactorTrivial.pypp.h"
+#include "DistributionGaussian.pypp.h"
+#include "Lattice1DIFParameters.pypp.h"
+#include "DistributionCosine.pypp.h"
+#include "IClusteredParticles.pypp.h"
 #include "FormFactorLorentz.pypp.h"
-#include "SpecularSimulation.pypp.h"
-#include "ISelectionRule.pypp.h"
-#include "FormFactorRipple2.pypp.h"
-#include "LayerRoughness.pypp.h"
 #include "Bin1DCVector.pypp.h"
-#include "FormFactorSphereGaussianRadius.pypp.h"
-#include "ParameterPool.pypp.h"
-#include "FormFactorPrism3.pypp.h"
+#include "InterferenceFunctionNone.pypp.h"
+#include "IFormFactorDecorator.pypp.h"
+#include "Bin1D.pypp.h"
+#include "ISampleBuilder.pypp.h"
+#include "IntensityDataIOFactory.pypp.h"
+#include "FTDistribution2DGate.pypp.h"
+#include "FormFactorInfLongRipple2.pypp.h"
 #include "IMaterial.pypp.h"
-#include "FTDistribution1DVoigt.pypp.h"
-#include "IntensityDataFunctions.pypp.h"
-#include "FormFactorPrism6.pypp.h"
-#include "IClusteredParticles.pypp.h"
-#include "VariableBinAxis.pypp.h"
+#include "InterferenceFunction2DParaCrystal.pypp.h"
+#include "Particle.pypp.h"
+#include "FormFactorTruncatedSpheroid.pypp.h"
+#include "MesoCrystal.pypp.h"
+#include "FTDistribution2DCauchy.pypp.h"
+#include "IInterferenceFunction.pypp.h"
+#include "FTDistribution1DCosine.pypp.h"
 #include "IParticle.pypp.h"
-#include "DistributionCosine.pypp.h"
-#include "FormFactorHemiEllipsoid.pypp.h"
+#include "FTDistribution1DTriangle.pypp.h"
+#include "FormFactorPrism6.pypp.h"
+#include "HomogeneousMagneticMaterial.pypp.h"
+#include "FormFactorFullSpheroid.pypp.h"
+#include "Transform3D.pypp.h"
+#include "IntensityData.pypp.h"
+#include "DistributionGate.pypp.h"
+#include "ThreadInfo.pypp.h"
 #include "IAxis.pypp.h"
-#include "vector_integer_t.pypp.h"
-#include "IntensityDataIOFactory.pypp.h"
-#include "ParameterDistribution.pypp.h"
-#include "Layer.pypp.h"
-#include "FormFactorPyramid.pypp.h"
 #include "CustomBinAxis.pypp.h"
-#include "FTDistribution2DCone.pypp.h"
-#include "IFTDistribution1D.pypp.h"
-#include "DistributionLorentz.pypp.h"
-#include "IDistribution1D.pypp.h"
-#include "HomogeneousMagneticMaterial.pypp.h"
-#include "FormFactorCuboctahedron.pypp.h"
-#include "cvector_t.pypp.h"
+#include "FormFactorCone6.pypp.h"
+#include "ICloneable.pypp.h"
 #include "PythonInterface_free_functions.pypp.h"
+#include "vector_longinteger_t.pypp.h"
 #include "FormFactorSphereLogNormalRadius.pypp.h"
-#include "FormFactorInfLongRipple1.pypp.h"
-#include "IResolutionFunction2D.pypp.h"
-#include "vector_IFormFactorPtr_t.pypp.h"
-#include "FormFactorFullSphere.pypp.h"
-#include "ParticleLayout.pypp.h"
-#include "FormFactorBox.pypp.h"
-#include "IParameterized.pypp.h"
-#include "Lattice2DIFParameters.pypp.h"
-#include "IFormFactorDecorator.pypp.h"
-#include "InterferenceFunction1DLattice.pypp.h"
-#include "ISample.pypp.h"
-#include "ISampleBuilder.pypp.h"
-#include "PythonInterface_global_variables.pypp.h"
-#include "Beam.pypp.h"
-#include "HomogeneousMaterial.pypp.h"
-#include "ICloneable.pypp.h"
-#include "ParticleCoreShell.pypp.h"
-#include "FormFactorDecoratorDebyeWaller.pypp.h"
-#include "MesoCrystal.pypp.h"
-#include "Lattice1DIFParameters.pypp.h"
-#include "IObserver.pypp.h"
-#include "IntensityData.pypp.h"
-#include "Lattice.pypp.h"
-#include "IInterferenceFunction.pypp.h"
-#include "ParticleInfo.pypp.h"
-#include "Instrument.pypp.h"
-#include "FormFactorInfLongBox.pypp.h"
-#include "FormFactorCone.pypp.h"
-#include "FTDistribution2DGauss.pypp.h"
+#include "FormFactorSphereGaussianRadius.pypp.h"
 #include "FormFactorTruncatedSphere.pypp.h"
-#include "FTDistribution2DVoigt.pypp.h"
-#include "FormFactorGauss.pypp.h"
-#include "InterferenceFunction2DParaCrystal.pypp.h"
-#include "Detector.pypp.h"
-#include "FormFactorInfLongRipple2.pypp.h"
-#include "LatticeBasis.pypp.h"
-#include "ICompositeSample.pypp.h"
-#include "Bin1D.pypp.h"
+#include "IntensityDataFunctions.pypp.h"
+#include "DistributionLorentz.pypp.h"
+#include "ConstKBinAxis.pypp.h"
+#include "ParameterDistribution.pypp.h"
+#include "ParameterPool.pypp.h"
 #include "vector_complex_t.pypp.h"
-#include "DistributionLogNormal.pypp.h"
+#include "SpecularSimulation.pypp.h"
+#include "FormFactorAnisoPyramid.pypp.h"
+#include "FormFactorCuboctahedron.pypp.h"
+#include "MultiLayer.pypp.h"
+#include "FormFactorCone.pypp.h"
+#include "IDistribution1D.pypp.h"
+#include "LayerRoughness.pypp.h"
+#include "VariableBinAxis.pypp.h"
+#include "SimpleSelectionRule.pypp.h"
+#include "FixedBinAxis.pypp.h"
+#include "InterferenceFunctionRadialParaCrystal.pypp.h"
+#include "Lattice2DIFParameters.pypp.h"
+#include "IFormFactor.pypp.h"
+#include "vdouble1d_t.pypp.h"
 #include "IFTDistribution2D.pypp.h"
+#include "Beam.pypp.h"
+#include "FormFactorSphereUniformRadius.pypp.h"
+#include "FTDistribution1DCauchy.pypp.h"
+#include "kvector_t.pypp.h"
+#include "Crystal.pypp.h"
+#include "ISelectionRule.pypp.h"
+#include "vector_integer_t.pypp.h"
 #include "__call_policies.pypp.hpp"
 #include "__convenience.pypp.hpp"
 #include "__call_policies.pypp.hpp"
diff --git a/Core/Tools/inc/IntensityDataFunctions.h b/Core/Tools/inc/IntensityDataFunctions.h
index 3aff9da9496260a8b6c2532082fecb55e21c225e..2e075d6a6ff04d6cab07b19f812a5179c517ebf5 100644
--- a/Core/Tools/inc/IntensityDataFunctions.h
+++ b/Core/Tools/inc/IntensityDataFunctions.h
@@ -27,15 +27,31 @@ class  BA_CORE_API_ IntensityDataFunctions
 {
 public:
 
-    //! Sets rectangular mask to IntensityData to exclude all points outside the
-    //! mask from analysis
+    //! @brief Sets rectangular mask to IntensityData to exclude all points outside the
+    //! mask from analysis. If masks alreay exists, they will be replaced.
+    //! @param data Intensity data object to set the mask
+    //! @param x1 x-cordinate of lower left corner of the rectangle
+    //! @param y1 y-cordinate of lower left corner of the rectangle
+    //! @param x2 x-cordinate of top right corner of the rectangle
+    //! @param y2 y-cordinate of top right corner of the rectangle
+    //! @param invert_flag if true the area will be included in the analysis
     static void setRectangularMask(OutputData<double>& data,
-        double x1, double y1, double x2, double y2);
+        double x1, double y1, double x2, double y2, bool invert_flag = false);
+
+    //! @brief Adds rectangular mask to IntensityData to exclude all points outside the
+    //! mask from analysis
+    static void addRectangularMask(OutputData<double>& data,
+        double x1, double y1, double x2, double y2, bool invert_flag = false);
 
     //! Sets elliptic mask to IntensityData to exclude all points outside the
     //! mask from analysis
     static void setEllipticMask(OutputData<double>& data,
-        double xc, double yc, double rx, double ry);
+        double xc, double yc, double rx, double ry, bool invert_flag = false);
+
+    //! Adds elliptic mask to IntensityData to exclude all points outside the
+    //! mask from analysis
+    static void addEllipticMask(OutputData<double>& data,
+        double xc, double yc, double rx, double ry, bool invert_flag = false);
 
     //! Returns relative difference between two data sets
     //! sum(result[i] - reference[i])/reference[i])
diff --git a/Core/Tools/inc/OutputDataFunctions.h b/Core/Tools/inc/OutputDataFunctions.h
index ebfd7ee230dc6753806b9ad7bfa64656b9637179..57c06fe46209a9bdf444e7d3a14b53d83d5fc4f0 100644
--- a/Core/Tools/inc/OutputDataFunctions.h
+++ b/Core/Tools/inc/OutputDataFunctions.h
@@ -81,18 +81,18 @@ namespace OutputDataFunctions
     //! limits
     BA_CORE_API_ Mask *CreateRectangularMask(
         const OutputData<double>& data,
-        const double *minima, const double *maxima);
+        const double *minima, const double *maxima, bool invert_flag = false);
     BA_CORE_API_ Mask *CreateRectangularMask(
         const OutputData<double>& data,
-        double x1, double y1, double x2, double y2);
+        double x1, double y1, double x2, double y2, bool invert_flag = false);
 
     //! Creates a elliptic mask based on the given OutputData object and limits
     BA_CORE_API_ Mask *CreateEllipticMask(
         const OutputData<double>& data,
-        const double *center, const double *radii);
+        const double *center, const double *radii, bool invert_flag = false);
     BA_CORE_API_ Mask *CreateEllipticMask(
         const OutputData<double>& data,
-        double xc, double yc, double rx, double ry);
+        double xc, double yc, double rx, double ry, bool invert_flag = false);
 
     // compare result with reference and return the difference
 //    BA_CORE_API_ double GetDifference(const OutputData<double> &result,
diff --git a/Core/Tools/src/FixedBinAxis.cpp b/Core/Tools/src/FixedBinAxis.cpp
index 11a0d4a79155fa17787a50425252bbe9891c0c8c..6655e98438ba867cfb21ab77acc4a9db3ba04031 100644
--- a/Core/Tools/src/FixedBinAxis.cpp
+++ b/Core/Tools/src/FixedBinAxis.cpp
@@ -61,11 +61,16 @@ double FixedBinAxis::getMax() const
 
 size_t FixedBinAxis::findClosestIndex(double value) const
 {
-    if (value < getMin() || value >= getMax()) {
-        std::ostringstream ostr;
-        ostr << "FixedBinAxis::findClosestIndex() -> Error! Given value not in any bin. ";
-        ostr << "value:" << value << " name:" << getName() << " min:" << getMin() << " max:" << getMax();
-        throw OutOfBoundsException(ostr.str());
+//    if (value < getMin() || value >= getMax()) {
+//        std::ostringstream ostr;
+//        ostr << "FixedBinAxis::findClosestIndex() -> Error! Given value not in any bin. ";
+//        ostr << "value:" << value << " name:" << getName() << " min:" << getMin() << " max:" << getMax();
+//        throw OutOfBoundsException(ostr.str());
+//    }
+    if( value < getMin()) {
+        return 0;
+    } else if(value >= getMax()) {
+        return m_nbins-1;
     }
 
     double step = (m_end - m_start)/m_nbins;
diff --git a/Core/Tools/src/IntensityDataFunctions.cpp b/Core/Tools/src/IntensityDataFunctions.cpp
index e688bf9428dddc8794b86bc41a5d1990a1af2c36..bf55d5b773ce7d4ab9bae08740028b46dbd24302 100644
--- a/Core/Tools/src/IntensityDataFunctions.cpp
+++ b/Core/Tools/src/IntensityDataFunctions.cpp
@@ -4,20 +4,32 @@
 #include <boost/scoped_ptr.hpp>
 
 void IntensityDataFunctions::setRectangularMask(OutputData<double>& data,
-    double x1, double y1, double x2, double y2)
+    double x1, double y1, double x2, double y2, bool invert_flag)
 {
-    Mask *mask1 = OutputDataFunctions::CreateRectangularMask(data, x1, y1, x2, y2);
+    Mask *mask1 = OutputDataFunctions::CreateRectangularMask(data, x1, y1, x2, y2, invert_flag);
     data.setMask(*mask1);
 }
 
+void IntensityDataFunctions::addRectangularMask(OutputData<double> &data, double x1, double y1, double x2, double y2, bool invert_flag)
+{
+    Mask *mask1 = OutputDataFunctions::CreateRectangularMask(data, x1, y1, x2, y2, invert_flag);
+    data.addMask(*mask1);
+}
+
 
 void IntensityDataFunctions::setEllipticMask(OutputData<double>& data,
-    double xc, double yc, double rx, double ry)
+    double xc, double yc, double rx, double ry, bool invert_flag)
 {
-    Mask *mask1 = OutputDataFunctions::CreateEllipticMask(data, xc, yc, rx, ry);
+    Mask *mask1 = OutputDataFunctions::CreateEllipticMask(data, xc, yc, rx, ry, invert_flag);
     data.setMask(*mask1);
 }
 
+void IntensityDataFunctions::addEllipticMask(OutputData<double> &data, double xc, double yc, double rx, double ry, bool invert_flag)
+{
+    Mask *mask1 = OutputDataFunctions::CreateEllipticMask(data, xc, yc, rx, ry, invert_flag);
+    data.addMask(*mask1);
+}
+
 double IntensityDataFunctions::getRelativeDifference(
         const OutputData<double> &result, const OutputData<double> &reference)
 {
diff --git a/Core/Tools/src/OutputDataFunctions.cpp b/Core/Tools/src/OutputDataFunctions.cpp
index 37ae3b50ffb7237f5f20281d14b139e3acc57198..cef497f0af01e94a7ebc961891849ab22059280d 100644
--- a/Core/Tools/src/OutputDataFunctions.cpp
+++ b/Core/Tools/src/OutputDataFunctions.cpp
@@ -398,7 +398,7 @@ void OutputDataFunctions::applyFunction(
 
 Mask* OutputDataFunctions::CreateRectangularMask(
     const OutputData<double>& data,
-    const double* minima, const double* maxima)
+    const double* minima, const double* maxima, bool invert_flag)
 {
     size_t rank = data.getRank();
     int *minima_i = new int[rank];
@@ -412,7 +412,7 @@ Mask* OutputDataFunctions::CreateRectangularMask(
     }
     MaskCoordinateRectangleFunction *p_rectangle_function =
             new MaskCoordinateRectangleFunction(rank, minima_i, maxima_i);
-    p_rectangle_function->setInvertFlag(true);
+    p_rectangle_function->setInvertFlag(invert_flag);
     delete[] minima_i;
     delete[] maxima_i;
     MaskCoordinates *p_result = new MaskCoordinates(rank, dims_i);
@@ -422,21 +422,20 @@ Mask* OutputDataFunctions::CreateRectangularMask(
 }
 
 
-Mask* OutputDataFunctions::CreateRectangularMask(
-    const OutputData<double>& data, double x1, double y1, double x2, double y2)
+Mask* OutputDataFunctions::CreateRectangularMask(const OutputData<double>& data, double x1, double y1, double x2, double y2, bool invert_flag)
 {
     if(data.getRank() != 2) throw LogicErrorException(
             "OutputDataFunctions::CreateRectangularMask2D()"
             " -> Error! Number of dimensions should be 2");
     const double minima[2]={x1, y1};
     const double maxima[2]={x2, y2};
-    return OutputDataFunctions::CreateRectangularMask(data, minima, maxima);
+    return OutputDataFunctions::CreateRectangularMask(data, minima, maxima, invert_flag);
 }
 
 
 Mask* OutputDataFunctions::CreateEllipticMask(
     const OutputData<double>& data,
-    const double* center, const double* radii)
+    const double* center, const double* radii, bool invert_flag)
 {
     size_t rank = data.getRank();
     int *center_i = new int[rank];
@@ -447,12 +446,15 @@ Mask* OutputDataFunctions::CreateEllipticMask(
         center_i[i] = (int)p_axis->findClosestIndex(center[i]);
         int lower_index = (int)p_axis->findClosestIndex(
                 (*p_axis)[center_i[i]] - radii[i] );
+//        int lower_index = (int)p_axis->findClosestIndex(center[i] - radii[i]);
+
         radii_i[i] = center_i[i] - lower_index;
+        std::cout << "XXX " << i << " center:" << center[i] << " center_i" << center_i[i] << " lower_index:" << lower_index << " radii_i[i]:" << radii_i[i]<< std::endl;
         dims_i[i] = (int)p_axis->getSize();
     }
     MaskCoordinateEllipseFunction *p_ellipse_function =
             new MaskCoordinateEllipseFunction(rank, center_i, radii_i);
-    p_ellipse_function->setInvertFlag(true);
+    p_ellipse_function->setInvertFlag(invert_flag);
     delete[] center_i;
     delete[] radii_i;
     MaskCoordinates *p_result = new MaskCoordinates(rank, dims_i);
@@ -462,15 +464,14 @@ Mask* OutputDataFunctions::CreateEllipticMask(
 }
 
 
-Mask* OutputDataFunctions::CreateEllipticMask(
-    const OutputData<double>& data, double xc, double yc, double rx, double ry)
+Mask* OutputDataFunctions::CreateEllipticMask(const OutputData<double>& data, double xc, double yc, double rx, double ry, bool invert_flag)
 {
     if(data.getRank() != 2) throw LogicErrorException(
             "OutputDataFunctions::CreateRectangularMask2D() -> Error! Number"
             " of dimensions should be 2");
     const double center[2]={xc, yc};
     const double radii[2]={rx, ry};
-    return OutputDataFunctions::CreateEllipticMask(data, center, radii);
+    return OutputDataFunctions::CreateEllipticMask(data, center, radii, invert_flag);
 }
 
 //double OutputDataFunctions::GetDifference(const OutputData<double> &result,
diff --git a/Core/Tools/src/VariableBinAxis.cpp b/Core/Tools/src/VariableBinAxis.cpp
index 9f8b3023e009d74d709e65399b0cf4cedee17429..155b7228f6584758a9e590c5276642a88b5ac66a 100644
--- a/Core/Tools/src/VariableBinAxis.cpp
+++ b/Core/Tools/src/VariableBinAxis.cpp
@@ -70,12 +70,18 @@ size_t VariableBinAxis::findClosestIndex(double value) const
     if(m_bin_boundaries.size()<2) {
         throw ClassInitializationException("VariableBinAxis::findClosestIndex() -> Error! VariableBinAxis not  correctly initialized");
     }
-    if (value < getMin() || value >= getMax()) {
-        std::ostringstream ostr;
-        ostr << "VariableBinAxis::findClosestIndex() -> Error! Given value not in any bin. ";
-        ostr << "value:" << value << " name:" << getName() << " min:" << getMin() << " max:" << getMax();
-        throw OutOfBoundsException(ostr.str());
+//    if (value < getMin() || value >= getMax()) {
+//        std::ostringstream ostr;
+//        ostr << "VariableBinAxis::findClosestIndex() -> Error! Given value not in any bin. ";
+//        ostr << "value:" << value << " name:" << getName() << " min:" << getMin() << " max:" << getMax();
+//        throw OutOfBoundsException(ostr.str());
+//    }
+    if( value < getMin()) {
+        return 0;
+    } else if(value >= getMax()) {
+        return m_nbins-1;
     }
+
     std::vector<double>::const_iterator top_limit = std::lower_bound(m_bin_boundaries.begin(), m_bin_boundaries.end(), value);
     if( *top_limit != value ) --top_limit;
     size_t nbin = top_limit - m_bin_boundaries.begin();
diff --git a/Examples/python/fitting/ex002_FitCylindersAndPrisms/FitCylindersPrisms_detailed.py b/Examples/python/fitting/ex002_FitCylindersAndPrisms/FitCylindersPrisms_detailed.py
index c1caa63ef9ff2f87763ab5035c7e51b6d0a208d7..63e33f6c24e4cabceada807d959a25a05e879546 100644
--- a/Examples/python/fitting/ex002_FitCylindersAndPrisms/FitCylindersPrisms_detailed.py
+++ b/Examples/python/fitting/ex002_FitCylindersAndPrisms/FitCylindersPrisms_detailed.py
@@ -145,6 +145,7 @@ class DrawObserver(IObserver):
                 pylab.text(0.01, 0.55 - i*0.1, str(fitpars[i].getName()) + " " + str(fitpars[i].getValue())[0:5] )
 
             pylab.draw()
+            pylab.pause(0.01)
 
 
 
diff --git a/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice.py b/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice.py
index adc8db22f6ff51893702b7ab88b2047a40ffa209..0845a4e457da35252c6c7b248c4e32e2fd134f0f 100644
--- a/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice.py
+++ b/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice.py
@@ -121,6 +121,7 @@ class DrawObserver(IObserver):
                 pylab.text(0.01, 0.55 - i*0.1, str(fitpars[i].getName()) + " " + str(fitpars[i].getValue())[0:5] )
 
             pylab.draw()
+            pylab.pause(0.01)
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py b/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py
index d521e7d840f81016b2b64132d37754d35bfc7f65..8a4e6d7a0b912fadc060d749f86d1992a5ffe228 100644
--- a/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py
+++ b/Examples/python/fitting/ex003_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py
@@ -141,6 +141,7 @@ class DrawObserver(IObserver):
                 pylab.text(0.01, 0.55 - i*0.1, str(fitpars[i].getName()) + " " + str(fitpars[i].getValue())[0:5] )
 
             pylab.draw()
+            pylab.pause(0.01)
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex004_FitScaleAndShift/FitScaleAndShift.py b/Examples/python/fitting/ex004_FitScaleAndShift/FitScaleAndShift.py
index 6691ee8537696e23d673a7df0c731c55e4faed6c..b6fbfab4b6f166c4a4f9fa8c9e4bdfd908c8ca9a 100644
--- a/Examples/python/fitting/ex004_FitScaleAndShift/FitScaleAndShift.py
+++ b/Examples/python/fitting/ex004_FitScaleAndShift/FitScaleAndShift.py
@@ -123,6 +123,7 @@ class DrawObserver(IObserver):
                 pylab.text(0.01, 0.55 - i*0.1, str(fitpars[i].getName()) + " " + str(fitpars[i].getValue())[0:5] )
 
             pylab.draw()
+            pylab.pause(0.01)
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex005_FitWithMasks/FitWithMasks.py b/Examples/python/fitting/ex005_FitWithMasks/FitWithMasks.py
index 37b0fc06c8c2ac1fdef561645139ecbd25856a93..d3518e2f882649b74788e9b71d00192d95cbceba 100644
--- a/Examples/python/fitting/ex005_FitWithMasks/FitWithMasks.py
+++ b/Examples/python/fitting/ex005_FitWithMasks/FitWithMasks.py
@@ -102,7 +102,8 @@ class DrawObserver(IObserver):
             pylab.colorbar(im)
             pylab.title('Simulated data')
             # plotting difference map
-            diff_map = (real_data - simulated_data)/(real_data + 1)
+            #diff_map = (real_data - simulated_data)/(real_data + 1)
+            diff_map= fit_suite.getFitObjects().getChiSquaredMap().getArray()
             pylab.subplot(2, 2, 3)
             im = pylab.imshow(numpy.rot90(diff_map, 1), norm=matplotlib.colors.LogNorm(), extent=[-1.0, 1.0, 0, 2.0], vmin = 0.001, vmax = 1.0)
             pylab.colorbar(im)
@@ -118,6 +119,7 @@ class DrawObserver(IObserver):
                 pylab.text(0.01, 0.55 - i*0.1, str(fitpars[i].getName()) + " " + str(fitpars[i].getValue())[0:5] )
 
             pylab.draw()
+            pylab.pause(0.01)
 
 
 def run_fitting():
@@ -132,7 +134,15 @@ def run_fitting():
 
     fit_suite = FitSuite()
 
-    IntensityDataFunctions.setRectangularMask(real_data, -0.1*degree, 0.0*degree, 0.1*degree, 2.0*degree) # x1,y1,x2,y2
+    # masking example: two excluded rectangles on the plot
+    IntensityDataFunctions.addRectangularMask(real_data, -0.1*degree, 0.1*degree, 0.1*degree, 0.2*degree)  # x1,y1,x2,y2
+    IntensityDataFunctions.addRectangularMask(real_data, -0.1*degree, 1.0*degree, 0.1*degree, 1.2*degree)  # x1,y1,x2,y2
+
+    # another mask example: one big square with two excluded areas on it
+    # IntensityDataFunctions.addRectangularMask(real_data, -0.6*degree, 0.0*degree, 0.6*degree, 1.5*degree, True)
+    # IntensityDataFunctions.addRectangularMask(real_data, -0.1*degree, 0.1*degree, 0.1*degree, 0.2*degree)
+    # IntensityDataFunctions.addEllipticMask(real_data, 0.0*degree, 1.2*degree, 0.3*degree, 0.2*degree)
+
     fit_suite.addSimulationAndRealData(simulation, real_data)
 
     fit_suite.initPrint(10)
@@ -141,8 +151,8 @@ def run_fitting():
     fit_suite.attachObserver(draw_observer)
 
     # setting fitting parameters with starting values
-    fit_suite.addFitParameter("*/FormFactorCylinder/height", 8.*nanometer, 0.01*nanometer, AttLimits.limited(4., 12.))
-    fit_suite.addFitParameter("*/FormFactorCylinder/radius", 8.*nanometer, 0.01*nanometer, AttLimits.limited(4., 12.))
+    fit_suite.addFitParameter("*/FormFactorCylinder/radius", 6.*nanometer, 0.01*nanometer, AttLimits.limited(4., 8.))
+    fit_suite.addFitParameter("*/FormFactorCylinder/height", 9.*nanometer, 0.01*nanometer, AttLimits.limited(8., 12.))
 
     # running fit
     fit_suite.runFit()
diff --git a/Fit/FitKernel/inc/FitSuiteObjects.h b/Fit/FitKernel/inc/FitSuiteObjects.h
index 7a6ce296109946b4bf0acaaa8538a8b3fc97ab78..4265455dc0ce01adda8bd477d8174f6f1618cf0a 100644
--- a/Fit/FitKernel/inc/FitSuiteObjects.h
+++ b/Fit/FitKernel/inc/FitSuiteObjects.h
@@ -101,6 +101,10 @@ class BA_CORE_API_  FitSuiteObjects : public IParameterized
     void setNfreeParameters(int nfree_parameters ) {
         m_nfree_parameters = nfree_parameters; }
 
+    //! Returns chi-squared map
+    OutputData<double> * getChiSquaredMap(size_t i_item = 0);
+
+
  protected:
     //! Registers some class members for later access via parameter pool
     virtual void init_parameters() {}
diff --git a/Fit/FitKernel/src/FitSuiteObjects.cpp b/Fit/FitKernel/src/FitSuiteObjects.cpp
index 7202caee59d6cf94b70e603e08a576b184c10161..8aac08d1d56451d122bb0a12fb1b9e7104f73138 100644
--- a/Fit/FitKernel/src/FitSuiteObjects.cpp
+++ b/Fit/FitKernel/src/FitSuiteObjects.cpp
@@ -152,3 +152,8 @@ std::string FitSuiteObjects::addParametersToExternalPool(
 
     return new_path;
 }
+
+OutputData<double> *FitSuiteObjects::getChiSquaredMap(size_t i_item)
+{
+    return m_fit_objects[check_index(i_item)]->getChiSquaredModule()->createChi2DifferenceMap();
+}
diff --git a/Fit/PythonAPI/src/FitSuiteObjects.pypp.cpp b/Fit/PythonAPI/src/FitSuiteObjects.pypp.cpp
index 0fbe73787eefc62bbbdf40ad3ef3aa0244135228..d664dbd5e7f9604718baa70b52cd30d95e253131 100644
--- a/Fit/PythonAPI/src/FitSuiteObjects.pypp.cpp
+++ b/Fit/PythonAPI/src/FitSuiteObjects.pypp.cpp
@@ -130,6 +130,17 @@ void register_FitSuiteObjects_class(){
                 "clear"
                 , clear_function_type( &::FitSuiteObjects::clear ) );
         
+        }
+        { //::FitSuiteObjects::getChiSquaredMap
+        
+            typedef ::OutputData< double > * ( ::FitSuiteObjects::*getChiSquaredMap_function_type)( ::std::size_t ) ;
+            
+            FitSuiteObjects_exposer.def( 
+                "getChiSquaredMap"
+                , getChiSquaredMap_function_type( &::FitSuiteObjects::getChiSquaredMap )
+                , ( bp::arg("i_item")=(::std::size_t)(0) )
+                , bp::return_value_policy< bp::manage_new_object >() );
+        
         }
         { //::FitSuiteObjects::getChiSquaredValue
         
diff --git a/Fit/PythonAPI/src/PythonModule.cpp b/Fit/PythonAPI/src/PythonModule.cpp
index ef8312b19c0e5f34a1dd036a15ac8689fdb68a4c..9e316599b4cf7992d66a9756e63a29a6c5020ac4 100644
--- a/Fit/PythonAPI/src/PythonModule.cpp
+++ b/Fit/PythonAPI/src/PythonModule.cpp
@@ -5,38 +5,38 @@ GCC_DIAG_OFF(missing-field-initializers)
 #include "boost/python.hpp"
 GCC_DIAG_ON(unused-parameter)
 GCC_DIAG_ON(missing-field-initializers)
-#include "IntensityFunctionSqrt.pypp.h"
-#include "MinimizerFactory.pypp.h"
-#include "IMinimizer.pypp.h"
-#include "vector_string_t.pypp.h"
-#include "SquaredFunctionSystematicError.pypp.h"
-#include "IntensityNormalizer.pypp.h"
-#include "IIntensityFunction.pypp.h"
-#include "INamed.pypp.h"
+#include "FitObject.pypp.h"
 #include "IntensityFunctionLog.pypp.h"
-#include "FitSuiteParameters.pypp.h"
-#include "AttFitting.pypp.h"
-#include "FitParameter.pypp.h"
-#include "IntensityScaleAndShiftNormalizer.pypp.h"
+#include "IntensityFunctionSqrt.pypp.h"
+#include "FitStrategyDefault.pypp.h"
 #include "IChiSquaredModule.pypp.h"
+#include "AttFitting.pypp.h"
+#include "FitStrategyAdjustParameters.pypp.h"
 #include "FitStrategyAdjustMinimizer.pypp.h"
-#include "IFitStrategy.pypp.h"
-#include "FitStrategyFixParameters.pypp.h"
+#include "ISquaredFunction.pypp.h"
 #include "SquaredFunctionGaussianError.pypp.h"
+#include "IIntensityFunction.pypp.h"
+#include "IntensityScaleAndShiftNormalizer.pypp.h"
+#include "SquaredFunctionMeanSquaredError.pypp.h"
+#include "FitSuiteParameters.pypp.h"
+#include "IMinimizer.pypp.h"
 #include "IIntensityNormalizer.pypp.h"
+#include "FitParameter.pypp.h"
+#include "INamed.pypp.h"
+#include "FitStrategyReleaseParameters.pypp.h"
 #include "FitSuite.pypp.h"
-#include "FitStrategyAdjustParameters.pypp.h"
+#include "IntensityNormalizer.pypp.h"
+#include "FitStrategyFixParameters.pypp.h"
+#include "FitSuiteObjects.pypp.h"
+#include "SquaredFunctionSystematicError.pypp.h"
+#include "MinimizerFactory.pypp.h"
 #include "ChiSquaredModule.pypp.h"
+#include "SquaredFunctionSimError.pypp.h"
 #include "MinimizerOptions.pypp.h"
-#include "SquaredFunctionDefault.pypp.h"
-#include "SquaredFunctionMeanSquaredError.pypp.h"
-#include "ISquaredFunction.pypp.h"
-#include "FitStrategyDefault.pypp.h"
 #include "AttLimits.pypp.h"
-#include "FitObject.pypp.h"
-#include "FitSuiteObjects.pypp.h"
-#include "SquaredFunctionSimError.pypp.h"
-#include "FitStrategyReleaseParameters.pypp.h"
+#include "SquaredFunctionDefault.pypp.h"
+#include "vector_string_t.pypp.h"
+#include "IFitStrategy.pypp.h"
 
 BOOST_PYTHON_MODULE(libBornAgainFit){
     boost::python::docstring_options doc_options(true, true, false);
diff --git a/Tests/UnitTests/TestCore/DWBASimulationTest.h b/Tests/UnitTests/TestCore/DWBASimulationTest.h
index d233aa7f388fc332abcc286187617d5d77c15a04..94042c88b677e091bc3159dd2a54e0d91c2283f2 100644
--- a/Tests/UnitTests/TestCore/DWBASimulationTest.h
+++ b/Tests/UnitTests/TestCore/DWBASimulationTest.h
@@ -83,7 +83,7 @@ TEST_F(DWBASimulationTest, MaskedThreadIterator)
         (*it) = double(index++);
     }
 
-    Mask *mask1 = OutputDataFunctions::CreateRectangularMask(m_data, 1.99, 0.99, 7.01, 3.01);
+    Mask *mask1 = OutputDataFunctions::CreateRectangularMask(m_data, 1.99, 0.99, 7.01, 3.01, true);
     m_data.setMask(*mask1);
 
     m_sim.setDetectorParameters(m_data);
diff --git a/Tests/UnitTests/TestCore/FixedBinAxisTest.h b/Tests/UnitTests/TestCore/FixedBinAxisTest.h
index bcd578b74ef09a89894f62fd0846a65215cb83e5..0464a4f20ddf62c8c5d972b7c43c2d53331b4754 100644
--- a/Tests/UnitTests/TestCore/FixedBinAxisTest.h
+++ b/Tests/UnitTests/TestCore/FixedBinAxisTest.h
@@ -49,7 +49,8 @@ TEST_F(FixedBinAxisTest, FindClosestIndex)
     EXPECT_EQ( size_t(0), v1.findClosestIndex(0.25));
     EXPECT_EQ( size_t(1), v1.findClosestIndex(0.5));
     EXPECT_EQ( size_t(1), v1.findClosestIndex(0.6));
-    ASSERT_THROW( v1.findClosestIndex(1.0), Exceptions::OutOfBoundsException);
+//    ASSERT_THROW( v1.findClosestIndex(1.0), Exceptions::OutOfBoundsException);
+    EXPECT_EQ( size_t(1), v1.findClosestIndex(1.0));
 
     FixedBinAxis v2("name", 3, -1.5, 1.5);
     EXPECT_EQ(size_t(0), v2.findClosestIndex(-1.5));
@@ -58,7 +59,8 @@ TEST_F(FixedBinAxisTest, FindClosestIndex)
     EXPECT_EQ(size_t(1), v2.findClosestIndex(0.0));
     EXPECT_EQ(size_t(2), v2.findClosestIndex(0.5));
     EXPECT_EQ(size_t(2), v2.findClosestIndex(1.499));
-    ASSERT_THROW( v2.findClosestIndex(1.5), Exceptions::OutOfBoundsException);
+//    ASSERT_THROW( v2.findClosestIndex(1.5), Exceptions::OutOfBoundsException);
+    EXPECT_EQ(size_t(2), v2.findClosestIndex(1.5));
 }
 
 TEST_F(FixedBinAxisTest, CheckBin)
diff --git a/Tests/UnitTests/TestCore/IntensityDataFunctionsTest.h b/Tests/UnitTests/TestCore/IntensityDataFunctionsTest.h
index 712fa0be668b631668a16b829047ca3d175985f3..26315a8935396f7a232d9cdcdd35d1ce1e0b3fbb 100644
--- a/Tests/UnitTests/TestCore/IntensityDataFunctionsTest.h
+++ b/Tests/UnitTests/TestCore/IntensityDataFunctionsTest.h
@@ -67,4 +67,83 @@ TEST_F(IntensityDataFunctionsTest, ClipDataSetVariable)
 
 }
 
+
+//-------------------------------------------------  2.5
+// 2 ||  2  |  5  |  8  | 11  | 14 | 17 | 20 | 23 |
+//------------------------------------------------
+// 1 ||  1  |  4  |  7  | 10  | 13 | 16 | 19 | 22 |
+//------------------------------------------------
+// 0 ||  0  |  3  |  6  |  9  | 12 | 15 | 18 | 21 |
+//================================================== -0.5
+//   ||  0  |  1  |  2  |  3  |  4 |  5 |  6 |  7 |
+// -4.5                                         3.5
+
+TEST_F(IntensityDataFunctionsTest, AddRectangularMask)
+{
+    OutputData<double > data;
+    data.addAxis("x", 8, -4.5, 3.5);
+    data.addAxis("y", 3, -0.5, 2.5);
+    data.setAllTo(0.0);
+    IntensityDataFunctions::addRectangularMask(data, -3.0, -0.5, 1.0, 1.49);
+    IntensityDataFunctions::addRectangularMask(data, 1.5, 0.5, 3.5, 2.5);
+
+    for(size_t i=0; i<data.getAllocatedSize(); ++i) {
+        data[i] = i;
+    }
+
+    int index(0);
+
+    std::vector<double> xref = boost::assign::list_of(-4.0)(-4.0)(-4.0)(-3.0)(-2.0)(-1.0)(0.0)(1.0)(2.0)(3.0);
+    std::vector<double> yref = boost::assign::list_of(0.0)(1.0)(2.0)(2.0)(2.0)(2.0)(2.0)(2.0)(0.0)(0.0);
+    std::vector<double> vref = boost::assign::list_of(0)(1)(2)(5)(8)(11)(14)(17)(18)(21);
+    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+        double x = data.getValueOfAxis("x", it.getIndex());
+        double y = data.getValueOfAxis("y", it.getIndex());
+        EXPECT_EQ(x, xref[index]);
+        EXPECT_EQ(y, yref[index]);
+        EXPECT_EQ(*it, vref[index]);
+        ++index;
+    }
+    data.removeAllMasks();
+    index=0;
+    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+        EXPECT_EQ( int(index++), int(it.getIndex()) );
+    }
+}
+
+
+//TEST_F(IntensityDataFunctionsTest, AddRectangularMask2)
+//{
+//    OutputData<double > data;
+//    data.addAxis("x", 8, -4.5, 3.5);
+//    data.addAxis("y", 3, -0.5, 2.5);
+//    data.setAllTo(0.0);
+//    IntensityDataFunctions::addRectangularMask(data, -3.0, -0.5, 1.0, 1.49, false);
+//    IntensityDataFunctions::addRectangularMask(data, 1.5, 0.5, 3.5, 2.5, false);
+
+//    for(size_t i=0; i<data.getAllocatedSize(); ++i) {
+//        data[i] = i;
+//    }
+
+//    int index(0);
+
+//    std::vector<double> xref = boost::assign::list_of(-4.0)(-4.0)(-4.0)(-3.0)(-2.0)(-1.0)(0.0)(1.0)(2.0)(3.0);
+//    std::vector<double> yref = boost::assign::list_of(0.0)(1.0)(2.0)(2.0)(2.0)(2.0)(2.0)(2.0)(0.0)(0.0);
+//    std::vector<double> vref = boost::assign::list_of(0)(1)(2)(5)(8)(11)(14)(17)(18)(21);
+//    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+//        double x = data.getValueOfAxis("x", it.getIndex());
+//        double y = data.getValueOfAxis("y", it.getIndex());
+//        EXPECT_EQ(x, xref[index]);
+//        EXPECT_EQ(y, yref[index]);
+//        EXPECT_EQ(*it, vref[index]);
+//        ++index;
+//    }
+//    data.removeAllMasks();
+//    index=0;
+//    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+//        EXPECT_EQ( int(index++), int(it.getIndex()) );
+//    }
+//}
+
+
 #endif
diff --git a/Tests/UnitTests/TestCore/OutputDataTest.h b/Tests/UnitTests/TestCore/OutputDataTest.h
index 54c7ee27e49c7f63fcc033c131952427ca9f6bcd..ef1f7f713f19a83b14e95f799ab8855fca317807 100644
--- a/Tests/UnitTests/TestCore/OutputDataTest.h
+++ b/Tests/UnitTests/TestCore/OutputDataTest.h
@@ -214,22 +214,22 @@ TEST_F(OutputDataTest, ValueOfAxis)
 
 // y |
 // --------------------------------------------
-// 1 | 5   11  17  23  29  35  41  47  53  59 |
-// 1 | 4   10  16  22  28  34  40  46  52  58 |
-// 1 | 3   9   15  21  27  33  39  45  51  57 |
-// 1 | 2   8   14  20  26  32  38  44  50  56 |
+// 5 | 5   11  17  23  29  35  41  47  53  59 |
+// 4 | 4   10  16  22  28  34  40  46  52  58 |
+// 3 | 3   9   15  21  27  33  39  45  51  57 |
+// 2 | 2   8   14  20  26  32  38  44  50  56 |
 // 1 | 1   7   13  19  25  31  37  43  49  55 |
 // 0 | 0   6   12  18  24  30  36  42  48  54 |
 // --------------------------------------------
 //   | 0   1   2   3   4   5   6   7   8   9  | x
 
-TEST_F(OutputDataTest, SetRectangularMask)
+TEST_F(OutputDataTest, SetInverseRectangularMask)
 {
     OutputData<double > data;
     data.addAxis("x", 10, 0., 10.);
     data.addAxis("y", 6, 0., 6.);
     data.setAllTo(0.0);
-    IntensityDataFunctions::setRectangularMask(data, 1.0, 1.0, 4.99, 2.99);
+    IntensityDataFunctions::setRectangularMask(data, 1.0, 1.0, 4.99, 2.99, true);
 
     for(size_t i=0; i<data.getAllocatedSize(); ++i) {
         data[i] = i;
@@ -256,6 +256,86 @@ TEST_F(OutputDataTest, SetRectangularMask)
 }
 
 
+//-------------------------------------------------  2.5
+// 2 ||  2  |  5  |  8  | 11  | 14 | 17 | 20 | 23 |
+//------------------------------------------------
+// 1 ||  1  |  4  |  7  | 10  | 13 | 16 | 19 | 22 |
+//------------------------------------------------
+// 0 ||  0  |  3  |  6  |  9  | 12 | 15 | 18 | 21 |
+//================================================== -0.5
+//   ||  0  |  1  |  2  |  3  |  4 |  5 |  6 |  7 |
+// -4.5                                         3.5
+
+TEST_F(OutputDataTest, SetRectangularMask)
+{
+    OutputData<double > data;
+    data.addAxis("x", 8, -4.5, 3.5);
+    data.addAxis("y", 3, -0.5, 2.5);
+    data.setAllTo(0.0);
+    IntensityDataFunctions::setRectangularMask(data, -3.0, -0.5, 1.0, 1.49);
+
+    for(size_t i=0; i<data.getAllocatedSize(); ++i) {
+        data[i] = i;
+    }
+
+    int index(0);
+
+    std::vector<double> xref = boost::assign::list_of(-4.0)(-4.0)(-4.0)(-3.0)(-2.0)(-1.0)(0.0)(1.0)(2.0)(2.0)(2.0)(3.0)(3.0)(3.0);
+    std::vector<double> yref = boost::assign::list_of(0.0)(1.0)(2.0)(2.0)(2.0)(2.0)(2.0)(2.0)(0.0)(1.0)(2.0)(0.0)(1.0)(2.0);
+    std::vector<double> vref = boost::assign::list_of(0)(1)(2)(5)(8)(11)(14)(17)(18)(19)(20)(21)(22)(23);
+    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+        double x = data.getValueOfAxis("x", it.getIndex());
+        double y = data.getValueOfAxis("y", it.getIndex());
+        EXPECT_EQ(x, xref[index]);
+        EXPECT_EQ(y, yref[index]);
+        EXPECT_EQ(*it, vref[index]);
+        ++index;
+    }
+    data.removeAllMasks();
+    index=0;
+    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+        EXPECT_EQ( int(index++), int(it.getIndex()) );
+    }
+}
+
+
+//-------------------------------------------------  2.5
+// 2 ||  2  |  5  |  8  | 11  | 14 | 17 | 20 | 23 |
+//------------------------------------------------
+// 1 ||  1  |  4  |  7  | 10  | 13 | 16 | 19 | 22 |
+//------------------------------------------------
+// 0 ||  0  |  3  |  6  |  9  | 12 | 15 | 18 | 21 |
+//================================================== -0.5
+//   ||  0  |  1  |  2  |  3  |  4 |  5 |  6 |  7 |
+// -4.5                                         3.5
+
+TEST_F(OutputDataTest, SetEllipticMask)
+{
+    OutputData<double > data;
+    data.addAxis("x", 8, -4.5, 3.5);
+    data.addAxis("y", 3, -0.5, 2.5);
+    data.setAllTo(0.0);
+    IntensityDataFunctions::setEllipticMask(data, 0.0, 1.0, 2.0, 1.0);
+
+    for(size_t i=0; i<data.getAllocatedSize(); ++i) {
+        data[i] = i;
+    }
+
+    int index(0);
+
+    std::vector<double> vref = boost::assign::list_of(0)(1)(2)(3)(4)(5)(6)(8)(9)(11)(15)(17)(18)(20)(21)(22)(23);
+    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+        EXPECT_EQ(*it, vref[index]);
+        ++index;
+    }
+    data.removeAllMasks();
+    index=0;
+    for(OutputData<double>::iterator it = data.begin(); it!=data.end(); ++it) {
+        EXPECT_EQ( int(index++), int(it.getIndex()) );
+    }
+}
+
+
 TEST_F(OutputDataTest, RectangularMaskVariableAxis)
 {
     static const double x_arr[] = {-5., -3., -2., 0.0, 0.5, 1.0, 3.0, 5.0, 6.0};
@@ -273,7 +353,7 @@ TEST_F(OutputDataTest, RectangularMaskVariableAxis)
         data[i] = i;
     }
 
-    IntensityDataFunctions::setRectangularMask(data, -2.5, 1.5, 0.99, 4.99);
+    IntensityDataFunctions::setRectangularMask(data, -2.5, 1.5, 0.99, 4.99, true);
     int index(0);
 
     std::vector<double> xref = boost::assign::list_of(-2.5)(-2.5)(-1.0)(-1.0)(0.25)(0.25)(0.75)(0.75);
@@ -388,7 +468,7 @@ TEST_F(OutputDataTest, ThreadInfoMaskedIterator)
         (*it) = double(index++);
     }
 
-    IntensityDataFunctions::setRectangularMask(data, 1.0, 1.0, 8.99, 2.99);
+    IntensityDataFunctions::setRectangularMask(data, 1.0, 1.0, 8.99, 2.99, true);
 
 
     const int nthreads = 4;
diff --git a/Tests/UnitTests/TestCore/VariableBinAxisTest.h b/Tests/UnitTests/TestCore/VariableBinAxisTest.h
index 4b8501bc46ea1c06f2a0fcc8a52f00a43afe5495..1a1e4b94ea10e7a7d83925fe947d091e217cb5f7 100644
--- a/Tests/UnitTests/TestCore/VariableBinAxisTest.h
+++ b/Tests/UnitTests/TestCore/VariableBinAxisTest.h
@@ -87,7 +87,8 @@ TEST_F(VariableBinAxisTest, FindClosestIndex)
     EXPECT_EQ( size_t(0), v1.findClosestIndex(0.25));
     EXPECT_EQ( size_t(1), v1.findClosestIndex(0.5));
     EXPECT_EQ( size_t(1), v1.findClosestIndex(0.6));
-    ASSERT_THROW( v1.findClosestIndex(1.0), Exceptions::OutOfBoundsException);
+//    ASSERT_THROW( v1.findClosestIndex(1.0), Exceptions::OutOfBoundsException);
+    EXPECT_EQ( size_t(1), v1.findClosestIndex(1.0));
 
     static const double arr2[] = {-1.5, -0.5, 0.5, 1.5};
     std::vector<double> values2 (arr2, arr2 + sizeof(arr2) / sizeof(arr2[0]) );
@@ -98,7 +99,8 @@ TEST_F(VariableBinAxisTest, FindClosestIndex)
     EXPECT_EQ(size_t(1), v2.findClosestIndex(0.0));
     EXPECT_EQ(size_t(2), v2.findClosestIndex(0.5));
     EXPECT_EQ(size_t(2), v2.findClosestIndex(1.499));
-    ASSERT_THROW( v2.findClosestIndex(1.5), Exceptions::OutOfBoundsException);
+//    ASSERT_THROW( v2.findClosestIndex(1.5), Exceptions::OutOfBoundsException);
+    EXPECT_EQ(size_t(2), v2.findClosestIndex(1.5));
 
     static const double arr3[] = {-1.0, -0.5, 0.5, 1.0, 2.0};
     std::vector<double> values3 (arr3, arr3 + sizeof(arr3) / sizeof(arr2[0]) );
diff --git a/dev-tools/python-bindings/settings_fit.py b/dev-tools/python-bindings/settings_fit.py
index 244175ae53bc2dac0e55a538fc9988611333623d..f38815a691674bd9fada2acb0319c41c4932ba41 100644
--- a/dev-tools/python-bindings/settings_fit.py
+++ b/dev-tools/python-bindings/settings_fit.py
@@ -172,6 +172,11 @@ def ManualClassTunings(mb):
         if "__gnu_cxx::__normal_iterator" in fun.decl_string:
             fun.exclude()
 
+    cl = mb.class_("FitSuiteObjects")
+    cl.member_function("getChiSquaredMap").call_policies = \
+        call_policies.return_value_policy(call_policies.manage_new_object)
+
+
 
 # excluding specific member functions
 def ManualExcludeMemberFunctions(mb):