C++ Interface to Tauola
VBF_distr.f
1 REAL*8 FUNCTION vbfdistr(ID1,ID2,ID3,ID4,HH1,HH2,PP,KEYIN)
2C***************************************************************************
3C* CALCULATES (matrix element)^2 for
4C I1 I2 -> I3 I4 TAU+TAU- with given tau polarizations H1 H2
5C for a given set of momenta P
6C I1,...,I4 are PDG particle codes
7C if I3 or I4=0, then summed over final jet flavors
8C H1 and H2 are tau+ and tau- helicities R: 1, L: -1, R+L: 0
9C KeyIn=0: no Higgs
10C KeyIn=1: SM Higgs
11C KeyIn=2: no Higgs
12C KeyIn=3: SM Higgs
13C KeyIn=4,5: non-standard state
14C***************************************************************************
15 IMPLICIT NONE
16
17 INTEGER I1,I2,I3,I4,ID1,ID2,ID3,ID4,H1,H2,HH1,HH2,BUF_H,KEY,KEYIN
18 real*8 p(0:3,6), pp(0:3,6),ans
19 INTEGER ICP
20 COMMON /cpstatus/ icp
21 real*8 buf(0:3)
22 INTEGER BUF_I,IGLU,I,J,K
23
24 LOGICAL INITIALIZED
25 DATA initialized/.false./
26 SAVE initialized
27 LOGICAL FLIPER,TESTUJEMY
28 SAVE keystored
29 INTEGER KEYSTORED
30 testujemy=id1.eq.-222.and.id2.eq.1.and.id3.eq.-1.and.id4.eq.1
31 icp=0
32
33 IF(.NOT.initialized) THEN
34 initialized = .true.
35 keystored=keyin
36 ENDIF
37
38 IF (keyin.NE.keystored) THEN
39 CALL vbf_reinit(keyin)
40 keystored=keyin
41 ENDIF
42
43 key=mod(keyin,2)
44 IF(keyin.GT.3) THEN
45 WRITE(*,*) 'non-standard state -- implementation not finished'
46 stop
47 ELSE IF(key.NE.0.AND.key.NE.1) THEN
48 WRITE(*,*) 'WRONG KEY'
49 stop
50 ENDIF
51 i1=id1
52 i2=id2
53 i3=id3
54 i4=id4
55 IF (testujemy) WRITE(*,*) 'idsy=',id1,id2,id3,id4
56C ------------
57C Fast track for cases where results must be zero.
58C ------------
59
60C bayron number conservation:
61C --------------------------
62
63 IF(i1+i2.EQ.42.AND.i3+i4.EQ.0) THEN
64C do nothing
65 ELSEIF(i3+i4.EQ.42.AND.i1+i2.EQ.0) THEN
66C do nothing
67 ELSEIF(i1*i2*i3*i4.LT.0) THEN
68 vbfdistr=0.0
69 RETURN
70 ENDIF
71
72 IF(mod(i1+i2+i3+i4,2).EQ.1) THEN
73 vbfdistr=0.0
74 RETURN
75 ENDIF
76
77 IF(i1.LT.0.and.i2.LT.0.AND.(i3.GT.0.OR.i4.GT.0)) THEN
78 vbfdistr=0.0
79 RETURN
80 ENDIF
81
82 IF(i3.LT.0.and.i4.LT.0.AND.(i1.GT.0.OR.i2.GT.0)) THEN
83 vbfdistr=0.0
84 RETURN
85 ENDIF
86
87C charge conservation
88C -------------------
89
90 IF(sign(mod(i1,2),i1)+sign(mod(i2,2),i2).NE.sign(mod(i3,2),i3)+sign(mod(i4,2),i4)) THEN
91 IF(i1+i2+i3+i4.LT.20) THEN
92 vbfdistr=0.0
93 RETURN
94 ENDIF
95 ENDIF
96
97C zero or two gluons
98C ------------------
99
100 iglu=0
101 IF(i1.EQ.21) iglu=iglu+1
102 IF(i2.EQ.21) iglu=iglu+1
103 IF(i3.EQ.21) iglu=iglu+1
104 IF(i4.EQ.21) iglu=iglu+1
105
106
107 IF(iglu.EQ.1.OR.iglu.GT.2) THEN
108 vbfdistr=0.0
109 RETURN
110 ENDIF
111
112C if gluons are present extra consitency check.
113C important, that this is before 5 to 1 replacement for configs with gluons.
114C --------------------------------------------------------------------------
115 IF(iglu.EQ.2.AND.i1+i2.NE.i3+i4) THEN
116
117 IF(.NOT.(i1+i2.EQ.0.OR.i3+i4.EQ.0)) THEN
118 vbfdistr=0.0
119 RETURN
120 ENDIF
121 ENDIF
122
123C =============================================
124C now we have a chance that result is non zero.
125C we copy configuration to local variables
126C =============================================
127 h1=hh1
128 h2=hh2
129 p(0:3,1) = pp(0:3,1)
130 p(0:3,2) = pp(0:3,2)
131 p(0:3,3) = pp(0:3,3)
132 p(0:3,4) = pp(0:3,4)
133 p(0:3,5) = pp(0:3,5)
134 p(0:3,6) = pp(0:3,6)
135
136
137C ===========================================
138C REARRANING order of 4-vectors and ID's
139C Parity reflection
140C ===========================================
141
142
143C treat quarks 5 as 1 if gluons present. It is possible
144C --------------------------------------
145 ! IF (IGLU.EQ.2) THEN
146 ! IF (ABS(I1).EQ.5) I1=I1/ABS(I1)
147 ! IF (ABS(I2).EQ.5) I2=I2/ABS(I2)
148 ! IF (ABS(I3).EQ.5) I3=I3/ABS(I3)
149 ! IF (ABS(I4).EQ.5) I4=I4/ABS(I4)
150 !ELSEIF((I1*I1.EQ.25.OR.I2*I2.EQ.25.OR.I3*I3.EQ.25.OR.I4*I4.EQ.25)) THEN
151 IF((i1*i1.EQ.25.OR.i2*i2.EQ.25.OR.i3*i3.EQ.25.OR.i4*i4.EQ.25)) THEN
152C in other cases processes with b-quarks are not yet installed.
153C -------------------------------------------------------------
154 vbfdistr=0.0
155 RETURN
156 ENDIF
157
158
159 if(testujemy) write(*,*) 'doszlimy do stepX',i1,i2, i3, i4
160
161C If I1 I2 particle/antiparticle (no gluon), |I1|>|I2| enforce
162C -------------------------------------------------------------
163 IF(i1*i2.LT.0.AND.i1+i2.LT.11.AND.i1**2.LT.i2**2) THEN
164 ! flip IDs
165 buf_i = i1
166 i1 = i2
167 i2 = buf_i
168
169 ! flip 4-vectors
170 buf(0:3) = p(0:3,1)
171 p(0:3,1) = p(0:3,2)
172 p(0:3,2) = buf(0:3)
173 ENDIF
174
175
176 if(testujemy) write(*,*) 'doszlimy do step2',i1,i2, i3, i4
177
178
179C C-reflection if all quarks antipartices
180C ---------------------------------------
181 IF(
182 $ (i1.LT.0.OR.i1.EQ.21).AND.
183 $ (i2.LT.0.OR.i2.EQ.21).AND.
184 $ (i3.LT.0.OR.i3.EQ.21).AND.
185 $ (i4.LT.0.OR.i4.EQ.21) ) THEN
186 IF(i1.NE.21) i1=-i1
187 IF(i2.NE.21) i2=-i2
188 IF(i3.NE.21) i3=-i3
189 IF(i4.NE.21) i4=-i4
190
191! C-parity on taus
192 buf_h=h2
193 h2=-h1
194 h1=-buf_h
195
196 buf(0:3) = p(0:3,5)
197 p(0:3,5) = p(0:3,6)
198 p(0:3,6) = buf(0:3)
199C and P-Parity
200 DO k=1,3
201 DO j=1,6
202 p(k,j)=-p(k,j)
203 enddo
204 enddo
205 icp=icp+1
206 ENDIF
207
208C for incoming particle antiparticle larger |id| must be first,
209C EXCEPTION: ID1 ID2 = (1,-2) or (-2,1) There is UDX.f but no DUX.f file
210C -----------------------------------------------------------------------
211
212
213 fliper=(i1*i2.LT.0.AND.i1+i2.LT.11)
214 IF(fliper) THEN
215 fliper=i2*i2.GT.i1*i1
216 IF(id1*id2.EQ.-2.AND.(id1.EQ.1.OR.id1.EQ.-2)) fliper=.NOT.fliper
217 ENDIF
218 IF(fliper) THEN
219 IF(i1.NE.21) i1=-i1
220 IF(i2.NE.21) i2=-i2
221 IF(i3.NE.21) i3=-i3
222 IF(i4.NE.21) i4=-i4
223
224! C-parity on taus
225 buf_h=h2
226 h2=-h1
227 h1=-buf_h
228 buf(0:3) = p(0:3,5)
229 p(0:3,5) = p(0:3,6)
230 p(0:3,6) = buf(0:3)
231
232C and P-Parity
233 DO k=1,3
234 DO j=1,6
235 p(k,j)=-p(k,j)
236 enddo
237 enddo
238
239 icp=icp+1
240 ENDIF
241
242 if(testujemy) write(*,*) 'doszlimy do step3',i1,i2,i3,i4
243
244C First incoming can not be antiparticle, second can not be lone gluon
245C ---------------------------------------------------------------------
246 IF(i1.LT.0.OR.(i2.EQ.21.AND.i1.NE.21)) THEN
247
248 ! flip IDs
249 buf_i = i1
250 i1 = i2
251 i2 = buf_i
252
253 ! flip 4-vectors
254 buf(0:3) = p(0:3,1)
255 p(0:3,1) = p(0:3,2)
256 p(0:3,2) = buf(0:3)
257 ENDIF
258
259C If both I1 I2 positive (no gluon) enforce that I1>I2
260 IF(i1.GT.0.AND.i2.GT.0.AND.i1+i2.LT.11.AND.i1.LT.i2) THEN
261 ! flip IDs
262 buf_i = i1
263 i1 = i2
264 i2 = buf_i
265
266 ! flip 4-vectors
267 buf(0:3) = p(0:3,1)
268 p(0:3,1) = p(0:3,2)
269 p(0:3,2) = buf(0:3)
270 ENDIF
271
272C ================
273C NOW FINAL STATES
274C ================
275
276C I3 never negative and I4 never alone 21
277C ---------------------------------------
278 IF(i3.LT.0.OR.(i4.EQ.21.AND.i3.NE.21)) THEN
279 ! flip IDs
280 buf_i = i3
281 i3 = i4
282 i4 = buf_i
283
284 ! flip 4-vectors
285 buf(0:3) = p(0:3,3)
286 p(0:3,3) = p(0:3,4)
287 p(0:3,4) = buf(0:3)
288 ENDIF
289
290C sole posibility <I3 even I4 odd> if both positive and non gluon
291C ---------------------------------------------------------------
292 IF(mod(i3,2).EQ.1.AND.mod(i4,2).EQ.0.AND.i3*i4.GT.0.AND.i3.NE.21) THEN
293 ! flip IDs
294 buf_i = i3
295 i3 = i4
296 i4 = buf_i
297
298 ! flip 4-vectors
299 buf(0:3) = p(0:3,3)
300 p(0:3,3) = p(0:3,4)
301 p(0:3,4) = buf(0:3)
302 ENDIF
303
304
305C if I3,I4 simultaneously odd I4 must be larger/equal I3
306C -------------------------------------------------------
307 IF(mod(i3,2).EQ.1.AND.mod(i4,2).EQ.1.AND.i3*i4.GT.0.AND.i3.GT.i4.AND.i3.NE.21) THEN
308 ! flip IDs
309 buf_i = i3
310 i3 = i4
311 i4 = buf_i
312
313 ! flip 4-vectors
314 buf(0:3) = p(0:3,3)
315 p(0:3,3) = p(0:3,4)
316 p(0:3,4) = buf(0:3)
317 ENDIF
318
319C
320C if I3,I4 simultaneously even I4 must be larger/equal I3
321C -------------------------------------------------------
322 IF(mod(i3,2).EQ.0.AND.mod(i4,2).EQ.0.AND.i3*i4.GT.0.AND.i3.GT.i4) THEN
323 ! flip IDs
324 buf_i = i3
325 i3 = i4
326 i4 = buf_i
327
328 ! flip 4-vectors
329 buf(0:3) = p(0:3,3)
330 p(0:3,3) = p(0:3,4)
331 p(0:3,4) = buf(0:3)
332 ENDIF
333
334
335 if(testujemy) write(*,*) 'doszlimy do case-a ',i1,i2,i3,i4
336
337 !
338 ! FINALLY select appropriate function
339
340C the ANS=0.0 is not needed unless there is something wrong in list for CASE below
341 ans=0.0
342
343 SELECT CASE(i1+i2)
344 CASE(0) ! UUX DDX CCX SSX INITIAL STATE
345 if(testujemy) write(*,*) 'doszlimy do 0 case-a ',i1,i2,i3,i4
346 IF(abs(i1).EQ.1) CALL ddx(p,i3,i4,h1,h2,key,ans)
347 IF(abs(i1).EQ.2) CALL uux(p,i3,i4,h1,h2,key,ans)
348 IF(abs(i1).EQ.3) CALL ssx(p,i3,i4,h1,h2,key,ans)
349 IF(abs(i1).EQ.4) CALL ccx(p,i3,i4,h1,h2,key,ans)
350 if(testujemy) write(*,*) 'doszlimy do 0 case-a ',i1,i2,i3,i4,ans
351 CASE(1) ! UDX INITIAL STATE
352 IF(abs(i1).EQ.2) CALL udx(p,i3,i4,h1,h2,key,ans)
353 IF(abs(i1).EQ.4) CALL csx(p,i3,i4,h1,h2,key,ans)
354 IF(abs(i1).EQ.3) CALL sux(p,i3,i4,h1,h2,key,ans)
355 CASE(2) ! DD INITIAL STATE
356 if(testujemy) write(*,*) 'doszlimy do 2 case-a ',i1,i2,i3,i4
357 IF(abs(i1).EQ.1) CALL dd(p,i3,i4,h1,h2,key,ans)
358 IF(abs(i1).EQ.4) CALL cux(p,i3,i4,h1,h2,key,ans)
359 IF(abs(i1).EQ.3) CALL sdx(p,i3,i4,h1,h2,key,ans)
360 CASE(3) ! UD INITIAL STATE
361 IF(abs(i1).EQ.2) CALL ud(p,i3,i4,h1,h2,key,ans)
362 IF(abs(i1).EQ.4) CALL cdx(p,i3,i4,h1,h2,key,ans)
363 CASE(4) ! UU INITIAL STATE
364 if(testujemy) write(*,*) 'doszlimy do 4 case-a ',i1,i2,i3,i4
365 IF(abs(i1).EQ.2) CALL uu(p,i3,i4,h1,h2,key,ans)
366 IF(abs(i1).EQ.1) CALL ds(p,i3,i4,h1,h2,key,ans)
367 IF(abs(i1).EQ.3) CALL sd(p,i3,i4,h1,h2,key,ans)
368 CASE(22) ! GLUON D INITIAL STATE
369 CALL gd(p,i3,i4,h1,h2,key,ans)
370 CASE(23) ! GLUON U INITIAL STATE
371 CALL gu(p,i3,i4,h1,h2,key,ans)
372 CASE(24) ! GLUON S INITIAL STATE
373 CALL gd(p,i3,i4,h1,h2,key,ans)
374 CASE(25) ! GLUON C INITIAL STATE
375 CALL gu(p,i3,i4,h1,h2,key,ans)
376 CASE(26) ! GLUON B INITIAL STATE
377 CALL gd(p,i3,i4,h1,h2,key,ans)
378 CASE(42) ! GLUON GLUON INITIAL STATE
379 CALL gg(p,i3,i4,h1,h2,key,ans)
380 CASE(8) ! CC INITIAL STATE
381 CALL cc(p,i3,i4,h1,h2,key,ans)
382 CASE(7) ! CS INITIAL STATE
383 CALL cs(p,i3,i4,h1,h2,key,ans)
384 CASE(5) ! DC INITIAL STATE
385 IF(abs(i1).EQ.1) CALL dc(p,i3,i4,h1,h2,key,ans)
386 IF(abs(i1).EQ.4) CALL cd(p,i3,i4,h1,h2,key,ans)
387 IF(abs(i1).EQ.3) CALL su(p,i3,i4,h1,h2,key,ans)
388 IF(abs(i1).EQ.2) CALL us(p,i3,i4,h1,h2,key,ans)
389 CASE(6) ! SS INITIAL STATE
390 if(testujemy) write(*,*) 'doszlismy do cu i1,i2=',i1,i2,i3,i4
391 IF(abs(i1).EQ.3) CALL ss(p,i3,i4,h1,h2,key,ans)
392 IF(abs(i1).EQ.4) CALL cu(p,i3,i4,h1,h2,key,ans)
393 CASE(-2) ! UCX INITIAL STATE
394 if(testujemy) write(*,*) 'doszlimy do -2 case-a ',i1,i2,i3,i4
395 IF(abs(i1).EQ.2) CALL ucx(p,i3,i4,h1,h2,key,ans)
396 IF(abs(i1).EQ.1) CALL dsx(p,i3,i4,h1,h2,key,ans)
397 CASE(-1) ! USX INITIAL STATE
398 if(testujemy) write(*,*) 'doszlimy do -1 case-a ',i1,i2,i3,i4
399 IF(abs(i1).EQ.2) CALL usx(p,i3,i4,h1,h2,key,ans)
400 IF(abs(i1).EQ.3) CALL scx(p,i3,i4,h1,h2,key,ans)
401 CASE(-3) ! DCX INITIAL STATE
402 if(testujemy) write(*,*) 'doszlimy do -3 case-a ',i1,i2,i3,i4
403 CALL dcx(p,i3,i4,h1,h2,key,ans)
404 CASE DEFAULT
405 ! WRITE(*,*) "VBFDISTR: UNSUPPORTED PROCESS:",I1,I2
406 ans=0.0
407 END SELECT
408 IF(i3.NE.i4) ans=ans/2.d0 ! we divide by 2 because we do not order I3,I4
409 ! but we will take both I3,I4 and I4,I3
410 vbfdistr=ans
411
412 END FUNCTION vbfdistr