Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pyinit_hijing.F
1 * $Id$
2       SUBROUTINE PYINIT_HIJING(FRAME,BEAM,TARGET,WIN)  
3     
4 C...Initializes the generation procedure; finds maxima of the   
5 C...differential cross-sections to be used for weighting.   
6 #include "ludat1_hijing.inc"
7 #include "ludat2_hijing.inc"
8 #include "ludat3_hijing.inc"
9 #include "ludat4_hijing.inc"
10 #include "pysubs_hijing.inc"
11 #include "pypars_hijing.inc"
12 #include "pyint1_hijing.inc"
13 #include "pyint2_hijing.inc"
14 #include "pyint5_hijing.inc"
15       CHARACTER*(*) FRAME,BEAM,TARGET   
16       CHARACTER CHFRAM*8,CHBEAM*8,CHTARG*8,CHMO(12)*3,CHLH(2)*6 
17       DATA CHMO/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep',  
18      &'Oct','Nov','Dec'/, CHLH/'lepton','hadron'/   
19     
20 C...Write headers.  
21 C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1000) MSTP(181),MSTP(182),  
22 C     &MSTP(185),CHMO(MSTP(184)),MSTP(183)   
23       CALL LULIST_HIJING(0)
24 C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1100)  
25     
26 C...Identify beam and target particles and initialize kinematics.   
27       CHFRAM=FRAME//' ' 
28       CHBEAM=BEAM//' '  
29       CHTARG=TARGET//' '    
30       CALL PYINKI_HIJING(CHFRAM,CHBEAM,CHTARG,WIN) 
31     
32 C...Select partonic subprocesses to be included in the simulation.  
33       IF(MSEL.NE.0) THEN    
34         DO 100 I=1,200  
35   100   MSUB(I)=0   
36       ENDIF 
37       IF(MINT(43).EQ.1.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN   
38 C...Lepton+lepton -> gamma/Z0 or W. 
39         IF(MINT(11)+MINT(12).EQ.0) MSUB(1)=1    
40         IF(MINT(11)+MINT(12).NE.0) MSUB(2)=1    
41       ELSEIF(MSEL.EQ.1) THEN    
42 C...High-pT QCD processes:  
43         MSUB(11)=1  
44         MSUB(12)=1  
45         MSUB(13)=1  
46         MSUB(28)=1  
47         MSUB(53)=1  
48         MSUB(68)=1  
49         IF(MSTP(82).LE.1.AND.CKIN(3).LT.PARP(81)) MSUB(95)=1    
50         IF(MSTP(82).GE.2.AND.CKIN(3).LT.PARP(82)) MSUB(95)=1    
51       ELSEIF(MSEL.EQ.2) THEN    
52 C...All QCD processes:  
53         MSUB(11)=1  
54         MSUB(12)=1  
55         MSUB(13)=1  
56         MSUB(28)=1  
57         MSUB(53)=1  
58         MSUB(68)=1  
59         MSUB(91)=1  
60         MSUB(92)=1  
61         MSUB(93)=1  
62         MSUB(95)=1  
63       ELSEIF(MSEL.GE.4.AND.MSEL.LE.8) THEN  
64 C...Heavy quark production. 
65         MSUB(81)=1  
66         MSUB(82)=1  
67         DO 110 J=1,MIN(8,MDCY(21,3))    
68   110   MDME(MDCY(21,2)+J-1,1)=0    
69         MDME(MDCY(21,2)+MSEL-1,1)=1 
70       ELSEIF(MSEL.EQ.10) THEN   
71 C...Prompt photon production:   
72         MSUB(14)=1  
73         MSUB(18)=1  
74         MSUB(29)=1  
75       ELSEIF(MSEL.EQ.11) THEN   
76 C...Z0/gamma* production:   
77         MSUB(1)=1   
78       ELSEIF(MSEL.EQ.12) THEN   
79 C...W+/- production:    
80         MSUB(2)=1   
81       ELSEIF(MSEL.EQ.13) THEN   
82 C...Z0 + jet:   
83         MSUB(15)=1  
84         MSUB(30)=1  
85       ELSEIF(MSEL.EQ.14) THEN   
86 C...W+/- + jet: 
87         MSUB(16)=1  
88         MSUB(31)=1  
89       ELSEIF(MSEL.EQ.15) THEN   
90 C...Z0 & W+/- pair production:  
91         MSUB(19)=1  
92         MSUB(20)=1  
93         MSUB(22)=1  
94         MSUB(23)=1  
95         MSUB(25)=1  
96       ELSEIF(MSEL.EQ.16) THEN   
97 C...H0 production:  
98         MSUB(3)=1   
99         MSUB(5)=1   
100         MSUB(8)=1   
101         MSUB(102)=1 
102       ELSEIF(MSEL.EQ.17) THEN   
103 C...H0 & Z0 or W+/- pair production:    
104         MSUB(24)=1  
105         MSUB(26)=1  
106       ELSEIF(MSEL.EQ.21) THEN   
107 C...Z'0 production: 
108         MSUB(141)=1 
109       ELSEIF(MSEL.EQ.22) THEN   
110 C...H+/- production:    
111         MSUB(142)=1 
112       ELSEIF(MSEL.EQ.23) THEN   
113 C...R production:   
114         MSUB(143)=1 
115       ENDIF 
116     
117 C...Count number of subprocesses on.    
118       MINT(44)=0    
119       DO 120 ISUB=1,200 
120       IF(MINT(43).LT.4.AND.ISUB.GE.91.AND.ISUB.LE.96.AND.   
121      &MSUB(ISUB).EQ.1) THEN 
122         WRITE(MSTU(11),1200) ISUB,CHLH(MINT(41)),CHLH(MINT(42)) 
123         STOP    
124       ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).EQ.-1) THEN 
125         WRITE(MSTU(11),1300) ISUB   
126         STOP    
127       ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).LE.-2) THEN 
128         WRITE(MSTU(11),1400) ISUB   
129         STOP    
130       ELSEIF(MSUB(ISUB).EQ.1) THEN  
131         MINT(44)=MINT(44)+1 
132       ENDIF 
133   120 CONTINUE  
134       IF(MINT(44).EQ.0) THEN    
135         WRITE(MSTU(11),1500)    
136         STOP    
137       ENDIF 
138       MINT(45)=MINT(44)-MSUB(91)-MSUB(92)-MSUB(93)-MSUB(94) 
139     
140 C...Maximum 4 generations; set maximum number of allowed flavours.  
141       MSTP(1)=MIN(4,MSTP(1))    
142       MSTU(114)=MIN(MSTU(114),2*MSTP(1))    
143       MSTP(54)=MIN(MSTP(54),2*MSTP(1))  
144     
145 C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton. 
146       DO 140 I=-20,20   
147       VINT(180+I)=0.    
148       IA=IABS(I)    
149       IF(IA.GE.1.AND.IA.LE.2*MSTP(1)) THEN  
150         DO 130 J=1,MSTP(1)  
151         IB=2*J-1+MOD(IA,2)  
152         IPM=(5-ISIGN(1,I))/2    
153         IDC=J+MDCY(IA,2)+2  
154   130   IF(MDME(IDC,1).EQ.1.OR.MDME(IDC,1).EQ.IPM) VINT(180+I)= 
155      &  VINT(180+I)+VCKM((IA+1)/2,(IB+1)/2) 
156       ELSEIF(IA.GE.11.AND.IA.LE.10+2*MSTP(1)) THEN  
157         VINT(180+I)=1.  
158       ENDIF 
159   140 CONTINUE  
160     
161 C...Choose Lambda value to use in alpha-strong. 
162       MSTU(111)=MSTP(2) 
163       IF(MSTP(3).GE.1) THEN 
164         ALAM=PARP(1)    
165         IF(MSTP(51).EQ.1) ALAM=0.2  
166         IF(MSTP(51).EQ.2) ALAM=0.29 
167         IF(MSTP(51).EQ.3) ALAM=0.2  
168         IF(MSTP(51).EQ.4) ALAM=0.4  
169         IF(MSTP(51).EQ.11) ALAM=0.16    
170         IF(MSTP(51).EQ.12) ALAM=0.26    
171         IF(MSTP(51).EQ.13) ALAM=0.36    
172         PARP(1)=ALAM    
173         PARP(61)=ALAM   
174         PARU(112)=ALAM  
175         PARJ(81)=ALAM   
176       ENDIF 
177     
178 C...Initialize widths and partial widths for resonances.    
179       CALL PYINRE_HIJING   
180     
181 C...Reset variables for cross-section calculation.  
182       DO 150 I=0,200    
183       DO 150 J=1,3  
184       NGEN(I,J)=0   
185   150 XSEC(I,J)=0.  
186       VINT(108)=0.  
187     
188 C...Find parametrized total cross-sections. 
189       IF(MINT(43).EQ.4) CALL PYXTOT_HIJING 
190     
191 C...Maxima of differential cross-sections.  
192       IF(MSTP(121).LE.0) CALL PYMAXI_HIJING    
193     
194 C...Initialize possibility of overlayed events. 
195       IF(MSTP(131).NE.0) CALL PYOVLY_HIJING(1) 
196     
197 C...Initialize multiple interactions with variable impact parameter.    
198       IF(MINT(43).EQ.4.AND.(MINT(45).NE.0.OR.MSTP(131).NE.0).AND.   
199      &MSTP(82).GE.2) CALL PYMULT_HIJING(1) 
200 C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1600)  
201     
202 C...Formats for initialization information. 
203  1000 FORMAT(///20X,'The Lund Monte Carlo - PYTHIA version ',I1,'.',I1/ 
204      &20X,'**  Last date of change:  ',I2,1X,A3,1X,I4,'  **'/)  
205  1100 FORMAT('1',18('*'),1X,'PYINIT_HIJING: initialization of PYTHIA ',    
206      &'routines',1X,17('*'))    
207  1200 FORMAT(1X,'Error: process number ',I3,' not meaningful for ',A6,  
208      &'-',A6,' interactions.'/1X,'Execution stopped!')  
209  1300 FORMAT(1X,'Error: requested subprocess',I4,' not implemented.'/   
210      &1X,'Execution stopped!')  
211  1400 FORMAT(1X,'Error: requested subprocess',I4,' not existing.'/  
212      &1X,'Execution stopped!')  
213  1500 FORMAT(1X,'Error: no subprocess switched on.'/    
214      &1X,'Execution stopped.')  
215  1600 FORMAT(/1X,22('*'),1X,'PYINIT_HIJING: initialization completed',1X
216      $     ,22('*'))  
217     
218       RETURN    
219       END