Millepede-II  V04-00-00_preview
 All Classes Files Functions Variables Enumerator Pages
mphistab.f90
Go to the documentation of this file.
1 
2 ! Code converted using TO_F90 by Alan Miller
3 ! Date: 2012-03-16 Time: 11:07:45
4 
70 
71 ! *************************** Histograms ******************************
72 
73 SUBROUTINE hmpdef(ih,xa,xb,text) ! book, reset histogram
74  IMPLICIT NONE
75  INTEGER:: i
76  INTEGER:: iha
77  INTEGER:: ihb
78  INTEGER:: ihc
79  INTEGER:: ix
80  INTEGER:: j
81  INTEGER:: lun
82  INTEGER:: lunw
83  INTEGER:: nbin
84  INTEGER:: nn
85  REAL:: x
86  REAL:: xcent
87  REAL:: xmean
88  REAL:: xsigm
89  ! book millepede histogram, 120 bins
90 
91  INTEGER, INTENT(IN) :: ih
92  REAL, INTENT(IN) :: xa
93  REAL, INTENT(IN) :: xb
94  CHARACTER (LEN=*), INTENT(IN) :: text
95  INTEGER, PARAMETER :: numhis=15
96  INTEGER:: inhist(120,numhis)
97  INTEGER::jnhist(5,numhis)
98  INTEGER::khist(numhis)
99  REAL:: fnhist(120,numhis)
100  equivalence(inhist(1,1),fnhist(1,1))
101  INTEGER:: kvers(numhis)
102  REAL:: xl(6,numhis)
103  DOUBLE PRECISION:: dl(2,numhis)
104  CHARACTER (LEN=60):: htext(numhis)
105  SAVE
106  DATA khist/numhis*0/,lun/7/
107  ! ...
108  IF(ih <= 0.OR.ih > numhis) return
109  ! IF(XA.EQ.XB) RETURN
110  DO i=1,120
111  inhist(i,ih)=0
112  END DO
113  DO j=1,5
114  jnhist(j,ih)=0
115  END DO
116  xl(1,ih)=xa
117  xl(2,ih)=xb
118  xl(3,ih)=0.0
119  IF(xa /= xb) xl(3,ih)=120.0/(xb-xa)
120  xl(6,ih)=0.5*(xa+xb) ! center
121  IF(khist(ih) == 0) THEN
122  kvers(ih)=0
123  ELSE
124  kvers(ih)=kvers(ih)+1
125  END IF
126  khist(ih)=1 ! flt.pt. (lin)
127  htext(ih)=text
128  dl(1,ih)=0.0d0
129  dl(2,ih)=0.0d0
130  return
131 
132  entry hmpldf(ih,text) ! book, reset log histogram
133  IF(ih <= 0.OR.ih > numhis) return
134  DO i=1,120
135  inhist(i,ih)=0
136  END DO
137  DO j=1,5
138  jnhist(j,ih)=0
139  END DO
140  IF(khist(ih) == 0) THEN
141  kvers(ih)=0
142  ELSE
143  kvers(ih)=kvers(ih)+1
144  END IF
145  khist(ih)=2 ! integer log
146  htext(ih)=text
147  xl(1,ih)=0.0
148  xl(2,ih)=6.0
149  return
150 
151  entry hmpent(ih,x) ! entry flt.pt.
152  IF(ih <= 0.OR.ih > numhis) return
153  IF(khist(ih) /= 1) return
154  IF(jnhist(4,ih) >= 2147483647) return
155  jnhist(4,ih)=jnhist(4,ih)+1 ! count
156  IF(jnhist(4,ih) <= 120) THEN
157  fnhist(jnhist(4,ih),ih)=x ! store value
158  IF(jnhist(4,ih) == 120) THEN
159  CALL hmpmak(inhist(1,ih),fnhist(1,ih),jnhist(1,ih), xl(1,ih),dl(1,ih))
160  END IF
161  return
162  END IF
163  ! IF(JNHIST(1,IH)+JNHIST(2,IH)+JNHIST(3,IH).EQ.0) THEN
164  ! XL(4,IH)=X
165  ! XL(5,IH)=X
166  ! END IF
167  i=int(1.0+xl(3,ih)*(x-xl(1,ih))) ! X - Xmin
168  j=2
169  IF(i < 1) j=1
170  IF(i > 120) j=3
171  jnhist(j,ih)=jnhist(j,ih)+1
172  xl(4,ih)=min(xl(4,ih),x)
173  xl(5,ih)=max(xl(5,ih),x)
174  IF(j /= 2) return
175  inhist(i,ih)=inhist(i,ih)+1
176  dl(1,ih)=dl(1,ih)+ x-xl(6,ih)
177  dl(2,ih)=dl(2,ih)+(x-xl(6,ih))**2
178  return
179 
180  entry hmplnt(ih,ix) ! entry integer
181  IF(ih <= 0.OR.ih > numhis) return
182  IF(khist(ih) /= 2) return
183  IF(jnhist(1,ih) >= 2147483647) return
184  IF(ix <= 0) THEN
185  jnhist(1,ih)=jnhist(1,ih)+1
186  ELSE
187  IF(jnhist(4,ih) == 0) jnhist(4,ih)=ix
188  IF(jnhist(5,ih) == 0) jnhist(5,ih)=ix
189  jnhist(4,ih)=min(jnhist(4,ih),ix)
190  jnhist(5,ih)=max(jnhist(5,ih),ix)
191  i=int(1.0+20.0*log10(REAL(ix)))
192  j=2
193  IF(i < 1) j=1
194  IF(i > 120) j=3
195  IF(j == 2) inhist(i,ih)=inhist(i,ih)+1
196  jnhist(j,ih)=jnhist(j,ih)+1
197  END IF
198  return
199 
200  entry hmprnt(ih) ! print, content vert
201  IF(ih == 0) THEN
202  iha=1
203  ihb=numhis
204  ELSE
205  IF(ih <= 0.OR.ih > numhis) return
206  iha=ih
207  ihb=ih
208  END IF
209  DO ihc=iha,ihb
210  IF(khist(ihc) /= 0) THEN
211  IF(khist(ihc) == 1) THEN
212  CALL hmpmak(inhist(1,ihc),fnhist(1,ihc),jnhist(1,ihc), &
213  xl(1,ihc),dl(1,ihc))
214  END IF
215  nn=jnhist(1,ihc)+jnhist(2,ihc)+jnhist(3,ihc)
216  IF(nn /= 0.OR.khist(ihc) == 3) THEN
217  WRITE(*,111)
218 111 format(' ______',2('______________________________'))
219  IF(kvers(ihc) == 1) THEN
220  WRITE(*,*) 'Histogram',ihc,': ',htext(ihc)
221  ELSE
222  WRITE(*,*) 'Histogram',ihc,'/',kvers(ihc),': ',htext(ihc)
223  END IF
224  IF(khist(ihc) == 1) THEN
225  WRITE(*,*) ' Out_low inside out_high = ', (jnhist(j,ihc),j=1,3)
226  ELSE IF(khist(ihc) == 2) THEN
227  WRITE(*,*) ' 0_or_negative inside above_10^6 = ', &
228  (jnhist(j,ihc),j=1,3)
229  END IF
230  IF(khist(ihc) == 3) THEN
231  CALL pfvert(120,fnhist(1,ihc))
232  END IF
233  IF(jnhist(2,ihc) /= 0) THEN ! integer content
234  CALL pivert(120,inhist(1,ihc))
235  IF(khist(ihc) == 1) THEN
236  CALL psvert(xl(1,ihc),xl(2,ihc))
237  ELSE IF(khist(ihc) == 2) THEN
238  CALL psvert(0.0,6.0)
239  END IF
240  END IF
241  IF(khist(ihc) == 1) THEN
242  WRITE(*,*) ' Min and Max are',xl(4,ihc),xl(5,ihc)
243  IF(jnhist(2,ihc) > 1) THEN
244  xmean=REAL(xl(6,ihc)+dl(1,ihc)/float(jnhist(2,ihc)))
245  xcent=0.5*(xl(1,ihc)+xl(2,ihc))
246  xsigm=REAL((dl(2,ihc)-dl(1,ihc)**2/float(jnhist(2,ihc))))
247  xsigm=sqrt(xsigm/float(jnhist(2,ihc)-1))
248  WRITE(*,*) ' Mean and sigma are', xmean,' +-',xsigm
249  END IF
250  ELSE IF(khist(ihc) == 2) THEN
251  WRITE(*,*) ' Plot of log10 of entries. Min and Max are', &
252  jnhist(4,ihc),jnhist(5,ihc)
253  END IF
254  END IF
255  END IF
256  END DO
257  return
258 
259  entry hmplun(lunw) ! unit for output
260  lun=lunw
261  return
262 
263  entry hmpwrt(ih) ! write histogram text file
264  IF(lun <= 0) return
265  IF(ih == 0) THEN
266  iha=1
267  ihb=numhis
268  ELSE
269  IF(ih <= 0.OR.ih > numhis) return
270  iha=ih
271  ihb=ih
272  END IF
273 
274  DO ihc=iha,ihb ! histogram loop
275  IF(khist(ihc) /= 0) THEN
276  IF(khist(ihc) == 1) THEN
277  CALL hmpmak(inhist(1,ihc),fnhist(1,ihc),jnhist(1,ihc), &
278  xl(1,ihc),dl(1,ihc))
279  END IF
280  nbin=120
281  WRITE(lun,204) ' '
282  WRITE(lun,201) ihc,kvers(ihc),khist(ihc)
283  WRITE(lun,204) htext(ihc)
284  IF (jnhist(1,ihc)+jnhist(2,ihc)+jnhist(3,ihc) == 0 &
285  .AND.xl(1,ihc) == xl(2,ihc)) THEN
286  ! hist is empty and hist range makes no sense
287  ! - cause: hist with 'variable edges' was never filled
288  ! - workaround: make lower and upper edge of hist differ in output
289  WRITE(lun,202) nbin,xl(1,ihc)-0.001,xl(2,ihc)+0.001
290  ELSE
291  WRITE(lun,202) nbin,xl(1,ihc),xl(2,ihc)
292  END IF
293  WRITE(lun,203) (jnhist(j,ihc),j=1,3)
294  WRITE(lun,204) 'bincontent'
295  IF(khist(ihc) == 1.OR.khist(ihc) == 2) THEN
296  CALL kprint(lun,inhist(1,ihc),nbin)
297  ELSE
298  WRITE(lun,219) (fnhist(i,ihc),i=1,nbin)
299  END IF
300 
301  IF(khist(ihc) == 1) THEN
302  WRITE(lun,205) xl(4,ihc),xl(5,ihc)
303  ELSE IF(khist(ihc) == 2) THEN
304  WRITE(lun,205) float(jnhist(4,ihc)),float(jnhist(5,ihc))
305  END IF
306  IF(khist(ihc) == 1) THEN
307  IF(jnhist(2,ihc) > 1) THEN
308  xmean=REAL(xl(6,ihc)+dl(1,ihc)/float(jnhist(2,ihc)))
309  xcent=0.5*(xl(1,ihc)+xl(2,ihc))
310  xsigm=REAL((dl(2,ihc)-dl(1,ihc)**2/float(jnhist(2,ihc))))
311  xsigm=sqrt(xsigm/float(jnhist(2,ihc)-1))
312  WRITE(lun,206) xmean,xsigm
313  END IF
314  END IF
315  WRITE(lun,204) 'end of histogram'
316  END IF
317  END DO
318 
319 201 format('Histogram ',i4,10x,'version ',i4,10x,'type',i2)
320 202 format(10x,' bins, limits ',i4,2g15.5)
321 203 format(10x,'out-low inside out-high ',3i10)
322 204 format(a)
323 205 format('minmax',2e15.7)
324 206 format('meansigma',2e15.7)
325 
326 219 format(4e15.7)
327 END SUBROUTINE hmpdef
328 
329 SUBROUTINE hmpmak(inhist,fnhist,jnhist,xl,dl) ! hist scale from data
330  IMPLICIT NONE
331  INTEGER:: i
332  INTEGER:: j
333  INTEGER:: k
334  INTEGER:: nn
335  REAL:: x
336  REAL:: xa
337  REAL:: xb
338 
339  INTEGER, INTENT(OUT) :: inhist(120)
340  REAL, INTENT(IN) :: fnhist(120)
341  INTEGER, INTENT(IN OUT) :: jnhist(5)
342  REAL, INTENT(IN OUT) :: xl(6)
343  DOUBLE PRECISION, INTENT(OUT) :: dl(2)
344  REAL:: cphist(120)
345 
346 
347 
348  SAVE
349  ! ...
350  nn=jnhist(4)
351  ! WRITE(*,*) 'HMPMAK: NN,JNHIST(5)',NN,JNHIST(5)
352  IF(nn == 0.OR.jnhist(5) /= 0) return
353  jnhist(5)=1
354  DO i=1,nn
355  ! WRITE(*,*) 'copy ',I,FNHIST(I)
356  cphist(i)=fnhist(i)
357  END DO
358  CALL heapf(cphist,nn)
359  IF(xl(3) == 0.0) THEN
360  CALL bintab(cphist,nn,xa,xb)
361  xl(1)=xa
362  xl(2)=xb
363  xl(3)=0.0
364  IF(xa /= xb) xl(3)=120.0/(xb-xa)
365  xl(6)=0.5*(xa+xb) ! center
366  END IF
367  xl(4)=cphist( 1)
368  xl(5)=cphist(nn)
369  ! WRITE(*,*) 'XL ',XL
370  DO i=1,nn
371  inhist(i)=0
372  END DO
373  DO k=1,nn
374  x=cphist(k)
375  i=int(1.0+xl(3)*(x-xl(1))) ! X - Xmin
376  ! WRITE(*,*) 'K,I,X ',K,I,X
377  j=2
378  IF(i < 1) j=1
379  IF(i > 120) j=3
380  jnhist(j)=jnhist(j)+1
381  IF(j == 2) THEN
382  inhist(i)=inhist(i)+1
383  dl(1)=dl(1)+ x-xl(6)
384  dl(2)=dl(2)+(x-xl(6))**2
385  END IF
386  END DO
387 END SUBROUTINE hmpmak
388 
389 
390 
391 
392 SUBROUTINE bintab(tab,n,xa,xb) ! hist scale from data
393  IMPLICIT NONE
394  REAL:: dd
395  REAL:: dx
396  INTEGER:: i
397  INTEGER:: iexp
398  INTEGER:: ii
399  INTEGER:: j
400  INTEGER:: m1
401  INTEGER:: m2
402  INTEGER:: n1
403  INTEGER:: n2
404  REAL:: rat
405  REAL:: x1
406  REAL:: x2
407  REAL:: xx
408  ! Bin limits XA and XB from TAB(N)
409 
410  REAL, INTENT(IN) :: tab(n)
411  INTEGER, INTENT(IN) :: n
412  REAL, INTENT(OUT) :: xa
413  REAL, INTENT(OUT) :: xb
414 
415  REAL:: bin(10)
416  DATA bin/1.0,1.5,2.0,3.0,4.0,5.0,8.0,10.0,15.0,20.0/
417  SAVE
418  ! ...
419 
420  CALL heapf(tab,n) ! reduced statistic
421  ! WRITE(*,*) ' '
422  ! WRITE(*,*) 'Sorted ',(TAB(I),I=1,N)
423  IF(n < 100) THEN
424  x1=tab(1)
425  x2=tab(n)
426  ! WRITE(*,*) 'reduced statistic X1 X2 ',X1,X2
427  ELSE ! large statistic
428  m1=int(1.0+0.05*float(n))
429  m2=int(1.0+0.16*float(n))
430  x1=tab(m1)-4.0*(tab(m2)-tab(m1))
431  IF(x1 < 0.0.AND.tab(1) >= 0.0) x1=tab(1)
432  x2=tab(n+1-m1)+4.0*(tab(n+1-m1)-tab(n+1-m2))
433  IF(x2 > 0.0.AND.tab(n) <= 0.0) x2=tab(n)
434  ! WRITE(*,*) 'large statistic ',X1,X2
435  ! WRITE(*,*) 'min und max ',TAB(1),TAB(N)
436  IF(x1*tab(1) <= 0.0) x1=0.0
437  IF(x2*tab(n) <= 0.0) x2=0.0
438  ! WRITE(*,*) 'large statistic zero ',X1,X2
439  IF(x1*x2 < 0.0.AND.min(-x1,x2) > 0.6*max(-x1,x2)) THEN
440  xx=max(-x1,x2) ! symmetry
441  x1=-xx
442  x2=+xx
443  ELSE IF(x1*x2 > 0.0.AND. & ! include zero ?
444  abs(min(x1,x2)) < 0.4*abs(max(x1,x2))) THEN
445  IF(x1 < 0.0) THEN
446  x2=0.0
447  ELSE
448  x1=0.0
449  END IF
450  END IF
451  ! WRITE(*,*) 'large statistic ',X1,X2
452  END IF
453  IF(x1 == x2) THEN
454  x1=x1-1.0
455  x2=x2+1.0
456  END IF
457  dx=x2-x1
458  ! WRITE(*,*) 'X1,X2,DX ',X1,X2,DX
459  rat=0.0
460  ii=1
461  DO j=1,11
462  i=j
463  IF(j == 11) i=ii
464  iexp=int(101.0+log10(dx)-log10(6.0*bin(i)))
465  iexp=iexp-100
466  dd=bin(i)*10.0**iexp
467 
468  n1=int(abs(x1)/dd)
469  IF(x1 < 0.0) n1=-n1
470  IF(float(n1)*dd > x1) n1=n1-1
471  ! WRITE(*,*) 'Bin ',I,N1,N1*DD,X1
472 
473  n2=int(abs(x2)/dd)
474  IF(x2 < 0.0) n2=-n2
475  IF(float(n2)*dd < x2) n2=n2+1
476  ! WRITE(*,*) 'Bin ',I,N2,N2*DD,X2
477 10 IF(n2-n1 < 6) THEN
478  IF(n1 /= 0) n1=n1-1
479  IF(n2-n1 < 6.AND.n2 /= 0) n2=n2+1
480  go to 10
481  END IF
482  ! WRITE(*,*) 'corrected N1 N2 ',N1,N2
483  xa=sign(float(n1)*dd,x1)
484  xb=sign(float(n2)*dd,x2)
485  ! WRITE(*,*) J,' resulting limits XA XB ',XA,XB
486  IF((x2-x1)/(xb-xa) > rat) THEN
487  ii=i
488  rat=(x2-x1)/(xb-xa)
489  END IF
490  END DO
491 ! WRITE(*,*) J,' resulting limits XA XB ',XA,XB
492 END SUBROUTINE bintab
493 
494 SUBROUTINE kprint(lun,list,n) ! print integer array
495  IMPLICIT NONE
496  INTEGER:: i
497  INTEGER:: ia
498  INTEGER:: ib
499  INTEGER:: k
500  INTEGER:: ln
501  INTEGER:: lp
502  INTEGER:: np
503  ! print integer array LIST(N)
504 
505  INTEGER, INTENT(IN OUT) :: lun
506  INTEGER, INTENT(IN) :: list(n)
507  INTEGER, INTENT(IN) :: n
508  INTEGER:: li(7)
509  DATA li/2,3,4,6,8,9,12/ ! number of characters
510  SAVE
511  ! ...
512  ib=0
513 10 ia=ib+1
514  IF(ia > n) return
515  DO k=1,7
516  np=72/li(k)
517  ib=min(ia-1+np,n)
518  IF(k <= 6) THEN
519  lp=10**(li(k)-1)-1 ! maximum positive
520  ln=-lp/10 ! minimum negative
521  DO i=ia,ib
522  IF(list(i) > lp.OR.list(i) < ln) go to 20
523  END DO
524  END IF
525  IF(k == 1) THEN
526  WRITE(lun,101) (list(i),i=ia,ib)
527  ELSE IF(k == 2) THEN
528  WRITE(lun,102) (list(i),i=ia,ib)
529  ELSE IF(k == 3) THEN
530  WRITE(lun,103) (list(i),i=ia,ib)
531  ELSE IF(k == 4) THEN
532  WRITE(lun,104) (list(i),i=ia,ib)
533  ELSE IF(k == 5) THEN
534  WRITE(lun,105) (list(i),i=ia,ib)
535  ELSE IF(k == 6) THEN
536  WRITE(lun,106) (list(i),i=ia,ib)
537  ELSE IF(k == 7) THEN
538  WRITE(lun,107) (list(i),i=ia,ib)
539  END IF
540  go to 10
541 20 continue
542  END DO
543 101 format(36i2)
544 102 format(24i3)
545 103 format(18i4)
546 104 format(12i6)
547 105 format( 9i8)
548 106 format( 8i9)
549 107 format( 6i12)
550 END SUBROUTINE kprint
551 
552 ! ***************************** XY data ****************************
553 
554 SUBROUTINE gmpdef(ig,ityp,text) ! book, reset XY storage
555  IMPLICIT NONE
556  REAL:: dx
557  REAL:: dy
558  INTEGER:: i
559  INTEGER:: iga
560  INTEGER:: igb
561  INTEGER:: igc
562  INTEGER:: j
563  INTEGER:: lun
564  INTEGER:: lunw
565  INTEGER:: n
566  INTEGER:: na
567  REAL:: wght
568  REAL:: x
569  REAL:: y
570  REAL:: y1
571  ! ITYP = 1 X,Y as dots
572  ! = 2 X,Y as line
573  ! = 3 X,Y as line and dots
574  ! = 4 X,Y, DX,DY symbols
575 
576  INTEGER, INTENT(IN) :: ig
577  INTEGER, INTENT(IN) :: ityp
578  CHARACTER (LEN=*), INTENT(IN) :: text
579  INTEGER, PARAMETER :: narr=1000
580  REAL:: array(2,narr)
581  REAL::array4(4,narr/2)
582  REAL::array1(narr+narr)
583  REAL::four(4)
584  equivalence(array(1,1),array4(1,1),array1(1))
585  INTEGER, PARAMETER :: numgxy=10
586  INTEGER, PARAMETER :: nlimit=500
587  INTEGER:: nstr(numgxy)
588  INTEGER::igtp(numgxy)
589  INTEGER::lvers(numgxy)
590  INTEGER::nst(3,numgxy)
591  REAL:: xyplws(10,numgxy)
592  INTEGER:: jflc(5,numgxy)
593  INTEGER::kflc(5,numgxy)
594  ! JFLC(1,.) = first used index
595  ! JFLC(2,.) = last used index
596  ! JFLC(3,.) = counter of used places
597  ! JFLC(4,.) = counter of ignored
598  ! JFLC(5,.) = limit for JFLC(3)
599  CHARACTER (LEN=60):: gtext(numgxy)
600 
601  LOGICAL:: start
602  SAVE
603  DATA start/.true./,lun/7/
604  DATA nstr/numgxy*0/ ! by GF
605  ! ...
606  IF(start) THEN
607  start=.false.
608  CALL stmars ! initialize storage
609  DO i=1,numgxy
610  DO j=1,5
611  jflc(j,i)=0
612  kflc(j,i)=0
613  END DO
614  END DO
615  END IF
616 
617  IF(ig < 1.OR.ig > numgxy) return
618  IF(ityp < 1.OR.ityp > 5) return
619  IF(nstr(ig) == 0) THEN
620  lvers(ig)=0
621  ELSE
622  lvers(ig)=lvers(ig)+1
623  END IF
624  nstr(ig)=1 ! by GF
625  ! remove stored elements
626  IF(jflc(1,ig) /= 0) CALL stmarm(jflc(1,ig))
627  IF(kflc(1,ig) /= 0) CALL stmarm(kflc(1,ig))
628  igtp(ig)=ityp
629  gtext(ig)=text
630  DO j=1,5
631  jflc(j,ig)=0
632  END DO
633  jflc(5,ig)=nlimit
634  IF(ityp == 5) THEN
635  DO j=1,5
636  kflc(j,ig)=0
637  END DO
638  jflc(5,ig)=128 ! maximum of 128 values
639  kflc(5,ig)=narr
640  nst(1,ig)=0
641  nst(2,ig)=0
642  nst(3,ig)=1
643  DO j=1,10
644  xyplws(j,ig)=0.0
645  END DO
646  END IF
647  return
648 
649  entry gmpxy(ig,x,y) ! add (X,Y) pair
650  IF(ig < 1.OR.ig > numgxy) return ! check argument IG
651  IF(igtp(ig) < 1.OR.igtp(ig) > 3) return ! check type
652  CALL stmapr(jflc(1,ig),x,y)
653  return
654 
655  entry gmpxyd(ig,x,y,dx,dy) ! add (X,Y,DX,DY)
656  IF(ig < 1.OR.ig > numgxy) return ! check argument IG
657  IF(igtp(ig) /= 4) return
658  four(1)=x
659  four(2)=y
660  four(3)=dx
661  four(4)=dy
662  CALL stmadp(jflc(1,ig),four)
663  return
664 
665  entry gmpms(ig,x,y) ! mean sigma(X) from Y
666  ! mean sigma from Y, as a function of X
667  ! WRITE(*,*) 'GMPMS ',IG,X,Y
668 
669  IF(ig < 1.OR.ig > numgxy) return ! check argument IG
670  IF(igtp(ig) /= 5) return
671 
672  xyplws(10,ig)=x ! last X coordinate
673  IF(nst(1,ig) == 0) THEN
674  y1=y
675  nst(1,ig)=1
676  IF(kflc(3,ig) == 0) xyplws(9,ig)=x ! start coordinate
677  ELSE
678  nst(1,ig)=0
679  CALL stmapr(kflc(1,ig),y1,y) ! store pair
680  IF(kflc(3,ig) >= kflc(5,ig)) THEN
681  CALL stmacp(kflc(1,ig),array,n) ! get data
682  CALL stmarm(kflc(1,ig)) ! remove data
683  n=n+n
684  CALL rmesig(array,n,xyplws(2,ig),xyplws(4,ig))
685  nst(2,ig)=nst(2,ig)+1
686  IF(nst(2,ig) == 1) xyplws(7,ig)=xyplws(9,ig)
687  xyplws(8,ig)=x ! end coordinate
688  xyplws(5,ig)=xyplws(5,ig)+xyplws(2,ig)
689  xyplws(6,ig)=xyplws(6,ig)+xyplws(4,ig)
690  IF(nst(2,ig) == nst(3,ig)) THEN
691  xyplws(1,ig)=0.5*(xyplws(7,ig)+xyplws(8,ig))
692  xyplws(2,ig)=xyplws(5,ig)/float(nst(3,ig))
693  xyplws(3,ig)=0.5*(xyplws(8,ig)-xyplws(7,ig))
694  xyplws(4,ig)=xyplws(6,ig)/float(nst(3,ig))
695  xyplws(5,ig)=0.0
696  xyplws(6,ig)=0.0
697  nst(2,ig)=0
698  CALL stmadp(jflc(1,ig),xyplws(1,ig))
699  IF(jflc(3,ig) >= jflc(5,ig)) THEN
700  CALL stmacp(jflc(1,ig),array4,n) ! get data
701  n=n/2
702  CALL stmarm(jflc(1,ig)) ! remove data
703  DO i=1,n,2 ! average
704  xyplws(7,ig)=array4(1,i )-array4(3, i)
705  xyplws(8,ig)=array4(1,i+1)+array4(3,i+1)
706  xyplws(1,ig)=0.5*(xyplws(7,ig)+xyplws(8,ig))
707  xyplws(2,ig)=0.5*(array4(2,i)+array4(2,i+1))
708  xyplws(3,ig)=0.5*(xyplws(8,ig)-xyplws(7,ig))
709  xyplws(4,ig)=0.5*(array4(4,i)+array4(4,i+1))
710  CALL stmadp(jflc(1,ig),xyplws(1,ig))
711  END DO
712  nst(3,ig)=nst(3,ig)+nst(3,ig)
713  END IF
714  END IF
715  END IF
716  END IF
717  return
718 
719  entry gmprnt(ig) ! print XY data
720  IF(ig == 0) THEN
721  iga=1
722  igb=numgxy
723  ELSE
724  IF(ig <= 0.OR.ig > numgxy) return
725  iga=ig
726  igb=ig
727  END IF
728  DO igc=iga,igb
729 
730  IF(igtp(igc) >= 1.AND.igtp(igc) <= 3) THEN
731  WRITE(*,*) ' '
732  WRITE(*,*) 'Store ',igc,': ',gtext(igc)
733  IF(jflc(4,igc) == 0) THEN
734  WRITE(*,*) ' stored n-tuples: ',jflc(3,igc)
735  ELSE
736  WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', &
737  jflc(3,igc),', ',jflc(4,igc)
738  END IF
739 
740  CALL stmacp(jflc(1,igc),array,na) ! get all data
741 
742  DO n=1,na
743  WRITE(*,102) n, array(1,n),array(2,n)
744  END DO
745 
746  ELSE IF(igtp(igc) == 4) THEN
747 
748  WRITE(*,*) ' '
749  WRITE(*,*) 'Store ',igc,': ',gtext(igc)
750  IF(jflc(4,igc) == 0) THEN
751  WRITE(*,*) ' stored n-tuples: ',jflc(3,igc)
752  ELSE
753  WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', &
754  jflc(3,igc),', ',jflc(4,igc)
755  END IF
756 
757  CALL stmacp(jflc(1,igc),array,na) ! get all data
758  na=na/2
759 
760  DO n=1,na
761  WRITE(*,102) n,(array4(j,n),j=1,4)
762  END DO
763 
764  ELSE IF(igtp(igc) == 5) THEN
765 
766  CALL stmacp(kflc(1,igc),array,n) ! get data
767  CALL stmarm(kflc(1,igc)) ! remove data
768  n=n+n
769  IF(nst(1,igc) == 1) THEN
770  n=n+1
771  array1(n)=y1
772  nst(1,igc)=0 ! reset
773  END IF
774  IF(n /= 0) THEN
775  xyplws(7,igc)=xyplws( 9,igc)
776  xyplws(8,igc)=xyplws(10,igc)
777  CALL rmesig(array1,n,xyplws(2,igc),xyplws(4,igc))
778  wght=float(n)/float(nst(3,igc)*kflc(5,igc))
779  xyplws(5,igc)=xyplws(5,igc)+xyplws(2,igc)*wght
780  xyplws(6,igc)=xyplws(6,igc)+xyplws(4,igc)*wght
781  xyplws(2,igc)=xyplws(5,igc)/(float(nst(2,igc))+wght)
782  xyplws(4,igc)=xyplws(6,igc)/(float(nst(2,igc))+wght)
783  xyplws(1,igc)=0.5*(xyplws(7,igc)+xyplws(8,igc))
784  xyplws(3,igc)=0.5*(xyplws(8,igc)-xyplws(7,igc))
785  CALL stmadp(jflc(1,igc),xyplws(1,igc))
786  END IF
787 
788  WRITE(*,*) ' '
789  WRITE(*,*) 'Store ',igc,': ',gtext(igc)
790  IF(jflc(4,igc) == 0) THEN
791  WRITE(*,*) ' stored n-tuples: ',jflc(3,igc)
792  ELSE
793  WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', &
794  jflc(3,igc),', ',jflc(4,igc)
795  END IF
796 
797  CALL stmacp(jflc(1,igc),array,na) ! get all data
798  na=na/2
799  DO n=1,na
800  WRITE(*,102) n,(array4(j,n),j=1,4)
801  END DO
802  END IF
803  END DO
804  return
805 
806  entry gmplun(lunw) ! unit for output
807  lun=lunw
808  return
809 
810  entry gmpwrt(ig) ! write XY text file
811  IF(lun <= 0) return
812  IF(ig == 0) THEN
813  iga=1
814  igb=numgxy
815  ELSE
816  IF(ig <= 0.OR.ig > numgxy) return
817  iga=ig
818  igb=ig
819  END IF
820  DO igc=iga,igb
821  IF(igtp(igc) == 5) THEN
822 
823  CALL stmacp(kflc(1,igc),array,n) ! get data
824  CALL stmarm(kflc(1,igc)) ! remove data
825  n=n+n
826  IF(nst(1,igc) == 1) THEN
827  n=n+1
828  array1(n)=y1
829  nst(1,igc)=0 ! reset
830  END IF
831  IF(n /= 0) THEN
832  xyplws(7,igc)=xyplws( 9,igc)
833  xyplws(8,igc)=xyplws(10,igc)
834  CALL rmesig(array1,n,xyplws(2,igc),xyplws(4,igc))
835  wght=float(n)/float(nst(3,igc)*kflc(5,igc))
836  xyplws(5,igc)=xyplws(5,igc)+xyplws(2,igc)*wght
837  xyplws(6,igc)=xyplws(6,igc)+xyplws(4,igc)*wght
838  xyplws(2,igc)=xyplws(5,igc)/(float(nst(2,igc))+wght)
839  xyplws(4,igc)=xyplws(6,igc)/(float(nst(2,igc))+wght)
840  xyplws(1,igc)=0.5*(xyplws(7,igc)+xyplws(8,igc))
841  xyplws(3,igc)=0.5*(xyplws(8,igc)-xyplws(7,igc))
842  CALL stmadp(jflc(1,igc),xyplws(1,igc))
843  END IF
844 
845  END IF
846  IF(jflc(3,igc)+jflc(4,igc) /= 0) THEN
847  WRITE(lun,204) ' '
848  WRITE(lun,201) igc,lvers(igc),igtp(igc)
849  WRITE(lun,204) gtext(igc)
850  WRITE(lun,203) jflc(3,igc)+jflc(4,igc)
851  CALL stmacp(jflc(1,igc),array,na) ! get all data
852  IF(igtp(igc) >= 1.AND.igtp(igc) <= 3) THEN
853  WRITE(lun,204) 'x-y'
854  DO n=1,na
855  WRITE(lun,205) array(1,n),array(2,n)
856  END DO
857  ELSE IF(igtp(igc) == 4.OR.igtp(igc) == 5) THEN
858  WRITE(lun,204) 'x-y-dx-dy'
859  na=na/2
860  DO n=1,na
861  WRITE(lun,205) (array4(j,n),j=1,4)
862  END DO
863  END IF
864  WRITE(lun,204) 'end of xy-data'
865  END IF
866  END DO
867 102 format(i12,4g15.7)
868  ! 103 FORMAT(' Index ___X___ ___Y___ '/
869  ! + ' ----- -------------- --------------')
870  ! 104 FORMAT(' Index ___X___ ___Y___ ',
871  ! + ' ___DX__ ___DY__ '/
872  ! + ' ----- -------------- --------------',
873  ! + ' -------------- --------------')
874 201 format('XY-Data ',i4,10x,'version ',i4,10x,'type',i2)
875 203 format(10x,'stored not-stored ',2i10)
876 204 format(a)
877 205 format(3x,4g15.7)
878 END SUBROUTINE gmpdef
879 
880 
881 SUBROUTINE stmars ! init/reset storage
882  IMPLICIT NONE
883  INTEGER:: i
884  INTEGER:: ifre
885  INTEGER:: ifrea
886  INTEGER:: ifreb
887  INTEGER:: ind
888  INTEGER:: j
889  INTEGER:: n
890  INTEGER:: ndim
891  REAL:: x
892  REAL:: y
893  parameter(ndim=5000) ! storage dimension, should be NUMGXY*NLIMIT
894  REAL:: tk(2,ndim) ! pair storage for data pairs
895  INTEGER:: next(ndim) ! pointer
896  INTEGER:: iflc1 ! first and last index of free pairs
897  INTEGER::iflc2 ! first and last index of free pairs
898  SAVE
899 
900  REAL:: four(4) ! double_pair, copy array
901  REAL::array(2,*) ! double_pair, copy array
902  INTEGER:: jflc(5) ! user array
903  ! JFLC(1) = first used index
904  ! JFLC(2) = last used index
905  ! JFLC(3) = counter of used places
906  ! JFLC(4) = counter of ignored
907  ! JFLC(5) = limit for JFLC(3)
908  ! ...
909  DO i=1,ndim
910  next(i)=i+1 ! pointer to next free location
911  tk(1,i)=0.0 ! reset
912  tk(2,i)=0.0
913  END DO
914  next(ndim)=0 ! ... and end pointer
915  iflc1=1 ! index first free pair
916  iflc2=ndim ! index last free pair
917  return
918 
919  entry stmapr(jflc,x,y) ! store pair (X,Y)
920  ifre=iflc1 ! index of free place
921  IF(ifre == 0.OR.jflc(3) >= jflc(5)) THEN ! overflow
922  jflc(4)=jflc(4)+1
923  ELSE
924  iflc1=next(ifre) ! pointer to new free location
925  IF(jflc(1) == 0) THEN ! first item
926  jflc(1)=ifre
927  ELSE
928  next(jflc(2))=ifre
929  END IF
930  next(ifre)=0
931  jflc(2)=ifre ! last index
932  jflc(3)=jflc(3)+1 ! counter
933  tk(1,ifre)=x
934  tk(2,ifre)=y
935  END IF
936  return
937 
938  entry stmadp(jflc,four) ! store double pair
939  ifrea=iflc1 ! index of 1. free place
940  IF(ifrea == 0) THEN ! overflow
941  jflc(4)=jflc(4)+1
942  ELSE
943  ifreb=next(iflc1) ! index of 2. free place
944  IF(ifreb == 0.OR.jflc(3) >= 2*jflc(5)) THEN ! overflow
945  jflc(4)=jflc(4)+1
946  ELSE
947  iflc1=next(ifreb) ! pointer to new free location
948  IF(jflc(1) == 0) THEN ! first item
949  jflc(1)=ifrea
950  ELSE
951  next(jflc(2))=ifrea
952  END IF
953  next(ifreb)=0
954  jflc(2)=ifreb ! last index
955  jflc(3)=jflc(3)+1 ! counter
956  tk(1,ifrea)=four(1)
957  tk(2,ifrea)=four(2)
958  tk(1,ifreb)=four(3)
959  tk(2,ifreb)=four(4)
960  END IF
961  END IF
962  return
963 
964  entry stmacp(jflc,array,n) ! copy (cp) all pairs to array
965  n=0
966  ind=jflc(1)
967 10 IF(ind == 0) return
968  n=n+1
969  array(1,n)=tk(1,ind)
970  array(2,n)=tk(2,ind)
971  ind=next(ind)
972  go to 10
973 
974  entry stmarm(jflc) ! remove (rm) stored paiirs
975  next(iflc2)=jflc(1) ! connect to free space
976  iflc2=jflc(2) ! new last free index
977  DO j=1,4
978  jflc(j)=0
979  END DO
980 END SUBROUTINE stmars ! init/
981 
982 SUBROUTINE rmesig(x,n,xloc,xsca) ! robust mean and sigma
983  IMPLICIT NONE
984  INTEGER:: i
985  ! robust determination of location and scale parameter,
986  ! for Gaussian data: location=mean and scale=standard deviation
987  ! XLOC = median of X_i (N values in array X(N))
988  ! XCSA = median of | X_i - XLOC |, times 1.4826
989 
990  REAL, INTENT(IN OUT) :: x(n) ! input array, modified
991  INTEGER, INTENT(IN) :: n
992  REAL, INTENT(OUT) :: xloc
993  REAL, INTENT(OUT) :: xsca
994  SAVE
995  ! ...
996  xloc=0.0
997  xsca=0.0
998  IF(n <= 0) return
999  CALL heapf(x,n) ! sort
1000  xloc=0.5*(x((n+1)/2)+x((n+2)/2)) ! location
1001  DO i=1,n
1002  x(i)=abs(x(i)-xloc)
1003  END DO
1004  CALL heapf(x,n) ! sort
1005  xsca=1.4826*0.5*(x((n+1)/2)+x((n+2)/2)) ! dispersion
1006 END SUBROUTINE rmesig
1007 
1008 
1009 
Generated on Wed May 23 2012 16:53:05 for Millepede-II by   doxygen 1.8.1