Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pymsin.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pymsin.f
1 
2 C*********************************************************************
3 
4 C...PYMSIN
5 C...Initializes supersymmetry: finds sparticle masses and
6 C...branching ratios and stores this information.
7 C...AUTHOR: STEPHEN MRENNA
8 C...Author: P. Skands (SLHA + RPV + ISASUSY Interface, NMSSM)
9 
10  SUBROUTINE pymsin
11 
12 C...Double precision and integer declarations.
13  IMPLICIT DOUBLE PRECISION(a-h, o-z)
14  IMPLICIT INTEGER(i-n)
15  INTEGER pyk,pychge,pycomp
16 C...Parameter statement to help give large particle numbers.
17  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
18  &kexcit=4000000,kdimen=5000000)
19 C...Commonblocks.
20  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
21  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
22  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
23  common/pydat4/chaf(500,2)
24  CHARACTER chaf*16
25  common/pypars/mstp(200),parp(200),msti(200),pari(200)
26  common/pyint4/mwid(500),wids(500,5)
27  common/pymssm/imss(0:99),rmss(0:99)
28  common/pymsrv/rvlam(3,3,3), rvlamp(3,3,3), rvlamb(3,3,3)
29  common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
30  &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
31  common/pyhtri/hhh(7)
32  SAVE /pydat1/,/pydat2/,/pydat3/,/pydat4/,/pypars/,/pyint4/,
33  &/pymssm/,/pymsrv/,/pyssmt/
34 
35 C...Local variables.
36  DOUBLE PRECISION alfa,beta
37  DOUBLE PRECISION tanb,al,be,cosa,cosb,sina,sinb,xw
38  INTEGER i,j,j1,i1,k1
39  INTEGER kc,lknt,idlam(400,3)
40  DOUBLE PRECISION xlam(0:400)
41  DOUBLE PRECISION wdtp(0:400),wdte(0:400,0:5)
42  DOUBLE PRECISION xarg,cos2b,xmw2,xmz2
43  DOUBLE PRECISION delm,xmdif
44  DOUBLE PRECISION dx,dy,ds,dmu2,dma2,dq2,du2,dd2,dl2,de2,dhu2,dhd2
45  DOUBLE PRECISION arg,sgnmu,r
46  INTEGER imssm
47  INTEGER irprty
48  INTEGER kfsusy(50),mwidsu(36),mdcysu(36)
49  SAVE mwidsu,mdcysu
50  DATA kfsusy/
51  &1000001,2000001,1000002,2000002,1000003,2000003,
52  &1000004,2000004,1000005,2000005,1000006,2000006,
53  &1000011,2000011,1000012,2000012,1000013,2000013,
54  &1000014,2000014,1000015,2000015,1000016,2000016,
55  &1000021,1000022,1000023,1000025,1000035,1000024,
56  &1000037,1000039, 25, 35, 36, 37,
57  & 6, 24, 45, 46,1000045, 9*0/
58  DATA init/0/
59 
60 C...Do nothing if SUSY not requested.
61  imssm=imss(1)
62  IF(imssm.EQ.0) RETURN
63 
64 C...Save copy of MWID(KC) and MDCY(KC,1) values before
65 C...they are set to zero for the LSP.
66  IF(init.EQ.0) THEN
67  init=1
68  DO 100 i=1,36
69  kf=kfsusy(i)
70  kc=pycomp(kf)
71  mwidsu(i)=mwid(kc)
72  mdcysu(i)=mdcy(kc,1)
73  100 CONTINUE
74  ENDIF
75 
76 C...Restore MWID(KC) and MDCY(KC,1) values previously zeroed for LSP.
77  DO 110 i=1,36
78  kf=kfsusy(i)
79  kc=pycomp(kf)
80  IF(mdcy(kc,1).EQ.0.AND.mdcysu(i).NE.0) THEN
81  mwid(kc)=mwidsu(i)
82  mdcy(kc,1)=mdcysu(i)
83  ENDIF
84  110 CONTINUE
85 
86 C...First part of routine: set masses and couplings.
87 
88 C...Reset mixing values in sfermion sector to pure left/right.
89  DO 120 i=1,16
90  sfmix(i,1)=1d0
91  sfmix(i,4)=1d0
92  sfmix(i,2)=0d0
93  sfmix(i,3)=0d0
94  120 CONTINUE
95 
96 C...Add NMSSM states if NMSSM switched on, and change old names.
97  IF (imss(13).NE.0) THEN
98 C... Switch on NMSSM
99  WRITE(mstu(11),*) '(PYMSIN:) switching on NMSSM'
100 
101  kfn=25
102  kcn=kfn
103  chaf(kcn,1)='H_10'
104  chaf(kcn,2)=' '
105 
106  kfn=35
107  kcn=kfn
108  chaf(kcn,1)='H_20'
109  chaf(kcn,2)=' '
110 
111  kfn=45
112  kcn=kfn
113  chaf(kcn,1)='H_30'
114  chaf(kcn,2)=' '
115 
116  kfn=36
117  kcn=kfn
118  chaf(kcn,1)='A_10'
119  chaf(kcn,2)=' '
120 
121  kfn=46
122  kcn=kfn
123  chaf(kcn,1)='A_20'
124  chaf(kcn,2)=' '
125 
126  kfn=1000045
127  kcn=pycomp(kfn)
128  IF (kcn.EQ.0) THEN
129  DO 123 kct=100,mstu(6)
130  IF(kchg(kct,4).GT.100) kcn=kct
131  123 CONTINUE
132  kcn=kcn+1
133  kchg(kcn,4)=kfn
134  mstu(20)=0
135  ENDIF
136 C... Set stable for now
137  pmas(kcn,2)=1d-6
138  mwid(kcn)=0
139  mdcy(kcn,1)=0
140  mdcy(kcn,2)=0
141  mdcy(kcn,3)=0
142  chaf(kcn,1)='~chi_50'
143  chaf(kcn,2)=' '
144  ENDIF
145 
146 C...Read spectrum from SLHA file.
147  IF (imssm.EQ.11.AND.imss(21).NE.0) THEN
148 C...First check for new states
149  CALL pyslha(0,0,ifail)
150 C...Then read spectrum
151  CALL pyslha(1,0,ifail)
152  ELSEIF (imss(21).NE.0) THEN
153 C...Check for new states but don't read spectrum
154  CALL pyslha(0,0,ifail)
155  ENDIF
156 
157 C...Common couplings.
158  tanb=rmss(5)
159  beta=atan(tanb)
160  cosb=cos(beta)
161  sinb=tanb*cosb
162  cos2b=cos(2d0*beta)
163  alfa=rmss(18)
164  xmw2=pmas(24,1)**2
165  xmz2=pmas(23,1)**2
166  xw=paru(102)
167 
168 C...Define sparticle masses for a general MSSM simulation.
169  IF(imssm.EQ.1) THEN
170  IF(imss(9).EQ.0) rmss(22)=rmss(9)
171  DO 130 i=1,5,2
172  kc=pycomp(ksusy1+i)
173  pmas(kc,1)=sqrt(rmss(8)**2-(2d0*xmw2+xmz2)*cos2b/6d0)
174  kc=pycomp(ksusy2+i)
175  pmas(kc,1)=sqrt(rmss(9)**2+(xmw2-xmz2)*cos2b/3d0)
176  kc=pycomp(ksusy1+i+1)
177  pmas(kc,1)=sqrt(rmss(8)**2+(4d0*xmw2-xmz2)*cos2b/6d0)
178  kc=pycomp(ksusy2+i+1)
179  pmas(kc,1)=sqrt(rmss(22)**2-(xmw2-xmz2)*cos2b*2d0/3d0)
180  130 CONTINUE
181  xarg=rmss(6)**2-pmas(24,1)**2*abs(cos(2d0*beta))
182  IF(xarg.LT.0d0) THEN
183  WRITE(mstu(11),*) ' SNEUTRINO MASS IS NEGATIVE'//
184  & ' FROM THE SUM RULE. '
185  WRITE(mstu(11),*) ' TRY A SMALLER VALUE OF TAN(BETA). '
186  RETURN
187  ELSE
188  xarg=sqrt(xarg)
189  ENDIF
190  DO 140 i=11,15,2
191  pmas(pycomp(ksusy1+i),1)=rmss(6)
192  pmas(pycomp(ksusy2+i),1)=rmss(7)
193  pmas(pycomp(ksusy1+i+1),1)=xarg
194  pmas(pycomp(ksusy2+i+1),1)=9999d0
195  140 CONTINUE
196  IF(imss(8).EQ.1) THEN
197  rmss(13)=rmss(6)
198  rmss(14)=rmss(7)
199  ENDIF
200 
201 C...Alternatively derive masses from SUGRA relations.
202  ELSEIF(imssm.EQ.2) THEN
203  rmss(36)=rmss(16)
204  CALL pyapps
205 C...Or use ISASUSY
206  ELSEIF(imssm.EQ.12.OR.imssm.EQ.13) THEN
207  rmss(36)=rmss(16)
208  CALL pysugi
209  alfa=rmss(18)
210  goto 170
211  ELSE
212  goto 170
213  ENDIF
214 
215 C...Add in extra D-term contributions.
216  IF(imss(7).EQ.1) THEN
217  r=0.43d0
218  dx=rmss(23)
219  dy=rmss(24)
220  ds=rmss(25)
221  WRITE(mstu(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
222  WRITE(mstu(11),*) 'C NEW DTERMS ADDED TO SCALAR MASSES '
223  WRITE(mstu(11),*) 'C IN A U(B-L) THEORY '
224  WRITE(mstu(11),*) 'C DX = ',dx
225  WRITE(mstu(11),*) 'C DY = ',dy
226  WRITE(mstu(11),*) 'C DS = ',ds
227  WRITE(mstu(11),*) 'C '
228  dy=r*dy-4d0/33d0*(1d0-r)*dx+(1d0-r)/33d0*ds
229  WRITE(mstu(11),*) 'C DY AT THE WEAK SCALE = ',dy
230  WRITE(mstu(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
231  dq2=dy/6d0-dx/3d0-ds/3d0
232  du2=-2d0*dy/3d0-dx/3d0-ds/3d0
233  dd2=dy/3d0+dx-2d0*ds/3d0
234  dl2=-dy/2d0+dx-2d0*ds/3d0
235  de2=dy-dx/3d0-ds/3d0
236  dhu2=dy/2d0+2d0*dx/3d0+2d0*ds/3d0
237  dhd2=-dy/2d0-2d0*dx/3d0+ds
238  dmu2=(-dy/2d0-2d0/3d0*dx+(cosb**2-2d0*sinb**2/3d0)*ds)
239  & /abs(cos2b)
240  dma2 = 2d0*dmu2+dhu2+dhd2
241  DO 150 i=1,5,2
242  kc=pycomp(ksusy1+i)
243  pmas(kc,1)=sqrt(pmas(kc,1)**2+dq2)
244  kc=pycomp(ksusy2+i)
245  pmas(kc,1)=sqrt(pmas(kc,1)**2+dd2)
246  kc=pycomp(ksusy1+i+1)
247  pmas(kc,1)=sqrt(pmas(kc,1)**2+dq2)
248  kc=pycomp(ksusy2+i+1)
249  pmas(kc,1)=sqrt(pmas(kc,1)**2+du2)
250  150 CONTINUE
251  DO 160 i=11,15,2
252  kc=pycomp(ksusy1+i)
253  pmas(kc,1)=sqrt(pmas(kc,1)**2+dl2)
254  kc=pycomp(ksusy2+i)
255  pmas(kc,1)=sqrt(pmas(kc,1)**2+de2)
256  kc=pycomp(ksusy1+i+1)
257  pmas(kc,1)=sqrt(pmas(kc,1)**2+dl2)
258  160 CONTINUE
259  IF(rmss(4)**2+dmu2.LT.0d0) THEN
260  WRITE(mstu(11),*) ' MU2 DRIVEN NEGATIVE '
261  CALL pystop(104)
262  ENDIF
263  sgnmu=sign(1d0,rmss(4))
264  rmss(4)=sgnmu*sqrt(rmss(4)**2+dmu2)
265  arg=rmss(10)**2*sign(1d0,rmss(10))+dq2
266  rmss(10)=sign(sqrt(abs(arg)),arg)
267  arg=rmss(11)**2*sign(1d0,rmss(11))+dd2
268  rmss(11)=sign(sqrt(abs(arg)),arg)
269  arg=rmss(12)**2*sign(1d0,rmss(12))+du2
270  rmss(12)=sign(sqrt(abs(arg)),arg)
271  arg=rmss(13)**2*sign(1d0,rmss(13))+dl2
272  rmss(13)=sign(sqrt(abs(arg)),arg)
273  arg=rmss(14)**2*sign(1d0,rmss(14))+de2
274  rmss(14)=sign(sqrt(abs(arg)),arg)
275  IF( rmss(19)**2 + dma2 .LE. 50d0 ) THEN
276  WRITE(mstu(11),*) ' MA DRIVEN TOO LOW '
277  CALL pystop(104)
278  ENDIF
279  rmss(19)=sqrt(rmss(19)**2+dma2)
280  rmss(6)=sqrt(rmss(6)**2+dl2)
281  rmss(7)=sqrt(rmss(7)**2+de2)
282  WRITE(mstu(11),*) ' MTL = ',rmss(10)
283  WRITE(mstu(11),*) ' MBR = ',rmss(11)
284  WRITE(mstu(11),*) ' MTR = ',rmss(12)
285  WRITE(mstu(11),*) ' SEL = ',rmss(6),rmss(13)
286  WRITE(mstu(11),*) ' SER = ',rmss(7),rmss(14)
287  ENDIF
288 
289 C...Fix the third generation sfermions.
290  CALL pythrg
291 
292 C...Fix the neutralino--chargino--gluino sector.
293  CALL pyinom
294 
295 C...Fix the Higgs sector.
296  CALL pyhggm(alfa)
297 
298 C...Choose the Gunion-Haber convention.
299  alfa=-alfa
300  rmss(18)=alfa
301 
302 C...Print information on mass parameters.
303  IF(imssm.EQ.2.AND.mstp(122).GT.0) THEN
304  WRITE(mstu(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
305  WRITE(mstu(11),*) ' USING APPROXIMATE SUGRA RELATIONS '
306  WRITE(mstu(11),*) ' M0 = ',rmss(8)
307  WRITE(mstu(11),*) ' M1/2=',rmss(1)
308  WRITE(mstu(11),*) ' TANB=',rmss(5)
309  WRITE(mstu(11),*) ' MU = ',rmss(4)
310  WRITE(mstu(11),*) ' AT = ',rmss(16)
311  WRITE(mstu(11),*) ' MA = ',rmss(19)
312  WRITE(mstu(11),*) ' MTOP=',pmas(6,1)
313  WRITE(mstu(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
314  ENDIF
315  IF(imss(20).EQ.1) THEN
316  WRITE(mstu(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
317  WRITE(mstu(11),*) ' DEBUG MODE '
318  WRITE(mstu(11),*) ' UMIX = ',umix(1,1),umix(1,2),
319  & umix(2,1),umix(2,2)
320  WRITE(mstu(11),*) ' UMIXI = ',umixi(1,1),umixi(1,2),
321  & umixi(2,1),umixi(2,2)
322  WRITE(mstu(11),*) ' VMIX = ',vmix(1,1),vmix(1,2),
323  & vmix(2,1),vmix(2,2)
324  WRITE(mstu(11),*) ' VMIXI = ',vmixi(1,1),vmixi(1,2),
325  & vmixi(2,1),vmixi(2,2)
326  WRITE(mstu(11),*) ' ZMIX = ',(zmix(1,i),i=1,4)
327  WRITE(mstu(11),*) ' ZMIXI = ',(zmixi(1,i),i=1,4)
328  WRITE(mstu(11),*) ' ZMIX = ',(zmix(2,i),i=1,4)
329  WRITE(mstu(11),*) ' ZMIXI = ',(zmixi(2,i),i=1,4)
330  WRITE(mstu(11),*) ' ZMIX = ',(zmix(3,i),i=1,4)
331  WRITE(mstu(11),*) ' ZMIXI = ',(zmixi(3,i),i=1,4)
332  WRITE(mstu(11),*) ' ZMIX = ',(zmix(4,i),i=1,4)
333  WRITE(mstu(11),*) ' ZMIXI = ',(zmixi(4,i),i=1,4)
334  WRITE(mstu(11),*) ' ALFA = ',alfa
335  WRITE(mstu(11),*) ' BETA = ',beta
336  WRITE(mstu(11),*) ' STOP = ',(sfmix(6,i),i=1,4)
337  WRITE(mstu(11),*) ' SBOT = ',(sfmix(5,i),i=1,4)
338  WRITE(mstu(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
339  ENDIF
340 
341 C...Set up the Higgs couplings - needed here since initialization
342 C...in PYINRE did not yet occur when PYWIDT is called below.
343  170 al=alfa
344  be=beta
345  sina=sin(al)
346  cosa=cos(al)
347  cosb=cos(be)
348  sinb=tanb*cosb
349  sbma=sin(be-al)
350  sapb=sin(al+be)
351  capb=cos(al+be)
352  cbma=cos(be-al)
353  c2a=cos(2d0*al)
354  c2b=cosb**2-sinb**2
355 C...tanb (used for H+)
356  paru(141)=tanb
357 
358 C...Firstly: h
359 C...Coupling to d-type quarks
360  paru(161)=sina/cosb
361 C...Coupling to u-type quarks
362  paru(162)=-cosa/sinb
363 C...Coupling to leptons
364  paru(163)=paru(161)
365 C...Coupling to Z
366  paru(164)=sbma
367 C...Coupling to W
368  paru(165)=paru(164)
369 
370 C...Secondly: H
371 C...Coupling to d-type quarks
372  paru(171)=-cosa/cosb
373 C...Coupling to u-type quarks
374  paru(172)=-sina/sinb
375 C...Coupling to leptons
376  paru(173)=paru(171)
377 C...Coupling to Z
378  paru(174)=cbma
379 C...Coupling to W
380  paru(175)=paru(174)
381 C...Coupling to h
382  IF(imss(4).GE.2) THEN
383  paru(176)=cos(2d0*al)*cos(be+al)-2d0*sin(2d0*al)*sin(be+al)
384  ELSE
385  hhh(3)=hhh(3)+hhh(4)+hhh(5)
386  paru(176)=-3d0/hhh(1)*(hhh(1)*sina**2*cosb*cosa+
387  1 hhh(2)*cosa**2*sinb*sina+hhh(3)*(sina**3*sinb+cosa**3*cosb-
388  2 2d0/3d0*cbma)-hhh(6)*sina*(cosb*c2a+cosa*capb)+
389  3 hhh(7)*cosa*(sinb*c2a+sina*capb))
390  ENDIF
391 C...Coupling to H+
392 C...Define later
393  IF(imss(4).GE.2) THEN
394  paru(168)=-sbma-cos(2d0*be)*sapb/2d0/(1d0-xw)
395  ELSE
396  paru(168)=1d0/hhh(1)*(hhh(1)*sinb**2*cosb*sina-
397  1 hhh(2)*cosb**2*sinb*cosa-hhh(3)*(sinb**3*cosa-cosb**3*sina)+
398  2 2d0*hhh(5)*sbma-hhh(6)*sinb*(cosb*sapb+sina*c2b)-
399  3 hhh(7)*cosb*(cosa*c2b-sinb*sapb)-(hhh(5)-hhh(4))*sbma)
400  ENDIF
401 C...Coupling to A
402  IF(imss(4).GE.2) THEN
403  paru(177)=cos(2d0*be)*cos(be+al)
404  ELSE
405  paru(177)=-1d0/hhh(1)*(hhh(1)*sinb**2*cosb*cosa+
406  1 hhh(2)*cosb**2*sinb*sina+hhh(3)*(sinb**3*sina+cosb**3*cosa)-
407  2 2d0*hhh(5)*cbma-hhh(6)*sinb*(cosb*capb+cosa*c2b)+
408  3 hhh(7)*cosb*(sinb*capb+sina*c2b))
409  ENDIF
410 C...Coupling to H+
411  IF(imss(4).GE.2) THEN
412  paru(178)=paru(177)
413  ELSE
414  paru(178)=paru(177)-(hhh(5)-hhh(4))/hhh(1)*cbma
415  ENDIF
416 C...Thirdly, A
417 C...Coupling to d-type quarks
418  paru(181)=tanb
419 C...Coupling to u-type quarks
420  paru(182)=1d0/paru(181)
421 C...Coupling to leptons
422  paru(183)=paru(181)
423  paru(184)=0d0
424  paru(185)=0d0
425 C...Coupling to Z h
426  paru(186)=cos(be-al)
427 C...Coupling to Z H
428  paru(187)=sin(be-al)
429  paru(188)=0d0
430  paru(189)=0d0
431  paru(190)=0d0
432 
433 C...Finally: H+
434 C...Coupling to W h
435  paru(195)=cos(be-al)
436 
437 C...Tell that all Higgs couplings have been set.
438  mstp(4)=1
439 
440 C...Set R-Violating couplings.
441 C...Set lambda couplings to common value or "natural values".
442  IF ((imss(51).NE.3).AND.(imss(51).NE.0)) THEN
443  vir3=1d0/(126d0)**3
444  DO 200 irk=1,3
445  DO 190 iri=1,3
446  DO 180 irj=1,3
447  IF (iri.NE.irj) THEN
448  IF (iri.LT.irj) THEN
449  rvlam(iri,irj,irk)=rmss(51)
450  IF (imss(51).EQ.2) rvlam(iri,irj,irk)=rmss(51)*
451  & sqrt(pmas(9+2*iri,1)*pmas(9+2*irj,1)*
452  & pmas(9+2*irk,1)*vir3)
453  ELSE
454  rvlam(iri,irj,irk)=-rvlam(irj,iri,irk)
455  ENDIF
456  ELSE
457  rvlam(iri,irj,irk)=0d0
458  ENDIF
459  180 CONTINUE
460  190 CONTINUE
461  200 CONTINUE
462  ENDIF
463 C...Set lambda' couplings to common value or "natural values".
464  IF ((imss(52).NE.3).AND.(imss(52).NE.0)) THEN
465  vir3=1d0/(126d0)**3
466  DO 230 iri=1,3
467  DO 220 irj=1,3
468  DO 210 irk=1,3
469  rvlamp(iri,irj,irk)=rmss(52)
470  IF (imss(52).EQ.2) rvlamp(iri,irj,irk)=rmss(52)*
471  & sqrt(pmas(9+2*iri,1)*0.5d0*(pmas(2*irj,1)+
472  & pmas(2*irj-1,1))*pmas(2*irk-1,1)*vir3)
473  210 CONTINUE
474  220 CONTINUE
475  230 CONTINUE
476  ENDIF
477 C...Set lambda'' couplings to common value or "natural values".
478  IF ((imss(53).NE.3).AND.(imss(53).NE.0)) THEN
479  vir3=1d0/(126d0)**3
480  DO 260 iri=1,3
481  DO 250 irj=1,3
482  DO 240 irk=1,3
483  IF (irj.NE.irk) THEN
484  IF (irj.LT.irk) THEN
485  rvlamb(iri,irj,irk)=rmss(53)
486  IF (imss(53).EQ.2) rvlamb(iri,irj,irk)=
487  & rmss(53)*sqrt(pmas(2*iri,1)*pmas(2*irj-1,1)*
488  & pmas(2*irk-1,1)*vir3)
489  ELSE
490  rvlamb(iri,irj,irk)=-rvlamb(iri,irk,irj)
491  ENDIF
492  ELSE
493  rvlamb(iri,irj,irk) = 0d0
494  ENDIF
495  240 CONTINUE
496  250 CONTINUE
497  260 CONTINUE
498  ENDIF
499 
500 C...Antisymmetrize couplings set by user
501  IF (imss(51).EQ.3.OR.imss(53).EQ.3) THEN
502  DO 290 iri=1,3
503  DO 280 irj=1,3
504  DO 270 irk=1,3
505  IF (rvlam(iri,irj,irk).NE.-rvlam(irj,iri,irk)) THEN
506  rvlam(irj,iri,irk)=-rvlam(iri,irj,irk)
507  IF (iri.EQ.irj) rvlam(iri,irj,irk)=0d0
508  ENDIF
509  IF (rvlamb(iri,irj,irk).NE.-rvlamb(iri,irk,irj)) THEN
510  rvlamb(iri,irk,irj)=-rvlamb(iri,irj,irk)
511  IF (irj.EQ.irk) rvlamb(iri,irj,irk)=0d0
512  ENDIF
513  270 CONTINUE
514  280 CONTINUE
515  290 CONTINUE
516  ENDIF
517 
518 C...Write spectrum to SLHA file
519  IF (imss(23).NE.0) THEN
520  ifail=0
521  CALL pyslha(3,0,ifail)
522  ENDIF
523 
524 C...Second part of routine: set decay modes and branching ratios.
525 
526 C...Allow chi10 -> gravitino + gamma or not.
527  kc=pycomp(ksusy1+39)
528  IF( imss(11) .NE. 0 ) THEN
529  pmas(kc,1)=rmss(21)/1d9
530  pmas(kc,2)=0d0
531  irprty=0
532  WRITE(mstu(11),*) ' ALLOWING DECAYS TO GRAVITINOS '
533  ELSE IF (imss(51).GE.1.OR.imss(52).GE.1.OR.imss(53).GE.1) THEN
534  irprty=0
535  IF (imss(51).GE.1) WRITE(mstu(11),*)
536  & ' ALLOWING SUSY LLE DECAYS'
537  IF (imss(52).GE.1) WRITE(mstu(11),*)
538  & ' ALLOWING SUSY LQD DECAYS'
539  IF (imss(53).GE.1) WRITE(mstu(11),*)
540  & ' ALLOWING SUSY UDD DECAYS'
541  IF (imss(53).GE.1.AND.imss(52).GE.1) WRITE(mstu(11),*)
542  & ' --- Warning: R-Violating couplings possibly',
543  & ' incompatible with proton decay'
544  ELSE
545  pmas(kc,1)=9999d0
546  irprty=1
547  ENDIF
548 
549 C...Loop over sparticle and Higgs species.
550  pmchi1=pmas(pycomp(ksusy1+22),1)
551 C...Find the LSP or NLSP for a gravitino LSP
552  ilsp=0
553  pmlsp=1d20
554  DO 300 i=1,36
555  kf=kfsusy(i)
556  IF(kf.EQ.1000039) goto 300
557  kc=pycomp(kf)
558  IF(pmas(kc,1).LT.pmlsp) THEN
559  ilsp=i
560  pmlsp=pmas(kc,1)
561  ENDIF
562  300 CONTINUE
563  DO 370 i=1,50
564  IF (i.GT.39.AND.imss(13).NE.1) goto 370
565  kf=kfsusy(i)
566  IF (kf.EQ.0) goto 370
567  kc=pycomp(kf)
568  lknt=0
569 
570 C...Check if there are any decays listed for this sparticle
571 C...in a file
572  IF (imss(22).NE.0) THEN
573  ifail=0
574 C...First look for MASS entry if not already done
575  IF (imss(1).NE.11.AND.imss(21).NE.0) CALL pyslha(5,kf,ifail)
576 C...Then look for decay info
577  ifail=0
578  CALL pyslha(2,kf,ifail)
579  IF (ifail.EQ.0.OR.kf.EQ.6.OR.kf.EQ.24) goto 370
580  ELSEIF (i.GE.37) THEN
581  goto 370
582  ENDIF
583 
584 C...Sfermion decays.
585  IF(i.LE.24) THEN
586 C...First check to see if sneutrino is lighter than chi10.
587  IF((i.EQ.15.OR.i.EQ.19.OR.i.EQ.23).AND.
588  & pmas(kc,1).LT.pmchi1) THEN
589  ELSE
590  CALL pysfdc(kf,xlam,idlam,lknt)
591  ENDIF
592 
593 C...Gluino decays.
594  ELSEIF(i.EQ.25) THEN
595  CALL pyglui(kf,xlam,idlam,lknt)
596  IF(i.EQ.ilsp.AND.irprty.EQ.1) lknt=0
597 
598 C...Neutralino decays.
599  ELSEIF(i.GE.26.AND.i.LE.29) THEN
600  CALL pynjdc(kf,xlam,idlam,lknt)
601 C...chi10 stable or chi10 -> gravitino + gamma.
602  IF(i.EQ.26.AND.irprty.EQ.1) THEN
603  pmas(kc,2)=1d-6
604  mdcy(kc,1)=0
605  mwid(kc)=0
606  ENDIF
607 
608 C...Chargino decays.
609  ELSEIF(i.GE.30.AND.i.LE.31) THEN
610  CALL pycjdc(kf,xlam,idlam,lknt)
611 
612 C...Gravitino is stable.
613  ELSEIF(i.EQ.32) THEN
614  mdcy(kc,1)=0
615  mwid(kc)=0
616 
617 C...Higgs decays.
618  ELSEIF(i.GE.33.AND.i.LE.36) THEN
619 C...Calculate decays to non-SUSY particles.
620  CALL pywidt(kf,pmas(kc,1)**2,wdtp,wdte)
621  lknt=0
622  DO 310 i1=0,100
623  xlam(i1)=0d0
624  310 CONTINUE
625  DO 330 i1=1,mdcy(kc,3)
626  k1=mdcy(kc,2)+i1-1
627  IF(iabs(kfdp(k1,1)).GT.ksusy1.OR.
628  & iabs(kfdp(k1,2)).GT.ksusy1) goto 330
629  xlam(i1)=wdtp(i1)
630  xlam(0)=xlam(0)+xlam(i1)
631  DO 320 j1=1,3
632  idlam(i1,j1)=kfdp(k1,j1)
633  320 CONTINUE
634  lknt=lknt+1
635  330 CONTINUE
636 C...Add the decays to SUSY particles.
637  CALL pyhext(kf,xlam,idlam,lknt)
638  ENDIF
639 C...Zero the branching ratios for use in loop mode
640 C...thanks to K. Matchev (FNAL)
641  DO 340 idc=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1
642  brat(idc)=0d0
643  340 CONTINUE
644 
645 C...Set stable particles.
646  IF(lknt.EQ.0) THEN
647  mdcy(kc,1)=0
648  mwid(kc)=0
649  pmas(kc,2)=1d-6
650  pmas(kc,3)=1d-5
651  pmas(kc,4)=0d0
652 
653 C...Store branching ratios in the standard tables.
654  ELSE
655  idc=mdcy(kc,2)+mdcy(kc,3)-1
656  delm=1d6
657  DO 360 il=1,lknt
658  idcsv=idc
659  350 idc=idc+1
660  brat(idc)=0d0
661  IF(idc.EQ.mdcy(kc,2)+mdcy(kc,3)) idc=mdcy(kc,2)
662  IF(idlam(il,1).EQ.kfdp(idc,1).AND.idlam(il,2).EQ.
663  & kfdp(idc,2).AND.idlam(il,3).EQ.kfdp(idc,3)) THEN
664  brat(idc)=xlam(il)/xlam(0)
665  xmdif=pmas(kc,1)
666  IF(mdme(idc,1).GE.1) THEN
667  xmdif=xmdif-pmas(pycomp(kfdp(idc,1)),1)-
668  & pmas(pycomp(kfdp(idc,2)),1)
669  IF(kfdp(idc,3).NE.0) xmdif=xmdif-
670  & pmas(pycomp(kfdp(idc,3)),1)
671  ENDIF
672  IF(i.LE.32) THEN
673  IF(xmdif.GE.0d0) THEN
674  delm=min(delm,xmdif)
675  ELSE
676  WRITE(mstu(11),*) ' ERROR WITH DELM ',delm,xmdif
677  WRITE(mstu(11),*) ' KF = ',kf
678  WRITE(mstu(11),*) ' KF(decay) = ',(kfdp(idc,j),j=1,3)
679  ENDIF
680  ENDIF
681  goto 360
682  ELSEIF(idc.EQ.idcsv) THEN
683  WRITE(mstu(11),*) ' Error in PYMSIN: SUSY decay ',
684  & 'channel not recognized:'
685  WRITE(mstu(11),*) kf,' -> ',(idlam(il,j),j=1,3)
686  goto 360
687  ELSE
688  goto 350
689  ENDIF
690  360 CONTINUE
691 
692 C...Store width, cutoff and lifetime.
693  pmas(kc,2)=xlam(0)
694  IF(pmas(kc,2).LT.0.1d0*delm) THEN
695  pmas(kc,3)=pmas(kc,2)*10d0
696  ELSE
697  pmas(kc,3)=0.95d0*delm
698  ENDIF
699  IF(pmas(kc,2).NE.0d0) THEN
700  pmas(kc,4)=paru(3)/pmas(kc,2)*1d-12
701  ENDIF
702 C...Write decays to SLHA file
703  IF (imss(24).NE.0) THEN
704  ifail=0
705  CALL pyslha(4,kf,ifail)
706  ENDIF
707 
708  ENDIF
709  370 CONTINUE
710 
711  RETURN
712  END