Changeset 8154 for trunk/MagicSoft

trunk/MagicSoft/Mars/Makefile

r8032	r8154
45	45	mhist \
46	46	manalysis \
	47	mextralgo \
47	48	msignal \
48	49	mbadpixels \

trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc

-              r8072
+              r8154
+}
+// -----------------------------------------------------------------------------
+//
+// 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
 …
     // 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++)
 …
+        }
+    }
+    /*
+     // 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
 …
         fTime -= timefineadjust;
+}
 #include <TH1.h>
 …
 #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;
 …
         cout << "        Shape: " << shape.GetNbinsX() << endl;
         cout << "        ACorr: " << autocorr.GetNbinsX() << endl;
         return kFALSE;
+        return -1;
+    }
 …
     // FIXME: DELETE!!!
     TH1 &derivative = *static_cast<TH1*>(shape.Clone());
+    derivative.SetDirectory(0);
     derivative.Reset();
 …
+        {
             cout << "ERROR - Division by zero: denom[0][0]==0 for i=" << i << "." << endl;
             return kFALSE;
+            return -1;
+        }
 …
+    }
     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;
 …
         cout << "        Shape: " << shape.GetNbinsX() << endl;
         cout << "        ACorr: " << autocorr.GetNbinsX() << endl;
         return;
+        return -1;
+    }
 …
     TSpline5 val("Signal", &gr);
+    TH1F derivative(shape);
+    // FIXME: DELETE!!!
+    TH1 &derivative = *static_cast<TH1*>(shape.Clone());
+    derivative.SetDirectory(0);
     derivative.Reset();
 …
         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);
 …
+        }
+    }
+    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

-              r8072
+              r8154
 #define MARS_MExtralgoDigitalFilter
 #ifndef ROOT_TROOT
 #include <TROOT.h>
+#ifndef ROOT_TMatrix
+#include <TMatrix.h>
 #endif
 …
 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
 …
     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
 …
         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++)
+        {
 …
 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; }
 …
     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

-              r7999
+              r8154
+}
+#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)
 …
     fSignalDev = -1;
     fTimeDev   = -1;
+/*
     //
     // Allow no saturated slice and
 …
     //
-    /*
     const Bool_t limlo = maxbin <      TMath::Ceil(fRiseTime);
     const Bool_t limup = maxbin > fNum-TMath::Ceil(fFallTime)-1;
 …
         return;
+    }
-    */
-    fTimeDev = fResolution;
     //
 …
     //
-    /*
     Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
 …
+        }
+    }
+  */
+    fTime      = maxpos;
+    fTimeDev   = fResolution;
+    fSignal    = CalcIntegral(maxpos);
+    fSignalDev = 0;  // means: is valid
+    return;
+*/
     // --- Start NEW ---
 …
     Float_t maxpos, maxval;
+    // FIXME: Check the dfeault if no maximum found!!!
     GetMaxAroundI(maxbin, maxpos, maxval);
 …
+    {
         fTime      = maxpos;
+        fTimeDev   = fResolution;
         fSignal    = maxval;
         fSignalDev = 0;  // means: is valid
 …
+    }
-    //
-    // 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
 …
     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;
+    }
+    */
     //
 …
     //
     //   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

-              r7999
+              r8154
     // Input
     Float_t *fVal;
     Int_t    fNum;
+    Float_t const *fVal;
+    const Int_t    fNum;
     Float_t *fDer1;
 …
         // Take a default in case no maximum is found
+        // FIXME: Check THIS!!!
         xmax=i;
         ymax=fVal[i];
 …
 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

-              r8141
+              r8154
 #include "MFTriggerPattern.h"
 #include "MGeomApply.h"
+#include "MPedestalSubtract.h"
 //#include "MMcPedestalCopy.h"
 #include "MPointingPosCalc.h"
 …
     MContinue conttp(&ftp, "ContTrigPattern");
+    // Create the pedestal subtracted raw-data
+    MPedestalSubtract pedsub;
+    pedsub.SetPedestalCam(&pedcamab);
     // Do signal and pedestal calculation
     MPedCalcFromLoGain     pedlo1("MPedCalcFundamental");
 …
         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());
 …
             pedlo2.SetExtractWindow(f, n);
             pedlo3.SetExtractWindow(f, n);
+            */
+            pedlo2.SetExtractWindow(f, w);
+            pedlo3.SetExtractWindow(f, w);
+        }
+    }
 …
     tlist2.AddToList(&apply);
     tlist2.AddToList(&merge);
+    tlist2.AddToList(&pedsub);
     tlist2.AddToList(&pedlo1);
     tlist2.AddToList(&pedlo2);

trunk/MagicSoft/Mars/mjobs/MJCalibration.cc

-              r8141
+              r8154
 #include "MRawFileRead.h"
 #include "MGeomApply.h"
+#include "MPedestalSubtract.h"
 #include "MTaskEnv.h"
 #include "MBadPixelsMerge.h"
 …
     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

-              r7739
+              r8154
 #include "MRawEvtData.h"
 #include "MGeomApply.h"
+#include "MPedestalSubtract.h"
 #include "MTriggerPatternDecode.h"
 #include "MBadPixelsMerge.h"
 …
     //
     if (!fBadPixelsFile.IsNull())
+      {
+    {
         *fLog << inf << "Excluding: " << fBadPixelsFile << endl;
         ifstream fin(fBadPixelsFile.Data());
         fBadPixels.AsciiRead((istream&)fin);
+      }
+        ifstream fin(fBadPixelsFile);
+        fBadPixels.AsciiRead(fin);
+    }
     MGeomApply geomapl;
 …
     MPedCalcPedRun pedcalc;
     pedcalc.SetPedestalUpdate(kFALSE);
+    //pedcalc.SetPedestalUpdate(kFALSE);
     MPedCalcFromLoGain pedlogain;
 …
         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);
     // ----------------------------------------------------------------------
 …
     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

r8153	r8154
1	1	/* ======================================================================== *\
	2	! $Name: not supported by cvs2svn $:$Id: MPedestalSubtract.cc,v 1.2 2006-10-24 08:18:07 tbretz Exp $
	3	! --------------------------------------------------------------------------
2	4	!
3	5	! *

trunk/MagicSoft/Mars/mpedestal/Makefile

-              r8151
+              r8154
            MHPedestalCor.cc \
            MPedestalSubtract.cc \
+           MPedestalSubtractedEvt.cc \
+           MPedCalcFromData.cc
+           MPedestalSubtractedEvt.cc
 ############################################################

trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h

r8151	r8154
19	19	#pragma link C++ class MPedestalSubtract+;
20	20	#pragma link C++ class MPedestalSubtractedEvt+;
21		~~#pragma link C++ class MPedCalcFromData+;~~
22	21
23	22	#pragma link C++ class MHPedestalCor+;

trunk/MagicSoft/Mars/msignal/MExtractTime.cc

-              r8144
+              r8154
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTime.cc,v 1.22 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
+!
 ! *
 …
+!
 !   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
+!
+!
 …
 // - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
 //                                      fLoGainFirst, fLoGainLast, fNumLoGainSamples);
 //
+/*
 Bool_t MExtractTime::ReInit(MParList *pList)
+{
 …
     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()
+{
 …
       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

-              r7091
+              r8154
   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

-              r7900
+              r8154
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.55 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
+!
 ! *
 …
+!
 !   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
+!
+!
 …
 //   MExtractedSignalCam
 //
+// For weired events check: Run 94127 Event 672, 1028
+//
 //////////////////////////////////////////////////////////////////////////////
 #include "MExtractTimeAndCharge.h"
 …
 #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"
 …
 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;
 // --------------------------------------------------------------------------
 …
 //
 // Sets:
-// - fLoGainFirstSave to 0
 // - fWindowSizeHiGain and fWindowSizeLoGain to 0
 // - fLoGainStartShift to fgLoGainStartShift+fgOffsetLoGain
 …
 //
 MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title)
     : fLoGainFirstSave(0), fWindowSizeHiGain(0), fWindowSizeLoGain(0)
+    : /*fLoGainFirstSave(0),*/ fWindowSizeHiGain(0), fWindowSizeLoGain(0)
+{
     fName  = name  ? name  : "MExtractTimeAndCharge";
 …
   if (fSignals)
     fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
+    fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples,
                                 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
   return kTRUE;
 …
 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;
+}
 …
 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

-              r7091
+              r8154
   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);
 …
   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,
 …
                                        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

-              r8129
+              r8154
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.73 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
+!
 ! *
 …
 #include "MCalibrationPattern.h"
 #include "MExtractedSignalCam.h"
+#include "MExtralgoDigitalFilter.h"
 #include "MPedestalPix.h"
 …
 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;
 …
 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;
 // --------------------------------------------------------------------------
 …
 //
 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";
 …
     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
+}
 …
 //
 // 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;
 …
   // 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);
 …
+}
 // --------------------------------------------------------------------------
 //
 …
 //
 // If new weights are set
-//  fTimeShiftHiGain
-//  fTimeShiftLoGain
 //  fNumHiGainSamples
 //  fNumLoGainSamples
 …
     //
+    // 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);
 …
     // From MExtractTimeAndCharge::ReInit
     if (fSignals)
         fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
+        fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples,
                                     fLoGainFirst, fLoGainLast, fNumLoGainSamples);
     return kTRUE;
 …
         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();
 …
+}
-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);
+}
 …
                                                                  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);
+}
+*/
 // --------------------------------------------------------------------------
 //
 …
+{
-  Byte_t hw = fWindowSizeHiGain;
-  Byte_t lw = fWindowSizeLoGain;
   Bool_t rc = kFALSE;
 …
       rc = kTRUE;
+  }
+   /*
+  Byte_t hw = fWindowSizeHiGain;
+  Byte_t lw = fWindowSizeLoGain;
   if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
 …
   if (rc)
     SetWindowSize(hw, lw);
   Bool_t rc2 = kFALSE;
   Int_t brh = fBinningResolutionHiGain;
 …
       rc = kTRUE;
+    }
+   */
   if (IsEnvDefined(env, prefix, "WeightsFile", print))
+    {
 …
+            }
             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;
 …
             fAmpWeightsHiGain .Set(len);
             fTimeWeightsHiGain.Set(len);
+            fPulseHiGain.Set(len);
             hi = kTRUE;
             continue;
 …
+            }
             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;
 …
             fAmpWeightsLoGain .Set(len);
             fTimeWeightsLoGain.Set(len);
+            fPulseLoGain.Set(len);
             lo = kTRUE;
             continue;
 …
             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;
 …
     *fLog << "with a resolution of " << fBinningResolutionLoGain << endl;
     CalcBinningResArrays();
+    //CalcBinningResArrays();
     switch (fWindowSizeHiGain)
 …
     return ReadWeightsFile(name, path);
+}
+/*
 //----------------------------------------------------------------------------
 //
 …
   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++)
+    {
 …
       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];
 …
       // 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++)
+        {
 …
           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];
 …
   return kTRUE;
+}
+*/
 void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
+{
 …
     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

-              r7194
+              r8154
+{
 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

-              r7930
+              r8154
 /* ======================================================================== *\
+! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeSpline.cc,v 1.61 2006-10-24 08:24:52 tbretz Exp $
+! --------------------------------------------------------------------------
+!
 ! *
 …
 ! *
+!
+!   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
+!
+!
 …
 #include "MExtractTimeAndChargeSpline.h"
+#include "MExtralgoSpline.h"
 #include "MPedestalPix.h"
 …
 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;
 // --------------------------------------------------------------------------
 …
 //
 MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title)
+    : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.),
+      fRandomIter(0), fExtractionType(kIntegral)
+    : fRandomIter(0), fExtractionType(kIntegral)
+{
 …
   SetLoGainStretch();
   SetOffsetLoGain(fgOffsetLoGain);
   SetLoGainStartShift(fgLoGainStartShift);
+//  SetLoGainStartShift(fgLoGainStartShift);
   SetRiseTimeHiGain();
 …
+{
   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);
 …
   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();
 …
+}
+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);
+}
 // --------------------------------------------------------------------------
 …
                                                           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);
+}
+*/
 // --------------------------------------------------------------------------
 …
   return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
+}
+}

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

-              r7491
+              r8154
 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

-              r7876
+              r8154
   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:
 …
   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

-              r7881
+              r8154
 #include "MPedestalCam.h"
 #include "MPedestalPix.h"
+//#include "MPedestalSubtractEvt.h"
 #include "MExtractedSignalCam.h"
 …
 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";
 …
 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";
 …
+    {
         *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;
+    }
 …
 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();
 …
+    }
+  return kTRUE;
+  */
+    return kTRUE;
+}
 …
     *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

-              r7804
+              r8154
 class MPedestalCam;
+class MPedestalSubtractedEvt;
 class MExtractedSignalCam;
 …
   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
 …
   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
 …
   Float_t GetResolutionPerPheHiGain() const { return fResolutionPerPheHiGain; }
   Float_t GetResolutionPerPheLoGain() const { return fResolutionPerPheLoGain; }
+  Byte_t  GetSaturationLimit() const { return fSaturationLimit; }
   Bool_t  IsNoiseCalculation () const { return fNoiseCalculation; }
 …
   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

r7897	r8154
21	21	INCLUDES = -I. -I../mbase -I../mgui -I../mraw -I../manalysis \
22	22	-I../mgeom -I../mtools -I../mpedestal -I../mbadpixels \
23		-I../mcalib
	23	-I../mcalib -I../mextralgo
24	24
25	25	# mgui (MCamEvent): MExtractSignalCam

trunk/MagicSoft/Mars/msignal/calibration_weights_UV46.dat

-              r6870
+              r8154
 # 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
 .0967702 1.89144
 .362204 1.98161
 .610111 2.06506
 .919902 2.20752
 .12282 1.94486
 .46305 1.70309
 .85728 1.58904
 .19155 1.21011
 .37154 0.603168
 .37187 0.000360131
 .25639 -0.535921
 .12495 -0.948388
 .97328 -1.286
 .72155 -1.70358
 .55218 -1.75991
 .21489 -1.86382
 .759117 -2.09403
 .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
 .00794832 0.102902
 .00785676 0.0396167
 .00459524 0.0111356
 .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
+.0967702 1.89144        0
+.362204 1.98161         0
+.610111 2.06506         0
+.919902 2.20752         0
+.12282 1.94486          0
+.46305 1.70309          0
+.85728 1.58904          0
+.19155 1.21011          0
+.37154 0.603168         0
+.37187 0.000360131      0
+.25639 -0.535921        0
+.12495 -0.948388        0
+.97328 -1.286           0
+.72155 -1.70358         0
+.55218 -1.75991         0
+.21489 -1.86382         0
+.759117 -2.09403        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
+.00794832 0.102902      0
+.00785676 0.0396167     0
+.00459524 0.0111356     0
+.00887564 0.0405        0
 # Low Gain Weights: 6 10
 # (Amplitude)  (Time)
 .0446612 -0.385273
 .038191 -0.00418687
 .0386966 0.0212324
 .0402881 0.0744799
 .0415794 0.229615
 .0598731 0.44332
 .0758477 0.661518
 .101509 1.10641
 .159323 1.64997
 .497256 2.83685
 .245087 3.27499
 .140546 2.46177
 .58086 2.2849
 .632721 2.45587
 .72819 2.52835
 .889583 2.48099
 .980812 2.50031
 .09885 2.55892
 .21374 2.78769
 .61928 3.08069
 .38544 1.95583
 .31998 1.1792
 .50633 0.591226
 .50916 0.0793899
 .5008 -0.33188
 .47339 -0.575386
 .45362 -0.915309
 .40214 -1.31593
 .34175 -1.77904
 .0661 -2.05471
 .31087 -1.49798
 .33793 -1.34758
 .10172 -1.21719
 .08133 -1.09356
 .04007 -0.981455
 .976745 -1.08299
 .930979 -1.14774
 .874203 -1.18348
 .816708 -1.20126
 .587354 -1.92869
 .783078 -1.89621
 .792771 -1.03439
 .622278 -0.781807
 .61184 -0.745831
 .578792 -0.683741
 .537336 -0.596328
 .51443 -0.592858
 .482294 -0.560586
 .462351 -0.827587
 .317989 -1.05649
 .459672 -0.775035
 .468287 -0.619961
 .374182 -0.31635
 .376946 -0.225242
 .367075 -0.347444
 .340737 -0.393231
 .321054 -0.187384
 .320654 -0.225558
 .302148 -0.399499
 .232954 -0.607578
+.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         0
+.09885 2.55892          0
+.21374 2.78769          0
+.61928 3.08069          0
+.38544 1.95583          0
+.31998 1.1792           0
+.50633 0.591226         0
+.50916 0.0793899        0
+.5008 -0.33188          0
+.47339 -0.575386        0
+.45362 -0.915309        0
+.40214 -1.31593         0
+.34175 -1.77904         0
+.0661 -2.05471          0
+.31087 -1.49798         0
+.33793 -1.34758         0
+.10172 -1.21719         0
+.08133 -1.09356         0
+.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

trunk/MagicSoft/Mars/msignal/cosmics_weights46.dat

-              r6870
+              r8154
 # 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
 .12767 1.83675
 .441087 1.70526
 .639658 1.50546
 .806192 1.58667
 .966865 1.41151
 .24058 1.34506
 .58293 1.39453
 .02092 1.46038
 .31508 0.8539
 .31705 0.051831
 .18497 -0.493294
 .01129 -0.749924
 .8675 -0.959329
 .75322 -1.25026
 .60962 -1.37231
 .3216 -1.58832
 .880364 -1.9924
 .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
 .00379863 0.393009
 .0805122 0.201843
 .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
+.12767 1.83675         0
+.441087 1.70526        0
+.639658 1.50546        0
+.806192 1.58667        0
+.966865 1.41151        0
+.24058 1.34506         0
+.58293 1.39453         0
+.02092 1.46038         0
+.31508 0.8539          0
+.31705 0.051831        0
+.18497 -0.493294       0
+.01129 -0.749924       0
+.8675 -0.959329        0
+.75322 -1.25026        0
+.60962 -1.37231        0
+.3216 -1.58832         0
+.880364 -1.9924        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
+.00379863 0.393009     0
+.0805122 0.201843      0
+.0443888 -0.0531952    0
+-0.0199606 -0.175625    0
+-0.0809393 -0.21419     0
 # Low Gain Weights: 6 10
 # (Amplitude)  (Time)
 .0446612 -0.385273
 .038191 -0.00418687
 .0386966 0.0212324
 .0402881 0.0744799
 .0415794 0.229615
 .0598731 0.44332
 .0758477 0.661518
 .101509 1.10641
 .159323 1.64997
 .497256 2.83685
 .245087 3.27499
 .140546 2.46177
 .58086 2.2849
 .632721 2.45587
 .72819 2.52835
 .889583 2.48099
 .980812 2.50031
 .09885 2.55892
 .21374 2.78769
 .61928 3.08069
 .38544 1.95583
 .31998 1.1792
 .50633 0.591226
 .50916 0.0793899
 .5008 -0.33188
 .47339 -0.575386
 .45362 -0.915309
 .40214 -1.31593
 .34175 -1.77904
 .0661 -2.05471
 .31087 -1.49798
 .33793 -1.34758
 .10172 -1.21719
 .08133 -1.09356
 .04007 -0.981455
 .976745 -1.08299
 .930979 -1.14774
 .874203 -1.18348
 .816708 -1.20126
 .587354 -1.92869
 .783078 -1.89621
 .792771 -1.03439
 .622278 -0.781807
 .61184 -0.745831
 .578792 -0.683741
 .537336 -0.596328
 .51443 -0.592858
 .482294 -0.560586
 .462351 -0.827587
 .317989 -1.05649
 .459672 -0.775035
 .468287 -0.619961
 .374182 -0.31635
 .376946 -0.225242
 .367075 -0.347444
 .340737 -0.393231
 .321054 -0.187384
 .320654 -0.225558
 .302148 -0.399499
 .232954 -0.607578
+.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        0
+.09885 2.55892         0
+.21374 2.78769         0
+.61928 3.08069         0
+.38544 1.95583         0
+.31998 1.1792          0
+.50633 0.591226        0
+.50916 0.0793899       0
+.5008 -0.33188         0
+.47339 -0.575386       0
+.45362 -0.915309       0
+.40214 -1.31593        0
+.34175 -1.77904        0
+.0661 -2.05471         0
+.31087 -1.49798        0
+.33793 -1.34758        0
+.10172 -1.21719        0
+.08133 -1.09356        0
+.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

Context Navigation

Legend:

trunk/MagicSoft/Mars/Makefile

trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc

trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc

trunk/MagicSoft/Mars/mjobs/MJCalibration.cc

trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

trunk/MagicSoft/Mars/mpedestal/MPedestalSubtract.cc

trunk/MagicSoft/Mars/mpedestal/Makefile

trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h

trunk/MagicSoft/Mars/msignal/MExtractTime.cc

trunk/MagicSoft/Mars/msignal/MExtractTime.h

trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc

trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

trunk/MagicSoft/Mars/msignal/MExtractedSignalPix.h

trunk/MagicSoft/Mars/msignal/MExtractor.cc

trunk/MagicSoft/Mars/msignal/MExtractor.h

trunk/MagicSoft/Mars/msignal/Makefile

trunk/MagicSoft/Mars/msignal/calibration_weights_UV46.dat

trunk/MagicSoft/Mars/msignal/cosmics_weights46.dat

Download in other formats: