IntegratePeakTimeSlices.h

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2009 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
/*
 * IntegratePeakTimeSlices.h
 *
 *  Created on: May 5, 2011
 *      Author: ruth
 */
#pragma once
 
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/IAlgorithm.h"
#include "MantidAPI/ITableWorkspace_fwd.h"
#include "MantidAPI/MatrixWorkspace_fwd.h"
#include "MantidAPI/SpectraDetectorTypes.h"
#include "MantidCrystal/DllConfig.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidGeometry/Crystal/IPeak.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidKernel/V3D.h"
#include "MantidKernel/cow_ptr.h"
 
#include <array>
 
// Forward Declarations
namespace Mantid {
namespace HistogramData {
class HistogramX;
}
} // namespace Mantid
namespace Mantid {
namespace Crystal {
class DataModeHandler {
public:
  DataModeHandler() { init(); }
 
  DataModeHandler(const DataModeHandler &handler);
  DataModeHandler(double baseRCRadius, double lastRCRadius, double lastRow, double lastCol, double CellWidth,
                  double CellHeight, bool CalcVariance, int MinCol, int MaxCol, int MinRow, int MaxRow) {
    init();
    this->baseRCRadius = baseRCRadius;
    this->lastRCRadius = lastRCRadius;
    this->lastRow = lastRow;
    this->lastCol = lastCol;
    this->CellWidth = CellWidth;
    this->CellHeight = CellHeight;
    this->CalcVariance = CalcVariance;
    this->MaxCol = MaxCol;
    this->MaxRow = MaxRow;
    this->MinCol = MinCol;
    this->MinRow = MinRow;
  }
 
  void setTime(double time) { this->time = time; }
 
  // Returns true if "Edge Peak" otherwise returns false
  bool setStatBase(std::vector<double> const &StatBase);
 
  bool isEdgePeak(const double *params, int nparams);
 
  void setHeightHalfWidthInfo(const MantidVec &xvals, const MantidVec &yvals, const MantidVec &counts);
 
  void setCurrentRadius(double radius) { currentRadius = radius; }
  void setCurrentCenter(const Kernel::V3D &newCenter) {
    Kernel::V3D XX(newCenter);
    currentPosition = XX;
  }
 
  double getCurrentRadius() { return currentRadius; }
  Kernel::V3D getCurrentCenter() { return currentPosition; }
  void updateEdgeXsize(double newsize) {
    if (EdgeX < 0)
 
      EdgeX = newsize;
 
    else if (newsize < EdgeX)
 
      EdgeX = newsize;
  }
 
  void updateEdgeYsize(double newsize) {
    if (EdgeY < 0)
 
      EdgeY = newsize;
 
    else if (newsize < EdgeY)
 
      EdgeY = newsize;
  }
 
  void CalcVariancesFromData(double background, double meanx, double meany, double &Varxx, double &Varxy, double &Varyy,
                             const std::vector<double> &StatBase);
 
  bool IsEnoughData(const double *ParameterValues, Kernel::Logger &);
 
  double getNewRCRadius();
 
  double getInitBackground() { return back_calc; }
 
  double getInitRow() { return row_calc; }
 
  double getInitCol() { return col_calc; }
 
  double getInitIntensity() { return Intensity_calc; }
 
  double getInitVarx() { return Vx_calc; }
 
  double getInitVary() { return Vy_calc; }
 
  double getInitVarxy() { return Vxy_calc; }
 
  std::string CalcConstraints(std::vector<std::pair<double, double>> &Bounds, bool CalcVariances);
  std::string getTies() { return ""; }
 
  bool CalcVariances();
 
  std::vector<double> GetParams(double b);
 
  double StatBaseVals(int index) {
    if (index < 0 || index >= (int)StatBase.size())
      return 0.0;
    return StatBase[index];
  }
 
  double CalcISAWIntensity(const double *params);
 
  double CalcISAWIntensityVariance(const double *params, const double *errs, double chiSqOvDOF);
 
  double CalcSampleIntensityMultiplier(const double *params) const;
 
  std::vector<double> InitValues(double Varx, double Vary, double b);
  double baseRCRadius;
  double lastRCRadius;
  double HalfWidthAtHalfHeightRadius;
  double calcNewRCRadius;
 
  double lastRow;
  int MinRow, MaxRow;
  double lastCol;
  int MinCol, MaxCol;
  double time;
  double CellWidth;
  double CellHeight;
 
  double VarxHW, VaryHW;
  double currentRadius;
  Kernel::V3D currentPosition;
  std::vector<double> StatBase;
 
  double EdgeX, EdgeY;
  double lastISAWIntensity, lastISAWVariance;
  bool CalcVariance;
  bool case4; // if true result of successful merge of dir =1 chan=0 and chan=1
  double back_calc, Intensity_calc, row_calc, col_calc, Vx_calc, Vy_calc, Vxy_calc;
 
private:
  void init() {
    baseRCRadius = -1;
    lastRCRadius = -1;
    lastRow = -1;
    lastCol = -1;
    EdgeX = EdgeY = -1;
    calcNewRCRadius = -1;
    MaxRow = -1;
    MaxCol = -1;
    MinRow = -1;
    MinCol = -1;
    time = -1;
    CalcVariance = true;
    CellWidth = CellHeight = 0;
    currentRadius = -1;
    lastISAWIntensity = -1;
    lastISAWVariance = -1;
    currentPosition = Kernel::V3D();
    HalfWidthAtHalfHeightRadius = -1;
    case4 = false;
 
    VarxHW = -1;
    VaryHW = -1;
    back_calc = Intensity_calc = row_calc = col_calc = Vx_calc = Vy_calc = Vxy_calc = -1;
  }
};
 
class MANTID_CRYSTAL_DLL IntegratePeakTimeSlices : public Mantid::API::Algorithm {
public:
  IntegratePeakTimeSlices();
 
  ~IntegratePeakTimeSlices() override;
 
  const std::string name() const override { return "IntegratePeakTimeSlices"; }
 
  const std::string summary() const override {
    return "The algorithm uses CurveFitting::BivariateNormal for fitting a "
           "time slice";
  }
 
  int version() const override { return 1; }
  const std::vector<std::string> seeAlso() const override { return {"PeakIntegration"}; }
 
  const std::string category() const override { return "Crystal\\Integration"; }
 
private:
  void init() override;
  void exec() override;
 
  bool m_EdgePeak;
 
  std::string m_AttributeNames[20];
 
  std::string m_ParameterNames[7];
 
  std::shared_ptr<DataModeHandler> m_AttributeValues;
  std::array<double, 7> m_ParameterValues;
 
  Mantid::detid2index_map m_wi_to_detid_map;
 
  int *m_NeighborIDs; // Stores IDs of nearest neighbors
  double m_R0;        
 
  Kernel::V3D m_center; 
  Kernel::V3D m_xvec;   
  Kernel::V3D m_yvec;   
  double m_ROW;         
  double m_COL;         
  double m_cellWidth;  
  double m_cellHeight; 
  int m_NROWS;
  int m_NCOLS;
 
  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, // from CentDet
                 double Radius, std::string &spec_idList);
 
  bool getNeighborPixIDs(const std::shared_ptr<Geometry::IComponent> &comp, const Kernel::V3D &Center, double &Radius,
                         int *&ArryofID);
 
  int CalculateTimeChannelSpan(Geometry::IPeak const &peak, const double dQ, const Mantid::HistogramData::HistogramX &X,
                               const int specNum, int &Centerchan);
 
  double CalculatePositionSpan(DataObjects::Peak const &peak, const double dQ);
 
  void InitializeColumnNamesInTableWorkspace(DataObjects::TableWorkspace_sptr &TabWS);
 
  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
 
  );
 
  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);
 
  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);
 
  std::string CalculateFunctionProperty_Fit();
 
  bool isGoodFit(std::vector<double> const &params, std::vector<double> const &errs,
                 std::vector<std::string> const &names, double chisqOverDOF);
  // returns last row added
  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);
 
  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);
 
  void updateStats(const double intensity, const double variance, const double row, const double col,
                   std::vector<double> &StatBase);
 
  int findNameInVector(std::string const &oneName, std::vector<std::string> const &nameList);
 
  double CalculateIsawIntegrateError(const double background, const double backError, const double ChiSqOverDOF,
                                     const double TotVariance, const int ncells);
 
  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;
 
  int findTimeChannel(const Mantid::HistogramData::HistogramX &X, const double time);
 
  // returns true if Neighborhood list is changed
  bool updateNeighbors(const std::shared_ptr<Geometry::IComponent> &comp, const Kernel::V3D &CentPos,
                       const Kernel::V3D &oldCenter, double NewRadius, double &neighborRadius);
};
} // namespace Crystal
} // namespace Mantid