C++ SpectrumPtr::getIntensityArray方法代码示例

PWIZ_API_DECL
msdata::SpectrumPtr SpectrumList_MZWindow::spectrum(size_t index, bool getBinaryData) const
{
    SpectrumPtr spectrum = inner_->spectrum(index, getBinaryData);
    if (!getBinaryData) return spectrum;

    vector<MZIntensityPair> data;
    spectrum->getMZIntensityPairs(data);

    vector<MZIntensityPair>::const_iterator begin = lower_bound(
        data.begin(), data.end(), MZIntensityPair(mzLow_,0), hasLowerMZ);

    vector<MZIntensityPair>::const_iterator end = upper_bound(
        data.begin(), data.end(), MZIntensityPair(mzHigh_,0), hasLowerMZ);

    vector<MZIntensityPair> newData;
    copy(begin, end, back_inserter(newData));

    BinaryDataArrayPtr intensityArray = spectrum->getIntensityArray();
    CVID intensityUnits = intensityArray.get() ? intensityArray->cvParam(MS_intensity_unit).cvid : CVID_Unknown;
    
    SpectrumPtr newSpectrum(new Spectrum(*spectrum));
    newSpectrum->binaryDataArrayPtrs.clear();
    newSpectrum->setMZIntensityPairs(newData, intensityUnits);

    return newSpectrum;
}

    void writeCompressedPeaks(SpectrumPtr s, ostream& os)
    {
        // Build arrays to hold peaks prior to compression
        int numPeaks = (int) s->defaultArrayLength;
        double *pD = new double[numPeaks];
        float *pF = new float[numPeaks];

        const BinaryDataArray& mzArray = *s->getMZArray();
        const BinaryDataArray& intensityArray = *s->getIntensityArray();
        for(int j = 0; j < numPeaks; j++)
        {
            pD[j] = mzArray.data[j];
            pF[j] = (float) intensityArray.data[j];
        }

        // compress mz
        uLong sizeM = (uLong) (numPeaks * sizeDoubleMSn);
        uLong comprLenM = compressBound(sizeM);
        Byte *comprM = (Byte*)calloc((uInt)comprLenM, 1);
        int retM = compress(comprM, &comprLenM, (const Bytef*)pD, sizeM);
        
        // compress intensity
        uLong sizeI = (uLong) (numPeaks * sizeFloatMSn);
        uLong comprLenI = compressBound(sizeI);
        Byte *comprI = (Byte*)calloc((uInt)comprLenI, 1);
        int retI = compress(comprI, &comprLenI, (const Bytef*)pF, sizeI);

        // Write the compressed peaks if all is well
        if ((Z_OK == retM) && (Z_OK == retI))
        {
            // write length of compressed array of m/z
            os.write(reinterpret_cast<char *>(&comprLenM), sizeIntMSn);

            // write length of compressed array of intensities
            os.write(reinterpret_cast<char *>(&comprLenI), sizeIntMSn);

            // write compressed array of m/z
            os.write(reinterpret_cast<char *>(comprM), comprLenM);

            // write compressed array of intensities
            os.write(reinterpret_cast<char *>(comprI), comprLenI);
        }

        // Clean up memory
        free(comprM);
        free(comprI);
        delete [] pD;
        delete [] pF;
        
        // In case of error, throw exception AFTER cleaning up memory
        if (Z_OK != retM || Z_OK != retI)
        {
            throw runtime_error("[Serializer_MSn::writeCompressedPeaks] Error compressing peaks.");
        }
    }

    void writeSpectrumText(SpectrumPtr s, ostream& os)
    {
        os << std::setprecision(7); // 124.4567
        
        // Write the scan numbers 
        os << "S\t";
        int scanNum = getScanNumber(s);
        os << scanNum <<  "\t" << scanNum << "\t";

        // Write the precursor mz
        Precursor& precur = s->precursors[0];
        SelectedIon& si = precur.selectedIons[0];
        double mz = si.cvParam(MS_selected_ion_m_z).valueAs<double>();
        os << mz << "\n";
        
        // Write the scan time, if available
        if( s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds() )
          os << "I\tRTime\t" << s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds() << "\n";

        // Collect charge and mass info
        vector<int> charges;
        vector<double> masses;
        int numChargeStates = 0;
        // for each selected ion
        BOOST_FOREACH(const SelectedIon& curIon, precur.selectedIons){
          numChargeStates += getChargeStates(curIon, charges, masses);
        }

        // Write EZ lines if accurate masses are available
        CVParam massParam = si.cvParam(MS_accurate_mass);
        if( !massParam.empty() ){
          for(int i=0; i < numChargeStates; i++){
            os << "I\tEZ\t" << charges[i] << "\t" << masses[i] << "\t0\t0" << endl; // pad last two fields with 0
          }
        }

        // For each charge, write the charge and mass
        for(int i = 0; i < numChargeStates; i++)
        {
            os << "Z\t" << charges[i] << "\t" << masses[i] << "\n"; 
        }

        
        // Write each mz, intensity pair
        const BinaryDataArray& mzArray = *s->getMZArray();
        const BinaryDataArray& intensityArray = *s->getIntensityArray();
        for (size_t p=0; p < s->defaultArrayLength; ++p)
        {
            os << mzArray.data[p] << " " << intensityArray.data[p] << "\n";
        }
    }

PWIZ_API_DECL SpectrumPtr SpectrumList_MetadataFixer::spectrum(size_t index, bool getBinaryData) const
{
    // always get binary data
    SpectrumPtr s = inner_->spectrum(index, true);

    BinaryDataArrayPtr mzArray = s->getMZArray();
    BinaryDataArrayPtr intensityArray = s->getIntensityArray();
    if (!mzArray.get() || !intensityArray.get())
        return s;

    vector<double>& mzs = mzArray->data;
    vector<double>& intensities = intensityArray->data;

    double tic = 0;
    if (!mzs.empty())
    {
        double bpmz, bpi = -1;
        for (size_t i=0, end=mzs.size(); i < end; ++i)
        {
            tic += intensities[i];
            if (bpi < intensities[i])
            {
                bpi = intensities[i];
                bpmz = mzs[i];
            }
        }

        replaceCvParam(*s, MS_base_peak_intensity, bpi, MS_number_of_counts);
        replaceCvParam(*s, MS_base_peak_m_z, bpmz, MS_m_z);
        replaceCvParam(*s, MS_lowest_observed_m_z, mzs.front(), MS_m_z);
        replaceCvParam(*s, MS_highest_observed_m_z, mzs.back(), MS_m_z);
    }

    replaceCvParam(*s, MS_TIC, tic, MS_number_of_counts);

    return s;
}

    void writeSpectrumBinary(SpectrumPtr s, int version, bool compress, ostream& os)
    {
      // todo
      // MSnScanInfo header(s);
      // header.writeSpectrumHeader(version, compress, os);

        int scanNum = getScanNumber(s);
        os.write(reinterpret_cast<char *>(&scanNum), sizeIntMSn);
        os.write(reinterpret_cast<char *>(&scanNum), sizeIntMSn); // Yes, there are two

        Precursor& precur = s->precursors[0];
        SelectedIon& si = precur.selectedIons[0];
        double mz = si.cvParam(MS_selected_ion_m_z).valueAs<double>();
        os.write(reinterpret_cast<char *>(&mz), sizeDoubleMSn);

        float rt = (float) s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds();
        os.write(reinterpret_cast<char *>(&rt), sizeFloatMSn);
        
        if (version >= 2)
        {
            float basePeakIntensity = s->cvParam(MS_base_peak_intensity).valueAs<float>();
            os.write(reinterpret_cast<char *>(&basePeakIntensity), sizeFloatMSn);

            double basePeakMZ = s->cvParam(MS_base_peak_m_z).valueAs<double>();
            os.write(reinterpret_cast<char *>(&basePeakMZ), sizeDoubleMSn);
            
            // We don't have this information, but we need to write something,
            // so pad with 0's. (version 2 specific data)
            double conversionFactorA = (double)0;
            os.write(reinterpret_cast<char *>(&conversionFactorA), sizeDoubleMSn); 
            double conversionFactorB = (double)0;
            os.write(reinterpret_cast<char *>(&conversionFactorB), sizeDoubleMSn); 

            double tic = s->cvParam(MS_total_ion_current).valueAs<double>();
            os.write(reinterpret_cast<char *>(&tic), sizeDoubleMSn);

            // TODO
            float ionInjectionTime = (float)0;
            os.write(reinterpret_cast<char *>(&ionInjectionTime), sizeFloatMSn);
        }
        
        vector<int> charges;
        vector<double> masses;
        int numChargeStates = 0;
        BOOST_FOREACH(const SelectedIon& curIon, precur.selectedIons)
        {
            numChargeStates += getChargeStates(curIon, charges, masses);
        }
        os.write(reinterpret_cast<char *>(&numChargeStates), sizeIntMSn);
        
        bool hasAccurateMass = false;
        if (version == 3)
        {
          int numEzStates = 0;
          CVParam massParam = si.cvParam(MS_accurate_mass);
          if (!massParam.empty())
          {
            numEzStates = numChargeStates;
            hasAccurateMass = true;
          }
          os.write(reinterpret_cast<char *>(&numEzStates), sizeIntMSn);
        }

        int numPeaks = (int) s->defaultArrayLength;
        os.write(reinterpret_cast<char *>(&numPeaks), sizeIntMSn);

        // end spectrum header info

        // Write out each charge state and corresponding mass
        for(int i = 0; i < numChargeStates; i++)
        {
            os.write(reinterpret_cast<char *>(&(charges[i])), sizeIntMSn);
            os.write(reinterpret_cast<char *>(&(masses[i])), sizeDoubleMSn);
        }
    
        // if there are accurate masses, write out EZ entries
        if( hasAccurateMass ){
          float blank = 0;  // we don't have rTime or area, pad with zeros
          for(int i=0; i < numChargeStates; i++){
            os.write(reinterpret_cast<char *>(&charges[i]), sizeIntMSn);
            os.write(reinterpret_cast<char *>(&masses[i]), sizeDoubleMSn);
            os.write(reinterpret_cast<char *>(&blank), sizeFloatMSn);
            os.write(reinterpret_cast<char *>(&blank), sizeFloatMSn);
          }
        }

        // Do we need to write compressed m/z, intensity arrays?
        if (compress)
        {
            writeCompressedPeaks(s, os);
        }
        else
        {
            // No need to compress, just write out the arrays
            const BinaryDataArray& mzArray = *s->getMZArray();
            const BinaryDataArray& intensityArray = *s->getIntensityArray();
            for(int i = 0; i < numPeaks; i++)
            {
                double mzPeak = mzArray.data[i];
                os.write(reinterpret_cast<char *>(&mzPeak), sizeDoubleMSn);
//.........这里部分代码省略.........

int test()
{
    int failedTests = 0;

    // create the spectrum list
    SpectrumListSimple* sl = new SpectrumListSimple;
    SpectrumListPtr originalList(sl);

    for (size_t i=0; i < testScanSummerSize; ++i)
    {
        TestScanSummerCalculator& t = testScanSummerCalculators[i];
        SpectrumPtr s(new Spectrum);
        s->set(MS_MSn_spectrum);
        s->set(MS_ms_level,2);
        s->index = i;
        s->scanList.scans.push_back(Scan());
        Scan& scanRef = s->scanList.scans[0];
        scanRef.set(MS_scan_start_time,t.rTime,UO_second);
        s->precursors.push_back(Precursor(t.inputPrecursorMZ));
        
        vector<double> inputMZArray = parseDoubleArray(t.inputMZArray);
        vector<double> inputIntensityArray = parseDoubleArray(t.inputIntensityArray);
        s->setMZIntensityArrays(inputMZArray, inputIntensityArray, MS_number_of_detector_counts);
        s->defaultArrayLength = inputMZArray.size();

        scanRef.scanWindows.push_back(ScanWindow());
        scanRef.scanWindows[0].set(MS_scan_window_lower_limit,inputMZArray[0]);
        scanRef.scanWindows[0].set(MS_scan_window_upper_limit,inputMZArray[inputMZArray.size()-1]);

        sl->spectra.push_back(s);
    }


    vector<double> goldMZArray = parseDoubleArray(goldStandard[0].inputMZArray);
    vector<double> goldIntensityArray = parseDoubleArray(goldStandard[0].inputIntensityArray);

    // run spectral summation
    try
    {

        SpectrumListPtr calculator(new SpectrumList_ScanSummer(originalList,0.05,10));

        for (size_t i=0; i < calculator->size(); ++i) 
        {
            SpectrumPtr s = calculator->spectrum(i,true);
            vector<double>& mzs = s->getMZArray()->data;
            vector<double>& intensities = s->getIntensityArray()->data;
            Precursor& precursor = s->precursors[0];
            SelectedIon& selectedIon = precursor.selectedIons[0];
            double precursorMZ = selectedIon.cvParam(MS_selected_ion_m_z).valueAs<double>();
            double rTime = s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds();

            vector<double> goldMZArray = parseDoubleArray(goldStandard[i].inputMZArray);
            vector<double> goldIntensityArray = parseDoubleArray(goldStandard[i].inputIntensityArray);

            unit_assert(mzs.size() == goldMZArray.size());
            unit_assert(intensities.size() == goldIntensityArray.size());
            unit_assert(precursorMZ == goldStandard[i].inputPrecursorMZ);
            unit_assert(rTime == goldStandard[i].rTime);

            for (size_t j=0; j < mzs.size(); ++j)
            {
                unit_assert_equal(mzs[j], goldMZArray[j], 1e-5);
                unit_assert_equal(intensities[j], goldIntensityArray[j], 1e-5);
            }

        }

        
    
            
    }
    catch (exception& e)
    {
        cerr << "Test failed:\n" << e.what() << endl;
        ++failedTests;
    }
    return failedTests;
}

//.........这里部分代码省略.........
    double startMz, stopMz;
    msExperiment->getAcquisitionMassRange(startMz, stopMz);
    scan.scanWindows.push_back(ScanWindow(startMz, stopMz, MS_m_z));

    // decide whether to use Points or Peaks to populate data arrays
    bool doCentroid = msLevelsToCentroid.contains(msLevel);

    bool continuousData = spectrum->getDataIsContinuous();
    if (continuousData && !doCentroid)
        result->set(MS_profile_spectrum);
    else
    {
        result->set(MS_centroid_spectrum);
        doCentroid = continuousData;
    }

    /*if (experimentType == MRM)
    {
        MRMTransitions^ transitions = msExperiment->MRMTransitions;
        double q1mz = transitions[ie.transition]->Q1Mass;//ie.transition->first;
        double q3mz = transitions[ie.transition]->Q3Mass;
        double intensity = points[ie.transition]->Y;
        result->defaultArrayLength = 1;//ie.transition->second.size();

        Precursor precursor;
        SelectedIon selectedIon;

        selectedIon.set(MS_selected_ion_m_z, q1mz, MS_m_z);

        precursor.activation.set(MS_CID); // assume CID

        precursor.selectedIons.push_back(selectedIon);
        result->precursors.push_back(precursor);

        if (getBinaryData)
        {
            mzArray.resize(result->defaultArrayLength, q3mz);
            intensityArray.resize(result->defaultArrayLength, intensity);
        }
    }
    else*/
    {
        if (spectrum->getHasPrecursorInfo())
        {
            double selectedMz, intensity;
            int charge;
            spectrum->getPrecursorInfo(selectedMz, intensity, charge);

            Precursor precursor;
            SelectedIon selectedIon;

            selectedIon.set(MS_selected_ion_m_z, selectedMz, MS_m_z);
            if (charge > 0)
                selectedIon.set(MS_charge_state, charge);

            precursor.activation.set(MS_beam_type_collision_induced_dissociation); // assume beam-type CID since all ABI instruments that write WIFFs are either QqTOF or QqLIT

            precursor.selectedIons.push_back(selectedIon);

            // fix isolation width
            if (spectrum->getHasIsolationInfo())
            {
                double centerMz, lowerLimit, upperLimit;
                spectrum->getIsolationInfo(centerMz, lowerLimit, upperLimit);

                precursor.set(MS_isolation_window_target_m_z, centerMz, MS_m_z);
                precursor.set(MS_isolation_window_upper_offset, upperLimit);
                precursor.set(MS_isolation_window_lower_offset, lowerLimit);
            }


            result->precursors.push_back(precursor);
        }

        //virtual bool getHasIsolationInfo() const;
        //virtual void getIsolationInfo(double& centerMz, double& lowerLimit, double& upperLimit) const;


        //result->set(MS_lowest_observed_m_z, spectrum->getMinX(), MS_m_z);
        //result->set(MS_highest_observed_m_z, spectrum->getMaxX(), MS_m_z);
        result->set(MS_base_peak_intensity, spectrum->getBasePeakY(), MS_number_of_detector_counts);
        result->set(MS_base_peak_m_z, spectrum->getBasePeakX(), MS_m_z);
        result->set(MS_total_ion_current, spectrum->getSumY(), MS_number_of_detector_counts);

        if (getBinaryData)
        {
            result->setMZIntensityArrays(std::vector<double>(), std::vector<double>(), MS_number_of_detector_counts);
            BinaryDataArrayPtr mzArray = result->getMZArray();
            BinaryDataArrayPtr intensityArray = result->getIntensityArray();

            spectrum->getData(doCentroid, mzArray->data, intensityArray->data);
            if (doCentroid)
                result->set(MS_profile_spectrum); // let SpectrumList_PeakPicker know this was a profile spectrum
        }

        result->defaultArrayLength = spectrum->getDataSize(doCentroid);
    }

    return result;
}

PWIZ_API_DECL SpectrumPtr SpectrumList_PeakPicker::spectrum(size_t index, bool getBinaryData) const
{
    SpectrumPtr s;
    
    switch (mode_)
    {
        case 1:
            s = dynamic_cast<detail::SpectrumList_Thermo*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 2:
            s = dynamic_cast<detail::SpectrumList_Bruker*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 3:
            s = dynamic_cast<detail::SpectrumList_ABI*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 4:
            s = dynamic_cast<detail::SpectrumList_Agilent*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 5:
            s = dynamic_cast<detail::SpectrumList_ABI_T2D*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 0:
        default:
            s = inner_->spectrum(index, true);
            break;
    }

    if (!msLevelsToPeakPick_.contains(s->cvParam(MS_ms_level).valueAs<int>()))
        return s;

    vector<CVParam>& cvParams = s->cvParams;
    vector<CVParam>::iterator itr = std::find(cvParams.begin(), cvParams.end(), MS_profile_spectrum);

    // return non-profile spectra as-is
    // (could have been acquired as centroid, or vendor may have done the centroiding)
    if (itr == cvParams.end())
        return s;

    // make sure the spectrum has binary data
    if (!s->getMZArray().get() || !s->getIntensityArray().get())
        s = inner_->spectrum(index, true);

    // replace profile term with centroid term
    vector<CVParam>& cvParams2 = s->cvParams;
    *(std::find(cvParams2.begin(), cvParams2.end(), MS_profile_spectrum)) = MS_centroid_spectrum;

    try
    {
        vector<double>& mzs = s->getMZArray()->data;
        vector<double>& intensities = s->getIntensityArray()->data;
        vector<double> xPeakValues, yPeakValues;
        algorithm_->detect(mzs, intensities, xPeakValues, yPeakValues);
        mzs.swap(xPeakValues);
        intensities.swap(yPeakValues);
        s->defaultArrayLength = mzs.size();
    }
    catch(std::exception& e)
    {
        throw std::runtime_error(std::string("[SpectrumList_PeakPicker::spectrum()] Error picking peaks: ") + e.what());
    }

    s->dataProcessingPtr = dp_;
    return s;
}

PWIZ_API_DECL void ThresholdFilter::operator () (const SpectrumPtr s) const
{
    if (!msLevelsToThreshold.contains(s->cvParam(MS_ms_level).valueAs<int>()))
        return;

    // do nothing to empty spectra
    if (s->defaultArrayLength == 0)
        return;

    if (byType == ThresholdingBy_Count ||
        byType == ThresholdingBy_CountAfterTies)
    {
        // if count threshold is greater than number of data points, return as is
        if (s->defaultArrayLength <= threshold)
            return;
        else if (threshold == 0)
        {
            s->getMZArray()->data.clear();
            s->getIntensityArray()->data.clear();
            s->defaultArrayLength = 0;
            return;
        }
    }

    vector<MZIntensityPair> mzIntensityPairs;
    s->getMZIntensityPairs(mzIntensityPairs);

    if (orientation == Orientation_MostIntense)
        sort(mzIntensityPairs.begin(), mzIntensityPairs.end(), orientationMore_Predicate);
    else if (orientation == Orientation_LeastIntense)
        sort(mzIntensityPairs.begin(), mzIntensityPairs.end(), orientationLess_Predicate);
    else
        throw runtime_error("[threshold()] invalid orientation type");

    double tic = accumulate(mzIntensityPairs.begin(), mzIntensityPairs.end(), 0.0, MZIntensityPairIntensitySum());

    if (tic == 0)
    {
        s->getMZArray()->data.clear();
        s->getIntensityArray()->data.clear();
        s->defaultArrayLength = 0;
        return;
    }

    double bpi = orientation == Orientation_MostIntense ? mzIntensityPairs.front().intensity
                                                        : mzIntensityPairs.back().intensity;

    // after the threshold is applied, thresholdItr should be set to the first data point to erase
    vector<MZIntensityPair>::iterator thresholdItr;

    switch (byType)
    {
        case ThresholdingBy_Count:
            // no need to check bounds on thresholdItr because it gets checked above
            thresholdItr = mzIntensityPairs.begin() + (size_t) threshold;

            // iterate backward until a non-tie is found
            while (true)
            {
                const double& i = thresholdItr->intensity;
                if (thresholdItr == mzIntensityPairs.begin())
                    break;
                else if (i != (--thresholdItr)->intensity)
                {
                    ++thresholdItr;
                    break;
                }
            }
            break;

        case ThresholdingBy_CountAfterTies:
            // no need to check bounds on thresholdItr because it gets checked above
            thresholdItr = mzIntensityPairs.begin() + ((size_t) threshold)-1;

            // iterate forward until a non-tie is found
            while (true)
            {
                const double& i = thresholdItr->intensity;
                if (++thresholdItr == mzIntensityPairs.end() ||
                    i != thresholdItr->intensity)
                    break;
            }
            break;

        case ThresholdingBy_AbsoluteIntensity:
            if (orientation == Orientation_MostIntense)
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
                                           mzIntensityPairs.end(),
                                           MZIntensityPair(0, threshold),
                                           orientationMore_Predicate);
            else
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
                                           mzIntensityPairs.end(),
                                           MZIntensityPair(0, threshold),
                                           orientationLess_Predicate);
            break;

        case ThresholdingBy_FractionOfBasePeakIntensity:
            if (orientation == Orientation_MostIntense)
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                //in some case (i suppose...) a previous peak could be larger
                //that is a reason to not reset the intensity to 0 too
                //if (! binnedSpectrum.getIntensityFromIndex(index_right))
                //binnedSpectrum.getOrInitIntensityFromIndex(index_right);
                /*} else {
                    binnedSpectrum.getOrInitIntensityFromIndex(index + 1);
                }*/
            }

            //pair<vector<double>, vector<double> > data = binnedSpectrum.getData();
            //printf("%d, %d\n", intData.size(), binnedSpectrum.mzData.size());
            ptr->setMZIntensityArrays(binnedSpectrum.mzData, intData,  CVID_Unknown);

            //here binning
            /*---comparing data intensities*/
            //printf("binnedSpectrumSize:%d, rawSize:%d\n", data.first.size(), rawSpectrum->getMZArray()->data.size());

            //bin
            //int nbBins = (int) (1700.0 / 0.1);
            map<int, vector<pair<double, double> > > rawIntensitiesSumByBinIndex, mzdbIntensitiesSumByIndex;
            vector<MZIntensityPair> rawMzIntensities;
            rawSpectrum->getMZIntensityPairs(rawMzIntensities);
            for (size_t j = 0; j < rawMzIntensities.size(); ++j) {
                const MZIntensityPair& pair = rawMzIntensities[j];
                int idx = (int) (pair.mz / 1);
                rawIntensitiesSumByBinIndex[idx].push_back(make_pair(pair.mz, pair.intensity));
            }

            for (size_t j= 0; j < intData.size(); ++j) {
                const double& mz = binnedSpectrum.mzData[j];
                const double& intensity = intData[j];
                int idx = (int) (mz / 1);
                mzdbIntensitiesSumByIndex[idx].push_back(make_pair(mz, intensity));
            }

            //iterate on raw map
            double counter = 0;
            double sum = 0;
            for (auto it = rawIntensitiesSumByBinIndex.begin(); it != rawIntensitiesSumByBinIndex.end(); ++it) {
                const int& idx = it->first;
                const vector<pair<double, double> >& vec = it->second;
                double integratedIntensity = 0;
                if (! vec.empty()) {
                    integrate(vec, integratedIntensity);
                }
                if (integratedIntensity) {
                    if (mzdbIntensitiesSumByIndex.find(idx) != mzdbIntensitiesSumByIndex.end()) {
                        const vector<pair<double, double> >& mzdbPoints = mzdbIntensitiesSumByIndex[idx];
                        //printf("v1 %f, %f\n", intensity, mzdbIntensity);
                        double mzDBIntegratedIntensity = 0;
                        integrate(mzdbPoints, mzDBIntegratedIntensity);
                        double error = (integratedIntensity - mzDBIntegratedIntensity) / integratedIntensity;
                        sum += error;
                        counter ++;
                    }
                }
            }
           //printf("%f, %f\n", sum, counter);
            /*if (scan->idMzDB == 5524) {
                printf("Scan 5524:%d, reconstructed profile:%d, raw profile:%d\n", scan->mz.size(), data.first.size(), rawSpectrum->getMZArray()->data.size());
            }*/
            fileHandle << (sum / counter) * 100.0 << endl;

            if (rawSpectrum->index == 1) {
                auto& rawData = rawSpectrum->getIntensityArray()->data;

                for (size_t kk = 0; kk < binnedSpectrum.mzData.size(); ++kk) {
                    //fileHandle_2 << binnedSpectrum.mzData[kk] << "\t" << rawData[kk] << "\t" << intData[kk] << endl;
                    fprintf(fileHandle_2, "%.4f\t%f\t%f\n", binnedSpectrum.mzData[kk] ,rawData[kk] , intData[kk] );
                }
            }


            //fileHandle_2 << ( (rawNonZeroCount - mzdbNonZeroCount) / rawNonZeroCount ) * 100 << endl;

        } else {
            ptr->setMZIntensityArrays(mz, intens, CVID_Unknown);
        }
        cache->addKeyValue(boost::lexical_cast<string>(i), ptr);

        //clear data
        //clearScanData(scan);

        int newPercent = (int) (((float) (i) / spec->size() * 100));
        if (newPercent != lastPercent) {
            printProgBar(newPercent);
            lastPercent = newPercent;
        }



    }//end iteration mzdb spectra;
    printProgBar(100);
    fileHandle.close(); fclose(fileHandle_2);/*---do not forget to close the fileHandle*/
    fclose(fileHandle_3);
    SpectrumListPtr spectrumListCached(new mzSpectrumListCache(f, mzdb.get(), cache));
    mzdb->run.spectrumListPtr = spectrumListCached;
    ChromatogramListPtr chromListPtr(new mzEmptyChromatogram);
    mzdb->run.chromatogramListPtr = chromListPtr;
}

      void mzDataCache::addKeyValue(const string& s, SpectrumPtr& spec) {
			vector<double>& mz = spec->getMZArray()->data;
			vector<double>& intensity = spec->getIntensityArray()->data;
			db.set((s + "mz").c_str(), (s + "mz").size() * sizeof (char), (char*) &mz[0], mz.size() * sizeof (double));
			db.set((s + "int").c_str(), (s + "int").size() * sizeof (char), (char*) &intensity[0], intensity.size() * sizeof (double));
		}