Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyinki.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyinki.f
1 
2 C*********************************************************************
3 
4 C...PYINKI
5 C...Sets up kinematics, including rotations and boosts to/from CM frame.
6 
7  SUBROUTINE pyinki(MODKI)
8 
9 C...Double precision and integer declarations.
10  IMPLICIT DOUBLE PRECISION(a-h, o-z)
11  IMPLICIT INTEGER(i-n)
12  INTEGER pyk,pychge,pycomp
13 
14 C...User process initialization commonblock.
15  INTEGER maxpup
16  parameter(maxpup=100)
17  INTEGER idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
18  DOUBLE PRECISION ebmup,xsecup,xerrup,xmaxup
19  common/heprup/idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
20  &idwtup,nprup,xsecup(maxpup),xerrup(maxpup),xmaxup(maxpup),
21  &lprup(maxpup)
22  SAVE /heprup/
23 
24 C...Commonblocks.
25  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
26  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
27  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
28  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
29  common/pypars/mstp(200),parp(200),msti(200),pari(200)
30  common/pyint1/mint(400),vint(400)
31  SAVE /pyjets/,/pydat1/,/pydat2/,/pysubs/,/pypars/,/pyint1/
32 
33 C...Set initial flavour state.
34  n=2
35  DO 100 i=1,2
36  k(i,1)=1
37  k(i,2)=mint(10+i)
38  IF(mint(140+i).NE.0) k(i,2)=mint(140+i)
39  100 CONTINUE
40 
41 C...Reset boost. Do kinematics for various cases.
42  DO 110 j=6,10
43  vint(j)=0d0
44  110 CONTINUE
45 
46 C...Set up kinematics for events defined in CM frame.
47  IF(mint(111).EQ.1) THEN
48  win=vint(290)
49  IF(modki.EQ.1) win=parp(171)*vint(290)
50  s=win**2
51  p(1,5)=vint(3)
52  p(2,5)=vint(4)
53  IF(mint(141).NE.0) p(1,5)=vint(303)
54  IF(mint(142).NE.0) p(2,5)=vint(304)
55  p(1,1)=0d0
56  p(1,2)=0d0
57  p(2,1)=0d0
58  p(2,2)=0d0
59  p(1,3)=sqrt(((s-p(1,5)**2-p(2,5)**2)**2-(2d0*p(1,5)*p(2,5))**2)/
60  & (4d0*s))
61  p(2,3)=-p(1,3)
62  p(1,4)=sqrt(p(1,3)**2+p(1,5)**2)
63  p(2,4)=sqrt(p(2,3)**2+p(2,5)**2)
64 
65 C...Set up kinematics for fixed target events.
66  ELSEIF(mint(111).EQ.2) THEN
67  win=vint(290)
68  IF(modki.EQ.1) win=parp(171)*vint(290)
69  p(1,5)=vint(3)
70  p(2,5)=vint(4)
71  IF(mint(141).NE.0) p(1,5)=vint(303)
72  IF(mint(142).NE.0) p(2,5)=vint(304)
73  p(1,1)=0d0
74  p(1,2)=0d0
75  p(2,1)=0d0
76  p(2,2)=0d0
77  p(1,3)=win
78  p(1,4)=sqrt(p(1,3)**2+p(1,5)**2)
79  p(2,3)=0d0
80  p(2,4)=p(2,5)
81  s=p(1,5)**2+p(2,5)**2+2d0*p(2,4)*p(1,4)
82  vint(10)=p(1,3)/(p(1,4)+p(2,4))
83  CALL pyrobo(0,0,0d0,0d0,0d0,0d0,-vint(10))
84 
85 C...Set up kinematics for events in user-defined frame.
86  ELSEIF(mint(111).EQ.3) THEN
87  p(1,5)=vint(3)
88  p(2,5)=vint(4)
89  IF(mint(141).NE.0) p(1,5)=vint(303)
90  IF(mint(142).NE.0) p(2,5)=vint(304)
91  p(1,4)=sqrt(p(1,1)**2+p(1,2)**2+p(1,3)**2+p(1,5)**2)
92  p(2,4)=sqrt(p(2,1)**2+p(2,2)**2+p(2,3)**2+p(2,5)**2)
93  DO 120 j=1,3
94  vint(7+j)=(p(1,j)+p(2,j))/(p(1,4)+p(2,4))
95  120 CONTINUE
96  CALL pyrobo(0,0,0d0,0d0,-vint(8),-vint(9),-vint(10))
97  vint(7)=pyangl(p(1,1),p(1,2))
98  CALL pyrobo(0,0,0d0,-vint(7),0d0,0d0,0d0)
99  vint(6)=pyangl(p(1,3),p(1,1))
100  CALL pyrobo(0,0,-vint(6),0d0,0d0,0d0,0d0)
101  s=p(1,5)**2+p(2,5)**2+2d0*(p(1,4)*p(2,4)-p(1,3)*p(2,3))
102 
103 C...Set up kinematics for events with user-defined four-vectors.
104  ELSEIF(mint(111).EQ.4) THEN
105  pms1=p(1,4)**2-p(1,1)**2-p(1,2)**2-p(1,3)**2
106  p(1,5)=sign(sqrt(abs(pms1)),pms1)
107  pms2=p(2,4)**2-p(2,1)**2-p(2,2)**2-p(2,3)**2
108  p(2,5)=sign(sqrt(abs(pms2)),pms2)
109  DO 130 j=1,3
110  vint(7+j)=(p(1,j)+p(2,j))/(p(1,4)+p(2,4))
111  130 CONTINUE
112  CALL pyrobo(0,0,0d0,0d0,-vint(8),-vint(9),-vint(10))
113  vint(7)=pyangl(p(1,1),p(1,2))
114  CALL pyrobo(0,0,0d0,-vint(7),0d0,0d0,0d0)
115  vint(6)=pyangl(p(1,3),p(1,1))
116  CALL pyrobo(0,0,-vint(6),0d0,0d0,0d0,0d0)
117  s=(p(1,4)+p(2,4))**2
118 
119 C...Set up kinematics for events with user-defined five-vectors.
120  ELSEIF(mint(111).EQ.5) THEN
121  DO 140 j=1,3
122  vint(7+j)=(p(1,j)+p(2,j))/(p(1,4)+p(2,4))
123  140 CONTINUE
124  CALL pyrobo(0,0,0d0,0d0,-vint(8),-vint(9),-vint(10))
125  vint(7)=pyangl(p(1,1),p(1,2))
126  CALL pyrobo(0,0,0d0,-vint(7),0d0,0d0,0d0)
127  vint(6)=pyangl(p(1,3),p(1,1))
128  CALL pyrobo(0,0,-vint(6),0d0,0d0,0d0,0d0)
129  s=(p(1,4)+p(2,4))**2
130 
131 C...Set up kinematics for events with external user processes.
132  ELSEIF(mint(111).GE.11) THEN
133  p(1,5)=vint(3)
134  p(2,5)=vint(4)
135  IF(mint(141).NE.0) p(1,5)=vint(303)
136  IF(mint(142).NE.0) p(2,5)=vint(304)
137  p(1,1)=0d0
138  p(1,2)=0d0
139  p(2,1)=0d0
140  p(2,2)=0d0
141  p(1,3)=sqrt(max(0d0,ebmup(1)**2-p(1,5)**2))
142  p(2,3)=-sqrt(max(0d0,ebmup(2)**2-p(2,5)**2))
143  p(1,4)=ebmup(1)
144  p(2,4)=ebmup(2)
145  vint(10)=(p(1,3)+p(2,3))/(p(1,4)+p(2,4))
146  CALL pyrobo(0,0,0d0,0d0,0d0,0d0,-vint(10))
147  s=(p(1,4)+p(2,4))**2
148  ENDIF
149 
150 C...Return or error for too low CM energy.
151  IF(modki.EQ.1.AND.s.LT.parp(2)**2) THEN
152  IF(mstp(172).LE.1) THEN
153  CALL pyerrm(23,
154  & '(PYINKI:) too low invariant mass in this event')
155  ELSE
156  msti(61)=1
157  RETURN
158  ENDIF
159  ENDIF
160 
161 C...Save information on incoming particles.
162  vint(1)=sqrt(s)
163  vint(2)=s
164  IF(mint(111).GE.4) THEN
165  IF(mint(141).EQ.0) THEN
166  vint(3)=p(1,5)
167  IF(mint(11).EQ.22.AND.p(1,5).LT.0) vint(307)=p(1,5)**2
168  ELSE
169  vint(303)=p(1,5)
170  ENDIF
171  IF(mint(142).EQ.0) THEN
172  vint(4)=p(2,5)
173  IF(mint(12).EQ.22.AND.p(2,5).LT.0) vint(308)=p(2,5)**2
174  ELSE
175  vint(304)=p(2,5)
176  ENDIF
177  ENDIF
178  vint(5)=p(1,3)
179  IF(modki.EQ.0) vint(289)=s
180  DO 150 j=1,5
181  v(1,j)=0d0
182  v(2,j)=0d0
183  vint(290+j)=p(1,j)
184  vint(295+j)=p(2,j)
185  150 CONTINUE
186 
187 C...Store pT cut-off and related constants to be used in generation.
188  IF(modki.EQ.0) vint(285)=ckin(3)
189  IF(mstp(82).LE.1) THEN
190  ptmn=parp(81)*(vint(1)/parp(89))**parp(90)
191  ELSE
192  ptmn=parp(82)*(vint(1)/parp(89))**parp(90)
193  ENDIF
194  vint(149)=4d0*ptmn**2/s
195  vint(154)=ptmn
196 
197  RETURN
198  END