-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathparini_30.f
8039 lines (7440 loc) · 268 KB
/
parini_30.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
subroutine parini(time_ini,parp2,win,psno,ijk,iMode,
c decpro,i_tune,time_par,time_had) ! 280524
! 260223 300623 Lei Added i_tune 220823 Lei Removed parp22
c210921 generate partonic initial state in relativistic
c pA,Ap,AA,lp, & lA collision based on 'PYTHIA' for C-framework
csa011223
c (iMode=3). generate intermediate final hadronic state before hadronic
c rescattering for A- and B-frameworks (iMode=1 and 2, respectively)
c it was composed by Ben-Hao Sa on 04/12/2003
c260223 input message is in 'pyjets'
c intermediate working arraies are in 'sa2'
c110123 'saf' is consistent with 'sa2'
c 'saf' to 'pyjets' after call 'scat' ! 220110
c110123 output message is in 'pyjets' (partons) and 'sbh' (hadrons) for case
c260223 of mstptj=0 (C-framework), but is in 'sbh' for mstptj=1 (A- and
csa011223 B-frameworks)
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
IMPLICIT INTEGER(I-N)
INTEGER PYK,PYCHGE,PYCOMP
PARAMETER (KSZJ=80000)
PARAMETER (NSIZE=280000)
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
c those variables in above common blocks are defined in 'jetset'
COMMON/PYSUBS/MSEL,MSUB(500),KFIN(2,-40:40),NON,CKIN(200)
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) ! 221203
c those variables in above common block are defined in 'PYTHIA'
COMMON/PYCIDAT2/KFMAXT,NONCI2,PARAM(20),WEIGH(600)
common/sa1/kjp21,non1,bp,iii,neve,nout,nosc
common/sa2/nsa,non2,ksa(kszj,5),psa(kszj,5),vsa(kszj,5)
common/sa4/tau(kszj),tlco(kszj,4)
common/sa5/kfmax,kfaco(100),numb(100),numbs(100),non5,
c disbe(100,100)
common/sa6/kfmaxi,nwhole
common/sa10/csnn,cspin,cskn,cspipi,cspsn,cspsm,rcsit,ifram,
& iabsb,iabsm,non10,ajpsi,csspn,csspm,csen ! 060813
c080104
common/sa14/ipyth(2000),idec(2000),iwide(2000)
common/sa21/pincl(5),pscal(5),pinch(5),vnu,fq2,w2l,yyl,zl,xb,pph
c ,vnlep ! 260314
common/sa24/adj1(40),nnstop,non24,zstop ! 140414
common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio ! 220110
common/sa27/itime,kjp22,gtime,astr,akapa(6),parj1,parj2,parj3,
c parj21,parj4,adiv,gpmax,nnc ! 070417 010518
common/sa30/vneump,vneumt,mstptj ! 241110 100821 230722
common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
common/saf/naf,nonaf,kaf(kszj,5),paf(kszj,5),vaf(kszj,5)
c080104
common/sbh/nbh,nonbh,kbh(kszj,5),pbh(kszj,5),vbh(kszj,5)
common/wz/c17(500,3),ishp(kszj),tp(500),coor(3),p17(500,4)
common/count/isinel(600)
common/papr/t0,sig,dep,ddt,edipi,epin,ecsnn,ekn,ecspsn,ecspsm
c ,rnt,rnp,rao,rou0,vneu,vneum,ecsspn,ecsspm,ecsen ! 060813
common/syspar/ipden,itden,suppm,suptm,suppc,suptc,r0p,r0t,
c nap,nat,nzp,nzt,pio
common/ctllist/nctl,noinel(600),nctl0,nctlm ! 180121 230121
common/sa12/ppsa(5),nchan,nsjp,sjp,taup,taujp
common/sa15/nps,npsi,pps(5000,5),ppsi(5000,5)
common/sa23/kpar,knn,kpp,knp,kep ! 200601 060813
common/sa33/smadel,ecce,secce,parecc,iparres ! 270312 240412 131212
common/schuds/schun,schudn,schudsn,sfra ! She and Lei
c iii : number of current event
c csen: e+p total x section in fm^2
c neve : total number of events
c bp : impact parameter
c 'sbe': store initial parton confiquration (with diquark)
c 'saf': store parton configuration after parton rescattering
c (w/o diquark)
c c17(i,1-3) : three position of i-th nucleon (origin is set at
c the center of overlap region) ! sa011223
c tp(i) : time of i-th nucleon counted since collision of two nuclei
c p17(i,1-4) : four momentum of i-th nucleon
c ishp(i)=1 if i-th particle inside the simulated volume
c =0 if i-th particle outside the simulated volume
c ecsen: largest collision distance between lepton and p ! 060813
c120214 note: incident lepton collides with nucleon in nucleus once
c only, because of very low total x-section. that collision is the one with
c lowest minimum approaching distance.
c sig (fm^2): cross section of pion + pion to kaon + kaon
c edipi: largest interaction distance between two pions.
c epin: largest interaction distance between pion and nucleon.
c ekn: largest interaction distance between kaon and nucleon.
c ecsnn: largest interaction distance between two nucleons.
c t0 : average proper formation time at rest.
c ddt : time accuracy used in parton initiation
c time accuracy used in parton cascade is dddt
c rou0 : normal nucleon density.
c rao : enlarged factor for the radius of simulated volume.
c nap and nzp (nat and nzt) are the mass and charge numbers of
c projectile (target) nucleus
c r0p=rnp : the standard radius of projectile
c r0t=rnt : the standard radius of target
c nctl: number of collision pairs in current collision list
c nctl0: number of collision pairs in initial collision list
c180121 nctlm: maxmimum number of collision pairs
c noinel(1): statistics of nn elas. colli.;
c noinel(592): statistics of nn colli. calling PYTHIA
c230121 noinel(593): statistics of nn colli. not calling PYTHIA
double precision bst(4),bzp,bzt,bbb(3),bb(3)
dimension peo(4),pi(4),pj(4),xi(4),xj(4)
dimension inoin(kszj)
dimension lc(nsize,5),tc(nsize),tw(nsize)
c-------------------------------------------------------------------------------
c----------------------- Local Variable Initializing -----------------------
c Counters of sub-collisions pairs (zero-pT, nn, pp, np/pn and lp/ln).
kpar=0
knn=0
kpp=0
knp=0
kep=0 ! 060813
c kep: to statistics of the # of times calling PYTHIA in
c case of lepton is projectile
c270312 initiation of x,y,xy,x^2,y^2 and sump (statistics of the number of
c nucleons in overlap region) ! 131212
sumx=0.
sumy=0.
sumxy=0. ! 131212
sumx2=0.
sumy2=0.
sump=0.
adj130=adj1(30) ! 121222 Lei
c270312
c----------------------- Local Variable Initializing -----------------------
c-------------------------------------------------------------------------------
c initiates pp (pA,Ap,AA,lp & lA) collision system
c241110
c creat the initial particle (nucleon) list
c-------------------------------------------------------------------------------
c-------------------------- Position Initializing --------------------------
c230311 in position phase space
c***************************** A+B Collisions ******************************
c191110
c A+B (nucleus-nucleus) ! 230311
if(ipden.eq.1 .and. itden.eq.1)then !! 230311
c distribute projectile nucleons by Woods-Saxon ! 060921
c------------------------ Nuclear Skin Depth Giving ------------------------
napt=nap
if(napt.lt.27)then
alpt=0.47
elseif(napt.gt.27.and.napt.lt.108)then
alpt=0.488
else
alpt=0.54
endif
if(napt.eq.27)then
alpt=0.478
elseif(napt.eq.28)then
alpt=0.48
elseif(napt.eq.32)then
alpt=0.49
elseif(napt.eq.56)then
alpt=0.49
elseif(napt.eq.64)then
alpt=0.49
elseif(napt.eq.108)then
alpt=0.495
elseif(napt.eq.184)then
alpt=0.53
elseif(napt.eq.197)then
alpt=0.54
elseif(napt.eq.207)then
alpt=0.545
elseif(napt.eq.238)then
alpt=0.55
endif
c------------------------ Nuclear Skin Depth Giving ------------------------
alp=alpt
r0=r0p
am=suppm ! upper bound in sampling the radius of projectile nucleon
ac=suppc ! maximum radius for projectile
ratps=vneump/nap ! ratio of projectile participant nucleons to total
do i1=1,nap
if( INT(adj130).eq.1 )then ! 121222 Lei
rann=pyr(1)
if(rann.lt.ratps)then
c sample position of projectile nucleon in overlap region of colliding
c nuclei
call arrove(i1,1,sumx,sumy,sumxy,sumx2,sumy2,sump,
c alp,r0,am,ac) ! 270312 131212 101014
else
c sample position of projectile nucleon according to Woods-Saxon
c distribution
call woodsax_samp(i1,1,alp,r0,am,ac,1) ! 230311
c230311 last argument here is 'iway', iway=1: particle i1 must be outside the
C230311 overlap region of colliding nuclei, iway=0: no more requirement
endif
elseif( INT(adj130).eq.0 )then ! 121222 Lei
call woodsax_samp(i1,1,alp,r0,am,ac,0) ! 230311 060921
endif
enddo
c230311
c distribute target nucleons by Woods-Saxon ! 060921
c------------------------ Nuclear Skin Depth Giving ------------------------
napt=nat
if(napt.lt.27)then
alpt=0.47
elseif(napt.gt.27.and.napt.lt.108)then
alpt=0.488
else
alpt=0.54
endif
if(napt.eq.27)then
alpt=0.478
elseif(napt.eq.28)then
alpt=0.48
elseif(napt.eq.32)then
alpt=0.49
elseif(napt.eq.56)then
alpt=0.49
elseif(napt.eq.64)then
alpt=0.49
elseif(napt.eq.108)then
alpt=0.495
elseif(napt.eq.184)then
alpt=0.53
elseif(napt.eq.197)then
alpt=0.54
elseif(napt.eq.207)then
alpt=0.545
elseif(napt.eq.238)then
alpt=0.55
endif
c------------------------ Nuclear Skin Depth Giving ------------------------
alp=alpt
r0=r0t
am=suptm ! upper bound in sampling the radius of target
ac=suptc ! maximum radius for target
ratps=vneumt/nat ! ratio of target participant nucleons to total
do i1=1,nat
i2=i1+nap
if( INT(adj130).eq.1 )then ! 121222 Lei
rann=pyr(1)
if(rann.lt.ratps)then
c sample position of target nucleon in overlap region of colliding nuclei
call arrove(i2,0,sumx,sumy,sumxy,sumx2,sumy2,sump,
c alp,r0,am,ac) ! 270312 131212 101014
else
c sample position of target nucleon according to Woods-Saxon
c distribution
call woodsax_samp(i2,0,alp,r0,am,ac,1)
endif
elseif( INT(adj130).eq.0 )then ! 121222 Lei
call woodsax_samp(i2,0,alp,r0,am,ac,0) ! 060921
endif
enddo
c------------- Impact-Parameter & Initial Geometry Calculating -------------
c191110
do i=1,nap
c050322 c17(i,1)=c17(i,1)+bp
c17(i,1)=c17(i,1)+0.5*bp ! 050322 move x-component of origin to 0.5*bp
c151222
c Calculates eccentricity correctly for both adj1(30)=0 and 1. ! 151222 Lei
x = c17(i,1)
y = c17(i,2)
z = c17(i,3)
c Relative distance between the projectile nucleon i and the target center (-bp/2., 0, 0)
rel_dist = SQRT( (x+bp/2.)**2 + y**2 + z**2 )
c The projectile nucleon i is inside the target, i.e. inside the overlap region.
if( rel_dist .le. r0t )then
sumx = sumx + x
sumy = sumy + y
sumxy = sumxy + x*y
sumx2 = sumx2 + x**2
sumy2 = sumy2 + y**2
sump = sump + 1.
end if
c151222
enddo
c050322
do i=nap+1,nap+nat
c17(i,1)=c17(i,1)-0.5*bp
c151222
c Calculates eccentricity correctly for both adj1(30)=0 and 1. ! 151222 Lei
x = c17(i,1)
y = c17(i,2)
z = c17(i,3)
c Relative distance between the target nucleon i and the projectile center (+bp/2., 0, 0)
rel_dist = SQRT( (x-bp/2.)**2 + y**2 + z**2 )
c The target nucleon i is inside the projectile, i.e. inside the overlap region.
if( rel_dist .le. r0p )then
sumx = sumx + x
sumy = sumy + y
sumxy = sumxy + x*y
sumx2 = sumx2 + x**2
sumy2 = sumy2 + y**2
sump = sump + 1.
end if
c151222
enddo
c050322
c191110
c------------- Impact-Parameter & Initial Geometry Calculating -------------
c***************************** A+B Collisions ******************************
c*************************** NA & lA Collisions ****************************
c p+A or lepton+A ! 060813 120214
elseif((ipden.eq.0.or.ipden.gt.1) .and. itden.eq.1)then !060813 120214
c100821 distribute projectile proton
do i=1,3
c17(1,i)=0.
if(i.eq.1)c17(1,i)=c17(1,i)+0.5*bp ! 050322 bp->0.5*bp
enddo
c distribute target nucleons by Woods-Saxon ! 180921
c distribute nat-vneumt target nucleons by Woods-Saxon ! 100821
c100821 vneumt: # of target participant nucleons
c180921 ineumt=int(vneumt) ! 100821
c------------------------ Nuclear Skin Depth Giving ------------------------
napt=nat ! -ineumt 180921
if(napt.lt.27)then
alpt=0.47
elseif(napt.gt.27.and.napt.lt.108)then
alpt=0.488
else
alpt=0.54
endif
if(napt.eq.27)then
alpt=0.478
elseif(napt.eq.28)then
alpt=0.48
elseif(napt.eq.32)then
alpt=0.49
elseif(napt.eq.56)then
alpt=0.49
elseif(napt.eq.64)then
alpt=0.49
elseif(napt.eq.108)then
alpt=0.495
elseif(napt.eq.184)then
alpt=0.53
elseif(napt.eq.197)then
alpt=0.54
elseif(napt.eq.207)then
alpt=0.545
elseif(napt.eq.238)then
alpt=0.55
endif
c------------------------ Nuclear Skin Depth Giving ------------------------
alp=alpt
r0=r0t
am=suptm ! upper bound in sampling the radius of target
ac=suptc ! maximum radius for target
do i1=1,napt ! 100821 nat->nat-ineumt=napt
i2=i1+nap
call woodsax_samp(i2,0,alp,r0,am,ac,0)
enddo
c240513
c050322
do i=nap+1,nap+nat
c17(i,1)=c17(i,1)-0.5*bp
enddo
c050322
c*************************** NA & lA Collisions ****************************
c****************************** AN Collisions ******************************
c A+p
elseif(ipden.eq.1 .and. itden.eq.0)then !!
c180921 distribute projectile nucleons by Woods-Saxon
c distribute nap-vneump projectile (spectator) nucleons by Woods-Saxon ! 100821
c100821 vneump: # of projectile participant nucleons
c180921 ineump=int(vneump)
c------------------------ Nuclear Skin Depth Giving ------------------------
napt=nap ! 180921 -ineump
if(napt.lt.27)then
alpt=0.47
elseif(napt.gt.27.and.napt.lt.108)then
alpt=0.488
else
alpt=0.54
endif
if(napt.eq.27)then
alpt=0.478
elseif(napt.eq.28)then
alpt=0.48
elseif(napt.eq.32)then
alpt=0.49
elseif(napt.eq.56)then
alpt=0.49
elseif(napt.eq.64)then
alpt=0.49
elseif(napt.eq.108)then
alpt=0.495
elseif(napt.eq.184)then
alpt=0.53
elseif(napt.eq.197)then
alpt=0.54
elseif(napt.eq.207)then
alpt=0.545
elseif(napt.eq.238)then
alpt=0.55
endif
c------------------------ Nuclear Skin Depth Giving ------------------------
alp=alpt
r0=r0p
am=suppm ! upper bound in sampling the radius of projectile nucleon
ac=suppc ! maximum radius for projectile
do i1=1,napt
call woodsax_samp(i1,1,alp,r0,am,ac,0)
enddo
c191110 100821
do i=1,napt
c17(i,1)=c17(i,1)+0.5*bp ! 050322 bp->0.5*bp
enddo
c191110 100821
c100821 move x-component of origin to 0.5*bp
do i=1,3
c17(nap+1,i)=0.
if(i.eq.1)c17(nap+1,i)=-0.5*bp ! 0.->-0.5*bp
enddo
c240513
c****************************** AN Collisions ******************************
c**************************** NN & lN Collisions ****************************
c p+p or lepton+p ! 070417
elseif((ipden.eq.0 .and. itden.eq.0) .or.
c (itden.eq.0 .and. ipden.ge.11))then !! 070417
do i=1,3
c17(1,i)=0.
c17(2,i)=0.
enddo
endif !!
c**************************** NN & lN Collisions ****************************
c230311
r0pt=r0p+r0t ! 191110
c********************** Initial Geometry Calculating ***********************
c270312
if(sump.gt.0.)then !!! 280224 .ne.0 -> .gt.0.
asumx=sumx/sump
sigmx2=sumx2/sump-asumx*asumx
asumy=sumy/sump
sigmy2=sumy2/sump-asumy*asumy
asumxy=sumxy/sump ! 131212
sigmxy=asumxy-asumx*asumy ! 131212
sigmsu=sigmy2+sigmx2 ! change from sigmxy to sigmsu 131212
sigmde=sigmy2-sigmx2 ! 131212
argu=sigmde*sigmde+4*sigmxy*sigmxy ! 131212
c131212
c participant eccentricity of participant nucleons
if(argu.gt.0. .and. sigmsu.gt.0.)
c ecce=sqrt(argu)/sigmsu !131212
c calculate transverse overlap area
argu1=sigmx2*sigmy2-sigmxy*sigmxy
if(argu1.gt.0.)secce=3.1416*sqrt(argu1) ! overlop area 250113
c131212
c--------------------- Transverse-Momentum Deformation ---------------------
c assuming ecce = geometric eccentricity of ellipsoid
c sqrt( 1-b^2/a^2 )
c with half major axis
c b = pt * ( 1 + smadel )
c and half minor axis
c a = pt * ( 1 - smadel ),
c the resulted
c smadel = -ecce*ecce/4,
c neglecting the samll term of
c ecce*ecce*( -2*smadel + smadel*smadel ).
ecc2=ecce*ecce ! 250113
smadel_a=parecc*ecc2/4. ! approximated deformation parameter 250113
c250113
delta1=(2.-ecc2+2.*(1.-ecc2)**0.5)/ecc2
delta2=(2.-ecc2-2.*(1.-ecc2)**0.5)/ecc2
if(delta1.le.1.)then
smadel=parecc*delta1 ! exact deformation parameter
elseif(delta2.le.1.)then
smadel=parecc*delta2 ! exact deformation parameter
endif
c250113
c here a sign change is introduced because of asymmetry of initial
c spatial space is oppsed to the final momentum space
c--------------------- Transverse-Momentum Deformation ---------------------
endif !!!
c270312
c********************** Initial Geometry Calculating ***********************
c021018 note: psno=0 (bp=0) for pp,lp and lA
c191110
c the beam direction is identified as the z axis
c the origin in position space is set on (0, 0, 0)
c and the origin of time is set at the moment of
c first nn colission assumed to be 1.e-5
c-------------------------- Position Initializing --------------------------
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c-------------------------- Momentum Initializing --------------------------
c230311 in momentum phase space
if(ifram.eq.1)then
ep1=0.5*win ! energy of projetile particle (if it is proton)
et1=ep1 ! energy of target particle (if proton)
ep2=0.5*win ! energy of projetile particle (if neutron)
et2=ep2 ! energy of target particle (if neutron)
pm2=pmas(pycomp(2212),1)**2 ! square mass of proton
pp1=dsqrt(ep1*ep1-pm2) ! momentum of projetile particle (if proton)
pt1=-dsqrt(et1*et1-pm2) ! momentum of target particle (if proton)
pm2=pmas(pycomp(2112),1)**2 ! square mass of nucleon
pp2=dsqrt(ep2*ep2-pm2) ! momentum of projetile particle (if neutron)
pt2=-dsqrt(et2*et2-pm2) ! momentum of target particle (if neutron)
c260314 set four momentum and mass for incident lepton
if(ipden.ge.11.and.ipden.le.16)then ! in cms
pincl(1)=0.
pincl(2)=0.
pincl(4)=0.5d0*win
pincl(5)=pmas(pycomp(ipden),1)
pincl3=pincl(4)*pincl(4)-pincl(5)*pincl(5)
pincl3=dmax1(pincl3,1.d-20)
pincl(3)=dsqrt(pincl3)
pinch(1)=0.
pinch(2)=0.
pinch(4)=0.5d0*win
pinch(5)=pmas(pycomp(2212),1)
pinch3=pinch(4)*pinch(4)-pinch(5)*pinch(5)
pinch3=dmax1(pinch3,1.d-20)
pinch(3)=dsqrt(pinch3)
endif
c260314
endif
if(ifram.eq.0)then
pp1=win ! momentum of projetile particle (if proton)
pt1=1.e-20 ! momentum of target particle (if proton)
pp2=win ! momentum of projetile particle (if neutron)
pt2=1.e-20 ! momentum of target particle (if neutron)
pm2=pmas(pycomp(2212),1)**2 ! square mass of proton
ep1=dsqrt(pp1*pp1+pm2) ! energy of projetile particle (if proton)
et1=dsqrt(pt1*pt1+pm2) ! energy of target particle (if proton)
pm2=pmas(pycomp(2112),1)**2 ! square mass of neutron
ep2=dsqrt(pp2*pp2+pm2) ! energy of projetile particle (if neutron)
et2=dsqrt(pt2*pt2+pm2) ! energy of target particle (if neutron)
c260314 set four momentum and mass for incident lepton
if(ipden.ge.11.and.ipden.le.16)then ! in lab
pincl(1)=0.
pincl(2)=0.
pincl(3)=win
pincl(5)=pmas(pycomp(ipden),1)
pincl4=pincl(3)*pincl(3)+pincl(5)*pincl(5)
pincl4=dmax1(pincl4,1.d-20)
pincl(4)=dsqrt(pincl4)
pinch(1)=0.
pinch(2)=0.
pinch(3)=0.
pinch(5)=pmas(pycomp(2212),1)
pinch(4)=pinch(5)
endif
c260314
endif
c260314
100 inzp=iabs(nzp)
inzt=iabs(nzt)
do i=1,nap
p17(i,1)=0. ! four momenta of projectile particle i
p17(i,2)=0.
if(i.le.inzp)then
p17(i,3)=pp1 ! projectile particle is proton
p17(i,4)=ep1
else
p17(i,3)=pp2 ! projectile particle is neutron
p17(i,4)=ep2
endif
enddo
napt=nap+nat
do i=nap+1,napt
p17(i,1)=0. ! four momenta of target particle i
p17(i,2)=0.
if(i.le.nap+inzt)then
p17(i,3)=pt1 ! target particle is proton
p17(i,4)=et1
else
p17(i,3)=pt2 ! target particle is neutron
p17(i,4)=et2
endif
enddo
c-------------------------- Momentum Initializing --------------------------
c-------------------------------------------------------------------------------
c calculate the velocity of the CM of collision system in LAB or
c in nucleon-nucleon CM system
bst(1)=p17(1,1)*nap+p17(nap+1,1)*nat
bst(2)=p17(1,2)*nap+p17(nap+1,2)*nat
bst(3)=p17(1,3)*nap+p17(nap+1,3)*nat
bst(4)=p17(1,4)*nap+p17(nap+1,4)*nat
bst(1)=-bst(1)/bst(4)
bst(2)=-bst(2)/bst(4)
bst(3)=-bst(3)/bst(4)
c-------------------------------------------------------------------------------
c----------------------- Local Variable Initializing -----------------------
do i=1,napt
tp(i)=0.
enddo
n=0
nbe=0 ! 080104
naf=0 ! 080104
nsa=0
nbh=0 ! 300623 Lei
idi=0
idio=0
ndiq=0
npt=0
ifcom=0 ! 220110
ishp=0
tau=0.
nctl=0
lc=0
tc=0.
tw=0.
numb=0
numbs=0
c----------------------- Local Variable Initializing -----------------------
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c----------------------- Particle Properties Giving ------------------------
c '1 -> |nzp|' are projectile protons or lepton, '|nzp|+1 -> nap'
c are projectile neutrons; 'nap+1 -> nap+nzt' are targer protons,
c the rest are target nuctrons in 'pyjets' after nuclear initiation above
n=napt
do i=1,n
k(i,1)=1
k(i,2)=2112
p(i,5)=pmas(pycomp(2112),1)
c160224 Lei For NN, NA(AN) and AA.
if( (i.le.abs(nzp).and.ipden.lt.2).or.(i.gt.nap .and. i.le.nap+
c nzt) )then ! 060813 120214
k(i,2)=2212
p(i,5)=pmas(pycomp(2212),1)
c061123 Lei For l+N & lbar + N
elseif( i.le.nap .AND.
& (ipden.ge.11 .AND. ipden.le.16 .AND. ABS(nzp).eq.1) )then
K(i,2) = SIGN( ipden, -nzp )
P(i,5) = PMAS( PYCOMP(ipden), 1 )
c061123 Lei
endif
do j=1,3
p(i,j)=p17(i,j)
v(i,j)=c17(i,j)
enddo
p(i,4)=p17(i,4)
v(i,4)=tp(i)
enddo
500 continue ! 031103
c v, vbh and vsa arraies are the position four vector
c note: for v etc., we do not take care of their fifth component
c for array k, we take care of only first three components
c----------------------- Particle Properties Giving ------------------------
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c--------------------- CMS Boost & Lorentz Contraction ---------------------
c boost PYJETS into cms of initial nucleus-nucleus collision system
c from lab or initial nucleon-nucleon cms system.
c call pyrobo(1,n,0.0,0.0,bst(1),bst(2),bst(3))
c Lorentz contract
bzp3=0.
bzp4=0.
bzt3=0.
bzt4=0.
do i=1,nap
bzp3=bzp3+p(i,3)
bzp4=bzp4+p(i,4)
enddo
do i=nap+1,napt
bzt3=bzt3+p(i,3)
bzt4=bzt4+p(i,4)
enddo
bzp=bzp3/bzp4
bzt=bzt3/bzt4
gamp=1.d0/dsqrt(dmax1(1.d-20,(1.0d0-bzp*bzp)))
c060813 120214 no Lorentz contraction for incident lepton
if(ipden.ge.2)gamp=1. ! 060813 120214
gamt=1.d0/dsqrt(dmax1(1.d-20,(1.0d0-bzt*bzt)))
c try no lorentz contract for target
c gamt=1.
do i=1,nap
c17(i,3)=c17(i,3)/gamp
v(i,3)=v(i,3)/gamp
enddo
do i=nap+1,napt
c17(i,3)=c17(i,3)/gamt
v(i,3)=v(i,3)/gamt
enddo
c180724 Lei
c Positions two particles at ( b/2, 0, -20 ) and ( -b/2, 0, 20 ).
c 20 fm is just a large enough distance.
do i=1,nap,1
V(i,3) = V(i,3) - 20D0
end do
do i = nap+1, napt, 1
V(i,3) = V(i,3) + 20D0
end do
c180724 Lei
c--------------------- CMS Boost & Lorentz Contraction ---------------------
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c---------------------- Particle Filtering & Ordering ----------------------
c260223
c if(iMode.eq.2.or.iMode.eq.3)then
c filter out those kind of particles wanted to study and make
c the order of proton, neutron, ... (cf. 'filt')
call filt
c060813 120214
c since lepton was moved to last position after calling filt, one has to
c remove it to the fist position
if(ipden.ge.2)call ltof(n)
c060813 120214
c endif
c260223
c161021 'pyjets' to 'sa2'
nsa=n
do j=1,5
do i=1,n
ksa(i,j)=k(i,j)
psa(i,j)=p(i,j)
vsa(i,j)=v(i,j)
enddo
enddo
do i=1,n
ishp(i)=1
enddo
c260223
if(iMode.eq.2.or.iMode.eq.3.or.iMode.eq.4)then ! 280524
do m=1,kfmax
numb(m)=numbs(m)
enddo
endif
c260223
c---------------------- Particle Filtering & Ordering ----------------------
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c note: particle list is composed of the arraies in common block
c 'sa2', the array 'ishp' in common block 'wz', the array 'tau' in
c common block 'sa4', and the array 'numb' in common block 'sa5'
time=time_ini ! 081010
irecon=0
csa011223
c upto here the initial nucleon list is available
c140223 Lei full_events_history of OSC1999A
call oscar(win,0)
c calculate the position for the center of mass of the
c non-freeze-out system. The distance of a particle, when checking
c is it freezing out or not, is measured with respect to this center
call copl(time)
c creat the initial collision list, note: be sure that the initial
c collision list must not be empty
call ctlcre(lc,tc,tw)
csa011223
c upto here the initial NN collision time list is available
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c------------------------ System Time Initializing -------------------------
csa time origin is set at the time of first NN collision
c find out colli. pair with least colli. time
call find(icp,tcp,lc,tc,tw,0)
if(icp.eq.0)stop 'initial collision list is empty' !
time=tcp
c070417 perform classical Newton motion in Lab. system for all particles
call his(time,lc,tc,tw,istop)
do ij=1,nsa
vsa(ij,4)=0.
enddo
c070417 move origin of time to collision time of first nucleon-nucleon collision
do ij=1,nctl
tc(ij)=tc(ij)-time+1.e-5
enddo
time=time_ini ! 081010
call copl(time)
c------------------------ System Time Initializing -------------------------
c-------------------------------------------------------------------------------
400 continue
c-------------------------------------------------------------------------------
c------------------------- Collision Implementing --------------------------
csa011223
c loop over implementing NN (hh) collision, updating hadron list, and
c updating collision time list untill the collision time list is empty.
c It is equivalent to implementing a nucleus-nucleus collision
call scat(time,lc,tc,tw,win,parp2,psno,ijk,ipau,irecon,
c gamt,iMode,decpro,i_tune,time_par,time_had) ! 021207 260223 280524
c300623 Lei Added i_tune 280823 Lei Removed parp22
if(ijk.eq.1)return
time_ini=time ! 081010
800 continue
c 'saf' to 'pyjets'
c180520 if(adj1(40).ne.5)call tran_saf ! 140414
if(mstptj.eq.0)call tran_saf ! 140414 180520 230722
n00=n ! 220110
c220110 n00: 'largest line number' in 'pyjets'
c220110 partons above n00 appear after inelastic collision
c 'sa2' to 'sbh'
c190224 'sbh' stores spectators, hadrons & leptons from the diffractive & the
c special sub-processes, e.g. NRQCD onia, or hadron beam remnants, or
c B-framework.
nbh=0
if(nsa.ge.1)then
nbh=nsa
do i2=1,5
do i1=1,nsa
kbh(i1,i2)=ksa(i1,i2)
pbh(i1,i2)=psa(i1,i2)
vbh(i1,i2)=vsa(i1,i2)
enddo
enddo
endif
c P(N,5)=SQRT(MAX(-P(N,1)**2-P(N,2)**2-P(N,3)**2+P(N,4)**2,0.0))
c P(N-1,5)=SQRT(MAX(-P(N-1,1)**2-P(N-1,2)**2-P(N-1,3)**2
c & +P(N-1,4)**2,0.0))
c call pyboro(1,n,0.0,0.0,-bst(1),-bst(2),-bst(3))
c boost PYJETS back to lab or nucleon-nucleon cms system.
c------------------------- Collision Implementing --------------------------
c-------------------------------------------------------------------------------
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine sysini(win) ! 060813
c give the initial values to quantities needed in calculation
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
IMPLICIT INTEGER(I-N)
INTEGER PYK,PYCHGE,PYCOMP
PARAMETER (KSZJ=80000)
COMMON/PYCIDAT1/KFACOT(100),DISDET(100),ISINELT(600)
common/sa5/kfmax,kfaco(100),numb(100),numbs(100),non5,
c disbe(100,100)
common/count/isinel(600)
COMMON/PYCIDAT2/KFMAXT,NONCI2,PARAM(20),WEIGH(600)
common/sa6/kfmaxi,nwhole
common/papr/t0,sig,dep,ddt,edipi,epin,ecsnn,ekn,ecspsn,ecspsm
c ,rnt,rnp,rao,rou0,vneu,vneum,ecsspn,ecsspm,ecsen ! 060813
common/syspar/ipden,itden,suppm,suptm,suppc,suptc,r0p,r0t,
c nap,nat,nzp,nzt,pio
common/sa10/csnn,cspin,cskn,cspipi,cspsn,cspsm,rcsit,ifram,
& iabsb,iabsm,non10,ajpsi,csspn,csspm,csen ! 060813
common/sa25/i_inel_proc,i_time_shower,para1_1,para1_2 ! 221203 250204
c cspsn : total cross section of J/psi (psi') + N
c cspsm : total cross section of J/psi (psi') + meson
c iabsb = 0 : without J/psi (psi') + baryon
c = 1 : with J/Psi (psi') + baryon
c iabsm = 0 : without J/psi (psi') + meson
c = 1 : with J/psi (psi') + meson
anat=nat
anap=nap
PARAM(1)=para1_1 ! 250204 200504
c rou0=PARAM(11)
c considering the nucleus as a sphere with radii rnt for target
c and rnp for projectile.
c rnt=(3.*anat/(4.*3.1415926*rou0))**(0.33333)
c rnp=(3.*anap/(4.*3.1415926*rou0))**(0.33333)
rp00=1.12 ! 1.05 to 1.12 070613
rt00=1.12 ! 1.05 to 1.12 070613
c070613 if(nap.gt.16)rp00=1.16*(1-1.16*anap**(-0.666666)) !rp00=1.122 (nat=208)
c070613 if(nat.gt.16)rt00=1.16*(1-1.16*anat**(-0.666666)) ! rt00=1.12 (nat=197)
c210924 Lei
c Uses PDG RPP2024 charge radius of proton 0.8409 +- 0.0004 fm.
c (PDG RPP2024 magnetic radius of neutron 0.864 +0.009 -0.008 fm)
c if(itden.eq.0)rnt=rt00*anat**(0.33333) ! 310805
if(itden.eq.0) rnt = 0.841
c210924 Lei
if(itden.eq.1)rnt=rt00*anat**(0.33333) ! +0.54 160511
if(nat.eq.2 .and. nzt.eq.1)rnt=4.0 ! 2.60 2.095 1.54 2603141
c060813 120214 if(ipden.eq.2)rnt=0.5 ! lepton
c210924 Lei
c Uses PDG RPP2024 charge radius of proton 0.8409 +- 0.0004 fm.
c (PDG RPP2024 magnetic radius of neutron 0.864 +0.009 -0.008 fm)
c if(ipden.eq.0)rnp=rp00*anap**(0.33333) ! 310805
if(ipden.eq.0) rnp = 0.841
c210924 Lei
if(ipden.eq.1)rnp=rp00*anap**(0.33333) ! +0.54 160511
if(ipden.ge.2)rnp=0.5 ! lepton ! 060813 120214
if(nap.eq.2 .and. nzp.eq.1)rnp=4.0 ! 2.60 2.095 1.54
rou0=3./4./3.1416*anat/(rnt*rnt*rnt) ! 310805
r0p=rnp
r0t=rnt
C set initial values to some quantities
c in the program the x-sections are given in a unit of fm^2 ! 060813
csnn=PARAM(1)*0.1
c250423 cspin=PARAM(2)*0.1
cspin=csnn*0.66666 ! 250423
c250423 0.66666=6/9, estimated by additive quark model (arXiv:2203.11061)
c250423 cskn=PARAM(3)*0.1
cskn=csnn*0.28444 ! 250423
c250423 0.8444=1.6*1.6/9, estimated by additive quark model
c250423 cspipi=PARAM(4)*0.1
cspipi=csnn*0.44444
c250423 0.44444=4/9, estimated by additive quark model
c250423 cspsn=PARAM(13)*0.1
cspsn=csnn*0.13333
c250423 0.13333=0.4*3/9, estimated by additive quark model
c250423 cspsm=PARAM(14)*0.1
cspsm=csnn*0.08888
c250423 0.08888=0.4*2/9, estimated by additive quark model
c250423 csspn=PARAM(15)*0.1
csspn=cspsn ! 250423
c250423 csspm=PARAM(16)*0.1
csspm=cspsm ! 250423
c060813 120214
if(ipden.ge.2)then
if(ifram.eq.0)then
ept=sqrt(win*win+0.938*0.938)
rots=sqrt((ept+0.938)*(ept+0.938)-win*win)
endif
if(ifram.eq.1)rots=win
call crosep(rots,csen) ! temporary using e^-p total x-section
c if(nzp.lt.0)call crosep(rots,csen) ! e^-p total x-section
c if(nzp.ge.0)call crosepp(rots,csen) ! e^+p total x-section
csen=csen*0.1
endif
c060813 120214
c largest collision distance between two colliding particles.
edipi=dsqrt(cspipi/3.1416)
epin=dsqrt(cspin/3.1416)
ekn=dsqrt(cskn/3.1416)
ecsnn=dsqrt(csnn/3.1416)
ecspsn=dsqrt(cspsn/3.1416)
ecspsm=dsqrt(cspsm/3.1416)
ecsspn=dsqrt(csspn/3.1416)
ecsspm=dsqrt(csspm/3.1416)
ecsen=sqrt(csen/3.1416) ! 060813
anp=nap**.3333
ant=nat**.3333
do ia=1,2
if(ia.eq.1)napt=nap
if(ia.eq.2)napt=nat
if(napt.lt.27)then
alpt=0.47
elseif(napt.ge.27.and.napt.le.108)then
alpt=0.488
else
alpt=0.54
endif
if(napt.eq.27)then
alpt=0.478
elseif(napt.eq.28)then
alpt=0.48
elseif(napt.eq.32)then
alpt=0.49
elseif(napt.eq.56)then
alpt=0.49
elseif(napt.eq.64)then
alpt=0.49
elseif(napt.eq.108)then
alpt=0.495
elseif(napt.eq.184)then
alpt=0.53
elseif(napt.eq.197)then
alpt=0.54
elseif(napt.eq.207)then
alpt=0.545
elseif(napt.eq.238)then
alpt=0.55