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