MSExperiment.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: Marc Sturm, Tom Waschischeck $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/KERNEL/AreaIterator.h>
#include <OpenMS/KERNEL/MSChromatogram.h>
#include <OpenMS/KERNEL/MSSpectrum.h>
#include <OpenMS/METADATA/ExperimentalSettings.h>
#include <OpenMS/DATASTRUCTURES/Matrix.h>
 
#include <vector>
 
 
namespace OpenMS
{
  class Peak1D;
  class ChromatogramPeak;
 
  using MSRun = MSExperiment;
  
  class OPENMS_DLLAPI MSExperiment final : public RangeManagerContainer<RangeRT, RangeMZ, RangeIntensity, RangeMobility>,
    public ExperimentalSettings
  {
 
public:
    typedef Peak1D PeakT;
    typedef ChromatogramPeak ChromatogramPeakT;
 
 
    typedef PeakT PeakType;
    typedef ChromatogramPeakT ChromatogramPeakType;
    typedef PeakType::CoordinateType CoordinateType;
    typedef PeakType::IntensityType IntensityType;
    typedef RangeManager<RangeRT, RangeMZ, RangeIntensity, RangeMobility> RangeManagerType;
    typedef RangeManagerContainer<RangeRT, RangeMZ, RangeIntensity, RangeMobility> RangeManagerContainerType;
    typedef MSSpectrum SpectrumType;
    typedef MSChromatogram ChromatogramType;
    typedef std::vector<SpectrumType> Base;
 
 
    typedef std::vector<SpectrumType>::iterator Iterator;
    typedef std::vector<SpectrumType>::const_iterator ConstIterator;
    typedef Internal::AreaIterator<PeakT, PeakT&, PeakT*, Iterator, SpectrumType::Iterator> AreaIterator;
    typedef Internal::AreaIterator<const PeakT, const PeakT&, const PeakT*, ConstIterator, SpectrumType::ConstIterator> ConstAreaIterator;
 
    // Attention: these refer to the spectra vector only!
    typedef Base::value_type value_type;
    typedef Base::iterator iterator;
    typedef Base::const_iterator const_iterator;
 
    MSExperiment();
 
    MSExperiment(const MSExperiment & source);
 
    MSExperiment(MSExperiment&&) = default;
 
    MSExperiment & operator=(const MSExperiment & source);
 
    MSExperiment& operator=(MSExperiment&&) & = default;
 
    MSExperiment & operator=(const ExperimentalSettings & source);
 
    ~MSExperiment() override;
 
    bool operator==(const MSExperiment & rhs) const;
 
    bool operator!=(const MSExperiment & rhs) const;
    
    inline Size size() const noexcept
    {
      return spectra_.size();
    }
 
    inline void resize(Size n)
    {
      spectra_.resize(n);
    }
 
    inline bool empty() const noexcept
    {
      return spectra_.empty();
    }
    
    inline void reserve(Size n)
    {
      spectra_.reserve(n);
    }
 
    inline SpectrumType& operator[](Size n)
    {
      return spectra_[n];
    }
 
    inline const SpectrumType& operator[](Size n) const
    {
      return spectra_[n];
    }
 
    inline Iterator begin() noexcept
    {
      return spectra_.begin();
    }
 
    inline ConstIterator begin() const noexcept
    {
      return spectra_.cbegin();
    }
 
    inline ConstIterator cbegin() const noexcept
    {
      return spectra_.cbegin();
    }
 
    inline Iterator end()
    {
      return spectra_.end();
    }
 
    inline ConstIterator end() const noexcept
    {
      return spectra_.cend();
    }
    
    inline ConstIterator cend() const noexcept
    {
      return spectra_.cend();
    }
 
    // Aliases / chromatograms
    void reserveSpaceSpectra(Size s);
    void reserveSpaceChromatograms(Size s);
 
 
 
    template <class Container>
    void get2DData(Container& cont) const
    {
      for (typename Base::const_iterator spec = spectra_.begin(); spec != spectra_.end(); ++spec)
      {
        if (spec->getMSLevel() != 1)
        {
          continue;
        }
        typename Container::value_type s; // explicit object here, since instantiation within push_back() fails on VS<12
        for (typename SpectrumType::const_iterator it = spec->begin(); it != spec->end(); ++it)
        {
          cont.push_back(s);
          cont.back().setRT(spec->getRT());
          cont.back().setMZ(it->getMZ());
          cont.back().setIntensity(it->getIntensity());
        }
      }
    }
 
    template <class Container>
    void set2DData(const Container& container)
    {
      set2DData<false, Container>(container);
    }
 
    template <class Container>
    void set2DData(const Container& container, const StringList& store_metadata_names)
    {
      // clean up the container first
      clear(true);
      SpectrumType* spectrum = nullptr;
      typename PeakType::CoordinateType current_rt = -std::numeric_limits<typename PeakType::CoordinateType>::max();
      for (typename Container::const_iterator iter = container.begin(); iter != container.end(); ++iter)
      {
        // check if the retention time has changed
        if (current_rt != iter->getRT() || spectrum == nullptr)
        {
          // append new spectrum
          if (current_rt > iter->getRT())
          {
            throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Input container is not sorted!");
          }
          current_rt =  iter->getRT();
          spectrum = createSpec_(current_rt, store_metadata_names);
        }
 
        // add either data point or mass traces (depending on template argument value)
        ContainerAdd_<typename Container::value_type, false>::addData_(spectrum, &(*iter), store_metadata_names);
      }
    }
 
    template <bool add_mass_traces, class Container>
    void set2DData(const Container& container)
    {
      // clean up the container first
      clear(true);
      SpectrumType* spectrum = nullptr;
      typename PeakType::CoordinateType current_rt = -std::numeric_limits<typename PeakType::CoordinateType>::max();
      for (typename Container::const_iterator iter = container.begin(); iter != container.end(); ++iter)
      {
        // check if the retention time has changed
        if (current_rt != iter->getRT() || spectrum == nullptr)
        {
          // append new spectrum
          if (current_rt > iter->getRT())
          {
            throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Input container is not sorted!");
          }
          current_rt =  iter->getRT();
          spectrum = createSpec_(current_rt);
        }
 
        // add either data point or mass traces (depending on template argument value)
        ContainerAdd_<typename Container::value_type, add_mass_traces>::addData_(spectrum, &(*iter));
      }
    }
 
 
 
 
    AreaIterator areaBegin(CoordinateType min_rt, CoordinateType max_rt, 
      CoordinateType min_mz, CoordinateType max_mz, UInt ms_level = 1);
 
    AreaIterator areaBegin(const RangeManagerType& range, UInt ms_level = 1);
 
    AreaIterator areaEnd();
 
    ConstAreaIterator areaBeginConst(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz, UInt ms_level = 1) const;
 
    ConstAreaIterator areaBeginConst(const RangeManagerType& range, UInt ms_level = 1) const;
 
    ConstAreaIterator areaEndConst() const;
 
    /* @brief Retrieves the peak data in the given mz-rt range and store data spectrum-wise in separate arrays.
     * 
     * For fast pyOpenMS access to peak data in format: [rt, [mz, intensity]]
     * 
     * @param min_rt The minimum retention time.
     * @param max_rt The maximum retention time.
     * @param min_mz The minimum m/z value.
     * @param max_mz The maximum m/z value.
     * @param ms_level The MS level of the spectra to consider.
     * @param rt The vector to store the retention times in.
     * @param mz The vector to store the m/z values in.
     * @param intensity The vector to store the intensities in.
     */
    void get2DPeakDataPerSpectrum(
      CoordinateType min_rt, 
      CoordinateType max_rt, 
      CoordinateType min_mz, 
      CoordinateType max_mz,
      Size ms_level,
      std::vector<float>& rt, 
      std::vector<std::vector<float>>& mz, 
      std::vector<std::vector<float>>& intensity) const
    {
      float t = -1.0;
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz, ms_level); it != areaEndConst(); ++it)
      {
        if (it.getRT() != t) 
        {
          t = (float)it.getRT();
          rt.push_back(t);
          mz.push_back(std::vector<float>());
          intensity.push_back(std::vector<float>());
        }
        mz.back().push_back((float)it->getMZ());
        intensity.back().push_back(it->getIntensity());
      }
    }
 
    /* @brief Retrieves the peak data in the given mz-rt range and store data spectrum-wise in separate arrays.
     * 
     * For fast pyOpenMS access to MS1 peak data in format: [rt, [mz, intensity, ion mobility]]
     * 
     * @param min_rt The minimum retention time.
     * @param max_rt The maximum retention time.
     * @param min_mz The minimum m/z value.
     * @param max_mz The maximum m/z value.
     * @param ms_level The MS level of the spectra to consider.
     * @param rt The vector to store the retention times in.
     * @param mz The vector to store the m/z values in.
     * @param intensity The vector to store the intensities in.
     * @param ion_mobility The vector to store the ion mobility values in.
    */
    void get2DPeakDataIMPerSpectrum(
      CoordinateType min_rt, 
      CoordinateType max_rt, 
      CoordinateType min_mz, 
      CoordinateType max_mz,
      Size ms_level,     
      std::vector<float>& rt, 
      std::vector<std::vector<float>>& mz,
      std::vector<std::vector<float>>& intensity, 
      std::vector<std::vector<float>>& ion_mobility) const
    {
      DriftTimeUnit unit;
      std::vector<float> im;
      float t = -1.0;
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz, ms_level); it != areaEndConst(); ++it)
      {
        if (it.getRT() != t)
        {
          t = (float)it.getRT();
          rt.push_back(t);
          std::tie(unit, im) = it.getSpectrum().maybeGetIMData();
          mz.push_back(std::vector<float>());
          intensity.push_back(std::vector<float>());
          ion_mobility.push_back(std::vector<float>());
        }
 
        if (unit != DriftTimeUnit::NONE)
        {
          const Size peak_index = it.getPeakIndex().peak;
          ion_mobility.back().push_back(im[peak_index]);
        }
        else
        {
          ion_mobility.back().push_back(-1.0);
        }
        mz.back().push_back((float)it->getMZ());
        intensity.back().push_back(it->getIntensity());
      }
    }
 
    /* @brief Retrieves the peak data in the given mz-rt range and store in separate arrays.
     * 
     * For fast pyOpenMS access to MS1 peak data in format: [rt, mz, intensity]
     * 
     * @param min_rt The minimum retention time.
     * @param max_rt The maximum retention time.
     * @param min_mz The minimum m/z value.
     * @param max_mz The maximum m/z value.
     * @param ms_level The MS level of the spectra to consider.
     * @param rt The vector to store the retention times in.
     * @param mz The vector to store the m/z values in.
     * @param intensity The vector to store the intensities in.
    */    
    void get2DPeakData(
      CoordinateType min_rt,
      CoordinateType max_rt,
      CoordinateType min_mz,
      CoordinateType max_mz,
      Size ms_level,
      std::vector<float>& rt,
      std::vector<float>& mz,
      std::vector<float>& intensity) 
      const
    {
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz, ms_level); it != areaEndConst(); ++it)
      {
        rt.push_back((float)it.getRT());
        mz.push_back((float)it->getMZ());
        intensity.push_back(it->getIntensity());
      }
    }
 
 
    /* @brief Retrieves the peak data in the given mz-rt range and store in separate arrays.
     * 
     * For fast pyOpenMS access to MS1 peak data in format: [rt, mz, intensity, ion mobility]
     * 
     * @param min_rt The minimum retention time.
     * @param max_rt The maximum retention time.
     * @param min_mz The minimum m/z value.
     * @param max_mz The maximum m/z value.
     * @param ms_level The MS level of the spectra to consider.
     * @param rt The vector to store the retention times in.
     * @param mz The vector to store the m/z values in.
     * @param intensity The vector to store the intensities in.
    */
    void get2DPeakDataIM(
      CoordinateType min_rt,
      CoordinateType max_rt,
      CoordinateType min_mz,
      CoordinateType max_mz,
      Size ms_level,
      std::vector<float>& rt,
      std::vector<float>& mz,
      std::vector<float>& intensity,
      std::vector<float>& ion_mobility) const
    {
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz, ms_level); it != areaEndConst(); ++it)
      {
        DriftTimeUnit unit = DriftTimeUnit::NONE;
        std::vector<float> im;
        float t = -1.0;
        if (it.getRT() != t)
        {
          t = (float)it.getRT();
          std::tie(unit, im) = it.getSpectrum().maybeGetIMData();
        }
        rt.push_back((float)it.getRT());
        mz.push_back((float)it->getMZ());
        intensity.push_back(it->getIntensity());
        if (unit != DriftTimeUnit::NONE)
        {
          const Size peak_index = it.getPeakIndex().peak;
          ion_mobility.push_back(im[peak_index]);
        }
        else
        {
          ion_mobility.push_back(-1.0);
        }
      }
    }
 
struct SumIntensityReduction {
 
  template <typename Iterator>
  auto operator()(Iterator begin, Iterator end) const {
    // Static assert to verify iterator type has intensity accessor
    using ValueType = typename std::iterator_traits<Iterator>::value_type;
    using IntensityType = decltype(std::declval<ValueType>().getIntensity());
    static_assert(std::is_member_function_pointer_v<decltype(&ValueType::getIntensity)>,
           "Iterator value type must have getIntensity() member function");
 
    IntensityType sum{};
    for (auto it = begin; it != end; ++it) {
      sum += it->getIntensity();
    }
    return sum;
  }
};
 
template<class MzReductionFunctionType>
std::vector<std::vector<MSExperiment::CoordinateType>> aggregate(
    const std::vector<std::pair<RangeMZ, RangeRT>>& mz_rt_ranges,
    unsigned int ms_level,
    MzReductionFunctionType func_mz_reduction) const
{
    // Early exit if there are no ranges
    if (mz_rt_ranges.empty()) 
    {
      // likely an error, but we return an empty vector instead of throwing an exception for now
      return {};
    }
 
    // Create a view of the spectra with given MS level
    std::vector<std::reference_wrapper<const MSSpectrum>> spectra_view;
    spectra_view.reserve(spectra_.size());
    std::copy_if(spectra_.begin(), spectra_.end(), 
                 std::back_inserter(spectra_view),
                 [ms_level](const auto& spec) { 
                     return spec.getMSLevel() == ms_level; 
                 });
 
    // Early exit if there are no spectra with the given MS level
    if (spectra_view.empty()) { // could be valid use or an error -> we return an empty vector
      return {};
    }
 
    // Get the indices of the spectra covered by the RT ranges by considering the MS level
    // If start and stop are the same, the range is empty
    auto getCoveredSpectra = [](
        const std::vector<std::reference_wrapper<const MSSpectrum>>& spectra_view, 
        const std::vector<std::pair<RangeMZ, RangeRT>>& mz_rt_ranges) 
        ->  std::vector<std::pair<size_t, size_t>>
      {
        std::vector<std::pair<size_t, size_t>> res;
        res.reserve(mz_rt_ranges.size());
        
        for (const auto & mz_rt : mz_rt_ranges) 
        {
          // std::cout << "rt range: " << mz_rt.second.getMin() << " - " << mz_rt.second.getMax() << std::endl;
          // std::cout << "specs start:" << spectra_view[0].get().getRT() << " specs end:" << spectra_view[spectra_view.size() - 1].get().getRT() << std::endl;
          auto start_it = std::lower_bound(spectra_view.begin(), spectra_view.end(), mz_rt.second.getMin(), 
            [](const auto& spec, double rt) 
            { return spec.get().getRT() < rt; });
 
          auto stop_it = std::upper_bound(spectra_view.begin(), spectra_view.end(), mz_rt.second.getMax(), 
            [](double rt, const auto& spec) 
            { return rt < spec.get().getRT(); });
 
            res.emplace_back(
                std::distance(spectra_view.begin(), start_it),
                std::distance(spectra_view.begin(), stop_it)
            );
            // std::cout << "start: " << std::distance(spectra_view.begin(), start_it) << " stop: " << std::distance(spectra_view.begin(), stop_it) << std::endl;
        }
        return res;
      };
 
    // For each range, gets (spectrum start index, spectrum stop index). The spectra covered by each RT range.
    const std::vector<std::pair<size_t, size_t>> rt_ranges_idcs = getCoveredSpectra(spectra_view, mz_rt_ranges);
 
    // Initialize result vector
    std::vector<std::vector<MSExperiment::CoordinateType>> result(mz_rt_ranges.size());
 
    // Initialize counts per spectrum index and total mappings
    std::vector<std::vector<size_t>> spec_idx_to_range_idx(spectra_view.size());
 
    // Build spectrum to range index mapping
    for (size_t i = 0; i < rt_ranges_idcs.size(); ++i) 
    {
        const auto& [start, stop] = rt_ranges_idcs[i];
        result[i].resize(stop - start);
        // std::cout << "start: " << start << " stop: " << stop << std::endl;
        for (size_t j = start; j < stop; ++j) 
        {
            spec_idx_to_range_idx[j].push_back(i);
        }
    }
 
   #pragma omp parallel for schedule(dynamic)
   for (Int64 i = 0; i < (Int64)spec_idx_to_range_idx.size(); ++i) // OpenMP on windows still requires signed loop variable
   {
      if (spec_idx_to_range_idx[i].empty()) continue; // no ranges for this spectrum? skip it
 
      const auto& spec = spectra_view[i].get();
      auto spec_begin = spec.cbegin();
      auto spec_end = spec.cend();
 
      for (size_t range_idx : spec_idx_to_range_idx[i]) 
      {
        const auto& mz_range = mz_rt_ranges[range_idx].first;
        
        // Find data points within MZ range
        auto start_it = spec.PosBegin(spec_begin, mz_range.getMinMZ(), spec_end);
        auto end_it = start_it;
        
        while (end_it != spec_end && end_it->getPosition() <= mz_range.getMaxMZ()) 
        {
          ++end_it;
        }
 
        // std::cout << "calculating reduction on range: " << range_idx << " for spectrum: " << i << " and peaks " << std::distance(spec.begin(), start_it) << " - " << std::distance(spec.begin(), end_it) << std::endl;
 
        // Calculate result using provided reduction function
        result[range_idx][i - rt_ranges_idcs[range_idx].first] = 
            func_mz_reduction(start_it, end_it);
      }
    }
    return result;
  }
 
// Overload without func_mz_reduction parameter (default to SumIntensityReduction). Needed because of template deduction issues
std::vector<std::vector<MSExperiment::CoordinateType>> aggregate(
    const std::vector<std::pair<RangeMZ, RangeRT>>& mz_rt_ranges,
    unsigned int ms_level) const
{
  return aggregate(mz_rt_ranges, ms_level, SumIntensityReduction());
}
 
template<class MzReductionFunctionType>
std::vector<MSChromatogram> extractXICs(
    const std::vector<std::pair<RangeMZ, RangeRT>>& mz_rt_ranges,
    unsigned int ms_level,
    MzReductionFunctionType func_mz_reduction) const
{
    // Early exit if there are no ranges
    if (mz_rt_ranges.empty()) 
    {
      // likely an error, but we return an empty vector instead of throwing an exception for now
      return {};
    }
 
    // Create a view of the spectra with given MS level
    std::vector<std::reference_wrapper<const MSSpectrum>> spectra_view;
    spectra_view.reserve(spectra_.size());
    std::copy_if(spectra_.begin(), spectra_.end(), 
                 std::back_inserter(spectra_view),
                 [ms_level](const auto& spec) { 
                     return spec.getMSLevel() == ms_level; 
                 });
 
    // Early exit if there are no spectra with the given MS level
    if (spectra_view.empty()) { // could be valid use or an error -> we return an empty vector
      return {};
    }
 
    // Get the indices of the spectra covered by the RT ranges by considering the MS level
    // If start and stop are the same, the range is empty
    auto getCoveredSpectra = [](
        const std::vector<std::reference_wrapper<const MSSpectrum>>& spectra_view, 
        const std::vector<std::pair<RangeMZ, RangeRT>>& mz_rt_ranges) 
        ->  std::vector<std::pair<size_t, size_t>>
      {
        std::vector<std::pair<size_t, size_t>> res;
        res.reserve(mz_rt_ranges.size());
 
        for (const auto & mz_rt : mz_rt_ranges) 
        {
          auto start_it = std::lower_bound(spectra_view.begin(), spectra_view.end(), mz_rt.second.getMin(), 
            [](const auto& spec, double rt) 
            { return spec.get().getRT() < rt; });
 
          auto stop_it = std::upper_bound(spectra_view.begin(), spectra_view.end(), mz_rt.second.getMax(), 
            [](double rt, const auto& spec) 
            { return rt < spec.get().getRT(); });
 
            res.emplace_back(
                std::distance(spectra_view.begin(), start_it),
                std::distance(spectra_view.begin(), stop_it)
            );
        }
        return res;
      };
 
    // For each range, gets (spectrum start index, spectrum stop index). The spectra covered by each RT range.
    const std::vector<std::pair<size_t, size_t>> rt_ranges_idcs = getCoveredSpectra(spectra_view, mz_rt_ranges);
 
    // Initialize result vector
    std::vector<MSChromatogram> result(mz_rt_ranges.size());
 
    // Initialize counts per spectrum index and total mappings
    std::vector<std::vector<size_t>> spec_idx_to_range_idx(spectra_view.size());
 
    // Build spectrum to range index mapping
    for (size_t i = 0; i < rt_ranges_idcs.size(); ++i) 
    {        
        const auto& [start, stop] = rt_ranges_idcs[i];
        result[i].resize(stop - start);
        result[i].getProduct().setMZ(
          (mz_rt_ranges[i].first.getMinMZ() + mz_rt_ranges[i].first.getMaxMZ()) / 2.0);
        for (size_t j = start; j < stop; ++j) 
        {
          spec_idx_to_range_idx[j].push_back(i);
        }
    }
 
   #pragma omp parallel for schedule(dynamic)
   for (Int64 i = 0; i < (Int64)spec_idx_to_range_idx.size(); ++i) // OpenMP on windows still requires signed loop variable
   {
      if (spec_idx_to_range_idx[i].empty()) continue; // no ranges for this spectrum? skip 
      
      const auto& spec = spectra_view[i].get();
      const double rt = spec.getRT();
      auto spec_begin = spec.cbegin();
      auto spec_end = spec.cend();
 
      for (size_t range_idx : spec_idx_to_range_idx[i]) 
      {
        const auto& mz_range = mz_rt_ranges[range_idx].first;
        
        // Find data points within MZ range
        auto start_it = spec.PosBegin(spec_begin, mz_range.getMinMZ(), spec_end);
        auto end_it = start_it;
        
        while (end_it != spec_end && end_it->getPosition() <= mz_range.getMaxMZ()) 
        {
          ++end_it;
        }
 
        // Calculate result using provided reduction function
        result[range_idx][i - rt_ranges_idcs[range_idx].first] = 
            ChromatogramPeak(rt, func_mz_reduction(start_it, end_it));
      }
    }
 
    for (auto& r : result) r.updateRanges(); // TODO: prob.. faster to look at first and last peaks as range is sorted
 
    return result;
  }
 
// Overload without func_mz_reduction parameter (needed because of template deduction issue)
std::vector<MSChromatogram> extractXICs(
    const std::vector<std::pair<RangeMZ, RangeRT>>& mz_rt_ranges,
    unsigned int ms_level) const
{
    return extractXICs(mz_rt_ranges, ms_level, SumIntensityReduction());
}
 
  std::vector<std::vector<MSExperiment::CoordinateType>> aggregateFromMatrix(
      const Matrix<double>& ranges,
      unsigned int ms_level,
      const std::string& mz_agg) const
  {
      // Check matrix dimensions
      if (ranges.cols() != 4)
      {
          throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
              "Range matrix must have 4 columns [mz_min, mz_max, rt_min, rt_max]");
      }
 
      // Convert matrix rows to vector of pairs
      std::vector<std::pair<RangeMZ, RangeRT>> mz_rt_ranges;
      mz_rt_ranges.reserve((Size)ranges.rows());
      
      for (Size i = 0; i < (Size)ranges.rows(); ++i)
      {
          mz_rt_ranges.emplace_back(
              RangeMZ(ranges(i, 0), ranges(i, 1)), // min max mz
              RangeRT(ranges(i, 2), ranges(i, 3))  // min max rt
          );
          // std::cout << "mz: " << ranges(i, 0) << " - " << ranges(i, 1) << " rt: " << ranges(i, 2) << " - " << ranges(i, 3) << std::endl;
      }
 
      // Call appropriate aggregation function based on mz_agg parameter
      if (mz_agg == "sum")
      {
          return aggregate(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)
              {
                  return std::accumulate(begin_it, end_it, 0.0,
                      [](double a, const Peak1D& b) { return a + b.getIntensity(); });
              });
      }
      else if (mz_agg == "max")
      {
          return aggregate(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)->double
              {
                  if (begin_it == end_it) return 0.0;
                  return std::max_element(begin_it, end_it,
                      [](const Peak1D& a, const Peak1D& b) { return a.getIntensity() < b.getIntensity(); }
                  )->getIntensity();
              });
      }
      else if (mz_agg == "min")
      {
          return aggregate(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)->double
              {
                  if (begin_it == end_it) return 0.0;
                  return std::min_element(begin_it, end_it,
                      [](const Peak1D& a, const Peak1D& b) { return a.getIntensity() < b.getIntensity(); }
                  )->getIntensity();
              });
      }
      else if (mz_agg == "mean")
      {
          return aggregate(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)
              {
                  if (begin_it == end_it) return 0.0;
                  double sum = std::accumulate(begin_it, end_it, 0.0,
                      [](double a, const Peak1D& b) { return a + b.getIntensity(); });
                  return sum / static_cast<double>(std::distance(begin_it, end_it));
              });
      }
      else
      {
          throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
              "Invalid aggregation function", mz_agg);
      }
  }
 
  std::vector<MSChromatogram> extractXICsFromMatrix(
      const Matrix<double>& ranges,
      unsigned int ms_level,
      const std::string& mz_agg) const
  {
      // Check matrix dimensions
      if (ranges.cols() != 4)
      {
          throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
              "Range matrix must have 4 columns [mz_min, mz_max, rt_min, rt_max]");
      }
 
      // Convert matrix rows to vector of pairs
      std::vector<std::pair<RangeMZ, RangeRT>> mz_rt_ranges;
      mz_rt_ranges.reserve((Size)ranges.rows());
      
      for (Size i = 0; i < (Size)ranges.rows(); ++i)
      {
          mz_rt_ranges.emplace_back(
              RangeMZ(ranges(i, 0), ranges(i, 1)),
              RangeRT(ranges(i, 2), ranges(i, 3))
          );
      }
 
      // Call appropriate extractXICs function based on mz_agg parameter
      if (mz_agg == "sum")
      {
          return extractXICs(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)
              {
                  return std::accumulate(begin_it, end_it, 0.0,
                      [](double a, const Peak1D& b) { return a + b.getIntensity(); });
              });
      }
      else if (mz_agg == "max")
      {
          return extractXICs(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)->double
              {
                  if (begin_it == end_it) return 0.0;
                  return std::max_element(begin_it, end_it,
                      [](const Peak1D& a, const Peak1D& b) { return a.getIntensity() < b.getIntensity(); }
                  )->getIntensity();
              });
      }
      else if (mz_agg == "min")
      {
          return extractXICs(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)->double
              {
                  if (begin_it == end_it) return 0.0;
                  return std::min_element(begin_it, end_it,
                      [](const Peak1D& a, const Peak1D& b) { return a.getIntensity() < b.getIntensity(); }
                  )->getIntensity();
              });
      }
      else if (mz_agg == "mean")
      {
          return extractXICs(mz_rt_ranges, ms_level,
              [](auto begin_it, auto end_it)
              {
                  if (begin_it == end_it) return 0.0;
                  double sum = std::accumulate(begin_it, end_it, 0.0,
                      [](double a, const Peak1D& b) { return a + b.getIntensity(); });
                  return sum / static_cast<double>(std::distance(begin_it, end_it));
              });
      }
      else
      {
          throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
              "Invalid aggregation function", mz_agg);
      }
  }
 
    ConstIterator RTBegin(CoordinateType rt) const;
 
    ConstIterator RTEnd(CoordinateType rt) const;
 
    Iterator RTBegin(CoordinateType rt);
 
    Iterator RTEnd(CoordinateType rt);
 
 
    ConstIterator IMBegin(CoordinateType im) const;
 
    ConstIterator IMEnd(CoordinateType im) const;
 
    // Docu in base class
    void updateRanges() override;
 
    void updateRanges(Int ms_level);
 
    UInt64 getSize() const;
 
    std::vector<UInt> getMSLevels() const;
 
 
    UInt64 getSqlRunID() const;
 
    void setSqlRunID(UInt64 id);
 
 
    void sortSpectra(bool sort_mz = true);
 
    void sortChromatograms(bool sort_rt = true);
 
    bool isSorted(bool check_mz = true) const;
 
 
    void reset();
 
    bool clearMetaDataArrays();
 
    const ExperimentalSettings& getExperimentalSettings() const;
 
    ExperimentalSettings& getExperimentalSettings();
 
    void getPrimaryMSRunPath(StringList& toFill) const;
 
    ConstIterator getPrecursorSpectrum(ConstIterator iterator) const;
 
    int getPrecursorSpectrum(int zero_based_index) const;
 
    ConstIterator getFirstProductSpectrum(ConstIterator iterator) const;
 
    int getFirstProductSpectrum(int zero_based_index) const;
 
    void swap(MSExperiment& from);
 
    void setSpectra(const std::vector<MSSpectrum>& spectra);
    void setSpectra(std::vector<MSSpectrum>&& spectra);
 
    void addSpectrum(const MSSpectrum& spectrum);
    void addSpectrum(MSSpectrum&& spectrum);
 
    const std::vector<MSSpectrum>& getSpectra() const;
 
    std::vector<MSSpectrum>& getSpectra();
 
    ConstIterator getClosestSpectrumInRT(const double RT) const;
    Iterator getClosestSpectrumInRT(const double RT);
 
    ConstIterator getClosestSpectrumInRT(const double RT, UInt ms_level) const;
    Iterator getClosestSpectrumInRT(const double RT, UInt ms_level);
 
    void setChromatograms(const std::vector<MSChromatogram>& chromatograms);
    void setChromatograms(std::vector<MSChromatogram>&& chromatograms);
 
    void addChromatogram(const MSChromatogram& chromatogram);
    void addChromatogram(MSChromatogram&& chrom);
 
    const std::vector<MSChromatogram>& getChromatograms() const;
 
    std::vector<MSChromatogram>& getChromatograms();
 
 
    MSChromatogram& getChromatogram(Size id);
 
    MSSpectrum& getSpectrum(Size id);
 
    Size getNrSpectra() const;
 
    Size getNrChromatograms() const;
 
    const MSChromatogram calculateTIC(float rt_bin_size = 0, UInt ms_level = 1) const;
 
    void clear(bool clear_meta_data);
 
    bool containsScanOfLevel(size_t ms_level) const;
 
    bool hasZeroIntensities(size_t ms_level) const;
 
    bool hasPeptideIdentifications() const;
 
    bool isIMFrame() const;
 
  protected:
    std::vector<MSChromatogram > chromatograms_;
    std::vector<SpectrumType> spectra_;
 
private:
 
    template<typename ContainerValueType, bool addMassTraces>
    struct ContainerAdd_
    {
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item);
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item, const StringList& store_metadata_names);
    };
 
    template<typename ContainerValueType>
    struct ContainerAdd_<ContainerValueType, false>
    {
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item)
      {
        // create temporary peak and insert it into spectrum
        spectrum->insert(spectrum->end(), PeakType());
        spectrum->back().setIntensity(item->getIntensity());
        spectrum->back().setPosition(item->getMZ());
      }
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item, const StringList& store_metadata_names)
      {
        addData_(spectrum, item);
        for (StringList::const_iterator itm = store_metadata_names.begin(); itm != store_metadata_names.end(); ++itm)
        {
          float val = std::numeric_limits<float>::quiet_NaN();
          if (item->metaValueExists(*itm)) val = item->getMetaValue(*itm);
          spectrum->getFloatDataArrays()[itm - store_metadata_names.begin()].push_back(val);
        }
      }
    };
 
    template<typename ContainerValueType>
    struct ContainerAdd_<ContainerValueType, true>
    {
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item)
      {
        if (item->metaValueExists("num_of_masstraces"))
        {
          Size mts = item->getMetaValue("num_of_masstraces");
          int charge = (item->getCharge()==0 ? 1 : item->getCharge()); // set to 1 if charge is 0, otherwise div/0 below
          for (Size i = 0; i < mts; ++i)
          {
            String meta_name = String("masstrace_intensity_") + i;
            if (!item->metaValueExists(meta_name))
            {
              throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, String("Meta value '") + meta_name + "' expected but not found in container.");
            }
            ContainerValueType p;
            p.setIntensity(item->getMetaValue(meta_name));
            p.setPosition(item->getMZ() + Constants::C13C12_MASSDIFF_U / charge * i);
            ContainerAdd_<ContainerValueType, false>::addData_(spectrum, &p);
          }
        }
        else ContainerAdd_<ContainerValueType, false>::addData_(spectrum, item);
      }
    };
 
    /*
      @brief Append a spectrum to current MSExperiment
 
      @param rt RT of new spectrum
      @return Pointer to newly created spectrum
    */
    SpectrumType* createSpec_(PeakType::CoordinateType rt);
 
    /*
      @brief Append a spectrum including floatdata arrays to current MSExperiment
 
      @param rt RT of new spectrum
      @param metadata_names Names of floatdata arrays attached to this spectrum
      @return Pointer to newly created spectrum
    */
    SpectrumType* createSpec_(PeakType::CoordinateType rt, const StringList& metadata_names);
 
  };
 
  OPENMS_DLLAPI std::ostream& operator<<(std::ostream& os, const MSExperiment& exp);
 
} // namespace OpenMS
 
#include <OpenMS/KERNEL/StandardTypes.h>