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 <QObject>
33 #include "../../../pappsoexception.h"
34 
35 
36 using 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 
76 {
77 }
78 
79 double
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) *
93 sqrt(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  std::vector<double> X;
122  X.push_back((m_mzCalibrationArr[0] - (double)tof));
123  X.push_back(m_mzCalibrationArr[1]);
124  if(m_mzCalibrationArr[2] != 0)
125  {
126  X.push_back(m_mzCalibrationArr[2]);
127  }
128  if(m_mzCalibrationArr[3] != 0)
129  {
130  X.push_back(m_mzCalibrationArr[3]);
131  // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
132  }
133  else
134  {
135  // qDebug() << "m_mzCalibrationArr[3]=" << m_mzCalibrationArr[3];
136  }
137  // qDebug() << "polynom_array :";
138  /*
139  for(double arg : X)
140  {
141  qDebug() << arg;
142  }
143  */
144  alglib::real_1d_array polynom_array;
145  try
146  {
147  polynom_array.setcontent(X.size(), &(X[0]));
148  }
149  catch(alglib::ap_error &error)
150  {
151  // PolynomialSolve: A[N]=0
153  QObject::tr("ERROR in alglib::polynom_array.setcontent :\n%1")
154  .arg(error.msg.c_str()));
155  }
156 
157 
158  /*
159  alglib::polynomialsolve(
160 real_1d_array a,
161 ae_int_t n,
162 complex_1d_array& x,
163 polynomialsolverreport& rep,
164 const xparams _params = alglib::xdefault);
165 */
166  alglib::complex_1d_array m;
167  alglib::polynomialsolverreport rep;
168 
169  alglib::xparams params = alglib::xdefault;
170  // qDebug();
171  try
172  {
173  alglib::polynomialsolve(polynom_array, X.size() - 1, m, rep, params);
174  }
175  catch(alglib::ap_error &error)
176  {
177  qDebug() << " X.size() - 1 = " << X.size() - 1;
178  qDebug() << m_mzCalibrationArr[0];
179  qDebug() << m_mzCalibrationArr[1];
180  qDebug() << m_mzCalibrationArr[2];
181  qDebug() << m_mzCalibrationArr[3];
182 
183  // PolynomialSolve: A[N]=0
185  QObject::tr("ERROR in MzCalibrationModel1::getMzFromTofIndex, "
186  "alglib::polynomialsolve :\n%1")
187  .arg(error.msg.c_str()));
188  }
189 
190 
191  // qDebug();
192 
193  if(m.length() == 0)
194  {
195  throw pappso::PappsoException(QObject::tr(
196  "ERROR in MzCalibrationModel1::getMzFromTofIndex m.size() == 0"));
197  }
198  // qDebug();
199  if(m[0].y != 0)
200  {
202  QObject::tr("ERROR in MzCalibrationModel1::getMzFromTofIndex m[0].y!= "
203  "0 for tof index=%1")
204  .arg(tof_index));
205  }
206 
207  // qDebug() << "m.length()=" << m.length();
208  // qDebug() << "m1=" << pow(m[0].x, 2);
209  // qDebug() << "m2=" << pow(m[1].x, 2);
210  return (pow(m[0].x, 2) - m_mzCalibrationArr[4]);
211 }
212 
213 quint32
215 {
216  // formula
217  // a = c0 = 1
218  // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
219  // c = c2, c2 = 207.775676931964 * m
220  // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
221  qDebug() << "mz=" << mz;
222 
223  mz = mz + m_mzCalibrationArr[4]; // mz_corr
224 
225  double tof = m_mzCalibrationArr[0];
226  qDebug() << "tof ( m_mzCalibrationArr[0])=" << tof;
227  // TODO cache value of std::sqrt((std::pow(10, 12) / m_mzCalibrationArr[1]))
228  tof += m_mzCalibrationArr[1] * std::sqrt(mz);
229  qDebug() << "tof=" << tof;
230  tof += m_mzCalibrationArr[2] * mz;
231  qDebug() << "tof=" << tof;
232  tof += m_mzCalibrationArr[3] * std::pow(mz, 1.5);
233  qDebug() << "tof=" << tof;
234  tof -= m_digitizerDelay;
235  qDebug() << "tof=" << tof;
236  tof = tof / m_digitizerTimebase;
237  qDebug() << "index=" << tof;
238  return (quint32)std::round(tof);
239 }
240 
242  double T1_frame,
243  double T2_frame,
244  double digitizerTimebase,
245  double digitizerDelay,
246  double C0,
247  double C1,
248  double C2,
249  double C3,
250  double C4,
251  double T1_ref,
252  double T2_ref,
253  double dC1,
254  double dC2)
255  : MzCalibrationModel1(T1_frame,
256  T2_frame,
257  digitizerTimebase,
258  digitizerDelay,
259  C0,
260  C1,
261  C2,
262  C3,
263  C4,
264  T1_ref,
265  T2_ref,
266  dC1,
267  dC2)
268 {
269 }
270 
272 {
273 }
274 
275 
276 double
278 {
279  if(m_max > tof_index)
280  {
281  if(m_arrMasses[tof_index] == 0)
282  {
283  m_arrMasses[tof_index] =
285  }
286  return m_arrMasses[tof_index];
287  }
288  else
289  {
290  return MzCalibrationModel1::getMzFromTofIndex(tof_index);
291  }
292 }
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