libpappsomspp
Library for mass spectrometry
timsframebase.cpp
Go to the documentation of this file.
1 /**
2  * \file pappsomspp/vendors/tims/timsframebase.cpp
3  * \date 16/12/2019
4  * \author Olivier Langella
5  * \brief handle a single Bruker's TimsTof frame without binary data
6  */
7 
8 /*******************************************************************************
9  * Copyright (c) 2019 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 #include "timsframebase.h"
28 #include "../../../pappsomspp/pappsoexception.h"
29 #include "../../../pappsomspp/exception/exceptionoutofrange.h"
30 #include <QDebug>
31 #include <cmath>
32 #include <solvers.h>
33 
34 namespace pappso
35 {
36 
37 TimsFrameBase::TimsFrameBase(std::size_t timsId, quint32 scanNum)
38 {
39  qDebug() << timsId;
40  m_timsId = timsId;
41 
42  m_scanNumber = scanNum;
43 }
44 
45 TimsFrameBase::TimsFrameBase([[maybe_unused]] const TimsFrameBase &other)
46 {
47 }
48 
50 {
51 }
52 
53 void
54 TimsFrameBase::setAccumulationTime(double accumulation_time_ms)
55 {
56  m_accumulationTime = accumulation_time_ms;
57 }
58 
59 
60 void
61 TimsFrameBase::setMzCalibration(double temperature_correction,
62  double digitizerTimebase,
63  double digitizerDelay,
64  double C0,
65  double C1,
66  double C2,
67  double C3)
68 {
69 
70  // temperature compensation
71  C1 = C1 * temperature_correction;
72  C2 = C2 / temperature_correction;
73 
74  m_digitizerDelay = digitizerDelay;
75  m_digitizerTimebase = digitizerTimebase;
76 
77  m_mzCalibrationArr.push_back(C0);
78  m_mzCalibrationArr.push_back(std::sqrt(std::pow(10, 12) / C1));
79  m_mzCalibrationArr.push_back(C2);
80  m_mzCalibrationArr.push_back(C3);
81 }
82 
83 bool
84 TimsFrameBase::checkScanNum(std::size_t scanNum) const
85 {
86  if(scanNum >= m_scanNumber)
87  {
89  QObject::tr("Invalid scan number : scanNum%1 > m_scanNumber")
90  .arg(scanNum));
91  }
92 
93  return true;
94 }
95 
96 std::size_t
97 TimsFrameBase::getNbrPeaks(std::size_t scanNum) const
98 {
99  throw PappsoException(
100  QObject::tr(
101  "ERROR unable to get number of peaks in TimsFrameBase for scan number %1")
102  .arg(scanNum));
103 }
104 
106 TimsFrameBase::getMassSpectrumSPtr(std::size_t scanNum) const
107 {
108  throw PappsoException(
109  QObject::tr(
110  "ERROR unable to getMassSpectrumSPtr in TimsFrameBase for scan number %1")
111  .arg(scanNum));
112 }
113 Trace
114 TimsFrameBase::cumulateScanToTrace(std::size_t scanNumBegin,
115  std::size_t scanNumEnd) const
116 {
117  throw PappsoException(
118  QObject::tr("ERROR unable to cumulateScanToTrace in TimsFrameBase for scan "
119  "number begin %1 end %2")
120  .arg(scanNumBegin)
121  .arg(scanNumEnd));
122 }
123 double
125 {
126  // mz calibration
127  return (index * m_digitizerTimebase) + m_digitizerDelay;
128 }
129 double
130 TimsFrameBase::getTofFromIndex(quint32 index) const
131 {
132  // mz calibration
133  return ((double)index * m_digitizerTimebase) + m_digitizerDelay;
134 }
135 double
136 TimsFrameBase::getMzFromTof(double tof) const
137 {
138  // http://www.alglib.net/equations/polynomial.php
139  // http://www.alglib.net/translator/man/manual.cpp.html#sub_polynomialsolve
140  // https://math.stackexchange.com/questions/1291208/number-of-roots-of-a-polynomial-of-non-integer-degree
141  // 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
142  // https://stackoverflow.com/questions/26091323/how-to-plot-a-function-curve-in-r
143  /*
144  * beware to put the function on a single line in R:
145 > eq <- function(m){ 1 + (sqrt((10^12)/670) * sqrt(m)) + (207.775676931964 * m)
146 + (59.2526676368822 * (m^1.5)) }
147 > eq <- function(m){ 313.577620892277 + (sqrt((10^12)/157424.07710945) *
148 sqrt(m)) + (0.000338743021989553 * m)
149 + (0 * (m^1.5)) }
150 > plot(eq(1:1000), type='l')
151 
152 
153 
154 > eq2 <- function(m2){ 1 + sqrt((10^12)/670) * m2 + 207.775676931964 * (m2^2)
155 + 59.2526676368822 * (m2^3) }
156 > plot(eq2(1:sqrt(1000)), type='l')
157 */
158  // How to Factor a Trinomial with Fractions as Coefficients
159 
160  // formula
161  // a = c0 = 1
162  // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
163  // c = c2, c2 = 207.775676931964 * m
164  // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
165  // double mz = 0;
166  std::vector<double> X;
167  X.push_back((m_mzCalibrationArr[0] - (double)tof));
168  X.push_back(m_mzCalibrationArr[1]);
169  if(m_mzCalibrationArr[2] != 0)
170  {
171  X.push_back(m_mzCalibrationArr[2]);
172  }
173  if(m_mzCalibrationArr[3] != 0)
174  {
175  X.push_back(m_mzCalibrationArr[3]);
176  }
177  else
178  {
179  qDebug() << m_mzCalibrationArr[3];
180  }
181 
182  alglib::real_1d_array polynom_array;
183  try
184  {
185  polynom_array.setcontent(X.size(), &(X[0]));
186  }
187  catch(alglib::ap_error &error)
188  {
189  // PolynomialSolve: A[N]=0
191  QObject::tr("ERROR in alglib::polynom_array.setcontent :\n%1")
192  .arg(error.msg.c_str()));
193  }
194 
195 
196  /*
197  alglib::polynomialsolve(
198 real_1d_array a,
199 ae_int_t n,
200 complex_1d_array& x,
201 polynomialsolverreport& rep,
202 const xparams _params = alglib::xdefault);
203 */
204  alglib::complex_1d_array m;
205  alglib::polynomialsolverreport rep;
206  // qDebug();
207  try
208  {
209  alglib::polynomialsolve(polynom_array, X.size() - 1, m, rep);
210  }
211  catch(alglib::ap_error &error)
212  {
213  qDebug() << " X.size() - 1 = " << X.size() - 1;
214  qDebug() << m_mzCalibrationArr[0];
215  qDebug() << m_mzCalibrationArr[1];
216  qDebug() << m_mzCalibrationArr[2];
217  qDebug() << m_mzCalibrationArr[3];
218 
219  // PolynomialSolve: A[N]=0
221  QObject::tr("ERROR in alglib::polynomialsolve :\n%1")
222  .arg(error.msg.c_str()));
223  }
224 
225 
226  // qDebug();
227 
228  if(m.length() == 0)
229  {
231  QObject::tr("ERROR in TimsFrame::getMzFromTof m.size() == 0"));
232  }
233  // qDebug();
234  if(m[0].y != 0)
235  {
237  QObject::tr("ERROR in TimsFrame::getMzFromTof m[0].y!= 0"));
238  }
239 
240  return pow(m[0].x, 2);
241 }
242 
243 quint32
244 TimsFrameBase::getRawIndexFromMz(double mz) const
245 {
246  // formula
247  // a = c0 = 1
248  // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
249  // c = c2, c2 = 207.775676931964 * m
250  // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
251  qDebug() << "mz=" << mz;
252 
253  double tof = m_mzCalibrationArr[0];
254  qDebug() << "tof ( m_mzCalibrationArr[0])=" << tof;
255  // TODO cache value of std::sqrt((std::pow(10, 12) / m_mzCalibrationArr[1]))
256  tof += m_mzCalibrationArr[1] * std::sqrt(mz);
257  qDebug() << "tof=" << tof;
258  tof += m_mzCalibrationArr[2] * mz;
259  qDebug() << "tof=" << tof;
260  tof += m_mzCalibrationArr[3] * std::pow(mz, 1.5);
261  qDebug() << "tof=" << tof;
263  qDebug() << "tof=" << tof;
264  tof = tof / m_digitizerTimebase;
265  qDebug() << "index=" << tof;
266  return (quint32)std::round(tof);
267 }
268 
269 void
270 TimsFrameBase::setTime(double time)
271 {
272  m_time = time;
273 }
274 
275 void
276 TimsFrameBase::setMsMsType(quint8 type)
277 {
278 
279  qDebug() << " m_msMsType=" << type;
280  m_msMsType = type;
281 }
282 
283 unsigned int
285 {
286  if(m_msMsType == 0)
287  return 1;
288  return 2;
289 }
290 
291 double
293 {
294  return m_time;
295 }
296 
297 std::size_t
298 TimsFrameBase::getId() const
299 {
300  return m_timsId;
301 }
302 void
303 TimsFrameBase::setTimsCalibration(int tims_model_type,
304  double C0,
305  double C1,
306  double C2,
307  double C3,
308  double C4,
309  [[maybe_unused]] double C5,
310  double C6,
311  double C7,
312  double C8,
313  double C9)
314 {
315  if(tims_model_type != 2)
316  {
317  throw pappso::PappsoException(QObject::tr(
318  "ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
319  }
320  m_timsDvStart = C2; // C2 from TimsCalibration
321  m_timsTtrans = C4; // C4 from TimsCalibration
322  m_timsNdelay = C0; // C0 from TimsCalibration
323  m_timsVmin = C8; // C8 from TimsCalibration
324  m_timsVmax = C9; // C9 from TimsCalibration
325  m_timsC6 = C6;
326  m_timsC7 = C7;
327 
328 
329  m_timsSlope =
330  (C3 - m_timsDvStart) / C1; // //C3 from TimsCalibration // C2 from
331  // TimsCalibration // C1 from TimsCalibration
332 }
333 double
334 TimsFrameBase::getVoltageTransformation(std::size_t scanNum) const
335 {
336  double v = m_timsDvStart +
337  m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
338 
339  if(v < m_timsVmin)
340  {
342  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
343  "calibration, v < m_timsVmin"));
344  }
345 
346 
347  if(v > m_timsVmax)
348  {
350  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
351  "calibration, v > m_timsVmax"));
352  }
353  return v;
354 }
355 double
356 TimsFrameBase::getDriftTime(std::size_t scanNum) const
357 {
358  return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
359 }
360 
361 double
362 TimsFrameBase::getOneOverK0Transformation(std::size_t scanNum) const
363 {
364  return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
365 }
366 
367 
368 std::size_t
369 TimsFrameBase::getScanNumFromOneOverK0(double one_over_k0) const
370 {
371  double temp = 1 / one_over_k0;
372  temp = temp - m_timsC6;
373  temp = temp / m_timsC7;
374  temp = 1 / temp;
375  temp = temp - m_timsDvStart;
376  temp = temp / m_timsSlope + m_timsTtrans + m_timsNdelay;
377  return (std::size_t)std::round(temp);
378 }
379 
380 } // namespace pappso
pappso::TimsFrameBase::m_digitizerTimebase
double m_digitizerTimebase
Definition: timsframebase.h:178
pappso::TimsFrameBase::checkScanNum
bool checkScanNum(std::size_t scanNum) const
Definition: timsframebase.cpp:102
pappso::TimsFrameBase::m_timsId
std::size_t m_timsId
Tims frame database id (the SQL identifier of this frame)
Definition: timsframebase.h:172
pappso::TimsFrameBase::setTimsCalibration
void setTimsCalibration(int tims_model_type, double C0, double C1, double C2, double C3, double C4, double C5, double C6, double C7, double C8, double C9)
Definition: timsframebase.cpp:321
pappso::TimsFrameBase::getMsLevel
unsigned int getMsLevel() const
Definition: timsframebase.cpp:302
pappso::TimsFrameBase::cumulateScanToTrace
virtual Trace cumulateScanToTrace(std::size_t scanNumBegin, std::size_t scanNumEnd) const
Definition: timsframebase.cpp:132
pappso::TimsFrameBase::m_timsSlope
double m_timsSlope
Definition: timsframebase.h:192
pappso::TimsFrameBase::m_timsC7
double m_timsC7
Definition: timsframebase.h:200
pappso
Definition: aa.cpp:38
pappso::TimsFrameBase::setAccumulationTime
void setAccumulationTime(double accumulation_time_ms)
Definition: timsframebase.cpp:72
pappso::TimsFrameBase::setMzCalibration
void setMzCalibration(double temperature_correction, double digitizerTimebase, double digitizerDelay, double C0, double C1, double C2, double C3)
Definition: timsframebase.cpp:79
pappso::TimsFrameBase::setTime
void setTime(double time)
Definition: timsframebase.cpp:288
pappso::TimsFrameBase::getNbrPeaks
virtual std::size_t getNbrPeaks(std::size_t scanNum) const
Definition: timsframebase.cpp:115
pappso::TimsFrameBase::getRawIndexFromMz
quint32 getRawIndexFromMz(double mz) const
get raw index of a given m/z
Definition: timsframebase.cpp:262
pappso::TimsFrameBase::getDriftTime
double getDriftTime(std::size_t scanNum) const
get drift time of a scan number in milliseconds
Definition: timsframebase.cpp:374
pappso::TimsFrameBase::m_accumulationTime
double m_accumulationTime
accumulation time in milliseconds
Definition: timsframebase.h:176
pappso::TimsFrameBase::m_time
double m_time
retention time
Definition: timsframebase.h:189
pappso::ExceptionOutOfRange
Definition: exceptionoutofrange.h:52
pappso::TimsFrameBase::getVoltageTransformation
double getVoltageTransformation(std::size_t scanNum) const
Definition: timsframebase.cpp:352
pappso::TimsFrameBase::getId
std::size_t getId() const
Definition: timsframebase.cpp:316
pappso::TimsFrameBase::m_timsDvStart
double m_timsDvStart
Definition: timsframebase.h:191
pappso::TimsFrameBase::getMassSpectrumSPtr
virtual MassSpectrumSPtr getMassSpectrumSPtr(std::size_t scanNum) const
Definition: timsframebase.cpp:124
pappso::TimsFrameBase::getScanNumFromOneOverK0
std::size_t getScanNumFromOneOverK0(double one_over_k0) const
get the scan number from a given 1/Ko mobility value
Definition: timsframebase.cpp:387
timsframebase.h
handle a single Bruker's TimsTof frame without binary data
pappso::TimsFrameBase::~TimsFrameBase
~TimsFrameBase()
Definition: timsframebase.cpp:67
pappso::TimsFrameBase::TimsFrameBase
TimsFrameBase(std::size_t timsId, quint32 scanNum)
constructor for binary independant tims frame
Definition: timsframebase.cpp:55
pappso::TimsFrameBase::setMsMsType
void setMsMsType(quint8 type)
Definition: timsframebase.cpp:294
pappso::TimsFrameBase::m_timsNdelay
double m_timsNdelay
Definition: timsframebase.h:196
pappso::TimsFrameBase::m_timsVmin
double m_timsVmin
Definition: timsframebase.h:197
pappso::TimsFrameBase::getOneOverK0Transformation
double getOneOverK0Transformation(std::size_t scanNum) const
get 1/K0 value of a given scan (mobility value)
Definition: timsframebase.cpp:380
pappso::PrecisionUnit::mz
@ mz
pappso::TimsFrameBase::getMzFromTof
double getMzFromTof(double tof) const
get m/z from time of flight
Definition: timsframebase.cpp:154
pappso::TimsFrameBase::m_scanNumber
quint32 m_scanNumber
total number of scans contained in this frame
Definition: timsframebase.h:166
pappso::TimsFrameBase::m_timsC6
double m_timsC6
Definition: timsframebase.h:199
pappso::TimsFrameBase::getTofFromIndex
double getTofFromIndex(quint32 index) const
get time of flight from raw index
Definition: timsframebase.cpp:148
pappso::TimsFrameBase::m_digitizerDelay
double m_digitizerDelay
Definition: timsframebase.h:179
pappso::TimsFrameBase::m_mzCalibrationArr
std::vector< double > m_mzCalibrationArr
MZ calibration parameters.
Definition: timsframebase.h:183
pappso::TimsFrameBase::m_msMsType
quint8 m_msMsType
Definition: timsframebase.h:185
pappso::TimsFrameBase::m_timsVmax
double m_timsVmax
Definition: timsframebase.h:198
pappso::TimsFrameBase::m_timsTtrans
double m_timsTtrans
Definition: timsframebase.h:195
pappso::TimsFrameBase::getTime
double getTime() const
Definition: timsframebase.cpp:310
pappso::MassSpectrumSPtr
std::shared_ptr< MassSpectrum > MassSpectrumSPtr
Definition: massspectrum.h:75
pappso::PappsoException
Definition: pappsoexception.h:62