Changeset 8154 for trunk/MagicSoft/Mars/mextralgo

Timestamp:

10/24/06 09:26:10 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mextralgo

Files:

• MExtralgoDigitalFilter.cc (modified) (13 diffs)
• MExtralgoDigitalFilter.h (modified) (6 diffs)
• MExtralgoSpline.cc (modified) (11 diffs)
• MExtralgoSpline.h (modified) (3 diffs)

trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc

@@ -38 +38 @@
 }
 
+// -----------------------------------------------------------------------------
+//
+// Calculates the chi2 of the fit, once the weights have been iterated.
+// Argument: time, obtained after a call to EvalDigitalFilterHiGain
+//
+Float_t MExtralgoDigitalFilter::GetChisq(const Int_t maxp, const Int_t frac, const Float_t sum) const
+{
+    /*
+    TMatrix g (windowh,1);
+    TMatrix gt(windowh,1);
+    TMatrix y (windowh,1);
+
+    const Float_t etau   = fFineAdjustHi*sumhi;
+    const Float_t t_fine = TMath::Abs(fFineAdjustHi)< 1./fBinningResolutionHiGain ? -fFineAdjustHi : 0.;
+
+    //   if (t_fine==0.)
+    //     return -1.;
+
+    if (fDebug)
+        gLog << all << " fMaxPHi: " << fMaxPHi << " fIterPHi " << fIterPHi << " t_fine: " << t_fine << endl;
+
+    //
+    // Slide with a window of size windowh over the sample
+    // and calculate the arrays by interpolating the pulse shape using the
+    // fine-tuned time information.
+    //
+    for (Int_t sample=0; sample &lt; windowh; sample++)
+    {
+        const Int_t idx = fArrBinningResHalfHiGain[sample] + fIterPHi;
+        const Int_t ids = fMaxPHi + sample;
+
+        y [sample][0]     = fHiGainSignalDF[ids];
+        g [sample][0]     = t_fine >= 0
+            ? (fPulseShapeHiGain[idx]    + t_fine*(fPulseShapeHiGain[idx+1]   -fPulseShapeHiGain[idx])  )*sumhi
+            : (fPulseShapeHiGain[idx]    + t_fine*(fPulseShapeHiGain[idx]     -fPulseShapeHiGain[idx-1]))*sumhi;
+        gt[sample][0]     = t_fine >= 0
+            ? (fPulseShapeDotHiGain[idx] + t_fine*(fPulseShapeDotHiGain[idx+1]-fPulseShapeDotHiGain[idx])   )*etau
+            : (fPulseShapeDotHiGain[idx] + t_fine*(fPulseShapeDotHiGain[idx]  -fPulseShapeDotHiGain[idx-1]) )*etau;
+    }
+
+    TMatrix second = y - g - gt;
+    TMatrix first(TMatrix::kTransposed,second);
+    TMatrix chisq = first * ((*fBHiInv)*second);
+
+    return chisq[0][0]/(windowh-2);
+    */
+/*
+
+    TMatrix S(fWindowSize, 1); // Signal (start==start of window)
+    for (int i=0; i<fWindowSize; i++)
+        S[i][0] = fVal[i+maxp];
+
+    TMatrix g(fWindowSize, 1);
+    //TMatrix gT(fWindowSize, 1);
+
+    for (int i=0; i<fWindowSize; i++)
+    {
+        Int_t idx = fWeightsPerBin*i + frac;
+
+        // FIXME: Maybe we could do an interpolation on time-fineadj?
+        //Float_t slope  = fPulseShapeHiGain[idx+1]   -fPulseShapeHiGain[idx];
+        //Float_t slopet = fPulseShapeDotHiGain[idx+1]-fPulseShapeDotHiGain[idx];
+
+        g[i][0] = fPulseShapeHiGain[idx]   *sumhi;
+        //gT[i][0] = fPulseShapeHiGainDot[idx]*tau;
+    }
+
+    TMatrix Ainv; // Autocorrelation Matrix (inverted)
+
+    TMatrix m = S - g;// - gT;
+    TMatrix mT(TMatrix::kTransposed, m);
+
+    TMatrix chisq = mT * (Ainv*m);
+
+    return chisq[0][0]/(fWindowSize-2);
+  */
+
+    Double_t sumc = 0;
+
+    TMatrix d(fWindowSize, 1); // Signal (start==start of window)
+    for (int i=0; i<fWindowSize; i++)
+    {
+        d[i][0] = fVal[i+maxp]/sum - fPulseShape[fWeightsPerBin*i + frac];
+        sumc += d[i][0]*d[i][0];
+    }
+
+/*
+    TMatrix Ainv; // Autocorrelation Matrix (inverted)
+
+    TMatrix dT(TMatrix::kTransposed, d);
+
+    TMatrix chisq = dT * (*fAinv*d);
+  */
+    return sumc;//chisq[0][0]/(fWindowSize-2);
+}
+
 #include <iostream>
-void MExtralgoDigitalFilter::Extract()
+void MExtralgoDigitalFilter::Extract(Int_t maxpos)
 {
     fSignal    =  0; // default is: no pulse found
@@ -54 +150 @@
     // Calculate the sum of the first fWindowSize slices
     //
-
-    // For the case of an even numberof weights/bin there is
+    // For the case of an even number of weights/bin there is
     // no central bin.So we create an artificial central bin.
     for (Int_t i=0; i<fNum-fWindowSize+1; i++)
@@ -66 +161 @@
         }
     }
+
+    /*
+     // This could be for a fast but less accurate extraction....
+     maxamp = Eval(fWeightsAmp, maxpos-fWindowSize/2);
+     maxp   = maxpos-fWindowSize/2;
+     */
+
     // The number of available slices were smaller than the
     // extraction window size of the extractor
@@ -136 +238 @@
         fTime -= timefineadjust;
 }
+
 
 #include <TH1.h>
@@ -145 +248 @@
 #include <TProfile.h>
 
-Bool_t MExtralgoDigitalFilter::CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb)
+Int_t MExtralgoDigitalFilter::CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb)
 {
     const Int_t weightsperbin = wpb<=0?shape.GetNbinsX()/autocorr.GetNbinsX():wpb;
@@ -154 +257 @@
         cout << "        Shape: " << shape.GetNbinsX() << endl;
         cout << "        ACorr: " << autocorr.GetNbinsX() << endl;
-        return kFALSE;
+        return -1;
     }
 
@@ -180 +283 @@
     // FIXME: DELETE!!!
     TH1 &derivative = *static_cast<TH1*>(shape.Clone());
+    derivative.SetDirectory(0);
     derivative.Reset();
 
@@ -233 +337 @@
         {
             cout << "ERROR - Division by zero: denom[0][0]==0 for i=" << i << "." << endl;
-            return kFALSE;
+            return -1;
         }
 
@@ -248 +352 @@
     }
 
-    return kTRUE;
+    return first*weightsperbin;
 }
 
-void MExtralgoDigitalFilter::CalculateWeights2(TH1F &shape, const TH2F &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb)
+Int_t MExtralgoDigitalFilter::CalculateWeights2(TH1 &shape, const TH2 &autocorr, TArrayF &weightsamp, TArrayF &weightstime, Int_t wpb)
 {
     const Int_t weightsperbin = wpb<=0?shape.GetNbinsX()/autocorr.GetNbinsX():wpb;
@@ -260 +364 @@
         cout << "        Shape: " << shape.GetNbinsX() << endl;
         cout << "        ACorr: " << autocorr.GetNbinsX() << endl;
-        return;
+        return -1;
     }
 
@@ -282 +386 @@
     TSpline5 val("Signal", &gr);
 
-    TH1F derivative(shape);
+    // FIXME: DELETE!!!
+    TH1 &derivative = *static_cast<TH1*>(shape.Clone());
+    derivative.SetDirectory(0);
     derivative.Reset();
 
@@ -331 +437 @@
         const TMatrixD denom  = (gT*(B*g))*(dT*(B*d)) - (dT*(B*g))*(dT*(B*g));
 
+        if (denom[0][0]==0)
+        {
+            cout << "ERROR - Division by zero: denom[0][0]==0 for i=" << i << "." << endl;
+            return -1;
+        }
+
         const TMatrixD w_amp  = (dT*(B*d))*(gT*B) - (gT*(B*d))*(dT*B);
         const TMatrixD w_time = (gT*(B*g))*(dT*B) - (gT*(B*d))*(gT*B);
@@ -342 +454 @@
         }
     }
+    return first*weightsperbin;
 }
-/*
-Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
-for (Int_t i = -fBinningResolutionHiGain/2+1; i<=fBinningResolutionHiGain/2; i++)
-{
-    // i = -4...5
-    for (Int_t count=0; count < fWindowSizeHiGain; count++)
-    {
-        // count=0: pos=10*(start+0) + 10 + i
-        // count=1: pos=10*(start+1) + 10 + i
-        // count=2: pos=10*(start+2) + 10 + i
-        // count=3: pos=10*(start+3) + 10 + i
-    }
-
-    for (Int_t count=0; count < fWindowSizeHiGain; count++)
-    {
-        // count=0: pos = 10*0 + 5 +i-1
-        // count=1: pos = 10*1 + 5 +i-1
-        // count=2: pos = 10*2 + 5 +i-1
-        // count=3: pos = 10*3 + 5 +i-1
-    }
-}
-
-Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionLoGain/2;
-for (Int_t i = -fBinningResolutionLoGain/2+1; i<=fBinningResolutionLoGain/2; i++)
-{
-    // i=-4..5
-    for (Int_t count=0; count < fWindowSizeLoGain; count++)
-    {
-        // count=0: pos = 10*(start+0) + 5 + i
-        // count=1: pos = 10*(start+1) + 5 + i
-        // count=2: pos = 10*(start+2) + 5 + i
-        // count=3: pos = 10*(start+3) + 5 + i
-    }
-
-    for (Int_t count=0; count < fWindowSizeLoGain; count++)
-    {
-        // count=0: pos = 10*0 + 5 +i-1
-        // count=1: pos = 10*1 + 5 +i-1
-        // count=2: pos = 10*2 + 5 +i-1
-        // count=3: pos = 10*3 + 5 +i-1
-    }
-}
-*/

trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h

@@ -2 +2 @@
 #define MARS_MExtralgoDigitalFilter
 
-#ifndef ROOT_TROOT
-#include <TROOT.h>
+#ifndef ROOT_TMatrix
+#include <TMatrix.h>
 #endif
 
@@ -12 +12 @@
 class TArrayF;
 
+//#include <TMatrix.h>
+
 class MExtralgoDigitalFilter
 {
 private:
     // Input
-    Float_t *fVal;
-    Int_t    fNum;
+    const Float_t *fVal;
+    Int_t fNum;
 
     Float_t const *fWeightsAmp;
     Float_t const *fWeightsTime;
+    Float_t const *fPulseShape;
+
+    const TMatrix *fAinv;
 
     const Int_t fWeightsPerBin; // Number of weights per data bin
@@ -30 +35 @@
     Float_t fSignal;
     Float_t fSignalDev;
+
+    Float_t GetChisq(const Int_t maxp, const Int_t frac, const Float_t sum) const;
+
+    inline Double_t ChiSq(const Double_t sum, const Int_t startv, const Int_t startw=0) const
+    {
+        //
+        // Slide with a window of size windowsize over the sample
+        // and multiply the entries with the corresponding weights
+        //
+        Double_t chisq = 0;
+
+        // Shift the start of the weight to the center of sample 0
+        Float_t const *w = fPulseShape + startw;
+
+        const Float_t *beg = fVal+startv;
+        for (Float_t const *pex=beg; pex<beg+fWindowSize; pex++)
+        {
+            const Double_t c = *w - *pex/sum;
+            chisq += c*c;
+            w += fWeightsPerBin;
+        }
+        return chisq;
+    }
 
     // Weights: Weights to evaluate
@@ -45 +73 @@
         Float_t const *w = weights + startw;
 
-        Float_t *const beg = fVal+startv;
+        const Float_t *beg = fVal+startv;
         for (Float_t const *pex=beg; pex<beg+fWindowSize; pex++)
         {
@@ -115 +143 @@
 
 public:
-    MExtralgoDigitalFilter(Int_t res, Int_t windowsize, Float_t *wa, Float_t *wt)
+    MExtralgoDigitalFilter(Int_t res, Int_t windowsize, Float_t *wa, Float_t *wt, Float_t *ps=0, TMatrix *ainv=0)
         : fVal(0), fNum(0), fWeightsAmp(wa+res/2), fWeightsTime(wt+res/2),
-        fWeightsPerBin(res), fWindowSize(windowsize),
+        fPulseShape(ps), fAinv(ainv), fWeightsPerBin(res), fWindowSize(windowsize),
         fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
     {
     }
 
-    void SetData(Int_t n, Float_t *val) { fNum=n; fVal=val; }
+    void SetData(Int_t n, Float_t const *val) { fNum=n; fVal=val; }
 
     Float_t GetTime() const          { return fTime; }
@@ -134 +162 @@
 
     Float_t ExtractNoise(Int_t iter) const;
-    void Extract();
+    void Extract(Int_t maxpos=-1);
 
-    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);
+    static Int_t CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
+    static Int_t CalculateWeights2(TH1 &shape, const TH2 &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
 };
 

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

@@ -207 +207 @@
 }
 
+#include <TRandom.h>
 Float_t MExtralgoSpline::ExtractNoise(Int_t iter)
 {
-    const Float_t nsx = iter * fResolution;
+    const Float_t nsx = gRandom->Uniform(); //iter * fResolution;
 
     if (fExtractionType == kAmplitude)
@@ -223 +224 @@
     fSignalDev = -1;
     fTimeDev   = -1;
-
+/*
     //
     // Allow no saturated slice and
@@ -229 +230 @@
     //
 
-    /*
     const Bool_t limlo = maxbin <      TMath::Ceil(fRiseTime);
     const Bool_t limup = maxbin > fNum-TMath::Ceil(fFallTime)-1;
@@ -248 +248 @@
         return;
     }
-    */
-
-    fTimeDev = fResolution;
 
     //
@@ -256 +253 @@
     //
 
-
-    /*
     Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
 
@@ -346 +341 @@
         }
     }
-  */
-
+
+    fTime      = maxpos;
+    fTimeDev   = fResolution;
+    fSignal    = CalcIntegral(maxpos);
+    fSignalDev = 0;  // means: is valid
+
+    return;
+*/
     // --- Start NEW ---
 
@@ -369 +370 @@
 
     Float_t maxpos, maxval;
+    // FIXME: Check the dfeault if no maximum found!!!
     GetMaxAroundI(maxbin, maxpos, maxval);
 
@@ -376 +378 @@
     {
         fTime      = maxpos;
+        fTimeDev   = fResolution;
         fSignal    = maxval;
         fSignalDev = 0;  // means: is valid
@@ -381 +384 @@
     }
 
-    //
-    // Loop from the beginning of the slice upwards to reach the maxhalf:
-    // With means of bisection:
-    //
-    /*
-    static const Float_t sqrt2 = TMath::Sqrt(2.);
-
-    step = sqrt2*3*0.061;//fRiseTime;
-    Float_t x = maxpos-0.86-3*0.061;//fRiseTime*1.25;
-
-//    step = sqrt2*0.5;//fRiseTime;
-//    Float_t x = maxpos-1.25;//fRiseTime*1.25;
-
-    Int_t  cnt  =0;
-    while (cnt++<30)
-    {
-        const Float_t y=EvalAt(x);
-
-        if (TMath::Abs(y-maxval/2)<fResolution)
-            break;
-
-        step /= sqrt2; // /2
-        x += y>maxval/2 ? -step : +step;
-    }
-    */
-
     // Search downwards for maxval/2
     // By doing also a search upwards we could extract the pulse width
@@ -412 +389 @@
 
     fTime      = x1;
+    fTimeDev   = fResolution;
     fSignal    = CalcIntegral(maxpos);
     fSignalDev = 0;  // means: is valid
 
-    //if (fSignal>100)
-    //    cout << "V=" << maxpos-x1 << endl;
+    //
+    // Loop from the beginning of the slice upwards to reach the maxhalf:
+    // With means of bisection:
+    //
+    /*
+    static const Float_t sqrt2 = TMath::Sqrt(2.);
+
+    step = sqrt2*3*0.061;//fRiseTime;
+    Float_t x = maxpos-0.86-3*0.061;//fRiseTime*1.25;
+
+//    step = sqrt2*0.5;//fRiseTime;
+//    Float_t x = maxpos-1.25;//fRiseTime*1.25;
+
+    Int_t  cnt  =0;
+    while (cnt++<30)
+    {
+        const Float_t y=EvalAt(x);
+
+        if (TMath::Abs(y-maxval/2)<fResolution)
+            break;
+
+        step /= sqrt2; // /2
+        x += y>maxval/2 ? -step : +step;
+    }
+    */
 
     //
@@ -422 +423 @@
     //
     //   fTime      = cnt==31 ? -1 : x;
+    //   fTimeDev   = fResolution;
     //   fSignal    = cnt==31 ? CalcIntegral(x) : CalcIntegral(maxpos);
     //   fSignalDev = 0;  // means: is valid

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

@@ -21 +21 @@
 
     // Input
-    Float_t *fVal;
-    Int_t    fNum;
+    Float_t const *fVal;
+    const Int_t    fNum;
 
     Float_t *fDer1;
@@ -238 +238 @@
 
         // Take a default in case no maximum is found
+        // FIXME: Check THIS!!!
         xmax=i;
         ymax=fVal[i];
@@ -438 +439 @@
 
 public:
-    MExtralgoSpline(Float_t *val, Int_t n, Float_t *der1, Float_t *der2)
+    MExtralgoSpline(const Float_t *val, Int_t n, Float_t *der1, Float_t *der2)
         : fExtractionType(kIntegral), fVal(val), fNum(n), fDer1(der1), fDer2(der2), fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
     {