Changeset 8072 for trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h

Timestamp:

10/16/06 18:02:31 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h (modified) (5 diffs)

trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h

@@ -5 +5 @@
 #include <TROOT.h>
 #endif
+
+class TH1;
+class TH2;
+class TH1F;
+class TH2F;
+class TArrayF;
 
 class MExtralgoDigitalFilter
@@ -13 +19 @@
     Int_t    fNum;
 
-    Float_t *fWeightsAmp;
-    Float_t *fWeightsTime;
+    Float_t const *fWeightsAmp;
+    Float_t const *fWeightsTime;
 
-    Int_t   fWeightsPerBin; // Number of weights per data bin
+    const Int_t fWeightsPerBin; // Number of weights per data bin
+    const Int_t fWindowSize;
 
     // Result
@@ -24 +31 @@
     Float_t fSignalDev;
 
-    inline void Eval(Double_t &sumamp, Double_t &sumtime, const Int_t windowsize, const Int_t startv, const Int_t startw=0) const
+    // Weights: Weights to evaluate
+    // Startv: Index of first bin of data
+    // startw: Offset on the weights
+    inline Double_t Eval(Float_t const *weights, const Int_t startv, const Int_t startw=0) const
     {
         //
@@ -30 +40 @@
         // and multiply the entries with the corresponding weights
         //
-        sumamp = 0;
-        sumtime = 0;
+        Double_t sum = 0;
 
         // Shift the start of the weight to the center of sample 0
-        Float_t const *wa  = fWeightsAmp  + fWeightsPerBin/2 + startw;
-        Float_t const *wt  = fWeightsTime + fWeightsPerBin/2 + startw;
+        Float_t const *w = weights + startw;
+
         Float_t *const beg = fVal+startv;
+        for (Float_t const *pex=beg; pex<beg+fWindowSize; pex++)
+        {
+            sum += *w * *pex;
+            w += fWeightsPerBin;
+        }
+        return sum;
+    }
 
-        for (Float_t const *pex=beg; pex<beg+windowsize; pex++)
+    inline void AlignIntoLimits(Int_t &maxp, Int_t &frac) const
+    {
+        // Align maxp into available range (TO BE CHECKED)
+        if (maxp < 0)
         {
-            sumamp  += *wa * *pex;
-            sumtime += *wt * *pex;
-
-            // Step forward by one bin
-            wa += fWeightsPerBin;
-            wt += fWeightsPerBin;
+            maxp = 0;
+            frac = fWeightsPerBin/2-1; // Assume peak at the end of the last slice
+        }
+        if (maxp > fNum-fWindowSize)
+        {
+            maxp = fNum-fWindowSize;
+            frac = -fWeightsPerBin/2; // Assume peak at the beginning of the first slice
         }
     }
 
+    inline Int_t AlignExtractionWindow(Int_t &maxp, Int_t &frac, const Double_t ampsum)
+    {
+        // Align extraction window to signal position
+
+        const Double_t timesum = Eval(fWeightsTime, maxp, frac);
+
+        // Because fWeightsPerBin/2 doesn't correspond to the center
+        // of a bin the time-values extracted are slightly positive.
+        // They are roughly between -0.45 and 0.55
+        const Double_t binoffset = TMath::Even(fWeightsPerBin) ? 0.5 : 0;
+
+        // This is the time offset from the extraction position
+        Double_t tmoffset = (frac+binoffset)/fWeightsPerBin + timesum/ampsum;
+
+        // Convert the residual fraction of one slice into an
+        // offset position in the extraction weights
+        const Int_t integ = TMath::FloorNint(tmoffset+0.5);
+
+        /*
+        if (integ>0)
+            tmoffset=0.49-0.05;
+        if (integ<0)
+            tmoffset=-0.49-0.05;
+        integ=0;
+        */
+
+        // move the extractor by an offset number of slices
+        // determined by the extracted time
+        maxp -= integ;
+
+        frac  = TMath::FloorNint((tmoffset-integ)*fWeightsPerBin);
+
+        // Align maxp into available range (TO BE CHECKED)
+        AlignIntoLimits(maxp, frac);
+
+        return integ;
+    }
+
+    inline void AlignExtractionWindow(Int_t &maxp, Int_t &frac)
+    {
+        const Double_t amp = Eval(fWeightsAmp, maxp, frac);
+        if (amp!=0)
+            AlignExtractionWindow(maxp, frac, amp);
+    }
+
 public:
-    MExtralgoDigitalFilter(Float_t *val, Int_t n, Float_t *wa, Float_t *wt)
-        : fVal(val), fNum(n),
-        fWeightsAmp(wa), fWeightsTime(wt), fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
+    MExtralgoDigitalFilter(Int_t res, Int_t windowsize, Float_t *wa, Float_t *wt)
+        : fVal(0), fNum(0), fWeightsAmp(wa+res/2), fWeightsTime(wt+res/2),
+        fWeightsPerBin(res), fWindowSize(windowsize),
+        fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
     {
     }
 
-    void SetWeightsPerBin(Int_t res) { fWeightsPerBin=res; }
+    void SetData(Int_t n, Float_t *val) { fNum=n; fVal=val; }
 
     Float_t GetTime() const          { return fTime; }
@@ -67 +133 @@
     void GetTime(Float_t &sig, Float_t &dsig) const   { sig=fTime; dsig=fTimeDev; }
 
-    Float_t ExtractNoise(Int_t iter, Int_t windowsize) const;
-    void Extract(Int_t windowsize, Float_t timeshift);
+    Float_t ExtractNoise(Int_t iter) const;
+    void Extract();
+
+    static Bool_t CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
+    static void CalculateWeights2(TH1F &shape, const TH2F &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
 };