FairRoot/PandaRoot
PndEmcFadcFilter.cxx
Go to the documentation of this file.
1 //-----------------------------------------------------------
2 //
3 // Description:
4 // FADC FIR (finite impulse response) filter
5 //
6 // 1) Parameters of the filter can be set directly with SetData() method
7 // one using one of the predefined filters (trapezoidal - SetupTrapez(),
8 // triangle - SetupTriangle(), bipolar triangle - SetupBipolarTriangle() etc.)
9 // 2) Two cosiquently applied filters can be convoluted with Convolute() method
10 // 3) Filtration is performed with Filter() method
11 // 4) Example of Moving Average Filter (average by 3 points)
12 // PndEmcFadcFilter *flt1 = new PndEmcFadcFilter();
13 // Double_t data[3]={1./3.,1./3.,1./3.};
14 // flt1->SetData(data,3,0);
15 //
16 // Author Dima Melnychuk
17 // (based on original code by A.Mykulyak)
18 //-----------------------------------------------------------
19 
20 #include <math.h>
21 #include <algorithm>
22 #include <iostream>
23 
24 #include "PndEmcFadcFilter.h"
25 #include "PndEmcAbsPulseshape.h"
26 
27 using namespace std;
28 
30  fCoeff(),
31  fOffset(0),
32  fType(arbitrary),
33  fIntegerize(kFALSE),
34  fNormFactor(1.),
35  fShiftCount(0)
36 {
37  fCoeff.clear();
38 }
39 
41 {}
42 
43 // Perform filtration of the input signal
44 void PndEmcFadcFilter::Filter(const std::vector<Double_t>& in, std::vector<Double_t>& out) const
45 {
46  Int_t i_size = in.size();
47  if (i_size == 0) {
48  out.clear();
49  return;
50  }
51 
52  Int_t i_nlpad = -MinIndex();
53  Int_t i_nrpad = MaxIndex();
54  Int_t i_ntab = fCoeff.size();
55  std::vector<Double_t> temp(i_size + i_nlpad + i_nrpad);
56 
57  Double_t* p_data = &temp[0];
58  for (Int_t i = 0; i < i_nlpad; i++) *p_data++ = in[0];
59  for (Int_t i = 0; i < i_size; i++) *p_data++ = in[i];
60  for (Int_t i = 0; i < i_nrpad; i++) *p_data++ = in[i_size-1];
61 
62  out.resize(i_size);
63 
64  p_data = &temp[0];
65  for (Int_t i = 0; i < i_size; i++) {
66  const Double_t* p_d = p_data;
67  const Double_t* p_c = &fCoeff[0];
68  Double_t d_sum = 0.;
69 
70  for (Int_t j = 0; j < i_ntab; j++)
71  d_sum += *p_c++ * *p_d++;
72  p_data ++;
73 
74  d_sum *= fNormFactor;
75  if (fIntegerize) d_sum = floor(d_sum);
76  out[i] = d_sum;
77  }
78 
79  return;
80 }
81 
82 // Convolute this filter with a filt .
84 {
85  PndEmcFadcFilter tmp(*this);
86  Convolute(tmp, filt);
87  return;
88 }
89 
90 // Convolute f1 with f2
92  const PndEmcFadcFilter& rhs)
93 {
94  const Double_t* p_lhs = &lhs.fCoeff[0];
95  const Double_t* p_rhs = &rhs.fCoeff[0];
96  Int_t i_nlhs = lhs.Size();
97  Int_t i_nrhs = rhs.Size();
98  Int_t i_size = i_nlhs + i_nrhs - 1;
99 
100  Clear();
101  fCoeff.resize(i_size);
102 
103  for (Int_t i = 0; i < i_size; i++) {
104  Double_t f_sum = 0.;
105  for (Int_t j = 0; j < i_nlhs; j++) {
106  Int_t k = j - i + i_nrhs - 1;
107  Double_t f_rhs = (k >= 0 && k < i_nrhs) ? p_rhs[k] : 0.;
108  f_sum += p_lhs[j] * f_rhs;
109  }
110  fCoeff[i] = f_sum;
111  }
112 
113  fType = lhs.fType * rhs.fType;
114  if (fType == arbitrary) {
115  fOffset = (i_size)/2; // ???
116  } else {
117  fOffset = (i_size)/2;
118  }
119 
120  fIntegerize = lhs.fIntegerize & rhs.fIntegerize;
121  fShiftCount = (fIntegerize) ? lhs.fShiftCount + rhs.fShiftCount : 0;
122  fNormFactor = lhs.fNormFactor * rhs.fNormFactor;
123 
124  return;
125 }
126 
127 // Set parameters of the filter
128 void PndEmcFadcFilter::SetData(Double_t data[], Int_t i_size, Int_t i_offset)
129 {
130  Clear();
131  fCoeff.resize(i_size);
132  for (Int_t i = 0; i < i_size; i++) fCoeff[i] = data[i];
133 
134  if (i_size == 1) {
135  fType = symmetric;
136  } else {
137  Int_t i_nsym = 0;
138  Int_t i_nasym = 0;
139  Int_t i_size2 = i_size/2;
140  for (Int_t i = 0; i < i_size2; i++) {
141  Double_t f_left = data[i];
142  Double_t f_right = data[i_size-i-1];
143  if (f_left == f_right) i_nsym += 1;
144  if (f_left == -f_right) i_nasym += 1;
145  }
146  fType = arbitrary;
147  if (i_nsym == i_size2) fType = symmetric;
148  if (i_nasym == i_size2) fType = antisymmetric;
149  }
150 
151  if (fType != arbitrary || i_offset < 0)
152  i_offset = i_size/2;
153 
154  fOffset = i_offset;
155 
156  return;
157 }
158 
159 // Setup moving average filter
160 void PndEmcFadcFilter::SetupMA(Int_t i_width)
161 {
162  Int_t i_size = i_width;
163  Clear();
164  fCoeff.resize(i_size);
165  fType = symmetric;
166 
167  for (Int_t i = 0; i < i_size; i++) {
168  fCoeff[i] = 1./i_width;
169  }
170 
171  return;
172 }
173 
174 // Setup moving window deconvolution filter
175 void PndEmcFadcFilter::SetupMWD(Int_t i_width, Double_t tau)
176 {
177  Int_t i_size = i_width;
178  Clear();
179  fCoeff.resize(i_size);
180  fType = arbitrary;
181 
182  fCoeff[0]=1;
183  fCoeff[i_size-1]=-(1-1./tau);
184 
185  for (Int_t i = 1; i < i_size-1; i++) {
186  fCoeff[i] = 1./tau;
187  }
188 
189  return;
190 }
191 
192 // Matched digital filter
193 void PndEmcFadcFilter::SetupMatchedFilter(Int_t i_width, PndEmcAbsPulseshape *pulseshape, Double_t sampleRate)
194 {
195  Int_t i_size = i_width;
196  Clear();
197  fCoeff.resize(i_size);
198  fType = arbitrary;
199 
200  Double_t amplitude=1.;
201  Double_t t;
202  Double_t val=0;
203 
204  for ( int i=0;i<i_size;i++)
205  {
206  t= i/sampleRate;
207  val=pulseshape->value(t,amplitude,0.);
208  fCoeff[i]=val;
209  }
210 
211  std::reverse(fCoeff.begin(), fCoeff.end());
212 
213  return;
214 
215 }
216 
217 void PndEmcFadcFilter::SetupBipolarTrapez(Int_t i_rise, Int_t i_flat, Int_t i_width)
218 {
219  Int_t i_size = 2*i_rise + 2*i_flat+i_width+1;
220  Clear();
221 
222  fCoeff.resize(i_size);
223  fType = symmetric;
224 
225  for (Int_t i = 0; i < i_rise; i++) {
226  fCoeff[i] = i;
227  fCoeff[i_size-i-1] = -i;
228  }
229 
230  for (Int_t i = i_rise; i < (i_rise+i_flat); i++)
231  {
232  fCoeff[i] = i_rise;
233  fCoeff[i_flat+i_width+i+1] = -i_rise;
234  }
235 
236  for (Int_t i = (i_rise+i_flat); i <= (i_rise+i_flat+i_width); i++) {
237  fCoeff[i] = i_rise-2*i_rise/i_width*(i-i_rise-i_width);
238  }
239 
240  return;
241 }
242 
243 // Setup filter with trapezoidal weighting function.
244 void PndEmcFadcFilter::SetupTrapez(Int_t i_rise, Int_t i_flat)
245 {
246  Int_t i_size = 2*i_rise + i_flat;
247  Clear();
248 
249  fOffset = (i_size)/2;
250  fCoeff.resize(i_size);
251  fType = symmetric;
252 
253  for (Int_t i = 0; i < i_size; i++) fCoeff[i] = i_rise+1;
254  for (Int_t i = 0; i < i_rise; i++) {
255  fCoeff[i] = i+1;
256  fCoeff[i_size-i-1] = i+1;
257  }
258 
259  return;
260 }
261 
262 
263 // Setup filter with bipolar triangular weighting function.
265 {
266  Int_t i_size = 4*i_rise + 1;
267  Clear();
268 
269  fOffset = 2*i_rise;
270  fCoeff.resize(i_size);
272 
273  for (Int_t i = 0; i <= i_rise; i++) {
274  fCoeff[i] = i;
275  fCoeff[fOffset-i] = i;
276  fCoeff[fOffset+i] = -i;
277  fCoeff[i_size-i-1] = -i;
278  }
279 
280  return;
281 }
282 
283 // Setup a differentiator
284 void PndEmcFadcFilter::SetupDifferentiator(Int_t i_lag, Int_t i_width)
285 {
286  Int_t i_size = i_lag + 2*i_width;
287  Clear();
288 
289  fOffset = (i_size)/2;
290  fCoeff.resize(i_size);
292 
293  for (Int_t i = 0; i < i_size; i++) fCoeff[i] = 0;
294  for (Int_t i = 0; i < i_width; i++) {
295  fCoeff[i] = -1.;
296  fCoeff[i_size-i-1] = 1.;
297  }
298 
299  return;
300 }
301 
302 // Setup a 2nd order differentiator.
303 void PndEmcFadcFilter::SetupDoubleDifferentiator(Int_t i_npos, Int_t i_nneg,
304  Int_t i_nzero)
305 {
306  Int_t i_size = i_npos + 2*(i_nzero + i_nneg);
307  Clear();
308 
309  fOffset = (i_size)/2;
310  fCoeff.resize(i_size);
311  fType = symmetric;
312 
313  Int_t i_apos = 2*i_nneg;
314  Int_t i_aneg = i_npos;
315  Int_t i_div = 1;
316  for (Int_t i = 2; i < max(i_apos,i_aneg); i++) {
317  if (i_apos%i == 0 && i_aneg%i == 0) i_div = i;
318  }
319  i_apos /= i_div;
320  i_aneg /= i_div;
321 
322  for (Int_t i = 0; i < i_size; i++) fCoeff[i] = i_apos;
323  for (Int_t i = 0; i < i_nneg; i++) {
324  fCoeff[i] = -i_aneg;
325  fCoeff[i_size-i-1] = -i_aneg;
326  }
327  for (Int_t i = i_nneg; i < i_nneg+i_nzero; i++) {
328  fCoeff[i] = 0.;
329  fCoeff[i_size-i-1] = 0.;
330  }
331 
332  return;
333 }
334 
335 // Setup a differentiator with Pole/Zero cancelation.
337 {
338  Int_t i_size = i_lag + 2;
339  Clear();
340 
341  fOffset = (i_size)/2;
342  fCoeff.resize(i_size);
343  fType = arbitrary;
344 
345  for (Int_t i = 0; i < i_size; i++) fCoeff[i] = 0.;
346  fCoeff[0] = -d_fac;
347  fCoeff[i_size-1] = 1.;
348 
349  return;
350 }
351 
352 // Use floating normalization with factor a d_norm .
354 {
355  fIntegerize = kFALSE;
356  fShiftCount = 0;
357  fNormFactor = d_norm;
358  return;
359 }
360 
361 // Use integer normalization with a shift count of a i_shift .
363 {
364  fIntegerize = kTRUE;
365  fShiftCount = i_shift;
366  fNormFactor = 1.;
367  for (Int_t i = 0; i < i_shift; i++) fNormFactor *= 0.5;
368  return;
369 }
370 
371 // Returns the lowest index.
373 {
374  return -fOffset;
375 }
376 
377 // Returns the highest index.
379 {
380  return fCoeff.size() - 1 - fOffset;
381 }
382 
383 // Returns size of filter (number of tabs).
385 {
386  return fCoeff.size();
387 }
388 
390 {
391  fCoeff.clear();
392  fOffset = 0;
393  fType = arbitrary;
394  return;
395 }
396 
void SetupPZDifferentiator(Int_t i_lag=0, Double_t d_fac=1.)
void SetupMA(Int_t i_width)
virtual void Filter(const std::vector< Double_t > &in, std::vector< Double_t > &out) const
Int_t i
Definition: run_full.C:25
void SetNormalizeFloating(Double_t d_norm=1.)
void SetupBipolarTriangle(Int_t i_rise)
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
void SetupBipolarTrapez(Int_t i_rise, Int_t i_flat, Int_t i_width)
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
FADC FIR (finite impulse response) filter.
void SetupMWD(Int_t i_width, Double_t tau)
void Convolute(const PndEmcFadcFilter &filt)
virtual double value(const double t, const double amp, const double toffset) const
void SetupMatchedFilter(Int_t i_width, PndEmcAbsPulseshape *pulseshape, Double_t sampleRate)
Double_t
Int_t MinIndex() const
void SetupDifferentiator(Int_t i_lag=0, Int_t i_width=1)
TFile * out
Definition: reco_muo.C:20
virtual ~PndEmcFadcFilter()
Int_t Size() const
void SetData(Double_t data[], Int_t i_size, Int_t i_offset)
void SetupDoubleDifferentiator(Int_t i_npos=1, Int_t i_nneg=1, Int_t i_nzero=0)
pulseshape interface
std::vector< Double_t > fCoeff
TTree * t
Definition: bump_analys.C:13
Int_t MaxIndex() const
void SetupTrapez(Int_t i_rise, Int_t i_flat)
ClassImp(PndEmcFadcFilter)
void SetNormalizeInteger(Int_t i_shift=0)