[u/mrichter/AliRoot.git] / HERWIG / src / hwccut.f


CDECK  ID>, HWCCUT.

*CMZ :-        -26/04/91  14.29.39  by  Federico Carminati

*-- Author :    Bryan Webber

C-----------------------------------------------------------------------

      SUBROUTINE HWCCUT(JHEP,KHEP,PCL,BTCLUS,SPLIT)

C-----------------------------------------------------------------------

C     Cuts into 2 the cluster, momentum PCL, made of partons JHEP & KHEP

C-----------------------------------------------------------------------

      INCLUDE 'HERWIG61.INC'

      DOUBLE PRECISION HWREXQ,HWUPCM,HWR,HWVDOT,EMC,QM1,QM2,EMX,EMY,

     & QM3,PXY,PCX,PCY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5),SKAPPA,DELTM,

     & VSCA,VTMP(4),RKAPPA,VCLUS

      INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,ID1,ID2,ID3,NTRY,NTRYMX,J,IB

      LOGICAL BTCLUS,SPLIT

      EXTERNAL HWREXQ,HWUPCM,HWR,HWVDOT,HWRINT

      COMMON/HWCFRM/VCLUS(4,NMXHEP)

      PARAMETER (SKAPPA=1.,NTRYMX=100)

      IF (IERROR.NE.0) RETURN

      EMC=PCL(5)

      ID1=IDHW(JHEP)

      ID2=IDHW(KHEP)

      QM1=RMASS(ID1)

      QM2=RMASS(ID2)

      SPLIT=.FALSE.

      NTRY=0

C Decide if cluster contains a b-(anti)quark

      IF (ID1.EQ.5.OR.ID1.EQ.11.OR.ID2.EQ.5.OR.ID2.EQ.11) THEN

        IB=2

      ELSE

        IB=1

      ENDIF

      IF (BTCLUS) THEN

C Split beam and target clusters as soft clusters

C Both (remnant) children treated like soft clusters if IOPREM=0(1)

  10    ID3=HWRINT(1,2)

        QM3=RMASS(ID3)

        IF (EMC.LE.QM1+QM2+2.*QM3) THEN

          ID3=3-ID3

          QM3=RMASS(ID3)

          IF (EMC.LE.QM1+QM2+2.*QM3) RETURN

        ENDIF

        PXY=EMC-QM1-QM2-TWO*QM3

        IF (ISTHEP(JHEP).EQ.153.OR.ISTHEP(JHEP).EQ.154.OR.

     &      IOPREM.EQ.0) THEN

          EMX=QM1+QM3+HWREXQ(BTCLM,PXY)

        ELSE

          EMX=QM1+QM3+PXY*HWR()**PSPLT(IB)

        ENDIF

        IF (ISTHEP(KHEP).EQ.153.OR.ISTHEP(KHEP).EQ.154.OR.

     &      IOPREM.EQ.0) THEN

          EMY=QM2+QM3+HWREXQ(BTCLM,PXY)

        ELSE

          EMY=QM2+QM3+PXY*HWR()**PSPLT(IB)

        ENDIF

        IF (EMX+EMY.GE.EMC) THEN

          NTRY=NTRY+1

          IF (NTRY.GT.NTRYMX) RETURN

          GOTO 10

        ENDIF

        PCX=HWUPCM(EMX,QM1,QM3)

        PCY=HWUPCM(EMY,QM2,QM3)

      ELSE

C Choose fragment masses for ordinary cluster

        PXY=EMC-QM1-QM2

  20    NTRY=NTRY+1

        IF (NTRY.GT.NTRYMX) RETURN

  30    EMX=QM1+PXY*HWR()**PSPLT(IB)

        EMY=QM2+PXY*HWR()**PSPLT(IB)

        IF (EMX+EMY.GE.EMC) GOTO 30

C u,d,s pair production with weights QWT

  40    ID3=HWRINT(1,3)

        IF (QWT(ID3).LT.HWR()) GOTO 40

        QM3=RMASS(ID3)

        PCX=HWUPCM(EMX,QM1,QM3)

        IF (PCX.LT.ZERO) GOTO 20

        PCY=HWUPCM(EMY,QM2,QM3)

        IF (PCY.LT.ZERO) GOTO 20

        SPLIT=.TRUE.

      ENDIF

C Boost antiquark to CoM frame to find axis

      CALL HWULOF(PCL,PHEP(1,KHEP),AX)

      RCM=1./SQRT(HWVDOT(3,AX,AX))

      CALL HWVSCA(3,RCM,AX,AX)

C Construct new CoM momenta (collinear)

      PXY=HWUPCM(EMC,EMX,EMY)

      CALL HWVSCA(3,PXY,AX,PC)

      PC(4)=SQRT(PXY**2+EMY**2)

      PC(5)=EMY

      CALL HWVSCA(3,PCY,AX,PA)

      PA(4)=SQRT(PCY**2+QM2**2)

      PA(5)=QM2

      CALL HWULOB(PC,PA,PB)

      CALL HWVDIF(4,PC,PB,PA)

      PA(5)=QM3

      LHEP=NHEP+1

      MHEP=NHEP+2

      CALL HWULOB(PCL,PB,PHEP(1,KHEP))

      CALL HWULOB(PCL,PA,PHEP(1,MHEP))

      CALL HWVSCA(3,-ONE,PC,PC)

      PC(4)=EMC-PC(4)

      PC(5)=EMX

      CALL HWVSCA(3,PCX,AX,PA)

      PA(4)=SQRT(PCX**2+QM3**2)

      CALL HWULOB(PC,PA,PB)

      CALL HWULOB(PCL,PB,PHEP(1,LHEP))

      DO 50 J=1,4

  50  PHEP(J,JHEP)=PCL(J)-PHEP(J,KHEP)-PHEP(J,LHEP)-PHEP(J,MHEP)

      PHEP(5,JHEP)=QM1

      CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))

C Construct new vertex positions

      RKAPPA=GEV2MM/SKAPPA

      CALL HWVSCA(3,RKAPPA,AX,AX)

      DELTM=(EMX-EMY)*(EMX+EMY)/(TWO*EMC)

      CALL HWVSCA(3,DELTM,AX,VTMP)

      VTMP(4)=(HALF*EMC-PXY)*RKAPPA

      CALL HWULB4(PCL,VTMP,VTMP)

      CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VHEP(1,LHEP))

      CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))

      VSCA=0.25*EMC+HALF*(PXY+DELTM)

      CALL HWVSCA(3,VSCA,AX,VTMP)

      VTMP(4)=(EMC-VSCA)*RKAPPA

      CALL HWULB4(PCL,VTMP,VTMP)

      CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,MHEP))

      VSCA=-0.25*EMC+HALF*(DELTM-PXY)

      CALL HWVSCA(3,VSCA,AX,VTMP)

      VTMP(4)=(EMC+VSCA)*RKAPPA

      CALL HWULB4(PCL,VTMP,VTMP)

      CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,JHEP))

C (Re-)label quarks

      IDHW(LHEP)=ID3+6

      IDHW(MHEP)=ID3

      IDHEP(MHEP)= IDPDG(ID3)

      IDHEP(LHEP)=-IDPDG(ID3)

      ISTHEP(LHEP)=151

      ISTHEP(MHEP)=151

      JMOHEP(2,JHEP)=LHEP

      JDAHEP(2,KHEP)=MHEP

      JMOHEP(1,LHEP)=JMOHEP(1,KHEP)

      JMOHEP(2,LHEP)=MHEP

      JDAHEP(1,LHEP)=0

      JDAHEP(2,LHEP)=JHEP

      JMOHEP(1,MHEP)=JMOHEP(1,JHEP)

      JMOHEP(2,MHEP)=KHEP

      JDAHEP(1,MHEP)=0

      JDAHEP(2,MHEP)=LHEP

      NHEP=NHEP+2

  999 END

CDECK  ID>, HWCDEC.

*CMZ :-        -26/04/91  10.18.56  by  Bryan Webber

*-- Author :    Bryan Webber

C-----------------------------------------------------------------------

      SUBROUTINE HWCDEC

C-----------------------------------------------------------------------

C     DECAYS CLUSTERS INTO PRIMARY HADRONS

C-----------------------------------------------------------------------

      INCLUDE 'HERWIG61.INC'

      INTEGER JCL,KCL,IP,JP,KP,IST,ID1,ID2,ID3

      IF (IERROR.NE.0) RETURN

      IF (IPROC/1000.EQ.9.OR.IPROC/1000.EQ.5) THEN

C---RELABEL CLUSTER CONNECTED TO REMNANT IN DIS

        DO 10 JCL=2,NHEP

        IF (ISTHEP(JCL).EQ.164) GOTO 20

        IF (ISTHEP(JCL).EQ.165) THEN

          IP=JMOHEP(1,JCL)

          JP=JMOHEP(2,JCL)

          KP=IP

          IF (ISTHEP(IP).EQ.162) THEN

            KP=JP

            JP=IP

          ENDIF

          IF (JMOHEP(2,KP).NE.JP) THEN

            IP=JMOHEP(2,KP)

          ELSE

            IP=JDAHEP(2,KP)

          ENDIF

          KCL=JDAHEP(1,IP)

          IF (ISTHEP(KCL)/10.NE.16) CALL HWWARN('HWCDEC',100,*999)

          ISTHEP(KCL)=164

          GOTO 20

        ENDIF

   10   CONTINUE

      ENDIF

   20 CONTINUE

      DO 30 JCL=1,NHEP

      IST=ISTHEP(JCL)

      IF (IST.GT.162.AND.IST.LT.166) THEN

C---DON'T HADRONIZE BEAM/TARGET CLUSTERS

        IF (IST.EQ.163.OR..NOT.GENSOF) THEN

C---SET UP FLAVOURS FOR CLUSTER DECAY

          CALL HWCFLA(IDHW(JMOHEP(1,JCL)),IDHW(JMOHEP(2,JCL)),ID1,ID3)

          CALL HWCHAD(JCL,ID1,ID3,ID2)

        ENDIF

      ENDIF

   30 CONTINUE

      ISTAT=50

  999 END

CDECK  ID>, HWCFLA.

*CMZ :-        -26/04/91  10.18.56  by  Bryan Webber

*-- Author :    Bryan Webber

C-----------------------------------------------------------------------

      SUBROUTINE HWCFLA(JD1,JD2,ID1,ID2)

C-----------------------------------------------------------------------

C     SETS UP FLAVOURS FOR CLUSTER DECAY

C-----------------------------------------------------------------------

      INTEGER JD1,JD2,ID1,ID2,JD,JDEC(12)

      DATA JDEC/1,2,3,10,11,12,4,5,6,7,8,9/

      JD=JD1

      IF (JD.GT.12) JD=JD-108

      ID1=JDEC(JD)

      JD=JD2

      IF (JD.GT.12) JD=JD-96

      ID2=JDEC(JD-6)

      END
Commit	Line	Data
3820ca8e	1
	2	CDECK ID>, HWCCUT.
	3
	4	*CMZ :- -26/04/91 14.29.39 by Federico Carminati
	5
	6	*-- Author : Bryan Webber
	7
	8	C-----------------------------------------------------------------------
	9
	10	SUBROUTINE HWCCUT(JHEP,KHEP,PCL,BTCLUS,SPLIT)
	11
	12	C-----------------------------------------------------------------------
	13
	14	C Cuts into 2 the cluster, momentum PCL, made of partons JHEP & KHEP
	15
	16	C-----------------------------------------------------------------------
	17
	18	INCLUDE 'HERWIG61.INC'
	19
	20	DOUBLE PRECISION HWREXQ,HWUPCM,HWR,HWVDOT,EMC,QM1,QM2,EMX,EMY,
	21
	22	& QM3,PXY,PCX,PCY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5),SKAPPA,DELTM,
	23
	24	& VSCA,VTMP(4),RKAPPA,VCLUS
	25
	26	INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,ID1,ID2,ID3,NTRY,NTRYMX,J,IB
	27
	28	LOGICAL BTCLUS,SPLIT
	29
	30	EXTERNAL HWREXQ,HWUPCM,HWR,HWVDOT,HWRINT
	31
	32	COMMON/HWCFRM/VCLUS(4,NMXHEP)
	33
	34	PARAMETER (SKAPPA=1.,NTRYMX=100)
	35
	36	IF (IERROR.NE.0) RETURN
	37
	38	EMC=PCL(5)
	39
	40	ID1=IDHW(JHEP)
	41
	42	ID2=IDHW(KHEP)
	43
	44	QM1=RMASS(ID1)
	45
	46	QM2=RMASS(ID2)
	47
	48	SPLIT=.FALSE.
	49
	50	NTRY=0
	51
	52	C Decide if cluster contains a b-(anti)quark
	53
	54	IF (ID1.EQ.5.OR.ID1.EQ.11.OR.ID2.EQ.5.OR.ID2.EQ.11) THEN
	55
	56	IB=2
	57
	58	ELSE
	59
	60	IB=1
	61
	62	ENDIF
	63
	64	IF (BTCLUS) THEN
65
66	C Split beam and target clusters as soft clusters
67
68	C Both (remnant) children treated like soft clusters if IOPREM=0(1)
69
70	10 ID3=HWRINT(1,2)
71
72	QM3=RMASS(ID3)
73
74	IF (EMC.LE.QM1+QM2+2.*QM3) THEN
75
76	ID3=3-ID3
77
78	QM3=RMASS(ID3)
79
80	IF (EMC.LE.QM1+QM2+2.*QM3) RETURN
81
82	ENDIF
83
84	PXY=EMC-QM1-QM2-TWO*QM3
85
86	IF (ISTHEP(JHEP).EQ.153.OR.ISTHEP(JHEP).EQ.154.OR.
87
88	& IOPREM.EQ.0) THEN
89
90	EMX=QM1+QM3+HWREXQ(BTCLM,PXY)
91
92	ELSE
93
94	EMX=QM1+QM3+PXYHWR()*PSPLT(IB)
95
96	ENDIF
97
98	IF (ISTHEP(KHEP).EQ.153.OR.ISTHEP(KHEP).EQ.154.OR.
99
100	& IOPREM.EQ.0) THEN
101
102	EMY=QM2+QM3+HWREXQ(BTCLM,PXY)
103
104	ELSE
105
106	EMY=QM2+QM3+PXYHWR()*PSPLT(IB)
107
108	ENDIF
109
110	IF (EMX+EMY.GE.EMC) THEN
111
112	NTRY=NTRY+1
113
114	IF (NTRY.GT.NTRYMX) RETURN
115
116	GOTO 10
117
118	ENDIF
119
120	PCX=HWUPCM(EMX,QM1,QM3)
121
122	PCY=HWUPCM(EMY,QM2,QM3)
123
124	ELSE
125
126	C Choose fragment masses for ordinary cluster
127
128	PXY=EMC-QM1-QM2
129
130	20 NTRY=NTRY+1
131
132	IF (NTRY.GT.NTRYMX) RETURN
133
134	30 EMX=QM1+PXYHWR()*PSPLT(IB)
135
136	EMY=QM2+PXYHWR()*PSPLT(IB)
137
138	IF (EMX+EMY.GE.EMC) GOTO 30
139
140	C u,d,s pair production with weights QWT
141
142	40 ID3=HWRINT(1,3)
143
144	IF (QWT(ID3).LT.HWR()) GOTO 40
145
146	QM3=RMASS(ID3)
147
148	PCX=HWUPCM(EMX,QM1,QM3)
149
150	IF (PCX.LT.ZERO) GOTO 20
151
152	PCY=HWUPCM(EMY,QM2,QM3)
153
154	IF (PCY.LT.ZERO) GOTO 20
155
156	SPLIT=.TRUE.
157
158	ENDIF
159
160	C Boost antiquark to CoM frame to find axis
161
162	CALL HWULOF(PCL,PHEP(1,KHEP),AX)
163
164	RCM=1./SQRT(HWVDOT(3,AX,AX))
165
166	CALL HWVSCA(3,RCM,AX,AX)
167
168	C Construct new CoM momenta (collinear)
169
170	PXY=HWUPCM(EMC,EMX,EMY)
171
172	CALL HWVSCA(3,PXY,AX,PC)
173
174	PC(4)=SQRT(PXY2+EMY2)
175
176	PC(5)=EMY
177
178	CALL HWVSCA(3,PCY,AX,PA)
179
180	PA(4)=SQRT(PCY2+QM22)
181
182	PA(5)=QM2
183
184	CALL HWULOB(PC,PA,PB)
185
186	CALL HWVDIF(4,PC,PB,PA)
187
188	PA(5)=QM3
189
190	LHEP=NHEP+1
191
192	MHEP=NHEP+2
193
194	CALL HWULOB(PCL,PB,PHEP(1,KHEP))
195
196	CALL HWULOB(PCL,PA,PHEP(1,MHEP))
197
198	CALL HWVSCA(3,-ONE,PC,PC)
199
200	PC(4)=EMC-PC(4)
201
202	PC(5)=EMX
203
204	CALL HWVSCA(3,PCX,AX,PA)
205
206	PA(4)=SQRT(PCX2+QM32)
207
208	CALL HWULOB(PC,PA,PB)
209
210	CALL HWULOB(PCL,PB,PHEP(1,LHEP))
211
212	DO 50 J=1,4
213
214	50 PHEP(J,JHEP)=PCL(J)-PHEP(J,KHEP)-PHEP(J,LHEP)-PHEP(J,MHEP)
215
216	PHEP(5,JHEP)=QM1
217
218	CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))
219
220	C Construct new vertex positions
221
222	RKAPPA=GEV2MM/SKAPPA
223
224	CALL HWVSCA(3,RKAPPA,AX,AX)
225
226	DELTM=(EMX-EMY)(EMX+EMY)/(TWOEMC)
227
228	CALL HWVSCA(3,DELTM,AX,VTMP)
229
230	VTMP(4)=(HALFEMC-PXY)RKAPPA
231
232	CALL HWULB4(PCL,VTMP,VTMP)
233
234	CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VHEP(1,LHEP))
235
236	CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))
237
238	VSCA=0.25EMC+HALF(PXY+DELTM)
239
240	CALL HWVSCA(3,VSCA,AX,VTMP)
241
242	VTMP(4)=(EMC-VSCA)*RKAPPA
243
244	CALL HWULB4(PCL,VTMP,VTMP)
245
246	CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,MHEP))
247
248	VSCA=-0.25EMC+HALF(DELTM-PXY)
249
250	CALL HWVSCA(3,VSCA,AX,VTMP)
251
252	VTMP(4)=(EMC+VSCA)*RKAPPA
253
254	CALL HWULB4(PCL,VTMP,VTMP)
255
256	CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,JHEP))
257
258	C (Re-)label quarks
259
260	IDHW(LHEP)=ID3+6
261
262	IDHW(MHEP)=ID3
263
264	IDHEP(MHEP)= IDPDG(ID3)
265
266	IDHEP(LHEP)=-IDPDG(ID3)
267
268	ISTHEP(LHEP)=151
269
270	ISTHEP(MHEP)=151
271
272	JMOHEP(2,JHEP)=LHEP
273
274	JDAHEP(2,KHEP)=MHEP
275
276	JMOHEP(1,LHEP)=JMOHEP(1,KHEP)
277
278	JMOHEP(2,LHEP)=MHEP
279
280	JDAHEP(1,LHEP)=0
281
282	JDAHEP(2,LHEP)=JHEP
283
284	JMOHEP(1,MHEP)=JMOHEP(1,JHEP)
285
286	JMOHEP(2,MHEP)=KHEP
287
288	JDAHEP(1,MHEP)=0
289
290	JDAHEP(2,MHEP)=LHEP
291
292	NHEP=NHEP+2
293
294	999 END
295
296	CDECK ID>, HWCDEC.
297
298	*CMZ :- -26/04/91 10.18.56 by Bryan Webber
299
300	*-- Author : Bryan Webber
301
302	C-----------------------------------------------------------------------
303
304	SUBROUTINE HWCDEC
305
306	C-----------------------------------------------------------------------
307
308	C DECAYS CLUSTERS INTO PRIMARY HADRONS
309
310	C-----------------------------------------------------------------------
311
312	INCLUDE 'HERWIG61.INC'
313
314	INTEGER JCL,KCL,IP,JP,KP,IST,ID1,ID2,ID3
315
316	IF (IERROR.NE.0) RETURN
317
318	IF (IPROC/1000.EQ.9.OR.IPROC/1000.EQ.5) THEN
319
320	C---RELABEL CLUSTER CONNECTED TO REMNANT IN DIS
321
322	DO 10 JCL=2,NHEP
323
324	IF (ISTHEP(JCL).EQ.164) GOTO 20
325
326	IF (ISTHEP(JCL).EQ.165) THEN
327
328	IP=JMOHEP(1,JCL)
329
330	JP=JMOHEP(2,JCL)
331
332	KP=IP
333
334	IF (ISTHEP(IP).EQ.162) THEN
335
336	KP=JP
337
338	JP=IP
339
340	ENDIF
341
342	IF (JMOHEP(2,KP).NE.JP) THEN
343
344	IP=JMOHEP(2,KP)
345
346	ELSE
347
348	IP=JDAHEP(2,KP)
349
350	ENDIF
351
352	KCL=JDAHEP(1,IP)
353
354	IF (ISTHEP(KCL)/10.NE.16) CALL HWWARN('HWCDEC',100,*999)
355
356	ISTHEP(KCL)=164
357
358	GOTO 20
359
360	ENDIF
361
362	10 CONTINUE
363
364	ENDIF
365
366	20 CONTINUE
367
368	DO 30 JCL=1,NHEP
369
370	IST=ISTHEP(JCL)
371
372	IF (IST.GT.162.AND.IST.LT.166) THEN
373
374	C---DON'T HADRONIZE BEAM/TARGET CLUSTERS
375
376	IF (IST.EQ.163.OR..NOT.GENSOF) THEN
377
378	C---SET UP FLAVOURS FOR CLUSTER DECAY
379
380	CALL HWCFLA(IDHW(JMOHEP(1,JCL)),IDHW(JMOHEP(2,JCL)),ID1,ID3)
381
382	CALL HWCHAD(JCL,ID1,ID3,ID2)
383
384	ENDIF
385
386	ENDIF
387
388	30 CONTINUE
389
390	ISTAT=50
391
392	999 END
393
394	CDECK ID>, HWCFLA.
395
396	*CMZ :- -26/04/91 10.18.56 by Bryan Webber
397
398	*-- Author : Bryan Webber
399
400	C-----------------------------------------------------------------------
401
402	SUBROUTINE HWCFLA(JD1,JD2,ID1,ID2)
403
404	C-----------------------------------------------------------------------
405
406	C SETS UP FLAVOURS FOR CLUSTER DECAY
407
408	C-----------------------------------------------------------------------
409
410	INTEGER JD1,JD2,ID1,ID2,JD,JDEC(12)
411
412	DATA JDEC/1,2,3,10,11,12,4,5,6,7,8,9/
413
414	JD=JD1
415
416	IF (JD.GT.12) JD=JD-108
417
418	ID1=JDEC(JD)
419
420	JD=JD2
421
422	IF (JD.GT.12) JD=JD-96
423
424	ID2=JDEC(JD-6)
425
426	END