Changeset 17862 for trunk/Mars/msignal


Ignore:
Timestamp:
05/15/14 19:09:49 (11 years ago)
Author:
tbretz
Message:
That's the best algorithm I could find within limited time. I think it is precise to a few percent level (5-15%)
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/msignal/MTreatSaturation.cc

    r17858 r17862  
    111111        // Find maximum
    112112        const Float_t *pmax = pbeg;
    113 
    114113        for (const Float_t *p=pbeg; p<pend; p++)
    115114            if (*p>*pmax)
    116115                pmax = p;
    117116
     117        // Is there any chance for saturation?
    118118        if (*pmax<1800)
    119119            continue;
    120120
    121         MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray());
    122         s.SetExtractionType(MExtralgoSpline::kAmplitudeRel);
    123 
    124         const double leading   = s.SearchYdn(1800., pmax-pbeg);
    125         const double falling   = s.SearchYup(1800., pmax-pbeg);
    126         const double width     = falling-leading;
    127         const double amplitude = 1800/(0.898417 - 0.0187633*width + 0.000163919*width*width - 6.87369e-7*width*width*width + 1.13264e-9*width*width*width*width);
    128         const double time      = leading-(4.46944 - 0.277272*width + 0.0035479*width*width);
    129 
     121        // Determine a rough estimate for the average baseline
    130122        double baseline = 0;
    131123        for (const Float_t *p=ptr+10; p<ptr+50; p++)
     
    133125        baseline /= 40;
    134126
    135         const UInt_t beg = TMath::CeilNint(leading);
    136         const UInt_t end = TMath::FloorNint(falling);
     127        // Do a spline interpolation to find the crossing of 1.8V
     128        // before and after the maximum
     129        MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray());
     130        s.SetExtractionType(MExtralgoSpline::kAmplitudeRel);
    137131
    138         ptr += fHiGainFirst+1;
    139         for (UInt_t i=0; i<end-beg; i++)
     132        const double leading = s.SearchYdn(pmax-pbeg, 1800);
     133        const double falling = s.SearchYup(pmax-pbeg, 1800);
     134        const double width   = falling-leading;
     135
     136        // If the width above 1.8V is below 3 samples... go ahead as usual.
     137        if (width>3)
    140138        {
    141             cout << i << " " << ptr[beg-i] << " ";
    142             ptr[beg+i] = amplitude*(1-1/(1+exp((i-time)/2.14)))*exp((time-i)/38.8)+baseline;
    143             cout << ptr[beg+i] << endl;
     139            // Estimate the amplitude and the arrival time from the width
     140            const double amplitude = (1800-baseline)/(0.898417 - 0.0187633*width + 0.000163919*width*width - 6.87369e-7*width*width*width + 1.13264e-9*width*width*width*width);
     141            const double deltat    = -1.41371-0.0525846*width + 93.2763/(width+13.196);
     142            const double time0     = leading-deltat-1;
     143
     144            const Int_t beg = TMath::FloorNint(leading);
     145            const Int_t end = TMath::CeilNint(falling);
     146
     147            ptr += fHiGainFirst+1;
     148            for (Int_t i=beg-5; i<end+5; i++)
     149            {
     150                const double x = i-time0;
     151                const double v = amplitude*(1-1/(1+exp(x/2.14)))*exp(-x/38.8)+baseline;
     152                if (v>ptr[i])
     153                    ptr[i] = v;
     154            }
    144155        }
    145 
    146         cout << endl;
    147156    }
    148157
Note: See TracChangeset for help on using the changeset viewer.