Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members  

/MbsAPI/f_find_peaks.c

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

Generated on Tue Nov 8 10:55:56 2005 for Go4-v2.10-5 by doxygen1.2.15