Mantid
Loading...
Searching...
No Matches
IntegratePeakTimeSlices.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7/*
8 * IntegratePeakTimeSlices.cpp
9 *
10 * Created on: May 5, 2011
11 * Author: ruth
12 *
13 *
14 *
15 */
22#include "MantidAPI/IFunction.h"
25
26#include "MantidHistogramData/BinEdges.h"
27
28#include <boost/math/special_functions/round.hpp>
29#include <utility>
30
31using namespace Mantid::Kernel;
32using namespace Mantid::API;
33using namespace Mantid::DataObjects;
34using namespace Mantid::Geometry;
35using namespace Mantid::HistogramData;
36using namespace std;
37namespace Mantid::Crystal {
38
39DECLARE_ALGORITHM(IntegratePeakTimeSlices)
40
41// Attr, m_AttributeValues, and StatBase indicies
42#define IStartRow 0
43#define IStartCol 1
44#define INRows 2
45#define INCol 3
46#define ISSIxx 4
47#define ISSIyy 5
48#define ISSIxy 6
49#define ISSxx 7
50#define ISSyy 8
51#define ISSxy 9
52#define ISSIx 10
53#define ISSIy 11
54#define ISSx 12
55#define ISSy 13
56#define IIntensities 14
57#define ISS1 15
58#define IVariance 16
59#define ITotBoundary 17
60#define INBoundary 18
61#define IVarBoundary 19
62#define NAttributes 20
63
64// Parameter indicies
65#define IBACK 0
66#define ITINTENS 1
67#define IXMEAN 2
68#define IYMEAN 3
69#define IVXX 4
70#define IVYY 5
71#define IVXY 6
72// TODO: Edge Peaks-Return NoPeak(Intensity=0,variance=0) if any center on any
73// slice is out of bounds
74// TODO: Calc ratio for edge peaks, should use the slant of bivariate normal
75// instead of assuming
76// distribution axes are lined up with the row and col vectors
77#define NParameters 7
78
79namespace {
80// # std sigs 0 .25 .5 .75 1 1.25 1.5 2 2.5
81const double probs[9] = {.5f, .5987f, .6915f, .7734f, .8413f, .8944f, .9322f, .9599f, .9772f};
82
83const int MinRowColSpan = 6;
84const int MaxRowColSpan = 36;
85const int MinTimeSpan = 3;
86const double NeighborhoodRadiusDivPeakRadius = 1.5;
87const double MaxNeighborhoodRadius = 10;
88const double NStdDevPeakSpan = 2;
89const double MaxGoodRatioFitvsExpIntenisites = 2.5;
90const double MinGoodRatioFitvsExpIntenisites = .25;
91const double MinGoodIoverSigI = 3.0;
92const double MinVariationInXYvalues = .6; // Peak spans one pixel only
93const double MaxCorrCoeffinXY = .9; // otherwise all data on one line
94} // namespace
95
97 : Algorithm(), m_R0(-1), m_ROW(0.), m_COL(0.), m_cellWidth(0.), m_cellHeight(0.), m_NROWS(0), m_NCOLS(0) {
98 m_EdgePeak = false;
99 m_NeighborIDs = new int[3];
100 m_NeighborIDs[0] = 3;
101 m_NeighborIDs[1] = 2;
102 m_AttributeNames[0] = "StartRow";
103 m_AttributeNames[1] = "StartCol";
104 m_AttributeNames[2] = "NRows";
105 m_AttributeNames[3] = "NCols";
106 m_AttributeNames[4] = "SSIxx";
107 m_AttributeNames[5] = "SSIyy";
108 m_AttributeNames[6] = "SSIxy";
109 m_AttributeNames[7] = "SSxx";
110 m_AttributeNames[8] = "SSyy";
111 m_AttributeNames[9] = "SSxy";
112 m_AttributeNames[10] = "SSIx";
113 m_AttributeNames[11] = "SSIy";
114 m_AttributeNames[12] = "SSx";
115 m_AttributeNames[13] = "SSy";
116 m_AttributeNames[14] = "Intensities";
117 m_AttributeNames[15] = " SS1";
118 m_AttributeNames[16] = "Variance";
119 m_AttributeNames[17] = "TotBoundary";
120 m_AttributeNames[18] = "NBoundary";
121 m_AttributeNames[19] = "VarianceBoundary";
122
123 m_ParameterNames[0] = "Background";
124 m_ParameterNames[1] = "Intensity";
125 m_ParameterNames[2] = "Mcol";
126 m_ParameterNames[3] = "Mrow";
127 m_ParameterNames[4] = "SScol";
128 m_ParameterNames[5] = "SSrow";
129 m_ParameterNames[6] = "SSrc";
130
131 std::fill(m_ParameterValues.begin(), m_ParameterValues.end(), 0.0);
132}
133
134double SQRT(double v) {
135 if (v < 0)
136 return -1;
137 return sqrt(v);
138}
141
143 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
144 "A 2D workspace with X values of time of flight");
145
146 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("OutputWorkspace", "", Direction::Output),
147 "Name of the output table workspace with Log info");
148
150 "Workspace of Peaks");
151
152 declareProperty("PeakIndex", 0, "Index of peak in PeaksWorkspace to integrate");
153
154 declareProperty("PeakQspan", .06, "Max magnitude of Q of Peak to Q of Peak Center, where mod(Q)=1/d");
155
156 declareProperty("CalculateVariances", true, "Calc (co)variances given parameter values versus fit (co)Variances ");
157
158 declareProperty("Ties", "", "Tie parameters(Background,Intensity, Mrow,...) to values/formulas.");
159
160 declareProperty("NBadEdgePixels", 0, "Number of bad Edge Pixels");
161
162 declareProperty("Intensity", 0.0, "Peak Integrated Intensity", Direction::Output);
163
164 declareProperty("SigmaIntensity", 0.0, "Peak Integrated Intensity Error", Direction::Output);
165}
166
177 time_t seconds1;
178
179 seconds1 = time(nullptr);
180
181 double dQ = getProperty("PeakQspan");
182
183 g_log.debug("------------------Start Peak Integrate-------------------");
184
185 if (dQ <= 0) {
186 g_log.error("Negative PeakQspans are not allowed. Use .17/G where G is the "
187 "max unit cell length");
188 throw std::runtime_error("Negative PeakQspans are not allowed in IntegratePeakTimeSlices");
189 }
190
191 MatrixWorkspace_const_sptr inpWkSpace = getProperty("InputWorkspace");
192 if (!inpWkSpace) {
193 g_log.error("Improper Input Workspace");
194 throw std::runtime_error("Improper Input Workspace in IntegratePeakTimeSlices");
195 }
196
197 PeaksWorkspace_sptr peaksW;
198 peaksW = getProperty("Peaks");
199 if (!peaksW) {
200 g_log.error("Improper Peaks Input");
201 throw std::runtime_error("Improper Peaks Input");
202 }
203
204 int indx = getProperty("PeakIndex");
205
206 Peak &peak = peaksW->getPeak(indx);
207
208 //------------------------------- Get Panel
209 //--------------------------------------
210 std::shared_ptr<const Geometry::IComponent> panel_const =
211 peak.getInstrument()->getComponentByName(peak.getBankName());
212
213 std::shared_ptr<Geometry::IComponent> panel = std::const_pointer_cast<Geometry::IComponent>(panel_const);
214
215 if (!panel || !panel_const) {
216 g_log.error("Cannot get panel for a peak");
217 throw std::runtime_error("Cannot get panel for a peak");
218 }
219
220 BoundingBox box;
221 panel->getBoundingBox(box);
222
223 int detID = peak.getDetectorID();
224 if (!box.isPointInside(peak.getDetPos())) {
225 g_log.error("Detector pixel is NOT inside the Peaks Bank");
226 throw std::runtime_error("Detector pixel is NOT inside the Peaks Bank");
227 }
228
230
231 g_log.debug() << " Peak Index " << indx << '\n';
232
233 double TotVariance = 0;
234 double TotIntensity = 0;
235 double lastRow = m_ROW;
236 double Row0 = lastRow;
237 double lastCol = m_COL;
238 double Col0 = lastCol;
239 string spec_idList;
240
241 // For quickly looking up workspace index from det id
242 m_wi_to_detid_map = inpWkSpace->getDetectorIDToWorkspaceIndexMap();
243
244 TableWorkspace_sptr TabWS = std::make_shared<TableWorkspace>(0);
245
246 //----------------------------- get Peak extents
247 //------------------------------
248 try {
249
250 // Find the workspace index for this detector ID
251 detid2index_map::const_iterator it = m_wi_to_detid_map.find(detID);
252 size_t wsIndx = (it->second);
253
254 double R = CalculatePositionSpan(peak, dQ) / 2;
255
256 R = min<double>(MaxRowColSpan * max<double>(m_cellWidth, m_cellHeight), R);
257 R = max<double>(MinRowColSpan * max<double>(m_cellWidth, m_cellHeight), R);
258
259 R = 2 * R; // Gets a few more background cells.
260 int Chan;
261
262 const auto &X = inpWkSpace->x(wsIndx);
263 int dChan = CalculateTimeChannelSpan(peak, dQ, X, int(wsIndx), Chan);
264
265 dChan = max<int>(dChan, MinTimeSpan);
266
267 double Centy = Row0;
268 double Centx = Col0;
269 IDetector_const_sptr CenterDet = peak.getDetector();
270
271 double neighborRadius; // last radius for finding neighbors
272
273 neighborRadius = min<double>(MaxNeighborhoodRadius, NeighborhoodRadiusDivPeakRadius * R);
274 auto Nneighbors = static_cast<int>(neighborRadius * neighborRadius / m_cellWidth / m_cellHeight * 4);
275
276 Nneighbors = min<int>(Nneighbors, static_cast<int>(inpWkSpace->getNumberHistograms()) - 2);
277 delete[] m_NeighborIDs;
278
279 m_NeighborIDs = new int[Nneighbors + 2];
280 m_NeighborIDs[0] = Nneighbors + 2;
281 m_NeighborIDs[1] = 2;
282 Kernel::V3D Cent = (m_center + m_xvec * (Centx - m_COL) + m_yvec * (Centy - m_ROW));
283
284 getNeighborPixIDs(panel, Cent, neighborRadius, m_NeighborIDs);
285
286 if (m_NeighborIDs[1] < 10) {
287 g_log.error("Not enough neighboring pixels to fit ");
288 throw std::runtime_error("Not enough neighboring pixels to fit ");
289 }
290 int NBadEdgeCells = getProperty("NBadEdgePixels");
291 int MaxChan = -1;
292 double MaxCounts = -1;
293
294 // --------------- Find Time Chan with max counts----------------
295 for (int dir = 1; dir > -2; dir -= 2) {
296 bool done = false;
297 for (int t = 0; t < dChan && !done; t++)
298 if (dir < 0 && t == 0) {
299 Centy = Row0;
300 Centx = Col0;
301 } else if (Chan + dir * t < 0 || Chan + dir * t >= static_cast<int>(X.size()))
302 done = true;
303 else {
304
305 int NN = m_NeighborIDs[1];
306 MatrixWorkspace_sptr Data = WorkspaceFactory::Instance().create(std::string("Workspace2D"), 3, NN, NN);
307
308 auto XXX = std::make_shared<DataModeHandler>(R, R, Centy, Centx, m_cellWidth, m_cellHeight,
309 getProperty("CalculateVariances"), NBadEdgeCells,
310 m_NCOLS - NBadEdgeCells, NBadEdgeCells, m_NROWS - NBadEdgeCells);
311 m_AttributeValues = XXX;
312 XXX->setCurrentRadius(R);
313
314 SetUpData1(Data, inpWkSpace, Chan + dir * t, Chan + dir * t, R, CenterDet->getPos(), spec_idList);
315
316 if (m_AttributeValues->StatBaseVals(ISSIxx) > 0) {
317 if (m_AttributeValues->StatBaseVals(IIntensities) > MaxCounts) {
318 MaxCounts = m_AttributeValues->StatBaseVals(IIntensities);
319 MaxChan = Chan + dir * t;
320 }
321 if (m_AttributeValues->StatBaseVals(IIntensities) > 0) {
322 Centx = m_AttributeValues->StatBaseVals(ISSIx) / m_AttributeValues->StatBaseVals(IIntensities);
323 Centy = m_AttributeValues->StatBaseVals(ISSIy) / m_AttributeValues->StatBaseVals(IIntensities);
324 } else
325 done = true;
326 } else
327 done = true;
328
329 if (t >= 3 && (m_AttributeValues->StatBaseVals(IIntensities) < MaxCounts / 2.0) && MaxCounts >= 0)
330 done = true;
331 }
332 }
333 if (MaxChan > 0)
334 Chan = MaxChan;
335
336 g_log.debug() << " largest Channel,Radius,m_cellWidth,m_cellHeight = " << Chan << " " << R << " " << m_cellWidth
337 << " " << m_cellHeight << '\n';
338
339 if (R < MinRowColSpan / 2 * max<double>(m_cellWidth, m_cellHeight) || dChan < MinTimeSpan) {
340 g_log.error("Not enough rows and cols or time channels ");
341 throw std::runtime_error("Not enough rows and cols or time channels ");
342 }
343
345
346 //------------------------------------- Start the Integrating
347 //-------------------------------
348 double time;
349 int ncells;
350
351 Mantid::API::Progress prog(this, 0.0, 1.0, dChan);
352
353 // Set from attributes replace by m_R0
354 m_R0 = -1;
355 int LastTableRow = -1;
356 auto origAttributeList = std::make_shared<DataModeHandler>();
357 auto lastAttributeList = std::make_shared<DataModeHandler>();
358
359 for (int dir = 1; dir >= -1; dir -= 2) {
360 bool done = false;
361
362 for (int chan = 0; chan < dChan && !done; chan++)
363 if (dir < 0 && chan == 0) {
364 lastRow = Row0;
365 lastCol = Col0;
366 lastAttributeList = origAttributeList;
367 if (TabWS->rowCount() > 0)
368 LastTableRow = 0;
369
370 } else if (Chan + dir * chan < 0 || Chan + dir * chan >= static_cast<int>(X.size()))
371 done = true;
372 else {
373
374 int xchan = Chan + dir * chan;
375
376 size_t topIndex = xchan + 1;
377 if (topIndex >= X.size())
378 topIndex = X.size() - 1;
379
380 time = (X[xchan] + X[topIndex]) / 2.0;
381
382 double Radius = R;
383
384 if (m_R0 > 0)
385 Radius = m_R0;
386
388 WorkspaceFactory::Instance().create(std::string("Workspace2D"), 3, m_NeighborIDs[1], m_NeighborIDs[1]);
389
390 g_log.debug() << " A:chan=" << xchan << " time=" << time << " Radius=" << Radius << "row= " << lastRow
391 << " col=" << lastCol << '\n';
392
393 SetUpData(Data, inpWkSpace, panel, xchan, xchan, lastCol, lastRow, Cent, neighborRadius, Radius, spec_idList);
394
395 m_AttributeValues->setTime(time);
396
397 // if( dir==1 && chan ==0)
398 // origAttributeList= m_AttributeValues;
399
400 ncells = static_cast<int>(m_AttributeValues->StatBaseVals(ISS1));
401
402 std::vector<double> params;
403 std::vector<double> errs;
404 std::vector<std::string> names;
405
406 if (m_AttributeValues->StatBaseVals(ISSIxx) > 0 &&
407 m_AttributeValues->IsEnoughData(m_ParameterValues.data(), g_log) && m_ParameterValues[ITINTENS] > 0) {
408 double chisqOverDOF;
409
410 Fit(Data, chisqOverDOF, done, names, params, errs, lastRow, lastCol, neighborRadius);
411
412 if (!done) // Bivariate error happened
413 {
414
415 if (isGoodFit(params, errs, names, chisqOverDOF)) {
416 LastTableRow = UpdateOutputWS(TabWS, dir, xchan, params, errs, names, chisqOverDOF,
417 m_AttributeValues->time, spec_idList);
418
419 double TotSliceIntensity = m_AttributeValues->StatBaseVals(IIntensities);
420 double TotSliceVariance = m_AttributeValues->StatBaseVals(IVariance);
421
422 updatePeakInformation(params, errs, names, TotVariance, TotIntensity, TotSliceIntensity,
423 TotSliceVariance, chisqOverDOF, ncells);
424
425 lastAttributeList = m_AttributeValues;
426
427 if (dir == 1 && chan == 0)
428 origAttributeList = lastAttributeList;
429 } else
430
431 done = true;
432 }
433
434 } else //(!IsEnoughData() || ParameterValues[ITINTENS] <= 0
435 {
436 done = true;
437 }
438
439 if (done) // try to merge
440 {
441 done = false;
442
443 int chanMin, chanMax;
444 if ((dir == 1 && chan == 0) || lastAttributeList->CellHeight <= 0) {
445 chanMin = xchan;
446 chanMax = xchan + 1;
447 if (dir < 0)
448 chanMax++;
449 auto XXX = std::make_shared<DataModeHandler>(*m_AttributeValues);
450 m_AttributeValues = XXX;
451 if (!X.empty())
452 m_AttributeValues->setTime((X[chanMax] + X[chanMin]) / 2.0);
453
454 } else // lastAttributeList exists
455
456 {
457 chanMin = std::min<int>(xchan, xchan - dir);
458 chanMax = chanMin + 1;
459 if (lastAttributeList->case4)
460 chanMax++;
461
462 auto XXX = std::make_shared<DataModeHandler>(*lastAttributeList);
463 m_AttributeValues = XXX;
464
465 m_AttributeValues->setTime((time + m_AttributeValues->time) / 2.0);
466 }
467
468 if (updateNeighbors(panel, m_AttributeValues->getCurrentCenter(), Cent,
469 m_AttributeValues->getCurrentRadius(), neighborRadius))
470 Cent = m_AttributeValues->getCurrentCenter();
471
472 Data =
473 WorkspaceFactory::Instance().create(std::string("Workspace2D"), 3, m_NeighborIDs[1], m_NeighborIDs[1]);
474
475 SetUpData1(Data, inpWkSpace, chanMin, chanMax, m_AttributeValues->getCurrentRadius(),
476 m_AttributeValues->getCurrentCenter(), spec_idList);
477
478 double chisqOverDOF;
479
480 g_log.debug("Try Merge 2 time slices");
481 if (m_AttributeValues->StatBaseVals(ISSIxx) >= 0 &&
482 m_AttributeValues->IsEnoughData(m_ParameterValues.data(), g_log))
483
484 Fit(Data, chisqOverDOF, done, names, params, errs, lastRow, lastCol, neighborRadius);
485 else
486 chisqOverDOF = -1;
487
488 if (!done && isGoodFit(params, errs, names, chisqOverDOF)) {
489
490 if (LastTableRow >= 0 && LastTableRow < static_cast<int>(TabWS->rowCount()))
491 TabWS->removeRow(LastTableRow);
492 else
493 LastTableRow = -1;
494
495 LastTableRow = UpdateOutputWS(TabWS, dir, (chanMin + chanMax) / 2.0, params, errs, names, chisqOverDOF,
496 m_AttributeValues->time, spec_idList);
497
498 if (lastAttributeList->lastISAWVariance > 0 && lastAttributeList->CellHeight > 0) {
499 TotIntensity -= lastAttributeList->lastISAWIntensity;
500 TotVariance -= lastAttributeList->lastISAWVariance;
501 }
502
503 double TotSliceIntensity = m_AttributeValues->StatBaseVals(IIntensities);
504
505 double TotSliceVariance = m_AttributeValues->StatBaseVals(IVariance);
506
507 updatePeakInformation(params, errs, names, TotVariance, TotIntensity, TotSliceIntensity, TotSliceVariance,
508 chisqOverDOF, static_cast<int>(m_AttributeValues->StatBaseVals(ISS1)));
509
510 // lastAttributeList= m_AttributeValues;
511
512 if (dir == 1 && (chan == 0 || chan == 1)) {
513 m_AttributeValues->case4 = true;
514 origAttributeList = m_AttributeValues;
515 } else
516 LastTableRow = -1;
517
518 } else {
519 auto XXX = std::make_shared<DataModeHandler>();
520 lastAttributeList = XXX;
521 }
522 done = true;
523 }
524
525 // Get ready for the next round
526 Data.reset();
527
528 if (!done) {
529
530 // Now set up the center for this peak
531 int i = findNameInVector("Mrow", names);
532 if (i < 0) {
533 throw std::runtime_error("Inconsistency found in algorithm "
534 "execution. The index for the parameter "
535 "Mrow is negative.");
536 }
537
538 lastRow = boost::math::iround(params[i]);
539 i = findNameInVector("Mcol", names);
540 if (i >= 0)
541 lastCol = boost::math::iround(params[i]);
542 prog.report();
543
544 } else if (dir > 0)
545 prog.report(dChan / 2);
546 else
547 prog.report(dChan);
548
549 params.clear();
550 errs.clear();
551 names.clear();
552 }
553 }
554
555 } catch (std::exception &EE1) {
556 std::cout << "Error in main reason=" << EE1.what() << '\n';
557
558 throw std::runtime_error(" Error IntegratePeakTimeSlices:" + std::string(EE1.what()));
559 } catch (std::string &mess) {
560 throw std::runtime_error("Error IntegratePeakTimeSlices:" + mess);
561
562 } catch (...) {
563 throw std::runtime_error("Error IntegratePeakTimeSlices:");
564 }
565
566 try {
567
568 setProperty("OutputWorkspace", TabWS);
569
570 setProperty("Intensity", TotIntensity);
571 setProperty("SigmaIntensity", SQRT(TotVariance));
572 time_t seconds2;
573
574 seconds2 = time(nullptr);
575 double dif = difftime(seconds2, seconds1);
576 g_log.debug() << "Finished Integr peak number " << indx << " in " << dif << " seconds\n";
577
578 } catch (std::exception &ss) {
579
580 std::cout << "Error occurred XX " << ss.what() << '\n';
581 throw std::runtime_error(ss.what());
582 }
583}
584
594bool IntegratePeakTimeSlices::getNeighborPixIDs(const std::shared_ptr<Geometry::IComponent> &comp,
595 const Kernel::V3D &Center, double &Radius, int *&ArryofID) {
596
597 int N = ArryofID[1];
598 int MaxN = ArryofID[0];
599
600 if (N >= MaxN)
601 return false;
602
604 comp->getBoundingBox(box);
605
606 double minx = Center.X() - Radius;
607 double miny = Center.Y() - Radius;
608 double minz = Center.Z() - Radius;
609 double maxx = Center.X() + Radius;
610 double maxy = Center.Y() + Radius;
611 double maxz = Center.Z() + Radius;
612
613 if (box.xMin() >= maxx)
614 return true;
615
616 if (box.xMax() <= minx)
617 return true;
618 ;
619
620 if (box.yMin() >= maxy)
621 return true;
622
623 if (box.yMax() <= miny)
624 return true;
625
626 if (box.zMin() >= maxz)
627 return true;
628
629 if (box.zMax() <= minz)
630 return true;
631 ;
632
633 auto det = std::dynamic_pointer_cast<Geometry::Detector>(comp);
634
635 if (det) {
636 V3D pos = det->getPos() - Center;
637 if (pos.X() * pos.X() + pos.Y() * pos.Y() + pos.Z() * pos.Z() < Radius * Radius) {
638 ArryofID[N] = det->getID();
639 N++;
640 ArryofID[1] = N;
641 }
642 return true;
643 }
644
645 auto Assembly = std::dynamic_pointer_cast<const Geometry::ICompAssembly>(comp);
646
647 if (!Assembly)
648 return true;
649
650 bool b = true;
651
652 for (int i = 0; i < Assembly->nelements() && b; i++)
653 b = getNeighborPixIDs(Assembly->getChild(i), Center, Radius, ArryofID);
654
655 return b;
656}
657
666bool IntegratePeakTimeSlices::updateNeighbors(const std::shared_ptr<Geometry::IComponent> &comp, const V3D &CentPos,
667 const V3D &oldCenter, double NewRadius, double &neighborRadius) {
668 double DD = (CentPos - oldCenter).norm();
669 bool changed = false;
670 if (DD + NewRadius > neighborRadius) {
671 auto NN = int(NStdDevPeakSpan * NeighborhoodRadiusDivPeakRadius * NewRadius / m_cellWidth * NStdDevPeakSpan *
672 NeighborhoodRadiusDivPeakRadius * NewRadius / m_cellHeight);
673 if (m_NeighborIDs[0] < NN) {
674 delete[] m_NeighborIDs;
675 m_NeighborIDs = new int[NN + 2];
676 m_NeighborIDs[0] = NN + 2;
677 }
678 m_NeighborIDs[1] = 2;
679 neighborRadius = NeighborhoodRadiusDivPeakRadius * NewRadius;
680
681 getNeighborPixIDs(comp, CentPos, neighborRadius, m_NeighborIDs);
682 changed = true;
683
684 } else // big enough neighborhood so
685 neighborRadius -= DD;
686
687 return changed;
688}
689
700double IntegratePeakTimeSlices::CalculatePositionSpan(Peak const &peak, const double dQ) {
701
702 try {
703 double Q = 0, ScatAngle = 0, dScatAngle = 0, DetSpan = 0;
704
705 Q = peak.getQLabFrame().norm();
707 const Geometry::IComponent_const_sptr sample = instr->getSample();
708 V3D pos = peak.getDetPos() - sample->getPos();
709
710 ScatAngle = acos(pos.Z() / pos.norm());
711
712 dScatAngle = 2 * dQ / Q * tan(ScatAngle / 2);
713
714 DetSpan = pos.norm() * dScatAngle; // s=r*theta
715
716 DetSpan = fabs(DetSpan);
717
718 // IDetector_sptr det = peak.getDetector();
719
720 return DetSpan;
721
722 } catch (std::exception &s) {
723 std::cout << "err in getNRowsCols, reason=" << s.what() << '\n';
724 return 0;
725 }
726}
727
742int IntegratePeakTimeSlices::CalculateTimeChannelSpan(Geometry::IPeak const &peak, const double dQ, const HistogramX &X,
743 const int specNum, int &Centerchan) {
744 UNUSED_ARG(specNum);
745 double Q = peak.getQLabFrame().norm(); // getQ( peak)/2/M_PI;
746
747 double time = peak.getTOF();
748 double dtime = dQ / Q * time;
749 int chanCenter = findTimeChannel(X, time);
750
751 Centerchan = chanCenter;
752 int chanLeft = findTimeChannel(X, time - dtime);
753 int chanRight = findTimeChannel(X, time + dtime);
754 int dchan = abs(chanCenter - chanLeft);
755
756 if (abs(chanRight - chanCenter) > dchan)
757 dchan = abs(chanRight - chanCenter);
758
759 dchan = max<int>(3, dchan);
760
761 return dchan + 5; // heuristic should be a lot more
762}
763
780void IntegratePeakTimeSlices::FindPlane(V3D &center, V3D &xvec, V3D &yvec, double &ROW, double &COL, int &NROWS,
781 int &NCOLS, double &pixWidthx, double &pixHeighty,
782 DataObjects::Peak const &peak) const {
783
784 NROWS = NCOLS = -1;
786 V3D detPos = det->getPos();
787
788 center.setX(detPos.X());
789 center.setY(detPos.Y());
790 center.setZ(detPos.Z());
791
792 std::shared_ptr<const Detector> dett = std::dynamic_pointer_cast<const Detector>(det);
793
794 pixWidthx = dett->getWidth();
795 pixHeighty = dett->getHeight();
796
797 Kernel::Quat Qt = dett->getRotation();
798 V3D yaxis0(0.0, 1.0, 0.0);
799
800 Qt.rotate(yaxis0);
801 yaxis0.normalize();
802
803 V3D xaxis0(1, 0, 0);
804 Qt.rotate(xaxis0);
805 xaxis0.normalize();
806
807 xvec.setX(xaxis0.X());
808 xvec.setY(xaxis0.Y());
809 xvec.setZ(xaxis0.Z());
810 yvec.setX(yaxis0.X());
811 yvec.setY(yaxis0.Y());
812 yvec.setZ(yaxis0.Z());
813 ROW = peak.getRow();
814 COL = peak.getCol();
816 if (!inst)
817 throw std::invalid_argument("No instrument for peak");
818 std::shared_ptr<const IComponent> panel = inst->getComponentByName(peak.getBankName());
819
820 std::shared_ptr<const RectangularDetector> ddet = std::dynamic_pointer_cast<const RectangularDetector>(panel);
821
822 if (ddet) {
823 std::pair<int, int> CR = ddet->getXYForDetectorID(det->getID());
824 ROW = CR.second;
825 COL = CR.first;
826 pixWidthx = ddet->xstep();
827 pixHeighty = ddet->ystep();
828
829 NROWS = ddet->ypixels();
830 NCOLS = ddet->xpixels();
831
832 return;
833 }
834 // Get NROWS and NCOLS for other panels
835 NROWS = NCOLS = -1;
836
837 if (!inst)
838 return;
839
840 if (!panel)
841 return;
842 std::shared_ptr<const Component> compPanel = std::dynamic_pointer_cast<const Component>(panel);
843 std::shared_ptr<IComponent> panel1(compPanel->base()->clone());
844 BoundingBox B;
845
846 Quat rot = panel1->getRotation();
847
848 rot.inverse();
849
850 panel1->rotate(rot);
851
852 panel1->getBoundingBox(B);
853
854 NROWS = boost::math::iround((B.yMax() - B.yMin()) / pixHeighty);
855 NCOLS = boost::math::iround((B.xMax() - B.xMin()) / pixWidthx);
856}
857
867void IntegratePeakTimeSlices::updateStats(const double intensity, const double variance, const double row,
868 const double col, std::vector<double> &StatBase) {
869
870 StatBase[ISSIxx] += col * col * intensity;
871 StatBase[ISSIyy] += intensity * row * row;
872 StatBase[ISSIxy] += intensity * row * col;
873 StatBase[ISSxx] += col * col;
874 StatBase[ISSyy] += row * row;
875 StatBase[ISSxy] += row * col;
876 StatBase[ISSIx] += intensity * col;
877 StatBase[ISSIy] += intensity * row;
878 StatBase[ISSx] += col;
879 StatBase[ISSy] += row;
880 StatBase[IIntensities] += intensity;
881 StatBase[IVariance] += variance;
882 StatBase[ISS1] += 1;
883}
884
895std::vector<double> DataModeHandler::InitValues(double Varx, double Vary, double b) {
896 std::vector<double> Res(7);
897
898 Res[IVXX] = Varx;
899 Res[IVYY] = Vary;
900 Res[IVXY] = 0;
901 auto nCells = static_cast<int>(StatBase[ISS1]);
902 double Den = StatBase[IIntensities] - b * nCells;
903 Res[IXMEAN] = (StatBase[ISSIx] - b * StatBase[ISSx]) / Den;
904 Res[IYMEAN] = (StatBase[ISSIy] - b * StatBase[ISSy]) / Den;
905 Res[IBACK] = b;
906 Res[ITINTENS] = StatBase[IIntensities] - b * nCells;
907
908 //---- Is Edge Cell ???-------
909 double NstdX = 4 * (currentRadius / CellWidth - EdgeX) / sqrt(Varx);
910 double NstdY = 4 * (currentRadius / CellHeight - EdgeY) / sqrt(Vary);
911 double sigx = 1;
912 double sigy = 1;
913 if (NstdX < 0)
914 sigx = -1;
915 if (NstdY < 0)
916 sigy = -1;
917
918 double x = 1;
919 if (sigy * NstdY < 7 && sigy * NstdY >= 0) // is close to row edge
920 {
921 x = probs[std::lround(sigy * NstdY)];
922 if (sigy < 0)
923 x = 1 - x;
924 double My2 = StatBase[IStartRow];
925 if (Res[IYMEAN] - My2 > My2 + StatBase[INRows] - Res[IYMEAN])
926 My2 += StatBase[INRows];
927 Res[IYMEAN] = Res[IYMEAN] * x + (1 - x) * My2;
928 }
929 double x1 = 1;
930 if (sigx * NstdX < 7 && sigx * NstdX > 0) // is close to x edge
931 {
932 x1 = probs[std::lround(sigx * NstdX)];
933 if (sigx < 0)
934 x1 = 1 - x1;
935 double Mx2 = StatBase[IStartCol];
936 if (Res[IXMEAN] - Mx2 > Mx2 + StatBase[INCol] - Res[IXMEAN])
937 Mx2 += StatBase[INCol];
938 Res[IXMEAN] = Res[IXMEAN] * x1 + (1 - x1) * Mx2;
939 }
940 Res[ITINTENS] /= x * x1;
941
942 return Res;
943}
944
952std::vector<double> DataModeHandler::GetParams(double b) {
953
954 auto nCells = static_cast<int>(StatBase[ISS1]);
955 double Den = StatBase[IIntensities] - b * nCells;
956 double Varx, Vary;
957
958 Varx = VarxHW;
959 Vary = VaryHW;
960
961 double Rx = lastRCRadius / CellWidth - EdgeX;
962 double Ry = lastRCRadius / CellHeight - EdgeY;
963 if (Varx <= 0)
965
966 if (Vary <= 0)
968 // Use
969 if (Rx * Rx < 4 * Varx || Ry * Ry < 4 * Vary) {
970 return InitValues(Varx, Vary, b);
971 }
972 if (Den < 0)
973 return std::vector<double>();
974
975 double Mx = StatBase[ISSIx] - b * StatBase[ISSx];
976 double My = StatBase[ISSIy] - b * StatBase[ISSy];
977
978 double Sxx = (StatBase[ISSIxx] - b * StatBase[ISSxx] - Mx * Mx / Den) / Den;
979 double Syy = (StatBase[ISSIyy] - b * StatBase[ISSyy] - My * My / Den) / Den;
980 double Sxy = (StatBase[ISSIxy] - b * StatBase[ISSxy] - Mx * My / Den) / Den;
981
982 double Intensity = StatBase[IIntensities] - b * nCells;
983 double col = Mx / Den;
984 double row = My / Den;
985 std::vector<double> Result(7);
986 Result[IBACK] = b;
987 Result[ITINTENS] = Intensity;
988 Result[IXMEAN] = col;
989 Result[IYMEAN] = row;
990 Result[IVXX] = Sxx;
991 Result[IVYY] = Syy;
992 Result[IVXY] = Sxy;
993
994 return Result;
995}
996
1003bool DataModeHandler::setStatBase(std::vector<double> const &statBase)
1004
1005{
1006 double TotBoundaryIntensities = statBase[ITotBoundary];
1007 auto nBoundaryCells = static_cast<int>(statBase[INBoundary]);
1008 this->StatBase = statBase;
1009 double b = 0;
1010 if (nBoundaryCells > 0)
1011 b = TotBoundaryIntensities / nBoundaryCells;
1012
1013 auto nCells = static_cast<int>(statBase[ISS1]);
1014 double Den = statBase[IIntensities] - b * nCells;
1015 int k = 0;
1016 while (Den <= 0 && b != 0) {
1017 b = b * .7;
1018 Den = statBase[IIntensities] - b * nCells;
1019
1020 if (k < 8)
1021 k++;
1022 else
1023 b = 0;
1024 }
1025
1026 double Varx, Vary;
1027 Varx = statBase[INCol] / 7; // Range = 3.5 standard deviations
1028 Vary = statBase[INRows] / 7;
1029 Varx *= Varx;
1030 Vary *= Vary;
1031
1032 double Rx = lastRCRadius / CellWidth - EdgeX;
1033 double Ry = lastRCRadius / CellHeight - EdgeY;
1034 if (CellWidth > 0 && currentRadius > 0 && lastCol > 0 && lastRow > 0)
1035 if (Rx * Rx < 4 * std::max(Varx, VarxHW) || HalfWidthAtHalfHeightRadius < 0 ||
1036 Ry * Ry < 4 * std::max(Vary, VaryHW)) // Edge peak so cannot use samples
1037 {
1038 Vx_calc = VarxHW;
1039 Vy_calc = VaryHW;
1040 Vxy_calc = 0;
1041 col_calc = lastCol;
1042 row_calc = lastRow;
1043 back_calc = b;
1044 Intensity_calc = statBase[IIntensities] - b * nCells;
1045 if (Vx_calc <= 0 || Vy_calc <= 0) // EdgePeak but not big enuf
1046 return true;
1047
1049 double r = CalcSampleIntensityMultiplier(params);
1050 Intensity_calc *= r;
1051 return true;
1052 }
1053 if (Den <= 0)
1054 Den = 1;
1055
1056 bool done = false;
1057 int ntimes = 0;
1058 double Mx, My, Sxx, Syy, Sxy;
1059
1060 double RangeX = statBase[INCol] / 2;
1061 double RangeY = statBase[INRows] / 2;
1062
1063 while (!done && ntimes < 29) {
1064 Mx = statBase[ISSIx] - b * statBase[ISSx];
1065 My = statBase[ISSIy] - b * statBase[ISSy];
1066 Sxx = (statBase[ISSIxx] - b * statBase[ISSxx] - Mx * Mx / Den) / Den;
1067 Syy = (statBase[ISSIyy] - b * statBase[ISSyy] - My * My / Den) / Den;
1068 Sxy = (statBase[ISSIxy] - b * statBase[ISSxy] - Mx * My / Den) / Den;
1069 ntimes++;
1070 done = false;
1071
1072 if (Sxx <= RangeX / 12 || Syy <= RangeY / 12 || Sxy * Sxy / Sxx / Syy > .9) {
1073 b = b * .95;
1074
1075 if (ntimes + 1 == 29)
1076 b = 0;
1077
1078 Den = statBase[IIntensities] - b * nCells;
1079 if (Den <= 1)
1080 Den = 1;
1081
1082 } else
1083 done = true;
1084 }
1085
1086 back_calc = b;
1087 Intensity_calc = statBase[IIntensities] - b * nCells;
1088 col_calc = Mx / Den;
1089 row_calc = My / Den;
1090 Vx_calc = Sxx;
1091 Vy_calc = Syy;
1092 Vxy_calc = Sxy;
1093 return false;
1094}
1095
1103 double Vx, Vy;
1104 Vx = VarxHW;
1105 Vy = VaryHW;
1106 if (Vx < 0)
1108 if (Vy < 0)
1110
1111 double Rx = lastRCRadius / CellWidth - EdgeX;
1112 double Ry = lastRCRadius / CellHeight - EdgeY;
1113 double mult = 1;
1114 if (Rx * Rx > 4 * Vx)
1115 Vx = std::max(VarxHW, Vx_calc);
1116 else
1117 mult = 1.35;
1118
1119 if (Ry * Ry > 4 * Vy)
1120 Vy = std::max(VaryHW, Vy_calc);
1121 else
1122 mult *= 1.35;
1123
1124 double DD = max<double>(sqrt(Vy) * CellHeight, sqrt(Vx) * CellWidth);
1125 double NewRadius = 1.4 * max<double>(MinRowColSpan * max<double>(CellWidth, CellHeight), 4.5 * DD);
1126 NewRadius = mult * min<double>(baseRCRadius, NewRadius);
1127 // 1.4 is needed to get more background cells. In rectangle the corners were
1128 // background
1129
1130 NewRadius = min<double>(MaxRowColSpan * max<double>(CellWidth, CellHeight), NewRadius);
1131
1132 return NewRadius;
1133}
1134
1144void DataModeHandler::setHeightHalfWidthInfo(const std::vector<double> &xvals, const std::vector<double> &yvals,
1145 const std::vector<double> &counts) {
1146 double minCount, maxCount;
1147 const auto &X = xvals;
1148 const auto &Y = yvals;
1149 const auto &C = counts;
1150 VarxHW = -1;
1151 VaryHW = -1;
1152 auto N = static_cast<int>(X.size());
1153
1155
1156 if (N <= 2)
1157 return;
1158
1159 minCount = maxCount = C[0];
1160 double MaxX = -1;
1161 double MaxY = -1;
1162 int nmax = 0;
1163 double lowX, lowY, highX, highY;
1164 lowX = highX = X[0];
1165 lowY = highY = Y[0];
1166
1167 for (int i = 1; i < N; i++) {
1168 if (X[i] < lowX)
1169 lowX = X[i];
1170 else if (X[i] > highX)
1171 highX = X[i];
1172
1173 if (Y[i] < lowY)
1174 lowY = Y[i];
1175 else if (Y[i] > highY)
1176 highY = Y[i];
1177
1178 if (C[i] > maxCount) {
1179 maxCount = C[i];
1180 MaxX = X[i];
1181 MaxY = Y[i];
1182 nmax = 1;
1183 } else if (C[i] < minCount) {
1184 minCount = C[i];
1185
1186 } else if (C[i] == maxCount) // Get a tolerance on this
1187 {
1188 MaxX += X[i];
1189 MaxY += Y[i];
1190 nmax++;
1191 }
1192 }
1193 if (minCount == maxCount)
1194 return;
1195
1196 MaxX /= nmax;
1197 MaxY /= nmax;
1198
1199 double dCount = std::max(.51, (maxCount - minCount) / 6.2);
1200 double CountUp = (maxCount + minCount) / 2 + dCount;
1201 double CountLow = (maxCount + minCount) / 2 - dCount;
1202 double dSpanx = (highX - lowX) / 6.;
1203 double dSpany = (highY - lowY) / 6.0;
1204
1205 int nMax = 0;
1206 int nMin = 0;
1207 double TotMax = 0;
1208 double TotMin = 0;
1209 double offset = std::max(.2, (maxCount - minCount) / 20);
1210 double TotR_max = 0;
1211 double TotR_min = 0;
1212 double TotRx0 = 0;
1213 double TotRy0 = 0;
1214 double TotCx = 0;
1215 double TotCy = 0;
1216 for (int i = 0; i < N; i++) {
1217 if (C[i] > maxCount - offset) {
1218 TotMax += C[i];
1219 nMax++;
1220 TotR_max += C[i] * sqrt((X[i] - MaxX) * (X[i] - MaxX) + (Y[i] - MaxY) * (Y[i] - MaxY));
1221 }
1222 if (C[i] < minCount + offset)
1223
1224 {
1225 TotMin += C[i];
1226 nMin++;
1227
1228 TotR_min += C[i] * sqrt((X[i] - MaxX) * (X[i] - MaxX) + (Y[i] - MaxY) * (Y[i] - MaxY));
1229 }
1230
1231 if (fabs(MaxY - Y[i]) < 1.2 && fabs(MaxX - X[i]) > 1.2 && C[i] >= CountLow && C[i] <= CountUp &&
1232 fabs(MaxX - X[i]) < dSpanx) {
1233 TotRx0 += (C[i] - minCount) * (X[i] - MaxX) * (X[i] - MaxX);
1234 TotCx += C[i] - minCount;
1235 }
1236
1237 if (fabs(MaxX - X[i]) < 1.2 && fabs(MaxY - Y[i]) > 1.2 && C[i] >= CountLow && C[i] <= CountUp &&
1238 fabs(MaxY - Y[i]) < dSpany) {
1239 TotRy0 += (C[i] - minCount) * (Y[i] - MaxY) * (Y[i] - MaxY);
1240 TotCy += C[i] - minCount;
1241 }
1242 }
1243
1244 if (nMax + nMin == N) // all data are on two levels essentially
1245 {
1246 if (TotMax <= 0)
1247 TotMax = 1;
1248 if (TotMin <= 0)
1249 TotMin = 1;
1250 double AvR = .5 * (TotR_max / TotMax + TotR_min / TotMin);
1251 HalfWidthAtHalfHeightRadius = AvR / .8326;
1252
1255 return;
1256 }
1257
1258 double TotR = 0, nR = -1, nRx = -1, nRy = -1;
1259 double MidVal = (TotMax / nMax + TotMin / nMin) / 2.0;
1260 double TotRx = 0, TotRy = 0;
1261 while ((nR <= 0 || nRy <= 0 || nRx <= 0) && offset < MidVal) {
1262 TotR = 0;
1263 nR = 0;
1264 TotRx = 0;
1265 TotRy = 0;
1266 nRx = 0;
1267 nRy = 0;
1268
1269 for (int i = 0; i < N; i++)
1270 if (C[i] < MidVal + offset && C[i] > MidVal - offset) {
1271 double X1 = X[i] - MaxX;
1272 double Y1 = Y[i] - MaxY;
1273 TotR += sqrt(X1 * X1 + Y1 * Y1);
1274 nR++;
1275 if ((X1 >= -1.2 && X1 <= 1.2) && fabs(Y1) > 1.2 && fabs(Y1) < dSpany) {
1276 nRy++;
1277 TotRy += abs(Y1);
1278 }
1279 if ((Y1 >= -1.2 && Y1 <= 1.2) && fabs(X1) > 1.2 && fabs(X1) < dSpanx) {
1280 nRx++;
1281 TotRx += fabs(X1);
1282 }
1283 }
1284 offset *= 1.1;
1285 }
1286
1287 double AvR = TotR / nR;
1288 HalfWidthAtHalfHeightRadius = AvR / .8326;
1289
1290 if (nRx > 0)
1291 VarxHW = (TotRx / nRx) * (TotRx / nRx) / .8326 / .8326;
1292 else if (TotCx > 0)
1293 VarxHW = TotRx0 * TotRx0 / TotCx / TotCx / .8326 / .8326;
1294 else if (HalfWidthAtHalfHeightRadius > 0)
1296 else
1297 VarxHW = -1;
1298
1299 if (nRy > 0)
1300 VaryHW = (TotRy / nRy) * (TotRy / nRy) / .8326 / .8326;
1301 else if (TotCy > 0)
1302 VaryHW = TotRy0 * TotRy0 / TotCy / TotCy / .8326 / .8326;
1303 else if (HalfWidthAtHalfHeightRadius > 0)
1305 else
1306 VaryHW = -1;
1307}
1328 const std::shared_ptr<Geometry::IComponent> &comp, const int chanMin,
1329 const int chanMax, double CentX, double CentY, Kernel::V3D &CentNghbr,
1330 double &neighborRadius, // from CentDetspec
1331 double Radius, string &spec_idList) {
1332
1333 Kernel::V3D CentPos1 = m_center + m_xvec * (CentX - m_COL) * m_cellWidth + m_yvec * (CentY - m_ROW) * m_cellHeight;
1334
1335 int NBadEdgeCells = getProperty("NBadEdgePixels");
1336
1337 auto X = std::make_shared<DataModeHandler>(Radius, Radius, CentY, CentX, m_cellWidth, m_cellHeight,
1338 getProperty("CalculateVariances"), NBadEdgeCells, m_NCOLS - NBadEdgeCells,
1339 NBadEdgeCells, m_NROWS - NBadEdgeCells);
1340
1342 m_AttributeValues->setCurrentRadius(Radius);
1343 m_AttributeValues->setCurrentCenter(CentPos1);
1344
1345 SetUpData1(Data, inpWkSpace, chanMin, chanMax, Radius, CentPos1, spec_idList);
1346
1347 if (m_AttributeValues->StatBaseVals(ISSIxx) < 0) // Not enough data
1348 return;
1349
1350 double NewRadius = m_AttributeValues->getNewRCRadius();
1351 if (m_R0 > 0) {
1352 NewRadius = m_R0;
1353 } else {
1354 m_R0 = NewRadius;
1355 }
1356
1357 CentX = m_ParameterValues[IXMEAN];
1358 CentY = m_ParameterValues[IYMEAN];
1359 Kernel::V3D CentPos = m_center + m_xvec * (CentX - m_COL) * m_cellWidth + m_yvec * (CentY - m_ROW) * m_cellHeight;
1360
1361 double DD = (CentPos - CentNghbr).norm();
1362
1363 if (DD + NewRadius > neighborRadius) {
1364 auto NN = int(NStdDevPeakSpan * NeighborhoodRadiusDivPeakRadius * NewRadius / m_cellWidth * NStdDevPeakSpan *
1365 NeighborhoodRadiusDivPeakRadius * NewRadius / m_cellHeight);
1366 if (m_NeighborIDs[0] < NN) {
1367 delete[] m_NeighborIDs;
1368 m_NeighborIDs = new int[NN + 2];
1369 m_NeighborIDs[0] = NN + 2;
1370 } // else
1371 // NN= m_NeighborIDs[0]-2;
1372 m_NeighborIDs[1] = 2;
1373 neighborRadius = NeighborhoodRadiusDivPeakRadius * NewRadius;
1374 CentNghbr = CentPos;
1375 getNeighborPixIDs(comp, CentPos, neighborRadius, m_NeighborIDs);
1376
1377 } else // big enough neighborhood so
1378 neighborRadius -= DD;
1379
1380 // if( changed) CentNghbr = CentPos.
1381 auto X1 = std::make_shared<DataModeHandler>(Radius, NewRadius, CentY, CentX, m_cellWidth, m_cellHeight,
1382 getProperty("CalculateVariances"), NBadEdgeCells, m_NCOLS - NBadEdgeCells,
1383 NBadEdgeCells, m_NROWS - NBadEdgeCells);
1384
1385 m_AttributeValues = X1;
1386 m_AttributeValues->setCurrentRadius(NewRadius);
1387 m_AttributeValues->setCurrentCenter(CentPos);
1388 SetUpData1(Data, inpWkSpace, chanMin, chanMax, NewRadius, CentPos, spec_idList);
1389}
1390
1403 API::MatrixWorkspace_const_sptr const &inpWkSpace, const int chanMin,
1404 const int chanMax, double Radius, const Kernel::V3D &CentPos,
1405 string &spec_idList) {
1407 if (m_NeighborIDs[1] < 10) {
1408 return;
1409 }
1410 std::vector<double> StatBase(NAttributes);
1411 std::shared_ptr<Workspace2D> ws = std::dynamic_pointer_cast<Workspace2D>(Data);
1412
1413 int NBadEdges = getProperty("NBadEdgePixels");
1414 spec_idList.clear();
1415
1416 for (int i = 0; i < NAttributes + 2; i++)
1417 StatBase.emplace_back(0);
1418
1419 std::vector<double> yvalB;
1420 std::vector<double> errB;
1421 std::vector<double> xvalB;
1422 std::vector<double> YvalB;
1423
1424 double TotBoundaryIntensities = 0;
1425 int nBoundaryCells = 0;
1426 double TotBoundaryVariances = 0;
1427
1428 double BoundaryRadius = min<double>(.90 * Radius, Radius - 1.5 * max<double>(m_cellWidth, m_cellHeight));
1429 double minRow = 20000, maxRow = -1, minCol = 20000, maxCol = -1;
1430
1431 int jj = 0;
1432
1433 std::vector<double> xRef;
1434 for (int i = 2; i < m_NeighborIDs[1]; i++) {
1435 int DetID = m_NeighborIDs[i];
1436
1437 size_t workspaceIndex;
1438 if (m_wi_to_detid_map.count(DetID) > 0)
1439 workspaceIndex = m_wi_to_detid_map.find(DetID)->second;
1440 else {
1441 throw std::runtime_error("No workspaceIndex for detID=" + std::to_string(DetID));
1442 }
1443
1444 IDetector_const_sptr Det = inpWkSpace->getDetector(workspaceIndex);
1445 V3D pixPos = Det->getPos();
1446
1447 if (i > 2)
1448 spec_idList += ",";
1449
1450 V3D dist = pixPos - CentPos;
1451 if (dist.scalar_prod(dist) < Radius * Radius)
1452
1453 {
1454 spec_idList += std::to_string(inpWkSpace->getSpectrum(workspaceIndex).getSpectrumNo());
1455
1456 double R1 = dist.scalar_prod(m_yvec);
1457 double R1a = R1 / m_cellHeight;
1458
1459 double row = m_ROW + R1a;
1460
1461 double C1 = dist.scalar_prod(m_xvec);
1462 double C1a = C1 / m_cellWidth;
1463
1464 double col = m_COL + C1a;
1465
1466 if (row > NBadEdges && col > NBadEdges && (m_NROWS < 0 || row < m_NROWS - NBadEdges) &&
1467 (m_NCOLS < 0 || col < m_NCOLS - NBadEdges)) {
1468 auto &histogram = inpWkSpace->y(workspaceIndex);
1469
1470 auto &histoerrs = inpWkSpace->e(workspaceIndex);
1471 double intensity = 0;
1472 double variance = 0;
1473 for (int chan = chanMin; chan <= chanMax; chan++) {
1474 intensity += histogram[chan];
1475 variance += histoerrs[chan] * histoerrs[chan];
1476 }
1477
1478 yvalB.emplace_back(intensity);
1479 double sigma = 1;
1480
1481 errB.emplace_back(sigma);
1482 xvalB.emplace_back(col);
1483 YvalB.emplace_back(row);
1484
1485 xRef.emplace_back(static_cast<double>(jj));
1486 jj++;
1487
1488 updateStats(intensity, variance, row, col, StatBase);
1489
1490 if ((pixPos - CentPos).norm() > BoundaryRadius) {
1491 TotBoundaryIntensities += intensity;
1492 nBoundaryCells++;
1493
1494 TotBoundaryVariances += variance;
1495 }
1496
1497 if (row < minRow)
1498 minRow = row;
1499 if (col < minCol)
1500 minCol = col;
1501 if (row > maxRow)
1502 maxRow = row;
1503 if (col > maxCol)
1504 maxCol = col;
1505
1506 } // if not bad edge
1507
1508 } // peak within radius
1509 } // for each neighbor
1510
1511 m_AttributeValues->EdgeY =
1512 max<double>(0.0, max<double>(-m_ROW + minRow + Radius / m_cellHeight, -maxRow + m_ROW + Radius / m_cellHeight));
1513 m_AttributeValues->EdgeX =
1514 max<double>(0.0, max<double>(-m_COL + minCol + Radius / m_cellWidth, -maxCol + m_COL + Radius / m_cellWidth));
1515 if (m_AttributeValues->EdgeY <= 1)
1516 m_AttributeValues->EdgeY = 0;
1517 if (m_AttributeValues->EdgeX <= 1)
1518 m_AttributeValues->EdgeX = 0;
1519
1520 auto pX = Kernel::make_cow<HistogramData::HistogramX>(std::move(xRef));
1521 Data->setX(0, pX);
1522 Data->setX(1, pX);
1523 Data->setX(2, pX);
1524
1525 ws->setCounts(0, yvalB);
1526 ws->setCountStandardDeviations(0, errB);
1527 ws->setCounts(1, xvalB);
1528 ws->setCounts(2, YvalB);
1529 m_AttributeValues->setHeightHalfWidthInfo(xvalB, YvalB, yvalB);
1530
1531 StatBase[IStartRow] = minRow;
1532 StatBase[IStartCol] = minCol;
1533 StatBase[INRows] = maxRow - minRow + 1;
1534 StatBase[INCol] = maxCol - minCol + 1;
1535
1536 StatBase[ITotBoundary] = TotBoundaryIntensities;
1537 StatBase[INBoundary] = nBoundaryCells;
1538 StatBase[IVarBoundary] = TotBoundaryVariances;
1539 m_EdgePeak = m_AttributeValues->setStatBase(StatBase);
1540
1541 m_ParameterValues[IBACK] = m_AttributeValues->getInitBackground();
1542 m_ParameterValues[ITINTENS] = m_AttributeValues->getInitIntensity();
1545 m_ParameterValues[IVXX] = m_AttributeValues->getInitVarx();
1546 m_ParameterValues[IVYY] = m_AttributeValues->getInitVary();
1547 m_ParameterValues[IVXY] = m_AttributeValues->getInitVarxy();
1548}
1549
1556int IntegratePeakTimeSlices::findTimeChannel(const HistogramX &X, const double time) {
1557 int sgn = 1;
1558
1559 if (X[0] > X[1])
1560 sgn = -1;
1561
1562 if (sgn * (X[0] - time) >= 0)
1563 return 0;
1564
1565 if (sgn * (time - X[X.size() - 1u]) >= 0)
1566 return static_cast<int>(X.size()) - 1;
1567
1568 size_t end = X.size() - 1u;
1569 for (size_t i = 0; i < end; i++) {
1570 if (sgn * (time - X[i]) >= 0 && sgn * (X[i + 1u] - time) >= 0)
1571 return static_cast<int>(i);
1572 }
1573
1574 return -1;
1575}
1576
1584bool DataModeHandler::isEdgePeak(const double *params, int nparams) {
1586 double Vary = Varx;
1587 if (nparams > 4) {
1588 Varx = params[IVXX];
1589 Vary = params[IVYY];
1590 }
1591
1592 if (Varx <= 0 || Vary <= 0 || HalfWidthAtHalfHeightRadius <= 0)
1593 return true;
1594
1595 double Rx = lastRCRadius / CellWidth - EdgeX; // span from center in x direction
1596 double Ry = lastRCRadius / CellHeight - EdgeY; // span from center in y direction
1597
1598 return Rx * Rx < NStdDevPeakSpan * NStdDevPeakSpan * std::max(Varx, VarxHW) ||
1599 Ry * Ry < NStdDevPeakSpan * NStdDevPeakSpan * std::max(Vary, VaryHW);
1600}
1601
1608
1609 std::ostringstream fun_str;
1610
1611 fun_str << "name=BivariateNormal,";
1612
1613 if (m_AttributeValues->CalcVariances())
1614 fun_str << "CalcVariances=1";
1615 else
1616 fun_str << "CalcVariances=-1";
1617
1618 int NN = NParameters;
1619 if (m_AttributeValues->CalcVariances())
1620 NN -= 3;
1621
1622 for (int i = 0; i < NN; i++) {
1623 fun_str << "," << m_ParameterNames[i] << "=" << m_ParameterValues[i];
1624 }
1625
1626 return fun_str.str();
1627}
1628
1637int IntegratePeakTimeSlices::findNameInVector(std::string const &oneName, std::vector<std::string> const &nameList)
1638
1639{
1640 for (size_t i = 0; i < nameList.size(); i++)
1641 if (oneName == nameList[i])
1642 return static_cast<int>(i);
1643
1644 return -1;
1645}
1646
1653
1654 this->baseRCRadius = handler.baseRCRadius;
1655 this->lastRCRadius = handler.lastRCRadius;
1657 this->calcNewRCRadius = handler.calcNewRCRadius;
1658 this->lastRow = handler.lastRow;
1659 this->lastCol = handler.lastCol;
1660 this->time = handler.time;
1661 this->CellWidth = handler.CellWidth;
1662 this->CellHeight = handler.CellHeight;
1663 this->currentRadius = handler.currentRadius;
1664 this->currentPosition = handler.currentPosition;
1665 this->StatBase = handler.StatBase;
1666 this->EdgeX = handler.EdgeX;
1667 this->EdgeY = handler.EdgeY;
1668 this->CalcVariance = handler.CalcVariance;
1669 this->VarxHW = handler.VarxHW;
1670 this->VaryHW = handler.VaryHW;
1671 this->MaxRow = handler.MaxRow;
1672 this->MaxCol = handler.MaxCol;
1673 this->MinRow = handler.MinRow;
1674 this->MinCol = handler.MinCol;
1675 this->lastISAWIntensity = handler.lastISAWIntensity;
1676 this->lastISAWVariance = handler.lastISAWIntensity;
1677 this->back_calc = handler.back_calc;
1678 this->Intensity_calc = handler.Intensity_calc;
1679 this->row_calc = handler.row_calc;
1680 this->col_calc = handler.col_calc;
1681 this->Vx_calc = handler.Vx_calc;
1682 this->Vy_calc = handler.Vy_calc;
1683 this->Vxy_calc = handler.Vxy_calc;
1684 this->case4 = handler.case4;
1685}
1700void DataModeHandler::CalcVariancesFromData(double background, double meanx, double meany, double &Varxx, double &Varxy,
1701 double &Varyy, const std::vector<double> &StatBase) {
1702
1703 double Den = StatBase[IIntensities] - background * StatBase[ISS1];
1704 Varxx = (StatBase[ISSIxx] - 2 * meanx * StatBase[ISSIx] + meanx * meanx * StatBase[IIntensities] -
1705 background * (StatBase[ISSxx] - 2 * meanx * StatBase[ISSx] + meanx * meanx * StatBase[ISS1])) /
1706 Den;
1707
1708 Varyy = (StatBase[ISSIyy] - 2 * meany * StatBase[ISSIy] + meany * meany * StatBase[IIntensities] -
1709 background * (StatBase[ISSyy] - 2 * meany * StatBase[ISSy] + meany * meany * StatBase[ISS1])) /
1710 Den;
1711
1712 Varxy =
1713 (StatBase[ISSIxy] - meanx * StatBase[ISSIy] - meany * StatBase[ISSIx] + meanx * meany * StatBase[IIntensities] -
1714 background *
1715 (StatBase[ISSxy] - meanx * StatBase[ISSy] - meany * StatBase[ISSx] + meanx * meany * StatBase[ISS1])) /
1716 Den;
1717
1718 if (CalcVariances()) {
1719
1720 Varxx = std::min(Varxx, 1.21 * getInitVarx()); // copied from BiVariateNormal
1721 Varxx = std::max(Varxx, .79 * getInitVarx());
1722 Varyy = std::min(Varyy, 1.21 * getInitVary());
1723 Varyy = std::max(Varyy, .79 * getInitVary());
1724 }
1725}
1732std::string DataModeHandler::CalcConstraints(std::vector<std::pair<double, double>> &Bounds, bool CalcVariances) {
1733 double TotIntensity = StatBase[IIntensities];
1734 double ncells = StatBase[ISS1];
1735 double Variance = StatBase[IVariance];
1736 double TotBoundaryIntensities = StatBase[ITotBoundary];
1737 double TotBoundaryVariances = StatBase[IVarBoundary];
1738
1739 double nBoundaryCells = StatBase[INBoundary];
1740 double back = TotBoundaryIntensities / nBoundaryCells;
1741 double backVar = std::max(nBoundaryCells / 50.0, TotBoundaryVariances) / nBoundaryCells / nBoundaryCells;
1742 double IntensVar = Variance + ncells * ncells * backVar;
1743
1744 double relError = .25;
1745
1746 if (back_calc != back)
1747 relError = .45;
1748
1749 int N = NParameters;
1750 if (CalcVariances)
1751 N = N - 3;
1752
1753 double NSigs = NStdDevPeakSpan;
1754 if (back_calc > 0)
1755 NSigs = std::max(NStdDevPeakSpan,
1756 7 - 5 * back_calc / back); // background too high
1757 ostringstream str;
1758
1759 NSigs *= max<double>(1.0, Intensity_calc / (TotIntensity - ncells * back_calc));
1760 str << max<double>(0.0, back_calc - NSigs * (1 + relError) * sqrt(backVar)) << "<Background<"
1761 << (back + NSigs * (1.8 + relError) * sqrt(backVar)) << ","
1762 << max<double>(0.0, Intensity_calc - NSigs * (1 + relError) * sqrt(IntensVar)) << "<Intensity<"
1763 << Intensity_calc + NSigs * (1 + relError) * sqrt(IntensVar);
1764
1765 double min = max<double>(0.0, back_calc - NSigs * (1 + relError) * sqrt(backVar));
1766 double maxx = back + NSigs * (1.8 + relError) * sqrt(backVar);
1767 Bounds.emplace_back(min, maxx);
1768 Bounds.emplace_back(max<double>(0.0, Intensity_calc - NSigs * (1 + relError) * sqrt(IntensVar)),
1769 Intensity_calc + NSigs * (1 + relError) * sqrt(IntensVar));
1770 double relErr1 = relError * .75;
1771 double val = col_calc;
1772 double minn = std::max(MinCol - .5, (1 - relErr1) * val);
1773 maxx = std::min((1 + relErr1) * val, MaxCol + .5);
1774
1775 str << "," << minn << "<"
1776 << "Mcol"
1777 << "<" << maxx;
1778 Bounds.emplace_back(minn, maxx);
1779
1780 val = row_calc;
1781
1782 minn = std::max(MinRow - .5, (1 - relErr1) * val);
1783 maxx = std::min((1 + relErr1) * val, MaxRow + .5);
1784 str << "," << minn << "<"
1785 << "Mrow"
1786 << "<" << maxx;
1787 Bounds.emplace_back(minn, maxx);
1788
1789 if (N >= 5) {
1790 val = Vx_calc;
1791 double valmin = val;
1792 double valmax = val;
1793 if (VarxHW > 0) {
1794 valmin = std::min(val, VarxHW);
1795 valmax = std::max(val, VarxHW);
1796 }
1797
1798 relErr1 *= .6; // Edge peaks: need to restrict sigmas.
1799 str << "," << (1 - relErr1) * valmin << "<"
1800 << "SScol"
1801 << "<" << (1 + relErr1) * valmax;
1802 Bounds.emplace_back((1 - relErr1) * valmin, (1 + relErr1) * valmax);
1803
1804 val = Vy_calc;
1805 valmin = val;
1806 valmax = val;
1807 if (VaryHW > 0) {
1808 valmin = std::min(val, VaryHW);
1809 valmax = std::max(val, VaryHW);
1810 }
1811 str << "," << (1 - relErr1) * valmin << "<"
1812 << "SSrow"
1813 << "<" << (1 + relErr1) * valmax;
1814 Bounds.emplace_back((1 - relErr1) * valmin, (1 + relErr1) * valmax);
1815 }
1816
1817 return str.str();
1818}
1819
1834void IntegratePeakTimeSlices::Fit(MatrixWorkspace_sptr &Data, double &chisqOverDOF, bool &done,
1835 std::vector<string> &names, std::vector<double> &params, std::vector<double> &errs,
1836 double lastRow, double lastCol, double neighborRadius) {
1837
1838 bool CalcVars = m_AttributeValues->CalcVariances();
1839 std::vector<std::pair<double, double>> Bounds;
1840 std::string Constraints = m_AttributeValues->CalcConstraints(Bounds, CalcVars);
1841 auto fit_alg = createChildAlgorithm("Fit");
1842 std::string fun_str = CalculateFunctionProperty_Fit();
1843
1844 std::string SSS(" Fit string ");
1845 SSS += fun_str;
1846 g_log.debug(SSS);
1847 g_log.debug() << " TotCount=" << m_AttributeValues->StatBase[IIntensities] << '\n';
1848
1849 fit_alg->setPropertyValue("Function", fun_str);
1850
1851 fit_alg->setProperty("InputWorkspace", Data);
1852 fit_alg->setProperty("WorkspaceIndex", 0);
1853 fit_alg->setProperty("StartX", 0.0);
1854 fit_alg->setProperty("EndX", 0.0 + static_cast<double>(m_NeighborIDs[1]));
1855 fit_alg->setProperty("MaxIterations", 5000);
1856 fit_alg->setProperty("CreateOutput", true);
1857
1858 fit_alg->setProperty("Output", "out");
1859
1860 fit_alg->setProperty("MaxIterations", 50);
1861
1862 std::string tie = getProperty("Ties");
1863 if (tie.length() > static_cast<size_t>(0))
1864 fit_alg->setProperty("Ties", tie);
1865 if (Constraints.length() > static_cast<size_t>(0))
1866 fit_alg->setProperty("Constraints", Constraints);
1867 try {
1868 fit_alg->executeAsChildAlg();
1869
1870 chisqOverDOF = fit_alg->getProperty("OutputChi2overDoF");
1871 std::string outputStatus = fit_alg->getProperty("OutputStatus");
1872 g_log.debug() << "Chisq/OutputStatus=" << chisqOverDOF << "/" << outputStatus << '\n';
1873
1874 names.clear();
1875 params.clear();
1876 errs.clear();
1877 ITableWorkspace_sptr RRes = fit_alg->getProperty("OutputParameters");
1878 for (int prm = 0; prm < static_cast<int>(RRes->rowCount()) - 1; prm++) {
1879 names.emplace_back(RRes->getRef<string>("Name", prm));
1880 params.emplace_back(RRes->getRef<double>("Value", prm));
1881 double error = RRes->getRef<double>("Error", prm);
1882 errs.emplace_back(error);
1883 }
1884 if (names.size() < 5) {
1885 names.emplace_back(m_ParameterNames[IVXX]);
1886 names.emplace_back(m_ParameterNames[IVYY]);
1887 names.emplace_back(m_ParameterNames[IVXY]);
1888 double Varxx, Varxy, Varyy;
1889 m_AttributeValues->CalcVariancesFromData(params[IBACK], params[IXMEAN], params[IYMEAN], Varxx, Varxy, Varyy,
1890 m_AttributeValues->StatBase);
1891 params.emplace_back(Varxx);
1892 params.emplace_back(Varyy);
1893 params.emplace_back(Varxy);
1894 errs.emplace_back(0);
1895 errs.emplace_back(0);
1896 errs.emplace_back(0);
1897 }
1898
1899 } catch (std::exception &Ex1) // ties or something else went wrong in BivariateNormal
1900 {
1901 done = true;
1902 g_log.error() << "Bivariate Error for PeakNum=" << static_cast<int>(getProperty("PeakIndex")) << ":"
1903 << std::string(Ex1.what()) << '\n';
1904 } catch (...) {
1905 done = true;
1906 g_log.error() << "Bivariate Error A for peakNum=" << static_cast<int>(getProperty("PeakIndex")) << '\n';
1907 }
1908 if (!done) // Bivariate error happened
1909 {
1910
1911 g_log.debug() << " Thru Algorithm: chiSq=" << setw(7) << chisqOverDOF << '\n';
1912 g_log.debug() << " Row,Col Radius=" << lastRow << "," << lastCol << "," << neighborRadius << '\n';
1913
1914 double sqrtChisq = -1;
1915 if (chisqOverDOF >= 0)
1916 sqrtChisq = (chisqOverDOF);
1917
1918 sqrtChisq =
1919 max<double>(sqrtChisq, m_AttributeValues->StatBaseVals(IIntensities) / m_AttributeValues->StatBaseVals(ISS1));
1920 sqrtChisq = SQRT(sqrtChisq);
1921
1922 for (size_t kk = 0; kk < params.size(); kk++) {
1923 g_log.debug() << " Parameter " << setw(8) << names[kk] << " " << setw(8) << params[kk];
1924 // if (names[kk].substr(0, 2) != string("SS"))
1925 g_log.debug() << "(" << setw(8) << (errs[kk] * sqrtChisq) << ")";
1926 if (Bounds.size() > kk) {
1927 pair<double, double> upLow = Bounds[kk];
1928 g_log.debug() << " Bounds(" << upLow.first << "," << upLow.second << ")";
1929 }
1930 g_log.debug() << '\n';
1931 }
1932
1933 double intensity = m_AttributeValues->CalcISAWIntensity(params.data());
1934 g_log.debug() << "IsawIntensity= " << intensity << '\n';
1935 }
1936}
1937
1954void IntegratePeakTimeSlices::PreFit(MatrixWorkspace_sptr &Data, double &chisqOverDOF, bool &done,
1955 std::vector<string> &names, std::vector<double> &params, std::vector<double> &errs,
1956 double lastRow, double lastCol, double neighborRadius) {
1957
1958 double background = m_ParameterValues[IBACK];
1959 int N = 3;
1960 if (background <= 0) {
1961 background = 0;
1962 N = 1;
1963 }
1964 bool CalcVars = m_AttributeValues->CalcVariances();
1965 int NParams = 4;
1966 if (!CalcVars)
1967 NParams += 3;
1968
1969 double minChiSqOverDOF = -1;
1970 double Bestparams[7];
1971 std::string Bestnames[7];
1972 for (int i = 0; i < N; i++) {
1973 Fit(Data, chisqOverDOF, done, names, params, errs, lastRow, lastCol, neighborRadius);
1974 g_log.debug() << "-----------------------" << i << "--------------------------\n";
1975 if ((minChiSqOverDOF < 0 || chisqOverDOF < minChiSqOverDOF) && (chisqOverDOF > 0) && !done) {
1976 for (int j = 0; j < NParams; j++) {
1977 Bestparams[j] = m_ParameterValues[j];
1978 Bestnames[j] = m_ParameterNames[j];
1979 }
1980 minChiSqOverDOF = chisqOverDOF;
1981 }
1982
1983 // Next round, reduce background
1984 background = background / 2;
1985 if (i + 1 == N - 1)
1986 background = 0;
1987
1988 std::vector<double> prms = m_AttributeValues->GetParams(background);
1989
1990 for (int j = 0; j < NParams; j++)
1991 m_ParameterValues[j] = prms[j];
1992 }
1993 vector<std::string> ParNames(m_ParameterNames, m_ParameterNames + NParams);
1994 for (int i = 0; i < NParams; i++) {
1995 int k = findNameInVector(Bestnames[i], ParNames);
1996 if (k >= 0 && k < NParams)
1997 m_ParameterValues[k] = Bestparams[k];
1998 }
1999
2000 Fit(Data, chisqOverDOF, done, names, params, errs, lastRow, lastCol, neighborRadius);
2001}
2002
2013bool IntegratePeakTimeSlices::isGoodFit(std::vector<double> const &params, std::vector<double> const &errs,
2014 std::vector<std::string> const &names, double chisqOverDOF) {
2015 int Ibk = findNameInVector("Background", names);
2016 if (Ibk < 0)
2017 throw std::runtime_error("Irrecoverable inconsistency found. The index for the "
2018 "parameter 'Background' is lower than zero.");
2019
2020 int IIntensity = findNameInVector("Intensity", names);
2021 if (IIntensity < 0)
2022 throw std::runtime_error("Irrecoverable inconsistency found. The index for the "
2023 "parameter 'Intensity' is lower than zero.");
2024
2025 if (chisqOverDOF < 0) {
2026
2027 g_log.debug() << " Bad Slice- negative chiSq= " << chisqOverDOF << '\n';
2028 ;
2029 return false;
2030 }
2031
2032 int NBadEdgeCells = getProperty("NBadEdgePixels");
2033 NBadEdgeCells = static_cast<int>(.6 * NBadEdgeCells);
2034 if (params[IXMEAN] < NBadEdgeCells || params[IYMEAN] < NBadEdgeCells || params[IXMEAN] > m_NCOLS - NBadEdgeCells ||
2035 params[IYMEAN] > m_NROWS - NBadEdgeCells)
2036 return false;
2037
2038 auto ncells = static_cast<int>(m_AttributeValues->StatBaseVals(ISS1));
2039
2040 if (m_AttributeValues->StatBaseVals(IIntensities) <= 0 ||
2041 (m_AttributeValues->StatBaseVals(IIntensities) - params[Ibk] * ncells) <= 0) {
2042
2043 g_log.debug() << " Bad Slice. Negative Counts= "
2044 << m_AttributeValues->StatBaseVals(IIntensities) - params[Ibk] * ncells << '\n';
2045 ;
2046 return false;
2047 }
2048
2049 double x = params[IIntensity] / (m_AttributeValues->StatBaseVals(IIntensities) - params[Ibk] * ncells);
2050
2051 if ((x < MinGoodRatioFitvsExpIntenisites || x > MaxGoodRatioFitvsExpIntenisites) &&
2052 !m_EdgePeak) // The fitted intensity should be close to tot intensity -
2053 // background
2054 {
2055 g_log.debug() << " Bad Slice. Fitted Intensity & Observed "
2056 "Intensity(-back) too different. ratio="
2057 << x << '\n';
2058
2059 return false;
2060 }
2061
2062 bool GoodNums = true;
2063 bool paramBad = false;
2064 auto BadParamNum = static_cast<size_t>(-1);
2065 for (size_t i = 0; i < errs.size(); i++)
2066 if (errs[i] != errs[i]) {
2067 GoodNums = false;
2068 paramBad = false;
2069 BadParamNum = i;
2070 } else if (errs[i] < 0) {
2071 GoodNums = false;
2072 paramBad = false;
2073 BadParamNum = i;
2074 } else if (params[i] != params[i]) {
2075 GoodNums = false;
2076 paramBad = true;
2077 BadParamNum = i;
2078 }
2079
2080 if (!GoodNums) {
2081 std::string obj = " parameter ";
2082 if (!paramBad)
2083 obj = " error ";
2084 g_log.debug() << " Bad Slice." << obj << BadParamNum << " is not a number\n";
2085 return false;
2086 }
2087
2088 GoodNums = true;
2089
2090 std::string Err("back ground is negative");
2091 if (params[Ibk] < -.002)
2092 GoodNums = false;
2093
2094 if (GoodNums)
2095 Err = "Intensity is negative";
2096 if (params[IIntensity] < 0)
2097 GoodNums = false;
2098
2099 double IsawIntensity = m_AttributeValues->CalcISAWIntensity(params.data());
2100 double IsawVariance = m_AttributeValues->CalcISAWIntensityVariance(params.data(), errs.data(), chisqOverDOF);
2101 if (GoodNums)
2102 Err = "Isaw Variance is negative";
2103 if (IsawVariance > 0) {
2104 if (GoodNums)
2105 Err = "I/sigI > 3";
2106 if (IsawIntensity * IsawIntensity / IsawVariance < MinGoodIoverSigI * MinGoodIoverSigI)
2107 GoodNums = false;
2108 } else
2109 GoodNums = false;
2110
2111 if (!GoodNums)
2112
2113 {
2114 g_log.debug() << Err << '\n';
2115
2116 return false;
2117 }
2118
2119 // Check weak peak. Max theoretical height should be more than 3
2120
2121 double maxPeakHeightTheoretical =
2122 params[ITINTENS] / 2 / M_PI / sqrt(params[IVXX] * params[IVYY] - params[IVXY] * params[IVXY]);
2123
2124 double AvHeight =
2125 m_AttributeValues->StatBaseVals(IIntensities) / m_AttributeValues->StatBaseVals(ISS1) - params[IBACK];
2126
2127 if (maxPeakHeightTheoretical < 2 * AvHeight || AvHeight < 0 || maxPeakHeightTheoretical < 0) {
2128
2129 g_log.debug() << " Bad Slice. Peak too small= " << maxPeakHeightTheoretical << "/" << AvHeight << '\n';
2130 return false;
2131 }
2132
2133 double Nrows = std::max(m_AttributeValues->StatBase[INRows], m_AttributeValues->StatBase[INCol]);
2134 if (maxPeakHeightTheoretical < 1 && (params[IVXX] > Nrows * Nrows / 4 || params[IVYY] > Nrows * Nrows / 4)) {
2135 g_log.debug() << "Peak is too flat \n";
2136 return false;
2137 }
2138
2139 // Exponential too steep, i.e. intensities at pixels 1 from center are <3*
2140 // intensity center
2141 if (params[IVXX] + params[IVYY] > 2.6 * (params[IVXX] * params[IVYY] - params[IVXY] * params[IVXY])) {
2142 g_log.debug() << " Bad Slice. Too steep of an exponential\n";
2143 return false;
2144 }
2145
2146 return true;
2147}
2148
2155bool DataModeHandler::IsEnoughData(const double *ParameterValues, Kernel::Logger & /*unused*/) {
2156 // Check if flat
2157 double Varx, Vary, Cov;
2158
2159 if (StatBase.empty())
2160 return false;
2161
2162 double ncells = static_cast<int>(StatBase[IIntensities]);
2163 if (ncells <= 0)
2164 return false;
2165
2166 double meanx = StatBase[ISSIx] / ncells;
2167 double meany = StatBase[ISSIy] / ncells;
2168
2169 if (!CalcVariances()) {
2170 Varx = ParameterValues[IVXX];
2171 Vary = ParameterValues[IVYY];
2172 Cov = ParameterValues[IVXY];
2173
2174 } else
2175 CalcVariancesFromData(ParameterValues[0], meanx, meany, Varx, Cov, Vary, StatBase);
2176
2177 if (Varx < MinVariationInXYvalues || Vary < MinVariationInXYvalues) // All data essentially the same.
2178 return false;
2179
2180 if (Cov * Cov > MaxCorrCoeffinXY * Varx * Vary) // All data on a obtuse line
2181 return false;
2182
2183 return true;
2184}
2185
2198double IntegratePeakTimeSlices::CalculateIsawIntegrateError(const double background, const double backError,
2199 const double ChiSqOverDOF, const double TotVariance,
2200 const int ncells) {
2201
2202 double B = TotVariance / ncells;
2203 if (B < ChiSqOverDOF)
2204 B = ChiSqOverDOF;
2205
2206 double Variance = TotVariance + (backError * backError * B) * ncells * ncells + background * ncells;
2207
2208 return SQRT(Variance);
2209}
2210
2217 // TabWS->setName("Log Table");
2218 TabWS->addColumn("double", "Time");
2219 TabWS->addColumn("double", "Channel");
2220 TabWS->addColumn("double", "Background");
2221 TabWS->addColumn("double", "Intensity");
2222 TabWS->addColumn("double", "Mcol");
2223 TabWS->addColumn("double", "Mrow");
2224 TabWS->addColumn("double", "SScol");
2225 TabWS->addColumn("double", "SSrow");
2226 TabWS->addColumn("double", "SSrc");
2227 TabWS->addColumn("double", "NCells");
2228 TabWS->addColumn("double", "ChiSqrOverDOF");
2229 TabWS->addColumn("double", "TotIntensity");
2230 TabWS->addColumn("double", "BackgroundError");
2231 TabWS->addColumn("double", "FitIntensityError");
2232 TabWS->addColumn("double", "ISAWIntensity");
2233 TabWS->addColumn("double", "ISAWIntensityError");
2234 TabWS->addColumn("double", "TotalBoundary");
2235 TabWS->addColumn("double", "NBoundaryCells");
2236 TabWS->addColumn("double", "Start Row");
2237 TabWS->addColumn("double", "End Row");
2238 TabWS->addColumn("double", "Start Col");
2239 TabWS->addColumn("double", "End Col");
2240 TabWS->addColumn("double", "TotIntensityError");
2241 TabWS->addColumn("str", "SpecIDs");
2242}
2243
2262 std::vector<double> const &params, std::vector<double> const &errs,
2263 std::vector<std::string> const &names, const double Chisq,
2264 const double time, string spec_idList) {
2265 int Ibk = findNameInVector("Background", names);
2266 int IIntensity = findNameInVector("Intensity", names);
2267 int IVx = findNameInVector("SScol", names);
2268 int IVy = findNameInVector("SSrow", names);
2269 int IVxy = findNameInVector("SSrc", names);
2270 int Irow = findNameInVector("Mrow", names);
2271 int Icol = findNameInVector("Mcol", names);
2272
2273 if (Ibk < 0 || IIntensity < 0 || IVx < 0 || IVy < 0 || IVxy < 0 || Irow < 0 || Icol < 0) {
2274 throw std::runtime_error("Inconsistency found when updating output "
2275 "workspace. None of the indices for the "
2276 "parameters 'Background', 'Intensity', 'SScol', "
2277 "'SSrow', 'SSrc', 'Mrow', 'Mcol' can be "
2278 "negative.");
2279 }
2280
2281 int newRowIndex = 0;
2282
2283 if (dir > 0)
2284 newRowIndex = static_cast<int>(TabWS->rowCount());
2285
2286 auto TableRow = static_cast<int>(TabWS->insertRow(newRowIndex));
2287
2288 auto ncells = static_cast<int>(m_AttributeValues->StatBaseVals(ISS1));
2289 double chisq = max<double>(Chisq, m_AttributeValues->StatBaseVals(IIntensities) / max<int>(ncells, 1));
2290
2291 TabWS->getRef<double>(std::string("Background"), TableRow) = params[Ibk];
2292 TabWS->getRef<double>(std::string("Channel"), TableRow) = chan;
2293
2294 TabWS->getRef<double>(std::string("Intensity"), TableRow) = params[IIntensity];
2295 TabWS->getRef<double>(std::string("FitIntensityError"), TableRow) = errs[IIntensity] * sqrt(chisq);
2296 TabWS->getRef<double>(std::string("Mcol"), TableRow) = params[Icol];
2297 TabWS->getRef<double>(std::string("Mrow"), TableRow) = params[Irow];
2298
2299 TabWS->getRef<double>(std::string("SScol"), TableRow) = params[IVx];
2300 TabWS->getRef<double>(std::string("SSrow"), TableRow) = params[IVy];
2301
2302 TabWS->getRef<double>(std::string("SSrc"), TableRow) = params[IVxy];
2303 TabWS->getRef<double>(std::string("NCells"), TableRow) = ncells;
2304 TabWS->getRef<double>(std::string("ChiSqrOverDOF"), TableRow) = chisq;
2305
2306 TabWS->getRef<double>(std::string("TotIntensity"), TableRow) = m_AttributeValues->StatBaseVals(IIntensities);
2307 TabWS->getRef<double>(std::string("BackgroundError"), TableRow) = errs[Ibk] * SQRT(chisq);
2308 TabWS->getRef<double>(std::string("ISAWIntensity"), TableRow) = m_AttributeValues->CalcISAWIntensity(params.data());
2309 // m_AttributeValues->StatBaseVals(IIntensities)
2310 // -
2311 // params[Ibk]
2312 // * ncells;
2313
2314 TabWS->getRef<double>(std::string("ISAWIntensityError"), TableRow) =
2315 sqrt(m_AttributeValues->CalcISAWIntensityVariance(params.data(), errs.data(), Chisq));
2316
2317 // CalculateIsawIntegrateError(
2318 // params[Ibk], errs[Ibk], chisq, m_AttributeValues->StatBaseVals(IVariance),
2319 // ncells);
2320
2321 TabWS->getRef<double>(std::string("Time"), TableRow) = time;
2322
2323 TabWS->getRef<double>(std::string("TotalBoundary"), TableRow) = m_AttributeValues->StatBaseVals(ITotBoundary);
2324 TabWS->getRef<double>(std::string("NBoundaryCells"), TableRow) = m_AttributeValues->StatBaseVals(INBoundary);
2325
2326 TabWS->getRef<double>(std::string("Start Row"), TableRow) = m_AttributeValues->StatBaseVals(IStartRow);
2327 TabWS->getRef<double>(std::string("End Row"), TableRow) =
2328 m_AttributeValues->StatBaseVals(IStartRow) + m_AttributeValues->StatBaseVals(INRows) - 1;
2329
2330 TabWS->getRef<double>(std::string("Start Col"), TableRow) = m_AttributeValues->StatBaseVals(IStartCol);
2331 TabWS->getRef<double>(std::string("End Col"), TableRow) =
2332 m_AttributeValues->StatBaseVals(IStartCol) + m_AttributeValues->StatBaseVals(INCol) - 1;
2333 TabWS->getRef<double>(std::string("TotIntensityError"), TableRow) = SQRT(m_AttributeValues->StatBaseVals(IVariance));
2334 TabWS->getRef<string>(std::string("SpecIDs"), TableRow) = std::move(spec_idList);
2335
2336 return newRowIndex;
2337}
2338
2353void IntegratePeakTimeSlices::updatePeakInformation(std::vector<double> const &params, std::vector<double> const &errs,
2354 std::vector<std::string> const &names, double &TotVariance,
2355 double &TotIntensity, double const TotSliceIntensity,
2356 double const TotSliceVariance, double const chisqdivDOF,
2357 const int ncells) {
2358 UNUSED_ARG(TotSliceIntensity);
2359 UNUSED_ARG(TotSliceVariance);
2360 UNUSED_ARG(names);
2361 UNUSED_ARG(ncells);
2362
2363 double err = 0;
2364 double intensity = 0;
2365
2366 err = m_AttributeValues->CalcISAWIntensityVariance(params.data(), errs.data(), chisqdivDOF);
2367
2368 intensity = m_AttributeValues->CalcISAWIntensity(params.data());
2369 TotIntensity += intensity;
2370
2371 TotVariance += err;
2372 g_log.debug() << "TotIntensity/TotVariance=" << TotIntensity << "/" << TotVariance << '\n';
2373}
2374
2384 if (!CalcVariance)
2385 return false;
2386
2388 return !isEdgePeak(param, 7);
2389}
2390
2400double DataModeHandler::CalcISAWIntensity(const double *params) {
2401
2402 double ExperimentalIntensity = StatBase[IIntensities] - params[IBACK] * StatBase[ISS1];
2403
2404 double r = CalcSampleIntensityMultiplier(params);
2405
2406 double alpha = .5 * (r - 1.0);
2407 alpha = std::min(1.0, alpha);
2408
2409 lastISAWIntensity = ExperimentalIntensity * r; //*( 1-alpha )+ alpha * FitIntensity;
2410 return lastISAWIntensity;
2411}
2422double DataModeHandler::CalcISAWIntensityVariance(const double *params, const double *errs, double chiSqOvDOF) {
2423
2424 auto ncells = static_cast<int>(StatBase[ISS1]);
2425 double B = StatBase[IVariance] / ncells;
2426 if (B < chiSqOvDOF)
2427 B = chiSqOvDOF;
2428
2429 double ExperimVar = StatBase[IVariance];
2430 double IntensityBackError = errs[IBACK] * sqrt(B);
2431
2432 ExperimVar += IntensityBackError * IntensityBackError * ncells * ncells + params[IBACK] * ncells;
2433
2434 double r = CalcSampleIntensityMultiplier(params);
2435 double alpha = .5 * (r - 1.0);
2436 alpha = std::min(1.0, alpha);
2437
2438 lastISAWVariance = ExperimVar * r * r; //*( 1 - alpha ) + alpha * FitVar;
2439 return lastISAWVariance;
2440}
2441
2450double DataModeHandler::CalcSampleIntensityMultiplier(const double *params) const {
2451 auto MinRow = static_cast<int>(StatBase[IStartRow]);
2452 int MaxRow = MinRow + static_cast<int>(StatBase[INRows]) - 1;
2453 auto MinCol = static_cast<int>(StatBase[IStartCol]);
2454 int MaxCol = MinCol + static_cast<int>(StatBase[INCol]) - 1;
2455 double r = 1;
2456
2457 if (params[IVXX] <= 0 || params[IVYY] <= 0)
2458 return 1.0;
2459
2460 // NstdX iand NstdY are the number of 1/4 standard deviations. Elements of
2461 // probs are in
2462 // 1/4 standard deviations
2463 double NstdX = 4 * min<double>(params[IXMEAN] - MinCol, MaxCol - params[IXMEAN]) / sqrt(params[IVXX]);
2464
2465 double NstdY = 4 * min<double>(params[IYMEAN] - MinRow, MaxRow - params[IYMEAN]) / sqrt(params[IVYY]);
2466
2467 double sgn = 1;
2468 if (NstdX < 0) {
2469 sgn = -1;
2470 }
2471 double P = 1;
2472 if (sgn * NstdX < 9) {
2473 auto xx = static_cast<int>(sgn * NstdX);
2474 double a = probs[xx];
2475 double b = 1;
2476 if (xx + 1 <= 8)
2477 b = probs[xx + 1];
2478 P = a + (b - a) * (sgn * NstdX - xx);
2479 }
2480 if (NstdX >= 7.5)
2481 r = 1.0;
2482 else if (sgn > 0)
2483 r = 1 / P;
2484 else
2485 r = 1 / (1 - P);
2486
2487 if (NstdY < 0) {
2488 sgn = -1;
2489 }
2490 P = 1;
2491 if (sgn * NstdY < 9) {
2492 auto xx = static_cast<int>(sgn * NstdY);
2493 double a = probs[xx];
2494 double b = 1;
2495 if (xx + 1 <= 8)
2496 b = probs[xx + 1];
2497 P = a + (b - a) * (sgn * NstdY - xx);
2498 }
2499 if (NstdY >= 7.5)
2500 r *= 1.0;
2501 else if (sgn > 0)
2502 r *= 1 / P;
2503 else
2504 r *= 1 / (1 - P);
2505
2506 r = std::max(r, 1.0);
2507 return r;
2508}
2509
2510} // namespace Mantid::Crystal
2511// Attr indicies
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double sigma
Definition: GetAllEi.cpp:156
double error
Definition: IndexPeaks.cpp:133
#define IVYY
#define IVarBoundary
#define IStartRow
#define ISS1
#define ISSIxx
#define ISSy
#define ISSyy
#define IBACK
#define INBoundary
#define ISSIx
#define INCol
#define IYMEAN
#define ISSIxy
#define ISSIy
#define ISSIyy
#define ITINTENS
#define NAttributes
#define IXMEAN
#define INRows
#define ITotBoundary
#define ISSx
#define IVariance
#define IVXX
#define ISSxy
#define IIntensities
#define IStartCol
#define ISSxx
#define NParameters
#define IVXY
#define fabs(x)
Definition: Matrix.cpp:22
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
double obj
the value of the quadratic function
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
Integrates each time slice using the BivariateNormal formula, adding the results to the peak object.
void setHeightHalfWidthInfo(const MantidVec &xvals, const MantidVec &yvals, const MantidVec &counts)
For edge peaks, the sample standard deviations do not work.
double CalcISAWIntensity(const double *params)
Calculates the Intensity designed for Edge Peaks.
double CalcISAWIntensityVariance(const double *params, const double *errs, double chiSqOvDOF)
Calculates the Error in the Intensity designed for Edge Peaks.
std::vector< double > InitValues(double Varx, double Vary, double b)
Returns init values with background and variances replaced by arguments.
bool isEdgePeak(const double *params, int nparams)
Determines if a Peak is an edge peak.
bool CalcVariances()
Determines whether the Variances can be calculated.
double getNewRCRadius()
Calculates the new radius for neighborhoods so as to include almost all of a peak.
void CalcVariancesFromData(double background, double meanx, double meany, double &Varxx, double &Varxy, double &Varyy, const std::vector< double > &StatBase)
Utility method to calculate variances from data given background and means.
double CalcSampleIntensityMultiplier(const double *params) const
For Edge Peaks.
std::vector< double > GetParams(double b)
Calculates the initial values of the parameters given background b.
bool IsEnoughData(const double *ParameterValues, Kernel::Logger &)
Calculates if there is enough data to for there to be a peak.
std::string CalcConstraints(std::vector< std::pair< double, double > > &Bounds, bool CalcVariances)
Calculates the string form of the constraints to be sent to the Fit Algorithm and also saves it in a ...
bool setStatBase(std::vector< double > const &StatBase)
Sets the Accumulated data values into this class, then updates other information like initial values.
void SetUpData1(API::MatrixWorkspace_sptr &Data, API::MatrixWorkspace_const_sptr const &inpWkSpace, const int chanMin, const int chanMax, double Radius, const Kernel::V3D &CentPos, std::string &spec_idList)
Prepares the data for futher analysis adding meta data and marking data on the edges of detectors.
void FindPlane(Kernel::V3D &center, Kernel::V3D &xvec, Kernel::V3D &yvec, double &ROW, double &COL, int &NROWS, int &NCOLS, double &pixWidthx, double &pixHeighty, DataObjects::Peak const &peak) const
For NonFlat banks, this determines the data of a small planar region approximating the instrument clo...
int CalculateTimeChannelSpan(Geometry::IPeak const &peak, const double dQ, const Mantid::HistogramData::HistogramX &X, const int specNum, int &Centerchan)
Calculates the span of channels needed to encompass all data around the peak with Q values within dQ ...
double m_COL
for Describing the Column(or 0) describing the center of the
double CalculatePositionSpan(DataObjects::Peak const &peak, const double dQ)
Calculates the span in rows and columns needed to include all data within dQ of the specified peak.
double CalculateIsawIntegrateError(const double background, const double backError, const double ChiSqOverDOF, const double TotVariance, const int ncells)
Calculates the error in integration closest to the latest ISAW calculations.
void exec() override
Executes this algorithm.
std::string CalculateFunctionProperty_Fit()
Calculates the string for the Function Property of the Fit Algorithm.
Kernel::V3D m_center
for Describing the Plane at the Peak
Kernel::V3D m_yvec
for Describing the Plane at the Peak
Kernel::V3D m_xvec
for Describing the Plane at the Peak
void SetUpData(API::MatrixWorkspace_sptr &Data, API::MatrixWorkspace_const_sptr const &inpWkSpace, const std::shared_ptr< Geometry::IComponent > &comp, const int chanMin, const int chanMax, double CentX, double CentY, Kernel::V3D &CentNghbr, double &neighborRadius, double Radius, std::string &spec_idList)
Initial phase at converting Detector data to workspace data that will be sent to the Fit Function,...
std::shared_ptr< DataModeHandler > m_AttributeValues
void init() override
Virtual method - must be overridden by concrete algorithm.
double m_ROW
for Describing the Row(or 0) describing the center of the Peak
double m_cellHeight
for Describing the Plane at the Peak
void Fit(API::MatrixWorkspace_sptr &Data, double &chisqOverDOF, bool &done, std::vector< std::string > &names, std::vector< double > &params, std::vector< double > &errs, double lastRow, double lastCol, double neighborRadius)
Sets up data for the Fit Algorithm call and invokes it.
bool getNeighborPixIDs(const std::shared_ptr< Geometry::IComponent > &comp, const Kernel::V3D &Center, double &Radius, int *&ArryofID)
Finds all neighbors within a given Radius of the Center on the given component.
int findNameInVector(std::string const &oneName, std::vector< std::string > const &nameList)
Utility to find a name in a vector of strings.
double m_R0
for Weak Peaks, these can be set using info from close
void updateStats(const double intensity, const double variance, const double row, const double col, std::vector< double > &StatBase)
Updates the cumulative statistics for the data being considered.
void InitializeColumnNamesInTableWorkspace(DataObjects::TableWorkspace_sptr &TabWS)
Initializes the column names in the output table workspace.
bool updateNeighbors(const std::shared_ptr< Geometry::IComponent > &comp, const Kernel::V3D &CentPos, const Kernel::V3D &oldCenter, double NewRadius, double &neighborRadius)
Checks and updates if needed the list of m_NeighborIDs.
void PreFit(API::MatrixWorkspace_sptr &Data, double &chisqOverDOF, bool &done, std::vector< std::string > &names, std::vector< double > &params, std::vector< double > &errs, double lastRow, double lastCol, double neighborRadius)
Tests several starting points in the Marquardt algorithm then calls Fit.
int UpdateOutputWS(DataObjects::TableWorkspace_sptr &TabWS, const int dir, const double chan, std::vector< double > const &params, std::vector< double > const &errs, std::vector< std::string > const &names, const double Chisq, const double time, std::string spec_idList)
Updates the information in the output OutputWorkspace.
bool isGoodFit(std::vector< double > const &params, std::vector< double > const &errs, std::vector< std::string > const &names, double chisqOverDOF)
Determines if the list of parameters and errors represent a "good" fit.
int findTimeChannel(const Mantid::HistogramData::HistogramX &X, const double time)
Finds the time channel with the given time in.
void updatePeakInformation(std::vector< double > const &params, std::vector< double > const &errs, std::vector< std::string > const &names, double &TotVariance, double &TotIntensity, double const TotSliceIntensity, double const TotSliceVariance, double const chisqdivDOF, const int ncells)
Updates m_AttributeValues with this peak information from this time slice.
A generic fitting algorithm.
Definition: Fit.h:78
Structure describing a single-crystal peak.
Definition: Peak.h:34
int getCol() const override
For RectangularDetectors only, returns the column (x) of the pixel of the detector or -1 if not found...
Definition: Peak.cpp:341
Geometry::Instrument_const_sptr getInstrument() const
Return a shared ptr to the instrument for this peak.
Definition: Peak.cpp:326
std::string getBankName() const
Find the name of the bank that is the parent of the detector.
Definition: Peak.cpp:353
int getRow() const override
For RectangularDetectors only, returns the row (y) of the pixel of the detector or -1 if not found.
Definition: Peak.cpp:332
Mantid::Kernel::V3D getQLabFrame() const override
Return the Q change (of the lattice, k_i - k_f) for this peak.
Definition: Peak.cpp:441
Geometry::IDetector_const_sptr getDetector() const
Return a shared ptr to the detector at center of peak.
Definition: Peak.cpp:323
virtual Mantid::Kernel::V3D getDetPos() const
Return the detector position vector.
Definition: Peak.cpp:710
int getDetectorID() const
Get the ID of the detector at the center of the peak
Definition: Peak.cpp:271
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
Definition: BoundingBox.h:34
bool isPointInside(const Kernel::V3D &point) const
Is the given point within the bounding box?
Definition: BoundingBox.cpp:28
double xMax() const
Return the maximum value of X.
Definition: BoundingBox.h:80
double zMin() const
Return the minimum value of Z.
Definition: BoundingBox.h:86
double zMax() const
Return the maximum value of Z.
Definition: BoundingBox.h:88
double yMax() const
Return the maximum value of Y.
Definition: BoundingBox.h:84
double xMin() const
Return the minimum value of X.
Definition: BoundingBox.h:78
double yMin() const
Return the minimum value of Y.
Definition: BoundingBox.h:82
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual double getTOF() const =0
virtual Mantid::Kernel::V3D getQLabFrame() const =0
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
Class for quaternions.
Definition: Quat.h:39
void inverse()
Inverse a quaternion (in the sense of rotation inversion)
Definition: Quat.cpp:376
void rotate(V3D &) const
Rotate a vector.
Definition: Quat.cpp:397
std::vector< double > getRotation(bool check_normalisation=false, bool throw_on_errors=false) const
returns the rotation matrix defined by this quaternion as an 9-point
Definition: Quat.cpp:453
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
Definition: V3D.h:274
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
void setZ(const double zz) noexcept
Set is z position.
Definition: V3D.h:230
double norm() const noexcept
Definition: V3D.h:263
void setX(const double xx) noexcept
Set is x position.
Definition: V3D.h:218
void setY(const double yy) noexcept
Set is y position.
Definition: V3D.h:224
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
Definition: IDetector.h:102
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54