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