libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
mzcalibrationmodel1.cpp
Go to the documentation of this file.
1/**
2 * \file pappsomspp/vendors/tims/mzcalibration/mzcalibrationmodel1.cpp
3 * \date 11/11/2020
4 * \author Olivier Langella
5 * \brief implement Bruker's model type 1 formula to compute m/z
6 */
7
8/*******************************************************************************
9 * Copyright (c) 2020 Olivier Langella <Olivier.Langella@u-psud.fr>.
10 *
11 * This file is part of the PAPPSOms++ library.
12 *
13 * PAPPSOms++ is free software: you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation, either version 3 of the License, or
16 * (at your option) any later version.
17 *
18 * PAPPSOms++ is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with PAPPSOms++. If not, see <http://www.gnu.org/licenses/>.
25 *
26 ******************************************************************************/
27
28#include "mzcalibrationmodel1.h"
29#include <cmath>
30#include <QDebug>
31#include <QObject>
32#include "../../../pappsoexception.h"
33#include "cardano.h"
34
35
36using namespace pappso;
37
39 double T2_frame,
40 double digitizerTimebase,
41 double digitizerDelay,
42 double C0,
43 double C1,
44 double C2,
45 double C3,
46 double C4,
47 double T1_ref,
48 double T2_ref,
49 double dC1,
50 double dC2)
51 : MzCalibrationInterface(digitizerTimebase, digitizerDelay)
52{
53
54 double temperature_correction =
55 dC1 * (T1_ref - T1_frame) + dC2 * (T2_ref - T2_frame);
56 temperature_correction = (double)1.0 + (temperature_correction / 1.0e6);
57
58 // temperature compensation
59 C1 = C1 * temperature_correction;
60 C2 = C2 / temperature_correction;
61
62
63 m_mzCalibrationArr.clear();
64
65 m_digitizerDelay = digitizerDelay;
66 m_digitizerTimebase = digitizerTimebase;
67
68 m_mzCalibrationArr.push_back(C0);
69 m_mzCalibrationArr.push_back(std::sqrt(std::pow(10, 12) / C1));
70 m_mzCalibrationArr.push_back(C2);
71 m_mzCalibrationArr.push_back(C3);
72 m_mzCalibrationArr.push_back(C4);
73}
74
78
79double
81{
82 double tof = ((double)tof_index * m_digitizerTimebase) + m_digitizerDelay;
83 // http://www.alglib.net/equations/polynomial.php
84 // http://www.alglib.net/translator/man/manual.cpp.html#sub_polynomialsolve
85 // https://math.stackexchange.com/questions/1291208/number-of-roots-of-a-polynomial-of-non-integer-degree
86 // 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
87 // https://stackoverflow.com/questions/26091323/how-to-plot-a-function-curve-in-r
88 /*
89 * beware to put the function on a single line in R:
90> eq <- function(m){ 1 + (sqrt((10^12)/670) * sqrt(m)) + (207.775676931964 * m)
91+ (59.2526676368822 * (m^1.5)) }
92> eq <- function(m){ 313.577620892277 + (sqrt((10^12)/157424.07710945) *
93sqrt(m)) + (0.000338743021989553 * m)
94+ (0 * (m^1.5)) }
95> plot(eq(1:1000), type='l')
96
97
98
99> eq2 <- function(m2){ 1 + sqrt((10^12)/670) * m2 + 207.775676931964 * (m2^2)
100+ 59.2526676368822 * (m2^3) }
101> plot(eq2(1:sqrt(1000)), type='l')
102*/
103 // How to Factor a Trinomial with Fractions as Coefficients
104
105 // formula
106 // a = c0 = 1
107 // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
108 // c = c2, c2 = 207.775676931964 * m
109 // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
110 // double mz = 0;
111
112
113 /* transformation formula given by Bruker 29/8/2019 :
114 * x = m + dm
115 *
116 * time = m_mzCalibrationArr[0]
117 * + sqrt ((10^12)/m_mzCalibrationArr[1]) * x^0.5
118 * + m_mzCalibrationArr[2] * x
119 * + m_mzCalibrationArr[3] * x^1.5
120 */
121
122
123 std::vector<double> X;
124 X.push_back((m_mzCalibrationArr[0] - (double)tof));
125 X.push_back(m_mzCalibrationArr[1]);
126 if(m_mzCalibrationArr[2] != 0)
127 {
128 X.push_back(m_mzCalibrationArr[2]);
129 }
130 if(m_mzCalibrationArr[3] != 0)
131 {
132 X.push_back(m_mzCalibrationArr[3]);
133 // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
134 }
135 else
136 {
137 // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
138 }
139
141
143 {
145 QObject::tr("ERROR in %1 %2 %3"
146 "inHousePolynomialSolve :\nresult is not valid")
147 .arg(__FILE__)
148 .arg(__FUNCTION__)
149 .arg(__LINE__));
150 }
151 return (pow(res.x1, 2) - m_mzCalibrationArr[4]);
152}
153
154quint32
156{
157 // formula
158 // a = c0 = 1
159 // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
160 // c = c2, c2 = 207.775676931964 * m
161 // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
162 qDebug() << "mz=" << mz;
163
164 mz = mz + m_mzCalibrationArr[4]; // mz_corr
165
166 double tof = m_mzCalibrationArr[0];
167 qDebug() << "tof ( m_mzCalibrationArr[0])=" << tof;
168 // TODO cache value of std::sqrt((std::pow(10, 12) / m_mzCalibrationArr[1]))
169 tof += m_mzCalibrationArr[1] * std::sqrt(mz);
170 qDebug() << "tof=" << tof;
171 tof += m_mzCalibrationArr[2] * mz;
172 qDebug() << "tof=" << tof;
173 tof += m_mzCalibrationArr[3] * std::pow(mz, 1.5);
174 qDebug() << "tof=" << tof;
175 tof -= m_digitizerDelay;
176 qDebug() << "tof=" << tof;
177 tof = tof / m_digitizerTimebase;
178 qDebug() << "index=" << tof;
179 return (quint32)std::round(tof);
180}
181
183 double T1_frame,
184 double T2_frame,
185 double digitizerTimebase,
186 double digitizerDelay,
187 double C0,
188 double C1,
189 double C2,
190 double C3,
191 double C4,
192 double T1_ref,
193 double T2_ref,
194 double dC1,
195 double dC2)
196 : MzCalibrationModel1(T1_frame,
197 T2_frame,
198 digitizerTimebase,
199 digitizerDelay,
200 C0,
201 C1,
202 C2,
203 C3,
204 C4,
205 T1_ref,
206 T2_ref,
207 dC1,
208 dC2)
209{
210}
211
215
216
217double
219{
220 if(m_max > tof_index)
221 {
222 if(m_arrMasses[tof_index] == 0)
223 {
224 m_arrMasses[tof_index] =
226 }
227 return m_arrMasses[tof_index];
228 }
229 else
230 {
232 }
233}
InHousePolynomialSolverResult inHousePolynomialSolve(const std::vector< double > &polynome)
Definition cardano.cpp:118
cubic solver adapted from https://www.codeproject.com/articles/798474/to-solve-a-cubic-equation thank...
std::vector< double > m_mzCalibrationArr
MZ calibration parameters.
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)
virtual double getMzFromTofIndex(quint32 tof_index) override
get m/z from time of flight raw index
virtual double getMzFromTofIndex(quint32 tof_index) override
get m/z from time of flight raw index
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)
virtual quint32 getTofIndexFromMz(double mz) override
get raw TOF index of a given m/z
implement Bruker's model type 1 formula to compute m/z
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39