A new version of PAW (2.05/20) has been released today 18 October together with CERNLIB 94B. The source files, binaries and libraries will be available, via anonymous ftp, on the asis system.
The PAW team is:
* NTUPLE/UWFUNC IDN FNAME [ CHOPT ]
IDN C 'Ntuple Identifier'
FNAME C 'File name'
CHOPT C 'Options' D=' '
Possible CHOPT values are:
' ' Generate the FORTRAN skeleton of a selection function.
E Present the selection function in the local editor.
P Code to print events is generated (not valid for new Ntuples).
T Names of the Ntuple variables are generated in DATA statements (not
valid for new Ntuples).
To generate the FORTRAN skeleton of a selection function or the INCLUDE
file with the columns declaration.
A FORTRAN function is generated if the FNAME is of the form, xxx.f,
xxx.for, xxx.fortran. Otherwise an INCLUDE file is generated. Example: If
Ntuple ID=30 has variable names [X,Y,Z,ETOT,EMISS,etc] then:
NTUPLE/UWFUNC 30 SELECT.FOR will generate the file SELECT.FOR with:
FUNCTION SELECT(XDUMMY)
COMMON/PAWIDN/IDNEVT,VIDN1,VIDN2,VIDN3,X,Y,Z,ETOT,EMISS,etc
SELECT=1.
END
Then using the command EDIT one can modify this file which could then
look something like (IDNEVT is the event number):
FUNCTION SELECT(XDUMMY)
COMMON/PAWIDN/IDNEVT,VIDN1,VIDN2,VIDN3,X,Y,Z,ETOT,EMISS,etc
IF(X**2+Y**2.GT.Z**2.OR.ETOT.GT.20.)THEN
SELECT=1.
ELSE
SELECT=0.
ENDIF
END
If in a subsequent command NTUPLE/PLOT, the selection function SELECT is
used, then:
If NTUPLE/PLOT 30.ETOT SELECT.FOR
VIDN1=ETOT
If NTUPLE/PLOT 30.SQRT(X**2+Y**2)%(ETOT-EMISS)
VIDN1=ETOT-EMISS
VIDN2=SQRT(X**2+Y**2)
NTUPLE/UWFUNC 30 SELECT.INC will generate an include file. This include
file may be referenced in a selection function in the following way:
FUNCTION SELECT(XDUMMY)
include 'select.inc'
SELECT=1.
IF(X.LE.Y)SELECT=0.
END
* NTUPLE/SCAN IDN [ UWFUNC NEVENT IFIRST OPTION VARLIS ]
IDN C 'Ntuple Identifier'
UWFUNC C 'User cut function' D='0'
NEVENT I 'Number of events' D=99999999
IFIRST I 'First event' D=1
OPTION C 'Options' D=' '
VARLIS C 'Names of the NVARS variables to scan' D=' ' Vararg
Possible OPTION values are:
' '
S Graphical scan (spider plot).
' ' Alphanumeric output of the Ntuple.
A Used with "S" it displays the average spider.
Scan the entries of an Ntuple subject to user cuts. Scan the variables
for NEVENT events starting at IFIRST, requiring that the events satisfy
cut UWFUNC. In the case of Alphanumeric output Up to 8 variables may be
scanned, the default is to scan the first 8 variables.
When the option S (Spider plot) is specified, each event is presented in
a graphical form (R versus PHI plot) to give a multi dimensional view of
the event. Each variable is represented on a separate axis with a scale
ranging from the minimum to the maximum value of the variable. A line
joins all the current points on every axis where each point corresponds
to the current value of the variable. When the HCOL parameter is
specified (eg SET HCOL 1002) a fill area is drawn.
VARLIS may contain a list of the original variables, expressions of the
original variables or/and ranges of variables. A range can be given in
the following form:
: means all variables (default).
var1:var2 means from variable var1 to variable var2 included.
var1: means from variable var1 to the last.
:var2 means from variable 1 to variable var2
For example, if IDN=30 has the 3 variables X,Y,Z,U,V,W one can do:
PAW > scan 30
PAW > scan 30 option=s
each event is drawn as a spider plot.
PAW > scan 30 option=sa
each event is drawn as a spider plot and the average spider
plot is also drawn.
PAW > scan 30 option=s X:Z W
PAW > scan 30 z>10
PAW > scan 30 z>10 ! ! ! z abs(x) y+z x func.for
where func.for is a COMIS function returning an expression
of the original variables. This function func.for may be
generated automatically by the PAW command:
PAW > uwfunc 30 func.for
* HISTOGRAM/FIT ID FUNC [ CHOPT NP PAR STEP PMIN PMAX ERRPAR ]
ID C 'Histogram Identifier'
FUNC C 'Function name' D=' '
CHOPT C 'Options' D=' '
NP I 'Number of parameters' D=0 R=0:34
PAR C 'Vector of parameters'
STEP C 'Vector of steps size'
PMIN C 'Vector of lower bounds'
PMAX C 'Vector of upper bounds'
ERRPAR C 'Vector of errors on parameters'
Possible CHOPT values are:
' ' Do the fit, plot the result and print the parameters.
0 Do not plot the result of the fit. By default the fitted function is
drawn unless the option 'N' below is specified.
N Do not store the result of the fit bin by bin with the histogram. By
default the function is calculated at the middle of each bin and the
fit results stored with the histogram data structure.
Q Quiet mode. No print
V Verbose mode. Results after each iteration are printed By default
only final results are printed.
B Some or all parameters are bounded. The vectors STEP,PMIN,PMAX must
be specified. Default is: All parameters vary freely.
L Use Log Likelihood. Default is chisquare method.
D The user is assumed to compute derivatives analytically using the
routine HDERIV. By default, derivatives are computed numerically.
W Sets weights equal to 1. Default weights taken from the square root
of the contents or from HPAKE/HBARX (PUT/ERRORS). If the L option is
given (Log Likelihood), bins with errors=0 are excluded of the fit.
M The interactive Minuit is invoked. (see Application HMINUIT below).
E Performs a better Error evaluation (MIGRAD + HESSE + MINOS).
U User function value is taken from /HCFITD/FITPAD(24),FITFUN.
K Keep the settings of Application HMINUIT for a subsequent command.
Fit a user defined (and parameter dependent) function to a histogram ID
(1-Dim or 2-Dim) in the specified range. FUNC may be:
A- The name of a file which contains the user defined
function to be minimized. Function name and file name
must be the same. For example file FUNC.FOR is:
FUNCTION FUNC(X) or FUNC(X,Y) for a 2-Dim histogram
COMMON/PAWPAR/PAR(2)
FUNC=PAR(1)*X +PAR(2)*EXP(-X)
END
Ex: His/fit 10 func.for ! 5 par
When the option U is given, the file FUNC.FOR should look like:
FUNCTION FUNC(X) or FUNC(X,Y) for a 2-Dim histogram
DOUBLE PRECISION FITPAD(24),FITFUN
COMMON/HCFITD/FITPAD,FITFUN
FITFUN=FITPAD(1)*X +FITPAD(2)*EXP(-X)
FUNC=FITFUN
END
B- One of the following keywords (1-Dim only):
G : to fit Func=par(1)*exp(-0.5*((x-par(2))/par(3))**2)
E : to fit Func=exp(par(1)+par(2)*x)
Pn: to fit Func=par(1)+par(2)*x+par(3)*x**2......+par(n+1)*x**n
Ex: His/fit 10 g
C- A combination of the keywords in B with the 2 operators + or *.
Ex: His/Fit 10 p4+g ! 8 par
His/Fit 10 p2*g+g ! 9 par
Note that in this case, the order of parameters in PAR must
correspond to the order of the basic functions.
For example, in the first case above, par(1:5) apply to
the polynomial of degree 4 and par(6:8) to the gaussian while
in the second case par(1:3) apply to the polynomial of degree 2,
par(4:6) to the first gaussian and par(7:9) to the second gaussian.
Blanks are not allowed in the expression.
For cases A and C, before the execution of this command, the vector PAR
must be filled (via Vector/Input) with the initial values. For case B,
if NP is set to 0, then the initial values of PAR will be calculated
automatically. After the fit, the vector PAR contains the new values of
parameters. If the vector ERRPAR is given, it will contain the errors on
the fitted parameters. A bin range may be specified with ID.
Ex. Histo/Fit 10(25:56).
When the Histo/it command is used in a macro, it might be convenient to
specify MINUIT directives in the macro itself via the Application HMINUIT
as described in this example:
Macro fit
Application HMINUIT exit
name 1 par_name1
name 2 par_name2
migrad
improve
exit
Histo/fit id fitfun.f M
Return
NOTE: The computation of the errors is based on a proposal by Stephane Coutu. below is a copy of the email with the proposal ONLY options 'S' and 'I' are implemented
I realized that there is another case where this kind of trouble occurs: if a bin has N data points all with the same value Y (especially possible when dealing with integers), the spread in Y for that bin is zero, and the uncertainty assigned is also zero, and the bin is ignored in making subsequent fits. If SQRT(Y) was the correct error in the case above, then SQRT(Y)/SQRT(N) would be the correct error here. In fact, any bin with non-zero number of entries N but with zero spread should have an uncertainty SQRT(Y)/SQRT(N).
Now, is SQRT(Y)/SQRT(N) really the correct uncertainty? I believe that it is only in the case where the Y variable is some sort of counting statistics, following a Poisson distribution. This should probably be set as the default case. However, Y can be any variable from an original NTUPLE, not necessarily distributed "Poissonly", and perhaps extra options could be offered with the command: PROFILE id title ncx xmin xmax ymin ymax [ chopt ] to allow the user to choose how errors are calculated. We could have, for example:
CHOPT
' ' (Default) Errors are Spread/SQRT(N) for Spread.ne.0. ,
" " SQRT(Y)/SQRT(N) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
'S' Errors are Spread for Spread.ne.0. ,
" " SQRT(Y) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
'I' Errors are Spread/SQRT(N) for Spread.ne.0. ,
" " 1./SQRT(12.*N) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
The third case above corresponds to Integer Y values for which the uncertainty is +-0.5, with the assumption that the probability that Y takes any value between Y-0.5 and Y+0.5 is uniform (the same argument goes for Y uniformly distributed between Y and Y+1); this would be useful if Y is an ADC measurement, for example. Other, fancier options would be possible, at the cost of adding one more parameter to the PROFILE command. For example, if all Y variables are distributed according to some known Gaussian of standard deviation Sigma, then:
'G' Errors are Spread/SQRT(N) for Spread.ne.0. ,
" " Sigma/SQRT(N) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
For example, this would be useful when all Y's are experimental quantities measured with the same instrument with precision Sigma.
Stephane Coutu
coutu@roo.physics.lsa.umich.edu
* HISTOGRAM/GET_VECT/REBIN ID X Y EX EY [ N IFIRST ILAST CHOPT ]
ID C 'Histogram Identifier'
X C 'Name of vector X'
Y C 'Name of vector Y'
EX C 'Name of vector EX'
EY C 'Name of vector EY'
N I 'Number of elements to fill' D=100
IFIRST I 'First bin' D=1
ILAST I 'Last bin' D=100
CHOPT C 'Option' D=' '
Possible CHOPT values are:
N Do not normalize values in Y
The specified channels of the 1-Dim histogram ID are cumulated (rebinned)
into new bins. The final contents of the new bin is the average of the
original bins by default. If the option N is given, the final contents of
the new bin is the sum of the original bins. Get contents and errors
into vectors, grouping bins. Bin width and centers are also extracted.
Allow to combine 2, 3 or more bins into one.
E.g.: REBIN 110 X Y EX EY 25 11 85
will group by 3 channels 11 to 85 and return
new abscissa, contents and errors.
Errors in X are equal to 1.5*BINWIDTH.
N.B.:
REBIN ID X Y EX EY is a convenient way to return in
one call abscissa, contents and errors for 1-Dim histogram.
In this case the errors in X are equal to 0.5*BINWIDTH.
* NTUPLE/HMERGE OUTFILE INFILES
OUTFILE C 'Output file name' D=' '
INFILES C 'Input file names' D=' ' Vararg
Merge HBOOK files containing histograms and/or ntuples. Ntuples are
merged and histograms with the same ID are added.
* KUIP/BUGREPORT
Send bug report to the PAW development team by email. The local editor is
automatically invoked with a template bug report. After the template has
been edited, full version information about the paw version and the
operating system is appended. The user is asked for a confirmation before
the report is send. ( The command is not yet implemented on all systems )
In PAW++ this command can be accessed from the `Help' menu of the
`Executive Window' (menu item `Mail Paw++ Develepers').
on Unix interpret command line by HOST_SHELL shell. For example if the current host shell is "ksh", the .kshrc script file is executed and all the aliases and variables defined in this file can be used. Before the command was passed to whatever shell was spawned by system().
SUBROUTINE HBAR2(ID2)
SUBROUTINE HRENAME(IDN,CHOLD,CHNEW)
CALL ISFACI(IFACI)
SUBROUTINE IGHTOR(RHI,RLI,RSI,R,G,B)
$CALL('fun.f(1.5)')
Rewrite of the KUMAC interpreter in C. The main difference is that EXEC inside a macro is treated like any other command.
ON ERROR CONTINUE
ON ERROR GOTO label
ON ERROR EXITM value
ON ERROR STOPM
Special thanks to Robert Franchisseur, Mike Kelsey, and Andrea Parri for beta-testing the new macro interpreter.
panel x.y command [label] [pixmap]
where:
x.y is the key position (column and row number),
command is the complete command (or list of commands) to
be executed when the corresponding button is pressed,
label (optional) is an alias-name for this command.
If specified, this alias-name is used for the button label
(when the appropriate ``View'' option is selected) instead of
the complete command (which is generally too long for a
``user-friendly'' button label).
pixmap (optional) has to be specified for graphical keys .
When ``pixmap'' is specified (for representing the key/button graphically with an icon) the graphical representation is displayed by default. It is possible at run time to vizualize the alphanumerical representation of the panel (with labels or commands) by selecting the appropriate entry in the menu_bar entry ``View'' at the top of the panel.
When the optional parameter ``label'' is missing but ``pixmap'' present, label should be replaced by a ``.''.
The application can define its own icons (pixmaps). This can be done by the application programmer in the command definition file (CDF) (following the directive ``>Icon_bitmaps'') or at run-time by the user himself with the command:
/MOTIF/ICON pixmap filename
where:
pixmap is the name given to the icon bitmap and used in the ``panel''
command described above,
filename is the name of the file where the icon bitmap datas are
stored.
To create a new icon bitmap (or pixmap) one can use the X11 standard bitmap editor ``bitmap''. E.g., to get a 20x20 pixel icon called m1 one can execute ``bitmap m1.bm 20x20''. The output file m1.bm containing ``#define m1_width 20 ...'' has to be referred in the command ``/MOTIF/ICON'' (e.g. /MOTIF/ICON m1 /user/.../.../m1.bm) or can be put in the command definition file (CDF).
It is possible to vizualize at run-time the panel in many ways by selecting the appropriate entry in the menu_bar menu ``View'' at the top of the panel:
MULTI_PANEL 'My Palette' '200x100+0+0'
will display a palette, or "multi_panel" widget, with the given title and geometry. After this command execution all panel definitions and executions will go into that widget. This can be done simply by executing the KUIP macro containing your panel definitions in the command line area, or, by selecting the "Add button" entry in the menu-bar entry "File" at the top of the palette.
The ``arrow buttons'' can be pressed either to reduce the panel to a label containing the panel title (arrow button is then turned left to right) or to display it (arrow button turned up to down).
To end a "multi-panel" setting one just have to issue the command:
MULTI_PANEL end
The next panel definitions and executions will be displayed as usual individual panels and will not go anymore into the palette of panels.
The following sequence of commands (which can be put inside a macro) can be used to set up a palette:
MULTI_PANEL
EXEC PANEL1.KUMAC
EXEC PANEL2.KUMAC
EXEC PANEL3.KUMAC
MULTI_PANEL end
where PANEL1.KUMAC, PANEL2.KUMAC, and PANEL3.KUMAC are KUIP macro files with ``usual'' panel setting and definition.
N.B. For users who are familiar with the ``UIMX'' User Interface Management System, the KUIP/Motif ``multi_panel'' is an emulation of the ``Palette'' widget which is built-in inside this program.
Other C callable routines have also been provided according to user request to facilitate, for instance, the implementation of a user-defined single selection list or file selection box (without having to be involved in direct Motif programming).
All information concerning this C callable interface can be found in the KUIP manual.
paw > call file.f ! compilation and interpretation with COMIS
paw > call file.f77 ! compilation of file.f with f77 and dynamic linking
paw > call file.c ! same with the C compiler
The possible ways to load shared objects:
Macro comp
appl comis quit
!file filename.f77
quit
in this case the COMIS compiler is invoked to process
filename.f COMIS fills a list of used variables for ntuple watched common
blocks, then COMIS creates a script file to invoke the native fortran compiler and loader to create a shared object, then the shared object is loaded.
Macro comp
appl comis quit
!file filename.c
quit
in this case COMIS creates a script file to invoke
the native C compiler and loader to create a shared object, then the shared
object is loaded.
Macro comp
appl comis quit
!file filename.sl
quit
in this case COMIS compiles a file 'filename.f' to
fill a list of used variables for ntuple watched common blocks, then the
shared object is loaded.
Macro comp
appl comis quit
!file filename.csl
quit
in this case the shared object (filename.sl) is
loaded The default script file looks something like this:
#! /bin/sh olddir=`pwd` cd CHPATH /bin/rm -f name.sl # for C compiler # cc -c .... name.c CHCC name.c #for f77 compiler # f77 -c .... name.f CHF77 name.f errno=$?' if [ $errno != 0 ] then exit $errno fi #for HPUX. ld -b -o name.sl name.o #for SUN. # ld -G -o name.sl name.o for solaris ld -o name.sl name.o #for SGI. ld -shared -o name.sl name.o # errno=$? if [ $errno != 0 ] then exit $errno fi /bin/chmod 555 name.sl /bin/rm -f name.o cd $olddir exit 0
This default script can be modified using COMIS directives:
!setopt 'string' path
The user can redefine CHPATH - the new directory where shared object and temporary files will be located( default is /tmp/ ),
!setopt 'string' f77
The user can redefine CHF77 - fortran compiler directive default is for HPUX
f77 -c +z +ppu -K -O
for SUN
f77 -c -pic
for SGI
f77 -c
!setopt 'string' cc
The user can redefine CHCC - C compiler directive default is for HPUX
cc -c +z -O
for SUN
cc -c -pic
for SGI
cc -cckr -c
PAW > H/PLOT 20 SURF,FB
plots the histogram 20 as a surface without the front box.
Zone 5 5
Do i=1,25
Set btyp [i]
Nul
Enddo
'Histogram_title ; XXX ; YYY ; ZZZ'
Where 'XXX', 'YYY' and 'ZZZ' are respectively the X, Y and Z axis titles. Example:
PAW > 1DHisto 10 'Histogram_title ; XXX ; YYY' 100 0 1
PAW > H/plot 10
PAW > 2DHisto 20 'Histogram_title ; XXX ; YYY ; ZZZ' 40 0 1 40 0 1
PAW > H/plot 20
Produce:
/ \
/ \
/ \
/ \
/ \
/ \
Y +--------------------------+ Z |\ /|
Y | | Z | \ / |
Y | | Z | \ / |
| | | \ / |
| | | \ / |
| | | \ / |
| | | | |
| | | | |
| | | | |
| | | | |
| | \ | / X
| | Y \ | / X
| | Y \ | / X
| | Y \ | /
| | \ | /
+--------------------------+ \|/
XXX
Histogram_title Histogram_title
In this model, one field can be blank. For example an histogram title like:
' ; XXX ; YYY'
will remove the Histogram_title but will leave the titles along the X and Y axis.
The command ATITLE has been updated according to this new scheme: it has an additional parameter: ZTIT.
* GRAPHICS/HPLOT/ATITLE [ XTIT YTIT ZTIT ]
XTIT C 'X Axis title' D=' '
YTIT C 'Y Axis title' D=' '
ZTIT C 'Z Axis title' D=' '
Draw axis titles on the axes of the present plot zone.
Thanks for all people who have contributed to the discussion, on the HEPLIB mailing list, about this new histogram title scheme.
This routine can be used to create the data structure to store errors for 2-D histograms (like HBARX for 1-Ds). By default, the errors are set to the sqrt(contents).
The HBOOK routine HPAKE (or the PAW command PUT/ERR) may be used to fill errors. It HBAR2 is not called before HPAKE, then HPAKE invokes HBAR2 automatically. When HBAR2 is called, routines HFILL or HF2 will accumulate the sum of square of weights. The errors for 2-Ds are used by the fit routines or the Histo/FIT command. The error bars are not drawn by the HPLOT routines.
This routine can be called from the PAW command line or from a COMIS routine
Subroutine to merge the NFILES HBOOK files with identical objects and directories into a new file FILOUT.
Same action as HMERGE, but the list of files to be merged and the output file are prompted interactively. The PAW command NTUPLE/HMERGE invokes this routine.
Duplicate definition of an Ntuple. A new Ntuple ID2 is created with the same definitions as ID1, but with 0 entries.
ID1 is the original id, ID2 the new
NEWBUF <0 use buffer size of ID1
NEWBUF =0 use current buffer size(10000 for RWN's)
NEWBUF >0 use NEWBUF as buffer size
CHTITL the new title, when blank use original,
CHOPT=' ' use original value (disk or memory resident)
CHOPT='M' to make the Ntuple memory resident
CHOPT='A' interactive mode: Sets addresses to /PAWCR4/
This routine can be called from the PAW command line or from a COMIS routine
Returns the variable definition as given in HBNAME for variable with index IVAR in N-tuple ID1. N-tuple must already be in memory.
To rename column CHOLD of ntuple ID into CHNEW. o New routine HCONVOL: CALL HCONVOL( D1,ID2,ID3,IERROR)
Action: Perform 1D convolution of ID2 with ID1 as kernel Result in ID3.
Setup: All histograms involved must exist before this call.
ID2 and ID3 histograms must be 1D.
ID1 can be 1-D or 2-D
Author: Per Steinar Iversen (PerSteinar.Iversen@fi.uib.no)
NB: This method scales badly for large histograms. The best general algorithm would be to unpack the histograms, add a suitable number of zeros, do the two FFTs, multiply the transforms, do yet another FFT and stuff the resulting histogram back into HBOOK. However, for small histograms, the naive method is probably faster, especially if recoded in terms of lower level calls. It will also work with 1D and 2D kernel-histograms that do not have matched coordinate systems - the FFT method implies equal binsize in X and Y for the kernel and the histogram to be folded; This simple method uses HBOOK to avoid this in one (X) dimension at least, corresponding to folding in a constant resolution term.
EXAMPLE of use
CALL HBOOK1(1,'Kernel 1 - 1D',100,-5.0,5.0,0.0)
CALL HBOOK2(2,'Kernel 2 - 2D',100,-5.0,5.0,100,0.0,100.0,0.0)
CALL HBPRO(2,0.0)
CALL HBOOK2(3,'Kernel 3 - 2D',100,-5.0,5.0,100,0.0,100.0,0.0)
CALL HBPRO(3,0.0)
CALL HBOOK1(4,'Function',100,0.0,100.0,0.0)
CALL HBOOK1(5,'Result 1',100,-10.0,110.0,0.0)
CALL HBOOK1(6,'Result 2',100,-10.0,110.0,0.0)
CALL HBOOK1(7,'Result 3',100,-10.0,110.0,0.0)
DO 10 I=1,100000
CALL RANNOR(A,B)
CALL HFILL(1,A,0.0,1.0)
CALL HFILL(1,B,0.0,1.0)
CALL RANNOR(A,B)
CALL HFILL(2,A,100.0*RNDM(I),1.0)
CALL HFILL(2,B,100.0*RNDM(I),1.0)
CALL RANNOR(A,B)
CALL HFILL(3,A,100.0*SQRT(RNDM(I+1)),1.0)
CALL HFILL(3,B,100.0*SQRT(RNDM(I+1)),1.0)
X = 30.0*(RNDM(I)-0.5)+50.0
CALL HFILL(4,X,0.0,1.0)
10 CONTINUE
CALL HCONVOL(1,4,5,IERROR)
CALL HCONVOL(2,4,6,IERROR)
CALL HCONVOL(3,4,7,IERROR)
A new option has been added in HLIMAP (Unix only). When the first parameter LIMIT=0, an existing shared memory is attached and become the current HBOOK directory. In this case HLIMIT must have been called before. The short example below illustrates how to use this option.
Program Example
Parameter (nwpaw=100000)
common/pawc/paw(nwpaw)
call hlimit(nwpaw)
*
call hlimap(0,'TEST')
*
1 read *,id
call hrin(id,9999,0)
call hprint(id)
call hdelet(id)
if(id.ne.0)go to 1
end
New option 'B' to be given with option 'E' to compute Binomial Errors in case of division.
CALL HOPERA(ID1,'/BE',ID2,ID3,1.,1.)
Compute errors for 2-Ds if option 'E' is given.
Options PROS or PROE are displayed with option SHOW
Can reset ntuples as well.
New option 'K' with option 'M'. When 'KM' is given, the parameters set in Application HMINUIT are not reset.
Profile histograms can be filled with weights.
The computation of the errors is based on a proposal by Stephane Coutu. Below is a copy of Coutu's proposal. ONLY options 'S' and 'I' are implemented now.
I realized that there is another case where this kind of trouble occurs: if a bin has N data points all with the same value Y (especially possible when dealing with integers), the spread in Y for that bin is zero, and the uncertainty assigned is also zero, and the bin is ignored in making subsequent fits. If SQRT(Y) was the correct error in the case above, then SQRT(Y)/SQRT(N) would be the correct error here. In fact, any bin with non-zero number of entries N but with zero spread should have an uncertainty SQRT(Y)/SQRT(N).
Now, is SQRT(Y)/SQRT(N) really the correct uncertainty? I believe that it is only in the case where the Y variable is some sort of counting statistics, following a Poisson distribution. This should probably be set as the default case. However, Y can be any variable from an original NTUPLE, not necessarily distributed "Poissonly", and perhaps extra options could be offered with the command: PROFILE id title ncx xmin xmax ymin ymax [ chopt ] to allow the user to choose how errors are calculated. We could have, for example:
CHOPT
' ' (Default) Errors are Spread/SQRT(N) for Spread.ne.0. ,
" " SQRT(Y)/SQRT(N) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
'S' Errors are Spread for Spread.ne.0. ,
" " SQRT(Y) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
'I' Errors are Spread/SQRT(N) for Spread.ne.0. ,
" " 1./SQRT(12.*N) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
The third case above corresponds to Integer Y values for which the uncertainty is +-0.5, with the assumption that the probability that Y takes any value between Y-0.5 and Y+0.5 is uniform (the same argument goes for Y uniformly distributed between Y and Y+1); this would be useful if Y is an ADC measurement, for example. Other, fancier options would be possible, at the cost of adding one more parameter to the PROFILE command. For example, if all Y variables are distributed according to some known Gaussian of standard deviation Sigma, then:
'G' Errors are Spread/SQRT(N) for Spread.ne.0. ,
" " Sigma/SQRT(N) for Spread.eq.0,N.gt.0 ,
" " 0. for N.eq.0
For example, this would be useful when all Y's are experimental quantities measured with the same instrument with precision Sigma.
Stephane Coutu
coutu@roo.physics.lsa.umich.edu
Routine HLNEXT supports RLOGIN directories (This means now all possible combinations of directories).
Mods in HMINUT to compute an equivalent chisquare in case of a log-likelihood fit.
2-D vectors may now be specified. Also apply to PAW command Vector/Fit.
HMCMLL already uses the new HMCLNL, old HMCLNL renamed as HMCLNO, but will be deleted in a couple of months.
HMCINI and the new and old versions of HMCLNL both contain a banner announcing the change, as it's not backwards compatible.
The weight histograms may be used for more than one of the fits if necessary - a check is made to make sure that they are not normalised more than once.
Version 3.00:
Remove check on RZ version in RZFILE Mods in RZOPEN and RZIODO to support PIAF file striping New routines RZSTRIP and RZSTRIR added Removal of 64K limitations (Sunanda Banerjee) 20/06/94