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