timsframe.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/timsframe.cpp
 * \date 23/08/2019
 * \author Olivier Langella
 * \brief handle a single Bruker's TimsTof frame
 */
 
/*******************************************************************************
 * Copyright (c) 2019 Olivier Langella <Olivier.Langella@u-psud.fr>.
 *
 * This file is part of the PAPPSOms++ library.
 *
 *     PAPPSOms++ is free software: you can redistribute it and/or modify
 *     it under the terms of the GNU General Public License as published by
 *     the Free Software Foundation, either version 3 of the License, or
 *     (at your option) any later version.
 *
 *     PAPPSOms++ is distributed in the hope that it will be useful,
 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *     GNU General Public License for more details.
 *
 *     You should have received a copy of the GNU General Public License
 *     along with PAPPSOms++.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
#include "timsframe.h"
#include "pappsomspp/core/pappsoexception.h"
#include <QDebug>
#include <QObject>
#include <QtEndian>
#include <memory>
 
namespace pappso
{
 
TimsFrame::XicComputeStructure::XicComputeStructure(
  const TimsFrame *fram_p, const XicCoordTims &xic_struct)
{
  xic_ptr = xic_struct.xicSptr.get();
 
  mobilityIndexBegin = xic_struct.scanNumBegin;
  mobilityIndexEnd   = xic_struct.scanNumEnd;
  mzIndexLowerBound =
    fram_p->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(
      xic_struct.mzRange.lower()); // convert mz to raw digitizer value
  mzIndexUpperBound =
    fram_p->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(
      xic_struct.mzRange.upper());
  tmpIntensity = 0;
}
 
TimsFrame::TimsFrame(std::size_t timsId, quint32 scanNum)
  : TimsFrameBase(timsId, scanNum)
{
  // m_timsDataFrame.resize(10);
}
 
TimsFrame::TimsFrame(std::size_t timsId,
                     quint32 scanNum,
                     char *p_bytes,
                     std::size_t len)
  : TimsFrameBase(timsId, scanNum)
{
  // langella@themis:~/developpement/git/bruker/cbuild$
  // ./src/sample/timsdataSamplePappso
  // /gorgone/pappso/fichiers_fabricants/Bruker/Demo_TimsTOF_juin2019/Samples/1922001/1922001-1_S-415_Pep_Pur-1ul_Slot1-10_1_2088.d/
  // qDebug() << timsId;
 
  m_binaryData.resize(len);
 
  if(p_bytes != nullptr)
    {
      unshufflePacket(p_bytes);
    }
  else
    {
      if(m_scanCount == 0)
        {
 
          throw pappso::PappsoException(
            QObject::tr("TimsFrame::TimsFrame(%1,%2,nullptr,%3) FAILED")
              .arg(m_frameId)
              .arg(m_scanCount)
              .arg(len));
        }
    }
}
 
TimsFrame::TimsFrame(const TimsFrame &other) : TimsFrameBase(other)
{
}
 
TimsFrame::~TimsFrame()
{
}
 
 
void
TimsFrame::unshufflePacket(const char *src)
{
  // qDebug();
  quint64 len = m_binaryData.size();
  if(len % 4 != 0)
    {
      throw pappso::PappsoException(
        QObject::tr("TimsFrame::unshufflePacket error: len % 4 != 0"));
    }
 
  quint64 nb_uint4 = len / 4;
 
  char *dest         = m_binaryData.data();
  quint64 src_offset = 0;
 
  for(quint64 j = 0; j < 4; j++)
    {
      for(quint64 i = 0; i < nb_uint4; i++)
        {
          dest[(i * 4) + j] = src[src_offset];
          src_offset++;
        }
    }
  // qDebug();
}
 
std::size_t
TimsFrame::getScanPeakCount(std::size_t scanNum) const
{
  if(m_binaryData.size() == 0)
    return 0;
  /*
    if(scanNum == 0)
      {
        quint32 res = (*(quint32 *)(m_timsDataFrame.constData() + 4)) -
                      (*(quint32 *)(m_timsDataFrame.constData()-4));
        return res / 2;
      }*/
  if(scanNum == (m_scanCount - 1))
    {
      auto nb_uint4 = m_binaryData.size() / 4;
 
      std::size_t cumul = 0;
      for(quint32 i = 0; i < m_scanCount; i++)
        {
          cumul += (*(quint32 *)(m_binaryData.constData() + (i * 4)));
        }
      return (nb_uint4 - cumul) / 2;
    }
  checkScanNum(scanNum);
 
  // quint32 *res = (quint32 *)(m_timsDataFrame.constData() + (scanNum * 4));
  // qDebug() << " res=" << *res;
  return (*(quint32 *)(m_binaryData.constData() + ((scanNum + 1) * 4))) / 2;
}
 
std::size_t
TimsFrame::getScanOffset(std::size_t scanNum) const
{
  std::size_t offset = 0;
  for(std::size_t i = 0; i < (scanNum + 1); i++)
    {
      offset += (*(quint32 *)(m_binaryData.constData() + (i * 4)));
    }
  return offset;
}
 
 
std::vector<quint32>
TimsFrame::getScanTofIndexList(std::size_t scanNum) const
{
  // qDebug();
  checkScanNum(scanNum);
  std::vector<quint32> scan_tof;
 
  if(m_binaryData.size() == 0)
    return scan_tof;
  scan_tof.resize(getScanPeakCount(scanNum));
 
  std::size_t offset = getScanOffset(scanNum);
 
  qint32 previous = -1;
  for(std::size_t i = 0; i < scan_tof.size(); i++)
    {
      scan_tof[i] =
        (*(quint32 *)(m_binaryData.constData() + (offset * 4) + (i * 8))) +
        previous;
      previous = scan_tof[i];
    }
  // qDebug();
  return scan_tof;
}
 
std::vector<quint32>
TimsFrame::getScanIntensityList(std::size_t scanNum) const
{
  // qDebug();
  checkScanNum(scanNum);
  std::vector<quint32> scan_intensities;
 
  if(m_binaryData.size() == 0)
    return scan_intensities;
 
  scan_intensities.resize(getScanPeakCount(scanNum));
 
  std::size_t offset = getScanOffset(scanNum);
 
  for(std::size_t i = 0; i < scan_intensities.size(); i++)
    {
      scan_intensities[i] =
        (*(quint32 *)(m_binaryData.constData() + (offset * 4) + (i * 8) + 4));
    }
  // qDebug();
  return scan_intensities;
}
 
 
quint64
TimsFrame::cumulateScanIntensities(std::size_t scanNum) const
{
  // qDebug();
 
  quint64 summed_intensities = 0;
 
  if(m_binaryData.size() == 0)
    return summed_intensities;
  // checkScanNum(scanNum);
 
  std::size_t size = getScanPeakCount(scanNum);
 
  std::size_t offset = getScanOffset(scanNum);
 
  qint32 previous = -1;
 
  for(std::size_t i = 0; i < size; i++)
    {
      quint32 x =
        (*(quint32 *)((m_binaryData.constData() + (offset * 4) + (i * 8))) +
         previous);
 
      quint32 y =
        (*(quint32 *)(m_binaryData.constData() + (offset * 4) + (i * 8) + 4));
 
      previous = x;
 
      summed_intensities += y;
    }
 
  // Normalization over the accumulation time for this frame.
  summed_intensities *= ((double)100.0 / m_acqDurationInMilliseconds);
 
  // qDebug();
 
  return summed_intensities;
}
 
 
/**
 * @brief ...
 *
 * @param mobility_scan_begin p_mobility_scan_begin:...
 * @param mobility_scan_end p_mobility_scan_end:...
 * @return quint64
 */
quint64
TimsFrame::cumulateScanRangeIntensities(std::size_t mobility_scan_begin,
                                        std::size_t mobility_scan_end) const
{
  quint64 summed_intensities = 0;
 
  // qDebug() << "begin mobility_scan_begin =" << mobility_scan_begin
  //<< "mobility_scan_end =" << mobility_scan_end;
 
  if(m_binaryData.size() == 0)
    return summed_intensities;
 
  try
    {
      std::size_t mobility_scan_max = mobility_scan_end + 1;
 
      for(std::size_t i = mobility_scan_begin; i < mobility_scan_max; i++)
        {
          // qDebug() << i;
          summed_intensities += cumulateScanIntensities(i);
          // qDebug() << i;
        }
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure in %1 %2 to %3 :\n %4")
                    .arg(__FUNCTION__)
                    .arg(mobility_scan_begin)
                    .arg(mobility_scan_end)
                    .arg(error.what());
    }
 
  // qDebug() << "end";
 
  return summed_intensities;
}
 
 
void
TimsFrame::cumulateScan(std::size_t scanNum,
                        TimsDataFastMap &accumulate_into) const
{
  // qDebug();
 
  if(m_binaryData.size() == 0)
    return;
  // qDebug();
 
  // checkScanNum(scanNum);
 
  std::size_t scan_size = getScanPeakCount(scanNum);
 
  std::size_t scan_offset = getScanOffset(scanNum);
 
  qint32 previous = -1;
  for(std::size_t i = 0; i < scan_size; i++)
    {
      quint32 x = (*(quint32 *)((m_binaryData.constData() + (scan_offset * 4) +
                                 (i * 8))) +
                   previous);
      quint32 y = (*(quint32 *)(m_binaryData.constData() + (scan_offset * 4) +
                                (i * 8) + 4));
 
      previous = x;
      // qDebug() << "x=" << x << " y=" << y;
      accumulate_into.accumulateIntensity(x, y);
    }
 
  qDebug();
}
 
// Proof of concept m/z range-conditions for accumulations
void
TimsFrame::cumulateScan2(std::size_t scanNum,
                         TimsDataFastMap &accumulate_into,
                         quint32 accepted_tof_index_range_begin,
                         quint32 accepted_tof_index_range_end) const
{
  // qDebug();
 
  if(m_binaryData.size() == 0)
    return;
 
  // checkScanNum(scanNum);
 
  std::size_t scan_size   = getScanPeakCount(scanNum);
  std::size_t scan_offset = getScanOffset(scanNum);
  qint32 previous         = -1;
 
  for(std::size_t i = 0; i < scan_size; i++)
    {
      quint32 x = (*(quint32 *)((m_binaryData.constData() + (scan_offset * 4) +
                                 (i * 8))) +
                   previous);
 
      previous = x;
 
      if(x < accepted_tof_index_range_begin)
        {
 
          // qDebug() << "TOF index still not in range, x:" << x;
          continue;
        }
      if(x > accepted_tof_index_range_end)
        {
          // qDebug() << "TOF index already out of range, x:" << x;
          break;
        }
 
      // qDebug() << "TOF index in range, x:" << x;
 
      quint32 y = (*(quint32 *)(m_binaryData.constData() + (scan_offset * 4) +
                                (i * 8) + 4));
 
      accumulate_into.accumulateIntensity(x, y);
    }
 
  // qDebug();
}
 
 
Trace
TimsFrame::cumulateScansToTrace(std::size_t mobility_scan_begin,
                                std::size_t mobility_scan_end) const
{
  // qDebug();
 
  Trace new_trace;
 
  try
    {
      if(m_binaryData.size() == 0)
        return new_trace;
      TimsDataFastMap &raw_spectrum =
        TimsDataFastMap::getTimsDataFastMapInstance();
      raw_spectrum.clear();
      // double local_accumulationTime = 0;
 
      std::size_t mobility_scan_max = mobility_scan_end + 1;
      qDebug();
      for(std::size_t i = mobility_scan_begin; i < mobility_scan_max; i++)
        {
          // qDebug() << i;
          cumulateScan(i, raw_spectrum);
          // qDebug() << i;
 
          // local_accumulationTime += m_accumulationTime;
        }
 
      // qDebug();
 
      pappso::DataPoint data_point_cumul;
 
 
      MzCalibrationInterface *mz_calibration_p =
        getMzCalibrationInterfaceSPtr().get();
 
 
      for(quint32 tof_index : raw_spectrum.getTofIndexList())
        {
          data_point_cumul.x = mz_calibration_p->getMzFromTofIndex(tof_index);
          // normalization
          data_point_cumul.y = raw_spectrum.readIntensity(tof_index) *
                               ((double)100.0 / m_acqDurationInMilliseconds);
          new_trace.push_back(data_point_cumul);
        }
      new_trace.sortX();
 
      // qDebug();
    }
 
  catch(std::exception &error)
    {
      qDebug() << QString(
                    "Failure in TimsFrame::cumulateScanToTrace %1 to %2 :\n %3")
                    .arg(mobility_scan_begin, mobility_scan_end)
                    .arg(error.what());
    }
  return new_trace;
}
 
Trace
TimsFrame::combineScansToTraceWithDowngradedMzResolution(
  std::size_t mz_index_merge_window,
  std::size_t mobility_scan_begin,
  std::size_t mobility_scan_end,
  quint32 &mz_minimum_index_out,
  quint32 &mz_maximum_index_out) const
{
  // qDebug();
 
  Trace new_trace;
 
  try
    {
      if(m_binaryData.size() == 0)
        {
          qDebug() << "The frame is empty, returning empty trace.";
          return new_trace;
        }
 
      // Allocate a map for (TOF,intensity) pairs to
      // accumulate ion mobility scans.
 
      TimsDataFastMap &raw_spectrum =
        TimsDataFastMap::getTimsDataFastMapInstance();
      raw_spectrum.clear();
      // double local_accumulationTime = 0;
 
      std::size_t mobility_scan_max = mobility_scan_end + 1;
 
      for(std::size_t i = mobility_scan_begin; i < mobility_scan_max; i++)
        {
          // qDebug() << "Going to cumulate currently iterated mobility scan:"
          // << i;
          cumulateScan(i, raw_spectrum);
          // qDebug() << "Done cumulating currently iterated mobility scan:" <<
          // i;
 
          // local_accumulationTime += m_accumulationTime;
        }
 
      // qDebug();
 
      pappso::DataPoint data_point_cumul;
 
      MzCalibrationInterface *mz_calibration_p =
        getMzCalibrationInterfaceSPtr().get();
 
      // If the caller asks that m/z values be binned larger than they are,
      // ask that the m/z raw map be reduced in resolution.
      if(mz_index_merge_window > 0)
        {
          raw_spectrum.downsizeMzRawMap(mz_index_merge_window);
        }
 
      // Store the first mz index and the last mz index of the current spectrum.
      // The values are set to the out parameters.
      mz_minimum_index_out = std::numeric_limits<quint32>::max();
      mz_maximum_index_out = 0;
 
      // FIXME: donc je comprends que les index sont contigus
      // puisqu'on utilise at(key) ?
      for(quint32 tof_index : raw_spectrum.getTofIndexList())
        {
          if(tof_index > mz_maximum_index_out)
            mz_maximum_index_out = tof_index;
          if(tof_index < mz_minimum_index_out)
            mz_minimum_index_out = tof_index;
 
          // Convert the TOF index to m/z
          data_point_cumul.x = mz_calibration_p->getMzFromTofIndex(tof_index);
 
          // Normalization
          data_point_cumul.y = raw_spectrum.readIntensity(tof_index) *
                               ((double)100.0 / m_acqDurationInMilliseconds);
 
          // Finally make the data point a new Trace point.
          new_trace.push_back(data_point_cumul);
        }
 
      // qDebug() << "At this point we have mz_minimum_index_out:"
      //          << mz_minimum_index_out
      //          << "and mz_maximum_index_out:" << mz_maximum_index_out;
 
      // FIXME: this does not seem to be necessary since raw_spectrum is a map
      // with auto-sorting on the keys which are quint32.
      // new_trace.sortX();
 
      // qDebug();
    }
  catch(std::exception &error)
    {
      qDebug() << QString(
                    "Failure in TimsFrame::cumulateScanToTrace %1 to %2 :\n %3")
                    .arg(mobility_scan_begin, mobility_scan_end)
                    .arg(error.what());
    }
 
  // qDebug() << "Returning new trace of size:" << new_trace.size();
 
  return new_trace;
}
 
 
Trace
TimsFrame::combineScansToTraceWithDowngradedMzResolution2(
  std::size_t mz_index_merge_window,
  double mz_range_begin,
  double mz_range_end,
  std::size_t mobility_scan_begin,
  std::size_t mobility_scan_end,
  quint32 &mz_minimum_index_out,
  quint32 &mz_maximum_index_out) const
{
  // qDebug() << "Calling cumulateScansToTraceMzDownResolution2 for both mz and "
  //             "im ranges accounting.";
 
  Trace new_trace;
 
  try
    {
      if(m_binaryData.size() == 0)
        {
          qDebug() << "The frame is empty, returning empty trace.";
          return new_trace;
        }
 
      // Allocate a map for (TOF index,intensity) pairs to
      // accumulate ion mobility scans.
 
      TimsDataFastMap &raw_spectrum =
        TimsDataFastMap::getTimsDataFastMapInstance();
      raw_spectrum.clear();
      // double local_accumulationTime = 0;
 
      std::size_t mobility_scan_max = mobility_scan_end + 1;
 
      quint32 tof_index_for_mz_range_begin =
        getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(
          mz_range_begin);
      quint32 tof_index_for_mz_range_end =
        getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(mz_range_end);
 
      // qDebug() << "TOF index for mz range begin:"
      //          << tof_index_for_mz_range_begin;
      // qDebug() << "TOF index for mz range end:"
      //          << tof_index_for_mz_range_end;
 
      for(std::size_t iter = mobility_scan_begin; iter < mobility_scan_max;
          iter++)
        {
          // qDebug() << "Going to cumulate currently iterated mobility scan:"
          // << iter;
          cumulateScan2(iter,
                        raw_spectrum,
                        tof_index_for_mz_range_begin,
                        tof_index_for_mz_range_end);
          // qDebug() << "Done cumulating currently iterated mobility scan:" <<
          // i;
 
          // local_accumulationTime += m_accumulationTime;
        }
 
      // qDebug();
 
      pappso::DataPoint data_point_cumul;
 
      MzCalibrationInterface *mz_calibration_p =
        getMzCalibrationInterfaceSPtr().get();
 
      // If the caller asks that m/z values be binned larger than they are,
      // ask that the m/z raw map be reduced in resolution.
      if(mz_index_merge_window > 0)
        {
          raw_spectrum.downsizeMzRawMap(mz_index_merge_window);
        }
 
      // Store the first mz index and the last mz index of the current spectrum.
      // The values are set to the out parameters.
      mz_minimum_index_out = std::numeric_limits<quint32>::max();
      mz_maximum_index_out = 0;
 
      // for(std::pair<quint32, quint32> pair_tof_intensity : raw_spectrum)
      for(quint32 tof_index : raw_spectrum.getTofIndexList())
        {
          std::size_t intensity = raw_spectrum.readIntensity(tof_index);
          if(tof_index > mz_maximum_index_out)
            mz_maximum_index_out = tof_index;
          if(tof_index < mz_minimum_index_out)
            mz_minimum_index_out = tof_index;
 
          // Convert the TOF index to m/z
          data_point_cumul.x = mz_calibration_p->getMzFromTofIndex(tof_index);
 
          // Normalization
          data_point_cumul.y =
            intensity * ((double)100.0 / m_acqDurationInMilliseconds);
 
          // Finally make the data point a new Trace point.
          new_trace.push_back(data_point_cumul);
        }
 
      // qDebug() << "At this point we have mz_minimum_index_out:"
      //          << mz_minimum_index_out
      //          << "and mz_maximum_index_out:" << mz_maximum_index_out;
 
      // FIXME: this does not seem to be necessary since raw_spectrum is a map
      // with auto-sorting on the keys which are quint32.
      // new_trace.sortX();
 
      // qDebug();
    }
  catch(std::exception &error)
    {
      qDebug() << QString(
                    "Failure in TimsFrame::cumulateScanToTrace %1 to %2 :\n %3")
                    .arg(mobility_scan_begin, mobility_scan_end)
                    .arg(error.what());
    }
 
  // qDebug() << "Returning new trace of size:" << new_trace.size();
 
  return new_trace;
}
 
 
void
TimsFrame::combineScansInTofIndexIntensityMap(TimsDataFastMap &rawSpectrum,
                                              std::size_t scan_index_begin,
                                              std::size_t scan_index_end) const
{
  // qDebug() << "begin mobility_scan_begin=" << mobility_scan_begin
  //<< " mobility_scan_end=" << mobility_scan_end;
 
  if(m_binaryData.size() == 0)
    return;
  try
    {
 
      std::size_t mobility_scan_max = scan_index_end + 1;
      // if (mobility_scan_max > m_scanCount) mobility_scan_max = m_scanCount;
      // qDebug();
      for(std::size_t i = scan_index_begin; i < mobility_scan_max; i++)
        {
          // qDebug() << i;
          cumulateScan(i, rawSpectrum);
          // qDebug() << i;
 
          // local_accumulationTime += m_accumulationTime;
        }
    }
 
  catch(std::exception &error)
    {
      qDebug() << QString("Failure in %1 %2 to %3 :\n %4")
                    .arg(__FUNCTION__)
                    .arg(scan_index_begin)
                    .arg(scan_index_end)
                    .arg(error.what());
    }
 
  // qDebug() << "end";
}
 
 
void
TimsFrame::combineScansInTofIndexIntensityMap(TimsDataFastMap &rawSpectrum,
                                              std::size_t scan_index_begin,
                                              std::size_t scan_index_end,
                                              quint32 tof_index_begin,
                                              quint32 tof_index_end) const
{
  // qDebug() << "tof_index_begin=" << tof_index_begin
  //          << " tof_index_end=" << tof_index_end;
 
  //<< " mobility_scan_end=" << mobility_scan_end;
 
  if(m_binaryData.size() == 0)
    return;
  try
    {
 
      std::size_t mobility_scan_max = scan_index_end + 1;
      // if (mobility_scan_max > m_scanCount) mobility_scan_max = m_scanCount;
      // qDebug();
      for(std::size_t i = scan_index_begin; i < mobility_scan_max; i++)
        {
          // qDebug() << i;
          cumulateScan2(i, rawSpectrum, tof_index_begin, tof_index_end);
          // qDebug() << i << rawSpectrum.getTofIndexList().size();
          // local_accumulationTime += m_accumulationTime;
        }
    }
 
  catch(std::exception &error)
    {
      qDebug() << QString("Failure in %1 %2 to %3 :\n %4")
                    .arg(__FUNCTION__)
                    .arg(scan_index_begin)
                    .arg(scan_index_end)
                    .arg(error.what());
    }
 
  // qDebug() << "end";
}
 
pappso::MassSpectrumSPtr
TimsFrame::getMassSpectrumSPtr(std::size_t scanNum) const
{
 
  // qDebug() << " scanNum=" << scanNum;
 
  checkScanNum(scanNum);
 
  // qDebug();
 
  pappso::MassSpectrumSPtr mass_spectrum_sptr =
    std::make_shared<pappso::MassSpectrum>();
  // std::vector<DataPoint>
 
  if(m_binaryData.size() == 0)
    return mass_spectrum_sptr;
 
  // qDebug();
 
  std::size_t size = getScanPeakCount(scanNum);
 
  std::size_t offset = getScanOffset(scanNum);
 
  MzCalibrationInterface *mz_calibration_p =
    getMzCalibrationInterfaceSPtr().get();
 
 
  qint32 previous = -1;
  qint32 tof_index;
  // std::vector<quint32> index_list;
  DataPoint data_point;
  for(std::size_t i = 0; i < size; i++)
    {
      tof_index =
        (*(quint32 *)((m_binaryData.constData() + (offset * 4) + (i * 8))) +
         previous);
      data_point.y =
        (*(quint32 *)(m_binaryData.constData() + (offset * 4) + (i * 8) + 4));
 
      // intensity normalization
      data_point.y *= 100.0 / m_acqDurationInMilliseconds;
 
      previous = tof_index;
 
 
      // mz calibration
      data_point.x = mz_calibration_p->getMzFromTofIndex(tof_index);
      mass_spectrum_sptr.get()->push_back(data_point);
    }
 
  // qDebug();
 
  return mass_spectrum_sptr;
}
 
 
Trace
TimsFrame::getMobilityScan(std::size_t scanNum,
                           std::size_t mz_index_merge_window,
                           double mz_range_begin,
                           double mz_range_end,
                           quint32 &mz_minimum_index_out,
                           quint32 &mz_maximum_index_out) const
{
  // qDebug() << "mz_range_begin:" << mz_range_begin
  //          << "mz_range_end:" << mz_range_end
  //          << "mz_index_merge_window:" << mz_index_merge_window;
 
  Trace spectrum;
 
  quint32 mz_index_begin = 0;
  quint32 mz_index_end   = std::numeric_limits<quint32>::max();
 
  if(mz_range_end > 0)
    {
      // qDebug() << "m/z range is requested.";
 
      mz_index_begin =
        msp_mzCalibration.get()->getTofIndexFromMz(mz_range_begin);
      mz_index_end = msp_mzCalibration.get()->getTofIndexFromMz(mz_range_end);
    }
 
  // qDebug() << "After conversion of mz indices, mz_index_begin:"
  //          << mz_index_begin << "mz_index_end;" << mz_index_end;
 
  auto raw_spectrum =
    getRawValuePairList(scanNum, mz_index_begin, mz_index_end);
 
  // qDebug() << " raw_spectrum.size();" << raw_spectrum.size();
 
  if(mz_index_merge_window > 0)
    {
      // qDebug() << "mz_index_merge_window;=" << mz_index_merge_window
      //          << " raw_spectrum.size()=" << raw_spectrum.size();
      raw_spectrum = downgradeResolutionOfTofIndexIntensityPairList(
        mz_index_merge_window, raw_spectrum);
    }
 
  if(raw_spectrum.size() > 0)
    {
      mz_minimum_index_out = raw_spectrum.front().tof_index;
      mz_maximum_index_out = raw_spectrum.back().tof_index;
 
      for(auto &&element : raw_spectrum)
        {
          spectrum.push_back(DataPoint(
            msp_mzCalibration.get()->getMzFromTofIndex(element.tof_index),
            static_cast<double>(element.intensity_index)));
 
          // intensity normalization
          spectrum.back().y *= 100.0 / m_acqDurationInMilliseconds;
        }
    }
 
  return spectrum;
}
 
void
TimsFrame::extractTimsXicListInRtRange(
  std::vector<XicCoordTims *>::iterator &itXicListbegin,
  std::vector<XicCoordTims *>::iterator &itXicListend,
  Enums::XicExtractMethod method) const
{
  // qDebug() << std::distance(itXicListbegin, itXicListend);
 
  std::vector<TimsFrame::XicComputeStructure> tmp_xic_list;
 
  for(auto it = itXicListbegin; it != itXicListend; it++)
    {
      tmp_xic_list.push_back(TimsFrame::XicComputeStructure(this, **it));
 
      // qDebug() << " tmp_xic_struct.mobilityIndexBegin="
      //          << tmp_xic_list.back().mobilityIndexBegin
      //          << " tmp_xic_struct.mobilityIndexEnd="
      //          << tmp_xic_list.back().mobilityIndexEnd;
 
      // qDebug() << " tmp_xic_struct.mzIndexLowerBound="
      //          << tmp_xic_list.back().mzIndexLowerBound
      //          << " tmp_xic_struct.mzIndexUpperBound="
      //          << tmp_xic_list.back().mzIndexUpperBound;
    }
  if(tmp_xic_list.size() == 0)
    return;
  /*
  std::sort(tmp_xic_list.begin(), tmp_xic_list.end(), [](const
  TimsXicStructure &a, const TimsXicStructure &b) { return
  a.mobilityIndexBegin < b.mobilityIndexBegin;
            });
            */
  std::vector<std::size_t> unique_scan_num_list;
  for(auto &&struct_xic : tmp_xic_list)
    {
      for(std::size_t scan = struct_xic.mobilityIndexBegin;
          (scan <= struct_xic.mobilityIndexEnd) && (scan < m_scanCount);
          scan++)
        {
          unique_scan_num_list.push_back(scan);
        }
    }
  std::sort(unique_scan_num_list.begin(), unique_scan_num_list.end());
  auto it_scan_num_end =
    std::unique(unique_scan_num_list.begin(), unique_scan_num_list.end());
  auto it_scan_num = unique_scan_num_list.begin();
 
  while(it_scan_num != it_scan_num_end)
    {
      TraceSPtr ms_spectrum = getRawTraceSPtr(*it_scan_num);
      // qDebug() << ms_spectrum.get()->toString();
      for(auto &&tmp_xic_struct : tmp_xic_list)
        {
          if(((*it_scan_num) >= tmp_xic_struct.mobilityIndexBegin) &&
             ((*it_scan_num) <= tmp_xic_struct.mobilityIndexEnd))
            {
              if(method == Enums::XicExtractMethod::max)
                {
                  tmp_xic_struct.tmpIntensity +=
                    ms_spectrum.get()->maxY(tmp_xic_struct.mzIndexLowerBound,
                                            tmp_xic_struct.mzIndexUpperBound);
 
                  // qDebug() << "tmp_xic_struct.tmpIntensity="
                  //          << tmp_xic_struct.tmpIntensity;
                }
              else
                {
                  // sum
                  tmp_xic_struct.tmpIntensity +=
                    ms_spectrum.get()->sumY(tmp_xic_struct.mzIndexLowerBound,
                                            tmp_xic_struct.mzIndexUpperBound);
 
                  // qDebug() << "tmp_xic_struct.tmpIntensity="
                  //          << tmp_xic_struct.tmpIntensity;
                }
            }
        }
      it_scan_num++;
    }
 
  for(auto &&tmp_xic_struct : tmp_xic_list)
    {
      if(tmp_xic_struct.tmpIntensity != 0)
        {
          // qDebug() << tmp_xic_struct.xic_ptr;
 
          tmp_xic_struct.xic_ptr->push_back(
            {m_rtInSeconds, tmp_xic_struct.tmpIntensity});
        }
    }
 
  // qDebug();
}
 
 
pappso::TraceSPtr
TimsFrame::getRawTraceSPtr(std::size_t scanNum) const
{
 
  // qDebug();
 
  pappso::TraceSPtr trace_sptr = std::make_shared<pappso::Trace>();
  // std::vector<DataPoint>
 
  if(m_binaryData.size() == 0)
    return trace_sptr;
  // qDebug();
 
  std::size_t size = getScanPeakCount(scanNum);
 
  std::size_t offset = getScanOffset(scanNum);
 
  qint32 previous = -1;
  std::vector<quint32> index_list;
  for(std::size_t i = 0; i < size; i++)
    {
      DataPoint data_point(
        (*(quint32 *)((m_binaryData.constData() + (offset * 4) + (i * 8))) +
         previous),
        (*(quint32 *)(m_binaryData.constData() + (offset * 4) + (i * 8) + 4)));
 
      // intensity normalization
      data_point.y *= 100.0 / m_acqDurationInMilliseconds;
 
      previous = data_point.x;
      trace_sptr.get()->push_back(data_point);
    }
  // qDebug();
  return trace_sptr;
}
 
 
std::vector<TimsFrame::TofIndexIntensityPair>
TimsFrame::getRawValuePairList(std::size_t scanNum,
                               quint32 accepted_tof_index_range_begin,
                               quint32 accepted_tof_index_range_end) const
{
  // qDebug() << accepted_tof_index_range_begin;
 
  std::vector<TimsFrame::TofIndexIntensityPair> raw_value_pairs;
  // std::vector<DataPoint>
 
  if(m_binaryData.size() == 0)
    return raw_value_pairs;
  // qDebug();
 
  std::size_t size = getScanPeakCount(scanNum);
 
  std::size_t offset = getScanOffset(scanNum);
 
  // qDebug() << size;
  qint32 previous = -1;
  std::vector<quint32> index_list;
  for(std::size_t i = 0; i < size; i++)
    {
 
      // qDebug() << i;
      TimsFrame::TofIndexIntensityPair raw_value_pair(
        {(*(quint32 *)((m_binaryData.constData() + (offset * 4) + (i * 8))) +
          previous),
         (*(quint32 *)(m_binaryData.constData() + (offset * 4) + (i * 8) +
                       4))});
 
      previous = raw_value_pair.tof_index;
      if(raw_value_pair.tof_index < accepted_tof_index_range_begin)
        {
          // qDebug() << "TOF index still not in range, x:" << x;
          continue;
        }
      if(raw_value_pair.tof_index > accepted_tof_index_range_end)
        {
          // qDebug() << "TOF index already out of range, x:" << x;
          break;
        }
 
      raw_value_pairs.push_back(raw_value_pair);
    }
  // qDebug() << raw_value_pairs.size();
  return raw_value_pairs;
}
 
} // namespace pappso