diff --git a/lib/XYSampleData.cpp b/lib/XYSampleData.cpp index de60a46..5ec3071 100644 --- a/lib/XYSampleData.cpp +++ b/lib/XYSampleData.cpp @@ -250,7 +250,8 @@ const XYStatData & XYSampleData::getData(void) } // fit ///////////////////////////////////////////////////////////////////////// -SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, +SampleFitResult XYSampleData::fit(std::vector &minimizer, + const DVec &init, const std::vector &v) { computeVarMat(); @@ -265,12 +266,17 @@ SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, FOR_STAT_ARRAY(result, s) { setDataToSample(s); - sampleResult = data_.fit(minimizer, initCopy, v); + if (s == central) + { + sampleResult = data_.fit(minimizer, initCopy, v); + } + else + { + sampleResult = data_.fit(*(minimizer.back()), initCopy, v); + } initCopy = sampleResult.segment(0, initCopy.size()); result[s] = sampleResult; result.chi2_[s] = sampleResult.getChi2(); - - for (unsigned int j = 0; j < v.size(); ++j) { result.model_[j].resize(nSample_); @@ -284,6 +290,15 @@ SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, return result; } +SampleFitResult XYSampleData::fit(Minimizer &minimizer, + const DVec &init, + const std::vector &v) +{ + vector mv{&minimizer}; + + return fit(mv, init, v); +} + // schedule data initilization from samples //////////////////////////////////// void XYSampleData::scheduleDataInit(void) { diff --git a/lib/XYSampleData.hpp b/lib/XYSampleData.hpp index 0e5342f..17f9842 100644 --- a/lib/XYSampleData.hpp +++ b/lib/XYSampleData.hpp @@ -92,11 +92,16 @@ public: // get internal XYStatData const XYStatData & getData(void); // fit + SampleFitResult fit(std::vector &minimizer, const DVec &init, + const std::vector &v); SampleFitResult fit(Minimizer &minimizer, const DVec &init, const std::vector &v); - template + template + SampleFitResult fit(std::vector &minimizer, const DVec &init, + const DoubleModel &model, const Ts... models); + template SampleFitResult fit(Minimizer &minimizer, const DVec &init, - const DoubleModel &model, const Mods... models); + const DoubleModel &model, const Ts... models); private: // schedule data initilization from samples void scheduleDataInit(void); @@ -118,17 +123,30 @@ private: * XYSampleData template implementation * ******************************************************************************/ template -SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, +SampleFitResult XYSampleData::fit(std::vector &minimizer, + const DVec &init, const DoubleModel &model, const Ts... models) { static_assert(static_or::value...>::value, - "model arguments are not compatible with DoubleModel &"); + "model arguments are not compatible with DoubleModel"); std::vector modelVector{&model, &models...}; return fit(minimizer, init, modelVector); } +template +SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, + const DoubleModel &model, const Ts... models) +{ + static_assert(static_or::value...>::value, + "model arguments are not compatible with DoubleModel"); + + std::vector mv{&minimizer}; + + return fit(mv, init, model, models...); +} + END_LATAN_NAMESPACE #endif // Latan_XYSampleData_hpp_ diff --git a/lib/XYStatData.cpp b/lib/XYStatData.cpp index 27878af..3793d53 100644 --- a/lib/XYStatData.cpp +++ b/lib/XYStatData.cpp @@ -24,6 +24,8 @@ using namespace std; using namespace Latan; +static constexpr double maxXsiDev = 10.; + /****************************************************************************** * FitResult implementation * ******************************************************************************/ @@ -220,7 +222,7 @@ const DMat & XYStatData::getFitVarMatPInv(void) } // fit ///////////////////////////////////////////////////////////////////////// -FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, +FitResult XYStatData::fit(vector &minimizer, const DVec &init, const vector &v) { // check model consistency @@ -257,14 +259,31 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, } // minimization - FitResult result; - DVec totalInit(totalNPar); + FitResult result; + DVec totalInit(totalNPar); + //// set total init vector totalInit.segment(0, nPar) = init; totalInit.segment(nPar, layout.totalXSize) = chi2DataVec_.segment(layout.totalYSize, layout.totalXSize); - minimizer.setInit(totalInit); - result = minimizer(chi2); + for (auto &m: minimizer) + { + m->setInit(totalInit); + //// do not allow more than maxXsiDev std. deviations on the x-axis + for (Index p = nPar; p < totalNPar; ++p) + { + double err; + + err = sqrt(fitVar_.diagonal()(layout.totalYSize + p - nPar)); + m->useLowLimit(p); + m->useHighLimit(p); + m->setLowLimit(p, totalInit(p) - maxXsiDev*err); + m->setHighLimit(p, totalInit(p) + maxXsiDev*err); + } + //// minimize and store results + result = (*m)(chi2); + totalInit = result; + } result.chi2_ = chi2(result); result.nPar_ = nPar; result.nDof_ = layout.totalYSize - nPar; @@ -281,6 +300,14 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, return result; } +FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, + const vector &v) +{ + vector mv{&minimizer}; + + return fit(mv, init, v); +} + // create data ///////////////////////////////////////////////////////////////// void XYStatData::createXData(const std::string name __dumb, const Index nData) { diff --git a/lib/XYStatData.hpp b/lib/XYStatData.hpp index 8db977f..226ab25 100644 --- a/lib/XYStatData.hpp +++ b/lib/XYStatData.hpp @@ -87,9 +87,14 @@ public: const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); // fit + FitResult fit(std::vector &minimizer, const DVec &init, + const std::vector &v); FitResult fit(Minimizer &minimizer, const DVec &init, const std::vector &v); template + FitResult fit(std::vector &minimizer, const DVec &init, + const DoubleModel &model, const Ts... models); + template FitResult fit(Minimizer &minimizer, const DVec &init, const DoubleModel &model, const Ts... models); protected: @@ -125,7 +130,7 @@ private: * XYStatData template implementation * ******************************************************************************/ template -FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, +FitResult XYStatData::fit(std::vector &minimizer, const DVec &init, const DoubleModel &model, const Ts... models) { static_assert(static_or::value...>::value, @@ -136,6 +141,18 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, return fit(minimizer, init, modelVector); } +template +FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, + const DoubleModel &model, const Ts... models) +{ + static_assert(static_or::value...>::value, + "model arguments are not compatible with DoubleModel"); + + std::vector mv{&minimizer}; + + return fit(mv, init, model, models...); +} + /****************************************************************************** * error check macros * ******************************************************************************/