ISAJET/isasusy/ssmhc.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE SSMHC(MHCNEG)
   3 C-----------------------------------------------------------------------
   4 C
   5 C         Calculates charged Higgs mass
   6 C         (scalar Higgs mixing angle) using radiative
   7 C         corrections calculated by M. Bisset
   8 C         and save results in /SSPAR/.
   9 C
  10 C         Both top and bottom couplings are now
  11 C         included.  Non-degenerate mixed squark
  12 C         masses and A-terms are also included.
  13 C         The D-terms from the squark mass matrix
  14 C         (terms prop. to g**2 * Yukawa coupling)
  15 C         are included by default.
  16 C
  17 C         There is an arbitrary mass scale that must
  18 C         chosen to avoid dimensionful logarithms.
  19 C         The choice does not matter if D-terms are
  20 C         not included, but it does matter if D-terms
  21 C         are included.
  22 C
  23 C         Arbitrary mass scale set to
  24 C              QQQ = HIGFRZ = SQRT(AMTLSS*AMTRSS)
  25 C         Updated to include running masses as 2-loop effect
  26 C
  27 C         It is assumed that the A-terms are real.
  28 C         (Complex A-terms are taken into account
  29 C         much of the subroutine; but, not in all
  30 C         cases.)
  31 C
  32 C-----------------------------------------------------------------------
  33 #if defined(CERNLIB_IMPNONE)
  34       IMPLICIT NONE
  35 #endif
  36 #include "isajet/sslun.inc"
  37 #include "isajet/sssm.inc"
  38 #include "isajet/sspar.inc"
  39 C
  40       INTEGER MHCNEG
  41       REAL PI,PI2,SR2,G2,GP2,GGP,GG1,GG2
  42       REAL GGPM,GG3,GG4
  43       REAL TANB,COTB,COSB,SINB,BE
  44       REAL SINB2,COSB2,COS2B,SIN2B
  45       REAL V2,VP2,V,VP,VVP,VPVM,VVPP,MT,MB
  46       REAL MT2,MB2,MT4,MB4,FT2,FB2,FT,FB
  47       REAL MW2,ZAP,QQQ2,EP,EP2,RR,MHP2
  48       REAL ATI,ABI,ATR,ABR,AT2,AB2
  49       REAL MSTL2,MSTR2,MSBL2,MSBR2
  50       REAL TLRM,BLRM,TLRP,BLRP
  51       REAL MST1,MST2,MSB1,MSB2
  52       REAL MST1SQ,MST2SQ,MSB1SQ,MSB2SQ
  53       REAL TTT1,TTT2,TTT3,BBB1,BBB2,BBB3
  54       REAL TEMPT,TEMPB,ROOTT,ROOTB,TM1B
  55 C
  56 C      For non-degenerate squarks
  57 C
  58       REAL TERMT,TERMB,DHRT1,DHRT2,DHRB1,DHRB2
  59       REAL DHRPT1,DHRPT2,DHRPB1,DHRPB2
  60       REAL DHPST1,DHPST2,DHPSB1,DHPSB2
  61       REAL DHMST1,DHMST2,DHMSB1,DHMSB2
  62       REAL ATA1,ATA2,BTA1,BTA2,ABA1,ABA2,BBA1,BBA2
  63       REAL ATA1SQ,ATA2SQ,BTA1SQ,BTA2SQ
  64       REAL ABA1SQ,ABA2SQ,BBA1SQ,BBA2SQ
  65       REAL NABT1,NABT2,NABB1,NABB2
  66       REAL FTG,FBG,BABA,PT1B1,PT1B2,PT2B1,PT2B2
  67       REAL PDPST1,PDPST2,PDPSB1,PDPSB2
  68       REAL PDMST1,PDMST2,PDMSB1,PDMSB2
  69       REAL PDPMT1,PDPMT2,PDPMB1,PDPMB2
  70       REAL LMST1,LMST2,LMSB1,LMSB2
  71       REAL EMI1,EMI2,EM3,TEMPF
  72       REAL DVRR,DVRPRP,DVRRP,TRACE,DETV
  73       REAL TERMSQ,GOLD2,MHC2
  74       REAL HIGFRZ,ASMB,MBMB,MBQ,ASMT,MTMT,MTQ,SUALFS
  75       DOUBLE PRECISION SSMQCD
  76 C
  77       MHCNEG=0
  78       PI=4.*ATAN(1.)
  79       PI2=PI**2
  80       SR2=SQRT(2.)
  81       G2=4.*PI*ALFAEM/SN2THW
  82       GP2=G2*SN2THW/(1.-SN2THW)
  83       HIGFRZ=SQRT(AMTLSS*AMTRSS)
  84       ASMB=SUALFS(AMBT**2,.36,AMTP,3)
  85       MBMB=AMBT*(1.-4*ASMB/3./PI)
  86       MBQ=SSMQCD(DBLE(MBMB),DBLE(HIGFRZ))
  87       ASMT=SUALFS(AMTP**2,.36,AMTP,3)
  88       MTMT=AMTP/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/AMTP))*
  89      $(ASMT/PI)**2)
  90       MTQ=SSMQCD(DBLE(MTMT),DBLE(HIGFRZ))
  91       MT=MTQ
  92       MB=MBQ
  93       MT2=MT**2
  94       MB2=MB**2
  95       MT4=MT2**2
  96       MB4=MB2**2
  97       MW2=AMW**2
  98       EP=TWOM1
  99       EP2=EP**2
 100       RR=RV2V1
 101       MHP2=AMHA**2
 102       TANB=1.0/RR
 103       COTB=RR
 104       BE=ATAN(1./RV2V1)
 105       SINB=SIN(BE)
 106       COSB=COS(BE)
 107       SINB2=SINB**2
 108       COSB2=COSB**2
 109       SIN2B=SIN(2.*BE)
 110       COS2B=COS(2.*BE)
 111       V2=2.0*MW2*SINB2/G2
 112       VP2=2.0*MW2*COSB2/G2
 113       V=SQRT(V2)
 114       VP=SQRT(VP2)
 115       VVP=SQRT(V2*VP2)
 116       VPVM=VP2-V2
 117       GGP=G2+GP2
 118       GGPM=G2-GP2
 119       GG1=G2-5.0*GP2/3.0
 120       GG2=G2-GP2/3.0
 121       GG3=G2+5.0*GP2/3.0
 122       GG4=G2+GP2/3.0
 123 C
 124       VVPP=2.0*AMZ**2/GGP
 125       FT2=MT2/V2
 126       FB2=MB2/VP2
 127       FT=SQRT(FT2)
 128       FB=SQRT(FB2)
 129 C
 130       MSTL2=AMTLSS**2
 131       MSTR2=AMTRSS**2
 132       MSBL2=AMBLSS**2
 133       MSBR2=AMBRSS**2
 134       TLRM=MSTL2-MSTR2
 135       BLRM=MSBL2-MSBR2
 136       TLRP=MSTL2+MSTR2
 137       BLRP=MSBL2+MSBR2
 138 C
 139 C          Charged Higgs mass calculation
 140 C          (AAT and AAB are also assumed to be real)
 141 C
 142       ATR=AAT
 143       ABR=AAB
 144       ATI=0.0
 145       ABI=0.0
 146       AT2=ATR**2+ATI**2
 147       AB2=ABR**2+ABI**2
 148 C
 149 C      UNFORTUNATELY, I HAVE USED MY OLD CONVENTION
 150 C      FOR THE STOP AND SBOTTOM EIGENVALUES HERE
 151 C      (T1 <==> T2 OF NOTATION IN X. TATA'S AND OTHER
 152 C      PEOPLE'S NOTATION).  SO THE NEXT EIGHT LINES ARE
 153 C      A FIX-UP UNTIL I GO THROUGH AND CHANGE THE
 154 C      NOTATION THROUGHOUT THIS SUBROUTINE.
 155 C
 156       MST2=AMT1SS
 157       MST1=AMT2SS
 158       MSB2=AMB1SS
 159       MSB1=AMB2SS
 160       MST2SQ=AMT1SS**2
 161       MST1SQ=AMT2SS**2
 162       MSB2SQ=AMB1SS**2
 163       MSB1SQ=AMB2SS**2
 164 C
 165       QQQ2=HIGFRZ**2
 166       ZAP=1
 167 C
 168 C          Non-degenerate squarks and/or D-terms. Since D-terms are
 169 C          always included, old dead code has been deleted. FEP
 170 C
 171 C          ROOTT recast as a sum of squares. Note that ROOTT=0
 172 C          could happen accidently and causes an error.
 173       TTT1=0.5*(MSTL2+MSTR2)+MT2+ZAP*VPVM*GGP/8.0
 174       TTT2=TLRM+ZAP*GG1*VPVM/4.0
 175       TTT3=4.0*FT2*(EP2*VP2+2.0*EP*VVP*ATR+AT2*V2)
 176       ROOTT=4*MT2*(COSB*EP + AAT*SINB)**2/SINB**2 +
 177      $(AMTLSS**2 - AMTRSS**2 +
 178      $AMW**2*(-5*GP2/3 + G2)*(COSB**2 - SINB**2)/(2*G2))**2
 179       ROOTT=0.5*SQRT(ROOTT)
 180       IF(ROOTT.EQ.0.0) THEN
 181         WRITE(LOUT,*) 'SSMHC: ERROR: ROOTT = 0,',
 182      $  '  CANNOT CALCULATE H+ MASS FOR THIS CASE.'
 183         STOP99
 184       ENDIF
 185 C
 186       BBB1=0.5*(MSBL2+MSBR2)+MB2-ZAP*VPVM*GGP/8.0
 187       BBB2=BLRM-ZAP*GG2*VPVM/4.0
 188       BBB3=4.0*FB2*(EP2*V2+2.0*EP*VVP*ABR+AB2*VP2)
 189 C          ROOTB recast as a sum of squares.
 190       ROOTB=4*MB2*(AAB*COSB + EP*SINB)**2/COSB**2 +
 191      $(AMBLSS**2 - AMBRSS**2 -
 192      $AMW**2*(-GP2/3 + G2)*(COSB**2 - SINB**2)/(2*G2))**2
 193       ROOTB=0.5*SQRT(ROOTB)
 194       IF(ROOTB.EQ.0.0) THEN
 195         WRITE(LOUT,*) 'SSMHC: ERROR: ROOTB = 0,',
 196      $  '  CANNOT CALCULATE H+ MASS FOR THIS CASE.'
 197         STOP99
 198       ENDIF
 199 C
 200 C      Calculate 2M1*B term
 201 C
 202       TEMPT=EP*FT2*VVP*ATI**2/(ROOTT**2)
 203       TEMPB=EP*FB2*VVP*ABI**2/(ROOTB**2)
 204       TM1B=-FT2*(TEMPT+ATR)*TTT1*LOG(MST1SQ/MST2SQ)/ROOTT
 205       TM1B=TM1B-FT2*ATR*LOG(MST1SQ*MST2SQ/QQQ2/QQQ2)
 206       TM1B=TM1B+FT2*(2.0*TEMPT-ATR)
 207       TM1B=TM1B-FB2*(TEMPB+ABR)*BBB1*LOG(MSB1SQ/MSB2SQ)/ROOTB
 208       TM1B=TM1B-FB2*ABR*LOG(MSB1SQ*MSB2SQ/QQQ2/QQQ2)
 209       TM1B=TM1B+FB2*(2.0*TEMPB-ABR)
 210       TM1B=3.0*EP*TM1B/32.0/PI2
 211       TM1B=TM1B-VVP*MHP2/VVPP
 212 C
 213 C       Calculate derivatives w.r.t. H_R divided by v*sqrt(2)
 214 C
 215       TEMPT=ZAP*GG1*(TLRM+VPVM*GG1/4.0)/8.0
 216       TERMT=-TEMPT+FT2*(EP*COTB*ATR+AT2)
 217       TERMT=TERMT/(2.0*ROOTT)
 218       DHRT1=FT2-ZAP*GGP/8.0 + TERMT
 219       DHRT2=FT2-ZAP*GGP/8.0 - TERMT
 220 C
 221       TEMPB=ZAP*GG2*(BLRM-VPVM*GG2/4.0)/8.0
 222       TERMB=TEMPB+FB2*EP*(EP+COTB*ABR)
 223       TERMB=TERMB/(2.0*ROOTB)
 224       DHRB1=ZAP*GGP/8.0 + TERMB
 225       DHRB2=ZAP*GGP/8.0 - TERMB
 226 C
 227 C       Calculate derivatives w.r.t. H_R' divided by v'*sqrt(2)
 228 C
 229       TERMT=TEMPT+FT2*EP*(EP+TANB*ATR)
 230       TERMT=TERMT/(2.0*ROOTT)
 231       DHRPT1=ZAP*GGP/8.0 + TERMT
 232       DHRPT2=ZAP*GGP/8.0 - TERMT
 233 C
 234       TERMB=-TEMPB+FB2*(EP*TANB*ABR+AB2)
 235       TERMB=TERMB/(2.0*ROOTB)
 236       DHRPB1=FB2-ZAP*GGP/8.0 + TERMB
 237       DHRPB2=FB2-ZAP*GGP/8.0 - TERMB
 238 C
 239 C       Calculate second derivatives w.r.t. H_R^+
 240 C
 241       TEMPT=(TLRM+ZAP*VPVM*GG1/4.0)/(4.0*ROOTT)
 242       TERMT=TEMPT*(-FT2+ZAP*GG3/4.0)
 243       DHPST1=FT2/2.0 + ZAP*GGPM/8.0 + TERMT
 244       DHPST2=FT2/2.0 + ZAP*GGPM/8.0 - TERMT
 245 C
 246       TEMPB=(BLRM-ZAP*VPVM*GG2/4.0)/(4.0*ROOTB)
 247       TERMB=TEMPB*(FT2-ZAP*GG4/4.0)
 248       DHPSB1=FT2/2.0 - ZAP*GGPM/8.0 + TERMB
 249       DHPSB2=FT2/2.0 - ZAP*GGPM/8.0 - TERMB
 250 C
 251 C       Calculate second derivatives w.r.t. H_R'^-
 252 C
 253       TERMT=TEMPT*(FB2-ZAP*GG3/4.0)
 254       DHMST1=FB2/2.0 - ZAP*GGPM/8.0 + TERMT
 255       DHMST2=FB2/2.0 - ZAP*GGPM/8.0 - TERMT
 256 C
 257       TERMB=TEMPB*(-FB2+ZAP*GG4/4.0)
 258       DHMSB1=FB2/2.0 + ZAP*GGPM/8.0 + TERMB
 259       DHMSB2=FB2/2.0 + ZAP*GGPM/8.0 - TERMB
 260 C
 261 C       From perturbative terms
 262 C       Here I assume A_t and A_b are real
 263 C       and therefore the eigenvectors are real
 264 C
 265 C       Find stop eigenvector factors
 266 C
 267       TEMPT=-TLRM/2.0 - ZAP*VPVM*GG1/8.0
 268       ATA1=TEMPT+ROOTT
 269       ATA2=TEMPT-ROOTT
 270       BTA1=-MT*(EP*COTB + ATR)
 271       BTA2=BTA1
 272       IF(ATA1.EQ.0.0 .AND. BTA1.EQ.0.0) THEN
 273         ATA1=-BTA1
 274         BTA1 = TEMPT - ROOTT
 275         IF(ATA1.EQ.0.0 .AND. BTA1.EQ.0.0) THEN
 276           WRITE(LOUT,*) 'SSMHC: ERROR: ZERO EIGENVECTOR FOR MST1SQ'
 277           STOP99
 278         ENDIF
 279       ENDIF
 280       IF(ATA2.EQ.0.0 .AND. BTA2.EQ.0.0) THEN
 281         ATA2=-BTA2
 282         BTA2=TEMPT+ROOTT
 283         IF(ATA2.EQ.0.0 .AND. BTA2.EQ.0.0) THEN
 284           WRITE(LOUT,*) 'SSMHC: ERROR:  ZERO EIGENVECTOR FOR MST2SQ'
 285           STOP99
 286         ENDIF
 287       ENDIF
 288       ATA1SQ=ATA1**2
 289       BTA1SQ=BTA1**2
 290       ATA2SQ=ATA2**2
 291       BTA2SQ=BTA2**2
 292       NABT1=1.0/(ATA1SQ+BTA1SQ)
 293       NABT2=1.0/(ATA2SQ+BTA2SQ)
 294 C
 295 C       Find sbottom eigenvector factors
 296 C
 297       TEMPB=-BLRM/2.0 + ZAP*VPVM*GG2/8.0
 298       ABA1=TEMPB+ROOTB
 299       ABA2=TEMPB-ROOTB
 300       BBA1=-MB*(EP*TANB + ABR)
 301       BBA2=BBA1
 302       IF(ABA1.EQ.0.0 .AND. BBA1.EQ.0.0) THEN
 303         ABA1=-BBA1
 304         BBA1=TEMPB-ROOTB
 305         IF(ABA1.EQ.0.0 .AND. BBA1.EQ.0.0) THEN
 306           WRITE(LOUT,*) 'SSMHC: ERROR: ZERO EIGENVECTOR FOR MSB1SQ'
 307           STOP99
 308         ENDIF
 309       ENDIF
 310       IF(ABA2.EQ.0.0 .AND. BBA2.EQ.0.0) THEN
 311         ABA2=-BBA2
 312         BBA2=TEMPB+ROOTB
 313         IF(ABA2.EQ.0.0 .AND. BBA2.EQ.0.0) THEN
 314           WRITE(LOUT,*) 'SSMHC: ERROR ZERO EIGENVECTOR FOR MSB2SQ'
 315           STOP99
 316         ENDIF
 317       ENDIF
 318       ABA1SQ=ABA1**2
 319       BBA1SQ=BBA1**2
 320       ABA2SQ=ABA2**2
 321       BBA2SQ=BBA2**2
 322       NABB1=1.0/(ABA1SQ+BBA1SQ)
 323       NABB2=1.0/(ABA2SQ+BBA2SQ)
 324 C
 325 C       Calculate perturbative terms
 326 C        from H_R^+2 derivative terms
 327 C
 328       FTG=FT2-ZAP*G2/2.0
 329       BABA=FT*FB*(VVP*FTG-EP*ATR)
 330       PT1B1=V2*(FTG**2)*BTA1SQ*BBA1SQ
 331       PT1B1=PT1B1+2.0*EP*FB*V*FTG*BTA1SQ*BBA1*ABA1
 332       PT1B1=PT1B1-2.0*ATR*MT*FTG*BTA1*ATA1*BBA1SQ
 333       PT1B1=PT1B1+2.0*BABA*BTA1*ATA1*BBA1*ABA1
 334       PT1B1=PT1B1-2.0*ATR*FT2*MB*ATA1SQ*BBA1*ABA1
 335       PT1B1=PT1B1+2.0*EP*FT*FB*MB*BTA1*ATA1*ABA1SQ
 336       PT1B1=PT1B1+FT2*AT2*ATA1SQ*BBA1SQ
 337       PT1B1=PT1B1+FB2*EP2*BTA1SQ*ABA1SQ
 338       PT1B1=PT1B1+FT2*MB2*ATA1SQ*ABA1SQ
 339       PT1B1=PT1B1*NABT1*NABB1
 340 C
 341       PT1B2=V2*(FTG**2)*BTA1SQ*BBA2SQ
 342       PT1B2=PT1B2+2.0*EP*FB*V*FTG*BTA1SQ*BBA2*ABA2
 343       PT1B2=PT1B2-2.0*ATR*MT*FTG*BTA1*ATA1*BBA2SQ
 344       PT1B2=PT1B2+2.0*BABA*BTA1*ATA1*BBA2*ABA2
 345       PT1B2=PT1B2-2.0*ATR*FT2*MB*ATA1SQ*BBA2*ABA2
 346       PT1B2=PT1B2+2.0*EP*FT*FB*MB*BTA1*ATA1*ABA2SQ
 347       PT1B2=PT1B2+FT2*AT2*ATA1SQ*BBA2SQ
 348       PT1B2=PT1B2+FB2*EP2*BTA1SQ*ABA2SQ
 349       PT1B2=PT1B2+FT2*MB2*ATA1SQ*ABA2SQ
 350       PT1B2=PT1B2*NABT1*NABB2
 351 C
 352       PT2B1=V2*(FTG**2)*BTA2SQ*BBA1SQ
 353       PT2B1=PT2B1+2.0*EP*FB*V*FTG*BTA2SQ*BBA1*ABA1
 354       PT2B1=PT2B1-2.0*ATR*MT*FTG*BTA2*ATA2*BBA1SQ
 355       PT2B1=PT2B1+2.0*BABA*BTA2*ATA2*BBA1*ABA1
 356       PT2B1=PT2B1-2.0*ATR*FT2*MB*ATA2SQ*BBA1*ABA1
 357       PT2B1=PT2B1+2.0*EP*FT*FB*MB*BTA2*ATA2*ABA1SQ
 358       PT2B1=PT2B1+FT2*AT2*ATA2SQ*BBA1SQ
 359       PT2B1=PT2B1+FB2*EP2*BTA2SQ*ABA1SQ
 360       PT2B1=PT2B1+FT2*MB2*ATA2SQ*ABA1SQ
 361       PT2B1=PT2B1*NABT2*NABB1
 362 C
 363       PT2B2=V2*(FTG**2)*BTA2SQ*BBA2SQ
 364       PT2B2=PT2B2+2.0*EP*FB*V*FTG*BTA2SQ*BBA2*ABA2
 365       PT2B2=PT2B2-2.0*ATR*MT*FTG*BTA2*ATA2*BBA2SQ
 366       PT2B2=PT2B2+2.0*BABA*BTA2*ATA2*BBA2*ABA2
 367       PT2B2=PT2B2-2.0*ATR*FT2*MB*ATA2SQ*BBA2*ABA2
 368       PT2B2=PT2B2+2.0*EP*FT*FB*MB*BTA2*ATA2*ABA2SQ
 369       PT2B2=PT2B2+FT2*AT2*ATA2SQ*BBA2SQ
 370       PT2B2=PT2B2+FB2*EP2*BTA2SQ*ABA2SQ
 371       PT2B2=PT2B2+FT2*MB2*ATA2SQ*ABA2SQ
 372       PT2B2=PT2B2*NABT2*NABB2
 373 C          The following used to add 1e-8, but this may be less than
 374 C          machine precision. Multiply by 1.001 instead.
 375       IF(MST1SQ.EQ.MSB1SQ) THEN
 376         WRITE(LOUT,*) 'SSMHC: WARNING: MST1 = MSB1',
 377      $  '  MST1 RAISED BY 0.0001'
 378         MST1SQ = MST1SQ*1.0001
 379       ENDIF
 380       IF(MST1SQ.EQ.MSB2SQ) THEN
 381         WRITE(LOUT,*) 'SSMHC: WARNING: MST1 = MSB2',
 382      $  '  MST1 RAISED BY 0.0001'
 383         MST1SQ = MST1SQ*1.0001
 384       ENDIF
 385       IF(MST2SQ.EQ.MSB1SQ) THEN
 386         WRITE(LOUT,*) 'SSMHC: WARNING: MST2 = MSB1',
 387      $  '  MST2 RAISED BY 0.0001'
 388         MST2SQ = MST2SQ*1.0001
 389       ENDIF
 390       IF(MST2SQ.EQ.MSB2SQ) THEN
 391         WRITE(LOUT,*) 'SSMHC: WARNING: MST2 = MSB2',
 392      $  '  MST2 RAISED BY 0.0001'
 393         MST2SQ = MST2SQ*1.0001
 394       ENDIF
 395 C
 396       PDPST1=PT1B1/(MST1SQ-MSB1SQ) +PT1B2/(MST1SQ-MSB2SQ)
 397       PDPST2=PT2B1/(MST2SQ-MSB1SQ) +PT2B2/(MST2SQ-MSB2SQ)
 398       PDPSB1=PT1B1/(MSB1SQ-MST1SQ) +PT2B1/(MSB1SQ-MST2SQ)
 399       PDPSB2=PT1B2/(MSB2SQ-MST1SQ) +PT2B2/(MSB2SQ-MST2SQ)
 400 C
 401 C       Calculate perturbative terms
 402 C        from H_R'^-2 derivative terms
 403 C
 404       FBG=FB2-ZAP*G2/2.0
 405       BABA=FT*FB*(VVP*FBG-EP*ABR)
 406       PT1B1=VP2*(FBG**2)*BTA1SQ*BBA1SQ
 407       PT1B1=PT1B1-2.0*ABR*MB*FBG*BTA1SQ*BBA1*ABA1
 408       PT1B1=PT1B1+2.0*EP*FT*VP*FBG*BTA1*ATA1*BBA1SQ
 409       PT1B1=PT1B1+2.0*BABA*BTA1*ATA1*BBA1*ABA1
 410       PT1B1=PT1B1+2.0*EP*FT*FB*MT*ATA1SQ*BBA1*ABA1
 411       PT1B1=PT1B1-2.0*ABR*FB2*MT*BTA1*ATA1*ABA1SQ
 412       PT1B1=PT1B1+FT2*EP2*ATA1SQ*BBA1SQ
 413       PT1B1=PT1B1+FB2*AB2*BTA1SQ*ABA1SQ
 414       PT1B1=PT1B1+FB2*MT2*ATA1SQ*ABA1SQ
 415       PT1B1=PT1B1*NABT1*NABB1
 416 C
 417       PT1B2=VP2*(FBG**2)*BTA1SQ*BBA2SQ
 418       PT1B2=PT1B2-2.0*ABR*MB*FBG*BTA1SQ*BBA2*ABA2
 419       PT1B2=PT1B2+2.0*EP*FT*VP*FBG*BTA1*ATA1*BBA2SQ
 420       PT1B2=PT1B2+2.0*BABA*BTA1*ATA1*BBA2*ABA2
 421       PT1B2=PT1B2+2.0*EP*FT*FB*MT*ATA1SQ*BBA2*ABA2
 422       PT1B2=PT1B2-2.0*ABR*FB2*MT*BTA1*ATA1*ABA2SQ
 423       PT1B2=PT1B2+FT2*EP2*ATA1SQ*BBA2SQ
 424       PT1B2=PT1B2+FB2*AB2*BTA1SQ*ABA2SQ
 425       PT1B2=PT1B2+FB2*MT2*ATA1SQ*ABA2SQ
 426       PT1B2=PT1B2*NABT1*NABB2
 427 C
 428       PT2B1=VP2*(FBG**2)*BTA2SQ*BBA1SQ
 429       PT2B1=PT2B1-2.0*ABR*MB*FBG*BTA2SQ*BBA1*ABA1
 430       PT2B1=PT2B1+2.0*EP*FT*VP*FBG*BTA2*ATA2*BBA1SQ
 431       PT2B1=PT2B1+2.0*BABA*BTA2*ATA2*BBA1*ABA1
 432       PT2B1=PT2B1+2.0*EP*FT*FB*MT*ATA2SQ*BBA1*ABA1
 433       PT2B1=PT2B1-2.0*ABR*FB2*MT*BTA2*ATA2*ABA1SQ
 434       PT2B1=PT2B1+FT2*EP2*ATA2SQ*BBA1SQ
 435       PT2B1=PT2B1+FB2*AB2*BTA2SQ*ABA1SQ
 436       PT2B1=PT2B1+FB2*MT2*ATA2SQ*ABA1SQ
 437       PT2B1=PT2B1*NABT2*NABB1
 438 C
 439       PT2B2=VP2*(FBG**2)*BTA2SQ*BBA2SQ
 440       PT2B2=PT2B2-2.0*ABR*MB*FBG*BTA2SQ*BBA2*ABA2
 441       PT2B2=PT2B2+2.0*EP*FT*VP*FBG*BTA2*ATA2*BBA2SQ
 442       PT2B2=PT2B2+2.0*BABA*BTA2*ATA2*BBA2*ABA2
 443       PT2B2=PT2B2+2.0*EP*FT*FB*MT*ATA2SQ*BBA2*ABA2
 444       PT2B2=PT2B2-2.0*ABR*FB2*MT*BTA2*ATA2*ABA2SQ
 445       PT2B2=PT2B2+FT2*EP2*ATA2SQ*BBA2SQ
 446       PT2B2=PT2B2+FB2*AB2*BTA2SQ*ABA2SQ
 447       PT2B2=PT2B2+FB2*MT2*ATA2SQ*ABA2SQ
 448       PT2B2=PT2B2*NABT2*NABB2
 449 C
 450       PDMST1=PT1B1/(MST1SQ-MSB1SQ) +PT1B2/(MST1SQ-MSB2SQ)
 451       PDMST2=PT2B1/(MST2SQ-MSB1SQ) +PT2B2/(MST2SQ-MSB2SQ)
 452       PDMSB1=PT1B1/(MSB1SQ-MST1SQ) +PT2B1/(MSB1SQ-MST2SQ)
 453       PDMSB2=PT1B2/(MSB2SQ-MST1SQ) +PT2B2/(MSB2SQ-MST2SQ)
 454 C
 455 C        Calculate perturbative terms
 456 C        from  H_R^+ H_R'^- derivative terms
 457 C
 458       BABA=FT*FB*(V2*FTG+VP2*FBG+EP2+ATR*ABR)
 459       PT1B1=-VVP*FTG*FBG*BTA1SQ*BBA1SQ
 460       PT1B1=PT1B1+FB*(ABR*V*FTG-EP*VP*FBG)*BTA1SQ*BBA1*ABA1
 461       PT1B1=PT1B1-FT*(EP*V*FTG-ATR*VP*FBG)*BTA1*ATA1*BBA1SQ
 462       PT1B1=PT1B1-BABA*BTA1*ATA1*BBA1*ABA1
 463       PT1B1=PT1B1+FT2*FB*(ATR*V-EP*VP)*ATA1SQ*BBA1*ABA1
 464       PT1B1=PT1B1+FT*FB2*(ABR*VP-EP*V)*BTA1*ATA1*ABA1SQ
 465       PT1B1=PT1B1+FT2*EP*ATR*ATA1SQ*BBA1SQ
 466       PT1B1=PT1B1+FB2*EP*ABR*BTA1SQ*ABA1SQ
 467       PT1B1=PT1B1-FT*FB*MT*MB*ATA1SQ*ABA1SQ
 468       PT1B1=PT1B1*NABT1*NABB1
 469 C
 470       PT1B2=-VVP*FTG*FBG*BTA1SQ*BBA2SQ
 471       PT1B2=PT1B2+FB*(ABR*V*FTG-EP*VP*FBG)*BTA1SQ*BBA2*ABA2
 472       PT1B2=PT1B2-FT*(EP*V*FTG-ATR*VP*FBG)*BTA1*ATA1*BBA2SQ
 473       PT1B2=PT1B2-BABA*BTA1*ATA1*BBA2*ABA2
 474       PT1B2=PT1B2+FT2*FB*(ATR*V-EP*VP)*ATA1SQ*BBA2*ABA2
 475       PT1B2=PT1B2+FT*FB2*(ABR*VP-EP*V)*BTA1*ATA1*ABA2SQ
 476       PT1B2=PT1B2+FT2*EP*ATR*ATA1SQ*BBA2SQ
 477       PT1B2=PT1B2+FB2*EP*ABR*BTA1SQ*ABA2SQ
 478       PT1B2=PT1B2-FT*FB*MT*MB*ATA1SQ*ABA2SQ
 479       PT1B2=PT1B2*NABT1*NABB2
 480 C
 481       PT2B1= -VVP*FTG*FBG*BTA2SQ*BBA1SQ
 482       PT2B1= PT2B1 + FB*(ABR*V*FTG - EP*VP*FBG)*BTA2SQ*BBA1*ABA1
 483       PT2B1= PT2B1 - FT*(EP*V*FTG - ATR*VP*FBG)*BTA2*ATA2*BBA1SQ
 484       PT2B1= PT2B1 - BABA*BTA2*ATA2*BBA1*ABA1
 485       PT2B1= PT2B1 + FT2*FB*(ATR*V - EP*VP)*ATA2SQ*BBA1*ABA1
 486       PT2B1= PT2B1 + FT*FB2*(ABR*VP - EP*V)*BTA2*ATA2*ABA1SQ
 487       PT2B1= PT2B1 + FT2*EP*ATR*ATA2SQ*BBA1SQ
 488       PT2B1= PT2B1 + FB2*EP*ABR*BTA2SQ*ABA1SQ
 489       PT2B1= PT2B1 - FT*FB*MT*MB*ATA2SQ*ABA1SQ
 490       PT2B1= PT2B1*NABT2*NABB1
 491 C
 492       PT2B2=-VVP*FTG*FBG*BTA2SQ*BBA2SQ
 493       PT2B2=PT2B2+FB*(ABR*V*FTG-EP*VP*FBG)*BTA2SQ*BBA2*ABA2
 494       PT2B2=PT2B2-FT*(EP*V*FTG-ATR*VP*FBG)*BTA2*ATA2*BBA2SQ
 495       PT2B2=PT2B2-BABA*BTA2*ATA2*BBA2*ABA2
 496       PT2B2=PT2B2+FT2*FB*(ATR*V-EP*VP)*ATA2SQ*BBA2*ABA2
 497       PT2B2=PT2B2+FT*FB2*(ABR*VP-EP*V)*BTA2*ATA2*ABA2SQ
 498       PT2B2=PT2B2+FT2*EP*ATR*ATA2SQ*BBA2SQ
 499       PT2B2=PT2B2+FB2*EP*ABR*BTA2SQ*ABA2SQ
 500       PT2B2=PT2B2-FT*FB*MT*MB*ATA2SQ*ABA2SQ
 501       PT2B2=PT2B2*NABT2*NABB2
 502 C
 503       PDPMT1=PT1B1/(MST1SQ-MSB1SQ) +PT1B2/(MST1SQ-MSB2SQ)
 504       PDPMT2=PT2B1/(MST2SQ-MSB1SQ) +PT2B2/(MST2SQ-MSB2SQ)
 505       PDPMB1=PT1B1/(MSB1SQ-MST1SQ) +PT2B1/(MSB1SQ-MST2SQ)
 506       PDPMB2=PT1B2/(MSB2SQ-MST1SQ) +PT2B2/(MSB2SQ-MST2SQ)
 507 C
 508       MST1SQ=MST1**2
 509       MST2SQ=MST2**2
 510       LMST1=MST1SQ*(2.0*LOG(MST1SQ/QQQ2)+1.0)
 511       LMST2=MST2SQ*(2.0*LOG(MST2SQ/QQQ2)+1.0)
 512       LMSB1=MSB1SQ*(2.0*LOG(MSB1SQ/QQQ2)+1.0)
 513       LMSB2=MSB2SQ*(2.0*LOG(MSB2SQ/QQQ2)+1.0)
 514 C
 515       EMI1=LMST1*(PDPST1+DHPST1-DHRT1)
 516       EMI1=EMI1+LMST2*(PDPST2+DHPST2-DHRT2)
 517       EMI1=EMI1+LMSB1*(PDPSB1+DHPSB1-DHRB1)
 518       EMI1=EMI1+LMSB2*(PDPSB2+DHPSB2-DHRB2)
 519 C
 520       EMI2=LMST1*(PDMST1+DHMST1-DHRPT1)
 521       EMI2=EMI2+LMST2*(PDMST2+DHMST2-DHRPT2)
 522       EMI2=EMI2+LMSB1*(PDMSB1+DHMSB1-DHRPB1)
 523       EMI2=EMI2+LMSB2*(PDMSB2+DHMSB2-DHRPB2)
 524 C
 525       EM3=LMST1*PDPMT1+LMST2*PDPMT2
 526       EM3=EM3+LMSB1*PDPMB1+LMSB2*PDPMB2
 527 C
 528       TEMPF=MT2*LOG(MT2/QQQ2)-MB2*LOG(MB2/QQQ2)
 529       DVRR=4.0*FT2*MB2*TEMPF/(MT2-MB2)
 530       DVRR=3.0*(EMI1-DVRR-2.0*FT2*MB2)/32.0/PI2
 531       DVRR=-TM1B*COTB +VP2*G2/2.0 +DVRR
 532 C
 533       DVRPRP=4.0*FB2*MT2*TEMPF/(MT2-MB2)
 534       DVRPRP=3.0*(EMI2-DVRPRP-2.0*FB2*MT2)/32.0/PI2
 535       DVRPRP=-TM1B*TANB +V2*G2/2.0 +DVRPRP
 536 C
 537       DVRRP=1.0 +(MT2+MB2)*LOG(MT2/MB2)/(MT2-MB2)
 538       DVRRP=2.0*FT*FB*MT*MB*(DVRRP+LOG(MT2*MB2/(QQQ2**2)))
 539       DVRRP=3.0*(EM3+DVRRP)/32.0/PI2
 540       DVRRP=TM1B -G2*VVP/2.0 +DVRRP
 541 C
 542       TRACE=DVRR+DVRPRP
 543       DETV=4.0*(DVRR*DVRPRP-DVRRP**2)
 544 C          Rewrite TERMSQ=TRACE**2-DETV
 545       TERMSQ=(DVRR-DVRPRP)**2+4*DVRRP**2
 546       TERMSQ=SQRT(TERMSQ)/2.0
 547       GOLD2=TRACE/2.0 -TERMSQ
 548       MHC2=TRACE/2.0 +TERMSQ
 549 C
 550       IF(MHC2.LT.0.0) THEN
 551         MHCNEG=1
 552         AMHC=0.
 553         GO TO 1000
 554       ENDIF
 555       AMHC=SQRT(MHC2)
 556 1000  RETURN
 557       END