Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members
Go4Fit

f_find_peaks.c

Go to the documentation of this file.
00001 //-------------------------------------------------------------
00002 //        Go4 Release Package v3.04-01 (build 30401)
00003 //                      28-November-2008
00004 //---------------------------------------------------------------
00005 //   The GSI Online Offline Object Oriented (Go4) Project
00006 //   Experiment Data Processing at EE 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 #include <stdlib.h>
00017 
00018 #define TYPE__SHORT  0
00019 #define TYPE__INT    1
00020 #define TYPE__FLOAT  2
00021 #define TYPE__DOUBLE 3
00022 
00023 #define UP 1
00024 #define DOWN -1
00025 #define HORI 0
00026 #define PRINT 0
00027 
00028 int go4fit_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum);
00029 
00030 void go4fit_find_peaks(
00031         void   *pfData,     // pointer to data
00032         int     lType,      // data type: 0=w 1=l 2=f 3=d
00033         int     lFirstChan, // first channel
00034         int     lPoints,    // number of points
00035         int     lAver,      // channels to sum up
00036         double  dDeltaFact, // noise factor
00037         double  dDeltaMin,  // noise minimum
00038         double *dNoise,     // noise array
00039         int     lPeaksMax,  // Maximum number of peaks
00040         int    *plPeaks,    // returns number of peaks found
00041         int    *plMinima,   // pointer to list of minima
00042         double *pdMinima,   // pointer to list of minima values
00043         int    *plMean,   // pointer to list of maxima
00044         int    *plWidth,    // pointer to list of width
00045         double *pdIntegral, // pointer to list of integrals
00046         int    *plLeft,     // pointer to array of left RMS index
00047         int    *plMaximum,     // pointer to array of maximum index
00048         int    *plRight,     // pointer to array of right RMS index
00049         float  *pfLowNoise,   // pointer to data array of low noise
00050         float  *pfHighNoise)  // pointer to data array of high noise
00051 {
00052   short *pw;
00053   int lMaxIndex,lMinIndex,lIndexMax,lDirection,lIndex,lLastUp,lLastDown,ii,i,lMax,lLeft,lRight,
00054     *plMinimaL,*plMeanL,*plMinimaR,*plMeanR,lPrint,*pl;
00055   float *pf;
00056   double dValue,dLow,dHigh,dDelta,dAver,*pdData,*pdNoise,dMax,dMin,dPos,dSig,dSum,*pd;
00057 
00058   *plPeaks=0;
00059   pdNoise=NULL;
00060   pdData = (double *)malloc(lPoints*8);
00061   pdNoise = (double *)malloc(lPoints*8);
00062   memset(pdData,0,lPoints*8);
00063   memset(pdNoise,0,lPoints*8);
00064   lPrint=PRINT;
00065   dMax=0.;
00066   dMin=1e20;
00067   pw = ((short  *)pfData) + lFirstChan;
00068   pl = ((int    *)pfData) + lFirstChan;
00069   pf = ((float  *)pfData) + lFirstChan;
00070   pd = ((double *)pfData) + lFirstChan;
00071   if(lAver == 0)lAver=1;
00072   lPoints=lPoints/lAver;
00073   dAver = (double) lAver;
00074   // copy the data field into local double field
00075   switch(lType)
00076     {
00077     case TYPE__SHORT:  for(ii=0;ii<lPoints;ii++)
00078       {  *(pdData+ii) = (double)*(pw+lAver*ii);
00079       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pw+lAver*ii+i);
00080       *(pdData+ii) = *(pdData+ii)/dAver;
00081       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00082       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00083     case TYPE__INT:    for(ii=0;ii<lPoints;ii++)
00084       {  *(pdData+ii) = (double)*(pl+lAver*ii);
00085       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pl+lAver*ii+i);
00086       *(pdData+ii) = *(pdData+ii)/dAver;
00087       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00088       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00089     case TYPE__FLOAT:  for(ii=0;ii<lPoints;ii++)
00090       {  *(pdData+ii) = (double)*(pf+lAver*ii);
00091       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pf+lAver*ii+i);
00092       *(pdData+ii) = *(pdData+ii)/dAver;
00093       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00094       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00095     case TYPE__DOUBLE: for(ii=0;ii<lPoints;ii++)
00096       {  *(pdData+ii) = (double)*(pd+lAver*ii);
00097       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pd+lAver*ii+i);
00098       *(pdData+ii) = *(pdData+ii)/dAver;
00099       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00100       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00101     default: printf("Data type %d not supported!\n",lType);
00102       printf("Type 0: short, 1: int, 2: float, 3: double\n");
00103       free(pdData);
00104       free(pdNoise);
00105       return;
00106     }
00107   if(dNoise == NULL)
00108     {// calculate from square root
00109       for(ii=0;ii<lPoints;ii++)
00110    {
00111      if(*(pdData+ii) == 0.0) *(pdNoise+ii) = dDeltaFact;
00112      else                    *(pdNoise+ii) = sqrt(*(pdData+ii)) * dDeltaFact;
00113      if(dDeltaMin > *(pdNoise+ii)) *(pdNoise+ii)=dDeltaMin;
00114    }
00115     }else {// sum up bins from noise input array
00116    for(ii=0;ii<lPoints;ii++)
00117      {  *(pdNoise+ii) = *(dNoise+lFirstChan+lAver*ii);
00118         for(i=1;i<lAver;i++)*(pdNoise+ii) += *(dNoise+lFirstChan+lAver*ii+i);
00119         *(pdNoise+ii) = *(pdNoise+ii)/dAver;
00120      }
00121       }
00122   lMaxIndex=0;
00123   lMinIndex=0;
00124   lIndexMax=lPoints-1;
00125   plMinimaL = (int *)malloc(lPeaksMax*16);
00126   plMeanL = (int *)(plMinimaL+lPeaksMax);
00127   plMinimaR = (int *)(plMinimaL+2*lPeaksMax);
00128   plMeanR = (int *)(plMinimaL+3*lPeaksMax);
00129   memset(plMinimaL,0,lPeaksMax*16);
00130 
00131   dDelta= *pdNoise;
00132   dHigh = *pdData + dDelta;
00133   dLow  = *pdData - dDelta;
00134   if(pfLowNoise)*pfLowNoise=dLow;
00135   if(pfHighNoise)*pfHighNoise=dHigh;
00136   //-----------------------------------------------------
00137   // loop: get points of minima and maxima
00138   lDirection=HORI;
00139   lIndex=0;
00140   lLastDown=0;
00141   lLastUp=0;
00142   while ((lIndex <= lIndexMax) & (lMaxIndex < lPeaksMax-2))
00143     {
00144       dValue = *(pdData+lIndex);
00145       if(dValue > dHigh)
00146    {// direction goes upstairs
00147      if (lDirection != UP)
00148        { // direction was downstairs: minimum found
00149          if(lPrint)printf("%5d - %5d - %5d d %8.1f\n",
00150                 lMinIndex,
00151                 lLastDown,
00152                 lIndex,
00153                 dDelta);
00154          *(plMinimaR+lMinIndex)=lIndex;
00155          *(plMinimaL+lMinIndex)=lLastDown;
00156          lMinIndex++;
00157        }
00158      lLastUp=lIndex;
00159      lDirection=UP;// direction goes upstairs
00160    }
00161       if(dValue < dLow)
00162    {// direction goes downstairs
00163      if (lDirection == UP)
00164        { // direction was upstairs: maximum found
00165          if(lPrint)printf("%5d + %5d - %5d d %8.1f\n",
00166                 lMaxIndex,
00167                 lLastUp,
00168                 lIndex,
00169                 dDelta);
00170          *(plMeanR+lMaxIndex)=lIndex;
00171          *(plMeanL+lMaxIndex)=lLastUp;
00172          lMaxIndex++;
00173        }
00174      lLastDown=lIndex;
00175      lDirection=DOWN;// direction goes downstairs
00176    }
00177       // calculate new limits if noise band has been left
00178       if((dValue < dLow)|(dValue > dHigh))
00179    {
00180           dDelta= *(pdNoise+lIndex);
00181      dHigh = dValue + dDelta;
00182      dLow  = dValue - dDelta;
00183    }
00184       lIndex++;
00185       if(pfLowNoise)*(pfLowNoise+lIndex)=dLow;
00186       if(pfHighNoise)*(pfHighNoise+lIndex)=dHigh;
00187     }
00188   if(lDirection == DOWN)
00189     {
00190       *(plMinimaL+lMinIndex)=lLastDown;
00191       *(plMinimaR+lMinIndex)=lIndexMax;
00192       if(lPrint)printf("%5d - %5d - %5d\n",
00193              lMinIndex,
00194              lLastDown,
00195              lIndexMax);
00196       lMinIndex++;
00197     }
00198   //-----------------------------------------------------
00199   if(lPrint)printf("-------------------------\n");
00200   if(lPrint)printf("Minima %d Maxima %d  \n",lMinIndex,lMaxIndex);
00201   // for each peak there must be a minimum left and right
00202   if(lMinIndex != (lMaxIndex+1))
00203     {
00204       free(plMinimaL);
00205       free(pdData);
00206       free(pdNoise);
00207       *plPeaks=0;
00208       if(lPrint)printf("Error, wrong minima\n");
00209       return;
00210     }
00211   // calculate mean values of minima L&R, mean content, peaks
00212   // start with first minimum
00213   if(lPrint)printf("  [ Left   Right ]  \n");
00214   *plMinima = (*plMinimaL + *plMinimaR)/2;
00215   *pdMinima=0.; // average content at mean minimum channel
00216   for(ii = *plMinimaL;ii <= *plMinimaR;ii++) *pdMinima += *(pdData+ii);
00217   *pdMinima = *pdMinima/(double)(*plMinimaR - *plMinimaL + 1);
00218   if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
00219          *plMinimaL,
00220          *plMinimaR,
00221          *plMinima,
00222          *pdMinima);
00223   for(lIndex=0;lIndex < lMaxIndex;lIndex++)
00224     {
00225       // simple mean peak value, overwritten later
00226       *(plMean+lIndex)=(*(plMeanL+lIndex)+*(plMeanR+lIndex))/2;
00227       if(lPrint)printf("+ [%6d %6d] %6d\n",
00228              *(plMeanL+lIndex),
00229              *(plMeanR+lIndex),
00230              *(plMean+lIndex));
00231       // mean minima values (channel and content)
00232       *(plMinima+lIndex+1) = (*(plMinimaL+lIndex+1) + *(plMinimaR+lIndex+1))/2;
00233       *(pdMinima+lIndex+1)=0.;
00234       for(ii = *(plMinimaL+lIndex+1);ii <= *(plMinimaR+lIndex+1);ii++)
00235    *(pdMinima+lIndex+1) += *(pdData+ii);
00236       *(pdMinima+lIndex+1) = *(pdMinima+lIndex+1)/(double)(*(plMinimaR+lIndex+1) - *(plMinimaL+lIndex+1)+1);
00237       if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
00238              *(plMinimaL+lIndex+1),
00239              *(plMinimaR+lIndex+1),
00240              *(plMinima+lIndex+1),
00241              *(pdMinima+lIndex+1));
00242       // Peak positions, width and sum by momenta
00243       ii = *(plMinimaL+lIndex+1) - *(plMinimaR+lIndex) + 1;
00244       if(ii > 1)
00245    {
00246      go4fit_position(ii,pdData + *(plMinimaR+lIndex),&dPos,&dSig,&dSum);
00247 
00248      *(plMean+*plPeaks)=(int)((dPos+(double)(*(plMinimaR+lIndex)))*dAver)+lFirstChan;
00249      *(plWidth+*plPeaks) =(int)(dSig*dAver + 0.5);
00250      if(*(pdMinima+lIndex) < *(pdMinima+lIndex+1))
00251        dSum=dSum-((double)ii * *(pdMinima+lIndex)+ (double)ii * (*(pdMinima+lIndex+1)-*(pdMinima+lIndex))/2.);
00252           else
00253        dSum=dSum-((double)ii * *(pdMinima+lIndex)- (double)ii * (*(pdMinima+lIndex)-*(pdMinima+lIndex+1))/2.);
00254      *(pdIntegral+*plPeaks) = dSum*dAver;
00255           (*plPeaks)++;
00256    }
00257       else
00258    {
00259      free(plMinimaL);
00260      free(pdData);
00261       free(pdNoise);
00262       *plPeaks=0;
00263      if(lPrint)printf("Error, negative interval\n");
00264      return;
00265    }
00266     }
00267   for(i=0;i<*plPeaks;i++)
00268     {
00269       dMax=0; // search local maximum in interval
00270       dMin=1000000000; // local minimum
00271       pd=pdData+*(plMinima+i);
00272       for(ii=*(plMinima+i);ii<*(plMinima+i+1);ii++){if(dMax < *pd)dMax = *pd; if(dMin > *pd)dMin = *pd; pd++;}
00273       dMin=dMin+(dMax-dMin)/2.; // threashold for RMS
00274       pd=pdData+*(plMinima+i); // search channel of local maximum
00275       for(lMax=*(plMinima+i);lMax<*(plMinima+i+1);lMax++){if(dMax == *pd)break; pd++;} // l_max is now channel of local maximum
00276       lLeft=*(plMinima+i);
00277       lRight=*(plMinima+i+1);
00278       *(plMaximum+i)=lMax;
00279       pd=pdData+lMax;
00280       for(ii=lMax;ii<lRight;ii++){if(dMin > *pd){lRight=ii;break;}pd++;}
00281       pd=pdData+lMax;
00282       for(ii=lMax;ii>lLeft;ii--){if(dMin > *pd){lLeft=ii;break;}pd--;}
00283       *(plLeft+i)=lLeft;
00284       *(plRight+i)=lRight;
00285     }
00286   *(plMinima) = (int)((double)(*(plMinima))*dAver + dAver/2.)+lFirstChan;
00287   for(i=0;i<*plPeaks;i++)
00288     {
00289       // calculate channel numbers of original data
00290       *(plMaximum+i)  = (int)((double)(*(plMaximum+i))*dAver + dAver/2.)+lFirstChan;
00291       *(plLeft+i)     = (int)((double)(*(plLeft+i))*dAver + dAver/2.)+lFirstChan;
00292       *(plRight+i)    = (int)((double)(*(plRight+i))*dAver + dAver/2.)+lFirstChan;
00293       *(plMinima+i+1) = (int)((double)(*(plMinima+i+1))*dAver + dAver/2.)+lFirstChan;
00294     }
00295   if(lPrint)printf("peaks found %d\n",*plPeaks);
00296 
00297   free(plMinimaL);
00298   free(pdData);
00299   free(pdNoise);
00300 }
00301 //**********************************************************************************
00302 int go4fit_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum)
00303 {
00304 
00305 #define LOW  2
00306 #define WIDTH  8
00307 
00308   double r_sig_f     ;
00309 //SL  int I,J,K,L     ;
00310   int J;
00311   double d_sum_prod  ;
00312 //SL  double l_max_chan  ;
00313   double *pl_data,d_max    ;
00314 
00315   /* 2.0E2*SQRT(2.0E0*LOG(2.0E0))/100. */;
00316   r_sig_f      = 55.4518;
00317   r_sig_f      = 2.;
00318   *pr_sum   = 0;
00319   d_sum_prod   = 0;
00320   *pr_sig      = 0.;
00321   *pr_pos      = 0.;
00322   d_max=0.0;
00323 
00324   /* get maximum channel content and integral */
00325   /*  pl_data = pa_data;
00326       for(J = 0; J < l_len; J++)
00327       {
00328       if(d_max < *pl_data)
00329       {
00330       d_max = *pl_data;
00331       l_max_chan = J;
00332       }
00333       pl_data++;
00334       }
00335   */
00336   /* Calculate first momentum */
00337   pl_data = pa_data;
00338   for(J = 1; J <= l_len; J++)
00339     {
00340       *pr_sum = *pr_sum +     *pl_data;
00341       d_sum_prod = d_sum_prod + J * *pl_data;
00342       pl_data++;
00343     }
00344   /* Calculate second momentum */
00345   if(*pr_sum > 0)
00346     {
00347       *pr_pos = (d_sum_prod/ *pr_sum + 0.5);
00348       pl_data = pa_data;
00349       for(J = 1; J <= l_len; J++)
00350    {
00351      *pr_sig = *pr_sig + ((double)J - *pr_pos) * ((double)J - *pr_pos) * *pl_data;
00352      pl_data++;
00353    }
00354       *pr_sig = r_sig_f * sqrt(*pr_sig/ *pr_sum);
00355     }
00356   else return(1);
00357   *pr_pos = *pr_pos -1.0;
00358   return(0);
00359 }
00362 //----------------------------END OF GO4 SOURCE FILE ---------------------
Generated on Fri Nov 28 12:59:11 2008 for Go4-v3.04-1 by  doxygen 1.4.2