.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/01-mass_calibration.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_01-mass_calibration.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_01-mass_calibration.py:


Mass calibration
================

This tutorial demonstrates how to estimate a contrast-to-mass calibration from a
calibration movie:

1. Load a calibration movie.
2. Run detection + PSF fitting to obtain a table of event contrasts.
3. Fit peaks in the contrast histogram.
4. Fit a (linear) contrast-to-mass calibration and plot the result.

.. GENERATED FROM PYTHON SOURCE LINES 13-14

.. code-block:: Python
   :dedent: 1









.. GENERATED FROM PYTHON SOURCE LINES 16-19

Imports
-------
We first import necessary packages and modules.

.. GENERATED FROM PYTHON SOURCE LINES 19-23

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np








.. GENERATED FROM PYTHON SOURCE LINES 24-30

Download example data
---------------------

This tutorial expects a calibration movie. For now, the example downloads a
small calibration file from Zenodo. You can replace this block with a direct
file path once you have a local calibration movie.

.. GENERATED FROM PYTHON SOURCE LINES 30-47

.. code-block:: Python

    from zenodozipdownloader import ZenodoZipDownloader

    import massphoto as mp

    files = ZenodoZipDownloader(
        doi="10.5281/zenodo.15800729",
        download_dir="example_data",
    ).download(inner_pattern="005_mfp1_500x.mp")
    if not files:
        raise RuntimeError(
            "Failed to download example data. Please download it manually from https://zenodo.org/record/15800729 "
            "or set path_calib to your own calibration movie."
        )
    path_calib = files[0]

    print(f"Using calibration movie: {path_calib}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Using calibration movie: /home/runner/work/massphoto/massphoto/docs/source/_examples/example_data/005_mfp1_500x.mp




.. GENERATED FROM PYTHON SOURCE LINES 48-54

Inputs
------

`expected_masses_kDa` should match the ladder present in the calibration movie.
`contrast_slope_guess` is only used to choose sensible histogram bins and
initial peak locations (kDa/contrast).

.. GENERATED FROM PYTHON SOURCE LINES 54-62

.. code-block:: Python


    # Expected masses for the first N peaks (kDa). Adjust to match your ladder.
    expected_masses_kDa = np.arange(1, 8) * 86.0

    # A rough contrast-to-mass slope estimate used to choose sensible binning and
    # initial peak locations (kDa/contrast). If unsure, keep the default.
    contrast_slope_guess = 3e4








.. GENERATED FROM PYTHON SOURCE LINES 63-68

Run detection + fitting
-----------------------

This is the main pipeline step. The output is a DataFrame with fitted event
positions and a `contrasts` column.

.. GENERATED FROM PYTHON SOURCE LINES 68-77

.. code-block:: Python



    movie, _metadata = mp.read_movie(path_calib)
    fitres_calib = mp.fit(movie, movie_correct_drift=True)

    print(f"Detected+fit events used for calibration: n={len(fitres_calib)}")
    fitres_calib[["frame", "x", "y", "contrasts"]].head()






.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    process movie chunks:   0%|          | 0/36 [00:00<?, ?it/s]
    process movie chunks: 100%|██████████| 36/36 [00:00<00:00, 2060.44it/s]

    process movie chunks:   0%|          | 0/31 [00:00<?, ?it/s]
    process movie chunks:   3%|▎         | 1/31 [00:00<00:10,  2.97it/s]
    process movie chunks:  29%|██▉       | 9/31 [00:00<00:01, 15.02it/s]
    process movie chunks:  55%|█████▍    | 17/31 [00:01<00:00, 19.08it/s]
    process movie chunks:  74%|███████▍  | 23/31 [00:01<00:00, 25.04it/s]
    process movie chunks:  87%|████████▋ | 27/31 [00:01<00:00, 22.07it/s]
    process movie chunks: 100%|██████████| 31/31 [00:01<00:00, 22.53it/s]

    process movie chunks:   0%|          | 0/289 [00:00<?, ?it/s]
    process movie chunks:   4%|▍         | 11/289 [00:00<00:11, 24.53it/s]
    process movie chunks:  16%|█▋        | 47/289 [00:00<00:04, 60.50it/s]
    process movie chunks:  19%|█▊        | 54/289 [00:01<00:06, 34.86it/s]
    process movie chunks:  26%|██▋       | 76/289 [00:01<00:04, 44.03it/s]
    process movie chunks:  29%|██▉       | 85/289 [00:01<00:04, 48.96it/s]
    process movie chunks:  32%|███▏      | 92/289 [00:02<00:08, 23.26it/s]
    process movie chunks:  36%|███▌      | 104/289 [00:07<00:26,  7.02it/s]
    process movie chunks:  37%|███▋      | 108/289 [00:07<00:25,  7.15it/s]
    process movie chunks:  51%|█████     | 148/289 [00:08<00:08, 15.77it/s]
    process movie chunks:  61%|██████    | 176/289 [00:09<00:05, 21.65it/s]
    process movie chunks:  62%|██████▏   | 180/289 [00:10<00:06, 16.40it/s]
    process movie chunks:  67%|██████▋   | 195/289 [00:10<00:04, 19.83it/s]
    process movie chunks:  69%|██████▉   | 199/289 [00:10<00:04, 20.69it/s]
    process movie chunks:  70%|███████   | 203/289 [00:11<00:04, 17.43it/s]
    process movie chunks:  76%|███████▌  | 220/289 [00:11<00:03, 18.89it/s]
    process movie chunks:  77%|███████▋  | 223/289 [00:12<00:03, 18.66it/s]
    process movie chunks:  82%|████████▏ | 237/289 [00:12<00:01, 26.54it/s]
    process movie chunks:  83%|████████▎ | 241/289 [00:12<00:02, 20.54it/s]
    process movie chunks:  89%|████████▉ | 258/289 [00:13<00:01, 23.46it/s]
    process movie chunks:  91%|█████████▏| 264/289 [00:13<00:01, 21.95it/s]
    process movie chunks:  94%|█████████▍| 273/289 [00:13<00:00, 25.44it/s]
    process movie chunks:  96%|█████████▌| 277/289 [00:14<00:00, 14.99it/s]
    process movie chunks:  99%|█████████▉| 286/289 [00:14<00:00, 19.02it/s]
    process movie chunks: 100%|██████████| 289/289 [00:16<00:00,  9.40it/s]
    process movie chunks: 100%|██████████| 289/289 [00:16<00:00, 17.78it/s]
    Detected+fit events used for calibration: n=2281


.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>frame</th>
          <th>x</th>
          <th>y</th>
          <th>contrasts</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>20</td>
          <td>94.167944</td>
          <td>32.324291</td>
          <td>0.002724</td>
        </tr>
        <tr>
          <th>1</th>
          <td>38</td>
          <td>123.995974</td>
          <td>20.961641</td>
          <td>0.010899</td>
        </tr>
        <tr>
          <th>2</th>
          <td>40</td>
          <td>93.441404</td>
          <td>21.917409</td>
          <td>0.002605</td>
        </tr>
        <tr>
          <th>3</th>
          <td>42</td>
          <td>141.535097</td>
          <td>40.692372</td>
          <td>0.005699</td>
        </tr>
        <tr>
          <th>4</th>
          <td>42</td>
          <td>143.770046</td>
          <td>52.605493</td>
          <td>0.008198</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 78-83

Fit contrast histogram peaks
----------------------------

We model the histogram as a sum of bin-integrated Gaussians. The fitted peak
centres (`mu`) are used as contrast values for the calibration fit.

.. GENERATED FROM PYTHON SOURCE LINES 83-96

.. code-block:: Python

    bins = np.arange(-20, 660, 4) / contrast_slope_guess
    initial_peaks = expected_masses_kDa / contrast_slope_guess

    mpfit_calib = mp.distfit.fit(
        fitres_calib["contrasts"],
        bins=bins,
        initial_peaks=initial_peaks,
        include_background=False,
    )

    print("\nFitted peak centres (contrast):")
    mpfit_calib.peaks[["mu", "sigma", "A"]]





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Fitted peak centres (contrast):


.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>mu</th>
          <th>sigma</th>
          <th>A</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.002924</td>
          <td>0.000348</td>
          <td>803.369926</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.005812</td>
          <td>0.000433</td>
          <td>518.591835</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.008627</td>
          <td>0.000440</td>
          <td>282.832731</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.011308</td>
          <td>0.000563</td>
          <td>169.690019</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.014191</td>
          <td>0.000526</td>
          <td>93.389237</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.016910</td>
          <td>0.000717</td>
          <td>70.160655</td>
        </tr>
        <tr>
          <th>6</th>
          <td>0.019726</td>
          <td>0.000333</td>
          <td>26.722994</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 97-102

Plot contrast histogram
-----------------------

This is a useful sanity check: you should see distinct peaks and a reasonable
fit overlay.

.. GENERATED FROM PYTHON SOURCE LINES 102-111

.. code-block:: Python

    plt.figure(figsize=(6, 3))
    plt.hist(fitres_calib["contrasts"], bins=bins, edgecolor="none")
    plt.plot(mpfit_calib.model_x, mpfit_calib.model_y, color="black", lw=1, label="fit")
    plt.xlabel("Contrast c")
    plt.ylabel("Counts")
    plt.title("Contrast histogram")
    plt.legend()
    plt.tight_layout()




.. image-sg:: /auto_examples/images/sphx_glr_01-mass_calibration_001.png
   :alt: Contrast histogram
   :srcset: /auto_examples/images/sphx_glr_01-mass_calibration_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 112-117

Fit contrast-to-mass calibration
--------------------------------

We fit a polynomial (default: linear) mapping contrast to mass (kDa). Often the
lower-mass peaks are the cleanest, so we only use the first few for the fit.

.. GENERATED FROM PYTHON SOURCE LINES 117-125

.. code-block:: Python


    n_values_for_fit = min(5, len(expected_masses_kDa), len(mpfit_calib.peaks))
    calibration = mp.mass_calibration.fit(
        mpfit_calib.peaks["mu"][:n_values_for_fit],
        expected_masses_kDa[:n_values_for_fit],
    )
    print("\nFitted calibration (contrast -> mass, kDa):", calibration)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Fitted calibration (contrast -> mass, kDa): -4.98688693 + 30677.65217447·x




.. GENERATED FROM PYTHON SOURCE LINES 126-133

Update from existing calibration (mass-domain peaks)
----------------------------------------------------

If you already have a previous calibration, measured peak positions may be in
mass units already. In that case, pass those values to ``fit``
with ``input_unit='auto'`` (or ``'mass'``) and provide the previous calibration
for inversion to contrast before fitting.

.. GENERATED FROM PYTHON SOURCE LINES 133-147

.. code-block:: Python


    previous_calibration = (-5.2, 3.1e4)
    measured_peak_masses_kDa = mp.mass_calibration.apply(
        mpfit_calib.peaks["mu"][:n_values_for_fit].to_numpy(),
        previous_calibration,
    )
    calibration_updated = mp.mass_calibration.fit(
        measured_peak_masses_kDa,
        expected_masses_kDa[:n_values_for_fit],
        input_unit="auto",
        reference_calibration=previous_calibration,
    )
    print("Updated calibration from mass-domain peak positions:", calibration_updated)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Updated calibration from mass-domain peak positions: -4.98688693 + 30677.65217447·x




.. GENERATED FROM PYTHON SOURCE LINES 148-149

Optionally, update the package wide defaults (`massphoto.defaults.mass_calibration`) with the fitted calibration.

.. GENERATED FROM PYTHON SOURCE LINES 149-154

.. code-block:: Python

    calibration = mp.mass_calibration.fit(
        mpfit_calib.peaks["mu"][:n_values_for_fit],
        expected_masses_kDa[:n_values_for_fit],
        update_defaults=True,
    )







.. GENERATED FROM PYTHON SOURCE LINES 155-160

Plot mass calibration
---------------------

This plot shows the fitted peak positions (with width estimates) and the
resulting calibration curve.

.. GENERATED FROM PYTHON SOURCE LINES 160-181

.. code-block:: Python



    plt.figure(figsize=(6, 3))
    plt.errorbar(
        x=mpfit_calib.peaks["mu"][:n_values_for_fit],
        y=expected_masses_kDa[:n_values_for_fit],
        xerr=mpfit_calib.peaks["sigma"][:n_values_for_fit],
        fmt="|",
        color="black",
        capsize=0,
        linewidth=2,
        label="peaks used for fit",
    )
    x_fit = np.array([bins[0], bins[-1]], dtype=float)
    plt.plot(x_fit, calibration(x_fit), lw=1, label="linear fit")
    plt.xlabel("Contrast c")
    plt.ylabel("Expected mass m (kDa)")
    plt.title("Mass calibration")
    plt.legend()
    plt.tight_layout()




.. image-sg:: /auto_examples/images/sphx_glr_01-mass_calibration_002.png
   :alt: Mass calibration
   :srcset: /auto_examples/images/sphx_glr_01-mass_calibration_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 182-186

Apply calibration
-----------------

Apply the mass calibration to add a ``masses_kDa`` column to the fit results.

.. GENERATED FROM PYTHON SOURCE LINES 186-191

.. code-block:: Python


    fitres_calib = mp.mass_calibration.apply(fitres_calib, calibration)

    fitres_calib[["frame", "x", "y", "contrasts", "masses_kDa"]].head()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>frame</th>
          <th>x</th>
          <th>y</th>
          <th>contrasts</th>
          <th>masses_kDa</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>20</td>
          <td>94.167944</td>
          <td>32.324291</td>
          <td>0.002724</td>
          <td>78.586279</td>
        </tr>
        <tr>
          <th>1</th>
          <td>38</td>
          <td>123.995974</td>
          <td>20.961641</td>
          <td>0.010899</td>
          <td>329.372803</td>
        </tr>
        <tr>
          <th>2</th>
          <td>40</td>
          <td>93.441404</td>
          <td>21.917409</td>
          <td>0.002605</td>
          <td>74.915870</td>
        </tr>
        <tr>
          <th>3</th>
          <td>42</td>
          <td>141.535097</td>
          <td>40.692372</td>
          <td>0.005699</td>
          <td>169.849288</td>
        </tr>
        <tr>
          <th>4</th>
          <td>42</td>
          <td>143.770046</td>
          <td>52.605493</td>
          <td>0.008198</td>
          <td>246.513887</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 192-197

Mass distribution after calibration
-----------------------------------

Finally, inspect the calibrated mass distribution. The expected ladder masses
are shown as dashed vertical lines.

.. GENERATED FROM PYTHON SOURCE LINES 197-218

.. code-block:: Python

    plt.figure(figsize=(6, 3))
    plt.hist(
        fitres_calib["masses_kDa"],
        bins=np.arange(-100, 700, 2.5),
        edgecolor="none",
        label="Single-protein events",
    )
    for mass in expected_masses_kDa:
        plt.axvline(
            mass,
            color="red",
            ls="--",
            lw=1,
            label="Expected masses" if mass == expected_masses_kDa[0] else None,
        )
    plt.xlabel("Mass (kDa)")
    plt.ylabel("Counts")
    plt.title("Mass distribution after calibration")
    plt.legend()
    plt.tight_layout()




.. image-sg:: /auto_examples/images/sphx_glr_01-mass_calibration_003.png
   :alt: Mass distribution after calibration
   :srcset: /auto_examples/images/sphx_glr_01-mass_calibration_003.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 36.692 seconds)


.. _sphx_glr_download_auto_examples_01-mass_calibration.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: 01-mass_calibration.ipynb <01-mass_calibration.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: 01-mass_calibration.py <01-mass_calibration.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: 01-mass_calibration.zip <01-mass_calibration.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_