IntegratePeakTimeSlices.cpp 88.6 KB
Newer Older
1
2
3
4
5
/*
 * IntegratePeakTimeSlices.cpp
 *
 *  Created on: May 5, 2011
 *      Author: ruth
6
7
 *
 *
8
9
 *
 */
10
#include "MantidCrystal/IntegratePeakTimeSlices.h"
11
#include "MantidAPI/ConstraintFactory.h"
12
#include "MantidAPI/IFunction.h"
13
14
15
16
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/IFuncMinimizer.h"
#include "MantidAPI/FuncMinimizerFactory.h"
#include "MantidAPI/FunctionFactory.h"
17
#include "MantidAPI/WorkspaceFactory.h"
18
#include "MantidDataObjects/Workspace2D.h"
19
//#include "MantidGeometry/Surfaces/Surface.h"
20

21
//#include <boost/random/poisson_distribution.hpp>
22
23
24
25
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
26
using namespace std;
27
28
namespace Mantid {
namespace Crystal {
29

30
DECLARE_ALGORITHM(IntegratePeakTimeSlices)
31

32
// Attr, m_AttributeValues, and StatBase indicies
33
34
35
36
37
38
39
40
41
42
43
44
#define IStartRow 0
#define IStartCol 1
#define INRows 2
#define INCol 3
#define ISSIxx 4
#define ISSIyy 5
#define ISSIxy 6
#define ISSxx 7
#define ISSyy 8
#define ISSxy 9
#define ISSIx 10
#define ISSIy 11
45
46
#define ISSx 12
#define ISSy 13
47
#define IIntensities 14
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#define ISS1 15
#define IVariance 16
#define ITotBoundary 17
#define INBoundary 18
#define IVarBoundary 19
#define NAttributes 20

// Parameter indicies
#define IBACK 0
#define ITINTENS 1
#define IXMEAN 2
#define IYMEAN 3
#define IVXX 4
#define IVYY 5
#define IVXY 6
// TODO: Edge Peaks-Return NoPeak(Intensity=0,variance=0) if any center on any
// slice is out of bounds
// TODO: Calc ratio for edge peaks, should use the slant of bivariate normal
// instead of assuming
67
//        distribution axes are lined up with the row and col vectors
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#define NParameters 7

namespace {
//           # std sigs  0    .25     .5     .75    1     1.25   1.5    2 2.5
const double probs[9] = {.5f,    .5987f, .6915f, .7734f, .8413f,
                         .8944f, .9322f, .9599f, .9772f};

const int MinRowColSpan = 6;
const int MaxRowColSpan = 36;
const int MinTimeSpan = 3;
const double NeighborhoodRadiusDivPeakRadius = 1.5;
const double MaxNeighborhoodRadius = 10;
const double NStdDevPeakSpan = 2;
const double MaxGoodRatioFitvsExpIntenisites = 2.5;
const double MinGoodRatioFitvsExpIntenisites = .25;
const double MinGoodIoverSigI = 3.0;
const double MinVariationInXYvalues = .6; // Peak spans one pixel only
const double MaxCorrCoeffinXY = .9;       // otherwise all data on one line
}

88
89
90
91
IntegratePeakTimeSlices::IntegratePeakTimeSlices()
    : Algorithm(), m_R0(-1), m_ROW(0.), m_COL(0.), m_cellWidth(0.),
      m_cellHeight(0.), m_NROWS(0), m_NCOLS(0) {
  m_debug = false;
92

93
  if (m_debug)
94
    g_log.setLevel(7);
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
  m_EdgePeak = false;
  m_NeighborIDs = new int[3];
  m_NeighborIDs[0] = 3;
  m_NeighborIDs[1] = 2;
  m_AttributeNames[0] = "StartRow";
  m_AttributeNames[1] = "StartCol";
  m_AttributeNames[2] = "NRows";
  m_AttributeNames[3] = "NCols";
  m_AttributeNames[4] = "SSIxx";
  m_AttributeNames[5] = "SSIyy";
  m_AttributeNames[6] = "SSIxy";
  m_AttributeNames[7] = "SSxx";
  m_AttributeNames[8] = "SSyy";
  m_AttributeNames[9] = "SSxy";
  m_AttributeNames[10] = "SSIx";
  m_AttributeNames[11] = "SSIy";
  m_AttributeNames[12] = "SSx";
  m_AttributeNames[13] = "SSy";
  m_AttributeNames[14] = "Intensities";
  m_AttributeNames[15] = " SS1";
  m_AttributeNames[16] = "Variance";
  m_AttributeNames[17] = "TotBoundary";
  m_AttributeNames[18] = "NBoundary";
  m_AttributeNames[19] = "VarianceBoundary";

  m_ParameterNames[0] = "Background";
  m_ParameterNames[1] = "Intensity";
  m_ParameterNames[2] = "Mcol";
  m_ParameterNames[3] = "Mrow";
  m_ParameterNames[4] = "SScol";
  m_ParameterNames[5] = "SSrow";
  m_ParameterNames[6] = "SSrc";
127
128

  // for (int i = 0; i < NAttributes; i++)
129
  //   m_AttributeValues[i] = 0;
130

131
132
  for (double &m_ParameterValue : m_ParameterValues)
    m_ParameterValue = 0;
Hahn, Steven's avatar
Hahn, Steven committed
133
  std::fill(m_ParameterValues.begin(), m_ParameterValues.end(), 0.0);
134
135
136
137
138
139
140
141
}

double SQRT(double v) {
  if (v < 0)
    return -1;
  return sqrt(v);
}
/// Destructor
142
IntegratePeakTimeSlices::~IntegratePeakTimeSlices() { delete[] m_NeighborIDs; }
143
144

void IntegratePeakTimeSlices::init() {
145
146
  declareProperty(Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      "InputWorkspace", "", Direction::Input),
147
148
                  "A 2D workspace with X values of time of flight");

149
150
  declareProperty(Kernel::make_unique<WorkspaceProperty<TableWorkspace>>(
                      "OutputWorkspace", "", Direction::Output),
151
152
                  "Name of the output table workspace with Log info");

153
154
155
  declareProperty(Kernel::make_unique<WorkspaceProperty<PeaksWorkspace>>(
                      "Peaks", "", Direction::Input),
                  "Workspace of Peaks");
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191

  declareProperty("PeakIndex", 0,
                  "Index of peak in PeaksWorkspace to integrate");

  declareProperty(
      "PeakQspan", .06,
      "Max magnitude of Q of Peak to Q of Peak Center, where mod(Q)=1/d");

  declareProperty(
      "CalculateVariances", true,
      "Calc (co)variances given parameter values versus fit (co)Variances ");

  declareProperty(
      "Ties", "",
      "Tie parameters(Background,Intensity, Mrow,...) to values/formulas.");

  declareProperty("NBadEdgePixels", 0, "Number of  bad Edge Pixels");

  declareProperty("Intensity", 0.0, "Peak Integrated Intensity",
                  Direction::Output);

  declareProperty("SigmaIntensity", 0.0, "Peak Integrated Intensity Error",
                  Direction::Output);
}

/**
 * Executes this algorithm
 *
 * Integrates one peak
 * -First attempts to find row/col extents, time extents to fully get most of
 *the peak
 * -Integrate each time slice.
 * -Report Results
 */
void IntegratePeakTimeSlices::exec() {
  time_t seconds1;
192

193
  seconds1 = time(nullptr);
194

195
  double dQ = getProperty("PeakQspan");
196

197
  g_log.debug("------------------Start Peak Integrate-------------------");
198

199
200
201
202
203
204
  if (dQ <= 0) {
    g_log.error("Negative PeakQspans are not allowed. Use .17/G where G is the "
                "max unit cell length");
    throw std::runtime_error(
        "Negative PeakQspans are not allowed in IntegratePeakTimeSlices");
  }
205

206
207
208
209
210
211
  MatrixWorkspace_const_sptr inpWkSpace = getProperty("InputWorkspace");
  if (!inpWkSpace) {
    g_log.error("Improper Input Workspace");
    throw std::runtime_error(
        "Improper Input Workspace in IntegratePeakTimeSlices");
  }
212

213
214
215
216
217
218
  PeaksWorkspace_sptr peaksW;
  peaksW = getProperty("Peaks");
  if (!peaksW) {
    g_log.error("Improper Peaks Input");
    throw std::runtime_error("Improper Peaks Input");
  }
219

220
  int indx = getProperty("PeakIndex");
221

222
  IPeak &peak = peaksW->getPeak(indx);
223

224
225
226
227
  //------------------------------- Get Panel
  //--------------------------------------
  boost::shared_ptr<const Geometry::IComponent> panel_const =
      peak.getInstrument()->getComponentByName(peak.getBankName());
228

229
230
  boost::shared_ptr<Geometry::IComponent> panel =
      boost::const_pointer_cast<Geometry::IComponent>(panel_const);
231

232
233
234
235
  if (!panel || !panel_const) {
    g_log.error("Cannot get panel for a peak");
    throw std::runtime_error("Cannot get panel for a peak");
  }
236

237
238
  BoundingBox box;
  panel->getBoundingBox(box);
239

240
241
242
243
244
  int detID = peak.getDetectorID();
  if (!box.isPointInside(peak.getDetPos())) {
    g_log.error("Detector pixel is NOT inside the Peaks Bank");
    throw std::runtime_error("Detector pixel is NOT inside the Peaks Bank");
  }
245

246
247
  FindPlane(m_center, m_xvec, m_yvec, m_ROW, m_COL, m_NROWS, m_NCOLS,
            m_cellWidth, m_cellHeight, peak);
248

249
  g_log.debug() << "   Peak Index " << indx << std::endl;
250

251
252
  double TotVariance = 0;
  double TotIntensity = 0;
253
  double lastRow = m_ROW;
254
  double Row0 = lastRow;
255
  double lastCol = m_COL;
256
257
  double Col0 = lastCol;
  string spec_idList = "";
258

259
  // For quickly looking up workspace index from det id
260
  m_wi_to_detid_map = inpWkSpace->getDetectorIDToWorkspaceIndexMap();
261

Peterson, Peter's avatar
Peterson, Peter committed
262
  TableWorkspace_sptr TabWS = boost::make_shared<TableWorkspace>(0);
263

264
265
266
  //----------------------------- get Peak extents
  //------------------------------
  try {
267

268
    // Find the workspace index for this detector ID
269
    detid2index_map::const_iterator it = m_wi_to_detid_map.find(detID);
270
    size_t wsIndx = (it->second);
271

272
    double R = CalculatePositionSpan(peak, dQ) / 2;
273

274
275
    R = min<double>(MaxRowColSpan * max<double>(m_cellWidth, m_cellHeight), R);
    R = max<double>(MinRowColSpan * max<double>(m_cellWidth, m_cellHeight), R);
276

277
278
    R = 2 * R; // Gets a few more background cells.
    int Chan;
279

280
281
    const MantidVec &X = inpWkSpace->dataX(wsIndx);
    int dChan = CalculateTimeChannelSpan(peak, dQ, X, int(wsIndx), Chan);
282

283
    dChan = max<int>(dChan, MinTimeSpan);
284

285
286
287
    double Centy = Row0;
    double Centx = Col0;
    IDetector_const_sptr CenterDet = peak.getDetector();
288

289
    double neighborRadius; // last radius for finding neighbors
290

291
292
    neighborRadius =
        min<double>(MaxNeighborhoodRadius, NeighborhoodRadiusDivPeakRadius * R);
Hahn, Steven's avatar
Hahn, Steven committed
293
294
    int Nneighbors = static_cast<int>(neighborRadius * neighborRadius /
                                      m_cellWidth / m_cellHeight * 4);
295

Hahn, Steven's avatar
Hahn, Steven committed
296
297
    Nneighbors = min<int>(
        Nneighbors, static_cast<int>(inpWkSpace->getNumberHistograms()) - 2);
298
    delete[] m_NeighborIDs;
299

300
301
302
303
304
    m_NeighborIDs = new int[Nneighbors + 2];
    m_NeighborIDs[0] = Nneighbors + 2;
    m_NeighborIDs[1] = 2;
    Kernel::V3D Cent =
        (m_center + m_xvec * (Centx - m_COL) + m_yvec * (Centy - m_ROW));
305

306
    getNeighborPixIDs(panel, Cent, neighborRadius, m_NeighborIDs);
307

308
    if (m_NeighborIDs[1] < 10) {
309
310
311
312
313
314
      g_log.error("Not enough neighboring pixels to fit ");
      throw std::runtime_error("Not enough neighboring pixels to fit ");
    }
    int NBadEdgeCells = getProperty("NBadEdgePixels");
    int MaxChan = -1;
    double MaxCounts = -1;
315

316
317
318
319
320
321
322
    //     --------------- Find Time Chan with max counts----------------
    for (int dir = 1; dir > -2; dir -= 2) {
      bool done = false;
      for (int t = 0; t < dChan && !done; t++)
        if (dir < 0 && t == 0) {
          Centy = Row0;
          Centx = Col0;
Hahn, Steven's avatar
Hahn, Steven committed
323
324
        } else if (Chan + dir * t < 0 ||
                   Chan + dir * t >= static_cast<int>(X.size()))
325
326
          done = true;
        else {
327

328
          int NN = m_NeighborIDs[1];
329
330
          MatrixWorkspace_sptr Data = WorkspaceFactory::Instance().create(
              std::string("Workspace2D"), 3, NN, NN);
331

Peterson, Peter's avatar
Peterson, Peter committed
332
          auto XXX = boost::make_shared<DataModeHandler>(
333
              R, R, Centy, Centx, m_cellWidth, m_cellHeight,
334
              getProperty("CalculateVariances"), NBadEdgeCells,
Peterson, Peter's avatar
Peterson, Peter committed
335
              m_NCOLS - NBadEdgeCells, NBadEdgeCells, m_NROWS - NBadEdgeCells);
336
          m_AttributeValues = XXX;
337
          XXX->setCurrentRadius(R);
338

339
340
          SetUpData1(Data, inpWkSpace, Chan + dir * t, Chan + dir * t, R,
                     CenterDet->getPos(), spec_idList);
341

342
343
344
          if (m_AttributeValues->StatBaseVals(ISSIxx) > 0) {
            if (m_AttributeValues->StatBaseVals(IIntensities) > MaxCounts) {
              MaxCounts = m_AttributeValues->StatBaseVals(IIntensities);
345
              MaxChan = Chan + dir * t;
346
            }
347
348
349
350
351
            if (m_AttributeValues->StatBaseVals(IIntensities) > 0) {
              Centx = m_AttributeValues->StatBaseVals(ISSIx) /
                      m_AttributeValues->StatBaseVals(IIntensities);
              Centy = m_AttributeValues->StatBaseVals(ISSIy) /
                      m_AttributeValues->StatBaseVals(IIntensities);
352
            } else
353
              done = true;
354
355
          } else
            done = true;
356

357
358
          if (t >= 3 && (m_AttributeValues->StatBaseVals(IIntensities) <
                         MaxCounts / 2.0) &&
359
360
              MaxCounts >= 0)
            done = true;
361
        }
362
363
364
    }
    if (MaxChan > 0)
      Chan = MaxChan;
365

366
367
368
    g_log.debug() << "   largest Channel,Radius,m_cellWidth,m_cellHeight = "
                  << Chan << " " << R << " " << m_cellWidth << " "
                  << m_cellHeight << std::endl;
369

370
    if (R < MinRowColSpan / 2 * max<double>(m_cellWidth, m_cellHeight) ||
371
372
373
374
        dChan < MinTimeSpan) {
      g_log.error("Not enough rows and cols or time channels ");
      throw std::runtime_error("Not enough rows and cols or time channels ");
    }
375

376
    InitializeColumnNamesInTableWorkspace(TabWS);
377

378
379
380
381
    //------------------------------------- Start the Integrating
    //-------------------------------
    double time;
    int ncells;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
382

383
    Mantid::API::Progress prog(this, 0.0, 100.0, dChan);
384

385
386
    // Set from attributes replace by m_R0
    m_R0 = -1;
387
    int LastTableRow = -1;
Peterson, Peter's avatar
Peterson, Peter committed
388
389
    auto origAttributeList = boost::make_shared<DataModeHandler>();
    auto lastAttributeList = boost::make_shared<DataModeHandler>();
390

391
392
    for (int dir = 1; dir >= -1; dir -= 2) {
      bool done = false;
393

394
395
396
397
398
399
400
      for (int chan = 0; chan < dChan && !done; chan++)
        if (dir < 0 && chan == 0) {
          lastRow = Row0;
          lastCol = Col0;
          lastAttributeList = origAttributeList;
          if (TabWS->rowCount() > 0)
            LastTableRow = 0;
401

Hahn, Steven's avatar
Hahn, Steven committed
402
403
        } else if (Chan + dir * chan < 0 ||
                   Chan + dir * chan >= static_cast<int>(X.size()))
404
405
          done = true;
        else {
406

407
          int xchan = Chan + dir * chan;
408

409
410
411
          size_t topIndex = xchan + 1;
          if (topIndex >= X.size())
            topIndex = X.size() - 1;
412

413
          time = (X[xchan] + X[topIndex]) / 2.0;
414

415
          double Radius = R;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
416

417
418
          if (m_R0 > 0)
            Radius = m_R0;
419

420
          int NN = m_NeighborIDs[1];
421

422
423
424
          MatrixWorkspace_sptr Data = WorkspaceFactory::Instance().create(
              std::string("Workspace2D"), 3, NN,
              NN); // neighbors.size(), neighbors.size());
425

426
427
428
          g_log.debug() << " A:chan=" << xchan << "  time=" << time
                        << "   Radius=" << Radius << "row= " << lastRow
                        << "  col=" << lastCol << std::endl;
429

430
431
          SetUpData(Data, inpWkSpace, panel, xchan, xchan, lastCol, lastRow,
                    Cent, neighborRadius, Radius, spec_idList);
432

433
          m_AttributeValues->setTime(time);
434

435
          // if( dir==1 && chan ==0)
436
          //    origAttributeList= m_AttributeValues;
437

438
          ncells = static_cast<int>(m_AttributeValues->StatBaseVals(ISS1));
439

440
441
442
          std::vector<double> params;
          std::vector<double> errs;
          std::vector<std::string> names;
Roman Tolchenov's avatar
Roman Tolchenov committed
443

444
          if (m_AttributeValues->StatBaseVals(ISSIxx) > 0 &&
Hahn, Steven's avatar
Hahn, Steven committed
445
446
              m_AttributeValues->IsEnoughData(m_ParameterValues.data(),
                                              g_log) &&
447
              m_ParameterValues[ITINTENS] > 0) {
448
            double chisqOverDOF;
449

450
451
            Fit(Data, chisqOverDOF, done, names, params, errs, lastRow, lastCol,
                neighborRadius);
452

453
454
            if (!done) // Bivariate error happened
            {
455

456
457
458
              if (isGoodFit(params, errs, names, chisqOverDOF)) {
                LastTableRow = UpdateOutputWS(
                    TabWS, dir, xchan, params, errs, names, chisqOverDOF,
459
                    m_AttributeValues->time, spec_idList);
460

461
                double TotSliceIntensity =
462
                    m_AttributeValues->StatBaseVals(IIntensities);
463
                double TotSliceVariance =
464
                    m_AttributeValues->StatBaseVals(IVariance);
465

466
467
468
                updatePeakInformation(params, errs, names, TotVariance,
                                      TotIntensity, TotSliceIntensity,
                                      TotSliceVariance, chisqOverDOF, ncells);
469

470
                lastAttributeList = m_AttributeValues;
471

472
473
474
                if (dir == 1 && chan == 0)
                  origAttributeList = lastAttributeList;
              } else
475

476
477
                done = true;
            }
478

479
480
481
482
          } else //(!IsEnoughData() || ParameterValues[ITINTENS] <= 0
          {
            done = true;
          }
483

484
485
486
487
488
          if (done) // try to merge
          {
            done = false;

            int chanMin, chanMax;
489
            if ((dir == 1 && chan == 0) || lastAttributeList->CellHeight <= 0) {
490
491
492
493
              chanMin = xchan;
              chanMax = xchan + 1;
              if (dir < 0)
                chanMax++;
Peterson, Peter's avatar
Peterson, Peter committed
494
495
              auto XXX =
                  boost::make_shared<DataModeHandler>(*m_AttributeValues);
496
              m_AttributeValues = XXX;
497
              if (!X.empty())
498
                m_AttributeValues->setTime((X[chanMax] + X[chanMin]) / 2.0);
499

500
            } else // lastAttributeList exists
501

502
503
504
505
506
            {
              chanMin = std::min<int>(xchan, xchan - dir);
              chanMax = chanMin + 1;
              if (lastAttributeList->case4)
                chanMax++;
507

Peterson, Peter's avatar
Peterson, Peter committed
508
509
              auto XXX =
                  boost::make_shared<DataModeHandler>(*lastAttributeList);
510
              m_AttributeValues = XXX;
511

512
513
              m_AttributeValues->setTime((time + m_AttributeValues->time) /
                                         2.0);
514
            }
515

516
517
            if (updateNeighbors(panel, m_AttributeValues->getCurrentCenter(),
                                Cent, m_AttributeValues->getCurrentRadius(),
518
                                neighborRadius))
519
              Cent = m_AttributeValues->getCurrentCenter();
520

521
            int NN = m_NeighborIDs[1];
522
523
            MatrixWorkspace_sptr Data = WorkspaceFactory::Instance().create(
                std::string("Workspace2D"), 3, NN, NN);
524

525
            SetUpData1(Data, inpWkSpace, chanMin, chanMax,
526
527
                       m_AttributeValues->getCurrentRadius(),
                       m_AttributeValues->getCurrentCenter(), spec_idList);
528

529
            double chisqOverDOF;
530

531
            g_log.debug("Try Merge 2 time slices");
532
            if (m_AttributeValues->StatBaseVals(ISSIxx) >= 0 &&
Hahn, Steven's avatar
Hahn, Steven committed
533
534
                m_AttributeValues->IsEnoughData(m_ParameterValues.data(),
                                                g_log))
535

536
537
538
539
              Fit(Data, chisqOverDOF, done, names, params, errs, lastRow,
                  lastCol, neighborRadius);
            else
              chisqOverDOF = -1;
540

541
            if (!done && isGoodFit(params, errs, names, chisqOverDOF)) {
542

Hahn, Steven's avatar
Hahn, Steven committed
543
544
              if (LastTableRow >= 0 &&
                  LastTableRow < static_cast<int>(TabWS->rowCount()))
545
546
547
                TabWS->removeRow(LastTableRow);
              else
                LastTableRow = -1;
548

549
550
              LastTableRow = UpdateOutputWS(
                  TabWS, dir, (chanMin + chanMax) / 2.0, params, errs, names,
551
                  chisqOverDOF, m_AttributeValues->time, spec_idList);
552

553
554
555
556
              if (lastAttributeList->lastISAWVariance > 0 &&
                  lastAttributeList->CellHeight > 0) {
                TotIntensity -= lastAttributeList->lastISAWIntensity;
                TotVariance -= lastAttributeList->lastISAWVariance;
557
558
              }

559
              double TotSliceIntensity =
560
                  m_AttributeValues->StatBaseVals(IIntensities);
561

562
              double TotSliceVariance =
563
                  m_AttributeValues->StatBaseVals(IVariance);
564

Hahn, Steven's avatar
Hahn, Steven committed
565
566
567
568
              updatePeakInformation(
                  params, errs, names, TotVariance, TotIntensity,
                  TotSliceIntensity, TotSliceVariance, chisqOverDOF,
                  static_cast<int>(m_AttributeValues->StatBaseVals(ISS1)));
569

570
              // lastAttributeList= m_AttributeValues;
571

572
              if (dir == 1 && (chan == 0 || chan == 1)) {
573
574
                m_AttributeValues->case4 = true;
                origAttributeList = m_AttributeValues;
575
576
              } else
                LastTableRow = -1;
577

578
            } else {
Peterson, Peter's avatar
Peterson, Peter committed
579
              auto XXX = boost::make_shared<DataModeHandler>();
580
581
582
583
              lastAttributeList = XXX;
            }
            done = true;
          }
584

585
586
          // Get ready for the next round
          Data.reset();
587

588
          if (!done) {
589

590
591
            // Now set up the center for this peak
            int i = find("Mrow", names);
592
593
594
595
596
597
            if (i < 0) {
              throw std::runtime_error("Inconsistency found in algorithm "
                                       "execution. The index for the parameter "
                                       "Mrow is negative.");
            }

598
            lastRow = static_cast<int>(params[i] + .5);
599
            i = find("Mcol", names);
600
            if (i >= 0)
601
              lastCol = static_cast<int>(params[i] + .5);
602
            prog.report();
603

604
605
606
607
          } else if (dir > 0)
            prog.report(dChan / 2);
          else
            prog.report(dChan);
608

609
610
611
          params.clear();
          errs.clear();
          names.clear();
612
        }
613
    }
614

615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
  } catch (std::exception &EE1) {
    std::cout << "Error in main reason=" << EE1.what() << std::endl;

    throw std::runtime_error(" Error IntegratePeakTimeSlices:" +
                             std::string(EE1.what()));
  } catch (std::string &mess) {
    throw std::runtime_error("Error IntegratePeakTimeSlices:" + mess);

  } catch (...) {
    throw std::runtime_error("Error IntegratePeakTimeSlices:");
  }

  try {

    setProperty("OutputWorkspace", TabWS);

    setProperty("Intensity", TotIntensity);
    setProperty("SigmaIntensity", SQRT(TotVariance));
    time_t seconds2;

635
    seconds2 = time(nullptr);
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
    double dif = difftime(seconds2, seconds1);
    g_log.debug() << "Finished Integr peak number " << indx << " in " << dif
                  << " seconds" << std::endl;

  } catch (std::exception &ss) {

    std::cout << "Error occurred XX " << ss.what() << std::endl;
    throw std::runtime_error(ss.what());
  }
}

/**
   * Finds all neighbors within a given Radius of the Center on the given
 * component.
   * @param comp -The component of interest
   * @param Center- the center of the neighbors
   * @param Radius - The radius from the center of neighbors to be included
   * @param ArryofID -The detector ID's of the neighbors. The id of the pixel at
 * the center may be included.
   */
bool IntegratePeakTimeSlices::getNeighborPixIDs(
    boost::shared_ptr<Geometry::IComponent> comp, Kernel::V3D &Center,
    double &Radius, int *&ArryofID) {

  int N = ArryofID[1];
  int MaxN = ArryofID[0];

  if (N >= MaxN)
    return false;

  Geometry::BoundingBox box;
  comp->getBoundingBox(box);

  double minx = Center.X() - Radius;
  double miny = Center.Y() - Radius;
  double minz = Center.Z() - Radius;
  double maxx = Center.X() + Radius;
  double maxy = Center.Y() + Radius;
  double maxz = Center.Z() + Radius;

  if (box.xMin() >= maxx)
    return true;

  if (box.xMax() <= minx)
    return true;
  ;

  if (box.yMin() >= maxy)
    return true;

  if (box.yMax() <= miny)
    return true;

  if (box.zMin() >= maxz)
    return true;

  if (box.zMax() <= minz)
    return true;
  ;

  boost::shared_ptr<Geometry::Detector> det =
      boost::dynamic_pointer_cast<Geometry::Detector>(comp);
  // if( comp->type().compare(0,8,"Detector")==0 ||
  // comp->type().compare("RectangularDetectorPixel")==0)
  if (det) {
    V3D pos = det->getPos() - Center;
    if (pos.X() * pos.X() + pos.Y() * pos.Y() + pos.Z() * pos.Z() <
        Radius * Radius) {
      ArryofID[N] = det->getID();
      N++;
      ArryofID[1] = N;
    }
    return true;
  }
710

711
712
  boost::shared_ptr<const Geometry::ICompAssembly> Assembly =
      boost::dynamic_pointer_cast<const Geometry::ICompAssembly>(comp);
713

714
715
  if (!Assembly)
    return true;
716

717
  bool b = true;
718

719
720
  for (int i = 0; i < Assembly->nelements() && b; i++)
    b = getNeighborPixIDs(Assembly->getChild(i), Center, Radius, ArryofID);
721

722
723
  return b;
}
724

725
/**
726
 * Checks and updates if needed the list of m_NeighborIDs
727
728
729
730
731
732
733
734
735
736
737
738
739
 * @param comp  Component with the neighboring pixels
 * @param CentPos    new Center
 * @param oldCenter   old Center
 * @param NewRadius   new Radius
 * @param neighborRadius  old the new neighborhood radius
 */
bool IntegratePeakTimeSlices::updateNeighbors(
    boost::shared_ptr<Geometry::IComponent> &comp, V3D CentPos, V3D oldCenter,
    double NewRadius, double &neighborRadius) {
  double DD = (CentPos - oldCenter).norm();
  bool changed = false;
  if (DD + NewRadius > neighborRadius) {
    int NN = int(NStdDevPeakSpan * NeighborhoodRadiusDivPeakRadius * NewRadius /
740
741
742
743
744
745
                 m_cellWidth * NStdDevPeakSpan *
                 NeighborhoodRadiusDivPeakRadius * NewRadius / m_cellHeight);
    if (m_NeighborIDs[0] < NN) {
      delete[] m_NeighborIDs;
      m_NeighborIDs = new int[NN + 2];
      m_NeighborIDs[0] = NN + 2;
746
    }
747
    m_NeighborIDs[1] = 2;
748
    neighborRadius = NeighborhoodRadiusDivPeakRadius * NewRadius;
749

750
    getNeighborPixIDs(comp, CentPos, neighborRadius, m_NeighborIDs);
751
    changed = true;
752

753
754
  } else // big enough neighborhood so
    neighborRadius -= DD;
755

756
757
  return changed;
}
758

759
760
761
762
763
764
765
766
767
768
/**
 * Calculates the span in rows and columns needed to include all data within dQ
 *of the
 * specified peak
 * @param peak  The peak of interest
 * @param dQ    The offset from the peak's Q value for the data of interest
 *
 * NOTE: differentials of Q =mv*sin(scatAng/2)/2 were used to calculate this
 *  Also s=r*theta was used to transfer d ScatAng to distance on a bank.
 */
769
770
771
double
IntegratePeakTimeSlices::CalculatePositionSpan(Geometry::IPeak const &peak,
                                               const double dQ) {
772

773
774
  try {
    double Q = 0, ScatAngle = 0, dScatAngle = 0, DetSpan = 0;
775

776
777
778
779
    Q = peak.getQLabFrame().norm();
    Geometry::Instrument_const_sptr instr = peak.getInstrument();
    const Geometry::IComponent_const_sptr sample = instr->getSample();
    V3D pos = peak.getDetPos() - sample->getPos();
780

781
    ScatAngle = acos(pos.Z() / pos.norm());
782

783
    dScatAngle = 2 * dQ / Q * tan(ScatAngle / 2);
784

785
    DetSpan = pos.norm() * dScatAngle; // s=r*theta
786

787
    DetSpan = fabs(DetSpan);
788

789
    // IDetector_sptr det = peak.getDetector();
790

791
    return DetSpan;
792

793
794
795
796
797
  } catch (std::exception &s) {
    std::cout << "err in getNRowsCols, reason=" << s.what() << std::endl;
    return 0;
  }
}
798

799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
/**
 * Calculates the span of channels needed to encompass all data around the peak
 *with Q values within dQ
 * of this peak's Q value
 *
 * @param peak     The peak of interest
 * @param dQ       The offset of peak's Q value whose data is considered part of
 *the peak
 * @param X        The list of time channel values.
 * @param specNum  The spectral number for the pixel(Not Currently Used)
 * @param Centerchan  The center time channel number( from X)
 *
 * @return The number of time channels around Centerchan to use
 */
int IntegratePeakTimeSlices::CalculateTimeChannelSpan(
Owen Arnold's avatar
Owen Arnold committed
814
    Geometry::IPeak const &peak, const double dQ, Mantid::MantidVec const &X,
815
816
817
    const int specNum, int &Centerchan) {
  UNUSED_ARG(specNum);
  double Q = peak.getQLabFrame().norm(); // getQ( peak)/2/M_PI;
818

819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
  double time = peak.getTOF();
  double dtime = dQ / Q * time;
  int chanCenter = find(X, time);

  Centerchan = chanCenter;
  int chanLeft = find(X, time - dtime);
  int chanRight = find(X, time + dtime);
  int dchan = abs(chanCenter - chanLeft);

  if (abs(chanRight - chanCenter) > dchan)
    dchan = abs(chanRight - chanCenter);

  dchan = max<int>(3, dchan);

  return dchan + 5; // heuristic should be a lot more
}

/**
 * For NonFlat banks, this determines the data of a small planar region
 *approximating the instrument
 * close to the peak
 *
 * @param center  The position of the center of this plane
 * @param xvec    The direction the column values increase
 * @param yvec    The direction the row values increase
 * @param ROW     The row for this peak( 0 if undefined)
 * @param COL     The col for this peak( 0 if undefined)
 * @param NROWS   The number? of rows for this bANK
 * @param NCOLS    The number of columns for this bank
 * @param pixWidthx   The width of a pixel
 * @param pixHeighty  The height of a pixel
 * @param peak
 */
void IntegratePeakTimeSlices::FindPlane(V3D &center, V3D &xvec, V3D &yvec,
                                        double &ROW, double &COL, int &NROWS,
                                        int &NCOLS, double &pixWidthx,
                                        double &pixHeighty,
Owen Arnold's avatar
Owen Arnold committed
856
                                        Geometry::IPeak const &peak) const {
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931

  NROWS = NCOLS = -1;
  IDetector_const_sptr det = peak.getDetector();
  V3D detPos = det->getPos();

  center.setX(detPos.X());
  center.setY(detPos.Y());
  center.setZ(detPos.Z());

  boost::shared_ptr<const Detector> dett =
      boost::dynamic_pointer_cast<const Detector>(det);

  pixWidthx = dett->getWidth();
  pixHeighty = dett->getHeight();

  Kernel::Quat Qt = dett->getRotation();
  V3D yaxis0(0.0, 1.0, 0.0);

  Qt.rotate(yaxis0);
  yaxis0.normalize();

  V3D xaxis0(1, 0, 0);
  Qt.rotate(xaxis0);
  xaxis0.normalize();

  xvec.setX(xaxis0.X());
  xvec.setY(xaxis0.Y());
  xvec.setZ(xaxis0.Z());
  yvec.setX(yaxis0.X());
  yvec.setY(yaxis0.Y());
  yvec.setZ(yaxis0.Z());
  ROW = peak.getRow();
  COL = peak.getCol();
  Geometry::Instrument_const_sptr inst = peak.getInstrument();
  if (!inst)
    throw std::invalid_argument("No instrument for peak");
  boost::shared_ptr<const IComponent> panel =
      inst->getComponentByName(peak.getBankName());

  boost::shared_ptr<const RectangularDetector> ddet =
      boost::dynamic_pointer_cast<const RectangularDetector>(panel);

  if (ddet) {
    std::pair<int, int> CR = ddet->getXYForDetectorID(det->getID());
    ROW = CR.second;
    COL = CR.first;
    pixWidthx = ddet->xstep();
    pixHeighty = ddet->ystep();

    NROWS = ddet->ypixels();
    NCOLS = ddet->xpixels();

    return;
  }
  // Get NROWS and NCOLS for other panels
  NROWS = NCOLS = -1;

  if (!inst)
    return;

  if (!panel)
    return;
  boost::shared_ptr<const Component> compPanel =
      boost::dynamic_pointer_cast<const Component>(panel);
  boost::shared_ptr<IComponent> panel1(compPanel->base()->clone());
  BoundingBox B;

  Quat rot = panel1->getRotation();

  rot.inverse();

  panel1->rotate(rot);

  panel1->getBoundingBox(B);

932
933
  NROWS = static_cast<int>((B.yMax() - B.yMin()) / pixHeighty + .5);
  NCOLS = static_cast<int>((B.xMax() - B.xMin()) / pixWidthx + .5);
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
}

/**
 * Updates the cumulative statistics for the data being considered
 *
 * @param intensity -The experimental intensity at the given pixel
 * @param variance  -The square of the errors in the above intensity
 * @param row       -The row of the given pixel
 * @param col       -The column of the given pixel
 * @param StatBase  -The data accumulator
 */
void IntegratePeakTimeSlices::updateStats(const double intensity,
                                          const double variance,
                                          const double row, const double col,
                                          std::vector<double> &StatBase) {

  StatBase[ISSIxx] += col * col * intensity;
  StatBase[ISSIyy] += intensity * row * row;
  StatBase[ISSIxy] += intensity * row * col;
  StatBase[ISSxx] += col * col;
  StatBase[ISSyy] += row * row;
  StatBase[ISSxy] += row * col;
  StatBase[ISSIx] += intensity * col;
  StatBase[ISSIy] += intensity * row;
  StatBase[ISSx] += col;
  StatBase[ISSy] += row;
  StatBase[IIntensities] += intensity;
  StatBase[IVariance] += variance;
  StatBase[ISS1] += 1;
}

/**
 * Finds and saves the initial values and modes( like isEdge) for this data.
 *
 * @param Varx  The Variance of the x(column) values
 * @param Vary  The Variance of the y(row)  values
 * @param b     The background
 *
 * NOTE:Used only in GetParams which is used only in PreFit which is not
 * currently used.
 */
std::vector<double> DataModeHandler::InitValues(double Varx, double Vary,
                                                double b) {
  std::vector<double> Res(7);

  Res[IVXX] = Varx;
  Res[IVYY] = Vary;
  Res[IVXY] = 0;
982
  int nCells = static_cast<int>(StatBase[ISS1]);
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
  double Den = StatBase[IIntensities] - b * nCells;
  Res[IXMEAN] = (StatBase[ISSIx] - b * StatBase[ISSx]) / Den;
  Res[IYMEAN] = (StatBase[ISSIy] - b * StatBase[ISSy]) / Den;
  Res[IBACK] = b;
  Res[ITINTENS] = StatBase[IIntensities] - b * nCells;

  //---- Is Edge Cell ???-------
  double NstdX = 4 * (currentRadius / CellWidth - EdgeX) / sqrt(Varx);
  double NstdY = 4 * (currentRadius / CellHeight - EdgeY) / sqrt(Vary);
  double sigx = 1;
  double sigy = 1;
  if (NstdX < 0)
    sigx = -1;
  if (NstdY < 0)
    sigy = -1;

  double x = 1;
  if (sigy * NstdY < 7 && sigy * NstdY >= 0) // is close to row edge
For faster browsing, not all history is shown. View entire blame