Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / luexec_hijing.F
1 * $Id$
2     
3 C*********************************************************************  
4     
5       SUBROUTINE LUEXEC_HIJING 
6     
7 C...Purpose: to administrate the fragmentation and decay chain. 
8 #include "lujets_hijing.inc"
9 #include "ludat1_hijing.inc"
10 #include "ludat2_hijing.inc"
11 #include "ludat3_hijing.inc"
12       DIMENSION PS(2,6) 
13     
14 C...Initialize and reset.   
15       MSTU(24)=0    
16       IF(MSTU(12).GE.1) CALL LULIST_HIJING(0)  
17       MSTU(31)=MSTU(31)+1   
18       MSTU(1)=0 
19       MSTU(2)=0 
20       MSTU(3)=0 
21       MCONS=1   
22     
23 C...Sum up momentum, energy and charge for starting entries.    
24       NSAV=N    
25       DO 100 I=1,2  
26       DO 100 J=1,6  
27   100 PS(I,J)=0.    
28       DO 120 I=1,N  
29       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120  
30       DO 110 J=1,4  
31   110 PS(1,J)=PS(1,J)+P(I,J)    
32       PS(1,6)=PS(1,6)+LUCHGE_HIJING(K(I,2))    
33   120 CONTINUE  
34       PARU(21)=PS(1,4)  
35     
36 C...Prepare system for subsequent fragmentation/decay.  
37       CALL LUPREP_HIJING(0)    
38     
39 C...Loop through jet fragmentation and particle decays. 
40       MBE=0 
41   130 MBE=MBE+1 
42       IP=0  
43   140 IP=IP+1   
44       KC=0  
45       IF(K(IP,1).GT.0.AND.K(IP,1).LE.10) KC=LUCOMP_HIJING(K(IP,2)) 
46       IF(KC.EQ.0) THEN  
47     
48 C...Particle decay if unstable and allowed. Save long-lived particle    
49 C...decays until second pass after Bose-Einstein effects.   
50       ELSEIF(KCHG(KC,2).EQ.0) THEN  
51         IF(MSTJ(21).GE.1.AND.MDCY(KC,1).GE.1.AND.(MSTJ(51).LE.0.OR.MBE. 
52      &  EQ.2.OR.PMAS(KC,2).GE.PARJ(91).OR.IABS(K(IP,2)).EQ.311))    
53      &  CALL LUDECY_HIJING(IP) 
54     
55 C...Decay products may develop a shower.    
56         IF(MSTJ(92).GT.0) THEN  
57           IP1=MSTJ(92)  
58           QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,  
59      &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))    
60           CALL LUSHOW_HIJING(IP1,IP1+1,QMAX)   
61           CALL LUPREP_HIJING(IP1)  
62           MSTJ(92)=0    
63         ELSEIF(MSTJ(92).LT.0) THEN  
64           IP1=-MSTJ(92) 
65           CALL LUSHOW_HIJING(IP1,-3,P(IP,5))   
66           CALL LUPREP_HIJING(IP1)  
67           MSTJ(92)=0    
68         ENDIF   
69     
70 C...Jet fragmentation: string or independent fragmentation. 
71       ELSEIF(K(IP,1).EQ.1.OR.K(IP,1).EQ.2) THEN 
72         MFRAG=MSTJ(1)   
73         IF(MFRAG.GE.1.AND.K(IP,1).EQ.1) MFRAG=2 
74         IF(MSTJ(21).GE.2.AND.K(IP,1).EQ.2.AND.N.GT.IP) THEN 
75           IF(K(IP+1,1).EQ.1.AND.K(IP+1,3).EQ.K(IP,3).AND.   
76      &    K(IP,3).GT.0.AND.K(IP,3).LT.IP) THEN  
77              IF(KCHG(LUCOMP_HIJING(K(K(IP,3),2)),2).EQ.0) MFRAG=MIN(1
78      $            ,MFRAG)    
79           ENDIF 
80         ENDIF   
81         IF(MFRAG.EQ.1) CALL LUSTRF_HIJING(IP)  
82         IF(MFRAG.EQ.2) CALL LUINDF_HIJING(IP)  
83         IF(MFRAG.EQ.2.AND.K(IP,1).EQ.1) MCONS=0 
84         IF(MFRAG.EQ.2.AND.(MSTJ(3).LE.0.OR.MOD(MSTJ(3),5).EQ.0)) MCONS=0    
85       ENDIF 
86     
87 C...Loop back if enough space left in LUJETS_HIJING and no error abort.    
88       IF(MSTU(24).NE.0.AND.MSTU(21).GE.2) THEN  
89       ELSEIF(IP.LT.N.AND.N.LT.MSTU(4)-20-MSTU(32)) THEN 
90         GOTO 140    
91       ELSEIF(IP.LT.N) THEN  
92          CALL LUERRM_HIJING(11
93      $        ,'(LUEXEC_HIJING:) no more memory left in LUJETS_HIJING')   
94       ENDIF 
95     
96 C...Include simple Bose-Einstein effect parametrization if desired. 
97       IF(MBE.EQ.1.AND.MSTJ(51).GE.1) THEN   
98         CALL LUBOEI_HIJING(NSAV)   
99         GOTO 130    
100       ENDIF 
101     
102 C...Check that momentum, energy and charge were conserved.  
103       DO 160 I=1,N  
104       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 160  
105       DO 150 J=1,4  
106   150 PS(2,J)=PS(2,J)+P(I,J)    
107       PS(2,6)=PS(2,6)+LUCHGE_HIJING(K(I,2))    
108   160 CONTINUE  
109       PDEV=(ABS(PS(2,1)-PS(1,1))+ABS(PS(2,2)-PS(1,2))+ABS(PS(2,3)-  
110      &PS(1,3))+ABS(PS(2,4)-PS(1,4)))/(1.+ABS(PS(2,4))+ABS(PS(1,4))) 
111       IF(MCONS.EQ.1.AND.PDEV.GT.PARU(11)) CALL LUERRM_HIJING(15,   
112      &'(LUEXEC_HIJING:) four-momentum was not conserved')  
113       IF(MCONS.EQ.1.AND.ABS(PS(2,6)-PS(1,6)).GT.0.1) CALL
114      $     LUERRM_HIJING(15,'(LUEXEC_HIJING:) charge was not conserved')
115      $     
116     
117       RETURN    
118       END