mzcalibrationmodel1.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/mzcalibration/mzcalibrationmodel1.cpp
 * \date 11/11/2020
 * \author Olivier Langella
 * \brief implement Bruker's model type 1 formula to compute m/z
 */
 
/*******************************************************************************
 * Copyright (c) 2020 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 "mzcalibrationmodel1.h"
#include <cmath>
#include <QDebug>
#include <QObject>
#include "pappsomspp/core/pappsoexception.h"
#include "cardano.h"
 
 
using namespace pappso;
 
MzCalibrationModel1::MzCalibrationModel1(double T1_frame,
                                         double T2_frame,
                                         double digitizerTimebase,
                                         double digitizerDelay,
                                         double C0,
                                         double C1,
                                         double C2,
                                         double C3,
                                         double C4,
                                         double T1_ref,
                                         double T2_ref,
                                         double dC1,
                                         double dC2)
  : MzCalibrationInterface(digitizerTimebase, digitizerDelay)
{
 
  double temperature_correction = dC1 * (T1_ref - T1_frame) + dC2 * (T2_ref - T2_frame);
  temperature_correction        = (double)1.0 + (temperature_correction / 1.0e6);
 
  // temperature compensation
  C1 = C1 * temperature_correction;
  C2 = C2 / temperature_correction;
 
 
  m_mzCalibrationArr.clear();
 
  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);
  m_mzCalibrationArr.push_back(C4);
}
 
MzCalibrationModel1::~MzCalibrationModel1()
{
}
 
double
MzCalibrationModel1::getMzFromTofIndex(quint32 tof_index)
{
  double tof = ((double)tof_index * m_digitizerTimebase) + m_digitizerDelay;
  // 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;
 
 
  /* transformation formula given by Bruker 29/8/2019 :
   * x = m + dm
   *
   * time = m_mzCalibrationArr[0]
   * + sqrt ((10^12)/m_mzCalibrationArr[1]) * x^0.5
   * + m_mzCalibrationArr[2] * x
   * + m_mzCalibrationArr[3] * x^1.5
   */
 
 
  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]);
      // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
    }
  else
    {
      // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
    }
 
  InHousePolynomialSolverResult res = inHousePolynomialSolve(X);
 
  if(res.type == CardanoResultCase::notvalid)
    {
      throw pappso::PappsoException(QObject::tr("ERROR in %1 %2 %3"
                                                "inHousePolynomialSolve :\nresult is not valid")
                                      .arg(__FILE__)
                                      .arg(__FUNCTION__)
                                      .arg(__LINE__));
    }
  return (pow(res.x1, 2) - m_mzCalibrationArr[4]);
}
 
quint32
MzCalibrationModel1::getTofIndexFromMz(double mz)
{
  // 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;
 
  mz = mz + m_mzCalibrationArr[4]; // mz_corr
 
  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);
}
 
pappso::MzCalibrationModel1Cached::MzCalibrationModel1Cached(double T1_frame,
                                                             double T2_frame,
                                                             double digitizerTimebase,
                                                             double digitizerDelay,
                                                             double C0,
                                                             double C1,
                                                             double C2,
                                                             double C3,
                                                             double C4,
                                                             double T1_ref,
                                                             double T2_ref,
                                                             double dC1,
                                                             double dC2)
  : MzCalibrationModel1(T1_frame,
                        T2_frame,
                        digitizerTimebase,
                        digitizerDelay,
                        C0,
                        C1,
                        C2,
                        C3,
                        C4,
                        T1_ref,
                        T2_ref,
                        dC1,
                        dC2)
{
}
 
pappso::MzCalibrationModel1Cached::~MzCalibrationModel1Cached()
{
}
 
 
double
MzCalibrationModel1Cached::getMzFromTofIndex(quint32 tof_index)
{
  if(m_max > tof_index)
    {
      if(m_arrMasses[tof_index] == 0)
        {
          m_arrMasses[tof_index] = MzCalibrationModel1::getMzFromTofIndex(tof_index);
        }
      return m_arrMasses[tof_index];
    }
  else
    {
      return MzCalibrationModel1::getMzFromTofIndex(tof_index);
    }
}