]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/pythia/pykmap.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New version of MUON from M.Bondila & A.Morsch
[u/mrichter/AliRoot.git] / PYTHIA / pythia / pykmap.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE PYKMAP(IVAR,MVAR,VVAR)
5
6C...Maps a uniform distribution into a distribution of a kinematical
7C...variable according to one of the possibilities allowed. It is
8C...assumed that kinematical limits have been set by a PYKLIM call.
9 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
10 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
11 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
12 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
13 COMMON/PYINT1/MINT(400),VINT(400)
14 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
15 SAVE /LUDAT1/,/LUDAT2/
16 SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/
17
18C...Convert VVAR to tau variable.
19 ISUB=MINT(1)
20 ISTSB=ISET(ISUB)
21 IF(IVAR.EQ.1) THEN
22 TAUMIN=VINT(11)
23 TAUMAX=VINT(31)
24 IF(MVAR.EQ.3.OR.MVAR.EQ.4) THEN
25 TAURE=VINT(73)
26 GAMRE=VINT(74)
27 ELSEIF(MVAR.EQ.5.OR.MVAR.EQ.6) THEN
28 TAURE=VINT(75)
29 GAMRE=VINT(76)
30 ENDIF
31 IF(MINT(47).EQ.1.AND.(ISTSB.EQ.1.OR.ISTSB.EQ.2.OR.ISTSB.EQ.6))
32 & THEN
33 TAU=1.
34 ELSEIF(MVAR.EQ.1) THEN
35 TAU=TAUMIN*(TAUMAX/TAUMIN)**VVAR
36 ELSEIF(MVAR.EQ.2) THEN
37 TAU=TAUMAX*TAUMIN/(TAUMIN+(TAUMAX-TAUMIN)*VVAR)
38 ELSEIF(MVAR.EQ.3.OR.MVAR.EQ.5) THEN
39 RATGEN=(TAURE+TAUMAX)/(TAURE+TAUMIN)*TAUMIN/TAUMAX
40 TAU=TAURE*TAUMIN/((TAURE+TAUMIN)*RATGEN**VVAR-TAUMIN)
41 ELSEIF(MVAR.EQ.4.OR.MVAR.EQ.6) THEN
42 AUPP=ATAN((TAUMAX-TAURE)/GAMRE)
43 ALOW=ATAN((TAUMIN-TAURE)/GAMRE)
44 TAU=TAURE+GAMRE*TAN(ALOW+(AUPP-ALOW)*VVAR)
45 ELSE
46 AUPP=LOG(MAX(2E-6,1.-TAUMAX))
47 ALOW=LOG(MAX(2E-6,1.-TAUMIN))
48 TAU=1.-EXP(AUPP+VVAR*(ALOW-AUPP))
49 ENDIF
50 VINT(21)=MIN(TAUMAX,MAX(TAUMIN,TAU))
51
52C...Convert VVAR to y* variable.
53 ELSEIF(IVAR.EQ.2) THEN
54 YSTMIN=VINT(12)
55 YSTMAX=VINT(32)
56 TAUE=VINT(21)
57 IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=VINT(26)
58 IF(MINT(47).EQ.1) THEN
59 YST=0.
60 ELSEIF(MINT(47).EQ.2) THEN
61 YST=-0.5*LOG(TAUE)
62 ELSEIF(MINT(47).EQ.3) THEN
63 YST=0.5*LOG(TAUE)
64 ELSEIF(MVAR.EQ.1) THEN
65 YST=YSTMIN+(YSTMAX-YSTMIN)*SQRT(VVAR)
66 ELSEIF(MVAR.EQ.2) THEN
67 YST=YSTMAX-(YSTMAX-YSTMIN)*SQRT(1.-VVAR)
68 ELSEIF(MVAR.EQ.3) THEN
69 AUPP=ATAN(EXP(YSTMAX))
70 ALOW=ATAN(EXP(YSTMIN))
71 YST=LOG(TAN(ALOW+(AUPP-ALOW)*VVAR))
72 ELSEIF(MVAR.EQ.4) THEN
73 YST0=-0.5*LOG(TAUE)
74 AUPP=LOG(MAX(1E-6,EXP(YST0-YSTMIN)-1.))
75 ALOW=LOG(MAX(1E-6,EXP(YST0-YSTMAX)-1.))
76 YST=YST0-LOG(1.+EXP(ALOW+VVAR*(AUPP-ALOW)))
77 ELSE
78 YST0=-0.5*LOG(TAUE)
79 AUPP=LOG(MAX(1E-6,EXP(YST0+YSTMIN)-1.))
80 ALOW=LOG(MAX(1E-6,EXP(YST0+YSTMAX)-1.))
81 YST=LOG(1.+EXP(AUPP+VVAR*(ALOW-AUPP)))-YST0
82 ENDIF
83 VINT(22)=MIN(YSTMAX,MAX(YSTMIN,YST))
84
85C...Convert VVAR to cos(theta-hat) variable.
86 ELSEIF(IVAR.EQ.3) THEN
87 RM34=MAX(1E-20,2.*VINT(63)*VINT(64)/(VINT(21)*VINT(2))**2)
88 RSQM=1.+RM34
89 IF(2.*VINT(71)**2/(VINT(21)*VINT(2)).LT.0.0001) RM34=MAX(RM34,
90 & 2.*VINT(71)**2/(VINT(21)*VINT(2)))
91 CTNMIN=VINT(13)
92 CTNMAX=VINT(33)
93 CTPMIN=VINT(14)
94 CTPMAX=VINT(34)
95 IF(MVAR.EQ.1) THEN
96 ANEG=CTNMAX-CTNMIN
97 APOS=CTPMAX-CTPMIN
98 IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
99 VCTN=VVAR*(ANEG+APOS)/ANEG
100 CTH=CTNMIN+(CTNMAX-CTNMIN)*VCTN
101 ELSE
102 VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
103 CTH=CTPMIN+(CTPMAX-CTPMIN)*VCTP
104 ENDIF
105 ELSEIF(MVAR.EQ.2) THEN
106 RMNMIN=MAX(RM34,RSQM-CTNMIN)
107 RMNMAX=MAX(RM34,RSQM-CTNMAX)
108 RMPMIN=MAX(RM34,RSQM-CTPMIN)
109 RMPMAX=MAX(RM34,RSQM-CTPMAX)
110 ANEG=LOG(RMNMIN/RMNMAX)
111 APOS=LOG(RMPMIN/RMPMAX)
112 IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
113 VCTN=VVAR*(ANEG+APOS)/ANEG
114 CTH=RSQM-RMNMIN*(RMNMAX/RMNMIN)**VCTN
115 ELSE
116 VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
117 CTH=RSQM-RMPMIN*(RMPMAX/RMPMIN)**VCTP
118 ENDIF
119 ELSEIF(MVAR.EQ.3) THEN
120 RMNMIN=MAX(RM34,RSQM+CTNMIN)
121 RMNMAX=MAX(RM34,RSQM+CTNMAX)
122 RMPMIN=MAX(RM34,RSQM+CTPMIN)
123 RMPMAX=MAX(RM34,RSQM+CTPMAX)
124 ANEG=LOG(RMNMAX/RMNMIN)
125 APOS=LOG(RMPMAX/RMPMIN)
126 IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
127 VCTN=VVAR*(ANEG+APOS)/ANEG
128 CTH=RMNMIN*(RMNMAX/RMNMIN)**VCTN-RSQM
129 ELSE
130 VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
131 CTH=RMPMIN*(RMPMAX/RMPMIN)**VCTP-RSQM
132 ENDIF
133 ELSEIF(MVAR.EQ.4) THEN
134 RMNMIN=MAX(RM34,RSQM-CTNMIN)
135 RMNMAX=MAX(RM34,RSQM-CTNMAX)
136 RMPMIN=MAX(RM34,RSQM-CTPMIN)
137 RMPMAX=MAX(RM34,RSQM-CTPMAX)
138 ANEG=1./RMNMAX-1./RMNMIN
139 APOS=1./RMPMAX-1./RMPMIN
140 IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
141 VCTN=VVAR*(ANEG+APOS)/ANEG
142 CTH=RSQM-1./(1./RMNMIN+ANEG*VCTN)
143 ELSE
144 VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
145 CTH=RSQM-1./(1./RMPMIN+APOS*VCTP)
146 ENDIF
147 ELSEIF(MVAR.EQ.5) THEN
148 RMNMIN=MAX(RM34,RSQM+CTNMIN)
149 RMNMAX=MAX(RM34,RSQM+CTNMAX)
150 RMPMIN=MAX(RM34,RSQM+CTPMIN)
151 RMPMAX=MAX(RM34,RSQM+CTPMAX)
152 ANEG=1./RMNMIN-1./RMNMAX
153 APOS=1./RMPMIN-1./RMPMAX
154 IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
155 VCTN=VVAR*(ANEG+APOS)/ANEG
156 CTH=1./(1./RMNMIN-ANEG*VCTN)-RSQM
157 ELSE
158 VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
159 CTH=1./(1./RMPMIN-APOS*VCTP)-RSQM
160 ENDIF
161 ENDIF
162 IF(CTH.LT.0.) CTH=MIN(CTNMAX,MAX(CTNMIN,CTH))
163 IF(CTH.GT.0.) CTH=MIN(CTPMAX,MAX(CTPMIN,CTH))
164 VINT(23)=CTH
165
166C...Convert VVAR to tau' variable.
167 ELSEIF(IVAR.EQ.4) THEN
168 TAU=VINT(21)
169 TAUPMN=VINT(16)
170 TAUPMX=VINT(36)
171 IF(MINT(47).EQ.1) THEN
172 TAUP=1.
173 ELSEIF(MVAR.EQ.1) THEN
174 TAUP=TAUPMN*(TAUPMX/TAUPMN)**VVAR
175 ELSEIF(MVAR.EQ.2) THEN
176 AUPP=(1.-TAU/TAUPMX)**4
177 ALOW=(1.-TAU/TAUPMN)**4
178 TAUP=TAU/MAX(1E-7,1.-(ALOW+(AUPP-ALOW)*VVAR)**0.25)
179 ELSE
180 AUPP=LOG(MAX(2E-6,1.-TAUPMX))
181 ALOW=LOG(MAX(2E-6,1.-TAUPMN))
182 TAUP=1.-EXP(AUPP+VVAR*(ALOW-AUPP))
183 ENDIF
184 VINT(26)=MIN(TAUPMX,MAX(TAUPMN,TAUP))
185
186C...Selection of extra variables needed in 2 -> 3 process:
187C...pT1, pT2, phi1, phi2, y3 for three outgoing particles.
188C...Since no options are available, the functions of PYKLIM
189C...and PYKMAP are joint for these choices.
190 ELSEIF(IVAR.EQ.5) THEN
191
192C...Read out total energy and particle masses.
193 MINT(51)=0
194 MPTPK=1
195 IF(ISUB.EQ.123.OR.ISUB.EQ.124.OR.ISUB.EQ.173.OR.ISUB.EQ.174
196 & .OR.ISUB.EQ.178.OR.ISUB.EQ.179) MPTPK=2
197 SHP=VINT(26)*VINT(2)
198 SHPR=SQRT(SHP)
199 PM1=VINT(201)
200 PM2=VINT(206)
201 PM3=SQRT(VINT(21))*VINT(1)
202 IF(PM1+PM2+PM3.GT.0.9999*SHPR) THEN
203 MINT(51)=1
204 RETURN
205 ENDIF
206 PMRS1=VINT(204)**2
207 PMRS2=VINT(209)**2
208
209C...Specify coefficients of pT choice; upper and lower limits.
210 IF(MPTPK.EQ.1) THEN
211 HWT1=0.4
212 HWT2=0.4
213 ELSE
214 HWT1=0.05
215 HWT2=0.05
216 ENDIF
217 HWT3=1.-HWT1-HWT2
218 PTSMX1=((SHP-PM1**2-(PM2+PM3)**2)**2-(2.*PM1*(PM2+PM3))**2)/
219 & (4.*SHP)
220 IF(CKIN(52).GT.0.) PTSMX1=MIN(PTSMX1,CKIN(52)**2)
221 PTSMN1=CKIN(51)**2
222 PTSMX2=((SHP-PM2**2-(PM1+PM3)**2)**2-(2.*PM2*(PM1+PM3))**2)/
223 & (4.*SHP)
224 IF(CKIN(54).GT.0.) PTSMX2=MIN(PTSMX2,CKIN(54)**2)
225 PTSMN2=CKIN(53)**2
226
227C...Select transverse momenta according to
228C...dp_T^2 * (a + b/(M^2 + p_T^2) + c/(M^2 + p_T^2)^2).
229 HMX=PMRS1+PTSMX1
230 HMN=PMRS1+PTSMN1
231 IF(HMX.LT.1.0001*HMN) THEN
232 MINT(51)=1
233 RETURN
234 ENDIF
235 HDE=PTSMX1-PTSMN1
236 RPT=RLU(0)
237 IF(RPT.LT.HWT1) THEN
238 PTS1=PTSMN1+RLU(0)*HDE
239 ELSEIF(RPT.LT.HWT1+HWT2) THEN
240 PTS1=MAX(PTSMN1,HMN*(HMX/HMN)**RLU(0)-PMRS1)
241 ELSE
242 PTS1=MAX(PTSMN1,HMN*HMX/(HMN+RLU(0)*HDE)-PMRS1)
243 ENDIF
244 WTPTS1=HDE/(HWT1+HWT2*HDE/(LOG(HMX/HMN)*(PMRS1+PTS1))+
245 & HWT3*HMN*HMX/(PMRS1+PTS1)**2)
246 HMX=PMRS2+PTSMX2
247 HMN=PMRS2+PTSMN2
248 IF(HMX.LT.1.0001*HMN) THEN
249 MINT(51)=1
250 RETURN
251 ENDIF
252 HDE=PTSMX2-PTSMN2
253 RPT=RLU(0)
254 IF(RPT.LT.HWT1) THEN
255 PTS2=PTSMN2+RLU(0)*HDE
256 ELSEIF(RPT.LT.HWT1+HWT2) THEN
257 PTS2=MAX(PTSMN2,HMN*(HMX/HMN)**RLU(0)-PMRS2)
258 ELSE
259 PTS2=MAX(PTSMN2,HMN*HMX/(HMN+RLU(0)*HDE)-PMRS2)
260 ENDIF
261 WTPTS2=HDE/(HWT1+HWT2*HDE/(LOG(HMX/HMN)*(PMRS2+PTS2))+
262 & HWT3*HMN*HMX/(PMRS2+PTS2)**2)
263
264C...Select azimuthal angles and check pT choice.
265 PHI1=PARU(2)*RLU(0)
266 PHI2=PARU(2)*RLU(0)
267 PHIR=PHI2-PHI1
268 PTS3=MAX(0.,PTS1+PTS2+2.*SQRT(PTS1*PTS2)*COS(PHIR))
269 IF(PTS3.LT.CKIN(55)**2.OR.(CKIN(56).GT.0..AND.PTS3.GT.
270 & CKIN(56)**2)) THEN
271 MINT(51)=1
272 RETURN
273 ENDIF
274
275C...Calculate transverse masses and check phase space not closed.
276 PMS1=PM1**2+PTS1
277 PMS2=PM2**2+PTS2
278 PMS3=PM3**2+PTS3
279 PMT1=SQRT(PMS1)
280 PMT2=SQRT(PMS2)
281 PMT3=SQRT(PMS3)
282 PM12=(PMT1+PMT2)**2
283 IF(PMT1+PMT2+PMT3.GT.0.9999*SHPR) THEN
284 MINT(51)=1
285 RETURN
286 ENDIF
287
288C...Select rapidity for particle 3 and check phase space not closed.
289 Y3MAX=LOG((SHP+PMS3-PM12+SQRT(MAX(0.,(SHP-PMS3-PM12)**2-
290 & 4.*PMS3*PM12)))/(2.*SHPR*PMT3))
291 IF(Y3MAX.LT.1E-6) THEN
292 MINT(51)=1
293 RETURN
294 ENDIF
295 Y3=(2.*RLU(0)-1.)*0.999999*Y3MAX
296 PZ3=PMT3*SINH(Y3)
297 PE3=PMT3*COSH(Y3)
298
299C...Find momentum transfers in two mirror solutions (in 1-2 frame).
300 PZ12=-PZ3
301 PE12=SHPR-PE3
302 PMS12=PE12**2-PZ12**2
303 SQL12=SQRT(MAX(0.,(PMS12-PMS1-PMS2)**2-4.*PMS1*PMS2))
304 IF(SQL12.LT.1E-6*SHP) THEN
305 MINT(51)=1
306 RETURN
307 ENDIF
308 PMM1=PMS12+PMS1-PMS2
309 PMM2=PMS12+PMS2-PMS1
310 TFAC=-SHPR/(2.*PMS12)
311 T1P=TFAC*(PE12-PZ12)*(PMM1-SQL12)
312 T1N=TFAC*(PE12-PZ12)*(PMM1+SQL12)
313 T2P=TFAC*(PE12+PZ12)*(PMM2-SQL12)
314 T2N=TFAC*(PE12+PZ12)*(PMM2+SQL12)
315
316C...Construct relative mirror weights and make choice.
317 IF(MPTPK.EQ.1) THEN
318 WTPU=1.
319 WTNU=1.
320 ELSE
321 WTPU=1./((T1P-PMRS1)*(T2P-PMRS2))**2
322 WTNU=1./((T1N-PMRS1)*(T2N-PMRS2))**2
323 ENDIF
324 WTP=WTPU/(WTPU+WTNU)
325 WTN=WTNU/(WTPU+WTNU)
326 EPS=1.
327 IF(WTN.GT.RLU(0)) EPS=-1.
328
329C...Store result of variable choice and associated weights.
330 VINT(202)=PTS1
331 VINT(207)=PTS2
332 VINT(203)=PHI1
333 VINT(208)=PHI2
334 VINT(205)=WTPTS1
335 VINT(210)=WTPTS2
336 VINT(211)=Y3
337 VINT(212)=Y3MAX
338 VINT(213)=EPS
339 IF(EPS.GT.0.) THEN
340 VINT(214)=1./WTP
341 VINT(215)=T1P
342 VINT(216)=T2P
343 ELSE
344 VINT(214)=1./WTN
345 VINT(215)=T1N
346 VINT(216)=T2N
347 ENDIF
348 VINT(217)=-0.5*TFAC*(PE12-PZ12)*(PMM2+EPS*SQL12)
349 VINT(218)=-0.5*TFAC*(PE12+PZ12)*(PMM1+EPS*SQL12)
350 VINT(219)=0.5*(PMS12-PTS3)
351 VINT(220)=SQL12
352 ENDIF
353
354 RETURN
355 END