]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HERWIG/src/hwissp.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Additional protection in Digitize, which was moved to the implementation file
[u/mrichter/AliRoot.git] / HERWIG / src / hwissp.f
CommitLineData
3820ca8e 1
2CDECK ID>, HWISSP.
3
4*CMZ :- -20/10/99 09:46:43 by Peter Richardson
5
6*-- Author : Bryan Webber, modified by Kosuke Odagiri
7
8C-----------------------------------------------------------------------
9
10 SUBROUTINE HWISSP(LRSUSY)
11
12C-----------------------------------------------------------------------
13
14C Reads in SUSY particle properties and decays,
15
16C in format generated by ISAWIG
17
18C-----------------------------------------------------------------------
19
20 INCLUDE 'HERWIG61.INC'
21
22 INTEGER I,J,K,IH,IHW,NSSP,NDEC,MDKYS,LRSUSY
23
24 DOUBLE PRECISION BETAH, WEINCOS,WEINSIN, MW,MZ, RMMAX
25
26 DOUBLE PRECISION FTM,FTMUU(4),FTMDD(4),FTMTT(4),FTMBB(4),FTMU,FTMD
27
28 DOUBLE PRECISION YTM,YTM1,DTERM(4), SQHF,SNBCSB,MZSW2
29
30 LOGICAL FIRST
31
32 EQUIVALENCE (MW,RMASS(198)), (MZ,RMASS(200))
33
34 DATA FIRST/.TRUE./
35
36 SAVE MDKYS
37
38 IF (FIRST) THEN
39
40 MDKYS=NDKYS
41
42 FIRST=.FALSE.
43
44 ELSE
45
46 NDKYS=MDKYS
47
48 ENDIF
49
50C--reset susy input flag
51
52 IF (LRSUSY.LT.0) CALL HWWARN('HWISSP',500,*999)
53
54 SUSYIN = .TRUE.
55
56C
57
58C Input SUSY particle + top quark table
59
60C
61
62 WRITE (6,10) LRSUSY
63
64 10 FORMAT (10X,'Reading in SUSY data from unit',I3)
65
66 READ (LRSUSY,'(I4)') NSSP
67
68 IF (NSSP.LE.0) RETURN
69
70 RMMAX=SQRT(HALF*(EBEAM1*EBEAM2+PBEAM1*PBEAM2))
71
72 RMMNSS=RMMAX
73
74 DO I=1,NSSP
75
76 READ (LRSUSY,1) IHW,RMASS(IHW),RLTIM(IHW)
77
78C Negative gaugino mass means physical field is gamma_5*psi
79
80C Store the signs
81
82 IF ((IHW.GE.450).AND.(IHW.LE.457)) THEN
83
84 IF (IHW.LE.453) THEN
85
86 J=IHW-449
87
88 ZSGNSS(J)=RMASS(IHW)/ABS(RMASS(IHW))
89
90 ELSEIF (IHW.LE.455) THEN
91
92 J=IHW-453
93
94 WSGNSS(J)=RMASS(IHW)/ABS(RMASS(IHW))
95
96 ENDIF
97
98 RMASS(IHW)=ABS(RMASS(IHW))
99
100 ENDIF
101
102 IF (ABS(IDPDG(IHW)).GT.1000000.AND.(RMASS(IHW).NE.ZERO))
103
104 & RMMNSS=MIN(RMMNSS,RMASS(IHW))
105
106 IF (IHW.GT.NRES) THEN
107
108 IF (IHW.GT.NMXRES) CALL HWWARN('HWISSP',501,*999)
109
110 NRES=IHW
111
112 ENDIF
113
114 ENDDO
115
116 XLMNSS=TWO*LOG(RMMNSS/RMMAX)
117
118 1 FORMAT(I5,F12.4,E15.5)
119
120C
121
122C Input decay modes
123
124C
125
126 DO I=1,NSSP
127
128 READ (LRSUSY,'(I4)') NDEC
129
130 IF (NDEC.GT.0) THEN
131
132 DO J=1,NDEC
133
134 NDKYS=NDKYS+1
135
136 IF (NDKYS.GT.NMXDKS) CALL HWWARN('HWISSP',100,*999)
137
138 READ (LRSUSY,11) IDK(NDKYS),BRFRAC(NDKYS),NME(NDKYS),
139
140 & (IDKPRD(K,NDKYS),K=1,5)
141
142 11 FORMAT(I6,F16.8,6I6)
143
144 ENDDO
145
146 ENDIF
147
148 ENDDO
149
150C
151
152C Mixings and other SUSY parameters
153
154C
155
156 READ (LRSUSY,'(2F16.8)') TANB,ALPHAH
157
158 DO I=1,4
159
160 READ (LRSUSY,13) ZMXNSS(I,1),ZMXNSS(I,2),ZMXNSS(I,3),ZMXNSS(I,4)
161
162 END DO
163
164 WEINSIN = SQRT(SWEIN)
165
166 WEINCOS = SQRT(1.-SWEIN)
167
168 DO I=1,4
169
170 ZMIXSS(I,1) = WEINCOS*ZMXNSS(I,1)+WEINSIN*ZMXNSS(I,2)
171
172 ZMIXSS(I,2) = -WEINSIN*ZMXNSS(I,1)+WEINCOS*ZMXNSS(I,2)
173
174 ZMIXSS(I,3) = ZMXNSS(I,3)
175
176 ZMIXSS(I,4) = ZMXNSS(I,4)
177
178 END DO
179
180 DO J=1,16
181
182 IF ((J.LE.6).OR.(J.GE.11)) THEN
183
184 LFCH(J)=VFCH(J,1)+AFCH(J,1)
185
186 RFCH(J)=VFCH(J,1)-AFCH(J,1)
187
188 DO I=1,4
189
190 SLFCH(J,I)= ZMIXSS(I,1)*QFCH(J)+ZMIXSS(I,2)*LFCH(J)
191
192 SRFCH(J,I)=-ZMIXSS(I,1)*QFCH(J)-ZMIXSS(I,2)*RFCH(J)
193
194 END DO
195
196 ENDIF
197
198 END DO
199
200 READ (LRSUSY,13) WMXVSS(1,1),WMXVSS(1,2), WMXVSS(2,1),WMXVSS(2,2)
201
202 READ (LRSUSY,13) WMXUSS(1,1),WMXUSS(1,2), WMXUSS(2,1),WMXUSS(2,2)
203
204 READ (LRSUSY,'(3F16.8)') THETAT,THETAB,THETAL
205
206 READ (LRSUSY,'(3F16.8)') ATSS,ABSS,ALSS
207
208 READ (LRSUSY,'( F16.8)') MUSS
209
210 DO I=1,6
211
212 QMIXSS(I,1,1)=1.
213
214 QMIXSS(I,1,2)=0.
215
216 QMIXSS(I,2,1)=0.
217
218 QMIXSS(I,2,2)=1.
219
220 LMIXSS(I,1,1)=1.
221
222 LMIXSS(I,1,2)=0.
223
224 LMIXSS(I,2,1)=0.
225
226 LMIXSS(I,2,2)=1.
227
228 END DO
229
230 QMIXSS(6,1,1)= COS(THETAT)
231
232 QMIXSS(6,1,2)= SIN(THETAT)
233
234 QMIXSS(6,2,1)=-QMIXSS(6,1,2)
235
236 QMIXSS(6,2,2)= QMIXSS(6,1,1)
237
238 QMIXSS(5,1,1)= COS(THETAB)
239
240 QMIXSS(5,1,2)= SIN(THETAB)
241
242 QMIXSS(5,2,1)=-QMIXSS(5,1,2)
243
244 QMIXSS(5,2,2)= QMIXSS(5,1,1)
245
246 LMIXSS(5,1,1)= COS(THETAL)
247
248 LMIXSS(5,1,2)= SIN(THETAL)
249
250 LMIXSS(5,2,1)=-LMIXSS(5,1,2)
251
252 LMIXSS(5,2,2)= LMIXSS(5,1,1)
253
254C--Evaluating Higgs parameters and couplings
255
256 BETAH=ATAN(TANB)
257
258 COTB=ONE/TANB
259
260 COSBPA=COS(BETAH+ALPHAH)
261
262 SINBPA=SIN(BETAH+ALPHAH)
263
264 COSBMA=COS(BETAH-ALPHAH)
265
266 SINBMA=SIN(BETAH-ALPHAH)
267
268 COSA=COS(ALPHAH)
269
270 SINA=SIN(ALPHAH)
271
272 COSB=COS(BETAH)
273
274 SINB=SIN(BETAH)
275
276 GHWWSS(1)=SINBMA
277
278 GHWWSS(2)=COSBMA
279
280 GHWWSS(3)=ZERO
281
282 DO 30 I=1,3
283
284 GHZZSS(I)=GHWWSS(I)
285
286 30 CONTINUE
287
288 GHDDSS(1)=-SINA/COSB
289
290 GHDDSS(2)= COSA/COSB
291
292 GHDDSS(3)= TANB
293
294 GHUUSS(1)= COSA/SINB
295
296 GHUUSS(2)= SINA/SINB
297
298 GHUUSS(3)= COTB
299
300 GHWHSS(1)= COSBMA
301
302 GHWHSS(2)= SINBMA
303
304 GHWHSS(3)= ONE
305
306 MZSW2 = MZ**2 * SQRT(SWEIN*(ONE-SWEIN))
307
308 DTERM(1) =-SINBPA*MZSW2
309
310 DTERM(2) = COSBPA*MZSW2
311
312 DTERM(3) = ZERO
313
314 FTMUU(1) =-MUSS*SINA/SINB
315
316 FTMUU(2) = MUSS*COSA/SINB
317
318 FTMUU(3) = MUSS
319
320 FTMUU(4) = MUSS
321
322 FTMTT(1) =-ATSS*COSA/SINB
323
324 FTMTT(2) =-ATSS*SINA/SINB
325
326 FTMTT(3) = ATSS*COTB
327
328 FTMTT(4) = ATSS*COTB
329
330 FTMDD(1) = MUSS*COSA/COSB
331
332 FTMDD(2) = MUSS*SINA/COSB
333
334 FTMDD(3) = MUSS
335
336 FTMDD(4) = MUSS
337
338 FTMBB(1) = ABSS*SINA/COSB
339
340 FTMBB(2) =-ABSS*COSA/COSB
341
342 FTMBB(3) = ABSS*TANB
343
344 FTMBB(4) = ABSS*TANB
345
346 DO 40 IH=1,4
347
348 FTMU=FTMUU(IH)
349
350 FTMD=FTMDD(IH)
351
352 DO 50 I=1,6
353
354 IF (I.EQ.5) FTMU=FTMU+FTMTT(IH)
355
356 IF (I.EQ.5) FTMD=FTMD+FTMBB(IH)
357
358 IF (MOD(I,2).EQ.0) THEN
359
360 YTM = GHUUSS(IH)
361
362 FTM = FTMU
363
364 ELSE
365
366 YTM = GHDDSS(IH)
367
368 FTM = FTMD
369
370 END IF
371
372 IF (IH.EQ.3) THEN
373
374 GHSQSS(IH,I,1,1) = ZERO
375
376 GHSQSS(IH,I,2,2) = ZERO
377
378 GHSQSS(IH,I,1,2) = FTM*HALF*RMASS(I)/MW
379
380 GHSQSS(IH,I,2,1) = - GHSQSS(IH,I,1,2)
381
382 GOTO 50
383
384 ELSEIF (IH.EQ.4) THEN
385
386 SQHF=SQRT(HALF)
387
388 SNBCSB=SINB*COSB
389
390 DO 60 J=1,2
391
392 DO 70 K=1,2
393
394 IF (MOD(I,2).EQ.1) THEN
395
396 GHSQSS(IH,I,J,K)=SQHF*(
397
398 & RMASS(I )*FTMD*QMIXSS(I,2,J)*QMIXSS(I+1,1,K)
399
400 & +RMASS(I+1)*FTMU*QMIXSS(I,1,J)*QMIXSS(I+1,2,K)
401
402 & +( MW**2*TWO*SNBCSB-RMASS(I+1)**2*COTB
403
404 & -RMASS(I )**2*TANB )*QMIXSS(I,1,J)*QMIXSS(I+1,1,K)
405
406 & -RMASS(I)*RMASS(I+1)/SNBCSB
407
408 & *QMIXSS(I,2,J)*QMIXSS(I+1,2,K) ) / MW
409
410 ELSE
411
412 GHSQSS(IH,I,J,K)=GHSQSS(IH,I-1,K,J)
413
414 END IF
415
416 70 END DO
417
418 60 END DO
419
420 ELSE
421
422 DO 80 J=1,2
423
424 DO 90 K=1,2
425
426 YTM1=ZERO
427
428 IF (J.EQ.K) YTM1=YTM*RMASS(I)**2
429
430 GHSQSS(IH,I,J,K)=( YTM1
431
432 & +( LFCH(I)*QMIXSS(I,1,J)*QMIXSS(I,1,K)
433
434 & -RFCH(I)*QMIXSS(I,2,J)*QMIXSS(I,2,K) )*DTERM(IH)
435
436 & +FTM*HALF*RMASS(I)*(QMIXSS(I,1,J)*QMIXSS(I,2,K)
437
438 & +QMIXSS(I,2,J)*QMIXSS(I,1,K)) ) / MW
439
440 90 CONTINUE
441
442 80 CONTINUE
443
444 END IF
445
446 50 CONTINUE
447
448 40 CONTINUE
449
450C--Rparity violation
451
452 READ (LRSUSY,'(L5)') RPARTY
453
454 IF(.NOT.RPARTY) THEN
455
456 READ(LRSUSY,20) (((LAMDA1(I,J,K),K=1,3),J=1,3),I=1,3)
457
458 READ(LRSUSY,20) (((LAMDA2(I,J,K),K=1,3),J=1,3),I=1,3)
459
460 READ(LRSUSY,20) (((LAMDA3(I,J,K),K=1,3),J=1,3),I=1,3)
461
462 ENDIF
463
464 13 FORMAT(4F16.8)
465
466 20 FORMAT(27E16.8)
467
468 CLOSE(LRSUSY)
469
470 999 END
471
472CDECK ID>, HWMEVT.
473
474*CMZ :- -04/05/99 14.28.59 by Bryan Webber
475
476*-- Author : Bryan Webber
477
478C-----------------------------------------------------------------------
479
480 SUBROUTINE HWMEVT
481
482C-----------------------------------------------------------------------
483
484C IPROC = 1000,... ADDS SOFT UNDERLYING EVENT
485
486C = 8000: CREATES MINIMUM-BIAS EVENT
487
488C SUPPRESSED BY ADDING 10000 TO IPROC
489
490C-----------------------------------------------------------------------
491
492 INCLUDE 'HERWIG61.INC'
493
494 DOUBLE PRECISION HWREXP,ENFAC,TECM,SECM,SUMM,EMCL,BMP(5),BMR(3,3)
495
496 INTEGER HWRINT,NETC,IBT,IDBT,ID1,ID2,ID3,KHEP,LHEP,NTRY,ICMS,
497
498 & NPPBAR,MCHT,JCL,JD1,JD2,JD3,ICH,MODC,NCHT,INHEP(2),
499
500 & INID(2,2),JBT
501
502 EXTERNAL HWREXP,HWRINT
503
504 IF (IERROR.NE.0) RETURN
505
506 IF (.NOT.GENSOF) GOTO 990
507
508 IF (IPROC.EQ.8000) THEN
509
510C---SET UP BEAM AND TARGET CLUSTERS
511
512 5 NETC=0
513
514 DO 10 IBT=1,2
515
516 JBT=IBT
517
518 IF (JDAHEP(1,IBT).NE.0) JBT=JDAHEP(1,IBT)
519
520 IDBT=IDHW(JBT)
521
522 IF (IDBT.EQ.73.OR.IDBT.EQ.75) THEN
523
524 INID(1,IBT)=HWRINT(1,2)
525
526 INID(2,IBT)=110
527
528 ELSEIF (IDBT.EQ.91.OR.IDBT.EQ.93) THEN
529
530 INID(1,IBT)=116
531
532 INID(2,IBT)=HWRINT(7,8)
533
534 ELSEIF (IDBT.EQ.30) THEN
535
536 INID(1,IBT)=HWRINT(1,2)
537
538 INID(2,IBT)=8
539
540 ELSEIF (IDBT.EQ.38) THEN
541
542 INID(1,IBT)=2
543
544 INID(2,IBT)=HWRINT(7,8)
545
546 ELSEIF (IDBT.EQ.34) THEN
547
548 INID(1,IBT)=3
549
550 INID(2,IBT)=HWRINT(7,8)
551
552 ELSEIF (IDBT.EQ.46) THEN
553
554 INID(1,IBT)=HWRINT(1,2)
555
556 INID(2,IBT)=9
557
558 ELSEIF (IDBT.EQ.59) THEN
559
560 INID(1,IBT)=HWRINT(1,2)
561
562 INID(2,IBT)=HWRINT(7,8)
563
564 ELSE
565
566 CALL HWWARN('HWMEVT',100,*999)
567
568 ENDIF
569
570 NETC=NETC+ICHRG(IDBT)
571
572 & -(ICHRG(INID(1,IBT))+ICHRG(INID(2,IBT)))/3
573
574 ENFAC=1.
575
576 IDHW(NHEP+IBT)=19
577
578 IDHEP(NHEP+IBT)=91
579
580 ISTHEP(NHEP+IBT)=163+IBT
581
582 JMOHEP(1,NHEP+IBT)=JBT
583
584 10 CONTINUE
585
586 IF (NETC.EQ.0) THEN
587
588 ID3=HWRINT(1,2)
589
590 ELSEIF (NETC.EQ.-1) THEN
591
592 ID3=1
593
594 ELSEIF (NETC.EQ.1) THEN
595
596 ID3=2
597
598 ELSE
599
600 GOTO 5
601
602 ENDIF
603
604 DO 12 IBT=1,2
605
606 NHEP=NHEP+1
607
608 JBT=IBT
609
610 IF (JDAHEP(1,IBT).NE.0) JBT=JDAHEP(1,IBT)
611
612 CALL HWVEQU(5,PHEP(1,JBT),PHEP(1,NHEP))
613
614 12 INHEP(IBT)=NHEP
615
616 ELSE
617
618C---FIND BEAM AND TARGET CLUSTERS
619
620 DO 20 IBT=1,2
621
622 DO 15 KHEP=1,NHEP
623
624 IF (ISTHEP(KHEP).EQ.163+IBT) THEN
625
626 INHEP(IBT)=KHEP
627
628 INID(1,IBT)=IDHW(JMOHEP(1,KHEP))
629
630 INID(2,IBT)=IDHW(JMOHEP(2,KHEP))
631
632 GOTO 20
633
634 ENDIF
635
636 15 CONTINUE
637
638C---COULDN'T FIND ONE
639
640 INHEP(IBT)=0
641
642 20 CONTINUE
643
644 JCL=-1
645
646C---TEST FOR BOTH FOUND
647
648 IF (INHEP(1).EQ.0) JCL=INHEP(2)
649
650 IF (INHEP(2).EQ.0) JCL=INHEP(1)
651
652 IF (JCL.EQ.0) CALL HWWARN('HWMEVT',101,*999)
653
654 IF (JCL.GT.0) THEN
655
656 ISTHEP(JCL)=163
657
658 CALL HWCFOR
659
660 CALL HWCDEC
661
662 CALL HWDHAD
663
664 CALL HWDHVY
665
666 GOTO 90
667
668 ENDIF
669
670 ID3=HWRINT(1,2)
671
672 ENFAC=ENSOF
673
674 NETC=0
675
676 ENDIF
677
678C---FIND SOFT CM MOMENTUM AND MULTIPLICITY
679
680 NTRY=0
681
682 NHEP=NHEP+1
683
684 IF (NHEP.GT.NMXHEP) CALL HWWARN('HWMEVT',102,*999)
685
686 ICMS=NHEP
687
688 IDHW(NHEP)=16
689
690 IDHEP(NHEP)=0
691
692 ISTHEP(NHEP)=170
693
694 CALL HWVSUM(4,PHEP(1,INHEP(1)),PHEP(1,INHEP(2)),PHEP(1,NHEP))
695
696 CALL HWUMAS(PHEP(1,NHEP))
697
698 TECM=PHEP(5,NHEP)
699
700 IF (IPROC/1000.EQ.9.OR.IPROC/1000.EQ.5) THEN
701
702 SECM=TECM*ENFAC
703
704 ELSE
705
706 SECM=PHEP(5,3)*ENFAC
707
708 ENDIF
709
710C---CHOOSE MULTIPLICITY
711
712 25 CALL HWMULT(SECM,NPPBAR)
713
714 30 NCL=0
715
716 MCHT=0
717
718 IERROR=0
719
720 NHEP =ICMS
721
722 SUMM=0.
723
724 NTRY=NTRY+1
725
726C---CREATE CLUSTERS
727
728 35 NCL=NCL+1
729
730 NHEP=NHEP+1
731
732 IF (NHEP.GT.NMXHEP) CALL HWWARN('HWMEVT',103,*999)
733
734 JCL=NHEP
735
736 IDHW(JCL)=19
737
738 IDHEP(JCL)=91
739
740 IF (NCL.LT.3) THEN
741
742 ISTHEP(JCL)=170+NCL
743
744 ID1=INID(1,NCL)
745
746 ID2=INID(2,NCL)
747
748 ELSE
749
750 ID1=ID2-6
751
752 IF (NCL.EQ.3) ID1=ID3
753
754 ID2=HWRINT(7,8)
755
756 ISTHEP(JCL)=173
757
758 ENDIF
759
760 JMOHEP(1,JCL)=ICMS
761
762 JMOHEP(2,JCL)=0
763
764 CALL HWVZRO(3,PHEP(1,JCL))
765
766 PHEP(4,JCL)=RMASS(ID1)+RMASS(ID2)+PMBM1+HWREXP(TWO/PMBM2)
767
768 PHEP(5,JCL)=PHEP(4,JCL)
769
770C---HADRONIZE AND DECAY CLUSTERS
771
772 CALL HWCFLA(ID1,ID2,JD1,JD2)
773
774 CALL HWCHAD(JCL,JD1,JD2,JD3)
775
776 IF (IERROR.NE.0) RETURN
777
778 IF (JD3.EQ.0) THEN
779
780 EMCL=RMASS(IDHW(NHEP))
781
782 IF (PHEP(4,JCL).NE.EMCL) THEN
783
784 PHEP(4,JCL)=EMCL
785
786 PHEP(5,JCL)=EMCL
787
788 PHEP(4,NHEP)=EMCL
789
790 PHEP(5,NHEP)=EMCL
791
792 ENDIF
793
794 ELSE
795
796 EMCL=PHEP(5,JCL)
797
798 ENDIF
799
800 IDCL(NCL)=JD3
801
802 PPCL(5,NCL)=EMCL
803
804 SUMM=SUMM +EMCL
805
806 CALL HWDHAD
807
808 CALL HWDHVY
809
810 IF (IERROR.NE.0) RETURN
811
812C---CHECK CHARGED MULTIPLICITY
813
814 MODC=0
815
816 DO 50 KHEP=JCL,NHEP
817
818 IF (ISTHEP(KHEP).EQ.1) THEN
819
820 ICH=ICHRG(IDHW(KHEP))
821
822 IF (ICH.NE.0) THEN
823
824 MCHT=MCHT+ABS(ICH)
825
826 MODC=MODC+ICH
827
828 ENDIF
829
830 ENDIF
831
832 50 CONTINUE
833
834 IF (NCL.EQ.1) THEN
835
836 NCHT=NPPBAR+NETC+ABS(MODC)
837
838 GOTO 35
839
840 ELSEIF (NCL.EQ.2) THEN
841
842 NCHT=NCHT+ABS(MODC)
843
844 IF (NCHT.LT.0) NCHT=NCHT+2
845
846 ENDIF
847
848 IF (MCHT.LT.NCHT) THEN
849
850 GOTO 35
851
852 ELSEIF (MCHT.GT.NCHT) THEN
853
854 IF (MOD(NTRY,50).EQ.0) GOTO 25
855
856 IF (NTRY.LT.NSTRY) GOTO 30
857
858C---NO PHASE SPACE FOR SOFT EVENT
859
860 NHEP=ICMS-1
861
862 IF (IPROC.EQ.8000) THEN
863
864C---MINIMUM BIAS: RELABEL BEAM AND TARGET CLUSTERS
865
866 DO 60 IBT=1,2
867
868 KHEP=INHEP(IBT)
869
870 LHEP=JMOHEP(1,KHEP)
871
872 ISTHEP(KHEP)=1
873
874 IDHEP(KHEP)=IDHEP(LHEP)
875
876 IDHW(KHEP)=IDHW(LHEP)
877
878 60 CONTINUE
879
880 ELSE
881
882C---UNDERLYING EVENT: DECAY THEM
883
884 ISTHEP(INHEP(1))=163
885
886 ISTHEP(INHEP(2))=163
887
888 CALL HWCFOR
889
890 CALL HWCDEC
891
892 CALL HWDHAD
893
894 CALL HWDHVY
895
896 ENDIF
897
898 GOTO 90
899
900 ENDIF
901
902C---GENERATE CLUSTER MOMENTA IN CLUSTER CM
903
904C FRAME. N.B. SECOND CLUSTER IS TARGET
905
906 IF (SUMM.GT.TECM) GOTO 25
907
908 CALL HWMLPS(TECM)
909
910 IF (NCL.EQ.0) GOTO 25
911
912 JCL=0
913
914C---ROTATE & BOOST CLUSTERS & DECAY PRODUCTS
915
916 CALL HWULOF(PHEP(1,ICMS),PHEP(1,INHEP(1)),BMP)
917
918 CALL HWUROT(BMP, ONE,ZERO,BMR)
919
920C---BMR PUTS BEAM ALONG Z AXIS (WE WANT INVERSE)
921
922 DO 70 KHEP=ICMS+1,NHEP
923
924 IF (ISTHEP(KHEP).GT.180.AND.ISTHEP(KHEP).LT.190
925
926 $ .AND.JMOHEP(1,KHEP).EQ.ICMS) THEN
927
928 ISTHEP(KHEP)=ISTHEP(KHEP)+3
929
930 LHEP=KHEP
931
932 JCL=JCL+1
933
934 CALL HWUROB(BMR,PPCL(1,JCL),PPCL(1,JCL))
935
936 CALL HWULOB(PHEP(1,ICMS),PPCL(1,JCL),PPCL(1,JCL))
937
938C---NOW PPCL(*,JCL) IS LAB MOMENTUM OF JTH CLUSTER
939
940 ENDIF
941
942 CALL HWULOB(PPCL(1,JCL),PHEP(1,KHEP),PHEP(1,KHEP))
943
944 70 CONTINUE
945
946 ISTHEP(INHEP(1))=167
947
948 ISTHEP(INHEP(2))=168
949
950 JMOHEP(1,ICMS)=INHEP(1)
951
952 JMOHEP(2,ICMS)=INHEP(2)
953
954 JDAHEP(1,INHEP(1))=ICMS
955
956 JDAHEP(2,INHEP(1))=0
957
958 JDAHEP(1,INHEP(2))=ICMS
959
960 JDAHEP(2,INHEP(2))=0
961
962 JDAHEP(1,ICMS)=ICMS+1
963
964 JDAHEP(2,ICMS)=LHEP
965
966 90 CONTINUE
967
968 990 ISTAT=100
969
970 999 END
971
972CDECK ID>, HWMLPS.
973
974*CMZ :- -04/05/99 14.17.04 by Bryan Webber
975
976*-- Author : David Ward, modified by Bryan Webber
977
978C-----------------------------------------------------------------------
979
980 SUBROUTINE HWMLPS(TECM)
981
982C-----------------------------------------------------------------------
983
984C GENERATES CYLINDRICAL PHASE SPACE USING THE METHOD OF JADACH
985
986C RETURNS WITH NCL=0 IF UNSUCCESSFUL
987
988C-----------------------------------------------------------------------
989
990 INCLUDE 'HERWIG61.INC'
991
992 DOUBLE PRECISION HWREXT,HWRUNG,HWUSQR,TECM,ESS,ALOGS,EPS,SUMX,
993
994 & SUMY,PT,PX,PY,PT2,SUMPT2,SUMTM,XIMIN,XIMAX,YY,SUM1,SUM2,SUM3,
995
996 & SUM4,EX,FY,DD,DYY,ZZ,E1,TM,SLOP,XI(NMXCL)
997
998 INTEGER NTRY,I,NIT,IY(NMXCL),IDP
999
1000 EXTERNAL HWREXT,HWRUNG,HWUSQR
1001
1002 IF (NCL.GT.NMXCL) THEN
1003
1004 CALL HWWARN('HWMLPS',1,*999)
1005
1006 NCL=NMXCL
1007
1008 ENDIF
1009
1010 ESS=TECM**2
1011
1012 ALOGS=LOG(ESS)
1013
1014 EPS=1D-10/NCL
1015
1016 NTRY=0
1017
1018 11 NTRY=NTRY+1
1019
1020 IF (NTRY.GT.NSTRY) THEN
1021
1022 NCL=0
1023
1024 RETURN
1025
1026 ENDIF
1027
1028 SUMX=0.
1029
1030 SUMY=0.
1031
1032 DO 12 I=1,NCL
1033
1034C---Pt distribution of form exp(-b*Mt)
1035
1036C---Factors for pt slopes to fit data. IDCL contains the type of
1037
1038C q-qbar pair produced in this cluster (0 if 1-particle cluster).
1039
1040 IDP=IDCL(I)
1041
1042 IF (IDP.LE.2) THEN
1043
1044 SLOP=PMBP1
1045
1046 ELSEIF(IDP.EQ.3.OR.IDP.EQ.10) THEN
1047
1048 SLOP=PMBP2
1049
1050 ELSEIF(IDP.GT.3.AND.IDP.LE.9) THEN
1051
1052 SLOP=PMBP3
1053
1054 ELSE
1055
1056 CALL HWWARN('HWMLPS',200,*999)
1057
1058 ENDIF
1059
1060 PT=HWREXT(PPCL(5,I),SLOP)
1061
1062 PT=HWUSQR(PT**2-PPCL(5,I)**2)
1063
1064 CALL HWRAZM(PT,PX,PY)
1065
1066 PPCL(1,I)=PX
1067
1068 PPCL(2,I)=PY
1069
1070 SUMX=SUMX+PPCL(1,I)
1071
1072 12 SUMY=SUMY+PPCL(2,I)
1073
1074 SUMX=SUMX/NCL
1075
1076 SUMY=SUMY/NCL
1077
1078 SUMPT2=0.
1079
1080 SUMTM=0.
1081
1082 DO 13 I=1,NCL
1083
1084 PPCL(1,I)=PPCL(1,I)-SUMX
1085
1086 PPCL(2,I)=PPCL(2,I)-SUMY
1087
1088 PT2=PPCL(1,I)**2+PPCL(2,I)**2
1089
1090 SUMPT2=SUMPT2+PT2
1091
1092C---STORE TRANSVERSE MASS IN PPCL(3,I) TEMPORARILY
1093
1094 PPCL(3,I)=SQRT(PT2+PPCL(5,I)**2)
1095
1096 13 SUMTM=SUMTM+PPCL(3,I)
1097
1098 IF (SUMTM.GT.TECM) GOTO 11
1099
1100 DO 14 I=1,NCL
1101
1102C---Form of "reduced rapidity" distribution
1103
1104 XI(I)=HWRUNG(0.6*ONE,ONE)
1105
1106 14 CONTINUE
1107
1108 CALL HWUSOR(XI,NCL,IY,1)
1109
1110 XIMIN=XI(1)
1111
1112 XIMAX=XI(NCL)-XI(1)
1113
1114C---N.B. TARGET CLUSTER IS SECOND
1115
1116 XI(1)=0.
1117
1118 DO 16 I=NCL-1,2,-1
1119
1120 XI(I+1)=(XI(I)-XIMIN)/XIMAX
1121
1122 16 CONTINUE
1123
1124 XI(2)=1.
1125
1126 YY=LOG(ESS/(PPCL(3,1)*PPCL(3,2)))
1127
1128 DO 18 NIT=1,10
1129
1130 SUM1=0.
1131
1132 SUM2=0.
1133
1134 SUM3=0.
1135
1136 SUM4=0.
1137
1138 DO 19 I=1,NCL
1139
1140 TM=PPCL(3,I)
1141
1142 EX=EXP(YY*XI(I))
1143
1144 SUM1=SUM1+(TM*EX)
1145
1146 SUM2=SUM2+(TM/EX)
1147
1148 SUM3=SUM3+(TM*EX)*XI(I)
1149
1150 19 SUM4=SUM4+(TM/EX)*XI(I)
1151
1152 FY=ALOGS-LOG(SUM1*SUM2)
1153
1154 DD=(SUM3*SUM2-SUM1*SUM4)/(SUM1*SUM2)
1155
1156 DYY=FY/DD
1157
1158 IF(ABS(DYY/YY).LT.EPS) GOTO 20
1159
1160 18 YY=YY+DYY
1161
1162C---Y ITERATIONS EXCEEDED - TRY AGAIN
1163
1164 IF (NTRY.LT.100) GOTO 11
1165
1166 EPS=10.*EPS
1167
1168 IF (EPS.GT.ONE) CALL HWWARN('HWMLPS',100,*999)
1169
1170 CALL HWWARN('HWMLPS',50,*11)
1171
1172 20 YY=YY+DYY
1173
1174 ZZ=LOG(TECM/SUM1)
1175
1176 DO 22 I=1,NCL
1177
1178 TM=PPCL(3,I)
1179
1180 E1=EXP(ZZ+YY*XI(I))
1181
1182 PPCL(3,I)=(0.5*TM)*((1./E1)-E1)
1183
1184 PPCL(4,I)=(0.5*TM)*((1./E1)+E1)
1185
1186 22 CONTINUE
1187
1188 999 END
1189
1190CDECK ID>, HWMNBI.
1191
1192*CMZ :- -26/04/91 11.11.55 by Bryan Webber
1193
1194*-- Author : David Ward, modified by Bryan Webber
1195
1196C-----------------------------------------------------------------------
1197
1198 FUNCTION HWMNBI(N,AVNCH,EK)
1199
1200C-----------------------------------------------------------------------
1201
1202C---Computes negative binomial probability
1203
1204C-----------------------------------------------------------------------
1205
1206 DOUBLE PRECISION HWMNBI,AVNCH,EK,R
1207
1208 INTEGER N,I
1209
1210 IF(N.LE.0) THEN
1211
1212 HWMNBI=0
1213
1214 ELSE
1215
1216 R=AVNCH/EK
1217
1218 HWMNBI=(1.+R)**(-EK)
1219
1220 R=R/(1.+R)
1221
1222 DO 1 I=1,N
1223
1224 HWMNBI=HWMNBI*R*(EK+I-1)/I
1225
1226 1 CONTINUE
1227
1228 ENDIF
1229
1230 END