diff --git a/App/inc/TestFresnelCoeff.h b/App/inc/TestFresnelCoeff.h
index 0f37d3e6fd06c8368ad6dc2936557b997f105a36..3a34bfdb7a52c8c43bee487ae2b374ab542203af 100644
--- a/App/inc/TestFresnelCoeff.h
+++ b/App/inc/TestFresnelCoeff.h
@@ -43,15 +43,8 @@ private:
     //! draw result of the test
     void draw_roughness_set();
 
-    MultiLayer *m_sample; //!< pointer to multilayer sample
-    OutputData<OpticalFresnel::MultiLayerCoeff_t  > *m_coeffs; //!< output data structure
+    MultiLayer *mp_sample; //!< pointer to multilayer sample
+    OutputData<OpticalFresnel::MultiLayerCoeff_t  > *mp_coeffs; //!< output data structure
 };
 
-
-
-
-
-
-
-
 #endif // TESTFRESNELCOEFF_H
diff --git a/App/inc/TestRootTree.h b/App/inc/TestRootTree.h
index 4d63a7a65e436b91ea4c622dbf5230647369363b..ef8af430a6305afe819fc8fa6e0ce82c087c3c99 100644
--- a/App/inc/TestRootTree.h
+++ b/App/inc/TestRootTree.h
@@ -56,9 +56,9 @@ private:
 
     void initializeMesoCrystal(double meso_alpha, double meso_phi, double nanopart_radius);
 
-    ISample *m_sample;
-    GISASExperiment *m_experiment;
-    OutputData<double> *m_data;
+    ISample *mp_sample;
+    GISASExperiment *mp_experiment;
+    OutputData<double> *mp_data;
 };
 
 
diff --git a/App/src/FitSuiteHelper.cpp b/App/src/FitSuiteHelper.cpp
index 00ab605df7a2e018f904d1a62a7f1f588c40fc89..a8ca6843c77401aef8f48f69adaf71b0c0d9137b 100644
--- a/App/src/FitSuiteHelper.cpp
+++ b/App/src/FitSuiteHelper.cpp
@@ -156,25 +156,25 @@ void FitSuiteObserverDraw::update(IObservable *subject)
 // return output data which contains relative difference between simulation and real data
 OutputData<double > *FitSuiteObserverDraw::getRelativeDifferenceMap(const IChiSquaredModule *chi_module)
 {
-    const OutputData<double> *simu_data = chi_module->getSimulationData();
-    const OutputData<double> *real_data = chi_module->getRealData();
-
-    OutputData<double > *difference_map = simu_data->clone();
-    difference_map->setAllTo(0.0);
-
-    simu_data->resetIndex();
-    real_data->resetIndex();
-    difference_map->resetIndex();
-    while (real_data->hasNext()) {
-        double value_simu = simu_data->currentValue();
-        double value_real = real_data->currentValue();
+    const OutputData<double> *p_simu_data = chi_module->getSimulationData();
+    const OutputData<double> *p_real_data = chi_module->getRealData();
+
+    OutputData<double > *p_difference_map = p_simu_data->clone();
+    p_difference_map->setAllTo(0.0);
+
+    OutputData<double>::const_iterator it_sim = p_simu_data->begin();
+    OutputData<double>::const_iterator it_real = p_real_data->begin();
+    OutputData<double>::iterator it_diff = p_difference_map->begin();
+
+    while (it_sim != p_simu_data->end()) {
+        double value_simu = *it_sim++;
+        double value_real = *it_real++;
         double value_diff(0);
         if( value_real > 0) value_diff = std::abs(value_real - value_simu)/value_real;
-        difference_map->next() = value_diff;
-        real_data->next(); simu_data->next();
+        *it_diff = value_diff;
+        ++it_diff;
     }
-
-    return difference_map;
+    return p_difference_map;
 }
 
 
diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index 33db0464333949d25f3d9ade27719d8d585ffeed..d6ae41db89c1621d30513967107a42a3b9d3d341 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -151,12 +151,12 @@ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const s
     h2->GetXaxis()->SetTitle( axis0->getName().c_str() );
     h2->GetYaxis()->SetTitle( axis1->getName().c_str() );
 
-    output.resetIndex();
-    while (output.hasNext())
+    OutputData<double>::const_iterator it = output.begin();
+    while (it != output.end())
     {
-        double x_value = output.getCurrentValueOfAxis<double>( axis0->getName().c_str() );
-        double y_value = output.getCurrentValueOfAxis<double>( axis1->getName().c_str() );
-        double z_value = output.next();
+        double x_value = output.getValueOfAxis<double>( axis0->getName().c_str(), it.getIndex() );
+        double y_value = output.getValueOfAxis<double>( axis1->getName().c_str(), it.getIndex() );
+        double z_value = *it++;
         h2->Fill(x_value, y_value, z_value);
     }
     h2->SetContour(50);
@@ -188,10 +188,10 @@ void IsGISAXSTools::drawOutputDataDistribution1D(const OutputData<double> &outpu
     TH1::SetDefaultBufferSize(5000);
     TH1D h1_spect("h1_spect", histo_name.c_str(), 200, 1.0, -1.0); // xmin > xmax as a sign of automatic binning
 
-    output.resetIndex();
-    while (output.hasNext())
+    OutputData<double>::const_iterator it = output.begin();
+    while (it != output.end())
     {
-        h1_spect.Fill( output.next() );
+        h1_spect.Fill( *it++ );
     }
 
     h1_spect.DrawCopy(draw_options.c_str());
@@ -221,14 +221,14 @@ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, c
     TH1D h1_spect("difference", histo_name.c_str(), 40, -20.0, 20.0);
     h1_spect.GetXaxis()->SetTitle("log10( (we-isgi)/isgi )");
 
-    left_clone->resetIndex();
-    while (left_clone->hasNext())
+    OutputData<double>::const_iterator it = left_clone->begin();
+    while (it != left_clone->end())
     {
-        double x = left_clone->next();
+        double x = *it++;
         if(x!=0) {
             x = log10(fabs(x));
         } else {
-            // lets put the cases then the difference is exactly 0 to underflow bin
+            // lets put the cases when the difference is exactly 0 to underflow bin
             x = -21.;
         }
         h1_spect.Fill( x );
@@ -308,16 +308,14 @@ void IsGISAXSTools::writeOutputDataToFile(const OutputData<double>& output,
         return;
         //throw FileNotIsOpenException("IsGISAXSTools::writeOutputDataToFile() -> Error. Can't open file '"+filename+"' for writing.");
     }
-    output.resetIndex();
     size_t row_length = output.getAxes()[1]->getSize();
-    int counter = 1;
-    while(output.hasNext()) {
-        double z_value = output.next();
+    OutputData<double>::const_iterator it = output.begin();
+    while(it != output.end()) {
+        double z_value = *it++;
         file << std::scientific << std::setprecision(precision) << z_value << "    ";
-        if(counter%row_length==0) {
+        if(it.getIndex()%row_length==0) {
             file << std::endl;
         }
-        ++counter;
     }
     if ( file.bad() ) {
         throw FileIsBadException("IsGISAXSTools::writeOutputDataToFile() -> Error! File is bad, probably it is a directory.");
@@ -377,21 +375,21 @@ OutputData<double> *IsGISAXSTools::readOutputDataFromFile(const std::string &fil
     // creating new OutputData and filling it with values from buffer_2d
     int y_size = buff_2d.size();
     int x_size = buff_2d.size() ? buff_2d[0].size() : 0;
-    OutputData<double> *output = new OutputData<double>;
-    output->addAxis(std::string("x-axis"), 0.0, double(x_size), x_size);
-    output->addAxis(std::string("y-axis"), 0.0, double(y_size), y_size);
-    output->setAllTo(0.0);
+    OutputData<double> *p_result = new OutputData<double>;
+    p_result->addAxis(std::string("x-axis"), 0.0, double(x_size), x_size);
+    p_result->addAxis(std::string("y-axis"), 0.0, double(y_size), y_size);
+    p_result->setAllTo(0.0);
 
-    output->resetIndex();
-    while (output->hasNext())
+    OutputData<double>::iterator it = p_result->begin();
+    while (it != p_result->end())
     {
-        size_t index_x = output->getCurrentIndexOfAxis("x-axis");
-        size_t index_y = output->getCurrentIndexOfAxis("y-axis");
-        //output->next() = buff_2d[index_y][index_x];
-        output->next() = buff_2d[index_x][index_y]; // strange
+        size_t index_x = p_result->toCoordinates(it.getIndex())[0];
+        size_t index_y = p_result->toCoordinates(it.getIndex())[1];
+        *it = buff_2d[index_x][index_y];
+        ++it;
     }
 
-    return output;
+    return p_result;
 }
 
 
@@ -402,7 +400,6 @@ void IsGISAXSTools::exportOutputDataInVectors2D(const OutputData<double> &output
 {
     if (output_data.getRank() != 2) return;
 
-    output_data.resetIndex();
     const NamedVector<double> *p_axis0 = dynamic_cast<const NamedVector<double>*>(output_data.getAxes()[0]);
     const NamedVector<double> *p_axis1 = dynamic_cast<const NamedVector<double>*>(output_data.getAxes()[1]);
     std::string axis0_name = p_axis0->getName();
@@ -424,13 +421,14 @@ void IsGISAXSTools::exportOutputDataInVectors2D(const OutputData<double> &output
         v_axis1[i].resize(axis1_size,0.0);
     }
 
-    while (output_data.hasNext())
+    OutputData<double>::const_iterator it = output_data.begin();
+    while (it != output_data.end())
     {
-        size_t axis0_index = output_data.getCurrentIndexOfAxis(axis0_name.c_str());
-        size_t axis1_index = output_data.getCurrentIndexOfAxis(axis1_name.c_str());
+        size_t axis0_index = output_data.toCoordinates(it.getIndex())[0];
+        size_t axis1_index = output_data.toCoordinates(it.getIndex())[1];
         double axis0_value = (*p_axis0)[axis0_index];
         double axis1_value = (*p_axis1)[axis1_index];
-        double intensity = output_data.next();
+        double intensity = *it++;
 
         v_data[axis0_index][axis1_index] = intensity;
         v_axis0[axis0_index][axis1_index] = axis0_value;
diff --git a/App/src/TestDiffuseReflection.cpp b/App/src/TestDiffuseReflection.cpp
index 078caaef358dc634a3486062a4e82ef026f4a418..ec856012fd2a846380f6c06dbdcf537c1e5c589c 100644
--- a/App/src/TestDiffuseReflection.cpp
+++ b/App/src/TestDiffuseReflection.cpp
@@ -81,12 +81,12 @@ void TestDiffuseReflection::execute()
         m_data_offspec->addAxis(std::string("alpha_i"), m_alphaMin, m_alphaMax, m_npoints);
         m_data_offspec->addAxis(std::string("alpha_f"), m_alphaMin, m_alphaMax, m_npoints);
 
-        m_data_offspec->resetIndex();
-        while (m_data_offspec->hasNext()) {
-            double alpha_i = m_data_offspec->getCurrentValueOfAxis<double>("alpha_i");
-            double alpha_f = m_data_offspec->getCurrentValueOfAxis<double>("alpha_f");
-            size_t index_alpha_i = m_data_offspec->getCurrentIndexOfAxis("alpha_i");
-            size_t index_alpha_f = m_data_offspec->getCurrentIndexOfAxis("alpha_f");
+        OutputData<double>::iterator it = m_data_offspec->begin();
+        while (it != m_data_offspec->end()) {
+            double alpha_i = m_data_offspec->getValueOfAxis<double>("alpha_i", it.getIndex());
+            double alpha_f = m_data_offspec->getValueOfAxis<double>("alpha_f", it.getIndex());
+            size_t index_alpha_i = m_data_offspec->getIndexOfAxis("alpha_i", it.getIndex());
+            size_t index_alpha_f = m_data_offspec->getIndexOfAxis("alpha_f", it.getIndex());
             ki.setLambdaAlphaPhi(1.54*Units::angstrom, -alpha_i, 0.0);
             kf.setLambdaAlphaPhi(1.54*Units::angstrom, alpha_f, 0.0);
             calc.execute(*m_sample, ki, kf);
@@ -98,7 +98,8 @@ void TestDiffuseReflection::execute()
             }
             //double intensity = rdwba.evaluate(ki, kf);
             //std::cout << "alpha_i " << alpha_i << " alpha_f " << alpha_f << " phi_f " << 0.0 << " inten " << intensity << std::endl;
-            m_data_offspec->next() = intensity;
+            *it = intensity;
+            ++it;
         }
 
         draw();
@@ -173,13 +174,14 @@ void TestDiffuseReflection::draw()
     h2.SetContour(50);
 //    h2.SetMinimum(0.001);
     h2.SetStats(0);
-    m_data_offspec->resetIndex();
-    while (m_data_offspec->hasNext()) {
-        double alpha_i = m_data_offspec->getCurrentValueOfAxis<double>("alpha_i");
-        double alpha_f = m_data_offspec->getCurrentValueOfAxis<double>("alpha_f");
-        size_t index_alpha_i = m_data_offspec->getCurrentIndexOfAxis("alpha_i");
-        size_t index_alpha_f = m_data_offspec->getCurrentIndexOfAxis("alpha_f");
-        double r = m_data_offspec->next();
+
+    OutputData<double>::const_iterator it = m_data_offspec->begin();
+    while (it != m_data_offspec->end()) {
+        double alpha_i = m_data_offspec->getValueOfAxis<double>("alpha_i", it.getIndex());
+        double alpha_f = m_data_offspec->getValueOfAxis<double>("alpha_f", it.getIndex());
+        size_t index_alpha_i = m_data_offspec->getIndexOfAxis("alpha_i", it.getIndex());
+        size_t index_alpha_f = m_data_offspec->getIndexOfAxis("alpha_f", it.getIndex());
+        double r = *it++;
         hspect.Fill(r);
         if(index_alpha_i==5) {
 
diff --git a/App/src/TestFittingModule1.cpp b/App/src/TestFittingModule1.cpp
index f11dc976f0816f390c9e48f92b87aeae4ff22b24..8718369a41777650237bb96dca9823301a4b98e6 100644
--- a/App/src/TestFittingModule1.cpp
+++ b/App/src/TestFittingModule1.cpp
@@ -200,13 +200,14 @@ void TestFittingModule1::generateRealData(double noise_factor)
     if (mp_real_data) delete mp_real_data;
 
     mp_real_data = mp_exact_data->clone();
-    mp_real_data->resetIndex();
-    while (mp_real_data->hasNext()) {
-        double current = mp_real_data->currentValue();
+    OutputData<double>::iterator it = mp_real_data->begin();
+    while (it != mp_real_data->end()) {
+        double current = *it;
         double sigma = noise_factor*std::sqrt(current);
         double random = MathFunctions::GenerateNormalRandom(current, sigma);
         if (random<0.0) random = 0.0;
-        mp_real_data->next() = random;
+        *it = random;
+        ++it;
     }
 }
 
diff --git a/App/src/TestFittingModule2.cpp b/App/src/TestFittingModule2.cpp
index a5c19dd7f64e351b8a2c87ed55cebf3cb095b2b8..2eccaa8468ca0e221d3afcd4e4e05a6873f56673 100644
--- a/App/src/TestFittingModule2.cpp
+++ b/App/src/TestFittingModule2.cpp
@@ -203,15 +203,15 @@ void TestFittingModule2::generateRealData(double noise_factor)
     if (mp_real_data) delete mp_real_data;
 
     mp_real_data = mp_exact_data->clone();
-    mp_real_data->resetIndex();
-    while (mp_real_data->hasNext()) {
-        double current = mp_real_data->currentValue();
+    OutputData<double>::iterator it = mp_real_data->begin();
+    while (it != mp_real_data->end()) {
+        double current = *it;
         double sigma = noise_factor*std::sqrt(current);
         double random = MathFunctions::GenerateNormalRandom(current, sigma);
         if (random<0.0) random = 0.0;
-        mp_real_data->next() = random;
+        *it = random;
+        ++it;
     }
-
 }
 
 
diff --git a/App/src/TestFormFactor.cpp b/App/src/TestFormFactor.cpp
index 2d8d4cb0614f40342258659c538d828b859f6f6a..b289ea429ca0f89bf4793caf2cc746707d541f5a 100644
--- a/App/src/TestFormFactor.cpp
+++ b/App/src/TestFormFactor.cpp
@@ -26,23 +26,23 @@ TestFormFactor::~TestFormFactor()
 
 void TestFormFactor::execute()
 {
-    MultiIndex& index = mp_intensity_output->getIndex();
+    OutputData<double>::iterator it = mp_intensity_output->begin();
     const NamedVector<double> *p_y_axis = dynamic_cast<const  NamedVector<double>*>(mp_intensity_output->getAxis("detector y-axis"));
     const NamedVector<double> *p_z_axis = dynamic_cast<const  NamedVector<double>*>(mp_intensity_output->getAxis("detector z-axis"));
     double lambda = 1.0;
     double alpha_i = 0.2*M_PI/180.0;
     cvector_t k_i;
     k_i.setLambdaAlphaPhi(lambda, -alpha_i, 0.0);
-    while (!index.endPassed())
+    while (it != mp_intensity_output->end())
     {
-        size_t index_y = index.getCurrentIndexOfAxis("detector y-axis");
-        size_t index_z = index.getCurrentIndexOfAxis("detector z-axis");
+        size_t index_y = mp_intensity_output->getIndexOfAxis("detector y-axis", it.getIndex());
+        size_t index_z = mp_intensity_output->getIndexOfAxis("detector z-axis", it.getIndex());
         double phi_f = M_PI*(*p_y_axis)[index_y]/180.0;
         double alpha_f = M_PI*(*p_z_axis)[index_z]/180.0;
         cvector_t k_f;
         k_f.setLambdaAlphaPhi(lambda, alpha_f, phi_f);
-        mp_intensity_output->currentValue() = std::pow(std::abs(m_ff.evaluate(k_i, k_f, alpha_i, alpha_f)),2);
-        ++index;
+        *it = std::pow(std::abs(m_ff.evaluate(k_i, k_f, alpha_i, alpha_f)),2);
+        ++it;
     }
     draw();
 }
@@ -53,8 +53,6 @@ void TestFormFactor::draw()
     TCanvas *c1 = new TCanvas("c1_test_formfactor", "Cylinder Formfactor", 0, 0, 1024, 768);
     (void)c1;
 
-    MultiIndex& index = mp_intensity_output->getIndex();
-    index.reset();
     const NamedVector<double> *p_y_axis = dynamic_cast<const NamedVector<double>*>(mp_intensity_output->getAxis("detector y-axis"));
     const NamedVector<double> *p_z_axis = dynamic_cast<const NamedVector<double>*>(mp_intensity_output->getAxis("detector z-axis"));
     size_t y_size = p_y_axis->getSize();
@@ -68,15 +66,16 @@ void TestFormFactor::draw()
     p_hist2D->GetXaxis()->SetTitle("phi_f");
     p_hist2D->GetYaxis()->SetTitle("alpha_f");
 
-    while (!index.endPassed())
+    OutputData<double>::const_iterator it = mp_intensity_output->begin();
+    while (it != mp_intensity_output->end())
     {
-        size_t index_y = index.getCurrentIndexOfAxis("detector y-axis");
-        size_t index_z = index.getCurrentIndexOfAxis("detector z-axis");
+        size_t index_y = mp_intensity_output->getIndexOfAxis("detector y-axis", it.getIndex());
+        size_t index_z = mp_intensity_output->getIndexOfAxis("detector z-axis", it.getIndex());
         double x_value = (*p_y_axis)[index_y];
         double y_value = (*p_z_axis)[index_z];
-        double z_value = std::log(mp_intensity_output->currentValue());
+        double z_value = std::log(*it);
         p_hist2D->Fill(x_value, y_value, z_value);
-        ++index;
+        ++it;
     }
     p_hist2D->SetContour(50);
     gStyle->SetPalette(51);
diff --git a/App/src/TestFresnelCoeff.cpp b/App/src/TestFresnelCoeff.cpp
index 033228cd38ba6feba212d98ee810c6414d993eca..1956ccc097a3df8b46d21ca30c1bcf81c5bf7345 100644
--- a/App/src/TestFresnelCoeff.cpp
+++ b/App/src/TestFresnelCoeff.cpp
@@ -27,6 +27,8 @@
 
 
 TestFresnelCoeff::TestFresnelCoeff()
+: mp_sample(0)
+, mp_coeffs(0)
 {
     std::cout << "TestFresnelCoeff::TestFresnelCoeff() -> Info." << std::endl;
 }
@@ -62,29 +64,29 @@ void TestFresnelCoeff::test_standard_samples()
 
     // loop over standard samples defined in SampleFactory and StandardSamples
     for(size_t i_sample=0; i_sample<snames.size(); i_sample++){
-        m_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample(snames[i_sample]));
+        mp_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample(snames[i_sample]));
 
-        m_coeffs = new OutputData<OpticalFresnel::MultiLayerCoeff_t >;
-        m_coeffs->addAxis(std::string("alpha_i"), 0.0*Units::degree, 2.0*Units::degree, 2000);
-        m_coeffs->resetIndex();
-        while (m_coeffs->hasNext()) {
-            double alpha_i = m_coeffs->getCurrentValueOfAxis<double>("alpha_i");
+        mp_coeffs = new OutputData<OpticalFresnel::MultiLayerCoeff_t >;
+        mp_coeffs->addAxis(std::string("alpha_i"), 0.0*Units::degree, 2.0*Units::degree, 2000);
+        OutputData<OpticalFresnel::MultiLayerCoeff_t >::iterator it = mp_coeffs->begin();
+        while (it != mp_coeffs->end()) {
+            double alpha_i = mp_coeffs->getValueOfAxis<double>("alpha_i", it.getIndex());
             kvector_t kvec;
             kvec.setLambdaAlphaPhi(1.54*Units::angstrom, -alpha_i, 0.0);
 
             OpticalFresnel::MultiLayerCoeff_t coeffs;
             OpticalFresnel fresnelCalculator;
-            fresnelCalculator.execute(*m_sample, kvec, coeffs);
+            fresnelCalculator.execute(*mp_sample, kvec, coeffs);
 
-            m_coeffs->next() = coeffs;
+            *it = coeffs;
+            ++it;
 
         } // alpha_i
 
         draw_standard_samples();
 
-        delete m_sample;
-        delete m_coeffs;
-        //break;
+        delete mp_sample;
+        delete mp_coeffs;
     } // i_sample
 }
 
@@ -96,7 +98,7 @@ void TestFresnelCoeff::draw_standard_samples()
 {
     static int ncall = 0;
 
-    size_t nlayers = m_sample->getNumberOfLayers();
+    size_t nlayers = mp_sample->getNumberOfLayers();
 
     // graphics for R,T coefficients in layers as a function of alpha_i
     size_t ncoeffs = 2;
@@ -114,32 +116,11 @@ void TestFresnelCoeff::draw_standard_samples()
     }
     TGraph *gr_absSum = new TGraph(); // |R_top|+|T_bottom|
 
-    m_coeffs->resetIndex();
+    OutputData<OpticalFresnel::MultiLayerCoeff_t >::const_iterator it = mp_coeffs->begin();
     int i_point = 0;
-    while (m_coeffs->hasNext())
-    {
-        double alpha_i = m_coeffs->getCurrentValueOfAxis<double>("alpha_i");
-        OpticalFresnel::MultiLayerCoeff_t coeffs = m_coeffs->next();
-
-        // debug printing
-//        size_t index_alpha = i_point;
-//        if( index_alpha%100==0 ) {
-//            std::cout << "alpha_i: " << index_alpha << " " <<std::setprecision(20) << alpha_i << std::endl;
-//            for(size_t i_layer=0; i_layer<nlayers; ++i_layer ) {
-//                std::cout << std::setprecision(12) << " L:" << i_layer
-//                          << " kz:" << coeffs[i_layer].kz/10.
-//                          << " rt:"
-//                          << coeffs[i_layer].r
-//                          << coeffs[i_layer].t
-//                          << coeffs[i_layer].rb
-//                          << coeffs[i_layer].tb
-//                          << "X:" << coeffs[i_layer].X
-//                          << "R:" << coeffs[i_layer].R
-//                          << "T:" << coeffs[i_layer].T
-//                          << std::endl;
-//                //if(index_alpha%100==0) std::cout << " L:" << i_layer << " R:" << coeffs[i_layer].R << " T:" << coeffs[i_layer].T << std::endl;
-//            }
-//        }
+    while (it != mp_coeffs->end()) {
+        double alpha_i = mp_coeffs->getValueOfAxis<double>("alpha_i", it.getIndex());
+        const OpticalFresnel::MultiLayerCoeff_t coeffs = *it++;
 
         // Filling graphics for R,T as a function of alpha_i
         for(size_t i_layer=0; i_layer<nlayers; ++i_layer ) {
@@ -158,7 +139,6 @@ void TestFresnelCoeff::draw_standard_samples()
             std::cout << " and Re(k_{z,N+1}) = " << coeffs[nlast].kz.real() << std::endl;
         }
         gr_absSum->SetPoint(i_point++, Units::rad2deg(alpha_i), sum);
-
     }
 
     // create name of canvas different for each new call of this method
@@ -230,58 +210,47 @@ void TestFresnelCoeff::draw_standard_samples()
 
     // drawing sample geometry
     c1->cd(nlayers+1);
-    DrawHelper::instance().DrawMultilayer(m_sample);
+    DrawHelper::instance().DrawMultilayer(mp_sample);
 }
 
-
-
 /* ************************************************************************* */
 //! calculate fresnel coefficients .vs. alpha_i for set of roughnesses
 /* ************************************************************************* */
 void TestFresnelCoeff::test_roughness_set()
 {
-    m_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample("MultilayerOffspecTestcase1a"));
-//    m_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample("SubstrateOnSubstrate"));
+    mp_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample("MultilayerOffspecTestcase1a"));
 
-    std::cout << *m_sample << std::endl;
+    std::cout << *mp_sample << std::endl;
     std::cout << "-----" << std::endl;
-    ParameterPool *newpool = m_sample->createParameterTree();
+    ParameterPool *newpool = mp_sample->createParameterTree();
     std::cout << *newpool << std::endl;
 
-//    std::cout << "-----" << std::endl;
     FitMultiParameter multipar;
-//    double x = 9.0;
-//    FitMultiParameter::parameter_t p(&x);
-//    std::cout << p  << " " << &x << std::endl;
-//    multipar.addParameter(p);
 
-//    multipar.addMatchedParametersFromPool("/multilayer/interface0/roughness/sigma",newpool);
     multipar.addMatchedParametersFromPool("/*/*/*/sigma",newpool);
     std::cout << "multipar: " << multipar << std::endl;
-//    multipar.setValue(0.0);
 
-    std::cout << *m_sample << std::endl;
+    std::cout << *mp_sample << std::endl;
 
-    m_coeffs = new OutputData<OpticalFresnel::MultiLayerCoeff_t >;
-    m_coeffs->addAxis(std::string("alpha_i"), 0.0*Units::degree, 2.0*Units::degree, 1000);
-    m_coeffs->addAxis(std::string("roughness"), 0.0, 12.0*Units::nanometer, 6);
-    m_coeffs->resetIndex();
-    while (m_coeffs->hasNext()) {
-        double alpha_i = m_coeffs->getCurrentValueOfAxis<double>("alpha_i");
-        double roughness = m_coeffs->getCurrentValueOfAxis<double>("roughness");
+    mp_coeffs = new OutputData<OpticalFresnel::MultiLayerCoeff_t >;
+    mp_coeffs->addAxis(std::string("alpha_i"), 0.0*Units::degree, 2.0*Units::degree, 1000);
+    mp_coeffs->addAxis(std::string("roughness"), 0.0, 12.0*Units::nanometer, 6);
+    OutputData<OpticalFresnel::MultiLayerCoeff_t >::iterator it = mp_coeffs->begin();
+    while (it != mp_coeffs->end()) {
+        double alpha_i = mp_coeffs->getValueOfAxis<double>("alpha_i", it.getIndex());
+        double roughness = mp_coeffs->getValueOfAxis<double>("roughness", it.getIndex());
         multipar.setValue(roughness);
 
         kvector_t kvec;
         kvec.setLambdaAlphaPhi(1.54*Units::angstrom, -alpha_i, 0.0);
         OpticalFresnel::MultiLayerCoeff_t coeffs;
         OpticalFresnel fresnelCalculator;
-        fresnelCalculator.execute(*m_sample, kvec, coeffs);
-        m_coeffs->next() = coeffs;
-        //std::cout << alpha_i << " " << roughness << " " << coeffs[0].R << std::endl;
+        fresnelCalculator.execute(*mp_sample, kvec, coeffs);
+        *it = coeffs;
+        ++it;
     }
 
     draw_roughness_set();
-
 }
 
 
@@ -293,16 +262,13 @@ void TestFresnelCoeff::draw_roughness_set()
 {
     static int ncall = 0;
 
-    size_t nlayers = m_sample->getNumberOfLayers();
+    size_t nlayers = mp_sample->getNumberOfLayers();
 
     // graphics for R,T coefficients in layers as a function of alpha_i
     size_t ncoeffs = 2;
     enum key_coeffs { kCoeffR, kCoeffT};
 
-    m_coeffs->resetIndex();
-    MultiIndex& index = m_coeffs->getIndex();
-    //NamedVector<double> *p_alpha = dynamic_cast<NamedVector<double>*>(m_coeffs->getAxis("alpha_i"));
-    const NamedVector<double> *p_rough = dynamic_cast<const NamedVector<double>*>(m_coeffs->getAxis("roughness"));
+    const NamedVector<double> *p_rough = dynamic_cast<const NamedVector<double>*>(mp_coeffs->getAxis("roughness"));
     size_t nroughness = p_rough->getSize();
 
 
@@ -317,36 +283,21 @@ void TestFresnelCoeff::draw_roughness_set()
             }
         }
     }
-//    TGraph *gr_absSum = new TGraph(); // |R_top|+|T_bottom|
-
 
-    while (!index.endPassed())
-    {
-        double alpha_i = m_coeffs->getCurrentValueOfAxis<double>("alpha_i");
-        //double roughness = m_coeffs->getCurrentValueOfAxis<double>("roughness");
-        size_t i_alpha = index.getCurrentIndexOfAxis("alpha_i");
-        size_t i_rough = index.getCurrentIndexOfAxis("roughness");
+    OutputData<OpticalFresnel::MultiLayerCoeff_t >::const_iterator it = mp_coeffs->begin();
+    while (it != mp_coeffs->end()) {
+        double alpha_i = mp_coeffs->getValueOfAxis<double>("alpha_i", it.getIndex());
+        size_t i_alpha = mp_coeffs->getIndexOfAxis("alpha_i", it.getIndex());
+        size_t i_rough = mp_coeffs->getIndexOfAxis("roughness", it.getIndex());
 
-        OpticalFresnel::MultiLayerCoeff_t coeffs = m_coeffs->currentValue();
+        OpticalFresnel::MultiLayerCoeff_t coeffs = *it;
 
         // Filling graphics for R,T as a function of alpha_i
         for(size_t i_layer=0; i_layer<nlayers; ++i_layer ) {
             gr_coeff[i_layer][kCoeffR][i_rough]->SetPoint(i_alpha, Units::rad2deg(alpha_i), std::abs(coeffs[i_layer].R) );
             gr_coeff[i_layer][kCoeffT][i_rough]->SetPoint(i_alpha, Units::rad2deg(alpha_i), std::abs(coeffs[i_layer].T) );
         }
-
-//        // Filling graphics for |R|+|T| as a function of alpha_i taking R from the top and T from the bottom layers
-//        int nlast = nlayers - 1;
-//        double sum;
-//        if(coeffs[0].kz.real()!=0.0) {
-//            sum = std::norm(coeffs[0].R) + std::norm(coeffs[nlast].T)*coeffs[nlast].kz.real()/coeffs[0].kz.real();
-//        } else {
-//            sum = 1.0;
-//            std::cout << "Re(k_{z,0}) = 0 for alpha_i = " << alpha_i << std::endl;
-//            std::cout << " and Re(k_{z,N+1}) = " << coeffs[nlast].kz.real() << std::endl;
-//        }
-//        gr_absSum->SetPoint(i_point++, Units::rad2deg(alpha_i), sum);
-        ++index;
+        ++it;
     }
 
     // create name of canvas different for each new call of this method
@@ -385,7 +336,6 @@ void TestFresnelCoeff::draw_roughness_set()
         h1ref.GetYaxis()->SetTitle("|R|, |T|");
         h1ref.DrawCopy();
 
-//        TLegend *leg = new TLegend(0.725,0.7,0.89,0.88);
         TLegend *leg = new TLegend(0.18,0.20,0.28,0.69);
         leg->SetTextSize(0.04);
         leg->SetBorderSize(1);
@@ -407,21 +357,9 @@ void TestFresnelCoeff::draw_roughness_set()
         }
         leg->Draw();
     }
-////    TGraph *gr = gr_absSum;
-////    gr->SetMarkerStyle(21);
-////    gr->SetMarkerSize(0.2);
-////    gr->SetLineColor(kMagenta);
-////    gr->SetMarkerColor(kMagenta);
-////    gr->Draw("pl same");
-//    TLegend *leg = new TLegend(0.625,0.6,0.89,0.69);
-//    leg->SetTextSize(0.04);
-//    leg->SetBorderSize(0);
-//    leg->SetFillStyle(0);
-//    leg->AddEntry(gr, "|R_top|+|T_bottom|","pl");
-//    leg->Draw();
 
     // drawing sample geometry
     c1->cd(nlayers+1);
-    DrawHelper::instance().DrawMultilayer(m_sample);
+    DrawHelper::instance().DrawMultilayer(mp_sample);
 }
 
diff --git a/App/src/TestMiscellaneous.cpp b/App/src/TestMiscellaneous.cpp
index c67a27d951fbfbeb12421d8c21d200778bf2d463..c4c614e067aa2da534a00d915dae8fee374a7220 100644
--- a/App/src/TestMiscellaneous.cpp
+++ b/App/src/TestMiscellaneous.cpp
@@ -52,19 +52,20 @@ void TestMiscellaneous::test_OutputDataTo2DArray()
     // [axis0][axis1]
     int axis0_size = 2;
     int axis1_size = 4;
-    OutputData<double> *output = new OutputData<double>;
-    output->addAxis(std::string("axis0"), 0.0, double(axis0_size), axis0_size);
-    output->addAxis(std::string("axis1"), 0.0, double(axis1_size), axis1_size);
-    output->setAllTo(0.0);
+    OutputData<double> *p_output = new OutputData<double>;
+    p_output->addAxis(std::string("axis0"), 0.0, double(axis0_size), axis0_size);
+    p_output->addAxis(std::string("axis1"), 0.0, double(axis1_size), axis1_size);
+    p_output->setAllTo(0.0);
 
-    output->resetIndex();
+    OutputData<double>::iterator it = p_output->begin();
     int nn=0;
-    while (output->hasNext())
+    while (it != p_output->end())
     {
-        size_t index0 = output->getCurrentIndexOfAxis("x-axis");
-        size_t index1 = output->getCurrentIndexOfAxis("y-axis");
+        size_t index0 = p_output->getIndexOfAxis("x-axis", it.getIndex());
+        size_t index1 = p_output->getIndexOfAxis("y-axis", it.getIndex());
         std::cout << " index0:" << index0 << " index1:" << index1 << std::endl;
-        output->next() = nn++;
+        *it = nn++;
+        ++it;
     }
 
 
@@ -217,28 +218,24 @@ void TestMiscellaneous::test_FormFactor()
     }
 
 
-    OutputData<double > *data1 = new OutputData<double >;
-    data1->addAxis(std::string("qx"), qmin, qmax, nbins);
-    data1->addAxis(std::string("qy"), qmin, qmax, nbins);
-    data1->addAxis(std::string("qz"), qmin, qmax, nbins);
-    data1->resetIndex();
-    while (data1->hasNext()) {
-        double x = data1->getCurrentValueOfAxis<double>("qx");
-        double y = data1->getCurrentValueOfAxis<double>("qy");
-        double z = data1->getCurrentValueOfAxis<double>("qz");
+    OutputData<double> *p_data = new OutputData<double>();
+    p_data->addAxis(std::string("qx"), qmin, qmax, nbins);
+    p_data->addAxis(std::string("qy"), qmin, qmax, nbins);
+    p_data->addAxis(std::string("qz"), qmin, qmax, nbins);
+    OutputData<double>::const_iterator it = p_data->begin();
+    while (it != p_data->end()) {
+        double x = p_data->getValueOfAxis<double>("qx", it.getIndex());
+        double y = p_data->getValueOfAxis<double>("qy", it.getIndex());
+        double z = p_data->getValueOfAxis<double>("qz", it.getIndex());
 
-        int ix = data1->getCurrentIndexOfAxis("qx");
-        int iy = data1->getCurrentIndexOfAxis("qy");
-        int iz = data1->getCurrentIndexOfAxis("qz");
+        int ix = p_data->getIndexOfAxis("qx", it.getIndex());
+        int iy = p_data->getIndexOfAxis("qy", it.getIndex());
+        int iz = p_data->getIndexOfAxis("qz", it.getIndex());
 
         cvector_t q(x,y,z);
         cvector_t q0(0,0,0);
         double value = std::abs(ff.evaluate(q,q0, 0.0, 0.0));
-//        complex_t qz(z,0.0);
-//        std::cout << q << " " << std::abs(ff.evaluate(q,q0)) << std::endl;
         if(iz==50) h2->Fill(x,y,std::abs(ff.evaluate(q,q0, 0.0, 0.0)));
-        //if(iy==0) h2->Fill(x,z,std::abs(ff.evaluate(q,q0)));
-        //if(ix==0) h2->Fill(z,y,std::abs(ff.evaluate(q,q0)));
 
         h3->Fill(x,y,z,std::abs(ff.evaluate(q,q0, 0.0, 0.0)));
 
@@ -255,7 +252,7 @@ void TestMiscellaneous::test_FormFactor()
         vh2_xz[iy] ->Fill(x,z,value);
         vh2_yz[ix] ->Fill(y,z,value);
 
-        data1->next();
+        ++it;
     }
 
     TCanvas *c1_xy = new TCanvas("c1_xy","c1_xy",1024,768);
diff --git a/App/src/TestMultiLayerRoughness.cpp b/App/src/TestMultiLayerRoughness.cpp
index e02ff2b5601ee22679c47915287c4b451df14c09..81f58393dde3d049538cbd44378c1abb7e491f21 100644
--- a/App/src/TestMultiLayerRoughness.cpp
+++ b/App/src/TestMultiLayerRoughness.cpp
@@ -33,12 +33,12 @@ void TestMultiLayerRoughness::execute()
     TH2D *h2 = new TH2D ("h2","h2", npoints, alphaMin-dalpha/2., alphaMax+dalpha/2., npoints, alphaMin-dalpha/2., alphaMax+dalpha/2. );
     h2->SetContour(50);
 
-    OutputData<double > *data_alpha_i = new OutputData<double >;
-    data_alpha_i->addAxis(std::string("alpha_i"), 0.0*Units::degree, 2.0*Units::degree, npoints);
-    data_alpha_i->resetIndex();
-    while (data_alpha_i->hasNext()) {
-        double alpha_i = data_alpha_i->getCurrentValueOfAxis<double>("alpha_i");
-        size_t index_alpha_i = data_alpha_i->getCurrentIndexOfAxis("alpha_i");
+    OutputData<double> *p_data_alpha_i = new OutputData<double>();
+    p_data_alpha_i->addAxis(std::string("alpha_i"), 0.0*Units::degree, 2.0*Units::degree, npoints);
+    OutputData<double>::const_iterator it_alpha_i = p_data_alpha_i->begin();
+    while (it_alpha_i != p_data_alpha_i->end()) {
+        double alpha_i = p_data_alpha_i->getValueOfAxis<double>("alpha_i", it_alpha_i.getIndex());
+        size_t index_alpha_i = p_data_alpha_i->getIndexOfAxis("alpha_i", it_alpha_i.getIndex());
         if(index_alpha_i%10 == 0) std::cout << index_alpha_i << " of " << npoints << std::endl;
 
         // setting experiment
@@ -50,19 +50,16 @@ void TestMultiLayerRoughness::execute()
         experiment.runSimulation();
 
         const OutputData<double> *output = experiment.getOutputData();
-        output->resetIndex();
-        while (output->hasNext()) {
-            double phi_f = output->getCurrentValueOfAxis<double>("phi_f");
-            double alpha_f = output->getCurrentValueOfAxis<double>("alpha_f");
-            //size_t index_phi = output->getCurrentIndexOfAxis("phi_f");
-            double intensity = output->next();
-            //std::cout << "alpha_i " << alpha_i << " alpha_f " << alpha_f << " phi_f " << phi_f << " inten " << intensity << std::endl;
+        OutputData<double>::const_iterator it_output = output->begin();
+        while (it_output != output->end()) {
+            double phi_f = output->getValueOfAxis<double>("phi_f", it_output.getIndex());
+            double alpha_f = output->getValueOfAxis<double>("alpha_f", it_output.getIndex());
+            double intensity = *it_output++;
             if(phi_f == 0) {
                 h2->Fill(Units::rad2deg(alpha_i), Units::rad2deg(alpha_f), intensity);
             }
         }
-
-        data_alpha_i->next();
+        ++it_alpha_i;
     }
 
     gStyle->SetPalette(1);
diff --git a/App/src/TestRootTree.cpp b/App/src/TestRootTree.cpp
index e8c7e3333c717a2a48346186fd3bd2977eca4275..4d20ca4190247331cd1d527614cc009ce5ed67cd 100644
--- a/App/src/TestRootTree.cpp
+++ b/App/src/TestRootTree.cpp
@@ -21,7 +21,7 @@
 
 #include <vector>
 
-TestRootTree::TestRootTree() : m_sample(0), m_experiment(0), m_data(0)
+TestRootTree::TestRootTree() : mp_sample(0), mp_experiment(0), mp_data(0)
 {
 
 }
@@ -29,9 +29,9 @@ TestRootTree::TestRootTree() : m_sample(0), m_experiment(0), m_data(0)
 
 TestRootTree::~TestRootTree()
 {
-    delete m_sample;
-    delete m_experiment;
-    delete m_data;
+    delete mp_sample;
+    delete mp_experiment;
+    delete mp_data;
 }
 
 
@@ -107,7 +107,7 @@ void TestRootTree::complex_write()
         double alpha_f_min(-0.4*Units::degree), alpha_f_max(0.066);
 
         GISASExperiment experiment(mp_options);
-        experiment.setSample(*m_sample);
+        experiment.setSample(*mp_sample);
         experiment.setDetectorParameters(nphi_f, phi_f_min, phi_f_max, nalpha_f , alpha_f_min, alpha_f_max);
         experiment.setBeamParameters(1.77*Units::angstrom, -alpha_i, phi_i);
         experiment.setBeamIntensity(1e7);
@@ -128,9 +128,9 @@ void TestRootTree::complex_write()
         event->mphi = Units::rad2deg(meso_phi);
         event->malpha = Units::rad2deg(meso_alpha);
         // copying output data into event frame
-        delete m_data;
-        m_data = experiment.getOutputDataClone();
-        IsGISAXSTools::exportOutputDataInVectors2D(*m_data, event->vi, event->vphi_f, event->valpha_f);
+        delete mp_data;
+        mp_data = experiment.getOutputDataClone();
+        IsGISAXSTools::exportOutputDataInVectors2D(*mp_data, event->vi, event->vphi_f, event->valpha_f);
 
         // lets switch to degtrees
         for(size_t i=0; i<event->vphi_f.size(); i++){
@@ -144,7 +144,7 @@ void TestRootTree::complex_write()
         c1->cd(); gPad->SetLogz();
         c1->Clear();
         IsGISAXSTools::setMinimum(1.);
-        IsGISAXSTools::drawOutputDataInPad(*m_data, "CONT4 Z", "IsGisaxs pyramid FF");
+        IsGISAXSTools::drawOutputDataInPad(*mp_data, "CONT4 Z", "IsGisaxs pyramid FF");
         c1->Modified();
         c1->Update();
 
@@ -202,17 +202,17 @@ void TestRootTree::complex_read()
 /* ************************************************************************* */
 void TestRootTree::simple_write()
 {
-    delete m_experiment;
-    delete m_sample;
-    delete m_data;
+    delete mp_experiment;
+    delete mp_sample;
+    delete mp_data;
 
-    m_sample = dynamic_cast<MultiLayer *>(SampleFactory::instance().createItem("IsGISAXS9_Pyramid"));
+    mp_sample = dynamic_cast<MultiLayer *>(SampleFactory::instance().createItem("IsGISAXS9_Pyramid"));
 
     // setting experiment
-    m_experiment = new GISASExperiment(mp_options);
-    m_experiment->setDetectorParameters(100, 0.0*Units::degree, 2.0*Units::degree, 100, 0.0*Units::degree, 2.0*Units::degree, true);
-    m_experiment->setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
-    m_experiment->setSample(*m_sample);
+    mp_experiment = new GISASExperiment(mp_options);
+    mp_experiment->setDetectorParameters(100, 0.0*Units::degree, 2.0*Units::degree, 100, 0.0*Units::degree, 2.0*Units::degree, true);
+    mp_experiment->setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
+    mp_experiment->setSample(*mp_sample);
 
     // variables below will be written in the tree
     double intens1(0), intens2(0), alpha_i(0), phi_i(0), alpha_f(0), phi_f(0);
@@ -254,38 +254,36 @@ void TestRootTree::simple_write()
         phi_i = 0;
         nev = i_ev;
 
-        m_experiment->setBeamParameters(1.0*Units::angstrom, alpha_i*Units::degree, phi_i);
-        m_experiment->runSimulation();
+        mp_experiment->setBeamParameters(1.0*Units::angstrom, alpha_i*Units::degree, phi_i);
+        mp_experiment->runSimulation();
 
-        m_data = m_experiment->getOutputDataClone();
+        mp_data = mp_experiment->getOutputDataClone();
         // accessing to scattering data
-        NamedVector<double> *axis0 = dynamic_cast<NamedVector<double>*>(m_data->getAxes()[0]);
-        NamedVector<double> *axis1 = dynamic_cast<NamedVector<double>*>(m_data->getAxes()[1]);
+        NamedVector<double> *axis0 = dynamic_cast<NamedVector<double>*>(mp_data->getAxes()[0]);
+        NamedVector<double> *axis1 = dynamic_cast<NamedVector<double>*>(mp_data->getAxes()[1]);
         std::string axis0_name = axis0->getName();
         std::string axis1_name = axis1->getName();
 
         c1->cd(); gPad->SetLogz();
         c1->Clear();
         IsGISAXSTools::setMinimum(1.);
-        IsGISAXSTools::drawOutputDataInPad(*m_data, "CONT4 Z", "IsGisaxs pyramid FF");
+        IsGISAXSTools::drawOutputDataInPad(*mp_data, "CONT4 Z", "IsGisaxs pyramid FF");
         c1->Modified();
         c1->Update();
         c1->Write();
 
-        m_data->resetIndex();
-        while (m_data->hasNext())
+        OutputData<double>::const_iterator it = mp_data->begin();
+        while (it != mp_data->end())
         {
-            size_t index_phi_f =  m_data->getCurrentIndexOfAxis(axis0_name.c_str());
-            size_t index_alpha_f = m_data->getCurrentIndexOfAxis(axis1_name.c_str());
+            size_t index_phi_f =  mp_data->getIndexOfAxis(axis0_name.c_str(), it.getIndex());
+            size_t index_alpha_f = mp_data->getIndexOfAxis(axis1_name.c_str(), it.getIndex());
             phi_f = Units::rad2deg( (*axis0)[index_phi_f]);
             alpha_f = Units::rad2deg( (*axis1)[index_alpha_f] );
             //std::cout << phi_f << " " << alpha_f << std::endl;
-            intens1 =m_data->next();
+            intens1 = *it++;
             tree->Fill();
-
         }
-
-        delete m_data;
+        delete mp_data;
     }
 
     tree->Write();
@@ -367,7 +365,7 @@ void TestRootTree::simple_read()
 void TestRootTree::initializeMesoCrystal(double meso_alpha, double meso_phi, double nanopart_radius)
 {
     (void)nanopart_radius;
-    delete m_sample;
+    delete mp_sample;
     // create mesocrystal
     double meso_radius = 300*Units::nanometer;
     double surface_filling_ratio = 0.25;
@@ -409,6 +407,6 @@ void TestRootTree::initializeMesoCrystal(double meso_alpha, double meso_phi, dou
     p_multi_layer->addLayer(air_layer);
     p_multi_layer->addLayer(avg_layer_decorator);
     p_multi_layer->addLayer(substrate_layer);
-    m_sample = p_multi_layer;
+    mp_sample = p_multi_layer;
 }
 
diff --git a/Core/Algorithms/inc/DWBASimulation.h b/Core/Algorithms/inc/DWBASimulation.h
index 1e78adc9a88c798a334c3d93976e563faa02e605..48d81a5593abc83e3fdea1872820d30ac1915fcd 100644
--- a/Core/Algorithms/inc/DWBASimulation.h
+++ b/Core/Algorithms/inc/DWBASimulation.h
@@ -41,12 +41,6 @@ public:
     virtual DWBASimulation *clone();
 
 protected:
-
-    virtual void resetIndex();
-    virtual bool hasNext() const;
-    virtual const double &next() const;
-    virtual double &next();
-
     OutputData<double> m_dwba_intensity;
     OutputData<double> m_output_data_mask;
     cvector_t m_ki;
diff --git a/Core/Algorithms/inc/GISASExperiment.h b/Core/Algorithms/inc/GISASExperiment.h
index 51074af1756af23ad1b3bfa932542efe0892ec0d..498f6ff81075f9b20453cc3cf312cdcdf638e9f3 100644
--- a/Core/Algorithms/inc/GISASExperiment.h
+++ b/Core/Algorithms/inc/GISASExperiment.h
@@ -44,12 +44,12 @@ private:
     virtual void init_parameters();
 
 	void initializeAnglesIsgisaxs(NamedVector<double> *p_axis, double start, double end, size_t size);
-	double getCurrentSolidAngle() const;
+	double getSolidAngle(size_t index) const;
 	double deltaAlpha(double alpha, double zeta) const;
 	double deltaPhi(double alpha, double phi, double zeta);
 	void createZetaAndProbVectors(std::vector<double> &zetas, std::vector<double> &probs, size_t nbr_zetas, double zeta_sigma);
 	void addToIntensityMap(double alpha, double phi, double value);
-	size_t findClosestIndex(const NamedVector<double> *p_axis, double value);
+	int findClosestIndex(const NamedVector<double> *p_axis, double value);
 };
 
 #endif /* GISASEXPERIMENT_H_ */
diff --git a/Core/Algorithms/src/ChiSquaredFrequency.cpp b/Core/Algorithms/src/ChiSquaredFrequency.cpp
index 838f76211f3e4987b18d4a43eac6a4c25fe20bb5..b108684f354c08931aa7709a69c13214014cc469 100644
--- a/Core/Algorithms/src/ChiSquaredFrequency.cpp
+++ b/Core/Algorithms/src/ChiSquaredFrequency.cpp
@@ -25,10 +25,10 @@ double ChiSquaredFrequency::calculateChiSquared(
     initWeights();
     size_t data_size = mp_weights->getAllocatedSize();
     OutputData<double> *p_difference = createChi2DifferenceMap();
-    mp_weights->resetIndex();
-    p_difference->resetIndex();
-    while(p_difference->hasNext()) {
-        result += p_difference->next()*mp_weights->next();
+    OutputData<double>::const_iterator it_weights = mp_weights->begin();
+    OutputData<double>::const_iterator it_diff = p_difference->begin();
+    while(it_diff != p_difference->end()) {
+        result += (*it_diff++)*(*it_weights++);
     }
     delete p_difference;
     m_chi2_value = result/data_size;
@@ -38,17 +38,18 @@ OutputData<double>* ChiSquaredFrequency::createChi2DifferenceMap() const
 {
     OutputData<double> *p_difference = mp_weights->clone();
 
-    p_difference->resetIndex();
-    mp_simulation_ft->resetIndex();
-    mp_real_ft->resetIndex();
+    OutputData<double>::iterator it_diff = p_difference->begin();
+    OutputData<complex_t>::const_iterator it_sim = mp_simulation_ft->begin();
+    OutputData<complex_t>::const_iterator it_real = mp_real_ft->begin();
 
-    while (p_difference->hasNext()) {
-        complex_t real = mp_real_ft->next();
-        complex_t simu = mp_simulation_ft->next();
+    while (it_diff != p_difference->end()) {
+        complex_t real = *it_real++;
+        complex_t simu = *it_sim++;
         complex_t diff = real - simu;
 //        double sum_norms = std::max(std::norm(real), 1.0);
         double squared_difference = std::norm(diff); // /std::sqrt(sum_norms);
-        p_difference->next() = squared_difference;
+        *it_diff = squared_difference;
+        ++it_diff;
     }
 
     return p_difference;
@@ -74,21 +75,21 @@ void ChiSquaredFrequency::initWeights()
     }
     mp_weights->setAxisSizes(rank, n_dims);
     delete n_dims;
-    mp_weights->resetIndex();
+    OutputData<double>::iterator it_weights = mp_weights->begin();
     size_t nbr_rows = mp_weights->getAllSizes()[0];
     size_t row_length = mp_weights->getAllSizes()[1];
     size_t row_number = 0;
-    size_t counter = 0;
-    while (mp_weights->hasNext()) {
+    while (it_weights != mp_weights->end()) {
         double weight = 0.0;
-        size_t column_index = counter%row_length;
+        size_t linear_index = it_weights.getIndex();
+        size_t column_index = linear_index%row_length;
         int shift = (column_index>=row_length/2 ? -(int)row_length : 0);
         int centered_column_index = (int)column_index - shift;
         if (row_number<m_cutoff*nbr_rows && centered_column_index < m_cutoff*row_length/2.0) {
             weight = 1.0;
         }
-        mp_weights->next() = weight;
-        ++counter;
-        if (counter%row_length==0) ++row_number;
+        *it_weights = weight;
+        if (linear_index%row_length==0) ++row_number;
+        ++it_weights;
     }
 }
diff --git a/Core/Algorithms/src/ChiSquaredModule.cpp b/Core/Algorithms/src/ChiSquaredModule.cpp
index 4ca4656d94c6fdbfd7222dea2880f21407c1596a..c7f54863ad15b6e7c43fd3081b458e8e33435a6f 100644
--- a/Core/Algorithms/src/ChiSquaredModule.cpp
+++ b/Core/Algorithms/src/ChiSquaredModule.cpp
@@ -22,10 +22,10 @@ double ChiSquaredModule::calculateChiSquared(
     size_t data_size = mp_real_data->getAllocatedSize();
     initWeights();
     OutputData<double> *p_difference = createChi2DifferenceMap();
-    mp_weights->resetIndex();
-    p_difference->resetIndex();
-    while(p_difference->hasNext()) {
-        result += p_difference->next()*mp_weights->next();
+    OutputData<double>::const_iterator it_weights = mp_weights->begin();
+    OutputData<double>::const_iterator it_diff = p_difference->begin();
+    while(it_diff != p_difference->end()) {
+        result += (*it_diff++)*(*it_weights++);
     }
     delete p_difference;
     m_chi2_value = result/data_size;
@@ -37,14 +37,16 @@ OutputData<double>* ChiSquaredModule::createChi2DifferenceMap() const
     OutputData<double > *p_difference = mp_simulation_data->clone();
     p_difference->setAllTo(0.0);
 
-    mp_simulation_data->resetIndex();
-    mp_real_data->resetIndex();
-    p_difference->resetIndex();
-    while (mp_real_data->hasNext()) {
-        double value_simu = mp_simulation_data->next();
-        double value_real = mp_real_data->next();
+    OutputData<double>::iterator it_diff = p_difference->begin();
+    OutputData<double>::const_iterator it_sim = mp_simulation_data->begin();
+    OutputData<double>::const_iterator it_real = mp_real_data->begin();
+
+    while (it_diff != p_difference->end()) {
+        double value_simu = *it_sim++;
+        double value_real = *it_real++;
         double squared_difference = mp_squared_function->calculateSquaredDifference(value_real, value_simu);
-        p_difference->next() = squared_difference;
+        *it_diff = squared_difference;
+        ++it_diff;
     }
 
     return p_difference;
diff --git a/Core/Algorithms/src/DWBASimulation.cpp b/Core/Algorithms/src/DWBASimulation.cpp
index f4f9cb56297e71e1f420d88857588b20b0d8f455..8db2845e7f28a36c4ee18f3d612cbef8800972ea 100644
--- a/Core/Algorithms/src/DWBASimulation.cpp
+++ b/Core/Algorithms/src/DWBASimulation.cpp
@@ -26,8 +26,6 @@ void DWBASimulation::init(const Experiment& experiment)
     m_output_data_mask.copyFrom(*experiment.getOutputDataMask());
 }
 
-
-
 DWBASimulation *DWBASimulation::clone()
 {
     DWBASimulation *sim = new DWBASimulation();
@@ -44,30 +42,3 @@ double DWBASimulation::getWaveLength() const
     kvector_t real_ki(m_ki.x().real(), m_ki.y().real(), m_ki.z().real());
     return 2.0*M_PI/real_ki.mag();
 }
-
-
-void DWBASimulation::resetIndex()
-{
-    m_output_data_mask.resetIndex();
-    m_dwba_intensity.resetIndex();
-}
-
-
-bool DWBASimulation::hasNext() const
-{
-    return m_dwba_intensity.hasNext();
-}
-
-
-const double &DWBASimulation::next() const
-{
-    m_output_data_mask.next();
-    return m_dwba_intensity.next();
-}
-
-double &DWBASimulation::next()
-{
-    m_output_data_mask.next();
-    return m_dwba_intensity.next();
-}
-
diff --git a/Core/Algorithms/src/DiffuseDWBASimulation.cpp b/Core/Algorithms/src/DiffuseDWBASimulation.cpp
index 50b6433b0ad3af29922895ae0d42b5e7b3348a76..3c272dab54f0d4c232d8bff6597b7c6fd4966d0c 100644
--- a/Core/Algorithms/src/DiffuseDWBASimulation.cpp
+++ b/Core/Algorithms/src/DiffuseDWBASimulation.cpp
@@ -22,17 +22,18 @@ void DiffuseDWBASimulation::run()
     initDiffuseFormFactorTerms(diffuse_terms, nbr_heights, samples_per_particle);
     double wavevector_scattering_factor = M_PI/getWaveLength()/getWaveLength();
 
-    resetIndex();
-    while ( hasNext() ) {
-        if( !m_output_data_mask.currentValue() ) {
-            next();
+    OutputData<double>::const_iterator it_mask = m_output_data_mask.begin();
+    OutputData<double>::iterator it_intensity = m_dwba_intensity.begin();
+    while ( it_intensity != m_dwba_intensity.end() ) {
+        if( !(*it_mask) ) {
+            ++it_mask, ++it_intensity;
             continue;
         }
 
-        double phi_f = getDWBAIntensity().getCurrentValueOfAxis<double>("phi_f");
-        double alpha_f = getDWBAIntensity().getCurrentValueOfAxis<double>("alpha_f");
+        double phi_f = getDWBAIntensity().getValueOfAxis<double>("phi_f", it_intensity.getIndex());
+        double alpha_f = getDWBAIntensity().getValueOfAxis<double>("alpha_f", it_intensity.getIndex());
         if (alpha_f<0) {
-            next();
+            ++it_mask, ++it_intensity;
             continue;
         }
         cvector_t k_f;
@@ -51,7 +52,8 @@ void DiffuseDWBASimulation::run()
             }
             total_intensity += p_diffuse_term->m_factor*(intensity - std::norm(amplitude));
         }
-        next() = total_intensity*wavevector_scattering_factor*wavevector_scattering_factor;
+        *it_intensity = total_intensity*wavevector_scattering_factor*wavevector_scattering_factor;
+        ++it_intensity, ++it_mask;
     }
 
     for (size_t i=0; i<diffuse_terms.size(); ++i) delete diffuse_terms[i];
diff --git a/Core/Algorithms/src/Experiment.cpp b/Core/Algorithms/src/Experiment.cpp
index 713319a705558fcf7c67521f32f4ee14d5f7ffc7..98a5dcf096e65f4cda6c4a9876293a96cc93175e 100644
--- a/Core/Algorithms/src/Experiment.cpp
+++ b/Core/Algorithms/src/Experiment.cpp
@@ -48,10 +48,10 @@ void Experiment::normalize()
 {
     double incident_intensity = m_beam.getIntensity();
     if (!m_is_normalized && incident_intensity!=1.0) {
-        m_intensity_map.resetIndex();
-        while (m_intensity_map.hasNext()) {
-            double old_value = m_intensity_map.currentValue();
-            m_intensity_map.next() = incident_intensity*old_value;
+        OutputData<double>::iterator it = m_intensity_map.begin();
+        while (it != m_intensity_map.end()) {
+            *it *= incident_intensity;
+            ++it;
         }
         m_is_normalized = true;
     }
@@ -110,13 +110,14 @@ void Experiment::setOutputDataMask(size_t n_chunks_total, size_t n_chunk )
     // copying topology from intensity data
     m_current_output_data_mask.copyFrom(m_intensity_map);
     // setting mask
-    m_current_output_data_mask.resetIndex();
-    while(m_current_output_data_mask.hasNext()) {
-        if(m_current_output_data_mask.getIndex().getPosition() % n_chunks_total == n_chunk) {
-            m_current_output_data_mask.next() = 1.0;
+    OutputData<double>::iterator it_mask = m_current_output_data_mask.begin();
+    while(it_mask != m_current_output_data_mask.end()) {
+        if(it_mask.getIndex() % n_chunks_total == n_chunk) {
+            *it_mask = 1.0;
         } else {
-            m_current_output_data_mask.next() = 0.0;
+            *it_mask = 0.0;
         }
+        ++it_mask;
     }
 
 }
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index f5040a5e9f57e80fe773d86a11d0468b7eb56271..ae150a4befd95b4ee92524cf30b96817702f1f79 100644
--- a/Core/Algorithms/src/GISASExperiment.cpp
+++ b/Core/Algorithms/src/GISASExperiment.cpp
@@ -95,7 +95,7 @@ void GISASExperiment::normalize()
         double incident_intensity = m_beam.getIntensity(); // Actually, this is the total number of neutrons hitting the sample
         double sin_alpha_i = std::abs(m_beam.getCentralK().cosTheta());
         for (OutputData<double>::iterator it = m_intensity_map.begin(); it != m_intensity_map.end(); ++it) {
-            double factor = incident_intensity*getCurrentSolidAngle()/sin_alpha_i;
+            double factor = incident_intensity*getSolidAngle(it.getIndex())/sin_alpha_i;
             (*it) *= factor;
         }
         m_is_normalized = true;
@@ -135,19 +135,18 @@ void GISASExperiment::smearIntensityFromZAxisTilting()
     std::vector<double> probs;
     createZetaAndProbVectors(zetas, probs, nbr_zetas, zeta_sigma);
 
-    OutputData<double> *p_old = m_intensity_map.clone();
+    OutputData<double> *p_clone = m_intensity_map.clone();
     m_intensity_map.setAllTo(0.0);
-    p_old->resetIndex();
-    while (p_old->hasNext()) {
-        double old_phi = p_old->getCurrentValueOfAxis<double>("phi_f");
-        double old_alpha = p_old->getCurrentValueOfAxis<double>("alpha_f");
+    OutputData<double>::const_iterator it_clone = p_clone->begin();
+    while (it_clone != p_clone->end()) {
+        double old_phi = p_clone->getValueOfAxis<double>("phi_f", it_clone.getIndex());
+        double old_alpha = p_clone->getValueOfAxis<double>("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]);
             double prob = probs[zeta_index];
-            addToIntensityMap(newalpha, newphi, prob*p_old->currentValue());
+            addToIntensityMap(newalpha, newphi, prob*(*it_clone++));
         }
-        p_old->next();
     }
 }
 
@@ -165,22 +164,22 @@ void GISASExperiment::initializeAnglesIsgisaxs(NamedVector<double> *p_axis, doub
     return;
 }
 
-double GISASExperiment::getCurrentSolidAngle() const
+double GISASExperiment::getSolidAngle(size_t index) const
 {
     const std::string s_alpha_f("alpha_f");
     const std::string s_phi_f("phi_f");
 
     const NamedVector<double> *p_alpha_axis = dynamic_cast<const NamedVector<double>* >(m_intensity_map.getAxis(s_alpha_f));
     const NamedVector<double> *p_phi_axis = dynamic_cast<const NamedVector<double>* >(m_intensity_map.getAxis(s_phi_f));
-    size_t alpha_index = m_intensity_map.getCurrentIndexOfAxis(s_alpha_f);
+    size_t alpha_index = m_intensity_map.getIndexOfAxis(s_alpha_f, index);
     size_t alpha_size = p_alpha_axis->getSize();
-    size_t phi_index = m_intensity_map.getCurrentIndexOfAxis(s_phi_f);
+    size_t phi_index = m_intensity_map.getIndexOfAxis(s_phi_f, index);
     size_t phi_size = p_phi_axis->getSize();
     if (alpha_size<2 || phi_size<2) {
         // Cannot determine detector cell size!
         return 0.0;
     }
-    double alpha_f = m_intensity_map.getCurrentValueOfAxis<double>(s_alpha_f);
+    double alpha_f = m_intensity_map.getValueOfAxis<double>(s_alpha_f, index);
     double cos_alpha_f = std::cos(alpha_f);
     double dalpha, dphi;
     if (alpha_index==0) {
@@ -243,19 +242,18 @@ void GISASExperiment::createZetaAndProbVectors(std::vector<double>& zetas,
 
 void GISASExperiment::addToIntensityMap(double alpha, double phi, double value)
 {
+    OutputData<double>::iterator it = m_intensity_map.begin();
     const NamedVector<double> *p_alpha_axis = dynamic_cast<const NamedVector<double> *>(m_intensity_map.getAxis("alpha_f"));
     const NamedVector<double> *p_phi_axis = dynamic_cast<const NamedVector<double> *>(m_intensity_map.getAxis("phi_f"));
-    MultiIndex &index = m_intensity_map.getIndex();
-    size_t i_alpha = findClosestIndex(p_alpha_axis, alpha);
-    size_t i_phi = findClosestIndex(p_phi_axis, phi);
-    index.setIndexOfAxis("alpha_f", i_alpha);
-    index.setIndexOfAxis("phi_f", i_phi);
-    m_intensity_map.currentValue() += value;
+    std::vector<int> coordinates;
+    coordinates.push_back(findClosestIndex(p_alpha_axis, alpha));
+    coordinates.push_back(findClosestIndex(p_phi_axis, phi));
+    it[m_intensity_map.toIndex(coordinates)] += value;
 }
 
-size_t GISASExperiment::findClosestIndex(const NamedVector<double> *p_axis, double value)
+int GISASExperiment::findClosestIndex(const NamedVector<double> *p_axis, double value)
 {
-    size_t result = 0;
+    int result = 0;
     double smallest_diff = std::abs(value-(*p_axis)[0]);
     for (size_t i=1; i<p_axis->getSize(); ++i) {
         double new_diff = std::abs(value-(*p_axis)[i]);
diff --git a/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp b/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
index b56872f253c81bfd090c4b615b16ca7210740503..c9c451500b939235f60f84a15f2df2b142923188 100644
--- a/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
+++ b/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
@@ -85,23 +85,26 @@ void LayerDecoratorDWBASimulation::calculateCoherentIntensity(IInterferenceFunct
     //std::cout << "Calculating coherent scattering..." << std::endl;
     double wavelength = getWaveLength();
     double total_surface_density = mp_layer_decorator->getTotalParticleSurfaceDensity();
-    resetIndex();
-    while ( hasNext() )
+
+    OutputData<double>::const_iterator it_mask = m_output_data_mask.begin();
+    OutputData<double>::iterator it_intensity = m_dwba_intensity.begin();
+    while ( it_intensity != m_dwba_intensity.end() )
     {
-        if( !m_output_data_mask.currentValue() ) {
-            next();
+        if( !(*it_mask) ) {
+            ++it_mask, ++it_intensity;
             continue;
         }
-        double phi_f = getDWBAIntensity().getCurrentValueOfAxis<double>("phi_f");
-        double alpha_f = getDWBAIntensity().getCurrentValueOfAxis<double>("alpha_f");
+        double phi_f = getDWBAIntensity().getValueOfAxis<double>("phi_f", it_intensity.getIndex());
+        double alpha_f = getDWBAIntensity().getValueOfAxis<double>("alpha_f", it_intensity.getIndex());
         if (alpha_f<0) {
-            next();
+            ++it_mask, ++it_intensity;
             continue;
         }
         cvector_t k_f;
         k_f.setLambdaAlphaPhi(wavelength, alpha_f, phi_f);
         k_f.setZ(mp_kz_function->evaluate(alpha_f));
-        next() = p_strategy->evaluate(m_ki, k_f, -m_alpha_i, alpha_f)*total_surface_density;
+        *it_intensity = p_strategy->evaluate(m_ki, k_f, -m_alpha_i, alpha_f)*total_surface_density;
+        ++it_mask, ++it_intensity;
     }
 }
 
diff --git a/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp b/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
index a0a7992da9f2c9fce538db0ade9fd42c0ed62b1f..a242746a9d168f060f80c5623f937261e91c8598 100644
--- a/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
+++ b/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
@@ -45,19 +45,22 @@ void MultiLayerRoughnessDWBASimulation::run()
 {
     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();
-    resetIndex();
-    while ( hasNext() )
+
+    OutputData<double>::const_iterator it_mask = m_output_data_mask.begin();
+    OutputData<double>::iterator it_intensity = m_dwba_intensity.begin();
+    while ( it_intensity != m_dwba_intensity.end() )
     {
-        if( !m_output_data_mask.currentValue() ) {
-            next();
+        if( !(*it_mask) ) {
+            ++it_mask, ++it_intensity;
             continue;
         }
 
-        double phi_f = getDWBAIntensity().getCurrentValueOfAxis<double>("phi_f");
-        double alpha_f = getDWBAIntensity().getCurrentValueOfAxis<double>("alpha_f");
+        double phi_f = getDWBAIntensity().getValueOfAxis<double>("phi_f", it_intensity.getIndex());
+        double alpha_f = getDWBAIntensity().getValueOfAxis<double>("alpha_f", it_intensity.getIndex());
         cvector_t k_f;
         k_f.setLambdaAlphaPhi(lambda, alpha_f, phi_f);
-        next() = evaluate(m_ki, k_f, -m_alpha_i, alpha_f);
+        *it_intensity = evaluate(m_ki, k_f, -m_alpha_i, alpha_f);
+        ++it_mask, ++it_intensity;
     }
 
 }
diff --git a/Core/Tools/inc/LLData.h b/Core/Tools/inc/LLData.h
index 822ae227e6b37e21247c0a0dd1cfc25073a70e44..63cba17c3e47f5bfbc8f39cd4a3025188ac384c3 100644
--- a/Core/Tools/inc/LLData.h
+++ b/Core/Tools/inc/LLData.h
@@ -15,6 +15,7 @@
 //! @date   Nov 7, 2012
 
 #include "Exceptions.h"
+#include "Numeric.h"
 
 #include <algorithm>
 
@@ -42,6 +43,7 @@ public:
     LLData<T> &operator+=(const LLData<T> &right);
     LLData<T> &operator-=(const LLData<T> &right);
     LLData<T> &operator*=(const LLData<T> &right);
+    LLData<T> &operator/=(const LLData<T> &right);
 
     // initialization, scaling
     void setAll(const T &value);
@@ -68,6 +70,7 @@ private:
 template <class T> LLData<T> &operator+(const LLData<T> &left, const LLData<T> &right);
 template <class T> LLData<T> &operator-(const LLData<T> &left, const LLData<T> &right);
 template <class T> LLData<T> &operator*(const LLData<T> &left, const LLData<T> &right);
+template <class T> LLData<T> &operator/(const LLData<T> &left, const LLData<T> &right);
 
 // Global helper functions for comparison
 template <class T> bool HaveSameDimensions(const LLData<T> &left, const LLData<T> &right);
@@ -162,6 +165,25 @@ template<class T> LLData<T>& LLData<T>::operator*=(const LLData& right)
     return *this;
 }
 
+template<class T> LLData<T>& LLData<T>::operator/=(const LLData& right)
+{
+    if (!HaveSameDimensions(*this, right)) {
+        throw RuntimeErrorException("Operation /= on LLData requires both operands to have the same dimensions");
+    }
+    for (size_t i=0; i<getTotalSize(); ++i) {
+        double ratio(0);
+        if( std::abs(m_data_array[i]) <= Numeric::double_epsilon && std::abs(right[i]) <= Numeric::double_epsilon) {
+            ratio = 0.0;
+        } else if (std::abs(right[i]) <= Numeric::double_epsilon) {
+            ratio = m_data_array[i]/Numeric::double_epsilon;
+        } else {
+            ratio = m_data_array[i]/right[i];
+        }
+        m_data_array[i] /= ratio;
+    }
+    return *this;
+}
+
 template<class T> void LLData<T>::setAll(const T& value)
 {
     for (size_t i=0; i<getTotalSize(); ++i) {
@@ -271,6 +293,13 @@ template<class T> LLData<T> &operator*(const LLData<T>& left, const LLData<T>& r
     return *p_result;
 }
 
+template<class T> LLData<T> &operator/(const LLData<T>& left, const LLData<T>& right)
+{
+    LLData<T> *p_result = new LLData<T>(left);
+    (*p_result) /= right;
+    return *p_result;
+}
+
 template<class T> bool HaveSameDimensions(const LLData<T>& left, const LLData<T>& right)
 {
     if (left.getRank() != right.getRank()) {
diff --git a/Core/Tools/inc/OutputData.h b/Core/Tools/inc/OutputData.h
index cbe1db7996b2291c2c0504c257e613618ea90208..87f30554bd54aac612081589aff6ad15de3bc62c 100644
--- a/Core/Tools/inc/OutputData.h
+++ b/Core/Tools/inc/OutputData.h
@@ -24,54 +24,6 @@
 #include <string>
 #include <sstream>
 
-template <class T> class OutputData;
-
-
-//- -------------------------------------------------------------------
-//! @class MultiIndex
-//! @brief Definition of MultiIndex class to navigate trough OutputData
-//- -------------------------------------------------------------------
-class MultiIndex
-{
-    template <class T> friend class OutputData;
-public:
-    std::vector<std::string> getLabels() { return m_labels; }
-    std::vector<size_t> getCurrentIndices();
-    size_t getCurrentIndexOfAxis(std::string axis_name);
-    size_t getPosition() const { return m_current_position; }
-    size_t getSize() const { return m_total_size; }
-    std::vector<size_t> getAllSizes() { return m_axis_sizes; }
-    void reset();
-    void setIndexOfAxis(std::string axis_name, size_t value);
-    void setAllIndices(std::vector<size_t> indices);
-    void incrementIndexOfAxis(std::string axis_name);
-    void decrementIndexOfAxis(std::string axis_name);
-    MultiIndex& operator++();
-    bool endPassed() { return m_end_passed; }
-private:
-    MultiIndex();
-    ~MultiIndex();
-    // Disabling copy constructor and assignment
-    MultiIndex(const MultiIndex& source);
-    MultiIndex operator=(const MultiIndex& source);
-    void updateCurrentIndices();
-    void updateCurrentPosition();
-    void setPosition(size_t position);
-
-    void init(const std::vector<NamedVectorBase*>& value_axes);
-    void clear();
-    std::vector<std::string> m_labels;
-    std::map<std::string, size_t> m_label_index_map;
-    size_t m_dimension;
-    size_t m_total_size;
-    size_t m_current_position;
-    std::vector<size_t> m_axis_sizes;
-    std::vector<size_t> m_current_coordinate;
-    std::vector<size_t> m_steps;
-    bool m_end_passed;
-};
-
-
 //- -------------------------------------------------------------------
 //! @class OutputData
 //! @brief Definition of OutputData class to store data of any type in
@@ -85,7 +37,6 @@ public:
     //! make object clone
     OutputData* clone() const;
 
-//    void copyInto(OutputData<T> &x) const;
     void copyFrom(const OutputData<T> &x);
 
     void addAxis(NamedVectorBase* p_new_axis);
@@ -95,6 +46,10 @@ public:
     const NamedVectorBase *getAxis(std::string label) const;
     const NamedVectorBase *getAxis(size_t index) const;
 
+    // ---------------------------------
+    // retrieve basic info
+    // ---------------------------------
+
     //! return number of dimensions
     size_t getRank() const { return m_value_axes.size(); }
 
@@ -116,8 +71,8 @@ public:
     //! fill raw array with data
     void fillRawDataArray(T *destination) const;
 
-    //! return multi index
-    MultiIndex& getIndex() const;
+    //! get sum of all values in the data structure
+    T totalSum() const;
 
     // ---------------------------------
     // external iterators
@@ -139,41 +94,26 @@ public:
     const_iterator begin() const { return const_iterator(this); }
 
     //! return a read/write iterator that points to the one past last element
-    iterator end() { return iterator(this, getAllocatedSize()); }
+    const iterator end() { return iterator(this, getAllocatedSize()); }
 
     //! return a read-only iterator that points to the one past last element
-    const_iterator end() const  { return const_iterator(this, getAllocatedSize()); }
+    const const_iterator end() const  { return const_iterator(this, getAllocatedSize()); }
 
     // ---------------------------------
-    // navigation and access to stored elements
+    // coordinate and index functions
     // ---------------------------------
 
-    //! reset index to point to the beginning of data structure
-    void resetIndex() const;
-
-    //! return true if next element exists
-    bool hasNext() const;
+    //! return vector of coordinates for given index
+    std::vector<int> toCoordinates(size_t index) const;
 
-    //! return const reference to current value and increment index
-    const T& next() const;
+    //! return index for specified coordinates
+    size_t toIndex(std::vector<int> coordinates) const;
 
-    //! return reference to current value and increment index
-    T& next();
+    //! return index of axis with given name for given total index
+    size_t getIndexOfAxis(std::string axis_name, size_t total_index) const;
 
-    //! return const reference to current value
-    const T& currentValue() const;
-
-    //! return reference to current value
-    T& currentValue();
-
-    //! return current index of axis with given name
-    size_t getCurrentIndexOfAxis(std::string axis_name) const;
-
-    //! return current value of axis with given name
-    template <class U> U getCurrentValueOfAxis(std::string axis_name) const;
-
-    //! get sum of all values in the data structure
-    T totalSum() const;
+    //! return value of axis with given name at given index
+    template <class U> U getValueOfAxis(std::string axis_name, size_t index) const;
 
     // ---------
     // modifiers
@@ -197,6 +137,19 @@ public:
     //! set new values to raw data array
     void setRawDataArray(const T *source);
 
+    //! addition-assignment operator for two output data
+    const OutputData<T> &operator+=(const OutputData<T> &right);
+
+    //! substraction-assignment operator for two output data
+    const OutputData<T> &operator-=(const OutputData<T> &right);
+
+    //! division-assignment operator for two output data
+    const OutputData<T> &operator/=(const OutputData<T> &right);
+
+    //! multiplication-assignment operator for two output data
+    const OutputData<T> &operator*=(const OutputData<T> &right);
+
+
 private:
     //! hidden copy constructor and assignment operators
     OutputData(const OutputData& source);
@@ -205,14 +158,19 @@ private:
     //! memory allocation for current dimensions configuration
     void allocate();
 
-    //! accessors for iterators
+    //! accessor for iterators
     T &operator[](size_t index) {
         if (mp_ll_data) return (*mp_ll_data)[index];
         throw ClassInitializationException("Low-level data objects was not yet initialized");
     }
 
+    //! constant accessor for iterators
+    const T &operator[](size_t index) const {
+        if (mp_ll_data) return (*mp_ll_data)[index];
+        throw ClassInitializationException("Low-level data objects was not yet initialized");
+    }
+
     std::vector<NamedVectorBase*> m_value_axes;
-    mutable MultiIndex m_index;
     LLData<T> *mp_ll_data;
 };
 
@@ -221,23 +179,12 @@ private:
 // specialized OutputData functions: global arithmetics
 /* **************** */
 
-//! addition-assignment operator for two output data
-const OutputData<double> &operator+=(OutputData<double> &left, const OutputData<double> &right);
-
-//! substraction-assignment operator for two output data
-const OutputData<double> &operator-=(OutputData<double> &left, const OutputData<double> &right);
-
-//! division-assignment operator for two output data
-const OutputData<double> &operator/=(OutputData<double> &left, const OutputData<double> &right);
-
-//! multiplication-assignment operator for two output data
-const OutputData<double> &operator*=(OutputData<double> &left, const OutputData<double> &right);
-
 //! double the bin size for each dimension
 OutputData<double> *doubleBinSize(const OutputData<double> &source);
 
-//! unnormalized Fourier transformations for real data
+//! unnormalized Fourier transformation for real data
 void fourierTransform(const OutputData<double> &source, OutputData<complex_t> *p_destination);
+//! unnormalized reverse Fourier transformation for real data
 void fourierTransformR(const OutputData<complex_t> &source, OutputData<double> *p_destination);
 
 OutputData<double> *getRealPart(const OutputData<complex_t> &source);
@@ -248,20 +195,17 @@ OutputData<double> *getModulusPart(const OutputData<complex_t> &source);
 // definitions
 /* ***************************************************************************/
 
-// default constructor
 template <class T> OutputData<T>::OutputData()
 : mp_ll_data(0)
 {
     allocate();
 }
 
-// destructor
 template <class T> OutputData<T>::~OutputData()
 {
     clear();
 }
 
-// make object clone
 template <class T> OutputData<T>* OutputData<T>::clone() const
 {
 	OutputData<T>* p_result = new OutputData<T>();
@@ -340,7 +284,6 @@ inline std::vector<T> OutputData<T>::getRawDataVector() const
     return result;
 }
 
-//! fill raw array with data
 template <class T> void OutputData<T>::fillRawDataArray(T *destination) const
 {
     for (size_t i=0; i<getAllocatedSize(); ++i) {
@@ -349,83 +292,74 @@ template <class T> void OutputData<T>::fillRawDataArray(T *destination) const
     return;
 }
 
-
-// return multi index
-template <class T> inline MultiIndex& OutputData<T>::getIndex() const
-{
-    return m_index;
-}
-
-
-// reset index to point to the beginning of data structure
-template <class T> inline void OutputData<T>::resetIndex() const
-{
-    m_index.reset();
-}
-
-
-// return true if next element exists
-template <class T> inline bool OutputData<T>::hasNext() const
-{
-    return !m_index.endPassed();
-}
-
-
-// return const reference to current value and increment index
-template <class T> inline const T& OutputData<T>::next() const
+template<class T> std::vector<int> OutputData<T>::toCoordinates(size_t index) const
 {
-    const T &temp = currentValue();
-    ++m_index;
-    return temp;
-}
-
-
-// return reference to current value and increment index
-template <class T> inline T& OutputData<T>::next()
-{
-    T &temp = currentValue();
-    ++m_index;
-    return temp;
+    size_t remainder = index;
+    std::vector<int> result;
+    result.resize(mp_ll_data->getRank());
+    for (size_t i=0; i<mp_ll_data->getRank(); ++i)
+    {
+        result[mp_ll_data->getRank()-1-i] = remainder % m_value_axes[mp_ll_data->getRank()-1-i]->getSize();
+        remainder /= m_value_axes[mp_ll_data->getRank()-1-i]->getSize();
+    }
+    return result;
 }
 
-
-// return const reference to current value
-template <class T> inline const T& OutputData<T>::currentValue() const
+template <class T> size_t OutputData<T>::toIndex(std::vector<int> coordinates) const
 {
-    return (*mp_ll_data)[m_index.m_current_position];
+    if (coordinates.size() != mp_ll_data->getRank()) {
+        throw LogicErrorException("Number of coordinates must match rank of data structure");
+    }
+    size_t result = 0;
+    int step_size = 1;
+    for (size_t i=mp_ll_data->getRank(); i>0; --i)
+    {
+        result += coordinates[i-1]*step_size;
+        step_size *= m_value_axes[i-1]->getSize();
+    }
+    return result;
 }
 
-
-// return reference to current value
-template <class T> inline T& OutputData<T>::currentValue()
+template <class T> size_t OutputData<T>::getIndexOfAxis(std::string axis_name, size_t total_index) const
 {
-    return (*mp_ll_data)[m_index.m_current_position];
-}
-
-
-template <class T> inline size_t OutputData<T>::getCurrentIndexOfAxis(std::string axis_name) const
-{
-    return m_index.getCurrentIndexOfAxis(axis_name);
+    std::vector<int> coordinates = toCoordinates(total_index);
+    const NamedVectorBase *p_axis;
+    for (size_t i=0; i<m_value_axes.size(); ++i) {
+        p_axis = m_value_axes[i];
+        if (p_axis->getName() == axis_name) {
+            return coordinates[i];
+        }
+    }
+    throw LogicErrorException("Axis with given name not found");
 }
 
 
-// return current value of axis with given name
 template <class T>
-template <class U> inline U OutputData<T>::getCurrentValueOfAxis(std::string axis_name) const
+template<class U> U OutputData<T>::getValueOfAxis(std::string axis_name, size_t index) const
 {
-    const NamedVector<U> *p_axis = dynamic_cast<const NamedVector<U>*>(getAxis(axis_name));
-    size_t index = getCurrentIndexOfAxis(axis_name);
-    return (*p_axis)[index];
+    std::vector<int> coordinates = toCoordinates(index);
+    const NamedVectorBase *p_axis;
+    const NamedVector<U> *p_derived = 0;
+    size_t axis_index = 0;
+    for (size_t i=0; i<m_value_axes.size(); ++i) {
+        p_axis = m_value_axes[i];
+        if (p_axis->getName() == axis_name) {
+            p_derived = dynamic_cast<const NamedVector<U>*>(getAxis(axis_name));
+            axis_index = i;
+        }
+    }
+    if (p_derived == 0) {
+        throw LogicErrorException("Axis with given name not found");
+    }
+    return (*p_derived)[coordinates[axis_index]];
 }
 
-
 template<class T>
 inline T OutputData<T>::totalSum() const
 {
     return mp_ll_data->getTotalSum();
 }
 
-// set object into initial state (no dimensions, data)
 template <class T> void OutputData<T>::clear()
 {
     for (size_t i=0; i<getRank(); ++i)
@@ -437,21 +371,16 @@ template <class T> void OutputData<T>::clear()
     mp_ll_data = 0;
 }
 
-
-// set content of output data to specific value
 template <class T> void OutputData<T>::setAllTo(const T &value)
 {
     mp_ll_data->setAll(value);
 }
 
-
-// multiply every item of this output data by value
 template <class T> void OutputData<T>::scaleAll(const T &factor)
 {
     mp_ll_data->scaleAll(factor);
 }
 
-// add <rank> axes with indicated sizes
 template <class T> void OutputData<T>::setAxisSizes(size_t rank, int *n_dims)
 {
     clear();
@@ -463,6 +392,29 @@ template <class T> void OutputData<T>::setAxisSizes(size_t rank, int *n_dims)
     }
 }
 
+template<class T> const OutputData<T> &OutputData<T>::operator+=(const OutputData<T>& right)
+{
+    *this->mp_ll_data += *right.mp_ll_data;
+    return *this;
+}
+
+template<class T> const OutputData<T> &OutputData<T>::operator-=(const OutputData<T>& right)
+{
+    *this->mp_ll_data -= *right.mp_ll_data;
+    return *this;
+}
+
+template<class T> const OutputData<T> &OutputData<T>::operator*=(const OutputData<T>& right)
+{
+    *this->mp_ll_data *= *right.mp_ll_data;
+    return *this;
+}
+
+template<class T> const OutputData<T> &OutputData<T>::operator/=(const OutputData<T>& right)
+{
+    *this->mp_ll_data /= *right.mp_ll_data;
+    return *this;
+}
 
 template <class T> void OutputData<T>::allocate()
 {
@@ -476,11 +428,8 @@ template <class T> void OutputData<T>::allocate()
     T default_value = T();
     mp_ll_data->setAll(default_value);
     delete[] dims;
-    m_index.init(m_value_axes);
 }
 
-
-// set new values to raw data vector
 template<class T> inline void OutputData<T>::setRawDataVector(const std::vector<T> &data_vector)
 {
     if (data_vector.size() != getAllocatedSize()) {
@@ -491,7 +440,6 @@ template<class T> inline void OutputData<T>::setRawDataVector(const std::vector<
     }
 }
 
-// set new values to raw data array
 template<class T> inline void OutputData<T>::setRawDataArray(const T *source)
 {
     for (size_t i=0; i<getAllocatedSize(); ++i) {
diff --git a/Core/Tools/inc/OutputDataIterator.h b/Core/Tools/inc/OutputDataIterator.h
index 40ea569dfd9d5e94325d452483dc5a2100161115..cde0a27a4aeb56520be1d4ecf95fe4d899e24324 100644
--- a/Core/Tools/inc/OutputDataIterator.h
+++ b/Core/Tools/inc/OutputDataIterator.h
@@ -14,8 +14,6 @@
 //! @author Scientific Computing Group at FRM II
 //! @date   Nov 12, 2012
 
-template <class T> class OutputData;
-
 template <class TValue> class IIterator
 {
 public:
@@ -27,14 +25,6 @@ public:
         TValue *p_new = new TValue();
         return *p_new;
     }
-
-    // comparison
-    bool operator==(const IIterator<TValue> &other) const {
-        return typeid(*this)==typeid(other) && equal(other);
-    }
-    bool operator!=(const IIterator<TValue> &other) const { return !(*this == other); }
-protected:
-    virtual bool equal(const IIterator<TValue> &other) const { (void)other; return true; }
 };
 
 template <class TValue, class TContainer> class OutputDataIterator : public IIterator<TValue>
@@ -45,24 +35,43 @@ public:
         m_current_index(start_at_index), mp_output_data(p_output_data) {}
 
     //! templated copy construction
-    template<class TValue2, class TContainer2> OutputDataIterator(const OutputDataIterator<TValue2, TContainer2> &other);
+    template<class TValue2, class TContainer2> OutputDataIterator(
+            const OutputDataIterator<TValue2, TContainer2> &other);
 
     //! templated copy assignment
-    template<class TValue2, class TContainer2> OutputDataIterator<TValue, TContainer> &operator=(const OutputDataIterator<TValue2, TContainer2> &right);
+    template<class TValue2, class TContainer2> OutputDataIterator<TValue, TContainer> &operator=(
+            const OutputDataIterator<TValue2, TContainer2> &right);
 
     virtual ~OutputDataIterator();
 
-    //! unary increment
+    //! prefix increment
     OutputDataIterator<TValue, TContainer> &operator++();
 
+    //! postfix increment
+    OutputDataIterator<const TValue, TContainer> operator++(int);
+
     //! retrieve current element
     virtual TValue &operator*() const;
 
+    //! pointer access
+    virtual TValue* operator->() const;
+
     //! retrieve indexed element
     TValue &operator[](size_t index) const;
 
+    //! get current index
+    const size_t getIndex() const { return m_current_index; }
+
+    //! get container pointer
+    TContainer *const getContainer() const { return mp_output_data; }
+
+    //! comparison
+    template <class TValue2, class TContainer2> bool operator==(
+            const OutputDataIterator<TValue2, TContainer2> &other);
+    template <class TValue2, class TContainer2> bool operator!=(
+            const OutputDataIterator<TValue2, TContainer2> &other) { return !(*this == other); }
+
 protected:
-    virtual bool equal(const IIterator<TValue> &other) const;
     virtual void swapContents(OutputDataIterator<TValue, TContainer> &other);
     size_t m_current_index;
     TContainer *mp_output_data;
@@ -74,8 +83,8 @@ OutputDataIterator<TValue, TContainer>::OutputDataIterator(const OutputDataItera
 : m_current_index(0)
 , mp_output_data(0)
 {
-    mp_output_data = other.mp_output_data;
-    m_current_index = other.m_current_index;
+    mp_output_data = other.getContainer();
+    m_current_index = other.getIndex();
 }
 
 template<class TValue, class TContainer>
@@ -97,20 +106,35 @@ template<class TValue, class TContainer> OutputDataIterator<TValue, TContainer>
     return *this;
 }
 
+template<class TValue, class TContainer> OutputDataIterator<const TValue, TContainer> OutputDataIterator<TValue, TContainer>::operator ++(int dummy)
+{
+    (void)dummy;
+    OutputDataIterator<const TValue, TContainer> result = *this;
+    ++(*this);
+    return result;
+}
+
 template<class TValue, class TContainer> TValue &OutputDataIterator<TValue, TContainer>::operator*() const
 {
     return (*this)[m_current_index];
 }
 
+template<class TValue, class TContainer> TValue* OutputDataIterator<TValue, TContainer>::operator ->() const
+{
+    return &((*this)[m_current_index]);
+}
+
 template<class TValue, class TContainer> TValue &OutputDataIterator<TValue, TContainer>::operator[](size_t index) const
 {
     return (*mp_output_data)[index];
 }
 
-template<class TValue, class TContainer> bool OutputDataIterator<TValue, TContainer>::equal(const IIterator<TValue>& other) const
+template<class TValue, class TContainer>
+template<class TValue2, class TContainer2>
+bool OutputDataIterator<TValue, TContainer>::operator ==(
+        const OutputDataIterator<TValue2, TContainer2>& other)
 {
-    const OutputDataIterator<TValue, TContainer> &other_cast = static_cast<const OutputDataIterator<TValue, TContainer> &>(other);
-    return m_current_index==other_cast.m_current_index;
+    return mp_output_data==other.getContainer() && m_current_index==other.getIndex();
 }
 
 template<class TValue, class TContainer> void OutputDataIterator<TValue, TContainer>::swapContents(OutputDataIterator<TValue, TContainer>& other)
diff --git a/Core/Tools/src/FitSuiteStrategy.cpp b/Core/Tools/src/FitSuiteStrategy.cpp
index b5cb929d90ba34e561bdb7171d5464de5ea038dc..fe49de623aa6d712ea15f453b88d84c49cbd97bb 100644
--- a/Core/Tools/src/FitSuiteStrategy.cpp
+++ b/Core/Tools/src/FitSuiteStrategy.cpp
@@ -253,15 +253,15 @@ void FitSuiteStrategyBootstrap::setFitSuiteParameterValues(const std::vector<dou
 OutputData<double> *FitSuiteStrategyBootstrap::generateNoisyData(double noise_factor, const OutputData<double> &source)
 {
     OutputData<double> *p_result = source.clone();
-    p_result->resetIndex();
-    while (p_result->hasNext()) {
-        double current = p_result->currentValue();
+    OutputData<double>::iterator it = p_result->begin();
+    while (it != p_result->end()) {
+        double current = *it;
         double sigma = noise_factor*std::sqrt(current);
         double random = MathFunctions::GenerateNormalRandom(current, sigma);
         if (random<0.0) random = 0.0;
-        p_result->next() = random;
+        *it = random;
+        ++it;
     }
-
     return p_result;
 }
 
diff --git a/Core/Tools/src/OutputData.cpp b/Core/Tools/src/OutputData.cpp
index 6930867b0ee1eaf55ea4cfc113071aef135104dd..4780abcb6b1f33b5daa1e6f3deeb5ded1e13a4a2 100644
--- a/Core/Tools/src/OutputData.cpp
+++ b/Core/Tools/src/OutputData.cpp
@@ -8,316 +8,6 @@
 void toFftw3Array(complex_t *source, size_t length, fftw_complex *destination);
 void fromFftw3Array(fftw_complex *source, size_t length, complex_t *destination);
 
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-MultiIndex::MultiIndex()
-    : m_dimension(0)
-    , m_total_size(1)
-    , m_current_position(0)
-    , m_end_passed(false)
-{
-}
-
-MultiIndex::~MultiIndex()
-{
-    clear();
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-MultiIndex& MultiIndex::operator++()
-{
-    if (m_current_position<m_total_size-1)
-    {
-        ++m_current_position;
-    }
-    else
-    {
-        m_end_passed = true;
-    }
-    return *this;
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-std::vector<size_t> MultiIndex::getCurrentIndices()
-{
-    updateCurrentIndices();
-    return m_current_coordinate;
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-size_t MultiIndex::getCurrentIndexOfAxis(std::string axis_name)
-{
-    if (m_label_index_map.count(axis_name) == 0)
-    {
-        std::string message = "MultiIndex::getCurrentIndexOfAxis() -> Error! Label not found: " + axis_name;
-        throw NullPointerException(message);
-    }
-    return getCurrentIndices()[m_label_index_map.find(axis_name)->second];
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::reset()
-{
-    m_current_position = 0;
-    m_end_passed = false;
-    updateCurrentIndices();
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::updateCurrentIndices()
-{
-    size_t remainder = m_current_position;
-    for (size_t i=0; i<m_dimension; ++i)
-    {
-        m_current_coordinate[m_dimension-1-i] = remainder % m_axis_sizes[m_dimension-1-i];
-        remainder /= m_axis_sizes[m_dimension-1-i];
-    }
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::updateCurrentPosition()
-{
-    m_end_passed = false;
-    m_current_position = 0;
-    for (size_t i=0; i<m_dimension; ++i)
-    {
-        m_current_position += m_current_coordinate[i]*m_steps[i];
-    }
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::setPosition(size_t position)
-{
-    if (position>=m_total_size) {
-        throw OutOfBoundsException("MultiIndex::setPosition() -> Error! Position value out of bounds!");
-    }
-    m_current_position = position;
-    updateCurrentIndices();
-}
-
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::setIndexOfAxis(std::string axis_name, size_t value)
-{
-    if (m_label_index_map.count(axis_name) == 0) return;
-    size_t index = m_label_index_map[axis_name];
-    if (value >= m_axis_sizes[index])
-    {
-        throw OutOfBoundsException(" MultiIndex::setIndexOfAxis() -> Coordinate value out of bounds!");
-    }
-    m_current_coordinate[index] = value;
-    updateCurrentPosition();
-}
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::setAllIndices(std::vector<size_t> indices)
-{
-    if (indices.size()!=m_dimension) {
-        throw OutOfBoundsException(" MultiIndex::setAllIndices() -> Indices vector not of the right size!");
-    }
-    m_current_coordinate = indices;
-    updateCurrentPosition();
-}
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::incrementIndexOfAxis(std::string axis_name)
-{
-    if (m_label_index_map.count(axis_name) == 0) return;
-    size_t index = m_label_index_map[axis_name];
-    if (m_current_coordinate[index] < m_axis_sizes[index]-2)
-    {
-        ++m_current_coordinate[index];
-        updateCurrentPosition();
-        return;
-    }
-    throw OutOfBoundsException("MultiIndex::incrementIndexOfAxis() -> Error! Coordinate value out of bounds!");
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::decrementIndexOfAxis(std::string axis_name)
-{
-    if (m_label_index_map.count(axis_name) == 0) return;
-    size_t index = m_label_index_map[axis_name];
-    if (m_current_coordinate[index] > 1)
-    {
-        --m_current_coordinate[index];
-        updateCurrentPosition();
-        return;
-    }
-    throw OutOfBoundsException("MultiIndex::decrementIndexOfAxis() -> Error! Coordinate value out of bounds!");
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::init(const std::vector<NamedVectorBase*>& value_axes)
-{
-    clear();
-    m_dimension = value_axes.size();
-    if (m_dimension==0) return;
-    m_axis_sizes.resize(m_dimension);
-    m_current_coordinate.resize(m_dimension);
-    m_steps.resize(m_dimension);
-    for (size_t i=0; i<m_dimension; ++i)
-    {
-        NamedVectorBase* p_axis = value_axes[i];
-        if (m_label_index_map.count(p_axis->getName()) != 0)
-        {
-            throw ClassInitializationException("MultiIndex::init() -> Error! Axis labels for OutputData object should be unique!");
-        }
-        m_label_index_map[p_axis->getName()] = i;
-        m_labels.push_back(p_axis->getName());
-        m_axis_sizes[i] = p_axis->getSize();
-        m_total_size *= p_axis->getSize();
-        m_current_coordinate[i] = 0;
-    }
-    size_t remaining_size = m_total_size;
-    for (size_t i=0; i<m_dimension; ++i)
-    {
-        m_steps[i] = remaining_size / m_axis_sizes[i];
-        remaining_size = m_steps[i];
-    }
-}
-
-
-/* ************************************************************************* */
-//
-/* ************************************************************************* */
-void MultiIndex::clear()
-{
-    m_axis_sizes.clear();
-    m_current_coordinate.clear();
-    m_labels.clear();
-    m_label_index_map.clear();
-    m_steps.clear();
-    m_dimension = 0;
-    m_total_size = 1;
-    m_current_position = 0;
-    m_end_passed = false;
-}
-
-
-/* ************************************************************************* */
-// addition-assignment operator for two output data
-/* ************************************************************************* */
-const OutputData<double> &operator+=(OutputData<double> &left, const OutputData<double> &right)
-{
-    if( &left == &right) throw LogicErrorException("OutputData &operator+=() -> Error. Left And right can't be the same");
-
-    size_t total_size = left.getAllocatedSize();
-    if (right.getAllocatedSize()!= total_size) {
-        throw LogicErrorException("OutputData<double> &operator+=() -> Error! Cannot add OutputData objects of different size.");
-    }
-    left.resetIndex();
-    right.resetIndex();
-    while (left.hasNext()) {
-        left.next() += right.next();
-    }
-    return left;
-}
-
-
-/* ************************************************************************* */
-// substraction-assignment operator for two output data
-/* ************************************************************************* */
-const OutputData<double> &operator-=(OutputData<double> &left, const OutputData<double> &right)
-{
-    if( &left == &right) throw LogicErrorException("OutputData &operator-=() -> Error. Left And right can't be the same");
-
-    size_t total_size = left.getAllocatedSize();
-    if (right.getAllocatedSize()!= total_size) {
-        throw LogicErrorException("OutputData<double> &operator-=() -> Error! Cannot substract OutputData objects of different size.");
-    }
-    left.resetIndex();
-    right.resetIndex();
-    while (right.hasNext()) {
-        left.next() -= right.next();
-    }
-    return left;
-}
-
-
-/* ************************************************************************* */
-// division-assignment operator for two output data
-/* ************************************************************************* */
-const OutputData<double> &operator/=(OutputData<double> &left, const OutputData<double> &right)
-{
-    if( &left == &right) throw LogicErrorException("OutputData &operator/=() -> Error. Left And right can't be the same");
-
-    size_t total_size = left.getAllocatedSize();
-    if (right.getAllocatedSize()!= total_size) {
-        throw LogicErrorException("OutputData<double> &operator/=() -> Error! Cannot substract OutputData objects of different size.");
-    }
-    left.resetIndex();
-    right.resetIndex();
-    while ( right.hasNext() ) {
-        double xleft = left.currentValue();
-        double xright = right.next();
-        double ratio(0);
-        if( fabs(xleft) <= Numeric::double_epsilon && fabs(xright) <= Numeric::double_epsilon) {
-            ratio = 0.0;
-        } else if (fabs(xright) <= Numeric::double_epsilon) {
-            ratio = xleft/Numeric::double_epsilon;
-        } else {
-            ratio = xleft/xright;
-        }
-        left.next() = ratio;
-    }
-    return left;
-}
-
-/* ************************************************************************* */
-// multiplication-assignment operator for two output data
-/* ************************************************************************* */
-const OutputData<double> &operator*=(OutputData<double> &left, const OutputData<double> &right)
-{
-    if( &left == &right) throw LogicErrorException("OutputData &operator*=() -> Error. Left And right can't be the same");
-
-    size_t total_size = left.getAllocatedSize();
-    if (right.getAllocatedSize()!= total_size) {
-        throw LogicErrorException("OutputData<double> &operator/=() -> Error! Cannot substract OutputData objects of different size.");
-    }
-    left.resetIndex();
-    right.resetIndex();
-    while ( right.hasNext() ) {
-        left.next() *= right.next();
-    }
-    return left;
-}
-
 /* ************************************************************************* */
 // double the bin size for each dimension
 /* ************************************************************************* */
@@ -334,23 +24,21 @@ OutputData<double> *doubleBinSize(const OutputData<double> &source)
         p_result->addAxis(source_axis->createDoubleBinSize());
     }
     // calculate new data content
-    source.resetIndex();
-    MultiIndex &source_index = source.getIndex();
-    MultiIndex &result_index = p_result->getIndex();
-    while (source.hasNext()) {
-        std::vector<size_t> source_indices = source_index.getCurrentIndices();
-        std::vector<size_t> dest_indices;
+    OutputData<double>::const_iterator it_source = source.begin();
+    OutputData<double>::iterator it_result = p_result->begin();
+    while (it_source != source.end()) {
+        std::vector<int> source_indices = source.toCoordinates(it_source.getIndex());
+        std::vector<int> dest_indices;
         double boundary_factor = 1.0;
         for (size_t i=0; i<source_indices.size(); ++i) {
             dest_indices.push_back(source_indices[i]/2);
             if (needs_resizing[i] &&
                     source_sizes[i]%2 == 1 &&
-                    source_indices[i] == source_sizes[i]-1) {
+                    source_indices[i] == (int)source_sizes[i]-1) {
                 boundary_factor *= 2.0;
             }
         }
-        result_index.setAllIndices(dest_indices);
-        p_result->currentValue() = boundary_factor*source.next();
+        it_result[p_result->toIndex(dest_indices)] = boundary_factor*(*it_source++);
     }
     return p_result;
 }
@@ -443,10 +131,11 @@ OutputData<double> *getRealPart(const OutputData<complex_t> &source)
     for (size_t i=0; i<source.getRank(); ++i) {
         p_result->addAxis(source.getAxis(i)->clone());
     }
-    source.resetIndex();
-    p_result->resetIndex();
-    while (source.hasNext()) {
-        p_result->next() = source.next().real();
+    OutputData<complex_t>::const_iterator it_source = source.begin();
+    OutputData<double>::iterator it_result = p_result->begin();
+    while (it_source != source.end()) {
+        *it_result = it_source->real();
+        ++it_source, ++it_result;
     }
     return p_result;
 }
@@ -457,10 +146,11 @@ OutputData<double>* getImagPart(const OutputData<complex_t>& source)
     for (size_t i=0; i<source.getRank(); ++i) {
         p_result->addAxis(source.getAxis(i)->clone());
     }
-    source.resetIndex();
-    p_result->resetIndex();
-    while (source.hasNext()) {
-        p_result->next() = source.next().imag();
+    OutputData<complex_t>::const_iterator it_source = source.begin();
+    OutputData<double>::iterator it_result = p_result->begin();
+    while (it_source != source.end()) {
+        *it_result = it_source->real();
+        ++it_source, ++it_result;
     }
     return p_result;
 }
@@ -471,10 +161,11 @@ OutputData<double>* getModulusPart(const OutputData<complex_t>& source)
     for (size_t i=0; i<source.getRank(); ++i) {
         p_result->addAxis(source.getAxis(i)->clone());
     }
-    source.resetIndex();
-    p_result->resetIndex();
-    while (source.hasNext()) {
-        p_result->next() = std::abs(source.next());
+    OutputData<complex_t>::const_iterator it_source = source.begin();
+    OutputData<double>::iterator it_result = p_result->begin();
+    while (it_source != source.end()) {
+        *it_result = std::abs(*it_source);
+        ++it_source, ++it_result;
     }
     return p_result;
 }
diff --git a/Core/Tools/src/OutputDataReader.cpp b/Core/Tools/src/OutputDataReader.cpp
index 9d4c9f16dcae4b06f143385c563d7dd019192706..98cbaa8ff15c3e6e40f01518d7e244f6b3df1d5f 100644
--- a/Core/Tools/src/OutputDataReader.cpp
+++ b/Core/Tools/src/OutputDataReader.cpp
@@ -96,19 +96,19 @@ OutputData<double > *OutputDataReadStreamV1::readOutputData(std::istream &input_
     NamedVector<double> *yaxis = new NamedVector<double>("y-axis");
     for(size_t i=0; i<buff_yaxis.size(); ++i) yaxis->push_back(buff_yaxis[i]);
 
-    OutputData<double > *output = new OutputData<double >;
-    output->addAxis(xaxis);
-    output->addAxis(yaxis);
-    output->setAllTo(0.0);
-    output->resetIndex();
-    while (output->hasNext())
+    OutputData<double > *p_result = new OutputData<double>;
+    p_result->addAxis(xaxis);
+    p_result->addAxis(yaxis);
+    p_result->setAllTo(0.0);
+    OutputData<double>::iterator it = p_result->begin();
+    while (it != p_result->end())
     {
-        size_t index_x = output->getCurrentIndexOfAxis("x-axis");
-        size_t index_y = output->getCurrentIndexOfAxis("y-axis");
-        output->next() = buff_data[index_y][index_x];
+        size_t index_x = p_result->toCoordinates(it.getIndex())[0];
+        size_t index_y = p_result->toCoordinates(it.getIndex())[1];
+        *it = buff_data[index_y][index_x];
+        ++it;
     }
-
-    return output;
+    return p_result;
 }
 
 
diff --git a/UnitTests/TestCore/OutputDataTest.h b/UnitTests/TestCore/OutputDataTest.h
index 918f7b2b93dd30a712bdf8dfcc58981c01b6bea9..a42966899aed1c86c0afbfd562c2676ea7f1606a 100644
--- a/UnitTests/TestCore/OutputDataTest.h
+++ b/UnitTests/TestCore/OutputDataTest.h
@@ -25,11 +25,11 @@ OutputDataTest::OutputDataTest()
     db_data_3d.addAxis(new NamedVector<double>("angle", 0.0, 0.1, 20));
     db_data_3d.addAxis(new NamedVector<double>("length", 0.0, 0.5, 10));
     db_data_3d.addAxis(new NamedVector<int>("index", 10, 1, 10));
-    MultiIndex &db_data_index = db_data_3d.getIndex();
+    OutputData<double>::iterator it = db_data_3d.begin();
     for (size_t i=0; i<2000; ++i)
     {
-        db_data_3d.currentValue() = (double)i;
-        ++db_data_index;
+        *it = (double)i;
+        ++it;
     }
 }
 
diff --git a/UnitTests/TestCore/main.cpp b/UnitTests/TestCore/main.cpp
index 74b6ea1291a5f305e0b42271e5414bcfdb7191b6..a03d74453a9a975bc5c68152b3743f53643e13b6 100644
--- a/UnitTests/TestCore/main.cpp
+++ b/UnitTests/TestCore/main.cpp
@@ -37,11 +37,12 @@ TEST_F(OutputDataTest, SizeAfterAddingAxes)
 
 TEST_F(OutputDataTest, DataInitialization)
 {
-    MultiIndex &db_data_index = db_data_3d.getIndex();
-    db_data_index.setIndexOfAxis("angle", 11);
-    db_data_index.setIndexOfAxis("length", 4);
-    db_data_index.setIndexOfAxis("index", 3);
-    EXPECT_DOUBLE_EQ((double)1143, db_data_3d.currentValue());
+    OutputData<double>::const_iterator it = db_data_3d.begin();
+    std::vector<int> coordinates;
+    coordinates.push_back(11);
+    coordinates.push_back(4);
+    coordinates.push_back(3);
+    EXPECT_DOUBLE_EQ((double)1143, it[db_data_3d.toIndex(coordinates)]);
 }
 
 int main(int argc, char** argv)