-
Notifications
You must be signed in to change notification settings - Fork 7
/
linana.f
executable file
·3536 lines (3501 loc) · 128 KB
/
linana.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
c.................... linana.f
c=======================================================================
c 10/26/01 update freezeout positions in case of interactions:
clin-3/2009 Note: freezeout spacetime values cannot be trusted for K0S & K0L
c as K0S/K0L are converted from K+/K- by hand at the end of hadron cascade.
subroutine hbtout(nnew,nt,ntmax)
c
PARAMETER (MAXSTR=150001,MAXR=1)
clin-5/2008 give tolerance to regular particles (perturbative probability 1):
PARAMETER (oneminus=0.99999,oneplus=1.00001)
dimension lastkp(MAXSTR), newkp(MAXSTR),xnew(3)
common /para7/ ioscar,nsmbbbar,nsmmeson
cc SAVE /para7/
COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
cc SAVE /hbt/
common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
cc SAVE /input1/
COMMON /AA/ R(3,MAXSTR)
cc SAVE /AA/
COMMON /BB/ P(3,MAXSTR)
cc SAVE /BB/
COMMON /CC/ E(MAXSTR)
cc SAVE /CC/
COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
cc SAVE /EE/
common /lastt/itimeh,bimp
cc SAVE /lastt/
COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
cc SAVE /tdecay/
COMMON /AREVT/ IAEVT, IARUN, MISS
cc SAVE /AREVT/
common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
cc SAVE /snn/
COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
cc SAVE /HJGLBR/
COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
clin-12/14/03:
COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
EXTERNAL IARFLV, INVFLV
common /para8/ idpert,npertd,idxsec
clin-2/2012:
common /phiHJ/iphirp,phiRP
SAVE
c
do 1001 i=1,max0(nlast,nnew)
lastkp(i)=0
1001 continue
do 1002 i=1,nnew
newkp(i)=0
1002 continue
c for each of the particles, search the freezeout record (common /hbt/)
c to find & keep those which do not have interactions during this timestep:
do 100 ip=1,nnew
do 1004 iplast=1,nlast
if(p(1,ip).eq.plast(1,iplast).and.
1 p(2,ip).eq.plast(2,iplast).and.
2 p(3,ip).eq.plast(3,iplast).and.
3 e(ip).eq.plast(4,iplast).and.
4 lb(ip).eq.lblast(iplast).and.
5 dpertp(ip).eq.dplast(iplast).and.lastkp(iplast).eq.0) then
clin-5/2008 modified below to the above in case we have perturbative particles:
c 5 lastkp(iplast).eq.0) then
deltat=nt*dt-xlast(4,iplast)
ene=sqrt(plast(1,iplast)**2+plast(2,iplast)**2
1 +plast(3,iplast)**2+plast(4,iplast)**2)
c xnew gives the coordinate if a particle free-streams to current time:
do 1003 ii=1,3
xnew(ii)=xlast(ii,iplast)+plast(ii,iplast)/ene*deltat
1003 continue
dr=sqrt((r(1,ip)-xnew(1))**2+(r(2,ip)-xnew(2))**2
1 +(r(3,ip)-xnew(3))**2)
c find particles with dp=0 and dr<0.01, considered to be those
c without any interactions during this timestep,
c thus keep their last positions and time:
if(dr.le.0.01) then
lastkp(iplast)=1
newkp(ip)=1
c if(lb(ip).eq.41) then
c write(95,*) 'nt,ip,px,x=',nt,ip,p(1,ip),r(1,ip),ftsv(ip)
c write(95,*) 'xnew=',xnew(1),xnew(2),xnew(3),xlast(4,ip)
c endif
clin-5/2009 Take care of formation time of particles read in at nt=ntmax-1:
if(nt.eq.ntmax.and.ftsv(ip).gt.((ntmax-1)*dt))
1 xlast(4,iplast)=ftsv(ip)
goto 100
endif
endif
1004 continue
100 continue
c for current particles with interactions, fill their current info in
c the freezeout record (if that record entry needs not to be kept):
do 150 ip=1,nnew
if(newkp(ip).eq.0) then
do 1005 iplast=1,nnew
if(lastkp(iplast).eq.0) then
ctest off: write collision info
c if(lb(ip).eq.41) then
c write(95,*) 'nt,lb(ip)=',nt,lb(ip)
c write(95,*) ' last p=',plast(1,iplast),
c 1 plast(2,iplast),plast(3,iplast),plast(4,iplast)
c write(95,*) ' after p=',p(1,ip),p(2,ip),p(3,ip),e(ip)
c write(95,*) 'after x=',r(1,ip),r(2,ip),r(3,ip),ftsv(ip)
c endif
c
xlast(1,iplast)=r(1,ip)
xlast(2,iplast)=r(2,ip)
xlast(3,iplast)=r(3,ip)
xlast(4,iplast)=nt*dt
c
if(nt.eq.ntmax) then
c freezeout time for decay daughters at the last timestep
c needs to include the decay time of the parent:
if(tfdcy(ip).gt.(ntmax*dt+0.001)) then
xlast(4,iplast)=tfdcy(ip)
c freezeout time for particles unformed at the next-to-last timestep
c needs to be their formation time instead of (ntmax*dt):
elseif(ftsv(ip).gt.((ntmax-1)*dt)) then
xlast(4,iplast)=ftsv(ip)
endif
endif
plast(1,iplast)=p(1,ip)
plast(2,iplast)=p(2,ip)
plast(3,iplast)=p(3,ip)
plast(4,iplast)=e(ip)
lblast(iplast)=lb(ip)
lastkp(iplast)=1
clin-5/2008:
dplast(iplast)=dpertp(ip)
goto 150
endif
1005 continue
endif
150 continue
c if the current particle list is shorter than the freezeout record,
c condense the last-collision record by filling new record from 1 to nnew,
c and label these entries as keep:
if(nnew.lt.nlast) then
do 170 iplast=1,nlast
if(lastkp(iplast).eq.0) then
do 1006 ip2=iplast+1,nlast
if(lastkp(ip2).eq.1) then
xlast(1,iplast)=xlast(1,ip2)
xlast(2,iplast)=xlast(2,ip2)
xlast(3,iplast)=xlast(3,ip2)
xlast(4,iplast)=xlast(4,ip2)
plast(1,iplast)=plast(1,ip2)
plast(2,iplast)=plast(2,ip2)
plast(3,iplast)=plast(3,ip2)
plast(4,iplast)=plast(4,ip2)
lblast(iplast)=lblast(ip2)
lastkp(iplast)=1
clin-5/2008:
dplast(iplast)=dplast(ip2)
goto 170
endif
1006 continue
endif
170 continue
endif
nlast=nnew
ctest off look inside each NT timestep (for debugging purpose):
c do ip=1,nlast
c write(99,*) ' p ',nt,ip,lblast(ip),plast(1,ip),
c 1 plast(2,ip),plast(3,ip),plast(4,ip),dplast(ip)
c write(99,*) ' x ',nt,ip,lblast(ip),xlast(1,ip),
c 1 xlast(2,ip),xlast(3,ip),xlast(4,ip),dplast(ip)
c enddo
c
if(nt.eq.ntmax) then
clin-5/2008 find final number of perturbative particles (deuterons only):
ndpert=0
do ip=1,nlast
if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
else
ndpert=ndpert+1
endif
enddo
c
c write(16,190) IAEVT,IARUN,nlast,bimp,npart1,npart2,
c 1 NELP,NINP,NELT,NINTHJ
clin-2/2012:
c write(16,190) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2,
c 1 NELP,NINP,NELT,NINTHJ
write(16,191) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2,
1 NELP,NINP,NELT,NINTHJ,phiRP
clin-5/2008 write out perturbatively-produced particles (deuterons only):
if(idpert.eq.1.or.idpert.eq.2)
1 write(90,190) IAEVT,IARUN,ndpert,bimp,npart1,npart2,
2 NELP,NINP,NELT,NINTHJ
do 1007 ip=1,nlast
clin-12/14/03 No formation time for spectator projectile or target nucleons,
c see ARINI1 in 'amptsub.f':
clin-3/2009 To be consistent with new particles produced in hadron cascade
c that are limited by the time-resolution (DT) of the hadron cascade,
c freezeout time of spectator projectile or target nucleons is written as
c DT as they are read at the 1st timestep and then propagated to time DT:
c
clin-9/2011 determine spectator nucleons consistently
c if(plast(1,ip).eq.0.and.plast(2,ip).eq.0
c 1 .and.(sqrt(plast(3,ip)**2+plast(4,ip)**2)*2/HINT1(1))
c 2 .gt.0.99.and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then
if(abs(plast(1,ip)).le.epsiPt.and.abs(plast(2,ip)).le.epsiPt
1 .and.(plast(3,ip).gt.amax1(0.,PZPROJ-epsiPz)
2 .or.plast(3,ip).lt.(-PZTARG+epsiPz))
3 .and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then
clin-5/2008 perturbatively-produced particles (currently only deuterons)
c are written to ana/ampt_pert.dat (without the column for the mass);
c ana/ampt.dat has regularly-produced particles (including deuterons);
c these two sets of deuteron data are close to each other(but not the same
c because of the bias from triggering the perturbative production);
c ONLY use one data set for analysis to avoid double-counting:
if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
write(16,200) INVFLV(lblast(ip)), plast(1,ip),
1 plast(2,ip),plast(3,ip),plast(4,ip),
2 xlast(1,ip),xlast(2,ip),xlast(3,ip),
3 xlast(4,ip)
clin-12/14/03-end
else
if(idpert.eq.1.or.idpert.eq.2) then
write(90,250) INVFLV(lblast(ip)), plast(1,ip),
1 plast(2,ip),plast(3,ip),
2 xlast(1,ip),xlast(2,ip),xlast(3,ip),
3 xlast(4,ip)
else
write(99,*) 'Unexpected perturbative particles'
endif
endif
elseif(amax1(abs(xlast(1,ip)),abs(xlast(2,ip)),
1 abs(xlast(3,ip)),abs(xlast(4,ip))).lt.9999) then
if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
write(16,200) INVFLV(lblast(ip)), plast(1,ip),
1 plast(2,ip),plast(3,ip),plast(4,ip),
2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
else
if(idpert.eq.1.or.idpert.eq.2) then
write(90,250) INVFLV(lblast(ip)),plast(1,ip),
1 plast(2,ip),plast(3,ip),
2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
3 dplast(ip)
else
write(99,*) 'Unexpected perturbative particles'
endif
endif
else
c change format for large numbers:
if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
write(16,201) INVFLV(lblast(ip)), plast(1,ip),
1 plast(2,ip),plast(3,ip),plast(4,ip),
2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
else
if(idpert.eq.1.or.idpert.eq.2) then
write(90,251) INVFLV(lblast(ip)), plast(1,ip),
1 plast(2,ip),plast(3,ip),
2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
3 dplast(ip)
else
write(99,*) 'Unexpected perturbative particles'
endif
endif
endif
1007 continue
if(ioscar.eq.1) call hoscar
endif
190 format(3(i7),f10.4,5x,6(i4))
191 format(3(i7),f10.4,5x,6(i4),5x,f7.4)
clin-3/2009 improve the output accuracy of Pz
200 format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,f8.2))
201 format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,e8.2))
250 format(I5,2(1x,f8.3),1x,f10.3,2(1x,f7.1),1x,f8.2,1x,f7.2,1x,e10.4)
251 format(I5,2(1x,f8.3),1x,f10.3,4(1x,e8.2),1x,e10.4)
c
return
end
c=======================================================================
SUBROUTINE decomp(px0,py0,pz0,xm0,i,itq1)
c
IMPLICIT DOUBLE PRECISION(D)
DOUBLE PRECISION enenew, pxnew, pynew, pznew
clin-8/2015 changed ptwo(2,5) and related variables to double precision
c to avoid IEEE_DIVIDE_BY_ZERO or IEEE_INVALID or IEEE_OVERFLOW_FLAG:
DOUBLE PRECISION de0, beta2, gam, ptwo, px0, py0, pz0, xm0
common /lor/ enenew, pxnew, pynew, pznew
cc SAVE /lor/
COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
cc SAVE /HPARNT/
common /decom/ptwo(2,5)
cc SAVE /decom/
COMMON/RNDF77/NSEED
cc SAVE /RNDF77/
COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd,
1 psembd,tmaxembd,phidecomp
SAVE
c
dcth=dble(RANART(NSEED))*2.d0-1.d0
dPHI=dble(RANART(NSEED)*HIPR1(40))*2.d0
clin-6/2009 Added if embedding a high-Pt quark pair after string melting:
if(iembed.ge.1.and.iembed.le.4) then
c Decompose the parent high-Pt pion to q and qbar with an internal momentum
c parallel to the pion direction so that one parton has ~the same hight Pt
c and the other parton has a very soft Pt:
c Note: htop() decomposes a meson to q as it(1) followed by qbar as it(2):
if(i.eq.(natt-2*nsembd).or.i.eq.(natt-2*nsembd-1)) then
dcth=0.d0
dphi=dble(phidecomp)
endif
endif
c
ds=xm0**2
dpcm=dsqrt((ds-(ptwo(1,5)+ptwo(2,5))**2)
1 *(ds-(ptwo(1,5)-ptwo(2,5))**2)/ds/4d0)
dpz=dpcm*dcth
dpx=dpcm*dsqrt(1.d0-dcth**2)*dcos(dphi)
dpy=dpcm*dsqrt(1.d0-dcth**2)*dsin(dphi)
de1=dsqrt(ptwo(1,5)**2+dpcm**2)
de2=dsqrt(ptwo(2,5)**2+dpcm**2)
c
de0=dsqrt(px0**2+py0**2+pz0**2+xm0**2)
dbex=px0/de0
dbey=py0/de0
dbez=pz0/de0
c boost the reference frame up by beta (pznew=gam(pz+beta e)):
beta2 = dbex ** 2 + dbey ** 2 + dbez ** 2
gam = 1.d0 / dsqrt(1.d0 - beta2)
if(beta2.ge.0.9999999999999d0) then
write(6,*) '1',dbex,dbey,dbez,beta2,gam
endif
c
call lorenz(de1,dpx,dpy,dpz,-dbex,-dbey,-dbez)
ptwo(1,1)=pxnew
ptwo(1,2)=pynew
ptwo(1,3)=pznew
ptwo(1,4)=enenew
call lorenz(de2,-dpx,-dpy,-dpz,-dbex,-dbey,-dbez)
ptwo(2,1)=pxnew
ptwo(2,2)=pynew
ptwo(2,3)=pznew
ptwo(2,4)=enenew
c
RETURN
END
c=======================================================================
SUBROUTINE HTOP
c
PARAMETER (MAXSTR=150001)
PARAMETER (MAXPTN=400001)
PARAMETER (MAXIDL=4001)
DOUBLE PRECISION GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
1 GXSGS,GYSGS,GZSGS,FTSGS, ptwo, xmdq, ptwox, ptwoy, ptwoz
dimension it(4)
COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
cc SAVE /HMAIN2/
COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
cc SAVE /HMAIN1/
COMMON /PARA1/ MUL
cc SAVE /PARA1/
COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
& PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
& XMASS0(MAXPTN), ITYP0(MAXPTN)
cc SAVE /prec1/
COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
cc SAVE /ilist7/
COMMON /ARPRC/ ITYPAR(MAXSTR),
& GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
& PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
& XMAR(MAXSTR)
cc SAVE /ARPRC/
common /decom/ptwo(2,5)
cc SAVE /decom/
COMMON/RNDF77/NSEED
cc SAVE /RNDF77/
COMMON /NOPREC/ NNOZPC, ITYPN(MAXIDL),
& GXN(MAXIDL), GYN(MAXIDL), GZN(MAXIDL), FTN(MAXIDL),
& PXN(MAXIDL), PYN(MAXIDL), PZN(MAXIDL), EEN(MAXIDL),
& XMN(MAXIDL)
cc SAVE /NOPREC/
COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
cc SAVE /HPARNT/
c 7/20/01: use double precision
c otherwise sometimes beta>1 and gamma diverge in lorenz():
COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
& PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
& GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
& K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
cc SAVE /SOFT/
common/anim/nevent,isoft,isflag,izpc
cc SAVE /anim/
clin-8/2015:
DOUBLE PRECISION vxp0,vyp0,vzp0,xstrg0,ystrg0,xstrg,ystrg
common /precpa/vxp0(MAXPTN),vyp0(MAXPTN),vzp0(MAXPTN),
1 xstrg0(MAXPTN),ystrg0(MAXPTN),
2 xstrg(MAXPTN),ystrg(MAXPTN),istrg0(MAXPTN),istrg(MAXPTN)
c DOUBLE PRECISION vxp0,vyp0,vzp0
c common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN)
cc SAVE /precpa/
common /para7/ ioscar,nsmbbbar,nsmmeson
COMMON /AREVT/ IAEVT, IARUN, MISS
common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
SAVE
c
npar=0
nnozpc=0
clin-5b/2008 calculate the number of hadrons to be converted to q/qbar:
if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
1 then
nsmbbbar=0
nsmmeson=0
do i=1,natt
id=ITYPAR(i)
idabs=iabs(id)
i2=MOD(idabs/10,10)
clin-9/2011 determine spectator nucleons consistently
c if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
c 1 .ge.(HINT1(1)/2*0.99).and.
c 2 .and.(id.eq.2112.or.id.eq.2212)) then
if(abs(PXAR(i)).le.epsiPt.and.abs(PYAR(i)).le.epsiPt
1 .and.(PZAR(i).gt.amax1(0.,PZPROJ-epsiPz)
2 .or.PZAR(i).lt.(-PZTARG+epsiPz))
3 .and.(id.eq.2112.or.id.eq.2212)) then
c spectator proj or targ nucleons without interactions, do not enter ZPC:
elseif(idabs.gt.1000.and.i2.ne.0) then
c baryons to be converted to q/qbar:
nsmbbbar=nsmbbbar+1
elseif((idabs.gt.100.and.idabs.lt.1000)
1 .or.idabs.gt.10000) then
c mesons to be converted to q/qbar:
nsmmeson=nsmmeson+1
endif
enddo
clin-6/2009:
if(ioscar.eq.2.or.ioscar.eq.3) then
write(92,*) iaevt,miss,3*nsmbbbar+2*nsmmeson,
1 nsmbbbar,nsmmeson,natt,natt-nsmbbbar-nsmmeson
endif
c write(92,*) iaevt, 3*nsmbbbar+2*nsmmeson
c write(92,*) ' event#, total # of initial partons after string
c 1 melting'
c write(92,*) 'String melting converts ',nsmbbbar, ' baryons &'
c 1, nsmmeson, ' mesons'
c write(92,*) 'Total # of initial particles= ',natt
c write(92,*) 'Total # of initial particles (gamma,e,muon,...)
c 1 not entering ZPC= ',natt-nsmbbbar-nsmmeson
endif
clin-5b/2008-over
do 100 i=1,natt
id=ITYPAR(i)
idabs=iabs(id)
i4=MOD(idabs/1000,10)
i3=MOD(idabs/100,10)
i2=MOD(idabs/10,10)
i1=MOD(idabs,10)
rnum=RANART(NSEED)
ftime=0.197*PEAR(i)/(PXAR(i)**2+PYAR(i)**2+XMAR(i)**2)
inozpc=0
it(1)=0
it(2)=0
it(3)=0
it(4)=0
c
clin-9/2011 determine spectator nucleons consistently
c if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
c 1 .ge.(HINT1(1)/2*0.99).and.((id.eq.2112).or.(id.eq.2212))) then
if(abs(PXAR(i)).le.epsiPt.and.abs(PYAR(i)).le.epsiPt
1 .and.(PZAR(i).gt.amax1(0.,PZPROJ-epsiPz)
2 .or.PZAR(i).lt.(-PZTARG+epsiPz))
3 .and.(id.eq.2112.or.id.eq.2212)) then
c spectator proj or targ nucleons without interactions, do not enter ZPC:
inozpc=1
elseif(idabs.gt.1000.and.i2.ne.0) then
c baryons:
if(((i4.eq.1.or.i4.eq.2).and.i4.eq.i3)
1 .or.(i4.eq.3.and.i3.eq.3)) then
if(i1.eq.2) then
if(rnum.le.(1./2.)) then
it(1)=i4
it(2)=i3*1000+i2*100+1
elseif(rnum.le.(2./3.)) then
it(1)=i4
it(2)=i3*1000+i2*100+3
else
it(1)=i2
it(2)=i4*1000+i3*100+3
endif
elseif(i1.eq.4) then
if(rnum.le.(2./3.)) then
it(1)=i4
it(2)=i3*1000+i2*100+3
else
it(1)=i2
it(2)=i4*1000+i3*100+3
endif
endif
elseif(i4.eq.1.or.i4.eq.2) then
if(i1.eq.2) then
if(rnum.le.(1./2.)) then
it(1)=i2
it(2)=i4*1000+i3*100+1
elseif(rnum.le.(2./3.)) then
it(1)=i2
it(2)=i4*1000+i3*100+3
else
it(1)=i4
it(2)=i3*1000+i2*100+3
endif
elseif(i1.eq.4) then
if(rnum.le.(2./3.)) then
it(1)=i2
it(2)=i4*1000+i3*100+3
else
it(1)=i4
it(2)=i3*1000+i2*100+3
endif
endif
elseif(i4.ge.3) then
it(1)=i4
if(i3.lt.i2) then
it(2)=i2*1000+i3*100+1
else
it(2)=i3*1000+i2*100+3
endif
endif
c antibaryons:
if(id.lt.0) then
it(1)=-it(1)
it(2)=-it(2)
endif
c isoft=4or5 decompose diquark flavor it(2) to two quarks it(3)&(4):
if(isoft.eq.4.or.isoft.eq.5) then
it(3)=MOD(it(2)/1000,10)
it(4)=MOD(it(2)/100,10)
endif
elseif((idabs.gt.100.and.idabs.lt.1000)
1 .or.idabs.gt.10000) then
c mesons:
if(i3.eq.i2) then
if(i3.eq.1.or.i3.eq.2) then
if(rnum.le.0.5) then
it(1)=1
it(2)=-1
else
it(1)=2
it(2)=-2
endif
else
it(1)=i3
it(2)=-i3
endif
else
if((isign(1,id)*(-1)**i3).eq.1) then
it(1)=i3
it(2)=-i2
else
it(1)=i2
it(2)=-i3
endif
endif
else
c save other particles (leptons and photons) outside of ZPC:
inozpc=1
endif
c
if(inozpc.eq.1) then
NJSGS(i)=0
nnozpc=nnozpc+1
itypn(nnozpc)=ITYPAR(i)
pxn(nnozpc)=PXAR(i)
pyn(nnozpc)=PYAR(i)
pzn(nnozpc)=PZAR(i)
een(nnozpc)=PEAR(i)
xmn(nnozpc)=XMAR(i)
gxn(nnozpc)=GXAR(i)
gyn(nnozpc)=GYAR(i)
gzn(nnozpc)=GZAR(i)
ftn(nnozpc)=FTAR(i)
else
NJSGS(i)=2
ptwo(1,5)=dble(ulmass(it(1)))
ptwo(2,5)=dble(ulmass(it(2)))
call decomp(dble(patt(i,1)),dble(patt(i,2)),
1 dble(patt(i,3)),dble(XMAR(i)),i,it(1))
ipamax=2
if((isoft.eq.4.or.isoft.eq.5)
1 .and.iabs(it(2)).gt.1000) ipamax=1
do 1001 ipar=1,ipamax
npar=npar+1
ityp0(npar)=it(ipar)
px0(npar)=ptwo(ipar,1)
py0(npar)=ptwo(ipar,2)
pz0(npar)=ptwo(ipar,3)
e0(npar)=ptwo(ipar,4)
xmass0(npar)=ptwo(ipar,5)
gx0(npar)=dble(GXAR(i))
gy0(npar)=dble(GYAR(i))
gz0(npar)=dble(GZAR(i))
ft0(npar)=dble(ftime)
lstrg0(npar)=i
lpart0(npar)=ipar
vxp0(npar)=dble(patt(i,1)/patt(i,4))
vyp0(npar)=dble(patt(i,2)/patt(i,4))
vzp0(npar)=dble(patt(i,3)/patt(i,4))
clin-8/2015: set parent string information for this parton:
xstrg(npar)=xstrg0(i)
ystrg(npar)=ystrg0(i)
istrg(npar)=istrg0(i)
1001 continue
200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
c
if((isoft.eq.4.or.isoft.eq.5)
1 .and.iabs(it(2)).gt.1000) then
NJSGS(i)=3
xmdq=ptwo(2,5)
ptwo(1,5)=dble(ulmass(it(3)))
ptwo(2,5)=dble(ulmass(it(4)))
c 8/19/02 avoid actual argument in common blocks of DECOMP:
c call decomp(ptwo(2,1),ptwo(2,2),ptwo(2,3),xmdq)
ptwox=ptwo(2,1)
ptwoy=ptwo(2,2)
ptwoz=ptwo(2,3)
call decomp(ptwox,ptwoy,ptwoz,xmdq,i,it(1))
c
do 1002 ipar=1,2
npar=npar+1
ityp0(npar)=it(ipar+2)
px0(npar)=ptwo(ipar,1)
py0(npar)=ptwo(ipar,2)
pz0(npar)=ptwo(ipar,3)
e0(npar)=ptwo(ipar,4)
xmass0(npar)=ptwo(ipar,5)
gx0(npar)=dble(GXAR(i))
gy0(npar)=dble(GYAR(i))
gz0(npar)=dble(GZAR(i))
ft0(npar)=dble(ftime)
lstrg0(npar)=i
lpart0(npar)=ipar+1
vxp0(npar)=dble(patt(i,1)/patt(i,4))
vyp0(npar)=dble(patt(i,2)/patt(i,4))
vzp0(npar)=dble(patt(i,3)/patt(i,4))
clin-8/2015: set parent string information for this parton:
xstrg(npar)=xstrg0(i)
ystrg(npar)=ystrg0(i)
istrg(npar)=istrg0(i)
1002 continue
endif
c
endif
100 continue
MUL=NPAR
c
clin-5b/2008:
if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
1 then
if((natt-nsmbbbar-nsmmeson).ne.nnozpc)
1 write(92,*) 'Problem with the total # of initial particles
2 (gamma,e,muon,...) not entering ZPC'
if((3*nsmbbbar+2*nsmmeson).ne.npar)
1 write(92,*) 'Problem with the total # of initial partons
2 after string melting'
endif
c
RETURN
END
c=======================================================================
SUBROUTINE PTOH
c
PARAMETER (MAXSTR=150001)
DOUBLE PRECISION gxp,gyp,gzp,ftp,pxp,pyp,pzp,pep,pmp
DOUBLE PRECISION gxp0,gyp0,gzp0,ft0fom,drlocl
DOUBLE PRECISION enenew, pxnew, pynew, pznew, beta2, gam
DOUBLE PRECISION ftavg0,gxavg0,gyavg0,gzavg0,bex,bey,bez
DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
1 GXSGS,GYSGS,GZSGS,FTSGS
DOUBLE PRECISION xmdiag,px1,py1,pz1,e1,px2,py2,pz2,e2,
1 px3,py3,pz3,e3,xmpair,etot
clin-9/2012: improve precision for argument in sqrt():
DOUBLE PRECISION p1,p2,p3
common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
1 pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
cc SAVE /loclco/
COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
cc SAVE /HMAIN1/
COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
cc SAVE /HMAIN2/
COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
& K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
& PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
cc SAVE /HJJET2/
COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
cc SAVE /ARPRNT/
COMMON /ARPRC/ ITYPAR(MAXSTR),
& GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
& PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
& XMAR(MAXSTR)
cc SAVE /ARPRC/
COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
& PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
& GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
& K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
cc SAVE /SOFT/
COMMON/RNDF77/NSEED
cc SAVE /RNDF77/
common/anim/nevent,isoft,isflag,izpc
cc SAVE /anim/
common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom
cc SAVE /prtn23/
common /nzpc/nattzp
cc SAVE /nzpc/
common /lor/ enenew, pxnew, pynew, pznew
cc SAVE /lor/
COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
cc SAVE /LUDAT1/
clin 4/19/2006
common /lastt/itimeh,bimp
COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
COMMON /AREVT/ IAEVT, IARUN, MISS
common /para7/ ioscar,nsmbbbar,nsmmeson
clin-5/2011
common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
c
dimension xmdiag(MAXSTR),indx(MAXSTR),ndiag(MAXSTR)
SAVE
c
call coales
c obtain particle mass here without broadening by Breit-Wigner width:
mstj24=MSTJ(24)
MSTJ(24)=0
nuudd=0
npich=0
nrhoch=0
ppi0=1.
prho0=0.
c determine hadron flavor except (pi0,rho0,eta,omega):
DO 1001 ISG = 1, NSG
if(NJSGS(ISG).ne.0) then
NATT=NATT+1
K1=K2SGS(ISG,1)
k1abs=iabs(k1)
PX1=PXSGS(ISG,1)
PY1=PYSGS(ISG,1)
PZ1=PZSGS(ISG,1)
K2=K2SGS(ISG,2)
k2abs=iabs(k2)
PX2=PXSGS(ISG,2)
PY2=PYSGS(ISG,2)
PZ2=PZSGS(ISG,2)
c 5/02/01 try lowest spin states as first choices,
c i.e. octet baryons and pseudoscalar mesons (ibs=2*baryonspin+1):
e1=PESGS(ISG,1)
e2=PESGS(ISG,2)
xmpair=dsqrt((e1+e2)**2-(px1+px2)**2-(py1+py2)**2
1 -(pz1+pz2)**2)
ibs=2
imspin=0
if(k1.eq.-k2.and.iabs(k1).le.2.
1 and.NJSGS(ISG).eq.2) then
nuudd=nuudd+1
xmdiag(nuudd)=xmpair
ndiag(nuudd)=natt
endif
K3=0
if((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3) then
K3=K2SGS(ISG,3)
k3abs=iabs(k3)
PX3=PXSGS(ISG,3)
PY3=PYSGS(ISG,3)
PZ3=PZSGS(ISG,3)
e3=PESGS(ISG,3)
xmpair=dsqrt((e1+e2+e3)**2-(px1+px2+px3)**2
1 -(py1+py2+py3)**2-(pz1+pz2+pz3)**2)
endif
c***** isoft=3 baryon decomposition is different:
if(isoft.eq.3.and.
1 (k1abs.gt.1000.or.k2abs.gt.1000)) then
if(k1abs.gt.1000) then
kdq=k1abs
kk=k2abs
else
kdq=k2abs
kk=k1abs
endif
ki=MOD(kdq/1000,10)
kj=MOD(kdq/100,10)
if(MOD(kdq,10).eq.1) then
idqspn=0
else
idqspn=1
endif
c
if(kk.gt.ki) then
ktemp=kk
kk=kj
kj=ki
ki=ktemp
elseif(kk.gt.kj) then
ktemp=kk
kk=kj
kj=ktemp
endif
c
if(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
if(idqspn.eq.0) then
kf=1000*ki+100*kk+10*kj+ibs
else
kf=1000*ki+100*kj+10*kk+ibs
endif
elseif(ki.eq.kj.and.ki.eq.kk) then
c can only be decuplet baryons:
kf=1000*ki+100*kj+10*kk+4
else
kf=1000*ki+100*kj+10*kk+ibs
endif
c form a decuplet baryon if the q+diquark mass is closer to its mass
c (and if the diquark has spin 1):
cc for now only include Delta, which is present in ART:
cc if(idqspn.eq.1.and.MOD(kf,10).eq.2) then
if(kf.eq.2112.or.kf.eq.2212) then
if(abs(sngl(xmpair)-ULMASS(kf)).gt.
1 abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
endif
if(k1.lt.0) kf=-kf
clin-6/22/01 isoft=4or5 baryons:
elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
1 then
if(k1abs.gt.k2abs) then
ki=k1abs
kk=k2abs
else
ki=k2abs
kk=k1abs
endif
if(k3abs.gt.ki) then
kj=ki
ki=k3abs
elseif(k3abs.lt.kk) then
kj=kk
kk=k3abs
else
kj=k3abs
endif
c
if(ki.eq.kj.and.ki.eq.kk) then
c can only be decuplet baryons (Delta-,++, Omega):
ibs=4
kf=1000*ki+100*kj+10*kk+ibs
elseif(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
c form Lambda or Sigma according to 3-quark mass,
c for now neglect decuplet (Sigma*0 etc) which is absent in ART:
ibs=2
kf1=1000*ki+100*kj+10*kk+ibs
kf2=1000*ki+100*kk+10*kj+ibs
kf=kf1
if(abs(sngl(xmpair)-ULMASS(kf1)).gt.
1 abs(sngl(xmpair)-ULMASS(kf2))) kf=kf2
else
ibs=2
kf=1000*ki+100*kj+10*kk+ibs
cc for now only include Delta0,+ as decuplets, which are present in ART:
if(kf.eq.2112.or.kf.eq.2212) then
if(abs(sngl(xmpair)-ULMASS(kf)).gt.
1 abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
endif
endif
if(k1.lt.0) kf=-kf
c***** mesons:
else
if(k1abs.eq.k2abs) then
if(k1abs.le.2) then
c treat diagonal mesons later in the subroutine:
kf=0
elseif(k1abs.le.3) then
c do not form eta', only form phi from s-sbar, since no eta' in ART:
kf=333
else
kf=100*k1abs+10*k1abs+2*imspin+1
endif
else
if(k1abs.gt.k2abs) then
kmax=k1abs
kmin=k2abs
elseif(k1abs.lt.k2abs) then
kmax=k2abs
kmin=k1abs
endif
kf=(100*kmax+10*kmin+2*imspin+1)
1 *isign(1,k1+k2)*(-1)**kmax
c form a vector meson if the q+qbar mass is closer to its mass:
if(MOD(iabs(kf),10).eq.1) then
if(abs(sngl(xmpair)-ULMASS(iabs(kf))).gt.
1 abs(sngl(xmpair)-ULMASS(iabs(kf)+2)))
2 kf=(iabs(kf)+2)*isign(1,kf)
endif
endif
endif
ITYPAR(NATT)=kf
KATT(NATT,1)=kf
if(iabs(kf).eq.211) then
npich=npich+1
elseif(iabs(kf).eq.213) then
nrhoch=nrhoch+1
endif
endif
clin-7/2011-check charm hadron flavors:
c if(k1abs.eq.4.or.k2abs.eq.4) then
c if(k3.eq.0) then
c write(99,*) iaevt,k1,k2,kf,xmpair,
c 1 ULMASS(iabs(kf)),ULMASS(iabs(kf)+2),isg
c else
c write(99,*) iaevt,k1,k2,k3,kf,xmpair,
c 1 ULMASS(iabs(kf)),ULMASS(iabs(kf)+2),isg
c endif
c endif
clin-7/2011-end
1001 CONTINUE
c assume Npi0=(Npi+ + Npi-)/2, Nrho0=(Nrho+ + Nrho-)/2 on the average:
if(nuudd.ne.0) then
ppi0=float(npich/2)/float(nuudd)
prho0=float(nrhoch/2)/float(nuudd)
endif
c determine diagonal mesons (pi0,rho0,eta and omega) from uubar/ddbar:
npi0=0
DO 1002 ISG = 1, NSG
if(K2SGS(ISG,1).eq.-K2SGS(ISG,2)
1 .and.iabs(K2SGS(ISG,1)).le.2.and.NJSGS(ISG).eq.2) then
if(RANART(NSEED).le.ppi0) npi0=npi0+1
endif
1002 CONTINUE
c
if(nuudd.gt.1) then
call index1(MAXSTR,nuudd,xmdiag,indx)
else
indx(1)=1
end if
c
DO 1003 ix=1,nuudd
iuudd=indx(ix)
inatt=ndiag(iuudd)
if(ix.le.npi0) then
kf=111
elseif(RANART(NSEED).le.(prho0/(1-ppi0+0.00001))) then
kf=113
else
c at T=150MeV, thermal weights for eta and omega(spin1) are about the same:
if(RANART(NSEED).le.0.5) then
kf=221
else
kf=223
endif
endif
ITYPAR(inatt)=kf
KATT(inatt,1)=kf
1003 CONTINUE
c determine hadron formation time, position and momentum:
inatt=0
clin-6/2009 write out parton info after coalescence:
if(ioscar.eq.3) then
WRITE (85, 395) IAEVT, 3*nsmbbbar+2*nsmmeson,nsmbbbar,nsmmeson,
1 bimp, NELP,NINP,NELT,NINTHJ,MISS
endif
395 format(4I8,f10.4,5I5)
c
DO 1006 ISG = 1, NSG
if(NJSGS(ISG).ne.0) then
inatt=inatt+1
K1=K2SGS(ISG,1)
k1abs=iabs(k1)
PX1=PXSGS(ISG,1)
PY1=PYSGS(ISG,1)
PZ1=PZSGS(ISG,1)
K2=K2SGS(ISG,2)
k2abs=iabs(k2)
PX2=PXSGS(ISG,2)
PY2=PYSGS(ISG,2)
PZ2=PZSGS(ISG,2)
e1=PESGS(ISG,1)
e2=PESGS(ISG,2)
c
if(NJSGS(ISG).eq.2) then
PXAR(inatt)=sngl(px1+px2)
PYAR(inatt)=sngl(py1+py2)
PZAR(inatt)=sngl(pz1+pz2)
PATT(inatt,1)=PXAR(inatt)
PATT(inatt,2)=PYAR(inatt)
PATT(inatt,3)=PZAR(inatt)
etot=e1+e2
clin-9/2012: improve precision for argument in sqrt():
p1=px1+px2
p2=py1+py2
p3=pz1+pz2
c
elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
1 then