| 1 | #ifndef MARS_MExtralgoDigitalFilter
|
|---|
| 2 | #define MARS_MExtralgoDigitalFilter
|
|---|
| 3 |
|
|---|
| 4 | #ifndef ROOT_TMatrix
|
|---|
| 5 | #include <TMatrix.h>
|
|---|
| 6 | #endif
|
|---|
| 7 |
|
|---|
| 8 | class TH1;
|
|---|
| 9 | class TH2;
|
|---|
| 10 | class TH1F;
|
|---|
| 11 | class TH2F;
|
|---|
| 12 | class TArrayF;
|
|---|
| 13 |
|
|---|
| 14 | class MExtralgoDigitalFilter
|
|---|
| 15 | {
|
|---|
| 16 | private:
|
|---|
| 17 | // Input
|
|---|
| 18 | const Float_t *fVal;
|
|---|
| 19 | Int_t fNum;
|
|---|
| 20 |
|
|---|
| 21 | Float_t const *fWeightsAmp;
|
|---|
| 22 | Float_t const *fWeightsTime;
|
|---|
| 23 | Float_t const *fPulseShape;
|
|---|
| 24 |
|
|---|
| 25 | const TMatrix *fAinv;
|
|---|
| 26 |
|
|---|
| 27 | const Int_t fWeightsPerBin; // Number of weights per data bin
|
|---|
| 28 | const Int_t fWindowSize;
|
|---|
| 29 |
|
|---|
| 30 | // Result
|
|---|
| 31 | Float_t fTime;
|
|---|
| 32 | Float_t fTimeDev;
|
|---|
| 33 | Float_t fSignal;
|
|---|
| 34 | Float_t fSignalDev;
|
|---|
| 35 |
|
|---|
| 36 | Float_t GetChisq(const Int_t maxp, const Int_t frac, const Float_t sum) const;
|
|---|
| 37 |
|
|---|
| 38 | inline Double_t ChiSq(const Double_t sum, const Int_t startv, const Int_t startw=0) const
|
|---|
| 39 | {
|
|---|
| 40 | //
|
|---|
| 41 | // Slide with a window of size windowsize over the sample
|
|---|
| 42 | // and multiply the entries with the corresponding weights
|
|---|
| 43 | //
|
|---|
| 44 | Double_t chisq = 0;
|
|---|
| 45 |
|
|---|
| 46 | // Shift the start of the weight to the center of sample 0
|
|---|
| 47 | Float_t const *w = fPulseShape + startw;
|
|---|
| 48 |
|
|---|
| 49 | const Float_t *beg = fVal+startv;
|
|---|
| 50 | for (Float_t const *pex=beg; pex<beg+fWindowSize; pex++)
|
|---|
| 51 | {
|
|---|
| 52 | const Double_t c = *w - *pex/sum;
|
|---|
| 53 | chisq += c*c;
|
|---|
| 54 | w += fWeightsPerBin;
|
|---|
| 55 | }
|
|---|
| 56 | return chisq;
|
|---|
| 57 | }
|
|---|
| 58 |
|
|---|
| 59 | // Weights: Weights to evaluate
|
|---|
| 60 | // Startv: Index of first bin of data
|
|---|
| 61 | // startw: Offset on the weights
|
|---|
| 62 | inline Double_t Eval(Float_t const *weights, const Int_t startv, const Int_t startw=0) const
|
|---|
| 63 | {
|
|---|
| 64 | //
|
|---|
| 65 | // Slide with a window of size windowsize over the sample
|
|---|
| 66 | // and multiply the entries with the corresponding weights
|
|---|
| 67 | //
|
|---|
| 68 | Double_t sum = 0;
|
|---|
| 69 |
|
|---|
| 70 | // Shift the start of the weight to the center of sample 0
|
|---|
| 71 | Float_t const *w = weights + startw;
|
|---|
| 72 |
|
|---|
| 73 | const Float_t *beg = fVal+startv;
|
|---|
| 74 | for (Float_t const *pex=beg; pex<beg+fWindowSize; pex++)
|
|---|
| 75 | {
|
|---|
| 76 | sum += *w * *pex;
|
|---|
| 77 | w += fWeightsPerBin;
|
|---|
| 78 | }
|
|---|
| 79 | return sum;
|
|---|
| 80 | }
|
|---|
| 81 |
|
|---|
| 82 | inline void AlignIntoLimits(Int_t &maxp, Int_t &frac) const
|
|---|
| 83 | {
|
|---|
| 84 | // Align maxp into available range (TO BE CHECKED)
|
|---|
| 85 | if (maxp < 0)
|
|---|
| 86 | {
|
|---|
| 87 | maxp = 0;
|
|---|
| 88 | frac = fWeightsPerBin/2-1; // Assume peak at the end of the last slice
|
|---|
| 89 | }
|
|---|
| 90 | if (maxp > fNum-fWindowSize)
|
|---|
| 91 | {
|
|---|
| 92 | maxp = fNum-fWindowSize;
|
|---|
| 93 | frac = -fWeightsPerBin/2; // Assume peak at the beginning of the first slice
|
|---|
| 94 | }
|
|---|
| 95 | }
|
|---|
| 96 |
|
|---|
| 97 | Int_t AlignExtractionWindow(Int_t &maxp, Int_t &frac, const Double_t ampsum);
|
|---|
| 98 | void AlignExtractionWindow(Int_t &maxp, Int_t &frac)
|
|---|
| 99 | {
|
|---|
| 100 | const Double_t amp = Eval(fWeightsAmp, maxp, frac);
|
|---|
| 101 | if (amp!=0)
|
|---|
| 102 | AlignExtractionWindow(maxp, frac, amp);
|
|---|
| 103 | }
|
|---|
| 104 |
|
|---|
| 105 | public:
|
|---|
| 106 | MExtralgoDigitalFilter(Int_t res, Int_t windowsize, Float_t *wa, Float_t *wt, Float_t *ps=0, TMatrix *ainv=0)
|
|---|
| 107 | : fVal(0), fNum(0), fWeightsAmp(wa+res/2), fWeightsTime(wt+res/2),
|
|---|
| 108 | fPulseShape(ps), fAinv(ainv), fWeightsPerBin(res), fWindowSize(windowsize),
|
|---|
| 109 | fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
|
|---|
| 110 | {
|
|---|
| 111 | }
|
|---|
| 112 |
|
|---|
| 113 | void SetData(Int_t n, Float_t const *val) { fNum=n; fVal=val; }
|
|---|
| 114 |
|
|---|
| 115 | Float_t GetTime() const { return fTime; }
|
|---|
| 116 | Float_t GetSignal() const { return fSignal; }
|
|---|
| 117 |
|
|---|
| 118 | Float_t GetTimeDev() const { return fTimeDev; }
|
|---|
| 119 | Float_t GetSignalDev() const { return fSignalDev; }
|
|---|
| 120 |
|
|---|
| 121 | void GetSignal(Float_t &sig, Float_t &dsig) const { sig=fSignal; dsig=fSignalDev; }
|
|---|
| 122 | void GetTime(Float_t &sig, Float_t &dsig) const { sig=fTime; dsig=fTimeDev; }
|
|---|
| 123 |
|
|---|
| 124 | Float_t ExtractNoise() const;
|
|---|
| 125 | void Extract(Int_t maxpos=-1);
|
|---|
| 126 |
|
|---|
| 127 | static Int_t CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
|
|---|
| 128 | static Int_t CalculateWeights2(TH1 &shape, const TH2 &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
|
|---|
| 129 | };
|
|---|
| 130 |
|
|---|
| 131 |
|
|---|
| 132 | #endif
|
|---|
