C++ Interface to Tauola
tauola-BBB/jetset-F/tauface-jetset.f
1 SUBROUTINE tauola(MODE,KEYPOL)
2C *************************************
3C general tauola interface, should work in every case until
4C hepevt is OK, does not check if hepevt is 'clean'
5C in particular will decay decayed taus...
6C only longitudinal spin effects are included.
7C in W decay v-a vertex is assumed
8C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001
9C this is the hepevt class in old style. No d_h_ class pre-name
10 INTEGER NMXHEP
11 parameter(nmxhep=4000)
12 real*8 phep, vhep ! to be real*4/ *8 depending on host
13 INTEGER nevhep,nhep,isthep,idhep,jmohep,
14 $ jdahep
15 COMMON /hepevt/
16 $ nevhep, ! serial number
17 $ nhep, ! number of particles
18 $ isthep(nmxhep), ! status code
19 $ idhep(nmxhep), ! particle ident KF
20 $ jmohep(2,nmxhep), ! parent particles
21 $ jdahep(2,nmxhep), ! childreen particles
22 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
23 $ vhep(4,nmxhep) ! vertex [mm]
24* ----------------------------------------------------------------------
25 LOGICAL qedrad
26 COMMON /phoqed/
27 $ qedrad(nmxhep) ! Photos flag
28* ----------------------------------------------------------------------
29 SAVE hepevt,phoqed
30
31
32 COMMON /taupos/ np1, np2
33 real*4 phoi(4),phof(4)
34 double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4)
35 COMMON / momdec / q1,q2,p1,p2,p3,p4
36* tauola, photos and jetset overall switches
37 COMMON /libra/ jak1,jak2,itdkrc,ifphot,ifhadm,ifhadp
38 real*4 rrr(1)
39 LOGICAL IFPSEUDO
40 common /pseudocoup/ csc,ssc
41 real*4 csc,ssc
42 save pseudocoup
43 COMMON / inout / inut,iout
44
45C to switch tau polarization OFF in taus
46 dimension pol1(4), pol2(4)
47 double precision POL1x(4), POL2x(4)
48 INTEGER ION(3)
49 DATA pol1 /0.0,0.0,0.0,0.0/
50 DATA pol2 /0.0,0.0,0.0,0.0/
51 DATA pi /3.141592653589793238462643d0/
52
53C store decay vertexes
54 dimension imother(20)
55 INTEGER KFHIGGS(3)
56
57C store daughter pointers
58 INTEGER ISON
59 COMMON /isons_tau/ison(2)
60 SAVE /isons_tau/
61
62 IF(mode.EQ.-1) THEN
63C ***********************
64
65 jak1 = 0 ! decay mode first tau
66 jak2 = 0 ! decay mode second tau
67 itdkrc=1.0 ! switch of radiative corrections in decay
68 ifphot=1.0 ! PHOTOS switch
69 ifhadm=1.0
70 ifhadp=1.0
71 pol=1.0 ! tau polarization dipswitch must be 1 or 0
72
73 kfhiggs(1) = 25
74 kfhiggs(2) = 35
75 kfhiggs(3) = 36
76 kfhigch = 37
77 kfz0 = 23
78 kfgam = 22
79 kftau = 15
80 kfnue = 16
81C couplings of the 'pseudoscalar higgs' as in CERN-TH/2003-166
82 psi=0.5*pi ! 0.15*PI
83 xmtau=1.777 ! tau mass
84 xmh=120 ! higgs boson mass
85 betah=sqrt(1d0-4*xmtau**2/xmh**2)
86 csc=cos(psi)*betah
87 ssc=sin(psi)
88C write(*,*) ' scalar component=',csc,' pseudo-scalar component=',ssc
89 IF (ifphot.EQ.1) CALL phoini ! this if PHOTOS was not initialized earlier
90 CALL inietc(jak1,jak2,itdkrc,ifphot)
91 CALL inimas
92 CALL iniphx(0.01d0)
93 CALL initdk
94C activation of pi0 and eta decays: (-1,1) means on, (-1,0) off
95 ion(1)=1
96 ion(2)=1
97 ion(3)=1
98 CALL taupi0(-1,1,ion)
99 CALL dekay(-1,pol1x)
100 WRITE(iout,7001) pol,psi,ion(1),ion(2),ion(3)
101 ELSEIF(mode.EQ.0) THEN
102C ***********************
103C
104C..... find tau-s and fill common block /TAUPOS/
105C this is to avoid LUND history fillings. This call is optional
106 CALL phyfix(nstop,nstart)
107C clear mothers of the previous event
108 DO ii=1,20
109 imother(ii)=0
110 ENDDO
111
112 DO ii=1,2
113 ison(ii)=0
114 ENDDO
115C ... and to find mothers giving taus.
116 ndec = 0
117C(BPK)--> LOOK FOR MOTHER, CHECK THAT IT IS NOT THE HISTORY ENTRY (E.G. MSTP(128)=0)
118 DO i=nstart,nhep
119 IF(abs(idhep(i)).EQ.kftau.AND.isthep(i).EQ.1.AND.
120 $ (isthep(i).GE.125.OR.isthep(i).LT.120)) THEN
121 imoth=jmohep(1,i)
122 DO WHILE (abs(idhep(imoth)).EQ.kftau) ! KEEP WALKING UP
123 imoth=jmohep(1,imoth)
124 ENDDO
125 IF (isthep(imoth).EQ.3.OR.
126 $ (isthep(imoth).GE.120.AND.isthep(imoth).LE.125)) THEN
127 DO j=nstart,nhep ! WE HAVE WALKED INTO HARD RECORD
128 IF (idhep(j).EQ.idhep(imoth).AND.
129 $ jmohep(1,j).EQ.imoth.AND.
130 $ isthep(j).EQ.2) THEN
131 jmoth=j
132 GOTO 66
133 ENDIF
134 ENDDO
135 ELSE
136 jmoth=imoth
137 ENDIF
138 66 CONTINUE
139 DO ii=1,ndec
140 IF(jmoth.EQ.imother(ii)) GOTO 9999
141 ENDDO
142C(BPK)--<
143 ndec=ndec+1
144 imother(ndec)= jmoth
145 ENDIF
146 9999 CONTINUE
147 ENDDO
148
149C ... taus of every mother are treated in this main loop
150 DO ii=1,ndec
151 im=imother(ii)
152 ncount=0
153 np1=0
154 np2=0
155
156
157C(BPK)-->
158C CORRECTING HEPEVT IS OUT OF QUESTION AT THIS POINT..
159 im0=im
160 IF (idhep(jmohep(1,im0)).EQ.idhep(im0)) im0=jmohep(1,im0)
161 isel=-1
162 DO i=nstart,nhep
163 IF (isthep(i).EQ.3.OR.
164 $ (isthep(i).GE.120.AND.isthep(i).LE.125)) THEN ! HARD RECORD
165 GOTO 76
166 ENDIF
167 imoth=jmohep(1,i)
168 DO WHILE (idhep(imoth).EQ.idhep(i).OR.abs(idhep(imoth)).EQ.kftau) ! KEEP WALKING UP
169 imoth=jmohep(1,imoth)
170 ENDDO
171 IF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.-1) THEN
172 ison(1)=i
173 isel=0
174 ELSEIF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.0) THEN
175 ison(2)=i
176 ELSEIF ((imoth.NE.im0.AND.imoth.NE.im).AND.isel.EQ.0) THEN
177 isel=1
178 GOTO 77
179 ENDIF
180 76 CONTINUE
181 ENDDO
182 77 CONTINUE
183C(BPK)--<
184
185
186C ... we correct HEPEVT (fix developped with Catherine BISCARAT)
187c IF (JDAHEP(2,IM).EQ.0) THEN ! ID of second daughter was missing
188c ISECU=1
189c DO I=JDAHEP(1,IM)+1,NHEP ! OK lets look for it
190c IF (JMOHEP(1,I).EQ.IM.AND.ISECU.EQ.1) THEN ! we have found one
191c JDAHEP(2,IM)=I
192c ELSEIF (JMOHEP(1,I).EQ.IM.AND.ISECU.NE.1) THEN ! we have found one after there
193c JDAHEP(2,IM)=0 ! was something else, lets kill game
194c ENDIF
195c IF (JMOHEP(1,I).NE.IM) ISECU=0 ! other stuff starts
196c ENDDO
197c ENDIF
198
199C ... we check whether there are just two or more tau-likes
200 DO i=ison(1),ison(2)
201 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) ncount=ncount+1
202 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) ncount=ncount+1
203 ENDDO
204
205C ... if there will be more we will come here again
206 666 CONTINUE
207
208C(BPK)-->
209 DO i=max(np1+1,ison(1)),ison(2)
210C(BPK)--<
211 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) np1=i
212 ENDDO
213C(BPK)-->
214 DO i=max(np2+1,ison(1)),ison(2)
215C(BPK)--<
216 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) np2=i
217 ENDDO
218 DO i=1,4
219 p1(i)= phep(i,np1) !momentum of tau+
220 p2(i)= phep(i,np2) !momentum of tau-
221 q1(i)= p1(i)+p2(i)
222 ENDDO
223
224 pol1(3)= 0d0
225 pol2(3)= 0d0
226
227 IF(keypol.EQ.1) THEN
228c.....include polarisation effect
229 CALL ranmar(rrr,1)
230
231 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
232 $ idhep(im).EQ.kfhiggs(3)) THEN ! case of Higgs
233 IF(rrr(1).LT.0.5) THEN
234 pol1(3)= pol
235 pol2(3)=-pol
236 ELSE
237 pol1(3)=-pol
238 pol2(3)= pol
239 ENDIF
240 ELSEIF((idhep(im).EQ.kfz0).OR.(idhep(im).EQ.kfgam)) THEN ! case of gamma/Z
241C there is no angular dependence in gamma/Z polarization
242C there is no s-dependence in gamma/Z polarization at all
243C there is even no Z polarization in any form
244C main reason is that nobody asked ...
245C but it is prepared and longitudinal correlations
246C can be included up to KORALZ standards
247
248 polz0=plzapx(.true.,im,np1,np2)
249 IF(rrr(1).LT.polz0) THEN
250 pol1(3)= pol
251 pol2(3)= pol
252 ELSE
253 pol1(3)=-pol
254 pol2(3)=-pol
255 ENDIF
256 ELSEIF(idhep(np1).EQ.-idhep(np2))THEN ! undef orig. only s-dep poss.
257 polz0=plzapx(.true.,im,np1,np2)
258 IF(rrr(1).LT.polz0) THEN
259 pol1(3)= pol
260 pol2(3)= pol
261 ELSE
262 pol1(3)=-pol
263 pol2(3)=-pol
264 ENDIF
265 ELSEIF(abs(idhep(im)).EQ.kfhigch) THEN ! case of charged Higgs
266 pol1(3)= pol
267 pol2(3)= pol
268 ELSE ! case of W+ or W-
269 pol1(3)= -pol
270 pol2(3)= -pol
271 ENDIF
272c.....include polarisation effect
273 ENDIF
274
275 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
276 $ idhep(im).EQ.kfhiggs(3)) THEN
277 IF(idhep(np1).EQ.-kftau .AND.
278 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep) .AND.
279 $ idhep(np2).EQ. kftau .AND.
280 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)
281 $ ) THEN
282 IF (idhep(im).EQ.kfhiggs(1)) THEN
283 ifpseudo= .false.
284 ELSEIF (idhep(im).EQ.kfhiggs(2)) THEN
285 ifpseudo= .false.
286 ELSEIF (idhep(im).EQ.kfhiggs(3)) THEN
287 ifpseudo= .true.
288 ELSE
289 WRITE(*,*) 'Warning from TAUOLA:'
290 WRITE(*,*) 'I stop this run, wrong IDHEP(IM)=',idhep(im)
291 stop
292 ENDIF
293 CALL spinhiggs(im,np1,np2,ifpseudo,pol1,pol2)
294 IF (ifphot.EQ.1) CALL photos(im) ! Bremsstrahlung in Higgs decay
295 ! AFTER adding taus !!
296
297
298 ENDIF
299 ELSE
300 IF(idhep(np1).EQ.-kftau.AND.
301 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep)) THEN
302C here check on if NP1 was not decayed should be verified
303 CALL dexay(1,pol1)
304 IF (ifphot.EQ.1) CALL photos(np1)
305 CALL taupi0(0,1,ion)
306 ENDIF
307
308 IF(idhep(np2).EQ. kftau.AND.
309 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)) THEN
310C here check on if NP2 was not decayed should be added
311 CALL dexay(2,pol2)
312 IF (ifphot.EQ.1) CALL photos(np2)
313 CALL taupi0(0,2,ion)
314 ENDIF
315 ENDIF
316 ncount=ncount-2
317 IF (ncount.GT.0) GOTO 666
318 ENDDO
319
320 ELSEIF(mode.EQ.1) THEN
321C ***********************
322C
323 CALL dexay(100,pol1)
324 CALL dekay(100,pol1x)
325 WRITE(iout,7002)
326 ENDIF
327C *****
328 7001 FORMAT(///1x,15(5h*****)
329 $ /,' *', 25x,'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
330 $ /,' *', 25x,'*****VERSION 1.21, September 2005******',9x,1h*,
331 $ /,' *', 25x,'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
332 $ /,' *', 25x,'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
333 $ /,' *', 25x,'****** Z. Was, M. Worek ***************',9x,1h*,
334 $ /,' *', 25x,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
335 $ /,' *', 25x,'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
336 $ /,' *', 25x,'****** are warmly acknowledged ********',9x,1h*,
337 $ /,' *', 25x,' ',9x,1h*,
338 $ /,' *', 25x,'********** INITIALIZATION ************',9x,1h*,
339 $ /,' *',f20.5,5x,'tau polarization switch must be 1 or 0 ',9x,1h*,
340 $ /,' *',f20.5,5x,'Higs scalar/pseudo mix CERN-TH/2003-166',9x,1h*,
341 $ /,' *',i10, 15x,'PI0 decay switch must be 1 or 0 ',9x,1h*,
342 $ /,' *',i10, 15x,'ETA decay switch must be 1 or 0 ',9x,1h*,
343 $ /,' *',i10, 15x,'K0S decay switch must be 1 or 0 ',9x,1h*,
344 $ /,1x,15(5h*****)/)
345
346 7002 FORMAT(///1x,15(5h*****)
347 $ /,' *', 25x,'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
348 $ /,' *', 25x,'*****VERSION 1.21, September2005 ******',9x,1h*,
349 $ /,' *', 25x,'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
350 $ /,' *', 25x,'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
351 $ /,' *', 25x,'****** Z. Was, M. Worek ***************',9x,1h*,
352 $ /,' *', 25x,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
353 $ /,' *', 25x,'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
354 $ /,' *', 25x,'****** are warmly acknowledged ********',9x,1h*,
355 $ /,' *', 25x,'****** END OF MODULE OPERATION ********',9x,1h*,
356 $ /,1x,15(5h*****)/)
357
358 END
359
360 SUBROUTINE spinhiggs(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
361 LOGICAL IFPSEUDO
362 real*8 hh1,hh2,wthiggs
363 dimension pol1(4), pol2(4),hh1(4),hh2(4), rrr(1)
364! CALL DEXAY(1,POL1) ! Kept for tests
365! CALL DEXAY(2,POL2) ! Kept for tests
366 INTEGER ION(3)
367 10 CONTINUE
368 CALL ranmar(rrr,1)
369 CALL dekay(1,hh1)
370 CALL dekay(2,hh2)
371 wt=wthiggs(ifpseudo,hh1,hh2)
372 IF (rrr(1).GT.wt) GOTO 10
373 CALL dekay(1+10,hh1)
374 CALL taupi0(0,1,ion)
375 CALL dekay(2+10,hh2)
376 CALL taupi0(0,2,ion)
377 END
378 FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
379 LOGICAL IFPSEUDO
380 common /pseudocoup/ csc,ssc
381 real*4 csc,ssc
382 save pseudocoup
383 real*8 hh1(4),hh2(4),r(4,4),wthiggs
384 DO k=1,4
385 DO l=1,4
386 r(k,l)=0
387 ENDDO
388 ENDDO
389 wthiggs=0d0
390
391 r(4,4)= 1d0 ! unpolarized part
392 r(3,3)=-1d0 ! longitudinal
393 ! other missing
394 IF (ifpseudo) THEN
395 r(1,1)=-1
396 r(2,2)= -1
397 r(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
398 r(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
399 r(1,2)=2*csc*ssc/(csc**2+ssc**2)
400 r(2,1)=-2*csc*ssc/(csc**2+ssc**2)
401 ELSE
402 r(1,1)=1
403 r(2,2)=1
404 ENDIF
405
406
407
408 DO k=1,4
409 DO l=1,4
410 wthiggs=wthiggs+r(k,l)*hh1(k)*hh2(l)
411 ENDDO
412 ENDDO
413 wthiggs=wthiggs/4d0
414 END
415
416 FUNCTION plzapx(HOPEin,IM0,NP1,NP2)
417C IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block.
418C the purpose of this routine is to calculate polarization of Z along tau direction.
419C this is highly non-trivial due to necessity of reading infromation from hard process
420C history in HEPEVT, which is often written not up to the gramatic rules.
421 real*8 plzapx,plzap0,svar,costhe,sini,sfin,zprop2,
422 $ p1(4),p2(4),q1(4),q2(4),qq(4),ph(4),pd1(4),pd2(4),
423 $ pq1(4),pq2(4),pb(4),pa(4)
424 INTEGER IM
425 LOGICAL HOPE,HOPEin
426C this is the hepevt class in old style. No d_h_ class pre-name
427 INTEGER NMXHEP
428 parameter(nmxhep=4000)
429 real*8 phep, vhep ! to be real*4/ *8 depending on host
430 INTEGER nevhep,nhep,isthep,idhep,jmohep,
431 $ jdahep
432 COMMON /hepevt/
433 $ nevhep, ! serial number
434 $ nhep, ! number of particles
435 $ isthep(nmxhep), ! status code
436 $ idhep(nmxhep), ! particle ident KF
437 $ jmohep(2,nmxhep), ! parent particles
438 $ jdahep(2,nmxhep), ! childreen particles
439 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
440 $ vhep(4,nmxhep) ! vertex [mm]
441* ----------------------------------------------------------------------
442 LOGICAL qedrad
443 COMMON /phoqed/
444 $ qedrad(nmxhep) ! Photos flag
445* ----------------------------------------------------------------------
446 SAVE hepevt,phoqed
447
448
449
450C(BPK)--> BROTHERS OF TAU ALREADY FOUND
451 INTEGER ISON
452 COMMON /isons_tau/ison(2)
453C(BPK)--<
454C >>
455C >> STEP 1: find where are particles in hepevent and pick them up
456C >>
457 hope=hopein
458C sometimes shade Z of Z is its mother ...
459 im=im0
460 im00=jmohep(1,im0)
461C to protect against check on mother of beam particles.
462 IF (im00.GT.0) THEN
463 IF (idhep(im0).EQ.idhep(im00)) im=jmohep(1,im0)
464 ENDIF
465C
466C find (host generator-level) incoming beam-bare-particles which form Z and co.
467 imo1=jmohep(1,im)
468 imo2=jmohep(2,im)
469
470C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS
471 im00=imo1
472 IF (isthep(im00).EQ.120) THEN
473 imo1=jmohep(1,im00)
474 imo2=jmohep(2,im00)
475 ENDIF
476C(BPK)--<
477
478 iffull=0
479C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first' in hepevt.
480 IF (imo1.EQ.0.AND.imo2.EQ.0) THEN
481 imo1=jmohep(1,np1)
482C(BPK)-->
483 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1) ! PROTECT AGAINST COPIES
484C(BPK)--<
485 imo2=jmohep(2,np1)
486C(BPK)-->
487 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2) ! PROTECT AGAINST COPIES
488C(BPK)--<
489 iffull=1
490C case when it was like qq --> tau+tau- gammas and qq were NOT 'first' in hepevt.
491
492 ELSEIF (idhep(im).NE.22.AND.idhep(im).NE.23) THEN
493 imo1=jmohep(1,np1)
494C(BPK)-->
495 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1) ! PROTECT AGAINST COPIES
496C(BPK)--<
497 imo2=jmohep(2,np1)
498C(BPK)-->
499 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2) ! PROTECT AGAINST COPIES
500C(BPK)--<
501 iffull=1
502 ENDIF
503
504
505C and check if it really happened
506 IF (imo1.EQ.0) hope=.false.
507 IF (imo2.EQ.0) hope=.false.
508 IF (imo1.EQ.imo2) hope=.false.
509
510C
511 DO i=1,4
512 q1(i)= phep(i,np1) !momentum of tau+
513 q2(i)= phep(i,np2) !momentum of tau-
514 ENDDO
515
516C corrections due to possible differences in 4-momentum of shadow vs true Z.
517C(BPK)-->
518 IF (im.EQ.jmohep(1,im0).AND.
519 $ (idhep(im).EQ.22.OR.idhep(im).EQ.23)) THEN
520 DO k=1,4
521 pb(k)=phep(k,im)
522 pa(k)=phep(k,im0)
523 ENDDO
524C(BPK)--<
525
526 CALL bostdq( 1,pa, q1, q1)
527 CALL bostdq( 1,pa, q2, q2)
528 CALL bostdq(-1,pb, q1, q1)
529 CALL bostdq(-1,pb, q2, q2)
530
531 ENDIF
532
533 DO i=1,4
534 qq(i)= q1(i)+q2(i) !momentum of Z
535 IF (hope) p1(i)=phep(i,imo1) !momentum of beam1
536 IF (hope) p2(i)=phep(i,imo2) !momentum of beam2
537 ph(i)=0d0
538 pd1(i)=0d0
539 pd2(i)=0d0
540 ENDDO
541! These momenta correspond to quarks, gluons photons or taus
542 idfq1=idhep(np1)
543 idfq2=idhep(np2)
544 IF (hope) idfp1=idhep(imo1)
545 IF (hope) idfp2=idhep(imo2)
546
547 svar=qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2
548 IF (.NOT.hope) THEN
549C options which may be useful in some cases of two heavy boson production
550C need individual considerations. To be developed.
551C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
552C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
553 plzapx=0.5 ! pure gamma
554 RETURN
555 ENDIF
556C >>
557C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles
558C >>
559C let us define beginning and end of particles which are produced in parallel to Z
560C in principle following should work
561
562C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS
563 nx1=jdahep(1,im00)
564 nx2=jdahep(2,im00)
565C but ...
566 inbr=im ! OK, HARD RECORD Z/GAMMA
567 IF (iffull.EQ.1) inbr=np1 ! OK, NO Z/GAMMA
568 IF (idhep(jmohep(1,inbr)).EQ.idhep(inbr)) inbr=jmohep(1,inbr) ! FORCE HARD RECORD
569C(BPK)--<
570 IF(nx1.EQ.0.OR.nx2.EQ.0) THEN
571 nx1=inbr
572 nx2=inbr
573 DO k=1,inbr-1
574 IF(jmohep(1,inbr-k).EQ.jmohep(1,inbr)) THEN
575 nx1=inbr-k
576 ELSE
577 GOTO 7
578 ENDIF
579 ENDDO
580 7 CONTINUE
581
582 DO k=inbr+1,nhep
583 IF(jmohep(1,k).EQ.jmohep(1,inbr)) THEN
584 nx2=k
585 ELSE
586 GOTO 8
587 ENDIF
588 ENDDO
589 8 CONTINUE
590 ENDIF
591
592C case of annihilation of two bosons is hopeles
593 IF (abs(idfp1).GE.20.AND.abs(idfp2).GE.20) hope=.false.
594C case of annihilation of non-matching flavors is hopeless
595 IF (abs(idfp1).LE.20.AND.abs(idfp2).LE.20.AND.idfp1+idfp2.NE.0)
596 $ hope=.false.
597 IF (.NOT.hope) THEN
598C options which may be useful in some cases of two heavy boson production
599C need individual considerations. To be developed.
600C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
601C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
602 plzapx=0.5 ! pure gamma
603 RETURN
604 ENDIF
605 IF (abs(idfp1).LT.20) ide= idfp1
606 IF (abs(idfp2).LT.20) ide=-idfp2
607
608
609C >>
610C >> STEP 3 we combine gluons, photons into incoming effective beams
611C >>
612
613C in the following we will ignore the possibility of photon emission from taus
614C however at certain moment it will be necessary to take care of
615
616 DO l=1,4
617 pd1(l)=p1(l)
618 pd2(l)=p2(l)
619 ENDDO
620
621 DO l=1,4
622 pq1(l)=q1(l)
623 pq2(l)=q2(l)
624 ENDDO
625
626 iflav=min(abs(idfp1),abs(idfp2))
627
628*--------------------------------------------------------------------------
629c IFLAV=min(ABS(IDFP1),ABS(IDFP2))
630c that means that always quark or lepton i.e. process like
631c f g(gamma) --> f Z0 --> tau tau
632c we glue fermions to effective beams that is f f --> Z0 --> tau tau
633c with gamma/g emission from initial fermion.
634*---------------------------------------------------------------------------
635
636 IF (abs(idfp1).GE.20) THEN
637 DO k=nx1,nx2
638 idp=idhep(k)
639 IF (abs(idp).EQ.iflav) THEN
640 DO l=1,4
641 pd1(l)=-phep(l,k)
642 ENDDO
643 ENDIF
644 ENDDO
645 ENDIF
646
647 IF (abs(idfp2).GE.20) THEN
648 DO k=nx1,nx2
649 idp=idhep(k)
650 IF (abs(idp).EQ.iflav) THEN
651 DO l=1,4
652 pd2(l)=-phep(l,k)
653 ENDDO
654 ENDIF
655 ENDDO
656 ENDIF
657
658C if first beam was boson: gluing
659
660 IF (abs(idfp1).GE.20) THEN
661 DO l=1,4
662 ph(l)=p1(l)
663 ENDDO
664 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
665 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
666 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
667 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
668 IF (xm1.LT.xm2) THEN
669 DO l=1,4
670 pd1(l)=pd1(l)+p1(l)
671 ENDDO
672 ELSE
673 DO l=1,4
674 pd2(l)=pd2(l)+p1(l)
675 ENDDO
676 ENDIF
677 ENDIF
678
679C if second beam was boson: gluing
680
681
682 IF (abs(idfp2).GE.20) THEN
683 DO l=1,4
684 ph(l)=p2(l)
685 ENDDO
686 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
687 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
688 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
689 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
690 IF (xm1.LT.xm2) THEN
691 DO l=1,4
692 pd1(l)=pd1(l)+p2(l)
693 ENDDO
694 ELSE
695 DO l=1,4
696 pd2(l)=pd2(l)+p2(l)
697 ENDDO
698 ENDIF
699 ENDIF
700
701C now spectators ...
702
703C(BPK)-->
704 nph1=np1
705 nph2=np2
706 IF (idhep(jmohep(1,np1)).EQ.idhep(np1)) nph1=jmohep(1,np1) ! SHOULD PUT US IN HARD REC.
707 IF (idhep(jmohep(1,np2)).EQ.idhep(np2)) nph2=jmohep(1,np2) ! SHOULD PUT US IN HARD REC.
708C(BPK)--<
709
710 DO k=nx1,nx2
711 IF (abs(idhep(k)).NE.iflav.AND.k.NE.im.AND.
712C(BPK)-->
713 $ k.NE.nph1.AND.k.NE.nph2) THEN
714C(BPK)--<
715 IF(idhep(k).EQ.22.AND.iffull.EQ.1) THEN
716 DO l=1,4
717 ph(l)=phep(l,k)
718 ENDDO
719 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
720 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
721 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
722 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
723 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
724 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
725 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
726 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
727
728
729 sini=abs((pd1(4)+pd2(4)-ph(4))**2-(pd1(3)+pd2(3)-ph(3))**2
730 $ -(pd1(2)+pd2(2)-ph(2))**2-(pd1(1)+pd2(1)-ph(1))**2)
731 sfin=abs((pd1(4)+pd2(4) )**2-(pd1(3)+pd2(3) )**2
732 $ -(pd1(2)+pd2(2) )**2-(pd1(1)+pd2(1) )**2)
733
734 facini=zprop2(sini)
735 facfin=zprop2(sfin)
736
737 xm1=xm1/facini
738 xm2=xm2/facini
739 xm3=xm3/facfin
740 xm4=xm4/facfin
741
742 xm=min(xm1,xm2,xm3,xm4)
743 IF (xm1.EQ.xm) THEN
744 DO l=1,4
745 pd1(l)=pd1(l)-ph(l)
746 ENDDO
747 ELSEIF (xm2.EQ.xm) THEN
748 DO l=1,4
749 pd2(l)=pd2(l)-ph(l)
750 ENDDO
751 ELSEIF (xm3.EQ.xm) THEN
752 DO l=1,4
753 q1(l)=pq1(l)+ph(l)
754 ENDDO
755 ELSE
756 DO l=1,4
757 q2(l)=pq2(l)+ph(l)
758 ENDDO
759 ENDIF
760 ELSE
761 DO l=1,4
762 ph(l)=phep(l,k)
763 ENDDO
764 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
765 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
766 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
767 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
768 IF (xm1.LT.xm2) THEN
769 DO l=1,4
770 pd1(l)=pd1(l)-ph(l)
771 ENDDO
772 ELSE
773 DO l=1,4
774 pd2(l)=pd2(l)-ph(l)
775 ENDDO
776 ENDIF
777 ENDIF
778 ENDIF
779 ENDDO
780
781
782C >>
783C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in
784c >> effective outcoming taus
785C >>
786C let us define beginning and end of particles which are produced in
787c parallel to tau
788
789
790
791C find outcoming particles which come from Z
792
793
794
795
796C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER
797 IF (abs(idhep(im0)).EQ.22.OR.abs(idhep(im0)).EQ.23) THEN
798 DO k=ison(1),ison(2)
799 IF(abs(idhep(k)).EQ.22) THEN
800C(BPK)--<
801
802 do l=1,4
803 ph(l)=phep(l,k)
804 enddo
805
806 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
807 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
808 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
809 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
810
811 xm=min(xm3,xm4)
812
813 IF (xm3.EQ.xm) THEN
814 DO l=1,4
815 q1(l)=pq1(l)+ph(l)
816 ENDDO
817 ELSE
818 DO l=1,4
819 q2(l)=pq2(l)+ph(l)
820 ENDDO
821 ENDIF
822 endif
823 enddo
824 ENDIF
825
826
827
828*------------------------------------------------------------------------
829
830
831C out of effective momenta we calculate COSTHE and later polarization
832 CALL angulu(pd1,pd2,q1,q2,costhe)
833
834 plzapx=plzap0(ide,idfq1,svar,costhe)
835 END
836
837 SUBROUTINE angulu(PD1,PD2,Q1,Q2,COSTHE)
838 real*8 pd1(4),pd2(4),q1(4),q2(4),costhe,p(4),qq(4),qt(4)
839C take effective beam which is less massive, it should be irrelevant
840C but in case HEPEVT is particulary dirty may help.
841C this routine calculate reduced system transver and cosine of scattering
842C angle.
843
844 xm1=abs(pd1(4)**2-pd1(3)**2-pd1(2)**2-pd1(1)**2)
845 xm2=abs(pd2(4)**2-pd2(3)**2-pd2(2)**2-pd2(1)**2)
846 IF (xm1.LT.xm2) THEN
847 sign=1d0
848 DO k=1,4
849 p(k)=pd1(k)
850 ENDDO
851 ELSE
852 sign=-1d0
853 DO k=1,4
854 p(k)=pd2(k)
855 ENDDO
856 ENDIF
857C calculate space like part of P (in Z restframe)
858 DO k=1,4
859 qq(k)=q1(k)+q2(k)
860 qt(k)=q1(k)-q2(k)
861 ENDDO
862
863 xmqq=sqrt(qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2)
864
865 qtxqq=qt(4)*qq(4)-qt(3)*qq(3)-qt(2)*qq(2)-qt(1)*qq(1)
866 DO k=1,4
867 qt(k)=qt(k)-qq(k)*qtxqq/xmqq**2
868 ENDDO
869
870 pxqq=p(4)*qq(4)-p(3)*qq(3)-p(2)*qq(2)-p(1)*qq(1)
871 DO k=1,4
872 p(k)=p(k)-qq(k)*pxqq/xmqq**2
873 ENDDO
874C calculate costhe
875 pxp =sqrt(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
876 qtxqt=sqrt(qt(3)**2+qt(2)**2+qt(1)**2-qt(4)**2)
877 pxqt =p(3)*qt(3)+p(2)*qt(2)+p(1)*qt(1)-p(4)*qt(4)
878 costhe=pxqt/pxp/qtxqt
879 costhe=costhe*sign
880 END
881
882 FUNCTION plzap0(IDE,IDF,SVAR,COSTH0)
883C this function calculates probability for the helicity +1 +1 configuration
884C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
885 real*8 plzap0,svar,costhe,costh0
886
887 costhe=costh0
888C >>>>> IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
889C >>>>> of first beam is used by T_GIVIZ0 including sign
890
891 IF (idf.GT.0) THEN
892 CALL initwk(ide,idf,svar)
893 ELSE
894 CALL initwk(-ide,-idf,svar)
895 ENDIF
896 plzap0=t_born(0,svar,costhe,1d0,1d0)
897 $ /(t_born(0,svar,costhe,1d0,1d0)+t_born(0,svar,costhe,-1d0,-1d0))
898
899! PLZAP0=0.5
900 END
901 FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
902C ----------------------------------------------------------------------
903C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
904C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
905C EXAMPLE OF THE METHOD APPLIED THERE
906C INPUT PARAMETERS ARE: SVAR -- transfer
907C COSTHE -- cosine of angle between tau+ and 1st beam
908C TA,TB -- helicity states of tau+ tau-
909C
910C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
911C ----------------------------------------------------------------------
912 IMPLICIT REAL*8(a-h,o-z)
913 COMMON / t_beampm / ene ,amin,amfin,ide,idf
914 real*8 ene ,amin,amfin
915 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
916 & ,xupgi ,xupzi ,xupgf ,xupzf
917 & ,ndiag0,ndiaga,keya,keyz
918 & ,itce,jtce,itcf,jtcf,kolor
919 real*8 ss,poln,t3e,qe,t3f,qf
920 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
921 real*8 seps1,seps2
922C=====================================================================
923 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
924 real*8 swsq,amw,amz,amh,amtop,gammz
925C SWSQ = sin2 (theta Weinberg)
926C AMW,AMZ = W & Z boson masses respectively
927C AMH = the Higgs mass
928C AMTOP = the top mass
929C GAMMZ = Z0 width
930 COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
931 COMPLEX*16 XUPZFP(2),XUPZIP(2)
932 COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
933 COMPLEX*16 PROPA,PROPZ
934 COMPLEX*16 XR,XI
935 COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF
936 COMPLEX*16 XTHING,XVE,XVF,XVEF
937 DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
938 DATA mode0 /-5/
939 DATA ide0 /-55/
940 DATA svar0,cost0 /-5.d0,-6.d0/
941 DATA pi /3.141592653589793238462643d0/
942 DATA seps1,seps2 /0d0,0d0/
943C
944C MEMORIZATION =========================================================
945 IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
946 $ .OR.ide0.NE.ide)THEN
947C
948 keygsw=1
949C ** PROPAGATORS
950 ide0=ide
951 mode0=mode
952 svar0=svar
953 cost0=costhe
954 sinthe=sqrt(1.d0-costhe**2)
955 beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
956C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
957 xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5*beta*(xupzf(1)-xupzf(2))
958 xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5*beta*(xupzf(1)-xupzf(2))
959 xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5*(xupzi(1)-xupzi(2))
960 xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5*(xupzi(1)-xupzi(2))
961C FINAL STATE VECTOR COUPLING
962 xupf =0.5d0*(xupzf(1)+xupzf(2))
963 xupi =0.5d0*(xupzi(1)+xupzi(2))
964 xthing =0d0
965
966 propa =1d0/svar
967 propz =1d0/dcmplx(svar-amz**2,svar/amz*gammz)
968 IF (keygsw.EQ.0) propz=0.d0
969 DO 50 i=1,2
970 DO 50 j=1,2
971 regula= (3-2*i)*(3-2*j) + costhe
972 regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
973 aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
974 azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
975 aborn(i,j)=aphot(i,j)+azett(i,j)
976 aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
977 azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
978 abornm(i,j)=aphotm(i,j)+azettm(i,j)
979 50 CONTINUE
980 ENDIF
981C
982C******************
983C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
984C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
985C* HELICITY CONSERVATION EXPLICITLY OBEYED
986 polar1= (seps1)
987 polar2= (-seps2)
988 born=0d0
989 DO 150 i=1,2
990 helic= 3-2*i
991 DO 150 j=1,2
992 helit=3-2*j
993 factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
994 factom=factor*(1+helit*ta)*(1-helit*tb)
995 factor=factor*(1+helit*ta)*(1+helit*tb)
996
997 born=born+cdabs(aborn(i,j))**2*factor
998C MASS TERM IN BORN
999 IF (mode.GE.1) THEN
1000 born=born+cdabs(abornm(i,j))**2*factom
1001 ENDIF
1002
1003 150 CONTINUE
1004C************
1005 funt=born
1006 IF(funt.LT.0.d0) funt=born
1007
1008C
1009 IF (svar.GT.4d0*amfin**2) THEN
1010C PHASE SPACE THRESHOLD FACTOR
1011 thresh=sqrt(1-4d0*amfin**2/svar)
1012 t_born= funt*svar**2*thresh
1013 ELSE
1014 thresh=0.d0
1015 t_born=0.d0
1016 ENDIF
1017C ZW HERE WAS AN ERROR 19. 05. 1989
1018! write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
1019! write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
1020 END
1021
1022 SUBROUTINE initwk(IDEX,IDFX,SVAR)
1023! initialization routine coupling masses etc.
1024 IMPLICIT REAL*8 (a-h,o-z)
1025 COMMON / t_beampm / ene ,amin,amfin,ide,idf
1026 real*8 ene ,amin,amfin
1027 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
1028 & ,xupgi ,xupzi ,xupgf ,xupzf
1029 & ,ndiag0,ndiaga,keya,keyz
1030 & ,itce,jtce,itcf,jtcf,kolor
1031 real*8 ss,poln,t3e,qe,t3f,qf
1032 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
1033 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
1034 real*8 swsq,amw,amz,amh,amtop,gammz
1035C SWSQ = sin2 (theta Weinberg)
1036C AMW,AMZ = W & Z boson masses respectively
1037C AMH = the Higgs mass
1038C AMTOP = the top mass
1039C GAMMZ = Z0 width
1040C
1041 ene=sqrt(svar)/2
1042 amin=0.511d-3
1043 swsq=0.23147
1044 amz=91.1882
1045 gammz=2.4952
1046 IF (idfx.EQ. 15) then
1047 idf=2 ! denotes tau +2 tau-
1048 amfin=1.77703 !this mass is irrelevant if small, used in ME only
1049 ELSEIF (idfx.EQ.-15) then
1050 idf=-2 ! denotes tau -2 tau-
1051 amfin=1.77703 !this mass is irrelevant if small, used in ME only
1052 ELSE
1053 WRITE(*,*) 'INITWK: WRONG IDFX'
1054 stop
1055 ENDIF
1056
1057 IF (idex.EQ. 11) then !electron
1058 ide= 2
1059 amin=0.511d-3
1060 ELSEIF (idex.EQ.-11) then !positron
1061 ide=-2
1062 amin=0.511d-3
1063 ELSEIF (idex.EQ. 13) then !mu+
1064 ide= 2
1065 amin=0.105659
1066 ELSEIF (idex.EQ.-13) then !mu-
1067 ide=-2
1068 amin=0.105659
1069 ELSEIF (idex.EQ. 1) then !d
1070 ide= 4
1071 amin=0.05
1072 ELSEIF (idex.EQ.- 1) then !d~
1073 ide=-4
1074 amin=0.05
1075 ELSEIF (idex.EQ. 2) then !u
1076 ide= 3
1077 amin=0.02
1078 ELSEIF (idex.EQ.- 2) then !u~
1079 ide=-3
1080 amin=0.02
1081 ELSEIF (idex.EQ. 3) then !s
1082 ide= 4
1083 amin=0.3
1084 ELSEIF (idex.EQ.- 3) then !s~
1085 ide=-4
1086 amin=0.3
1087 ELSEIF (idex.EQ. 4) then !c
1088 ide= 3
1089 amin=1.3
1090 ELSEIF (idex.EQ.- 4) then !c~
1091 ide=-3
1092 amin=1.3
1093 ELSEIF (idex.EQ. 5) then !b
1094 ide= 4
1095 amin=4.5
1096 ELSEIF (idex.EQ.- 5) then !b~
1097 ide=-4
1098 amin=4.5
1099 ELSEIF (idex.EQ. 12) then !nu_e
1100 ide= 1
1101 amin=0.1d-3
1102 ELSEIF (idex.EQ.- 12) then !nu_e~
1103 ide=-1
1104 amin=0.1d-3
1105 ELSEIF (idex.EQ. 14) then !nu_mu
1106 ide= 1
1107 amin=0.1d-3
1108 ELSEIF (idex.EQ.- 14) then !nu_mu~
1109 ide=-1
1110 amin=0.1d-3
1111 ELSEIF (idex.EQ. 16) then !nu_tau
1112 ide= 1
1113 amin=0.1d-3
1114 ELSEIF (idex.EQ.- 16) then !nu_tau~
1115 ide=-1
1116 amin=0.1d-3
1117
1118 ELSE
1119 WRITE(*,*) 'INITWK: WRONG IDEX'
1120 stop
1121 ENDIF
1122
1123C ----------------------------------------------------------------------
1124C
1125C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
1126C
1127C called by : KORALZ
1128C ----------------------------------------------------------------------
1129 itce=ide/iabs(ide)
1130 jtce=(1-itce)/2
1131 itcf=idf/iabs(idf)
1132 jtcf=(1-itcf)/2
1133 CALL t_givizo( ide, 1,aizor,qe,kdumm)
1134 CALL t_givizo( ide,-1,aizol,qe,kdumm)
1135 xupgi(1)=qe
1136 xupgi(2)=qe
1137 t3e = aizol+aizor
1138 xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
1139 xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
1140 CALL t_givizo( idf, 1,aizor,qf,kolor)
1141 CALL t_givizo( idf,-1,aizol,qf,kolor)
1142 xupgf(1)=qf
1143 xupgf(2)=qf
1144 t3f = aizol+aizor
1145 xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
1146 xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
1147C
1148 ndiag0=2
1149 ndiaga=11
1150 keya = 1
1151 keyz = 1
1152C
1153C
1154 RETURN
1155 END
1156
1157 SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1158C ----------------------------------------------------------------------
1159C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
1160C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
1161C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
1162C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
1163C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
1164C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
1165C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
1166C
1167C called by : EVENTE, EVENTM, FUNTIH, .....
1168C ----------------------------------------------------------------------
1169 IMPLICIT REAL*8(a-h,o-z)
1170C
1171 IF(idferm.EQ.0.OR.iabs(idferm).GT.4) GOTO 901
1172 IF(iabs(ihelic).NE.1) GOTO 901
1173 ih =ihelic
1174 idtype =iabs(idferm)
1175 ic =idferm/idtype
1176 lepqua=int(idtype*0.4999999d0)
1177 iupdow=idtype-2*lepqua-1
1178 charge =(-iupdow+2d0/3d0*lepqua)*ic
1179 sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
1180 kolor=1+2*lepqua
1181C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
1182C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1183 RETURN
1184 901 print *,' STOP IN GIVIZO: WRONG PARAMS.'
1185 stop
1186 END
1187 SUBROUTINE phyfix(NSTOP,NSTART)
1188 common/lujets/n,k(4000,5),p(4000,5),v(4000,5)
1189 SAVE /lujets/
1190C NSTOP NSTART : when PHYTIA history ends and event starts.
1191 nstop=0
1192 nstart=1
1193 DO i=1, n
1194 IF(k(i,1).NE.21) THEN
1195 nstop = i-1
1196 nstart= i
1197 GOTO 500
1198 ENDIF
1199 ENDDO
1200 500 CONTINUE
1201 END
1202 SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1203C ----------------------------------------------------------------------
1204C this subroutine fills one entry into the HEPEVT common
1205C and updates the information for affected mother entries
1206C
1207C written by Martin W. Gruenewald (91/01/28)
1208C
1209C called by : ZTOHEP,BTOHEP,DWLUxy
1210C ----------------------------------------------------------------------
1211C
1212C this is the hepevt class in old style. No d_h_ class pre-name
1213C this is the hepevt class in old style. No d_h_ class pre-name
1214 INTEGER NMXHEP
1215 parameter(nmxhep=4000)
1216 real*8 phep, vhep ! to be real*4/ *8 depending on host
1217 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1218 $ jdahep
1219 COMMON /hepevt/
1220 $ nevhep, ! serial number
1221 $ nhep, ! number of particles
1222 $ isthep(nmxhep), ! status code
1223 $ idhep(nmxhep), ! particle ident KF
1224 $ jmohep(2,nmxhep), ! parent particles
1225 $ jdahep(2,nmxhep), ! childreen particles
1226 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1227 $ vhep(4,nmxhep) ! vertex [mm]
1228* ----------------------------------------------------------------------
1229 LOGICAL qedrad
1230 COMMON /phoqed/
1231 $ qedrad(nmxhep) ! Photos flag
1232* ----------------------------------------------------------------------
1233 SAVE hepevt,phoqed
1234
1235
1236 LOGICAL PHFLAG
1237C
1238 real*4 p4(4)
1239C
1240C check address mode
1241 IF (n.EQ.0) THEN
1242C
1243C append mode
1244 ihep=nhep+1
1245 ELSE IF (n.GT.0) THEN
1246C
1247C absolute position
1248 ihep=n
1249 ELSE
1250C
1251C relative position
1252 ihep=nhep+n
1253 END IF
1254C
1255C check on IHEP
1256 IF ((ihep.LE.0).OR.(ihep.GT.nmxhep)) RETURN
1257C
1258C add entry
1259 nhep=ihep
1260 isthep(ihep)=ist
1261 idhep(ihep)=id
1262 jmohep(1,ihep)=jmo1
1263 IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
1264 jmohep(2,ihep)=jmo2
1265 IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
1266 jdahep(1,ihep)=jda1
1267 jdahep(2,ihep)=jda2
1268C
1269 DO i=1,4
1270 phep(i,ihep)=p4(i)
1271C
1272C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
1273 vhep(i,ihep)=0.0
1274 END DO
1275 phep(5,ihep)=pinv
1276C FLAG FOR PHOTOS...
1277 qedrad(ihep)=phflag
1278C
1279C update process:
1280 DO ip=jmohep(1,ihep),jmohep(2,ihep)
1281 IF(ip.GT.0)THEN
1282C
1283C if there is a daughter at IHEP, mother entry at IP has decayed
1284 IF(isthep(ip).EQ.1)isthep(ip)=2
1285C
1286C and daughter pointers of mother entry must be updated
1287 IF(jdahep(1,ip).EQ.0)THEN
1288 jdahep(1,ip)=ihep
1289 jdahep(2,ip)=ihep
1290 ELSE
1291 jdahep(2,ip)=max(ihep,jdahep(2,ip))
1292 END IF
1293 END IF
1294 END DO
1295C
1296 RETURN
1297 END
1298
1299
1300 FUNCTION ihepdim(DUM)
1301C this is the hepevt class in old style. No d_h_ class pre-name
1302C this is the hepevt class in old style. No d_h_ class pre-name
1303 INTEGER NMXHEP
1304 parameter(nmxhep=4000)
1305 real*8 phep, vhep ! to be real*4/ *8 depending on host
1306 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1307 $ jdahep
1308 COMMON /hepevt/
1309 $ nevhep, ! serial number
1310 $ nhep, ! number of particles
1311 $ isthep(nmxhep), ! status code
1312 $ idhep(nmxhep), ! particle ident KF
1313 $ jmohep(2,nmxhep), ! parent particles
1314 $ jdahep(2,nmxhep), ! childreen particles
1315 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1316 $ vhep(4,nmxhep) ! vertex [mm]
1317* ----------------------------------------------------------------------
1318 LOGICAL qedrad
1319 COMMON /phoqed/
1320 $ qedrad(nmxhep) ! Photos flag
1321* ----------------------------------------------------------------------
1322 SAVE hepevt,phoqed
1323
1324
1325 ihepdim=nhep
1326 END
1327 FUNCTION zprop2(S)
1328 IMPLICIT REAL*8(a-h,o-z)
1329 COMPLEX*16 CPRZ0,CPRZ0M
1330 amz=91.1882
1331 gammz=2.49
1332 cprz0=dcmplx((s-amz**2),s/amz*gammz)
1333 cprz0m=1/cprz0
1334 zprop2=(abs(cprz0m))**2
1335 END
1336
1337 SUBROUTINE taupi0(MODE,JAK,ION)
1338C no initialization required. Must be called once after every:
1339C 1) CALL DEKAY(1+10,...)
1340C 2) CALL DEKAY(2+10,...)
1341C 3) CALL DEXAY(1,...)
1342C 4) CALL DEXAY(2,...)
1343C subroutine to decay originating from TAUOLA's taus:
1344C 1) etas (with CALL TAUETA(JAK))
1345C 2) later pi0's from taus.
1346C 3) extensions to other applications possible.
1347C this routine belongs to >tauola universal interface<, but uses
1348C routines from >tauola< utilities as well. 25.08.2005
1349C this is the hepevt class in old style. No d_h_ class pre-name
1350 INTEGER NMXHEP
1351 parameter(nmxhep=4000)
1352 real*8 phep, vhep ! to be real*4/ *8 depending on host
1353 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1354 $ jdahep
1355 COMMON /hepevt/
1356 $ nevhep, ! serial number
1357 $ nhep, ! number of particles
1358 $ isthep(nmxhep), ! status code
1359 $ idhep(nmxhep), ! particle ident KF
1360 $ jmohep(2,nmxhep), ! parent particles
1361 $ jdahep(2,nmxhep), ! childreen particles
1362 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1363 $ vhep(4,nmxhep) ! vertex [mm]
1364* ----------------------------------------------------------------------
1365 LOGICAL qedrad
1366 COMMON /phoqed/
1367 $ qedrad(nmxhep) ! Photos flag
1368* ----------------------------------------------------------------------
1369 SAVE hepevt,phoqed
1370
1371
1372
1373C position of taus, must be defined by host program:
1374 COMMON /taupos/ np1,np2
1375c
1376 REAL PHOT1(4),PHOT2(4)
1377 real*8 r,x(4),y(4),pi0(4)
1378 INTEGER JEZELI(3),ION(3)
1379 DATA jezeli /0,0,0/
1380 SAVE jezeli
1381 IF (mode.EQ.-1) THEN
1382 jezeli(1)=ion(1)
1383 jezeli(2)=ion(2)
1384 jezeli(3)=ion(3)
1385 RETURN
1386 ENDIF
1387 IF (jezeli(1).EQ.0) RETURN
1388 IF (jezeli(2).EQ.1) CALL taueta(jak)
1389 IF (jezeli(3).EQ.1) CALL tauk0s(jak)
1390C position of decaying particle:
1391 IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
1392 nps=np1
1393 ELSE
1394 nps=np2
1395 ENDIF
1396 nhepm=nhep ! to avoid infinite loop
1397 DO k=jdahep(1,nps),nhepm ! we search for pi0's from tau till eor.
1398 IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k) THEN ! IF we found pi0
1399 DO l=1,4
1400 pi0(l)= phep(l,k)
1401 ENDDO
1402! random 3 vector on the sphere, masless
1403 r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1404 CALL spherd(r,x)
1405 x(4)=r
1406 y(4)=r
1407
1408 y(1)=-x(1)
1409 y(2)=-x(2)
1410 y(3)=-x(3)
1411! boost to lab and to real*4
1412 CALL bostdq(-1,pi0,x,x)
1413 CALL bostdq(-1,pi0,y,y)
1414 DO l=1,4
1415 phot1(l)=x(l)
1416 phot2(l)=y(l)
1417 ENDDO
1418C to hepevt
1419 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1420 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1421 ENDIF
1422 ENDDO
1423C
1424 END
1425 SUBROUTINE taueta(JAK)
1426C subroutine to decay etas's from taus.
1427C this routine belongs to tauola universal interface, but uses
1428C routines from tauola utilities. Just flat phase space, but 4 channels.
1429C it is called at the beginning of SUBR. TAUPI0(JAK)
1430C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1431C this is the hepevt class in old style. No d_h_ class pre-name
1432 INTEGER NMXHEP
1433 parameter(nmxhep=4000)
1434 real*8 phep, vhep ! to be real*4/ *8 depending on host
1435 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1436 $ jdahep
1437 COMMON /hepevt/
1438 $ nevhep, ! serial number
1439 $ nhep, ! number of particles
1440 $ isthep(nmxhep), ! status code
1441 $ idhep(nmxhep), ! particle ident KF
1442 $ jmohep(2,nmxhep), ! parent particles
1443 $ jdahep(2,nmxhep), ! childreen particles
1444 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1445 $ vhep(4,nmxhep) ! vertex [mm]
1446* ----------------------------------------------------------------------
1447 LOGICAL qedrad
1448 COMMON /phoqed/
1449 $ qedrad(nmxhep) ! Photos flag
1450* ----------------------------------------------------------------------
1451 SAVE hepevt,phoqed
1452
1453
1454 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1455 * ,ampiz,ampi,amro,gamro,ama1,gama1
1456 * ,amk,amkz,amkst,gamkst
1457*
1458 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1459 * ,ampiz,ampi,amro,gamro,ama1,gama1
1460 * ,amk,amkz,amkst,gamkst
1461
1462C position of taus, must be defined by host program:
1463 COMMON /taupos/ np1,np2
1464c
1465 REAL RRR(1),BRSUM(3), RR(2)
1466 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1467 real*8 x(4), y(4), z(4)
1468 REAL YM1,YM2,YM3
1469 real*8 r,ru,peta(4),xm1,xm2,xm3,xm,xlam
1470 real*8 a,b,c
1471 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1472C position of decaying particle:
1473 IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
1474 nps=np1
1475 ELSE
1476 nps=np2
1477 ENDIF
1478 nhepm=nhep ! to avoid infinite loop
1479 DO k=jdahep(1,nps),nhepm ! we search for etas's from tau till eor.
1480 IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k) THEN ! IF we found eta
1481 DO l=1,4
1482 peta(l)= phep(l,k) ! eta 4 momentum
1483 ENDDO
1484C eta cumulated branching ratios:
1485 brsum(1)=0.389 ! gamma gamma
1486 brsum(2)=brsum(1)+0.319 ! 3 pi0
1487 brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1488 CALL ranmar(rrr,1)
1489
1490 IF (rrr(1).LT.brsum(1)) THEN ! gamma gamma channel exactly like pi0
1491! random 3 vector on the sphere, masless
1492 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1493 CALL spherd(r,x)
1494 x(4)=r
1495 y(4)=r
1496
1497 y(1)=-x(1)
1498 y(2)=-x(2)
1499 y(3)=-x(3)
1500! boost to lab and to real*4
1501 CALL bostdq(-1,peta,x,x)
1502 CALL bostdq(-1,peta,y,y)
1503 DO l=1,4
1504 phot1(l)=x(l)
1505 phot2(l)=y(l)
1506 ENDDO
1507C to hepevt
1508 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1509 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1510 ELSE ! 3 body channels
1511 IF(rrr(1).LT.brsum(2)) THEN ! 3 pi0
1512 id1= 111
1513 id2= 111
1514 id3= 111
1515 xm1=ampiz ! masses
1516 xm2=ampiz
1517 xm3=ampiz
1518 ELSEIF(rrr(1).LT.brsum(3)) THEN ! pi+ pi- pi0
1519 id1= 211
1520 id2=-211
1521 id3= 111
1522 xm1=ampi ! masses
1523 xm2=ampi
1524 xm3=ampiz
1525 ELSE ! pi+ pi- gamma
1526 id1= 211
1527 id2=-211
1528 id3= 22
1529 xm1=ampi ! masses
1530 xm2=ampi
1531 xm3=0.0
1532 ENDIF
1533 7 CONTINUE ! we generate mass of the first pair:
1534 CALL ranmar(rr,2)
1535 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1536 amin=xm1+xm2
1537 amax=r-xm3
1538 am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1539C weight for flat phase space
1540 wt=xlam(r**2,am2**2,xm3**2)*xlam(am2**2,xm1**2,xm2**2)
1541 & /r**2 /am2**2
1542 IF (rr(2).GT.wt) GOTO 7
1543
1544 ru=xlam(am2**2,xm1**2,xm2**2)/am2/2 ! momenta of the
1545 ! first two products
1546 ! in the rest frame of that pair
1547 CALL spherd(ru,x)
1548 x(4)=sqrt(ru**2+xm1**2)
1549 y(4)=sqrt(ru**2+xm2**2)
1550
1551 y(1)=-x(1)
1552 y(2)=-x(2)
1553 y(3)=-x(3)
1554C generate momentum of that pair in rest frame of eta:
1555 ru=xlam(r**2,am2**2,xm3**2)/r/2
1556 CALL spherd(ru,z)
1557 z(4)=sqrt(ru**2+am2**2)
1558C and boost first two decay products to rest frame of eta.
1559 CALL bostdq(-1,z,x,x)
1560 CALL bostdq(-1,z,y,y)
1561C redefine Z(4) to 4-momentum of the last decay product:
1562 z(1)=-z(1)
1563 z(2)=-z(2)
1564 z(3)=-z(3)
1565 z(4)=sqrt(ru**2+xm3**2)
1566C boost all to lab and move to real*4; also masses
1567 CALL bostdq(-1,peta,x,x)
1568 CALL bostdq(-1,peta,y,y)
1569 CALL bostdq(-1,peta,z,z)
1570 DO l=1,4
1571 phot1(l)=x(l)
1572 phot2(l)=y(l)
1573 phot3(l)=z(l)
1574 ENDDO
1575 ym1=xm1
1576 ym2=xm2
1577 ym3=xm3
1578C to hepevt
1579 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1580 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1581 CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1582 ENDIF
1583
1584 ENDIF
1585 ENDDO
1586C
1587 END
1588 SUBROUTINE tauk0s(JAK)
1589C subroutine to decay K0S's from taus.
1590C this routine belongs to tauola universal interface, but uses
1591C routines from tauola utilities. Just flat phase space, but 4 channels.
1592C it is called at the beginning of SUBR. TAUPI0(JAK)
1593C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1594C this is the hepevt class in old style. No d_h_ class pre-name
1595 INTEGER NMXHEP
1596 parameter(nmxhep=4000)
1597 real*8 phep, vhep ! to be real*4/ *8 depending on host
1598 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1599 $ jdahep
1600 COMMON /hepevt/
1601 $ nevhep, ! serial number
1602 $ nhep, ! number of particles
1603 $ isthep(nmxhep), ! status code
1604 $ idhep(nmxhep), ! particle ident KF
1605 $ jmohep(2,nmxhep), ! parent particles
1606 $ jdahep(2,nmxhep), ! childreen particles
1607 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1608 $ vhep(4,nmxhep) ! vertex [mm]
1609* ----------------------------------------------------------------------
1610 LOGICAL qedrad
1611 COMMON /phoqed/
1612 $ qedrad(nmxhep) ! Photos flag
1613* ----------------------------------------------------------------------
1614 SAVE hepevt,phoqed
1615
1616
1617
1618 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1619 * ,ampiz,ampi,amro,gamro,ama1,gama1
1620 * ,amk,amkz,amkst,gamkst
1621*
1622 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1623 * ,ampiz,ampi,amro,gamro,ama1,gama1
1624 * ,amk,amkz,amkst,gamkst
1625
1626C position of taus, must be defined by host program:
1627 COMMON /taupos/ np1,np2
1628c
1629 REAL RRR(1),BRSUM(3), RR(2)
1630 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1631 real*8 x(4), y(4), z(4)
1632 REAL YM1,YM2,YM3
1633 real*8 r,ru,peta(4),xm1,xm2,xm3,xm,xlam
1634 real*8 a,b,c
1635 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1636C position of decaying particle:
1637 IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
1638 nps=np1
1639 ELSE
1640 nps=np2
1641 ENDIF
1642 nhepm=nhep ! to avoid infinite loop
1643 DO k=jdahep(1,nps),nhepm ! we search for K0S's from tau till eor.
1644 IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k) THEN ! IF we found K0S
1645
1646
1647 DO l=1,4
1648 peta(l)= phep(l,k) ! K0S 4 momentum (this is cloned from eta decay)
1649 ENDDO
1650C K0S cumulated branching ratios:
1651 brsum(1)=0.313 ! 2 PI0
1652 brsum(2)=1.0 ! BRSUM(1)+0.319 ! Pi+ PI-
1653 brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1654 CALL ranmar(rrr,1)
1655
1656 IF(rrr(1).LT.brsum(1)) THEN ! 2 pi0
1657 id1= 111
1658 id2= 111
1659 xm1=ampiz ! masses
1660 xm2=ampiz
1661 ELSEIF(rrr(1).LT.brsum(2)) THEN ! pi+ pi-
1662 id1= 211
1663 id2=-211
1664 xm1=ampi ! masses
1665 xm2=ampi
1666 ELSE ! gamma gamma unused !!!
1667 id1= 22
1668 id2= 22
1669 xm1= 0.0 ! masses
1670 xm2= 0.0
1671 ENDIF
1672
1673! random 3 vector on the sphere, of equal mass !!
1674 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1675 r4=r
1676 r=sqrt(abs(r**2-xm1**2))
1677 CALL spherd(r,x)
1678 x(4)=r4
1679 y(4)=r4
1680
1681 y(1)=-x(1)
1682 y(2)=-x(2)
1683 y(3)=-x(3)
1684! boost to lab and to real*4
1685 CALL bostdq(-1,peta,x,x)
1686 CALL bostdq(-1,peta,y,y)
1687 DO l=1,4
1688 phot1(l)=x(l)
1689 phot2(l)=y(l)
1690 ENDDO
1691
1692 ym1=xm1
1693 ym2=xm2
1694C to hepevt
1695 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1696 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1697
1698C
1699 ENDIF
1700 ENDDO
1701
1702 END