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