Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyupda.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyupda.f
1 
2 C*********************************************************************
3 
4 C...PYUPDA
5 C...Facilitates the updating of particle and decay data
6 C...by allowing it to be done in an external file.
7 
8  SUBROUTINE pyupda(MUPDA,LFN)
9 
10 C...Double precision and integer declarations.
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13  INTEGER pyk,pychge,pycomp
14 C...Commonblocks.
15  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
16  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
17  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
18  common/pydat4/chaf(500,2)
19  CHARACTER chaf*16
20  common/pyint4/mwid(500),wids(500,5)
21  SAVE /pydat1/,/pydat2/,/pydat3/,/pydat4/,/pyint4/
22 C...Local arrays, character variables and data.
23  CHARACTER chinl*120,chkf*9,chvar(22)*9,chlin*72,
24  &chblk(20)*72,chold*16,chtmp*16,chnew*16,chcom*24
25  DATA chvar/ 'KCHG(I,1)','KCHG(I,2)','KCHG(I,3)','KCHG(I,4)',
26  &'PMAS(I,1)','PMAS(I,2)','PMAS(I,3)','PMAS(I,4)','MDCY(I,1)',
27  &'MDCY(I,2)','MDCY(I,3)','MDME(I,1)','MDME(I,2)','BRAT(I) ',
28  &'KFDP(I,1)','KFDP(I,2)','KFDP(I,3)','KFDP(I,4)','KFDP(I,5)',
29  &'CHAF(I,1)','CHAF(I,2)','MWID(I) '/
30 
31 C...Write header if not yet done.
32  IF(mstu(12).NE.12345) CALL pylist(0)
33 
34 C...Write information on file for editing.
35  IF(mupda.EQ.1) THEN
36  DO 110 kc=1,500
37  WRITE(lfn,5000) kchg(kc,4),(chaf(kc,j1),j1=1,2),
38  & (kchg(kc,j2),j2=1,3),(pmas(kc,j3),j3=1,4),
39  & mwid(kc),mdcy(kc,1)
40  DO 100 idc=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1
41  WRITE(lfn,5100) mdme(idc,1),mdme(idc,2),brat(idc),
42  & (kfdp(idc,j),j=1,5)
43  100 CONTINUE
44  110 CONTINUE
45 
46 C...Read complete set of information from edited file or
47 C...read partial set of new or updated information from edited file.
48  ELSEIF(mupda.EQ.2.OR.mupda.EQ.3) THEN
49 
50 C...Reset counters.
51  kcc=100
52  ndc=0
53  chkf=' '
54  IF(mupda.EQ.2) THEN
55  DO 120 i=1,mstu(6)
56  kchg(i,4)=0
57  120 CONTINUE
58  ELSE
59  DO 130 kc=1,mstu(6)
60  IF(kc.GT.100.AND.kchg(kc,4).GT.100) kcc=kc
61  ndc=max(ndc,mdcy(kc,2)+mdcy(kc,3)-1)
62  130 CONTINUE
63  ENDIF
64 
65 C...Begin of loop: read new line; unknown whether particle or
66 C...decay data.
67  140 READ(lfn,5200,end=190) chinl
68 
69 C...Identify particle code and whether already defined (for MUPDA=3).
70  IF(chinl(2:10).NE.' ') THEN
71  chkf=chinl(2:10)
72  READ(chkf,5300) kf
73  IF(mupda.EQ.2) THEN
74  IF(kf.LE.100) THEN
75  kc=kf
76  ELSE
77  kcc=kcc+1
78  kc=kcc
79  ENDIF
80  ELSE
81  kcrep=0
82  IF(kf.LE.100) THEN
83  kcrep=kf
84  ELSE
85  DO 150 kcr=101,kcc
86  IF(kchg(kcr,4).EQ.kf) kcrep=kcr
87  150 CONTINUE
88  ENDIF
89 C...Remove duplicate old decay data.
90  IF(kcrep.NE.0.AND.mdcy(kcrep,3).GT.0) THEN
91  idcrep=mdcy(kcrep,2)
92  ndcrep=mdcy(kcrep,3)
93  DO 160 i=1,kcc
94  IF(mdcy(i,2).GT.idcrep) mdcy(i,2)=mdcy(i,2)-ndcrep
95  160 CONTINUE
96  DO 180 i=idcrep,ndc-ndcrep
97  mdme(i,1)=mdme(i+ndcrep,1)
98  mdme(i,2)=mdme(i+ndcrep,2)
99  brat(i)=brat(i+ndcrep)
100  DO 170 j=1,5
101  kfdp(i,j)=kfdp(i+ndcrep,j)
102  170 CONTINUE
103  180 CONTINUE
104  ndc=ndc-ndcrep
105  kc=kcrep
106  ELSEIF(kcrep.NE.0) THEN
107  kc=kcrep
108  ELSE
109  kcc=kcc+1
110  kc=kcc
111  ENDIF
112  ENDIF
113 
114 C...Study line with particle data.
115  IF(kc.GT.mstu(6)) CALL pyerrm(27,
116  & '(PYUPDA:) Particle arrays full by KF ='//chkf)
117  READ(chinl,5000) kchg(kc,4),(chaf(kc,j1),j1=1,2),
118  & (kchg(kc,j2),j2=1,3),(pmas(kc,j3),j3=1,4),
119  & mwid(kc),mdcy(kc,1)
120  mdcy(kc,2)=0
121  mdcy(kc,3)=0
122 
123 C...Study line with decay data.
124  ELSE
125  ndc=ndc+1
126  IF(ndc.GT.mstu(7)) CALL pyerrm(27,
127  & '(PYUPDA:) Decay data arrays full by KF ='//chkf)
128  IF(mdcy(kc,2).EQ.0) mdcy(kc,2)=ndc
129  mdcy(kc,3)=mdcy(kc,3)+1
130  READ(chinl,5100) mdme(ndc,1),mdme(ndc,2),brat(ndc),
131  & (kfdp(ndc,j),j=1,5)
132  ENDIF
133 
134 C...End of loop; ensure that PYCOMP tables are updated.
135  goto 140
136  190 CONTINUE
137  mstu(20)=0
138 
139 C...Perform possible tests that new information is consistent.
140  DO 220 kc=1,mstu(6)
141  kf=kchg(kc,4)
142  IF(kf.EQ.0) goto 220
143  WRITE(chkf,5300) kf
144  IF(min(pmas(kc,1),pmas(kc,2),pmas(kc,3),pmas(kc,1)-pmas(kc,3),
145  & pmas(kc,4)).LT.0d0.OR.mdcy(kc,3).LT.0) CALL pyerrm(17,
146  & '(PYUPDA:) Mass/width/life/(# channels) wrong for KF ='//chkf)
147  brsum=0d0
148  DO 210 idc=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1
149  IF(mdme(idc,2).GT.80) goto 210
150  kq=kchg(kc,1)
151  pms=pmas(kc,1)-pmas(kc,3)-parj(64)
152  merr=0
153  DO 200 j=1,5
154  kp=kfdp(idc,j)
155  IF(kp.EQ.0.OR.kp.EQ.81.OR.iabs(kp).EQ.82) THEN
156  IF(kp.EQ.81) kq=0
157  ELSEIF(pycomp(kp).EQ.0) THEN
158  merr=3
159  ELSE
160  kq=kq-pychge(kp)
161  kpc=pycomp(kp)
162  pms=pms-pmas(kpc,1)
163  IF(mstj(24).GT.0) pms=pms+0.5d0*min(pmas(kpc,2),
164  & pmas(kpc,3))
165  ENDIF
166  200 CONTINUE
167  IF(kq.NE.0) merr=max(2,merr)
168  IF(mwid(kc).EQ.0.AND.kf.NE.311.AND.pms.LT.0d0)
169  & merr=max(1,merr)
170  IF(merr.EQ.3) CALL pyerrm(17,
171  & '(PYUPDA:) Unknown particle code in decay of KF ='//chkf)
172  IF(merr.EQ.2) CALL pyerrm(17,
173  & '(PYUPDA:) Charge not conserved in decay of KF ='//chkf)
174  IF(merr.EQ.1) CALL pyerrm(7,
175  & '(PYUPDA:) Kinematically unallowed decay of KF ='//chkf)
176  brsum=brsum+brat(idc)
177  210 CONTINUE
178  WRITE(chtmp,5500) brsum
179  IF(abs(brsum).GT.0.0005d0.AND.abs(brsum-1d0).GT.0.0005d0)
180  & CALL pyerrm(7,'(PYUPDA:) Sum of branching ratios is '//
181  & chtmp(9:16)//' for KF ='//chkf)
182  220 CONTINUE
183 
184 C...Write DATA statements for inclusion in program.
185  ELSEIF(mupda.EQ.4) THEN
186 
187 C...Find out how many codes and decay channels are actually used.
188  kcc=0
189  ndc=0
190  DO 230 i=1,mstu(6)
191  IF(kchg(i,4).NE.0) THEN
192  kcc=i
193  ndc=max(ndc,mdcy(i,2)+mdcy(i,3)-1)
194  ENDIF
195  230 CONTINUE
196 
197 C...Initialize writing of DATA statements for inclusion in program.
198  DO 300 ivar=1,22
199  ndim=mstu(6)
200  IF(ivar.GE.12.AND.ivar.LE.19) ndim=mstu(7)
201  nlin=1
202  chlin=' '
203  chlin(7:35)='DATA ('//chvar(ivar)//',I= 1, )/'
204  llin=35
205  chold='START'
206 
207 C...Loop through variables for conversion to characters.
208  DO 280 idim=1,ndim
209  IF(ivar.EQ.1) WRITE(chtmp,5400) kchg(idim,1)
210  IF(ivar.EQ.2) WRITE(chtmp,5400) kchg(idim,2)
211  IF(ivar.EQ.3) WRITE(chtmp,5400) kchg(idim,3)
212  IF(ivar.EQ.4) WRITE(chtmp,5400) kchg(idim,4)
213  IF(ivar.EQ.5) WRITE(chtmp,5500) pmas(idim,1)
214  IF(ivar.EQ.6) WRITE(chtmp,5500) pmas(idim,2)
215  IF(ivar.EQ.7) WRITE(chtmp,5500) pmas(idim,3)
216  IF(ivar.EQ.8) WRITE(chtmp,5500) pmas(idim,4)
217  IF(ivar.EQ.9) WRITE(chtmp,5400) mdcy(idim,1)
218  IF(ivar.EQ.10) WRITE(chtmp,5400) mdcy(idim,2)
219  IF(ivar.EQ.11) WRITE(chtmp,5400) mdcy(idim,3)
220  IF(ivar.EQ.12) WRITE(chtmp,5400) mdme(idim,1)
221  IF(ivar.EQ.13) WRITE(chtmp,5400) mdme(idim,2)
222  IF(ivar.EQ.14) WRITE(chtmp,5600) brat(idim)
223  IF(ivar.EQ.15) WRITE(chtmp,5400) kfdp(idim,1)
224  IF(ivar.EQ.16) WRITE(chtmp,5400) kfdp(idim,2)
225  IF(ivar.EQ.17) WRITE(chtmp,5400) kfdp(idim,3)
226  IF(ivar.EQ.18) WRITE(chtmp,5400) kfdp(idim,4)
227  IF(ivar.EQ.19) WRITE(chtmp,5400) kfdp(idim,5)
228  IF(ivar.EQ.20) chtmp=chaf(idim,1)
229  IF(ivar.EQ.21) chtmp=chaf(idim,2)
230  IF(ivar.EQ.22) WRITE(chtmp,5400) mwid(idim)
231 
232 C...Replace variables beyond what is properly defined.
233  IF(ivar.LE.4) THEN
234  IF(idim.GT.kcc) chtmp=' 0'
235  ELSEIF(ivar.LE.8) THEN
236  IF(idim.GT.kcc) chtmp=' 0.0'
237  ELSEIF(ivar.LE.11) THEN
238  IF(idim.GT.kcc) chtmp=' 0'
239  ELSEIF(ivar.LE.13) THEN
240  IF(idim.GT.ndc) chtmp=' 0'
241  ELSEIF(ivar.LE.14) THEN
242  IF(idim.GT.ndc) chtmp=' 0.0'
243  ELSEIF(ivar.LE.19) THEN
244  IF(idim.GT.ndc) chtmp=' 0'
245  ELSEIF(ivar.LE.21) THEN
246  IF(idim.GT.kcc) chtmp=' '
247  ELSE
248  IF(idim.GT.kcc) chtmp=' 0'
249  ENDIF
250 
251 C...Length of variable, trailing decimal zeros, quotation marks.
252  llow=1
253  lhig=1
254  DO 240 ll=1,16
255  IF(chtmp(17-ll:17-ll).NE.' ') llow=17-ll
256  IF(chtmp(ll:ll).NE.' ') lhig=ll
257  240 CONTINUE
258  chnew=chtmp(llow:lhig)//' '
259  lnew=1+lhig-llow
260  IF((ivar.GE.5.AND.ivar.LE.8).OR.ivar.EQ.14) THEN
261  lnew=lnew+1
262  250 lnew=lnew-1
263  IF(lnew.GE.2.AND.chnew(lnew:lnew).EQ.'0') goto 250
264  IF(chnew(lnew:lnew).EQ.'.') lnew=lnew-1
265  IF(lnew.EQ.0) THEN
266  chnew(1:3)='0D0'
267  lnew=3
268  ELSE
269  chnew(lnew+1:lnew+2)='D0'
270  lnew=lnew+2
271  ENDIF
272  ELSEIF(ivar.EQ.20.OR.ivar.EQ.21) THEN
273  DO 260 ll=lnew,1,-1
274  IF(chnew(ll:ll).EQ.'''') THEN
275  chtmp=chnew
276  chnew=chtmp(1:ll)//''''//chtmp(ll+1:11)
277  lnew=lnew+1
278  ENDIF
279  260 CONTINUE
280  lnew=min(14,lnew)
281  chtmp=chnew
282  chnew(1:lnew+2)=''''//chtmp(1:lnew)//''''
283  lnew=lnew+2
284  ENDIF
285 
286 C...Form composite character string, often including repetition counter.
287  IF(chnew.NE.chold) THEN
288  nrpt=1
289  chold=chnew
290  chcom=chnew
291  lcom=lnew
292  ELSE
293  lrpt=lnew+1
294  IF(nrpt.GE.2) lrpt=lnew+3
295  IF(nrpt.GE.10) lrpt=lnew+4
296  IF(nrpt.GE.100) lrpt=lnew+5
297  IF(nrpt.GE.1000) lrpt=lnew+6
298  llin=llin-lrpt
299  nrpt=nrpt+1
300  WRITE(chtmp,5400) nrpt
301  lrpt=1
302  IF(nrpt.GE.10) lrpt=2
303  IF(nrpt.GE.100) lrpt=3
304  IF(nrpt.GE.1000) lrpt=4
305  chcom(1:lrpt+1+lnew)=chtmp(17-lrpt:16)//'*'//chnew(1:lnew)
306  lcom=lrpt+1+lnew
307  ENDIF
308 
309 C...Add characters to end of line, to new line (after storing old line),
310 C...or to new block of lines (after writing old block).
311  IF(llin+lcom.LE.70) THEN
312  chlin(llin+1:llin+lcom+1)=chcom(1:lcom)//','
313  llin=llin+lcom+1
314  ELSEIF(nlin.LE.19) THEN
315  chlin(llin+1:72)=' '
316  chblk(nlin)=chlin
317  nlin=nlin+1
318  chlin(6:6+lcom+1)='&'//chcom(1:lcom)//','
319  llin=6+lcom+1
320  ELSE
321  chlin(llin:72)='/'//' '
322  chblk(nlin)=chlin
323  WRITE(chtmp,5400) idim-nrpt
324  chblk(1)(30:33)=chtmp(13:16)
325  DO 270 ilin=1,nlin
326  WRITE(lfn,5700) chblk(ilin)
327  270 CONTINUE
328  nlin=1
329  chlin=' '
330  chlin(7:35+lcom+1)='DATA ('//chvar(ivar)//
331  & ',I= , )/'//chcom(1:lcom)//','
332  WRITE(chtmp,5400) idim-nrpt+1
333  chlin(25:28)=chtmp(13:16)
334  llin=35+lcom+1
335  ENDIF
336  280 CONTINUE
337 
338 C...Write final block of lines.
339  chlin(llin:72)='/'//' '
340  chblk(nlin)=chlin
341  WRITE(chtmp,5400) ndim
342  chblk(1)(30:33)=chtmp(13:16)
343  DO 290 ilin=1,nlin
344  WRITE(lfn,5700) chblk(ilin)
345  290 CONTINUE
346  300 CONTINUE
347  ENDIF
348 
349 C...Formats for reading and writing particle data.
350  5000 FORMAT(1x,i9,2x,a16,2x,a16,3i3,3f12.5,1p,e13.5,2i3)
351  5100 FORMAT(10x,2i5,f12.6,5i10)
352  5200 FORMAT(a120)
353  5300 FORMAT(i9)
354  5400 FORMAT(i16)
355  5500 FORMAT(f16.5)
356  5600 FORMAT(f16.6)
357  5700 FORMAT(a72)
358 
359  RETURN
360  END