StRoot: StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.cxx Source File

00001 #include <stdio.h>
00002 #include <string.h>
00003 #include <math.h>
00004 /*********************************************************************
00005  * $Id: L2Histo.cxx,v 1.6 2010/04/18 06:05:32 pibero Exp $
00006  * \author Jan Balewski, IUCF, 2006 
00007  *********************************************************************
00008  * Descripion:
00009  *  primitive 1-D or 2-D histograming I/O
00010  *  all indexes run [0,nb-1]
00011  *********************************************************************
00012  */
00013 
00014 #include "L2Histo.h"
00015 
00016 //=====================================
00017 //=====================================
00018 void
00019 L2Histo::set( int idx, const char *tit, int nbx, int nby) {
00020   if (!(idx>0 && nbx >0 && nby >0 && tit))
00021     {
00022       head.hid=1;
00023       head.nBin=1;
00024       head.nBinX=1;
00025       head.nBinY=1;
00026       head.dataSize=head.nBin*sizeof(int);
00027       data=new int [head.nBin];
00028       memset(data,0,head.dataSize);
00029       head.title[0]=0;
00030       strncpy(head.title,"set(): HISTOGRAM IS BROKEN IN. CONTACT EXPERT.",mxTx);
00031       return;
00032     }
00033 
00034   head.nOver=head.nUnder=0;
00035   head.ver=version;
00036   head.hid=idx;
00037   head.nBin=nbx*nby;
00038   head.nBinX=nbx;
00039   head.nBinY=nby;
00040   head.dataSize=head.nBin*sizeof(int);
00041   data=new int [head.nBin];
00042   memset(data,0,head.dataSize);
00043   head.title[0]=0;
00044   strncpy(head.title,tit,mxTx);
00045 }
00046 
00047 //=====================================
00048 //=====================================
00049 void
00050 L2Histo::setTitle(const char *tit){
00051   head.title[0]=0;
00052   if (!tit)
00053     {
00054       strncpy(head.title,"setTitle():  HISTOGRAM IS BROKEN. CONTACT EXPERT.",mxTx);
00055       return;
00056     }
00057   strncpy(head.title,tit,mxTx);
00058 }
00059 
00060 //=====================================
00061 //=====================================
00062 L2Histo::L2Histo(){
00063   head.nOver=head.nUnder=-1;
00064   head.ver=version;
00065   head.hid=0;
00066   head.nBin=0;
00067   head.nBinX=head.nBin;
00068   head.nBinY=1;
00069   head.dataSize=0;
00070   data=0;
00071   head.title[0]=0;
00072 }
00073 
00074 //=====================================
00075 //=====================================
00076 void
00077 L2Histo::reset() {
00078   head.nOver=head.nUnder=0;
00079   if(head.dataSize>0) memset(data,0,head.dataSize);
00080 }
00081 
00082 //=====================================
00083 //=====================================
00084 void
00085 L2Histo::print( int flag, FILE *fd){
00086   int dim=1;
00087   if(head.nBinY>1) dim=2;
00088   fprintf(fd," L2Histo::print h%d, %dD, tit='%s' head.nBin=%d nBinX=%d ver=%d\n", head.hid,  dim,  head.title,  head.nBin,  head.nBinX,head.ver);
00089 
00090   int i=0;
00091   int sum=0;
00092   for(i=0;i< head.nBin;i++) {
00093     sum+= data[i];
00094     if(flag) {
00095       if(dim==1) {
00096         fprintf(fd," bin=%2d value=%4d sum=%d\n",i, data[i],sum);
00097       } else {
00098         int ix=i%head.nBinX;
00099         int iy=i/head.nBinX;
00100         fprintf(fd," bin=%d  ix=%d iy=%d value=%d \n",i,ix,iy, data[i]);
00101       }
00102     }
00103   }
00104   fprintf(fd,"   Counts: inRange=%d  nUnder=%d, nOver=%d\n",sum,  head.nUnder,  head.nOver);
00105 }
00106 
00107 
00108 //=====================================
00109 //=====================================
00110 void
00111 L2Histo::printCSV(  FILE *fd){
00112   if (!fd) {
00113     printf("printCSV() passed bad output file!\n");
00114     return;
00115   }
00116   fprintf(fd,"#L2H%d,%d,%s,", head.hid,  head.nBin,  head.title);
00117   if(head.nBinY>1) {
00118     fprintf(fd,"ny=%d,CVS format not supported for 2D histos,\n",head.nBinY);
00119     return;
00120   }
00121   int i=0;
00122 
00123   for(i=0;i< head.nBin;i++)  fprintf(fd,"%d,",data[i]);
00124   fprintf(fd,"\n");
00125   return;
00126 }
00127 
00128 
00129 //=====================================
00130 //=====================================
00131 void
00132 L2Histo::fill( int bin){
00133   /* should be:   assert(head.nBinY==1) 
00134      but it is an online code - proceed, good luck
00135   */ 
00136   if(bin<0) { head.nUnder++; return;}
00137   if(bin>=head.nBin) { head.nOver++; return;}
00138   data[bin]++;
00139 }
00140 
00141 
00142 //=====================================
00143 //=====================================
00144 void
00145 L2Histo::fillW( int bin, int w){
00146   //Fill data[bin] with weight w.
00147   /* should be:   assert(head.nBinY==1) 
00148      but it is an online code - proceed, good luck
00149   */ 
00150   if(bin<0) { head.nUnder+=w; return;}
00151   if(bin>=head.nBin) { head.nOver+=w; return;}
00152   data[bin]+=w;
00153 }
00154 
00155 //=====================================
00156 //=====================================
00157 void
00158 L2Histo::fill( int binX, int binY){
00159   /* will work fine for 1-D, unless binY is negative;
00160      just computes too much
00161      inverse is not true !
00162   */ 
00163   if(binX<0 || binY <0) { head.nUnder++; return;}
00164   if(binX>=head.nBinX || binY>=head.nBinY) { head.nOver++; return;}
00165   data[binY*head.nBinX+binX]++;
00166 }
00167 
00168 //=====================================
00169 //=====================================
00170 void
00171 L2Histo::fillW( int binX, int binY, int w){
00172   /* will work fine for 1-D, unless binY is negative;
00173      just computes too much
00174      inverse is not true !
00175   */ 
00176   if(binX<0 || binY <0) { head.nUnder+=w; return;}
00177   if(binX>=head.nBinX || binY>=head.nBinY) { head.nOver+=w; return;}
00178   data[binY*head.nBinX+binX]+=w;
00179 }
00180 
00181 //=====================================
00182 //=====================================
00183 void
00184 L2Histo::write(FILE *fd, int dbg) {
00185   if(fd==0) return;
00186   if( head.nBin<=0) return;
00187   if(dbg) print(dbg-1);
00188   char *c;
00189   unsigned int i;
00190   //.... save histo header
00191   c=(char *)&head;
00192   for(i=0;i<sizeof(L2Histo::Head);i++)  fputc(c[i],fd);
00193  
00194    //.... save histo data
00195   c=(char *) data;
00196   for(i=0;i< head.dataSize;i++)  fputc(c[i],fd);
00197 
00198 } 
00199 
00200 
00201 //=====================================
00202 //=====================================
00203 int
00204 L2Histo::read(FILE *fd, int dbg) {
00205   if (!fd)
00206     {
00207       printf("L2Histo::read called with no file.  Aborting\n");
00208       return -1;
00209     }
00210   char *c;
00211   unsigned int i;
00212   
00213   // printf("aa=%d\n",sizeof(*h));
00214   //.... read histo header
00215   c=(char *)&head;
00216   for(i=0;i<sizeof(L2Histo::Head);i++) {
00217     int val=fgetc(fd);
00218     if(i==0 && val==EOF) return 1;
00219     if(!(val!=EOF))
00220       {
00221         printf("L2Histo::read  val==EOF.  Aborting\n");
00222         return -1;
00223       }
00224     c[i]=val;
00225   } 
00226   //  printf("rd ver %d %d \n", head.ver,version);
00227   if (!(  head.ver==version)) // change it in the future if header ever changes
00228     {
00229       printf("L2Histo::read head.ver!=version.  Aborting\n");
00230       return -1;
00231     }
00232   if (!( head.nBin>0))
00233     {
00234       printf("L2Histo::read  head.nBin<=0.  Aborting\n");
00235       return -1;
00236     }
00237   if (!( head.dataSize== head.nBin*sizeof(int)))
00238     {
00239       printf("L2Histo::read  head.dataSize!= head.nBin*sizeof(int).  Aborting\n");
00240       return -1;
00241     }
00242 
00243   data=new int [head.dataSize];
00244   memset( data,0, head.dataSize); 
00245   
00246   //.... read histo data
00247   c=(char *) data;
00248   for(i=0;i< head.dataSize;i++) {
00249     int val=fgetc(fd);
00250     if(dbg)printf("i=%d val=%d\n",i,val);
00251     if(i==0 && val==EOF) return 1;
00252     if(!(val!=EOF))
00253     {
00254       printf("L2Histo::read  val==EOF.  Aborting\n");
00255       return -1;
00256     }
00257     c[i]=val;
00258   }
00259 
00260  if(dbg) printf("got histo hid=%d  , head.nBin=%d, tit='%s'\n", head.hid,  head.nBin, head.title);
00261 
00262   if(dbg) print(1);
00263   return 0;
00264 }
00265 
00266 
00267 /*========================================
00268   ======================================== */
00269 bool 
00270 L2Histo::findMax( int *iMax, int *iFWHM) {
00271   // returns true if both max & FWHM make sense
00272   *iMax=*iFWHM=-1;  
00273 
00274   /* finds pedestal & FWHM */
00275   int maxVal=-1, maxJ=-1;
00276   int j;
00277   for(j=0; j<head.nBin;j++) {
00278     // printf("j=%d  maxVal=%d maxJ=%d  yield=%d\n",j,maxVal, maxJ,data[j]);
00279     if(data[j]<=maxVal) continue;
00280     maxVal=data[j];
00281     maxJ=j;
00282   }
00283 
00284   if(maxJ<0) return false; /* no maximum found */
00285 
00286   *iMax=maxJ; // maximum was found
00287 
00288   if(head.nBinY>1) { // FWHM makes no sense for 2D plots
00289    *iMax=maxJ;
00290    *iFWHM=999;
00291    return true;
00292   }
00293 
00294 
00295   // search for FWHM
00296   int halfMax=maxVal/2;
00297    int j1=-1, j2=-1;
00298    for(j=maxJ+1; j<head.nBin;j++) {// forward
00299      if(data[j]>halfMax) continue;
00300      j2=j;
00301      break;
00302    }
00303 
00304    for(j=maxJ-1; j>=0;j--) { // backward
00305      if(data[j]>halfMax) continue;
00306      j1=j;
00307      break;
00308    }
00309    if(j1<0 || j2<0) return false; /* FWHM is 0 ? */
00310 
00311 
00312    *iFWHM=(j2-j1);
00313  
00314    return true;
00315 }
00316 
00317 /*========================================
00318   ======================================== */
00319 bool 
00320 L2Histo::findMean( int *iMean, int *iRMS) {
00321   // returns true if both mean & RMS make sense
00322 
00323   *iMean=*iRMS=-1;  
00324   if (head.nBinY>1)
00325     return false; //this method doesn't work for 2d plots
00326   /* finds mean & RMS without counting overflows */
00327   int totalVal=0, weightedVal=0, meanJ=-1, rms=-1;
00328   int j;
00329   for(j=0; j<head.nBin;j++) {
00330     // printf("j=%d  maxVal=%d maxJ=%d  yield=%d\n",j,maxVal, maxJ,data[j]);
00331     weightedVal+=data[j]*j;
00332     totalVal+=data[j];
00333   }
00334   if (totalVal>0) meanJ=weightedVal/totalVal;
00335   if(meanJ<0) return false; /* no maximum found */
00336 
00337   *iMean=meanJ; // mean was found
00338 
00339   // search for RMS
00340   weightedVal=0;
00341   for(j=0; j<head.nBin;j++) {
00342     // printf("j=%d  maxVal=%d maxJ=%d  yield=%d\n",j,maxVal, maxJ,data[j]);
00343     weightedVal+=data[j]*(j-meanJ)*(j-meanJ);
00344   }
00345   rms=weightedVal/totalVal;
00346 
00347   *iRMS=rms; 
00348    return true;
00349 }
00350 
00351 //=====================================
00352 //=====================================
00353 void
00354 L2Histo::printPed( FILE *fd, int x0, int x1,char term){
00355   if(head.nBinY>1) {
00356     fprintf(fd," printped NOT implemented for 2D histos: hid=%d tit=%s\n",head.hid, head.title);
00357     return;
00358   }
00359      int j;
00360       for(j=0; j<head.nBin;j++) {
00361         int adc=x0+j;
00362         if(adc>=x1) break;
00363         float val=data[j];
00364         if(val<=0) { // print grid for empty bins 
00365           if(adc==0) 
00366             fprintf(fd,"*");
00367           else if(adc==10 ) 
00368             fprintf(fd,">");
00369           else if(adc==-10 ) 
00370             fprintf(fd,"<");
00371           else if(adc%10==0) 
00372             fprintf(fd,",");
00373           else
00374             fprintf(fd,".");
00375           continue;
00376         }
00377         // print yield for populated bins
00378 
00379         fprintf(fd,"%c",y2c(val));
00380       }
00381 
00382       // sum higher bins t print overflow
00383       int sum=head.nOver;
00384       for(j=x1; j<head.nBin;j++) sum+=data[j];
00385       
00386       fprintf(fd,"ov=%d %c",sum,term);
00387 }
00388 
00389 //=====================================
00390 //=====================================
00391 char
00392 L2Histo::y2c(float val){
00393   float valLog=log(val);
00394   char k='?';
00395   if(valLog<0.) 
00396     k='-';
00397   else if(valLog<9) 
00398     k='1'+(int)valLog;
00399   else  
00400     k='a'+(int)(valLog-9);
00401   return k;
00402 }
00403 
00404 
00405 
00406 /*
00407 *********************************************************************
00408   $Log: L2Histo.cxx,v $
00409   Revision 1.6  2010/04/18 06:05:32  pibero
00410   Address compiler warnings.
00411 
00412   Revision 1.5  2010/04/17 16:42:09  pibero
00413   *** empty log message ***
00414 
00415   Revision 1.4  2008/01/16 23:32:33  balewski
00416   toward token dependent compute()
00417 
00418   Revision 1.3  2008/01/10 22:45:39  balewski
00419   reduce memory allocation for every histo
00420 
00421   Revision 1.2  2007/11/14 03:58:07  balewski
00422   cleanup of common timing measurement
00423 
00424   Revision 1.1  2007/10/11 00:33:14  balewski
00425   L2algo added
00426 
00427   Revision 1.3  2006/03/28 19:33:23  balewski
00428   ver16b , in L2new
00429 
00430   Revision 1.2  2006/03/11 17:08:32  balewski
00431   now CVS comments should work
00432 
00433 */
00434