MbsAPI/f_find_peaks.c
Go to the documentation of this file.
00001 // $Id: f_find_peaks.c 1196 2014-04-10 14:06:24Z linev $ 00002 //----------------------------------------------------------------------- 00003 // The GSI Online Offline Object Oriented (Go4) Project 00004 // Experiment Data Processing at EE department, GSI 00005 //----------------------------------------------------------------------- 00006 // Copyright (C) 2000- GSI Helmholtzzentrum für Schwerionenforschung GmbH 00007 // Planckstr. 1, 64291 Darmstadt, Germany 00008 // Contact: http://go4.gsi.de 00009 //----------------------------------------------------------------------- 00010 // This software can be used under the license agreements as stated 00011 // in Go4License.txt file which is part of the distribution. 00012 //----------------------------------------------------------------------- 00013 00014 00015 /**************************************************************** 00016 * 00017 * Copyright (C) 00018 * 00019 * Gesellschaft f. Schwerionenforschung, GSI 00020 * Planckstr. 1 00021 * 64291 Darmstadt 00022 * Germany 00023 * 00024 * Author: H.Essel(at)gsi.de 00025 * 00026 * This program is free software; you can redistribute it and/or 00027 * modify it under the terms of the 00028 * GNU General Public License 00029 * as published by the Free Software Foundation; either version 2 00030 * of the License, or (at your option) any later version. 00031 * 00032 * This program is distributed in the hope that it will be useful, 00033 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00034 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 00035 * 00036 * See the GNU General Public License for more details 00037 * (http://www.gnu.org). 00038 * 00039 *****************************************************************/ 00040 00041 #include <stdio.h> 00042 #include <stdlib.h> 00043 #include <string.h> 00044 #include <math.h> 00045 00046 #define UP 1 00047 #define DOWN -1 00048 #define HORI 0 00049 #define PRINT 0 00050 #include "f_find_peaks.h" 00051 00052 void f_find_peaks( 00053 void *pfData, // pointer to data 00054 int lType, // data type: 0=w 1=l 2=f 3=d 00055 int lFirstChan, // first channel 00056 int lPoints, // number of points 00057 int lAver, // channels to sum up 00058 double dDeltaFact, // noise factor 00059 double dDeltaMin, // noise minimum 00060 double *dNoise, // noise array 00061 int lPeaksMax, // Maximum number of peaks 00062 int *plPeaks, // returns number of peaks found 00063 int *plMinima, // pointer to list of minima 00064 double *pdMinima, // pointer to list of minima values 00065 int *plMaxima, // pointer to list of maxima 00066 int *plWidth, // pointer to list of width 00067 double *pdIntegral) // pointer to list of integrals 00068 { 00069 short *pw; 00070 int lMaxIndex,lMinIndex,lIndexMax,lDirection,lIndex,lLastUp,lLastDown,ii,i, 00071 *plMinimaL,*plMaximaL,*plMinimaR,*plMaximaR,lPrint,*pl; 00072 float *pf; 00073 double dValue,dLow,dHigh,dDelta,dAver,*pdData,*pdNoise,dMax,dMin,dPos,dSig,dSum,*pd; 00074 00075 *plPeaks=0; 00076 pdNoise=NULL; 00077 pdData = (double *)malloc(lPoints*8); 00078 pdNoise = (double *)malloc(lPoints*8); 00079 memset(pdData,0,lPoints*8); 00080 memset(pdNoise,0,lPoints*8); 00081 lPrint=PRINT; 00082 dMax=0.; 00083 dMin=1e20; 00084 pw = ((short *)pfData) + lFirstChan; 00085 pl = ((int *)pfData) + lFirstChan; 00086 pf = ((float *)pfData) + lFirstChan; 00087 pd = ((double *)pfData) + lFirstChan; 00088 if(lAver == 0)lAver=1; 00089 lPoints=lPoints/lAver; 00090 dAver = (double) lAver; 00091 // copy the data field into local double field 00092 switch(lType) 00093 { 00094 case TYPE__SHORT: for(ii=0;ii<lPoints;ii++) 00095 { *(pdData+ii) = (double)*(pw+lAver*ii); 00096 for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pw+lAver*ii+i); 00097 *(pdData+ii) = *(pdData+ii)/dAver; 00098 if(*(pdData+ii) > dMax) dMax=*(pdData+ii); 00099 if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break; 00100 case TYPE__INT: for(ii=0;ii<lPoints;ii++) 00101 { *(pdData+ii) = (double)*(pl+lAver*ii); 00102 for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pl+lAver*ii+i); 00103 *(pdData+ii) = *(pdData+ii)/dAver; 00104 if(*(pdData+ii) > dMax) dMax=*(pdData+ii); 00105 if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break; 00106 case TYPE__FLOAT: for(ii=0;ii<lPoints;ii++) 00107 { *(pdData+ii) = (double)*(pf+lAver*ii); 00108 for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pf+lAver*ii+i); 00109 *(pdData+ii) = *(pdData+ii)/dAver; 00110 if(*(pdData+ii) > dMax) dMax=*(pdData+ii); 00111 if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break; 00112 case TYPE__DOUBLE: for(ii=0;ii<lPoints;ii++) 00113 { *(pdData+ii) = (double)*(pd+lAver*ii); 00114 for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pd+lAver*ii+i); 00115 *(pdData+ii) = *(pdData+ii)/dAver; 00116 if(*(pdData+ii) > dMax) dMax=*(pdData+ii); 00117 if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break; 00118 default: printf("Data type %d not supported!\n",lType); 00119 printf("Type 0: short, 1: int, 2: float, 3: double\n"); 00120 free(pdData); 00121 free(pdNoise); 00122 return; 00123 } 00124 if(dNoise == NULL) 00125 {// calculate from square root 00126 for(ii=0;ii<lPoints;ii++) 00127 { 00128 if(*(pdData+ii) == 0.0) *(pdNoise+ii) = dDeltaFact; 00129 else *(pdNoise+ii) = sqrt(*(pdData+ii)) * dDeltaFact; 00130 if(dDeltaMin > *(pdNoise+ii)) *(pdNoise+ii)=dDeltaMin; 00131 } 00132 }else {// sum up bins from noise input array 00133 for(ii=0;ii<lPoints;ii++) 00134 { *(pdNoise+ii) = *(dNoise+lFirstChan+lAver*ii); 00135 for(i=1;i<lAver;i++)*(pdNoise+ii) += *(dNoise+lFirstChan+lAver*ii+i); 00136 *(pdNoise+ii) = *(pdNoise+ii)/dAver; 00137 } 00138 } 00139 lMaxIndex=0; 00140 lMinIndex=0; 00141 lIndexMax=lPoints-1; 00142 plMinimaL = (int *)malloc(lPeaksMax*16); 00143 plMaximaL = (int *)(plMinimaL+lPeaksMax); 00144 plMinimaR = (int *)(plMinimaL+2*lPeaksMax); 00145 plMaximaR = (int *)(plMinimaL+3*lPeaksMax); 00146 memset(plMinimaL,0,lPeaksMax*16); 00147 00148 dDelta= *pdNoise; 00149 dHigh = *pdData + dDelta; 00150 dLow = *pdData - dDelta; 00151 //----------------------------------------------------- 00152 // loop: get points of minima and maxima 00153 lDirection=HORI; 00154 lIndex=0; 00155 lLastDown=0; 00156 lLastUp=0; 00157 while ((lIndex <= lIndexMax) & (lMaxIndex < lPeaksMax-2)) 00158 { 00159 dValue = *(pdData+lIndex); 00160 if(dValue > dHigh) 00161 {// direction goes upstairs 00162 if (lDirection != UP) 00163 { // direction was downstairs: minimum found 00164 if(lPrint)printf("%5d - %5d - %5d d %8.1f\n", 00165 lMinIndex, 00166 lLastDown, 00167 lIndex, 00168 dDelta); 00169 *(plMinimaR+lMinIndex)=lIndex; 00170 *(plMinimaL+lMinIndex)=lLastDown; 00171 lMinIndex++; 00172 } 00173 lLastUp=lIndex; 00174 lDirection=UP;// direction goes upstairs 00175 } 00176 if(dValue < dLow) 00177 {// direction goes downstairs 00178 if (lDirection == UP) 00179 { // direction was upstairs: maximum found 00180 if(lPrint)printf("%5d + %5d - %5d d %8.1f\n", 00181 lMaxIndex, 00182 lLastUp, 00183 lIndex, 00184 dDelta); 00185 *(plMaximaR+lMaxIndex)=lIndex; 00186 *(plMaximaL+lMaxIndex)=lLastUp; 00187 lMaxIndex++; 00188 } 00189 lLastDown=lIndex; 00190 lDirection=DOWN;// direction goes downstairs 00191 } 00192 // calculate new limits if noise band has been left 00193 if((dValue < dLow)|(dValue > dHigh)) 00194 { 00195 dDelta= *(pdNoise+lIndex); 00196 dHigh = dValue + dDelta; 00197 dLow = dValue - dDelta; 00198 } 00199 lIndex++; 00200 } 00201 if(lDirection == DOWN) 00202 { 00203 *(plMinimaL+lMinIndex)=lLastDown; 00204 *(plMinimaR+lMinIndex)=lIndexMax; 00205 if(lPrint)printf("%5d - %5d - %5d\n", 00206 lMinIndex, 00207 lLastDown, 00208 lIndexMax); 00209 lMinIndex++; 00210 } 00211 //----------------------------------------------------- 00212 if(lPrint)printf("-------------------------\n"); 00213 if(lPrint)printf("Minima %d Maxima %d \n",lMinIndex,lMaxIndex); 00214 // for each peak there must be a minimum left and right 00215 if(lMinIndex != (lMaxIndex+1)) 00216 { 00217 free(plMinimaL); 00218 free(pdData); 00219 free(pdNoise); 00220 printf("Error, wrong minima\n"); 00221 return; 00222 } 00223 // calculate mean values of minima L&R, mean content, peaks 00224 // start with first minimum 00225 if(lPrint)printf(" [ Left Right ] \n"); 00226 *plMinima = (*plMinimaL + *plMinimaR)/2; 00227 *plMinima = (int)((double)(*plMinima)*dAver + dAver/2.)+lFirstChan; 00228 *pdMinima=0.; // average content at mean minimum channel 00229 for(ii = *plMinimaL;ii <= *plMinimaR;ii++) *pdMinima += *(pdData+ii); 00230 *pdMinima = *pdMinima/(double)(*plMinimaR - *plMinimaL + 1); 00231 if(lPrint)printf("- [%6d %6d] %6d %6.1f\n", 00232 *plMinimaL, 00233 *plMinimaR, 00234 *plMinima, 00235 *pdMinima); 00236 for(lIndex=0;lIndex < lMaxIndex;lIndex++) 00237 { 00238 // simple mean peak value, overwritten later 00239 *(plMaxima+lIndex)=(*(plMaximaL+lIndex)+*(plMaximaR+lIndex))/2; 00240 if(lPrint)printf("+ [%6d %6d] %6d\n", 00241 *(plMaximaL+lIndex), 00242 *(plMaximaR+lIndex), 00243 *(plMaxima+lIndex)); 00244 // mean minima values (channel and content) 00245 *(plMinima+lIndex+1) = (*(plMinimaL+lIndex+1) + *(plMinimaR+lIndex+1))/2; 00246 *(plMinima+lIndex+1) = (int)((double)(*(plMinima+lIndex+1))*dAver + dAver/2.)+lFirstChan; 00247 *(pdMinima+lIndex+1)=0.; 00248 for(ii = *(plMinimaL+lIndex+1);ii <= *(plMinimaR+lIndex+1);ii++) 00249 *(pdMinima+lIndex+1) += *(pdData+ii); 00250 *(pdMinima+lIndex+1) = *(pdMinima+lIndex+1)/ 00251 (double)(*(plMinimaR+lIndex+1) - *(plMinimaL+lIndex+1)+1); 00252 if(lPrint)printf("- [%6d %6d] %6d %6.1f\n", 00253 *(plMinimaL+lIndex+1), 00254 *(plMinimaR+lIndex+1), 00255 *(plMinima+lIndex+1), 00256 *(pdMinima+lIndex+1)); 00257 // Peak positions, width and sum by momenta 00258 ii = *(plMinimaL+lIndex+1) - *(plMinimaR+lIndex) + 1; 00259 if(ii > 1) 00260 { 00261 f_position(ii,pdData + *(plMinimaR+lIndex),&dPos,&dSig,&dSum); 00262 00263 *(plMaxima+*plPeaks)=(int)((dPos+(double)(*(plMinimaR+lIndex)))*dAver)+lFirstChan; 00264 *(plWidth+*plPeaks) =(int)(dSig*dAver + 0.5); 00265 if(*(pdMinima+lIndex) < *(pdMinima+lIndex+1)) 00266 dSum=dSum-((double)ii * *(pdMinima+lIndex)+ (double)ii * (*(pdMinima+lIndex+1)-*(pdMinima+lIndex))/2.); 00267 else 00268 dSum=dSum-((double)ii * *(pdMinima+lIndex)- (double)ii * (*(pdMinima+lIndex)-*(pdMinima+lIndex+1))/2.); 00269 *(pdIntegral+*plPeaks) = dSum*dAver; 00270 (*plPeaks)++; 00271 } 00272 else 00273 { 00274 free(plMinimaL); 00275 free(pdData); 00276 free(pdNoise); 00277 printf("Error, negative interval\n"); 00278 return; 00279 } 00280 } 00281 if(lPrint)printf("peaks found %d\n",*plPeaks); 00282 00283 free(plMinimaL); 00284 free(pdData); 00285 free(pdNoise); 00286 } 00287 //********************************************************************************** 00288 int f_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum) 00289 { 00290 00291 #define LOW 2 00292 #define WIDTH 8 00293 00294 double r_sig_f ; 00295 int I,J,K,L ; 00296 double d_sum_prod ; 00297 double l_max_chan ; 00298 double *pl_data,d_max ; 00299 00300 /* 2.0E2*SQRT(2.0E0*LOG(2.0E0))/100. */; 00301 r_sig_f = 55.4518; 00302 r_sig_f = 2.; 00303 *pr_sum = 0; 00304 d_sum_prod = 0; 00305 *pr_sig = 0.; 00306 *pr_pos = 0.; 00307 d_max=0.0; 00308 00309 /* get maximum channel content and integral */ 00310 /* pl_data = pa_data; 00311 for(J = 0; J < l_len; J++) 00312 { 00313 if(d_max < *pl_data) 00314 { 00315 d_max = *pl_data; 00316 l_max_chan = J; 00317 } 00318 pl_data++; 00319 } 00320 */ 00321 /* Calculate first momentum */ 00322 pl_data = pa_data; 00323 for(J = 1; J <= l_len; J++) 00324 { 00325 *pr_sum = *pr_sum + *pl_data; 00326 d_sum_prod = d_sum_prod + J * *pl_data; 00327 pl_data++; 00328 } 00329 /* Calculate second momentum */ 00330 if(*pr_sum > 0) 00331 { 00332 *pr_pos = (d_sum_prod/ *pr_sum + 0.5); 00333 pl_data = pa_data; 00334 for(J = 1; J <= l_len; J++) 00335 { 00336 *pr_sig = *pr_sig + ((double)J - *pr_pos) * ((double)J - *pr_pos) * *pl_data; 00337 pl_data++; 00338 } 00339 *pr_sig = r_sig_f * sqrt(*pr_sig/ *pr_sum); 00340 } 00341 else return(1); 00342 *pr_pos = *pr_pos -1.0; 00343 return(0); 00344 }
