/** Loads, checks and passes back the values passed to the algorithm
 * @param whiteBeam1 :: A white beam vanadium spectrum that will be used to
 * check detector efficiency variations
 * @param whiteBeam2 :: The other white beam vanadium spectrum from the same
 * instrument to use for comparison
 * @param variation :: The maximum fractional variation above the median that is
 * allowed for god detectors
 * @param startWsIndex :: Index number of the first spectrum to use
 * @param endWsIndex :: Index number of the last spectrum to use
 * @throw invalid_argument if there is an incapatible property value and so the
 * algorithm can't continue
 */
void DetectorEfficiencyVariation::retrieveProperties(
    API::MatrixWorkspace_sptr &whiteBeam1,
    API::MatrixWorkspace_sptr &whiteBeam2, double &variation, int &startWsIndex,
    int &endWsIndex) {
  whiteBeam1 = getProperty("WhiteBeamBase");
  whiteBeam2 = getProperty("WhiteBeamCompare");
  if (whiteBeam1->getInstrument()->getName() !=
      whiteBeam2->getInstrument()->getName()) {
    throw std::invalid_argument("The two input white beam vanadium workspaces "
                                "must be from the same instrument");
  }
  int maxWsIndex = static_cast<int>(whiteBeam1->getNumberHistograms()) - 1;
  if (maxWsIndex !=
      static_cast<int>(whiteBeam2->getNumberHistograms()) -
          1) { // we would get a crash later on if this were not true
    throw std::invalid_argument("The input white beam vanadium workspaces must "
                                "be have the same number of histograms");
  }

  variation = getProperty("Variation");

  startWsIndex = getProperty("StartWorkspaceIndex");
  if ((startWsIndex < 0) || (startWsIndex > maxWsIndex)) {
    g_log.warning("StartWorkspaceIndex out of range, changed to 0");
    startWsIndex = 0;
  }
  endWsIndex = getProperty("EndWorkspaceIndex");
  if (endWsIndex == Mantid::EMPTY_INT())
    endWsIndex = maxWsIndex;
  if ((endWsIndex < 0) || (endWsIndex > maxWsIndex)) {
    g_log.warning(
        "EndWorkspaceIndex out of range, changed to max Workspace number");
    endWsIndex = maxWsIndex;
  }
  if ((endWsIndex < startWsIndex)) {
    g_log.warning(
        "EndWorkspaceIndex can not be less than the StartWorkspaceIndex, "
        "changed to max Workspace number");
    endWsIndex = maxWsIndex;
  }
}

/// Validate the some of the algorithm's input properties.
std::map<std::string, std::string> ReflectometrySumInQ::validateInputs() {
  std::map<std::string, std::string> issues;
  API::MatrixWorkspace_sptr inWS;
  Indexing::SpectrumIndexSet indices;

  // validateInputs is called on the individual workspaces when the algorithm
  // is executed, it but may get called on a group from AlgorithmDialog. This
  // isn't handled in getWorkspaceAndIndices. We should fix this properly but
  // for now skip validation for groups to avoid an exception. See #22933
  try {
    std::tie(inWS, indices) =
        getWorkspaceAndIndices<API::MatrixWorkspace>(Prop::INPUT_WS);
  } catch (std::runtime_error &) {
    return issues;
  }

  const auto &spectrumInfo = inWS->spectrumInfo();
  const double beamCentre = getProperty(Prop::BEAM_CENTRE);
  const size_t beamCentreIndex = static_cast<size_t>(beamCentre);
  bool beamCentreFound{false};
  for (const auto i : indices) {
    if (spectrumInfo.isMonitor(i)) {
      issues["InputWorkspaceIndexSet"] = "Index set cannot include monitors.";
      break;
    } else if ((i > 0 && spectrumInfo.isMonitor(i - 1)) ||
               (i < spectrumInfo.size() - 1 && spectrumInfo.isMonitor(i + 1))) {
      issues["InputWorkspaceIndexSet"] =
          "A neighbour to any detector in the index set cannot be a monitor";
      break;
    }
    if (i == beamCentreIndex) {
      beamCentreFound = true;
      break;
    }
  }
  if (!beamCentreFound) {
    issues[Prop::BEAM_CENTRE] =
        "Beam centre is not included in InputWorkspaceIndexSet.";
  }
  return issues;
}

/*
 * Convert DAS log to a vector of absolute time
 * @param  orderedtofs: tofs with abstimevec
 */
void ProcessDasNexusLog::convertToAbsoluteTime(
    API::MatrixWorkspace_sptr ws, std::string logname,
    std::vector<Kernel::DateAndTime> &abstimevec,
    std::vector<double> &orderedtofs) {
  // 1. Get log
  Kernel::Property *log = ws->run().getProperty(logname);
  Kernel::TimeSeriesProperty<double> *tslog =
      dynamic_cast<Kernel::TimeSeriesProperty<double> *>(log);
  if (!tslog)
    throw std::runtime_error("Invalid time series log: it could not be cast "
                             "(interpreted) as a time series property");
  std::vector<Kernel::DateAndTime> times = tslog->timesAsVector();
  std::vector<double> values = tslog->valuesAsVector();

  // 2. Get converted
  size_t numsamepulses = 0;
  std::vector<double> tofs;
  Kernel::DateAndTime prevtime(0);

  for (size_t i = 0; i < times.size(); i++) {
    Kernel::DateAndTime tnow = times[i];
    if (tnow > prevtime) {
      // (a) Process previous logs
      std::sort(tofs.begin(), tofs.end());
      for (double tof : tofs) {
        Kernel::DateAndTime temptime =
            prevtime + static_cast<int64_t>(tof * 100);
        abstimevec.push_back(temptime);
        orderedtofs.push_back(tof);
      }
      // (b) Clear
      tofs.clear();
      // (c) Update time
      prevtime = tnow;
    } else {
      numsamepulses++;
    }
    // (d) Push the current value
    tofs.push_back(values[i]);
  } // ENDFOR
  // Clear the last
  if (!tofs.empty()) {
    // (a) Process previous logs: note value is in unit of 100 nano-second
    std::sort(tofs.begin(), tofs.end());
    for (double tof : tofs) {
      Kernel::DateAndTime temptime = prevtime + static_cast<int64_t>(tof * 100);
      abstimevec.push_back(temptime);
      orderedtofs.push_back(tof);
    }
  } else {
    throw std::runtime_error("Impossible for this to happen!");
  }
} // END Function

    std::vector<std::vector<size_t> > DetectorDiagnostic::makeInstrumentMap(API::MatrixWorkspace_sptr countsWS)
    {
      std::vector<std::vector<size_t> > mymap;
      std::vector<size_t> single;

      for(size_t i=0;i < countsWS->getNumberHistograms();i++)
      {
        single.push_back(i);
      }
      mymap.push_back(single);
      return mymap;
    }

void ConvertEmptyToTof::setTofInWS(const std::vector<double> &tofAxis,
                                   API::MatrixWorkspace_sptr outputWS) {

  const size_t numberOfSpectra = m_inputWS->getNumberHistograms();

  g_log.debug() << "Setting the TOF X Axis for numberOfSpectra="
                << numberOfSpectra << '\n';

  auto axisPtr = Kernel::make_cow<HistogramData::HistogramX>(tofAxis);
  HistogramData::BinEdges edges(tofAxis);
  Progress prog(this, 0.0, 0.2, numberOfSpectra);

  for (size_t i = 0; i < numberOfSpectra; ++i) {
    // Replace bin edges with tof axis
    outputWS->setBinEdges(i, edges);

    prog.report();
  } // end for i

  outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
}

/** Smoothing using Butterworth filter.
 *  @param n ::     The cutoff frequency control parameter.
 *               Cutoff frequency = my/n where my is the
 *               number of sample points in the data.
 *               As with the "Zeroing" case, the cutoff
 *               frequency is truncated to an integer value
 *               and set to 1 if the truncated value was zero.
 *  @param order :: The order of the Butterworth filter, 1, 2, etc.
 *               This must be a positive integer.
 *  @param unfilteredWS :: workspace for storing the unfiltered Fourier
 * transform of the input spectrum
 *  @param filteredWS :: workspace for storing the filtered spectrum
 */
void FFTSmooth2::Butterworth(int n, int order,
                             API::MatrixWorkspace_sptr &unfilteredWS,
                             API::MatrixWorkspace_sptr &filteredWS) {
  int mx = static_cast<int>(unfilteredWS->readX(0).size());
  int my = static_cast<int>(unfilteredWS->readY(0).size());
  int ny = my / n;

  if (ny == 0)
    ny = 1;

  filteredWS =
      API::WorkspaceFactory::Instance().create(unfilteredWS, 2, mx, my);

  const Mantid::MantidVec &Yr = unfilteredWS->readY(0);
  const Mantid::MantidVec &Yi = unfilteredWS->readY(1);
  const Mantid::MantidVec &X = unfilteredWS->readX(0);

  Mantid::MantidVec &yr = filteredWS->dataY(0);
  Mantid::MantidVec &yi = filteredWS->dataY(1);
  Mantid::MantidVec &xr = filteredWS->dataX(0);
  Mantid::MantidVec &xi = filteredWS->dataX(1);

  xr.assign(X.begin(), X.end());
  xi.assign(X.begin(), X.end());
  yr.assign(Yr.size(), 0);
  yi.assign(Yr.size(), 0);

  double cutoff = ny;

  for (int i = 0; i < my; i++) {
    double scale = 1.0 / (1.0 + pow(i / cutoff, 2 * order));
    yr[i] = scale * Yr[i];
    yi[i] = scale * Yi[i];
  }
}

/** Smoothing by zeroing.
 *  @param n :: The order of truncation
 *  @param unfilteredWS :: workspace for storing the unfiltered Fourier
 * transform of the input spectrum
 *  @param filteredWS :: workspace for storing the filtered spectrum
 */
void FFTSmooth2::zero(int n, API::MatrixWorkspace_sptr &unfilteredWS,
                      API::MatrixWorkspace_sptr &filteredWS) {
  int mx = static_cast<int>(unfilteredWS->readX(0).size());
  int my = static_cast<int>(unfilteredWS->readY(0).size());
  int ny = my / n;

  if (ny == 0)
    ny = 1;

  filteredWS =
      API::WorkspaceFactory::Instance().create(unfilteredWS, 2, mx, my);

  const Mantid::MantidVec &Yr = unfilteredWS->readY(0);
  const Mantid::MantidVec &Yi = unfilteredWS->readY(1);
  const Mantid::MantidVec &X = unfilteredWS->readX(0);

  Mantid::MantidVec &yr = filteredWS->dataY(0);
  Mantid::MantidVec &yi = filteredWS->dataY(1);
  Mantid::MantidVec &xr = filteredWS->dataX(0);
  Mantid::MantidVec &xi = filteredWS->dataX(1);

  xr.assign(X.begin(), X.end());
  xi.assign(X.begin(), X.end());
  yr.assign(Yr.size(), 0);
  yi.assign(Yr.size(), 0);

  for (int i = 0; i < ny; i++) {
    yr[i] = Yr[i];
    yi[i] = Yi[i];
  }
}

/***
 * This will ensure the spectrum numbers do not overlap by starting the second
 *on at the first + 1
 *
 * @param ws1 The first workspace supplied to the algorithm.
 * @param ws2 The second workspace supplied to the algorithm.
 * @param output The workspace that is going to be returned by the algorithm.
 */
void ConjoinWorkspaces::fixSpectrumNumbers(API::MatrixWorkspace_const_sptr ws1,
                                           API::MatrixWorkspace_const_sptr ws2,
                                           API::MatrixWorkspace_sptr output) {
  bool needsFix(false);

  if (this->getProperty("CheckOverlapping")) {
    // If CheckOverlapping is required, then either skip fixing spectrum number
    // or get stopped by an exception
    if (!m_overlapChecked)
      checkForOverlap(ws1, ws2, true);
    needsFix = false;
  } else {
    // It will be determined later whether spectrum number needs to be fixed.
    needsFix = true;
  }
  if (!needsFix)
    return;

  // is everything possibly ok?
  specid_t min;
  specid_t max;
  getMinMax(output, min, max);
  if (max - min >= static_cast<specid_t>(
                       output->getNumberHistograms())) // nothing to do then
    return;

  // information for remapping the spectra numbers
  specid_t ws1min;
  specid_t ws1max;
  getMinMax(ws1, ws1min, ws1max);

  // change the axis by adding the maximum existing spectrum number to the
  // current value
  for (size_t i = ws1->getNumberHistograms(); i < output->getNumberHistograms();
       i++) {
    specid_t origid;
    origid = output->getSpectrum(i)->getSpectrumNo();
    output->getSpectrum(i)->setSpectrumNo(origid + ws1max);
  }
}

/** Checks and retrieves the requested spectrum out of the input workspace
 *  @param inputWorkspace The input workspace.
 *  @param spectra_num The spectra number.
 *  @returns A workspace containing the monitor spectrum only.
 *  @returns spectra number (WS ID) which is used to normalize by.
 *  @throw std::runtime_error If the properties are invalid
 */
API::MatrixWorkspace_sptr NormaliseToMonitor::getInWSMonitorSpectrum(API::MatrixWorkspace_sptr inputWorkspace,int &spectra_num)
{
 // this is the index of the spectra within the workspace and we need to indetnify it either from DetID or fron SpecID
 // size_t spectra_num(-1);
// try monitor spectrum. If it is specified, it overides everything
  int monitorSpec = getProperty("MonitorSpectrum");
  if(monitorSpec<0){ 
      // Get hold of the monitor spectrum through detector ID
    int monitorID = getProperty("MonitorID");
    if (monitorID < 0) {
        throw std::runtime_error("Both MonitorSpectrum and MonitorID can not be negative");
    }
    // set spectra of detector's ID of one selected monitor ID
    std::vector<detid_t> detID(1,monitorID);
    // got the index of correspondent spectras (should be only one). 
    std::vector<size_t>  indexList;
    inputWorkspace->getIndicesFromDetectorIDs(detID,indexList);
    if(indexList.empty()){
         throw std::runtime_error("Can not find spectra, coorespoinding to the requested monitor ID");
    }
    if(indexList.size()>1){
         throw std::runtime_error("More then one spectra coorespods to the requested monitor ID, which is unheard of");
    }
    spectra_num = (int)indexList[0];
  }else{ // monitor spectrum is specified.
        spec2index_map specs;
        const SpectraAxis* axis = dynamic_cast<const SpectraAxis*>(inputWorkspace->getAxis(1));
        if ( ! axis)
        {
            throw std::runtime_error("Cannot retrieve monitor spectrum - spectrum numbers not attached to workspace");
        }
        axis->getSpectraIndexMap(specs);
        if ( ! specs.count(monitorSpec) )
        {
            throw std::runtime_error("Input workspace does not contain spectrum number given for MonitorSpectrum");
        }
        spectra_num = (int)specs[monitorSpec];
  }
  return this->extractMonitorSpectrum(inputWorkspace,spectra_num);
}

/**  Load logs from Nexus file. Logs are expected to be in
*   /run/sample group of the file.
*   @param ws :: The workspace to load the logs to.
*   @param entry :: The Nexus entry
*   @param period :: The period of this workspace
*/
void LoadMuonNexus2::loadLogs(API::MatrixWorkspace_sptr ws, NXEntry &entry,
                              int period) {
  // Avoid compiler warning
  (void)period;

  std::string start_time = entry.getString("start_time");

  std::string sampleName = entry.getString("sample/name");
  NXMainClass runlogs = entry.openNXClass<NXMainClass>("sample");
  ws->mutableSample().setName(sampleName);

  for (std::vector<NXClassInfo>::const_iterator it = runlogs.groups().begin();
       it != runlogs.groups().end(); ++it) {
    NXLog nxLog = runlogs.openNXLog(it->nxname);
    Kernel::Property *logv = nxLog.createTimeSeries(start_time);
    if (!logv)
      continue;
    ws->mutableRun().addLogData(logv);
  }

  ws->setTitle(entry.getString("title"));

  if (entry.containsDataSet("notes")) {
    ws->setComment(entry.getString("notes"));
  }

  std::string run_num = std::to_string(entry.getInt("run_number"));
  // The sample is left to delete the property
  ws->mutableRun().addLogData(
      new PropertyWithValue<std::string>("run_number", run_num));

  ws->populateInstrumentParameters();
}

/** Execute the algorithm.
 */
void CalculateDIFC::exec() {

  DataObjects::OffsetsWorkspace_const_sptr offsetsWs =
      getProperty("OffsetsWorkspace");
  API::ITableWorkspace_const_sptr calibWs = getProperty("CalibrationWorkspace");
  API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
  API::MatrixWorkspace_sptr outputWs = getProperty("OutputWorkspace");

  if ((!bool(inputWs == outputWs)) ||
      (!bool(boost::dynamic_pointer_cast<SpecialWorkspace2D>(outputWs)))) {
    outputWs = boost::dynamic_pointer_cast<MatrixWorkspace>(
        boost::make_shared<SpecialWorkspace2D>(inputWs->getInstrument()));
    outputWs->setTitle("DIFC workspace");
  }

  // convert to actual type being used
  DataObjects::SpecialWorkspace2D_sptr outputSpecialWs =
      boost::dynamic_pointer_cast<DataObjects::SpecialWorkspace2D>(outputWs);

  API::Progress progress(this, 0.0, 1.0, inputWs->getNumberHistograms());
  if (bool(calibWs)) {
    calculateFromTable(progress, *outputSpecialWs, *calibWs);
  } else {
    // this method handles calculating from instrument geometry as well
    const auto &detectorInfo = inputWs->detectorInfo();
    calculateFromOffset(progress, *outputSpecialWs, offsetsWs.get(),
                        detectorInfo);
  }

  setProperty("OutputWorkspace", outputWs);
}

/**
 * Here is the main logic to perform the transformation, to calculate the bin position in degree for each spectrum.
 * 
 * The first part of the method is to check if the pixel position is inside the ring defined as minRadio and maxRadio. 
 * 
 * To do this, it deducts the pixel position. This deduction follows the followin assumption: 
 * 
 *  - the spectrum_index == row number 
 *  - the position in the 'Y' direction is given by getAxis(1)[spectrum_index]
 *  - the position in the 'X' direction is the central point of the bin (dataX[column] + dataX[column+1])/2
 *
 * Having the position of the pixel, as defined above, if the distance is outside the ring defined by minRadio, maxRadio, 
 * it defines the bin position as -1.
 * 
 * If the pixel is inside the ring, it calculates the angle of the pixel and calls fromAngleToBin to define the bin 
 * position. 
 * @param ws: pointer to the workspace
 * @param spectrum_index: index of the spectrum
 * @param bins_pos: bin positions (for each column inside the spectrum, the correspondent bin_pos)
 */
void RingProfile::getBinForPixel(const API::MatrixWorkspace_sptr ws, 
                                 int spectrum_index, std::vector<int> & bins_pos ){
  
  if (bins_pos.size() != ws->dataY(spectrum_index).size())
    throw std::runtime_error("Invalid bin positions vector"); 
  
  API::NumericAxis *oldAxis2 = dynamic_cast<API::NumericAxis*>(ws->getAxis(1));
  // assumption y position is the ws->getAxis(1)(spectrum_index)

  // calculate ypos, the difference of y - centre and  the square of this difference
  double ypos = (*oldAxis2)(spectrum_index);
  double diffy = ypos-centre_y;
  double diffy_quad = pow(diffy, 2.0);

  // the reference to X bins (the limits for each pixel in the horizontal direction)
  auto xvec = ws->dataX(spectrum_index);

  // for each pixel inside this row
  for (size_t i = 0; i< xvec.size()-1; i++){

    double xpos = (xvec[i] + xvec[i+1])/2.0; // the x position is the centre of the bins boundaries
    double diffx = xpos - centre_x; 
    // calculate the distance => norm of pixel position - centre
    double distance = sqrt(pow(diffx, 2.0) + diffy_quad);
    
    // check if the distance is inside the ring
    if (distance < min_radius || distance > max_radius || distance == 0){
      bins_pos[i] = -1; 
      continue;
    }
    
    double angle = atan2(diffy, diffx);

    // call fromAngleToBin (radians)
    bins_pos[i] = fromAngleToBin(angle, false); 

  } 
  
  }; 

/*
 * Get instrument property as double
 * @s - input property name
 *
 */
double
LoadHelper::getInstrumentProperty(const API::MatrixWorkspace_sptr &workspace,
                                  std::string s) {
    std::vector<std::string> prop =
        workspace->getInstrument()->getStringParameter(s);
    if (prop.empty()) {
        g_log.debug("Property <" + s + "> doesn't exist!");
        return EMPTY_DBL();
    } else {
        g_log.debug() << "Property <" + s + "> = " << prop[0] << '\n';
        return boost::lexical_cast<double>(prop[0]);
    }
}

 V3D LoadHelper::getComponentPosition(API::MatrixWorkspace_sptr ws, const std::string &componentName)
 {
   try
   {
     Geometry::Instrument_const_sptr instrument = ws->getInstrument();
     Geometry::IComponent_const_sptr component = instrument->getComponentByName(componentName);
     V3D pos = component->getPos();
     return pos;
   } catch (Mantid::Kernel::Exception::NotFoundError&)
   {
     throw std::runtime_error("Error when trying to move the " + componentName + " : NotFoundError");
   }
 }

/** Convert the workspace units according to a simple output = a * (input^b) relationship
 *  @param outputWS :: the output workspace
 *  @param factor :: the conversion factor a to apply
 *  @param power :: the Power b to apply to the conversion
 */
void ConvertUnits::convertQuickly(API::MatrixWorkspace_sptr outputWS, const double& factor, const double& power)
{
  Progress prog(this,0.2,1.0,m_numberOfSpectra);
  int64_t numberOfSpectra_i = static_cast<int64_t>(m_numberOfSpectra); // cast to make openmp happy

  // See if the workspace has common bins - if so the X vector can be common
  // First a quick check using the validator
  CommonBinsValidator sameBins;
  bool commonBoundaries = false;
  if ( sameBins.isValid(outputWS) == "" )
  {
    commonBoundaries =  WorkspaceHelpers::commonBoundaries(outputWS);
    // Only do the full check if the quick one passes
    if (commonBoundaries)
    {
      // Calculate the new (common) X values
      MantidVec::iterator iter;
      for (iter = outputWS->dataX(0).begin(); iter != outputWS->dataX(0).end(); ++iter)
      {
        *iter = factor * std::pow(*iter,power);
      }

      MantidVecPtr xVals;
      xVals.access() = outputWS->dataX(0);

      PARALLEL_FOR1(outputWS)
      for (int64_t j = 1; j < numberOfSpectra_i; ++j)
      {
        PARALLEL_START_INTERUPT_REGION
        outputWS->setX(j,xVals);
        prog.report("Convert to " + m_outputUnit->unitID());
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
      if (!m_inputEvents) // if in event mode the work is done
        return;
    }
  }

/**
 * The main method to calculate the ring profile for workspaces based on instruments. 
 * 
 * It will iterate over all the spectrum inside the workspace. 
 * For each spectrum, it will use the RingProfile::getBinForPixel method to identify 
 * where, in the output_bins, the sum of all the spectrum values should be placed in. 
 * 
 * @param inputWS: pointer to the input workspace
 * @param output_bins: the reference to the vector to be filled with the integration values
 */
void RingProfile::processInstrumentRingProfile(const API::MatrixWorkspace_sptr inputWS, 
                                               std::vector<double> & output_bins){

  for (int i= 0; i< (int) inputWS->getNumberHistograms(); i++){
    m_progress->report("Computing ring bins positions for detectors"); 
    // for the detector based, the positions will be taken from the detector itself.
    try{
      Mantid::Geometry::IDetector_const_sptr det = inputWS->getDetector(i);
      
      // skip monitors  
      if (det->isMonitor()){
        continue;
      }
      
      // this part will be executed if the instrument is attached to the workspace
      
      // get the bin position
      int bin_n = getBinForPixel(det); 
        
      if (bin_n < 0) // -1 is the agreement for an invalid bin, or outside the ring being integrated
        continue;       

      g_log.debug() << "Bin for the index " << i << " = " << bin_n << " Pos = " << det->getPos() << std::endl; 

      // get the reference to the spectrum
      auto spectrum_pt = inputWS->getSpectrum(i); 
      const MantidVec & refY = spectrum_pt->dataY(); 
      // accumulate the values of this spectrum inside this bin
      for (size_t sp_ind = 0; sp_ind < inputWS->blocksize(); sp_ind ++)
        output_bins[bin_n] += refY[sp_ind];         

    }catch(Kernel::Exception::NotFoundError & ex){
      g_log.information() << "It found that detector for " << i << " is not valid. " << ex.what() << std::endl; 
      continue;
    }
    
  }
}

/** Performs the Debye-Bueche transformation: 1/sqrt(I) v Q^2
 *  The output is set to zero for negative input Y values
 *  @param ws The workspace to be transformed
 */
void IQTransform::debyeBueche(API::MatrixWorkspace_sptr ws)
{
  MantidVec& X = ws->dataX(0);
  MantidVec& Y = ws->dataY(0);
  MantidVec& E = ws->dataE(0);
  std::transform(X.begin(),X.end(),X.begin(),VectorHelper::Squares<double>());
  for(size_t i = 0; i < Y.size(); ++i)
  {
    if ( Y[i] > 0.0 )
    {
      Y[i] = 1.0 / std::sqrt(Y[i]);
      E[i] *= std::pow(Y[i],3);
    }
    else
    {
      Y[i] = 0.0;
      E[i] = 0.0;
    }
  }

  ws->setYUnitLabel("1/sqrt(I)");
  m_label->setLabel("Q^2");
}

/** Extracts a single spectrum from a Workspace2D into a new workspaces. Uses CropWorkspace to do this.
 *  @param WS ::    The workspace containing the spectrum to extract
 *  @param index :: The workspace index of the spectrum to extract
 *  @return A Workspace2D containing the extracted spectrum
 */
API::MatrixWorkspace_sptr CalculateTransmissionBeamSpreader::extractSpectrum(API::MatrixWorkspace_sptr WS, const size_t index)
{
  // Check that given spectra are monitors
  if ( !WS->getDetector(index)->isMonitor() )
  {
    g_log.information("The Incident Beam Monitor UDET provided is not marked as a monitor");
  }

  Algorithm_sptr childAlg = createSubAlgorithm("ExtractSingleSpectrum",0.0,0.4);
  childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS);
  childAlg->setProperty<int>("WorkspaceIndex", static_cast<int>(index));
  childAlg->executeAsSubAlg();
  return childAlg->getProperty("OutputWorkspace");
}

/**
 * Perform a call to nxgetslab, via the NexusClasses wrapped methods for a given blocksize. This assumes that the
 * xbins have alread been cached
 * @param data :: The NXDataSet object of y values
 * @param errors :: The NXDataSet object of error values
 * @param blocksize :: The blocksize to use
 * @param nchannels :: The number of channels for the block
 * @param hist :: The workspace index to start reading into
 * @param local_workspace :: A pointer to the workspace
 */
void LoadNexusProcessed::loadBlock(NXDataSetTyped<double> & data, NXDataSetTyped<double> & errors,
    int64_t blocksize, int64_t nchannels, int64_t &hist,
    API::MatrixWorkspace_sptr local_workspace)
{
  data.load(static_cast<int>(blocksize),static_cast<int>(hist));
  errors.load(static_cast<int>(blocksize),static_cast<int>(hist));
  double *data_start = data();
  double *data_end = data_start + nchannels;
  double *err_start = errors();
  double *err_end = err_start + nchannels;
  int64_t final(hist + blocksize);
  while( hist < final )
  {
    MantidVec& Y = local_workspace->dataY(hist);
    Y.assign(data_start, data_end);
    data_start += nchannels; data_end += nchannels;
    MantidVec& E = local_workspace->dataE(hist);
    E.assign(err_start, err_end);
    err_start += nchannels; err_end += nchannels;
    local_workspace->setX(hist, m_xbins);
    ++hist;
  }
}

/** Initialization method:
@param bkgWS    -- shared pointer to the workspace which contains background
@param sourceWS -- shared pointer to the workspace to remove background from
@param emode    -- energy conversion mode used during internal units conversion
(0 -- elastic, 1-direct, 2 indirect, as defined in Units conversion
@param pLog     -- pointer to the logger class which would report errors
@param nThreads -- number of threads to be used for background removal
@param inPlace  -- if the background removal occurs from the existing workspace
or target workspace has to be cloned.
*/
void BackgroundHelper::initialize(const API::MatrixWorkspace_const_sptr &bkgWS,
                                  const API::MatrixWorkspace_sptr &sourceWS,
                                  int emode, Kernel::Logger *pLog, int nThreads,
                                  bool inPlace) {
  m_bgWs = bkgWS;
  m_wkWS = sourceWS;
  m_Emode = emode;
  m_pgLog = pLog;
  m_inPlace = inPlace;

  std::string bgUnits = bkgWS->getAxis(0)->unit()->unitID();
  if (bgUnits != "TOF")
    throw std::invalid_argument(" Background Workspace: " + bkgWS->getName() +
                                " should be in the units of TOF");

  if (!(bkgWS->getNumberHistograms() == 1 ||
        sourceWS->getNumberHistograms() == bkgWS->getNumberHistograms()))
    throw std::invalid_argument(" Background Workspace: " + bkgWS->getName() +
                                " should have the same number of spectra as "
                                "source workspace or be a single histogram "
                                "workspace");

  auto WSUnit = sourceWS->getAxis(0)->unit();
  if (!WSUnit)
    throw std::invalid_argument(" Source Workspace: " + sourceWS->getName() +
                                " should have units");

  Geometry::IComponent_const_sptr source =
      sourceWS->getInstrument()->getSource();
  m_Sample = sourceWS->getInstrument()->getSample();
  if ((!source) || (!m_Sample))
    throw std::invalid_argument(
        "Instrument on Source workspace:" + sourceWS->getName() +
        "is not sufficiently defined: failed to get source and/or sample");
  m_L1 = source->getDistance(*m_Sample);

  // just in case.
  this->deleteUnitsConverters();
  // allocate the array of units converters to avoid units reallocation within a
  // loop
  m_WSUnit.assign(nThreads, NULL);
  for (int i = 0; i < nThreads; i++) {
    m_WSUnit[i] = WSUnit->clone();
  }

  m_singleValueBackground = false;
  if (bkgWS->getNumberHistograms() == 0)
    m_singleValueBackground = true;
  const MantidVec &dataX = bkgWS->dataX(0);
  const MantidVec &dataY = bkgWS->dataY(0);
  // const MantidVec& dataE = bkgWS->dataE(0);
  m_NBg = dataY[0];
  m_dtBg = dataX[1] - dataX[0];
  // m_ErrSq  = dataE[0]*dataE[0]; // needs further clarification

  m_Efix = this->getEi(sourceWS);
}