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