skyline: Stopped inferring a Waters lockmass function from waters_connect mzML

pwiz/analysis/spectrum_processing/SpectrumList_Filter.cpp

@@ -318,12 +318,25 @@ PWIZ_API_DECL SpectrumList_FilterPredicate_MSLevelSet::SpectrumList_FilterPredic
 
 PWIZ_API_DECL boost::logic::tribool SpectrumList_FilterPredicate_MSLevelSet::accept(const msdata::Spectrum& spectrum) const
 {
-    CVParam param = spectrum.cvParamChild(MS_spectrum_type);
-    if (param.cvid == CVID_Unknown) return boost::logic::indeterminate;
-    if (!cvIsA(param.cvid, MS_mass_spectrum))
-        return msLevelSet_.contains(0); // non-MS spectra are considered ms level 0
-    param = spectrum.cvParam(MS_ms_level);
-    if (param.cvid == CVID_Unknown) return boost::logic::indeterminate;
+    if (!spectrum.hasCVParamChild(MS_spectrum_type))
+        return boost::logic::indeterminate;
+
+    // A declared ms level decides this, whatever the spectrum type says. Not every child of
+    // "spectrum type" is under "mass spectrum" - "calibration spectrum" is not, and the UIMF reader
+    // writes it as the sole type on frames that carry a perfectly good ms level - so keying off the
+    // type first discarded that level and called such a spectrum ms level 0.
+    //
+    // Ask by child rather than for the first one: a spectrum may carry more than one child of
+    // "spectrum type", and no rule fixes the order a writer emits them in.
+    CVParam param = spectrum.cvParam(MS_ms_level);
+    if (param.cvid == CVID_Unknown)
+    {
+        // Nothing declared, so the type is all there is to go on. A spectrum that is not a mass
+        // spectrum at all - an emission spectrum, say - is considered ms level 0.
+        if (!spectrum.hasCVParamChild(MS_mass_spectrum))
+            return msLevelSet_.contains(0);
+        return boost::logic::indeterminate;
+    }
     int msLevel = param.valueAs<int>();
     bool result = msLevelSet_.contains(msLevel);
     return result;

pwiz/analysis/spectrum_processing/SpectrumList_FilterTest.cpp

@@ -518,6 +518,64 @@ void testMSLevelSet(SpectrumListPtr sl)
     }
 }
 
+
+// "calibration spectrum" is a child of "spectrum type" but not of "mass spectrum", and the UIMF
+// reader writes it as the sole type on a frame that still carries an ms level. Such a spectrum must
+// be filtered on the level it declares, not treated as ms level 0 the way a type-only spectrum is.
+// Uses its own list, since the shared one is pinned by exact sizes and ids throughout this file.
+void testMSLevelSetCalibrationSpectrum()
+{
+    if (os_) *os_ << "testMSLevelSetCalibrationSpectrum:\n";
+
+    SpectrumListSimplePtr sl(new SpectrumListSimple);
+
+    SpectrumPtr ms1(new Spectrum);
+    ms1->index = 0;
+    ms1->id = "scan=1";
+    ms1->set(MS_ms_level, 1);
+    ms1->set(MS_MS1_spectrum);
+    sl->spectra.push_back(ms1);
+
+    SpectrumPtr calibration(new Spectrum);
+    calibration->index = 1;
+    calibration->id = "scan=2";
+    calibration->set(MS_ms_level, 1);
+    calibration->set(MS_calibration_spectrum);
+    sl->spectra.push_back(calibration);
+
+    // A spectrum type with no ms level at all is what the level 0 rule is for
+    SpectrumPtr emission(new Spectrum);
+    emission->index = 2;
+    emission->id = "scan=3";
+    emission->set(MS_emission_spectrum);
+    sl->spectra.push_back(emission);
+
+    {
+        SpectrumList_Filter filter(sl, SpectrumList_FilterPredicate_MSLevelSet(IntegerSet(1)));
+        if (os_)
+        {
+            printSpectrumList(filter, *os_);
+            *os_ << endl;
+        }
+        // The calibration spectrum declares ms level 1, so asking for level 1 must keep it
+        unit_assert_operator_equal(2, filter.size());
+        unit_assert(filter.spectrumIdentity(0).id == "scan=1");
+        unit_assert(filter.spectrumIdentity(1).id == "scan=2");
+    }
+
+    {
+        SpectrumList_Filter filter(sl, SpectrumList_FilterPredicate_MSLevelSet(IntegerSet(0)));
+        if (os_)
+        {
+            printSpectrumList(filter, *os_);
+            *os_ << endl;
+        }
+        // And asking for level 0 must not, leaving only the spectrum that declares no level
+        unit_assert_operator_equal(1, filter.size());
+        unit_assert(filter.spectrumIdentity(0).id == "scan=3");
+    }
+}
+
 void testMS2Activation(SpectrumListPtr sl)
 {
     if (os_) *os_ << "testMS2Activation:\n";
@@ -840,6 +898,7 @@ void test()
     testScanEventSet(sl);
     testScanTimeRange(sl);
     testMSLevelSet(sl);
+    testMSLevelSetCalibrationSpectrum();
     testMS2Activation(sl);
     testMassAnalyzerFilter(sl);
     testMZPresentFilter(sl);

pwiz/data/msdata/SpectrumListWrapper.hpp

@@ -55,6 +55,12 @@ class PWIZ_API_DECL SpectrumListWrapper : public SpectrumListBase
     virtual IndexList findNameValue(const std::string& name, const std::string& value) const {return size() == inner_->size() ? inner_->findNameValue(name, value) : SpectrumList::findNameValue(name, value);}
     virtual IndexList findSpotID(const std::string& spotID) const {return size() == inner_->size() ? inner_->findSpotID(spotID) : SpectrumList::findSpotID(spotID);}
 
+    // A wrapper only ever filters or transforms what the inner list presents, so if the inner list is
+    // omitting calibration spectra then so is the wrapper. Without forwarding this, the base
+    // implementation answers false and callers cannot tell a wrapped vendor list from one with no
+    // calibration spectra at all.
+    virtual bool calibrationSpectraAreOmitted() const {return inner_->calibrationSpectraAreOmitted();}
+
     // no default implementation, because otherwise subclasses could override the DetailLevel overload and the getBinaryData overload would be inconsistent
     virtual SpectrumPtr spectrum(size_t index, bool getBinaryData = false) const = 0;
 

pwiz/data/vendor_readers/Waters/ChromatogramList_Waters.cpp

@@ -124,8 +124,19 @@ PWIZ_API_DECL ChromatogramPtr ChromatogramList_Waters::chromatogram(size_t index
 
             multimap<double, pair<int, double>> fullFileTIC;
 
+            // When ignoreCalibrationScans is set the lockmass function is kept out of the spectrum list,
+            // so it must not be summed into the global TIC either - otherwise the TIC carries points that
+            // no spectrum in the file accounts for.
+            int lockmassFunction = -1;
+            bool excludeLockmassFunction = config_.ignoreCalibrationScans &&
+                                           rawdata_->Info.TryGetLockMassFunction(lockmassFunction) &&
+                                           lockmassFunction >= 0;
+
             for(int function : rawdata_->FunctionIndexList())
             {
+                if (excludeLockmassFunction && function == lockmassFunction)
+                    continue;
+
                 if (config_.globalChromatogramsAreMs1Only)
                 {
                     int msLevel;

pwiz/data/vendor_readers/Waters/Reader_Waters_Test.cpp

@@ -130,6 +130,15 @@ int main(int argc, char* argv[])
             result += pwiz::util::testReader(reader, testArgs, testAcceptOnly, requireUnicodeSupport, pwiz::util::IsNamedRawFile({ "HDDDA_Short_noLM.raw" }), newConfig);
         }
 
+        // test ignoreCalibrationScans on a file that actually has a lockspray function. HDDDA_Short_noLM
+        // above has none, so it cannot show that the calibration spectra are dropped, nor that the
+        // global TIC then stops summing them - MSe_Short has function 3.
+        {
+            auto newConfig = config;
+            newConfig.ignoreCalibrationScans = true;
+            result += pwiz::util::testReader(reader, testArgs, testAcceptOnly, requireUnicodeSupport, pwiz::util::IsNamedRawFile({ "MSe_Short.raw" }), newConfig);
+        }
+
         // test ddaProcessing
         {
             auto newConfig = config;

pwiz/data/vendor_readers/Waters/SpectrumList_Waters.cpp

@@ -42,11 +42,21 @@ namespace detail {
 using namespace Waters;
 
 SpectrumList_Waters::SpectrumList_Waters(MSData& msd, RawDataPtr rawdata, const Reader::Config& config)
-    : msd_(msd), rawdata_(rawdata), config_(config), lockmassFunction_(LOCKMASS_FUNCTION_UNINIT)
+    : msd_(msd), rawdata_(rawdata), config_(config), lockmassFunction_(LOCKMASS_FUNCTION_UNKNOWN)
 {
     useDDAProcessor_ = config_.ddaProcessing;
     rawdata_->EnableProcessing(useDDAProcessor_);
 
+    // Resolve the lockmass function once, here, while still single threaded. Doing it lazily would let
+    // spectrum() (which holds readMutex) and calibrationSpectraAreOmitted() (which does not) race on
+    // the member and enter the MassLynx SDK unserialized. Note the out-parameter goes to a local: the
+    // SDK wrapper leaves it untouched when the call fails, so handing it the member directly would
+    // cache whatever happened to be there.
+    int lockmassFunction = LOCKMASS_FUNCTION_UNKNOWN;
+    if (!rawdata_->Info.TryGetLockMassFunction(lockmassFunction) || lockmassFunction < 0)
+        lockmassFunction = LOCKMASS_FUNCTION_UNKNOWN;
+    lockmassFunction_ = lockmassFunction;
+
     if (useDDAProcessor_)
     {       
         createDDAIndex();
@@ -107,16 +117,14 @@ PWIZ_API_DECL SpectrumPtr SpectrumList_Waters::spectrum(size_t index, bool getBi
     return spectrum(index, getBinaryData ? DetailLevel_FullData : DetailLevel_FullMetadata, lockmassMzPosScans, lockmassMzNegScans, lockmassTolerance, msLevelsToCentroid);
 }
 
+PWIZ_API_DECL int SpectrumList_Waters::lockMassFunction() const
+{
+    return lockmassFunction_; // resolved in the constructor, so this is safe to call from any thread
+}
+
 PWIZ_API_DECL bool SpectrumList_Waters::isLockMassFunction(int function) const
 {
-    if (lockmassFunction_ == LOCKMASS_FUNCTION_UNINIT)
-    {
-        if (!rawdata_->Info.TryGetLockMassFunction(lockmassFunction_))
-        {
-            lockmassFunction_ = LOCKMASS_FUNCTION_UNKNOWN;
-        }
-    }
-    return function == lockmassFunction_;
+    return function == lockMassFunction();
 }
 
 PWIZ_API_DECL SpectrumPtr SpectrumList_Waters::spectrum(size_t index, DetailLevel detailLevel, double lockmassMzPosScans, double lockmassMzNegScans, double lockmassTolerance, const pwiz::util::IntegerSet& msLevelsToCentroid) const
@@ -629,7 +637,7 @@ PWIZ_API_DECL void SpectrumList_Waters::calculatePeakMetadata(SpectrumPtr& spect
 
 PWIZ_API_DECL bool SpectrumList_Waters::calibrationSpectraAreOmitted() const
 {
-    return config_.ignoreCalibrationScans && lockmassFunction_ >= 0;
+    return config_.ignoreCalibrationScans && lockMassFunction() >= 0;
 }
 
 PWIZ_API_DECL void SpectrumList_Waters::createIndex()
@@ -834,6 +842,7 @@ bool SpectrumList_Waters::hasCombinedIonMobility() const {return false;}
 bool SpectrumList_Waters::canConvertIonMobilityAndCCS() const {return false;}
 double SpectrumList_Waters::ionMobilityToCCS(double ionMobility, double mz, int charge) const {return 0;}
 double SpectrumList_Waters::ccsToIonMobility(double ccs, double mz, int charge) const {return 0;}
+int SpectrumList_Waters::lockMassFunction() const {return LOCKMASS_FUNCTION_UNKNOWN;}
 bool SpectrumList_Waters::isLockMassFunction(int function) const {return false;}
 bool SpectrumList_Waters::calibrationSpectraAreOmitted() const {return false;}
 SpectrumPtr SpectrumList_Waters::spectrum(size_t index, bool getBinaryData) const {return SpectrumPtr();}

pwiz/data/vendor_readers/Waters/SpectrumList_Waters.hpp

@@ -120,9 +120,15 @@ class PWIZ_API_DECL SpectrumList_Waters : public SpectrumListIonMobilityBase
 
 #endif // PWIZ_READER_WATERS
 
-    mutable int lockmassFunction_; // 0-based. Special values: -1=uninitialized -2=unknown
+    int lockmassFunction_; // 0-based, or LOCKMASS_FUNCTION_UNKNOWN. Resolved by the constructor.
 #define LOCKMASS_FUNCTION_UNKNOWN -2
-#define LOCKMASS_FUNCTION_UNINIT -1
+
+    /// The 0-based lockmass function number, or LOCKMASS_FUNCTION_UNKNOWN if the source has none.
+    /// Resolved once by the constructor, so this is a plain read and is safe from any thread. Do not
+    /// make it lazy again: spectrum() calls it holding readMutex while calibrationSpectraAreOmitted()
+    /// calls it without, so a lazy write would race and would also let two threads into the MassLynx
+    /// SDK at once.
+    int lockMassFunction() const;
 };
 
 } // detail

pwiz_tools/BiblioSpec/src/Jamfile.jam

@@ -348,6 +348,17 @@ exe LibToSqlite3
     <library>blib_v_1
  ;
 
+unit-test-if-exists SpectrumTest
+  : # sources
+    SpectrumTest.cpp
+    blib
+    $(PWIZ_ROOT_PATH)/libraries/SQLite//sqlite3
+  : # requirements
+    <link>static
+    <threading>multi
+    <library>/ext/zlib//z
+ ;
+
 # for make targets below, <relevant> features make <location> be calculated properly
 
 alias all-relevant : : <relevant>cxxstd <relevant>cxxstd-dialect <relevant>target-os <relevant>toolset <relevant>variant <relevant>debug-symbols <relevant>runtime-link <relevant>link <relevant>address-model <relevant>architecture <relevant>asynch-exceptions <relevant>threading <relevant>debug-store <relevant>embed-manifest <relevant>exception-handling <relevant>extern-c-nothrow <relevant>inlining <relevant>instruction-set <relevant>optimization <relevant>pch <relevant>rtti <relevant>runtime-debugging <relevant>threading <relevant>warnings <relevant>warnings-as-errors <relevant>windows-api ;

pwiz_tools/BiblioSpec/src/PwizReader.cpp

@@ -163,6 +163,7 @@ bool PwizReader::getSpectrum(int identifier,
     }
     unique_ptr<SpectrumInfo> specInfo(new SpectrumInfo());
     specInfo->SpectrumInfo::update(*foundSpec, getPeaks);
+    ensureMzAscending(*specInfo);
 
     // confirm that it's an ms/ms spectrum
     if( getPeaks && specInfo->msLevel != 2 ){
@@ -303,6 +304,7 @@ bool PwizReader::getSpectrum(int identifier,
     }
     unique_ptr<SpectrumInfo> specInfo(new SpectrumInfo());
     specInfo->SpectrumInfo::update(*foundSpec, getPeaks);
+    ensureMzAscending(*specInfo);
 
     // confirm that it's an ms/ms spectrum
     if( specInfo->msLevel != 2 ){
@@ -426,11 +428,28 @@ void PwizReader::addCharges(BiblioSpec::Spectrum& returnSpectrum,
     } // next precursor
 }
 
+/**
+ * Put a spectrum's peaks in ascending m/z order if its writer did not. Ascending m/z is nowhere
+ * required of a writer, but everything downstream assumes it - and this has to happen here, on the
+ * SpectrumInfo, because it is the only point both transfer paths share. BlibBuild never builds a
+ * BiblioSpec::Spectrum at all: it goes BuildParser -> getSpectrum(int, SpecData&) -> transferSpec,
+ * and hands the raw arrays to insertPeaks, so a library built from a writer that presented some
+ * other order would be stored in that order and stay that way for every consumer of the .blib.
+ */
+void PwizReader::ensureMzAscending(SpectrumInfo& specInfo)
+{
+    if (is_sorted(specInfo.data.begin(), specInfo.data.end(),
+                  [](const MZIntensityPair& a, const MZIntensityPair& b) { return a.mz < b.mz; }))
+        return;
+    sort(specInfo.data.begin(), specInfo.data.end(),
+         [](const MZIntensityPair& a, const MZIntensityPair& b) { return a.mz < b.mz; });
+}
+
 /**
  * Copy the information from the Pwiz spectrum to the BiblioSpec
  * SpecData.
  */
-void PwizReader::transferSpec(BiblioSpec::SpecData& returnData, 
+void PwizReader::transferSpec(BiblioSpec::SpecData& returnData,
                               unique_ptr<SpectrumInfo>& specInfo){
 
     returnData.id = specInfo->scanNumber;

pwiz_tools/BiblioSpec/src/PwizReader.h

@@ -150,7 +150,11 @@ class PwizReader : public BiblioSpec::SpecFileReader {
      * Copy the information from the Pwiz spectrum to the BiblioSpec
      * SpecData.
      */
-    void transferSpec(BiblioSpec::SpecData& returnData, 
+    // Ascending m/z is nowhere required of a writer but is assumed everywhere downstream; both
+    // transfer paths below share the SpectrumInfo, so it is corrected there.
+    static void ensureMzAscending(SpectrumInfo& specInfo);
+
+    void transferSpec(BiblioSpec::SpecData& returnData,
                       unique_ptr<SpectrumInfo>& specInfo);
 
     /**

pwiz_tools/BiblioSpec/src/Spectrum.cpp

@@ -274,7 +274,23 @@ double Spectrum::getSignalToNoise() {
 
 // deletes existing peaks
 void Spectrum::setRawPeaks(const vector<PEAK_T>& newpeaks) {
-    rawPeaks_.assign(newpeaks.begin(), newpeaks.end()); 
+    rawPeaks_.assign(newpeaks.begin(), newpeaks.end());
+
+    // Ascending m/z is nowhere required of a writer, but everything downstream assumes it:
+    // binPeaks() merges a peak into the previous bin only when the two are adjacent in the array,
+    // and removePrecursorPeaks() binary searches the m/z axis. Peaks presented in some other order
+    // therefore leave same-bin peaks unmerged and erase the wrong range around the precursor, with
+    // nothing reported either way.
+    // This covers spectra arriving from a .blib and from Spectrum copies. Spectra read through
+    // PwizReader are ordered there instead, on the SpectrumInfo, because BlibBuild never builds a
+    // Spectrum at all - see PwizReader::ensureMzAscending.
+    // The check is a single pass; the sort only runs for a writer that did not present m/z order.
+    if (!is_sorted(rawPeaks_.begin(), rawPeaks_.end(),
+                   [](const PEAK_T& a, const PEAK_T& b) { return a.mz < b.mz; }))
+    {
+        sort(rawPeaks_.begin(), rawPeaks_.end(),
+             [](const PEAK_T& a, const PEAK_T& b) { return a.mz < b.mz; });
+    }
 }
 
 void Spectrum::setProcessedPeaks(const vector<PEAK_T>& newpeaks) {