corems.mass_spectrum.calc.Calibration

  1# -*- coding: utf-8 -*-
  2"""
  3Created on Wed May 13 02:16:09 2020
  4
  5@author: Will Kew
  6"""
  7
  8# import modules
  9import csv
 10import warnings
 11from io import BytesIO, StringIO
 12from pathlib import Path
 13
 14import numpy as np
 15import pandas as pd
 16from s3path import S3Path
 17
 18# import scipy modules for calibration
 19from scipy.optimize import minimize
 20
 21
 22class MzDomainCalibration:
 23    """MzDomainCalibration class for recalibrating mass spectra
 24
 25    Parameters
 26    ----------
 27    mass_spectrum : CoreMS MassSpectrum Object
 28        The mass spectrum to be calibrated.
 29    ref_masslist : str
 30        The path to a reference mass list.
 31    mzsegment : tuple of floats, optional
 32        The mz range to recalibrate, or None. Used for calibration of specific parts of the mz domain at a time.
 33        Future work - allow multiple mzsegments to be passed.
 34
 35    Attributes
 36    ----------
 37    mass_spectrum : CoreMS MassSpectrum Object
 38        The mass spectrum to be calibrated.
 39    mzsegment : tuple of floats or None
 40        The mz range to recalibrate, or None.
 41    ref_mass_list_path : str or Path
 42        The path to the reference mass list.
 43
 44    Methods
 45    -------
 46    * run().
 47        Main function to run this class.
 48    * load_ref_mass_list().
 49        Load reference mass list (Bruker format).
 50    * gen_ref_mass_list_from_assigned(min_conf=0.7).
 51        Generate reference mass list from assigned masses.
 52    * find_calibration_points(df_ref, calib_ppm_error_threshold=(-1, 1), calib_snr_threshold=5).
 53        Find calibration points in the mass spectrum based on the reference mass list.
 54    * robust_calib(param, cal_peaks_mz, cal_refs_mz, order=1).
 55        Recalibration function.
 56    * recalibrate_mass_spectrum(cal_peaks_mz, cal_refs_mz, order=1, diagnostic=False).
 57        Main recalibration function which uses a robust linear regression.
 58
 59
 60    """
 61
 62    def __init__(self, mass_spectrum, ref_masslist, mzsegment=None):
 63        self.mass_spectrum = mass_spectrum
 64        self.mzsegment = mzsegment
 65
 66        # define reference mass list - bruker .ref format
 67        self.ref_mass_list_path = ref_masslist
 68        if self.mass_spectrum.percentile_assigned(mute_output=True)[0] != 0:
 69            warnings.warn(
 70                "Warning: calibrating spectra which have already been assigned may yield erroneous results"
 71            )
 72        self.mass_spectrum.mz_cal = self.mass_spectrum.mz_exp
 73        self.mass_spectrum.mz_cal_profile = self.mass_spectrum._mz_exp
 74
 75        if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
 76            print(
 77                "MS Obj loaded - " + str(len(mass_spectrum.mspeaks)) + " peaks found."
 78            )
 79     
 80    def load_ref_mass_list(self):
 81        """
 82        Load reference mass list (Bruker format).
 83
 84        Loads in a reference mass list from a .ref file. Some versions of
 85        Bruker's software produce .ref files with a different format where
 86        the header lines (starting with '#' or '##') and delimiters may vary.
 87        The file may be located locally or on S3 and will be handled accordingly.
 88        
 89        Returns
 90        -------
 91        df_ref : Pandas DataFrame
 92            Reference mass list object.
 93        """
 94        # Get a Path-like object from the input path string or S3Path
 95        refmasslist = (Path(self.ref_mass_list_path)
 96                    if isinstance(self.ref_mass_list_path, str)
 97                    else self.ref_mass_list_path)
 98        
 99        # Make sure the file exists
100        if not refmasslist.exists():
101            raise FileExistsError("File does not exist: %s" % refmasslist)
102        
103        # Read all lines from the file (handling S3 vs local differently)
104        if isinstance(refmasslist, S3Path):
105            # For S3, read the file in binary, then decode to string and split into lines.
106            content = refmasslist.open("rb").read()
107            all_lines = content.decode("utf-8").splitlines(keepends=True)
108        else:
109            # For a local file, open in text mode and read lines.
110            with refmasslist.open("r") as f:
111                all_lines = f.readlines()
112        
113        # Identify the index of the first line of the actual data.
114        # We assume header lines start with '#' (or '##') and ignore blank lines.
115        data_start_idx = 0
116        for idx, line in enumerate(all_lines):
117            if line.strip() and not line.lstrip().startswith("#"):
118                data_start_idx = idx
119                break
120        
121        # If there are not enough lines to guess the dialect, throw an error
122        if data_start_idx >= len(all_lines):
123            raise ValueError("The file does not appear to contain any data lines.")
124        
125        # Use a couple of the data lines to let csv.Sniffer detect the delimiter
126        sample_lines = "".join(all_lines[data_start_idx:data_start_idx+2])
127        try:
128            dialect = csv.Sniffer().sniff(sample_lines)
129            delimiter = dialect.delimiter
130        except csv.Error:
131            # If csv.Sniffer fails, default to a common delimiter (e.g., comma)
132            delimiter = ","
133        
134        # Join the lines from the beginning of data (this might include a blank line) 
135        joined_data = "".join(all_lines[data_start_idx:])
136        
137        # Depending on whether the file is S3 or local, wrap the data as needed for pandas
138        if isinstance(refmasslist, S3Path):
139            data_buffer = BytesIO(joined_data.encode("utf-8"))
140        else:
141            data_buffer = StringIO(joined_data)
142        
143        # Read data into a DataFrame.
144        # Adjust columns and names as needed – here we assume at least 2 columns:
145        df_ref = pd.read_csv(data_buffer,
146                            sep=delimiter,
147                            header=None,
148                            usecols=[0, 1],   # Modify if more columns are required.
149                            names=["Formula", "m/z"])
150        
151        df_ref.sort_values(by="m/z", ascending=True, inplace=True)
152        
153        if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
154            print("Reference mass list loaded - {} calibration masses loaded.".format(len(df_ref)))
155        
156        return df_ref
157
158    def gen_ref_mass_list_from_assigned(self, min_conf: float = 0.7):
159        """Generate reference mass list from assigned masses
160
161        This function will generate a ref mass dataframe object from an assigned corems mass spec obj
162        using assigned masses above a certain minimum confidence threshold.
163
164        This function needs to be retested and check it is covered in the unit tests.
165
166        Parameters
167        ----------
168        min_conf : float, optional
169            minimum confidence score. The default is 0.7.
170
171        Returns
172        -------
173        df_ref : Pandas DataFrame
174            reference mass list - based on calculated masses.
175
176        """
177        # TODO this function needs to be retested and check it is covered in the unit tests
178        df = self.mass_spectrum.to_dataframe()
179        df = df[df["Confidence Score"] > min_conf]
180        df_ref = pd.DataFrame(columns=["m/z"])
181        df_ref["m/z"] = df["Calculated m/z"]
182        if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
183            print(
184                "Reference mass list generated - "
185                + str(len(df_ref))
186                + " calibration masses."
187            )
188        return df_ref
189
190    def find_calibration_points(
191        self,
192        df_ref,
193        calib_ppm_error_threshold: tuple[float, float] = (-1, 1),
194        calib_snr_threshold: float = 5,
195        calibration_ref_match_method: str = "merged",
196        calibration_ref_match_tolerance: float = 0.003,
197        calibration_ref_match_std_raw_error_limit: float = 1.5,
198    ):
199        """Function to find calibration points in the mass spectrum
200
201        Based on the reference mass list.
202
203        Parameters
204        ----------
205        df_ref : Pandas DataFrame
206            reference mass list for recalibration.
207        calib_ppm_error_threshold : tuple of floats, optional
208            ppm error for finding calibration masses in the spectrum. The default is -1,1.
209            Note: This is based on the calculation of ppm = ((mz_measure - mz_theoretical)/mz_theoretical)*1e6.
210                Some software does this the other way around and value signs must be inverted for that to work.
211        calib_snr_threshold : float, optional
212            snr threshold for finding calibration masses in the spectrum. The default is 5.
213            If SNR data is unavailable, peaks are filtered by intensity percentile using the formula:
214            percentile = max(5, 100 - calib_snr_threshold)
215        calibration_ref_match_method : str, optional
216            method for matching calibration references. The default is "merged".
217        calibration_ref_match_tolerance : float, optional
218            tolerance for matching calibration references. The default is 0.003.
219        calibration_ref_match_std_raw_error_limit : float, optional
220            standard deviation raw error limit for calibration references. The default is 1.5.
221
222        Returns
223        -------
224        cal_peaks_mz : list of floats
225            masses of measured ions to use in calibration routine
226        cal_refs_mz : list of floats
227            reference mz values of found calibration points.
228
229        """
230
231        # Check if SNR data is available by testing the first peak
232        use_snr = False
233        if len(self.mass_spectrum.mspeaks) > 0:
234            first_peak = self.mass_spectrum.mspeaks[0]
235            if (hasattr(first_peak, 'signal_to_noise') and 
236                first_peak.signal_to_noise is not None and 
237                not np.isnan(first_peak.signal_to_noise) and
238                first_peak.signal_to_noise > 0):
239                use_snr = True
240
241        # This approach is much more efficient and expedient than the original implementation.
242        peaks_mz = []
243        peaks_intensity = []
244        
245        if use_snr:
246            # Use SNR filtering
247            for x in self.mass_spectrum.mspeaks:
248                if x.signal_to_noise > calib_snr_threshold:
249                    if self.mzsegment:
250                        if min(self.mzsegment) <= x.mz_exp <= max(self.mzsegment):
251                            peaks_mz.append(x.mz_exp)
252                    else:
253                        peaks_mz.append(x.mz_exp)
254        else:
255            # Fallback to intensity percentile filtering
256            intensity_percentile = max(5, 100 - calib_snr_threshold)
257            warnings.warn(
258                f"SNR data unavailable for calibration. Using intensity-based filtering instead. "
259                f"SNR threshold of {calib_snr_threshold} corresponds to intensity percentile >= {intensity_percentile}%."
260            )
261            
262            # Collect all peaks and their intensities
263            all_peaks_data = []
264            for x in self.mass_spectrum.mspeaks:
265                if self.mzsegment:
266                    if min(self.mzsegment) <= x.mz_exp <= max(self.mzsegment):
267                        all_peaks_data.append((x.mz_exp, x.abundance))
268                else:
269                    all_peaks_data.append((x.mz_exp, x.abundance))
270            
271            if all_peaks_data:
272                peaks_mz_list, intensities = zip(*all_peaks_data)
273                intensity_threshold = np.percentile(intensities, intensity_percentile)
274                
275                for mz, intensity in all_peaks_data:
276                    if intensity >= intensity_threshold:
277                        peaks_mz.append(mz)
278        
279        peaks_mz = np.asarray(peaks_mz)
280
281        if calibration_ref_match_method == "legacy":
282            # This legacy approach iterates through each reference match and finds the entries within 1 mz and within the user defined PPM error threshold
283            # Then it removes ambiguities - which means the calibration threshold hasto be very tight.
284            cal_peaks_mz = []
285            cal_refs_mz = []
286            for mzref in df_ref["m/z"]:
287                tmp_peaks_mz = peaks_mz[abs(peaks_mz - mzref) < 1]
288                for mzmeas in tmp_peaks_mz:
289                    delta_mass = ((mzmeas - mzref) / mzref) * 1e6
290                    if delta_mass < max(calib_ppm_error_threshold):
291                        if delta_mass > min(calib_ppm_error_threshold):
292                            cal_peaks_mz.append(mzmeas)
293                            cal_refs_mz.append(mzref)
294
295            # To remove entries with duplicated indices (reference masses matching multiple peaks)
296            tmpdf = pd.Series(index=cal_refs_mz, data=cal_peaks_mz, dtype=float)
297            tmpdf = tmpdf[~tmpdf.index.duplicated(keep=False)]
298
299            cal_peaks_mz = list(tmpdf.values)
300            cal_refs_mz = list(tmpdf.index)
301        elif calibration_ref_match_method == "merged":
302            # This is a new approach (August 2024) which uses Pandas 'merged_asof' to find the peaks closest in m/z between
303            # reference and measured masses. This is a quicker way to match, and seems to get more matches.
304            # It may not work as well when the data are far from correc initial mass
305            # e.g. if the correct peak is further from the reference than an incorrect peak.
306            meas_df = pd.DataFrame(columns=["meas_m/z"], data=peaks_mz)
307            tolerance = calibration_ref_match_tolerance
308            merged_df = pd.merge_asof(
309                df_ref,
310                meas_df,
311                left_on="m/z",
312                right_on="meas_m/z",
313                tolerance=tolerance,
314                direction="nearest",
315            )
316            merged_df.dropna(how="any", inplace=True)
317            merged_df["Error_ppm"] = (
318                (merged_df["meas_m/z"] - merged_df["m/z"]) / merged_df["m/z"]
319            ) * 1e6
320            median_raw_error = merged_df["Error_ppm"].median()
321            std_raw_error = merged_df["Error_ppm"].std()
322            if std_raw_error > calibration_ref_match_std_raw_error_limit:
323                std_raw_error = calibration_ref_match_std_raw_error_limit
324            self.mass_spectrum.calibration_raw_error_median = median_raw_error
325            self.mass_spectrum.calibration_raw_error_stdev = std_raw_error
326            merged_df = merged_df[
327                (merged_df["Error_ppm"] > (median_raw_error - 1.5 * std_raw_error))
328                & (merged_df["Error_ppm"] < (median_raw_error + 1.5 * std_raw_error))
329            ]
330            # merged_df= merged_df[(merged_df['Error_ppm']>min(calib_ppm_error_threshold))&(merged_df['Error_ppm']<max(calib_ppm_error_threshold))]
331            cal_peaks_mz = list(merged_df["meas_m/z"])
332            cal_refs_mz = list(merged_df["m/z"])
333        else:
334            raise ValueError(f"{calibration_ref_match_method} not allowed.")
335
336        # it is crucial the mass lists are in same order
337        # corems likes to do masses from high to low.
338        cal_refs_mz.sort(reverse=False)
339        cal_peaks_mz.sort(reverse=False)
340        if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
341            print(
342                str(len(cal_peaks_mz))
343                + " calibration points matched within thresholds."
344            )
345        return cal_peaks_mz, cal_refs_mz
346
347    def robust_calib(
348        self,
349        param: list[float],
350        cal_peaks_mz: list[float],
351        cal_refs_mz: list[float],
352        order: int = 1,
353    ):
354        """Recalibration function
355
356        Computes the rms of m/z errors to minimize when calibrating.
357        This is adapted from from spike.
358
359        Parameters
360        ----------
361        param : list of floats
362            generated by minimize function from scipy optimize.
363        cal_peaks_mz : list of floats
364            masses of measured peaks to use in mass calibration.
365        cal_peaks_mz : list of floats
366            reference mz values of found calibration points.
367        order : int, optional
368            order of the recalibration function. 1 = linear, 2 = quadratic. The default is 1.
369
370        Returns
371        -------
372        rmserror : float
373            root mean square mass error for calibration points.
374
375        """
376        Aterm = param[0]
377        Bterm = param[1]
378        try:
379            Cterm = param[2]
380        except IndexError:
381            pass
382
383        # get the mspeaks from the mass spectrum object which were calibration points
384        # mspeaks = [self.mass_spectrum.mspeaks[x] for x in imzmeas]
385        # get their calibrated mass values
386        # mspeakmzs = [x.mz_cal for x in mspeaks]
387        cal_peaks_mz = np.asarray(cal_peaks_mz)
388
389        # linearz
390        if order == 1:
391            ref_recal_points = (Aterm * cal_peaks_mz) + Bterm
392        # quadratic
393        elif order == 2:
394            ref_recal_points = (Aterm * (cal_peaks_mz)) + (
395                Bterm * np.power((cal_peaks_mz), 2) + Cterm
396            )
397
398        # sort both the calibration points (measured, recalibrated)
399        ref_recal_points.sort()
400        # and sort the calibration points (theoretical, predefined)
401        cal_refs_mz.sort()
402
403        # calculate the ppm error for each calibration point
404        error = ((ref_recal_points - cal_refs_mz) / cal_refs_mz) * 1e6
405        # calculate the root mean square error - this is our target to minimize
406        rmserror = np.sqrt(np.mean(error**2))
407        return rmserror
408
409    def recalibrate_mass_spectrum(
410        self,
411        cal_peaks_mz: list[float],
412        cal_refs_mz: list[float],
413        order: int = 1,
414        diagnostic: bool = False,
415    ):
416        """Main recalibration function which uses a robust linear regression
417
418        This function performs the recalibration of the mass spectrum object.
419        It iteratively applies
420
421        Parameters
422        ----------
423        cal_peaks_mz : list of float
424            masses of measured peaks to use in mass calibration.
425        cal_refs_mz : list of float
426            reference mz values of found calibration points.
427        order : int, optional
428            order of the recalibration function. 1 = linear, 2 = quadratic. The default is 1.
429
430        Returns
431        -------
432        mass_spectrum : CoreMS mass spectrum object
433            Calibrated mass spectrum object
434
435
436        Notes
437        -----
438        This function is adapted, in part, from the SPIKE project [1,2] and is based on the robust linear regression method.
439
440        References
441        ----------
442        1.  Chiron L., Coutouly M-A., Starck J-P., Rolando C., Delsuc M-A.
443            SPIKE a Processing Software dedicated to Fourier Spectroscopies
444            https://arxiv.org/abs/1608.06777 (2016)
445        2.  SPIKE - https://github.com/spike-project/spike
446
447        """
448        # initialise parameters for recalibration
449        # these are the 'Aterm, Bterm, Cterm'
450        # as spectra are already freq->mz calibrated, these terms are very small
451        # may be beneficial to formally separate them from the freq->mz terms
452        if order == 1:
453            Po = [1, 0]
454        elif order == 2:
455            Po = [1, 0, 0]
456
457        if len(cal_peaks_mz) >= 2:
458            if self.mzsegment:  # If only part of the spectrum is to be recalibrated
459                mz_exp_peaks = np.array(
460                    [mspeak.mz_exp for mspeak in self.mass_spectrum]
461                )
462                # Split the array into two parts - one to recailbrate, one to keep unchanged.
463                mz_exp_peaks_tocal = mz_exp_peaks[
464                    (mz_exp_peaks >= min(self.mzsegment))
465                    & (mz_exp_peaks <= max(self.mzsegment))
466                ]
467                mz_exp_peaks_unchanged = mz_exp_peaks[
468                    ~(mz_exp_peaks >= min(self.mzsegment))
469                    | ~(mz_exp_peaks <= max(self.mzsegment))
470                ]
471                # TODO: - segmented calibration needs a way to better track the calibration args/values...
472                if not self.mass_spectrum.is_centroid:
473                    mz_exp_profile = np.array(self.mass_spectrum.mz_exp_profile)
474                    # Split the array into two parts - one to recailbrate, one to keep unchanged.
475                    mz_exp_profile_tocal = mz_exp_profile[
476                        (mz_exp_profile >= min(self.mzsegment))
477                        & (mz_exp_profile <= max(self.mzsegment))
478                    ]
479                    mz_exp_profile_unchanged = mz_exp_profile[
480                        ~(mz_exp_profile >= min(self.mzsegment))
481                        | ~(mz_exp_profile <= max(self.mzsegment))
482                    ]
483            else:  # if just recalibrating the whole spectrum
484                mz_exp_peaks_tocal = np.array(
485                    [mspeak.mz_exp for mspeak in self.mass_spectrum]
486                )
487                if not self.mass_spectrum.is_centroid:
488                    mz_exp_profile_tocal = np.array(self.mass_spectrum.mz_exp_profile)
489
490            minimize_method = self.mass_spectrum.settings.calib_minimize_method
491            res = minimize(
492                self.robust_calib,
493                Po,
494                args=(cal_peaks_mz, cal_refs_mz, order),
495                method=minimize_method,
496            )
497            if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
498                print(
499                    "minimize function completed with RMS error of: {:0.3f} ppm".format(
500                        res["fun"]
501                    )
502                )
503                print(
504                    "minimize function performed {:1d} fn evals and {:1d} iterations".format(
505                        res["nfev"], res["nit"]
506                    )
507                )
508            Pn = res.x
509
510            # mz_exp_ms = np.array([mspeak.mz_exp for mspeak in self.mass_spectrum])
511
512            if order == 1:
513                mz_domain = (Pn[0] * mz_exp_peaks_tocal) + Pn[1]
514                if not self.mass_spectrum.is_centroid:
515                    mz_profile_calc = (Pn[0] * mz_exp_profile_tocal) + Pn[1]
516
517            elif order == 2:
518                mz_domain = (Pn[0] * (mz_exp_peaks_tocal)) + (
519                    Pn[1] * np.power((mz_exp_peaks_tocal), 2) + Pn[2]
520                )
521
522                if not self.mass_spectrum.is_centroid:
523                    mz_profile_calc = (Pn[0] * (mz_exp_profile_tocal)) + (
524                        Pn[1] * np.power((mz_exp_profile_tocal), 2) + Pn[2]
525                    )
526
527            if self.mzsegment:
528                # Recombine the mass domains
529                mz_domain = np.concatenate([mz_domain, mz_exp_peaks_unchanged])
530                mz_domain.sort()
531                if not self.mass_spectrum.is_centroid:
532                    mz_profile_calc = np.concatenate(
533                        [mz_profile_calc, mz_exp_profile_unchanged]
534                    )
535                    mz_profile_calc.sort()
536                # Sort them
537                if (
538                    mz_exp_peaks[0] > mz_exp_peaks[1]
539                ):  # If originally descending mass order
540                    mz_domain = mz_domain[::-1]
541                    if not self.mass_spectrum.is_centroid:
542                        mz_profile_calc = mz_profile_calc[::-1]
543
544            self.mass_spectrum.mz_cal = mz_domain
545            if not self.mass_spectrum.is_centroid:
546                self.mass_spectrum.mz_cal_profile = mz_profile_calc
547
548            self.mass_spectrum.calibration_order = order
549            self.mass_spectrum.calibration_RMS = float(res["fun"])
550            self.mass_spectrum.calibration_points = int(len(cal_refs_mz))
551            self.mass_spectrum.calibration_ref_mzs = cal_refs_mz
552            self.mass_spectrum.calibration_meas_mzs = cal_peaks_mz
553
554            self.mass_spectrum.calibration_segment = self.mzsegment
555
556            if diagnostic:
557                return self.mass_spectrum, res
558            return self.mass_spectrum
559        else:
560            warnings.warn("Too few calibration points - aborting.")
561            return self.mass_spectrum
562
563    def run(self):
564        """Run the calibration routine
565
566        This function runs the calibration routine.
567
568        """
569        calib_snr_threshold = self.mass_spectrum.settings.calib_sn_threshold
570        max_calib_ppm_error = self.mass_spectrum.settings.max_calib_ppm_error
571        min_calib_ppm_error = self.mass_spectrum.settings.min_calib_ppm_error
572        calib_pol_order = self.mass_spectrum.settings.calib_pol_order
573        calibration_ref_match_method = (
574            self.mass_spectrum.settings.calibration_ref_match_method
575        )
576        calibration_ref_match_tolerance = (
577            self.mass_spectrum.settings.calibration_ref_match_tolerance
578        )
579        calibration_ref_match_std_raw_error_limit = (
580            self.mass_spectrum.settings.calibration_ref_match_std_raw_error_limit
581        )
582
583        # load reference mass list
584        df_ref = self.load_ref_mass_list()
585
586        # find calibration points
587        cal_peaks_mz, cal_refs_mz = self.find_calibration_points(
588            df_ref,
589            calib_ppm_error_threshold=(min_calib_ppm_error, max_calib_ppm_error),
590            calib_snr_threshold=calib_snr_threshold,
591            calibration_ref_match_method=calibration_ref_match_method,
592            calibration_ref_match_tolerance=calibration_ref_match_tolerance,
593            calibration_ref_match_std_raw_error_limit=calibration_ref_match_std_raw_error_limit,
594        )
595        if len(cal_peaks_mz) == 2:
596            self.mass_spectrum.settings.calib_pol_order = 1
597            calib_pol_order = 1
598            if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
599                print("Only 2 calibration points found, forcing a linear recalibration")
600        elif len(cal_peaks_mz) < 2:
601            warnings.warn("Too few calibration points found, function will fail")
602        self.recalibrate_mass_spectrum(cal_peaks_mz, cal_refs_mz, order=calib_pol_order)