]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HERWIG/src/hwegam.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
splitting of simulation and reconstruction code (T.Kuhr)
[u/mrichter/AliRoot.git] / HERWIG / src / hwegam.f
CommitLineData
3820ca8e 1
2CDECK ID>, HWEGAM.
3
4*CMZ :- -26/04/91 11.11.55 by Bryan Webber
5
6*-- Author : Bryan Webber & Luca Stanco
7
8C-----------------------------------------------------------------------
9
10 SUBROUTINE HWEGAM(IHEP,ZMI,ZMA,WWA)
11
12C-----------------------------------------------------------------------
13
14C GENERATES A PHOTON IN WEIZSACKER-WILLIAMS (WWA=.TRUE.) OR
15
16C ELSE EQUIVALENT PHOTON APPROX FROM INCOMING E+, E-, MU+ OR MU-
17
18C-----------------------------------------------------------------------
19
20 INCLUDE 'HERWIG61.INC'
21
22 DOUBLE PRECISION HWR,HWRUNI,EGMIN,ZMIN,ZMAX,ZGAM,SS,ZMI,ZMA,
23
24 & PPL,PMI,QT2,Q2,QQMIN,QQMAX,S0,A,RPM(2)
25
26 INTEGER IHEP,IHADIS,HQ,I
27
28 LOGICAL WWA
29
30 EXTERNAL HWR,HWRUNI
31
32 DATA EGMIN/5.D0/
33
34 IF (IERROR.NE.0) RETURN
35
36 IF (IHEP.LT.1.OR.IHEP.GT.2) CALL HWWARN('HWEGAM',500,*999)
37
38 SS=PHEP(5,3)
39
40 IF (IHEP.EQ.1) THEN
41
42 IHADIS=2
43
44 ELSE
45
46 IHADIS=1
47
48 IF (JDAHEP(1,IHADIS).NE.0) IHADIS=JDAHEP(1,IHADIS)
49
50 ENDIF
51
52C---DEFINE LIMITS FOR GAMMA MOMENTUM FRACTION
53
54 IF (ZMI.LE.ZERO .OR. ZMA.GT.ONE) THEN
55
56 IF (IPRO.EQ.13.OR.IPRO.EQ.14) THEN
57
58 S0 = EMMIN**2
59
60 ELSEIF(IPRO.EQ.15.OR.IPRO.EQ.18.OR.IPRO.EQ.22.OR.IPRO.EQ.24.OR.
61
62 & IPRO.EQ.50.OR.IPRO.EQ.53.OR.IPRO.EQ.55)THEN
63
64 S0 = 4.D0*PTMIN**2
65
66 ELSEIF (IPRO.EQ.17.OR.IPRO.EQ.51) THEN
67
68 HQ = MOD(IPROC,100)
69
70 S0 = 4.D0*(PTMIN**2+RMASS(HQ)**2)
71
72 ELSEIF (IPRO.EQ.16.OR.IPRO.EQ.19.OR.IPRO.EQ.95) THEN
73
74 S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
75
76 ELSEIF (IPRO.EQ.23) THEN
77
78 S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
79
80 S0 = (PTMIN+SQRT(PTMIN**2+S0))**2
81
82 ELSEIF (IPRO.EQ.20) THEN
83
84 S0 = RMASS(201)**2
85
86 ELSEIF (IPRO.EQ.21) THEN
87
88 S0 = (PTMIN+SQRT(PTMIN**2+RMASS(198)**2))**2
89
90C--PR MOD 7/7/99
91
92 ELSEIF(IPRO.EQ.30) THEN
93
94 S0 = 4.0D0*(PTMIN**2+RMMNSS**2)
95
96 ELSEIF(IPRO.EQ.40.OR.IPRO.EQ.41) THEN
97
98 HQ = IPROC-100*IPRO
99
100 RPM(1) = RMMNSS
101
102 RPM(2) = ZERO
103
104 IF(HQ.GE.10.AND.HQ.LT.20) THEN
105
106 RPM(1) = ABS(RMASS(450))
107
108 IF(HQ.GT.10) RPM(1) = ABS(RMASS(449+MOD(HQ,10)))
109
110 ELSEIF(HQ.GE.20.AND.HQ.LT.30) THEN
111
112 RPM(1) = ABS(RMASS(454))
113
114 IF(HQ.GT.20) RPM(1) = ABS(RMASS(453+MOD(HQ,20)))
115
116 ELSEIF(HQ.EQ.30) THEN
117
118 RPM(1) = RMASS(449)
119
120 ELSEIF(HQ.EQ.40) THEN
121
122 IF(IPRO.EQ.40) THEN
123
124 RPM(1) = RMASS(425)
125
126 DO I=1,5
127
128 RPM(1) = MIN(RPM(1),RMASS(425+I))
129
130 ENDDO
131
132 ELSE
133
134 RPM(1) = MIN(RMASS(405),RMASS(406))
135
136 ENDIF
137
138 RPM(2) = RMASS(198)
139
140 ELSEIF(HQ.EQ.50) THEN
141
142 IF(IPRO.EQ.40) THEN
143
144 RPM(1) = RMASS(425)
145
146 DO I=1,5
147
148 RPM(1) = MIN(RPM(1),RMASS(425+I))
149
150 ENDDO
151
152 DO I=1,3
153
154 RPM(2) = MIN(RPM(1),RMASS(433+2*I))
155
156 ENDDO
157
158 RPM(1) = MIN(RPM(1),RPM(2))
159
160 RPM(2) = RMASS(203)
161
162 DO I=1,2
163
164 RPM(2) = MIN(RPM(2),RMASS(204+I))
165
166 ENDDO
167
168 ELSE
169
170 RPM(1) = RMASS(401)
171
172 RPM(2) = RMASS(413)
173
174 DO I=1,5
175
176 RPM(1) = MIN(RPM(1),RMASS(401+I))
177
178 RPM(2) = MIN(RPM(2),RMASS(413+I))
179
180 ENDDO
181
182 RPM(1) = MIN(RPM(1),RPM(2))
183
184 RPM(2) = RMASS(203)
185
186 DO I=1,2
187
188 RPM(2) = MIN(RPM(2),RMASS(204+I))
189
190 ENDDO
191
192 ENDIF
193
194 RPM(2) = RMASS(203)
195
196 DO I=1,2
197
198 RPM(2) = MIN(RPM(2),RMASS(204+I))
199
200 ENDDO
201
202 ELSEIF(HQ.GE.60) THEN
203
204 RPM(1) = ZERO
205
206 ENDIF
207
208 RPM(1) = RPM(1)**2
209
210 RPM(2) = RPM(2)**2
211
212 S0 = RPM(1)+RPM(2)+TWO*(PTMIN**2+
213
214 & SQRT(RPM(1)*RPM(2)+PTMIN**2*(RPM(1)+RPM(2)+PTMIN**2)))
215
216C--end of mod
217
218 ELSEIF (IPRO.EQ.52) THEN
219
220 HQ = MOD(IPROC,100)
221
222 S0 = (PTMIN+SQRT(PTMIN**2+RMASS(HQ)**2))**2
223
224 ELSEIF (IPRO.EQ.60) THEN
225
226 HQ = MOD(IPROC,100)
227
228 IF (HQ.EQ.0) THEN
229
230 S0 = 4.D0*PTMIN**2
231
232 ELSE
233
234 IF (HQ.GT.6) HQ=2*HQ+107
235
236 IF (HQ.EQ.127) HQ=198
237
238 S0 = 4.D0*(PTMIN**2+RMASS(HQ)**2)
239
240 ENDIF
241
242 ELSEIF (IPRO.EQ.80) THEN
243
244 S0 = WHMIN**2
245
246 ELSEIF (IPRO.EQ.90) THEN
247
248 S0 = Q2MIN
249
250 ELSEIF (IPRO.EQ.91.OR.IPRO.EQ.92) THEN
251
252 S0 = Q2MIN+4.D0*PTMIN**2
253
254 HQ = MOD(IPROC,100)
255
256 IF (HQ.GT.0) S0 = S0+4.D0*RMASS(HQ)**2
257
258 IF (IPRO.EQ.91) S0 = MAX(S0,EMMIN**2)
259
260 ELSE
261
262 S0 = 0
263
264 ENDIF
265
266 IF (S0.GT.ZERO) THEN
267
268 S0 = (SQRT(S0)+ABS(PHEP(5,IHADIS)))**2-PHEP(5,IHADIS)**2
269
270 S0 = MAX(S0,WHMIN**2)
271
272 ZMIN = S0 / (SS**2 - PHEP(5,IHEP)**2 - PHEP(5,IHADIS)**2)
273
274 ZMAX = ONE
275
276 ELSE
277
278C---UNKNOWN PROCESS: USE ENERGY CUTOFF, AND WARN USER
279
280 IF (FSTWGT) CALL HWWARN('HWEGAM',1,*999)
281
282 ZMIN = EGMIN / PHEP(4,IHEP)
283
284 ZMAX = ONE
285
286 ENDIF
287
288 ELSE
289
290 ZMIN=ZMI
291
292 ZMAX=ZMA
293
294 ENDIF
295
296C---APPLY USER DEFINED CUTS YWWMIN,YWWMAX AND INDIRECT LIMITS ON Z
297
298 IF (.NOT.WWA) THEN
299
300 ZMIN=MAX(ZMIN,YWWMIN,SQRT(Q2WWMN)/ABS(PHEP(3,IHEP)))
301
302 ZMAX=MIN(ZMAX,YWWMAX)
303
304 IF (ZMIN.GT.ZMAX) THEN
305
306 GAMWT=ZERO
307
308 RETURN
309
310 ENDIF
311
312 ENDIF
313
314C---GENERATE GAMMA MOMENTUM FRACTION
315
316 A=HALF
317
318 10 IF (HWR().LT.A) THEN
319
320 ZGAM=(ZMIN/ZMAX)**HWR()*ZMAX
321
322 ELSE
323
324 ZGAM=(ZMAX-ZMIN)*HWR()+ZMIN
325
326 ENDIF
327
328 GAMWT = GAMWT * .5*ALPHEM/PIFAC *
329
330 + (1+(1-ZGAM)**2)/(A/LOG(ZMAX/ZMIN)+(1-A)/(ZMAX-ZMIN)*ZGAM)
331
332 IF (WWA) THEN
333
334 GAMWT = GAMWT * LOG((ONE-ZGAM)/ZGAM*(SS/PHEP(5,IHEP))**2)
335
336 ELSE
337
338C---Q2WWMN AND Q2WWMX ARE USER-DEFINED LIMITS IN THE Q**2 INTEGRATION
339
340 QQMAX=MIN(Q2WWMX,(ZGAM*PHEP(3,IHEP))**2)
341
342 QQMIN=MAX(Q2WWMN,(PHEP(5,IHEP)*ZGAM)**2/(1.-ZGAM))
343
344 IF (QQMIN.GT.QQMAX) CALL HWWARN('HWEGAM',50,*10)
345
346 Q2=EXP(HWRUNI(0,LOG(QQMIN),LOG(QQMAX)))
347
348 GAMWT = GAMWT * LOG(QQMAX/QQMIN)
349
350 ENDIF
351
352 IF (GAMWT.LT.ZERO) GAMWT=ZERO
353
354C---FILL PHOTON
355
356 NHEP=NHEP+1
357
358 IDHW(NHEP)=59
359
360 ISTHEP(NHEP)=3
361
362 IDHEP(NHEP)=22
363
364 JMOHEP(1,NHEP)=IHEP
365
366 JMOHEP(2,NHEP)=0
367
368 JDAHEP(1,NHEP)=0
369
370 JDAHEP(2,NHEP)=0
371
372 JDAHEP(1,IHEP)=NHEP
373
374 IF (WWA) THEN
375
376C---FOR COLLINEAR KINEMATICS, ZGAM IS THE ENERGY FRACTION
377
378 PHEP(4,NHEP)=PHEP(4,IHEP)*ZGAM
379
380 PHEP(3,NHEP)=PHEP(3,IHEP)-SIGN(SQRT(
381
382 & (PHEP(4,IHEP)-PHEP(4,NHEP))**2-PHEP(5,IHEP)**2),PHEP(3,IHEP))
383
384 PHEP(2,NHEP)=0
385
386 PHEP(1,NHEP)=0
387
388 CALL HWUMAS(PHEP(1,NHEP))
389
390 ELSE
391
392C---FOR EXACT KINEMATICS, ZGAM IS TAKEN TO BE FRACTION OF (E+PZ)
393
394 PPL=ZGAM*(ABS(PHEP(3,IHEP))+PHEP(4,IHEP))
395
396 QT2=(ONE-ZGAM)*Q2-(ZGAM*PHEP(5,IHEP))**2
397
398 PMI=(QT2-Q2)/PPL
399
400 PHEP(5,NHEP)=-SQRT(Q2)
401
402 PHEP(4,NHEP)=(PPL+PMI)/TWO
403
404 PHEP(3,NHEP)=SIGN((PPL-PMI)/TWO,PHEP(3,IHEP))
405
406 CALL HWRAZM(SQRT(QT2),PHEP(1,NHEP),PHEP(2,NHEP))
407
408 ENDIF
409
410C---UPDATE OVERALL CM FRAME
411
412 JMOHEP(IHEP,3)=NHEP
413
414 CALL HWVDIF(4,PHEP(1,3),PHEP(1,IHEP),PHEP(1,3))
415
416 CALL HWVSUM(4,PHEP(1,NHEP),PHEP(1,3),PHEP(1,3))
417
418 CALL HWUMAS(PHEP(1,3))
419
420C---FILL OUTGOING LEPTON
421
422 NHEP=NHEP+1
423
424 IDHW(NHEP)=IDHW(IHEP)
425
426 ISTHEP(NHEP)=1
427
428 IDHEP(NHEP)=IDHEP(IHEP)
429
430 JMOHEP(1,NHEP)=IHEP
431
432 JMOHEP(2,NHEP)=0
433
434 JDAHEP(1,NHEP)=0
435
436 JDAHEP(2,NHEP)=0
437
438 JDAHEP(2,IHEP)=NHEP
439
440 CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,NHEP-1),PHEP(1,NHEP))
441
442 PHEP(5,NHEP)=PHEP(5,IHEP)
443
444 999 END
445
446CDECK ID>, HWEINI.
447
448*CMZ :- -26/04/91 12.42.30 by Federico Carminati
449
450*-- Author : Bryan Webber
451
452C-----------------------------------------------------------------------
453
454 SUBROUTINE HWEINI
455
456C-----------------------------------------------------------------------
457
458C INITIALISES ELEMENTARY PROCESS
459
460C-----------------------------------------------------------------------
461
462 INCLUDE 'HERWIG61.INC'
463
464 DOUBLE PRECISION HWRSET,DUMMY,SAFETY
465
466 EXTERNAL HWRSET
467
468 PARAMETER (SAFETY=1.001)
469
470 INTEGER NBSH,I
471
472C---NO OF WEIGHT GENERATED
473
474 NWGTS=0
475
476C---ACCUMULATED WEIGHTS
477
478 WGTSUM=0.
479
480C---ACCUMULATED WEIGHT-SQUARED
481
482 WSQSUM=0.
483
484C---CURRENT MAX WEIGHT
485
486 WBIGST=0.
487
488C---LAST VALUE OF SCALE
489
490 EMLST=0.
491
492C---NUMBER OF ERRORS REPORTED
493
494 NUMER=0
495
496C---NUMBER OF ERRORS UNREPORTED
497
498 NUMERU=0
499
500C---FIND MAXIMUM EVENT WEIGHT IN CASES WHERE THIS IS REQUIRED
501
502 IF (NOWGT) THEN
503
504 IF (WGTMAX.EQ.ZERO) THEN
505
506 NBSH=IBSH
507
508 DUMMY = HWRSET(IBRN)
509
510 WRITE(6,10) IPROC,IBRN,NBSH
511
512 10 FORMAT(/10X,'INITIAL SEARCH FOR MAX WEIGHT'//
513
514 & 10X,'PROCESS CODE IPROC = ',I11/
515
516 & 10X,'RANDOM NO. SEED 1 = ',I11/
517
518 & 10X,' SEED 2 = ',I11/
519
520 & 10X,'NUMBER OF SHOTS = ',I11)
521
522 NEVHEP=0
523
524 DO 11 I=1,NBSH
525
526 CALL HWEPRO
527
528 11 CONTINUE
529
530 WRITE(6,20)
531
532 20 FORMAT(/10X,'INITIAL SEARCH FINISHED')
533
534 IF (WBIGST*NWGTS.LT.SAFETY*WGTSUM)
535
536 & WGTMAX=SAFETY*WBIGST
537
538 CALL HWEFIN
539
540 NWGTS=0
541
542 WGTSUM=0.
543
544 WSQSUM=0.
545
546 WBIGST=0.
547
548 ELSE
549
550 WRITE(6,21) AVWGT,WGTMAX
551
552 21 FORMAT(/1P,10X,'INPUT EVT WEIGHT =',E12.4/
553
554 & 10X,'INPUT MAX WEIGHT =',E12.4)
555
556 ENDIF
557
558 ENDIF
559
560C---RESET RANDOM NUMBER
561
562 DUMMY = HWRSET(NRN)
563
564 ISTAT=5
565
566 999 END
567
568CDECK ID>, HWEISR.
569
570*CMZ :- -01/04/99 19.55.17 by Mike Seymour
571
572*-- Author : Mike Seymour
573
574C-----------------------------------------------------------------------
575
576 SUBROUTINE HWEISR(IHEP)
577
578C-----------------------------------------------------------------------
579
580C GENERATES AN ISR PHOTON FROM INCOMING E+, E-, MU+ OR MU-
581
582C-----------------------------------------------------------------------
583
584 INCLUDE 'HERWIG61.INC'
585
586 DOUBLE PRECISION HWR,QSQMAX,QSQMIN,A,B,B1,B2,B3,B4,B5,B6,B7,B8,
587
588 $ R,AA,T0,T1,C1,C2,T,Z(2),QSQ(2),PHI(2),C
589
590 INTEGER IHEP,I,J
591
592 EXTERNAL HWR
593
594 SAVE Z,QSQ,PHI
595
596C---IF ZMXISR IS ZERO, THERE CAN BE NO ISR
597
598 IF (ZMXISR.EQ.ZERO .OR. (IPRO.GT.3.AND.IPRO.NE.6)) RETURN
599
600C---CHECK CONSISTENCY OF TMNISR AND ZMXISR
601
602 IF (ZMXISR**2.LT.TMNISR) CALL HWWARN('HWEISR',200,*999)
603
604C---CALCULATE VIRTUALITY LIMITS
605
606 QSQMAX=4*PHEP(4,IHEP)**2
607
608 QSQMIN=PHEP(5,IHEP)**2
609
610C---AND THEREFORE THE Z DEPENDENCE
611
612 A=ALPHEM/PIFAC
613
614 B=A*(LOG(QSQMAX/QSQMIN)-1)
615
616C---DECIDE HOW MUCH WEIGHT TO GIVE THE Z RESONANCE
617
618 IF (IHEP.EQ.1) THEN
619
620 IF (IPRO.EQ.1.OR.IPRO.EQ.6) THEN
621
622 AA=10
623
624 ELSEIF (IPRO.EQ.2) THEN
625
626 AA=0
627
628 ELSEIF (IPRO.EQ.3) THEN
629
630 AA=1
631
632 ELSE
633
634 RETURN
635
636 ENDIF
637
638 T0=RMASS(200)**2/QSQMAX
639
640 T1=GAMZ*RMASS(200)/QSQMAX
641
642 IF (T0.GT.ONE) THEN
643
644 T0=0
645
646 AA=0
647
648 ENDIF
649
650 AA=AA*(1-T0)
651
652C---GENERATE A T VALUE BETWEEN TMNISR AND 1 ACCORDING TO:
653
654C ( b**2*log(zmxisr**2/t)/t + 2*b*(1-(1-zmxisr)**b)*((1-t)**(2*b-1)+1/t
655
656C +(1-t0)**(2b-1)*aa*t1/((t-t0)**2+t1**2)) ) *theta(zmxisr**2-t)
657
658C +( 2*b*(1-zmxisr)**b*((1-t)**(b-1)+1/t
659
660C +(1-t0)**(b-1)*aa*t1/((t-t0)**2+t1**2)) ) *theta(zmxisr-t)
661
662C +( (1-zmxisr)**(2*b) ) *delta(1-t)
663
664 B1=(1-ZMXISR)**(2*B)
665
666 B2=B1+2*(1-ZMXISR)**B*((1-TMNISR)**B-(1-ZMXISR)**B)
667
668 B3=B2+2*B*(1-ZMXISR)**B*LOG(ZMXISR/TMNISR)
669
670 B4=B3+2*B*(1-ZMXISR)**B*AA*(1-T0)**(B-1)
671
672 $ *(ATAN((ZMXISR-T0)/T1)-ATAN((TMNISR-T0)/T1))
673
674 B5=B4+(1-(1-ZMXISR)**B)*((1-TMNISR)**(2*B)-(1-ZMXISR**2)**(2*B))
675
676 B6=B5+2*B*(1-(1-ZMXISR)**B)*LOG(ZMXISR**2/TMNISR)
677
678 B7=B6+B**2*LOG(ZMXISR**2/TMNISR)**2/2
679
680 B8=B7+2*B*(1-(1-ZMXISR)**B)*AA*(1-T0)**(2*B-1)
681
682 $ *(ATAN((ZMXISR**2-T0)/T1)-ATAN((TMNISR-T0)/T1))
683
684 R=B8*HWR()
685
686 IF (R.LE.B1) THEN
687
688C---NEITHER EMITS
689
690 T=1
691
692 GAMWT=GAMWT*B8/B1
693
694 Z(1)=1
695
696 ELSEIF (R.LE.B4) THEN
697
698C---ONE EMITS
699
700 IF (R.LE.B2) THEN
701
702 R=(R-B1)/(B2-B1)
703
704 T=1-(1-TMNISR)*(1-R*(1-((1-ZMXISR)/(1-TMNISR))**B))**(1/B)
705
706 ELSEIF (R.LE.B3) THEN
707
708 R=(R-B2)/(B3-B2)
709
710 T=(TMNISR/ZMXISR)**R*ZMXISR
711
712 ELSE
713
714 R=(R-B3)/(B4-B3)
715
716 T=T0+T1*TAN(
717
718 $ ATAN((ZMXISR-T0)/T1)*R+ATAN((TMNISR-T0)/T1)*(1-R))
719
720 ENDIF
721
722 GAMWT=GAMWT*B8/(2*B*(1-ZMXISR)**B*((1-T)**(B-1)+1/T+
723
724 $ (1-T0)**(B-1)*AA*T1/((T-T0)**2+T1**2)))
725
726 Z(1)=1
727
728 IF (HWR().GT.HALF) Z(1)=T
729
730 GAMWT=GAMWT*2
731
732 ELSE
733
734C---BOTH EMIT
735
736 IF (R.LE.B5) THEN
737
738 R=(R-B4)/(B5-B4)
739
740 T=1-(1-TMNISR)*
741
742 $ (1-R*(1-((1-ZMXISR**2)/(1-TMNISR))**(2*B)))**(.5/B)
743
744 ELSEIF (R.LE.B6) THEN
745
746 R=(R-B5)/(B6-B5)
747
748 T=(TMNISR/ZMXISR**2)**R*ZMXISR**2
749
750 ELSEIF (R.LE.B7) THEN
751
752 R=(R-B6)/(B7-B6)
753
754 T=(TMNISR/ZMXISR**2)**SQRT(R)*ZMXISR**2
755
756 ELSE
757
758 R=(R-B7)/(B8-B7)
759
760 T=T0+T1*TAN(
761
762 $ ATAN((ZMXISR**2-T0)/T1)*R+ATAN((TMNISR-T0)/T1)*(1-R))
763
764 ENDIF
765
766 GAMWT=GAMWT*B8/(B**2*LOG(ZMXISR**2/T)/T
767
768 $ + 2*B*(1-(1-ZMXISR)**B)*((1-T)**(2*B-1)+1/T+
769
770 $ (1-T0)**(B-1)*AA*T1/((T-T0)**2+T1**2)))
771
772C---GENERATE A Z VALUE BETWEEN T/ZMXISR AND ZMXISR ACCORDING TO:
773
774C 1/z+(1-z)**(b-1)+t/z**2*(1-t/z)**(b-1)
775
776 C1=LOG(ZMXISR**2/T)
777
778 C2=C1+2/B*((1-T/ZMXISR)**B-(1-ZMXISR)**B)
779
780 IF (C2.GT.ZERO) THEN
781
782 R=C2*HWR()
783
784 IF (R.LE.C1) THEN
785
786 Z(1)=(T/ZMXISR**2)**HWR()*ZMXISR
787
788 ELSE
789
790 Z(1)=1-(1-T/ZMXISR)*
791
792 $ (1-HWR()*(1-((1-ZMXISR)/(1-T/ZMXISR))**B))**(1/B)
793
794 IF (2*R.LE.C2+C1) Z(1)=T/Z(1)
795
796 ENDIF
797
798 ELSE
799
800 Z(1)=SQRT(T)
801
802 ENDIF
803
804 GAMWT=GAMWT*C2/Z(1)
805
806 $ /(1/Z(1)+(1-Z(1))**(B-1)+T/Z(1)**2*(1-T/Z(1))**(B-1))
807
808 ENDIF
809
810C---INCLUDE DISTRIBUTION FUNCTIONS
811
812 Z(2)=T/Z(1)
813
814 DO 10 I=1,2
815
816 IF (Z(I).GT.ZMXISR) THEN
817
818 Z(I)=1
819
820 GAMWT=GAMWT*(1-ZMXISR)**B*EXP(3*B/4)*(1-B**2*PIFAC**2/12)
821
822 ELSE
823
824 GAMWT=GAMWT*(B*(1-Z(I))**(B-1)*(1+Z(I)**2)/2
825
826 $ *EXP(B*Z(I)/2*(1+Z(I)/2))*(1-B**2*PIFAC**2/12)
827
828 $ +B**2/8*((1+Z(I))*((1+Z(I))**2+3*LOG(Z(I)))
829
830 $ -4*LOG(Z(I))/(1-Z(I))))
831
832 ENDIF
833
834 10 CONTINUE
835
836C---CHOOSE BOTH QSQ VALUES
837
838 DO 30 I=1,2
839
840 IF (Z(I).GT.ZMXISR .OR. COLISR) THEN
841
842 QSQ(I)=0
843
844 ELSE
845
846 J=3-I
847
848C---ACCORDING TO 1/(QSQ+QSQMIN) FROM 0 TO (1-Z)*(T/(Z+T))*QSQMAX
849
850 20 QSQ(I)=(((1-Z(I))*(T/(Z(I)+T))
851
852 $ *QSQMAX/QSQMIN+1)**HWR()-1)*QSQMIN
853
854C---AND REJECT TO QSQ/(QSQ+QSQMIN)**2
855
856 IF (HWR()*(QSQ(I)+QSQMIN).GT.QSQ(I)) GOTO 20
857
858 ENDIF
859
860 30 CONTINUE
861
862C---CHOOSE BOTH AZIMUTHS
863
864 PHI(1)=HWR()*2*PIFAC
865
866 PHI(2)=HWR()*2*PIFAC
867
868C---USE S-HAT PRESCRIPTION TO MODIFY Z VALUES
869
870 I=0
871
872 IF ((1-Z(1))*QSQ(1).GT.(1-Z(2))*QSQ(2)) I=1
873
874 IF ((1-Z(2))*QSQ(2).GT.(1-Z(1))*QSQ(1)) I=2
875
876 IF (I.GT.0) THEN
877
878 J=3-I
879
880 Z(I)=Z(I)+QSQ(I)/QSQMAX
881
882 IF (QSQ(J).GT.ZERO) THEN
883
884 Z(J)=((QSQ(I)*QSQMAX+QSQ(J)*QSQMAX
885
886 $ -QSQ(I)*QSQ(J))/QSQMAX**2+T)/Z(I)
887
888 C=COS(PHI(1)-PHI(2))*SQRT(QSQ(1)*QSQ(2))/QSQMAX
889
890 Z(J)=Z(J)+(-2*C**2*(1-Z(I))+2*C*SQRT((1-Z(I))
891
892 $ *(C**2*(1-Z(I))+Z(I)**2*(1-Z(J)))))/Z(I)**2
893
894 ENDIF
895
896 ENDIF
897
898 ELSEIF (IHEP.EQ.2) THEN
899
900C---EVERYTHING WAS GENERATED LAST TIME
901
902 ELSE
903
904C---ROUTINE CALLED UNEXPECTEDLY
905
906 CALL HWWARN('HWEISR',201,*999)
907
908 ENDIF
909
910C---IF Z IS TOO LARGE THERE IS NO EMISSION
911
912 IF (Z(IHEP).GT.ZMXISR) RETURN
913
914C---PUT NEW LEPTON IN EVENT RECORD
915
916 NHEP=NHEP+1
917
918 IDHW(NHEP)=IDHW(IHEP)
919
920 IDHEP(NHEP)=IDHEP(IHEP)
921
922 ISTHEP(NHEP)=3
923
924 JMOHEP(1,NHEP)=IHEP
925
926 JMOHEP(2,NHEP)=0
927
928 JDAHEP(1,NHEP)=0
929
930 JDAHEP(2,NHEP)=0
931
932 JDAHEP(1,IHEP)=NHEP
933
934C---AND OUTGOING PHOTON
935
936 NHEP=NHEP+1
937
938 IDHW(NHEP)=59
939
940 IDHEP(NHEP)=22
941
942 ISTHEP(NHEP)=1
943
944 JMOHEP(1,NHEP)=IHEP
945
946 JMOHEP(2,NHEP)=0
947
948 JDAHEP(1,NHEP)=0
949
950 JDAHEP(2,NHEP)=0
951
952 JDAHEP(2,IHEP)=NHEP
953
954C---RECONSTRUCT PHOTON KINEMATICS (Z IS LIGHT-CONE MOMENTUM FRACTION)
955
956 PHEP(1,NHEP)=SQRT(QSQ(IHEP)*(1-Z(IHEP)))*COS(PHI(IHEP))
957
958 PHEP(2,NHEP)=SQRT(QSQ(IHEP)*(1-Z(IHEP)))*SIN(PHI(IHEP))
959
960 PHEP(3,NHEP)=(1-Z(IHEP))*PHEP(4,IHEP)-QSQ(IHEP)/(4*PHEP(4,IHEP))
961
962 IF (IHEP.EQ.2) PHEP(3,NHEP)=-PHEP(3,NHEP)
963
964 PHEP(4,NHEP)=(1-Z(IHEP))*PHEP(4,IHEP)+QSQ(IHEP)/(4*PHEP(4,IHEP))
965
966 PHEP(5,NHEP)=0
967
968C---AND LEPTON
969
970 CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,NHEP),PHEP(1,NHEP-1))
971
972 CALL HWUMAS(PHEP(1,NHEP-1))
973
974C---UPDATE OVERALL CM FRAME
975
976 JMOHEP(IHEP,3)=NHEP-1
977
978 CALL HWVDIF(4,PHEP(1,3),PHEP(1,IHEP),PHEP(1,3))
979
980 CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,3),PHEP(1,3))
981
982 CALL HWUMAS(PHEP(1,3))
983
984 999 END
985
986CDECK ID>, HWEONE.
987
988*CMZ :- -26/04/91 11.11.55 by Bryan Webber
989
990*-- Author : Bryan Webber
991
992C-----------------------------------------------------------------------
993
994 SUBROUTINE HWEONE
995
996C-----------------------------------------------------------------------
997
998C SETS UP 2->1 (COLOUR SINGLET) HARD SUBPROCESS
999
1000C-----------------------------------------------------------------------
1001
1002 INCLUDE 'HERWIG61.INC'
1003
1004 DOUBLE PRECISION PA
1005
1006 INTEGER ICMF,I,IBM,IHEP
1007
1008C---INCOMING LINES
1009
1010 ICMF=NHEP+3
1011
1012 DO 15 I=1,2
1013
1014 IBM=I
1015
1016C---FIND BEAM AND TARGET
1017
1018 IF (JDAHEP(1,I).NE.0) IBM=JDAHEP(1,I)
1019
1020 IHEP=NHEP+I
1021
1022 IDHW(IHEP)=IDN(I)
1023
1024 IDHEP(IHEP)=IDPDG(IDN(I))
1025
1026 ISTHEP(IHEP)=110+I
1027
1028 JMOHEP(1,IHEP)=ICMF
1029
1030 JMOHEP(I,ICMF)=IHEP
1031
1032 JDAHEP(1,IHEP)=ICMF
1033
1034C---SPECIAL - IF INCOMING PARTON IS INCOMING BEAM THEN COPY IT
1035
1036 IF (XX(I).EQ.ONE.AND.IDHW(IBM).EQ.IDN(I)) THEN
1037
1038 CALL HWVEQU(5,PHEP(1,IBM),PHEP(1,IHEP))
1039
1040 IF (I.EQ.2) PHEP(3,IHEP)=-PHEP(3,IHEP)
1041
1042 ELSE
1043
1044 PHEP(1,IHEP)=0.
1045
1046 PHEP(2,IHEP)=0.
1047
1048 PHEP(5,IHEP)=RMASS(IDN(I))
1049
1050 PA=XX(I)*(PHEP(4,IBM)+ABS(PHEP(3,IBM)))
1051
1052 PHEP(4,IHEP)=0.5*(PA+PHEP(5,IHEP)**2/PA)
1053
1054 PHEP(3,IHEP)=PA-PHEP(4,IHEP)
1055
1056 ENDIF
1057
1058 15 CONTINUE
1059
1060 PHEP(3,NHEP+2)=-PHEP(3,NHEP+2)
1061
1062C---HARD CENTRE OF MASS
1063
1064 IDHW(ICMF)=IDCMF
1065
1066 IDHEP(ICMF)=IDPDG(IDCMF)
1067
1068 ISTHEP(ICMF)=110
1069
1070 CALL HWVSUM(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,ICMF))
1071
1072 CALL HWUMAS(PHEP(1,ICMF))
1073
1074C---SET UP COLOUR STRUCTURE LABELS
1075
1076 JMOHEP(2,NHEP+1)=NHEP+2
1077
1078 JDAHEP(2,NHEP+1)=NHEP+2
1079
1080 JMOHEP(2,NHEP+2)=NHEP+1
1081
1082 JDAHEP(2,NHEP+2)=NHEP+1
1083
1084 JDAHEP(1,NHEP+3)=NHEP+3
1085
1086 JDAHEP(2,NHEP+3)=NHEP+3
1087
1088 NHEP=NHEP+3
1089
1090 999 END
1091
1092CDECK ID>, HWEPRO.
1093
1094*CMZ :- -01/04/99 19.41.18 by Mike Seymour
1095
1096*-- Author : Bryan Webber
1097
1098C-----------------------------------------------------------------------
1099
1100 SUBROUTINE HWEPRO
1101
1102C-----------------------------------------------------------------------
1103
1104C WHEN NEVHEP=0, CHOOSES X VALUES AND FINDS WEIGHT FOR PROCESS IPROC
1105
1106C OTHERWISE, CHOOSES AND LOADS ALL VARIABLES FOR HARD PROCESS
1107
1108C-----------------------------------------------------------------------
1109
1110 INCLUDE 'HERWIG61.INC'
1111
1112 DOUBLE PRECISION HWR
1113
1114 EXTERNAL HWR
1115
1116 IF (IERROR.NE.0) RETURN
1117
1118C---ROUTINE LOOPS BACK TO HERE IF GENERATED WEIGHT WAS NOT ACCEPTED
1119
1120 10 GENEV=.FALSE.
1121
1122C---FSTWGT IS .TRUE. DURING FIRST CALL TO HARD PROCESS ROUTINE
1123
1124 FSTWGT=NWGTS.EQ.0
1125
1126C---FSTEVT IS .TRUE. THROUGHOUT THE FIRST EVENT
1127
1128 FSTEVT=NEVHEP.EQ.1
1129
1130C---SET COLOUR CORRECTION TO FALSE
1131
1132 COLUPD = .FALSE.
1133
1134 HRDCOL(1,1)=0
1135
1136 HRDCOL(1,3)=0
1137
1138C---SET UP INITIAL STATE
1139
1140 NHEP=1
1141
1142 ISTHEP(NHEP)=101
1143
1144 PHEP(1,NHEP)=0.
1145
1146 PHEP(2,NHEP)=0.
1147
1148 PHEP(3,NHEP)=PBEAM1
1149
1150 PHEP(4,NHEP)=EBEAM1
1151
1152 PHEP(5,NHEP)=RMASS(IPART1)
1153
1154 JMOHEP(1,NHEP)=0
1155
1156 JMOHEP(2,NHEP)=0
1157
1158 JDAHEP(1,NHEP)=0
1159
1160 JDAHEP(2,NHEP)=0
1161
1162 IDHW(NHEP)=IPART1
1163
1164 IDHEP(NHEP)=IDPDG(IPART1)
1165
1166 NHEP=NHEP+1
1167
1168 ISTHEP(NHEP)=102
1169
1170 PHEP(1,NHEP)=0.
1171
1172 PHEP(2,NHEP)=0.
1173
1174 PHEP(3,NHEP)=-PBEAM2
1175
1176 PHEP(4,NHEP)=EBEAM2
1177
1178 PHEP(5,NHEP)=RMASS(IPART2)
1179
1180 JMOHEP(1,NHEP)=0
1181
1182 JMOHEP(2,NHEP)=0
1183
1184 JDAHEP(1,NHEP)=0
1185
1186 JDAHEP(2,NHEP)=0
1187
1188 IDHW(NHEP)=IPART2
1189
1190 IDHEP(NHEP)=IDPDG(IPART2)
1191
1192C---NEXT ENTRY IS OVERALL CM FRAME
1193
1194 NHEP=NHEP+1
1195
1196 IDHW(NHEP)=14
1197
1198 IDHEP(NHEP)=0
1199
1200 ISTHEP(NHEP)=103
1201
1202 JMOHEP(1,NHEP)=NHEP-2
1203
1204 JMOHEP(2,NHEP)=NHEP-1
1205
1206 JDAHEP(1,NHEP)=0
1207
1208 JDAHEP(2,NHEP)=0
1209
1210 CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP-2),PHEP(1,NHEP))
1211
1212 CALL HWUMAS(PHEP(1,NHEP))
1213
1214C Select a primary interaction point
1215
1216 IF (PIPSMR) THEN
1217
1218 CALL HWRPIP
1219
1220 ELSE
1221
1222 CALL HWVZRO(4,VTXPIP)
1223
1224 ENDIF
1225
1226 CALL HWVEQU(3,VTXPIP,VHEP(1,NHEP))
1227
1228 VHEP(4,NHEP)=0.0
1229
1230C---GENERATE PHOTONS (WEIZSACKER-WILLIAMS APPROX)
1231
1232C FOR HADRONIC PROCESSES WITH LEPTON BEAMS
1233
1234 GAMWT=ONE
1235
1236 IF (IPRO.GT.10.AND.IPRO.LT.90) THEN
1237
1238 IF (ABS(IDHEP(1)).EQ.11.OR.ABS(IDHEP(1)).EQ.13)
1239
1240 & CALL HWEGAM(1,ZERO, ONE,.FALSE.)
1241
1242 IF (ABS(IDHEP(2)).EQ.11.OR.ABS(IDHEP(2)).EQ.13)
1243
1244 & CALL HWEGAM(2,ZERO, ONE,.FALSE.)
1245
1246 ELSEIF (IPRO.GE.90) THEN
1247
1248 IF (ABS(IDHEP(2)).EQ.11.OR.ABS(IDHEP(2)).EQ.13)
1249
1250 & CALL HWEGAM(2,ZERO, ONE,.FALSE.)
1251
1252 ENDIF
1253
1254C---GENERATE ISR PHOTONS FOR LEPTONIC PROCESSES
1255
1256 IF (IPRO.LT.10) THEN
1257
1258 CALL HWEISR(1)
1259
1260 CALL HWEISR(2)
1261
1262 ENDIF
1263
1264C---IF USER LIMITS WERE TOO TIGHT, MIGHT NOT BE ANY PHASE-SPACE
1265
1266 IF (GAMWT.LE.ZERO) GOTO 30
1267
1268C---IF CMF HAS ACQUIRED A TRANSVERSE BOOST, OR USER REQUESTS IT ANYWAY,
1269
1270C BOOST EVENT RECORD BACK TO CMF
1271
1272 IF (PHEP(1,3)**2+PHEP(2,3)**2.GT.ZERO .OR. USECMF) CALL HWUBST(1)
1273
1274C---ROUTINE LOOPS BACK TO HERE IF GENERATED WEIGHT WAS ACCEPTED
1275
1276 20 CONTINUE
1277
1278C---IPRO=MOD(IPROC/100,100)
1279
1280 IF (IPRO.EQ.1) THEN
1281
1282 IF (IPROC.LT.110.OR.IPROC.GE.120) THEN
1283
1284C--- E+E- -> Q-QBAR OR L-LBAR
1285
1286 CALL HWHEPA
1287
1288 ELSE
1289
1290C--- E+E- -> Q-QBAR-GLUON
1291
1292 CALL HWHEPG
1293
1294 ENDIF
1295
1296 ELSEIF (IPRO.EQ.2) THEN
1297
1298C--- E+E- -> W+ W-
1299
1300 CALL HWHEWW
1301
1302 ELSEIF (IPRO.EQ.3) THEN
1303
1304C---E+E- -> Z H
1305
1306 CALL HWHIGZ
1307
1308 ELSEIF (IPRO.EQ.4) THEN
1309
1310C---E+E- -> NUEB NUE H
1311
1312 CALL HWHIGW
1313
1314 ELSEIF (IPRO.EQ.5 .AND. IPROC.LT.550) THEN
1315
1316C---EE -> EE GAMGAM -> EE FFBAR/WW
1317
1318 CALL HWHEGG
1319
1320 ELSEIF (IPRO.EQ.5) THEN
1321
1322C---EE -> ENU GAMW -> ENU FF'BAR/WZ
1323
1324 CALL HWHEGW
1325
1326 ELSEIF (IPRO.EQ.6) THEN
1327
1328C---EE -> FOUR JETS
1329
1330 CALL HWH4JT
1331
1332 ELSEIF (IPRO.EQ.13) THEN
1333
1334C---GAMMA/Z0/Z' DRELL-YAN PROCESS
1335
1336 CALL HWHDYP
1337
1338 ELSEIF (IPRO.EQ.14) THEN
1339
1340C---W+/- PRODUCTION VIA DRELL-YAN PROCESS
1341
1342 CALL HWHWPR
1343
1344 ELSEIF (IPRO.EQ.15) THEN
1345
1346C---QCD HARD 2->2 PROCESSES
1347
1348 CALL HWHQCD
1349
1350 ELSEIF (IPRO.EQ.16) THEN
1351
1352C---HIGGS PRODUCTION VIA GLUON FUSION
1353
1354 CALL HWHIGS
1355
1356 ELSEIF (IPRO.EQ.17) THEN
1357
1358C---QCD HEAVY FLAVOUR PRODUCTION
1359
1360 CALL HWHHVY
1361
1362 ELSEIF (IPRO.EQ.18) THEN
1363
1364C---QCD DIRECT PHOTON + JET PRODUCTION
1365
1366 CALL HWHPHO
1367
1368 ELSEIF (IPRO.EQ.19) THEN
1369
1370C---HIGGS PRODUCTION VIA W FUSION
1371
1372 CALL HWHIGW
1373
1374 ELSEIF (IPRO.EQ.20) THEN
1375
1376C---TOP PRODUCTION FROM W EXCHANGE
1377
1378 CALL HWHWEX
1379
1380 ELSEIF (IPRO.EQ.21) THEN
1381
1382C---VECTOR BOSON + JET PRODUCTION
1383
1384 CALL HWHV1J
1385
1386 ELSEIF (IPRO.EQ.22) THEN
1387
1388C QCD direct photon pair production
1389
1390 CALL HWHPH2
1391
1392 ELSEIF (IPRO.EQ.23) THEN
1393
1394C QCD Higgs plus jet production
1395
1396 CALL HWHIGJ
1397
1398 ELSEIF (IPRO.EQ.24) THEN
1399
1400C---COLOUR-SINGLET EXCHANGE
1401
1402 CALL HWHSNG
1403
1404 ELSEIF (IPRO.EQ.30) THEN
1405
1406C---HADRON-HADRON SUSY PROCESSES
1407
1408 CALL HWHSSP
1409
1410 ELSEIF(IPRO.EQ.40.OR.IPRO.EQ.41) THEN
1411
1412C---HADRON-HADRON R-PARITY VIOLATING SUSY PROCESSES
1413
1414 CALL HWHRSP
1415
1416 ELSEIF (IPRO.EQ.50) THEN
1417
1418C Point-like photon two-jet production
1419
1420 CALL HWHPPT
1421
1422 ELSEIF (IPRO.EQ.51) THEN
1423
1424C Point-like photon/QCD heavy flavour pair production
1425
1426 CALL HWHPPH
1427
1428 ELSEIF (IPRO.EQ.52) THEN
1429
1430C Point-like photon/QCD heavy flavour single excitation
1431
1432 CALL HWHPPE
1433
1434 ELSEIF (IPRO.EQ.53) THEN
1435
1436C Compton scattering of point-like photon and (anti)quark
1437
1438 CALL HWHPQS
1439
1440 ELSEIF (IPRO.EQ.55) THEN
1441
1442C Point-like photon/higher twist meson production
1443
1444 CALL HWHPPM
1445
1446 ELSEIF (IPRO.EQ.60) THEN
1447
1448C---QPM GAMMA-GAMMA-->QQBAR
1449
1450 CALL HWHQPM
1451
1452 ELSEIF (IPRO.GE.70.AND.IPRO.LE.79) THEN
1453
1454C---BARYON-NUMBER VIOLATION, AND OTHER MULTI-W PRODUCTION PROCESSES
1455
1456 CALL HVHBVI
1457
1458 ELSEIF (IPRO.EQ.80) THEN
1459
1460C---MINIMUM-BIAS: NO HARD SUBPROCESS
1461
1462C FIND WEIGHT
1463
1464 CALL HWMWGT
1465
1466 ELSEIF (IPRO.EQ.90) THEN
1467
1468C---DEEP INELASTIC
1469
1470 CALL HWHDIS
1471
1472 ELSEIF(IPRO.EQ.91) THEN
1473
1474C---BOSON - GLUON(QUARK) FUSION --> ANTIQUARK(GLUON) + QUARK
1475
1476 CALL HWHBGF
1477
1478 ELSEIF(IPRO.EQ.92) THEN
1479
1480C---DEEP INELASTIC WITH EXTRA JET: OBSOLETE PROCESS
1481
1482 WRITE (6,40)
1483
1484 40 FORMAT (1X,' IPROC=92** is no longer supported.'
1485
1486 & /1X,' Please use IPROC=91** instead.')
1487
1488 CALL HWWARN('HWEPRO',500,*999)
1489
1490 ELSEIF(IPRO.EQ.95) THEN
1491
1492C---HIGGS PRODUCTION VIA W FUSION IN E P
1493
1494 CALL HWHIGW
1495
1496 ELSE
1497
1498C---UNKNOWN PROCESS
1499
1500 CALL HWWARN('HWEPRO',102,*999)
1501
1502 ENDIF
1503
1504 30 IF (GENEV) THEN
1505
1506 IF (NOWGT) EVWGT=AVWGT
1507
1508 ISTAT=10
1509
1510 RETURN
1511
1512 ELSE
1513
1514C---IF AN EVENT IS CANCELLED BEFORE IT IS GENERATED, GIVE IT ZERO WEIGHT
1515
1516 IF (IERROR.NE.0) THEN
1517
1518 EVWGT=ZERO
1519
1520 IERROR=0
1521
1522 ENDIF
1523
1524 EVWGT=EVWGT*GAMWT
1525
1526 NWGTS=NWGTS+1
1527
1528 WGTSUM=WGTSUM+EVWGT
1529
1530 WSQSUM=WSQSUM+EVWGT**2
1531
1532 IF (EVWGT.GT.WBIGST) THEN
1533
1534 WBIGST=EVWGT
1535
1536 IF (NOWGT.AND.WBIGST.GT.WGTMAX) THEN
1537
1538 IF (NEVHEP.NE.0) CALL HWWARN('HWEPRO',1,*999)
1539
1540 WGTMAX=WBIGST*1.1
1541
1542 WRITE (6,99) WGTMAX
1543
1544 ENDIF
1545
1546 ELSEIF (EVWGT.LT.ZERO) THEN
1547
1548 IF (EVWGT.LT.-1.D-9) CALL HWWARN('HWEPRO',3,*999)
1549
1550 EVWGT=0.
1551
1552 ENDIF
1553
1554 IF (NEVHEP.NE.0) THEN
1555
1556C---LOW EFFICIENCY WARNINGS:
1557
1558C WARN AT 10*EFFMIN, STOP AT EFFMIN
1559
1560 IF (10*EFFMIN*NWGTS.GT.NEVHEP) THEN
1561
1562 IF (EFFMIN*NWGTS.GT.NEVHEP) CALL HWWARN('HWEPRO',200,*999)
1563
1564 IF (EFFMIN.GT.ZERO) THEN
1565
1566 IF (MOD(NWGTS,INT(10/EFFMIN)).EQ.0) THEN
1567
1568 CALL HWWARN('HWEPRO',2,*999)
1569
1570 WRITE (6,98) WGTMAX
1571
1572 ENDIF
1573
1574 ENDIF
1575
1576 ENDIF
1577
1578 IF (NOWGT) THEN
1579
1580 GENEV=EVWGT.GT.WGTMAX*HWR()
1581
1582 ELSE
1583
1584 GENEV=EVWGT.NE.ZERO
1585
1586 ENDIF
1587
1588 IF (GENEV) GOTO 20
1589
1590 GOTO 10
1591
1592 ENDIF
1593
1594 ENDIF
1595
1596 98 FORMAT(10X,' MAXIMUM WEIGHT =',1PG24.16)
1597
1598 99 FORMAT(10X,'NEW MAXIMUM WEIGHT =',1PG24.16)
1599
1600 999 END
1601
1602CDECK ID>, HWETWO.
1603
1604*CMZ :- -26/04/91 11.11.55 by Bryan Webber
1605
1606*-- Author : Bryan Webber
1607
1608C-----------------------------------------------------------------------
1609
1610 SUBROUTINE HWETWO
1611
1612C-----------------------------------------------------------------------
1613
1614C SETS UP 2->2 HARD SUBPROCESS
1615
1616C-----------------------------------------------------------------------
1617
1618 INCLUDE 'HERWIG61.INC'
1619
1620 DOUBLE PRECISION HWUMBW,HWUPCM,PA,PCM
1621
1622 INTEGER ICMF,IBM,I,J,K,IHEP,NTRY
1623
1624 EXTERNAL HWUPCM
1625
1626C---INCOMING LINES
1627
1628 ICMF=NHEP+3
1629
1630 DO 15 I=1,2
1631
1632 IBM=I
1633
1634C---FIND BEAM AND TARGET
1635
1636 IF (JDAHEP(1,I).NE.0) IBM=JDAHEP(1,I)
1637
1638 IHEP=NHEP+I
1639
1640 IDHW(IHEP)=IDN(I)
1641
1642 IDHEP(IHEP)=IDPDG(IDN(I))
1643
1644 ISTHEP(IHEP)=110+I
1645
1646 JMOHEP(1,IHEP)=ICMF
1647
1648 JMOHEP(I,ICMF)=IHEP
1649
1650 JDAHEP(1,IHEP)=ICMF
1651
1652C---SPECIAL - IF INCOMING PARTON IS INCOMING BEAM THEN COPY IT
1653
1654 IF (XX(I).EQ.ONE.AND.IDHW(IBM).EQ.IDN(I)) THEN
1655
1656 CALL HWVEQU(5,PHEP(1,IBM),PHEP(1,IHEP))
1657
1658 IF (I.EQ.2) PHEP(3,IHEP)=-PHEP(3,IHEP)
1659
1660 ELSE
1661
1662 PHEP(1,IHEP)=0.
1663
1664 PHEP(2,IHEP)=0.
1665
1666 PHEP(5,IHEP)=RMASS(IDN(I))
1667
1668 PA=XX(I)*(PHEP(4,IBM)+ABS(PHEP(3,IBM)))
1669
1670 PHEP(4,IHEP)=0.5*(PA+PHEP(5,IHEP)**2/PA)
1671
1672 PHEP(3,IHEP)=PA-PHEP(4,IHEP)
1673
1674 ENDIF
1675
1676 15 CONTINUE
1677
1678 PHEP(3,NHEP+2)=-PHEP(3,NHEP+2)
1679
1680C---HARD CENTRE OF MASS
1681
1682 IDHW(ICMF)=IDCMF
1683
1684 IDHEP(ICMF)=IDPDG(IDCMF)
1685
1686 ISTHEP(ICMF)=110
1687
1688 CALL HWVSUM(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,ICMF))
1689
1690 CALL HWUMAS(PHEP(1,ICMF))
1691
1692C---OUTGOING LINES
1693
1694 NTRY=0
1695
1696 19 CONTINUE
1697
1698 DO 20 I=3,4
1699
1700 IHEP=NHEP+I+1
1701
1702 IDHW(IHEP)=IDN(I)
1703
1704 IDHEP(IHEP)=IDPDG(IDN(I))
1705
1706 ISTHEP(IHEP)=110+I
1707
1708 JMOHEP(1,IHEP)=ICMF
1709
1710 JDAHEP(I-2,ICMF)=IHEP
1711
1712 20 PHEP(5,IHEP)=HWUMBW(IDN(I))
1713
1714 PCM=HWUPCM(PHEP(5,NHEP+3),PHEP(5,NHEP+4),PHEP(5,NHEP+5))
1715
1716 IF (PCM.LT.ZERO) THEN
1717
1718 NTRY=NTRY+1
1719
1720 IF (NTRY.LE.NETRY) GO TO 19
1721
1722 CALL HWWARN('HWETWO',103,*999)
1723
1724 ENDIF
1725
1726 IHEP=NHEP+4
1727
1728 PHEP(4,IHEP)=SQRT(PCM**2+PHEP(5,IHEP)**2)
1729
1730 PHEP(3,IHEP)=PCM*COSTH
1731
1732 PHEP(1,IHEP)=SQRT((PCM+PHEP(3,IHEP))*(PCM-PHEP(3,IHEP)))
1733
1734 CALL HWRAZM(PHEP(1,IHEP),PHEP(1,IHEP),PHEP(2,IHEP))
1735
1736 CALL HWULOB(PHEP(1,NHEP+3),PHEP(1,IHEP),PHEP(1,IHEP))
1737
1738 CALL HWVDIF(4,PHEP(1,NHEP+3),PHEP(1,IHEP),PHEP(1,NHEP+5))
1739
1740C---SET UP COLOUR STRUCTURE LABELS
1741
1742 DO 30 I=1,4
1743
1744 J=I
1745
1746 IF (J.GT.2) J=J+1
1747
1748 K=ICO(I)
1749
1750 IF (K.GT.2) K=K+1
1751
1752 JMOHEP(2,NHEP+J)=NHEP+K
1753
1754 30 JDAHEP(2,NHEP+K)=NHEP+J
1755
1756 NHEP=NHEP+5
1757
1758 999 END
1759
1760CDECK ID>, HWH4JT.
1761
1762*CMZ :- -01/04/99 19.47.55 by Mike Seymour
1763
1764*-- Author : Ian Knowles
1765
1766C-----------------------------------------------------------------------
1767
1768 SUBROUTINE HWH4JT
1769
1770C-----------------------------------------------------------------------
1771
1772C Four jet production in e^+e^- annihilation: qqbar+gg & qqbar+qqbar
1773
1774C IOP4JT controls the treatment of the colour flow interference term
1775
1776C qqbar-gg case:
1777
1778C IOP4JT(1)=0 neglect, =1 extreme 2341; =2 extreme 3421
1779
1780C qqbar-qqbar (identical quark flavour) case:
1781
1782C IOP4JT(2)=0 neglect, =1 extreme 4123; =2 extreme 2143
1783
1784C
1785
1786C Matrix elements based on Ellis Ross & Terrano and Catani & Seymour
1787
1788C
1789
1790C WARNING: Phase space factor inaccurate for JADE y_cut > 0.14.
1791
1792C-----------------------------------------------------------------------
1793
1794 INCLUDE 'HERWIG61.INC'
1795
1796 INTEGER LM,LP,IQK,I,J,IDMN,IDMX,ID1,ID2,IST(4)
1797
1798 DOUBLE PRECISION HWR,HWUALF,HWUAEM,HWULDO,HWH4J1,HWH4J2,
1799
1800 & HWH4J4,HWH4J5,HWH4J6,HWH4J7,QNOW,Q2NOW,QLST,SCUT,PSFAC,FACT,X,
1801
1802 & COLA,COLB,COLC,CLF(7,6),P12,P13,P14,P23,P24,P34,FACTR,EP1,EP2,
1803
1804 & EP3,EP4,GG1,GG2,GG12,GG3,GG13,GG23,GGINT,WTGG,QQ,QP,QQINT,QQ1,
1805
1806 & QQ2,WTQQ,WTQP,HCS,WTAB,WTBA,WTOT,RCS
1807
1808 $ ,EF,QF,E(4)
1809
1810 LOGICAL INCLQG(6),INCLQQ(6,6),ORIENT
1811
1812 EXTERNAL HWR,HWUALF,HWUAEM,HWULDO,HWH4J1,HWH4J2,HWH4J4,
1813
1814 & HWH4J5,HWH4J6,HWH4J7
1815
1816 SAVE HCS,QLST,WTQP,WTQQ,WTGG,FACTR,COLA,COLB,COLC,IDMN,IDMX,
1817
1818 & CLF,GG1,GG2,GGINT,INCLQG,INCLQQ,LM,LP,QQ1,QQ2,QQINT,FACT,ORIENT,
1819
1820 & Q2NOW,SCUT
1821
1822 DATA QLST,IST/-1.,113,114,114,114/
1823
1824C
1825
1826 IF (GENEV) THEN
1827
1828 RCS=HCS*HWR()
1829
1830 ELSE
1831
1832 IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWH4JT',100,*999)
1833
1834 QNOW=PHEP(5,3)
1835
1836 IF (QNOW.NE.QLST) THEN
1837
1838 QLST=QNOW
1839
1840 Q2NOW=QNOW**2
1841
1842 SCUT=Y4JT*Q2NOW
1843
1844C Calculate allowed fraction of Phase Space using parameterization
1845
1846 IF (DURHAM) THEN
1847
1848 PSFAC=(1.-6.*Y4JT)**5.50*(1.-173.3*Y4JT*(1.-247.3*Y4JT
1849
1850 & *(1.+148.3*Y4JT*(1.+3.913*Y4JT))))
1851
1852 & /(1.-8.352*Y4JT*(1.-1102.*Y4JT
1853
1854 & *(1.+1603.*Y4JT*(1.+22.99*Y4JT))))
1855
1856 ELSE
1857
1858 PSFAC=(1.-6.*Y4JT)**4.62*(1.-44.72*Y4JT*(1.-176.0*Y4JT
1859
1860 & *(1.+102.9*Y4JT*(1.-6.579*Y4JT))))
1861
1862 & /(1.-3.392*Y4JT*(1.-946.5*Y4JT
1863
1864 & *(1.+423.4*Y4JT*(1.-3.971*Y4JT))))
1865
1866 ENDIF
1867
1868 FACT=GEV2NB*HWUAEM(Q2NOW)**2*CFFAC*FLOAT(NCOLO)*PSFAC
1869
1870 & /(THREE*16*PIFAC)
1871
1872 COLA=CFFAC
1873
1874 COLB=CFFAC-HALF*CAFAC
1875
1876 COLC=HALF
1877
1878 LM=1
1879
1880 IF (JDAHEP(1,LM).NE.0) LM=JDAHEP(1,LM)
1881
1882 LP=2
1883
1884 IF (JDAHEP(1,LP).NE.0) LP=JDAHEP(1,LP)
1885
1886 IQK=MOD(IPROC,10)
1887
1888 IF (IQK.NE.0) THEN
1889
1890 IDMN=IQK
1891
1892 IDMX=IQK
1893
1894 ELSE
1895
1896 IDMN=1
1897
1898 IDMX=6
1899
1900 ENDIF
1901
1902 DO 10 I=1,6
1903
1904 CALL HWUCFF(11,I,Q2NOW,CLF(1,I))
1905
1906 IF (QNOW.GT.TWO*(RMASS(I)+RMASS(13))) THEN
1907
1908 INCLQG(I)=.TRUE.
1909
1910 ELSE
1911
1912 INCLQG(I)=.FALSE.
1913
1914 ENDIF
1915
1916 DO 10 J=I,6
1917
1918 IF (QNOW.GT.TWO*(RMASS(I)+RMASS(J ))) THEN
1919
1920 INCLQQ(I,J)=.TRUE.
1921
1922 INCLQQ(J,I)=.TRUE.
1923
1924 ELSE
1925
1926 INCLQQ(I,J)=.FALSE.
1927
1928 INCLQQ(J,I)=.FALSE.
1929
1930 ENDIF
1931
1932 10 CONTINUE
1933
1934 IF (MOD(IPROC/10,10).EQ.5) THEN
1935
1936 ORIENT=.FALSE.
1937
1938 ELSE
1939
1940 ORIENT=.TRUE.
1941
1942 ENDIF
1943
1944 ENDIF
1945
1946C Generate phase space point and check it passes cuts
1947
1948 CALL HWVEQU(5,PHEP(1,3),PHEP(1,NHEP+1))
1949
1950 DO 20 I=2,5
1951
1952 20 PHEP(5,NHEP+I)=0.
1953
1954 30 CALL HWDFOR(PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,NHEP+3),
1955
1956 & PHEP(1,NHEP+4),PHEP(1,NHEP+5))
1957
1958 IF (DURHAM) THEN
1959
1960 P12=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+3))
1961
1962 X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+3),
1963
1964 & PHEP(4,NHEP+3)/PHEP(4,NHEP+2))*P12
1965
1966 IF (X.GT.SCUT) THEN
1967
1968 P13=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+4))
1969
1970 X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+4),
1971
1972 & PHEP(4,NHEP+4)/PHEP(4,NHEP+2))*P13
1973
1974 IF (X.GT.SCUT) THEN
1975
1976 P14=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+5))
1977
1978 X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+5),
1979
1980 & PHEP(4,NHEP+5)/PHEP(4,NHEP+2))*P14
1981
1982 IF (X.GT.SCUT) THEN
1983
1984 P23=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+4))
1985
1986 X=MIN(PHEP(4,NHEP+3)/PHEP(4,NHEP+4),
1987
1988 & PHEP(4,NHEP+4)/PHEP(4,NHEP+3))*P23
1989
1990 IF (X.GT.SCUT) THEN
1991
1992 P24=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+5))
1993
1994 X=MIN(PHEP(4,NHEP+3)/PHEP(4,NHEP+5),
1995
1996 & PHEP(4,NHEP+5)/PHEP(4,NHEP+3))*P24
1997
1998 IF (X.GT.SCUT) THEN
1999
2000 P34=2*HWULDO(PHEP(1,NHEP+4),PHEP(1,NHEP+5))
2001
2002 X=MIN(PHEP(4,NHEP+4)/PHEP(4,NHEP+5),
2003
2004 & PHEP(4,NHEP+5)/PHEP(4,NHEP+4))*P34
2005
2006 IF (X.GT.SCUT) GOTO 40
2007
2008 ENDIF
2009
2010 ENDIF
2011
2012 ENDIF
2013
2014 ENDIF
2015
2016 ENDIF
2017
2018 ELSE
2019
2020 P12=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+3))
2021
2022 IF (P12.GT.SCUT) THEN
2023
2024 P13=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+4))
2025
2026 IF (P13.GT.SCUT) THEN
2027
2028 P14=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+5))
2029
2030 IF (P14.GT.SCUT) THEN
2031
2032 P23=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+4))
2033
2034 IF (P23.GT.SCUT) THEN
2035
2036 P24=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+5))
2037
2038 IF (P24.GT.SCUT) THEN
2039
2040 P34=2*HWULDO(PHEP(1,NHEP+4),PHEP(1,NHEP+5))
2041
2042 IF (P34.GT.SCUT) GOTO 40
2043
2044 ENDIF
2045
2046 ENDIF
2047
2048 ENDIF
2049
2050 ENDIF
2051
2052 ENDIF
2053
2054 ENDIF
2055
2056C Failed cuts retry
2057
2058 GOTO 30
2059
2060C Passed cuts: calculate contributions to Matrix Elements
2061
2062 40 EMSCA=SQRT(MIN(P12,P13,P14,P23,P24,P34))
2063
2064 FACTR=FACT*HWUALF(1,EMSCA)**2
2065
2066 IF (ORIENT) THEN
2067
2068 QF=HWULDO(PHEP(1,LP),PHEP(1,3))
2069
2070 EF=Q2NOW/(2*SQRT(QF**2-HWULDO(PHEP(1,LP),PHEP(1,LP))*Q2NOW))
2071
2072 QF=HALF-EF*QF/Q2NOW
2073
2074 DO I=1,4
2075
2076 E(I)=EF*PHEP(I,LP)+QF*PHEP(I,3)
2077
2078 ENDDO
2079
2080 EP1=HWULDO(E,PHEP(1,NHEP+2))
2081
2082 EP2=HWULDO(E,PHEP(1,NHEP+3))
2083
2084 EP3=HWULDO(E,PHEP(1,NHEP+4))
2085
2086 EP4=HWULDO(E,PHEP(1,NHEP+5))
2087
2088 ENDIF
2089
2090C q-qbar-g-g
2091
2092 GG1=HWH4J1(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2093
2094 & +HWH4J1(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2095
2096 GG2=HWH4J1(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2097
2098 & +HWH4J1(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2099
2100 GG12=HWH4J2(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2101
2102 & +HWH4J2(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2103
2104 & +HWH4J2(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2105
2106 & +HWH4J2(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2107
2108 GG3=HWH4J4(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2109
2110 & +HWH4J4(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2111
2112 GG13=GG3+HWH4J5(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2113
2114 & +HWH4J5(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2115
2116 GG23=GG3+HWH4J5(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2117
2118 & +HWH4J5(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2119
2120C Add up weights
2121
2122 GG1 =COLA*(GG1 +GG13)
2123
2124 GG2 =COLA*(GG2 +GG23)
2125
2126 GGINT=COLB*(GG12-GG13-GG23)
2127
2128 WTGG=FACTR*(GG1+GG2+GGINT)
2129
2130C q-qbar-q-qbar
2131
2132 QP=HWH4J6(P13,P12,P14,P23,P34,P24,EP1,EP3,EP2,EP4,ORIENT)
2133
2134 & +HWH4J6(P24,P12,P23,P14,P34,P13,EP2,EP4,EP1,EP3,ORIENT)
2135
2136 & +HWH4J6(P13,P34,P23,P14,P12,P24,EP3,EP1,EP4,EP2,ORIENT)
2137
2138 & +HWH4J6(P24,P34,P14,P23,P12,P13,EP4,EP2,EP3,EP1,ORIENT)
2139
2140 QQ=HWH4J6(P13,P23,P34,P12,P14,P24,EP3,EP1,EP2,EP4,ORIENT)
2141
2142 & +HWH4J6(P24,P23,P12,P34,P14,P13,EP2,EP4,EP3,EP1,ORIENT)
2143
2144 & +HWH4J6(P13,P14,P12,P34,P23,P24,EP1,EP3,EP4,EP2,ORIENT)
2145
2146 & +HWH4J6(P24,P14,P34,P12,P23,P13,EP4,EP2,EP1,EP3,ORIENT)
2147
2148 QQINT=HWH4J7(P13,P12,P14,P23,P34,P24,EP1,EP3,EP2,EP4,ORIENT)
2149
2150 & +HWH4J7(P24,P12,P23,P14,P34,P13,EP2,EP4,EP1,EP3,ORIENT)
2151
2152 & +HWH4J7(P13,P23,P34,P12,P14,P24,EP3,EP1,EP2,EP4,ORIENT)
2153
2154 & +HWH4J7(P24,P23,P12,P34,P14,P13,EP2,EP4,EP3,EP1,ORIENT)
2155
2156 & +HWH4J7(P13,P14,P12,P34,P23,P24,EP1,EP3,EP4,EP2,ORIENT)
2157
2158 & +HWH4J7(P24,P14,P34,P12,P23,P13,EP4,EP2,EP1,EP3,ORIENT)
2159
2160 & +HWH4J7(P13,P34,P23,P14,P12,P24,EP3,EP1,EP4,EP2,ORIENT)
2161
2162 & +HWH4J7(P24,P34,P14,P23,P12,P13,EP4,EP2,EP3,EP1,ORIENT)
2163
2164C Add up weights
2165
2166 WTQP=FACTR*COLC*QP/TWO
2167
2168 QQ1 =COLC*QP
2169
2170 QQ2 =COLC*QQ
2171
2172 QQINT=COLB*QQINT
2173
2174 WTQQ=FACTR*(QQ1+QQ2+QQINT)/2
2175
2176 ENDIF
2177
2178C
2179
2180 HCS=0.
2181
2182 DO 60 ID1=IDMN,IDMX
2183
2184 IF (INCLQG(ID1)) THEN
2185
2186C Gluon channel
2187
2188 HCS=HCS+CLF(1,ID1)*WTGG
2189
2190 IF (GENEV.AND.HCS.GT.RCS) THEN
2191
2192C Select colour flow
2193
2194 WTAB=GG1
2195
2196 WTBA=GG2
2197
2198 IF (IOP4JT(1).EQ.1) THEN
2199
2200 IF (GGINT.GE.ZERO) THEN
2201
2202 WTAB=WTAB+GGINT
2203
2204 ELSE
2205
2206 WTBA=MAX(WTBA,WTBA+GGINT)
2207
2208 ENDIF
2209
2210 ELSEIF (IOP4JT(1).EQ.2) THEN
2211
2212 IF (GGINT.GE.ZERO) THEN
2213
2214 WTBA=WTBA+GGINT
2215
2216 ELSE
2217
2218 WTAB=MAX(WTAB,WTAB+GGINT)
2219
2220 ENDIF
2221
2222 ELSEIF (IOP4JT(1).NE.0) THEN
2223
2224 CALL HWWARN('HWH4JT',101,*999)
2225
2226 ENDIF
2227
2228 WTOT=WTAB+WTBA
2229
2230 IF (WTAB.GT.HWR()*WTOT) THEN
2231
2232 CALL HWHQCP( 13, 13,3142,91,*99)
2233
2234 ELSE
2235
2236 CALL HWHQCP( 13, 13,4123,92,*99)
2237
2238 ENDIF
2239
2240 ENDIF
2241
2242 ENDIF
2243
2244C Quark channels
2245
2246 DO 50 ID2=1,6
2247
2248C Identical quark pairs
2249
2250 IF (ID1.EQ.ID2.AND.INCLQQ(ID1,ID1)) THEN
2251
2252 HCS=HCS+CLF(1,ID1)*WTQQ
2253
2254 IF (GENEV.AND.HCS.GT.RCS) THEN
2255
2256C Select colour flow
2257
2258 WTAB=QQ1
2259
2260 WTBA=QQ2
2261
2262 IF (IOP4JT(2).EQ.1) THEN
2263
2264 IF (QQINT.GE.ZERO) THEN
2265
2266 WTAB=WTAB+QQINT
2267
2268 ELSE
2269
2270 WTBA=MAX(WTBA,WTBA+QQINT)
2271
2272 ENDIF
2273
2274 ELSEIF (IOP4JT(2).EQ.2) THEN
2275
2276 IF (QQINT.GE.ZERO) THEN
2277
2278 WTBA=WTBA+QQINT
2279
2280 ELSE
2281
2282 WTAB=MAX(WTAB,WTAB+QQINT)
2283
2284 ENDIF
2285
2286 ELSEIF (IOP4JT(2).NE.0) THEN
2287
2288 CALL HWWARN('HWH4JT',102,*999)
2289
2290 ENDIF
2291
2292 WTOT=WTAB+WTBA
2293
2294 IF (WTAB.GT.HWR()*WTOT) THEN
2295
2296 CALL HWHQCP(ID1,ID1+6,4123,93,*99)
2297
2298 ELSE
2299
2300 CALL HWHQCP(ID1,ID1+6,2143,94,*99)
2301
2302 ENDIF
2303
2304 ENDIF
2305
2306C Unlike quark pairs
2307
2308 ELSEIF (INCLQQ(ID1,ID2)) THEN
2309
2310 HCS=HCS+(CLF(1,ID1)+CLF(1,ID2))*WTQP
2311
2312 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID2,ID2+6,4123,95,*99)
2313
2314 ENDIF
2315
2316 50 CONTINUE
2317
2318 60 CONTINUE
2319
2320 EVWGT=HCS
2321
2322 RETURN
2323
2324C Set up labels for selected final state
2325
2326 99 IDN(1)=ID1
2327
2328 IDN(2)=ID1+6
2329
2330 J=NHEP+1
2331
2332 IDHW(J)=200
2333
2334 IDHEP(J)=23
2335
2336 ISTHEP(J)=110
2337
2338 JMOHEP(1,J)=LM
2339
2340 JMOHEP(2,J)=LP
2341
2342 JDAHEP(1,J)=NHEP+2
2343
2344 JDAHEP(2,J)=NHEP+5
2345
2346 DO 100 I=1,4
2347
2348 J=NHEP+1+I
2349
2350 IDHW(J)=IDN(I)
2351
2352 IDHEP(J)=IDPDG(IDN(I))
2353
2354 ISTHEP(J)=IST(I)
2355
2356 JMOHEP(1,J)=NHEP+1
2357
2358 100 JDAHEP(1,J)=0
2359
2360C And colour structure pointers
2361
2362 DO 110 I=1,4
2363
2364 J=ICO(I)
2365
2366 JMOHEP(2,NHEP+1+I)=NHEP+1+J
2367
2368 110 JDAHEP(2,NHEP+1+J)=NHEP+1+I
2369
2370 NHEP=NHEP+5
2371
2372 999 END
2373
2374CDECK ID>, HWH4J1.
2375
2376*CMZ :- -01/04/99 19.47.55 by Mike Seymour
2377
2378*-- Author : Ian Knowles
2379
2380C-----------------------------------------------------------------------
2381
2382 FUNCTION HWH4J1(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2383
2384C-----------------------------------------------------------------------
2385
2386C Evaluate `ERT' functions A, B, C, D, E; S12=(p1+p2)^2 etc.
2387
2388C-----------------------------------------------------------------------
2389
2390 IMPLICIT NONE
2391
2392 DOUBLE PRECISION HWH4J1,HWH4J2,HWH4J4,HWH4J5,HWH4J6,HWH4J7,
2393
2394 & S12,S13,S14,S23,S24,S34,S123,S124,S134,S234,S,EP1,EP2,EP3,EP4,
2395
2396 & SUM
2397
2398 LOGICAL ORIENT
2399
2400 S123=S12+S13+S23
2401
2402 S124=S12+S14+S24
2403
2404 S134=S13+S14+S34
2405
2406 S234=S23+S24+S34
2407
2408 S=S12+S13+S14+S23+S24+S34
2409
2410 HWH4J1=(S12*((S12+S14+S23+S34)**2+S13*(S12+S14-S24)+S24*(S12+S23))
2411
2412 & +(S14*S23-S12*S34-S13*S24)*(S14+S23+S34)/2)
2413
2414 & /(S13*S24*S134*S234)
2415
2416 & +((S12+S24)*(S13+S34)-S14*S23)/(S13*S134**2)
2417
2418 & +2*S23*(S-S13)/(S13*S134*S24) + S34/(2*S13*S24)
2419
2420 IF (ORIENT) THEN
2421
2422 HWH4J1=HWH4J1
2423
2424 & +4*((EP1*EP1*((S-S13)*(S23+S24)-S24*S34)
2425
2426 & -EP1*EP2*(S12*(S123+S124)+(S+S12)*(S14+S23)+2*S14*S23
2427
2428 & +S24*S134+S234*(S13+2*S234))
2429
2430 & +EP1*EP3*(S*(S24-S12)+S12*S13+(S14+2*S234-S34)*S24)
2431
2432 & -EP1*EP4*(S12*S124+S23*(S+S12+S14))
2433
2434 & +EP2*EP2*((S-S24)*(S13+S14)+2*(S13+S34)*S234-S13*S34)
2435
2436 & -EP2*EP3*((S+S23)*(S12+S14)+(S12+2*(S23+S234))*S234)
2437
2438 & +EP2*EP4*(S12*(S24-S)+S13*(S+S23-S34)+2*(S13+S34-S234)*S234)
2439
2440 & +EP3*EP3*(S14+2*S234)*S24
2441
2442 & +EP3*EP4*(-S234*(2*(S12+S23)+S134)+S12*S34-S13*S24-S14*S23)
2443
2444 & +EP4*EP4*S13*S23)*S134
2445
2446 & +EP2*(EP1+EP3+EP4)*2*S14*S24*S234)/(S*S13*S24*S134**2*S234)
2447
2448 ELSE
2449
2450 HWH4J1=2*HWH4J1/3
2451
2452 ENDIF
2453
2454 RETURN
2455
2456C-----------------------------------------------------------------------
2457
2458 ENTRY HWH4J2(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2459
2460C-----------------------------------------------------------------------
2461
2462 S123=S12+S13+S23
2463
2464 S124=S12+S14+S24
2465
2466 S134=S13+S14+S34
2467
2468 S234=S23+S24+S34
2469
2470 S=S12+S13+S14+S23+S24+S34
2471
2472 HWH4J2=(S12*S14*(S24+S34)+S24*(S12*(S14+S34)+S13*(S14-S24)))
2473
2474 & /(S14*S23*S13*S134)
2475
2476 & +S12*(S+S34)*S124/(S24*S234*S14*S134)
2477
2478 & -(S13*(2*(S12+S24)+S23)+S14**2)/(S134*S13*S14)
2479
2480 & +S12*S123*S124/(2*S13*S24*S14*S23)
2481
2482 IF (ORIENT) THEN
2483
2484 HWH4J2=HWH4J2
2485
2486 & +4*((EP1*EP1*(S12*S134*S234-4*S23*S24*S34)
2487
2488 & +EP1*EP2*(2*(2*S13*S234+S14*S123)*S24-S12*S134*(S+S12+S34))
2489
2490 & +EP1*EP3*(S12*(4*S24*S34-S134*(S12+S14-S24))
2491
2492 & -4*(S13*S24-S14*S23)*S24)
2493
2494 & +EP1*EP4*(4*(S13+S14)*S23*S24-S12*S134*(S12+S13-S23))
2495
2496 & +EP2*EP2*(S12*S134-4*S13*S24)*S134
2497
2498 & +EP2*EP3*(4*S13*(S12+S23+S24)*S24-S12*S134*(S12-S14+S24))
2499
2500 & -EP2*EP4*(4*(S12*(S14+S134)+S13*(S134-S234))*S24
2501
2502 & +S12*(S12-S13+S23)*S134)
2503
2504 & -EP3*EP3*4*S12*S14*S24
2505
2506 & -EP3*EP4*2*S12*(2*S14*S24+S12*S134))*S234
2507
2508 & +(EP1*(EP1*(S23+S24)+EP2*(S134-2*S))
2509
2510 & -(EP1+EP2)*(EP3+EP4)*S12+EP2*EP2*(S13+S14))*2*S14*S24*S123)
2511
2512 & /(2*S*S13*S14*S234*S23*S24*S134)
2513
2514 ELSE
2515
2516 HWH4J2=2*HWH4J2/3
2517
2518 ENDIF
2519
2520 RETURN
2521
2522C-----------------------------------------------------------------------
2523
2524 ENTRY HWH4J4(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2525
2526C-----------------------------------------------------------------------
2527
2528 S134=S13+S14+S34
2529
2530 S234=S23+S24+S34
2531
2532 S=S12+S13+S14+S23+S24+S34
2533
2534 HWH4J4=-(S12*(S34*(3*(S+S34)+S12)-S134*S234-2*(S13*S24+S14*S23))
2535
2536 & +(S14*S23-S13*S24)*(S13-S14+S24-S23))/(2*S134*S234*S34**2)
2537
2538 & -(S12*(S134**2/2+2*S13*S14+S34*(S13+S14-S34))
2539
2540 & +S34*((S13+S14)*(S23+S24)+S14*S24+S13*S23)
2541
2542 & +(S13*S24-S14*S23)*(S14-S13))/(S34*S134)**2
2543
2544 IF (ORIENT) THEN
2545
2546 HWH4J4=HWH4J4
2547
2548 & +4*((-EP1*EP1*2*(S23+S24)*S34
2549
2550 & -EP1*EP2*(S13*(S23+3*S24)+S14*(3*S23+S24)-(4*S12-S34)*S34)
2551
2552 & +EP1*EP3*((2*S12-S24)*S34-(S13-S14)*S24)
2553
2554 & +EP1*EP4*((2*S12-S23)*S34+(S13-S14)*S23)
2555
2556 & -EP2*EP2*2*(S13+S14)*S34
2557
2558 & +EP2*EP3*(2*S12*S34-S14*(S23-S24+S34))
2559
2560 & +EP2*EP4*(2*S12*S34+S13*(S23-S24-S34))
2561
2562 & +EP3*EP3*2*S14*S24
2563
2564 & +EP3*EP4*2*(S12*S34-S13*S24-S14*S23)
2565
2566 & +EP4*EP4*2*S13*S23)/(S*S134*S234*S34**2)
2567
2568 & +(EP1*EP2*(S134*(S134+2*S34)+4*(S13*S14-S34**2))
2569
2570 & +EP2*EP3*2*(2*S13*S34+S14*(S13-S14+S34))
2571
2572 & +EP2*EP4*2*(2*S14*S34-S13*(S13-S14-S34)))
2573
2574 & /(S*(S134*S34)**2))
2575
2576 ELSE
2577
2578 HWH4J4=2*HWH4J4/3
2579
2580 ENDIF
2581
2582 RETURN
2583
2584C-----------------------------------------------------------------------
2585
2586 ENTRY HWH4J5(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2587
2588C-----------------------------------------------------------------------
2589
2590 S123=S12+S13+S23
2591
2592 S124=S12+S14+S24
2593
2594 S134=S13+S14+S34
2595
2596 S234=S23+S24+S34
2597
2598 S=S12+S13+S14+S23+S24+S34
2599
2600 HWH4J5=(3*S12*S34**2-3*S13*S24*S34+3*S12*S24*S34+3*S14*S23*S34-
2601
2602 $ S13*S24**2-S12*S23*S34+6*S12*S14*S34+2*S12*S13*S34-
2603
2604 $ 2*S12**2*S34+S14*S23*S24-3*S13*S23*S24-2*S13*S14*S24+
2605
2606 $ 4*S12*S14*S24+2*S12*S13*S24+3*S14*S23**2+2*S14**2*S23+
2607
2608 $ 2*S14**2*S12+2*S12**2*S14+6*S12*S14*S23-2*S12*S13**2-
2609
2610 $ 2*S12**2*S13)/(2*S13*S134*S234*S34)+
2611
2612 $ (2*S12*S34**2-2*S13*S24*S34+S12*S24*S34+4*S13*S23*S34+
2613
2614 $ 4*S12*S14*S34+2*S12*S13*S34+2*S12**2*S34-S13*S24**2+
2615
2616 $ 3*S14*S23*S24+4*S13*S23*S24-2*S13*S14*S24+4*S12*S14*S24+
2617
2618 $ 2*S12*S13*S24+2*S14*S23**2+4*S13*S23**2+2*S13*S14*S23+
2619
2620 $ 2*S12*S14*S23+4*S12*S13*S23+2*S12*S14**2+4*S12**2*S13+
2621
2622 $ 4*S12*S13*S14+2*S12**2*S14)/(2*S13*S134*S24*S34)-
2623
2624 $ (S12*S34**2-2*S14*S24*S34-2*S13*S24*S34-S14*S23*S34+
2625
2626 $ S13*S23*S34+S12*S14*S34+2*S12*S13*S34-2*S14**2*S24-
2627
2628 $ 4*S13*S14*S24-4*S13**2*S24-S14**2*S23-S13**2*S23+
2629
2630 $ S12*S13*S14-S12*S13**2)/(S13*S34*S134**2)
2631
2632 IF (ORIENT) THEN
2633
2634 SUM=
2635
2636 & +EP1*EP1*((S13-S14+S23-3*S24)*S34+(S134+S14+2*S34)*S234)
2637
2638 & *S24*S134
2639
2640 & +EP1*EP2*((2*(S12-S24)+S34)*S134-S14*(4*S12+S14+3*S23)
2641
2642 & +S13*(S13+S23)+S24*S34 )*S24*S134
2643
2644 & -EP1*EP2*(((2*S12*S134+S13*(2*(S12+S14+S23)-S24+S34)
2645
2646 & +S14*(S14-S23)+(2*S14-S34)*S234)*S234)*S134
2647
2648 & + 4*S13**2*S24*S234)
2649
2650 & +EP1*EP3*(S12*(2*S13-S134)+S13*(S24+2*S234)+S14*(3*S24-S234)
2651
2652 & +S34*(S234-3*S24))*S24*S134
2653
2654 & +EP1*EP4*((S12*(S13-S14+3*S34)-S23*(S13+3*S14-S34))*S24
2655
2656 & -(S12*(S13+S134+2*S34)+2*S13*S24
2657
2658 & +(S13-2*S14)*S23)*S234)*S134
2659
2660 & +EP2*EP2*(S13*((2*S13+S34)*S234+S24*(S134-2*S34))
2661
2662 & +2*S14*S134*(S24+S234))*S134
2663
2664 SUM=SUM
2665
2666 & -EP2*EP3*(((S12*(S13+2*S14-S34)+S14*(S+2*S23-S34))*S24
2667
2668 & +(S12*(S13+S134)+(S13+S24+2*S234)*S14
2669
2670 & +2*S13*(2*S23+S34))*S234)*S134
2671
2672 & +4*S13**2*S24*S234)
2673
2674 & +EP2*EP4*(((S12*(S13-2*S134)+S13*(S+2*S23-3*S34))*S24
2675
2676 & -((S-3*S13+S23+2*S24)*S13+2*S12*S14
2677
2678 & +2*S14*(S23+2*S24))*S234)*S134-4*S13**2*S24*S234)
2679
2680 & +EP3*EP3*2*(S13*S234+S14*S24)*S24*S134
2681
2682 & +EP3*EP4*(2*(S12*S34-S13*S24-S14*S23)*S24
2683
2684 & -(S12*S134+2*S13*S23)*S234)*S134
2685
2686 & +EP4*EP4*2*(S12*S234+S23*S24)*S13*S134
2687
2688 HWH4J5=HWH4J5+4*SUM/(S*S234*S134**2*S13*S34*S24)
2689
2690 ELSE
2691
2692 HWH4J5=2*HWH4J5/3
2693
2694 ENDIF
2695
2696 RETURN
2697
2698C-----------------------------------------------------------------------
2699
2700 ENTRY HWH4J6(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2701
2702C-----------------------------------------------------------------------
2703
2704 S123=S12+S13+S23
2705
2706 S124=S12+S14+S24
2707
2708 S134=S13+S14+S34
2709
2710 S234=S23+S24+S34
2711
2712 S=S12+S13+S14+S23+S24+S34
2713
2714 HWH4J6=(S23*(S123*S234-S*S23)+S12*(S123*S124-S*S12))/(S13*S123)**2
2715
2716 & -(S12*S34*(S234-2*S23)+S14*S23*(S234-2*S34)
2717
2718 & -S13*S24*(S234+S13))/(S13**2*S123*S134)
2719
2720 IF (ORIENT) THEN
2721
2722 HWH4J6=HWH4J6
2723
2724 & +4*(-EP1*EP1*2*S23*S34
2725
2726 & +EP1*EP2*((S12-S23)*S34-S13*(S24-S34))
2727
2728 & +(EP1*EP3+EP2*EP4)*2*(S12*S34-S13*S24+S14*S23)
2729
2730 & -EP1*EP4*(S13*S24-(3*(S13+S14)+S34)*S23)
2731
2732 & -(EP1+EP2+EP3)*EP4*2
2733
2734 & *(S12*(S13+S23)+(S12+S13)*S23)*S134/S123
2735
2736 & +EP2*EP2*S13*(S14+S34)
2737
2738 & +EP2*EP3*(S13*(S14-S24)-(S12-S23)*S14)
2739
2740 & -EP3*EP3*2*S12*S14
2741
2742 & -EP3*EP4*(S13*S24-(3*(S13+S34)+S14)*S12)
2743
2744 & +EP4*EP4*(S12+S23)*S13)/(S*S134*S123*S13**2)
2745
2746 ELSE
2747
2748 HWH4J6=2*HWH4J6/3
2749
2750 ENDIF
2751
2752 RETURN
2753
2754C-----------------------------------------------------------------------
2755
2756 ENTRY HWH4J7(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2757
2758C-----------------------------------------------------------------------
2759
2760 S123=S12+S13+S23
2761
2762 S124=S12+S14+S24
2763
2764 S134=S13+S14+S34
2765
2766 S234=S23+S24+S34
2767
2768 S=S12+S13+S14+S23+S24+S34
2769
2770 HWH4J7=((S12*S34+S13*S24-S14*S23)*(S13+S14+S23+S24)-2*S12*S24*S34)
2771
2772 & /(S13*S134*S23*S123)
2773
2774 & -S12*(S12*S-S123*S124)/(S123**2*S13*S23)
2775
2776 & -(S13+S14)*(S23+S24)*S34/(S13*S134*S23*S234)
2777
2778 IF (ORIENT) THEN
2779
2780 HWH4J7=HWH4J7
2781
2782 & +4*(+2*(EP1+EP2)*(S23*EP1-S13*EP2)*S34*S134
2783
2784 & -EP1*EP2*2*S34**2*S123
2785
2786 & +EP1*EP3*(S123*(S23+S24)*S34+2*S134*(S13*S24-S14*S23))
2787
2788 & +EP1*EP4*(S123*(S23+S24)*S34+2*S12**2*S134*S234/S123
2789
2790 & +2*S134*(S24*(S13-S12)-S23*(S12+S14)))
2791
2792 & +EP2*EP3*(2*(S12*S34+S13*S24-S14*S23)*S134
2793
2794 & +S123*(S13+S14)*S34)
2795
2796 & +EP2*EP4*(S123*(S13+S14)*S34+2*S12**2*S234*S134/S123
2797
2798 & -2*S134*(S12*S234-S13*S24+S14*S23))
2799
2800 & -EP3*EP3*S12*(2*S24*S134+S123*S34)
2801
2802 & +EP3*EP4*2*S12*(S134*(S23-S24)-S34*S123+S12*S134*S234/S123)
2803
2804 & +EP4*EP4*S12*(2*S23*S134-S123*S34))
2805
2806 & /(S*S13*S23*S123*S134*S234)
2807
2808 ELSE
2809
2810 HWH4J7=2*HWH4J7/3
2811
2812 ENDIF
2813
2814 RETURN
2815
2816 END