Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pypole.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pypole.f
1 
2 C*********************************************************************
3 
4 C...PYPOLE
5 C...This subroutine computes the CP-even higgs and CP-odd pole
6 c...Higgs masses and mixing angles.
7 
8 C...Program based on the work by M. Carena, M. Quiros
9 C...and C.E.M. Wagner, "Effective potential methods and
10 C...the Higgs mass spectrum in the MSSM", CERN-TH/95-157
11 
12 C...Inputs: IHIGGS(explained below),MCHI,MA,TANB,MQ,MUR,MDR,MTOP,
13 C...AT,AB,MU
14 C...where MCHI is the largest chargino mass, MA is the running
15 C...CP-odd higgs mass, TANB is the value of the ratio of vacuum
16 C...expectaion values at the scale MTOP, MQ is the third generation
17 C...left handed squark mass parameter, MUR is the third generation
18 C...right handed stop mass parameter, MDR is the third generation
19 C...right handed sbottom mass parameter, MTOP is the pole top quark
20 C...mass; AT,AB are the soft supersymmetry breaking trilinear
21 C...couplings of the stop and sbottoms, respectively, and MU is the
22 C...supersymmetric mass parameter
23 
24 C...The parameter IHIGGS=0,1,2,3 corresponds to the number of
25 C...Higgses whose pole mass is computed. If IHIGGS=0 only running
26 C...masses are given, what makes the running of the program
27 c...much faster and it is quite generally a good approximation
28 c...(for a theoretical discussion see ref. above). If IHIGGS=1,
29 C...only the pole mass for H is computed. If IHIGGS=2, then h and H,
30 c...and if IHIGGS=3, then h,H,A polarizations are computed
31 
32 C...Output: MH and MHP which are the lightest CP-even Higgs running
33 C...and pole masses, respectively; HM and HMP are the heaviest CP-even
34 C...Higgs running and pole masses, repectively; SA and CA are the
35 C...SIN(ALPHA) and COS(ALPHA) where ALPHA is the Higgs mixing angle
36 C...AMP is the CP-odd Higgs pole mass. STOP1,STOP2,SBOT1 and SBOT2
37 C...are the stop and sbottom mass eigenvalues. Finally, TANBA is
38 C...the value of TANB at the CP-odd Higgs mass scale
39 
40 C...This subroutine makes use of CERN library subroutine
41 C...integration package, which makes the computation of the
42 C...pole Higgs masses somewhat faster. We thank P. Janot for this
43 C...improvement. Those who are not able to call the CERN
44 C...libraries, please use the subroutine SUBHPOLE2.F, which
45 C...although somewhat slower, gives identical results
46 
47  SUBROUTINE pypole(IHIGGS,XMC,XMA,TANB,XMQ,XMUR,XMDR,XMT,AT,AB,XMU,
48  &xmh,xmhp,hm,hmp,amp,sa,ca,stop1,stop2,sbot1,sbot2,tanba,xmg,dt,db)
49 
50 C...Double precision and integer declarations.
51  IMPLICIT DOUBLE PRECISION(a-h, o-z)
52  IMPLICIT INTEGER(i-n)
53 
54 C...Parameters.
55  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
56  SAVE /pydat1/
57  INTEGER pyk,pychge,pycomp
58 
59 C...Local variables.
60  dimension delta(2,2),coupt(2,2),t(2,2),sstop2(2),
61  &ssbot2(2),b(2,2),coupb(2,2),
62  &hcoupt(2,2),hcoupb(2,2),
63  &acoupt(2,2),acoupb(2,2),pr(3), polar(3)
64 
65  delta(1,1) = 1d0
66  delta(2,2) = 1d0
67  delta(1,2) = 0d0
68  delta(2,1) = 0d0
69  v = 174.1d0
70  xmz=91.18d0
71  pi=paru(1)
72  rxmt=pymrun(6,xmt**2)
73  CALL pyrghm(xmc,xma,tanb,xmq,xmur,xmdr,xmt,at,ab,
74  &xmu,xmh,hm,xmch,sa,ca,sab,cab,tanba,xmg,dt,db)
75 
76  sinb = tanb/(tanb**2+1d0)**0.5d0
77  cosb = 1d0/(tanb**2+1d0)**0.5d0
78  cos2b = sinb**2 - cosb**2
79  sinbpa = sinb*ca + cosb*sa
80  cosbpa = cosb*ca - sinb*sa
81  rmbot = pymrun(5,xmt**2)
82  xmq2 = xmq**2
83  xmur2 = xmur**2
84  IF(xmur.LT.0d0) xmur2=-xmur2
85  xmdr2 = xmdr**2
86  xmst11 = rxmt**2 + xmq2 - 0.35d0*xmz**2*cos2b
87  xmst22 = rxmt**2 + xmur2 - 0.15d0*xmz**2*cos2b
88  IF(xmst11.LT.0d0) goto 500
89  IF(xmst22.LT.0d0) goto 500
90  xmsb11 = rmbot**2 + xmq2 + 0.42d0*xmz**2*cos2b
91  xmsb22 = rmbot**2 + xmdr2 + 0.08d0*xmz**2*cos2b
92  IF(xmsb11.LT.0d0) goto 500
93  IF(xmsb22.LT.0d0) goto 500
94 C WMST11 = RXMT**2 + XMQ2
95 C WMST22 = RXMT**2 + XMUR2
96  xmst12 = rxmt*(at - xmu/tanb)
97  xmsb12 = rmbot*(ab - xmu*tanb)
98 
99 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
100 C...STOP EIGENVALUES CALCULATION
101 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
102 
103  stop12 = 0.5d0*(xmst11+xmst22) +
104  &0.5d0*((xmst11+xmst22)**2 -
105  &4d0*(xmst11*xmst22 - xmst12**2))**0.5d0
106  stop22 = 0.5d0*(xmst11+xmst22) -
107  &0.5d0*((xmst11+xmst22)**2 - 4d0*(xmst11*xmst22 -
108  &xmst12**2))**0.5d0
109 
110  IF(stop22.LT.0d0) goto 500
111  sstop2(1) = stop12
112  sstop2(2) = stop22
113  stop1 = stop12**0.5d0
114  stop2 = stop22**0.5d0
115 C STOP1W = STOP1
116 C STOP2W = STOP2
117 
118  IF(xmst12.EQ.0d0) xst11 = 1d0
119  IF(xmst12.EQ.0d0) xst12 = 0d0
120  IF(xmst12.EQ.0d0) xst21 = 0d0
121  IF(xmst12.EQ.0d0) xst22 = 1d0
122 
123  IF(xmst12.EQ.0d0) goto 110
124 
125  100 xst11 = xmst12/(xmst12**2+(xmst11-stop12)**2)**0.5d0
126  xst12 = - (xmst11-stop12)/(xmst12**2+(xmst11-stop12)**2)**0.5d0
127  xst21 = xmst12/(xmst12**2+(xmst11-stop22)**2)**0.5d0
128  xst22 = - (xmst11-stop22)/(xmst12**2+(xmst11-stop22)**2)**0.5d0
129 
130  110 t(1,1) = xst11
131  t(2,2) = xst22
132  t(1,2) = xst12
133  t(2,1) = xst21
134 
135  sbot12 = 0.5d0*(xmsb11+xmsb22) +
136  &0.5d0*((xmsb11+xmsb22)**2 -
137  &4d0*(xmsb11*xmsb22 - xmsb12**2))**0.5d0
138  sbot22 = 0.5d0*(xmsb11+xmsb22) -
139  &0.5d0*((xmsb11+xmsb22)**2 - 4d0*(xmsb11*xmsb22 -
140  &xmsb12**2))**0.5d0
141  IF(sbot22.LT.0d0) goto 500
142  sbot1 = sbot12**0.5d0
143  sbot2 = sbot22**0.5d0
144 
145  ssbot2(1) = sbot12
146  ssbot2(2) = sbot22
147 
148  IF(xmsb12.EQ.0d0) xsb11 = 1d0
149  IF(xmsb12.EQ.0d0) xsb12 = 0d0
150  IF(xmsb12.EQ.0d0) xsb21 = 0d0
151  IF(xmsb12.EQ.0d0) xsb22 = 1d0
152 
153  IF(xmsb12.EQ.0d0) goto 130
154 
155  120 xsb11 = xmsb12/(xmsb12**2+(xmsb11-sbot12)**2)**0.5d0
156  xsb12 = - (xmsb11-sbot12)/(xmsb12**2+(xmsb11-sbot12)**2)**0.5d0
157  xsb21 = xmsb12/(xmsb12**2+(xmsb11-sbot22)**2)**0.5d0
158  xsb22 = - (xmsb11-sbot22)/(xmsb12**2+(xmsb11-sbot22)**2)**0.5d0
159 
160  130 b(1,1) = xsb11
161  b(2,2) = xsb22
162  b(1,2) = xsb12
163  b(2,1) = xsb21
164 
165 
166  sint = 0.2320d0
167  sqr = dsqrt(2d0)
168  vp = 174.1d0*sqr
169 
170 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
171 C...STARTING OF LIGHT HIGGS
172 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
173 
174  IF(ihiggs.EQ.0) goto 490
175 
176  DO 150 i = 1,2
177  DO 140 j = 1,2
178  coupt(i,j) =
179  & sint*xmz**2*2d0*sqr/174.1d0/3d0*sinbpa*(delta(i,j) +
180  & (3d0 - 8d0*sint)/4d0/sint*t(1,i)*t(1,j))
181  & -rxmt**2/174.1d0**2*vp/sinb*ca*delta(i,j)
182  & -rxmt/vp/sinb*(at*ca + xmu*sa)*(t(1,i)*t(2,j) +
183  & t(1,j)*t(2,i))
184  140 CONTINUE
185  150 CONTINUE
186 
187 
188  DO 170 i = 1,2
189  DO 160 j = 1,2
190  coupb(i,j) =
191  & -sint*xmz**2*2d0*sqr/174.1d0/6d0*sinbpa*(delta(i,j) +
192  & (3d0 - 4d0*sint)/2d0/sint*b(1,i)*b(1,j))
193  & +rmbot**2/174.1d0**2*vp/cosb*sa*delta(i,j)
194  & +rmbot/vp/cosb*(ab*sa + xmu*ca)*(b(1,i)*b(2,j) +
195  & b(1,j)*b(2,i))
196  160 CONTINUE
197  170 CONTINUE
198 
199  prun = xmh
200  eps = 1d-4*prun
201  iter = 0
202  180 iter = iter + 1
203  DO 230 i3 = 1,3
204 
205  pr(i3)=prun+(i3-2)*eps/2
206  p2=pr(i3)**2
207  polt = 0d0
208  DO 200 i = 1,2
209  DO 190 j = 1,2
210  polt = polt + coupt(i,j)**2*3d0*
211  & pyfint(p2,sstop2(i),sstop2(j))/16d0/pi**2
212  190 CONTINUE
213  200 CONTINUE
214 
215  polb = 0d0
216  DO 220 i = 1,2
217  DO 210 j = 1,2
218  polb = polb + coupb(i,j)**2*3d0*
219  & pyfint(p2,ssbot2(i),ssbot2(j))/16d0/pi**2
220  210 CONTINUE
221  220 CONTINUE
222 C RXMT2 = RXMT**2
223  xmt2=xmt**2
224 
225  poltt =
226  & 3d0*rxmt**2/8d0/pi**2/ v **2*
227  & ca**2/sinb**2 *
228  & (-2d0*xmt**2+0.5d0*p2)*
229  & pyfint(p2,xmt2,xmt2)
230 
231  pol = polt + polb + poltt
232  polar(i3) = p2 - xmh**2 - pol
233  230 CONTINUE
234  deriv = (polar(3)-polar(1))/eps
235  drun = - polar(2)/deriv
236  prun = prun + drun
237  p2 = prun**2
238  IF( abs(drun) .LT. 1d-4 .OR.iter.GT.500) goto 240
239  goto 180
240  240 CONTINUE
241 
242  xmhp = dsqrt(p2)
243 
244 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
245 C...END OF LIGHT HIGGS
246 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
247 
248  250 IF(ihiggs.EQ.1) goto 490
249 
250 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
251 C... STARTING OF HEAVY HIGGS
252 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
253 
254  DO 270 i = 1,2
255  DO 260 j = 1,2
256  hcoupt(i,j) =
257  & -sint*xmz**2*2d0*sqr/174.1d0/3d0*cosbpa*(delta(i,j) +
258  & (3d0 - 8d0*sint)/4d0/sint*t(1,i)*t(1,j))
259  & -rxmt**2/174.1d0**2*vp/sinb*sa*delta(i,j)
260  & -rxmt/vp/sinb*(at*sa - xmu*ca)*(t(1,i)*t(2,j) +
261  & t(1,j)*t(2,i))
262  260 CONTINUE
263  270 CONTINUE
264 
265  DO 290 i = 1,2
266  DO 280 j = 1,2
267  hcoupb(i,j) =
268  & sint*xmz**2*2d0*sqr/174.1d0/6d0*cosbpa*(delta(i,j) +
269  & (3d0 - 4d0*sint)/2d0/sint*b(1,i)*b(1,j))
270  & -rmbot**2/174.1d0**2*vp/cosb*ca*delta(i,j)
271  & -rmbot/vp/cosb*(ab*ca - xmu*sa)*(b(1,i)*b(2,j) +
272  & b(1,j)*b(2,i))
273  hcoupb(i,j)=0d0
274  280 CONTINUE
275  290 CONTINUE
276 
277  prun = hm
278  eps = 1d-4*prun
279  iter = 0
280  300 iter = iter + 1
281  DO 350 i3 = 1,3
282  pr(i3)=prun+(i3-2)*eps/2
283  hp2=pr(i3)**2
284 
285  hpolt = 0d0
286  DO 320 i = 1,2
287  DO 310 j = 1,2
288  hpolt = hpolt + hcoupt(i,j)**2*3d0*
289  & pyfint(hp2,sstop2(i),sstop2(j))/16d0/pi**2
290  310 CONTINUE
291  320 CONTINUE
292 
293  hpolb = 0d0
294  DO 340 i = 1,2
295  DO 330 j = 1,2
296  hpolb = hpolb + hcoupb(i,j)**2*3d0*
297  & pyfint(hp2,ssbot2(i),ssbot2(j))/16d0/pi**2
298  330 CONTINUE
299  340 CONTINUE
300 
301 C RXMT2 = RXMT**2
302  xmt2 = xmt**2
303 
304  hpoltt =
305  & 3d0*rxmt**2/8d0/pi**2/ v **2*
306  & sa**2/sinb**2 *
307  & (-2d0*xmt**2+0.5d0*hp2)*
308  & pyfint(hp2,xmt2,xmt2)
309 
310  hpol = hpolt + hpolb + hpoltt
311  polar(i3) =hp2-hm**2-hpol
312  350 CONTINUE
313  deriv = (polar(3)-polar(1))/eps
314  drun = - polar(2)/deriv
315  prun = prun + drun
316  hp2 = prun**2
317  IF( abs(drun) .LT. 1d-4 .OR.iter.GT.500) goto 360
318  goto 300
319  360 CONTINUE
320 
321 
322  370 CONTINUE
323  hmp = hp2**0.5d0
324 
325 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
326 C... END OF HEAVY HIGGS
327 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
328 
329  IF(ihiggs.EQ.2) goto 490
330 
331 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
332 C...BEGINNING OF PSEUDOSCALAR HIGGS
333 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
334 
335  DO 390 i = 1,2
336  DO 380 j = 1,2
337  acoupt(i,j) =
338  & -rxmt/vp/sinb*(at*cosb + xmu*sinb)*
339  & (t(1,i)*t(2,j) -t(1,j)*t(2,i))
340  380 CONTINUE
341  390 CONTINUE
342  DO 410 i = 1,2
343  DO 400 j = 1,2
344  acoupb(i,j) =
345  & rmbot/vp/cosb*(ab*sinb + xmu*cosb)*
346  & (b(1,i)*b(2,j) -b(1,j)*b(2,i))
347  400 CONTINUE
348  410 CONTINUE
349 
350  prun = xma
351  eps = 1d-4*prun
352  iter = 0
353  420 iter = iter + 1
354  DO 470 i3 = 1,3
355  pr(i3)=prun+(i3-2)*eps/2
356  ap2=pr(i3)**2
357  apolt = 0d0
358  DO 440 i = 1,2
359  DO 430 j = 1,2
360  apolt = apolt + acoupt(i,j)**2*3d0*
361  & pyfint(ap2,sstop2(i),sstop2(j))/16d0/pi**2
362  430 CONTINUE
363  440 CONTINUE
364  apolb = 0d0
365  DO 460 i = 1,2
366  DO 450 j = 1,2
367  apolb = apolb + acoupb(i,j)**2*3d0*
368  & pyfint(ap2,ssbot2(i),ssbot2(j))/16d0/pi**2
369  450 CONTINUE
370  460 CONTINUE
371 C RXMT2 = RXMT**2
372  xmt2=xmt**2
373  apoltt =
374  & 3d0*rxmt**2/8d0/pi**2/ v **2*
375  & cosb**2/sinb**2 *
376  & (-0.5d0*ap2)*
377  & pyfint(ap2,xmt2,xmt2)
378  apol = apolt + apolb + apoltt
379  polar(i3) = ap2 - xma**2 -apol
380  470 CONTINUE
381  deriv = (polar(3)-polar(1))/eps
382  drun = - polar(2)/deriv
383  prun = prun + drun
384  ap2 = prun**2
385  IF( abs(drun) .LT. 1d-4 .OR.iter.GT.500) goto 480
386  goto 420
387  480 CONTINUE
388 
389  amp = dsqrt(ap2)
390 
391 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
392 C...END OF PSEUDOSCALAR HIGGS
393 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
394 
395  IF(ihiggs.EQ.3) goto 490
396 
397  490 CONTINUE
398  RETURN
399  500 CONTINUE
400  WRITE(mstu(11),*) ' EXITING IN PYPOLE '
401  WRITE(mstu(11),*) ' XMST11,XMST22 = ',xmst11,xmst22
402  WRITE(mstu(11),*) ' XMSB11,XMSB22 = ',xmsb11,xmsb22
403  WRITE(mstu(11),*) ' STOP22,SBOT22 = ',stop22,sbot22
404  CALL pystop(107)
405  END