luboei.f

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

    
C*********************************************************************  
    
      SUBROUTINE luboei(NSAV)   
    
C...Purpose: to modify event so as to approximately take into account   
C...Bose-Einstein effects according to a simple phenomenological    
C...parametrization.    
      IMPLICIT DOUBLE PRECISION(d)  
      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/ 
      dimension dps(4),kfbe(9),nbe(0:9),bei(100)    
      DATA kfbe/211,-211,111,321,-321,130,310,221,331/  
    
C...Boost event to overall CM frame. Calculate CM energy.   
      IF((mstj(51).NE.1.AND.mstj(51).NE.2).OR.n-nsav.LE.1) RETURN   
      DO 100 j=1,4  
  100 dps(j)=0. 
      DO 120 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 120  
      DO 110 j=1,4  
  110 dps(j)=dps(j)+p(i,j)  
  120 CONTINUE  
      CALL ludbrb(0,0,0.,0.,-dps(1)/dps(4),-dps(2)/dps(4),  
     &-dps(3)/dps(4))   
      pecm=0.   
      DO 130 i=1,n  
  130 IF(k(i,1).GE.1.AND.k(i,1).LE.10) pecm=pecm+p(i,4) 
    
C...Reserve copy of particles by species at end of record.  
      nbe(0)=n+mstu(3)  
      DO 160 ibe=1,min(9,mstj(51))  
      nbe(ibe)=nbe(ibe-1)   
      DO 150 i=nsav+1,n 
      IF(k(i,2).NE.kfbe(ibe)) goto 150  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 150  
      IF(nbe(ibe).GE.mstu(4)-mstu(32)-5) THEN   
        CALL luerrm(11,'(LUBOEI:) no more memory left in LUJETS')   
        RETURN  
      ENDIF 
      nbe(ibe)=nbe(ibe)+1   
      k(nbe(ibe),1)=i   
      DO 140 j=1,3  
  140 p(nbe(ibe),j)=0.  
  150 CONTINUE  
  160 CONTINUE  
    
C...Tabulate integral for subsequent momentum shift.    
      DO 210 ibe=1,min(9,mstj(51))  
      IF(ibe.NE.1.AND.ibe.NE.4.AND.ibe.LE.7) goto 180   
      IF(ibe.EQ.1.AND.max(nbe(1)-nbe(0),nbe(2)-nbe(1),nbe(3)-nbe(2)).   
     &le.1) goto 180    
      IF(ibe.EQ.4.AND.max(nbe(4)-nbe(3),nbe(5)-nbe(4),nbe(6)-nbe(5),    
     &nbe(7)-nbe(6)).LE.1) goto 180 
      IF(ibe.GE.8.AND.nbe(ibe)-nbe(ibe-1).LE.1) goto 180    
      IF(ibe.EQ.1) pmhq=2.*ulmass(211)  
      IF(ibe.EQ.4) pmhq=2.*ulmass(321)  
      IF(ibe.EQ.8) pmhq=2.*ulmass(221)  
      IF(ibe.EQ.9) pmhq=2.*ulmass(331)  
      qdel=0.1*min(pmhq,parj(93))   
      IF(mstj(51).EQ.1) THEN    
        nbin=min(100,nint(9.*parj(93)/qdel))    
        beex=exp(0.5*qdel/parj(93)) 
        bert=exp(-qdel/parj(93))    
      ELSE  
        nbin=min(100,nint(3.*parj(93)/qdel))    
      ENDIF 
      DO 170 ibin=1,nbin    
      qbin=qdel*(ibin-0.5)  
      bei(ibin)=qdel*(qbin**2+qdel**2/12.)/sqrt(qbin**2+pmhq**2)    
      IF(mstj(51).EQ.1) THEN    
        beex=beex*bert  
        bei(ibin)=bei(ibin)*beex    
      ELSE  
        bei(ibin)=bei(ibin)*exp(-(qbin/parj(93))**2)    
      ENDIF 
  170 IF(ibin.GE.2) bei(ibin)=bei(ibin)+bei(ibin-1) 
    
C...Loop through particle pairs and find old relative momentum. 
  180 DO 200 i1m=nbe(ibe-1)+1,nbe(ibe)-1    
      i1=k(i1m,1)   
      DO 200 i2m=i1m+1,nbe(ibe) 
      i2=k(i2m,1)   
      q2old=max(0.,(p(i1,4)+p(i2,4))**2-(p(i1,1)+p(i2,1))**2-(p(i1,2)+  
     &p(i2,2))**2-(p(i1,3)+p(i2,3))**2-(p(i1,5)+p(i2,5))**2)    
      qold=sqrt(q2old)  
    
C...Calculate new relative momentum.    
      IF(qold.LT.0.5*qdel) THEN 
        qmov=qold/3.    
      ELSEIF(qold.LT.(nbin-0.1)*qdel) THEN  
        rbin=qold/qdel  
        ibin=rbin   
        rinp=(rbin**3-ibin**3)/(3*ibin*(ibin+1)+1)  
        qmov=(bei(ibin)+rinp*(bei(ibin+1)-bei(ibin)))*  
     &  sqrt(q2old+pmhq**2)/q2old   
      ELSE  
        qmov=bei(nbin)*sqrt(q2old+pmhq**2)/q2old    
      ENDIF 
      q2new=q2old*(qold/(qold+3.*parj(92)*qmov))**(2./3.)   
    
C...Calculate and save shift to be performed on three-momenta.  
      hc1=(p(i1,4)+p(i2,4))**2-(q2old-q2new)    
      hc2=(q2old-q2new)*(p(i1,4)-p(i2,4))**2    
      ha=0.5*(1.-sqrt(hc1*q2new/(hc1*q2old-hc2)))   
      DO 190 j=1,3  
      pd=ha*(p(i2,j)-p(i1,j))   
      p(i1m,j)=p(i1m,j)+pd  
  190 p(i2m,j)=p(i2m,j)-pd  
  200 CONTINUE  
  210 CONTINUE  
    
C...Shift momenta and recalculate energies. 
      DO 230 im=nbe(0)+1,nbe(min(9,mstj(51)))   
      i=k(im,1) 
      DO 220 j=1,3  
  220 p(i,j)=p(i,j)+p(im,j) 
  230 p(i,4)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)  
    
C...Rescale all momenta for energy conservation.    
      pes=0.    
      pqs=0.    
      DO 240 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 240  
      pes=pes+p(i,4)    
      pqs=pqs+p(i,5)**2/p(i,4)  
  240 CONTINUE  
      fac=(pecm-pqs)/(pes-pqs)  
      DO 260 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 260  
      DO 250 j=1,3  
  250 p(i,j)=fac*p(i,j) 
      p(i,4)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)  
  260 CONTINUE  
    
C...Boost back to correct reference frame.  
      CALL ludbrb(0,0,0.,0.,dps(1)/dps(4),dps(2)/dps(4),dps(3)/dps(4))  
    
      RETURN    
      END