C++InterfacetoTauola
VBF_distr.f
1  REAL*8 FUNCTION vbfdistr(ID1,ID2,ID3,ID4,HH1,HH2,PP,KEYIN)
2 C***************************************************************************
3 C* CALCULATES (matrix element)^2 for
4 C I1 I2 -> I3 I4 TAU+TAU- with given tau polarizations H1 H2
5 C for a given set of momenta P
6 C I1,...,I4 are PDG particle codes
7 C if I3 or I4=0, then summed over final jet flavors
8 C H1 and H2 are tau+ and tau- helicities R: 1, L: -1, R+L: 0
9 C KeyIn=0: no Higgs
10 C KeyIn=1: SM Higgs
11 C KeyIn=2: no Higgs
12 C KeyIn=3: SM Higgs
13 C KeyIn=4,5: non-standard state
14 C***************************************************************************
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
56 C ------------
57 C Fast track for cases where results must be zero.
58 C ------------
59 
60 C bayron number conservation:
61 C --------------------------
62 
63  IF(i1+i2.EQ.42.AND.i3+i4.EQ.0) THEN
64 C do nothing
65  ELSEIF(i3+i4.EQ.42.AND.i1+i2.EQ.0) THEN
66 C 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 
87 C charge conservation
88 C -------------------
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 
97 C zero or two gluons
98 C ------------------
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 
112 C if gluons are present extra consitency check.
113 C important, that this is before 5 to 1 replacement for configs with gluons.
114 C --------------------------------------------------------------------------
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 
123 C =============================================
124 C now we have a chance that result is non zero.
125 C we copy configuration to local variables
126 C =============================================
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 
137 C ===========================================
138 C REARRANING order of 4-vectors and ID's
139 C Parity reflection
140 C ===========================================
141 
142 
143 C treat quarks 5 as 1 if gluons present. It is possible
144 C --------------------------------------
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
152 C in other cases processes with b-quarks are not yet installed.
153 C -------------------------------------------------------------
154  vbfdistr=0.0
155  RETURN
156  ENDIF
157 
158 
159  if(testujemy) write(*,*) 'doszlimy do stepX',i1,i2, i3, i4
160 
161 C If I1 I2 particle/antiparticle (no gluon), |I1|>|I2| enforce
162 C -------------------------------------------------------------
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 
179 C C-reflection if all quarks antipartices
180 C ---------------------------------------
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)
199 C 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 
208 C for incoming particle antiparticle larger |id| must be first,
209 C EXCEPTION: ID1 ID2 = (1,-2) or (-2,1) There is UDX.f but no DUX.f file
210 C -----------------------------------------------------------------------
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 
232 C 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 
244 C First incoming can not be antiparticle, second can not be lone gluon
245 C ---------------------------------------------------------------------
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 
259 C 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 
272 C ================
273 C NOW FINAL STATES
274 C ================
275 
276 C I3 never negative and I4 never alone 21
277 C ---------------------------------------
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 
290 C sole posibility <I3 even I4 odd> if both positive and non gluon
291 C ---------------------------------------------------------------
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 
305 C if I3,I4 simultaneously odd I4 must be larger/equal I3
306 C -------------------------------------------------------
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 
319 C
320 C if I3,I4 simultaneously even I4 must be larger/equal I3
321 C -------------------------------------------------------
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 
340 C 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
Generated on Tue Mar 12 2024 22:57:51 for C++InterfacetoTauola by   doxygen 1.8.7