src/RcppPwiz.cpp
84a4b166
 #include "RcppPwiz.h"
 
 RcppPwiz::RcppPwiz()
 {
8a535503
   msd = NULL;
2aed7135
   nativeIdFormat = CVID_Unknown;
8a535503
   instrumentInfo = Rcpp::List::create();
   chromatogramsInfo = Rcpp::DataFrame::create();
   isInCacheInstrumentInfo = FALSE;
   allScanHeaderInfo = Rcpp::List::create();
   isInCacheAllScanHeaderInfo = FALSE;
84a4b166
 }
d2f35232
 /* Destructor*/
 RcppPwiz::~RcppPwiz()
 {
   RcppPwiz::close();
 }
84a4b166
 
2aed7135
 // void RcppPwiz::open(const string& fileName)
 void RcppPwiz::open(Rcpp::StringVector fileName)
84a4b166
 {
 
2aed7135
   filename = Rcpp::as<std::string>(fileName(0));
   msd = new MSDataFile(filename);
4c03c487
   nativeIdFormat = id::getDefaultNativeIDFormat(*msd);
84a4b166
 }
d2f35232
 
 /* Release all memory on close. */
 void RcppPwiz::close()
 {
   if (msd != NULL)
     {
       delete msd;
       msd = NULL;
2aed7135
       nativeIdFormat = CVID_Unknown;
d2f35232
       instrumentInfo = Rcpp::List::create();
       chromatogramsInfo = Rcpp::DataFrame::create();
       isInCacheInstrumentInfo = FALSE;
       allScanHeaderInfo = Rcpp::List::create();
       isInCacheAllScanHeaderInfo = FALSE;
     }
 }
 
96ea799f
 string RcppPwiz::getFilename() {
8a535503
   return filename;
84a4b166
 }
 
96ea799f
 int RcppPwiz::getLastScan() const {
8a535503
   if (msd != NULL) {
     SpectrumListPtr slp = msd->run.spectrumListPtr;
     return slp->size();
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
8a535503
   return -1;
96ea799f
 }
 
 int RcppPwiz::getLastChrom() const {
8a535503
   if (msd != NULL) {
     ChromatogramListPtr clp = msd->run.chromatogramListPtr;
     return clp->size();
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
8a535503
   return -1;
84a4b166
 }
 
 Rcpp::List RcppPwiz::getInstrumentInfo ( )
 {
8a535503
   if (msd != NULL)
84a4b166
     {
8a535503
       if (!isInCacheInstrumentInfo)
84a4b166
         {
 
45834cf0
 	  vector<InstrumentConfigurationPtr> icp = msd->instrumentConfigurationPtrs; // NULL for mzData
8a535503
 	  if (icp.size() != 0)
84a4b166
             {
8a535503
 	      CVTranslator cvTranslator;
 	      LegacyAdapter_Instrument adapter(*icp[0], cvTranslator);
 	      vector<SoftwarePtr> sp = msd->softwarePtrs;
 	      std::vector<SamplePtr> sample = msd->samplePtrs;
 	      std::vector<ScanSettingsPtr> scansetting = msd->scanSettingsPtrs;
 	      std::string ionisation = "";
 	      std::string analyzer = "";
 	      std::string detector = "";
 	      // Fix issue #113
 	      // if (icp[0]->componentList.size() > 0)
 	      // That does still not mean we have a ionisation available.
 	      // Could be that either analyzer or detector or ionisation is
 	      // defined.
 	      // Have to use try-catch
 	      try {
45834cf0
 		ionisation = std::string(adapter.ionisation());
8a535503
 	      } catch(...) {}
 	      try {
45834cf0
 		analyzer = std::string(adapter.analyzer());
8a535503
 	      } catch(...) {}
 	      try {
45834cf0
 		detector = std::string(adapter.detector());
8a535503
 	      } catch(...) {}
 	      instrumentInfo = Rcpp::List::create(
 						  Rcpp::_["manufacturer"]  = std::string(adapter.manufacturer()),
 						  Rcpp::_["model"]         = std::string(adapter.model()),
 						  Rcpp::_["ionisation"]    = ionisation,
 						  Rcpp::_["analyzer"]      = analyzer,
 						  Rcpp::_["detector"]      = detector,
 						  Rcpp::_["software"]      = (sp.size()>0?sp[0]->id + " " + sp[0]->version:""),
 						  Rcpp::_["sample"]		  = (sample.size()>0?sample[0]->name+sample[0]->id:""),
 						  Rcpp::_["source"]        = (scansetting.size()>0?scansetting[0]->sourceFilePtrs[0]->location:"")
 						  ) ;
84a4b166
 
             }
8a535503
 	  else
84a4b166
             {
8a535503
 	      instrumentInfo = Rcpp::List::create(
 						  Rcpp::_["manufacturer"]  = "",
 						  Rcpp::_["model"]         = "",
 						  Rcpp::_["ionisation"]    = "",
 						  Rcpp::_["analyzer"]      = "",
 						  Rcpp::_["detector"]      = "",
 						  Rcpp::_["software"]      = "",
 						  Rcpp::_["sample"]		  = "",
 						  Rcpp::_["source"]		  = ""
 						  ) ;
84a4b166
             }
 
8a535503
 	  isInCacheInstrumentInfo = TRUE;
84a4b166
         }
8a535503
       return(instrumentInfo);
84a4b166
     }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
8a535503
   return instrumentInfo;
84a4b166
 }
 
2aed7135
 int RcppPwiz::getAcquisitionNumber(string id, size_t index) const
 {
4c03c487
   // const SpectrumIdentity& si = msd->run.spectrumListPtr->spectrumIdentity(index);
   string scanNumber = id::translateNativeIDToScanNumber(nativeIdFormat, id);
   if (scanNumber.empty())
     return static_cast<int>(index) + 1;
   else
     return lexical_cast<int>(scanNumber);
2aed7135
 }
 
4c03c487
 // Using this function instead of pwiz translateNativeIDToScanNumber because
 // it randomly causes segfaults on macOS.
 // int RcppPwiz::getAcquisitionNumber(string id, size_t index) const
 // {
 //   if (id.find("controllerType") != std::string::npos) {
 //     if (id.find("controllerType=0 controllerNumber=1") == std::string::npos)
 //       return static_cast<int>(index) + 1;
 //   }
 //   string e;
 //   std::smatch match;
 //   if (id.find("scan=") != std::string::npos)
 //     e ="scan=(\\d+)";
 //   else if (id.find("index=") != std::string::npos)
 //     e = "index=(\\d+)";
 //   else if (id.find("spectrum=") != std::string::npos)
 //     e = "spectrum=(\\d+)";
 //   else if (id.find("scanId=") != std::string::npos)
 //     e = "scanId=(\\d+)";
 //   else return static_cast<int>(index) + 1;
 //   if (std::regex_search(id, match, std::regex(e)))
 //     return lexical_cast<int>(match[1]);
 //   else return static_cast<int>(index) + 1;
 // }
 
32246a88
 Rcpp::DataFrame RcppPwiz::getScanHeaderInfo (Rcpp::IntegerVector whichScan) {
   if (msd != NULL) {
     SpectrumListPtr slp = msd->run.spectrumListPtr;
2aed7135
     size_t N = slp->size();
     size_t N_scans = whichScan.size();
32246a88
     Rcpp::IntegerVector seqNum(N_scans); // number in sequence observed file (1-based)
     Rcpp::IntegerVector acquisitionNum(N_scans); // scan number as declared in File (may be gaps)
     Rcpp::IntegerVector msLevel(N_scans);
     Rcpp::IntegerVector polarity(N_scans);
     Rcpp::IntegerVector peaksCount(N_scans);
     Rcpp::NumericVector totIonCurrent(N_scans);
     Rcpp::NumericVector retentionTime(N_scans);        /* in seconds */
     Rcpp::NumericVector basePeakMZ(N_scans);
     Rcpp::NumericVector basePeakIntensity(N_scans);
     Rcpp::NumericVector collisionEnergy(N_scans);
     Rcpp::NumericVector ionisationEnergy(N_scans);
     Rcpp::NumericVector lowMZ(N_scans);
     Rcpp::NumericVector highMZ(N_scans);
     Rcpp::IntegerVector precursorScanNum(N_scans); /* only if MS level > 1 */
     Rcpp::NumericVector precursorMZ(N_scans);  /* only if MS level > 1 */
     Rcpp::IntegerVector precursorCharge(N_scans);  /* only if MS level > 1 */
     Rcpp::NumericVector precursorIntensity(N_scans);  /* only if MS level > 1 */
     Rcpp::IntegerVector mergedScan(N_scans);  /* only if MS level > 1 */
     Rcpp::IntegerVector mergedResultScanNum(N_scans); /* scan number of the resultant merged scan */
     Rcpp::IntegerVector mergedResultStartScanNum(N_scans); /* smallest scan number of the scanOrigin for merged scan */
     Rcpp::IntegerVector mergedResultEndScanNum(N_scans); /* largest scan number of the scanOrigin for merged scan */
     Rcpp::NumericVector ionInjectionTime(N_scans); /* The time spent filling an ion trapping device*/
     Rcpp::StringVector filterString(N_scans);
     Rcpp::StringVector spectrumId(N_scans);
     Rcpp::LogicalVector centroided(N_scans);
     Rcpp::NumericVector ionMobilityDriftTime(N_scans);
     Rcpp::NumericVector isolationWindowTargetMZ(N_scans);
     Rcpp::NumericVector isolationWindowLowerOffset(N_scans);
     Rcpp::NumericVector isolationWindowUpperOffset(N_scans);
     Rcpp::NumericVector scanWindowLowerLimit(N_scans);
     Rcpp::NumericVector scanWindowUpperLimit(N_scans);
45834cf0
 
2aed7135
     for (size_t i = 0; i < N_scans; i++) {
32246a88
       int current_scan = whichScan[i];
2aed7135
       size_t current_index = static_cast<size_t>(current_scan - 1);
       // SpectrumPtr sp = slp->spectrum(current_index, false);
       SpectrumPtr sp = slp->spectrum(current_index, DetailLevel_FullMetadata);
32246a88
       Scan dummy;
       Scan& scan = sp->scanList.scans.empty() ? dummy : sp->scanList.scans[0];
2aed7135
       if (scan.empty())
 	Rprintf("Scan with index %d empty.\n", current_scan);
32246a88
       // seqNum
       seqNum[i] = current_scan;
2aed7135
       acquisitionNum[i] = getAcquisitionNumber(sp->id, current_index);
32246a88
       // spectrumId
2aed7135
       spectrumId[i] = Rcpp::String(sp->id);
32246a88
       // msLevel
       msLevel[i] = sp->cvParam(MS_ms_level).valueAs<int>();
       // peaksCount
       peaksCount[i] = static_cast<int>(sp->defaultArrayLength);
       // totIonCurrent
       totIonCurrent[i] = sp->cvParam(MS_total_ion_current).valueAs<double>();
       // basePeakMZ
       basePeakMZ[i] = sp->cvParam(MS_base_peak_m_z).valueAs<double>();
       // basePeakIntensity
       basePeakIntensity[i] = sp->cvParam(MS_base_peak_intensity).valueAs<double>();
       // ionisationEnerty
       ionisationEnergy[i] = sp->cvParam(MS_ionization_energy_OBSOLETE).valueAs<double>();
       // lowMZ
       lowMZ[i] = sp->cvParam(MS_lowest_observed_m_z).valueAs<double>();
       // highMZ
       highMZ[i] = sp->cvParam(MS_highest_observed_m_z).valueAs<double>();
       // polarity
       CVParam param = sp->cvParamChild(MS_scan_polarity);
       polarity[i] = (param.cvid==MS_negative_scan ? 0 : (param.cvid==MS_positive_scan ? +1 : -1 ) );
       // centroided
       param = sp->cvParamChild(MS_spectrum_representation);
45834cf0
       centroided[i] = (param.cvid==MS_centroid_spectrum ? TRUE : (param.cvid==MS_profile_spectrum ? FALSE : NA_LOGICAL));
32246a88
       // retentionTime
       retentionTime[i] = scan.cvParam(MS_scan_start_time).timeInSeconds();
       // ionInjectionTime
       ionInjectionTime[i] = (scan.cvParam(MS_ion_injection_time).timeInSeconds() * 1000);
       // filterString
       filterString[i] = scan.cvParam(MS_filter_string).value.empty() ? NA_STRING : Rcpp::String(scan.cvParam(MS_filter_string).value);
       // ionMobilityDriftTime
       ionMobilityDriftTime[i] = scan.cvParam(MS_ion_mobility_drift_time).value.empty() ? NA_REAL : (scan.cvParam(MS_ion_mobility_drift_time).timeInSeconds() * 1000);
       // scanWindowLowerLimit and scanWindowUpperLimit
       if (!scan.scanWindows.empty()) {
 	scanWindowLowerLimit[i] = scan.scanWindows[0].cvParam(MS_scan_window_lower_limit).valueAs<double>();
 	scanWindowUpperLimit[i] = scan.scanWindows[0].cvParam(MS_scan_window_upper_limit).valueAs<double>();
       } else {
 	scanWindowLowerLimit[i] = NA_REAL;
 	scanWindowUpperLimit[i] = NA_REAL;
       }
       // mergedScan - also not supported by RAMPAdapter
       mergedScan[i] = NA_INTEGER;
       mergedResultScanNum[i] = NA_INTEGER;
       mergedResultStartScanNum[i] = NA_INTEGER;
       mergedResultEndScanNum[i] = NA_INTEGER;
       if (!sp->precursors.empty()) {
 	const Precursor& precursor = sp->precursors[0];
 	// collisionEnergy
 	collisionEnergy[i] = precursor.activation.cvParam(MS_collision_energy).valueAs<double>();
 	// precursorScanNum
 	size_t precursorIndex = slp->find(precursor.spectrumID);
2aed7135
 	if (precursorIndex < N) {
 	  precursorScanNum[i] = getAcquisitionNumber(precursor.spectrumID, precursorIndex);
32246a88
 	} else {
 	  precursorScanNum[i] = NA_INTEGER;
606b1d91
 	}
32246a88
 	// precursorMZ, precursorCharge, precursorIntensity
 	if (!precursor.selectedIons.empty()) {
 	  precursorMZ[i] = precursor.selectedIons[0].cvParam(MS_selected_ion_m_z).value.empty() ? precursor.selectedIons[0].cvParam(MS_m_z).valueAs<double>() : precursor.selectedIons[0].cvParam(MS_selected_ion_m_z).valueAs<double>();
 	  precursorCharge[i] = precursor.selectedIons[0].cvParam(MS_charge_state).valueAs<int>();
 	  precursorIntensity[i] = precursor.selectedIons[0].cvParam(MS_peak_intensity).valueAs<double>();
 	}
 	// isolationWindowTargetMZ, ...
 	IsolationWindow iwin = sp->precursors[0].isolationWindow;
 	if (!iwin.empty()) {
 	  isolationWindowTargetMZ[i] = iwin.cvParam(MS_isolation_window_target_m_z).value.empty() ? NA_REAL : iwin.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
 	  isolationWindowLowerOffset[i] = iwin.cvParam(MS_isolation_window_lower_offset).value.empty() ? NA_REAL : iwin.cvParam(MS_isolation_window_lower_offset).valueAs<double>();
 	  isolationWindowUpperOffset[i] = iwin.cvParam(MS_isolation_window_upper_offset).value.empty() ? NA_REAL : iwin.cvParam(MS_isolation_window_upper_offset).valueAs<double>();
 	} else {
 	  isolationWindowTargetMZ[i] = NA_REAL;
 	  isolationWindowLowerOffset[i] = NA_REAL;
 	  isolationWindowUpperOffset[i] = NA_REAL;
 	}
       } else {
 	collisionEnergy[i] = NA_REAL;
 	precursorScanNum[i] = NA_INTEGER;
 	precursorMZ[i] = NA_REAL;
 	precursorCharge[i] = NA_INTEGER;
 	precursorIntensity[i] = NA_REAL;
 	mergedScan[i] = NA_INTEGER;
 	mergedResultScanNum[i] = NA_INTEGER;
 	mergedResultStartScanNum[i] = NA_INTEGER;
 	mergedResultEndScanNum[i] = NA_INTEGER;
 	isolationWindowTargetMZ[i] = NA_REAL;
 	isolationWindowLowerOffset[i] = NA_REAL;
 	isolationWindowUpperOffset[i] = NA_REAL;
       }
84a4b166
     }
45834cf0
 
32246a88
     Rcpp::List header(31);
     std::vector<std::string> names;
2aed7135
     size_t i = 0;
32246a88
     names.push_back("seqNum");
     header[i++] = Rcpp::wrap(seqNum);
     names.push_back("acquisitionNum");
     header[i++] = Rcpp::wrap(acquisitionNum);
     names.push_back("msLevel");
     header[i++] = Rcpp::wrap(msLevel);
     names.push_back("polarity");
     header[i++] = Rcpp::wrap(polarity);
     names.push_back("peaksCount");
     header[i++] = Rcpp::wrap(peaksCount);
     names.push_back("totIonCurrent");
     header[i++] = Rcpp::wrap(totIonCurrent);
     names.push_back("retentionTime");
     header[i++] = Rcpp::wrap(retentionTime);
     names.push_back("basePeakMZ");
     header[i++] = Rcpp::wrap(basePeakMZ);
     names.push_back("basePeakIntensity");
     header[i++] = Rcpp::wrap(basePeakIntensity);
     names.push_back("collisionEnergy");
     header[i++] = Rcpp::wrap(collisionEnergy);
     names.push_back("ionisationEnergy");
     header[i++] = Rcpp::wrap(ionisationEnergy);
     names.push_back("lowMZ");
     header[i++] = Rcpp::wrap(lowMZ);
     names.push_back("highMZ");
     header[i++] = Rcpp::wrap(highMZ);
     names.push_back("precursorScanNum");
     header[i++] = Rcpp::wrap(precursorScanNum);
     names.push_back("precursorMZ");
     header[i++] = Rcpp::wrap(precursorMZ);
     names.push_back("precursorCharge");
     header[i++] = Rcpp::wrap(precursorCharge);
     names.push_back("precursorIntensity");
     header[i++] = Rcpp::wrap(precursorIntensity);
     names.push_back("mergedScan");
     header[i++] = Rcpp::wrap(mergedScan);
     names.push_back("mergedResultScanNum");
     header[i++] = Rcpp::wrap(mergedResultScanNum);
     names.push_back("mergedResultStartScanNum");
     header[i++] = Rcpp::wrap(mergedResultStartScanNum);
     names.push_back("mergedResultEndScanNum");
     header[i++] = Rcpp::wrap(mergedResultEndScanNum);
     names.push_back("injectionTime");
     header[i++] = Rcpp::wrap(ionInjectionTime);
     names.push_back("filterString");
     header[i++] = Rcpp::wrap(filterString);
     names.push_back("spectrumId");
     header[i++] = Rcpp::wrap(spectrumId);
     names.push_back("centroided");
     header[i++] = Rcpp::wrap(centroided);
     names.push_back("ionMobilityDriftTime");
45834cf0
     header[i++] = Rcpp::wrap(ionMobilityDriftTime);
32246a88
     names.push_back("isolationWindowTargetMZ");
45834cf0
     header[i++] = Rcpp::wrap(isolationWindowTargetMZ);
32246a88
     names.push_back("isolationWindowLowerOffset");
45834cf0
     header[i++] = Rcpp::wrap(isolationWindowLowerOffset);
32246a88
     names.push_back("isolationWindowUpperOffset");
45834cf0
     header[i++] = Rcpp::wrap(isolationWindowUpperOffset);
32246a88
     names.push_back("scanWindowLowerLimit");
45834cf0
     header[i++] = Rcpp::wrap(scanWindowLowerLimit);
32246a88
     names.push_back("scanWindowUpperLimit");
45834cf0
     header[i++] = Rcpp::wrap(scanWindowUpperLimit);
32246a88
     header.attr("names") = names;
45834cf0
 
32246a88
     return header;
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
8a535503
   return Rcpp::DataFrame::create( );
84a4b166
 }
 
32246a88
 Rcpp::DataFrame RcppPwiz::getAllScanHeaderInfo ( ) {
   if (msd != NULL) {
     if (!isInCacheAllScanHeaderInfo) {
       SpectrumListPtr slp = msd->run.spectrumListPtr;
2aed7135
       size_t N = slp->size();
45834cf0
 
2c7479f5
       allScanHeaderInfo = getScanHeaderInfo(Rcpp::seq(1, N));
45834cf0
       isInCacheAllScanHeaderInfo = TRUE;
84a4b166
     }
32246a88
     return allScanHeaderInfo ;
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
8a535503
   return Rcpp::DataFrame::create( );
84a4b166
 }
 
5bcb9d2b
 Rcpp::List RcppPwiz::getPeakList(Rcpp::IntegerVector whichScan) {
32246a88
   if (msd != NULL) {
     SpectrumListPtr slp = msd->run.spectrumListPtr;
2aed7135
     size_t n_scans = slp->size();
     size_t n_want = whichScan.size();
5bcb9d2b
     int current_scan;
     SpectrumPtr sp;
     BinaryDataArrayPtr mzs,ints;
     std::vector<double> data;
     Rcpp::NumericVector data_matrix;
     Rcpp::List res(n_want);
c2ab2b8c
     Rcpp::CharacterVector cn = Rcpp::CharacterVector::create("mz", "intensity");
2aed7135
     for (size_t i = 0; i < n_want; i++) {
5bcb9d2b
       current_scan = whichScan[i];
       if (current_scan < 1 || current_scan > n_scans) {
 	Rprintf("Index whichScan out of bounds [1 ... %d].\n", n_scans);
 	return Rcpp::List::create( );
32246a88
       }
2aed7135
       size_t current_index = static_cast<size_t>(current_scan - 1);
       // sp = slp->spectrum(current_index, true);
       sp = slp->spectrum(current_index, DetailLevel_FullData);
5bcb9d2b
       mzs = sp->getMZArray();
       ints = sp->getIntensityArray();
       if (!mzs.get() || !ints.get()) {
 	Rcpp::NumericMatrix pks(0, 2);
c2ab2b8c
 	Rcpp::colnames(pks) = cn;
5bcb9d2b
 	res[i] = pks;
 	continue;
       }
       if (mzs->data.size() != ints->data.size())
 	Rcpp::Rcerr << "Sizes of mz and intensity arrays don't match." << std::endl;
       data = mzs->data;
       data.insert(data.end(), ints->data.begin(), ints->data.end());
       data_matrix = Rcpp::wrap(data);
       data_matrix.attr("dim") = Rcpp::Dimension(ints->data.size(), 2);
c2ab2b8c
       Rcpp::colnames(data_matrix) = cn;
5bcb9d2b
       res[i] = data_matrix;
84a4b166
     }
5bcb9d2b
     return res;
32246a88
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
5bcb9d2b
   return Rcpp::List::create();
84a4b166
 }
 
19dc1af2
 /**
  * copyWriteMSFile copies (general) content from the originating MS file and
  * replaces the Spectrum list with the new data provided with arguments
  * spctr_header and spctr_data.
  * TODO: add strings that describe what processings have been done in R.
  * We're copying:
  * o fileDescription (adding also the originating MS file)
  * o softwareList (adding also additional info)
  * o instrumentConfigurationList
  * o dataProcessingList (adding also additional info)
  * o run: all "general" info from the run.
  * Potential additional parameters:
  * - centroided: whether spectra data is centroided: NA, TRUE, FALSE.
  * - processings: string/vector with all processing steps.
  * - software(s): name (and version?) of software.
5afbbd4f
  * software_processing: has to be a list of character vectors.
19dc1af2
  **/
e501fc3e
 void RcppPwiz::copyWriteMSfile(const string& file, const string& format,
 			       const string& originalFile,
 			       Rcpp::DataFrame spctr_header,
 			       Rcpp::List spctr_data,
5afbbd4f
 			       bool rtime_seconds,
 			       Rcpp::List software_processing) {
e501fc3e
   MSDataFile *msd;
   msd = new MSDataFile(originalFile);
   MSData newmsd;
   newmsd.cvs = defaultCVList();
 
   Rcpp::IntegerVector msLevel = spctr_header["msLevel"];
   // Copy data from the original file.
   // o fileDescription with: fileContent, sourceFileList
d4bcbe13
   //   NOTE: don't copy the file description for mzXML export - somehow the
   //   spectra data will then not be written.
   if (format != "mzxml")
     newmsd.fileDescription = msd->fileDescription;
e501fc3e
   bool is_ms1 = false;
   bool is_msn = false;
   for (int i = 0; i < msLevel.size(); i++) {
     if (msLevel[i] == 1)
       is_ms1 = true;
     if (msLevel[i] > 1)
       is_msn = true;
   }
   if (is_ms1)
     newmsd.fileDescription.fileContent.set(MS_MS1_spectrum);
   if (is_msn)
     newmsd.fileDescription.fileContent.set(MS_MSn_spectrum);
 
5afbbd4f
   // o paramGroupList
d4bcbe13
   if (format != "mzxml")
     newmsd.paramGroupPtrs = msd->paramGroupPtrs;
19dc1af2
   // o sampleList
e501fc3e
   newmsd.samplePtrs = msd->samplePtrs;
5afbbd4f
   // o instrumentConfigurationList
   newmsd.instrumentConfigurationPtrs = msd->instrumentConfigurationPtrs;
   // o softwareList
   newmsd.softwarePtrs = msd->softwarePtrs;
e501fc3e
   // o dataProcessingList
   newmsd.dataProcessingPtrs = msd->dataProcessingPtrs;
5afbbd4f
   // Add new software and processings:
   if (software_processing.size() > 0) {
     for (int sp = 0; sp < software_processing.size(); sp++) {
       addDataProcessing(newmsd, Rcpp::as<Rcpp::StringVector>(software_processing(sp)));
     }
19dc1af2
   }
45834cf0
 
5afbbd4f
   // o run
e501fc3e
   // Initialize the run and fill with data from the original file.
   Run &original_run = msd->run;
   newmsd.run.id = original_run.id;
d4bcbe13
   if (format != "mzxml") {
     newmsd.run.defaultInstrumentConfigurationPtr =
       original_run.defaultInstrumentConfigurationPtr;
     newmsd.run.samplePtr = original_run.samplePtr;
     newmsd.run.startTimeStamp = original_run.startTimeStamp;
     newmsd.run.defaultSourceFilePtr = original_run.defaultSourceFilePtr;
   }
e501fc3e
   // Now filling with new data
   addSpectrumList(newmsd, spctr_header, spctr_data, rtime_seconds);
45834cf0
 
d4bcbe13
   if (format == "mgf") {
5afbbd4f
     std::ofstream* mgfOutFileP = new std::ofstream(file.c_str());
     Serializer_MGF serializerMGF;
     serializerMGF.write(*mgfOutFileP, newmsd);
     mgfOutFileP->flush();
     mgfOutFileP->close();
d4bcbe13
   } else if (format == "mzxml") {
5afbbd4f
     std::ofstream mzXMLOutFileP(file.c_str());
     Serializer_mzXML::Config config;
     config.binaryDataEncoderConfig.compression = BinaryDataEncoder::Compression_Zlib;
     Serializer_mzXML serializerMzXML(config);
     serializerMzXML.write(mzXMLOutFileP, newmsd);
     mzXMLOutFileP.flush();
     mzXMLOutFileP.close();
d4bcbe13
   } else if (format == "mzml") {
5afbbd4f
     std::ofstream mzXMLOutFileP(file.c_str());
     Serializer_mzML::Config config;
     config.binaryDataEncoderConfig.compression = BinaryDataEncoder::Compression_Zlib;
     Serializer_mzML mzmlSerializer(config);
     mzmlSerializer.write(mzXMLOutFileP, newmsd);
     mzXMLOutFileP.flush();
     mzXMLOutFileP.close();
   }
e501fc3e
   else
     Rcpp::Rcerr << format << " format not supported! Please try mgf, mzML, mzXML or mz5." << std::endl;
45834cf0
 
e501fc3e
   // Cleanup.
   delete msd;
 }
 
 // writeSpectrumList: writes the provided spectrum data to a file.
 void RcppPwiz::writeSpectrumList(const string& file, const string& format,
 				 Rcpp::DataFrame spctr_header,
 				 Rcpp::List spctr_data,
5afbbd4f
 				 bool rtime_seconds,
 				 Rcpp::List software_processing) {
e501fc3e
   MSData newmsd;
   newmsd.cvs = defaultCVList();
 
   Rcpp::IntegerVector msLevel = spctr_header["msLevel"];
   bool is_ms1 = false;
   bool is_msn = false;
   for (int i = 0; i < msLevel.size(); i++) {
     if (msLevel[i] == 1)
       is_ms1 = true;
     if (msLevel[i] > 1)
       is_msn = true;
   }
   if (is_ms1)
     newmsd.fileDescription.fileContent.set(MS_MS1_spectrum);
   if (is_msn)
     newmsd.fileDescription.fileContent.set(MS_MSn_spectrum);
5afbbd4f
 
   // Add software_processing:
   if (software_processing.size() > 0) {
     for (int sp = 0; sp < software_processing.size(); sp++) {
       addDataProcessing(newmsd, Rcpp::as<Rcpp::StringVector>(software_processing(sp)));
     }
19dc1af2
   }
45834cf0
 
d4bcbe13
   newmsd.run.id = "Experiment_1";
e501fc3e
 
   // Now filling with new data
   addSpectrumList(newmsd, spctr_header, spctr_data, rtime_seconds);
 
d4bcbe13
   if (format == "mgf") {
5afbbd4f
     std::ofstream* mgfOutFileP = new std::ofstream(file.c_str());
     Serializer_MGF serializerMGF;
     serializerMGF.write(*mgfOutFileP, newmsd);
     mgfOutFileP->flush();
     mgfOutFileP->close();
d4bcbe13
   } else if (format == "mzxml") {
5afbbd4f
     std::ofstream mzXMLOutFileP(file.c_str());
     Serializer_mzXML::Config config;
     config.binaryDataEncoderConfig.compression = BinaryDataEncoder::Compression_Zlib;
     Serializer_mzXML serializerMzXML(config);
     serializerMzXML.write(mzXMLOutFileP, newmsd);
     mzXMLOutFileP.flush();
     mzXMLOutFileP.close();
d4bcbe13
   } else if (format == "mzml") {
5afbbd4f
     std::ofstream mzXMLOutFileP(file.c_str());
     Serializer_mzML::Config config;
     config.binaryDataEncoderConfig.compression = BinaryDataEncoder::Compression_Zlib;
     Serializer_mzML mzmlSerializer(config);
     mzmlSerializer.write(mzXMLOutFileP, newmsd);
     mzXMLOutFileP.flush();
     mzXMLOutFileP.close();
   }
e501fc3e
   else
     Rcpp::Rcerr << format << " format not supported! Please try mgf, mzML, mzXML or mz5." << std::endl;
 }
 
5afbbd4f
 /*
  * o soft_proc: is supposed to be a character vector of length >= 2:
  *   soft_proc[0]: The software name (required).
  *   soft_proc[1]: The software version (required).
45834cf0
  *   soft_proc[2]: The CV ID of the software. Use "MS:-1" if not known, in
82debcef
  *                 which case we are NOT writing the corresponding CV element.
45834cf0
  *   soft_proc[3-length]: CV IDs of the processing steps (optional).
5afbbd4f
  */
 void RcppPwiz::addDataProcessing(MSData& msd, Rcpp::StringVector soft_proc) {
19dc1af2
   SoftwarePtr new_soft(new Software);
5afbbd4f
   new_soft->id = soft_proc(0);
   new_soft->version = soft_proc(1);
   int soft_proc_size = soft_proc.size();
   if (soft_proc_size > 2) {
82debcef
     if (soft_proc(2) != "MS:-1") {
       CVTermInfo cv_term = cvTermInfo(soft_proc(2));
       new_soft->set(cv_term.cvid);
     }
19dc1af2
   }
5afbbd4f
   // Order: get the number of already present dataProcessingPtrs and
   // increment
   int order = msd.dataProcessingPtrs.size() + 1;
   DataProcessingPtr data_processing(new DataProcessing);
   std::ostringstream oss;
   oss << soft_proc[0] << "_processing";
   data_processing->id = oss.str();
   ProcessingMethod proc_meth;
   proc_meth.order = order;
   proc_meth.softwarePtr = new_soft;
   if (soft_proc_size > 3) {
     // Got also processing steps.
     for (int i = 3; i < soft_proc_size; i++) {
       CVTermInfo cv_term = cvTermInfo(soft_proc(i));
       proc_meth.set(cv_term.cvid);
19dc1af2
     }
ed65fd49
   }
5afbbd4f
   data_processing->processingMethods.push_back(proc_meth);
   msd.softwarePtrs.push_back(new_soft);
   msd.dataProcessingPtrs.push_back(data_processing);
ed65fd49
 }
5afbbd4f
 
ed65fd49
 /** Adds information provided in the header and spectra data to the spectrumList
  *  content of the MSData.
  *  TODO: OPEN QUESTION: what to use as spectrum ID? See issue #105
  *       For now: use scan=acquisitionNum[i]. According to the code of the
  *       RAMPAdapter.cpp this seems to be correct - the scan number (i.e.
  *       acquisitionNum) is extracted/guessed from the id of the spectrum.
  *       Need to test: what if the acquisitionNum has gaps? Are MSn spectra
  *       still linked correctly to their precursor?
  *       Alternative: scan=seqNum[i].
  **/
e501fc3e
 void RcppPwiz::addSpectrumList(MSData& msd,
 			       Rcpp::DataFrame& spctr_header,
 			       Rcpp::List& spctr_data,
 			       bool rtime_seconds) {
630000b7
   int precursor_idx;
e501fc3e
   // Break the header down into its elements/columns:
   Rcpp::IntegerVector seqNum = spctr_header["seqNum"];
   Rcpp::IntegerVector acquisitionNum = spctr_header["acquisitionNum"];
   Rcpp::IntegerVector msLevel = spctr_header["msLevel"];
   Rcpp::IntegerVector polarity = spctr_header["polarity"];
   Rcpp::IntegerVector peaksCount = spctr_header["peaksCount"];
5afbbd4f
   Rcpp::NumericVector totIonCurrent = spctr_header["totIonCurrent"];
e501fc3e
   Rcpp::NumericVector retentionTime = spctr_header["retentionTime"];
   Rcpp::NumericVector basePeakMZ = spctr_header["basePeakMZ"];
   Rcpp::NumericVector basePeakIntensity = spctr_header["basePeakIntensity"];
   Rcpp::NumericVector collisionEnergy = spctr_header["collisionEnergy"];
   Rcpp::NumericVector ionisationEnergy = spctr_header["ionisationEnergy"];
   Rcpp::NumericVector lowMZ = spctr_header["lowMZ"];
   Rcpp::NumericVector highMZ = spctr_header["highMZ"];
   Rcpp::IntegerVector precursorScanNum = spctr_header["precursorScanNum"];
   Rcpp::NumericVector precursorMZ = spctr_header["precursorMZ"];
   Rcpp::IntegerVector precursorCharge = spctr_header["precursorCharge"];
   Rcpp::NumericVector precursorIntensity = spctr_header["precursorIntensity"];
5afbbd4f
   Rcpp::IntegerVector mergedScan = spctr_header["mergedScan"];
e501fc3e
   // Skipping mergedResultScanNum, mergedResultStartScanNum and mergedResultEndScanNum
5afbbd4f
   Rcpp::NumericVector ionInjectionTime = spctr_header["injectionTime"];
1ce0075c
   Rcpp::StringVector filterString = spctr_header["filterString"];
630000b7
   Rcpp::StringVector spectrumId = spctr_header["spectrumId"];
a4b3e08f
   Rcpp::LogicalVector centroided = spctr_header["centroided"];
817d1420
   Rcpp::NumericVector ionMobilityDriftTime = spctr_header["ionMobilityDriftTime"];
2edb9032
   Rcpp::NumericVector isolationWindowTargetMZ = spctr_header["isolationWindowTargetMZ"];
   Rcpp::NumericVector isolationWindowLowerOffset = spctr_header["isolationWindowLowerOffset"];
   Rcpp::NumericVector isolationWindowUpperOffset = spctr_header["isolationWindowUpperOffset"];
52172a4d
   Rcpp::NumericVector scanWindowLowerLimit = spctr_header["scanWindowLowerLimit"];
   Rcpp::NumericVector scanWindowUpperLimit = spctr_header["scanWindowUpperLimit"];
45834cf0
 
e501fc3e
   // From MSnbase::Spectrum        Column in the header
   // msLevel integer               $msLevel
   // peaksCount integer
   // rt numeric
   // acquisitionNum integer        $acquisitionNum
   // scanIndex integer             $seqNum
   // tic numeric                   $totIonCurrent
   // mz numeric                    peaks()[, 1]
   // intensity numeric             peaks()[, 2]
   // fromFile integer
   // centroided logical
   // smoothed logical
   // polarity integer              $polarity: 0 negative, 1 positive, -1 unknown
   // Spectrum2
   // merged numeric                $mergedScan
   // precScanNum integer           $precursorScanNum
   // precursorMz numeric           $precursorMz
   // precursorIntensity numeric    $precursorIntensity
   // precursorCharge integer       $precursorCharge
   // collisionEnergy numeric       $collisionEnergy
45834cf0
 
e501fc3e
   // Now filling with new data
   shared_ptr<SpectrumListSimple> spectrumList(new SpectrumListSimple);
30e00142
   // Add the default Processing pointer (fix issue #151
   spectrumList->dp = msd.dataProcessingPtrs[(msd.dataProcessingPtrs.size() - 1)];
e501fc3e
   msd.run.spectrumListPtr = spectrumList;
   // TODO add also eventual processings.
   for (int i = 0; i < spctr_data.size(); i++) {
     spectrumList->spectra.push_back(SpectrumPtr(new Spectrum));
     Spectrum& spct = *spectrumList->spectra[i];
     spct.set(MS_ms_level, msLevel[i]);
8a535503
     // centroided
3b091f3a
     if (!Rcpp::LogicalVector::is_na(centroided[i]) && centroided[i] == TRUE)
8a535503
       spct.set(MS_centroid_spectrum);
3b091f3a
     if (!Rcpp::LogicalVector::is_na(centroided[i]) && centroided[i] == FALSE)
8a535503
       spct.set(MS_profile_spectrum);
e501fc3e
     // [X] polarity
     if (polarity[i] == 0)
       spct.set(MS_negative_scan);
     if (polarity[i] == 1)
       spct.set(MS_positive_scan);
     if (msLevel[i] == 1)
       spct.set(MS_MS1_spectrum);
     else
       spct.set(MS_MSn_spectrum);
8a535503
     spct.set(MS_lowest_observed_m_z, lowMZ[i], MS_m_z);
     spct.set(MS_highest_observed_m_z, highMZ[i], MS_m_z);
     spct.set(MS_base_peak_m_z, basePeakMZ[i], MS_m_z);
     spct.set(MS_base_peak_intensity, basePeakIntensity[i],
 	     MS_number_of_detector_counts);
e501fc3e
     spct.set(MS_total_ion_current, totIonCurrent[i]);
     // TODO:
     // [X] seqNum: number observed in file.
630000b7
     spct.index = seqNum[i] - 1;	// Or just i?
e501fc3e
     // [X] acquisitionNum: number as reported (there might be gaps).
630000b7
     // spct.id = "scan=" + boost::lexical_cast<std::string>(acquisitionNum[i]);
     spct.id = spectrumId[i];	// Use the provided ID instead
e501fc3e
     // [ ] peaksCount: no need to set this?
     // [X] retentionTime
     spct.scanList.scans.push_back(Scan());
     spct.scanList.set(MS_no_combination);
     Scan &spct_scan = spct.scanList.scans.back();
8a535503
     if (rtime_seconds)
e501fc3e
       spct_scan.set(MS_scan_start_time, retentionTime[i], UO_second);
8a535503
     else
e501fc3e
       spct_scan.set(MS_scan_start_time, retentionTime[i], UO_minute);
8a535503
     if (ionInjectionTime[i] > 0)
       spct_scan.set(MS_ion_injection_time, ionInjectionTime[i], UO_millisecond);
b5ca4753
 
     if (!Rcpp::StringVector::is_na(filterString[i]))
       spct_scan.set(MS_filter_string, filterString[i]);
 
3b091f3a
     if (!Rcpp::NumericVector::is_na(ionMobilityDriftTime[i]))
817d1420
       spct_scan.set(MS_ion_mobility_drift_time, ionMobilityDriftTime[i],
 		    UO_millisecond);
52172a4d
 
     // scanWindow
3b091f3a
     if (!Rcpp::NumericVector::is_na(scanWindowLowerLimit[i]) &&
 	!Rcpp::NumericVector::is_na(scanWindowUpperLimit[i])) {
52172a4d
       spct_scan.scanWindows.push_back(ScanWindow(scanWindowLowerLimit[i], scanWindowUpperLimit[i], MS_m_z));
     }
e501fc3e
     // MSn - precursor:
45834cf0
     if ( (precursorScanNum[i] > 0) | (precursorMZ[i] > 0) ) {
6e9d806c
       // Fill precursor data. This preserves the precursor data even if the
       // precursor scan is not available (e.g. after MS level filtering).
       spct.precursors.resize(1);
       Precursor& prec = spct.precursors.front();
       if (collisionEnergy[i] != 0) {
 	prec.activation.set(MS_collision_induced_dissociation);
 	prec.activation.set(MS_collision_energy, collisionEnergy[i],
 			    UO_electronvolt);
       }
       prec.selectedIons.resize(1);
       prec.selectedIons[0].set(MS_selected_ion_m_z, precursorMZ[i], MS_m_z);
       prec.selectedIons[0].set(MS_peak_intensity, precursorIntensity[i],
 			       MS_number_of_detector_counts);
       prec.selectedIons[0].set(MS_charge_state, precursorCharge[i]);
630000b7
       // Get the spectrumId of the precursor. Assuming that precursorScanNum is
       // linked to the acquisitionNum of the precursor.
32234431
       // This seems to be correct, since both the acquisitionNum and the
       // precursorNum are extracted from the respective spectrum's ID.
       precursor_idx = -1;
630000b7
       for (int j = 0; j < spctr_data.size(); j++) {
 	if (precursorScanNum[i] == acquisitionNum[j]) {
 	  precursor_idx = j;
 	  break;
 	}
       }
32234431
       if (precursor_idx >= 0) {
 	prec.spectrumID = spectrumId[precursor_idx];
5afbbd4f
       }
2edb9032
       // isolation window
3b091f3a
       if (!Rcpp::NumericVector::is_na(isolationWindowTargetMZ[i])) {
2edb9032
 	prec.isolationWindow.set(MS_isolation_window_target_m_z, isolationWindowTargetMZ[i]);
 	prec.isolationWindow.set(MS_isolation_window_lower_offset, isolationWindowLowerOffset[i]);
 	prec.isolationWindow.set(MS_isolation_window_upper_offset, isolationWindowUpperOffset[i]);
       }
e501fc3e
     }
     // [X] collisionEnergy
     // [ ] ionisationEnergy
     // [X] precursorScanNum
     // [X] precursorMZ
     // [X] precursorCharge
     // [X] precursorIntensity
     // [ ] mergedScan
45834cf0
 
e501fc3e
     Rcpp::NumericMatrix spct_vals = spctr_data[i];
     // mz values
     Rcpp::NumericVector mz_vals = spct_vals( Rcpp::_, 0);
     BinaryDataArrayPtr spct_mz(new BinaryDataArray);
     spct_mz->set(MS_m_z_array, "", MS_m_z);
     spct_mz->data.resize(mz_vals.size());
     for (int j = 0; j < mz_vals.size(); j++)
       spct_mz->data[j] = mz_vals[j];
     spct.binaryDataArrayPtrs.push_back(spct_mz);
     // intensity values
     Rcpp::NumericVector ints_vals = spct_vals( Rcpp::_, 1);
     BinaryDataArrayPtr spct_ints(new BinaryDataArray);
     spct_ints->set(MS_intensity_array, "", MS_number_of_detector_counts);
     spct_ints->data.resize(ints_vals.size());
     for (int j = 0; j < ints_vals.size(); j++)
       spct_ints->data[j] = ints_vals[j];
     spct.binaryDataArrayPtrs.push_back(spct_ints);
     spct.defaultArrayLength = spct_mz->data.size();
   }
 }
 
96ea799f
 Rcpp::DataFrame RcppPwiz::getChromatogramsInfo( int whichChrom )
84a4b166
 {
8a535503
   if (msd != NULL) {
     ChromatogramListPtr clp = msd->run.chromatogramListPtr;
     if (clp.get() == 0) {
f8520e7a
       Rf_warningcall(R_NilValue, "The direct support for chromatogram info is only available in mzML format.");
8a535503
       return Rcpp::DataFrame::create();
     } else if (clp->size() == 0) {
f8520e7a
       Rf_warningcall(R_NilValue, "No available chromatogram info.");
8a535503
       return Rcpp::DataFrame::create();
     } else if ( (whichChrom < 0) || (whichChrom > clp->size()) ) {
       Rprintf("Index whichChrom out of bounds [0 ... %d].\n", (clp->size())-1);
       return Rcpp::DataFrame::create( );
     } else {
       std::vector<double> time;
       std::vector<double> intensity;
       ChromatogramPtr c = clp->chromatogram(whichChrom, true);
       vector<TimeIntensityPair> pairs;
       c->getTimeIntensityPairs (pairs);
 
       for (int i = 0; i < pairs.size(); i++) {
 	TimeIntensityPair p = pairs.at(i);
 	time.push_back(p.time);
 	intensity.push_back(p.intensity);
       }
96ea799f
 
8a535503
       chromatogramsInfo = Rcpp::DataFrame::create(Rcpp::_["time"] = time,
 						  Rcpp::_[c->id]  = intensity);
96ea799f
 
84a4b166
     }
8a535503
     return(chromatogramsInfo);
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
8a535503
   return Rcpp::DataFrame::create( );
84a4b166
 }
 
3555a9aa
 // get the header info for chromatograms.
 Rcpp::DataFrame RcppPwiz::getChromatogramHeaderInfo (Rcpp::IntegerVector whichChrom)
 {
   if (msd != NULL) {
2aed7135
     // CVID nativeIdFormat_ = id::getDefaultNativeIDFormat(*msd);
3555a9aa
     ChromatogramListPtr clp = msd->run.chromatogramListPtr;
     if (clp.get() == 0) {
f8520e7a
       Rf_warningcall(R_NilValue, "The direct support for chromatogram info is only available in mzML format.");
3555a9aa
       return Rcpp::DataFrame::create();
     } else if (clp->size() == 0) {
f8520e7a
       Rf_warningcall(R_NilValue, "No available chromatogram info.");
3555a9aa
       return Rcpp::DataFrame::create();
     }
 
45834cf0
     int N = clp->size();
3555a9aa
     int N_chrom = whichChrom.size();
 
     Rcpp::StringVector chromatogramId(N_chrom); // the ID from the chrom
     Rcpp::IntegerVector chromatogramIndex(N_chrom);  // In contrast to the acquisitionNum we report here the index (1 based) of the chromatogram within the file.
     Rcpp::IntegerVector polarity(N_chrom);
     // MS:1000827: isolation window target m/z
     // MS:1000828: isolation window lower offset
     // MS:1000829: isolation window upper offset
     Rcpp::NumericVector precursorIsolationWindowTargetMZ(N_chrom);
     Rcpp::NumericVector precursorIsolationWindowLowerOffset(N_chrom);
     Rcpp::NumericVector precursorIsolationWindowUpperOffset(N_chrom);
     Rcpp::NumericVector precursorCollisionEnergy(N_chrom);
     Rcpp::NumericVector productIsolationWindowTargetMZ(N_chrom);
     Rcpp::NumericVector productIsolationWindowLowerOffset(N_chrom);
     Rcpp::NumericVector productIsolationWindowUpperOffset(N_chrom);
45834cf0
 
3555a9aa
     for (int i = 0; i < N_chrom; i++) {
       int current_chrom = whichChrom[i];
       if (current_chrom < 0 || current_chrom > N) {
f8520e7a
 	Rf_warningcall(R_NilValue, "Provided index out of bounds.");
3555a9aa
 	Rcpp::Rcerr << "Provided index out of bounds" << std::endl;
       }
       ChromatogramPtr ch = clp->chromatogram(current_chrom - 1, false);
       chromatogramId[i] = ch->id;
       chromatogramIndex[i] = current_chrom;
       CVParam param = ch->cvParamChild(MS_scan_polarity);
       polarity[i] = (param.cvid==MS_negative_scan ? 0 : (param.cvid==MS_positive_scan ? +1 : -1 ) );
       if (!ch->precursor.empty()) {
 	precursorIsolationWindowTargetMZ[i] = ch->precursor.isolationWindow.cvParam(MS_isolation_window_target_m_z).value.empty() ? NA_REAL : ch->precursor.isolationWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
a9633868
 	precursorIsolationWindowLowerOffset[i] = ch->precursor.isolationWindow.cvParam(MS_isolation_window_lower_offset).value.empty() ? NA_REAL : ch->precursor.isolationWindow.cvParam(MS_isolation_window_lower_offset).valueAs<double>();
 	precursorIsolationWindowUpperOffset[i] = ch->precursor.isolationWindow.cvParam(MS_isolation_window_upper_offset).value.empty() ? NA_REAL : ch->precursor.isolationWindow.cvParam(MS_isolation_window_upper_offset).valueAs<double>();
45834cf0
 	precursorCollisionEnergy[i] = ch->precursor.activation.cvParam(MS_collision_energy).value.empty() ? NA_REAL : ch->precursor.activation.cvParam(MS_collision_energy).valueAs<double>();
3555a9aa
       } else {
 	precursorIsolationWindowTargetMZ[i] = NA_REAL;
 	precursorIsolationWindowLowerOffset[i] = NA_REAL;
 	precursorIsolationWindowUpperOffset[i] = NA_REAL;
 	precursorCollisionEnergy[i] = NA_REAL;
       }
       if (!ch->product.empty()) {
 	productIsolationWindowTargetMZ[i] = ch->product.isolationWindow.cvParam(MS_isolation_window_target_m_z).value.empty() ? NA_REAL : ch->product.isolationWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
a9633868
 	productIsolationWindowLowerOffset[i] = ch->product.isolationWindow.cvParam(MS_isolation_window_lower_offset).value.empty() ? NA_REAL : ch->product.isolationWindow.cvParam(MS_isolation_window_lower_offset).valueAs<double>();
 	productIsolationWindowUpperOffset[i] = ch->product.isolationWindow.cvParam(MS_isolation_window_upper_offset).value.empty() ? NA_REAL : ch->product.isolationWindow.cvParam(MS_isolation_window_upper_offset).valueAs<double>();
3555a9aa
       } else {
 	productIsolationWindowTargetMZ[i] = NA_REAL;
 	productIsolationWindowLowerOffset[i] = NA_REAL;
 	productIsolationWindowUpperOffset[i] = NA_REAL;
       }
     }
     Rcpp::List chromHeader(10);
     std::vector<std::string> names;
     int i = 0;
     names.push_back("chromatogramId");
     chromHeader[i++] = Rcpp::wrap(chromatogramId);
     names.push_back("chromatogramIndex");
     chromHeader[i++] = Rcpp::wrap(chromatogramIndex);
     names.push_back("polarity");
     chromHeader[i++] = Rcpp::wrap(polarity);
     names.push_back("precursorIsolationWindowTargetMZ");
     chromHeader[i++] = Rcpp::wrap(precursorIsolationWindowTargetMZ);
     names.push_back("precursorIsolationWindowLowerOffset");
     chromHeader[i++] = Rcpp::wrap(precursorIsolationWindowLowerOffset);
     names.push_back("precursorIsolationWindowUpperOffset");
     chromHeader[i++] = Rcpp::wrap(precursorIsolationWindowUpperOffset);
     names.push_back("precursorCollisionEnergy");
     chromHeader[i++] = Rcpp::wrap(precursorCollisionEnergy);
     names.push_back("productIsolationWindowTargetMZ");
     chromHeader[i++] = Rcpp::wrap(productIsolationWindowTargetMZ);
     names.push_back("productIsolationWindowLowerOffset");
     chromHeader[i++] = Rcpp::wrap(productIsolationWindowLowerOffset);
     names.push_back("productIsolationWindowUpperOffset");
     chromHeader[i++] = Rcpp::wrap(productIsolationWindowUpperOffset);
45834cf0
 
3555a9aa
     chromHeader.attr("names") = names;
     return chromHeader;
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
3555a9aa
   return Rcpp::DataFrame::create( );
 }
 
 Rcpp::DataFrame RcppPwiz::getAllChromatogramHeaderInfo ( ) {
   if (msd != NULL) {
     ChromatogramListPtr clp = msd->run.chromatogramListPtr;
4b41c6ae
     if (clp.get() == 0) {
f8520e7a
       Rf_warningcall(R_NilValue, "The direct support for chromatogram info is only available in mzML format.");
4b41c6ae
       return Rcpp::DataFrame::create();
     }
3555a9aa
     int N = clp->size();
4b41c6ae
     if (N > 0) {
b61d82bf
       return getChromatogramHeaderInfo(Rcpp::seq(1, N));
4b41c6ae
     } else {
f8520e7a
       Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
4b41c6ae
     }
3555a9aa
   }
   return Rcpp::DataFrame::create( );
 }
 
84a4b166
 Rcpp::NumericMatrix RcppPwiz::get3DMap ( std::vector<int> scanNumbers, double whichMzLow, double whichMzHigh, double resMz )
 {
8a535503
   if (msd != NULL)
84a4b166
     {
 
8a535503
       SpectrumListPtr slp = msd->run.spectrumListPtr;
       double f = 1 / resMz;
       int low = round(whichMzLow * f);
       int high = round(whichMzHigh * f);
       int dmz = high - low + 1;
       int drt = scanNumbers.size();
84a4b166
 
8a535503
       Rcpp::NumericMatrix map3d(drt, dmz);
84a4b166
 
8a535503
       for (int i = 0; i < drt; i++)
84a4b166
         {
8a535503
 	  for (int j = 0; j < dmz; j++)
84a4b166
             {
8a535503
 	      map3d(i,j) = 0.0;
84a4b166
             }
         }
 
8a535503
       int j=0;
f8520e7a
       //Rprintf("%d\n",1);
8a535503
       for (int i = 0; i < scanNumbers.size(); i++)
84a4b166
         {
2aed7135
 	  SpectrumPtr s = slp->spectrum(scanNumbers[i] - 1, DetailLevel_FullMetadata);
8a535503
 	  vector<MZIntensityPair> pairs;
 	  s->getMZIntensityPairs(pairs);
84a4b166
 
8a535503
 	  for (int k=0; k < pairs.size(); k++)
84a4b166
             {
8a535503
 	      MZIntensityPair p = pairs.at(k);
 	      j = round(p.mz * f) - low;
 	      if ((j >= 0) & (j < dmz))
84a4b166
                 {
8a535503
 		  if (p.intensity > map3d(i,j))
84a4b166
                     {
8a535503
 		      map3d(i,j) = p.intensity;
84a4b166
                     }
                 }
             }
 
         }
8a535503
       return(map3d);
84a4b166
     }
 
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
8a535503
   return Rcpp::NumericMatrix(0,0);
84a4b166
 }
e501fc3e
 
e7c9b847
 string RcppPwiz::getRunStartTimeStamp() {
   if (msd != NULL) {
     return msd->run.startTimeStamp;
   }
f8520e7a
   Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
e7c9b847
   return "";
 }