pyhikmap.f

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

    
C*********************************************************************  
    
      SUBROUTINE pyhikmap(IVAR,MVAR,VVAR) 
    
C...Maps a uniform distribution into a distribution of a kinematical    
C...variable according to one of the possibilities allowed. It is   
C...assumed that kinematical limits have been set by a PYKLIM call. 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
      common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2) 
      SAVE /pyhiint2/ 
    
C...Convert VVAR to tau variable.   
      isub=mint(1)  
      IF(ivar.EQ.1) THEN    
        taumin=vint(11) 
        taumax=vint(31) 
        IF(mvar.EQ.3.OR.mvar.EQ.4) THEN 
          taure=vint(73)    
          gamre=vint(74)    
        ELSEIF(mvar.EQ.5.OR.mvar.EQ.6) THEN 
          taure=vint(75)    
          gamre=vint(76)    
        ENDIF   
        IF(mint(43).EQ.1.AND.(iset(isub).EQ.1.OR.iset(isub).EQ.2)) THEN 
          tau=1.    
        ELSEIF(mvar.EQ.1) THEN  
          tau=taumin*(taumax/taumin)**vvar  
        ELSEIF(mvar.EQ.2) THEN  
          tau=taumax*taumin/(taumin+(taumax-taumin)*vvar)   
        ELSEIF(mvar.EQ.3.OR.mvar.EQ.5) THEN 
          ratgen=(taure+taumax)/(taure+taumin)*taumin/taumax    
          tau=taure*taumin/((taure+taumin)*ratgen**vvar-taumin) 
        ELSE    
          aupp=atan((taumax-taure)/gamre)   
          alow=atan((taumin-taure)/gamre)   
          tau=taure+gamre*tan(alow+(aupp-alow)*vvar)    
        ENDIF   
        vint(21)=min(taumax,max(taumin,tau))    
    
C...Convert VVAR to y* variable.    
      ELSEIF(ivar.EQ.2) THEN    
        ystmin=vint(12) 
        ystmax=vint(32) 
        IF(mint(43).EQ.1) THEN  
          yst=0.    
        ELSEIF(mint(43).EQ.2) THEN  
          IF(iset(isub).LE.2) yst=-0.5*log(vint(21))    
          IF(iset(isub).GE.3) yst=-0.5*log(vint(26))    
        ELSEIF(mint(43).EQ.3) THEN  
          IF(iset(isub).LE.2) yst=0.5*log(vint(21)) 
          IF(iset(isub).GE.3) yst=0.5*log(vint(26)) 
        ELSEIF(mvar.EQ.1) THEN  
          yst=ystmin+(ystmax-ystmin)*sqrt(vvar) 
        ELSEIF(mvar.EQ.2) THEN  
          yst=ystmax-(ystmax-ystmin)*sqrt(1.-vvar)  
        ELSE    
          aupp=atan(exp(ystmax))    
          alow=atan(exp(ystmin))    
          yst=log(tan(alow+(aupp-alow)*vvar))   
        ENDIF   
        vint(22)=min(ystmax,max(ystmin,yst))    
    
C...Convert VVAR to cos(theta-hat) variable.    
      ELSEIF(ivar.EQ.3) THEN    
        rm34=2.*vint(63)*vint(64)/(vint(21)*vint(2))**2 
        rsqm=1.+rm34    
        IF(2.*vint(71)**2/(vint(21)*vint(2)).LT.0.0001) rm34=max(rm34,  
     &  2.*vint(71)**2/(vint(21)*vint(2)))  
        ctnmin=vint(13) 
        ctnmax=vint(33) 
        ctpmin=vint(14) 
        ctpmax=vint(34) 
        IF(mvar.EQ.1) THEN  
          aneg=ctnmax-ctnmin    
          apos=ctpmax-ctpmin    
          IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg) THEN  
            vctn=vvar*(aneg+apos)/aneg  
            cth=ctnmin+(ctnmax-ctnmin)*vctn 
          ELSE  
            vctp=(vvar*(aneg+apos)-aneg)/apos   
            cth=ctpmin+(ctpmax-ctpmin)*vctp 
          ENDIF 
        ELSEIF(mvar.EQ.2) THEN  
          rmnmin=max(rm34,rsqm-ctnmin)  
          rmnmax=max(rm34,rsqm-ctnmax)  
          rmpmin=max(rm34,rsqm-ctpmin)  
          rmpmax=max(rm34,rsqm-ctpmax)  
          aneg=log(rmnmin/rmnmax)   
          apos=log(rmpmin/rmpmax)   
          IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg) THEN  
            vctn=vvar*(aneg+apos)/aneg  
            cth=rsqm-rmnmin*(rmnmax/rmnmin)**vctn   
          ELSE  
            vctp=(vvar*(aneg+apos)-aneg)/apos   
            cth=rsqm-rmpmin*(rmpmax/rmpmin)**vctp   
          ENDIF 
        ELSEIF(mvar.EQ.3) THEN  
          rmnmin=max(rm34,rsqm+ctnmin)  
          rmnmax=max(rm34,rsqm+ctnmax)  
          rmpmin=max(rm34,rsqm+ctpmin)  
          rmpmax=max(rm34,rsqm+ctpmax)  
          aneg=log(rmnmax/rmnmin)   
          apos=log(rmpmax/rmpmin)   
          IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg) THEN  
            vctn=vvar*(aneg+apos)/aneg  
            cth=rmnmin*(rmnmax/rmnmin)**vctn-rsqm   
          ELSE  
            vctp=(vvar*(aneg+apos)-aneg)/apos   
            cth=rmpmin*(rmpmax/rmpmin)**vctp-rsqm   
          ENDIF 
        ELSEIF(mvar.EQ.4) THEN  
          rmnmin=max(rm34,rsqm-ctnmin)  
          rmnmax=max(rm34,rsqm-ctnmax)  
          rmpmin=max(rm34,rsqm-ctpmin)  
          rmpmax=max(rm34,rsqm-ctpmax)  
          aneg=1./rmnmax-1./rmnmin  
          apos=1./rmpmax-1./rmpmin  
          IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg) THEN  
            vctn=vvar*(aneg+apos)/aneg  
            cth=rsqm-1./(1./rmnmin+aneg*vctn)   
          ELSE  
            vctp=(vvar*(aneg+apos)-aneg)/apos   
            cth=rsqm-1./(1./rmpmin+apos*vctp)   
          ENDIF 
        ELSEIF(mvar.EQ.5) THEN  
          rmnmin=max(rm34,rsqm+ctnmin)  
          rmnmax=max(rm34,rsqm+ctnmax)  
          rmpmin=max(rm34,rsqm+ctpmin)  
          rmpmax=max(rm34,rsqm+ctpmax)  
          aneg=1./rmnmin-1./rmnmax  
          apos=1./rmpmin-1./rmpmax  
          IF(aneg.GT.0..AND.vvar*(aneg+apos).LE.aneg) THEN  
            vctn=vvar*(aneg+apos)/aneg  
            cth=1./(1./rmnmin-aneg*vctn)-rsqm   
          ELSE  
            vctp=(vvar*(aneg+apos)-aneg)/apos   
            cth=1./(1./rmpmin-apos*vctp)-rsqm   
          ENDIF 
        ENDIF   
        IF(cth.LT.0.) cth=min(ctnmax,max(ctnmin,cth))   
        IF(cth.GT.0.) cth=min(ctpmax,max(ctpmin,cth))   
        vint(23)=cth    
    
C...Convert VVAR to tau' variable.  
      ELSEIF(ivar.EQ.4) THEN    
        tau=vint(11)    
        taupmn=vint(16) 
        taupmx=vint(36) 
        IF(mint(43).EQ.1) THEN  
          taup=1.   
        ELSEIF(mvar.EQ.1) THEN  
          taup=taupmn*(taupmx/taupmn)**vvar 
        ELSE    
          aupp=(1.-tau/taupmx)**4   
          alow=(1.-tau/taupmn)**4   
          taup=tau/(1.-(alow+(aupp-alow)*vvar)**0.25)   
        ENDIF   
        vint(26)=min(taupmx,max(taupmn,taup))   
      ENDIF 
    
      RETURN    
      END