1 #include "isajet/pilot.h"
2 SUBROUTINE SSMHC(MHCNEG)
3 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/.
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.
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
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
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
32 C-----------------------------------------------------------------------
33 #if defined(CERNLIB_IMPNONE)
36 #include "isajet/sslun.inc"
37 #include "isajet/sssm.inc"
38 #include "isajet/sspar.inc"
41 REAL PI,PI2,SR2,G2,GP2,GGP,GG1,GG2
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
56 C For non-degenerate squarks
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
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))*
90 MTQ=SSMQCD(DBLE(MTMT),DBLE(HIGFRZ))
139 C Charged Higgs mass calculation
140 C (AAT and AAB are also assumed to be real)
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.
168 C Non-degenerate squarks and/or D-terms. Since D-terms are
169 C always included, old dead code has been deleted. FEP
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.'
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.'
200 C Calculate 2M1*B term
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
213 C Calculate derivatives w.r.t. H_R divided by v*sqrt(2)
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
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
227 C Calculate derivatives w.r.t. H_R' divided by v'*sqrt(2)
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
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
239 C Calculate second derivatives w.r.t. H_R^+
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
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
251 C Calculate second derivatives w.r.t. H_R'^-
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
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
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
265 C Find stop eigenvector factors
267 TEMPT=-TLRM/2.0 - ZAP*VPVM*GG1/8.0
270 BTA1=-MT*(EP*COTB + ATR)
272 IF(ATA1.EQ.0.0 .AND. BTA1.EQ.0.0) THEN
275 IF(ATA1.EQ.0.0 .AND. BTA1.EQ.0.0) THEN
276 WRITE(LOUT,*) 'SSMHC: ERROR: ZERO EIGENVECTOR FOR MST1SQ'
280 IF(ATA2.EQ.0.0 .AND. BTA2.EQ.0.0) THEN
283 IF(ATA2.EQ.0.0 .AND. BTA2.EQ.0.0) THEN
284 WRITE(LOUT,*) 'SSMHC: ERROR: ZERO EIGENVECTOR FOR MST2SQ'
292 NABT1=1.0/(ATA1SQ+BTA1SQ)
293 NABT2=1.0/(ATA2SQ+BTA2SQ)
295 C Find sbottom eigenvector factors
297 TEMPB=-BLRM/2.0 + ZAP*VPVM*GG2/8.0
300 BBA1=-MB*(EP*TANB + ABR)
302 IF(ABA1.EQ.0.0 .AND. BBA1.EQ.0.0) THEN
305 IF(ABA1.EQ.0.0 .AND. BBA1.EQ.0.0) THEN
306 WRITE(LOUT,*) 'SSMHC: ERROR: ZERO EIGENVECTOR FOR MSB1SQ'
310 IF(ABA2.EQ.0.0 .AND. BBA2.EQ.0.0) THEN
313 IF(ABA2.EQ.0.0 .AND. BBA2.EQ.0.0) THEN
314 WRITE(LOUT,*) 'SSMHC: ERROR ZERO EIGENVECTOR FOR MSB2SQ'
322 NABB1=1.0/(ABA1SQ+BBA1SQ)
323 NABB2=1.0/(ABA2SQ+BBA2SQ)
325 C Calculate perturbative terms
326 C from H_R^+2 derivative terms
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
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
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
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
380 IF(MST1SQ.EQ.MSB2SQ) THEN
381 WRITE(LOUT,*) 'SSMHC: WARNING: MST1 = MSB2',
382 $ ' MST1 RAISED BY 0.0001'
383 MST1SQ = MST1SQ*1.0001
385 IF(MST2SQ.EQ.MSB1SQ) THEN
386 WRITE(LOUT,*) 'SSMHC: WARNING: MST2 = MSB1',
387 $ ' MST2 RAISED BY 0.0001'
388 MST2SQ = MST2SQ*1.0001
390 IF(MST2SQ.EQ.MSB2SQ) THEN
391 WRITE(LOUT,*) 'SSMHC: WARNING: MST2 = MSB2',
392 $ ' MST2 RAISED BY 0.0001'
393 MST2SQ = MST2SQ*1.0001
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)
401 C Calculate perturbative terms
402 C from H_R'^-2 derivative terms
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
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
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
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
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)
455 C Calculate perturbative terms
456 C from H_R^+ H_R'^- derivative terms
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
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
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
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
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)
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)
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)
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)
525 EM3=LMST1*PDPMT1+LMST2*PDPMT2
526 EM3=EM3+LMSB1*PDPMB1+LMSB2*PDPMB2
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
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
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
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