source: trunk/Mars/mcore/DrsCalib.h@ 12278

class CalibData
{
protected:
    uint64_t fNumEntries;

    size_t fNumSamples;
    size_t fNumChannels;

    vector<int64_t> fSum;
    vector<int64_t> fSum2;

    /*
    vector<map<int16_t, uint8_t> > fMedian;
    vector<int16_t> fResult;
     */
public:
    CalibData() : fNumEntries(0), fNumSamples(0), fNumChannels(0) { }
    void Reset()
    {
        fNumEntries  = 0;
        fNumSamples  = 0;
        fNumChannels = 0;

        fSum.clear();
        fSum2.clear();
    }

    void InitSize(uint16_t channels, uint16_t samples)
    {
        fNumChannels = channels;
        fNumSamples  = samples;

        fSum.resize(samples*channels);
        fSum2.resize(samples*channels);

        //fMedian.resize(40*9*1024);
    }
/*
                continue;

                if (ch<40*9)
                {
                    fMedian[abs][val[rel]]++;

                    int n= 8;
                    if (fNumEntries>0 && fNumEntries%(n*100)==0)
                    {
                        fResult.resize(40*9*1024);

                        uint16_t sum = 0;
                        for (map<int16_t, uint8_t>::const_iterator it=fMedian[abs].begin();
                             it!=fMedian[abs].end(); it++)
                        {
                            map<int16_t, uint8_t>::const_iterator is = it;
                            is++;
                            if (is==fMedian[abs].end())
                                break;

                            sum += it->second;

                            double med = 0;

                            med += int64_t(it->second)*int64_t(it->first);
                            med += int64_t(is->second)*int64_t(is->first);
                            med /= int64_t(it->second)+int64_t(is->second);

                            if (sum==n*50)
                            {
                                fResult[abs] = it->first;
                                cout << ch << " MED+(50): " << it->first << endl;
                            }
                            if (sum<n*50 && sum+is->second>n*50)
                            {
                                fResult[abs] = med;
                                cout << ch << " MED-(50): " << med << endl;
                            }
                        }
                        cout << ch << " AVG=" << float(fSum[abs])/fNumEntries << endl;
                        fMedian[abs].clear();
                    }
                } // -2029 -2012 -2003 -1996 -1990 // -1955
                */

    void AddRel(const int16_t *val, const int16_t *start)
    {
        for (size_t ch=0; ch<fNumChannels; ch++)
        {
            const int16_t spos = start[ch];
            if (spos<0)
                continue;

            const size_t pos = ch*fNumSamples;

            for (size_t i=0; i<fNumSamples; i++)
            {
                // Value is relative to trigger
                // Abs is corresponding index relative to DRS pipeline
                const size_t rel = pos +  i;
                const size_t abs = pos + (spos+i)%fNumSamples;

                const int64_t v = val[rel];

                fSum[abs]  += v;
                fSum2[abs] += v*v;
            }
        }

        fNumEntries++;
    }

    void AddRel(const int16_t *val,    const int16_t *start,
                const int32_t *offset, const uint32_t scale)
    {
        for (size_t ch=0; ch<fNumChannels; ch++)
        {
            const int16_t spos = start[ch];
            if (spos<0)
                continue;

            const size_t pos = ch*fNumSamples;

            for (size_t i=0; i<fNumSamples; i++)
            {
                // Value is relative to trigger
                // Offset is relative to DRS pipeline
                // Abs is corresponding index relative to DRS pipeline
                const size_t rel = pos +  i;
                const size_t abs = pos + (spos+i)%fNumSamples;

                const int64_t v = int64_t(val[rel])*scale-offset[abs];

                fSum[abs]  += v;
                fSum2[abs] += v*v;
            }
        }

        fNumEntries++;
    }

    void AddAbs(const int16_t *val,    const  int16_t *start,
                const int32_t *offset, const uint32_t scale)
    {
        for (size_t ch=0; ch<fNumChannels; ch++)
        {
            const int16_t spos = start[ch];
            if (spos<0)
                continue;

            const size_t pos = ch*fNumSamples;

            for (size_t i=0; i<fNumSamples; i++)
            {
                // Value is relative to trigger
                // Offset is relative to DRS pipeline
                // Abs is corresponding index relative to DRS pipeline
                const size_t rel = pos +  i;
                const size_t abs = pos + (spos+i)%fNumSamples;

                const int64_t v = int64_t(val[rel])*scale-offset[abs];

                fSum[rel]  += v;
                fSum2[rel] += v*v;
            }
        }

        fNumEntries++;
    }

    static void Apply(int16_t *val, const int16_t *start, uint32_t roi,
                      const int32_t *offset, const uint32_t scaleabs,
                      const int64_t *gain,   const uint32_t scalegain,
                      const int64_t *trgoff, const uint32_t scalerel)
    {
        return ;// "FIXME: scalegain exceeds valid range..."
        for (size_t ch=0; ch<1440; ch++)
        {
            const int16_t spos = start[ch];
            if (spos<0)
                continue;

            const size_t pos = ch*roi;

            for (size_t i=0; i<roi; i++)
            {
                // Value is relative to trigger
                // Offset is relative to DRS pipeline
                // Abs is corresponding index relative to DRS pipeline
                const size_t rel = pos +  i;
                const size_t abs = pos + (spos+i)%1024;

                const int64_t v =
                    + int64_t(val[rel])   *scaleabs*scalerel
                    - int64_t(trgoff[rel])*scaleabs
                    - int64_t(offset[abs])*scalerel
                    ;

                // The gain value is multiplied by fNumOffset and fNumGain
                const int64_t div = int64_t(gain[abs])*scaleabs*scalerel;
                val[rel] = (v*scalegain)/div;
             }
        }
    }

    static void Apply(float *vec, const int16_t *val, const int16_t *start, uint32_t roi,
                      const int32_t *offset, const uint32_t scaleabs,
                      const int64_t *gain,   /*const*/ uint64_t scalegain,
                      const int64_t *trgoff, /*const*/ uint64_t scalerel)
    {
        /*
        scalegain *= scaleabs;
        scalerel  *= scaleabs;

        for (size_t ch=0; ch<1440; ch++)
        {
            const size_t pos = ch*roi;
            const int16_t spos = start[ch];
            if (spos<0)
            {
                memset(vec+pos, 0, roi);
                continue;
            }

            for (size_t i=0; i<roi; i++)
            {
                // Value is relative to trigger
                // Offset is relative to DRS pipeline
                // Abs is corresponding index relative to DRS pipeline
                const size_t rel = pos +  i;
                const size_t abs = pos + (spos+i)%1024;

                const int64_t v =
                    + int64_t(val[rel])   *scaleabs*scalerel
                    - int64_t(trgoff[rel])*scaleabs
                    - int64_t(offset[abs])*scalerel
                    ;

                const int64_t div = int64_t(gain[abs])*scaleabs*scalerel;

                vec[rel] = double(v)*scalegain/div;
            }
        }*/

        for (size_t ch=0; ch<1440; ch++)
        {
            const size_t pos = ch*roi;

            const int16_t spos = start[ch];
            if (spos<0)
            {
                memset(vec+pos, 0, roi);
                continue;
            }

            for (size_t i=0; i<roi; i++)
            {
                // Value is relative to trigger
                // Offset is relative to DRS pipeline
                // Abs is corresponding index relative to DRS pipeline
                const size_t rel = pos +  i;
                const size_t abs = pos + (spos+i)%1024;

                const int64_t v =
                    + (int64_t(val[rel])*scaleabs-offset[abs])*scalerel
                    - trgoff[rel]
                    ;

                const int64_t div = gain[abs]*scalerel;
                vec[rel] = double(v)*scalegain/div;
            }
        }
    }

    static void RemoveSpikes(float *vec, uint32_t roi)
    {
        if (roi<4)
            return;

        cout << "X" << endl;

        for (size_t ch=0; ch<1440; ch++)
        {
            float *p = vec + ch*roi;

            for (size_t i=1; i<roi-2; i++)
            {
                if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
                {
                    p[i] = (p[i-1]+p[i+1])/2;
                }

                if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
                {
                    p[i] = (p[i-1]+p[i+2])/2;
                    p[i+1] = p[i];
                }
            }
        }
    }

    pair<vector<double>,vector<double> > GetSampleStats() const
    {
        if (fNumEntries==0)
            return make_pair(vector<double>(),vector<double>());

        vector<double> mean(fSum.size());
        vector<double> error(fSum.size());

        vector<int64_t>::const_iterator it = fSum.begin();
        vector<int64_t>::const_iterator i2 = fSum2.begin();
        vector<double>::iterator        im = mean.begin();
        vector<double>::iterator        ie = error.begin();

        while (it!=fSum.end())
        {
            *im = /*cnt<fResult.size() ? fResult[cnt] :*/ double(*it)/fNumEntries;
            *ie = sqrt(double(*i2*int64_t(fNumEntries) - *it * *it))/fNumEntries;

            im++;
            ie++;
            it++;
            i2++;
        }


        /*
         valarray<double> ...

         mean /= fNumEntries;
         error = sqrt(error/fNumEntries - mean*mean);
         */

        return make_pair(mean, error);
    }

    void GetSampleStats(float *ptr, float scale) const
    {
        if (fNumEntries==0)
        {
            memset(ptr, 0, sizeof(float)*1024*1440*2);
            return;
        }

        vector<int64_t>::const_iterator it = fSum.begin();
        vector<int64_t>::const_iterator i2 = fSum2.begin();

        while (it!=fSum.end())
        {
            *ptr             = scale*double(*it)/fNumEntries;
            *(ptr+1024*1440) = scale*sqrt(double(*i2*int64_t(fNumEntries) - *it * *it))/fNumEntries;

            ptr++;
            it++;
            i2++;
        }
    }

    static void GetPixelStats(float *ptr, const float *data, uint16_t roi)
    {
        if (roi==0)
            return;

        for (int i=0; i<1440; i++)
        {
            const float *vec = data+i*roi;

            int    pos  = 0;
            double sum  = vec[0];
            double sum2 = vec[0]*vec[0];
            for (int j=1; j<roi; j++)
            {
                sum  += vec[j];
                sum2 += vec[j]*vec[j];

                if (vec[j]>vec[pos])
                    pos = j;
            }
            sum  /= roi;
            sum2 /= roi;

            *(ptr+0*1440+i) = sum;
            *(ptr+1*1440+i) = sqrt(sum2 - sum * sum);
            *(ptr+2*1440+i) = vec[pos];
            *(ptr+3*1440+i) = pos;
        }
    }

    static void GetPixelMax(float *max, const float *data, uint16_t roi, int32_t first, int32_t last)
    {
        if (roi==0 || first<0 || last<0 || first>=roi || last>=roi || last<first)
            return;

        for (int i=0; i<1440; i++)
        {
            const float *beg = data+i*roi+first;
            const float *end = data+i*roi+last;

            const float *pmax = beg;

            for (const float *ptr=beg+1; ptr<=end; ptr++)
                if (*ptr>*pmax)
                    pmax = ptr;

            max[i] = *pmax;
        }

    }

    const vector<int64_t> &GetSum() const { return fSum; }

    uint64_t GetNumEntries() const { return fNumEntries; }
};