[u/mrichter/AliRoot.git] / GEANT321 / gheisha / phasp.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:01  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.39  by  S.Giani
*-- Author :
      SUBROUTINE PHASP
C
C *** NVE 29-MAR-1988 CERN GENEVA ***
C
C CALLED BY : NUCREC TWOCLU GENXPT
C ORIGIN : H.FESEFELDT (02-DEC-1986)
C
#include "geant321/s_prntfl.inc"
#include "geant321/s_genio.inc"
#include "geant321/limits.inc"
C
#if !defined(CERNLIB_SINGLE)
      DOUBLE PRECISION WTMAX,WTMAXQ,WTFC,TWGT,ONE,TEXPXL,TEXPXU
      PARAMETER (ONE=1.D0)
#endif
#if defined(CERNLIB_SINGLE)
      PARAMETER (ONE=1.0)
#endif
      LOGICAL LZERO
      DIMENSION EMM(18)
      DIMENSION RNO(50)
      DIMENSION EM(18),PD(18),EMS(18),SM(18),FFQ(18),PCM1(90)
      EQUIVALENCE (NT,NPG),(AMASS(1),EM(1)),(PCM1(1),PCM(1,1))
      SAVE KNT
C
      DATA  FFQ/0.,3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
     $                       256.3704, 268.4705, 240.9780, 189.2637,
     $                       132.1308,  83.0202,  47.4210,  24.8295,
     $                        12.0006,   5.3858,   2.2560,   0.8859/
      DATA  KNT , TWOPI /  1 , 6.2831853073 /
C
C --- Initialise local arrays and the result array PCM ---
      DO 10 JZERO=1,18
         PCM(1,JZERO)=0.
         PCM(2,JZERO)=0.
         PCM(3,JZERO)=0.
         PCM(4,JZERO)=0.
         PCM(5,JZERO)=0.
         EMM(JZERO)  =0.
         PD(JZERO)   =0.
         EMS(JZERO)  =0.
         SM(JZERO)   =0.
  10  CONTINUE
C
      KNT = KNT + 1
      IF (.NOT.NPRT(3).AND..NOT.NPRT(4)) GOTO 100
      WRITE(NEWBCD,1200) NPG,TECM,(AMASS(JK),JK=1,NPG)
  100 CONTINUE
  150 IF (NT .LT. 2)  GO TO 1001
      IF (NT .GT. 18)  GO TO 1002
      NTM1=NT-1
      NTM2=NT-2
      NTNM4 = 3*NT - 4
      EMM(1)=EM(1)
      TM=0.0
      DO 200 I=1,NT
      EMS(I)=EM(I)**2
      TM=TM+EM(I)
 200  SM(I)=TM
      WGT=1.
 210  TECMTM=TECM-TM
      IF (TECMTM .LE. 0.0)  GO TO 1000
      EMM(NT)=TECM
      IF (KGENEV.GT.1) GO TO 400
      EMMAX=TECMTM+EM(1)
      EMMIN=0.0
C
C          For weight calculation, form sum of log's of terms
C          instead of product of terms. Note that thereby WTMAX
C          and WTMAXQ are changed in their contents; they are
C          currently not used outside the range from here to
C          label 531. We also need to check for zero factors now.
C          Negative values cannot appear as GPDK always returns a
C          nonnegative number. As coded, even the exotic cases
C          NT<2 (first loop not executed) and NTM1<1 (second loop
C          not executed) should be safe and give the same result
C          for WTG in the end as the old code.
C
      WTMAX=0.0
      LZERO=.TRUE.
      DO 350 I=2,NT
      EMMIN=EMMIN+EM(I-1)
      EMMAX=EMMAX+EM(I)
      WTFC=GPDK(EMMAX,EMMIN,EM(I))
      IF(WTFC.LE.0.) THEN
      LZERO=.FALSE.
      GOTO 351
      ENDIF
      WTMAX=WTMAX+LOG(WTFC)
 350  CONTINUE
 351  WTMAXQ= EXPXU
      IF(LZERO) WTMAXQ= -WTMAX
      GO TO 455
  400 WTMAXQ=LOG(ONE*TECMTM**NTM2*FFQ(NT) / TECM)
  455 CONTINUE
      CALL GRNDM(RNO,NTNM4)
      IF(NTM2) 900,509,460
  460 CONTINUE
      CALL FLPSOR(RNO,NTM2)
      DO 508 J=2,NTM1
  508 EMM(J)=RNO(J-1)*TECMTM+SM(J)
  509 TWGT=WTMAXQ
      IR=NTM2
      LZERO=.TRUE.
      DO 530 I=1,NTM1
      PD(I)=GPDK(EMM(I+1),EMM(I),EM(I+1))
      IF(PD(I).LE.0.0) THEN
      LZERO=.FALSE.
      ELSE
      TWGT=TWGT+LOG(ONE*PD(I))
      ENDIF
  530 CONTINUE
  531 WGT=0.0
      IF(LZERO) THEN
      TEXPXU=EXPXU
      TEXPXL=EXPXL
      WGT=EXP(MAX(MIN(TWGT,TEXPXU),TEXPXL))
      ENDIF
      PCM(1,1)=0.0
      PCM(2,1)=PD(1)
      PCM(3,1)=0.0
      DO 570 I=2,NT
      PCM(1,I)=0.0
      PCM(2,I) = -PD(I-1)
      PCM(3,I)=0.0
      IR=IR+1
      BANG=TWOPI*RNO(IR)
      CB=COS(BANG)
      SB=SIN(BANG)
      IR=IR+1
      C=2.0*RNO(IR)-1.0
      S=SQRT(ABS(1.0-C*C))
      IF(I.EQ.NT) GO TO 1567
      ESYS=SQRT(PD(I)**2+EMM(I)**2)
      BETA=PD(I)/ESYS
      GAMA=ESYS/EMM(I)
      DO 568 J=1,I
      NDX = 5*J - 5
      AA= PCM1(NDX+1)**2 + PCM1(NDX+2)**2 + PCM1(NDX+3)**2
      PCM1(NDX+5) = SQRT(AA)
      PCM1(NDX+4) = SQRT(AA+EMS(J))
      CALL ROTES2(C,S,CB,SB,PCM,J)
      PSAVE = GAMA*(PCM(2,J)+BETA*PCM(4,J))
  568 PCM(2,J)=PSAVE
      GO TO 570
 1567 DO 1568 J=1,I
      AA=PCM(1,J)**2 + PCM(2,J)**2 + PCM(3,J)**2
      PCM(5,J)=SQRT(AA)
      PCM(4,J)=SQRT(AA+EMS(J))
      CALL ROTES2(C,S,CB,SB,PCM,J)
 1568 CONTINUE
  570 CONTINUE
  900 CONTINUE
      RETURN
 1000 DO 212 I=1,NPG
      PCM(1,I)=0.
      PCM(2,I)=0.
      PCM(3,I)=0.
      PCM(4,I)=AMASS(I)
  212 PCM(5,I)=AMASS(I)
      WGT=0.
      RETURN
 1001 IF(NPRT(3).OR.NPRT(4)) WRITE(NEWBCD,1101)
      GO TO 1050
 1002 WRITE(NEWBCD,1102)
 1050 WRITE(NEWBCD,1150) KNT
      WRITE(NEWBCD,1200) NPG,TECM,(AMASS(JK),JK=1,NPG)
      RETURN
 1100 FORMAT (1H0,'AVAILABLE ENERGY NEGATIVE')
 1101 FORMAT (1H0,'LESS THAN 2 OUTGOING PARTICLES')
 1102 FORMAT (1H0,'MORE THAN 18 OUTGOING PARTICLES')
 1150 FORMAT (1H0,'ABOVE ERROR DETECTED IN PHASP AT CALL NUMBER',I7)
 1200 FORMAT (1H0,'INPUT DATA TO PHASP.         NPG= ' ,I6/
     +2X,9H   TECM=  ,D15.7   ,18H  PARTICLE MASSES=,5D15.7/(42X,5D15.7)
     +)
      END
Commit	Line	Data
fe4da5cc	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(ONETECMTMNTM2FFQ(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)+BETAPCM(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