Changeset 8154 for trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc (modified) (13 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
-    }
-}
-*/