Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyinom.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyinom.f
1 
2 C*********************************************************************
3 
4 C...PYINOM
5 C...Finds the mass eigenstates and mixing matrices for neutralinos
6 C...and charginos.
7 
8  SUBROUTINE pyinom
9 
10 C...Double precision and integer declarations.
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13  INTEGER pycomp
14 C...Parameter statement to help give large particle numbers.
15  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
16  &kexcit=4000000,kdimen=5000000)
17 C...Commonblocks.
18  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
19  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
20  common/pymssm/imss(0:99),rmss(0:99)
21  common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
22  &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
23  SAVE /pydat1/,/pydat2/,/pymssm/,/pyssmt/
24 
25 C...Local variables.
26  DOUBLE PRECISION xmw,xmz,xm(4)
27  DOUBLE PRECISION ar(4,4),wr(4),zr(4,4),zi(4,4),ai(4,4)
28  DOUBLE PRECISION wi(4),fv1(4),fv2(4),fv3(4)
29  DOUBLE PRECISION cosw,sinw
30  DOUBLE PRECISION xmu
31  DOUBLE PRECISION tanb,cosb,sinb
32  DOUBLE PRECISION xm1,xm2,xm3,beta
33  DOUBLE PRECISION q2,aem,a1,a2,aq,rm1,rm2
34  DOUBLE PRECISION arg,x0,x1,ax0,ax1,at,bt
35  DOUBLE PRECISION y0,y1,amgx0,am1x0,amgx1,am1x1
36  DOUBLE PRECISION argx0,ar1x0,argx1,ar1x1
37  DOUBLE PRECISION pyalps,pyalem
38  DOUBLE PRECISION pyrnm3
39  COMPLEX*16 car(4,4),cai(4,4),ca1,ca2
40  INTEGER ierr,index(4),i,j,k,iopt,ilr,kfnchi(4)
41  DATA kfnchi/1000022,1000023,1000025,1000035/
42 
43  iopt=imss(2)
44  IF(imss(1).EQ.2) THEN
45  iopt=1
46  ENDIF
47 C...M1, M2, AND M3 ARE INDEPENDENT
48  IF(iopt.EQ.0) THEN
49  xm1=rmss(1)
50  xm2=rmss(2)
51  xm3=rmss(3)
52  ELSEIF(iopt.GE.1) THEN
53  q2=pmas(23,1)**2
54  aem=pyalem(q2)
55  a2=aem/paru(102)
56  a1=aem/(1d0-paru(102))
57  xm1=rmss(1)
58  xm2=rmss(2)
59  IF(imss(1).EQ.2) xm1=rmss(1)/rmss(20)*a1*5d0/3d0
60  IF(iopt.EQ.1) THEN
61  xm2=xm1*a2/a1*3d0/5d0
62  rmss(2)=xm2
63  ELSEIF(iopt.EQ.3) THEN
64  xm1=xm2*5d0/3d0*a1/a2
65  rmss(1)=xm1
66  ENDIF
67  xm3=pyrnm3(xm2/a2)
68  rmss(3)=xm3
69  IF(xm3.LE.0d0) THEN
70  WRITE(mstu(11),*) ' ERROR WITH M3 = ',xm3
71  CALL pystop(105)
72  ENDIF
73  ENDIF
74 
75 C...GLUINO MASS
76  IF(imss(3).EQ.1) THEN
77  pmas(pycomp(ksusy1+21),1)=abs(xm3)
78  ELSE
79  aq=0d0
80  DO 110 i=1,4
81  DO 100 ilr=1,2
82  rm1=pmas(pycomp(ilr*ksusy1+i),1)**2/xm3**2
83  aq=aq+0.5d0*((2d0-rm1)*(rm1*log(rm1)-1d0)
84  & +(1d0-rm1)**2*log(abs(1d0-rm1)))
85  100 CONTINUE
86  110 CONTINUE
87 
88  DO 130 i=5,6
89  DO 120 ilr=1,2
90  rm1=pmas(pycomp(ilr*ksusy1+i),1)**2/xm3**2
91  rm2=pmas(i,1)**2/xm3**2
92  arg=(rm1-rm2-1d0)**2-4d0*rm2**2
93  IF(arg.GE.0d0) THEN
94  x0=0.5d0*(1d0+rm2-rm1-sqrt(arg))
95  ax0=abs(x0)
96  x1=0.5d0*(1d0+rm2-rm1+sqrt(arg))
97  ax1=abs(x1)
98  IF(x0.EQ.1d0) THEN
99  at=-1d0
100  bt=0.25d0
101  ELSEIF(x0.EQ.0d0) THEN
102  at=0d0
103  bt=-0.25d0
104  ELSE
105  at=0.5d0*log(abs(1d0-x0))*(1d0-x0**2)+
106  & 0.5d0*x0**2*log(ax0)
107  bt=(-1d0-2d0*x0)/4d0
108  ENDIF
109  IF(x1.EQ.1d0) THEN
110  at=-1d0+at
111  bt=0.25d0+bt
112  ELSEIF(x1.EQ.0d0) THEN
113  at=0d0+at
114  bt=-0.25d0+bt
115  ELSE
116  at=0.5d0*log(abs(1d0-x1))*(1d0-x1**2)+0.5d0*
117  & x1**2*log(ax1)+at
118  bt=(-1d0-2d0*x1)/4d0+bt
119  ENDIF
120  aq=aq+at+bt
121  ELSE
122  x0=0.5d0*(1d0+rm2-rm1)
123  y0=-0.5d0*sqrt(-arg)
124  amgx0=sqrt(x0**2+y0**2)
125  am1x0=sqrt((1d0-x0)**2+y0**2)
126  argx0=atan2(-x0,-y0)
127  ar1x0=atan2(1d0-x0,y0)
128  x1=x0
129  y1=-y0
130  amgx1=amgx0
131  am1x1=am1x0
132  argx1=atan2(-x1,-y1)
133  ar1x1=atan2(1d0-x1,y1)
134  at=0.5d0*log(am1x0)*(1d0-x0**2+3d0*y0**2)
135  & +0.5d0*(x0**2-y0**2)*log(amgx0)
136  bt=(-1d0-2d0*x0)/4d0+x0*y0*( ar1x0-argx0 )
137  at=at+0.5d0*log(am1x1)*(1d0-x1**2+3d0*y1**2)
138  & +0.5d0*(x1**2-y1**2)*log(amgx1)
139  bt=bt+(-1d0-2d0*x1)/4d0+x1*y1*( ar1x1-argx1 )
140  aq=aq+at+bt
141  ENDIF
142  120 CONTINUE
143  130 CONTINUE
144  pmas(pycomp(ksusy1+21),1)=abs(xm3)*(1d0+pyalps(xm3**2)
145  & /(2d0*paru(2))*(15d0+aq))
146  ENDIF
147 
148 C...NEUTRALINO MASSES
149  DO 150 i=1,4
150  DO 140 j=1,4
151  ai(i,j)=0d0
152  140 CONTINUE
153  150 CONTINUE
154  xmz=pmas(23,1)
155  xmw=pmas(24,1)
156  xmu=rmss(4)
157  sinw=sqrt(paru(102))
158  cosw=sqrt(1d0-paru(102))
159  tanb=rmss(5)
160  beta=atan(tanb)
161  cosb=cos(beta)
162  sinb=tanb*cosb
163 
164 C... Definitions:
165 C... psi^0 =(-i bino^0, -i wino^0, h_d^0(=H_1^0), h_u^0(=H_2^0))
166 C... => L_neutralino = -1/2*(psi^0)^T * [AR] * psi^0 + h.c.
167  ar(1,1) = xm1*cos(rmss(30))
168  ai(1,1) = xm1*sin(rmss(30))
169  ar(2,2) = xm2*cos(rmss(31))
170  ai(2,2) = xm2*sin(rmss(31))
171  ar(3,3) = 0d0
172  ar(4,4) = 0d0
173  ar(1,2) = 0d0
174  ar(2,1) = 0d0
175  ar(1,3) = -xmz*sinw*cosb
176  ar(3,1) = ar(1,3)
177  ar(1,4) = xmz*sinw*sinb
178  ar(4,1) = ar(1,4)
179  ar(2,3) = xmz*cosw*cosb
180  ar(3,2) = ar(2,3)
181  ar(2,4) = -xmz*cosw*sinb
182  ar(4,2) = ar(2,4)
183  ar(3,4) = -xmu*cos(rmss(33))
184  ai(3,4) = -xmu*sin(rmss(33))
185  ar(4,3) = -xmu*cos(rmss(33))
186  ai(4,3) = -xmu*sin(rmss(33))
187 C CALL PYEIG4(AR,WR,ZR)
188  CALL pyeicg(4,4,ar,ai,wr,wi,1,zr,zi,fv1,fv2,fv3,ierr)
189  IF(ierr.NE.0) THEN
190  WRITE(mstu(11),*) ' PROBLEM WITH PYEICG IN PYINOM '
191  ENDIF
192  DO 160 i=1,4
193  index(i)=i
194  xm(i)=abs(wr(i))
195  160 CONTINUE
196  DO 180 i=2,4
197  k=i
198  DO 170 j=i-1,1,-1
199  IF(xm(k).LT.xm(j)) THEN
200  itmp=index(j)
201  xtmp=xm(j)
202  index(j)=index(k)
203  xm(j)=xm(k)
204  index(k)=itmp
205  xm(k)=xtmp
206  k=k-1
207  ELSE
208  goto 180
209  ENDIF
210  170 CONTINUE
211  180 CONTINUE
212 
213 
214  DO 210 i=1,4
215  k=index(i)
216  smz(i)=wr(k)
217  pmas(pycomp(kfnchi(i)),1)=abs(smz(i))
218  s=0d0
219  DO 190 j=1,4
220  s=s+zr(j,k)**2+zi(j,k)**2
221  190 CONTINUE
222  DO 200 j=1,4
223  zmix(i,j)=zr(j,k)/sqrt(s)
224  zmixi(i,j)=zi(j,k)/sqrt(s)
225  IF(abs(zmix(i,j)).LT.1d-6) zmix(i,j)=0d0
226  IF(abs(zmixi(i,j)).LT.1d-6) zmixi(i,j)=0d0
227  200 CONTINUE
228  210 CONTINUE
229 
230 C...CHARGINO MASSES
231 C.....Find eigenvectors of X X^*
232  ai(1,1) = 0d0
233  ai(2,2) = 0d0
234  ar(1,1) = xm2**2+2d0*xmw**2*sinb**2
235  ar(2,2) = xmu**2+2d0*xmw**2*cosb**2
236  ar(1,2) = sqrt(2d0)*xmw*(xm2*cos(rmss(31))*cosb+
237  &xmu*cos(rmss(33))*sinb)
238  ai(1,2) = sqrt(2d0)*xmw*(xm2*sin(rmss(31))*cosb-
239  &xmu*sin(rmss(33))*sinb)
240  ar(2,1) = sqrt(2d0)*xmw*(xm2*cos(rmss(31))*cosb+
241  &xmu*cos(rmss(33))*sinb)
242  ai(2,1) = sqrt(2d0)*xmw*(-xm2*sin(rmss(31))*cosb+
243  &xmu*sin(rmss(33))*sinb)
244  CALL pyeicg(4,2,ar,ai,wr,wi,1,zr,zi,fv1,fv2,fv3,ierr)
245  IF(ierr.NE.0) THEN
246  WRITE(mstu(11),*) ' PROBLEM WITH PYEICG IN PYINOM '
247  ENDIF
248  index(1)=1
249  index(2)=2
250  IF(wr(2).LT.wr(1)) THEN
251  index(1)=2
252  index(2)=1
253  ENDIF
254 
255  DO 240 i=1,2
256  k=index(i)
257  smw(i)=sqrt(wr(k))
258  s=0d0
259  DO 220 j=1,2
260  s=s+zr(j,k)**2+zi(j,k)**2
261  220 CONTINUE
262  DO 230 j=1,2
263  umix(i,j)=zr(j,k)/sqrt(s)
264  umixi(i,j)=-zi(j,k)/sqrt(s)
265  IF(abs(umix(i,j)).LT.1d-6) umix(i,j)=0d0
266  IF(abs(umixi(i,j)).LT.1d-6) umixi(i,j)=0d0
267  230 CONTINUE
268  240 CONTINUE
269 C...Force chargino mass > neutralino mass
270  IF(abs(smw(1)).LT.abs(smz(1))+2d0*pmas(pycomp(111),1)) THEN
271  CALL pyerrm(18,'(PYINOM:) '//
272  & 'forcing m(~chi+_1) > m(~chi0_1) + 2m(pi0)')
273  smw(1)=sign(abs(smz(1))+2d0*pmas(pycomp(111),1),smw(1))
274  ENDIF
275  pmas(pycomp(ksusy1+24),1)=smw(1)
276  pmas(pycomp(ksusy1+37),1)=smw(2)
277 
278 C.....Find eigenvectors of X^* X
279  ai(1,1) = 0d0
280  ai(2,2) = 0d0
281  ar(1,1) = xm2**2+2d0*xmw**2*cosb**2
282  ar(2,2) = xmu**2+2d0*xmw**2*sinb**2
283  ar(1,2) = sqrt(2d0)*xmw*(xm2*cos(rmss(31))*sinb+
284  &xmu*cos(rmss(33))*cosb)
285  ai(1,2) = sqrt(2d0)*xmw*(-xm2*sin(rmss(31))*sinb+
286  &xmu*sin(rmss(33))*cosb)
287  ar(2,1) = sqrt(2d0)*xmw*(xm2*cos(rmss(31))*sinb+
288  &xmu*cos(rmss(33))*cosb)
289  ai(2,1) = sqrt(2d0)*xmw*(xm2*sin(rmss(31))*sinb-
290  &xmu*sin(rmss(33))*cosb)
291  CALL pyeicg(4,2,ar,ai,wr,wi,1,zr,zi,fv1,fv2,fv3,ierr)
292  IF(ierr.NE.0) THEN
293  WRITE(mstu(11),*) ' PROBLEM WITH PYEICG IN PYINOM '
294  ENDIF
295  index(1)=1
296  index(2)=2
297  IF(wr(2).LT.wr(1)) THEN
298  index(1)=2
299  index(2)=1
300  ENDIF
301 
302  DO 270 i=1,2
303  k=index(i)
304  s=0d0
305  DO 250 j=1,2
306  s=s+zr(j,k)**2+zi(j,k)**2
307  250 CONTINUE
308  DO 260 j=1,2
309  vmix(i,j)=zr(j,k)/sqrt(s)
310  vmixi(i,j)=-zi(j,k)/sqrt(s)
311  IF(abs(vmix(i,j)).LT.1d-6) vmix(i,j)=0d0
312  IF(abs(vmixi(i,j)).LT.1d-6) vmixi(i,j)=0d0
313  260 CONTINUE
314  270 CONTINUE
315 
316 
317  RETURN
318  END