diff --git a/App/inc/TestIsGISAXS12.h b/App/inc/TestIsGISAXS12.h
index c4a5c010bd9d3d4b5b21531f03e7a1f243612025..a7db5cc02005af6ff2df625a29bc2715e5d55be2 100644
--- a/App/inc/TestIsGISAXS12.h
+++ b/App/inc/TestIsGISAXS12.h
@@ -51,10 +51,11 @@ 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_dispersion_radius1;
double m_height_aspect_ratio1;
+ double m_particle_probability2;
double m_particle_radius2;
double m_dispersion_radius2;
double m_height_aspect_ratio2;
@@ -125,8 +126,8 @@ private:
//! convert isgisaxs 1d scan to output data 2d object
OutputData<double> *convert_isgi_scan(std::vector<IsgiData > &isgi_data);
- //! run test fit
- void run_test_fit();
+ //! run test minimizer to check whole chain
+ void run_test_minimizer();
void print_axes(DataSet &data);
diff --git a/App/src/DrawHelper.cpp b/App/src/DrawHelper.cpp
index 4c189752f1dc185334c1c8773bce6daa6e2e839b..8b8387ab48ddef7c568fbdec332175de406a530f 100644
--- a/App/src/DrawHelper.cpp
+++ b/App/src/DrawHelper.cpp
@@ -62,7 +62,7 @@ void DrawHelper::SetMagnifier(TCanvas *canvas)
void DrawHelper::ExecuteMagnifier(int event, int px, int py, TObject *sel)
{
(void)sel;
- if ( (event == kButton1Double) ) {
+ if ( event == kButton1Double ) {
TCanvas *c = (TCanvas *) gTQSender;
char cname[256];
sprintf(cname,"%s_%d",c->GetTitle(),(int)time(0));
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/TestIsGISAXS12.cpp b/App/src/TestIsGISAXS12.cpp
index 976df02a3934c6fe962a607f6723e4510de2617b..5713c5f8db21d6a57fe2b009b43d8bb88958fdcd 100644
--- a/App/src/TestIsGISAXS12.cpp
+++ b/App/src/TestIsGISAXS12.cpp
@@ -36,6 +36,7 @@
#include "TH2D.h"
#include "TLine.h"
#include "TROOT.h"
+#include "TLegend.h"
@@ -77,11 +78,12 @@ void TestIsGISAXS12::execute()
// plot isgisaxs data
plot_isgisaxs_data();
+ // run test fit
+ run_test_minimizer();
+
// run isgisaxs ex-12 style fit
- //run_isgisaxs_fit();
+ run_isgisaxs_fit();
- // run test fit
- //run_test_fit();
}
@@ -94,7 +96,6 @@ void TestIsGISAXS12::run_isgisaxs_comparison()
{
// run simulation for default sample parameters
m_experiment->runSimulation();
- m_experiment->normalize();
IsGISAXSTools::writeOutputDataToFile(*(m_experiment->getOutputData()), m_data_path+"this_fitconstraints.ima");
// plotting results of comparison we/isgisaxs for the sample with default parameters
@@ -104,15 +105,10 @@ void TestIsGISAXS12::run_isgisaxs_comparison()
// -------------
// plot results
// -------------
- OutputData<double> *isgi_data_orig = IsGISAXSTools::readOutputDataFromFile(isgi_file);
- OutputData<double> *our_data_orig = IsGISAXSTools::readOutputDataFromFile(this_file);
-
- OutputDataNormalizerScaleAndShift mr(1.0, 0.0);
- OutputData<double> *isgi_data = mr.createNormalizedData(*isgi_data_orig);
- OutputData<double> *our_data = mr.createNormalizedData(*our_data_orig);
+ OutputData<double> *isgi_data = IsGISAXSTools::readOutputDataFromFile(isgi_file);
+ 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.);
@@ -127,8 +123,8 @@ void TestIsGISAXS12::run_isgisaxs_comparison()
// difference
c1->cd(3);
-// IsGISAXSTools::setMinimum(-0.0001);
-// IsGISAXSTools::setMaximum(0.0001);
+ IsGISAXSTools::setMinimum(-0.0001);
+ IsGISAXSTools::setMaximum(0.0001);
IsGISAXSTools::drawOutputDataRelativeDifference2D(*our_data, *isgi_data, "CONT4 Z", "2D Difference map");
// difference
@@ -137,10 +133,10 @@ void TestIsGISAXS12::run_isgisaxs_comparison()
IsGISAXSTools::setMinimum(1);
IsGISAXSTools::drawOutputDataDifference1D(*our_data, *isgi_data, "", "Difference spectra");
-
delete isgi_data;
delete our_data;
+ c1->Update();
}
@@ -151,13 +147,7 @@ void TestIsGISAXS12::run_isgisaxs_fit()
{
// reading info about two 1D scans defined in isgisaxs example
DataSet isgi_scans;
- read_isgisaxs_datfile(m_data_path+"isgi_fitconstraints.dat", isgi_scans);
-
-
-
- // making some plots together with test simulation
-// plot_isgisaxs_data();
-// return;
+ read_isgisaxs_outfile(m_data_path+"isgi_fitconstraints.out", isgi_scans);
// creating fit suite
m_fitSuite = new FitSuite();
@@ -170,37 +160,38 @@ 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("*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("*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_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("*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/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_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("*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("*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("*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_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("*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/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_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;
@@ -214,22 +205,51 @@ void TestIsGISAXS12::run_isgisaxs_fit()
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(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();
}
@@ -252,161 +272,62 @@ void TestIsGISAXS12::plot_isgisaxs_data()
print_axes(isgi_scans_smoothed);
print_axes(isgi_results);
- TCanvas *c1 = DrawHelper::instance().createAndRegisterCanvas("c1_isgisaxs_data", "c1_isgisaxs_data", 768, 1024);
+ TCanvas *c1 = DrawHelper::instance().createAndRegisterCanvas("c1_isgisaxs_data", "Looking on IsGISAXS data and fit results", 768, 1024);
c1->Divide(2,3);
// 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");
}
+ 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");
}
+ c1->cd(3); leg2->Draw();
+ c1->cd(4); leg2->Draw();
-
- /*
- *
- *
- 1 1 # Global scale factor : 1.31728E+05
- 2 2 # Global shift factor : 5.90612E-01
- 3 3 # ( 1) Particle probability : 6.60819E-01
- 4 4 # ( 1) Particle radius (nm) : 4.92536E+00
- 5 5 # ( 1) Dispersion of radius: 1.86101E-01
- 6 6 # ( 1) Height aspect ratio : 9.84164E-01
- 7 7 # ( 2) Particle probability : 3.44720E-01
- 8 8 # ( 2) Particle radius (nm) : 5.26120E+00
- 9 9 # ( 2) Dispersion of radius: 2.16915E-01
- 10 10 # ( 2) Height aspect ratio : 1.00898E+00
- 11 11 # D (nm) : 1.49763E+01
- 12 12 # w (nm) : 3.02653E+00
- */
-
- // 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.31728E+05, 100, AttLimits::limited(1e4, 2e5));
-// m_fitSuite->addFitParameter("*Normalizer/shift", 5.90612E-01, 1, AttLimits::limited(1., 20.));
-// m_fitSuite->addFitParameter("*SampleBuilder/particle_probability", 6.60819E-01, 0.1, AttLimits::limited(0.01, 1.0) );
-
-// m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 4.92536E+00, 1*Units::nanometer, AttLimits::limited(1., 10.) );
-// m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 1.86101E-01, 0.1, AttLimits::limited(0.01, 1.) );
-// m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 9.84164E-01, 0.1, AttLimits::limited(0.01, 10.) );
-
-// m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 5.26120E+00, 1*Units::nanometer, AttLimits::limited(1., 10.) );
-// m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 2.16915E-01, 0.1, AttLimits::limited(0.01, 1.) );
-// m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 1.00898E+00, 0.1, AttLimits::limited(0.01, 10.) );
-
-// m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 1.49763E+01, 1*Units::nanometer, AttLimits::limited(0.01, 50.0) );
-// m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 3.02653E+00, 1*Units::nanometer, AttLimits::limited(0.01, 10.) );
-
-// // setting up fitSuite
-// ChiSquaredModule chiModule;
-// chiModule.setOutputDataNormalizer( OutputDataNormalizerScaleAndShift() );
-
-// for(DataSet::iterator it=isgi_scans_smoothed.begin(); it!= isgi_scans_smoothed.end(); ++it) {
-// m_fitSuite->addExperimentAndRealData(*m_experiment, *(*it), chiModule);
-// }
-
-// m_fitSuite->runFit();
-
- // drawing isgsaxs fit results on top of isgisaxs real data
+ // drawing ration 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); gPad->SetLogy();
- //TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*isgi_scans_smoothed[i_set], "data");
- TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*isgi_results[i_set], "data");
+ 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();
-// const OutputData<double > *data = m_fitSuite->getFitObjects()->getObject(i_set)->getChiSquaredModule()->getSimulationData();
-// TH1D *hdata_smoothed = IsGISAXSTools::getOutputDataScanHist(*data, "data_from_module");
-// hdata_smoothed->SetLineColor(kBlue);
-// hdata_smoothed->DrawCopy("same");
+ if(i_set==0) leg3->AddEntry(hdata,"isgisaxs results/data","lp");
}
+ c1->cd(5); leg3->Draw();
+ c1->cd(6); leg3->Draw();
+ c1->Update();
-
- m_experiment->setDetectorParameters(*isgi_results[0]);
- m_experiment->runSimulation();
- //m_experiment->normalize();
- c1->cd(5);
- OutputDataNormalizerScaleAndShift mr(1.31728E+05, 5.90612E-01);
- OutputData<double > *xxx = mr.createNormalizedData(*m_experiment->getOutputData());
- TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*xxx, "data");
- hdata->DrawCopy("same");
-
-
-// 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 (phi_f axis from one scan (with alpha_f fixed), alpha_f axis from another scan (with phi_f fixed)
-// 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();
-
-// // 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 ));
-
-// // drawing data and scans
-// TCanvas *c1 = new TCanvas("c1","c1",1024, 768);
-// 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;
-// }
-// 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;
-// }
-
-// // 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->Update();
-
-// for(DataScan_t::iterator it=data_scans.begin(); it!= data_scans.end(); ++it) delete (*it);
}
@@ -482,6 +403,9 @@ void TestIsGISAXS12::read_isgisaxs_datfile(const std::string &filename, DataSet
/* ************************************************************************* */
// 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)
{
@@ -517,7 +441,7 @@ void TestIsGISAXS12::read_isgisaxs_outfile(const std::string &filename, DataSet
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");
+ 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]);
@@ -607,137 +531,110 @@ OutputData<double> *TestIsGISAXS12::convert_isgi_scan(std::vector<IsgiData > &is
++it;
}
-// OutputData<double > *data_half = OutputDataFunctions::doubleBinSize(*data);
-// OutputData<double > *data_quarter = OutputDataFunctions::doubleBinSize(*data_half);
-
-// std::cout << data->getAllocatedSize() << std::endl;
-// std::cout << data_half->getAllocatedSize() << std::endl;
-// std::cout << data_quarter->getAllocatedSize() << std::endl;
-// delete data;
-// delete data_half;
-// TCanvas *c1 = dynamic_cast<TCanvas *>( gROOT->FindObject("c1_data_interpol") );
-// if( !c1) {
-// c1 = new TCanvas("c1_data_interpol","c1_data_interpol",800, 600);
-// c1->Divide(2,2);
-// c1->cd(1);
-// } else {
-// c1->cd(2);
-// }
-// TH1D *hdata = IsGISAXSTools::getOutputDataScanHist(*data,"data");
-// TH1D *hdata_half = IsGISAXSTools::getOutputDataScanHist(*data_half,"data_half");
-// TH1D *hdata_quarter = IsGISAXSTools::getOutputDataScanHist(*data_quarter,"data_quarter");
-// hdata_half->SetLineColor(kBlue);
-// hdata_quarter->SetLineColor(kRed);
-// hdata->Draw();
-// hdata_half->Draw("same");
-// hdata_quarter->Draw("same");
-
return data;
}
/* ************************************************************************* */
-// run test fitting (no isgisaxs comparison here)
+// run test minimizer to check the whole chain
/* ************************************************************************* */
-void TestIsGISAXS12::run_test_fit()
+void TestIsGISAXS12::run_test_minimizer()
{
- // creating "real" data
- //m_experiment->setDetectorResolutionFunction(new ResolutionFunction2DSimple(0.0002, 0.0002));
- m_experiment->setBeamIntensity(1e10);
- m_experiment->runSimulation();
- m_experiment->normalize();
+ // 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);
- OutputData<double > *real_data = IsGISAXSTools::createNoisyData( *m_experiment->getOutputData() );
+ // 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
- // setting up 1d scans by making slices on real data
- DataSet 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) );
+ m_fitSuite = new FitSuite();
+ m_fitSuite->setMinimizer( new TestMinimizer() );
- // 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(DataSet::const_iterator it=data_scans.begin(); it!= data_scans.end(); ++it) {
- TLine *line = IsGISAXSTools::getOutputDataScanLine(*(*it));
- line->DrawClone();
- delete line;
- }
+ 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.));
- int npad(2);
- for(DataSet::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();
+ 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.) );
- // creating fit suite
- m_fitSuite = new FitSuite();
- m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
-// m_fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiFit", "") );
- 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(DataSet::iterator it=data_scans.begin(); it!= data_scans.end(); ++it) {
- m_fitSuite->addExperimentAndRealData(*m_experiment, *(*it));
- }
+ 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_scans_smoothed.begin(); it!= isgi_scans_smoothed.end(); ++it) {
+ m_fitSuite->addExperimentAndRealData(*m_experiment, *(*it), chiModule);
+ }
m_fitSuite->runFit();
- delete real_data;
- for(DataSet::iterator it=data_scans.begin(); it!= data_scans.end(); ++it) delete (*it);
-}
+ 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();
+}
/* ************************************************************************* */
//
/* ************************************************************************* */
TestIsGISAXS12::TestSampleBuilder::TestSampleBuilder()
-// starting parameters
-// : m_particle_probability(0.4)
-// , m_particle_radius1(4*Units::nanometer)
-// , m_dispersion_radius1(0.2)
-// , m_height_aspect_ratio1(0.8)
-// , m_particle_radius2(4*Units::nanometer)
-// , m_dispersion_radius2(0.2)
-// , m_height_aspect_ratio2(0.8)
-// , m_interf_distance(12.0*Units::nanometer)
-// , m_interf_width(6*Units::nanometer)
-// optimal parameters
- : m_particle_probability(6.60819E-01)
- , m_particle_radius1(4.92536E+00*Units::nanometer)
- , m_dispersion_radius1(1.86101E-01)
- , m_height_aspect_ratio1(9.84164E-01)
- , m_particle_radius2(5.26120E+00*Units::nanometer)
- , m_dispersion_radius2(2.16915E-01)
- , m_height_aspect_ratio2(1.00898E+00)
- , m_interf_distance(1.49763E+01*Units::nanometer)
- , m_interf_width(3.02653E+00*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();
@@ -746,10 +643,11 @@ 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("dispersion_radius1", &m_dispersion_radius1);
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);
@@ -772,6 +670,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;
@@ -796,9 +698,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
diff --git a/Core/Algorithms/inc/IChiSquaredModule.h b/Core/Algorithms/inc/IChiSquaredModule.h
index 169647ce196b62ce1a25262dfe7af39b42170b9e..5b6299f835d72ff6396f4b626313259d080822a0 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 free parameters
+ virtual void setNfreeParameters(int nfree_parameters) { m_nfree_parameters = nfree_parameters; }
+
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_nfree_parameters;
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..82789b15c9a5bd9420f7916bbefc6f4fbeb411dc 100644
--- a/Core/Algorithms/src/ChiSquaredModule.cpp
+++ b/Core/Algorithms/src/ChiSquaredModule.cpp
@@ -24,20 +24,20 @@ double ChiSquaredModule::calculateChiSquared()
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();
+ size_t data_size = mp_real_data->getAllocatedSize() - m_nfree_parameters;
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();
@@ -45,7 +45,7 @@ double ChiSquaredModule::calculateChiSquared()
result += (*it_diff++)*(*it_weights++);
}
delete p_difference;
- m_chi2_value = result/data_size;
+ m_chi2_value = result/(double)data_size;
return m_chi2_value;
}
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index a21fdabe29a7a3eb9c070db3805499d77b1ae96f..04ae42bec538421eb5c8b5106dcf94c5ad2e6271 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..b9d19768123f47279baff909c4a104623ba0879c 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_nfree_parameters(0)
, m_chi2_value(0)
{
mp_squared_function = new SquaredFunctionDefault();
@@ -24,6 +25,7 @@ IChiSquaredModule::IChiSquaredModule(const IChiSquaredModule &other)
, mp_data_selector(0)
, mp_data_normalizer(0)
, mp_intensity_function(0)
+ , m_nfree_parameters(0)
, m_chi2_value(0)
{
if(other.mp_real_data) mp_real_data = other.mp_real_data;
@@ -33,6 +35,7 @@ IChiSquaredModule::IChiSquaredModule(const IChiSquaredModule &other)
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_nfree_parameters = other.m_nfree_parameters;
}
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/Tools/inc/FitSuiteParameters.h b/Core/Tools/inc/FitSuiteParameters.h
index 860fbc5fb935e48087060906caa7578974c2c0b2..e95b891f303f3d5560d239e9aa9eb54239c3458c 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/src/FitSuite.cpp b/Core/Tools/src/FitSuite.cpp
index 8dbd0206fe5544d28da850cd85955ff6c4afd6a0..b11f9a6c0943c313abe5e2050f823dd687a7ec00 100644
--- a/Core/Tools/src/FitSuite.cpp
+++ b/Core/Tools/src/FitSuite.cpp
@@ -121,6 +121,12 @@ void FitSuite::runFit()
// linking fit parameters with parameters defined in the experiment
link_fit_parameters();
+ // FIXME (find better place) propagating number of free parameters
+ for(size_t i=0; i<m_fit_objects.size(); ++i) {
+ m_fit_objects.getChiSquaredModule(i)->setNfreeParameters(m_fit_parameters.getNfreeParameters());
+ }
+
+
// running minimizer
if( m_fit_strategies.empty() ) {
// running single minimization round
diff --git a/Core/Tools/src/FitSuiteObjects.cpp b/Core/Tools/src/FitSuiteObjects.cpp
index 99916ff302de11cfa69ebb9c562bf19ffb43e90a..c96c2443454be75d0c89d6bf18309203fa28907d 100644
--- a/Core/Tools/src/FitSuiteObjects.cpp
+++ b/Core/Tools/src/FitSuiteObjects.cpp
@@ -35,7 +35,7 @@ void FitSuiteObjects::runSimulation()
{
for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
(*it)->getExperiment()->runSimulation();
- (*it)->getExperiment()->normalize();
+ //(*it)->getExperiment()->normalize();
}
}
@@ -45,11 +45,27 @@ void FitSuiteObjects::runSimulation()
/* ************************************************************************* */
double FitSuiteObjects::getChiSquaredValue()
{
- double chi_squared(0);
+ // determining maximum intensity in all datasets
+ 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;
+
+ double chi_sum(0);
+ for(FitObjects_t::iterator it = m_fit_objects.begin(); it!= m_fit_objects.end(); ++it) {
+
+ // normalizing datasets to the maximum intensity and scale, if Chi squared module have normalizer module defined
+ OutputDataNormalizerScaleAndShift *data_normalizer = dynamic_cast<OutputDataNormalizerScaleAndShift *>((*it)->getChiSquaredModule()->getOutputDataNormalizer());
+ if( data_normalizer) {
+ data_normalizer->setMaximumIntensity(max_intensity);
+ }
+ double chi_squared = (*it)->calculateChiSquared();
+ chi_sum += chi_squared;
+ std::cout << " chi " << chi_squared << " chi_sum:" << chi_sum << std::endl;
+ }
+ return chi_sum;
}
diff --git a/Core/Tools/src/FitSuiteParameters.cpp b/Core/Tools/src/FitSuiteParameters.cpp
index 404e3afd031b6da0404f2964fc959bd9d708464b..ebec4c4fcc9143e095cde25b7105e62c77c67fc1 100644
--- a/Core/Tools/src/FitSuiteParameters.cpp
+++ b/Core/Tools/src/FitSuiteParameters.cpp
@@ -79,6 +79,18 @@ 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++;
+ }
+ std::cout << "XXX FitSuiteParameters::getNfreeParameters() " << 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)