Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyedit.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyedit.f
1 
2 C*********************************************************************
3 
4 C...PYEDIT
5 C...Performs global manipulations on the event record, in particular
6 C...to exclude unstable or undetectable partons/particles.
7 
8  SUBROUTINE pyedit(MEDIT)
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...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/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
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/pyctag/nct,mct(4000,2)
22  SAVE /pyjets/,/pydat1/,/pydat2/,/pyctag/
23 C...Local arrays.
24  dimension ns(2),pts(2),pls(2)
25 
26 C...Remove unwanted partons/particles.
27  IF((medit.GE.0.AND.medit.LE.3).OR.medit.EQ.5) THEN
28  imax=n
29  IF(mstu(2).GT.0) imax=mstu(2)
30  i1=max(1,mstu(1))-1
31  DO 110 i=max(1,mstu(1)),imax
32  IF(k(i,1).EQ.0.OR.(k(i,1).GE.21.AND.k(i,1).LE.40)) goto 110
33  IF(medit.EQ.1) THEN
34  IF(k(i,1).GT.10.AND.k(i,1).NE.41.AND.k(i,1).NE.42) goto 110
35  ELSEIF(medit.EQ.2) THEN
36  IF(k(i,1).GT.10.AND.k(i,1).NE.41.AND.k(i,1).NE.42) goto 110
37  kc=pycomp(k(i,2))
38  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
39  & kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
40  & k(i,2).EQ.ksusy1+39) goto 110
41  ELSEIF(medit.EQ.3) THEN
42  IF(k(i,1).GT.10.AND.k(i,1).NE.41.AND.k(i,1).NE.42) goto 110
43  kc=pycomp(k(i,2))
44  IF(kc.EQ.0) goto 110
45  IF(kchg(kc,2).EQ.0.AND.pychge(k(i,2)).EQ.0) goto 110
46  ELSEIF(medit.EQ.5) THEN
47  IF(k(i,1).EQ.13.OR.k(i,1).EQ.14.OR.k(i,1).EQ.52) goto 110
48  kc=pycomp(k(i,2))
49  IF(kc.EQ.0) goto 110
50  IF(k(i,1).GT.10.AND.k(i,1).NE.41.AND.k(i,1).NE.42.AND.
51  & kchg(kc,2).EQ.0) goto 110
52  ENDIF
53 
54 C...Pack remaining partons/particles. Origin no longer known.
55  i1=i1+1
56  DO 100 j=1,5
57  k(i1,j)=k(i,j)
58  p(i1,j)=p(i,j)
59  v(i1,j)=v(i,j)
60  100 CONTINUE
61  k(i1,3)=0
62  110 CONTINUE
63  IF(i1.LT.n) mstu(3)=0
64  IF(i1.LT.n) mstu(70)=0
65  n=i1
66 
67 C...Selective removal of class of entries. New position of retained.
68  ELSEIF(medit.GE.11.AND.medit.LE.15) THEN
69  i1=0
70  DO 120 i=1,n
71  k(i,3)=mod(k(i,3),mstu(5))
72  IF(medit.EQ.11.AND.k(i,1).LT.0) goto 120
73  IF(medit.EQ.12.AND.k(i,1).EQ.0) goto 120
74  IF(medit.EQ.13.AND.(k(i,1).EQ.11.OR.k(i,1).EQ.12.OR.
75  & k(i,1).EQ.15.OR.k(i,1).EQ.51).AND.k(i,2).NE.94) goto 120
76  IF(medit.EQ.14.AND.(k(i,1).EQ.13.OR.k(i,1).EQ.14.OR.
77  & k(i,1).EQ.52.OR.k(i,2).EQ.94)) goto 120
78  IF(medit.EQ.15.AND.k(i,1).GE.21.AND.k(i,1).LE.40) goto 120
79  i1=i1+1
80  k(i,3)=k(i,3)+mstu(5)*i1
81  120 CONTINUE
82 
83 C...Find new event history information and replace old.
84  DO 140 i=1,n
85  IF(k(i,1).LE.0.OR.(k(i,1).GE.21.AND.k(i,1).LE.40).OR.
86  & k(i,3)/mstu(5).EQ.0) goto 140
87  id=i
88  130 im=mod(k(id,3),mstu(5))
89  IF(medit.EQ.13.AND.im.GT.0.AND.im.LE.n) THEN
90  IF((k(im,1).EQ.11.OR.k(im,1).EQ.12.OR.k(im,1).EQ.15.OR.
91  & k(im,1).EQ.51).AND.k(im,2).NE.94) THEN
92  id=im
93  goto 130
94  ENDIF
95  ELSEIF(medit.EQ.14.AND.im.GT.0.AND.im.LE.n) THEN
96  IF(k(im,1).EQ.13.OR.k(im,1).EQ.14.OR.k(im,1).EQ.52.OR.
97  & k(im,2).EQ.94) THEN
98  id=im
99  goto 130
100  ENDIF
101  ENDIF
102  k(i,3)=mstu(5)*(k(i,3)/mstu(5))
103  IF(im.NE.0) k(i,3)=k(i,3)+k(im,3)/mstu(5)
104  IF(k(i,1).NE.3.AND.k(i,1).NE.13.AND.k(i,1).NE.14.AND.
105  & k(i,1).NE.42.AND.k(i,1).NE.52) THEN
106  IF(k(i,4).GT.0.AND.k(i,4).LE.mstu(4)) k(i,4)=
107  & k(k(i,4),3)/mstu(5)
108  IF(k(i,5).GT.0.AND.k(i,5).LE.mstu(4)) k(i,5)=
109  & k(k(i,5),3)/mstu(5)
110  ELSE
111  kcm=mod(k(i,4)/mstu(5),mstu(5))
112  IF(kcm.GT.0.AND.kcm.LE.mstu(4).AND.k(i,1).NE.42.AND.
113  & k(i,1).NE.52) kcm=k(kcm,3)/mstu(5)
114  kcd=mod(k(i,4),mstu(5))
115  IF(kcd.GT.0.AND.kcd.LE.mstu(4)) kcd=k(kcd,3)/mstu(5)
116  k(i,4)=mstu(5)**2*(k(i,4)/mstu(5)**2)+mstu(5)*kcm+kcd
117  kcm=mod(k(i,5)/mstu(5),mstu(5))
118  IF(kcm.GT.0.AND.kcm.LE.mstu(4)) kcm=k(kcm,3)/mstu(5)
119  kcd=mod(k(i,5),mstu(5))
120  IF(kcd.GT.0.AND.kcd.LE.mstu(4)) kcd=k(kcd,3)/mstu(5)
121  k(i,5)=mstu(5)**2*(k(i,5)/mstu(5)**2)+mstu(5)*kcm+kcd
122  ENDIF
123  140 CONTINUE
124 
125 C...Pack remaining entries.
126  i1=0
127  mstu90=mstu(90)
128  mstu(90)=0
129  DO 170 i=1,n
130  IF(k(i,3)/mstu(5).EQ.0) goto 170
131  i1=i1+1
132  DO 150 j=1,5
133  k(i1,j)=k(i,j)
134  p(i1,j)=p(i,j)
135  v(i1,j)=v(i,j)
136  150 CONTINUE
137 C...Also update LHA1 colour tags
138  mct(i1,1)=mct(i,1)
139  mct(i1,2)=mct(i,2)
140  k(i1,3)=mod(k(i1,3),mstu(5))
141  DO 160 iz=1,mstu90
142  IF(i.EQ.mstu(90+iz)) THEN
143  mstu(90)=mstu(90)+1
144  mstu(90+mstu(90))=i1
145  paru(90+mstu(90))=paru(90+iz)
146  ENDIF
147  160 CONTINUE
148  170 CONTINUE
149  IF(i1.LT.n) mstu(3)=0
150  IF(i1.LT.n) mstu(70)=0
151  n=i1
152 
153 C...Fill in some missing daughter pointers (lost in colour flow).
154  ELSEIF(medit.EQ.16) THEN
155  DO 220 i=1,n
156  IF(k(i,1).LE.10.OR.(k(i,1).GE.21.AND.k(i,1).LE.50)) goto 220
157  IF(k(i,4).NE.0.OR.k(i,5).NE.0) goto 220
158 C...Find daughters who point to mother.
159  DO 180 i1=i+1,n
160  IF(k(i1,3).NE.i) THEN
161  ELSEIF(k(i,4).EQ.0) THEN
162  k(i,4)=i1
163  ELSE
164  k(i,5)=i1
165  ENDIF
166  180 CONTINUE
167  IF(k(i,5).EQ.0) k(i,5)=k(i,4)
168  IF(k(i,4).NE.0) goto 220
169 C...Find daughters who point to documentation version of mother.
170  im=k(i,3)
171  IF(im.LE.0.OR.im.GE.i) goto 220
172  IF(k(im,1).LE.20.OR.k(im,1).GT.30) goto 220
173  IF(k(im,2).NE.k(i,2).OR.abs(p(im,5)-p(i,5)).GT.1d-2) goto 220
174  DO 190 i1=i+1,n
175  IF(k(i1,3).NE.im) THEN
176  ELSEIF(k(i,4).EQ.0) THEN
177  k(i,4)=i1
178  ELSE
179  k(i,5)=i1
180  ENDIF
181  190 CONTINUE
182  IF(k(i,5).EQ.0) k(i,5)=k(i,4)
183  IF(k(i,4).NE.0) goto 220
184 C...Find daughters who point to documentation daughters who,
185 C...in their turn, point to documentation mother.
186  id1=im
187  id2=im
188  DO 200 i1=im+1,i-1
189  IF(k(i1,3).EQ.im.AND.k(i1,1).GE.21.AND.k(i1,1).LE.30) THEN
190  id2=i1
191  IF(id1.EQ.im) id1=i1
192  ENDIF
193  200 CONTINUE
194  DO 210 i1=i+1,n
195  IF(k(i1,3).NE.id1.AND.k(i1,3).NE.id2) THEN
196  ELSEIF(k(i,4).EQ.0) THEN
197  k(i,4)=i1
198  ELSE
199  k(i,5)=i1
200  ENDIF
201  210 CONTINUE
202  IF(k(i,5).EQ.0) k(i,5)=k(i,4)
203  220 CONTINUE
204 
205 C...Save top entries at bottom of PYJETS commonblock.
206  ELSEIF(medit.EQ.21) THEN
207  IF(2*n.GE.mstu(4)) THEN
208  CALL pyerrm(11,'(PYEDIT:) no more memory left in PYJETS')
209  RETURN
210  ENDIF
211  DO 240 i=1,n
212  DO 230 j=1,5
213  k(mstu(4)-i,j)=k(i,j)
214  p(mstu(4)-i,j)=p(i,j)
215  v(mstu(4)-i,j)=v(i,j)
216  230 CONTINUE
217  240 CONTINUE
218  mstu(32)=n
219 
220 C...Restore bottom entries of commonblock PYJETS to top.
221  ELSEIF(medit.EQ.22) THEN
222  DO 260 i=1,mstu(32)
223  DO 250 j=1,5
224  k(i,j)=k(mstu(4)-i,j)
225  p(i,j)=p(mstu(4)-i,j)
226  v(i,j)=v(mstu(4)-i,j)
227  250 CONTINUE
228  260 CONTINUE
229  n=mstu(32)
230 
231 C...Mark primary entries at top of commonblock PYJETS as untreated.
232  ELSEIF(medit.EQ.23) THEN
233  i1=0
234  DO 270 i=1,n
235  kh=k(i,3)
236  IF(kh.GE.1) THEN
237  IF(k(kh,1).GE.21.AND.k(kh,1).LE.30) kh=0
238  ENDIF
239  IF(kh.NE.0) goto 280
240  i1=i1+1
241  IF(k(i,1).GE.11.AND.k(i,1).LE.20) k(i,1)=k(i,1)-10
242  IF(k(i,1).GE.51.AND.k(i,1).LE.60) k(i,1)=k(i,1)-10
243  270 CONTINUE
244  280 n=i1
245 
246 C...Place largest axis along z axis and second largest in xy plane.
247  ELSEIF(medit.EQ.31.OR.medit.EQ.32) THEN
248  CALL pyrobo(1,n+mstu(3),0d0,-pyangl(p(mstu(61),1),
249  & p(mstu(61),2)),0d0,0d0,0d0)
250  CALL pyrobo(1,n+mstu(3),-pyangl(p(mstu(61),3),
251  & p(mstu(61),1)),0d0,0d0,0d0,0d0)
252  CALL pyrobo(1,n+mstu(3),0d0,-pyangl(p(mstu(61)+1,1),
253  & p(mstu(61)+1,2)),0d0,0d0,0d0)
254  IF(medit.EQ.31) RETURN
255 
256 C...Rotate to put slim jet along +z axis.
257  DO 290 is=1,2
258  ns(is)=0
259  pts(is)=0d0
260  pls(is)=0d0
261  290 CONTINUE
262  DO 300 i=1,n
263  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 300
264  IF(mstu(41).GE.2) THEN
265  kc=pycomp(k(i,2))
266  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
267  & kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
268  & k(i,2).EQ.ksusy1+39) goto 300
269  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.pychge(k(i,2))
270  & .EQ.0) goto 300
271  ENDIF
272  is=2d0-sign(0.5d0,p(i,3))
273  ns(is)=ns(is)+1
274  pts(is)=pts(is)+sqrt(p(i,1)**2+p(i,2)**2)
275  300 CONTINUE
276  IF(ns(1)*pts(2)**2.LT.ns(2)*pts(1)**2)
277  & CALL pyrobo(1,n+mstu(3),paru(1),0d0,0d0,0d0,0d0)
278 
279 C...Rotate to put second largest jet into -z,+x quadrant.
280  DO 310 i=1,n
281  IF(p(i,3).GE.0d0) goto 310
282  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 310
283  IF(mstu(41).GE.2) THEN
284  kc=pycomp(k(i,2))
285  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
286  & kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
287  & k(i,2).EQ.ksusy1+39) goto 310
288  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.pychge(k(i,2))
289  & .EQ.0) goto 310
290  ENDIF
291  is=2d0-sign(0.5d0,p(i,1))
292  pls(is)=pls(is)-p(i,3)
293  310 CONTINUE
294  IF(pls(2).GT.pls(1)) CALL pyrobo(1,n+mstu(3),0d0,paru(1),
295  & 0d0,0d0,0d0)
296  ENDIF
297 
298  RETURN
299  END