Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyevol.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyevol.f
1 
2 C***********************************************************************
3 
4 C...PYEVOL
5 C...Handles intertwined pT-ordered spacelike initial-state parton
6 C...and multiple interactions.
7 
8  SUBROUTINE pyevol(MODE,PT2MAX,PT2MIN)
9 C...Mode = -1 : Initialize first time. Determine MAX and MIN scales.
10 C...MODE = 0 : (Re-)initialize ISR/MI evolution.
11 C...Mode = 1 : Evolve event from PT2MAX to PT2MIN.
12 
13 C...Double precision and integer declarations.
14  IMPLICIT DOUBLE PRECISION(a-h, o-z)
15  IMPLICIT INTEGER(i-n)
16  INTEGER pyk,pychge,pycomp
17 C...External
18  EXTERNAL pyalps
19  DOUBLE PRECISION pyalps
20 C...Parameter statement for maximum size of showers.
21  parameter(maxnur=1000)
22 C...Commonblocks.
23  common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
24  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
25  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
26  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
27  common/pypars/mstp(200),parp(200),msti(200),pari(200)
28  common/pyint1/mint(400),vint(400)
29  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
30  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
31  common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
32  & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
33  & xmi(2,240),pt2mi(240),imisep(0:240)
34  common/pyctag/nct,mct(4000,2)
35  common/pyismx/mimx,jsmx,kflamx,kflcmx,kfbeam(2),nisgen(2,240),
36  & pt2mx,pt2amx,zmx,rm2cmx,q2bmx,phimx
37  common/pyisjn/mjn1mx,mjn2mx,mjoind(2,240)
38 C...Local arrays and saved variables.
39  dimension vintsv(11:80),ksav(4,5),psav(4,5),vsav(4,5),shat(240)
40  SAVE nsav,nparts,m15sv,m16sv,m21sv,m22sv,vintsv,shat,isubhd,alam3
41  & ,psav,ksav,vsav
42 
43  SAVE /pypart/,/pyjets/,/pydat1/,/pydat2/,/pypars/,/pyint1/,
44  & /pyint2/,/pyint3/,/pyintm/,/pyctag/,/pyismx/,/pyisjn/
45 
46 C----------------------------------------------------------------------
47 C...MODE=-1: Pre-initialization. Store info on hard scattering etc,
48 C...done only once per event, while MODE=0 is repeated each time the
49 C...evolution needs to be restarted.
50  IF (mode.EQ.-1) THEN
51  isubhd=mint(1)
52  nsav=n
53  nparts=npart
54 C...Store hard scattering variables
55  m15sv=mint(15)
56  m16sv=mint(16)
57  m21sv=mint(21)
58  m22sv=mint(22)
59  DO 100 j=11,80
60  vintsv(j)=vint(j)
61  100 CONTINUE
62  DO 120 j=1,5
63  DO 110 is=1,4
64  i=is+mint(84)
65  psav(is,j)=p(i,j)
66  ksav(is,j)=k(i,j)
67  vsav(is,j)=v(i,j)
68  110 CONTINUE
69  120 CONTINUE
70 
71 C...Set shat for hardest scattering
72  shat(1)=vint(44)
73  IF(iset(isubhd).GE.3.AND.iset(isubhd).LE.5) shat(1)=vint(26)
74  & *vint(2)
75 
76 C...Compute 3-Flavour Lambda_QCD (sets absolute lowest PT scale below)
77  rmc=pmas(4,1)
78  rmb=pmas(5,1)
79  alam4=parp(61)
80  IF(mstu(112).LT.4) alam4=parp(61)*(parp(61)/rmc)**(2d0/25d0)
81  IF(mstu(112).GT.4) alam4=parp(61)*(rmb/parp(61))**(2d0/25d0)
82  alam3=alam4*(rmc/alam4)**(2d0/27d0)
83 
84 C----------------------------------------------------------------------
85 C...MODE= 0: Initialize ISR/MI evolution, i.e. begin from hardest
86 C...interaction initiators, with no previous evolution. Check the input
87 C...PT2MAX and PT2MIN and impose extra constraints on minimum PT2 (e.g.
88 C...must be larger than Lambda_QCD) and maximum PT2 (e.g. must be
89 C...smaller than the CM energy / 2.)
90  ELSEIF (mode.EQ.0) THEN
91 C...Reset counters and switches
92  n=nsav
93  npart=nparts
94  mint(30)=0
95  mint(31)=1
96  mint(36)=1
97 C...Reset hard scattering variables
98  mint(1)=isubhd
99  DO 130 j=11,80
100  vint(j)=vintsv(j)
101  130 CONTINUE
102  DO 150 j=1,5
103  DO 140 is=1,4
104  i=is+mint(84)
105  p(i,j)=psav(is,j)
106  k(i,j)=ksav(is,j)
107  v(i,j)=vsav(is,j)
108  p(mint(83)+4+is,j)=psav(is,j)
109  v(mint(83)+4+is,j)=vsav(is,j)
110  140 CONTINUE
111  150 CONTINUE
112 C...Reset statistics on activity in event.
113  DO 160 j=351,359
114  mint(j)=0
115  vint(j)=0d0
116  160 CONTINUE
117 C...Reset extra companion reweighting factor
118  vint(140)=1d0
119 
120 C...We do not generate MI for soft process (ISUB=95), but the
121 C...initialization must be done regardless, for later purposes.
122  mint(36)=1
123 
124 C...Initialize multiple interactions.
125  CALL pyptmi(-1,ptdum1,ptdum2,ptdum3,idum)
126  IF(mint(51).NE.0) RETURN
127 
128 C...Decide whether quarks in hard scattering were valence or sea
129  pt2hd=vint(54)
130  DO 170 js=1,2
131  mint(30)=js
132  CALL pyptmi(2,pt2hd,ptdum2,ptdum3,idum)
133  IF(mint(51).NE.0) RETURN
134  170 CONTINUE
135 
136 C...Set lower cutoff for PT2 iteration and colour interference PT2 scale
137  vint(18)=0d0
138  IF(mstp(70).EQ.0) THEN
139  pt20=parp(62)**2
140  pt2min=max(pt2min,pt20,(1.1d0*alam3)**2)
141  ELSEIF(mstp(70).EQ.1) THEN
142  pt20=(parp(81)*(vint(1)/parp(89))**parp(90))**2
143  pt2min=max(pt2min,pt20,(1.1d0*alam3)**2)
144  ELSE
145  vint(18)=(parp(82)*(vint(1)/parp(89))**parp(90))**2
146  pt2min=max(pt2min,(1.1d0*alam3)**2)
147  ENDIF
148 C...Also store PT2MIN in VINT(17).
149  180 vint(17)=pt2min
150 
151 C...Set FS masses zero now.
152  vint(63)=0d0
153  vint(64)=0d0
154 
155 C...Initialize IS showers with VINT(56) as max scale.
156  pt2isr=vint(56)
157  CALL pyptis(-1,pt2isr,pt2min,pt2dum,ifail)
158  IF(mint(51).NE.0) RETURN
159 
160  RETURN
161 
162 C----------------------------------------------------------------------
163 C...MODE= 1: Evolve event from PTMAX to PTMIN.
164  ELSEIF (mode.EQ.1) THEN
165 
166 C...Skip if no phase space.
167  190 IF (pt2max.LE.pt2min) goto 330
168 
169 C...Starting pT2 max scale (to be udpated successively).
170  pt2cmx=pt2max
171 
172 C...Evolve two sides of the event to find which branches at highest pT.
173  200 jsmx=-1
174  mimx=0
175  pt2mx=0d0
176 
177 C...Loop over current shower initiators.
178  IF (mstp(61).GE.1) THEN
179  DO 230 mi=1,mint(31)
180  IF (mi.GE.2.AND.mstp(84).LE.0) goto 230
181  isub=96
182  IF (mi.EQ.1) isub=isubhd
183  mint(1)=isub
184  mint(36)=mi
185 C...Set up shat, initiator x values, and x remaining in BR.
186  vint(44)=shat(mi)
187  vint(141)=xmi(1,mi)
188  vint(142)=xmi(2,mi)
189  vint(143)=1d0
190  vint(144)=1d0
191  DO 210 ji=1,mint(31)
192  IF (ji.EQ.mint(36)) goto 210
193  vint(143)=vint(143)-xmi(1,ji)
194  vint(144)=vint(144)-xmi(2,ji)
195  210 CONTINUE
196 C...Loop over sides.
197 C...Generate trial branchings for this interaction. The hardest
198 C...branching so far is automatically updated if necessary in /PYISMX/.
199  DO 220 js=1,2
200  mint(30)=js
201  CALL pyptis(0,pt2cmx,pt2min,pt2new,ifail)
202  IF (mint(51).NE.0) RETURN
203  220 CONTINUE
204  230 CONTINUE
205  ENDIF
206 
207 C...Generate trial additional interaction.
208  mint(36)=mint(31)+1
209  240 IF (mod(mstp(81),10).GE.1) THEN
210  mint(1)=96
211 C...Set up X remaining in BR.
212  vint(143)=1d0
213  vint(144)=1d0
214  DO 250 ji=1,mint(31)
215  vint(143)=vint(143)-xmi(1,ji)
216  vint(144)=vint(144)-xmi(2,ji)
217  250 CONTINUE
218 C...Generate trial interaction
219  260 CALL pyptmi(0,pt2cmx,pt2min,pt2new,ifail)
220  IF (mint(51).EQ.1) RETURN
221  ENDIF
222 
223 C...And the winner is:
224  IF (pt2mx.LT.pt2min) THEN
225  goto 330
226  ELSEIF (jsmx.EQ.0) THEN
227 C...Accept additional interaction (may still fail).
228  CALL pyptmi(1,pt2new,pt2min,pt2dum,ifail)
229  IF(mint(51).NE.0) RETURN
230  IF (ifail.EQ.0) THEN
231  shat(mint(36))=vint(44)
232 C...Decide on flavours (valence/sea/companion).
233  DO 270 js=1,2
234  mint(30)=js
235  CALL pyptmi(2,pt2new,pt2min,pt2dum,ifail)
236  IF(mint(51).NE.0) RETURN
237  270 CONTINUE
238  ENDIF
239  ELSEIF (jsmx.EQ.1.OR.jsmx.EQ.2) THEN
240 C...Reconstruct kinematics of acceptable ISR branching.
241 C...Set up shat, initiator x values, and x remaining in BR.
242  mint(30)=jsmx
243  mint(36)=mimx
244  vint(44)=shat(mint(36))
245  vint(141)=xmi(1,mint(36))
246  vint(142)=xmi(2,mint(36))
247  vint(143)=1d0
248  vint(144)=1d0
249  DO 280 ji=1,mint(31)
250  IF (ji.EQ.mint(36)) goto 280
251  vint(143)=vint(143)-xmi(1,ji)
252  vint(144)=vint(144)-xmi(2,ji)
253  280 CONTINUE
254  pt2new=pt2mx
255  CALL pyptis(1,pt2new,pt2dm1,pt2dm2,ifail)
256  IF (mint(51).EQ.1) RETURN
257  ELSEIF (jsmx.EQ.3.OR.jsmx.EQ.4) THEN
258 C...Bookeep joining. Cannot (yet) be constructed kinematically.
259  mint(354)=mint(354)+1
260  vint(354)=vint(354)+sqrt(pt2mx)
261  IF (mint(354).EQ.1) vint(359)=sqrt(pt2mx)
262  mjoind(jsmx-2,mjn1mx)=mjn2mx
263  mjoind(jsmx-2,mjn2mx)=mjn1mx
264  ENDIF
265 
266 C...Update PT2 iteration scale.
267  pt2cmx=pt2mx
268 
269 C...Loop back to continue evolution.
270  IF(n.GT.mstu(4)-mstu(32)-10) THEN
271  CALL pyerrm(11,'(PYEVOL:) no more memory left in PYJETS')
272  ELSE
273  IF (jsmx.GE.0.AND.pt2cmx.GE.pt2min) goto 200
274  ENDIF
275 
276 C----------------------------------------------------------------------
277 C...MODE= 2: (Re-)store user information on hardest interaction etc.
278  ELSEIF (mode.EQ.2) THEN
279 
280 C...Revert to "ordinary" meanings of some parameters.
281  290 DO 310 js=1,2
282  mint(12+js)=k(imi(js,1,1),2)
283  vint(140+js)=xmi(js,1)
284  IF(mint(18+js).EQ.1) vint(140+js)=vint(154+js)*xmi(js,1)
285  vint(142+js)=1d0
286  DO 300 mi=1,mint(31)
287  vint(142+js)=vint(142+js)-xmi(js,mi)
288  300 CONTINUE
289  310 CONTINUE
290 
291 C...Restore saved quantities for hardest interaction.
292  mint(1)=isubhd
293  mint(15)=m15sv
294  mint(16)=m16sv
295  mint(21)=m21sv
296  mint(22)=m22sv
297  DO 320 j=11,80
298  vint(j)=vintsv(j)
299  320 CONTINUE
300 
301  ENDIF
302 
303  330 RETURN
304  END