~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook_gfortran8.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
C####
 
2
C#### YOU WILL NOT NEED WHAT FOLLOWS IN YOUR OWN ANALYSIS
 
3
C####
 
4
C************************************************************************
 
5
C
 
6
C                           HISTOGRAMMING PACKAGE 
 
7
C                                 M. MANGANO
 
8
C
 
9
C************************************************************************
 
10
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,
 
17
C                               MTFILL,...)
 
18
C**********************************************************************
 
19
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).
 
24
C
 
25
C INITIALIZATION:
 
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
 
33
C Example:
 
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
 
38
C
 
39
C FILLING:
 
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. 
 
42
C
 
43
C PLAYING AROUND:
 
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.
 
71
C
 
72
C FINAL ACCOUNTING:
 
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.
 
78
C
 
79
C OUTPUT:
 
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--------------------------------------------------------------------------
 
97
C
 
98
C  COMMON/HISTO/  Histogram N
 
99
C                           
 
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
 
103
C
 
104
C   HMIN(N),      Min value of x range
 
105
C   HMAX(N),      Max value of x range
 
106
C   HDEL(N),      Bin width
 
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
 
120
C
 
121
 
 
122
 
 
123
      SUBROUTINE INIHIST
 
124
      IMPLICIT NONE
 
125
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
126
      INTEGER I
 
127
      NHIST=0
 
128
      NHIST2=0
 
129
      DO 1, I=1,NMB             
 
130
   1  BOOK(I)=' NO'
 
131
      DO 2 I=1,50
 
132
   2  BOOK2(I)=' NO'
 
133
      END  
 
134
    
 
135
      SUBROUTINE MBOOK(N,TIT,DEL,XMIN,XMAX)
 
136
      IMPLICIT NONE
 
137
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
138
      INTEGER N
 
139
      CHARACTER*(*) TIT
 
140
      REAL*8 DEL,XMIN,XMAX
 
141
      INTEGER I
 
142
      NHIST = MAX(N,NHIST)
 
143
      IF(BOOK(N)(1:1).EQ.'Y') THEN
 
144
         CALL MBWARN('MBOOK')
 
145
         WRITE(*,*) 'Histogram',N,TITLE(N),' already in use. '
 
146
         WRITE(*,*) 'superseded by ',TIT
 
147
      ENDIF
 
148
      BOOK(N) = 'YES'
 
149
      TITLE(N) = ' '//TIT
 
150
1     HDEL(N) = DEL
 
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,' ):'
 
155
        READ(*,*) DEL
 
156
        GO TO 1
 
157
      ENDIF
 
158
      HMIN(N) = XMIN
 
159
      HMAX(N) = NBIN(N)*DEL+XMIN
 
160
      IF(ABS(HMAX(N)-XMAX).GT.1.D-3*DEL) THEN
 
161
         CALL MBWARN('MBOOK')
 
162
         WRITE(*,*)
 
163
     #'Histogram ', TIT, ' Change of upper limit:',xmax,'-->',HMAX(N)
 
164
      ENDIF
 
165
      IENT(N) = 0
 
166
      IUSCORE(N) = 0
 
167
      IOSCORE(N) = 0
 
168
      USCORE(N) = 0
 
169
      OSCORE(N) = 0
 
170
      HAVG(N) = 0
 
171
      HINT(N) = 0
 
172
      HSIG(N) = 0
 
173
      DO I=1,NBIN(N)
 
174
         XHIS(N,I)=HMIN(N)+HDEL(N)*(DFLOAT(I)-0.5)
 
175
         IHIS(N,I)=0
 
176
         HIST(N,I)=0
 
177
      ENDDO
 
178
      END
 
179
 
 
180
      SUBROUTINE MFILL(N,X,Y)
 
181
      IMPLICIT NONE
 
182
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
183
      INTEGER N,I
 
184
      REAL*8 X,Y,SUM,XI
 
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
 
191
      ELSE
 
192
         XI=((X-HMIN(N))/HDEL(N))+1
 
193
         I=INT(XI)
 
194
         IENT(N)=IENT(N)+1
 
195
         IHIS(N,I)=IHIS(N,I)+1
 
196
C SF, 10/3/2010 -- avoids numerical problems
 
197
         SUM=HIST(N,I)+Y
 
198
         IF(ABS(SUM).LT.1.d-19)THEN
 
199
           HIST(N,I)=0.d0
 
200
         ELSE
 
201
           HIST(N,I)=SUM
 
202
         ENDIF
 
203
      ENDIF
 
204
      END
 
205
 
 
206
 
 
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.
 
215
      IMPLICIT NONE
 
216
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
217
      INTEGER NIN,NOUT,IDIR,IPOW,I,L,M
 
218
      REAL*8 SUMPOW
 
219
      CHARACTER*14  C
 
220
      DIMENSION C(2) 
 
221
      DATA C/' INTEG BELOW X',' INTEG ABOVE X'/
 
222
      IF(BOOK(NIN)(1:1).NE.'Y') RETURN
 
223
      M = NBIN(NIN)                                           
 
224
      I = (IDIR + 3)/2
 
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))
 
228
      ENDIF
 
229
      IF(IDIR.EQ.1) THEN
 
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
 
233
         DO L=2,M                      
 
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
 
237
         ENDDO
 
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
 
244
         DO L=M-1,1,-1                        
 
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
 
248
         ENDDO
 
249
         USCORE(NOUT) = SUMPOW(USCORE(NIN),HIST(NIN,1),IPOW)
 
250
         IUSCORE(NOUT) = IUSCORE(NIN)+IHIS(NIN,1)
 
251
      ELSE                                 
 
252
         CALL MBWARN('MINTEG')
 
253
         WRITE(*,*) 'Wrong idir in minteg: OPERATION NOT PERFORMED'
 
254
         RETURN
 
255
      ENDIF
 
256
      END
 
257
 
 
258
      FUNCTION SUMPOW(X,Y,IPOW)
 
259
      IMPLICIT NONE
 
260
      REAL*8 SUMPOW,X,Y
 
261
      INTEGER IPOW
 
262
      IF(IPOW.EQ.1) THEN
 
263
         SUMPOW = X + Y
 
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'
 
269
      ELSE
 
270
         SUMPOW = (X**IPOW+Y**IPOW)**(1./IPOW)
 
271
      ENDIF
 
272
      END
 
273
 
 
274
      SUBROUTINE MOPERA(I,OPER,J,K,X,Y)
 
275
      IMPLICIT NONE
 
276
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
277
      INTEGER I,J,K
 
278
      REAL*8 X,Y
 
279
      INTEGER L
 
280
      REAL*8 XAVG,XNORM,XSQAVG,XSUM,XSUMSQ,XXX
 
281
      CHARACTER OPER*1
 
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')
 
285
          WRITE(*,*) I,J                               
 
286
  20      FORMAT(' ****** INCOMPATIBLE OPERATION HIST ',I2,' &',I2,
 
287
     &    '*******'/)
 
288
          RETURN
 
289
      ENDIF
 
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))
 
293
      IF(OPER.EQ.'E') THEN
 
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))
 
299
         ELSE
 
300
            XXX = 0
 
301
         ENDIF
 
302
         DO L=1,NBIN(I)
 
303
            XSUM   = HIST(I,L)
 
304
            XSUMSQ = HIST(J,L)
 
305
            HIST(K,L)=HIST(K,L) + XXX*HIST(I,L)
 
306
            HIST(J,L)=XXX*SQRT(ABS(XSUMSQ-XSUM**2*XXX))
 
307
         ENDDO
 
308
         IENT(K)=IENT(K)+IENT(I)
 
309
         XSUM = USCORE(I)
 
310
         XSUMSQ = USCORE(J)
 
311
         USCORE(K) = USCORE(K)+XXX*XSUM
 
312
         USCORE(J) = XXX*SQRT(ABS(XSUMSQ-XSUM**2*XXX))
 
313
         XSUM = OSCORE(I)
 
314
         XSUMSQ = OSCORE(J)
 
315
         OSCORE(K) = OSCORE(K)+XXX*XSUM
 
316
         OSCORE(J) = XXX*SQRT(ABS(XSUMSQ-XSUM**2*XXX))
 
317
      ELSEIF(OPER.EQ.'Q') THEN
 
318
         DO L=1,NBIN(I)
 
319
            HIST(K,L) = SQRT(HIST(J,L)**2+HIST(I,L)**2)
 
320
         ENDDO
 
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
 
324
         DO L=1,NBIN(I)
 
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
 
329
            IENT(K) = IENT(K)+1
 
330
            HIST(I,L) = 0
 
331
            IHIS(I,L) = 0
 
332
         ENDDO
 
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
 
342
         IENT(I) = 0
 
343
         IUSCORE(I) = 0
 
344
         IOSCORE(I) = 0
 
345
         USCORE(I) = 0
 
346
         OSCORE(I) = 0
 
347
      ELSE
 
348
        DO L=1,NBIN(I)
 
349
        IF(OPER.EQ.'+') THEN
 
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
 
357
            HIST(K,L)=0.
 
358
          ELSE
 
359
            HIST(K,L)=X*HIST(I,L) / (Y*HIST(J,L))
 
360
          ENDIF
 
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))
 
366
          ELSE                           
 
367
            HIST(K,L)=0.
 
368
          ENDIF
 
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
 
373
             HIST(K,L)=0.
 
374
           ELSE
 
375
             HIST(K,L)=X*LOG10(Y*HIST(I,L))
 
376
           ENDIF
 
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))
 
380
           IF(XNORM.NE.0.) THEN
 
381
             XAVG=HIST(J,L)/XNORM
 
382
             HIST(K,L)=
 
383
     &       SQRT(ABS(-XAVG**2+HIST(K,L)/XNORM)/DFLOAT(IHIS(I,L)))
 
384
             HIST(J,L)=XAVG 
 
385
           ELSE 
 
386
             HIST(K,L)=0.
 
387
             HIST(J,L)=0.
 
388
           ENDIF
 
389
        ELSEIF(OPER.EQ.'V') THEN                 
 
390
           XAVG=HIST(I,L)*X
 
391
           XSQAVG=HIST(J,L)*X
 
392
           XNORM=DFLOAT(IHIS(I,L))*X
 
393
           IF(XNORM.NE.0.) THEN
 
394
              HIST(K,L)=SQRT(ABS(XSQAVG-XAVG**2)/XNORM)
 
395
              HIST(I,L)=XAVG
 
396
           ELSE  
 
397
              HIST(K,L)=0.
 
398
           ENDIF 
 
399
        ELSE 
 
400
         CALL MBWARN('MOPERA')
 
401
         WRITE(*,*) OPER
 
402
   5     FORMAT(' ****** OPERATION ="',A1,'" UNKNOWN ********'/)
 
403
         RETURN
 
404
        ENDIF
 
405
        ENDDO
 
406
        IF(OPER.EQ.'+') THEN
 
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
 
417
            USCORE(K)=0.
 
418
          ELSE
 
419
            USCORE(K)=X*USCORE(I) / (Y*USCORE(J))
 
420
          ENDIF
 
421
          IF(Y.EQ.0..OR.OSCORE(J).EQ.0.) THEN
 
422
            OSCORE(K)=0.
 
423
          ELSE
 
424
            OSCORE(K)=X*OSCORE(I) / (Y*OSCORE(J))
 
425
          ENDIF
 
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))
 
432
          ELSE                           
 
433
            USCORE(K)=0.
 
434
          ENDIF     
 
435
          IF(OSCORE(I).GT.0.) THEN
 
436
            OSCORE(K)=X*SQRT(OSCORE(I))
 
437
          ELSE                           
 
438
            OSCORE(K)=0.
 
439
          ENDIF     
 
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
 
445
             USCORE(K)=0.
 
446
           ELSE
 
447
             USCORE(K)=X*LOG10(Y*USCORE(I))
 
448
           ENDIF                         
 
449
          IF(OSCORE(I).EQ.0..OR.J.EQ.0.) THEN
 
450
             OSCORE(K)=0.
 
451
           ELSE
 
452
             OSCORE(K)=X*LOG10(Y*OSCORE(I))
 
453
           ENDIF                         
 
454
        ENDIF
 
455
      ENDIF
 
456
      RETURN
 
457
      END
 
458
 
 
459
      SUBROUTINE MCOPY(NIN,NOUT)
 
460
      IMPLICIT NONE
 
461
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
462
      INTEGER NIN,NOUT,I
 
463
C
 
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'
 
472
      HDEL(NOUT)=HDEL(NIN)
 
473
      HMIN(NOUT)=HMIN(NIN)
 
474
      HMAX(NOUT)=HMAX(NIN)
 
475
      NBIN(NOUT)=NBIN(NIN)
 
476
      IENT(NOUT)=IENT(NIN)
 
477
      IUSCORE(NOUT)=IUSCORE(NIN)
 
478
      IOSCORE(NOUT)=IOSCORE(NIN)
 
479
      USCORE(NOUT)=USCORE(NIN)
 
480
      OSCORE(NOUT)=OSCORE(NIN)
 
481
      HAVG(NOUT)=HAVG(NIN)
 
482
      HINT(NOUT)=HINT(NIN)
 
483
      HSIG(NOUT)=HSIG(NIN)
 
484
      DO I=1,NBIN(NOUT)
 
485
        HIST(NOUT,I)=HIST(NIN,I)
 
486
        XHIS(NOUT,I)=XHIS(NIN,I)
 
487
        IHIS(NOUT,I)=IHIS(NIN,I)
 
488
      ENDDO
 
489
      END
 
490
 
 
491
      SUBROUTINE MZERO(N)
 
492
      IMPLICIT NONE
 
493
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
494
      INTEGER N,I
 
495
      BOOK(N)='RES'
 
496
      IENT(N)=0
 
497
      IUSCORE(N)=0
 
498
      IOSCORE(N)=0
 
499
      HAVG(N)=0.
 
500
      HINT(N)=0.
 
501
      DO 1 I=1,NBIN(N)
 
502
   1  HIST(N,I)=0.
 
503
      END                   
 
504
 
 
505
      SUBROUTINE MDUMP(N)
 
506
      IMPLICIT NONE
 
507
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
508
      INTEGER N
 
509
      REAL*8 HISTL(100),XHISL(100),HDELL,HMINL,HMAXL,USCOREL,OSCOREL,
 
510
     &HAVGL,HINTL,HSIGL
 
511
      INTEGER IHISL(100),NBINL,IUSCOREL,IOSCOREL,IENTL,I
 
512
      CHARACTER*50 TITLEL
 
513
      CHARACTER*3 BOOKL
 
514
      SAVE HISTL,XHISL,HDELL,HMINL,HMAXL,USCOREL,OSCOREL,
 
515
     &NBINL,IHISL,IUSCOREL,IOSCOREL,IENTL,HAVGL,HINTL,HSIGL,
 
516
     &BOOKL,TITLEL
 
517
      BOOKL=BOOK(N)                                                           
 
518
      TITLEL=TITLE(N)                   
 
519
      HDELL=HDEL(N)  
 
520
      HMINL=HMIN(N)
 
521
      HMAXL=HMAX(N)
 
522
      USCOREL=USCORE(N)
 
523
      OSCOREL=OSCORE(N)
 
524
      NBINL=NBIN(N)
 
525
      HSIGL=HSIG(N)
 
526
      IENTL=IENT(N)     
 
527
      IUSCOREL=IUSCORE(N)
 
528
      IOSCOREL=IOSCORE(N)      
 
529
      HAVGL=HAVG(N)         
 
530
      HINTL=HINT(N)   
 
531
      DO 1 I=1,NBIN(N)
 
532
      XHISL(I)=XHIS(N,I)
 
533
      IHISL(I)=IHIS(N,I)
 
534
   1  HISTL(I)=HIST(N,I)
 
535
      RETURN
 
536
      ENTRY MFETCH(N)
 
537
      BOOK(N)=BOOKL                   
 
538
      TITLE(N)=TITLEL                   
 
539
      HDEL(N)=HDELL  
 
540
      HMIN(N)=HMINL
 
541
      HMAX(N)=HMAXL
 
542
      USCORE(N)=USCOREL
 
543
      OSCORE(N)=OSCOREL
 
544
      NBIN(N)=NBINL
 
545
      HSIG(N)=HSIGL
 
546
      IENT(N)=IENTL     
 
547
      IUSCORE(N)=IUSCOREL
 
548
      IOSCORE(N)=IOSCOREL      
 
549
      HAVG(N)=HAVGL         
 
550
      HINT(N)=HINTL   
 
551
      DO 2 I=1,NBINL
 
552
      XHIS(N,I)=XHISL(I)
 
553
      IHIS(N,I)=IHISL(I)
 
554
   2  HIST(N,I)=HISTL(I)
 
555
      END               
 
556
 
 
557
      SUBROUTINE MRESET(N)
 
558
      IMPLICIT NONE
 
559
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
560
      INTEGER N
 
561
      BOOK(N)='RES'
 
562
      END
 
563
 
 
564
      SUBROUTINE PUTTAG(J,NAME)
 
565
      IMPLICIT NONE
 
566
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
567
      INTEGER J
 
568
c Per marcare un istogramma
 
569
      CHARACTER * (*) NAME, TAG
 
570
      BOOK(J) = NAME
 
571
      RETURN
 
572
      ENTRY GETTAG(J,TAG)
 
573
      TAG = BOOK(J)
 
574
      END
 
575
 
 
576
      SUBROUTINE MFINAL(N)
 
577
      IMPLICIT NONE
 
578
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
579
      INTEGER N,J
 
580
      REAL*8 AVG,XIN,SIG,IF,X
 
581
c
 
582
      AVG=0
 
583
      XIN=0                                  
 
584
      SIG=0
 
585
      IF=0
 
586
      DO J=1,NBIN(N)
 
587
         X=HIST(N,J)
 
588
         AVG=AVG+X*XHIS(N,J)
 
589
         XIN=XIN+X
 
590
         IF(X.NE.0) IF=1
 
591
      ENDDO             
 
592
      IF(XIN.EQ.0) GO TO 10
 
593
      AVG = AVG/XIN
 
594
      DO J=1,NBIN(N)
 
595
         SIG=HIST(N,J)*(XHIS(N,J)-AVG)**2+SIG
 
596
      ENDDO
 
597
      SIG=SQRT(ABS(SIG/XIN))
 
598
 10   CONTINUE
 
599
      HINT(N) = XIN
 
600
      HAVG(N) = AVG
 
601
      HSIG(N) = SIG
 
602
      IF(IF.EQ.0) BOOK(N)='RES'
 
603
      END               
 
604
 
 
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
 
608
      IMPLICIT NONE
 
609
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
610
      INTEGER N,J
 
611
      REAL*8 AVG,XIN,SIG,IF,X
 
612
      IF(BOOK(N)(1:1).NE.'Y') RETURN
 
613
      AVG=0
 
614
      XIN=0                                  
 
615
      SIG=0
 
616
      IF=0
 
617
      IF(HMIN(N).NE.HMAX(N))IF=1
 
618
      DO J=1,NBIN(N)
 
619
         X=HIST(N,J)
 
620
         AVG=AVG+X*XHIS(N,J)
 
621
         XIN=XIN+X
 
622
         IF(X.NE.0) IF=1
 
623
      ENDDO             
 
624
      IF(XIN.EQ.0) GO TO 10
 
625
      AVG = AVG/XIN
 
626
      DO J=1,NBIN(N)
 
627
         SIG=HIST(N,J)*(XHIS(N,J)-AVG)**2+SIG
 
628
      ENDDO
 
629
      SIG=SQRT(ABS(SIG/XIN))
 
630
 10   CONTINUE
 
631
      HINT(N) = XIN
 
632
      HAVG(N) = AVG
 
633
      HSIG(N) = SIG
 
634
      IF(IF.EQ.0) BOOK(N)='RES'
 
635
      END               
 
636
 
 
637
      SUBROUTINE MNORM(N,X)
 
638
      IMPLICIT NONE
 
639
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
640
      INTEGER N,I
 
641
      REAL*8 X,Y
 
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'
 
647
        RETURN
 
648
      ELSE
 
649
        Y=X/HINT(N)
 
650
      ENDIF
 
651
      DO 1, I=1,NBIN(N)
 
652
    1 HIST(N,I)=HIST(N,I)*Y
 
653
      HINT(N)=X
 
654
      OSCORE(N)=OSCORE(N)*Y
 
655
      USCORE(N)=USCORE(N)*Y
 
656
      END
 
657
 
 
658
      SUBROUTINE MPRINT(N)
 
659
      IMPLICIT NONE
 
660
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
661
      INTEGER N,J
 
662
      CHARACTER CTIME*10,CDATE*10
 
663
      INTEGER INI,IMON,IDAY,IYEAR
 
664
      DATA INI/0/
 
665
      IF(INI.EQ.0) THEN
 
666
      CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
 
667
      INI=1
 
668
      ENDIF
 
669
      IF(BOOK(N)(:1).NE.'Y') RETURN
 
670
      WRITE(98,7) N,IYEAR,IMON,IDAY,CTIME(1:5)
 
671
      WRITE(98,*) TITLE(N)
 
672
      DO 1 J=1,NBIN(N)
 
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)
 
676
    1 CONTINUE
 
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/)
 
680
   10 FORMAT(4X,2E10.3)
 
681
   15 FORMAT(/' AVG =',E10.3,4X,' RMS =',E10.3,' INTEGRAL =',E10.3,/)
 
682
   20 FORMAT(' ENTRIES=',I5,4X,'UNDERSCORE=',I5,4x,'OVERSCORE=',I5,//)
 
683
      END
 
684
 
 
685
      SUBROUTINE MTOP(N,M,BTIT,LTIT,SCALE)
 
686
      IMPLICIT NONE
 
687
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
688
      INTEGER N,M,J
 
689
      CHARACTER CTIME*10,CDATE*10
 
690
      INTEGER INI,IMON,IDAY,IYEAR
 
691
      CHARACTER*(*) LTIT,BTIT,SCALE
 
692
      DATA INI/0/
 
693
      IF(INI.EQ.0) THEN
 
694
      CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
 
695
      INI=1
 
696
      ENDIF
 
697
      IF(BOOK(N)(:1).NE.'Y') RETURN
 
698
      WRITE(99,100) TITLE(N),BTIT,LTIT,SCALE,HMIN(N),HMAX(N)
 
699
  100 FORMAT( /1x,                               
 
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 ')
 
711
      DO 1 J=1,NBIN(N)
 
712
      WRITE(99,'(3X,F10.4,2(2X,E15.4))')  
 
713
     &                            XHIS(N,J),HIST(N,J),HIST(M,J)
 
714
    1 CONTINUE
 
715
      WRITE(99,200)
 
716
  200 FORMAT('   PLOT')
 
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)
 
719
  300 FORMAT( /1x,                               
 
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=',
 
729
     &      1PE10.3,'"',/1X,
 
730
     &' TITLE 7.5 0.6 "',I2,'-',I2,', 19',I2,2X,A5,'"',/1X,                
 
731
     &' SET TITLE SIZE -2')                            
 
732
      WRITE(99,400)
 
733
  400 FORMAT('   NEW PLOT')
 
734
      END
 
735
 
 
736
      SUBROUTINE MNRTOP(N,M,BTIT,LTIT,SCALE)
 
737
      IMPLICIT NONE
 
738
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
739
      INTEGER N,M,J
 
740
      CHARACTER CTIME*10,CDATE*10
 
741
      INTEGER INI,IMON,IDAY,IYEAR
 
742
      CHARACTER*(*) LTIT,BTIT,SCALE
 
743
      DATA INI/0/
 
744
      IF(INI.EQ.0) THEN
 
745
         CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
 
746
         INI=1
 
747
      ENDIF
 
748
      CALL MFINAL(N)
 
749
      WRITE(99,100) TITLE(N)(1:20),BTIT,LTIT,SCALE,HMIN(N),HMAX(N)
 
750
  100 FORMAT( /1x,                               
 
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 ')
 
762
      DO 1 J=1,NBIN(N)
 
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)
 
766
    1 CONTINUE
 
767
      WRITE(99,200)
 
768
  200 FORMAT('   PLOT')
 
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)
 
771
  300 FORMAT( /1x,                               
 
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=',
 
781
     &      1PE10.3,'"',/1X,
 
782
     &' TITLE 7.5 0.6 "',I2,'-',I2,', 19',I2,2X,A5,'"',/1X,                
 
783
     &' SET TITLE SIZE -2')
 
784
      WRITE(99,400)
 
785
  400 FORMAT('   NEW PLOT')
 
786
      END
 
787
 
 
788
 
 
789
      SUBROUTINE MBOOK2(N,TIT,DELX,XMIN,XMAX,DELY,YMIN,YMAX)
 
790
      IMPLICIT NONE
 
791
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
792
      INTEGER N
 
793
      CHARACTER*(*) TIT
 
794
      REAL*8 DELX,XMIN,XMAX,DELY,YMIN,YMAX
 
795
      INTEGER I,J
 
796
      NHIST2=MAX(N,NHIST2)
 
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
 
801
      ENDIF
 
802
      TITLE2(N)='   '//TIT                     
 
803
      BOOK2(N)='YES'
 
804
C-- setup x-boundaries
 
805
      HDEL2(1,N)=DELX
 
806
      HMIN2(1,N)=XMIN  
 
807
      HMAX2(1,N)=XMAX
 
808
      NBIN2(1,N)=INT((XMAX-XMIN)/(DELX*0.9999))
 
809
C-- setup y-boundaries
 
810
      HDEL2(2,N)=DELY
 
811
      HMIN2(2,N)=YMIN  
 
812
      HMAX2(2,N)=YMAX
 
813
      NBIN2(2,N)=INT((YMAX-YMIN)/(DELY*0.9999))
 
814
      IENT2(N)=0                      
 
815
      IOSCORE2(N)=0
 
816
      HAVG2(1,N)=0.
 
817
      HAVG2(2,N)=0.
 
818
      HINT2(N)=0.
 
819
      DO 1 I=1,NBIN2(1,N)
 
820
      XHIS2(N,I)=HMIN2(1,N)+HDEL2(1,N)*(DFLOAT(I)-0.5)
 
821
   1  CONTINUE
 
822
      DO 2 I=1,NBIN2(2,N)
 
823
      YHIS2(N,I)=HMIN2(2,N)+HDEL2(2,N)*(DFLOAT(I)-0.5)
 
824
   2  CONTINUE
 
825
      DO 3 I=1,NBIN2(1,N)
 
826
      DO 3 J=1,NBIN2(2,N)
 
827
      HIST2(N,I,J)=0.                   
 
828
   3  CONTINUE
 
829
      END        
 
830
 
 
831
      SUBROUTINE MFILL2(N,X,Y,WGT)
 
832
      IMPLICIT NONE
 
833
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
834
      INTEGER N
 
835
      REAL*8 X,Y,WGT
 
836
      INTEGER I,J
 
837
      REAL*8 XI,YI
 
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 
 
842
      ELSE
 
843
         XI=((X-HMIN2(1,N))/HDEL2(1,N))+1.
 
844
         YI=((Y-HMIN2(2,N))/HDEL2(2,N))+1.
 
845
         I=INT(XI)
 
846
         J=INT(YI)
 
847
         IF( I.GT.0.AND.I.LE.NBIN2(1,N) .AND.
 
848
     #       J.GT.0.AND.J.LE.NBIN2(2,N) )THEN
 
849
           IENT2(N)=IENT2(N)+1
 
850
           IHIS2(N,I,J)=IHIS2(N,I,J)+1
 
851
           HIST2(N,I,J)=HIST2(N,I,J)+WGT
 
852
         ELSE
 
853
           IOSCORE2(N)=IOSCORE2(N)+1 
 
854
         ENDIF
 
855
      ENDIF
 
856
      END
 
857
 
 
858
 
 
859
      SUBROUTINE MOPERA2(I,OPER,J,K,X,Y)
 
860
      IMPLICIT NONE
 
861
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
862
      INTEGER I,J,K
 
863
      REAL*8 X,Y
 
864
      CHARACTER OPER*1
 
865
      INTEGER L1,L2
 
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
 
869
      DO L1=1,NBIN2(1,I)
 
870
      DO L2=1,NBIN2(2,I)
 
871
      IF(OPER.EQ.'+') THEN
 
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
 
879
          HIST2(K,L1,L2)=0.
 
880
          ELSE
 
881
          HIST2(K,L1,L2)=X*HIST2(I,L1,L2) / (Y*HIST2(J,L1,L2))
 
882
        ENDIF
 
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))
 
888
        ELSE
 
889
        HIST2(K,L1,L2)=0.
 
890
        ENDIF
 
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
 
895
             HIST2(K,L1,L2)=0.
 
896
             ELSE
 
897
             HIST2(K,L1,L2)=X*LOG10(Y*HIST2(I,L1,L2))
 
898
        ENDIF
 
899
      ELSE
 
900
      WRITE(98,5) OPER
 
901
   5  FORMAT(' ****** OPERATION ="',A1,'" UNKNOWN ********'/)
 
902
      RETURN
 
903
      ENDIF
 
904
      END DO
 
905
      ENDDO
 
906
      RETURN
 
907
  10  WRITE(98,20) I,J
 
908
  20  FORMAT(' ****** INCOMPATIBLE OPERATION HIST2 ',I2,' &',I2,
 
909
     &                                                   '*******'/)
 
910
      END
 
911
     
 
912
 
 
913
      SUBROUTINE MCOPY2(NIN,NOUT)
 
914
      IMPLICIT NONE
 
915
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
916
      INTEGER NIN,NOUT,I,J
 
917
C
 
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)
 
944
      DO I=1,NBIN2(1,NOUT)
 
945
        XHIS2(NOUT,I)=XHIS2(NIN,I)
 
946
        XPROJ(NOUT,I)=XPROJ(NIN,I)
 
947
      ENDDO
 
948
      DO I=1,NBIN2(2,NOUT)
 
949
        YHIS2(NOUT,I)=YHIS2(NIN,I)
 
950
        YPROJ(NOUT,I)=YPROJ(NIN,I)
 
951
      ENDDO
 
952
      DO I=1,NBIN2(1,NOUT)
 
953
        DO J=1,NBIN2(2,NOUT)
 
954
          HIST2(NOUT,I,J)=HIST2(NIN,I,J)
 
955
          IHIS2(NOUT,I,J)=IHIS2(NIN,I,J)
 
956
        ENDDO
 
957
      ENDDO
 
958
      END
 
959
 
 
960
 
 
961
      SUBROUTINE MFINAL2(N)
 
962
      IMPLICIT NONE
 
963
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
964
      INTEGER N,I,J
 
965
      REAL*8 XIN
 
966
      IF(BOOK2(N)(:1).NE.'Y') RETURN
 
967
      XIN=0.                                  
 
968
C-- projection on the x-axis
 
969
      DO 2 I=1,NBIN2(1,N)
 
970
        DO 1 J=1,NBIN2(2,N)
 
971
   1    XPROJ(N,I)=XPROJ(N,I)+HIST2(N,I,J)
 
972
   2  XIN=XIN+XPROJ(N,I)
 
973
      IF(XIN.EQ.0.) GO TO 10
 
974
C-- projection on the y-axis
 
975
      DO 3 J=1,NBIN2(2,N)
 
976
        DO 3 I=1,NBIN2(1,N)
 
977
   3    YPROJ(N,J)=YPROJ(N,J)+HIST2(N,I,J)
 
978
      HINT2(N)=XIN
 
979
      RETURN
 
980
  10  BOOK2(N)=' NO'
 
981
      END               
 
982
 
 
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)
 
989
      IMPLICIT NONE
 
990
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
991
      INTEGER N,I,J
 
992
      REAL*8 XIN
 
993
      IF(BOOK2(N)(1:1).NE.'Y') RETURN
 
994
      XIN=0.                                  
 
995
      DO I=1,NBIN2(1,N)
 
996
        XPROJ(N,I)=0.
 
997
      ENDDO
 
998
      DO J=1,NBIN2(2,N)
 
999
        YPROJ(N,J)=0.
 
1000
      ENDDO
 
1001
C-- projection on the x-axis
 
1002
      DO 2 I=1,NBIN2(1,N)
 
1003
        DO 1 J=1,NBIN2(2,N)
 
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
 
1008
      DO 3 J=1,NBIN2(2,N)
 
1009
        DO 3 I=1,NBIN2(1,N)
 
1010
   3    YPROJ(N,J)=YPROJ(N,J)+HIST2(N,I,J)
 
1011
      HINT2(N)=XIN
 
1012
      RETURN
 
1013
  10  BOOK2(N)=' NO'
 
1014
      END               
 
1015
 
 
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.
 
1019
      IMPLICIT NONE
 
1020
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
1021
      INTEGER N,M,IAX,I
 
1022
      BOOK(M)='YES'      
 
1023
      NBIN(M)=NBIN2(IAX,N)
 
1024
      HDEL(M)=HDEL2(IAX,N)
 
1025
      HMIN(M)=HMIN2(IAX,N)
 
1026
      HMAX(M)=HMAX2(IAX,N)
 
1027
      NHIST=MAX(NHIST,M)
 
1028
      TITLE(M)=TITLE2(N)//'(PROJ)'
 
1029
      DO I=1,NBIN(M)
 
1030
      IF(IAX.EQ.1)      THEN
 
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)
 
1036
      ENDIF
 
1037
      ENDDO
 
1038
      END                             
 
1039
 
 
1040
 
 
1041
      SUBROUTINE MNORM2(N,X)    
 
1042
      IMPLICIT NONE
 
1043
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
1044
      INTEGER N,I,J
 
1045
      REAL*8 X
 
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
 
1050
      HINT2(N)=X                      
 
1051
      END  
 
1052
 
 
1053
      SUBROUTINE MPRINT2(N)
 
1054
      IMPLICIT NONE
 
1055
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
1056
      INTEGER N,I,J
 
1057
      CHARACTER CTIME*10,CDATE*10
 
1058
      INTEGER INI,IMON,IDAY,IYEAR
 
1059
      DATA INI/0/               
 
1060
      IF(INI.EQ.0) THEN
 
1061
      CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
 
1062
      INI=1
 
1063
      ENDIF
 
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))
 
1068
     *                            ,I=1,NBIN2(1,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,//)
 
1074
      END
 
1075
 
 
1076
 
 
1077
      SUBROUTINE MTOP2(N,BTIT,LTIT,PLOPT)
 
1078
      IMPLICIT NONE
 
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
 
1086
      DATA INI/0/                  
 
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
 
1096
      IF(INI.EQ.0) THEN
 
1097
      ISEED=2143567+2*INT(SECNDS(0.0))
 
1098
      CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
 
1099
      INI=1
 
1100
      ENDIF
 
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)
 
1104
  100 FORMAT( /1x,                               
 
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,
 
1114
     &' PLOT AXES',/1X,
 
1115
     &' SET ORDER X Y ')
 
1116
 
 
1117
C-- squares with area proportional to the bin weight
 
1118
      IF(PLOPT.EQ.'BOX') THEN
 
1119
      SCMAX=0
 
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)))
 
1123
  1   CONTINUE
 
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))
 
1132
      ELSE
 
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))
 
1136
      ENDIF
 
1137
  2   CONTINUE
 
1138
 211  CONTINUE
 
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)
 
1141
      WRITE(99,300)
 
1142
  300 FORMAT('   PLOT')
 
1143
      WRITE(99,400) HINT2(N),IENT2(N)
 
1144
     &   ,IOSCORE2(N),IMON,IDAY,IYEAR,CTIME(1:5)
 
1145
  400 FORMAT( /1x,                               
 
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')
 
1154
      WRITE(99,500)
 
1155
  500 FORMAT('   NEW PLOT')
 
1156
 
 
1157
C-- dots - as many as the value of HINT2(n), distributed according to wgt
 
1158
      ELSEIF(PLOPT.EQ.'DOT') THEN                                        
 
1159
      XINT=0.
 
1160
      DO 11, I=1,NBIN2(1,N)
 
1161
      DO 11, J=1,NBIN2(2,N)
 
1162
      XINT=XINT+HIST2(N,I,J)    
 
1163
 11   CONTINUE
 
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)
 
1169
C      RN=(RANXX(ISEED))
 
1170
C      REST=HIST2(N,I,J)-FLOAT(NDOTS)
 
1171
C      IF(RN.LT.REST) NDOTS=NDOTS+1
 
1172
      DO ND=1,NDOTS
 
1173
      X=XHIS2(N,I)+(2.*RANXX(ISEED)-1.)*HDEL2(1,N)
 
1174
      Y=YHIS2(N,J)+(2.*RANXX(ISEED)-1.)*HDEL2(2,N)
 
1175
      WRITE(99,*) X,Y
 
1176
      ENDDO
 
1177
      NTOT=NTOT+NDOTS
 
1178
  3   CONTINUE       
 
1179
      WRITE(99,*) '  SET INTENSITY 1'
 
1180
      WRITE(99,300)                 
 
1181
      WRITE(99,600) HINT2(N),IENT2(N),NTOT
 
1182
     &   ,IMON,IDAY,IYEAR,CTIME(1:5)
 
1183
  600 FORMAT( /1x,                               
 
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')
 
1192
      WRITE(99,500)
 
1193
 
 
1194
      ENDIF
 
1195
      END
 
1196
 
 
1197
            
 
1198
      FUNCTION RANXX(SEED)
 
1199
*     -----------------
 
1200
* Ref.: K. Park and K.W. Miller, Comm. of the ACM 31 (1988) p.1192
 
1201
* Use seed = 1 as first value.
 
1202
*
 
1203
      IMPLICIT INTEGER(A-Z)
 
1204
      DOUBLE PRECISION MINV,RANXX
 
1205
      SAVE
 
1206
      PARAMETER(M=2147483647,A=16807,Q=127773,R=2836)
 
1207
      PARAMETER(MINV=0.46566128752458d-09)
 
1208
      HI = SEED/Q
 
1209
      LO = MOD(SEED,Q)
 
1210
      SEED = A*LO - R*HI
 
1211
      IF(SEED.LE.0) SEED = SEED + M
 
1212
      RANXX = SEED*MINV
 
1213
      END
 
1214
 
 
1215
 
 
1216
      SUBROUTINE MULTITOP(NH,NE,N,M,BTIT,LTIT,SCA)
 
1217
      IMPLICIT NONE
 
1218
      INCLUDE 'mcatnlo_hbook_gfortran8.inc'
 
1219
      INTEGER NH,NE,N,M
 
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,
 
1224
     # NOLD,NS
 
1225
      CHARACTER*(*) LTIT,BTIT,SCA
 
1226
      CHARACTER*7  PLOT(4)
 
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/
 
1234
C  DEFAULT SIZES
 
1235
      REAL*8 TIT0,LAB0,TIC0
 
1236
      DATA TIT0,LAB0,TIC0/3.,3.,0.06/
 
1237
      DATA INI/0/
 
1238
      IF(INI.EQ.0) THEN
 
1239
      CALL MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
 
1240
      IFRAME=0
 
1241
      WRITE(99,71) CDATE,CTIME(1:5)
 
1242
   71 FORMAT(4X,' ( ',A10,1X,A5/)
 
1243
      INI=1
 
1244
      ENDIF
 
1245
      IF(SCA.EQ.'REF') THEN
 
1246
        IFRAME=0 
 
1247
        RETURN
 
1248
      ENDIF
 
1249
      IF(BOOK(NH)(:1).NE.'Y') RETURN
 
1250
      IFRMAX=N*M
 
1251
      IFRAME=IFRAME+1
 
1252
      IF(IFRAME.GT.IFRMAX.OR.N.NE.NOLD.OR.M.NE.MOLD) THEN
 
1253
        IFRAME=1
 
1254
        WRITE(99,202)
 
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,'"')
 
1258
      ENDIF
 
1259
      IF(IFRAME.EQ.1) THEN
 
1260
        I=1
 
1261
        J=1
 
1262
      ELSEIF(IFRAME.LE.IFRMAX) THEN
 
1263
        IF(I.LE.N) I=I+1
 
1264
        IF(I.GT.N) THEN
 
1265
                I=1
 
1266
                J=J+1
 
1267
        ENDIF
 
1268
      ENDIF
 
1269
      IF(N.EQ.NOLD) GO TO 10
 
1270
      NS=N-1
 
1271
      XD=WIDTH/DFLOAT(N)
 
1272
      SRED=SQRT(DFLOAT(N*M))
 
1273
      TITS=TIT0/SRED
 
1274
      LABS=LAB0/SRED
 
1275
      TICS=TIC0/SRED
 
1276
      XTIT0=0.55*XPFRAC*XD
 
1277
      NOLD=N
 
1278
10    IF(M.EQ.MOLD) GO TO 20
 
1279
      YD=HEIGHT/DFLOAT(M)
 
1280
      YTIT0=0.06*YD
 
1281
      MOLD=M
 
1282
20    CONTINUE
 
1283
      XL=(I-1)*XD + XCORN
 
1284
      YL=(M-J)*YD + YCORN
 
1285
      XU=XL+XD*XPFRAC
 
1286
      YU=YL+YD*YPFRAC
 
1287
      IP=0
 
1288
c inizio modifiche
 
1289
      XMX=0.
 
1290
      DO IBIN=1,NBIN(NH)
 
1291
        X=HIST(NH,IBIN)
 
1292
        IF(X.NE.0.) THEN
 
1293
           IF(XMX.EQ.0.) THEN
 
1294
              FMX = X + HIST(NE,IBIN)
 
1295
              FMN = X - HIST(NE,IBIN)
 
1296
           ELSE
 
1297
              FMX=MAX(FMX,X + HIST(NE,IBIN))
 
1298
              FMN=MIN(FMN,X - HIST(NE,IBIN))
 
1299
           ENDIF
 
1300
        ENDIF
 
1301
        XMX=MAX(XMX,ABS(X)+ HIST(NE,IBIN))
 
1302
      ENDDO
 
1303
      SCALE=SCA
 
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
 
1310
cRF
 
1311
c$$$        IF(FMN.LE.0.) THEN
 
1312
c$$$                SCALE='LIN'
 
1313
c$$$                GO TO 50
 
1314
c$$$        ENDIF
 
1315
        FMAX=10.**( AINT(LOG10(ABS(FMX))+1000001) - 1000000 )
 
1316
        FMIN=10.**( AINT(LOG10(ABS(FMN))+1000000) - 1000000 )
 
1317
      ENDIF
 
1318
C fine modifiche
 
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)
 
1326
      XTIT=XL+XTIT0
 
1327
      YTIT=YU+YTIT0
 
1328
      WRITE(99,101) XL,YTIT,TITLE(NH)(1:40)
 
1329
101   FORMAT('  TITLE ',2(F8.4,1X),'"',A,'"')
 
1330
      YTIT=YTIT-2.*YTIT0
 
1331
      WRITE(99,102) XTIT,YTIT,HINT(NH)
 
1332
102   FORMAT('  TITLE ',2(F8.4,1X),'" INT=',1PE10.3,'"')
 
1333
      YTIT=YTIT-YTIT0
 
1334
      WRITE(99,103) XTIT,YTIT,IENT(NH)
 
1335
103   FORMAT('  TITLE ',2(F8.4,1X),'" ENT=',I9,'"')
 
1336
      YTIT=YTIT-YTIT0
 
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,'"')
 
1340
        YTIT=YTIT-YTIT0
 
1341
      ENDIF
 
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,'"')
 
1345
        YTIT=YTIT-YTIT0
 
1346
      ENDIF
 
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),/,
 
1349
     *       2X,'JOIN TEXT',/,
 
1350
     *       2X,2(F8.4,1X),/,2(F8.4,1X),/,
 
1351
     *       2X,'JOIN TEXT')
 
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
 
1356
107   FORMAT(
 
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')
 
1364
C
 
1365
C  END HEADER , FILL TOPDRAWER WITH DATA
 
1366
C
 
1367
      ENTRY MTFILL(NH,NE,N,M,BTIT,LTIT,SCA)
 
1368
      IP=IP+1
 
1369
      IF(IP.GT.4) IP=1
 
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)
 
1375
200   CONTINUE
 
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',/,/)
 
1380
203   RETURN
 
1381
      END
 
1382
 
 
1383
      SUBROUTINE NEWPLOT
 
1384
      WRITE(99,202) 
 
1385
202   FORMAT('   NEW PLOT',/,/)
 
1386
      CALL MULTITOP(1,1,1,1,' ',' ','REF')
 
1387
      END                         
 
1388
 
 
1389
C*******************************************************************
 
1390
C     END OF THE HISTOGRAMMING PACKAGE
 
1391
C*******************************************************************
 
1392
 
 
1393
 
 
1394
C***** GENERIC UTILITIES *********************************
 
1395
C
 
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
 
1403
 
 
1404
      SUBROUTINE VSUM(P1,P2,Q)
 
1405
      IMPLICIT NONE
 
1406
      REAL*8 P1(4),P2(4),Q(4)
 
1407
      INTEGER I
 
1408
      DO I=1,4
 
1409
      Q(I)=P1(I)+P2(I)
 
1410
      END DO
 
1411
      END
 
1412
 
 
1413
      SUBROUTINE RNDINT(N,INDX)
 
1414
C   returns a random ordering of the string of integers (1,2,..,N)
 
1415
      IMPLICIT NONE
 
1416
      INTEGER N
 
1417
      INTEGER INDX(N),ISEED,I
 
1418
      COMMON/RNDM/ISEED
 
1419
      DOUBLE PRECISION RANXX
 
1420
      REAL*8 X(500)
 
1421
      IF(N.GT.500) WRITE(6,*) 'WARNING**RNDINT'
 
1422
      DO I=1,N
 
1423
      X(I)=RANXX(ISEED)
 
1424
      ENDDO
 
1425
      CALL UTSORT(X,N,INDX,2)
 
1426
      END
 
1427
 
 
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------------------------------------------------------------------------
 
1433
      IMPLICIT NONE
 
1434
      INTEGER N,IOPT
 
1435
      REAL*8 A(N),B(500)
 
1436
      INTEGER K(N),IL(500),IR(500),I,J
 
1437
      IF (N.GT.500) WRITE(6,*) 'WARNING**UTSORT'
 
1438
      IL(1)=0
 
1439
      IR(1)=0
 
1440
      DO 10 I=2,N
 
1441
      IL(I)=0
 
1442
      IR(I)=0
 
1443
      J=1
 
1444
   2  IF(A(I).GT.A(J)) GO TO 5
 
1445
   3  IF(IL(J).EQ.0) GO TO 4
 
1446
      J=IL(J)
 
1447
      GO TO 2
 
1448
   4  IR(I)=-J
 
1449
      IL(J)=I
 
1450
      GO TO 10
 
1451
   5  IF(IR(J).LE.0) GO TO 6
 
1452
      J=IR(J)
 
1453
      GO TO 2
 
1454
   6  IR(I)=IR(J)
 
1455
      IR(J)=I
 
1456
  10  CONTINUE
 
1457
      I=1
 
1458
      J=1
 
1459
      GO TO 8
 
1460
  20  J=IL(J)
 
1461
   8  IF(IL(J).GT.0) GO TO 20
 
1462
   9  K(I)=J
 
1463
      B(I)=A(J)
 
1464
      I=I+1
 
1465
      IF(IR(J)) 12,30,13
 
1466
  13  J=IR(J)
 
1467
      GO TO 8
 
1468
  12  J=-IR(J)
 
1469
      GO TO 9
 
1470
  30  IF(IOPT.EQ.2) RETURN
 
1471
      DO 31 I=1,N
 
1472
  31  A(I)=B(I)
 
1473
 999  END
 
1474
 
 
1475
      SUBROUTINE MBWARN(ROUT)
 
1476
      CHARACTER*(*) ROUT
 
1477
      WRITE(*,*) '***********************************************'
 
1478
      WRITE(*,*) '***** WARNING CALLED FROM ROUTINE ',ROUT,':'
 
1479
      END
 
1480
 
 
1481
 
 
1482
      SUBROUTINE MYDATE(IMON,IDAY,IYEAR,CTIME,CDATE)
 
1483
      IMPLICIT NONE
 
1484
      INTEGER IMON,IDAY,IYEAR,IVAL(3)
 
1485
      CHARACTER*10 CTIME,CDATE
 
1486
      CHARACTER*10 CTIME2,CDATE2
 
1487
      CALL IDATE(IVAL)
 
1488
      IDAY=IVAL(1)
 
1489
      IMON=IVAL(2)
 
1490
      IYEAR=IVAL(3)
 
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)
 
1495
      END