Changeset 17862 for trunk/Mars

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:

• trunk/Mars/msignal/MTreatSaturation.cc (modified) (2 diffs)

trunk/Mars/msignal/MTreatSaturation.cc

@@ -111 +111 @@
         // Find maximum
         const Float_t *pmax = pbeg;
-
         for (const Float_t *p=pbeg; p<pend; p++)
             if (*p>*pmax)
                 pmax = p;
 
+        // Is there any chance for saturation?
         if (*pmax<1800)
             continue;
 
-        MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray());
-        s.SetExtractionType(MExtralgoSpline::kAmplitudeRel);
-
-        const double leading   = s.SearchYdn(1800., pmax-pbeg);
-        const double falling   = s.SearchYup(1800., pmax-pbeg);
-        const double width     = falling-leading;
-        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);
-        const double time      = leading-(4.46944 - 0.277272*width + 0.0035479*width*width);
-
+        // Determine a rough estimate for the average baseline
         double baseline = 0;
         for (const Float_t *p=ptr+10; p<ptr+50; p++)
@@ -133 +125 @@
         baseline /= 40;
 
-        const UInt_t beg = TMath::CeilNint(leading);
-        const UInt_t end = TMath::FloorNint(falling);
+        // Do a spline interpolation to find the crossing of 1.8V
+        // before and after the maximum
+        MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray());
+        s.SetExtractionType(MExtralgoSpline::kAmplitudeRel);
 
-        ptr += fHiGainFirst+1;
-        for (UInt_t i=0; i<end-beg; i++)
+        const double leading = s.SearchYdn(pmax-pbeg, 1800);
+        const double falling = s.SearchYup(pmax-pbeg, 1800);
+        const double width   = falling-leading;
+
+        // If the width above 1.8V is below 3 samples... go ahead as usual.
+        if (width>3)
         {
-            cout << i << " " << ptr[beg-i] << " ";
-            ptr[beg+i] = amplitude*(1-1/(1+exp((i-time)/2.14)))*exp((time-i)/38.8)+baseline;
-            cout << ptr[beg+i] << endl;
+            // Estimate the amplitude and the arrival time from the width
+            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);
+            const double deltat    = -1.41371-0.0525846*width + 93.2763/(width+13.196);
+            const double time0     = leading-deltat-1;
+
+            const Int_t beg = TMath::FloorNint(leading);
+            const Int_t end = TMath::CeilNint(falling);
+
+            ptr += fHiGainFirst+1;
+            for (Int_t i=beg-5; i<end+5; i++)
+            {
+                const double x = i-time0;
+                const double v = amplitude*(1-1/(1+exp(x/2.14)))*exp(-x/38.8)+baseline;
+                if (v>ptr[i])
+                    ptr[i] = v;
+            }
         }
-
-        cout << endl;
     }