GEANT321/gheisha/phpnuc.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:06  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.40  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE PHPNUC
  13 C
  14 C *** DOUBLE PRECISION VERSION OF THE PHASE SPACE ROUTINE "PHASP"
  15 C *** THIS ROUTINE MUST BE CALLED BY THE NUCLEAR INTERACTION ROUTINE
  16 C *** "NUCREC" (SEE ALSO COMMENTS THEREIN). THE REASON IS SIMPLY THAT
  17 C *** ENERGY-MOMENTUM CALCULATIONS ARE NOT POSSIBLE WITHIN ONLY
  18 C *** 6 DIGITS OF ACCURACY FOR TOTAL ENERGIES
  19 C *** IN THE ORDER OF HUNDREDS OF GEV (URANIUM NUCLEUS), COMPARED WITH
  20 C *** KINETIC ENERGIES IN THE ORDER OF MEV (NEUTRONS, PROTONS AND
  21 C *** PHOTONS IN THE REACTIONS A(X,Y(GAMMA,GAMMA))A'). IN THE ORIGINAL
  22 C *** GHEISHA8 CODE ALL THESE CALCULATIONS ARE DONE IN DOUBLE PRECISION
  23 C *** HMF 29-AUG-1989 RWTH AACHEN
  24 C
  25 C CALLED BY : NUCREC
  26 C ORIGIN    : H.FESEFELDT
  27 C
  28 #if !defined(CERNLIB_SINGLE)
  29       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  30 #endif
  31       REAL RNDM(1)
  32 C
  33 #include "geant321/s_prntfl.inc"
  34 #include "geant321/s_nucio.inc"
  35 C
  36       SAVE  KNT, TWOPI, FFQ
  37       DIMENSION EMM(18)
  38       DIMENSION RNO(50)
  39       DIMENSION EM(18),PD(18),EMS(18),SM(18),FFQ(18),PCM1(90)
  40       EQUIVALENCE (NT,NPG),(AMASS(1),EM(1)),(PCM1(1),PCM(1,1))
  41       DATA  FFQ/0.,3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
  42      2                       256.3704, 268.4705, 240.9780, 189.2637,
  43      3                       132.1308,  83.0202,  47.4210,  24.8295,
  44      4                        12.0006,   5.3858,   2.2560,   0.8859/
  45       DATA  KNT , TWOPI /  1 , 6.2831853073 /
  46 C
  47 C --- Initialise local arrays and the result array PCM ---
  48       CALL VZERO(PCM,90)
  49       DO 80 JZERO=1,18
  50          EMM(JZERO)=0.
  51          PD(JZERO) =0.
  52          EMS(JZERO)=0.
  53          SM(JZERO) =0.
  54   80  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       NTP1=NT+1
  65       NTNM4 = 3*NT - 4
  66       EMM(1)=EM(1)
  67       TM=0.0
  68       DO 200 I=1,NT
  69       EMS(I)=EM(I)**2
  70       TM=TM+EM(I)
  71  200  SM(I)=TM
  72       WGT=1.
  73  210  TECMTM=TECM-TM
  74       IF (TECMTM .LE. 0.0)  GO TO 1000
  75       EMM(NT)=TECM
  76       IF (KGENEV.GT.1) GO TO 400
  77       EMMAX=TECMTM+EM(1)
  78       EMMIN=0.0
  79       WTMAX=1.0
  80       DO 350 I=2,NT
  81       EMMIN=EMMIN+EM(I-1)
  82       EMMAX=EMMAX+EM(I)
  83   350 WTMAX=WTMAX*DPDNUC(EMMAX,EMMIN,EM(I))
  84       WTMAXQ=1.0/WTMAX
  85       GO TO 455
  86   400 WTMAXQ=TECMTM**NTM2*FFQ(NT) / TECM
  87   455 CONTINUE
  88       DO 457 I= 1, NTNM4
  89       CALL GRNDM(RNDM,1)
  90 #if defined(CERNLIB_SINGLE)
  91   457 RNO(I) = RNDM(1)
  92 #endif
  93 #if !defined(CERNLIB_SINGLE)
  94   457 RNO(I) = DBLE(RNDM(1))
  95 #endif
  96       IF(NTM2) 900,509,460
  97   460 CONTINUE
  98       CALL DLPNUC(RNO,NTM2)
  99       DO 508 J=2,NTM1
 100   508 EMM(J)=RNO(J-1)*(TECMTM)+SM(J)
 101   509 WGT=WTMAXQ
 102       IR=NTM2
 103       DO 530 I=1,NTM1
 104       PD(I)=DPDNUC(EMM(I+1),EMM(I),EM(I+1))
 105   530 WGT=WGT*PD(I)
 106       PCM(1,1)=0.0
 107       PCM(2,1)=PD(1)
 108       PCM(3,1)=0.0
 109       DO 570 I=2,NT
 110       PCM(1,I)=0.0
 111       PCM(2,I) = -PD(I-1)
 112       PCM(3,I)=0.0
 113       IR=IR+1
 114       BANG=TWOPI*RNO(IR)
 115       CB=COS(BANG)
 116       SB=SIN(BANG)
 117       IR=IR+1
 118       C=2.0*RNO(IR)-1.0
 119       S=SQRT(1.0-C*C)
 120       IF(I.EQ.NT) GO TO 1567
 121       ESYS=SQRT(PD(I)**2+EMM(I)**2)
 122       BETA=PD(I)/ESYS
 123       GAMA=ESYS/EMM(I)
 124       DO 568 J=1,I
 125       NDX = 5*J - 5
 126       AA= PCM1(NDX+1)**2 + PCM1(NDX+2)**2 + PCM1(NDX+3)**2
 127       PCM1(NDX+5) = SQRT(AA)
 128       PCM1(NDX+4) = SQRT(AA+EMS(J))
 129       CALL DOTNUC(C,S,CB,SB,PCM,J)
 130       PSAVE = GAMA*(PCM(2,J)+BETA*PCM(4,J))
 131   568 PCM(2,J)=PSAVE
 132       GO TO 570
 133  1567 DO 1568 J=1,I
 134       AA=PCM(1,J)**2 + PCM(2,J)**2 + PCM(3,J)**2
 135       PCM(5,J)=SQRT(AA)
 136       PCM(4,J)=SQRT(AA+EMS(J))
 137       CALL DOTNUC(C,S,CB,SB,PCM,J)
 138  1568 CONTINUE
 139   570 CONTINUE
 140   900 CONTINUE
 141       RETURN
 142  1000 DO 212 I=1,NPG
 143       PCM(1,I)=0.
 144       PCM(2,I)=0.
 145       PCM(3,I)=0.
 146       PCM(4,I)=AMASS(I)
 147   212 PCM(5,I)=AMASS(I)
 148       WGT=0.
 149       RETURN
 150  1001 IF(NPRT(3).OR.NPRT(4)) WRITE(NEWBCD,1101)
 151       GO TO 1050
 152  1002 WRITE(NEWBCD,1102)
 153  1050 WRITE(NEWBCD,1150) KNT
 154       WRITE(NEWBCD,1200) NPG,TECM,(AMASS(JK),JK=1,NPG)
 155       STOP
 156  1100 FORMAT ('0*PHPNUC* AVAILABLE ENERGY NEGATIVE')
 157  1101 FORMAT ('0*PHPNUC* LESS THAN 2 OUTGOING PARTICLES')
 158  1102 FORMAT ('0*PHPNUC* MORE THAN 18 OUTGOING PARTICLES')
 159  1150 FORMAT ('0*PHPNUC* ABOVE ERROR DETECTED IN PHASP',
 160      $ ' AT CALL NUMBER ',I7)
 161  1200 FORMAT ('0*PHPNUC* INPUT DATA TO PHPNUC. NPG = ',I6/
 162      $ ' TECM = ',E15.7,' PARTICLE MASSES = ',5E15.7/(42X,5E15.7))
 163       END