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