~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/util/hist_funcs.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 #include <stdio.h>
  2 #include <math.h>
  3 #include <values.h>
  4 #include <tina/hist_funcs.h>
  5 
  6 double (*mfunc)(int order,  double *params, float x) = NULL;
  7 double (*merror_func)(shistogram *ph, float x)=NULL;
  8 shistogram *mph;
  9 
 10 static double **maketarray(int l_array, int w_array, float gl)
 11 {
 12     double **array;
 13     short   i, j;
 14 
 15     if ((array = (double **) ralloc((unsigned) l_array * sizeof(double *)))
 16                                == NULL)
 17         return (NULL);
 18     if ((array[0] = (double *) ralloc((unsigned) l_array * w_array * sizeof(double)))
 19                                == NULL)
 20         return (void *) (NULL);
 21     for (i = 1; i < l_array; ++i)
 22         array[i] = array[i - 1] + w_array;
 23     for (i = 0; i < l_array; ++i)
 24         for (j = 0; j < w_array; j++)
 25             array[i][j] = gl;
 26     return (array);
 27 }
 28 
 29 static void    freetarray(double **array)
 30 {
 31     if (array == NULL) return;
 32     if (array[0] != NULL) rfree((void *) array[0]);
 33     if (array != NULL) rfree((void *) array);
 34 }
 35 
 36 static float **makeharray(int l_array,int w_array,float gl)
 37 {
 38      float **array;
 39      short i,j;
 40      if((array=(float **)ralloc(l_array*sizeof(float *)))==NULL)
 41            herror("no room for array");
 42      if((array[0]=(float *)ralloc(l_array*w_array*sizeof(float)))==NULL)
 43            herror("no room for array");
 44      for (i=1;i<l_array;++i)
 45           array[i]=array[i-1]+w_array;
 46      for (i=0;i<l_array;++i)
 47            for (j=0;j<w_array;j++)
 48                array[i][j] = gl;
 49      return(array);
 50 }         
 51 
 52 static void freeharray(float **array)
 53 {
 54     int i;
 55     if (array == NULL) return;
 56     if (array[0]!=NULL) rfree(array[0]);        
 57     if (array!=NULL) rfree(array);      
 58 }
 59 
 60 shistogram *hbook1(int id, char *title,float xmin,float xmax,int xbins)
 61 {
 62      shistogram *ph;
 63      int i;
 64 
 65      if (id >= NHIST) 
 66      {
 67         herror("cannot book histogram \n");
 68         return(NULL);
 69      }
 70      ph = (shistogram *)ralloc(sizeof(shistogram));
 71      ph->id = id;
 72      ph->title = title;
 73      ph->type = 1.0;
 74      ph->super = NULL;
 75      ph->xmax = xmax;
 76      ph->xmin = xmin;
 77      ph->contents = 0.0;
 78      ph->mean = 0.0;
 79      ph->mean2 = 0.0;
 80      ph->under = 0.0;
 81      ph->over = 0.0; 
 82      ph->entries = 0.0;
 83      ph->xbins = xbins;
 84      ph->ybins=1;
 85      ph->xincr = (xmax-xmin)/(float)xbins;
 86      ph->array=makeharray(1,xbins,0.0);
 87      return(ph);
 88 }
 89 
 90 shistogram *hfree(shistogram *ph)
 91 {
 92    if (ph!=NULL)
 93    {
 94       freeharray(ph->array);
 95       if (ph->par != NULL) rfree(ph->par);
 96       rfree(ph);
 97    }
 98    return(NULL);
 99 }
100 
101 shistogram *hbook2(int id,char *title,float xmin,float xmax,int xbins,
102                                       float ymin,float ymax,int ybins)
103 {
104      shistogram *ph;
105      int i;
106 
107      if (id >= NHIST) 
108      {
109         herror("cannot book histogram \n");
110         return(NULL);
111      }
112      ph = (shistogram *)ralloc(sizeof(shistogram));
113      ph->id = id;
114      ph->title = title;
115      ph->type = 2;
116      ph->super= NULL;
117      ph->xmax = xmax;
118      ph->xmin = xmin;
119      ph->ymin = ymin;
120      ph->ymax = ymax;
121      ph->contents = 0.0;
122      ph->entries = 0.0;
123      ph->mean = 0.0;
124      ph->under = 0.0;
125      ph->over =0.0; 
126      ph->above = 0.0;
127      ph->below = 0.0;
128      ph->xbins = xbins;
129      ph->ybins = ybins;
130      ph->xincr = (xmax-xmin)/(float)xbins;
131      ph->yincr = (ymax-ymin)/(float)ybins;
132      ph->array = makeharray(ybins,xbins,0.0);
133      return(ph);
134 }
135 
136 float hfill1(shistogram *ph,float x,float w)
137 {
138    int i;
139    if (ph==NULL) return(0.0);
140    i = floor((x-ph->xmin)/ph->xincr);
141    if (x<ph->xmax && x>=ph->xmin )
142    {
143       if (w==0.0) return(ph->array[0][i]);
144       ph->contents+=w;
145       ph->entries++;
146       ph->mean += w*(x-ph->mean)/ph->contents;
147       ph->mean2 += w*(x*x-ph->mean2)/ph->contents;
148 /* correct numerical instability on first time through */
149       if (ph->mean*ph->mean>ph->mean2)
150           ph->mean2 = ph->mean*ph->mean;
151       return(ph->array[0][i]+=w);
152    }
153    else
154    {
155       if (x>ph->xmax) ph->over += w;
156       else ph->under += w;
157       return(0.0);
158    }
159 }
160 
161 float hfill1s(shistogram *ph,float x,float w)
162 {
163    int i;
164    double fac;
165    if (ph==NULL) return(0.0);
166    i = floor((x-ph->xmin)/ph->xincr + 0.5);
167    if (x<ph->xmax && x>=ph->xmin )
168    {
169       if (w==0.0) return(ph->array[0][i]);
170       ph->contents+=w;
171       ph->entries++;
172       ph->mean += w*(x-ph->mean)/ph->contents;
173       ph->mean2 += w*(x*x-ph->mean2)/ph->contents;
174 /* correct numerical instability on first time through */
175       if (ph->mean*ph->mean>ph->mean2)
176           ph->mean2 = ph->mean*ph->mean;
177       fac = (x-i*ph->xincr-ph->xmin)/ph->xincr;
178       if (i>0) ph->array[0][i-1]+=w*(0.5-fac); 
179       if (i<ph->xbins) ph->array[0][i]+=w*(0.5+fac); 
180       return(ph->array[0][i]);
181    }
182    else
183    {
184       if (x>ph->xmax) ph->over += w;
185       else ph->under += w;
186       return(0.0);
187    }
188 }
189 
190 float hfill2(shistogram *ph,float x,float y,float w)
191 {
192    int i,j;
193    if (ph==NULL||w==0.0) return(0.0);
194    i = floor((x-ph->xmin)/ph->xincr);
195    j = floor((y-ph->ymin)/ph->yincr);
196    if (x<ph->xmax && x>=ph->xmin && y<ph->ymax && y>=ph->ymin)
197    {
198       ph->contents+=w;
199       ph->entries++;
200       ph->mean += w*(x-ph->mean)/ph->contents;
201       return(ph->array[j][i]+=w);
202    }
203    else 
204    {
205       if(x<ph->xmin) return(ph->under +=w);
206       if(x>=ph->xmax) return(ph->over +=w);
207       if(y<ph->ymin) return(ph->below += w);
208       if(y>=ph->ymax) return(ph->above += w);
209    }
210 }
211 
212 void hprint(FILE *fp,shistogram *ph)
213 {
214    int i,j;
215    float scale=1.0,max,variance,y;
216    char temp[20];
217    if(ph==NULL)
218    {
219       fprintf(fp,"histogram empty! \n");
220       return;
221    }
222 
223    if(ph->entries == 0)
224       fprintf(fp," h id %d %s empty \n\n",ph->id,ph->title);
225    else
226       fprintf(fp," h id %d %s\n\n",ph->id,ph->title);
227 
228    if (ph->type==1)
229    {
230       if(ph->entries!=0)
231       {
232          max = 0.0;
233          for (i=0;i<ph->xbins;i++)
234          {
235             if (ph->array[0][i]>max) max=ph->array[0][i];
236          }
237          scale =  max/100.0;
238          scale = pow(10.0,(double)((int)(log10((double)scale))));
239          if(max>50.0*scale) scale *= 2.0;
240          if(max>50.0*scale) scale *= 2.5;
241          if(max>50.0*scale) scale *= 2.0;
242          if(max>50.0*scale) scale *= 2.0;
243 
244          for (i=0;i<ph->xbins;i++)
245          {
246             fprintf(fp," %6.2f ",ph->xmin+i*ph->xincr);
247             if(ph->super==NULL)
248             {
249                for (j=0;(float)j<ph->array[0][i]/scale;j++)
250                {
251                   fprintf(fp,"x");
252                }
253             }
254             else
255             {
256                y = ph->super(ph->npar,ph->par,(ph->xmin+(i+0.5)*ph->xincr));
257                if (y/scale>60.0) y = scale*60.0;
258 
259                for (j=0;(float)j<ph->array[0][i]/scale
260                       ||(float)j<=y/scale;j++)
261                { 
262                    if ((int)(y/scale)==j) fprintf(fp,"*");
263                    else if ((float)j+0.5<ph->array[0][i]/scale) fprintf(fp,"x");
264                    else fprintf(fp," ");
265                }
266             }
267             fprintf(fp,"\n");
268          }
269          fprintf(fp,"\n");
270       }
271       fprintf(fp," h id %d entries %d contents %5.3f underflows %5.3f overflows %5.3f\n",
272               ph->id,ph->entries,ph->contents,ph->under,ph->over);
273       variance = (float)sqrt(fabs(ph->mean2-ph->mean*ph->mean));
274       fprintf(fp," h id %d mean  %5.3f var(**0.5) %5.3f scale %5.3f res %5.3f run %5.3f grad %5.3f \n",
275               ph->id,ph->mean,variance,scale,hlnorm(ph,2),hrunstat(ph),hgradstat(ph));
276       fprintf(fp,"\n");
277       fprintf(fp,"\n");
278    }
279    else if(ph->type==2)
280    {
281       if(ph->entries!=0)
282       {
283          for (j=ph->ybins-1;j>=0;j--)
284             {
285             fprintf(fp," %6.2f ",ph->ymin+j*ph->yincr);
286             for (i=0;i<ph->xbins;i++)
287             {
288                 if(ph->array[j][i]>=1.0 && ph->array[j][i]<10.0)
289                     fprintf(fp,"%d",(int)ph->array[j][i]);
290                 else if(ph->array[j][i]==0.0)
291                     fprintf(fp," ");
292                 else if(ph->array[j][i]<1.0 && ph->array[j][i]>0.0)
293                     fprintf(fp,".");
294                 else if(ph->array[j][i]>=200.0)
295                     fprintf(fp,"*");
296                 else if(ph->array[j][i]>=150.0)
297                     fprintf(fp,"H");
298                 else if(ph->array[j][i]>=100.0)
299                     fprintf(fp,"G");
300                 else if(ph->array[j][i]>=65.0)
301                     fprintf(fp,"F");
302                 else if(ph->array[j][i]>=45.0)
303                     fprintf(fp,"E");
304                 else if(ph->array[j][i]>=30.0)
305                     fprintf(fp,"D");
306                 else if(ph->array[j][i]>=20.0)
307                     fprintf(fp,"C");
308                 else if(ph->array[j][i]>=15.0)
309                     fprintf(fp,"B");
310                 else if(ph->array[j][i]>=10.0)
311                     fprintf(fp,"A");
312                 else if(ph->array[j][i]<0.0)
313                     fprintf(fp,"-");
314                 else 
315                     fprintf(fp,"?");
316             }
317             fprintf(fp,"\n");
318          }
319          sprintf(temp,"%4.2f",ph->xmin);
320          for (j=0;temp[j]!='\0';j++)
321          {
322             fprintf(fp,"        ");
323             for (i=0;i<ph->xbins;i++)
324             {
325                sprintf(temp,"%4.2f",ph->xmin+i*ph->xincr);
326                fprintf(fp,"%c",temp[j]);
327             }
328             fprintf(fp,"\n");
329          }
330          fprintf(fp,"\n");
331          fprintf(fp," Key: .<1.0,  A<15.0,  B<20.0,  C<30.0,  D<45.0,\n");
332          fprintf(fp,"      E<65.0  F<100.0, G<150.0, H<200.0, *<Infty.\n");
333       }
334 
335       fprintf(fp," entries %d contents %6.2f\n", ph->entries,ph->contents);
336       fprintf(fp," under%6.2f over %6.2f above %6.2f below %6.2f\n",
337               ph->under,ph->over,ph->above,ph->below);
338       fprintf(fp,"\n");
339    }
340 }
341 
342 double hresidual(shistogram *ph, int i)
343 {
344     double y, residual;
345     y = ph->super(ph->npar,ph->par,(ph->xmin+(i+0.5)*ph->xincr));
346     residual = y-ph->array[0][i];
347     return(residual);
348 }
349 
350 double hlnorm(shistogram *ph, int norm)
351 {
352     int i;
353     double residual = 0.0;
354     if (ph->super==NULL) return(0.0);
355     for (i=1;i<ph->xbins-1;i++)
356     {
357         if (norm == 1) residual += fabs(hresidual(ph,i));
358         else if (norm==2) residual += hresidual(ph,i)*hresidual(ph,i);
359     }
360     residual /= (float)(ph->xbins-2);
361  
362     if (norm == 1) return(residual);
363     else return(sqrt(residual));
364 }
365 
366 double hgradstat(shistogram *ph)
367 {
368     int i;
369     double y, old_y, data, old_data;
370     double grad, gradtest;
371     int count=0;
372     if (ph->super==NULL) return(0.0);
373     old_y = ph->super(ph->npar,ph->par,(ph->xmin+(0.5)*ph->xincr));
374     old_data = ph->array[0][0];
375     
376     for (i=1;i<ph->xbins;i++)
377     {
378         y = ph->super(ph->npar,ph->par,(ph->xmin+(i+0.5)*ph->xincr));
379         data = ph->array[0][i];
380         gradtest = (y-old_y )*(data - old_data);
381         if ( gradtest>0.0 )
382         {
383            count++;
384         }
385         old_y = y;
386         old_data = data;
387     }
388 
389     return((double)count/((float)ph->xbins -1));
390 }
391 
392 double hrunstat(shistogram *ph)
393 {
394     int i,count[8]={0,0,0,0,0,0,0,0};
395     int ncount = 0;
396     double new_residual,residual=0.0,y;
397     double bhatt = 0.0, temp=1.0;
398 
399     if (ph->super==NULL) return(0.0);
400 /* analyse run count */
401     for (i=0;i<ph->xbins;i++)
402     {
403         y = ph->super(ph->npar,ph->par,(ph->xmin+(i+0.5)*ph->xincr));
404         new_residual = y-ph->array[0][i];
405         if ((residual <0.0 && new_residual<0.0)
406           ||(residual >=0.0 && new_residual>=0.0))
407         {
408            if (i!=0) ncount++;
409         }
410         else
411         {
412            ncount = 0;
413         }
414         if(ncount>7) ncount = 7;
415         count[ncount]++;
416         residual = new_residual;
417     }
418 /* compare frequency distribution against expected using B statistic */
419     for (i=0;i<7;i++)
420     {
421         temp *= 0.5;
422         bhatt += sqrt(temp)*sqrt((float)count[i]);
423     }
424 /* finish off last bin */
425     bhatt += sqrt(temp)*sqrt((float)count[7]); 
426     return(bhatt/sqrt((float)ph->xbins)); 
427 }
428 
429 void histdo(FILE *fp)
430 {
431     int i;
432     if (fp==NULL) return;
433     for (i=0;i<NHIST;i++)
434     {
435        if(hist[i]!=NULL) hprint(fp,hist[i]);
436     }
437 }
438 
439 void hreset(shistogram *ph)
440 {
441     int i,j;
442 
443     if(ph==NULL)
444     {
445        for (i=0;i<NHIST;i++)
446        {
447           if(hist[i]!=NULL) hreset(hist[i]);
448        }
449     }
450 
451     else
452     {
453        for(j=0;j<ph->ybins;j++)
454        {
455          for(i=0;i<ph->xbins;i++)
456          {
457            ph->array[j][i] = 0.0;
458          }
459        }
460        ph->contents = 0.0;
461        ph->entries = 0.0;
462        ph->under=0.0;
463        ph->over=0.0;
464        ph->above=0.0;
465        ph->below=0.0;
466        ph->mean=0.0;
467        ph->mean2=0.0;
468     }
469 }
470 
471 void hpxprint(FILE *fp, shistogram *ph)
472 {
473     int i,j;
474     char *title;
475     shistogram *temp=NULL;
476 
477     if(ph==NULL)
478     {
479         return;
480     }
481     else
482     {
483        title = (char *)ralloc(256*sizeof(char));
484        sprintf(title," X projection > %s ",ph->title);
485        temp = hbook1(NHIST-1,title,ph->xmin,ph->xmax,ph->xbins);
486        for(j=0;j<ph->ybins;j++)
487        {
488          for(i=0;i<ph->xbins;i++)
489          {
490              hfill1(temp,ph->xmin+(i+0.5)*ph->xincr,ph->array[j][i]); 
491          }
492        }
493     }
494     hprint(fp,temp);
495     freeharray(temp->array);
496     rfree(temp);
497 }
498 
499 void hpyprint(FILE *fp, shistogram *ph)
500 {
501     int i,j;
502     char *title;
503     shistogram *temp=NULL;
504 
505     if(ph==NULL)
506     {
507         return;
508     }
509     else
510     {
511        title = (char *)ralloc(256*sizeof(char));
512        sprintf(title," Y projection > %s ",ph->title);
513        temp = hbook1(NHIST-1,title,ph->ymin,ph->ymax,ph->ybins);
514        for(j=0;j<ph->ybins;j++)
515        {
516          for(i=0;i<ph->xbins;i++)
517          {
518              hfill1(temp,ph->ymin+(j+0.5)*ph->yincr,ph->array[j][i]);
519          }
520        }
521     }
522     hprint(fp,temp);
523     freeharray(temp->array);
524     rfree(temp);
525 }
526 
527 double hdiag(shistogram *ph)
528 {
529     int i,j;
530     double diag = 0.0;
531     if (ph && ph->type == 2)
532     {
533        for(j=0;j<ph->ybins;j++)
534        {
535          if(j<ph->xbins)
536          {
537            diag += ph->array[j][j];
538          }
539        }
540     }
541     return(diag);
542 }
543 
544 void hintegf(shistogram *ph)
545 {
546     int i,j;
547     char *title;
548 
549     title = (char *)ralloc(256*sizeof(char));
550 
551     if (ph->type == 1)
552     {
553        for (i=1;i<ph->xbins;i++)
554            ph->array[0][i] +=  ph->array[0][i-1];
555        sprintf(title," Integrated > %s ",ph->title);
556        ph->title = title;
557     }
558 }
559 
560 void hintegb(shistogram *ph)
561 {
562     int i,j;
563     char *title;
564 
565     title = (char *)ralloc(256*sizeof(char));
566 
567     if (ph->type == 1)
568     {
569        for (i=ph->xbins;i>0;i--)
570            ph->array[0][i-1] +=  ph->array[0][i];
571        sprintf(title," Integrated < %s ",ph->title);
572        ph->title = title;
573     }
574 }
575 
576 void hsmoof(shistogram *ph)
577 {
578     int i,j;
579 
580     if (ph->type == 1)
581     {
582          for (i=1;i<ph->xbins-1;i++)
583              ph->array[0][i-1] += 0.25*ph->array[0][i]-0.125*ph->array[0][i-1];
584              ph->array[0][i+1] += 0.25*ph->array[0][i];
585              ph->array[0][i] = 0.5*ph->array[0][i];
586     }
587     else if (ph->type==2)
588     {
589         for(j=1;j<ph->ybins-1;j++)
590         {
591             for (i=1;i<ph->xbins-1;i++)
592                 ph->array[j][i] = 0.0625*ph->array[j-1][i-1]
593                                 + 0.0625*ph->array[j-1][i]
594                                 + 0.0625*ph->array[j-1][i+1]
595                                 + 0.0625*ph->array[j][i-1]
596                                 + 0.5*ph->array[j][i]
597                                 + 0.0625*ph->array[j][i+1]
598                                 + 0.0625*ph->array[j+1][i-1]
599                                 + 0.0625*ph->array[j+1][i]
600                                 + 0.0625*ph->array[j+1][i+1];
601         }
602     }
603 }
604 
605 float hmax1(shistogram *ph, float *maxx)
606 {
607     int i;
608     float hmax = ph->array[0][0];
609     int imax = 0;
610     double centre = 0.0;
611     double h0,h1,hm1;
612 
613     if (ph->type == 1)
614     {
615          for (i=0;i<ph->xbins;i++)
616          {
617              if (ph->array[0][i]>hmax)
618              {
619                  hmax = ph->array[0][i];
620                  imax = i;
621              }
622          }
623     }
624     else
625        herror("wrong histogram type to hmax1");
626     h0  = hmax;
627     h1  = ph->array[0][imax+1];
628     hm1 = ph->array[0][imax-1];
629     /* assuming parabola fit to 3 centre bins */
630     centre = (hm1-h1)/(2*(h1+hm1)-4*h0);
631     *maxx = ph->xmin+(imax+0.5)*ph->xincr;
632     *maxx = centre*(ph->xincr) + *maxx;
633 
634     return(hmax);
635 }
636 
637 float hnearmax1(shistogram *ph, float *maxx)
638 {
639     int i;
640     float hmax = ph->array[0][0];
641     int imax = 0, tmax;
642     double centre = 0.0;
643     double h0,h1,hm1,newmaxx,olddiff=MAXFLOAT;
644 
645     hmax = ph->array[0][0];
646     imax = 0;
647 
648     if (ph->type == 1)
649     {
650          for (i=1;i<ph->xbins;i++)
651          {
652              if (ph->array[0][i]>ph->array[0][i-1] &&
653                  ph->array[0][i]>ph->array[0][i+1])
654              {
655                 if (ph->array[0][i]>hmax)
656                 {
657                    h0  = ph->array[0][i];
658                    h1  = ph->array[0][i+1];
659                    hm1 = ph->array[0][i-1];
660                    /* assuming parabola fit to 3 centre bins */
661                    centre = (hm1-h1)/(2*(h1+hm1)-4*h0);
662                    newmaxx = ph->xmin+(i+0.5)*ph->xincr;
663                    newmaxx = centre*(ph->xincr) + newmaxx;
664                    if (fabs(newmaxx-*maxx) < olddiff)
665                    {
666                       olddiff = fabs(newmaxx-*maxx);
667                       hmax = ph->array[0][i];
668                       imax = i;
669                    }
670                 }
671              }
672          }
673     }
674     else
675        herror("wrong histogram type to hmax1");
676     h0  = hmax;
677     h1  = ph->array[0][imax+1];
678     hm1 = ph->array[0][imax-1];
679     /* assuming parabola fit to 3 centre bins */
680     centre = (hm1-h1)/(2*(h1+hm1)-4*h0);
681     *maxx = ph->xmin+(imax+0.5)*ph->xincr;
682     *maxx = centre*(ph->xincr) + *maxx;
683 
684     return(hmax);
685 }
686 
687 float hmax2(shistogram *ph, float *maxx, float *maxy)
688 {
689     int i,j;
690     float hmax = ph->array[0][0];
691     int imax = 0, jmax=0;
692 
693     if (ph->type == 1)
694        herror("wrong histogram type to hmax2");
695     else
696     {
697         for (j=0;j<ph->ybins;j++)
698         {
699             for (i=0;i<ph->xbins;i++)
700             {
701                 if (ph->array[j][i]>hmax)
702                 {
703                     hmax = ph->array[j][i];
704                     imax = i;
705                     jmax = j;
706                 }
707             }
708         }
709     }
710     *maxx = ph->xmin+(imax+0.5)*ph->xincr;
711     *maxy = ph->ymin+(jmax+0.5)*ph->yincr;
712 
713     return(hmax);
714 }
715 
716 void hstore(FILE *fp,shistogram *ph)
717 /* warning not yet tested */
718 {
719     int i,j;
720 
721     fprintf(fp,"%d \n",ph->id);
722     fprintf(fp,"%s \n",ph->title);
723     fprintf(fp,"%d \n",ph->type);
724     fprintf(fp,"%3.2f %3.2f %d \n",ph->xmin,ph->xmax,ph->xbins);
725     fprintf(fp,"%3.2f %3.2f %d \n",ph->ymin,ph->ymax,ph->ybins);
726     if (ph->type == 1)
727     {
728          for (i=0;i<ph->xbins;i++)
729              fprintf(fp,"%6.2f %3.2f \n ",ph->xmin+(i+0.5)*ph->xincr,ph->array[0][i]);
730     }
731     else if (ph->type==2)
732     {
733         for(j=0;j<ph->ybins;j++)
734         {
735             for (i=0;i<ph->xbins;i++)
736                 fprintf(fp,"%6.2f %6.2f %3.2f \n ",
737                     ph->xmin+(i+0.5)*ph->xincr,ph->ymin+(j+0.5)*ph->yincr,ph->array[j][i]);
738         }
739         fprintf(fp,"\n");
740     }
741 }
742 
743 void hfetch(FILE *fp)
744 /* warning not yet tested */
745 {
746      shistogram *ph;
747      int id;
748      char title[1024];
749      int type;
750      float xmin,xmax;
751      int xbins;
752      float ymin,ymax;
753      int ybins;
754      int i,j;
755      float entry;
756 
757      fscanf(fp,"%d",&id);
758      fscanf(fp,"%s",title);
759      fscanf(fp,"%d",&type);
760      fscanf(fp,"%f %f %d",&xmin,&xmax,&xbins);
761      fscanf(fp,"%f %f %d",&ymin,&ymax,&ybins);
762 
763      if (hist[id] == NULL && type ==1)
764          ph = hist[id] = hbook1(id,title,xmin,xmax,xbins);
765      if (hist[id] == NULL && type ==2)
766          ph = hist[id] = hbook2(id,title,xmin,xmax,xbins,ymin,ymax,ybins);
767      for (j=0;j<ph->ybins;j++)
768      {
769          for (i=0;i<ph->xbins;i++)
770          {
771              if (type == 1)
772              {
773                  fscanf(fp,"%*f %f",&entry);
774                  hfill1(ph,ph->xmin+i*ph->xincr,entry);
775              }
776              if (type ==2)
777              {
778                  fscanf(fp,"%*f %*f %f",&entry);
779                  hfill2(ph,ph->xmin+i*ph->xincr,ph->ymin+j*ph->yincr,entry);
780              }
781          }
782      }
783 }
784 
785 void hopera()
786 {
787 }
788 
789 void herror(char *str)
790 {
791      printf("error **** %s\n",str);
792      exit();
793 }
794 
795 static double total_chisq(int npar, double *a)
796 {
797       double bh,btot = 0.0;
798       float x;
799       int j;
800 
801       btot = 0.0;
802       for (x=mph->xmin+0.5*(mph->xmax-mph->xmin)/(double)mph->xbins ,j=0;
803            j<mph->xbins;
804            x+=(mph->xmax-mph->xmin)/(double)mph->xbins,j++)
805       {
806            bh = mph->array[0][j]-mfunc(npar,a,x);
807            if(merror_func != NULL) 
808                bh *= 0.5*bh/(merror_func(mph,x)*merror_func(mph,x));
809            else
810                bh *= 0.5*bh;
811            btot += bh;
812       }
813       return(btot);
814 }
815 
816 void hfit(FILE *fp,shistogram *ph,int n,double *a,double *w,
817           double (*fitfunc)(int, double *,float) ,
818           double(*errfunc)(shistogram *,float))
819 {
820        double min;
821        double  simplexmin2(int n, double *x, double *lambda,
822                           double (*funk) ( /* ??? */ ),
823                           double ftol, void (*text) ( /* ??? */ ));
824 
825        mfunc = fitfunc;
826        merror_func = errfunc;
827        mph = ph;
828        min = (float)simplexmin2(n, a, w, total_chisq, 0.000001, NULL);
829 }
830 
831 double **herror_analysis(FILE *fp,shistogram *ph,int n,double *a,
832           double (*fitfunc)(int, double *,float),
833           double (*derfunc)(double *,float,int),
834           double(*errfunc)(shistogram *,float),
835           int exclude)
836 /* a single data point may be selected for exclusion */
837 {
838     double **c,**cinv, *placebo;
839     int i,j,k;
840     double d,x;
841 
842     c = maketarray(n,n,0.0);
843     for (i=0;i<n;i++)
844     {
845         for (j=0;j<n;j++)
846         {
847             c[i][j] = 0.0;
848             for (x=ph->xmin+ph->xincr/2.0,k=0;x<ph->xmax;x+=ph->xincr,k++)
849             {
850                 if (k!=exclude)
851                 {
852                     d = (double)(*errfunc)(ph,x);
853                     d *= d;
854                     c[i][j] += (double)(*derfunc)(a,x,i)*(double)(*derfunc)(a,x,j)/d;
855                 }
856             }
857         }
858     }
859 
860     cinv = hmatinv(n,a,c, placebo);
861     freetarray(c);
862    
863     if (fp!=(FILE *)0)
864     {
865         for (i=0;i<n;i++)
866         {
867             fprintf(fp,"     parameter %d = %f +- %f\n",i,a[i],sqrt(cinv[i][i]));
868         }
869     }
870     return(cinv);
871 }
872 
873 void hsuper(shistogram *ph,int n,double (*func)(int, double*, float),double *a)
874 {
875     int i;
876 
877     if (ph->npar!=n&&ph->par!=NULL)
878     {
879         rfree(ph->par);
880         ph->par = NULL;
881     }
882     if (ph->par==NULL) ph->par = (double *)ralloc(n*sizeof(double));
883     ph->super = func;
884     ph->npar = n;
885     for (i=0;i<n;i++) ph->par[i] = a[i];
886 }
887 
888 void hfit_gauss(FILE *fp,shistogram *ph)
889 {
890     double a[3];
891     double w[3];
892     int i;
893 
894     if (ph==NULL||fp==NULL) return;
895 
896     a[1] = ph->mean;
897     a[2] = (double)sqrt((double)(ph->mean2-ph->mean*ph->mean));
898 
899     a[0] = 0.0;
900 
901     for (i=0;i<ph->xbins;i++)
902        if (ph->array[0][i]>a[0]) a[0]=ph->array[0][i];
903 
904     w[0] = a[0];
905     w[1] = a[2];
906     w[2] = a[2];
907 
908     if (a[2]==0.0) return;
909 
910     hfit(fp,ph,3,a,w,hgaussian,NULL);
911 }
912 
913 double herrdum(shistogram *ph,float x)
914 {
915     double d;
916     double y;
917 
918     y = hfill1(ph,x,0.0);
919 
920     if (y<=0.0) y = 1.0;
921     d = sqrt(y);
922 
923     return(d);
924 }
925 
926 double hgaussian(int n,double *a,float x)
927 {
928    return(a[0]*exp(-(a[1]-x)*(a[1]-x)/(2.0*a[2]*a[2])));
929 }
930 
931 double hgaussder(double *a,float x,int i)
932 {
933     double y;
934     y = hgaussian(3,a,x);
935     if (i==0) return(exp(-(a[1]-x)*(a[1]-x)/(2.0*a[2]*a[2])));
936 
937     if (i==1) return((x-a[1])*y/(a[2]*a[2]));
938 
939     if(i==2) return((x-a[1])*(x-a[1])/(a[2]*a[2]*a[2])*y);
940 
941 }
942 
943 double **hmatinv(int n, double *params, double **a, double *nresult)
944 {
945     double **y, **p, **q, *d, *nparams;
946                 double temp;
947 /*
948     double condition = 0.000001;
949 */
950     double condition = 1.0;
951     double  wmax, thresh;
952     int i,j,k;
953 
954     y = maketarray(n,n,0.0);
955     p = maketarray(n,n,0.0);
956     q = maketarray(n,n,0.0);
957     d = (double *)ralloc(n*sizeof(double));
958     nparams = (double *)ralloc(n*sizeof(double));
959     svd(n,n,a,p,d,q);
960 
961 /*
962     wmax = 0.0;
963     for (j = 0; j < n; j++)
964         if (d[j] > wmax)
965             wmax = d[j];
966     thresh = condition * wmax;
967     for (j = 0; j < n; j++)
968         if (d[j] < thresh)
969             d[j] / thresh*thresh;
970         else
971             d[j] = 1.0/d[j];
972 */
973 /* compute rotated parameter set */
974     for (i = 0; i < n; i++)
975     {
976          nparams[i] = 0.0;
977          for (j=0;j<n;j++)
978          {
979               nparams[i] += params[j]*p[j][i];
980          }
981          nparams[i] *= nparams[i];
982     }
983 /* test new parameter set for significant non-zero */
984                 *nresult = 0.0;
985     for (j = 0; j < n; j++)
986     {
987         temp = d[j];
988         if (d[j]*nparams[j] < condition)
989             d[j] *= nparams[j]*nparams[j]/(condition*condition);
990         else
991             d[j] = 1.0/d[j];
992         *nresult += temp*d[j];
993     }
994 
995     for (i = 0; i < n; i++)
996     {
997          for (j=0;j<n;j++)
998          {
999              for (k=0;k<n;k++)
1000              {
1001                  y[j][k] += d[i]*p[j][i]*q[k][i];
1002              }
1003          }
1004     }
1005 
1006     freetarray(p);
1007     freetarray(q);
1008     rfree(d);
1009     return(y); 
1010 }
1011 
1012 
1013 /*
1014         a.lacey additions
1015 */
1016 
1017 
1018 
1019 void          hbin_dump_plain(FILE *out, shistogram *ph)
1020 {
1021         int     i;
1022         float   c;
1023  
1024         if ((ph == NULL) || (out == NULL)) return;
1025  
1026         c = ph->xmin;
1027         for (i = 0; i < ph->xbins; i++)
1028                 {
1029                         fprintf(out, "%f\t%f\n", c, ph->array[0][i]);
1030                         c += ph->xincr;
1031                 }
1032 }
1033 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.