This is hdata.c in view mode; [Download] [Up]
/*
* hippodata.c - Process data in ntuple into the data which will
* be plotted.
*
* Copyright (C) 1992 The Board of Trustees of The Leland Stanford
* Junior University. All Rights Reserved.
*
* $Id: hdata.c,v 5.0 1993/08/17 21:54:32 rensing Exp $
*
* by Paul Rensing, June 4, 1992
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include "hippo.h"
#include "hutil.h"
#include "hdata.h"
GLOB_QUAL const char hippodata_c_rcsid[] =
"$Id: hdata.c,v 5.0 1993/08/17 21:54:32 rensing Exp $";
/*
* private functions
*/
static void bin1D(display disp, int strtpt);
static void bin2D(display disp, int strtpt);
/* binning for weighted plots */
static void bin1D_w(display disp, int strtpt);
static void bin2D_w(display disp, int strtpt);
int h_cum(display disp, double x, double w, ... )
{
int twoD = 0;
double y = 0.0;
va_list argPtr;
int xbin, ybin;
float *xpt, *last_xpt;
bins_t *bins;
if (disp == NULL) return -1;
if (disp->tuple)
{
h_error("h_cum : This display has an ntuple bound to it.");
return -1;
}
twoD = disp->graphtype == COLORPLOT || disp->graphtype == LEGOPLOT;
if (!twoD && disp->graphtype != HISTOGRAM) return -1;
if (!disp->bins.flags.fixed) disp->bins.flags.fixed = 1;
bins = &disp->bins;
if (twoD)
{
y = w;
va_start(argPtr, w);
w = va_arg(argPtr, double);
va_end(argPtr);
}
xbin = (x - disp->xAxis.low) * ((float) bins->xAxis.nBins)
/ (disp->xAxis.high - disp->xAxis.low);
if (!twoD)
{
if (xbin > bins->xAxis.nBins)
bins->totals[2][0] += w;
else if (xbin < 0)
bins->totals[0][0] += w;
else
{
bins->data[xbin] += w;
bins->variance[xbin] += w * w;
bins->totals[1][0] += w;
bins->xAxis.moments[0] += w;
bins->xAxis.moments[1] += w * x;
bins->xAxis.moments[2] += w * x * x;
}
}
else
{
ybin = (y - disp->yAxis.low) * ((float) bins->yAxis.nBins)
/ (disp->yAxis.high - disp->yAxis.low);
if (xbin > bins->xAxis.nBins)
{
if (ybin > bins->yAxis.nBins)
bins->totals[2][2] += w;
else if (ybin < 0)
bins->totals[2][0] += w;
else
bins->totals[2][1] += w;
}
else if (xbin < 0)
{
if (ybin > bins->yAxis.nBins)
bins->totals[0][2] += w;
else if (ybin < 0)
bins->totals[0][0] += w;
else
bins->totals[0][1] += w;
}
else
{
if (ybin > bins->yAxis.nBins)
bins->totals[1][2] += w;
else if (ybin < 0)
bins->totals[1][0] += w;
else
{
bins->data[xbin * bins->yAxis.nBins + ybin] += w;
bins->variance[xbin * bins->yAxis.nBins + ybin] += w * w;
bins->totals[1][1] += w;
bins->xAxis.moments[0] += w;
bins->xAxis.moments[1] += w * x;
bins->xAxis.moments[2] += w * x * x;
bins->yAxis.moments[0] += w;
bins->yAxis.moments[1] += w * y;
bins->yAxis.moments[2] += w * y * y;
}
}
}
bins->binMin = FLT_MAX;
bins->binMax = -FLT_MAX;
if (twoD)
last_xpt = &bins->data[bins->yAxis.nBins * bins->xAxis.nBins];
else
last_xpt = &bins->data[bins->xAxis.nBins];
for (xpt = bins->data; xpt<last_xpt; xpt++)
{
if (*xpt < bins->binMin) bins->binMin = *xpt;
if (*xpt > bins->binMax) bins->binMax = *xpt;
}
return 0;
}
int h_bin(display disp)
{
int memNeeded = 1, clear=0, spt;
if (disp == NULL) return -1;
if (disp->tuple == NULL) return 0;
if (disp->bins.flags.fixed) return 0;
if (disp->tuple->ndata == 0)
{
disp->bins.ndata = 0;
disp->nt_rev = disp->tuple->rev;
return 0;
}
/* do autoscaling, if necessary */
h_autoScale( disp );
/* scatterplots don't have bins */
if (disp->graphtype == SCATTERPLOT
|| disp->graphtype == XYPLOT
|| disp->graphtype == STRIPCHART
|| disp->graphtype == THREEDSCATTER)
return 0;
/* check if re-binning is actually needed. */
if (!disp->bins.flags.dirty
&& (disp->bins.ndata == disp->tuple->ndata)
&& (disp->nt_rev == disp->tuple->rev)
&& (disp->bins.data != NULL) )
return 0;
/*
* if display parameters have been changed, must reset bins.
*/
if (disp->bins.flags.dirty) clear = 1;
/*
* if revision number has changed, clear bins.
*/
if ( (disp->nt_rev == disp->tuple->rev) ||
(disp->bins.ndata > disp->tuple->ndata) )
clear = 1;
/*
* check how many bins are needed and reallocate if necessary
*/
switch (disp->dim)
{
case 2:
memNeeded *= disp->bins.yAxis.nBins;
case 1:
memNeeded *= disp->bins.xAxis.nBins;
}
if (memNeeded <= 0) return -1;
if (disp->bins.binAlloc < memNeeded)
{
if (disp->bins.data == NULL || disp->bins.binAlloc == 0)
{
disp->bins.data =
(float *) malloc( memNeeded * sizeof(float));
disp->bins.variance =
(float *) malloc( memNeeded * sizeof(float));
}
else
{
disp->bins.data =
(float *) realloc(disp->bins.data,
memNeeded * sizeof(float));
disp->bins.variance =
(float *) realloc(disp->bins.variance,
memNeeded * sizeof(float));
}
if ( (disp->bins.data == NULL) || (disp->bins.variance == NULL) )
{
free(disp->bins.data);
free(disp->bins.variance);
disp->bins.data = disp->bins.variance = NULL;
disp->bins.binAlloc = 0;
return -1;
}
disp->bins.binAlloc = memNeeded;
clear = 1;
}
/*
* clear the bins (if necessary).
*/
if (clear)
{
spt = 0;
/* use memset for speed */
memset((char *)disp->bins.data, (char)0,
memNeeded*sizeof(float) );
memset((char *)disp->bins.variance, (char)0,
memNeeded*sizeof(float) );
disp->bins.binMin = FLT_MAX;
disp->bins.binMax = -FLT_MAX;
memset((char *)disp->bins.totals, (char)0,
sizeof(disp->bins.totals) );
memset((char *)disp->bins.xAxis.moments, (char)0,
sizeof(disp->bins.xAxis.moments) );
memset((char *)disp->bins.yAxis.moments, (char)0,
sizeof(disp->bins.yAxis.moments) );
}
else
spt = disp->bins.ndata;
switch (disp->dim)
{
case 1:
if (disp->binding.weight < 0) bin1D(disp, spt);
else bin1D_w(disp, spt);
break;
case 2:
if (disp->binding.weight < 0) bin2D(disp, spt);
else bin2D_w(disp, spt);
break;
default:
break;
}
disp->bins.ndata = disp->tuple->ndata;
disp->bins.flags.dirty = 0;
disp->nt_rev = disp->tuple->rev;
return 0;
}
static void bin1D(display disp, int strtpt)
{
float xlow = disp->xAxis.low;
float xhigh = disp->xAxis.high;
float bwin = ((float) disp->bins.xAxis.nBins) / (xhigh - xlow);
float *ept;
func_id cutlist = disp->nt_cut;
bins_t *bins = &disp->bins;
/*
* set xpt to point to first nt point to accumulate.
* set last_data to point to last point in nt data array.
*/
float *last_xpt =
&disp->tuple->data[disp->binding.x][disp->tuple->ndata];
register float *xpt = &disp->tuple->data[disp->binding.x][strtpt];
register int i_nt = strtpt;
register int bin;
for ( ; xpt < last_xpt; xpt++, i_nt++ )
{
if (cutlist != NULL)
if (doCuts(disp->tuple->data, i_nt, cutlist)) continue;
if (*xpt >= xhigh)
bins->totals[2][0] += 1.0;
else if (*xpt < xlow)
bins->totals[0][0] += 1.0;
else
{
bin = (*xpt - xlow) * bwin;
bins->data[bin]++;
bins->totals[1][0] += 1.0;
bins->xAxis.moments[0] += 1.0;
bins->xAxis.moments[1] += *xpt;
bins->xAxis.moments[2] += *xpt * *xpt;
}
}
/*
* Look for min and max AFTER accumulating bins.
* Minimum can only be found after all points accumulated, since
* first accumulation gives min. eg. all bins start at 0; first point
* make a bin = 1; this would be the min.
* Also, generally faster. Usually # data points << number of bins,
* so much fewer compares.
*/
last_xpt = &bins->data[disp->bins.xAxis.nBins];
bins->binMin = FLT_MAX;
bins->binMax = -FLT_MAX;
for (xpt = bins->data, ept = bins->variance;
xpt<last_xpt;
xpt++, ept++)
{
if (*xpt < bins->binMin) bins->binMin = *xpt;
if (*xpt > bins->binMax) bins->binMax = *xpt;
*ept = *xpt;
}
}
/*
* binning for 1D hist. with weight. (keep separate for speed)
*/
static void bin1D_w(display disp, int strtpt)
{
float xlow = disp->xAxis.low;
float xhigh = disp->xAxis.high;
float bwin = ((float) disp->bins.xAxis.nBins) / (xhigh - xlow);
func_id cutlist = disp->nt_cut;
bins_t *bins = &disp->bins;
/*
* set xpt to point to first nt point to accumulate.
* set last_data to point to last point in nt data array.
*/
float *last_xpt =
&disp->tuple->data[disp->binding.x][disp->tuple->ndata];
register float *xpt = &disp->tuple->data[disp->binding.x][strtpt];
register int i_nt = strtpt;
register float *wpt =
&disp->tuple->data[disp->binding.weight][strtpt];
register int bin;
for ( ; xpt < last_xpt; xpt++, wpt++, i_nt++ )
{
if (cutlist != NULL)
if (doCuts(disp->tuple->data, i_nt, cutlist)) continue;
if (*xpt >= xhigh) bins->totals[2][0] += *wpt;
else if (*xpt < xlow) bins->totals[0][0] += *wpt;
else
{
bin = (*xpt - xlow) * bwin;
bins->data[bin] += *wpt;
bins->variance[bin] += (*wpt) * (*wpt);
bins->totals[1][0] += *wpt;
bins->xAxis.moments[0] += *wpt;
bins->xAxis.moments[1] += *wpt * *xpt;
bins->xAxis.moments[2] += *wpt * *xpt * *xpt;
}
}
bins->binMin = FLT_MAX;
bins->binMax = -FLT_MAX;
last_xpt = &bins->data[disp->bins.xAxis.nBins];
for (xpt = bins->data; xpt<last_xpt; xpt++)
{
if (*xpt < bins->binMin) bins->binMin = *xpt;
if (*xpt > bins->binMax) bins->binMax = *xpt;
}
return;
}
static void bin2D(display disp, int strtpt)
{
int xbin, ybin;
int nYBins = disp->bins.yAxis.nBins;
float xlow = disp->xAxis.low;
float xhigh = disp->xAxis.high;
float xbwin = ((float) disp->bins.xAxis.nBins) / (xhigh - xlow);
float ylow = disp->yAxis.low;
float yhigh = disp->yAxis.high;
float ybwin = ((float) nYBins) / (yhigh - ylow);
func_id cutlist = disp->nt_cut;
bins_t *bins = &disp->bins;
/*
* set xpt to point to first nt point to accumulate.
* set last_data to point beyond last point in nt data array.
*/
float *last_xpt =
&disp->tuple->data[disp->binding.x][disp->tuple->ndata];
register float *xpt = &disp->tuple->data[disp->binding.x][strtpt];
register float *ypt = &disp->tuple->data[disp->binding.y][strtpt];
register int i_nt = strtpt;
for ( ; xpt<last_xpt; xpt++, ypt++, i_nt++ )
{
if (cutlist != NULL)
if (doCuts(disp->tuple->data, i_nt, cutlist)) continue;
if (*xpt >= xhigh)
{
if (*ypt >= yhigh)
bins->totals[2][2] += 1.0;
else if (*ypt < ylow)
bins->totals[2][0] += 1.0;
else
bins->totals[2][1] += 1.0;
}
else if (*xpt < xlow)
{
if (*ypt >= yhigh)
bins->totals[0][2] += 1.0;
else if (*ypt < ylow)
bins->totals[0][0] += 1.0;
else
bins->totals[0][1] += 1.0;
}
else
{
if (*ypt >= yhigh)
bins->totals[1][2] += 1.0;
else if (*ypt < ylow)
bins->totals[1][0] += 1.0;
else
{
xbin = (*xpt - xlow) * xbwin;
ybin = (*ypt - ylow) * ybwin;
bins->data[xbin * nYBins + ybin]++;
bins->totals[1][1] += 1.0;
bins->xAxis.moments[0] += 1.0;
bins->xAxis.moments[1] += *xpt;
bins->xAxis.moments[2] += *xpt * *xpt;
bins->yAxis.moments[0] += 1.0;
bins->yAxis.moments[1] += *ypt;
bins->yAxis.moments[2] += *ypt * *ypt;
}
}
}
bins->binMin = FLT_MAX;
bins->binMax = -FLT_MAX;
last_xpt = &bins->data[nYBins * disp->bins.xAxis.nBins];
for (xpt = bins->data, ypt=bins->variance; xpt<last_xpt; xpt++, ypt++)
{
if (*xpt < bins->binMin) bins->binMin = *xpt;
if (*xpt > bins->binMax) bins->binMax = *xpt;
*ypt = *xpt;
}
return;
}
static void bin2D_w(display disp, int strtpt)
{
int xbin, ybin;
int nYBins = disp->bins.yAxis.nBins;
float xlow = disp->xAxis.low;
float xhigh = disp->xAxis.high;
float xbwin = ((float) disp->bins.xAxis.nBins) / (xhigh - xlow);
float ylow = disp->yAxis.low;
float yhigh = disp->yAxis.high;
float ybwin = ((float) nYBins) / (yhigh - ylow);
func_id cutlist = disp->nt_cut;
bins_t *bins = &disp->bins;
/*
* set xpt to point to first nt point to accumulate.
* set last_data to point to last point in nt data array.
*/
float *last_xpt =
&disp->tuple->data[disp->binding.x][disp->tuple->ndata];
register float *xpt = &disp->tuple->data[disp->binding.x][strtpt];
register float *ypt = &disp->tuple->data[disp->binding.y][strtpt];
register float *wpt = &disp->tuple->data[disp->binding.weight][strtpt];
register int i_nt = strtpt;
for ( ; xpt<last_xpt; xpt++, ypt++, wpt++, i_nt++ )
{
if (cutlist != NULL)
if (doCuts(disp->tuple->data, i_nt, cutlist)) continue;
if (*xpt >= xhigh)
{
if (*ypt >= yhigh)
bins->totals[2][2] += *wpt;
else if (*ypt < ylow)
bins->totals[2][0] += *wpt;
else
bins->totals[2][1] += *wpt;
}
else if (*xpt < xlow)
{
if (*ypt >= yhigh)
bins->totals[0][2] += *wpt;
else if (*ypt < ylow)
bins->totals[0][0] += *wpt;
else
bins->totals[0][1] += *wpt;
}
else
{
if (*ypt >= yhigh)
bins->totals[1][2] += *wpt;
else if (*ypt < ylow)
bins->totals[1][0] += *wpt;
else
{
xbin = (*xpt - xlow) * xbwin;
ybin = (*ypt - ylow) * ybwin;
bins->data[xbin * nYBins + ybin] += *wpt;
bins->variance[xbin * nYBins + ybin] += (*wpt)*(*wpt);
bins->totals[1][1] += *wpt;
bins->xAxis.moments[0] += *wpt;
bins->xAxis.moments[1] += *wpt * *xpt;
bins->xAxis.moments[2] += *wpt * *xpt * *xpt;
bins->yAxis.moments[0] += *wpt;
bins->yAxis.moments[1] += *wpt * *ypt;
bins->yAxis.moments[2] += *wpt * *ypt * *ypt;
}
}
}
bins->binMin = FLT_MAX;
bins->binMax = -FLT_MAX;
last_xpt = &bins->data[nYBins * disp->bins.xAxis.nBins];
for (xpt = bins->data; xpt<last_xpt; xpt++)
{
if (*xpt < bins->binMin) bins->binMin = *xpt;
if (*xpt > bins->binMax) bins->binMax = *xpt;
}
return;
}
int getXYData( display disp, float **xy, float **xerr,
float **yerr, int **over, struct maxs *m, int *xmin_pt)
{
int i_nt,i_pdata;
float *xp,*yp;
float *last_xpt;
float xlow, xhigh;
float ylow, yhigh;
float *xerrorp = NULL, *yerrorp = NULL;
float *xerrp, *yerrp;
float *xyp;
int *overp, isover;
int ylog, xlog;
func_id cutlist = disp->nt_cut;
if (disp->binding.x < 0 || disp->binding.y < 0 ||
disp->binding.x >= disp->tuple->ndim ||
disp->binding.y >= disp->tuple->ndim )
return -1;
m->xl = FLT_MAX;
m->xh = -FLT_MAX;
m->yl = FLT_MAX;
m->yh = -FLT_MAX;
m->zl = FLT_MAX;
m->zh = -FLT_MAX;
m->xsum = 0.0;
m->x2sum = 0.0;
m->ysum = 0.0;
m->y2sum = 0.0;
m->xysum = 0.0;
ylog = disp->yAxis.flags.log;
xlog = disp->xAxis.flags.log;
if (disp->xAxis.flags.autoScale)
{
xlow = -FLT_MAX;
xhigh = FLT_MAX;
}
else
{
xlow = disp->xAxis.low;
xhigh = disp->xAxis.high;
}
if (disp->yAxis.flags.autoScale)
{
ylow = -FLT_MAX;
yhigh = FLT_MAX;
}
else
{
ylow = disp->yAxis.low;
yhigh = disp->yAxis.high;
}
memset((char *)disp->bins.totals, (char)0, sizeof(disp->bins.totals) );
if ( (*xy = (float *)malloc( 2*disp->tuple->ndata*sizeof(float)))
== NULL )
{
h_error("h_plot has run out of memory.");
return -1;
}
if ( disp->binding.xerror >=0 )
{
if ((*xerr = (float *)malloc( 2*disp->tuple->ndata*sizeof(float)))
== NULL )
{
free(*xy);
h_error("h_plot has run out of memory.");
return -1;
}
}
else *xerr = NULL;
if ( disp->binding.yerror >= 0 )
{
if ( (*yerr = (float *)malloc( 2*disp->tuple->ndata*sizeof(float)))
== NULL )
{
free(xy);
free(xerr);
h_error("h_plot has run out of memory.");
return -1;
}
}
else *yerr = NULL;
if ( (*over = (int *)malloc( (disp->tuple->ndata+1)*sizeof(int)))
== NULL )
{
free(xy);
free(xerr);
free(yerr);
h_error("h_plot has run out of memory.");
return -1;
}
last_xpt = &disp->tuple->data[disp->binding.x][disp->tuple->ndata];
xp = &disp->tuple->data[disp->binding.x][0];
yp = &disp->tuple->data[disp->binding.y][0];
if (disp->binding.xerror >=0 )
xerrorp = &disp->tuple->data[disp->binding.xerror][0];
if (disp->binding.yerror >=0 )
yerrorp = &disp->tuple->data[disp->binding.yerror][0];
xerrp = *xerr;
yerrp = *yerr;
xyp = *xy;
overp = *over;
for (i_pdata=i_nt=0;
xp<last_xpt;
xp++, yp++, i_nt++, xerrorp++, yerrorp++ )
{
if (cutlist != NULL)
if (doCuts(disp->tuple->data, i_nt, cutlist)) continue;
isover = 0;
if (*xp > xhigh)
{
if (*yp > yhigh)
disp->bins.totals[2][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[2][0] += 1.0;
else
disp->bins.totals[2][1] += 1.0;
}
else if (*xp < xlow)
{
if (*yp > yhigh)
disp->bins.totals[0][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[0][0] += 1.0;
else
disp->bins.totals[0][1] += 1.0;
}
else
{
if (*yp > yhigh)
disp->bins.totals[1][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[1][0] += 1.0;
else
disp->bins.totals[1][1] += 1.0;
}
if (!xlog)
*xyp++ = *xp;
else
{
if (*xp > 0.0)
*xyp++ = log10(*xp);
else
{
isover = 1;
*xyp++ = -FLT_MAX;
}
}
if (*xp > xhigh || *xp < xlow) isover = 1;
if (!ylog)
*xyp++ = *yp;
else
{
if (*yp > 0.0)
*xyp++ = log10(*yp);
else
{
isover = 1;
*xyp++ = -FLT_MAX;
}
}
if (*yp > yhigh || *yp < ylow ) isover = 1;
if (! isover)
{
if (*xp > m->xh) m->xh = *xp;
if (*xp < m->xl)
{
m->xl = *xp;
*xmin_pt = i_pdata;
}
if (*yp > m->yh) m->yh = *yp;
if (*yp < m->yl) m->yl = *yp;
m->xsum += *xp;
m->x2sum += (*xp) * (*xp);
m->ysum += *yp;
m->y2sum += (*yp) * (*yp);
m->xysum += (*xp) * (*yp);
}
if (xerrp != NULL)
{
*xerrp = *xp + *xerrorp;
if (!isover && *xerrp > m->xh) m->xh = *xerrp;
if (xlog)
{
if (*xerrp > 0.0)
*xerrp = log10(*xerrp);
else
{
*xerrp = -FLT_MAX;
isover = 1;
}
}
if (*xerrp > xhigh || *xerrp < xlow) isover = 1;
xerrp++;
*xerrp = *xp - *xerrorp;
if (!isover && *xerrp < m->xl &&
(!xlog || (xlog && *xerrp>0.0)) ) m->xl = *xerrp;
if (xlog)
{
if (*xerrp > 0.0)
*xerrp = log10(*xerrp);
else
{
*xerrp = -FLT_MAX;
isover = 1;
}
}
if (*xerrp > xhigh || *xerrp < xlow) isover = 1;
xerrp++;
}
if (yerrp != NULL)
{
*yerrp = *yp + *yerrorp;
if (!isover && *yerrp > m->yh) m->yh = *yerrp;
if (ylog)
{
if (*yerrp > 0.0)
*yerrp = log10(*yerrp);
else
{
*yerrp = -FLT_MAX;
isover = 1;
}
}
if (*yerrp > yhigh || *yerrp < ylow ) isover = 1;
yerrp++;
*yerrp = *yp - *yerrorp;
if (!isover && *yerrp < m->yl &&
(!ylog || (ylog && *yerrp > 0.0)) )
m->yl = *yerrp;
if (ylog)
{
if (*yerrp > 0.0)
*yerrp = log10(*yerrp);
else
{
*yerrp = -FLT_MAX;
isover = 1;
}
}
if (*yerrp > yhigh || *yerrp < ylow ) isover = 1;
yerrp++;
}
if (isover)
{
*overp++ = i_pdata;
isover = 1;
}
i_pdata++; /* gets incremented if row not cut */
}
/* indicate last item in over */
*overp = i_pdata;
return i_pdata;
}
int getHistoData( display disp, float **xy, float **xylego,
float **xerr, float **yerr, int **over,
struct maxs *m)
{
int i;
float binwidth, binside;
float *yp;
float *last;
float *xyp;
int *overp;
float xlow,xhigh;
float ylow, yhigh;
float *yerrorp, *xerrp, *yerrp;
float *xylegop;
float def_bin;
int ylog, isover;
if (disp->bins.data == NULL) return -1;
ylog = disp->yAxis.flags.log;
if (disp->xAxis.flags.autoScale)
{
xlow = -FLT_MAX;
xhigh = FLT_MAX;
}
else
{
xlow = disp->xAxis.low;
xhigh = disp->xAxis.high;
}
if (disp->yAxis.flags.autoScale)
{
ylow = -FLT_MAX;
yhigh = FLT_MAX;
}
else
{
ylow = disp->yAxis.low;
yhigh = disp->yAxis.high;
}
m->xl = FLT_MAX;
m->xh = -FLT_MAX;
m->yl = FLT_MAX;
m->yh = -FLT_MAX;
m->zl = FLT_MAX;
m->zh = -FLT_MAX;
m->xsum = 0.0;
m->x2sum = 0.0;
m->ysum = 0.0;
m->y2sum = 0.0;
m->xysum = 0.0;
if ( (*xy = (float *)malloc( 2*disp->bins.xAxis.nBins*sizeof(float)))
== NULL )
{
h_error("h_plot has run out of memory.");
return -1;
}
if ( (*xylego = (float *)malloc( 4*(disp->bins.xAxis.nBins+2)*sizeof(float)))
== NULL )
{
free(*xy);
h_error("h_plot has run out of memory.");
return -1;
}
if ( disp->drawtype&ERRBAR )
{
if ( (*xerr = (float *)malloc( 2*disp->bins.xAxis.nBins*sizeof(float)))
== NULL )
{
free(*xy);
free(*xylego);
h_error("h_plot has run out of memory.");
return -1;
}
}
else *xerr = NULL;
if ( disp->drawtype&ERRBAR )
{
if ( (*yerr = (float *)malloc( 2*disp->bins.xAxis.nBins*sizeof(float)))
== NULL )
{
free(*xy);
free(*xylego);
free(*xerr);
h_error("h_plot has run out of memory.");
return -1;
}
}
else *yerr = NULL;
if ( (*over = (int *)malloc( (disp->bins.xAxis.nBins+1)*sizeof(int)))
== NULL )
{
free(*xy);
free(*xylego);
free(*xerr);
free(*yerr);
h_error("h_plot has run out of memory.");
return -1;
}
last = &( disp->bins.data[disp->bins.xAxis.nBins] );
yp = disp->bins.data;
yerrorp = disp->bins.variance;
i = 0;
xerrp = *xerr;
yerrp = *yerr;
xyp = *xy;
xylegop = *xylego;
overp = *over;
binside = disp->xAxis.low;
binwidth = (disp->xAxis.high - disp->xAxis.low)/disp->bins.xAxis.nBins;
/*
* connect the bins to the walls of the plot.
* this is only useful if we allow plot axis to be bigger
* than the bin limits.
*/
if (disp->yAxis.low > 0.0) def_bin = disp->yAxis.low;
else def_bin = 0.0;
if (ylog)
if (def_bin > 0.0)
def_bin = log10(def_bin);
else
def_bin = -FLT_MAX;
*xylegop++ = disp->xAxis.low;
*xylegop++ = def_bin;
*xylegop++ = binside;
*xylegop++ = def_bin;
for ( ; yp<last; yp++, yerrorp++ )
{
isover = 0;
*xyp++ = binside + binwidth/2.0;
if (!ylog)
*xyp++ = *yp;
else
{
if (*yp > 0.0)
*xyp++ = log10(*yp);
else
{
isover = 1;
*xyp++ = -FLT_MAX;
}
}
if (*yp > yhigh || *yp < ylow) isover = 1;
if (!isover)
{
if (*yp > m->yh) m->yh = *yp;
if (*yp < m->yl && (!ylog || (ylog && *yp>0.0)) )
m->yl = *yp;
}
*xylegop++ = binside;
if (xerrp != NULL) *xerrp++ = binside;
*xylegop = *yp;
if (ylog)
{
if (*xylegop > 0.0)
*xylegop = log10(*xylegop);
else
*xylegop = -FLT_MAX;
}
xylegop++;
binside += binwidth;
*xylegop++ = binside;
*xylegop = *(xylegop-2);
xylegop++;
if (xerrp != NULL) *xerrp++ = binside;
if (yerrp != NULL)
{
*yerrp = *yp + sqrt(*yerrorp);
if (!isover && *yerrp > m->yh) m->yh = *yerrp;
if (ylog)
{
if (*yerrp > 0.0)
*yerrp = log10(*yerrp);
else
{
*yerrp = -FLT_MAX;
isover = 1;
}
}
if (*yerrp > yhigh || *yerrp < ylow) isover = 1;
yerrp++;
*yerrp = *yp - sqrt(*yerrorp);
if (!isover && *yerrp < m->yl &&
!( ylog && *yerrp<=0.0))
m->yl = *yerrp;
if (ylog)
{
if (*yerrp > 0.0)
*yerrp = log10(*yerrp);
else
{
*yerrp = -FLT_MAX;
isover = 1;
}
}
if (*yerrp > yhigh || *yerrp < ylow) isover = 1;
yerrp++;
}
if (isover)
{
isover = 1;
*overp++ = i;
}
i++;
}
/* connect to the walls */
*xylegop++ = binside;
*xylegop++ = def_bin;
*xylegop++ = disp->xAxis.high;
*xylegop++ = def_bin;
/* indicate last item in over. */
*overp = i;
return i;
}
int getScatData( display disp, float **xy, struct maxs *m )
{
int i_nt;
float *xp,*yp;
float *last_xpt;
float *xyp;
float xlow, xhigh;
float ylow, yhigh;
int xlog, ylog;
func_id cutlist = disp->nt_cut;
if (disp->binding.x < 0 || disp->binding.y < 0 ||
disp->binding.x >= disp->tuple->ndim ||
disp->binding.y >= disp->tuple->ndim )
return -1;
xlog = disp->xAxis.flags.log;
ylog = disp->yAxis.flags.log;
m->xl = FLT_MAX;
m->xh = -FLT_MAX;
m->yl = FLT_MAX;
m->yh = -FLT_MAX;
m->zl = FLT_MAX;
m->zh = -FLT_MAX;
m->xsum = 0.0;
m->x2sum = 0.0;
m->ysum = 0.0;
m->y2sum = 0.0;
m->xysum = 0.0;
if (disp->xAxis.flags.autoScale)
{
xlow = -FLT_MAX;
xhigh = FLT_MAX;
}
else
{
xlow = disp->xAxis.low;
xhigh = disp->xAxis.high;
}
if (disp->yAxis.flags.autoScale)
{
ylow = -FLT_MAX;
yhigh = FLT_MAX;
}
else
{
ylow = disp->yAxis.low;
yhigh = disp->yAxis.high;
}
if (xlog && xlow < FLT_MIN) xlow = FLT_MIN;
if (ylog && ylow < FLT_MIN) ylow = FLT_MIN;
if (xhigh <= xlow) xhigh = xlow + 1.0;
if (yhigh <= ylow) yhigh = ylow + 1.0;
if ( (*xy = (float *)malloc( 2*disp->tuple->ndata*sizeof(float)))
== NULL )
{
h_error("h_plot has run out of memory.");
return -1;
}
memset((char *)disp->bins.totals, (char)0, sizeof(disp->bins.totals) );
last_xpt = &disp->tuple->data[disp->binding.x][disp->tuple->ndata];
xp = &disp->tuple->data[disp->binding.x][0];
yp = &disp->tuple->data[disp->binding.y][0];
i_nt = 0;
xyp = *xy;
for ( ; xp<last_xpt; xp++, yp++, i_nt++ )
{
if (cutlist != NULL)
if (doCuts(disp->tuple->data, i_nt, cutlist)) continue;
if (*xp >= xhigh)
{
if (*yp > yhigh)
disp->bins.totals[2][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[2][0] += 1.0;
else
disp->bins.totals[2][1] += 1.0;
}
else if (*xp < xlow)
{
if (*yp >= yhigh)
disp->bins.totals[0][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[0][0] += 1.0;
else
disp->bins.totals[0][1] += 1.0;
}
else
{
if (*yp > yhigh)
disp->bins.totals[1][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[1][0] += 1.0;
else
{
disp->bins.totals[1][1] += 1.0;
if (!xlog)
*xyp++ = *xp;
else
*xyp++ = log10(*xp);
if (!ylog)
*xyp++ = *yp;
else
*xyp++ = log10(*yp);
if (*xp > m->xh) m->xh = *xp;
if (*xp < m->xl) m->xl = *xp;
if (*yp > m->yh) m->yh = *yp;
if (*yp < m->yl) m->yl = *yp;
m->xsum += *xp;
m->x2sum += (*xp) * (*xp);
m->ysum += *yp;
m->y2sum += (*yp) * (*yp);
m->xysum += (*xp) * (*yp);
}
}
}
return (xyp - *xy)/2;
}
int getScatData3D( display disp, float **xyz, struct maxs *m )
{
int i_nt;
float *xp,*yp,*zp;
float *last_xpt;
float *xyzp;
float xlow, xhigh;
float ylow, yhigh;
float zlow, zhigh;
int xlog, ylog, zlog;
func_id cutlist = disp->nt_cut;
if (disp->binding.x < 0 || disp->binding.y < 0 || disp->binding.z < 0 ||
disp->binding.x >= disp->tuple->ndim ||
disp->binding.y >= disp->tuple->ndim ||
disp->binding.z >= disp->tuple->ndim )
return -1;
xlog = disp->xAxis.flags.log;
ylog = disp->yAxis.flags.log;
zlog = disp->zAxis.flags.log;
m->xl = FLT_MAX;
m->xh = -FLT_MAX;
m->yl = FLT_MAX;
m->yh = -FLT_MAX;
m->zl = FLT_MAX;
m->zh = -FLT_MAX;
m->xsum = 0.0;
m->x2sum = 0.0;
m->ysum = 0.0;
m->y2sum = 0.0;
m->xysum = 0.0;
if (disp->xAxis.flags.autoScale)
{
xlow = -FLT_MAX;
xhigh = FLT_MAX;
}
else
{
xlow = disp->xAxis.low;
xhigh = disp->xAxis.high;
}
if (disp->yAxis.flags.autoScale)
{
ylow = -FLT_MAX;
yhigh = FLT_MAX;
}
else
{
ylow = disp->yAxis.low;
yhigh = disp->yAxis.high;
}
if (disp->zAxis.flags.autoScale)
{
zlow = -FLT_MAX;
zhigh = FLT_MAX;
}
else
{
zlow = disp->zAxis.low;
zhigh = disp->zAxis.high;
}
if (xlog && xlow < FLT_MIN) xlow = FLT_MIN;
if (ylog && ylow < FLT_MIN) ylow = FLT_MIN;
if (zlog && zlow < FLT_MIN) zlow = FLT_MIN;
if (xhigh <= xlow) xhigh = xlow + 1.0;
if (yhigh <= ylow) yhigh = ylow + 1.0;
if (zhigh <= zlow) zhigh = zlow + 1.0;
if ( (*xyz = (float *)malloc( 3*disp->tuple->ndata*sizeof(float)))
== NULL )
{
h_error("h_plot has run out of memory.");
return -1;
}
memset((char *)disp->bins.totals, (char)0, sizeof(disp->bins.totals) );
last_xpt = &disp->tuple->data[disp->binding.x][disp->tuple->ndata];
xp = &disp->tuple->data[disp->binding.x][0];
yp = &disp->tuple->data[disp->binding.y][0];
zp = &disp->tuple->data[disp->binding.z][0];
i_nt = 0;
xyzp = *xyz;
for ( ; xp<last_xpt; xp++, yp++, zp++, i_nt++ )
{
if (cutlist != NULL)
if (doCuts(disp->tuple->data, i_nt, cutlist)) continue;
if (*xp >= xhigh)
{
if (*yp > yhigh)
disp->bins.totals[2][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[2][0] += 1.0;
else
disp->bins.totals[2][1] += 1.0;
}
else if (*xp < xlow)
{
if (*yp >= yhigh)
disp->bins.totals[0][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[0][0] += 1.0;
else
disp->bins.totals[0][1] += 1.0;
}
else /* x in range ... */
{
if (*yp > yhigh)
disp->bins.totals[1][2] += 1.0;
else if (*yp < ylow)
disp->bins.totals[1][0] += 1.0;
else /* y in range ... */
{
/*** Note: we should expand this to include the z totals....someday ***/
if (*zp > zhigh) {}
/****disp->bins.totals[...][...] += 1.0;****/
else if (*zp < zlow) {}
/****disp->bins.totals[...][...] += 1.0;****/
else /* z in range ... */
{
disp->bins.totals[1][1] += 1.0;
if (!xlog)
*xyzp++ = *xp;
else
*xyzp++ = log10(*xp);
if (!ylog)
*xyzp++ = *yp;
else
*xyzp++ = log10(*yp);
if (!zlog)
*xyzp++ = *zp;
else
*xyzp++ = log10(*zp);
if (*xp > m->xh) m->xh = *xp;
if (*xp < m->xl) m->xl = *xp;
if (*yp > m->yh) m->yh = *yp;
if (*yp < m->yl) m->yl = *yp;
if (*zp > m->zh) m->zh = *zp;
if (*zp < m->zl) m->zl = *zp;
}
}
}
}
return (xyzp - *xyz)/3;
}
int h_fillWithFitData( display dsp, ntuple nt, int unbinned )
{
int rc = 0;
if (dsp == NULL || nt == NULL) return -1;
if (nt->ndim < 4 ||
((dsp->graphtype == LEGOPLOT || dsp->graphtype == COLORPLOT ||
dsp->graphtype == THREEDSCATTER) && nt->ndim < 6))
return -1;
h_clrNt( nt );
switch (dsp->graphtype)
{
case (int)HISTOGRAM:
if (unbinned)
{
float *xd, *wd = NULL;
int i;
int useWeight = dsp->binding.weight >= 0;
func_id cutlist = dsp->nt_cut;
float w = 1.0;
xd = dsp->tuple->data[dsp->binding.x];
if (useWeight)
wd = dsp->tuple->data[dsp->binding.weight];
for (i=0; i<dsp->tuple->ndata && rc==0; i++, xd++, wd++)
{
if (cutlist != NULL)
if (doCuts(dsp->tuple->data, i, cutlist)) continue;
if (*xd < dsp->xAxis.low
|| *xd > dsp->xAxis.high) continue;
if (useWeight) w = *wd;
rc = h_fill( nt, *xd, w );
}
break;
}
else
{
float *yd, *yed;
float x, xe;
float *last;
h_bin( dsp );
yd = dsp->bins.data;
yed = dsp->bins.variance;
xe = (dsp->xAxis.high - dsp->xAxis.low)
/ (2.0 * dsp->bins.xAxis.nBins);
x = dsp->xAxis.low + xe;
last = &dsp->bins.data[dsp->bins.xAxis.nBins];
for (; yd < last && rc == 0; yd++, yed++, x+=2.0*xe )
rc = h_fill( nt, x, xe, *yd, sqrt(*yed) );
break;
}
case (int)XYPLOT:
case (int)STRIPCHART:
{
float *xd, *yd, *xed=NULL, *yed=NULL;
float xerr=1.0, yerr=1.0;
int useXerr, useYerr;
int i;
func_id cutlist = dsp->nt_cut;
xd = dsp->tuple->data[dsp->binding.x];
yd = dsp->tuple->data[dsp->binding.y];
useXerr = (dsp->binding.xerror >= 0);
useYerr = (dsp->binding.yerror >= 0);
if (useXerr)
xed = dsp->tuple->data[dsp->binding.xerror];
if (useYerr)
yed = dsp->tuple->data[dsp->binding.yerror];
for (i=0; i<dsp->tuple->ndata && rc == 0;
i++, xd++, yd++, xed++, yed++)
{
if (cutlist != NULL)
if (doCuts(dsp->tuple->data, i, cutlist)) continue;
if (*xd < dsp->xAxis.low || *xd > dsp->xAxis.high) continue;
if (useXerr) xerr = *xed;
if (useYerr) yerr = *yed;
rc = h_fill( nt, *xd, xerr, *yd, yerr );
}
break;
}
case COLORPLOT:
case LEGOPLOT:
/* if unbinned, fall through to scatterplot */
if (!unbinned)
{
float *zd, *zed;
float x, xe;
float y, ye;
float *last, *last2;
h_bin( dsp );
zd = dsp->bins.data;
zed = dsp->bins.variance;
xe = (dsp->bins.xAxis.high - dsp->bins.xAxis.low)
/ (2.0 * dsp->bins.xAxis.nBins);
ye = (dsp->bins.yAxis.high - dsp->bins.yAxis.low)
/ (2.0 * dsp->bins.yAxis.nBins);
x = dsp->bins.xAxis.low + xe;
last = &dsp->bins.data[dsp->bins.xAxis.nBins
* dsp->bins.yAxis.nBins];
for (; zd < last && rc == 0; x+=2.0*xe)
{
y = dsp->bins.yAxis.low + ye;
last2 = zd + dsp->bins.yAxis.nBins;
for (; zd<last2 && rc==0; zd++, zed++, y+=2.0*ye )
rc = h_fill( nt, x, xe, y, ye, *zd, *zed );
}
break;
}
case SCATTERPLOT:
{
float *xd, *yd;
int i;
func_id cutlist = dsp->nt_cut;
xd = dsp->tuple->data[dsp->binding.x];
yd = dsp->tuple->data[dsp->binding.y];
for (i=0; i<dsp->tuple->ndata && rc==0; i++, xd++, yd++)
{
if (cutlist != NULL)
if (doCuts(dsp->tuple->data, i, cutlist)) continue;
if (*xd < dsp->xAxis.low || *xd > dsp->xAxis.high) continue;
if (*yd < dsp->yAxis.low || *yd > dsp->yAxis.high) continue;
rc = h_fill( nt, *xd, 1.0, *yd, 1.0 );
}
break;
}
case THREEDSCATTER:
{
float *xd, *yd, *zd;
int i;
func_id cutlist = dsp->nt_cut;
xd = dsp->tuple->data[dsp->binding.x];
yd = dsp->tuple->data[dsp->binding.y];
zd = dsp->tuple->data[dsp->binding.z];
for (i=0; i<dsp->tuple->ndata && rc==0; i++, xd++, yd++, zd++)
{
if (cutlist != NULL)
if (doCuts(dsp->tuple->data, i, cutlist)) continue;
if (*xd < dsp->xAxis.low || *xd > dsp->xAxis.high) continue;
if (*yd < dsp->yAxis.low || *yd > dsp->yAxis.high) continue;
if (*zd < dsp->zAxis.low || *zd > dsp->zAxis.high) continue;
rc = h_fill( nt, *xd, 1.0, *yd, 1.0, *zd, 1.0 );
}
break;
}
}
return rc;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.