HIJING/hijing1_36/vegas.F

   1 * $Id$
   2 C*******************************************************************
   3 C
   4 C
   5 C
   6 C
   7 C*******************************************************************
   8 C   SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N
   9 C      - BY G.P. LEPAGE   SEPT 1976/(REV)APR 1978
  10 C*******************************************************************
  11 C
  12       SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A)
  13       IMPLICIT REAL*8(A-H,O-Z)
  14 #include "bveg1_hijing.inc"
  15 #include "bveg2_hijing.inc"
  16 #include "bveg3_hijing.inc"
  17       EXTERNAL FXN
  18       DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10)
  19      1   ,KG(10),IA(10)
  20       REAL*4 QRAN(10)
  21       DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/-1/
  22       SAVE
  23 C
  24       NDO=1
  25       DO 1 J=1,NDIM
  26 1     XI(1,J)=ONE
  27 C
  28       ENTRY VEGAS1(FXN,AVGI,SD,CHI2A)
  29 C         - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
  30       IT=0
  31       SI=0.
  32       SI2=SI
  33       SWGT=SI
  34       SCHI=SI
  35 C
  36       ENTRY VEGAS2(FXN,AVGI,SD,CHI2A)
  37 C         - NO INITIALIZATION
  38       ND=NDMX
  39       NG=1
  40       IF(MDS.EQ.0) GO TO 2
  41       NG=(NCALL/2.)**(1./NDIM)
  42       MDS=1
  43       IF((2*NG-NDMX).LT.0) GO TO 2
  44       MDS=-1
  45       NPG=NG/NDMX+1
  46       ND=NG/NPG
  47       NG=NPG*ND
  48 2     K=NG**NDIM
  49       NPG=NCALL/K
  50       IF(NPG.LT.2) NPG=2
  51       CALLS=NPG*K
  52       DXG=ONE/NG
  53       DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE)
  54       XND=ND
  55       NDM=ND-1
  56       DXG=DXG*XND
  57       XJAC=ONE/CALLS
  58       DO 3 J=1,NDIM
  59 c***this is the line 50
  60       DX(J)=XU(J)-XL(J)
  61 3     XJAC=XJAC*DX(J)
  62 C
  63 C   REBIN PRESERVING BIN DENSITY
  64 C
  65       IF(ND.EQ.NDO) GO TO 8
  66       RC=NDO/XND
  67       DO 7 J=1,NDIM
  68       K=0
  69       XN=0.
  70       DR=XN
  71       I=K
  72 4     K=K+1
  73       DR=DR+ONE
  74       XO=XN
  75       XN=XI(K,J)
  76 5     IF(RC.GT.DR) GO TO 4
  77       I=I+1
  78       DR=DR-RC
  79       XIN(I)=XN-(XN-XO)*DR
  80       IF(I.LT.NDM) GO TO 5
  81       DO 6 I=1,NDM
  82 6     XI(I,J)=XIN(I)
  83 7     XI(ND,J)=ONE
  84       NDO=ND
  85 C
  86 8     CONTINUE
  87       IF(NPRN.NE.0) WRITE(16,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND
  88      1                           ,(XL(J),XU(J),J=1,NDIM)
  89 C
  90       ENTRY VEGAS3(FXN,AVGI,SD,CHI2A)
  91 C         - MAIN INTEGRATION LOOP
  92 9     IT=IT+1
  93       TI=0.
  94       TSI=TI
  95       DO 10 J=1,NDIM
  96       KG(J)=1
  97       DO 10 I=1,ND
  98       D(I,J)=TI
  99 10    DI(I,J)=TI
 100 C
 101 11    FB=0.
 102       F2B=FB
 103       K=0
 104 12    K=K+1
 105       CALL ARAN9(QRAN,NDIM)
 106       WGT=XJAC
 107       DO 15 J=1,NDIM
 108       XN=(KG(J)-QRAN(J))*DXG+ONE
 109 c*****this is the line 100
 110       IA(J)=XN
 111       IF(IA(J).GT.1) GO TO 13
 112       XO=XI(IA(J),J)
 113       RC=(XN-IA(J))*XO
 114       GO TO 14
 115 13    XO=XI(IA(J),J)-XI(IA(J)-1,J)
 116       RC=XI(IA(J)-1,J)+(XN-IA(J))*XO
 117 14    X(J)=XL(J)+RC*DX(J)
 118       WGT=WGT*XO*XND
 119 15    CONTINUE
 120 C
 121       F=WGT
 122       F=F*FXN(X,WGT)
 123       F2=F*F
 124       FB=FB+F
 125       F2B=F2B+F2
 126       DO 16 J=1,NDIM
 127       DI(IA(J),J)=DI(IA(J),J)+F
 128 16    IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2
 129       IF(K.LT.NPG) GO TO 12
 130 C
 131       F2B=DSQRT(F2B*NPG)
 132       F2B=(F2B-FB)*(F2B+FB)
 133       TI=TI+FB
 134       TSI=TSI+F2B
 135       IF(MDS.GE.0) GO TO 18
 136       DO 17 J=1,NDIM
 137 17    D(IA(J),J)=D(IA(J),J)+F2B
 138 18    K=NDIM
 139 19    KG(K)=MOD(KG(K),NG)+1
 140       IF(KG(K).NE.1) GO TO 11
 141       K=K-1
 142       IF(K.GT.0) GO TO 19
 143 C
 144 C   FINAL RESULTS FOR THIS ITERATION
 145 C
 146       TSI=TSI*DV2G
 147       TI2=TI*TI
 148       WGT=TI2/(TSI+1.0d-37)
 149       SI=SI+TI*WGT
 150       SI2=SI2+TI2
 151       SWGT=SWGT+WGT
 152       SWGT=SWGT+1.0D-37
 153       SI2=SI2+1.0D-37
 154       SCHI=SCHI+TI2*WGT
 155       AVGI=SI/(SWGT)
 156       SD=SWGT*IT/(SI2)
 157       CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(IT-.999)
 158       SD=DSQRT(ONE/SD)
 159 C****this is the line 150
 160       IF(NPRN.EQ.0) GO TO 21
 161       TSI=DSQRT(TSI)
 162       WRITE(16,201) IT,TI,TSI,AVGI,SD,CHI2A
 163       IF(NPRN.GE.0) GO TO 21
 164       DO 20 J=1,NDIM
 165 20    WRITE(16,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND)
 166 C
 167 C   REFINE GRID
 168 C
 169 21    DO 23 J=1,NDIM
 170       XO=D(1,J)
 171       XN=D(2,J)
 172       D(1,J)=(XO+XN)/2.
 173       DT(J)=D(1,J)
 174       DO 22 I=2,NDM
 175       D(I,J)=XO+XN
 176       XO=XN
 177       XN=D(I+1,J)
 178       D(I,J)=(D(I,J)+XN)/3.
 179 22    DT(J)=DT(J)+D(I,J)
 180       D(ND,J)=(XN+XO)/2.
 181 23    DT(J)=DT(J)+D(ND,J)
 182 C
 183       DO 28 J=1,NDIM
 184       RC=0.
 185       DO 24 I=1,ND
 186       R(I)=0.
 187       IF (DT(J).GE.1.0D18) THEN
 188        WRITE(6,*) '************** A SINGULARITY >1.0D18'
 189 C      WRITE(5,1111)
 190 C1111  FORMAT(1X,'**************IMPORTANT NOTICE***************')
 191 C      WRITE(5,1112)
 192 C1112  FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18')
 193 C      WRITE(5,1113)
 194 C1113  FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS')
 195 C      WRITE(5,1114)
 196 C1114  FORMAT(1X,'**************END NOTICE*************')
 197       END IF
 198       IF(D(I,J).LE.1.0D-18) GO TO 24
 199       XO=DT(J)/D(I,J)
 200       R(I)=((XO-ONE)/XO/DLOG(XO))**ALPH
 201 24    RC=RC+R(I)
 202       RC=RC/XND
 203       K=0
 204       XN=0.
 205       DR=XN
 206       I=K
 207 25    K=K+1
 208       DR=DR+R(K)
 209       XO=XN
 210 c****this is the line 200
 211       XN=XI(K,J)
 212 26    IF(RC.GT.DR) GO TO 25
 213       I=I+1
 214       DR=DR-RC
 215       XIN(I)=XN-(XN-XO)*DR/(R(K)+1.0d-30)
 216       IF(I.LT.NDM) GO TO 26
 217       DO 27 I=1,NDM
 218 27    XI(I,J)=XIN(I)
 219 28    XI(ND,J)=ONE
 220 C
 221       IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9
 222 200   FORMAT('0INPUT PARAMETERS FOR VEGAS:  NDIM=',I3,'  NCALL=',F8.0
 223      1    /28X,'  IT=',I5,'  ITMX=',I5/28X,'  ACC=',G9.3
 224      2    /28X,'  MDS=',I3,'   ND=',I4/28X,'  (XL,XU)=',
 225      3    (T40,'( ',G12.6,' , ',G12.6,' )'))
 226 201   FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3,
 227      1    ':   INTEGRAL =',G14.8/21X,'STD DEV  =',G10.4 /
 228      2    ' ACCUMULATED RESULTS:   INTEGRAL =',G14.8 /
 229      3    24X,'STD DEV  =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4)
 230 202   FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,'  DELT I  ',
 231      1    2X,' CONV''CE  ',11X,'X',7X,'  DELT I  ',2X,' CONV''CE  '
 232      2   ,11X,'X',7X,'  DELT I  ',2X,' CONV''CE  ' /
 233      2    (' ',3G12.4,5X,3G12.4,5X,3G12.4))
 234       RETURN
 235       END