diff --git a/App/inc/ROOTMinimizer.h b/App/inc/ROOTMinimizer.h
index 96d5d1793f1626a0bb5b4a76b5e292072405d12b..34ffc8809e9389708371193ff383897b7e507d94 100644
--- a/App/inc/ROOTMinimizer.h
+++ b/App/inc/ROOTMinimizer.h
@@ -17,6 +17,7 @@
 #include "IMinimizer.h"
 #include "OutputData.h"
 #include "Exceptions.h"
+#include "ROOTMinimizerFunction.h"
 #include <string>
 // from ROOT
 #include "Math/Minimizer.h"
@@ -31,12 +32,14 @@
 class ROOTMinimizer : public IMinimizer
 {
 public:
-    ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type);
+    ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type=std::string());
     virtual ~ROOTMinimizer();
 
     virtual void setVariable(int index, const FitParameter *par) ;
-    virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1);
-//    virtual void setFunctionAndGradient(boost::function<double(const double *)> fcn, boost::function<double(const double *, int)> fcn_deriv, int ndim=1);
+
+    virtual void setFunction(function_t fcn, int ndims, element_function_t element_fcn, int nelements);
+//    virtual void setElementFunction(element_function_t element_fcn, int ndims, int nelements) { m_element_fcn = element_fcn, m_ndims = ndims; m_nelements = nelements; }
+
     virtual void minimize();
 
     //! return pointer to created minimizer
@@ -71,8 +74,16 @@ public:
 
 private:
     ROOT::Math::Minimizer *m_root_minimizer;
-    ROOT::Math::Functor *m_fcn; //! function to minimize
+    ROOTMinimizerFunction * m_minfunc;
+    ROOTMinimizerElementFunction * m_minfunc_element;
+
+//    IMinimizerFunction *m_minimizer_function;
+//    ROOT::Math::Functor *m_fcn; //! function to minimize
 //    ROOT::Math::GradFunctor *m_fcn_grad; //! gradient of function to minimize
+//    function_t m_fcn;
+//    element_function_t m_element_fcn;
+//    int m_ndims;
+//    int m_nelements;
 };
 
 #endif // ROOTMINIMIZER_H
diff --git a/App/inc/ROOTMinimizerFunction.h b/App/inc/ROOTMinimizerFunction.h
index 09489e3095f4dce01fbffe1e0cd0f864ed6d7ad9..357cdf3adb11c1377341c0cea7fc4af72e6ef663 100644
--- a/App/inc/ROOTMinimizerFunction.h
+++ b/App/inc/ROOTMinimizerFunction.h
@@ -15,64 +15,58 @@
 //! @date   Dec 10, 2012
 
 
+#include "IMinimizer.h"
+#include "Math/Functor.h"
+#include "Math/FitMethodFunction.h"
+
+//- -------------------------------------------------------------------
+//! @class ROOTMinimizerFunction
+//! @brief Basic minimizer function
+//- -------------------------------------------------------------------
+class ROOTMinimizerFunction : public ROOT::Math::Functor
+{
+public:
+    ROOTMinimizerFunction(IMinimizer::function_t fcn, int ndims ) : ROOT::Math::Functor(fcn, ndims), m_fcn(fcn) {}
+    virtual ~ROOTMinimizerFunction(){}
+    IMinimizer::function_t m_fcn;
+};
+
+
+//- -------------------------------------------------------------------
+//! @class ROOTMinimizerElementFunction
+//! @brief Minimizer function with access to single data element residuals.
+//! Required by Fumili, Fumili2 and GSLMultiMin minimizers
+//- -------------------------------------------------------------------
+class ROOTMinimizerElementFunction : public ROOT::Math::FitMethodFunction
+{
+public:
+    typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t  Type_t;
+
+    ROOTMinimizerElementFunction(IMinimizer::function_t fcn, int ndims, IMinimizer::element_function_t element_fcn, int nelements)
+        : ROOT::Math::FitMethodFunction(ndims, nelements)
+        , m_fcn(fcn)
+        , m_element_fcn(element_fcn)
+        , m_ndims(ndims)
+        , m_nelements(nelements) { }
+
+    virtual ~ROOTMinimizerElementFunction(){}
+
+    Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; }
+    ROOT::Math::IMultiGenFunction * Clone() const { return new ROOTMinimizerElementFunction(m_fcn, m_ndims, m_element_fcn, m_nelements); }
+
+    //! evaluation of chi2
+    double DoEval(const double * par) const { return m_fcn(par); }
+
+    //! evaluation of single data element residual
+    double DataElement(const double *par, unsigned int i, double *g = 0) const { return m_element_fcn(par,i,g); }
+
+private:
+    IMinimizer::function_t m_fcn;
+    IMinimizer::element_function_t m_element_fcn;
+    int m_ndims;
+    int m_nelements;
+};
 
-//#include "Math/Minimizer.h"
-//#include "Math/Factory.h"
-//#include "Math/FitMethodFunction.h"
-
-
-////- -------------------------------------------------------------------
-////! @class ROOTMinimizerFunction
-////! @brief
-////- -------------------------------------------------------------------
-
-//class ROOTMinimizerFunction : public ROOT::Math::FitMethodFunction
-//{
-//public:
-//    typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t  Type_t;
-//    typedef boost::function<double(const double *)> function_eval_t;
-//    typedef boost::function<double(const double *, unsigned int, double *)> function_element_t;
-
-//    ROOTMinimizerFunction();
-//    virtual ~ROOTMinimizerFunction(){}
-
-//    Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; }
-//    ROOT::Math::IMultiGenFunction * Clone() const { return new MyChi2Function(m_test); }
-
-//    double DoEval(const double * par) const {
-
-
-//};
-
-
-//class MyChi2Function : public ROOT::Math::FitMethodFunction
-//{
-//public:
-//    typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t  Type_t;
-
-//    MyChi2Function(TestFumiliLMA *test) : ROOT::Math::FitMethodFunction(test->m_ndim, test->m_real_data->getAllocatedSize()), m_test(test) {}
-//    virtual ~MyChi2Function(){}
-
-//    Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; }
-//    ROOT::Math::IMultiGenFunction * Clone() const { return new MyChi2Function(m_test); }
-
-//    // evaluation of the all chi2
-//    double DoEval(const double * par) const {
-//       int ndata = NPoints();
-//       std::cout << "DoEval: " << ndata << std::endl;
-//       double chi2 = 0;
-//       for (int i = 0; i <  ndata; ++i) {
-//          double res = DataElement( par, i);
-//          chi2 += res*res;
-//       }
-//       //std::cout << "DoEval: chi" << chi2/double(ndata) << std::endl;
-//       return chi2/double(ndata);
-//    }
-
-//    double DataElement(const double *par, unsigned int i, double *g = 0) const;
-
-//    TestFumiliLMA *m_test;
-//};
 
 
 
diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index 606db559ab9e2a0f90f25b83469130adb429ce4c..50f2192ec813c4acaaaee4abddb618d203887780 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -123,8 +123,8 @@ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const s
     OutputData<double>::const_iterator it = output.begin();
     while (it != output.end())
     {
-        double x = output.getValueOfAxis( haxises[0].name, it.getIndex() );
-        double y = output.getValueOfAxis( haxises[1].name, it.getIndex() );
+        double x = output.getValueOfAxis( haxises[0].name.c_str(), it.getIndex() );
+        double y = output.getValueOfAxis( haxises[1].name.c_str(), it.getIndex() );
         double value = *it++;
         hist2->Fill(x, y, value);
     }
@@ -216,7 +216,7 @@ TH1 *IsGISAXSTools::getOutputDataTH123D(const OutputData<double>& output, const
     {
         std::vector<double > xyz;
         for(size_t i_axis=0; i_axis<haxises.size(); ++i_axis) {
-            xyz.push_back(output.getValueOfAxis( haxises[i_axis].name, it.getIndex() ) );
+            xyz.push_back(output.getValueOfAxis( haxises[i_axis].name.c_str(), it.getIndex() ) );
         }
         double value = *it++;
         if(hist1) hist1->Fill(xyz[0], value);
diff --git a/App/src/ROOTMinimizer.cpp b/App/src/ROOTMinimizer.cpp
index 5ae2f91574e8a7f7df78cc91fd29d4aa4529f2bc..848a60876d1cf0e80e146ec10104eb901abec518 100644
--- a/App/src/ROOTMinimizer.cpp
+++ b/App/src/ROOTMinimizer.cpp
@@ -1,10 +1,13 @@
 #include "ROOTMinimizer.h"
 #include "Exceptions.h"
 #include "Utils.h"
+#include "ROOTMinimizerFunction.h"
 #include <iomanip>
 #include <sstream>
 
-ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type) : m_fcn(0)
+ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type)
+    : m_minfunc(0)
+    , m_minfunc_element(0)
 {
     if(minimizer_name == "Minuit2" && algo_type == "Fumili2") {
         throw LogicErrorException("ROOTMinimizer::ROOTMinimizer() -> Error! Use word Fumili instead of 'Fumili2'");
@@ -22,7 +25,8 @@ ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::strin
 ROOTMinimizer::~ROOTMinimizer()
 {
     delete m_root_minimizer;
-    delete m_fcn;
+    delete m_minfunc;
+    delete m_minfunc_element;
 }
 
 
@@ -61,19 +65,23 @@ void ROOTMinimizer::minimize()
 /* ************************************************************************* */
 // set fcn function for minimizer
 /* ************************************************************************* */
-void ROOTMinimizer::setFunction(boost::function<double(const double *)> fcn, int ndim)
+void ROOTMinimizer::setFunction(function_t fcn, int ndims, element_function_t element_fcn, int nelements)
 {
-    m_fcn = new ROOT::Math::Functor(fcn, ndim);
-    m_root_minimizer->SetFunction(*m_fcn);
-}
-
+    if( fcn && element_fcn ) {
+        std::cout << " ROOTMinimizer::setFunction() -> XXX 1.1 making ROOTMinimizerElementFunction " << std::endl;
+        delete m_minfunc;
+        m_minfunc_element = new ROOTMinimizerElementFunction(fcn, ndims, element_fcn, nelements);
+        m_root_minimizer->SetFunction(*m_minfunc_element);
+    } else if( fcn ) {
+        std::cout << " ROOTMinimizer::setFunction() -> XXX 1.2 making ROOTMinimizerFunction" << std::endl;
+        delete m_minfunc;
+        m_minfunc = new ROOTMinimizerFunction(fcn, ndims);
+        m_root_minimizer->SetFunction(*m_minfunc);
+    } else {
+        throw LogicErrorException("ROOTMinimizer::minimize() -> Error! Can't guess minimization function type");
+    }
 
-//void ROOTMinimizer::setFunctionAndGradient(boost::function<double(const double *)> fcn, boost::function<double(const double *, int)> fcn_deriv, int ndim)
-//{
-//    m_fcn_grad = new ROOT::Math::GradFunctor(fcn, fcn_deriv, ndim);
-//    m_root_minimizer->SetFunction(*m_fcn_grad);
-//    m_root_minimizer->SetFunction(*m_fcn_grad);
-//}
+}
 
 
 /* ************************************************************************* */
diff --git a/App/src/TestFumiliLMA.cpp b/App/src/TestFumiliLMA.cpp
index cf6fc7be325d57f67877b85c107d0e1a5ecff25c..97e6b24d1bd88e1413098fac2375c78622e327d1 100644
--- a/App/src/TestFumiliLMA.cpp
+++ b/App/src/TestFumiliLMA.cpp
@@ -184,7 +184,7 @@ double MyChi2Function::DataElement(const double *pars, unsigned int i, double *g
 {
     double xx[2];
     std::cout << " DataElement() -> " << g << " " << i;
-    for(size_t i=0; i<m_test->m_func->GetNpar(); ++i) std::cout << " " << pars[i];
+    for(int i=0; i<m_test->m_func->GetNpar(); ++i) std::cout << " " << pars[i];
     std::cout << std::endl;
     const AxisDouble *xaxis = m_test->m_real_data->getAxis(0);
     const AxisDouble *yaxis = m_test->m_real_data->getAxis(1);
diff --git a/App/src/TestToyExperiment.cpp b/App/src/TestToyExperiment.cpp
index f32e2afa090ed732bcee0386ca665fd626e8dba2..ba5f6276adbdc784609aa553986e9d4c1c6b0049 100644
--- a/App/src/TestToyExperiment.cpp
+++ b/App/src/TestToyExperiment.cpp
@@ -17,15 +17,18 @@ void ToyExperiment::runSimulation()
     m_func->SetParameters(&pars[0]);
     m_intensity_map.setAllTo(0.0);
     OutputData<double>::iterator it = m_intensity_map.begin();
+    const std::string s_phi_f("phi_f");
+    const std::string s_alpha_f("alpha_f");
     while( it!= m_intensity_map.end() ) {
-        double phi_f = m_intensity_map.getValueOfAxis("phi_f", it.getIndex());
-        double alpha_f = m_intensity_map.getValueOfAxis("alpha_f", it.getIndex());
+        double phi_f = m_intensity_map.getValueOfAxis(s_phi_f, it.getIndex());
+        double alpha_f = m_intensity_map.getValueOfAxis(s_alpha_f, it.getIndex());
         double value = m_func->Eval(phi_f, alpha_f);
         (*it) = value;
         ++it;
     }
 }
 
+
 void ToyExperiment::init_parameters()
 {
     getParameterPool()->clear();
@@ -72,7 +75,8 @@ void TestToyExperiment::execute()
 
     // setting up fitSuite
     FitSuite *m_fitSuite = new FitSuite();
-    m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
+    //m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
+    m_fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiFit") );
     m_fitSuite->attachObserver( new FitSuiteObserverPrint() );
     m_fitSuite->attachObserver( new FitSuiteObserverDraw() );
 
diff --git a/Core/Algorithms/inc/ChiSquaredModule.h b/Core/Algorithms/inc/ChiSquaredModule.h
index d7b1c1234c6ca731df4b81262de5096d746bee37..9acbb22f28eef9852a169b0a48980370a64cffdb 100644
--- a/Core/Algorithms/inc/ChiSquaredModule.h
+++ b/Core/Algorithms/inc/ChiSquaredModule.h
@@ -31,6 +31,9 @@ public:
     //! return output data which contains chi^2 values
     virtual OutputData<double > *createChi2DifferenceMap() const;
 
+    //! return residual between data and simulation for single element
+    virtual double getResidualValue(int index ) const;
+
 private:
     // disabling assignment operator
     ChiSquaredModule &operator=(const ChiSquaredModule &);
diff --git a/Core/Algorithms/inc/IChiSquaredModule.h b/Core/Algorithms/inc/IChiSquaredModule.h
index 37d8a3b6ab78bae7ea078470d5da90bb939aec7b..e9730d15c96965df46eb105751770158fc4804bf 100644
--- a/Core/Algorithms/inc/IChiSquaredModule.h
+++ b/Core/Algorithms/inc/IChiSquaredModule.h
@@ -73,6 +73,9 @@ public:
     //! set number of degree of freedom
     void setNdegreeOfFreedom(int ndegree_of_freedom) { m_ndegree_of_freedom = ndegree_of_freedom; }
 
+    //! return residual between data and simulation for single element
+    virtual double getResidualValue(int /* index */ ) const { throw NotImplementedException("IChiSquaredModule::getResidualValue() -> Error! Not implemented."); }
+
 protected:
     // disabling assignment operator
     IChiSquaredModule &operator=(const IChiSquaredModule &);
diff --git a/Core/Algorithms/src/ChiSquaredModule.cpp b/Core/Algorithms/src/ChiSquaredModule.cpp
index bf2764c33d985c1fd55ce6bad7e9b29a7a92573f..c1c7f688ba84d8f24b9b98e5b0786270f2318242 100644
--- a/Core/Algorithms/src/ChiSquaredModule.cpp
+++ b/Core/Algorithms/src/ChiSquaredModule.cpp
@@ -52,6 +52,16 @@ double ChiSquaredModule::calculateChiSquared()
 }
 
 
+// FIXME: ChiSquaredModule::getResidualValue() is implemented very quickly and very dirty (no weights, no protection against missed call to calculateChiSquared, against wrong index)
+double ChiSquaredModule::getResidualValue(int index ) const
+{
+    double value_real = (*mp_real_data)[index];
+    double value_simu  = (*mp_simulation_data)[index];
+    double squared_difference = mp_squared_function->calculateSquaredDifference(value_real, value_simu);
+    return (squared_difference > 0 ? std::sqrt(squared_difference) : 0.0);
+}
+
+
 OutputData<double>* ChiSquaredModule::createChi2DifferenceMap() const
 {
     if( !mp_real_data ) throw NullPointerException("ChiSquaredModule::calculateChiSquared() -> Error! No real data has been set");
diff --git a/Core/Algorithms/src/DiffuseDWBASimulation.cpp b/Core/Algorithms/src/DiffuseDWBASimulation.cpp
index b3932d95eed398c3a464e8b7432a47327fa09c23..d75154c648ff2e4873a4786242b54e8b195c18fd 100644
--- a/Core/Algorithms/src/DiffuseDWBASimulation.cpp
+++ b/Core/Algorithms/src/DiffuseDWBASimulation.cpp
@@ -16,6 +16,9 @@ DiffuseDWBASimulation::~DiffuseDWBASimulation()
 
 void DiffuseDWBASimulation::run()
 {
+    const std::string s_phi_f("phi_f");
+    const std::string s_alpha_f("alpha_f");
+
     std::vector<DiffuseFormFactorTerm *> diffuse_terms;
     size_t nbr_heights = 50;
     size_t samples_per_particle = 9;
@@ -24,8 +27,8 @@ void DiffuseDWBASimulation::run()
 
     DWBASimulation::iterator it_intensity = begin();
     while ( it_intensity != end() ) {
-        double phi_f = getDWBAIntensity().getValueOfAxis("phi_f", it_intensity.getIndex());
-        double alpha_f = getDWBAIntensity().getValueOfAxis("alpha_f", it_intensity.getIndex());
+        double phi_f = getDWBAIntensity().getValueOfAxis(s_phi_f, it_intensity.getIndex());
+        double alpha_f = getDWBAIntensity().getValueOfAxis(s_alpha_f, it_intensity.getIndex());
         if (alpha_f<0) {
             ++it_intensity;
             continue;
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index f400c35ad38ce6267668f00a8071e6fc8f36a906..b900691d727734c738e38f64d305936bd7e77745 100644
--- a/Core/Algorithms/src/GISASExperiment.cpp
+++ b/Core/Algorithms/src/GISASExperiment.cpp
@@ -129,9 +129,12 @@ void GISASExperiment::normalize()
 void GISASExperiment::setDetectorParameters(size_t n_phi, double phi_f_min, double phi_f_max,
                                             size_t n_alpha, double alpha_f_min, double alpha_f_max, bool isgisaxs_style)
 {
+    const std::string s_phi_f("phi_f");
+    const std::string s_alpha_f("alpha_f");
+
     m_detector.clear();
-    AxisDouble phi_axis("phi_f");
-    AxisDouble alpha_axis("alpha_f");
+    AxisDouble phi_axis(s_phi_f);
+    AxisDouble alpha_axis(s_alpha_f);
     if (isgisaxs_style) {
         initializeAnglesIsgisaxs(&phi_axis, phi_f_min, phi_f_max, n_phi);
         initializeAnglesIsgisaxs(&alpha_axis, alpha_f_min, alpha_f_max, n_alpha);
@@ -153,6 +156,9 @@ void GISASExperiment::setDetectorResolutionFunction(IResolutionFunction2D *p_res
 
 void GISASExperiment::smearIntensityFromZAxisTilting()
 {
+    const std::string s_phi_f("phi_f");
+    const std::string s_alpha_f("alpha_f");
+
     size_t nbr_zetas = 5;
     double zeta_sigma = 45.0*Units::degree;
     std::vector<double> zetas;
@@ -163,8 +169,8 @@ void GISASExperiment::smearIntensityFromZAxisTilting()
     m_intensity_map.setAllTo(0.0);
     OutputData<double>::const_iterator it_clone = p_clone->begin();
     while (it_clone != p_clone->end()) {
-        double old_phi = p_clone->getValueOfAxis("phi_f", it_clone.getIndex());
-        double old_alpha = p_clone->getValueOfAxis("alpha_f", it_clone.getIndex());
+        double old_phi = p_clone->getValueOfAxis(s_phi_f, it_clone.getIndex());
+        double old_alpha = p_clone->getValueOfAxis(s_alpha_f, it_clone.getIndex());
         for (size_t zeta_index=0; zeta_index<zetas.size(); ++zeta_index) {
             double newphi = old_phi + deltaPhi(old_alpha, old_phi, zetas[zeta_index]);
             double newalpha = old_alpha + deltaAlpha(old_alpha, zetas[zeta_index]);
@@ -191,8 +197,8 @@ void GISASExperiment::initializeAnglesIsgisaxs(AxisDouble *p_axis, double start,
 
 double GISASExperiment::getSolidAngle(size_t index) const
 {
-    const std::string s_alpha_f("alpha_f");
     const std::string s_phi_f("phi_f");
+    const std::string s_alpha_f("alpha_f");
 
     const AxisDouble *p_alpha_axis = m_intensity_map.getAxis(s_alpha_f);
     const AxisDouble *p_phi_axis = m_intensity_map.getAxis(s_phi_f);
diff --git a/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp b/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
index b529436d54576ee1c7b85967a1016bd8e79da39e..45e04ad8a349ac3e941a6abe705fc68ac550e044 100644
--- a/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
+++ b/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
@@ -80,6 +80,8 @@ std::vector<IFormFactor *> LayerDecoratorDWBASimulation::createDWBAFormFactors()
 
 void LayerDecoratorDWBASimulation::calculateCoherentIntensity(IInterferenceFunctionStrategy *p_strategy)
 {
+    const std::string s_phi_f("phi_f");
+    const std::string s_alpha_f("alpha_f");
     //std::cout << "Calculating coherent scattering..." << std::endl;
     double wavelength = getWaveLength();
     double total_surface_density = mp_layer_decorator->getTotalParticleSurfaceDensity();
@@ -87,8 +89,8 @@ void LayerDecoratorDWBASimulation::calculateCoherentIntensity(IInterferenceFunct
     DWBASimulation::iterator it_intensity = begin();
     while ( it_intensity != end() )
     {
-        double phi_f = getDWBAIntensity().getValueOfAxis("phi_f", it_intensity.getIndex());
-        double alpha_f = getDWBAIntensity().getValueOfAxis("alpha_f", it_intensity.getIndex());
+        double phi_f = getDWBAIntensity().getValueOfAxis(s_phi_f, it_intensity.getIndex());
+        double alpha_f = getDWBAIntensity().getValueOfAxis(s_alpha_f, it_intensity.getIndex());
         if (alpha_f<0) {
             ++it_intensity;
             continue;
diff --git a/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp b/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
index 181e5c7983ac95939b2f6d2ed50387cdb3162edf..a1f714c9557daa619acdbad56f95f296da0f7dbb 100644
--- a/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
+++ b/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
@@ -23,14 +23,17 @@ void MultiLayerRoughnessDWBASimulation::setReflectionTransmissionFunction(size_t
 
 void MultiLayerRoughnessDWBASimulation::run()
 {
+    const std::string s_phi_f("phi_f");
+    const std::string s_alpha_f("alpha_f");
+
     kvector_t m_ki_real(m_ki.x().real(), m_ki.y().real(), m_ki.z().real());
     double lambda = 2.0*M_PI/m_ki_real.mag();
 
     DWBASimulation::iterator it_intensity = begin();
     while ( it_intensity != m_dwba_intensity.end() )
     {
-        double phi_f = getDWBAIntensity().getValueOfAxis("phi_f", it_intensity.getIndex());
-        double alpha_f = getDWBAIntensity().getValueOfAxis("alpha_f", it_intensity.getIndex());
+        double phi_f = getDWBAIntensity().getValueOfAxis(s_phi_f, it_intensity.getIndex());
+        double alpha_f = getDWBAIntensity().getValueOfAxis(s_alpha_f, it_intensity.getIndex());
         cvector_t k_f;
         k_f.setLambdaAlphaPhi(lambda, alpha_f, phi_f);
         *it_intensity = evaluate(m_ki, k_f, -m_alpha_i, alpha_f);
diff --git a/Core/Tools/inc/AxisDouble.h b/Core/Tools/inc/AxisDouble.h
index dfa456f9cc5597499c3332679abf64c9da2f4f13..3107cf34b75996bd145c5e27afbe76699f3c2fea 100644
--- a/Core/Tools/inc/AxisDouble.h
+++ b/Core/Tools/inc/AxisDouble.h
@@ -39,28 +39,28 @@ public:
     ~AxisDouble() {}
 
     //! retrieve the number of bins
-    size_t getSize() const { return m_value_vector.size(); }
+    inline size_t getSize() const { return m_value_vector.size(); }
 
     //! retrieve the label of the axis
-    std::string getName() const { return m_name; }
+    inline std::string getName() const { return m_name; }
 
     //! set the axis label
-    void setName(std::string name) { m_name = name; }
+    inline void setName(std::string name) { m_name = name; }
 
     //! add point to the end of the axis
-    void push_back(double element) { m_value_vector.push_back(element); }
+    inline void push_back(double element) { m_value_vector.push_back(element); }
 
     //! indexed accessor
-    double& operator[](size_t index) { return m_value_vector.at(index); }
+    inline double& operator[](size_t index) { return m_value_vector.at(index); }
 
     //! const indexed accessor
-    const double& operator[](size_t index) const { return m_value_vector.at(index); }
+    inline const double& operator[](size_t index) const { return m_value_vector.at(index); }
 
     //! get value of first point of axis
-    double getMin() const { return m_value_vector.front(); }
+    inline double getMin() const { return m_value_vector.front(); }
 
     //! get value of last point of axis
-    double getMax() const { return m_value_vector.back(); }
+    inline double getMax() const { return m_value_vector.back(); }
 
     //! initialize axis points
     void initElements(size_t size, double start, double end);
diff --git a/Core/Tools/inc/FitSuite.h b/Core/Tools/inc/FitSuite.h
index 49833193e842da525d6e98c4b5c28911dbdea4cf..f55fba07feba3153632092514906b8b80f337ec9 100644
--- a/Core/Tools/inc/FitSuite.h
+++ b/Core/Tools/inc/FitSuite.h
@@ -72,6 +72,9 @@ public:
     //! function to minimize
     double functionToMinimize(const double *pars_current_values);
 
+    //! provides minimizer with gradients wrt parameters for single data element
+    double elementFunction(const double *pars_current_values, unsigned int index, double *deriv);
+
     //! return reference to the kit with data
     FitSuiteObjects *getFitObjects() { return &m_fit_objects; }
 
diff --git a/Core/Tools/inc/FitSuiteObjects.h b/Core/Tools/inc/FitSuiteObjects.h
index 02a648ad80c665f7862d8591d246ce94e87e4fd8..d9aaf10b66a3fe23cc399745f50d3ae231c629ee 100644
--- a/Core/Tools/inc/FitSuiteObjects.h
+++ b/Core/Tools/inc/FitSuiteObjects.h
@@ -50,6 +50,9 @@ public:
     //! get sum of chi squared values for all fit objects
     double getChiSquaredValue(int n_free_fit_parameters = 0);
 
+    //! get residuals for single data element
+    double getResidualValue(int index);
+
     //! get experiment
     const Experiment *getExperiment(size_t i_item = 0) const { return m_fit_objects[check_index(i_item)]->getExperiment(); }
 
diff --git a/Core/Tools/inc/IMinimizer.h b/Core/Tools/inc/IMinimizer.h
index f22fd3b63ae86449c9abc5e96b749ee7cb53b0be..b25f3cabaa354a1834a5bffee965f4b2a342c85a 100644
--- a/Core/Tools/inc/IMinimizer.h
+++ b/Core/Tools/inc/IMinimizer.h
@@ -28,6 +28,11 @@
 class IMinimizer
 {
 public:
+    //! signature of function to minimize
+    typedef boost::function<double(const double *)> function_t;
+    //! signature of function to minimize with acess to single element residual
+    typedef boost::function<double(const double *, unsigned int, double *)> element_function_t;
+
     IMinimizer(){}
     virtual ~IMinimizer(){}
 
@@ -35,7 +40,9 @@ public:
     virtual void setVariable(int index, const FitParameter *par) = 0;
 
     //! set function to minimize
-    virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1) = 0;
+//    virtual void setFunction(function_t fcn, int ndims) = 0;
+//    virtual void setElementFunction(element_function_t fcn, int ndims, int nelements) = 0;
+    virtual void setFunction(function_t fcn, int ndims, element_function_t element_fcn = element_function_t(), int nelements = 0) = 0;
 
     //! run minimization
     virtual void minimize() = 0;
@@ -76,7 +83,10 @@ public:
     virtual void setVariable(int index, const FitParameter *par) { m_values[index] = par->getValue(); }
 
     //! set function to minimize
-    virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1) { m_fcn = fcn, m_ndim = ndim; }
+    //virtual void setFunction(function_t fcn, int ndims) { m_fcn = fcn; (void)ndims; }
+    virtual void setFunction(function_t fcn, int /* ndims */, element_function_t /* element_fcn */, int /* nelements */ ) { m_fcn = fcn; }
+
+//    virtual void setElementFunction(element_function_t /* fcn */, int /* ndims */, int /* nelements */) { throw NotImplementedException("TestMinimizer::setGradientFunction"); }
 
     //! run minimization
     virtual void minimize()
@@ -119,8 +129,7 @@ public:
 
 private:
     std::map<int, double > m_values;
-    boost::function<double(const double *)> m_fcn;
-    int m_ndim;
+    function_t m_fcn;
 };
 
 
diff --git a/Core/Tools/inc/LLData.h b/Core/Tools/inc/LLData.h
index 98b70a512833c545fcf8602afd2a53c12150c422..243f8714e534e559f74ae5651e61808c3ea7fb9e 100644
--- a/Core/Tools/inc/LLData.h
+++ b/Core/Tools/inc/LLData.h
@@ -51,7 +51,7 @@ public:
 
     // retrieve basic info
     size_t getTotalSize() const;
-    size_t getRank() const { return m_rank; }
+    inline size_t getRank() const { return m_rank; }
     const int *getDimensions() const { return m_dims; }
     T getTotalSum() const;
 
diff --git a/Core/Tools/inc/OutputData.h b/Core/Tools/inc/OutputData.h
index bfcdd8c61bd9853377fad7996f77a7e9c3acba30..37fdf160bc9272d93f1d2c568453fb247f719bf9 100644
--- a/Core/Tools/inc/OutputData.h
+++ b/Core/Tools/inc/OutputData.h
@@ -124,14 +124,17 @@ public:
     //! return vector of coordinates for given index
     std::vector<int> toCoordinates(size_t index) const;
 
+    //! return coordinate for given index and axis number
+    int toCoordinate(size_t index, size_t i_selected_axis) const;
+
     //! return index for specified coordinates
     size_t toIndex(std::vector<int> coordinates) const;
 
     //! return index of axis with given name for given total index
-    size_t getIndexOfAxis(std::string axis_name, size_t total_index) const;
+    size_t getIndexOfAxis(const std::string &axis_name, size_t total_index) const;
 
     //! return value of axis with given name at given index
-    double getValueOfAxis(std::string axis_name, size_t index) const;
+    double getValueOfAxis(const std::string &axis_name, size_t index) const;
 
     // ---------
     // modifiers
@@ -380,6 +383,20 @@ template<class T> std::vector<int> OutputData<T>::toCoordinates(size_t index) co
     return result;
 }
 
+template<class T> int OutputData<T>::toCoordinate(size_t index, size_t i_selected_axis) const
+{
+    size_t remainder(index);
+    for (size_t i=0; i<mp_ll_data->getRank(); ++i)
+    {
+        size_t i_axis = mp_ll_data->getRank()-1-i;
+        int result = (int)(remainder % m_value_axes[i_axis].getSize());
+        if(i_selected_axis == i_axis ) return result;
+        remainder /= m_value_axes[i_axis].getSize();
+    }
+    throw LogicErrorException("OutputData<T>::toCoordinate() -> Error! No axis with given number");
+}
+
+
 template <class T> size_t OutputData<T>::toIndex(std::vector<int> coordinates) const
 {
     if (coordinates.size() != mp_ll_data->getRank()) {
@@ -395,7 +412,7 @@ template <class T> size_t OutputData<T>::toIndex(std::vector<int> coordinates) c
     return result;
 }
 
-template <class T> size_t OutputData<T>::getIndexOfAxis(std::string axis_name, size_t total_index) const
+template <class T> size_t OutputData<T>::getIndexOfAxis(const std::string &axis_name, size_t total_index) const
 {
     std::vector<int> coordinates = toCoordinates(total_index);
     for (size_t i=0; i<m_value_axes.size(); ++i) {
@@ -403,22 +420,36 @@ template <class T> size_t OutputData<T>::getIndexOfAxis(std::string axis_name, s
             return coordinates[i];
         }
     }
-    throw LogicErrorException("OutputData<T>::getIndexOfAxis() -> Error! Axis with given name not found '"+axis_name+std::string("'"));
+    throw LogicErrorException("OutputData<T>::getIndexOfAxis() -> Error! Axis with given name not found '"+std::string(axis_name)+std::string("'"));
 }
 
+//template <class T>
+//double OutputData<T>::getValueOfAxis(const char *axis_name, size_t index) const
+////double OutputData<T>::getValueOfAxis(std::string axis_name, size_t index) const
+//{
+//    std::vector<int> coordinates = toCoordinates(index);
+//    for (size_t i=0; i<m_value_axes.size(); ++i) {
+//        if (m_value_axes[i].getName() == axis_name) {
+//            return m_value_axes[i][coordinates[i]];
+//        }
+//    }
+//    throw LogicErrorException("OutputData<T>::getValueOfAxis() -> Error! Axis with given name not found '"+std::string(axis_name)+std::string("'"));
+//}
+
 
 template <class T>
-double OutputData<T>::getValueOfAxis(std::string axis_name, size_t index) const
+double OutputData<T>::getValueOfAxis(const std::string &axis_name, size_t index) const
 {
-    std::vector<int> coordinates = toCoordinates(index);
     for (size_t i=0; i<m_value_axes.size(); ++i) {
         if (m_value_axes[i].getName() == axis_name) {
-            return m_value_axes[i][coordinates[i]];
+            int axis_index = toCoordinate(index, i);
+            return m_value_axes[i][axis_index];
         }
     }
-    throw LogicErrorException("OutputData<T>::getValueOfAxis() -> Error! Axis with given name not found '"+axis_name+std::string("'"));
+    throw LogicErrorException("OutputData<T>::getValueOfAxis() -> Error! Axis with given name not found '"+std::string(axis_name)+std::string("'"));
 }
 
+
 template<class T>
 inline T OutputData<T>::totalSum() const
 {
diff --git a/Core/Tools/src/FitSuite.cpp b/Core/Tools/src/FitSuite.cpp
index bde80ad67448bb644b458968543870cf9abf1832..74b045b5c56322e0e28c8bb46fb06e2bd23058ba 100644
--- a/Core/Tools/src/FitSuite.cpp
+++ b/Core/Tools/src/FitSuite.cpp
@@ -84,13 +84,24 @@ void FitSuite::link_fit_parameters()
 void FitSuite::minimize()
 {
     // initializing minimizer with fcn function belonging to given class
-    m_minimizer->setFunction( boost::bind(&FitSuite::functionToMinimize, this, _1), (int)m_fit_parameters.size() );
-
-    // propagating local fit parameters to the minimizer's internal list of parameters
+    IMinimizer::function_t fcn = boost::bind(&FitSuite::functionToMinimize, this, _1);
+    //m_minimizer->setFunction( fcn, (int)m_fit_parameters.size() );
+    // FIXME: FitSuite::minimize() where to get number of elements?
+    int nelements = m_fit_objects.getRealData()->getAllocatedSize();
+    IMinimizer::element_function_t element_fcn = boost::bind(&FitSuite::elementFunction, this, _1, _2, _3);
+    m_minimizer->setFunction( fcn, (int)m_fit_parameters.size(), element_fcn, nelements );
+
+    // propagating values of local fit parameters to the minimizer's internal parameters
     for(size_t i_par = 0; i_par<m_fit_parameters.size(); i_par++) {
+        std::cout << " i_par " << i_par << std::endl;
         m_minimizer->setVariable((int)i_par, m_fit_parameters[i_par] );
     }
-    if( m_fit_parameters.size() != m_minimizer->getNumberOfVariables())  std::cout << "FitSuite::minimize() -> Warning. Something unexpected" << std::endl;
+    if( m_fit_parameters.size() != m_minimizer->getNumberOfVariables())  {
+        std::ostringstream ostr;
+        ostr << "FitSuite::minimize() -> Error! Number of variables defined in minimizer (" << m_minimizer->getNumberOfVariables() << ") ";
+        ostr << "doesn't coincide with number of FitSuite's parameters (" << m_fit_parameters.size() << ")";
+        throw LogicErrorException(ostr.str());
+    }
 
     // minimizing
     m_minimizer->minimize();
@@ -166,3 +177,31 @@ double FitSuite::functionToMinimize(const double *pars_current_values)
 }
 
 
+/* ************************************************************************* */
+// provides minimizer with gradients wrt parameters for single data element
+/* ************************************************************************* */
+double FitSuite::elementFunction(const double *pars_current_values, unsigned int index, double *deriv)
+{
+    if(index%10 == 0) std::cout << "elementFunction "  << index << std::endl;
+
+    // set fitting parameters to values suggested by the minimizer
+    m_fit_parameters.setValues(pars_current_values);
+
+    // run simulations
+    m_fit_objects.runSimulation();
+
+    // caclulate residual value
+//    m_fit_objects.getChiSquaredValue(m_fit_parameters.getNfreeParameters());
+    double residual = m_fit_objects.getResidualValue(index);
+
+    if( deriv ) {
+        throw NotImplementedException("FitSuite::elementFunction() -> Error! Calculation of derivatives is not implemented");
+    }
+
+    if(index == 100 ) assert(0);
+
+    return residual;
+}
+
+
+
diff --git a/Core/Tools/src/FitSuiteObjects.cpp b/Core/Tools/src/FitSuiteObjects.cpp
index 2e1493edd87cd0bdba43b1c773cb5092e9d5893f..90a01a1f9b92b191de042d59027e5d7c69ce7864 100644
--- a/Core/Tools/src/FitSuiteObjects.cpp
+++ b/Core/Tools/src/FitSuiteObjects.cpp
@@ -43,10 +43,11 @@ void FitSuiteObjects::runSimulation()
 
 /* ************************************************************************* */
 // get sum of chi squared values for all fit objects
-// FIXME: refactor FitSuiteObjects::getChiSquaredValue()
+// FIXME: refactor FitSuiteObjects::getChiSquaredValue() (the main problem is duplication calculateChiSquared() for ChiSquaredModule and FitObject)
 /* ************************************************************************* */
 double FitSuiteObjects::getChiSquaredValue(int n_free_fit_parameters)
 {
+    double max_intensity = getSimulationMaxIntensity();
     double chi_sum(0);
     for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
         IChiSquaredModule *chi = (*it)->getChiSquaredModule();
@@ -54,7 +55,7 @@ double FitSuiteObjects::getChiSquaredValue(int n_free_fit_parameters)
         chi->setNdegreeOfFreedom( m_fit_objects.size() * (*it)->getRealData()->getAllocatedSize() - n_free_fit_parameters);
         // normalizing datasets to the maximum intensity over all fit objects defined
         OutputDataNormalizerScaleAndShift *data_normalizer =  dynamic_cast<OutputDataNormalizerScaleAndShift *>(chi->getOutputDataNormalizer());
-        if( data_normalizer) data_normalizer->setMaximumIntensity( getSimulationMaxIntensity() );
+        if( data_normalizer) data_normalizer->setMaximumIntensity( max_intensity );
 
         double weight = (*it)->getWeight()/m_total_weight;
         double chi_squared = (weight*weight) * (*it)->calculateChiSquared();
@@ -65,8 +66,28 @@ double FitSuiteObjects::getChiSquaredValue(int n_free_fit_parameters)
 }
 
 
+double FitSuiteObjects::getResidualValue(int index)
+{
+    double residual_sum(0);
+    for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
+        IChiSquaredModule *chi = (*it)->getChiSquaredModule();
+        const OutputData<double> *data_real = (*it)->getRealData();
+        const OutputData<double> *data_simu = (*it)->getSimulationData();
+        double value_real = (*data_real)[index];
+        double value_simu = (*data_simu)[index];
+        double squared_difference = chi->getSquaredFunction()->calculateSquaredDifference(value_real, value_simu);
+        double weight = (*it)->getWeight()/m_total_weight;
+        double residual(0);
+        (squared_difference > 0 ? residual = weight*std::sqrt(squared_difference) : 0.0);
+        residual_sum += residual;
+    }
+    return residual_sum;
+}
+
+
+
 /* ************************************************************************* */
-// calculate maximum intensity in simulated data over all fir objects defined
+// calculate maximum intensity in simulated data over all fit objects defined
 /* ************************************************************************* */
 double FitSuiteObjects::getSimulationMaxIntensity()
 {
diff --git a/UnitTests/TestCore/OutputDataTest.h b/UnitTests/TestCore/OutputDataTest.h
index c6bd0714fe30aca37fae8918eef77dc651aea7db..90c7e83cf49f4f6441fc1e634dabdef1c1a88e63 100644
--- a/UnitTests/TestCore/OutputDataTest.h
+++ b/UnitTests/TestCore/OutputDataTest.h
@@ -114,4 +114,20 @@ TEST_F(OutputDataTest, MaxElement)
     EXPECT_EQ( double(10.0), (*cit));
 }
 
+TEST_F(OutputDataTest, ValueOfAxis)
+{
+    OutputData<double > data;
+    data.addAxis("axis1", 10, 0., 9.);
+    data.addAxis("axis2", 2, 0., 10.);
+    EXPECT_EQ( double(0.0), data.getValueOfAxis("axis1", 0));
+    EXPECT_EQ( double(0.0), data.getValueOfAxis("axis2", 0));
+    EXPECT_EQ( double(0.0), data.getValueOfAxis("axis1", 1));
+    EXPECT_EQ( double(10.0), data.getValueOfAxis("axis2", 1));
+    EXPECT_EQ( double(1.0), data.getValueOfAxis("axis1", 2));
+    EXPECT_EQ( double(0.0), data.getValueOfAxis("axis2", 2));
+    EXPECT_EQ( double(9.0), data.getValueOfAxis("axis1", 19));
+    EXPECT_EQ( double(10.0), data.getValueOfAxis("axis2", 19));
+}
+
+
 #endif // OUTPUTDATATEST_H