]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HERWIG/src/hwhigw.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
use TMath::Abs() instead of ambiguous abs().
[u/mrichter/AliRoot.git] / HERWIG / src / hwhigw.f
CommitLineData
3820ca8e 1
2CDECK ID>, HWHIGW.
3
4*CMZ :- -26/04/91 14.55.44 by Federico Carminati
5
6*-- Author : Mike Seymour
7
8C-----------------------------------------------------------------------
9
10 SUBROUTINE HWHIGW
11
12C-----------------------------------------------------------------------
13
14C HIGGS PRODUCTION VIA W BOSON FUSION
15
16C MEAN EVWGT = HIGGS PRODN C-S * BRANCHING FRACTION IN NB
17
18C-----------------------------------------------------------------------
19
20 INCLUDE 'HERWIG61.INC'
21
22 DOUBLE PRECISION HWULDO,HWRUNI,HWR,HWUAEM,K1MAX2,K1MIN2,K12,
23
24 & K2MAX2,K2MIN2,K22,EMW2,EMW,ROOTS,EMH2,EMH,ROOTS2,P1,PHI1,PHI2,
25
26 & COSPHI,COSTH1,SINTH1,COSTH2,SINTH2,P2,WEIGHT,TAU,TAULN,CSFAC,
27
28 & PSUM,PROB,Q1(5),Q2(5),H(5),A,B,C,TERM2,BRHIGQ,G1WW,G2WW,G1ZZ(6),
29
30 & G2ZZ(6),AWW,AZZ(6),PWW,PZZ(6),EMZ,EMZ2,RSUM,GLUSQ,GRUSQ,GLDSQ,
31
32 & GRDSQ,GLESQ,GRESQ,CW,CZ,EMFAC,CV,CA,BR,X2,ETA,P1JAC,FACTR,EH2
33
34 INTEGER HWRINT,IDEC,I,ID1,ID2,IHAD
35
36 LOGICAL EE,EP
37
38 EXTERNAL HWULDO,HWRUNI,HWR,HWUAEM,HWRINT
39
40 SAVE EMW2,EMZ2,EE,GLUSQ,GRUSQ,GLDSQ,GRDSQ,GLESQ,GRESQ,G1ZZ,G2ZZ,
41
42 & G1WW,G2WW,CW,CZ,PSUM,AWW,PWW,AZZ,PZZ,ROOTS,Q1,Q2,H,FACTR
43
44 EQUIVALENCE (EMW,RMASS(198)),(EMZ,RMASS(200))
45
46 IHAD=2
47
48 IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
49
50 IF (FSTWGT) THEN
51
52 EMW2=EMW**2
53
54 EMZ2=EMZ**2
55
56 GLUSQ=(VFCH(2,1)+AFCH(2,1))**2
57
58 GRUSQ=(VFCH(2,1)-AFCH(2,1))**2
59
60 GLDSQ=(VFCH(1,1)+AFCH(1,1))**2
61
62 GRDSQ=(VFCH(1,1)-AFCH(1,1))**2
63
64 GLESQ=(VFCH(11,1)+AFCH(11,1))**2
65
66 GRESQ=(VFCH(11,1)-AFCH(11,1))**2
67
68 G1ZZ(1)=GLUSQ*GLUSQ+GRUSQ*GRUSQ
69
70 G2ZZ(1)=GLUSQ*GRUSQ+GRUSQ*GLUSQ
71
72 G1ZZ(2)=GLUSQ*GLDSQ+GRUSQ*GRDSQ
73
74 G2ZZ(2)=GLUSQ*GRDSQ+GRUSQ*GLDSQ
75
76 G1ZZ(3)=GLDSQ*GLDSQ+GRDSQ*GRDSQ
77
78 G2ZZ(3)=GLDSQ*GRDSQ+GRDSQ*GLDSQ
79
80 G1ZZ(4)=GLESQ*GLESQ+GRESQ*GRESQ
81
82 G2ZZ(4)=GLESQ*GRESQ+GRESQ*GLESQ
83
84 G1ZZ(5)=GLESQ*GLUSQ+GRESQ*GRUSQ
85
86 G2ZZ(5)=GLESQ*GRUSQ+GRESQ*GLUSQ
87
88 G1ZZ(6)=GLESQ*GLDSQ+GRESQ*GRDSQ
89
90 G2ZZ(6)=GLESQ*GRDSQ+GRESQ*GLDSQ
91
92 G1WW=0.25
93
94 G2WW=0
95
96 FACTR=GEV2NB/(128.*PIFAC**3)
97
98 EH2=RMASS(201)**2
99
100 CW=256*(PIFAC*HWUAEM(EH2)/SWEIN)**3*EMW2
101
102 CZ=256.*(PIFAC*HWUAEM(EH2))**3*EMZ2/(SWEIN*(1.-SWEIN))
103
104 ENDIF
105
106 EE=IPRO.LT.10
107
108 EP=IPRO.GE.90
109
110 IF (.NOT.GENEV) THEN
111
112C---CHOOSE PARAMETERS
113
114 EVWGT=0.
115
116 CALL HWHIGM(EMH,EMFAC)
117
118 IF (EMH.LE.ZERO .OR. EMH.GE.PHEP(5,3)) RETURN
119
120 EMSCA=EMH
121
122 IF (EE) THEN
123
124 ROOTS=PHEP(5,3)
125
126 ELSE
127
128 TAU=(EMH/PHEP(5,3))**2
129
130 TAULN=LOG(TAU)
131
132 ROOTS=PHEP(5,3)*SQRT(EXP(HWRUNI(0,-1D-10,TAULN)))
133
134 ENDIF
135
136 EMH2=EMH**2
137
138 ROOTS2=ROOTS**2
139
140C---CHOOSE P1 ACCORDING TO (1-ETA)*(ETA-X2)/ETA**2
141
142C WHERE ETA=1-2P1/ROOTS AND X2=EMH**2/S
143
144 X2=EMH2/ROOTS2
145
146 1 ETA=X2**HWR()
147
148 IF (HWR()*(1-EMH/ROOTS)**2*ETA.GT.(1-ETA)*(ETA-X2))GOTO 1
149
150 P1JAC=0.5*ROOTS*ETA**2/((1-ETA)*(ETA-X2))
151
152 & *(-LOG(X2)*(1+X2)-2*(1-X2))
153
154 P1=0.5*ROOTS*(1-ETA)
155
156C---CHOOSE PHI1,2 UNIFORMLY
157
158 PHI1=2*PIFAC*HWR()
159
160 PHI2=2*PIFAC*HWR()
161
162 COSPHI=COS(PHI2-PHI1)
163
164C---CHOOSE K1^2, ON PROPAGATOR FACTOR
165
166 K1MAX2=2*P1*ROOTS
167
168 K1MIN2=0
169
170 K12=EMW2-(EMW2+K1MAX2)*(EMW2+K1MIN2)/
171
172 & ((K1MAX2-K1MIN2)*HWR()+(EMW2+K1MIN2))
173
174C---CALCULATE COSTH1 FROM K1^2
175
176 COSTH1=1+K12/(P1*ROOTS)
177
178 SINTH1=SQRT(1-COSTH1**2)
179
180C---CHOOSE K2^2
181
182 K2MAX2=ROOTS*(ROOTS2-EMH2-2*ROOTS*P1)*(ROOTS-P1-P1*COSTH1)
183
184 & /((ROOTS-P1)**2-(P1*COSTH1)**2-(P1*SINTH1*COSPHI)**2)
185
186 K2MIN2=0
187
188 K22=EMW2-(EMW2+K2MAX2)*(EMW2+K2MIN2)/
189
190 & ((K2MAX2-K2MIN2)*HWR()+(EMW2+K2MIN2))
191
192C---CALCULATE A,B,C FACTORS, AND...
193
194 A=-2*K22*P1*COSTH1 - ROOTS*(ROOTS2-EMH2-2*ROOTS*P1)
195
196 B=-2*K22*P1*SINTH1*COSPHI
197
198 C=+2*K22*P1 - 2*ROOTS*K22 - ROOTS*(ROOTS2-EMH2-2*ROOTS*P1)
199
200C---SOLVE A*COSTH2 + B*SINTH2 + C = 0 FOR COSTH2
201
202 TERM2=B**2 + A**2 - C**2
203
204 IF (TERM2.LT.ZERO) RETURN
205
206 TERM2=B*SQRT(TERM2)
207
208 IF (A.GE.ZERO) RETURN
209
210 COSTH2=(-A*C + TERM2)/(A**2+B**2)
211
212 SINTH2=SQRT(1-COSTH2**2)
213
214C---FINALLY, GET P2
215
216 IF (COSTH2.EQ.-ONE) RETURN
217
218 P2=-K22/(ROOTS*(1+COSTH2))
219
220C---LOAD UP CMF MOMENTA
221
222 Q1(1)=P1*SINTH1*COS(PHI1)
223
224 Q1(2)=P1*SINTH1*SIN(PHI1)
225
226 Q1(3)=P1*COSTH1
227
228 Q1(4)=P1
229
230 Q1(5)=0
231
232 Q2(1)=P2*SINTH2*COS(PHI2)
233
234 Q2(2)=P2*SINTH2*SIN(PHI2)
235
236 Q2(3)=P2*COSTH2
237
238 Q2(4)=P2
239
240 Q2(5)=0
241
242 H(1)=-Q1(1)-Q2(1)
243
244 H(2)=-Q1(2)-Q2(2)
245
246 H(3)=-Q1(3)-Q2(3)
247
248 H(4)=-Q1(4)-Q2(4)+ROOTS
249
250 CALL HWUMAS(H)
251
252C---CALCULATE MATRIX ELEMENTS SQUARED
253
254 AWW=ENHANC(10)**2 * CW*(ROOTS2/2*HWULDO(Q1,Q2)*G1WW
255
256 & +ROOTS2/4*P1*P2*(1+COSTH1)*(1-COSTH2)*G2WW)
257
258 DO 10 I=1,6
259
260 AZZ(I)=ENHANC(11)**2 * CZ*(ROOTS2/2*HWULDO(Q1,Q2)*G1ZZ(I)
261
262 & +ROOTS2/4*P1*P2*(1+COSTH1)*(1-COSTH2)*G2ZZ(I))
263
264 & *((K12-EMW2)/(K12-EMZ2)*(K22-EMW2)/(K22-EMZ2))**2
265
266 10 CONTINUE
267
268C---CALCULATE WEIGHT IN INTEGRAL
269
270 WEIGHT=FACTR*P2*P1JAC/(ROOTS2**2*HWULDO(H,Q2))
271
272 & *(K1MAX2-K1MIN2)/((K1MAX2+EMW2)*(K1MIN2+EMW2))
273
274 & *(K2MAX2-K2MIN2)/((K2MAX2+EMW2)*(K2MIN2+EMW2))
275
276 & * EMFAC
277
278 EMSCA=EMW
279
280 XXMIN=(ROOTS/PHEP(5,3))**2
281
282 XLMIN=LOG(XXMIN)
283
284C---INCLUDE BRANCHING RATIO OF HIGGS
285
286 IDEC=MOD(IPROC,100)
287
288 IF (IDEC.GT.0.AND.IDEC.LE.12) WEIGHT=WEIGHT*BRHIG(IDEC)
289
290 IF (IDEC.EQ.0) THEN
291
292 BRHIGQ=0
293
294 DO 20 I=1,6
295
296 20 BRHIGQ=BRHIGQ+BRHIG(I)
297
298 WEIGHT=WEIGHT*BRHIGQ
299
300 ENDIF
301
302 IF (IDEC.EQ.10) THEN
303
304 CALL HWDBOZ(198,ID1,ID2,CV,CA,BR,1)
305
306 CALL HWDBOZ(199,ID1,ID2,CV,CA,BR,1)
307
308 WEIGHT=WEIGHT*BR
309
310 ELSEIF (IDEC.EQ.11) THEN
311
312 CALL HWDBOZ(200,ID1,ID2,CV,CA,BR,1)
313
314 CALL HWDBOZ(200,ID1,ID2,CV,CA,BR,1)
315
316 WEIGHT=WEIGHT*BR
317
318 ENDIF
319
320 IF (EE) THEN
321
322 CSFAC=WEIGHT
323
324 PSUM=AWW+AZZ(4)
325
326 EVWGT=CSFAC*PSUM
327
328 ELSEIF (EP) THEN
329
330 CSFAC=-WEIGHT*TAULN
331
332 XX(1)=ONE
333
334 XX(2)=XXMIN
335
336 CALL HWSFUN(XX(2),EMSCA,IDHW(IHAD),NSTRU,DISF(1,2),2)
337
338 IF (IDHW(1).LE.126) THEN
339
340 PWW=(DISF(2,2)+DISF(4,2)+DISF(7,2)+DISF( 9,2))*AWW
341
342 ELSE
343
344 PWW=(DISF(1,2)+DISF(3,2)+DISF(8,2)+DISF(10,2))*AWW
345
346 ENDIF
347
348 PZZ(5)=(DISF(2,2)+DISF(4,2)+DISF(8,2)+DISF(10,2))*AZZ(5)
349
350 PZZ(6)=(DISF(1,2)+DISF(3,2)+DISF(7,2)+DISF( 9,2))*AZZ(6)
351
352 PSUM=PWW+PZZ(5)+PZZ(6)
353
354 EVWGT=CSFAC*PSUM
355
356 ELSE
357
358 CSFAC=WEIGHT*TAULN*XLMIN
359
360 CALL HWSGEN(.TRUE.)
361
362 PWW=((DISF(2,1)+DISF(4, 1)+DISF(7,1)+DISF(9,1))
363
364 & *(DISF(8,2)+DISF(10,2)+DISF(1,2)+DISF(3,2))
365
366 & +(DISF(8,1)+DISF(10,1)+DISF(1,1)+DISF(3,1))
367
368 & *(DISF(2,2)+DISF(4, 2)+DISF(7,2)+DISF(9,2)))
369
370 & *AWW
371
372 PZZ(1)=((DISF(2,1)+DISF(4,1)+DISF(8,1)+DISF(10,1))
373
374 & *(DISF(2,2)+DISF(4,2)+DISF(8,2)+DISF(10,2)))
375
376 & *AZZ(1)
377
378 PZZ(2)=((DISF(2,1)+DISF(4,1)+DISF(8,1)+DISF(10,1))
379
380 & *(DISF(1,2)+DISF(3,2)+DISF(7,2)+DISF(9, 2))
381
382 & +(DISF(1,1)+DISF(3,1)+DISF(7,1)+DISF(9, 1))
383
384 & *(DISF(2,2)+DISF(4,2)+DISF(8,2)+DISF(10,2)))
385
386 & *AZZ(2)
387
388 PZZ(3)=((DISF(1,1)+DISF(3,1)+DISF(7,1)+DISF(9,1))
389
390 & *(DISF(1,2)+DISF(3,2)+DISF(7,2)+DISF(9,2)))
391
392 & *AZZ(3)
393
394 PSUM=PWW+PZZ(1)+PZZ(2)+PZZ(3)
395
396C---EVENT WEIGHT IS SUM OVER ALL COMBINATIONS
397
398 EVWGT=CSFAC*PSUM
399
400 ENDIF
401
402 ELSE
403
404C---GENERATE EVENT
405
406C---CHOOSE EVENT TYPE
407
408 RSUM=PSUM*HWR()
409
410C---ELECTRON BEAMS?
411
412 IF (EE) THEN
413
414 IDN(1)=IDHW(1)
415
416 IDN(2)=IDHW(2)
417
418C---WW FUSION?
419
420 IF (RSUM.LT.AWW) THEN
421
422 IDN(3)=IDN(1)+1
423
424 IDN(4)=IDN(2)+1
425
426C---ZZ FUSION?
427
428 ELSE
429
430 IDN(3)=IDN(1)
431
432 IDN(4)=IDN(2)
433
434 ENDIF
435
436C---LEPTON-HADRON COLISION?
437
438 ELSEIF (EP) THEN
439
440C---WW FUSION?
441
442 IDN(1)=IDHW(1)
443
444 IF (RSUM.LT.PWW) THEN
445
446 24 IDN(2)=HWRINT(1,8)
447
448 IF (IDN(2).GE.5) IDN(2)=IDN(2)+2
449
450 IF (ICHRG(IDN(1))*ICHRG(IDN(2)).GT.0) GOTO 24
451
452 PROB=DISF(IDN(2),2)*AWW/PWW
453
454 IF (HWR().GT.PROB) GOTO 24
455
456 IDN(3)=IDN(1)+1
457
458 IF (HWR().GT.SCABI) THEN
459
460 IDN(4)= 4*INT((IDN(2)-1)/2)-IDN(2)+3
461
462 ELSE
463
464 IDN(4)=12*INT((IDN(2)-1)/6)-IDN(2)+5
465
466 ENDIF
467
468C---ZZ FUSION FROM U-TYPE QUARK?
469
470 ELSEIF (RSUM.LT.PWW+PZZ(5)) THEN
471
472 26 IDN(2)=2*HWRINT(1,4)
473
474 IF (IDN(2).GE.5) IDN(2)=IDN(2)+2
475
476 PROB=DISF(IDN(2),2)*AZZ(5)/PZZ(5)
477
478 IF (HWR().GT.PROB) GOTO 26
479
480 IDN(3)=IDN(1)
481
482 IDN(4)=IDN(2)
483
484C---ZZ FUSION FROM D-TYPE QUARK?
485
486 ELSE
487
488 28 IDN(2)=2*HWRINT(1,4)-1
489
490 IF (IDN(2).GE.5) IDN(2)=IDN(2)+2
491
492 PROB=DISF(IDN(2),2)*AZZ(6)/PZZ(6)
493
494 IF (HWR().GT.PROB) GOTO 28
495
496 IDN(3)=IDN(1)
497
498 IDN(4)=IDN(2)
499
500 ENDIF
501
502C---HADRON BEAMS?
503
504 ELSE
505
506C---WW FUSION?
507
508 IF (RSUM.LT.PWW) THEN
509
510 31 DO 32 I=1,2
511
512 IDN(I)=HWRINT(1,8)
513
514 IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
515
516 32 CONTINUE
517
518 IF (ICHRG(IDN(1))*ICHRG(IDN(2)).GT.0) GOTO 31
519
520 PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AWW/PWW
521
522 IF (HWR().GT.PROB) GOTO 31
523
524C---CHOOSE OUTGOING QUARKS
525
526 DO 33 I=1,2
527
528 IF (HWR().GT.SCABI) THEN
529
530 IDN(I+2)=4*INT((IDN(I)-1)/2)-IDN(I)+3
531
532 ELSE
533
534 IDN(I+2)=12*INT((IDN(I)-1)/6)-IDN(I)+5
535
536 ENDIF
537
538 33 CONTINUE
539
540C---ZZ FUSION FROM U-TYPE QUARKS?
541
542 ELSEIF (RSUM.LT.PWW+PZZ(1)) THEN
543
544 41 DO 42 I=1,2
545
546 IDN(I)=2*HWRINT(1,4)
547
548 IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
549
550 42 CONTINUE
551
552 PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AZZ(1)/PZZ(1)
553
554 IF (HWR().GT.PROB) GOTO 41
555
556 IDN(3)=IDN(1)
557
558 IDN(4)=IDN(2)
559
560C---ZZ FUSION FROM D-TYPE QUARKS?
561
562 ELSEIF (RSUM.LT.PWW+PZZ(1)+PZZ(3)) THEN
563
564 51 DO 52 I=1,2
565
566 IDN(I)=2*HWRINT(1,4)-1
567
568 IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
569
570 52 CONTINUE
571
572 PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AZZ(3)/PZZ(3)
573
574 IF (HWR().GT.PROB) GOTO 51
575
576 IDN(3)=IDN(1)
577
578 IDN(4)=IDN(2)
579
580C---ZZ FUSION FROM UD-TYPE PAIRS?
581
582 ELSE
583
584 61 IF (HWR().GT.HALF) THEN
585
586 IDN(1)=2*HWRINT(1,4)-1
587
588 IDN(2)=2*HWRINT(1,4)
589
590 ELSE
591
592 IDN(1)=2*HWRINT(1,4)
593
594 IDN(2)=2*HWRINT(1,4)-1
595
596 ENDIF
597
598 DO 62 I=1,2
599
600 62 IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
601
602 PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AZZ(2)/PZZ(2)
603
604 IF (HWR().GT.PROB) GOTO 61
605
606 IDN(3)=IDN(1)
607
608 IDN(4)=IDN(2)
609
610 ENDIF
611
612 ENDIF
613
614C---NOW BOOST TO LAB, AND SET UP STATUS CODES etc
615
616 IDCMF=15
617
618C---INCOMING
619
620 IF (.NOT.EE) CALL HWEONE
621
622C---CMF POINTERS
623
624 JDAHEP(1,NHEP)=NHEP+1
625
626 JDAHEP(2,NHEP)=NHEP+3
627
628 JMOHEP(1,NHEP+1)=NHEP
629
630 JMOHEP(1,NHEP+2)=NHEP
631
632 JMOHEP(1,NHEP+3)=NHEP
633
634C---OUTGOING MOMENTA (GIVE QUARKS MASS NON-COVARIANTLY!)
635
636 Q1(5)=RMASS(IDN(1))
637
638 Q1(4)=SQRT(Q1(4)**2+Q1(5)**2)
639
640 Q2(5)=RMASS(IDN(2))
641
642 Q2(4)=SQRT(Q2(4)**2+Q2(5)**2)
643
644 H(4)=-Q1(4)-Q2(4)+PHEP(5,NHEP)
645
646 CALL HWUMAS(H)
647
648 CALL HWULOB(PHEP(1,NHEP),Q1,PHEP(1,NHEP+1))
649
650 CALL HWULOB(PHEP(1,NHEP),Q2,PHEP(1,NHEP+2))
651
652 CALL HWULOB(PHEP(1,NHEP),H,PHEP(1,NHEP+3))
653
654C---STATUS AND IDs
655
656 ISTHEP(NHEP+1)=113
657
658 ISTHEP(NHEP+2)=114
659
660 ISTHEP(NHEP+3)=114
661
662 IDHW(NHEP+1)=IDN(3)
663
664 IDHEP(NHEP+1)=IDPDG(IDN(3))
665
666 IDHW(NHEP+2)=IDN(4)
667
668 IDHEP(NHEP+2)=IDPDG(IDN(4))
669
670 IDHW(NHEP+3)=201
671
672 IDHEP(NHEP+3)=IDPDG(201)
673
674C---COLOUR LABELS
675
676 JMOHEP(2,NHEP+1)=NHEP-2
677
678 JMOHEP(2,NHEP+2)=NHEP-1
679
680 JMOHEP(2,NHEP-1)=NHEP+2
681
682 JMOHEP(2,NHEP-2)=NHEP+1
683
684 JMOHEP(2,NHEP+3)=NHEP+3
685
686 JDAHEP(2,NHEP+1)=NHEP-2
687
688 JDAHEP(2,NHEP+2)=NHEP-1
689
690 JDAHEP(2,NHEP-1)=NHEP+2
691
692 JDAHEP(2,NHEP-2)=NHEP+1
693
694 JDAHEP(2,NHEP+3)=NHEP+3
695
696 NHEP=NHEP+3
697
698 ENDIF
699
700 999 END
701
702CDECK ID>, HWHIGY.
703
704*CMZ :- -26/04/91 13.37.37 by Federico Carminati
705
706*-- Author : Mike Seymour
707
708C-----------------------------------------------------------------------
709
710 FUNCTION HWHIGY(A,B,XP)
711
712C-----------------------------------------------------------------------
713
714C CALCULATE THE INTEGRAL OF BERENDS AND KLEISS APPENDIX B
715
716C-----------------------------------------------------------------------
717
718 IMPLICIT NONE
719
720 COMPLEX XQ,Z1,Z2,Z3,Z4,C0,C1,C2,C3,C4,C5,C6,C7,C8,FUN,Z
721
722 DOUBLE PRECISION HWHIGY,TWO,A,B,XP,Y
723
724 PARAMETER (TWO=2.D0)
725
726C---DECLARE ALL THE STATEMENT-FUNCTION DEFINITIONS
727
728 C0(Z,A)=(Z**2-A)**2*((Z**2+A)**2-24*Z*(Z**2+A)+8*Z**2*(A+6))/Z**4
729
730 C1(Z,A)=A**4/(3*Z)
731
732 C2(Z,A)=-A**3*(24*Z-A)/(2*Z**2)
733
734 C3(Z,A)=A**2*(8*Z**2*(A+6)-24*A*Z+A**2)/Z**3
735
736 C4(Z,A)=-A**2*(24*Z**3+8*Z**2*(A+6)-24*A*Z+A**2)/Z**4
737
738 C5(Z,A)=Z**3-24*Z**2+8*Z*(A+6)+24*A
739
740 C6(Z,A)=0.5*Z**2-12*Z+4*(A+6)
741
742 C7(Z,A)=Z/3-8
743
744 C8(Z,A)=0.25
745
746 FUN(Z,Y,A)=C0(Z,A)*LOG(Y-Z)
747
748 & +C1(Z,A)/Y**3
749
750 & +C2(Z,A)/Y**2
751
752 & +C3(Z,A)/Y
753
754 & +C4(Z,A)*LOG(Y)
755
756 & +C5(Z,A)*Y
757
758 & +C6(Z,A)*Y**2
759
760 & +C7(Z,A)*Y**3
761
762 & +C8(Z,A)*Y**4
763
764C---NOW EVALUATE THE INTEGRAL
765
766 HWHIGY=0
767
768 IF (A.GT.4) RETURN
769
770 XQ=CMPLX(XP,B)
771
772 Z1=XQ+SQRT(XQ**2-A)
773
774 Z2=XQ-SQRT(XQ**2-A)
775
776 Z3=FUN(Z1,TWO,A)-FUN(Z1,SQRT(A),A)
777
778 Z4=FUN(Z2,TWO,A)-FUN(Z2,SQRT(A),A)
779
780 HWHIGY=AIMAG((Z3-Z4)/(Z1-Z2))/(8*B)
781
782 END