Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyinit.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyinit.f
1 
2 C*********************************************************************
3 
4 C...PYINIT
5 C...Initializes the generation procedure; finds maxima of the
6 C...differential cross-sections to be used for weighting.
7 
8  SUBROUTINE pyinit(FRAME,BEAM,TARGET,WIN)
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/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
21  common/pypars/mstp(200),parp(200),msti(200),pari(200)
22  common/pyint1/mint(400),vint(400)
23  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
24  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
25  SAVE /pydat1/,/pydat2/,/pydat3/,/pydat4/,/pysubs/,/pypars/,
26  &/pyint1/,/pyint2/,/pyint5/
27 C...Local arrays and character variables.
28  dimension alamin(20),nfin(20)
29  CHARACTER*(*) frame,beam,target
30  CHARACTER chfram*12,chbeam*12,chtarg*12,chlh(2)*6
31 
32 C...Interface to PDFLIB.
33  common/w50511/nptype,ngroup,nset,mode,nfl,lo,tmas
34  common/w50512/qcdl4,qcdl5
35  SAVE /w50511/,/w50512/
36  DOUBLE PRECISION value(20),tmas,qcdl4,qcdl5
37  CHARACTER*20 parm(20)
38  DATA value/20*0d0/,parm/20*' '/
39 
40 C...Data:Lambda and n_f values for parton distributions..
41  DATA alamin/0.177d0,0.239d0,0.247d0,0.2322d0,0.248d0,0.248d0,
42  &0.192d0,0.326d0,2*0.2d0,0.2d0,0.2d0,0.29d0,0.2d0,0.4d0,5*0.2d0/,
43  &nfin/20*4/
44  DATA chlh/'lepton','hadron'/
45 
46 C...Check that BLOCK DATA PYDATA has been loaded.
47  CALL pyckbd
48 
49 C...Reset MINT and VINT arrays. Write headers.
50  msti(53)=0
51  DO 100 j=1,400
52  mint(j)=0
53  vint(j)=0d0
54  100 CONTINUE
55  IF(mstu(12).NE.12345) CALL pylist(0)
56  IF(mstp(122).GE.1) WRITE(mstu(11),5100)
57 
58 C...Reset error counters.
59  mstu(23)=0
60  mstu(27)=0
61  mstu(30)=0
62 
63 C...Reset processes that should not be on.
64  msub(96)=0
65  msub(97)=0
66 
67 C...Select global FSR/ISR/UE parameter set = 'tune'
68 C...See routine PYTUNE for details
69  IF (mstp(5).NE.0) THEN
70  mstp5=mstp(5)
71  CALL pytune(mstp5)
72  ENDIF
73 
74 C...Call user process initialization routine.
75  IF(frame(1:1).EQ.'u'.OR.frame(1:1).EQ.'U') THEN
76  msel=0
77  CALL upinit
78  msel=0
79  ENDIF
80 
81 C...Maximum 4 generations; set maximum number of allowed flavours.
82  mstp(1)=min(4,mstp(1))
83  mstu(114)=min(mstu(114),2*mstp(1))
84  mstp(58)=min(mstp(58),2*mstp(1))
85 
86 C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton.
87  DO 120 i=-20,20
88  vint(180+i)=0d0
89  ia=iabs(i)
90  IF(ia.GE.1.AND.ia.LE.2*mstp(1)) THEN
91  DO 110 j=1,mstp(1)
92  ib=2*j-1+mod(ia,2)
93  IF(ib.GE.6.AND.mstp(9).EQ.0) goto 110
94  ipm=(5-isign(1,i))/2
95  idc=j+mdcy(ia,2)+2
96  IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.ipm) vint(180+i)=
97  & vint(180+i)+vckm((ia+1)/2,(ib+1)/2)
98  110 CONTINUE
99  ELSEIF(ia.GE.11.AND.ia.LE.10+2*mstp(1)) THEN
100  vint(180+i)=1d0
101  ENDIF
102  120 CONTINUE
103 
104 C...Initialize parton distributions: PDFLIB.
105  IF(mstp(52).EQ.2) THEN
106  parm(1)='NPTYPE'
107  value(1)=1
108  parm(2)='NGROUP'
109  value(2)=mstp(51)/1000
110  parm(3)='NSET'
111  value(3)=mod(mstp(51),1000)
112  parm(4)='TMAS'
113  value(4)=pmas(6,1)
114  CALL pdfset(parm,value)
115  mint(93)=1000000+mstp(51)
116  ENDIF
117 
118 C...Choose Lambda value to use in alpha-strong.
119  mstu(111)=mstp(2)
120  IF(mstp(3).GE.2) THEN
121  alam=0.2d0
122  nf=4
123  IF(mstp(52).EQ.1.AND.mstp(51).GE.1.AND.mstp(51).LE.20) THEN
124  alam=alamin(mstp(51))
125  nf=nfin(mstp(51))
126  ELSEIF(mstp(52).EQ.2.AND.nfl.EQ.5) THEN
127  alam=qcdl5
128  nf=5
129  ELSEIF(mstp(52).EQ.2) THEN
130  alam=qcdl4
131  nf=4
132  ENDIF
133  parp(1)=alam
134  parp(61)=alam
135  parp(72)=alam
136  paru(112)=alam
137  mstu(112)=nf
138  IF(mstp(3).EQ.3) parj(81)=alam
139  ENDIF
140 
141 C...Initialize the SUSY generation: couplings, masses,
142 C...decay modes, branching ratios, and so on.
143  CALL pymsin
144 C...Initialize widths and partial widths for resonances.
145  CALL pyinre
146 C...Set Z0 mass and width for e+e- routines.
147  parj(123)=pmas(23,1)
148  parj(124)=pmas(23,2)
149 
150 C...Identify beam and target particles and frame of process.
151  chfram=frame//' '
152  chbeam=beam//' '
153  chtarg=TARGET//' '
154  CALL pyinbm(chfram,chbeam,chtarg,win)
155  IF(mint(65).EQ.1) goto 170
156 
157 C...For gamma-p or gamma-gamma allow many (3 or 6) alternatives.
158 C...For e-gamma allow 2 alternatives.
159  mint(121)=1
160  IF(mstp(14).EQ.10.AND.(msel.EQ.1.OR.msel.EQ.2)) THEN
161  IF((mint(11).EQ.22.OR.mint(12).EQ.22).AND.
162  & (iabs(mint(11)).GT.100.OR.iabs(mint(12)).GT.100)) mint(121)=3
163  IF(mint(11).EQ.22.AND.mint(12).EQ.22) mint(121)=6
164  IF((mint(11).EQ.22.OR.mint(12).EQ.22).AND.
165  & (iabs(mint(11)).EQ.11.OR.iabs(mint(12)).EQ.11)) mint(121)=2
166  ELSEIF(mstp(14).EQ.20.AND.(msel.EQ.1.OR.msel.EQ.2)) THEN
167  IF((mint(11).EQ.22.OR.mint(12).EQ.22).AND.
168  & (iabs(mint(11)).GT.100.OR.iabs(mint(12)).GT.100)) mint(121)=3
169  IF(mint(11).EQ.22.AND.mint(12).EQ.22) mint(121)=9
170  ELSEIF(mstp(14).EQ.25.AND.(msel.EQ.1.OR.msel.EQ.2)) THEN
171  IF((mint(11).EQ.22.OR.mint(12).EQ.22).AND.
172  & (iabs(mint(11)).GT.100.OR.iabs(mint(12)).GT.100)) mint(121)=2
173  IF(mint(11).EQ.22.AND.mint(12).EQ.22) mint(121)=4
174  ELSEIF(mstp(14).EQ.30.AND.(msel.EQ.1.OR.msel.EQ.2)) THEN
175  IF((mint(11).EQ.22.OR.mint(12).EQ.22).AND.
176  & (iabs(mint(11)).GT.100.OR.iabs(mint(12)).GT.100)) mint(121)=4
177  IF(mint(11).EQ.22.AND.mint(12).EQ.22) mint(121)=13
178  ENDIF
179  mint(123)=mstp(14)
180  IF((mstp(14).EQ.10.OR.mstp(14).EQ.20.OR.mstp(14).EQ.25.OR.
181  &mstp(14).EQ.30).AND.msel.NE.1.AND.msel.NE.2) mint(123)=0
182  IF(mstp(14).GE.11.AND.mstp(14).LE.19) THEN
183  IF(mstp(14).EQ.11) mint(123)=0
184  IF(mstp(14).EQ.12.OR.mstp(14).EQ.14) mint(123)=5
185  IF(mstp(14).EQ.13.OR.mstp(14).EQ.17) mint(123)=6
186  IF(mstp(14).EQ.15) mint(123)=2
187  IF(mstp(14).EQ.16.OR.mstp(14).EQ.18) mint(123)=7
188  IF(mstp(14).EQ.19) mint(123)=3
189  ELSEIF(mstp(14).GE.21.AND.mstp(14).LE.24) THEN
190  IF(mstp(14).EQ.21) mint(123)=0
191  IF(mstp(14).EQ.22.OR.mstp(14).EQ.23) mint(123)=4
192  IF(mstp(14).EQ.24) mint(123)=1
193  ELSEIF(mstp(14).GE.26.AND.mstp(14).LE.29) THEN
194  IF(mstp(14).EQ.26.OR.mstp(14).EQ.28) mint(123)=8
195  IF(mstp(14).EQ.27.OR.mstp(14).EQ.29) mint(123)=9
196  ENDIF
197 
198 C...Set up kinematics of process.
199  CALL pyinki(0)
200 
201 C...Set up kinematics for photons inside leptons.
202  IF(mint(141).NE.0.OR.mint(142).NE.0) CALL pygaga(1,wtgaga)
203 
204 C...Precalculate flavour selection weights.
205  CALL pykfin
206 
207 C...Loop over gamma-p or gamma-gamma alternatives.
208  ckin3=ckin(3)
209  msav48=0
210  DO 160 iga=1,mint(121)
211  ckin(3)=ckin3
212  mint(122)=iga
213 
214 C...Select partonic subprocesses to be included in the simulation.
215  CALL pyinpr
216  mint(101)=1
217  mint(102)=1
218  mint(103)=mint(11)
219  mint(104)=mint(12)
220 
221 C...Count number of subprocesses on.
222  mint(48)=0
223  DO 130 isub=1,500
224  IF(mint(50).EQ.0.AND.isub.GE.91.AND.isub.LE.96.AND.
225  & msub(isub).EQ.1.AND.mint(121).GT.1) THEN
226  msub(isub)=0
227  ELSEIF(mint(50).EQ.0.AND.isub.GE.91.AND.isub.LE.96.AND.
228  & msub(isub).EQ.1) THEN
229  WRITE(mstu(11),5200) isub,chlh(mint(41)),chlh(mint(42))
230  CALL pystop(1)
231  ELSEIF(msub(isub).EQ.1.AND.iset(isub).EQ.-1) THEN
232  WRITE(mstu(11),5300) isub
233  CALL pystop(1)
234  ELSEIF(msub(isub).EQ.1.AND.iset(isub).LE.-2) THEN
235  WRITE(mstu(11),5400) isub
236  CALL pystop(1)
237  ELSEIF(msub(isub).EQ.1) THEN
238  mint(48)=mint(48)+1
239  ENDIF
240  130 CONTINUE
241 
242 C...Stop or raise warning flag if no subprocesses on.
243  IF(mint(121).EQ.1.AND.mint(48).EQ.0) THEN
244  IF(mstp(127).NE.1) THEN
245  WRITE(mstu(11),5500)
246  CALL pystop(1)
247  ELSE
248  WRITE(mstu(11),5700)
249  msti(53)=1
250  ENDIF
251  ENDIF
252  mint(49)=mint(48)-msub(91)-msub(92)-msub(93)-msub(94)
253  msav48=msav48+mint(48)
254 
255 C...Reset variables for cross-section calculation.
256  DO 150 i=0,500
257  DO 140 j=1,3
258  ngen(i,j)=0
259  xsec(i,j)=0d0
260  140 CONTINUE
261  150 CONTINUE
262 
263 C...Find parametrized total cross-sections.
264  CALL pyxtot
265  vint(318)=vint(317)
266 
267 C...Maxima of differential cross-sections.
268  IF(mstp(121).LE.1) CALL pymaxi
269 
270 C...Initialize possibility of pileup events.
271  IF(mint(121).GT.1) mstp(131)=0
272  IF(mstp(131).NE.0) CALL pypile(1)
273 
274 C...Initialize multiple interactions with variable impact parameter.
275  IF(mint(50).EQ.1) THEN
276  ptmn=parp(82)*(vint(1)/parp(89))**parp(90)
277  IF(mod(mstp(81),10).EQ.0.AND.(ckin(3).GT.ptmn.OR.
278  & ((msel.NE.1.AND.msel.NE.2)))) mstp(82)=min(1,mstp(82))
279  IF((mint(49).NE.0.OR.mstp(131).NE.0).AND.mstp(82).GE.2) THEN
280  mint(35)=1
281  CALL pymult(1)
282  mint(35)=3
283  CALL pymign(1)
284  ENDIF
285  ENDIF
286 
287 C...Save results for gamma-p and gamma-gamma alternatives.
288  IF(mint(121).GT.1) CALL pysave(1,iga)
289  160 CONTINUE
290 
291 C...Initialization finished.
292  IF(msav48.EQ.0) THEN
293  IF(mstp(127).NE.1) THEN
294  WRITE(mstu(11),5500)
295  CALL pystop(1)
296  ELSE
297  WRITE(mstu(11),5700)
298  msti(53)=1
299  ENDIF
300  ENDIF
301  170 IF(mstp(122).GE.1) WRITE(mstu(11),5600)
302 
303 C...Formats for initialization information.
304  5100 FORMAT('1',18('*'),1x,'PYINIT: initialization of PYTHIA ',
305  &'routines',1x,17('*'))
306  5200 FORMAT(1x,'Error: process number ',i3,' not meaningful for ',a6,
307  &'-',a6,' interactions.'/1x,'Execution stopped!')
308  5300 FORMAT(1x,'Error: requested subprocess',i4,' not implemented.'/
309  &1x,'Execution stopped!')
310  5400 FORMAT(1x,'Error: requested subprocess',i4,' not existing.'/
311  &1x,'Execution stopped!')
312  5500 FORMAT(1x,'Error: no subprocess switched on.'/
313  &1x,'Execution stopped.')
314  5600 FORMAT(/1x,22('*'),1x,'PYINIT: initialization completed',1x,
315  &22('*'))
316  5700 FORMAT(1x,'Error: no subprocess switched on.'/
317  &1x,'Execution will stop if you try to generate events.')
318 
319  RETURN
320  END