Changeset 8154

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Makefile (modified) (1 diff)
• mextralgo/MExtralgoDigitalFilter.cc (modified) (13 diffs)
• mextralgo/MExtralgoDigitalFilter.h (modified) (6 diffs)
• mextralgo/MExtralgoSpline.cc (modified) (11 diffs)
• mextralgo/MExtralgoSpline.h (modified) (3 diffs)
• mjobs/MJCalibrateSignal.cc (modified) (5 diffs)
• mjobs/MJCalibration.cc (modified) (2 diffs)
• mjobs/MJPedestal.cc (modified) (5 diffs)
• mpedestal/MPedestalSubtract.cc (modified) (1 diff)
• mpedestal/MPedestalSubtract.h (added)
• mpedestal/MPedestalSubtractedEvt.cc (added)
• mpedestal/MPedestalSubtractedEvt.h (added)
• mpedestal/Makefile (modified) (1 diff)
• mpedestal/PedestalLinkDef.h (modified) (1 diff)
• msignal/MExtractTime.cc (modified) (5 diffs)
• msignal/MExtractTime.h (modified) (1 diff)
• msignal/MExtractTimeAndCharge.cc (modified) (10 diffs)
• msignal/MExtractTimeAndCharge.h (modified) (3 diffs)
• msignal/MExtractTimeAndChargeDigitalFilter.cc (modified) (32 diffs)
• msignal/MExtractTimeAndChargeDigitalFilter.h (modified) (1 diff)
• msignal/MExtractTimeAndChargeSpline.cc (modified) (11 diffs)
• msignal/MExtractTimeAndChargeSpline.h (modified) (1 diff)
• msignal/MExtractedSignalPix.h (modified) (2 diffs)
• msignal/MExtractor.cc (modified) (7 diffs)
• msignal/MExtractor.h (modified) (5 diffs)
• msignal/Makefile (modified) (1 diff)
• msignal/calibration_weights_UV46.dat (modified) (1 diff)
• msignal/cosmics_weights46.dat (modified) (1 diff)

trunk/MagicSoft/Mars/Makefile

@@ -45 +45 @@
           mhist \
           manalysis \
+          mextralgo \
           msignal \
           mbadpixels \

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)
     {

trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc

@@ -94 +94 @@
 #include "MFTriggerPattern.h"
 #include "MGeomApply.h"
+#include "MPedestalSubtract.h"
 //#include "MMcPedestalCopy.h"
 #include "MPointingPosCalc.h"
@@ -450 +451 @@
     MContinue conttp(&ftp, "ContTrigPattern");
 
+    // Create the pedestal subtracted raw-data
+    MPedestalSubtract pedsub;
+    pedsub.SetPedestalCam(&pedcamab);
+
     // Do signal and pedestal calculation
     MPedCalcFromLoGain     pedlo1("MPedCalcFundamental");
@@ -471 +476 @@
         extractor1->SetPedestals(&pedcamab);
 
+        // Setup to use the hi-gain extraction window in the lo-gain
+        // range (the start of the lo-gain range is added automatically
+        // by MPedCalcFromLoGain)
+        //
+        // The window size of the extractor is not yet initialized,
+        // so we have to stick to the extraction range
+        const Int_t f = extractor1->GetHiGainFirst();
+        const Int_t l = extractor1->GetHiGainLast();
+        const Int_t w = (l-f+1)&~1;
+
+        pedlo1.SetExtractWindow(f, w);
+
         if (extractor1->InheritsFrom("MExtractTimeAndCharge"))
         {
             pedlo2.SetExtractor((MExtractTimeAndCharge*)extractor1);
             pedlo3.SetExtractor((MExtractTimeAndCharge*)extractor1);
+            /*
             const Int_t win = ((MExtractTimeAndCharge*)extractor1)->GetWindowSizeHiGain();
             pedlo1.SetExtractWindow(15, win);
-            pedlo2.SetExtractWindow(15, win/*obsolete*/);
-            pedlo3.SetExtractWindow(15, win/*obsolete*/);
+            pedlo2.SetExtractWindow(15, win//obsolete//);
+            pedlo3.SetExtractWindow(15, win//obsolete//);
+            */
         }
         else
         {
+            /*
             // FIXME: How to get the fixed value 15 automatically?
             const Int_t f = (Int_t)(15.5+extractor1->GetHiGainFirst());
@@ -488 +508 @@
             pedlo2.SetExtractWindow(f, n);
             pedlo3.SetExtractWindow(f, n);
+            */
+            pedlo2.SetExtractWindow(f, w);
+            pedlo3.SetExtractWindow(f, w);
+
         }
     }
@@ -657 +681 @@
     tlist2.AddToList(&apply);
     tlist2.AddToList(&merge);
+    tlist2.AddToList(&pedsub);
     tlist2.AddToList(&pedlo1);
     tlist2.AddToList(&pedlo2);

trunk/MagicSoft/Mars/mjobs/MJCalibration.cc

@@ -154 +154 @@
 #include "MRawFileRead.h"
 #include "MGeomApply.h"
+#include "MPedestalSubtract.h"
 #include "MTaskEnv.h"
 #include "MBadPixelsMerge.h"
@@ -1907 +1908 @@
     tlist.AddToList(&apply);
 
-    MPedCalcPedRun           pedcalc;
-    pedcalc.SetExtractWindow(fExtractor->GetHiGainFirst(),TMath::Nint(fExtractor->GetNumHiGainSamples()));
+    // Produce pedestal subtracted raw-data
+    MPedestalSubtract pedsub;
+    pedsub.SetPedestalCam(fExtractor->GetPedestals());
+    tlist.AddToList(&pedsub);
+
+    // Setup to use the hi-gain extraction window in the lo-gain
+    // range (the start of the lo-gain range is added automatically
+    // by MPedCalcFromLoGain)
+    //
+    // The window size of the extractor is not yet initialized,
+    // so we have to stick to the extraction range
+    const Int_t f = fExtractor->GetHiGainFirst();
+    const Int_t l = fExtractor->GetHiGainLast();
+    const Int_t w = (l-f+1)&~1;
+
+    MPedCalcPedRun pedcalc;
+    pedcalc.SetExtractWindow(f, w);
 
     if (IsIntensity())

trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

@@ -82 +82 @@
 #include "MRawEvtData.h"
 #include "MGeomApply.h"
+#include "MPedestalSubtract.h"
 #include "MTriggerPatternDecode.h"
 #include "MBadPixelsMerge.h"
@@ -1064 +1065 @@
     // 
     if (!fBadPixelsFile.IsNull())
-      {
+    {
         *fLog << inf << "Excluding: " << fBadPixelsFile << endl;
-        ifstream fin(fBadPixelsFile.Data());
-        fBadPixels.AsciiRead((istream&)fin);
-      }
+        ifstream fin(fBadPixelsFile);
+        fBadPixels.AsciiRead(fin);
+    }
 
     MGeomApply geomapl;
@@ -1074 +1075 @@
 
     MPedCalcPedRun pedcalc;
-    pedcalc.SetPedestalUpdate(kFALSE);
+    //pedcalc.SetPedestalUpdate(kFALSE);
 
     MPedCalcFromLoGain pedlogain;
@@ -1132 +1133 @@
         tlist.AddToList(&fillpul);
     }
+
+    // produce pedestal subtracted raw-data
+    MPedestalSubtract pedsub;
+    if (fExtractor && fExtractionType!=kFundamental)
+        pedsub.SetPedestalCam(&fPedestalCamIn);
+    else
+        pedsub.SetNamePedestalCam(""); // Only copy hi- and lo-gain together!
+    tlist.AddToList(&pedsub);
 
     // ----------------------------------------------------------------------
@@ -1186 +1195 @@
     if (fExtractor)
     {
-      fExtractor->SetPedestals(&fPedestalCamIn);
-
-      if (fExtractionType!=kFundamental)
+        fExtractor->SetPedestals(&fPedestalCamIn);
+
+        if (fExtractionType!=kFundamental)
         {
-          pedcalc.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
-          pedlogain.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
-          
-          pedcalc.SetExtractor((MExtractTimeAndCharge*)fExtractor);
-          pedlogain.SetExtractor((MExtractTimeAndCharge*)fExtractor);
+            pedcalc.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
+            pedlogain.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
+
+            pedcalc.SetExtractor((MExtractTimeAndCharge*)fExtractor);
+            pedlogain.SetExtractor((MExtractTimeAndCharge*)fExtractor);
         }
-      
-      if (fExtractor->InheritsFrom("MExtractTimeAndCharge"))
+        else
         {
-
-          const Float_t f = 0.1+fExtractor->GetHiGainFirst();
-          const Int_t win = ((MExtractTimeAndCharge*)fExtractor)->GetWindowSizeHiGain();
-          pedcalc.SetExtractWindow((Int_t)f, win);
-          pedlogain.SetExtractWindow((Int_t)(15+f), win);
-
+            // The window size of the extractor is not yet initialized,
+            // so we have to stick to the extraction range
+            const Int_t f = fExtractor->GetHiGainFirst();
+            const Int_t l = fExtractor->GetHiGainLast();
+            const Int_t w = (l-f+1)&~1;
+
+            // Setup to use the hi-gain extraction window in the lo-gain
+            // range (the start of the lo-gain range is added automatically
+            // by MPedCalcFromLoGain)
+            pedcalc.SetExtractWindow(f, w);
+            pedlogain.SetExtractWindow(f, w);
         }
-      else
+
+        if (!fExtractor->InheritsFrom("MExtractTimeAndCharge") && fExtractionType!=kFundamental)
         {
-          const Float_t f = 0.1+fExtractor->GetHiGainFirst();
-          const Float_t n = 0.1+fExtractor->GetNumHiGainSamples();
-          pedcalc.SetExtractWindow((Int_t)f, (Int_t)n);
-          pedlogain.SetExtractWindow((Int_t)(15+f), (Int_t)n);
-
-          if (fExtractionType!=kFundamental)
-            {
-              *fLog << inf;
-              *fLog << "Signal extractor doesn't inherit from MExtractTimeAndCharge..." << endl;
-              *fLog << " --> falling back to fundamental pedestal extraction." << endl;
-              fExtractionType=kFundamental;
-            }
+            *fLog << inf;
+            *fLog << "Signal extractor doesn't inherit from MExtractTimeAndCharge..." << endl;
+            *fLog << " --> falling back to fundamental pedestal extraction." << endl;
+            fExtractionType=kFundamental;
         }
     }

trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.cc

@@ -1 +1 @@
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MPedestalSubtract.cc,v 1.2 2006-10-24 08:18:07 tbretz Exp $
+! --------------------------------------------------------------------------
 !
 ! *

trunk/MagicSoft/Mars/mpedestal/Makefile

@@ -42 +42 @@
            MHPedestalCor.cc \
            MPedestalSubtract.cc \
-           MPedestalSubtractedEvt.cc \
-           MPedCalcFromData.cc
+           MPedestalSubtractedEvt.cc
 
 ############################################################

trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h

@@ -19 +19 @@
 #pragma link C++ class MPedestalSubtract+;
 #pragma link C++ class MPedestalSubtractedEvt+;
-#pragma link C++ class MPedCalcFromData+;
 
 #pragma link C++ class MHPedestalCor+;

trunk/MagicSoft/Mars/msignal/MExtractTime.cc

@@ -1 +1 @@
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTime.cc,v 1.22 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
 !
 ! *
@@ -17 +19 @@
 !
 !   Author(s): Markus Gaug, 04/2004 <mailto:markus@ifae.es>
-!
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Author(s): Thomas Bretz, 04/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2006
 !
 !
@@ -137 +140 @@
 // - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
 //                                      fLoGainFirst, fLoGainLast, fNumLoGainSamples);
-//
+/*
 Bool_t MExtractTime::ReInit(MParList *pList)
 {
@@ -148 +151 @@
     return kTRUE;
 }
-
+*/
 // --------------------------------------------------------------------------
 //
 // Calculate the integral of the FADC time slices and store them as a new
 // pixel in the MArrivalTimeCam container.
-//
+/*
 Int_t MExtractTime::Process()
 {
@@ -183 +186 @@
       pix.SetGainSaturation(sathi, satlo);
  
-    } /* while (pixel.Next()) */
+    }
 
     fArrTime->SetReadyToSave();
 
     return kTRUE;
-}
+}*/
 
 void MExtractTime::Print(Option_t *o) const
 {
-    if (IsA()==MExtractTime::Class())
-        *fLog << GetDescriptor() << ":" << endl;
+//    if (IsA()==MExtractTime::Class())
+//        *fLog << GetDescriptor() << ":" << endl;
+    MExtractor::Print(o);
     *fLog << " Offset Lo-Gain:     " << fOffsetLoGain << endl;
-    MExtractor::Print(o);
-}
+}

trunk/MagicSoft/Mars/msignal/MExtractTime.h

@@ -17 +17 @@
   
   MArrivalTimeCam *fArrTime;          //! Container with the photons arrival times
-  
+/*
   virtual void FindTimeHiGain(Byte_t *firstused, Float_t &time, Float_t &dtime,
                               Byte_t &sat, const MPedestalPix &ped) const {}
   virtual void FindTimeLoGain(Byte_t *firstused, Float_t &time, Float_t &dtime,
                               Byte_t &sat, const MPedestalPix &ped) const {}
-
+  */
   Int_t  PreProcess( MParList *pList );
-  Bool_t ReInit    ( MParList *pList );
-  Int_t  Process   ();
+//  Bool_t ReInit    ( MParList *pList );
+//  Int_t  Process   ();
 
 public:

trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc

@@ -1 +1 @@
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.55 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
 !
 ! *
@@ -17 +19 @@
 !
 !   Author(s): Markus Gaug, 05/2004 <mailto:markus@ifae.es>
+!   Author(s): Thomas Bretz, 05/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2006
 !
 !
@@ -60 +63 @@
 //   MExtractedSignalCam
 //
+// For weired events check: Run 94127 Event 672, 1028
+//
 //////////////////////////////////////////////////////////////////////////////
 #include "MExtractTimeAndCharge.h"
@@ -70 +75 @@
 #include "MParList.h"
 
-#include "MRawEvtData.h"
+//#include "MArrayB.h"
+//#include "MArrayF.h"
+
+//#include "MRawEvtData.h"
 #include "MRawEvtPixelIter.h"
-#include "MRawRunHeader.h"
-
-#include "MPedestalCam.h"
-#include "MPedestalPix.h"
+//#include "MRawRunHeader.h"
+
+//#include "MPedestalCam.h"
+//#include "MPedestalPix.h"
+#include "MPedestalSubtractedEvt.h"
 
 #include "MArrivalTimeCam.h"
@@ -87 +96 @@
 using namespace std;
 
-const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -3.5; 
-const Byte_t  MExtractTimeAndCharge::fgLoGainSwitch = 120; 
+const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -6.0;
+const Byte_t  MExtractTimeAndCharge::fgLoGainSwitch     =  120;
 
 // --------------------------------------------------------------------------
@@ -95 +104 @@
 //
 // Sets: 
-// - fLoGainFirstSave to 0
 // - fWindowSizeHiGain and fWindowSizeLoGain to 0
 // - fLoGainStartShift to fgLoGainStartShift+fgOffsetLoGain
@@ -101 +109 @@
 //
 MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title)
-    : fLoGainFirstSave(0), fWindowSizeHiGain(0), fWindowSizeLoGain(0)
+    : /*fLoGainFirstSave(0),*/ fWindowSizeHiGain(0), fWindowSizeLoGain(0)
 {
     fName  = name  ? name  : "MExtractTimeAndCharge";
@@ -157 +165 @@
   
   if (fSignals)
-    fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
+    fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples,
                                 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
-
 
   return kTRUE;
@@ -175 +182 @@
 Int_t MExtractTimeAndCharge::Process()
 {
-
-  MRawEvtPixelIter pixel(fRawEvt);
-
-  while (pixel.Next())
-  {
-      // COPY HERE PRODUCING ARRAY WITH SAMPLES
-
-      /*
-       const MPedestalPix  &ped = (*fPedestals)[pixidx];
-
-       const Float_t pedes  = ped.GetPedestal();
-       const Float_t ABoffs = ped.GetPedestalABoffset();
-
-       const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
-
-       const Int_t num = pixel.GetNumHiGainSamples()+pixel.GetNumLoGainSamples();
-
-       MArrayF sample(num);
-
-       const Int_t abflag = pixel.HasABFlag() ? 1 : 0;
-       const Int_t ids0 = fHiGainFirst + abflag;
-
-       Int_t ids  = ids0;
-
-       Float_t  null   = 0;      // Starting value for maxpos
-       Float_t *maxpos = &null;  // Position of maximum
-       Float_t *sat1   = 0;      // First saturating slice
-       Float_t *sat2   = 0;      // Last saturating slice
-
-       const Float_t *beg = sample.GetArray();
-       const Float_t *end = beg + num;
-
-       Float_t *ptr = beg;
-       while (ptr<end)
-       {
-          *sample++ = (Float_t)*ptr - pedmean[ids++&0x1];
-
-          // This extraction must be done independant for lo- and hi-gain
-          // if (*ptr > *maxpos)
-          //     maxpos = ptr;
-          //
-          // if (*ptr >= fSaturationLimit)
-          // {
-          //    sat2 = ptr;
-          //    if (!sat1)
-          //       sat1 = ptr;
-          // }
-          //
-          // ptr++;
+    MRawEvtPixelIter pixel(fRawEvt);
+
+    while (pixel.Next())
+    {
+        const Int_t pixidx = pixel.GetPixelId();
+
+        //
+        // Preparations for the pedestal subtraction (with AB-noise correction)
+        //
+        Float_t *sig = fSignal->GetSamples(pixidx);
+
+        // ====================================================================
+        //    MOVE THIS TO MExtractTimeAndCharge
+        // ====================================================================
+        // number of hi- and lo-gains
+        const Int_t numh = pixel.GetNumHiGainSamples();
+        const Int_t numl = pixel.GetNumLoGainSamples();
+/*
+        // COPY HERE PRODUCING ARRAY WITH SAMPLES
+        static MArrayB sample;
+        sample.Set(numh+numl);
+
+        if (numl>0)
+        {
+            memcpy(sample.GetArray(),      pixel.GetHiGainSamples(), numh);
+            memcpy(sample.GetArray()+numh, pixel.GetLoGainSamples(), numl);
         }
-       */
-
-      //
-      // Find signal in hi- and lo-gain
-      //
-      Float_t sumhi =0., deltasumhi =0; // Set hi-gain of MExtractedSignalPix valid
-      Float_t timehi=0., deltatimehi=0; // Set hi-gain of MArrivalTimePix valid
-      Byte_t sathi=0;
-
-      // Initialize fMaxBinContent for the case, it gets not set by the derived class
-      fMaxBinContent = fLoGainSwitch + 1; 
-
-      const Int_t pixidx = pixel.GetPixelId();
-      const MPedestalPix  &ped = (*fPedestals)[pixidx];
-      const Bool_t higainabflag = pixel.HasABFlag();
-
-      FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), 
-                              sumhi, deltasumhi, timehi, deltatimehi,
-                              sathi, ped, higainabflag);
-
-      // Make sure that in cases the time couldn't be correctly determined
-      // more meaningfull default values are assigned
-      if (timehi<=0 || timehi>pixel.GetNumHiGainSamples())
-          timehi = gRandom->Uniform(pixel.GetNumHiGainSamples());
-      
-      Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
-      Float_t timelo=0., deltatimelo=-1;  // invalidate logain of MArrivalTimePix
-      Byte_t satlo=0;
-      
-      //
-      // Adapt the low-gain extraction range from the obtained high-gain time
-      //
-
-      // IN THIS CASE THE PIXEL SHOULD BE MARKED BAD!!!!
-      // MEANS: Hi gain has saturated, but it is too early to extract
-      // the lo-gain properly
-      // THIS produces pulse positions ~= -1
-      // The signal might be handled in MCalibrateData, but hwat's about
-      // the arrival times in MCalibrateRelTime
-      if (sathi && fMaxBinContent<=fLoGainSwitch)
-          deltasumlo=deltasumhi=deltatimelo=deltatimehi=-1;
-
-      // FIXME: What to do with the pixel if it saturates too early???
-      if (pixel.HasLoGain() && (fMaxBinContent > fLoGainSwitch /*|| sathi>0*/) )
-      {
-          fLoGainFirstSave = fLoGainFirst;
-
-          // sathi is the number (not index!) of the first saturating slice
-          // 0 indicates that no saturation was found
-          // FIMXME: Is 3 an accurate value?
-
-          const Float_t pos = sathi==0 ? timehi : (int)(sathi)-3;
-
-          if (pos+fLoGainStartShift>0)
-              fLoGainFirst = (Byte_t)(pos + fLoGainStartShift);
-
-          if (fLoGainFirst<fLoGainFirstSave)
-              fLoGainFirst = fLoGainFirstSave;
-
-          if (fLoGainLast-fLoGainFirst >= fWindowSizeLoGain)
-          {
-              deltasumlo  = 0; // make logain of MExtractedSignalPix valid
-              deltatimelo = 0; // make logain of MArrivalTimePix valid
-
-              const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
-              FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
-                                      sumlo, deltasumlo, timelo, deltatimelo,
-                                      satlo, ped, logainabflag);
-
-              if (satlo>6)
-                  deltasumlo=deltatimelo=-1;
-          }
-          fLoGainFirst = fLoGainFirstSave;
-
-          // Make sure that in cases the time couldn't be correctly determined
-          // more meaningfull default values are assigned
-          if (timelo<=0 || timelo>pixel.GetNumLoGainSamples())
-              timelo = gRandom->Uniform(pixel.GetNumLoGainSamples());
-      }
-
-      MExtractedSignalPix &pix = (*fSignals)[pixidx];
-      MArrivalTimePix     &tix = (*fArrTime)[pixidx];
-      pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo);
-      pix.SetGainSaturation(sathi, satlo);
-
-      tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
-      tix.SetGainSaturation(sathi, satlo);
- 
+
+        // get pedestal information
+        const MPedestalPix &pedpix = (*fPedestals)[pixidx];
+
+        // pedestal information
+        const Int_t   ab  = pixel.HasABFlag() ? 1 : 0;
+        const Float_t ped = pedpix.GetPedestal();
+
+        // destination array
+        static MArrayF sig;
+        sig.Set(numh+numl);
+
+        // start of destination array, end of hi-gain destination array
+        // and start of hi-gain samples
+        Float_t *beg = sig.GetArray();
+        Float_t *end = beg + sig.GetSize();
+
+        // determine with which pedestal (+/- AB offset) to start
+        const Bool_t  noswap  = (ab&1)==((beg-beg)&1);
+        const Float_t offh    = noswap ? pedpix.GetPedestalABoffset() : -pedpix.GetPedestalABoffset();
+        const Float_t mean[2] = { ped + offh, ped - offh };
+
+        //const Int_t rangeh = fHiGainLast - fHiGainFirst + 1  + fHiLoLast;
+        //const Int_t rangel = fLoGainLast - fLoGainFirst + 1;
+
+        const Byte_t *src = numl>0 ? sample.GetArray() : pixel.GetHiGainSamples();
+        //const Byte_t *src = sample.GetArray();
+
+        // Copy hi-gains into array and substract pedestal
+        // FIXME: Shell we really subtract the pedestal from saturating slices???
+        for (Float_t *ptr=beg; ptr<end; ptr++)
+            *ptr = (Float_t)*src++ - mean[(ptr-beg)&1];
+
+        // Determin saturation of hi-gains
+        Byte_t *p0 = fSignal->GetSamplesRaw(pixidx); //numl>0 ? sample.GetArray() : pixel.GetHiGainSamples();
+//        Byte_t *p0 = sample.GetArray();
+
+        Byte_t maxcont  =  0;
+        Int_t  maxposhi = -1;
+
+        Byte_t *sathi0 = 0; // first saturating hi-gain slice
+        Byte_t *sathi1 = 0; // last saturating hi-gain slice
+        for (Byte_t *ptr=p0+fHiGainFirst; ptr<p0+fHiGainLast+fHiLoLast+1; ptr++)
+        {
+            // Get hi-gain maxbincontent
+            // Put this into its own function and loop?
+            if (*ptr>maxcont)
+            {
+                // ORGINAL:
+                //Int_t range = (fHiGainLast - fHiGainFirst + 1 + fHiLoLast) - fWindowSizeHiGain + 1;
+                // maxpos>1 && maxpos<=range
+
+                // This is a strange workaround to fit the previous
+                // Spline setup: TO BE CHECKED!
+                const Int_t pos = ptr-p0;
+                if (pos<fHiGainLast+1)
+                {
+                    maxposhi = pos-fHiGainFirst;
+
+                    if (maxposhi>1 && maxposhi<fHiGainLast-fHiGainFirst*//*+1 -fWindowSizeHiGain+1*//*+fHiLoLast*//*)
+                        maxcont = *ptr;
+                }
+            }
+
+            // FIXME: Do we also have to really COUNT the number
+            // of saturating slices, eg. for crosschecks?
+            if (*ptr>=fSaturationLimit)
+            {
+                sathi1 = ptr;
+                if (!sathi0)
+                    sathi0 = ptr;
+            }
+        }
+*/
+        Int_t sathi0 = fHiGainFirst;          // First slice to extract and first saturating slice
+        Int_t sathi1 = fHiGainLast/*+fHiLoLast*/; // Last  slice to extract and last saturating slice
+
+        Int_t maxcont;
+        Int_t maxposhi = fSignal->GetMax(pixidx, sathi0, sathi1, maxcont);
+        // Would it be better to take lastsat-firstsat?
+        Int_t numsathi = fSignal->GetSaturation(pixidx, fSaturationLimit, sathi0, sathi1);
+/*
+        // sathi2 is the number (not the index) of first saturating slice
+//        const Byte_t sathi2   = sathi0<0 ? 0 : sathi0+1;
+
+        // Number (not index) of first saturating slice
+ //       const Byte_t sathi2 = sathi0==0 ? 0 : sathi0-p0+1;
+//        const Byte_t numsathi = sathi0==0 ? 0 : sathi1-sathi0+1;
+
+//        Int_t maxb = maxcont;
+
+        // ====================================================================
+        // FIXME: Move range out of the extractors...
+
+        // Initialize maxcont for the case, it gets not set by the derived class
+        Float_t sumhi2 =0., deltasumhi2 =-1; // Set hi-gain of MExtractedSignalPix valid
+        Float_t timehi2=0., deltatimehi2=-1; // Set hi-gain of MArrivalTimePix valid
+*/
+        Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
+        Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
+        if (numsathi<1)
+        {
+
+        const Int_t rangehi = fHiGainLast - fHiGainFirst + 1/* + fHiLoLast*/;
+        FindTimeAndChargeHiGain2(sig/*.GetArray()*/+fHiGainFirst, rangehi,
+                                 sumhi, deltasumhi, timehi, deltatimehi,
+                                 numsathi, maxposhi);
+
+        timehi += fHiGainFirst;
+        }
+/*
+        Float_t sumhi = sumhi2;
+        Float_t deltasumhi = deltasumhi2;
+        Float_t timehi = timehi2;
+        Float_t deltatimehi=deltatimehi2;
+//        Byte_t sathi = sathi2;
+
+
+        //
+        // Find signal in hi- and lo-gain
+        //
+        Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
+        Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
+        Byte_t sathi=0;
+
+        // Initialize maxcont for the case, it gets not set by the derived class
+        maxcont = fLoGainSwitch + 1;
+
+        //const Int_t pixidx = pixel.GetPixelId();
+        //const MPedestalPix  &ped = (*fPedestals)[pixidx];
+        const Bool_t higainabflag = pixel.HasABFlag();
+        FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(),
+                                sumhi, deltasumhi, timehi, deltatimehi,
+                                sathi, pedpix, higainabflag);
+        timehi += fHiGainFirst;
+
+  */
+        // Make sure that in cases the time couldn't be correctly determined
+        // more meaningfull default values are assigned
+        if (timehi<=fHiGainFirst || timehi>=fHiGainLast)
+            timehi = gRandom->Uniform(fHiGainLast-fHiGainFirst)+fHiGainFirst;
+
+        Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
+        Float_t timelo=0., deltatimelo=-1;  // invalidate logain of MArrivalTimePix
+        Byte_t satlo=0;
+        Int_t numsatlo=0;
+
+        //
+        // Adapt the low-gain extraction range from the obtained high-gain time
+        //
+
+        // IN THIS CASE THE PIXEL SHOULD BE MARKED BAD!!!!
+        // MEANS: Hi gain has saturated, but the signal is to dim
+        // to extract the lo-gain properly
+        // THIS produces pulse positions ~= -1
+        // The signal might be handled in MCalibrateData, but hwat's about
+        // the arrival times in MCalibrateRelTime
+        if (numsathi>0 && maxcont<=fLoGainSwitch)
+            deltasumlo=deltasumhi=deltatimelo=deltatimehi=-1;
+
+        // If more than 8 hi-gain slices have saturated this is
+        // no physical event. We just assume that something with
+        // the extraction is wrong
+        if (numsathi>8) // FIXME: Should be something like 2?
+            deltasumhi=deltatimehi=-1;
+
+        // FIXME: What to do with the pixel if it saturates too early???
+        if (pixel.HasLoGain() && (maxcont > fLoGainSwitch /*|| sathi>0*/) )
+        {
+            // To determin the window in which the lo-gain is extracted
+            // clearly more information about the relation between the
+            // extraction window and the reslting time is necessary.
+            //fLoGainStartShift = -6.0;
+
+            const Byte_t fLoGainFirstSave = fLoGainFirst;
+
+            // sathi is the number (not index!) of the first saturating slice
+            // 0 indicates that no saturation was found
+            // FIXME: Is 0.5 should be expressed by the rise time of
+            //        the hi-gain signal!
+            const Float_t pos = numsathi==0 ? timehi : sathi0-0.5;
+
+            const Int_t lostart = TMath::FloorNint(pos+6);
+
+            const Int_t newfirst = TMath::FloorNint(pos+fLoGainStartShift);
+            fLoGainFirst = newfirst>fLoGainFirstSave ? newfirst : fLoGainFirstSave;
+
+//            if (fLoGainLast-fLoGainFirst >= fWindowSizeLoGain)
+            {
+/*
+                Float_t sumlo2 =0., deltasumlo2 =-1.; // invalidate logain of MExtractedSignalPix
+                Float_t timelo2=0., deltatimelo2=-1.;  // invalidate logain of MArrivalTimePix
+
+                // Determin saturation in lo-gains
+                Byte_t *satlo0 = 0; // first saturating hi-gain slice
+                Byte_t *satlo1 = 0; // last saturating hi-gain slice
+                Int_t maxposlo = -1;
+                Byte_t maxlo = 0;
+
+                Byte_t *p0 = fSignal->GetSamplesRaw(pixidx);
+                for (Byte_t *ptr=p0+numh+fLoGainFirst; ptr<p0+numh+fLoGainLast+1; ptr++)
+                {
+                    if (*ptr>maxlo)
+                    {
+                        // This is a strange workaround to fit the previous
+                        // Spline setup: TO BE CHECKED!
+                        const Int_t ipos = ptr-p0-numh;
+                        if (ipos>=fLoGainFirst && ipos<=fLoGainLast)
+                        {
+                            maxposlo = ipos-fLoGainFirst;
+                            maxlo = *ptr;
+                        }
+                    }
+                    if (*ptr>=fSaturationLimit)
+                    {
+                        satlo1 = ptr;
+                        if (!satlo0)
+                            satlo0 = ptr;
+                    }
+                    // FIXME: Do we also have to really COUNT the number
+                    // of saturating slices, eg. for crosschecks?
+                }
+
+                // Number of saturating lo-gain slices (this is not necessary
+                // identical to a counter counting the number of saturating
+                // slices!)
+                const Byte_t numsatlo = satlo0==0 ? 0 : satlo1-satlo0+1;
+*/
+
+                Int_t satlo0 = numh+fLoGainFirst;   // First slice to extract and first saturating slice
+                Int_t satlo1 = numh+fLoGainLast;    // Last  slice to extract and last saturating slice
+
+                // Would it be better to take lastsat-firstsat?
+                Int_t maxlo;
+                Int_t maxposlo = fSignal->GetMax(pixidx, satlo0, satlo1, maxlo);
+                numsatlo = fSignal->GetSaturation(pixidx, fSaturationLimit, satlo0, satlo1);
+
+                //                const Byte_t satlo2 = numsatlo;
+
+                const Int_t rangelo = fLoGainLast - fLoGainFirst + 1;
+                FindTimeAndChargeLoGain2(sig/*.GetArray()*/+numh+fLoGainFirst, rangelo,
+                                         sumlo, deltasumlo, timelo, deltatimelo,
+                                         numsatlo, maxposlo);
+                timelo += fLoGainFirst;
+
+//                sumlo =sumlo2, deltasumlo=deltasumlo2; // invalidate logain of MExtractedSignalPix
+//                timelo=timelo2, deltatimelo=deltatimelo2;  // invalidate logain of MArrivalTimePix
+//                satlo = satlo2;
+
+                /*
+                // THIS IS NOW THE JOB OF THE EXTRACTOR!
+                //deltasumlo  = 0; // make logain of MExtractedSignalPix valid
+                //deltatimelo = 0; // make logain of MArrivalTimePix valid
+
+                const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
+                FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
+                                        sumlo, deltasumlo, timelo, deltatimelo,
+                                        satlo, pedpix, logainabflag);
+                timelo += fLoGainFirst;
+                */
+
+
+                // If more than 6 lo-gain slices have saturated this is
+                // no physical event. We just assume that something with
+                // the extraction is wrong
+                if (numsatlo>6)
+                    deltasumlo=deltatimelo=-1;
+            }
+            fLoGainFirst = fLoGainFirstSave;
+
+            // Make sure that in cases the time couldn't be correctly determined
+            // more meaningfull default values are assigned
+            if (timelo<=fLoGainFirst || timelo>=fLoGainLast)
+                timelo = gRandom->Uniform(fLoGainLast-fLoGainFirst)+fLoGainFirst;
+       }
+
+        MExtractedSignalPix &pix = (*fSignals)[pixidx];
+        MArrivalTimePix     &tix = (*fArrTime)[pixidx];
+        pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo);
+        pix.SetGainSaturation(numsathi, numsatlo);
+//        pix.SetGainSaturation(sathi, satlo);
+
+        tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
+        tix.SetGainSaturation(numsathi, numsatlo);
+//        tix.SetGainSaturation(sathi, satlo);
     } /* while (pixel.Next()) */
 
-  fArrTime->SetReadyToSave();
-  fSignals->SetReadyToSave();
-
-  return kTRUE;
+    fArrTime->SetReadyToSave();
+    fSignals->SetReadyToSave();
+
+    return kTRUE;
 }
 
@@ -352 +533 @@
 void MExtractTimeAndCharge::Print(Option_t *o) const
 {
-  if (IsA()==Class())
-    *fLog << GetDescriptor() << ":" << endl;
+  MExtractTime::Print(o);
+//  if (IsA()==Class())
+//    *fLog << GetDescriptor() << ":" << endl;
 
   *fLog << dec;
-  *fLog << " Taking " << fWindowSizeHiGain
-        << " HiGain samples from slice " << (Int_t)fHiGainFirst
-        << " to " << (Int_t)(fHiGainLast+fHiLoLast) << " incl" << endl;
-  *fLog << " Taking " << fWindowSizeLoGain
-        << " LoGain samples from slice " << (Int_t)fLoGainFirst
-        << " to " << (Int_t)fLoGainLast << " incl" << endl;
-  
+  //*fLog << " Taking " << fWindowSizeHiGain
+  //      << " HiGain samples from slice " << (Int_t)fHiGainFirst
+  //      << " to " << (Int_t)(fHiGainLast/*+fHiLoLast*/) << " incl" << endl;
+  //*fLog << " Taking " << fWindowSizeLoGain
+  //      << " LoGain samples from slice " << (Int_t)fLoGainFirst
+  //      << " to " << (Int_t)fLoGainLast << " incl" << endl;
+  //
   *fLog << " LoGainStartShift:   " << fLoGainStartShift << endl;
   *fLog << " LoGainSwitch:       " << (Int_t)fLoGainSwitch << endl;
-  MExtractTime::Print(o);
 }

trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h

@@ -14 +14 @@
   static const Byte_t  fgLoGainSwitch;     //! Default for fLoGainSwitch     (now set to: 100)
   
-protected:
+private:
 
-  Byte_t  fLoGainFirstSave;       //! Temporary variable to store the original position of low-gain start slice
+//  Byte_t  fLoGainFirstSave;       //! Temporary variable to store the original position of low-gain start slice
   Float_t fLoGainStartShift;      // Shift to start searching the low-gain signal obtained from the high-gain times.
 
   Byte_t  fLoGainSwitch;          // Limit for max. bin content before the low-gain gets extracted
 
+protected:
+
   Int_t  fWindowSizeHiGain;       //  Window Size High-Gain
   Int_t  fWindowSizeLoGain;       //  Window Size Low-Gain
 
-  Byte_t fMaxBinContent;          //  Maximum bin content
-
   Int_t  PreProcess(MParList *pList);
   Int_t  Process();
-  Bool_t ReInit(MParList *pList);
-
   Int_t  ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 
@@ -47 +45 @@
   virtual void SetWindowSize(Int_t windowh, Int_t windowl) { fWindowSizeHiGain = windowh;
                                                            fWindowSizeLoGain = windowl;  }
+
+  Bool_t ReInit(MParList *pList);
   
   virtual Bool_t InitArrays() { return kTRUE; }
-
+/*
   virtual void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum,
                                        Float_t &time, Float_t &dtime,
@@ -57 +57 @@
                                        Float_t &time, Float_t &dtime,
                                        Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { } 
+  */
+  virtual void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum,
+                                       Float_t &time, Float_t &dtime,
+                                       Byte_t sat, Int_t maxpos) { }
+  
+  virtual void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum,  Float_t &dsum,
+                                       Float_t &time, Float_t &dtime,
+                                       Byte_t sat, Int_t maxpos) { }
+
+  // For MExtractPedestal
 
   ClassDef(MExtractTimeAndCharge, 2)   // Time And Charge Extractor Base Class

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc

@@ -1 +1 @@
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.73 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
 !
 ! *
@@ -72 +74 @@
 #include "MCalibrationPattern.h"
 #include "MExtractedSignalCam.h"
+#include "MExtralgoDigitalFilter.h"
 
 #include "MPedestalPix.h"
@@ -83 +86 @@
 const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst             =  1;
 const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast              = 14;
-const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain        =  4;
-const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain        =  6;
+//const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain        =  4;
+//const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain        =  6;
 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
@@ -90 +93 @@
 const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain    =  4;
 const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain           =  0.95;
-const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift       = -1.8;
+//const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift       = -1.8;
 
 // --------------------------------------------------------------------------
@@ -104 +107 @@
 //
 MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title) 
-    : fTimeShiftHiGain(0.), fTimeShiftLoGain(0.), fAutomaticWeights(kTRUE), fRandomIter(0)
+    : fBinningResolutionHiGain(fgBinningResolutionHiGain),
+    fBinningResolutionLoGain(fgBinningResolutionLoGain),
+    fAutomaticWeights(kTRUE), fRandomIter(0)
 {
     fName  = name  ? name  : "MExtractTimeAndChargeDigitalFilter";
@@ -110 +115 @@
 
     SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
-    SetWindowSize();
-    SetBinningResolution();
-    SetSignalStartBin();
+    SetWindowSize(3, 5);
+//    SetBinningResolution();
+//    SetSignalStartBin();
 
     SetOffsetLoGain(fgOffsetLoGain);          // Order between both
-    SetLoGainStartShift(fgLoGainStartShift);  // is important
+//    SetLoGainStartShift(fgLoGainStartShift);  // is important
 }
 
@@ -121 +126 @@
 //
 // Checks:
-// - if a window is bigger than the one defined by the ranges, set it to the available range
+// - if a window is bigger than the one defined by the ranges, set it
+// to the available range
 // 
 // Sets:
 // - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
 // - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
-//  
+//
+// This function might be used to turn the digital filter into a
+// sliding window extractor by setting the filename to NULL
+//
 void MExtractTimeAndChargeDigitalFilter::SetWindowSize(Int_t windowh, Int_t windowl)
 {
-
-  if (windowh != fgWindowSizeHiGain)
-    *fLog << warn << GetDescriptor() 
-          << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default." << endl;
-  if (windowl != fgWindowSizeLoGain)
-    *fLog << warn << GetDescriptor() 
-          << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default" << endl;
-
   fWindowSizeHiGain = windowh;
   fWindowSizeLoGain = windowl;
@@ -184 +185 @@
   // size of 6. The exact numbers have to be found still.
   //
-  fNumHiGainSamples = fAmpWeightsHiGain.GetSum()/fBinningResolutionHiGain;
-  fNumLoGainSamples = fAmpWeightsLoGain.GetSum()/fBinningResolutionLoGain;
+  fNumHiGainSamples  = fWindowSizeHiGain;
+  fNumLoGainSamples  = fWindowSizeLoGain;
   fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
   fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
@@ -191 +192 @@
 }
 
+
 // --------------------------------------------------------------------------
 //
@@ -213 +215 @@
 //
 // If new weights are set
-//  fTimeShiftHiGain
-//  fTimeShiftLoGain
 //  fNumHiGainSamples
 //  fNumLoGainSamples
@@ -234 +234 @@
 
     //
-    // shift the times back to the right reference (start counting from 0)
+    // We need here the effective number of samples. In pricipal the number
+    // is different depending on the weights used and must be set
+    // event by event.
     //
-    // The high-gain has the extraction offset (fHiGainFirst) already included
-    // here for speed reason. The low-gain has a changing extraction offset,
-    // so it will be added at every event (in FindTimeAndChargeLoGain)
-    fTimeShiftHiGain = 0.5 + 1./fBinningResolutionHiGain;
-    fTimeShiftLoGain = 0.5 + 1./fBinningResolutionLoGain;
-
-    //
-    // We need here the effective number of samples which is about 2.5 in the case of a window
-    // size of 6. The exact numbers have to be found still.
-    //
-    fNumHiGainSamples = (Float_t)fWindowSizeHiGain/2.4;
-    fNumLoGainSamples = (Float_t)fWindowSizeLoGain/2.4;
+    fNumHiGainSamples  = fAmpWeightsHiGain.GetSum()/fBinningResolutionHiGain;
+    fNumLoGainSamples  = fAmpWeightsLoGain.GetSum()/fBinningResolutionLoGain;
     fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
     fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
@@ -253 +245 @@
     // From MExtractTimeAndCharge::ReInit
     if (fSignals)
-        fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
+        fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples,
                                     fLoGainFirst, fLoGainLast, fNumLoGainSamples);
     return kTRUE;
@@ -269 +261 @@
         return kFALSE;
 
-    const Int_t rangehi = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
-    const Int_t rangelo = (Int_t)(fLoGainLast - fLoGainFirst + 1);
-
-    fHiGainSignal.Set(rangehi);
-    fLoGainSignal.Set(rangelo);
+    //const Int_t rangehi = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
+    //const Int_t rangelo = (Int_t)(fLoGainLast - fLoGainFirst + 1);
+
+    //cout << "ARRAYS INITIALIZED TO : " << rangehi << " " << rangelo << endl;
+
+    //fHiGainSignal.Set(rangehi);
+    //fLoGainSignal.Set(rangelo);
 
     return GetWeights();
@@ -295 +289 @@
 }
 
-void MExtractTimeAndChargeDigitalFilter::CalcBinningResArrays()
-{
-
-  fArrBinningResHiGain.Set(fWindowSizeHiGain);
-  fArrBinningResHalfHiGain.Set(fWindowSizeHiGain);
-
-  for (int i=0; i<fWindowSizeHiGain; i++)
-    {
-      fArrBinningResHiGain[i] = fBinningResolutionHiGain*i;
-      fArrBinningResHalfHiGain[i] = fArrBinningResHiGain[i] + fBinningResolutionHalfHiGain;
-    }
-
-  fArrBinningResLoGain.Set(fWindowSizeLoGain);
-  fArrBinningResHalfLoGain.Set(fWindowSizeLoGain);
-  
-  for (int i=0; i<fWindowSizeLoGain; i++)
-    {
-      fArrBinningResLoGain[i] = fBinningResolutionLoGain*i;
-      fArrBinningResHalfLoGain[i] = fArrBinningResLoGain[i] + fBinningResolutionHalfLoGain;
-    }
-}
-
 // --------------------------------------------------------------------------
 //
 // Apply the digital filter algorithm to the high-gain slices.
 //
-void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Float_t &dsum, 
+void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain2(const Float_t *ptr, Int_t num,
+                                                                  Float_t &sum, Float_t &dsum,
+                                                                  Float_t &time, Float_t &dtime,
+                                                                  Byte_t sat, Int_t maxpos)
+{
+    // Do some handling if maxpos is last slice!
+
+    MExtralgoDigitalFilter df(fBinningResolutionHiGain, fWindowSizeHiGain,
+                              fAmpWeightsHiGain.GetArray(),
+                              fTimeWeightsHiGain.GetArray(),
+                              fPulseHiGain.GetArray());
+    df.SetData(num, ptr);
+
+    if (IsNoiseCalculation())
+    {
+        if (fRandomIter == fBinningResolutionHiGain)
+            fRandomIter = 0;
+
+        sum = df.ExtractNoise(fRandomIter);
+        fRandomIter++;
+        return;
+    }
+
+    df.Extract(/*maxpos*/);
+    df.GetSignal(sum, dsum);
+    df.GetTime(time, dtime);
+}
+
+void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain2(const Float_t *ptr, Int_t num,
+                                                                  Float_t &sum, Float_t &dsum,
+                                                                  Float_t &time, Float_t &dtime,
+                                                                  Byte_t sat, Int_t maxpos)
+{
+    MExtralgoDigitalFilter df(fBinningResolutionLoGain, fWindowSizeLoGain,
+                              fAmpWeightsLoGain.GetArray(),
+                              fTimeWeightsLoGain.GetArray(),
+                              fPulseLoGain.GetArray());
+
+    df.SetData(num, ptr);
+
+    if (IsNoiseCalculation())
+    {
+        if (fRandomIter == fBinningResolutionLoGain)
+            fRandomIter = 0;
+
+        sum = df.ExtractNoise(fRandomIter);
+        return;
+    }
+
+    df.Extract(/*maxpos*/);
+    df.GetSignal(sum, dsum);
+    df.GetTime(time, dtime);
+}
+
+/*
+// --------------------------------------------------------------------------
+//
+// Apply the digital filter algorithm to the high-gain slices.
+//
+void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
                                                                  Float_t &time, Float_t &dtime,
                                                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
 {
-
-  Int_t range = fHiGainLast - fHiGainFirst + 1;
-  
-  const Byte_t *end = ptr + range;
-  Byte_t *p     = ptr;
-
-  //
-  // Preparations for the pedestal subtraction (with AB-noise correction)
-  //
-  const Float_t pedes  = ped.GetPedestal();
-  const Float_t ABoffs = ped.GetPedestalABoffset();
-
-  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
-
-  range += fHiLoLast;
-  fMaxBinContent = 0;
-  //
-  // Check for saturation in all other slices
-  //
-  Int_t ids = fHiGainFirst;
-  Float_t *sample = fHiGainSignal.GetArray();
-  while (p<end)
-    {
-      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
-
-      if (*p > fMaxBinContent)
-        {
-          Byte_t maxpos = p-ptr;
-
-          // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
-          if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
-            fMaxBinContent = *p;
-        }
-
-      if (*p++ >= fSaturationLimit)
-        if (!sat)
-          sat = ids-fHiGainFirst;
-    }
-
-  //
-  // Slide with a window of size fWindowSizeHiGain over the sample 
-  // and multiply the entries with the corresponding weights
-  //
-  if (IsNoiseCalculation())
-    {
-      if (fRandomIter == fBinningResolutionHiGain)
-        fRandomIter = 0;
-      for (Int_t ids=0; ids < fWindowSizeHiGain; ids++)
+    Int_t range = fHiGainLast - fHiGainFirst + 1;
+    const Byte_t *end = first + range;
+    Byte_t       *p  = first;
+
+    const Float_t pedes  = ped.GetPedestal();
+    const Float_t ABoffs = ped.GetPedestalABoffset();
+
+    const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
+
+    Byte_t maxcont = 0;
+    Int_t  maxpos  = 0;  // obsolete for Digital Filter (used in spline!)
+    sat = 0;
+
+    //
+    // Check for saturation in all other slices
+    //
+    Int_t ids = fHiGainFirst;
+    Float_t *sample = fHiGainSignal.GetArray();
+
+    while (p<end)
+    {
+        *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
+
+        if (*p > maxcont)
         {
-          const Int_t   idx = fArrBinningResHiGain[ids] + fRandomIter;
-          sum              += fAmpWeightsHiGain [idx]*fHiGainSignal[ids]; 
+            maxpos = ids-fHiGainFirst-1;
+            // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
+            //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
+            if (maxpos > 1 && maxpos < (range-1*//* - fWindowSizeHiGain + 1*//*))
+                maxcont = *p;
         }
-      fRandomIter++;
-      return;
-    }
-  
-  if (fHiLoLast != 0)
-    {
-
-      end = logain + fHiLoLast;
-
-      while (logain<end)
+
+        if (*p++ >= fSaturationLimit)
+            if (!sat)
+                sat = ids-fHiGainFirst;
+    }
+
+    if (fHiLoLast != 0)
+    {
+        end = logain + fHiLoLast;
+
+        while (logain<end)
         {
-
-          *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
-
-          if (*logain++ >= fSaturationLimit)
-            if (!sat)
-              sat = ids-fHiGainFirst;
+            *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
+
+            if (*logain > maxcont)
+            {
+                maxpos = ids-fHiGainFirst-1;
+
+                // FIXME!!! MUST BE THERE!
+                // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
+                //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
+                //    maxcont = *logain;
+            }
+
+            if (*logain++ >= fSaturationLimit)
+                if (!sat)
+                    sat = ids-fHiGainFirst;
+
+            range++;
         }
     }
 
-  //
-  // allow no saturated slice 
-  //
-  //  if (sat > 0)
-  //    return;
-  
-  Float_t time_sum  = 0.;
-  Float_t fmax      = -FLT_MAX;
-  Float_t ftime_max = 0.;
-  Int_t   max_p     = -1;
-  
-  //
-  // Calculate the sum of the first fWindowSize slices
-  //
-  for (Int_t i=0;i<range-fWindowSizeHiGain+1;i++)
-    {
-      sum      = 0.;
-      time_sum = 0.;
-
-      //
-      // Slide with a window of size fWindowSizeHiGain over the sample
-      // and multiply the entries with the corresponding weights
-      //
-      for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
-      {
-        const Int_t   idx = fBinningResolutionHiGain*sample+fBinningResolutionHalfHiGain;
-        const Float_t pex = fHiGainSignal[sample+i];
-        sum              += fAmpWeightsHiGain [idx]*pex;
-        time_sum         += fTimeWeightsHiGain[idx]*pex;
-      }
-
-      if (sum>fmax)
-      {
-        fmax      = sum;
-        ftime_max = time_sum;
-        max_p     = i;
-      }
-    } /*   for (Int_t i=0;i<range-fWindowSizeHiGain+1;i++) */
-
-  if (max_p<0)
-  {
-      sum = 0;
-      time = -1;
-      dsum = -1;
-      dtime = -1;
-      return;
-  }
-
-  if (fmax==0)
-  {
-      sum = 0;
-      time = -1;
-      dsum = 0;
-      dtime = 0;
-      return;
-  }
-
-  ftime_max        /= fmax;
-  Int_t t_iter      = Int_t(ftime_max*fBinningResolutionHiGain);
-  Int_t sample_iter = 0;
-
-  while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain )
-    {
-      if (t_iter > fBinningResolutionHalfHiGain-1)
-        {
-          t_iter -= fBinningResolutionHiGain;
-          max_p--; 
-          sample_iter--;
-        }
-      if (t_iter < -fBinningResolutionHalfHiGain)
-        {
-          t_iter += fBinningResolutionHiGain;
-          max_p++; 
-          sample_iter++;
-        }
-    }
-
-  sum = 0.;
-  time_sum = 0.;
-  if (max_p < 0)
-    max_p = 0;
-  if (max_p+fWindowSizeHiGain > range)
-    max_p = range-fWindowSizeHiGain;
-  //
-  // Slide with a window of size fWindowSizeHiGain over the sample 
-  // and multiply the entries with the corresponding weights
-  //
-  for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
-  {
-    const Int_t   idx = fArrBinningResHalfHiGain[sample] + t_iter;
-    const Int_t   ids = max_p + sample;
-    const Float_t pex = fHiGainSignal[ids];
-    sum              += fAmpWeightsHiGain [idx]*pex; 
-    time_sum         += fTimeWeightsHiGain[idx]*pex;
-  }
-
-  if (sum == 0)
-  {
-      time = -1;
-      dsum = 0;
-      dtime = 0;
-      return;
-  }
-
-  // here, the first high-gain slice is already included in the fTimeShiftHiGain
-//  time = fTimeShiftHiGain + max_p - Float_t(t_iter)/fBinningResolutionHiGain;
-  time = max_p + fTimeShiftHiGain + (Float_t)fHiGainFirst /* this shifts the time to the start of the rising edge */
-      - ((Float_t)t_iter)/fBinningResolutionHiGain;
-
-  const Float_t timefineadjust = time_sum/sum;
-  
-  if (TMath::Abs(timefineadjust) < 4./fBinningResolutionHiGain)
-    time -= timefineadjust;
-
+    FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range,
+                             sum, dsum, time, dtime,
+                             sat, 0);
 }
 
@@ -521 +432 @@
                                                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
 {
-
-  const Int_t range = fLoGainLast - fLoGainFirst + 1;
-
-  const Byte_t *end = ptr + range;
-  Byte_t *p     = ptr;
-  //
-  // Prepare the low-gain pedestal
-  //
-  const Float_t pedes  = ped.GetPedestal();
-  const Float_t ABoffs = ped.GetPedestalABoffset();
-        
-  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
-
-  //
-  // Check for saturation in all other slices
-  //
-  Float_t *sample = fLoGainSignal.GetArray();
-  Int_t    ids    = fLoGainFirst;
-  while (p<end)
-    {
-      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
-
-      if (*p++ >= fSaturationLimit)
-        sat++;
-    }
-
-  //
-  // Slide with a window of size fWindowSizeLoGain over the sample 
-  // and multiply the entries with the corresponding weights
-  //
-  if (IsNoiseCalculation())
-    {
-      if (fRandomIter == fBinningResolutionLoGain)
-        fRandomIter = 0;
-      for (Int_t ids=0; ids < fWindowSizeLoGain; ids++)
-        {
-          const Int_t   idx = fArrBinningResLoGain[ids] + fRandomIter;
-          sum              += fAmpWeightsLoGain [idx]*fLoGainSignal[ids]; 
-        }
-      return;
-    }
-
-  Float_t time_sum  = 0.;
-  Float_t fmax      = -FLT_MAX;
-  Float_t ftime_max = 0.;
-  Int_t   max_p     = -1;
-  
-  //
-  // Calculate the sum of the first fWindowSize slices
-  //
-  for (Int_t i=0;i<range-fWindowSizeLoGain+1;i++)
-    {
-      sum      = 0.;
-      time_sum = 0.;
-      
-      //
-      // Slide with a window of size fWindowSizeLoGain over the sample 
-      // and multiply the entries with the corresponding weights
-      //
-      for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
-      {
-        const Int_t   idx = fArrBinningResHalfLoGain[sample];
-        const Float_t pex = fLoGainSignal[sample+i];
-        sum              += fAmpWeightsLoGain [idx]*pex; 
-        time_sum         += fTimeWeightsLoGain[idx]*pex;
-      }
-
-      if (sum>fmax)
-      {
-        fmax      = sum;
-        ftime_max = time_sum;
-        max_p     = i;
-      }
-    } /*   for (Int_t i=0;i<range-fWindowSizeLoGain+1;i++) */
-
-  if (max_p<0)
-  {
-      sum = 0;
-      time = -1;
-      dsum = -1;
-      dtime = -1;
-      return;
-  }
-
-  if (fmax==0)
-  {
-      sum = 0;
-      time = -1;
-      dsum = 0;
-      dtime = 0;
-      return;
-  }
-
-  ftime_max        /= fmax;
-  Int_t t_iter      = Int_t(ftime_max*fBinningResolutionLoGain);
-  Int_t sample_iter = 0;
-
-  while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain )
-    {
-      if (t_iter > fBinningResolutionHalfLoGain-1)
-        {
-          t_iter -= fBinningResolutionLoGain;
-          max_p--; 
-          sample_iter--;
-        }
-      if (t_iter < -fBinningResolutionHalfLoGain)
-        {
-          t_iter += fBinningResolutionLoGain;
-          max_p++; 
-          sample_iter++;
-        }
-    }
-  
-  sum = 0.;
-  time_sum = 0.;
-
-  //
-  // Slide with a window of size fWindowSizeLoGain over the sample 
-  // and multiply the entries with the corresponding weights
-  //
-  for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
-  {
-    const Int_t   idx = fArrBinningResHalfLoGain[sample] + t_iter;
-    const Int_t   ids = max_p + sample;
-    const Float_t pex = ids < 0 ? 0. : ( ids >= range ? 0. : fLoGainSignal[ids]);
-    sum              += fAmpWeightsLoGain [idx]*pex; 
-    time_sum         += fTimeWeightsLoGain[idx]*pex;
-  }
-
-  if (sum == 0)
-  {
-      time = -1;
-      dsum = 0;
-      dtime = 0;
-      return;
-  }
-
-  time = max_p + fTimeShiftLoGain + (Float_t)fLoGainFirst /* this shifts the time to the start of the rising edge */
-      - ((Float_t)t_iter)/fBinningResolutionLoGain;
-
-  const Float_t timefineadjust = time_sum/sum;
-  
-  if (TMath::Abs(timefineadjust) < 4./fBinningResolutionLoGain)
-    time -= timefineadjust;
-
-}
-
+    const Int_t range = fLoGainLast - fLoGainFirst + 1;
+
+    const Byte_t *end = ptr + range;
+    Byte_t *p     = ptr;
+
+    //
+    // Prepare the low-gain pedestal
+    //
+    const Float_t pedes  = ped.GetPedestal();
+    const Float_t ABoffs = ped.GetPedestalABoffset();
+
+    const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
+
+    //
+    // Check for saturation in all other slices
+    //
+    Float_t *sample = fLoGainSignal.GetArray();
+    Int_t    ids    = fLoGainFirst;
+
+    while (p<end)
+    {
+        *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
+
+        if (*p++ >= fSaturationLimit)
+            sat++;
+    }
+
+    FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range,
+                             sum, dsum, time, dtime,
+                             sat, 0);
+
+}
+*/
 // --------------------------------------------------------------------------
 //
@@ -681 +478 @@
 {
 
-  Byte_t hw = fWindowSizeHiGain;
-  Byte_t lw = fWindowSizeLoGain;
   Bool_t rc = kFALSE;
 
@@ -690 +485 @@
       rc = kTRUE;
   }
+   /*
+  Byte_t hw = fWindowSizeHiGain;
+  Byte_t lw = fWindowSizeLoGain;
 
   if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
@@ -704 +502 @@
   if (rc)
     SetWindowSize(hw, lw);
-  
+
   Bool_t rc2 = kFALSE;
   Int_t brh = fBinningResolutionHiGain;
@@ -725 +523 @@
       rc = kTRUE;
     }
-  
+   */
   if (IsEnvDefined(env, prefix, "WeightsFile", print))
     {
@@ -823 +621 @@
             }
 
-            if (2!=sscanf(str.Data(), "# High Gain Weights:%2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain))
+            if (2!=sscanf(str.Data(), "# High Gain Weights: %2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain))
             {
                 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
@@ -833 +631 @@
             fAmpWeightsHiGain .Set(len);
             fTimeWeightsHiGain.Set(len);
+            fPulseHiGain.Set(len);
             hi = kTRUE;
             continue;
@@ -845 +644 @@
             }
 
-            if (2!=sscanf(str.Data(),"# Low Gain Weights:%2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain))
+            if (2!=sscanf(str.Data(),"# Low Gain Weights: %2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain))
             {
                 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
@@ -855 +654 @@
             fAmpWeightsLoGain .Set(len);
             fTimeWeightsLoGain.Set(len);
+            fPulseLoGain.Set(len);
             lo = kTRUE;
             continue;
@@ -867 +667 @@
             continue;
 
-        if (2!=sscanf(str.Data(), "%f %f",
+        if (3!=sscanf(str.Data(), "%f %f %f",
                       lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],
-                      lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]))
+                      lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt],
+                      lo ? &fPulseLoGain[cnt] : &fPulseHiGain[cnt]))
         {
             *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
@@ -903 +704 @@
     *fLog << "with a resolution of " << fBinningResolutionLoGain << endl;
 
-    CalcBinningResArrays();
+    //CalcBinningResArrays();
 
     switch (fWindowSizeHiGain)
@@ -1058 +859 @@
     return ReadWeightsFile(name, path);
 }
-
+/*
 //----------------------------------------------------------------------------
 //
@@ -1125 +926 @@
   Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
 
-  for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++)
+  for (Int_t i = -fBinningResolutionHiGain/2+1; i<=fBinningResolutionHiGain/2; i++)
     {
   
@@ -1178 +979 @@
       for (Int_t count=0; count < fWindowSizeHiGain; count++)
         {
-          const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1;
+          const Int_t idx = i+fBinningResolutionHiGain/2+fBinningResolutionHiGain*count-1;
           fAmpWeightsHiGain [idx] = w_amp [0][count];
           fTimeWeightsHiGain[idx] = w_time[0][count];
@@ -1246 +1047 @@
       // Loop over relative time in one BinningResolution interval
       //
-      Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain;
-      
-      for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++)
+      Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionLoGain/2;
+      
+      for (Int_t i = -fBinningResolutionLoGain/2+1; i<=fBinningResolutionLoGain/2; i++)
         {
           
@@ -1301 +1102 @@
           for (Int_t count=0; count < fWindowSizeLoGain; count++)
             {
-              const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1;
+              const Int_t idx = i+fBinningResolutionLoGain/2+fBinningResolutionLoGain*count-1;
               fAmpWeightsLoGain [idx] = w_amp [0][count];
               fTimeWeightsLoGain[idx] = w_time[0][count];
@@ -1329 +1130 @@
   return kTRUE;
 }
-
+*/
 void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
 {
@@ -1336 +1137 @@
 
     MExtractTimeAndCharge::Print(o);
-    *fLog << " Time Shift  HiGain: " << fTimeShiftHiGain         << "  LoGain: " << fTimeShiftLoGain << endl;
     *fLog << " Window Size HiGain: " << fWindowSizeHiGain        << "  LoGain: " << fWindowSizeLoGain << endl;
     *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << "  LoGain: " << fBinningResolutionHiGain << endl;

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h

@@ -22 +22 @@
 {
 private:
-  static const Byte_t fgHiGainFirst;             //! Default for fHiGainFirst       (now set to: 0)
-  static const Byte_t fgHiGainLast;              //! Default for fHiGainLast        (now set to:14)
-  static const Byte_t fgLoGainFirst;             //! Default for fLoGainFirst       (now set to: 3)
-  static const Byte_t fgLoGainLast;              //! Default for fLoGainLast        (now set to:14)
-  static const Int_t  fgWindowSizeHiGain;        //! Default for fWindowSizeHiGain  (now set to: 6)
-  static const Int_t  fgWindowSizeLoGain;        //! Default for fWindowSizeLoGain  (now set to: 6)
-  static const Int_t  fgBinningResolutionHiGain; //! Default for fBinningResolutionHiGain (now set to: 10)
-  static const Int_t  fgBinningResolutionLoGain; //! Default for fBinningResolutionLoGain (now set to: 10)
-  static const Int_t  fgSignalStartBinHiGain;    //! Default for fSignalStartBinHiGain (now set to: 4)
-  static const Int_t  fgSignalStartBinLoGain;    //! Default for fSignalStartBinLoGain (now set to: 4)
-  static const TString fgNameWeightsFile;        //! "cosmics_weights.dat"
-  static const Float_t fgOffsetLoGain;           //! Default for fOffsetLoGain (now set to 1.7)
-  static const Float_t fgLoGainStartShift;       //! Default for fLoGainStartShift (now set to -1.8)
+    static const Byte_t fgHiGainFirst;             //! Default for fHiGainFirst       (now set to: 0)
+    static const Byte_t fgHiGainLast;              //! Default for fHiGainLast        (now set to:14)
+    static const Byte_t fgLoGainFirst;             //! Default for fLoGainFirst       (now set to: 3)
+    static const Byte_t fgLoGainLast;              //! Default for fLoGainLast        (now set to:14)
+//    static const Int_t  fgWindowSizeHiGain;        //! Default for fWindowSizeHiGain  (now set to: 6)
+//    static const Int_t  fgWindowSizeLoGain;        //! Default for fWindowSizeLoGain  (now set to: 6)
+    static const Int_t  fgBinningResolutionHiGain; //! Default for fBinningResolutionHiGain (now set to: 10)
+    static const Int_t  fgBinningResolutionLoGain; //! Default for fBinningResolutionLoGain (now set to: 10)
+    static const Int_t  fgSignalStartBinHiGain;    //! Default for fSignalStartBinHiGain (now set to: 4)
+    static const Int_t  fgSignalStartBinLoGain;    //! Default for fSignalStartBinLoGain (now set to: 4)
+    static const TString fgNameWeightsFile;        //! "cosmics_weights.dat"
+    static const Float_t fgOffsetLoGain;           //! Default for fOffsetLoGain (now set to 1.7)
+//    static const Float_t fgLoGainStartShift;       //! Default for fLoGainStartShift (now set to -1.8)
 
-  MCalibrationPattern  *fCalibPattern;          //! Calibration DM pattern
-  
-  MArrayF fHiGainSignal;                        //! Need fast access to the signals in a float way
-  MArrayF fLoGainSignal;                        //! Store them in separate arrays
+    MCalibrationPattern  *fCalibPattern;          //! Calibration DM pattern
 
-  MArrayI fArrBinningResHiGain;                 //! helping arrays to hold binningres * i 
-  MArrayI fArrBinningResLoGain;                 //! helping arrays to hold binningres * i 
-  MArrayI fArrBinningResHalfHiGain;             //! helping arrays to hold binningres * i + binningreshalf
-  MArrayI fArrBinningResHalfLoGain;             //! helping arrays to hold binningres * i + binningreshalf
-  
-  Float_t fTimeShiftHiGain;                     //  Time shift from when on to apply the filter
-  Float_t fTimeShiftLoGain;                     //  Time shift from when on to apply the filter
-  
-  Int_t   fSignalStartBinHiGain;                //! Start bin from when on to apply weights
-  Int_t   fSignalStartBinLoGain;                //! Start bin from when on to apply weights
+//    MArrayF fHiGainSignal;                        //! Need fast access to the signals in a float way
+//    MArrayF fLoGainSignal;                        //! Store them in separate arrays
 
-  Int_t   fBinningResolutionHiGain;             //  Number of weights per bin High-Gain
-  Int_t   fBinningResolutionHalfHiGain;         //  Half Number of weights per bin High-Gain
-  Int_t   fBinningResolutionLoGain;             //  Number of weights per bin Low-Gain
-  Int_t   fBinningResolutionHalfLoGain;         //  Half Number of weights per bin Low-Gain
-  
-  MArrayF fAmpWeightsHiGain;                    //! Amplitude weights High-Gain (from weights file)
-  MArrayF fTimeWeightsHiGain;                   //! Time weights High-Gain (from weights file)
-  MArrayF fAmpWeightsLoGain;                    //! Amplitude weights Low-Gain (from weights file)
-  MArrayF fTimeWeightsLoGain;                   //! Time weights Low-Gain (from weights file)
+//    Int_t   fSignalStartBinHiGain;                //! Start bin from when on to apply weights
+//    Int_t   fSignalStartBinLoGain;                //! Start bin from when on to apply weights
 
-  TString fNameWeightsFile;                     // Name of the weights file
-  Bool_t  fAutomaticWeights;                    // Flag whether weight should be determined automatically
-  TString fNameWeightsFileSet;                  //! Flag if weights have alreayd been set
-  Int_t   fRandomIter;                          //! Counter used to randomize weights for noise calculation
+    Int_t   fBinningResolutionHiGain;             //  Number of weights per bin High-Gain
+    Int_t   fBinningResolutionLoGain;             //  Number of weights per bin Low-Gain
 
-  // MExtractTimeAndChargeDigitalFilter
-  void    CalcBinningResArrays();
-  Int_t   GetAutomaticWeights();
-  Bool_t  GetWeights();
-  Int_t   ReadWeightsFile(TString filename, TString path="");
-  TString CompileWeightFileName(TString path, const TString &name) const;
+    MArrayF fAmpWeightsHiGain;                    //! Amplitude weights High-Gain (from weights file)
+    MArrayF fTimeWeightsHiGain;                   //! Time weights High-Gain (from weights file)
+    MArrayF fAmpWeightsLoGain;                    //! Amplitude weights Low-Gain (from weights file)
+    MArrayF fTimeWeightsLoGain;                   //! Time weights Low-Gain (from weights file)
 
-  // MExtractTimeAndCharge
-  Bool_t InitArrays();
+    MArrayF fPulseHiGain;                    //! 
+    MArrayF fPulseLoGain;                    //! 
 
-  // MTask
-  Int_t PreProcess(MParList *pList);
-  Int_t Process();
+    TString fNameWeightsFile;                     // Name of the weights file
+    Bool_t  fAutomaticWeights;                    // Flag whether weight should be determined automatically
+    TString fNameWeightsFileSet;                  //! Flag if weights have alreayd been set
+    Int_t   fRandomIter;                          //! Counter used to randomize weights for noise calculation
+
+    // MExtractTimeAndChargeDigitalFilter
+    void    CalcBinningResArrays();
+    Int_t   GetAutomaticWeights();
+    Bool_t  GetWeights();
+    Int_t   ReadWeightsFile(TString filename, TString path="");
+    TString CompileWeightFileName(TString path, const TString &name) const;
+
+
+    // MExtractTimeAndCharge
+    Bool_t InitArrays();
+
+    // MTask
+    Int_t PreProcess(MParList *pList);
+    Int_t Process();
 
 protected:
-  // MParContainer
-  Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+    // MParContainer
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 
 public:
-  MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL);
-  ~MExtractTimeAndChargeDigitalFilter() { }
+    MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL);
+    ~MExtractTimeAndChargeDigitalFilter() { }
+/*
+    Bool_t WriteWeightsFile(TString filename,
+                            TH1F *shapehi, TH2F *autocorrhi,
+                            TH1F *shapelo=NULL, TH2F *autocorrlo=NULL );
 
-  Bool_t WriteWeightsFile(TString filename,
-                          TH1F *shapehi, TH2F *autocorrhi,
-                          TH1F *shapelo=NULL, TH2F *autocorrlo=NULL );
+  */
+    void SetNameWeightsFile(TString s="")
+    {
+        s.ReplaceAll("\015", ""); // This is a fix for TEnv files edited with windows editors
+        fNameWeightsFile = s;
+        fNameWeightsFileSet="";
+    }
 
+    void EnableAutomaticWeights(Bool_t b=kTRUE) { fAutomaticWeights = b; }
 
-  void SetNameWeightsFile(TString s="") { s.ReplaceAll("\015", ""); // This is a fix for TEnv files edited with windows editors
-      fNameWeightsFile = s; fNameWeightsFileSet=""; }
-  void EnableAutomaticWeights(Bool_t b=kTRUE) { fAutomaticWeights = b; }
+    void SetBinningResolution(const Int_t rh=fgBinningResolutionHiGain, const Int_t rl=fgBinningResolutionLoGain)
+    {
+        fBinningResolutionHiGain = rh & ~1;
+        fBinningResolutionLoGain = rl & ~1;
+    }
+/*
+    void SetSignalStartBin( const Int_t sh=fgSignalStartBinHiGain, const Int_t sl=fgSignalStartBinLoGain)
+    {
+        fSignalStartBinHiGain = sh;
+        fSignalStartBinLoGain = sl;
+    }
 
-  void SetBinningResolution(const Int_t rh=fgBinningResolutionHiGain, const Int_t rl=fgBinningResolutionLoGain)  {
-    fBinningResolutionHiGain     = rh & ~1;
-    fBinningResolutionHalfHiGain = fBinningResolutionHiGain/2;
-    fBinningResolutionLoGain     = rl & ~1;
-    fBinningResolutionHalfLoGain = fBinningResolutionLoGain/2;
-  }
-  
-  void SetSignalStartBin( const Int_t sh=fgSignalStartBinHiGain, const Int_t sl=fgSignalStartBinLoGain) {
-    fSignalStartBinHiGain = sh;
-    fSignalStartBinLoGain = sl; 
-  }
+ */
+    void SetWindowSize( Int_t windowh, Int_t windowl);
+    const char* GetNameWeightsFile() const  { return fNameWeightsFile.Data(); }
 
-  void SetWindowSize( Int_t windowh=fgWindowSizeHiGain, Int_t windowl=fgWindowSizeLoGain);
+    void Print(Option_t *o="") const; //*MENU*
+/*
+    void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum,
+                                 Float_t &time, Float_t &dtime,
+                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
+    void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum,  Float_t &dsum,
+                                 Float_t &time, Float_t &dtime,
+                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
+ */
+    void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum,
+                                  Float_t &time, Float_t &dtime,
+                                  Byte_t sat, Int_t maxpos);
 
-  const char* GetNameWeightsFile() const  { return fNameWeightsFile.Data(); }
+    void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum,  Float_t &dsum,
+                                  Float_t &time, Float_t &dtime,
+                                  Byte_t sat, Int_t maxpos);
 
-  void Print(Option_t *o="") const; //*MENU*
-
-  void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum,
-                               Float_t &time, Float_t &dtime,
-                               Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
-  void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum,  Float_t &dsum,
-                               Float_t &time, Float_t &dtime,
-                               Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
-  
-  ClassDef(MExtractTimeAndChargeDigitalFilter, 2)   // Hendrik's digital filter
+    ClassDef(MExtractTimeAndChargeDigitalFilter, 2)   // Hendrik's digital filter
 };
 

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

@@ -1 +1 @@
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeSpline.cc,v 1.61 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
 !
 ! *
@@ -15 +17 @@
 ! *
 !
-!   Author(s): Markus Gaug       09/2004 <mailto:markus@ifae.es> 
+!   Author(s): Thomas Bretz <mailto:tbretz@astro.uni-wuerzbrug.de>
+!   Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
 !
-!   Copyright: MAGIC Software Development, 2002-2004
+!   Copyright: MAGIC Software Development, 2002-2006
 !
 !
@@ -142 +145 @@
 #include "MExtractTimeAndChargeSpline.h"
 
+#include "MExtralgoSpline.h"
+
 #include "MPedestalPix.h"
 
@@ -156 +161 @@
 const Byte_t  MExtractTimeAndChargeSpline::fgLoGainLast       = 14;
 const Float_t MExtractTimeAndChargeSpline::fgResolution       = 0.05;
-const Float_t MExtractTimeAndChargeSpline::fgRiseTimeHiGain   = 0.5;
-const Float_t MExtractTimeAndChargeSpline::fgFallTimeHiGain   = 0.5;
+const Float_t MExtractTimeAndChargeSpline::fgRiseTimeHiGain   = 0.7;
+const Float_t MExtractTimeAndChargeSpline::fgFallTimeHiGain   = 1.0;
 const Float_t MExtractTimeAndChargeSpline::fgLoGainStretch    = 1.5;
 const Float_t MExtractTimeAndChargeSpline::fgOffsetLoGain     = 1.3;
-const Float_t MExtractTimeAndChargeSpline::fgLoGainStartShift = -2.4;
+//const Float_t MExtractTimeAndChargeSpline::fgLoGainStartShift = -2.4;
 
 // --------------------------------------------------------------------------
@@ -177 +182 @@
 //
 MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title) 
-    : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), 
-      fRandomIter(0), fExtractionType(kIntegral)
+    : fRandomIter(0), fExtractionType(kIntegral)
 {
 
@@ -187 +191 @@
   SetLoGainStretch();
   SetOffsetLoGain(fgOffsetLoGain);
-  SetLoGainStartShift(fgLoGainStartShift);
+//  SetLoGainStartShift(fgLoGainStartShift);
 
   SetRiseTimeHiGain();
@@ -310 +314 @@
 {
 
-  Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
-
-  fHiGainSignal     .Set(range);
+  Int_t range = fHiGainLast - fHiGainFirst + 1;// + fHiLoLast;
+
+//  fHiGainSignal     .Set(range);
   fHiGainFirstDeriv .Set(range);
   fHiGainSecondDeriv.Set(range);
@@ -318 +322 @@
   range = fLoGainLast - fLoGainFirst + 1;
 
-  fLoGainSignal     .Set(range);
+//  fLoGainSignal     .Set(range);
   fLoGainFirstDeriv .Set(range);
   fLoGainSecondDeriv.Set(range);
 
-  fHiGainSignal     .Reset();
+//  fHiGainSignal     .Reset();
   fHiGainFirstDeriv .Reset();
   fHiGainSecondDeriv.Reset();
 
-  fLoGainSignal     .Reset();
+//  fLoGainSignal     .Reset();
   fLoGainFirstDeriv .Reset();
   fLoGainSecondDeriv.Reset();
@@ -361 +365 @@
 }
 
+void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain2(const Float_t *ptr, Int_t num,
+                                                           Float_t &sum, Float_t &dsum,
+                                                           Float_t &time, Float_t &dtime,
+                                                           Byte_t sat, Int_t maxpos)
+{
+    //    const Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
+
+    // Do some handling if maxpos is last slice!
+
+    MExtralgoSpline s(ptr, num, fHiGainFirstDeriv.GetArray(), fHiGainSecondDeriv.GetArray());
+
+    s.SetRiseFallTime(fRiseTimeHiGain, fFallTimeHiGain);
+    s.SetResolution(fResolution);
+
+    if (IsNoiseCalculation())
+    {
+        if (fRandomIter == int(1./fResolution))
+            fRandomIter = 0;
+
+        sum = s.ExtractNoise(fRandomIter);
+        fRandomIter++;
+        return;
+    }
+
+    s.Extract(sat, maxpos);
+    s.GetTime(time, dtime);
+    s.GetSignal(sum, dsum);
+}
+
+void MExtractTimeAndChargeSpline::FindTimeAndChargeLoGain2(const Float_t *ptr, Int_t num,
+                                                           Float_t &sum,  Float_t &dsum,
+                                                           Float_t &time, Float_t &dtime,
+                                                           Byte_t sat, Int_t maxpos)
+{
+    MExtralgoSpline s(ptr, num, fLoGainFirstDeriv.GetArray(), fLoGainSecondDeriv.GetArray());
+
+    s.SetRiseFallTime(fRiseTimeLoGain, fFallTimeLoGain);
+    s.SetResolution(fResolution);
+
+    if (IsNoiseCalculation())
+    {
+        if (fRandomIter == int(1./fResolution))
+            fRandomIter = 0;
+
+        sum = s.ExtractNoise(fRandomIter);
+        return;
+    }
+
+    s.Extract(sat, maxpos);
+    s.GetTime(time, dtime);
+    s.GetSignal(sum, dsum);
+}
+
 // --------------------------------------------------------------------------
 //
-// Calculates the arrival time and charge for each pixel 
-//
-void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum, 
+// Calculates the arrival time and charge for each pixel
+//
+/*
+void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
                                                           Float_t &time, Float_t &dtime, 
                                                           Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
 {
-  Int_t range = fHiGainLast - fHiGainFirst + 1;
-  const Byte_t *end = first + range;
-  Byte_t       *p  = first;
-
-  sat = 0;
-  dsum = 0; // In all cases the extracted signal is valid
-
-  const Float_t pedes  = ped.GetPedestal();
-  const Float_t ABoffs = ped.GetPedestalABoffset();
-
-  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
-
-  fAbMax         = 0.;
-  fAbMaxPos      = 0.;
-  fHalfMax       = 0.;
-  fMaxBinContent = 0;
-  Int_t  maxpos  = 0;
-
-  //
-  // Check for saturation in all other slices
-  //
-  Int_t ids = fHiGainFirst;
-  Float_t *sample = fHiGainSignal.GetArray();
-  while (p<end)
-    {
-
-      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
-
-      if (*p > fMaxBinContent)
+    Int_t range = fHiGainLast - fHiGainFirst + 1;
+    const Byte_t *end = first + range;
+    Byte_t       *p  = first;
+
+    const Float_t pedes  = ped.GetPedestal();
+    const Float_t ABoffs = ped.GetPedestalABoffset();
+
+    const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
+
+    Byte_t maxcont = 0;
+    Int_t  maxpos  = 0;
+    sat = 0;
+
+    //
+    // Check for saturation in all other slices
+    //
+    Int_t ids = fHiGainFirst;
+    Float_t *sample = fHiGainSignal.GetArray();
+    while (p<end)
+    {
+        *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
+
+        if (*p > maxcont)
         {
-          maxpos = ids-fHiGainFirst-1;
-          // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
-          if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
-              fMaxBinContent = *p;
+            maxpos = ids-fHiGainFirst-1;
+            // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
+            //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
+            if (maxpos > 1 && maxpos < (range-1*//* - fWindowSizeHiGain + 1*//*))
+                maxcont = *p;
         }
 
-      if (*p++ >= fSaturationLimit)
-        if (!sat)
-          sat = ids-fHiGainFirst;
-      
-    }
-  
-  if (fHiLoLast != 0)
-    {
-
-      end = logain + fHiLoLast;
-
-      while (logain<end)
+        if (*p++ >= fSaturationLimit)
+            if (!sat)
+                sat = ids-fHiGainFirst;
+    }
+
+    if (fHiLoLast != 0)
+    {
+        end = logain + fHiLoLast;
+
+        while (logain<end)
         {
-
-          *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
-
-          if (*logain > fMaxBinContent)
-            { 
+            *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
+
+            if (*logain > maxcont)
+            {
                 maxpos = ids-fHiGainFirst-1;
+
+                // FIXME!!! MUST BE THERE!
                 // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
                 //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
-                //    fMaxBinContent = *logain;
+                //    maxcont = *logain;
             }
-          
-          if (*logain++ >= fSaturationLimit)
-            if (!sat)
-              sat = ids-fHiGainFirst;
-
-          range++;
+
+            if (*logain++ >= fSaturationLimit)
+                if (!sat)
+                    sat = ids-fHiGainFirst;
+
+            range++;
         }
     }
 
-  fAbMax = fHiGainSignal[maxpos];
-
-  fHiGainSecondDeriv[0] = 0.;
-  fHiGainFirstDeriv[0]  = 0.;
-
-  for (Int_t i=1;i<range-1;i++)
-    {
-      const Float_t pp = fHiGainSecondDeriv[i-1] + 4.;
-      fHiGainSecondDeriv[i] = -1.0/pp;
-      fHiGainFirstDeriv [i] = fHiGainSignal[i+1] - 2*fHiGainSignal[i] + fHiGainSignal[i-1];
-      fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
-    }
-
-  fHiGainSecondDeriv[range-1] = 0.;
-
-  for (Int_t k=range-2;k>=0;k--)
-    fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
-  for (Int_t k=range-2;k>=0;k--)
-    fHiGainSecondDeriv[k] /= 6.;
-  
-  if (IsNoiseCalculation())
-    {
-
-      if (fRandomIter == int(1./fResolution))
-        fRandomIter = 0;
-      
-      const Float_t nsx = fRandomIter * fResolution;
-
-      if (fExtractionType == kAmplitude)
-        {
-          const Float_t b = nsx;
-          const Float_t a = 1. - nsx;
-          
-          sum = a*fHiGainSignal[1]
-            + b*fHiGainSignal[2]
-            + (a*a*a-a)*fHiGainSecondDeriv[1] 
-            + (b*b*b-b)*fHiGainSecondDeriv[2];
-        }
-      else
-          sum = CalcIntegralHiGain(2. + nsx, range);
-
-      fRandomIter++;
-      return;
-    }
-
-  //
-  // Allow no saturated slice and
-  // Don't start if the maxpos is too close to the limits.
-  //
-  const Bool_t limlo = maxpos <       TMath::Ceil(fRiseTimeHiGain);
-  const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeHiGain)-1;
-  if (sat || limlo || limup)
-    {
-      dtime = 1.0;
-      if (fExtractionType == kAmplitude)
-        {
-          sum  = fAbMax;
-          time = (Float_t)(fHiGainFirst + maxpos);
-          return;
-        }
-
-      sum  = CalcIntegralHiGain(limlo ? 0 : range, range);
-      time = (Float_t)(fHiGainFirst + maxpos - 1);
-      return;
-    }
-      
-  dtime = fResolution;
-
-  //
-  // Now find the maximum  
-  //
-  Float_t step    = 0.2; // start with step size of 1ns and loop again with the smaller one
-  Float_t lower   = -1. + maxpos;
-  Float_t upper   = (Float_t)maxpos;
-  fAbMaxPos       = upper;
-  Float_t x       = lower;
-  Float_t y       = 0.;
-  Float_t a       = 1.;
-  Float_t b       = 0.;
-  Int_t   klo     = maxpos-1;
-  Int_t   khi     = maxpos;
-
-  //
-  // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2.
-  // If no maximum is found, go to interval maxpos+1.
-  //
-  while ( x < upper - 0.3 )
-    {
-
-      x += step;
-      a -= step;
-      b += step;
-
-      y = a*fHiGainSignal[klo]
-        + b*fHiGainSignal[khi]
-        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
-        + (b*b*b-b)*fHiGainSecondDeriv[khi];
-
-      if (y > fAbMax)
-        {
-          fAbMax    = y;
-          fAbMaxPos = x;
-        }
-
-    }
-
-  //
-  // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2
-  //
-  if (fAbMaxPos > upper-0.1)
-    {
-      upper   = 1. + maxpos;
-      lower   = (Float_t)maxpos;
-      x       = lower;
-      a       = 1.;
-      b       = 0.;
-      khi     = maxpos+1;
-      klo     = maxpos;
-
-      while (x<upper-0.3)
-        {
-
-          x += step;
-          a -= step;
-          b += step;
-          
-          y = a*fHiGainSignal[klo]
-            + b*fHiGainSignal[khi]
-            + (a*a*a-a)*fHiGainSecondDeriv[klo] 
-            + (b*b*b-b)*fHiGainSecondDeriv[khi];
-
-          if (y > fAbMax)
-            {
-              fAbMax    = y;
-              fAbMaxPos = x;
-            }
-        }
-  }
-  //
-  // Now, the time, abmax and khicont and klocont are set correctly within the previous precision.
-  // Try a better precision. 
-  //
-  const Float_t up = fAbMaxPos+step - 3.0*fResolution;
-  const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
-  const Float_t maxpossave = fAbMaxPos;
-  
-  x     = fAbMaxPos;
-  a     = upper - x;
-  b     = x - lower;
- 
-  step  = 2.*fResolution; // step size of 0.1 FADC slices
- 
-  while (x<up)
-    {
-
-      x += step;
-      a -= step;
-      b += step;
-      
-      y = a*fHiGainSignal[klo]
-        + b*fHiGainSignal[khi]
-        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
-        + (b*b*b-b)*fHiGainSecondDeriv[khi];
-      
-      if (y > fAbMax)
-        {
-          fAbMax    = y;
-          fAbMaxPos = x;
-        }
-    }
-
-  //
-  // Second, try from time down to time-0.2 in steps of fResolution.
-  //
-  x     = maxpossave;
-
-  //
-  // Test the possibility that the absolute maximum has not been found between
-  // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos
-  // which requires new setting of klocont and khicont
-  //
-  if (x < lower + fResolution)
-    {
-      klo--;
-      khi--;
-      upper -= 1.;
-      lower -= 1.;
-    }
-
-  a     = upper - x;
-  b     = x - lower;
-  
-  while (x>lo)
-    {
-
-      x -= step;
-      a += step;
-      b -= step;
-      
-      y = a*fHiGainSignal[klo]
-        + b*fHiGainSignal[khi]
-        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
-        + (b*b*b-b)*fHiGainSecondDeriv[khi];
-      
-      if (y > fAbMax)
-        {
-          fAbMax    = y;
-          fAbMaxPos = x;
-        }
-    }
-
-  if (fExtractionType == kAmplitude)
-    {
-      time  = fAbMaxPos + (Int_t)fHiGainFirst;
-      sum   = fAbMax;
-      return;
-    }
-
-  fHalfMax = fAbMax/2.;
-  
-  //
-  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
-  // First, find the right FADC slice:
-  // 
-  klo  = maxpos;
-  while (klo > 0)
-    {
-      if (fHiGainSignal[--klo] < fHalfMax)
-        break;
-    }
-  
-  khi = klo+1;
-  //
-  // Loop from the beginning of the slice upwards to reach the fHalfMax:
-  // With means of bisection:
-  // 
-  x     = (Float_t)klo;
-  a     = 1.;
-  b     = 0.;
-  
-  step = 0.5;
-  Bool_t back = kFALSE;
-  
-  Int_t maxcnt = 20;
-  Int_t cnt    = 0;
-  
-  while (TMath::Abs(y-fHalfMax) > fResolution)
-    {
-      
-      if (back)
-        {
-          x -= step;
-          a += step;
-          b -= step;
-        }
-      else
-        {
-          x += step;
-          a -= step;
-          b += step;
-        }
-      
-      y = a*fHiGainSignal[klo]
-        + b*fHiGainSignal[khi]
-        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
-        + (b*b*b-b)*fHiGainSecondDeriv[khi];
-      
-      back = y > fHalfMax;
-      
-      if (++cnt > maxcnt)
-        break;
-      
-      step /= 2.;
-    }
-  
-  //
-  // Now integrate the whole thing!
-  // 
-  time = (Float_t)fHiGainFirst + x;
-  sum  = CalcIntegralHiGain(fAbMaxPos - fRiseTimeHiGain, range);
-}
-
+    FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range,
+                             sum, dsum, time, dtime,
+                             sat, maxpos);
+}
 
 // --------------------------------------------------------------------------
@@ -728 +502 @@
                                                           Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) 
 {
-  Int_t range = fLoGainLast - fLoGainFirst + 1;
-  const Byte_t *end = first + range;
-  Byte_t       *p  = first;
-
-  const Float_t pedes  = ped.GetPedestal();
-  const Float_t ABoffs = ped.GetPedestalABoffset();
-
-  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
-
-  fAbMax        = 0.;
-  fAbMaxPos     = 0.;
-  Int_t  maxpos = 0;
-  Int_t  max    = -9999;
-
-  dsum = 0; // In all cases the extracted signal is valid
-
-  //
-  // Check for saturation in all other slices
-  //
-  Int_t    ids    = fLoGainFirst;
-  Float_t *sample = fLoGainSignal.GetArray();
-  while (p<end)
-    {
-
-      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
-
-      if (*p > max)
+    Int_t range = fLoGainLast - fLoGainFirst + 1;
+    const Byte_t *end = first + range;
+    Byte_t       *p  = first;
+
+    const Float_t pedes  = ped.GetPedestal();
+    const Float_t ABoffs = ped.GetPedestalABoffset();
+
+    const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
+
+    Int_t  maxpos = 0;
+    Int_t  max    = -9999;
+
+    //
+    // Check for saturation in all other slices
+    //
+    Int_t    ids    = fLoGainFirst;
+    Float_t *sample = fLoGainSignal.GetArray();
+    while (p<end)
+    {
+        *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
+
+        if (*p > max)
         {
-          maxpos = ids-fLoGainFirst-1;
-          max    = *p;
+            maxpos = ids-fLoGainFirst-1;
+            max    = *p;
         }
 
-      if (*p++ >= fSaturationLimit)
-        sat++;
-    }
-  
-  fAbMax = fLoGainSignal[maxpos];
-  
-  fLoGainSecondDeriv[0] = 0.;
-  fLoGainFirstDeriv[0]  = 0.;
-
-  for (Int_t i=1;i<range-1;i++)
-    {
-      const Float_t pp = fLoGainSecondDeriv[i-1] + 4.;
-      fLoGainSecondDeriv[i] = -1.0/pp;
-      fLoGainFirstDeriv [i] = fLoGainSignal[i+1] - 2*fLoGainSignal[i] + fLoGainSignal[i-1];
-      fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
-    }
-
-  fLoGainSecondDeriv[range-1] = 0.;
-
-  for (Int_t k=range-2;k>=0;k--)
-    fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
-  for (Int_t k=range-2;k>=0;k--)
-    fLoGainSecondDeriv[k] /= 6.;
-  
-  if (IsNoiseCalculation())
-    {
-      if (fRandomIter == int(1./fResolution))
-        fRandomIter = 0;
-      
-      const Float_t nsx = fRandomIter * fResolution;
-
-      if (fExtractionType == kAmplitude)
-        {
-          const Float_t b = nsx;
-          const Float_t a = 1. - nsx;
-          
-          sum = a*fLoGainSignal[1]
-            + b*fLoGainSignal[2]
-            + (a*a*a-a)*fLoGainSecondDeriv[1] 
-            + (b*b*b-b)*fLoGainSecondDeriv[2];
-        }
-      else
-          sum = CalcIntegralLoGain(2. + nsx, range);
-
-      fRandomIter++;
-      return;
-    }
-  //
-  // Allow no saturated slice and
-  // Don't start if the maxpos is too close to the limits.
-  //
-  const Bool_t limlo = maxpos <       TMath::Ceil(fRiseTimeLoGain);
-  const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeLoGain)-1;
-  if (sat || limlo || limup)
-    {
-      dtime = 1.0;
-      if (fExtractionType == kAmplitude)
-        {
-          time = (Float_t)(fLoGainFirst + maxpos);
-          sum = fAbMax;
-          return;
-        }
-
-      sum  = CalcIntegralLoGain(limlo ? 0 : range, range);
-      time = (Float_t)(fLoGainFirst + maxpos - 1);
-      return;
-    }
-
-  dtime = fResolution;
-
-  //
-  // Now find the maximum  
-  //
-  Float_t step    = 0.2; // start with step size of 1ns and loop again with the smaller one
-  Float_t lower   = -1. + maxpos;
-  Float_t upper   = (Float_t)maxpos;
-  fAbMaxPos       = upper;
-  Float_t x       = lower;
-  Float_t y       = 0.;
-  Float_t a       = 1.;
-  Float_t b       = 0.;
-  Int_t   klo     = maxpos-1;
-  Int_t   khi     = maxpos;
-
-  //
-  // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2.
-  // If no maximum is found, go to interval maxpos+1.
-  //
-  while ( x < upper - 0.3 )
-    {
-
-      x += step;
-      a -= step;
-      b += step;
-
-      y = a*fLoGainSignal[klo]
-        + b*fLoGainSignal[khi]
-        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
-        + (b*b*b-b)*fLoGainSecondDeriv[khi];
-
-      if (y > fAbMax)
-        {
-          fAbMax    = y;
-          fAbMaxPos = x;
-        }
-
-    }
-
-  //
-  // Test the possibility that the absolute maximum has not been found before the 
-  // maxpos and search from maxpos to maxpos+1 in steps of 0.2
-  //
-  if (fAbMaxPos > upper-0.1)
-    {
-
-      upper   = 1. + maxpos;
-      lower   = (Float_t)maxpos;
-      x       = lower;
-      a       = 1.;
-      b       = 0.;
-      khi     = maxpos+1;
-      klo     = maxpos;
-
-      while (x<upper-0.3)
-        {
-
-          x += step;
-          a -= step;
-          b += step;
-          
-          y = a*fLoGainSignal[klo]
-            + b*fLoGainSignal[khi]
-            + (a*a*a-a)*fLoGainSecondDeriv[klo] 
-            + (b*b*b-b)*fLoGainSecondDeriv[khi];
-
-          if (y > fAbMax)
-            {
-              fAbMax    = y;
-              fAbMaxPos = x;
-            }
-        }
-  }
-
-
-  //
-  // Now, the time, abmax and khicont and klocont are set correctly within the previous precision.
-  // Try a better precision. 
-  //
-  const Float_t up = fAbMaxPos+step - 3.0*fResolution;
-  const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
-  const Float_t maxpossave = fAbMaxPos;
-  
-  x     = fAbMaxPos;
-  a     = upper - x;
-  b     = x - lower;
- 
-  step  = 2.*fResolution; // step size of 0.1 FADC slice
- 
-  while (x<up)
-    {
-
-      x += step;
-      a -= step;
-      b += step;
-      
-      y = a*fLoGainSignal[klo]
-        + b*fLoGainSignal[khi]
-        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
-        + (b*b*b-b)*fLoGainSecondDeriv[khi];
-      
-      if (y > fAbMax)
-        {
-          fAbMax    = y;
-          fAbMaxPos = x;
-        }
-    }
-
-  //
-  // Second, try from time down to time-0.2 in steps of 0.025.
-  //
-  x     = maxpossave;
-
-  //
-  // Test the possibility that the absolute maximum has not been found between
-  // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos
-  // which requires new setting of klocont and khicont
-  //
-  if (x < lower + fResolution)
-    {
-      klo--;
-      khi--;
-      upper -= 1.;
-      lower -= 1.;
-    }
-
-  a     = upper - x;
-  b     = x - lower;
-  
-  while (x>lo)
-    {
-
-      x -= step;
-      a += step;
-      b -= step;
-      
-      y = a*fLoGainSignal[klo]
-        + b*fLoGainSignal[khi]
-        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
-        + (b*b*b-b)*fLoGainSecondDeriv[khi];
-      
-      if (y > fAbMax)
-        {
-          fAbMax    = y;
-          fAbMaxPos = x;
-        }
-    }
-
-  if (fExtractionType == kAmplitude)
-    {
-      time = fAbMaxPos + (Int_t)fLoGainFirst;
-      sum  = fAbMax;
-      return;
-    }
-  
-  fHalfMax = fAbMax/2.;
-  
-  //
-  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
-  // First, find the right FADC slice:
-  // 
-  klo  = maxpos;
-  while (klo > 0)
-    {
-      klo--;
-      if (fLoGainSignal[klo] < fHalfMax)
-        break;
-    }
-  
-  khi = klo+1;
-  //
-  // Loop from the beginning of the slice upwards to reach the fHalfMax:
-  // With means of bisection:
-  // 
-  x     = (Float_t)klo;
-  a     = 1.;
-  b     = 0.;
-  
-  step = 0.5;
-  Bool_t back = kFALSE;
-  
-  Int_t maxcnt = 20;
-  Int_t cnt    = 0;
-  
-  while (TMath::Abs(y-fHalfMax) > fResolution)
-    {
-      
-      if (back)
-        {
-          x -= step;
-          a += step;
-          b -= step;
-        }
-      else
-        {
-          x += step;
-          a -= step;
-          b += step;
-        }
-      
-      y = a*fLoGainSignal[klo]
-        + b*fLoGainSignal[khi]
-        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
-        + (b*b*b-b)*fLoGainSecondDeriv[khi];
-      
-      back = y > fHalfMax;
-      
-      if (++cnt > maxcnt)
-        break;
-      
-      step /= 2.;
-    }
-  
-  //
-  // Now integrate the whole thing!
-  // 
-  time = x + (Int_t)fLoGainFirst;
-  sum  = CalcIntegralLoGain(fAbMaxPos - fRiseTimeLoGain, range);
-}
-
-Float_t MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t start, Float_t range) const
-{
-    // The number of steps is calculated directly from the integration
-    // window. This is the only way to ensure we are not dealing with
-    // numerical rounding uncertanties, because we always get the same
-    // value under the same conditions -- it might still be different on
-    // other machines!
-    const Float_t step  = 0.2;
-    const Float_t width = fRiseTimeHiGain+fFallTimeHiGain;
-    const Float_t max   = range-1 - (width+step);
-    const Int_t   num   = TMath::Nint(width/step);
-
-    // The order is important. In some cases (limlo-/limup-check) it can
-    // happen than max<0. In this case we start at 0
-    if (start > max)
-        start = max;
-    if (start < 0)
-        start = 0;
-
-    start += step/2;
-
-    Double_t sum = 0.;
-    for (Int_t i=0; i<num; i++)
-    {
-        const Float_t x = start+i*step;
-        const Int_t klo = (Int_t)TMath::Floor(x);
-        const Int_t khi = klo + 1;
-        // Note: if x is close to one integer number (= a FADC sample)
-        // we get the same result by using that sample as klo, and the
-        // next one as khi, or using the sample as khi and the previous
-        // one as klo (the spline is of course continuous). So we do not
-        // expect problems from rounding issues in the argument of
-        // Floor() above (we have noticed differences in roundings
-        // depending on the compilation options).
-
-        const Float_t a = khi - x; // Distance from x to next FADC sample
-        const Float_t b = x - klo; // Distance from x to previous FADC sample
-
-        sum += a*fHiGainSignal[klo]
-            +  b*fHiGainSignal[khi]
-            + (a*a*a-a)*fHiGainSecondDeriv[klo]
-            + (b*b*b-b)*fHiGainSecondDeriv[khi];
-
-        // FIXME? Perhaps the integral should be done analitically
-        // between every two FADC slices, instead of numerically
-    }
-
-    sum *= step; // Transform sum in integral
-    return sum;
-}
-
-Float_t MExtractTimeAndChargeSpline::CalcIntegralLoGain(Float_t start, Float_t range) const
-{
-    // The number of steps is calculated directly from the integration
-    // window. This is the only way to ensure we are not dealing with
-    // numerical rounding uncertanties, because we always get the same
-    // value under the same conditions -- it might still be different on
-    // other machines!
-    const Float_t step  = 0.2;
-    const Float_t width = fRiseTimeLoGain+fFallTimeLoGain;
-    const Float_t max   = range-1 - (width+step);
-    const Int_t   num   = TMath::Nint(width/step);
-
-    // The order is important. In some cases (limlo-/limup-check) it can
-    // happen that max<0. In this case we start at 0
-    if (start > max)
-        start = max;
-    if (start < 0)
-        start = 0;
-
-    start += step/2;
-
-    Double_t sum = 0.;
-    for (Int_t i=0; i<num; i++)
-    {
-        const Float_t x = start+i*step;
-        const Int_t klo = (Int_t)TMath::Floor(x);
-        const Int_t khi = klo + 1;
-        // Note: if x is close to one integer number (= a FADC sample)
-        // we get the same result by using that sample as klo, and the
-        // next one as khi, or using the sample as khi and the previous
-        // one as klo (the spline is of course continuous). So we do not
-        // expect problems from rounding issues in the argument of
-        // Floor() above (we have noticed differences in roundings
-        // depending on the compilation options).
-
-        const Float_t a = khi - x; // Distance from x to next FADC sample
-        const Float_t b = x - klo; // Distance from x to previous FADC sample
-
-        sum += a*fLoGainSignal[klo]
-            +  b*fLoGainSignal[khi]
-            + (a*a*a-a)*fLoGainSecondDeriv[klo]
-            + (b*b*b-b)*fLoGainSecondDeriv[khi];
-
-        // FIXME? Perhaps the integral should be done analitically
-        // between every two FADC slices, instead of numerically
-    }
-    sum *= step; // Transform sum in integral
-    return sum;
-}
+        if (*p++ >= fSaturationLimit)
+            sat++;
+    }
+
+    FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range,
+                             sum, dsum, time, dtime,
+                             sat, maxpos);
+}
+*/
 
 // --------------------------------------------------------------------------
@@ -1198 +587 @@
 
   return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
-
-}
-
-
+}

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

@@ -11 +11 @@
 
 class MPedestalPix;
+
 class MExtractTimeAndChargeSpline : public MExtractTimeAndCharge
 {
+private:
+    static const Byte_t  fgHiGainFirst;      //! Default for fHiGainFirst  (now set to: 2)
+    static const Byte_t  fgHiGainLast;       //! Default for fHiGainLast   (now set to: 14)
+    static const Byte_t  fgLoGainFirst;      //! Default for fLoGainFirst  (now set to: 2)
+    static const Byte_t  fgLoGainLast;       //! Default for fLoGainLast   (now set to: 14)
+    static const Float_t fgResolution;       //! Default for fResolution   (now set to: 0.003)
+    static const Float_t fgRiseTimeHiGain;   //! Default for fRiseTime     (now set to: 1.5)
+    static const Float_t fgFallTimeHiGain;   //! Default for fFallTime     (now set to: 4.5)
+    static const Float_t fgLoGainStretch;    //! Default for fLoGainStretch    (now set to: 1.5)
+    static const Float_t fgOffsetLoGain;     //! Default for fOffsetLoGain     (now set to 1.7)
+//    static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to -1.6)
+
+//    MArrayF fHiGainSignal;                   //! Need fast access to the signals in a float way
+//    MArrayF fLoGainSignal;                   //! Store them in separate arrays
+    MArrayF fHiGainFirstDeriv;               //! High-gain discretized first derivatives
+    MArrayF fLoGainFirstDeriv;               //! Low-gain discretized first derivatives
+    MArrayF fHiGainSecondDeriv;              //! High-gain discretized second derivatives
+    MArrayF fLoGainSecondDeriv;              //! Low-gain discretized second derivatives
+
+    Float_t fResolution;                     // The time resolution in FADC units
+
+    Float_t fRiseTimeHiGain;                 // The usual rise time of the pulse in the high-gain
+    Float_t fFallTimeHiGain;                 // The usual fall time of the pulse in the high-gain
+    Float_t fRiseTimeLoGain;                 // The usual rise time of the pulse in the low-gain
+    Float_t fFallTimeLoGain;                 // The usual fall time of the pulse in the low-gain
+
+    Float_t fLoGainStretch;                  // The stretch of the low-gain w.r.t. the high-gain pulse
+
+    Int_t   fRandomIter;                     //! Counter used to randomize weights for noise calculation
+
+    Int_t   ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+    Bool_t  InitArrays();
+
+public:
+    enum ExtractionType_t { kAmplitude, kIntegral };    //! Possible time and charge extraction types
 
 private:
-  
-  static const Byte_t  fgHiGainFirst;      //! Default for fHiGainFirst  (now set to: 2)
-  static const Byte_t  fgHiGainLast;       //! Default for fHiGainLast   (now set to: 14)
-  static const Byte_t  fgLoGainFirst;      //! Default for fLoGainFirst  (now set to: 2)
-  static const Byte_t  fgLoGainLast;       //! Default for fLoGainLast   (now set to: 14)
-  static const Float_t fgResolution;       //! Default for fResolution   (now set to: 0.003)
-  static const Float_t fgRiseTimeHiGain;   //! Default for fRiseTime     (now set to: 1.5)
-  static const Float_t fgFallTimeHiGain;   //! Default for fFallTime     (now set to: 4.5)
-  static const Float_t fgLoGainStretch;    //! Default for fLoGainStretch    (now set to: 1.5)
-  static const Float_t fgOffsetLoGain;     //! Default for fOffsetLoGain     (now set to 1.7)
-  static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to -1.6)
-  
-  MArrayF fHiGainSignal;                //! Need fast access to the signals in a float way
-  MArrayF fLoGainSignal;                //! Store them in separate arrays
-  MArrayF fHiGainFirstDeriv;            //! High-gain discretized first derivatives
-  MArrayF fLoGainFirstDeriv;            //! Low-gain discretized first derivatives
-  MArrayF fHiGainSecondDeriv;           //! High-gain discretized second derivatives
-  MArrayF fLoGainSecondDeriv;           //! Low-gain discretized second derivatives
-
-  Float_t fAbMax;                       //! Current maximum of the spline
-  Float_t fAbMaxPos;                    //! Current position of the maximum of the spline
-  Float_t fHalfMax;                     //! Current half maximum of the spline
-
-  Float_t fResolution;                  // The time resolution in FADC units
-
-  Float_t fRiseTimeHiGain;              // The usual rise time of the pulse in the high-gain
-  Float_t fFallTimeHiGain;              // The usual fall time of the pulse in the high-gain
-  Float_t fRiseTimeLoGain;              // The usual rise time of the pulse in the low-gain
-  Float_t fFallTimeLoGain;              // The usual fall time of the pulse in the low-gain
-
-  Float_t fLoGainStretch;               // The stretch of the low-gain w.r.t. the high-gain pulse
-
-  Int_t   fRandomIter;                  //! Counter used to randomize weights for noise calculation
-  
-  Int_t   ReadEnv(const TEnv &env, TString prefix, Bool_t print);
-
-  Bool_t  InitArrays();
-  
-  Float_t CalcIntegralHiGain(Float_t start, Float_t range) const;
-  Float_t CalcIntegralLoGain(Float_t start, Float_t range) const;
-
-public:  
-  enum ExtractionType_t { kAmplitude, kIntegral };    //! Possible time and charge extraction types
-
-private:
-  
-  ExtractionType_t fExtractionType;
+    ExtractionType_t fExtractionType;
 
 public:
+    MExtractTimeAndChargeSpline(const char *name=NULL, const char *title=NULL);
+    ~MExtractTimeAndChargeSpline() {}
 
-  MExtractTimeAndChargeSpline(const char *name=NULL, const char *title=NULL);
-  ~MExtractTimeAndChargeSpline() {}
+    Float_t GetRiseTimeHiGain() const { return fRiseTimeHiGain; }
+    Float_t GetFallTimeHiGain() const { return fFallTimeHiGain; }
 
-  Float_t GetRiseTimeHiGain() const { return fRiseTimeHiGain; }
-  Float_t GetFallTimeHiGain() const { return fFallTimeHiGain; }
+    void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0 );
 
-  void SetRange      ( Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0 );  
-  void SetResolution ( const Float_t f=fgResolution  )  { fResolution  = f;  }
-  void SetRiseTimeHiGain( const Float_t f=fgRiseTimeHiGain    )
+    void SetResolution(const Float_t f=fgResolution)  { fResolution  = f;  }
+
+    void SetRiseTimeHiGain(const Float_t f=fgRiseTimeHiGain)
     {
-      fRiseTimeHiGain    = f;
-      fRiseTimeLoGain    = f*fLoGainStretch;
-      fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
-      fWindowSizeHiGain  = TMath::Nint(fRiseTimeHiGain + fFallTimeHiGain);
+        fRiseTimeHiGain    = f;
+        fRiseTimeLoGain    = f*fLoGainStretch;
+        fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
+        fWindowSizeHiGain  = TMath::Nint(fRiseTimeHiGain + fFallTimeHiGain);
     }
-  void SetFallTimeHiGain( const Float_t f=fgFallTimeHiGain    )
+    void SetFallTimeHiGain(const Float_t f=fgFallTimeHiGain)
     {
-      fFallTimeHiGain    = f;
-      fFallTimeLoGain    = f*fLoGainStretch;
-      fNumHiGainSamples  = fRiseTimeHiGain + fFallTimeHiGain;
-      fNumLoGainSamples  = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.;
-      fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
-      fWindowSizeLoGain  = TMath::Nint(fRiseTimeLoGain + fFallTimeLoGain);
+        fFallTimeHiGain    = f;
+        fFallTimeLoGain    = f*fLoGainStretch;
+        fNumHiGainSamples  = fRiseTimeHiGain + fFallTimeHiGain;
+        fNumLoGainSamples  = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.;
+        fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
+        fWindowSizeLoGain  = TMath::Nint(fRiseTimeLoGain + fFallTimeLoGain);
     }
 
-  void SetLoGainStretch ( const Float_t f=fgLoGainStretch    )  { fLoGainStretch = f;   }
-  
-  void SetChargeType ( const ExtractionType_t typ=kIntegral);
-  
-  void FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
-                               Float_t &time, Float_t &dtime,
-                               Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
-  void FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum,  Float_t &dsum,
-                               Float_t &time, Float_t &dtime,
-                               Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
+    void SetLoGainStretch(const Float_t f=fgLoGainStretch) { fLoGainStretch = f;   }
 
-  ClassDef(MExtractTimeAndChargeSpline, 4)   // Task to Extract Arrival Times and Charges using a Cubic Spline
+    void SetChargeType(const ExtractionType_t typ=kIntegral);
+/*
+    void FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
+                                 Float_t &time, Float_t &dtime,
+                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
+    void FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum,  Float_t &dsum,
+                                 Float_t &time, Float_t &dtime,
+                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
+*/
+    void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum,
+                                  Float_t &time, Float_t &dtime,
+                                  Byte_t sat, Int_t maxpos);
+
+    void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum,  Float_t &dsum,
+                                  Float_t &time, Float_t &dtime,
+                                  Byte_t sat, Int_t maxpos);
+
+    ClassDef(MExtractTimeAndChargeSpline, 4)   // Task to Extract Arrival Times and Charges using a Cubic Spline
 };
 
 #endif
-
-
-

trunk/MagicSoft/Mars/msignal/MExtractedSignalPix.h

@@ -15 +15 @@
   Float_t fExtractedSignalLoGainError; // error of the mean value of the extracted signal
 
-  Byte_t fNumHiGainSaturated;          // Number of first hi-gain slice which has saturated (could be negative if already the first slice saturates)
-  Byte_t fNumLoGainSaturated;          // Number of first lo-gain slices which have saturated
+  Byte_t fNumHiGainSaturated;          // Number of how many hi-gain slices saturated
+  Byte_t fNumLoGainSaturated;          // Number of how many lo-gain slices saturated
+
+  Byte_t GetNumHiGainSaturated() const { return fNumHiGainSaturated; }
+  Byte_t GetNumLoGainSaturated() const { return fNumLoGainSaturated; }
 
 public:
@@ -36 +39 @@
   Float_t GetExtractedSignalLoGainError() const { return fExtractedSignalLoGainError;  }
 
-  Byte_t GetNumHiGainSaturated()          const { return fNumHiGainSaturated;          }
-  Byte_t GetNumLoGainSaturated()          const { return fNumLoGainSaturated;          }
-
   Bool_t IsHiGainSaturated()              const { return fNumHiGainSaturated>0;        }
   Bool_t IsLoGainSaturated()              const { return fNumLoGainSaturated>0;        }

trunk/MagicSoft/Mars/msignal/MExtractor.cc

@@ -93 +93 @@
 #include "MPedestalCam.h"
 #include "MPedestalPix.h"
+//#include "MPedestalSubtractEvt.h"
 
 #include "MExtractedSignalCam.h"
@@ -101 +102 @@
 using namespace std;
 
-const Byte_t  MExtractor::fgSaturationLimit = 250;
+const Byte_t  MExtractor::fgSaturationLimit = 245;
 const TString MExtractor::fgNamePedestalCam = "MPedestalCam";
 const TString MExtractor::fgNameSignalCam   = "MExtractedSignalCam";
@@ -123 +124 @@
 MExtractor::MExtractor(const char *name, const char *title)
     : fResolutionPerPheHiGain(0), fResolutionPerPheLoGain(0),
-      fPedestals(NULL), fSignals(NULL), fRawEvt(NULL), fRunHeader(NULL),
-      fHiLoLast(0), fNumHiGainSamples(0), fNumLoGainSamples(0)
+    fPedestals(NULL), fSignals(NULL), fRawEvt(NULL), fRunHeader(NULL),
+    fSignal(NULL),
+      /*fHiLoLast(0),*/ fNumHiGainSamples(0), fNumLoGainSamples(0)
 {
     fName  = name  ? name  : "MExtractor";
@@ -172 +174 @@
     {
         *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fSignal = (MPedestalSubtractedEvt*)pList->FindObject(AddSerialNumber("MPedestalSubtractedEvt"));
+    if (!fSignal)
+    {
+        *fLog << err << AddSerialNumber("MPedestalSubtractedEvt") << " not found... aborting." << endl;
         return kFALSE;
     }
@@ -231 +240 @@
 Bool_t MExtractor::ReInit(MParList *pList)
 {
-
+    const Int_t numl = fRunHeader->GetNumSamplesLoGain();
+    const Int_t numh = fRunHeader->GetNumSamplesHiGain();
+    const Int_t num  = numh+numl;
+
+    if (fHiGainLast>=num)
+    {
+        *fLog << err << "ERROR - Last hi-gain slice must not exceed " << num-1 << endl;
+        return kFALSE;
+    }
+    if (fLoGainLast>=num)
+    {
+        *fLog << err << "ERROR - Last lo-gain slice must not exceed " << num-1 << endl;
+        return kFALSE;
+    }
+
+/*
   const Int_t logainsamples = fRunHeader->GetNumSamplesLoGain();
 
@@ -289 +313 @@
     }
   
-  return kTRUE;
+  */
+    return kTRUE;
 }
 
@@ -431 +456 @@
     *fLog << " Hi Gain Range:      " << Form("%2d %2d", fHiGainFirst, fHiGainLast) << endl;
     *fLog << " Lo Gain Range:      " << Form("%2d %2d", fLoGainFirst, fLoGainLast) << endl;
-    *fLog << " Gain Overlap to Lo: " << Form("%2d", fHiLoLast)        << endl;
+//    *fLog << " Gain Overlap to Lo: " << Form("%2d", fHiLoLast)        << endl;
     *fLog << " Saturation Lim:     " << Form("%3d", fSaturationLimit) << endl;
     *fLog << " Num Samples Hi/Lo:  " << Form("%2.1f %2.1f", fNumHiGainSamples, fNumLoGainSamples) << endl;

trunk/MagicSoft/Mars/msignal/MExtractor.h

@@ -18 +18 @@
 
 class MPedestalCam;
+class MPedestalSubtractedEvt;
 class MExtractedSignalCam;
 
@@ -42 +43 @@
   MRawEvtData         *fRawEvt;            //! Raw event data (time slices)
   MRawRunHeader       *fRunHeader;         //! RunHeader information
+
+  MPedestalSubtractedEvt *fSignal;         //!
                                        
   Byte_t   fHiGainFirst;                   // First FADC slice nr. to extract the High Gain signal
@@ -48 +51 @@
   Byte_t   fLoGainLast;                    // Last FADC slice nr. to extract the Low Gain signal
                                            
-  Byte_t   fHiLoLast;                      // Number of slices in fLoGainSamples counted for the High-Gain signal
+//  Byte_t   fHiLoLast;                      // Number of slices in fLoGainSamples counted for the High-Gain signal
                                            
   Float_t  fNumHiGainSamples;              // Number High Gain FADC slices used to extract the signal
@@ -91 +94 @@
   Float_t GetResolutionPerPheHiGain() const { return fResolutionPerPheHiGain; }
   Float_t GetResolutionPerPheLoGain() const { return fResolutionPerPheLoGain; }
+  Byte_t  GetSaturationLimit() const { return fSaturationLimit; }
 
   Bool_t  IsNoiseCalculation () const { return fNoiseCalculation; }
@@ -107 +111 @@
 
   void SetPedestals (MPedestalCam *pedcam)   { fPedestals = pedcam; }
+  MPedestalCam *GetPedestals()   { return fPedestals; }
 
   // TObject
   void Clear(Option_t *o="") 
     {
-      fHiGainFirst = fHiGainLast = fLoGainFirst = fLoGainLast = fHiLoLast = 0;
+      fHiGainFirst = fHiGainLast = fLoGainFirst = fLoGainLast = /*fHiLoLast =*/ 0;
     }
 

trunk/MagicSoft/Mars/msignal/Makefile

@@ -21 +21 @@
 INCLUDES =  -I. -I../mbase -I../mgui -I../mraw -I../manalysis \
             -I../mgeom -I../mtools -I../mpedestal -I../mbadpixels \
-            -I../mcalib
+            -I../mcalib -I../mextralgo
 
 # mgui (MCamEvent):         MExtractSignalCam

trunk/MagicSoft/Mars/msignal/calibration_weights_UV46.dat

@@ -1 +1 @@
 # High Gain Weights: 4 10
 # (Amplitude)  (Time)
--0.735581 -0.673266
--0.800386 -0.47606
--0.902801 -0.192719
--0.945945 0.237423
--0.851618 0.770424
--0.598107 1.28444
--0.290663 1.69546
-0.0967702 1.89144
-0.362204 1.98161
-0.610111 2.06506
-0.919902 2.20752
-1.12282 1.94486
-1.46305 1.70309
-1.85728 1.58904
-2.19155 1.21011
-2.37154 0.603168
-2.37187 0.000360131
-2.25639 -0.535921
-2.12495 -0.948388
-1.97328 -1.286
-1.72155 -1.70358
-1.55218 -1.75991
-1.21489 -1.86382
-0.759117 -2.09403
-0.266577 -2.03167
--0.102906 -1.60487
--0.284545 -1.16641
--0.367327 -0.679947
--0.376403 -0.3943
--0.379568 -0.229929
--0.218954 0.0443019
--0.22659 0.163101
--0.195144 0.249649
--0.136182 0.328464
--0.0541723 0.334035
--0.00249719 0.207255
-0.00794832 0.102902
-0.00785676 0.0396167
-0.00459524 0.0111356
-0.00887564 0.0405
+-0.735581 -0.673266      0
+-0.800386 -0.47606       0
+-0.902801 -0.192719      0
+-0.945945 0.237423       0
+-0.851618 0.770424       0
+-0.598107 1.28444        0
+-0.290663 1.69546        0
+0.0967702 1.89144        0
+0.362204 1.98161         0
+0.610111 2.06506         0
+0.919902 2.20752         0
+1.12282 1.94486          0
+1.46305 1.70309          0
+1.85728 1.58904          0
+2.19155 1.21011          0
+2.37154 0.603168         0
+2.37187 0.000360131      0
+2.25639 -0.535921        0
+2.12495 -0.948388        0
+1.97328 -1.286           0
+1.72155 -1.70358         0
+1.55218 -1.75991         0
+1.21489 -1.86382         0
+0.759117 -2.09403        0
+0.266577 -2.03167        0
+-0.102906 -1.60487       0
+-0.284545 -1.16641       0
+-0.367327 -0.679947      0
+-0.376403 -0.3943        0
+-0.379568 -0.229929      0
+-0.218954 0.0443019      0
+-0.22659 0.163101        0
+-0.195144 0.249649       0
+-0.136182 0.328464       0
+-0.0541723 0.334035      0
+-0.00249719 0.207255     0
+0.00794832 0.102902      0
+0.00785676 0.0396167     0
+0.00459524 0.0111356     0
+0.00887564 0.0405        0
 # Low Gain Weights: 6 10
 # (Amplitude)  (Time)
-0.0446612 -0.385273
-0.038191 -0.00418687
-0.0386966 0.0212324
-0.0402881 0.0744799
-0.0415794 0.229615
-0.0598731 0.44332
-0.0758477 0.661518
-0.101509 1.10641
-0.159323 1.64997
-0.497256 2.83685
-0.245087 3.27499
-0.140546 2.46177
-0.58086 2.2849
-0.632721 2.45587
-0.72819 2.52835
-0.889583 2.48099
-0.980812 2.50031
-1.09885 2.55892
-1.21374 2.78769
-1.61928 3.08069
-1.38544 1.95583
-1.31998 1.1792
-1.50633 0.591226
-1.50916 0.0793899
-1.5008 -0.33188
-1.47339 -0.575386
-1.45362 -0.915309
-1.40214 -1.31593
-1.34175 -1.77904
-1.0661 -2.05471
-1.31087 -1.49798
-1.33793 -1.34758
-1.10172 -1.21719
-1.08133 -1.09356
-1.04007 -0.981455
-0.976745 -1.08299
-0.930979 -1.14774
-0.874203 -1.18348
-0.816708 -1.20126
-0.587354 -1.92869
-0.783078 -1.89621
-0.792771 -1.03439
-0.622278 -0.781807
-0.61184 -0.745831
-0.578792 -0.683741
-0.537336 -0.596328
-0.51443 -0.592858
-0.482294 -0.560586
-0.462351 -0.827587
-0.317989 -1.05649
-0.459672 -0.775035
-0.468287 -0.619961
-0.374182 -0.31635
-0.376946 -0.225242
-0.367075 -0.347444
-0.340737 -0.393231
-0.321054 -0.187384
-0.320654 -0.225558
-0.302148 -0.399499
-0.232954 -0.607578
+0.0446612 -0.385273      0
+0.038191 -0.00418687     0
+0.0386966 0.0212324      0
+0.0402881 0.0744799      0
+0.0415794 0.229615       0
+0.0598731 0.44332        0
+0.0758477 0.661518       0
+0.101509 1.10641         0
+0.159323 1.64997         0
+0.497256 2.83685         0
+0.245087 3.27499         0
+0.140546 2.46177         0
+0.58086 2.2849           0
+0.632721 2.45587         0
+0.72819 2.52835          0
+0.889583 2.48099         0
+0.980812 2.50031         0
+1.09885 2.55892          0
+1.21374 2.78769          0
+1.61928 3.08069          0
+1.38544 1.95583          0
+1.31998 1.1792           0
+1.50633 0.591226         0
+1.50916 0.0793899        0
+1.5008 -0.33188          0
+1.47339 -0.575386        0
+1.45362 -0.915309        0
+1.40214 -1.31593         0
+1.34175 -1.77904         0
+1.0661 -2.05471          0
+1.31087 -1.49798         0
+1.33793 -1.34758         0
+1.10172 -1.21719         0
+1.08133 -1.09356         0
+1.04007 -0.981455        0
+0.976745 -1.08299        0
+0.930979 -1.14774        0
+0.874203 -1.18348        0
+0.816708 -1.20126        0
+0.587354 -1.92869        0
+0.783078 -1.89621        0
+0.792771 -1.03439        0
+0.622278 -0.781807       0
+0.61184 -0.745831        0
+0.578792 -0.683741       0
+0.537336 -0.596328       0
+0.51443 -0.592858        0
+0.482294 -0.560586       0
+0.462351 -0.827587       0
+0.317989 -1.05649        0
+0.459672 -0.775035       0
+0.468287 -0.619961       0
+0.374182 -0.31635        0
+0.376946 -0.225242       0
+0.367075 -0.347444       0
+0.340737 -0.393231       0
+0.321054 -0.187384       0
+0.320654 -0.225558       0
+0.302148 -0.399499       0
+0.232954 -0.607578       0

trunk/MagicSoft/Mars/msignal/cosmics_weights46.dat

@@ -1 +1 @@
 # High Gain Weights: 4 10
 # (Amplitude)  (Time)
--0.637443 -0.648763
--0.707555 -0.610204
--0.772465 -0.469694
--0.896081 -0.305143
--1.02783 0.00696361
--1.03684 0.731979
--0.73698 1.49428
--0.23419 1.57509
-0.12767 1.83675
-0.441087 1.70526
-0.639658 1.50546
-0.806192 1.58667
-0.966865 1.41151
-1.24058 1.34506
-1.58293 1.39453
-2.02092 1.46038
-2.31508 0.8539
-2.31705 0.051831
-2.18497 -0.493294
-2.01129 -0.749924
-1.8675 -0.959329
-1.75322 -1.25026
-1.60962 -1.37231
-1.3216 -1.58832
-0.880364 -1.9924
-0.219661 -2.50675
--0.35077 -1.96955
--0.553454 -0.855454
--0.491715 -0.387673
--0.38701 -0.0297507
--0.293549 0.092153
--0.278646 0.194742
--0.263716 0.25195
--0.221238 0.323789
--0.121941 0.381038
-0.00379863 0.393009
-0.0805122 0.201843
-0.0443888 -0.0531952
--0.0199606 -0.175625
--0.0809393 -0.21419
+-0.637443 -0.648763     0
+-0.707555 -0.610204     0
+-0.772465 -0.469694     0
+-0.896081 -0.305143     0
+-1.02783 0.00696361     0
+-1.03684 0.731979       0
+-0.73698 1.49428        0
+-0.23419 1.57509        0
+0.12767 1.83675         0
+0.441087 1.70526        0
+0.639658 1.50546        0
+0.806192 1.58667        0
+0.966865 1.41151        0
+1.24058 1.34506         0
+1.58293 1.39453         0
+2.02092 1.46038         0
+2.31508 0.8539          0
+2.31705 0.051831        0
+2.18497 -0.493294       0
+2.01129 -0.749924       0
+1.8675 -0.959329        0
+1.75322 -1.25026        0
+1.60962 -1.37231        0
+1.3216 -1.58832         0
+0.880364 -1.9924        0
+0.219661 -2.50675       0
+-0.35077 -1.96955       0
+-0.553454 -0.855454     0
+-0.491715 -0.387673     0
+-0.38701 -0.0297507     0
+-0.293549 0.092153      0
+-0.278646 0.194742      0
+-0.263716 0.25195       0
+-0.221238 0.323789      0
+-0.121941 0.381038      0
+0.00379863 0.393009     0
+0.0805122 0.201843      0
+0.0443888 -0.0531952    0
+-0.0199606 -0.175625    0
+-0.0809393 -0.21419     0
 # Low Gain Weights: 6 10
 # (Amplitude)  (Time)
-0.0446612 -0.385273
-0.038191 -0.00418687
-0.0386966 0.0212324
-0.0402881 0.0744799
-0.0415794 0.229615
-0.0598731 0.44332
-0.0758477 0.661518
-0.101509 1.10641
-0.159323 1.64997
-0.497256 2.83685
-0.245087 3.27499
-0.140546 2.46177
-0.58086 2.2849
-0.632721 2.45587
-0.72819 2.52835
-0.889583 2.48099
-0.980812 2.50031
-1.09885 2.55892
-1.21374 2.78769
-1.61928 3.08069
-1.38544 1.95583
-1.31998 1.1792
-1.50633 0.591226
-1.50916 0.0793899
-1.5008 -0.33188
-1.47339 -0.575386
-1.45362 -0.915309
-1.40214 -1.31593
-1.34175 -1.77904
-1.0661 -2.05471
-1.31087 -1.49798
-1.33793 -1.34758
-1.10172 -1.21719
-1.08133 -1.09356
-1.04007 -0.981455
-0.976745 -1.08299
-0.930979 -1.14774
-0.874203 -1.18348
-0.816708 -1.20126
-0.587354 -1.92869
-0.783078 -1.89621
-0.792771 -1.03439
-0.622278 -0.781807
-0.61184 -0.745831
-0.578792 -0.683741
-0.537336 -0.596328
-0.51443 -0.592858
-0.482294 -0.560586
-0.462351 -0.827587
-0.317989 -1.05649
-0.459672 -0.775035
-0.468287 -0.619961
-0.374182 -0.31635
-0.376946 -0.225242
-0.367075 -0.347444
-0.340737 -0.393231
-0.321054 -0.187384
-0.320654 -0.225558
-0.302148 -0.399499
-0.232954 -0.607578
+0.0446612 -0.385273     0
+0.038191 -0.00418687    0
+0.0386966 0.0212324     0
+0.0402881 0.0744799     0
+0.0415794 0.229615      0
+0.0598731 0.44332       0
+0.0758477 0.661518      0
+0.101509 1.10641        0
+0.159323 1.64997        0
+0.497256 2.83685        0
+0.245087 3.27499        0
+0.140546 2.46177        0
+0.58086 2.2849          0
+0.632721 2.45587        0
+0.72819 2.52835         0
+0.889583 2.48099        0
+0.980812 2.50031        0
+1.09885 2.55892         0
+1.21374 2.78769         0
+1.61928 3.08069         0
+1.38544 1.95583         0
+1.31998 1.1792          0
+1.50633 0.591226        0
+1.50916 0.0793899       0
+1.5008 -0.33188         0
+1.47339 -0.575386       0
+1.45362 -0.915309       0
+1.40214 -1.31593        0
+1.34175 -1.77904        0
+1.0661 -2.05471         0
+1.31087 -1.49798        0
+1.33793 -1.34758        0
+1.10172 -1.21719        0
+1.08133 -1.09356        0
+1.04007 -0.981455       0
+0.976745 -1.08299       0
+0.930979 -1.14774       0
+0.874203 -1.18348       0
+0.816708 -1.20126       0
+0.587354 -1.92869       0
+0.783078 -1.89621       0
+0.792771 -1.03439       0
+0.622278 -0.781807      0
+0.61184 -0.745831       0
+0.578792 -0.683741      0
+0.537336 -0.596328      0
+0.51443 -0.592858       0
+0.482294 -0.560586      0
+0.462351 -0.827587      0
+0.317989 -1.05649       0
+0.459672 -0.775035      0
+0.468287 -0.619961      0
+0.374182 -0.31635       0
+0.376946 -0.225242      0
+0.367075 -0.347444      0
+0.340737 -0.393231      0
+0.321054 -0.187384      0
+0.320654 -0.225558      0
+0.302148 -0.399499      0
+0.232954 -0.607578      0