pyhiinki.f

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

    
C*********************************************************************  
    
      SUBROUTINE pyhiinki(CHFRAM,CHBEAM,CHTARG,WIN)   
    
C...Identifies the two incoming particles and sets up kinematics,   
C...including rotations and boosts to/from CM frame.    
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200) 
      SAVE /pyhisubs/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
      CHARACTER chfram*8,chbeam*8,chtarg*8,chcom(3)*8,chalp(2)*26,  
     &chidnt(3)*8,chtemp*8,chcde(18)*8,chinit*76    
      dimension len(3),kcde(18) 
      DATA chalp/'abcdefghijklmnopqrstuvwxyz',  
     &'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ 
      DATA chcde/'e-      ','e+      ','nue     ','nue~    ',   
     &'mu-     ','mu+     ','numu    ','numu~   ','tau-    ',   
     &'tau+    ','nutau   ','nutau~  ','pi+     ','pi-     ',   
     &'n       ','n~      ','p       ','p~      '/  
      DATA kcde/11,-11,12,-12,13,-13,14,-14,15,-15,16,-16,  
     &211,-211,2112,-2112,2212,-2212/   
    
C...Convert character variables to lowercase and find their length. 
      chcom(1)=chfram   
      chcom(2)=chbeam   
      chcom(3)=chtarg   
      DO 120 i=1,3  
      len(i)=8  
      DO 100 ll=8,1,-1  
      IF(len(i).EQ.ll.AND.chcom(i)(ll:ll).EQ.' ') len(i)=ll-1   
      DO 100 la=1,26    
  100 IF(chcom(i)(ll:ll).EQ.chalp(2)(la:la)) chcom(i)(ll:ll)=   
     &chalp(1)(la:la)   
      chidnt(i)=chcom(i)    
      DO 110 ll=1,6 
      IF(chidnt(i)(ll:ll+2).EQ.'bar') THEN  
        chtemp=chidnt(i)    
        chidnt(i)=chtemp(1:ll-1)//'~'//chtemp(ll+3:8)//'  ' 
      ENDIF 
  110 CONTINUE  
      DO 120 ll=1,8 
      IF(chidnt(i)(ll:ll).EQ.'_') THEN  
        chtemp=chidnt(i)    
        chidnt(i)=chtemp(1:ll-1)//chtemp(ll+1:8)//' '   
      ENDIF 
  120 CONTINUE  
    
C...Set initial state. Error for unknown codes. Reset variables.    
      n=2   
      DO 140 i=1,2  
      k(i,2)=0  
      DO 130 j=1,18 
  130 IF(chidnt(i+1).EQ.chcde(j)) k(i,2)=kcde(j)    
      p(i,5)=ulmass(k(i,2)) 
      mint(40+i)=1  
      IF(iabs(k(i,2)).GT.100) mint(40+i)=2  
      DO 140 j=1,5  
  140 v(i,j)=0. 
      IF(k(1,2).EQ.0) WRITE(mstu(11),1000) chbeam(1:len(2)) 
      IF(k(2,2).EQ.0) WRITE(mstu(11),1100) chtarg(1:len(3)) 
      IF(k(1,2).EQ.0.OR.k(2,2).EQ.0) stop   
      DO 150 j=6,10 
  150 vint(j)=0.    
      chinit=' '    
    
C...Set up kinematics for events defined in CM frame.   
      IF(chcom(1)(1:2).EQ.'cm') THEN    
        IF(chcom(2)(1:1).NE.'e') THEN   
          loffs=(34-(len(2)+len(3)))/2  
          chinit(loffs+1:76)='PYTHIA will be initialized for a '//  
     &    chcom(2)(1:len(2))//'-'//chcom(3)(1:len(3))//' collider'//' ' 
        ELSE    
          loffs=(33-(len(2)+len(3)))/2  
          chinit(loffs+1:76)='PYTHIA will be initialized for an '// 
     &    chcom(2)(1:len(2))//'-'//chcom(3)(1:len(3))//' collider'//' ' 
        ENDIF   
C        WRITE(MSTU(11),1200) CHINIT 
C        WRITE(MSTU(11),1300) WIN    
        s=win**2    
        p(1,1)=0.   
        p(1,2)=0.   
        p(2,1)=0.   
        p(2,2)=0.   
        p(1,3)=sqrt(((s-p(1,5)**2-p(2,5)**2)**2-(2.*p(1,5)*p(2,5))**2)/ 
     &  (4.*s)) 
        p(2,3)=-p(1,3)  
        p(1,4)=sqrt(p(1,3)**2+p(1,5)**2)    
        p(2,4)=sqrt(p(2,3)**2+p(2,5)**2)    
    
C...Set up kinematics for fixed target events.  
      ELSEIF(chcom(1)(1:3).EQ.'fix') THEN   
        loffs=(29-(len(2)+len(3)))/2    
        chinit(loffs+1:76)='PYTHIA will be initialized for '//  
     &  chcom(2)(1:len(2))//' on '//chcom(3)(1:len(3))//    
     &  ' fixed target'//' '    
C        WRITE(MSTU(11),1200) CHINIT 
C        WRITE(MSTU(11),1400) WIN    
        p(1,1)=0.   
        p(1,2)=0.   
        p(2,1)=0.   
        p(2,2)=0.   
        p(1,3)=win  
        p(1,4)=sqrt(p(1,3)**2+p(1,5)**2)    
        p(2,3)=0.   
        p(2,4)=p(2,5)   
        s=p(1,5)**2+p(2,5)**2+2.*p(2,4)*p(1,4)  
        vint(10)=p(1,3)/(p(1,4)+p(2,4)) 
        CALL lurobo(0.,0.,0.,0.,-vint(10))  
C        WRITE(MSTU(11),1500) SQRT(S)    
    
C...Set up kinematics for events in user-defined frame. 
      ELSEIF(chcom(1)(1:3).EQ.'use') THEN   
        loffs=(13-(len(1)+len(2)))/2    
        chinit(loffs+1:76)='PYTHIA will be initialized for '//  
     &  chcom(2)(1:len(2))//' on '//chcom(3)(1:len(3))//    
     &  'user-specified configuration'//' ' 
C        WRITE(MSTU(11),1200) CHINIT 
C        WRITE(MSTU(11),1600)    
C        WRITE(MSTU(11),1700) CHCOM(2),P(1,1),P(1,2),P(1,3)  
C        WRITE(MSTU(11),1700) CHCOM(3),P(2,1),P(2,2),P(2,3)  
        p(1,4)=sqrt(p(1,1)**2+p(1,2)**2+p(1,3)**2+p(1,5)**2)    
        p(2,4)=sqrt(p(2,1)**2+p(2,2)**2+p(2,3)**2+p(2,5)**2)    
        DO 160 j=1,3    
  160   vint(7+j)=(dble(p(1,j))+dble(p(2,j)))/dble(p(1,4)+p(2,4))   
        CALL lurobo(0.,0.,-vint(8),-vint(9),-vint(10))  
        vint(7)=ulangl(p(1,1),p(1,2))   
        CALL lurobo(0.,-vint(7),0.,0.,0.)   
        vint(6)=ulangl(p(1,3),p(1,1))   
        CALL lurobo(-vint(6),0.,0.,0.,0.)   
        s=p(1,5)**2+p(2,5)**2+2.*(p(1,4)*p(2,4)-p(1,3)*p(2,3))  
C        WRITE(MSTU(11),1500) SQRT(S)    
    
C...Unknown frame. Error for too low CM energy. 
      ELSE  
        WRITE(mstu(11),1800) chfram(1:len(1))   
        stop    
      ENDIF 
      IF(s.LT.parp(2)**2) THEN  
        WRITE(mstu(11),1900) sqrt(s)    
        stop    
      ENDIF 
    
C...Save information on incoming particles. 
      mint(11)=k(1,2)   
      mint(12)=k(2,2)   
      mint(43)=2*mint(41)+mint(42)-2    
      vint(1)=sqrt(s)   
      vint(2)=s 
      vint(3)=p(1,5)    
      vint(4)=p(2,5)    
      vint(5)=p(1,3)    
    
C...Store constants to be used in generation.   
      IF(mstp(82).LE.1) vint(149)=4.*parp(81)**2/s  
      IF(mstp(82).GE.2) vint(149)=4.*parp(82)**2/s  
    
C...Formats for initialization and error information.   
 1000 FORMAT(1x,'Error: unrecognized beam particle ''',a,'''.'/ 
     &1x,'Execution stopped!')  
 1100 FORMAT(1x,'Error: unrecognized target particle ''',a,'''.'/   
     &1x,'Execution stopped!')  
 1200 FORMAT(/1x,78('=')/1x,'I',76x,'I'/1x,'I',a76,'I') 
 1300 FORMAT(1x,'I',18x,'at',1x,f10.3,1x,'GeV center-of-mass energy',   
     &19x,'I'/1x,'I',76x,'I'/1x,78('='))    
 1400 FORMAT(1x,'I',22x,'at',1x,f10.3,1x,'GeV/c lab-momentum',22x,'I')  
 1500 FORMAT(1x,'I',76x,'I'/1x,'I',11x,'corresponding to',1x,f10.3,1x,  
     &'GeV center-of-mass energy',12x,'I'/1x,'I',76x,'I'/1x,78('='))    
 1600 FORMAT(1x,'I',76x,'I'/1x,'I',24x,'px (GeV/c)',3x,'py (GeV/c)',3x, 
     &'pz (GeV/c)',16x,'I') 
 1700 FORMAT(1x,'I',15x,a8,3(2x,f10.3,1x),15x,'I')  
 1800 FORMAT(1x,'Error: unrecognized coordinate frame ''',a,'''.'/  
     &1x,'Execution stopped!')  
 1900 FORMAT(1x,'Error: too low CM energy,',f8.3,' GeV for event ', 
     &'generation.'/1x,'Execution stopped!')    
    
      RETURN    
      END