Removing MC position of the primary vertex (Jouri)

[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / luonia_hijing.F

* $Id$
    
C*********************************************************************  
    
      SUBROUTINE LUONIA_HIJING(KFL,ECM)    
    
C...Purpose: to generate Upsilon and toponium decays into three 
C...gluons or two gluons and a photon.  
#include "lujets_hijing.inc"
#include "ludat1_hijing.inc"
#include "ludat2_hijing.inc"
    
C...Printout. Check input parameters.   
      IF(MSTU(12).GE.1) CALL LULIST_HIJING(0)  
      IF(KFL.LT.0.OR.KFL.GT.8) THEN 
         CALL LUERRM_HIJING(16
     $        ,'(LUONIA_HIJING:) called with unknown flavour code')    
        IF(MSTU(21).GE.1) RETURN    
      ENDIF 
      IF(ECM.LT.PARJ(127)+2.02*PARF(101)) THEN  
         CALL LUERRM_HIJING(16
     $        ,'(LUONIA_HIJING:) called with too small CM energy') 
        IF(MSTU(21).GE.1) RETURN    
      ENDIF 
    
C...Initial e+e- and onium state (optional).    
      NC=0  
      IF(MSTJ(115).GE.2) THEN   
        NC=NC+2 
        CALL LU1ENT_HIJING(NC-1,11,0.5*ECM,0.,0.)  
        K(NC-1,1)=21    
        CALL LU1ENT_HIJING(NC,-11,0.5*ECM,PARU(1),0.)  
        K(NC,1)=21  
      ENDIF 
      KFLC=IABS(KFL)    
      IF(MSTJ(115).GE.3.AND.KFLC.GE.5) THEN 
        NC=NC+1 
        KF=110*KFLC+3   
        MSTU10=MSTU(10) 
        MSTU(10)=1  
        P(NC,5)=ECM 
        CALL LU1ENT_HIJING(NC,KF,ECM,0.,0.)    
        K(NC,1)=21  
        K(NC,3)=1   
        MSTU(10)=MSTU10 
      ENDIF 
    
C...Choose x1 and x2 according to matrix element.   
      NTRY=0    
  100 X1=RLU_HIJING(0) 
      X2=RLU_HIJING(0) 
      X3=2.-X1-X2   
      IF(X3.GE.1..OR.((1.-X1)/(X2*X3))**2+((1.-X2)/(X1*X3))**2+ 
     &((1.-X3)/(X1*X2))**2.LE.2.*RLU_HIJING(0)) GOTO 100   
      NTRY=NTRY+1   
      NJET=3    
      IF(MSTJ(101).LE.4) CALL LU3ENT_HIJING(NC+1,21,21,21,ECM,X1,X3)   
      IF(MSTJ(101).GE.5) CALL LU3ENT_HIJING(-(NC+1),21,21,21,ECM,X1,X3)    
    
C...Photon-gluon-gluon events. Small system modifications. Jet origin.  
      MSTU(111)=MSTJ(108)   
      IF(MSTJ(108).EQ.2.AND.(MSTJ(101).EQ.0.OR.MSTJ(101).EQ.1)) 
     &MSTU(111)=1   
      PARU(112)=PARJ(121)   
      IF(MSTU(111).EQ.2) PARU(112)=PARJ(122)    
      QF=0. 
      IF(KFLC.NE.0) QF=KCHG(KFLC,1)/3.  
      RGAM=7.2*QF**2*PARU(101)/ULALPS_HIJING(ECM**2)   
      MK=0  
      ECMC=ECM  
      IF(RLU_HIJING(0).GT.RGAM/(1.+RGAM)) THEN 
        IF(1.-MAX(X1,X2,X3).LE.MAX((PARJ(126)/ECM)**2,PARJ(125)))   
     &  NJET=2  
        IF(NJET.EQ.2.AND.MSTJ(101).LE.4) CALL LU2ENT_HIJING(NC+1,21,21
     $       ,ECM)    
        IF(NJET.EQ.2.AND.MSTJ(101).GE.5) CALL LU2ENT_HIJING(-(NC+1),21
     $       ,21,ECM) 
      ELSE  
        MK=1    
        ECMC=SQRT(1.-X1)*ECM    
        IF(ECMC.LT.2.*PARJ(127)) GOTO 100   
        K(NC+1,1)=1 
        K(NC+1,2)=22    
        K(NC+1,4)=0 
        K(NC+1,5)=0 
        IF(MSTJ(101).GE.5) K(NC+2,4)=MSTU(5)*(NC+3) 
        IF(MSTJ(101).GE.5) K(NC+2,5)=MSTU(5)*(NC+3) 
        IF(MSTJ(101).GE.5) K(NC+3,4)=MSTU(5)*(NC+2) 
        IF(MSTJ(101).GE.5) K(NC+3,5)=MSTU(5)*(NC+2) 
        NJET=2  
        IF(ECMC.LT.4.*PARJ(127)) THEN   
          MSTU10=MSTU(10)   
          MSTU(10)=1    
          P(NC+2,5)=ECMC    
          CALL LU1ENT_HIJING(NC+2,83,0.5*(X2+X3)*ECM,PARU(1),0.)   
          MSTU(10)=MSTU10   
          NJET=0    
        ENDIF   
      ENDIF 
      DO 110 IP=NC+1,N  
  110 K(IP,3)=K(IP,3)+(MSTJ(115)/2)+(KFLC/5)*(MSTJ(115)/3)*(NC-1)   
    
C...Differential cross-sections. Upper limit for cross-section. 
      IF(MSTJ(106).EQ.1) THEN   
        SQ2=SQRT(2.)    
        HF1=1.-PARJ(131)*PARJ(132)  
        HF3=PARJ(133)**2    
        CT13=(X1*X3-2.*X1-2.*X3+2.)/(X1*X3) 
        ST13=SQRT(1.-CT13**2)   
        SIGL=0.5*X3**2*((1.-X2)**2+(1.-X3)**2)*ST13**2  
        SIGU=(X1*(1.-X1))**2+(X2*(1.-X2))**2+(X3*(1.-X3))**2-SIGL   
        SIGT=0.5*SIGL   
        SIGI=(SIGL*CT13/ST13+0.5*X1*X3*(1.-X2)**2*ST13)/SQ2 
        SIGMAX=(2.*HF1+HF3)*ABS(SIGU)+2.*(HF1+HF3)*ABS(SIGL)+2.*(HF1+   
     &  2.*HF3)*ABS(SIGT)+2.*SQ2*(HF1+2.*HF3)*ABS(SIGI) 
    
C...Angular orientation of event.   
  120   CHI=PARU(2)*RLU_HIJING(0)  
        CTHE=2.*RLU_HIJING(0)-1.   
        PHI=PARU(2)*RLU_HIJING(0)  
        CCHI=COS(CHI)   
        SCHI=SIN(CHI)   
        C2CHI=COS(2.*CHI)   
        S2CHI=SIN(2.*CHI)   
        THE=ACOS(CTHE)  
        STHE=SIN(THE)   
        C2PHI=COS(2.*(PHI-PARJ(134)))   
        S2PHI=SIN(2.*(PHI-PARJ(134)))   
        SIG=((1.+CTHE**2)*HF1+STHE**2*C2PHI*HF3)*SIGU+2.*(STHE**2*HF1-  
     &  STHE**2*C2PHI*HF3)*SIGL+2.*(STHE**2*C2CHI*HF1+((1.+CTHE**2)*    
     &  C2CHI*C2PHI-2.*CTHE*S2CHI*S2PHI)*HF3)*SIGT-2.*SQ2*(2.*STHE*CTHE*    
     &  CCHI*HF1-2.*STHE*(CTHE*CCHI*C2PHI-SCHI*S2PHI)*HF3)*SIGI 
        IF(SIG.LT.SIGMAX*RLU_HIJING(0)) GOTO 120   
        CALL LUDBRB_HIJING(NC+1,N,0.,CHI,0D0,0D0,0D0)  
        CALL LUDBRB_HIJING(NC+1,N,THE,PHI,0D0,0D0,0D0) 
      ENDIF 
    
C...Generate parton shower. Rearrange along strings and check.  
      IF(MSTJ(101).GE.5.AND.NJET.GE.2) THEN 
        CALL LUSHOW_HIJING(NC+MK+1,-NJET,ECMC) 
        MSTJ14=MSTJ(14) 
        IF(MSTJ(105).EQ.-1) MSTJ(14)=0  
        IF(MSTJ(105).GE.0) MSTU(28)=0   
        CALL LUPREP_HIJING(0)  
        MSTJ(14)=MSTJ14 
        IF(MSTJ(105).GE.0.AND.MSTU(28).NE.0) GOTO 100   
      ENDIF 
    
C...Generate fragmentation. Information for LUTABU_HIJING: 
      IF(MSTJ(105).EQ.1) CALL LUEXEC_HIJING    
      MSTU(161)=110*KFLC+3  
      MSTU(162)=0   
    
      RETURN    
      END