Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
example_3.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file example_3.f
1 C****************************************************************************
2 C Program # 2 from Comp. Phys. Commun. 83 (1994) 307
3 C by M. Gyulassy and X-.N. Wang
4 C Modified by V.Uzhinsky, CERN, Oct. 2003
5 C***************************************************************************
6 
7  CHARACTER frame*8,proj*8,targ*8
8 
9  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
10  SAVE /hiparnt/
11 
12  common/himain1/natt,eatt,jatt,nt,np,n0,n01,n10,n11
13  SAVE /himain1/
14 
15  common/himain2/katt(130000,4),patt(130000,4)
16  SAVE /himain2/
17 C
18  dimension gb(101), xb(101), dndp(50)
19 
20  common/ranseed/nseed
21  SAVE /ranseed/
22 
23  nseed=0
24 
25  do i=1,50
26  dndp(i)=0.
27  enddo
28 
29  frame='CMS'
30 
31  write(6,*)'===================================================='
32  write(6,*)' Calculation of transverse momentum distribution '
33  write(6,*)' of charged pions in minimum bias '
34  write(6,*)' hh-, hA- and AA-collisions '
35  write(6,*)' Calculation will be performed in CM System '
36  write(6,*)'===================================================='
37 
38  write(6,*)
39  write(6,*)'Enter the energy per NN-collision (GeV)'
40  read(5,*) efrm
41 
42  write(6,*)
43  write(6,*)'Enter a type of the "projectile" particle'
44  write(6,*)
45  write(6,*)' P proton, PBAR anti-proton,'
46  write(6,*)' N neutron, NBAR anti-neutron,'
47  write(6,*)' PI+ - positive pion, PI- negative pion,'
48  write(6,*)' K+ positive kaon, K- negative kaon'
49  write(6,*)
50  write(6,*)' A - nucleus --------------------------'
51 
52  read(5,1) proj
53 1 format(a8)
54 
55  if(proj.ne.'A') then
56  iap=1
57  if(proj.eq.'P' ) izp= 1
58  if(proj.eq.'PBAR') izp=-1
59  if(proj.eq.'N' ) izp= 0
60  if(proj.eq.'NBAR') izp= 0
61  if(proj.eq.'PI+' ) izp= 1
62  if(proj.eq.'PI-' ) izp=-1
63  if(proj.eq.'K+' ) izp= 1
64  if(proj.eq.'K-' ) izp=-1
65  else
66  write(6,*)
67  write(6,*)'Enter mass number and charge of the proj. nucleus'
68  read(5,*) iap, izp
69  endif
70 
71  write(6,*)
72  write(6,*)'Enter a type of the "target" particle (same notations)'
73  read(5,1) targ
74 
75  if(targ.ne.'A') then
76  iat=1
77  if(targ.eq.'P' ) izt= 1
78  if(targ.eq.'PBAR') izt=-1
79  if(targ.eq.'N' ) izt= 0
80  if(targ.eq.'NBAR') izt= 0
81  if(targ.eq.'PI+' ) izt= 1
82  if(targ.eq.'PI-' ) izt=-1
83  if(targ.eq.'K+' ) izt= 1
84  if(targ.eq.'K-' ) izt=-1
85  else
86  write(6,*)
87  write(6,*)'Enter mass number and charge of the target nucleus'
88  read(5,*) iat, izt
89  endif
90 
91  write(6,*)
92  write(6,*)'Enter number of events per each value of impact',
93  ,' parameter (e.a. 10)'
94 
95  read(5,*) n_event
96 
97 C....initialize HIJING for requested interactions
98  CALL hijset(efrm,frame,proj,targ,iap,izp,iat,izt)
99 
100  write(6,*)' Simulation of interactions with'
101  write(6,*)
102  write(6,*)' Proj = ',proj,' and Targ = ',targ
103  write(6,*)' IAP =',iap ,' IAT =',iat
104  write(6,*)' IZP =',izp ,' IZT =',izt
105  write(6,*)
106  write(6,*)' Reference frame - ',frame
107  write(6,*)' ENERGY ',efrm,' GeV'
108  write(6,*)' Number of generated events per B interval -',n_event
109  write(6,*)
110 
111 C....set BMIN=0 and BMAX=R_A+R_B
112 
113  bmin=0.0
114  bmax=hipr1(34)+hipr1(35)
115 
116 C....calculate the Glauber probability and its integrated value:
117 
118  dip=(bmax-bmin)/100.0
119  gbtot=0.0
120  DO 100 i=1,101
121  xb(i)=bmin+(i-1)*dip
122  ov=profile(xb(i))
123  gb(i)=xb(i)*(1.0-exp(-hint1(12)*ov))
124  gbtot=gbtot+gb(i)
125 100 CONTINUE
126  write(6,*)'Inelastic X-section (mb) ',gbtot*dip*10.*6.28
127 
128  pause
129 C....generating N_EVENT for each of 100 impact parameter intervals:
130 
131  nont=0
132  gnorm=gbtot
133 
134  DO 300 ib=1,100
135  b1=xb(ib)
136  b2=xb(ib+1)
137 
138 C....normalized Glauber probability:
139 
140  w_gb=(gb(ib)+gb(ib+1))/2.0/gbtot
141 
142  DO 200 ie=1,n_event
143  CALL hijing(frame,b1,b2)
144 
145 C....count number of events without any interaction
146 C....and renormalize the total Glauber prabability:
147 
148  IF(np+nt .EQ. 0) THEN
149  nont=nont+1
150  gnorm=gnorm-gb(ib)/float(n_event)
151  go to 200
152  ENDIF
153 
154 C....calculate pt distribution of charged pions:
155 
156  DO 150 k=1,natt
157 C....select charged pions only:
158  IF(abs(katt(k,1)) .NE. 211) go to 150
159 
160 C....calculate pt:
161  ptr=sqrt(patt(k,1)**2+patt(k,2)**2)
162 
163 C....calculate pt distribution and weight with normalized
164 C....Glauber probability to get minimum bias results:
165 
166  IF(ptr .GE. 10.0) go to 150
167  ipt=1+ptr/0.2
168  dndp(ipt)=dndp(ipt)+w_gb/float(n_event)/0.2
169 150 CONTINUE
170 200 CONTINUE
171 300 CONTINUE
172 
173 C....renormalize the distribution by the renormalized Glauber
174 C....probability which excludes the events without any interaction:
175 
176  IF(nont.NE.0) THEN
177  DO 400 i=1,50
178  dndp(i)=dndp(i)*gbtot/gnorm
179  if(dndp(i).ne.0.) write(6,350)0.2*(i-1),dndp(i)
180 350 format(1x,f5.1,2x,e11.4)
181 400 CONTINUE
182  ENDIF
183 
184  stop
185  END
186 
187  FUNCTION ran(NSEED)
188  ran=rlu(nseed)
189  RETURN
190  END