timsddaprecursors.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/timsddaprecursors.h
 * \date 30/06/2024
 * \brief handle specific data for DDA MS runs
 */
 
/*******************************************************************************
 * Copyright (c) 2024 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 "timsddaprecursors.h"
#include "pappsomspp/core/exception/exceptionnotfound.h"
#include "pappsomspp/core/exception/exceptioninterrupted.h"
#include <QSqlError>
#include "pappsomspp/core/processing/combiners/tracepluscombiner.h"
#include "pappsomspp/core/processing/filters/filtertriangle.h"
#include "pappsomspp/core/processing/filters/filtersuitestring.h"
#include <QtConcurrent>
namespace pappso
{
 
TimsDdaPrecursors::TimsDdaPrecursors(QSqlQuery &query, TimsData *tims_data_origin)
  : mp_timsDataOrigin(tims_data_origin)
{
 
  // get number of precursors
  qDebug();
  m_totalPrecursorCount = 0;
  if(!query.exec("SELECT COUNT( DISTINCT Id) FROM Precursors;"))
    {
      qDebug();
      throw pappso::ExceptionNotFound(QObject::tr("ERROR : no Precursors in SqlLite database"));
    }
  else
    {
      if(query.next())
        {
          m_totalPrecursorCount = query.value(0).toLongLong();
        }
    }
 
  qDebug();
 
  mcsp_ms2Filter =
    std::make_shared<FilterSuiteString>("chargeDeconvolution|0.02dalton mzExclusion|0.01dalton");
 
 
  std::shared_ptr<FilterTriangle> ms1filter = std::make_shared<FilterTriangle>();
  ms1filter.get()->setTriangleSlope(50, 0.01);
  mcsp_ms1Filter = ms1filter;
}
 
TimsDdaPrecursors::~TimsDdaPrecursors()
{
}
 
std::vector<TimsDdaPrecursors::SpectrumDescr>
TimsDdaPrecursors::getSpectrumDescrListByFrameId(std::size_t frame_id) const
{
  std::vector<TimsDdaPrecursors::SpectrumDescr> spectrum_descr_list;
 
  try
    {
      // QMutexLocker lock(&m_mutex);
      //  Go get records!
 
 
      QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
 
      QSqlQuery q(qdb);
      q.exec(QString("SELECT PasefFrameMsMsInfo.Frame, "    // 0
                     "PasefFrameMsMsInfo.ScanNumBegin, "    // 1
                     "PasefFrameMsMsInfo.ScanNumEnd, "      // 2
                     "PasefFrameMsMsInfo.IsolationMz, "     // 3
                     "PasefFrameMsMsInfo.IsolationWidth, "  // 4
                     "PasefFrameMsMsInfo.CollisionEnergy, " // 5
                     "PasefFrameMsMsInfo.Precursor, "       // 6
                     "Precursors.Id, "                      // 7
                     "Precursors.LargestPeakMz, "           // 8
                     "Precursors.AverageMz, "               // 9
                     "Precursors.MonoisotopicMz, "          // 10
                     "Precursors.Charge, "                  // 11
                     "Precursors.ScanNumber, "              // 12
                     "Precursors.Intensity, "               // 13
                     "Precursors.Parent "                   // 14
                     "FROM PasefFrameMsMsInfo "
                     "INNER JOIN Precursors ON "
 
 
                     "PasefFrameMsMsInfo.Precursor=Precursors.Id "
                     "WHERE PasefFrameMsMsInfo.Frame=%1;")
               .arg(frame_id));
      if(q.lastError().isValid())
        {
          qDebug();
          throw PappsoException(
            QObject::tr("ERROR in TIMS sqlite database file %1, executing SQL "
                        "command %2:\n%3\n%4\n%5")
              .arg(mp_timsDataOrigin->m_timsDataDirectory.absoluteFilePath("analysis.tdf"))
              .arg(q.lastQuery())
              .arg(qdb.lastError().databaseText())
              .arg(qdb.lastError().driverText())
              .arg(qdb.lastError().nativeErrorCode()));
        }
 
      q.last(); // strange bug : get the last sql record and get back,
      // otherwise it will not retrieve all records.
      q.first();
      // std::size_t i = 0;
      do
        {
          TimsDdaPrecursors::SpectrumDescr spectrum_descr;
          spectrum_descr.tims_frame_list.clear();
          // mass_spectrum.setPrecursorCharge(q.value(11).toInt());
          // mass_spectrum.setPrecursorMz(q.value(10).toDouble());
          // mass_spectrum.setPrecursorIntensity(q.value(13).toDouble());
          spectrum_descr.precursor_ion_data =
            PrecursorIonData(q.value(10).toDouble(), q.value(11).toInt(), q.value(13).toDouble());
 
          spectrum_descr.precursor_id = q.value(6).toLongLong();
          spectrum_descr.ms2_index    = (spectrum_descr.precursor_id * 2) - 1;
          spectrum_descr.ms1_index    = (spectrum_descr.precursor_id * 2) - 2;
 
          spectrum_descr.scan_mobility_start = q.value(1).toLongLong();
          spectrum_descr.scan_mobility_end   = q.value(2).toLongLong();
 
          spectrum_descr.isolationMz     = q.value(3).toDouble();
          spectrum_descr.isolationWidth  = q.value(4).toDouble();
          spectrum_descr.collisionEnergy = q.value(5).toFloat();
          spectrum_descr.parent_frame    = q.value(14).toLongLong();
 
          spectrum_descr_list.push_back(spectrum_descr);
        }
      while(q.next());
    }
  catch(PappsoException &error)
    {
      throw error;
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
  return spectrum_descr_list;
}
 
TimsDataFastMap &
TimsDdaPrecursors::getCombinedMs2ScansByPrecursorId(std::size_t precursor_id)
{
 
  qDebug();
  TimsDataFastMap &raw_spectrum = TimsDataFastMap::getTimsDataFastMapInstance();
  raw_spectrum.clear();
  try
    {
      QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
 
      QSqlQuery q(qdb);
      q.exec(QString("SELECT PasefFrameMsMsInfo.*, Precursors.* FROM "
                     "PasefFrameMsMsInfo INNER JOIN Precursors ON "
                     "PasefFrameMsMsInfo.Precursor=Precursors.Id where "
                     "Precursors.Id=%1;")
               .arg(precursor_id));
      if(q.lastError().isValid())
        {
          qDebug();
          throw PappsoException(
            QObject::tr("ERROR in TIMS sqlite database file %1, executing SQL "
                        "command %2:\n%3\n%4\n%5")
              .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf"))
              .arg(q.lastQuery())
              .arg(qdb.lastError().databaseText())
              .arg(qdb.lastError().driverText())
              .arg(qdb.lastError().nativeErrorCode()));
        }
      qDebug();
      // m_mutex.unlock();
      if(q.size() == 0)
        {
 
          throw ExceptionNotFound(
            QObject::tr("ERROR in getQualifiedMassSpectrumByPrecursorId, precursor "
                        "id=%1 not found")
              .arg(precursor_id));
        }
      else
        {
          //  qDebug() << " q.size()="<< q.size();
          qDebug();
          bool first                      = true;
          std::size_t scan_mobility_start = 0;
          std::size_t scan_mobility_end   = 0;
          std::vector<std::size_t> tims_frame_list;
 
          while(q.next())
            {
              tims_frame_list.push_back(q.value(0).toLongLong());
              if(first)
                {
 
                  scan_mobility_start = q.value(1).toLongLong();
                  scan_mobility_end   = q.value(2).toLongLong();
 
                  first = false;
                }
            }
          // QMutexLocker locker(&m_mutex_spectrum);
          qDebug();
          TimsFrameCstSPtr tims_frame, previous_frame;
          // TracePlusCombiner combiner;
          // MapTrace combiner_result;
          for(std::size_t tims_id : tims_frame_list)
            {
              tims_frame = mp_timsDataOrigin->getTimsFrameCstSPtrCached(tims_id);
              qDebug();
              /*combiner.combine(combiner_result,
                tims_frame.get()->cumulateScanToTrace(
                scan_mobility_start, scan_mobility_end));*/
              if(previous_frame.get() != nullptr)
                {
                  if(previous_frame.get()->hasSameCalibrationData(*tims_frame.get()))
                    {
                    }
                  else
                    {
                      throw ExceptionNotFound(
                        QObject::tr("ERROR in %1 %2, different calibration data "
                                    "between frame id %3 and frame id %4")
                          .arg(__FILE__)
                          .arg(__FUNCTION__)
                          .arg(previous_frame.get()->getId())
                          .arg(tims_frame.get()->getId()));
                    }
                }
              tims_frame.get()->combineScansInTofIndexIntensityMap(
                raw_spectrum, scan_mobility_start, scan_mobility_end);
              qDebug();
 
              previous_frame = tims_frame;
            }
          qDebug() << " precursor_index=" << precursor_id << " num_rows=" << tims_frame_list.size()
                   << " sql=" << q.lastQuery() << " " << (std::size_t)QThread::currentThreadId();
          if(first == true)
            {
              throw ExceptionNotFound(
                QObject::tr("ERROR in getQualifiedMassSpectrumByPrecursorId, precursor "
                            "id=%1 not found")
                  .arg(precursor_id));
            }
          qDebug();
        }
    }
 
  catch(PappsoException &error)
    {
      throw PappsoException(QObject::tr("ERROR in %1 (precursor_index=%2):\n%3")
                              .arg(__FUNCTION__)
                              .arg(precursor_id)
                              .arg(error.qwhat()));
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
  return raw_spectrum;
  qDebug();
}
 
std::size_t
TimsDdaPrecursors::getTotalPrecursorCount() const
{
  return m_totalPrecursorCount;
}
 
XicCoordTims
TimsDdaPrecursors::getXicCoordTimsFromPrecursorId(std::size_t precursor_id,
                                                  PrecisionPtr precision_ptr)
{
 
  qDebug();
  XicCoordTims xic_coord_tims_struct;
 
  try
    {
      if(m_mapXicCoordRecord.size() == 0)
        {
          QMutexLocker lock(&m_mutex);
          // Go get records!
 
          // We proceed in this way:
 
          // 1. For each Precursor reference to the Precursors table's ID
          // found in the PasefFrameMsMsInfo table, store the precursor ID for
          // step 2.
 
          // 2. From the Precursors table's ID from step 1, get the
          // MonoisotopicMz.
 
          // 3. From the PasefFrameMsMsInfo table, for the Precursors table's
          // ID reference, get a reference to the Frames table's ID. Thanks to
          // the Frames ID, look for the Time value (acquisition retention
          // time) for the MS/MS spectrum. The Time value in the Frames tables
          // always corresponds to a Frame of MsMsType 8 (that is, MS/MS),
          // which is expected since we are looking into MS/MS data.
 
          // 4. From the PasefFrameMsMsInfo table, associate the values
          // ScanNumBegin and ScanNumEnd, the mobility bins in which the
          // precursor was found.
 
 
          QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
          QSqlQuery q(qdb);
          q.exec(
            QString("SELECT Precursors.id, "
                    "min(Frames.Time), "
                    "min(PasefFrameMsMsInfo.ScanNumBegin), "
                    "max(PasefFrameMsMsInfo.ScanNumEnd), "
                    "Precursors.MonoisotopicMz "
                    "FROM "
                    "PasefFrameMsMsInfo INNER JOIN Precursors ON "
                    "PasefFrameMsMsInfo.Precursor=Precursors.Id "
                    "INNER JOIN Frames ON "
                    "PasefFrameMsMsInfo.Frame=Frames.Id "
                    "GROUP BY Precursors.id;"));
          if(q.lastError().isValid())
            {
              qDebug();
              throw PappsoException(
                QObject::tr("ERROR in TIMS sqlite database file %1, executing SQL "
                            "command %2:\n%3\n%4\n%5")
                  .arg(mp_timsDataOrigin->m_timsDataDirectory.absoluteFilePath("analysis.tdf"))
                  .arg(q.lastQuery())
                  .arg(qdb.lastError().databaseText())
                  .arg(qdb.lastError().driverText())
                  .arg(qdb.lastError().nativeErrorCode()));
            }
 
          q.last(); // strange bug : get the last sql record and get back,
          // otherwise it will not retrieve all records.
          q.first();
          // std::size_t i = 0;
          do
            {
              QSqlRecord record = q.record();
              m_mapXicCoordRecord.insert(std::pair<std::size_t, QSqlRecord>(
                (std::size_t)record.value(0).toULongLong(), record));
            }
          while(q.next());
        }
 
 
      auto it_map_xiccoord = m_mapXicCoordRecord.find(precursor_id);
      if(it_map_xiccoord == m_mapXicCoordRecord.end())
        {
 
          throw ExceptionNotFound(
            QObject::tr("ERROR Precursors database id %1 not found").arg(precursor_id));
        }
 
      auto &q                            = it_map_xiccoord->second;
      xic_coord_tims_struct.mzRange      = MzRange(q.value(4).toDouble(), precision_ptr);
      xic_coord_tims_struct.scanNumBegin = q.value(2).toUInt();
      xic_coord_tims_struct.scanNumEnd   = q.value(3).toUInt();
      xic_coord_tims_struct.rtTarget     = q.value(1).toDouble();
      // xic_structure.charge       = q.value(5).toUInt();
      xic_coord_tims_struct.xicSptr = std::make_shared<Xic>();
    }
  catch(PappsoException &error)
    {
      throw error;
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
  return xic_coord_tims_struct;
}
 
 
TimsDdaPrecursors::SpectrumDescr
TimsDdaPrecursors::getSpectrumDescrWithPrecursorId(std::size_t precursor_id) const
{
 
  SpectrumDescr spectrum_descr;
  QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
  QSqlQuery q(qdb);
  q.exec(QString("SELECT PasefFrameMsMsInfo.Frame, "    // 0
                 "PasefFrameMsMsInfo.ScanNumBegin, "    // 1
                 "PasefFrameMsMsInfo.ScanNumEnd, "      // 2
                 "PasefFrameMsMsInfo.IsolationMz, "     // 3
                 "PasefFrameMsMsInfo.IsolationWidth, "  // 4
                 "PasefFrameMsMsInfo.CollisionEnergy, " // 5
                 "PasefFrameMsMsInfo.Precursor, "       // 6
                 "Precursors.Id, "                      // 7
                 "Precursors.LargestPeakMz, "           // 8
                 "Precursors.AverageMz, "               // 9
                 "Precursors.MonoisotopicMz, "          // 10
                 "Precursors.Charge, "                  // 11
                 "Precursors.ScanNumber, "              // 12
                 "Precursors.Intensity, "               // 13
                 "Precursors.Parent "                   // 14
                 "FROM PasefFrameMsMsInfo "
                 "INNER JOIN Precursors ON "
 
 
                 "PasefFrameMsMsInfo.Precursor=Precursors.Id "
                 "WHERE Precursors.Id=%1;")
           .arg(precursor_id));
  if(q.lastError().isValid())
    {
 
      throw PappsoException(
        QObject::tr("ERROR in TIMS sqlite database file %1, executing SQL "
                    "command %2:\n%3\n%4\n%5")
          .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf"))
          .arg(q.lastQuery())
          .arg(qdb.lastError().databaseText())
          .arg(qdb.lastError().driverText())
          .arg(qdb.lastError().nativeErrorCode()));
    }
 
 
  bool first = true;
  while(q.next())
    {
 
      qDebug() << " cumul tims frame:" << q.value(0).toLongLong();
      spectrum_descr.tims_frame_list.push_back(q.value(0).toLongLong());
      if(first)
        {
          // mass_spectrum.setPrecursorCharge(q.value(11).toInt());
          // mass_spectrum.setPrecursorMz(q.value(10).toDouble());
          // mass_spectrum.setPrecursorIntensity(q.value(13).toDouble());
          spectrum_descr.precursor_ion_data =
            PrecursorIonData(q.value(10).toDouble(), q.value(11).toInt(), q.value(13).toDouble());
 
          spectrum_descr.precursor_id = q.value(6).toLongLong();
          spectrum_descr.ms2_index    = (spectrum_descr.precursor_id * 2) - 1;
          spectrum_descr.ms1_index    = (spectrum_descr.precursor_id * 2) - 2;
 
          spectrum_descr.scan_mobility_start = q.value(1).toLongLong();
          spectrum_descr.scan_mobility_end   = q.value(2).toLongLong();
 
          spectrum_descr.isolationMz     = q.value(3).toDouble();
          spectrum_descr.isolationWidth  = q.value(4).toDouble();
          spectrum_descr.collisionEnergy = q.value(5).toFloat();
          spectrum_descr.parent_frame    = q.value(14).toLongLong();
 
 
          first = false;
        }
    }
  if(spectrum_descr.precursor_id == 0)
    {
      throw ExceptionNotFound(QObject::tr("ERROR in %1 %2 : precursor id (%3) NOT FOUND ")
                                .arg(__FILE__)
                                .arg(__FUNCTION__)
                                .arg(precursor_id));
    }
  return spectrum_descr;
}
 
TimsDdaPrecursors::SpectrumDescr
TimsDdaPrecursors::getSpectrumDescrWithScanCoordinates(
  const std::pair<std::size_t, std::size_t> &scan_coordinates)
{
 
  SpectrumDescr spectrum_descr;
  QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
  QSqlQuery q(qdb);
  q.exec(QString("SELECT PasefFrameMsMsInfo.Frame, "    // 0
                 "PasefFrameMsMsInfo.ScanNumBegin, "    // 1
                 "PasefFrameMsMsInfo.ScanNumEnd, "      // 2
                 "PasefFrameMsMsInfo.IsolationMz, "     // 3
                 "PasefFrameMsMsInfo.IsolationWidth, "  // 4
                 "PasefFrameMsMsInfo.CollisionEnergy, " // 5
                 "PasefFrameMsMsInfo.Precursor, "       // 6
                 "Precursors.Id, "                      // 7
                 "Precursors.LargestPeakMz, "           // 8
                 "Precursors.AverageMz, "               // 9
                 "Precursors.MonoisotopicMz, "          // 10
                 "Precursors.Charge, "                  // 11
                 "Precursors.ScanNumber, "              // 12
                 "Precursors.Intensity, "               // 13
                 "Precursors.Parent "                   // 14
                 "FROM PasefFrameMsMsInfo "
                 "INNER JOIN Precursors ON "
                 "PasefFrameMsMsInfo.Precursor=Precursors.Id "
                 "WHERE "
                 "PasefFrameMsMsInfo.Frame=%1 and "
                 "(PasefFrameMsMsInfo.ScanNumBegin "
                 "<= %2 and PasefFrameMsMsInfo.ScanNumEnd >= %2);")
           .arg(scan_coordinates.first)
           .arg(scan_coordinates.second));
  if(q.lastError().isValid())
    {
 
      throw PappsoException(
        QObject::tr("ERROR in TIMS sqlite database file %1, executing SQL "
                    "command %2:\n%3\n%4\n%5")
          .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf"))
          .arg(q.lastQuery())
          .arg(qdb.lastError().databaseText())
          .arg(qdb.lastError().driverText())
          .arg(qdb.lastError().nativeErrorCode()));
    }
 
  if(q.next())
    {
 
      qDebug() << " cumul tims frame:" << q.value(0).toLongLong();
      spectrum_descr.tims_frame_list.push_back(q.value(0).toLongLong());
      // mass_spectrum.setPrecursorCharge(q.value(11).toInt());
      // mass_spectrum.setPrecursorMz(q.value(10).toDouble());
      // mass_spectrum.setPrecursorIntensity(q.value(13).toDouble());
      spectrum_descr.precursor_ion_data =
        PrecursorIonData(q.value(10).toDouble(), q.value(11).toInt(), q.value(13).toDouble());
 
      spectrum_descr.precursor_id = q.value(6).toLongLong();
      spectrum_descr.ms2_index    = (spectrum_descr.precursor_id * 2) - 1;
      spectrum_descr.ms1_index    = (spectrum_descr.precursor_id * 2) - 2;
 
      spectrum_descr.scan_mobility_start = q.value(1).toLongLong();
      spectrum_descr.scan_mobility_end   = q.value(2).toLongLong();
 
      spectrum_descr.isolationMz     = q.value(3).toDouble();
      spectrum_descr.isolationWidth  = q.value(4).toDouble();
      spectrum_descr.collisionEnergy = q.value(5).toFloat();
      spectrum_descr.parent_frame    = q.value(14).toLongLong();
    }
  return spectrum_descr;
}
 
 
void
TimsDdaPrecursors::fillSpectrumDescriptionWithSqlRecord(
  TimsDdaPrecursors::SpectrumDescr &spectrum_descr, QSqlQuery &qprecursor_list)
{
 
  spectrum_descr.tims_frame_list.clear();
  spectrum_descr.tims_frame_list.push_back(qprecursor_list.value(0).toLongLong());
  // mass_spectrum.setPrecursorCharge(q.value(11).toInt());
  // mass_spectrum.setPrecursorMz(q.value(10).toDouble());
  // mass_spectrum.setPrecursorIntensity(q.value(13).toDouble());
  spectrum_descr.precursor_ion_data = PrecursorIonData(qprecursor_list.value(10).toDouble(),
                                                       qprecursor_list.value(11).toInt(),
                                                       qprecursor_list.value(13).toDouble());
 
  spectrum_descr.precursor_id = qprecursor_list.value(6).toLongLong();
  spectrum_descr.ms2_index    = (spectrum_descr.precursor_id * 2) - 1;
  spectrum_descr.ms1_index    = (spectrum_descr.precursor_id * 2) - 2;
 
  spectrum_descr.scan_mobility_start = qprecursor_list.value(1).toLongLong();
  spectrum_descr.scan_mobility_end   = qprecursor_list.value(2).toLongLong();
 
  spectrum_descr.isolationMz     = qprecursor_list.value(3).toDouble();
  spectrum_descr.isolationWidth  = qprecursor_list.value(4).toDouble();
  spectrum_descr.collisionEnergy = qprecursor_list.value(5).toFloat();
  spectrum_descr.parent_frame    = qprecursor_list.value(14).toLongLong();
}
 
void
TimsDdaPrecursors::getQualifiedMs1MassSpectrumBySpectrumDescr(const MsRunIdCstSPtr &msrun_id,
                                                              QualifiedMassSpectrum &mass_spectrum,
                                                              const SpectrumDescr &spectrum_descr,
                                                              bool want_binary_data)
{
 
 
  qDebug() << " ms2_index=" << spectrum_descr.ms2_index
           << " precursor_index=" << spectrum_descr.precursor_id;
 
  TracePlusCombiner combiner;
  MapTrace combiner_result;
 
  try
    {
      mass_spectrum.setMsLevel(1);
      mass_spectrum.setPrecursorSpectrumIndex(0);
      mass_spectrum.setEmptyMassSpectrum(true);
 
      MassSpectrumId spectrum_id;
      spectrum_id.setSpectrumIndex(spectrum_descr.ms1_index);
      spectrum_id.setNativeId(QString("frame_id=%1 begin=%2 end=%3 precursor=%4 idxms1=%5")
                                .arg(spectrum_descr.parent_frame)
                                .arg(spectrum_descr.scan_mobility_start)
                                .arg(spectrum_descr.scan_mobility_end)
                                .arg(spectrum_descr.precursor_id)
                                .arg(spectrum_descr.ms1_index));
 
      spectrum_id.setMsRunId(msrun_id);
 
      mass_spectrum.setMassSpectrumId(spectrum_id);
 
 
      TimsFrameBaseCstSPtr tims_frame;
      if(want_binary_data)
        {
          qDebug() << "bindec";
          tims_frame = mp_timsDataOrigin->getTimsFrameCstSPtrCached(spectrum_descr.parent_frame);
        }
      else
        {
          tims_frame =
            mp_timsDataOrigin->getTimsFrameBaseCstSPtrCached(spectrum_descr.parent_frame);
        }
      mass_spectrum.setRtInSeconds(tims_frame.get()->getRtInSeconds());
 
      mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::IonMobOneOverK0Begin,
        tims_frame.get()->getOneOverK0Transformation(spectrum_descr.scan_mobility_start));
 
      mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::IonMobOneOverK0End,
        tims_frame.get()->getOneOverK0Transformation(spectrum_descr.scan_mobility_end));
 
 
      if(want_binary_data)
        {
          combiner.combine(combiner_result,
                           tims_frame.get()->cumulateScansToTrace(
                             spectrum_descr.scan_mobility_start, spectrum_descr.scan_mobility_end));
 
          Trace trace(combiner_result);
          qDebug();
 
          if(trace.size() > 0)
            {
              if(mcsp_ms1Filter != nullptr)
                {
                  mcsp_ms1Filter->filter(trace);
                }
 
              qDebug();
              mass_spectrum.setMassSpectrumSPtr(MassSpectrum(trace).makeMassSpectrumSPtr());
              mass_spectrum.setEmptyMassSpectrum(false);
            }
          else
            {
              mass_spectrum.setMassSpectrumSPtr(nullptr);
              mass_spectrum.setEmptyMassSpectrum(true);
            }
        }
      qDebug();
    }
 
  catch(PappsoException &error)
    {
      throw error;
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
}
 
 
void
TimsDdaPrecursors::setMs2FilterCstSPtr(FilterInterfaceCstSPtr &filter)
{
  mcsp_ms2Filter = filter;
}
void
TimsDdaPrecursors::setMs1FilterCstSPtr(FilterInterfaceCstSPtr &filter)
{
  mcsp_ms1Filter = filter;
}
 
void
TimsDdaPrecursors::getQualifiedMs2MassSpectrumBySpectrumDescr(const MsRunIdCstSPtr &msrun_id,
                                                              QualifiedMassSpectrum &mass_spectrum,
                                                              const SpectrumDescr &spectrum_descr,
                                                              bool want_binary_data)
{
 
 
  try
    {
      qDebug();
      MassSpectrumId spectrum_id;
 
      spectrum_id.setSpectrumIndex(spectrum_descr.ms2_index);
      spectrum_id.setNativeId(QString("precursor=%1 idxms2=%2")
                                .arg(spectrum_descr.precursor_id)
                                .arg(spectrum_descr.ms2_index));
      spectrum_id.setMsRunId(msrun_id);
 
      mass_spectrum.setMassSpectrumId(spectrum_id);
 
      mass_spectrum.setMsLevel(2);
      qDebug() << "spectrum_descr.precursor_id=" << spectrum_descr.precursor_id
               << " spectrum_descr.ms1_index=" << spectrum_descr.ms1_index
               << " spectrum_descr.ms2_index=" << spectrum_descr.ms2_index;
      mass_spectrum.setPrecursorSpectrumIndex(spectrum_descr.ms1_index);
 
      mass_spectrum.setEmptyMassSpectrum(true);
 
      qDebug();
 
 
      mass_spectrum.appendPrecursorIonData(spectrum_descr.precursor_ion_data);
 
      mass_spectrum.setPrecursorNativeId(
        QString("frame_id=%1 begin=%2 end=%3 precursor=%4 idxms1=%5")
          .arg(spectrum_descr.parent_frame)
          .arg(spectrum_descr.scan_mobility_start)
          .arg(spectrum_descr.scan_mobility_end)
          .arg(spectrum_descr.precursor_id)
          .arg(spectrum_descr.ms1_index));
 
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IsolationMz,
                                      spectrum_descr.isolationMz);
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IsolationMzWidth,
                                      spectrum_descr.isolationWidth);
 
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::CollisionEnergy,
                                      spectrum_descr.collisionEnergy);
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::BrukerPrecursorIndex,
                                      (quint64)spectrum_descr.precursor_id);
 
      // QMutexLocker locker(&m_mutex_spectrum);
      qDebug();
      TimsFrameBaseCstSPtr tims_frame, previous_frame;
      // TracePlusCombiner combiner;
      // MapTrace combiner_result;
      TimsDataFastMap &raw_spectrum = TimsDataFastMap::getTimsDataFastMapInstance();
      raw_spectrum.clear();
      bool first = true;
      for(std::size_t tims_id : spectrum_descr.tims_frame_list)
        {
          qDebug() << " precursor_index=" << spectrum_descr.precursor_id << " tims_id=" << tims_id
                   << (std::size_t)QThread::currentThreadId();
 
          if(want_binary_data)
            {
              qDebug() << "bindec";
              tims_frame = mp_timsDataOrigin->getTimsFrameCstSPtrCached(tims_id);
            }
          else
            {
              tims_frame = mp_timsDataOrigin->getTimsFrameBaseCstSPtrCached(tims_id);
            }
          qDebug() << (std::size_t)QThread::currentThreadId();
 
          if(first)
            {
              mass_spectrum.setRtInSeconds(tims_frame.get()->getRtInSeconds());
 
              mass_spectrum.setParameterValue(
                QualifiedMassSpectrumParameter::IonMobOneOverK0Begin,
                tims_frame.get()->getOneOverK0Transformation(spectrum_descr.scan_mobility_start));
 
              mass_spectrum.setParameterValue(
                QualifiedMassSpectrumParameter::IonMobOneOverK0End,
                tims_frame.get()->getOneOverK0Transformation(spectrum_descr.scan_mobility_end));
 
              first = false;
            }
 
 
          if(want_binary_data)
            {
              qDebug();
              /*combiner.combine(combiner_result,
                tims_frame.get()->cumulateScanToTrace(
                scan_mobility_start, scan_mobility_end));*/
              if(previous_frame.get() != nullptr)
                {
                  if(previous_frame.get()->hasSameCalibrationData(*tims_frame.get()))
                    {
                    }
                  else
                    {
                      throw ExceptionNotFound(
                        QObject::tr("ERROR in %1 %2, different calibration data "
                                    "between frame id %3 and frame id %4")
                          .arg(__FILE__)
                          .arg(__FUNCTION__)
                          .arg(previous_frame.get()->getId())
                          .arg(tims_frame.get()->getId()));
                    }
                }
              qDebug() << (std::size_t)QThread::currentThreadId();
 
              tims_frame.get()->combineScansInTofIndexIntensityMap(
                raw_spectrum, spectrum_descr.scan_mobility_start, spectrum_descr.scan_mobility_end);
              qDebug() << (std::size_t)QThread::currentThreadId();
            }
          previous_frame = tims_frame;
        }
      qDebug() << " precursor_index=" << spectrum_descr.precursor_id
               << " num_rows=" << spectrum_descr.tims_frame_list.size()
               << (std::size_t)QThread::currentThreadId();
      if(first == true)
        {
          throw ExceptionNotFound(
            QObject::tr("ERROR in getQualifiedMassSpectrumByPrecursorId, precursor "
                        "id=%1 not found")
              .arg(spectrum_descr.precursor_id));
        }
      if(want_binary_data)
        {
          qDebug() << " precursor_index=" << spectrum_descr.precursor_id;
          // peak_pick.filter(trace);
          Trace trace;
          if(m_builtinMs2Centroid)
            {
              // raw_spectrum.removeArtefactualSpike();
              raw_spectrum.builtInCentroid();
            }
 
          trace = tims_frame.get()->getTraceFromTofIndexIntensityMap(raw_spectrum);
 
          if(trace.size() > 0)
            {
              qDebug() << " precursor_index=" << spectrum_descr.precursor_id << " " << trace.size()
                       << " " << (std::size_t)QThread::currentThreadId();
 
              if(mcsp_ms2Filter != nullptr)
                {
                  // FilterTriangle filter;
                  // filter.setTriangleSlope(50, 0.02);
                  // filter.filter(trace);
                  // trace.filter(FilterHighPass(10));
                  mcsp_ms2Filter->filter(trace);
                }
 
              // FilterScaleFactorY filter_scale((double)1 /
              // (double)tims_frame_list.size());
              // filter_scale.filter(trace);
              qDebug() << " precursor_index=" << spectrum_descr.precursor_id;
              mass_spectrum.setMassSpectrumSPtr(MassSpectrum(trace).makeMassSpectrumSPtr());
              mass_spectrum.setEmptyMassSpectrum(false);
            }
          else
            {
              mass_spectrum.setMassSpectrumSPtr(nullptr);
              mass_spectrum.setEmptyMassSpectrum(true);
            }
 
          qDebug();
        }
      qDebug();
    }
 
  catch(PappsoException &error)
    {
      throw PappsoException(QObject::tr("ERROR in %1 (ms2_index=%2 precursor_index=%3):\n%4")
                              .arg(__FUNCTION__)
                              .arg(spectrum_descr.ms2_index)
                              .arg(spectrum_descr.precursor_id)
                              .arg(error.qwhat()));
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
  qDebug();
}
 
void
TimsDdaPrecursors::setMs2BuiltinCentroid(bool centroid)
{
  m_builtinMs2Centroid = centroid;
}
 
bool
TimsDdaPrecursors::getMs2BuiltinCentroid() const
{
  return m_builtinMs2Centroid;
}
 
void
TimsDdaPrecursors::ms2ReaderSpectrumCollectionByMsLevel(const MsRunIdCstSPtr &msrun_id,
                                                        SpectrumCollectionHandlerInterface &handler,
                                                        unsigned int ms_level)
{
  qDebug() << " ms_level=" << ms_level;
 
  QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
  QSqlQuery qprecursor_list(qdb);
  qprecursor_list.exec(
    QString("SELECT PasefFrameMsMsInfo.Frame, "    // 0
            "PasefFrameMsMsInfo.ScanNumBegin, "    // 1
            "PasefFrameMsMsInfo.ScanNumEnd, "      // 2
            "PasefFrameMsMsInfo.IsolationMz, "     // 3
            "PasefFrameMsMsInfo.IsolationWidth, "  // 4
            "PasefFrameMsMsInfo.CollisionEnergy, " // 5
            "PasefFrameMsMsInfo.Precursor, "       // 6
            "Precursors.Id, "                      // 7
            "Precursors.LargestPeakMz, "           // 8
            "Precursors.AverageMz, "               // 9
            "Precursors.MonoisotopicMz, "          // 10
            "Precursors.Charge, "                  // 11
            "Precursors.ScanNumber, "              // 12
            "Precursors.Intensity, "               // 13
            "Precursors.Parent "                   // 14
            "FROM PasefFrameMsMsInfo "
            "INNER JOIN Precursors ON "
            "PasefFrameMsMsInfo.Precursor=Precursors.Id "
            "ORDER BY PasefFrameMsMsInfo.Precursor, PasefFrameMsMsInfo.Frame ;"));
  if(qprecursor_list.lastError().isValid())
    {
 
      throw PappsoException(
        QObject::tr("ERROR in TIMS sqlite database file %1, executing SQL "
                    "command %2:\n%3\n%4\n%5")
          .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf"))
          .arg(qprecursor_list.lastQuery())
          .arg(qdb.lastError().databaseText())
          .arg(qdb.lastError().driverText())
          .arg(qdb.lastError().nativeErrorCode()));
    }
 
 
  qDebug() << "qprecursor_list.size()=" << qprecursor_list.size();
  qDebug() << QObject::tr(
                "TIMS sqlite database file %1, executing SQL "
                "command %2:\n%3\n%4\n%5")
                .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf"))
                .arg(qprecursor_list.lastQuery())
                .arg(qdb.lastError().databaseText())
                .arg(qdb.lastError().driverText())
                .arg(qdb.lastError().nativeErrorCode());
 
  qDebug() << "qprecursor_list.isActive()=" << qprecursor_list.isActive();
  qDebug() << "qprecursor_list.isSelect()=" << qprecursor_list.isSelect();
  bool first = true;
  SpectrumDescr spectrum_descr;
 
  /*
std::size_t i = 0;
while(qprecursor_list.next())
                         {
                         qDebug() << "i=" << i;
                         i++;
                         }*/
 
  qprecursor_list.last(); // strange bug : get the last sql record and get
  // back, otherwise it will not retrieve all records.
 
  qDebug() << "qprecursor_list.at()=" << qprecursor_list.at();
  qprecursor_list.first();
  std::vector<TimsDdaPrecursors::SpectrumDescr> spectrum_description_list;
  spectrum_descr.precursor_id = 0;
  // std::size_t i = 0;
 
  do
    {
 
      if(spectrum_descr.precursor_id != (std::size_t)qprecursor_list.value(6).toLongLong())
        {
          // new precursor
          if(spectrum_descr.precursor_id > 0)
            {
              spectrum_description_list.push_back(spectrum_descr);
            }
 
          spectrum_descr.tims_frame_list.clear();
          first = true;
        }
      qDebug() << " qprecursor_list.value(6).toLongLong() ="
               << qprecursor_list.value(6).toLongLong();
      spectrum_descr.precursor_id = (std::size_t)qprecursor_list.value(6).toLongLong();
      qDebug() << " spectrum_descr.precursor_id =" << spectrum_descr.precursor_id;
      qDebug() << " cumul tims frame:" << qprecursor_list.value(0).toLongLong();
      spectrum_descr.tims_frame_list.push_back(qprecursor_list.value(0).toLongLong());
      qDebug() << " first =" << first;
      if(first)
        {
          qDebug();
          // mass_spectrum.setPrecursorCharge(q.value(11).toInt());
          // mass_spectrum.setPrecursorMz(q.value(10).toDouble());
          // mass_spectrum.setPrecursorIntensity(q.value(13).toDouble());
          spectrum_descr.precursor_ion_data =
            PrecursorIonData(qprecursor_list.value(10).toDouble(),
                             qprecursor_list.value(11).toInt(),
                             qprecursor_list.value(13).toDouble());
 
          // spectrum_descr.precursor_id = q.value(6).toLongLong();
          spectrum_descr.ms2_index = (spectrum_descr.precursor_id * 2) - 1;
          spectrum_descr.ms1_index = (spectrum_descr.precursor_id * 2) - 2;
 
          spectrum_descr.scan_mobility_start = qprecursor_list.value(1).toLongLong();
          spectrum_descr.scan_mobility_end   = qprecursor_list.value(2).toLongLong();
 
          spectrum_descr.isolationMz     = qprecursor_list.value(3).toDouble();
          spectrum_descr.isolationWidth  = qprecursor_list.value(4).toDouble();
          spectrum_descr.collisionEnergy = qprecursor_list.value(5).toFloat();
          spectrum_descr.parent_frame    = qprecursor_list.value(14).toLongLong();
 
 
          first = false;
        }
      // qDebug() << "qprecursor_list.executedQuery()="
      //         << qprecursor_list.executedQuery();
      // qDebug() << "qprecursor_list.last()=" << qprecursor_list.last();
      // i++;
    }
  while(qprecursor_list.next());
 
  // last One
 
  // new precursor
  if(spectrum_descr.precursor_id > 0)
    {
      spectrum_description_list.push_back(spectrum_descr);
    }
 
 
  QString local_filepath =
    mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf");
 
  if(m_isMonoThread)
    {
      for(SpectrumDescr &spectrum_descr : spectrum_description_list)
        {
 
          std::vector<QualifiedMassSpectrum> mass_spectrum_list;
          ms2ReaderGenerateMS1MS2Spectrum(
            msrun_id, mass_spectrum_list, handler, spectrum_descr, ms_level);
 
          for(auto &qualified_spectrum : mass_spectrum_list)
            {
              handler.setQualifiedMassSpectrum(qualified_spectrum);
            }
 
          if(handler.shouldStop())
            {
              qDebug() << "The operation was cancelled. Breaking the loop.";
              throw ExceptionInterrupted(
                QObject::tr("reading TimsTOF job cancelled by the user :\n%1").arg(local_filepath));
            }
        }
    }
  else
    {
 
 
      TimsDdaPrecursors *itself                           = this;
      SpectrumCollectionHandlerInterface *pointer_handler = &handler;
 
 
      std::function<std::vector<QualifiedMassSpectrum>(const TimsDdaPrecursors::SpectrumDescr &)>
        map_function_generate_spectrum = [itself, msrun_id, pointer_handler, ms_level](
                                           const TimsDdaPrecursors::SpectrumDescr &spectrum_descr)
        -> std::vector<QualifiedMassSpectrum> {
        std::vector<QualifiedMassSpectrum> mass_spectrum_list;
        itself->ms2ReaderGenerateMS1MS2Spectrum(
          msrun_id, mass_spectrum_list, *pointer_handler, spectrum_descr, ms_level);
 
 
        return mass_spectrum_list;
      };
 
      std::function<void(std::size_t,
                         const std::vector<QualifiedMassSpectrum> &qualified_spectrum_list)>
        reduce_function_spectrum_list =
          [pointer_handler, local_filepath](
            std::size_t res, const std::vector<QualifiedMassSpectrum> &qualified_spectrum_list) {
            for(auto &qualified_spectrum : qualified_spectrum_list)
              {
                pointer_handler->setQualifiedMassSpectrum(qualified_spectrum);
              }
 
            if(pointer_handler->shouldStop())
              {
                qDebug() << "The operation was cancelled. Breaking the loop.";
                throw ExceptionInterrupted(
                  QObject::tr("reading TimsTOF job on %1 cancelled by the user")
                    .arg(local_filepath));
              }
            res++;
          };
 
 
      QFuture<std::size_t> res;
      res = QtConcurrent::mappedReduced<std::size_t>(spectrum_description_list.begin(),
                                                     spectrum_description_list.end(),
                                                     map_function_generate_spectrum,
                                                     reduce_function_spectrum_list,
                                                     QtConcurrent::OrderedReduce);
      res.waitForFinished();
    }
  handler.loadingEnded();
  mp_timsDataOrigin->getTimsBinDecPtr()->closeLinearRead();
}
 
void
TimsDdaPrecursors::setMonoThread(bool is_mono_thread)
{
  m_isMonoThread = is_mono_thread;
}
 
void
TimsDdaPrecursors::ms2ReaderGenerateMS1MS2Spectrum(
  const MsRunIdCstSPtr &msrun_id,
  std::vector<QualifiedMassSpectrum> &qualified_mass_spectrum_list,
  SpectrumCollectionHandlerInterface &handler,
  const TimsDdaPrecursors::SpectrumDescr &spectrum_descr,
  unsigned int ms_level)
{
 
  qDebug() << " ms_level=" << ms_level;
  // The handler will receive the index of the mass spectrum in the
  // current run via the mass spectrum id member datum.
  if((ms_level == 0) || (ms_level == 1))
    {
      qualified_mass_spectrum_list.push_back(QualifiedMassSpectrum());
      getQualifiedMs1MassSpectrumBySpectrumDescr(msrun_id,
                                                 qualified_mass_spectrum_list.back(),
                                                 spectrum_descr,
                                                 handler.needMsLevelPeakList(1));
    }
  if((ms_level == 0) || (ms_level == 2))
    {
      qualified_mass_spectrum_list.push_back(QualifiedMassSpectrum());
      getQualifiedMs2MassSpectrumBySpectrumDescr(msrun_id,
                                                 qualified_mass_spectrum_list.back(),
                                                 spectrum_descr,
                                                 handler.needMsLevelPeakList(2));
    }
  qDebug();
}
 
void
TimsDdaPrecursors::rawReaderSpectrumCollectionByMsLevel(const MsRunIdCstSPtr &msrun_id,
                                                        SpectrumCollectionHandlerInterface &handler,
                                                        unsigned int ms_level)
{
 
 
  // We'll need it to perform the looping in the spectrum list.
  std::size_t spectrum_list_size = mp_timsDataOrigin->getTotalScanCount();
 
  //   qDebug() << "The spectrum list has size:" << spectrum_list_size;
 
  // Inform the handler of the spectrum list so that it can handle feedback to
  // the user.
  handler.spectrumListHasSize(spectrum_list_size);
 
  QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
  QSqlQuery qprecursor_list(qdb);
  qprecursor_list.exec(
    QString("SELECT DISTINCT "
            "PasefFrameMsMsInfo.Frame, "           // 0
            "PasefFrameMsMsInfo.ScanNumBegin, "    // 1
            "PasefFrameMsMsInfo.ScanNumEnd, "      // 2
            "PasefFrameMsMsInfo.IsolationMz, "     // 3
            "PasefFrameMsMsInfo.IsolationWidth, "  // 4
            "PasefFrameMsMsInfo.CollisionEnergy, " // 5
            "PasefFrameMsMsInfo.Precursor, "       // 6
            "Precursors.Id, "                      // 7
            "Precursors.LargestPeakMz, "           // 8
            "Precursors.AverageMz, "               // 9
            "Precursors.MonoisotopicMz, "          // 10
            "Precursors.Charge, "                  // 11
            "Precursors.ScanNumber, "              // 12
            "Precursors.Intensity, "               // 13
            "Precursors.Parent "                   // 14
            "FROM PasefFrameMsMsInfo "
            "INNER JOIN Precursors ON "
            "PasefFrameMsMsInfo.Precursor=Precursors.Id "
            "ORDER BY PasefFrameMsMsInfo.Frame, PasefFrameMsMsInfo.ScanNumBegin ;"));
  if(qprecursor_list.lastError().isValid())
    {
      throw PappsoException(
        QObject::tr("ERROR in TIMS sqlite database file %1, executing SQL "
                    "command %2:\n%3\n%4\n%5")
          .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf"))
          .arg(qprecursor_list.lastQuery())
          .arg(qdb.lastError().databaseText())
          .arg(qdb.lastError().driverText())
          .arg(qdb.lastError().nativeErrorCode()));
    }
 
 
  std::size_t i = 0; // iterate on each Spectrum
 
  qprecursor_list.last(); // strange bug : get the last sql record and get
  // back, unless it will not retrieve all records.
 
  qDebug() << "qprecursor_list.at()=" << qprecursor_list.at();
  qprecursor_list.first();
 
  TimsFrameBaseCstSPtr tims_frame;
  SpectrumDescr spectrum_descr;
 
  for(const FrameIdDescr &current_frame : mp_timsDataOrigin->getFrameIdDescrList())
    {
 
      // If the user of this reader instance wants to stop reading the
      // spectra, then break this loop.
      if(handler.shouldStop())
        {
          qDebug() << "The operation was cancelled. Breaking the loop.";
          throw ExceptionInterrupted(
            QObject::tr("reading TimsTOF job cancelled by the user :\n%1")
              .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf")));
        }
 
      tims_frame = mp_timsDataOrigin->getTimsFrameBaseCstSPtrCached(current_frame.m_frameId);
      unsigned int tims_ms_level = tims_frame.get()->getMsLevel();
 
      if((ms_level != 0) && (ms_level != tims_ms_level))
        { // bypass
          i += current_frame.m_scanCount;
        }
      else
        {
          bool want_binary_data = handler.needMsLevelPeakList(tims_ms_level);
          qDebug() << "want_binary_data=" << want_binary_data;
          if(want_binary_data)
            {
              qDebug() << "bindec";
              tims_frame = mp_timsDataOrigin->getTimsFrameCstSPtrCached(current_frame.m_frameId);
            }
 
          bool possible_precursor = false;
          if(tims_ms_level == 2)
            {
              // seek the precursor record:
              while(qprecursor_list.value(0).toULongLong() < current_frame.m_frameId)
                {
                  qprecursor_list.next();
 
                  if(qprecursor_list.value(0).toULongLong() == current_frame.m_frameId)
                    {
                      possible_precursor = true;
                    }
                  fillSpectrumDescriptionWithSqlRecord(spectrum_descr, qprecursor_list);
                }
            }
 
          for(std::size_t scan_num = 0; scan_num < current_frame.m_scanCount; scan_num++)
            {
              bool has_a_precursor = false;
              if(possible_precursor)
                {
                  if(spectrum_descr.scan_mobility_end < scan_num)
                    {
                      // seek the precursor record:
                      while(qprecursor_list.value(0).toULongLong() < current_frame.m_frameId)
                        {
                          qprecursor_list.next();
 
                          if(qprecursor_list.value(0).toULongLong() != current_frame.m_frameId)
                            {
                              possible_precursor = false;
                            }
                          fillSpectrumDescriptionWithSqlRecord(spectrum_descr, qprecursor_list);
                        }
                    }
 
                  if(possible_precursor && (spectrum_descr.scan_mobility_start < scan_num))
                    {
                      // we are in
                      has_a_precursor = true;
                    }
                } // end to determine if we are in a precursor for this
                  // spectrum
 
              QualifiedMassSpectrum mass_spectrum;
 
 
              MassSpectrumId spectrum_id;
 
              spectrum_id.setSpectrumIndex(i);
              spectrum_id.setMsRunId(msrun_id);
              spectrum_id.setNativeId(QString("frame_id=%1 scan_index=%2 global_scan_index=%3")
                                        .arg(current_frame.m_frameId)
                                        .arg(scan_num)
                                        .arg(i));
 
              mass_spectrum.setMassSpectrumId(spectrum_id);
 
              mass_spectrum.setMsLevel(tims_frame.get()->getMsLevel());
              mass_spectrum.setRtInSeconds(tims_frame.get()->getRtInSeconds());
 
              mass_spectrum.setDtInMilliSeconds(
                tims_frame.get()->getDriftTimeInMilliseconds(scan_num));
              // 1/K0
              mass_spectrum.setParameterValue(
                QualifiedMassSpectrumParameter::IonMobOneOverK0,
                tims_frame.get()->getOneOverK0Transformation(scan_num));
 
              mass_spectrum.setEmptyMassSpectrum(true);
              if(want_binary_data)
                {
                  try
                    {
                      mass_spectrum.setMassSpectrumSPtr(
                        tims_frame.get()->getMassSpectrumSPtr(scan_num));
                    }
                  catch(PappsoException &error)
                    {
                      throw PappsoException(
                        QObject::tr("ERROR in %1 (scan_num=%2 spectrum_index=%3):\n%4")
                          .arg(__FUNCTION__)
                          .arg(scan_num)
                          .arg(spectrum_id.getSpectrumIndex())
                          .arg(error.qwhat()));
                    }
                  if(mass_spectrum.size() > 0)
                    {
                      mass_spectrum.setEmptyMassSpectrum(false);
                    }
                }
              else
                {
                  // if(tims_frame.get()->getNbrPeaks(coordinate.second) > 0)
                  //{
                  mass_spectrum.setEmptyMassSpectrum(false);
                  // }
                }
              if(has_a_precursor)
                {
                  if(spectrum_descr.precursor_id > 0)
                    {
 
                      mass_spectrum.appendPrecursorIonData(spectrum_descr.precursor_ion_data);
 
                      std::size_t prec_spectrum_index =
                        mp_timsDataOrigin->getGlobalScanIndexByScanCoordinates(
                          spectrum_descr.parent_frame, scan_num);
 
                      mass_spectrum.setPrecursorSpectrumIndex(prec_spectrum_index);
                      mass_spectrum.setPrecursorNativeId(
                        QString("frame_id=%1 scan_index=%2 global_scan_index=%3")
                          .arg(spectrum_descr.parent_frame)
                          .arg(scan_num)
                          .arg(prec_spectrum_index));
 
                      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IsolationMz,
                                                      spectrum_descr.isolationMz);
                      mass_spectrum.setParameterValue(
                        QualifiedMassSpectrumParameter::IsolationMzWidth,
                        spectrum_descr.isolationWidth);
 
                      mass_spectrum.setParameterValue(
                        QualifiedMassSpectrumParameter::CollisionEnergy,
                        spectrum_descr.collisionEnergy);
                      mass_spectrum.setParameterValue(
                        QualifiedMassSpectrumParameter::BrukerPrecursorIndex,
                        (quint64)spectrum_descr.precursor_id);
                    }
                }
 
              handler.setQualifiedMassSpectrum(mass_spectrum);
              i++;
            }
        }
    }
}
 
 
std::vector<std::size_t>
TimsDdaPrecursors::getPrecursorsByMzRtCharge(int charge, double mz_val, double rt_sec, double k0)
{
 
  std::vector<std::size_t> precursor_ids;
  std::vector<std::vector<double>> ids;
 
  QSqlDatabase qdb = mp_timsDataOrigin->openDatabaseConnection();
  QSqlQuery q(QString("SELECT Frames.Time, Precursors.MonoisotopicMz, Precursors.Charge, "
                      "Precursors.Id, Frames.Id, PasefFrameMsMsInfo.ScanNumBegin, "
                      "PasefFrameMsMsInfo.scanNumEnd "
                      "FROM Frames "
                      "INNER JOIN PasefFrameMsMsInfo ON Frames.Id = "
                      "PasefFrameMsMsInfo.Frame "
                      "INNER JOIN Precursors ON "
                      "PasefFrameMsMsInfo.Precursor= Precursors.Id "
                      "WHERE Precursors.Charge == %1 "
                      "AND Precursors.MonoisotopicMz > %2 -0.01 "
                      "AND Precursors.MonoisotopicMz < %2 +0.01 "
                      "AND Frames.Time >= %3 -1 "
                      "AND Frames.Time < %3 +1; ")
                .arg(charge)
                .arg(mz_val)
                .arg(rt_sec));
  q.exec();
  if(q.lastError().isValid())
    {
 
      throw PappsoException(
        QObject::tr("ERROR in TIMS sqlite database file %1, database name %2, "
                    "executing SQL "
                    "command %3:\n%4\n%5\n%6")
          .arg(mp_timsDataOrigin->getTimsDataDirectory().absoluteFilePath("analysis.tdf"))
          .arg(qdb.databaseName())
          .arg(q.lastQuery())
          .arg(qdb.lastError().databaseText())
          .arg(qdb.lastError().driverText())
          .arg(qdb.lastError().nativeErrorCode()));
    }
  while(q.next())
    {
      // qInfo() << q.value(0).toDouble() << q.value(1).toDouble()
      //       << q.value(2).toDouble() << q.value(3).toDouble();
 
      std::vector<double> sql_values;
      sql_values.push_back(q.value(4).toDouble()); // frame id
      sql_values.push_back(q.value(3).toDouble()); // precursor id
      sql_values.push_back(q.value(5).toDouble()); // scan num begin
      sql_values.push_back(q.value(6).toDouble()); // scan num end
      sql_values.push_back(q.value(1).toDouble()); // mz_value
 
      ids.push_back(sql_values);
 
 
      if(std::find(precursor_ids.begin(), precursor_ids.end(), q.value(3).toDouble()) ==
         precursor_ids.end())
        {
          precursor_ids.push_back(q.value(3).toDouble());
        }
    }
 
  if(precursor_ids.size() > 1)
    {
      // std::vector<std::size_t> precursor_ids_ko =
      // getMatchPrecursorIdByKo(ids, values[3]);
      if(precursor_ids.size() > 1)
        {
          precursor_ids = getClosestPrecursorIdByMz(ids, k0);
        }
      return precursor_ids;
    }
  else
    {
      return precursor_ids;
    }
}
 
std::vector<std::size_t>
TimsDdaPrecursors::getMatchPrecursorIdByKo(std::vector<std::vector<double>> ids, double ko_value)
{
  std::vector<std::size_t> precursor_id;
  for(std::vector<double> index : ids)
    {
      auto coordinate = mp_timsDataOrigin->getScanCoordinatesByGlobalScanIndex(index[0]);
 
      TimsFrameBaseCstSPtr tims_frame;
      tims_frame = mp_timsDataOrigin->getTimsFrameBaseCstSPtrCached(coordinate.first);
 
      double bko = tims_frame.get()->getOneOverK0Transformation(index[2]);
      double eko = tims_frame.get()->getOneOverK0Transformation(index[3]);
 
      // qInfo() << "diff" << (bko + eko) / 2;
      double mean_ko = (bko + eko) / 2;
 
      if(mean_ko > ko_value - 0.1 && mean_ko < ko_value + 0.1)
        {
          precursor_id.push_back(index[1]);
        }
    }
  return precursor_id;
}
 
std::vector<std::size_t>
TimsDdaPrecursors::getClosestPrecursorIdByMz(std::vector<std::vector<double>> ids, double mz_value)
{
  std::vector<std::size_t> best_precursor;
  double best_value     = 1;
  int count             = 1;
  int best_val_position = 0;
 
  for(std::vector<double> values : ids)
    {
      double new_val = abs(mz_value - values[4]);
      if(new_val < best_value)
        {
          best_value        = new_val;
          best_val_position = count;
        }
      count++;
    }
  best_precursor.push_back(ids[best_val_position][1]);
  return best_precursor;
}
 
 
} // namespace pappso