| 1 | #ifndef MARS_MExtralgoDigitalFilter
|
|---|
| 2 | #define MARS_MExtralgoDigitalFilter
|
|---|
| 3 |
|
|---|
| 4 | #ifndef ROOT_TROOT
|
|---|
| 5 | #include <TROOT.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 | Float_t *fVal;
|
|---|
| 19 | Int_t fNum;
|
|---|
| 20 |
|
|---|
| 21 | Float_t const *fWeightsAmp;
|
|---|
| 22 | Float_t const *fWeightsTime;
|
|---|
| 23 |
|
|---|
| 24 | const Int_t fWeightsPerBin; // Number of weights per data bin
|
|---|
| 25 | const Int_t fWindowSize;
|
|---|
| 26 |
|
|---|
| 27 | // Result
|
|---|
| 28 | Float_t fTime;
|
|---|
| 29 | Float_t fTimeDev;
|
|---|
| 30 | Float_t fSignal;
|
|---|
| 31 | Float_t fSignalDev;
|
|---|
| 32 |
|
|---|
| 33 | // Weights: Weights to evaluate
|
|---|
| 34 | // Startv: Index of first bin of data
|
|---|
| 35 | // startw: Offset on the weights
|
|---|
| 36 | inline Double_t Eval(Float_t const *weights, const Int_t startv, const Int_t startw=0) const
|
|---|
| 37 | {
|
|---|
| 38 | //
|
|---|
| 39 | // Slide with a window of size windowsize over the sample
|
|---|
| 40 | // and multiply the entries with the corresponding weights
|
|---|
| 41 | //
|
|---|
| 42 | Double_t sum = 0;
|
|---|
| 43 |
|
|---|
| 44 | // Shift the start of the weight to the center of sample 0
|
|---|
| 45 | Float_t const *w = weights + startw;
|
|---|
| 46 |
|
|---|
| 47 | Float_t *const beg = fVal+startv;
|
|---|
| 48 | for (Float_t const *pex=beg; pex<beg+fWindowSize; pex++)
|
|---|
| 49 | {
|
|---|
| 50 | sum += *w * *pex;
|
|---|
| 51 | w += fWeightsPerBin;
|
|---|
| 52 | }
|
|---|
| 53 | return sum;
|
|---|
| 54 | }
|
|---|
| 55 |
|
|---|
| 56 | inline void AlignIntoLimits(Int_t &maxp, Int_t &frac) const
|
|---|
| 57 | {
|
|---|
| 58 | // Align maxp into available range (TO BE CHECKED)
|
|---|
| 59 | if (maxp < 0)
|
|---|
| 60 | {
|
|---|
| 61 | maxp = 0;
|
|---|
| 62 | frac = fWeightsPerBin/2-1; // Assume peak at the end of the last slice
|
|---|
| 63 | }
|
|---|
| 64 | if (maxp > fNum-fWindowSize)
|
|---|
| 65 | {
|
|---|
| 66 | maxp = fNum-fWindowSize;
|
|---|
| 67 | frac = -fWeightsPerBin/2; // Assume peak at the beginning of the first slice
|
|---|
| 68 | }
|
|---|
| 69 | }
|
|---|
| 70 |
|
|---|
| 71 | inline Int_t AlignExtractionWindow(Int_t &maxp, Int_t &frac, const Double_t ampsum)
|
|---|
| 72 | {
|
|---|
| 73 | // Align extraction window to signal position
|
|---|
| 74 |
|
|---|
| 75 | const Double_t timesum = Eval(fWeightsTime, maxp, frac);
|
|---|
| 76 |
|
|---|
| 77 | // Because fWeightsPerBin/2 doesn't correspond to the center
|
|---|
| 78 | // of a bin the time-values extracted are slightly positive.
|
|---|
| 79 | // They are roughly between -0.45 and 0.55
|
|---|
| 80 | const Double_t binoffset = TMath::Even(fWeightsPerBin) ? 0.5 : 0;
|
|---|
| 81 |
|
|---|
| 82 | // This is the time offset from the extraction position
|
|---|
| 83 | Double_t tmoffset = (frac+binoffset)/fWeightsPerBin + timesum/ampsum;
|
|---|
| 84 |
|
|---|
| 85 | // Convert the residual fraction of one slice into an
|
|---|
| 86 | // offset position in the extraction weights
|
|---|
| 87 | const Int_t integ = TMath::FloorNint(tmoffset+0.5);
|
|---|
| 88 |
|
|---|
| 89 | /*
|
|---|
| 90 | if (integ>0)
|
|---|
| 91 | tmoffset=0.49-0.05;
|
|---|
| 92 | if (integ<0)
|
|---|
| 93 | tmoffset=-0.49-0.05;
|
|---|
| 94 | integ=0;
|
|---|
| 95 | */
|
|---|
| 96 |
|
|---|
| 97 | // move the extractor by an offset number of slices
|
|---|
| 98 | // determined by the extracted time
|
|---|
| 99 | maxp -= integ;
|
|---|
| 100 |
|
|---|
| 101 | frac = TMath::FloorNint((tmoffset-integ)*fWeightsPerBin);
|
|---|
| 102 |
|
|---|
| 103 | // Align maxp into available range (TO BE CHECKED)
|
|---|
| 104 | AlignIntoLimits(maxp, frac);
|
|---|
| 105 |
|
|---|
| 106 | return integ;
|
|---|
| 107 | }
|
|---|
| 108 |
|
|---|
| 109 | inline void AlignExtractionWindow(Int_t &maxp, Int_t &frac)
|
|---|
| 110 | {
|
|---|
| 111 | const Double_t amp = Eval(fWeightsAmp, maxp, frac);
|
|---|
| 112 | if (amp!=0)
|
|---|
| 113 | AlignExtractionWindow(maxp, frac, amp);
|
|---|
| 114 | }
|
|---|
| 115 |
|
|---|
| 116 | public:
|
|---|
| 117 | MExtralgoDigitalFilter(Int_t res, Int_t windowsize, Float_t *wa, Float_t *wt)
|
|---|
| 118 | : fVal(0), fNum(0), fWeightsAmp(wa+res/2), fWeightsTime(wt+res/2),
|
|---|
| 119 | fWeightsPerBin(res), fWindowSize(windowsize),
|
|---|
| 120 | fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
|
|---|
| 121 | {
|
|---|
| 122 | }
|
|---|
| 123 |
|
|---|
| 124 | void SetData(Int_t n, Float_t *val) { fNum=n; fVal=val; }
|
|---|
| 125 |
|
|---|
| 126 | Float_t GetTime() const { return fTime; }
|
|---|
| 127 | Float_t GetSignal() const { return fSignal; }
|
|---|
| 128 |
|
|---|
| 129 | Float_t GetTimeDev() const { return fTimeDev; }
|
|---|
| 130 | Float_t GetSignalDev() const { return fSignalDev; }
|
|---|
| 131 |
|
|---|
| 132 | void GetSignal(Float_t &sig, Float_t &dsig) const { sig=fSignal; dsig=fSignalDev; }
|
|---|
| 133 | void GetTime(Float_t &sig, Float_t &dsig) const { sig=fTime; dsig=fTimeDev; }
|
|---|
| 134 |
|
|---|
| 135 | Float_t ExtractNoise(Int_t iter) const;
|
|---|
| 136 | void Extract();
|
|---|
| 137 |
|
|---|
| 138 | static Bool_t CalculateWeights(TH1 &shape, const TH2 &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
|
|---|
| 139 | static void CalculateWeights2(TH1F &shape, const TH2F &autocorr, TArrayF &wa, TArrayF &wt, Int_t wpb=-1);
|
|---|
| 140 | };
|
|---|
| 141 |
|
|---|
| 142 |
|
|---|
| 143 | #endif
|
|---|
