GSI Object Oriented Online Offline (Go4)  GO4-5.3.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
f_find_peaks.c
Go to the documentation of this file.
1 // $Id: f_find_peaks.c 1956 2016-11-18 16:46:37Z linev $
2 //-----------------------------------------------------------------------
3 // The GSI Online Offline Object Oriented (Go4) Project
4 // Experiment Data Processing at EE department, GSI
5 //-----------------------------------------------------------------------
6 // Copyright (C) 2000- GSI Helmholtzzentrum f�r Schwerionenforschung GmbH
7 // Planckstr. 1, 64291 Darmstadt, Germany
8 // Contact: http://go4.gsi.de
9 //-----------------------------------------------------------------------
10 // This software can be used under the license agreements as stated
11 // in Go4License.txt file which is part of the distribution.
12 //-----------------------------------------------------------------------
13 
14 
15 /****************************************************************
16  *
17  * Copyright (C)
18  *
19  * Gesellschaft f. Schwerionenforschung, GSI
20  * Planckstr. 1
21  * 64291 Darmstadt
22  * Germany
23  *
24  * Author: H.Essel(at)gsi.de
25  *
26  * This program is free software; you can redistribute it and/or
27  * modify it under the terms of the
28  * GNU General Public License
29  * as published by the Free Software Foundation; either version 2
30  * of the License, or (at your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful,
33  * but WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
35  *
36  * See the GNU General Public License for more details
37  * (http://www.gnu.org).
38  *
39  *****************************************************************/
40 
41 #include <stdio.h>
42 #include <stdlib.h>
43 #include <string.h>
44 #include <math.h>
45 
46 #define UP 1
47 #define DOWN -1
48 #define HORI 0
49 #define PRINT 0
50 #include "f_find_peaks.h"
51 
52 int f_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum)
53 {
54 
55 #define LOW 2
56 #define WIDTH 8
57 
58  double r_sig_f ;
59  int I,J,K,L ;
60  double d_sum_prod ;
61  double l_max_chan ;
62  double *pl_data,d_max ;
63 
64  /* 2.0E2*SQRT(2.0E0*LOG(2.0E0))/100. */;
65  r_sig_f = 55.4518;
66  r_sig_f = 2.;
67  *pr_sum = 0;
68  d_sum_prod = 0;
69  *pr_sig = 0.;
70  *pr_pos = 0.;
71  d_max=0.0;
72 
73  /* get maximum channel content and integral */
74  /* pl_data = pa_data;
75  for(J = 0; J < l_len; J++)
76  {
77  if(d_max < *pl_data)
78  {
79  d_max = *pl_data;
80  l_max_chan = J;
81  }
82  pl_data++;
83  }
84  */
85  /* Calculate first momentum */
86  pl_data = pa_data;
87  for(J = 1; J <= l_len; J++)
88  {
89  *pr_sum = *pr_sum + *pl_data;
90  d_sum_prod = d_sum_prod + J * *pl_data;
91  pl_data++;
92  }
93  /* Calculate second momentum */
94  if(*pr_sum > 0)
95  {
96  *pr_pos = (d_sum_prod/ *pr_sum + 0.5);
97  pl_data = pa_data;
98  for(J = 1; J <= l_len; J++)
99  {
100  *pr_sig = *pr_sig + ((double)J - *pr_pos) * ((double)J - *pr_pos) * *pl_data;
101  pl_data++;
102  }
103  *pr_sig = r_sig_f * sqrt(*pr_sig/ *pr_sum);
104  }
105  else return(1);
106  *pr_pos = *pr_pos -1.0;
107  return(0);
108 }
109 
110 //**********************************************************************************
111 
112 void f_find_peaks(
113  void *pfData, // pointer to data
114  int lType, // data type: 0=w 1=l 2=f 3=d
115  int lFirstChan, // first channel
116  int lPoints, // number of points
117  int lAver, // channels to sum up
118  double dDeltaFact, // noise factor
119  double dDeltaMin, // noise minimum
120  double *dNoise, // noise array
121  int lPeaksMax, // Maximum number of peaks
122  int *plPeaks, // returns number of peaks found
123  int *plMinima, // pointer to list of minima
124  double *pdMinima, // pointer to list of minima values
125  int *plMaxima, // pointer to list of maxima
126  int *plWidth, // pointer to list of width
127  double *pdIntegral) // pointer to list of integrals
128 {
129  short *pw;
130  int lMaxIndex,lMinIndex,lIndexMax,lDirection,lIndex,lLastUp,lLastDown,ii,i,
131  *plMinimaL,*plMaximaL,*plMinimaR,*plMaximaR,lPrint,*pl;
132  float *pf;
133  double dValue,dLow,dHigh,dDelta,dAver,*pdData,*pdNoise,dMax,dMin,dPos,dSig,dSum,*pd;
134 
135  *plPeaks=0;
136  pdNoise=NULL;
137  pdData = (double *)malloc(lPoints*8);
138  pdNoise = (double *)malloc(lPoints*8);
139  memset(pdData,0,lPoints*8);
140  memset(pdNoise,0,lPoints*8);
141  lPrint=PRINT;
142  dMax=0.;
143  dMin=1e20;
144  pw = ((short *)pfData) + lFirstChan;
145  pl = ((int *)pfData) + lFirstChan;
146  pf = ((float *)pfData) + lFirstChan;
147  pd = ((double *)pfData) + lFirstChan;
148  if(lAver == 0)lAver=1;
149  lPoints=lPoints/lAver;
150  dAver = (double) lAver;
151  // copy the data field into local double field
152  switch(lType)
153  {
154  case TYPE__SHORT: for(ii=0;ii<lPoints;ii++)
155  { *(pdData+ii) = (double)*(pw+lAver*ii);
156  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pw+lAver*ii+i);
157  *(pdData+ii) = *(pdData+ii)/dAver;
158  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
159  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
160  case TYPE__INT: for(ii=0;ii<lPoints;ii++)
161  { *(pdData+ii) = (double)*(pl+lAver*ii);
162  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pl+lAver*ii+i);
163  *(pdData+ii) = *(pdData+ii)/dAver;
164  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
165  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
166  case TYPE__FLOAT: for(ii=0;ii<lPoints;ii++)
167  { *(pdData+ii) = (double)*(pf+lAver*ii);
168  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pf+lAver*ii+i);
169  *(pdData+ii) = *(pdData+ii)/dAver;
170  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
171  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
172  case TYPE__DOUBLE: for(ii=0;ii<lPoints;ii++)
173  { *(pdData+ii) = (double)*(pd+lAver*ii);
174  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pd+lAver*ii+i);
175  *(pdData+ii) = *(pdData+ii)/dAver;
176  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
177  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
178  default: printf("Data type %d not supported!\n",lType);
179  printf("Type 0: short, 1: int, 2: float, 3: double\n");
180  free(pdData);
181  free(pdNoise);
182  return;
183  }
184  if(dNoise == NULL)
185  {// calculate from square root
186  for(ii=0;ii<lPoints;ii++)
187  {
188  if(*(pdData+ii) == 0.0) *(pdNoise+ii) = dDeltaFact;
189  else *(pdNoise+ii) = sqrt(*(pdData+ii)) * dDeltaFact;
190  if(dDeltaMin > *(pdNoise+ii)) *(pdNoise+ii)=dDeltaMin;
191  }
192  }else {// sum up bins from noise input array
193  for(ii=0;ii<lPoints;ii++)
194  { *(pdNoise+ii) = *(dNoise+lFirstChan+lAver*ii);
195  for(i=1;i<lAver;i++)*(pdNoise+ii) += *(dNoise+lFirstChan+lAver*ii+i);
196  *(pdNoise+ii) = *(pdNoise+ii)/dAver;
197  }
198  }
199  lMaxIndex=0;
200  lMinIndex=0;
201  lIndexMax=lPoints-1;
202  plMinimaL = (int *)malloc(lPeaksMax*16);
203  plMaximaL = (int *)(plMinimaL+lPeaksMax);
204  plMinimaR = (int *)(plMinimaL+2*lPeaksMax);
205  plMaximaR = (int *)(plMinimaL+3*lPeaksMax);
206  memset(plMinimaL,0,lPeaksMax*16);
207 
208  dDelta= *pdNoise;
209  dHigh = *pdData + dDelta;
210  dLow = *pdData - dDelta;
211  //-----------------------------------------------------
212  // loop: get points of minima and maxima
213  lDirection=HORI;
214  lIndex=0;
215  lLastDown=0;
216  lLastUp=0;
217  while ((lIndex <= lIndexMax) & (lMaxIndex < lPeaksMax-2))
218  {
219  dValue = *(pdData+lIndex);
220  if(dValue > dHigh)
221  {// direction goes upstairs
222  if (lDirection != UP)
223  { // direction was downstairs: minimum found
224  if(lPrint)printf("%5d - %5d - %5d d %8.1f\n",
225  lMinIndex,
226  lLastDown,
227  lIndex,
228  dDelta);
229  *(plMinimaR+lMinIndex)=lIndex;
230  *(plMinimaL+lMinIndex)=lLastDown;
231  lMinIndex++;
232  }
233  lLastUp=lIndex;
234  lDirection=UP;// direction goes upstairs
235  }
236  if(dValue < dLow)
237  {// direction goes downstairs
238  if (lDirection == UP)
239  { // direction was upstairs: maximum found
240  if(lPrint)printf("%5d + %5d - %5d d %8.1f\n",
241  lMaxIndex,
242  lLastUp,
243  lIndex,
244  dDelta);
245  *(plMaximaR+lMaxIndex)=lIndex;
246  *(plMaximaL+lMaxIndex)=lLastUp;
247  lMaxIndex++;
248  }
249  lLastDown=lIndex;
250  lDirection=DOWN;// direction goes downstairs
251  }
252  // calculate new limits if noise band has been left
253  if((dValue < dLow)|(dValue > dHigh))
254  {
255  dDelta= *(pdNoise+lIndex);
256  dHigh = dValue + dDelta;
257  dLow = dValue - dDelta;
258  }
259  lIndex++;
260  }
261  if(lDirection == DOWN)
262  {
263  *(plMinimaL+lMinIndex)=lLastDown;
264  *(plMinimaR+lMinIndex)=lIndexMax;
265  if(lPrint)printf("%5d - %5d - %5d\n",
266  lMinIndex,
267  lLastDown,
268  lIndexMax);
269  lMinIndex++;
270  }
271  //-----------------------------------------------------
272  if(lPrint)printf("-------------------------\n");
273  if(lPrint)printf("Minima %d Maxima %d \n",lMinIndex,lMaxIndex);
274  // for each peak there must be a minimum left and right
275  if(lMinIndex != (lMaxIndex+1))
276  {
277  free(plMinimaL);
278  free(pdData);
279  free(pdNoise);
280  printf("Error, wrong minima\n");
281  return;
282  }
283  // calculate mean values of minima L&R, mean content, peaks
284  // start with first minimum
285  if(lPrint)printf(" [ Left Right ] \n");
286  *plMinima = (*plMinimaL + *plMinimaR)/2;
287  *plMinima = (int)((double)(*plMinima)*dAver + dAver/2.)+lFirstChan;
288  *pdMinima=0.; // average content at mean minimum channel
289  for(ii = *plMinimaL;ii <= *plMinimaR;ii++) *pdMinima += *(pdData+ii);
290  *pdMinima = *pdMinima/(double)(*plMinimaR - *plMinimaL + 1);
291  if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
292  *plMinimaL,
293  *plMinimaR,
294  *plMinima,
295  *pdMinima);
296  for(lIndex=0;lIndex < lMaxIndex;lIndex++)
297  {
298  // simple mean peak value, overwritten later
299  *(plMaxima+lIndex)=(*(plMaximaL+lIndex)+*(plMaximaR+lIndex))/2;
300  if(lPrint)printf("+ [%6d %6d] %6d\n",
301  *(plMaximaL+lIndex),
302  *(plMaximaR+lIndex),
303  *(plMaxima+lIndex));
304  // mean minima values (channel and content)
305  *(plMinima+lIndex+1) = (*(plMinimaL+lIndex+1) + *(plMinimaR+lIndex+1))/2;
306  *(plMinima+lIndex+1) = (int)((double)(*(plMinima+lIndex+1))*dAver + dAver/2.)+lFirstChan;
307  *(pdMinima+lIndex+1)=0.;
308  for(ii = *(plMinimaL+lIndex+1);ii <= *(plMinimaR+lIndex+1);ii++)
309  *(pdMinima+lIndex+1) += *(pdData+ii);
310  *(pdMinima+lIndex+1) = *(pdMinima+lIndex+1)/
311  (double)(*(plMinimaR+lIndex+1) - *(plMinimaL+lIndex+1)+1);
312  if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
313  *(plMinimaL+lIndex+1),
314  *(plMinimaR+lIndex+1),
315  *(plMinima+lIndex+1),
316  *(pdMinima+lIndex+1));
317  // Peak positions, width and sum by momenta
318  ii = *(plMinimaL+lIndex+1) - *(plMinimaR+lIndex) + 1;
319  if(ii > 1)
320  {
321  f_position(ii,pdData + *(plMinimaR+lIndex),&dPos,&dSig,&dSum);
322 
323  *(plMaxima+*plPeaks)=(int)((dPos+(double)(*(plMinimaR+lIndex)))*dAver)+lFirstChan;
324  *(plWidth+*plPeaks) =(int)(dSig*dAver + 0.5);
325  if(*(pdMinima+lIndex) < *(pdMinima+lIndex+1))
326  dSum=dSum-((double)ii * *(pdMinima+lIndex)+ (double)ii * (*(pdMinima+lIndex+1)-*(pdMinima+lIndex))/2.);
327  else
328  dSum=dSum-((double)ii * *(pdMinima+lIndex)- (double)ii * (*(pdMinima+lIndex)-*(pdMinima+lIndex+1))/2.);
329  *(pdIntegral+*plPeaks) = dSum*dAver;
330  (*plPeaks)++;
331  }
332  else
333  {
334  free(plMinimaL);
335  free(pdData);
336  free(pdNoise);
337  printf("Error, negative interval\n");
338  return;
339  }
340  }
341  if(lPrint)printf("peaks found %d\n",*plPeaks);
342 
343  free(plMinimaL);
344  free(pdData);
345  free(pdNoise);
346 }
TYPE__FLOAT
#define TYPE__FLOAT
Definition: f_find_peaks.c:18
PRINT
#define PRINT
Definition: f_find_peaks.c:49
TYPE__INT
#define TYPE__INT
Definition: f_find_peaks.c:17
HORI
#define HORI
Definition: f_find_peaks.c:48
DOWN
#define DOWN
Definition: f_find_peaks.c:47
f_find_peaks.h
TYPE__DOUBLE
#define TYPE__DOUBLE
Definition: f_find_peaks.c:19
f_position
int f_position(int l_len, double *pa_data, double *pr_pos, double *pr_sig, double *pr_sum)
Definition: f_find_peaks.c:52
TYPE__SHORT
#define TYPE__SHORT
Definition: f_find_peaks.c:16
UP
#define UP
Definition: f_find_peaks.c:46
f_find_peaks
void f_find_peaks(void *pfData, int lType, int lFirstChan, int lPoints, int lAver, double dDeltaFact, double dDeltaMin, double *dNoise, int lPeaksMax, int *plPeaks, int *plMinima, double *pdMinima, int *plMaxima, int *plWidth, double *pdIntegral)
Definition: f_find_peaks.c:112
Generated on Fri Oct 26 2018 08:34:51 for GSI Object Oriented Online Offline (Go4) by   doxygen 1.8.8