1 #include "isajet/pilot.h"
2 C--------------------------------------------------------------------
3 SUBROUTINE SUGRA(M0,MHF,A0,TANB,SGNMU,MT,IMODEL)
4 C--------------------------------------------------------------------
6 C Calculate supergravity spectra for ISAJET using as inputs
7 C M0 = M_0 = common scalar mass at GUT scale
8 C MHF = M_(1/2) = common gaugino mass at GUT scale
9 C A0 = A_0 = trilinear soft breaking parameter at GUT scale
10 C TANB = tan(beta) = ratio of vacuum expectation values v_1/v_2
11 C SGNMU = sgn(mu) = +-1 = sign of Higgsino mass term
12 C MT = M_t = mass of t quark
13 C M0 = Lambda = ratio of vevs <F>/<S>
14 C MHF = M_Mes = messenger scale
15 C A0 = n_5 = number of messenger fields
16 C IMODEL = 1 for SUGRA model
20 C Uses Runge-Kutta method to integrate RGE's from M_Z to M_GUT
21 C and back, putting in correct thresholds. For the first iteration
22 C only the first 6 couplings are included and a common threshold
25 C See /SUGMG/ for definitions of couplings and masses.
27 #if defined(CERNLIB_IMPNONE)
30 #include "isajet/sslun.inc"
31 #include "isajet/sspar.inc"
32 #include "isajet/sssm.inc"
33 #include "isajet/sugxin.inc"
34 #include "isajet/sugmg.inc"
35 #include "isajet/sugpas.inc"
36 #include "isajet/sugnu.inc"
37 #include "isajet/ssinf.inc"
38 REAL GY(7),W1(21),G(29),W2(87)
41 DOUBLE PRECISION DDILOG,XLM
43 EXTERNAL SURG06,SURG26
44 REAL M0,MHF,A0,TANB,SGNMU,MT,XLAMGM,XMESGM,XN5GM
46 REAL M2,SUALFE,SUALFS,Q,T,A1I,AGUT,A3I,A2I,MTMT,ASMT,DT,
47 $TGUT,TZ,GGUT,SIG2,SIG1,MH1S,MH2S,AGUTI,
48 $MUS,MBMZ,MB,MTAU,MZ,MW,SR2,PI,ALEM,MTAMZ,
49 $MTAMB,MTAMTA,MBMB,ASMB,BETA,COTB,SINB,COS2B,COSB,XC,
50 $MSN,MG,MT1,MT2,MB1,MB2,MW1,MW2,AMU,BTHAT,BBHAT,BLHAT,AM2
52 REAL G0SAVE(26),DELG0,DELLIM,THRF,THRG,DY,QOLD
57 DATA MZ/91.187/,MTAU/1.777/,MB/4.9/,ALEM/.0078186/
58 C This choice is a compromise between precision and speed:
59 DATA MXITER/20/,NSTEP0/200/,DELLIM/2.E-2/
61 C Define REAL(COMPLEX*16) for g77. This might need to be
62 C changed for 64-bit machines?
66 C Save input parameters
83 IF (XGMIN(12).EQ.0.) XGMIN(12)=XN5GM
84 IF (XGMIN(13).EQ.0.) XGMIN(13)=XN5GM
85 IF (XGMIN(14).EQ.0.) XGMIN(14)=XN5GM
87 C Compute gauge mediated threshold functions
91 THRF=((1.D0+XLM)*(LOG(1.D0+XLM)-2*DDILOG(XLM/(1.D0+XLM))+
92 , .5*DDILOG(2*XLM/(1.D0+XLM)))+
93 , (1.D0-XLM)*(LOG(1.D0-XLM)-2*DDILOG(-XLM/(1.D0-XLM))+
94 , .5*DDILOG(-2*XLM/(1.D0-XLM))))/XLM**2
95 THRG=((1.D0+XLM)*LOG(1.D0+XLM)+(1.D0-XLM)*LOG(1.D0-XLM))/XLM**2
98 C Initialize standard model parameters in /SSSM/:
123 XW=.2324-1.03E-7*(MT**2-138.**2)
126 A1MZ=5*ALEM/3./(1.-XW)
129 GP=SQRT(3./5.*A1MZ*4.*PI)
137 IF (IMODEL.EQ.1) THEN
138 MSUSY=SQRT(M0**2+4*MHF**2)
139 ELSE IF (IMODEL.EQ.2) THEN
141 ELSE IF (IMODEL.EQ.7) THEN
142 MSUSY=SQRT(M0**2+(.01*MHF)**2)
144 C USE PIERCE PRESCRIPTION FOR MAGNITUDE OF VEV
145 C VEV=SR2*(248.6+0.9*LOG(MSUSY/AMZ)
146 C V=SQRT(VEV**2/(1.+COTB))
147 C PREVIOUS PRESCRIPTION
148 V=SQRT(2*MW**2/G2**2/(1.+COTB**2))
152 C Compute m(tau), m(b) at z scale using qcd, qed
154 MTAMTA=MTAU*(1.-SUALFE(MTAU**2)/PI)
155 MTAMB=MTAMTA*(SUALFE(MB**2)/SUALFE(MTAU**2))**(-27./76.)
156 MTAMZ=MTAMB*(SUALFE(MZ**2)/SUALFE(MB**2))**(-27./80.)
158 ASMB=SUALFS(MB**2,.36,MT,3)
159 MBMB=MB*(1.-4*ASMB/3./PI)
160 ASMZ=SUALFS(MZ**2,.36,MT,3)
161 MBMZ=MBMB*(ASMZ/ASMB)**(12./23.)*
162 $ (SUALFE(MZ**2)/SUALFE(MB**2))**(-3./80.)
164 ASMT=SUALFS(MT**2,.36,MT,3)
165 MTMT=MT/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/MT))*(ASMT/PI)**2)
167 FNMZ=SQRT(XNRIN(2)*XNRIN(1)/(SINB*VEV)**2)
170 C Run the 3 gauge and 3 Yukawa's up to find M_GUT ,A_GUT and
175 GY(1)=SQRT(4*PI*A1MZ)
176 GY(2)=SQRT(4*PI*A2MZ)
177 GY(3)=SQRT(4*PI*ALFA3)
182 IF (IMODEL.EQ.1.OR.IMODEL.EQ.7) THEN
183 IF (XSUGIN(7).EQ.0.) THEN
188 ELSE IF (IMODEL.EQ.2) THEN
193 DT=(TGUT-TZ)/FLOAT(NSTEP)
195 T=TZ+(TGUT-TZ)*FLOAT(II-1)/FLOAT(NSTEP)
197 IF (Q.GT.MT.AND.GY(6).EQ.0.) GY(6)=FTMT
198 IF (Q.GT.XNRIN(2).AND.GY(7).EQ.0.) GY(7)=FNMZ
199 CALL RKSTP(7,DT,T,GY,SURG06,W1)
203 IF (GY(5).GT.10..OR.GY(6).GT.10..OR.GY(7).GT.10.) THEN
207 IF (A1I.LT.A2I.AND.XSUGIN(7).EQ.0.) GO TO 10
209 IF (MGUT.EQ.1.E19) THEN
210 WRITE(LOUT,*) 'SUGRA: NO UNIFICATION FOUND'
213 10 IF (XSUGIN(7).EQ.0.) THEN
218 AGUT=(GY(1)**2/4./PI+GY(2)**2/4./PI)/2.
224 IF (XNRIN(1).EQ.0..AND.XNRIN(2).LT.1.E19) THEN
231 C Define parameters at GUT scale
234 IF (IMODEL.EQ.1) THEN
238 ELSE IF (IMODEL.EQ.2) THEN
240 G(J+6)=XGMIN(11+J)*XGMIN(8)*THRG*(GY(J)/4./PI)**2*XLAMGM
244 C OVERWRITE ALFA_3 UNIFICATION TO GET ALFA_3(MZ) RIGHT
245 IF (IMODEL.EQ.1.AND.IAL3UN.NE.0) G(3)=GGUT
249 C IF NR MAJORANA MASS EXISTS, SET EXTRA NR RGE PARAMETERS
250 IF (XNRIN(2).LT.1.E19) THEN
259 IF (IMODEL.EQ.1) THEN
263 C Set possible non-universal boundary conditions
265 IF (XNUSUG(J).LT.1.E19) THEN
270 IF (XNUSUG(J).LT.1.E19) THEN
274 ELSE IF (IMODEL.EQ.2) THEN
276 DY=SQRT(3./5.)*GY(1)*XGMIN(11)
277 G(13)=XC*(.75*XGMIN(13)*(GY(2)/4./PI)**4+.6*.25*
278 ,XGMIN(12)*(GY(1)/4./PI)**4)+XGMIN(9)-DY
279 G(14)=XC*(.75*XGMIN(13)*(GY(2)/4./PI)**4+.6*.25*
280 ,XGMIN(12)*(GY(1)/4./PI)**4)+XGMIN(10)+DY
281 G(15)=XC*(.6*XGMIN(12)*(GY(1)/4./PI)**4)+2*DY
282 G(16)=XC*(.75*XGMIN(13)*(GY(2)/4./PI)**4+.6*.25*
283 ,XGMIN(12)*(GY(1)/4./PI)**4)-DY
284 G(17)=XC*(4*XGMIN(14)*(GY(3)/4./PI)**4/3.+.6*XGMIN(12)*
285 ,(GY(1)/4./PI)**4/9.)+2*DY/3.
286 G(18)=XC*(4*XGMIN(14)*(GY(3)/4./PI)**4/3.+.6*4*XGMIN(12)*
287 ,(GY(1)/4./PI)**4/9.)-4*DY/3.
288 G(19)=XC*(4*XGMIN(14)*(GY(3)/4./PI)**4/3.+.75*XGMIN(13)*
289 ,(GY(2)/4./PI)**4+.6*XGMIN(12)*(GY(1)/4./PI)**4/36.)+DY/3.
295 ELSE IF (IMODEL.EQ.7) THEN
299 BLHAT=G(4)*(-9*G(1)**2/5.-3*G(2)**2+3*G(5)**2+4*G(4)**2)
300 BBHAT=G(5)*(-7*G(1)**2/15.-3*G(2)**2-16*G(3)**2/3.+
301 , G(6)**2+6*G(5)**2+G(4)**2)
302 BTHAT=G(6)*(-13*G(1)**2/15.-3*G(2)**2-16*G(3)**2/3.+
304 G(7)=-33*MHF*G(1)**2/5./16./PI**2
305 G(8)=-MHF*G(2)**2/16./PI**2
306 G(9)=3*MHF*G(3)**2/16./PI**2
307 G(10)=BLHAT*MHF/G(4)/16./PI**2
308 G(11)=BBHAT*MHF/G(5)/16./PI**2
309 G(12)=BTHAT*MHF/G(6)/16./PI**2
310 G(13)=(-99*G(1)**4/50.-3*G(2)**4/2.+3*G(5)*BBHAT+G(4)*BLHAT)*
311 , MHF**2/(16*PI**2)**2
312 G(14)=(-99*G(1)**4/50.-3*G(2)**4/2.+3*G(6)*BTHAT)*
313 , MHF**2/(16*PI**2)**2
314 G(15)=(-198*G(1)**4/25.)*MHF**2/(16*PI**2)**2
315 G(16)=(-99*G(1)**4/50.-3*G(2)**4/2.)*MHF**2/(16*PI**2)**2
316 G(17)=(-22*G(1)**4/25.+8*G(3)**4)*MHF**2/(16*PI**2)**2
317 G(18)=(-88*G(1)**4/25.+8*G(3)**4)*MHF**2/(16*PI**2)**2
318 G(19)=(-11*G(1)**4/50.-3*G(2)**4/2.+8*G(3)**4)*
319 , MHF**2/(16*PI**2)**2
320 G(20)=(-198*G(1)**4/25.+2*G(4)*BLHAT)*MHF**2/(16*PI**2)**2
321 G(21)=(-99*G(1)**4/50.-3*G(2)**4/2.+G(4)*BLHAT)*
322 , MHF**2/(16*PI**2)**2
323 G(22)=(-22*G(1)**4/25.+8*G(3)**4+2*G(5)*BBHAT)*
324 ,MHF**2/(16*PI**2)**2
325 G(23)=(-88*G(1)**4/25.+8*G(3)**4+2*G(6)*BTHAT)*
326 ,MHF**2/(16*PI**2)**2
327 G(24)=(-11*G(1)**4/50.-3*G(2)**4/2.+8*G(3)**4+G(5)*BBHAT+
328 , G(6)*BTHAT)*MHF**2/(16*PI**2)**2
337 C Check for tachyonic sleptons at GUT scale
338 IF (G(15).LT.0..OR.G(16).LT.0.) THEN
342 C Initialize thresholds
351 C Evolve parameters from mgut to mz
355 DT=(TZ-TGUT)/FLOAT(NSTEP)
356 C Freeze Higgs parameters at HIGFRZ = Drees' value
357 C AMTLSS, AMTRSS initialized to 0 for later use in HIGFRZ
358 IF (IMODEL.EQ.1) THEN
359 HIGFRZ=SQRT(M0**2+3*MHF**2)
360 ELSE IF (IMODEL.EQ.2) THEN
362 ELSE IF (IMODEL.EQ.7) THEN
363 HIGFRZ=SQRT(M0**2+(.01*MHF)**2)
368 T=TGUT+(TZ-TGUT)*FLOAT(II-1)/FLOAT(NSTEP)
371 CALL RKSTP(29,DT,T,G,SURG26,W2)
372 IF (Q.LT.AMNRMJ.AND.QOLD.GE.AMNRMJ.AND.FNMZ.EQ.0.) THEN
375 IF (Q.LT.AMNRMJ) THEN
380 CALL SUGFRZ(Q,G,G0,IG)
381 IF (NOGOOD.NE.0) GO TO 100
382 IF (Q.LT.MZ) GO TO 20
386 C Electroweak breaking constraints; tree level
387 MUS=(G0(13)-G0(14)*TANB**2)/(TANB**2-1.)-MZ**2/2.
392 MU=SQRT(MUS)*SIGN(1.,SGNMU)
393 B=(G0(13)+G0(14)+2*MUS)*SIN2B/MU/2.
394 C Compute tree level masses
395 CALL SUGMAS(G0,0,IMODEL)
396 IF (NOGOOD.NE.0) GO TO 100
397 C Compute effective potential corrections
398 CALL SUGEFF(G0,SIG1,SIG2)
401 MUS=(MH1S-MH2S*TANB**2)/(TANB**2-1.)-MZ**2/2.
406 MU=SQRT(MUS)*SIGN(1.,SGNMU)
407 B=(MH1S+MH2S+2*MUS)*SIN2B/MU/2.
409 C Recompute weak scale Yukawa couplings including SUSY loops
410 C Follow formulae of Pierce et al. NPB491, 3 (1997)
424 C Be careful in using our convention vs Pierce et al.
425 IF (ABS(COS(THETAT)).LT..707107) THEN
426 MTMT=MT/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/MT))*(ASMT/PI)**2
427 $ -ASMT/3./PI*(REAL8(SSB1(MT**2,MG,MT1))+
428 $ REAL8(SSB1(MT**2,MG,MT2))-SIN(2*THETAT)*MG/MT*
429 $ (REAL8(SSB0(MT**2,MG,MT1))-REAL8(SSB0(MT**2,MG,MT2)))))
431 MTMT=MT/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/MT))*(ASMT/PI)**2
432 $ -ASMT/3./PI*(REAL8(SSB1(MT**2,MG,MT2))+
433 $ REAL8(SSB1(MT**2,MG,MT1))+SIN(2*THETAT)*MG/MT*
434 $ (REAL8(SSB0(MT**2,MG,MT2))-REAL8(SSB0(MT**2,MG,MT1)))))
438 IF (ABS(COS(THETAB)).LT..707107) THEN
439 MBMZ=MBMZ*(1.+ASMZ/3./PI*(REAL8(SSB1(MZ**2,MG,MB1))+
440 $ REAL8(SSB1(MZ**2,MG,MB2))-SIN(2*THETAB)*MG/MB*
441 $ (REAL8(SSB0(MZ**2,MG,MB1))-REAL8(SSB0(MZ**2,MG,MB2))))
442 $ -FTMT**2*MU*(-AAT*TANB+MU)/16./PI**2/(MT1**2-MT2**2)*
443 $ (REAL8(SSB0(MZ**2,AMU,MT1))-REAL8(SSB0(MZ**2,AMU,MT2)))+
444 $ G2**2*MU*M2*TANB/16./PI**2/(AMU**2-M2**2)*
445 $ (SIN(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT1))-
446 $ REAL8(SSB0(MZ**2,AMU,MT1)))+
447 $ COS(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT2))-
448 $ REAL8(SSB0(MZ**2,AMU,MT2)))))
450 MBMZ=MBMZ*(1.+ASMZ/3./PI*(REAL8(SSB1(MZ**2,MG,MB2))+
451 $ REAL8(SSB1(MZ**2,MG,MB1))+SIN(2*THETAB)*MG/MB*
452 $ (REAL8(SSB0(MZ**2,MG,MB2))-REAL8(SSB0(MZ**2,MG,MB1))))
453 $ -FTMT**2*MU*(-AAT*TANB+MU)/16./PI**2/(MT2**2-MT1**2)*
454 $ (REAL8(SSB0(MZ**2,AMU,MT2))-REAL8(SSB0(MZ**2,AMU,MT1)))+
455 $ G2**2*MU*M2*TANB/16./PI**2/(AMU**2-M2**2)*
456 $ (COS(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT2))-
457 $ REAL8(SSB0(MZ**2,AMU,MT2)))+
458 $ SIN(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT1))-
459 $ REAL8(SSB0(MZ**2,AMU,MT1)))))
462 MTAMZ=MTAMZ*(1.+G2**2*MU*M2*TANB/16./PI**2/(MUS-M2**2)*
463 $(REAL8(SSB0(MZ**2,AM2,MSN))-REAL8(SSB0(MZ**2,AMU,MSN))))
466 C Iterate entire process, increasing NSTEP each time
467 C This time, freeze out parameters at sqrt(t_l t_r)
469 HIGFRZ=MAX(AMZ,(G0(23)*G0(24))**0.25)
474 CALL SUGRGE(M0,MHF,A0,TANB,SGNMU,MT,G,G0,IG,W2,NSTEP,IMODEL)
475 IF(NOGOOD.NE.0) GO TO 100
478 320 DELG0=MAX(DELG0,ABS((G0(J)-G0SAVE(J))/G0(J)))
479 IF(DELG0.LT.DELLIM) GO TO 400
481 WRITE(LOUT,1000) MXITER
482 1000 FORMAT(/' SUGRA WARNING: NO CONVERGENCE IN',I4,' ITERATIONS')
493 C Fill XISAIN common block
499 XISAIN(5)=SQRT(G0(19))
500 XISAIN(6)=SQRT(G0(17))
501 XISAIN(7)=SQRT(G0(18))
502 XISAIN(8)=SQRT(G0(16))
503 XISAIN(9)=SQRT(G0(15))
509 XISAIN(15)=SQRT(G0(24))
510 XISAIN(16)=SQRT(G0(22))
511 XISAIN(17)=SQRT(G0(23))
512 XISAIN(18)=SQRT(G0(21))
513 XISAIN(19)=SQRT(G0(20))