Public Member Functions | Private Member Functions | Private Attributes | List of all members
WireCell::SigProc::PeakFinding Class Reference

#include <PeakFinding.h>

Public Member Functions

 PeakFinding (int fMaxPeaks=200, double sigma=1, double threshold=0.05, bool backgroundRemove=false, int deconIterations=3, bool markov=true, int averWindow=3)
 
 ~PeakFinding ()
 
int find_peak (Waveform::realseq_t &signal)
 
void Clear ()
 
int GetNPeaks ()
 
double * GetPositionX ()
 
double * GetPositionY ()
 

Private Member Functions

int SearchHighRes ()
 

Private Attributes

int fMaxPeaks
 
double sigma
 
double threshold
 
bool backgroundRemove
 
int deconIterations
 
bool markov
 
int averWindow
 
double * source
 
int ssize
 
double * destVector
 
double * fPositionX
 
double * fPositionY
 
int npeaks
 
Log::logptr_t log
 

Detailed Description

Definition at line 10 of file PeakFinding.h.

Constructor & Destructor Documentation

PeakFinding::PeakFinding ( int  fMaxPeaks = 200,
double  sigma = 1,
double  threshold = 0.05,
bool  backgroundRemove = false,
int  deconIterations = 3,
bool  markov = true,
int  averWindow = 3 
)

Definition at line 9 of file PeakFinding.cxx.

13  : fMaxPeaks(fMaxPeaks)
14  , sigma(sigma)
15  , threshold(threshold)
16  , backgroundRemove(backgroundRemove)
17  , deconIterations(deconIterations)
18  , markov(markov)
19  , averWindow(averWindow)
20  , log(Log::logger("sigproc"))
21 {
22 }
WireCell::Log::logger
logptr_t logger(std::string name)
Definition: Logging.cxx:71
PeakFinding::~PeakFinding ( )

Definition at line 24 of file PeakFinding.cxx.

24  {
25  Clear();
26 }

Member Function Documentation

void PeakFinding::Clear ( void  )

Definition at line 28 of file PeakFinding.cxx.

28  {
29 
30 
31 
32  delete [] source;
33  delete [] destVector;
34  delete [] fPositionX;
35  delete [] fPositionY;
36 }
int PeakFinding::find_peak ( Waveform::realseq_t &  signal)

Definition at line 38 of file PeakFinding.cxx.

38  {
39  ssize = int(signal.size());
40  source = new double[ssize];
41  for (size_t i = 0; i!=signal.size();i++){
42  *(source+i) = signal.at(i);
43  }
44 
45  destVector = new double[ssize];
46  fPositionX = new double[ssize];
47  fPositionY = new double[ssize];
48 
49  if (ssize >= 2 * sigma + 1){ // limit the size ...
50  npeaks = SearchHighRes();
51  }else{
52  npeaks = 0;
53  }
54 
55  //fill fPositionY
56  for (int i=0; i<npeaks; i++){
57  int peak_pos1 = std::round(fPositionX[i]);
58  fPositionY[i] = source[peak_pos1];
59  }
60 
61 
62  return npeaks;
63 }
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
int WireCell::SigProc::PeakFinding::GetNPeaks ( )
inline

Definition at line 22 of file PeakFinding.h.

22 {return npeaks;};
double* WireCell::SigProc::PeakFinding::GetPositionX ( )
inline

Definition at line 23 of file PeakFinding.h.

23 {return fPositionX;};
double* WireCell::SigProc::PeakFinding::GetPositionY ( )
inline

Definition at line 24 of file PeakFinding.h.

24 {return fPositionY;};
int PeakFinding::SearchHighRes ( )
private

Definition at line 67 of file PeakFinding.cxx.

68 {
69  const int PEAK_WINDOW = 1024;
70 
71  int i=0, j=0, numberIterations = (int)(7 * sigma + 0.5);
72  double a=0, b=0, c=0;
73  int k=0, lindex=0, posit=0, imin=0, imax=0, jmin=0, jmax=0, lh_gold=0, priz=0;
74  double lda=0, ldb=0, ldc=0, area=0, maximum=0, maximum_decon=0;
75  int xmin=0, xmax=0, l=0, peak_index = 0, w=0;
76  int size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2;
77  double maxch=0;
78  double nom=0, nip=0, nim=0, sp=0, sm=0, plocha = 0;
79  double m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow=0,av=0,men=0;
80  if (sigma < 1) {
81  log->error("SearchHighRes: invalid sigma {}, must be greater than or equal to 1", sigma);
82  return 0;
83  }
84 
85  if(threshold<=0 || threshold>=100){
86  log->error("SearchHighRes: invalid threshold {}, must be positive and less than 100", threshold);
87  return 0;
88  }
89 
90  j = (int) (5.0 * sigma + 0.5);
91  if (j >= PEAK_WINDOW / 2) {
92  log->error("SearchHighRes: too large sigma: j={}",j);
93  return 0;
94  }
95 
96  if (markov == true) {
97  if (averWindow <= 0) {
98  log->error("SearchHighRes: averaging window must be positive, got {}", averWindow);
99  return 0;
100  }
101  }
102 
103  if(backgroundRemove == true){
104  if(ssize < 2 * numberIterations + 1){
105  log->error("SearchHighRes: too large clipping window: {}", ssize);
106  return 0;
107  }
108  }
109 
110  k = int(2 * sigma+0.5);
111  if(k >= 2){
112  for(i = 0;i < k;i++){
113  a = i,b = source[i];
114  m0low += 1,m1low += a,m2low += a * a,l0low += b,l1low += a * b;
115  }
116  detlow = m0low * m2low - m1low * m1low;
117  if(detlow != 0)
118  l1low = (-l0low * m1low + l1low * m0low) / detlow;
119 
120  else
121  l1low = 0;
122  if(l1low > 0)
123  l1low=0;
124  }
125 
126  else{
127  l1low = 0;
128  }
129 
130  i = (int)(7 * sigma + 0.5);
131  i = 2 * i;
132  double *working_space = new double [7 * (ssize + i)];
133  for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
134  for(i = 0; i < size_ext; i++){
135  if(i < shift){
136  a = i - shift;
137  working_space[i + size_ext] = source[0] + l1low * a;
138  if(working_space[i + size_ext] < 0)
139  working_space[i + size_ext]=0;
140  }
141 
142  else if(i >= ssize + shift){
143  a = i - (ssize - 1 + shift);
144  working_space[i + size_ext] = source[ssize - 1];
145  if(working_space[i + size_ext] < 0)
146  working_space[i + size_ext]=0;
147  }
148 
149  else
150  working_space[i + size_ext] = source[i - shift];
151  }
152 
153  if(backgroundRemove == true){
154  for(i = 1; i <= numberIterations; i++){
155  for(j = i; j < size_ext - i; j++){
156  if(markov == false){
157  a = working_space[size_ext + j];
158  b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
159  if(b < a)
160  a = b;
161 
162  working_space[j]=a;
163  }
164 
165  else{
166  a = working_space[size_ext + j];
167  av = 0;
168  men = 0;
169  for (w = j - bw; w <= j + bw; w++){
170  if ( w >= 0 && w < size_ext){
171  av += working_space[size_ext + w];
172  men +=1;
173  }
174  }
175  av = av / men;
176  b = 0;
177  men = 0;
178  for (w = j - i - bw; w <= j - i + bw; w++){
179  if ( w >= 0 && w < size_ext){
180  b += working_space[size_ext + w];
181  men +=1;
182  }
183  }
184  b = b / men;
185  c = 0;
186  men = 0;
187  for (w = j + i - bw; w <= j + i + bw; w++){
188  if ( w >= 0 && w < size_ext){
189  c += working_space[size_ext + w];
190  men +=1;
191  }
192  }
193  c = c / men;
194  b = (b + c) / 2;
195  if (b < a)
196  av = b;
197  working_space[j]=av;
198  }
199  }
200  for(j = i; j < size_ext - i; j++)
201  working_space[size_ext + j] = working_space[j];
202  }
203  for(j = 0;j < size_ext; j++){
204  if(j < shift){
205  a = j - shift;
206  b = source[0] + l1low * a;
207  if (b < 0) b = 0;
208  working_space[size_ext + j] = b - working_space[size_ext + j];
209  }
210 
211  else if(j >= ssize + shift){
212  a = j - (ssize - 1 + shift);
213  b = source[ssize - 1];
214  if (b < 0) b = 0;
215  working_space[size_ext + j] = b - working_space[size_ext + j];
216  }
217 
218  else{
219  working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
220  }
221  }
222  for(j = 0;j < size_ext; j++){
223  if(working_space[size_ext + j] < 0) working_space[size_ext + j] = 0;
224  }
225  }
226 
227  for(i = 0; i < size_ext; i++){
228  working_space[i + 6*size_ext] = working_space[i + size_ext];
229  }
230 
231  if(markov == true){
232  for(j = 0; j < size_ext; j++)
233  working_space[2 * size_ext + j] = working_space[size_ext + j];
234  xmin = 0,xmax = size_ext - 1;
235  for(i = 0, maxch = 0; i < size_ext; i++){
236  working_space[i] = 0;
237  if(maxch < working_space[2 * size_ext + i])
238  maxch = working_space[2 * size_ext + i];
239  plocha += working_space[2 * size_ext + i];
240  }
241  if(maxch == 0) {
242  delete [] working_space;
243  return 0;
244  }
245 
246  nom = 1;
247  working_space[xmin] = 1;
248  for(i = xmin; i < xmax; i++){
249  nip = working_space[2 * size_ext + i] / maxch;
250  nim = working_space[2 * size_ext + i + 1] / maxch;
251  sp = 0,sm = 0;
252  for(l = 1; l <= averWindow; l++){
253  if((i + l) > xmax)
254  a = working_space[2 * size_ext + xmax] / maxch;
255 
256  else
257  a = working_space[2 * size_ext + i + l] / maxch;
258 
259  b = a - nip;
260  if(a + nip <= 0)
261  a=1;
262 
263  else
264  a = sqrt(a + nip);
265 
266  b = b / a;
267  b = exp(b);
268  sp = sp + b;
269  if((i - l + 1) < xmin)
270  a = working_space[2 * size_ext + xmin] / maxch;
271 
272  else
273  a = working_space[2 * size_ext + i - l + 1] / maxch;
274 
275  b = a - nim;
276  if(a + nim <= 0)
277  a = 1;
278 
279  else
280  a = sqrt(a + nim);
281 
282  b = b / a;
283  b = exp(b);
284  sm = sm + b;
285  }
286  a = sp / sm;
287  a = working_space[i + 1] = working_space[i] * a;
288  nom = nom + a;
289  }
290  for(i = xmin; i <= xmax; i++){
291  working_space[i] = working_space[i] / nom;
292  }
293  for(j = 0; j < size_ext; j++)
294  working_space[size_ext + j] = working_space[j] * plocha;
295  for(j = 0; j < size_ext; j++){
296  working_space[2 * size_ext + j] = working_space[size_ext + j];
297  }
298  if(backgroundRemove == true){
299  for(i = 1; i <= numberIterations; i++){
300  for(j = i; j < size_ext - i; j++){
301  a = working_space[size_ext + j];
302  b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
303  if(b < a)
304  a = b;
305  working_space[j] = a;
306  }
307  for(j = i; j < size_ext - i; j++)
308  working_space[size_ext + j] = working_space[j];
309  }
310  for(j = 0; j < size_ext; j++){
311  working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
312  }
313  }
314  }
315 //deconvolution starts
316  area = 0;
317  lh_gold = -1;
318  posit = 0;
319  maximum = 0;
320 //generate response vector
321  for(i = 0; i < size_ext; i++){
322  lda = (double)i - 3 * sigma;
323  lda = lda * lda / (2 * sigma * sigma);
324  j = (int)(1000 * exp(-lda));
325  lda = j;
326  if(lda != 0)
327  lh_gold = i + 1;
328 
329  working_space[i] = lda;
330  area = area + lda;
331  if(lda > maximum){
332  maximum = lda;
333  posit = i;
334  }
335  }
336 //read source vector
337  for(i = 0; i < size_ext; i++)
338  working_space[2 * size_ext + i] = fabs(working_space[size_ext + i]);
339 //create matrix at*a(vector b)
340  i = lh_gold - 1;
341  if(i > size_ext)
342  i = size_ext;
343 
344  imin = -i,imax = i;
345  for(i = imin; i <= imax; i++){
346  lda = 0;
347  jmin = 0;
348  if(i < 0)
349  jmin = -i;
350  jmax = lh_gold - 1 - i;
351  if(jmax > (lh_gold - 1))
352  jmax = lh_gold - 1;
353 
354  for(j = jmin;j <= jmax; j++){
355  ldb = working_space[j];
356  ldc = working_space[i + j];
357  lda = lda + ldb * ldc;
358  }
359  working_space[size_ext + i - imin] = lda;
360  }
361 //create vector p
362  i = lh_gold - 1;
363  imin = -i,imax = size_ext + i - 1;
364  for(i = imin; i <= imax; i++){
365  lda = 0;
366  for(j = 0; j <= (lh_gold - 1); j++){
367  ldb = working_space[j];
368  k = i + j;
369  if(k >= 0 && k < size_ext){
370  ldc = working_space[2 * size_ext + k];
371  lda = lda + ldb * ldc;
372  }
373 
374  }
375  working_space[4 * size_ext + i - imin] = lda;
376  }
377 //move vector p
378  for(i = imin; i <= imax; i++)
379  working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
380 //initialization of resulting vector
381  for(i = 0; i < size_ext; i++)
382  working_space[i] = 1;
383 //START OF ITERATIONS
384  for(lindex = 0; lindex < deconIterations; lindex++){
385  for(i = 0; i < size_ext; i++){
386  if(fabs(working_space[2 * size_ext + i]) > 0.00001 && fabs(working_space[i]) > 0.00001){
387  lda=0;
388  jmin = lh_gold - 1;
389  if(jmin > i)
390  jmin = i;
391 
392  jmin = -jmin;
393  jmax = lh_gold - 1;
394  if(jmax > (size_ext - 1 - i))
395  jmax=size_ext-1-i;
396 
397  for(j = jmin; j <= jmax; j++){
398  ldb = working_space[j + lh_gold - 1 + size_ext];
399  ldc = working_space[i + j];
400  lda = lda + ldb * ldc;
401  }
402  ldb = working_space[2 * size_ext + i];
403  if(lda != 0)
404  lda = ldb / lda;
405 
406  else
407  lda = 0;
408 
409  ldb = working_space[i];
410  lda = lda * ldb;
411  working_space[3 * size_ext + i] = lda;
412  }
413  }
414  for(i = 0; i < size_ext; i++){
415  working_space[i] = working_space[3 * size_ext + i];
416  }
417  }
418 //shift resulting spectrum
419  for(i=0;i<size_ext;i++){
420  lda = working_space[i];
421  j = i + posit;
422  j = j % size_ext;
423  working_space[size_ext + j] = lda;
424  }
425 //write back resulting spectrum
426  maximum = 0, maximum_decon = 0;
427  j = lh_gold - 1;
428  for(i = 0; i < size_ext - j; i++){
429  if(i >= shift && i < ssize + shift){
430  working_space[i] = area * working_space[size_ext + i + j];
431  if(maximum_decon < working_space[i])
432  maximum_decon = working_space[i];
433  if(maximum < working_space[6 * size_ext + i])
434  maximum = working_space[6 * size_ext + i];
435  }
436 
437  else
438  working_space[i] = 0;
439  }
440  lda=1;
441  if(lda>threshold)
442  lda=threshold;
443  lda=lda/100;
444 
445 //searching for peaks in deconvolved spectrum
446  for(i = 1; i < size_ext - 1; i++){
447  if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
448  if(i >= shift && i < ssize + shift){
449  if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
450  for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
451  a += (double)(j - shift) * working_space[j];
452  b += working_space[j];
453  }
454  a = a / b;
455  if(a < 0)
456  a = 0;
457 
458  if(a >= ssize)
459  a = ssize - 1;
460  if(peak_index == 0){
461  fPositionX[0] = a;
462  peak_index = 1;
463  }
464 
465  else{
466  for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
467  if(working_space[6 * size_ext + shift + (int)a] > working_space[6 * size_ext + shift + (int)fPositionX[j]])
468  priz = 1;
469  }
470  if(priz == 0){
471  if(j < fMaxPeaks){
472  fPositionX[j] = a;
473  }
474  }
475 
476  else{
477  for(k = peak_index; k >= j; k--){
478  if(k < fMaxPeaks){
479  fPositionX[k] = fPositionX[k - 1];
480  }
481  }
482  fPositionX[j - 1] = a;
483  }
484  if(peak_index < fMaxPeaks)
485  peak_index += 1;
486  }
487  }
488  }
489  }
490  }
491 
492  for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
493  delete [] working_space;
494  // fNPeaks = peak_index;
495  // if(peak_index == fMaxPeaks){
496  // log->debug("SearchHighRes: peak buffer full with {} peaks, nticks: {}, theshold: {}",
497  // fMaxPeaks, ssize, threshold);
498 
499  // }
500  return peak_index;
501 }
l
static QStrList * l
Definition: config.cpp:1044
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
a
const GenericPointer< typename T::ValueType > T2 T::AllocatorType & a
Definition: pointer.h:1124
b
static bool * b
Definition: config.cpp:1043
maximum
float maximum(const std::vector< float > &wave)
Definition: test_misconfigure.cxx:38

Member Data Documentation

int WireCell::SigProc::PeakFinding::averWindow
private

Definition at line 34 of file PeakFinding.h.

bool WireCell::SigProc::PeakFinding::backgroundRemove
private

Definition at line 31 of file PeakFinding.h.

int WireCell::SigProc::PeakFinding::deconIterations
private

Definition at line 32 of file PeakFinding.h.

double* WireCell::SigProc::PeakFinding::destVector
private

Definition at line 40 of file PeakFinding.h.

int WireCell::SigProc::PeakFinding::fMaxPeaks
private

Definition at line 24 of file PeakFinding.h.

double* WireCell::SigProc::PeakFinding::fPositionX
private

Definition at line 41 of file PeakFinding.h.

double* WireCell::SigProc::PeakFinding::fPositionY
private

Definition at line 42 of file PeakFinding.h.

Log::logptr_t WireCell::SigProc::PeakFinding::log
private

Definition at line 49 of file PeakFinding.h.

bool WireCell::SigProc::PeakFinding::markov
private

Definition at line 33 of file PeakFinding.h.

int WireCell::SigProc::PeakFinding::npeaks
private

Definition at line 44 of file PeakFinding.h.

double WireCell::SigProc::PeakFinding::sigma
private

Definition at line 29 of file PeakFinding.h.

double* WireCell::SigProc::PeakFinding::source
private

Definition at line 37 of file PeakFinding.h.

int WireCell::SigProc::PeakFinding::ssize
private

Definition at line 38 of file PeakFinding.h.

double WireCell::SigProc::PeakFinding::threshold
private

Definition at line 30 of file PeakFinding.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11