#include "RcppPwiz.h"

RcppPwiz::RcppPwiz()
{
  msd = NULL;
  nativeIdFormat = CVID_Unknown;
  instrumentInfo = Rcpp::List::create();
  chromatogramsInfo = Rcpp::DataFrame::create();
  isInCacheInstrumentInfo = FALSE;
  allScanHeaderInfo = Rcpp::List::create();
  isInCacheAllScanHeaderInfo = FALSE;
}
/* Destructor*/
RcppPwiz::~RcppPwiz()
{
  RcppPwiz::close();
}

// void RcppPwiz::open(const string& fileName)
void RcppPwiz::open(Rcpp::StringVector fileName)
{

  filename = Rcpp::as<std::string>(fileName(0));
  msd = new MSDataFile(filename);
  nativeIdFormat = id::getDefaultNativeIDFormat(*msd);
}

/* Release all memory on close. */
void RcppPwiz::close()
{
  if (msd != NULL)
    {
      delete msd;
      msd = NULL;
      nativeIdFormat = CVID_Unknown;
      instrumentInfo = Rcpp::List::create();
      chromatogramsInfo = Rcpp::DataFrame::create();
      isInCacheInstrumentInfo = FALSE;
      allScanHeaderInfo = Rcpp::List::create();
      isInCacheAllScanHeaderInfo = FALSE;
    }
}

string RcppPwiz::getFilename() {
  return filename;
}

int RcppPwiz::getLastScan() const {
  if (msd != NULL) {
    SpectrumListPtr slp = msd->run.spectrumListPtr;
    return slp->size();
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return -1;
}

int RcppPwiz::getLastChrom() const {
  if (msd != NULL) {
    ChromatogramListPtr clp = msd->run.chromatogramListPtr;
    return clp->size();
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return -1;
}

Rcpp::List RcppPwiz::getInstrumentInfo ( )
{
  if (msd != NULL)
    {
      if (!isInCacheInstrumentInfo)
        {

	  vector<InstrumentConfigurationPtr> icp = msd->instrumentConfigurationPtrs; // NULL for mzData
	  if (icp.size() != 0)
            {
	      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 {
		ionisation = std::string(adapter.ionisation());
	      } catch(...) {}
	      try {
		analyzer = std::string(adapter.analyzer());
	      } catch(...) {}
	      try {
		detector = std::string(adapter.detector());
	      } 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:"")
						  ) ;

            }
	  else
            {
	      instrumentInfo = Rcpp::List::create(
						  Rcpp::_["manufacturer"]  = "",
						  Rcpp::_["model"]         = "",
						  Rcpp::_["ionisation"]    = "",
						  Rcpp::_["analyzer"]      = "",
						  Rcpp::_["detector"]      = "",
						  Rcpp::_["software"]      = "",
						  Rcpp::_["sample"]		  = "",
						  Rcpp::_["source"]		  = ""
						  ) ;
            }

	  isInCacheInstrumentInfo = TRUE;
        }
      return(instrumentInfo);
    }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return instrumentInfo;
}

int RcppPwiz::getAcquisitionNumber(string id, size_t index) const
{
  // 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);
}

// 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;
// }

Rcpp::DataFrame RcppPwiz::getScanHeaderInfo (Rcpp::IntegerVector whichScan) {
  if (msd != NULL) {
    SpectrumListPtr slp = msd->run.spectrumListPtr;
    size_t N = slp->size();
    size_t N_scans = whichScan.size();
    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);

    for (size_t i = 0; i < N_scans; i++) {
      int current_scan = whichScan[i];
      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);
      Scan dummy;
      Scan& scan = sp->scanList.scans.empty() ? dummy : sp->scanList.scans[0];
      if (scan.empty())
	Rprintf("Scan with index %d empty.\n", current_scan);
      // seqNum
      seqNum[i] = current_scan;
      acquisitionNum[i] = getAcquisitionNumber(sp->id, current_index);
      // spectrumId
      spectrumId[i] = Rcpp::String(sp->id);
      // 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);
      centroided[i] = (param.cvid==MS_centroid_spectrum ? TRUE : (param.cvid==MS_profile_spectrum ? FALSE : NA_LOGICAL));
      // 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);
	if (precursorIndex < N) {
	  precursorScanNum[i] = getAcquisitionNumber(precursor.spectrumID, precursorIndex);
	} else {
	  precursorScanNum[i] = NA_INTEGER;
	}
	// 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;
      }
    }

    Rcpp::List header(31);
    std::vector<std::string> names;
    size_t i = 0;
    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");
    header[i++] = Rcpp::wrap(ionMobilityDriftTime);
    names.push_back("isolationWindowTargetMZ");
    header[i++] = Rcpp::wrap(isolationWindowTargetMZ);
    names.push_back("isolationWindowLowerOffset");
    header[i++] = Rcpp::wrap(isolationWindowLowerOffset);
    names.push_back("isolationWindowUpperOffset");
    header[i++] = Rcpp::wrap(isolationWindowUpperOffset);
    names.push_back("scanWindowLowerLimit");
    header[i++] = Rcpp::wrap(scanWindowLowerLimit);
    names.push_back("scanWindowUpperLimit");
    header[i++] = Rcpp::wrap(scanWindowUpperLimit);
    header.attr("names") = names;

    return header;
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return Rcpp::DataFrame::create( );
}

Rcpp::DataFrame RcppPwiz::getAllScanHeaderInfo ( ) {
  if (msd != NULL) {
    if (!isInCacheAllScanHeaderInfo) {
      SpectrumListPtr slp = msd->run.spectrumListPtr;
      size_t N = slp->size();

      allScanHeaderInfo = getScanHeaderInfo(Rcpp::seq(1, N));
      isInCacheAllScanHeaderInfo = TRUE;
    }
    return allScanHeaderInfo ;
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return Rcpp::DataFrame::create( );
}

Rcpp::List RcppPwiz::getPeakList(Rcpp::IntegerVector whichScan) {
  if (msd != NULL) {
    SpectrumListPtr slp = msd->run.spectrumListPtr;
    size_t n_scans = slp->size();
    size_t n_want = whichScan.size();
    int current_scan;
    SpectrumPtr sp;
    BinaryDataArrayPtr mzs,ints;
    std::vector<double> data;
    Rcpp::NumericVector data_matrix;
    Rcpp::List res(n_want);
    for (size_t i = 0; i < n_want; i++) {
      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( );
      }
      size_t current_index = static_cast<size_t>(current_scan - 1);
      // sp = slp->spectrum(current_index, true);
      sp = slp->spectrum(current_index, DetailLevel_FullData);
      mzs = sp->getMZArray();
      ints = sp->getIntensityArray();
      if (!mzs.get() || !ints.get()) {
	Rcpp::NumericMatrix pks(0, 2);
	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);
      res[i] = data_matrix;
    }
    return res;
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return Rcpp::List::create();
}

/**
 * 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.
 * software_processing: has to be a list of character vectors.
 **/
void RcppPwiz::copyWriteMSfile(const string& file, const string& format,
			       const string& originalFile,
			       Rcpp::DataFrame spctr_header,
			       Rcpp::List spctr_data,
			       bool rtime_seconds,
			       Rcpp::List software_processing) {
  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
  //   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;
  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);

  // o paramGroupList
  if (format != "mzxml")
    newmsd.paramGroupPtrs = msd->paramGroupPtrs;
  // o sampleList
  newmsd.samplePtrs = msd->samplePtrs;
  // o instrumentConfigurationList
  newmsd.instrumentConfigurationPtrs = msd->instrumentConfigurationPtrs;
  // o softwareList
  newmsd.softwarePtrs = msd->softwarePtrs;
  // o dataProcessingList
  newmsd.dataProcessingPtrs = msd->dataProcessingPtrs;
  // 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)));
    }
  }

  // o run
  // Initialize the run and fill with data from the original file.
  Run &original_run = msd->run;
  newmsd.run.id = original_run.id;
  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;
  }
  // Now filling with new data
  addSpectrumList(newmsd, spctr_header, spctr_data, rtime_seconds);

  if (format == "mgf") {
    std::ofstream* mgfOutFileP = new std::ofstream(file.c_str());
    Serializer_MGF serializerMGF;
    serializerMGF.write(*mgfOutFileP, newmsd);
    mgfOutFileP->flush();
    mgfOutFileP->close();
  } else if (format == "mzxml") {
    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();
  } else if (format == "mzml") {
    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();
  }
  else
    Rcpp::Rcerr << format << " format not supported! Please try mgf, mzML, mzXML or mz5." << std::endl;

  // 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,
				 bool rtime_seconds,
				 Rcpp::List software_processing) {
  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);

  // 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)));
    }
  }

  newmsd.run.id = "Experiment_1";

  // Now filling with new data
  addSpectrumList(newmsd, spctr_header, spctr_data, rtime_seconds);

  if (format == "mgf") {
    std::ofstream* mgfOutFileP = new std::ofstream(file.c_str());
    Serializer_MGF serializerMGF;
    serializerMGF.write(*mgfOutFileP, newmsd);
    mgfOutFileP->flush();
    mgfOutFileP->close();
  } else if (format == "mzxml") {
    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();
  } else if (format == "mzml") {
    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();
  }
  else
    Rcpp::Rcerr << format << " format not supported! Please try mgf, mzML, mzXML or mz5." << std::endl;
}

/*
 * 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).
 *   soft_proc[2]: The CV ID of the software. Use "MS:-1" if not known, in
 *                 which case we are NOT writing the corresponding CV element.
 *   soft_proc[3-length]: CV IDs of the processing steps (optional).
 */
void RcppPwiz::addDataProcessing(MSData& msd, Rcpp::StringVector soft_proc) {
  SoftwarePtr new_soft(new Software);
  new_soft->id = soft_proc(0);
  new_soft->version = soft_proc(1);
  int soft_proc_size = soft_proc.size();
  if (soft_proc_size > 2) {
    if (soft_proc(2) != "MS:-1") {
      CVTermInfo cv_term = cvTermInfo(soft_proc(2));
      new_soft->set(cv_term.cvid);
    }
  }
  // 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);
    }
  }
  data_processing->processingMethods.push_back(proc_meth);
  msd.softwarePtrs.push_back(new_soft);
  msd.dataProcessingPtrs.push_back(data_processing);
}

/** 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].
 **/
void RcppPwiz::addSpectrumList(MSData& msd,
			       Rcpp::DataFrame& spctr_header,
			       Rcpp::List& spctr_data,
			       bool rtime_seconds) {
  int precursor_idx;
  // 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"];
  Rcpp::NumericVector totIonCurrent = spctr_header["totIonCurrent"];
  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"];
  Rcpp::IntegerVector mergedScan = spctr_header["mergedScan"];
  // Skipping mergedResultScanNum, mergedResultStartScanNum and mergedResultEndScanNum
  Rcpp::NumericVector ionInjectionTime = spctr_header["injectionTime"];
  Rcpp::StringVector filterString = spctr_header["filterString"];
  Rcpp::StringVector spectrumId = spctr_header["spectrumId"];
  Rcpp::LogicalVector centroided = spctr_header["centroided"];
  Rcpp::NumericVector ionMobilityDriftTime = spctr_header["ionMobilityDriftTime"];
  Rcpp::NumericVector isolationWindowTargetMZ = spctr_header["isolationWindowTargetMZ"];
  Rcpp::NumericVector isolationWindowLowerOffset = spctr_header["isolationWindowLowerOffset"];
  Rcpp::NumericVector isolationWindowUpperOffset = spctr_header["isolationWindowUpperOffset"];
  Rcpp::NumericVector scanWindowLowerLimit = spctr_header["scanWindowLowerLimit"];
  Rcpp::NumericVector scanWindowUpperLimit = spctr_header["scanWindowUpperLimit"];

  // 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

  // Now filling with new data
  shared_ptr<SpectrumListSimple> spectrumList(new SpectrumListSimple);
  // Add the default Processing pointer (fix issue #151
  spectrumList->dp = msd.dataProcessingPtrs[(msd.dataProcessingPtrs.size() - 1)];
  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]);
    // centroided
    if (centroided[i] != NA_LOGICAL && centroided[i] == TRUE)
      spct.set(MS_centroid_spectrum);
    if (centroided[i] != NA_LOGICAL && centroided[i] == FALSE)
      spct.set(MS_profile_spectrum);
    // [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);
    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);
    spct.set(MS_total_ion_current, totIonCurrent[i]);
    // TODO:
    // [X] seqNum: number observed in file.
    spct.index = seqNum[i] - 1;	// Or just i?
    // [X] acquisitionNum: number as reported (there might be gaps).
    // spct.id = "scan=" + boost::lexical_cast<std::string>(acquisitionNum[i]);
    spct.id = spectrumId[i];	// Use the provided ID instead
    // [ ] 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();
    if (rtime_seconds)
      spct_scan.set(MS_scan_start_time, retentionTime[i], UO_second);
    else
      spct_scan.set(MS_scan_start_time, retentionTime[i], UO_minute);
    if (ionInjectionTime[i] > 0)
      spct_scan.set(MS_ion_injection_time, ionInjectionTime[i], UO_millisecond);

    if (!Rcpp::StringVector::is_na(filterString[i]))
      spct_scan.set(MS_filter_string, filterString[i]);

    if (ionMobilityDriftTime[i] != NA_REAL)
      spct_scan.set(MS_ion_mobility_drift_time, ionMobilityDriftTime[i],
		    UO_millisecond);

    // scanWindow
    if (scanWindowLowerLimit[i] != NA_REAL && scanWindowUpperLimit[i] != NA_REAL) {
      spct_scan.scanWindows.push_back(ScanWindow(scanWindowLowerLimit[i], scanWindowUpperLimit[i], MS_m_z));
    }
    // MSn - precursor:
    if ( (precursorScanNum[i] > 0) | (precursorMZ[i] > 0) ) {
      // 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]);
      // Get the spectrumId of the precursor. Assuming that precursorScanNum is
      // linked to the acquisitionNum of the precursor.
      // This seems to be correct, since both the acquisitionNum and the
      // precursorNum are extracted from the respective spectrum's ID.
      precursor_idx = -1;
      for (int j = 0; j < spctr_data.size(); j++) {
	if (precursorScanNum[i] == acquisitionNum[j]) {
	  precursor_idx = j;
	  break;
	}
      }
      if (precursor_idx >= 0) {
	prec.spectrumID = spectrumId[precursor_idx];
      }
      // isolation window
      if (isolationWindowTargetMZ[i] != NA_REAL) {
	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]);
      }
    }
    // [X] collisionEnergy
    // [ ] ionisationEnergy
    // [X] precursorScanNum
    // [X] precursorMZ
    // [X] precursorCharge
    // [X] precursorIntensity
    // [ ] mergedScan

    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();
  }
}

Rcpp::DataFrame RcppPwiz::getChromatogramsInfo( int whichChrom )
{
  if (msd != NULL) {
    ChromatogramListPtr clp = msd->run.chromatogramListPtr;
    if (clp.get() == 0) {
      Rf_warningcall(R_NilValue, "The direct support for chromatogram info is only available in mzML format.");
      return Rcpp::DataFrame::create();
    } else if (clp->size() == 0) {
      Rf_warningcall(R_NilValue, "No available chromatogram info.");
      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);
      }

      chromatogramsInfo = Rcpp::DataFrame::create(Rcpp::_["time"] = time,
						  Rcpp::_[c->id]  = intensity);

    }
    return(chromatogramsInfo);
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return Rcpp::DataFrame::create( );
}

// get the header info for chromatograms.
Rcpp::DataFrame RcppPwiz::getChromatogramHeaderInfo (Rcpp::IntegerVector whichChrom)
{
  if (msd != NULL) {
    // CVID nativeIdFormat_ = id::getDefaultNativeIDFormat(*msd);
    ChromatogramListPtr clp = msd->run.chromatogramListPtr;
    if (clp.get() == 0) {
      Rf_warningcall(R_NilValue, "The direct support for chromatogram info is only available in mzML format.");
      return Rcpp::DataFrame::create();
    } else if (clp->size() == 0) {
      Rf_warningcall(R_NilValue, "No available chromatogram info.");
      return Rcpp::DataFrame::create();
    }

    int N = clp->size();
    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);

    for (int i = 0; i < N_chrom; i++) {
      int current_chrom = whichChrom[i];
      if (current_chrom < 0 || current_chrom > N) {
	Rf_warningcall(R_NilValue, "Provided index out of bounds.");
	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>();
	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>();
	precursorCollisionEnergy[i] = ch->precursor.activation.cvParam(MS_collision_energy).value.empty() ? NA_REAL : ch->precursor.activation.cvParam(MS_collision_energy).valueAs<double>();
      } 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>();
	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>();
      } 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);

    chromHeader.attr("names") = names;
    return chromHeader;
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return Rcpp::DataFrame::create( );
}

Rcpp::DataFrame RcppPwiz::getAllChromatogramHeaderInfo ( ) {
  if (msd != NULL) {
    ChromatogramListPtr clp = msd->run.chromatogramListPtr;
    if (clp.get() == 0) {
      Rf_warningcall(R_NilValue, "The direct support for chromatogram info is only available in mzML format.");
      return Rcpp::DataFrame::create();
    }
    int N = clp->size();
    if (N > 0) {
      return getChromatogramHeaderInfo(Rcpp::seq(1, N));
    } else {
      Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
    }
  }
  return Rcpp::DataFrame::create( );
}

Rcpp::NumericMatrix RcppPwiz::get3DMap ( std::vector<int> scanNumbers, double whichMzLow, double whichMzHigh, double resMz )
{
  if (msd != NULL)
    {

      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();

      Rcpp::NumericMatrix map3d(drt, dmz);

      for (int i = 0; i < drt; i++)
        {
	  for (int j = 0; j < dmz; j++)
            {
	      map3d(i,j) = 0.0;
            }
        }

      int j=0;
      //Rprintf("%d\n",1);
      for (int i = 0; i < scanNumbers.size(); i++)
        {
	  SpectrumPtr s = slp->spectrum(scanNumbers[i] - 1, DetailLevel_FullMetadata);
	  vector<MZIntensityPair> pairs;
	  s->getMZIntensityPairs(pairs);

	  for (int k=0; k < pairs.size(); k++)
            {
	      MZIntensityPair p = pairs.at(k);
	      j = round(p.mz * f) - low;
	      if ((j >= 0) & (j < dmz))
                {
		  if (p.intensity > map3d(i,j))
                    {
		      map3d(i,j) = p.intensity;
                    }
                }
            }

        }
      return(map3d);
    }

  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return Rcpp::NumericMatrix(0,0);
}

string RcppPwiz::getRunStartTimeStamp() {
  if (msd != NULL) {
    return msd->run.startTimeStamp;
  }
  Rf_warningcall(R_NilValue, "pwiz not yet initialized.");
  return "";
}