diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetDetOffsetsMultiPeaks.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetDetOffsetsMultiPeaks.h index 2e392a0fd83d..d34349967783 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetDetOffsetsMultiPeaks.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetDetOffsetsMultiPeaks.h @@ -4,6 +4,7 @@ #include "MantidAPI/Algorithm.h" #include "MantidKernel/System.h" #include "MantidDataObjects/OffsetsWorkspace.h" +#include "MantidDataObjects/EventWorkspace.h" #include "MantidDataObjects/TableWorkspace.h" #include "MantidAPI/MatrixWorkspace.h" #include @@ -15,6 +16,24 @@ namespace Mantid { namespace Algorithms { + +struct FitPeakOffsetResult +{ + double mask; + double offset; + double chi2; + /// ??? + double fitSum; + /// summation of chi-square + double chisqSum; + /// + double peakPosFittedSize; + int numpeakstofit; + int numpeaksfitted; + int numpeaksindrange; + std::string fitoffsetstatus; +}; + /** Find the offsets for each detector @@ -65,10 +84,30 @@ class DLLExport GetDetOffsetsMultiPeaks: public API::Algorithm // Overridden Algorithm methods void init(); void exec(); + + /// Generate output information table workspace + Mantid::DataObjects::TableWorkspace_sptr createOutputInfoTable(); + + /// Generate output peak information table workspace + Mantid::DataObjects::TableWorkspace_sptr createOutputPeakOffsetTable(); + + FitPeakOffsetResult calculatePeakOffset(const int wi, std::vector& fittedpeakpositions, std::vector& tofitpeakpositions); + + void makeFitSummary(const DataObjects::TableWorkspace_sptr& infows); + + API::MatrixWorkspace_sptr inputW; + DataObjects::EventWorkspace_const_sptr eventW; + bool isEvent; + std::string m_backType; std::string m_peakType; double m_maxChiSq; double m_minPeakHeight; + + double maxOffset; + + std::vector peakPositions; + std::vector fitWindows; }; } // namespace Algorithm diff --git a/Code/Mantid/Framework/Algorithms/src/GetDetOffsetsMultiPeaks.cpp b/Code/Mantid/Framework/Algorithms/src/GetDetOffsetsMultiPeaks.cpp index 42e78e110615..abf17df0c222 100644 --- a/Code/Mantid/Framework/Algorithms/src/GetDetOffsetsMultiPeaks.cpp +++ b/Code/Mantid/Framework/Algorithms/src/GetDetOffsetsMultiPeaks.cpp @@ -6,6 +6,12 @@ The algorithm iterates over each spectrum in the workspace and fits a [[Gaussian This is then written into a [[CalFile|.cal file]] for every detector that contributes to that spectrum. All of the entries in the cal file are initially set to both be included, but also to all group into a single group on [[DiffractionFocussing]]. The [[CreateCalFileByNames]] algorithm can be used to alter the grouping in the cal file. +== Fit for peak offset == +The algorithm to calculate offset of peaks' positions is to +minimize a cost function as + \sum_{pi} |X_{0, pi} - (1+offset)X_{0, pi}|/chi^2_{pi} +, which pi is the index of a peak whose position is within MinD and MaxD. + == Usage == '''Python''' @@ -44,6 +50,7 @@ namespace Algorithms { /// Factor to convert full width half max to sigma for calculations of I/sigma. const double FWHM_TO_SIGMA = 2.0*sqrt(2.0*std::log(2.0)); + const double BAD_OFFSET(1000.); // mark things that didn't work with this /** * Helper function for calculating costs in gsl. @@ -53,6 +60,7 @@ namespace Algorithms */ double gsl_costFunction(const gsl_vector *v, void *params) { + // FIXME - there is no need to use vectors peakPosToFit, peakPosFitted and chisq double *p = (double *)params; size_t n = static_cast(p[0]); std::vector peakPosToFit(n); @@ -119,7 +127,7 @@ namespace Algorithms { declareProperty(new WorkspaceProperty<>("InputWorkspace","",Direction::Input, - boost::make_shared("dSpacing")),"A 2D workspace with X values of d-spacing"); + boost::make_shared("dSpacing")),"A 2D matrix workspace with X values of d-spacing"); declareProperty(new ArrayProperty("DReference"),"Enter a comma-separated list of the expected X-position of the centre of the peaks. Only peaks near these positions will be fitted." ); declareProperty("FitWindowMaxWidth", 0., @@ -162,8 +170,11 @@ namespace Algorithms auto offsetwsprop = new WorkspaceProperty("PeaksOffsetTableWorkspace", "PeakOffsetTable", Direction::Output); declareProperty(offsetwsprop, "Name of an output table workspace containing peaks' offset data."); + declareProperty("OutputFitSummary", false, "If specified, a summary on fitting peak summary is made as log's notice level. "); + } + //---------------------------------------------------------------------------------------------- /** * The windows should be half of the distance between the peaks of maxWidth, whichever is smaller. * @param dmin The minimum d-spacing for the workspace @@ -215,10 +226,13 @@ namespace Algorithms */ void GetDetOffsetsMultiPeaks::exec() { - const double BAD_OFFSET(1000.); // mark things that didn't work with this - MatrixWorkspace_sptr inputW=getProperty("InputWorkspace"); - double maxOffset=getProperty("MaxOffset"); + //************************************************************************* + // Process input information + //************************************************************************* + // MatrixWorkspace_sptr inputW=getProperty("InputWorkspace"); + inputW=getProperty("InputWorkspace"); + maxOffset=getProperty("MaxOffset"); int nspec=static_cast(inputW->getNumberHistograms()); // Create the output OffsetsWorkspace OffsetsWorkspace_sptr outputW(new OffsetsWorkspace(inputW->getInstrument())); @@ -232,9 +246,11 @@ namespace Algorithms //To get the workspace index from the detector ID const detid2index_map pixel_to_wi = maskWS->getDetectorIDToWorkspaceIndexMap(true); // the peak positions and where to fit - std::vector peakPositions = getProperty("DReference"); + // std::vector peakPositions = getProperty("DReference"); + peakPositions = getProperty("DReference"); std::sort(peakPositions.begin(), peakPositions.end()); - std::vector fitWindows = generateWindows(wkspDmin, wkspDmax, peakPositions, this->getProperty("FitWindowMaxWidth")); + // std::vector fitWindows = generateWindows(wkspDmin, wkspDmax, peakPositions, this->getProperty("FitWindowMaxWidth")); + fitWindows = generateWindows(wkspDmin, wkspDmax, peakPositions, this->getProperty("FitWindowMaxWidth")); g_log.information() << "windows : "; if (fitWindows.empty()) { @@ -248,8 +264,10 @@ namespace Algorithms } // some shortcuts for event workspaces - EventWorkspace_const_sptr eventW = boost::dynamic_pointer_cast( inputW ); - bool isEvent = false; + // EventWorkspace_const_sptr eventW = boost::dynamic_pointer_cast( inputW ); + eventW = boost::dynamic_pointer_cast( inputW ); + // bool isEvent = false; + isEvent = false; if (eventW) isEvent = true; @@ -260,27 +278,20 @@ namespace Algorithms m_maxChiSq = this->getProperty("MaxChiSq"); m_minPeakHeight = this->getProperty("MinimumPeakHeight"); - // Create output information tableworkspace - auto m_infoTableWS = boost::make_shared(); - m_infoTableWS->addColumn("int", "WorkspaceIndex"); - m_infoTableWS->addColumn("int", "NumberPeaksFitted"); - m_infoTableWS->addColumn("int", "NumberPeaksInRange"); - m_infoTableWS->addColumn("str", "OffsetFitStatus"); - m_infoTableWS->addColumn("double", "ChiSquare"); + //************************************************************************* + // Process output workspaces + //************************************************************************* + // Output: + TableWorkspace_sptr m_infoTableWS = createOutputInfoTable(); setProperty("SpectraFitInfoTableWorkspace", m_infoTableWS); // Create output peak offset workspace - auto m_peakOffsetTableWS = boost::make_shared(); - m_peakOffsetTableWS->addColumn("int", "WorkspaceIndex"); - for (size_t i = 0; i < peakPositions.size(); ++i) - { - std::stringstream namess; - namess << "@" << std::setprecision(5) << peakPositions[i]; - m_peakOffsetTableWS->addColumn("str", namess.str()); - } - m_peakOffsetTableWS->addColumn("double", "OffsetDeviation"); + TableWorkspace_sptr m_peakOffsetTableWS = createOutputPeakOffsetTable(); setProperty("PeaksOffsetTableWorkspace", m_peakOffsetTableWS); + //************************************************************************* + // Calculate offset of each detector + //************************************************************************* // Fit all the spectra with a gaussian Progress prog(this, 0, 1.0, nspec); @@ -289,6 +300,8 @@ namespace Algorithms for (int wi=0;wi fittedpeakpositions, tofitpeakpositions; + FitPeakOffsetResult offsetresult = calculatePeakOffset(wi, fittedpeakpositions, tofitpeakpositions); + +#endif // Get the list of detectors in this pixel const std::set & dets = inputW->getSpectrum(wi)->getDetectorIDs(); @@ -443,36 +461,44 @@ namespace Algorithms std::set::iterator it; for (it = dets.begin(); it != dets.end(); ++it) { - outputW->setValue(*it, offset, fitSum); - outputNP->setValue(*it, peakPosFittedSize, chisqSum); + // Set value to output peak offset workspace + outputW->setValue(*it, offsetresult.offset, offsetresult.fitSum); + // Set value to output peak number workspace + outputNP->setValue(*it, offsetresult.peakPosFittedSize, offsetresult.chisqSum); + // Set value to mask workspace const auto mapEntry = pixel_to_wi.find(*it); - if ( mapEntry == pixel_to_wi.end() ) continue; + if ( mapEntry == pixel_to_wi.end() ) + continue; const size_t workspaceIndex = mapEntry->second; - if (mask == 1.) + if (offsetresult.mask == 1.) { // Being masked maskWS->maskWorkspaceIndex(workspaceIndex); - maskWS->dataY(workspaceIndex)[0] = mask; + maskWS->dataY(workspaceIndex)[0] = offsetresult.mask; } else { // Using the detector - maskWS->dataY(workspaceIndex)[0] = mask; + maskWS->dataY(workspaceIndex)[0] = offsetresult.mask; } - } + } // ENDFOR (detectors) // Add report to TableWorkspace TableRow newrow = m_infoTableWS->appendRow(); - newrow << wi << numpeaksfitted << numpeaksindrange << fitoffsetstatus << chi2; + newrow << wi << offsetresult.numpeaksfitted << offsetresult.numpeaksindrange + << offsetresult.fitoffsetstatus << offsetresult.chi2; // Add report to offset info table + // FIXME - This does not make much sense... + // TODO - Best way is to re-fit the spectra with given offset and see how good + // the peak positions are matching the theoretical value TableRow newrow2 = m_peakOffsetTableWS->appendRow(); newrow2 << wi; - if (numpeaksfitted > 0) + if (offsetresult.numpeaksfitted > 0) { - std::vector haspeakvec(numpeakstofit, false); - std::vector deltavec(numpeakstofit, 0.0); - for (int i = 0; i < numpeaksfitted; ++i) + std::vector haspeakvec(offsetresult.numpeakstofit, false); + std::vector deltavec(offsetresult.numpeakstofit, 0.0); + for (int i = 0; i < offsetresult.numpeaksfitted; ++i) { double peakcentre = tofitpeakpositions[i]; int index = static_cast( @@ -489,7 +515,7 @@ namespace Algorithms double sumdelta1 = 0.0; double sumdelta2 = 0.0; double numdelta = 0.0; - for (int i = 0; i < numpeakstofit; ++i) + for (int i = 0; i < offsetresult.numpeakstofit; ++i) { if (haspeakvec[i]) { @@ -511,15 +537,11 @@ namespace Algorithms double stddev = sumdelta2/numdelta - (sumdelta1/numdelta)*(sumdelta1/numdelta); newrow2 << stddev; - } - else - { - // Do nothing - } + } // ENDIF (numpeaksfitted > 0) + } // End of critical region - } prog.report(); PARALLEL_END_INTERUPT_REGION } @@ -542,6 +564,12 @@ namespace Algorithms childAlg->executeAsChildAlg(); } + // Make summary + bool dosummary = getProperty("OutputFitSummary"); + if (dosummary) + makeFitSummary(m_infoTableWS); + + return; } //----------------------------------------------------------------------------------------- @@ -591,8 +619,8 @@ namespace Algorithms * @param nparams :: Number of parameters. * @param minD :: Min distance. * @param maxD :: Max distance. - * @param peakPosToFit :: Peak positions to fit. - * @param peakPosFitted :: Peak positions fitted. + * @param peakPosToFit :: Peak positions to fit (output). + * @param peakPosFitted :: Peak positions fitted (output). * @param chisq :: chisq. * @return The number of peaks in range */ @@ -809,5 +837,237 @@ namespace Algorithms return numPeaksInRange; } + //---------------------------------------------------------------------------------------------- + /** Create information table workspace for output + */ + TableWorkspace_sptr GetDetOffsetsMultiPeaks::createOutputInfoTable() + { + // Create output information tableworkspace + auto m_infoTableWS = boost::make_shared(); + m_infoTableWS->addColumn("int", "WorkspaceIndex"); + m_infoTableWS->addColumn("int", "NumberPeaksFitted"); + m_infoTableWS->addColumn("int", "NumberPeaksInRange"); + m_infoTableWS->addColumn("str", "OffsetFitStatus"); + m_infoTableWS->addColumn("double", "ChiSquare"); + + return m_infoTableWS; + } + + //---------------------------------------------------------------------------------------------- + /** Create peak offset table workspace for output + */ + TableWorkspace_sptr GetDetOffsetsMultiPeaks::createOutputPeakOffsetTable() + { + auto m_peakOffsetTableWS = boost::make_shared(); + m_peakOffsetTableWS->addColumn("int", "WorkspaceIndex"); + for (size_t i = 0; i < peakPositions.size(); ++i) + { + std::stringstream namess; + namess << "@" << std::setprecision(5) << peakPositions[i]; + m_peakOffsetTableWS->addColumn("str", namess.str()); + } + m_peakOffsetTableWS->addColumn("double", "OffsetDeviation"); + + return m_peakOffsetTableWS; + } + + //---------------------------------------------------------------------------------------------- + /** Calculate offset for one spectrum + */ + FitPeakOffsetResult GetDetOffsetsMultiPeaks::calculatePeakOffset(const int wi, std::vector& fittedpeakpositions, std::vector& tofitpeakpositions) + { + // Initialize the structure to return + FitPeakOffsetResult fr; + + fr.offset = 0.0; + fr.fitoffsetstatus = "N/A"; + fr.chi2 = -1; + + fr.fitSum = 0.0; + fr.chisqSum = 0.0; + + // Number of peak fitted of this spectrum + fr.peakPosFittedSize = 0.0; + + fr.numpeaksfitted = 0; + fr.numpeakstofit = 0; + fr.numpeaksindrange = 0; + + // checks for dead detectors + if ((isEvent) && (eventW->getEventList(wi).empty())) + { + // empty detector will be masked + fr.offset = BAD_OFFSET; + fr.fitoffsetstatus = "empty det"; + } + else + { + // dead detector will be masked + const MantidVec& Y = inputW->readY(wi); + const int YLength = static_cast(Y.size()); + double sumY = 0.0; + for (int i = 0; i < YLength; i++) sumY += Y[i]; + if (sumY < 1.e-30) + { + // Dead detector will be masked + fr.offset=BAD_OFFSET; + fr.fitoffsetstatus = "dead det"; + } + } + + if (fr.offset < 10.) + { + // Fit the peak + std::vector peakPosToFit, peakPosFitted, chisq; + size_t nparams; + double minD, maxD; + fr.numpeaksindrange = + fitSpectra(wi, inputW, peakPositions, fitWindows, nparams, minD, maxD, + peakPosToFit, peakPosFitted, chisq); + fr.numpeakstofit = static_cast(peakPositions.size()); + fr.numpeaksfitted = static_cast(peakPosFitted.size()); + fittedpeakpositions = peakPosFitted; + tofitpeakpositions = peakPosToFit; + + // Fit offset + if (nparams > 0 && fr.numpeaksindrange > 0) + { + //double * params = new double[2*nparams+1]; + double params[153]; + if(nparams > 50) nparams = 50; + params[0] = static_cast(nparams); + params[1] = minD; + params[2] = maxD; + for (size_t i = 0; i < nparams; i++) + { + params[i+3] = peakPosToFit[i]; + } + for (size_t i = 0; i < nparams; i++) + { + params[i+3+nparams] = peakPosFitted[i]; + } + fr.peakPosFittedSize = static_cast(peakPosFitted.size()); + for (size_t i = 0; i < nparams; i++) + { + params[i+3+2*nparams] = chisq[i]; + fr.chisqSum += chisq[i]; + } + + const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; + gsl_multimin_fminimizer *s = NULL; + gsl_vector *ss, *x; + gsl_multimin_function minex_func; + + // finally do the fitting + size_t nopt = 1; + size_t iter = 0; + int status = 0; + double size; + + /* Starting point */ + x = gsl_vector_alloc (nopt); + gsl_vector_set_all (x, 0.0); + + /* Set initial step sizes to 0.001 */ + ss = gsl_vector_alloc (nopt); + gsl_vector_set_all (ss, 0.001); + + /* Initialize method and iterate */ + minex_func.n = nopt; + minex_func.f = &gsl_costFunction; + minex_func.params = ¶ms; + + s = gsl_multimin_fminimizer_alloc (T, nopt); + gsl_multimin_fminimizer_set (s, &minex_func, x, ss); + + do + { + iter++; + status = gsl_multimin_fminimizer_iterate(s); + if (status) + break; + + size = gsl_multimin_fminimizer_size (s); + status = gsl_multimin_test_size (size, 1e-4); + + } + while (status == GSL_CONTINUE && iter < 50); + + // Output summary to log file + std::string reportOfDiffractionEventCalibrateDetectors = gsl_strerror(status); + g_log.debug() << " Workspace Index = " << wi << + " Method used = " << " Simplex" << + " Iteration = " << iter << + " Status = " << reportOfDiffractionEventCalibrateDetectors << + " Minimize Sum = " << s->fval << + " Offset = " << gsl_vector_get (s->x, 0) << " \n"; + fr.offset = gsl_vector_get (s->x, 0); + fr.fitSum = s->fval; + gsl_vector_free(x); + gsl_vector_free(ss); + gsl_multimin_fminimizer_free (s); + + fr.fitoffsetstatus = reportOfDiffractionEventCalibrateDetectors; + fr.chi2 = s->fval; + // FIXME - Whether chi2 is normalized by number of peaks entering minimization + //delete [] params; + } + else + { + // Not enough peaks have been found. + // Output warning + g_log.debug() << "Spectra " << wi << " has 0 parameter for it. Set to bad_offset." << ".\n"; + fr.offset = BAD_OFFSET; + fr.fitoffsetstatus = "no peaks"; + } + } + + // Final check offset + fr.mask=0.0; + if (std::abs(fr.offset) > maxOffset) + { + fr.offset = 0.0; + fr.mask = 1.0; + } + + return fr; + } /// ENDFUNCTION: GetDetOffsetsMultiPeaks + + //---------------------------------------------------------------------------------------------- + /** Make a summary of fitting + */ + void GetDetOffsetsMultiPeaks::makeFitSummary(const TableWorkspace_sptr& infows) + { + size_t numrows = infows->rowCount(); + double sumchi2 = 0; + double weight_sumchi2 = 0; + size_t weight_numfittedpeaks = 0; + size_t numunmasked = 0; + for (size_t i = 0; i < numrows; ++i) + { + TableRow row = infows->getRow(i); + int wi, numpeakfitted, numpeakinrange; + double chi2, offset; + std::string fitofsetstatus; + row >> wi >> numpeakfitted >> numpeakinrange >> fitofsetstatus >> chi2 >> offset; + if (numpeakfitted * numpeakinrange > 0) + { + ++ numunmasked; + sumchi2 += chi2; + weight_sumchi2 += static_cast(numpeakfitted*chi2); + weight_numfittedpeaks += numpeakfitted; + } + } + + double avgchi2 = sumchi2/static_cast(numunmasked); + double wtavgchi2 = weight_sumchi2/static_cast(weight_numfittedpeaks); + + g_log.notice() << "Unmasked spectra = " << numunmasked << " \t" + << "Average chi-sq = " << avgchi2 << " \t" + << "Weighted avg. chi-sq = " << wtavgchi2 << "\n"; + + } + + } // namespace Algorithm } // namespace Mantid