luupda.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luupda.f

    
C*********************************************************************  
    
      SUBROUTINE luupda(MUPDA,LFN)  
    
C...Purpose: to facilitate the updating of particle and decay data. 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      common/ludat4/chaf(500)   
      CHARACTER chaf*8  
      SAVE /ludat4/ 
      CHARACTER chinl*80,chkc*4,chvar(19)*9,chlin*72,   
     &chblk(20)*72,chold*12,chtmp*12,chnew*12,chcom*12  
      DATA chvar/ 'KCHG(I,1)','KCHG(I,2)','KCHG(I,3)','PMAS(I,1)',  
     &'PMAS(I,2)','PMAS(I,3)','PMAS(I,4)','MDCY(I,1)','MDCY(I,2)',  
     &'MDCY(I,3)','MDME(I,1)','MDME(I,2)','BRAT(I)  ','KFDP(I,1)',  
     &'KFDP(I,2)','KFDP(I,3)','KFDP(I,4)','KFDP(I,5)','CHAF(I)  '/  
    
C...Write information on file for editing.  
      IF(mstu(12).GE.1) CALL lulist(0)  
      IF(mupda.EQ.1) THEN   
        DO 110 kc=1,mstu(6) 
        WRITE(lfn,1000) kc,chaf(kc),(kchg(kc,j1),j1=1,3),   
     &  (pmas(kc,j2),j2=1,4),mdcy(kc,1) 
        DO 100 idc=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1   
  100   WRITE(lfn,1100) mdme(idc,1),mdme(idc,2),brat(idc),  
     &  (kfdp(idc,j),j=1,5) 
  110   CONTINUE    
    
C...Reset variables and read information from edited file.  
      ELSEIF(mupda.EQ.2) THEN   
        DO 120 i=1,mstu(7)  
        mdme(i,1)=1 
        mdme(i,2)=0 
        brat(i)=0.  
        DO 120 j=1,5    
  120   kfdp(i,j)=0 
        kc=0    
        idc=0   
        ndc=0   
  130   READ(lfn,1200,end=140) chinl    
        IF(chinl(2:5).NE.'    ') THEN   
          chkc=chinl(2:5)   
          IF(kc.NE.0) THEN  
            mdcy(kc,2)=0    
            IF(ndc.NE.0) mdcy(kc,2)=idc+1-ndc   
            mdcy(kc,3)=ndc  
          ENDIF 
          READ(chkc,1300) kc    
          IF(kc.LE.0.OR.kc.GT.mstu(6)) CALL luerrm(27,  
     &    '(LUUPDA:) Read KC code illegal, KC ='//chkc) 
          READ(chinl,1000) kcr,chaf(kc),(kchg(kc,j1),j1=1,3),   
     &    (pmas(kc,j2),j2=1,4),mdcy(kc,1)   
          ndc=0 
        ELSE    
          idc=idc+1 
          ndc=ndc+1 
          IF(idc.GE.mstu(7)) CALL luerrm(27,    
     &    '(LUUPDA:) Decay data arrays full by KC ='//chkc) 
          READ(chinl,1100) mdme(idc,1),mdme(idc,2),brat(idc),   
     &    (kfdp(idc,j),j=1,5)   
        ENDIF   
        goto 130    
  140   mdcy(kc,2)=0    
        IF(ndc.NE.0) mdcy(kc,2)=idc+1-ndc   
        mdcy(kc,3)=ndc  
    
C...Perform possible tests that new information is consistent.  
        mstj24=mstj(24) 
        mstj(24)=0  
        DO 170 kc=1,mstu(6) 
        WRITE(chkc,1300) kc 
        IF(min(pmas(kc,1),pmas(kc,2),pmas(kc,3),pmas(kc,1)-pmas(kc,3),  
     &  pmas(kc,4)).LT.0..OR.mdcy(kc,3).LT.0) CALL luerrm(17,   
     &  '(LUUPDA:) Mass/width/life/(# channels) wrong for KC ='//chkc)  
        brsum=0.    
        DO 160 idc=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1   
        IF(mdme(idc,2).GT.80) goto 160  
        kq=kchg(kc,1)   
        pms=pmas(kc,1)-pmas(kc,3)-parj(64)  
        merr=0  
        DO 150 j=1,5    
        kp=kfdp(idc,j)  
        IF(kp.EQ.0.OR.kp.EQ.81.OR.iabs(kp).EQ.82) THEN  
        ELSEIF(lucomp(kp).EQ.0) THEN    
          merr=3    
        ELSE    
          kq=kq-luchge(kp)  
          pms=pms-ulmass(kp)    
        ENDIF   
  150   CONTINUE    
        IF(kq.NE.0) merr=max(2,merr)    
        IF(kfdp(idc,2).NE.0.AND.(kc.LE.20.OR.kc.GT.40).AND. 
     &  (kc.LE.80.OR.kc.GT.100).AND.mdme(idc,2).NE.34.AND.  
     &  mdme(idc,2).NE.61.AND.pms.LT.0.) merr=max(1,merr)   
        IF(merr.EQ.3) CALL luerrm(17,   
     &  '(LUUPDA:) Unknown particle code in decay of KC ='//chkc)   
        IF(merr.EQ.2) CALL luerrm(17,   
     &  '(LUUPDA:) Charge not conserved in decay of KC ='//chkc)    
        IF(merr.EQ.1) CALL luerrm(7,    
     &  '(LUUPDA:) Kinematically unallowed decay of KC ='//chkc)    
        brsum=brsum+brat(idc)   
  160   CONTINUE    
        WRITE(chtmp,1500) brsum 
        IF(abs(brsum).GT.0.0005.AND.abs(brsum-1.).GT.0.0005) CALL   
     &  luerrm(7,'(LUUPDA:) Sum of branching ratios is '//chtmp(5:12)// 
     &  ' for KC ='//chkc)  
  170   CONTINUE    
        mstj(24)=mstj24 
    
C...Initialize writing of DATA statements for inclusion in program. 
      ELSEIF(mupda.EQ.3) THEN   
        DO 240 ivar=1,19    
        ndim=mstu(6)    
        IF(ivar.GE.11.AND.ivar.LE.18) ndim=mstu(7)  
        nlin=1  
        chlin=' '   
        chlin(7:35)='DATA ('//chvar(ivar)//',I=   1,    )/' 
        llin=35 
        chold='START'   
    
C...Loop through variables for conversion to characters.    
        DO 220 idim=1,ndim  
        IF(ivar.EQ.1) WRITE(chtmp,1400) kchg(idim,1)    
        IF(ivar.EQ.2) WRITE(chtmp,1400) kchg(idim,2)    
        IF(ivar.EQ.3) WRITE(chtmp,1400) kchg(idim,3)    
        IF(ivar.EQ.4) WRITE(chtmp,1500) pmas(idim,1)    
        IF(ivar.EQ.5) WRITE(chtmp,1500) pmas(idim,2)    
        IF(ivar.EQ.6) WRITE(chtmp,1500) pmas(idim,3)    
        IF(ivar.EQ.7) WRITE(chtmp,1500) pmas(idim,4)    
        IF(ivar.EQ.8) WRITE(chtmp,1400) mdcy(idim,1)    
        IF(ivar.EQ.9) WRITE(chtmp,1400) mdcy(idim,2)    
        IF(ivar.EQ.10) WRITE(chtmp,1400) mdcy(idim,3)   
        IF(ivar.EQ.11) WRITE(chtmp,1400) mdme(idim,1)   
        IF(ivar.EQ.12) WRITE(chtmp,1400) mdme(idim,2)   
        IF(ivar.EQ.13) WRITE(chtmp,1500) brat(idim) 
        IF(ivar.EQ.14) WRITE(chtmp,1400) kfdp(idim,1)   
        IF(ivar.EQ.15) WRITE(chtmp,1400) kfdp(idim,2)   
        IF(ivar.EQ.16) WRITE(chtmp,1400) kfdp(idim,3)   
        IF(ivar.EQ.17) WRITE(chtmp,1400) kfdp(idim,4)   
        IF(ivar.EQ.18) WRITE(chtmp,1400) kfdp(idim,5)   
        IF(ivar.EQ.19) chtmp=chaf(idim) 
    
C...Length of variable, trailing decimal zeros, quotation marks.    
        llow=1  
        lhig=1  
        DO 180 ll=1,12  
        IF(chtmp(13-ll:13-ll).NE.' ') llow=13-ll    
  180   IF(chtmp(ll:ll).NE.' ') lhig=ll 
        chnew=chtmp(llow:lhig)//' ' 
        lnew=1+lhig-llow    
        IF((ivar.GE.4.AND.ivar.LE.7).OR.ivar.EQ.13) THEN    
          lnew=lnew+1   
  190     lnew=lnew-1   
          IF(chnew(lnew:lnew).EQ.'0') goto 190  
          IF(lnew.EQ.1) chnew(1:2)='0.' 
          IF(lnew.EQ.1) lnew=2  
        ELSEIF(ivar.EQ.19) THEN 
          DO 200 ll=lnew,1,-1   
          IF(chnew(ll:ll).EQ.'''') THEN 
            chtmp=chnew 
            chnew=chtmp(1:ll)//''''//chtmp(ll+1:11) 
            lnew=lnew+1 
          ENDIF 
  200     CONTINUE  
          chtmp=chnew   
          chnew(1:lnew+2)=''''//chtmp(1:lnew)//'''' 
          lnew=lnew+2   
        ENDIF   
    
C...Form composite character string, often including repetition counter.    
        IF(chnew.NE.chold) THEN 
          nrpt=1    
          chold=chnew   
          chcom=chnew   
          lcom=lnew 
        ELSE    
          lrpt=lnew+1   
          IF(nrpt.GE.2) lrpt=lnew+3 
          IF(nrpt.GE.10) lrpt=lnew+4    
          IF(nrpt.GE.100) lrpt=lnew+5   
          IF(nrpt.GE.1000) lrpt=lnew+6  
          llin=llin-lrpt    
          nrpt=nrpt+1   
          WRITE(chtmp,1400) nrpt    
          lrpt=1    
          IF(nrpt.GE.10) lrpt=2 
          IF(nrpt.GE.100) lrpt=3    
          IF(nrpt.GE.1000) lrpt=4   
          chcom(1:lrpt+1+lnew)=chtmp(13-lrpt:12)//'*'//chnew(1:lnew)    
          lcom=lrpt+1+lnew  
        ENDIF   
    
C...Add characters to end of line, to new line (after storing old line),    
C...or to new block of lines (after writing old block). 
        IF(llin+lcom.LE.70) THEN    
          chlin(llin+1:llin+lcom+1)=chcom(1:lcom)//','  
          llin=llin+lcom+1  
        ELSEIF(nlin.LE.19) THEN 
          chlin(llin+1:72)=' '  
          chblk(nlin)=chlin 
          nlin=nlin+1   
          chlin(6:6+lcom+1)='&'//chcom(1:lcom)//',' 
          llin=6+lcom+1 
        ELSE    
          chlin(llin:72)='/'//' '   
          chblk(nlin)=chlin 
          WRITE(chtmp,1400) idim-nrpt   
          chblk(1)(30:33)=chtmp(9:12)   
          DO 210 ilin=1,nlin    
  210     WRITE(lfn,1600) chblk(ilin)   
          nlin=1    
          chlin=' ' 
          chlin(7:35+lcom+1)='DATA ('//chvar(ivar)//',I=    ,    )/'//  
     &    chcom(1:lcom)//','    
          WRITE(chtmp,1400) idim-nrpt+1 
          chlin(25:28)=chtmp(9:12)  
          llin=35+lcom+1    
        ENDIF   
  220   CONTINUE    
    
C...Write final block of lines. 
        chlin(llin:72)='/'//' ' 
        chblk(nlin)=chlin   
        WRITE(chtmp,1400) ndim  
        chblk(1)(30:33)=chtmp(9:12) 
        DO 230 ilin=1,nlin  
  230   WRITE(lfn,1600) chblk(ilin) 
  240   CONTINUE    
      ENDIF 
    
C...Formats for reading and writing particle data.  
 1000 FORMAT(1x,i4,2x,a8,3i3,3f12.5,2x,f12.5,i3)    
 1100 FORMAT(5x,2i5,f12.5,5i8)  
 1200 FORMAT(a80)   
 1300 FORMAT(i4)    
 1400 FORMAT(i12)   
 1500 FORMAT(f12.5) 
 1600 FORMAT(a72)   
    
      RETURN    
      END