C++ PeakSpectrum类代码示例

  void CompNovoIdentificationBase::getCIDSpectrumLight_(PeakSpectrum & spec, const String & sequence, DoubleReal prefix, DoubleReal suffix)
  {
    static DoubleReal h2o_mass = EmpiricalFormula("H2O").getMonoWeight();
    Peak1D p;
    DoubleReal b_pos(0.0 + prefix);
    DoubleReal y_pos(h2o_mass + suffix);
    for (Size i = 0; i != sequence.size() - 1; ++i)
    {
      char aa(sequence[i]);
      b_pos += aa_to_weight_[aa];

      char aa2(sequence[sequence.size() - i - 1]);
      y_pos += aa_to_weight_[aa2];

      if (b_pos > min_mz_ && b_pos < max_mz_)
      {
        p.setPosition(b_pos + Constants::PROTON_MASS_U);
        p.setIntensity(1.0f);
        spec.push_back(p);
      }

      if (y_pos > min_mz_ && y_pos < max_mz_)
      {
        p.setPosition(y_pos + Constants::PROTON_MASS_U);
        p.setIntensity(1.0f);
        spec.push_back(p);
      }
    }

    spec.sortByPosition();
    return;
  }

 std::vector<PeakSpectrum> AScore::peakPickingPerWindowsInSpectrum_(PeakSpectrum &real_spectrum) const
 {
   vector<PeakSpectrum> windows_top10;
   
   double spect_lower_bound = floor(real_spectrum.front().getMZ() / 100) * 100;
   double spect_upper_bound = ceil(real_spectrum.back().getMZ() / 100) * 100;
   
   Size number_of_windows = static_cast<Size>(ceil((spect_upper_bound - spect_lower_bound) / 100));
   windows_top10.resize(number_of_windows);
   
   PeakSpectrum::Iterator it_current_peak = real_spectrum.begin();
   Size window_upper_bound(spect_lower_bound + 100);
   
   for (Size current_window = 0; current_window < number_of_windows; ++current_window)
   {
     PeakSpectrum real_window;
     while (((*it_current_peak).getMZ() <= window_upper_bound) && (it_current_peak < real_spectrum.end()))
     {
       real_window.push_back(*it_current_peak);
       ++it_current_peak;
     }
     
     real_window.sortByIntensity(true);
     for (Size i = 0; (i < 10) & (i < real_window.size()); ++i)
     {
       windows_top10[current_window].push_back(real_window[i]);
     }
     
     window_upper_bound += 100;
   }
   return windows_top10;
 }

  double SpectrumAlignmentScore::operator()(const PeakSpectrum & s1, const PeakSpectrum & s2) const
  {
    const double tolerance = (double)param_.getValue("tolerance");
    bool is_relative_tolerance = param_.getValue("is_relative_tolerance").toBool();
    bool use_linear_factor = param_.getValue("use_linear_factor").toBool();
    bool use_gaussian_factor = param_.getValue("use_gaussian_factor").toBool();

    if (use_linear_factor && use_gaussian_factor)
    {
      cerr << "Warning: SpectrumAlignmentScore, use either 'use_linear_factor' or 'use_gaussian_factor'!" << endl;
    }

    SpectrumAlignment aligner;
    Param p;
    p.setValue("tolerance", tolerance);
    p.setValue("is_relative_tolerance", (String)param_.getValue("is_relative_tolerance"));
    aligner.setParameters(p);

    vector<pair<Size, Size> > alignment;
    aligner.getSpectrumAlignment(alignment, s1, s2);

    double score(0), sum(0), sum1(0), sum2(0);
    for (PeakSpectrum::ConstIterator it1 = s1.begin(); it1 != s1.end(); ++it1)
    {
      sum1 += it1->getIntensity() * it1->getIntensity();
    }

    for (PeakSpectrum::ConstIterator it1 = s2.begin(); it1 != s2.end(); ++it1)
    {
      sum2 += it1->getIntensity() * it1->getIntensity();
    }

    for (vector<pair<Size, Size> >::const_iterator it = alignment.begin(); it != alignment.end(); ++it)
    {
      //double factor(0.0);
      //factor = (epsilon - fabs(s1[it->first].getPosition()[0] - s2[it->second].getPosition()[0])) / epsilon;
      double mz_tolerance(tolerance);

      if (is_relative_tolerance)
      {
        mz_tolerance = mz_tolerance * s1[it->first].getPosition()[0] / 1e6;
      }

      double mz_difference(fabs(s1[it->first].getPosition()[0] - s2[it->second].getPosition()[0]));
      double factor = 1.0;

      if (use_linear_factor || use_gaussian_factor)
      {
        factor = getFactor_(mz_tolerance, mz_difference, use_gaussian_factor);
      }
      sum += sqrt(s1[it->first].getIntensity() * s2[it->second].getIntensity() * factor);
    }

    score = sum / (sqrt(sum1 * sum2));

    return score;
  }

  double XQuestScores::xCorrelationPrescore(const PeakSpectrum & spec1, const PeakSpectrum & spec2, double tolerance)
  {
    // return 0 = no correlation, when one of the spectra is empty
    if (spec1.size() == 0 || spec2.size() == 0) {
      return 0.0;
    }

    double maxionsize = std::max(spec1[spec1.size()-1].getMZ(), spec2[spec2.size()-1].getMZ());
    Int table_size = ceil(maxionsize / tolerance)+1;
    std::vector< double > ion_table1(table_size, 0);
    std::vector< double > ion_table2(table_size, 0);

    // Build tables of the same size, each bin has the size of the tolerance
    for (Size i = 0; i < spec1.size(); ++i)
    {
      Size pos = static_cast<Size>(ceil(spec1[i].getMZ() / tolerance));
      ion_table1[pos] = 1;
    }
    for (Size i = 0; i < spec2.size(); ++i)
    {
      Size pos =static_cast<Size>(ceil(spec2[i].getMZ() / tolerance));
      ion_table2[pos] = 1;

    }

    double dot_product = 0.0;
    for (Size i = 0; i < ion_table1.size(); ++i)
    {
      dot_product += ion_table1[i] * ion_table2[i];
    }

    // determine the smaller spectrum and normalize by the number of peaks in it
    double peaks = std::min(spec1.size(), spec2.size());
    return dot_product / peaks;
  }

 void SuffixArrayPeptideFinder::getCandidates(vector<vector<pair<SuffixArrayPeptideFinder::FASTAEntry, String> > >& candidates, const String& DTA_file)
 {
   DTAFile dta_file;
   PeakSpectrum s;
   dta_file.load(DTA_file, s);
   s.sortByPosition();
   PeakSpectrum::ConstIterator it(s.begin());
   vector<double> spec;
   for (; it != s.end(); ++it)
   {
     spec.push_back(it->getPosition()[0]);
   }
   const vector<double> specc(spec);
   getCandidates(candidates, specc);
   return;
 }

  map<Size, PeakSpectrum > PScore::calculatePeakLevelSpectra(const PeakSpectrum& spec, const vector<Size>& ranks, Size min_level, Size max_level)
  {
    map<Size, MSSpectrum<Peak1D> > peak_level_spectra;

    if (spec.empty()) return peak_level_spectra;

    // loop over all peaks and associated (zero-based) ranks
    for (Size i = 0; i != ranks.size(); ++i)
    {
      // start at the highest (less restrictive) level
      for (int j = static_cast<int>(max_level); j >= static_cast<int>(min_level); --j)
      {
        // if the current peak is annotated to have lower or equal rank then allowed for this peak level add it
        if (static_cast<int>(ranks[i]) <= j)
        {
          peak_level_spectra[j].push_back(spec[i]);
        }
        else
        {
          // if the current peak has higher rank than the current level then all it is also to high for the lower levels
          break;
        }
      }
    }
    return peak_level_spectra;
  }

  void CompNovoIdentificationBase::windowMower_(PeakSpectrum & spec, DoubleReal windowsize, Size no_peaks)
  {
    PeakSpectrum copy(spec);
    vector<Peak1D> to_be_deleted;
    for (Size i = 0; i < spec.size(); ++i)
    {
      PeakSpectrum sub_spec;
      bool end(false);
      for (Size j = i; spec[j].getPosition()[0] - spec[i].getPosition()[0] < windowsize; )
      {
        sub_spec.push_back(spec[j]);
        if (++j == spec.size())
        {
          end = true;
          break;
        }
      }

      sub_spec.sortByIntensity(true);

      for (Size k = no_peaks; k < sub_spec.size(); ++k)
      {
        Peak1D p(sub_spec[k]);
        to_be_deleted.push_back(p);
      }

      if (end)
      {
        break;
      }
    }

    spec.clear(false);
    for (PeakSpectrum::ConstIterator it = copy.begin(); it != copy.end(); ++it)
    {
      if (find(to_be_deleted.begin(), to_be_deleted.end(), *it) == to_be_deleted.end())
      {
        spec.push_back(*it);
      }
    }

    spec.sortByPosition();

  }

  double XQuestScores::logOccupancyProb(const PeakSpectrum& theoretical_spec,  const Size matched_size, double fragment_mass_tolerance, bool fragment_mass_tolerance_unit_ppm)
  {
    using boost::math::binomial;
    Size theo_size = theoretical_spec.size();

    if (matched_size < 1 || theo_size < 1)
    {
      return 0;
    }

    double range;
    double used_tolerance;

    if (fragment_mass_tolerance_unit_ppm)
    {
      range = std::log(theoretical_spec.back().getMZ()) - std::log(theoretical_spec[0].getMZ());
      used_tolerance = fragment_mass_tolerance / 1e6;
    }
    else
    {
      range = theoretical_spec.back().getMZ() - theoretical_spec[0].getMZ();
      used_tolerance = fragment_mass_tolerance;
    }

    // A priori probability of a random match given info about the theoretical spectrum
    double a_priori_p = 0;
    a_priori_p = 1 - pow(1 - 2 * used_tolerance / range,  static_cast<double>(theo_size));

    double log_occu_prob = 0;
    binomial flip(theo_size, a_priori_p);
    // min double number to avoid 0 values, causing scores with the value "inf"
    log_occu_prob = -log(1 - cdf(flip, matched_size) + std::numeric_limits<double>::min());

    // score lower than 0 does not make sense, but can happen, if cfd = 0, then -log( 1 + <double>::min() ) < 0
    if (log_occu_prob >= 0.0)
    {
      return log_occu_prob;
    }
    else // underflow warning?
    {
      return 0;
    }
  }

  std::vector< double > XQuestScores::xCorrelation(const PeakSpectrum & spec1, const PeakSpectrum & spec2, Int maxshift, double tolerance)
  {
    // generate vector of results, filled with zeroes
    std::vector< double > results(maxshift * 2 + 1, 0);

    // return 0 = no correlation, when one of the spectra is empty
    if (spec1.size() == 0 || spec2.size() == 0) {
      return results;
    }

    double maxionsize = std::max(spec1[spec1.size()-1].getMZ(), spec2[spec2.size()-1].getMZ());
    Int table_size = ceil(maxionsize / tolerance)+1;
    std::vector< double > ion_table1(table_size, 0);
    std::vector< double > ion_table2(table_size, 0);

    // Build tables of the same size, each bin has the size of the tolerance
    for (Size i = 0; i < spec1.size(); ++i)
    {
      Size pos = static_cast<Size>(ceil(spec1[i].getMZ() / tolerance));
      ion_table1[pos] = 10.0;
    }
    for (Size i = 0; i < spec2.size(); ++i)
    {
      Size pos =static_cast<Size>(ceil(spec2[i].getMZ() / tolerance));
      ion_table2[pos] = 10.0;
    }

    // Compute means
    double mean1 = (std::accumulate(ion_table1.begin(), ion_table1.end(), 0.0)) / table_size;
    double mean2 = (std::accumulate(ion_table2.begin(), ion_table2.end(), 0.0)) / table_size;

    // Compute denominator
    double s1 = 0;
    double s2 = 0;
    for (Int i = 0; i < table_size; ++i)
    {
      s1 += pow((ion_table1[i] - mean1), 2);
      s2 += pow((ion_table2[i] - mean2), 2);
    }
    double denom = sqrt(s1 * s2);

    // Calculate correlation for each shift
    for (Int shift = -maxshift; shift <= maxshift; ++shift)
    {
      double s = 0;
      for (Int i = 0; i < table_size; ++i)
      {
        Int j = i + shift;
        if ( (j >= 0) && (j < table_size))
        {
          s += (ion_table1[i] - mean1) * (ion_table2[j] - mean2);
        }
      }
      if (denom > 0)
      {
        results[shift + maxshift] = s / denom;
      }
    }
    return results;
  }

  /**
  @brief Similarity pairwise score

  This function return the similarity score of two spectra based on SteinScott.

  @param s1  const PeakSpectrum Spectrum 1
  @param s2  const PeakSpectrum Spectrum 2
  @see SteinScottImproveScore()
  */
  double SteinScottImproveScore::operator()(const PeakSpectrum & s1, const PeakSpectrum & s2) const
  {
    const double epsilon = (double)param_.getValue("tolerance");
    const double constant = epsilon / 10000;

    //const double c(0.0004);
    double score(0), sum(0), sum1(0), sum2(0), sum3(0), sum4(0);
    /* std::cout << s1 << std::endl;
    std::cout << std::endl;
    std::cout << s2 << std::endl;*/

    for (PeakSpectrum::ConstIterator it1 = s1.begin(); it1 != s1.end(); ++it1)
    {
      double temp = it1->getIntensity();
      sum1 += temp * temp;
      sum3 += temp;
    }

    for (PeakSpectrum::ConstIterator it1 = s2.begin(); it1 != s2.end(); ++it1)
    {
      double temp = it1->getIntensity();
      sum2 += temp * temp;
      sum4 += temp;
    }
    double z = constant * (sum3 * sum4);
    Size j_left(0);
    for (Size i = 0; i != s1.size(); ++i)
    {
      for (Size j = j_left; j != s2.size(); ++j)
      {
        double pos1(s1[i].getMZ()), pos2(s2[j].getMZ());
        if (std::abs(pos1 - pos2) <= 2 * epsilon)
        {
          sum += s1[i].getIntensity() * s2[j].getIntensity();
        }
        else
        {
          if (pos2 > pos1)
          {
            break;
          }
          else
          {
            j_left = j;
          }
        }
      }
    }
    //std::cout<< sum << " Sum " << z << " z " << std::endl;
    score = (sum - z) / (std::sqrt((sum1 * sum2)));
    // std::cout<<score<< " score" << std::endl;
    if (score < (float)param_.getValue("threshold"))
    {
      score = 0;
    }

    return score;
  }

  double XQuestScores::matchOddsScore(const PeakSpectrum& theoretical_spec,  const Size matched_size, double fragment_mass_tolerance, bool fragment_mass_tolerance_unit_ppm, bool is_xlink_spectrum, Size n_charges)
  {
    using boost::math::binomial;
    Size theo_size = theoretical_spec.size();

    if (matched_size < 1 || theo_size < 1)
    {
      return 0;
    }

    double range = theoretical_spec[theo_size-1].getMZ() -  theoretical_spec[0].getMZ();

    // Compute fragment tolerance in Da for the mean of MZ values, if tolerance in ppm (rough approximation)
    double mean = 0.0;
    for (Size i = 0; i < theo_size; ++i)
    {
      mean += theoretical_spec[i].getMZ();
    }
    mean = mean / theo_size;
    double tolerance_Th = fragment_mass_tolerance_unit_ppm ? mean * 1e-6 * fragment_mass_tolerance : fragment_mass_tolerance;

    // A priori probability of a random match given info about the theoretical spectrum
    double a_priori_p = 0;

    if (is_xlink_spectrum)
    {
      a_priori_p = (1 - ( pow( (1 - 2 * tolerance_Th / (0.5 * range)),  (static_cast<double>(theo_size) / static_cast<double>(n_charges)))));
    }
    else
    {
      a_priori_p = (1 - ( pow( (1 - 2 * tolerance_Th / (0.5 * range)),  static_cast<int>(theo_size))));
    }

    double match_odds = 0;

    binomial flip(theo_size, a_priori_p);
    // min double number to avoid 0 values, causing scores with the value "inf"
    match_odds = -log(1 - cdf(flip, matched_size) + std::numeric_limits<double>::min());

    // score lower than 0 does not make sense, but can happen if cfd = 0, -log( 1 + min() ) < 0
    if (match_odds >= 0.0)
    {
      return match_odds;
    }
    else
    {
      return 0;
    }
  }

  // s1 should be the original spectrum
  DoubleReal CompNovoIdentificationBase::compareSpectra_(const PeakSpectrum & s1, const PeakSpectrum & s2)
  {
    DoubleReal score(0.0);

    PeakSpectrum::ConstIterator it1 = s1.begin();
    PeakSpectrum::ConstIterator it2 = s2.begin();

    Size num_matches(0);
    while (it1 != s1.end() && it2 != s2.end())
    {
      DoubleReal pos1(it1->getPosition()[0]), pos2(it2->getPosition()[0]);
      if (fabs(pos1 - pos2) < fragment_mass_tolerance_)
      {
        score += it1->getIntensity();
        ++num_matches;
      }

      if (pos1 <= pos2)
      {
        ++it1;
      }
      else
      {
        ++it2;
      }
    }

    if (num_matches == 0)
    {
      return 0;
    }

    score /= sqrt((DoubleReal)num_matches);

    return score;
  }

 Size AScore::numberOfMatchedIons_(const PeakSpectrum & th, const PeakSpectrum & window, Size depth, double fragment_mass_tolerance, bool fragment_mass_tolerance_ppm) const
 {
   PeakSpectrum window_reduced = window;
   if (window_reduced.size() > depth)
   {
     window_reduced.resize(depth);
   }
   
   window_reduced.sortByPosition();
   Size n = 0;
   for (Size i = 0; i < th.size(); ++i)
   {
     Size nearest_peak = -1;
     try
     {
       nearest_peak = window_reduced.findNearest(th[i].getMZ());
     }
     catch (Exception::Precondition) {}
     
     if (nearest_peak < window_reduced.size())
     {
       double window_mz = window_reduced[nearest_peak].getMZ();
       double error = abs(window_mz - th[i].getMZ());
       
       if (fragment_mass_tolerance_ppm)
       {
         error = error / window_mz * 1e6;
       }
       if (error < fragment_mass_tolerance)
       {
         ++n;
       }
     }      
   }
   return n;
 }

  void CompNovoIonScoring::scoreETDFeatures_(Size /*charge*/, double precursor_weight, Map<double, IonScore> & ion_scores, const PeakSpectrum & CID_spec, const PeakSpectrum & ETD_spec)
  {
    //double fragment_mass_tolerance((double)param_.getValue("fragment_mass_tolerance"));
    Size max_isotope_to_score(param_.getValue("max_isotope_to_score"));

    for (PeakSpectrum::ConstIterator it1 = CID_spec.begin(); it1 != CID_spec.end(); ++it1)
    {
      double pos1(it1->getPosition()[0]);
      double b_sum(0.0), y_sum(0.0);

      // score a-ions
      for (PeakSpectrum::ConstIterator it2 = CID_spec.begin(); it2 != CID_spec.end(); ++it2)
      {
        double pos2(it2->getPosition()[0]);
        if (fabs(pos1 - pos2 - 28.0) < fragment_mass_tolerance_)
        {
          double factor((fragment_mass_tolerance_ - fabs(pos1 - pos2 - 28.0)) / fragment_mass_tolerance_);
#ifdef SCORE_ETDFEATURES_DEBUG
          cerr << "scoreETDFeatures: found a-ion " << pos1 << " (" << pos2 << ") (factor=" << factor << ") " << b_sum << " -> ";
#endif
          b_sum += it2->getIntensity() * factor;
#ifdef SCORE_ETDFEATURES_DEBUG
          cerr << endl;
#endif
        }
      }

      for (PeakSpectrum::ConstIterator it2 = ETD_spec.begin(); it2 != ETD_spec.end(); ++it2)
      {
        double pos2(it2->getPosition()[0]);

        // check if pos2 is precursor doubly charged, which has not fragmented
        double pre_diff_lower = (precursor_weight + Constants::PROTON_MASS_U) / 2.0 - fragment_mass_tolerance_;
        double pre_diff_upper = (precursor_weight + 4.0 * Constants::PROTON_MASS_U) / 2.0 + fragment_mass_tolerance_;
        if (pos2 > pre_diff_lower && pos2 < pre_diff_upper)
        {
#ifdef SCORE_ETDFEATURES_DEBUG
          cerr << "scoreETDFeatures: pre-range: " << pos2 << " is in precursor peak range: " << pre_diff_lower << " <-> " << pre_diff_upper << endl;
#endif
          continue;
        }

        //double diff(pos2 - pos1);

        // pos1 is CID ion; pos2 is ETD ion
        // pos1 b-ion, pos2 c-ion
        if (fabs(pos1 + 17.0 - pos2) < fragment_mass_tolerance_)
        {
          // now test if the ETD peak has "isotope" pattern
          double factor((fragment_mass_tolerance_ - fabs(pos1 + 17.0 - pos2)) / fragment_mass_tolerance_);
#ifdef SCORE_ETDFEATURES_DEBUG
          cerr << "scoreETDFeatures: is b-ion: " << pos1 << " (" << pos2 << ") (factor=" << factor << ") " << b_sum << " -> ";
#endif
          vector<double> iso_pattern;
          iso_pattern.push_back(it1->getIntensity());
          double actual_pos = it1->getPosition()[0];
          for (PeakSpectrum::ConstIterator it3 = it2; it3 != ETD_spec.end(); ++it3)
          {
            double it3_pos(it3->getPosition()[0]);
            if (fabs(fabs(actual_pos - it3_pos) - Constants::NEUTRON_MASS_U) < fragment_mass_tolerance_)
            {
              iso_pattern.push_back(it3->getIntensity());
              actual_pos = it3_pos;
            }
            if (iso_pattern.size() == max_isotope_to_score)
            {
              break;
            }
          }

          if (ion_scores[it1->getPosition()[0]].is_isotope_1_mono != -1)
          {
            b_sum += it2->getIntensity() * iso_pattern.size() * factor;
          }
#ifdef SCORE_ETDFEATURES_DEBUG
          cerr << b_sum << endl;
#endif
        }



        // pos1 z-ion, pos2 y-ion
        if (fabs(pos2 + 16.0 - pos1) < fragment_mass_tolerance_)
        {
          double factor((fragment_mass_tolerance_ - fabs(pos2 + 16.0 - pos1)) / fragment_mass_tolerance_);
          // now test if the ETD peak has "isotope" pattern
#ifdef SCORE_ETDFEATURES_DEBUG
          cerr << "scoreETDFeatures: is y-ion: " << pos1 << " (" << pos2 << ") (factor=" << factor << ") " << y_sum << " -> ";
#endif
          vector<double> iso_pattern;
          iso_pattern.push_back(it1->getIntensity());
          double actual_pos = it1->getPosition()[0];
          for (PeakSpectrum::ConstIterator it3 = it2; it3 != ETD_spec.end(); ++it3)
          {
            double it3_pos(it3->getPosition()[0]);
            if (fabs(fabs(actual_pos - it3_pos) - Constants::NEUTRON_MASS_U) < fragment_mass_tolerance_)
            {
              iso_pattern.push_back(it3->getIntensity());
              actual_pos = it3_pos;
            }
//.........这里部分代码省略.........

  void CompNovoIonScoring::scoreWitnessSet_(Size charge, double precursor_weight, Map<double, IonScore> & ion_scores, const PeakSpectrum & CID_spec)
  {
    vector<double> diffs;
    //diffs.push_back(28.0);
    diffs.push_back(17.0);
    diffs.push_back(18.0);

    // witnesses of CID spec (diffs)
    for (PeakSpectrum::ConstIterator it1 = CID_spec.begin(); it1 != CID_spec.end(); ++it1)
    {
      //Size num_wit(0);
      double wit_score(0.0);
      double pos1(it1->getPosition()[0]);
      wit_score += it1->getIntensity();
      for (PeakSpectrum::ConstIterator it2 = CID_spec.begin(); it2 != CID_spec.end(); ++it2)
      {
        double pos2(it2->getPosition()[0]);

        // direct ++
        if (charge > 1)
        {
          if (fabs(pos2 * 2 - Constants::PROTON_MASS_U - pos1) < fragment_mass_tolerance_)
          {
            double factor((fragment_mass_tolerance_ - fabs(pos2 * 2 - Constants::PROTON_MASS_U - pos1)) / fragment_mass_tolerance_);
            // pos1 is ion, pos2 is ++ion
#ifdef SCORE_WITNESSSET_DEBUG
            cerr << "scoreWitnessSet: ++ion " << pos1 << " " << pos2 << " (factor=" << factor << ") " << wit_score << " -> ";
#endif
            if (ion_scores[it2->getPosition()[0]].s_isotope_pattern_2 < 0.2)
            {
              wit_score += it2->getIntensity() * /* 0.2 */ factor;
            }
            else
            {
              wit_score += it2->getIntensity() * ion_scores[it2->getPosition()[0]].s_isotope_pattern_2 * factor;
            }
#ifdef SCORE_WITNESSSET_DEBUG
            cerr << wit_score << endl;
#endif
          }
        }

        // diffs?
        for (vector<double>::const_iterator it = diffs.begin(); it != diffs.end(); ++it)
        {
          // pos1 is ion, pos2 loss peak
          if (fabs(pos1 - pos2 - *it) < fragment_mass_tolerance_)
          {
            double factor((fragment_mass_tolerance_ - fabs(pos1 - pos2 - *it)) / fragment_mass_tolerance_);
#ifdef SCORE_WITNESSSET_DEBUG
            cerr << "scoreWitnessSet: diff " << pos1 << " (" << pos2 << ") " << *it << " (factor=" << factor << ") " << wit_score << " -> ";
#endif
            wit_score += it2->getIntensity() /* / 5.0*/ * factor;
#ifdef SCORE_WITNESSSET_DEBUG
            cerr << wit_score << endl;
#endif
          }
        }

        // is there a b-ion?; pos1 is ion, pos2 complementary ion
        if (fabs(pos1 + pos2 - 1 * Constants::PROTON_MASS_U - precursor_weight) < fragment_mass_tolerance_)
        {
          double factor((fragment_mass_tolerance_ - fabs(pos1 + pos2 - Constants::PROTON_MASS_U - precursor_weight)) / fragment_mass_tolerance_);
          /*factor *= 0.2;*/
#ifdef SCORE_WITNESSSET_DEBUG
          cerr << "scoreWitnessSet: complementary " << pos1 << " (" << pos2 << ") (factor=" << factor << ") " << wit_score << " -> ";
#endif
          // found complementary ion
          if (ion_scores[it2->getPosition()[0]].s_isotope_pattern_1 < 0.5 || ion_scores[it2->getPosition()[0]].is_isotope_1_mono != 1)
          {
            wit_score += it2->getIntensity() /* * 0.5*/ * factor;
          }
          else
          {
            wit_score += it2->getIntensity() * ion_scores[it2->getPosition()[0]].s_isotope_pattern_1 * factor;
          }
#ifdef SCORE_WITNESSSET_DEBUG
          cerr << wit_score << endl;
#endif

          if (ion_scores[it2->getPosition()[0]].s_bion != 0)
          {
#ifdef SCORE_WITNESSSET_DEBUG
            cerr << "scoreWitnessSet: complementary is b-ion " << pos1 << "(" << pos2 << ")" << wit_score << " -> ";
#endif
            wit_score += ion_scores[it2->getPosition()[0]].s_bion * factor;
#ifdef SCORE_WITNESSSET_DEBUG
            cerr << wit_score << endl;
#endif
          }

        }
      }

      // isotope pattern ok?
      if (ion_scores[it1->getPosition()[0]].s_isotope_pattern_1 > 0 && ion_scores[it1->getPosition()[0]].is_isotope_1_mono == 1)
      {
#ifdef SCORE_WITNESSSET_DEBUG
        cerr << "scoreWitnessSet: isotope pattern: " << pos1 << " " << wit_score << " -> ";
#endif
//.........这里部分代码省略.........