diff --git a/App/inc/ROOTMinimizer.h b/App/inc/ROOTMinimizer.h
index 5878c288e7f41d08e9a132bb41a04923a2652935..79d1b5de1501433146b9db3c62fcd15884340aa9 100644
--- a/App/inc/ROOTMinimizer.h
+++ b/App/inc/ROOTMinimizer.h
@@ -71,6 +71,7 @@ public:
private:
ROOT::Math::Minimizer *m_root_minimizer;
ROOT::Math::Functor *m_fcn;
+ //ROOT::Math::GradFunctor *m_fcn;
};
#endif // ROOTMINIMIZER_H
diff --git a/App/inc/TestIsGISAXS12.h b/App/inc/TestIsGISAXS12.h
index 9faf27bb1f60099d3b4a1058430947eadb302fed..dbb1d2a62f1f0eea2c3dbf34c7f346964806afb3 100644
--- a/App/inc/TestIsGISAXS12.h
+++ b/App/inc/TestIsGISAXS12.h
@@ -35,7 +35,10 @@ public:
virtual void execute();
private:
- typedef std::vector<OutputData<double > *> DataScan_t;
+// typedef std::vector<OutputData<double > *> DataScan_t;
+// enum isgisaxs_scans {kFixedAlphaf, kFixedPhif };
+
+
//! builds sample for fitter testing
class TestSampleBuilder : public ISampleBuilder
@@ -48,17 +51,49 @@ private:
protected:
//! initialize pool parameters, i.e. register some of class members for later access via parameter pool
virtual void init_parameters();
- double m_particle_probability;
+ double m_particle_probability1;
double m_particle_radius1;
- double m_particle_radius2;
double m_dispersion_radius1;
- double m_dispersion_radius2;
double m_height_aspect_ratio1;
+ double m_particle_probability2;
+ double m_particle_radius2;
+ double m_dispersion_radius2;
double m_height_aspect_ratio2;
double m_interf_distance;
double m_interf_width;
};
+ //! represent collection of OutputData object's
+ class DataSet {
+ public:
+ typedef std::vector<OutputData<double > *> dataset_t;
+ typedef std::vector<OutputData<double > *>::iterator iterator;
+ typedef std::vector<OutputData<double > *>::const_iterator const_iterator;
+ DataSet(){}
+ ~DataSet() { clear(); }
+ void clear() {for(size_t i=0; i<m_data.size(); ++i) delete m_data[i]; m_data.clear(); }
+ DataSet(const DataSet &other)
+ {
+ for(size_t i=0; i<other.m_data.size(); ++i) m_data.push_back(other.m_data[i]->clone());
+ }
+ DataSet &operator=(const DataSet &other)
+ {
+ if( this != &other ) {
+ clear();
+ for(size_t i=0; i<other.m_data.size(); ++i) m_data.push_back(other.m_data[i]->clone());
+ }
+ return *this;
+ }
+ void push_back(OutputData<double > *data) { m_data.push_back(data); }
+ dataset_t::iterator begin() { return m_data.begin(); }
+ dataset_t::iterator end() { return m_data.end(); }
+ size_t size() { return m_data.size(); }
+ const OutputData<double> *operator[](int index) const { return m_data[index]; }
+ OutputData<double> *operator[](int index) { return m_data[index]; }
+
+ dataset_t m_data;
+ };
+
//! represent single line stored in isgisaxs *.dat file with data to fit
class IsgiData {
public:
@@ -76,9 +111,6 @@ private:
//! run standard isgisaxs comparison for the given sample
void run_isgisaxs_comparison();
- //! run test fit
- void run_test_fit();
-
//! run isgisaxs ex-12 style fit
void run_isgisaxs_fit();
@@ -86,17 +118,29 @@ private:
void plot_isgisaxs_data();
//! read special isgisaxs *.dat file with data to fit
- void read_isgisaxs_datfile(const std::string &filename);
+ void read_isgisaxs_datfile(const std::string &filename, DataSet &dataset);
+
+ //! read special isgisaxs *.out file with isgisaxs results
+ void read_isgisaxs_outfile(const std::string &filename, DataSet &dataset, bool read_fit_results = false);
//! convert isgisaxs 1d scan to output data 2d object
OutputData<double> *convert_isgi_scan(std::vector<IsgiData > &isgi_data);
+ //! run test minimizer to check whole chain
+ void run_test_minimizer();
+
+ //! run chi module test on isgisaxs data/result pair to check module numericaly
+ void run_test_chimodule();
+
+ void print_axes(DataSet &data);
+
+
std::string m_data_path;
GISASExperiment *m_experiment;
ISampleBuilder *m_sample_builder;
FitSuite *m_fitSuite;
- DataScan_t m_isgi_scans; //! vector of OutputData's representing isgisaxs scans
+// DataScan_t m_isgi_scans; //! vector of OutputData's representing isgisaxs scans
double m_isgi_fixed_alphaf;
double m_isgi_fixed_phif;
};
diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index 3c81357910dc039cdd75a51fce1618326d780e34..cdb08ccc697f26f5766f343be8f9dc59e5df4172 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -299,7 +299,7 @@ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, c
h1_spect.Fill( x );
}
- gPad->SetLogy();
+// gPad->SetLogy();
gPad->SetRightMargin(0.115);
gPad->SetLeftMargin(0.115);
h1_spect.SetStats(1);
diff --git a/App/src/ROOTMinimizer.cpp b/App/src/ROOTMinimizer.cpp
index 829726f6304cdaa2be186fd4ba4650a99c62ff31..8f57e6468091d8ffe1bea488492fd72b3b000e98 100644
--- a/App/src/ROOTMinimizer.cpp
+++ b/App/src/ROOTMinimizer.cpp
@@ -61,6 +61,7 @@ void ROOTMinimizer::minimize()
void ROOTMinimizer::setFunction(boost::function<double(const double *)> fcn, int ndim)
{
m_fcn = new ROOT::Math::Functor(fcn, ndim);
+// m_fcn = new ROOT::Math::GradFunctor(fcn, ndim);
m_root_minimizer->SetFunction(*m_fcn);
}
diff --git a/App/src/TestDetectorResolution.cpp b/App/src/TestDetectorResolution.cpp
index f7f3a5382c7947288843ab77cb2648937fc0b604..a0e7129b66c2efdb8a875111434da103cabc8b95 100644
--- a/App/src/TestDetectorResolution.cpp
+++ b/App/src/TestDetectorResolution.cpp
@@ -13,13 +13,13 @@
namespace
{
-double testResolutionFunction(double u, double v)
-{
- double sigma_u = 0.001;
- double sigma_v = 0.001;
- return MathFunctions::IntegratedGaussian(u, 0.0, sigma_u)
- * MathFunctions::IntegratedGaussian(v, 0.0, sigma_v);
-}
+//double testResolutionFunction(double u, double v)
+//{
+// double sigma_u = 0.001;
+// double sigma_v = 0.001;
+// return MathFunctions::IntegratedGaussian(u, 0.0, sigma_u)
+// * MathFunctions::IntegratedGaussian(v, 0.0, sigma_v);
+//}
}
TestDetectorResolution::TestDetectorResolution()
diff --git a/App/src/TestFittingModule1.cpp b/App/src/TestFittingModule1.cpp
index b2ad0a236d2cab90b87af14c5259510ba090a6f6..12b05872d642e101a3fd17b7dd93c217539e9294 100644
--- a/App/src/TestFittingModule1.cpp
+++ b/App/src/TestFittingModule1.cpp
@@ -186,7 +186,7 @@ void TestFittingModule1::initializeRealData()
if( !mp_experiment ) throw NullPointerException("TestFittingModule2::initializeRealData() -> Error! No experiment o sample defined ");
mp_experiment->runSimulation();
- mp_experiment->normalize();
+ //mp_experiment->normalize();
delete mp_real_data;
mp_real_data = IsGISAXSTools::createNoisyData(*mp_experiment->getOutputData());
}
diff --git a/App/src/TestFittingModule2.cpp b/App/src/TestFittingModule2.cpp
index e761bb4341e0a846ac83e46c1d66a7c34344a3e1..627a9d88578c7d7b1e60176865b9780ba211dc81 100644
--- a/App/src/TestFittingModule2.cpp
+++ b/App/src/TestFittingModule2.cpp
@@ -101,21 +101,26 @@ void TestFittingModule2::fit_example_chimodule()
//mp_experiment->setDetectorResolutionFunction(new ResolutionFunction2DSimple(0.0002, 0.0002));
initializeRealData();
+// m_fitSuite->addFitParameter("*SampleBuilder/m_cylinder_height", 12*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+// m_fitSuite->addFitParameter("*SampleBuilder/m_cylinder_radius", 2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+// m_fitSuite->addFitParameter("*SampleBuilder/m_prism3_half_side", 12*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+// m_fitSuite->addFitParameter("*SampleBuilder/m_prism3_height", 2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/m_cylinder_height", 12*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/m_cylinder_radius", 2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/m_prism3_half_side", 12*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/m_prism3_height", 2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*Normalizer/scale", 1e10, 1e3, AttLimits::limited(1e8, 1e12));
+ m_fitSuite->addFitParameter("*SampleBuilder/m_cylinder_height", 12*Units::nanometer, 0.01*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/m_cylinder_radius", 2*Units::nanometer, 0.01*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/m_prism3_half_side", 12*Units::nanometer, 0.01*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/m_prism3_height", 2*Units::nanometer, 0.01*Units::nanometer, AttLimits::lowerLimited(0.01) );
+// m_fitSuite->addFitParameter("*Normalizer/scale", 1e10, 10., AttLimits::limited(1e9, 2*1e10));
// setting up fitSuite
ChiSquaredModule chiModule;
//chiModule.setChiSquaredFunction( SquaredFunctionDefault() );
+// chiModule.setChiSquaredFunction( SquaredFunctionWhichOnlyWorks() ); // it works only with resolution function, without it fit doesn't converge
chiModule.setChiSquaredFunction( SquaredFunctionWithSystematicError() );
- //chiModule.setChiSquaredFunction( SquaredFunctionWhichOnlyWorks() ); // it works only with resolution function, without it fit doesn't converge
- //chiModule.setOutputDataNormalizer( OutputDataNormalizerScaleAndShift(1e10,0) );
- chiModule.setIntensityFunction( IntensityFunctionLog() );
+// chiModule.setOutputDataNormalizer( OutputDataNormalizerScaleAndShift(1e10,0) );
+ //chiModule.setIntensityFunction( IntensityFunctionLog() );
m_fitSuite->addExperimentAndRealData(*mp_experiment, *mp_real_data, chiModule);
+ m_fitSuite->getFitObjects()->setExperimentNormalize(true);
m_fitSuite->runFit();
diff --git a/App/src/TestIsGISAXS12.cpp b/App/src/TestIsGISAXS12.cpp
index 4d580c079e0f1adb0def705ab7ed52ac0ac50b7d..e74ff6746242b11c1da170f86b3f130cf086038b 100644
--- a/App/src/TestIsGISAXS12.cpp
+++ b/App/src/TestIsGISAXS12.cpp
@@ -24,6 +24,7 @@
#include "ROOTMinimizer.h"
#include "OutputDataFunctions.h"
#include "NamedVector.h"
+#include "Types.h"
#include <iostream>
#include <fstream>
@@ -34,6 +35,8 @@
#include "TH1D.h"
#include "TH2D.h"
#include "TLine.h"
+#include "TROOT.h"
+#include "TLegend.h"
@@ -57,9 +60,6 @@ TestIsGISAXS12::~TestIsGISAXS12()
delete m_experiment;
delete m_sample_builder;
delete m_fitSuite;
- for(DataScan_t::iterator it=m_isgi_scans.begin(); it!=m_isgi_scans.end(); ++it) {
- delete (*it);
- }
}
@@ -73,13 +73,18 @@ void TestIsGISAXS12::execute()
initialiseExperiment();
// run our standard isgisaxs comparison for given sample
- run_isgisaxs_comparison();
+// run_isgisaxs_comparison();
+
+ // plot isgisaxs data
+// plot_isgisaxs_data();
// run test fit
- //run_test_fit();
+// run_test_minimizer();
// run isgisaxs ex-12 style fit
- run_isgisaxs_fit();
+ run_isgisaxs_fit();
+
+// run_test_chimodule();
}
@@ -95,7 +100,7 @@ void TestIsGISAXS12::run_isgisaxs_comparison()
IsGISAXSTools::writeOutputDataToFile(*(m_experiment->getOutputData()), m_data_path+"this_fitconstraints.ima");
// plotting results of comparison we/isgisaxs for the sample with default parameters
- std::string isgi_file(m_data_path+"isgi_fitconstraints.ima");
+ std::string isgi_file(m_data_path+"isgi_fitconstraints_optimal.ima");
std::string this_file(m_data_path+"this_fitconstraints.ima");
// -------------
@@ -105,12 +110,12 @@ void TestIsGISAXS12::run_isgisaxs_comparison()
OutputData<double> *our_data = IsGISAXSTools::readOutputDataFromFile(this_file);
TCanvas *c1 = DrawHelper::instance().createAndRegisterCanvas("TestIsGISAXS12_c1", "ex-12: Mixture of cylindrical particles with different size distribution");
-
c1->Divide(2,2);
- IsGISAXSTools::setMinimum(1.);
+// IsGISAXSTools::setMinimum(1.);
// our calculations
c1->cd(1); gPad->SetLogz();
+// IsGISAXSTools::drawOutputDataInPad(*(m_experiment->getOutputData()), "CONT4 Z", "this");
IsGISAXSTools::drawOutputDataInPad(*our_data, "CONT4 Z", "this");
// isgisaxs data
@@ -132,82 +137,7 @@ void TestIsGISAXS12::run_isgisaxs_comparison()
delete isgi_data;
delete our_data;
-}
-
-
-/* ************************************************************************* */
-// run test fitting
-/* ************************************************************************* */
-void TestIsGISAXS12::run_test_fit()
-{
- // creating "real" data
- //m_experiment->setDetectorResolutionFunction(new ResolutionFunction2DSimple(0.0002, 0.0002));
- m_experiment->setBeamIntensity(1e10);
- m_experiment->runSimulation();
- m_experiment->normalize();
-
- OutputData<double > *real_data = IsGISAXSTools::createNoisyData( *m_experiment->getOutputData() );
-
- // setting up 1d scans by making slices on real data
- DataScan_t data_scans;
- data_scans.push_back( OutputDataFunctions::selectRangeOnOneAxis(*real_data, "alpha_f", 0.012, 0.012) );
- data_scans.push_back( OutputDataFunctions::selectRangeOnOneAxis(*real_data, "phi_f", 0.011, 0.011) );
-
- // drawing data and scans
- TCanvas *c1 = new TCanvas("c1","c1",1024, 768);
- c1->Divide(2,2);
- c1->cd(1); gPad->SetLogz();
- TH2D *hist1 = dynamic_cast<TH2D *>(IsGISAXSTools::getOutputDataTH123D( *real_data, "real_data"));
- hist1->Draw("COLZ");
- for(DataScan_t::const_iterator it=data_scans.begin(); it!= data_scans.end(); ++it) {
- TLine *line = IsGISAXSTools::getOutputDataScanLine(*(*it));
- line->DrawClone();
- delete line;
- }
-
- int npad(2);
- for(DataScan_t::iterator it=data_scans.begin(); it!= data_scans.end(); ++it, ++npad) {
- c1->cd(npad);
- TH1D *hist = IsGISAXSTools::getOutputDataScanHist(*(*it));
- hist->DrawCopy();
- delete hist;
- }
c1->Update();
-
- // creating fit suite
- m_fitSuite = new FitSuite();
- m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
- m_fitSuite->attachObserver( new FitSuiteObserverPrint() );
- m_fitSuite->attachObserver( new FitSuiteObserverDraw() );
-
-// m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 0.2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 0.2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 0.8*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 0.8*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 12*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 6*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/particle_probability", 0.4, 0.1, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 4*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-// m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 4*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-
- m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 0.3*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 0.3*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 0.7*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 0.7*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 10*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 5*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/particle_probability", 0.4, 0.1, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 3*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
- m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 3*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
-
- for(DataScan_t::iterator it=data_scans.begin(); it!= data_scans.end(); ++it) {
- m_fitSuite->addExperimentAndRealData(*m_experiment, *(*it));
- }
-
- m_fitSuite->runFit();
-
- delete real_data;
- for(DataScan_t::iterator it=data_scans.begin(); it!= data_scans.end(); ++it) delete (*it);
}
@@ -217,16 +147,13 @@ void TestIsGISAXS12::run_test_fit()
void TestIsGISAXS12::run_isgisaxs_fit()
{
// reading info about two 1D scans defined in isgisaxs example
- read_isgisaxs_datfile(m_data_path+"isgi_fitconstraints.dat");
-
- // making some plots together with test simulation
- plot_isgisaxs_data();
-// return;
+ DataSet isgi_scans;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_scans);
// creating fit suite
m_fitSuite = new FitSuite();
m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
- //m_fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiMin", "SteepestDescent") );
+ //m_fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiFit", "") );
m_fitSuite->attachObserver( new FitSuiteObserverPrint(10) );
m_fitSuite->attachObserver( new FitSuiteObserverDraw(50) );
@@ -234,115 +161,305 @@ void TestIsGISAXS12::run_isgisaxs_fit()
minim->SetStrategy(1);
// minim->SetPrecision(1.);
- m_fitSuite->addFitParameter("*Normalizer/scale", 1e5, 100, AttLimits::limited(1e4, 2e5));
- m_fitSuite->addFitParameter("*Normalizer/shift", 10, 1, AttLimits::limited(1., 20.));
- m_fitSuite->addFitParameter("*SampleBuilder/particle_probability", 0.4, 0.1, AttLimits::limited(0.01, 1.0) );
+ m_fitSuite->addFitParameter("*Normalizer/scale", 1e5, 1, AttLimits::limited(1e4, 2e5));
+ m_fitSuite->addFitParameter("*Normalizer/shift", 10, 0.01, AttLimits::limited(1., 20.));
+
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_probability1", 0.4, 0.01, AttLimits::limited(0.01, 1.0) );
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 4*Units::nanometer, 0.01*Units::nanometer, AttLimits::limited(1., 10.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 0.2, 0.01, AttLimits::limited(0.01, 1.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 0.8, 0.01, AttLimits::limited(0.01, 10.) );
+
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_probability2", 0.6, 0.01, AttLimits::limited(0.01, 1.0) );
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 4*Units::nanometer, 0.01*Units::nanometer, AttLimits::limited(1., 10.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 0.2, 0.01, AttLimits::limited(0.01, 1.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 0.8, 0.01, AttLimits::limited(0.01, 10.) );
+
+ m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 12*Units::nanometer, 0.01*Units::nanometer, AttLimits::limited(0.01, 50.0) );
+ m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 6*Units::nanometer, 0.01*Units::nanometer, AttLimits::limited(0.01, 10.) );
- m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 4*Units::nanometer, 1*Units::nanometer, AttLimits::limited(1., 10.) );
- m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 0.2, 0.1, AttLimits::limited(0.01, 1.) );
- m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 0.8, 0.1, AttLimits::limited(0.01, 10.) );
- m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 4*Units::nanometer, 1*Units::nanometer, AttLimits::limited(1., 10.) );
- m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 0.2, 0.1, AttLimits::limited(0.01, 1.) );
- m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 0.8, 0.1, AttLimits::limited(0.01, 10.) );
+// m_fitSuite->addFitParameter("*Normalizer/scale", 1.31159E+05, 100, AttLimits::limited(1e4, 2e5));
+// m_fitSuite->addFitParameter("*Normalizer/shift", -8.10009E-02, 1, AttLimits::limited(-10., 20.));
- m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 12*Units::nanometer, 1*Units::nanometer, AttLimits::limited(0.01, 50.0) );
- m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 6*Units::nanometer, 1*Units::nanometer, AttLimits::limited(0.01, 10.) );
+// m_fitSuite->addFitParameter("*SampleBuilder/particle_probability1", 5.34055E-01, 0.1, AttLimits::limited(0.01, 1.0) );
+// m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 4.90801E+00, 1*Units::nanometer, AttLimits::limited(1., 10.) );
+// m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 1.90651E-01, 0.1, AttLimits::limited(0.01, 1.) );
+// m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 1.00193E+00, 0.1, AttLimits::limited(0.01, 10.) );
+
+// m_fitSuite->addFitParameter("*SampleBuilder/particle_probability2", 4.70783E-01, 0.1, AttLimits::limited(0.01, 1.0) );
+// m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 5.16801E+00, 1*Units::nanometer, AttLimits::limited(1., 10.) );
+// m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 2.03908E-01, 0.1, AttLimits::limited(0.01, 1.) );
+// m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 9.77402E-01, 0.1, AttLimits::limited(0.01, 10.) );
+
+// m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 1.49681E+01, 1*Units::nanometer, AttLimits::limited(0.01, 50.0) );
+// m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 3.03315E+00, 1*Units::nanometer, AttLimits::limited(0.01, 10.) );
// setting up fitSuite
ChiSquaredModule chiModule;
chiModule.setChiSquaredFunction( SquaredFunctionWithSystematicError(0.08) );
chiModule.setOutputDataNormalizer( OutputDataNormalizerScaleAndShift() );
- for(DataScan_t::iterator it=m_isgi_scans.begin(); it!= m_isgi_scans.end(); ++it) {
+ for(DataSet::iterator it=isgi_scans.begin(); it!= isgi_scans.end(); ++it) {
m_fitSuite->addExperimentAndRealData(*m_experiment, *(*it), chiModule);
}
m_fitSuite->runFit();
// drawing results
- TCanvas *c2 = new TCanvas("c2","c2",800,600);
+ TCanvas *c2 = new TCanvas("c2","GISASFW fit results",800,600);
c2->Divide(2,2);
- for(size_t i=0; i<m_fitSuite->getFitObjects()->size(); ++i) {
- c2->cd((int)i+1);
- const FitObject *obj = m_fitSuite->getFitObjects()->getObject(i);
- TH1D *hreal = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getRealData(),"real");
- TH1D *hsimul = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getSimulationData(),"simul");
+ TLegend *leg1 = new TLegend(0.5,0.6,0.85,0.85);
+ leg1->SetBorderSize(1);
+ leg1->SetFillStyle(0);
+ for(size_t i_set=0; i_set<m_fitSuite->getFitObjects()->size(); ++i_set) {
+ c2->cd(i_set+1);
+ const FitObject *obj = m_fitSuite->getFitObjects()->getObject(i_set);
+ TH1D *hreal = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getRealData(),"gisasfw_real");
+ TH1D *hsimul = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getSimulationData(),"gisasfw_simul");
hreal->SetLineColor(kBlue);
gPad->SetLogy();
hreal->DrawCopy();
hsimul->DrawCopy("same");
- delete hreal;
- delete hsimul;
+ if(i_set==0) leg1->AddEntry(hreal,"GISASFW data","lp");
+ if(i_set==0) leg1->AddEntry(hsimul,"GISASFW simul","lp");
+ }
+ c2->cd(1); leg1->Draw();
+ c2->cd(2); leg1->Draw();
+
+ // drawing ratio
+ TLegend *leg2 = new TLegend(0.5,0.6,0.85,0.85);
+ leg2->SetBorderSize(1);
+ leg2->SetFillStyle(0);
+ for(size_t i_set=0; i_set<m_fitSuite->getFitObjects()->size(); ++i_set) {
+ c2->cd(3+1);
+ const FitObject *obj = m_fitSuite->getFitObjects()->getObject(i_set);
+ OutputData<double > *real = obj->getChiSquaredModule()->getRealData()->clone();
+ OutputData<double > *simul = obj->getChiSquaredModule()->getSimulationData()->clone();
+
+ c2->cd(i_set+3);
+ *simul /= *real;
+ TH1D *hratio = IsGISAXSTools::getOutputDataScanHist(*simul,"gisasfw_real_simul_ratio");
+ hratio->DrawCopy();
+ if(i_set==0) {
+ leg2->AddEntry(hratio,"GISASFW simul/real","lp");
+ }
+ delete real;
+ delete simul;
}
+ c2->cd(3); leg2->Draw();
+ c2->cd(4); leg2->Draw();
+ c2->Update();
}
+/* ************************************************************************* */
+// run chi module test on isgisaxs data/result pair to check module numericaly
+/* ************************************************************************* */
+void TestIsGISAXS12::run_test_chimodule()
+{
+ DataSet isgi_scans;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_scans);
+
+ DataSet isgi_results;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_results, true);
+
+ // setting up fitSuite
+ ChiSquaredModule chiModule;
+ chiModule.setChiSquaredFunction( SquaredFunctionWithSystematicError(0.08) );
+ chiModule.setOutputDataNormalizer( OutputDataNormalizerScaleAndShift() );
+
+ OutputDataNormalizerScaleAndShift normalizer(1.31159E+05, -8.10009E-02);
+
+ double max_intensity(0);
+ for(size_t i=0; i<isgi_results.size(); ++i) {
+ OutputData<double>::const_iterator cit = std::max_element(isgi_results[i]->begin(), isgi_results[i]->end());
+ max_intensity = std::max(max_intensity, *cit);
+ }
+ std::cout << "XXX " << max_intensity << std::endl;
+ normalizer.setMaximumIntensity(max_intensity);
+
+ chiModule.setOutputDataNormalizer( normalizer );
+
+ double chi_sum(0);
+ for(size_t i=0; i<isgi_scans.size(); ++i) {
+ chiModule.setRealAndSimulatedData(*isgi_scans[i], *isgi_results[i]);
+ std::cout << " AAA " << isgi_scans.size()*isgi_results[i]->getAllocatedSize() - 12 << std::endl;
+ chiModule.setNdegreeOfFreedom(isgi_scans.size()*isgi_results[i]->getAllocatedSize() - 12);
+ double chi = 0.5*0.5*chiModule.calculateChiSquared();
+ chi_sum += chi;
+ std::cout << "chi : " << chi << " chi_sum:" << chi_sum << std::endl;
+ }
+ std::cout << "chi_sum " << chi_sum << std::endl;
+
+ return;
+}
+
/* ************************************************************************* */
// plot isgisaxs data together with test simulation
/* ************************************************************************* */
void TestIsGISAXS12::plot_isgisaxs_data()
{
- if(m_isgi_scans.size() != 2) throw LogicErrorException("TestIsGISAXS12::plot_isgisaxs_data() -> Load isgisaxs scans first");
-
- enum isgisaxs_scans {kFixedAlphaf, kFixedPhif };
-
- // making 2d OutputData with axes like in two isgisaxs scans
- OutputData<double > *data_as_isgisaxs = new OutputData<double>;
- const NamedVector<double > *axis_phif_scan = dynamic_cast<const NamedVector<double> *>(m_isgi_scans[kFixedAlphaf]->getAxis("phi_f"));
- const NamedVector<double > *axis_alphaf_scan = dynamic_cast<const NamedVector<double> *>(m_isgi_scans[kFixedPhif]->getAxis("alpha_f"));
- if( !axis_phif_scan || !axis_phif_scan ) throw NullPointerException("TestIsGISAXS12::plot_isgisaxs_data() -> Error. Can'get get axis");
- data_as_isgisaxs->addAxis(axis_phif_scan->clone());
- data_as_isgisaxs->addAxis(axis_alphaf_scan->clone());
- data_as_isgisaxs->setAllTo(0.0);
- std::cout << axis_phif_scan->getSize() << " " << axis_alphaf_scan->getSize() << " " << data_as_isgisaxs->getAllocatedSize() << std::endl;
-
- // running 2d simulation
- m_experiment->setDetectorParameters(*data_as_isgisaxs);
- delete data_as_isgisaxs;
-// m_experiment->setBeamIntensity(1e10);
- m_experiment->runSimulation();
-// m_experiment->normalize();
+ // reading two 1D scans defined in isgisaxs example (399 points/scan)
+ DataSet isgi_scans;
+ read_isgisaxs_datfile(m_data_path+"isgi_fitconstraints.dat", isgi_scans);
+
+ // reading isgisaxs scans which actually have been used for a fit together with fit results (100 points/scan)
+ DataSet isgi_scans_smoothed;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_scans_smoothed);
+ DataSet isgi_results;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_results, true);
- // setting up 1d scans by making slices on just simulated data
- DataScan_t data_scans;
- data_scans.push_back( OutputDataFunctions::selectRangeOnOneAxis(*m_experiment->getOutputData(), "alpha_f", m_isgi_fixed_alphaf, m_isgi_fixed_alphaf ));
- data_scans.push_back( OutputDataFunctions::selectRangeOnOneAxis(*m_experiment->getOutputData(), "phi_f", m_isgi_fixed_phif, m_isgi_fixed_phif ));
+ print_axes(isgi_scans);
+ print_axes(isgi_scans_smoothed);
+ print_axes(isgi_results);
- // drawing data and scans
- TCanvas *c1 = new TCanvas("c1","c1",1024, 768);
+ TCanvas *c1 = DrawHelper::instance().createAndRegisterCanvas("c1_isgisaxs_data", "Looking on IsGISAXS data and fit results", 768, 1024);
c1->Divide(2,3);
- c1->cd(1); gPad->SetLogz();
- TH2D *hist1 = dynamic_cast<TH2D *>(IsGISAXSTools::getOutputDataTH123D( *m_experiment->getOutputData(), "real_data"));
- hist1->Draw("COLZ");
- for(DataScan_t::const_iterator it=data_scans.begin(); it!= data_scans.end(); ++it) {
- TLine *line = IsGISAXSTools::getOutputDataScanLine(*(*it));
- line->DrawClone();
- delete line;
+
+ // drawing real data with fine and coars granularity on top of each other
+ TLegend *leg1 = new TLegend(0.5,0.6,0.85,0.85);
+ leg1->SetBorderSize(1);
+ leg1->SetFillStyle(0);
+
+ for(size_t i_set=0; i_set<isgi_scans.size(); ++i_set) {
+ c1->cd(1+i_set); gPad->SetLogy();
+ TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*isgi_scans[i_set], "data");
+ hdata->DrawCopy();
+ if(i_set==0) leg1->AddEntry(hdata,"isgisaxs data fine","lp");
+ TH1D *hdata_smoothed = IsGISAXSTools::getOutputDataScanHist(*isgi_scans_smoothed[i_set], "data_smoothed");
+ hdata_smoothed->SetLineColor(kBlue);
+ hdata_smoothed->DrawCopy("same");
+ if(i_set==0) leg1->AddEntry(hdata_smoothed,"isgisaxs data","lp");
}
- int npad(3);
- for(DataScan_t::iterator it=data_scans.begin(); it!= data_scans.end(); ++it, ++npad) {
- c1->cd(npad);
- TH1D *hist = IsGISAXSTools::getOutputDataScanHist(*(*it), "sim");
- hist->DrawCopy();
- delete hist;
+ c1->cd(1); leg1->Draw();
+ c1->cd(2); leg1->Draw();
+
+ // drawing isgsaxs fit results on top of isgisaxs real data
+ TLegend *leg2 = new TLegend(0.5,0.6,0.85,0.85);
+ leg2->SetBorderSize(1);
+ leg2->SetFillStyle(0);
+ for(size_t i_set=0; i_set<isgi_scans.size(); ++i_set) {
+ c1->cd(3+i_set); gPad->SetLogy();
+ TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*isgi_results[i_set], "data");
+ hdata->SetLineColor(kRed);
+ hdata->DrawCopy();
+ if(i_set==0) leg2->AddEntry(hdata,"isgisaxs results","lp");
+ TH1D *hdata_smoothed = IsGISAXSTools::getOutputDataScanHist(*isgi_scans_smoothed[i_set], "data_smoothed");
+ hdata_smoothed->SetLineColor(kBlue);
+ hdata_smoothed->DrawCopy("same");
+ if(i_set==0) leg2->AddEntry(hdata_smoothed,"isgisaxs data","lp");
}
-
- // drawing isgi scans
- npad=5;
- for(DataScan_t::iterator it=m_isgi_scans.begin(); it!= m_isgi_scans.end(); ++it, ++npad) {
- c1->cd(npad);
- TH1D *hist = IsGISAXSTools::getOutputDataScanHist(*(*it),"isgi");
- hist->SetLineColor(kBlue);
- hist->DrawCopy();
- delete hist;
+ c1->cd(3); leg2->Draw();
+ c1->cd(4); leg2->Draw();
+
+ // drawing ratio of isgisaxs data/results
+ TLegend *leg3 = new TLegend(0.5,0.6,0.85,0.85);
+ leg3->SetBorderSize(1);
+ leg3->SetFillStyle(0);
+ for(size_t i_set=0; i_set<isgi_scans.size(); ++i_set) {
+ c1->cd(5+i_set);
+ *isgi_results[i_set] /= *isgi_scans_smoothed[i_set];
+ TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*isgi_results[i_set], "isgisaxs_results_data");
+ hdata->SetLineColor(kRed);
+ hdata->DrawCopy();
+ if(i_set==0) leg3->AddEntry(hdata,"isgisaxs results/data","lp");
}
+ c1->cd(5); leg3->Draw();
+ c1->cd(6); leg3->Draw();
+
c1->Update();
- for(DataScan_t::iterator it=data_scans.begin(); it!= data_scans.end(); ++it) delete (*it);
+}
+
+
+/* ************************************************************************* */
+// run test minimizer to check the whole chain
+/* ************************************************************************* */
+void TestIsGISAXS12::run_test_minimizer()
+{
+ // reading isgisaxs real data
+ DataSet isgi_scans_smoothed;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_scans_smoothed);
+ // isgisaxs fit results
+ DataSet isgi_results;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_results, true);
+
+ // Putting parameters found by isgisaxs into our sample and run FitSuite once with the help of TestMinimizer to see if
+ // our simulation produces numerically same results
+
+ m_fitSuite = new FitSuite();
+ m_fitSuite->setMinimizer( new TestMinimizer() );
+
+ m_fitSuite->addFitParameter("*Normalizer/scale", 1.31159E+05, 100, AttLimits::limited(1e4, 2e5));
+ m_fitSuite->addFitParameter("*Normalizer/shift", -8.10009E-02, 1, AttLimits::limited(-10., 20.));
+
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_probability1", 5.34055E-01, 0.1, AttLimits::limited(0.01, 1.0) );
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 4.90801E+00, 1*Units::nanometer, AttLimits::limited(1., 10.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 1.90651E-01, 0.1, AttLimits::limited(0.01, 1.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 1.00193E+00, 0.1, AttLimits::limited(0.01, 10.) );
+
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_probability2", 4.70783E-01, 0.1, AttLimits::limited(0.01, 1.0) );
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 5.16801E+00, 1*Units::nanometer, AttLimits::limited(1., 10.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 2.03908E-01, 0.1, AttLimits::limited(0.01, 1.) );
+ m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 9.77402E-01, 0.1, AttLimits::limited(0.01, 10.) );
+
+ m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 1.49681E+01, 1*Units::nanometer, AttLimits::limited(0.01, 50.0) );
+ m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 3.03315E+00, 1*Units::nanometer, AttLimits::limited(0.01, 10.) );
+
+ // setting up fitSuite
+ ChiSquaredModule chiModule;
+ chiModule.setChiSquaredFunction( SquaredFunctionWithSystematicError(0.08) );
+ chiModule.setOutputDataNormalizer( OutputDataNormalizerScaleAndShift() );
+// for(DataSet::iterator it=isgi_results.begin(); it!= isgi_results.end(); ++it) {
+// m_fitSuite->addExperimentAndRealData(*m_experiment, *(*it), chiModule);
+// }
+ for(DataSet::iterator it=isgi_scans_smoothed.begin(); it!= isgi_scans_smoothed.end(); ++it) {
+ m_fitSuite->addExperimentAndRealData(*m_experiment, *(*it), chiModule);
+ }
+ m_fitSuite->runFit();
+
+ TCanvas *c1 = new TCanvas("c1_test_minimizer","TestMinimizer", 800, 600);
+ c1->Divide(2,2);
+
+ // drawing GISASFW simul on top of isgisaxs simul
+ TLegend *leg1 = new TLegend(0.5,0.6,0.85,0.85);
+ leg1->SetBorderSize(1);
+ leg1->SetFillStyle(0);
+ for(size_t i_set=0; i_set<isgi_results.size(); ++i_set) {
+ c1->cd(1+i_set);
+ gPad->SetLogy();
+ TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*isgi_results[i_set], "data");
+ hdata->SetLineColor(kRed);
+ hdata->DrawCopy();
+ const OutputData<double > *data = m_fitSuite->getFitObjects()->getObject(i_set)->getChiSquaredModule()->getSimulationData();
+ TH1D *simul_data = IsGISAXSTools::getOutputDataScanHist(*data, "data_from_module");
+ simul_data->SetLineColor(kBlue);
+ simul_data->DrawCopy("same");
+
+ if(i_set==0) leg1->AddEntry(hdata,"isgisaxs results","lp");
+ if(i_set==0) leg1->AddEntry(simul_data,"gisasfw simul","lp");
+ }
+ c1->cd(1); leg1->Draw();
+ c1->cd(2); leg1->Draw();
+
+ TLegend *leg2 = new TLegend(0.5,0.6,0.85,0.85);
+ leg2->SetBorderSize(1);
+ leg2->SetFillStyle(0);
+ for(size_t i_set=0; i_set<isgi_results.size(); ++i_set) {
+ c1->cd(3+i_set);
+ OutputData<double > *data = m_fitSuite->getFitObjects()->getObject(i_set)->getChiSquaredModule()->getSimulationData()->clone();
+ *data /= *isgi_results[i_set];
+ TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*data, "gisasfw_isgisaxs_simul");
+ hdata->SetLineColor(kRed);
+ hdata->DrawCopy();
+ delete data;
+ if(i_set==0) leg2->AddEntry(hdata,"gisasfw/isgisaxs simul","lp");
+ }
+ c1->cd(3); leg1->Draw();
+ c1->cd(4); leg1->Draw();
}
@@ -366,8 +483,10 @@ void TestIsGISAXS12::initialiseExperiment()
/* ************************************************************************* */
// read special isgisaxs *.dat file with data to fit
/* ************************************************************************* */
-void TestIsGISAXS12::read_isgisaxs_datfile(const std::string &filename)
+void TestIsGISAXS12::read_isgisaxs_datfile(const std::string &filename, DataSet &dataset)
{
+ dataset.clear();
+
// opening ASCII file
std::ifstream fin;
fin.open(filename.c_str(), std::ios::in);
@@ -387,8 +506,7 @@ void TestIsGISAXS12::read_isgisaxs_datfile(const std::string &filename)
n_dataset_line = 0; // it's a beginning of new data set ("cross-section" in isgisaxs terminology)
if( !isgiScan.empty() ) {
OutputData<double > *data = convert_isgi_scan(isgiScan);
- m_isgi_scans.push_back(data);
- std::cout << "YYY " << isgiScan.size() << std::endl;
+ dataset.push_back(data);
isgiScan.clear();
}
}
@@ -400,11 +518,9 @@ void TestIsGISAXS12::read_isgisaxs_datfile(const std::string &filename)
if ( !(iss >> ctmp >> isgiData.phif >> isgiData.alphaf >> isgiData.intensity) ) throw FormatErrorException("TestIsGISAXS12::read_isgisaxs_datfile() -> Error!");
iss >> isgiData.err; // column with errors can be absent in file, so no check for success here
ctmp == 'T' ? isgiData.use_it = true : isgiData.use_it = false;
- if(isgiData.phif > .0001) {
- isgiData.phif = std::asin(isgiData.phif); // because isgisax in fact stores in *.dat file sin(phif), and sin(alphaf) instead of phif, alphaf
- isgiData.alphaf = std::asin(isgiData.alphaf); // because isgisax in fact stores in *.dat file sin(phif), and sin(alphaf) instead of phif, alphaf
- isgiScan.push_back(isgiData);
- }
+ isgiData.phif = std::asin(isgiData.phif); // because isgisax in fact stores in *.dat file sin(phif), and sin(alphaf) instead of phif, alphaf
+ isgiData.alphaf = std::asin(isgiData.alphaf); // because isgisax in fact stores in *.dat file sin(phif), and sin(alphaf) instead of phif, alphaf
+ isgiScan.push_back(isgiData);
}
n_dataset_line++;
@@ -417,6 +533,71 @@ void TestIsGISAXS12::read_isgisaxs_datfile(const std::string &filename)
}
+/* ************************************************************************* */
+// read special isgisaxs *.out file with isgisaxs adjusted data and fit results
+//
+// if read_fit_results == false, then it loads isgisaxs data to fit
+// if read_fit_results == true, then it loads isgisaxs fit results
+/* ************************************************************************* */
+void TestIsGISAXS12::read_isgisaxs_outfile(const std::string &filename, DataSet &dataset, bool read_fit_results)
+{
+ enum isgisaxs_keys_outfile {kSin_twotheta, kSin_alphaf, kQx, kQy, kQz, kGISAXS, kData2fit, kErrorbar, kIobs_Icalc, kFitted };
+
+ dataset.clear();
+
+ // opening ASCII file
+ std::ifstream fin;
+ fin.open(filename.c_str(), std::ios::in);
+ if( !fin.is_open() ) {
+ throw FileNotIsOpenException("TestIsGISAXS12::read_isgisaxs_datfile() -> Error. Can't open file '"+filename+"' for reading.");
+ }
+
+ std::vector<IsgiData > isgiScan;
+
+ std::string sline;
+ while( std::getline(fin, sline))
+ {
+ std::string::size_type pos=sline.find("# Cut # =");
+ if( pos!= std::string::npos ) {
+ //std::cout << "beginning of one dataset " << isgiScan.size() << std::endl;
+ if( !isgiScan.empty() ) {
+ OutputData<double > *data = convert_isgi_scan(isgiScan);
+ dataset.push_back(data);
+ isgiScan.clear();
+ }
+ } else if (sline[0] != '#' && sline.size() >2){
+ // not a header, parsing records
+ std::istringstream iss(sline.c_str());
+ IsgiData isgiData;
+
+ vdouble1d_t buff;
+ std::copy(std::istream_iterator<double>(iss), std::istream_iterator<double>(), back_inserter(buff));
+ if( buff.size() != kFitted+1) {
+ throw LogicErrorException("TestIsGISAXS12::read_isgisaxs_outfile -> Error! Line doesn't have enough double numbers. Not an *.out file? Line '"+sline+"'");
+ }
+
+ isgiData.phif = std::asin(buff[kSin_twotheta]);
+ isgiData.alphaf = std::asin(buff[kSin_alphaf]);
+ isgiData.intensity = buff[kData2fit];
+ if( read_fit_results ) isgiData.intensity = buff[kGISAXS];
+
+ isgiScan.push_back(isgiData);
+ }
+ }
+ if ( fin.bad() ) {
+ throw FileIsBadException("TestIsGISAXS12::read_isgisaxs_datfile() -> Error! File is bad after readline(), probably it is a directory.");
+ }
+ fin.close();
+
+ if( !isgiScan.empty() ) {
+ OutputData<double > *data = convert_isgi_scan(isgiScan);
+ dataset.push_back(data);
+ isgiScan.clear();
+ }
+}
+
+
+
/* ************************************************************************* */
// convert isgisaxs 1d scan to output data 2d object
/* ************************************************************************* */
@@ -424,6 +605,14 @@ OutputData<double> *TestIsGISAXS12::convert_isgi_scan(std::vector<IsgiData > &is
{
if(isgi_data.size() <2 ) throw LogicErrorException("TestIsGISAXS12::convert_isgi_scan() -> Error! Too short vector.");
+// // adjust data to get rid from the point with phi_f~0.0
+// std::vector<IsgiData > adjusted_data;
+// for(size_t i_point=0; i_point<isgi_data.size(); ++i_point) {
+// if(isgi_data[i_point].phif > 0.001) {
+// adjusted_data.push_back(isgi_data[i_point]);
+// }
+// }
+
// check if it is scan with fixed phi_f or with fixed alpha_f
bool fixed_phif(true);
bool fixed_alphaf(true);
@@ -482,15 +671,18 @@ OutputData<double> *TestIsGISAXS12::convert_isgi_scan(std::vector<IsgiData > &is
//
/* ************************************************************************* */
TestIsGISAXS12::TestSampleBuilder::TestSampleBuilder()
- : m_particle_probability(0.4)
- , m_particle_radius1(4*Units::nanometer)
- , m_particle_radius2(4*Units::nanometer)
- , m_dispersion_radius1(0.2)
- , m_dispersion_radius2(0.2)
- , m_height_aspect_ratio1(0.8)
- , m_height_aspect_ratio2(0.8)
- , m_interf_distance(12.0*Units::nanometer)
- , m_interf_width(6*Units::nanometer)
+// optimal parameters as found by isgisaxs
+ : m_particle_probability1(5.34055E-01)
+ , m_particle_radius1(4.90801E+00*Units::nanometer)
+ , m_dispersion_radius1(1.90651E-01)
+ , m_height_aspect_ratio1(1.00193E+00)
+ , m_particle_probability2(4.70783E-01)
+ , m_particle_radius2(5.16801E+00*Units::nanometer)
+ , m_dispersion_radius2(2.03908E-01)
+ , m_height_aspect_ratio2(9.77402E-01)
+ , m_interf_distance(1.49681E+01*Units::nanometer)
+ , m_interf_width(3.03315E+00*Units::nanometer)
+
{
init_parameters();
}
@@ -498,12 +690,13 @@ TestIsGISAXS12::TestSampleBuilder::TestSampleBuilder()
void TestIsGISAXS12::TestSampleBuilder::init_parameters()
{
getParameterPool()->clear();
- getParameterPool()->registerParameter("particle_probability", &m_particle_probability);
+ getParameterPool()->registerParameter("particle_probability1", &m_particle_probability1);
getParameterPool()->registerParameter("particle_radius1", &m_particle_radius1);
- getParameterPool()->registerParameter("particle_radius2", &m_particle_radius2);
getParameterPool()->registerParameter("dispersion_radius1", &m_dispersion_radius1);
- getParameterPool()->registerParameter("dispersion_radius2", &m_dispersion_radius2);
getParameterPool()->registerParameter("height_aspect_ratio1", &m_height_aspect_ratio1);
+ getParameterPool()->registerParameter("particle_probability2", &m_particle_probability2);
+ getParameterPool()->registerParameter("particle_radius2", &m_particle_radius2);
+ getParameterPool()->registerParameter("dispersion_radius2", &m_dispersion_radius2);
getParameterPool()->registerParameter("height_aspect_ratio2", &m_height_aspect_ratio2);
getParameterPool()->registerParameter("interf_distance", &m_interf_distance);
getParameterPool()->registerParameter("interf_width", &m_interf_width);
@@ -524,6 +717,10 @@ ISample *TestIsGISAXS12::TestSampleBuilder::buildSample() const
Layer air_layer(air_material);
// preparing nano particles prototypes for seeding layer's particle_decoration
+ double particle_probability1 = m_particle_probability1;
+// double particle_probability2 = 1. - m_particle_probability1;
+ double particle_probability2 = m_particle_probability2;
+
double radius1 = m_particle_radius1;
double radius2 = m_particle_radius2;
double height1 = m_height_aspect_ratio1*radius1;
@@ -548,9 +745,9 @@ ISample *TestIsGISAXS12::TestSampleBuilder::buildSample() const
// building nano particles
ParticleBuilder builder;
- builder.setPrototype(cylinder1,"/Particle/FormFactorCylinder/radius", par1, m_particle_probability);
+ builder.setPrototype(cylinder1,"/Particle/FormFactorCylinder/radius", par1, particle_probability1);
builder.plantParticles(particle_decoration);
- builder.setPrototype(cylinder2,"/Particle/FormFactorCylinder/radius", par2, 1-m_particle_probability);
+ builder.setPrototype(cylinder2,"/Particle/FormFactorCylinder/radius", par2, particle_probability2);
builder.plantParticles(particle_decoration);
// making layer holding all whose nano particles
@@ -565,6 +762,18 @@ ISample *TestIsGISAXS12::TestSampleBuilder::buildSample() const
}
+void TestIsGISAXS12::print_axes(DataSet &data)
+{
+ for(size_t i_set=0; i_set<data.size(); ++i_set) {
+ std::cout << "scan #" << i_set << " ";
+ for(size_t i_axis=0; i_axis<data[i_set]->getNdimensions(); ++i_axis) {
+ const NamedVector<double > *axis = dynamic_cast<const NamedVector<double > *>(data[i_set]->getAxis(i_axis));
+ if( !axis) throw NullPointerException("TestIsGISAXS12::plot_isgisaxs_data() -> Error! Can't get axis");
+ std::cout << "( " << axis->getName() << ", " << axis->getSize() << ", " << axis->getMin() << ", " << axis->getMax() << " ) ";
+ }
+ std::cout << std::endl;
+ }
+}
diff --git a/Core/Algorithms/inc/IChiSquaredModule.h b/Core/Algorithms/inc/IChiSquaredModule.h
index 169647ce196b62ce1a25262dfe7af39b42170b9e..37d8a3b6ab78bae7ea078470d5da90bb939aec7b 100644
--- a/Core/Algorithms/inc/IChiSquaredModule.h
+++ b/Core/Algorithms/inc/IChiSquaredModule.h
@@ -58,6 +58,7 @@ public:
//! get data normalizer
virtual const IOutputDataNormalizer *getOutputDataNormalizer() const {return mp_data_normalizer; }
+ virtual IOutputDataNormalizer *getOutputDataNormalizer() {return mp_data_normalizer; }
//! set data normalizer
virtual void setOutputDataNormalizer(const IOutputDataNormalizer &data_normalizer);
@@ -69,6 +70,9 @@ public:
//! return last calculated chi squared value
virtual double getValue() const { return m_chi2_value; }
+ //! set number of degree of freedom
+ void setNdegreeOfFreedom(int ndegree_of_freedom) { m_ndegree_of_freedom = ndegree_of_freedom; }
+
protected:
// disabling assignment operator
IChiSquaredModule &operator=(const IChiSquaredModule &);
@@ -81,6 +85,7 @@ protected:
IFittingDataSelector *mp_data_selector;
IOutputDataNormalizer *mp_data_normalizer;
IIntensityFunction *mp_intensity_function;
+ int m_ndegree_of_freedom;
double m_chi2_value;
};
diff --git a/Core/Algorithms/inc/IOutputDataNormalizer.h b/Core/Algorithms/inc/IOutputDataNormalizer.h
index 0221b67879ebda3ed3b3c361e5f526918c7dd2d3..a456b63db081a16624bbe44e1ab60ee86a73cb5f 100644
--- a/Core/Algorithms/inc/IOutputDataNormalizer.h
+++ b/Core/Algorithms/inc/IOutputDataNormalizer.h
@@ -44,6 +44,8 @@ public:
virtual OutputDataNormalizerScaleAndShift *clone() const;
+ void setMaximumIntensity(double max_intensity) { m_max_intensity = max_intensity; }
+
protected:
//! initialize pool parameters, i.e. register some of class members for later access via parameter pool
virtual void init_parameters();
@@ -54,6 +56,7 @@ private:
double m_scale;
double m_shift;
+ double m_max_intensity;
};
diff --git a/Core/Algorithms/src/ChiSquaredModule.cpp b/Core/Algorithms/src/ChiSquaredModule.cpp
index 71005602b55b61829603a7add9730298b33aca95..bf2764c33d985c1fd55ce6bad7e9b29a7a92573f 100644
--- a/Core/Algorithms/src/ChiSquaredModule.cpp
+++ b/Core/Algorithms/src/ChiSquaredModule.cpp
@@ -23,35 +23,40 @@ double ChiSquaredModule::calculateChiSquared()
if( !mp_real_data ) throw NullPointerException("ChiSquaredModule::calculateChiSquared() -> Error! No real data has been set");
if( !mp_simulation_data ) throw NullPointerException("ChiSquaredModule::calculateChiSquared() -> Error! No simulated data has been set");
- double result = 0.0;
- size_t data_size = mp_real_data->getAllocatedSize();
initWeights();
- if( mp_intensity_function ) {
- OutputDataFunctions::applyFunction(*mp_simulation_data, mp_intensity_function);
- OutputDataFunctions::applyFunction(*mp_real_data, mp_intensity_function);
- }
-
if(mp_data_normalizer) {
OutputData<double > *normalized_simulation = mp_data_normalizer->createNormalizedData(*mp_simulation_data);
delete mp_simulation_data;
mp_simulation_data = normalized_simulation;
}
+ if( mp_intensity_function ) {
+ OutputDataFunctions::applyFunction(*mp_simulation_data, mp_intensity_function);
+ OutputDataFunctions::applyFunction(*mp_real_data, mp_intensity_function);
+ }
+
OutputData<double> *p_difference = createChi2DifferenceMap();
OutputData<double>::const_iterator it_weights = mp_weights->begin();
OutputData<double>::const_iterator it_diff = p_difference->begin();
+ double result(0);
while(it_diff != p_difference->end()) {
result += (*it_diff++)*(*it_weights++);
}
delete p_difference;
- m_chi2_value = result/data_size;
+
+ double fnorm(0);
+ m_ndegree_of_freedom > 0 ? fnorm=double(m_ndegree_of_freedom) : fnorm = double(mp_real_data->getAllocatedSize());
+ m_chi2_value = result/fnorm;
return m_chi2_value;
}
OutputData<double>* ChiSquaredModule::createChi2DifferenceMap() const
{
+ if( !mp_real_data ) throw NullPointerException("ChiSquaredModule::calculateChiSquared() -> Error! No real data has been set");
+ if( !mp_simulation_data ) throw NullPointerException("ChiSquaredModule::calculateChiSquared() -> Error! No simulated data has been set");
+
OutputData<double > *p_difference = mp_simulation_data->clone();
p_difference->setAllTo(0.0);
diff --git a/Core/Algorithms/src/Experiment.cpp b/Core/Algorithms/src/Experiment.cpp
index 9196375c2243389a97a17c0c431ecb5dcd48a232..7e0b0063697b89998e9f8fa6b3bfdc3ddb093f60 100644
--- a/Core/Algorithms/src/Experiment.cpp
+++ b/Core/Algorithms/src/Experiment.cpp
@@ -181,7 +181,7 @@ void Experiment::updateSample()
void Experiment::setDetectorParameters(const OutputData<double > &output_data)
{
-// std::cout << "Experiment::setDetectorParameters() -> Info. Adjusting detector to have shape as in given output data" << std::endl;
+ //std::cout << "Experiment::setDetectorParameters() -> Info. Adjusting detector to have shape as in given output data" << std::endl;
m_detector.clear();
for(size_t i_axis=0; i_axis<output_data.getNdimensions(); ++i_axis) {
const NamedVector<double> *axis = reinterpret_cast<const NamedVector<double>*>(output_data.getAxes()[i_axis]);
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index 5209c36d9e8fd08567cbd5697cb41735c3e745da..de8c12ffa50e8137344fd9e10f78eb758d16a1ed 100644
--- a/Core/Algorithms/src/GISASExperiment.cpp
+++ b/Core/Algorithms/src/GISASExperiment.cpp
@@ -237,9 +237,12 @@ double GISASExperiment::getSolidAngle(size_t index) const
} else {
//std::cout << "GISASExperiment::getSolidAngle() -> Warning! Only one bin on phi_f axis, size of the bin will be taken from alpha_f axis" << std::endl;
}
- if(!dalpha) dalpha=dphi;
- if(!dphi) dphi=dalpha;
- return cos_alpha_f*dalpha*dphi;
+ if(!dalpha || !dphi) {
+ std::cout << "GISASExperiment::getSolidAngle() -> Warning! Not defined normalization" << std::endl;
+ return 1;
+ } else {
+ return cos_alpha_f*dalpha*dphi;
+ }
}
diff --git a/Core/Algorithms/src/IChiSquaredModule.cpp b/Core/Algorithms/src/IChiSquaredModule.cpp
index 592f83d7c4a0fa0e7e81a65c01efdf21baa1d8b1..112d3db7e5a9b177e3b52406c7f5a21058175fc0 100644
--- a/Core/Algorithms/src/IChiSquaredModule.cpp
+++ b/Core/Algorithms/src/IChiSquaredModule.cpp
@@ -9,6 +9,7 @@ IChiSquaredModule::IChiSquaredModule()
, mp_data_selector(0)
, mp_data_normalizer(0)
, mp_intensity_function(0)
+ , m_ndegree_of_freedom(0)
, m_chi2_value(0)
{
mp_squared_function = new SquaredFunctionDefault();
@@ -24,15 +25,17 @@ IChiSquaredModule::IChiSquaredModule(const IChiSquaredModule &other)
, mp_data_selector(0)
, mp_data_normalizer(0)
, mp_intensity_function(0)
+ , m_ndegree_of_freedom(0)
, m_chi2_value(0)
{
- if(other.mp_real_data) mp_real_data = other.mp_real_data;
- if(other.mp_simulation_data) mp_simulation_data = other.mp_simulation_data;
+ if(other.mp_real_data) mp_real_data = other.mp_real_data->clone();
+ if(other.mp_simulation_data) mp_simulation_data = other.mp_simulation_data->clone();
if(other.mp_weights) mp_weights = other.mp_weights->clone();
if(other.mp_squared_function) mp_squared_function = other.mp_squared_function->clone();
if(other.mp_data_selector) mp_data_selector = other.mp_data_selector->clone();
if(other.mp_data_normalizer) mp_data_normalizer = other.mp_data_normalizer->clone();
if(other.mp_intensity_function) mp_intensity_function = other.mp_intensity_function->clone();
+ m_ndegree_of_freedom = other.m_ndegree_of_freedom;
}
diff --git a/Core/Algorithms/src/IOutputDataNormalizer.cpp b/Core/Algorithms/src/IOutputDataNormalizer.cpp
index e5bbd0ed7e29f4022098520b0f0e946ee553601f..6ad14b8f60148f0d0f455d94d7616cae70758951 100644
--- a/Core/Algorithms/src/IOutputDataNormalizer.cpp
+++ b/Core/Algorithms/src/IOutputDataNormalizer.cpp
@@ -6,6 +6,7 @@
OutputDataNormalizerScaleAndShift::OutputDataNormalizerScaleAndShift()
: m_scale(1.0)
, m_shift(0)
+ , m_max_intensity(0)
{
setName("Normalizer");
init_parameters();
@@ -14,6 +15,7 @@ OutputDataNormalizerScaleAndShift::OutputDataNormalizerScaleAndShift()
OutputDataNormalizerScaleAndShift::OutputDataNormalizerScaleAndShift(double scale, double shift)
: m_scale(scale)
, m_shift(shift)
+ , m_max_intensity(0)
{
setName("Normalizer");
init_parameters();
@@ -23,6 +25,7 @@ OutputDataNormalizerScaleAndShift::OutputDataNormalizerScaleAndShift(const Outpu
{
m_scale = other.m_scale;
m_shift = other.m_shift;
+ m_max_intensity = other.m_max_intensity;
init_parameters();
}
@@ -45,17 +48,17 @@ OutputData<double> *OutputDataNormalizerScaleAndShift::createNormalizedData(cons
{
OutputData<double > *normalized_data = data.clone();
- OutputData<double >::const_iterator cit = std::max_element(data.begin(), data.end());
- double max_intensity = (*cit);
- if(max_intensity) {
- OutputData<double >::iterator it = normalized_data->begin();
- while(it!=normalized_data->end()) {
- double value = (*it);
- (*it) = m_scale*(value/max_intensity) + m_shift;
- ++it;
- }
- } else {
- std::cout << "OutputDataNormalizerScaleAndShift::createNormalizedData() -> Warning! Zero maximum intensity" << std::endl;
+ double max_intensity = m_max_intensity;
+ if( max_intensity == 0 ) {
+ OutputData<double >::const_iterator cit = std::max_element(data.begin(), data.end());
+ max_intensity = (*cit);
+ }
+ if(max_intensity == 0) throw LogicErrorException("OutputDataNormalizerScaleAndShift::createNormalizedData() -> Error! Zero maximum intensity");
+ OutputData<double >::iterator it = normalized_data->begin();
+ while(it!=normalized_data->end()) {
+ double value = (*it);
+ (*it) = m_scale*(value/max_intensity) + m_shift;
+ ++it;
}
return normalized_data;
diff --git a/Core/Core.pro b/Core/Core.pro
index dc0a46e0e385f179ef6b046ebc52db1c263d3e80..a09bd6f7e305b26376250bed1377560be9ce0472 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -5,7 +5,7 @@ TARGET = ScattCore
TEMPLATE = lib
CONFIG += plugin # to remove versions from file name
QT -= core gui
-CONFIG += BUILD_PYTHON_BOOST_MODULE # to generate python interface
+#CONFIG += BUILD_PYTHON_BOOST_MODULE # to generate python interface
# including common project properties
include($$PWD/../shared.pri)
diff --git a/Core/Tools/inc/FitObject.h b/Core/Tools/inc/FitObject.h
index 054082cd91da312a0263151231207f64a063652c..16516287facbded668cca87aac964a8ee7b01d20 100644
--- a/Core/Tools/inc/FitObject.h
+++ b/Core/Tools/inc/FitObject.h
@@ -28,7 +28,7 @@
class FitObject : public IParameterized
{
public:
- FitObject(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module=ChiSquaredModule());
+ FitObject(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module=ChiSquaredModule(), double weight = 1.0);
~FitObject();
//! get experiment
@@ -57,6 +57,9 @@ public:
//! add parameters from local pool to external pool and call recursion over direct children
virtual std::string addParametersToExternalPool(std::string path, ParameterPool *external_pool, int copy_number=-1) const;
+ //! return weight of data set in chi2 calculations
+ double getWeight() const { return m_weight; }
+
protected:
//! initialize pool parameters, i.e. register some of class members for later access via parameter pool
virtual void init_parameters();
@@ -68,6 +71,7 @@ private:
Experiment *m_experiment; //! external experiment (not owned by this)
OutputData<double > *m_real_data; //! real data
IChiSquaredModule *m_chi2_module; //! chi2 module
+ double m_weight; //! weight of data set in chi2 calculations
};
diff --git a/Core/Tools/inc/FitSuiteObjects.h b/Core/Tools/inc/FitSuiteObjects.h
index 88e7c526629b47b979d8b5f2d88df8448247374c..02a648ad80c665f7862d8591d246ce94e87e4fd8 100644
--- a/Core/Tools/inc/FitSuiteObjects.h
+++ b/Core/Tools/inc/FitSuiteObjects.h
@@ -42,13 +42,13 @@ public:
size_t size() const { return m_fit_objects.size(); }
//! add to kit pair of (experiment, real data) for consecutive simulation and chi2 module
- void add(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module);
+ void add(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module, double weight = 1.0);
//! loop through all defined experiments and run they simulation
void runSimulation();
- //! get chi squared value calculated for all pairs of (experiment, real data)
- double getChiSquaredValue();
+ //! get sum of chi squared values for all fit objects
+ double getChiSquaredValue(int n_free_fit_parameters = 0);
//! get experiment
const Experiment *getExperiment(size_t i_item = 0) const { return m_fit_objects[check_index(i_item)]->getExperiment(); }
@@ -72,10 +72,16 @@ public:
//! add parameters from local pool to external pool and call recursion over direct children
virtual std::string addParametersToExternalPool(std::string path, ParameterPool *external_pool, int copy_number=-1) const;
+ //! set experiment normalize flag
+ void setExperimentNormalize(bool experiment_normalize) { m_experiment_normalize = experiment_normalize; }
+
protected:
//! initialize pool parameters, i.e. register some of class members for later access via parameter pool
virtual void init_parameters();
+ //! calculate maximum intensity in simulated data
+ double getSimulationMaxIntensity();
+
private:
//! disabled copy constructor and assignment operator
FitSuiteObjects &operator=(const FitSuiteObjects &);
@@ -84,7 +90,9 @@ private:
//! check if index inside vector bounds
inline size_t check_index(size_t index) const { return index < m_fit_objects.size() ? index : throw OutOfBoundsException("FitSuiteKit::check() -> Index outside of range"); }
- FitObjects_t m_fit_objects; // set of experiments and corresponding real data
+ FitObjects_t m_fit_objects; //! set of experiments and corresponding real data
+ double m_total_weight; //! sum of weights of fit sets
+ bool m_experiment_normalize;
};
diff --git a/Core/Tools/inc/FitSuiteParameters.h b/Core/Tools/inc/FitSuiteParameters.h
index f1b90aee5c87a8755720bce251394bfbe454afa3..6a0a3711da7e963b1132c2bd46d64d66f4a9ef48 100644
--- a/Core/Tools/inc/FitSuiteParameters.h
+++ b/Core/Tools/inc/FitSuiteParameters.h
@@ -72,6 +72,9 @@ public:
//! linking fit parameters with pool parameters
void link_to_pool(const ParameterPool *pool);
+ //! return number of free parameters
+ size_t getNfreeParameters() const;
+
private:
//! disabled copy constructor and assignment operator
FitSuiteParameters &operator=(const FitSuiteParameters &other);
diff --git a/Core/Tools/inc/IMinimizer.h b/Core/Tools/inc/IMinimizer.h
index a4015075d3671b786b004b7824e5d3642a21e6b8..f22fd3b63ae86449c9abc5e96b749ee7cb53b0be 100644
--- a/Core/Tools/inc/IMinimizer.h
+++ b/Core/Tools/inc/IMinimizer.h
@@ -17,6 +17,8 @@
#include "FitParameter.h"
#include <boost/function.hpp>
+#include <map>
+#include "Exceptions.h"
//- -------------------------------------------------------------------
@@ -58,4 +60,69 @@ public:
};
+
+
+//- -------------------------------------------------------------------
+//! @class TestMinimizer
+//! @brief Minimizer which calls minimization function once to test whole chain
+//- -------------------------------------------------------------------
+class TestMinimizer : public IMinimizer
+{
+public:
+ TestMinimizer(){}
+ virtual ~TestMinimizer(){}
+
+ //! set variable
+ virtual void setVariable(int index, const FitParameter *par) { m_values[index] = par->getValue(); }
+
+ //! set function to minimize
+ virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1) { m_fcn = fcn, m_ndim = ndim; }
+
+ //! run minimization
+ virtual void minimize()
+ {
+ std::vector<double > buffer;
+ buffer.resize(m_values.size(), 0.0);
+ for(std::map<int, double >::iterator it=m_values.begin(); it!= m_values.end(); ++it ) {
+ buffer[it->first] = it->second;
+ std::cout << " minimize(): " << it->first << " " << it->second << std::endl;
+ }
+ std::cout << "TestMinimizer::minimize() -> Info. Calling fcn" << std::endl;
+ m_fcn(&buffer[0]);
+ }
+
+ //! get number of variables to fit
+ virtual size_t getNumberOfVariables() const { return m_values.size(); }
+
+ //! return minimum function value
+ virtual double getMinValue() const { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); return 0.0; }
+
+ //! return pointer to the parameters values at the minimum
+ virtual double getValueOfVariableAtMinimum(size_t i) const
+ {
+ std::map<int, double >::const_iterator pos = m_values.find(i);
+ if(pos != m_values.end()){
+ return pos->second;
+ } else {
+ throw LogicErrorException("TestMinimizer::getValueOfVariableAtMinimum() -> Not found!");
+ }
+ }
+
+ //! return pointer to the parameters values at the minimum
+ virtual double getErrorOfVariable(size_t /* i*/) const { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); return 0.0; }
+
+ //! clear resources (parameters) for consecutives minimizations
+ virtual void clear() { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); }
+
+ //! print fit results
+ virtual void printResults() const { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); }
+
+private:
+ std::map<int, double > m_values;
+ boost::function<double(const double *)> m_fcn;
+ int m_ndim;
+};
+
+
+
#endif // IMINIMIZER_H
diff --git a/Core/Tools/src/Exceptions.cpp b/Core/Tools/src/Exceptions.cpp
index 2a8de5a244885285a27ce84ac13ecd7a44c6a867..f9a68eea06953ff2bed20bfa24e8c9f85a8dbc14 100644
--- a/Core/Tools/src/Exceptions.cpp
+++ b/Core/Tools/src/Exceptions.cpp
@@ -5,7 +5,7 @@ namespace Exceptions {
void LogExceptionMessage(const std::string &message)
{
- std::cout << message << std::endl;
+ std::cerr << message << std::endl;
}
NotImplementedException::NotImplementedException(const std::string &message)
diff --git a/Core/Tools/src/FitObject.cpp b/Core/Tools/src/FitObject.cpp
index ee0335df06554e6aefbd1ea48b2710e5d30049ac..1ceb4844f05156705f496c949ed480929df16a46 100644
--- a/Core/Tools/src/FitObject.cpp
+++ b/Core/Tools/src/FitObject.cpp
@@ -5,10 +5,11 @@
/* ************************************************************************* */
// FitObject c-tors
/* ************************************************************************* */
-FitObject::FitObject(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module)
+FitObject::FitObject(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module, double weight)
: m_experiment(experiment.clone())
, m_real_data(real_data.clone())
, m_chi2_module(chi2_module.clone())
+ , m_weight(weight)
{
setName("FitObject");
if( !m_real_data->hasSameShape(*m_experiment->getOutputData()) ) {
@@ -33,9 +34,9 @@ void FitObject::setRealData(const OutputData<double > &real_data)
m_real_data = real_data.clone();
if( m_experiment) {
if( !m_real_data->hasSameShape(*m_experiment->getOutputData()) ) {
- std::cout << "FitSuiteKit::KitItem::setRealData() -> Real data and the detector have different shape. Adjusting detector..." << std::endl;
+ std::cout << "FitObject::setRealData() -> Real data and the detector have different shape. Adjusting detector..." << std::endl;
} else {
- std::cout << "FitSuiteKit::KitItem::setRealData() -> Real data and the detector have same shape. No nedd to adjust detector." << std::endl;
+ std::cout << "FitObject::setRealData() -> Real data and the detector have same shape. No need to adjust detector." << std::endl;
}
m_experiment->setDetectorParameters(*m_real_data);
}
diff --git a/Core/Tools/src/FitSuite.cpp b/Core/Tools/src/FitSuite.cpp
index b6e93cebb8bf1b9c2afaae887cb2e22c00b24232..e08ec9cc2abf49064404dc1ce93a277e6f6c05ae 100644
--- a/Core/Tools/src/FitSuite.cpp
+++ b/Core/Tools/src/FitSuite.cpp
@@ -70,7 +70,7 @@ void FitSuite::link_fit_parameters()
{
ParameterPool *pool = m_fit_objects.createParameterTree();
m_fit_parameters.link_to_pool(pool);
- std::cout << "XXXXXX FitSuite::link_fit_parameters() -> " << std::endl;
+ std::cout << "FitSuite::link_fit_parameters() -> Parameter pool:" << std::endl;
std::cout << *pool << std::endl;
std::cout << "----------------" << std::endl;
delete pool;
@@ -155,7 +155,7 @@ double FitSuite::functionToMinimize(const double *pars_current_values)
m_fit_objects.runSimulation();
// caclulate chi2 value
- double chi_squared = m_fit_objects.getChiSquaredValue();
+ double chi_squared = m_fit_objects.getChiSquaredValue(m_fit_parameters.getNfreeParameters());
notifyObservers();
m_n_call++;
diff --git a/Core/Tools/src/FitSuiteObjects.cpp b/Core/Tools/src/FitSuiteObjects.cpp
index 99916ff302de11cfa69ebb9c562bf19ffb43e90a..2e1493edd87cd0bdba43b1c773cb5092e9d5893f 100644
--- a/Core/Tools/src/FitSuiteObjects.cpp
+++ b/Core/Tools/src/FitSuiteObjects.cpp
@@ -2,7 +2,7 @@
-FitSuiteObjects::FitSuiteObjects()
+FitSuiteObjects::FitSuiteObjects() : m_total_weight(0), m_experiment_normalize(false)
{
setName("FitSuiteObjects");
init_parameters();
@@ -22,9 +22,10 @@ void FitSuiteObjects::clear()
/* ************************************************************************* */
// add to kit pair of (experiment, real data) for consecutive simulation and chi2 module
/* ************************************************************************* */
-void FitSuiteObjects::add(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module)
+void FitSuiteObjects::add(const Experiment &experiment, const OutputData<double > &real_data, const IChiSquaredModule &chi2_module, double weight)
{
- m_fit_objects.push_back(new FitObject(experiment, real_data, chi2_module));
+ m_total_weight += weight;
+ m_fit_objects.push_back(new FitObject(experiment, real_data, chi2_module, weight));
}
@@ -35,21 +36,47 @@ void FitSuiteObjects::runSimulation()
{
for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
(*it)->getExperiment()->runSimulation();
- (*it)->getExperiment()->normalize();
+ if(m_experiment_normalize) (*it)->getExperiment()->normalize();
}
}
/* ************************************************************************* */
-// loop through all defined experiments and run they simulation
+// get sum of chi squared values for all fit objects
+// FIXME: refactor FitSuiteObjects::getChiSquaredValue()
+/* ************************************************************************* */
+double FitSuiteObjects::getChiSquaredValue(int n_free_fit_parameters)
+{
+ double chi_sum(0);
+ for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
+ IChiSquaredModule *chi = (*it)->getChiSquaredModule();
+
+ chi->setNdegreeOfFreedom( m_fit_objects.size() * (*it)->getRealData()->getAllocatedSize() - n_free_fit_parameters);
+ // normalizing datasets to the maximum intensity over all fit objects defined
+ OutputDataNormalizerScaleAndShift *data_normalizer = dynamic_cast<OutputDataNormalizerScaleAndShift *>(chi->getOutputDataNormalizer());
+ if( data_normalizer) data_normalizer->setMaximumIntensity( getSimulationMaxIntensity() );
+
+ double weight = (*it)->getWeight()/m_total_weight;
+ double chi_squared = (weight*weight) * (*it)->calculateChiSquared();
+ chi_sum += chi_squared;
+// std::cout << " chi " << chi_squared << " chi_sum:" << chi_sum << std::endl;
+ }
+ return chi_sum;
+}
+
+
+/* ************************************************************************* */
+// calculate maximum intensity in simulated data over all fir objects defined
/* ************************************************************************* */
-double FitSuiteObjects::getChiSquaredValue()
+double FitSuiteObjects::getSimulationMaxIntensity()
{
- double chi_squared(0);
+ double max_intensity(0);
for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
- chi_squared += (*it)->calculateChiSquared();
+ const OutputData<double > *data = (*it)->getExperiment()->getOutputData();
+ OutputData<double >::const_iterator cit = std::max_element(data->begin(), data->end());
+ max_intensity = std::max(max_intensity, *cit);
}
- return chi_squared;
+ return max_intensity;
}
diff --git a/Core/Tools/src/FitSuiteParameters.cpp b/Core/Tools/src/FitSuiteParameters.cpp
index 404e3afd031b6da0404f2964fc959bd9d708464b..7543a1f52a9c2403216ce960f1eca4a490e8dee7 100644
--- a/Core/Tools/src/FitSuiteParameters.cpp
+++ b/Core/Tools/src/FitSuiteParameters.cpp
@@ -79,6 +79,17 @@ std::vector<double > FitSuiteParameters::getValues() const
}
+size_t FitSuiteParameters::getNfreeParameters() const
+{
+ size_t n_free(0);
+ for(parameters_t::const_iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it)
+ {
+ if( !(*it)->isFixed() ) n_free++;
+ }
+ return n_free;
+}
+
+
/* ************************************************************************* */
// linking fit parameters with pool parameters
diff --git a/Core/Tools/src/OutputDataReader.cpp b/Core/Tools/src/OutputDataReader.cpp
index 98cbaa8ff15c3e6e40f01518d7e244f6b3df1d5f..64de72be79953095e9d95db9c9a1fa7b5854d147 100644
--- a/Core/Tools/src/OutputDataReader.cpp
+++ b/Core/Tools/src/OutputDataReader.cpp
@@ -14,7 +14,12 @@
#include <boost/iostreams/filtering_streambuf.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/copy.hpp>
+
+#include "Macros.h"
+GCC_DIAG_OFF(unused-parameter);
#include <boost/iostreams/filter/gzip.hpp>
+GCC_DIAG_ON(unused-parameter);
+
OutputData<double > *OutputDataReader::getOutputData(const std::string &file_name)
diff --git a/UnitTests/TestCore/ChiSquaredModuleTest.h b/UnitTests/TestCore/ChiSquaredModuleTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..6c7cb3a82aeeeb8e0006afe87daec7c0e256af10
--- /dev/null
+++ b/UnitTests/TestCore/ChiSquaredModuleTest.h
@@ -0,0 +1,126 @@
+#ifndef CHISQUAREDMODULETEST_H
+#define CHISQUAREDMODULETEST_H
+
+#include "ChiSquaredModule.h"
+#include "IFittingDataSelector.h"
+#include "ISquaredFunction.h"
+
+#include "gtest/gtest.h"
+
+
+class ChiSquaredModuleTest : public ::testing::Test
+{
+protected:
+ ChiSquaredModuleTest();
+ virtual ~ChiSquaredModuleTest();
+
+ ChiSquaredModule m_chi_empty;
+ ChiSquaredModule m_chi_default;
+ OutputData<double > m_real_data;
+ OutputData<double > m_simul_data;
+};
+
+
+ChiSquaredModuleTest::ChiSquaredModuleTest()
+{
+ m_real_data.addAxis(new NamedVector<double>("phi_f", 0.0, 10.0, 10));
+ m_real_data.addAxis(new NamedVector<double>("alpha_f", 0.0, 10.0, 10));
+ m_real_data.setAllTo(1.0);
+ m_simul_data.copyFrom(m_real_data);
+ m_simul_data.setAllTo(1.1);
+
+ // default module with set of real and simulated data
+ m_chi_default.setRealAndSimulatedData(m_real_data, m_simul_data);
+
+}
+
+
+ChiSquaredModuleTest::~ChiSquaredModuleTest()
+{
+
+}
+
+
+TEST_F(ChiSquaredModuleTest, InitialState)
+{
+ EXPECT_EQ( NULL, m_chi_empty.getRealData());
+ EXPECT_EQ( NULL, m_chi_empty.getSimulationData());
+ EXPECT_TRUE( dynamic_cast<const SquaredFunctionDefault *>(m_chi_empty.getSquaredFunction()));
+ EXPECT_TRUE( dynamic_cast<const DefaultAllDataSelector *>(m_chi_empty.getFittingDataSelector()));
+ EXPECT_EQ( NULL, m_chi_empty.getOutputDataNormalizer());
+ EXPECT_EQ( NULL, m_chi_empty.getIntensityFunction());
+ EXPECT_EQ( double(0), m_chi_empty.getValue());
+ ASSERT_THROW(m_chi_empty.calculateChiSquared(), NullPointerException);
+ ASSERT_THROW(m_chi_empty.createChi2DifferenceMap(), NullPointerException);
+}
+
+TEST_F(ChiSquaredModuleTest, CloneOfEmpty)
+{
+ ChiSquaredModule *clone_of_empty = m_chi_empty.clone();
+ EXPECT_EQ( NULL, clone_of_empty->getRealData());
+ EXPECT_EQ( NULL, clone_of_empty->getSimulationData());
+ EXPECT_TRUE( dynamic_cast<const SquaredFunctionDefault *>(clone_of_empty->getSquaredFunction()));
+ EXPECT_TRUE( dynamic_cast<const DefaultAllDataSelector *>(clone_of_empty->getFittingDataSelector()));
+ EXPECT_EQ( NULL, clone_of_empty->getOutputDataNormalizer());
+ EXPECT_EQ( NULL, clone_of_empty->getIntensityFunction());
+ EXPECT_EQ( double(0), clone_of_empty->getValue());
+ ASSERT_THROW(clone_of_empty->calculateChiSquared(), NullPointerException);
+ ASSERT_THROW(clone_of_empty->createChi2DifferenceMap(), NullPointerException);
+ delete clone_of_empty;
+}
+
+TEST_F(ChiSquaredModuleTest, DefaultModule)
+{
+ EXPECT_FLOAT_EQ( double(0.01), m_chi_default.calculateChiSquared());
+ EXPECT_FLOAT_EQ( double(0.01), m_chi_default.getValue());
+}
+
+TEST_F(ChiSquaredModuleTest, CloneOfDefault)
+{
+ ChiSquaredModule *clone_of_default = m_chi_default.clone();
+ EXPECT_FLOAT_EQ( double(0.01), clone_of_default->calculateChiSquared());
+ EXPECT_FLOAT_EQ( double(0.01), clone_of_default->getValue());
+ clone_of_default->setNdegreeOfFreedom(1);
+ EXPECT_FLOAT_EQ( double(1.0), clone_of_default->calculateChiSquared());
+ delete clone_of_default;
+}
+
+
+TEST_F(ChiSquaredModuleTest, IsGISAXSLikeModule)
+{
+ ChiSquaredModule chi_isgisaxs;
+ OutputData<double > real_data;
+ OutputData<double > simul_data;
+ const size_t nbins(5);
+ real_data.addAxis(new NamedVector<double>("phi_f", 0.0, 1.0, nbins));
+ simul_data.addAxis(new NamedVector<double>("phi_f", 0.0, 1.0, nbins));
+ const double a_real_data[nbins] = {1., 10., 100., 10., 1. };
+ const double a_simul_data[nbins] = {10., 100., 1000., 100., 10. };
+ OutputData<double >::iterator it_real = real_data.begin();
+ OutputData<double >::iterator it_simul = simul_data.begin();
+ int index(0);
+ while(it_real != real_data.end()) {
+ *it_real = a_real_data[index];
+ *it_simul = a_simul_data[index];
+ ++index; ++it_real; ++it_simul;
+ }
+ chi_isgisaxs.setRealAndSimulatedData(real_data, simul_data);
+ OutputDataNormalizerScaleAndShift normalizer(100., 0.0);
+ chi_isgisaxs.setOutputDataNormalizer( normalizer );
+ EXPECT_FLOAT_EQ( double(0.0), chi_isgisaxs.calculateChiSquared());
+
+// m_chi_isgisaxs.setChiSquaredFunction( SquaredFunctionWithSystematicError(1.0) );
+// EXPECT_FLOAT_EQ( double(0.005), m_chi_isgisaxs.calculateChiSquared());
+
+// m_chi_isgisaxs.setOutputDataNormalizer( OutputDataNormalizerScaleAndShift() );
+// EXPECT_FLOAT_EQ( double(0.005), m_chi_isgisaxs.calculateChiSquared());
+
+
+}
+
+
+#endif // CHISQUAREDMODULETEST_H
+
+
+
+
diff --git a/UnitTests/TestCore/TestCore.pro b/UnitTests/TestCore/TestCore.pro
index 2ab675f7229c86f3cfa27a3a3f095fdfbdf454fc..f8677b63241d0ede8d53de32c41a8d2a24526e8d 100644
--- a/UnitTests/TestCore/TestCore.pro
+++ b/UnitTests/TestCore/TestCore.pro
@@ -19,7 +19,8 @@ HEADERS += \
NamedVectorTest.h \
OutputDataTest.h \
OutputDataIteratorTest.h \
- GISASExperimentTest.h
+ GISASExperimentTest.h \
+ ChiSquaredModuleTest.h
OBJECTS_DIR = obj
@@ -68,4 +69,4 @@ for(dep, MY_DEPENDENCY_LIB) {
###############################################################################
# runs automatically tests right after linking
###############################################################################
-QMAKE_POST_LINK = ./$(TARGET)
+QMAKE_POST_LINK = ./$(TARGET) 2> /dev/null
diff --git a/UnitTests/TestCore/main.cpp b/UnitTests/TestCore/main.cpp
index 10af4c472ed6cab820bcc8bdad023c139f5683f2..17c20c1cc91f6c5b196c96f90665c19bb398a082 100644
--- a/UnitTests/TestCore/main.cpp
+++ b/UnitTests/TestCore/main.cpp
@@ -1,6 +1,7 @@
#include "gtest/gtest.h"
#include "BeamTest.h"
+#include "ChiSquaredModuleTest.h"
#include "DetectorTest.h"
#include "ExperimentTest.h"
#include "GISASExperimentTest.h"