Loading...
Searching...
No Matches
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);
}
}
inHousePolynomialSolve
InHousePolynomialSolverResult inHousePolynomialSolve(const std::vector< double > &polynome)
Definition cardano.cpp:115
cardano.h
CardanoResultCase::notvalid
@ notvalid
Definition cardano.h:37
pappso::MzCalibrationInterface::MzCalibrationInterface
MzCalibrationInterface(double digitizerTimebase, double digitizerDelay)
Definition mzcalibrationinterface.cpp:34
pappso::MzCalibrationInterface::m_mzCalibrationArr
std::vector< double > m_mzCalibrationArr
MZ calibration parameters.
Definition mzcalibrationinterface.h:101
pappso::MzCalibrationInterface::m_digitizerTimebase
double m_digitizerTimebase
Definition mzcalibrationinterface.h:96
pappso::MzCalibrationInterface::m_digitizerDelay
double m_digitizerDelay
Definition mzcalibrationinterface.h:97
pappso::MzCalibrationModel1Cached::m_max
quint32 m_max
Definition mzcalibrationmodel1.h:105
pappso::MzCalibrationModel1Cached::~MzCalibrationModel1Cached
virtual ~MzCalibrationModel1Cached()
Definition mzcalibrationmodel1.cpp:209
pappso::MzCalibrationModel1Cached::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)
Definition mzcalibrationmodel1.cpp:180
pappso::MzCalibrationModel1Cached::m_arrMasses
double m_arrMasses[600000]
Definition mzcalibrationmodel1.h:104
pappso::MzCalibrationModel1Cached::getMzFromTofIndex
virtual double getMzFromTofIndex(quint32 tof_index) override
get m/z from time of flight raw index
Definition mzcalibrationmodel1.cpp:215
pappso::MzCalibrationModel1::getMzFromTofIndex
virtual double getMzFromTofIndex(quint32 tof_index) override
get m/z from time of flight raw index
Definition mzcalibrationmodel1.cpp:79
pappso::MzCalibrationModel1::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)
Definition mzcalibrationmodel1.cpp:38
pappso::MzCalibrationModel1::getTofIndexFromMz
virtual quint32 getTofIndexFromMz(double mz) override
get raw TOF index of a given m/z
Definition mzcalibrationmodel1.cpp:153
pappso::MzCalibrationModel1::~MzCalibrationModel1
virtual ~MzCalibrationModel1()
Definition mzcalibrationmodel1.cpp:74
pappso::PappsoException
Definition pappsoexception.h:41
mzcalibrationmodel1.h
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
pappsoexception.h
InHousePolynomialSolverResult
Definition cardano.h:47