This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gheisha / phasp.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:01  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.39  by  S.Giani
11 *-- Author :
12       SUBROUTINE PHASP
13 C
14 C *** NVE 29-MAR-1988 CERN GENEVA ***
15 C
16 C CALLED BY : NUCREC TWOCLU GENXPT
17 C ORIGIN : H.FESEFELDT (02-DEC-1986)
18 C
19 #include "geant321/s_prntfl.inc"
20 #include "geant321/s_genio.inc"
21 #include "geant321/limits.inc"
22 C
23 #if !defined(CERNLIB_SINGLE)
24       DOUBLE PRECISION WTMAX,WTMAXQ,WTFC,TWGT,ONE,TEXPXL,TEXPXU
25       PARAMETER (ONE=1.D0)
26 #endif
27 #if defined(CERNLIB_SINGLE)
28       PARAMETER (ONE=1.0)
29 #endif
30       LOGICAL LZERO
31       DIMENSION EMM(18)
32       DIMENSION RNO(50)
33       DIMENSION EM(18),PD(18),EMS(18),SM(18),FFQ(18),PCM1(90)
34       EQUIVALENCE (NT,NPG),(AMASS(1),EM(1)),(PCM1(1),PCM(1,1))
35       SAVE KNT
36 C
37       DATA  FFQ/0.,3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
38      $                       256.3704, 268.4705, 240.9780, 189.2637,
39      $                       132.1308,  83.0202,  47.4210,  24.8295,
40      $                        12.0006,   5.3858,   2.2560,   0.8859/
41       DATA  KNT , TWOPI /  1 , 6.2831853073 /
42 C
43 C --- Initialise local arrays and the result array PCM ---
44       DO 10 JZERO=1,18
45          PCM(1,JZERO)=0.
46          PCM(2,JZERO)=0.
47          PCM(3,JZERO)=0.
48          PCM(4,JZERO)=0.
49          PCM(5,JZERO)=0.
50          EMM(JZERO)  =0.
51          PD(JZERO)   =0.
52          EMS(JZERO)  =0.
53          SM(JZERO)   =0.
54   10  CONTINUE
55 C
56       KNT = KNT + 1
57       IF (.NOT.NPRT(3).AND..NOT.NPRT(4)) GOTO 100
58       WRITE(NEWBCD,1200) NPG,TECM,(AMASS(JK),JK=1,NPG)
59   100 CONTINUE
60   150 IF (NT .LT. 2)  GO TO 1001
61       IF (NT .GT. 18)  GO TO 1002
62       NTM1=NT-1
63       NTM2=NT-2
64       NTNM4 = 3*NT - 4
65       EMM(1)=EM(1)
66       TM=0.0
67       DO 200 I=1,NT
68       EMS(I)=EM(I)**2
69       TM=TM+EM(I)
70  200  SM(I)=TM
71       WGT=1.
72  210  TECMTM=TECM-TM
73       IF (TECMTM .LE. 0.0)  GO TO 1000
74       EMM(NT)=TECM
75       IF (KGENEV.GT.1) GO TO 400
76       EMMAX=TECMTM+EM(1)
77       EMMIN=0.0
78 C
79 C          For weight calculation, form sum of log's of terms
80 C          instead of product of terms. Note that thereby WTMAX
81 C          and WTMAXQ are changed in their contents; they are
82 C          currently not used outside the range from here to
83 C          label 531. We also need to check for zero factors now.
84 C          Negative values cannot appear as GPDK always returns a
85 C          nonnegative number. As coded, even the exotic cases
86 C          NT<2 (first loop not executed) and NTM1<1 (second loop
87 C          not executed) should be safe and give the same result
88 C          for WTG in the end as the old code.
89 C
90       WTMAX=0.0
91       LZERO=.TRUE.
92       DO 350 I=2,NT
93       EMMIN=EMMIN+EM(I-1)
94       EMMAX=EMMAX+EM(I)
95       WTFC=GPDK(EMMAX,EMMIN,EM(I))
96       IF(WTFC.LE.0.) THEN
97       LZERO=.FALSE.
98       GOTO 351
99       ENDIF
100       WTMAX=WTMAX+LOG(WTFC)
101  350  CONTINUE
102  351  WTMAXQ= EXPXU
103       IF(LZERO) WTMAXQ= -WTMAX
104       GO TO 455
105   400 WTMAXQ=LOG(ONE*TECMTM**NTM2*FFQ(NT) / TECM)
106   455 CONTINUE
107       CALL GRNDM(RNO,NTNM4)
108       IF(NTM2) 900,509,460
109   460 CONTINUE
110       CALL FLPSOR(RNO,NTM2)
111       DO 508 J=2,NTM1
112   508 EMM(J)=RNO(J-1)*TECMTM+SM(J)
113   509 TWGT=WTMAXQ
114       IR=NTM2
115       LZERO=.TRUE.
116       DO 530 I=1,NTM1
117       PD(I)=GPDK(EMM(I+1),EMM(I),EM(I+1))
118       IF(PD(I).LE.0.0) THEN
119       LZERO=.FALSE.
120       ELSE
121       TWGT=TWGT+LOG(ONE*PD(I))
122       ENDIF
123   530 CONTINUE
124   531 WGT=0.0
125       IF(LZERO) THEN
126       TEXPXU=EXPXU
127       TEXPXL=EXPXL
128       WGT=EXP(MAX(MIN(TWGT,TEXPXU),TEXPXL))
129       ENDIF
130       PCM(1,1)=0.0
131       PCM(2,1)=PD(1)
132       PCM(3,1)=0.0
133       DO 570 I=2,NT
134       PCM(1,I)=0.0
135       PCM(2,I) = -PD(I-1)
136       PCM(3,I)=0.0
137       IR=IR+1
138       BANG=TWOPI*RNO(IR)
139       CB=COS(BANG)
140       SB=SIN(BANG)
141       IR=IR+1
142       C=2.0*RNO(IR)-1.0
143       S=SQRT(ABS(1.0-C*C))
144       IF(I.EQ.NT) GO TO 1567
145       ESYS=SQRT(PD(I)**2+EMM(I)**2)
146       BETA=PD(I)/ESYS
147       GAMA=ESYS/EMM(I)
148       DO 568 J=1,I
149       NDX = 5*J - 5
150       AA= PCM1(NDX+1)**2 + PCM1(NDX+2)**2 + PCM1(NDX+3)**2
151       PCM1(NDX+5) = SQRT(AA)
152       PCM1(NDX+4) = SQRT(AA+EMS(J))
153       CALL ROTES2(C,S,CB,SB,PCM,J)
154       PSAVE = GAMA*(PCM(2,J)+BETA*PCM(4,J))
155   568 PCM(2,J)=PSAVE
156       GO TO 570
157  1567 DO 1568 J=1,I
158       AA=PCM(1,J)**2 + PCM(2,J)**2 + PCM(3,J)**2
159       PCM(5,J)=SQRT(AA)
160       PCM(4,J)=SQRT(AA+EMS(J))
161       CALL ROTES2(C,S,CB,SB,PCM,J)
162  1568 CONTINUE
163   570 CONTINUE
164   900 CONTINUE
165       RETURN
166  1000 DO 212 I=1,NPG
167       PCM(1,I)=0.
168       PCM(2,I)=0.
169       PCM(3,I)=0.
170       PCM(4,I)=AMASS(I)
171   212 PCM(5,I)=AMASS(I)
172       WGT=0.
173       RETURN
174  1001 IF(NPRT(3).OR.NPRT(4)) WRITE(NEWBCD,1101)
175       GO TO 1050
176  1002 WRITE(NEWBCD,1102)
177  1050 WRITE(NEWBCD,1150) KNT
178       WRITE(NEWBCD,1200) NPG,TECM,(AMASS(JK),JK=1,NPG)
179       RETURN
180  1100 FORMAT (1H0,'AVAILABLE ENERGY NEGATIVE')
181  1101 FORMAT (1H0,'LESS THAN 2 OUTGOING PARTICLES')
182  1102 FORMAT (1H0,'MORE THAN 18 OUTGOING PARTICLES')
183  1150 FORMAT (1H0,'ABOVE ERROR DETECTED IN PHASP AT CALL NUMBER',I7)
184  1200 FORMAT (1H0,'INPUT DATA TO PHASP.         NPG= ' ,I6/
185      +2X,9H   TECM=  ,D15.7   ,18H  PARTICLE MASSES=,5D15.7/(42X,5D15.7)
186      +)
187       END