WorkspaceCreationHelper.cpp

  pAxis0->setUnit("DeltaE");

  ws->replaceAxis(0, pAxis0);

  // define oriented lattice which requested for processed ws
  OrientedLattice *latt = new OrientedLattice(1, 1, 1, 90., 90., 90.);
  ws->mutableSample().setOrientedLattice(latt);
  delete latt;

  // TODO: clarify if this property indeed goes there;
  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->setName(wsname);
  outputWS->initialize(1, 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
  std::vector<std::size_t> pixel_to_wkspindex;
  pixel_to_wkspindex.reserve(
      detid_max + 1); // starting at zero up to and including detid_max
  pixel_to_wkspindex.assign(detid_max + 1, 0);
  size_t workspaceIndex = 0;
  for (it = detector_map.begin(); it != detector_map.end(); ++it) {
    if (!it->second->isMonitor()) {
      pixel_to_wkspindex[it->first] = workspaceIndex;
      DataObjects::EventList &spec =
          outputWS->getOrAddEventList(workspaceIndex);
      spec.addDetectorID(it->first);
      // Start the spectrum number at 1
      spec.setSpectrumNo(specid_t(workspaceIndex + 1));
      workspaceIndex += 1;
    }
  }

  // Clear
  pixel_to_wkspindex.clear();

  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
  MantidVec qbins;
  qbins.push_back(0.0);
  qbins.push_back(1.0);
  qbins.push_back(4.0);
  MantidVec 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 x-axis for histogramming.
  MantidVecPtr x1;
  MantidVec &xRef = x1.access();
  double x0 = -3;
  xRef.resize(numX);
  for (int i = 0; i < numX; ++i) {
    xRef[i] = x0 + i;
  }

  // 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->setX(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
  outputWS->dataY(0)[1] = 2.0;
  outputWS->dataY(0)[2] = 3.0;
  outputWS->dataY(0)[3] = 3.0;
  outputWS->dataY(0)[4] = 2.0;
  outputWS->dataE(0)[1] = 2.0;
  outputWS->dataE(0)[2] = 3.0;
  outputWS->dataE(0)[3] = 3.0;
  outputWS->dataE(0)[4] = 2.0;
  outputWS->dataF(0)[1] = 2.0;
  outputWS->dataF(0)[2] = 3.0;
  outputWS->dataF(0)[3] = 3.0;
  outputWS->dataF(0)[4] = 1.0;
  // Q bin #2
  outputWS->dataY(1)[1] = 1.0;
  outputWS->dataY(1)[2] = 3.0;
  outputWS->dataY(1)[3] = 3.0;
  outputWS->dataY(1)[4] = 2.0;
  outputWS->dataY(1)[5] = 2.0;
  outputWS->dataE(1)[1] = 1.0;
  outputWS->dataE(1)[2] = 3.0;
  outputWS->dataE(1)[3] = 3.0;
  outputWS->dataE(1)[4] = 2.0;
  outputWS->dataE(1)[5] = 2.0;
  outputWS->dataF(1)[1] = 1.0;
  outputWS->dataF(1)[2] = 3.0;
  outputWS->dataF(1)[3] = 3.0;
  outputWS->dataF(1)[4] = 1.0;
  outputWS->dataF(1)[5] = 2.0;
  // Q bin #3
  outputWS->dataY(2)[1] = 1.0;
  outputWS->dataY(2)[2] = 2.0;
  outputWS->dataY(2)[3] = 3.0;
  outputWS->dataY(2)[4] = 1.0;
  outputWS->dataE(2)[1] = 1.0;
  outputWS->dataE(2)[2] = 2.0;
  outputWS->dataE(2)[3] = 3.0;
  outputWS->dataE(2)[4] = 1.0;
  outputWS->dataF(2)[1] = 1.0;
  outputWS->dataF(2)[2] = 2.0;
  outputWS->dataF(2)[3] = 2.0;
  outputWS->dataF(2)[4] = 1.0;
  // Q bin #4
  outputWS->dataY(3)[0] = 1.0;
  outputWS->dataY(3)[1] = 2.0;
  outputWS->dataY(3)[2] = 3.0;
  outputWS->dataY(3)[3] = 2.0;
  outputWS->dataY(3)[4] = 1.0;
  outputWS->dataE(3)[0] = 1.0;
  outputWS->dataE(3)[1] = 2.0;
  outputWS->dataE(3)[2] = 3.0;
  outputWS->dataE(3)[3] = 2.0;
  outputWS->dataE(3)[4] = 1.0;
  outputWS->dataF(3)[0] = 1.0;
  outputWS->dataF(3)[1] = 2.0;
  outputWS->dataF(3)[2] = 3.0;
  outputWS->dataF(3)[3] = 2.0;
  outputWS->dataF(3)[4] = 1.0;
  outputWS->dataF(3)[5] = 1.0;

  // Set representation
  outputWS->finalize();

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

  return outputWS;
}

Mantid::DataObjects::PeaksWorkspace_sptr
createPeaksWorkspace(const int numPeaks, const bool createOrientedLattice) {
  PeaksWorkspace_sptr peaksWS(new 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::shared_ptr<TableWorkspace>(new 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);
  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();

    // get detector or detector group which corresponds to the spectra i
    Geometry::IDetector_const_sptr spDet;
    try {
      spDet = inputWS->getDetector(i);
    } catch (Kernel::Exception::NotFoundError &) {
      continue;
    }

    // Check that we aren't dealing with monitor...
    if (spDet->isMonitor())
      continue;

    // calculate the requested values;
    sp2detMap[i] = liveDetectorsCount;
    detId[liveDetectorsCount] = int32_t(spDet->getID());
    detIDMap[liveDetectorsCount] = i;
    L2[liveDetectorsCount] = spDet->getDistance(*sample);

    double polar = inputWS->detectorTwoTheta(spDet);
    double azim = spDet->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);
    }
  }
}
}