libpappsomspp
Library for mass spectrometry
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 <solvers.h>
30#include <cmath>
31#include <QDebug>
32#include "../../../pappsoexception.h"
33
34
35using namespace pappso;
36
37MzCalibrationModel1::MzCalibrationModel1(double T1_frame,
38 double T2_frame,
39 double digitizerTimebase,
40 double digitizerDelay,
41 double C0,
42 double C1,
43 double C2,
44 double C3,
45 double C4,
46 double T1_ref,
47 double T2_ref,
48 double dC1,
49 double dC2)
50 : MzCalibrationInterface(digitizerTimebase, digitizerDelay)
51{
52
53 double temperature_correction =
54 dC1 * (T1_ref - T1_frame) + dC2 * (T2_ref - T2_frame);
55 temperature_correction = (double)1.0 + (temperature_correction / 1.0e6);
56
57 // temperature compensation
58 C1 = C1 * temperature_correction;
59 C2 = C2 / temperature_correction;
60
61
62 m_mzCalibrationArr.clear();
63
64 m_digitizerDelay = digitizerDelay;
65 m_digitizerTimebase = digitizerTimebase;
66
67 m_mzCalibrationArr.push_back(C0);
68 m_mzCalibrationArr.push_back(std::sqrt(std::pow(10, 12) / C1));
69 m_mzCalibrationArr.push_back(C2);
70 m_mzCalibrationArr.push_back(C3);
71 m_mzCalibrationArr.push_back(C4);
72}
73
75{
76}
77
78double
80{
81 double tof = ((double)tof_index * m_digitizerTimebase) + m_digitizerDelay;
82 // http://www.alglib.net/equations/polynomial.php
83 // http://www.alglib.net/translator/man/manual.cpp.html#sub_polynomialsolve
84 // https://math.stackexchange.com/questions/1291208/number-of-roots-of-a-polynomial-of-non-integer-degree
85 // 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
86 // https://stackoverflow.com/questions/26091323/how-to-plot-a-function-curve-in-r
87 /*
88 * beware to put the function on a single line in R:
89> eq <- function(m){ 1 + (sqrt((10^12)/670) * sqrt(m)) + (207.775676931964 * m)
90+ (59.2526676368822 * (m^1.5)) }
91> eq <- function(m){ 313.577620892277 + (sqrt((10^12)/157424.07710945) *
92sqrt(m)) + (0.000338743021989553 * m)
93+ (0 * (m^1.5)) }
94> plot(eq(1:1000), type='l')
95
96
97
98> eq2 <- function(m2){ 1 + sqrt((10^12)/670) * m2 + 207.775676931964 * (m2^2)
99+ 59.2526676368822 * (m2^3) }
100> plot(eq2(1:sqrt(1000)), type='l')
101*/
102 // How to Factor a Trinomial with Fractions as Coefficients
103
104 // formula
105 // a = c0 = 1
106 // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
107 // c = c2, c2 = 207.775676931964 * m
108 // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
109 // double mz = 0;
110
111
112 /* transformation formula given by Bruker 29/8/2019 :
113 * x = m + dm
114 *
115 * time = m_mzCalibrationArr[0]
116 * + sqrt ((10^12)/m_mzCalibrationArr[1]) * x^0.5
117 * + m_mzCalibrationArr[2] * x
118 * + m_mzCalibrationArr[3] * x^1.5
119 */
120 std::vector<double> X;
121 X.push_back((m_mzCalibrationArr[0] - (double)tof));
122 X.push_back(m_mzCalibrationArr[1]);
123 if(m_mzCalibrationArr[2] != 0)
124 {
125 X.push_back(m_mzCalibrationArr[2]);
126 }
127 if(m_mzCalibrationArr[3] != 0)
128 {
129 X.push_back(m_mzCalibrationArr[3]);
130 // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
131 }
132 else
133 {
134 // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
135 }
136 // qDebug() << "polynom_array :";
137 /*
138 for(double arg : X)
139 {
140 qDebug() << arg;
141 }
142 */
143 alglib::real_1d_array polynom_array;
144 try
145 {
146 polynom_array.setcontent(X.size(), &(X[0]));
147 }
148 catch(alglib::ap_error &error)
149 {
150 // PolynomialSolve: A[N]=0
152 QObject::tr("ERROR in alglib::polynom_array.setcontent :\n%1")
153 .arg(error.msg.c_str()));
154 }
155
156
157 /*
158 alglib::polynomialsolve(
159real_1d_array a,
160ae_int_t n,
161complex_1d_array& x,
162polynomialsolverreport& rep,
163const xparams _params = alglib::xdefault);
164*/
165 alglib::complex_1d_array m;
166 alglib::polynomialsolverreport rep;
167 // qDebug();
168 try
169 {
170 alglib::polynomialsolve(polynom_array, X.size() - 1, m, rep);
171 }
172 catch(alglib::ap_error &error)
173 {
174 qDebug() << " X.size() - 1 = " << X.size() - 1;
175 qDebug() << m_mzCalibrationArr[0];
176 qDebug() << m_mzCalibrationArr[1];
177 qDebug() << m_mzCalibrationArr[2];
178 qDebug() << m_mzCalibrationArr[3];
179
180 // PolynomialSolve: A[N]=0
182 QObject::tr("ERROR in MzCalibrationModel1::getMzFromTofIndex, "
183 "alglib::polynomialsolve :\n%1")
184 .arg(error.msg.c_str()));
185 }
186
187
188 // qDebug();
189
190 if(m.length() == 0)
191 {
192 throw pappso::PappsoException(QObject::tr(
193 "ERROR in MzCalibrationModel1::getMzFromTofIndex m.size() == 0"));
194 }
195 // qDebug();
196 if(m[0].y != 0)
197 {
199 QObject::tr("ERROR in MzCalibrationModel1::getMzFromTofIndex m[0].y!= "
200 "0 for tof index=%1")
201 .arg(tof_index));
202 }
203
204 // qDebug() << "m.length()=" << m.length();
205 // qDebug() << "m1=" << pow(m[0].x, 2);
206 // qDebug() << "m2=" << pow(m[1].x, 2);
207 return (pow(m[0].x, 2) - m_mzCalibrationArr[4]);
208}
209
210quint32
212{
213 // formula
214 // a = c0 = 1
215 // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
216 // c = c2, c2 = 207.775676931964 * m
217 // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
218 qDebug() << "mz=" << mz;
219
220 mz = mz + m_mzCalibrationArr[4]; // mz_corr
221
222 double tof = m_mzCalibrationArr[0];
223 qDebug() << "tof ( m_mzCalibrationArr[0])=" << tof;
224 // TODO cache value of std::sqrt((std::pow(10, 12) / m_mzCalibrationArr[1]))
225 tof += m_mzCalibrationArr[1] * std::sqrt(mz);
226 qDebug() << "tof=" << tof;
227 tof += m_mzCalibrationArr[2] * mz;
228 qDebug() << "tof=" << tof;
229 tof += m_mzCalibrationArr[3] * std::pow(mz, 1.5);
230 qDebug() << "tof=" << tof;
231 tof -= m_digitizerDelay;
232 qDebug() << "tof=" << tof;
233 tof = tof / m_digitizerTimebase;
234 qDebug() << "index=" << tof;
235 return (quint32)std::round(tof);
236}
237
239 double T1_frame,
240 double T2_frame,
241 double digitizerTimebase,
242 double digitizerDelay,
243 double C0,
244 double C1,
245 double C2,
246 double C3,
247 double C4,
248 double T1_ref,
249 double T2_ref,
250 double dC1,
251 double dC2)
252 : MzCalibrationModel1(T1_frame,
253 T2_frame,
254 digitizerTimebase,
255 digitizerDelay,
256 C0,
257 C1,
258 C2,
259 C3,
260 C4,
261 T1_ref,
262 T2_ref,
263 dC1,
264 dC2)
265{
266}
267
269{
270}
271
272
273double
275{
276 if(m_max > tof_index)
277 {
278 if(m_arrMasses[tof_index] == 0)
279 {
280 m_arrMasses[tof_index] =
282 }
283 return m_arrMasses[tof_index];
284 }
285 else
286 {
288 }
289}
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
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