c----------------------------------------------------------------------- subroutine bjinta(ier) c----------------------------------------------------------------------- c fin. state interactions and decays c----------------------------------------------------------------------- include 'epos.inc' double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctimel/ntc common/col3/ncol,kolpt double precision ttaun,ttau0,rcproj,rctarg common/cttaun/ttaun /cttau0/ttau0 /geom1/rcproj,rctarg logical go,lclean call utpri('bjinta',ish,ishini,4) ier=0 if(ncol.eq.0.and.iappl.eq.2)goto1000 if(nevt.ne.1.or.ifrade.eq.0)goto1000 if(iappl.eq.4.or.iappl.eq.9)then goto5000 endif !if(iappl.eq.1)then ! tauxx=0.7+0.94*max(radnuc(maproj),radnuc(matarg))/(0.5*engy) !else ! tauxx=0. !endif !tauzz=max(taumin,tauxx) !print*,'====',taumin,tauxx,tauzz !ttaus=dble(tauzz) ttaus=taumin ttau0=dsqrt(rcproj*rctarg) call jtauin ! initialize hyperbola if(iappl.ne.1)goto 5000 c no-secondary-interactions or parton-ladder-fusion c ------------------------------------------------- if(iorsce.eq.0.and.iorsdf.eq.0.and.iorshh.eq.0 & .or.iorsdf.eq.3)then if(iorsdf.eq.3)then lclean=.false. if(nclean.gt.0.and.nptl.gt.mxptl/5)then ! if nptl already very big, clean up useless particles in cptl list. !(do not use it when gakstr() is called (some information lost) nptli=maproj+matarg+1 do iii=nptli,nptl go=.true. if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false. if(go.and.mod(istptl(iii),10).ne.0)istptl(iii)=99 enddo nptl0=nptl call utclea(nptli,nptl0) lclean=.true. endif nptlbpo=nptl call jintpo(lclean,iret) !parton-ladder-fusion if(ish.ge.2)call alist('parton-ladder-fusion&',nptlbpo+1,nptl) if(iret.eq.1)goto 1001 endif goto 5000 else stop'bjinta: not supported any more (310305). ' endif 5000 continue nptlbd=nptl call xSpaceTime if(ifrade.eq.0)goto779 !skip decay if(idecay.eq.0)goto779 !skip decay if(ish.ge.2)call alist('final decay&',0,0) if(iappl.eq.4.or.iappl.eq.7.or.iappl.eq.9)then nptli=1 else nptli=maproj+matarg+1 endif np1=nptli 41 np2=nptl nptli=np1 ip=np1-1 do while (ip.lt.np2) ip=ip+1 if(istptl(ip).eq.0)then call hdecas(ip,iret) if(iret.eq.1)goto 1001 if(iret.eq.-1)goto 42 c remove useless particles if not enough space if(nclean.gt.0.and.nptl.gt.mxptl/2)then nnnpt=0 do iii=nptli,ip go=.true. if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false. if(go.and.mod(istptl(iii),10).ne.0)then istptl(iii)=99 nnnpt=nnnpt+1 endif enddo if(nnnpt.gt.mxptl-nptl)then nptl0=nptl call utclea(nptli,nptl0) np2=np2-nnnpt ip=ip-nnnpt nptli=ip endif endif endif 42 continue enddo nptli=max(nptli,np1) np1=np2+1 if(np1.le.nptl)then if(ish.ge.2)then if(ish.ge.3)call alist('partial list&',0,0) do 6 ip=np1,nptl call alist('&',ip,ip) 6 continue endif goto 41 endif 779 continue c if(ish.ge.2)call alist('complete list&',1,nptl) c on shell check c -------------- c if(iappl.eq.1)call jresc 1000 continue call utprix('bjinta',ish,ishini,4) return 1001 continue ier=1 goto 1000 end cc---------------------------------------------------------------------- c subroutine jintcs(i,j,ecm,bij,nq,jc,ics) cc---------------------------------------------------------------------- cc compare hadron distance with energy dependent cross section cc data taken from particle data group, durham and juelich cc input: cc i,j: particle indices cc ecm: center-of-mass energy cc bij: impact parameter cc nq: net quark number of fused object cc jc: jc of fused object cc output: cc ics=0 if distance larger than sqrt(sig(E_CMS)/pi) cc ics=1 else cc The data are from HEPDATA, cc the formulas from Rev. Particle Properties 1995 cc---------------------------------------------------------------------- c include 'epos.inc' c integer jci(nflav,2),jcj(nflav,2),jc(nflav,2),kc(nflav) c *,kci(nflav),kcj(nflav) c common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl) c *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3) c parameter(npp=249,napp=205,npn=411,napn=31,npip=441) c parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91) cc parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91,nlp=35) c parameter(npi1=12,npi2=12,npi3=18,npi4=21,npi5=9) cc real ppecm(npp) c real ppbmx(npp) c real appecm(napp),appbmx(napp) c real pnecm(npn),pnbmx(npn) c real apnecm(napn),apnbmx(napn) c real pipecm(npip),pipbmx(npip) c real pimecm(npim),pimbmx(npim) c real kmpecm(nkmp),kmpbmx(nkmp) c real kmnecm(nkmn),kmnbmx(nkmn) c real kppecm(nkpp),kppbmx(nkpp) c real kpnecm(nkpn),kpnbmx(nkpn) cc real lpecm(nlp),lpbmx(nlp) c real pi1ecm(npi1),pi1bmx(npi1) c real pi2ecm(npi2),pi2bmx(npi2) c real pi3ecm(npi3),pi3bmx(npi3) c real pi4ecm(npi4),pi4bmx(npi4) c real pi5ecm(npi5),pi5bmx(npi5) c cc data ppecm/ cc * 1.8812, 1.8855, 1.8910, 1.8963, 1.9073 cc *, 1.9108, 1.9145, 1.9224, 1.9244, 1.9352 cc *, 1.9466, 1.9468, 1.9542, 1.9592, 1.9636 cc *, 1.9772, 1.9860, 1.9945, 2.0032, 2.0052 cc *, 2.0070, 2.0272, 2.0275, 2.0302, 2.0333 cc *, 2.0402, 2.0427, 2.0586, 2.0608, 2.0692 cc *, 2.0702, 2.0708, 2.0715, 2.0718, 2.0751 cc *, 2.0797, 2.0813, 2.0829, 2.0843, 2.0846 cc *, 2.0935, 2.1062, 2.1113, 2.1123, 2.1170 cc *, 2.1173, 2.1184, 2.1268, 2.1289, 2.1357 cc *, 2.1467, 2.1511, 2.1522, 2.1553, 2.1618 cc *, 2.1639, 2.1726, 2.1771, 2.1795, 2.1799 cc *, 2.1802, 2.1813, 2.1868, 2.2152, 2.2184 cc *, 2.2212, 2.2254, 2.2395, 2.2405, 2.2532 cc *, 2.2606, 2.2613, 2.2861, 2.2889, 2.2914 cc *, 2.2988, 2.3101, 2.3136, 2.3228, 2.3348 cc *, 2.3490, 2.3525, 2.3529, 2.3800, 2.3842 cc *, 2.3912, 2.4088, 2.4130, 2.4193, 2.4298 cc *, 2.4315, 2.4392, 2.4472, 2.4573, 2.5034 cc *, 2.5131, 2.5268, 2.5743, 2.5848, 2.5916 cc *, 2.6327, 2.6620, 2.6700, 2.6984, 2.7080 cc *, 2.7205, 2.7534, 2.7573, 2.7651, 2.7670 cc *, 2.7844, 2.8024, 2.8092, 2.8127, 2.8533 cc *, 2.8556, 2.8638, 2.9079, 2.9395, 2.9500 cc *, 2.9776, 2.9961, 3.0495, 3.0769, 3.0879 cc *, 3.1358, 3.1547, 3.2251, 3.2371, 3.3020 cc *, 3.3527, 3.3620, 3.4221, 3.4967, 3.5019 cc *, 3.5035, 3.5814, 3.5829, 3.6266, 3.8549 cc *, 3.8742, 4.0503, 4.0663, 4.0698, 4.0732 cc *, 4.0778, 4.1074, 4.1300, 4.4976, 4.5183 cc *, 4.5389, 4.5410, 4.5615, 4.6808, 4.9146 cc *, 4.9336, 5.0088, 5.2993, 5.3011, 5.4731 cc *, 5.6083, 5.6416, 5.6465, 5.9171, 5.9502 cc *, 5.9644, 6.1697, 6.1803, 6.2706, 6.3034 cc *, 6.3390, 6.5627, 6.7040, 6.8424, 6.9105 cc *, 6.9780, 7.1111, 7.1662, 7.6202, 7.8624 cc *, 8.2124, 8.7647, 9.0282, 9.2843, 9.5825 cc *, 9.7763, 9.9851, 10.2447, 10.6927, 10.8926 cc *, 11.4549, 11.5365, 11.7779, 11.8519, 13.6241 cc *, 13.6883, 13.7611, 13.8968, 15.0628, 15.1868 cc *, 16.6595, 16.8275, 17.9077, 18.0121, 18.1677 cc *, 19.2213, 19.4156, 19.6556, 19.7002, 21.2604 cc *, 22.9574, 23.3624, 23.4057, 23.4965, 23.4967 cc *, 23.5964, 23.7605, 23.8787, 24.1521, 25.2904 cc *, 26.3796, 27.5960, 30.5240, 30.5954, 30.5957 cc *, 30.6555, 30.6954, 30.7954, 35.1947, 44.6933 cc *, 44.6937, 44.6938, 44.7706, 44.8228, 44.8933 cc *, 45.1933, 52.6798, 52.7090, 52.7921, 52.7921 cc *, 52.7927, 52.8927, 53.1921, 62.2907, 62.3914 cc *, 62.4907, 62.6907, 62.6913, 62.7906 cc */ c data ppbmx/ c * 3.1615, 2.2212, 1.7113, 1.4927, 1.1631 c *, 1.0911, 1.0388, .9525, .9390, .8885 c *, .8019, .8956, .9115, .8593, .8840 c *, .9062, .8444, .8444, .8482, .8713 c *, .8575, .8795, .8795, .8759, .8630 c *, .8822, .8593, .8813, .9036, .9150 c *, .8740, .9219, .8740, .9253, .8722 c *, .8795, .9288, .8704, .9398, .9271 c *, .9407, .9373, .9756, .9926, 1.0357 c *, 1.0037, 1.0408, .9739, 1.0108, 1.0479 c *, 1.0645, 1.0705, 1.1113, 1.0794, 1.0955 c *, 1.1085, 1.1256, 1.1535, 1.1731, 1.1424 c *, 1.0794, 1.1480, 1.1590, 1.1888, 1.2218 c *, 1.2164, 1.2087, 1.2296, 1.2231, 1.2335 c *, 1.2322, 1.2309, 1.2114, 1.2296, 1.2293 c *, 1.2489, 1.2303, 1.2322, 1.2322, 1.2127 c *, 1.2192, 1.2295, 1.2399, 1.2290, 1.2374 c *, 1.2205, 1.2278, 1.2284, 1.2061, 1.2296 c *, 1.2296, 1.2540, 1.2008, 1.2260, 1.2229 c *, 1.2257, 1.2188, 1.2118, 1.2078, 1.1982 c *, 1.2039, 1.2012, 1.1991, 1.1808, 1.1969 c *, 1.1959, 1.1922, 1.1902, 1.1897, 1.1878 c *, 1.1888, 1.1860, 1.1856, 1.1850, 1.2244 c *, 1.1782, 1.1790, 1.1718, 1.1696, 1.1726 c *, 1.1576, 1.1656, 1.1606, 1.1603, 1.1581 c *, 1.1590, 1.1530, 1.1576, 1.1487, 1.1476 c *, 1.1447, 1.1794, 1.1448, 1.1507, 1.1507 c *, 1.1407, 1.1403, 1.1507, 1.1581, 1.1645 c *, 1.1616, 1.1507, 1.1332, 1.1294, 1.1284 c *, 1.1227, 1.1284, 1.1298, 1.1261, 1.1199 c *, 1.1365, 1.1438, 1.1284, 1.1424, 1.1230 c *, 1.1218, 1.1142, 1.1156, 1.1202, 1.1198 c *, 1.1099, 1.1099, 1.1175, 1.1241, 1.1168 c *, 1.1099, 1.1128, 1.1241, 1.1103, 1.1149 c *, 1.1155, 1.1083, 1.1197, 1.1127, 1.1185 c *, 1.1113, 1.1128, 1.1113, 1.1083, 1.1056 c *, 1.1067, 1.1070, 1.1059, 1.1063, 1.1070 c *, 1.1028, 1.0979, 1.1060, 1.1062, 1.1774 c *, 1.0805, 1.1039, 1.0940, 1.0955, 1.0940 c *, 1.0955, 1.1082, 1.1128, 1.1082, 1.0852 c *, 1.1003, 1.1097, 1.1118, 1.0911, 1.1227 c *, 1.1070, 1.1206, 1.1142, 1.1213, 1.1174 c *, 1.1199, 1.1142, 1.1185, 1.1180, 1.1113 c *, 1.1099, 1.1297, 1.1142, 1.1226, 1.1240 c *, 1.1251, 1.1368, 1.1403, 1.1319, 1.1296 c *, 1.1333, 1.1303, 1.1284, 1.1343, 1.1516 c *, 1.1521, 1.1549, 1.1641, 1.1631, 1.1562 c *, 1.1631, 1.1697, 1.1726, 1.1756, 1.1659 c *, 1.1660, 1.1617, 1.1686, 1.1774, 1.1713 c *, 1.1744, 1.1810, 1.1694, 1.1848 c */ c data appecm/ c * 1.9002, 1.9050, 1.9072, 1.9078, 1.9091 c *, 1.9129, 1.9157, 1.9162, 1.9174, 1.9176 c *, 1.9180, 1.9195, 1.9201, 1.9224, 1.9226 c *, 1.9246, 1.9252, 1.9255, 1.9257, 1.9271 c *, 1.9282, 1.9293, 1.9301, 1.9310, 1.9319 c *, 1.9328, 1.9334, 1.9345, 1.9359, 1.9370 c *, 1.9372, 1.9384, 1.9393, 1.9398, 1.9407 c *, 1.9426, 1.9430, 1.9433, 1.9452, 1.9454 c *, 1.9473, 1.9485, 1.9495, 1.9500, 1.9510 c *, 1.9515, 1.9547, 1.9559, 1.9562, 1.9579 c *, 1.9610, 1.9615, 1.9644, 1.9680, 1.9718 c *, 1.9755, 1.9788, 1.9829, 1.9871, 1.9911 c *, 1.9954, 1.9994, 2.0813, 2.0979, 2.1146 c *, 2.1180, 2.1316, 2.1487, 2.1660, 2.1834 c *, 2.1868, 2.1938, 2.1991, 2.2008, 2.2184 c *, 2.2226, 2.2359, 2.2500, 2.2606, 2.2712 c *, 2.2889, 2.2995, 2.3066, 2.3243, 2.3419 c *, 2.3490, 2.3511, 2.3596, 2.3701, 2.3772 c *, 2.3860, 2.3877, 2.3947, 2.4035, 2.4123 c *, 2.4298, 2.4472, 2.4629, 2.4820, 2.4993 c *, 2.5165, 2.5268, 2.5337, 2.5508, 2.5678 c *, 2.5848, 2.5950, 2.6017, 2.6186, 2.6353 c *, 2.6377, 2.6520, 2.6654, 2.6687, 2.6852 c *, 2.7017, 2.7182, 2.7345, 2.7508, 2.7670 c *, 2.7832, 2.7896, 2.7992, 2.8152, 2.8312 c *, 2.8439, 2.8470, 2.8565, 2.8628, 2.9377 c *, 2.9561, 2.9745, 3.0351, 3.0769, 3.0828 c *, 3.1648, 3.2507, 3.2788, 3.3620, 3.4568 c *, 3.5492, 3.6266, 3.7893, 3.8501, 3.8597 c *, 3.8742, 3.9455, 4.1074, 4.3499, 4.3926 c *, 4.5389, 4.5799, 4.6808, 4.7991, 4.9336 c *, 4.9901, 5.1742, 5.2993, 5.3520, 5.4731 c *, 5.6416, 5.6911, 5.8534, 5.9644, 6.0113 c *, 6.2706, 6.3153, 6.7040, 6.9780, 7.3061 c *, 7.6202, 7.7422, 7.8624, 7.8743, 8.0393 c *, 8.2124, 8.4930, 8.7647, 9.0282, 9.2843 c *, 9.5335, 9.7763, 11.5365, 13.7611, 13.8630 c *, 15.0628, 16.8275, 17.9077, 19.4156, 21.2604 c *, 22.9574, 30.4098, 30.5943, 30.6861, 52.5843 c *, 52.7979, 52.7979, 62.2853, 62.4957, 62.6905 c *, 539.9198, 546.9191, 899.8658, 900.0000,1803.0007 c */ c data appbmx/ c * 2.7553, 2.6857, 2.6732, 2.6517, 2.6306 c *, 2.5376, 2.4398, 2.4779, 2.5212, 2.5143 c *, 2.4573, 2.4799, 2.4495, 2.4286, 2.4521 c *, 2.3683, 2.4463, 2.4240, 2.4489, 2.3743 c *, 2.4162, 2.3561, 2.3903, 2.3796, 2.3221 c *, 2.3689, 2.4181, 2.3221, 2.3823, 2.2987 c *, 2.3337, 2.3568, 2.3487, 2.2911, 2.3344 c *, 2.2701, 2.3097, 2.3207, 2.2631, 2.2603 c *, 2.2666, 2.2412, 2.2319, 2.3097, 2.2327 c *, 2.2104, 2.2104, 2.2645, 2.2369, 2.1880 c *, 2.1756, 2.2162, 2.1491, 2.1484, 2.1365 c *, 2.1320, 2.1027, 2.1065, 2.0913, 2.0776 c *, 2.0668, 2.0883, 1.9333, 1.9098, 1.8851 c *, 1.8932, 1.8721, 1.8623, 1.8520, 1.8409 c *, 1.8678, 1.8575, 1.8325, 1.8429, 1.8088 c *, 1.7698, 1.7941, 1.7864, 1.7814, 1.7736 c *, 1.7645, 1.7623, 1.7576, 1.7523, 1.7445 c *, 1.7419, 1.7019, 1.7342, 1.7311, 1.7271 c *, 1.7210, 1.7182, 1.7161, 1.7119, 1.7054 c *, 1.6947, 1.6816, 1.6780, 1.6678, 1.6599 c *, 1.6509, 1.6328, 1.6449, 1.6396, 1.6319 c *, 1.6270, 1.5761, 1.6043, 1.6120, 1.6069 c *, 1.5978, 1.6018, 1.5891, 1.5948, 1.5905 c *, 1.5849, 1.5778, 1.5736, 1.5679, 1.5785 c *, 1.5583, 1.5476, 1.5519, 1.5467, 1.5416 c *, 1.5233, 1.5368, 1.5492, 1.5313, 1.4895 c *, 1.5579, 1.5107, 1.4680, 1.4725, 1.4586 c *, 1.3889, 1.4799, 1.4472, 1.4371, 1.3576 c *, 1.4228, 1.3922, 1.3762, 1.3854, 1.3669 c *, 1.4161, 1.3623, 1.3638, 1.3530, 1.3762 c *, 1.3273, 1.2866, 1.2989, 1.3159, 1.2828 c *, 1.3086, 1.2797, 1.2704, 1.2940, 1.2741 c *, 1.2514, 1.2540, 1.2803, 1.2653, 1.2603 c *, 1.2386, 1.2192, 1.2283, 1.2231, 1.2140 c *, 1.2048, 1.2114, 1.2099, 1.2052, 1.2048 c *, 1.1984, 1.1928, 1.1901, 1.1902, 1.1888 c *, 1.1848, 1.1769, 1.1706, 1.1579, 1.1608 c *, 1.1521, 1.1534, 1.1520, 1.1495, 1.1549 c *, 1.1550, 1.1580, 1.1672, 1.1562, 1.1743 c *, 1.1859, 1.1950, 1.1851, 1.1821, 1.1987 c *, 1.4740, 1.4037, 1.4483, 1.4417, 1.5149 c */ c data pnecm/ c * 1.87867, 1.87869, 1.87872, 1.87875, 1.87878 c *, 1.87881, 1.87884, 1.87887, 1.87890, 1.87893 c *, 1.87896, 1.87899, 1.87902, 1.87906, 1.87909 c *, 1.87913, 1.87916, 1.87920, 1.87923, 1.87927 c *, 1.87931, 1.87934, 1.87938, 1.87942, 1.87946 c *, 1.87950, 1.87954, 1.87958, 1.87962, 1.87966 c *, 1.87970, 1.87975, 1.87979, 1.87983, 1.87988 c *, 1.87992, 1.87997, 1.88001, 1.88006, 1.88011 c *, 1.88015, 1.88020, 1.88025, 1.88030, 1.88035 c *, 1.88040, 1.88045, 1.88050, 1.88055, 1.88061 c *, 1.88066, 1.88071, 1.88077, 1.88082, 1.88087 c *, 1.88093, 1.88099, 1.88104, 1.88110, 1.88116 c *, 1.88121, 1.88127, 1.88133, 1.88139, 1.88145 c *, 1.88151, 1.88157, 1.88163, 1.88170, 1.88176 c *, 1.88182, 1.88189, 1.88195, 1.88202, 1.88208 c *, 1.88215, 1.88221, 1.88228, 1.88235, 1.88241 c *, 1.88248, 1.88255, 1.88262, 1.88269, 1.88276 c *, 1.88283, 1.88290, 1.88297, 1.88305, 1.88312 c *, 1.88319, 1.88327, 1.88334, 1.88342, 1.88349 c *, 1.88357, 1.88364, 1.88372, 1.88380, 1.88388 c *, 1.88396, 1.88403, 1.88411, 1.88419, 1.88428 c *, 1.88436, 1.88444, 1.88452, 1.88460, 1.88469 c *, 1.88477, 1.88485, 1.88494, 1.88502, 1.88511 c *, 1.88519, 1.88528, 1.88537, 1.88546, 1.88554 c *, 1.88563, 1.88572, 1.88581, 1.88590, 1.88599 c *, 1.88608, 1.88618, 1.88627, 1.88636, 1.88645 c *, 1.88655, 1.88664, 1.88674, 1.88683, 1.88693 c *, 1.88702, 1.88712, 1.88722, 1.88731, 1.88741 c *, 1.88751, 1.88761, 1.88771, 1.88781, 1.88791 c *, 1.88801, 1.88811, 1.88822, 1.88832, 1.88842 c *, 1.88852, 1.88863, 1.88873, 1.88884, 1.88894 c *, 1.88905, 1.88916, 1.88926, 1.88937, 1.88948 c *, 1.88959, 1.88970, 1.88980, 1.88991, 1.89003 c *, 1.89014, 1.89025, 1.89036, 1.89047, 1.89058 c *, 1.89070, 1.89081, 1.89093, 1.89104, 1.89116 c *, 1.89127, 1.89139, 1.89150, 1.89162, 1.89174 c *, 1.89186, 1.89198, 1.89209, 1.89221, 1.89233 c *, 1.89245, 1.89258, 1.89270, 1.89282, 1.89294 c *, 1.89306, 1.89319, 1.89331, 1.89344, 1.89356 c *, 1.89369, 1.89381, 1.89394, 1.89406, 1.89419 c *, 1.89432, 1.89445, 1.89458, 1.89496, 1.89549 c *, 1.89602, 1.89615, 1.89656, 1.89697, 1.89724 c *, 1.89766, 1.89780, 1.89850, 1.89893, 1.89922 c *, 1.89995, 1.90068, 1.90098, 1.90143, 1.90174 c *, 1.90235, 1.90250, 1.90312, 1.90407, 1.90503 c *, 1.90519, 1.90616, 1.90732, 1.90749, 1.90833 c *, 1.90936, 1.91216, 1.91379, 1.91545, 1.91714 c *, 1.91905, 1.91924, 1.92100, 1.92160, 1.92179 c *, 1.92259, 1.92319, 1.92421, 1.92502, 1.92543 c *, 1.92605, 1.92646, 1.92792, 1.93046, 1.93241 c *, 1.93306, 1.93570, 1.93704, 1.93953, 1.94021 c *, 1.94183, 1.94699, 1.94723, 1.95206, 1.95329 c *, 1.95452, 1.95651, 1.96029, 1.96080, 1.97927 c *, 1.98533, 1.99806, 2.00884, 2.01269, 2.02779 c *, 2.02963, 2.04329, 2.04643, 2.05915, 2.05947 c *, 2.05979, 2.07207, 2.07337, 2.07533, 2.09178 c *, 2.10847, 2.11385, 2.12061, 2.12536, 2.14243 c *, 2.15308, 2.15343, 2.15964, 2.17072, 2.17698 c *, 2.18185, 2.19442, 2.21193, 2.21369, 2.22353 c *, 2.22846, 2.22952, 2.24715, 2.25597, 2.28851 c *, 2.29134, 2.29382, 2.31293, 2.31328, 2.32673 c *, 2.34935, 2.35501, 2.36207, 2.38251, 2.39729 c *, 2.41135, 2.41591, 2.42817, 2.44006, 2.45995 c *, 2.48220, 2.49952, 2.50643, 2.52951, 2.56750 c *, 2.58756, 2.63546, 2.64719, 2.66487, 2.67285 c *, 2.71088, 2.72336, 2.75178, 2.75634, 2.76802 c *, 2.76867, 2.76996, 2.78741, 2.80542, 2.81567 c *, 2.85639, 2.85860, 2.86681, 2.91102, 2.94265 c *, 2.96543, 3.05273, 3.08729, 3.09115, 3.15805 c *, 3.16821, 3.22857, 3.24053, 3.29551, 3.35628 c *, 3.42580, 3.50727, 3.55297, 3.58520, 3.58675 c *, 3.59580, 3.63049, 3.72425, 3.75633, 4.05460 c *, 4.06219, 4.07412, 4.11175, 4.37426, 4.52314 c *, 4.54378, 4.74537, 4.93885, 5.12513, 5.30495 c *, 5.40999, 5.47892, 5.61424, 5.64757, 5.93444 c *, 5.97071, 6.27732, 6.42517, 6.51227, 6.56970 c *, 6.71114, 6.87703, 6.98545, 7.18446, 7.23904 c *, 7.24942, 7.62830, 8.10603, 8.22114, 8.54948 c *, 8.77406, 9.29418, 9.78673, 10.15712, 10.25566 c *, 10.70407, 11.13446, 11.54882, 12.33587, 13.77577 c *, 15.07880, 15.74960, 16.84544, 17.92675, 18.18703 c *, 18.44365, 19.43624, 20.14863, 21.28302, 22.69375 c *, 22.98187 c */ c data pnbmx/ c * 10.7613, 10.6195, 10.5290, 10.4010, 10.2901 c *, 10.1970, 10.1029, 10.0088, 9.9009, 9.7973 c *, 9.7125, 9.6198, 9.5613, 9.4281, 9.3628 c *, 9.2801, 9.2068, 9.1103, 9.0239, 8.9609 c *, 8.8562, 8.8011, 8.7574, 8.6332, 8.5487 c *, 8.4788, 8.4157, 8.3606, 8.2746, 8.2604 c *, 8.1583, 8.0778, 8.0127, 7.9567, 7.8691 c *, 7.8176, 7.8190, 7.7028, 7.6634, 7.5956 c *, 7.5808, 7.4956, 7.4144, 7.3600, 7.3071 c *, 7.2947, 7.1947, 7.1389, 7.1333, 7.0404 c *, 6.9856, 6.9468, 6.9171, 6.8369, 6.7869 c *, 6.7829, 6.6804, 6.6483, 6.5982, 6.5464 c *, 6.5060, 6.4541, 6.4539, 6.3684, 6.3219 c *, 6.2777, 6.2832, 6.1881, 6.1494, 6.1048 c *, 6.0642, 6.0207, 6.0255, 5.9756, 5.9161 c *, 5.8717, 5.8276, 5.7878, 5.7481, 5.7319 c *, 5.6977, 5.6492, 5.6142, 5.5796, 5.5428 c *, 5.4939, 5.4659, 5.4549, 5.3933, 5.3601 c *, 5.3311, 5.3037, 5.2760, 5.2429, 5.2381 c *, 5.2384, 5.1470, 5.1112, 5.0812, 5.0483 c *, 5.0350, 5.0330, 4.9640, 4.9410, 4.9050 c *, 4.8714, 4.8408, 4.8482, 4.8518, 4.7611 c *, 4.7157, 4.7049, 4.6834, 4.6939, 4.6206 c *, 4.5931, 4.5723, 4.5498, 4.5210, 4.5337 c *, 4.4726, 4.4898, 4.4251, 4.4013, 4.3784 c *, 4.3514, 4.3264, 4.3530, 4.2846, 4.2567 c *, 4.2339, 4.2171, 4.1919, 4.1757, 4.1503 c *, 4.1329, 4.0924, 4.0797, 4.0624, 4.0446 c *, 4.0210, 4.0010, 3.9812, 3.9627, 3.9402 c *, 3.9221, 3.9029, 3.8925, 3.8680, 3.8439 c *, 3.8244, 3.8094, 3.7816, 3.7719, 3.7396 c *, 3.7271, 3.7011, 3.6849, 3.6804, 3.6640 c *, 3.6496, 3.5954, 3.6042, 3.5697, 3.5552 c *, 3.5463, 3.5265, 3.5212, 3.4992, 3.1432 c *, 3.4630, 3.4483, 3.4484, 3.4240, 3.4127 c *, 3.3595, 3.3767, 3.3582, 3.3059, 3.3258 c *, 3.3182, 3.2582, 3.2740, 3.2621, 3.1984 c *, 3.2408, 3.2043, 3.2087, 3.1652, 3.1720 c *, 3.1476, 3.0754, 3.1236, 3.1494, 3.0653 c *, 3.0757, 3.1030, 3.0692, 3.0466, 3.0392 c *, 3.0346, 2.9933, 3.0146, 2.8774, 2.7972 c *, 2.6463, 2.8068, 2.6869, 2.6643, 2.6673 c *, 2.6684, 2.6409, 2.5520, 2.4341, 2.4978 c *, 2.4534, 2.4417, 2.4293, 2.3283, 2.4978 c *, 2.3042, 2.3111, 2.1996, 2.1996, 2.1960 c *, 2.1469, 2.1223, 2.0791, 2.0599, 1.9891 c *, 2.0027, 1.8678, 1.7850, 1.7868, 1.7362 c *, 1.6478, 1.6254, 1.6166, 1.6478, 1.6516 c *, 1.5554, 1.5727, 1.5656, 1.5290, 1.5149 c *, 1.5554, 1.5615, 1.5522, 1.4863, 1.4745 c *, 1.4745, 1.3991, 1.4139, 1.1686, 1.3458 c *, 1.3587, 1.2425, 1.3171, 1.2828, 1.2153 c *, 1.2679, 1.2766, 1.1410, 1.2514, 1.0852 c *, 1.1445, 1.1672, 1.0624, 1.1089, 1.0899 c *, 1.0171, 1.0399, 1.0645, 1.0417, .9934 c *, 1.0403, 1.0029, 1.0357, 1.0388, 1.0337 c *, 1.0428, 1.0555, 1.0663, 1.0546, 1.0626 c *, 1.0013, 1.0705, 1.0711, 1.0852, 1.0830 c *, 1.1119, 1.0940, 1.0976, 1.0675, 1.1205 c *, 1.0464, 1.0995, .9508, 1.0972, 1.1170 c *, 1.1015, 1.1251, 1.1296, 1.1027, 1.1028 c *, .9271, 1.1363, 1.1120, 1.1455, 1.1192 c *, 1.1498, 1.1491, 1.0108, 1.1282, 1.1550 c *, 1.1617, 1.1354, 1.1586, 1.1631, 1.1378 c *, 1.1656, 1.1684, 1.1389, 1.1694, 1.1690 c *, 1.1702, 1.1705, 1.1390, 1.1714, 1.1716 c *, 1.1326, 1.1517, 1.1697, 1.1731, 1.1716 c *, 1.0867, 1.1537, 1.1698, 1.1642, 1.1638 c *, 1.1199, 1.1635, 1.1713, 1.1631, 1.1601 c *, 1.1340, 1.0823, 1.1598, 1.1754, 1.1572 c *, 1.1565, 1.1567, 1.1631, 1.1538, 1.0852 c *, 1.0342, 1.1645, 1.1452, 1.1099, 1.0940 c *, 1.1185, 1.1470, 1.1388, 1.1424, 1.0705 c *, 1.1374, 1.1199, 1.1262, 1.1142, 1.1205 c *, 1.0867, 1.1095, 1.0734, 1.1234, 1.0925 c *, 1.1135, 1.1076, 1.1070, 1.0955, 1.1027 c *, 1.1073, 1.0630, 1.1007, 1.1185, 1.1133 c *, 1.1128, 1.1028, 1.1033, 1.1086, 1.1120 c *, 1.1013, 1.0991, 1.1086, 1.1028, 1.1027 c *, 1.1064, 1.1036, 1.1135, 1.1139, 1.1136 c *, 1.1145, 1.1166, 1.1152, 1.1168, 1.1172 c *, 1.1230, 1.1187, 1.1254, 1.1224, 1.1329 c *, 1.1216 c */ c data apnecm/ c * 2.1288, 2.2007, 2.2500, 2.3044, 2.3529 c *, 2.9284, 3.5136, 3.6305, 3.7933, 4.1117 c *, 4.5438, 4.9389, 5.1797, 5.3049, 5.4789 c *, 5.6476, 6.2773, 6.9855, 7.6283, 8.2211 c *, 8.7741, 9.2942, 9.7867, 11.5488, 13.7758 c *, 15.0788, 16.8454, 17.9267, 19.4362, 21.2830 c *, 22.9819 c */ c data apnbmx/ c * 1.9462, 1.7481, 1.8881, 1.8019, 1.8627 c *, 1.4745, 1.3231, 1.3762, 1.3351, 1.3505 c *, 1.2841, 1.3086, 1.2565, 1.3038, 1.2218 c *, 1.2952, 1.2127, 1.2035, 1.1968, 1.1951 c *, 1.1604, 1.1794, 1.1747, 1.1724, 1.1595 c *, 1.1637, 1.1477, 1.1520, 1.1510, 1.1478 c *, 1.1466 c */ c data pipecm/ c * 1.1050, 1.1154, 1.1165, 1.1256, 1.1273 c *, 1.1333, 1.1370, 1.1382, 1.1394, 1.1438 c *, 1.1495, 1.1579, 1.1592, 1.1677, 1.1691 c *, 1.1697, 1.1757, 1.1771, 1.1777, 1.1784 c *, 1.1798, 1.1838, 1.1892, 1.1906, 1.1926 c *, 1.1933, 1.1953, 1.1960, 1.1967, 1.1981 c *, 1.1994, 1.2008, 1.2015, 1.2022, 1.2028 c *, 1.2056, 1.2063, 1.2069, 1.2097, 1.2104 c *, 1.2111, 1.2124, 1.2131, 1.2138, 1.2152 c *, 1.2166, 1.2172, 1.2193, 1.2200, 1.2214 c *, 1.2221, 1.2255, 1.2262, 1.2269, 1.2276 c *, 1.2283, 1.2303, 1.2310, 1.2317, 1.2352 c *, 1.2358, 1.2365, 1.2372, 1.2400, 1.2407 c *, 1.2421, 1.2434, 1.2462, 1.2476, 1.2503 c *, 1.2517, 1.2524, 1.2538, 1.2545, 1.2565 c *, 1.2613, 1.2627, 1.2696, 1.2716, 1.2730 c *, 1.2751, 1.2799, 1.2819, 1.2833, 1.2860 c *, 1.2867, 1.2915, 1.2922, 1.2929, 1.2990 c *, 1.3010, 1.3024, 1.3119, 1.3126, 1.3180 c *, 1.3200, 1.3220, 1.3254, 1.3261, 1.3321 c *, 1.3382, 1.3388, 1.3415, 1.3422, 1.3449 c *, 1.3469, 1.3522, 1.3608, 1.3615, 1.3622 c *, 1.3655, 1.3661, 1.3767, 1.3786, 1.3813 c *, 1.3898, 1.3917, 1.3950, 1.4047, 1.4080 c *, 1.4112, 1.4157, 1.4163, 1.4176, 1.4208 c *, 1.4215, 1.4285, 1.4298, 1.4304, 1.4336 c *, 1.4349, 1.4425, 1.4470, 1.4495, 1.4526 c *, 1.4539, 1.4652, 1.4677, 1.4702, 1.4764 c *, 1.4777, 1.4857, 1.4863, 1.4870, 1.4882 c *, 1.4919, 1.4993, 1.4999, 1.5030, 1.5097 c *, 1.5103, 1.5121, 1.5146, 1.5212, 1.5249 c *, 1.5285, 1.5327, 1.5351, 1.5429, 1.5435 c *, 1.5465, 1.5513, 1.5531, 1.5548, 1.5566 c *, 1.5637, 1.5655, 1.5673, 1.5708, 1.5720 c *, 1.5732, 1.5761, 1.5790, 1.5849, 1.5866 c *, 1.5965, 1.5977, 1.5994, 1.6023, 1.6063 c *, 1.6121, 1.6133, 1.6144, 1.6264, 1.6287 c *, 1.6327, 1.6344, 1.6406, 1.6417, 1.6429 c *, 1.6502, 1.6519, 1.6536, 1.6625, 1.6664 c *, 1.6670, 1.6676, 1.6687, 1.6692, 1.6714 c *, 1.6792, 1.6825, 1.6853, 1.6913, 1.6935 c *, 1.6941, 1.6963, 1.6985, 1.6990, 1.7012 c *, 1.7170, 1.7175, 1.7208, 1.7240, 1.7267 c *, 1.7391, 1.7417, 1.7423, 1.7428, 1.7503 c *, 1.7535, 1.7609, 1.7635, 1.7651, 1.7714 c *, 1.7735, 1.7767, 1.7793, 1.7798, 1.7908 c *, 1.7923, 1.7949, 1.7960, 1.8032, 1.8063 c *, 1.8099, 1.8197, 1.8207, 1.8238, 1.8253 c *, 1.8289, 1.8304, 1.8320, 1.8386, 1.8436 c *, 1.8497, 1.8507, 1.8537, 1.8643, 1.8673 c *, 1.8678, 1.8703, 1.8713, 1.8738, 1.8777 c *, 1.8792, 1.8812, 1.8827, 1.8881, 1.8911 c *, 1.8916, 1.8936, 1.8980, 1.9039, 1.9073 c *, 1.9108, 1.9151, 1.9161, 1.9181, 1.9186 c *, 1.9249, 1.9283, 1.9317, 1.9331, 1.9350 c *, 1.9384, 1.9423, 1.9519, 1.9571, 1.9614 c *, 1.9652, 1.9681, 1.9756, 1.9780, 1.9804 c *, 1.9893, 1.9921, 1.9968, 1.9987, 2.0034 c *, 2.0085, 2.0196, 2.0344, 2.0358, 2.0385 c *, 2.0545, 2.0636, 2.0654, 2.0704, 2.0816 c *, 2.0839, 2.0933, 2.1106, 2.1151, 2.1257 c *, 2.1279, 2.1388, 2.1545, 2.1606, 2.1761 c *, 2.1770, 2.1804, 2.1817, 2.1920, 2.2022 c *, 2.2060, 2.2187, 2.2480, 2.2646, 2.2934 c *, 2.3056, 2.3339, 2.3499, 2.3535, 2.3538 c *, 2.3737, 2.3894, 2.4050, 2.4128, 2.4283 c *, 2.4398, 2.4513, 2.4703, 2.4892, 2.5080 c *, 2.5266, 2.5451, 2.5561, 2.5634, 2.5816 c *, 2.5997, 2.6069, 2.6177, 2.6355, 2.6568 c *, 2.6708, 2.6987, 2.7057, 2.7195, 2.7401 c *, 2.7469, 2.7606, 2.7741, 2.7977, 2.8010 c *, 2.8077, 2.8343, 2.8409, 2.8672, 2.8769 c *, 2.8997, 2.9093, 2.9350, 2.9414, 2.9636 c *, 2.9731, 3.0013, 3.0045, 3.0107, 3.0355 c *, 3.0570, 3.0662, 3.0876, 3.0967, 3.1268 c *, 3.1328, 3.1566, 3.1862, 3.2067, 3.2155 c *, 3.2445, 3.2733, 3.3018, 3.3216, 3.3329 c *, 3.3609, 3.3887, 3.4163, 3.4190, 3.4327 c *, 3.4436, 3.4707, 3.4869, 3.4976, 3.5244 c *, 3.5509, 3.6033, 3.6550, 3.7059, 3.7462 c *, 3.9247, 3.9887, 3.9981, 4.4001, 4.4341 c *, 4.8193, 4.8387, 5.2120, 5.2246, 5.3889 c *, 5.5535, 5.5603, 5.7265, 5.8880, 5.8912 c *, 5.9671, 6.1984, 6.2421, 6.5084, 6.6369 c *, 6.9138, 7.5617, 8.1584, 8.7144, 8.9794 c *, 9.2369, 9.7314, 9.9412, 10.2020, 10.6517 c *, 11.4987, 13.7295, 15.0339, 16.6334, 16.8018 c *, 17.8835, 18.1439, 21.2400, 22.9386, 24.1342 c *, 25.2733 c */ c data pipbmx/ c * .4424, .5585, .6180, .7485, .7092 c *, .8058, .7777, .9113, .8974, .9934 c *, 1.0896, 1.2653, 1.3078, 1.5076, 1.5472 c *, 1.5656, 1.7019, 1.7504, 1.7934, 1.7626 c *, 1.8168, 1.9706, 2.0163, 2.1140, 2.1749 c *, 2.0576, 2.1851, 2.1148, 2.1924, 2.1851 c *, 2.3104, 2.3480, 2.2883, 2.1924, 2.3602 c *, 2.4570, 2.3262, 2.2708, 2.4722, 2.3296 c *, 2.5124, 2.3194, 2.4476, 2.4033, 2.5497 c *, 2.5313, 2.5514, 2.5181, 2.5125, 2.5156 c *, 2.5105, 2.4398, 2.4201, 2.5074, 2.4978 c *, 2.4463, 2.4863, 2.4856, 2.4469, 2.5231 c *, 2.3534, 2.5357, 2.4482, 2.3796, 2.3823 c *, 2.3582, 2.3650, 2.3870, 2.1705, 2.1185 c *, 2.2341, 2.2057, 2.2284, 2.1178, 2.1705 c *, 2.0622, 2.0490, 2.0122, 1.9041, 1.8873 c *, 1.9091, 1.7897, 1.7499, 1.7535, 1.8797 c *, 1.8455, 1.6087, 1.6381, 1.6361, 1.6737 c *, 1.5329, 1.5327, 1.4879, 1.4351, 1.4461 c *, 1.3854, 1.3517, 1.3334, 1.3695, 1.2741 c *, 1.2989, 1.2183, 1.2361, 1.2074, 1.1456 c *, 1.1997, 1.1439, 1.1381, 1.0786, 1.1185 c *, 1.0693, 1.1339, .9958, 1.0077, 1.0491 c *, .9997, .9611, .9452, .9141, .9190 c *, .8806, .8885, .9200, .8618, .9680 c *, .8702, .8643, .8263, .8234, .7150 c *, .8412, .8101, .8704, .7767, .8766 c *, .7722, .7694, .7956, .7590, .7128 c *, .7159, .6894, .6894, .7139, .6898 c *, .7412, .6784, .7181, .6947, .7695 c *, .6704, .7110, .6910, .6699, .7199 c *, .6935, .6759, .6843, .6933, .7190 c *, .7159, .6865, .6958, .7017, .6794 c *, .7081, .7004, .7128, .7569, .7269 c *, .7356, .7230, .7332, .7548, .7464 c *, .7695, .7703, .7707, .8288, .8117 c *, .8170, .8054, .7990, .7878, .8556 c *, .8423, .8917, .8664, .8253, .8649 c *, .8627, .8813, .8461, .8938, .9288 c *, .8759, .8345, .8972, .8625, .8818 c *, .8983, .8953, .9080, .9277, .9117 c *, .9080, .8579, .8974, .8970, .9184 c *, .9267, .9184, .9092, .8704, .9146 c *, .9184, .9591, .8649, .9296, .9071 c *, .9449, .9407, .9536, .9721, .9608 c *, .9837, .9233, .9542, 1.0004, .9853 c *, .9934, 1.0005, 1.0061, .9358, 1.0180 c *, 1.0280, 1.0573, .9982, 1.0357, 1.0555 c *, .9877, 1.0682, 1.0694, 1.0817, 1.0720 c *, 1.1046, 1.1041, 1.0596, 1.1013, 1.1217 c *, 1.1217, 1.1291, 1.1113, 1.1270, 1.1354 c *, 1.0799, 1.1199, 1.1439, 1.1459, 1.1213 c *, 1.0907, 1.1368, 1.1488, 1.1008, 1.1500 c *, 1.1156, 1.1502, 1.0947, 1.1410, 1.1427 c *, 1.1480, 1.0729, 1.1634, 1.1319, 1.1480 c *, 1.1295, 1.1213, 1.0783, 1.1044, 1.1027 c *, 1.1217, 1.0998, 1.0418, 1.1031, 1.0712 c *, 1.0720, 1.0600, 1.0540, 1.1180, 1.0490 c *, 1.0464, 1.0380, 1.0124, 1.0202, .9821 c *, 1.0013, .9945, .9997, 1.0187, .9805 c *, .9873, .9796, .9608, 1.0045, .9756 c *, .9695, .9674, .9646, .9662, .9659 c *, .9654, .9491, .9643, .9566, .9654 c *, .9636, .9619, .9676, .9774, .9747 c *, .9916, .9843, .9950, .9932, .9641 c *, .9895, .9861, .9871, .9905, .9916 c *, .9657, .9872, .9770, .9796, .9751 c *, .9637, .9689, .9633, .9635, .9601 c *, .9580, .9707, .9581, .9540, .9505 c *, .9509, .9422, .9494, .9657, .9478 c *, .9513, .9541, .9472, .9439, .9640 c *, .9452, .9397, .9440, .9390, .9416 c *, .9424, .9394, .9334, .9366, .9348 c *, .9338, .9414, .9312, .9115, .9286 c *, .9373, .9287, .9190, .9248, .9232 c *, .9339, .9214, .9201, .9183, .9181 c *, .9170, .9150, .9138, .9071, .9117 c *, .9106, .9092, .9084, .9253, .9053 c *, .9061, .9052, .9132, .9032, .9029 c *, .9008, .8973, .8954, .8928, .9115 c *, .9021, .8938, .8997, .8907, .8921 c *, .8834, .8831, .8795, .8774, .8755 c *, .8745, .8682, .8849, .8649, .8705 c *, .8616, .8684, .8691, .8634, .8687 c *, .8636, .8616, .8581, .8571, .8575 c *, .8582, .8551, .8575, .8582, .8618 c *, .8597, .8619, .8634, .8693, .8645 c *, .8682, .8538, .8759, .8818, .8831 c *, .8853/ c data (pimecm(i),i=1,400)/ c * 1.1046, 1.1133, 1.1394, 1.1425, 1.1495 c *, 1.1579, 1.1585, 1.1592, 1.1598, 1.1677 c *, 1.1691, 1.1731, 1.1777, 1.1784, 1.1798 c *, 1.1831, 1.1858, 1.1879, 1.1892, 1.1906 c *, 1.1940, 1.1967, 1.1994, 1.2008, 1.2015 c *, 1.2028, 1.2069, 1.2076, 1.2090, 1.2097 c *, 1.2111, 1.2124, 1.2131, 1.2159, 1.2166 c *, 1.2179, 1.2186, 1.2200, 1.2214, 1.2234 c *, 1.2269, 1.2283, 1.2296, 1.2303, 1.2317 c *, 1.2324, 1.2352, 1.2358, 1.2365, 1.2372 c *, 1.2407, 1.2421, 1.2441, 1.2462, 1.2476 c *, 1.2510, 1.2517, 1.2524, 1.2545, 1.2551 c *, 1.2579, 1.2586, 1.2593, 1.2607, 1.2613 c *, 1.2627, 1.2634, 1.2641, 1.2662, 1.2668 c *, 1.2682, 1.2696, 1.2710, 1.2716, 1.2730 c *, 1.2758, 1.2778, 1.2806, 1.2819, 1.2826 c *, 1.2833, 1.2847, 1.2888, 1.2915, 1.2922 c *, 1.2929, 1.2963, 1.3004, 1.3024, 1.3038 c *, 1.3058, 1.3065, 1.3072, 1.3112, 1.3119 c *, 1.3126, 1.3146, 1.3180, 1.3187, 1.3200 c *, 1.3207, 1.3220, 1.3227, 1.3254, 1.3261 c *, 1.3301, 1.3308, 1.3321, 1.3335, 1.3362 c *, 1.3375, 1.3388, 1.3415, 1.3422, 1.3455 c *, 1.3482, 1.3515, 1.3522, 1.3555, 1.3562 c *, 1.3615, 1.3628, 1.3641, 1.3655, 1.3681 c *, 1.3694, 1.3760, 1.3773, 1.3786, 1.3898 c *, 1.3917, 1.3963, 1.3995, 1.4002, 1.4015 c *, 1.4047, 1.4099, 1.4112, 1.4163, 1.4176 c *, 1.4183, 1.4196, 1.4208, 1.4228, 1.4234 c *, 1.4272, 1.4292, 1.4298, 1.4304, 1.4393 c *, 1.4425, 1.4470, 1.4476, 1.4489, 1.4507 c *, 1.4539, 1.4589, 1.4596, 1.4608, 1.4652 c *, 1.4658, 1.4683, 1.4727, 1.4739, 1.4758 c *, 1.4764, 1.4777, 1.4795, 1.4808, 1.4833 c *, 1.4851, 1.4857, 1.4876, 1.4882, 1.4894 c *, 1.4901, 1.4919, 1.4925, 1.4938, 1.4968 c *, 1.4974, 1.4981, 1.4993, 1.5030, 1.5054 c *, 1.5060, 1.5066, 1.5072, 1.5097, 1.5103 c *, 1.5109, 1.5146, 1.5152, 1.5182, 1.5206 c *, 1.5212, 1.5224, 1.5231, 1.5249, 1.5309 c *, 1.5327, 1.5345, 1.5357, 1.5387, 1.5405 c *, 1.5411, 1.5429, 1.5435, 1.5441, 1.5465 c *, 1.5489, 1.5507, 1.5513, 1.5519, 1.5548 c *, 1.5590, 1.5614, 1.5637, 1.5655, 1.5673 c *, 1.5702, 1.5708, 1.5761, 1.5802, 1.5825 c *, 1.5831, 1.5843, 1.5849, 1.5866, 1.5890 c *, 1.5936, 1.5948, 1.5977, 1.5994, 1.6017 c *, 1.6052, 1.6058, 1.6087, 1.6092, 1.6138 c *, 1.6173, 1.6178, 1.6190, 1.6236, 1.6258 c *, 1.6276, 1.6287, 1.6298, 1.6315, 1.6327 c *, 1.6367, 1.6372, 1.6378, 1.6400, 1.6406 c *, 1.6423, 1.6480, 1.6513, 1.6519, 1.6525 c *, 1.6547, 1.6553, 1.6614, 1.6631, 1.6637 c *, 1.6648, 1.6653, 1.6659, 1.6681, 1.6692 c *, 1.6703, 1.6720, 1.6726, 1.6731, 1.6742 c *, 1.6792, 1.6798, 1.6803, 1.6820, 1.6825 c *, 1.6853, 1.6858, 1.6864, 1.6875, 1.6935 c *, 1.6957, 1.6963, 1.6979, 1.6990, 1.7067 c *, 1.7072, 1.7121, 1.7127, 1.7132, 1.7186 c *, 1.7202, 1.7208, 1.7213, 1.7224, 1.7235 c *, 1.7240, 1.7267, 1.7278, 1.7294, 1.7342 c *, 1.7348, 1.7391, 1.7407, 1.7428, 1.7476 c *, 1.7503, 1.7535, 1.7545, 1.7609, 1.7614 c *, 1.7667, 1.7714, 1.7735, 1.7741, 1.7767 c *, 1.7772, 1.7798, 1.7824, 1.7835, 1.7871 c *, 1.7929, 1.7960, 1.8001, 1.8032, 1.8058 c *, 1.8130, 1.8182, 1.8212, 1.8218, 1.8238 c *, 1.8253, 1.8258, 1.8289, 1.8320, 1.8370 c *, 1.8386, 1.8436, 1.8446, 1.8507, 1.8512 c *, 1.8557, 1.8588, 1.8638, 1.8678, 1.8713 c *, 1.8762, 1.8777, 1.8792, 1.8807, 1.8827 c *, 1.8857, 1.8881, 1.8886, 1.8936, 1.8975 c *, 1.9010, 1.9063, 1.9132, 1.9186, 1.9205 c *, 1.9249, 1.9254, 1.9273, 1.9283, 1.9317 c *, 1.9326, 1.9374, 1.9423, 1.9451, 1.9495 c *, 1.9519, 1.9538, 1.9614, 1.9709, 1.9733 c *, 1.9747, 1.9799, 1.9884, 1.9903, 1.9987 c *, 2.0034, 2.0117, 2.0164, 2.0182, 2.0196 c *, 2.0330, 2.0335, 2.0358, 2.0540, 2.0581 c *, 2.0636, 2.0654, 2.0798, 2.0812, 2.0816 c */ c data (pimecm(i),i=401,578)/ c * 2.0933, 2.1013, 2.1075, 2.1221, 2.1261 c *, 2.1305, 2.1375, 2.1445, 2.1588, 2.1658 c *, 2.1740, 2.1804, 2.1877, 2.2026, 2.2119 c *, 2.2166, 2.2187, 2.2305, 2.2510, 2.2626 c *, 2.2688, 2.2717, 2.3040, 2.3121, 2.3315 c *, 2.3483, 2.3538, 2.3737, 2.3744, 2.3925 c *, 2.4050, 2.4105, 2.4128, 2.4267, 2.4302 c *, 2.4513, 2.4627, 2.4892, 2.5069, 2.5266 c *, 2.5561, 2.5634, 2.5802, 2.5925, 2.5997 c *, 2.6284, 2.6355, 2.6557, 2.6708, 2.6987 c *, 2.7057, 2.7401, 2.7673, 2.7741, 2.7966 c *, 2.8077, 2.8343, 2.8409, 2.8672, 2.8769 c *, 2.8997, 2.9093, 2.9318, 2.9341, 2.9414 c *, 2.9636, 2.9731, 3.0045, 3.0262, 3.0355 c *, 3.0570, 3.0656, 3.0662, 3.0876, 3.0967 c *, 3.1268, 3.1566, 3.1774, 3.1862, 3.2067 c *, 3.2155, 3.2445, 3.2733, 3.3018, 3.3216 c *, 3.3329, 3.3609, 3.3887, 3.4163, 3.4190 c *, 3.4436, 3.4675, 3.4707, 3.4869, 3.4976 c *, 3.5509, 3.6033, 3.6550, 3.6575, 3.7059 c *, 3.7312, 3.7462, 3.8402, 3.8935, 3.9887 c *, 4.0004, 4.3873, 4.4341, 4.6811, 4.7210 c *, 4.7408, 4.8387, 4.8406, 5.0288, 5.0844 c *, 5.2120, 5.2353, 5.3889, 5.4254, 5.5603 c *, 5.6123, 5.7265, 5.7787, 5.8880, 5.9451 c *, 5.9671, 5.9953, 6.0792, 6.1984, 6.2241 c *, 6.3479, 6.5070, 6.6369, 6.8140, 6.9138 c *, 7.0734, 7.2450, 7.3962, 7.5617, 7.7092 c *, 7.9841, 7.9865, 8.1584, 8.1814, 8.3628 c *, 8.4410, 8.6582, 8.6593, 8.7144, 8.9794 c *, 9.2369, 9.2845, 9.4874, 9.7314, 9.7775 c *, 9.9694, 10.2020, 10.3580, 10.4293, 10.5988 c *, 10.6517, 10.8697, 11.0833, 11.4987, 13.7295 c *, 15.0339, 16.6334, 16.8018, 17.8835, 18.1439 c *, 19.3933, 19.6336, 21.2400, 22.9386, 24.1342 c *, 25.2733, 26.0050, 26.3632 c */ c data (pimbmx(i),i=1,400)/ c * .5314, .5314, .6580, .7092, .7506 c *, .8215, .8579, .8406, .8349, .9441 c *, .9657, 1.0363, 1.0686, 1.1085, 1.1185 c *, 1.1899, 1.1658, 1.2218, 1.2322, 1.2653 c *, 1.3273, 1.2374, 1.3611, 1.3716, 1.3267 c *, 1.3986, 1.4150, 1.3399, 1.4701, 1.4516 c *, 1.4723, 1.4955, 1.4494, 1.4161, 1.4558 c *, 1.4625, 1.5006, 1.5128, 1.5160, 1.4584 c *, 1.4991, 1.4791, 1.4217, 1.4799, 1.4799 c *, 1.4588, 1.4844, 1.4172, 1.4273, 1.4439 c *, 1.4166, 1.4127, 1.3739, 1.3624, 1.3971 c *, 1.3297, 1.3231, 1.3285, 1.2878, 1.2952 c *, 1.2679, 1.2641, 1.2976, 1.2386, 1.2463 c *, 1.3025, 1.2489, 1.2239, 1.1902, 1.2114 c *, 1.1955, 1.2083, 1.1658, 1.1686, 1.1603 c *, 1.1424, 1.1185, 1.1160, 1.1013, 1.1070 c *, 1.1041, 1.0823, 1.0645, 1.0779, 1.0295 c *, 1.0311, .9950, 1.0155, .9950, 1.0029 c *, .9767, 1.0249, .9853, .9657, .9906 c *, .9608, .9591, .9805, .9458, .9575 c *, .9422, .9407, .9558, .9473, .9271 c *, .9132, .9236, .9219, .9219, .9575 c *, .9097, .9224, .9094, .9202, .8903 c *, .9224, .9160, .9122, .9094, .8956 c *, .9167, .9271, .9398, .9184, .9277 c *, .9310, .9184, .9354, .9219, .9300 c *, .9303, .9508, .9451, .9680, .9399 c *, .9420, .9387, .9667, .9248, .9915 c *, .9698, .9626, 1.0187, .9659, .9441 c *, .9766, .9694, 1.0061, .9789, .9757 c *, .9803, 1.0495, 1.0218, 1.0026, 1.0706 c *, 1.0277, 1.0302, 1.0690, 1.0552, 1.0406 c *, 1.0801, 1.0522, 1.1284, 1.1226, 1.1270 c *, 1.0934, 1.0968, 1.1381, 1.0915, 1.1543 c *, 1.1775, 1.2101, 1.1247, 1.1506, 1.1781 c *, 1.1546, 1.1968, 1.1604, 1.1697, 1.2072 c *, 1.2058, 1.2489, 1.1912, 1.2141, 1.1928 c *, 1.1986, 1.1978, 1.2256, 1.2565, 1.2029 c *, 1.2127, 1.1997, 1.2035, 1.1791, 1.1853 c *, 1.1914, 1.1927, 1.2616, 1.2303, 1.1848 c *, 1.1507, 1.2399, 1.1510, 1.1898, 1.1499 c *, 1.1337, 1.1433, 1.1089, 1.1517, 1.2008 c *, 1.1598, 1.1030, 1.0798, 1.1165, 1.1213 c *, 1.1190, 1.0861, 1.0776, 1.0653, 1.0730 c *, 1.1106, 1.0908, 1.0945, 1.0650, 1.0882 c *, 1.0720, 1.0540, 1.0925, 1.0718, 1.1377 c *, 1.0867, 1.0802, 1.1090, 1.0936, 1.0915 c *, 1.1316, 1.1575, 1.1202, 1.1275, 1.0908 c *, 1.1142, 1.1507, 1.1890, 1.1863, 1.2114 c *, 1.1142, 1.2001, 1.1312, 1.1930, 1.2563 c *, 1.2489, 1.2919, 1.2652, 1.2425, 1.2149 c *, 1.2386, 1.2616, 1.3207, 1.3299, 1.2816 c *, 1.3296, 1.3562, 1.2957, 1.3831, 1.3028 c *, 1.3806, 1.3704, 1.3202, 1.4082, 1.3482 c *, 1.3730, 1.3351, 1.3886, 1.3348, 1.3957 c *, 1.3929, 1.3790, 1.3791, 1.3946, 1.3641 c *, 1.3877, 1.3566, 1.3564, 1.3762, 1.3253 c *, 1.3739, 1.3252, 1.3234, 1.3423, 1.3159 c *, 1.2643, 1.2700, 1.2584, 1.2466, 1.2841 c *, 1.2191, 1.2052, 1.2794, 1.2835, 1.2374 c *, 1.2061, 1.1936, 1.1947, 1.2476, 1.1481 c *, 1.1496, 1.2187, 1.1427, 1.2101, 1.1192 c *, 1.1445, 1.1371, 1.1613, 1.0994, 1.1020 c *, 1.1233, 1.1041, 1.0896, 1.0890, 1.0880 c *, 1.0802, 1.0946, 1.0940, 1.0670, 1.0761 c *, 1.0817, 1.0702, 1.0723, 1.0779, 1.0789 c *, 1.0733, 1.0771, 1.0767, 1.0710, 1.0705 c *, 1.0720, 1.0729, 1.0472, 1.0959, 1.0788 c *, 1.0739, 1.0645, 1.0681, 1.0749, 1.0789 c *, 1.0799, 1.0795, 1.0763, 1.0588, 1.0805 c *, 1.0724, 1.0660, 1.0663, 1.0801, 1.0988 c *, 1.0797, 1.0816, 1.0687, 1.0690, 1.0724 c *, 1.0721, 1.0777, 1.0650, 1.0715, 1.0617 c *, 1.0690, 1.0678, 1.0714, 1.0323, 1.0672 c *, 1.0630, 1.0564, 1.0600, 1.0636, 1.0499 c *, 1.0479, 1.0481, 1.0463, 1.0510, 1.0450 c *, 1.0171, 1.0517, 1.0498, 1.0510, 1.0412 c *, 1.0374, 1.0501, 1.0388, 1.0449, 1.0449 c *, 1.0311, 1.0510, 1.0505, 1.0449, 1.0567 c *, 1.0572, 1.0495, 1.0632, 1.0461, 1.0403 c */ c data (pimbmx(i),i=401,578)/ c * 1.0629, 1.0642, 1.0525, 1.0709, 1.0585 c *, 1.0612, 1.0717, 1.0731, 1.0660, 1.0761 c *, 1.0696, 1.0767, 1.0755, 1.0720, 1.0615 c *, 1.0665, 1.0660, 1.0714, 1.0670, 1.0650 c *, 1.0627, 1.0621, 1.0499, 1.0499, 1.0457 c *, 1.0405, 1.0372, 1.0412, 1.0331, 1.0299 c *, 1.0299, 1.0241, 1.0321, 1.0306, 1.0232 c *, 1.0252, 1.0171, 1.0211, 1.0226, 1.0178 c *, 1.0083, 1.0157, 1.0138, .9918, 1.0140 c *, 1.0028, 1.0120, 1.0077, 1.0088, .9979 c *, 1.0052, 1.0025, .9876, .9987, .9990 c *, .9944, .9845, .9918, .9777, .9892 c *, .9754, .9856, .9812, .9816, .9831 c *, .9674, .9800, .9781, .9679, .9756 c *, .9805, .9759, .9739, .9672, .9718 c *, .9701, .9682, .9632, .9664, .9538 c *, .9647, .9628, .9613, .9596, .9449 c *, .9585, .9572, .9553, .9540, .9624 c *, .9527, .9468, .9513, .9525, .9503 c *, .9479, .9459, .9438, .9424, .9422 c *, .9407, .9457, .9399, .9385, .9368 c *, .9409, .9248, .9209, .9034, .8974 c *, .9150, .9089, .9145, .9146, .9080 c *, .9044, .9082, .9069, .9062, .8938 c *, .9034, .9045, .9011, .8921, .8979 c *, .8925, .8982, .8975, .8976, .8947 c *, .8962, .8932, .8916, .8913, .8888 c *, .8889, .8892, .8880, .8881, .8842 c *, .8815, .8826, .8832, .8813, .8826 c *, .8791, .8835, .8831, .8833, .8782 c *, .8789, .8813, .8797, .8780, .8744 c *, .8782, .8842, .8777, .8777, .8789 c *, .8815, .8800, .8839, .8740, .8732 c *, .8751, .8784, .8757, .8779, .8630 c *, .8802, .8775, .8851, .8874, .8903 c *, .8935, .8965, .8965 c */ c data kmpecm/ c * 1.4691, 1.4720, 1.4750, 1.4780, 1.4811 c *, 1.4837, 1.4843, 1.4860, 1.4876, 1.4910 c *, 1.4944, 1.4979, 1.5014, 1.5032, 1.5050 c *, 1.5087, 1.5091, 1.5124, 1.5132, 1.5162 c *, 1.5170, 1.5200, 1.5220, 1.5239, 1.5278 c *, 1.5318, 1.5354, 1.5358, 1.5362, 1.5378 c *, 1.5399, 1.5440, 1.5523, 1.5607, 1.5654 c *, 1.5688, 1.5775, 1.5784, 1.5863, 1.5916 c *, 1.5947, 1.6023, 1.6050, 1.6055, 1.6086 c *, 1.6145, 1.6159, 1.6172, 1.6191, 1.6236 c *, 1.6319, 1.6328, 1.6332, 1.6420, 1.6461 c *, 1.6466, 1.6522, 1.6563, 1.6582, 1.6614 c *, 1.6642, 1.6694, 1.6712, 1.6717, 1.6768 c *, 1.6806, 1.6811, 1.6839, 1.6843, 1.6867 c *, 1.6885, 1.6961, 1.6965, 1.7003, 1.7022 c *, 1.7083, 1.7087, 1.7172, 1.7177, 1.7181 c *, 1.7229, 1.7243, 1.7276, 1.7342, 1.7374 c *, 1.7436, 1.7459, 1.7483, 1.7539, 1.7610 c *, 1.7629, 1.7633, 1.7718, 1.7770, 1.7789 c *, 1.7793, 1.7817, 1.7840, 1.7873, 1.7892 c *, 1.8028, 1.8037, 1.8075, 1.8135, 1.8140 c *, 1.8219, 1.8247, 1.8261, 1.8308, 1.8369 c *, 1.8401, 1.8406, 1.8410, 1.8475, 1.8480 c *, 1.8489, 1.8540, 1.8559, 1.8605, 1.8647 c *, 1.8721, 1.8744, 1.8767, 1.8785, 1.8836 c *, 1.8951, 1.8955, 1.8983, 1.9001, 1.9029 c *, 1.9065, 1.9184, 1.9202, 1.9243, 1.9348 c *, 1.9393, 1.9411, 1.9434, 1.9483, 1.9547 c *, 1.9628, 1.9637, 1.9659, 1.9699, 1.9798 c *, 1.9838, 1.9923, 1.9958, 2.0047, 2.0127 c *, 2.0149, 2.0162, 2.0255, 2.0277, 2.0308 c *, 2.0417, 2.0430, 2.0487, 2.0579, 2.0653 c *, 2.0679, 2.0713, 2.0813, 2.0882, 2.0925 c *, 2.1028, 2.1105, 2.1126, 2.1211, 2.1232 c *, 2.1258, 2.1351, 2.1444, 2.1507, 2.1528 c *, 2.1654, 2.1675, 2.1687, 2.1737, 2.1837 c *, 2.1862, 2.1937, 2.2044, 2.2131, 2.2143 c *, 2.2274, 2.2356, 2.2478, 2.2546, 2.2659 c *, 2.2756, 2.2836, 2.2976, 2.2996, 2.3162 c *, 2.3166, 2.3296, 2.3335, 2.3363, 2.3535 c *, 2.3562, 2.3725, 2.3729, 2.3748, 2.3887 c *, 2.3918, 2.3933, 2.4006, 2.4109, 2.4174 c *, 2.4223, 2.4299, 2.4352, 2.4488, 2.4518 c *, 2.4675, 2.4705, 2.4861, 2.4887, 2.4936 c *, 2.5046, 2.5230, 2.5412, 2.5593, 2.5701 c *, 2.5773, 2.5952, 2.6023, 2.6059, 2.6130 c *, 2.6306, 2.6447, 2.6482, 2.6656, 2.6795 c *, 2.6829, 2.7002, 2.7173, 2.7848, 2.8540 c *, 2.9216, 2.9407, 2.9878, 3.0096, 3.0250 c *, 3.0526, 3.1783, 3.2191, 3.2993, 3.3887 c *, 3.5239, 3.6924, 3.7800, 3.9967, 4.0200 c *, 4.1348, 4.4617, 4.4826, 4.7663, 4.8636 c *, 4.9779, 5.1080, 5.1263, 5.2349, 5.3237 c *, 5.4110, 5.5479, 5.5816, 5.8281, 5.9080 c *, 6.0801, 6.2173, 6.3664, 6.6545, 6.9306 c *, 7.2610, 7.5770, 7.7605, 7.8207, 7.9985 c *, 8.1725, 8.2297, 8.4546, 8.7275, 8.9922 c *, 9.2493, 9.4994, 9.7431, 9.9809, 10.2131 c *, 11.5086, 13.7378, 15.0415, 16.8085, 17.8898 c *, 19.3991, 21.2453, 22.9435, 24.1389 c */ c data kmpbmx/ c * 1.9033, 1.7662, 1.7298, 1.7544, 1.5461 c *, 1.6991, 1.6205, 1.5898, 1.5817, 1.5023 c *, 1.5554, 1.5086, 1.5065, 1.4948, 1.4799 c *, 1.4927, 1.4854, 1.6136, 1.6017, 1.7312 c *, 1.4884, 1.7075, 1.5574, 1.5260, 1.5002 c *, 1.4571, 1.3991, 1.3219, 1.3434, 1.3635 c *, 1.3315, 1.2966, 1.2213, 1.1604, 1.1930 c *, 1.1382, 1.1354, 1.1583, 1.1185, 1.1156 c *, 1.0645, 1.0247, 1.0969, 1.0779, 1.1083 c *, 1.0749, 1.0155, 1.0464, 1.0754, 1.0202 c *, 1.0470, 1.0615, 1.0540, 1.0357, 1.0605 c *, 1.0428, 1.0265, 1.0187, 1.0611, 1.0615 c *, 1.0295, 1.0651, 1.0823, 1.0077, 1.0764 c *, 1.1298, 1.1013, 1.1080, 1.1070, 1.1368 c *, 1.1382, 1.1410, 1.1418, 1.1312, 1.1466 c *, 1.1378, 1.1241, 1.0896, 1.1264, 1.0720 c *, 1.1368, 1.1013, 1.1442, 1.1466, 1.1612 c *, 1.1726, 1.1755, 1.1576, 1.1726, 1.1848 c *, 1.1617, 1.2029, 1.1686, 1.2274, 1.2252 c *, 1.1726, 1.2256, 1.2244, 1.2008, 1.2332 c *, 1.2828, 1.2548, 1.2412, 1.2887, 1.2853 c *, 1.2653, 1.2540, 1.2527, 1.2449, 1.2118 c *, 1.1781, 1.1902, 1.1767, 1.1859, 1.1410 c *, 1.1594, 1.1800, 1.1185, 1.1256, 1.1354 c *, 1.1131, 1.1562, 1.1143, 1.1382, 1.0841 c *, 1.0587, 1.0632, 1.0311, 1.0617, 1.0150 c *, 1.0309, 1.0174, 1.0110, 1.0171, .9961 c *, .9719, .9938, .9869, .9953, .9966 c *, 1.0189, .9977, .9918, .9948, 1.0034 c *, .9737, 1.0066, 1.0137, 1.0041, 1.0171 c *, 1.0303, 1.0223, 1.0322, .9751, 1.0331 c *, 1.0118, 1.0403, 1.0399, 1.0429, 1.0171 c *, 1.0279, 1.0467, 1.0414, 1.0217, 1.0432 c *, 1.0379, 1.0280, 1.0351, 1.0171, 1.0278 c *, 1.0314, 1.0213, 1.0133, 1.0240, 1.0080 c *, 1.0018, .9964, .9984, .9889, .9910 c *, .9903, .9801, .9852, .9853, .9847 c *, .9800, .9832, .9770, .9749, .9754 c *, .9723, .9741, .9744, .9738, .9751 c *, .9780, .9738, .9738, .9684, .9712 c *, .9669, .9680, .9671, .9576, .9619 c *, .9624, .9585, .9588, .9586, .9525 c *, .9253, .9518, .9503, .9491, .9463 c *, .9476, .9420, .9458, .9434, .9341 c *, .9444, .9412, .9393, .9395, .9210 c *, .9370, .9358, .8974, .9400, .9342 c *, .9305, .9141, .9271, .9267, .9228 c *, .9233, .9219, .9247, .9296, .9253 c *, .9089, .8992, .8946, .8992, .9062 c *, .9069, .8874, .9219, .8746, .8795 c *, .8740, .8704, .8785, .8795, .8667 c *, .8849, .8510, .8463, .8612, .8292 c *, .8292, .8387, .8273, .8273, .8292 c *, .8292, .8349, .8234, .8349, .8176 c *, .8292, .8180, .8193, .8154, .8130 c *, .8121, .8145, .8078, .7959, .8088 c *, .8075, .8064, .8056, .8086, .8048 c *, .8080, .8068, .8050, .8042, .8050 c *, .8054, .8064, .8096, .8095, .8107 c *, .8135, .8234, .8238, .8263 c */ c data kmnecm/ c * 1.6212, 1.6781, 1.7053, 1.7863, 1.7882 c *, 1.8271, 1.9025, 2.0152, 2.1237, 2.2157 c *, 2.4252, 2.6054, 2.6160, 2.9441, 3.5279 c *, 3.6965, 4.0245, 4.4666, 4.8690, 5.1136 c *, 5.2406, 5.4169, 5.5877, 5.9144, 6.2241 c *, 6.9381, 7.5852, 8.1813, 8.7369, 9.2592 c *, 9.7536, 10.2241, 11.5209, 13.7525, 15.0575 c *, 16.8264, 17.9089, 19.4198, 21.2680, 22.9680 c *, 24.1646 c */ c data kmnbmx/ c * .9132, .9966, .9772, .9966, 1.1382 c *, 1.0794, .9674, .9167, .8795, .8500 c *, .8482, .8444, .8330, .8078, .8349 c *, .8368, .7919, .8146, .8019, .8156 c *, .7999, .7961, .8038, .8038, .7949 c *, .7868, .7864, .7910, .7850, .7870 c *, .7906, .7885, .7945, .7971, .8011 c *, .8003, .8029, .8082, .8088, .8156 c *, .8189 c */ c data kppecm/ c * 1.4447, 1.4516, 1.4522, 1.4585, 1.4663 c *, 1.4750, 1.5036, 1.5050, 1.5091, 1.5239 c *, 1.5378, 1.5523, 1.5654, 1.5784, 1.5916 c *, 1.5929, 1.6014, 1.6037, 1.6050, 1.6150 c *, 1.6159, 1.6191, 1.6259, 1.6264, 1.6268 c *, 1.6328, 1.6378, 1.6461, 1.6517, 1.6587 c *, 1.6605, 1.6652, 1.6792, 1.6843, 1.6853 c *, 1.6928, 1.7073, 1.7101, 1.7210, 1.7294 c *, 1.7374, 1.7422, 1.7483, 1.7539, 1.7643 c *, 1.7662, 1.7704, 1.7789, 1.7793, 1.7864 c *, 1.7897, 1.8028, 1.8070, 1.8135, 1.8191 c *, 1.8327, 1.8355, 1.8373, 1.8517, 1.8587 c *, 1.8605, 1.8679, 1.8725, 1.8813, 1.8836 c *, 1.9038, 1.9070, 1.9289, 1.9298, 1.9321 c *, 1.9524, 1.9533, 1.9749, 1.9807, 1.9950 c *, 1.9972, 1.9994, 2.0074, 2.0193, 2.0435 c *, 2.0492, 2.0635, 2.0653, 2.0851, 2.1040 c *, 2.1066, 2.1083, 2.1279, 2.1296, 2.1490 c *, 2.1507, 2.1717, 2.1908, 2.1924, 2.2110 c *, 2.2131, 2.2213, 2.2319, 2.2335, 2.2538 c *, 2.2724, 2.2740, 2.2940, 2.3123, 2.3138 c *, 2.3375, 2.3531, 2.3725, 2.3903, 2.3918 c *, 2.4109, 2.4197, 2.4299, 2.4413, 2.4488 c *, 2.4675, 2.4861, 2.5046, 2.5230, 2.5266 c *, 2.5412, 2.5521, 2.5593, 2.5773, 2.5952 c *, 2.6130, 2.6306, 2.6482, 2.6656, 2.6829 c *, 2.7002, 2.7173, 3.0096, 3.0556, 3.1754 c *, 3.3887, 3.5239, 3.7800, 4.0200, 4.1348 c *, 4.4617, 4.6469, 4.7663, 4.8636, 4.9591 c *, 5.1263, 5.2349, 5.4110, 5.5816, 5.7308 c *, 5.9080, 6.0646, 6.2173, 6.9306, 7.5770 c *, 7.8207, 8.1725, 8.7275, 9.2493, 9.7431 c *, 10.2131, 11.5086, 13.7378, 15.0415, 16.6402 c *, 16.8085, 17.8898, 18.1501, 19.3991, 21.2453 c *, 22.9435, 24.1389 c */ c data kppbmx/ c * .5412, .6308, .6024, .6050, .5971 c *, .6037, .6232, .6103, .6482, .6601 c *, .6386, .6466, .6438, .6204, .6482 c *, .6358, .6333, .6445, .6443, .6346 c *, .6410, .6227, .6283, .6308, .6403 c *, .6290, .6457, .5984, .6333, .5955 c *, .5944, .6295, .6346, .6090, .6433 c *, .6383, .6482, .6425, .6543, .6588 c *, .6652, .6768, .6730, .6723, .6815 c *, .7040, .6898, .7014, .7001, .7181 c *, .7130, .7159, .7067, .7440, .7345 c *, .7365, .7485, .7382, .7474, .7574 c *, .7588, .7559, .7590, .7582, .7668 c *, .7592, .7682, .7661, .7697, .7548 c *, .7661, .7626, .7626, .7563, .7590 c *, .7578, .7611, .7557, .7555, .7506 c *, .7498, .7517, .7508, .7540, .7464 c *, .7538, .7512, .7527, .7534, .7527 c *, .7565, .7521, .7529, .7525, .7444 c *, .7517, .7334, .7485, .7491, .7510 c *, .7466, .7476, .7478, .7472, .7485 c *, .7378, .7451, .7468, .7474, .7476 c *, .7459, .7410, .7461, .7457, .7414 c *, .7464, .7457, .7444, .7444, .7444 c *, .7442, .7291, .7421, .7429, .7421 c *, .7397, .7386, .7373, .7389, .7384 c *, .7384, .7386, .7378, .7569, .7878 c *, .7548, .7356, .7526, .7421, .7715 c *, .7568, .7590, .7777, .7421, .7632 c *, .7464, .7442, .7548, .7367, .7736 c *, .7378, .7421, .7455, .7502, .7510 c *, .7653, .7529, .7580, .7544, .7614 c *, .7605, .7661, .7703, .7805, .7883 c *, .7850, .7907, .7611, .7961, .8023 c *, .8068, .8111 c */ c data kpnecm/ c * 1.6875, 1.7430, 1.7670, 1.7816, 1.7905 c *, 1.8144, 1.8383, 1.8615, 1.8749, 1.8846 c *, 1.9080, 1.9308, 1.9345, 1.9535, 1.9760 c *, 1.9974, 1.9983, 2.0205, 2.0460, 2.0647 c *, 2.0678, 2.0864, 2.1066, 2.1079, 2.1109 c *, 2.1292, 2.1322, 2.1504, 2.1533, 2.1743 c *, 2.1922, 2.1951, 2.2157, 2.2239, 2.2334 c *, 2.2362, 2.2565, 2.2739, 2.2767, 2.2967 c *, 2.3138, 2.3166, 2.3402, 2.3559, 2.3753 c *, 2.3919, 2.3946, 2.4138, 2.4328, 2.4517 c *, 2.4704, 2.4891, 2.5076, 2.5259, 2.5442 c *, 2.5551, 2.5623, 2.5803, 2.5982, 2.6160 c *, 2.6337, 2.6513, 2.6687, 2.6861, 2.7033 c *, 2.7204, 3.5279, 4.0245, 4.4666, 4.8690 c *, 5.2406, 5.4169, 5.5877, 5.9144, 6.2241 c *, 6.9381, 7.5852, 8.1813, 8.7369, 9.2592 c *, 9.7536, 10.2241, 11.5209, 13.7525, 15.0575 c *, 16.8264, 17.9089, 19.4198, 21.2680, 22.9680 c *, 24.1646 c */ c data kpnbmx/ c * .7024, .7324, .7485, .7527, .7680 c *, .7758, .8100, .8224, .7611, .8151 c *, .8031, .7915, .7674, .7842, .7822 c *, .7590, .7791, .7767, .7758, .7734 c *, .7754, .7709, .7674, .7713, .7742 c *, .7752, .7748, .7721, .7680, .7707 c *, .7674, .7713, .7715, .7695, .7684 c *, .7734, .7682, .7709, .7672, .7659 c *, .7653, .7653, .7506, .7626, .7624 c *, .7701, .7588, .7622, .7592, .7491 c *, .7588, .7574, .7592, .7582, .7571 c *, .7464, .7559, .7538, .7529, .7529 c *, .7534, .7538, .7487, .7487, .7498 c *, .7474, .7464, .7485, .7464, .7485 c *, .7464, .7565, .7442, .7485, .7545 c *, .7555, .7542, .7634, .7653, .7686 c *, .7671, .7723, .7695, .7785, .7824 c *, .7905, .7927, .7923, .8052, .8100 c *, .8137 c */ cc data lpecm/ cc * 2.0577, 2.0583, 2.0593, 2.0595, 2.0609 cc *, 2.0617, 2.0629, 2.0642, 2.0643, 2.0647 cc *, 2.0666, 2.0671, 2.0683, 2.0692, 2.0709 cc *, 2.0720, 2.0783, 2.0818, 2.0935, 2.1022 cc *, 2.1058, 2.1117, 2.1326, 2.1559, 2.1731 cc *, 2.1811, 2.2079, 2.2360, 2.2505, 2.3107 cc *, 2.3733, 2.4534, 2.4695, 2.9278, 5.2476 cc */ cc data lpbmx/ cc * 2.5793, 2.3937, 1.9381, 2.3736, 2.0342 cc *, 2.2068, 1.9381, 1.8797, 1.7841, 1.7930 cc *, 1.2676, 1.6641, 1.4604, 1.0403, 1.3470 cc *, 1.0420, .8722, .8368, .5323, .5352 cc *, .8867, .5382, .8106, .8043, .8058 cc *, .8630, .7031, .7777, .9441, 1.0403 cc *, .9772, .6628, 1.1968, 1.0555, 1.0495 cc */ c data pi1ecm/ c * .45000, .55000, .62500, .68750, .73750 c *, .78750, .83750, .88750, .95000, 1.02500 c *, 1.10000, 1.17500 c */ c data pi1bmx/ c * .48533, .47873, .55852, .58087, .55279 c *, .49185, .47203, .46181, .46181, .47203 c *, .41459, .40684 c */ c data pi2ecm/ c * .45000, .55000, .62500, .68750, .73750 c *, .78750, .83750, .88750, .95000, 1.02500 c *, 1.10000, 1.17500 c */ c data pi2bmx/ c * .45135, .80385, .87767, 1.28779, 1.81771 c *, 2.07296, 1.55536, .95912, .94407, .77768 c *, .70467, .54408 c */ c data pi3ecm/ c * .45000, .55000, .61250, .65000, .67500 c *, .70000, .72500, .75000, .77500, .80000 c *, .82500, .85000, .87500, .90000, .93750 c *, 1.02500, 1.10000, 1.17500 c */ c data pi3bmx/ c * .50463, .79988, 1.13541, 1.30497, 1.34462 c *, 1.52957, 1.63809, 1.94216, 1.98511, 1.66316 c *, 1.69163, 1.62933, 1.06152, .96574, .83302 c *, .71365, .66517, .59173 c */ c data pi4ecm/ c * .28500, .29500, .30500, .31500, .37000 c *, .46000, .54375, .61250, .65000, .67500 c *, .70000, .72500, .75000, .77500, .80000 c *, .82500, .85000, .87500, .90125, .92750 c *, .96000 c */ c data pi4bmx/ c * .29854, .46181, .54115, .54115, .51709 c *, .72031, .80976, .96078, .90798, 1.07788 c *, 1.37389, 1.70101, 2.02716, 1.83687, 1.52644 c *, 1.57469, 1.35640, .93390, .98531, .84062 c *, .79188 c */ c data pi5ecm/0.75,0.80,0.85,0.90,0.95,1.00,1.05,1.10,1.15/ c data pi5bmx/0.4,0.5,0.7,0.9,1.1,1.1,0.7,0.6,0.5/ c c call utpri('jintcs',ish,ishini,7) c c ics=0 c c if(idptl(i).gt.10000)then c call idquad(i,nqi,nai,jci) c if(nqi.eq.0)then c idi=999 c else c idi=9999 c endif c else c idi=idptl(i) c endif c c if(idptl(j).gt.10000)then c call idquad(j,nqj,naj,jcj) c if(nqj.eq.0)then c idj=999 c else c idj=9999 c endif c else c idj=idptl(j) c endif c c do k=1,nflav c kc(k)=jc(k,1)-jc(k,2) c if(nq.lt.0)kc(k)=-kc(k) c enddo c if(nq.lt.0)nq=-nq c c if(ish.ge.7)then c write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j) c write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6) c write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6) c endif c cc check minimal kinetic energy c ekin=ecm-pptl(5,i)-pptl(5,j) c if(ekin.lt.amimel)goto1000 c cc ------------------------------------------------------------------------- c if(iabs(idi).gt.1000.or.iabs(idj).gt.1000)then !baryon involved cc ------------------------------------------------------------------------- c c if(nq.eq.6)then !------------baryon-baryon -----> c c if(kc(1).eq.kc(2))then !isospin_z=0 cc pn c if(ish.ge.7)write(ifch,*)'sig_pn chosen' c if(ecm.lt.2)then c call utindx(npn,pnecm,ecm,iecm) c bmx=pnbmx(iecm) c else c p=ecm*sqrt((ecm**2/4./.94**2)-1.) c if(p.lt.1.)then c sig=33.+196.*(abs(0.95-p))**2.5 c elseif(p.lt.2)then c sig=24.2+8.9*p c else c sig=47.3+.513*(alog(p)**2)-4.27*alog(p) c endif c bmx=sqrt(sig/10./pi) c endif c else !isospin_z=+-1 cc pp c if(ish.ge.7)write(ifch,*)'sig_pp chosen' c if(ecm.le.1.882)then c bmx=ppbmx(1) c if(ish.ge.7)write(ifch,*)'b_mx:',bmx c elseif(ecm.le.1.887)then c bmx=ppbmx(2) c if(ish.ge.7)write(ifch,*)'b_mx:',bmx c elseif(ecm.le.1.893)then c bmx=ppbmx(3) c if(ish.ge.7)write(ifch,*)'b_mx:',bmx c elseif(ecm.le.1.899)then c bmx=ppbmx(4) c elseif(ecm.le.1.909)then c bmx=ppbmx(5) c elseif(ecm.le.1.913)then c bmx=ppbmx(6) c elseif(ecm.le.1.918)then c bmx=ppbmx(7) c else c p=ecm*sqrt((ecm**2/4./.94**2)-1.) c if(p.lt.0.8)then c sig=23.5+1000.*(0.7-p)**4 c elseif(p.lt.1.5)then c sig=23.5+24.6/(1+exp(-(p-1.2)/0.1)) c elseif(p.lt.5)then c sig=41.+60.*(p-0.9)*exp(-1.2*p) c else c sig=48+.522*(alog(p)**2)-4.51*alog(p) c endif c bmx=sqrt(sig/10./pi) c endif c endif c c elseif(nq.eq.0)then !-----------baryon-antibaryon ----> c c if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0)then !p-ap, n-an cc app c if(ish.ge.7)write(ifch,*)'sig_app chosen' c if(ecm.lt.2)then c call utindx(napp,appecm,ecm,iecm) c bmx=appbmx(iecm) c else c p=ecm*sqrt((ecm**2/4./.94**2)-1.) c sig=38.4+77.6*p**(-.64)+.26*(alog(p)**2)-1.2*alog(p) c bmx=sqrt(sig/10./pi) c endif c else !p-an, n-ap cc apn c if(ish.ge.7)write(ifch,*)'sig_apn chosen' c if(ecm.lt.2)then c call utindx(napn,apnecm,ecm,iecm) c bmx=apnbmx(iecm) c else c p=ecm*sqrt((ecm**2/4./.94**2)-1.) c sig=133.6*p**(-.7)-1.22*(alog(p)**2)+13.7*alog(p) c bmx=sqrt(sig/10./pi) c endif c endif c c elseif(nq.eq.3)then !----------baryon-meson ----> c c if(kc(3).eq.0)then cc no kaons involved(except for K-L) c if(kc(1).eq.0.or.kc(2).eq.0)then cc pip c if(ish.ge.7)write(ifch,*)'sig_pip chosen' c if(ecm.lt.2.5)then c call utindx(npip,pipecm,ecm,iecm) c bmx=pipbmx(iecm) c else c p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2) c sig=16.4+19.3*p**(-.42)+.19*(alog(p)**2) c bmx=sqrt(sig/10./pi) c endif c else cc pim c if(ish.ge.7)write(ifch,*)'sig_pim chosen' c if(ecm.lt.2.5)then c call utindx(npim,pimecm,ecm,iecm) c bmx=pimbmx(iecm) c else c p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2) c sig=33.0+14.0*p**(-1.36)+.456*(alog(p)**2)-4.03*alog(p) c bmx=sqrt(sig/10./pi) c endif c endif c elseif(kc(3).eq.1)then cc strange particles involved c if(kc(1).eq.0.or.kc(2).eq.0)then cc kmn c if(ish.ge.7)write(ifch,*)'sig_kmn chosen' c if(ecm.lt.2.5)then c call utindx(nkmn,kmnecm,ecm,iecm) c bmx=kmnbmx(iecm) c else c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) c sig=25.2+.38*(alog(p)**2)-2.9*alog(p) c bmx=sqrt(sig/10./pi) c endif c else cc kmp c if(ish.ge.7)write(ifch,*)'sig_kmp chosen' c if(ecm.lt.2.5)then c call utindx(nkmp,kmpecm,ecm,iecm) c bmx=kmpbmx(iecm) c else c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) c sig=32.1+.66*(alog(p)**2)-5.6*alog(p) c bmx=sqrt(sig/10./pi) c endif c endif c elseif(kc(3).eq.-1)then cc strange particles involved c if(kc(1).eq.3.or.kc(2).eq.3)then cc kpp c if(ish.ge.7)write(ifch,*)'sig_kpp chosen' c if(ecm.lt.2.5)then c call utindx(nkpp,kppecm,ecm,iecm) c bmx=kppbmx(iecm) c else c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) c sig=18.1+.26*(alog(p)**2)-alog(p) c bmx=sqrt(sig/10./pi) c endif c else cc kpn c if(ish.ge.7)write(ifch,*)'sig_kpn chosen' c if(ecm.lt.2.5)then c call utindx(nkpn,kpnecm,ecm,iecm) c bmx=kpnbmx(iecm) c else c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) c sig=18.7+.21*(alog(p)**2)-.89*alog(p) c bmx=sqrt(sig/10./pi) c endif c endif c elseif(kc(3).ge.2)then cc two strange particles involved c bmx=1. c elseif(kc(3).le.-2)then cc two strange particles involved c bmx=1. c endif c c else !-----------baryon_cluster-anything----> c c write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j) c write(ifch,*)'r_i:',radptl(i),' r_j:',radptl(j) c write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6) c write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6) c bmx=radptl(i)+radptl(j) c c endif ! <-------------- c cc ------------------------------------- c else ! meson-meson cc ------------------------------------- c c call idquad(i,nqi,nai,jci) c call idquad(j,nqj,naj,jcj) c do l=1,nflav c kci(l)=jci(l,1)-jci(l,2) c kcj(l)=jcj(l,1)-jcj(l,2) c enddo c c if(kci(3).eq.0.and.pptl(5,i).le.1.0 c * .and.kcj(3).eq.0.and.pptl(5,j).le.1.0)then c c if(kci(1).eq.0.and.kci(2).eq.0.and.pptl(5,i).le.0.5)then c if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then cc pi0 pi0 c goto104 c elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then cc pi0 pi+ c goto102 c elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then cc pi0 pi- c goto103 c else cc pi0 eta c goto105 c endif c elseif(kci(1).eq.1.and.kci(2).eq.-1)then c if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then cc pi+ pi0 c goto102 c elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then cc pi+ pi+ c goto101 c elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then cc pi+ pi- c goto104 c else cc pi+ eta c goto105 c endif c elseif(kci(1).eq.-1.and.kci(2).eq.1)then c if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then cc pi- pi0 c goto103 c elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then cc pi- pi+ c goto104 c elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then cc pi- pi- c goto101 c else cc pi- eta c goto105 c endif c else c if(pptl(5,j).le.0.5)then cc eta pi c goto105 c else cc eta eta c goto106 c endif c endif c101 continue c if(ish.ge.7)write(ifch,*)'sig_pi+pi+ chosen' c call utindx(npi1,pi1ecm,ecm,iecm) c bmx=pi1bmx(iecm) c goto110 c102 continue c if(ish.ge.7)write(ifch,*)'sig_pi+pi0 chosen' c call utindx(npi2,pi2ecm,ecm,iecm) c bmx=pi2bmx(iecm) c goto110 c103 continue c if(ish.ge.7)write(ifch,*)'sig_pi-pi0 chosen' c call utindx(npi3,pi3ecm,ecm,iecm) c bmx=pi3bmx(iecm) c goto110 c104 continue c if(ish.ge.7)write(ifch,*)'sig_pi-pi+ chosen' c call utindx(npi4,pi4ecm,ecm,iecm) c bmx=pi4bmx(iecm) c goto110 c105 continue c if(ish.ge.7)write(ifch,*)'sig_pi_eta chosen' c call utindx(npi5,pi5ecm,ecm,iecm) c bmx=pi5bmx(iecm) c goto110 c106 continue c if(ish.ge.7)write(ifch,*)'sig_eta_eta chosen' c bmx=0.4 ! approx. sig=5mb c110 continue c c else !something else involved, strange etc. c c bmx=0.6 ! approx. sig=10mb c c endif c cc -------------------------------- c endif cc -------------------------------- c c if(bij.le.bmx)ics=1 c if(ish.ge.7)then c write(ifch,*)'b_mx:',bmx,' b_ij:',bij,' ics:',ics c endif c c1000 continue c call utprix('jintcs',ish,ishini,7) c return c end c cc---------------------------------------------------------------------- c subroutine jintel(i,j,p,amim,xaver) cc---------------------------------------------------------------------- cc elastic scattering of ptls i,j cc---------------------------------------------------------------------- c include 'epos.inc' c real xaver(4) c real p(5),u(5),pei(5),pej(5) c cc determine momenta of outgoing particles (pei,pej) cc ------------------------------------------------- c if(p(5).le.(pptl(5,i)+pptl(5,j))*.99)then c if(ish.ge.1)then c call utmsg('jintel') c write(ifch,132)p(5),pptl(5,i)+pptl(5,j) c132 format(1x,'***** m_fus < m_i+m_j ---> qcm set zero ( ' c *,2f10.3,' )') c write(ifch,133)'p_i: ',(pptl(k,i),k=1,5) c write(ifch,133)'p_j: ',(pptl(k,j),k=1,5) c write(ifch,133)'p_fus:',(p(k),k=1,5) c133 format(1x,a6,3x,5f10.3) c call utmsgf c endif c qcm=0. c elseif(p(5).le.pptl(5,i)+pptl(5,j))then c qcm=0. c else c qcm=utpcm(p(5),pptl(5,i),pptl(5,j)) c endif c cc isotropic c u(3)=2.*rangen()-1. c phi=2.*pi*rangen() c u(1)=sqrt(1.-u(3)**2)*cos(phi) c u(2)=sqrt(1.-u(3)**2)*sin(phi) c do 47 k=1,3 c pei(k)= qcm*u(k) c47 pej(k)=-qcm*u(k) c cc nonisotropic cc-c pt=sqrt(2./pi)*ptq*sqrt(-2*alog(rangen()) !2-dim Gauss cc-c if(pt.ge.qcm)pt=rangen()*qcm cc-c qpl=sqrt(qcm**2-pt**2) cc-c u(3)=qpl cc-c phi=2.*pi*rangen() cc-c u(1)=pt*cos(phi) cc-c u(2)=pt*sin(phi) cc-c call utaxis(i,j,a1,a2,a3) cc-c ivt=1 cc-c if(a3.lt.0.)then cc-c a1=-a1 cc-c a2=-a2 cc-c a3=-a3 cc-c ivt=-1 cc-c endif cc-c call utrota(-1,a1,a2,a3,u(1),u(2),u(3)) cc-c do 47 k=1,3 cc-c pei(k)= u(k)*ivt cc-c47 pej(k)=-u(k)*ivt c c pei(4)=sqrt(qcm**2+pptl(5,i)**2) c pej(4)=sqrt(qcm**2+pptl(5,j)**2) c pei(5)=pptl(5,i) c pej(5)=pptl(5,j) c call utlobo(-1,p(1),p(2),p(3),p(4),p(5) c *,pei(1),pei(2),pei(3),pei(4)) c call utlobo(-1,p(1),p(2),p(3),p(4),p(5) c *,pej(1),pej(2),pej(3),pej(4)) c cc fill /cptl/ cc ----------- c do 49 lo=1,2 c nptl=nptl+1 c if(lo.eq.1)ij=i c if(lo.eq.2)ij=j c do 48 k=1,5 c if(lo.eq.1)pptl(k,nptl)=pei(k) c if(lo.eq.2)pptl(k,nptl)=pej(k) c48 continue c istptl(nptl)=0 c idptl(nptl)=idptl(ij) c ibptl(1,nptl)=ibptl(1,ij) c ibptl(2,nptl)=ibptl(2,ij) c ibptl(3,nptl)=ibptl(3,ij) c ibptl(4,nptl)=ibptl(4,ij) c xorptl(1,nptl)=xaver(1) c xorptl(2,nptl)=xaver(2) c xorptl(3,nptl)=xaver(3) c xorptl(4,nptl)=xaver(4) c iorptl(nptl)=i c jorptl(nptl)=j c tivptl(1,nptl)=xaver(4) c call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm) c tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen())) c ifrptl(1,nptl)=0 c ifrptl(2,nptl)=0 c ityptl(nptl)=ityptl(ij) c49 continue c c1000 return c end c cc---------------------------------------------------------------------- c subroutine jintfs(i,j,p,nq,jc,amim,iret) cc---------------------------------------------------------------------- cc input: cc i,j: particle indices cc output: cc p: 5-momentum of fused ptl cc nq: net quark number of fused ptl cc jc: jc of fused ptl cc amim : minimum mass of fused ptl cc iret: return code cc 0=ok cc 1=mass p(5) less than amim cc---------------------------------------------------------------------- c include 'epos.inc' c integer jci(nflav,2),jcj(nflav,2),jc(nflav,2) c real p(5) c double precision ppfus(4),pp52 c c iret=0 c c do 35 k=1,4 c p(k)=pptl(k,i)+pptl(k,j) c35 ppfus(k)=dble(p(k)) c pp52=ppfus(4)**2-ppfus(3)**2-ppfus(2)**2-ppfus(1)**2 c if(pp52.le.0.)then c if(ish.ge.1)then c call utmsg('jintfs') c write(ifch,*)'***** mfus**2 < 0 (',pp52,' )' c write(ifch,*)(ppfus(m),m=1,4) c call utmsgf c endif c goto1001 c endif c p(5)=sngl(dsqrt(pp52)) c c call idquad(i,nqi,nai,jci) c call idquad(j,nqj,naj,jcj) c do 29 n=1,nflav c do 29 k=1,2 c29 jc(n,k)=jci(n,k)+jcj(n,k) c nq=0 c do 54 n=1,nflav c54 nq=nq+jc(n,1)-jc(n,2) c c call idcomj(jc) c amim=utamnz(jc,5)+amimfs c if(p(5).lt.amim)goto1001 c goto1000 c1001 iret=1 c1000 return c end c cc---------------------------------------------------------------------- c subroutine jintfu(i,j,p,jc) cc---------------------------------------------------------------------- cc input: cc i,j: particle indices cc p,jc: momentum and jc of fused object cc outout: cc id, ib, ist, ior, jor, ifr, ity of fused particle cc written to /cptl/ after increasing nptl by 1 cc---------------------------------------------------------------------- c c include 'epos.inc' c parameter (mxlook=10000,mxdky=2000) c common/dkytab/look(mxlook),cbr(mxdky),mode(5,mxdky) c real p(5) c integer jc(nflav,2),ic(2),ib(4) c c nptl=nptl+1 c cc momentum c do k=1,5 c pptl(k,nptl)=p(k) c enddo c amf=p(5) c cc determine idr, ib(1-4) c idr=0 c do 40 nf=1,nflav c do 40 ij=1,2 c if(jc(nf,ij).ge.10)idr=7*10**8 c40 continue c if(idr/10**8.ne.7)then c call idenco(jc,ic,ireten) c if(ireten.eq.1)call utstop('jintfu: idenco ret code = 1&') c id=idtra(ic,0,0,3) c43 amc=amf c call idres(id,amc,idr,iadj) c if(idr.ne.0)then c lid=look(iabs(idr)) c if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0) c *.and.pptl(5,nptl).gt.amc+1e-3)then c amf=amf+0.010 c goto43 c endif c if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0) c *.and.abs(amc-pptl(5,nptl)).gt.1e-3)then c if(ish.ge.1)then c call utmsg('jintfu') c write(ifch,*)'***** not on mass shell after fusion: ' c *,pptl(5,nptl),amc c call utmsgf c endif c endif c endif c if(idr.eq.0)then c if(mod(ic(1),100).ne.0.or.mod(ic(2),100).ne.0)then c idr=9*10**8 c else c idr=8*10**8+ic(1)*100+ic(2)/100 c endif c endif c else c call idtrbi(jc,ib(1),ib(2),ib(3),ib(4)) c idr=idr c *+mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4 c *+mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4) c if(ish.ge.7)write(ifch,*)'ib:',(ib(kk),kk=1,4) c endif c c n=nptl c cc fill /cptl/ cc ----------- c idptl(n)=idr c do k=1,4 c ibptl(k,n)=ib(k) c enddo c istptl(n)=0 c iorptl(n)=i c jorptl(n)=j c ifrptl(1,n)=0 c ifrptl(2,n)=0 c ityptl(n)=50 c c return c end c c---------------------------------------------------------------------- subroutine jintpo(lclean,iret) c---------------------------------------------------------------------- c parton-ladder-fusion -- 3D grid -- based on string segments c---------------------------------------------------------------------- include 'epos.inc' include 'epos.incems' include 'epos.incico' common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl) *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3) parameter(m1grid=10,kgrid=3,kegrid=3) parameter(m3xgrid=21*kgrid*kegrid) parameter(mxcl=4000,mxcli=1000) integer idropgrid(m1grid,m1grid,m3xgrid) & ,jdropgrid(m1grid,m1grid,m3xgrid) & ,jclu(m3xgrid),nsegp4(mxcl) & ,irep(mxcl),mmji(mxcl,mxcli) & ,jccl(mxcl,nflav,2),nseg(mxcl),mseg(mxcl),kclu(mxcl) & ,naseg(0:mxcl),nfseg(mxcl),nsegmx(mxcl) & ,jc(nflav,2),ke(6),jcjj(nflav,2) & ,nsegmt(mxptl),kmean(2,1001:1000+mxcl) & ,ncluj(1001:1000+mxptl) !,ic(2),nclk(m3xgrid) double precision tpro,zpro,ttar,ztar,ttaus,detap,detat & ,pptld(5,mxcl),p4tmp double precision enp,pzp,ena,pza,ccc,sss,enew,pznew,taa2,pta2 .,pold,pnew,pnew2,p5a2,eeta,enf,pzf common /cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat common/cdelzet/delzet,delsgr /cvocell/vocell common/cranphi/ranphi parameter(maxp=40*mxcl) dimension yrad(maxp),phirad(maxp),pe(5),yrad2(maxp),nfrag(mxptl) logical first,lnew(m1grid,m1grid),lold(m1grid,m1grid),lcont,lpass &,lclean double precision ptest(5),ttest,p52,p4mean(4,mxcl),zor,tor,xmxms &,ppp(5),be(4),energ,bp,rmean(2,mxcl) !,qqq(5),www(5) dimension mpair(mamx,mamx) save ncntpo!,ntrr,ntrt data ncntpo /0/!,ntrr/0/,ntrt/0/ ncntpo=ncntpo+1 call utpri('jintpo',ish,ishini,4) tauico=ttaus if(kigrid.gt.kgrid)then write(ifmt,*)'kigrid,kgrid:',kigrid,kgrid stop'jintpo: kigrid too big. ' endif fegrid=yhaha/5.36 !rapidity range/rap range at RHIC if(fegrid.gt.kegrid)then write(ifmt,*)'fegrid,kegrid:',fegrid,kegrid stop'jintpo: kegrid too small. ' endif m3grid=m3xgrid /kgrid *kigrid /kegrid *fegrid c if(iLHC.eq.1)m3grid=min(m3xgrid,int(m3grid*0.75)) if(mod(m3grid,2).eq.0)m3grid=m3grid+1 xgrid=8 sgrid=fsgrid *ttaus *fegrid *6.0 delxgr=2*xgrid/m1grid delsgr=2*sgrid/m3grid shiftx=(1.-2.*rangen())*rangen()*delxgr/2. shifts=(1.-2.*rangen())*rangen()*delsgr/2. vocell= delxgr * delxgr * delsgr delzet=delsgr/ttaus nptla=nptl iflhc=1 if(iLHC.eq.1)iflhc=min(3,max(1,nsegce-1)) nsegsuj=max(nsegce,nint(float(nsegsu/iflhc)*fegrid)) xmxms=150d0 !maximum mass for a subcluster if(ncntpo.eq.1)then write(ifmt,*)'EPOS used with FUSION option' if(ish.ge.1)then print*,'+',('-',i=1,57) print*,'| ttaus sgrid nsegce:',ttaus,fsgrid,' ',nsegce print*,'| sgrid delxgr delsgr vocell:',sgrid,delxgr,delsgr,vocell print*,'+',('-',i=1,57) endif endif c print*,'+',('-',i=1,57) c...compute x,y,z if(ish.ge.6)write(ifch,*)'compute x,y,z' do n=1,nptla iaaptl(n)=1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! the following check is important for ioclude= 4 or 5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(ioclude.ne.3.and.istptl(n).eq.10)then stop'\n\n remnant cluster detected in jintpo\n\n' endif if(istptl(n).eq.0.or.istptl(n).eq.10)then call jtain(n,xptl(n),yptl(n),zptl(n),tptl(n),nnn,0) call idtau(idptl(n),pptl(4,n),pptl(5,n),taurem) if(abs(taurem).gt.1e-6.and. . ityptl(n).ge.40.and.ityptl(n).le.59.and.nnn.eq.1)then iaaptl(n)=0 !!!!nnn=1: ptl lives later c print*,n,ityptl(n),idptl(iorptl(n)),istptl(iorptl(n)),nnn endif call jtaus(zptl(n),tz,sptl(n)) strap=1e10 xpl=tptl(n)+zptl(n) xmi=tptl(n)-zptl(n) if(xmi.gt.0.0.and.xpl.gt.0.0)then strap=0.5*log(xpl/xmi) else !particle at eta=inf iaaptl(n)=0 endif dezptl(n)=strap !space-time-rapidity c write(ifch,*)'dez ???',n,nnn,taurem,ityptl(n) c & ,pptl(4,n),pptl(3,n),idptl(n),xptl(n),yptl(n),zptl(n),sptl(n) else iaaptl(n)=0 xptl(n)=0. yptl(n)=0. zptl(n)=0. sptl(n)=0. dezptl(n)=1e10 endif enddo ntry=0 c...valid particles c print *,'start -----------------------------------------' ptmax=ptclu**2 nptlo=1 if(iLHC.eq.1)nptlo=maproj+matarg+1 do n=1,nptla !random number calls before if's !to keep random number sequence for testing rdm=rangen() if(istptl(n).eq.0)then call idquac(n,idum1,idum2,idum3,jc) else jc(4,1)=0 jc(4,2)=0 endif if(iaaptl(n).ne.0)then if(istptl(n).ne.10)then pt2=pptl(1,n)**2+pptl(2,n)**2 am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2 if(n.le.nptlo)then !spectators iaaptl(n)=0 elseif(ityptl(n).eq.47.or.ityptl(n).eq.57)then iaaptl(n)=0 !no active spectators elseif(istptl(n).ne.0.or.ityptl(n).ge.60)then iaaptl(n)=0 elseif(abs(idptl(n)).le.100)then iaaptl(n)=0 !no fundamental particle (electron...) !elseif(ityptl(n).eq.41.or.ityptl(n).eq.51)then ! iaaptl(n)=0 !avoid particles coming from remnant droplet decay ! !(already droplet decay products !) elseif(abs(am2tmp-pptl(5,n)*pptl(5,n)).gt.30.)then iaaptl(n)=0 !to discard off shell particles elseif(jc(4,1).ne.0.or.jc(4,2).ne.0)then iaaptl(n)=0 !to discard charmed particles elseif(iLHC.eq.1)then if(pt2.lt.1.e-3)then iaaptl(n)=0 !to discard slow proton (spectators) elseif(ptclu.gt.0.)then if(pt2.gt.(1.5*ptclu)**2)then ptmax=max(ptmax,pt2) if(ioquen.eq.0.or.fploss.le.0.)then !high pt particle escape iaaptl(n)=0 elseif(fploss.gt.1.e10)then !high pt particle absorbed iaaptl(n)=1 else iaaptl(n)=-1 !high pt particle lose energy endif elseif(pt2.gt.(0.5*ptclu)**2)then if(rdm.lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then ptmax=max(ptmax,pt2) if(ioquen.eq.0.or.fploss.le.0.)then !high pt particle escape iaaptl(n)=0 elseif(fploss.gt.1.e10)then !high pt particle absorbed iaaptl(n)=1 else iaaptl(n)=-1 !high pt particle lose energy endif endif endif else if(ioquen.eq.0.or.fploss.le.0.)then !high pt particle escape iaaptl(n)=0 elseif(fploss.gt.1.e10)then !high pt particle absorbed iaaptl(n)=1 else iaaptl(n)=-1 !high pt particle lose energy endif endif elseif(pt2.gt.(1.5*ptclu)**2)then ptmax=max(ptmax,pt2) if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)then iaaptl(n)=0 endif elseif(pt2.gt.(0.5*ptclu)**2)then if(rdm.lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)iaaptl(n)=0 endif endif endif endif c if(iaaptl(n).eq.0.and.abs(idptl(n)).le.10000.and.abs(idptl(n)) c & .ge.100.and.mod(abs(idptl(n)),100).ne.0 c & .and.ityptl(n).ne.0.and.ityptl(n).lt.50) c & print*,'not selected',idptl(n),ityptl(n),pptl(4,n) enddo c...Start cluster formation 8888 continue nsegsuj=max(nsegce,nint(nsegsuj/1.08**ntry)) ntry=ntry+1 if(ntry.gt.90) !do not put more than 100 or change limit for p4mean &call utstop('jintpo: cluster formation impossible ! &') nptl=nptla do 1 k=1,m3grid do 1 j=1,m1grid do 1 i=1,m1grid idropgrid(i,j,k)=0 1 continue c...count string segments in cell if(ish.ge.6)write(ifch,*)'count string segments in cell' do n=1,nptla if(iaaptl(n).ne.0)then i=1+(xptl(n)+xgrid+shiftx)/delxgr j=1+(yptl(n)+xgrid+shiftx)/delxgr k=1+(sptl(n)+sgrid+shifts)/delsgr if( i.ge.1.and.i.le.m1grid & .and.j.ge.1.and.j.le.m1grid & .and.k.ge.1.and.k.le.m3grid)then nfrag(n)=1 if(iLHC.eq.1)then !count number of quarks to absorbe more baryon than mesons call idquac(n,idum1,idum2,idum3,jc) nfrag(n)=0 do nj=1,2 do ni=1,nflav nfrag(n)=nfrag(n)+ni*jc(ni,nj) enddo enddo endif c print *,idptl(n),nfrag(n) idropgrid(i,j,k)=idropgrid(i,j,k)+nfrag(n) endif endif enddo if(iLHC.eq.1)then cc...check high pt segments c !to use this part one has to define: c !... 1 = valid c !... -1 = valid but high pt c !... 0 = not valid c c esu=0 c do i=1,nptla c if(istptl(i).eq.0.or.istptl(i).eq.3)esu=esu+pptl(4,i) c enddo c write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu if(ish.ge.6)write(ifch,*)'check high pt segments' if(fploss.gt.0.0)then ein=0 elo=0 do n=1,nptla delen=0. if(iaaptl(n).eq.-1)then i=1+(xptl(n)+xgrid+shiftx)/delxgr j=1+(yptl(n)+xgrid+shiftx)/delxgr k=1+(sptl(n)+sgrid+shifts)/delsgr if( i.ge.1.and.i.le.m1grid & .and.j.ge.1.and.j.le.m1grid & .and.k.ge.1.and.k.le.m3grid)then iescape=0 xa=xptl(n) ya=yptl(n) eta=dezptl(n) pz=pptl(3,n) en=pptl(4,n) taa2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2 pza=pz ena=sqrt(taa2+pza**2) enaxx=en p5a2=pptl(5,n)**2 pta2=pptl(1,n)**2+pptl(2,n)**2 !transform to bjo eeta=eta ccc=cosh(eeta) sss=sinh(eeta) enp= ena*ccc-pza*sss pzp=-ena*sss+pza*ccc vz=pzp/enp vx=pptl(1,n)/enp vy=pptl(2,n)/enp !print*,'+++++++',eta,vz,pz/en if(vx.ne.0.0.or.vy.ne.0.0)then if(abs(vx).ge.abs(vy))then ica=1 rat=vy/vx is=sign(1.,vx) l=i va=xa wa=ya else ica=2 rat=vx/vy is=sign(1.,vy) l=j va=ya wa=xa endif if(is.eq.-1)then imax=1 else !if(is.eq. 1) imax=m1grid endif vr=sqrt(vx**2+vy**2) ratx=vz/vr qu=fploss*delxgr*sqrt(1+rat**2)*sqrt(1+ratx**2) do lu=l,imax,is vu=-(xgrid+shiftx)+(lu-0.5)*delxgr wu=wa+rat*(vu-va) mu=1+(wu+xgrid+shiftx)/delxgr if(mu.ge.1.and.mu.le.m1grid)then if(ica.eq.1)then ix=lu jx=mu else ix=mu jx=lu endif if(idropgrid(i,j,k).ge.nsegce)then dens=float(idropgrid(ix,jx,k))/float(iflhc) else dens=0 endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! DeltaE ~ sqrt(T**3*E) * dL (BDMPS2008) ! -> DeltaE ~ dens^3/8 * sqrt(E) *dL !del=dens**(3./8.)*max(1.,sqrt(pptl(4,n))) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ escale=6 xen=enp/escale del=0.02*dens**(3./8.)*max(1.,sqrt(xen)) delen=delen+del*escale*qu c if(k.ge.m3grid/2-2.and.k.le.m3grid/2+4) c & write(ifch,'(2i3,4x,2i3,4x,2i3,4x,i3,3f7.3)') c & k,ica,i,j,ix,jx,idropgrid(ix,jx,k),qu,delen c & ,pptl(4,n)-delen endif enddo endif pold=sqrt(pta2+pzp**2) enew=enp-delen pnew2=enew**2-p5a2 if(enew.gt.0..and.pnew2.gt.0.)iescape=1 if(iescape.eq.1)then iaaptl(n)=0 idropgrid(i,j,k)=idropgrid(i,j,k)-nfrag(n) if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible. ' pnew=sqrt(pnew2) p1=pnew*pptl(1,n)/pold p2=pnew*pptl(2,n)/pold pznew=pnew*pzp /pold enf= enew*ccc+pznew*sss pzf= enew*sss+pznew*ccc if(ish.ge.1.and.enf.gt.ena.or.enf.ne.enf)then if(ish.ge.1.and.(enf-ena.gt.0.1d0.or.enf.ne.enf))then write(ifmt,*)'WWWW en gain for eta =',eta write(ifmt,*)'e pz bf (lab)',ena,pza,enaxx write(ifmt,*)'e pz bf (bjo)',enp,pzp ,enp**2-pzp**2 write(ifmt,*)'e pz af (bjo)',enew,pznew,enew**2-pznew**2 write(ifmt,*)'e pz af (lab)',enf,pzf ,enf**2-pzf**2 c stop endif elseif(delen.gt.1e-3)then c create fake particle with energy lost in cluster nptl=nptl+1 nptla=nptla+1 istptl(nptl)=3 !daughter part of cluster iaaptl(nptl)=1 pptl(1,nptl)=pptl(1,n)-p1 pptl(2,nptl)=pptl(2,n)-p2 pptl(3,nptl)=pptl(3,n)-pzf pptl(4,nptl)=pptl(4,n)-enf do l=1,3 xorptl(l,nptl)=xorptl(l,n) enddo ! mother pptl(1,n)=p1 pptl(2,n)=p2 pptl(3,n)=pzf tivptl(1,n)=xorptl(4,n) pptl(4,n)=sqrt(pptl(1,n)**2+pptl(2,n)**2 & +pptl(3,n)**2+pptl(5,n)**2) if(abs(pptl(4,n)-enf).gt.0.01*pptl(4,n)) & print *,'jintpo eloss',pptl(4,n),enf,(pptl(4,n)-enf)/pptl(4,n) call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm) tivptl(2,n)=tivptl(1,n)+taugm*(-alog(rangen())) ifrptl(1,n)=nptl ! daughter pptl(5,nptl)=0. pptl(4,nptl)=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2 & +pptl(3,nptl)**2+pptl(5,nptl)**2) xorptl(4,nptl)=xorptl(4,n) tivptl(1,nptl)=tivptl(1,n) iorptl(nptl)=n jorptl(nptl)=0 idptl(nptl)=idptl(n)*100+sign(99,idptl(n)) ityptl(nptl)=ityptl(n)+2 xptl(nptl)=xptl(n) yptl(nptl)=yptl(n) zptl(nptl)=zptl(n) sptl(nptl)=sptl(n) dezptl(nptl)=dezptl(n) if(ish.ge.6)write(ifch,*)'--> Virtual particle : ',nptl & ,idptl(nptl),iorptl(nptl),ityptl(nptl),(pptl(kk,nptl),kk=1,5) endif elo= elo+ena-pptl(4,n) !write(ifch,*)n,' escape' else iaaptl(n)=1 ein=ein+pptl(4,n) !write(ifch,*)n,'core' endif endif endif enddo else do n=1,nptla if(iaaptl(n).eq.-1)then iaaptl(n)=0 endif enddo endif c eloss=esu c esu=0 c do i=1,nptla c if(istptl(i).eq.0.or.istptl(i).eq.3)esu=esu+pptl(4,i) c enddo c write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu c eloss=eloss-esu c write(ifch,*)'+++++ein,elo,eloss ',ein,elo,eloss c do n=1,nptla c if(iaaptl(n).eq.-1)then c i=1+(xptl(n)+xgrid+shiftx)/delxgr c j=1+(yptl(n)+xgrid+shiftx)/delxgr c k=1+(sptl(n)+sgrid+shifts)/delsgr c if( i.ge.1.and.i.le.m1grid c & .and.j.ge.1.and.j.le.m1grid c & .and.k.ge.1.and.k.le.m3grid)then c if(idropgrid(i,j,k).lt.2*nsegce)then !surface segments c iaaptl(n)=0 c idropgrid(i,j,k)=idropgrid(i,j,k)-1 c if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible. ' cc else cc iaaptl(n)=1 c endif c endif c endif c enddo endif c...identify clusters ifirst=0 if(ish.ge.6)write(ifch,*)'identify clusters' do k=1,m3grid !~~~~~~k-loop jjj=0 do j=1,m1grid first=.true. do i=1,m1grid if(idropgrid(i,j,k).ge.nsegce)then if(first)then ifirst=i jjj=jjj+1 if(jjj.gt.mxcl)stop'jintpo: mxcl too small. ' irep(jjj)=0 jj=jjj first=.false. endif jdropgrid(i,j,k)=jj if(j.gt.1)then if(jdropgrid(i,j-1,k).ne.0)then jjo=jdropgrid(i,j-1,k) if(jjo.lt.jj)then if(jj.eq.jjj)jjj=jjj-1 jjx=jj jj=jjo do ii=ifirst,i jdropgrid(ii,j,k)=jj if(jdropgrid(ii,j-1,k).eq.jjx)then if(jjx.gt.jjj)jjj=jjj+1 jja=jjx jjb=jj 90 continue if(irep(jja).eq.0.or.irep(jja).eq.jjb)then irep(jja)=jjb else mn=min(irep(jja),jjb) mx=max(irep(jja),jjb) irep(jja)=mn jja=mx jjb=mn goto 90 endif endif enddo elseif(jdropgrid(i,j-1,k).gt.jj)then irep(jjo)=jj endif endif endif else jdropgrid(i,j,k)=0 first=.true. endif enddo enddo !~~~~cluster jj ---> cluster irep(jj) do jj=jjj,1,-1 if(irep(jj).ne.0)then do j=1,m1grid do i=1,m1grid if(jdropgrid(i,j,k).eq.jj)jdropgrid(i,j,k)=irep(jj) enddo enddo endif enddo !~~~~~remove empty cluster indices jjjx=jjj jjj=0 jj=0 do jjx=1,jjjx if(irep(jjx).eq.0)then jj=jj+1 jjj=jjj+1 else do j=1,m1grid do i=1,m1grid if(jdropgrid(i,j,k).gt.jj) & jdropgrid(i,j,k)=jdropgrid(i,j,k)-1 enddo enddo endif enddo !~~~~~ jclu(k)=jjj enddo !~~~~~~~~~~~~~~~~~ END k-loop c...add segments to avoid holes if(ish.ge.6)write(ifch,*)'add segments from holes' if(iohole.eq.1.and.iLHC.eq.1)then do 83 n=1,nptla if(iaaptl(n).eq.0)goto 83 ihit=0 i=1+(xptl(n)+xgrid+shiftx)/delxgr j=1+(yptl(n)+xgrid+shiftx)/delxgr k=1+(sptl(n)+sgrid+shifts)/delsgr if( i.ge.1.and.i.le.m1grid & .and.j.ge.1.and.j.le.m1grid & .and.k.ge.1.and.k.le.m3grid)then jgr=jdropgrid(i,j,k) nplus=idropgrid(i,j,k) if(jgr.eq.0)then isgi=1 if(rangen().gt.0.5)isgi=-1 isgj=1 if(rangen().gt.0.5)isgj=-1 isgk=1 if(rangen().gt.0.5)isgk=-1 do ii=i-isgi,i+isgi,isgi do jj=j-isgj,j+isgj,isgj do kk=k-isgk,k+isgk,2*isgk if( ii.ge.1.and.ii.le.m1grid . .and.jj.ge.1.and.jj.le.m1grid . .and.kk.ge.1.and.kk.le.m3grid)then if(jdropgrid(ii,jj,kk).gt.0)nplus=nplus+1 c nplus=idropgrid(ii,jj,kk)+iflhc if(nplus.ge.nsegce)then ihit=1 c goto 84 endif endif enddo enddo enddo endif endif c 84 continue if(ihit.eq.1)then idropgrid(i,j,k)=nplus!idropgrid(i,j,k)+iflhc c add cell randomly to adjacent cluster if any isgi=1 if(rangen().gt.0.5)isgi=-1 isgj=1 if(rangen().gt.0.5)isgj=-1 do ii=i-isgi,i+isgi,2*isgi do jj=j-isgj,j+isgj,2*isgj if( ii.ge.1.and.ii.le.m1grid . .and.jj.ge.1.and.jj.le.m1grid)then jjj=jdropgrid(ii,jj,k) if(jjj.gt.0)then jdropgrid(i,j,k)=jjj goto 83 endif endif enddo enddo c if no adjacent cluster, create a new one jclu(k)=jclu(k)+1 jdropgrid(i,j,k)=jclu(k) endif 83 continue endif if(ish.ge.8)then write(ifch,*)'segment positions' do k=1,m3grid write(ifch,*)'k=',k,' jclu=',jclu(k) & ,' s=',(k-1)*delsgr-sgrid-shifts do j=m1grid,1,-1 write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,m1grid) & ,(jdropgrid(i,j,k),i=1,m1grid) enddo enddo endif c...absolute clusters numbering (for all k) if(ish.ge.6)write(ifch,*)'absolute clusters numbering' if(iLHC.eq.-1)then !fuse touching clusters along k (not used) c midbin=m3grid/2 c if(rangen().gt.0.5)midbin=midbin+1 jjj=1000 k=1 jclu(k)=0 do j=1,m1grid do i=1,m1grid if(jdropgrid(i,j,k).gt.0)then jclu(k)=jclu(k)+1 jjj=jjj+1 ncluj(jjj)=0 kmean(1,jjj)=k kmean(2,jjj)=k jjx=jdropgrid(i,j,k) do jj=j,m1grid do ii=i,m1grid if(jdropgrid(ii,jj,k).eq.jjx)then jdropgrid(ii,jj,k)=jjj ncluj(jjj)=ncluj(jjj)+idropgrid(ii,jj,k) endif enddo enddo endif enddo enddo do k=2,m3grid do j=1,m1grid do i=1,m1grid lnew(i,j)=.true. lold(i,j)=.true. enddo enddo jclu(k)=0 do j=1,m1grid do i=1,m1grid if(lnew(i,j).and.jdropgrid(i,j,k).gt.0)then lcont=.false. if(jdropgrid(i,j,k-1).gt.0)then jj0=jdropgrid(i,j,k-1) nlim=0 if(jj0.gt.1000)nlim=ncluj(jj0) if(nlim.lt.nsegsuj)then jclu(k)=jclu(k)+1 jjx=jdropgrid(i,j,k) if(jjx.gt.jj0)then do jj=jjx,jjj ncluj(jj)=ncluj(jj+1) kmean(1,jj)=kmean(1,jj+1) kmean(2,jj)=kmean(2,jj+1) enddo jjj=jjj-1 jclu(k)=jclu(k)-1 endif kmean(2,jj0)=k do jj=1,m1grid do ii=1,m1grid if(lnew(ii,jj))then if(jdropgrid(ii,jj,k).eq.jjx)then jdropgrid(ii,jj,k)=jj0 ncluj(jj0)=ncluj(jj0)+idropgrid(ii,jj,k) lnew(ii,jj)=.false. elseif(jdropgrid(ii,jj,k).gt.jjx.and.jjx.gt.jj0)then jdropgrid(ii,jj,k)=jdropgrid(ii,jj,k)-1 endif endif enddo enddo else lcont=.true. endif else lcont=.true. endif if(lcont.and.lold(i,j))then jclu(k)=jclu(k)+1 jjj=jjj+1 ncluj(jjj)=0 kmean(1,jjj)=k kmean(2,jjj)=k jjx=jdropgrid(i,j,k) do jj=j,m1grid do ii=i,m1grid if(jdropgrid(ii,jj,k).eq.jjx)then jdropgrid(ii,jj,k)=jjj ncluj(jjj)=ncluj(jjj)+idropgrid(ii,jj,k) lold(ii,jj)=.false. endif enddo enddo endif endif enddo enddo enddo jjs=0 do jj=1001,jjj jjs=jjs+1 kclu(jjs)=(kmean(1,jj)+kmean(2,jj))/2 nseg(jjs)=0 nsegp4(jjs)=0 p4mean(1,jjs)=0d0 p4mean(2,jjs)=0d0 p4mean(3,jjs)=0d0 p4mean(4,jjs)=0d0 rmean(1,jjs)=0d0 !cluster distance to center rmean(2,jjs)=0d0 !cluster maximal radius enddo do k=1,m3grid do j=1,m1grid do i=1,m1grid if(jdropgrid(i,j,k).gt.0)then jdropgrid(i,j,k)=jdropgrid(i,j,k)-1000 endif enddo enddo enddo jjj=jjj-1000 else jjj=jclu(1) do k=2,m3grid do j=1,m1grid do i=1,m1grid if(jdropgrid(i,j,k).gt.0) & jdropgrid(i,j,k)=jdropgrid(i,j,k)+jjj enddo enddo jjj=jjj+jclu(k) enddo jjs=0 do k=1,m3grid do jj=1,jclu(k) jjs=jjs+1 kclu(jjs)=k nseg(jjs)=0 nsegp4(jjs)=0 p4mean(1,jjs)=0d0 p4mean(2,jjs)=0d0 p4mean(3,jjs)=0d0 p4mean(4,jjs)=0d0 rmean(1,jjs)=0d0 !cluster distance to center rmean(2,jjs)=0d0 !cluster maximal radius enddo enddo endif do 20 i=1,maproj do 20 j=1,matarg 20 mpair(j,i)=0 c...calculate mean energy of segments going into in clusters if(ish.ge.6)write(ifch,*) &'calculate mean energy of segments going into in clusters' rmax=0. do 95 n=1,nptla if(iaaptl(n).eq.0)goto 95 i=1+(xptl(n)+xgrid+shiftx)/delxgr j=1+(yptl(n)+xgrid+shiftx)/delxgr k=1+(sptl(n)+sgrid+shifts)/delsgr if( i.ge.1.and.i.le.m1grid & .and.j.ge.1.and.j.le.m1grid & .and.k.ge.1.and.k.le.m3grid)then jj=jdropgrid(i,j,k) if(jj.gt.0)then nsegp4(jj)=nsegp4(jj)+1 io=iorptl(n) if(iLHC.eq.1.and.ityptl(n).lt.40.and.io.gt.0)then !look for corresponding pair of nucleons if(iorptl(io).gt.0)then do while(iorptl(iorptl(io)).ge.0) io=iorptl(io) enddo ip=iorptl(io) it=jorptl(io)-maproj if(ip.gt.0.and.it.gt.0)mpair(it,ip)=mpair(it,ip)+1 endif endif do kk=1,4 p4mean(kk,jj)=p4mean(kk,jj)+dble(pptl(kk,n)) enddo rr=sqrt(xptl(n)**2+yptl(n)**2) rmean(1,jj)=rmean(1,jj)+rr rmax=max(rr,rmax) c if(jj.eq.9)write(ifch,*)'after',n c & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k elseif(istptl(n).eq.3)then !virtual particle is lost (no cluster here) iaaptl(n) = 0 !don't use this particle next time istptl(n) = 5 ior=iorptl(n) ifrptl(1,ior) = 0 ifrptl(2,ior) = 0 do k=1,3 !restore energy to mother particle pptl(k,ior)=pptl(k,ior)+pptl(k,n) enddo pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2 & +pptl(3,ior)**2+pptl(5,ior)**2) call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm) tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen())) endif endif 95 continue ectot=0. amctot=0. maptot=0 mapmax=0 npair=0 do jj=1,jjj ectot=ectot+p4mean(4,jj) amctot=amctot+(p4mean(4,jj)+p4mean(3,jj)) & *(p4mean(4,jj)-p4mean(3,jj))-p4mean(1,jj)**2-p4mean(2,jj)**2 if(nsegp4(jj).ge.nsegce/iflhc)then p4mean(4,jj)=p4mean(4,jj)/dble(nsegp4(jj)) rmean(1,jj)=rmean(1,jj)/dble(nsegp4(jj)) else p4mean(4,jj)=1d50 rmean(1,jj)=-1d0 endif enddo yco=0. ycor=1. if(iLHC.eq.1.and.amctot.gt.0.)then do ip=1,maproj do it=1,matarg if(mpair(it,ip).gt.0)npair=npair+1 mapmax=max(mapmax,mpair(it,ip)) maptot=maptot+mpair(it,ip) enddo enddo amctot=sqrt(amctot) if(amctot.gt.amuseg**2)then visco=1. yrmaxi=yradmx c if(fvisco.gt.0)yrmaxi=yrmaxi*(1.-visco) yrmaxi=yrmaxi*log(exp(yradpi)+amctot**2) c yrmaxi=yrmaxi*amctot**yradpi if(yrmaxi.gt.1e-2)then yyrmax=dble(yrmaxi) fradflii=sngl(1d0/ & ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0))) else yrmaxi=0. fradflii=1. endif else amctot=1. visco=0. fradflii=1. yrmaxi=ainfin !to define ityptl as 19 if mass is too low (aumin=amuseg+yrmaxi in epos-hnb.f) endif if(ylongmx.lt.0.)then delzet=abs(ylongmx)*log(exp(yradmi)+amctot**2) c delzet=abs(ylongmx)*log(exp(yradmi)+amctot) c delzet=abs(ylongmx)*amctot**yradmi yco=delzet !* 1.75 else yco=ylongmx endif ycor=yco if(fvisco.gt.0.)ycor=0. if(ycor.gt.1.e-2)then ycor=sqrt(sinh(ycor)/ycor)/fradflii else ycor=1. endif else if(ylongmx.lt.0)delzet=1.75*delzet endif c...mark segments going into in clusters, count them if(ish.ge.6)write(ifch,*)'mark segments going into in clusters' do 96 n=1,nptla if(iaaptl(n).eq.0)goto 96 i=1+(xptl(n)+xgrid+shiftx)/delxgr j=1+(yptl(n)+xgrid+shiftx)/delxgr k=1+(sptl(n)+sgrid+shifts)/delsgr if( i.ge.1.and.i.le.m1grid & .and.j.ge.1.and.j.le.m1grid & .and.k.ge.1.and.k.le.m3grid)then jj=jdropgrid(i,j,k) if(jj.gt.0)then ! count only particles with reasonnable energy pt2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2 am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2 if(abs(am2tmp).lt.0.1.or. & dble(pptl(4,n)).le.100.d0/dble(ntry)*p4mean(4,jj))then istptl(n) = 3 nseg(jj)=nseg(jj)+1 if(rmean(1,jj).gt.0d0)then rmean(2,jj)=max(rmean(2,jj) & ,abs(rmean(1,jj)-sqrt(xptl(n)**2+yptl(n)**2))) c print *,'r',jj,rmean(1,jj),sqrt(xptl(n)**2+yptl(n)**2),rmean(2,jj) endif c if(jj.eq.9)write(ifch,*)'after',n c & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k else c write(ifmt,*)'Rejected particle : ',n,idptl(n) if(ish.ge.1)write(ifch,*)'Rejected particle : ',n,idptl(n) & ,ityptl(n),(pptl(kk,n),kk=1,5),am2tmp & ,dble(pptl(4,n))/dble(nsegp4(jj)),p4mean(4,jj),nsegp4(jj),jj nsegp4(jj)=nsegp4(jj)-1 iaaptl(n)=0 c if(iaaptl(n).eq.0)print*,'rejected',idptl(n),ityptl(n) c & ,pptl(4,n) if(nsegp4(jj).ge.nsegce/iflhc)then p4mean(4,jj)=(p4mean(4,jj)*dble(nsegp4(jj)+1)-pptl(4,n)) & /dble(nsegp4(jj)) else p4mean(4,jj)=1d50 endif endif endif endif 96 continue c...add segments moving towards clusters if(ish.ge.6)write(ifch,*)'add segments moving towards clusters' if(iocluin.eq.1)then do 93 n=1,nptla if(iaaptl(n).eq.0)goto 93 ihit=0 i=1+(xptl(n)+xgrid+shiftx)/delxgr j=1+(yptl(n)+xgrid+shiftx)/delxgr k=1+(sptl(n)+sgrid+shifts)/delsgr if( i.ge.1.and.i.le.m1grid & .and.j.ge.1.and.j.le.m1grid & .and.k.ge.1.and.k.le.m3grid)then jgr=jdropgrid(i,j,k) if(jgr.eq.0)then if(i.ge.m1grid/2)then isi=-1 else isi=1 endif if(j.ge.m1grid/2)then jsi=-1 else jsi=1 endif do ii=i,i+2*isi,isi do jj=j,j+2*jsi,jsi if(.not.(ii.eq.i.and.jj.eq.j))then if(ii.ge.1.and.ii.le.m1grid . .and.jj.ge.1.and.jj.le.m1grid)then jg=jdropgrid(ii,jj,k) if(jg.gt.0)then c if(nseg(jg).gt.50)then x=xptl(n) y=yptl(n) vrad=( x*pptl(1,n)/pptl(4,n)+y*pptl(2,n)/pptl(4,n)) if(vrad.lt.0.)then ihit=1 goto 94 endif c endif endif endif endif enddo enddo endif endif 94 continue if(ihit.eq.1)then delx=delxgr*(ii-i) dely=delxgr*(jj-j) xn=xptl(n)+delx yn=yptl(n)+dely ix=1+(xn+xgrid+shiftx)/delxgr jx=1+(yn+xgrid+shiftx)/delxgr jgrx=jdropgrid(ix,jx,k) if(jgrx.gt.0)then xptl(n)=xn yptl(n)=yn istptl(n) = 3 nseg(jgrx)=nseg(jgrx)+1 if(rmean(1,jg).gt.0d0) & rmean(2,jgrx)=max(rmean(2,jgrx) & ,abs(rmean(1,jgrx)-sqrt(xptl(n)**2+yptl(n)**2))) endif endif 93 continue endif c Decay Strings not included in clusters taking into account Z c do it here to have this particles before the one from cluster c (to define maxfra properly) if(ish.ge.6)write(ifch,*)'decay Strings not included in clusters' if(ifrade.ne.0.and.ntry.eq.1.and..not.lclean)then nptlx=nptl+1 c copy decayed strings into new strings if no fragments are used for cluster do n=1,nptl if(istptl(n).eq.29.and.ityptl(n).lt.40)then !fragmented central strings iused=0 do nn=ifrptl(1,n),ifrptl(2,n) if(istptl(nn).eq.3.or.ifrptl(1,nn).ne.0)iused=iused+1 enddo if(iused.eq.0)then !if no fragment used, reset string istptl(n)=28 do nn=ifrptl(1,n),ifrptl(2,n) !suppress product istptl(nn)=4 iaaptl(nn)=0 enddo do nn=iorptl(n),jorptl(n) !reinitialize string ends istptl(nn)=20 enddo endif endif enddo call gakfra(0,iret) !fragmentation using Z if(iret.gt.0)goto 1000 do nn=nptlx,nptl !exclude new particle from cluster formation iaaptl(nn)=0 call jtain(nn,xptl(nn),yptl(nn),zptl(nn),tptl(nn),nnn,0) call jtaus(zptl(nn),tz,sptl(nn)) strap=1e10 xpl=tptl(nn)+zptl(nn) xmi=tptl(nn)-zptl(nn) if(xmi.gt.0.0.and.xpl.gt.0.0)strap=0.5*log(xpl/xmi) dezptl(nn)=strap !space-time-rapidity enddo maxfra=nptl nptla=nptl if(ish.ge.6.and.nptl.ne.nptlx-1) & call alist('list after second fragmentation&',nptlx,nptl) if(irescl.eq.1)then call utghost(iret) if(iret.gt.0)goto 1000 endif endif c if(iLHC.eq.1)then cc...split high pt segment into one part going out of the cluster and one part included in cluster (number of particle in cluster conserved) c if(ish.ge.6) c & write(ifch,*)'mark high pt segments going into in clusters' c do 100 n=1,nptla c if(iaaptl(n).ge.0)goto 100 c i=1+(xptl(n)+xgrid+shiftx)/delxgr c j=1+(yptl(n)+xgrid+shiftx)/delxgr c k=1+(sptl(n)+sgrid+shifts)/delsgr c if( i.ge.1.and.i.le.m1grid c & .and.j.ge.1.and.j.le.m1grid c & .and.k.ge.1.and.k.le.m3grid)then c jj=jdropgrid(i,j,k) c if(jj.gt.0)then c if(rmean(1,jj).gt.0d0)then c ! path through cluster c rp=sngl(rmean(2,jj) c & -abs(sqrt(xptl(n)**2+yptl(n)**2)-rmean(1,jj))) c ! energy density (cluster mass approximate as sqrt(E)) c rh=sngl(sqrt(p4mean(4,jj)*nsegp4(jj)) c & /(4d0/3d0*pi*rmean(2,jj)**3)) c ! Et used as Energy (transverse energy) c pt2=pptl(1,n)**2+pptl(2,n)**2 cc ept=sqrt(pt2+pptl(3,n)**2) c ept=sqrt(pt2+pptl(5,n)**2) c ! energy loss based on BDMPS (peigne et al.) c eloss=fploss*rh**(3./8.)*max(1.,sqrt(ept))*rp c ! scaling factor c fscal=eloss/ept cc print *,ept,eloss,nsegp4(jj),rmean(2,jj),rh,rp,fscal c c if(fscal.le.0.)then !exclude fragment from cluster c nseg(jj)=nseg(jj)-1 c istptl(n)=0 c iaaptl(n)=0 c ! update cluster mean energy c nsegp4(jj)=nsegp4(jj)-1 c if(nsegp4(jj).ge.nsegce/iflhc)then c p4mean(4,jj)=(p4mean(4,jj)*dble(nsegp4(jj)+1) c & -pptl(4,n))/dble(nsegp4(jj)) c else c p4mean(4,jj)=1d50 c rmean(1,jj)=-1d0 c endif c elseif(fscal.lt.1..and.(1.-fscal)**2*pt2.gt. c & (0.5-log(rangen()))**2*(ptclu**2))then cc & (0.5-log(rangen()))**2*(pptl(5,n)**2+ptclu**2))then cc print *,ept,eloss,nsegp4(jj),rmean(2,jj),rh,rp,fscal c if(ish.ge.6)write(ifch,*)'Quenched particle : ',n c & ,idptl(n),ityptl(n),(pptl(kk,n),kk=1,5),ept,eloss,nsegp4(jj),jj c iaaptl(n)=0 !mother particle leave cluster c istptl(n)=0 c nptl=nptl+1 !create fake particle with energy lost in cluster c nptla=nptla+1 c istptl(nptl)=3 !daughter part of cluster c iaaptl(nptl)=1 c do l=1,3 c pptl(l,nptl)=fscal*pptl(l,n) c pptl(l,n)=(1.-fscal)*pptl(l,n) c xorptl(l,nptl)=xorptl(l,n) c enddo c! mother c tivptl(1,n)=xorptl(4,n) c pptl(4,n)=sqrt(pptl(1,n)**2+pptl(2,n)**2 c & +pptl(3,n)**2+pptl(5,n)**2) c call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm) c tivptl(2,n)=tivptl(1,n)+taugm*(-alog(rangen())) c ! daughter c pptl(5,nptl)=0. c pptl(4,nptl)=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2 c & +pptl(3,nptl)**2+pptl(5,nptl)**2) c xorptl(4,nptl)=xorptl(4,n) c tivptl(1,nptl)=tivptl(1,n) c iorptl(nptl)=n c idptl(nptl)=idptl(n)*100+sign(99,idptl(n)) c call idquac(nptl,idum1,idum2,idum3,jc) c ityptl(nptl)=ityptl(n)+2 c xptl(nptl)=xptl(n) c yptl(nptl)=yptl(n) c zptl(nptl)=zptl(n) c sptl(nptl)=sptl(n) c dezptl(nptl)=dezptl(n) c if(ish.ge.6)write(ifch,*)'--> Virtual particle : ',nptl c & ,idptl(nptl),ityptl(nptl),(pptl(kk,nptl),kk=1,5) c ! update cluster mean energy c p4mean(4,jj)=(p4mean(4,jj)*dble(nsegp4(jj)) c & -pptl(4,n)+pptl(4,nptl))/dble(nsegp4(jj)) c else c istptl(n)=3 c iaaptl(n)=1 c if(ish.ge.6)write(ifch,*)'Absorbed particle : ',n c & ,idptl(n),ityptl(n),(pptl(kk,n),kk=1,5),ept,eloss,nsegp4(jj),jj c endif c endif c endif c endif c 100 continue c c c c endif c c nptlb0=nptl nptlb=nptl+jjj c...prepare /cptl/ for clusters if(ish.ge.6)write(ifch,*)'prepare /cptl/ for clusters' do jj=1,jjj nptl=nptl+1 istptl(nptl)=12 do l=1,4 pptl(l,nptl)=0. xorptl(l,nptl)=0 enddo sptl(nptl)=0 uptl(nptl)=0 optl(nptl)=0 desptl(nptl)=0 iorptl(nptl)=0 jorptl(nptl)=0 c limit the maximum number of subcluster to half the number of particle c (not to have empty subclusters) nsegmx(jj)=max(1,nint(float(nseg(jj))/float(nsegsuj))) if(ish.ge.6)write(ifch,*)'nsegmx',jj,nseg(jj),nsegmx(jj) & ,nsegsuj enddo c...prepare /cptl/ for subclusters if(ish.ge.6)write(ifch,*)'prepare /cptl/ for subclusters' mm=0 do jj=1,jjj do ii=1,nsegmx(jj) mm=mm+1 if(mm.gt.mxcl)stop'jintpo: mxcl too small. ' mmji(jj,ii)=mm mseg(mm)=0 nptl=nptl+1 istptl(nptl)=10 do l=1,4 pptld(l,mm)=0d0 pptl(l,nptl)=0. xorptl(l,nptl)=0 enddo sptl(nptl)=0 uptl(nptl)=0 optl(nptl)=0 desptl(nptl)=0 do l=1,nflav jccl(mm,l,1)=0 jccl(mm,l,2)=0 enddo iorptl(nptl)=nptla+jj jorptl(nptl)=0 if(ii.eq.1)ifrptl(1,nptla+jj)=nptlb+mm ifrptl(2,nptla+jj)=nptlb+mm enddo enddo c...separate string segments, add dense area segments to clusters if(ish.ge.6)write(ifch,*)'separate string segments' do 98 n=1,nptla if(istptl(n).ne.3)goto 98 i=1+(xptl(n)+xgrid+shiftx)/delxgr j=1+(yptl(n)+xgrid+shiftx)/delxgr k=1+(sptl(n)+sgrid+shifts)/delsgr if( i.ge.1.and.i.le.m1grid & .and.j.ge.1.and.j.le.m1grid & .and.k.ge.1.and.k.le.m3grid)then jj=jdropgrid(i,j,k) if(jj.gt.0)then iimx=nsegmx(jj) do ii=1,iimx mm=mmji(jj,ii) if(mseg(mm).eq.0)goto 10 !not to have an empty cluster enddo ii=1+rangen()*iimx ii=min(ii,iimx) iini=ii am2tmpmx=1e30 9 ntmp=mmji(jj,ii) am2tmp=(pptld(4,ntmp)+pptld(3,ntmp)) & *(pptld(4,ntmp)-pptld(3,ntmp)) if(am2tmp.gt.xmxms**2/2.)then if(am2tmp.lt.am2tmpmx)then mm=mmji(jj,ii) am2tmpmx=am2tmp endif ii=ii+1 if(ii.gt.iimx)ii=1 if(ii.ne.iini)then goto 9 else goto 10 endif endif mm=mmji(jj,ii) 10 continue mseg(mm)=mseg(mm)+1 ifrptl(1,n)=mm !local use of ifrptl c write(ifch,*)'end',n c & ,pptl(4,n),pptl(3,n),idptl(n),i,j,k,sptl(n),mm p4tmp=0d0 do l=1,3 pptld(l,mm)= pptld(l,mm) + dble(pptl(l,n)) p4tmp=p4tmp+dble(pptl(l,n))*dble(pptl(l,n)) enddo p4tmp=sqrt(p4tmp+dble(pptl(5,n)*pptl(5,n))) pptld(4,mm)= pptld(4,mm) + p4tmp c if(mm.eq.86)write(ifch,*)'other',n c & ,pptl(4,n),pptl(3,n),pptl(5,n),idptl(n),k,p4tmp c & ,pptld(4,mm),pptld(3,mm),pptld(2,mm),pptld(1,mm) c & ,(pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm)) xorptl(1,nptlb+mm)=xorptl(1,nptlb+mm)+xptl(n) xorptl(2,nptlb+mm)=xorptl(2,nptlb+mm)+yptl(n) xorptl(3,nptlb+mm)=xorptl(3,nptlb+mm)+zptl(n) xorptl(4,nptlb+mm)=xorptl(4,nptlb+mm)+tptl(n) sptl(nptlb+mm)=sptl(nptlb+mm)+sptl(n) aa=cos(phievt) bb=sin(phievt) cc=-sin(phievt) dd=cos(phievt) xrot=xptl(n)*aa+yptl(n)*bb yrot=xptl(n)*cc+yptl(n)*dd uptl(nptlb+mm)=uptl(nptlb+mm)+xrot**2 optl(nptlb+mm)=optl(nptlb+mm)+yrot**2 desptl(nptlb+mm)=desptl(nptlb+mm)+xrot*yrot call idquac(n,idum1,idum2,idum3,jc) c id=idptl(n) c ida=iabs(id/10) c ids=id/iabs(id) c if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10 c if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids c if(ida.eq.213)id=1230*ids c ic(1)=idtrai(1,id,1) c ic(2)=idtrai(2,id,1) c call iddeco(ic,jc) do l=1,nflav jccl(mm,l,1)=jccl(mm,l,1)+jc(l,1) jccl(mm,l,2)=jccl(mm,l,2)+jc(l,2) enddo else idropgrid(i,j,k)=0 endif endif 98 continue c...associate segments to clusters if(ish.ge.6)write(ifch,*)'associate segments to clusters' naseg(0)=0 do jj=1,jjj do ii=1,nsegmx(jj) mm=mmji(jj,ii) naseg(mm)=naseg(mm-1)+mseg(mm) nfseg(mm)=0 enddo enddo do 97 n=1,nptla if(istptl(n).ne.3)goto 97 istptl(n) = 2 c write(ifch,*)'final segments ',n c & ,ityptl(n),istptl(n),idptl(n),dezptl(n),pptl(3,n) mm=ifrptl(1,n) nfseg(mm)=nfseg(mm)+1 nsegmt(naseg(mm-1)+nfseg(mm))=n 97 continue do jj=1,jjj nst=0 do ii=1,nsegmx(jj) mm=mmji(jj,ii) if(mseg(mm).ne.nfseg(mm))stop'jintpo: mseg.ne.nfseg ' nst=nst+mseg(mm) enddo if(nst.ne.nseg(jj))stop'sum(mseg(mm)).ne.nseg(jj)' enddo c...finish cluster storage to /cptl/ if(ish.ge.6)write(ifch,*)'finish cluster storage to /cptl/' xx=0. yy=0. xy=0. mjjsegsum=0 do jj=1,jjj njj=nptla+jj mjjseg=0 do l=1,nflav jcjj(l,1)=0 jcjj(l,2)=0 enddo do ii=1,nsegmx(jj) mm=mmji(jj,ii) n=nptlb+mm do l=1,nflav jc(l,1)=jccl(mm,l,1) jc(l,2)=jccl(mm,l,2) ke(l)=jc(l,1)-jc(l,2) jcjj(l,1)=jcjj(l,1)+jc(l,1) jcjj(l,2)=jcjj(l,2)+jc(l,2) enddo call idenct(jc,idptl(n) * ,ibptl(1,n),ibptl(2,n),ibptl(3,n),ibptl(4,n)) ttest=0d0 do ji=1,4 ptest(ji)=0d0 do ns=naseg(mm-1)+1,naseg(mm) ni=nsegmt(ns) ptest(ji)=ptest(ji)+dble(pptl(ji,ni)) enddo ptest(ji)=abs(ptest(ji)-pptld(ji,mm)) ttest=ttest+ptest(ji) enddo amcmin=1.01*utamnu(ke(1),ke(2),ke(3),ke(4),ke(5),ke(6),4) p52=((pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm)) & -pptld(1,mm)**2-pptld(2,mm)**2) if(iLHC.eq.1)then amcmin=sqrt(max(amcmin**2,sngl(p52))) if(amcmin/ycor.gt.amuseg)amcmin=amcmin*ycor c & **max(0.,min(1.,pptld(4,mm)-2.*amcmin*ycor)) endif pptld(5,mm)=0d0 if(p52.gt.0d0)then pptld(5,mm)=sqrt(p52) jerr(2)=jerr(2)+1 if(iLHC.eq.1.and.pptld(5,mm).lt.amcmin & .and.pptld(4,mm).gt.amcmin)then pptld(5,mm)=dble(amcmin) bp=sqrt((pptld(4,mm)+pptld(5,mm))*(pptld(4,mm)-pptld(5,mm)) & /(pptld(3,mm)*pptld(3,mm)+pptld(2,mm)*pptld(2,mm) & +pptld(1,mm)*pptld(1,mm))) pptld(1,mm)=bp*pptld(1,mm) pptld(2,mm)=bp*pptld(2,mm) pptld(3,mm)=bp*pptld(3,mm) c write(ifch,*)"ici ",n,sqrt(p52),pptld(5,mm),bp,pptld(4,mm), c & sqrt(pptld(3,mm)*pptld(3,mm) c & +pptld(2,mm)*pptld(2,mm) c & +pptld(1,mm)*pptld(1,mm) c & +pptld(5,mm)*pptld(5,mm)) c write(ifch,*)'droplet uds=',ke(1),ke(2),ke(3),' E=',pptld(5,mm) endif elseif(p52.le.0d0)then jerr(3)=jerr(3)+1 pptld(5,mm)=dble(amcmin) pptld(4,mm)=sqrt(pptld(3,mm)*pptld(3,mm) & +pptld(2,mm)*pptld(2,mm) & +pptld(1,mm)*pptld(1,mm) & +pptld(5,mm)*pptld(5,mm)) endif if(ish.ge.1.and.(abs(ttest).gt.1.d0.or.pptld(5,mm).gt.xmxms)) & then call utmsg('jintpo&') write(ifmt,*)'***** Warning in jintpo !',ntry write(ifch,*)'***** jintpo: momenta messed up (ttest > 0)' write(ifch,*)'*****',mm,n,mseg(mm),p52,ttest write(ifch,*)'*****',jj,nsegp4(jj),p4mean(4,jj) write(ifch,'(a,10x,4f15.4)')'*****',(pptld(ji,mm),ji=1,4) do ns=naseg(mm-1)+1,naseg(mm) ni=nsegmt(ns) write(ifch,'(a,i5,i9,5f15.4,f12.4)')'*****',ni,idptl(ni) * ,(pptl(ji,ni),ji=1,4),pptl(5,ni)**2 * ,(pptl(4,ni)+pptl(3,ni))*(pptl(4,ni)-pptl(3,ni)) * -pptl(1,ni)**2-pptl(2,ni)**2 enddo endif if(pptld(5,mm).gt.xmxms)then p4max=0. nh=0 do ns=naseg(mm-1)+1,naseg(mm) ni=nsegmt(ns) if(pptl(4,ni).ge.p4max)then nh=ni p4max=pptl(4,ni) endif enddo if(nh.le.0)then stop'Cannot be in jintpo ...' else !put back nh as normal particle iaaptl(nh) = 0 !don't use this fragment next time c if(iaaptl(nh).eq.0)print*,'excluded',idptl(nh),ityptl(nh) c & ,pptl(4,nh) if(mod(abs(idptl(nh)),100).eq.99)then !restore lost energy istptl(nh) = 5 ior=iorptl(nh) iaaptl(ior)=0 !don't use this fragment next time ifrptl(1,ior) = 0 ifrptl(2,ior) = 0 if(istptl(ior).eq.0)then do k=1,3 pptl(k,ior)=pptl(k,ior)+pptl(k,n) enddo pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2 & +pptl(3,ior)**2+pptl(5,ior)**2) call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm) tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen())) else istptl(ior) = 0 endif elseif(idptl(nh).lt.1e4)then istptl(nh) = 0 !particle ifrptl(1,nh) = 0 ifrptl(2,nh) = 0 else istptl(nh) = 10 !droplet endif endif if(ish.ge.1) & write(ifch,*)'***** Redo cluster without heavy particle :' & ,nh,ntry goto 8888 endif do l=1,5 pptl(l,n)=sngl(pptld(l,mm)) enddo mjjseg=mjjseg+mseg(mm) do l=1,4 pptl(l,njj)=pptl(l,njj)+pptl(l,n) xorptl(l,njj)=xorptl(l,njj)+xorptl(l,n) xorptl(l,n)=xorptl(l,n)/float(mseg(mm)) enddo sptl(njj)=sptl(njj)+sptl(n) uptl(njj)=uptl(njj)+uptl(n) optl(njj)=optl(njj)+optl(n) desptl(njj)=desptl(njj)+desptl(n) sptl(n)=sptl(n)/float(mseg(mm)) uptl(n)=uptl(n)/float(mseg(mm)) optl(n)=optl(n)/float(mseg(mm)) desptl(n)=desptl(n)/float(mseg(mm)) radptl(n)=0 istptl(n)=10 ifrptl(1,n)=0 ifrptl(2,n)=0 tivptl(1,n)=xorptl(4,n) tivptl(2,n)=ainfin ityptl(n)=60 radptl(n)=0 dezptl(n)=0. enddo !ii do l=1,4 xorptl(l,njj)=xorptl(l,njj)/float(mjjseg) enddo mjjsegsum=mjjsegsum+mjjseg xx=xx+uptl(njj) yy=yy+optl(njj) xy=xy+desptl(njj) sptl(njj)=sptl(njj)/float(mjjseg) uptl(njj)=uptl(njj)/float(mjjseg) optl(njj)=optl(njj)/float(mjjseg) desptl(njj)=desptl(njj)/float(mjjseg) pjj52=(pptl(4,njj)+pptl(3,njj))*(pptl(4,njj)-pptl(3,njj)) & -pptl(1,njj)**2-pptl(2,njj)**2 pptl(5,njj)=0 if(pjj52.gt.0)then pptl(5,njj)=sqrt(pjj52) endif ityptl(njj)=60 call idenct(jc,idptl(njj) * ,ibptl(1,njj),ibptl(2,njj),ibptl(3,njj),ibptl(4,njj)) enddo !jj c...ranphi ranphi=0 rini=0. if(mjjsegsum.gt.0)then xx=xx/float(mjjsegsum) yy=yy/float(mjjsegsum) xy=xy/float(mjjsegsum) dta=0.5*(xx-yy) c eba=0.5*(xx+yy) c ww=-xy c !------------------------------------------------------- c ! inertia tensor: c !----------------------+-------------------------------+ c ! - ! with =uptl ! c ! - ! =optl =desptl ! c !----------------------+-------------------------------+ c ! Eigenvalues ev1, ev2 c !------------------------------------------------------- c ev1=eba+sqrt(dta**2+ww**2) c ev2=eba-sqrt(dta**2+ww**2) if(xy.lt.0..and.dta.ne.0.)then ranphi=0.5*atan(-xy/dta) elseif(xy.gt.0..and.dta.ne.0.)then ranphi=-0.5*atan(xy/dta) else ranphi=0 endif c if(dta.ne.0.)then c ranphi=0.5*atan(abs(ww)/abs(dta)) c if( ww.gt.0..and.dta.gt.0.)then c ranphi=ranphi c elseif(ww.lt.0..and.dta.gt.0.)then c ranphi=-ranphi c elseif(ww.gt.0..and.dta.lt.0.)then c ranphi=pi-ranphi c elseif(ww.lt.0..and.dta.lt.0.)then c ranphi=ranphi-pi c endif c else c ranphi=2.*pi*rangen() c endif rini=max(0.01,sqrt(5./3.*(xx+yy))) !=3/5*R**2 for sphere of radius R volu=(4./3.*pi*(xx+yy)**1.5) c rho=amctot/volu flowpp=0. flowaa=0. if(iLHC.eq.1.and.visco.gt.0.)then if(npair.gt.0)then fcorr=min(1.,float(mapmax)*abs(fvisco)/float(maptot)) visco=min(1.,2.*float(maptot)/float(npair)/float(mapmax))**2 c & *abs(fvisco)) elseif(lclean)then !large number of particles, npair can't be calculated visco=1e-6 fcorr=1. else !cluster from remnants only visco=1. fcorr=1. endif if(visco.ge.1.)yrmaxi=0. !yrmaxi*(1.-visco) c visco=exp(-min(50.,max(0., c & float(koievt)/log(amctot)-abs(fvisco))**yradpi)) !mix flow c & max(0.,rmax**2-abs(fvisco)))) !mix flow c yrmaxi=log(amctot**2) c yrmaxi=yradmx*yrmaxi*(1.-visco) c if(visco.lt.1.and.yrmaxi.gt.1e-2)then c yyrmax=dble(yrmaxi) c fradflii=sngl(1d0/ c & ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0))) c else c visco=1. c yrmaxi=0. c fradflii=1. c endif c if(visco.gt.1e-5)then c yrmaxi=yradmx*(yrmaxi c & +visco*(log(amctot)*yradpx/yradmx-yrmaxi)) c else c yrmaxi=yradmx*yrmaxi c endif fradflii=1. if(yrmaxi.gt.0)then flowpp=visco*log(fcorr*amctot)*yradpx flowaa=yrmaxi if(rangen().lt.flowaa/flowpp)then visco=0. yrmaxi=flowaa+max(0.,flowpp-flowaa) yyrmax=dble(yrmaxi) fradflii=sngl(1d0/ & ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0))) else yrmaxi=0. visco=log(fcorr*amctot)/log(amctot) endif elseif(fcorr.lt.1.)then visco=log(fcorr*amctot)/log(amctot) endif endif if(ish.ge.3)write(ifch,*)'yrmaxi,delzet=',yrmaxi,delzet c print *,'->',bimevt,yrmaxi,visco*yradpx*log(amctot),flowaa c & ,flowpp,maptot,log(ectot),visco,log(amctot),rho c & ,mapmax,npair c & ,min(1.,2.*float(maptot)/float(npair)/float(mapmax))**2 endif c...print if(ish.ge.5)then write(ifch,*)'print' do k=1,m3grid write(ifch,*)'k=',k,' jclu=',jclu(k) & ,' s=',(k-1)*delsgr-sgrid-shifts do j=m1grid,1,-1 write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,m1grid) & ,(jdropgrid(i,j,k),i=1,m1grid) enddo enddo write(ifch,'(a,a)') & ' k jj nseg mm mseg n mass' & ,' s y z t ' do jj=1,jjj do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj)) mm=mmji(jj,ii) n=nptlb+mm sg=pptl(3,n)/abs(pptl(3,n)) tm=sqrt(pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2) y=sg*alog((pptl(4,n)+sg*pptl(3,n))/tm) c if(kclu(jj).eq.44)print *,tm,pptl(4,n),pptl(3,n),iorptl(n) write(ifch,'(2i5,i6,i8,2i6,5f10.3)') & kclu(jj),jj,nseg(jj),mm,mseg(mm),n,pptl(5,n) & ,sptl(n),y,xorptl(3,n),xorptl(4,n) enddo enddo endif c...decay iret=0 if(jjj.gt.0)then !decay only if some cluster produced if(4-abs(typevt).gt.0.0001)typevt=-typevt !typevt < 0 if fusion but only if not SD (sign used for something else for SD ... and no fusion produced for SD events) if(ish.ge.5)write(ifch,*)'decay ...' if(ifrade.eq.0.or.ispherio.gt.0)goto1000 if(jdecay.eq.0)goto1000 nptlbcf=nptl nptl0=nptl if(hydt.ne.'---')then call HydroFO(ier) else nclu=0 ptest(1)=0d0 ptest(2)=0d0 ptest(3)=0d0 ptest(4)=0d0 ptest(5)=0d0 do jj=1,jjj do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj)) mm=mmji(jj,ii) np=nptlb+mm c print *,'decay',jj if(ioclude.eq.3)then call hnbaaa(np,iret) else call DropletDecay(np,iret) endif if(iret.eq.1)then istptl(np)=istptl(np)+2 do ns=naseg(mm-1)+1,naseg(mm) n=nsegmt(ns) if(mod(abs(idptl(n)),100).eq.99)then !restore lost energy istptl(n) = 5 ior=iorptl(n) ifrptl(1,ior) = 0 ifrptl(2,ior) = 0 if(istptl(ior).eq.0)then do k=1,3 pptl(k,ior)=pptl(k,ior)+pptl(k,n) enddo pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2 & +pptl(3,ior)**2+pptl(5,ior)**2) call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm) tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen())) else istptl(ior) = 0 endif elseif(idptl(n).lt.1e4)then istptl(n) = 0 !particle ifrptl(1,n) = 0 ifrptl(2,n) = 0 else istptl(n) = 10 !droplet endif enddo elseif(ioclude.eq.3)then do i=nptl0+1,nptl if(ityptl(i).eq.60)then nclu=nclu+1 ptest(1)=ptest(1)+dble(pptl(1,i)) ptest(2)=ptest(2)+dble(pptl(2,i)) ptest(3)=ptest(3)+dble(pptl(3,i)) ptest(4)=ptest(4)+dble(pptl(4,i)) endif enddo nptl0=nptl endif enddo enddo endif do jj=1,jjj do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj)) mm=mmji(jj,ii) np=nptlb+mm istptl(np)=istptl(np)+1 ifrptl(1,np)=nptlbcf+1 ifrptl(2,np)=nptl rinptl(np)=kclu(jj)-m3grid/2 enddo enddo c add global flow on all particles of all decayed cluster if(iLHC.eq.1)then yrmax=0. if(fvisco.gt.0.)yrmax=yradpx*visco else yrmax=yradpx endif c print *,bimevt,rini,yrmax,yrmaxi,delzet,np,ptmax,visco if(ioclude.eq.3.and.yrmax.gt.1e-3.and.nclu.gt.0)then c set angular informations fecc=0 aa=1 bb=0 cc=0 dd=1 dta=0.5*abs(xx-yy) ev1=(xx+yy)/2+sqrt(dta**2+xy**2) ev2=(xx+yy)/2-sqrt(dta**2+xy**2) ecc=(ev1-ev2)/(ev1+ev2) c fecc=facecc*ecc!/(1.+yrmax) c print*,'pp',ecc,ranphi fecc=min(facecc,ecc) !be careful : fecc change since it is the elliptical deformation of the sub cluster (give strength of v2) phiclu=mod(phievt+ranphi,2.*pi) !do not change otherwise v2 is gone if(phiclu.lt.-pi)phiclu=phiclu+2*pi if(phiclu.gt.pi)phiclu=phiclu-2*pi aa=cos(phiclu) bb=sin(phiclu) cc=-sin(phiclu) dd=cos(phiclu) errlim=0.00005 c loop on particles for each main cluster c ncl=nptlb0 npass=max(1,min(nclu/5,jjj)) !to have the same number of group of particles than original clusters but different repartition of particles npart=nclu/npass if(npart*npass-nclu.gt.max(5,npart/2))npass=npass+1 c print *,'ici',nclu,npass,npart c lcont=.true. c if(nclu.lt.50)lcont=.false. lcont=.false. ncl=0 nmin=nptlbcf+1 nmax=nptl idrc=-1+2.*int(0.5+rangen()) ntot=nclu c prepare debug output for flow nall=0 if(ish.ge.5)then nall=nmax-nmin+1 do ii=1,nall idptl(nptl+ii)=idptl(nptlbcf+ii) enddo iorptl(nptl+nall+1)=nptl+1 jorptl(nptl+nall+1)=nptl+nall do k=1,5 pptl(k,nptl+nall+1)=0 do ii=1,nall pptl(k,nptl+ii)=pptl(k,nptlbcf+ii) pptl(k,nptl+nall+1)=pptl(k,nptl+nall+1)+pptl(k,nptlbcf+ii) enddo enddo endif c initialization for rescaling on ipart particles ipart=0 tecm0=0. tecm=0. ini0=0 ifi0=0 c do 900 while (ncl.le.nptlb-1) do while (ntot.gt.0) ncl=ncl+1 if(iLHC.eq.1)then yrmax=yradpx*visco else yrmax=yradpx endif c cms frame of all particles from same cluster do k=1,5 ppp(k)=0d0 enddo ini=nptl ifi=1 if(idrc.gt.0)then imax=nmax imin=nmin-1 if(ipart.eq.0)ini0=nmin else imax=nmin imin=nmax+1 if(ipart.eq.0)ifi0=nmax endif c npclu=0 c 880 ncl=ncl+1 n=0 i=imin lpass=.true. do while ((n.lt.npart.or.ncl.eq.npass) & .and.idrc*i.lt.idrc*imax.and.lpass) i=i+idrc c if(jorptl(i).eq.ncl)then c if(jorptl(i).eq.ncl.and.ityptl(i).eq.60)then if(ityptl(i).eq.60)then n=n+1 ntot=ntot-1 ini=min(ini,i) ifi=max(ifi,i) c if(ityptl(i).eq.60)npclu=npclu+1 do k=1,4 ppp(k)=ppp(k)+dble(pptl(k,i)) enddo elseif(.not.lcont.and.n.gt.0)then lpass=.false. endif c print *,ityptl(i),i,imin,imax,idrc,n,npart,ncl,nclu,nptl enddo np=n if(idrc.gt.0)then nmin=nmin+i-imin ifi0=i else nmax=nmax+i-imin ini0=i endif c if(ncl.lt.nptlb.and.npclu.lt.int(0.2*(nptl-nptlbcf+1)))goto880 if(ifi.le.ini)goto 900 c record info for rescaling ipart=ipart+np c test mass ppp(5)=(ppp(4)-ppp(3))*(ppp(4)+ppp(3))-ppp(2)**2-ppp(1)**2 if(ppp(5).gt.0d0)then ppp(5)=sqrt(ppp(5)) else if(ish.ge.1)write(ifch,*)'Precision problem in jintpo, p:', & (ppp(k),k=1,5) ppp(5)=0d0 endif if(ish.ge.4) & write(ifch,*)'Group of particle: ', & idrc,ini,ifi,ncl,'/',npass,npart,nclu if(ppp(5).gt.0d0)then ! here all particle should have flow do i=ini,ifi if(ityptl(i).eq.60)then tecm=tecm+pptl(4,i) !use energy in collision frame as reference call utlob4(1,ppp(1),ppp(2),ppp(3),ppp(4),ppp(5) $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) endif enddo if(tecm.gt.0.)then yrmax=yrmax*log(amctot) !/rini**yradpi else yrmax=0. endif c print *,bimevt,rini,ppp(5),yrmax,yrmaxi,delzet,np,ptmax,visco if(ish.ge.4) & write(ifch,*)'Radial flow: ',yrmax,tecm,visco,yradpx,yradpp if(yrmax.gt.0.)then if(np.gt.maxp)stop'maxp too small in jintpo' i=0 if(ish.ge.8)call clist('list before flow&',ini,ifi,60,60) do ii=ini,ifi if(ityptl(ii).eq.60)then i=i+1 yrad(i)=sqrt(rangen()) phirad(i)=2.*pi*rangen() pt2=(pptl(1,ii)**2+pptl(2,ii)**2)!+pptl(5,ii)**2) bex=sinh(yrad(i)*yrmax)*cos(phirad(i)) & *(1+fecc/(1.+pt2)) bey=sinh(yrad(i)*yrmax)*sin(phirad(i)) & *(1-fecc/(1.+pt2)) be(1)=aa*bex+cc*bey be(2)=bb*bex+dd*bey be(3)=0d0 be(4)=sqrt(1+be(1)**2+be(2)**2) c call utlob4(1,be(1),be(2),be(3),be(4),1d0 c * , pptl(1,ii), pptl(2,ii), pptl(3,ii), pptl(4,ii)) c mimic boost transformation but protect against to high values of p(3) (p(3)~p(4)) pt2=pptl(1,ii)**2+pptl(2,ii)**2 bet=-(pptl(1,ii)*be(1)+pptl(2,ii)*be(2))/(be(4)+1.) pt=yradpp**max(1.,pptl(5,ii)) fac=1./(1.+sqrt(pt2/pptl(5,ii)**2))**pt bet=bet+sqrt(pt2+(pptl(3,ii)**2+pptl(5,ii)**2)*fac) c bet=bet+sqrt(pt2+max(yradpp,pptl(5,ii))**2) !simplified version with yradpp~proton mass pptl(1,ii)=pptl(1,ii)-bet*be(1) pptl(2,ii)=pptl(2,ii)-bet*be(2) pptl(4,ii)=sqrt(pptl(1,ii)**2+pptl(2,ii)**2 * +pptl(3,ii)**2+pptl(5,ii)**2) else yrad(i)=0. phirad(i)=0. endif enddo if(ish.ge.8)call clist('list after flow&',ini,ifi,60,60) c boost back pe(1)=0. pe(2)=0. pe(4)=0. do i=ini,ifi if(ityptl(i).eq.60)then call utlob4(-1,ppp(1),ppp(2),ppp(3),ppp(4),ppp(5) $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) pe(1)=pe(1)+pptl(1,i) pe(2)=pe(2)+pptl(2,i) pe(4)=pe(4)+pptl(4,i) endif enddo cc random rescaling c rescaling has to be done for different bins in eta to conserve c energy flow (if everything is done at the end on all particles, c the energy flow is concentred at mid-rapidity in contradiction c with ATLAS data c On the other hand, a sufficient number of particles is necessary c to have proper eta distribution in particular at high pt. c Since the clusters are ordered in rapidity we can do the rescaling on c group of clusters having enough particles for the rescaling but c not too much to keep the eta dependence of the energy flow. c set ipart.ge.1 to do rescaling for each subclusters if((ipart.ge.1.and.ntot.ge.ipart/2).or.ntot.eq.0)then ini=ini0 ifi=ifi0 if(ish.ge.8) & call clist('list before flow rescaling&',ini,ifi,60,60) if(fplmin.gt.0)then c ntrt=ntrt+1 niter=0 611 energ=0. energ0=0. niter=niter+1 i=0 ptmax=0. plmax=0. do j=ini,ifi if(ityptl(j).eq.60)then i=i+1 pt2=pptl(1,j)**2+pptl(2,j)**2 pz2=pptl(3,j)**2 pp2=pt2+pz2 c et2=(pp2+pptl(5,j)**2)*pt2/pp2 pt=sqrt(pt2) c pp=sqrt(pp2) c base necessary to avoid peak at pt or pl=0 c epsi change the shape of eta distributions and pt for c identified particles (shift the maximum of the distribution) c epsi=min(0.99,fplmin*rangen()**0.3) base=0. !sqrt(1./(1.-epsi)**2-1.)*pptl(5,j)/pp finc=2. if(energ0.lt.tecm)then finc=finc*sqrt(float(max(1,niter-300))) c base=base+min(tecm/pe(4), c & (max(0,niter-300))*0.3*pptl(5,j)/pp) else c base=base/log10(max(10.,float(niter-300))) finc=finc/sqrt(max(1.,float(niter-300))) endif c if(1.-pptl(5,j)/pptl(4,j).ge.epsi c & .and.1.+pptl(5,j)**2*(1./pp2-1./pt2).gt.0.)then if(1.+pptl(5,j)**2*(1./pp2-1./pt2).gt.0.)then ptmax=max(pt,ptmax) plmax=max(abs(pptl(3,j)),plmax) yrad2(i)=rangen() yrad2(i)=yrad2(i)* & max(0.,min(1.,finc*tecm/pe(4))-base) c necessary even with epsi cut to smooth eta distribution c and since sumEt has to be conserved, we can constrain scaling for pt c and pz : c ytmp=yrad2(i) yrad2(i)= yrad2(i)**((1.-(pptl(5,j) & /sqrt((min(1.,yrad2(i)*float(max(1,niter-100))))**2 & *(pt**2+pptl(3,j)**2)+pptl(5,j)**2))**1.) & *exp(-fplmin*max(0.,pt**2-(pe(4)/tecm)**2))) yrad2(i)=min(1.,base+yrad2(i)) !should be here to avoid peak at eta=0. else yrad2(i)=1. c ytmp=yrad2(i) endif be(1)= pptl(1,j)*yrad2(i) be(2)= pptl(2,j)*yrad2(i) be(3)= pptl(3,j)*yrad2(i) c print *,niter,i,ntry,yrad2(i),energ/tecm,pt,pptl(3,j),finc,base c * ,ytmp,(1.-(pptl(5,j) c & /sqrt((min(1.,yrad2(i)*float(max(1,niter-100))))**2 c & *(pt**2+pptl(3,j)**2)+pptl(5,j)**2))**0.1) energ=energ+sqrt(be(1)**2+be(2)**2 * +be(3)**2+pptl(5,j)**2) endif enddo energ0=energ c print *,'fin',niter,energ/tecm if(abs(energ-tecm)/tecm.gt.1..and.niter.lt.1000)then goto 611 elseif(niter.ge.1000)then c print *,'Rescaling failed:',pe(4),tecm,ptmax,plmax,ipart,ini0,ifi if(ish.ge.2)write(ifch,*)'Random rescaling failed:' & ,energ,tecm,ptmax,plmax,ipart,ini0,ifi c ntrr=ntrr+1 goto 200 endif c print *,'done',niter,energ/tecm,ipart,ini0,ifi,finc if(ish.ge.5)write(ifch,*)'Rescaling done:' & ,tecm,energ/tecm,niter,ptmax,plmax,ipart,ini0,ifi i=0 do j=ini,ifi if(ityptl(j).eq.60)then i=i+1 pptl(1,j)= pptl(1,j)*yrad2(i) pptl(2,j)= pptl(2,j)*yrad2(i) pptl(3,j)= pptl(3,j)*yrad2(i) pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2 * +pptl(3,j)**2+pptl(5,j)**2) endif enddo endif 200 continue if(ish.ge.8) & call clist('list after flow rescaling&',ini,ifi,60,60) tecm0=tecm0+tecm tecm=0. ipart=0 endif endif ! yrmax endif !p(5)>0 900 continue enddo c rescale momentum precisely and globaly to avoid artefacts for matrix c but only at 10% level to keep dependence of sumEt with eta if(tecm0.gt.0)then ini=nptlbcf+1 ifi=nptl esoll=tecm0 scal=1. do ipass=1,2000 sum=0. n=0 do j=ini,ifi if(ityptl(j).eq.60)then n=n+1 c this part is EXTREMELY important for the pseudorapidity shape at various pt c if nothing special a broad peak appear at eta=0 c to avoid that, scal has to be reduced when p3 or pt reach 0 if(scal.lt.1.)then scal0=scal pt=sqrt(pptl(1,j)**2+pptl(2,j)**2) if(fplmin.le.0.)then pt=abs(fplmin)/sqrt(pptl(5,j))*pt c pt=abs(fplmin)*pt else pt=5.*pt endif pow=sqrt(pptl(5,j)**2+scal**2*(pt**2 * +pptl(3,j)**2)) pow=(1.-(1./sqrt(pptl(5,j))+pptl(5,j)) * /(1./sqrt(pptl(5,j))+pow)) pow=rangen()*pow else pow=1.-2.*pptl(4,j)/engy !to avoid particle with energy larger than beam energy if(pow.lt.0.)then scal0=1./((1.-pow) * *exp(-0.25*max(-4.,log(rangen())))) pow=1. c print *,j,pptl(4,j),pow,scal0,scal c * ,pptl(3,j)*scal0**pow else scal0=scal pow=rangen()*pow endif endif do k=1,3 pptl(k,j)=pptl(k,j)*scal0**pow !to smooth distributions enddo pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2 * +pptl(3,j)**2+pptl(5,j)**2) sum=sum+pptl(4,j) endif enddo scal=esoll/sum c write(ifmt,*)'ipass,scal,e,esoll:' c $ ,ipass,scal,sum,esoll if(abs(scal-1.).le.errlim) goto 300 enddo 300 continue c write(ifmt,*)'ipass,scal,e,esoll:' c $ ,ipass,scal,sum,esoll c adjust pt to have pt conservation in cms of particles having flow if(nclu.gt.0)then ptest(5)=(ptest(4)-ptest(3))*(ptest(4)+ptest(3))-ptest(2)**2 & -ptest(1)**2 if(ptest(5).gt.0d0)then ptest(5)=sqrt(ptest(5)) else if(ish.ge.1)write(ifch,*)'Precision problem in jintpo, p:', & (ptest(k),k=1,5) ptest(5)=0d0 endif be(1)=0.d0 be(2)=0.d0 do i=ini,ifi if(ityptl(i).eq.60)then call utlob4(1,ptest(1),ptest(2),ptest(3),ptest(4),ptest(5) $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) be(1)=be(1)+dble(pptl(1,i)) be(2)=be(2)+dble(pptl(2,i)) endif enddo c shift nclu particles to have sum_pt=0. and boost back in global cms pt1shift=-sngl(be(1)/dble(nclu)) pt2shift=-sngl(be(2)/dble(nclu)) do i=ini,ifi if(ityptl(i).eq.60)then pptl(1,i)=pptl(1,i)+pt1shift pptl(2,i)=pptl(2,i)+pt2shift pptl(4,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2 & +pptl(3,i)**2+pptl(5,i)**2) call utlob4(-1,ptest(1),ptest(2),ptest(3),ptest(4),ptest(5) & ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) endif enddo endif endif c if(ntrt.gt.0)print *,"jintpo rescaling :",float(ntrr)/float(ntrt) c define life time n=0 do i=ini,ifi if(ityptl(i).eq.60)then n=n+1 r=1.15*rini*yrad(n) !yrad=y/ymax tau=2.25/sqrt(yrad(n)**2+0.04)-0.75 z=xorptl(3,i) t=xorptl(4,i) ! zeta=0.5*log((t+z)/(t-z))-0.5*delzet+2*0.5*delzet*rangen() test=(pptl(4,i)-pptl(3,i))*(pptl(4,i)+pptl(3,i)) if(test.gt.0.)then zeta=0.5*log((pptl(4,i)+pptl(3,i)) & /(pptl(4,i)-pptl(3,i))) else !in case of precision problem (but not always good neither for p<0 pt=sqrt(pptl(2,i)**2+pptl(1,i)**2) zeta=0.5*log(1+2*pptl(3,i)*(pptl(4,i)+pptl(3,i)) & /(pt*pt+pptl(5,i)**2)) endif z=tau*sinh(zeta) t=tau*cosh(zeta) xorptl(1,i)=xorptl(1,i)+r*cos(phirad(n)) xorptl(2,i)=xorptl(2,i)+r*sin(phirad(n)) xorptl(3,i)=z xorptl(4,i)=t endif enddo if(ish.ge.5)then do k=1,5 pptl(k,nptl+nall+2)=0 do ii=nptlbcf+1,nptl pptl(k,nptl+nall+2)=pptl(k,nptl+nall+2)+pptl(k,ii) enddo enddo iorptl(nptl+nall+2)=nptlbcf+1 jorptl(nptl+nall+2)=nptl call alist2('longitudinal and radial flow&',nptl+1 & ,nptl+nall,nptlbcf+1,nptl) call alist2('momentum sum&',nptl+nall+1,nptl+nall+1 & ,nptl+nall+2,nptl+nall+2) write(ifch,'(1x,50a1/)')('-',k=1,50) endif endif !ioclude=3 and flow do n=nptlbcf+1,nptl if(ioclude.ne.3)then iorptl(n)=nptlb+1 jorptl(n)=nptlbcf rinptl(n)=rinptl((iorptl(n)+jorptl(n))/2) else rinptl(n)=rinptl(iorptl(n)) endif istptl(n)=0 ifrptl(1,n)=0 ifrptl(2,n)=0 tivptl(1,n)=xorptl(4,n) call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm) r=rangen() tivptl(2,n)=tivptl(1,n)+taugm*(-alog(r)) radptl(n)=0. dezptl(n)=0. itsptl(n)=0 enddo endif c Decay droplets not included in clusters iret=0 do mm=1,nptla nptlb=nptl if(istptl(mm).eq.10)then if(ish.ge.5)write(ifch,*)'Decay remaining droplet :',mm if(nptlb.gt.mxptl-10) & call utstop('jintpo: mxptl too small (2)&') if(ioclude.eq.3)then call hnbaaa(mm,iret) else call DropletDecay(mm,iret) !Decay remn iret=0 endif if(iret.eq.0.and.nptl.ne.nptlb)then ! ---successful decay--- istptl(mm)=istptl(mm)+1 ifrptl(1,mm)=nptlb+1 ifrptl(2,mm)=nptl t=tivptl(2,mm) x=xorptl(1,mm)+(t-xorptl(4,mm))*pptl(1,mm)/pptl(4,mm) y=xorptl(2,mm)+(t-xorptl(4,mm))*pptl(2,mm)/pptl(4,mm) z=xorptl(3,mm)+(t-xorptl(4,mm))*pptl(3,mm)/pptl(4,mm) do 21 n=nptlb+1,nptl iorptl(n)=mm jorptl(n)=0 istptl(n)=0 ifrptl(1,n)=0 ifrptl(2,n)=0 radius=0.8*sqrt(rangen()) phi=2*pi*rangen() ti=t zi=z xorptl(1,n)=x + radius*cos(phi) xorptl(2,n)=y + radius*sin(phi) xorptl(3,n)=zi xorptl(4,n)=ti iioo=mm zor=dble(xorptl(3,iioo)) tor=dble(xorptl(4,iioo)) r=rangen() tauran=-taurea*alog(r) call jtaix(n,tauran,zor,tor,zis,tis) tivptl(1,n)=amax1(ti,tis) call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm) r=rangen() tivptl(2,n)=t+taugm*(-alog(r)) radptl(n)=0. dezptl(n)=0. itsptl(n)=0 rinptl(nptl)=-9999 21 continue else ! Unsuccessful decay if(ish.ge.1)write(ifch,*) * '***** Unsuccessful remnant cluster decay' * ,' --> redo event.' endif endif enddo 1000 continue call utprix('jintpo',ish,ishini,4) end cc----------------------------------------------------------------------- c subroutine jrad(i,nq,na,jc,rad) cc----------------------------------------------------------------------- cc return hadron radius (data taken from huefner and povh) cc----------------------------------------------------------------------- c include 'epos.inc' c integer jc(nflav,2),kc(nflav) c c id=iabs(idptl(i)) c am=pptl(5,i) c if(id.lt.10000)then c k=mod(id,10) c else c k=1 c endif c do l=1,nflav c kc(l)=iabs(jc(l,1)-jc(l,2)) c enddo c c if(nq.eq.0)then ! mesons c if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0.and.kc(4).eq.0)then c if(k.eq.0)then ! flavor singlet pseudoscalar mesons c if(am.ge.0.000)then c rad=0.64 ! pi0 c if(am.ge.0.500)then c rad=0.60 ! eta c if(am.ge.0.900)then c rad=0.40 ! eta prime c if(am.ge.2.900)then c rad=0.17 ! eta charm c endif c endif c endif c else c write(ifch,*) c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6) c call utstop('jrad: meson radius not defined&') c endif c else ! flavor singlet vector mesons c if(am.ge.0.000)then c rad=0.72 ! rho,omega c if(am.ge.1.000)then c rad=0.46 ! phi c if(am.ge.3.000)then c rad=0.20 ! J/psi c endif c endif c else c write(ifch,*) c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6) c call utstop('jrad: meson radius not defined&') c endif c endif c elseif(kc(3).eq.0.and.kc(4).eq.0)then ! nonstrange, noncharmed c if(k.eq.0)then c rad=0.64 ! pi c else c rad=0.72 ! resonances c endif c elseif(kc(3).ne.0.and.kc(4).eq.0)then ! strange c if(k.eq.0)then c rad=0.59 ! kaons c else c rad=0.68 ! kaon resonances c endif c else ! charmed c write(ifch,*)'i:',i,' id:',idptl(i) c call utstop('jrad: radius of meson not defined&') c endif c else !baryons c if(kc(4).gt.0)then ! charmed c write(ifch,*) c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am c write(ifch,*)'i:',i,' id:',idptl(i) c call utstop('jrad: radius of charmed baryon not defined&') c elseif(kc(3).eq.0)then ! nonstrange c if(k.eq.0)then c rad=0.82 !nucleons c else c rad=1.00 !resonances c endif c elseif(kc(3).eq.1)then ! strange c if(k.eq.0)then c rad=0.76 !lambda, sigma c else c rad=0.93 !resonances c endif c elseif(kc(3).eq.2)then ! double strange c if(k.eq.0)then c rad=0.71 !cascades c else c rad=0.87 !resonances c endif c elseif(kc(3).ge.3)then ! triple strange c rad=0.79 !omega c else c write(ifch,*) c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am c write(ifch,*) c * 'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6), c * ' |q-qbar|:',(kc(l),l=1,6) c call utstop('jrad: should not happen&') c endif cc string fragments with |#q|>3 c if(na.gt.3)then c a=(na/3.)**(1./3.) c if(ish.ge.7)then c call utmsg('jrad ') c write(ifch,*) c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am c write(ifch,*) c * 'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6), c * ' |q-qbar|:',(kc(l),l=1,6) c write(ifch,*)'nq:',nq,' na:',na,' r:',rad,' ar:',a*rad c call utmsgf c endif c rad=rad*a c endif c endif c c if(ish.ge.7)then c write(ifch,*) c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am,' rad:',rad c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6) c endif c c return c end c c----------------------------------------------------------------------- subroutine jresc c----------------------------------------------------------------------- include 'epos.inc' double precision pa(5),pj(5) integer ipptl(mxptl) call utpri('jresc ',ish,ishini,4) iret=0 nptlpt=maproj+matarg np=0 do i=nptlpt+1,nptl if(istptl(i).eq.0 * .and.idptl(i).lt.10000.and.pptl(5,i).gt.0.01)then np=np+1 ipptl(np)=i endif enddo if(np.lt.2)goto1001 do ii=1,np i=ipptl(ii) if(mod(iabs(idptl(i)),10).lt.2)then call idmass(idptl(i),ami) dm=abs(ami-pptl(5,i)) if(dm.gt.0.001)then ntry=0 1 continue ntry=ntry+1 2 jj=1+int(rangen()*np) j=ipptl(jj) if(ish.ge.4)write(ifch,*)i,j,istptl(j) if(j.eq.i)goto2 if(mod(iabs(idptl(j)),10).lt.2)then call idmass(idptl(j),amj) else amj=pptl(5,j) endif do l=1,5 pa(l)=dble(pptl(l,i)) pj(l)=dble(pptl(l,j)) enddo if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70) if(ish.ge.4)write(ifch,11)i,idptl(i),'before:',pa,'want:',ami if(ish.ge.4)write(ifch,11)j,idptl(j),'before:',pj,'want:',amj call jrescl(pa,dble(ami),pj,dble(amj),iret) if(iret.eq.1)then if(ntry.le.50)then goto1 else goto1001 endif endif if(ish.ge.4)write(ifch,11)i,idptl(i),' after:',pa if(ish.ge.4)write(ifch,11)j,idptl(j),' after:',pj if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70) do l=1,5 pptl(l,i)=sngl(pa(l)) pptl(l,j)=sngl(pj(l)) enddo endif endif enddo 11 format(i5,1x,i5,1x,a,1x,5(d8.2,1x),a,1x,e8.2) 1000 continue call utprix('jresc ',ish,ishini,4) return 1001 continue if(ish.ge.1)then write(ifmt,'(a)')'jresc: could not put on shell' endif goto1000 end c----------------------------------------------------------------------- subroutine jrescl(p1,am1,p2,am2,iret) c----------------------------------------------------------------------- c rescale momenta of two particles such that the masses assume given c values. c input: c p1, p2: momenta of the two particles c am1, am2: desired masses of the two particles c output: c p1, p2: new momenta of the two particles c----------------------------------------------------------------------- include 'epos.inc' double precision p1(5),p2(5) * ,p1n(5),p2n(5) * ,a1,a2,a12,am1,am2 * ,b1,b2,c,d,e,f,g,p,q,r call utpri('jrescl',ish,ishini,7) iret=0 a1=p1(5)**2 a2=p2(5)**2 a12=p1(4)*p2(4)-p1(3)*p2(3)-p1(2)*p2(2)-p1(1)*p2(1) if(a12.le.(a1+a2))then if(ish.ge.7)write(ifch,*)'a_12 < a_1 + a_2' if(ish.ge.7)write(ifch,*)a12,' < ',a1+a2 c goto1001 endif 11 format(5(d9.3,1x)) if(ish.ge.7)write(ifch,11)p1,a1 if(ish.ge.7)write(ifch,11)p2,a2 if(ish.ge.7)write(ifch,*)a12 c=(a1+a12)/(a2+a12) d=(a1-am1**2-a2+am2**2)/(a2+a12)*0.5d0 e=a1-2d0*a12*c+a2*c**2 f=2d0*(a1-a12*(c+d)+a2*c*d) g=a1-2d0*a12*d+a2*d**2-am1**2 p=f/e q=g/e r=p**2-4d0*q if(ish.ge.7)write(ifch,*)'c:',c,' d:',d if(ish.ge.7)write(ifch,*)'e:',e,' f:',f,' g:',g if(ish.ge.7)write(ifch,*)'p:',p,' q:',q,' r:',r if(r.lt.0d0)goto1001 b1=-0.5d0*(p-dsqrt(r)) b2=b1*c+d if(ish.ge.7)write(ifch,*)'b_1:',b1,' b_2:',b2 do i=1,4 p1n(i)=(1d0+b1)*p1(i)-b2*p2(i) p2n(i)=(1d0+b2)*p2(i)-b1*p1(i) enddo a1=p1n(4)**2-p1n(3)**2-p1n(2)**2-p1n(1)**2 a2=p2n(4)**2-p2n(3)**2-p2n(2)**2-p2n(1)**2 if(a1.gt.0d0.and.a2.gt.0d0)then do i=1,4 p1(i)=p1n(i) p2(i)=p2n(i) enddo p1(5)=dsqrt(a1) p2(5)=dsqrt(a2) if(ish.ge.7)write(ifch,11)p1,a1 if(ish.ge.7)write(ifch,11)p2,a2 else goto1001 endif if(p1(4).lt.0..or.p2(4).lt.0.)goto1001 1000 continue call utprix('jrescl',ish,ishini,7) return 1001 continue iret=1 goto1000 end c----------------------------------------------------------------------- subroutine jtain(i,x,y,z,t,n,iopt) c----------------------------------------------------------------------- c returns intersection (x,y,z,t) of ptl-i-trajectory with taus-line. c input: c i: particle number c iopt: formation time considered (0) or not (1) c output: c x,y,z,t: 4-vector of intersection point c n: exit code c n=0: ok c n=1: ptl lives later c n=2: ptl lives earlier c n=9: tiv1>tiv2 c----------------------------------------------------------------------- include 'epos.inc' double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat double precision vv,zza,zz,tt,xo3,xo4,ti1,ti2,derr,dd double precision ttp,zzp,ttt,zzt,vvt,vvp,spt2m2E,p4 common/ctfi/tin,tfi double precision ttau0 common/cttau0/ttau0 n=0 tin=0 tfi=0 derr=1d-2 ttp=tpro*ttaus zzp=zpro*ttaus ttt=ttar*ttaus zzt=ztar*ttaus vv=sign(min(1.d0,abs(dble(pptl(3,i)))/dble(pptl(4,i))) & ,dble(pptl(3,i))) if(abs(vv).ge.1.d0)then spt2m2E=dble(pptl(1,i)*pptl(1,i)+pptl(2,i)*pptl(2,i) & +pptl(5,i)*pptl(5,i)) c if(pptl(4,i).le.0.)then p4=sqrt(dble(pptl(3,i)*pptl(3,i))+spt2m2E) c else c p4=dble(pptl(4,i)) c endif ctp to avoid precision problem, replace abs(p3)/p4 by sqrt(1-(pt2+m2)/E2) spt2m2E=min(1.d0,sqrt(spt2m2E)/p4) vv=sign(sqrt((1d0+spt2m2E)*(1d0-spt2m2E)),dble(pptl(3,i))) endif xo3=dble(xorptl(3,i)) xo4=dble(xorptl(4,i)) zza=xo3-xo4*vv if(iopt.eq.0)then ti1=dble(tivptl(1,i)) elseif(iopt.eq.1)then ti1=dble(xo4) else ti1=0 call utstop("Wrong iopt in jtain !&") endif ti2=dble(tivptl(2,i)) if(ti1.gt.ti2)then n=9 goto1 endif zfi=sngl(xo3+(ti2-xo4)*vv) call jtaus(zfi,tzfi,szfi) tfi=tzfi if(tfi.ge.sngl(ti2))then n=2 goto1 endif zin=sngl(xo3+(ti1-xo4)*vv) call jtaus(zin,tzin,szin) tin=tzin if(tin.le.sngl(ti1))then n=1 goto1 endif 1 continue if(ttaus.le.ttau0)then tt=ttaus zz=xo3+(tt-xo4)*vv if(tt.lt.ti1.and.n.eq.0)n=1 if(tt.ge.ti2.and.n.eq.0)n=2 goto1000 else vvt=zzt/ttt vvp=zzp/ttp tt=(ttt+(zza-zzt)*vvt)/(1-vv*vvt) zz=xo3+(tt-xo4)*vv if(zz.le.zzt)then if(tt.lt.ti1.and.n.eq.0)n=1 if(tt.ge.ti2.and.n.eq.0)n=2 goto1000 endif tt=(ttp+(zza-zzp)*vvp)/(1-vv*vvp) zz=xo3+(tt-xo4)*vv if(zz.ge.zzp)then if(tt.lt.ti1.and.n.eq.0)n=1 if(tt.ge.ti2.and.n.eq.0)n=2 goto1000 endif dd=1-vv**2 if(sngl(dd).eq.0..and.vv.gt.0.)then tt=-(ttaus**2+zza**2)/2d0/zza elseif(sngl(dd).eq.0..and.vv.lt.0.)then tt=(ttaus**2+zza**2)/2d0/zza else tt=(zza*vv+dsqrt(zza**2+ttaus**2*dd))/dd endif zz=xo3+(tt-xo4)*vv if(tt.lt.ti1.and.n.eq.0)n=1 if(tt.ge.ti2.and.n.eq.0)n=2 if(dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr*ttaus**2.and. *dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr)then if(ish.ge.1)then call utmsg('jtain') write(ifch,*)'***** ttaus**2 .ne. (tt+zz)*(tt-zz)' write(ifch,*)sngl(ttaus**2),sngl((tt+zz)*(tt-zz)) call utmsgf endif goto1000 endif endif 1000 t=sngl(tt) z=sngl(zz) x=xorptl(1,i)+(t-xorptl(4,i))*pptl(1,i)/pptl(4,i) y=xorptl(2,i)+(t-xorptl(4,i))*pptl(2,i)/pptl(4,i) return end c----------------------------------------------------------------------- subroutine jtaix(i,tau,zor,tor,z,t) c----------------------------------------------------------------------- c returns intersection z,t of ptl-i-trajectory with hyperbola h. c h: (t-tor)**2-(z-zor)**2=tau**2 . c zor, tor double precision. c----------------------------------------------------------------------- include 'epos.inc' double precision tor,zor,tors,zors,vv,cc,dd,ttau,derr,tt,zz derr=1d-3 ttau=dble(tau) zors=dble(xorptl(3,i))-zor tors=dble(xorptl(4,i))-tor vv=dble(pptl(3,i)/pptl(4,i)) vv=dmin1(vv,1d0) vv=dmax1(vv,-1d0) cc=zors-tors*vv dd=1d0-vv**2 dd=dmax1(dd,0d0) if(dd.eq.0d0.and.cc.eq.0d0)then if(tau.eq.0.)tt=0d0 if(tau.ne.0.)tt=dble(ainfin) zz=tt goto1000 elseif(dd.eq.0d0)then tt=-(ttau**2+cc**2)/2d0/cc/vv elseif(dd.lt.1e-8)then tt=-(ttau**2+cc**2)/2d0/cc/vv call utmsg('jtaix') write(ifch,*)'***** dd = ',dd,' treated as zero' call utmsgf else tt=(cc*vv+dsqrt(cc**2+ttau**2*dd)) tt=tt/dd endif zz=cc+tt*vv if(dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr*ttau**2 *.and.dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr *.and.tors**2-zors**2.lt.1e6)then if(ish.ge.2)then call utmsg('jtaix') write(ifch,*)'***** ttau**2 .ne. (tt+zz)*(tt-zz)' write(ifch,*)sngl(ttau**2),sngl((tt+zz)*(tt-zz)) write(ifch,*)'tau,t,z:' write(ifch,*)tau,tt,zz write(ifch,*)'#,id(ptl):',i,idptl(i) write(ifch,*)'zor,tor(str):',zor,tor write(ifch,*)'zors,tors,p,e(ptl):' write(ifch,*)sngl(zors),sngl(tors),pptl(3,i),pptl(4,i) call utmsgf endif endif 1000 z=sngl(zz+zor) t=sngl(tt+tor) return end c----------------------------------------------------------------------- subroutine jtaug(su,so,g,y) c----------------------------------------------------------------------- c returns g factor and rapidity y for given su, so c----------------------------------------------------------------------- include 'epos.inc' double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat double precision ttp,zzp,ttt,zzt,ssp,sst,ssu,sso,ss1,ss2,gg *,ssav,yyav,hh double precision ttau0 common/cttau0/ttau0 ssu=dble(su) sso=dble(so) if(ssu.ge.sso)then sso=(ssu+sso)*0.5d0 + dble(abs(dezzer))*ttaus*0.5d0 ssu=(ssu+sso)*0.5d0 - dble(abs(dezzer))*ttaus*0.5d0 so=real(sso) su=real(ssu) endif if(ssu.ge.sso)then print*,ssu,sso,dble(abs(dezzer))*ttaus*0.5d0 stop'STOP: sr jtaug: ssu.ge.sso' endif g=1 if(ttaus.le.ttau0)return ttp=tpro*ttaus zzp=zpro*ttaus ttt=ttar*ttaus zzt=ztar*ttaus ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp)) sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt)) ssav=(ssu+sso)/2d0 yyav=ssav/ttaus if(ssav.ge.ssp)yyav=detap if(ssav.le.sst)yyav=detat gg=0 if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu) if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu)) if(ssu.lt.ssp.and.sso.gt.sst)then ss1=dmax1(ssu,sst) ss2=dmin1(sso,ssp) gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) ) endif gg=gg/(sso-ssu) hh=0 if(ssu.lt.sst)hh=hh + dsinh(detat-yyav) * (dmin1(sst,sso)-ssu) if(sso.gt.ssp)hh=hh + dsinh(detap-yyav) * (sso-dmax1(ssp,ssu)) if(ssu.lt.ssp.and.sso.gt.sst)then ss1=dmax1(ssu,sst) ss2=dmin1(sso,ssp) hh=hh+ttaus*( dcosh(ss2/ttaus-yyav)-dcosh(ss1/ttaus-yyav) ) endif hh=hh/(sso-ssu) yyav=yyav+0.5d0*dlog((gg+hh)/(gg-hh)) gg=0 if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu) if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu)) if(ssu.lt.ssp.and.sso.gt.sst)then ss1=dmax1(ssu,sst) ss2=dmin1(sso,ssp) gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) ) endif gg=gg/(sso-ssu) g=sngl(gg) y=sngl(yyav) return end c----------------------------------------------------------------------- subroutine jtaui(s,ts,zs) c----------------------------------------------------------------------- c returns time ts and coord zs corresponding to ttaus and inv. length s c----------------------------------------------------------------------- double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat double precision ttau0 common/cttau0/ttau0 double precision ttp,zzp,ttt,zzt,ssp,sst,ss,zeta zs=s ts=sngl(ttaus) if(ttaus.le.ttau0)return ttp=tpro*ttaus zzp=zpro*ttaus ttt=ttar*ttaus zzt=ztar*ttaus ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp)) sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt)) ss=dble(s) if(ss.le.sst)then zs=sngl(zzt+ttar*(ss-sst)) ts=sngl(ttt+(dble(zs)-zzt)*zzt/ttt) elseif(ss.ge.ssp)then zs=sngl(zzp+tpro*(ss-ssp)) ts=sngl(ttp+(dble(zs)-zzp)*zzp/ttp) else zeta=ss/ttaus ts=sngl(ttaus*dcosh(zeta)) zs=sngl(ttaus*dsinh(zeta)) endif return end c----------------------------------------------------------------------- subroutine jtauin c----------------------------------------------------------------------- c initializes equal time hyperbola at ttaus c----------------------------------------------------------------------- include 'epos.inc' double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat double precision ttau0,rcproj,rctarg common/geom1/rcproj,rctarg common/cttau0/ttau0 call utpri('jtauin',ish,ishini,6) if(ttaus.gt.ttau0)then if(rcproj.gt.1d-10)then detap=dmin1(dble((ypjtl-yhaha)*etafac),dlog(ttaus/rcproj)) else detap=dble((ypjtl-yhaha)*etafac) endif if(rctarg.gt.1d-10)then detat=dmax1(dble(-yhaha*etafac),dlog(rctarg/ttaus)) else detat=dble(-yhaha*etafac) endif tpro=dcosh(detap) zpro=dsinh(detap) ttar=dcosh(detat) ztar=dsinh(detat) else detap=0d0 detat=0d0 tpro=0d0 zpro=0d0 ttar=0d0 ztar=0d0 endif if(ish.ge.6)then write(ifch,*)'hyperbola at tau=',ttaus write(ifch,*)'r_p:',rcproj,' r_t:',rctarg write(ifch,*)'y_p:',detap,' y_t:',detat write(ifch,*)'t_p:',tpro,' z_p:',zpro write(ifch,*)'t_t:',ttar,' z_t:',ztar endif call utprix('jtauin',ish,ishini,6) return end c----------------------------------------------------------------------- subroutine jtaus(z,tz,sz) c----------------------------------------------------------------------- c returns time tz and inv length sz corresponding to ttaus and z c----------------------------------------------------------------------- include 'epos.inc' double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat double precision ttau0 common/cttau0/ttau0 double precision ttp,zzp,ttt,zzt,zz,tzz tz=sngl(ttaus) sz=z if(ttaus.le.ttau0)return ttp=tpro*ttaus zzp=zpro*ttaus ttt=ttar*ttaus zzt=ztar*ttaus zz=dble(z) if(zz.le.zzt)then tz=sngl(ttt+(zz-zzt)*zzt/ttt) sz=sngl(ttaus*detat+(zz-zzt)/ttar) elseif(zz.ge.zzp)then tz=sngl(ttp+(zz-zzp)*zzp/ttp) sz=sngl(ttaus*detap+(zz-zzp)/tpro) else if(sngl(ttaus).ge.ainfin)then tz=sngl(ttaus) sz=0. if(ish.ge.1)then call utmsg('jtaus') write(ifch,*)'***** large ttaus; set tz=ttaus, sz=0' write(ifch,*)'ttaus=',ttaus,'zz=',zz call utmsgf endif else tzz=dsqrt(ttaus**2+zz**2) tz=sngl(tzz) sz=sngl(ttaus*0.5d0*dlog((tzz+zz)/(tzz-zz))) endif endif return end c----------------------------------------------------------------------- subroutine jtaux(t,z,ttaux) c----------------------------------------------------------------------- c returns ttaux (-> tau-line) for given t and z. c ttaux: double precision. c----------------------------------------------------------------------- double precision ttaux,tt,zz,rcproj,rctarg,zt1,zp1,zt2,zp2,ttau0 common/geom1/rcproj,rctarg common/cttau0/ttau0 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat tt=dble(t) zz=dble(z) if(tt.gt.ttau0)then zt1=rctarg-tt zp1=tt-rcproj zt2=ztar/ttar*tt zp2=zpro/tpro*tt if(zz.le.dmax1(zt1,zt2))then if(zt1.gt.zt2)then ttaux=rctarg*dsqrt((tt-zz)/(2d0*rctarg-tt-zz)) else ttaux=(ttar*tt-ztar*zz)/(ttar**2-ztar**2) endif elseif(zz.ge.dmin1(zp1,zp2))then if(zp1.lt.zp2)then ttaux=rcproj*dsqrt((tt+zz)/(2d0*rcproj-tt+zz)) else ttaux=(tpro*tt-zpro*zz)/(tpro**2-zpro**2) endif else ttaux=dsqrt(tt**2-zz**2) endif else ttaux=tt endif return end c----------------------------------------------------------------------- subroutine xjden1(ii,itau,x,y,rad,o,u) c----------------------------------------------------------------------- c ii=0: initialization c ii=1: determining dense regions in space of individual events c x,y,rad: tranverse coordinates and radius of particle i c o,u: specifies long range: u < s < o (s: long coordinate) c ii=2: plot of individual event c xpar1: itau ; valid: 1,..,10 c xpar2: z-range: -xpar2 < z < xpar2 c xpar3, x-range: -xpar3 < x < xpar3 c xpar4, y-range: -xpar4 < y < xpar4 c----------------------------------------------------------------------- include "epos.inc" double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat if(idensi.ne.1)stop'STOP in xjden1: idensi must be set 1' dlcoox=0.5 dlcooy=0.5 if(ii.eq.0)then do i=1,nzeta do j=1,mxcoox do k=1,mxcooy kdensi(itau,i,j,k)=0 enddo enddo enddo elseif(ii.eq.1)then if(itau.lt.1.or.itau.gt.mxtau)return tau=sngl(ttaus) zu=u/tau zo=o/tau do 1 i=1,nzeta zi=-nzeta*dlzeta/2-dlzeta/2+i*dlzeta if(zu.gt.zi.or.zo.lt.zi)goto1 do 2 j=1,mxcoox xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox do 3 k=1,mxcooy yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy if((x-xj)**2+(y-yk)**2.gt.rad**2)goto3 kdensi(itau,i,j,k)=1 3 continue 2 continue 1 continue elseif(ii.eq.2)then itaux=nint(xpar1) if(itaux.gt.mxtau)stop'STOP in xjden1: itaux too large' iz=nint(xpar2/dlzeta) ix=nint(xpar3/dlcoox) iy=nint(xpar4/dlcooy) if(iz.gt.nzeta/2)stop'STOP in xjden1: zeta-range too large' if(ix.gt.mxcoox/2)stop'STOP in xjden1: x-range too large' if(iy.gt.mxcooy/2)stop'STOP in xjden1: y-range too large' do k=mxcooy/2+1-iy,mxcooy/2+iy write(ifhi,'(a,f7.2)') '! tau: ',tauv(itaux) write(ifhi,'(a)') 'openhisto' write(ifhi,'(a,2f7.2)')'xrange',-iz*dlzeta,iz*dlzeta write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox write(ifhi,'(a)') 'set ityp2d 3' write(ifhi,'(a)') 'txt "xaxis space-time rapidity [z]"' write(ifhi,'(a)') 'txt "yaxis transverse coordinate x (fm)"' write(ifhi,'(a,i4)') 'array2d',2*iz do j=mxcoox/2+1-ix,mxcoox/2+ix write(ifhi,'(40i2)') (kdensi(itaux,i,j,k), * i=nzeta/2+1-iz,nzeta/2+iz) enddo write(ifhi,'(a)') ' endarray' write(ifhi,'(a)') 'closehisto plot2d' enddo else stop'STOP in xjden1: wrong option' endif return end c----------------------------------------------------------------------- subroutine xjden2(ii,itau,x,y,rad,s) c----------------------------------------------------------------------- c ii=0: initialization c ii=1: determining dense regions in space of individual events c x,y,rad: tranverse coordinates and radius of particle i c s: long coordinate c ii=2: plot of individual event c xpar1: itau ; valid: 1,..,10 c xpar2: s-range: -xpar2 < s < xpar2 c xpar3, x-range: -xpar3 < x < xpar3 c xpar4, y-range: -xpar4 < y < xpar4 c----------------------------------------------------------------------- include "epos.inc" double precision tpro,zpro,ttar,ztar,ttaus,detap,detat common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat parameter (mxcoos=60) common/cdensh/kdensh(matau,mxcoos,mxcoox,mxcooy),ktot(matau) character cy*3 dlcoox=0.5 dlcooy=0.5 dlcoos=0.5 if(ii.eq.0)then do i=1,mxcoos do j=1,mxcoox do k=1,mxcooy kdensh(itau,i,j,k)=0 enddo enddo enddo ktot(itau)=0 elseif(ii.eq.1)then if(itau.lt.1.or.itau.gt.mxtau)return tau=sngl(ttaus) z=s/tau do 1 i=1,mxcoos si=-mxcoos*dlcoos/2-dlcoos/2+i*dlcoos do 2 j=1,mxcoox xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox do 3 k=1,mxcooy yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy if(((x-xj)**2+(y-yk)**2+(z-si)**2).gt.rad**2)goto3 kdensh(itau,i,j,k)=kdensh(itau,i,j,k)+1 ktot(itau)=ktot(itau)+1 3 continue 2 continue 1 continue elseif(ii.eq.2)then itaux=nint(xpar1) if(itaux.gt.mxtau)stop'STOP in xjden2: itaux too large' is=nint(xpar2/dlcoos) ix=nint(xpar3/dlcoox) iy=nint(xpar4/dlcooy) if(is.gt.mxcoos/2)stop'STOP in xjden2: s-range too large' if(ix.gt.mxcoox/2)stop'STOP in xjden2: x-range too large' if(iy.gt.mxcooy/2)stop'STOP in xjden2: y-range too large' do k=mxcooy/2+1-iy,mxcooy/2+iy write(cy,'(f3.1)')-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy write(ifhi,'(a)') 'openhisto' write(ifhi,'(a,2f7.2)')'xrange',-is*dlcoos,is*dlcoos write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox write(ifhi,'(a)') 'set ityp2d 5' write(ifhi,'(a)') 'txt "xaxis [z] "' write(ifhi,'(a)') *'txt "yaxis x (fm), y='//cy//' fm"' write(ifhi,'(a,i4)') 'array2d',2*is do j=mxcoox/2+1-ix,mxcoox/2+ix do i=mxcoos/2+1-is,mxcoos/2+is write(ifhi,'(e11.3)') *kdensh(itaux,i,j,k)/dlcooy/dlcoos/dlcoox/ktot(itaux) enddo enddo write(ifhi,'(a)') ' endarray' write(ifhi,'(a)') 'closehisto plot2d' enddo else stop'STOP in xjden2: wrong option' endif return end cc----------------------------------------------------------------------- c subroutine postscript(iii,ii,ic) cc----------------------------------------------------------------------- c include 'epos.inc' c character*10 color(10) c if(iii.eq.0)then c ifps=21 c open(unit=ifps,file='zzz.ps',status='unknown') c WRITE(ifps,'(a)') '%!PS-Adobe-2.0' c WRITE(ifps,'(a)') '%%Title: tt2.fig' c WRITE(ifps,'(a)') '%%Orientation: Portrait' c WRITE(ifps,'(a)') '%%BeginSetup' c WRITE(ifps,'(a)') '%%IncludeFeature: *PageSize A4' c WRITE(ifps,'(a)') '%%EndSetup' c WRITE(ifps,'(a)') '%%EndComments' c WRITE(ifps,*) '/l {lineto} bind def' c WRITE(ifps,*) '/rl {rlineto} bind def' c WRITE(ifps,*) '/m {moveto} bind def' c WRITE(ifps,*) '/rm {rmoveto} bind def' c WRITE(ifps,*) '/s {stroke} bind def' c WRITE(ifps,*) '/gr {grestore} bind def' c WRITE(ifps,*) '/gs {gsave} bind def' c WRITE(ifps,*) '/cp {closepath} bind def' c WRITE(ifps,*) '/tr {translate} bind def' c WRITE(ifps,*) '/sc {scale} bind def' c WRITE(ifps,*) '/sd {setdash} bind def' c WRITE(ifps,*) '/sdo {[.01 .05] 0 sd} bind def' c WRITE(ifps,*) '/sdf {[1 .0] 0 sd} bind def' c WRITE(ifps,*) '/n {newpath} bind def' c WRITE(ifps,*) '/slw {setlinewidth } bind def' c write(ifps,*) '/srgb {setrgbcolor} bind def' c write(ifps,*) '/lgrey { 0 0.95 0.95 srgb} bind def' c write(ifps,*) '/black { 0 0 0 srgb} bind def' c write(ifps,*) '/red { 1 0 0 srgb} bind def ' c write(ifps,*) '/green { 0 1 0 srgb} bind def ' c write(ifps,*) '/blue { 0 0 1 srgb} bind def ' c write(ifps,*) '/yellow { 1 0.5 0 srgb} bind def ' c write(ifps,*) '/turquoise { 0 1 1 srgb} bind def ' c write(ifps,*) '/purple { 1 0 1 srgb} bind def ' cc write(ifps,*) '/ { srgb} bind def ' cc write(ifps,*) '/ { srgb} bind def ' c write(ifps,*) '/ef {eofill} bind def' c WRITE(ifps,'(a)') '%%EndProlog' c WRITE(ifps,*) 'gsave' c WRITE(ifps,*) '/Helvetica findfont 10 scalefont setfont' c color(9)='lgrey ' c color(1)='black ' c color(2)='red ' c color(3)='green ' c color(4)='blue ' c color(7)='yellow ' c color(5)='turquoise ' c color(6)='purple ' c np=0 c elseif(iii.eq.1)then c np=np+1 c write(ifps,'(a,i4)') '%%Page: number ',np c write(ifps,'(a)') 'gsave' c WRITE(ifps,*) '100 700 tr' c scale=0.125 c WRITE(ifps,*) 1./scale,1./scale,' sc' c WRITE(ifps,*) scale/2.,' slw' c WRITE(ifps,*) '/Helvetica findfont ',15.*scale c & ,' scalefont setfont' c write(ifps,*) color(1),' n ',smin,xmin,' m ( tau:',tau,') show ' c c WRITE(ifps,*) '/Helvetica findfont ',5.*scale c & ,' scalefont setfont' c c c yb=-2. c dy=4./12. c yb=yb-dy/2 c do iyb=0,11 c yb=yb+dy c WRITE(ifps,*) 'gsave' c WRITE(ifps,*) (xmax-xmin)*1.1*float(int(iyb/4)) c & ,-(xmax-xmin)*1.1*mod(iyb,4),' tr' c write(ifps,*) ' n ',smin,xmin,' m ',smax,xmin,' l ' c & ,smax,xmax,' l ',smin,xmax,' l cp s ' cc.......particles in layer iyb............. c do i=1,nptl c if(ii.gt.0)then c write(ifps,*) color(mod(i,5)+2) c & ,' n ',u,x-r,' m ',o,x-r,' l ' c & ,o,x+r,' l ',u,x+r,' l cp s ' c write(ifps,*) ' n ',u,x-r,' m (',i,ior,') show ' c else c write(ifps,*) ' n ',s,x,r,0,360,' arc ',color(ic),' s ' c write(ifps,*) ' n ',s-r,x,' m (',i,io,') show ' c endif c 10 enddo c write(ifps,*) color(1),' n ',smin,xmin,' m (',yb,') show ' c WRITE(ifps,*) 'grestore' c enddo !yb bin c write(ifps,'(a)') 'grestore' c write(ifps,*) 'showpage' c elseif(iii.eq.2)then c write(ifps,*) 'gr' c c write(ifps,'(a)') '%%Trailer' c write(ifps,'(a,i4)') '%%Pages: ',np c write(ifps,'(a)') '%%EOF' c close(unit=ifps) c endif c c return c end c c------------------------------------------------------------------------------ subroutine xtauev(iii) c------------------------------------------------------------------------------ jdum=iii end c------------------------------------------------------------------------------ subroutine wimi c------------------------------------------------------------------------------ end c------------------------------------------------------------------------------ subroutine wimino c------------------------------------------------------------------------------ end c------------------------------------------------------------------------------ subroutine xspace(iii) c------------------------------------------------------------------------------ jdum=iii end c------------------------------------------------------------------------------ subroutine wclu c------------------------------------------------------------------------------ end c------------------------------------------------------------------------------ subroutine wclufi c------------------------------------------------------------------------------ end c------------------------------------------------------------------------------ subroutine wtime(iii) c------------------------------------------------------------------------------ jdum=iii end