EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
example_6.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file example_6.f
1 
2 C****************************************************************************
3 C The program for calculation of the general properties of
4 C interactions in 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  common/himain2/katt(130000,4),patt(130000,4)
12 
13 C ********information of produced particles
14 C
15 
16  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
17  SAVE /hiparnt/
18 
19  CHARACTER *1 key
20  CHARACTER *12 fname
21 
22  parameter(nwpawc=1000000)
23  common/pawc/hmemory(nwpawc)
24 
25  common/ranseed/nseed
26  SAVE /ranseed/
27  nseed=0
28 
29  write(6,*)'====================================================='
30  write(6,*)'The program for calculation of general properties of '
31  write(6,*)' interactions in the HIJING model. '
32  write(6,*)'====================================================='
33  write(6,*)
34  write(6,*)' You can work in Lab. or CM systems '
35  write(6,*)' Would you like to use CM system? Y - Yes, N - No? '
36 
37  read(5,1001) key
38 1001 format(a1)
39 
40  if(key.eq.'Y'.or.key.eq.'y') then
41  write(6,*)'CM system is used --------------------------------'
42  frame='CMS'
43  else
44  write(6,*)'LAB system is used -------------------------------'
45  frame='LAB'
46  endif
47 
48  write(6,*)'Enter the corresponding energy per NN collisions (GeV)'
49 
50  read(5,*) efrm
51 
52  write(6,*)
53  write(6,*)'Enter a type of the projectile particle'
54  write(6,*)
55  write(6,*)' P proton, PBAR anti-proton,'
56  write(6,*)' N neutron, NBAR anti-neutron,'
57  write(6,*)' PI+ - positive pion, PI- negative pion,'
58  write(6,*)' K+ positive kaon, K- negative kaon'
59  write(6,*)' A - nucleus --------------------------'
60 
61 
62  read(5,1002) proj
63 1002 format(a8)
64 
65  if(proj.ne.'A') then
66  iap=1
67  if(proj.eq.'P' ) izp= 1
68  if(proj.eq.'PBAR') izp=-1
69  if(proj.eq.'N' ) izp= 0
70  if(proj.eq.'NBAR') izp= 0
71  if(proj.eq.'PI+' ) izp= 1
72  if(proj.eq.'PI-' ) izp=-1
73  if(proj.eq.'K+' ) izp= 1
74  if(proj.eq.'K-' ) izp=-1
75  else
76  write(6,*)
77  write(6,*)'Enter mass number and charge of the proj. nucleus'
78  read(5,*) iap, izp
79  endif
80 
81  write(6,*)
82  write(6,*)'Enter a type of the target particle (same notations)'
83  read(5,1002) targ
84 
85  if(targ.ne.'A') then
86  iat=1
87  if(targ.eq.'P' ) izt= 1
88  if(targ.eq.'PBAR') izt=-1
89  if(targ.eq.'N' ) izt= 0
90  if(targ.eq.'NBAR') izt= 0
91  if(targ.eq.'PI+' ) izt= 1
92  if(targ.eq.'PI-' ) izt=-1
93  if(targ.eq.'K+' ) izt= 1
94  if(targ.eq.'K-' ) izt=-1
95  else
96  write(6,*)
97  write(6,*)'Enter mass number and charge of the target nucleus'
98  read(5,*) iat, izt
99  endif
100 
101  write(6,*)
102  write(6,*)'Enter number of events'
103  read(5,*) n_events
104 
105  write(6,*)'Enter FILENAME for HBOOK output'
106  read(5,1003) fname
107 1003 format(a12)
108 
109 c------------------------------------------------------------------
110  CALL hijset(efrm,frame,proj,targ,iap,izp,iat,izt)
111 C ********Initialize HIJING
112 
113  WRITE(6,*)' Simulation of interactions with'
114  WRITE(6,*)
115  WRITE(6,*)' Proj = ',proj,' and Targ = ',targ
116  WRITE(6,*)' IAP =',iap ,' IAT =',iat
117  WRITE(6,*)' IZP =',izp ,' IZT =',izt
118  WRITE(6,*)
119  WRITE(6,*)' Reference frame - ',frame
120  WRITE(6,*)' ENERGY ',efrm,' GeV'
121  WRITE(6,*)' Number of generated events -',n_events
122  WRITE(6,*)
123 
124  bmin=0.
125  bmax=hipr1(34)+hipr1(35)
126 
127 c------------------------------------------------------------------
128 c----------------- Charged particles ------------------------------
129 c------------------------------------------------------------------
130  CALL hlimit(nwpawc)
131 
132  CALL hbook1(2,'Impact parameter distribution',
133  ,100,bmin,bmax,0.)
134 
135  nu_max=min(iap*iat+1,2000)
136  CALL hbook1(4,'Nu-distribution',
137  ,100,-0.5,float(nu_max),0.)
138 
139  CALL hbook1(6,'Distribution on the number of wounded A nucleons',
140  ,iap+1,-0.5,float(iap)+0.5,0.)
141 
142  CALL hbook1(8,'Distribution on the number of wounded B nucleons',
143  ,iat+1,-0.5,float(iat)+0.5,0.)
144 
145  if(proj.eq.'A'.or.targ.eq.'A') then
146  ch_max=5*nu_max
147  m_neg=ch_max/2.
148  else
149  ch_max=99.5
150  m_neg=ch_max
151  endif
152  CALL hbook1(10,' Charged particle multiplicity distribution',
153  ,100,-0.5,ch_max,0.)
154 
155  if(frame.eq.'CMS') then
156  eta_l=-15.
157  eta_h=+15.
158  else
159  eta_l=-2.0
160  eta_h=alog(2.*efrm)+2.
161  endif
162  nbineta=(eta_h-eta_l)/0.2
163 
164  CALL hbook1(20,' Charged particle pseudo-rapidity distribution',
165  ,nbineta,eta_l,eta_h,0.)
166 
167  if(frame.eq.'CMS') then
168  y_l=-alog(efrm)-2.
169  y_h=+alog(efrm)+2.
170  else
171  y_l=-2.0
172  y_h=alog(2.*efrm)+2.
173  endif
174  nbiny=(y_h-y_l)/0.2
175 
176  CALL hbook1(30,' Charged particle rapidity distribution',
177  ,nbiny,y_l,y_h,0.)
178 
179  CALL hbook1(40,' Charged particle Pt distribution',
180  ,100,0.,10.,0.)
181 
182  if(frame.eq.'CMS') then
183  e_l= 0.
184  e_h=efrm/2.
185  else
186  e_l= 0.
187  e_h=efrm
188  endif
189  nbine=(e_h-e_l)/0.5
190 
191  CALL hbook1(50,' Charged particle energy distribution',
192  ,nbine,e_l,e_h,0.)
193 
194  CALL hbook1(60,' Charged particle Cos(Theta) distribution',
195  ,40,-1.,1.,0.)
196 
197  CALL hbook1(70,' Charged particle Phi distribution',
198  ,180,0.,180.,0.)
199 
200  CALL hbook2(80,' Phi - Eta correlation of charged particles',
201  ,nbineta,eta_l,eta_h,180,-180.,180.,0.)
202 
203  CALL hbook1(90,' Particle composition (IDs)',
204  ,6600,-3300.,3300.,0.)
205 
206 C----------------------------------------------------------------
207 c--------------- Negative charged particles ---------------------
208 c----------------------------------------------------------------
209  CALL hbook1(110,' Negative particle multiplicity distribution',
210  ,100,-0.5,float(m_neg),0.)
211 
212  CALL hbook1(120,' Negative particle pseudo-rapidity distribution',
213  ,nbineta,eta_l,eta_h,0.)
214 
215  CALL hbook1(130,' Negative particle rapidity distribution',
216  ,nbiny,y_l,y_h,0.)
217 
218  CALL hbook1(140,' Negative particle Pt distribution',
219  ,100,0.,10.,0.)
220 
221  CALL hbook1(150,' Negative particle energy distribution',
222  ,nbine,e_l,e_h,0.)
223 
224  CALL hbook1(160,' Negative particle Cos(Theta) distribution',
225  ,40,-1.,1.,0.)
226 
227  CALL hbook1(170,' Negative particle Phi distribution',
228  ,180,0.,180.,0.)
229 
230  CALL hbook2(180,' Phi - Eta correlation of negative particles',
231  ,nbineta,eta_l,eta_h,180,-180.,180.,0.)
232 
233 c----------------------------------------------------------------
234 c----------------Positive charged particles ---------------------
235 c----------------------------------------------------------------
236  CALL hbook1(210,' Positive particle multiplicity distribution',
237  ,100,-0.5,float(m_neg),0.)
238 
239  CALL hbook1(220,' Positive particle pseudo-rapidity distribution',
240  ,nbineta,eta_l,eta_h,0.)
241 
242  CALL hbook1(230,' Positive particle rapidity distribution',
243  ,nbiny,y_l,y_h,0.)
244 
245  CALL hbook1(240,' Positive particle Pt distribution',
246  ,100,0.,10.,0.)
247 
248  CALL hbook1(250,' Positive particle energy distribution',
249  ,nbine,e_l,e_h,0.)
250 
251  CALL hbook1(260,' Positive particle Cos(Theta) distribution',
252  ,40,-1.,1.,0.)
253 
254  CALL hbook1(270,' Positive particle Phi distribution',
255  ,180,0.,180.,0.)
256 
257  CALL hbook2(280,' Phi - Eta correlation of positive particles',
258  ,nbineta,eta_l,eta_h,180,-180.,180.,0.)
259 
260 c----------------------------------------------------------------
261 c------------------------- Protons ------------------------------
262 c----------------------------------------------------------------
263  CALL hbook1(310,' Proton multiplicity distribution',
264  ,100,-0.5,float(izp+izt)+10.,0.)
265 
266  CALL hbook1(320,' Proton pseudo-rapidity distribution',
267  ,nbineta,eta_l,eta_h,0.)
268 
269  CALL hbook1(330,' Proton rapidity distribution',
270  ,nbiny,y_l,y_h,0.)
271 
272  CALL hbook1(340,' Proton Pt distribution',
273  ,100,0.,10.,0.)
274 
275  CALL hbook1(350,' Proton energy distribution',
276  ,nbine,e_l,e_h,0.)
277 
278  CALL hbook1(360,' Proton Cos(Theta) distribution',
279  ,40,-1.,1.,0.)
280 
281  CALL hbook1(370,' Proton Phi distribution',
282  ,180,0.,180.,0.)
283  CALL hbook2(380,' Phi - Eta correlation of protons',
284  ,nbineta,eta_l,eta_h,180,-180.,180.,0.)
285 
286 
287 c----------------------------------------------------------------
288 c----------------------- Gamma quanta ---------------------------
289 c----------------------------------------------------------------
290  CALL hbook1(410,' Gamma multiplicity distribution',
291  ,100,-0.5,float(m_neg),0.)
292 
293  CALL hbook1(420,' Gamma pseudo-rapidity distribution',
294  ,nbineta,eta_l,eta_h,0.)
295 
296  CALL hbook1(430,' Gamma rapidity distribution',
297  ,nbiny,y_l,y_h,0.)
298 
299  CALL hbook1(440,' Gamma Pt distribution',
300  ,100,0.,10.,0.)
301 
302  CALL hbook1(450,' Gamma energy distribution',
303  ,nbine,e_l,e_h,0.)
304 
305  CALL hbook1(460,' Gamma Cos(Theta) distribution',
306  ,40,-1.,1.,0.)
307 
308  CALL hbook1(470,' Gamma Phi distribution',
309  ,180,0.,180.,0.)
310 
311  CALL hbook2(480,' Phi - Eta correlation of Gammas',
312  ,nbineta,eta_l,eta_h,180,-180.,180.,0.)
313 c----------------------------------------------------------------
314 
315  pi=4.*atan(1.) ! 3.14159
316 
317  DO 2000 i_event=1,n_events
318 
319  WRITE(6,*)' Event # ',i_event,' ------------------'
320 C
321  CALL hijing(frame,bmin,bmax)
322 C
323  b=hint1(19)
324  Call hf1(2,b,1.)
325 
326  nu=n0+n01+n10+n11
327  Call hf1(4,float(nu),1.)
328 
329  Call hf1(6,float(np),1.)
330  Call hf1(8,float(nt),1.)
331 
332  n_ch=0
333  n_neg=0
334  n_pos=0
335  n_protons=0
336  n_gammas=0
337 
338  DO 3000 i=1,natt
339 
340  if(katt(i,2).eq. 0) go to 3000 ! reject non-interacting projectile
341  if(katt(i,2).eq. 1) go to 3000 ! reject elastic scattering
342  if(katt(i,2).eq.10) go to 3000 ! reject non-interacting target
343  if(katt(i,2).eq.11) go to 3000 ! reject elastic scattering
344 
345  id=katt(i,1)
346 
347  ich=luchge(katt(i,1))/3
348 
349  amass=ulmass(katt(i,1))
350 
351  e=patt(i,4)
352  pz=patt(i,3)
353  pt=sqrt(patt(i,1)**2+patt(i,2)**2)
354  p=sqrt(pz**2+pt**2)
355 
356  if(p-pz.ge.1.0e-4) then
357  eta= 0.5*alog((p+pz)/(p-pz))
358  elseif(p+pz.ge.1.0e-4) then
359  eta=-0.5*alog((p-pz)/(p+pz))
360  else
361  eta=15.
362  endif
363 
364  if(e-pz.ge.1.0e-4) then
365  y= 0.5*alog((e+pz)/(e-pz))
366  elseif(e+pz.ge.1.0e-4) then
367  y=-0.5*alog((e-pz)/(e+pz))
368  else
369  y=15.
370  endif
371 
372 
373  if(p.ge.1.0e-4) then
374  cos_theta=pz/p
375  else
376  cos_theta=1.
377  endif
378 
379  phi=ulangl(patt(i,1),patt(i,2))*180./pi
380 c==========================================================
381 
382  if(ich.ne.0) then
383  n_ch=n_ch+1
384 
385  Call hf1(20,eta,1.)
386  Call hf1(30, y,1.)
387  Call hf1(40, pt,1.)
388  Call hf1(50, e,1.)
389  Call hf1(60,cos_theta,1.)
390  Call hf1(70,phi,1.)
391  Call hfill(80,eta,phi,1.)
392  Call hf1(90,float(id),1.)
393  endif
394 
395  if(ich.lt.0) then
396  n_neg=n_neg+1
397 
398  Call hf1(120,eta,1.)
399  Call hf1(130, y,1.)
400  Call hf1(140, pt,1.)
401  Call hf1(150, e,1.)
402  Call hf1(160,cos_theta,1.)
403  Call hf1(170,phi,1.)
404  Call hfill(180,eta,phi,1.)
405  endif
406 
407  if(ich.gt.0) then
408  n_pos=n_pos+1
409 
410  Call hf1(220,eta,1.)
411  Call hf1(230, y,1.)
412  Call hf1(240, pt,1.)
413  Call hf1(250, e,1.)
414  Call hf1(260,cos_theta,1.)
415  Call hf1(270,phi,1.)
416  Call hfill(280,eta,phi,1.)
417  endif
418 
419  if(id.eq.2212) then
420  n_protons=n_protons+1
421 
422  Call hf1(320,eta,1.)
423  Call hf1(330, y,1.)
424  Call hf1(340, pt,1.)
425  Call hf1(350, e,1.)
426  Call hf1(360,cos_theta,1.)
427  Call hf1(370,phi,1.)
428  Call hfill(380,eta,phi,1.)
429  endif
430 
431  if(id.eq. 22) then
432  n_gammas=n_gammas+1
433 
434  Call hf1(420,eta,1.)
435  Call hf1(430, y,1.)
436  Call hf1(440, pt,1.)
437  Call hf1(450, e,1.)
438  Call hf1(460,cos_theta,1.)
439  Call hf1(470,phi,1.)
440  Call hfill(480,eta,phi,1.)
441  endif
442 
443 
444 3000 CONTINUE
445 
446 
447  CALL hf1( 10,float(n_ch),1.)
448  CALL hf1(110,float(n_neg),1.)
449  CALL hf1(210,float(n_pos),1.)
450  CALL hf1(310,float(n_protons),1.)
451  CALL hf1(410,float(n_gammas ),1.)
452 
453 2000 CONTINUE
454 
455 c================ Normalization ===========================
456 
457  c1=1./float(n_events) ! Multiplicity distr.
458  c2=0.
459 
460  Call hopera( 10,'+', 10, 10,c1,c2)
461  Call hopera(110,'+',110,110,c1,c2)
462  Call hopera(210,'+',210,210,c1,c2)
463  Call hopera(310,'+',310,310,c1,c2)
464  Call hopera(410,'+',410,410,c1,c2)
465 
466  c1=1./float(n_events)/((eta_h-eta_l)/nbineta)! Eta distr.
467  c2=0.
468 
469  Call hopera( 20,'+', 20, 20,c1,c2)
470  Call hopera(120,'+',120,120,c1,c2)
471  Call hopera(220,'+',220,220,c1,c2)
472  Call hopera(320,'+',320,320,c1,c2)
473  Call hopera(420,'+',420,420,c1,c2)
474 
475  c1=1./float(n_events)/((y_h-y_l)/nbiny) ! Y distr.
476  c2=0.
477 
478  Call hopera( 30,'+', 30, 30,c1,c2)
479  Call hopera(130,'+',130,130,c1,c2)
480  Call hopera(230,'+',230,230,c1,c2)
481  Call hopera(330,'+',330,330,c1,c2)
482  Call hopera(430,'+',430,430,c1,c2)
483 
484  c1=1./float(n_events)/0.1 ! Pt distr.
485  c2=0.
486 
487  Call hopera( 40,'+', 40, 40,c1,c2)
488  Call hopera(140,'+',140,140,c1,c2)
489  Call hopera(240,'+',240,240,c1,c2)
490  Call hopera(340,'+',340,340,c1,c2)
491  Call hopera(440,'+',440,440,c1,c2)
492 
493  c1=1./float(n_events)/((e_h-e_l)/nbine) ! E distr.
494  c2=0.
495 
496  Call hopera( 50,'+', 50, 50,c1,c2)
497  Call hopera(150,'+',150,150,c1,c2)
498  Call hopera(250,'+',250,250,c1,c2)
499  Call hopera(350,'+',350,350,c1,c2)
500  Call hopera(450,'+',450,450,c1,c2)
501 
502  c1=1./float(n_events)/0.05 ! Cos Theta distr.
503  c2=0.
504 
505  Call hopera( 60,'+', 60, 60,c1,c2)
506  Call hopera(160,'+',160,160,c1,c2)
507  Call hopera(260,'+',260,260,c1,c2)
508  Call hopera(360,'+',360,360,c1,c2)
509  Call hopera(460,'+',460,460,c1,c2)
510 
511  c1=1./float(n_events) ! Phi distr.
512  c2=0.
513 
514  Call hopera( 70,'+', 70, 70,c1,c2)
515  Call hopera(170,'+',170,170,c1,c2)
516  Call hopera(270,'+',270,270,c1,c2)
517  Call hopera(370,'+',370,370,c1,c2)
518  Call hopera(470,'+',470,470,c1,c2)
519 
520 c================ Writing results =========================
521 
522  CALL hrput(0,fname,'N')
523  END
524 
525  FUNCTION ran(NSEED)
526  ran=rlu(nseed)
527  RETURN
528  END