Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyinre.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyinre.f
1 
2 C*********************************************************************
3 
4 C...PYINRE
5 C...Calculates full and effective widths of gauge bosons, stores
6 C...masses and widths, rescales coefficients to be used for
7 C...resonance production generation.
8 
9  SUBROUTINE pyinre
10 
11 C...Double precision and integer declarations.
12  IMPLICIT DOUBLE PRECISION(a-h, o-z)
13  IMPLICIT INTEGER(i-n)
14  INTEGER pyk,pychge,pycomp
15 C...Parameter statement to help give large particle numbers.
16  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
17  &kexcit=4000000,kdimen=5000000)
18 C...Commonblocks.
19  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
20  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
21  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
22  common/pydat4/chaf(500,2)
23  CHARACTER chaf*16
24  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
25  common/pypars/mstp(200),parp(200),msti(200),pari(200)
26  common/pyint1/mint(400),vint(400)
27  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
28  common/pyint4/mwid(500),wids(500,5)
29  common/pyint6/proc(0:500)
30  CHARACTER proc*28
31  common/pymssm/imss(0:99),rmss(0:99)
32  SAVE /pydat1/,/pydat2/,/pydat3/,/pydat4/,/pysubs/,/pypars/,
33  &/pyint1/,/pyint2/,/pyint4/,/pyint6/,/pymssm/
34 C...Local arrays and data.
35  dimension wdtp(0:400),wdte(0:400,0:5),wdtpm(0:400),
36  &wdtem(0:400,0:5),kcord(500),pmord(500)
37 
38 C...Born level couplings in MSSM Higgs doublet sector.
39  xw=paru(102)
40  xwv=xw
41  IF(mstp(8).GE.2) xw=1d0-(pmas(24,1)/pmas(23,1))**2
42  xw1=1d0-xw
43  IF(mstp(4).EQ.2) THEN
44  tanbe=paru(141)
45  ratbe=((1d0-tanbe**2)/(1d0+tanbe**2))**2
46  sqmz=pmas(23,1)**2
47  sqmw=pmas(24,1)**2
48  sqmh=pmas(25,1)**2
49  sqma=sqmh*(sqmz-sqmh)/(sqmz*ratbe-sqmh)
50  sqmhp=0.5d0*(sqma+sqmz+sqrt((sqma+sqmz)**2-4d0*sqma*sqmz*ratbe))
51  sqmhc=sqma+sqmw
52  IF(sqmh.GE.sqmz.OR.min(sqma,sqmhp,sqmhc).LE.0d0) THEN
53  WRITE(mstu(11),5000)
54  CALL pystop(101)
55  ENDIF
56  pmas(35,1)=sqrt(sqmhp)
57  pmas(36,1)=sqrt(sqma)
58  pmas(37,1)=sqrt(sqmhc)
59  alsu=0.5d0*atan(2d0*tanbe*(sqma+sqmz)/((1d0-tanbe**2)*
60  & (sqma-sqmz)))
61  besu=atan(tanbe)
62  paru(142)=1d0
63  paru(143)=1d0
64  paru(161)=-sin(alsu)/cos(besu)
65  paru(162)=cos(alsu)/sin(besu)
66  paru(163)=paru(161)
67  paru(164)=sin(besu-alsu)
68  paru(165)=paru(164)
69  paru(168)=sin(besu-alsu)+0.5d0*cos(2d0*besu)*sin(besu+alsu)/xw
70  paru(171)=cos(alsu)/cos(besu)
71  paru(172)=sin(alsu)/sin(besu)
72  paru(173)=paru(171)
73  paru(174)=cos(besu-alsu)
74  paru(175)=paru(174)
75  paru(176)=cos(2d0*alsu)*cos(besu+alsu)-2d0*sin(2d0*alsu)*
76  & sin(besu+alsu)
77  paru(177)=cos(2d0*besu)*cos(besu+alsu)
78  paru(178)=cos(besu-alsu)-0.5d0*cos(2d0*besu)*cos(besu+alsu)/xw
79  paru(181)=tanbe
80  paru(182)=1d0/tanbe
81  paru(183)=paru(181)
82  paru(184)=0d0
83  paru(185)=paru(184)
84  paru(186)=cos(besu-alsu)
85  paru(187)=sin(besu-alsu)
86  paru(188)=paru(186)
87  paru(189)=paru(187)
88  paru(190)=0d0
89  paru(195)=cos(besu-alsu)
90  ENDIF
91 
92 C...Reset effective widths of gauge bosons.
93  DO 110 i=1,500
94  DO 100 j=1,5
95  wids(i,j)=1d0
96  100 CONTINUE
97  110 CONTINUE
98 
99 C...Order resonances by increasing mass (except Z0 and W+/-).
100  nres=0
101  DO 140 kc=1,500
102  kf=kchg(kc,4)
103  IF(kf.EQ.0) goto 140
104  IF(mwid(kc).EQ.0) goto 140
105  IF(kc.EQ.7.OR.kc.EQ.8.OR.kc.EQ.17.OR.kc.EQ.18) THEN
106  IF(mstp(1).LE.3) goto 140
107  ENDIF
108  IF(kf/ksusy1.EQ.1.OR.kf/ksusy1.EQ.2) THEN
109  IF(imss(1).LE.0) goto 140
110  ENDIF
111  nres=nres+1
112  pmres=pmas(kc,1)
113  IF(kc.EQ.23.OR.kc.EQ.24) pmres=0d0
114  DO 120 i1=nres-1,1,-1
115  IF(pmres.GE.pmord(i1)) goto 130
116  kcord(i1+1)=kcord(i1)
117  pmord(i1+1)=pmord(i1)
118  120 CONTINUE
119  130 kcord(i1+1)=kc
120  pmord(i1+1)=pmres
121  140 CONTINUE
122 
123 C...Loop over possible resonances.
124  DO 180 i=1,nres
125  kc=kcord(i)
126  kf=kchg(kc,4)
127 
128 C...Check that no fourth generation channels on by mistake.
129  IF(mstp(1).LE.3) THEN
130  DO 150 j=1,mdcy(kc,3)
131  idc=j+mdcy(kc,2)-1
132  kfa1=iabs(kfdp(idc,1))
133  kfa2=iabs(kfdp(idc,2))
134  IF(kfa1.EQ.7.OR.kfa1.EQ.8.OR.kfa1.EQ.17.OR.kfa1.EQ.18.OR.
135  & kfa2.EQ.7.OR.kfa2.EQ.8.OR.kfa2.EQ.17.OR.kfa2.EQ.18)
136  & mdme(idc,1)=-1
137  150 CONTINUE
138  ENDIF
139 
140 C...Check that no supersymmetric channels on by mistake.
141  IF(imss(1).LE.0) THEN
142  DO 160 j=1,mdcy(kc,3)
143  idc=j+mdcy(kc,2)-1
144  kfa1s=iabs(kfdp(idc,1))/ksusy1
145  kfa2s=iabs(kfdp(idc,2))/ksusy1
146  IF(kfa1s.EQ.1.OR.kfa1s.EQ.2.OR.kfa2s.EQ.1.OR.kfa2s.EQ.2)
147  & mdme(idc,1)=-1
148  160 CONTINUE
149  ENDIF
150 
151 C...Find mass and evaluate width.
152  pmr=pmas(kc,1)
153  IF(kf.EQ.25.OR.kf.EQ.35.OR.kf.EQ.36) mint(62)=1
154  IF(mwid(kc).EQ.3) mint(63)=1
155  CALL pywidt(kf,pmr**2,wdtp,wdte)
156  mint(51)=0
157 
158 C...Evaluate suppression factors due to non-simulated channels.
159  IF(kchg(kc,3).EQ.0) THEN
160  wdtp0i=0d0
161  IF(wdtp(0).GT.0d0) wdtp0i=1d0/wdtp(0)
162  wids(kc,1)=((wdte(0,1)+wdte(0,2))**2+
163  & 2d0*(wdte(0,1)+wdte(0,2))*(wdte(0,4)+wdte(0,5))+
164  & 2d0*wdte(0,4)*wdte(0,5))*wdtp0i**2
165  wids(kc,2)=(wdte(0,1)+wdte(0,2)+wdte(0,4))*wdtp0i
166  wids(kc,3)=0d0
167  wids(kc,4)=0d0
168  wids(kc,5)=0d0
169  ELSE
170  IF(mwid(kc).EQ.3) mint(63)=1
171  CALL pywidt(-kf,pmr**2,wdtpm,wdtem)
172  mint(51)=0
173  wdtp0i=0d0
174  IF(wdtp(0).GT.0d0) wdtp0i=1d0/wdtp(0)
175  wids(kc,1)=((wdte(0,1)+wdte(0,2))*(wdtem(0,1)+wdtem(0,3))+
176  & (wdte(0,1)+wdte(0,2))*(wdtem(0,4)+wdtem(0,5))+
177  & (wdte(0,4)+wdte(0,5))*(wdtem(0,1)+wdtem(0,3))+
178  & wdte(0,4)*wdtem(0,5)+wdte(0,5)*wdtem(0,4))*wdtp0i**2
179  wids(kc,2)=(wdte(0,1)+wdte(0,2)+wdte(0,4))*wdtp0i
180  wids(kc,3)=(wdtem(0,1)+wdtem(0,3)+wdtem(0,4))*wdtp0i
181  wids(kc,4)=((wdte(0,1)+wdte(0,2))**2+
182  & 2d0*(wdte(0,1)+wdte(0,2))*(wdte(0,4)+wdte(0,5))+
183  & 2d0*wdte(0,4)*wdte(0,5))*wdtp0i**2
184  wids(kc,5)=((wdtem(0,1)+wdtem(0,3))**2+
185  & 2d0*(wdtem(0,1)+wdtem(0,3))*(wdtem(0,4)+wdtem(0,5))+
186  & 2d0*wdtem(0,4)*wdtem(0,5))*wdtp0i**2
187  ENDIF
188 
189 C...Set resonance widths and branching ratios;
190 C...also on/off switch for decays.
191  IF(mwid(kc).EQ.1.OR.mwid(kc).EQ.3) THEN
192  pmas(kc,2)=wdtp(0)
193  pmas(kc,3)=min(0.9d0*pmas(kc,1),10d0*pmas(kc,2))
194  IF(mstp(41).EQ.0.OR.mstp(41).EQ.1) mdcy(kc,1)=mstp(41)
195  DO 170 j=1,mdcy(kc,3)
196  idc=j+mdcy(kc,2)-1
197  brat(idc)=0d0
198  IF(wdtp(0).GT.0d0) brat(idc)=wdtp(j)/wdtp(0)
199  170 CONTINUE
200  ENDIF
201  180 CONTINUE
202 
203 C...Flavours of leptoquark: redefine charge and name.
204  kflqq=kfdp(mdcy(42,2),1)
205  kflql=kfdp(mdcy(42,2),2)
206  kchg(42,1)=kchg(pycomp(kflqq),1)*isign(1,kflqq)+
207  &kchg(pycomp(kflql),1)*isign(1,kflql)
208  ll=1
209  IF(iabs(kflql).EQ.13) ll=2
210  IF(iabs(kflql).EQ.15) ll=3
211  chaf(42,1)='LQ_'//chaf(iabs(kflqq),1)(1:1)//
212  &chaf(iabs(kflql),1)(1:ll)//' '
213  chaf(42,2)=chaf(42,2)(1:4+ll)//'bar '
214 
215 C...Special cases in treatment of gamma*/Z0: redefine process name.
216  IF(mstp(43).EQ.1) THEN
217  proc(1)='f + fbar -> gamma*'
218  proc(15)='f + fbar -> g + gamma*'
219  proc(19)='f + fbar -> gamma + gamma*'
220  proc(30)='f + g -> f + gamma*'
221  proc(35)='f + gamma -> f + gamma*'
222  ELSEIF(mstp(43).EQ.2) THEN
223  proc(1)='f + fbar -> Z0'
224  proc(15)='f + fbar -> g + Z0'
225  proc(19)='f + fbar -> gamma + Z0'
226  proc(30)='f + g -> f + Z0'
227  proc(35)='f + gamma -> f + Z0'
228  ELSEIF(mstp(43).EQ.3) THEN
229  proc(1)='f + fbar -> gamma*/Z0'
230  proc(15)='f + fbar -> g + gamma*/Z0'
231  proc(19)='f+ fbar -> gamma + gamma*/Z0'
232  proc(30)='f + g -> f + gamma*/Z0'
233  proc(35)='f + gamma -> f + gamma*/Z0'
234  ENDIF
235 
236 C...Special cases in treatment of gamma*/Z0/Z'0: redefine process name.
237  IF(mstp(44).EQ.1) THEN
238  proc(141)='f + fbar -> gamma*'
239  ELSEIF(mstp(44).EQ.2) THEN
240  proc(141)='f + fbar -> Z0'
241  ELSEIF(mstp(44).EQ.3) THEN
242  proc(141)='f + fbar -> Z''0'
243  ELSEIF(mstp(44).EQ.4) THEN
244  proc(141)='f + fbar -> gamma*/Z0'
245  ELSEIF(mstp(44).EQ.5) THEN
246  proc(141)='f + fbar -> gamma*/Z''0'
247  ELSEIF(mstp(44).EQ.6) THEN
248  proc(141)='f + fbar -> Z0/Z''0'
249  ELSEIF(mstp(44).EQ.7) THEN
250  proc(141)='f + fbar -> gamma*/Z0/Z''0'
251  ENDIF
252 
253 C...Special cases in treatment of WW -> WW: redefine process name.
254  IF(mstp(45).EQ.1) THEN
255  proc(77)='W+ + W+ -> W+ + W+'
256  ELSEIF(mstp(45).EQ.2) THEN
257  proc(77)='W+ + W- -> W+ + W-'
258  ELSEIF(mstp(45).EQ.3) THEN
259  proc(77)='W+/- + W+/- -> W+/- + W+/-'
260  ENDIF
261 
262 C...Format for error information.
263  5000 FORMAT(1x,'Error: unphysical input tan^2(beta) and m_H ',
264  &'combination'/1x,'Execution stopped!')
265 
266  RETURN
267  END