*
* $Id: genbod.F,v 1.1.1.1 1996/03/22 16:42:49 mclareni Exp $
*
* $Log: genbod.F,v $
* Revision 1.1.1.1  1996/03/22 16:42:49  mclareni
* Phtools
*
*
      SUBROUTINE GENBOD
C   SUBROUTINE TO GENERATE N-BODY EVENT
C   ACCORDING TO FERMI LORENTZ-INVARIANT PHASE SPACE
C   ADAPTED FROM FOWL (CERN W505) SEPT. 1974 BY F. JAMES
C   EVENTS ARE GENERATED IN THEIR OWN CENTER-OF-MASS,
C   BUT MAY BE TRANSFORMED TO ANY FRAME USING LOREN4
C
C   INPUT TO SUBROUTINE IS THRU COMMON BLOCK GENIN
C             NP=NUMBER OF OUTGOING PARTICLES (.LT. 19)
C             TECM=TOTAL ENERGY IN CENTER-OF-MASS
C             AMASS(I)=MASS OF ITH OUTGOING PARTICLE
C             KGENEV=1 FOR CONSTANT CROSS SECTION
C                      2 FOR FERMI ENERGY-DEPENDANCE
C
C   OUTPUT FROM SUBROUTINE IS THRU COMMON BLOCK GENOUT
C             PCM(1,I)=X-MOMENTUM IF ITH PARTICLE
C             PCM(2,I)=Y-MOMENTUM IF ITH PARTICLE
C             PCM(3,I)=Z-MOMENTUM IF ITH PARTICLE
C             PCM(4,I)=ENERGY OF ITH PARTICLE
C             PCM(5,I)=MOMENTUM OF ITH PARTICLE
C             WT=WEIGHT OF EVENT
      COMMON/GENIN / NP, TECM, AMASS(18), KGENEV
      COMMON/GENOUT/ PCM(5,18) , WT
      DIMENSION EMM(18)
      DIMENSION RNO(50)
C--       PCM1 IS LINEAR EQUIV. OF PCM TO AVOID DOUBLE INDICES
      DIMENSION EM(18),PD(18),EMS(18),SM(18),FFQ(18),PCM1(90)
      EQUIVALENCE (NT,NP),(AMASS(1),EM(1)),(PCM1(1),PCM(1,1))
C FFQ(N)=PI * (TWOPI)**(N-2) / (N-2)FACTORIAL
      DATA FFQ/0.,3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
     2                       256.3704, 268.4705, 240.9780, 189.2637,
     3                       132.1308,  83.0202,  47.4210,  24.8295,
     4                        12.0006,   5.3858,   2.2560,   0.8859/
      DATA KNT,TWOPI/0,6.2831853073/
#if defined(CERNLIB_CDC)
      DATA KGENEV/ 2 /
#endif
C        INITIALIZATION
      KNT=KNT + 1
      IF(KNT.GT.1) GOTO 100
      WRITE(6,1160)
      WRITE(6,1200) NP,TECM,(AMASS(JK),JK=1,NP)
  100 CONTINUE
      IF(NT.LT.2) GOTO 1001
      IF(NT.GT.18) GOTO 1002
      NTM1=NT-1
      NTM2=NT-2
      NTP1=NT+1
      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
C        CONSTANTS DEPENDING ON TECM
      TECMTM=TECM-TM
      IF(TECMTM.LE.0.0) GOTO 1000
      EMM(NT)=TECM
      IF(KGENEV.GT.1) GOTO 400
C        CONSTANT CROSS SECTION AS FUNCTION OF TECM
      EMMAX=TECMTM+EM(1)
      EMMIN=0.0
      WTMAX=1.0
      DO 350 I=2,NT
      EMMIN=EMMIN+EM(I-1)
      EMMAX=EMMAX+EM(I)
  350 WTMAX=WTMAX*PDK(EMMAX,EMMIN,EM(I))
      WTMAXQ=1.0/WTMAX
      GOTO 455
C--      FERMI ENERGY DEPENDENCE FOR CROSS SECTION
  400 WTMAXQ=TECMTM**NTM2*FFQ(NT) / TECM
C        CALCULATION OF WT BASED ON EFFECTIVE MASSES EMM
  455 CONTINUE
C--               FILL RNO WITH 3*NT-4 RANDOM NUMBERS,
C--            OF WHICH THE FIRST NT-2 ARE ORDERED.
      DO 457 I= 1, NTNM4
  457 RNO(I)=RNDM(I)
      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 WT=WTMAXQ
      IR=NTM2
      DO 530 I=1,NTM1
      PD(I)=PDK(EMM(I+1),EMM(I),EM(I+1))
  530 WT=WT*PD(I)
C--       COMPLETE SPECIFICATION OF EVENT (RAUBOLD-LYNCH METHOD)
      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(1.0-C*C)
      IF(I.EQ.NT) GOTO 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
      GOTO 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
C          ERROR RETURNS
 1000 WRITE(6,1100)
      GOTO 1050
 1001 WRITE(6,1101)
      GOTO 1050
 1002 WRITE(6,1102)
 1050 WRITE(6,1150) KNT
      WRITE(6,1200) NP,TECM,(AMASS(JK),JK=1,NP)
      STOP
 1100 FORMAT(28H0 AVAILABLE ENERGY NEGATIVE )
 1101 FORMAT(33H0 LESS THAN 2 OUTGOING PARTICLES )
 1102 FORMAT(34H0 MORE THAN 18 OUTGOING PARTICLES )
 1150 FORMAT(47H0 ABOVE ERROR DETECTED IN GENBOD AT CALL NUMBER,I7)
 1160 FORMAT(34H0 FIRST CALL TO SUBROUTINE GENBOD )
 1200 FORMAT(36H  INPUT DATA TO GENBOD.         NP=   ,I6/
     +  ,8H   TECM=,E16.7,18H  PARTICLE MASSES=,5E15.5/(42X,5E15.5)
     +)
      END
