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