GSI Object Oriented Online Offline (Go4)  GO4-6.1.4
 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 3064 2021-03-12 15:55:22Z 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 fuer 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 #include <stdlib.h>
15 
16 #define TYPE__SHORT 0
17 #define TYPE__INT 1
18 #define TYPE__FLOAT 2
19 #define TYPE__DOUBLE 3
20 
21 #define UP 1
22 #define DOWN -1
23 #define HORI 0
24 #define PRINT 0
25 
26 int go4fit_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum);
27 
29  void *pfData, // pointer to data
30  int lType, // data type: 0=w 1=l 2=f 3=d
31  int lFirstChan, // first channel
32  int lPoints, // number of points
33  int lAver, // channels to sum up
34  double dDeltaFact, // noise factor
35  double dDeltaMin, // noise minimum
36  double *dNoise, // noise array
37  int lPeaksMax, // Maximum number of peaks
38  int *plPeaks, // returns number of peaks found
39  int *plMinima, // pointer to list of minima
40  double *pdMinima, // pointer to list of minima values
41  int *plMean, // pointer to list of maxima
42  int *plWidth, // pointer to list of width
43  double *pdIntegral, // pointer to list of integrals
44  int *plLeft, // pointer to array of left RMS index
45  int *plMaximum, // pointer to array of maximum index
46  int *plRight, // pointer to array of right RMS index
47  float *pfLowNoise, // pointer to data array of low noise
48  float *pfHighNoise) // pointer to data array of high noise
49 {
50  short *pw;
51  int lMaxIndex,lMinIndex,lIndexMax,lDirection,lIndex,lLastUp,lLastDown,ii,i,lMax,lLeft,lRight,
52  *plMinimaL,*plMeanL,*plMinimaR,*plMeanR,lPrint,*pl;
53  float *pf;
54  double dValue,dLow,dHigh,dDelta,dAver,*pdData,*pdNoise,dMax,dMin,dPos,dSig,dSum,*pd;
55 
56  *plPeaks=0;
57  pdNoise=NULL;
58  plMinimaL = NULL;
59  pdData = (double *)malloc(lPoints*8);
60  pdNoise = (double *)malloc(lPoints*8);
61  if (!pdData || !pdNoise)
62  goto clear_data;
63 
64  memset(pdData,0,lPoints*8);
65  memset(pdNoise,0,lPoints*8);
66  lPrint=PRINT;
67  dMax=0.;
68  dMin=1e20;
69  pw = ((short *)pfData) + lFirstChan;
70  pl = ((int *)pfData) + lFirstChan;
71  pf = ((float *)pfData) + lFirstChan;
72  pd = ((double *)pfData) + lFirstChan;
73  if(lAver == 0)lAver=1;
74  lPoints=lPoints/lAver;
75  dAver = (double) lAver;
76  // copy the data field into local double field
77  switch(lType)
78  {
79  case TYPE__SHORT: for(ii=0;ii<lPoints;ii++)
80  { *(pdData+ii) = (double)*(pw+lAver*ii);
81  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pw+lAver*ii+i);
82  *(pdData+ii) = *(pdData+ii)/dAver;
83  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
84  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
85  case TYPE__INT: for(ii=0;ii<lPoints;ii++)
86  { *(pdData+ii) = (double)*(pl+lAver*ii);
87  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pl+lAver*ii+i);
88  *(pdData+ii) = *(pdData+ii)/dAver;
89  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
90  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
91  case TYPE__FLOAT: for(ii=0;ii<lPoints;ii++)
92  { *(pdData+ii) = (double)*(pf+lAver*ii);
93  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pf+lAver*ii+i);
94  *(pdData+ii) = *(pdData+ii)/dAver;
95  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
96  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
97  case TYPE__DOUBLE: for(ii=0;ii<lPoints;ii++)
98  { *(pdData+ii) = (double)*(pd+lAver*ii);
99  for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pd+lAver*ii+i);
100  *(pdData+ii) = *(pdData+ii)/dAver;
101  if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
102  if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
103  default: printf("Data type %d not supported!\n",lType);
104  printf("Type 0: short, 1: int, 2: float, 3: double\n");
105  free(pdData);
106  free(pdNoise);
107  return;
108  }
109  if(dNoise == NULL)
110  {// calculate from square root
111  for(ii=0;ii<lPoints;ii++)
112  {
113  if(*(pdData+ii) == 0.0) *(pdNoise+ii) = dDeltaFact;
114  else *(pdNoise+ii) = sqrt(*(pdData+ii)) * dDeltaFact;
115  if(dDeltaMin > *(pdNoise+ii)) *(pdNoise+ii)=dDeltaMin;
116  }
117  }else {// sum up bins from noise input array
118  for(ii=0;ii<lPoints;ii++)
119  { *(pdNoise+ii) = *(dNoise+lFirstChan+lAver*ii);
120  for(i=1;i<lAver;i++)*(pdNoise+ii) += *(dNoise+lFirstChan+lAver*ii+i);
121  *(pdNoise+ii) = *(pdNoise+ii)/dAver;
122  }
123  }
124  lMaxIndex=0;
125  lMinIndex=0;
126  lIndexMax=lPoints-1;
127  plMinimaL = (int *)malloc(lPeaksMax*16);
128  if (!plMinimaL)
129  goto clear_data;
130  plMeanL = (int *)(plMinimaL+lPeaksMax);
131  plMinimaR = (int *)(plMinimaL+2*lPeaksMax);
132  plMeanR = (int *)(plMinimaL+3*lPeaksMax);
133  memset(plMinimaL,0,lPeaksMax*16);
134 
135  dDelta= *pdNoise;
136  dHigh = *pdData + dDelta;
137  dLow = *pdData - dDelta;
138  if(pfLowNoise)*pfLowNoise=dLow;
139  if(pfHighNoise)*pfHighNoise=dHigh;
140  //-----------------------------------------------------
141  // loop: get points of minima and maxima
142  lDirection=HORI;
143  lIndex=0;
144  lLastDown=0;
145  lLastUp=0;
146  while ((lIndex <= lIndexMax) & (lMaxIndex < lPeaksMax-2))
147  {
148  dValue = *(pdData+lIndex);
149  if(dValue > dHigh)
150  {// direction goes upstairs
151  if (lDirection != UP)
152  { // direction was downstairs: minimum found
153  if(lPrint)printf("%5d - %5d - %5d d %8.1f\n",
154  lMinIndex,
155  lLastDown,
156  lIndex,
157  dDelta);
158  *(plMinimaR+lMinIndex)=lIndex;
159  *(plMinimaL+lMinIndex)=lLastDown;
160  lMinIndex++;
161  }
162  lLastUp=lIndex;
163  lDirection=UP;// direction goes upstairs
164  }
165  if(dValue < dLow)
166  {// direction goes downstairs
167  if (lDirection == UP)
168  { // direction was upstairs: maximum found
169  if(lPrint)printf("%5d + %5d - %5d d %8.1f\n",
170  lMaxIndex,
171  lLastUp,
172  lIndex,
173  dDelta);
174  *(plMeanR+lMaxIndex)=lIndex;
175  *(plMeanL+lMaxIndex)=lLastUp;
176  lMaxIndex++;
177  }
178  lLastDown=lIndex;
179  lDirection=DOWN;// direction goes downstairs
180  }
181  // calculate new limits if noise band has been left
182  if((dValue < dLow)|(dValue > dHigh))
183  {
184  dDelta= *(pdNoise+lIndex);
185  dHigh = dValue + dDelta;
186  dLow = dValue - dDelta;
187  }
188  lIndex++;
189  if(pfLowNoise)*(pfLowNoise+lIndex)=dLow;
190  if(pfHighNoise)*(pfHighNoise+lIndex)=dHigh;
191  }
192  if(lDirection == DOWN)
193  {
194  *(plMinimaL+lMinIndex)=lLastDown;
195  *(plMinimaR+lMinIndex)=lIndexMax;
196  if(lPrint)printf("%5d - %5d - %5d\n",
197  lMinIndex,
198  lLastDown,
199  lIndexMax);
200  lMinIndex++;
201  }
202  //-----------------------------------------------------
203  if(lPrint)printf("-------------------------\n");
204  if(lPrint)printf("Minima %d Maxima %d \n",lMinIndex,lMaxIndex);
205  // for each peak there must be a minimum left and right
206  if(lMinIndex != (lMaxIndex+1))
207  {
208  free(plMinimaL);
209  free(pdData);
210  free(pdNoise);
211  *plPeaks=0;
212  if(lPrint)printf("Error, wrong minima\n");
213  return;
214  }
215  // calculate mean values of minima L&R, mean content, peaks
216  // start with first minimum
217  if(lPrint)printf(" [ Left Right ] \n");
218  *plMinima = (*plMinimaL + *plMinimaR)/2;
219  *pdMinima=0.; // average content at mean minimum channel
220  for(ii = *plMinimaL;ii <= *plMinimaR;ii++) *pdMinima += *(pdData+ii);
221  *pdMinima = *pdMinima/(double)(*plMinimaR - *plMinimaL + 1);
222  if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
223  *plMinimaL,
224  *plMinimaR,
225  *plMinima,
226  *pdMinima);
227  for(lIndex=0;lIndex < lMaxIndex;lIndex++)
228  {
229  // simple mean peak value, overwritten later
230  *(plMean+lIndex)=(*(plMeanL+lIndex)+*(plMeanR+lIndex))/2;
231  if(lPrint)printf("+ [%6d %6d] %6d\n",
232  *(plMeanL+lIndex),
233  *(plMeanR+lIndex),
234  *(plMean+lIndex));
235  // mean minima values (channel and content)
236  *(plMinima+lIndex+1) = (*(plMinimaL+lIndex+1) + *(plMinimaR+lIndex+1))/2;
237  *(pdMinima+lIndex+1)=0.;
238  for(ii = *(plMinimaL+lIndex+1);ii <= *(plMinimaR+lIndex+1);ii++)
239  *(pdMinima+lIndex+1) += *(pdData+ii);
240  *(pdMinima+lIndex+1) = *(pdMinima+lIndex+1)/(double)(*(plMinimaR+lIndex+1) - *(plMinimaL+lIndex+1)+1);
241  if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
242  *(plMinimaL+lIndex+1),
243  *(plMinimaR+lIndex+1),
244  *(plMinima+lIndex+1),
245  *(pdMinima+lIndex+1));
246  // Peak positions, width and sum by momenta
247  ii = *(plMinimaL+lIndex+1) - *(plMinimaR+lIndex) + 1;
248  if(ii > 1)
249  {
250  go4fit_position(ii,pdData + *(plMinimaR+lIndex),&dPos,&dSig,&dSum);
251 
252  *(plMean+*plPeaks)=(int)((dPos+(double)(*(plMinimaR+lIndex)))*dAver)+lFirstChan;
253  *(plWidth+*plPeaks) =(int)(dSig*dAver + 0.5);
254  if(*(pdMinima+lIndex) < *(pdMinima+lIndex+1))
255  dSum=dSum-((double)ii * *(pdMinima+lIndex)+ (double)ii * (*(pdMinima+lIndex+1)-*(pdMinima+lIndex))/2.);
256  else
257  dSum=dSum-((double)ii * *(pdMinima+lIndex)- (double)ii * (*(pdMinima+lIndex)-*(pdMinima+lIndex+1))/2.);
258  *(pdIntegral+*plPeaks) = dSum*dAver;
259  (*plPeaks)++;
260  }
261  else
262  {
263  free(plMinimaL);
264  free(pdData);
265  free(pdNoise);
266  *plPeaks=0;
267  if(lPrint)printf("Error, negative interval\n");
268  return;
269  }
270  }
271  for(i=0;i<*plPeaks;i++)
272  {
273  dMax=0; // search local maximum in interval
274  dMin=1000000000; // local minimum
275  pd=pdData+*(plMinima+i);
276  for(ii=*(plMinima+i);ii<*(plMinima+i+1);ii++){if(dMax < *pd)dMax = *pd; if(dMin > *pd)dMin = *pd; pd++;}
277  dMin=dMin+(dMax-dMin)/2.; // threashold for RMS
278  pd=pdData+*(plMinima+i); // search channel of local maximum
279  for(lMax=*(plMinima+i);lMax<*(plMinima+i+1);lMax++){if(dMax == *pd)break; pd++;} // l_max is now channel of local maximum
280  lLeft=*(plMinima+i);
281  lRight=*(plMinima+i+1);
282  *(plMaximum+i)=lMax;
283  pd=pdData+lMax;
284  for(ii=lMax;ii<lRight;ii++){if(dMin > *pd){lRight=ii;break;}pd++;}
285  pd=pdData+lMax;
286  for(ii=lMax;ii>lLeft;ii--){if(dMin > *pd){lLeft=ii;break;}pd--;}
287  *(plLeft+i)=lLeft;
288  *(plRight+i)=lRight;
289  }
290  *(plMinima) = (int)((double)(*(plMinima))*dAver + dAver/2.)+lFirstChan;
291  for(i=0;i<*plPeaks;i++)
292  {
293  // calculate channel numbers of original data
294  *(plMaximum+i) = (int)((double)(*(plMaximum+i))*dAver + dAver/2.)+lFirstChan;
295  *(plLeft+i) = (int)((double)(*(plLeft+i))*dAver + dAver/2.)+lFirstChan;
296  *(plRight+i) = (int)((double)(*(plRight+i))*dAver + dAver/2.)+lFirstChan;
297  *(plMinima+i+1) = (int)((double)(*(plMinima+i+1))*dAver + dAver/2.)+lFirstChan;
298  }
299  if(lPrint)printf("peaks found %d\n",*plPeaks);
300 
301 clear_data:
302  free(plMinimaL);
303  free(pdData);
304  free(pdNoise);
305 }
306 //**********************************************************************************
307 int go4fit_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum)
308 {
309 
310 #define LOW 2
311 #define WIDTH 8
312 
313  double r_sig_f;
314  int J;
315  double d_sum_prod ;
316  double *pl_data;
317 
318  /* 2.0E2*SQRT(2.0E0*LOG(2.0E0))/100. */;
319  r_sig_f = 55.4518;
320  r_sig_f = 2.;
321  *pr_sum = 0;
322  d_sum_prod = 0;
323  *pr_sig = 0.;
324  *pr_pos = 0.;
325 
326  /* Calculate first momentum */
327  pl_data = pa_data;
328  for(J = 1; J <= l_len; J++)
329  {
330  *pr_sum = *pr_sum + *pl_data;
331  d_sum_prod = d_sum_prod + J * *pl_data;
332  pl_data++;
333  }
334  /* Calculate second momentum */
335  if(*pr_sum > 0)
336  {
337  *pr_pos = (d_sum_prod/ *pr_sum + 0.5);
338  pl_data = pa_data;
339  for(J = 1; J <= l_len; J++)
340  {
341  *pr_sig = *pr_sig + ((double)J - *pr_pos) * ((double)J - *pr_pos) * *pl_data;
342  pl_data++;
343  }
344  *pr_sig = r_sig_f * sqrt(*pr_sig/ *pr_sum);
345  }
346  else return(1);
347  *pr_pos = *pr_pos -1.0;
348  return(0);
349 }
#define UP
Definition: f_find_peaks.c:21
#define PRINT
Definition: f_find_peaks.c:24
void go4fit_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 *plMean, int *plWidth, double *pdIntegral, int *plLeft, int *plMaximum, int *plRight, float *pfLowNoise, float *pfHighNoise)
Definition: f_find_peaks.c:28
#define TYPE__INT
Definition: f_find_peaks.c:17
#define TYPE__DOUBLE
Definition: f_find_peaks.c:19
int go4fit_position(int l_len, double *pa_data, double *pr_pos, double *pr_sig, double *pr_sum)
Definition: f_find_peaks.c:307
#define DOWN
Definition: f_find_peaks.c:22
#define TYPE__SHORT
Definition: f_find_peaks.c:16
#define HORI
Definition: f_find_peaks.c:23
#define TYPE__FLOAT
Definition: f_find_peaks.c:18