pythiaeRHIC/pythia/pyupda.f

0001
0002 C*********************************************************************
0003
0004 C...PYUPDA
0005 C...Facilitates the updating of particle and decay data
0006 C...by allowing it to be done in an external file.
0007
0008       SUBROUTINE PYUPDA(MUPDA,LFN)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0017       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0018       COMMON/PYDAT4/CHAF(500,2)
0019       CHARACTER CHAF*16
0020       COMMON/PYINT4/MWID(500),WIDS(500,5)
0021       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYDAT4/,/PYINT4/
0022 C...Local arrays, character variables and data.
0023       CHARACTER CHINL*120,CHKF*9,CHVAR(22)*9,CHLIN*72,
0024      &CHBLK(20)*72,CHOLD*16,CHTMP*16,CHNEW*16,CHCOM*24
0025       DATA CHVAR/ 'KCHG(I,1)','KCHG(I,2)','KCHG(I,3)','KCHG(I,4)',
0026      &'PMAS(I,1)','PMAS(I,2)','PMAS(I,3)','PMAS(I,4)','MDCY(I,1)',
0027      &'MDCY(I,2)','MDCY(I,3)','MDME(I,1)','MDME(I,2)','BRAT(I)  ',
0028      &'KFDP(I,1)','KFDP(I,2)','KFDP(I,3)','KFDP(I,4)','KFDP(I,5)',
0029      &'CHAF(I,1)','CHAF(I,2)','MWID(I)  '/
0030
0031 C...Write header if not yet done.
0032       IF(MSTU(12).NE.12345) CALL PYLIST(0)
0033
0034 C...Write information on file for editing.
0035       IF(MUPDA.EQ.1) THEN
0036         DO 110 KC=1,500
0037           WRITE(LFN,5000) KCHG(KC,4),(CHAF(KC,J1),J1=1,2),
0038      &    (KCHG(KC,J2),J2=1,3),(PMAS(KC,J3),J3=1,4),
0039      &    MWID(KC),MDCY(KC,1)
0040           DO 100 IDC=MDCY(KC,2),MDCY(KC,2)+MDCY(KC,3)-1
0041             WRITE(LFN,5100) MDME(IDC,1),MDME(IDC,2),BRAT(IDC),
0042      &      (KFDP(IDC,J),J=1,5)
0043   100     CONTINUE
0044   110   CONTINUE
0045
0046 C...Read complete set of information from edited file or
0047 C...read partial set of new or updated information from edited file.
0048       ELSEIF(MUPDA.EQ.2.OR.MUPDA.EQ.3) THEN
0049
0050 C...Reset counters.
0051         KCC=100
0052         NDC=0
0053         CHKF='         '
0054         IF(MUPDA.EQ.2) THEN
0055           DO 120 I=1,MSTU(6)
0056             KCHG(I,4)=0
0057   120     CONTINUE
0058         ELSE
0059           DO 130 KC=1,MSTU(6)
0060             IF(KC.GT.100.AND.KCHG(KC,4).GT.100) KCC=KC
0061             NDC=MAX(NDC,MDCY(KC,2)+MDCY(KC,3)-1)
0062   130     CONTINUE
0063         ENDIF
0064
0065 C...Begin of loop: read new line; unknown whether particle or
0066 C...decay data.
0067   140   READ(LFN,5200,END=190) CHINL
0068
0069 C...Identify particle code and whether already defined  (for MUPDA=3).
0070         IF(CHINL(2:10).NE.'         ') THEN
0071           CHKF=CHINL(2:10)
0072           READ(CHKF,5300) KF
0073           IF(MUPDA.EQ.2) THEN
0074             IF(KF.LE.100) THEN
0075               KC=KF
0076             ELSE
0077               KCC=KCC+1
0078               KC=KCC
0079             ENDIF
0080           ELSE
0081             KCREP=0
0082             IF(KF.LE.100) THEN
0083               KCREP=KF
0084             ELSE
0085               DO 150 KCR=101,KCC
0086                 IF(KCHG(KCR,4).EQ.KF) KCREP=KCR
0087   150         CONTINUE
0088             ENDIF
0089 C...Remove duplicate old decay data.
0090             IF(KCREP.NE.0.AND.MDCY(KCREP,3).GT.0) THEN
0091               IDCREP=MDCY(KCREP,2)
0092               NDCREP=MDCY(KCREP,3)
0093               DO 160 I=1,KCC
0094                 IF(MDCY(I,2).GT.IDCREP) MDCY(I,2)=MDCY(I,2)-NDCREP
0095   160         CONTINUE
0096               DO 180 I=IDCREP,NDC-NDCREP
0097                 MDME(I,1)=MDME(I+NDCREP,1)
0098                 MDME(I,2)=MDME(I+NDCREP,2)
0099                 BRAT(I)=BRAT(I+NDCREP)
0100                 DO 170 J=1,5
0101                   KFDP(I,J)=KFDP(I+NDCREP,J)
0102   170           CONTINUE
0103   180         CONTINUE
0104               NDC=NDC-NDCREP
0105               KC=KCREP
0106             ELSEIF(KCREP.NE.0) THEN
0107               KC=KCREP
0108             ELSE
0109               KCC=KCC+1
0110               KC=KCC
0111             ENDIF
0112           ENDIF
0113
0114 C...Study line with particle data.
0115           IF(KC.GT.MSTU(6)) CALL PYERRM(27,
0116      &    '(PYUPDA:) Particle arrays full by KF ='//CHKF)
0117           READ(CHINL,5000) KCHG(KC,4),(CHAF(KC,J1),J1=1,2),
0118      &    (KCHG(KC,J2),J2=1,3),(PMAS(KC,J3),J3=1,4),
0119      &    MWID(KC),MDCY(KC,1)
0120           MDCY(KC,2)=0
0121           MDCY(KC,3)=0
0122
0123 C...Study line with decay data.
0124         ELSE
0125           NDC=NDC+1
0126           IF(NDC.GT.MSTU(7)) CALL PYERRM(27,
0127      &    '(PYUPDA:) Decay data arrays full by KF ='//CHKF)
0128           IF(MDCY(KC,2).EQ.0) MDCY(KC,2)=NDC
0129           MDCY(KC,3)=MDCY(KC,3)+1
0130           READ(CHINL,5100) MDME(NDC,1),MDME(NDC,2),BRAT(NDC),
0131      &    (KFDP(NDC,J),J=1,5)
0132         ENDIF
0133
0134 C...End of loop; ensure that PYCOMP tables are updated.
0135         GOTO 140
0136   190   CONTINUE
0137         MSTU(20)=0
0138
0139 C...Perform possible tests that new information is consistent.
0140         DO 220 KC=1,MSTU(6)
0141           KF=KCHG(KC,4)
0142           IF(KF.EQ.0) GOTO 220
0143           WRITE(CHKF,5300) KF
0144           IF(MIN(PMAS(KC,1),PMAS(KC,2),PMAS(KC,3),PMAS(KC,1)-PMAS(KC,3),
0145      &    PMAS(KC,4)).LT.0D0.OR.MDCY(KC,3).LT.0) CALL PYERRM(17,
0146      &    '(PYUPDA:) Mass/width/life/(# channels) wrong for KF ='//CHKF)
0147           BRSUM=0D0
0148           DO 210 IDC=MDCY(KC,2),MDCY(KC,2)+MDCY(KC,3)-1
0149             IF(MDME(IDC,2).GT.80) GOTO 210
0150             KQ=KCHG(KC,1)
0151             PMS=PMAS(KC,1)-PMAS(KC,3)-PARJ(64)
0152             MERR=0
0153             DO 200 J=1,5
0154               KP=KFDP(IDC,J)
0155               IF(KP.EQ.0.OR.KP.EQ.81.OR.IABS(KP).EQ.82) THEN
0156                 IF(KP.EQ.81) KQ=0
0157               ELSEIF(PYCOMP(KP).EQ.0) THEN
0158                 MERR=3
0159               ELSE
0160                 KQ=KQ-PYCHGE(KP)
0161                 KPC=PYCOMP(KP)
0162                 PMS=PMS-PMAS(KPC,1)
0163                 IF(MSTJ(24).GT.0) PMS=PMS+0.5D0*MIN(PMAS(KPC,2),
0164      &          PMAS(KPC,3))
0165               ENDIF
0166   200       CONTINUE
0167             IF(KQ.NE.0) MERR=MAX(2,MERR)
0168             IF(MWID(KC).EQ.0.AND.KF.NE.311.AND.PMS.LT.0D0)
0169      &      MERR=MAX(1,MERR)
0170             IF(MERR.EQ.3) CALL PYERRM(17,
0171      &      '(PYUPDA:) Unknown particle code in decay of KF ='//CHKF)
0172             IF(MERR.EQ.2) CALL PYERRM(17,
0173      &      '(PYUPDA:) Charge not conserved in decay of KF ='//CHKF)
0174             IF(MERR.EQ.1) CALL PYERRM(7,
0175      &      '(PYUPDA:) Kinematically unallowed decay of KF ='//CHKF)
0176             BRSUM=BRSUM+BRAT(IDC)
0177   210     CONTINUE
0178           WRITE(CHTMP,5500) BRSUM
0179           IF(ABS(BRSUM).GT.0.0005D0.AND.ABS(BRSUM-1D0).GT.0.0005D0)
0180      &    CALL PYERRM(7,'(PYUPDA:) Sum of branching ratios is '//
0181      &    CHTMP(9:16)//' for KF ='//CHKF)
0182   220   CONTINUE
0183
0184 C...Write DATA statements for inclusion in program.
0185       ELSEIF(MUPDA.EQ.4) THEN
0186
0187 C...Find out how many codes and decay channels are actually used.
0188         KCC=0
0189         NDC=0
0190         DO 230 I=1,MSTU(6)
0191           IF(KCHG(I,4).NE.0) THEN
0192             KCC=I
0193             NDC=MAX(NDC,MDCY(I,2)+MDCY(I,3)-1)
0194           ENDIF
0195   230   CONTINUE
0196
0197 C...Initialize writing of DATA statements for inclusion in program.
0198         DO 300 IVAR=1,22
0199           NDIM=MSTU(6)
0200           IF(IVAR.GE.12.AND.IVAR.LE.19) NDIM=MSTU(7)
0201           NLIN=1
0202           CHLIN=' '
0203           CHLIN(7:35)='DATA ('//CHVAR(IVAR)//',I=   1,    )/'
0204           LLIN=35
0205           CHOLD='START'
0206
0207 C...Loop through variables for conversion to characters.
0208           DO 280 IDIM=1,NDIM
0209             IF(IVAR.EQ.1) WRITE(CHTMP,5400) KCHG(IDIM,1)
0210             IF(IVAR.EQ.2) WRITE(CHTMP,5400) KCHG(IDIM,2)
0211             IF(IVAR.EQ.3) WRITE(CHTMP,5400) KCHG(IDIM,3)
0212             IF(IVAR.EQ.4) WRITE(CHTMP,5400) KCHG(IDIM,4)
0213             IF(IVAR.EQ.5) WRITE(CHTMP,5500) PMAS(IDIM,1)
0214             IF(IVAR.EQ.6) WRITE(CHTMP,5500) PMAS(IDIM,2)
0215             IF(IVAR.EQ.7) WRITE(CHTMP,5500) PMAS(IDIM,3)
0216             IF(IVAR.EQ.8) WRITE(CHTMP,5500) PMAS(IDIM,4)
0217             IF(IVAR.EQ.9) WRITE(CHTMP,5400) MDCY(IDIM,1)
0218             IF(IVAR.EQ.10) WRITE(CHTMP,5400) MDCY(IDIM,2)
0219             IF(IVAR.EQ.11) WRITE(CHTMP,5400) MDCY(IDIM,3)
0220             IF(IVAR.EQ.12) WRITE(CHTMP,5400) MDME(IDIM,1)
0221             IF(IVAR.EQ.13) WRITE(CHTMP,5400) MDME(IDIM,2)
0222             IF(IVAR.EQ.14) WRITE(CHTMP,5600) BRAT(IDIM)
0223             IF(IVAR.EQ.15) WRITE(CHTMP,5400) KFDP(IDIM,1)
0224             IF(IVAR.EQ.16) WRITE(CHTMP,5400) KFDP(IDIM,2)
0225             IF(IVAR.EQ.17) WRITE(CHTMP,5400) KFDP(IDIM,3)
0226             IF(IVAR.EQ.18) WRITE(CHTMP,5400) KFDP(IDIM,4)
0227             IF(IVAR.EQ.19) WRITE(CHTMP,5400) KFDP(IDIM,5)
0228             IF(IVAR.EQ.20) CHTMP=CHAF(IDIM,1)
0229             IF(IVAR.EQ.21) CHTMP=CHAF(IDIM,2)
0230             IF(IVAR.EQ.22) WRITE(CHTMP,5400) MWID(IDIM)
0231
0232 C...Replace variables beyond what is properly defined.
0233             IF(IVAR.LE.4) THEN
0234               IF(IDIM.GT.KCC) CHTMP='               0'
0235             ELSEIF(IVAR.LE.8) THEN
0236               IF(IDIM.GT.KCC) CHTMP='             0.0'
0237             ELSEIF(IVAR.LE.11) THEN
0238               IF(IDIM.GT.KCC) CHTMP='               0'
0239             ELSEIF(IVAR.LE.13) THEN
0240               IF(IDIM.GT.NDC) CHTMP='               0'
0241             ELSEIF(IVAR.LE.14) THEN
0242               IF(IDIM.GT.NDC) CHTMP='             0.0'
0243             ELSEIF(IVAR.LE.19) THEN
0244               IF(IDIM.GT.NDC) CHTMP='               0'
0245             ELSEIF(IVAR.LE.21) THEN
0246               IF(IDIM.GT.KCC) CHTMP='                '
0247             ELSE
0248               IF(IDIM.GT.KCC) CHTMP='               0'
0249             ENDIF
0250
0251 C...Length of variable, trailing decimal zeros, quotation marks.
0252             LLOW=1
0253             LHIG=1
0254             DO 240 LL=1,16
0255               IF(CHTMP(17-LL:17-LL).NE.' ') LLOW=17-LL
0256               IF(CHTMP(LL:LL).NE.' ') LHIG=LL
0257   240       CONTINUE
0258             CHNEW=CHTMP(LLOW:LHIG)//' '
0259             LNEW=1+LHIG-LLOW
0260             IF((IVAR.GE.5.AND.IVAR.LE.8).OR.IVAR.EQ.14) THEN
0261               LNEW=LNEW+1
0262   250         LNEW=LNEW-1
0263               IF(LNEW.GE.2.AND.CHNEW(LNEW:LNEW).EQ.'0') GOTO 250
0264               IF(CHNEW(LNEW:LNEW).EQ.'.') LNEW=LNEW-1
0265               IF(LNEW.EQ.0) THEN
0266                 CHNEW(1:3)='0D0'
0267                 LNEW=3
0268               ELSE
0269                 CHNEW(LNEW+1:LNEW+2)='D0'
0270                 LNEW=LNEW+2
0271               ENDIF
0272             ELSEIF(IVAR.EQ.20.OR.IVAR.EQ.21) THEN
0273               DO 260 LL=LNEW,1,-1
0274                 IF(CHNEW(LL:LL).EQ.'''') THEN
0275                   CHTMP=CHNEW
0276                   CHNEW=CHTMP(1:LL)//''''//CHTMP(LL+1:11)
0277                   LNEW=LNEW+1
0278                 ENDIF
0279   260         CONTINUE
0280               LNEW=MIN(14,LNEW)
0281               CHTMP=CHNEW
0282               CHNEW(1:LNEW+2)=''''//CHTMP(1:LNEW)//''''
0283               LNEW=LNEW+2
0284             ENDIF
0285
0286 C...Form composite character string, often including repetition counter.
0287             IF(CHNEW.NE.CHOLD) THEN
0288               NRPT=1
0289               CHOLD=CHNEW
0290               CHCOM=CHNEW
0291               LCOM=LNEW
0292             ELSE
0293               LRPT=LNEW+1
0294               IF(NRPT.GE.2) LRPT=LNEW+3
0295               IF(NRPT.GE.10) LRPT=LNEW+4
0296               IF(NRPT.GE.100) LRPT=LNEW+5
0297               IF(NRPT.GE.1000) LRPT=LNEW+6
0298               LLIN=LLIN-LRPT
0299               NRPT=NRPT+1
0300               WRITE(CHTMP,5400) NRPT
0301               LRPT=1
0302               IF(NRPT.GE.10) LRPT=2
0303               IF(NRPT.GE.100) LRPT=3
0304               IF(NRPT.GE.1000) LRPT=4
0305               CHCOM(1:LRPT+1+LNEW)=CHTMP(17-LRPT:16)//'*'//CHNEW(1:LNEW)
0306               LCOM=LRPT+1+LNEW
0307             ENDIF
0308
0309 C...Add characters to end of line, to new line (after storing old line),
0310 C...or to new block of lines (after writing old block).
0311             IF(LLIN+LCOM.LE.70) THEN
0312               CHLIN(LLIN+1:LLIN+LCOM+1)=CHCOM(1:LCOM)//','
0313               LLIN=LLIN+LCOM+1
0314             ELSEIF(NLIN.LE.19) THEN
0315               CHLIN(LLIN+1:72)=' '
0316               CHBLK(NLIN)=CHLIN
0317               NLIN=NLIN+1
0318               CHLIN(6:6+LCOM+1)='&'//CHCOM(1:LCOM)//','
0319               LLIN=6+LCOM+1
0320             ELSE
0321               CHLIN(LLIN:72)='/'//' '
0322               CHBLK(NLIN)=CHLIN
0323               WRITE(CHTMP,5400) IDIM-NRPT
0324               CHBLK(1)(30:33)=CHTMP(13:16)
0325               DO 270 ILIN=1,NLIN
0326                 WRITE(LFN,5700) CHBLK(ILIN)
0327   270         CONTINUE
0328               NLIN=1
0329               CHLIN=' '
0330               CHLIN(7:35+LCOM+1)='DATA ('//CHVAR(IVAR)//
0331      &        ',I=    ,    )/'//CHCOM(1:LCOM)//','
0332               WRITE(CHTMP,5400) IDIM-NRPT+1
0333               CHLIN(25:28)=CHTMP(13:16)
0334               LLIN=35+LCOM+1
0335             ENDIF
0336   280     CONTINUE
0337
0338 C...Write final block of lines.
0339           CHLIN(LLIN:72)='/'//' '
0340           CHBLK(NLIN)=CHLIN
0341           WRITE(CHTMP,5400) NDIM
0342           CHBLK(1)(30:33)=CHTMP(13:16)
0343           DO 290 ILIN=1,NLIN
0344             WRITE(LFN,5700) CHBLK(ILIN)
0345   290     CONTINUE
0346   300   CONTINUE
0347       ENDIF
0348
0349 C...Formats for reading and writing particle data.
0350  5000 FORMAT(1X,I9,2X,A16,2X,A16,3I3,3F12.5,1P,E13.5,2I3)
0351  5100 FORMAT(10X,2I5,F12.6,5I10)
0352  5200 FORMAT(A120)
0353  5300 FORMAT(I9)
0354  5400 FORMAT(I16)
0355  5500 FORMAT(F16.5)
0356  5600 FORMAT(F16.6)
0357  5700 FORMAT(A72)
0358
0359       RETURN
0360       END