Fix for #84624: Problem in TPC cluster map
[u/mrichter/AliRoot.git] / HERWIG / hwlhin.f
1 C Collects all of the Les Houches interface routines, plus utilities
2 C for colour codes
3 C
4 C----------------------------------------------------------------------
5       SUBROUTINE UPEVNT_GUP
6 C----------------------------------------------------------------------
7 C  Reads MC@NLO input files and fills Les Houches event common HEPEUP
8 C----------------------------------------------------------------------
9       INCLUDE 'herwig65.inc'
10 C---Les Houches Event Common Block
11       INTEGER MAXNUP
12       PARAMETER (MAXNUP=500)
13       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
14       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP,
15      & XMP2,XMA2,XMB2,BETA,VA,VB,SIGMA,DELTA,S2,XKA,XKB,PTF,E,PL
16       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
17      &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
18      &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
19      &              SPINUP(MAXNUP)
20       DOUBLE PRECISION PCM(5),PTR,XMTR,HWVDOT,HWULDO,PDB(5)
21       INTEGER I,J,IC,JPR,MQQ,NQQ,IUNIT,ISCALE,I1HPRO,IBOS,NP,IG,
22      & ILEP,ID,IA,IB,ICOL4(4,4),ICOL5(5,18),JJPROC,IVHVEC,IVHLEP,MUP
23       PARAMETER (IUNIT=61)
24       LOGICAL BOPRO,NODEC,REMIT
25       COMMON/NQQCOM/MQQ,NQQ
26       COMMON/VHLIN/IVHVEC,IVHLEP
27 C---Colour flows for heavy quark pair production
28       DATA ICOL4/
29      & 10,02,10,02,01,20,20,01,12,23,10,03,12,31,30,02/
30       DATA ICOL5/
31      & 10,02,13,30,02, 10,02,32,10,03,
32      & 10,21,30,20,03, 10,23,20,10,03,
33      & 01,20,23,30,01, 01,20,31,20,03,
34      & 01,23,03,20,01, 01,12,03,30,02,
35      & 12,20,30,10,03, 12,30,10,30,02,
36      & 12,03,02,10,03, 12,01,03,30,02,
37      & 12,23,14,40,03, 12,34,32,10,04,
38      & 12,23,43,10,04, 12,31,34,40,02,
39      & 12,34,14,30,02, 12,31,42,30,04/
40       IF (IERROR.NE.0) RETURN
41 C---READ AN EVENT
42       IF(NQQ.GE.MQQ)CALL HWWARN('UPEVNT',201,*999)
43       READ(IUNIT,901) I1HPRO,IC,NP
44       READ(IUNIT,902) (IDUP(I),I=1,NP)
45       READ(IUNIT,903) XWGTUP
46 C---Les Houches expects mean weight to be the cross section in pb
47       XWGTUP= XWGTUP*MQQ
48       READ(IUNIT,904) ((PUP(J,I),J=1,4),I=1,NP)
49       NQQ=NQQ+1
50 C---Input format is now (px,py,pz,m)
51       DO I=1,NP
52          E=SQRT(HWVDOT(4,PUP(1,I),PUP(1,I)))
53          PUP(5,I)=PUP(4,I)
54          PUP(4,I)=E
55       ENDDO
56       CALL HWVSUM(4,PUP(1,1),PUP(1,2),PCM)
57       CALL HWUMAS(PCM)
58 C---REMIT MEANS A REAL PARTON EMISSION OCCURRED
59       REMIT=PUP(4,3).NE.ZERO
60 C---NODEC MEANS DECAYS NOT YET DONE
61       NODEC=NP.EQ.5
62       NUP=NP
63 C---CHECK PROCESS CODE
64       JJPROC=MOD(ABS(IPROC),10000)
65       JPR=JJPROC/100
66       BOPRO=JPR.EQ.13.OR.JPR.EQ.14.OR.JPR.EQ.16.OR.JPR.EQ.36
67       IF (BOPRO) THEN
68 C----------------------------------------------------------------------
69 C   SINGLE GAUGE OR HIGGS BOSON PRODUCTION
70 C   B = Z/gamma, W or H (SM or any MSSM neutral Higgs)
71 C-----------------------------------------------------------------------
72 C I1HPRO IDENTIFIES THE PARTONIC SUBPROCESS, WITH THE FOLLOWING CONVENTIONS:
73 C   I1HPRO         PROCESS
74 C    401        q qbar -> g B
75 C    402        q g    -> q B
76 C    403        qbar q -> g B
77 C    404        qbar g -> qbar B
78 C    405        g q    -> q B
79 C    406        g qbar -> qbar B
80 C    407        g g    -> g B
81 C-----------------------------------------------------------------------
82 C---NODEC=.TRUE. FOR HIGGS AND UNDECAYED EW BOSON
83          NODEC=NP.EQ.4
84          IHPRO=I1HPRO-400
85          ISCALE=0
86          IF(JPR.EQ.16)ISCALE=2
87       ELSEIF (JPR.EQ.17.OR.JPR.EQ.20) THEN
88 C----------------------------------------------------------------------
89 C   HEAVY Q and/or QBAR PRODUCTION
90 C   IPROC=-1705,-1706 for Q=b,t
91 C   IPROC=-2000 for single top
92 C-----------------------------------------------------------------------
93 C I1HPRO IDENTIFIES THE PARTONIC SUBPROCESS, WITH THE FOLLOWING CONVENTIONS:
94 C   I1HPRO         PROCESS
95 C    401        q qbar -> g Q Qbar
96 C    402        q g    -> q Q Qbar
97 C    403        qbar q -> g Q Qbar
98 C    404        qbar g -> qbar Q Qbar
99 C    405        g q    -> q Q Qbar
100 C    406        g qbar -> qbar Q Qbar
101 C    407        g g    -> g Q Qbar
102 C    408        q q    -> g t q
103 C    409        qbar qbar -> g tbar qbar
104 C-----------------------------------------------------------------------
105 C IC SPECIFIES THE COLOUR CONNECTION (NOW IN INPUT FILE)
106 C-----------------------------------------------------------------------
107 C---SET IHPRO AS FOR DIRECT PHOTON (IPROC=1800)
108          IHPRO=I1HPRO-360
109          ISCALE=0
110          IF(ABS(IPROC).EQ.1705.OR.ABS(IPROC).EQ.11705)ISCALE=5
111       ELSEIF (JPR.EQ.28) THEN
112 C----------------------------------------------------------------------
113 C   GAUGE BOSON PAIR PRODUCTION
114 C   VV=WW,ZZ,ZW+,ZW- FOR IPROC=-2850,-2860,-2870,-2880
115 C-----------------------------------------------------------------------
116 C I1HPRO IDENTIFIES THE PARTONIC SUBPROCESS, WITH THE FOLLOWING CONVENTIONS:
117 C   I1HPRO         PROCESS
118 C    401        q qbar -> g V V
119 C    402        q g    -> q V V
120 C    403        qbar q -> g V V
121 C    404        qbar g -> qbar V V
122 C    405        g q    -> q V V
123 C    406        g qbar -> qbar V V
124 C-----------------------------------------------------------------------
125          IHPRO=I1HPRO-400
126          ISCALE=0
127       ELSEIF (JPR.EQ.26.OR.JPR.EQ.27) THEN
128 C----------------------------------------------------------------------
129 C   GAUGE BOSON PLUS HIGGS PRODUCTION
130 C   VH=WH,ZH FOR IPROC=-2600-ID,-2700-ID
131 C   WHERE ID CONTROLS HIGGS DECAY AS IN STANDARD HERWIG
132 C-----------------------------------------------------------------------
133          IHPRO=I1HPRO-400
134          ISCALE=0
135       ELSE
136          CALL HWWARN('UPEVNT',202,*999)
137       ENDIF
138 C---HARD SCALE
139       SCALUP=PCM(5)
140       IF (REMIT) THEN
141          IF (ISCALE.EQ.0) THEN
142             PTR=SQRT(PUP(1,3)**2+PUP(2,3)**2)
143             SCALUP=PCM(5)-2.*PTR
144          ELSEIF(ISCALE.EQ.1)THEN
145             SCALUP=PCM(5)
146          ELSEIF (ISCALE.EQ.2) THEN
147             SCALUP=SQRT(PUP(1,3)**2+PUP(2,3)**2)
148          ELSEIF (ISCALE.EQ.3.OR.ISCALE.EQ.4.OR.ISCALE.EQ.5) THEN
149             PTR=SQRT(PUP(1,3)**2+PUP(2,3)**2)
150             IA=4
151             IB=5
152             XMP2=PUP(5,3)**2
153             XMA2=PUP(5,IA)**2
154             XMB2=PUP(5,IB)**2
155             S2=XMA2+XMB2+2*HWULDO(PUP(1,IA),PUP(1,IB))
156             SIGMA=XMA2+XMB2
157             DELTA=XMA2-XMB2
158             BETA=SQRT(1-2*SIGMA/S2+(DELTA/S2)**2)
159             VA=BETA/(1+DELTA/S2)
160             VB=BETA/(1-DELTA/S2)
161             XKA=HWULDO(PUP(1,3),PUP(1,IA))
162             XKB=HWULDO(PUP(1,3),PUP(1,IB))
163             E=(XKA+XKB)/SQRT(S2)
164             PL=-2.0/((VA+VB)*BETA*SQRT(S2))*(VA*XKA-VB*XKB)
165             PTF=E**2-PL**2-XMP2
166             IF (PTF.LE.ZERO) CALL HWWARN('UPEVNT',103,*999)
167             PTF=SQRT(PTF)
168             IF(ISCALE.EQ.3)THEN
169               SCALUP=PCM(5)-2.*MIN(PTR,PTF)
170             ELSEIF(ISCALE.EQ.4)THEN
171               SCALUP=MIN(PTR,PTF)
172             ELSE
173               SCALUP=(MIN(PTR,PTF))**2+(XMA2+XMB2)/2.D0
174               SCALUP=SQRT(SCALUP)
175             ENDIF
176             IF (SCALUP.LE.ZERO) CALL HWWARN('UPEVNT',100,*999)
177          ELSEIF (ISCALE.EQ.6) THEN
178             XMTR=SQRT(PUP(5,4)**2+PUP(1,4)**2+PUP(2,4)**2)
179             PTR=SQRT(PUP(1,3)**2+PUP(2,3)**2)
180             SCALUP=PCM(5)-PTR-XMTR
181             IF (SCALUP.LE.ZERO) CALL HWWARN('UPEVNT',100,*999)
182          ELSEIF (ISCALE.EQ.7) THEN
183             SCALUP=SQRT(PUP(5,4)**2+PUP(1,4)**2+PUP(2,4)**2)
184          ELSE
185             CALL HWWARN('UPEVNT',501,*999)
186          ENDIF
187       ELSE
188          NUP=NUP-1
189       ENDIF
190 C---INITIAL STATE
191       DO I=1,2
192          ISTUP(I)=-1
193          MOTHUP(1,I)=0
194          MOTHUP(2,I)=0
195       ENDDO
196 C---FINAL STATE
197       DO I=3,NUP
198          ISTUP(I)=1
199          MOTHUP(1,I)=1
200          MOTHUP(2,I)=2
201       ENDDO
202       IF (BOPRO.AND.NODEC) THEN
203 C---SINGLE BOSON (UNDECAYED)
204          IF (REMIT) THEN
205 C---SET COLOUR CONNECTIONS
206             DO I=1,3
207                ICOLUP(1,I)=501
208                ICOLUP(2,I)=502
209             ENDDO
210             IF (IHPRO.EQ.1) THEN
211                ICOLUP(2,1)=0
212                ICOLUP(1,2)=0
213             ELSEIF (IHPRO.EQ.2) THEN
214                ICOLUP(1,1)=502
215                ICOLUP(2,1)=0
216                ICOLUP(2,3)=0
217             ELSEIF (IHPRO.EQ.3) THEN
218                ICOLUP(1,1)=0
219                ICOLUP(2,2)=0
220             ELSEIF (IHPRO.EQ.4) THEN
221                ICOLUP(1,1)=0
222                ICOLUP(2,1)=501
223                ICOLUP(1,3)=0
224             ELSEIF (IHPRO.EQ.5) THEN
225                ICOLUP(1,2)=502
226                ICOLUP(2,2)=0
227                ICOLUP(2,3)=0
228             ELSEIF (IHPRO.EQ.6) THEN
229                ICOLUP(1,2)=0
230                ICOLUP(2,2)=501
231                ICOLUP(1,3)=0
232             ELSEIF (IHPRO.EQ.7) THEN
233                ICOLUP(1,2)=502
234                ICOLUP(2,2)=503
235                ICOLUP(2,3)=503
236             ELSE
237                CALL HWWARN('UPEVNT',101,*999)
238             ENDIF
239          ELSE
240             CALL HWVEQU(5,PUP(1,4),PUP(1,3))
241 C---SET COLOUR CONNECTIONS
242             DO I=1,2
243                ICOLUP(1,I)=0
244                ICOLUP(2,I)=0
245             ENDDO
246             IF (IDUP(1).GT.0) THEN
247                ICOLUP(1,1)=501
248                ICOLUP(2,2)=501
249                IF (IDUP(1).GT.0) THEN
250 C---GG FUSION
251                   ICOLUP(2,1)=502
252                   ICOLUP(1,2)=502
253                ENDIF
254             ELSE
255 C---QBAR Q
256                ICOLUP(2,1)=501
257                ICOLUP(1,2)=501
258             ENDIF
259          ENDIF
260          ICOLUP(1,NUP)=0
261          ICOLUP(2,NUP)=0
262 C---LOAD BOSON ID
263          IF (JPR.EQ.13) THEN
264             IDUP(NUP)=23
265          ELSEIF (JPR.EQ.16) THEN
266             IDUP(NUP)=25
267          ELSEIF (JPR.EQ.36) THEN
268             IBOS=MOD(JJPROC,100)
269             IF (IBOS.EQ.10) THEN
270                IDUP(NUP)=26
271             ELSEIF (IBOS.EQ.20) THEN
272                IDUP(NUP)=35
273             ELSEIF (IBOS.EQ.30) THEN
274                IDUP(NUP)=36
275             ELSE
276                CALL HWWARN('UPEVNT',502,*999)
277             ENDIF
278          ELSEIF (JPR.EQ.14) THEN
279             IBOS=0
280             DO I=1,NUP-1
281                ID=IDUP(I)
282                IF (ID.EQ.21) THEN
283                   IC=0
284                ELSEIF (ID.GT.0) THEN
285                   IC=ICHRG(ID)
286                ELSE
287                   IC=ICHRG(6-ID)
288                ENDIF
289                IBOS=IBOS+IC
290             ENDDO
291             IF (REMIT) IBOS=IBOS-2*IC
292             IF (ABS(IBOS).NE.3) CALL  HWWARN('UPEVNT',503,*999)
293             IDUP(NUP)=8*IBOS
294          ENDIF
295       ELSEIF (JPR.EQ.17) THEN
296 C---HEAVY QUARKS
297          IF (REMIT) THEN
298 C---3-BODY FINAL STATE
299 C---SET COLOUR CONNECTIONS
300             IF (IC.LE.18) THEN
301                DO I=1,5
302                   CALL UPCODE(ICOL5(I,IC),ICOLUP(1,I))
303                ENDDO
304             ELSE
305                CALL HWWARN('UPEVNT',105,*999)
306             ENDIF
307          ELSE
308 C---2-BODY FINAL STATE
309             IDUP(3)=IDUP(4)
310             IDUP(4)=IDUP(5)
311             CALL HWVEQU(5,PUP(1,4),PUP(1,3))
312             CALL HWVEQU(5,PUP(1,5),PUP(1,4))
313 C---SET COLOUR CONNECTIONS
314             IF (IC.LE.4) THEN
315                DO I=1,4
316                   CALL UPCODE(ICOL4(I,IC),ICOLUP(1,I))
317                ENDDO
318             ELSE
319                CALL HWWARN('UPEVNT',104,*999)
320             ENDIF
321          ENDIF
322       ELSEIF (JPR.EQ.20) THEN
323 C---SINGLE TOP: IA,IB ARE THE QUARKS THAT ARE COLOUR CONNECTED
324 C   I.E. (FOR H EVENTS) THOSE THAT ARE NOT CONNECTED TO GLUON
325          IA=IC/10
326          IB=IC-10*IA
327          IF (IA.LT.1.OR.IA.GT.5) CALL HWWARN('UPEVNT',108,*999)
328          IF (IB.LT.1.OR.IB.GT.5) CALL HWWARN('UPEVNT',109,*999)
329          IF (IA.EQ.IB) CALL HWWARN('UPEVNT',110,*999)
330          DO I=1,5
331             IF (I.EQ.IA.OR.I.EQ.IB) THEN
332                IF (IDUP(I).GT.0) THEN
333                   ICOLUP(1,I)=501
334                   ICOLUP(2,I)=0
335                ELSE
336                   ICOLUP(1,I)=0
337                   ICOLUP(2,I)=501
338                ENDIF
339             ELSEIF (IDUP(I).EQ.21) THEN
340                IG=I
341                ICOLUP(1,I)=502
342                ICOLUP(2,I)=503
343             ELSEIF (IDUP(I).GT.0) THEN
344                ICOLUP(1,I)=502
345                ICOLUP(2,I)=0
346             ELSE
347                ICOLUP(1,I)=0
348                ICOLUP(2,I)=502
349             ENDIF
350          ENDDO
351          IF (REMIT) THEN
352 C---3-BODY FINAL STATE
353 C---COMPLETE GLUON COLOUR CONNECTIONS
354             DO I=1,5
355                IF (I.NE.IA.AND.I.NE.IB.AND.I.NE.IG) THEN
356                   IF (IDUP(I).GT.0) THEN
357                      IF((I.LT.3.AND.IG.LT.3)
358      &              .OR.(I.GT.2.AND.IG.GT.2)) ICOLUP(1,I)=503
359                   ELSE
360                      IF((I.LT.3.AND.IG.GT.2)
361      &              .OR.(I.GT.2.AND.IG.LT.3)) ICOLUP(2,I)=503
362                   ENDIF
363                ENDIF
364             ENDDO
365          ELSE
366 C---2-BODY FINAL STATE
367             IDUP(3)=IDUP(4)
368             IDUP(4)=IDUP(5)
369             ICOLUP(1,3)=ICOLUP(1,4)
370             ICOLUP(2,3)=ICOLUP(2,4)
371             ICOLUP(1,4)=ICOLUP(1,5)
372             ICOLUP(2,4)=ICOLUP(2,5)
373             CALL HWVEQU(5,PUP(1,4),PUP(1,3))
374             CALL HWVEQU(5,PUP(1,5),PUP(1,4))
375          ENDIF
376       ELSE
377 C---BOSON PAIR OR LEPTON PAIR
378          IF (BOPRO.OR.NODEC) THEN
379             NUP=6
380             DO I=6,5,-1
381                CALL HWVEQU(5,PUP(1,I-1),PUP(1,I))
382                IDUP(I)=IDUP(I-1)
383                ISTUP(I)=1
384             ENDDO
385          ELSE
386 C---BOSON PAIR: ONE OR BOTH DECAYED
387 C---ADD BOSON(S) TO EVENT RECORD
388             IF (ABS(IDUP(6)).LT.20) THEN
389                NUP=8
390                I=2
391                IF (ABS(IDUP(4)).LT.20) THEN
392                   NUP=10
393                   I=3
394                ENDIF
395                MUP=NUP-1
396                CALL HWVEQU(5,PUP(1,MUP-I),PUP(1,MUP))
397                CALL HWVEQU(5,PUP(1,NUP-I),PUP(1,NUP))
398                CALL HWVSUM(4,PUP(1,MUP),PUP(1,NUP),PUP(1,6))
399                CALL HWUMAS(PUP(1,6))
400                IDUP(MUP)=IDUP(MUP-I)
401                IDUP(NUP)=IDUP(NUP-I)
402                ISTUP(MUP)=1
403                ISTUP(NUP)=1
404                MOTHUP(1,MUP)=6
405                MOTHUP(2,MUP)=6
406                MOTHUP(1,NUP)=6
407                MOTHUP(2,NUP)=6
408                ISTUP(6)=2
409                ID=IDUP(MUP)+IDUP(NUP)
410                IF (ID.EQ.0) THEN
411                   IDUP(6)=23
412                ELSEIF (ABS(ID).EQ.1) THEN
413                   IDUP(6)=24*ID
414                ELSE
415                   CALL HWWARN('UPEVNT',106,*999)
416                ENDIF
417             ENDIF
418             IF (ABS(IDUP(4)).LT.20) THEN
419                CALL HWVZRO(4,PDB)
420                DO I=8,7,-1
421                   CALL HWVEQU(5,PUP(1,I-3),PUP(1,I))
422                   CALL HWVSUM(4,PUP(1,I),PDB,PDB)
423                   IDUP(I)=IDUP(I-3)
424                   ISTUP(I)=1
425                   MOTHUP(1,I)=5
426                   MOTHUP(2,I)=5
427                ENDDO
428                CALL HWUMAS(PDB)
429                CALL HWVEQU(5,PDB,PUP(1,5))
430                ISTUP(5)=2
431                ID=IDUP(7)+IDUP(8)
432                IF (ID.EQ.0) THEN
433                   IDUP(5)=23
434                ELSEIF (ABS(ID).EQ.1) THEN
435                   IDUP(5)=24*ID
436                ELSE
437                   CALL HWWARN('UPEVNT',107,*999)
438                ENDIF
439             ELSE
440                CALL HWVEQU(5,PUP(1,4),PUP(1,5))
441                IDUP(5)=IDUP(4)
442                ISTUP(5)=1
443                MOTHUP(1,5)=4
444                MOTHUP(2,5)=4
445             ENDIF
446          ENDIF
447 C---ADD DIBOSON OR DILEPTON TO EVENT RECORD (TO FIX ITS MASS)
448          CALL HWVZRO(4,PDB)
449          DO I=6,5,-1
450             CALL HWVSUM(4,PUP(1,I),PDB,PDB)
451             MOTHUP(1,I)=4
452             MOTHUP(2,I)=4
453          ENDDO
454          CALL HWUMAS(PDB)
455          CALL HWVEQU(5,PDB,PUP(1,4))
456          ISTUP(4)=2
457          IDUP(4)=0
458          IF (REMIT) THEN
459 C---SET COLOUR CONNECTIONS
460             DO I=1,3
461                ICOLUP(1,I)=501
462                ICOLUP(2,I)=502
463             ENDDO
464             IF (IHPRO.EQ.1) THEN
465                ICOLUP(2,1)=0
466                ICOLUP(1,2)=0
467             ELSEIF (IHPRO.EQ.2) THEN
468                ICOLUP(1,1)=502
469                ICOLUP(2,1)=0
470                ICOLUP(2,3)=0
471             ELSEIF (IHPRO.EQ.3) THEN
472                ICOLUP(1,1)=0
473                ICOLUP(2,2)=0
474             ELSEIF (IHPRO.EQ.4) THEN
475                ICOLUP(1,1)=0
476                ICOLUP(2,1)=501
477                ICOLUP(1,3)=0
478             ELSEIF (IHPRO.EQ.5) THEN
479                ICOLUP(1,2)=502
480                ICOLUP(2,2)=0
481                ICOLUP(2,3)=0
482             ELSEIF (IHPRO.EQ.6) THEN
483                ICOLUP(1,2)=0
484                ICOLUP(2,2)=501
485                ICOLUP(1,3)=0
486             ELSE
487                CALL HWWARN('UPEVNT',102,*999)
488             ENDIF
489             DO I=4,NUP
490                ICOLUP(1,I)=0
491                ICOLUP(2,I)=0
492             ENDDO
493          ELSE
494             DO I=5,NUP
495                CALL HWVEQU(5,PUP(1,I),PUP(1,I-2))
496                IDUP(I-2)=IDUP(I)
497                ISTUP(I-2)=ISTUP(I)
498                MOTHUP(1,I-2)=MOTHUP(1,I)-2
499                MOTHUP(2,I-2)=MOTHUP(1,I)-2
500             ENDDO
501             MOTHUP(1,3)=1
502             MOTHUP(1,4)=1
503             NUP=NUP-2
504 C---SET COLOUR CONNECTIONS
505             DO I=1,NUP
506                ICOLUP(1,I)=0
507                ICOLUP(2,I)=0
508             ENDDO
509             IF (IDUP(1).GT.0) THEN
510                ICOLUP(1,1)=501
511                ICOLUP(2,2)=501
512             ELSE
513                ICOLUP(2,1)=501
514                ICOLUP(1,2)=501
515             ENDIF
516          ENDIF
517          IF (BOPRO) THEN
518 C---DILEPTON PRODUCTION
519             IBOS=MOD(JJPROC,100)
520             ILEP=MOD(JJPROC,10)
521             IBOS=IBOS-ILEP
522 C---LOAD LEPTON AND BOSON ID
523             I=NUP-1
524             J=NUP
525             IF ( IBOS.EQ.50 .OR.
526      #          (IBOS.EQ.60.AND.JPR.EQ.13) .OR.
527      #          (IBOS.EQ.70.AND.JPR.EQ.13) ) THEN
528                IDUP(I)=-ILEP-10
529                IDUP(J)=-IDUP(I)
530                IF (REMIT) IDUP(4)=23
531             ELSEIF (IBOS.EQ.60.AND.JPR.EQ.14) THEN
532                IDUP(I)=-9-2*ILEP
533                IDUP(J)=1-IDUP(I)
534                IF (REMIT) IDUP(4)=24
535             ELSEIF (IBOS.EQ.70.AND.JPR.EQ.14) THEN
536                IDUP(I)=-10-2*ILEP
537                IDUP(J)=-1-IDUP(I)
538                IF (REMIT) IDUP(4)=-24
539             ELSE
540                CALL HWWARN('UPEVNT',504,*999)
541             ENDIF
542          ENDIF
543       ENDIF
544  999  CONTINUE
545       IF(IERROR.LT.100) RETURN
546       PRINT *
547       DO I=1,NUP
548          PRINT '(4I4,3F8.2)',IDUP(I),ISTUP(I),(ICOLUP(J,I),J=1,2),
549      &        (PUP(J,I),J=1,3)
550       ENDDO
551 c       IPR, IC, NP
552  901  FORMAT(1X,I3,2(1X,I2))
553 c      (ID(I),I=1,NP)
554  902  FORMAT(7(1X,I3))
555 c       XEVWGT
556  903  FORMAT(1X,D14.8)
557 c      ((P(J,I),J=1,4),I=1,NP)
558  904  FORMAT(28(1X,D14.8))
559 c 901  FORMAT(1X,I3,4(1X,I2))
560 c 902  FORMAT(1X,D14.8)
561 c 903  FORMAT(16(1X,D14.8))
562       END
563 C----------------------------------------------------------------------
564       SUBROUTINE UPCODE(ICODE,ICOL)
565 C--DECODES COLOUR CONNECTIONS
566 C----------------------------------------------------------------------
567       IMPLICIT NONE
568       INTEGER ICODE,ICOL(2)
569       ICOL(1)=ICODE/10
570       IF (ICOL(1).NE.0) ICOL(1)=ICOL(1)+500
571       ICOL(2)=MOD(ICODE,10)
572       IF (ICOL(2).NE.0) ICOL(2)=ICOL(2)+500
573       END
574 C----------------------------------------------------------------------
575       SUBROUTINE UPINIT_GUP
576 C----------------------------------------------------------------------
577 C  Reads MC@NLO input headers and fills Les Houches run common HEPRUP
578 C----------------------------------------------------------------------
579       INCLUDE 'herwig65.inc'
580 C--Les Houches Common Blocks
581       INTEGER MAXPUP
582       PARAMETER(MAXPUP=100)
583       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
584       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
585       COMMON /HEPRUP/ IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
586      &                IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),
587      &                XMAXUP(MAXPUP),LPRUP(MAXPUP)
588       INTEGER MAXNUP
589       PARAMETER (MAXNUP=500)
590       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
591       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
592       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
593      &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
594      &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
595      &              SPINUP(MAXNUP)
596       DOUBLE PRECISION XCKECM,XTMP1,XTMP2,XTMP3,XTMP4,XMT,XMW,XMZ,
597      & XMH,XMV,XM1,XM2,XM3,XM4,XM5,XM21,XLAM,GAH,TINY
598       DOUBLE PRECISION XMV1,GAV1,GAMAX1,XMV2,GAV2,GAMAX2
599       INTEGER IVVCODE,IFAIL,MQQ,NQQ,IHW,I,NDNS,JPR,JPR0,IH,
600      & IVHVEC,IVHLEP,IVLEP1,IVLEP2
601       CHARACTER*60 TMPSTR
602       CHARACTER*4 STRP1,STRP2
603       CHARACTER*8 STRGRP
604       CHARACTER*2 STRSCH
605       CHARACTER*3 STRFMT
606       CHARACTER*50 QQIN
607       LOGICAL FK88STRNOEQ
608       DATA TINY/1.D-3/
609       COMMON/NQQCOM/MQQ,NQQ
610       COMMON/VVJIN/QQIN
611       COMMON/VHLIN/IVHVEC,IVHLEP
612       COMMON/VVLIN/IVLEP1,IVLEP2
613 C
614       IF (IERROR.NE.0) RETURN
615 C--SET UP INPUT FILES
616       OPEN(UNIT=61,FILE=QQIN,STATUS='UNKNOWN')
617 C--READ HEADERS OF EVENT FILE
618       READ(61,801)XCKECM,XTMP1,XTMP2,XTMP3,XTMP4,TMPSTR
619       READ(61,802)IVVCODE,TMPSTR
620       IVVCODE=MOD(IVVCODE,10000)
621 C---CHECK PROCESS CODE
622       JPR0=MOD(ABS(IPROC),10000)
623       JPR=JPR0/100
624       IF (JPR.NE.IVVCODE/100) CALL HWWARN('UPINIT',500,*999)
625       IF ((JPR.EQ.17.OR.JPR.EQ.28.OR.JPR.EQ.36).AND.
626      & IVVCODE.NE.MOD(ABS(IPROC),10000)) CALL HWWARN('UPINIT',501,*999)
627       IF (JPR.EQ.13.OR.JPR.EQ.14) THEN
628          IF(JPR0.EQ.1396.OR.JPR0.EQ.1371.OR.
629      #      JPR0.EQ.1372.OR.JPR0.EQ.1373)THEN
630            READ(61,808)EMMIN,EMMAX,TMPSTR
631          ELSE
632            READ(61,809)XMV,GAH,GAMMAX,TMPSTR
633          ENDIF
634 C-- CHECK VECTOR BOSON MASS
635          IF( (IVVCODE.EQ.1397.AND.ABS(XMV-RMASS(200)).GT.TINY) .OR.
636      #       (IVVCODE.EQ.1497.AND.ABS(XMV-RMASS(198)).GT.TINY) .OR.
637      #       (IVVCODE.EQ.1498.AND.ABS(XMV-RMASS(199)).GT.TINY) )
638      #      CALL HWWARN('UPINIT',502,*999)
639       ELSEIF (JPR.EQ.26.OR.JPR.EQ.27) THEN
640          READ(61,810)IVHVEC,IVHLEP,TMPSTR
641          READ(61,809)XMV,GAH,GAMMAX,TMPSTR
642          READ(61,809)XMH,GAH,GAMMAX,TMPSTR
643          IF( (JPR.EQ.26.AND.ABS(XMV-RMASS(199)).GT.TINY) .OR.
644      #       (JPR.EQ.27.AND.ABS(XMV-RMASS(200)).GT.TINY) )
645      #      CALL HWWARN('UPINIT',508,*999)
646          IF(ABS(XMH-RMASS(201)).GT.TINY) CALL HWWARN('UPINIT',509,*999)
647       ELSEIF (JPR.EQ.28) THEN
648          READ(61,808)XMW,XMZ,TMPSTR
649 C-- CHECK VECTOR BOSON MASSES
650          IF(ABS(XMW-RMASS(198)).GT.TINY .OR.
651      #      ABS(XMZ-RMASS(200)).GT.TINY) CALL HWWARN('UPINIT',502,*999)
652          READ(61,810)IVLEP1,IVLEP2,TMPSTR
653          READ(61,809)XMV1,GAV1,GAMAX1,TMPSTR
654          READ(61,809)XMV2,GAV2,GAMAX2,TMPSTR
655       ELSEIF (JPR.EQ.16.OR.JPR.EQ.36) THEN
656          READ(61,809)XMH,GAH,XMT,TMPSTR
657 C-- CHECK HIGGS AND TOP MASSES
658          IH=201
659          IF (JPR.EQ.36) IH=IVVCODE/10-158
660          IF(ABS(XMH-RMASS(IH)).GT.TINY) CALL HWWARN('UPINIT',503,*999)
661          IF(ABS(XMT-RMASS(6)) .GT.TINY) CALL HWWARN('UPINIT',504,*999)
662       ELSEIF (JPR.EQ.17) THEN
663          READ(61,803)XMT,TMPSTR
664 C-- CHECK HEAVY QUARK MASS
665          IF( (IVVCODE.EQ.1706.AND.ABS(XMT-RMASS(6)).GT.TINY) .OR.
666      #       (IVVCODE.EQ.1705.AND.ABS(XMT-RMASS(5)).GT.TINY) .OR.
667      #       (IVVCODE.EQ.1704.AND.ABS(XMT-RMASS(4)).GT.TINY) )
668      #   CALL HWWARN('UPINIT',505,*999)
669       ELSEIF (JPR.EQ.20) THEN
670          READ(61,803)XMT,TMPSTR
671 C-- CHECK HEAVY QUARK MASS
672          IF(ABS(XMT-RMASS(6)).GT.TINY) CALL HWWARN('UPINIT',511,*999)
673       ELSE
674          CALL HWWARN('UPINIT',506,*999)
675       ENDIF
676       READ(61,804)XM1,XM2,XM3,XM4,XM5,XM21,TMPSTR
677       READ(61,805)STRP1,STRP2,TMPSTR
678       READ(61,806)STRGRP,NDNS,TMPSTR
679       READ(61,807)XLAM,STRSCH,TMPSTR
680 C--CHECK THAT EVENT FILE HAS BEEN GENERATED CONSISTENTLY WITH
681 C--HERWIG PARAMETERS ADOPTED HERE
682       IFAIL=0
683 C-- CM ENERGY
684       IF( ABS(XCKECM-PBEAM1-PBEAM2).GT.TINY .OR.
685 C-- QUARK AND GLUON MASSES
686      #     ABS(XM1-RMASS(1)).GT.TINY .OR.
687      #     ABS(XM2-RMASS(2)).GT.TINY .OR.
688      #     ABS(XM3-RMASS(3)).GT.TINY .OR.
689      #     ABS(XM4-RMASS(4)).GT.TINY .OR.
690      #     ABS(XM5-RMASS(5)).GT.TINY .OR.
691      #     ABS(XM21-RMASS(13)).GT.TINY .OR.
692 C-- LAMBDA_QCD: NOW REMOVED TO ALLOW MORE FLEXIBILITY (NNLO EFFECT ANYHOW)
693 C     #     ABS(XLAM-QCDLAM).GT.TINY .OR.
694 C-- REPLACE THE FOLLOWING WITH A CONDITION ON STRSCH, IF CONSISTENT
695 C-- INFORMATION ON PDF SCHEME WILL BE AVAILABLE FROM PDF LIBRARIES AND HERWIG
696 C-- COLLIDING PARTICLE TYPE
697      #     FK88STRNOEQ(STRP1,PART1) .OR.
698      #     FK88STRNOEQ(STRP2,PART2) ) IFAIL=1
699 C--IF PDF LIBRARY IS USED, CHECK PDF CONSISTENCY
700       IF( IFAIL.EQ.0 .AND. MODPDF(1).NE.-1)THEN
701          IF( ABS(NDNS-MODPDF(1)).GT.TINY .OR.
702      #       ABS(NDNS-MODPDF(2)).GT.TINY .OR.
703 C--IF LHAPDF IS USED, STRGRP AND AUTPDF ARE DIFFERENT
704      #       STRGRP.NE.'LHAPDF'.AND.FK88STRNOEQ(STRGRP,AUTPDF(1)))THEN
705                 IFAIL=0
706       ENDIF
707 C--WHEN LHAPDF IS LINKED, AUTPDF() IS A MC@NLO-DEFINED STRING
708          IF(AUTPDF(1).EQ.'HWLHAPDF'.OR.AUTPDF(1).EQ.'LHAEXT')THEN
709             AUTPDF(1)='DEFAULT'
710             AUTPDF(2)='DEFAULT'
711          ENDIF
712       ENDIF
713       IF(IFAIL.EQ.1) CALL HWWARN('UPINIT',507,*999)
714       CALL HWUIDT(3,IDBMUP(1),IHW,PART1)
715       CALL HWUIDT(3,IDBMUP(2),IHW,PART2)
716       DO I=1,2
717          EBMUP(I)=HALF*XCKECM
718          PDFGUP(I)=-1
719          PDFSUP(I)=-1
720       ENDDO
721       IDWTUP=-4
722       NPRUP=1
723       LPRUP(1)=IVVCODE
724 C-- TEST FOR NEW FORMAT INPUT MOMENTA: (PX,PY,PZ,M)
725       READ(61,811) STRFMT,TMPSTR
726       IF (STRFMT.NE.'P,M') CALL HWWARN('UPINIT',510,*999)
727       READ(61,900) MQQ
728       NQQ=0
729 C-- LARGEST EXPECTED NUMBER OF LEGS
730       NUP=10
731       AQEDUP=ZERO
732       AQCDUP=ZERO
733       DO I=1,NUP
734          VTIMUP(I)=ZERO
735          SPINUP(I)=9.
736       ENDDO
737  801  FORMAT(5(1X,D10.4),1X,A)
738  802  FORMAT(1X,I6,1X,A)
739  803  FORMAT(1X,D10.4,1X,A)
740  804  FORMAT(6(1X,D10.4),1X,A)
741  805  FORMAT(2(1X,A4),1X,A)
742  806  FORMAT(1X,A8,1X,I6,1X,A)
743  807  FORMAT(1X,D10.4,1X,A2,1X,A)
744  808  FORMAT(2(1X,D10.4),1X,A)
745  809  FORMAT(3(1X,D10.4),1X,A)
746  810  FORMAT(2(1X,I2),1X,A)
747  811  FORMAT(1X,A3,1X,A)
748  900  FORMAT(I9)
749  999  END