timsframebase.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/timsframebase.cpp
 * \date 16/12/2019
 * \author Olivier Langella
 * \brief handle a single Bruker's TimsTof frame without binary data
 */
 
/*******************************************************************************
 * 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 "timsframebase.h"
#include "../../../pappsomspp/pappsoexception.h"
#include "../../../pappsomspp/exception/exceptionoutofrange.h"
#include <QDebug>
#include <cmath>
#include <solvers.h>
 
namespace pappso
{
 
TimsFrameBase::TimsFrameBase(std::size_t timsId, quint32 scanNum)
{
  qDebug() << timsId;
  m_timsId = timsId;
 
  m_scanNumber = scanNum;
}
 
TimsFrameBase::TimsFrameBase(const TimsFrameBase &other)
{
}
 
TimsFrameBase::~TimsFrameBase()
{
}
 
void
TimsFrameBase::setAccumulationTime(double accumulation_time_ms)
{
  m_accumulationTime = accumulation_time_ms;
}
 
 
void
TimsFrameBase::setMzCalibration(double temperature_correction,
                                double digitizerTimebase,
                                double digitizerDelay,
                                double C0,
                                double C1,
                                double C2,
                                double C3)
{
 
  // temperature compensation
  C1 = C1 * temperature_correction;
  C2 = C2 / temperature_correction;
 
  m_digitizerDelay    = digitizerDelay;
  m_digitizerTimebase = digitizerTimebase;
 
  m_mzCalibrationArr.push_back(C0);
  m_mzCalibrationArr.push_back(std::sqrt(std::pow(10, 12) / C1));
  m_mzCalibrationArr.push_back(C2);
  m_mzCalibrationArr.push_back(C3);
}
 
bool
TimsFrameBase::checkScanNum(std::size_t scanNum) const
{
  if(scanNum < 0)
    {
      throw pappso::ExceptionOutOfRange(
        QObject::tr("Invalid scan number : scanNum%1 < 0").arg(scanNum));
    }
  if(scanNum >= m_scanNumber)
    {
      throw pappso::ExceptionOutOfRange(
        QObject::tr("Invalid scan number : scanNum%1  > m_scanNumber")
          .arg(scanNum));
    }
 
  return true;
}
 
std::size_t
TimsFrameBase::getNbrPeaks(std::size_t scanNum) const
{
  throw PappsoException(
    QObject::tr(
      "ERROR unable to get number of peaks in TimsFrameBase for scan number %1")
      .arg(scanNum));
}
 
MassSpectrumSPtr
TimsFrameBase::getMassSpectrumSPtr(std::size_t scanNum) const
{
  throw PappsoException(
    QObject::tr(
      "ERROR unable to getMassSpectrumSPtr in TimsFrameBase for scan number %1")
      .arg(scanNum));
}
Trace
TimsFrameBase::cumulateScanToTrace(std::size_t scanNumBegin,
                                   std::size_t scanNumEnd) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to cumulateScanToTrace in TimsFrameBase for scan "
                "number begin %1 end %2")
      .arg(scanNumBegin)
      .arg(scanNumEnd));
}
double
TimsFrameBase::getTofFromIndex(double index) const
{
  // mz calibration
  return (index * m_digitizerTimebase) + m_digitizerDelay;
}
double
TimsFrameBase::getTofFromIndex(quint32 index) const
{
  // mz calibration
  return ((double)index * m_digitizerTimebase) + m_digitizerDelay;
}
double
TimsFrameBase::getMzFromTof(double tof) const
{
  // http://www.alglib.net/equations/polynomial.php
  // http://www.alglib.net/translator/man/manual.cpp.html#sub_polynomialsolve
  // https://math.stackexchange.com/questions/1291208/number-of-roots-of-a-polynomial-of-non-integer-degree
  // https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=2&ved=2ahUKEwiWhLOFxqrkAhVLxYUKHVqqDFcQFjABegQIAxAB&url=https%3A%2F%2Fkluge.in-chemnitz.de%2Fopensource%2Fspline%2Fexample_alglib.cpp&usg=AOvVaw0guGejJGPmkOVg48_GJYR8
  // https://stackoverflow.com/questions/26091323/how-to-plot-a-function-curve-in-r
  /*
   * beware to put the function on a single line in R:
> eq <- function(m){ 1 + (sqrt((10^12)/670) * sqrt(m)) + (207.775676931964 * m)
+ (59.2526676368822 * (m^1.5)) }
> eq <- function(m){ 313.577620892277 + (sqrt((10^12)/157424.07710945) *
sqrt(m)) + (0.000338743021989553 * m)
+ (0 * (m^1.5)) }
> plot(eq(1:1000), type='l')
 
 
 
> eq2 <- function(m2){ 1 + sqrt((10^12)/670) * m2 + 207.775676931964 * (m2^2)
+ 59.2526676368822 * (m2^3) }
> plot(eq2(1:sqrt(1000)), type='l')
*/
  // How to Factor a Trinomial with Fractions as Coefficients
 
  // formula
  // a = c0 = 1
  // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
  // c = c2, c2 = 207.775676931964  * m
  // d = c3, c3 = 59.2526676368822  * m^1.5  (3/2)
  // double mz = 0;
  std::vector<double> X;
  X.push_back((m_mzCalibrationArr[0] - (double)tof));
  X.push_back(m_mzCalibrationArr[1]);
  if(m_mzCalibrationArr[2] != 0)
    {
      X.push_back(m_mzCalibrationArr[2]);
    }
  if(m_mzCalibrationArr[3] != 0)
    {
      X.push_back(m_mzCalibrationArr[3]);
    }
  else
    {
      qDebug() << m_mzCalibrationArr[3];
    }
 
  alglib::real_1d_array polynom_array;
  try
    {
      polynom_array.setcontent(X.size(), &(X[0]));
    }
  catch(alglib::ap_error &error)
    {
      // PolynomialSolve: A[N]=0
      throw pappso::PappsoException(
        QObject::tr("ERROR in alglib::polynom_array.setcontent :\n%1")
          .arg(error.msg.c_str()));
    }
 
 
  /*
  alglib::polynomialsolve(
real_1d_array a,
ae_int_t n,
complex_1d_array& x,
polynomialsolverreport& rep,
const xparams _params = alglib::xdefault);
*/
  alglib::complex_1d_array m;
  alglib::polynomialsolverreport rep;
  // qDebug();
  try
    {
      alglib::polynomialsolve(polynom_array, X.size() - 1, m, rep);
    }
  catch(alglib::ap_error &error)
    {
      qDebug() << " X.size() - 1 = " << X.size() - 1;
      qDebug() << m_mzCalibrationArr[0];
      qDebug() << m_mzCalibrationArr[1];
      qDebug() << m_mzCalibrationArr[2];
      qDebug() << m_mzCalibrationArr[3];
 
      // PolynomialSolve: A[N]=0
      throw pappso::PappsoException(
        QObject::tr("ERROR in alglib::polynomialsolve :\n%1")
          .arg(error.msg.c_str()));
    }
 
 
  // qDebug();
 
  if(m.length() == 0)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR in TimsFrame::getMzFromTof m.size() == 0"));
    }
  // qDebug();
  if(m[0].y != 0)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR in TimsFrame::getMzFromTof m[0].y!= 0"));
    }
 
  return pow(m[0].x, 2);
}
 
quint32
TimsFrameBase::getRawIndexFromMz(double mz) const
{
  // formula
  // a = c0 = 1
  // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
  // c = c2, c2 = 207.775676931964  * m
  // d = c3, c3 = 59.2526676368822  * m^1.5  (3/2)
  qDebug() << "mz=" << mz;
 
  double tof = m_mzCalibrationArr[0];
  qDebug() << "tof ( m_mzCalibrationArr[0])=" << tof;
  // TODO cache value of  std::sqrt((std::pow(10, 12) / m_mzCalibrationArr[1]))
  tof += m_mzCalibrationArr[1] * std::sqrt(mz);
  qDebug() << "tof=" << tof;
  tof += m_mzCalibrationArr[2] * mz;
  qDebug() << "tof=" << tof;
  tof += m_mzCalibrationArr[3] * std::pow(mz, 1.5);
  qDebug() << "tof=" << tof;
  tof -= m_digitizerDelay;
  qDebug() << "tof=" << tof;
  tof = tof / m_digitizerTimebase;
  qDebug() << "index=" << tof;
  return (quint32)std::round(tof);
}
 
void
TimsFrameBase::setTime(double time)
{
  m_time = time;
}
 
void
TimsFrameBase::setMsMsType(quint8 type)
{
 
  qDebug() << " m_msMsType=" << type;
  m_msMsType = type;
}
 
unsigned int
TimsFrameBase::getMsLevel() const
{
  if(m_msMsType == 0)
    return 1;
  return 2;
}
 
double
TimsFrameBase::getTime() const
{
  return m_time;
}
 
std::size_t
TimsFrameBase::getId() const
{
  return m_timsId;
}
void
TimsFrameBase::setTimsCalibration(int tims_model_type,
                                  double C0,
                                  double C1,
                                  double C2,
                                  double C3,
                                  double C4,
                                  double C5,
                                  double C6,
                                  double C7,
                                  double C8,
                                  double C9)
{
  if(tims_model_type != 2)
    {
      throw pappso::PappsoException(QObject::tr(
        "ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
    }
  m_timsDvStart = C2; // C2 from TimsCalibration
  m_timsTtrans  = C4; // C4 from TimsCalibration
  m_timsNdelay  = C0; // C0 from TimsCalibration
  m_timsVmin    = C8; // C8 from TimsCalibration
  m_timsVmax    = C9; // C9 from TimsCalibration
  m_timsC6      = C6;
  m_timsC7      = C7;
 
 
  m_timsSlope =
    (C3 - m_timsDvStart) / C1; //  //C3 from TimsCalibration // C2 from
                               //  TimsCalibration // C1 from TimsCalibration
}
double
TimsFrameBase::getVoltageTransformation(std::size_t scanNum) const
{
  double v = m_timsDvStart +
             m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
 
  if(v < m_timsVmin)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
                    "calibration, v < m_timsVmin"));
    }
 
 
  if(v > m_timsVmax)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
                    "calibration, v > m_timsVmax"));
    }
  return v;
}
double
TimsFrameBase::getDriftTime(std::size_t scanNum) const
{
  return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
}
 
double
TimsFrameBase::getOneOverK0Transformation(std::size_t scanNum) const
{
  return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
}
} // namespace pappso