Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyindf.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyindf.f
1 
2 C*********************************************************************
3 
4 C...PYINDF
5 C...Handles the fragmentation of a jet system (or a single
6 C...jet) according to independent fragmentation models.
7 
8  SUBROUTINE pyindf(IP)
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/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
16  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
17  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
18  SAVE /pyjets/,/pydat1/,/pydat2/
19 C...Local arrays.
20  dimension dps(5),psi(4),nfi(3),nfl(3),ifet(3),kflf(3),
21  &kflo(2),pxo(2),pyo(2),wo(2)
22 
23 C.. MOPS error message
24  IF(mstj(12).GT.3) CALL pyerrm(9,'(PYINDF:) MSTJ(12)>3 options'//
25  &' are not treated as expected in independent fragmentation')
26 
27 C...Reset counters. Identify parton system and take copy. Check flavour.
28  nsav=n
29  mstu90=mstu(90)
30  njet=0
31  kqsum=0
32  DO 100 j=1,5
33  dps(j)=0d0
34  100 CONTINUE
35  i=ip-1
36  110 i=i+1
37  IF(i.GT.min(n,mstu(4)-mstu(32))) THEN
38  CALL pyerrm(12,'(PYINDF:) failed to reconstruct jet system')
39  IF(mstu(21).GE.1) RETURN
40  ENDIF
41  IF(k(i,1).NE.1.AND.k(i,1).NE.2) goto 110
42  kc=pycomp(k(i,2))
43  IF(kc.EQ.0) goto 110
44  kq=kchg(kc,2)*isign(1,k(i,2))
45  IF(kq.EQ.0) goto 110
46  njet=njet+1
47  IF(kq.NE.2) kqsum=kqsum+kq
48  DO 120 j=1,5
49  k(nsav+njet,j)=k(i,j)
50  p(nsav+njet,j)=p(i,j)
51  dps(j)=dps(j)+p(i,j)
52  120 CONTINUE
53  k(nsav+njet,3)=i
54  IF(k(i,1).EQ.2.OR.(mstj(3).LE.5.AND.n.GT.i.AND.
55  &k(i+1,1).EQ.2)) goto 110
56  IF(njet.NE.1.AND.kqsum.NE.0) THEN
57  CALL pyerrm(12,'(PYINDF:) unphysical flavour combination')
58  IF(mstu(21).GE.1) RETURN
59  ENDIF
60 
61 C...Boost copied system to CM frame. Find CM energy and sum flavours.
62  IF(njet.NE.1) THEN
63  mstu(33)=1
64  CALL pyrobo(nsav+1,nsav+njet,0d0,0d0,-dps(1)/dps(4),
65  & -dps(2)/dps(4),-dps(3)/dps(4))
66  ENDIF
67  pecm=0d0
68  DO 130 j=1,3
69  nfi(j)=0
70  130 CONTINUE
71  DO 140 i=nsav+1,nsav+njet
72  pecm=pecm+p(i,4)
73  kfa=iabs(k(i,2))
74  IF(kfa.LE.3) THEN
75  nfi(kfa)=nfi(kfa)+isign(1,k(i,2))
76  ELSEIF(kfa.GT.1000) THEN
77  kfla=mod(kfa/1000,10)
78  kflb=mod(kfa/100,10)
79  IF(kfla.LE.3) nfi(kfla)=nfi(kfla)+isign(1,k(i,2))
80  IF(kflb.LE.3) nfi(kflb)=nfi(kflb)+isign(1,k(i,2))
81  ENDIF
82  140 CONTINUE
83 
84 C...Loop over attempts made. Reset counters.
85  ntry=0
86  150 ntry=ntry+1
87  IF(ntry.GT.200) THEN
88  CALL pyerrm(14,'(PYINDF:) caught in infinite loop')
89  IF(mstu(21).GE.1) RETURN
90  ENDIF
91  n=nsav+njet
92  mstu(90)=mstu90
93  DO 160 j=1,3
94  nfl(j)=nfi(j)
95  ifet(j)=0
96  kflf(j)=0
97  160 CONTINUE
98 
99 C...Loop over jets to be fragmented.
100  DO 230 ip1=nsav+1,nsav+njet
101  mstj(91)=0
102  nsav1=n
103  mstu91=mstu(90)
104 
105 C...Initial flavour and momentum values. Jet along +z axis.
106  kflh=iabs(k(ip1,2))
107  IF(kflh.GT.10) kflh=mod(kflh/1000,10)
108  kflo(2)=0
109  wf=p(ip1,4)+sqrt(p(ip1,1)**2+p(ip1,2)**2+p(ip1,3)**2)
110 
111 C...Initial values for quark or diquark jet.
112  170 IF(iabs(k(ip1,2)).NE.21) THEN
113  nstr=1
114  kflo(1)=k(ip1,2)
115  CALL pyptdi(0,pxo(1),pyo(1))
116  wo(1)=wf
117 
118 C...Initial values for gluon treated like random quark jet.
119  ELSEIF(mstj(2).LE.2) THEN
120  nstr=1
121  IF(mstj(2).EQ.2) mstj(91)=1
122  kflo(1)=int(1d0+(2d0+parj(2))*pyr(0))*(-1)**int(pyr(0)+0.5d0)
123  CALL pyptdi(0,pxo(1),pyo(1))
124  wo(1)=wf
125 
126 C...Initial values for gluon treated like quark-antiquark jet pair,
127 C...sharing energy according to Altarelli-Parisi splitting function.
128  ELSE
129  nstr=2
130  IF(mstj(2).EQ.4) mstj(91)=1
131  kflo(1)=int(1d0+(2d0+parj(2))*pyr(0))*(-1)**int(pyr(0)+0.5d0)
132  kflo(2)=-kflo(1)
133  CALL pyptdi(0,pxo(1),pyo(1))
134  pxo(2)=-pxo(1)
135  pyo(2)=-pyo(1)
136  wo(1)=wf*pyr(0)**(1d0/3d0)
137  wo(2)=wf-wo(1)
138  ENDIF
139 
140 C...Initial values for rank, flavour, pT and W+.
141  DO 220 istr=1,nstr
142  180 i=n
143  mstu(90)=mstu91
144  irank=0
145  kfl1=kflo(istr)
146  px1=pxo(istr)
147  py1=pyo(istr)
148  w=wo(istr)
149 
150 C...New hadron. Generate flavour and hadron species.
151  190 i=i+1
152  IF(i.GE.mstu(4)-mstu(32)-njet-5) THEN
153  CALL pyerrm(11,'(PYINDF:) no more memory left in PYJETS')
154  IF(mstu(21).GE.1) RETURN
155  ENDIF
156  irank=irank+1
157  k(i,1)=1
158  k(i,3)=ip1
159  k(i,4)=0
160  k(i,5)=0
161  200 CALL pykfdi(kfl1,0,kfl2,k(i,2))
162  IF(k(i,2).EQ.0) goto 180
163  IF(irank.EQ.1.AND.iabs(kfl1).LE.10.AND.iabs(kfl2).GT.10) THEN
164  IF(pyr(0).GT.parj(19)) goto 200
165  ENDIF
166 
167 C...Find hadron mass. Generate four-momentum.
168  p(i,5)=pymass(k(i,2))
169  CALL pyptdi(kfl1,px2,py2)
170  p(i,1)=px1+px2
171  p(i,2)=py1+py2
172  pr=p(i,5)**2+p(i,1)**2+p(i,2)**2
173  CALL pyzdis(kfl1,kfl2,pr,z)
174  mzsav=0
175  IF(iabs(kfl1).GE.4.AND.iabs(kfl1).LE.8.AND.mstu(90).LT.8) THEN
176  mzsav=1
177  mstu(90)=mstu(90)+1
178  mstu(90+mstu(90))=i
179  paru(90+mstu(90))=z
180  ENDIF
181  p(i,3)=0.5d0*(z*w-pr/max(1d-4,z*w))
182  p(i,4)=0.5d0*(z*w+pr/max(1d-4,z*w))
183  IF(mstj(3).GE.1.AND.irank.EQ.1.AND.kflh.GE.4.AND.
184  & p(i,3).LE.0.001d0) THEN
185  IF(w.GE.p(i,5)+0.5d0*parj(32)) goto 180
186  p(i,3)=0.0001d0
187  p(i,4)=sqrt(pr)
188  z=p(i,4)/w
189  ENDIF
190 
191 C...Remaining flavour and momentum.
192  kfl1=-kfl2
193  px1=-px2
194  py1=-py2
195  w=(1d0-z)*w
196  DO 210 j=1,5
197  v(i,j)=0d0
198  210 CONTINUE
199 
200 C...Check if pL acceptable. Go back for new hadron if enough energy.
201  IF(mstj(3).GE.0.AND.p(i,3).LT.0d0) THEN
202  i=i-1
203  IF(mzsav.EQ.1) mstu(90)=mstu(90)-1
204  ENDIF
205  IF(w.GT.parj(31)) goto 190
206  n=i
207  220 CONTINUE
208  IF(mod(mstj(3),5).EQ.4.AND.n.EQ.nsav1) wf=wf+0.1d0*parj(32)
209  IF(mod(mstj(3),5).EQ.4.AND.n.EQ.nsav1) goto 170
210 
211 C...Rotate jet to new direction.
212  the=pyangl(p(ip1,3),sqrt(p(ip1,1)**2+p(ip1,2)**2))
213  phi=pyangl(p(ip1,1),p(ip1,2))
214  mstu(33)=1
215  CALL pyrobo(nsav1+1,n,the,phi,0d0,0d0,0d0)
216  k(k(ip1,3),4)=nsav1+1
217  k(k(ip1,3),5)=n
218 
219 C...End of jet generation loop. Skip conservation in some cases.
220  230 CONTINUE
221  IF(njet.EQ.1.OR.mstj(3).LE.0) goto 490
222  IF(mod(mstj(3),5).NE.0.AND.n-nsav-njet.LT.2) goto 150
223 
224 C...Subtract off produced hadron flavours, finished if zero.
225  DO 240 i=nsav+njet+1,n
226  kfa=iabs(k(i,2))
227  kfla=mod(kfa/1000,10)
228  kflb=mod(kfa/100,10)
229  kflc=mod(kfa/10,10)
230  IF(kfla.EQ.0) THEN
231  IF(kflb.LE.3) nfl(kflb)=nfl(kflb)-isign(1,k(i,2))*(-1)**kflb
232  IF(kflc.LE.3) nfl(kflc)=nfl(kflc)+isign(1,k(i,2))*(-1)**kflb
233  ELSE
234  IF(kfla.LE.3) nfl(kfla)=nfl(kfla)-isign(1,k(i,2))
235  IF(kflb.LE.3) nfl(kflb)=nfl(kflb)-isign(1,k(i,2))
236  IF(kflc.LE.3) nfl(kflc)=nfl(kflc)-isign(1,k(i,2))
237  ENDIF
238  240 CONTINUE
239  nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+
240  &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3
241  IF(nreq.EQ.0) goto 320
242 
243 C...Take away flavour of low-momentum particles until enough freedom.
244  nrem=0
245  250 irem=0
246  p2min=pecm**2
247  DO 260 i=nsav+njet+1,n
248  p2=p(i,1)**2+p(i,2)**2+p(i,3)**2
249  IF(k(i,1).EQ.1.AND.p2.LT.p2min) irem=i
250  IF(k(i,1).EQ.1.AND.p2.LT.p2min) p2min=p2
251  260 CONTINUE
252  IF(irem.EQ.0) goto 150
253  k(irem,1)=7
254  kfa=iabs(k(irem,2))
255  kfla=mod(kfa/1000,10)
256  kflb=mod(kfa/100,10)
257  kflc=mod(kfa/10,10)
258  IF(kfla.GE.4.OR.kflb.GE.4) k(irem,1)=8
259  IF(k(irem,1).EQ.8) goto 250
260  IF(kfla.EQ.0) THEN
261  isgn=isign(1,k(irem,2))*(-1)**kflb
262  IF(kflb.LE.3) nfl(kflb)=nfl(kflb)+isgn
263  IF(kflc.LE.3) nfl(kflc)=nfl(kflc)-isgn
264  ELSE
265  IF(kfla.LE.3) nfl(kfla)=nfl(kfla)+isign(1,k(irem,2))
266  IF(kflb.LE.3) nfl(kflb)=nfl(kflb)+isign(1,k(irem,2))
267  IF(kflc.LE.3) nfl(kflc)=nfl(kflc)+isign(1,k(irem,2))
268  ENDIF
269  nrem=nrem+1
270  nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+
271  &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3
272  IF(nreq.GT.nrem) goto 250
273  DO 270 i=nsav+njet+1,n
274  IF(k(i,1).EQ.8) k(i,1)=1
275  270 CONTINUE
276 
277 C...Find combination of existing and new flavours for hadron.
278  280 nfet=2
279  IF(nfl(1)+nfl(2)+nfl(3).NE.0) nfet=3
280  IF(nreq.LT.nrem) nfet=1
281  IF(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3)).EQ.0) nfet=0
282  DO 290 j=1,nfet
283  ifet(j)=1+(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3)))*pyr(0)
284  kflf(j)=isign(1,nfl(1))
285  IF(ifet(j).GT.iabs(nfl(1))) kflf(j)=isign(2,nfl(2))
286  IF(ifet(j).GT.iabs(nfl(1))+iabs(nfl(2))) kflf(j)=isign(3,nfl(3))
287  290 CONTINUE
288  IF(nfet.EQ.2.AND.(ifet(1).EQ.ifet(2).OR.kflf(1)*kflf(2).GT.0))
289  &goto 280
290  IF(nfet.EQ.3.AND.(ifet(1).EQ.ifet(2).OR.ifet(1).EQ.ifet(3).OR.
291  &ifet(2).EQ.ifet(3).OR.kflf(1)*kflf(2).LT.0.OR.kflf(1)*kflf(3)
292  &.LT.0.OR.kflf(1)*(nfl(1)+nfl(2)+nfl(3)).LT.0)) goto 280
293  IF(nfet.EQ.0) kflf(1)=1+int((2d0+parj(2))*pyr(0))
294  IF(nfet.EQ.0) kflf(2)=-kflf(1)
295  IF(nfet.EQ.1) kflf(2)=isign(1+int((2d0+parj(2))*pyr(0)),-kflf(1))
296  IF(nfet.LE.2) kflf(3)=0
297  IF(kflf(3).NE.0) THEN
298  kflfc=isign(1000*max(iabs(kflf(1)),iabs(kflf(3)))+
299  & 100*min(iabs(kflf(1)),iabs(kflf(3)))+1,kflf(1))
300  IF(kflf(1).EQ.kflf(3).OR.(1d0+3d0*parj(4))*pyr(0).GT.1d0)
301  & kflfc=kflfc+isign(2,kflfc)
302  ELSE
303  kflfc=kflf(1)
304  ENDIF
305  CALL pykfdi(kflfc,kflf(2),kfldmp,kf)
306  IF(kf.EQ.0) goto 280
307  DO 300 j=1,max(2,nfet)
308  nfl(iabs(kflf(j)))=nfl(iabs(kflf(j)))-isign(1,kflf(j))
309  300 CONTINUE
310 
311 C...Store hadron at random among free positions.
312  npos=min(1+int(pyr(0)*nrem),nrem)
313  DO 310 i=nsav+njet+1,n
314  IF(k(i,1).EQ.7) npos=npos-1
315  IF(k(i,1).EQ.1.OR.npos.NE.0) goto 310
316  k(i,1)=1
317  k(i,2)=kf
318  p(i,5)=pymass(k(i,2))
319  p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)
320  310 CONTINUE
321  nrem=nrem-1
322  nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+
323  &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3
324  IF(nrem.GT.0) goto 280
325 
326 C...Compensate for missing momentum in global scheme (3 options).
327  320 IF(mod(mstj(3),5).NE.0.AND.mod(mstj(3),5).NE.4) THEN
328  DO 340 j=1,3
329  psi(j)=0d0
330  DO 330 i=nsav+njet+1,n
331  psi(j)=psi(j)+p(i,j)
332  330 CONTINUE
333  340 CONTINUE
334  psi(4)=psi(1)**2+psi(2)**2+psi(3)**2
335  pws=0d0
336  DO 350 i=nsav+njet+1,n
337  IF(mod(mstj(3),5).EQ.1) pws=pws+p(i,4)
338  IF(mod(mstj(3),5).EQ.2) pws=pws+sqrt(p(i,5)**2+(psi(1)*p(i,1)+
339  & psi(2)*p(i,2)+psi(3)*p(i,3))**2/psi(4))
340  IF(mod(mstj(3),5).EQ.3) pws=pws+1d0
341  350 CONTINUE
342  DO 370 i=nsav+njet+1,n
343  IF(mod(mstj(3),5).EQ.1) pw=p(i,4)
344  IF(mod(mstj(3),5).EQ.2) pw=sqrt(p(i,5)**2+(psi(1)*p(i,1)+
345  & psi(2)*p(i,2)+psi(3)*p(i,3))**2/psi(4))
346  IF(mod(mstj(3),5).EQ.3) pw=1d0
347  DO 360 j=1,3
348  p(i,j)=p(i,j)-psi(j)*pw/pws
349  360 CONTINUE
350  p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)
351  370 CONTINUE
352 
353 C...Compensate for missing momentum withing each jet separately.
354  ELSEIF(mod(mstj(3),5).EQ.4) THEN
355  DO 390 i=n+1,n+njet
356  k(i,1)=0
357  DO 380 j=1,5
358  p(i,j)=0d0
359  380 CONTINUE
360  390 CONTINUE
361  DO 410 i=nsav+njet+1,n
362  ir1=k(i,3)
363  ir2=n+ir1-nsav
364  k(ir2,1)=k(ir2,1)+1
365  pls=(p(i,1)*p(ir1,1)+p(i,2)*p(ir1,2)+p(i,3)*p(ir1,3))/
366  & (p(ir1,1)**2+p(ir1,2)**2+p(ir1,3)**2)
367  DO 400 j=1,3
368  p(ir2,j)=p(ir2,j)+p(i,j)-pls*p(ir1,j)
369  400 CONTINUE
370  p(ir2,4)=p(ir2,4)+p(i,4)
371  p(ir2,5)=p(ir2,5)+pls
372  410 CONTINUE
373  pss=0d0
374  DO 420 i=n+1,n+njet
375  IF(k(i,1).NE.0) pss=pss+p(i,4)/(pecm*(0.8d0*p(i,5)+0.2d0))
376  420 CONTINUE
377  DO 440 i=nsav+njet+1,n
378  ir1=k(i,3)
379  ir2=n+ir1-nsav
380  pls=(p(i,1)*p(ir1,1)+p(i,2)*p(ir1,2)+p(i,3)*p(ir1,3))/
381  & (p(ir1,1)**2+p(ir1,2)**2+p(ir1,3)**2)
382  DO 430 j=1,3
383  p(i,j)=p(i,j)-p(ir2,j)/k(ir2,1)+(1d0/(p(ir2,5)*pss)-1d0)*
384  & pls*p(ir1,j)
385  430 CONTINUE
386  p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)
387  440 CONTINUE
388  ENDIF
389 
390 C...Scale momenta for energy conservation.
391  IF(mod(mstj(3),5).NE.0) THEN
392  pms=0d0
393  pes=0d0
394  pqs=0d0
395  DO 450 i=nsav+njet+1,n
396  pms=pms+p(i,5)
397  pes=pes+p(i,4)
398  pqs=pqs+p(i,5)**2/p(i,4)
399  450 CONTINUE
400  IF(pms.GE.pecm) goto 150
401  neco=0
402  460 neco=neco+1
403  pfac=(pecm-pqs)/(pes-pqs)
404  pes=0d0
405  pqs=0d0
406  DO 480 i=nsav+njet+1,n
407  DO 470 j=1,3
408  p(i,j)=pfac*p(i,j)
409  470 CONTINUE
410  p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)
411  pes=pes+p(i,4)
412  pqs=pqs+p(i,5)**2/p(i,4)
413  480 CONTINUE
414  IF(neco.LT.10.AND.abs(pecm-pes).GT.2d-6*pecm) goto 460
415  ENDIF
416 
417 C...Origin of produced particles and parton daughter pointers.
418  490 DO 500 i=nsav+njet+1,n
419  IF(mstu(16).NE.2) k(i,3)=nsav+1
420  IF(mstu(16).EQ.2) k(i,3)=k(k(i,3),3)
421  500 CONTINUE
422  DO 510 i=nsav+1,nsav+njet
423  i1=k(i,3)
424  k(i1,1)=k(i1,1)+10
425  IF(mstu(16).NE.2) THEN
426  k(i1,4)=nsav+1
427  k(i1,5)=nsav+1
428  ELSE
429  k(i1,4)=k(i1,4)-njet+1
430  k(i1,5)=k(i1,5)-njet+1
431  IF(k(i1,5).LT.k(i1,4)) THEN
432  k(i1,4)=0
433  k(i1,5)=0
434  ENDIF
435  ENDIF
436  510 CONTINUE
437 
438 C...Document independent fragmentation system. Remove copy of jets.
439  nsav=nsav+1
440  k(nsav,1)=11
441  k(nsav,2)=93
442  k(nsav,3)=ip
443  k(nsav,4)=nsav+1
444  k(nsav,5)=n-njet+1
445  DO 520 j=1,4
446  p(nsav,j)=dps(j)
447  v(nsav,j)=v(ip,j)
448  520 CONTINUE
449  p(nsav,5)=sqrt(max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))
450  v(nsav,5)=0d0
451  DO 540 i=nsav+njet,n
452  DO 530 j=1,5
453  k(i-njet+1,j)=k(i,j)
454  p(i-njet+1,j)=p(i,j)
455  v(i-njet+1,j)=v(i,j)
456  530 CONTINUE
457  540 CONTINUE
458  n=n-njet+1
459  DO 550 iz=mstu90+1,mstu(90)
460  mstu(90+iz)=mstu(90+iz)-njet+1
461  550 CONTINUE
462 
463 C...Boost back particle system. Set production vertices.
464  IF(njet.NE.1) CALL pyrobo(nsav+1,n,0d0,0d0,dps(1)/dps(4),
465  &dps(2)/dps(4),dps(3)/dps(4))
466  DO 570 i=nsav+1,n
467  DO 560 j=1,4
468  v(i,j)=v(ip,j)
469  560 CONTINUE
470  570 CONTINUE
471 
472  RETURN
473  END