C++ Interface to Tauola
tauola-modified/tauola.f
1
2
3
4
5
6 SUBROUTINE jaker(JAK)
7C *********************
8C
9C **********************************************************************
10C *
11
12
13
14
15C *********TAUOLA LIBRARY: VERSION 2.9 ******** *
16C **************October 2011****************** *
17
18C ** AUTHORS: S.JADACH, Z.WAS ***** *
19C ** R. DECKER, M. JEZABEK, J.H.KUEHN, ***** *
20C ********AVAILABLE FROM: www.cern.ch/wasm **** *
21C *******PUBLISHED IN COMP. PHYS. COMM.******** *
22C *** 76 (1993) 361 **** *
23C *** 64 (1990) 275 **** *
24C *** 70 (1992) 69 **** *
25C *** CLEO initialization: **** *
26C *** Alain Weinstein www home page: **** *
27C *** http://www.cithep.caltech.edu/~ajw/ **** *
28C *** RChL initialization: **** *
29C *** O. Shekhovtsova, T. Przedzinski, **** *
30C *** P. Roig and Z. Was **** *
31C *** IFJPAN-2013-5, UAB-FT-731 **** *
32C **********************************************************************
33C
34C ----------------------------------------------------------------------
35c SUBROUTINE JAKER,
36C CHOOSES DECAY MODE ACCORDING TO LIST OF BRANCHING RATIOS
37C JAK=1 ELECTRON MODE
38C JAK=2 MUON MODE
39C JAK=3 PION MODE
40C JAK=4 RHO MODE
41C JAK=5 A1 MODE
42C JAK=6 K MODE
43C JAK=7 K* MODE
44
45
46
47
48
49C JAK=8 nPI MODE
50
51C
52C called by : DEXAY
53C ----------------------------------------------------------------------
54 COMMON / taubra / gamprt(30),jlist(30),nchan
55
56
57C REAL CUMUL(20)
58
59 REAL CUMUL(30),RRR(1)
60C
61 IF(nchan.LE.0.OR.nchan.GT.30) GOTO 902
62 CALL ranmar(rrr,1)
63 sum=0
64 DO 20 i=1,nchan
65 sum=sum+gamprt(i)
66 20 cumul(i)=sum
67 DO 25 i=nchan,1,-1
68 IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
69 25 CONTINUE
70 jak=jlist(ji)
71 RETURN
72 902 print 9020
73 9020 FORMAT(' ----- JAKER: WRONG NCHAN')
74 stop
75 END
76 SUBROUTINE dekay(KTO,HX)
77C ***********************
78C THIS DEKAY IS IN SPIRIT OF THE 'DECAY' WHICH
79C WAS INCLUDED IN KORAL-B PROGRAM, COMP. PHYS. COMMUN.
80C VOL. 36 (1985) 191, SEE COMMENTS ON GENERAL PHILOSOPHY THERE.
81C KTO=0 INITIALISATION (OBLIGATORY)
82C KTO=1,11 DENOTES TAU+ AND KTO=2,12 TAU-
83C DEKAY(1,H) AND DEKAY(2,H) IS CALLED INTERNALLY BY MC GENERATOR.
84C H DENOTES THE POLARIMETRIC VECTOR, USED BY THE HOST PROGRAM FOR
85C CALCULATION OF THE SPIN WEIGHT.
86C USER MAY OPTIONALLY CALL DEKAY(11,H) DEKAY(12,H) IN ORDER
87C TO TRANSFORM DECAY PRODUCTS TO CMS AND WRITE LUND RECORD IN /LUJETS/.
88C KTO=100, PRINT FINAL REPORT (OPTIONAL).
89C DECAY MODES:
90C JAK=1 ELECTRON DECAY
91C JAK=2 MU DECAY
92C JAK=3 PI DECAY
93C JAK=4 RHO DECAY
94C JAK=5 A1 DECAY
95C JAK=6 K DECAY
96C JAK=7 K* DECAY
97
98
99
100
101
102
103C JAK=8 NPI DECAY
104C JAK=0 INCLUSIVE: JAK=1,2,3,4,5,6,7,8
105
106 REAL H(4)
107 real*8 hx(4)
108 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
109
110
111
112 COMMON / idfc / idf
113
114 COMMON /taupos/ np1,np2
115 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
116 real*4 gampmc ,gamper
117 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
118
119
120
121 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
122
123 & ,names
124 CHARACTER NAMES(NMODE)*31
125 COMMON / inout / inut,iout
126 COMMON /ipcht/ iver
127 INTEGER IVER
128
129 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4),HDUM(4)
130 REAL PDUMX(4,9)
131 DATA iwarm/0/
132 ktom=kto
133
134
135
136 IF(kto.EQ.-1) THEN
137C ==================
138C INITIALISATION OR REINITIALISATION
139C first or second tau positions in HEPEVT as in KORALB/Z
140 np1=3
141 np2=4
142 ktom=1
143 IF (iwarm.EQ.1) x=5/(iwarm-1)
144 iwarm=1
145 WRITE(iout,7001) jak1,jak
146 WRITE(iout,7002) iver
147 nevtot=0
148 nev1=0
149 nev2=0
150 IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
151 CALL dadmel(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
152 CALL dadmmu(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
153 CALL dadmpi(-1,idum,hdum,pdum1,pdum2)
154 CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
155 CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
156 CALL dadmkk(-1,idum,hdum,pdum1,pdum2)
157 CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
158 CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
159 ENDIF
160 DO 21 i=1,30
161 nevdec(i)=0
162 gampmc(i)=0
163 21 gamper(i)=0
164 ELSEIF(kto.EQ.1) THEN
165C =====================
166C DECAY OF TAU+ IN THE TAU REST FRAME
167 nevtot=nevtot+1
168 IF(iwarm.EQ.0) GOTO 902
169 isgn= idf/iabs(idf)
170
171
172
173C AJWMOD to change BRs depending on sign:
174 CALL taurdf(kto)
175
176 CALL dekay1(0,h,isgn)
177 ELSEIF(kto.EQ.2) THEN
178C =================================
179C DECAY OF TAU- IN THE TAU REST FRAME
180 nevtot=nevtot+1
181 IF(iwarm.EQ.0) GOTO 902
182 isgn=-idf/iabs(idf)
183
184
185
186C AJWMOD to change BRs depending on sign:
187 CALL taurdf(kto)
188
189 CALL dekay2(0,h,isgn)
190 ELSEIF(kto.EQ.11) THEN
191C ======================
192C REST OF DECAY PROCEDURE FOR ACCEPTED TAU+ DECAY
193 nev1=nev1+1
194 isgn= idf/iabs(idf)
195 CALL dekay1(1,h,isgn)
196 ELSEIF(kto.EQ.12) THEN
197C ======================
198C REST OF DECAY PROCEDURE FOR ACCEPTED TAU- DECAY
199 nev2=nev2+1
200 isgn=-idf/iabs(idf)
201 CALL dekay2(1,h,isgn)
202 ELSEIF(kto.EQ.100) THEN
203C =======================
204 IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
205 CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
206 CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
207 CALL dadmpi( 1,idum,hdum,pdum1,pdum2)
208 CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
209 CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
210 CALL dadmkk( 1,idum,hdum,pdum1,pdum2)
211 CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
212 CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
213 WRITE(iout,7010) nev1,nev2,nevtot
214 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
215 WRITE(iout,7012)
216 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
217 WRITE(iout,7013)
218 ENDIF
219 ELSE
220C ====
221 GOTO 910
222 ENDIF
223C =====
224 DO 78 k=1,4
225 78 hx(k)=h(k)
226 RETURN
227 7001 FORMAT(///1x,15(5h*****)
228 $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
229 $ /,' *', 25x,'***********October 2011 ***************',9x,1h*,
230 $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
231 $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
232 $ /,' *', 25x,'**AVAILABLE FROM: www.cern.ch/wasm**** ',9x,1h*,
233 $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
234 $ /,' *', 25x,'0: Physics initialization CLEO collab ',9x,1h*,
235 $ /,' *', 25x,' see Alain Weinstein www home page: ',9x,1h*,
236 $ /,' *', 25x,'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
237 $ /,' *', 25x,'/korb_doc.html#files ',9x,1h*,
238 $ /,' *', 25x,'1: Physics initialization RChL of: ',9x,1h*,
239 $ /,' *', 25x,' O. Shekhovtsova, T. Przedzinski, ',9x,1h*,
240 $ /,' *', 25x,' P. Roig and Z. Was ',9x,1h*,
241 $ /,' *', 25x,' IFJPAN-2013-5, UAB-FT-731 ',9x,1h*,
242 $ /,' *', 25x,'*******CPC 76 (1993) 361 *****',9x,1h*,
243 $ /,' *', 25x,'**5 or more pi dec.: precision limited ',9x,1h*,
244 $ /,' *', 25x,'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
245 $ /,' *',i20 ,5x,'JAK1 = DECAY MODE TAU+ ',9x,1h*,
246 $ /,' *',i20 ,5x,'JAK2 = DECAY MODE TAU- ',9x,1h*,
247 $ /,1x,15(5h*****)/)
248 7002 FORMAT(' *',i20 ,5x,'IVER = hadronic current version ',9x,1h*)
249 7010 FORMAT(///1x,15(5h*****)
250 $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
251 $ /,' *', 25x,'***********October 2011 ***************',9x,1h*,
252 $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
253 $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
254 $ /,' *', 25x,'**AVAILABLE FROM: www.cern.ch/wasm ****',9x,1h*,
255 $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
256 $ /,' *', 25x,'******* 64 (1990) 275 *****',9x,1h*,
257 $ /,' *', 25x,'******* 70 (1992) 69 *****',9x,1h*,
258 $ /,' *', 25x,'******* 76 (1993) 361 *****',9x,1h*,
259 $ /,' *', 25x,'******* IFJPAN-2013-5, UAB-FT-731 **',9x,1h*,
260 $ /,' *', 25x,'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
261 $ /,' *',i20 ,5x,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
262 $ /,' *',i20 ,5x,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
263 $ /,' *',i20 ,5x,'NEVTOT = SUM ',9x,1h*,
264 $ /,' *',' NOEVTS ',
265 $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
266 7011 FORMAT(1x,'*'
267 $ ,i10,2f12.7 ,' DADMEL ELECTRON ',9x,1h*
268 $ /,' *',i10,2f12.7 ,' DADMMU MUON ',9x,1h*
269 $ /,' *',i10,2f12.7 ,' DADMPI PION ',9x,1h*
270 $ /,' *',i10,2f12.7, ' DADMRO RHO (->2PI) ',9x,1h*
271 $ /,' *',i10,2f12.7, ' DADMAA A1 (->3PI) ',9x,1h*
272 $ /,' *',i10,2f12.7, ' DADMKK KAON ',9x,1h*
273 $ /,' *',i10,2f12.7, ' DADMKS K* ',9x,1h*)
274 7012 FORMAT(1x,'*'
275 $ ,i10,2f12.7,a31 ,8x,1h*)
276 7013 FORMAT(1x,'*'
277 $ ,20x,'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
278 $ /,' *',20x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
279 $ /,1x,15(5h*****)/)
280 902 print 9020
281 9020 FORMAT(' ----- DEKAY: LACK OF INITIALISATION')
282 stop
283 910 print 9100
284 9100 FORMAT(' ----- DEKAY: WRONG VALUE OF KTO ')
285 stop
286 END
287 SUBROUTINE dekay1(IMOD,HH,ISGN)
288C *******************************
289C THIS ROUTINE SIMULATES TAU+ DECAY
290 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
291 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
292 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
293 real*4 gampmc ,gamper
294 REAL HH(4)
295 REAL HV(4),PNU(4),PPI(4)
296 REAL PWB(4),PMU(4),PNM(4)
297 REAL PRHO(4),PIC(4),PIZ(4)
298 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
299 REAL PKK(4),PKS(4)
300 REAL PNPI(4,9)
301 REAL PHOT(4)
302 REAL PDUM(4)
303 DATA nev,nprin/0,10/
304 kto=1
305 IF(jak1.EQ.-1) RETURN
306 imd=imod
307 IF(imd.EQ.0) THEN
308C =================
309 jak=jak1
310 IF(jak1.EQ.0) CALL jaker(jak)
311 IF(jak.EQ.1) THEN
312 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
313 ELSEIF(jak.EQ.2) THEN
314 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
315 ELSEIF(jak.EQ.3) THEN
316 CALL dadmpi(0, isgn,hv,ppi,pnu)
317 ELSEIF(jak.EQ.4) THEN
318 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
319 ELSEIF(jak.EQ.5) THEN
320 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
321 ELSEIF(jak.EQ.6) THEN
322 CALL dadmkk(0, isgn,hv,pkk,pnu)
323 ELSEIF(jak.EQ.7) THEN
324 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
325 ELSE
326 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
327 ENDIF
328 DO 33 i=1,3
329 33 hh(i)=hv(i)
330 hh(4)=1.0
331
332 ELSEIF(imd.EQ.1) THEN
333C =====================
334 nev=nev+1
335 IF (jak.LT.31) THEN
336 nevdec(jak)=nevdec(jak)+1
337 ENDIF
338 DO 34 i=1,4
339 34 pdum(i)=.0
340 IF(jak.EQ.1) THEN
341 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
342 CALL dwrph(ktom,phot)
343 DO 10 i=1,4
344 10 pp1(i)=pmu(i)
345
346 ELSEIF(jak.EQ.2) THEN
347 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
348 CALL dwrph(ktom,phot)
349 DO 20 i=1,4
350 20 pp1(i)=pmu(i)
351
352 ELSEIF(jak.EQ.3) THEN
353 CALL dwlupi(1,isgn,ppi,pnu)
354 DO 30 i=1,4
355 30 pp1(i)=ppi(i)
356
357 ELSEIF(jak.EQ.4) THEN
358 CALL dwluro(1,isgn,pnu,prho,pic,piz)
359 DO 40 i=1,4
360 40 pp1(i)=prho(i)
361
362 ELSEIF(jak.EQ.5) THEN
363 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
364 DO 50 i=1,4
365 50 pp1(i)=paa(i)
366 ELSEIF(jak.EQ.6) THEN
367 CALL dwlukk(1,isgn,pkk,pnu)
368 DO 60 i=1,4
369 60 pp1(i)=pkk(i)
370 ELSEIF(jak.EQ.7) THEN
371 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
372 DO 70 i=1,4
373 70 pp1(i)=pks(i)
374 ELSE
375CAM MULTIPION DECAY
376 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
377 DO 80 i=1,4
378 80 pp1(i)=pwb(i)
379 ENDIF
380
381 ENDIF
382C =====
383 END
384 SUBROUTINE dekay2(IMOD,HH,ISGN)
385C *******************************
386C THIS ROUTINE SIMULATES TAU- DECAY
387 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
388 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
389 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
390 real*4 gampmc ,gamper
391 REAL HH(4)
392 REAL HV(4),PNU(4),PPI(4)
393 REAL PWB(4),PMU(4),PNM(4)
394 REAL PRHO(4),PIC(4),PIZ(4)
395 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
396 REAL PKK(4),PKS(4)
397 REAL PNPI(4,9)
398 REAL PHOT(4)
399 REAL PDUM(4)
400 DATA nev,nprin/0,10/
401 kto=2
402 IF(jak2.EQ.-1) RETURN
403 imd=imod
404 IF(imd.EQ.0) THEN
405C =================
406 jak=jak2
407 IF(jak2.EQ.0) CALL jaker(jak)
408 IF(jak.EQ.1) THEN
409 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
410 ELSEIF(jak.EQ.2) THEN
411 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
412 ELSEIF(jak.EQ.3) THEN
413 CALL dadmpi(0, isgn,hv,ppi,pnu)
414 ELSEIF(jak.EQ.4) THEN
415 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
416 ELSEIF(jak.EQ.5) THEN
417 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
418 ELSEIF(jak.EQ.6) THEN
419 CALL dadmkk(0, isgn,hv,pkk,pnu)
420 ELSEIF(jak.EQ.7) THEN
421 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
422 ELSE
423 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
424 ENDIF
425 DO 33 i=1,3
426 33 hh(i)=hv(i)
427 hh(4)=1.0
428 ELSEIF(imd.EQ.1) THEN
429C =====================
430 nev=nev+1
431 IF (jak.LT.31) THEN
432 nevdec(jak)=nevdec(jak)+1
433 ENDIF
434 DO 34 i=1,4
435 34 pdum(i)=.0
436 IF(jak.EQ.1) THEN
437 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
438 CALL dwrph(ktom,phot)
439 DO 10 i=1,4
440 10 pp2(i)=pmu(i)
441
442 ELSEIF(jak.EQ.2) THEN
443 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
444 CALL dwrph(ktom,phot)
445 DO 20 i=1,4
446 20 pp2(i)=pmu(i)
447
448 ELSEIF(jak.EQ.3) THEN
449 CALL dwlupi(2,isgn,ppi,pnu)
450 DO 30 i=1,4
451 30 pp2(i)=ppi(i)
452
453 ELSEIF(jak.EQ.4) THEN
454 CALL dwluro(2,isgn,pnu,prho,pic,piz)
455 DO 40 i=1,4
456 40 pp2(i)=prho(i)
457
458 ELSEIF(jak.EQ.5) THEN
459 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
460 DO 50 i=1,4
461 50 pp2(i)=paa(i)
462 ELSEIF(jak.EQ.6) THEN
463 CALL dwlukk(2,isgn,pkk,pnu)
464 DO 60 i=1,4
465 60 pp1(i)=pkk(i)
466 ELSEIF(jak.EQ.7) THEN
467 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
468 DO 70 i=1,4
469 70 pp1(i)=pks(i)
470 ELSE
471CAM MULTIPION DECAY
472 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
473 DO 80 i=1,4
474 80 pp1(i)=pwb(i)
475 ENDIF
476C
477 ENDIF
478C =====
479 END
480 SUBROUTINE dexay(KTO,POL)
481C ----------------------------------------------------------------------
482C THIS 'DEXAY' IS A ROUTINE WHICH GENERATES DECAY OF THE SINGLE
483C POLARIZED TAU, POL IS A POLARIZATION VECTOR (NOT A POLARIMETER
484C VECTOR AS IN DEKAY) OF THE TAU AND IT IS AN INPUT PARAMETER.
485C KTO=0 INITIALISATION (OBLIGATORY)
486C KTO=1 DENOTES TAU+ AND KTO=2 TAU-
487C DEXAY(1,POL) AND DEXAY(2,POL) ARE CALLED INTERNALLY BY MC GENERATOR.
488C DECAY PRODUCTS ARE TRANSFORMED READILY
489C TO CMS AND WRITEN IN THE LUND RECORD IN /LUJETS/
490C KTO=100, PRINT FINAL REPORT (OPTIONAL).
491C
492C called by : KORALZ
493C ----------------------------------------------------------------------
494 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
495 real*4 gampmc ,gamper
496 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
497 COMMON / idfc / idff
498 COMMON /taupos/ np1,np2
499 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
500 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
501 & ,names
502 CHARACTER NAMES(NMODE)*31
503 COMMON / inout / inut,iout
504 COMMON /ipcht/ iver
505 INTEGER IVER
506
507 REAL POL(4)
508 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
509 REAL PDUM(4)
510 REAL PDUMI(4,9)
511 DATA iwarm/0/
512 ktom=kto
513C
514 IF(kto.EQ.-1) THEN
515C ==================
516
517C INITIALISATION OR REINITIALISATION
518C first or second tau positions in HEPEVT as in KORALB/Z
519 np1=3
520 np2=4
521 iwarm=1
522 WRITE(iout, 7001) jak1,jak2
523 WRITE(iout,7002) iver
524 nevtot=0
525 nev1=0
526 nev2=0
527 IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
528 CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
529 CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
530 CALL dexpi(-1,idum,pdum,pdum1,pdum2)
531 CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
532 CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
533 CALL dexkk(-1,idum,pdum,pdum1,pdum2)
534 CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
535 CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
536 ENDIF
537 DO 21 i=1,30
538 nevdec(i)=0
539 gampmc(i)=0
540 21 gamper(i)=0
541 ELSEIF(kto.EQ.1) THEN
542C =====================
543C DECAY OF TAU+ IN THE TAU REST FRAME
544 nevtot=nevtot+1
545 nev1=nev1+1
546 IF(iwarm.EQ.0) GOTO 902
547 isgn=idff/iabs(idff)
548CAM CALL DEXAY1(POL,ISGN)
549 CALL dexay1(kto,jak1,jakp,pol,isgn)
550 ELSEIF(kto.EQ.2) THEN
551C =================================
552C DECAY OF TAU- IN THE TAU REST FRAME
553 nevtot=nevtot+1
554 nev2=nev2+1
555 IF(iwarm.EQ.0) GOTO 902
556 isgn=-idff/iabs(idff)
557CAM CALL DEXAY2(POL,ISGN)
558 CALL dexay1(kto,jak2,jakm,pol,isgn)
559 ELSEIF(kto.EQ.100) THEN
560C =======================
561 IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
562 CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
563 CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
564 CALL dexpi( 1,idum,pdum,pdum1,pdum2)
565 CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
566 CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
567 CALL dexkk( 1,idum,pdum,pdum1,pdum2)
568 CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
569 CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
570 WRITE(iout,7010) nev1,nev2,nevtot
571 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
572 WRITE(iout,7012)
573 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
574 WRITE(iout,7013)
575 ENDIF
576 ELSE
577 GOTO 910
578 ENDIF
579 RETURN
580 7001 FORMAT(///1x,15(5h*****)
581 $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
582 $ /,' *', 25x,'***********October 2011***************',9x,1h*,
583 $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
584 $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
585 $ /,' *', 25x,'**AVAILABLE FROM: www.cern.ch/wasm****',9x,1h*,
586 $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
587 $ /,' *', 25x,'0: Physics initialization CLEO collab ',9x,1h*,
588 $ /,' *', 25x,' see Alain Weinstein www home page: ',9x,1h*,
589 $ /,' *', 25x,'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
590 $ /,' *', 25x,'/korb_doc.html#files ',9x,1h*,
591 $ /,' *', 25x,'1: Physics initialization RChL of: ',9x,1h*,
592 $ /,' *', 25x,' O. Shekhovtsova, T. Przedzinski, ',9x,1h*,
593 $ /,' *', 25x,' P. Roig and Z. Was ',9x,1h*,
594 $ /,' *', 25x,' IFJPAN-2013-5, UAB-FT-731 ',9x,1h*,
595 $ /,' *', 25x,'****** CPC 76 (1993) 361 ******',9x,1h*,
596 $ /,' *', 25x,'**5 or more pi dec.: precision limited ',9x,1h*,
597 $ /,' *', 25x,'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
598 $ /,' *',i20 ,5x,'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
599 $ /,' *',i20 ,5x,'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
600 $ /,1x,15(5h*****)/)
601 7002 FORMAT(' *',i20 ,5x,'IVER = hadronic current version ',9x,1h*)
602CHBU format 7010 had more than 19 continuation lines
603CHBU split into two
604 7010 FORMAT(///1x,15(5h*****)
605 $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
606 $ /,' *', 25x,'***********October 2011***************',9x,1h*,
607 $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
608 $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
609 $ /,' *', 25x,'**AVAILABLE FROM: www.cern.ch/wasm ****',9x,1h*,
610 $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
611 $ /,' *', 25x,'******* 64 (1990) 275 *****',9x,1h*,
612 $ /,' *', 25x,'******* 70 (1992) 69 *****',9x,1h*,
613 $ /,' *', 25x,'******* 76 (1993) 361 *****',9x,1h*,
614 $ /,' *', 25x,'******* IFJPAN-2013-5, UAB-FT-731 **',9x,1h*,
615 $ /,' *', 25x,'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
616 $ /,' *',i20 ,5x,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
617 $ /,' *',i20 ,5x,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
618 $ /,' *',i20 ,5x,'NEVTOT = SUM ',9x,1h*
619 $ /,' *',' NOEVTS ',
620 $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
621 7011 FORMAT(1x,'*'
622 $ ,i10,2f12.7 ,' DADMEL ELECTRON ',9x,1h*
623 $ /,' *',i10,2f12.7 ,' DADMMU MUON ',9x,1h*
624 $ /,' *',i10,2f12.7 ,' DADMPI PION ',9x,1h*
625 $ /,' *',i10,2f12.7, ' DADMRO RHO (->2PI) ',9x,1h*
626 $ /,' *',i10,2f12.7, ' DADMAA A1 (->3PI) ',9x,1h*
627 $ /,' *',i10,2f12.7, ' DADMKK KAON ',9x,1h*
628 $ /,' *',i10,2f12.7, ' DADMKS K* ',9x,1h*)
629 7012 FORMAT(1x,'*'
630 $ ,i10,2f12.7,a31 ,8x,1h*)
631 7013 FORMAT(1x,'*'
632 $ ,20x,'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
633 $ /,' *',20x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
634 $ /,1x,15(5h*****)/)
635 902 WRITE(iout, 9020)
636 9020 FORMAT(' ----- DEXAY: LACK OF INITIALISATION')
637 stop
638 910 WRITE(iout, 9100)
639 9100 FORMAT(' ----- DEXAY: WRONG VALUE OF KTO ')
640 stop
641 END
642 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
643C ---------------------------------------------------------------------
644C THIS ROUTINE SIMULATES TAU+- DECAY
645C
646C called by : DEXAY
647C ---------------------------------------------------------------------
648 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
649 real*4 gampmc ,gamper
650 COMMON / inout / inut,iout
651 REAL POL(4),POLAR(4)
652 REAL PNU(4),PPI(4)
653 REAL PRHO(4),PIC(4),PIZ(4)
654 REAL PWB(4),PMU(4),PNM(4)
655 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
656 REAL PKK(4),PKS(4)
657 REAL PNPI(4,9)
658 REAL PHOT(4)
659 REAL PDUM(4)
660C
661 IF(jakin.EQ.-1) RETURN
662 DO 33 i=1,3
663 33 polar(i)=pol(i)
664 polar(4)=0.
665 DO 34 i=1,4
666 34 pdum(i)=.0
667 jak=jakin
668 IF(jak.EQ.0) CALL jaker(jak)
669CAM
670 IF(jak.EQ.1) THEN
671 CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
672 CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
673 CALL dwrph(kto,phot )
674 ELSEIF(jak.EQ.2) THEN
675 CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
676 CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
677 CALL dwrph(kto,phot )
678 ELSEIF(jak.EQ.3) THEN
679 CALL dexpi(0, isgn,polar,ppi,pnu)
680 CALL dwlupi(kto,isgn,ppi,pnu)
681 ELSEIF(jak.EQ.4) THEN
682 CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
683 CALL dwluro(kto,isgn,pnu,prho,pic,piz)
684 ELSEIF(jak.EQ.5) THEN
685 CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
686 CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
687 ELSEIF(jak.EQ.6) THEN
688 CALL dexkk(0, isgn,polar,pkk,pnu)
689 CALL dwlukk(kto,isgn,pkk,pnu)
690 ELSEIF(jak.EQ.7) THEN
691 CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
692 CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
693 ELSE
694 jnpi=jak-7
695 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
696 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
697 ENDIF
698 nevdec(jak)=nevdec(jak)+1
699 END
700 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
701C ----------------------------------------------------------------------
702C THIS SIMULATES TAU DECAY IN TAU REST FRAME
703C INTO ELECTRON AND TWO NEUTRINOS
704C
705C called by : DEXAY,DEXAY1
706C ----------------------------------------------------------------------
707 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
708 DATA iwarm/0/
709C
710 IF(mode.EQ.-1) THEN
711C ===================
712 iwarm=1
713 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
714CC CALL HBOOK1(813,'WEIGHT DISTRIBUTION DEXEL $',100,0,2)
715C
716 ELSEIF(mode.EQ. 0) THEN
717C =======================
718300 CONTINUE
719 IF(iwarm.EQ.0) GOTO 902
720 CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
721 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
722CC CALL HFILL(813,WT)
723 CALL ranmar(rn,1)
724 IF(rn(1).GT.wt) GOTO 300
725C
726 ELSEIF(mode.EQ. 1) THEN
727C =======================
728 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
729CC CALL HPRINT(813)
730 ENDIF
731C =====
732 RETURN
733 902 print 9020
734 9020 FORMAT(' ----- DEXEL: LACK OF INITIALISATION')
735 stop
736 END
737 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
738C ----------------------------------------------------------------------
739C THIS SIMULATES TAU DECAY IN ITS REST FRAME
740C INTO MUON AND TWO NEUTRINOS
741C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
742C PWB W-BOSON
743C Q1 MUON
744C Q2 MUON-NEUTRINO
745C ----------------------------------------------------------------------
746 COMMON / inout / inut,iout
747 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
748 DATA iwarm/0/
749C
750 IF(mode.EQ.-1) THEN
751C ===================
752 iwarm=1
753 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
754CC CALL HBOOK1(814,'WEIGHT DISTRIBUTION DEXMU $',100,0,2)
755C
756 ELSEIF(mode.EQ. 0) THEN
757C =======================
758300 CONTINUE
759 IF(iwarm.EQ.0) GOTO 902
760 CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
761 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
762CC CALL HFILL(814,WT)
763 CALL ranmar(rn,1)
764 IF(rn(1).GT.wt) GOTO 300
765C
766 ELSEIF(mode.EQ. 1) THEN
767C =======================
768 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
769CC CALL HPRINT(814)
770 ENDIF
771C =====
772 RETURN
773 902 WRITE(iout, 9020)
774 9020 FORMAT(' ----- DEXMU: LACK OF INITIALISATION')
775 stop
776 END
777 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
778C ----------------------------------------------------------------------
779C
780C called by : DEXEL,(DEKAY,DEKAY1)
781C ----------------------------------------------------------------------
782 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
783 real*4 gfermi,gv,ga,ccabib,scabib,gamel
784 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
785 * ,ampiz,ampi,amro,gamro,ama1,gama1
786 * ,amk,amkz,amkst,gamkst
787C
788 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
789 * ,ampiz,ampi,amro,gamro,ama1,gama1
790 * ,amk,amkz,amkst,gamkst
791 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
792 real*4 gampmc ,gamper
793 real*4 phx(4)
794 COMMON / inout / inut,iout
795 REAL HHV(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4)
796 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
797 real*4 rrr(3)
798 real*8 swt, sswt
799 DATA pi /3.141592653589793238462643/
800 DATA iwarm/0/
801C
802 IF(mode.EQ.-1) THEN
803C ===================
804 iwarm=1
805 nevraw=0
806 nevacc=0
807 nevovr=0
808 swt=0
809 sswt=0
810 wtmax=1e-20
811 DO 15 i=1,500
812 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
813 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
81415 CONTINUE
815CC CALL HBOOK1(803,'WEIGHT DISTRIBUTION DADMEL $',100,0,2)
816C
817 ELSEIF(mode.EQ. 0) THEN
818C =======================
819300 CONTINUE
820 IF(iwarm.EQ.0) GOTO 902
821 nevraw=nevraw+1
822 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
823CC CALL HFILL(803,WT/WTMAX)
824 swt=swt+wt
825 sswt=sswt+wt**2
826 CALL ranmar(rrr,3)
827 rn=rrr(1)
828 IF(wt.GT.wtmax) nevovr=nevovr+1
829 IF(rn*wtmax.GT.wt) GOTO 300
830C ROTATIONS TO BASIC TAU REST FRAME
831 rr2=rrr(2)
832 costhe=-1.+2.*rr2
833 thet=acos(costhe)
834 rr3=rrr(3)
835 phi =2*pi*rr3
836 CALL rotor2(thet,pnu,pnu)
837 CALL rotor3( phi,pnu,pnu)
838 CALL rotor2(thet,pwb,pwb)
839 CALL rotor3( phi,pwb,pwb)
840 CALL rotor2(thet,q1,q1)
841 CALL rotor3( phi,q1,q1)
842 CALL rotor2(thet,q2,q2)
843 CALL rotor3( phi,q2,q2)
844 CALL rotor2(thet,hv,hv)
845 CALL rotor3( phi,hv,hv)
846 CALL rotor2(thet,phx,phx)
847 CALL rotor3( phi,phx,phx)
848 DO 44,i=1,3
849 44 hhv(i)=-isgn*hv(i)
850 nevacc=nevacc+1
851C
852 ELSEIF(mode.EQ. 1) THEN
853C =======================
854 IF(nevraw.EQ.0) RETURN
855 pargam=swt/float(nevraw+1)
856 error=0
857 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
858 rat=pargam/gamel
859 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
860CC CALL HPRINT(803)
861 gampmc(1)=rat
862 gamper(1)=error
863CAM NEVDEC(1)=NEVACC
864 ENDIF
865C =====
866 RETURN
867 7010 FORMAT(///1x,15(5h*****)
868 $ /,' *', 25x,'******** DADMEL FINAL REPORT ******** ',9x,1h*
869 $ /,' *',i20 ,5x,'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
870 $ /,' *',i20 ,5x,'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
871 $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
872 $ /,' *',e20.5,5x,'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
873 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
874 $ /,' *',f20.9,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
875 $ /,' *',25x, 'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
876 $ /,' *',25x, 'BUT ONLY V-A CUPLINGS ',9x,1h*
877 $ /,1x,15(5h*****)/)
878 902 WRITE(iout, 9020)
879 9020 FORMAT(' ----- DADMEL: LACK OF INITIALISATION')
880 stop
881 END
882 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
883C ----------------------------------------------------------------------
884 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
885 * ,ampiz,ampi,amro,gamro,ama1,gama1
886 * ,amk,amkz,amkst,gamkst
887C
888 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
889 * ,ampiz,ampi,amro,gamro,ama1,gama1
890 * ,amk,amkz,amkst,gamkst
891 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
892 real*4 gfermi,gv,ga,ccabib,scabib,gamel
893 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
894 real*4 gampmc ,gamper
895 COMMON / inout / inut,iout
896 real*4 phx(4)
897 REAL HHV(4),HV(4),PNU(4),PWB(4),Q1(4),Q2(4)
898 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
899 real*4 rrr(3)
900 real*8 swt, sswt
901 DATA pi /3.141592653589793238462643/
902 DATA iwarm /0/
903C
904 IF(mode.EQ.-1) THEN
905C ===================
906 iwarm=1
907 nevraw=0
908 nevacc=0
909 nevovr=0
910 swt=0
911 sswt=0
912 wtmax=1e-20
913 DO 15 i=1,500
914 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
915 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
91615 CONTINUE
917CC CALL HBOOK1(802,'WEIGHT DISTRIBUTION DADMMU $',100,0,2)
918C
919 ELSEIF(mode.EQ. 0) THEN
920C =======================
921300 CONTINUE
922 IF(iwarm.EQ.0) GOTO 902
923 nevraw=nevraw+1
924 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
925CC CALL HFILL(802,WT/WTMAX)
926 swt=swt+wt
927 sswt=sswt+wt**2
928 CALL ranmar(rrr,3)
929 rn=rrr(1)
930 IF(wt.GT.wtmax) nevovr=nevovr+1
931 IF(rn*wtmax.GT.wt) GOTO 300
932C ROTATIONS TO BASIC TAU REST FRAME
933 costhe=-1.+2.*rrr(2)
934 thet=acos(costhe)
935 phi =2*pi*rrr(3)
936 CALL rotor2(thet,pnu,pnu)
937 CALL rotor3( phi,pnu,pnu)
938 CALL rotor2(thet,pwb,pwb)
939 CALL rotor3( phi,pwb,pwb)
940 CALL rotor2(thet,q1,q1)
941 CALL rotor3( phi,q1,q1)
942 CALL rotor2(thet,q2,q2)
943 CALL rotor3( phi,q2,q2)
944 CALL rotor2(thet,hv,hv)
945 CALL rotor3( phi,hv,hv)
946 CALL rotor2(thet,phx,phx)
947 CALL rotor3( phi,phx,phx)
948 DO 44,i=1,3
949 44 hhv(i)=-isgn*hv(i)
950 nevacc=nevacc+1
951C
952 ELSEIF(mode.EQ. 1) THEN
953C =======================
954 IF(nevraw.EQ.0) RETURN
955 pargam=swt/float(nevraw+1)
956 error=0
957 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
958 rat=pargam/gamel
959 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
960CC CALL HPRINT(802)
961 gampmc(2)=rat
962 gamper(2)=error
963CAM NEVDEC(2)=NEVACC
964 ENDIF
965C =====
966 RETURN
967 7010 FORMAT(///1x,15(5h*****)
968 $ /,' *', 25x,'******** DADMMU FINAL REPORT ******** ',9x,1h*
969 $ /,' *',i20 ,5x,'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
970 $ /,' *',i20 ,5x,'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
971 $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
972 $ /,' *',e20.5,5x,'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
973 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
974 $ /,' *',f20.9,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
975 $ /,' *',25x, 'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
976 $ /,' *',25x, 'BUT ONLY V-A CUPLINGS ',9x,1h*
977 $ /,1x,15(5h*****)/)
978 902 WRITE(iout, 9020)
979 9020 FORMAT(' ----- DADMMU: LACK OF INITIALISATION')
980 stop
981 END
982 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
983C XNX,XNA was flipped in parameters of dphsel and dphsmu
984C *********************************************************************
985C * ELECTRON DECAY MODE *
986C *********************************************************************
987 real*4 phx(4)
988 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
989 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
990 real*8 dgamt
991 ielmu=1
992 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
993 DO 7 k=1,4
994 hvx(k)=hv(k)
995 phx(k)=ph(k)
996 paax(k)=paa(k)
997 xax(k)=xa(k)
998 qpx(k)=qp(k)
999 xnx(k)=xn(k)
1000 7 CONTINUE
1001 dgamx=dgamt
1002 END
1003 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1004C XNX,XNA was flipped in parameters of dphsel and dphsmu
1005C *********************************************************************
1006C * MUON DECAY MODE *
1007C *********************************************************************
1008 real*4 phx(4)
1009 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
1010 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1011 real*8 dgamt
1012 ielmu=2
1013 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1014 DO 7 k=1,4
1015 hvx(k)=hv(k)
1016 phx(k)=ph(k)
1017 paax(k)=paa(k)
1018 xax(k)=xa(k)
1019 qpx(k)=qp(k)
1020 xnx(k)=xn(k)
1021 7 CONTINUE
1022 dgamx=dgamt
1023 END
1024 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
1025 IMPLICIT REAL*8 (a-h,o-z)
1026C ----------------------------------------------------------------------
1027* IT SIMULATES E,MU CHANNELS OF TAU DECAY IN ITS REST FRAME WITH
1028* QED ORDER ALPHA CORRECTIONS
1029C ----------------------------------------------------------------------
1030 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1031 * ,ampiz,ampi,amro,gamro,ama1,gama1
1032 * ,amk,amkz,amkst,gamkst
1033C
1034 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1035 * ,ampiz,ampi,amro,gamro,ama1,gama1
1036 * ,amk,amkz,amkst,gamkst
1037 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1038 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1039 COMMON / inout / inut,iout
1040 COMMON / taurad / xk0dec,itdkrc
1041 real*8 xk0dec
1042 real*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1043 real*8 pr(4)
1044 real*4 rrr(6)
1045 LOGICAL IHARD
1046 DATA pi /3.141592653589793238462643d0/
1047C AJWMOD to satisfy compiler, comment out this unused function.
1048C AMRO, GAMRO IS ONLY A PARAMETER FOR GETING HIGHT EFFICIENCY
1049C
1050C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
1051C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P)
1052 phspac=1./2**17/pi**8
1053 amtax=amtau
1054C TAU MOMENTUM
1055 pt(1)=0.d0
1056 pt(2)=0.d0
1057 pt(3)=0.d0
1058 pt(4)=amtax
1059C
1060 CALL ranmar(rrr,6)
1061C
1062 IF (ielmu.EQ.1) THEN
1063 amu=amel
1064 ELSE
1065 amu=ammu
1066 ENDIF
1067C
1068 prhard=0.30d0
1069 IF ( itdkrc.EQ.0) prhard=0d0
1070 prsoft=1.-prhard
1071 IF(prsoft.LT.0.1) THEN
1072 print *, 'ERROR IN DRCMU; PRSOFT=',prsoft
1073 stop
1074 ENDIF
1075C
1076 rr5=rrr(5)
1077 ihard=(rr5.GT.prsoft)
1078 IF (ihard) THEN
1079C TAU DECAY TO 'TAU+photon'
1080 rr1=rrr(1)
1081 ams1=(amu+amnuta)**2
1082 ams2=(amtax)**2
1083 xk1=1-ams1/ams2
1084 xl1=log(xk1/2/xk0dec)
1085 xl0=log(2*xk0dec)
1086 xk=exp(xl1*rr1+xl0)
1087 am3sq=(1-xk)*ams2
1088 am3 =sqrt(am3sq)
1089 phspac=phspac*ams2*xl1*xk
1090 phspac=phspac/prhard
1091 ELSE
1092 am3=amtax
1093 phspac=phspac*2**6*pi**3
1094 phspac=phspac/prsoft
1095 ENDIF
1096C MASS OF NEUTRINA SYSTEM
1097 rr2=rrr(2)
1098 ams1=(amnuta)**2
1099 ams2=(am3-amu)**2
1100CAM
1101CAM
1102* FLAT PHASE SPACE;
1103 am2sq=ams1+ rr2*(ams2-ams1)
1104 am2 =sqrt(am2sq)
1105 phspac=phspac*(ams2-ams1)
1106* NEUTRINA REST FRAME, DEFINE XN AND XA
1107 enq1=(am2sq+amnuta**2)/(2*am2)
1108 enq2=(am2sq-amnuta**2)/(2*am2)
1109 ppi= enq1**2-amnuta**2
1110 pppi=sqrt(abs(enq1**2-amnuta**2))
1111 phspac=phspac*(4*pi)*(2*pppi/am2)
1112* NU TAU IN NUNU REST FRAME
1113 CALL spherd(pppi,xn)
1114 xn(4)=enq1
1115* NU LIGHT IN NUNU REST FRAME
1116 DO 30 i=1,3
1117 30 xa(i)=-xn(i)
1118 xa(4)=enq2
1119* TAU-prim REST FRAME, DEFINE QP (muon
1120* NUNU MOMENTUM
1121 pr(1)=0
1122 pr(2)=0
1123 pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1124 pr(3)= sqrt(abs(pr(4)**2-am2**2))
1125 ppi = pr(4)**2-am2**2
1126* MUON MOMENTUM
1127 qp(1)=0
1128 qp(2)=0
1129 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1130 qp(3)=-pr(3)
1131 phspac=phspac*(4*pi)*(2*pr(3)/am3)
1132* NEUTRINA BOOSTED FROM THEIR FRAME TO TAU-prim REST FRAME
1133 exe=(pr(4)+pr(3))/am2
1134 CALL bostd3(exe,xn,xn)
1135 CALL bostd3(exe,xa,xa)
1136 rr3=rrr(3)
1137 rr4=rrr(4)
1138 IF (ihard) THEN
1139 eps=4*(amu/amtax)**2
1140 xl1=log((2+eps)/eps)
1141 xl0=log(eps)
1142 eta =exp(xl1*rr3+xl0)
1143 cthet=1+eps-eta
1144 thet =acos(cthet)
1145 phspac=phspac*xl1/2*eta
1146 phi = 2*pi*rr4
1147 CALL rotpox(thet,phi,xn)
1148 CALL rotpox(thet,phi,xa)
1149 CALL rotpox(thet,phi,qp)
1150 CALL rotpox(thet,phi,pr)
1151C
1152* NOW TO THE TAU REST FRAME, DEFINE TAU-prim AND GAMMA MOMENTA
1153* tau-prim MOMENTUM
1154 paa(1)=0
1155 paa(2)=0
1156 paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1157 paa(3)= sqrt(abs(paa(4)**2-am3**2))
1158 ppi = paa(4)**2-am3**2
1159 phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1160* GAMMA MOMENTUM
1161 ph(1)=0
1162 ph(2)=0
1163 ph(4)=paa(3)
1164 ph(3)=-paa(3)
1165* ALL MOMENTA BOOSTED FROM TAU-prim REST FRAME TO TAU REST FRAME
1166* Z-AXIS ANTIPARALLEL TO PHOTON MOMENTUM
1167 exe=(paa(4)+paa(3))/am3
1168 CALL bostd3(exe,xn,xn)
1169 CALL bostd3(exe,xa,xa)
1170 CALL bostd3(exe,qp,qp)
1171 CALL bostd3(exe,pr,pr)
1172 ELSE
1173 thet =acos(-1.+2*rr3)
1174 phi = 2*pi*rr4
1175 CALL rotpox(thet,phi,xn)
1176 CALL rotpox(thet,phi,xa)
1177 CALL rotpox(thet,phi,qp)
1178 CALL rotpox(thet,phi,pr)
1179C
1180* NOW TO THE TAU REST FRAME, DEFINE TAU-prim AND GAMMA MOMENTA
1181* tau-prim MOMENTUM
1182 paa(1)=0
1183 paa(2)=0
1184 paa(4)=amtax
1185 paa(3)=0
1186* GAMMA MOMENTUM
1187 ph(1)=0
1188 ph(2)=0
1189 ph(4)=0
1190 ph(3)=0
1191 ENDIF
1192C PARTIAL WIDTH CONSISTS OF PHASE SPACE AND AMPLITUDE
1193 CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1194 dgamt=1/(2.*amtax)*amplit*phspac
1195 END
1196 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1197 IMPLICIT REAL*8 (a-h,o-z)
1198C ----------------------------------------------------------------------
1199C IT CALCULATES MATRIX ELEMENT FOR THE
1200C TAU --> MU(E) NU NUBAR DECAY MODE
1201C INCLUDING COMPLETE ORDER ALPHA QED CORRECTIONS.
1202C ----------------------------------------------------------------------
1203 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1204 * ,ampiz,ampi,amro,gamro,ama1,gama1
1205 * ,amk,amkz,amkst,gamkst
1206C
1207 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1208 * ,ampiz,ampi,amro,gamro,ama1,gama1
1209 * ,amk,amkz,amkst,gamkst
1210 real*8 hv(4),qp(4),xn(4),xa(4),xk(4)
1211C
1212 hv(4)=1.d0
1213 ak0=xk0dec*amtau
1214 IF(xk(4).LT.0.1d0*ak0) THEN
1215 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1216 ELSE
1217 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1218 ENDIF
1219 RETURN
1220 END
1221 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1222C
1223C **********************************************************************
1224C REAL PHOTON MATRIX ELEMENT SQUARED *
1225C PARAMETERS: *
1226C HV- POLARIMETRIC FOUR-VECTOR OF TAU *
1227C QP,XN,XA,XK - 4-momenta of electron (muon), NU, NUBAR and PHOTON *
1228C All four-vectors in TAU rest frame (in GeV) *
1229C AK0 - INFRARED CUTOFF, MINIMAL ENERGY OF HARD PHOTONS (GEV) *
1230C SQM2 - value for S=0 *
1231C see Eqs. (2.9)-(2.10) from CJK ( Nucl.Phys.B(1991) ) *
1232C **********************************************************************
1233C
1234 IMPLICIT REAL*8(a-h,o-z)
1235 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1236 * ,ampiz,ampi,amro,gamro,ama1,gama1
1237 * ,amk,amkz,amkst,gamkst
1238C
1239 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1240 * ,ampiz,ampi,amro,gamro,ama1,gama1
1241 * ,amk,amkz,amkst,gamkst
1242 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1243 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1244 COMMON / qedprm /alfinv,alfpi,xk0
1245 real*8 alfinv,alfpi,xk0
1246 real*8 qp(4),xn(4),xa(4),xk(4)
1247 real*8 r(4)
1248 real*8 hv(4)
1249 real*8 s0(3),rxa(3),rxk(3),rqp(3)
1250 DATA pi /3.141592653589793238462643d0/
1251C
1252 tmass=amtau
1253 gf=gfermi
1254 alphai=alfinv
1255 tmass2=tmass**2
1256 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1257 r(4)=tmass
1258C SCALAR PRODUCTS OF FOUR-MOMENTA
1259 DO 7 i=1,3
1260 r(1)=0.d0
1261 r(2)=0.d0
1262 r(3)=0.d0
1263 r(i)=tmass
1264 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1265C RXN(I)=R(4)*XN(4)-R(1)*XN(1)-R(2)*XN(2)-R(3)*XN(3)
1266 rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1267 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1268 7 CONTINUE
1269 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1270 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1271 qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1272c XNXA=XN(4)*XA(4)-XN(1)*XA(1)-XN(2)*XA(2)-XN(3)*XA(3)
1273 xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1274 xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1275 txn=tmass*xn(4)
1276 txa=tmass*xa(4)
1277 tqp=tmass*qp(4)
1278 txk=tmass*xk(4)
1279C
1280 x= xnxk/qpxn
1281 z= txk/tqp
1282 a= 1+x
1283 b= 1+ x*(1+z)/2+z/2
1284 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1285 $tmass2/txk**2) +
1286 $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1287 $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1288 const4=256*pi/alphai*gf**2
1289 IF (itdkrc.EQ.0) const4=0d0
1290 sqm2=s1*const4
1291 DO 5 i=1,3
1292 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1293 $ tmass2/txk**2) +
1294 $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1295 $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1296 5 hv(i)=s0(i)/s1-1.d0
1297 RETURN
1298 END
1299 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1300C
1301C **********************************************************************
1302C BORN +VIRTUAL+SOFT PHOTON MATRIX ELEMENT**2 O(ALPHA) *
1303C PARAMETERS: *
1304C HV- POLARIMETRIC FOUR-VECTOR OF TAU *
1305C QP,XN,XA - FOUR-MOMENTA OF ELECTRON (MUON), NU AND NUBAR IN GEV *
1306C ALL FOUR-VECTORS IN TAU REST FRAME *
1307C AK0 - INFRARED CUTOFF, MINIMAL ENERGY OF HARD PHOTONS *
1308C THB - VALUE FOR S=0 *
1309C SEE EQS. (2.2),(2.4)-(2.5) FROM CJK (NUCL.PHYS.B351(1991)70 *
1310C AND (C.2) FROM JK (NUCL.PHYS.B320(1991)20 ) *
1311C **********************************************************************
1312C
1313 IMPLICIT REAL*8(a-h,o-z)
1314 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1315 * ,ampiz,ampi,amro,gamro,ama1,gama1
1316 * ,amk,amkz,amkst,gamkst
1317C
1318 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1319 * ,ampiz,ampi,amro,gamro,ama1,gama1
1320 * ,amk,amkz,amkst,gamkst
1321 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1322 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1323 COMMON / qedprm /alfinv,alfpi,xk0
1324 real*8 alfinv,alfpi,xk0
1325 dimension qp(4),xn(4),xa(4)
1326 real*8 hv(4)
1327 dimension r(4)
1328 real*8 rxa(3),rxn(3),rqp(3)
1329 real*8 bornpl(3),am3pol(3),xm3pol(3)
1330 DATA pi /3.141592653589793238462643d0/
1331C
1332 tmass=amtau
1333 gf=gfermi
1334 alphai=alfinv
1335C
1336 tmass2=tmass**2
1337 r(4)=tmass
1338 DO 7 i=1,3
1339 r(1)=0.d0
1340 r(2)=0.d0
1341 r(3)=0.d0
1342 r(i)=tmass
1343 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1344 rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1345C RXK(I)=R(4)*XK(4)-R(1)*XK(1)-R(2)*XK(2)-R(3)*XK(3)
1346 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1347 7 CONTINUE
1348C QUASI TWO-BODY VARIABLES
1349 u0=qp(4)/tmass
1350 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1351 w3=u3
1352 w0=(xn(4)+xa(4))/tmass
1353 up=u0+u3
1354 um=u0-u3
1355 wp=w0+w3
1356 wm=w0-w3
1357 yu=log(up/um)/2
1358 yw=log(wp/wm)/2
1359 eps2=u0**2-u3**2
1360 eps=sqrt(eps2)
1361 y=w0**2-w3**2
1362 al=ak0/tmass
1363C FORMFACTORS
1364 f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1365 $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1366 $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1367 $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1368 fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1369 fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1370 f3= eps2*(fp+fm)/2
1371C SCALAR PRODUCTS OF FOUR-MOMENTA
1372 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1373 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1374 xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1375 txn=tmass*xn(4)
1376 txa=tmass*xa(4)
1377 tqp=tmass*qp(4)
1378C DECAY DIFFERENTIAL WIDTH WITHOUT AND WITH POLARIZATION
1379 const3=1/(2*alphai*pi)*64*gf**2
1380 IF (itdkrc.EQ.0) const3=0d0
1381 xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1382 $fm* qpxn*qpxa + f3* tmass2*xnxa )
1383 am3=xm3*const3
1384C V-A AND V+A COUPLINGS, BUT IN THE BORN PART ONLY
1385 brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1386 & -(gv**2-ga**2)*tmass*amnuta*qpxa
1387 born= 32*(gfermi**2/2.)*brak
1388 DO 5 i=1,3
1389 xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1390 $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1391 $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1392 am3pol(i)=xm3pol(i)*const3
1393C V-A AND V+A COUPLINGS, BUT IN THE BORN PART ONLY
1394 bornpl(i)=born+(
1395 & (gv+ga)**2*tmass*xnxa*qp(i)
1396 & -(gv-ga)**2*tmass*qpxn*xa(i)
1397 & +(gv**2-ga**2)*amnuta*txa*qp(i)
1398 & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1399 & 32*(gfermi**2/2.)
1400 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1401 thb=born+am3
1402 IF (thb/born.LT.0.1d0) THEN
1403 print *, 'ERROR IN THB, THB/BORN=',thb/born
1404 thb=0.d0
1405 ENDIF
1406 RETURN
1407 END
1408 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1409C ----------------------------------------------------------------------
1410C TAU DECAY INTO PION AND TAU-NEUTRINO
1411C IN TAU REST FRAME
1412C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
1413C PPI PION CHARGED
1414C ----------------------------------------------------------------------
1415 REAL POL(4),HV(4),PNU(4),PPI(4),RN(1)
1416CC
1417 IF(mode.EQ.-1) THEN
1418C ===================
1419 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1420CC CALL HBOOK1(815,'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
1421
1422 ELSEIF(mode.EQ. 0) THEN
1423C =======================
1424300 CONTINUE
1425 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1426 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1427CC CALL HFILL(815,WT)
1428 CALL ranmar(rn,1)
1429 IF(rn(1).GT.wt) GOTO 300
1430C
1431 ELSEIF(mode.EQ. 1) THEN
1432C =======================
1433 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1434CC CALL HPRINT(815)
1435 ENDIF
1436C =====
1437 RETURN
1438 END
1439 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1440C ----------------------------------------------------------------------
1441 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1442 * ,ampiz,ampi,amro,gamro,ama1,gama1
1443 * ,amk,amkz,amkst,gamkst
1444C
1445 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1446 * ,ampiz,ampi,amro,gamro,ama1,gama1
1447 * ,amk,amkz,amkst,gamkst
1448 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1449 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1450 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1451 real*4 gampmc ,gamper
1452 COMMON / inout / inut,iout
1453 REAL PPI(4),PNU(4),HV(4)
1454 DATA pi /3.141592653589793238462643/
1455C
1456 IF(mode.EQ.-1) THEN
1457C ===================
1458 nevtot=0
1459 ELSEIF(mode.EQ. 0) THEN
1460C =======================
1461 nevtot=nevtot+1
1462 epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1463 enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1464 xpi= sqrt(epi**2-ampi**2)
1465C PI MOMENTUM
1466 CALL sphera(xpi,ppi)
1467 ppi(4)=epi
1468C TAU-NEUTRINO MOMENTUM
1469 DO 30 i=1,3
147030 pnu(i)=-ppi(i)
1471 pnu(4)=enu
1472 pxq=amtau*epi
1473 pxn=amtau*enu
1474 qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1475 brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1476 & +(gv**2-ga**2)*amtau*amnuta*ampi**2
1477 DO 40 i=1,3
147840 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1479 hv(4)=1
1480C
1481 ELSEIF(mode.EQ. 1) THEN
1482C =======================
1483 IF(nevtot.EQ.0) RETURN
1484 fpi=0.1284
1485C GAMM=(GFERMI*FPI)**2/(16.*PI)*AMTAU**3*
1486C * (BRAK/AMTAU**4)**2
1487CZW 7.02.93 here was an error affecting non standard model
1488C configurations only
1489 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1490 $ (brak/amtau**4)*
1491 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1492 $ -4*ampi**2*amnuta**2 )/amtau**2
1493 error=0
1494 rat=gamm/gamel
1495 WRITE(iout, 7010) nevtot,gamm,rat,error
1496 gampmc(3)=rat
1497 gamper(3)=error
1498CAM NEVDEC(3)=NEVTOT
1499 ENDIF
1500C =====
1501 RETURN
1502 7010 FORMAT(///1x,15(5h*****)
1503 $ /,' *', 25x,'******** DADMPI FINAL REPORT ******** ',9x,1h*
1504 $ /,' *',i20 ,5x,'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1505 $ /,' *',e20.5,5x,'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1506 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1507 $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1508 $ /,1x,15(5h*****)/)
1509 END
1510 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1511C ----------------------------------------------------------------------
1512C THIS SIMULATES TAU DECAY IN TAU REST FRAME
1513C INTO NU RHO, NEXT RHO DECAYS INTO PION PAIR.
1514C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
1515C PRO RHO
1516C PIC PION CHARGED
1517C PIZ PION ZERO
1518C ----------------------------------------------------------------------
1519 COMMON / inout / inut,iout
1520 REAL POL(4),HV(4),PRO(4),PNU(4),PIC(4),PIZ(4),RN(1)
1521 DATA iwarm/0/
1522C
1523 IF(mode.EQ.-1) THEN
1524C ===================
1525 iwarm=1
1526 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1527CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXRO $',100,0,2)
1528CC CALL HBOOK1(916,'ABS2 OF HV IN ROUTINE DEXRO $',100,0,2)
1529C
1530 ELSEIF(mode.EQ. 0) THEN
1531C =======================
1532300 CONTINUE
1533 IF(iwarm.EQ.0) GOTO 902
1534 CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1535 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1536CC CALL HFILL(816,WT)
1537CC XHELP=HV(1)**2+HV(2)**2+HV(3)**2
1538CC CALL HFILL(916,XHELP)
1539 CALL ranmar(rn,1)
1540 IF(rn(1).GT.wt) GOTO 300
1541C
1542 ELSEIF(mode.EQ. 1) THEN
1543C =======================
1544 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1545CC CALL HPRINT(816)
1546CC CALL HPRINT(916)
1547 ENDIF
1548C =====
1549 RETURN
1550 902 WRITE(iout, 9020)
1551 9020 FORMAT(' ----- DEXRO: LACK OF INITIALISATION')
1552 stop
1553 END
1554 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1555C ----------------------------------------------------------------------
1556 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1557 * ,ampiz,ampi,amro,gamro,ama1,gama1
1558 * ,amk,amkz,amkst,gamkst
1559C
1560 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1561 * ,ampiz,ampi,amro,gamro,ama1,gama1
1562 * ,amk,amkz,amkst,gamkst
1563 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1564 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1565 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1566 real*4 gampmc ,gamper
1567 COMMON / inout / inut,iout
1568 REAL HHV(4)
1569 REAL HV(4),PRO(4),PNU(4),PIC(4),PIZ(4)
1570 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
1571 real*4 rrr(3)
1572 real*8 swt, sswt
1573 DATA pi /3.141592653589793238462643/
1574 DATA iwarm/0/
1575C
1576 IF(mode.EQ.-1) THEN
1577C ===================
1578 iwarm=1
1579 nevraw=0
1580 nevacc=0
1581 nevovr=0
1582 swt=0
1583 sswt=0
1584 wtmax=1e-20
1585 DO 15 i=1,500
1586 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1587 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
158815 CONTINUE
1589CC CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADMRO $',100,0,2)
1590CC PRINT 7003,WTMAX
1591C
1592 ELSEIF(mode.EQ. 0) THEN
1593C =======================
1594300 CONTINUE
1595 IF(iwarm.EQ.0) GOTO 902
1596 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1597CC CALL HFILL(801,WT/WTMAX)
1598 nevraw=nevraw+1
1599 swt=swt+wt
1600 sswt=sswt+wt**2
1601 CALL ranmar(rrr,3)
1602 rn=rrr(1)
1603 IF(wt.GT.wtmax) nevovr=nevovr+1
1604 IF(rn*wtmax.GT.wt) GOTO 300
1605C ROTATIONS TO BASIC TAU REST FRAME
1606 costhe=-1.+2.*rrr(2)
1607 thet=acos(costhe)
1608 phi =2*pi*rrr(3)
1609 CALL rotor2(thet,pnu,pnu)
1610 CALL rotor3( phi,pnu,pnu)
1611 CALL rotor2(thet,pro,pro)
1612 CALL rotor3( phi,pro,pro)
1613 CALL rotor2(thet,pic,pic)
1614 CALL rotor3( phi,pic,pic)
1615 CALL rotor2(thet,piz,piz)
1616 CALL rotor3( phi,piz,piz)
1617 CALL rotor2(thet,hv,hv)
1618 CALL rotor3( phi,hv,hv)
1619 DO 44 i=1,3
1620 44 hhv(i)=-isgn*hv(i)
1621 nevacc=nevacc+1
1622C
1623 ELSEIF(mode.EQ. 1) THEN
1624C =======================
1625 IF(nevraw.EQ.0) RETURN
1626 pargam=swt/float(nevraw+1)
1627 error=0
1628 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1629 rat=pargam/gamel
1630 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1631CC CALL HPRINT(801)
1632 gampmc(4)=rat
1633 gamper(4)=error
1634CAM NEVDEC(4)=NEVACC
1635 ENDIF
1636C =====
1637 RETURN
1638 7003 FORMAT(///1x,15(5h*****)
1639 $ /,' *', 25x,'******** DADMRO INITIALISATION ********',9x,1h*
1640 $ /,' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1641 $ /,1x,15(5h*****)/)
1642 7010 FORMAT(///1x,15(5h*****)
1643 $ /,' *', 25x,'******** DADMRO FINAL REPORT ******** ',9x,1h*
1644 $ /,' *',i20 ,5x,'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1645 $ /,' *',i20 ,5x,'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1646 $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1647 $ /,' *',e20.5,5x,'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1648 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1649 $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1650 $ /,1x,15(5h*****)/)
1651 902 WRITE(iout, 9020)
1652 9020 FORMAT(' ----- DADMRO: LACK OF INITIALISATION')
1653 stop
1654 END
1655 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1656C ----------------------------------------------------------------------
1657C IT SIMULATES RHO DECAY IN TAU REST FRAME WITH
1658C Z-AXIS ALONG RHO MOMENTUM
1659C ----------------------------------------------------------------------
1660 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1661 * ,ampiz,ampi,amro,gamro,ama1,gama1
1662 * ,amk,amkz,amkst,gamkst
1663C
1664 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1665 * ,ampiz,ampi,amro,gamro,ama1,gama1
1666 * ,amk,amkz,amkst,gamkst
1667 REAL HV(4),PT(4),PN(4),PR(4),PIC(4),PIZ(4),QQ(4),RR1(1)
1668 DATA pi /3.141592653589793238462643/
1669 DATA icont /0/
1670C
1671C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
1672 phspac=1./2**11/pi**5
1673C TAU MOMENTUM
1674 pt(1)=0.
1675 pt(2)=0.
1676 pt(3)=0.
1677 pt(4)=amtau
1678C MASS OF (REAL/VIRTUAL) RHO
1679 ams1=(ampi+ampiz)**2
1680 ams2=(amtau-amnuta)**2
1681C FLAT PHASE SPACE
1682C AMX2=AMS1+ RR1*(AMS2-AMS1)
1683C AMX=SQRT(AMX2)
1684C PHSPAC=PHSPAC*(AMS2-AMS1)
1685C PHASE SPACE WITH SAMPLING FOR RHO RESONANCE
1686 alp1=atan((ams1-amro**2)/amro/gamro)
1687 alp2=atan((ams2-amro**2)/amro/gamro)
1688CAM
1689 100 CONTINUE
1690 CALL ranmar(rr1,1)
1691 alp=alp1+rr1(1)*(alp2-alp1)
1692 amx2=amro**2+amro*gamro*tan(alp)
1693 amx=sqrt(amx2)
1694 IF(amx.LT.2.*ampi) GO TO 100
1695CAM
1696 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1697 phspac=phspac*(alp2-alp1)
1698C
1699C TAU-NEUTRINO MOMENTUM
1700 pn(1)=0
1701 pn(2)=0
1702 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1703 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
1704C RHO MOMENTUM
1705 pr(1)=0
1706 pr(2)=0
1707 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1708 pr(3)=-pn(3)
1709 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1710C
1711CAM
1712 enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1713 enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1714 pppi=sqrt((enq1-ampi)*(enq1+ampi))
1715 phspac=phspac*(4*pi)*(2*pppi/amx)
1716C CHARGED PI MOMENTUM IN RHO REST FRAME
1717 CALL sphera(pppi,pic)
1718 pic(4)=enq1
1719C NEUTRAL PI MOMENTUM IN RHO REST FRAME
1720 DO 20 i=1,3
172120 piz(i)=-pic(i)
1722 piz(4)=enq2
1723 exe=(pr(4)+pr(3))/amx
1724C PIONS BOOSTED FROM RHO REST FRAME TO TAU REST FRAME
1725 CALL bostr3(exe,pic,pic)
1726 CALL bostr3(exe,piz,piz)
1727
1728 CALL dam2pi(0,pt,pn,pic,piz,amplit,hv)
1729 dgamt=1/(2.*amtau)*amplit*phspac
1730
1731 RETURN
1732 END
1733
1734
1735C ----------------------------------------------------------------------
1736C ----------------------------------------------------------------------
1737C ----------------------------------------------------------------------
1738C RCHL UPDATE - NEW FUNCTIONS
1739C ----------------------------------------------------------------------
1740C ----------------------------------------------------------------------
1741C ----------------------------------------------------------------------
1742
1743
1744 SUBROUTINE dam2pi(MNUM,PT,PN,PIM1,PIM2,AMPLIT,HV)
1745C ----------------------------------------------------------------------
1746* CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
1747* FOR TAU DECAY INTO 2 scalar MODES
1748* ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
1749* CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
1750C MNUM DECAY MODE IDENTIFIER.
1751C
1752C called by : DPHSAA
1753C ----------------------------------------------------------------------
1754 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1755 * ,ampiz,ampi,amro,gamro,ama1,gama1
1756 * ,amk,amkz,amkst,gamkst
1757C
1758 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1759 * ,ampiz,ampi,amro,gamro,ama1,gama1
1760 * ,amk,amkz,amkst,gamkst
1761 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1762 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1763 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4)
1764 REAL PIVEC(4),PIAKS(4),HVM(4)
1765 COMPLEX HADCUR(4)
1766 DATA pi /3.141592653589793238462643/
1767 DATA icont /0/
1768C
1769 IF (mnum.EQ.0) THEN
1770 CALL curr_pipi0(pim1,pim2,hadcur)
1771 ELSEIF (mnum.EQ.1) THEN
1772 CALL curr_pik0(pim1,pim2,hadcur)
1773 ELSEIF (mnum.EQ.2) THEN
1774 CALL curr_kpi0(pim1,pim2,hadcur)
1775 ELSEIF (mnum.EQ.3) THEN
1776 CALL curr_kk0(pim1,pim2,hadcur)
1777 ELSE
1778 write(*,*) 'DAM2PI: wrong MNUM= ',mnum
1779 stop
1780 ENDIF
1781
1782C
1783* CALCULATE PI-VECTORS: VECTOR AND AXIAL
1784 CALL clvec(hadcur,pn,pivec)
1785 CALL claxi(hadcur,pn,piaks)
1786 CALL clnut(hadcur,brakm,hvm)
1787* SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
1788 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
1789 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
1790 IF (mnum.EQ.0.OR.mnum.EQ.3) THEN
1791 amplit=(ccabib*gfermi)**2*brak
1792 ELSE
1793 amplit=(scabib*gfermi)**2*brak
1794 ENDIF
1795C POLARIMETER VECTOR IN TAU REST FRAME
1796 DO 90 i=1,3
1797 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
1798 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
1799C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
1800 hv(i)=-hv(i)/brak
1801 90 CONTINUE
1802 END
1803
1804
1805 SUBROUTINE curr_pipi0(PC,PN,HADCUR)
1806C standard TAUOLA current for tau to pi pi0 nu decay
1807C now it has universal form eg. it is straighforward to add
1808C scalar part
1809C NOTE:
1810C PC 4-momentum of pi
1811C PN 4-momentum of pi0
1812C 06.08.2011
1813 IMPLICIT NONE
1814 COMPLEX BWIGS,HADCUR(4),FKPIPL,FRHO_PI
1815 COMPLEX*16 FPIBEL
1816 REAL PC(4),PN(4),QQ(4),PKS(4),FPIRHO
1817 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1818 * ,ampiz,ampi,amro,gamro,ama1,gama1
1819 * ,amk,amkz,amkst,gamkst
1820C
1821 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1822 * ,ampiz,ampi,amro,gamro,ama1,gama1
1823 * ,amk,amkz,amkst,gamkst
1824 COMMON /ipcht/ iver
1825 INTEGER IVER
1826 INTEGER FF2PIRHO
1827
1828
1829 REAL PKSD,QQPKS
1830 INTEGER IK,K
1831 DO ik=1,4
1832 pks(ik)=pc(ik)+ pn(ik)
1833 qq(ik)=pc(ik)- pn(ik)
1834 ENDDO
1835C QQ transverse to PKS
1836 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1837 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1838 DO 31 ik=1,4
1839 31 qq(ik)=qq(ik)-pks(ik)*qqpks/pksd
1840
1841 IF (iver.EQ.1) THEN
1842 CALL getff2pirho(ff2pirho)
1843 IF (ff2pirho.EQ.2) THEN ! Belle,
1844C ! all fit parameters, par(1...11), are free
1845 DO k=1,4
1846 hadcur(k)=qq(k)* fpibel(sqrt(pksd),0)
1847 ENDDO
1848 ELSEIF (ff2pirho.EQ.3) THEN ! Belle
1849c ! all fit parameter free except for
1850c ! par(1)=F_pi(0)=1-fixed
1851
1852 DO k=1,4
1853 hadcur(k)=qq(k)* fpibel(sqrt(pksd),1)
1854 ENDDO
1855 ELSE
1856 write(*,*) 'problem in 2-scalars current FF2PIRHO=',ff2pirho
1857 stop
1858 ENDIF
1859 ELSEIF (iver.EQ.0) THEN ! cleo
1860 DO k=1,4
1861 hadcur(k)=qq(k)* sqrt(fpirho(sqrt(pksd)))
1862 ENDDO
1863
1864 ELSE
1865 write(*,*) 'problem in 2-scalars current IVER=',iver
1866 stop
1867 ENDIF
1868
1869 END
1870
1871
1872 SUBROUTINE curr_pik0(PC,PN,HADCUR)
1873C standard TAUOLA current for tau to pi K0 nu decay
1874C now it has universal form eg. it is straighforward to add
1875C scalar part
1876C NOTE:
1877C PC 4-momentum of pi
1878C PN 4-momentum of K0
1879C 06.08.2011
1880 implicit none
1881 COMPLEX BWIGS,HADCUR(4),FKPIPL
1882 REAL PC(4),PN(4),QQ(4),PKS(4),FKPISC,PKSD,QQPKS
1883 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1884 & ,ampiz,ampi,amro,gamro,ama1,gama1
1885 & ,amk,amkz,amkst,gamkst
1886C
1887 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1888 & ,ampiz,ampi,amro,gamro,ama1,gama1
1889 & ,amk,amkz,amkst,gamkst,fact_k0pi
1890 COMMON / taukle / bra1,brk0,brk0b,brks
1891 real*4 bra1,brk0,brk0b,brks
1892 Integer I,K
1893
1894 DO i=1,4
1895 pks(i)=pc(i)+ pn(i)
1896 qq(i)=pc(i)- pn(i)
1897 ENDDO
1898
1899C QQ transverse to PKS
1900 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1901 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1902 DO 31 i=1,4
1903 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
1904
1905
1906 DO k=1,4
1907 hadcur(k)=qq(k)*bwigs(pksd,amkst,gamkst)
1908 ENDDO
1909
1910 END
1911
1912
1913 SUBROUTINE curr_kpi0(PC,PN,HADCUR)
1914C standard TAUOLA current for tau to pi pi0 nu decay
1915C now it has universal form eg. it is straighforward to add
1916C scalar part
1917C NOTE:
1918C PC 4-momentum of K
1919C PN 4-momentum of pi0
1920C 06.08.2011
1921 implicit none
1922 COMPLEX BWIGS,HADCUR(4),FKPIPL
1923 REAL PC(4),PN(4),QQ(4),PKS(4),FKPISC,PKSD,QQPKS
1924 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1925 * ,ampiz,ampi,amro,gamro,ama1,gama1
1926 * ,amk,amkz,amkst,gamkst
1927C
1928 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1929 * ,ampiz,ampi,amro,gamro,ama1,gama1
1930 * ,amk,amkz,amkst,gamkst,fact_kpi0
1931 COMMON / taukle / bra1,brk0,brk0b,brks
1932 real*4 bra1,brk0,brk0b,brks
1933 INTEGER I,K
1934
1935 DO 30 i=1,4
1936 pks(i)=pc(i)+ pn(i)
1937 30 qq(i)=pc(i)- pn(i)
1938C QQ transverse to PKS
1939 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1940 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1941 DO 31 i=1,4
1942 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
1943 DO k=1,4
1944 hadcur(k)=qq(k)*bwigs(pksd,amkst,gamkst)
1945 ENDDO
1946
1947 END
1948
1949
1950 SUBROUTINE curr_kk0(PC,PN,HADCUR)
1951C standard TAUOLA current for tau to K K0 nu decay
1952C now it has universal form eg. it is straighforward to add
1953C scalar part
1954C NOTE:
1955C PC 4-momentum of K
1956C PN 4-momentum of K0
1957C 06.08.2011
1958 IMPLICIT NONE
1959 COMPLEX BWIGS,HADCUR(4),FKK0_RCHT
1960 REAL PC(4),PN(4),QQ(4),PKS(4),PKSD,QQPKS,FPIRK
1961 INTEGER I,K
1962 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1963 * ,ampiz,ampi,amro,gamro,ama1,gama1
1964 * ,amk,amkz,amkst,gamkst
1965C
1966 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1967 * ,ampiz,ampi,amro,gamro,ama1,gama1
1968 * ,amk,amkz,amkst,gamkst
1969 DO i=1,4
1970 pks(i)=pc(i)+ pn(i)
1971 qq(i)=pc(i)- pn(i)
1972 ENDDO
1973C QQ transverse to PKS
1974 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1975 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1976 DO 31 i=1,4
1977 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
1978
1979 DO k=1,4
1980 hadcur(k)=qq(k)*sqrt(fpirk(sqrt(pksd)))
1981 ENDDO
1982
1983 END
1984
1985
1986 FUNCTION coef(I,J)
1987C clebsh gordan (or so ...) coefs for 3 scalar final states
1988 implicit none
1989C IVER=0 TAUOLA cleo COEF(I,J) = COEFc(I,J)
1990C IVER=1 TAUOLA RChL COEF(I,J) = COEFr(I,J)
1991 COMMON /ipcht/ iver
1992 INTEGER IVER
1993 REAL COEFc(1:5,0:7)
1994 REAL COEFr(1:5,0:7)
1995 REAL COEF,COEFrr
1996 DATA pi /3.141592653589793238462643/
1997 REAL PI
1998 DATA icont /0/
1999 INTEGER ICONT
2000 INTEGER I,J
2001 REAL FPIc,FPIr
2002
2003C initialization of FPI matrix defined in ...
2004C FPIc is to be used with cleo initialization
2005
2006C actual choice is made in ???
2007
2008 DATA fpic /93.3e-3/
2009
2010
2011C initialization of COEF matrix defined in ...
2012C COEFc is to be used with cleo initialization
2013
2014 IF (icont.EQ.0) THEN
2015 icont=1
2016C
2017C*****COEFc(I,J)
2018
2019 coefc(1,0)= 2.0*sqrt(2.)/3.0
2020 coefc(2,0)=-2.0*sqrt(2.)/3.0
2021C AJW 2/98: Add in the D-wave and I=0 3pi substructure:
2022 coefc(3,0)= 2.0*sqrt(2.)/3.0
2023 coefc(4,0)= fpic
2024 coefc(5,0)= 0.0
2025C
2026 coefc(1,1)=-sqrt(2.)/3.0
2027 coefc(2,1)= sqrt(2.)/3.0
2028 coefc(3,1)= 0.0
2029 coefc(4,1)= fpic
2030 coefc(5,1)= sqrt(2.)
2031
2032C
2033 coefc(1,2)=-sqrt(2.)/3.0
2034 coefc(2,2)= sqrt(2.)/3.0
2035 coefc(3,2)= 0.0
2036
2037 coefc(4,2)= 0.0
2038 coefc(5,2)=-sqrt(2.)
2039
2040
2041C AJW 11/97: Add in the K*-prim-s, ala Finkemeier&Mirkes
2042 coefc(1,3)= 1./3.
2043 coefc(2,3)=-2./3.
2044 coefc(3,3)= 2./3.
2045 coefc(4,3)= 0.0
2046 coefc(5,3)= 0.0
2047C
2048 coefc(1,4)= 1.0/sqrt(2.)/3.0
2049 coefc(2,4)=-1.0/sqrt(2.)/3.0
2050 coefc(3,4)= 0.0
2051 coefc(4,4)= 0.0
2052 coefc(5,4)= 0.0
2053C
2054 coefc(1,5)=-sqrt(2.)/3.0
2055 coefc(2,5)= sqrt(2.)/3.0
2056 coefc(3,5)= 0.0
2057 coefc(4,5)= 0.0
2058 coefc(5,5)=-sqrt(2.)
2059C
2060C AJW 11/97: Add in the K*-prim-s, ala Finkemeier&Mirkes
2061 coefc(1,6)= 1./3.
2062 coefc(2,6)=-2./3.
2063 coefc(3,6)= 2./3.
2064 coefc(4,6)= 0.0
2065 coefc(5,6)=-2.0
2066C
2067 coefc(1,7)= 0.0
2068 coefc(2,7)= 0.0
2069 coefc(3,7)= 0.0
2070 coefc(4,7)= 0.0
2071 coefc(5,7)=-sqrt(2.0/3.0)
2072
2073 ENDIF
2074 IF (iver.EQ.0.OR.j.NE.0) THEN ! so far rchl only for 3pi modes
2075 coef=coefc(i,j)
2076 ELSEIF (iver.EQ.1) THEN
2077 coef=coefrr(i,j)
2078 ELSE
2079 write(*,*) 'wrong IVER=',iver
2080 stop
2081 ENDIF
2082 END
2083
2084
2085 SUBROUTINE inirchl(IVERI)
2086C routine to set version no for the currents physics initialization
2087C IVER=0 TAUOLA cleo
2088C IVER=1 TAUOLA RChL
2089
2090 implicit none
2091 INTEGER IVERI
2092
2093 COMMON /ipcht/ iver
2094 INTEGER IVER
2095 iver=iveri
2096
2097 IF (iver.EQ.1) THEN
2098 CALL rchl_parameters(1)
2099 ENDIF
2100 end
2101
2102
2103 SUBROUTINE inirchlget(I)
2104C routine to get version no for the currents physics initialization
2105C IVER=0 TAUOLA cleo
2106C IVER=1 TAUOLA RChL
2107 COMMON /ipcht/ iver
2108 INTEGER IVER
2109
2110 i=iver
2111 end
2112
2113C ----------------------------------------------------------------------
2114C ----------------------------------------------------------------------
2115C ----------------------------------------------------------------------
2116C RCHL UPDATE - END OF NEW FUNCTIONS
2117C ----------------------------------------------------------------------
2118C ----------------------------------------------------------------------
2119C ----------------------------------------------------------------------
2120
2121
2122 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
2123C ----------------------------------------------------------------------
2124* THIS SIMULATES TAU DECAY IN TAU REST FRAME
2125* INTO NU A1, NEXT A1 DECAYS INTO RHO PI AND FINALLY RHO INTO PI PI.
2126* OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
2127* PAA A1
2128* PIM1 PION MINUS (OR PI0) 1 (FOR TAU MINUS)
2129* PIM2 PION MINUS (OR PI0) 2
2130* PIPL PION PLUS (OR PI-)
2131* (PIPL,PIM1) FORM A RHO
2132C ----------------------------------------------------------------------
2133 COMMON / inout / inut,iout
2134 REAL POL(4),HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4),RN(1)
2135 DATA iwarm/0/
2136C
2137 IF(mode.EQ.-1) THEN
2138C ===================
2139 iwarm=1
2140 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
2141CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
2142C
2143 ELSEIF(mode.EQ. 0) THEN
2144* =======================
2145 300 CONTINUE
2146 IF(iwarm.EQ.0) GOTO 902
2147 CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
2148 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2149CC CALL HFILL(816,WT)
2150 CALL ranmar(rn,1)
2151 IF(rn(1).GT.wt) GOTO 300
2152C
2153 ELSEIF(mode.EQ. 1) THEN
2154* =======================
2155 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
2156CC CALL HPRINT(816)
2157 ENDIF
2158C =====
2159 RETURN
2160 902 WRITE(iout, 9020)
2161 9020 FORMAT(' ----- DEXAA: LACK OF INITIALISATION')
2162 stop
2163 END
2164 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
2165C ----------------------------------------------------------------------
2166* A1 DECAY UNWEIGHTED EVENTS
2167C ----------------------------------------------------------------------
2168 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2169 * ,ampiz,ampi,amro,gamro,ama1,gama1
2170 * ,amk,amkz,amkst,gamkst
2171C
2172 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2173 * ,ampiz,ampi,amro,gamro,ama1,gama1
2174 * ,amk,amkz,amkst,gamkst
2175 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2176 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2177 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2178 real*4 gampmc ,gamper
2179 COMMON / inout / inut,iout
2180 REAL HHV(4)
2181 REAL HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4)
2182 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
2183 real*4 rrr(3)
2184 real*8 swt, sswt
2185 DATA pi /3.141592653589793238462643/
2186 DATA iwarm/0/
2187C
2188 IF(mode.EQ.-1) THEN
2189C ===================
2190 iwarm=1
2191 nevraw=0
2192 nevacc=0
2193 nevovr=0
2194 swt=0
2195 sswt=0
2196 wtmax=1e-20
2197 DO 15 i=1,500
2198 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
2199 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
220015 CONTINUE
2201CC CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADMAA $',100,0,2)
2202C
2203 ELSEIF(mode.EQ. 0) THEN
2204C =======================
2205300 CONTINUE
2206 IF(iwarm.EQ.0) GOTO 902
2207 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
2208CC CALL HFILL(801,WT/WTMAX)
2209 nevraw=nevraw+1
2210 swt=swt+wt
2211ccM.S.>>>>>>
2212cc SSWT=SSWT+WT**2
2213 sswt=sswt+dble(wt)**2
2214ccM.S.<<<<<<
2215 CALL ranmar(rrr,3)
2216 rn=rrr(1)
2217 IF(wt.GT.wtmax) nevovr=nevovr+1
2218 IF(rn*wtmax.GT.wt) GOTO 300
2219C ROTATIONS TO BASIC TAU REST FRAME
2220 costhe=-1.+2.*rrr(2)
2221 thet=acos(costhe)
2222 phi =2*pi*rrr(3)
2223 CALL rotpol(thet,phi,pnu)
2224 CALL rotpol(thet,phi,paa)
2225 CALL rotpol(thet,phi,pim1)
2226 CALL rotpol(thet,phi,pim2)
2227 CALL rotpol(thet,phi,pipl)
2228 CALL rotpol(thet,phi,hv)
2229 DO 44 i=1,3
2230 44 hhv(i)=-isgn*hv(i)
2231 nevacc=nevacc+1
2232C
2233 ELSEIF(mode.EQ. 1) THEN
2234C =======================
2235 IF(nevraw.EQ.0) RETURN
2236 pargam=swt/float(nevraw+1)
2237 error=0
2238 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2239 rat=pargam/gamel
2240 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2241CC CALL HPRINT(801)
2242 gampmc(5)=rat
2243 gamper(5)=error
2244CAM NEVDEC(5)=NEVACC
2245 ENDIF
2246C =====
2247 RETURN
2248 7003 FORMAT(///1x,15(5h*****)
2249 $ /,' *', 25x,'******** DADMAA INITIALISATION ********',9x,1h*
2250 $ /,' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2251 $ /,1x,15(5h*****)/)
2252 7010 FORMAT(///1x,15(5h*****)
2253 $ /,' *', 25x,'******** DADMAA FINAL REPORT ******** ',9x,1h*
2254 $ /,' *',i20 ,5x,'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
2255 $ /,' *',i20 ,5x,'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
2256 $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2257 $ /,' *',e20.5,5x,'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
2258 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2259 $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2260 $ /,1x,15(5h*****)/)
2261 902 WRITE(iout, 9020)
2262 9020 FORMAT(' ----- DADMAA: LACK OF INITIALISATION')
2263 stop
2264 END
2265 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
2266C ----------------------------------------------------------------------
2267* IT SIMULATES A1 DECAY IN TAU REST FRAME WITH
2268* Z-AXIS ALONG A1 MOMENTUM
2269C ----------------------------------------------------------------------
2270 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2271 * ,ampiz,ampi,amro,gamro,ama1,gama1
2272 * ,amk,amkz,amkst,gamkst
2273C
2274 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2275 * ,ampiz,ampi,amro,gamro,ama1,gama1
2276 * ,amk,amkz,amkst,gamkst
2277 COMMON / taukle / bra1,brk0,brk0b,brks
2278 real*4 bra1,brk0,brk0b,brks
2279 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
2280
2281
2282 real*4 rrr(1)
2283C MATRIX ELEMENT NUMBER:
2284 mnum=0
2285C TYPE OF THE GENERATION:
2286 keyt=1
2287 CALL ranmar(rrr,1)
2288 rmod=rrr(1)
2289 IF (rmod.LT.bra1) THEN
2290 jaa=1
2291 amp1=ampi
2292 amp2=ampi
2293 amp3=ampi
2294 ELSE
2295 jaa=2
2296 amp1=ampiz
2297 amp2=ampiz
2298 amp3=ampi
2299 ENDIF
2300 CALL ch3piset(jaa) ! information on sub-chanel passed for further use.
2301 call
2302 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2303 END
2304 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
2305C ----------------------------------------------------------------------
2306C TAU DECAY INTO KAON AND TAU-NEUTRINO
2307C IN TAU REST FRAME
2308C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
2309C PKK KAON CHARGED
2310C ----------------------------------------------------------------------
2311 REAL POL(4),HV(4),PNU(4),PKK(4),RN(1)
2312C
2313 IF(mode.EQ.-1) THEN
2314C ===================
2315 CALL dadmkk(-1,isgn,hv,pkk,pnu)
2316CC CALL HBOOK1(815,'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
2317C
2318 ELSEIF(mode.EQ. 0) THEN
2319C =======================
2320300 CONTINUE
2321 CALL dadmkk( 0,isgn,hv,pkk,pnu)
2322 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2323CC CALL HFILL(815,WT)
2324 CALL ranmar(rn,1)
2325 IF(rn(1).GT.wt) GOTO 300
2326C
2327 ELSEIF(mode.EQ. 1) THEN
2328C =======================
2329 CALL dadmkk( 1,isgn,hv,pkk,pnu)
2330CC CALL HPRINT(815)
2331 ENDIF
2332C =====
2333 RETURN
2334 END
2335 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
2336C ----------------------------------------------------------------------
2337C FZ
2338 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2339 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2340 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2341 * ,ampiz,ampi,amro,gamro,ama1,gama1
2342 * ,amk,amkz,amkst,gamkst
2343C
2344 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2345 * ,ampiz,ampi,amro,gamro,ama1,gama1
2346 * ,amk,amkz,amkst,gamkst
2347 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2348 real*4 gampmc ,gamper
2349 COMMON / inout / inut,iout
2350 REAL PKK(4),PNU(4),HV(4)
2351 DATA pi /3.141592653589793238462643/
2352C
2353 IF(mode.EQ.-1) THEN
2354C ===================
2355 nevtot=0
2356 ELSEIF(mode.EQ. 0) THEN
2357C =======================
2358 nevtot=nevtot+1
2359 ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
2360 enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
2361 xkk= sqrt(ekk**2-amk**2)
2362C K MOMENTUM
2363 CALL sphera(xkk,pkk)
2364 pkk(4)=ekk
2365C TAU-NEUTRINO MOMENTUM
2366 DO 30 i=1,3
236730 pnu(i)=-pkk(i)
2368 pnu(4)=enu
2369 pxq=amtau*ekk
2370 pxn=amtau*enu
2371 qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
2372 brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
2373 & +(gv**2-ga**2)*amtau*amnuta*amk**2
2374 DO 40 i=1,3
237540 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
2376 hv(4)=1
2377C
2378 ELSEIF(mode.EQ. 1) THEN
2379C =======================
2380 IF(nevtot.EQ.0) RETURN
2381 fkk=0.0354
2382CFZ THERE WAS BRAK/AMTAU**4 BEFORE
2383C GAMM=(GFERMI*FKK)**2/(16.*PI)*AMTAU**3*
2384C * (BRAK/AMTAU**4)**2
2385CZW 7.02.93 here was an error affecting non standard model
2386C configurations only
2387 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
2388 $ (brak/amtau**4)*
2389 $ sqrt((amtau**2-amk**2-amnuta**2)**2
2390 $ -4*amk**2*amnuta**2 )/amtau**2
2391 error=0
2392
2393 error=0
2394 rat=gamm/gamel
2395 WRITE(iout, 7010) nevtot,gamm,rat,error
2396 gampmc(6)=rat
2397 gamper(6)=error
2398CAM NEVDEC(6)=NEVTOT
2399 ENDIF
2400C =====
2401 RETURN
2402 7010 FORMAT(///1x,15(5h*****)
2403 $ /,' *', 25x,'******** DADMKK FINAL REPORT ********',9x,1h*
2404 $ /,' *',i20 ,5x,'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
2405 $ /,' *',e20.5,5x,'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
2406 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2407 $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
2408 $ /,1x,15(5h*****)/)
2409 END
2410 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2411C ----------------------------------------------------------------------
2412C THIS SIMULATES TAU DECAY IN TAU REST FRAME
2413C INTO NU K*, THEN K* DECAYS INTO PI0,K+-(JKST=20)
2414C OR PI+-,K0(JKST=10).
2415C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
2416C PKS K* CHARGED
2417C PK0 K ZERO
2418C PKC K CHARGED
2419C PIC PION CHARGED
2420C PIZ PION ZERO
2421C ----------------------------------------------------------------------
2422 COMMON / inout / inut,iout
2423 REAL POL(4),HV(4),PKS(4),PNU(4),PKK(4),PPI(4),RN(1)
2424 DATA iwarm/0/
2425C
2426 IF(mode.EQ.-1) THEN
2427C ===================
2428 iwarm=1
2429CFZ INITIALISATION DONE WITH THE GHARGED PION NEUTRAL KAON MODE(JKST=10
2430 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2431CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXKS $',100,0,2)
2432CC CALL HBOOK1(916,'ABS2 OF HV IN ROUTINE DEXKS $',100,0,2)
2433C
2434 ELSEIF(mode.EQ. 0) THEN
2435C =======================
2436300 CONTINUE
2437 IF(iwarm.EQ.0) GOTO 902
2438 CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2439 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2440CC CALL HFILL(816,WT)
2441CC XHELP=HV(1)**2+HV(2)**2+HV(3)**2
2442CC CALL HFILL(916,XHELP)
2443 CALL ranmar(rn,1)
2444 IF(rn(1).GT.wt) GOTO 300
2445C
2446 ELSEIF(mode.EQ. 1) THEN
2447C ======================================
2448 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2449CC CALL HPRINT(816)
2450CC CALL HPRINT(916)
2451 ENDIF
2452C =====
2453 RETURN
2454 902 WRITE(iout, 9020)
2455 9020 FORMAT(' ----- DEXKS: LACK OF INITIALISATION')
2456 stop
2457 END
2458 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2459C ----------------------------------------------------------------------
2460 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2461 * ,ampiz,ampi,amro,gamro,ama1,gama1
2462 * ,amk,amkz,amkst,gamkst
2463C
2464 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2465 * ,ampiz,ampi,amro,gamro,ama1,gama1
2466 * ,amk,amkz,amkst,gamkst
2467 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2468 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2469 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2470 real*4 gampmc ,gamper
2471 COMMON / taukle / bra1,brk0,brk0b,brks
2472 real*4 bra1,brk0,brk0b,brks
2473 COMMON / inout / inut,iout
2474 REAL HHV(4)
2475 REAL HV(4),PKS(4),PNU(4),PKK(4),PPI(4)
2476 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
2477 real*4 rrr(3),rmod(1)
2478 real*8 swt, sswt
2479 DATA pi /3.141592653589793238462643/
2480 DATA iwarm/0/
2481C
2482 IF(mode.EQ.-1) THEN
2483C ===================
2484 iwarm=1
2485 nevraw=0
2486 nevacc=0
2487 nevovr=0
2488 swt=0
2489 sswt=0
2490 wtmax=1e-20
2491 DO 15 i=1,5000
2492C THE INITIALISATION IS DONE WITH THE 66.7% MODE
2493 jkst=10
2494 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2495 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
249615 CONTINUE
2497CC CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADMKS $',100,0,2)
2498CC PRINT 7003,WTMAX
2499CC CALL HBOOK1(112,'-------- K* MASS -------- $',100,0.,2.)
2500 ELSEIF(mode.EQ. 0) THEN
2501C =====================================
2502 IF(iwarm.EQ.0) GOTO 902
2503C HERE WE CHOOSE RANDOMLY BETWEEN K0 PI+_ (66.7%)
2504C AND K+_ PI0 (33.3%)
2505 dec1=brks
2506400 CONTINUE
2507 CALL ranmar(rmod,1)
2508 IF(rmod(1).LT.dec1) THEN
2509 jkst=10
2510 ELSE
2511 jkst=20
2512 ENDIF
2513 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2514 CALL ranmar(rrr,3)
2515 rn=rrr(1)
2516 IF(wt.GT.wtmax) nevovr=nevovr+1
2517 nevraw=nevraw+1
2518 swt=swt+wt
2519 sswt=sswt+wt**2
2520 IF(rn*wtmax.GT.wt) GOTO 400
2521C ROTATIONS TO BASIC TAU REST FRAME
2522 costhe=-1.+2.*rrr(2)
2523 thet=acos(costhe)
2524 phi =2*pi*rrr(3)
2525 CALL rotor2(thet,pnu,pnu)
2526 CALL rotor3( phi,pnu,pnu)
2527 CALL rotor2(thet,pks,pks)
2528 CALL rotor3( phi,pks,pks)
2529 CALL rotor2(thet,pkk,pkk)
2530 CALL rotor3(phi,pkk,pkk)
2531 CALL rotor2(thet,ppi,ppi)
2532 CALL rotor3( phi,ppi,ppi)
2533 CALL rotor2(thet,hv,hv)
2534 CALL rotor3( phi,hv,hv)
2535 DO 44 i=1,3
2536 44 hhv(i)=-isgn*hv(i)
2537 nevacc=nevacc+1
2538C
2539 ELSEIF(mode.EQ. 1) THEN
2540C =======================
2541 IF(nevraw.EQ.0) RETURN
2542 pargam=swt/float(nevraw+1)
2543 error=0
2544 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2545 rat=pargam/gamel
2546 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2547CC CALL HPRINT(801)
2548 gampmc(7)=rat
2549 gamper(7)=error
2550CAM NEVDEC(7)=NEVACC
2551 ENDIF
2552C =====
2553 RETURN
2554 7003 FORMAT(///1x,15(5h*****)
2555 $ /,' *', 25x,'******** DADMKS INITIALISATION ********',9x,1h*
2556 $ /,' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2557 $ /,1x,15(5h*****)/)
2558 7010 FORMAT(///1x,15(5h*****)
2559 $ /,' *', 25x,'******** DADMKS FINAL REPORT ********',9x,1h*
2560 $ /,' *',i20 ,5x,'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2561 $ /,' *',i20 ,5x,'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2562 $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2563 $ /,' *',e20.5,5x,'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2564 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2565 $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2566 $ /,1x,15(5h*****)/)
2567 902 WRITE(iout, 9020)
2568 9020 FORMAT(' ----- DADMKS: LACK OF INITIALISATION')
2569 stop
2570 END
2571 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2572C ----------------------------------------------------------------------
2573C IT SIMULATES KAON* DECAY IN TAU REST FRAME WITH
2574C Z-AXIS ALONG KAON* MOMENTUM
2575C JKST=10 FOR K* --->K0 + PI+-
2576C JKST=20 FOR K* --->K+- + PI0
2577C ----------------------------------------------------------------------
2578 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2579 * ,ampiz,ampi,amro,gamro,ama1,gama1
2580 * ,amk,amkz,amkst,gamkst
2581C
2582 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2583 * ,ampiz,ampi,amro,gamro,ama1,gama1
2584 * ,amk,amkz,amkst,gamkst
2585 REAL HV(4),PT(4),PN(4),PKS(4),PKK(4),PPI(4),QQ(4),RR1(1),RR2(1)
2586 DATA pi /3.141592653589793238462643/
2587C
2588 DATA icont /0/
2589C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
2590 phspac=1./2**11/pi**5
2591C TAU MOMENTUM
2592 pt(1)=0.
2593 pt(2)=0.
2594 pt(3)=0.
2595 pt(4)=amtau
2596 CALL ranmar(rr1,1)
2597C HERE BEGIN THE K0,PI+_ DECAY
2598 IF(jkst.EQ.10)THEN
2599C ==================
2600C MASS OF (REAL/VIRTUAL) K*
2601 ams1=(ampi+amkz)**2
2602 ams2=(amtau-amnuta)**2
2603 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2604 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2605 CALL ranmar(rr2,1)
2606 prob1=0.2
2607 IF (rr2(1).LT.prob1) THEN
2608C FLAT PHASE SPACE
2609 amx2=ams1+ rr1(1)*(ams2-ams1)
2610 amx=sqrt(amx2)
2611 ELSE
2612C PHASE SPACE WITH SAMPLING FOR K* RESONANCE
2613 alp=alp1+rr1(1)*(alp2-alp1)
2614 amx2=amkst**2+amkst*gamkst*tan(alp)
2615 amx=sqrt(amx2)
2616 ENDIF
2617C merging of the two channels
2618 phspac1=(ams2-ams1)
2619 phspac2=((amx2-amkst**2)**2+(amkst*gamkst)**2)
2620 & /(amkst*gamkst)
2621 phspac2=phspac2*(alp2-alp1)
2622 a1=0.0
2623 a2=0.0
2624 IF (phspac1.NE.0.0) a1=prob1 /phspac1
2625 IF (phspac2.NE.0.0) a2=(1-prob1)/phspac2
2626
2627
2628
2629 IF (a1+a2.NE.0.0) THEN
2630 phspac=phspac/(a1+a2)
2631 ELSE
2632 phspac=0
2633 ENDIF
2634
2635C
2636C TAU-NEUTRINO MOMENTUM
2637 pn(1)=0
2638 pn(2)=0
2639 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2640 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
2641C
2642C K* MOMENTUM
2643 pks(1)=0
2644 pks(2)=0
2645 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2646 pks(3)=-pn(3)
2647 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2648C
2649CAM
2650 enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2651 pppi=sqrt(abs(enpi-ampi)*(enpi+ampi))
2652 phspac=phspac*(4*pi)*(2*pppi/amx)
2653C CHARGED PI MOMENTUM IN KAON* REST FRAME
2654 CALL sphera(pppi,ppi)
2655 ppi(4)=enpi
2656C NEUTRAL KAON MOMENTUM IN K* REST FRAME
2657 DO 20 i=1,3
265820 pkk(i)=-ppi(i)
2659 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2660 exe=(pks(4)+pks(3))/amx
2661C PION AND K BOOSTED FROM K* REST FRAME TO TAU REST FRAME
2662 CALL bostr3(exe,ppi,ppi)
2663 CALL bostr3(exe,pkk,pkk)
2664
2665 CALL dam2pi(1,pt,pn,ppi,pkk,amplit,hv)
2666 dgamt=1/(2.*amtau)*amplit*phspac
2667C
2668C HERE BEGIN THE K+-,PI0 DECAY
2669 ELSEIF(jkst.EQ.20)THEN
2670C ======================
2671C MASS OF (REAL/VIRTUAL) K*
2672 ams1=(ampiz+amk)**2
2673 ams2=(amtau-amnuta)**2
2674 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2675 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2676
2677 CALL ranmar(rr2,1)
2678 prob1=0.2
2679 IF (rr2(1).LT.prob1) THEN
2680C FLAT PHASE SPACE
2681 amx2=ams1+ rr1(1)*(ams2-ams1)
2682 amx=sqrt(amx2)
2683 ELSE
2684C PHASE SPACE WITH SAMPLING FOR K* RESONANCE
2685 alp=alp1+rr1(1)*(alp2-alp1)
2686 amx2=amkst**2+amkst*gamkst*tan(alp)
2687 amx=sqrt(amx2)
2688 ENDIF
2689C merging of the two channels
2690 phspac1=(ams2-ams1)
2691 phspac2=((amx2-amkst**2)**2+(amkst*gamkst)**2)
2692 & /(amkst*gamkst)
2693 phspac2=phspac2*(alp2-alp1)
2694 a1=0.0
2695 a2=0.0
2696 IF (phspac1.NE.0.0) a1=prob1 /phspac1
2697 IF (phspac2.NE.0.0) a2=(1-prob1)/phspac2
2698
2699 IF (a1+a2.NE.0.0) THEN
2700 phspac=phspac/(a1+a2)
2701 ELSE
2702 phspac=0
2703 ENDIF
2704
2705C
2706C TAU-NEUTRINO MOMENTUM
2707 pn(1)=0
2708 pn(2)=0
2709 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2710 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
2711C KAON* MOMENTUM
2712 pks(1)=0
2713 pks(2)=0
2714 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2715 pks(3)=-pn(3)
2716 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2717C
2718CAM
2719 enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2720 pppi=sqrt(abs(enpi-ampiz)*(enpi+ampiz))
2721 phspac=phspac*(4*pi)*(2*pppi/amx)
2722C NEUTRAL PI MOMENTUM IN K* REST FRAME
2723 CALL sphera(pppi,ppi)
2724 ppi(4)=enpi
2725C CHARGED KAON MOMENTUM IN K* REST FRAME
2726 DO 50 i=1,3
272750 pkk(i)=-ppi(i)
2728 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2729 exe=(pks(4)+pks(3))/amx
2730C PION AND K BOOSTED FROM K* REST FRAME TO TAU REST FRAME
2731 CALL bostr3(exe,ppi,ppi)
2732 CALL bostr3(exe,pkk,pkk)
2733
2734 CALL dam2pi(2,pt,pn,pkk,ppi,amplit,hv)
2735 dgamt=1/(2.*amtau)*amplit*phspac
2736
2737 ENDIF
2738
2739 RETURN
2740 END
2741
2742
2743
2744 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2745C ----------------------------------------------------------------------
2746C IT SIMULATES MULTIPI DECAY IN TAU REST FRAME WITH
2747C Z-AXIS OPPOSITE TO NEUTRINO MOMENTUM
2748C ----------------------------------------------------------------------
2749 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2750 * ,ampiz,ampi,amro,gamro,ama1,gama1
2751 * ,amk,amkz,amkst,gamkst
2752C
2753 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2754 * ,ampiz,ampi,amro,gamro,ama1,gama1
2755 * ,amk,amkz,amkst,gamkst
2756 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2757 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2758 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2759 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2760 & ,names
2761 CHARACTER NAMES(NMODE)*31
2762 real*8 wetmax(20)
2763C
2764 real*8 pn(4),pr(4),ppi(4,9),hv(4)
2765 real*4 pnx(4),prx(4),ppix(4,9),hvx(4)
2766 real*8 pv(5,9),pt(4),ue(3),be(3)
2767 real*8 pawt,amx,ams1,ams2,pa,phs,phsmax,pmin,pmax
2768!!! M.S. to fix underflow >>>
2769 real*8 phspac
2770!!! M.S. to fix underflow <<<
2771 real*8 gam,bep,phi,a,b,c
2772 real*8 ampik
2773 real*4 rrr(9),rrx(2),rn(1),rr2(1)
2774C
2775 DATA pi /3.141592653589793238462643/
2776 DATA wetmax /20*1d-15/
2777C
2778CC-- PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A)
2779C
2780 pawt(a,b,c)=
2781 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2782C
2783 ampik(i,j)=dcdmas(idffin(i,j))
2784C
2785C
2786 IF ((jnpi.LE.0).OR.jnpi.GT.20) THEN
2787 WRITE(6,*) 'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2788 stop
2789 ENDIF
2790
2791C TAU MOMENTUM
2792 pt(1)=0.
2793 pt(2)=0.
2794 pt(3)=0.
2795 pt(4)=amtau
2796C
2797 500 CONTINUE
2798C MASS OF VIRTUAL W
2799 nd=mulpik(jnpi)
2800 ps=0.
2801 phspac = 1./2.**5 /pi**2
2802 DO 4 i=1,nd
28034 ps =ps+ampik(i,jnpi)
2804 CALL ranmar(rr2,1)
2805 ams1=ps**2
2806 ams2=(amtau-amnuta)**2
2807C
2808C
2809 amx2=ams1+ rr2(1)*(ams2-ams1)
2810 amx =sqrt(amx2)
2811 amw =amx
2812 phspac=phspac * (ams2-ams1)
2813C
2814C TAU-NEUTRINO MOMENTUM
2815 pn(1)=0
2816 pn(2)=0
2817 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2818 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
2819C W MOMENTUM
2820 pr(1)=0
2821 pr(2)=0
2822 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2823 pr(3)=-pn(3)
2824 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2825C
2826C AMPLITUDE (cf YS.Tsai Phys.Rev.D4,2821(1971)
2827C or F.Gilman SH.Rhie Phys.Rev.D31,1066(1985)
2828C
2829 pxq=amtau*pr(4)
2830 pxn=amtau*pn(4)
2831 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2832C HERE WAS AN ERROR. 20.10.91 (ZW)
2833C BRAK=2*(GV**2+GA**2)*(2*PXQ*PXN+AMX2*QXN)
2834 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2835 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2836CAM Assume neutrino mass=0. and sum over final polarisation
2837C BRAK= 2*(AMTAU**2-AMX2) * (AMTAU**2+2.*AMX2)
2838 amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,jnpi)
2839 dgamt=1./(2.*amtau)*amplit*phspac
2840C
2841C ISOTROPIC W DECAY IN W REST FRAME
2842 phsmax = 1.
2843 DO 200 i=1,4
2844 200 pv(i,1)=pr(i)
2845 pv(5,1)=amw
2846 pv(5,nd)=ampik(nd,jnpi)
2847C COMPUTE MAX. PHASE SPACE FACTOR
2848 pmax=amw-ps+ampik(nd,jnpi)
2849 pmin=.0
2850 DO 220 il=nd-1,1,-1
2851 pmax=pmax+ampik(il,jnpi)
2852 pmin=pmin+ampik(il+1,jnpi)
2853 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2854
2855C --- 2.02.94 ZW 9 lines
2856 amx=amw
2857 DO 222 il=1,nd-2
2858 ams1=.0
2859 DO 223 jl=il+1,nd
2860 223 ams1=ams1+ampik(jl,jnpi)
2861 ams1=ams1**2
2862 amx =(amx-ampik(il,jnpi))
2863 ams2=(amx)**2
2864 phsmax=phsmax * (ams2-ams1)
2865 222 CONTINUE
2866 ncont=0
2867 100 CONTINUE
2868 ncont=ncont+1
2869CAM GENERATE ND-2 EFFECTIVE MASSES
2870 phs=1.d0
2871 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2872 amx=amw
2873 CALL ranmar(rrr,nd-2)
2874 DO 230 il=1,nd-2
2875 ams1=.0d0
2876 DO 231 jl=il+1,nd
2877 231 ams1=ams1+ampik(jl,jnpi)
2878 ams1=ams1**2
2879 ams2=(amx-ampik(il,jnpi))**2
2880 rr1=rrr(il)
2881 amx2=ams1+ rr1*(ams2-ams1)
2882 amx=sqrt(amx2)
2883 pv(5,il+1)=amx
2884 phspac=phspac * (ams2-ams1)
2885C --- 2.02.94 ZW 1 line
2886 phs=phs* (ams2-ams1)
2887 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2888 phs =phs *pa/pv(5,il)
2889 230 CONTINUE
2890 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2891 phs =phs *pa/pv(5,nd-1)
2892 CALL ranmar(rn,1)
2893
2894 IF(phsmax.NE.0.0) THEN ! TP 5.10.2011 due to rounding errs.
2895 ! PHSMAX may be zero, protect div. by it
2896 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2897 ELSE
2898 wetmax(jnpi)=1.2d0*wetmax(jnpi)/1.2d0
2899 ENDIF
2900
2901 IF (ncont.EQ.500 000) THEN
2902 xnpi=0.0
2903 DO kk=1,nd
2904 xnpi=xnpi+ampik(kk,jnpi)
2905 ENDDO
2906 WRITE(6,*) 'ROUNDING INSTABILITY IN DPHNPI ?'
2907 WRITE(6,*) 'AMW=',amw,'XNPI=',xnpi
2908 WRITE(6,*) 'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
2909 WRITE(6,*) 'PHS=',phs,'PHSMAX=',phsmax
2910 GOTO 500
2911 ENDIF
2912 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs) GO TO 100
2913C...PERFORM SUCCESSIVE TWO-PARTICLE DECAYS IN RESPECTIVE CM FRAME
2914 280 DO 300 il=1,nd-1
2915 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2916 CALL ranmar(rrx,2)
2917 ue(3)=2.*rrx(1)-1.
2918 phi=2.*pi*rrx(2)
2919 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2920 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2921 DO 290 j=1,3
2922 ppi(j,il)=pa*ue(j)
2923 290 pv(j,il+1)=-pa*ue(j)
2924 ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2925 pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2926 phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2927 300 CONTINUE
2928C...LORENTZ TRANSFORM DECAY PRODUCTS TO TAU FRAME
2929 DO 310 j=1,4
2930 310 ppi(j,nd)=pv(j,nd)
2931 DO 340 il=nd-1,1,-1
2932 DO 320 j=1,3
2933 320 be(j)=pv(j,il)/pv(4,il)
2934 gam=pv(4,il)/pv(5,il)
2935 DO 340 i=il,nd
2936 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2937 DO 330 j=1,3
2938 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2939 ppi(4,i)=gam*(ppi(4,i)+bep)
2940 340 CONTINUE
2941C
2942 hv(4)=1.
2943 hv(3)=0.
2944 hv(2)=0.
2945 hv(1)=0.
2946 DO k=1,4
2947 pnx(k)=pn(k)
2948 prx(k)=pr(k)
2949 hvx(k)=hv(k)
2950 DO l=1,nd
2951 ppix(k,l)=ppi(k,l)
2952 ENDDO
2953 ENDDO
2954 RETURN
2955 END
2956 FUNCTION sigee(Q2,JNP)
2957C ----------------------------------------------------------------------
2958C e+e- cross section in the (1.GEV2,AMTAU**2) region
2959C normalised to sig0 = 4/3 pi alfa2
2960C used in matrix element for multipion tau decays
2961C cf YS.Tsai Phys.Rev D4 ,2821(1971)
2962C F.Gilman et al Phys.Rev D17,1846(1978)
2963C C.Kiesling, to be pub. in High Energy e+e- Physics (1988)
2964C DATSIG(*,1) = e+e- -> pi+pi-2pi0
2965C DATSIG(*,2) = e+e- -> 2pi+2pi-
2966C DATSIG(*,3) = 5-pion contribution (a la TN.Pham et al)
2967C (Phys Lett 78B,623(1978)
2968C DATSIG(*,5) = e+e- -> 6pi
2969C
2970C 4- and 6-pion cross sections from data
2971C 5-pion contribution related to 4-pion cross section
2972C
2973C Called by DPHNPI
2974C ----------------------------------------------------------------------
2975 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2976 * ,ampiz,ampi,amro,gamro,ama1,gama1
2977 * ,amk,amkz,amkst,gamkst
2978C
2979 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2980 * ,ampiz,ampi,amro,gamro,ama1,gama1
2981 * ,amk,amkz,amkst,gamkst
2982 real*4 datsig(17,6)
2983C
2984 DATA datsig/
2985 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2986 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2987 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2988 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2989 5 17*.0,
2990 6 17*.0,
2991 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2992 8 17*.0/
2993 DATA sig0 / 86.8 /
2994 DATA pi /3.141592653589793238462643/
2995 DATA init / 0 /
2996C
2997 jnpi=jnp
2998 IF(jnp.EQ.4) jnpi=3
2999 IF(jnp.EQ.3) jnpi=4
3000 IF(init.EQ.0) THEN
3001 init=1
3002C AJWMOD: initialize if called from outside QQ:
3003C IF (AMPI.LT.0.139) AMPI = 0.1395675
3004 ampi2=ampi**2
3005 fpi = .943*ampi
3006 DO 100 i=1,17
3007 datsig(i,2) = datsig(i,2)/2.
3008 datsig(i,1) = datsig(i,1) + datsig(i,2)
3009 s = 1.025+(i-1)*.05
3010 fact=0.
3011 s2=s**2
3012 DO 200 j=1,17
3013 t= 1.025+(j-1)*.05
3014 IF(t . gt. s-ampi ) GO TO 201
3015 t2=t**2
3016 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
3017 fact = fact * (datsig(j,1)+datsig(j+1,1))
3018 200 datsig(i,3) = datsig(i,3) + fact
3019 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
3020 datsig(i,4) = datsig(i,3)
3021 datsig(i,6) = datsig(i,5)
3022 100 CONTINUE
3023C WRITE(6,1000) DATSIG
3024 1000 FORMAT(///1x,' EE SIGMA USED IN MULTIPI DECAYS'/
3025 % (17f7.2/))
3026 ENDIF
3027 q=sqrt(q2)
3028 qmin=1.
3029 IF(q.LT.qmin) THEN
3030 sigee=datsig(1,jnpi)+
3031 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
3032 ELSEIF(q.LT.1.8) THEN
3033 DO 1 i=1,16
3034 qmax = qmin + .05
3035 IF(q.LT.qmax) GO TO 2
3036 qmin = qmin + .05
3037 1 CONTINUE
3038 2 sigee=datsig(i,jnpi)+
3039 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
3040 ELSEIF(q.GT.1.8) THEN
3041 sigee=datsig(17,jnpi)+
3042 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
3043 ENDIF
3044 IF(sigee.LT..0) sigee=0.
3045C
3046 sigee = sigee/(6.*pi**2*sig0)
3047C
3048 RETURN
3049 END
3050
3051 FUNCTION sigold(Q2,JNPI)
3052C ----------------------------------------------------------------------
3053C e+e- cross section in the (1.GEV2,AMTAU**2) region
3054C normalised to sig0 = 4/3 pi alfa2
3055C used in matrix element for multipion tau decays
3056C cf YS.Tsai Phys.Rev D4 ,2821(1971)
3057C F.Gilman et al Phys.Rev D17,1846(1978)
3058C C.Kiesling, to be pub. in High Energy e+e- Physics (1988)
3059C DATSIG(*,1) = e+e- -> pi+pi-2pi0
3060C DATSIG(*,2) = e+e- -> 2pi+2pi-
3061C DATSIG(*,3) = 5-pion contribution (a la TN.Pham et al)
3062C (Phys Lett 78B,623(1978)
3063C DATSIG(*,4) = e+e- -> 6pi
3064C
3065C 4- and 6-pion cross sections from data
3066C 5-pion contribution related to 4-pion cross section
3067C
3068C Called by DPHNPI
3069C ----------------------------------------------------------------------
3070 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3071 * ,ampiz,ampi,amro,gamro,ama1,gama1
3072 * ,amk,amkz,amkst,gamkst
3073C
3074 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3075 * ,ampiz,ampi,amro,gamro,ama1,gama1
3076 * ,amk,amkz,amkst,gamkst
3077 real*4 datsig(17,4)
3078C
3079 DATA datsig/
3080 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
3081 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
3082 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
3083 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
3084 5 17*.0,
3085 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
3086 DATA sig0 / 86.8 /
3087 DATA pi /3.141592653589793238462643/
3088 DATA init / 0 /
3089C
3090 IF(init.EQ.0) THEN
3091 init=1
3092 ampi2=ampi**2
3093 fpi = .943*ampi
3094 DO 100 i=1,17
3095 datsig(i,2) = datsig(i,2)/2.
3096 datsig(i,1) = datsig(i,1) + datsig(i,2)
3097 s = 1.025+(i-1)*.05
3098 fact=0.
3099 s2=s**2
3100 DO 200 j=1,17
3101 t= 1.025+(j-1)*.05
3102 IF(t . gt. s-ampi ) GO TO 201
3103 t2=t**2
3104 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
3105 fact = fact * (datsig(j,1)+datsig(j+1,1))
3106 200 datsig(i,3) = datsig(i,3) + fact
3107 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
3108 100 CONTINUE
3109C WRITE(6,1000) DATSIG
3110 1000 FORMAT(///1x,' EE SIGMA USED IN MULTIPI DECAYS'/
3111 % (17f7.2/))
3112 ENDIF
3113 q=sqrt(q2)
3114 qmin=1.
3115 IF(q.LT.qmin) THEN
3116 sigee=datsig(1,jnpi)+
3117 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
3118 ELSEIF(q.LT.1.8) THEN
3119 DO 1 i=1,16
3120 qmax = qmin + .05
3121 IF(q.LT.qmax) GO TO 2
3122 qmin = qmin + .05
3123 1 CONTINUE
3124 2 sigee=datsig(i,jnpi)+
3125 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
3126 ELSEIF(q.GT.1.8) THEN
3127 sigee=datsig(17,jnpi)+
3128 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
3129 ENDIF
3130 IF(sigee.LT..0) sigee=0.
3131C
3132 sigee = sigee/(6.*pi**2*sig0)
3133 sigold=sigee
3134C
3135 RETURN
3136 END
3137 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
3138C ----------------------------------------------------------------------
3139* IT SIMULATES THREE PI (K) DECAY IN THE TAU REST FRAME
3140* Z-AXIS ALONG HADRONIC SYSTEM
3141C ----------------------------------------------------------------------
3142 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3143 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3144 & ,names
3145 CHARACTER NAMES(NMODE)*31
3146
3147 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PNPI(4,9)
3148C MATRIX ELEMENT NUMBER:
3149 mnum=jaa
3150C TYPE OF THE GENERATION:
3151 keyt=4
3152 IF(jaa.EQ.7) keyt=3
3153C --- MASSES OF THE DECAY PRODUCTS
3154 amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
3155 amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
3156 amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
3157 call
3158 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
3159 DO i=1,4
3160 pnpi(i,1)=pim1(i)
3161 pnpi(i,2)=pim2(i)
3162 pnpi(i,3)=pipl(i)
3163 ENDDO
3164 END
3165
3166
3167
3168
3169 subroutine
3170 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
3171C ----------------------------------------------------------------------
3172* IT SIMULATES A1 DECAY IN TAU REST FRAME WITH
3173* Z-AXIS ALONG A1 MOMENTUM
3174* it can be also used to generate K K pi and K pi pi tau decays.
3175* INPUT PARAMETERS
3176* KEYT - algorithm controlling switch
3177* 2 - flat phase space PIM1 PIM2 symmetrized statistical factor 1/2
3178* 1 - like 1 but peaked around a1 and rho (two channels) masses.
3179* 3 - peaked around omega, all particles different
3180* other- flat phase space, all particles different
3181* AMP1 - mass of first pi, etc. (1-3)
3182* MNUM - matrix element type
3183* 0 - a1 matrix element
3184* 1-6 - matrix element for K pi pi, K K pi decay modes
3185* 7 - pi- pi0 gamma matrix element
3186C ----------------------------------------------------------------------
3187 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3188 * ,ampiz,ampi,amro,gamro,ama1,gama1
3189 * ,amk,amkz,amkst,gamkst
3190C
3191 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3192 * ,ampiz,ampi,amro,gamro,ama1,gama1
3193 * ,amk,amkz,amkst,gamkst
3194 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3195 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3196 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
3197 REAL PR(4)
3198 real*4 rrr(5)
3199 DATA pi /3.141592653589793238462643/
3200 DATA icont /0/
3201 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3202C AMRO, GAMRO IS ONLY A PARAMETER FOR GETING HIGHT EFFICIENCY
3203C
3204C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
3205C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P)
3206 phspac=1./2**17/pi**8
3207C TAU MOMENTUM
3208 pt(1)=0.
3209 pt(2)=0.
3210 pt(3)=0.
3211 pt(4)=amtau
3212C
3213 CALL ranmar(rrr,5)
3214 rr=rrr(5)
3215C
3216 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
3217 $ amrx,gamrx,amra,gamra,amrb,gamrb)
3218 IF (ichan.EQ.1) THEN
3219 amp1=ampb
3220 amp2=ampa
3221 ELSEIF (ichan.EQ.2) THEN
3222 amp1=ampa
3223 amp2=ampb
3224 ELSE
3225 amp1=ampb
3226 amp2=ampa
3227 ENDIF
3228CAM
3229 rr1=rrr(1)
3230 ams1=(amp1+amp2+amp3)**2
3231 ams2=(amtau-amnuta)**2
3232* PHASE SPACE WITH SAMPLING FOR A1 RESONANCE
3233 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3234 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3235 alp=alp1+rr1*(alp2-alp1)
3236 am3sq =amrx**2+amrx*gamrx*tan(alp)
3237 am3 =sqrt(am3sq)
3238 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3239 phspac=phspac*(alp2-alp1)
3240C MASS OF (REAL/VIRTUAL) RHO -
3241 rr2=rrr(2)
3242 ams1=(amp2+amp3)**2
3243 ams2=(am3-amp1)**2
3244 IF (ichan.LE.2) THEN
3245* PHASE SPACE WITH SAMPLING FOR RHO RESONANCE,
3246 alp1=atan((ams1-amra**2)/amra/gamra)
3247 alp2=atan((ams2-amra**2)/amra/gamra)
3248 alp=alp1+rr2*(alp2-alp1)
3249 am2sq =amra**2+amra*gamra*tan(alp)
3250 am2 =sqrt(am2sq)
3251C --- THIS PART OF THE JACOBIAN WILL BE RECOVERED LATER ---------------
3252C PHSPAC=PHSPAC*(ALP2-ALP1)
3253C PHSPAC=PHSPAC*((AM2SQ-AMRA**2)**2+(AMRA*GAMRA)**2)/(AMRA*GAMRA)
3254C----------------------------------------------------------------------
3255 ELSE
3256* FLAT PHASE SPACE;
3257 am2sq=ams1+ rr2*(ams2-ams1)
3258 am2 =sqrt(am2sq)
3259 phf0=(ams2-ams1)
3260 ENDIF
3261* RHO RESTFRAME, DEFINE PIPL AND PIM1
3262 enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
3263 enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
3264 ppi= enq1**2-amp3**2
3265 pppi=sqrt(abs(enq1**2-amp3**2))
3266C --- this part of jacobian will be recovered later
3267 phf1=(4*pi)*(2*pppi/am2)
3268* PI MINUS MOMENTUM IN RHO REST FRAME
3269 CALL sphera(pppi,pipl)
3270 pipl(4)=enq1
3271* PI0 1 MOMENTUM IN RHO REST FRAME
3272 DO 30 i=1,3
3273 30 pim1(i)=-pipl(i)
3274 pim1(4)=enq2
3275* A1 REST FRAME, DEFINE PIM2
3276* RHO MOMENTUM
3277 pr(1)=0
3278 pr(2)=0
3279 pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
3280 pr(3)= sqrt(abs(pr(4)**2-am2**2))
3281 ppi = pr(4)**2-am2**2
3282* PI0 2 MOMENTUM
3283 pim2(1)=0
3284 pim2(2)=0
3285 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
3286 pim2(3)=-pr(3)
3287 phf2=(4*pi)*(2*pr(3)/am3)
3288* OLD PIONS BOOSTED FROM RHO REST FRAME TO A1 REST FRAME
3289 exe=(pr(4)+pr(3))/am2
3290 CALL bostr3(exe,pipl,pipl)
3291 CALL bostr3(exe,pim1,pim1)
3292 rr3=rrr(3)
3293 rr4=rrr(4)
3294CAM THET =PI*RR3
3295 thet =acos(-1.+2*rr3)
3296 phi = 2*pi*rr4
3297 CALL rotpol(thet,phi,pipl)
3298 CALL rotpol(thet,phi,pim1)
3299 CALL rotpol(thet,phi,pim2)
3300 CALL rotpol(thet,phi,pr)
3301C
3302* NOW TO THE TAU REST FRAME, DEFINE A1 AND NEUTRINO MOMENTA
3303* A1 MOMENTUM
3304 paa(1)=0
3305 paa(2)=0
3306 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
3307 paa(3)= sqrt(abs(paa(4)**2-am3**2))
3308 ppi = paa(4)**2-am3**2
3309 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
3310* TAU-NEUTRINO MOMENTUM
3311 pn(1)=0
3312 pn(2)=0
3313 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
3314 pn(3)=-paa(3)
3315C HERE WE CORRECT FOR THE JACOBIANS OF THE TWO CHAINS
3316C ---FIRST CHANNEL ------- PIM1+PIPL
3317 ams1=(amp2+amp3)**2
3318 ams2=(am3-amp1)**2
3319 alp1=atan((ams1-amra**2)/amra/gamra)
3320 alp2=atan((ams2-amra**2)/amra/gamra)
3321 xpro = (pim1(3)+pipl(3))**2
3322 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
3323 am2sq=-xpro+(pim1(4)+pipl(4))**2
3324C JACOBIAN OF SPEEDING
3325 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3326 ff1 =ff1 *(alp2-alp1)
3327C LAMBDA OF RHO DECAY
3328 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
3329C LAMBDA OF A1 DECAY
3330 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
3331 xjaje=gg1*(ams2-ams1)
3332C ---SECOND CHANNEL ------ PIM2+PIPL
3333 ams1=(amp1+amp3)**2
3334 ams2=(am3-amp2)**2
3335 alp1=atan((ams1-amrb**2)/amrb/gamrb)
3336 alp2=atan((ams2-amrb**2)/amrb/gamrb)
3337 xpro = (pim2(3)+pipl(3))**2
3338 $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
3339 am2sq=-xpro+(pim2(4)+pipl(4))**2
3340 ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
3341 ff2 =ff2 *(alp2-alp1)
3342 gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
3343 gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
3344 xjadw=gg2*(ams2-ams1)
3345C
3346 a1=0.0
3347 a2=0.0
3348 a3=0.0
3349 xjac1=ff1*gg1
3350 xjac2=ff2*gg2
3351 IF (ichan.EQ.2) THEN
3352 xjac3=xjadw
3353 ELSE
3354 xjac3=xjaje
3355 ENDIF
3356 IF (xjac1.NE.0.0) a1=prob1/xjac1
3357 IF (xjac2.NE.0.0) a2=prob2/xjac2
3358 IF (xjac3.NE.0.0) a3=prob3/xjac3
3359C
3360 IF (a1+a2+a3.NE.0.0) THEN
3361 phspac=phspac/(a1+a2+a3)
3362 ELSE
3363 phspac=0.0
3364 ENDIF
3365 IF(ichan.EQ.2) THEN
3366 DO 70 i=1,4
3367 x=pim1(i)
3368 pim1(i)=pim2(i)
3369 70 pim2(i)=x
3370 ENDIF
3371* ALL PIONS BOOSTED FROM A1 REST FRAME TO TAU REST FRAME
3372* Z-AXIS ANTIPARALLEL TO NEUTRINO MOMENTUM
3373 exe=(paa(4)+paa(3))/am3
3374 CALL bostr3(exe,pipl,pipl)
3375 CALL bostr3(exe,pim1,pim1)
3376 CALL bostr3(exe,pim2,pim2)
3377 CALL bostr3(exe,pr,pr)
3378C PARTIAL WIDTH CONSISTS OF PHASE SPACE AND AMPLITUDE
3379 IF (mnum.EQ.8) THEN
3380 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
3381C ELSEIF (MNUM.EQ.0) THEN
3382C CALL DAMPAA(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3383 ELSE
3384 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
3385 ENDIF
3386 IF (keyt.EQ.1.OR.keyt.EQ.2) THEN
3387C THE STATISTICAL FACTOR FOR IDENTICAL PI-S IS CANCELLED WITH
3388C TWO, FOR TWO MODES OF A1 DECAY NAMELLY PI+PI-PI- AND PI-PI0PI0
3389 phspac=phspac*2.0
3390 phspac=phspac/2.
3391 ENDIF
3392 dgamt=1/(2.*amtau)*amplit*phspac
3393 END
3394 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3395C ----------------------------------------------------------------------
3396* CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
3397* FOR TAU DECAY INTO A1, A1 DECAYS NEXT INTO RHO+PI AND RHO INTO PI+PI.
3398* ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
3399* CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
3400* THE ROUTINE IS WRITEN FOR ZERO NEUTRINO MASS.
3401C
3402C called by : DPHSAA
3403C ----------------------------------------------------------------------
3404 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3405 * ,ampiz,ampi,amro,gamro,ama1,gama1
3406 * ,amk,amkz,amkst,gamkst
3407C
3408 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3409 * ,ampiz,ampi,amro,gamro,ama1,gama1
3410 * ,amk,amkz,amkst,gamkst
3411 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3412 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3413 COMMON /testa1/ keya1
3414 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3415 REAL PAA(4),VEC1(4),VEC2(4)
3416 REAL PIVEC(4),PIAKS(4),HVM(4)
3417 COMPLEX BWIGN,HADCUR(4),FPIK
3418 DATA icont /1/
3419C
3420* F CONSTANTS FOR A1, A1-RHO-PI, AND RHO-PI-PI
3421*
3422 DATA fpi /93.3e-3/
3423* THIS INLINE FUNCT. CALCULATES THE SCALAR PART OF THE PROPAGATOR
3424 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3425C
3426* FOUR MOMENTUM OF A1
3427 DO 10 i=1,4
3428 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3429* MASSES OF A1, AND OF TWO PI-PAIRS WHICH MAY FORM RHO
3430 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3431 xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3432 $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3433 xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3434 $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3435* ELEMENTS OF HADRON CURRENT
3436 prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3437 $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3438 prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3439 $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3440 DO 40 i=1,4
3441 vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3442 40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3443* HADRON CURRENT SATURATED WITH A1 AND RHO RESONANCES
3444 IF (keya1.EQ.1) THEN
3445 fa1=9.87
3446 faropi=1.0
3447 fro2pi=1.0
3448 fnorm=fa1/sqrt(2.)*faropi*fro2pi
3449 DO 45 i=1,4
3450 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3451 $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3452 $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3453 45 CONTINUE
3454 ELSE
3455 fnorm=2.0*sqrt(2.)/3.0/fpi
3456 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3457 DO 46 i=1,4
3458 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3459 $ *(cmplx(vec1(i))*fpik(xmro1)
3460 $ +cmplx(vec2(i))*fpik(xmro2))
3461 46 CONTINUE
3462 ENDIF
3463C
3464* CALCULATE PI-VECTORS: VECTOR AND AXIAL
3465 CALL clvec(hadcur,pn,pivec)
3466 CALL claxi(hadcur,pn,piaks)
3467 CALL clnut(hadcur,brakm,hvm)
3468* SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
3469 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3470 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3471 amplit=(gfermi*ccabib)**2*brak/2.
3472C THE STATISTICAL FACTOR FOR IDENTICAL PI-S WAS CANCELLED WITH
3473C TWO, FOR TWO MODES OF A1 DECAY NAMELLY PI+PI-PI- AND PI-PI0PI0
3474C POLARIMETER VECTOR IN TAU REST FRAME
3475 DO 90 i=1,3
3476 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3477 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3478C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
3479 hv(i)=-hv(i)/brak
3480 90 CONTINUE
3481 END
3482
3483 FUNCTION gfun(QKWA)
3484C ****************************************************************
3485C G-FUNCTION USED TO INRODUCE ENERGY DEPENDENCE IN A1 WIDTH
3486C ****************************************************************
3487 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3488 * ,ampiz,ampi,amro,gamro,ama1,gama1
3489 * ,amk,amkz,amkst,gamkst
3490C
3491 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3492 * ,ampiz,ampi,amro,gamro,ama1,gama1
3493 * ,amk,amkz,amkst,gamkst
3494C
3495 IF (qkwa.LT.(amro+ampi)**2) THEN
3496 gfun=4.1*(qkwa-9*ampiz**2)**3
3497 $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3498 ELSE
3499 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3500 ENDIF
3501 END
3502 COMPLEX FUNCTION bwigs(S,M,G)
3503C **********************************************************
3504C P-WAVE BREIT-WIGNER FOR K*
3505C **********************************************************
3506 REAL S,M,G
3507 REAL PI,PIM,QS,QM,W,GS,MK
3508C AJW: add K*-prim possibility:
3509 REAL PM, PG, PBETA
3510 COMPLEX BW,BWP
3511 DATA init /0/
3512 p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3513 $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3514C ------------ PARAMETERS --------------------
3515 IF (init.EQ.0) THEN
3516 init=1
3517 pi=3.141592654
3518 pim=.139
3519 mk=.493667
3520C AJW: add K*-prim possibility:
3521 pm = pkorb(1,16)
3522 pg = pkorb(2,16)
3523 pbeta = pkorb(3,16)
3524C ------- BREIT-WIGNER -----------------------
3525 ENDIF
3526 qs=p(s,pim**2,mk**2)
3527 qm=p(m**2,pim**2,mk**2)
3528 w=sqrt(s)
3529 gs=g*(m/w)*(qs/qm)**3
3530 bw=m**2/cmplx(m**2-s,-m*gs)
3531 qpm=p(pm**2,pim**2,mk**2)
3532 g1=pg*(pm/w)*(qs/qpm)**3
3533 bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3534 bwigs= (bw+pbeta*bwp)/(1+pbeta)
3535 RETURN
3536 END
3537 COMPLEX FUNCTION bwig(S,M,G)
3538C **********************************************************
3539C P-WAVE BREIT-WIGNER FOR RHO
3540C **********************************************************
3541 REAL S,M,G
3542 REAL PI,PIM,QS,QM,W,GS
3543 DATA init /0/
3544C ------------ PARAMETERS --------------------
3545 IF (init.EQ.0) THEN
3546 init=1
3547 pi=3.141592654
3548 pim=.139
3549C ------- BREIT-WIGNER -----------------------
3550 ENDIF
3551 IF (s.GT.4.*pim**2) THEN
3552 qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3553 qm=sqrt(m**2/4.-pim**2)
3554 w=sqrt(s)
3555 gs=g*(m/w)*(qs/qm)**3
3556 ELSE
3557 gs=0.0
3558 ENDIF
3559 bwig=m**2/cmplx(m**2-s,-m*gs)
3560 RETURN
3561 END
3562 COMPLEX FUNCTION fpik(W)
3563C **********************************************************
3564C PION FORM FACTOR
3565C **********************************************************
3566 COMPLEX BWIG
3567 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
3568 EXTERNAL bwig
3569 DATA init /0/
3570C
3571C ------------ PARAMETERS --------------------
3572 IF (init.EQ.0 ) THEN
3573 init=1
3574 pi=3.141592654
3575 pim=.140
3576 rom=pkorb(1,9)
3577 rog=pkorb(2,9)
3578 rom1=pkorb(1,15)
3579 rog1=pkorb(2,15)
3580 beta1=pkorb(3,15)
3581 ENDIF
3582C -----------------------------------------------
3583 s=w**2
3584 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3585 & /(1+beta1)
3586 RETURN
3587 END
3588 FUNCTION fpirho(W)
3589C **********************************************************
3590C SQUARE OF PION FORM FACTOR
3591C **********************************************************
3592 COMPLEX FPIK
3593 fpirho=cabs(fpik(w))**2
3594 END
3595 SUBROUTINE clvec(HJ,PN,PIV)
3596C ----------------------------------------------------------------------
3597* CALCULATES THE "VECTOR TYPE" PI-VECTOR PIV
3598* NOTE THAT THE NEUTRINO MOM. PN IS ASSUMED TO BE ALONG Z-AXIS
3599C
3600C called by : DAMPAA
3601C ----------------------------------------------------------------------
3602 REAL PIV(4),PN(4)
3603 COMPLEX HJ(4),HN
3604C
3605 hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3606 hh= real(hj(4)*conjg(hj(4))-hj(3)*conjg(hj(3))
3607 $ -hj(2)*conjg(hj(2))-hj(1)*conjg(hj(1)))
3608 DO 10 i=1,4
3609 10 piv(i)=4.*real(hn*conjg(hj(i)))-2.*hh*pn(i)
3610 RETURN
3611 END
3612 SUBROUTINE claxi(HJ,PN,PIA)
3613C ----------------------------------------------------------------------
3614* CALCULATES THE "AXIAL TYPE" PI-VECTOR PIA
3615* NOTE THAT THE NEUTRINO MOM. PN IS ASSUMED TO BE ALONG Z-AXIS
3616C SIGN is chosen +/- for decay of TAU +/- respectively
3617C called by : DAMPAA, CLNUT
3618C ----------------------------------------------------------------------
3619 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3620 COMMON / idfc / idff
3621 REAL PIA(4),PN(4)
3622 COMPLEX HJ(4),HJC(4)
3623C DET2(I,J)=AIMAG(HJ(I)*HJC(J)-HJ(J)*HJC(I))
3624C -- here was an error (ZW, 21.11.1991)
3625 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3626C -- it was affecting sign of A_LR asymmetry in a1 decay.
3627C -- note also collision of notation of gamma_va as defined in
3628C -- TAUOLA paper and J.H. Kuhn and Santamaria Z. Phys C 48 (1990) 445
3629* -----------------------------------
3630 IF (ktom.EQ.1.OR.ktom.EQ.-1) THEN
3631 sign= idff/abs(idff)
3632 ELSEIF (ktom.EQ.2) THEN
3633 sign=-idff/abs(idff)
3634 ELSE
3635 print *, 'STOP IN CLAXI: KTOM=',ktom
3636 stop
3637 ENDIF
3638C
3639 DO 10 i=1,4
3640 10 hjc(i)=conjg(hj(i))
3641 pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3642 pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3643 pia(3)= 2.*pn(4)*det2(1,2)
3644 pia(4)= 2.*pn(3)*det2(1,2)
3645C ALL FOUR INDICES ARE UP SO PIA(3) AND PIA(4) HAVE SAME SIGN
3646 DO 20 i=1,4
3647 20 pia(i)=pia(i)*sign
3648 END
3649 SUBROUTINE clnut(HJ,B,HV)
3650C ----------------------------------------------------------------------
3651* CALCULATES THE CONTRIBUTION BY NEUTRINO MASS
3652* NOTE THE TAU IS ASSUMED TO BE AT REST
3653C
3654C called by : DAMPAA
3655C ----------------------------------------------------------------------
3656 COMPLEX HJ(4)
3657 REAL HV(4),P(4)
3658 DATA p /3*0.,1.0/
3659C
3660 CALL claxi(hj,p,hv)
3661 b=real( hj(4)*aimag(hj(4)) - hj(3)*aimag(hj(3))
3662 & - hj(2)*aimag(hj(2)) - hj(1)*aimag(hj(1)) )
3663 RETURN
3664 END
3665 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3666C ----------------------------------------------------------------------
3667* CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
3668* FOR TAU DECAY INTO A1, A1 DECAYS NEXT INTO RHO+PI AND RHO INTO PI+PI.
3669* ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
3670* CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
3671* THE ROUTINE IS WRITEN FOR ZERO NEUTRINO MASS.
3672C
3673C called by : DPHSAA
3674C ----------------------------------------------------------------------
3675 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3676 * ,ampiz,ampi,amro,gamro,ama1,gama1
3677 * ,amk,amkz,amkst,gamkst
3678C
3679 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3680 * ,ampiz,ampi,amro,gamro,ama1,gama1
3681 * ,amk,amkz,amkst,gamkst
3682 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3683 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3684 COMMON /testa1/ keya1
3685 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3686 REAL PAA(4),VEC1(4),VEC2(4)
3687 REAL PIVEC(4),PIAKS(4),HVM(4)
3688 COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3689 DATA icont /1/
3690* THIS INLINE FUNCT. CALCULATES THE SCALAR PART OF THE PROPAGATOR
3691C AJWMOD to satisfy compiler, comment out this unused function.
3692C
3693* FOUR MOMENTUM OF A1
3694 DO 10 i=1,4
3695 vec1(i)=0.0
3696 vec2(i)=0.0
3697 hv(i) =0.0
3698 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3699 vec1(1)=1.0
3700* MASSES OF A1, AND OF TWO PI-PAIRS WHICH MAY FORM RHO
3701 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3702 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3703 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3704 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3705* ELEMENTS OF HADRON CURRENT
3706 prod1 =vec1(1)*pipl(1)
3707 prod2 =vec2(2)*pipl(2)
3708 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3709 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3710 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3711 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3712 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3713 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3714 DO 40 i=1,3
3715 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3716 40 CONTINUE
3717 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3718 DO 41 i=1,3
3719 vec1(i)= vec1(i)/gnorm
3720 41 CONTINUE
3721 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3722 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3723 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3724 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3725 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3726 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3727 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3728 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3729 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3730 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3731 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3732* HADRON CURRENT
3733 fnorm=formom(xmaa,xmom)
3734 brak=0.0
3735 DO 120 jj=1,2
3736 DO 45 i=1,4
3737 IF (jj.EQ.1) THEN
3738 hadcur(i) = fnorm *(
3739 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3740 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3741 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3742 ELSE
3743 hadcur(i) = fnorm *(
3744 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3745 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3746 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3747 ENDIF
3748 45 CONTINUE
3749C
3750* CALCULATE PI-VECTORS: VECTOR AND AXIAL
3751 CALL clvec(hadcur,pn,pivec)
3752 CALL claxi(hadcur,pn,piaks)
3753 CALL clnut(hadcur,brakm,hvm)
3754* SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
3755 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3756 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3757 DO 90 i=1,3
3758 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3759 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3760 90 CONTINUE
3761C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
3762 120 CONTINUE
3763 amplit=(gfermi*ccabib)**2*brak/2.
3764C THE STATISTICAL FACTOR FOR IDENTICAL PI-S WAS CANCELLED WITH
3765C TWO, FOR TWO MODES OF A1 DECAY NAMELLY PI+PI-PI- AND PI-PI0PI0
3766C POLARIMETER VECTOR IN TAU REST FRAME
3767 DO 91 i=1,3
3768 hv(i)=-hv(i)/brak
3769 91 CONTINUE
3770
3771 END
3772 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3773C ----------------------------------------------------------------------
3774* CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
3775* FOR TAU DECAY INTO K K pi, K pi pi.
3776* ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
3777* CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
3778C MNUM DECAY MODE IDENTIFIER.
3779C
3780C called by : DPHSAA
3781C ----------------------------------------------------------------------
3782 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3783 * ,ampiz,ampi,amro,gamro,ama1,gama1
3784 * ,amk,amkz,amkst,gamkst
3785C
3786 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3787 * ,ampiz,ampi,amro,gamro,ama1,gama1
3788 * ,amk,amkz,amkst,gamkst
3789 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3790 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3791 COMMON /ipcht/ iver
3792 INTEGER IVER
3793 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3794 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3795 REAL PIVEC(4),PIAKS(4),HVM(4)
3796 REAL FNORM(0:7),COEF
3797 DOUBLE PRECISION GETFPIRPT
3798 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3799 COMPLEX F1,F2,F3,F4,F5
3800 EXTERNAL form1,form2,form3,form4,form5
3801 DATA pi /3.141592653589793238462643/
3802 DATA icont /0/
3803CC
3804 DATA fpic /93.3e-3/
3805 IF (icont.EQ.0) THEN
3806 icont=1
3807 uroj=cmplx(0.0,1.0)
3808 dwapi0=sqrt(2.0)
3809
3810 ENDIF
3811 IF (iver.EQ.0.OR.mnum.NE.0) THEN ! so far rchl only for 3pi modes
3812 fpi=fpic
3813 ELSEIF (iver.EQ.1) THEN
3814 fpi=getfpirpt(1) ! GET defined in in ffwid3pi.f of RChL-currents
3815 ELSE
3816 write(*,*) 'wrong IVER=',iver
3817 stop
3818 ENDIF
3819 fnorm(0)=ccabib/fpi
3820 fnorm(1)=ccabib/fpi
3821 fnorm(2)=ccabib/fpi
3822 fnorm(3)=ccabib/fpi
3823 fnorm(4)=scabib/fpi/dwapi0
3824 fnorm(5)=scabib/fpi
3825 fnorm(6)=scabib/fpi
3826 fnorm(7)=ccabib/fpi
3827
3828
3829C
3830 DO 10 i=1,4
3831 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3832 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3833 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3834 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3835 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3836 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3837 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3838 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3839* ELEMENTS OF HADRON CURRENT
3840 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3841 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3842 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3843 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3844 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3845 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3846 DO 40 i=1,4
3847 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3848 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3849 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3850 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3851 CALL prod5(pim1,pim2,pim3,vec5)
3852* HADRON CURRENT
3853C be aware that sign of vec2 is opposite to sign of vec1 in a1 case
3854C Rationalize this code:
3855 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3856 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3857 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3858 f4 = (-1.0*uroj)*
3859 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3860 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3861 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3862
3863 DO 45 i=1,4
3864 hadcur(i)= cmplx(fnorm(mnum)) * (
3865 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3866 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3867 45 CONTINUE
3868C
3869* CALCULATE PI-VECTORS: VECTOR AND AXIAL
3870 CALL clvec(hadcur,pn,pivec)
3871 CALL claxi(hadcur,pn,piaks)
3872 CALL clnut(hadcur,brakm,hvm)
3873* SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
3874 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3875 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3876 amplit=(gfermi)**2*brak/2.
3877 IF (mnum.GE.9) THEN
3878 print *, 'MNUM=',mnum
3879 znak=-1.0
3880 xm1=0.0
3881 xm2=0.0
3882 xm3=0.0
3883 DO 77 k=1,4
3884 IF (k.EQ.4) znak=1.0
3885 xm1=znak*pim1(k)**2+xm1
3886 xm2=znak*pim2(k)**2+xm2
3887 xm3=znak*pim3(k)**2+xm3
3888 77 print *, 'PIM1=',pim1(k),'PIM2=',pim2(k),'PIM3=',pim3(k)
3889 print *, 'XM1=',sqrt(xm1),'XM2=',sqrt(xm2),'XM3=',sqrt(xm3)
3890 print *, '************************************************'
3891 ENDIF
3892C POLARIMETER VECTOR IN TAU REST FRAME
3893 DO 90 i=1,3
3894 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3895 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3896C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
3897 hv(i)=-hv(i)/brak
3898 90 CONTINUE
3899 END
3900 SUBROUTINE prod5(P1,P2,P3,PIA)
3901C ----------------------------------------------------------------------
3902C external product of P1, P2, P3 4-momenta.
3903C SIGN is chosen +/- for decay of TAU +/- respectively
3904C called by : DAMPAA, CLNUT
3905C ----------------------------------------------------------------------
3906 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3907 COMMON / idfc / idff
3908 REAL PIA(4),P1(4),P2(4),P3(4)
3909 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3910* -----------------------------------
3911 IF (ktom.EQ.1.OR.ktom.EQ.-1) THEN
3912 sign= idff/abs(idff)
3913 ELSEIF (ktom.EQ.2) THEN
3914 sign=-idff/abs(idff)
3915 ELSE
3916 print *, 'STOP IN PROD5: KTOM=',ktom
3917 stop
3918 ENDIF
3919C
3920C EPSILON( p1(1), p2(2), p3(3), (4) ) = 1
3921C
3922 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3923 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3924 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3925 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3926C ALL FOUR INDICES ARE UP SO PIA(3) AND PIA(4) HAVE SAME SIGN
3927 DO 20 i=1,4
3928 20 pia(i)=pia(i)*sign
3929 END
3930
3931 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3932C ----------------------------------------------------------------------
3933* THIS SIMULATES TAU DECAY IN TAU REST FRAME
3934* INTO NU A1, NEXT A1 DECAYS INTO RHO PI AND FINALLY RHO INTO PI PI.
3935* OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
3936* PAA A1
3937* PIM1 PION MINUS (OR PI0) 1 (FOR TAU MINUS)
3938* PIM2 PION MINUS (OR PI0) 2
3939* PIPL PION PLUS (OR PI-)
3940* (PIPL,PIM1) FORM A RHO
3941C ----------------------------------------------------------------------
3942 COMMON / inout / inut,iout
3943 REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3944 DATA iwarm/0/
3945C
3946 IF(mode.EQ.-1) THEN
3947C ===================
3948 iwarm=1
3949 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3950CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
3951C
3952 ELSEIF(mode.EQ. 0) THEN
3953* =======================
3954 300 CONTINUE
3955 IF(iwarm.EQ.0) GOTO 902
3956 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3957 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3958CC CALL HFILL(816,WT)
3959 CALL ranmar(rn,1)
3960 IF(rn(1).GT.wt) GOTO 300
3961C
3962 ELSEIF(mode.EQ. 1) THEN
3963* =======================
3964 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3965CC CALL HPRINT(816)
3966 ENDIF
3967C =====
3968 RETURN
3969 902 WRITE(iout, 9020)
3970 9020 FORMAT(' ----- DEXNEW: LACK OF INITIALISATION')
3971 stop
3972 END
3973 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3974C ----------------------------------------------------------------------
3975 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3976 * ,ampiz,ampi,amro,gamro,ama1,gama1
3977 * ,amk,amkz,amkst,gamkst
3978C
3979 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3980 * ,ampiz,ampi,amro,gamro,ama1,gama1
3981 * ,amk,amkz,amkst,gamkst
3982 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3983 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3984 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3985 real*4 gampmc ,gamper
3986 COMMON / inout / inut,iout
3987 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3988 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3989 & ,names
3990 CHARACTER NAMES(NMODE)*31
3991
3992 real*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3993 real*4 pdum1(4),pdum2(4),pdumi(4,9)
3994 real*4 rrr(3)
3995 real*4 wtmax(nmode)
3996 real*8 swt(nmode),sswt(nmode)
3997 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3998C
3999 DATA pi /3.141592653589793238462643/
4000 DATA iwarm/0/
4001C
4002 IF(mode.EQ.-1) THEN
4003C ===================
4004C -- AT THE MOMENT ONLY TWO DECAY MODES OF MULTIPIONS HAVE M. ELEM
4005 nmod=nmode
4006 iwarm=1
4007C PRINT 7003
4008 DO 1 jnpi=1,nmod
4009 nevraw(jnpi)=0
4010 nevacc(jnpi)=0
4011 nevovr(jnpi)=0
4012 swt(jnpi)=0
4013 sswt(jnpi)=0
4014 wtmax(jnpi)=-1.
4015C for 4pi phase space, need lots more trials at initialization,
4016C or use the WTMAX determined with many trials for default model:
4017 ntrials = 5000
4018 IF (jnpi.LE.nm4) THEN
4019C 11.Oct.11: fix for BINP and KARLSRUHE currents added
4020 wtmax(jnpi) = pkorb(3,37+jnpi)
4021 ntrials = 20000
4022 END IF
4023 DO i=1,ntrials
4024 IF (jnpi.LE.0) THEN
4025 GOTO 903
4026 ELSEIF(jnpi.LE.nm4) THEN
4027 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4028 ELSEIF(jnpi.LE.nm4+nm5) THEN
4029 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4030 ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
4031 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4032 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
4033 inum=jnpi-nm4-nm5-nm6
4034 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
4035 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
4036 inum=jnpi-nm4-nm5-nm6-nm3
4037 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
4038 ELSE
4039 GOTO 903
4040 ENDIF
4041 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
4042 ENDDO
4043C PRINT *,' DADNEW JNPI,NTRIALS,WTMAX =',JNPI,NTRIALS,WTMAX(JNPI)
4044C CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0)
4045C PRINT 7004,WTMAX(JNPI)
40461 CONTINUE
4047 WRITE(iout,7005)
4048C
4049 ELSEIF(mode.EQ. 0) THEN
4050C =======================
4051 IF(iwarm.EQ.0) GOTO 902
4052C
4053300 CONTINUE
4054 IF (jnpi.LE.0) THEN
4055 GOTO 903
4056 ELSEIF(jnpi.LE.nm4) THEN
4057 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4058 ELSEIF(jnpi.LE.nm4+nm5) THEN
4059 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4060 ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
4061 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
4062 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
4063 inum=jnpi-nm4-nm5-nm6
4064 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
4065 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
4066 inum=jnpi-nm4-nm5-nm6-nm3
4067 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
4068 ELSE
4069 GOTO 903
4070 ENDIF
4071 DO i=1,4
4072 hv(i)=-isgn*hhv(i)
4073 ENDDO
4074C CALL HFILL(801,WT/WTMAX(JNPI))
4075 nevraw(jnpi)=nevraw(jnpi)+1
4076 swt(jnpi)=swt(jnpi)+wt
4077cccM.S.>>>>>>
4078cc SSWT(JNPI)=SSWT(JNPI)+WT**2
4079 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
4080cccM.S.<<<<<<
4081 CALL ranmar(rrr,3)
4082 rn=rrr(1)
4083 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
4084 IF(rn*wtmax(jnpi).GT.wt) GOTO 300
4085C ROTATIONS TO BASIC TAU REST FRAME
4086 costhe=-1.+2.*rrr(2)
4087 thet=acos(costhe)
4088 phi =2*pi*rrr(3)
4089 CALL rotor2(thet,pnu,pnu)
4090 CALL rotor3( phi,pnu,pnu)
4091 CALL rotor2(thet,pwb,pwb)
4092 CALL rotor3( phi,pwb,pwb)
4093 CALL rotor2(thet,hv,hv)
4094 CALL rotor3( phi,hv,hv)
4095 nd=mulpik(jnpi)
4096 DO 301 i=1,nd
4097 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
4098 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
4099301 CONTINUE
4100 nevacc(jnpi)=nevacc(jnpi)+1
4101C
4102 ELSEIF(mode.EQ. 1) THEN
4103C =======================
4104 DO 500 jnpi=1,nmod
4105 IF(nevraw(jnpi).EQ.0) GOTO 500
4106 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
4107 error=0
4108 IF(nevraw(jnpi).NE.0)
4109 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
4110 rat=pargam/gamel
4111 WRITE(iout, 7010) names(jnpi),
4112 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
4113CC CALL HPRINT(801)
4114 gampmc(8+jnpi-1)=rat
4115 gamper(8+jnpi-1)=error
4116CAM NEVDEC(8+JNPI-1)=NEVACC(JNPI)
4117 500 CONTINUE
4118 ENDIF
4119C =====
4120 RETURN
4121 7003 FORMAT(///1x,15(5h*****)
4122 $ /,' *', 25x,'******** DADNEW INITIALISATION ********',9x,1h*
4123 $ )
4124 7004 FORMAT(' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
4125 7005 FORMAT(
4126 $ /,1x,15(5h*****)/)
4127 7010 FORMAT(///1x,15(5h*****)
4128 $ /,' *', 25x,'******** DADNEW FINAL REPORT ******** ',9x,1h*
4129 $ /,' *', 25x,'CHANNEL:',a31 ,9x,1h*
4130 $ /,' *',i20 ,5x,'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
4131 $ /,' *',i20 ,5x,'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
4132 $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
4133 $ /,' *',e20.5,5x,'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
4134 $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
4135 $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
4136 $ /,1x,15(5h*****)/)
4137 902 WRITE(iout, 9020)
4138 9020 FORMAT(' ----- DADNEW: LACK OF INITIALISATION')
4139 stop
4140 903 WRITE(iout, 9030) jnpi,mode
4141 9030 FORMAT(' ----- DADNEW: WRONG JNPI',2i5)
4142 stop
4143 END
4144
4145
4146 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4147C ----------------------------------------------------------------------
4148* IT SIMULATES A1 DECAY IN TAU REST FRAME WITH
4149* Z-AXIS ALONG A1 MOMENTUM
4150C ----------------------------------------------------------------------
4151 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4152 * ,ampiz,ampi,amro,gamro,ama1,gama1
4153 * ,amk,amkz,amkst,gamkst
4154C
4155 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4156 * ,ampiz,ampi,amro,gamro,ama1,gama1
4157 * ,amk,amkz,amkst,gamkst
4158 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4159 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4160 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
4161 REAL PR(4),PIZ(4)
4162 real*4 rrr(9)
4163 real*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
4164 DATA pi /3.141592653589793238462643/
4165 DATA icont /0/
4166 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
4167C AMRO, GAMRO IS ONLY A PARAMETER FOR GETING HIGHT EFFICIENCY
4168C
4169C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
4170C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P)
4171 phspac=1./2**23/pi**11
4172 phsp=1./2**5/pi**2
4173 IF (jnpi.EQ.1) THEN
4174 prez=0.7
4175 amp1=ampi
4176 amp2=ampi
4177 amp3=ampi
4178 amp4=ampiz
4179 amrx=pkorb(1,14)
4180 gamrx=pkorb(2,14)
4181C AJW: cant simply change AMROP, etc, here!
4182C CHOICE is a by-hand tuning/optimization, no simple relationship
4183C to actual resonance masses (accd to Z.Was).
4184C What matters in the end is what you put in formf/curr .
4185 amrop =1.2
4186 gamrop=.46
4187 ELSE
4188 prez=0.0
4189 amp1=ampiz
4190 amp2=ampiz
4191 amp3=ampiz
4192 amp4=ampi
4193 amrx=1.4
4194 gamrx=.6
4195 amrop =amrx
4196 gamrop=gamrx
4197
4198 ENDIF
4199 rrb=0.3
4200 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
4201 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4202 prez=prob1+prob2
4203C TAU MOMENTUM
4204 pt(1)=0.
4205 pt(2)=0.
4206 pt(3)=0.
4207 pt(4)=amtau
4208C
4209 CALL ranmar(rrr,9)
4210C
4211* MASSES OF 4, 3 AND 2 PI SYSTEMS
4212C 3 PI WITH SAMPLING FOR RESONANCE
4213CAM
4214 rr1=rrr(6)
4215 ams1=(amp1+amp2+amp3+amp4)**2
4216 ams2=(amtau-amnuta)**2
4217 alp1=atan((ams1-amrop**2)/amrop/gamrop)
4218 alp2=atan((ams2-amrop**2)/amrop/gamrop)
4219 alp=alp1+rr1*(alp2-alp1)
4220 am4sq =amrop**2+amrop*gamrop*tan(alp)
4221 am4 =sqrt(am4sq)
4222 phspac=phspac*
4223 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4224 phspac=phspac*(alp2-alp1)
4225
4226C
4227 rr1=rrr(1)
4228 ams1=(amp2+amp3+amp4)**2
4229 ams2=(am4-amp1)**2
4230 IF (rrr(9).GT.prez) THEN
4231 am3sq=ams1+ rr1*(ams2-ams1)
4232 am3 =sqrt(am3sq)
4233C --- this part of jacobian will be recovered later
4234 ff1=ams2-ams1
4235 ELSE
4236* PHASE SPACE WITH SAMPLING FOR OMEGA RESONANCE,
4237 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4238 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4239 alp=alp1+rr1*(alp2-alp1)
4240 am3sq =amrx**2+amrx*gamrx*tan(alp)
4241 am3 =sqrt(am3sq)
4242C --- THIS PART OF THE JACOBIAN WILL BE RECOVERED LATER ---------------
4243 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4244 ff1=ff1*(alp2-alp1)
4245 ENDIF
4246C MASS OF 2
4247 rr2=rrr(2)
4248 ams1=(amp3+amp4)**2
4249 ams2=(am3-amp2)**2
4250* FLAT PHASE SPACE;
4251 am2sq=ams1+ rr2*(ams2-ams1)
4252 am2 =sqrt(am2sq)
4253C --- this part of jacobian will be recovered later
4254 ff2=(ams2-ams1)
4255* 2 RESTFRAME, DEFINE PIZ AND PIPL
4256 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4257 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4258 ppi= enq1**2-amp4**2
4259 pppi=sqrt(abs(enq1**2-amp4**2))
4260 phspac=phspac*(4*pi)*(2*pppi/am2)
4261* PIZ MOMENTUM IN 2 REST FRAME
4262 CALL sphera(pppi,piz)
4263 piz(4)=enq1
4264* PIPL MOMENTUM IN 2 REST FRAME
4265 DO 30 i=1,3
4266 30 pipl(i)=-piz(i)
4267 pipl(4)=enq2
4268* 3 REST FRAME, DEFINE PIM1
4269* PR MOMENTUM
4270 pr(1)=0
4271 pr(2)=0
4272 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4273 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4274 ppi = pr(4)**2-am2**2
4275* PIM1 MOMENTUM
4276 pim1(1)=0
4277 pim1(2)=0
4278 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4279 pim1(3)=-pr(3)
4280C --- this part of jacobian will be recovered later
4281 ff3=(4*pi)*(2*pr(3)/am3)
4282* OLD PIONS BOOSTED FROM 2 REST FRAME TO 3 REST FRAME
4283 exe=(pr(4)+pr(3))/am2
4284 CALL bostr3(exe,piz,piz)
4285 CALL bostr3(exe,pipl,pipl)
4286 rr3=rrr(3)
4287 rr4=rrr(4)
4288 thet =acos(-1.+2*rr3)
4289 phi = 2*pi*rr4
4290 CALL rotpol(thet,phi,pipl)
4291 CALL rotpol(thet,phi,pim1)
4292 CALL rotpol(thet,phi,piz)
4293 CALL rotpol(thet,phi,pr)
4294* 4 REST FRAME, DEFINE PIM2
4295* PR MOMENTUM
4296 pr(1)=0
4297 pr(2)=0
4298 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4299 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4300 ppi = pr(4)**2-am3**2
4301* PIM2 MOMENTUM
4302 pim2(1)=0
4303 pim2(2)=0
4304 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4305 pim2(3)=-pr(3)
4306C --- this part of jacobian will be recovered later
4307 ff4=(4*pi)*(2*pr(3)/am4)
4308* OLD PIONS BOOSTED FROM 3 REST FRAME TO 4 REST FRAME
4309 exe=(pr(4)+pr(3))/am3
4310 CALL bostr3(exe,piz,piz)
4311 CALL bostr3(exe,pipl,pipl)
4312 CALL bostr3(exe,pim1,pim1)
4313 rr3=rrr(7)
4314 rr4=rrr(8)
4315 thet =acos(-1.+2*rr3)
4316 phi = 2*pi*rr4
4317 CALL rotpol(thet,phi,pipl)
4318 CALL rotpol(thet,phi,pim1)
4319 CALL rotpol(thet,phi,pim2)
4320 CALL rotpol(thet,phi,piz)
4321 CALL rotpol(thet,phi,pr)
4322C
4323* NOW TO THE TAU REST FRAME, DEFINE PAA AND NEUTRINO MOMENTA
4324* PAA MOMENTUM
4325 paa(1)=0
4326 paa(2)=0
4327 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4328 paa(3)= sqrt(abs(paa(4)**2-am4**2))
4329 ppi = paa(4)**2-am4**2
4330 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4331 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4332* TAU-NEUTRINO MOMENTUM
4333 pn(1)=0
4334 pn(2)=0
4335 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4336 pn(3)=-paa(3)
4337C ZBW 20.12.2002 bug fix
4338 IF(rrr(9).LE.0.5*prez) THEN
4339 DO 72 i=1,4
4340 x=pim1(i)
4341 pim1(i)=pim2(i)
4342 72 pim2(i)=x
4343 ENDIF
4344C end of bug fix
4345C WE INCLUDE REMAINING PART OF THE JACOBIAN
4346C --- FLAT CHANNEL
4347 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4348 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4349 ams2=(am4-amp2)**2
4350 ams1=(amp1+amp3+amp4)**2
4351 ff1=(ams2-ams1)
4352 ams1=(amp3+amp4)**2
4353 ams2=(sqrt(am3sq)-amp1)**2
4354 ff2=ams2-ams1
4355 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4356 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4357 uu=ff1*ff2*ff3*ff4
4358C --- FIRST CHANNEL
4359 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4360 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4361 ams2=(am4-amp2)**2
4362 ams1=(amp1+amp3+amp4)**2
4363 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4364 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4365 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4366 ff1=ff1*(alp2-alp1)
4367 ams1=(amp3+amp4)**2
4368 ams2=(sqrt(am3sq)-amp1)**2
4369 ff2=ams2-ams1
4370 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4371 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4372 ff=ff1*ff2*ff3*ff4
4373C --- SECOND CHANNEL
4374 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4375 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4376 ams2=(am4-amp1)**2
4377 ams1=(amp2+amp3+amp4)**2
4378 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4379 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4380 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4381 gg1=gg1*(alp2-alp1)
4382 ams1=(amp3+amp4)**2
4383 ams2=(sqrt(am3sq)-amp2)**2
4384 gg2=ams2-ams1
4385 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4386 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4387 gg=gg1*gg2*gg3*gg4
4388C --- JACOBIAN AVERAGED OVER THE TWO
4389 ! 05.10.2011 missing factor in IF( (0.5*PREZ* ... (1-PREZ)* ...) added
4390 IF ( (0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg).GT.0.0d0) THEN
4391 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4392 phspac=phspac*rr
4393 ELSE
4394 phspac=0.0
4395 ENDIF
4396* MOMENTA OF THE TWO PI-MINUS ARE RANDOMLY SYMMETRISED
4397 IF (jnpi.EQ.1) THEN
4398 rr5= rrr(5)
4399 IF(rr5.LE.0.5) THEN
4400 DO 70 i=1,4
4401 x=pim1(i)
4402 pim1(i)=pim2(i)
4403 70 pim2(i)=x
4404 ENDIF
4405 phspac=phspac/2.
4406 ELSE
4407C MOMENTA OF PI0-S ARE GENERATED UNIFORMLY ONLY IF PREZ=0.0
4408 rr5= rrr(5)
4409 IF(rr5.LE.0.5) THEN
4410 DO 71 i=1,4
4411 x=pim1(i)
4412 pim1(i)=pim2(i)
4413 71 pim2(i)=x
4414 ENDIF
4415 phspac=phspac/6.
4416 ENDIF
4417* ALL PIONS BOOSTED FROM 4 REST FRAME TO TAU REST FRAME
4418* Z-AXIS ANTIPARALLEL TO NEUTRINO MOMENTUM
4419 exe=(paa(4)+paa(3))/am4
4420 CALL bostr3(exe,piz,piz)
4421 CALL bostr3(exe,pipl,pipl)
4422 CALL bostr3(exe,pim1,pim1)
4423 CALL bostr3(exe,pim2,pim2)
4424 CALL bostr3(exe,pr,pr)
4425C PARTIAL WIDTH CONSISTS OF PHASE SPACE AND AMPLITUDE
4426C CHECK ON CONSISTENCY WITH DADNPI, THEN, CODE BREAKES UNIFORM PION
4427C DISTRIBUTION IN HADRONIC SYSTEM
4428CAM Assume neutrino mass=0. and sum over final polarisation
4429C AMX2=AM4**2
4430C BRAK= 2*(AMTAU**2-AMX2) * (AMTAU**2+2.*AMX2)
4431C AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AMX2*SIGEE(AMX2,1)
4432 IF (jnpi.EQ.1) THEN
4433 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4434 ELSEIF (jnpi.EQ.2) THEN
4435 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4436 ENDIF
4437 dgamt=1/(2.*amtau)*amplit*phspac
4438C PHASE SPACE CHECK
4439C DGAMT=PHSPAC
4440 DO 77 k=1,4
4441 pmult(k,1)=pim1(k)
4442 pmult(k,2)=pim2(k)
4443 pmult(k,3)=pipl(k)
4444 pmult(k,4)=piz(k)
4445 77 CONTINUE
4446 END
4447 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4448C ----------------------------------------------------------------------
4449* CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
4450* FOR TAU DECAY INTO 4 PI MODES
4451* ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
4452* CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
4453C MNUM DECAY MODE IDENTIFIER.
4454C
4455C called by : DPHSAA
4456C ----------------------------------------------------------------------
4457 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4458 * ,ampiz,ampi,amro,gamro,ama1,gama1
4459 * ,amk,amkz,amkst,gamkst
4460C
4461 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4462 * ,ampiz,ampi,amro,gamro,ama1,gama1
4463 * ,amk,amkz,amkst,gamkst
4464 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4465 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4466 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4467 REAL PIVEC(4),PIAKS(4),HVM(4)
4468 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4469 EXTERNAL form1,form2,form3,form4,form5
4470 DATA pi /3.141592653589793238462643/
4471 DATA icont /0/
4472C
4473 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4474C
4475* CALCULATE PI-VECTORS: VECTOR AND AXIAL
4476 CALL clvec(hadcur,pn,pivec)
4477 CALL claxi(hadcur,pn,piaks)
4478 CALL clnut(hadcur,brakm,hvm)
4479* SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
4480 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4481 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4482 amplit=(ccabib*gfermi)**2*brak/2.
4483C POLARIMETER VECTOR IN TAU REST FRAME
4484 DO 90 i=1,3
4485 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4486 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4487C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
4488 hv(i)=-hv(i)/brak
4489 90 CONTINUE
4490 END
4491 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4492C ----------------------------------------------------------------------
4493* IT SIMULATES 5pi DECAY IN TAU REST FRAME WITH
4494* Z-AXIS ALONG 5pi MOMENTUM
4495C ----------------------------------------------------------------------
4496 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4497 * ,ampiz,ampi,amro,gamro,ama1,gama1
4498 * ,amk,amkz,amkst,gamkst
4499C
4500 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4501 * ,ampiz,ampi,amro,gamro,ama1,gama1
4502
4503
4504 * ,amk,amkz,amkst,gamkst
4505 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4506 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4507 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4508 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4509 & ,names
4510 CHARACTER NAMES(NMODE)*31
4511 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4512 real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4513 real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4514 real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4515 real*4 rrr(10)
4516 real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4517 real*8 xm,am,gamma
4518ccM.S.>>>>>>
4519 real*8 phspac
4520ccM.S.<<<<<<
4521 DATA pi /3.141592653589793238462643/
4522 DATA icont /0/
4523 data fpi /93.3e-3/
4524c
4525 COMPLEX BWIGN
4526C
4527 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4528
4529C
4530 amom=.782
4531 gamom=0.0085
4532c
4533C 6 BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
4534C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P)
4535 phspac=1./2**29/pi**14
4536c PHSPAC=1./2**5/PI**2
4537C init 5pi decay mode (JNPI)
4538 amp1=dcdmas(idffin(1,jnpi))
4539 amp2=dcdmas(idffin(2,jnpi))
4540 amp3=dcdmas(idffin(3,jnpi))
4541 amp4=dcdmas(idffin(4,jnpi))
4542 amp5=dcdmas(idffin(5,jnpi))
4543c
4544C TAU MOMENTUM
4545 pt(1)=0.
4546 pt(2)=0.
4547 pt(3)=0.
4548 pt(4)=amtau
4549C
4550 CALL ranmar(rrr,10)
4551C
4552c masses of 5, 4, 3 and 2 pi systems
4553c 3 pi with sampling for omega resonance
4554cam
4555c mass of 5 (12345)
4556 rr1=rrr(10)
4557 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4558 ams2=(amtau-amnuta)**2
4559 am5sq=ams1+ rr1*(ams2-ams1)
4560 am5 =sqrt(am5sq)
4561 phspac=phspac*(ams2-ams1)
4562c
4563c mass of 4 (2345)
4564c flat phase space
4565 rr1=rrr(9)
4566 ams1=(amp2+amp3+amp4+amp5)**2
4567 ams2=(am5-amp1)**2
4568 am4sq=ams1+ rr1*(ams2-ams1)
4569 am4 =sqrt(am4sq)
4570 gg1=ams2-ams1
4571c
4572c mass of 3 (234)
4573C phase space with sampling for omega resonance
4574 rr1=rrr(1)
4575 ams1=(amp2+amp3+amp4)**2
4576 ams2=(am4-amp5)**2
4577 alp1=atan((ams1-amom**2)/amom/gamom)
4578 alp2=atan((ams2-amom**2)/amom/gamom)
4579 alp=alp1+rr1*(alp2-alp1)
4580 am3sq =amom**2+amom*gamom*tan(alp)
4581 am3 =sqrt(am3sq)
4582c --- this part of the jacobian will be recovered later ---------------
4583 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4584 gg2=gg2*(alp2-alp1)
4585c flat phase space;
4586C am3sq=ams1+ rr1*(ams2-ams1)
4587C am3 =sqrt(am3sq)
4588c --- this part of jacobian will be recovered later
4589C gg2=ams2-ams1
4590c
4591C mass of 2 (34)
4592 rr2=rrr(2)
4593 ams1=(amp3+amp4)**2
4594 ams2=(am3-amp2)**2
4595c flat phase space;
4596 am2sq=ams1+ rr2*(ams2-ams1)
4597 am2 =sqrt(am2sq)
4598c --- this part of jacobian will be recovered later
4599 gg3=ams2-ams1
4600c
4601c (34) restframe, define pi3 and pi4
4602 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4603 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4604 ppi= enq1**2-amp3**2
4605 pppi=sqrt(abs(enq1**2-amp3**2))
4606 ff1=(4*pi)*(2*pppi/am2)
4607c pi3 momentum in (34) rest frame
4608 call sphera(pppi,pi3)
4609 pi3(4)=enq1
4610c pi4 momentum in (34) rest frame
4611 do 30 i=1,3
4612 30 pi4(i)=-pi3(i)
4613 pi4(4)=enq2
4614c
4615c (234) rest frame, define pi2
4616c pr momentum
4617 pr(1)=0
4618 pr(2)=0
4619 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4620 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4621 ppi = pr(4)**2-am2**2
4622c pi2 momentum
4623 pi2(1)=0
4624 pi2(2)=0
4625 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4626 pi2(3)=-pr(3)
4627c --- this part of jacobian will be recovered later
4628 ff2=(4*pi)*(2*pr(3)/am3)
4629c old pions boosted from 2 rest frame to 3 rest frame
4630 exe=(pr(4)+pr(3))/am2
4631 call bostr3(exe,pi3,pi3)
4632 call bostr3(exe,pi4,pi4)
4633 rr3=rrr(3)
4634 rr4=rrr(4)
4635 thet =acos(-1.+2*rr3)
4636 phi = 2*pi*rr4
4637 call rotpol(thet,phi,pi2)
4638 call rotpol(thet,phi,pi3)
4639 call rotpol(thet,phi,pi4)
4640C
4641C (2345) rest frame, define pi5
4642c pr momentum
4643 pr(1)=0
4644 pr(2)=0
4645 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4646 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4647 ppi = pr(4)**2-am3**2
4648c pi5 momentum
4649 pi5(1)=0
4650 pi5(2)=0
4651 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4652 pi5(3)=-pr(3)
4653c --- this part of jacobian will be recovered later
4654 ff3=(4*pi)*(2*pr(3)/am4)
4655c old pions boosted from 3 rest frame to 4 rest frame
4656 exe=(pr(4)+pr(3))/am3
4657 call bostr3(exe,pi2,pi2)
4658 call bostr3(exe,pi3,pi3)
4659 call bostr3(exe,pi4,pi4)
4660 rr3=rrr(5)
4661 rr4=rrr(6)
4662 thet =acos(-1.+2*rr3)
4663 phi = 2*pi*rr4
4664 call rotpol(thet,phi,pi2)
4665 call rotpol(thet,phi,pi3)
4666 call rotpol(thet,phi,pi4)
4667 call rotpol(thet,phi,pi5)
4668C
4669C (12345) rest frame, define pi1
4670c pr momentum
4671 pr(1)=0
4672 pr(2)=0
4673 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4674 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4675 ppi = pr(4)**2-am4**2
4676c pi1 momentum
4677 pi1(1)=0
4678 pi1(2)=0
4679 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4680 pi1(3)=-pr(3)
4681c --- this part of jacobian will be recovered later
4682 ff4=(4*pi)*(2*pr(3)/am5)
4683c old pions boosted from 4 rest frame to 5 rest frame
4684 exe=(pr(4)+pr(3))/am4
4685 call bostr3(exe,pi2,pi2)
4686 call bostr3(exe,pi3,pi3)
4687 call bostr3(exe,pi4,pi4)
4688 call bostr3(exe,pi5,pi5)
4689 rr3=rrr(7)
4690 rr4=rrr(8)
4691 thet =acos(-1.+2*rr3)
4692 phi = 2*pi*rr4
4693 call rotpol(thet,phi,pi1)
4694 call rotpol(thet,phi,pi2)
4695 call rotpol(thet,phi,pi3)
4696 call rotpol(thet,phi,pi4)
4697 call rotpol(thet,phi,pi5)
4698c
4699* now to the tau rest frame, define paa and neutrino momenta
4700* paa momentum
4701 paa(1)=0
4702 paa(2)=0
4703c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4704c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4705c ppi = paa(4)**2-am5**2
4706 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4707 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4708 ppi = paa(4)**2-am5sq
4709 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4710* tau-neutrino momentum
4711 pn(1)=0
4712 pn(2)=0
4713 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4714 pn(3)=-paa(3)
4715c
4716 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4717c
4718C all pions boosted from 5 rest frame to tau rest frame
4719C z-axis antiparallel to neutrino momentum
4720 exe=(paa(4)+paa(3))/am5
4721 call bostr3(exe,pi1,pi1)
4722 call bostr3(exe,pi2,pi2)
4723 call bostr3(exe,pi3,pi3)
4724 call bostr3(exe,pi4,pi4)
4725 call bostr3(exe,pi5,pi5)
4726c
4727C partial width consists of phase space and amplitude
4728C AMPLITUDE (cf YS.Tsai Phys.Rev.D4,2821(1971)
4729C or F.Gilman SH.Rhie Phys.Rev.D31,1066(1985)
4730C
4731 pxq=amtau*paa(4)
4732 pxn=amtau*pn(4)
4733 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4734 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4735 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4736 fompp = cabs(bwign(am3,amom,gamom))**2
4737c normalisation factor (to some numerical undimensioned factor;
4738c cf R.Fischer et al ZPhys C3, 313 (1980))
4739 fnorm = 1/fpi**6
4740c AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AM5SQ*SIGEE(AM5SQ,JNPI)
4741 amplit=ccabib**2*gfermi**2/2. * brak
4742 amplit = amplit * fompp * fnorm
4743c phase space test
4744c amplit = amplit * fnorm
4745 dgamt=1/(2.*amtau)*amplit*phspac
4746c ignore spin terms
4747 DO 40 i=1,3
4748 40 hv(i)=0.
4749c
4750 do 77 k=1,4
4751 pmult(k,1)=pi1(k)
4752 pmult(k,2)=pi2(k)
4753 pmult(k,3)=pi3(k)
4754 pmult(k,4)=pi4(k)
4755 pmult(k,5)=pi5(k)
4756 77 continue
4757 return
4758C missing: transposition of identical particles, startistical factors
4759C for identical matrices, polarimetric vector. Matrix element rather naive.
4760C flat phase space in pion system + with breit wigner for omega
4761C anyway it is better than nothing, and code is improvable.
4762 end
4763 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4764C ----------------------------------------------------------------------
4765C IT SIMULATES RHO DECAY IN TAU REST FRAME WITH
4766C Z-AXIS ALONG RHO MOMENTUM
4767C Rho decays to K Kbar
4768C ----------------------------------------------------------------------
4769 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4770 * ,ampiz,ampi,amro,gamro,ama1,gama1
4771 * ,amk,amkz,amkst,gamkst
4772C
4773 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4774 * ,ampiz,ampi,amro,gamro,ama1,gama1
4775 * ,amk,amkz,amkst,gamkst
4776
4777 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4778 REAL RR1(1)
4779 DATA pi /3.141592653589793238462643/
4780 DATA icont /0/
4781C
4782C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
4783 phspac=1./2**11/pi**5
4784C TAU MOMENTUM
4785 pt(1)=0.
4786 pt(2)=0.
4787 pt(3)=0.
4788 pt(4)=amtau
4789C MASS OF (REAL/VIRTUAL) RHO
4790 ams1=(amk+amkz)**2
4791 ams2=(amtau-amnuta)**2
4792C FLAT PHASE SPACE
4793 CALL ranmar(rr1,1)
4794 amx2=ams1+ rr1(1)*(ams2-ams1)
4795 amx=sqrt(amx2)
4796 phspac=phspac*(ams2-ams1)
4797C PHASE SPACE WITH SAMPLING FOR RHO RESONANCE
4798c ALP1=ATAN((AMS1-AMRO**2)/AMRO/GAMRO)
4799c ALP2=ATAN((AMS2-AMRO**2)/AMRO/GAMRO)
4800CAM
4801 100 CONTINUE
4802c CALL RANMAR(RR1,1)
4803c ALP=ALP1+RR1(1)*(ALP2-ALP1)
4804c AMX2=AMRO**2+AMRO*GAMRO*TAN(ALP)
4805c AMX=SQRT(AMX2)
4806c IF(AMX.LT.(AMK+AMKZ)) GO TO 100
4807CAM
4808c PHSPAC=PHSPAC*((AMX2-AMRO**2)**2+(AMRO*GAMRO)**2)/(AMRO*GAMRO)
4809c PHSPAC=PHSPAC*(ALP2-ALP1)
4810C
4811C TAU-NEUTRINO MOMENTUM
4812 pn(1)=0
4813 pn(2)=0
4814 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4815 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
4816C RHO MOMENTUM
4817 pr(1)=0
4818 pr(2)=0
4819 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4820 pr(3)=-pn(3)
4821 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4822C
4823CAM
4824 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4825 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4826 pppi=sqrt(abs(enq1-amk)*(enq1+amk))
4827 phspac=phspac*(4*pi)*(2*pppi/amx)
4828C CHARGED PI MOMENTUM IN RHO REST FRAME
4829 CALL sphera(pppi,pkc)
4830 pkc(4)=enq1
4831C NEUTRAL PI MOMENTUM IN RHO REST FRAME
4832 DO 20 i=1,3
483320 pkz(i)=-pkc(i)
4834 pkz(4)=enq2
4835 exe=(pr(4)+pr(3))/amx
4836C PIONS BOOSTED FROM RHO REST FRAME TO TAU REST FRAME
4837 CALL bostr3(exe,pkc,pkc)
4838 CALL bostr3(exe,pkz,pkz)
4839
4840
4841 CALL dam2pi(3,pt,pn,pkc,pkz,amplit,hv)
4842 dgamt=1/(2.*amtau)*amplit*phspac
4843
4844 do 77 k=1,4
4845 pmult(k,1)=pkc(k)
4846 pmult(k,2)=pkz(k)
4847 77 continue
4848 RETURN
4849 END
4850 FUNCTION fpirk(W)
4851C ----------------------------------------------------------
4852c square of pion form factor
4853C ----------------------------------------------------------
4854 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4855 * ,ampiz,ampi,amro,gamro,ama1,gama1
4856 * ,amk,amkz,amkst,gamkst
4857C
4858 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4859 * ,ampiz,ampi,amro,gamro,ama1,gama1
4860 * ,amk,amkz,amkst,gamkst
4861c COMPLEX FPIKMK
4862 COMPLEX FPIKM
4863 fpirk=cabs(fpikm(w,amk,amkz))**2
4864c FPIRK=CABS(FPIKMK(W,AMK,AMKZ))**2
4865 END
4866 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4867C **********************************************************
4868C Kaon form factor
4869C **********************************************************
4870 COMPLEX BWIGM
4871 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4872 EXTERNAL bwig
4873 DATA init /0/
4874C
4875C ------------ PARAMETERS --------------------
4876 IF (init.EQ.0 ) THEN
4877 init=1
4878 pi=3.141592654
4879 pim=.140
4880 rom=0.773
4881 rog=0.145
4882 rom1=1.570
4883 rog1=0.510
4884c BETA1=-0.111
4885 beta1=-0.221
4886 ENDIF
4887C -----------------------------------------------
4888 s=w**2
4889 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4890 & /(1+beta1)
4891 RETURN
4892 END
4893 SUBROUTINE reslux
4894C ****************
4895C INITIALIZE LUND COMMON
4896 nhep=0
4897 END
4898 SUBROUTINE dwrph(KTO,PHX)
4899C
4900C -------------------------
4901C
4902 IMPLICIT REAL*8 (a-h,o-z)
4903 real*4 phx(4)
4904 real*4 qhot(4)
4905C
4906 DO 9 k=1,4
4907 qhot(k) =0.0
4908 9 CONTINUE
4909C CASE OF TAU RADIATIVE DECAYS.
4910C FILLING OF THE LUND COMMON BLOCK.
4911 DO 1002 i=1,4
4912 1002 qhot(i)=phx(i)
4913 IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
4914 RETURN
4915 END
4916 SUBROUTINE dwluph(KTO,PHOT)
4917C---------------------------------------------------------------------
4918C Lorentz transformation to CMsystem and
4919C Updating of HEPEVT record
4920C
4921C called by : DEXAY1,(DEKAY1,DEKAY2)
4922C
4923C used when radiative corrections in decays are generated
4924C---------------------------------------------------------------------
4925C
4926 REAL PHOT(4)
4927 COMMON /taupos/ np1,np2
4928C
4929C check energy
4930 IF (phot(4).LE.0.0) RETURN
4931C
4932C position of decaying particle:
4933 IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
4934 nps=np1
4935 ELSE
4936 nps=np2
4937 ENDIF
4938C
4939 ktos=kto
4940 IF(ktos.GT.10) ktos=ktos-10
4941C boost and append photon (gamma is 22)
4942 CALL tralo4(ktos,phot,phot,am)
4943 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4944C
4945 RETURN
4946 END
4947
4948 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4949C ----------------------------------------------------------------------
4950C Lorentz transformation to CMsystem and
4951C Updating of HEPEVT record
4952C
4953C ISGN = 1/-1 for tau-/tau+
4954C
4955C called by : DEXAY,(DEKAY1,DEKAY2)
4956C ----------------------------------------------------------------------
4957C
4958 REAL PNU(4),PWB(4),PEL(4),PNE(4)
4959 COMMON /taupos/ np1,np2
4960C
4961C position of decaying particle:
4962 IF(kto.EQ. 1) THEN
4963 nps=np1
4964 ELSE
4965 nps=np2
4966 ENDIF
4967C
4968C tau neutrino (nu_tau is 16)
4969 CALL tralo4(kto,pnu,pnu,am)
4970 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4971C
4972C W boson (W+ is 24)
4973 CALL tralo4(kto,pwb,pwb,am)
4974C CALL FILHEP(0,2,-24*ISGN,NPS,NPS,0,0,PWB,AM,.TRUE.)
4975C
4976C electron (e- is 11)
4977 CALL tralo4(kto,pel,pel,am)
4978 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4979C
4980C anti electron neutrino (nu_e is 12)
4981 CALL tralo4(kto,pne,pne,am)
4982 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4983C
4984 RETURN
4985 END
4986 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4987C ----------------------------------------------------------------------
4988C Lorentz transformation to CMsystem and
4989C Updating of HEPEVT record
4990C
4991C ISGN = 1/-1 for tau-/tau+
4992C
4993C called by : DEXAY,(DEKAY1,DEKAY2)
4994C ----------------------------------------------------------------------
4995C
4996 REAL PNU(4),PWB(4),PMU(4),PNM(4)
4997 COMMON /taupos/ np1,np2
4998C
4999C position of decaying particle:
5000 IF(kto.EQ. 1) THEN
5001 nps=np1
5002 ELSE
5003 nps=np2
5004 ENDIF
5005C
5006C tau neutrino (nu_tau is 16)
5007 CALL tralo4(kto,pnu,pnu,am)
5008 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5009C
5010C W boson (W+ is 24)
5011 CALL tralo4(kto,pwb,pwb,am)
5012C CALL FILHEP(0,2,-24*ISGN,NPS,NPS,0,0,PWB,AM,.TRUE.)
5013C
5014C muon (mu- is 13)
5015 CALL tralo4(kto,pmu,pmu,am)
5016 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5017C
5018C anti muon neutrino (nu_mu is 14)
5019 CALL tralo4(kto,pnm,pnm,am)
5020 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5021C
5022 RETURN
5023 END
5024 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5025C ----------------------------------------------------------------------
5026C Lorentz transformation to CMsystem and
5027C Updating of HEPEVT record
5028C
5029C ISGN = 1/-1 for tau-/tau+
5030C
5031C called by : DEXAY,(DEKAY1,DEKAY2)
5032C ----------------------------------------------------------------------
5033C
5034 REAL PNU(4),PPI(4)
5035 COMMON /taupos/ np1,np2
5036C
5037C position of decaying particle:
5038 IF(kto.EQ. 1) THEN
5039 nps=np1
5040 ELSE
5041 nps=np2
5042 ENDIF
5043C
5044C tau neutrino (nu_tau is 16)
5045 CALL tralo4(kto,pnu,pnu,am)
5046 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5047C
5048C charged pi meson (pi+ is 211)
5049 CALL tralo4(kto,ppi,ppi,am)
5050 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5051C
5052 RETURN
5053 END
5054 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5055C ----------------------------------------------------------------------
5056C Lorentz transformation to CMsystem and
5057C Updating of HEPEVT record
5058C
5059C ISGN = 1/-1 for tau-/tau+
5060C
5061C called by : DEXAY,(DEKAY1,DEKAY2)
5062C ----------------------------------------------------------------------
5063C
5064 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5065 COMMON /taupos/ np1,np2
5066C
5067C position of decaying particle:
5068 IF(kto.EQ. 1) THEN
5069 nps=np1
5070 ELSE
5071 nps=np2
5072 ENDIF
5073C
5074C tau neutrino (nu_tau is 16)
5075 CALL tralo4(kto,pnu,pnu,am)
5076 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5077C
5078C charged rho meson (rho+ is 213)
5079 CALL tralo4(kto,prho,prho,am)
5080 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5081C
5082C charged pi meson (pi+ is 211)
5083 CALL tralo4(kto,pic,pic,am)
5084 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5085C
5086C pi0 meson (pi0 is 111)
5087 CALL tralo4(kto,piz,piz,am)
5088 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5089C
5090 RETURN
5091 END
5092 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5093C ----------------------------------------------------------------------
5094C Lorentz transformation to CMsystem and
5095C Updating of HEPEVT record
5096C
5097C ISGN = 1/-1 for tau-/tau+
5098C JAA = 1 (2) FOR A_1- DECAY TO PI+ 2PI- (PI- 2PI0)
5099C
5100C called by : DEXAY,(DEKAY1,DEKAY2)
5101C ----------------------------------------------------------------------
5102C
5103 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5104 COMMON /taupos/ np1,np2
5105C
5106C position of decaying particle:
5107 IF(kto.EQ. 1) THEN
5108 nps=np1
5109 ELSE
5110 nps=np2
5111 ENDIF
5112C
5113C tau neutrino (nu_tau is 16)
5114 CALL tralo4(kto,pnu,pnu,am)
5115 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5116C
5117C charged a_1 meson (a_1+ is 20213)
5118 CALL tralo4(kto,paa,paa,am)
5119 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5120C
5121C two possible decays of the charged a1 meson
5122 IF(jaa.EQ.1) THEN
5123C
5124C A1 --> PI+ PI- PI- (or charged conjugate)
5125C
5126C pi minus (or c.c.) (pi+ is 211)
5127 CALL tralo4(kto,pim2,pim2,am)
5128 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5129C
5130C pi minus (or c.c.) (pi+ is 211)
5131 CALL tralo4(kto,pim1,pim1,am)
5132 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5133C
5134C pi plus (or c.c.) (pi+ is 211)
5135 CALL tralo4(kto,pipl,pipl,am)
5136 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5137C
5138 ELSE IF (jaa.EQ.2) THEN
5139C
5140C A1 --> PI- PI0 PI0 (or charged conjugate)
5141C
5142C pi zero (pi0 is 111)
5143 CALL tralo4(kto,pim2,pim2,am)
5144 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5145C
5146C pi zero (pi0 is 111)
5147 CALL tralo4(kto,pim1,pim1,am)
5148 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5149C
5150C pi minus (or c.c.) (pi+ is 211)
5151 CALL tralo4(kto,pipl,pipl,am)
5152 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5153C
5154 ENDIF
5155C
5156 RETURN
5157 END
5158 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5159C ----------------------------------------------------------------------
5160C Lorentz transformation to CMsystem and
5161C Updating of HEPEVT record
5162C
5163C ISGN = 1/-1 for tau-/tau+
5164C
5165C ----------------------------------------------------------------------
5166C
5167 REAL PKK(4),PNU(4)
5168 COMMON /taupos/ np1,np2
5169C
5170C position of decaying particle
5171 IF (kto.EQ.1) THEN
5172 nps=np1
5173 ELSE
5174 nps=np2
5175 ENDIF
5176C
5177C tau neutrino (nu_tau is 16)
5178 CALL tralo4 (kto,pnu,pnu,am)
5179 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5180C
5181C K meson (K+ is 321)
5182 CALL tralo4 (kto,pkk,pkk,am)
5183 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5184C
5185 RETURN
5186 END
5187 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5188 COMMON / taukle / bra1,brk0,brk0b,brks
5189 real*4 bra1,brk0,brk0b,brks
5190C ----------------------------------------------------------------------
5191C Lorentz transformation to CMsystem and
5192C Updating of HEPEVT record
5193C
5194C ISGN = 1/-1 for tau-/tau+
5195C JKST=10 (20) corresponds to K0B pi- (K- pi0) decay
5196C
5197C ----------------------------------------------------------------------
5198C
5199 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5200 COMMON /taupos/ np1,np2
5201C
5202C position of decaying particle
5203 IF(kto.EQ. 1) THEN
5204 nps=np1
5205 ELSE
5206 nps=np2
5207 ENDIF
5208C
5209C tau neutrino (nu_tau is 16)
5210 CALL tralo4(kto,pnu,pnu,am)
5211 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5212C
5213C charged K* meson (K*+ is 323)
5214 CALL tralo4(kto,pks,pks,am)
5215 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5216C
5217C two possible decay modes of charged K*
5218 IF(jkst.EQ.10) THEN
5219C
5220C K*- --> pi- K0B (or charged conjugate)
5221C
5222C charged pi meson (pi+ is 211)
5223 CALL tralo4(kto,ppi,ppi,am)
5224 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5225C
5226 bran=brk0b
5227 IF (isgn.EQ.-1) bran=brk0
5228C K0 --> K0_long (is 130) / K0_short (is 310) = 1/1
5229 CALL ranmar(xio,1)
5230 IF(xio(1).GT.bran) THEN
5231 k0type = 130
5232 ELSE
5233 k0type = 310
5234 ENDIF
5235C
5236 CALL tralo4(kto,pkk,pkk,am)
5237 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5238C
5239 ELSE IF(jkst.EQ.20) THEN
5240C
5241C K*- --> pi0 K-
5242C
5243C pi zero (pi0 is 111)
5244 CALL tralo4(kto,ppi,ppi,am)
5245 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5246C
5247C charged K meson (K+ is 321)
5248 CALL tralo4(kto,pkk,pkk,am)
5249 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5250C
5251 ENDIF
5252C
5253 RETURN
5254 END
5255 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5256C ----------------------------------------------------------------------
5257C Lorentz transformation to CMsystem and
5258C Updating of HEPEVT record
5259C
5260C ISGN = 1/-1 for tau-/tau+
5261C
5262C called by : DEXAY,(DEKAY1,DEKAY2)
5263C ----------------------------------------------------------------------
5264C
5265 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5266 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5267 & ,names
5268 COMMON /taupos/ np1,np2
5269 CHARACTER NAMES(NMODE)*31
5270 REAL PNU(4),PWB(4),PNPI(4,9)
5271 REAL PPI(4)
5272C
5273 jnpi=mode-7
5274C position of decaying particle
5275 IF(kto.EQ. 1) THEN
5276 nps=np1
5277 ELSE
5278 nps=np2
5279 ENDIF
5280C
5281C tau neutrino (nu_tau is 16)
5282 CALL tralo4(kto,pnu,pnu,am)
5283 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5284C
5285C W boson (W+ is 24)
5286 CALL tralo4(kto,pwb,pwb,am)
5287 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5288C
5289C multi pi mode JNPI
5290C
5291C get multiplicity of mode JNPI
5292 nd=mulpik(jnpi)
5293 DO i=1,nd
5294 kfpi=lunpik(idffin(i,jnpi),-isgn)
5295C for charged conjugate case, change charged pions only
5296C IF(KFPI.NE.111)KFPI=KFPI*ISGN
5297 DO j=1,4
5298 ppi(j)=pnpi(j,i)
5299 END DO
5300 CALL tralo4(kto,ppi,ppi,am)
5301 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5302 END DO
5303C
5304 RETURN
5305 END
5306 FUNCTION amast(PP)
5307C ----------------------------------------------------------------------
5308C CALCULATES MASS OF PP (DOUBLE PRECISION)
5309C
5310C USED BY : RADKOR
5311C ----------------------------------------------------------------------
5312 IMPLICIT REAL*8 (a-h,o-z)
5313 real*8 pp(4)
5314 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5315C
5316 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5317 amast=aaa
5318 RETURN
5319 END
5320 FUNCTION amas4(PP)
5321C ******************
5322C ----------------------------------------------------------------------
5323C CALCULATES MASS OF PP
5324C
5325C USED BY :
5326C ----------------------------------------------------------------------
5327 REAL PP(4)
5328 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5329 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5330 amas4=aaa
5331 RETURN
5332 END
5333 FUNCTION angxy(X,Y)
5334C ----------------------------------------------------------------------
5335C
5336C USED BY : KORALZ RADKOR
5337C ----------------------------------------------------------------------
5338 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5339 DATA pi /3.141592653589793238462643d0/
5340C
5341 IF(abs(y).LT.abs(x)) THEN
5342 the=atan(abs(y/x))
5343 IF(x.LE.0d0) the=pi-the
5344 ELSE
5345 the=acos(x/sqrt(x**2+y**2))
5346 ENDIF
5347 angxy=the
5348 RETURN
5349 END
5350 FUNCTION angfi(X,Y)
5351C ----------------------------------------------------------------------
5352* CALCULATES ANGLE IN (0,2*PI) RANGE OUT OF X-Y
5353C
5354C USED BY : KORALZ RADKOR
5355C ----------------------------------------------------------------------
5356 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5357 DATA pi /3.141592653589793238462643d0/
5358C
5359 IF(abs(y).LT.abs(x)) THEN
5360 the=atan(abs(y/x))
5361 IF(x.LE.0d0) the=pi-the
5362 ELSE
5363 the=acos(x/sqrt(x**2+y**2))
5364 ENDIF
5365 IF(y.LT.0d0) the=2d0*pi-the
5366 angfi=the
5367 END
5368 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5369C ----------------------------------------------------------------------
5370C
5371C USED BY : KORALZ
5372C ----------------------------------------------------------------------
5373 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5374 dimension pvec(4),qvec(4),rvec(4)
5375C
5376 phi=ph1
5377 cs=cos(phi)
5378 sn=sin(phi)
5379 DO 10 i=1,4
5380 10 rvec(i)=pvec(i)
5381 qvec(1)=rvec(1)
5382 qvec(2)= cs*rvec(2)-sn*rvec(3)
5383 qvec(3)= sn*rvec(2)+cs*rvec(3)
5384 qvec(4)=rvec(4)
5385 RETURN
5386 END
5387 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5388C ----------------------------------------------------------------------
5389C
5390C USED BY : KORALZ RADKOR
5391C ----------------------------------------------------------------------
5392 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5393 dimension pvec(4),qvec(4),rvec(4)
5394C
5395 phi=ph1
5396 cs=cos(phi)
5397 sn=sin(phi)
5398 DO 10 i=1,4
5399 10 rvec(i)=pvec(i)
5400 qvec(1)= cs*rvec(1)+sn*rvec(3)
5401 qvec(2)=rvec(2)
5402 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5403 qvec(4)=rvec(4)
5404 RETURN
5405 END
5406 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5407C ----------------------------------------------------------------------
5408C
5409C USED BY : KORALZ RADKOR
5410C ----------------------------------------------------------------------
5411 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5412C
5413 dimension pvec(4),qvec(4),rvec(4)
5414 phi=ph1
5415 cs=cos(phi)
5416 sn=sin(phi)
5417 DO 10 i=1,4
5418 10 rvec(i)=pvec(i)
5419 qvec(1)= cs*rvec(1)-sn*rvec(2)
5420 qvec(2)= sn*rvec(1)+cs*rvec(2)
5421 qvec(3)=rvec(3)
5422 qvec(4)=rvec(4)
5423 END
5424 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5425C ----------------------------------------------------------------------
5426C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
5427C
5428C USED BY : TAUOLA KORALZ (?)
5429C ----------------------------------------------------------------------
5430 real*4 pvec(4),qvec(4),rvec(4)
5431C
5432 DO 10 i=1,4
5433 10 rvec(i)=pvec(i)
5434 rpl=rvec(4)+rvec(3)
5435 rmi=rvec(4)-rvec(3)
5436 qpl=rpl*exe
5437 qmi=rmi/exe
5438 qvec(1)=rvec(1)
5439 qvec(2)=rvec(2)
5440 qvec(3)=(qpl-qmi)/2
5441 qvec(4)=(qpl+qmi)/2
5442 END
5443 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5444C ----------------------------------------------------------------------
5445C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
5446C
5447C USED BY : KORALZ RADKOR
5448C ----------------------------------------------------------------------
5449 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5450 dimension pvec(4),qvec(4),rvec(4)
5451C
5452 DO 10 i=1,4
5453 10 rvec(i)=pvec(i)
5454 rpl=rvec(4)+rvec(3)
5455 rmi=rvec(4)-rvec(3)
5456 qpl=rpl*exe
5457 qmi=rmi/exe
5458 qvec(1)=rvec(1)
5459 qvec(2)=rvec(2)
5460 qvec(3)=(qpl-qmi)/2
5461 qvec(4)=(qpl+qmi)/2
5462 RETURN
5463 END
5464 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5465C ----------------------------------------------------------------------
5466C
5467C called by :
5468C ----------------------------------------------------------------------
5469 real*4 pvec(4),qvec(4),rvec(4)
5470C
5471 phi=ph1
5472 cs=cos(phi)
5473 sn=sin(phi)
5474 DO 10 i=1,4
5475 10 rvec(i)=pvec(i)
5476 qvec(1)=rvec(1)
5477 qvec(2)= cs*rvec(2)-sn*rvec(3)
5478 qvec(3)= sn*rvec(2)+cs*rvec(3)
5479 qvec(4)=rvec(4)
5480 END
5481 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5482C ----------------------------------------------------------------------
5483C
5484C USED BY : TAUOLA
5485C ----------------------------------------------------------------------
5486 IMPLICIT REAL*4(a-h,o-z)
5487 real*4 pvec(4),qvec(4),rvec(4)
5488C
5489 phi=ph1
5490 cs=cos(phi)
5491 sn=sin(phi)
5492 DO 10 i=1,4
5493 10 rvec(i)=pvec(i)
5494 qvec(1)= cs*rvec(1)+sn*rvec(3)
5495 qvec(2)=rvec(2)
5496 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5497 qvec(4)=rvec(4)
5498 END
5499 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5500C ----------------------------------------------------------------------
5501C
5502C USED BY : TAUOLA
5503C ----------------------------------------------------------------------
5504 real*4 pvec(4),qvec(4),rvec(4)
5505C
5506 cs=cos(phi)
5507 sn=sin(phi)
5508 DO 10 i=1,4
5509 10 rvec(i)=pvec(i)
5510 qvec(1)= cs*rvec(1)-sn*rvec(2)
5511 qvec(2)= sn*rvec(1)+cs*rvec(2)
5512 qvec(3)=rvec(3)
5513 qvec(4)=rvec(4)
5514 END
5515 SUBROUTINE spherd(R,X)
5516C ----------------------------------------------------------------------
5517C GENERATES UNIFORMLY THREE-VECTOR X ON SPHERE OF RADIUS R
5518C DOUBLE PRECISON VERSION OF SPHERA
5519C ----------------------------------------------------------------------
5520 real*8 r,x(4),pi,costh,sinth
5521 real*4 rrr(2)
5522 DATA pi /3.141592653589793238462643d0/
5523C
5524 CALL ranmar(rrr,2)
5525 costh=-1+2*rrr(1)
5526 sinth=sqrt(1 -costh**2)
5527 x(1)=r*sinth*cos(2*pi*rrr(2))
5528 x(2)=r*sinth*sin(2*pi*rrr(2))
5529 x(3)=r*costh
5530 RETURN
5531 END
5532 SUBROUTINE rotpox(THET,PHI,PP)
5533 IMPLICIT REAL*8 (a-h,o-z)
5534C ----------------------------------------------------------------------
5535C
5536C ----------------------------------------------------------------------
5537 dimension pp(4)
5538C
5539 CALL rotod2(thet,pp,pp)
5540 CALL rotod3( phi,pp,pp)
5541 RETURN
5542 END
5543 SUBROUTINE sphera(R,X)
5544C ----------------------------------------------------------------------
5545C GENERATES UNIFORMLY THREE-VECTOR X ON SPHERE OF RADIUS R
5546C
5547C called by : DPHSxx,DADMPI,DADMKK
5548C ----------------------------------------------------------------------
5549 REAL X(4)
5550 real*4 rrr(2)
5551 DATA pi /3.141592653589793238462643/
5552C
5553 CALL ranmar(rrr,2)
5554 costh=-1.+2.*rrr(1)
5555 sinth=sqrt(1.-costh**2)
5556 x(1)=r*sinth*cos(2*pi*rrr(2))
5557 x(2)=r*sinth*sin(2*pi*rrr(2))
5558 x(3)=r*costh
5559 RETURN
5560 END
5561 SUBROUTINE rotpol(THET,PHI,PP)
5562C ----------------------------------------------------------------------
5563C
5564C called by : DADMAA,DPHSAA
5565C ----------------------------------------------------------------------
5566 REAL PP(4)
5567C
5568 CALL rotor2(thet,pp,pp)
5569 CALL rotor3( phi,pp,pp)
5570 RETURN
5571 END
5572 SUBROUTINE ranmar(RVEC,LENV)
5573C ----------------------------------------------------------------------
5574C<<<<<FUNCTION RANMAR(IDUMM)
5575C CERNLIB V113, VERSION WITH AUTOMATIC DEFAULT INITIALIZATION
5576C Transformed to SUBROUTINE to be as in CERNLIB
5577C AM.Lutz November 1988, Feb. 1989
5578C
5579C!Universal random number generator proposed by Marsaglia and Zaman
5580C in report FSU-SCRI-87-50
5581C modified by F. James, 1988 and 1989, to generate a vector
5582C of pseudorandom numbers RVEC of length LENV, and to put in
5583C the COMMON block everything needed to specify currrent state,
5584C and to add input and output entry points RMARIN, RMARUT.
5585C
5586C Unique random number used in the program
5587C ----------------------------------------------------------------------
5588 COMMON / inout / inut,iout
5589 dimension rvec(*)
5590 common/raset1/u(97),c,i97,j97
5591 parameter(modcns=1000000000)
5592 DATA ntot,ntot2,ijkl/-1,0,0/
5593C
5594 IF (ntot .GE. 0) GO TO 50
5595C
5596C Default initialization. User has called RANMAR without RMARIN.
5597 ijkl = 54217137
5598 ntot = 0
5599 ntot2 = 0
5600 kalled = 0
5601 GO TO 1
5602C
5603 entry rmarin(ijklin, ntotin,ntot2n)
5604C Initializing routine for RANMAR, may be called before
5605C generating pseudorandom numbers with RANMAR. The input
5606C values should be in the ranges: 0<=IJKLIN<=900 OOO OOO
5607C 0<=NTOTIN<=999 999 999
5608C 0<=NTOT2N<<999 999 999!
5609C To get the standard values in Marsaglia-s paper, IJKLIN=54217137
5610C NTOTIN,NTOT2N=0
5611 ijkl = ijklin
5612 ntot = max(ntotin,0)
5613 ntot2= max(ntot2n,0)
5614 kalled = 1
5615C always come here to initialize
5616 1 CONTINUE
5617 ij = ijkl/30082
5618 kl = ijkl - 30082*ij
5619 i = mod(ij/177, 177) + 2
5620 j = mod(ij, 177) + 2
5621 k = mod(kl/169, 178) + 1
5622 l = mod(kl, 169)
5623 WRITE(iout,201) ijkl,ntot,ntot2
5624 201 FORMAT(1x,' RANMAR INITIALIZED: ',i10,2x,2i10)
5625 DO 2 ii= 1, 97
5626 s = 0.
5627 t = .5
5628 DO 3 jj= 1, 24
5629 m = mod(mod(i*j,179)*k, 179)
5630 i = j
5631 j = k
5632 k = m
5633 l = mod(53*l+1, 169)
5634 IF (mod(l*m,64) .GE. 32) s = s+t
5635 3 t = 0.5*t
5636 2 u(ii) = s
5637 twom24 = 1.0
5638 DO 4 i24= 1, 24
5639 4 twom24 = 0.5*twom24
5640 c = 362436.*twom24
5641 cd = 7654321.*twom24
5642 cm = 16777213.*twom24
5643 i97 = 97
5644 j97 = 33
5645C Complete initialization by skipping
5646C (NTOT2*MODCNS + NTOT) random numbers
5647 DO 45 loop2= 1, ntot2+1
5648 now = modcns
5649 IF (loop2 .EQ. ntot2+1) now=ntot
5650 IF (now .GT. 0) THEN
5651 WRITE (iout,'(A,I15)') ' RMARIN SKIPPING OVER ',now
5652 DO 40 idum = 1, ntot
5653 uni = u(i97)-u(j97)
5654 IF (uni .LT. 0.) uni=uni+1.
5655 u(i97) = uni
5656 i97 = i97-1
5657 IF (i97 .EQ. 0) i97=97
5658 j97 = j97-1
5659 IF (j97 .EQ. 0) j97=97
5660 c = c - cd
5661 IF (c .LT. 0.) c=c+cm
5662 40 CONTINUE
5663 ENDIF
5664 45 CONTINUE
5665 IF (kalled .EQ. 1) RETURN
5666C
5667C Normal entry to generate LENV random numbers
5668 50 CONTINUE
5669 DO 100 ivec= 1, lenv
5670 uni = u(i97)-u(j97)
5671 IF (uni .LT. 0.) uni=uni+1.
5672 u(i97) = uni
5673 i97 = i97-1
5674 IF (i97 .EQ. 0) i97=97
5675 j97 = j97-1
5676 IF (j97 .EQ. 0) j97=97
5677 c = c - cd
5678 IF (c .LT. 0.) c=c+cm
5679 uni = uni-c
5680 IF (uni .LT. 0.) uni=uni+1.
5681C Replace exact zeroes by uniform distr. *2**-24
5682 IF (uni .EQ. 0.) THEN
5683 uni = twom24*u(2)
5684C An exact zero here is very unlikely, but lets be safe.
5685 IF (uni .EQ. 0.) uni= twom24*twom24
5686 ENDIF
5687 rvec(ivec) = uni
5688 100 CONTINUE
5689 ntot = ntot + lenv
5690 IF (ntot .GE. modcns) THEN
5691 ntot2 = ntot2 + 1
5692 ntot = ntot - modcns
5693 ENDIF
5694 RETURN
5695C Entry to output current status
5696 entry rmarut(ijklut,ntotut,ntot2t)
5697 ijklut = ijkl
5698 ntotut = ntot
5699 ntot2t = ntot2
5700 RETURN
5701 END
5702 FUNCTION dilogt(X)
5703C *****************
5704 IMPLICIT REAL*8(a-h,o-z)
5705CERN C304 VERSION 29/07/71 DILOG 59 C
5706 z=-1.64493406684822
5707 IF(x .LT.-1.0) GO TO 1
5708 IF(x .LE. 0.5) GO TO 2
5709 IF(x .EQ. 1.0) GO TO 3
5710 IF(x .LE. 2.0) GO TO 4
5711 z=3.2898681336964
5712 1 t=1.0/x
5713 s=-0.5
5714 z=z-0.5* log(abs(x))**2
5715 GO TO 5
5716 2 t=x
5717 s=0.5
5718 z=0.
5719 GO TO 5
5720 3 dilogt=1.64493406684822
5721 RETURN
5722 4 t=1.0-x
5723 s=-0.5
5724 z=1.64493406684822 - log(x)* log(abs(t))
5725 5 y=2.66666666666666 *t+0.66666666666666
5726 b= 0.00000 00000 00001
5727 a=y*b +0.00000 00000 00004
5728 b=y*a-b+0.00000 00000 00011
5729 a=y*b-a+0.00000 00000 00037
5730 b=y*a-b+0.00000 00000 00121
5731 a=y*b-a+0.00000 00000 00398
5732 b=y*a-b+0.00000 00000 01312
5733 a=y*b-a+0.00000 00000 04342
5734 b=y*a-b+0.00000 00000 14437
5735 a=y*b-a+0.00000 00000 48274
5736 b=y*a-b+0.00000 00001 62421
5737 a=y*b-a+0.00000 00005 50291
5738 b=y*a-b+0.00000 00018 79117
5739 a=y*b-a+0.00000 00064 74338
5740 b=y*a-b+0.00000 00225 36705
5741 a=y*b-a+0.00000 00793 87055
5742 b=y*a-b+0.00000 02835 75385
5743 a=y*b-a+0.00000 10299 04264
5744 b=y*a-b+0.00000 38163 29463
5745 a=y*b-a+0.00001 44963 00557
5746 b=y*a-b+0.00005 68178 22718
5747 a=y*b-a+0.00023 20021 96094
5748 b=y*a-b+0.00100 16274 96164
5749 a=y*b-a+0.00468 63619 59447
5750 b=y*a-b+0.02487 93229 24228
5751 a=y*b-a+0.16607 30329 27855
5752 a=y*a-b+1.93506 43008 6996
5753 dilogt=s*t*(a-b)+z
5754 RETURN
5755C=======================================================================
5756C===================END OF CPC PART ====================================
5757C=======================================================================
5758 END