diff --git a/App/inc/FitSuiteHelper.h b/App/inc/FitSuiteHelper.h
index c807db9ad6416abc4e1302e4a5d5353678e64529..004893497ee3b2d669fe81799e904e38e330a3e8 100644
--- a/App/inc/FitSuiteHelper.h
+++ b/App/inc/FitSuiteHelper.h
@@ -16,6 +16,8 @@
#include "IObserver.h"
+#include "OutputData.h"
+#include "ChiSquaredModule.h"
#include <string>
@@ -40,6 +42,9 @@ class FitSuiteObserverDraw : public IObserver
public:
FitSuiteObserverDraw( const std::string &canvas_name = std::string("FitSuiteObserverDraw_c1") ) : m_canvas_name(canvas_name) {}
void update(IObservable *subject);
+
+ //! return output data which contains chi2 values from ChisSquaredModule of FitSuite
+ OutputData<double > *getDifferenceMap(const ChiSquaredModule *chi_module);
private:
std::string m_canvas_name; //! canvas name were to draw
};
diff --git a/App/inc/TestFittingModule.h b/App/inc/TestFittingModule.h
index 215a75c324b9fa58b7fef129504390c4d9b21551..a41d7b02d4fca29ef7525d5f8d1ed4b74d152e82 100644
--- a/App/inc/TestFittingModule.h
+++ b/App/inc/TestFittingModule.h
@@ -10,7 +10,7 @@
// * mollis quis. Mauris commodo rhoncus porttitor. *
// ********************************************************************
//! @file TestFittingModule.h
-//! @brief Definition of
+//! @brief Definition of TestFittingModule class
//! @author Scientific Computing Group at FRM II
//! @date Jul 20, 2012
@@ -21,6 +21,7 @@
#include "FitMultiParameter.h"
#include "GISASExperiment.h"
+
class TestFittingModule : public IFunctionalTest
{
public:
diff --git a/App/src/FitSuiteHelper.cpp b/App/src/FitSuiteHelper.cpp
index a2416de7bfd28abd1d8587d7ddc77f077b0d1d91..84a5e0fa99eb7abdb6d86aeade9f5ed6d67f4c37 100644
--- a/App/src/FitSuiteHelper.cpp
+++ b/App/src/FitSuiteHelper.cpp
@@ -74,16 +74,20 @@ void FitSuiteObserverDraw::update(IObservable *subject)
gPad->SetLogz();
gPad->SetLeftMargin(0.12);
gPad->SetRightMargin(0.12);
- IsGISAXSTools::setMinimum(10.);
+// IsGISAXSTools::setMinimum(10.);
IsGISAXSTools::drawOutputDataRelativeDifference2D(*fitSuite->getChiSquaredModule()->getSimulationData(), *fitSuite->getChiSquaredModule()->getRealData(), "COLZ", "relative difference");
// chi2 difference
c1->cd(4);
gPad->SetLogz();
gPad->SetLeftMargin(0.12);
+ OutputData<double > *diff_map = getDifferenceMap(fitSuite->getChiSquaredModule());
gPad->SetRightMargin(0.12);
- IsGISAXSTools::setMinimum(0.0000001);
- IsGISAXSTools::setMaximum(1.0);
- IsGISAXSTools::drawOutputDataChi2Difference2D(*fitSuite->getChiSquaredModule()->getSimulationData(), *fitSuite->getChiSquaredModule()->getRealData(), "COLZ", "relative chi2 difference");
+// IsGISAXSTools::setMinimum(0.0000001);
+// IsGISAXSTools::resetMinimumAndMaximum();
+// IsGISAXSTools::setMaximum(1.0);
+// IsGISAXSTools::setMinimum(0.000001);
+ IsGISAXSTools::drawOutputDataInPad(*diff_map, "COLZ", "chi2 difference map");
+ delete diff_map;
c1->cd(5);
TPaveText *pt = new TPaveText(.05,.1,.95,.8);
@@ -101,6 +105,38 @@ void FitSuiteObserverDraw::update(IObservable *subject)
}
+// return output data which contains chi2 values from ChisSquaredModule of FitSuite
+OutputData<double > *FitSuiteObserverDraw::getDifferenceMap(const ChiSquaredModule *chi_module)
+{
+ const ISquaredFunction *squared_function = chi_module->getSquaredFunction();
+
+ const OutputData<double> *simu_data = chi_module->getSimulationData();
+ const OutputData<double> *real_data = chi_module->getRealData();
+
+ OutputData<double > *difference = simu_data->clone();
+ difference->setAllTo(0.0);
+
+ double norm(0);
+
+ simu_data->resetIndex();
+ real_data->resetIndex();
+ difference->resetIndex();
+ while (real_data->hasNext()) {
+ double value_simu = simu_data->currentValue();
+ double value_real = real_data->currentValue();
+ double squared_difference = squared_function->calculateSquaredDifference(value_real, value_simu);
+ difference->next() = squared_difference;
+ norm += squared_difference;
+ real_data->next(); simu_data->next();
+ }
+ //norm *= real_data->getAllocatedSize();
+ difference->scaleAll(1./norm);
+// std::cout << "XXX " << aaa/real_data->getAllocatedSize() << " " << chi_module->getValue() << std::endl;
+
+ return difference;
+}
+
+
/* ************************************************************************* */
// Save results of each fit iteration to the disk in the form of ROOT tree
@@ -112,12 +148,6 @@ void FitSuiteObserverWriteTree::update(IObservable *subject)
FitSuite *fitSuite = dynamic_cast<FitSuite *>(subject);
if( !fitSuite ) throw NullPointerException("FitSuiteObserverWriteTree::update() -> Error! Can't cast FitSuite");
- // printing parameter values
- for(FitSuite::fitparameters_t::iterator it = fitSuite->fitparams_begin(); it!=fitSuite->fitparams_end(); ++it) {
- std::cout << " " << (*it)->getName() << " " << (*it)->getValue();
- // std::cout << *(*it);
- }
- std::cout << std::endl;
// preparing root file for writing
// if it is first call the file will be opened in 'recreate' mode, otherwise in 'update' mode
TFile *top(0);
diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index ddedfbf3ebb9e3f517a3a18935ef7b831e226996..d10a7d9cc816c9a142939ff56545705a15214788 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -97,8 +97,8 @@ void IsGISAXSTools::drawOutputDataInPad(const OutputData<double>& output,
TH2D *h2 = IsGISAXSTools::getOutputDataTH2D(output, "p_hist2D");
gStyle->SetPalette(1);
gStyle->SetOptStat(0);
- gPad->SetRightMargin(0.115);
- gPad->SetLeftMargin(0.115);
+// gPad->SetRightMargin(0.115);
+// gPad->SetLeftMargin(0.115);
if( hasMinimum() ) h2->SetMinimum(m_hist_min);
if( hasMaximum() ) h2->SetMaximum(m_hist_max);
h2->SetTitle(histogram_title.c_str());
diff --git a/App/src/ROOTMinimizer.cpp b/App/src/ROOTMinimizer.cpp
index 0970467a40646da4b980c388b691d9f802c3a414..fc7a13eb556c107a674530878f9617cfa728d909 100644
--- a/App/src/ROOTMinimizer.cpp
+++ b/App/src/ROOTMinimizer.cpp
@@ -64,7 +64,7 @@ void ROOTMinimizer::printResults() const
for(size_t i=0; i<getNumberOfVariables(); ++i) {
std::cout << " #" << i
- << " '" << Utils::AdjustStringLength(m_root_minimizer->VariableName(i), 30) << "' "
+ << " " << Utils::AdjustStringLength(m_root_minimizer->VariableName(i), 30) << " "
<< " value:" << getValueOfVariableAtMinimum(i)
<< " error:" << getErrorOfVariable(i) << std::endl;
}
diff --git a/App/src/TestFittingModule.cpp b/App/src/TestFittingModule.cpp
index 4d1b8b620a6cd2f4ccc6d42e81e1ca5e1b7dbfc4..e9afbd357b32e6312dd2ff673b66b87209af1de7 100644
--- a/App/src/TestFittingModule.cpp
+++ b/App/src/TestFittingModule.cpp
@@ -13,6 +13,7 @@
#include "Exceptions.h"
#include "DrawHelper.h"
#include "FitSuiteHelper.h"
+#include "ResolutionFunction2DSimple.h"
#include "IObserver.h"
#include "FitSuite.h"
@@ -57,18 +58,27 @@ void TestFittingModule::execute()
TCanvas *c1 = new TCanvas(canvas_name.c_str(), "Test of the fitting suite", 800, 600);
c1->cd(); gPad->SetLogz();
IsGISAXSTools::drawOutputDataInPad(*mp_exact_data, "CONT4 Z", "exact data");
+ c1->Update();
\
// setting fitSuite
FitSuite *fitSuite = new FitSuite();
fitSuite->setExperiment(mp_experiment);
fitSuite->setRealData(*mp_real_data);
- fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
+ //fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
+
+// ROOT::Math::Minimizer *m_minim = (dynamic_cast<ROOTMinimizer *>(fitSuite->getMinimizer()))->getROOTMinimizer();
+// std::cout << m_minim->SetTolerance(0.00001) << std::endl;
+
+ //fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiMin", "ConjugateFR") );
+ //fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiMin", "BFGS") );
+ //fitSuite->setMinimizer( new ROOTMinimizer("GSLSimAn", "") );
+
fitSuite->addFitParameter("*/MultiLayer/Layer0/thickness", 12*Units::nanometer, 2*Units::nanometer, TRange<double>(1.0, 20.0) );
fitSuite->addFitParameter("*/FormFactorCylinder/radius", 2*Units::nanometer, 2*Units::nanometer, TRange<double>(1.0, 20.0) );
fitSuite->attachObserver( new FitSuiteObserverPrint() );
fitSuite->attachObserver( new FitSuiteObserverDraw() );
- fitSuite->attachObserver( new FitSuiteObserverWriteTree() );
+ //fitSuite->attachObserver( new FitSuiteObserverWriteTree() );
fitSuite->runFit();
@@ -115,6 +125,7 @@ void TestFittingModule::initializeExperiment()
mp_experiment->setSample(*mp_sample);
mp_experiment->setDetectorParameters(100, 0.0*Units::degree, 2.0*Units::degree,100 , 0.0*Units::degree, 2.0*Units::degree);
mp_experiment->setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
+ mp_experiment->setDetectorResolutionFunction(new ResolutionFunction2DSimple(0.0002, 0.0002));
mp_experiment->setBeamIntensity(1e10);
}
@@ -169,7 +180,7 @@ void TestFittingModule::generateRealData(double noise_factor)
mp_real_data->resetIndex();
while (mp_real_data->hasNext()) {
double current = mp_real_data->currentValue();
- double sigma = noise_factor*current;
+ 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;
diff --git a/Core/Algorithms/inc/ChiSquaredModule.h b/Core/Algorithms/inc/ChiSquaredModule.h
index 586f9e28f165ba142554904e58d7ba1a0e8f9cf4..36f8f400d87c3e5f6590201dd548f77e5b0c13af 100644
--- a/Core/Algorithms/inc/ChiSquaredModule.h
+++ b/Core/Algorithms/inc/ChiSquaredModule.h
@@ -43,6 +43,9 @@ public:
//! return chi2 value (should be called after calculateChiSquared)
virtual double getValue() const { return m_chi2_value; }
+ //! return squared function
+ const ISquaredFunction *getSquaredFunction() const { return mp_squared_function; }
+
protected:
OutputData<double> *mp_real_data;
OutputData<double> *mp_simulation_data;
diff --git a/Core/Algorithms/inc/ISquaredFunction.h b/Core/Algorithms/inc/ISquaredFunction.h
index fb97e0206e0bb347f3860b181f419aebd26bef1a..0a466a7411177d4ab5b70725a493bc13f2fdcc48 100644
--- a/Core/Algorithms/inc/ISquaredFunction.h
+++ b/Core/Algorithms/inc/ISquaredFunction.h
@@ -17,6 +17,7 @@
#include "Numeric.h"
#include <cmath>
+#include <iostream>
class ISquaredFunction
{
@@ -41,9 +42,12 @@ inline double DefaultSquaredFunction::calculateSquaredDifference(
double real_value, double simulated_value) const
{
double diff_squared = (simulated_value-real_value)*(simulated_value-real_value);
- if (diff_squared < Numeric::double_epsilon) return 0.0;
- //double normalization = std::max(std::abs(real_value), 1.0);
- double normalization = 1.0;
+ if (diff_squared < Numeric::double_epsilon) {
+ return 0.0;
+ }
+ double normalization = std::max(std::abs(real_value), 1.0);
+// double normalization = real_value;
+ normalization = 1.0;
return diff_squared/normalization;
}