C++ Spectrum::getScanNumber方法代码示例

//#ifdef _MSC_VER 
//#ifndef MSSINGLESCAN_MAIN
//int mssinglescan_main(int argc, char * argv[] ) {
//#else
//int main(int argc, char *argv[]){
//#endif
//#else
int main(int argc, char *argv[]){
//#endif


	//Here are all the variable we are going to need
	MSReader r;
	Spectrum s;
  int j;

  if(argc==1){
    printf("DESCRIPTION: Reads an MS/MS spectrum from any MSToolkit supported file type and outputs to screen in MS2 format.\n\n");
    printf("USAGE: MSSingleScan [scan number] [file]\n");
    exit(0);
  };

	r.readFile(argv[2],s,atoi(argv[1]));
  if(s.getScanNumber()==0) exit(-1);
  printf("S\t%d\t%d\t%.*f\n",s.getScanNumber(),s.getScanNumber(),2,s.getMZ());
	if(s.getRTime()>0) printf("I\tRTime\t%.*f\n",4,s.getRTime());
	for(j=0;j<s.sizeZ();j++){
		printf("Z\t%d\t%.*f\n",s.atZ(j).z,2,s.atZ(j).mz);
	};

	for(j=0;j<s.size();j++){
		printf("%.4f %.4f\n",s.at(j).mz,s.at(j).intensity);
	};

  return 0;

};

// assumes at least one spectrum in allRefs
void SearchLibrary::runSearch(Spectrum& s)
{
    scoreMatches(s, cachedSpectra_, targetMatches_);
    scoreMatches(s, cachedDecoySpectra_, decoyMatches_);

    // keep scores from all target psms for estimating Weibull parameters
    vector<double> allScores;
    if(compute_pvalues_){
        for(size_t i=0; i < targetMatches_.size(); i++){
            double dotp = targetMatches_[i].getScore(DOTP);
            allScores.push_back(dotp);
        }
    }

    // there may have been spectra in cachedSpectra_ but none at the
    // correct charge state.  Check again
    if( targetMatches_.size() == 0 ){
        Verbosity::warn("No library spectra found for query %d "
                        "(precursor m/z %.2f).", s.getScanNumber(), s.getMz());
        return;
    }
    if( compute_pvalues_ ){
        addNullScores(s, allScores);
    }

    // sort the matches descending
    sort(targetMatches_.begin(), targetMatches_.end(), compMatchDotScore); 
    sort(decoyMatches_.begin(), decoyMatches_.end(), compMatchDotScore); 

    setRank();
    
    if( printAll_ ){
        cout << "spec " << s.getScanNumber() << endl;
    }

    if( compute_pvalues_ ){
        weibullEstimator_.estimateParams(allScores);
        
        // print params to file
        if( weibullParamFile_.is_open() ){
            weibullParamFile_ << s.getScanNumber() << "\t"
                              << weibullEstimator_.getEta() << "\t"
                              << weibullEstimator_.getBeta() << "\t"
                              << weibullEstimator_.getShift() << "\t"
                              << weibullEstimator_.getCorrelation() << "\t"
                              << weibullEstimator_.getNumPointsFit() 
                //(int)(allScores.size() * fraction_to_fit_)
                              << endl;
        }
        setMatchesPvalues(allScores.size());
    }
}

/**
 *  Generate more scores by creating decoy spectra and comparing them
 *  to query.  Create one decoy for each target until there are the
 *  minimum number of scores.  Do not save decoy spectrum or its
 *  Match.  Makes no changes to cachedSpectra_.
 */
void SearchLibrary::addNullScores(Spectrum s, vector<double>& allScores){
    int shiftAmount = 5;

    while((int)allScores.size() < minWeibullScores_) {
        int specAdded = 0; // make sure spectra were successfully added
        
        //loop through all candidate refs, create shifted spectrum, compare
        for(size_t i=0; i < targetMatches_.size(); i++) {
            const RefSpectrum* targetSpec = targetMatches_.at(i).getRefSpec();
            RefSpectrum* decoySpec = targetSpec->newDecoy(shiftAmount,
                                                          shiftRawSpectra_);
            if( decoySpec == NULL ){
                continue;
            }
            specAdded++;
            Match thisMatch(&s, decoySpec); // spectrum deleted with match?
            thisMatch.setMatchLibID(0);
            Verbosity::comment(V_ALL, 
                      "Comparing query spec %d and shifted library spec %d",
                      s.getScanNumber(), decoySpec->getLibSpecID() );
            
            DotProduct::compare(thisMatch);
            allScores.push_back(thisMatch.getScore(DOTP));
            
        } // next ref spectrum

        shiftAmount += 5;

        if( specAdded == 0 ){
            break;  // avoid infinite loop if number of scores not increasing
        }
    } // next pass through all ref spectra
}

/**
 * Compare the given query spectrum to all library spectra.  Create a
 * match for each and add to matches.
 */
void SearchLibrary::scoreMatches(Spectrum& s, deque<RefSpectrum*>& spectra,
                                 vector<Match>& matches ){
    Verbosity::debug("Scoring %d matches.", spectra.size());
    // get the charge states we will search
    const vector<int>& charges = s.getPossibleCharges();
    
    // compare all ref spec to query, create match for each
    for(size_t i=0; i< spectra.size(); i++) {
        // is there a better place to check this?
        if(spectra.at(i)->getNumProcessedPeaks() == 0 ){ 
            Verbosity::debug("Skipping library spectrum %d.  No peaks.", 
                             spectra.at(i)->getLibSpecID());
            continue;
        }
        
        if( ! checkCharge(charges, spectra.at(i)->getCharge()) ){
            continue;
        }
        
        Match thisMatch(&s, spectra.at(i));  
        
        thisMatch.setMatchLibID(spectra.at(i)->getLibID());
        
        Verbosity::comment(V_ALL, "Comparing query spec %d and library spec %d",
                           s.getScanNumber(), spectra.at(i)->getLibSpecID());
        
        DotProduct::compare(thisMatch); //static method
        
        // save match for reporting
        matches.push_back(thisMatch);
    }  
}

bool CPeptideDatabase::readMascot(char* mascotFile, char* dataFile){
  MSReader r;
  MascotParser p;
  Spectrum s;
  MascotLite m;

  int i;
  int percent=0;

  p.clear();
  p.readText(mascotFile);

  if(p.isDistiller()){
    cout << mascotFile << "\t" << p.size() << " peptides from Distiller, no need for RT lookup." << endl;
  } else if(p.hasRTime()){
    cout << mascotFile << "\t" << p.size() << " peptides have RT already, no need for lookup." << endl;
  } else {
    cout << mascotFile << "\tLooking up retention times for " << p.size() << " Mascot peptides..." << endl;
    cout << percent;
    r.readFile(dataFile,s);
  }

  for(i=0;i<p.size();i++){
    if(p.isDistiller()||p.hasRTime()){ 
      m.rTime=p[i].rTime;
    } else {
      if( (int)i*100/p.size() > percent){
        percent=(int)i*100/p.size();
        cout << "\b\b\b" << percent;
      }

      r.readFile(NULL,s,p[i].scanNum);
      if(s.getScanNumber()==0){
        cout << "Scan not found. Please check that correct data file is used. Stopping Mascot parsing." << endl;
        return false;
      }
      m.rTime=s.getRTime();
    }
    m.charge=p[i].charge;
    m.fileID=fileID;
    m.monoMass=p[i].zeroMass;
    m.scanNum=p[i].scanNum;
    strcpy(m.peptide,p[i].sequence_Long);
    strcpy(m.protein,p[i].gene);
    mascot.push_back(m);
  }
  cout << "  Done!" << endl;

  fileID++;

  return true;
}

bool CNoiseReduction::DeNoiseD(Spectrum& sp){

  Spectrum tmpSpec;  
  vector<int> v;

  sp.clear();

  if(pos==0){
    if(cs.scan.iLower>0) {
      if(!NewScanAverage(sp,cs.inFile,(int)(cs.boxcar/2),0.1f,cs.scan.iLower)) return false;
    } else {
      NewScanAverage(sp,cs.inFile,(int)(cs.boxcar/2),0.1f);
    }
    pos=1;
  } else {
    NewScanAverage(sp,NULL,(int)(cs.boxcar/2),0.1f);
  }

  if(sp.getScanNumber()==0) return false;
  FirstDerivativePeaks(sp,1);
  return true;

}

bool CNoiseReduction::NewScanAverage(Spectrum& sp, char* file, int width, float cutoff, int scanNum){
  
  Spectrum ts;

  vector<int> vPos;
   
  int i;
  int j;
  int k;
  int m;
  int widthCount=0;
  int numScans=1;
  double dif;
  double prec;
  double dt;
  double c;

  bool bLeft=true;
  int posLeft;
  int posRight;
  int index;
  char cFilter1[256];
  //char cFilter2[256];

  //double slope;
  //double intercept;

  sp.clear();
  Spectrum* specs;
  specs = new Spectrum[width*2+1];

  //if file is not null, create new buffer
  if(file!=NULL){
    strcpy(lastFile,file);
    bs.clear();
    if(scanNum>0) r->readFile(file,ts,scanNum);
    else r->readFile(file,ts);
    if(ts.getScanNumber()==0) {
      delete [] specs;
      return false;
    }
    bs.push_back(ts);
    specs[0]=bs[0];
    c=CParam(specs[0],3);
    posA=0;
  } else {
    posA++;
    if(posA>=(int)bs.size()) { //end of buffer, no more data
      delete [] specs;
      return false; 
    }
    specs[0]=bs[posA];
    c=CParam(specs[0],3);
  }

  specs[0].getRawFilter(cFilter1,256);

  posLeft=posA;
  posRight=posA;
  while(widthCount<(width*2)){

    index=-1;

    //Alternate looking left and right
    if(bLeft){
      bLeft=false;
      widthCount++;
      
      while(true){
        posLeft--;
        if(posLeft<0) { //buffer is too short on left, add spectra
          i=bs[0].getScanNumber();
          while(true){
            i--;
            if(i==0) break;
            r->readFile(lastFile,ts,i);
            if(ts.getScanNumber()==0) continue;
            else break;
          }
          if(i==0) break;
          bs.push_front(ts);
          posA++;
          posRight++;
          posLeft=0;
          //ts.getRawFilter(cFilter2,256);
          if(ts.getMsLevel()==cs.msLevel) {
            index=posLeft;
            break;
          }
        } else {
          //bs[posLeft].getRawFilter(cFilter2,256);
          if(bs[posLeft].getMsLevel()==cs.msLevel) {
            index=posLeft;
            break;
          }
        }
      }

    } else {
      bLeft=true;
//.........这里部分代码省略.........

bool CNoiseReduction::ScanAverage(Spectrum& sp, char* file, int width, float cutoff){
  
  Spectrum ts;
  Spectrum ps=sp;
  MSReader r;
 
  int i;
  int j;
  int k;
  int widthCount=0;
  int numScans=1;
  double dif;
  double prec;
  double dt;
  double c=CParam(ps,3);

  bool bLeft=true;
  int posLeft=ps.getScanNumber()-1;
  int posRight=ps.getScanNumber()+1;
  char cFilter1[256];
  char cFilter2[256];

  ps.getRawFilter(cFilter1,256);

  while(widthCount<(width*2)){

    //Alternate looking left and right
    if(bLeft){
      bLeft=false;
      widthCount++;
      while(posLeft>0){
        r.readFile(file,ts,posLeft);
        if(ts.getScanNumber()==0) break;
        ts.getRawFilter(cFilter2,256);
        if(strcmp(cFilter1,cFilter2)==0) break;
        posLeft--;
      }
    } else {
      bLeft=true;
      widthCount++;
      while(true){
        r.readFile(file,ts,posRight);
        if(ts.getScanNumber()==0) break;
        ts.getRawFilter(cFilter2,256);
        if(strcmp(cFilter1,cFilter2)==0) break;
        posRight++;
      }
    }
    if(ts.getScanNumber()==0) continue;

    numScans++;

    //Match peaks between pivot scan and temp scan
    k=0;
    for(i=0;i<ps.size();i++){
      dif=100000.0;
      prec = c * ps.at(i).mz * ps.at(i).mz / 2;

      for(j=k;j<ts.size();j++){
        dt=fabs(ps.at(i).mz-ts.at(j).mz);
        if(dt<=dif) {
          if(dt<prec) {
            ps.at(i).intensity+=ts.at(j).intensity;
            ts.at(j).mz=-1.0;
            k=j+1;
            break;
          }
          dif=dt;
        } else {
          k=j-1;
          break;
        }
      }
    }

    //Add unmatched peaks from temp scan
    for(i=0;i<ts.size();i++){
      if(ts.at(i).mz>-1.0) ps.add(ts.at(i));
    }

    //Sort pivot scan peaks for fast traversal against next temp scan
    ps.sortMZ();

  }

  //Average points and apply cutoff
  sp.clear();
  for(i=0;i<ps.size();i++) {
    ps.at(i).intensity/=numScans;
    if(ps.at(i).intensity>=cutoff) sp.add(ps.at(i));
  }

  sp.setScanNumber(ps.getScanNumber());
  sp.setScanNumber(ps.getScanNumber(true),true);
  sp.setRTime(ps.getRTime());

  return true;
}

bool CNoiseReduction::DeNoise(Spectrum& sp){

  double ppm;
  int i,j,k;
  int index;
  int matchCount;
  char cFilter1[256];
  char cFilter2[256];

  Spectrum tmpSpec;  
  vector<int> v;

  sp.clear();

  if(pos==0){
    if((cs.scan.iLower>0)) {
      k=cs.scan.iLower; 
      r->readFile(&cs.inFile[0],tmpSpec,cs.scan.iLower);
    } else {
      r->readFile(&cs.inFile[0],tmpSpec);
      k=tmpSpec.getScanNumber();
    }
    if(tmpSpec.getScanNumber()==0) return false;
    tmpSpec.getRawFilter(cFilter1,256);
      
    //Gather left side of scan
    i=1;
    j=0;
    while( k-i > 0){
      r->readFile(&cs.inFile[0],tmpSpec,cs.scan.iLower-i);
      if(tmpSpec.getScanNumber()==0) {
        i++;
        continue;
      }
      tmpSpec.getRawFilter(cFilter2,256);
      
      //Assume High resolution data at all times
      if(!cs.centroid) {
        FirstDerivativePeaks(tmpSpec,1);
        tmpSpec.setRawFilter(cFilter2);
      }

      s.push_front(tmpSpec);
      if(strcmp(cFilter1,cFilter2)==0){
        j++;
				if(j==(int)(cs.boxcar/2)) break;
      }
      i++;
    }

    //cout << "Done left " << s.size() << " " << cs.rawAvgWidth << endl;

    //Get our position again
    r->readFile(&cs.inFile[0],tmpSpec,k);

    //Assume High resolution data at all times
    if(!cs.centroid) {
      FirstDerivativePeaks(tmpSpec,1);
      tmpSpec.setRawFilter(cFilter1);
    }
    
    //Add first target scan
    s.push_back(tmpSpec);
    pos=s.size()-1; 
    
    //Add right side of scan
    i=1;
    j=0;
    while(true){
      r->readFile(NULL,tmpSpec);
      if(tmpSpec.getScanNumber()==0) break;
      tmpSpec.getRawFilter(cFilter2,256);

      //Assume High resolution data at all times
      if(!cs.centroid) {
        FirstDerivativePeaks(tmpSpec,1);
        tmpSpec.setRawFilter(cFilter2);
      }

      s.push_back(tmpSpec);
      //cout << tmpSpec.getScanNumber() << " " << cFilter1 << " xx " << cFilter2 << endl;
      if(strcmp(cFilter1,cFilter2)==0){  
        j++;
				if(j==(int)(cs.boxcar/2)) break;
      }
      i++;
    }

  }

  //return false if we reached the end
  if(pos>=(int)s.size()) return false;

  //Because Deque may have mixed spectra, create list of Deque indexes to compare
  //Shrink Deque on the left if some spectra are not needed
  //Expand Deque to the right if needed and possible
  //cout << "Checking scan " << s[pos].getScanNumber() << endl;
  v.clear();
  s[pos].getRawFilter(cFilter1,256);

//.........这里部分代码省略.........

//First derivative method taken from CSpecAnalyze, returns base peak intensity of the set
void CNoiseReduction::FirstDerivativePeaks(Spectrum& sp, int winSize){
  int i,j;
  float maxIntensity;
  int bestPeak;
  bool bLastPos;
  Spectrum gp;

	int nextBest;
	double FWHM;
	Peak_T centroid;

  bLastPos=false;
  for(i=0;i<sp.size()-winSize;i++){

    if(sp.at(i).intensity<sp.at(i+winSize).intensity) {
      bLastPos=true;
      continue;
    } else {
      if(bLastPos){
				bLastPos=false;
	
				//find max and add peak
				maxIntensity=0;
				for(j=i;j<i+winSize;j++){
				  if (sp.at(j).intensity>maxIntensity){
				    maxIntensity=sp.at(j).intensity;
				    bestPeak = j;
				  }
				}

				//Best estimate of Gaussian centroid
				//Get 2nd highest point of peak
				if(bestPeak==sp.size()-1){
					nextBest=bestPeak-1;
				} else if(sp.at(bestPeak-1).intensity > sp.at(bestPeak+1).intensity){
					nextBest=bestPeak-1;
				} else {
					nextBest=bestPeak+1;
				}

				//Get FWHM
				FWHM = calcFWHM(sp.at(bestPeak).mz);

				//Calc centroid MZ (in three lines for easy reading)
				centroid.mz = pow(FWHM,2)*log(sp.at(bestPeak).intensity/sp.at(nextBest).intensity);
				centroid.mz /= GC*(sp.at(bestPeak).mz-sp.at(nextBest).mz);
				centroid.mz += (sp.at(bestPeak).mz+sp.at(nextBest).mz)/2;

				//Calc centroid intensity
				centroid.intensity=(float)(sp.at(bestPeak).intensity/exp(-pow((sp.at(bestPeak).mz-centroid.mz)/FWHM,2)*GC));

				//some peaks are funny shaped and have bad gaussian fit.
				//if error is more than 10%, keep existing intensity
				if( fabs((sp.at(bestPeak).intensity - centroid.intensity) / centroid.intensity * 100) > 10 ||
            //not a good check for infinity
            centroid.intensity>999999999999.9 ||
            centroid.intensity < 0 ) {
					centroid.intensity=sp.at(bestPeak).intensity;
				}

				//Hack until I put in mass ranges
				if(centroid.mz<0 || centroid.mz>2000) {
					//do nothing if invalid mz
				} else {
					gp.add(centroid);
				}
				i+=winSize-1;
      }

    }
  }
  
  int scanNumber=sp.getScanNumber();
  int scanNumber2=sp.getScanNumber(false);
  float rTime=sp.getRTime();
  sp = gp;
  sp.setRTime(rTime);
  sp.setScanNumber(scanNumber);
  sp.setScanNumber(scanNumber2,true);

}

bool CNoiseReduction::NewScanAveragePlusDeNoise(Spectrum& sp, char* file, int width, float cutoff, int scanNum){
  
  Spectrum ts;

  vector<int> vPos;
 
  int i;
  int j;
  int k;
  int m;
  int widthCount=0;
  int numScans=1;
  int match;
  double dif;
  double prec;
  double dt;
  double c;

  bool bLeft=true;
  int posLeft;
  int posRight;
  int index;
  char cFilter1[256];
  char cFilter2[256];

  sp.clear();

  Spectrum* specs;
  specs = new Spectrum[width*2+1];

  //if file is not null, create new buffer
  if(file!=NULL){
    strcpy(lastFile,file);
    bs.clear();
    if(scanNum>0) r->readFile(file,ts,scanNum);
    else r->readFile(file,ts);
    if(ts.getScanNumber()==0) {
      delete [] specs;
      return false;
    }
    bs.push_back(ts);
    specs[0]=bs[0];
    c=CParam(specs[0],3);
    posA=0;
  } else {
    posA++;
    if(posA>=(int)bs.size()) { //end of buffer, no more data
      delete [] specs;
      return false; 
    }
    specs[0]=bs[posA];
    c=CParam(specs[0],3);
  }

  //set our pivot spectrum
  specs[0].getRawFilter(cFilter1,256);

  posLeft=posA;
  posRight=posA;
  while(widthCount<(width*2)){

    index=-1;

    //Alternate looking left and right
    if(bLeft){
      bLeft=false;
      widthCount++;
      
      while(true){
        posLeft--;
        if(posLeft<0) { //buffer is too short on left, add spectra
          i=bs[0].getScanNumber();
          while(true){
            i--;
            if(i==0) break;
            r->readFile(lastFile,ts,i);
            if(ts.getScanNumber()==0) continue;
            else break;
          }
          if(i==0) break;
          bs.push_front(ts);
          posA++;
          posRight++;
          posLeft=0;
          ts.getRawFilter(cFilter2,256);
          if(strcmp(cFilter1,cFilter2)==0) {
            index=posLeft;
            break;
          }
        } else {
          bs[posLeft].getRawFilter(cFilter2,256);
          if(strcmp(cFilter1,cFilter2)==0) {
            index=posLeft;
            break;
          }
        }
      }

    } else {
      bLeft=true;
//.........这里部分代码省略.........

bool CNoiseReduction::ScanAveragePlusDeNoise(Spectrum& sp, char* file, int width, float cutoff, int scanNum){
  
  Spectrum ts;
  Spectrum ps;
  //MSReader r;

  vector<int> v;
  vector<int> vPos;
 
  int i;
  int j;
  int k;
  int widthCount=0;
  int numScans=1;
  int match;
  double dif;
  double prec;
  double dt;
  double c;

  bool bLeft=true;
  int posLeft=ps.getScanNumber()-1;
  int posRight=ps.getScanNumber()+1;
  int index;
  char cFilter1[256];
  //char cFilter2[256];

  sp.clear();

  //if file is not null, create new buffer
  if(file!=NULL){
    strcpy(lastFile,file);
    bs.clear();
    if(scanNum>0) r->readFile(file,ts,scanNum);
    else r->readFile(file,ts);
    if(ts.getScanNumber()==0) return false;
    bs.push_back(ts);
    ps=bs[0];
    c=CParam(ps,3);
    posA=0;
  } else {
    posA++;
    //cout << "ER: " << posA << " " << bs.size() << endl;
    if(posA>=(int)bs.size()) return false; //end of buffer, no more data
    ps=bs[posA];
    c=CParam(ps,3);
  }

  //set our pivot spectrum
  //ps=bs[posA];
  ps.getRawFilter(cFilter1,256);
  //cout << "Averaging: " << ps.getScanNumber() << endl;

  posLeft=posA;
  posRight=posA;
  while(widthCount<(width*2)){

    index=-1;

    //Alternate looking left and right
    if(bLeft){
      bLeft=false;
      widthCount++;
      
      while(true){
        posLeft--;
        //cout << posLeft << endl;
        if(posLeft<0) { //buffer is too short on left, add spectra
          i=bs[0].getScanNumber();
          while(true){
            i--;
            //cout << "I: " << i << endl;
            if(i==0) break;
            r->readFile(lastFile,ts,i);
            if(ts.getScanNumber()==0) continue;
            else break;
          }
          if(i==0) break;
          bs.push_front(ts);
          for(i=0;i<(int)v.size();i++)v[i]++;
          posA++;
          posRight++;
          posLeft=0;
          //ts.getRawFilter(cFilter2,256);
          if(ts.getMsLevel()==cs.msLevel) {
            index=posLeft;
            break;
          }
        } else {
          //bs[posLeft].getRawFilter(cFilter2,256);
          if(bs[posLeft].getMsLevel()==cs.msLevel) {
            index=posLeft;
            break;
          }
        }
      }

    } else {
      bLeft=true;
      widthCount++;
//.........这里部分代码省略.........

// Function reads in the spectra from the data stored for each PSMClass object
void PepXMLClass::readInSpectra() {
	string curSpectrumFilePath;

	deque<PSMClass>::iterator curPSM;
        SpecStruct *spec = NULL;
        bool status;
        string spectrumFileName;
        int scanNum;
        int ctr = 0;
        int N = (signed) PSMvec->size();
        
        if(g_ext == "mgf") {
            parseMGF();
            return;
        }
        
        // mstoolkit to read MS2 spectra
        MSReader *reader = new MSReader();
        reader->setFilter(MS2);

       // mstoolkit spectrum object
       Spectrum *S = NULL;

	// Extract from each PSM, it's parent spectrum file and scan number
	for(curPSM = PSMvec->begin(); curPSM != PSMvec->end(); curPSM++) {
            spectrumFileName = curPSM->getSpectrumFileName();
            scanNum = curPSM->getScanNumber();
            
            
            // code to get the full path to the input spectrum file
            filesystem::path curFile( spectrumFileName.c_str() );
            filesystem::path spectral_dir( g_srcDir.c_str() );
            filesystem::path curFilePath( spectral_dir/curFile );
            curSpectrumFilePath = curFilePath.file_string();

            // if the source file is not located, drop the spectrum
            if( !boost::filesystem::exists(curSpectrumFilePath) ) {
                curPSM = PSMvec->erase(curPSM);
                continue;
            }
            
            
            // Read the spectrum for 'scanNum' into 'S'
            S = new Spectrum();
            status = reader->readFile(curSpectrumFilePath.c_str(), *S, scanNum);
            if(S->getScanNumber() == 0) {
                cerr << "Failed to get " << curSpectrumFilePath.c_str() << " scan#: " << scanNum << endl;
                exit(0);
            }
            
            spec = new SpecStruct();
            for(int j = 0; j < S->size(); j++) {
                spec->mz.push_back( S->at(j).mz );
                spec->intensity.push_back( S->at(j).intensity );
            }
            delete(S); S = NULL;
            
            // Assign the data in 'S' to curPSM
            curPSM->recordSpectrum(*spec);
            delete(spec); spec = NULL;
            ctr++;
            
            printProgress("Reading in spectra (please be patient)...", ctr, N);
	}
       delete(reader); reader = NULL;
       cerr << endl; // prettier stderr
}