2
C#### YOU WILL NOT NEED WHAT FOLLOWS IN YOUR OWN ANALYSIS
4
C************************************************************************
6
C HISTOGRAMMING PACKAGE
9
C************************************************************************
11
C**********************************************************************
12
C SIMPLE HISTOGRAMMING PACKAGE -- SIMPLIFIED VERSION OF HBOOK
13
C BY Michelangelo Mangano NOVEMBER 1988
14
C LAST REVISED NOVEMBER 9, 1988
15
C LAST REVISED JUNE 12, 1989 (ADD SCATTER PLOTS)
16
C LAST REVISED oct 1990 (Add multi-plots on one page, routines MULTITOP,
18
C**********************************************************************
20
C Fills up to 100 histograms with up to 100 bins.
21
C Gives a data file (to be specified in the calling program by assigning
22
C a file name to unit 98) and a topdrawer file (to be specified in the
23
C calling program by assigning a file name to unit 99).
26
C Call once INIHIST; this just resets a few counters and logicals
27
C Call MBOOK(N,'TITLE',DEL,XMIN,XMAX) for each histogram to be booked.
28
C N (an integer) is the label of the histogram;
29
C 'TITLE' is the name of the histogram (no more then 100 characters);
30
C DEL (real*8) is the bin size;
31
C XMIN (real*8) is the lower limit of the first bin;
32
C XMAX (real*8)is the upper limit of the last bin
34
C call mbook(2,'pt distribution',1.,10,70)
35
C This call initializes histogram number 2, called 'pt distribution';
36
C The bin size will be 1. (possibly GeV, if that's what you want), the
37
C first bin being 10<x<11. and the last one being 69.<x<70
40
C When it's time, call MFILL(N,X,Y); this will add Y (real*8) to the bin
41
C in which X (real*8) happens to be, within histogram N.
44
C At the end of the day you may want to sum, divide, cancel, etc.etc.
45
C various histograms (bin by bin). Then you call MOPERA(I,'O',J,K,X,Y).
46
C The 1-character string O can take the following values:
47
C + : sums X*(hist I) with Y*(hist J) and puts the result in hist K;
48
C - : subtracts X*(hist I) with Y*(hist J) and puts the result in hist K;
49
C * : multiplies X*(hist I) with Y*(hist J) and puts the result in hist K;
50
C / : divides X*(hist I) with Y*(hist J) and puts the result in hist K;
51
C F : multiplies hist I by the factor X, and puts the result in hist K;
52
C R : takes the square root of hist I, and puts the result in hist K;if
53
C the value at a given bin is less than or equal to 0, puts 0 in K
54
C S : takes the square of hist I, and puts the result in hist K;
55
C L : takes the log_10 of hist I, and puts the result in hist K; if the
56
C value at a given bin is less than or equal to 0, puts 0 in K
57
C M : statistical analysis; if I contains the weights (let's say WGT),
58
C J contains variable times weight (F*WGT) and K contains the
59
C variable squared times the weight (F**2*WGT), then, after using 'M',
60
C J will contain the average value of the variable <F> and K will
61
C contain the sigma of the average: sigma=sqrt(<F**2>-<F>**2).
62
C If WGT=1. for all the entries, then it is enough to put I=J, and
63
C it is not necessary to book a hist with the weights.
64
C V : estimates errors for vegas evaluation of differential distributions.
65
C Fill I with the values of
66
C the functions do integrate times the Vegas weight (fun*wgt); fill
67
C J with fun**2*wgt; then K will contain an estimate of the error
68
C of the integration. Putting X=1/(#of iterations) performs the
69
C avegare over the iterations, and gives the right normalization to
70
C the differential distribution, I, and to the errors, K. J stays the same.
73
C Now we can finalize our histograms; MFINAL(N) will calculate the integral
74
C of the histogram N, the mean value of the X variable and its RMS.
75
C If we now want to renormalize the hist's, we can call MNORM(N,X), which
76
C will normalize the integral to X -- CAUTION: do not call MNORM before
77
C MFINAL, it will blow up.
80
C To get a .dat file containing the values of the histograms, together with
81
C some information (like integral, mean values, etc.etc.) call MPRINT(N),
82
C for each hist N that you want in the .dat file. Before the call to MPRINT
83
C you want to open unit 98 and give it a name:
84
C OPEN(UNIT=98,NAME='NAME.DAT',STATUS='NEW')
85
C If you want a topdrawer file with a plot of the hist values, call
86
C MTOP(N,M,'X','Y','SCALE'). The points of the plot will be taken from histogram
87
C N, the error bars from histogram M. 'SCALE', character*(*), determines
88
C the scale for y, logarithmic or linear (SCALE=LOG,LIN).
89
C If you do not want error bars, keep
90
C a histogram of zeros, or just call a hist that had not been booked.
91
C X will appear as a 'bottom title', and Y will appear as a 'left title'.
92
C The top title is by default the name of the histogram itself.
93
C A little box below the plot will contain some information on the plot
94
C itself. Before calling MTOP,
95
C OPEN(UNIT=99,NAME='NAME.TOP',STATUS='NEW')
96
C--------------------------------------------------------------------------
98
C COMMON/HISTO/ Histogram N
100
C BOOK(N), Three-letter character-string: 'NO' if histogram was not
101
C Booked, 'RES' if it was resetted, 'YES' otherwise
102
C TITLE(N), Title of the histogram
104
C HMIN(N), Min value of x range
105
C HMAX(N), Max value of x range
107
C NBIN(N), Total number of bins
108
C USCORE(N), Total integral of underscores with x < HMIN(N)
109
C OSCORE(N), Total integral of onderscores with x > HMAX(N)
110
C IUSCORE(N), Number of entries with x < HMIN(N)
111
C IOSCORE(N), Number of entries with x > HMAX(N)
112
C IENT(N), Total number of entries within x range HMIN(N)<x<HMAX(N)
113
C HINT(N), Integral of the histogram within HMIN(N)<x<HMAX(N)
114
C HAVG(N), Average value of x, weighted over the x range of the histo
115
C HSIG(N), Quadratic dispersion of x around the average
116
C HIST(N,L), Value of bin L-th
117
C XHIS(N,L), Central x value of bin L-th
118
C IHIS(N,L), Number of entries within bin L-th
119
C NHIST Total number of booked histograms
125
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
135
SUBROUTINE MBOOK(N,TIT,DEL,XMIN,XMAX)
137
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
143
IF(BOOK(N)(1:1).EQ.'Y') THEN
145
WRITE(*,*) 'Histogram',N,TITLE(N),' already in use. '
146
WRITE(*,*) 'superseded by ',TIT
151
NBIN(N) = INT((XMAX-XMIN)/(DEL*0.999999d0))
152
IF(NBIN(N).GT.100) THEN
153
WRITE(*,*) 'TOO MANY BINS (',NBIN(N),') REQUIRED IN HIST ',N
154
WRITE(*,*) 'RE-ENTER BIN SIZE DELTA (OLD BIN = ',DEL,' ):'
159
HMAX(N) = NBIN(N)*DEL+XMIN
160
IF(ABS(HMAX(N)-XMAX).GT.1.D-3*DEL) THEN
163
#'Histogram ', TIT, ' Change of upper limit:',xmax,'-->',HMAX(N)
174
XHIS(N,I)=HMIN(N)+HDEL(N)*(DFLOAT(I)-0.5)
180
SUBROUTINE MFILL(N,X,Y)
182
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
185
IF(X.LT.HMIN(N)) THEN
186
USCORE(N) = USCORE(N) + Y
187
IUSCORE(N) = IUSCORE(N) + 1
188
ELSEIF(X.GT.HMAX(N)) THEN
189
OSCORE(N) = OSCORE(N) + Y
190
IOSCORE(N) = IOSCORE(N) + 1
192
XI=((X-HMIN(N))/HDEL(N))+1
195
IHIS(N,I)=IHIS(N,I)+1
196
C SF, 10/3/2010 -- avoids numerical problems
198
IF(ABS(SUM).LT.1.d-19)THEN
207
SUBROUTINE MINTEG(NIN,NOUT,IDIR,IPOW)
208
C If IPOW=1 performs the integral of the distribution contained in histogram
209
C NIN up to the value specified by the abscissa (if IDIR=1) or from this
210
C value on (if IDIR=-1). The resulting integral distribution is put into
211
C NOUT, which is automatically booked if NOUT.ne.NIN . Choosing IPOW=2
212
C the routine will return the square root of the integral of the squares,
213
C as is required, for example, for the propagation of the mean quadratic error
214
C of a given distribution. Overscores and underscores are included.
216
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
217
INTEGER NIN,NOUT,IDIR,IPOW,I,L,M
221
DATA C/' INTEG BELOW X',' INTEG ABOVE X'/
222
IF(BOOK(NIN)(1:1).NE.'Y') RETURN
225
IF(NOUT.NE.NIN.AND.BOOK(NOUT)(1:1).NE.'Y') THEN
226
CALL MBOOK(NOUT,TITLE(NIN)//C(I),
227
& HDEL(NIN),HMIN(NIN),HMAX(NIN))
230
HIST(NOUT,1) = SUMPOW(HIST(NIN,1),USCORE(NIN),IPOW)
231
IHIS(NOUT,1) = IHIS(NIN,1) + IUSCORE(NIN)
232
XHIS(NOUT,1) = XHIS(NIN,1) + HDEL(NIN)/2
234
HIST(NOUT,L) = SUMPOW(HIST(NIN,L),HIST(NOUT,L-1),IPOW)
235
IHIS(NOUT,L) = IHIS(NIN,L) + IHIS(NOUT,L-1)
236
XHIS(NOUT,L) = XHIS(NIN,L) + HDEL(NIN)/2
238
OSCORE(NOUT) = SUMPOW(OSCORE(NIN),HIST(NIN,M),IPOW)
239
IOSCORE(NOUT) = IOSCORE(NIN) + IHIS(NIN,M)
240
ELSEIF(IDIR.EQ.-1) THEN
241
HIST(NOUT,M) = SUMPOW(HIST(NIN,M),OSCORE(NIN),IPOW)
242
IHIS(NOUT,M) = IHIS(NIN,M) + IOSCORE(NIN)
243
XHIS(NOUT,M) = XHIS(NIN,M) - HDEL(NIN)/2
245
HIST(NOUT,L) = SUMPOW(HIST(NIN,L),HIST(NOUT,L+1),IPOW)
246
IHIS(NOUT,L) = IHIS(NIN,L) + IHIS(NOUT,L+1)
247
XHIS(NOUT,L) = XHIS(NIN,L) - HDEL(NIN)/2
249
USCORE(NOUT) = SUMPOW(USCORE(NIN),HIST(NIN,1),IPOW)
250
IUSCORE(NOUT) = IUSCORE(NIN)+IHIS(NIN,1)
252
CALL MBWARN('MINTEG')
253
WRITE(*,*) 'Wrong idir in minteg: OPERATION NOT PERFORMED'
258
FUNCTION SUMPOW(X,Y,IPOW)
264
ELSEIF(IPOW.EQ.2) THEN
265
SUMPOW = SQRT(X**2+Y**2)
266
ELSEIF(IPOW.EQ.0) THEN
267
CALL MBWARN('SUMPOW')
268
WRITE(*,*)'Error: IPOW=0 not allowed in SUMPOW'
270
SUMPOW = (X**IPOW+Y**IPOW)**(1./IPOW)
274
SUBROUTINE MOPERA(I,OPER,J,K,X,Y)
276
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
280
REAL*8 XAVG,XNORM,XSQAVG,XSUM,XSUMSQ,XXX
282
IF(NBIN(I).NE.NBIN(J).AND.(OPER.EQ.'+'.OR.OPER.EQ.'-'.OR.OPER.EQ.
283
& '*'.OR.OPER.EQ.'/'.OR.OPER.EQ.'M'.OR.OPER.EQ.'A')) THEN
284
CALL MBWARN('MOPERA')
286
20 FORMAT(' ****** INCOMPATIBLE OPERATION HIST ',I2,' &',I2,
290
IF(BOOK(I)(1:1).NE.'Y'.OR.BOOK(J)(1:1).NE.'Y') RETURN
291
IF(BOOK(K)(1:1).NE.'Y')
292
& CALL MBOOK(K,TITLE(I),HDEL(I),HMIN(I),HMAX(I))
294
c If I contains the accumulated weights, J the accumulated squares of the
295
c weights and IHIS(J,1) the number of accumulated entries, 'E' will add
296
c the average value of I to K and will put in J the quadratic dispersion.
297
IF(IHIS(J,1).NE.0) THEN
298
XXX = 1.D0/DFLOAT(IHIS(J,1))
305
HIST(K,L)=HIST(K,L) + XXX*HIST(I,L)
306
HIST(J,L)=XXX*SQRT(ABS(XSUMSQ-XSUM**2*XXX))
308
IENT(K)=IENT(K)+IENT(I)
311
USCORE(K) = USCORE(K)+XXX*XSUM
312
USCORE(J) = XXX*SQRT(ABS(XSUMSQ-XSUM**2*XXX))
315
OSCORE(K) = OSCORE(K)+XXX*XSUM
316
OSCORE(J) = XXX*SQRT(ABS(XSUMSQ-XSUM**2*XXX))
317
ELSEIF(OPER.EQ.'Q') THEN
319
HIST(K,L) = SQRT(HIST(J,L)**2+HIST(I,L)**2)
321
USCORE(K) = SQRT(USCORE(J)**2+USCORE(I)**2)
322
OSCORE(K) = SQRT(OSCORE(J)**2+OSCORE(I)**2)
323
ELSEIF(OPER.EQ.'A') THEN
325
HIST(J,L) = HIST(J,L) + HIST(I,L)
326
IHIS(J,L) = IHIS(J,L) + IHIS(I,L)
327
HIST(K,L) = HIST(K,L) + HIST(I,L)**2
328
IHIS(K,L) = IHIS(K,L) + 1
333
IENT(J) = IENT(J)+IENT(I)
334
IUSCORE(J) = IUSCORE(J) + IUSCORE(I)
335
USCORE(J) = USCORE(J) + USCORE(I)
336
IOSCORE(J) = IOSCORE(J) + IOSCORE(I)
337
OSCORE(J) = OSCORE(J) + OSCORE(I)
338
IUSCORE(K) = IUSCORE(K) + 1
339
USCORE(K) = USCORE(K) + USCORE(I)**2
340
IOSCORE(K) = IOSCORE(K) + 1
341
OSCORE(K) = OSCORE(K) + OSCORE(I)**2
350
HIST(K,L)=X*HIST(I,L) + Y*HIST(J,L)
351
ELSEIF(OPER.EQ.'-') THEN
352
HIST(K,L)=X*HIST(I,L) - Y*HIST(J,L)
353
ELSEIF(OPER.EQ.'*') THEN
354
HIST(K,L)=X*HIST(I,L) * Y*HIST(J,L)
355
ELSEIF(OPER.EQ.'/') THEN
356
IF(Y.EQ.0..OR.HIST(J,L).EQ.0.) THEN
359
HIST(K,L)=X*HIST(I,L) / (Y*HIST(J,L))
361
ELSEIF(OPER.EQ.'F') THEN
362
HIST(K,L)=X*HIST(I,L)
363
ELSEIF(OPER.EQ.'R') THEN
364
IF(HIST(I,L).GT.0.) THEN
365
HIST(K,L)=X*SQRT(HIST(I,L))
369
ELSEIF(OPER.EQ.'S') THEN
370
HIST(K,L)=X*HIST(I,L)**2
371
ELSEIF(OPER.EQ.'L') THEN
372
IF(HIST(I,L).EQ.0..OR.J.EQ.0.) THEN
375
HIST(K,L)=X*LOG10(Y*HIST(I,L))
377
ELSEIF(OPER.EQ.'M') THEN
378
IF(I.NE.J) XNORM=HIST(I,L)
379
IF(I.EQ.J) XNORM=DFLOAT(IHIS(J,L))
383
& SQRT(ABS(-XAVG**2+HIST(K,L)/XNORM)/DFLOAT(IHIS(I,L)))
389
ELSEIF(OPER.EQ.'V') THEN
392
XNORM=DFLOAT(IHIS(I,L))*X
394
HIST(K,L)=SQRT(ABS(XSQAVG-XAVG**2)/XNORM)
400
CALL MBWARN('MOPERA')
402
5 FORMAT(' ****** OPERATION ="',A1,'" UNKNOWN ********'/)
407
USCORE(K)=X*USCORE(I) + Y*USCORE(J)
408
OSCORE(K)=X*OSCORE(I) + Y*OSCORE(J)
409
ELSEIF(OPER.EQ.'-') THEN
410
USCORE(K)=X*USCORE(I) - Y*USCORE(J)
411
OSCORE(K)=X*OSCORE(I) - Y*OSCORE(J)
412
ELSEIF(OPER.EQ.'*') THEN
413
USCORE(K)=X*USCORE(I) * Y*USCORE(J)
414
OSCORE(K)=X*OSCORE(I) * Y*OSCORE(J)
415
ELSEIF(OPER.EQ.'/') THEN
416
IF(Y.EQ.0..OR.USCORE(J).EQ.0.) THEN
419
USCORE(K)=X*USCORE(I) / (Y*USCORE(J))
421
IF(Y.EQ.0..OR.OSCORE(J).EQ.0.) THEN
424
OSCORE(K)=X*OSCORE(I) / (Y*OSCORE(J))
426
ELSEIF(OPER.EQ.'F') THEN
427
USCORE(K)=X*USCORE(I)
428
OSCORE(K)=X*OSCORE(I)
429
ELSEIF(OPER.EQ.'R') THEN
430
IF(USCORE(I).GT.0.) THEN
431
USCORE(K)=X*SQRT(USCORE(I))
435
IF(OSCORE(I).GT.0.) THEN
436
OSCORE(K)=X*SQRT(OSCORE(I))
440
ELSEIF(OPER.EQ.'S') THEN
441
USCORE(K)=X*USCORE(I)**2
442
OSCORE(K)=X*OSCORE(I)**2
443
ELSEIF(OPER.EQ.'L') THEN
444
IF(USCORE(I).EQ.0..OR.J.EQ.0.) THEN
447
USCORE(K)=X*LOG10(Y*USCORE(I))
449
IF(OSCORE(I).EQ.0..OR.J.EQ.0.) THEN
452
OSCORE(K)=X*LOG10(Y*OSCORE(I))
459
SUBROUTINE MCOPY(NIN,NOUT)
461
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
464
C If the histogram is not booked or resetted, do not copy it
465
IF(BOOK(NIN)(1:1).NE.'Y'.AND.BOOK(NIN)(1:1).NE.'R') RETURN
466
C If the histogram was resetted without having being booked, do not copy it.
467
C This may happen if MFINAL is called prior to MCOPY
468
IF(BOOK(NIN)(1:1).EQ.'R'.AND.HMIN(NIN).EQ.HMAX(NIN)) RETURN
469
IF(BOOK(NOUT)(1:1).NE.'Y'.AND.BOOK(NOUT)(1:1).NE.'R')
470
& CALL MBOOK(NOUT,TITLE(NIN)(2:),HDEL(NIN),HMIN(NIN),HMAX(NIN))
471
IF(BOOK(NIN).EQ.'RES')BOOK(NOUT)='RES'
477
IUSCORE(NOUT)=IUSCORE(NIN)
478
IOSCORE(NOUT)=IOSCORE(NIN)
479
USCORE(NOUT)=USCORE(NIN)
480
OSCORE(NOUT)=OSCORE(NIN)
485
HIST(NOUT,I)=HIST(NIN,I)
486
XHIS(NOUT,I)=XHIS(NIN,I)
487
IHIS(NOUT,I)=IHIS(NIN,I)
493
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
507
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
509
REAL*8 HISTL(100),XHISL(100),HDELL,HMINL,HMAXL,USCOREL,OSCOREL,
511
INTEGER IHISL(100),NBINL,IUSCOREL,IOSCOREL,IENTL,I
514
SAVE HISTL,XHISL,HDELL,HMINL,HMAXL,USCOREL,OSCOREL,
515
&NBINL,IHISL,IUSCOREL,IOSCOREL,IENTL,HAVGL,HINTL,HSIGL,
559
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
564
SUBROUTINE PUTTAG(J,NAME)
566
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
568
c Per marcare un istogramma
569
CHARACTER * (*) NAME, TAG
578
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
580
REAL*8 AVG,XIN,SIG,IF,X
592
IF(XIN.EQ.0) GO TO 10
595
SIG=HIST(N,J)*(XHIS(N,J)-AVG)**2+SIG
597
SIG=SQRT(ABS(SIG/XIN))
602
IF(IF.EQ.0) BOOK(N)='RES'
605
SUBROUTINE MFINAL3(N)
606
C Identical to MFINAL, except for the fact that if the histogram is void
607
C but was previously booked, BOOK is left as it is and not set to RES
609
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
611
REAL*8 AVG,XIN,SIG,IF,X
612
IF(BOOK(N)(1:1).NE.'Y') RETURN
617
IF(HMIN(N).NE.HMAX(N))IF=1
624
IF(XIN.EQ.0) GO TO 10
627
SIG=HIST(N,J)*(XHIS(N,J)-AVG)**2+SIG
629
SIG=SQRT(ABS(SIG/XIN))
634
IF(IF.EQ.0) BOOK(N)='RES'
637
SUBROUTINE MNORM(N,X)
639
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
642
CHARACTER CTIME*10,CDATE*10
643
IF(BOOK(N).NE.'YES')RETURN
644
IF(HINT(N).EQ.0.) THEN
645
c CALL MBKWRN('MNORM')
646
WRITE(*,*)' INTEGRAL HIST ',N,' IS ZERO: CANNOT RENORMALIZE'
652
1 HIST(N,I)=HIST(N,I)*Y
654
OSCORE(N)=OSCORE(N)*Y
655
USCORE(N)=USCORE(N)*Y
660
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
662
CHARACTER CTIME*10,CDATE*10
663
INTEGER INI,IMON,IDAY,IYEAR
666
CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
669
IF(BOOK(N)(:1).NE.'Y') RETURN
670
WRITE(98,7) N,IYEAR,IMON,IDAY,CTIME(1:5)
673
IF(HIST(N,J).EQ.0.) GO TO 1
674
WRITE(98,'(3X,F10.4,2X,E15.4)')
675
& XHIS(N,J),HIST(N,J)
677
WRITE(98,15) HAVG(N),HSIG(N),HINT(N)
678
WRITE(98,20) IENT(N),IUSCORE(N),IOSCORE(N)
679
7 FORMAT(4X,'HIST = ',I3,' 19',I2,'-',I2,'-',I2,1X,A5/)
681
15 FORMAT(/' AVG =',E10.3,4X,' RMS =',E10.3,' INTEGRAL =',E10.3,/)
682
20 FORMAT(' ENTRIES=',I5,4X,'UNDERSCORE=',I5,4x,'OVERSCORE=',I5,//)
685
SUBROUTINE MTOP(N,M,BTIT,LTIT,SCALE)
687
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
689
CHARACTER CTIME*10,CDATE*10
690
INTEGER INI,IMON,IDAY,IYEAR
691
CHARACTER*(*) LTIT,BTIT,SCALE
694
CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
697
IF(BOOK(N)(:1).NE.'Y') RETURN
698
WRITE(99,100) TITLE(N),BTIT,LTIT,SCALE,HMIN(N),HMAX(N)
700
&' SET WINDOW Y 2.5 TO 9.'/,1X,
701
&' SET WINDOW X 2.5 TO 10.'/,1X,
702
&' SET FONT DUPLEX '/1X,
703
&' SET SYMBOL 5O SIZE 1.8'/,1X,
704
&' TITLE TOP ','"',A,'"',/1X,
705
&' TITLE BOTTOM ','"',A,'"',/1X,
706
&' TITLE LEFT ','"',A,'"',/1X,
707
&' SET SCALE Y ',A,/1X,
708
&' (SET TICKS TOP OFF) '/1x,
709
&' SET LIMITS X ',F10.5,' ',F10.5,/1X,
710
&' SET ORDER X Y DY ')
712
WRITE(99,'(3X,F10.4,2(2X,E15.4))')
713
& XHIS(N,J),HIST(N,J),HIST(M,J)
717
WRITE(99,300) HINT(N),HAVG(N),HSIG(N),IENT(N),IUSCORE(N)
718
& ,IOSCORE(N),USCORE(N),OSCORE(N),IMON,IDAY,IYEAR,CTIME(1:5)
720
&' BOX 6.25 0.9 SIZE 7.5 1.2'/,1X,
721
&' SET WINDOW Y 0. TO 2.'/,1X,
722
&' SET TITLE SIZE -1.5'/1X,
723
&' SET FONT DUPLEX '/1X,
724
&' TITLE 2.8 1.2 "INT =',1PE10.3,' AVG =',1PE10.3,
725
& ' RMS =',1PE10.3,'"',/1X,
726
&' TITLE 2.8 0.9 "Entries =',I8,2x,'Undersc =',I6,2X
727
& ,'Oversc =',I6,'"',/1X,
728
&' TITLE 2.8 0.6 "Intgr ufloat=',1PE10.3,' Intgr ofloat=',
730
&' TITLE 7.5 0.6 "',I2,'-',I2,', 19',I2,2X,A5,'"',/1X,
731
&' SET TITLE SIZE -2')
733
400 FORMAT(' NEW PLOT')
736
SUBROUTINE MNRTOP(N,M,BTIT,LTIT,SCALE)
738
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
740
CHARACTER CTIME*10,CDATE*10
741
INTEGER INI,IMON,IDAY,IYEAR
742
CHARACTER*(*) LTIT,BTIT,SCALE
745
CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
749
WRITE(99,100) TITLE(N)(1:20),BTIT,LTIT,SCALE,HMIN(N),HMAX(N)
751
&' SET WINDOW Y 2.5 TO 7.'/,1X,
752
&' SET WINDOW X 2.5 TO 10.'/,1X,
753
&' SET FONT DUPLEX '/1X,
754
&' SET SYMBOL 5O SIZE 1.8'/,1X,
755
&' TITLE TOP ','"',A,'"',/1X,
756
&' TITLE BOTTOM ','"',A,'"',/1X,
757
&' TITLE LEFT ','"',A,'"',/1X,
758
&' SET SCALE Y ',A,/1X,
759
&' (SET TICKS TOP OFF) '/1x,
760
&' SET LIMITS X ',F10.5,' ',F10.5,/1X,
761
&' SET ORDER X Y DY ')
763
IF(HIST(N,J).EQ.0) GO TO 1
764
WRITE(99,'(3X,F10.4,2(2X,E10.3))')
765
& XHIS(N,J),HIST(N,J),HIST(M,J)
769
WRITE(99,300) HINT(N),HAVG(N),HSIG(N),IENT(N),IUSCORE(N)
770
& ,IOSCORE(N),USCORE(N),OSCORE(N),IMON,IDAY,IYEAR,CTIME(1:5)
772
&' BOX 6.25 0.9 SIZE 7.5 1.2'/,1X,
773
&' SET WINDOW Y 0. TO 2.'/,1X,
774
&' SET TITLE SIZE -1.5'/1X,
775
&' SET FONT DUPLEX '/1X,
776
&' TITLE 2.8 1.2 "INT =',1PE10.3,' AVG =',1PE10.3,
777
& ' RMS =',1PE10.3,'"',/1X,
778
&' TITLE 2.8 0.9 "Entries =',I8,2x,'Undersc =',I6,2X
779
& ,'Oversc =',I6,'"',/1X,
780
&' TITLE 2.8 0.6 "Intgr ufloat=',1PE10.3,' Intgr ofloat=',
782
&' TITLE 7.5 0.6 "',I2,'-',I2,', 19',I2,2X,A5,'"',/1X,
783
&' SET TITLE SIZE -2')
785
400 FORMAT(' NEW PLOT')
789
SUBROUTINE MBOOK2(N,TIT,DELX,XMIN,XMAX,DELY,YMIN,YMAX)
791
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
794
REAL*8 DELX,XMIN,XMAX,DELY,YMIN,YMAX
797
IF(BOOK2(N)(1:1).EQ.'Y') THEN
798
CALL MBWARN('MBOOK2')
799
WRITE(*,*) 'Histogram',N,TITLE2(N),' already in use. '
800
WRITE(*,*) 'superseded by ',TIT
804
C-- setup x-boundaries
808
NBIN2(1,N)=INT((XMAX-XMIN)/(DELX*0.9999))
809
C-- setup y-boundaries
813
NBIN2(2,N)=INT((YMAX-YMIN)/(DELY*0.9999))
820
XHIS2(N,I)=HMIN2(1,N)+HDEL2(1,N)*(DFLOAT(I)-0.5)
823
YHIS2(N,I)=HMIN2(2,N)+HDEL2(2,N)*(DFLOAT(I)-0.5)
831
SUBROUTINE MFILL2(N,X,Y,WGT)
833
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
838
C Modified by SF on 10/8/2005; now similar to mfill()
839
IF( X.LT.HMIN2(1,N) .OR. X.GE.HMAX2(1,N) .OR.
840
# Y.LT.HMIN2(2,N) .OR. Y.GE.HMAX2(2,N) )THEN
841
IOSCORE2(N)=IOSCORE2(N)+1
843
XI=((X-HMIN2(1,N))/HDEL2(1,N))+1.
844
YI=((Y-HMIN2(2,N))/HDEL2(2,N))+1.
847
IF( I.GT.0.AND.I.LE.NBIN2(1,N) .AND.
848
# J.GT.0.AND.J.LE.NBIN2(2,N) )THEN
850
IHIS2(N,I,J)=IHIS2(N,I,J)+1
851
HIST2(N,I,J)=HIST2(N,I,J)+WGT
853
IOSCORE2(N)=IOSCORE2(N)+1
859
SUBROUTINE MOPERA2(I,OPER,J,K,X,Y)
861
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
866
IF(NBIN2(1,I).NE.NBIN2(1,J).OR.NBIN2(2,I).NE.NBIN2(2,J).
867
&AND.(OPER.EQ.'+'.OR.OPER.EQ.'-'.OR.OPER.EQ.
868
&'*'.OR.OPER.EQ.'/'.OR.OPER.EQ.'M')) GO TO 10
872
HIST2(K,L1,L2)=X*HIST2(I,L1,L2) + Y*HIST2(J,L1,L2)
873
ELSEIF(OPER.EQ.'-') THEN
874
HIST2(K,L1,L2)=X*HIST2(I,L1,L2) - Y*HIST2(J,L1,L2)
875
ELSEIF(OPER.EQ.'*') THEN
876
HIST2(K,L1,L2)=X*HIST2(I,L1,L2) * Y*HIST2(J,L1,L2)
877
ELSEIF(OPER.EQ.'/') THEN
878
IF(Y.EQ.0..OR.HIST2(J,L1,L2).EQ.0.) THEN
881
HIST2(K,L1,L2)=X*HIST2(I,L1,L2) / (Y*HIST2(J,L1,L2))
883
ELSEIF(OPER.EQ.'F') THEN
884
HIST2(K,L1,L2)=X*HIST2(I,L1,L2)
885
ELSEIF(OPER.EQ.'R') THEN
886
IF(HIST2(I,L1,L2).GT.0.) THEN
887
HIST2(K,L1,L2)=X*SQRT(HIST2(I,L1,L2))
891
ELSEIF(OPER.EQ.'S') THEN
892
HIST2(K,L1,L2)=X*HIST2(I,L1,L2)**2
893
ELSEIF(OPER.EQ.'L') THEN
894
IF(HIST2(I,L1,L2).EQ.0..OR.J.EQ.0.) THEN
897
HIST2(K,L1,L2)=X*LOG10(Y*HIST2(I,L1,L2))
901
5 FORMAT(' ****** OPERATION ="',A1,'" UNKNOWN ********'/)
908
20 FORMAT(' ****** INCOMPATIBLE OPERATION HIST2 ',I2,' &',I2,
913
SUBROUTINE MCOPY2(NIN,NOUT)
915
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
918
C If the histogram is not booked or resetted, do not copy it
919
IF(BOOK2(NIN)(1:1).NE.'Y'.AND.BOOK2(NIN)(1:1).NE.'R') RETURN
920
C If the histogram was resetted without having being booked, do not copy it.
921
C This may happen if MFINAL2 is called prior to MCOPY2
922
IF(BOOK2(NIN)(1:1).EQ.'R'.AND.
923
& HMIN2(1,NIN).EQ.HMAX2(1,NIN)) RETURN
924
IF(BOOK2(NOUT)(1:1).NE.'Y'.AND.BOOK2(NOUT)(1:1).NE.'R')
925
& CALL MBOOK2(NOUT,TITLE2(NIN),
926
& HDEL2(1,NIN),HMIN2(1,NIN),HMAX2(1,NIN),
927
& HDEL2(2,NIN),HMIN2(2,NIN),HMAX2(2,NIN))
928
IF(BOOK2(NIN).EQ.'RES')BOOK2(NOUT)='RES'
929
HDEL2(1,NOUT)=HDEL2(1,NIN)
930
HMIN2(1,NOUT)=HMIN2(1,NIN)
931
HMAX2(1,NOUT)=HMAX2(1,NIN)
932
NBIN2(1,NOUT)=NBIN2(1,NIN)
933
HDEL2(2,NOUT)=HDEL2(2,NIN)
934
HMIN2(2,NOUT)=HMIN2(2,NIN)
935
HMAX2(2,NOUT)=HMAX2(2,NIN)
936
NBIN2(2,NOUT)=NBIN2(2,NIN)
937
IENT2(NOUT)=IENT2(NIN)
938
IOSCORE2(NOUT)=IOSCORE2(NIN)
939
HAVG2(1,NOUT)=HAVG2(1,NIN)
940
HAVG2(2,NOUT)=HAVG2(2,NIN)
941
HINT2(NOUT)=HINT2(NIN)
942
HSIG2(1,NOUT)=HSIG2(1,NIN)
943
HSIG2(2,NOUT)=HSIG2(2,NIN)
945
XHIS2(NOUT,I)=XHIS2(NIN,I)
946
XPROJ(NOUT,I)=XPROJ(NIN,I)
949
YHIS2(NOUT,I)=YHIS2(NIN,I)
950
YPROJ(NOUT,I)=YPROJ(NIN,I)
954
HIST2(NOUT,I,J)=HIST2(NIN,I,J)
955
IHIS2(NOUT,I,J)=IHIS2(NIN,I,J)
961
SUBROUTINE MFINAL2(N)
963
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
966
IF(BOOK2(N)(:1).NE.'Y') RETURN
968
C-- projection on the x-axis
971
1 XPROJ(N,I)=XPROJ(N,I)+HIST2(N,I,J)
973
IF(XIN.EQ.0.) GO TO 10
974
C-- projection on the y-axis
977
3 YPROJ(N,J)=YPROJ(N,J)+HIST2(N,I,J)
983
SUBROUTINE MFINAL32(N)
984
C Identical to MFINAL3, except for the fact that if the histogram is void
985
C but was previously booked, BOOK is left as it is and not set to RES,
986
C and for the fact that XPROJ and YPROJ are initialized to zero, which
987
C might prevent errors in the computation of the integral in the case that
988
C the histograms are plotted several times during a single run (via MCOPY2)
990
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
993
IF(BOOK2(N)(1:1).NE.'Y') RETURN
1001
C-- projection on the x-axis
1004
1 XPROJ(N,I)=XPROJ(N,I)+HIST2(N,I,J)
1005
2 XIN=XIN+XPROJ(N,I)
1006
IF(XIN.EQ.0.) GO TO 10
1007
C-- projection on the y-axis
1010
3 YPROJ(N,J)=YPROJ(N,J)+HIST2(N,I,J)
1016
SUBROUTINE PROHIS(N,M,IAX)
1017
C-- projects the scatter plot N onto the IAX axis (x,y->IAX=1,2) and
1018
C put the contents in histogram M, after automatically booking it.
1020
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
1023
NBIN(M)=NBIN2(IAX,N)
1024
HDEL(M)=HDEL2(IAX,N)
1025
HMIN(M)=HMIN2(IAX,N)
1026
HMAX(M)=HMAX2(IAX,N)
1028
TITLE(M)=TITLE2(N)//'(PROJ)'
1031
HIST(M,I)=XPROJ(N,I)
1032
XHIS(M,I)=XHIS2(N,I)
1033
ELSEIF(IAX.EQ.2) THEN
1034
HIST(M,I)=YPROJ(N,I)
1035
XHIS(M,I)=YHIS2(N,I)
1041
SUBROUTINE MNORM2(N,X)
1043
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
1046
IF(BOOK2(N)(:1).NE.'Y')RETURN
1047
DO 1, I=1,NBIN2(1,N)
1048
DO 1, J=1,NBIN2(2,N)
1049
1 HIST2(N,I,J)=HIST2(N,I,J)/HINT2(N)*X
1053
SUBROUTINE MPRINT2(N)
1055
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
1057
CHARACTER CTIME*10,CDATE*10
1058
INTEGER INI,IMON,IDAY,IYEAR
1061
CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
1064
IF(BOOK2(N)(:1).NE.'Y') RETURN
1065
WRITE(98,7) N,IYEAR,IMON,IDAY,CTIME(1:5)
1066
WRITE(98,*) TITLE2(N)
1067
WRITE(98,10) ((XHIS2(N,I),YHIS2(N,J),HIST2(N,I,J),J=1,NBIN2(2,N))
1069
WRITE(98,20) IENT2(N),HINT2(N),IOSCORE2(N)
1070
7 FORMAT(4X,'HIST = ',I3,' 19',I2,'-',I2,'-',I2,1X,A5/)
1071
10 FORMAT(4X,3E10.3)
1072
20 FORMAT(' ENTRIES=',I5,4X,' INTEGRAL =',E10.3,4x,
1073
* 'OVERSCORE=',I5,//)
1077
SUBROUTINE MTOP2(N,BTIT,LTIT,PLOPT)
1079
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
1080
INTEGER N,I,J,ISEED,ND,NDOTS,NTOT
1081
REAL*8 SCMAX,X,Y,XINT
1082
CHARACTER CTIME*10,CDATE*10
1083
INTEGER INI,IMON,IDAY,IYEAR
1084
CHARACTER*(*) LTIT,BTIT,PLOPT
1085
DOUBLE PRECISION RANXX
1087
C Modified by SF on 26/10/2007, for a better treatment of negative entries.
1088
C When BOX is chosen, negative bins are represented by circles rather
1089
C than boxes. SCMAX need be redefined as the maximum of the absolute
1090
C values of all entries. When DOT is chosen, the number of dots is now
1091
C proportional to the absolute value of the bin entry.
1092
C The layout of the BOX option has been changed. The size of boxes (or squares)
1093
C is now such that, if the cross section is flat, the boxes approximately
1094
C cover the whole plane without overlapping and without leaving empty spaces.
1095
C This assumes that the length of the axes is about 5
1097
ISEED=2143567+2*INT(SECNDS(0.0))
1098
CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
1101
IF(BOOK2(N)(:1).NE.'Y') RETURN
1102
WRITE(99,100) TITLE2(N),BTIT,LTIT,HMIN2(1,N),HMAX2(1,N),
1103
* HMIN2(2,N),HMAX2(2,N)
1105
&' SET WINDOW Y 2.5 TO 9.'/,1X,
1106
&' SET WINDOW X 2.5 TO 10.'/,1X,
1107
&' SET FONT DUPLEX '/1X,
1108
&' TITLE TOP ','"',A,'"',/1X,
1109
&' TITLE BOTTOM ','"',A,'"',/1X,
1110
&' TITLE LEFT ','"',A,'"',/1X,
1111
&' (SET TICKS TOP OFF) '/1x,
1112
&' SET LIMITS X ',F10.5,' ',F10.5,/1X,
1113
&' SET LIMITS Y ',F10.5,' ',F10.5,/1X,
1117
C-- squares with area proportional to the bin weight
1118
IF(PLOPT.EQ.'BOX') THEN
1120
DO 1, I=1,NBIN2(1,N)
1121
DO 1, J=1,NBIN2(2,N)
1122
SCMAX=MAX(SCMAX,ABS(HIST2(N,I,J)))
1124
IF(SCMAX.EQ.0.) GO TO 211
1125
DO 2, I=1,NBIN2(1,N)
1126
DO 2, J=1,NBIN2(2,N)
1127
IF(HIST2(N,I,J).EQ.0.) GO TO 2
1128
IF(HIST2(N,I,J).GT.0.)THEN
1129
WRITE(99,200) XHIS2(N,I),YHIS2(N,J),
1130
# HIST2(N,I,J)/SCMAX*5/NBIN2(1,N)*ABS(HMAX2(1,N)-HMIN2(1,N)),
1131
# HIST2(N,I,J)/SCMAX*5/NBIN2(2,N)*ABS(HMAX2(2,N)-HMIN2(2,N))
1133
WRITE(99,210) XHIS2(N,I),YHIS2(N,J),
1134
# -HIST2(N,I,J)/SCMAX*5/NBIN2(1,N)*ABS(HMAX2(1,N)-HMIN2(1,N)),
1135
# -HIST2(N,I,J)/SCMAX*5/NBIN2(2,N)*ABS(HMAX2(2,N)-HMIN2(2,N))
1139
200 FORMAT(2X,'BOX',F14.8,2X,F14.8,3X,'DATA SIZE',F14.8,1X,F14.8)
1140
210 FORMAT(2X,'CIRCLE',F14.8,2X,F14.8,3X,'DATA SIZE',F14.8,1X,F14.8)
1143
WRITE(99,400) HINT2(N),IENT2(N)
1144
& ,IOSCORE2(N),IMON,IDAY,IYEAR,CTIME(1:5)
1146
&' BOX 6.25 1.0 SIZE 7.5 1.'/,1X,
1147
&' SET WINDOW Y 0. TO 2.'/,1X,
1148
&' SET TITLE SIZE -1.5'/1X,
1149
&' SET FONT DUPLEX '/1X,
1150
&' TITLE 2.8 1.2 "INT =',1PE10.3,' Entries =',I8,2x
1151
& ,'Oversc =',I6,'"',/1X,
1152
&' TITLE 2.8 0.8 "',I2,'-',I2,', 19',I2,2X,A5,'"'/1X,
1153
&' SET TITLE SIZE -2')
1155
500 FORMAT(' NEW PLOT')
1157
C-- dots - as many as the value of HINT2(n), distributed according to wgt
1158
ELSEIF(PLOPT.EQ.'DOT') THEN
1160
DO 11, I=1,NBIN2(1,N)
1161
DO 11, J=1,NBIN2(2,N)
1162
XINT=XINT+HIST2(N,I,J)
1164
WRITE(99,*) ' SET INTENSITY 3'
1165
DO 3, I=1,NBIN2(1,N)
1166
DO 3, J=1,NBIN2(2,N)
1167
IF(HIST2(N,I,J).EQ.0..OR.XINT.EQ.0.) GO TO 3
1168
NDOTS=INT(500*ABS(HIST2(N,I,J))/XINT)
1170
C REST=HIST2(N,I,J)-FLOAT(NDOTS)
1171
C IF(RN.LT.REST) NDOTS=NDOTS+1
1173
X=XHIS2(N,I)+(2.*RANXX(ISEED)-1.)*HDEL2(1,N)
1174
Y=YHIS2(N,J)+(2.*RANXX(ISEED)-1.)*HDEL2(2,N)
1179
WRITE(99,*) ' SET INTENSITY 1'
1181
WRITE(99,600) HINT2(N),IENT2(N),NTOT
1182
& ,IMON,IDAY,IYEAR,CTIME(1:5)
1184
&' BOX 6.25 1.0 SIZE 7.5 1.'/,1X,
1185
&' SET WINDOW Y 0. TO 2.'/,1X,
1186
&' SET TITLE SIZE -1.5'/1X,
1187
&' SET FONT DUPLEX '/1X,
1188
&' TITLE 2.8 1.2 "INT =',1PE10.3,' Entries =',I8,2x
1189
& ,'Points =',I6,'"',/1X,
1190
&' TITLE 2.8 0.8 "',I2,'-',I2,', 19',I2,2X,A5,'"'/1X,
1191
&' SET TITLE SIZE -2')
1198
FUNCTION RANXX(SEED)
1200
* Ref.: K. Park and K.W. Miller, Comm. of the ACM 31 (1988) p.1192
1201
* Use seed = 1 as first value.
1203
IMPLICIT INTEGER(A-Z)
1204
DOUBLE PRECISION MINV,RANXX
1206
PARAMETER(M=2147483647,A=16807,Q=127773,R=2836)
1207
PARAMETER(MINV=0.46566128752458d-09)
1211
IF(SEED.LE.0) SEED = SEED + M
1216
SUBROUTINE MULTITOP(NH,NE,N,M,BTIT,LTIT,SCA)
1218
INCLUDE 'mcatnlo_hbook_gfortran8.inc'
1220
REAL*8 LABS,FMX,FMN,FMAX,FMIN,SRED,TICS,TITS,X,XD,YD,XL,YL,
1221
# XU,YU,XMX,XTIT,YTIT,XTIT0,YTIT0
1222
CHARACTER CTIME*10,CDATE*10,SCALE*3
1223
INTEGER IMON,IDAY,IYEAR,I,J,IBIN,IFRAME,IFRMAX,INI,IP,MOLD,
1225
CHARACTER*(*) LTIT,BTIT,SCA
1227
DATA PLOT/'SOLID','DASHES','DOTS','DOTDASH'/
1228
C PLOT SIZE, CORNERS
1229
REAL*8 WIDTH,HEIGHT,XCORN,YCORN
1230
DATA WIDTH,HEIGHT/11.5,8.5/,XCORN,YCORN/1.5,1./
1231
C PLOT VERSUS TEXT FRACTION
1232
REAL*8 XPFRAC,YPFRAC,XTFRAC,YTFRAC
1233
DATA XPFRAC,YPFRAC/0.75,0.75/,XTFRAC,YTFRAC/0.25,0.25/
1235
REAL*8 TIT0,LAB0,TIC0
1236
DATA TIT0,LAB0,TIC0/3.,3.,0.06/
1239
CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
1241
WRITE(99,71) CDATE,CTIME(1:5)
1242
71 FORMAT(4X,' ( ',A10,1X,A5/)
1245
IF(SCA.EQ.'REF') THEN
1249
IF(BOOK(NH)(:1).NE.'Y') RETURN
1252
IF(IFRAME.GT.IFRMAX.OR.N.NE.NOLD.OR.M.NE.MOLD) THEN
1255
WRITE(99,1) CDATE,CTIME(1:5)
1256
1 FORMAT(' ( SET FONT DUPLEX',/,' SET TITLE SIZE 2',/,
1257
+ ' TITLE 12.8 9 ANGLE -90 ','" MLM ',A10,1X,A5,'"')
1259
IF(IFRAME.EQ.1) THEN
1262
ELSEIF(IFRAME.LE.IFRMAX) THEN
1269
IF(N.EQ.NOLD) GO TO 10
1272
SRED=SQRT(DFLOAT(N*M))
1276
XTIT0=0.55*XPFRAC*XD
1278
10 IF(M.EQ.MOLD) GO TO 20
1294
FMX = X + HIST(NE,IBIN)
1295
FMN = X - HIST(NE,IBIN)
1297
FMX=MAX(FMX,X + HIST(NE,IBIN))
1298
FMN=MIN(FMN,X - HIST(NE,IBIN))
1301
XMX=MAX(XMX,ABS(X)+ HIST(NE,IBIN))
1304
50 IF(SCALE.EQ.'LIN') THEN
1305
IF(FMN.GE.0.) FMIN=0.
1306
IF(FMN.LT.0.) FMIN=FMN*1.3
1307
IF(FMX.GT.0.) FMAX=FMX*1.3
1308
IF(FMX.LT.0.) FMAX=0.
1309
ELSEIF(SCALE.EQ.'LOG') THEN
1311
c$$$ IF(FMN.LE.0.) THEN
1315
FMAX=10.**( AINT(LOG10(ABS(FMX))+1000001) - 1000000 )
1316
FMIN=10.**( AINT(LOG10(ABS(FMN))+1000000) - 1000000 )
1319
WRITE(99,100) TITS,LABS,TICS,XL,XU,YL,YU
1320
100 FORMAT(2X,'( SET FONT DUPLEX',/,
1321
* 2X,'SET TITLE SIZE ',F8.4,/,
1322
* 2X,'SET LABEL SIZE ',F8.4,/,
1323
* 2X,'SET TICKS TOP OFF SIZE ',F8.4,/,
1324
* 2X,'SET WINDOW X ',F8.4,' TO ',F8.4,/,
1325
* 2X,'SET WINDOW Y ',F8.4,' TO ',F8.4)
1328
WRITE(99,101) XL,YTIT,TITLE(NH)(1:40)
1329
101 FORMAT(' TITLE ',2(F8.4,1X),'"',A,'"')
1331
WRITE(99,102) XTIT,YTIT,HINT(NH)
1332
102 FORMAT(' TITLE ',2(F8.4,1X),'" INT=',1PE10.3,'"')
1334
WRITE(99,103) XTIT,YTIT,IENT(NH)
1335
103 FORMAT(' TITLE ',2(F8.4,1X),'" ENT=',I9,'"')
1337
IF(USCORE(NH).NE.0.) THEN
1338
WRITE(99,104) XTIT,YTIT,USCORE(NH)
1339
104 FORMAT(' TITLE ',2(F8.4,1X),'" UFL=',1PE10.3,'"')
1342
IF(OSCORE(NH).NE.0.) THEN
1343
WRITE(99,105) XTIT,YTIT,OSCORE(NH)
1344
105 FORMAT(' TITLE ',2(F8.4,1X),'" OFL=',1PE10.3,'"')
1347
WRITE(99,106) XTIT,YTIT,XU,YTIT,XTIT,YTIT,XTIT,YU
1348
106 FORMAT(2X,'SET ORD X Y ',/,2(F8.4,1X),/,2(F8.4,1X),/,
1350
* 2X,2(F8.4,1X),/,2(F8.4,1X),/,
1352
WRITE(99,108) TITS*1.5
1353
108 FORMAT(2X,'SET TITLE SIZE ',F8.4)
1354
WRITE(99,107) BTIT,XL-0.75*XD*XTFRAC,YL+(YU-YL)/3.,LTIT,SCALE,
1355
* HMIN(NH),HMAX(NH),FMIN,FMAX
1357
&' TITLE BOTTOM ','"',A,'"',/1X,
1358
&' TITLE ',f10.5,f10.5,' ANGLE 90 ','"',A,'"',/1X,
1359
&' SET SCALE Y ',A,/1X,
1360
&' SET TICKS TOP OFF '/1x,
1361
&' SET LIMITS X ',F10.5,' ',F10.5,/1X,
1362
&' SET LIMITS Y ',1PE10.3,' ',1PE10.3,/1X,
1363
&' SET ORDER X Y DY')
1365
C END HEADER , FILL TOPDRAWER WITH DATA
1367
ENTRY MTFILL(NH,NE,N,M,BTIT,LTIT,SCA)
1370
WRITE(99,110) TITLE(NH),HINT(NH),IENT(NH)
1371
110 FORMAT(' ( ',A,/,' ( INT=',1PE10.3,' ENTRIES=',I12)
1372
DO 200 IBIN=1,NBIN(NH)
1373
WRITE(99,'(3X,F10.4,2(2X,E15.4))')
1374
& XHIS(NH,IBIN),HIST(NH,IBIN),HIST(NE,IBIN)
1376
WRITE(99,201) PLOT(IP)
1377
IF(BOOK(NE).NE.'NO') WRITE(99,*) ' PLOT'
1378
201 FORMAT(2X,'HIST ',A)
1379
202 FORMAT(' NEW PLOT',/,/)
1385
202 FORMAT(' NEW PLOT',/,/)
1386
CALL MULTITOP(1,1,1,1,' ',' ','REF')
1389
C*******************************************************************
1390
C END OF THE HISTOGRAMMING PACKAGE
1391
C*******************************************************************
1394
C***** GENERIC UTILITIES *********************************
1396
C VSUM(P,Q,K) returns K(I)=P(I)+Q(I) , I=1,4
1397
C DOT(P,Q) Lorentz scalar product (+---)
1398
C RNDINT(N,INDX) returns a random ordering of 1,2,..,N
1399
C UTSORT(A,N,K,IOPT):
1400
C Sort A(N) into ascending order
1401
C IOPT = 1 : return sorted A and index array K
1402
C IOPT = 2 : return index array K only
1404
SUBROUTINE VSUM(P1,P2,Q)
1406
REAL*8 P1(4),P2(4),Q(4)
1413
SUBROUTINE RNDINT(N,INDX)
1414
C returns a random ordering of the string of integers (1,2,..,N)
1417
INTEGER INDX(N),ISEED,I
1419
DOUBLE PRECISION RANXX
1421
IF(N.GT.500) WRITE(6,*) 'WARNING**RNDINT'
1425
CALL UTSORT(X,N,INDX,2)
1428
SUBROUTINE UTSORT(A,N,K,IOPT)
1429
C Sort A(N) into ascending order
1430
C IOPT = 1 : return sorted A and index array K
1431
C IOPT = 2 : return index array K only
1432
C------------------------------------------------------------------------
1436
INTEGER K(N),IL(500),IR(500),I,J
1437
IF (N.GT.500) WRITE(6,*) 'WARNING**UTSORT'
1444
2 IF(A(I).GT.A(J)) GO TO 5
1445
3 IF(IL(J).EQ.0) GO TO 4
1451
5 IF(IR(J).LE.0) GO TO 6
1461
8 IF(IL(J).GT.0) GO TO 20
1470
30 IF(IOPT.EQ.2) RETURN
1475
SUBROUTINE MBWARN(ROUT)
1477
WRITE(*,*) '***********************************************'
1478
WRITE(*,*) '***** WARNING CALLED FROM ROUTINE ',ROUT,':'
1482
SUBROUTINE MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
1484
INTEGER IMON,IDAY,IYEAR,IVAL(3)
1485
CHARACTER*10 CTIME,CDATE
1486
CHARACTER*10 CTIME2,CDATE2
1491
CALL DATE_AND_TIME(TIME=CTIME2)
1492
CALL DATE_AND_TIME(DATE=CDATE2)
1493
CTIME=CTIME2(1:2)//":"//CTIME2(3:4)
1494
CDATE=CDATE2(1:4)//"-"//CDATE2(5:6)//"-"//CDATE2(7:8)