More volume overlaps corrected
[u/mrichter/AliRoot.git] / ISAJET / code / decps1.F
1 #include "isajet/pilot.h"
2       SUBROUTINE DECPS1(IP,NADD,PGEN)
3 C
4 C          Generate masses for uniform NADD-body phase space in DECPS2.
5 C          Auxiliary routine for DECAY.
6 C
7 #if defined(CERNLIB_IMPNONE)
8       IMPLICIT NONE
9 #endif
10 C
11 #include "isajet/itapes.inc"
12 #include "isajet/partcl.inc"
13 C
14       INTEGER IP,NADD
15       REAL PGEN(5,5)
16       REAL REDUCE(5),RND(5)
17       REAL RANF,PCM,DBLPCM
18       REAL WTMAX,SUM1,SUM2,SUM,RNEW,WT,A,B,C
19       INTEGER I,NADD1,J,I1,JJ1,JSAVE
20 C
21 C          Function definitions.
22 C
23 #if defined(CERNLIB_SINGLE)
24       PCM(A,B,C)=SQRT((A-B-C)*(A+B+C)*(A-B+C)*(A+B-C))/(2.*A)
25 #endif
26 #if defined(CERNLIB_DOUBLE)
27       PCM(A,B,C)=DBLPCM(A,B,C)
28 #endif
29 C
30       DATA REDUCE/1.,1.,2.,5.,15./
31 C
32 C          Calculate maximum phase-space weight.
33 C
34       IF(NADD.LE.2) RETURN
35       NADD1=NADD-1
36       WTMAX=1./REDUCE(NADD)
37       SUM=0
38       DO 100 I=1,NADD
39         SUM=SUM+PPTCL(5,NPTCL+I)
40 100   CONTINUE
41       SUM1=PGEN(5,1)
42       SUM2=SUM-PPTCL(5,NPTCL+1)
43       DO 110 I=1,NADD1
44         WTMAX=WTMAX*PCM(SUM1,SUM2,PPTCL(5,NPTCL+I))
45         SUM1=SUM1-PPTCL(5,NPTCL+I)
46         SUM2=SUM2-PPTCL(5,NPTCL+I+1)
47 110   CONTINUE
48 C
49 C          Generate masses for uniform NADD-body phase space.
50 C
51 200   CONTINUE
52       RND(1)=1.
53       DO 210 I=2,NADD1
54         RNEW=RANF()
55         I1=I-1
56         DO 220 JJ1=1,I1
57           J=I-JJ1
58           JSAVE=J+1
59           IF(RNEW.LE.RND(J)) GO TO 210
60           RND(JSAVE)=RND(J)
61 220     CONTINUE
62 210   RND(JSAVE)=RNEW
63       RND(NADD)=0.
64       WT=1.
65       SUM1=SUM
66       DO 230 I=2,NADD
67         SUM1=SUM1-PPTCL(5,NPTCL+I-1)
68         PGEN(5,I)=SUM1+RND(I)*(PGEN(5,1)-SUM)
69         IF(PGEN(5,I-1).LE.PGEN(5,I)+PPTCL(5,NPTCL+I-1)) GO TO 200
70         WT=WT*PCM(PGEN(5,I-1),PGEN(5,I),PPTCL(5,NPTCL+I-1))
71 230   CONTINUE
72       IF(WT.LT.RANF()*WTMAX) GO TO 200
73 C
74       RETURN
75       END