WorkspaceCreationHelper.cpp

 * @param name :: property name
 * @param val :: value
 */
void addTSPEntry(Run &runInfo, std::string name, double val) {
  TimeSeriesProperty<double> *tsp;
  tsp = new TimeSeriesProperty<double>(name);
  tsp->addValue("2011-05-24T00:00:00", val);
  runInfo.addProperty(tsp);
}

// =====================================================================================
/** Sets the OrientedLattice in the crystal as an crystal with given lattice
 *lengths, angles of 90 deg
 *
 * @param ws :: workspace to set
 * @param a :: lattice length
 * @param b :: lattice length
 * @param c :: lattice length
 */
void setOrientedLattice(Mantid::API::MatrixWorkspace_sptr ws, double a,
                        double b, double c) {
  auto latt =
      Mantid::Kernel::make_unique<OrientedLattice>(a, b, c, 90., 90., 90.);
  ws->mutableSample().setOrientedLattice(latt.release());
}

// =====================================================================================
/** Create a default universal goniometer and set its angles
 *
 * @param ws :: workspace to set
 * @param phi :: +Y rotation angle (deg)
 * @param chi :: +X rotation angle (deg)
 * @param omega :: +Y rotation angle (deg)
 */
void setGoniometer(Mantid::API::MatrixWorkspace_sptr ws, double phi, double chi,
                   double omega) {
  addTSPEntry(ws->mutableRun(), "phi", phi);
  addTSPEntry(ws->mutableRun(), "chi", chi);
  addTSPEntry(ws->mutableRun(), "omega", omega);
  Mantid::Geometry::Goniometer gm;
  gm.makeUniversalGoniometer();
  ws->mutableRun().setGoniometer(gm, true);
}

//
Mantid::API::MatrixWorkspace_sptr
createProcessedWorkspaceWithCylComplexInstrument(size_t numPixels,
                                                 size_t numBins,
                                                 bool has_oriented_lattice) {
  size_t rHist = static_cast<size_t>(std::sqrt(static_cast<double>(numPixels)));
  while (rHist * rHist < numPixels)
    rHist++;

  Mantid::API::MatrixWorkspace_sptr ws =
      createGroupedWorkspace2DWithRingsAndBoxes(rHist, 10, 0.1);
  auto pAxis0 = new NumericAxis(numBins);
  for (size_t i = 0; i < numBins; i++) {
    double dE = -1.0 + static_cast<double>(i) * 0.8;
    pAxis0->setValue(i, dE);
  }
  pAxis0->setUnit("DeltaE");
  ws->replaceAxis(0, pAxis0);
  if (has_oriented_lattice) {
    auto latt =
        Mantid::Kernel::make_unique<OrientedLattice>(1, 1, 1, 90., 90., 90.);
    ws->mutableSample().setOrientedLattice(latt.release());

    addTSPEntry(ws->mutableRun(), "phi", 0);
    addTSPEntry(ws->mutableRun(), "chi", 0);
    addTSPEntry(ws->mutableRun(), "omega", 0);
    Mantid::Geometry::Goniometer gm;
    gm.makeUniversalGoniometer();
    ws->mutableRun().setGoniometer(gm, true);
  }

  return ws;
}

/// Create a workspace with all components needed for inelastic analysis and 3
/// detectors in specific places
/// @param L2        -- the sample to detector flight path
/// @param polar     -- the detector polar angle
/// @param azimutal  -- the detector azimuthal
/// @param numBins   -- the number of histogram bins for the workspace
/// @param Emin      -- minimal energy transfer
/// @param Emax      -- maxinal energy transfer
/// @param Ei        -- input beam energy
Mantid::API::MatrixWorkspace_sptr
createProcessedInelasticWS(const std::vector<double> &L2,
                           const std::vector<double> &polar,
                           const std::vector<double> &azimutal, size_t numBins,
                           double Emin, double Emax, double Ei) {
  // not used but interface needs it
  std::set<int64_t> maskedWorkspaceIndices;
  size_t numPixels = L2.size();

  Mantid::API::MatrixWorkspace_sptr ws =
      create2DWorkspaceWithValues(uint64_t(numPixels), uint64_t(numBins), true,
                                  maskedWorkspaceIndices, 0, 1, 0.1);

  // detectors at L2, sample at 0 and source at -L2_min
  ws->setInstrument(
      ComponentCreationHelper::createCylInstrumentWithDetInGivenPositions(
          L2, polar, azimutal));

  for (int g = 0; g < static_cast<int>(numPixels); g++) {
    auto &spec = ws->getSpectrum(g);
    // we just made (in createCylInstrumentWithDetInGivenPosisions) det ID-s to
    // start from 1
    spec.setDetectorID(g + 1);
    // and this is absolutely different nummer, corresponding to det ID just by
    // chance ? -- some uncertainties remain
    spec.setSpectrumNo(g + 1);
    // spec->setSpectrumNo(g+1);
    //   spec->addDetectorID(g*9);
    //   spec->setSpectrumNo(g+1); // Match detector ID and spec NO
  }

  const double dE = (Emax - Emin) / static_cast<double>(numBins);
  for (size_t j = 0; j < numPixels; j++) {
    std::vector<double> E_transfer;
    E_transfer.reserve(numBins);
    for (size_t i = 0; i <= numBins; i++) {
      E_transfer.push_back(Emin + static_cast<double>(i) * dE);
    }
    ws->mutableX(j) = E_transfer;
  }

  // set axis, correspondent to the X-values
  auto pAxis0 = new NumericAxis(numBins);
  const auto &E_transfer = ws->x(0);
  for (size_t i = 0; i < numBins; i++) {
    double E = 0.5 * (E_transfer[i] + E_transfer[i + 1]);
    pAxis0->setValue(i, E);
  }

  pAxis0->setUnit("DeltaE");

  ws->replaceAxis(0, pAxis0);

  // define oriented lattice which requested for processed ws
  auto latt =
      Mantid::Kernel::make_unique<OrientedLattice>(1, 1, 1, 90., 90., 90.);
  ws->mutableSample().setOrientedLattice(latt.release());

  ws->mutableRun().addProperty(
      new PropertyWithValue<std::string>("deltaE-mode", "Direct"), true);
  ws->mutableRun().addProperty(new PropertyWithValue<double>("Ei", Ei), true);
  // these properties have to be different -> specific for processed ws, as time
  // now should be reconciled
  addTSPEntry(ws->mutableRun(), "phi", 0);
  addTSPEntry(ws->mutableRun(), "chi", 0);
  addTSPEntry(ws->mutableRun(), "omega", 0);
  Mantid::Geometry::Goniometer gm;
  gm.makeUniversalGoniometer();
  ws->mutableRun().setGoniometer(gm, true);

  return ws;
}

/*
 * Create an EventWorkspace from a source EventWorkspace.
 * The new workspace should be exactly the same as the source workspace but
 * without any events
 */
Mantid::DataObjects::EventWorkspace_sptr
createEventWorkspace3(Mantid::DataObjects::EventWorkspace_const_sptr sourceWS,
                      std::string wsname, API::Algorithm *alg) {
  UNUSED_ARG(wsname);
  // 1. Initialize:use dummy numbers for arguments, for event workspace it
  // doesn't matter
  Mantid::DataObjects::EventWorkspace_sptr outputWS =
      Mantid::DataObjects::EventWorkspace_sptr(
          new DataObjects::EventWorkspace());
  outputWS->initialize(sourceWS->getInstrument()->getDetectorIDs(true).size(),
                       1, 1);

  // 2. Set the units
  outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
  outputWS->setYUnit("Counts");
  outputWS->setTitle("Empty_Title");

  // 3. Add the run_start property:
  int runnumber = sourceWS->getRunNumber();
  outputWS->mutableRun().addProperty("run_number", runnumber);

  std::string runstartstr = sourceWS->run().getProperty("run_start")->value();
  outputWS->mutableRun().addProperty("run_start", runstartstr);

  // 4. Instrument
  Mantid::API::Algorithm_sptr loadInst =
      alg->createChildAlgorithm("LoadInstrument");
  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  loadInst->setPropertyValue("InstrumentName",
                             sourceWS->getInstrument()->getName());
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", outputWS);
  loadInst->setProperty("RewriteSpectraMap",
                        Mantid::Kernel::OptionalBool(true));
  loadInst->executeAsChildAlg();
  // Populate the instrument parameters in this workspace - this works around a
  // bug
  outputWS->populateInstrumentParameters();

  // 6. Build spectrum and event list
  // a) We want to pad out empty pixels.
  detid2det_map detector_map;
  outputWS->getInstrument()->getDetectors(detector_map);

  // b) determine maximum pixel id
  detid2det_map::iterator it;
  detid_t detid_max = 0; // seems like a safe lower bound
  for (it = detector_map.begin(); it != detector_map.end(); ++it)
    if (it->first > detid_max)
      detid_max = it->first;

  // c) Pad all the pixels and Set to zero
  size_t workspaceIndex = 0;
  const auto &detectorInfo = outputWS->detectorInfo();
  for (it = detector_map.begin(); it != detector_map.end(); ++it) {
    if (!detectorInfo.isMonitor(detectorInfo.indexOf(it->first))) {
      auto &spec = outputWS->getSpectrum(workspaceIndex);
      spec.addDetectorID(it->first);
      // Start the spectrum number at 1
      spec.setSpectrumNo(specnum_t(workspaceIndex + 1));
      workspaceIndex += 1;
    }
  }

  return outputWS;
}

RebinnedOutput_sptr createRebinnedOutputWorkspace() {
  RebinnedOutput_sptr outputWS =
      Mantid::DataObjects::RebinnedOutput_sptr(new RebinnedOutput());
  // outputWS->setName("rebinTest");
  Mantid::API::AnalysisDataService::Instance().add("rebinTest", outputWS);

  // Set Q ('y') axis binning
  std::vector<double> qbins{0.0, 1.0, 4.0};
  std::vector<double> qaxis;
  const int numY =
      static_cast<int>(VectorHelper::createAxisFromRebinParams(qbins, qaxis));

  // Initialize the workspace
  const int numHist = numY - 1;
  const int numX = 7;
  outputWS->initialize(numHist, numX, numX - 1);

  // Set the normal units
  outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("DeltaE");
  outputWS->setYUnit("Counts");
  outputWS->setTitle("Empty_Title");

  // Create the i-axis for histogramming.
  HistogramData::BinEdges x1{-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0};

  // Create a numeric axis to replace the default vertical one
  Axis *const verticalAxis = new NumericAxis(numY);
  outputWS->replaceAxis(1, verticalAxis);

  // Now set the axis values
  for (int i = 0; i < numHist; ++i) {
    outputWS->setBinEdges(i, x1);
    verticalAxis->setValue(i, qaxis[i]);
  }
  // One more to set on the 'y' axis
  verticalAxis->setValue(numHist, qaxis[numHist]);

  // Set the 'y' axis units
  verticalAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
  verticalAxis->title() = "|Q|";

  // Set the X axis title (for conversion to MD)
  outputWS->getAxis(0)->title() = "Energy transfer";

  // Now, setup the data
  // Q bin #1

  // Populates destination starting at index 1 with the following data
  // e.g. y(0)[1] = 2.0, y(0)[2] = 3.0 ..etc. as the starting index is 1
  // if you change the values in the line below please update this comment!
  populateWsWithInitList(outputWS->mutableY(0), 1, {2.0, 3.0, 3.0, 2.0});
  populateWsWithInitList(outputWS->mutableE(0), 1, {2.0, 3.0, 3.0, 2.0});
  populateWsWithInitList(outputWS->dataF(0), 1, {2.0, 3.0, 3.0, 1.0});

  // Q bin #2
  populateWsWithInitList(outputWS->mutableY(1), 1, {1.0, 3.0, 3.0, 2.0, 2.0});
  populateWsWithInitList(outputWS->mutableE(1), 1, {1.0, 3.0, 3.0, 2.0, 2.0});
  populateWsWithInitList(outputWS->dataF(1), 1, {1.0, 3.0, 3.0, 1.0, 2.0});

  // Q bin #3
  populateWsWithInitList(outputWS->mutableY(2), 1, {1.0, 2.0, 3.0, 1.0});
  populateWsWithInitList(outputWS->mutableE(2), 1, {1.0, 2.0, 3.0, 1.0});
  populateWsWithInitList(outputWS->dataF(2), 1, {1.0, 2.0, 2.0, 1.0});

  // Q bin #4
  populateWsWithInitList(outputWS->mutableY(3), 0, {1.0, 2.0, 3.0, 2.0, 1.0});
  populateWsWithInitList(outputWS->mutableE(3), 0, {1.0, 2.0, 3.0, 2.0, 1.0});
  populateWsWithInitList(outputWS->dataF(3), 0, {1.0, 2.0, 3.0, 2.0, 1.0, 1.0});

  // Set representation
  outputWS->finalize();

  // Make errors squared rooted
  for (int i = 0; i < numHist; ++i) {
    auto &mutableE = outputWS->mutableE(i);
    for (int j = 0; j < numX - 1; ++j) {
      mutableE[j] = std::sqrt(mutableE[j]);
    }
  }

  return outputWS;
}

/**
  * Populates the destination array (usually a mutable histogram)
  * starting at the index specified with the doubles provided in an
  * initializer list. Note the caller is responsible for ensuring
  * the destination has capacity for startingIndex + size(initializer list)
  * number of values
  *
  * @param destination :: The array to populate with data
  * @param startingIndex :: The index to start populating data at
  * @param values :: The initializer list to populate the array with
  * starting at the index specified
  */
template <typename T>
void populateWsWithInitList(T &destination, size_t startingIndex,
                            const std::initializer_list<double> &values) {
  size_t index = 0;
  for (const double val : values) {
    destination[startingIndex + index] = val;
    index++;
  }
}

Mantid::DataObjects::PeaksWorkspace_sptr
createPeaksWorkspace(const int numPeaks, const bool createOrientedLattice) {
  auto peaksWS = boost::make_shared<PeaksWorkspace>();
  Instrument_sptr inst =
      ComponentCreationHelper::createTestInstrumentRectangular2(1, 10);
  peaksWS->setInstrument(inst);

  for (int i = 0; i < numPeaks; ++i) {
    Peak peak(inst, i, i + 0.5);
    peaksWS->addPeak(peak);
  }

  if (createOrientedLattice) {
    Mantid::Geometry::OrientedLattice lattice;
    peaksWS->mutableSample().setOrientedLattice(&lattice);
  }
  return peaksWS;
}

/** helper method to create preprocessed detector's table workspace */
boost::shared_ptr<DataObjects::TableWorkspace>
createTableWorkspace(const API::MatrixWorkspace_const_sptr &inputWS) {
  const size_t nHist = inputWS->getNumberHistograms();

  // set the target workspace
  auto targWS = boost::make_shared<TableWorkspace>(nHist);
  // detectors positions
  if (!targWS->addColumn("V3D", "DetDirections"))
    throw(std::runtime_error("Can not add column DetDirectrions"));
  // sample-detector distance;
  if (!targWS->addColumn("double", "L2"))
    throw(std::runtime_error("Can not add column L2"));
  // Diffraction angle
  if (!targWS->addColumn("double", "TwoTheta"))
    throw(std::runtime_error("Can not add column TwoTheta"));
  if (!targWS->addColumn("double", "Azimuthal"))
    throw(std::runtime_error("Can not add column Azimuthal"));
  // the detector ID;
  if (!targWS->addColumn("int", "DetectorID"))
    throw(std::runtime_error("Can not add column DetectorID"));
  // stores spectra index which corresponds to a valid detector index;
  if (!targWS->addColumn("size_t", "detIDMap"))
    throw(std::runtime_error("Can not add column detIDMap"));
  // stores detector index which corresponds to the workspace index;
  if (!targWS->addColumn("size_t", "spec2detMap"))
    throw(std::runtime_error("Can not add column spec2detMap"));

  // will see about that
  // sin^2(Theta)
  //    std::vector<double>      SinThetaSq;

  //,"If the detectors were actually processed from real instrument or generated
  // for some fake one ");
  return targWS;
}

/** method does preliminary calculations of the detectors positions to convert
 results into k-dE space ;
 and places the resutls into static cash to be used in subsequent calls to this
 algorithm */
void processDetectorsPositions(const API::MatrixWorkspace_const_sptr &inputWS,
                               DataObjects::TableWorkspace_sptr &targWS,
                               double Ei) {
  Geometry::Instrument_const_sptr instrument = inputWS->getInstrument();
  //
  Geometry::IComponent_const_sptr source = instrument->getSource();
  Geometry::IComponent_const_sptr sample = instrument->getSample();
  if ((!source) || (!sample)) {

    throw Kernel::Exception::InstrumentDefinitionError(
        "Instrubment not sufficiently defined: failed to get source and/or "
        "sample");
  }

  // L1
  try {
    double L1 = source->getDistance(*sample);
    targWS->logs()->addProperty<double>("L1", L1, true);
  } catch (Kernel::Exception::NotFoundError &) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Unable to calculate source-sample distance for workspace",
        inputWS->getTitle());
  }
  // Instrument name
  std::string InstrName = instrument->getName();
  targWS->logs()->addProperty<std::string>("InstrumentName", InstrName, true);
  targWS->logs()->addProperty<bool>("FakeDetectors", false, true);
  targWS->logs()->addProperty<double>("Ei", Ei, true); //"Incident energy for
  // Direct or Analysis
  // energy for indirect
  // instrument");

  // get access to the workspace memory
  auto &sp2detMap = targWS->getColVector<size_t>("spec2detMap");
  auto &detId = targWS->getColVector<int32_t>("DetectorID");
  auto &detIDMap = targWS->getColVector<size_t>("detIDMap");
  auto &L2 = targWS->getColVector<double>("L2");
  auto &TwoTheta = targWS->getColVector<double>("TwoTheta");
  auto &Azimuthal = targWS->getColVector<double>("Azimuthal");
  auto &detDir = targWS->getColVector<Kernel::V3D>("DetDirections");

  //// progress messave appearence
  size_t nHist = targWS->rowCount();
  //// Loop over the spectra
  uint32_t liveDetectorsCount(0);
  const auto &spectrumInfo = inputWS->spectrumInfo();
  for (size_t i = 0; i < nHist; i++) {
    sp2detMap[i] = std::numeric_limits<size_t>::quiet_NaN();
    detId[i] = std::numeric_limits<int32_t>::quiet_NaN();
    detIDMap[i] = std::numeric_limits<size_t>::quiet_NaN();
    L2[i] = std::numeric_limits<double>::quiet_NaN();
    TwoTheta[i] = std::numeric_limits<double>::quiet_NaN();
    Azimuthal[i] = std::numeric_limits<double>::quiet_NaN();

    if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i))
      continue;

    // calculate the requested values;
    sp2detMap[i] = liveDetectorsCount;
    detId[liveDetectorsCount] = int32_t(spectrumInfo.detector(i).getID());
    detIDMap[liveDetectorsCount] = i;
    L2[liveDetectorsCount] = spectrumInfo.l2(i);

    double polar = spectrumInfo.twoTheta(i);
    double azim = spectrumInfo.detector(i).getPhi();
    TwoTheta[liveDetectorsCount] = polar;
    Azimuthal[liveDetectorsCount] = azim;

    double sPhi = sin(polar);
    double ez = cos(polar);
    double ex = sPhi * cos(azim);
    double ey = sPhi * sin(azim);

    detDir[liveDetectorsCount].setX(ex);
    detDir[liveDetectorsCount].setY(ey);
    detDir[liveDetectorsCount].setZ(ez);

    // double sinTheta=sin(0.5*polar);
    // this->SinThetaSq[liveDetectorsCount]  = sinTheta*sinTheta;

    liveDetectorsCount++;
  }
  targWS->logs()->addProperty<uint32_t>(
      "ActualDetectorsNum", liveDetectorsCount,
      true); //,"The actual number of detectors receivinv signal");
}

boost::shared_ptr<Mantid::DataObjects::TableWorkspace>
buildPreprocessedDetectorsWorkspace(Mantid::API::MatrixWorkspace_sptr ws) {
  Mantid::DataObjects::TableWorkspace_sptr DetPos = createTableWorkspace(ws);
  double Ei = ws->run().getPropertyValueAsType<double>("Ei");
  processDetectorsPositions(ws, DetPos, Ei);

  return DetPos;
}
void create2DAngles(std::vector<double> &L2, std::vector<double> &polar,
                    std::vector<double> &azim, size_t nPolar, size_t nAzim,
                    double polStart, double polEnd, double azimStart,
                    double azimEnd) {
  size_t nDet = nPolar * nAzim;
  L2.resize(nDet, 10);
  polar.resize(nDet);
  azim.resize(nDet);

  double dPolar = (polEnd - polStart) / static_cast<double>(nDet - 1);
  double dAzim = (azimEnd - azimEnd) / static_cast<double>(nDet - 1);
  for (size_t i = 0; i < nPolar; i++) {
    for (size_t j = 0; j < nAzim; j++) {
      polar[i * nPolar + j] = polStart + dPolar * static_cast<double>(i);
      azim[i * nPolar + j] = azimStart + dAzim * static_cast<double>(j);
    }
  }
}

ITableWorkspace_sptr
createEPPTableWorkspace(const std::vector<EPPTableRow> &rows) {
  ITableWorkspace_sptr ws = boost::make_shared<TableWorkspace>(rows.size());
  auto wsIndexColumn = ws->addColumn("int", "WorkspaceIndex");
  auto centreColumn = ws->addColumn("double", "PeakCentre");
  auto centreErrorColumn = ws->addColumn("double", "PeakCentreError");
  auto sigmaColumn = ws->addColumn("double", "Sigma");
  auto sigmaErrorColumn = ws->addColumn("double", "SigmaError");
  auto heightColumn = ws->addColumn("double", "Height");
  auto heightErrorColumn = ws->addColumn("double", "HeightError");
  auto chiSqColumn = ws->addColumn("double", "chiSq");
  auto statusColumn = ws->addColumn("str", "FitStatus");
  for (size_t i = 0; i != rows.size(); ++i) {
    const auto &row = rows[i];
    if (row.workspaceIndex < 0) {
      wsIndexColumn->cell<int>(i) = static_cast<int>(i);
    } else {
      wsIndexColumn->cell<int>(i) = row.workspaceIndex;
    }
    centreColumn->cell<double>(i) = row.peakCentre;
    centreErrorColumn->cell<double>(i) = row.peakCentreError;
    sigmaColumn->cell<double>(i) = row.sigma;
    sigmaErrorColumn->cell<double>(i) = row.sigmaError;
    heightColumn->cell<double>(i) = row.height;
    heightErrorColumn->cell<double>(i) = row.heightError;
    chiSqColumn->cell<double>(i) = row.chiSq;
    statusColumn->cell<std::string>(i) =
        row.fitStatus == EPPTableRow::FitStatus::SUCCESS ? "success" : "failed";
  }
  return ws;
}
} // namespace WorkspaceCreationHelper