Method extractXICs

src/topp/MetaProSIP.cpp:1825–1877

{
public:
  static vector<vector<double> > extractXICs(double seed_rt, vector<double> xic_mzs, double mz_toelrance_ppm, double rt_tolerance_s, const PeakMap& peak_map)
  {
    // point on first spectrum in tolerance window
    PeakMap::ConstIterator rt_begin = peak_map.RTBegin(seed_rt - rt_tolerance_s);

    // point on after last spectrum in tolerance window
    PeakMap::ConstIterator rt_end = peak_map.RTBegin(seed_rt + rt_tolerance_s);

    // create set containing all rts of spectra in tolerance window
    set<double> all_rts;
    for (PeakMap::ConstIterator rt_it = rt_begin; rt_it != rt_end; ++rt_it)
    {
      all_rts.insert(rt_it->getRT());
    }

    vector<vector<double> > xics(xic_mzs.size(), vector<double>());

    for (Size i = 0; i < xic_mzs.size(); ++i)
    {
      // create and initialize xic to contain values for all rts
      map<double, double> xic; // rt to summed intensity
      for (set<double>::const_iterator sit = all_rts.begin(); sit != all_rts.end(); ++sit)
      {
        xic[*sit] = 0;
      }

      double mz_da = mz_toelrance_ppm * xic_mzs[i] * 1e-6; // mz tolerance in Dalton
      PeakMap::ConstAreaIterator it = peak_map.areaBeginConst(seed_rt - rt_tolerance_s, seed_rt + rt_tolerance_s, xic_mzs[i] - mz_da, xic_mzs[i] + mz_da);

      for (; it != peak_map.areaEndConst(); ++it)
      {
        double rt = it.getRT();
        if (xic.find(rt) != xic.end())
        {
          xic[rt] += it->getIntensity();
        }
        else
        {
          OPENMS_LOG_WARN << "RT: " << rt << " not contained in rt set." << endl;
        }
      }

      // copy map to vector for easier processing
      vector<double> v;
      for (map<double, double>::const_iterator xic_it = xic.begin(); xic_it != xic.end(); ++xic_it)
      {
        v.push_back(xic_it->second);
      }

      xics[i] = v;
    }
    return xics;
  }

  static vector<double> correlateXICsToMono(const vector<vector<double> >& xics)
  {