EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
example_5.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file example_5.f
1 
2 C****************************************************************************
3 C The program for checking conservation laws fulfillment in
4 C the HIJING model.
5 C Written by V.Uzhinsky, CERN, Oct. 2003
6 C***************************************************************************
7 
8  CHARACTER frame*8,proj*8,targ*8
9 
10  common/himain1/ natt,eatt,jatt,nt,np,n0,n01,n10,n11
11  SAVE /himain1/
12 
13  common/himain2/katt(130000,4),patt(130000,4)
14  SAVE /himain2/
15 
16 C ********information of produced particles
17 C
18  CHARACTER *1 key
19  CHARACTER *12 fname
20 
21  parameter(nwpawc=50000)
22  common/pawc/hmemory(nwpawc)
23 
24  common/ranseed/nseed
25  SAVE /ranseed/
26  nseed=0
27 
28  write(6,*)'====================================================='
29  write(6,*)'The program for checking conservation laws fulfillment'
30  write(6,*)' in the HIJING model '
31  write(6,*)' Only hadronic (hh) interactions are considered '
32  write(6,*)'====================================================='
33  write(6,*)' You can work in Lab. or CM systems '
34  write(6,*)' Would you like to use CM system? Y - Yes, N - No? '
35 
36  read(5,1001) key
37 1001 format(a1)
38 
39  if(key.eq.'Y'.or.key.eq.'y') then
40  write(6,*)'CM system is used --------------------------------'
41  frame='CMS'
42  else
43  write(6,*)'LAB system is used -------------------------------'
44  frame='LAB'
45  endif
46 
47  write(6,*)'Enter the corresponding energy (GeV)'
48 
49  read(5,*) efrm
50 
51  write(6,*)
52  write(6,*)'Enter a type of the projectile particle'
53  write(6,*)
54  write(6,*)' P proton, PBAR anti-proton,'
55  write(6,*)' N neutron, NBAR anti-neutron,'
56  write(6,*)' PI+ - positive pion, PI- negative pion,'
57  write(6,*)' K+ positive kaon, K- negative kaon'
58 
59  read(5,1002) proj
60 1002 format(a8)
61 
62  if(proj.ne.'A') then
63  iap=1
64  if(proj.eq.'P' ) izp= 1
65  if(proj.eq.'PBAR') izp=-1
66  if(proj.eq.'N' ) izp= 0
67  if(proj.eq.'NBAR') izp= 0
68  if(proj.eq.'PI+' ) izp= 1
69  if(proj.eq.'PI-' ) izp=-1
70  if(proj.eq.'K+' ) izp= 1
71  if(proj.eq.'K-' ) izp=-1
72  else
73  write(6,*)
74  write(6,*)'Enter mass number and charge of the proj. nucleus'
75  read(5,*) iap, izp
76  endif
77 
78  write(6,*)
79  write(6,*)'Enter a type of the target particle (same notations)'
80  read(5,1002) targ
81 
82  if(targ.ne.'A') then
83  iat=1
84  if(targ.eq.'P' ) izt= 1
85  if(targ.eq.'PBAR') izt=-1
86  if(targ.eq.'N' ) izt= 0
87  if(targ.eq.'NBAR') izt= 0
88  if(targ.eq.'PI+' ) izt= 1
89  if(targ.eq.'PI-' ) izt=-1
90  if(targ.eq.'K+' ) izt= 1
91  if(targ.eq.'K-' ) izt=-1
92  else
93  write(6,*)
94  write(6,*)'Enter mass number and charge of the target nucleus'
95  read(5,*) iat, izt
96  endif
97 
98  write(6,*)
99  write(6,*)'Enter number of events'
100  read(5,*) n_events
101 
102  write(6,*)'Enter FILENAME for HBOOK output'
103  read(5,1003) fname
104 1003 format(a12)
105 
106 c------------------------------------------------------------------
107  CALL hlimit(nwpawc)
108 
109  efrml=efrm*(9./10.)
110  efrmh=efrm*(11./10.)
111  CALL hbook1(10,' Energy conservation',100,efrml,efrmh,0.)
112 
113  pxl=-5.
114  pxh= 5.
115  CALL hbook1(20,' Px momentum conservation',100,pxl,pxh,0.)
116 
117  pyl=-5.
118  pyh= 5.
119  CALL hbook1(30,' Py momentum conservation',100,pyl,pyh,0.)
120 
121  if(frame.eq.'CMS') then
122  pzl=-5.
123  pzh= 5.
124  else
125  if(proj.eq.'P'.or.proj.eq.'PBAR'.or.
126  & proj.eq.'N'.or.proj.eq.'NBAR'.or.proj.eq.'A') then
127  plab=sqrt(efrm**2-0.88)
128  endif
129 
130  if(proj.eq.'PI+'.or.proj.eq.'PI-') then
131  plab=sqrt(efrm**2-0.0196)
132  endif
133 
134  if(proj.eq.'K+'.or.proj.eq.'K-') then
135  plab=sqrt(efrm-0.25)
136  endif
137 
138  pzl=plab*(3./4.)
139  pzh=plab*(5./4.)
140  endif
141 
142  CALL hbook1(40,' Pz momentum conservation',100,pzl,pzh,0.)
143 
144  CALL hbook1(50,' Charge conservation' ,100,-5.,5.,0.)
145 
146  CALL hbook1(60,' Baryon number conservation',100,-5.,5.,0.)
147 
148  CALL hbook1(70,' Lepton number conservation',100,-5.,5.,0.)
149 
150  CALL hbook1(80,' Azimuthal isotropy' ,100,-3.2,3.2,0.)
151 C----------------------------------------------------------------
152 
153  CALL hijset(efrm,frame,proj,targ,iap,izp,iat,izt)
154 C ********Initialize HIJING
155 
156  WRITE(6,*)' Simulation of interactions with'
157  WRITE(6,*)
158  WRITE(6,*)' Proj = ',proj,' and Targ = ',targ
159  WRITE(6,*)' IAP =',iap ,' IAT =',iat
160  WRITE(6,*)' IZP =',izp ,' IZT =',izt
161  WRITE(6,*)
162  WRITE(6,*)' Reference frame - ',frame
163  WRITE(6,*)' ENERGY ',efrm,' GeV'
164  WRITE(6,*)' Number of generated events -',n_events
165  WRITE(6,*)
166 
167 
168 
169  bmin=0.0
170  bmax=0.0
171 
172  DO 2000 i_event=1,n_events
173 C
174  WRITE(6,*)' Event # ',i_event,' ------------------'
175 C
176  CALL hijing(frame,bmin,bmax)
177 C
178  sum_e =0.
179  sum_px=0.
180  sum_py=0.
181  sum_pz=0.
182  sum_q =0.
183  sum_b =0.
184  sum_l =0.
185 
186  DO 3000 i=1,natt
187 
188  if(katt(i,2).eq. 0) go to 3000 ! reject non-interacting projectile
189  if(katt(i,2).eq. 1) go to 3000 ! reject elastic scattering
190  if(katt(i,2).eq.10) go to 3000 ! reject non-interacting target
191  if(katt(i,2).eq.11) go to 3000 ! reject elastic scattering
192 
193  ich=luchge(katt(i,1))/3
194 
195  amass=ulmass(katt(i,1))
196 
197  sum_e =sum_e + patt(i,4)
198  sum_px=sum_px+ patt(i,1)
199  sum_py=sum_py+ patt(i,2)
200  sum_pz=sum_pz+ patt(i,3)
201  sum_q =sum_q + ich
202 
203  if(iabs(katt(i,1)).gt.2000) then
204  sum_b=sum_b+iabs(katt(i,1))/katt(i,1)
205  endif
206 
207  if((iabs(katt(i,1)).gt.10.).and.
208  & (iabs(katt(i,1)).lt.20.)) then
209  sum_l=sum_l+iabs(katt(i,1))/katt(i,1)
210  endif
211 
212  phi=ulangl(patt(i,1),patt(i,2))
213  CALL hf1(80,phi,1.)
214 
215 3000 CONTINUE
216 
217  CALL hf1(10,sum_e ,1.)
218  CALL hf1(20,sum_px,1.)
219  CALL hf1(30,sum_py,1.)
220  CALL hf1(40,sum_pz,1.)
221  CALL hf1(50,sum_q ,1.)
222  CALL hf1(60,sum_b ,1.)
223  CALL hf1(70,sum_l ,1.)
224 2000 CONTINUE
225 
226  CALL hrput(0,fname,'N')
227  END
228 
229  FUNCTION ran(NSEED)
230  ran=rlu(nseed)
231  RETURN
232  END