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