diff --git a/App/App.pro b/App/App.pro
index 3f99c8d12c60f9c687f2a6f3a6c9a60eef109773..d37c5cd02c97b4c45bdcef77727eb2cfb5360901 100644
--- a/App/App.pro
+++ b/App/App.pro
@@ -33,9 +33,10 @@ SOURCES += \
     src/TestIsGISAXS9.cpp \
     src/TestIsGISAXS10.cpp \
     src/TestMesoCrystal.cpp \
+    src/TestMultiLayerRoughness.cpp \
+    src/TestPerformance.cpp \
     src/TestRootTree.cpp \
-    src/TestRoughness.cpp \
-    src/TestPerformance.cpp
+    src/TestRoughness.cpp
 
 HEADERS += \
     inc/App.h \
@@ -59,9 +60,10 @@ HEADERS += \
     inc/TestIsGISAXS9.h \
     inc/TestIsGISAXS10.h \
     inc/TestMesoCrystal.h \
+    inc/TestMultiLayerRoughness.h \
+    inc/TestPerformance.h \
     inc/TestRootTree.h \
-    inc/TestRoughness.h \
-    inc/TestPerformance.h
+    inc/TestRoughness.h
 
 INCLUDEPATH += ./inc ../Core/Algorithms/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc ../Core/PythonAPI/inc
 DEPENDPATH  += ./inc ../Core/Algorithms/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc ../Core/PythonAPI/inc
diff --git a/App/inc/TestMultiLayerRoughness.h b/App/inc/TestMultiLayerRoughness.h
new file mode 100644
index 0000000000000000000000000000000000000000..39af242402d2de78b28db08f5f1c8e631dffb418
--- /dev/null
+++ b/App/inc/TestMultiLayerRoughness.h
@@ -0,0 +1,37 @@
+#ifndef TESTMULTILAYERROUGHNESS_H
+#define TESTMULTILAYERROUGHNESS_H
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   TestMultiLayerRoughness.h
+//! @brief  Tests of specular and diffuse reflection from rough multilayer in DWBA.
+//! @author Scientific Computing Group at FRM II
+//! @date   01.08.2012
+
+#include "IFunctionalTest.h"
+#include "OutputData.h"
+#include "ISample.h"
+
+
+//- -------------------------------------------------------------------
+//! @class TestMultiLayerRoughness
+//! @brief Tests of specular and diffuse reflection from rough multilayer in DWBA.
+//! (this is replacement for TestDiffuseReflection test)
+//- -------------------------------------------------------------------
+class TestMultiLayerRoughness : public IFunctionalTest
+{
+public:
+    TestMultiLayerRoughness();
+    virtual ~TestMultiLayerRoughness(){}
+
+    virtual void execute();
+
+};
+
+#endif // TESTMULTILAYERROUGHNESS_H
diff --git a/App/src/CommandLine.cpp b/App/src/CommandLine.cpp
index 9355e749b683058ad60ab0a3617f7eb40ae148ae..1345cdbbf9844d7224f40bb20f0d83b7dc5a65de 100644
--- a/App/src/CommandLine.cpp
+++ b/App/src/CommandLine.cpp
@@ -37,19 +37,20 @@ CommandLine::CommandLine(int argc, char **argv) : m_argc(argc), m_argv(argv)
     m_description["profile"]      = "profile specified test";
     // default description of functional tests
     for(arguments_t::iterator it = m_defined_functional_tests.begin(); it!=m_defined_functional_tests.end(); ++it ) {
-        m_description[ (*it) ]    = "functional test: no description";
+        m_description[ (*it) ]        = "functional test: no description";
     }
     // additional description of functional tests (overwrites previous default)
-    m_description["isgisaxs1"]    = "functional test: isgisaxs ex-1 (2 types of particles without inteference on top of substrate)";
-    m_description["isgisaxs3"]    = "functional test: isgisaxs ex-3 (cylinder FF)";
-    m_description["isgisaxs9"]    = "functional test: isgisaxs ex-9 (pyramid FF)";
-    m_description["isgisaxs10"]   = "functional test: isgisaxs ex-10 (cylinders with interference on top of substrate)";
-    m_description["convolution"]  = "functional test: test of convolution via fft";
-    m_description["diffuse"]      = "functional test: diffuse scattering from multi layer with roughness";
-    m_description["formfactor"]   = "functional test: some formfactor";
-    m_description["roughness"]    = "functional test: roughness parameters";
-    m_description["roottree"]     = "functional test: using root trees to read/write data from/to disk";
-    m_description["performance"]  = "functional test: run performance test for several predefined tasks";
+    m_description["isgisaxs1"]        = "functional test: isgisaxs ex-1 (2 types of particles without inteference on top of substrate)";
+    m_description["isgisaxs3"]        = "functional test: isgisaxs ex-3 (cylinder FF)";
+    m_description["isgisaxs9"]        = "functional test: isgisaxs ex-9 (rotated pyramid FF)";
+    m_description["isgisaxs10"]       = "functional test: isgisaxs ex-10 (cylinders with interference on top of substrate)";
+    m_description["convolution"]      = "functional test: test of convolution via fft";
+    m_description["diffuse"]          = "functional test: diffuse scattering from multi layer with roughness";
+    m_description["formfactor"]       = "functional test: some formfactor";
+    m_description["roughness"]        = "functional test: roughness parameters";
+    m_description["roottree"]         = "functional test: using root trees to read/write data from/to disk";
+    m_description["performance"]      = "functional test: run performance test for several predefined tasks";
+    m_description["roughdwba"]        = "functional test: diffuse scattering from multi layer with roughness";
 
 }
 
diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index cb678bbf15bc42a4121f06eaaa84ac302be6969d..40383534dc35787e9865a6c6702fe7775138e2bc 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -39,6 +39,10 @@ void IsGISAXSTools::drawLogOutputData(const OutputData<double>& output,
 void IsGISAXSTools::drawOutputDataInPad(const OutputData<double>& output,
         const std::string& draw_options, const std::string &histogram_title)
 {
+    if(!gPad) {
+        throw NullPointerException("IsGISAXSTools::drawOutputDataInPad() -> Error! No canvas exists.");
+    }
+
     const OutputData<double> *p_output = &output;
     if (p_output->getDimension() != 2) return;
     // creation of 2D histogram from calculated intensities
@@ -92,6 +96,10 @@ void IsGISAXSTools::drawOutputDataInPad(const OutputData<double>& output,
 /* ************************************************************************* */
 void IsGISAXSTools::drawOutputDataDistribution1D(const OutputData<double> &output, const std::string &draw_options, const std::string &histogram_title)
 {
+    if(!gPad) {
+        throw NullPointerException("IsGISAXSTools::drawOutputDataDistribution1D() -> Error! No canvas exists.");
+    }
+
     std::string histo_name = histogram_title;
     if (histo_name.empty()) {
         histo_name = gPad->GetTitle();
@@ -116,6 +124,10 @@ void IsGISAXSTools::drawOutputDataDistribution1D(const OutputData<double> &outpu
 /* ************************************************************************* */
 void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, const OutputData<double> &right, const std::string &draw_options, const std::string &histogram_title)
 {
+    if(!gPad) {
+        throw NullPointerException("IsGISAXSTools::drawOutputDataDifference1D() -> Error! No canvas exists.");
+    }
+
     OutputData<double> *left_clone = left.clone();
     OutputData<double> *right_clone = right.clone();
 
@@ -158,6 +170,9 @@ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, c
 /* ************************************************************************* */
 void IsGISAXSTools::drawOutputDataDifference2D(const OutputData<double> &left, const OutputData<double> &right, const std::string &draw_options, const std::string &histogram_title)
 {
+    if(!gPad) {
+        throw NullPointerException("IsGISAXSTools::drawOutputDataDifference2D -> Error! No canvas exists.");
+    }
     OutputData<double> *left_clone = left.clone();
     OutputData<double> *right_clone = right.clone();
 
diff --git a/App/src/StandardSamples.cpp b/App/src/StandardSamples.cpp
index 9f26fb5c5f127f400bca0ee432adfd93ff08c8ab..484b3b14dfc9f39e56ae3692f1b2c29608fb2e61 100644
--- a/App/src/StandardSamples.cpp
+++ b/App/src/StandardSamples.cpp
@@ -140,7 +140,7 @@ ISample *StandardSamples::MultilayerOffspecTestcase1a()
     lSubstrate.setMaterial(mSubstrate, 0);
 
     LayerRoughness roughness;
-    roughness.setSigma(0.0*Units::nanometer);
+    roughness.setSigma(1.0*Units::nanometer);
     roughness.setHurstParameter(0.3);
     roughness.setLatteralCorrLength(500*Units::nanometer);
 
diff --git a/App/src/TestConvolution.cpp b/App/src/TestConvolution.cpp
index 3fcec4317663578c05959bae4f9c5a2e0f855450..d85927f521b2b8c5e6ada1ec7db3f58645666e78 100644
--- a/App/src/TestConvolution.cpp
+++ b/App/src/TestConvolution.cpp
@@ -177,7 +177,7 @@ void TestConvolution::test_convolve1d()
         href->DrawCopy();
         gr_result->Draw("pl same");
     }
-    float_t rp, cp;
+    Float_t rp, cp;
     std::cout << "--- TestConvolution::test_convolve_1d() -> Peformance." << std::endl;
     benchmark.Summary(rp, cp);
 
@@ -320,7 +320,7 @@ void TestConvolution::test_convolve2d()
     }
 
     // benchmark summary
-    float_t rp, cp;
+    Float_t rp, cp;
     std::cout << "--- TestConvolution::test_convolve_2d() -> Peformance." << std::endl;
     benchmark.Summary(rp, cp);
 
diff --git a/App/src/TestDiffuseReflection.cpp b/App/src/TestDiffuseReflection.cpp
index b8c5907af915899118bf10faee34a27b4fbfd397..9a15084a3bd3e823d07bfb1823c54102fbf73def 100644
--- a/App/src/TestDiffuseReflection.cpp
+++ b/App/src/TestDiffuseReflection.cpp
@@ -53,6 +53,8 @@ void TestDiffuseReflection::execute()
     kvector_t ki, kf;
     for(size_t i_sample=0; i_sample<samples.size(); i_sample++){
         m_sample = samples[i_sample];
+        std::cout << *m_sample << std::endl;
+        m_sample->walk_and_print();
 
         // specular reflectivity alpha_i = alpha_f
         m_data_spec = new OutputData<double >;
@@ -75,7 +77,7 @@ void TestDiffuseReflection::execute()
         m_data_offspec = new OutputData<double >;
         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_spec->resetIndex();
+        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");
@@ -84,6 +86,7 @@ void TestDiffuseReflection::execute()
             //double r = calc.execute0(*m_sample, ki, kf);
             calc.execute(*m_sample, ki, kf);
             m_data_offspec->next() = calc.getDiffuseAutocorr() + calc.getDiffuseCrosscorr();
+//            std::cout << alpha_i << " " << alpha_f << " " << calc.getDiffuseAutocorr() << std::endl;
         }
 
         draw();
@@ -149,7 +152,7 @@ void TestDiffuseReflection::draw()
     double dalpha = (m_alphaMax - m_alphaMin) / (m_npoints-1);
     TH2D h2("h2","h2", m_npoints, Units::rad2deg(m_alphaMin-dalpha/2.), Units::rad2deg(m_alphaMax+dalpha/2.), m_npoints, Units::rad2deg(m_alphaMin-dalpha/2.), Units::rad2deg(m_alphaMax+dalpha/2.) );
     h2.SetContour(50);
-    h2.SetMinimum(0.001);
+//    h2.SetMinimum(0.001);
     h2.SetStats(0);
     m_data_offspec->resetIndex();
     while (m_data_offspec->hasNext()) {
diff --git a/App/src/TestFactory.cpp b/App/src/TestFactory.cpp
index 240689f136be91df0b884a8fa682c99ef9f7fbcb..527090533580c028cb129c30e2266ab920c81381 100644
--- a/App/src/TestFactory.cpp
+++ b/App/src/TestFactory.cpp
@@ -13,6 +13,7 @@
 #include "TestRootTree.h"
 #include "TestFittingModule.h"
 #include "TestPerformance.h"
+#include "TestMultiLayerRoughness.h"
 
 
 #include "TBenchmark.h"
@@ -37,6 +38,7 @@ TestFactory::TestFactory() : m_benchmark(0)
     registerItem("roottree",    IFactoryCreateFunction<TestRootTree, IFunctionalTest> );
     registerItem("fitting",     IFactoryCreateFunction<TestFittingModule, IFunctionalTest> );
     registerItem("performance", IFactoryCreateFunction<TestPerformance, IFunctionalTest> );
+    registerItem("roughdwba",   IFactoryCreateFunction<TestMultiLayerRoughness, IFunctionalTest> );
 
     m_benchmark = new TBenchmark();
 }
@@ -54,7 +56,7 @@ TestFactory::~TestFactory()
 void TestFactory::print_benchmarks()
 {
     std::cout << "--- TestFactory::print_benchmarks() ---" << std::endl;
-    float_t rp, cp;
+    Float_t rp, cp;
     m_benchmark->Summary(rp, cp);
 }
 
diff --git a/App/src/TestIsGISAXS9.cpp b/App/src/TestIsGISAXS9.cpp
index b67a8fd48ebe542f1ed625de7a996bdc51e387c0..17f9a88000a52768cf61ab4aef230a9a7f003f3c 100644
--- a/App/src/TestIsGISAXS9.cpp
+++ b/App/src/TestIsGISAXS9.cpp
@@ -47,24 +47,18 @@ void TestIsGISAXS9::clear()
 void TestIsGISAXS9::execute()
 {
 
-    //MultiLayer *ml = dynamic_cast<MultiLayer *>(SampleFactory::createSample("IsGISAXS9_RotatedPyramid"));
-    MultiLayer *ml = dynamic_cast<MultiLayer *>(SampleFactory::createSample("MesoCrystal1"));
+//    //MultiLayer *ml = dynamic_cast<MultiLayer *>(SampleFactory::createSample("IsGISAXS9_RotatedPyramid"));
+//    MultiLayer *ml = dynamic_cast<MultiLayer *>(SampleFactory::createSample("MesoCrystal1"));
+//    std::cout << *ml << std::endl;
+//    std::cout << "------------------------------------" << std::endl;
+//    ml->walk_and_print();
+//    std::cout << "------------------------------------" << std::endl;
+//    ParameterPool *pool = ml->createParameterTree();
+//    std::cout << *pool << std::endl;
+//    std::cout << pool->setMatchedParametersValue("/*NanoParticleInfo*/depth",99);
+//    ml->walk_and_print();
+//    return;
 
-    std::cout << *ml << std::endl;
-
-    std::cout << "------------------------------------" << std::endl;
-    ml->walk_and_print();
-
-    std::cout << "------------------------------------" << std::endl;
-    ParameterPool *pool = ml->createParameterTree();
-    std::cout << *pool << std::endl;
-
-    std::cout << pool->setMatchedParametersValue("/*NanoParticleInfo*/depth",99);
-
-    ml->walk_and_print();
-
-
-    return;
     clear();
 
     // setting experiment
diff --git a/App/src/TestMultiLayerRoughness.cpp b/App/src/TestMultiLayerRoughness.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..97bc1064f1994e31ea43bc6f254179b874a22645
--- /dev/null
+++ b/App/src/TestMultiLayerRoughness.cpp
@@ -0,0 +1,47 @@
+#include "TestMultiLayerRoughness.h"
+#include "GISASExperiment.h"
+#include "Units.h"
+#include "IsGISAXSTools.h"
+#include "SampleFactory.h"
+
+#include "TCanvas.h"
+
+TestMultiLayerRoughness::TestMultiLayerRoughness()
+{
+
+}
+
+
+
+
+void TestMultiLayerRoughness::execute()
+{
+
+//    std::vector<std::string > snames;
+//    snames.push_back("MultilayerOffspecTestcase1a");
+//    snames.push_back("MultilayerOffspecTestcase1b");
+////    snames.push_back("MultilayerOffspecTestcase2a");
+////    snames.push_back("MultilayerOffspecTestcase2b");
+
+//    std::vector<MultiLayer *> samples;
+//    for(size_t i_sample=0; i_sample<snames.size(); i_sample++){
+//        samples.push_back( dynamic_cast<MultiLayer *>(SampleFactory::instance().createItem(snames[i_sample])) );
+//    }
+
+
+    ISample *sample = SampleFactory::instance().createItem("MultilayerOffspecTestcase1a");
+
+    // setting experiment
+    GISASExperiment experiment;
+    experiment.setDetectorParameters(100, 0.0*Units::degree, 2.0*Units::degree, 100, 0.0*Units::degree, 2.0*Units::degree, true);
+    experiment.setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
+
+    experiment.setSample(sample);
+    experiment.runSimulation();
+
+    const OutputData<double> *output = experiment.getOutputData();
+
+    TCanvas *c1 = new TCanvas("c1","c1",1024, 768);
+    c1->cd();
+    IsGISAXSTools::drawOutputDataInPad(*output, "CONT4 Z", "Calculated pyramid FF");
+}
diff --git a/Core/Algorithms/inc/MultiLayerDWBASimulation.h b/Core/Algorithms/inc/MultiLayerDWBASimulation.h
index e9c7c26265fe91ac4514762f3361a24c26ba8bc0..88b113246cc11b41b60866f3ad48a3423f46bc29 100644
--- a/Core/Algorithms/inc/MultiLayerDWBASimulation.h
+++ b/Core/Algorithms/inc/MultiLayerDWBASimulation.h
@@ -19,6 +19,8 @@
 #include <map>
 
 class MultiLayer;
+class MultiLayerRoughnessDWBASimulation;
+
 
 class MultiLayerDWBASimulation : public DWBASimulation
 {
@@ -32,6 +34,7 @@ public:
 protected:
     std::map<size_t, LayerDWBASimulation*> m_layer_dwba_simulation_map;
     MultiLayer *mp_multi_layer;
+    MultiLayerRoughnessDWBASimulation *m_roughness_dwba_simulation;
 };
 
 
diff --git a/Core/Algorithms/inc/MultiLayerRoughnessDWBASimulation.h b/Core/Algorithms/inc/MultiLayerRoughnessDWBASimulation.h
new file mode 100644
index 0000000000000000000000000000000000000000..15d2a3da43cc3ab536c3de8604c978f3f3b40799
--- /dev/null
+++ b/Core/Algorithms/inc/MultiLayerRoughnessDWBASimulation.h
@@ -0,0 +1,58 @@
+#ifndef MULTILAYERROUGHNESSDWBASIMULATION_H
+#define MULTILAYERROUGHNESSDWBASIMULATION_H
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   MultiLayerRoughnessDWBASimulation.h
+//! @brief  Calculation of diffuse reflection from multilayer with rough interfaces
+//! @author Scientific Computing Group at FRM II
+//! @date   02.08.2012
+
+#include "DWBASimulation.h"
+#include <vector>
+
+class MultiLayer;
+#include "IDoubleToComplexFunction.h"
+
+//- -------------------------------------------------------------------
+//! @class MultiLayerRoughnessDWBASimulation
+//! @brief Calculation of diffuse reflection from multilayer with rough interfaces
+//- -------------------------------------------------------------------
+class MultiLayerRoughnessDWBASimulation : public DWBASimulation
+{
+public:
+    MultiLayerRoughnessDWBASimulation(const MultiLayer *p_multi_layer);
+    virtual ~MultiLayerRoughnessDWBASimulation();
+
+    virtual void run();
+
+    // set Kz, T and R functions for given layer
+    void setKzTAndRFunctions(int i, const IDoubleToComplexFunction &kz_function, const IDoubleToComplexFunction &T_function, const IDoubleToComplexFunction &R_function);
+
+    // evaluate
+    virtual double evaluate(kvector_t k_i, kvector_t k_f);
+
+protected:
+    MultiLayerRoughnessDWBASimulation(const MultiLayerRoughnessDWBASimulation &);
+    MultiLayerRoughnessDWBASimulation &operator=(const MultiLayerRoughnessDWBASimulation &);
+
+    complex_t get_refractive_term(int ilayer);
+    complex_t get_sum4terms(int ilayer);
+
+    std::vector<IDoubleToComplexFunction *> mp_kz_function;
+    std::vector<IDoubleToComplexFunction *> mp_T_function;
+    std::vector<IDoubleToComplexFunction *> mp_R_function;
+
+    MultiLayer *mp_multi_layer;
+
+    double m_alpha_f;
+    kvector_t m_kf;
+};
+
+#endif // MULTILAYERROUGHNESSDWBASIMULATION_H
diff --git a/Core/Algorithms/inc/OpticalFresnel.h b/Core/Algorithms/inc/OpticalFresnel.h
index b589258df102857be74c05dd335b94385c400395..c92ddbbaf9efef9cc2ef9cacef12c910cd7fe150 100644
--- a/Core/Algorithms/inc/OpticalFresnel.h
+++ b/Core/Algorithms/inc/OpticalFresnel.h
@@ -67,11 +67,4 @@ private:
 };
 
 
-// class to help pyplusplus to expose MultiLayerCoeff_t in python during automatic code generation
-//class PyHelperOpticalFresnel
-//{
-//    size_t python_boost_helper() { return sizeof(OpticalFresnel::MultiLayerCoeff_t); }
-//};
-
-
 #endif // OPTICALFRESNEL_H
diff --git a/Core/Algorithms/src/LayerDWBASimulation.cpp b/Core/Algorithms/src/LayerDWBASimulation.cpp
index 17f045e25c36201e719962d8aee0591b1b426f8b..01ee05c7d3cc6c8559a8af9940544eb77828b4a3 100644
--- a/Core/Algorithms/src/LayerDWBASimulation.cpp
+++ b/Core/Algorithms/src/LayerDWBASimulation.cpp
@@ -17,18 +17,21 @@ LayerDWBASimulation::~LayerDWBASimulation()
 void LayerDWBASimulation::setKzFunction(
         const IDoubleToComplexFunction& kz_function)
 {
+    delete mp_kz_function;
     mp_kz_function = kz_function.clone();
 }
 
 void LayerDWBASimulation::setTFunction(
         const IDoubleToComplexFunction& T_function)
 {
+    delete mp_T_function;
     mp_T_function = T_function.clone();
 }
 
 void LayerDWBASimulation::setRFunction(
         const IDoubleToComplexFunction& R_function)
 {
+    delete mp_R_function;
     mp_R_function = R_function.clone();
 }
 
diff --git a/Core/Algorithms/src/MultiLayerDWBASimulation.cpp b/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
index bbb18f889c5f81216ccbe17edf01fba4cbcecd08..a23655a12889f883ee60c766054ba6f364b5a069 100644
--- a/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
+++ b/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
@@ -1,9 +1,11 @@
 #include "MultiLayerDWBASimulation.h"
 #include "OpticalFresnel.h"
 #include "DoubleToComplexInterpolatingFunction.h"
+#include "MultiLayerRoughnessDWBASimulation.h"
+
 
 MultiLayerDWBASimulation::MultiLayerDWBASimulation(
-        const MultiLayer* p_multi_layer)
+        const MultiLayer* p_multi_layer) : m_roughness_dwba_simulation(0)
 {
     mp_multi_layer = p_multi_layer->clone();
 }
@@ -15,6 +17,7 @@ MultiLayerDWBASimulation::~MultiLayerDWBASimulation()
     {
         delete (*it).second;
     }
+    delete m_roughness_dwba_simulation;
 }
 
 void MultiLayerDWBASimulation::init(const Experiment& experiment)
@@ -27,6 +30,14 @@ void MultiLayerDWBASimulation::init(const Experiment& experiment)
             p_layer_dwba_sim->init(experiment);
         }
     }
+    // scattering from rough surfaces in DWBA
+    for (size_t i=0; i<mp_multi_layer->getNumberOfInterfaces(); ++i) {
+        if(mp_multi_layer->getLayerInterface(i)->getRoughness() ) {
+            m_roughness_dwba_simulation = new MultiLayerRoughnessDWBASimulation(mp_multi_layer);
+            m_roughness_dwba_simulation->init(experiment);
+            break;
+        }
+    }
 }
 
 void MultiLayerDWBASimulation::run()
@@ -39,7 +50,6 @@ void MultiLayerDWBASimulation::run()
     std::map<double, OpticalFresnel::MultiLayerCoeff_t> fresnel_coeff_map;
     for (size_t i=0; i<p_alpha_axis->getSize(); ++i) {
         double angle = (*p_alpha_axis)[i];
-        //kvector_t kvec = kvector_t::LambdaAlphaPhi(lambda, -angle, 0.0);
         kvector_t kvec;
         kvec.setLambdaAlphaPhi(lambda, -angle, 0.0);
         OpticalFresnel::MultiLayerCoeff_t coeffs;
@@ -50,21 +60,47 @@ void MultiLayerDWBASimulation::run()
     OpticalFresnel::MultiLayerCoeff_t coeffs;
     fresnelCalculator.execute(*mp_multi_layer, m_ki, coeffs);
     fresnel_coeff_map[-m_alpha_i] = coeffs;
-    // Run through layers for DWBA corrections
-    for (std::map<size_t, LayerDWBASimulation*>::const_iterator it=m_layer_dwba_simulation_map.begin();
-            it!=m_layer_dwba_simulation_map.end(); ++it) {
-        size_t i = (*it).first;
-        LayerDWBASimulation *p_layer_dwba_sim = (*it).second;
-        // Construct kz, T and R functions for use inside layer i
+
+//    // Run through layers for DWBA corrections
+//    for (std::map<size_t, LayerDWBASimulation*>::const_iterator it=m_layer_dwba_simulation_map.begin();
+//            it!=m_layer_dwba_simulation_map.end(); ++it) {
+//        size_t i = (*it).first;
+//        LayerDWBASimulation *p_layer_dwba_sim = (*it).second;
+//        // Construct kz, T and R functions for use inside layer i
+//        std::map<double, complex_t> kz_map;
+//        std::map<double, complex_t> T_map;
+//        std::map<double, complex_t> R_map;
+//        for (std::map<double, OpticalFresnel::MultiLayerCoeff_t>::const_iterator it=fresnel_coeff_map.begin();
+//                it!=fresnel_coeff_map.end(); ++it) {
+//            double angle = (*it).first;
+//            complex_t kz = (*it).second[i].kz;
+//            complex_t T = (*it).second[i].T;
+//            complex_t R = (*it).second[i].R;
+//            kz_map[angle] = kz;
+//            T_map[angle] = T;
+//            R_map[angle] = R;
+//        }
+//        DoubleToComplexInterpolatingFunction kz_function(kz_map);
+//        DoubleToComplexInterpolatingFunction T_function(T_map);
+//        DoubleToComplexInterpolatingFunction R_function(R_map);
+//        p_layer_dwba_sim->setKzTAndRFunctions(kz_function, T_function, R_function);
+//        p_layer_dwba_sim->run();
+//        m_dwba_intensity += p_layer_dwba_sim->getDWBAIntensity();
+
+
+//    }
+
+    // run through layers and construct T,R functions
+    for(size_t i_layer=0; i_layer<mp_multi_layer->getNumberOfLayers(); ++i_layer) {
         std::map<double, complex_t> kz_map;
         std::map<double, complex_t> T_map;
         std::map<double, complex_t> R_map;
         for (std::map<double, OpticalFresnel::MultiLayerCoeff_t>::const_iterator it=fresnel_coeff_map.begin();
                 it!=fresnel_coeff_map.end(); ++it) {
             double angle = (*it).first;
-            complex_t kz = (*it).second[i].kz;
-            complex_t T = (*it).second[i].T;
-            complex_t R = (*it).second[i].R;
+            complex_t kz = (*it).second[i_layer].kz;
+            complex_t T = (*it).second[i_layer].T;
+            complex_t R = (*it).second[i_layer].R;
             kz_map[angle] = kz;
             T_map[angle] = T;
             R_map[angle] = R;
@@ -72,10 +108,28 @@ void MultiLayerDWBASimulation::run()
         DoubleToComplexInterpolatingFunction kz_function(kz_map);
         DoubleToComplexInterpolatingFunction T_function(T_map);
         DoubleToComplexInterpolatingFunction R_function(R_map);
-        p_layer_dwba_sim->setKzTAndRFunctions(kz_function, T_function, R_function);
-        p_layer_dwba_sim->run();
-        // TODO: this must be a += statement when this is implemented in OutputData<double>
-        //m_dwba_intensity += (*p_layer_dwba_sim->getDWBAIntensity());
-        m_dwba_intensity += p_layer_dwba_sim->getDWBAIntensity();
+
+        // layer DWBA simulation
+        std::map<size_t, LayerDWBASimulation*>::const_iterator pos = m_layer_dwba_simulation_map.find(i_layer);
+        if(pos != m_layer_dwba_simulation_map.end() ) {
+            LayerDWBASimulation *p_layer_dwba_sim = pos->second;
+            p_layer_dwba_sim->setKzTAndRFunctions(kz_function, T_function, R_function);
+            p_layer_dwba_sim->run();
+            m_dwba_intensity += p_layer_dwba_sim->getDWBAIntensity();
+        }
+
+        if(m_roughness_dwba_simulation) {
+            m_roughness_dwba_simulation->setKzTAndRFunctions(i_layer, kz_function, T_function, R_function);
+        }
+
+    } // i_layer
+
+
+    if(m_roughness_dwba_simulation) {
+        m_roughness_dwba_simulation->run();
+        m_dwba_intensity += m_roughness_dwba_simulation->getDWBAIntensity();
     }
+
 }
+
+
diff --git a/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp b/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1ddc4ff4acf2d76172705ae62b1074880dd96d93
--- /dev/null
+++ b/Core/Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp
@@ -0,0 +1,251 @@
+#include "MultiLayerRoughnessDWBASimulation.h"
+#include "MultiLayer.h"
+
+
+MultiLayerRoughnessDWBASimulation::MultiLayerRoughnessDWBASimulation(const MultiLayer *p_multi_layer)
+{
+    mp_multi_layer = p_multi_layer->clone();
+    mp_kz_function.resize(mp_multi_layer->getNumberOfLayers(), 0);
+    mp_T_function.resize(mp_multi_layer->getNumberOfLayers(), 0);
+    mp_R_function.resize(mp_multi_layer->getNumberOfLayers(), 0);
+    std::cout << "XXX " << std::endl;
+    std::cout << *mp_multi_layer << std::endl;
+    std::cout << "XXX XXX" << std::endl;
+    mp_multi_layer->walk_and_print();
+}
+
+
+MultiLayerRoughnessDWBASimulation::~MultiLayerRoughnessDWBASimulation()
+{
+    for(size_t i=0; i<mp_kz_function.size(); ++i) delete mp_kz_function[i];
+    for(size_t i=0; i<mp_T_function.size(); ++i) delete mp_T_function[i];
+    for(size_t i=0; i<mp_R_function.size(); ++i) delete mp_R_function[i];
+    delete mp_multi_layer;
+}
+
+
+void MultiLayerRoughnessDWBASimulation::setKzTAndRFunctions(int i_layer,
+        const IDoubleToComplexFunction& kz_function,
+        const IDoubleToComplexFunction& T_function,
+        const IDoubleToComplexFunction& R_function)
+{
+    std::cout << " XXX ks functions:" << std::endl;
+    delete mp_kz_function[i_layer];
+    mp_kz_function[i_layer] = kz_function.clone();
+
+    delete mp_T_function[i_layer];
+    mp_T_function[i_layer] = T_function.clone();
+
+    delete mp_R_function[i_layer];
+    mp_R_function[i_layer] = R_function.clone();
+}
+
+
+
+void MultiLayerRoughnessDWBASimulation::run()
+{
+
+    m_dwba_intensity.resetIndex();
+    double lambda = 2.0*M_PI/m_ki.mag();
+
+//    complex_t k_iz = -mp_kz_function->evaluate(-m_alpha_i);
+//    complex_t k_fz = mp_kz_function->evaluate(alpha_f);
+
+    while (m_dwba_intensity.hasNext())
+    {
+        double phi_f = m_dwba_intensity.getCurrentValueOfAxis<double>("phi_f");
+        m_alpha_f = m_dwba_intensity.getCurrentValueOfAxis<double>("alpha_f");
+        m_kf.setLambdaAlphaPhi(lambda, m_alpha_f, phi_f);
+        std::cout << m_alpha_i << " " << m_alpha_f << std::endl;
+        m_dwba_intensity.next() = evaluate(m_ki, m_kf);
+    }
+
+}
+
+
+double MultiLayerRoughnessDWBASimulation::evaluate(kvector_t k_i, kvector_t k_f)
+{
+    kvector_t q = k_f - k_i;
+    double autocorr(0);
+    for(size_t i=0; i<mp_multi_layer->getNumberOfLayers()-1; i++){
+        const LayerRoughness *rough = mp_multi_layer->getLayerBottomInterface(i)->getRoughness();
+        if(rough) {
+            autocorr += std::norm( get_refractive_term(i)) * std::norm(get_sum4terms(i) ) * rough->getSpectralFun(q);
+        } else {
+            std::cout << "MultiLayerRoughnessDWBASimulation::evaluate() -> No roughness for layer " << i << std::endl;
+        }
+    }
+    return autocorr*m_ki.mag2()/16./M_PI;
+
+}
+
+
+complex_t MultiLayerRoughnessDWBASimulation::get_refractive_term(int ilayer)
+{
+    complex_t n1 = mp_multi_layer->getLayer(ilayer)->getRefractiveIndex();
+    complex_t n2 = mp_multi_layer->getLayer(ilayer+1)->getRefractiveIndex();
+    return n1*n1-n2*n2;
+}
+
+
+complex_t MultiLayerRoughnessDWBASimulation::get_sum4terms(int ilayer)
+{
+    double qz1 =  m_ki.z() + m_kf.z();
+    double qz2 =  m_ki.z() - m_kf.z();
+    double qz3 = -m_ki.z() + m_kf.z();
+    double qz4 = -m_ki.z() - m_kf.z();
+
+    std::cout << " ZZZ " << mp_T_function.size() << " " << std::endl;
+    for(size_t i=0; i<mp_T_function.size(); ++i) {
+        std::cout << mp_T_function[i] << std::endl;
+    }
+
+    complex_t Ti = mp_T_function[ilayer+1]->evaluate(m_alpha_i);
+    complex_t Tf = mp_T_function[ilayer+1]->evaluate(m_alpha_f);
+    complex_t Ri = mp_R_function[ilayer+1]->evaluate(m_alpha_i);
+    complex_t Rf = mp_R_function[ilayer+1]->evaluate(m_alpha_f);
+    double sigma2 = -0.5*std::pow(mp_multi_layer->getLayerBottomInterface(ilayer)->getRoughness()->getSigma(), 2);
+    complex_t term1 = Ti * Tf * std::exp( sigma2*qz1*qz1 );
+    complex_t term2 = Ti * Rf * std::exp( sigma2*qz2*qz2 );
+    complex_t term3 = Ri * Tf * std::exp( sigma2*qz3*qz3 );
+    complex_t term4 = Ri * Rf * std::exp( sigma2*qz4*qz4 );
+
+    return term1 + term2 + term3 + term4;
+}
+
+
+
+/*
+
+void DWBADiffuseReflection::execute(const MultiLayer &sample, const kvector_t &ki, const kvector_t &kf)
+{
+    setSample(sample);
+    setKvectors(ki, kf);
+
+    diffuse_autocorr();
+
+    if( m_sample->getCrossCorrLength() != 0) {
+        diffuse_crosscorr();
+    } else {
+        m_diffuse_crosscorr = 0;
+    }
+
+}
+
+
+void DWBADiffuseReflection::setKvectors(const kvector_t &ki, const kvector_t &kf)
+{
+    m_ki = ki;
+    m_q = kf - ki;
+    m_qz1 = std::pow(ki.z() + kf.z(), 2);
+    m_qz2 = std::pow(ki.z() - kf.z(), 2);
+    m_qz3 = std::pow(-ki.z() + kf.z(), 2);
+    m_qz4 = std::pow(-ki.z() - kf.z(), 2);
+    OpticalFresnel calculator;
+    calculator.execute(*m_sample, ki, m_fcoeff_i);
+    calculator.execute(*m_sample, kf, m_fcoeff_f);
+}
+
+
+void DWBADiffuseReflection::diffuse_autocorr()
+{
+    double autocorr(0);
+    for(size_t i=0; i<m_sample->getNumberOfLayers()-1; i++){
+        autocorr += std::norm( get_refractive_term(i)) * std::norm(get_sum4terms(i) ) * m_sample->getLayerBottomInterface(i)->getRoughness()->getSpectralFun(m_q);
+    }
+    m_diffuse_autocorr = autocorr*m_ki.mag2()/16./M_PI;
+}
+
+
+void DWBADiffuseReflection::diffuse_crosscorr()
+{
+    complex_t crosscorr(0);
+
+    std::vector<complex_t > rterm;
+    rterm.resize(m_sample->getNumberOfLayers());
+    std::vector<complex_t > sterm;
+    sterm.resize(m_sample->getNumberOfLayers());
+
+    for(size_t j=0; j<m_sample->getNumberOfLayers()-1; j++){
+        rterm[j] = get_refractive_term(j);
+        sterm[j] = get_sum4terms(j);
+    }
+
+    for(size_t j=0; j<m_sample->getNumberOfLayers()-1; j++){
+        for(size_t k=0; k<m_sample->getNumberOfLayers()-1; k++) {
+            if(j==k) continue;
+            crosscorr += rterm[j]*sterm[j]*rterm[k]*m_sample->getCrossCorrSpectralFun(m_q,j,k)*std::conj(sterm[k]);
+        }
+    }
+    m_diffuse_crosscorr = crosscorr.real()*m_ki.mag2()/16./M_PI;
+}
+
+
+complex_t DWBADiffuseReflection::get_refractive_term(int ilayer)
+{
+    complex_t n1 = m_sample->getLayer(ilayer)->getRefractiveIndex();
+    complex_t n2 = m_sample->getLayer(ilayer+1)->getRefractiveIndex();
+    return n1*n1-n2*n2;
+}
+
+
+complex_t DWBADiffuseReflection::get_sum4terms(int ilayer)
+{
+//    double qz1 = m_ki.z() + m_kf.z();
+//    double qz2 = m_ki.z() - m_kf.z();
+//    double qz3 = -m_ki.z() + m_kf.z();
+//    double qz4 = -m_ki.z() - m_kf.z();
+//    complex_t term1 = m_fcoeff_i[ilayer+1].T * m_fcoeff_f[ilayer+1].T * std::exp( -0.5*std::pow(sigma * qz1 ,2) );
+//    complex_t term2 = m_fcoeff_i[ilayer+1].T * m_fcoeff_f[ilayer+1].R * std::exp( -0.5*std::pow(sigma * qz2 ,2) );
+//    complex_t term3 = m_fcoeff_i[ilayer+1].R * m_fcoeff_f[ilayer+1].T * std::exp( -0.5*std::pow(sigma * qz3 ,2) );
+//    complex_t term4 = m_fcoeff_i[ilayer+1].R * m_fcoeff_f[ilayer+1].R * std::exp( -0.5*std::pow(sigma * qz4 ,2) );
+//
+//    complex_t term1 = m_fcoeff_i[ilayer+1].T * m_fcoeff_f[ilayer+1].T * std::exp( -0.5*std::pow(sigma * m_qz1 ,2) );
+//    complex_t term2 = m_fcoeff_i[ilayer+1].T * m_fcoeff_f[ilayer+1].R * std::exp( -0.5*std::pow(sigma * m_qz2 ,2) );
+//    complex_t term3 = m_fcoeff_i[ilayer+1].R * m_fcoeff_f[ilayer+1].T * std::exp( -0.5*std::pow(sigma * m_qz3 ,2) );
+//    complex_t term4 = m_fcoeff_i[ilayer+1].R * m_fcoeff_f[ilayer+1].R * std::exp( -0.5*std::pow(sigma * m_qz4 ,2) );
+
+    double sigma2 = -0.5*std::pow(m_sample->getLayerBottomInterface(ilayer)->getRoughness()->getSigma(), 2);
+    complex_t term1 = m_fcoeff_i[ilayer+1].T * m_fcoeff_f[ilayer+1].T * std::exp( sigma2*m_qz1 );
+    complex_t term2 = m_fcoeff_i[ilayer+1].T * m_fcoeff_f[ilayer+1].R * std::exp( sigma2*m_qz2 );
+    complex_t term3 = m_fcoeff_i[ilayer+1].R * m_fcoeff_f[ilayer+1].T * std::exp( sigma2*m_qz3 );
+    complex_t term4 = m_fcoeff_i[ilayer+1].R * m_fcoeff_f[ilayer+1].R * std::exp( sigma2*m_qz4 );
+
+    return term1 + term2 + term3 + term4;
+}
+
+*/
+
+//double DWBADiffuseReflection::execute0(const MultiLayer &sample, const kvector_t &ki, const kvector_t &kf)
+//{
+
+//    OpticalFresnel::MultiLayerCoeff_t c_f;
+//    OpticalFresnel::execute(sample, ki, c_f);
+
+//    OpticalFresnel::MultiLayerCoeff_t c_i;
+//    OpticalFresnel::execute(sample, kf, c_i);
+
+//    double autosum(0);
+//    for(size_t i=0; i<sample.getNumberOfLayers()-1; i++){
+//        complex_t n1 = sample.getLayer(i)->getRefractiveIndex();
+//        complex_t n2 = sample.getLayer(i+1)->getRefractiveIndex();
+//        const LayerRoughness &rough = sample.getLayerBottomInterface(i)->getRoughness();
+//        double sigma = rough.getSigma();
+
+//        kvector_t q = ki - kf;
+//        double qz1 = ki.z() + kf.z();
+//        double qz2 = -ki.z() - kf.z();
+//        double qz3 = ki.z() - kf.z();
+//        double qz4 = -ki.z() + kf.z();
+//        complex_t term1 = c_i[i+1].T*c_f[i+1].T * std::exp( -0.5*std::pow(sigma * qz1 ,2) );
+//        complex_t term2 = c_i[i+1].R*c_f[i+1].T * std::exp( -0.5*std::pow(sigma * qz2 ,2) );
+//        complex_t term3 = c_i[i+1].T*c_f[i+1].R * std::exp( -0.5*std::pow(sigma * qz3 ,2) );
+//        complex_t term4 = c_i[i+1].R*c_f[i+1].R * std::exp( -0.5*std::pow(sigma * qz4 ,2) );
+//        double x1 = std::norm( std::pow(n1,2)-std::pow(n2,2) ) * std::norm(term1 + term2 + term3 + term4) * rough.getSpectralFun(q);
+//        autosum += x1;
+//    }
+
+//    const double S = 1.0;
+//    double CS = autosum*S*ki.mag2()/16./M_PI;
+//    return CS;
+//}
diff --git a/Core/Core.pro b/Core/Core.pro
index 0e4c3cc725a3f97f7cb8d9d12af25ed210f0dce9..c2705bdb5787b3cdef627150b9bb7c3f8547c101 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -36,6 +36,7 @@ SOURCES += \
     Algorithms/src/LayerDWBASimulation.cpp \
     Algorithms/src/LocalMonodisperseApproximationStrategy.cpp \
     Algorithms/src/MultiLayerDWBASimulation.cpp \
+    Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp \
     Algorithms/src/OpticalFresnel.cpp \
     \
     Geometry/src/BasicVector3D.cpp \
@@ -121,6 +122,7 @@ HEADERS += \
     Algorithms/inc/LayerDWBASimulation.h \
     Algorithms/inc/LocalMonodisperseApproximationStrategy.h \
     Algorithms/inc/MultiLayerDWBASimulation.h \
+    Algorithms/inc/MultiLayerRoughnessDWBASimulation.h \
     Algorithms/inc/OpticalFresnel.h \
     \
     Geometry/inc/BasicVector3D.h \
diff --git a/Core/Samples/inc/Layer.h b/Core/Samples/inc/Layer.h
index 388c2db90fb7f0ce120f296e038a8adcf73c29bc..8620741a59e63e8f004c03d87a28c6738a49a39e 100644
--- a/Core/Samples/inc/Layer.h
+++ b/Core/Samples/inc/Layer.h
@@ -54,6 +54,9 @@ public:
     //! return refractive index of the layer's material
     virtual complex_t getRefractiveIndex() const { return (dynamic_cast<const HomogeneousMaterial *>(mp_material))->getRefractiveIndex(); }
 
+    //! return false (override is important for polymorphism of LayerDecorator)
+    virtual bool hasDWBASimulation() const { return false; }
+
     //! return zero pointer (override is important for polymorphism of LayerDecorator)
     virtual LayerDWBASimulation *createDWBASimulation() const { return 0; }
 
diff --git a/Core/Samples/inc/LayerDecorator.h b/Core/Samples/inc/LayerDecorator.h
index f492a2dde94f59e980aae4ab40d32cbabdb4b0ad..8da51a71c9407d6b23e19117c6b9e3b2296ddc1e 100644
--- a/Core/Samples/inc/LayerDecorator.h
+++ b/Core/Samples/inc/LayerDecorator.h
@@ -65,6 +65,8 @@ public:
     const NanoParticleDecoration* getDecoration() const { return mp_decoration; }
     void setDecoration(NanoParticleDecoration* mpDecoration) { mp_decoration = mpDecoration; }
 
+    virtual bool hasDWBASimulation() const { return true; }
+
     virtual LayerDecoratorDWBASimulation *createDWBASimulation() const {
         return new LayerDecoratorDWBASimulation(this);
     }
diff --git a/Core/Samples/inc/LayerInterface.h b/Core/Samples/inc/LayerInterface.h
index 0c24318399f6f3752ecfe13be2c03f3ffdee9a76..a78d5b641dbe64a8d9ba7383becd740b1d4cede2 100644
--- a/Core/Samples/inc/LayerInterface.h
+++ b/Core/Samples/inc/LayerInterface.h
@@ -76,4 +76,5 @@ private:
     const Layer *m_LayerBottom;    //!< pointer to the layer below interface
 };
 
+
 #endif // LAYERINTERFACE_H
diff --git a/Core/Samples/inc/MultiLayer.h b/Core/Samples/inc/MultiLayer.h
index b4972643032568d3a37f1046e1ffe001dcd05673..6b50bb2a7bd0fe698f6f6c759526d0ffea5935cc 100644
--- a/Core/Samples/inc/MultiLayer.h
+++ b/Core/Samples/inc/MultiLayer.h
@@ -46,6 +46,9 @@ public:
     //! return number of layers in multilayer
     inline size_t getNumberOfLayers() const { return m_layers.size(); }
 
+    //! return number of interfaces in multilayer
+    inline size_t getNumberOfInterfaces() const { return m_interfaces.size(); }
+
     //! add object to multilayer, overrides from ISample
     void addLayer(const Layer &p_child);
 
@@ -55,6 +58,9 @@ public:
     //! return layer with given index
     inline const Layer *getLayer(size_t i_layer) const { return m_layers[ check_layer_index(i_layer) ]; }
 
+    //! return layer with given index
+    inline const LayerInterface *getLayerInterface(size_t i_interface) const { return m_interfaces[ check_interface_index(i_interface) ]; }
+
     //! return z-coordinate of the layer's bottom
     inline double getLayerBottomZ(size_t i_layer) const { return m_layers_z[ check_layer_index(i_layer) ]; }
 
diff --git a/Core/Samples/src/MultiLayer.cpp b/Core/Samples/src/MultiLayer.cpp
index f80f29763f9c037548da1a99089c4969df2d9624..d6965952c6f3934816795be0423e56f2ed50b9da 100644
--- a/Core/Samples/src/MultiLayer.cpp
+++ b/Core/Samples/src/MultiLayer.cpp
@@ -67,7 +67,12 @@ MultiLayer *MultiLayer::clone() const
         const Layer *topLayer = newMultiLayer->m_layers[i];
         const Layer *bottomLayer = newMultiLayer->m_layers[i+1];
 
-        LayerInterface *newInterface = LayerInterface::createRoughInterface(topLayer, bottomLayer, *m_interfaces[i]->getRoughness() );
+        LayerInterface *newInterface(0);
+        if(m_interfaces[i]->getRoughness()) {
+            newInterface = LayerInterface::createRoughInterface(topLayer, bottomLayer, *m_interfaces[i]->getRoughness() );
+        } else {
+            newInterface = LayerInterface::createSmoothInterface(topLayer, bottomLayer );
+        }
         newMultiLayer->addAndRegisterInterface( newInterface );
     }
 
@@ -108,7 +113,12 @@ void MultiLayer::addLayerWithTopRoughness(const Layer &layer, const LayerRoughne
     if (getNumberOfLayers())
     {
         const Layer *p_last_layer = m_layers.back();
-        LayerInterface *interface = LayerInterface::createRoughInterface( p_last_layer, p_new_layer, roughness);
+        LayerInterface *interface(0);
+        if(roughness.getSigma() != 0) {
+            interface = LayerInterface::createRoughInterface( p_last_layer, p_new_layer, roughness);
+        } else {
+            interface = LayerInterface::createSmoothInterface( p_last_layer, p_new_layer);
+        }
         addAndRegisterLayer(p_new_layer);
         addAndRegisterInterface(interface);
         m_layers_z.push_back(m_layers_z.back() - layer.getThickness() );
@@ -123,9 +133,8 @@ void MultiLayer::addLayerWithTopRoughness(const Layer &layer, const LayerRoughne
 /* ************************************************************************* */
 // add layer with default (zero) roughness
 /* ************************************************************************* */
-void MultiLayer::addLayer(const Layer &p_child)
+void MultiLayer::addLayer(const Layer &layer)
 {
-    const Layer &layer = dynamic_cast<const Layer &>(p_child);
     addLayerWithTopRoughness(layer, LayerRoughness());
 }
 
@@ -153,15 +162,15 @@ void MultiLayer::addLayer(const Layer &p_child)
 /* ************************************************************************* */
 double MultiLayer::getCrossCorrSpectralFun(const kvector_t &kvec, int j, int k) const
 {
+    if(m_crossCorrLength == 0) return 0.0;
     double z_j = getLayerBottomZ(j);
     double z_k = getLayerBottomZ(k);
     const LayerRoughness *rough_j = getLayerBottomInterface(j)->getRoughness();
     const LayerRoughness *rough_k = getLayerBottomInterface(k)->getRoughness();
+    if( !rough_j || !rough_k ) return 0.0;
     double sigma_j = rough_j->getSigma();
     double sigma_k = rough_k->getSigma();
-    if(sigma_j == 0 || sigma_k ==0 || m_crossCorrLength==0) {
-        return 0.0;
-    }
+    if(sigma_j == 0 || sigma_k == 0) return 0.0;
     double corr = 0.5*( (sigma_k/sigma_j)*rough_j->getSpectralFun(kvec) + (sigma_j/sigma_k)*rough_k->getSpectralFun(kvec) ) * std::exp( -1*std::abs(z_j-z_k)/m_crossCorrLength );
     return corr;
 }
@@ -188,12 +197,18 @@ void MultiLayer::setLayerThickness(size_t i_layer, double thickness)
 MultiLayerDWBASimulation* MultiLayer::createDWBASimulation() const
 {
     for (size_t i=0; i<getNumberOfLayers(); ++i) {
-        LayerDWBASimulation *p_layer_dwba_sim = getLayer(i)->createDWBASimulation();
-        if ( p_layer_dwba_sim ) {
-            delete p_layer_dwba_sim;
-            return new MultiLayerDWBASimulation(this);
-        }
+        if( getLayer(i)->hasDWBASimulation() ) return new MultiLayerDWBASimulation(this);
+//        LayerDWBASimulation *p_layer_dwba_sim = getLayer(i)->createDWBASimulation();
+//        if ( p_layer_dwba_sim ) {
+//            delete p_layer_dwba_sim;
+//            return new MultiLayerDWBASimulation(this);
+//        }
+    }
+
+    for(size_t i=0; i<getNumberOfInterfaces(); ++i) {
+        if( getLayerInterface(i)->getRoughness() ) return new MultiLayerDWBASimulation(this);
     }
+
     return 0;
 }
 
diff --git a/Core/Tools/inc/Exceptions.h b/Core/Tools/inc/Exceptions.h
index 59f46b2faabd708e43a6e11e9f4eb8fefe584fe8..ed5f5b2d32f2d05977b3e84c95ea72c5eda3cffc 100644
--- a/Core/Tools/inc/Exceptions.h
+++ b/Core/Tools/inc/Exceptions.h
@@ -18,85 +18,90 @@
 #include <string>
 
 
-class NotImplementedException : public std::logic_error
-{
-public:
-    NotImplementedException(const std::string &message);
-};
-
-class NullPointerException : public std::logic_error
-{
-public:
-    NullPointerException(const std::string& message);
-};
-
-class OutOfBoundsException : public std::logic_error
-{
-public:
-    OutOfBoundsException(const std::string& message);
-};
-
-class ClassInitializationException : public std::runtime_error
-{
-public:
-    ClassInitializationException(const std::string& message);
-};
-
-class SelfReferenceException : public std::logic_error
-{
-public:
-    SelfReferenceException(const std::string &message);
-};
-
-class DeadReferenceException : public std::runtime_error
-{
-public:
-    DeadReferenceException(const std::string& message);
-};
-
-class UnknownClassRegistrationException : public std::runtime_error
-{
-public:
-    UnknownClassRegistrationException(const std::string& message);
-};
-
-class ExistingClassRegistrationException : public std::runtime_error
-{
-public:
-    ExistingClassRegistrationException(const std::string& message);
-};
-
-class LogicErrorException : public std::logic_error
-{
-public:
-    LogicErrorException(const std::string& message);
-};
-
-class RuntimeErrorException : public std::runtime_error
-{
-public:
-    RuntimeErrorException(const std::string& message);
-};
-
-class DivisionByZeroException : public std::runtime_error
-{
-public:
-    DivisionByZeroException(const std::string& message);
-};
-
-class DomainErrorException : public std::domain_error
-{
-public:
-    DomainErrorException(const std::string& message);
-};
-
-class FileNotIsOpenException : public std::runtime_error
-{
-public:
-    FileNotIsOpenException(const std::string& message);
-};
-
-void LogExceptionMessage(const std::string &message);
-
+namespace Exceptions {
+
+    class NotImplementedException : public std::logic_error
+    {
+    public:
+        NotImplementedException(const std::string &message);
+    };
+
+    class NullPointerException : public std::logic_error
+    {
+    public:
+        NullPointerException(const std::string& message);
+    };
+
+    class OutOfBoundsException : public std::logic_error
+    {
+    public:
+        OutOfBoundsException(const std::string& message);
+    };
+
+    class ClassInitializationException : public std::runtime_error
+    {
+    public:
+        ClassInitializationException(const std::string& message);
+    };
+
+    class SelfReferenceException : public std::logic_error
+    {
+    public:
+        SelfReferenceException(const std::string &message);
+    };
+
+    class DeadReferenceException : public std::runtime_error
+    {
+    public:
+        DeadReferenceException(const std::string& message);
+    };
+
+    class UnknownClassRegistrationException : public std::runtime_error
+    {
+    public:
+        UnknownClassRegistrationException(const std::string& message);
+    };
+
+    class ExistingClassRegistrationException : public std::runtime_error
+    {
+    public:
+        ExistingClassRegistrationException(const std::string& message);
+    };
+
+    class LogicErrorException : public std::logic_error
+    {
+    public:
+        LogicErrorException(const std::string& message);
+    };
+
+    class RuntimeErrorException : public std::runtime_error
+    {
+    public:
+        RuntimeErrorException(const std::string& message);
+    };
+
+    class DivisionByZeroException : public std::runtime_error
+    {
+    public:
+        DivisionByZeroException(const std::string& message);
+    };
+
+    class DomainErrorException : public std::domain_error
+    {
+    public:
+        DomainErrorException(const std::string& message);
+    };
+
+    class FileNotIsOpenException : public std::runtime_error
+    {
+    public:
+        FileNotIsOpenException(const std::string& message);
+    };
+
+    void LogExceptionMessage(const std::string &message);
+
+}
+
+using namespace Exceptions;
 
 #endif // EXCEPTIONS_H
diff --git a/Core/Tools/src/Exceptions.cpp b/Core/Tools/src/Exceptions.cpp
index 8fca8756a8f9d20a5f13a6da901ec3021b8fdb49..2fb30024c07d40cf587c2f64e330d790c6f541ec 100644
--- a/Core/Tools/src/Exceptions.cpp
+++ b/Core/Tools/src/Exceptions.cpp
@@ -1,6 +1,8 @@
 #include "Exceptions.h"
 #include <iostream>
 
+namespace Exceptions {
+
 void LogExceptionMessage(const std::string &message)
 {
     std::cout << message << std::endl;
@@ -84,5 +86,5 @@ FileNotIsOpenException::FileNotIsOpenException(const std::string &message)
     LogExceptionMessage(message);
 }
 
-
+}
 
diff --git a/Doc/Doxygen/Doxyfile b/Doc/Doxygen/Doxyfile
index f8357d017d67dd79bb65e849a5d7d281a1d18125..31fd7034bd17b4f1f6a6153bc4ac15c8ce574118 100644
--- a/Doc/Doxygen/Doxyfile
+++ b/Doc/Doxygen/Doxyfile
@@ -38,7 +38,7 @@ PROJECT_NUMBER         =
 # for a project that appears at the top of each page and should give viewer
 # a quick idea about the purpose of the project. Keep the description short.
 
-PROJECT_BRIEF          =
+PROJECT_BRIEF          = "Simulation of neutron and x-ray scattering at grazing incidence"
 
 # With the PROJECT_LOGO tag one can specify an logo or icon that is
 # included in the documentation. The maximum height of the logo should not
@@ -395,7 +395,7 @@ HIDE_IN_BODY_DOCS      = NO
 # to NO (the default) then the documentation will be excluded.
 # Set it to YES to include the internal documentation.
 
-INTERNAL_DOCS          = NO
+INTERNAL_DOCS          = YES
 
 # If the CASE_SENSE_NAMES tag is set to NO then Doxygen will only generate
 # file names in lower-case letters. If set to YES upper-case letters are also
@@ -674,7 +674,7 @@ EXCLUDE_SYMLINKS       = NO
 # against the file with absolute path, so to exclude all test directories
 # for example use the pattern */test/*
 
-EXCLUDE_PATTERNS       =
+EXCLUDE_PATTERNS       = */PythonAPI/*
 
 # The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names
 # (namespaces, classes, functions, etc.) that should be excluded from the
diff --git a/Examples/Performance/perf_history.txt b/Examples/Performance/perf_history.txt
index 8c2d403f2ccf46b61cd3f4754fabf5609d749731..482f239ba3de8e1b9945afff7559257b99a335d7 100644
--- a/Examples/Performance/perf_history.txt
+++ b/Examples/Performance/perf_history.txt
@@ -25,11 +25,9 @@
 2012-07-30 16:19:46 | jcnsopc73  | macosx64, 2800 MHz      | 81967.2 | 3.8835  | 3.83142 | 0.19175 | 
 
 # multiple changes to get parameters from MesoCrystal, LayerDecorator etc into parameter pool
-2012-07-31 11:21:26 | jcnsopc73  | macosx64, 2800 MHz      | 80321.3 | 3.96825 | 3.94477 | 0.34305 | 
-2012-07-31 11:25:57 | jcnsopc73  | macosx64, 2800 MHz      | 82304.5 | 3.89864 | 3.861   | 0.33444 | 
-2012-07-31 11:27:52 | jcnsopc73  | macosx64, 2800 MHz      | 82987.5 | 3.98406 | 3.861   | 0.34188 | 
-
-# switching off debug
+# switching off debug, it looks like MesoCrystal is now slower (0.45Hz vs 0.70 Hz), origin is not understood
 2012-07-31 11:39:09 | jcnsopc73  | macosx64, 2800 MHz      | 263158  | 13.1579 | 13.0719 | 0.45351 | 
 2012-07-31 11:39:23 | jcnsopc73  | macosx64, 2800 MHz      | 266667  | 12.9032 | 13.0719 | 0.45351 | 
 2012-07-31 11:40:12 | jcnsopc73  | macosx64, 2800 MHz      | 270270  | 13.0719 | 12.9032 | 0.45351 | 
+
+2012-08-03 09:21:48 | jcnsopc73  | macosx64, 2800 MHz      | 77519.4 | 4.04858 | 3.89864 | 0.17574 |