File indexing completed on 2024-04-06 12:14:03
0001
0002 subroutine bjinta(ier)
0003
0004
0005
0006 include 'epos.inc'
0007 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
0008 common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctimel/ntc
0009 common/col3/ncol,kolpt
0010 double precision ttaun,ttau0,rcproj,rctarg
0011 common/cttaun/ttaun /cttau0/ttau0 /geom1/rcproj,rctarg
0012 logical go,lclean
0013
0014 call utpri('bjinta',ish,ishini,4)
0015
0016 ier=0
0017
0018 if(ncol.eq.0.and.iappl.eq.2)goto1000
0019 if(nevt.ne.1.or.ifrade.eq.0)goto1000
0020
0021 if(iappl.eq.4.or.iappl.eq.9)then
0022 goto5000
0023 endif
0024
0025 !if(iappl.eq.1)then
0026 ! tauxx=0.7+0.94*max(radnuc(maproj),radnuc(matarg))/(0.5*engy)
0027 !else
0028 ! tauxx=0.
0029 !endif
0030 !tauzz=max(taumin,tauxx)
0031 !print*,'====',taumin,tauxx,tauzz
0032 !ttaus=dble(tauzz)
0033 ttaus=taumin
0034 ttau0=dsqrt(rcproj*rctarg)
0035 call jtauin ! initialize hyperbola
0036
0037 if(iappl.ne.1)goto 5000
0038
0039
0040
0041
0042 if(iorsce.eq.0.and.iorsdf.eq.0.and.iorshh.eq.0
0043 & .or.iorsdf.eq.3)then
0044 if(iorsdf.eq.3)then
0045 lclean=.false.
0046 if(nclean.gt.0.and.nptl.gt.mxptl/5)then
0047 ! if nptl already very big, clean up useless particles in cptl list.
0048 !(do not use it when gakstr() is called (some information lost)
0049 nptli=maproj+matarg+1
0050 do iii=nptli,nptl
0051 go=.true.
0052 if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
0053 if(go.and.mod(istptl(iii),10).ne.0)istptl(iii)=99
0054 enddo
0055 nptl0=nptl
0056 call utclea(nptli,nptl0)
0057 lclean=.true.
0058 endif
0059 nptlbpo=nptl
0060 call jintpo(lclean,iret) !parton-ladder-fusion
0061 if(ish.ge.2)call alist('parton-ladder-fusion&',nptlbpo+1,nptl)
0062 if(iret.eq.1)goto 1001
0063 endif
0064 goto 5000
0065 else
0066 stop'bjinta: not supported any more (310305). '
0067 endif
0068
0069 5000 continue
0070
0071 nptlbd=nptl
0072
0073 call xSpaceTime
0074
0075 if(ifrade.eq.0)goto779 !skip decay
0076 if(idecay.eq.0)goto779 !skip decay
0077
0078
0079 if(ish.ge.2)call alist('final decay&',0,0)
0080 if(iappl.eq.4.or.iappl.eq.7.or.iappl.eq.9)then
0081 nptli=1
0082 else
0083 nptli=maproj+matarg+1
0084 endif
0085 np1=nptli
0086 41 np2=nptl
0087 nptli=np1
0088 ip=np1-1
0089 do while (ip.lt.np2)
0090 ip=ip+1
0091 if(istptl(ip).eq.0)then
0092 call hdecas(ip,iret)
0093 if(iret.eq.1)goto 1001
0094 if(iret.eq.-1)goto 42
0095
0096 if(nclean.gt.0.and.nptl.gt.mxptl/2)then
0097 nnnpt=0
0098 do iii=nptli,ip
0099 go=.true.
0100 if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
0101 if(go.and.mod(istptl(iii),10).ne.0)then
0102 istptl(iii)=99
0103 nnnpt=nnnpt+1
0104 endif
0105 enddo
0106 if(nnnpt.gt.mxptl-nptl)then
0107 nptl0=nptl
0108 call utclea(nptli,nptl0)
0109 np2=np2-nnnpt
0110 ip=ip-nnnpt
0111 nptli=ip
0112 endif
0113 endif
0114 endif
0115 42 continue
0116 enddo
0117 nptli=max(nptli,np1)
0118 np1=np2+1
0119 if(np1.le.nptl)then
0120 if(ish.ge.2)then
0121 if(ish.ge.3)call alist('partial list&',0,0)
0122 do 6 ip=np1,nptl
0123 call alist('&',ip,ip)
0124 6 continue
0125 endif
0126 goto 41
0127 endif
0128 779 continue
0129
0130
0131
0132
0133
0134
0135
0136 1000 continue
0137 call utprix('bjinta',ish,ishini,4)
0138 return
0139
0140 1001 continue
0141 ier=1
0142 goto 1000
0143
0144 end
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528
0529
0530
0531
0532
0533
0534
0535
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633
0634
0635
0636
0637
0638
0639
0640
0641
0642
0643
0644
0645
0646
0647
0648
0649
0650
0651
0652
0653
0654
0655
0656
0657
0658
0659
0660
0661
0662
0663
0664
0665
0666
0667
0668
0669
0670
0671
0672
0673
0674
0675
0676
0677
0678
0679
0680
0681
0682
0683
0684
0685
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813
0814
0815
0816
0817
0818
0819
0820
0821
0822
0823
0824
0825
0826
0827
0828
0829
0830
0831
0832
0833
0834
0835
0836
0837
0838
0839
0840
0841
0842
0843
0844
0845
0846
0847
0848
0849
0850
0851
0852
0853
0854
0855
0856
0857
0858
0859
0860
0861
0862
0863
0864
0865
0866
0867
0868
0869
0870
0871
0872
0873
0874
0875
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899
0900
0901
0902
0903
0904
0905
0906
0907
0908
0909
0910
0911
0912
0913
0914
0915
0916
0917
0918
0919
0920
0921
0922
0923
0924
0925
0926
0927
0928
0929
0930
0931
0932
0933
0934
0935
0936
0937
0938
0939
0940
0941
0942
0943
0944
0945
0946
0947
0948
0949
0950
0951
0952
0953
0954
0955
0956
0957
0958
0959
0960
0961
0962
0963
0964
0965
0966
0967
0968
0969
0970
0971
0972
0973
0974
0975
0976
0977
0978
0979
0980
0981
0982
0983
0984
0985
0986
0987
0988
0989
0990
0991
0992
0993
0994
0995
0996
0997
0998
0999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897 subroutine jintpo(lclean,iret)
1898
1899
1900
1901 include 'epos.inc'
1902 include 'epos.incems'
1903 include 'epos.incico'
1904 common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
1905 *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
1906 parameter(m1grid=10,kgrid=3,kegrid=3)
1907 parameter(m3xgrid=21*kgrid*kegrid)
1908 parameter(mxcl=4000,mxcli=1000)
1909 integer idropgrid(m1grid,m1grid,m3xgrid)
1910 & ,jdropgrid(m1grid,m1grid,m3xgrid)
1911 & ,jclu(m3xgrid),nsegp4(mxcl)
1912 & ,irep(mxcl),mmji(mxcl,mxcli)
1913 & ,jccl(mxcl,nflav,2),nseg(mxcl),mseg(mxcl),kclu(mxcl)
1914 & ,naseg(0:mxcl),nfseg(mxcl),nsegmx(mxcl)
1915 & ,jc(nflav,2),ke(6),jcjj(nflav,2)
1916 & ,nsegmt(mxptl),kmean(2,1001:1000+mxcl)
1917 & ,ncluj(1001:1000+mxptl) !,ic(2),nclk(m3xgrid)
1918 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
1919 & ,pptld(5,mxcl),p4tmp
1920 double precision enp,pzp,ena,pza,ccc,sss,enew,pznew,taa2,pta2
1921 .,pold,pnew,pnew2,p5a2,eeta,enf,pzf
1922 common /cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
1923 common/cdelzet/delzet,delsgr /cvocell/vocell
1924 common/cranphi/ranphi
1925 parameter(maxp=40*mxcl)
1926 dimension yrad(maxp),phirad(maxp),pe(5),yrad2(maxp),nfrag(mxptl)
1927 logical first,lnew(m1grid,m1grid),lold(m1grid,m1grid),lcont,lpass
1928 &,lclean
1929 double precision ptest(5),ttest,p52,p4mean(4,mxcl),zor,tor,xmxms
1930 &,ppp(5),be(4),energ,bp,rmean(2,mxcl) !,qqq(5),www(5)
1931 dimension mpair(mamx,mamx)
1932 save ncntpo!,ntrr,ntrt
1933 data ncntpo /0/!,ntrr/0/,ntrt/0/
1934 ncntpo=ncntpo+1
1935 call utpri('jintpo',ish,ishini,4)
1936 tauico=ttaus
1937
1938 if(kigrid.gt.kgrid)then
1939 write(ifmt,*)'kigrid,kgrid:',kigrid,kgrid
1940 stop'jintpo: kigrid too big. '
1941 endif
1942
1943 fegrid=yhaha/5.36 !rapidity range/rap range at RHIC
1944 if(fegrid.gt.kegrid)then
1945 write(ifmt,*)'fegrid,kegrid:',fegrid,kegrid
1946 stop'jintpo: kegrid too small. '
1947 endif
1948
1949 m3grid=m3xgrid /kgrid *kigrid /kegrid *fegrid
1950
1951 if(mod(m3grid,2).eq.0)m3grid=m3grid+1
1952
1953 xgrid=8
1954 sgrid=fsgrid *ttaus *fegrid *6.0
1955 delxgr=2*xgrid/m1grid
1956 delsgr=2*sgrid/m3grid
1957 shiftx=(1.-2.*rangen())*rangen()*delxgr/2.
1958 shifts=(1.-2.*rangen())*rangen()*delsgr/2.
1959 vocell= delxgr * delxgr * delsgr
1960 delzet=delsgr/ttaus
1961 nptla=nptl
1962 iflhc=1
1963 if(iLHC.eq.1)iflhc=min(3,max(1,nsegce-1))
1964 nsegsuj=max(nsegce,nint(float(nsegsu/iflhc)*fegrid))
1965 xmxms=150d0 !maximum mass for a subcluster
1966
1967 if(ncntpo.eq.1)then
1968 write(ifmt,*)'EPOS used with FUSION option'
1969 if(ish.ge.1)then
1970 print*,'+',('-',i=1,57)
1971 print*,'| ttaus sgrid nsegce:',ttaus,fsgrid,' ',nsegce
1972 print*,'| sgrid delxgr delsgr vocell:',sgrid,delxgr,delsgr,vocell
1973 print*,'+',('-',i=1,57)
1974 endif
1975 endif
1976
1977
1978
1979
1980 if(ish.ge.6)write(ifch,*)'compute x,y,z'
1981 do n=1,nptla
1982 iaaptl(n)=1
1983 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1984 ! the following check is important for ioclude= 4 or 5
1985 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1986 if(ioclude.ne.3.and.istptl(n).eq.10)then
1987 stop'\n\n remnant cluster detected in jintpo\n\n'
1988 endif
1989 if(istptl(n).eq.0.or.istptl(n).eq.10)then
1990 call jtain(n,xptl(n),yptl(n),zptl(n),tptl(n),nnn,0)
1991 call idtau(idptl(n),pptl(4,n),pptl(5,n),taurem)
1992 if(abs(taurem).gt.1e-6.and.
1993 . ityptl(n).ge.40.and.ityptl(n).le.59.and.nnn.eq.1)then
1994 iaaptl(n)=0 !!!!nnn=1: ptl lives later
1995
1996 endif
1997 call jtaus(zptl(n),tz,sptl(n))
1998 strap=1e10
1999 xpl=tptl(n)+zptl(n)
2000 xmi=tptl(n)-zptl(n)
2001 if(xmi.gt.0.0.and.xpl.gt.0.0)then
2002 strap=0.5*log(xpl/xmi)
2003 else !particle at eta=inf
2004 iaaptl(n)=0
2005 endif
2006 dezptl(n)=strap !space-time-rapidity
2007
2008
2009 else
2010 iaaptl(n)=0
2011 xptl(n)=0.
2012 yptl(n)=0.
2013 zptl(n)=0.
2014 sptl(n)=0.
2015 dezptl(n)=1e10
2016 endif
2017 enddo
2018 ntry=0
2019
2020
2021
2022
2023 ptmax=ptclu**2
2024 nptlo=1
2025 if(iLHC.eq.1)nptlo=maproj+matarg+1
2026 do n=1,nptla
2027 !random number calls before if's
2028 !to keep random number sequence for testing
2029 rdm=rangen()
2030 if(istptl(n).eq.0)then
2031 call idquac(n,idum1,idum2,idum3,jc)
2032 else
2033 jc(4,1)=0
2034 jc(4,2)=0
2035 endif
2036 if(iaaptl(n).ne.0)then
2037 if(istptl(n).ne.10)then
2038 pt2=pptl(1,n)**2+pptl(2,n)**2
2039 am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2
2040 if(n.le.nptlo)then !spectators
2041 iaaptl(n)=0
2042 elseif(ityptl(n).eq.47.or.ityptl(n).eq.57)then
2043 iaaptl(n)=0 !no active spectators
2044 elseif(istptl(n).ne.0.or.ityptl(n).ge.60)then
2045 iaaptl(n)=0
2046 elseif(abs(idptl(n)).le.100)then
2047 iaaptl(n)=0 !no fundamental particle (electron...)
2048 !elseif(ityptl(n).eq.41.or.ityptl(n).eq.51)then
2049 ! iaaptl(n)=0 !avoid particles coming from remnant droplet decay
2050 ! !(already droplet decay products !)
2051 elseif(abs(am2tmp-pptl(5,n)*pptl(5,n)).gt.30.)then
2052 iaaptl(n)=0 !to discard off shell particles
2053 elseif(jc(4,1).ne.0.or.jc(4,2).ne.0)then
2054 iaaptl(n)=0 !to discard charmed particles
2055 elseif(iLHC.eq.1)then
2056 if(pt2.lt.1.e-3)then
2057 iaaptl(n)=0 !to discard slow proton (spectators)
2058 elseif(ptclu.gt.0.)then
2059 if(pt2.gt.(1.5*ptclu)**2)then
2060 ptmax=max(ptmax,pt2)
2061 if(ioquen.eq.0.or.fploss.le.0.)then !high pt particle escape
2062 iaaptl(n)=0
2063 elseif(fploss.gt.1.e10)then !high pt particle absorbed
2064 iaaptl(n)=1
2065 else
2066 iaaptl(n)=-1 !high pt particle lose energy
2067 endif
2068 elseif(pt2.gt.(0.5*ptclu)**2)then
2069 if(rdm.lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then
2070 ptmax=max(ptmax,pt2)
2071 if(ioquen.eq.0.or.fploss.le.0.)then !high pt particle escape
2072 iaaptl(n)=0
2073 elseif(fploss.gt.1.e10)then !high pt particle absorbed
2074 iaaptl(n)=1
2075 else
2076 iaaptl(n)=-1 !high pt particle lose energy
2077 endif
2078 endif
2079 endif
2080 else
2081 if(ioquen.eq.0.or.fploss.le.0.)then !high pt particle escape
2082 iaaptl(n)=0
2083 elseif(fploss.gt.1.e10)then !high pt particle absorbed
2084 iaaptl(n)=1
2085 else
2086 iaaptl(n)=-1 !high pt particle lose energy
2087 endif
2088 endif
2089 elseif(pt2.gt.(1.5*ptclu)**2)then
2090 ptmax=max(ptmax,pt2)
2091 if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)then
2092 iaaptl(n)=0
2093 endif
2094 elseif(pt2.gt.(0.5*ptclu)**2)then
2095 if(rdm.lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then
2096 if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)iaaptl(n)=0
2097 endif
2098 endif
2099 endif
2100 endif
2101 c if(iaaptl(n).eq.0.and.abs(idptl(n)).le.10000.and.abs(idptl(n))
2102 c & .ge.100.and.mod(abs(idptl(n)),100).ne.0
2103 c & .and.ityptl(n).ne.0.and.ityptl(n).lt.50)
2104 c & print*,'not selected',idptl(n),ityptl(n),pptl(4,n)
2105 enddo
2106
2107 c...Start cluster formation
2108
2109 8888 continue
2110
2111 nsegsuj=max(nsegce,nint(nsegsuj/1.08**ntry))
2112 ntry=ntry+1
2113 if(ntry.gt.90) !do not put more than 100 or change limit for p4mean
2114 &call utstop('jintpo: cluster formation impossible ! &')
2115 nptl=nptla
2116
2117 do 1 k=1,m3grid
2118 do 1 j=1,m1grid
2119 do 1 i=1,m1grid
2120 idropgrid(i,j,k)=0
2121 1 continue
2122
2123 c...count string segments in cell
2124
2125 if(ish.ge.6)write(ifch,*)'count string segments in cell'
2126 do n=1,nptla
2127 if(iaaptl(n).ne.0)then
2128 i=1+(xptl(n)+xgrid+shiftx)/delxgr
2129 j=1+(yptl(n)+xgrid+shiftx)/delxgr
2130 k=1+(sptl(n)+sgrid+shifts)/delsgr
2131 if( i.ge.1.and.i.le.m1grid
2132 & .and.j.ge.1.and.j.le.m1grid
2133 & .and.k.ge.1.and.k.le.m3grid)then
2134 nfrag(n)=1
2135 if(iLHC.eq.1)then !count number of quarks to absorbe more baryon than mesons
2136 call idquac(n,idum1,idum2,idum3,jc)
2137 nfrag(n)=0
2138 do nj=1,2
2139 do ni=1,nflav
2140 nfrag(n)=nfrag(n)+ni*jc(ni,nj)
2141 enddo
2142 enddo
2143
2144 endif
2145 c print *,idptl(n),nfrag(n)
2146 idropgrid(i,j,k)=idropgrid(i,j,k)+nfrag(n)
2147 endif
2148 endif
2149 enddo
2150
2151
2152 if(iLHC.eq.1)then
2153 cc...check high pt segments
2154 c !to use this part one has to define:
2155 c !... 1 = valid
2156 c !... -1 = valid but high pt
2157 c !... 0 = not valid
2158 c
2159
2160 c esu=0
2161 c do i=1,nptla
2162 c if(istptl(i).eq.0.or.istptl(i).eq.3)esu=esu+pptl(4,i)
2163 c enddo
2164 c write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu
2165
2166 if(ish.ge.6)write(ifch,*)'check high pt segments'
2167 if(fploss.gt.0.0)then
2168 ein=0
2169 elo=0
2170 do n=1,nptla
2171 delen=0.
2172 if(iaaptl(n).eq.-1)then
2173 i=1+(xptl(n)+xgrid+shiftx)/delxgr
2174 j=1+(yptl(n)+xgrid+shiftx)/delxgr
2175 k=1+(sptl(n)+sgrid+shifts)/delsgr
2176 if( i.ge.1.and.i.le.m1grid
2177 & .and.j.ge.1.and.j.le.m1grid
2178 & .and.k.ge.1.and.k.le.m3grid)then
2179 iescape=0
2180 xa=xptl(n)
2181 ya=yptl(n)
2182 eta=dezptl(n)
2183 pz=pptl(3,n)
2184 en=pptl(4,n)
2185 taa2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2
2186 pza=pz
2187 ena=sqrt(taa2+pza**2)
2188 enaxx=en
2189 p5a2=pptl(5,n)**2
2190 pta2=pptl(1,n)**2+pptl(2,n)**2
2191 !transform to bjo
2192 eeta=eta
2193 ccc=cosh(eeta)
2194 sss=sinh(eeta)
2195 enp= ena*ccc-pza*sss
2196 pzp=-ena*sss+pza*ccc
2197 vz=pzp/enp
2198 vx=pptl(1,n)/enp
2199 vy=pptl(2,n)/enp
2200 !print*,'+++++++',eta,vz,pz/en
2201 if(vx.ne.0.0.or.vy.ne.0.0)then
2202 if(abs(vx).ge.abs(vy))then
2203 ica=1
2204 rat=vy/vx
2205 is=sign(1.,vx)
2206 l=i
2207 va=xa
2208 wa=ya
2209 else
2210 ica=2
2211 rat=vx/vy
2212 is=sign(1.,vy)
2213 l=j
2214 va=ya
2215 wa=xa
2216 endif
2217 if(is.eq.-1)then
2218 imax=1
2219 else !if(is.eq. 1)
2220 imax=m1grid
2221 endif
2222 vr=sqrt(vx**2+vy**2)
2223 ratx=vz/vr
2224 qu=fploss*delxgr*sqrt(1+rat**2)*sqrt(1+ratx**2)
2225 do lu=l,imax,is
2226 vu=-(xgrid+shiftx)+(lu-0.5)*delxgr
2227 wu=wa+rat*(vu-va)
2228 mu=1+(wu+xgrid+shiftx)/delxgr
2229 if(mu.ge.1.and.mu.le.m1grid)then
2230 if(ica.eq.1)then
2231 ix=lu
2232 jx=mu
2233 else
2234 ix=mu
2235 jx=lu
2236 endif
2237 if(idropgrid(i,j,k).ge.nsegce)then
2238 dens=float(idropgrid(ix,jx,k))/float(iflhc)
2239 else
2240 dens=0
2241 endif
2242 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2243 ! DeltaE ~ sqrt(T**3*E) * dL (BDMPS2008)
2244 ! -> DeltaE ~ dens^3/8 * sqrt(E) *dL
2245 !del=dens**(3./8.)*max(1.,sqrt(pptl(4,n)))
2246 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2247 escale=6
2248 xen=enp/escale
2249 del=0.02*dens**(3./8.)*max(1.,sqrt(xen))
2250 delen=delen+del*escale*qu
2251 c if(k.ge.m3grid/2-2.and.k.le.m3grid/2+4)
2252 c & write(ifch,'(2i3,4x,2i3,4x,2i3,4x,i3,3f7.3)')
2253 c & k,ica,i,j,ix,jx,idropgrid(ix,jx,k),qu,delen
2254 c & ,pptl(4,n)-delen
2255 endif
2256 enddo
2257 endif
2258 pold=sqrt(pta2+pzp**2)
2259 enew=enp-delen
2260 pnew2=enew**2-p5a2
2261 if(enew.gt.0..and.pnew2.gt.0.)iescape=1
2262 if(iescape.eq.1)then
2263 iaaptl(n)=0
2264 idropgrid(i,j,k)=idropgrid(i,j,k)-nfrag(n)
2265 if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible. '
2266 pnew=sqrt(pnew2)
2267 p1=pnew*pptl(1,n)/pold
2268 p2=pnew*pptl(2,n)/pold
2269 pznew=pnew*pzp /pold
2270 enf= enew*ccc+pznew*sss
2271 pzf= enew*sss+pznew*ccc
2272 if(ish.ge.1.and.enf.gt.ena.or.enf.ne.enf)then
2273 if(ish.ge.1.and.(enf-ena.gt.0.1d0.or.enf.ne.enf))then
2274 write(ifmt,*)'WWWW en gain for eta =',eta
2275 write(ifmt,*)'e pz bf (lab)',ena,pza,enaxx
2276 write(ifmt,*)'e pz bf (bjo)',enp,pzp ,enp**2-pzp**2
2277 write(ifmt,*)'e pz af (bjo)',enew,pznew,enew**2-pznew**2
2278 write(ifmt,*)'e pz af (lab)',enf,pzf ,enf**2-pzf**2
2279 c stop
2280 endif
2281 elseif(delen.gt.1e-3)then
2282 c create fake particle with energy lost in cluster
2283 nptl=nptl+1
2284 nptla=nptla+1
2285 istptl(nptl)=3 !daughter part of cluster
2286 iaaptl(nptl)=1
2287 pptl(1,nptl)=pptl(1,n)-p1
2288 pptl(2,nptl)=pptl(2,n)-p2
2289 pptl(3,nptl)=pptl(3,n)-pzf
2290 pptl(4,nptl)=pptl(4,n)-enf
2291 do l=1,3
2292 xorptl(l,nptl)=xorptl(l,n)
2293 enddo
2294 ! mother
2295 pptl(1,n)=p1
2296 pptl(2,n)=p2
2297 pptl(3,n)=pzf
2298 tivptl(1,n)=xorptl(4,n)
2299 pptl(4,n)=sqrt(pptl(1,n)**2+pptl(2,n)**2
2300 & +pptl(3,n)**2+pptl(5,n)**2)
2301 if(abs(pptl(4,n)-enf).gt.0.01*pptl(4,n))
2302 & print *,'jintpo eloss',pptl(4,n),enf,(pptl(4,n)-enf)/pptl(4,n)
2303 call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
2304 tivptl(2,n)=tivptl(1,n)+taugm*(-alog(rangen()))
2305 ifrptl(1,n)=nptl
2306 ! daughter
2307 pptl(5,nptl)=0.
2308 pptl(4,nptl)=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2
2309 & +pptl(3,nptl)**2+pptl(5,nptl)**2)
2310 xorptl(4,nptl)=xorptl(4,n)
2311 tivptl(1,nptl)=tivptl(1,n)
2312 iorptl(nptl)=n
2313 jorptl(nptl)=0
2314 idptl(nptl)=idptl(n)*100+sign(99,idptl(n))
2315 ityptl(nptl)=ityptl(n)+2
2316 xptl(nptl)=xptl(n)
2317 yptl(nptl)=yptl(n)
2318 zptl(nptl)=zptl(n)
2319 sptl(nptl)=sptl(n)
2320 dezptl(nptl)=dezptl(n)
2321 if(ish.ge.6)write(ifch,*)'--> Virtual particle : ',nptl
2322 & ,idptl(nptl),iorptl(nptl),ityptl(nptl),(pptl(kk,nptl),kk=1,5)
2323 endif
2324 elo= elo+ena-pptl(4,n)
2325 !write(ifch,*)n,' escape'
2326 else
2327 iaaptl(n)=1
2328 ein=ein+pptl(4,n)
2329 !write(ifch,*)n,'core'
2330 endif
2331 endif
2332 endif
2333 enddo
2334 else
2335 do n=1,nptla
2336 if(iaaptl(n).eq.-1)then
2337 iaaptl(n)=0
2338 endif
2339 enddo
2340 endif
2341
2342 c eloss=esu
2343 c esu=0
2344 c do i=1,nptla
2345 c if(istptl(i).eq.0.or.istptl(i).eq.3)esu=esu+pptl(4,i)
2346 c enddo
2347 c write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu
2348 c eloss=eloss-esu
2349 c write(ifch,*)'+++++ein,elo,eloss ',ein,elo,eloss
2350
2351
2352 c do n=1,nptla
2353 c if(iaaptl(n).eq.-1)then
2354 c i=1+(xptl(n)+xgrid+shiftx)/delxgr
2355 c j=1+(yptl(n)+xgrid+shiftx)/delxgr
2356 c k=1+(sptl(n)+sgrid+shifts)/delsgr
2357 c if( i.ge.1.and.i.le.m1grid
2358 c & .and.j.ge.1.and.j.le.m1grid
2359 c & .and.k.ge.1.and.k.le.m3grid)then
2360 c if(idropgrid(i,j,k).lt.2*nsegce)then !surface segments
2361 c iaaptl(n)=0
2362 c idropgrid(i,j,k)=idropgrid(i,j,k)-1
2363 c if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible. '
2364 cc else
2365 cc iaaptl(n)=1
2366 c endif
2367 c endif
2368 c endif
2369 c enddo
2370
2371 endif
2372
2373 c...identify clusters
2374
2375 ifirst=0
2376 if(ish.ge.6)write(ifch,*)'identify clusters'
2377 do k=1,m3grid !~~~~~~k-loop
2378 jjj=0
2379 do j=1,m1grid
2380 first=.true.
2381 do i=1,m1grid
2382 if(idropgrid(i,j,k).ge.nsegce)then
2383 if(first)then
2384 ifirst=i
2385 jjj=jjj+1
2386 if(jjj.gt.mxcl)stop'jintpo: mxcl too small. '
2387 irep(jjj)=0
2388 jj=jjj
2389 first=.false.
2390 endif
2391 jdropgrid(i,j,k)=jj
2392 if(j.gt.1)then
2393 if(jdropgrid(i,j-1,k).ne.0)then
2394 jjo=jdropgrid(i,j-1,k)
2395 if(jjo.lt.jj)then
2396 if(jj.eq.jjj)jjj=jjj-1
2397 jjx=jj
2398 jj=jjo
2399 do ii=ifirst,i
2400 jdropgrid(ii,j,k)=jj
2401 if(jdropgrid(ii,j-1,k).eq.jjx)then
2402 if(jjx.gt.jjj)jjj=jjj+1
2403 jja=jjx
2404 jjb=jj
2405 90 continue
2406 if(irep(jja).eq.0.or.irep(jja).eq.jjb)then
2407 irep(jja)=jjb
2408 else
2409 mn=min(irep(jja),jjb)
2410 mx=max(irep(jja),jjb)
2411 irep(jja)=mn
2412 jja=mx
2413 jjb=mn
2414 goto 90
2415 endif
2416 endif
2417 enddo
2418 elseif(jdropgrid(i,j-1,k).gt.jj)then
2419 irep(jjo)=jj
2420 endif
2421 endif
2422 endif
2423 else
2424 jdropgrid(i,j,k)=0
2425 first=.true.
2426 endif
2427 enddo
2428 enddo
2429 !~~~~cluster jj ---> cluster irep(jj)
2430 do jj=jjj,1,-1
2431 if(irep(jj).ne.0)then
2432 do j=1,m1grid
2433 do i=1,m1grid
2434 if(jdropgrid(i,j,k).eq.jj)jdropgrid(i,j,k)=irep(jj)
2435 enddo
2436 enddo
2437 endif
2438 enddo
2439 !~~~~~remove empty cluster indices
2440 jjjx=jjj
2441 jjj=0
2442 jj=0
2443 do jjx=1,jjjx
2444 if(irep(jjx).eq.0)then
2445 jj=jj+1
2446 jjj=jjj+1
2447 else
2448 do j=1,m1grid
2449 do i=1,m1grid
2450 if(jdropgrid(i,j,k).gt.jj)
2451 & jdropgrid(i,j,k)=jdropgrid(i,j,k)-1
2452 enddo
2453 enddo
2454 endif
2455 enddo
2456 !~~~~~
2457 jclu(k)=jjj
2458 enddo !~~~~~~~~~~~~~~~~~ END k-loop
2459
2460 c...add segments to avoid holes
2461
2462 if(ish.ge.6)write(ifch,*)'add segments from holes'
2463 if(iohole.eq.1.and.iLHC.eq.1)then
2464 do 83 n=1,nptla
2465 if(iaaptl(n).eq.0)goto 83
2466 ihit=0
2467 i=1+(xptl(n)+xgrid+shiftx)/delxgr
2468 j=1+(yptl(n)+xgrid+shiftx)/delxgr
2469 k=1+(sptl(n)+sgrid+shifts)/delsgr
2470 if( i.ge.1.and.i.le.m1grid
2471 & .and.j.ge.1.and.j.le.m1grid
2472 & .and.k.ge.1.and.k.le.m3grid)then
2473 jgr=jdropgrid(i,j,k)
2474 nplus=idropgrid(i,j,k)
2475 if(jgr.eq.0)then
2476 isgi=1
2477 if(rangen().gt.0.5)isgi=-1
2478 isgj=1
2479 if(rangen().gt.0.5)isgj=-1
2480 isgk=1
2481 if(rangen().gt.0.5)isgk=-1
2482 do ii=i-isgi,i+isgi,isgi
2483 do jj=j-isgj,j+isgj,isgj
2484 do kk=k-isgk,k+isgk,2*isgk
2485 if( ii.ge.1.and.ii.le.m1grid
2486 . .and.jj.ge.1.and.jj.le.m1grid
2487 . .and.kk.ge.1.and.kk.le.m3grid)then
2488 if(jdropgrid(ii,jj,kk).gt.0)nplus=nplus+1
2489 c nplus=idropgrid(ii,jj,kk)+iflhc
2490 if(nplus.ge.nsegce)then
2491 ihit=1
2492 c goto 84
2493 endif
2494 endif
2495 enddo
2496 enddo
2497 enddo
2498 endif
2499 endif
2500 c 84 continue
2501 if(ihit.eq.1)then
2502 idropgrid(i,j,k)=nplus!idropgrid(i,j,k)+iflhc
2503 c add cell randomly to adjacent cluster if any
2504 isgi=1
2505 if(rangen().gt.0.5)isgi=-1
2506 isgj=1
2507 if(rangen().gt.0.5)isgj=-1
2508 do ii=i-isgi,i+isgi,2*isgi
2509 do jj=j-isgj,j+isgj,2*isgj
2510 if( ii.ge.1.and.ii.le.m1grid
2511 . .and.jj.ge.1.and.jj.le.m1grid)then
2512 jjj=jdropgrid(ii,jj,k)
2513 if(jjj.gt.0)then
2514 jdropgrid(i,j,k)=jjj
2515 goto 83
2516 endif
2517 endif
2518 enddo
2519 enddo
2520 c if no adjacent cluster, create a new one
2521 jclu(k)=jclu(k)+1
2522 jdropgrid(i,j,k)=jclu(k)
2523 endif
2524 83 continue
2525 endif
2526
2527 if(ish.ge.8)then
2528 write(ifch,*)'segment positions'
2529 do k=1,m3grid
2530 write(ifch,*)'k=',k,' jclu=',jclu(k)
2531 & ,' s=',(k-1)*delsgr-sgrid-shifts
2532 do j=m1grid,1,-1
2533 write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,m1grid)
2534 & ,(jdropgrid(i,j,k),i=1,m1grid)
2535 enddo
2536 enddo
2537 endif
2538
2539 c...absolute clusters numbering (for all k)
2540
2541 if(ish.ge.6)write(ifch,*)'absolute clusters numbering'
2542 if(iLHC.eq.-1)then !fuse touching clusters along k (not used)
2543 c midbin=m3grid/2
2544 c if(rangen().gt.0.5)midbin=midbin+1
2545 jjj=1000
2546 k=1
2547 jclu(k)=0
2548 do j=1,m1grid
2549 do i=1,m1grid
2550 if(jdropgrid(i,j,k).gt.0)then
2551 jclu(k)=jclu(k)+1
2552 jjj=jjj+1
2553 ncluj(jjj)=0
2554 kmean(1,jjj)=k
2555 kmean(2,jjj)=k
2556 jjx=jdropgrid(i,j,k)
2557 do jj=j,m1grid
2558 do ii=i,m1grid
2559 if(jdropgrid(ii,jj,k).eq.jjx)then
2560 jdropgrid(ii,jj,k)=jjj
2561 ncluj(jjj)=ncluj(jjj)+idropgrid(ii,jj,k)
2562 endif
2563 enddo
2564 enddo
2565 endif
2566 enddo
2567 enddo
2568 do k=2,m3grid
2569 do j=1,m1grid
2570 do i=1,m1grid
2571 lnew(i,j)=.true.
2572 lold(i,j)=.true.
2573 enddo
2574 enddo
2575 jclu(k)=0
2576 do j=1,m1grid
2577 do i=1,m1grid
2578 if(lnew(i,j).and.jdropgrid(i,j,k).gt.0)then
2579 lcont=.false.
2580 if(jdropgrid(i,j,k-1).gt.0)then
2581 jj0=jdropgrid(i,j,k-1)
2582 nlim=0
2583 if(jj0.gt.1000)nlim=ncluj(jj0)
2584 if(nlim.lt.nsegsuj)then
2585 jclu(k)=jclu(k)+1
2586 jjx=jdropgrid(i,j,k)
2587 if(jjx.gt.jj0)then
2588 do jj=jjx,jjj
2589 ncluj(jj)=ncluj(jj+1)
2590 kmean(1,jj)=kmean(1,jj+1)
2591 kmean(2,jj)=kmean(2,jj+1)
2592 enddo
2593 jjj=jjj-1
2594 jclu(k)=jclu(k)-1
2595 endif
2596 kmean(2,jj0)=k
2597 do jj=1,m1grid
2598 do ii=1,m1grid
2599 if(lnew(ii,jj))then
2600 if(jdropgrid(ii,jj,k).eq.jjx)then
2601 jdropgrid(ii,jj,k)=jj0
2602 ncluj(jj0)=ncluj(jj0)+idropgrid(ii,jj,k)
2603 lnew(ii,jj)=.false.
2604 elseif(jdropgrid(ii,jj,k).gt.jjx.and.jjx.gt.jj0)then
2605 jdropgrid(ii,jj,k)=jdropgrid(ii,jj,k)-1
2606 endif
2607 endif
2608 enddo
2609 enddo
2610 else
2611 lcont=.true.
2612 endif
2613 else
2614 lcont=.true.
2615 endif
2616 if(lcont.and.lold(i,j))then
2617 jclu(k)=jclu(k)+1
2618 jjj=jjj+1
2619 ncluj(jjj)=0
2620 kmean(1,jjj)=k
2621 kmean(2,jjj)=k
2622 jjx=jdropgrid(i,j,k)
2623 do jj=j,m1grid
2624 do ii=i,m1grid
2625 if(jdropgrid(ii,jj,k).eq.jjx)then
2626 jdropgrid(ii,jj,k)=jjj
2627 ncluj(jjj)=ncluj(jjj)+idropgrid(ii,jj,k)
2628 lold(ii,jj)=.false.
2629 endif
2630 enddo
2631 enddo
2632 endif
2633 endif
2634 enddo
2635 enddo
2636 enddo
2637 jjs=0
2638 do jj=1001,jjj
2639 jjs=jjs+1
2640 kclu(jjs)=(kmean(1,jj)+kmean(2,jj))/2
2641 nseg(jjs)=0
2642 nsegp4(jjs)=0
2643 p4mean(1,jjs)=0d0
2644 p4mean(2,jjs)=0d0
2645 p4mean(3,jjs)=0d0
2646 p4mean(4,jjs)=0d0
2647 rmean(1,jjs)=0d0 !cluster distance to center
2648 rmean(2,jjs)=0d0 !cluster maximal radius
2649 enddo
2650 do k=1,m3grid
2651 do j=1,m1grid
2652 do i=1,m1grid
2653 if(jdropgrid(i,j,k).gt.0)then
2654 jdropgrid(i,j,k)=jdropgrid(i,j,k)-1000
2655 endif
2656 enddo
2657 enddo
2658 enddo
2659 jjj=jjj-1000
2660 else
2661 jjj=jclu(1)
2662 do k=2,m3grid
2663 do j=1,m1grid
2664 do i=1,m1grid
2665 if(jdropgrid(i,j,k).gt.0)
2666 & jdropgrid(i,j,k)=jdropgrid(i,j,k)+jjj
2667 enddo
2668 enddo
2669 jjj=jjj+jclu(k)
2670 enddo
2671 jjs=0
2672 do k=1,m3grid
2673 do jj=1,jclu(k)
2674 jjs=jjs+1
2675 kclu(jjs)=k
2676 nseg(jjs)=0
2677 nsegp4(jjs)=0
2678 p4mean(1,jjs)=0d0
2679 p4mean(2,jjs)=0d0
2680 p4mean(3,jjs)=0d0
2681 p4mean(4,jjs)=0d0
2682 rmean(1,jjs)=0d0 !cluster distance to center
2683 rmean(2,jjs)=0d0 !cluster maximal radius
2684 enddo
2685 enddo
2686 endif
2687 do 20 i=1,maproj
2688 do 20 j=1,matarg
2689 20 mpair(j,i)=0
2690
2691
2692 c...calculate mean energy of segments going into in clusters
2693
2694 if(ish.ge.6)write(ifch,*)
2695 &'calculate mean energy of segments going into in clusters'
2696 rmax=0.
2697 do 95 n=1,nptla
2698 if(iaaptl(n).eq.0)goto 95
2699 i=1+(xptl(n)+xgrid+shiftx)/delxgr
2700 j=1+(yptl(n)+xgrid+shiftx)/delxgr
2701 k=1+(sptl(n)+sgrid+shifts)/delsgr
2702 if( i.ge.1.and.i.le.m1grid
2703 & .and.j.ge.1.and.j.le.m1grid
2704 & .and.k.ge.1.and.k.le.m3grid)then
2705 jj=jdropgrid(i,j,k)
2706 if(jj.gt.0)then
2707 nsegp4(jj)=nsegp4(jj)+1
2708 io=iorptl(n)
2709 if(iLHC.eq.1.and.ityptl(n).lt.40.and.io.gt.0)then !look for corresponding pair of nucleons
2710 if(iorptl(io).gt.0)then
2711 do while(iorptl(iorptl(io)).ge.0)
2712 io=iorptl(io)
2713 enddo
2714 ip=iorptl(io)
2715 it=jorptl(io)-maproj
2716 if(ip.gt.0.and.it.gt.0)mpair(it,ip)=mpair(it,ip)+1
2717 endif
2718 endif
2719 do kk=1,4
2720 p4mean(kk,jj)=p4mean(kk,jj)+dble(pptl(kk,n))
2721 enddo
2722 rr=sqrt(xptl(n)**2+yptl(n)**2)
2723 rmean(1,jj)=rmean(1,jj)+rr
2724 rmax=max(rr,rmax)
2725 c if(jj.eq.9)write(ifch,*)'after',n
2726 c & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k
2727 elseif(istptl(n).eq.3)then !virtual particle is lost (no cluster here)
2728 iaaptl(n) = 0 !don't use this particle next time
2729 istptl(n) = 5
2730 ior=iorptl(n)
2731 ifrptl(1,ior) = 0
2732 ifrptl(2,ior) = 0
2733 do k=1,3 !restore energy to mother particle
2734 pptl(k,ior)=pptl(k,ior)+pptl(k,n)
2735 enddo
2736 pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2
2737 & +pptl(3,ior)**2+pptl(5,ior)**2)
2738 call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm)
2739 tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen()))
2740 endif
2741 endif
2742 95 continue
2743
2744 ectot=0.
2745 amctot=0.
2746 maptot=0
2747 mapmax=0
2748 npair=0
2749 do jj=1,jjj
2750 ectot=ectot+p4mean(4,jj)
2751 amctot=amctot+(p4mean(4,jj)+p4mean(3,jj))
2752 & *(p4mean(4,jj)-p4mean(3,jj))-p4mean(1,jj)**2-p4mean(2,jj)**2
2753
2754 if(nsegp4(jj).ge.nsegce/iflhc)then
2755 p4mean(4,jj)=p4mean(4,jj)/dble(nsegp4(jj))
2756 rmean(1,jj)=rmean(1,jj)/dble(nsegp4(jj))
2757 else
2758 p4mean(4,jj)=1d50
2759 rmean(1,jj)=-1d0
2760 endif
2761 enddo
2762 yco=0.
2763 ycor=1.
2764 if(iLHC.eq.1.and.amctot.gt.0.)then
2765 do ip=1,maproj
2766 do it=1,matarg
2767 if(mpair(it,ip).gt.0)npair=npair+1
2768 mapmax=max(mapmax,mpair(it,ip))
2769 maptot=maptot+mpair(it,ip)
2770 enddo
2771 enddo
2772 amctot=sqrt(amctot)
2773 if(amctot.gt.amuseg**2)then
2774 visco=1.
2775 yrmaxi=yradmx
2776
2777 yrmaxi=yrmaxi*log(exp(yradpi)+amctot**2)
2778
2779 if(yrmaxi.gt.1e-2)then
2780 yyrmax=dble(yrmaxi)
2781 fradflii=sngl(1d0/
2782 & ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
2783 else
2784 yrmaxi=0.
2785 fradflii=1.
2786 endif
2787 else
2788 amctot=1.
2789 visco=0.
2790 fradflii=1.
2791 yrmaxi=ainfin !to define ityptl as 19 if mass is too low (aumin=amuseg+yrmaxi in epos-hnb.f)
2792 endif
2793 if(ylongmx.lt.0.)then
2794 delzet=abs(ylongmx)*log(exp(yradmi)+amctot**2)
2795
2796
2797 yco=delzet !* 1.75
2798 else
2799 yco=ylongmx
2800 endif
2801 ycor=yco
2802 if(fvisco.gt.0.)ycor=0.
2803 if(ycor.gt.1.e-2)then
2804 ycor=sqrt(sinh(ycor)/ycor)/fradflii
2805 else
2806 ycor=1.
2807 endif
2808 else
2809 if(ylongmx.lt.0)delzet=1.75*delzet
2810 endif
2811
2812
2813
2814 if(ish.ge.6)write(ifch,*)'mark segments going into in clusters'
2815 do 96 n=1,nptla
2816 if(iaaptl(n).eq.0)goto 96
2817 i=1+(xptl(n)+xgrid+shiftx)/delxgr
2818 j=1+(yptl(n)+xgrid+shiftx)/delxgr
2819 k=1+(sptl(n)+sgrid+shifts)/delsgr
2820 if( i.ge.1.and.i.le.m1grid
2821 & .and.j.ge.1.and.j.le.m1grid
2822 & .and.k.ge.1.and.k.le.m3grid)then
2823 jj=jdropgrid(i,j,k)
2824 if(jj.gt.0)then
2825 ! count only particles with reasonnable energy
2826 pt2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2
2827 am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2
2828 if(abs(am2tmp).lt.0.1.or.
2829 & dble(pptl(4,n)).le.100.d0/dble(ntry)*p4mean(4,jj))then
2830 istptl(n) = 3
2831 nseg(jj)=nseg(jj)+1
2832 if(rmean(1,jj).gt.0d0)then
2833 rmean(2,jj)=max(rmean(2,jj)
2834 & ,abs(rmean(1,jj)-sqrt(xptl(n)**2+yptl(n)**2)))
2835
2836 endif
2837
2838
2839 else
2840
2841 if(ish.ge.1)write(ifch,*)'Rejected particle : ',n,idptl(n)
2842 & ,ityptl(n),(pptl(kk,n),kk=1,5),am2tmp
2843 & ,dble(pptl(4,n))/dble(nsegp4(jj)),p4mean(4,jj),nsegp4(jj),jj
2844 nsegp4(jj)=nsegp4(jj)-1
2845 iaaptl(n)=0
2846
2847
2848 if(nsegp4(jj).ge.nsegce/iflhc)then
2849 p4mean(4,jj)=(p4mean(4,jj)*dble(nsegp4(jj)+1)-pptl(4,n))
2850 & /dble(nsegp4(jj))
2851 else
2852 p4mean(4,jj)=1d50
2853 endif
2854 endif
2855 endif
2856 endif
2857 96 continue
2858
2859
2860
2861
2862 if(ish.ge.6)write(ifch,*)'add segments moving towards clusters'
2863 if(iocluin.eq.1)then
2864 do 93 n=1,nptla
2865 if(iaaptl(n).eq.0)goto 93
2866 ihit=0
2867 i=1+(xptl(n)+xgrid+shiftx)/delxgr
2868 j=1+(yptl(n)+xgrid+shiftx)/delxgr
2869 k=1+(sptl(n)+sgrid+shifts)/delsgr
2870 if( i.ge.1.and.i.le.m1grid
2871 & .and.j.ge.1.and.j.le.m1grid
2872 & .and.k.ge.1.and.k.le.m3grid)then
2873 jgr=jdropgrid(i,j,k)
2874 if(jgr.eq.0)then
2875 if(i.ge.m1grid/2)then
2876 isi=-1
2877 else
2878 isi=1
2879 endif
2880 if(j.ge.m1grid/2)then
2881 jsi=-1
2882 else
2883 jsi=1
2884 endif
2885 do ii=i,i+2*isi,isi
2886 do jj=j,j+2*jsi,jsi
2887 if(.not.(ii.eq.i.and.jj.eq.j))then
2888 if(ii.ge.1.and.ii.le.m1grid
2889 . .and.jj.ge.1.and.jj.le.m1grid)then
2890 jg=jdropgrid(ii,jj,k)
2891 if(jg.gt.0)then
2892
2893 x=xptl(n)
2894 y=yptl(n)
2895 vrad=( x*pptl(1,n)/pptl(4,n)+y*pptl(2,n)/pptl(4,n))
2896 if(vrad.lt.0.)then
2897 ihit=1
2898 goto 94
2899 endif
2900
2901 endif
2902 endif
2903 endif
2904 enddo
2905 enddo
2906 endif
2907 endif
2908 94 continue
2909 if(ihit.eq.1)then
2910 delx=delxgr*(ii-i)
2911 dely=delxgr*(jj-j)
2912 xn=xptl(n)+delx
2913 yn=yptl(n)+dely
2914 ix=1+(xn+xgrid+shiftx)/delxgr
2915 jx=1+(yn+xgrid+shiftx)/delxgr
2916 jgrx=jdropgrid(ix,jx,k)
2917 if(jgrx.gt.0)then
2918 xptl(n)=xn
2919 yptl(n)=yn
2920 istptl(n) = 3
2921 nseg(jgrx)=nseg(jgrx)+1
2922 if(rmean(1,jg).gt.0d0)
2923 & rmean(2,jgrx)=max(rmean(2,jgrx)
2924 & ,abs(rmean(1,jgrx)-sqrt(xptl(n)**2+yptl(n)**2)))
2925 endif
2926 endif
2927 93 continue
2928 endif
2929
2930
2931
2932
2933
2934 if(ish.ge.6)write(ifch,*)'decay Strings not included in clusters'
2935 if(ifrade.ne.0.and.ntry.eq.1.and..not.lclean)then
2936 nptlx=nptl+1
2937
2938 do n=1,nptl
2939 if(istptl(n).eq.29.and.ityptl(n).lt.40)then !fragmented central strings
2940 iused=0
2941 do nn=ifrptl(1,n),ifrptl(2,n)
2942 if(istptl(nn).eq.3.or.ifrptl(1,nn).ne.0)iused=iused+1
2943 enddo
2944 if(iused.eq.0)then !if no fragment used, reset string
2945 istptl(n)=28
2946 do nn=ifrptl(1,n),ifrptl(2,n) !suppress product
2947 istptl(nn)=4
2948 iaaptl(nn)=0
2949 enddo
2950 do nn=iorptl(n),jorptl(n) !reinitialize string ends
2951 istptl(nn)=20
2952 enddo
2953 endif
2954 endif
2955 enddo
2956 call gakfra(0,iret) !fragmentation using Z
2957 if(iret.gt.0)goto 1000
2958 do nn=nptlx,nptl !exclude new particle from cluster formation
2959 iaaptl(nn)=0
2960 call jtain(nn,xptl(nn),yptl(nn),zptl(nn),tptl(nn),nnn,0)
2961 call jtaus(zptl(nn),tz,sptl(nn))
2962 strap=1e10
2963 xpl=tptl(nn)+zptl(nn)
2964 xmi=tptl(nn)-zptl(nn)
2965 if(xmi.gt.0.0.and.xpl.gt.0.0)strap=0.5*log(xpl/xmi)
2966 dezptl(nn)=strap !space-time-rapidity
2967 enddo
2968 maxfra=nptl
2969 nptla=nptl
2970 if(ish.ge.6.and.nptl.ne.nptlx-1)
2971 & call alist('list after second fragmentation&',nptlx,nptl)
2972 if(irescl.eq.1)then
2973 call utghost(iret)
2974 if(iret.gt.0)goto 1000
2975 endif
2976 endif
2977
2978
2979
2980
2981
2982
2983
2984
2985
2986
2987
2988
2989
2990
2991
2992
2993
2994
2995
2996
2997
2998
2999
3000
3001
3002
3003
3004
3005
3006
3007
3008
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
3024
3025
3026
3027
3028
3029
3030
3031
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
3051
3052
3053
3054
3055
3056
3057
3058
3059
3060
3061
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084 nptlb=nptl+jjj
3085
3086
3087 if(ish.ge.6)write(ifch,*)'prepare /cptl/ for clusters'
3088 do jj=1,jjj
3089 nptl=nptl+1
3090 istptl(nptl)=12
3091 do l=1,4
3092 pptl(l,nptl)=0.
3093 xorptl(l,nptl)=0
3094 enddo
3095 sptl(nptl)=0
3096 uptl(nptl)=0
3097 optl(nptl)=0
3098 desptl(nptl)=0
3099 iorptl(nptl)=0
3100 jorptl(nptl)=0
3101
3102
3103 nsegmx(jj)=max(1,nint(float(nseg(jj))/float(nsegsuj)))
3104 if(ish.ge.6)write(ifch,*)'nsegmx',jj,nseg(jj),nsegmx(jj)
3105 & ,nsegsuj
3106 enddo
3107
3108
3109
3110 if(ish.ge.6)write(ifch,*)'prepare /cptl/ for subclusters'
3111 mm=0
3112 do jj=1,jjj
3113 do ii=1,nsegmx(jj)
3114 mm=mm+1
3115 if(mm.gt.mxcl)stop'jintpo: mxcl too small. '
3116 mmji(jj,ii)=mm
3117 mseg(mm)=0
3118 nptl=nptl+1
3119 istptl(nptl)=10
3120 do l=1,4
3121 pptld(l,mm)=0d0
3122 pptl(l,nptl)=0.
3123 xorptl(l,nptl)=0
3124 enddo
3125 sptl(nptl)=0
3126 uptl(nptl)=0
3127 optl(nptl)=0
3128 desptl(nptl)=0
3129 do l=1,nflav
3130 jccl(mm,l,1)=0
3131 jccl(mm,l,2)=0
3132 enddo
3133 iorptl(nptl)=nptla+jj
3134 jorptl(nptl)=0
3135 if(ii.eq.1)ifrptl(1,nptla+jj)=nptlb+mm
3136 ifrptl(2,nptla+jj)=nptlb+mm
3137 enddo
3138 enddo
3139
3140
3141
3142 if(ish.ge.6)write(ifch,*)'separate string segments'
3143 do 98 n=1,nptla
3144 if(istptl(n).ne.3)goto 98
3145 i=1+(xptl(n)+xgrid+shiftx)/delxgr
3146 j=1+(yptl(n)+xgrid+shiftx)/delxgr
3147 k=1+(sptl(n)+sgrid+shifts)/delsgr
3148 if( i.ge.1.and.i.le.m1grid
3149 & .and.j.ge.1.and.j.le.m1grid
3150 & .and.k.ge.1.and.k.le.m3grid)then
3151 jj=jdropgrid(i,j,k)
3152 if(jj.gt.0)then
3153 iimx=nsegmx(jj)
3154 do ii=1,iimx
3155 mm=mmji(jj,ii)
3156 if(mseg(mm).eq.0)goto 10 !not to have an empty cluster
3157 enddo
3158 ii=1+rangen()*iimx
3159 ii=min(ii,iimx)
3160 iini=ii
3161 am2tmpmx=1e30
3162 9 ntmp=mmji(jj,ii)
3163 am2tmp=(pptld(4,ntmp)+pptld(3,ntmp))
3164 & *(pptld(4,ntmp)-pptld(3,ntmp))
3165 if(am2tmp.gt.xmxms**2/2.)then
3166 if(am2tmp.lt.am2tmpmx)then
3167 mm=mmji(jj,ii)
3168 am2tmpmx=am2tmp
3169 endif
3170 ii=ii+1
3171 if(ii.gt.iimx)ii=1
3172 if(ii.ne.iini)then
3173 goto 9
3174 else
3175 goto 10
3176 endif
3177 endif
3178 mm=mmji(jj,ii)
3179 10 continue
3180 mseg(mm)=mseg(mm)+1
3181 ifrptl(1,n)=mm !local use of ifrptl
3182
3183
3184 p4tmp=0d0
3185 do l=1,3
3186 pptld(l,mm)= pptld(l,mm) + dble(pptl(l,n))
3187 p4tmp=p4tmp+dble(pptl(l,n))*dble(pptl(l,n))
3188 enddo
3189 p4tmp=sqrt(p4tmp+dble(pptl(5,n)*pptl(5,n)))
3190 pptld(4,mm)= pptld(4,mm) + p4tmp
3191
3192
3193
3194
3195 xorptl(1,nptlb+mm)=xorptl(1,nptlb+mm)+xptl(n)
3196 xorptl(2,nptlb+mm)=xorptl(2,nptlb+mm)+yptl(n)
3197 xorptl(3,nptlb+mm)=xorptl(3,nptlb+mm)+zptl(n)
3198 xorptl(4,nptlb+mm)=xorptl(4,nptlb+mm)+tptl(n)
3199 sptl(nptlb+mm)=sptl(nptlb+mm)+sptl(n)
3200 aa=cos(phievt)
3201 bb=sin(phievt)
3202 cc=-sin(phievt)
3203 dd=cos(phievt)
3204 xrot=xptl(n)*aa+yptl(n)*bb
3205 yrot=xptl(n)*cc+yptl(n)*dd
3206 uptl(nptlb+mm)=uptl(nptlb+mm)+xrot**2
3207 optl(nptlb+mm)=optl(nptlb+mm)+yrot**2
3208 desptl(nptlb+mm)=desptl(nptlb+mm)+xrot*yrot
3209 call idquac(n,idum1,idum2,idum3,jc)
3210
3211
3212
3213
3214
3215
3216
3217
3218
3219 do l=1,nflav
3220 jccl(mm,l,1)=jccl(mm,l,1)+jc(l,1)
3221 jccl(mm,l,2)=jccl(mm,l,2)+jc(l,2)
3222 enddo
3223 else
3224 idropgrid(i,j,k)=0
3225 endif
3226 endif
3227 98 continue
3228
3229
3230
3231 if(ish.ge.6)write(ifch,*)'associate segments to clusters'
3232 naseg(0)=0
3233 do jj=1,jjj
3234 do ii=1,nsegmx(jj)
3235 mm=mmji(jj,ii)
3236 naseg(mm)=naseg(mm-1)+mseg(mm)
3237 nfseg(mm)=0
3238 enddo
3239 enddo
3240 do 97 n=1,nptla
3241 if(istptl(n).ne.3)goto 97
3242 istptl(n) = 2
3243
3244
3245 mm=ifrptl(1,n)
3246 nfseg(mm)=nfseg(mm)+1
3247 nsegmt(naseg(mm-1)+nfseg(mm))=n
3248 97 continue
3249 do jj=1,jjj
3250 nst=0
3251 do ii=1,nsegmx(jj)
3252 mm=mmji(jj,ii)
3253 if(mseg(mm).ne.nfseg(mm))stop'jintpo: mseg.ne.nfseg '
3254 nst=nst+mseg(mm)
3255 enddo
3256 if(nst.ne.nseg(jj))stop'sum(mseg(mm)).ne.nseg(jj)'
3257 enddo
3258
3259
3260
3261 if(ish.ge.6)write(ifch,*)'finish cluster storage to /cptl/'
3262 xx=0.
3263 yy=0.
3264 xy=0.
3265 mjjsegsum=0
3266 do jj=1,jjj
3267 njj=nptla+jj
3268 mjjseg=0
3269 do l=1,nflav
3270 jcjj(l,1)=0
3271 jcjj(l,2)=0
3272 enddo
3273 do ii=1,nsegmx(jj)
3274 mm=mmji(jj,ii)
3275 n=nptlb+mm
3276
3277 do l=1,nflav
3278 jc(l,1)=jccl(mm,l,1)
3279 jc(l,2)=jccl(mm,l,2)
3280 ke(l)=jc(l,1)-jc(l,2)
3281 jcjj(l,1)=jcjj(l,1)+jc(l,1)
3282 jcjj(l,2)=jcjj(l,2)+jc(l,2)
3283 enddo
3284 call idenct(jc,idptl(n)
3285 * ,ibptl(1,n),ibptl(2,n),ibptl(3,n),ibptl(4,n))
3286 ttest=0d0
3287 do ji=1,4
3288 ptest(ji)=0d0
3289 do ns=naseg(mm-1)+1,naseg(mm)
3290 ni=nsegmt(ns)
3291 ptest(ji)=ptest(ji)+dble(pptl(ji,ni))
3292 enddo
3293 ptest(ji)=abs(ptest(ji)-pptld(ji,mm))
3294 ttest=ttest+ptest(ji)
3295 enddo
3296 amcmin=1.01*utamnu(ke(1),ke(2),ke(3),ke(4),ke(5),ke(6),4)
3297 p52=((pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm))
3298 & -pptld(1,mm)**2-pptld(2,mm)**2)
3299 if(iLHC.eq.1)then
3300 amcmin=sqrt(max(amcmin**2,sngl(p52)))
3301 if(amcmin/ycor.gt.amuseg)amcmin=amcmin*ycor
3302
3303 endif
3304
3305 pptld(5,mm)=0d0
3306 if(p52.gt.0d0)then
3307 pptld(5,mm)=sqrt(p52)
3308 jerr(2)=jerr(2)+1
3309 if(iLHC.eq.1.and.pptld(5,mm).lt.amcmin
3310 & .and.pptld(4,mm).gt.amcmin)then
3311 pptld(5,mm)=dble(amcmin)
3312 bp=sqrt((pptld(4,mm)+pptld(5,mm))*(pptld(4,mm)-pptld(5,mm))
3313 & /(pptld(3,mm)*pptld(3,mm)+pptld(2,mm)*pptld(2,mm)
3314 & +pptld(1,mm)*pptld(1,mm)))
3315 pptld(1,mm)=bp*pptld(1,mm)
3316 pptld(2,mm)=bp*pptld(2,mm)
3317 pptld(3,mm)=bp*pptld(3,mm)
3318
3319
3320
3321
3322
3323
3324 endif
3325 elseif(p52.le.0d0)then
3326 jerr(3)=jerr(3)+1
3327 pptld(5,mm)=dble(amcmin)
3328 pptld(4,mm)=sqrt(pptld(3,mm)*pptld(3,mm)
3329 & +pptld(2,mm)*pptld(2,mm)
3330 & +pptld(1,mm)*pptld(1,mm)
3331 & +pptld(5,mm)*pptld(5,mm))
3332 endif
3333 if(ish.ge.1.and.(abs(ttest).gt.1.d0.or.pptld(5,mm).gt.xmxms))
3334 & then
3335 call utmsg('jintpo&')
3336 write(ifmt,*)'***** Warning in jintpo !',ntry
3337 write(ifch,*)'***** jintpo: momenta messed up (ttest > 0)'
3338 write(ifch,*)'*****',mm,n,mseg(mm),p52,ttest
3339 write(ifch,*)'*****',jj,nsegp4(jj),p4mean(4,jj)
3340 write(ifch,'(a,10x,4f15.4)')'*****',(pptld(ji,mm),ji=1,4)
3341 do ns=naseg(mm-1)+1,naseg(mm)
3342 ni=nsegmt(ns)
3343 write(ifch,'(a,i5,i9,5f15.4,f12.4)')'*****',ni,idptl(ni)
3344 * ,(pptl(ji,ni),ji=1,4),pptl(5,ni)**2
3345 * ,(pptl(4,ni)+pptl(3,ni))*(pptl(4,ni)-pptl(3,ni))
3346 * -pptl(1,ni)**2-pptl(2,ni)**2
3347 enddo
3348 endif
3349 if(pptld(5,mm).gt.xmxms)then
3350 p4max=0.
3351 nh=0
3352 do ns=naseg(mm-1)+1,naseg(mm)
3353 ni=nsegmt(ns)
3354 if(pptl(4,ni).ge.p4max)then
3355 nh=ni
3356 p4max=pptl(4,ni)
3357 endif
3358 enddo
3359 if(nh.le.0)then
3360 stop'Cannot be in jintpo ...'
3361 else !put back nh as normal particle
3362 iaaptl(nh) = 0 !don't use this fragment next time
3363 c if(iaaptl(nh).eq.0)print*,'excluded',idptl(nh),ityptl(nh)
3364 c & ,pptl(4,nh)
3365 if(mod(abs(idptl(nh)),100).eq.99)then !restore lost energy
3366 istptl(nh) = 5
3367 ior=iorptl(nh)
3368 iaaptl(ior)=0 !don't use this fragment next time
3369 ifrptl(1,ior) = 0
3370 ifrptl(2,ior) = 0
3371 if(istptl(ior).eq.0)then
3372 do k=1,3
3373 pptl(k,ior)=pptl(k,ior)+pptl(k,n)
3374 enddo
3375 pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2
3376 & +pptl(3,ior)**2+pptl(5,ior)**2)
3377 call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm)
3378 tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen()))
3379 else
3380 istptl(ior) = 0
3381 endif
3382 elseif(idptl(nh).lt.1e4)then
3383 istptl(nh) = 0 !particle
3384 ifrptl(1,nh) = 0
3385 ifrptl(2,nh) = 0
3386 else
3387 istptl(nh) = 10 !droplet
3388 endif
3389 endif
3390 if(ish.ge.1)
3391 & write(ifch,*)'***** Redo cluster without heavy particle :'
3392 & ,nh,ntry
3393 goto 8888
3394 endif
3395 do l=1,5
3396 pptl(l,n)=sngl(pptld(l,mm))
3397 enddo
3398 mjjseg=mjjseg+mseg(mm)
3399 do l=1,4
3400 pptl(l,njj)=pptl(l,njj)+pptl(l,n)
3401 xorptl(l,njj)=xorptl(l,njj)+xorptl(l,n)
3402 xorptl(l,n)=xorptl(l,n)/float(mseg(mm))
3403 enddo
3404 sptl(njj)=sptl(njj)+sptl(n)
3405 uptl(njj)=uptl(njj)+uptl(n)
3406 optl(njj)=optl(njj)+optl(n)
3407 desptl(njj)=desptl(njj)+desptl(n)
3408 sptl(n)=sptl(n)/float(mseg(mm))
3409 uptl(n)=uptl(n)/float(mseg(mm))
3410 optl(n)=optl(n)/float(mseg(mm))
3411 desptl(n)=desptl(n)/float(mseg(mm))
3412 radptl(n)=0
3413 istptl(n)=10
3414 ifrptl(1,n)=0
3415 ifrptl(2,n)=0
3416 tivptl(1,n)=xorptl(4,n)
3417 tivptl(2,n)=ainfin
3418 ityptl(n)=60
3419 radptl(n)=0
3420 dezptl(n)=0.
3421 enddo !ii
3422 do l=1,4
3423 xorptl(l,njj)=xorptl(l,njj)/float(mjjseg)
3424 enddo
3425 mjjsegsum=mjjsegsum+mjjseg
3426 xx=xx+uptl(njj)
3427 yy=yy+optl(njj)
3428 xy=xy+desptl(njj)
3429 sptl(njj)=sptl(njj)/float(mjjseg)
3430 uptl(njj)=uptl(njj)/float(mjjseg)
3431 optl(njj)=optl(njj)/float(mjjseg)
3432 desptl(njj)=desptl(njj)/float(mjjseg)
3433 pjj52=(pptl(4,njj)+pptl(3,njj))*(pptl(4,njj)-pptl(3,njj))
3434 & -pptl(1,njj)**2-pptl(2,njj)**2
3435 pptl(5,njj)=0
3436 if(pjj52.gt.0)then
3437 pptl(5,njj)=sqrt(pjj52)
3438 endif
3439 ityptl(njj)=60
3440 call idenct(jc,idptl(njj)
3441 * ,ibptl(1,njj),ibptl(2,njj),ibptl(3,njj),ibptl(4,njj))
3442 enddo !jj
3443
3444
3445
3446
3447 ranphi=0
3448 rini=0.
3449 if(mjjsegsum.gt.0)then
3450 xx=xx/float(mjjsegsum)
3451 yy=yy/float(mjjsegsum)
3452 xy=xy/float(mjjsegsum)
3453 dta=0.5*(xx-yy)
3454
3455
3456
3457
3458
3459
3460
3461
3462
3463
3464
3465
3466 if(xy.lt.0..and.dta.ne.0.)then
3467 ranphi=0.5*atan(-xy/dta)
3468 elseif(xy.gt.0..and.dta.ne.0.)then
3469 ranphi=-0.5*atan(xy/dta)
3470 else
3471 ranphi=0
3472 endif
3473
3474
3475
3476
3477
3478
3479
3480
3481
3482
3483
3484
3485
3486
3487 rini=max(0.01,sqrt(5./3.*(xx+yy))) !<r**2>=3/5*R**2 for sphere of radius R
3488 volu=(4./3.*pi*(xx+yy)**1.5)
3489
3490 flowpp=0.
3491 flowaa=0.
3492
3493 if(iLHC.eq.1.and.visco.gt.0.)then
3494 if(npair.gt.0)then
3495 fcorr=min(1.,float(mapmax)*abs(fvisco)/float(maptot))
3496 visco=min(1.,2.*float(maptot)/float(npair)/float(mapmax))**2
3497
3498 elseif(lclean)then !large number of particles, npair can't be calculated
3499 visco=1e-6
3500 fcorr=1.
3501 else !cluster from remnants only
3502 visco=1.
3503 fcorr=1.
3504 endif
3505 if(visco.ge.1.)yrmaxi=0. !yrmaxi*(1.-visco)
3506 c visco=exp(-min(50.,max(0.,
3507 c & float(koievt)/log(amctot)-abs(fvisco))**yradpi)) !mix flow
3508 c & max(0.,rmax**2-abs(fvisco)))) !mix flow
3509 c yrmaxi=log(amctot**2)
3510 c yrmaxi=yradmx*yrmaxi*(1.-visco)
3511 c if(visco.lt.1.and.yrmaxi.gt.1e-2)then
3512 c yyrmax=dble(yrmaxi)
3513 c fradflii=sngl(1d0/
3514 c & ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
3515 c else
3516 c visco=1.
3517 c yrmaxi=0.
3518 c fradflii=1.
3519 c endif
3520 c if(visco.gt.1e-5)then
3521 c yrmaxi=yradmx*(yrmaxi
3522 c & +visco*(log(amctot)*yradpx/yradmx-yrmaxi))
3523 c else
3524 c yrmaxi=yradmx*yrmaxi
3525 c endif
3526 fradflii=1.
3527 if(yrmaxi.gt.0)then
3528 flowpp=visco*log(fcorr*amctot)*yradpx
3529 flowaa=yrmaxi
3530 if(rangen().lt.flowaa/flowpp)then
3531 visco=0.
3532 yrmaxi=flowaa+max(0.,flowpp-flowaa)
3533 yyrmax=dble(yrmaxi)
3534 fradflii=sngl(1d0/
3535 & ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
3536 else
3537 yrmaxi=0.
3538 visco=log(fcorr*amctot)/log(amctot)
3539 endif
3540 elseif(fcorr.lt.1.)then
3541 visco=log(fcorr*amctot)/log(amctot)
3542 endif
3543 endif
3544
3545 if(ish.ge.3)write(ifch,*)'yrmaxi,delzet=',yrmaxi,delzet
3546 c print *,'->',bimevt,yrmaxi,visco*yradpx*log(amctot),flowaa
3547 c & ,flowpp,maptot,log(ectot),visco,log(amctot),rho
3548 c & ,mapmax,npair
3549 c & ,min(1.,2.*float(maptot)/float(npair)/float(mapmax))**2
3550 endif
3551
3552 c...print
3553
3554 if(ish.ge.5)then
3555 write(ifch,*)'print'
3556 do k=1,m3grid
3557 write(ifch,*)'k=',k,' jclu=',jclu(k)
3558 & ,' s=',(k-1)*delsgr-sgrid-shifts
3559 do j=m1grid,1,-1
3560 write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,m1grid)
3561 & ,(jdropgrid(i,j,k),i=1,m1grid)
3562 enddo
3563 enddo
3564 write(ifch,'(a,a)')
3565 & ' k jj nseg mm mseg n mass'
3566 & ,' s y z t '
3567 do jj=1,jjj
3568 do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
3569 mm=mmji(jj,ii)
3570 n=nptlb+mm
3571 sg=pptl(3,n)/abs(pptl(3,n))
3572 tm=sqrt(pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2)
3573 y=sg*alog((pptl(4,n)+sg*pptl(3,n))/tm)
3574 c if(kclu(jj).eq.44)print *,tm,pptl(4,n),pptl(3,n),iorptl(n)
3575 write(ifch,'(2i5,i6,i8,2i6,5f10.3)')
3576 & kclu(jj),jj,nseg(jj),mm,mseg(mm),n,pptl(5,n)
3577 & ,sptl(n),y,xorptl(3,n),xorptl(4,n)
3578 enddo
3579 enddo
3580 endif
3581
3582 c...decay
3583 iret=0
3584 if(jjj.gt.0)then !decay only if some cluster produced
3585 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)
3586 if(ish.ge.5)write(ifch,*)'decay ...'
3587 if(ifrade.eq.0.or.ispherio.gt.0)goto1000
3588 if(jdecay.eq.0)goto1000
3589 nptlbcf=nptl
3590 nptl0=nptl
3591 if(hydt.ne.'---')then
3592 call HydroFO(ier)
3593 else
3594 nclu=0
3595 ptest(1)=0d0
3596 ptest(2)=0d0
3597 ptest(3)=0d0
3598 ptest(4)=0d0
3599 ptest(5)=0d0
3600 do jj=1,jjj
3601 do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
3602 mm=mmji(jj,ii)
3603 np=nptlb+mm
3604 c print *,'decay',jj
3605 if(ioclude.eq.3)then
3606 call hnbaaa(np,iret)
3607 else
3608 call DropletDecay(np,iret)
3609 endif
3610 if(iret.eq.1)then
3611 istptl(np)=istptl(np)+2
3612 do ns=naseg(mm-1)+1,naseg(mm)
3613 n=nsegmt(ns)
3614 if(mod(abs(idptl(n)),100).eq.99)then !restore lost energy
3615 istptl(n) = 5
3616 ior=iorptl(n)
3617 ifrptl(1,ior) = 0
3618 ifrptl(2,ior) = 0
3619 if(istptl(ior).eq.0)then
3620 do k=1,3
3621 pptl(k,ior)=pptl(k,ior)+pptl(k,n)
3622 enddo
3623 pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2
3624 & +pptl(3,ior)**2+pptl(5,ior)**2)
3625 call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm)
3626 tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen()))
3627 else
3628 istptl(ior) = 0
3629 endif
3630 elseif(idptl(n).lt.1e4)then
3631 istptl(n) = 0 !particle
3632 ifrptl(1,n) = 0
3633 ifrptl(2,n) = 0
3634 else
3635 istptl(n) = 10 !droplet
3636 endif
3637 enddo
3638 elseif(ioclude.eq.3)then
3639 do i=nptl0+1,nptl
3640 if(ityptl(i).eq.60)then
3641 nclu=nclu+1
3642 ptest(1)=ptest(1)+dble(pptl(1,i))
3643 ptest(2)=ptest(2)+dble(pptl(2,i))
3644 ptest(3)=ptest(3)+dble(pptl(3,i))
3645 ptest(4)=ptest(4)+dble(pptl(4,i))
3646 endif
3647 enddo
3648 nptl0=nptl
3649 endif
3650 enddo
3651 enddo
3652 endif
3653 do jj=1,jjj
3654 do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
3655 mm=mmji(jj,ii)
3656 np=nptlb+mm
3657 istptl(np)=istptl(np)+1
3658 ifrptl(1,np)=nptlbcf+1
3659 ifrptl(2,np)=nptl
3660 rinptl(np)=kclu(jj)-m3grid/2
3661 enddo
3662 enddo
3663 c add global flow on all particles of all decayed cluster
3664 if(iLHC.eq.1)then
3665 yrmax=0.
3666 if(fvisco.gt.0.)yrmax=yradpx*visco
3667 else
3668 yrmax=yradpx
3669 endif
3670 c print *,bimevt,rini,yrmax,yrmaxi,delzet,np,ptmax,visco
3671 if(ioclude.eq.3.and.yrmax.gt.1e-3.and.nclu.gt.0)then
3672 c set angular informations
3673 fecc=0
3674 aa=1
3675 bb=0
3676 cc=0
3677 dd=1
3678 dta=0.5*abs(xx-yy)
3679 ev1=(xx+yy)/2+sqrt(dta**2+xy**2)
3680 ev2=(xx+yy)/2-sqrt(dta**2+xy**2)
3681 ecc=(ev1-ev2)/(ev1+ev2)
3682 c fecc=facecc*ecc!/(1.+yrmax)
3683 c print*,'pp',ecc,ranphi
3684 fecc=min(facecc,ecc) !be careful : fecc change <pt> since it is the elliptical deformation of the sub cluster (give strength of v2)
3685
3686 phiclu=mod(phievt+ranphi,2.*pi) !do not change otherwise v2 is gone
3687 if(phiclu.lt.-pi)phiclu=phiclu+2*pi
3688 if(phiclu.gt.pi)phiclu=phiclu-2*pi
3689 aa=cos(phiclu)
3690 bb=sin(phiclu)
3691 cc=-sin(phiclu)
3692 dd=cos(phiclu)
3693 errlim=0.00005
3694 c loop on particles for each main cluster
3695 c ncl=nptlb0
3696 npass=max(1,min(nclu/5,jjj)) !to have the same number of group of particles than original clusters but different repartition of particles
3697 npart=nclu/npass
3698 if(npart*npass-nclu.gt.max(5,npart/2))npass=npass+1
3699 c print *,'ici',nclu,npass,npart
3700 c lcont=.true.
3701 c if(nclu.lt.50)lcont=.false.
3702 lcont=.false.
3703 ncl=0
3704 nmin=nptlbcf+1
3705 nmax=nptl
3706 idrc=-1+2.*int(0.5+rangen())
3707 ntot=nclu
3708
3709 c prepare debug output for flow
3710 nall=0
3711 if(ish.ge.5)then
3712 nall=nmax-nmin+1
3713 do ii=1,nall
3714 idptl(nptl+ii)=idptl(nptlbcf+ii)
3715 enddo
3716 iorptl(nptl+nall+1)=nptl+1
3717 jorptl(nptl+nall+1)=nptl+nall
3718 do k=1,5
3719 pptl(k,nptl+nall+1)=0
3720 do ii=1,nall
3721 pptl(k,nptl+ii)=pptl(k,nptlbcf+ii)
3722 pptl(k,nptl+nall+1)=pptl(k,nptl+nall+1)+pptl(k,nptlbcf+ii)
3723 enddo
3724 enddo
3725 endif
3726
3727 c initialization for rescaling on ipart particles
3728 ipart=0
3729 tecm0=0.
3730 tecm=0.
3731 ini0=0
3732 ifi0=0
3733
3734 c do 900 while (ncl.le.nptlb-1)
3735 do while (ntot.gt.0)
3736 ncl=ncl+1
3737 if(iLHC.eq.1)then
3738 yrmax=yradpx*visco
3739 else
3740 yrmax=yradpx
3741 endif
3742 c cms frame of all particles from same cluster
3743 do k=1,5
3744 ppp(k)=0d0
3745 enddo
3746 ini=nptl
3747 ifi=1
3748 if(idrc.gt.0)then
3749 imax=nmax
3750 imin=nmin-1
3751 if(ipart.eq.0)ini0=nmin
3752 else
3753 imax=nmin
3754 imin=nmax+1
3755 if(ipart.eq.0)ifi0=nmax
3756 endif
3757 c npclu=0
3758 c 880 ncl=ncl+1
3759 n=0
3760 i=imin
3761 lpass=.true.
3762 do while ((n.lt.npart.or.ncl.eq.npass)
3763 & .and.idrc*i.lt.idrc*imax.and.lpass)
3764 i=i+idrc
3765 c if(jorptl(i).eq.ncl)then
3766 c if(jorptl(i).eq.ncl.and.ityptl(i).eq.60)then
3767
3768 if(ityptl(i).eq.60)then
3769 n=n+1
3770 ntot=ntot-1
3771 ini=min(ini,i)
3772 ifi=max(ifi,i)
3773 c if(ityptl(i).eq.60)npclu=npclu+1
3774 do k=1,4
3775 ppp(k)=ppp(k)+dble(pptl(k,i))
3776 enddo
3777 elseif(.not.lcont.and.n.gt.0)then
3778 lpass=.false.
3779 endif
3780 c print *,ityptl(i),i,imin,imax,idrc,n,npart,ncl,nclu,nptl
3781 enddo
3782 np=n
3783 if(idrc.gt.0)then
3784 nmin=nmin+i-imin
3785 ifi0=i
3786 else
3787 nmax=nmax+i-imin
3788 ini0=i
3789 endif
3790
3791 c if(ncl.lt.nptlb.and.npclu.lt.int(0.2*(nptl-nptlbcf+1)))goto880
3792 if(ifi.le.ini)goto 900
3793
3794 c record info for rescaling
3795 ipart=ipart+np
3796
3797 c test mass
3798 ppp(5)=(ppp(4)-ppp(3))*(ppp(4)+ppp(3))-ppp(2)**2-ppp(1)**2
3799 if(ppp(5).gt.0d0)then
3800 ppp(5)=sqrt(ppp(5))
3801 else
3802 if(ish.ge.1)write(ifch,*)'Precision problem in jintpo, p:',
3803 & (ppp(k),k=1,5)
3804 ppp(5)=0d0
3805 endif
3806 if(ish.ge.4)
3807 & write(ifch,*)'Group of particle: ',
3808 & idrc,ini,ifi,ncl,'/',npass,npart,nclu
3809 if(ppp(5).gt.0d0)then ! here all particle should have flow
3810 do i=ini,ifi
3811 if(ityptl(i).eq.60)then
3812 tecm=tecm+pptl(4,i) !use energy in collision frame as reference
3813 call utlob4(1,ppp(1),ppp(2),ppp(3),ppp(4),ppp(5)
3814 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
3815 endif
3816 enddo
3817 if(tecm.gt.0.)then
3818 yrmax=yrmax*log(amctot) !/rini**yradpi
3819 else
3820 yrmax=0.
3821 endif
3822 c print *,bimevt,rini,ppp(5),yrmax,yrmaxi,delzet,np,ptmax,visco
3823 if(ish.ge.4)
3824 & write(ifch,*)'Radial flow: ',yrmax,tecm,visco,yradpx,yradpp
3825 if(yrmax.gt.0.)then
3826 if(np.gt.maxp)stop'maxp too small in jintpo'
3827 i=0
3828 if(ish.ge.8)call clist('list before flow&',ini,ifi,60,60)
3829 do ii=ini,ifi
3830 if(ityptl(ii).eq.60)then
3831 i=i+1
3832 yrad(i)=sqrt(rangen())
3833 phirad(i)=2.*pi*rangen()
3834 pt2=(pptl(1,ii)**2+pptl(2,ii)**2)!+pptl(5,ii)**2)
3835 bex=sinh(yrad(i)*yrmax)*cos(phirad(i))
3836 & *(1+fecc/(1.+pt2))
3837 bey=sinh(yrad(i)*yrmax)*sin(phirad(i))
3838 & *(1-fecc/(1.+pt2))
3839 be(1)=aa*bex+cc*bey
3840 be(2)=bb*bex+dd*bey
3841 be(3)=0d0
3842 be(4)=sqrt(1+be(1)**2+be(2)**2)
3843 c call utlob4(1,be(1),be(2),be(3),be(4),1d0
3844 c * , pptl(1,ii), pptl(2,ii), pptl(3,ii), pptl(4,ii))
3845 c mimic boost transformation but protect against to high values of p(3) (p(3)~p(4))
3846 pt2=pptl(1,ii)**2+pptl(2,ii)**2
3847 bet=-(pptl(1,ii)*be(1)+pptl(2,ii)*be(2))/(be(4)+1.)
3848
3849 pt=yradpp**max(1.,pptl(5,ii))
3850 fac=1./(1.+sqrt(pt2/pptl(5,ii)**2))**pt
3851 bet=bet+sqrt(pt2+(pptl(3,ii)**2+pptl(5,ii)**2)*fac)
3852
3853 pptl(1,ii)=pptl(1,ii)-bet*be(1)
3854 pptl(2,ii)=pptl(2,ii)-bet*be(2)
3855 pptl(4,ii)=sqrt(pptl(1,ii)**2+pptl(2,ii)**2
3856 * +pptl(3,ii)**2+pptl(5,ii)**2)
3857 else
3858 yrad(i)=0.
3859 phirad(i)=0.
3860 endif
3861 enddo
3862 if(ish.ge.8)call clist('list after flow&',ini,ifi,60,60)
3863
3864
3865 c boost back
3866 pe(1)=0.
3867 pe(2)=0.
3868 pe(4)=0.
3869 do i=ini,ifi
3870 if(ityptl(i).eq.60)then
3871 call utlob4(-1,ppp(1),ppp(2),ppp(3),ppp(4),ppp(5)
3872 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
3873 pe(1)=pe(1)+pptl(1,i)
3874 pe(2)=pe(2)+pptl(2,i)
3875 pe(4)=pe(4)+pptl(4,i)
3876 endif
3877 enddo
3878
3879
3880
3881 cc random rescaling
3882 c rescaling has to be done for different bins in eta to conserve
3883 c energy flow (if everything is done at the end on all particles,
3884 c the energy flow is concentred at mid-rapidity in contradiction
3885 c with ATLAS data
3886 c On the other hand, a sufficient number of particles is necessary
3887 c to have proper eta distribution in particular at high pt.
3888 c Since the clusters are ordered in rapidity we can do the rescaling on
3889 c group of clusters having enough particles for the rescaling but
3890 c not too much to keep the eta dependence of the energy flow.
3891
3892 c set ipart.ge.1 to do rescaling for each subclusters
3893
3894 if((ipart.ge.1.and.ntot.ge.ipart/2).or.ntot.eq.0)then
3895
3896 ini=ini0
3897 ifi=ifi0
3898
3899 if(ish.ge.8)
3900 & call clist('list before flow rescaling&',ini,ifi,60,60)
3901
3902 if(fplmin.gt.0)then
3903
3904 c ntrt=ntrt+1
3905
3906 niter=0
3907 611 energ=0.
3908 energ0=0.
3909 niter=niter+1
3910 i=0
3911 ptmax=0.
3912 plmax=0.
3913 do j=ini,ifi
3914 if(ityptl(j).eq.60)then
3915 i=i+1
3916 pt2=pptl(1,j)**2+pptl(2,j)**2
3917 pz2=pptl(3,j)**2
3918 pp2=pt2+pz2
3919 c et2=(pp2+pptl(5,j)**2)*pt2/pp2
3920 pt=sqrt(pt2)
3921 c pp=sqrt(pp2)
3922 c base necessary to avoid peak at pt or pl=0
3923 c epsi change the shape of eta distributions and pt for
3924 c identified particles (shift the maximum of the distribution)
3925 c epsi=min(0.99,fplmin*rangen()**0.3)
3926 base=0. !sqrt(1./(1.-epsi)**2-1.)*pptl(5,j)/pp
3927 finc=2.
3928 if(energ0.lt.tecm)then
3929 finc=finc*sqrt(float(max(1,niter-300)))
3930 c base=base+min(tecm/pe(4),
3931 c & (max(0,niter-300))*0.3*pptl(5,j)/pp)
3932 else
3933 c base=base/log10(max(10.,float(niter-300)))
3934 finc=finc/sqrt(max(1.,float(niter-300)))
3935 endif
3936 c if(1.-pptl(5,j)/pptl(4,j).ge.epsi
3937 c & .and.1.+pptl(5,j)**2*(1./pp2-1./pt2).gt.0.)then
3938 if(1.+pptl(5,j)**2*(1./pp2-1./pt2).gt.0.)then
3939 ptmax=max(pt,ptmax)
3940 plmax=max(abs(pptl(3,j)),plmax)
3941 yrad2(i)=rangen()
3942 yrad2(i)=yrad2(i)*
3943 & max(0.,min(1.,finc*tecm/pe(4))-base)
3944
3945
3946 c necessary even with epsi cut to smooth eta distribution
3947 c and since sumEt has to be conserved, we can constrain scaling for pt
3948 c and pz :
3949 c ytmp=yrad2(i)
3950 yrad2(i)= yrad2(i)**((1.-(pptl(5,j)
3951 & /sqrt((min(1.,yrad2(i)*float(max(1,niter-100))))**2
3952 & *(pt**2+pptl(3,j)**2)+pptl(5,j)**2))**1.)
3953 & *exp(-fplmin*max(0.,pt**2-(pe(4)/tecm)**2)))
3954
3955 yrad2(i)=min(1.,base+yrad2(i)) !should be here to avoid peak at eta=0.
3956 else
3957 yrad2(i)=1.
3958 c ytmp=yrad2(i)
3959 endif
3960 be(1)= pptl(1,j)*yrad2(i)
3961 be(2)= pptl(2,j)*yrad2(i)
3962 be(3)= pptl(3,j)*yrad2(i)
3963
3964 c print *,niter,i,ntry,yrad2(i),energ/tecm,pt,pptl(3,j),finc,base
3965 c * ,ytmp,(1.-(pptl(5,j)
3966 c & /sqrt((min(1.,yrad2(i)*float(max(1,niter-100))))**2
3967 c & *(pt**2+pptl(3,j)**2)+pptl(5,j)**2))**0.1)
3968
3969 energ=energ+sqrt(be(1)**2+be(2)**2
3970 * +be(3)**2+pptl(5,j)**2)
3971 endif
3972 enddo
3973 energ0=energ
3974 c print *,'fin',niter,energ/tecm
3975 if(abs(energ-tecm)/tecm.gt.1..and.niter.lt.1000)then
3976 goto 611
3977 elseif(niter.ge.1000)then
3978 c print *,'Rescaling failed:',pe(4),tecm,ptmax,plmax,ipart,ini0,ifi
3979 if(ish.ge.2)write(ifch,*)'Random rescaling failed:'
3980 & ,energ,tecm,ptmax,plmax,ipart,ini0,ifi
3981 c ntrr=ntrr+1
3982 goto 200
3983 endif
3984 c print *,'done',niter,energ/tecm,ipart,ini0,ifi,finc
3985 if(ish.ge.5)write(ifch,*)'Rescaling done:'
3986 & ,tecm,energ/tecm,niter,ptmax,plmax,ipart,ini0,ifi
3987 i=0
3988 do j=ini,ifi
3989 if(ityptl(j).eq.60)then
3990 i=i+1
3991 pptl(1,j)= pptl(1,j)*yrad2(i)
3992 pptl(2,j)= pptl(2,j)*yrad2(i)
3993 pptl(3,j)= pptl(3,j)*yrad2(i)
3994 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
3995 * +pptl(3,j)**2+pptl(5,j)**2)
3996 endif
3997 enddo
3998 endif
3999
4000 200 continue
4001
4002
4003 if(ish.ge.8)
4004 & call clist('list after flow rescaling&',ini,ifi,60,60)
4005
4006 tecm0=tecm0+tecm
4007 tecm=0.
4008 ipart=0
4009
4010 endif
4011
4012 endif ! yrmax
4013
4014 endif !p(5)>0
4015
4016
4017 900 continue
4018 enddo
4019
4020 c rescale momentum precisely and globaly to avoid artefacts for matrix
4021 c but only at 10% level to keep dependence of sumEt with eta
4022 if(tecm0.gt.0)then
4023 ini=nptlbcf+1
4024 ifi=nptl
4025 esoll=tecm0
4026 scal=1.
4027 do ipass=1,2000
4028 sum=0.
4029 n=0
4030 do j=ini,ifi
4031 if(ityptl(j).eq.60)then
4032 n=n+1
4033 c this part is EXTREMELY important for the pseudorapidity shape at various pt
4034 c if nothing special a broad peak appear at eta=0
4035 c to avoid that, scal has to be reduced when p3 or pt reach 0
4036 if(scal.lt.1.)then
4037 scal0=scal
4038 pt=sqrt(pptl(1,j)**2+pptl(2,j)**2)
4039 if(fplmin.le.0.)then
4040 pt=abs(fplmin)/sqrt(pptl(5,j))*pt
4041 c pt=abs(fplmin)*pt
4042 else
4043 pt=5.*pt
4044 endif
4045 pow=sqrt(pptl(5,j)**2+scal**2*(pt**2
4046 * +pptl(3,j)**2))
4047 pow=(1.-(1./sqrt(pptl(5,j))+pptl(5,j))
4048 * /(1./sqrt(pptl(5,j))+pow))
4049 pow=rangen()*pow
4050 else
4051 pow=1.-2.*pptl(4,j)/engy !to avoid particle with energy larger than beam energy
4052 if(pow.lt.0.)then
4053 scal0=1./((1.-pow)
4054 * *exp(-0.25*max(-4.,log(rangen()))))
4055 pow=1.
4056 c print *,j,pptl(4,j),pow,scal0,scal
4057 c * ,pptl(3,j)*scal0**pow
4058 else
4059 scal0=scal
4060 pow=rangen()*pow
4061 endif
4062 endif
4063
4064 do k=1,3
4065 pptl(k,j)=pptl(k,j)*scal0**pow !to smooth distributions
4066 enddo
4067 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
4068 * +pptl(3,j)**2+pptl(5,j)**2)
4069 sum=sum+pptl(4,j)
4070 endif
4071 enddo
4072 scal=esoll/sum
4073 c write(ifmt,*)'ipass,scal,e,esoll:'
4074 c $ ,ipass,scal,sum,esoll
4075 if(abs(scal-1.).le.errlim) goto 300
4076 enddo
4077 300 continue
4078 c write(ifmt,*)'ipass,scal,e,esoll:'
4079 c $ ,ipass,scal,sum,esoll
4080
4081 c adjust pt to have pt conservation in cms of particles having flow
4082 if(nclu.gt.0)then
4083 ptest(5)=(ptest(4)-ptest(3))*(ptest(4)+ptest(3))-ptest(2)**2
4084 & -ptest(1)**2
4085 if(ptest(5).gt.0d0)then
4086 ptest(5)=sqrt(ptest(5))
4087 else
4088 if(ish.ge.1)write(ifch,*)'Precision problem in jintpo, p:',
4089 & (ptest(k),k=1,5)
4090 ptest(5)=0d0
4091 endif
4092 be(1)=0.d0
4093 be(2)=0.d0
4094 do i=ini,ifi
4095 if(ityptl(i).eq.60)then
4096 call utlob4(1,ptest(1),ptest(2),ptest(3),ptest(4),ptest(5)
4097 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
4098 be(1)=be(1)+dble(pptl(1,i))
4099 be(2)=be(2)+dble(pptl(2,i))
4100 endif
4101 enddo
4102 c shift nclu particles to have sum_pt=0. and boost back in global cms
4103 pt1shift=-sngl(be(1)/dble(nclu))
4104 pt2shift=-sngl(be(2)/dble(nclu))
4105 do i=ini,ifi
4106 if(ityptl(i).eq.60)then
4107 pptl(1,i)=pptl(1,i)+pt1shift
4108 pptl(2,i)=pptl(2,i)+pt2shift
4109 pptl(4,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2
4110 & +pptl(3,i)**2+pptl(5,i)**2)
4111 call utlob4(-1,ptest(1),ptest(2),ptest(3),ptest(4),ptest(5)
4112 & ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
4113 endif
4114 enddo
4115 endif
4116
4117 endif
4118
4119
4120 c if(ntrt.gt.0)print *,"jintpo rescaling :",float(ntrr)/float(ntrt)
4121 c define life time
4122 n=0
4123 do i=ini,ifi
4124 if(ityptl(i).eq.60)then
4125 n=n+1
4126 r=1.15*rini*yrad(n) !yrad=y/ymax
4127 tau=2.25/sqrt(yrad(n)**2+0.04)-0.75
4128 z=xorptl(3,i)
4129 t=xorptl(4,i)
4130 ! zeta=0.5*log((t+z)/(t-z))-0.5*delzet+2*0.5*delzet*rangen()
4131 test=(pptl(4,i)-pptl(3,i))*(pptl(4,i)+pptl(3,i))
4132 if(test.gt.0.)then
4133 zeta=0.5*log((pptl(4,i)+pptl(3,i))
4134 & /(pptl(4,i)-pptl(3,i)))
4135 else !in case of precision problem (but not always good neither for p<0
4136 pt=sqrt(pptl(2,i)**2+pptl(1,i)**2)
4137 zeta=0.5*log(1+2*pptl(3,i)*(pptl(4,i)+pptl(3,i))
4138 & /(pt*pt+pptl(5,i)**2))
4139 endif
4140 z=tau*sinh(zeta)
4141 t=tau*cosh(zeta)
4142 xorptl(1,i)=xorptl(1,i)+r*cos(phirad(n))
4143 xorptl(2,i)=xorptl(2,i)+r*sin(phirad(n))
4144 xorptl(3,i)=z
4145 xorptl(4,i)=t
4146 endif
4147 enddo
4148
4149
4150 if(ish.ge.5)then
4151 do k=1,5
4152 pptl(k,nptl+nall+2)=0
4153 do ii=nptlbcf+1,nptl
4154 pptl(k,nptl+nall+2)=pptl(k,nptl+nall+2)+pptl(k,ii)
4155 enddo
4156 enddo
4157 iorptl(nptl+nall+2)=nptlbcf+1
4158 jorptl(nptl+nall+2)=nptl
4159 call alist2('longitudinal and radial flow&',nptl+1
4160 & ,nptl+nall,nptlbcf+1,nptl)
4161 call alist2('momentum sum&',nptl+nall+1,nptl+nall+1
4162 & ,nptl+nall+2,nptl+nall+2)
4163 write(ifch,'(1x,50a1/)')('-',k=1,50)
4164 endif
4165
4166
4167 endif !ioclude=3 and flow
4168
4169
4170 do n=nptlbcf+1,nptl
4171 if(ioclude.ne.3)then
4172 iorptl(n)=nptlb+1
4173 jorptl(n)=nptlbcf
4174 rinptl(n)=rinptl((iorptl(n)+jorptl(n))/2)
4175 else
4176 rinptl(n)=rinptl(iorptl(n))
4177 endif
4178 istptl(n)=0
4179 ifrptl(1,n)=0
4180 ifrptl(2,n)=0
4181 tivptl(1,n)=xorptl(4,n)
4182 call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
4183 r=rangen()
4184 tivptl(2,n)=tivptl(1,n)+taugm*(-alog(r))
4185 radptl(n)=0.
4186 dezptl(n)=0.
4187 itsptl(n)=0
4188 enddo
4189
4190 endif
4191
4192 c Decay droplets not included in clusters
4193 iret=0
4194 do mm=1,nptla
4195 nptlb=nptl
4196 if(istptl(mm).eq.10)then
4197 if(ish.ge.5)write(ifch,*)'Decay remaining droplet :',mm
4198 if(nptlb.gt.mxptl-10)
4199 & call utstop('jintpo: mxptl too small (2)&')
4200 if(ioclude.eq.3)then
4201 call hnbaaa(mm,iret)
4202 else
4203 call DropletDecay(mm,iret) !Decay remn
4204 iret=0
4205 endif
4206 if(iret.eq.0.and.nptl.ne.nptlb)then ! ---successful decay---
4207 istptl(mm)=istptl(mm)+1
4208 ifrptl(1,mm)=nptlb+1
4209 ifrptl(2,mm)=nptl
4210 t=tivptl(2,mm)
4211 x=xorptl(1,mm)+(t-xorptl(4,mm))*pptl(1,mm)/pptl(4,mm)
4212 y=xorptl(2,mm)+(t-xorptl(4,mm))*pptl(2,mm)/pptl(4,mm)
4213 z=xorptl(3,mm)+(t-xorptl(4,mm))*pptl(3,mm)/pptl(4,mm)
4214 do 21 n=nptlb+1,nptl
4215 iorptl(n)=mm
4216 jorptl(n)=0
4217 istptl(n)=0
4218 ifrptl(1,n)=0
4219 ifrptl(2,n)=0
4220 radius=0.8*sqrt(rangen())
4221 phi=2*pi*rangen()
4222 ti=t
4223 zi=z
4224 xorptl(1,n)=x + radius*cos(phi)
4225 xorptl(2,n)=y + radius*sin(phi)
4226 xorptl(3,n)=zi
4227 xorptl(4,n)=ti
4228 iioo=mm
4229 zor=dble(xorptl(3,iioo))
4230 tor=dble(xorptl(4,iioo))
4231 r=rangen()
4232 tauran=-taurea*alog(r)
4233 call jtaix(n,tauran,zor,tor,zis,tis)
4234 tivptl(1,n)=amax1(ti,tis)
4235 call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
4236 r=rangen()
4237 tivptl(2,n)=t+taugm*(-alog(r))
4238 radptl(n)=0.
4239 dezptl(n)=0.
4240 itsptl(n)=0
4241 rinptl(nptl)=-9999
4242 21 continue
4243 else ! Unsuccessful decay
4244 if(ish.ge.1)write(ifch,*)
4245 * '***** Unsuccessful remnant cluster decay'
4246 * ,' --> redo event.'
4247 endif
4248 endif
4249 enddo
4250
4251
4252
4253 1000 continue
4254 call utprix('jintpo',ish,ishini,4)
4255 end
4256
4257 cc-----------------------------------------------------------------------
4258 c subroutine jrad(i,nq,na,jc,rad)
4259 cc-----------------------------------------------------------------------
4260 cc return hadron radius (data taken from huefner and povh)
4261 cc-----------------------------------------------------------------------
4262 c include 'epos.inc'
4263 c integer jc(nflav,2),kc(nflav)
4264 c
4265 c id=iabs(idptl(i))
4266 c am=pptl(5,i)
4267 c if(id.lt.10000)then
4268 c k=mod(id,10)
4269 c else
4270 c k=1
4271 c endif
4272 c do l=1,nflav
4273 c kc(l)=iabs(jc(l,1)-jc(l,2))
4274 c enddo
4275 c
4276 c if(nq.eq.0)then ! mesons
4277 c if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0.and.kc(4).eq.0)then
4278 c if(k.eq.0)then ! flavor singlet pseudoscalar mesons
4279 c if(am.ge.0.000)then
4280 c rad=0.64 ! pi0
4281 c if(am.ge.0.500)then
4282 c rad=0.60 ! eta
4283 c if(am.ge.0.900)then
4284 c rad=0.40 ! eta prime
4285 c if(am.ge.2.900)then
4286 c rad=0.17 ! eta charm
4287 c endif
4288 c endif
4289 c endif
4290 c else
4291 c write(ifch,*)
4292 c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am
4293 c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
4294 c call utstop('jrad: meson radius not defined&')
4295 c endif
4296 c else ! flavor singlet vector mesons
4297 c if(am.ge.0.000)then
4298 c rad=0.72 ! rho,omega
4299 c if(am.ge.1.000)then
4300 c rad=0.46 ! phi
4301 c if(am.ge.3.000)then
4302 c rad=0.20 ! J/psi
4303 c endif
4304 c endif
4305 c else
4306 c write(ifch,*)
4307 c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am
4308 c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
4309 c call utstop('jrad: meson radius not defined&')
4310 c endif
4311 c endif
4312 c elseif(kc(3).eq.0.and.kc(4).eq.0)then ! nonstrange, noncharmed
4313 c if(k.eq.0)then
4314 c rad=0.64 ! pi
4315 c else
4316 c rad=0.72 ! resonances
4317 c endif
4318 c elseif(kc(3).ne.0.and.kc(4).eq.0)then ! strange
4319 c if(k.eq.0)then
4320 c rad=0.59 ! kaons
4321 c else
4322 c rad=0.68 ! kaon resonances
4323 c endif
4324 c else ! charmed
4325 c write(ifch,*)'i:',i,' id:',idptl(i)
4326 c call utstop('jrad: radius of meson not defined&')
4327 c endif
4328 c else !baryons
4329 c if(kc(4).gt.0)then ! charmed
4330 c write(ifch,*)
4331 c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am
4332 c write(ifch,*)'i:',i,' id:',idptl(i)
4333 c call utstop('jrad: radius of charmed baryon not defined&')
4334 c elseif(kc(3).eq.0)then ! nonstrange
4335 c if(k.eq.0)then
4336 c rad=0.82 !nucleons
4337 c else
4338 c rad=1.00 !resonances
4339 c endif
4340 c elseif(kc(3).eq.1)then ! strange
4341 c if(k.eq.0)then
4342 c rad=0.76 !lambda, sigma
4343 c else
4344 c rad=0.93 !resonances
4345 c endif
4346 c elseif(kc(3).eq.2)then ! double strange
4347 c if(k.eq.0)then
4348 c rad=0.71 !cascades
4349 c else
4350 c rad=0.87 !resonances
4351 c endif
4352 c elseif(kc(3).ge.3)then ! triple strange
4353 c rad=0.79 !omega
4354 c else
4355 c write(ifch,*)
4356 c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am
4357 c write(ifch,*)
4358 c * 'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6),
4359 c * ' |q-qbar|:',(kc(l),l=1,6)
4360 c call utstop('jrad: should not happen&')
4361 c endif
4362 cc string fragments with |#q|>3
4363 c if(na.gt.3)then
4364 c a=(na/3.)**(1./3.)
4365 c if(ish.ge.7)then
4366 c call utmsg('jrad ')
4367 c write(ifch,*)
4368 c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am
4369 c write(ifch,*)
4370 c * 'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6),
4371 c * ' |q-qbar|:',(kc(l),l=1,6)
4372 c write(ifch,*)'nq:',nq,' na:',na,' r:',rad,' ar:',a*rad
4373 c call utmsgf
4374 c endif
4375 c rad=rad*a
4376 c endif
4377 c endif
4378 c
4379 c if(ish.ge.7)then
4380 c write(ifch,*)
4381 c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am,' rad:',rad
4382 c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
4383 c endif
4384 c
4385 c return
4386 c end
4387 c
4388 c-----------------------------------------------------------------------
4389 subroutine jresc
4390 c-----------------------------------------------------------------------
4391 include 'epos.inc'
4392 double precision pa(5),pj(5)
4393 integer ipptl(mxptl)
4394
4395 call utpri('jresc ',ish,ishini,4)
4396
4397 iret=0
4398 nptlpt=maproj+matarg
4399 np=0
4400 do i=nptlpt+1,nptl
4401 if(istptl(i).eq.0
4402 * .and.idptl(i).lt.10000.and.pptl(5,i).gt.0.01)then
4403 np=np+1
4404 ipptl(np)=i
4405 endif
4406 enddo
4407 if(np.lt.2)goto1001
4408 do ii=1,np
4409 i=ipptl(ii)
4410 if(mod(iabs(idptl(i)),10).lt.2)then
4411 call idmass(idptl(i),ami)
4412 dm=abs(ami-pptl(5,i))
4413 if(dm.gt.0.001)then
4414 ntry=0
4415 1 continue
4416 ntry=ntry+1
4417 2 jj=1+int(rangen()*np)
4418 j=ipptl(jj)
4419 if(ish.ge.4)write(ifch,*)i,j,istptl(j)
4420 if(j.eq.i)goto2
4421 if(mod(iabs(idptl(j)),10).lt.2)then
4422 call idmass(idptl(j),amj)
4423 else
4424 amj=pptl(5,j)
4425 endif
4426 do l=1,5
4427 pa(l)=dble(pptl(l,i))
4428 pj(l)=dble(pptl(l,j))
4429 enddo
4430 if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70)
4431 if(ish.ge.4)write(ifch,11)i,idptl(i),'before:',pa,'want:',ami
4432 if(ish.ge.4)write(ifch,11)j,idptl(j),'before:',pj,'want:',amj
4433 call jrescl(pa,dble(ami),pj,dble(amj),iret)
4434 if(iret.eq.1)then
4435 if(ntry.le.50)then
4436 goto1
4437 else
4438 goto1001
4439 endif
4440 endif
4441 if(ish.ge.4)write(ifch,11)i,idptl(i),' after:',pa
4442 if(ish.ge.4)write(ifch,11)j,idptl(j),' after:',pj
4443 if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70)
4444 do l=1,5
4445 pptl(l,i)=sngl(pa(l))
4446 pptl(l,j)=sngl(pj(l))
4447 enddo
4448 endif
4449 endif
4450 enddo
4451 11 format(i5,1x,i5,1x,a,1x,5(d8.2,1x),a,1x,e8.2)
4452
4453 1000 continue
4454 call utprix('jresc ',ish,ishini,4)
4455 return
4456
4457 1001 continue
4458 if(ish.ge.1)then
4459 write(ifmt,'(a)')'jresc: could not put on shell'
4460 endif
4461 goto1000
4462
4463 end
4464
4465 c-----------------------------------------------------------------------
4466 subroutine jrescl(p1,am1,p2,am2,iret)
4467 c-----------------------------------------------------------------------
4468 c rescale momenta of two particles such that the masses assume given
4469 c values.
4470 c input:
4471 c p1, p2: momenta of the two particles
4472 c am1, am2: desired masses of the two particles
4473 c output:
4474 c p1, p2: new momenta of the two particles
4475 c-----------------------------------------------------------------------
4476 include 'epos.inc'
4477 double precision p1(5),p2(5)
4478 * ,p1n(5),p2n(5)
4479 * ,a1,a2,a12,am1,am2
4480 * ,b1,b2,c,d,e,f,g,p,q,r
4481
4482 call utpri('jrescl',ish,ishini,7)
4483
4484 iret=0
4485 a1=p1(5)**2
4486 a2=p2(5)**2
4487 a12=p1(4)*p2(4)-p1(3)*p2(3)-p1(2)*p2(2)-p1(1)*p2(1)
4488 if(a12.le.(a1+a2))then
4489 if(ish.ge.7)write(ifch,*)'a_12 < a_1 + a_2'
4490 if(ish.ge.7)write(ifch,*)a12,' < ',a1+a2
4491 c goto1001
4492 endif
4493
4494 11 format(5(d9.3,1x))
4495 if(ish.ge.7)write(ifch,11)p1,a1
4496 if(ish.ge.7)write(ifch,11)p2,a2
4497 if(ish.ge.7)write(ifch,*)a12
4498
4499 c=(a1+a12)/(a2+a12)
4500 d=(a1-am1**2-a2+am2**2)/(a2+a12)*0.5d0
4501
4502 e=a1-2d0*a12*c+a2*c**2
4503 f=2d0*(a1-a12*(c+d)+a2*c*d)
4504 g=a1-2d0*a12*d+a2*d**2-am1**2
4505
4506 p=f/e
4507 q=g/e
4508 r=p**2-4d0*q
4509
4510 if(ish.ge.7)write(ifch,*)'c:',c,' d:',d
4511 if(ish.ge.7)write(ifch,*)'e:',e,' f:',f,' g:',g
4512 if(ish.ge.7)write(ifch,*)'p:',p,' q:',q,' r:',r
4513 if(r.lt.0d0)goto1001
4514
4515 b1=-0.5d0*(p-dsqrt(r))
4516
4517 b2=b1*c+d
4518
4519 if(ish.ge.7)write(ifch,*)'b_1:',b1,' b_2:',b2
4520
4521 do i=1,4
4522 p1n(i)=(1d0+b1)*p1(i)-b2*p2(i)
4523 p2n(i)=(1d0+b2)*p2(i)-b1*p1(i)
4524 enddo
4525
4526 a1=p1n(4)**2-p1n(3)**2-p1n(2)**2-p1n(1)**2
4527 a2=p2n(4)**2-p2n(3)**2-p2n(2)**2-p2n(1)**2
4528 if(a1.gt.0d0.and.a2.gt.0d0)then
4529 do i=1,4
4530 p1(i)=p1n(i)
4531 p2(i)=p2n(i)
4532 enddo
4533 p1(5)=dsqrt(a1)
4534 p2(5)=dsqrt(a2)
4535 if(ish.ge.7)write(ifch,11)p1,a1
4536 if(ish.ge.7)write(ifch,11)p2,a2
4537 else
4538 goto1001
4539 endif
4540
4541 if(p1(4).lt.0..or.p2(4).lt.0.)goto1001
4542
4543 1000 continue
4544 call utprix('jrescl',ish,ishini,7)
4545 return
4546
4547 1001 continue
4548 iret=1
4549 goto1000
4550 end
4551
4552 c-----------------------------------------------------------------------
4553 subroutine jtain(i,x,y,z,t,n,iopt)
4554 c-----------------------------------------------------------------------
4555 c returns intersection (x,y,z,t) of ptl-i-trajectory with taus-line.
4556 c input:
4557 c i: particle number
4558 c iopt: formation time considered (0) or not (1)
4559 c output:
4560 c x,y,z,t: 4-vector of intersection point
4561 c n: exit code
4562 c n=0: ok
4563 c n=1: ptl lives later
4564 c n=2: ptl lives earlier
4565 c n=9: tiv1>tiv2
4566 c-----------------------------------------------------------------------
4567 include 'epos.inc'
4568 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
4569 common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
4570 double precision vv,zza,zz,tt,xo3,xo4,ti1,ti2,derr,dd
4571 double precision ttp,zzp,ttt,zzt,vvt,vvp,spt2m2E,p4
4572 common/ctfi/tin,tfi
4573 double precision ttau0
4574 common/cttau0/ttau0
4575
4576 n=0
4577
4578 tin=0
4579 tfi=0
4580
4581 derr=1d-2
4582 ttp=tpro*ttaus
4583 zzp=zpro*ttaus
4584 ttt=ttar*ttaus
4585 zzt=ztar*ttaus
4586 vv=sign(min(1.d0,abs(dble(pptl(3,i)))/dble(pptl(4,i)))
4587 & ,dble(pptl(3,i)))
4588
4589
4590 if(abs(vv).ge.1.d0)then
4591 spt2m2E=dble(pptl(1,i)*pptl(1,i)+pptl(2,i)*pptl(2,i)
4592 & +pptl(5,i)*pptl(5,i))
4593 c if(pptl(4,i).le.0.)then
4594 p4=sqrt(dble(pptl(3,i)*pptl(3,i))+spt2m2E)
4595 c else
4596 c p4=dble(pptl(4,i))
4597 c endif
4598 ctp to avoid precision problem, replace abs(p3)/p4 by sqrt(1-(pt2+m2)/E2)
4599 spt2m2E=min(1.d0,sqrt(spt2m2E)/p4)
4600 vv=sign(sqrt((1d0+spt2m2E)*(1d0-spt2m2E)),dble(pptl(3,i)))
4601 endif
4602 xo3=dble(xorptl(3,i))
4603 xo4=dble(xorptl(4,i))
4604 zza=xo3-xo4*vv
4605 if(iopt.eq.0)then
4606 ti1=dble(tivptl(1,i))
4607 elseif(iopt.eq.1)then
4608 ti1=dble(xo4)
4609 else
4610 ti1=0
4611 call utstop("Wrong iopt in jtain !&")
4612 endif
4613 ti2=dble(tivptl(2,i))
4614
4615 if(ti1.gt.ti2)then
4616 n=9
4617 goto1
4618 endif
4619
4620 zfi=sngl(xo3+(ti2-xo4)*vv)
4621 call jtaus(zfi,tzfi,szfi)
4622 tfi=tzfi
4623 if(tfi.ge.sngl(ti2))then
4624 n=2
4625 goto1
4626 endif
4627
4628 zin=sngl(xo3+(ti1-xo4)*vv)
4629 call jtaus(zin,tzin,szin)
4630 tin=tzin
4631 if(tin.le.sngl(ti1))then
4632 n=1
4633 goto1
4634 endif
4635
4636
4637 1 continue
4638
4639 if(ttaus.le.ttau0)then
4640 tt=ttaus
4641 zz=xo3+(tt-xo4)*vv
4642 if(tt.lt.ti1.and.n.eq.0)n=1
4643 if(tt.ge.ti2.and.n.eq.0)n=2
4644 goto1000
4645 else
4646 vvt=zzt/ttt
4647 vvp=zzp/ttp
4648 tt=(ttt+(zza-zzt)*vvt)/(1-vv*vvt)
4649 zz=xo3+(tt-xo4)*vv
4650 if(zz.le.zzt)then
4651 if(tt.lt.ti1.and.n.eq.0)n=1
4652 if(tt.ge.ti2.and.n.eq.0)n=2
4653 goto1000
4654 endif
4655 tt=(ttp+(zza-zzp)*vvp)/(1-vv*vvp)
4656 zz=xo3+(tt-xo4)*vv
4657 if(zz.ge.zzp)then
4658 if(tt.lt.ti1.and.n.eq.0)n=1
4659 if(tt.ge.ti2.and.n.eq.0)n=2
4660 goto1000
4661 endif
4662 dd=1-vv**2
4663 if(sngl(dd).eq.0..and.vv.gt.0.)then
4664 tt=-(ttaus**2+zza**2)/2d0/zza
4665 elseif(sngl(dd).eq.0..and.vv.lt.0.)then
4666 tt=(ttaus**2+zza**2)/2d0/zza
4667 else
4668 tt=(zza*vv+dsqrt(zza**2+ttaus**2*dd))/dd
4669 endif
4670 zz=xo3+(tt-xo4)*vv
4671 if(tt.lt.ti1.and.n.eq.0)n=1
4672 if(tt.ge.ti2.and.n.eq.0)n=2
4673 if(dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr*ttaus**2.and.
4674 *dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr)then
4675 if(ish.ge.1)then
4676 call utmsg('jtain')
4677 write(ifch,*)'***** ttaus**2 .ne. (tt+zz)*(tt-zz)'
4678 write(ifch,*)sngl(ttaus**2),sngl((tt+zz)*(tt-zz))
4679 call utmsgf
4680 endif
4681 goto1000
4682 endif
4683 endif
4684
4685 1000 t=sngl(tt)
4686 z=sngl(zz)
4687 x=xorptl(1,i)+(t-xorptl(4,i))*pptl(1,i)/pptl(4,i)
4688 y=xorptl(2,i)+(t-xorptl(4,i))*pptl(2,i)/pptl(4,i)
4689 return
4690 end
4691
4692 c-----------------------------------------------------------------------
4693 subroutine jtaix(i,tau,zor,tor,z,t)
4694 c-----------------------------------------------------------------------
4695 c returns intersection z,t of ptl-i-trajectory with hyperbola h.
4696 c h: (t-tor)**2-(z-zor)**2=tau**2 .
4697 c zor, tor double precision.
4698 c-----------------------------------------------------------------------
4699 include 'epos.inc'
4700 double precision tor,zor,tors,zors,vv,cc,dd,ttau,derr,tt,zz
4701 derr=1d-3
4702 ttau=dble(tau)
4703 zors=dble(xorptl(3,i))-zor
4704 tors=dble(xorptl(4,i))-tor
4705 vv=dble(pptl(3,i)/pptl(4,i))
4706 vv=dmin1(vv,1d0)
4707 vv=dmax1(vv,-1d0)
4708 cc=zors-tors*vv
4709 dd=1d0-vv**2
4710 dd=dmax1(dd,0d0)
4711 if(dd.eq.0d0.and.cc.eq.0d0)then
4712 if(tau.eq.0.)tt=0d0
4713 if(tau.ne.0.)tt=dble(ainfin)
4714 zz=tt
4715 goto1000
4716 elseif(dd.eq.0d0)then
4717 tt=-(ttau**2+cc**2)/2d0/cc/vv
4718 elseif(dd.lt.1e-8)then
4719 tt=-(ttau**2+cc**2)/2d0/cc/vv
4720 call utmsg('jtaix')
4721 write(ifch,*)'***** dd = ',dd,' treated as zero'
4722 call utmsgf
4723 else
4724 tt=(cc*vv+dsqrt(cc**2+ttau**2*dd))
4725 tt=tt/dd
4726 endif
4727 zz=cc+tt*vv
4728 if(dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr*ttau**2
4729 *.and.dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr
4730 *.and.tors**2-zors**2.lt.1e6)then
4731 if(ish.ge.2)then
4732 call utmsg('jtaix')
4733 write(ifch,*)'***** ttau**2 .ne. (tt+zz)*(tt-zz)'
4734 write(ifch,*)sngl(ttau**2),sngl((tt+zz)*(tt-zz))
4735 write(ifch,*)'tau,t,z:'
4736 write(ifch,*)tau,tt,zz
4737 write(ifch,*)'#,id(ptl):',i,idptl(i)
4738 write(ifch,*)'zor,tor(str):',zor,tor
4739 write(ifch,*)'zors,tors,p,e(ptl):'
4740 write(ifch,*)sngl(zors),sngl(tors),pptl(3,i),pptl(4,i)
4741 call utmsgf
4742 endif
4743 endif
4744 1000 z=sngl(zz+zor)
4745 t=sngl(tt+tor)
4746 return
4747 end
4748
4749 c-----------------------------------------------------------------------
4750 subroutine jtaug(su,so,g,y)
4751 c-----------------------------------------------------------------------
4752 c returns g factor and rapidity y for given su, so
4753 c-----------------------------------------------------------------------
4754 include 'epos.inc'
4755 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
4756 common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
4757 double precision ttp,zzp,ttt,zzt,ssp,sst,ssu,sso,ss1,ss2,gg
4758 *,ssav,yyav,hh
4759 double precision ttau0
4760 common/cttau0/ttau0
4761
4762 ssu=dble(su)
4763 sso=dble(so)
4764
4765 if(ssu.ge.sso)then
4766 sso=(ssu+sso)*0.5d0 + dble(abs(dezzer))*ttaus*0.5d0
4767 ssu=(ssu+sso)*0.5d0 - dble(abs(dezzer))*ttaus*0.5d0
4768 so=real(sso)
4769 su=real(ssu)
4770 endif
4771 if(ssu.ge.sso)then
4772 print*,ssu,sso,dble(abs(dezzer))*ttaus*0.5d0
4773 stop'STOP: sr jtaug: ssu.ge.sso'
4774 endif
4775
4776 g=1
4777
4778 if(ttaus.le.ttau0)return
4779
4780 ttp=tpro*ttaus
4781 zzp=zpro*ttaus
4782 ttt=ttar*ttaus
4783 zzt=ztar*ttaus
4784 ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp))
4785 sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt))
4786
4787 ssav=(ssu+sso)/2d0
4788 yyav=ssav/ttaus
4789 if(ssav.ge.ssp)yyav=detap
4790 if(ssav.le.sst)yyav=detat
4791
4792 gg=0
4793 if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu)
4794 if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu))
4795 if(ssu.lt.ssp.and.sso.gt.sst)then
4796 ss1=dmax1(ssu,sst)
4797 ss2=dmin1(sso,ssp)
4798 gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) )
4799 endif
4800 gg=gg/(sso-ssu)
4801
4802 hh=0
4803 if(ssu.lt.sst)hh=hh + dsinh(detat-yyav) * (dmin1(sst,sso)-ssu)
4804 if(sso.gt.ssp)hh=hh + dsinh(detap-yyav) * (sso-dmax1(ssp,ssu))
4805 if(ssu.lt.ssp.and.sso.gt.sst)then
4806 ss1=dmax1(ssu,sst)
4807 ss2=dmin1(sso,ssp)
4808 hh=hh+ttaus*( dcosh(ss2/ttaus-yyav)-dcosh(ss1/ttaus-yyav) )
4809 endif
4810 hh=hh/(sso-ssu)
4811
4812 yyav=yyav+0.5d0*dlog((gg+hh)/(gg-hh))
4813
4814 gg=0
4815 if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu)
4816 if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu))
4817 if(ssu.lt.ssp.and.sso.gt.sst)then
4818 ss1=dmax1(ssu,sst)
4819 ss2=dmin1(sso,ssp)
4820 gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) )
4821 endif
4822 gg=gg/(sso-ssu)
4823
4824 g=sngl(gg)
4825 y=sngl(yyav)
4826 return
4827 end
4828
4829 c-----------------------------------------------------------------------
4830 subroutine jtaui(s,ts,zs)
4831 c-----------------------------------------------------------------------
4832 c returns time ts and coord zs corresponding to ttaus and inv. length s
4833 c-----------------------------------------------------------------------
4834
4835 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
4836 common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
4837 double precision ttau0
4838 common/cttau0/ttau0
4839
4840 double precision ttp,zzp,ttt,zzt,ssp,sst,ss,zeta
4841
4842 zs=s
4843 ts=sngl(ttaus)
4844
4845 if(ttaus.le.ttau0)return
4846
4847 ttp=tpro*ttaus
4848 zzp=zpro*ttaus
4849 ttt=ttar*ttaus
4850 zzt=ztar*ttaus
4851 ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp))
4852 sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt))
4853 ss=dble(s)
4854
4855 if(ss.le.sst)then
4856 zs=sngl(zzt+ttar*(ss-sst))
4857 ts=sngl(ttt+(dble(zs)-zzt)*zzt/ttt)
4858 elseif(ss.ge.ssp)then
4859 zs=sngl(zzp+tpro*(ss-ssp))
4860 ts=sngl(ttp+(dble(zs)-zzp)*zzp/ttp)
4861 else
4862 zeta=ss/ttaus
4863 ts=sngl(ttaus*dcosh(zeta))
4864 zs=sngl(ttaus*dsinh(zeta))
4865 endif
4866
4867 return
4868 end
4869
4870 c-----------------------------------------------------------------------
4871 subroutine jtauin
4872 c-----------------------------------------------------------------------
4873 c initializes equal time hyperbola at ttaus
4874 c-----------------------------------------------------------------------
4875 include 'epos.inc'
4876 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
4877 common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
4878 double precision ttau0,rcproj,rctarg
4879 common/geom1/rcproj,rctarg
4880 common/cttau0/ttau0
4881
4882 call utpri('jtauin',ish,ishini,6)
4883
4884 if(ttaus.gt.ttau0)then
4885 if(rcproj.gt.1d-10)then
4886 detap=dmin1(dble((ypjtl-yhaha)*etafac),dlog(ttaus/rcproj))
4887 else
4888 detap=dble((ypjtl-yhaha)*etafac)
4889 endif
4890 if(rctarg.gt.1d-10)then
4891 detat=dmax1(dble(-yhaha*etafac),dlog(rctarg/ttaus))
4892 else
4893 detat=dble(-yhaha*etafac)
4894 endif
4895 tpro=dcosh(detap)
4896 zpro=dsinh(detap)
4897 ttar=dcosh(detat)
4898 ztar=dsinh(detat)
4899 else
4900 detap=0d0
4901 detat=0d0
4902 tpro=0d0
4903 zpro=0d0
4904 ttar=0d0
4905 ztar=0d0
4906 endif
4907
4908 if(ish.ge.6)then
4909 write(ifch,*)'hyperbola at tau=',ttaus
4910 write(ifch,*)'r_p:',rcproj,' r_t:',rctarg
4911 write(ifch,*)'y_p:',detap,' y_t:',detat
4912 write(ifch,*)'t_p:',tpro,' z_p:',zpro
4913 write(ifch,*)'t_t:',ttar,' z_t:',ztar
4914 endif
4915
4916 call utprix('jtauin',ish,ishini,6)
4917 return
4918 end
4919
4920 c-----------------------------------------------------------------------
4921 subroutine jtaus(z,tz,sz)
4922 c-----------------------------------------------------------------------
4923 c returns time tz and inv length sz corresponding to ttaus and z
4924 c-----------------------------------------------------------------------
4925 include 'epos.inc'
4926 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
4927 common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
4928 double precision ttau0
4929 common/cttau0/ttau0
4930
4931 double precision ttp,zzp,ttt,zzt,zz,tzz
4932
4933 tz=sngl(ttaus)
4934 sz=z
4935
4936 if(ttaus.le.ttau0)return
4937
4938 ttp=tpro*ttaus
4939 zzp=zpro*ttaus
4940 ttt=ttar*ttaus
4941 zzt=ztar*ttaus
4942 zz=dble(z)
4943
4944 if(zz.le.zzt)then
4945 tz=sngl(ttt+(zz-zzt)*zzt/ttt)
4946 sz=sngl(ttaus*detat+(zz-zzt)/ttar)
4947 elseif(zz.ge.zzp)then
4948 tz=sngl(ttp+(zz-zzp)*zzp/ttp)
4949 sz=sngl(ttaus*detap+(zz-zzp)/tpro)
4950 else
4951 if(sngl(ttaus).ge.ainfin)then
4952 tz=sngl(ttaus)
4953 sz=0.
4954 if(ish.ge.1)then
4955 call utmsg('jtaus')
4956 write(ifch,*)'***** large ttaus; set tz=ttaus, sz=0'
4957 write(ifch,*)'ttaus=',ttaus,'zz=',zz
4958 call utmsgf
4959 endif
4960 else
4961 tzz=dsqrt(ttaus**2+zz**2)
4962 tz=sngl(tzz)
4963 sz=sngl(ttaus*0.5d0*dlog((tzz+zz)/(tzz-zz)))
4964 endif
4965 endif
4966
4967 return
4968 end
4969
4970 c-----------------------------------------------------------------------
4971 subroutine jtaux(t,z,ttaux)
4972 c-----------------------------------------------------------------------
4973 c returns ttaux (-> tau-line) for given t and z.
4974 c ttaux: double precision.
4975 c-----------------------------------------------------------------------
4976 double precision ttaux,tt,zz,rcproj,rctarg,zt1,zp1,zt2,zp2,ttau0
4977 common/geom1/rcproj,rctarg
4978 common/cttau0/ttau0
4979 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
4980 common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat
4981
4982 tt=dble(t)
4983 zz=dble(z)
4984
4985 if(tt.gt.ttau0)then
4986 zt1=rctarg-tt
4987 zp1=tt-rcproj
4988 zt2=ztar/ttar*tt
4989 zp2=zpro/tpro*tt
4990 if(zz.le.dmax1(zt1,zt2))then
4991 if(zt1.gt.zt2)then
4992 ttaux=rctarg*dsqrt((tt-zz)/(2d0*rctarg-tt-zz))
4993 else
4994 ttaux=(ttar*tt-ztar*zz)/(ttar**2-ztar**2)
4995 endif
4996 elseif(zz.ge.dmin1(zp1,zp2))then
4997 if(zp1.lt.zp2)then
4998 ttaux=rcproj*dsqrt((tt+zz)/(2d0*rcproj-tt+zz))
4999 else
5000 ttaux=(tpro*tt-zpro*zz)/(tpro**2-zpro**2)
5001 endif
5002 else
5003 ttaux=dsqrt(tt**2-zz**2)
5004 endif
5005 else
5006 ttaux=tt
5007 endif
5008
5009 return
5010 end
5011
5012 c-----------------------------------------------------------------------
5013 subroutine xjden1(ii,itau,x,y,rad,o,u)
5014 c-----------------------------------------------------------------------
5015 c ii=0: initialization
5016 c ii=1: determining dense regions in space of individual events
5017 c x,y,rad: tranverse coordinates and radius of particle i
5018 c o,u: specifies long range: u < s < o (s: long coordinate)
5019 c ii=2: plot of individual event
5020 c xpar1: itau ; valid: 1,..,10
5021 c xpar2: z-range: -xpar2 < z < xpar2
5022 c xpar3, x-range: -xpar3 < x < xpar3
5023 c xpar4, y-range: -xpar4 < y < xpar4
5024 c-----------------------------------------------------------------------
5025 include "epos.inc"
5026 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
5027 common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat
5028
5029 if(idensi.ne.1)stop'STOP in xjden1: idensi must be set 1'
5030
5031 dlcoox=0.5
5032 dlcooy=0.5
5033
5034 if(ii.eq.0)then
5035
5036 do i=1,nzeta
5037 do j=1,mxcoox
5038 do k=1,mxcooy
5039 kdensi(itau,i,j,k)=0
5040 enddo
5041 enddo
5042 enddo
5043
5044 elseif(ii.eq.1)then
5045
5046 if(itau.lt.1.or.itau.gt.mxtau)return
5047
5048 tau=sngl(ttaus)
5049 zu=u/tau
5050 zo=o/tau
5051
5052 do 1 i=1,nzeta
5053 zi=-nzeta*dlzeta/2-dlzeta/2+i*dlzeta
5054 if(zu.gt.zi.or.zo.lt.zi)goto1
5055 do 2 j=1,mxcoox
5056 xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox
5057 do 3 k=1,mxcooy
5058 yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
5059 if((x-xj)**2+(y-yk)**2.gt.rad**2)goto3
5060 kdensi(itau,i,j,k)=1
5061 3 continue
5062 2 continue
5063 1 continue
5064
5065 elseif(ii.eq.2)then
5066
5067 itaux=nint(xpar1)
5068 if(itaux.gt.mxtau)stop'STOP in xjden1: itaux too large'
5069
5070 iz=nint(xpar2/dlzeta)
5071 ix=nint(xpar3/dlcoox)
5072 iy=nint(xpar4/dlcooy)
5073 if(iz.gt.nzeta/2)stop'STOP in xjden1: zeta-range too large'
5074 if(ix.gt.mxcoox/2)stop'STOP in xjden1: x-range too large'
5075 if(iy.gt.mxcooy/2)stop'STOP in xjden1: y-range too large'
5076
5077 do k=mxcooy/2+1-iy,mxcooy/2+iy
5078 write(ifhi,'(a,f7.2)') '! tau: ',tauv(itaux)
5079 write(ifhi,'(a)') 'openhisto'
5080 write(ifhi,'(a,2f7.2)')'xrange',-iz*dlzeta,iz*dlzeta
5081 write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox
5082 write(ifhi,'(a)') 'set ityp2d 3'
5083 write(ifhi,'(a)') 'txt "xaxis space-time rapidity [z]"'
5084 write(ifhi,'(a)') 'txt "yaxis transverse coordinate x (fm)"'
5085 write(ifhi,'(a,i4)') 'array2d',2*iz
5086 do j=mxcoox/2+1-ix,mxcoox/2+ix
5087 write(ifhi,'(40i2)') (kdensi(itaux,i,j,k),
5088 * i=nzeta/2+1-iz,nzeta/2+iz)
5089 enddo
5090 write(ifhi,'(a)') ' endarray'
5091 write(ifhi,'(a)') 'closehisto plot2d'
5092 enddo
5093
5094 else
5095
5096 stop'STOP in xjden1: wrong option'
5097
5098 endif
5099
5100 return
5101 end
5102
5103 c-----------------------------------------------------------------------
5104 subroutine xjden2(ii,itau,x,y,rad,s)
5105 c-----------------------------------------------------------------------
5106 c ii=0: initialization
5107 c ii=1: determining dense regions in space of individual events
5108 c x,y,rad: tranverse coordinates and radius of particle i
5109 c s: long coordinate
5110 c ii=2: plot of individual event
5111 c xpar1: itau ; valid: 1,..,10
5112 c xpar2: s-range: -xpar2 < s < xpar2
5113 c xpar3, x-range: -xpar3 < x < xpar3
5114 c xpar4, y-range: -xpar4 < y < xpar4
5115 c-----------------------------------------------------------------------
5116 include "epos.inc"
5117 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
5118 common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat
5119 parameter (mxcoos=60)
5120 common/cdensh/kdensh(matau,mxcoos,mxcoox,mxcooy),ktot(matau)
5121 character cy*3
5122
5123 dlcoox=0.5
5124 dlcooy=0.5
5125 dlcoos=0.5
5126
5127 if(ii.eq.0)then
5128
5129 do i=1,mxcoos
5130 do j=1,mxcoox
5131 do k=1,mxcooy
5132 kdensh(itau,i,j,k)=0
5133 enddo
5134 enddo
5135 enddo
5136 ktot(itau)=0
5137
5138 elseif(ii.eq.1)then
5139
5140 if(itau.lt.1.or.itau.gt.mxtau)return
5141
5142 tau=sngl(ttaus)
5143 z=s/tau
5144
5145 do 1 i=1,mxcoos
5146 si=-mxcoos*dlcoos/2-dlcoos/2+i*dlcoos
5147 do 2 j=1,mxcoox
5148 xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox
5149 do 3 k=1,mxcooy
5150 yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
5151 if(((x-xj)**2+(y-yk)**2+(z-si)**2).gt.rad**2)goto3
5152 kdensh(itau,i,j,k)=kdensh(itau,i,j,k)+1
5153 ktot(itau)=ktot(itau)+1
5154 3 continue
5155 2 continue
5156 1 continue
5157
5158 elseif(ii.eq.2)then
5159
5160 itaux=nint(xpar1)
5161 if(itaux.gt.mxtau)stop'STOP in xjden2: itaux too large'
5162
5163 is=nint(xpar2/dlcoos)
5164 ix=nint(xpar3/dlcoox)
5165 iy=nint(xpar4/dlcooy)
5166 if(is.gt.mxcoos/2)stop'STOP in xjden2: s-range too large'
5167 if(ix.gt.mxcoox/2)stop'STOP in xjden2: x-range too large'
5168 if(iy.gt.mxcooy/2)stop'STOP in xjden2: y-range too large'
5169
5170 do k=mxcooy/2+1-iy,mxcooy/2+iy
5171 write(cy,'(f3.1)')-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
5172 write(ifhi,'(a)') 'openhisto'
5173 write(ifhi,'(a,2f7.2)')'xrange',-is*dlcoos,is*dlcoos
5174 write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox
5175 write(ifhi,'(a)') 'set ityp2d 5'
5176 write(ifhi,'(a)') 'txt "xaxis [z] "'
5177 write(ifhi,'(a)')
5178 *'txt "yaxis x (fm), y='//cy//' fm"'
5179 write(ifhi,'(a,i4)') 'array2d',2*is
5180 do j=mxcoox/2+1-ix,mxcoox/2+ix
5181 do i=mxcoos/2+1-is,mxcoos/2+is
5182 write(ifhi,'(e11.3)')
5183 *kdensh(itaux,i,j,k)/dlcooy/dlcoos/dlcoox/ktot(itaux)
5184 enddo
5185 enddo
5186 write(ifhi,'(a)') ' endarray'
5187 write(ifhi,'(a)') 'closehisto plot2d'
5188 enddo
5189
5190 else
5191
5192 stop'STOP in xjden2: wrong option'
5193
5194 endif
5195
5196 return
5197 end
5198
5199 cc-----------------------------------------------------------------------
5200 c subroutine postscript(iii,ii,ic)
5201 cc-----------------------------------------------------------------------
5202 c include 'epos.inc'
5203 c character*10 color(10)
5204 c if(iii.eq.0)then
5205 c ifps=21
5206 c open(unit=ifps,file='zzz.ps',status='unknown')
5207 c WRITE(ifps,'(a)') '%!PS-Adobe-2.0'
5208 c WRITE(ifps,'(a)') '%%Title: tt2.fig'
5209 c WRITE(ifps,'(a)') '%%Orientation: Portrait'
5210 c WRITE(ifps,'(a)') '%%BeginSetup'
5211 c WRITE(ifps,'(a)') '%%IncludeFeature: *PageSize A4'
5212 c WRITE(ifps,'(a)') '%%EndSetup'
5213 c WRITE(ifps,'(a)') '%%EndComments'
5214 c WRITE(ifps,*) '/l {lineto} bind def'
5215 c WRITE(ifps,*) '/rl {rlineto} bind def'
5216 c WRITE(ifps,*) '/m {moveto} bind def'
5217 c WRITE(ifps,*) '/rm {rmoveto} bind def'
5218 c WRITE(ifps,*) '/s {stroke} bind def'
5219 c WRITE(ifps,*) '/gr {grestore} bind def'
5220 c WRITE(ifps,*) '/gs {gsave} bind def'
5221 c WRITE(ifps,*) '/cp {closepath} bind def'
5222 c WRITE(ifps,*) '/tr {translate} bind def'
5223 c WRITE(ifps,*) '/sc {scale} bind def'
5224 c WRITE(ifps,*) '/sd {setdash} bind def'
5225 c WRITE(ifps,*) '/sdo {[.01 .05] 0 sd} bind def'
5226 c WRITE(ifps,*) '/sdf {[1 .0] 0 sd} bind def'
5227 c WRITE(ifps,*) '/n {newpath} bind def'
5228 c WRITE(ifps,*) '/slw {setlinewidth } bind def'
5229 c write(ifps,*) '/srgb {setrgbcolor} bind def'
5230 c write(ifps,*) '/lgrey { 0 0.95 0.95 srgb} bind def'
5231 c write(ifps,*) '/black { 0 0 0 srgb} bind def'
5232 c write(ifps,*) '/red { 1 0 0 srgb} bind def '
5233 c write(ifps,*) '/green { 0 1 0 srgb} bind def '
5234 c write(ifps,*) '/blue { 0 0 1 srgb} bind def '
5235 c write(ifps,*) '/yellow { 1 0.5 0 srgb} bind def '
5236 c write(ifps,*) '/turquoise { 0 1 1 srgb} bind def '
5237 c write(ifps,*) '/purple { 1 0 1 srgb} bind def '
5238 cc write(ifps,*) '/ { srgb} bind def '
5239 cc write(ifps,*) '/ { srgb} bind def '
5240 c write(ifps,*) '/ef {eofill} bind def'
5241 c WRITE(ifps,'(a)') '%%EndProlog'
5242 c WRITE(ifps,*) 'gsave'
5243 c WRITE(ifps,*) '/Helvetica findfont 10 scalefont setfont'
5244 c color(9)='lgrey '
5245 c color(1)='black '
5246 c color(2)='red '
5247 c color(3)='green '
5248 c color(4)='blue '
5249 c color(7)='yellow '
5250 c color(5)='turquoise '
5251 c color(6)='purple '
5252 c np=0
5253 c elseif(iii.eq.1)then
5254 c np=np+1
5255 c write(ifps,'(a,i4)') '%%Page: number ',np
5256 c write(ifps,'(a)') 'gsave'
5257 c WRITE(ifps,*) '100 700 tr'
5258 c scale=0.125
5259 c WRITE(ifps,*) 1./scale,1./scale,' sc'
5260 c WRITE(ifps,*) scale/2.,' slw'
5261 c WRITE(ifps,*) '/Helvetica findfont ',15.*scale
5262 c & ,' scalefont setfont'
5263 c write(ifps,*) color(1),' n ',smin,xmin,' m ( tau:',tau,') show '
5264 c
5265 c WRITE(ifps,*) '/Helvetica findfont ',5.*scale
5266 c & ,' scalefont setfont'
5267 c
5268 c
5269 c yb=-2.
5270 c dy=4./12.
5271 c yb=yb-dy/2
5272 c do iyb=0,11
5273 c yb=yb+dy
5274 c WRITE(ifps,*) 'gsave'
5275 c WRITE(ifps,*) (xmax-xmin)*1.1*float(int(iyb/4))
5276 c & ,-(xmax-xmin)*1.1*mod(iyb,4),' tr'
5277 c write(ifps,*) ' n ',smin,xmin,' m ',smax,xmin,' l '
5278 c & ,smax,xmax,' l ',smin,xmax,' l cp s '
5279 cc.......particles in layer iyb.............
5280 c do i=1,nptl
5281 c if(ii.gt.0)then
5282 c write(ifps,*) color(mod(i,5)+2)
5283 c & ,' n ',u,x-r,' m ',o,x-r,' l '
5284 c & ,o,x+r,' l ',u,x+r,' l cp s '
5285 c write(ifps,*) ' n ',u,x-r,' m (',i,ior,') show '
5286 c else
5287 c write(ifps,*) ' n ',s,x,r,0,360,' arc ',color(ic),' s '
5288 c write(ifps,*) ' n ',s-r,x,' m (',i,io,') show '
5289 c endif
5290 c 10 enddo
5291 c write(ifps,*) color(1),' n ',smin,xmin,' m (',yb,') show '
5292 c WRITE(ifps,*) 'grestore'
5293 c enddo !yb bin
5294 c write(ifps,'(a)') 'grestore'
5295 c write(ifps,*) 'showpage'
5296 c elseif(iii.eq.2)then
5297 c write(ifps,*) 'gr'
5298 c
5299 c write(ifps,'(a)') '%%Trailer'
5300 c write(ifps,'(a,i4)') '%%Pages: ',np
5301 c write(ifps,'(a)') '%%EOF'
5302 c close(unit=ifps)
5303 c endif
5304 c
5305 c return
5306 c end
5307 c
5308
5309 c------------------------------------------------------------------------------
5310 subroutine xtauev(iii)
5311 c------------------------------------------------------------------------------
5312 jdum=iii
5313 end
5314 c------------------------------------------------------------------------------
5315 subroutine wimi
5316 c------------------------------------------------------------------------------
5317 end
5318 c------------------------------------------------------------------------------
5319 subroutine wimino
5320 c------------------------------------------------------------------------------
5321 end
5322 c------------------------------------------------------------------------------
5323 subroutine xspace(iii)
5324 c------------------------------------------------------------------------------
5325 jdum=iii
5326 end
5327 c------------------------------------------------------------------------------
5328 subroutine wclu
5329 c------------------------------------------------------------------------------
5330 end
5331 c------------------------------------------------------------------------------
5332 subroutine wclufi
5333 c------------------------------------------------------------------------------
5334 end
5335 c------------------------------------------------------------------------------
5336 subroutine wtime(iii)
5337 c------------------------------------------------------------------------------
5338 jdum=iii
5339 end