Geant4_10
dpm25hist.f
Go to the documentation of this file.
1 *-- Author :
2 C-------------------------------------------------------------------
3 C
4 C FILE DNDISTRM FORTRAN
5 C
6 C------------------------------------------------------------------
7 C DEBUG SUBCHK
8 C END DEBUG
9  SUBROUTINE distr(IOP,NHKKH1,PO,IGENER)
10  IMPLICIT DOUBLE PRECISION (a-h,o-z)
11 C***
12 C REDUCED VERSION FOR SCORING RESULTS FROM PRIMARY AA INTERACTION
13 C FROM /HKKEVT/
14 C kinematics in lab (target rest system)
15 C - MULTIPLICITIES
16 C - number of inc generations inside the nucleus
17 C - RAPIDITY DISTRIBUTION (without cuts)
18 C - pseudorapidity distribution (without cuts)
19 C***
20 C conventions for rapidity distributions:
21 C 1 Proton 11 NEGATIVES
22 C 2 Neutron 12 NBAR=9
23 C 3 PI+=13 13 LAMBDA=17
24 C 4 PI-=14 14 SIGMA==20,21,22
25 C 5 K+ =15 15 LAMBDABAR=18
26 C 6 K- =16 16 SIGMABAR=99,100,101
27 C 7 neutral kaons=12,19,24,25 17 THETA=97,98
28 C 8 pbar=2 18 THETABAR=102,103
29 C 9 charged hadrons
30 C 10 total hadrons
31 C-----------------------------------------------------------------------
32  dimension p4p4p4(4)
33 *KEEP,HKKEVT.
34 c INCLUDE (HKKEVT)
35  parameter(nmxhkk= 89998)
36 c PARAMETER (NMXHKK=25000)
37  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
38  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
39  +(4,nmxhkk)
40 C
41 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
42 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
43 C THE POSITIONS OF THE PROJECTILE NUCLEONS
44 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
45 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
46 C COMPLETELY CONSISTENT. THE TIMES IN THE
47 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
48 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
49 C
50 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
51 C
52 C NMXHKK: maximum numbers of entries (partons/particles) that can be
53 C stored in the commonblock.
54 C
55 C NHKK: the actual number of entries stored in current event. These are
56 C found in the first NHKK positions of the respective arrays below.
57 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
58 C entry.
59 C
60 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
61 C = 0 : null entry.
62 C = 1 : an existing entry, which has not decayed or fragmented.
63 C This is the main class of entries which represents the
64 C "final state" given by the generator.
65 C = 2 : an entry which has decayed or fragmented and therefore
66 C is not appearing in the final state, but is retained for
67 C event history information.
68 C = 3 : a documentation line, defined separately from the event
69 C history. (incoming reacting
70 C particles, etc.)
71 C = 4 - 10 : undefined, but reserved for future standards.
72 C = 11 - 20 : at the disposal of each model builder for constructs
73 C specific to his program, but equivalent to a null line in the
74 C context of any other program. One example is the cone defining
75 C vector of HERWIG, another cluster or event axes of the JETSET
76 C analysis routines.
77 C = 21 - : at the disposal of users, in particular for event tracking
78 C in the detector.
79 C
80 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
81 C standard.
82 C
83 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
84 C The value is 0 for initial entries.
85 C
86 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
87 C one mother exist, in which case the value 0 is used. In cluster
88 C fragmentation models, the two mothers would correspond to the q
89 C and qbar which join to form a cluster. In string fragmentation,
90 C the two mothers of a particle produced in the fragmentation would
91 C be the two endpoints of the string (with the range in between
92 C implied).
93 C
94 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
95 C entry has not decayed, this is 0.
96 C
97 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
98 C entry has not decayed, this is 0. It is assumed that the daughters
99 C of a particle (or cluster or string) are stored sequentially, so
100 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
101 C daughters. Even in cases where only one daughter is defined (e.g.
102 C K0 -> K0S) both values should be defined, to make for a uniform
103 C approach in terms of loop constructions.
104 C
105 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
106 C
107 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
108 C
109 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
110 C
111 C PHKK(4,IHKK) : energy, in GeV.
112 C
113 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
114 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
115 C
116 C VHKK(1,IHKK) : production vertex x position, in mm.
117 C
118 C VHKK(2,IHKK) : production vertex y position, in mm.
119 C
120 C VHKK(3,IHKK) : production vertex z position, in mm.
121 C
122 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
123 C********************************************************************
124 *KEEP,NSHMAK.
125  COMMON /nshmak/ nnshma,npshma,ntshma,nshmac
126 *KEEP,DPAR.
127 C /DPAR/ CONTAINS PARTICLE PROPERTIES
128 C ANAME = LITERAL NAME OF THE PARTICLE
129 C AAM = PARTICLE MASS IN GEV
130 C GA = DECAY WIDTH
131 C TAU = LIFE TIME OF INSTABLE PARTICLES
132 C IICH = ELECTRIC CHARGE OF THE PARTICLE
133 C IIBAR = BARYON NUMBER
134 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
135 C
136  CHARACTER*8 aname
137  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
138  +iibar(210),k1(210),k2(210)
139 C------------------
140 *KEEP,NUCC.
141  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
142 *KEEP,DPRIN.
143  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
144 *KEND.
145 C modified DPMJET
146  COMMON /bufueh/ annvv,annss,annsv,annvs,anncc,
147  * anndv,annvd,annds,annsd,
148  * annhh,annzz,
149  * ptvv,ptss,ptsv,ptvs,ptcc,ptdv,ptvd,ptds,ptsd,
150  * pthh,ptzz,
151  * eevv,eess,eesv,eevs,eecc,eedv,eevd,eeds,eesd,
152  * eehh,eezz
153  * ,anndi,ptdi,eedi
154  * ,annzd,anndz,ptzd,ptdz,eezd,eedz
155 C---------------------
156 C---------------------
157  dimension xyl(50,20),yyl(50,20),yylps(50,20),indx(28)
158  dimension disgen(50),xgen(50)
159  parameter(numtyp=160)
160  dimension avmult(numtyp),ave(numtyp),ake(numtyp),avept(numtyp)
161  DATA disgen /50*0.0/
162  DATA indx/ 1, 8,-1,-1,-1,-1,-1, 2,12,-1,-1,7,3,4,5,6,13,15,7,14,
163  * 14,14,19, 7, 7,-1,-1,-1/
164 C----------------------------------------------------------------------
165  oneone=1
166  zero=0
167  twotwo=2
168  hundth=1.d-2
169 C CALL DISPT(IOP,NHKKH1,PO)
170 C CALL DISEVA(IOP,NHKKH1,PO,IGENER)
171  go to(1,2,3),iop
172 * initialization call
173 1 CONTINUE
174  anchsq=0.
175  delrap=1.
176  IF(ip.EQ.1)delrap=0.1
177  nhkkh2=nhkkh1
178  IF (nhkkh1.EQ.0)nhkkh2=1
179  eeo=sqrt(po**2+aam(nhkkh2)**2)
180  WRITE(6, 1001)eeo,po,nhkkh2,aam(nhkkh2)
181  1001 FORMAT (' EEO',f10.2,f10.2,i10,f10.2)
182  dy=0.4
183 C
184  DO 11 i=1,20
185  DO 13 j=1,50
186  xyl(j,i)=-2.0 + (j-1)*dy
187  yyl(j,i)=1.e-18
188  yylps(j,i)=1.e-18
189 13 CONTINUE
190 11 CONTINUE
191  DO 61 i=1,numtyp
192  ave(i) =1.e-18
193  avept(i) =1.e-18
194  avmult(i) =1.e-18
195  61 CONTINUE
196  RETURN
197 C-------------------------------------------------------------------
198  2 CONTINUE
199 C SEWEW secondary interactions
200  CALL sewew(1,nhkkh1)
201 * scoring of results from individual events
202 C NUMBER OF GENERATIONS INSIDE THE NUCLEUS
203  IF(igener.LT.1.OR.igener.GT.50) igener=50
204  disgen(igener)=disgen(igener) + 1.0
205 C
206 C WRITE (18,1711)
207 C1711 FORMAT (' event px,py,pz,e,id ')
208  nhad=0
209  avmulc=0.
210 C HBOOK HISTOGRAMS
211 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
212 C WRITE(6,'(A)')' before plomb'
213  IF(ihbook.EQ.1)CALL plomb(2,p4p4p4,ccchrg,xfxfxf,1,ijproj)
214 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
215 C Here we include the decayed resonances into the
216 C multiplicity Table
217  ipriev=0
218  DO 521 i=nhkkh1,nhkk
219  IF (isthkk(i).EQ.2)THEN
220 C WRITE(18,1712)PHKK(1,I),PHKK(2,I),PHKK(3,I),Phkk(4,I),IDHKK(I)
221 C NHAD=NHAD+1
222  nrhkk=mcihad(idhkk(i))
223  IF (nrhkk.LE.0.OR.nrhkk.GT.210)THEN
224  WRITE(6,1389)nrhkk,i,idhkk(i),nhkkh1,nhkk
225  nrhkk=1
226  ENDIF
227  nre=nrhkk
228  nrem=nre
229  ichhkk=iich(nrhkk)
230  IF(nre.GT.30)THEN
231  IF(nre.GT.160)go to 521
232  ave(nrem)=ave(nrem) + phkk(4,i)
233  avept(nrem)=avept(nrem) + pt
234  avmult(nrem)=avmult(nrem) + 1.
235  ENDIF
236  ENDIF
237  521 CONTINUE
238  DO 21 i=nhkkh1,nhkk
239  IF (isthkk(i).EQ.1)THEN
240 C WRITE(18,1712)PHKK(1,I),PHKK(2,I),PHKK(3,I),Phkk(4,I),IDHKK(I)
241  1712 FORMAT (4e14.5,i8)
242  nhad=nhad+1
243  nrhkk=mcihad(idhkk(i))
244  IF (nrhkk.LE.0.OR.nrhkk.GT.210)THEN
245  WRITE(6,1389)nrhkk,i,idhkk(i),nhkkh1,nhkk
246  1389 FORMAT (' DISTR: NRHKK ERROR ',5i10)
247  nrhkk=1
248  ENDIF
249  nre=nrhkk
250  nrem=nre
251  ichhkk=iich(nrhkk)
252 * rapidity
253  ptt=phkk(1,i)**2+phkk(2,i)**2+0.000001
254  pt=sqrt(ptt)
255 C IF(PT.LT.0.5)GO TO 21
256  amt=sqrt(ptt+phkk(5,i)**2)
257 C*** AMT=SQRT(PTT+AMPI**2)
258 C ETOT=SQRT(AMT**2 + PHKK(3,I)**2)
259 C YL=LOG((ETOT + PHKK(3,I))/AMT+1.E-18)
260  yl=log((abs(phkk(3,i) + phkk(4,i)))/amt+1.e-18)
261 C IF(YL.LT.0.75.OR.YL.GT.3.) GO TO 21
262  IF (nre.GT.25) nre=28
263  IF (nre.LT. 1) nre=28
264  IF(nrem.GT.numtyp) nrem=28
265  IF(nrem.LT.1) nrem=28
266  ni=indx(nre)
267  IF (nrhkk.LE.101.AND.nrhkk.GE.99) ni=16
268  IF (nrhkk.EQ.97.OR.nrhkk.EQ.98) ni=17
269  IF (nrhkk.EQ.102.OR.nrhkk.EQ.103) ni=18
270  ave(nrem)=ave(nrem) + phkk(4,i)
271  avept(nrem)=avept(nrem) + pt
272  avmult(nrem)=avmult(nrem) + 1.
273 C TOTAL=30
274  ave(30)=ave(30) + phkk(4,i)
275  avept(30)=avept(30) + pt
276 C CHARGED=27
277  IF (ichhkk.NE.0) THEN
278  ave(27)=ave(27) + phkk(4,i)
279  avept(27)=avept(27) + pt
280  avmulc=avmulc + 1.
281  avmult(27)=avmult(27) + 1.
282  ENDIF
283  iyl=(yl+2.0)/dy + 1
284  IF (iyl.LT.1) iyl=1
285  IF (iyl.GT.50) iyl=50
286  IF (ichhkk.NE.0) THEN
287  yyl(iyl,9)=yyl(iyl,9)+1.
288  ENDIF
289  IF (ichhkk.LT.0) yyl(iyl,11)=yyl(iyl,11)+1.
290  IF(ni.GT.0) yyl(iyl,ni)=yyl(iyl,ni)+1.
291  yyl(iyl,10)=yyl(iyl,10)+1.
292 * pseudorapidity
293  pt=sqrt(ptt)
294  ptot=sqrt(ptt+phkk(3,i)**2)
295  ylps=log((ptot+phkk(3,i))/pt)
296  iyl=(ylps+2.0)/dy + 1
297  IF (iyl.LT.1) iyl=1
298  IF (iyl.GT.50) iyl=50
299  IF (ichhkk.NE.0) yylps(iyl,9)=yylps(iyl,9)+1.
300  IF (ichhkk.LT.0) yylps(iyl,11)=yylps(iyl,11)+1.
301  IF(ni.GT.0) yylps(iyl,ni)=yylps(iyl,ni)+1.
302  yylps(iyl,10)=yylps(iyl,10)+1.
303 C HBOOK HISTOGRAMS
304 C HBOOK HISTOGRAMS
305 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
306  p4p4p4(1) = phkk(1,i)
307  p4p4p4(2) = phkk(2,i)
308  p4p4p4(3) = phkk(3,i)
309  p4p4p4(4) = phkk(4,i)
310  itif = nrhkk
311  xfl=phkk(4,i)/eeo
312  xfxfxf=xfl
313  ccchrg = ichhkk
314  IF(ihbook.EQ.1)CALL plomb(3,p4p4p4,
315  * ccchrg,xfxfxf,itif,ijproj)
316 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
317  ENDIF
318 21 CONTINUE
319 * total multiplicity
320  avmult(30)=avmult(30) + nhad
321  anchsq=anchsq+avmulc**2
322  anhad=nhad
323 C HBOOK HISTOGRAMS
324 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
325  IF(ihbook.EQ.1)CALL plomb(4,p4p4p4,
326  * ccchrg,xfxfxf,1,ijproj)
327 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
328  RETURN
329 C
330 C--------------------------------------------------------------------
331 C
332  3 CONTINUE
333 * output of final results
334 C WRITE(6,'(1H1,50(1H*))')
335 C WRITE(6,'(/10X,A/)')'OUTP.FR.DISTR 0.75.lt.y.lt.3., pt.gt.o.5'
336 C WRITE(6,'(50(1H*))')
337 C NORMALIZE AND PRINT COUNTERS
338  IF (annvv.GT.1.)THEN
339  ptvv=ptvv/annvv
340  eevv=eevv/annvv
341  ENDIF
342  IF (annss.GT.1.)THEN
343  ptss=ptss/annss
344  eess=eess/annss
345  ENDIF
346  IF (annsv.GT.1.)THEN
347  ptsv=ptsv/annsv
348  eesv=eesv/annsv
349  ENDIF
350  IF (annvs.GT.1.)THEN
351  ptvs=ptvs/annvs
352  eevs=eevs/annvs
353  ENDIF
354  IF (anncc.GE.1.)THEN
355  ptcc=ptcc/anncc
356  eecc=eecc/anncc
357  ENDIF
358  IF (annvd.GT.1.)THEN
359  ptvd=ptvd/annvd
360  eevd=eevd/annvd
361  ENDIF
362  IF (anndv.GT.1.)THEN
363  ptdv=ptdv/anndv
364  eedv=eedv/anndv
365  ENDIF
366  IF (annsd.GT.1.)THEN
367  ptsd=ptsd/annsd
368  eesd=eesd/annsd
369  ENDIF
370  IF (annds.GT.1.)THEN
371  ptds=ptds/annds
372  eeds=eeds/annds
373  ENDIF
374  IF (annhh.GE.1.)THEN
375  pthh=pthh/annhh
376  eehh=eehh/annhh
377  ENDIF
378  IF (annzz.GE.1.)THEN
379  ptzz=ptzz/annzz
380  eezz=eezz/annzz
381  ENDIF
382  IF (nhkkh1.GT.1.)THEN
383  annvv=annvv/nhkkh1
384  annss=annss/nhkkh1
385  annsv=annsv/nhkkh1
386  annvs=annvs/nhkkh1
387  anncc=anncc/nhkkh1
388  annhh=annhh/nhkkh1
389  annzz=annzz/nhkkh1
390  anndv=anndv/nhkkh1
391  annvd=annvd/nhkkh1
392  annds=annds/nhkkh1
393  annsd=annsd/nhkkh1
394  anchsq=anchsq/nhkkh1
395 C WRITE (6,7431)ANNVV,PTVV,EEVV,
396 C * ANNSS,PTSS,EESS,
397 C * ANNSV,PTSV,EESV,
398 C * ANNVS,PTVS,EEVS,
399 C * ANNDV,PTDV,EEDV,
400 C * ANNVD,PTVD,EEVD,
401 C * ANNDS,PTDS,EEDS,
402 C * ANNSD,PTSD,EESD,
403 C * ANNHH,PTHH,EEHH,
404 C * ANNZZ,PTZZ,EEZZ,
405 C * ANNCC,PTCC,EECC
406 C7431 FORMAT (' VV CHAINS NN,PT ECM: ',3F12.4/
407 C * ' SS CHAINS NN,PT ECM: ',3F12.4/
408 C * ' SV CHAINS NN,PT ECM: ',3F12.4/
409 C * ' VS CHAINS NN,PT ECM: ',3F12.4/
410 C * ' DV CHAINS NN,PT ECM: ',3F12.4/
411 C * ' VD CHAINS NN,PT ECM: ',3F12.4/
412 C * ' DS CHAINS NN,PT ECM: ',3F12.4/
413 C * ' SD CHAINS NN,PT ECM: ',3F12.4/
414 C * ' HH CHAINS NN,PT ECM: ',3F12.4/
415 C * ' ZZ CHAINS NN,PT ECM: ',3F12.4/
416 C * ' CC CHAINS NN,PT ECM: ',3F12.4)
417  ENDIF
418 C
419 C DISTRIBUTION FOR
420 C NUMBER OF GENERATIONS INSIDE THE NUCLEUS
421 C (NORMALIZED TO 1)
422  avgnor=0.0
423  avgen=0.0
424  DO 801 ige=1,50
425  avgnor=avgnor + disgen(ige)
426  avgen=avgen + ige*disgen(ige)
427  801 CONTINUE
428  DO 802 ige=1,50
429  xgen(ige)=float(ige) - 0.5
430  disgen(ige)=disgen(ige)/avgnor
431  802 CONTINUE
432  avgen=avgen/avgnor
433 C WRITE(*,'(A,1PE10.2//A)')
434 C & ' AVERAGE NUMBER OF GENERATIONS INSIDE THE NUCLEUS =',
435 C & AVGEN,
436 C & ' CORRESPONDING DISTRIBUTION (NORMALIZED TO 1.0)'
437 C CALL PLOT(XGEN,DISGEN,50,1,50,ZERO,ONEONE,ZERO,HUNDTH)
438 C
439  DO 310 i=1,numtyp
440  avmult(i)=avmult(i)/nhkkh1
441  ave(i)=ave(i)/nhkkh1
442  avept(i)=avept(i)/(nhkkh1*avmult(i))
443 310 CONTINUE
444  ffff2=sqrt(anchsq-avmult(27)**2)/avmult(27)
445 C WRITE(6,7772)FFFF2,ANCHSQ,AVMULT(27)
446 C7772 FORMAT(' FFFF2,ANCHSQ,AVMULT(27):',3E15.5)
447  WRITE(6, 64)
448  64 FORMAT(' PARTICLE REF,CHAR,IBAR, MASS AVERAGE',
449  *' ENERGY, MULTIPLICITY, INELASTICITY')
450  DO 62 i=1,numtyp
451  ake(i)=ave(i)/eeo
452  WRITE(6, 63) aname(i),i,iich(i),iibar(i),aam(i),
453  * ave(i),avmult(i),ake(i),avept(i)
454  63 FORMAT (' ',a8,3i5,f10.3,4f15.6)
455  62 CONTINUE
456 C
457  DO 33 i=1,20
458  DO 35 j=1,50
459  yyl(j,i) =yyl(j,i) /(nhkkh1*dy)
460  yylps(j,i)=yylps(j,i)/(nhkkh1*dy)
461 35 CONTINUE
462 33 CONTINUE
463  WRITE(6,'(1H1,11(A/))')
464  & ' Conventions for rapidity distributions: ',
465  & ' 1 Proton 11 NEGATIVES ',
466  & ' 2 Neutron 12 NBAR=9 ',
467  & ' 3 PI+=13 13 LAMBDA=17 ',
468  & ' 4 PI-=14 14 SIGMA==20,21,22 ',
469  & ' 5 K+ =15 15 LAMBDABAR=18 ',
470  & ' 6 K- =16 16 SIGMABAR=99,100,101',
471  & ' 7 neutral kaons=12,19,24,25 17 THETA=97,98 ',
472  & ' 8 pbar=2 18 THETABAR=102,103 ',
473  & ' 9 charged hadrons ',
474  & ' 10 total hadrons '
475  WRITE(6, 66)
476  66 FORMAT(' RAPIDITY DISTRIBUTION')
477  WRITE(6,302)
478  302 FORMAT (' (first number gives the lower bin limit)')
479  DO 36 j=1,50
480  WRITE(6, 37) xyl(j,1),(yyl(j,i),i=1,10)
481 37 FORMAT (f10.2,10e11.3)
482 36 CONTINUE
483  WRITE(6, 66)
484  WRITE(6,302)
485  DO 306 j=1,50
486  WRITE(6, 37) xyl(j,1),(yyl(j,i),i=11,20)
487 306 CONTINUE
488 C
489  WRITE(6, 66 )
490  CALL plot(xyl,yyl,1000,20,50,-twotwo,dy,zero,delrap)
491 * rescaled rapidity distribution
492  CALL plot(xyl,yyl,1000,20,50,-twotwo,dy,zero,5.*delrap)
493 * logarithmic rapidity distribution
494  DO 38 i=1,20
495  DO 38 j=1,50
496  yyl(j,i)=log10(abs(yyl(j,i))+1.e-8)
497  38 CONTINUE
498  WRITE(6, 66)
499  CALL plot(xyl,yyl,1000,20,50,-twotwo,dy,-twotwo,5.*hundth)
500 C
501  WRITE(6,301)
502  301 FORMAT ('1 PSEUDORAPIDITY DISTRIBUTION')
503  WRITE(6,302)
504  DO 303 j=1,50
505  WRITE(6,37) xyl(j,1),(yylps(j,i),i=1,10)
506  303 CONTINUE
507  WRITE(6,301)
508  WRITE(6,302)
509  DO 304 j=1,50
510  WRITE(6,37) xyl(j,1),(yylps(j,i),i=11,20)
511  304 CONTINUE
512  WRITE(6,301)
513  CALL plot(xyl,yylps,1000,20,50,-twotwo,dy,zero,delrap)
514  CALL plot(xyl,yylps,1000,20,50,-twotwo,dy,zero,5.*delrap)
515  RETURN
516  END
517 *CMZ : 1.00/02 08/01/92 17.24.06 by H.-J. Moehring & J. Ranft
518 *CMZ : 1.00/01 25/11/91 15.30.00 by H.-J. M¶hring
519 *-- Author :
520 C
521 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
522 C
523 C
524 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
525  SUBROUTINE distrc(IOP,NHKKH1,PO,IGENER)
526  IMPLICIT DOUBLE PRECISION (a-h,o-z)
527 C***
528 C REDUCED VERSION FOR SCORING RESULTS FROM PRIMARY AA INTERACTION
529 C FROM /HKKEVT/
530 C kinematics in cms
531 C - MULTIPLICITIES
532 C - number of inc generations inside the nucleus
533 C - RAPIDITY DISTRIBUTION (without cuts)
534 C - pseudorapidity distribution (without cuts)
535 C***
536 C conventions for rapidity distributions:
537 C 1 Proton 11 NEGATIVES
538 C 2 Neutron 12 NBAR=9
539 C 3 PI+=13 13 LAMBDA=17
540 C 4 PI-=14 14 SIGMA==20,21,22
541 C 5 K+ =15 15 LAMBDABAR=18
542 C 6 K- =16 16 SIGMABAR=99,100,101
543 C 7 neutral kaons=12,19,24,25 17 THETA=97,98
544 C 8 pbar=2 18 THETABAR=102,103
545 C 9 charged hadrons
546 C 10 total hadrons
547 C-----------------------------------------------------------------------
548 *KEEP,HKKEVT.
549 c INCLUDE (HKKEVT)
550 C PARAMETER (NMXHKK=89998)
551  parameter(nmxhkk= 89998)
552 c PARAMETER (NMXHKK=25000)
553  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
554  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
555  +(4,nmxhkk)
556 C
557 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
558 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
559 C THE POSITIONS OF THE PROJECTILE NUCLEONS
560 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
561 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
562 C COMPLETELY CONSISTENT. THE TIMES IN THE
563 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
564 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
565 C
566 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
567 C
568 C NMXHKK: maximum numbers of entries (partons/particles) that can be
569 C stored in the commonblock.
570 C
571 C NHKK: the actual number of entries stored in current event. These are
572 C found in the first NHKK positions of the respective arrays below.
573 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
574 C entry.
575 C
576 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
577 C = 0 : null entry.
578 C = 1 : an existing entry, which has not decayed or fragmented.
579 C This is the main class of entries which represents the
580 C "final state" given by the generator.
581 C = 2 : an entry which has decayed or fragmented and therefore
582 C is not appearing in the final state, but is retained for
583 C event history information.
584 C = 3 : a documentation line, defined separately from the event
585 C history. (incoming reacting
586 C particles, etc.)
587 C = 4 - 10 : undefined, but reserved for future standards.
588 C = 11 - 20 : at the disposal of each model builder for constructs
589 C specific to his program, but equivalent to a null line in the
590 C context of any other program. One example is the cone defining
591 C vector of HERWIG, another cluster or event axes of the JETSET
592 C analysis routines.
593 C = 21 - : at the disposal of users, in particular for event tracking
594 C in the detector.
595 C
596 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
597 C standard.
598 C
599 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
600 C The value is 0 for initial entries.
601 C
602 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
603 C one mother exist, in which case the value 0 is used. In cluster
604 C fragmentation models, the two mothers would correspond to the q
605 C and qbar which join to form a cluster. In string fragmentation,
606 C the two mothers of a particle produced in the fragmentation would
607 C be the two endpoints of the string (with the range in between
608 C implied).
609 C
610 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
611 C entry has not decayed, this is 0.
612 C
613 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
614 C entry has not decayed, this is 0. It is assumed that the daughters
615 C of a particle (or cluster or string) are stored sequentially, so
616 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
617 C daughters. Even in cases where only one daughter is defined (e.g.
618 C K0 -> K0S) both values should be defined, to make for a uniform
619 C approach in terms of loop constructions.
620 C
621 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
622 C
623 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
624 C
625 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
626 C
627 C PHKK(4,IHKK) : energy, in GeV.
628 C
629 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
630 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
631 C
632 C VHKK(1,IHKK) : production vertex x position, in mm.
633 C
634 C VHKK(2,IHKK) : production vertex y position, in mm.
635 C
636 C VHKK(3,IHKK) : production vertex z position, in mm.
637 C
638 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
639 C********************************************************************
640 *KEEP,NSHMAK.
641  COMMON /nshmak/ nnshma,npshma,ntshma,nshmac
642 *KEEP,DPAR.
643 C /DPAR/ CONTAINS PARTICLE PROPERTIES
644 C ANAME = LITERAL NAME OF THE PARTICLE
645 C AAM = PARTICLE MASS IN GEV
646 C GA = DECAY WIDTH
647 C TAU = LIFE TIME OF INSTABLE PARTICLES
648 C IICH = ELECTRIC CHARGE OF THE PARTICLE
649 C IIBAR = BARYON NUMBER
650 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
651 C
652  CHARACTER*8 aname
653  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
654  +iibar(210),k1(210),k2(210)
655 C------------------
656 *KEEP,NUCC.
657  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
658 *KEEP,DPRIN.
659  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
660 *KEND.
661 C modified DPMJET
662  COMMON /bufueh/ annvv,annss,annsv,annvs,anncc,
663  * anndv,annvd,annds,annsd,
664  * annhh,annzz,
665  * ptvv,ptss,ptsv,ptvs,ptcc,ptdv,ptvd,ptds,ptsd,
666  * pthh,ptzz,
667  * eevv,eess,eesv,eevs,eecc,eedv,eevd,eeds,eesd,
668  * eehh,eezz
669  * ,anndi,ptdi,eedi
670  * ,annzd,anndz,ptzd,ptdz,eezd,eedz
671 C---------------------
672  CHARACTER*80 title
673  CHARACTER*8 projty,targty
674 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
675 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
676  COMMON /user1/title,projty,targty
677  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
678  dimension xyl(50,20),yyl(50,20),yylps(50,20),indx(28)
679  dimension disgen(50),xgen(50)
680  parameter(numtyp=160)
681  dimension avmult(numtyp),ave(numtyp),ake(numtyp),avept(numtyp)
682  DATA disgen /50*0.0/
683  DATA indx/ 1, 8,-1,-1,-1,-1,-1, 2,12,-1,-1,7,3,4,5,6,13,15,7,14,
684  * 14,14,19, 7, 7,-1,-1,-1/
685  DATA niprie /0/
686 C----------------------------------------------------------------------
687 C CALL DISPT(IOP,NHKKH1,PO)
688  go to(1,2,3),iop
689 * initialization call
690 1 CONTINUE
691 C initialize ALICE output
692 C CALL SHINIT
693 C
694  anchsq=0.
695  delrap=1.
696  IF(ip.EQ.1)delrap=0.1
697  nhkkh2=nhkkh1
698  IF (nhkkh1.EQ.0)nhkkh2=1
699  eeo=sqrt(po**2+aam(nhkkh2)**2)
700  WRITE(6, 1001)eeo,po,nhkkh2,aam(nhkkh2)
701  1001 FORMAT (' EEO',f10.2,f10.2,i10,f10.2)
702  dy=0.4
703 C
704  DO 11 i=1,20
705  DO 13 j=1,50
706  xyl(j,i)=-10.0 + (j-1)*dy
707  yyl(j,i)=1.e-18
708  yylps(j,i)=1.e-18
709 13 CONTINUE
710 11 CONTINUE
711  DO 61 i=1,numtyp
712  ave(i) =1.e-18
713  avept(i) =1.e-18
714  avmult(i) =1.e-18
715  61 CONTINUE
716  RETURN
717 C-------------------------------------------------------------------
718  2 CONTINUE
719 C SEWEW secondary interactions
720  CALL sewew(1,nhkkh1)
721  ievt=ievt+1
722 C output for ALICE
723 C CALL SHEVUT(IEVT)
724 C
725 * scoring of results from individual events
726 C NUMBER OF GENERATIONS INSIDE THE NUCLEUS
727  IF(igener.LT.1.OR.igener.GT.50) igener=50
728  disgen(igener)=disgen(igener) + 1.0
729 C
730 C WRITE (18,1711)
731 C1711 FORMAT (' event px,py,pz,e,id ')
732  nhad=0
733  avmulc=0.
734 C Here we include the decayed resonances into the
735 C multiplicity Table
736  ipriev=0
737  DO 521 i=nhkkh1,nhkk
738  IF (isthkk(i).EQ.2)THEN
739 C WRITE(18,1712)PHKK(1,I),PHKK(2,I),PHKK(3,I),Phkk(4,I),IDHKK(I)
740 C NHAD=NHAD+1
741  nrhkk=mcihad(idhkk(i))
742  IF (nrhkk.LE.0.OR.nrhkk.GT.210)THEN
743  WRITE(6,1389)nrhkk,i,idhkk(i),nhkkh1,nhkk
744  nrhkk=1
745  ENDIF
746  nre=nrhkk
747  nrem=nre
748  ichhkk=iich(nrhkk)
749  IF(nre.GT.30)THEN
750  IF(nre.GT.160)go to 521
751  avept(nrem)=avept(nrem) + pt
752  avmult(nrem)=avmult(nrem) + 1.
753  ENDIF
754  ENDIF
755  521 CONTINUE
756  DO 21 i=nhkkh1,nhkk
757  IF (isthkk(i).EQ.1)THEN
758 C WRITE(18,1712)PHKK(1,I),PHKK(2,I),PHKK(3,I),Phkk(4,I),IDHKK(I)
759  1712 FORMAT (4e14.5,i8)
760  IF(phkk(4,i).GT.cmener/2.d0)THEN
761  ipriev=1
762  niprie=niprie+1
763  ENDIF
764  nhad=nhad+1
765  nrhkk=mcihad(idhkk(i))
766  IF (nrhkk.LE.0.OR.nrhkk.GT.210)THEN
767  WRITE(6,1389)nrhkk,i,idhkk(i),nhkkh1,nhkk
768  1389 FORMAT (' DISTR: NRHKK ERROR ',5i10)
769  nrhkk=1
770  ENDIF
771  nre=nrhkk
772  nrem=nre
773  ichhkk=iich(nrhkk)
774 * rapidity
775  ptt=phkk(1,i)**2+phkk(2,i)**2+0.000001
776  pt=sqrt(ptt)
777 C IF(PT.LT.0.5)GO TO 21
778  amt=sqrt(ptt+phkk(5,i)**2)
779 C*** AMT=SQRT(PTT+AMPI**2)
780 C ETOT=SQRT(AMT**2 + PHKK(3,I)**2)
781 C YL=LOG((ETOT + PHKK(3,I))/AMT+1.E-18)
782  yl=log((abs(phkk(3,i) + phkk(4,i)))/amt+1.e-18)
783 C IF(YL.LT.0.75.OR.YL.GT.3.) GO TO 21
784  IF (nre.GT.25) nre=28
785  IF (nre.LT. 1) nre=28
786  IF(nrem.GT.numtyp) nrem=28
787  IF(nrem.LT.1) nrem=28
788  ni=indx(nre)
789  IF (nrhkk.LE.101.AND.nrhkk.GE.99) ni=16
790  IF (nrhkk.EQ.97.OR.nrhkk.EQ.98) ni=17
791  IF (nrhkk.EQ.102.OR.nrhkk.EQ.103) ni=18
792  ave(nrem)=ave(nrem) + phkk(4,i)
793  avept(nrem)=avept(nrem) + pt
794  avmult(nrem)=avmult(nrem) + 1.
795 C TOTAL=30
796  ave(30)=ave(30) + phkk(4,i)
797  avept(30)=avept(30) + pt
798 C CHARGED=27
799  IF (ichhkk.NE.0) THEN
800  ave(27)=ave(27) + phkk(4,i)
801  avept(27)=avept(27) + pt
802  avmulc=avmulc + 1.
803  avmult(27)=avmult(27) + 1.
804  ENDIF
805  iyl=(yl+10.0)/dy + 1
806  IF (iyl.LT.1) iyl=1
807  IF (iyl.GT.50) iyl=50
808  IF (ichhkk.NE.0) THEN
809  yyl(iyl,9)=yyl(iyl,9)+1.
810  ENDIF
811  IF (ichhkk.LT.0) yyl(iyl,11)=yyl(iyl,11)+1.
812  IF(ni.GT.0) yyl(iyl,ni)=yyl(iyl,ni)+1.
813  yyl(iyl,10)=yyl(iyl,10)+1.
814 * pseudorapidity
815  pt=sqrt(ptt)
816  ptot=sqrt(ptt+phkk(3,i)**2)
817  ylps=log((ptot+phkk(3,i))/pt)
818  iyl=(ylps+10.0)/dy + 1
819  IF (iyl.LT.1) iyl=1
820  IF (iyl.GT.50) iyl=50
821  IF (ichhkk.NE.0) yylps(iyl,9)=yylps(iyl,9)+1.
822  IF (ichhkk.LT.0) yylps(iyl,11)=yylps(iyl,11)+1.
823  IF(ni.GT.0) yylps(iyl,ni)=yylps(iyl,ni)+1.
824  yylps(iyl,10)=yylps(iyl,10)+1.
825  ENDIF
826 21 CONTINUE
827  IF(ipriev.EQ.-1)THEN
828  IF(niprie.LT.20)THEN
829  WRITE(6,*)' IPRIEV = 1 '
830  DO 120 ihkk=1,nhkk
831  WRITE(6,1050) ihkk,isthkk(ihkk),idhkk(ihkk),jmohkk(1,ihkk),
832  + jmohkk(2,ihkk), jdahkk(1,ihkk),jdahkk(2,ihkk),
833  + (phkk(khkk,ihkk),khkk=1,5)
834 120 CONTINUE
835  1050 FORMAT (i6,i4,5i6,5e16.8)
836  ENDIF
837  ENDIF
838 * total multiplicity
839  avmult(30)=avmult(30) + nhad
840  anchsq=anchsq+avmulc**2
841  anhad=nhad
842  RETURN
843 C
844 C--------------------------------------------------------------------
845 C
846  3 CONTINUE
847 C output for ALICE
848 C CALL SHEND
849 C
850 * output of final results
851 C WRITE(6,'(1H1,50(1H*))')
852 C WRITE(6,'(/10X,A/)')'OUTP.FR.DISTR 0.75.lt.y.lt.3., pt.gt.o.5'
853 C WRITE(6,'(50(1H*))')
854 C NORMALIZE AND PRINT COUNTERS
855  IF (annvv.GT.1.)THEN
856  ptvv=ptvv/annvv
857  eevv=eevv/annvv
858  ENDIF
859  IF (annss.GT.1.)THEN
860  ptss=ptss/annss
861  eess=eess/annss
862  ENDIF
863  IF (annsv.GT.1.)THEN
864  ptsv=ptsv/annsv
865  eesv=eesv/annsv
866  ENDIF
867  IF (annvs.GT.1.)THEN
868  ptvs=ptvs/annvs
869  eevs=eevs/annvs
870  ENDIF
871  IF (anncc.GE.1.)THEN
872  ptcc=ptcc/anncc
873  eecc=eecc/anncc
874  ENDIF
875  IF (nhkkh1.GT.1.)THEN
876  annvv=annvv/nhkkh1
877  annss=annss/nhkkh1
878  annsv=annsv/nhkkh1
879  annvs=annvs/nhkkh1
880  anncc=anncc/nhkkh1
881  anchsq=anchsq/nhkkh1
882  WRITE (6,7431)annvv,ptvv,eevv,
883  * annss,ptss,eess,
884  * annsv,ptsv,eesv,
885  * annvs,ptvs,eevs,
886  * anncc,ptcc,eecc
887  7431 FORMAT (' VV CHAINS NN,PT ECM: ',3f12.4/
888  * ' SS CHAINS NN,PT ECM: ',3f12.4/
889  * ' SV CHAINS NN,PT ECM: ',3f12.4/
890  * ' VS CHAINS NN,PT ECM: ',3f12.4/
891  * ' CC CHAINS NN,PT ECM: ',3f12.4)
892  ENDIF
893 C
894 C DISTRIBUTION FOR
895 C NUMBER OF GENERATIONS INSIDE THE NUCLEUS
896 C (NORMALIZED TO 1)
897  avgnor=0.0
898  avgen=0.0
899  DO 801 ige=1,50
900  avgnor=avgnor + disgen(ige)
901  avgen=avgen + ige*disgen(ige)
902  801 CONTINUE
903  DO 802 ige=1,50
904  xgen(ige)=float(ige) - 0.5
905  disgen(ige)=disgen(ige)/avgnor
906  802 CONTINUE
907  avgen=avgen/avgnor
908 C WRITE(*,'(A,1PE10.2//A)')
909 C & ' AVERAGE NUMBER OF GENERATIONS INSIDE THE NUCLEUS =',
910 C & AVGEN,
911 C & ' CORRESPONDING DISTRIBUTION (NORMALIZED TO 1.0)'
912 C CALL PLOT(XGEN,DISGEN,50,1,50,0.0,1.0,0.0,0.01)
913 C
914  DO 310 i=1,numtyp
915  avmult(i)=avmult(i)/nhkkh1
916  ave(i)=ave(i)/nhkkh1
917  avept(i)=avept(i)/(nhkkh1*avmult(i))
918 310 CONTINUE
919  ffff2=sqrt(anchsq-avmult(27)**2)/avmult(27)
920  WRITE(6,7772)ffff2,anchsq,avmult(27)
921  7772 FORMAT(' FFFF2,ANCHSQ,AVMULT(27):',3e15.5)
922  WRITE(6, 64)
923  64 FORMAT(' PARTICLE REF,CHAR,IBAR, MASS AVERAGE',
924  *' ENERGY, MULTIPLICITY, INELASTICITY')
925  DO 62 i=1,numtyp
926  ake(i)=ave(i)/cmener
927  WRITE(6, 63) aname(i),i,iich(i),iibar(i),aam(i),
928  * ave(i),avmult(i),ake(i),avept(i)
929  63 FORMAT (' ',a8,3i5,f10.3,4e15.5)
930  62 CONTINUE
931 C
932  DO 33 i=1,20
933  DO 35 j=1,50
934  yyl(j,i) =yyl(j,i) /(nhkkh1*dy)
935  yylps(j,i)=yylps(j,i)/(nhkkh1*dy)
936 35 CONTINUE
937 33 CONTINUE
938  WRITE(6,'(1H1,11(A/))')
939  & ' Conventions for rapidity distributions: ',
940  & ' 1 Proton 11 NEGATIVES ',
941  & ' 2 Neutron 12 NBAR=9 ',
942  & ' 3 PI+=13 13 LAMBDA=17 ',
943  & ' 4 PI-=14 14 SIGMA==20,21,22 ',
944  & ' 5 K+ =15 15 LAMBDABAR=18 ',
945  & ' 6 K- =16 16 SIGMABAR=99,100,101',
946  & ' 7 neutral kaons=12,19,24,25 17 THETA=97,98 ',
947  & ' 8 pbar=2 18 THETABAR=102,103 ',
948  & ' 9 charged hadrons ',
949  & ' 10 total hadrons '
950  WRITE(6, 66)
951  66 FORMAT(' RAPIDITY DISTRIBUTION')
952  WRITE(6,302)
953  302 FORMAT (' (first number gives the lower bin limit)')
954  DO 36 j=1,50
955  WRITE(6, 37) xyl(j,1),(yyl(j,i),i=1,10)
956 37 FORMAT (f10.2,10e11.3)
957 36 CONTINUE
958  WRITE(6, 66)
959  WRITE(6,302)
960  DO 306 j=1,50
961  WRITE(6, 37) xyl(j,1),(yyl(j,i),i=11,20)
962 306 CONTINUE
963 C
964  WRITE(6, 66 )
965  CALL plot(xyl,yyl,1000,20,50,-10.,dy,0.,delrap)
966 * rescaled rapidity distribution
967  CALL plot(xyl,yyl,1000,20,50,-10.,dy,0.,5.*delrap)
968 * logarithmic rapidity distribution
969  DO 38 i=1,20
970  DO 38 j=1,50
971  yyl(j,i)=log10(abs(yyl(j,i))+1.e-8)
972  38 CONTINUE
973  WRITE(6, 66)
974  CALL plot(xyl,yyl,1000,20,50,-10.,dy,-2.0,0.05)
975 C
976  WRITE(6,301)
977  301 FORMAT ('1 PSEUDORAPIDITY DISTRIBUTION')
978  WRITE(6,302)
979  DO 303 j=1,50
980  WRITE(6,37) xyl(j,1),(yylps(j,i),i=1,10)
981  303 CONTINUE
982  WRITE(6,301)
983  WRITE(6,302)
984  DO 304 j=1,50
985  WRITE(6,37) xyl(j,1),(yylps(j,i),i=11,20)
986  304 CONTINUE
987  WRITE(6,301)
988  CALL plot(xyl,yylps,1000,20,50,-10.,dy,0.,delrap)
989  CALL plot(xyl,yylps,1000,20,50,-10.,dy,0.,5.*delrap)
990  RETURN
991  END
992 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
993 C
994  SUBROUTINE distrp(IOP,NHKKH1,PO)
995  IMPLICIT DOUBLE PRECISION (a-h,o-z)
996  RETURN
997  END
998 C--------------------------------------------Dummies----------
999  SUBROUTINE distco(IOP,IJPROJ,PPN,IDUMMY)
1000  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1001  RETURN
1002  END
1003  SUBROUTINE distpa(IOP)
1004  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1005  RETURN
1006  END
1007  SUBROUTINE disres(IOP,IJPROJ,PPN)
1008  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1009  RETURN
1010  END
1011  SUBROUTINE edensi(I,J)
1012  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1013  RETURN
1014  END
1015  SUBROUTINE plomb(I,PP,CHAR,XF,ITIF,IJPROJ)
1016  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1017  RETURN
1018  END
1019 
1020  SUBROUTINE plombc(I,PP,CHAR,XF,ITIF,IJPROJ)
1021  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1022  RETURN
1023  END
1024 
1025  SUBROUTINE dispt(IOP,NHKKH1,PO)
1026  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1027 C
1028 C*** 1=P, 2=N, 3=PI+, 4=PI-, 5=PIO, 6=GAM+HYP, 7=K, 8=ANUC, 9=CHARGED
1029 C*** 10=TOT, 11=TOTHAD
1030 C
1031 *KEEP,HKKEVT.
1032 c INCLUDE (HKKEVT)
1033  parameter(nmxhkk= 89998)
1034 c PARAMETER (NMXHKK=25000)
1035  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
1036  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
1037  +(4,nmxhkk)
1038 C
1039 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
1040 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
1041 C THE POSITIONS OF THE PROJECTILE NUCLEONS
1042 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
1043 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
1044 C COMPLETELY CONSISTENT. THE TIMES IN THE
1045 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
1046 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
1047 C
1048 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
1049 C
1050 C NMXHKK: maximum numbers of entries (partons/particles) that can be
1051 C stored in the commonblock.
1052 C
1053 C NHKK: the actual number of entries stored in current event. These are
1054 C found in the first NHKK positions of the respective arrays below.
1055 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
1056 C entry.
1057 C
1058 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
1059 C = 0 : null entry.
1060 C = 1 : an existing entry, which has not decayed or fragmented.
1061 C This is the main class of entries which represents the
1062 C "final state" given by the generator.
1063 C = 2 : an entry which has decayed or fragmented and therefore
1064 C is not appearing in the final state, but is retained for
1065 C event history information.
1066 C = 3 : a documentation line, defined separately from the event
1067 C history. (incoming reacting
1068 C particles, etc.)
1069 C = 4 - 10 : undefined, but reserved for future standards.
1070 C = 11 - 20 : at the disposal of each model builder for constructs
1071 C specific to his program, but equivalent to a null line in the
1072 C context of any other program. One example is the cone defining
1073 C vector of HERWIG, another cluster or event axes of the JETSET
1074 C analysis routines.
1075 C = 21 - : at the disposal of users, in particular for event tracking
1076 C in the detector.
1077 C
1078 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
1079 C standard.
1080 C
1081 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
1082 C The value is 0 for initial entries.
1083 C
1084 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
1085 C one mother exist, in which case the value 0 is used. In cluster
1086 C fragmentation models, the two mothers would correspond to the q
1087 C and qbar which join to form a cluster. In string fragmentation,
1088 C the two mothers of a particle produced in the fragmentation would
1089 C be the two endpoints of the string (with the range in between
1090 C implied).
1091 C
1092 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
1093 C entry has not decayed, this is 0.
1094 C
1095 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
1096 C entry has not decayed, this is 0. It is assumed that the daughters
1097 C of a particle (or cluster or string) are stored sequentially, so
1098 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
1099 C daughters. Even in cases where only one daughter is defined (e.g.
1100 C K0 -> K0S) both values should be defined, to make for a uniform
1101 C approach in terms of loop constructions.
1102 C
1103 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
1104 C
1105 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
1106 C
1107 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
1108 C
1109 C PHKK(4,IHKK) : energy, in GeV.
1110 C
1111 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
1112 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
1113 C
1114 C VHKK(1,IHKK) : production vertex x position, in mm.
1115 C
1116 C VHKK(2,IHKK) : production vertex y position, in mm.
1117 C
1118 C VHKK(3,IHKK) : production vertex z position, in mm.
1119 C
1120 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
1121 C********************************************************************
1122 *KEEP,NSHMAK.
1123  COMMON /nshmak/ nnshma,npshma,ntshma,nshmac
1124 *KEEP,DPAR.
1125 C /DPAR/ CONTAINS PARTICLE PROPERTIES
1126 C ANAME = LITERAL NAME OF THE PARTICLE
1127 C AAM = PARTICLE MASS IN GEV
1128 C GA = DECAY WIDTH
1129 C TAU = LIFE TIME OF INSTABLE PARTICLES
1130 C IICH = ELECTRIC CHARGE OF THE PARTICLE
1131 C IIBAR = BARYON NUMBER
1132 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
1133 C
1134  CHARACTER*8 aname
1135  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
1136  +iibar(210),k1(210),k2(210)
1137 C------------------
1138 *KEEP,NUCC.
1139  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
1140 *KEEP,DPRIN.
1141  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1142 *KEND.
1143 
1144 C modified DPMJET
1145  COMMON /bufues/ bnnvv,bnnss,bnnsv,bnnvs,bnncc,
1146  * bnndv,bnnvd,bnnds,bnnsd,
1147  * bnnhh,bnnzz,
1148  * bptvv,bptss,bptsv,bptvs,bptcc,bptdv,
1149  * bptvd,bptds,bptsd,
1150  * bpthh,bptzz,
1151  * beevv,beess,beesv,beevs,beecc,beedv,
1152  * beevd,beeds,beesd,
1153  * beehh,beezz
1154  * ,bnndi,bptdi,beedi
1155  * ,bnnzd,bnndz,bptzd,bptdz,beezd,beedz
1156  COMMON /ncoucs/ bcouvv,bcouss,bcousv,bcouvs,
1157  * bcouzz,bcouhh,bcouds,bcousd,
1158  * bcoudz,bcouzd,bcoudi,
1159  * bcoudv,bcouvd,bcoucc
1160  COMMON /bufueh/ annvv,annss,annsv,annvs,anncc,
1161  * anndv,annvd,annds,annsd,
1162  * annhh,annzz,
1163  * ptvv,ptss,ptsv,ptvs,ptcc,ptdv,ptvd,ptds,ptsd,
1164  * pthh,ptzz,
1165  * eevv,eess,eesv,eevs,eecc,eedv,eevd,eeds,eesd,
1166  * eehh,eezz
1167  * ,anndi,ptdi,eedi
1168  * ,annzd,anndz,ptzd,ptdz,eezd,eedz
1169  COMMON /ncouch/ acouvv,acouss,acousv,acouvs,
1170  * acouzz,acouhh,acouds,acousd,
1171  * acoudz,acouzd,acoudi,
1172  * acoudv,acouvd,acoucc
1173 C---------------------
1174  COMMON /eventa/idumtp
1175 C
1176  dimension ptp(50,20),pty(50,20),emty(50,20),indx(28)
1177  COMMON /sigla/siglau
1178  COMMON /final/ifinal
1179 C------------------------------------------------------------------
1180 C------------------------------------------------------------------
1181 
1182  DATA indx/1,8,10,10,10,10,7,2,7,10,10,7,3,4,5,6,
1183  * 7,7,7,7,7,7,7,7,7,7,7,7/
1184 C kinematics in lab (target rest system)
1185 C***
1186 C conventions for rapidity distributions:
1187 C 1 Proton 11 NEGATIVES
1188 C 2 Neutron 12 NBAR=9
1189 C 3 PI+=13 13 LAMBDA=17
1190 C 4 PI-=14 14 SIGMA==20,21,22
1191 C 5 K+ =15 15 LAMBDABAR=18
1192 C 6 K- =16 16 SIGMABAR=99,100,101
1193 C 7 neutral kaons=12,19,24,25 17 THETA=97,98
1194 C 8 pbar=2 18 THETABAR=102,103
1195 C 9 charged hadrons 19 neutral Kaons 12,19,24,25
1196 C 10 total hadrons
1197 C-----------------------------------------------------------------------
1198  DATA kpl /1/
1199  DATA ievl /0/
1200 C----------------------------------------------------------------------
1201  zero=0
1202  oneone=1
1203  twotwo=2
1204  nip=ip
1205  aip=ip
1206  delrap=1.
1207  IF(ip.EQ.1)delrap=0.1
1208 C
1209 C-----------------------------
1210  go to(1,2,3),iop
1211 1 CONTINUE
1212  dxfl=0.02
1213  dy=0.4
1214  dpt=0.10
1215 C DIMENSION PTP(50,20),PTY(50,20),EMTY(50,20)
1216  DO 102 i=1,20
1217  DO 102 j=1,50
1218  ptp(j,i)=(j-1)*dpt
1219  pty(j,i)=1.d-12
1220  emty(j,i)=1.d-12
1221  102 CONTINUE
1222  nhkkh2=nhkkh1
1223  IF (nhkkh1.EQ.0)nhkkh2=1
1224  eeo=sqrt(po**2+aam(nhkkh2)**2)
1225  WRITE(6, 1001)eeo,po,nhkkh2,aam(nhkkh2)
1226  1001 FORMAT (' EEO,PO ',f13.2,f13.2,i10,f10.2)
1227 2 CONTINUE
1228  nhad=0
1229 C HBOOK HISTOGRAMS
1230 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1231 C WRITE(6,'(A)')' before plomb'
1232  IF(ihbook.EQ.1)CALL plomb(2,p4p4p4,ccchrg,xfxfxf,1,ijproj)
1233 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1234  DO 21 i=nhkkh1,nhkk
1235  IF (isthkk(i).EQ.1)THEN
1236 C j.r. temp for s-s
1237  ptt=phkk(1,i)**2+phkk(2,i)**2+0.00001
1238  pt=sqrt(ptt)+0.001
1239 C IF(PT.LT.0.5)GO TO 21
1240 C
1241  nhad=nhad+1
1242  nrhkk=mcihad(idhkk(i))
1243  IF (nrhkk.LE.0.OR.nrhkk.GT.210)THEN
1244  nrhkk=1
1245  ENDIF
1246  nre=nrhkk
1247  ichhkk=iich(nrhkk)
1248  IF (nre.GT.160) nre=28
1249  IF (nre.LT. 1) nre=28
1250  nrex=nre
1251  IF(nrex.GE.28)nrex=28
1252  ni=indx(nrex)
1253  nix=ni
1254  IF (nrhkk.EQ.9)nix=12
1255  IF (nrhkk.EQ.17.OR.nrhkk.EQ.22)nix=13
1256  IF (nrhkk.LE.22.AND.nrhkk.GE.20)nix=14
1257  IF (nrhkk.EQ.18.OR.nrhkk.EQ.100)nix=15
1258  IF (nrhkk.LE.101.AND.nrhkk.GE.99)nix=16
1259  IF (nrhkk.EQ.98)nix=17
1260  IF (nrhkk.EQ.103)nix=18
1261  IF (nrhkk.EQ.12.OR.nrhkk.EQ.19)nix=19
1262  IF (nrhkk.EQ.24.OR.nrhkk.EQ.25)nix=19
1263  IF(nre.EQ.28) ni=7
1264  ptt=phkk(1,i)**2+phkk(2,i)**2+0.00001
1265  pt=sqrt(ptt)+0.001
1266  amt=sqrt(ptt+phkk(5,i)**2)
1267  yl=log((abs(phkk(3,i)+sqrt(phkk(3,i)**2+amt**2)))/amt+1.e-18)
1268  yllps=log((abs(phkk(3,i)+sqrt(phkk(3,i)**2+ptt)))/sqrt(ptt)
1269  * +1.e-18)
1270  kpl=1
1271  ipt=pt/dpt+1.
1272  IF (ipt.LT.1)ipt=1
1273  IF (ipt.GT.50) ipt=50
1274  IF (ichhkk.NE.0)pty(ipt,9)=pty(ipt,9)+1./pt
1275  IF (ichhkk.EQ.-1)pty(ipt,11)=pty(ipt,11)+1./pt
1276  pty(ipt,nix)=pty(ipt,nix)+1./pt
1277  pty(ipt,10)=pty(ipt,10)+1./pt
1278  IF(yl.GT.2.3.AND.yl.LE.3.)THEN
1279  iamt=amt/dpt+1.
1280  IF (iamt.LT.1)iamt=1
1281  IF (iamt.GT.48) iamt=48
1282  IF (ichhkk.NE.0)emty(iamt,9)=emty(iamt,9)+1./amt
1283  IF (ichhkk.EQ.-1)emty(iamt,11)=emty(iamt,11)+1./amt
1284  emty(iamt,nix)=emty(iamt,nix)+1./amt
1285  emty(iamt,10)=emty(iamt,10)+1./amt
1286  IF(pt.GT.1.0.AND.pt.LT.2.0)THEN
1287  iamt=49
1288  IF (ichhkk.NE.0)emty(iamt,9)=emty(iamt,9)+dpt
1289  IF (ichhkk.EQ.-1)emty(iamt,11)=emty(iamt,11)+dpt
1290  emty(iamt,nix)=emty(iamt,nix)+dpt
1291  emty(iamt,10)=emty(iamt,10)+dpt
1292  ENDIF
1293  IF(amt.GT.1.72.AND.yl.LE.2.6)THEN
1294  iamt=50
1295  IF (ichhkk.NE.0)emty(iamt,9)=emty(iamt,9)+dpt
1296  IF (ichhkk.EQ.-1)emty(iamt,11)=emty(iamt,11)+dpt
1297  emty(iamt,nix)=emty(iamt,nix)+dpt
1298  emty(iamt,10)=emty(iamt,10)+dpt
1299  ENDIF
1300  ENDIF
1301  ENDIF
1302 21 CONTINUE
1303  RETURN
1304 C
1305 C--------------------------------------------------------------------
1306 C
1307  3 CONTINUE
1308 C Normalization and Printing
1309 C------------------------------------------------------------------
1310 C------------------------------------------------------------------
1311  kkk=13
1312 C Normalization
1313 C
1314 C NORMALIZE AND PRINT COUNTERS
1315  IF (bnndi.GT.1.)THEN
1316  bptdi=bptdi/bnndi
1317  beedi=beedi/bnndi
1318  ENDIF
1319  IF (bnnvv.GT.1.)THEN
1320  bptvv=bptvv/bnnvv
1321  beevv=beevv/bnnvv
1322  ENDIF
1323  IF (bnnss.GT.1.)THEN
1324  bptss=bptss/bnnss
1325  beess=beess/bnnss
1326  ENDIF
1327  IF (bnnsv.GT.1.)THEN
1328  bptsv=bptsv/bnnsv
1329  beesv=beesv/bnnsv
1330  ENDIF
1331  IF (bnnvs.GT.1.)THEN
1332  bptvs=bptvs/bnnvs
1333  beevs=beevs/bnnvs
1334  ENDIF
1335  IF (bnncc.GE.1.)THEN
1336  bptcc=bptcc/bnncc
1337  beecc=beecc/bnncc
1338  ENDIF
1339  IF (bnndv.GE.1.)THEN
1340  bptdv=bptdv/bnndv
1341  beedv=beedv/bnndv
1342  ENDIF
1343  IF (bnnvd.GE.1.)THEN
1344  bptvd=bptvd/bnnvd
1345  beevd=beevd/bnnvd
1346  ENDIF
1347  IF (bnnds.GE.1.)THEN
1348  bptds=bptds/bnnds
1349  beeds=beeds/bnnds
1350  ENDIF
1351  IF (bnndz.GE.1.)THEN
1352  bptdz=bptdz/bnndz
1353  beedz=beedz/bnndz
1354  ENDIF
1355  IF (bnnhh.GE.1.)THEN
1356  bpthh=bpthh/bnnhh
1357  beehh=beehh/bnnhh
1358  ENDIF
1359  IF (bnnzz.GE.1.)THEN
1360  bptzz=bptzz/bnnzz
1361  beezz=beezz/bnnzz
1362  ENDIF
1363  IF (bnnsd.GE.1.)THEN
1364  bptsd=bptsd/bnnsd
1365  beesd=beesd/bnnsd
1366  ENDIF
1367  IF (bnnzd.GE.1.)THEN
1368  bptzd=bptzd/bnnzd
1369  beezd=beezd/bnnzd
1370  ENDIF
1371  IF (nhkkh1.GT.1.)THEN
1372  bnnvv=bnnvv/nhkkh1
1373  bnndi=bnndi/nhkkh1
1374  bnnss=bnnss/nhkkh1
1375  bnnsv=bnnsv/nhkkh1
1376  bnnvs=bnnvs/nhkkh1
1377  bnncc=bnncc/nhkkh1
1378  bnnvd=bnnvd/nhkkh1
1379  bnndv=bnndv/nhkkh1
1380  bnnds=bnnds/nhkkh1
1381  bnnsd=bnnsd/nhkkh1
1382  bnndz=bnndz/nhkkh1
1383  bnnzd=bnnzd/nhkkh1
1384  bnnhh=bnnhh/nhkkh1
1385  bnnzz=bnnzz/nhkkh1
1386  bcouvv=bcouvv/nhkkh1
1387  bcouss=bcouss/nhkkh1
1388  bcousv=bcousv/nhkkh1
1389  bcouvs=bcouvs/nhkkh1
1390  bcouzz=bcouzz/nhkkh1
1391  bcouhh=bcouhh/nhkkh1
1392  bcouds=bcouds/nhkkh1
1393  bcousd=bcousd/nhkkh1
1394  bcoudz=bcoudz/nhkkh1
1395  bcouzd=bcouzd/nhkkh1
1396  bcoudi=bcoudi/nhkkh1
1397  bcoudv=bcoudv/nhkkh1
1398  bcouvd=bcouvd/nhkkh1
1399  bcoucc=bcoucc/nhkkh1
1400  WRITE (6,7431)bnnvv,bptvv,beevv,bcouvv,
1401  * bnnss,bptss,beess,bcouss,
1402  * bnnsv,bptsv,beesv,bcousv,
1403  * bnnvs,bptvs,beevs,bcouvs,
1404  * bnncc,bptcc,beecc,bcoucc,
1405  * bnndv,bptdv,beedv,bcoudv,
1406  * bnnvd,bptvd,beevd,bcouvd,
1407  * bnnds,bptds,beeds,bcouds,
1408  * bnnsd,bptsd,beesd,bcousd,
1409  * bnndz,bptdz,beedz,bcoudz,
1410  * bnnzd,bptzd,beezd,bcouzd,
1411  * bnnhh,bpthh,beehh,bcouhh,
1412  * bnndi,bptdi,beedi,bcoudi,
1413  * bnnzz,bptzz,beezz,bcouzz
1414  7431 FORMAT (' VV CHAINS NN,PT ECM: ',4f12.4/
1415  * ' SS CHAINS NN,PT ECM: ',4f12.4/
1416  * ' SV CHAINS NN,PT ECM: ',4f12.4/
1417  * ' VS CHAINS NN,PT ECM: ',4f12.4/
1418  * ' CC CHAINS NN,PT ECM: ',4f12.4/
1419  * ' DV CHAINS NN,PT ECM: ',4f12.4/
1420  * ' VD CHAINS NN,PT ECM: ',4f12.4/
1421  * ' DS CHAINS NN,PT ECM: ',4f12.4/
1422  * ' SD CHAINS NN,PT ECM: ',4f12.4/
1423  * ' DZ CHAINS NN,PT ECM: ',4f12.4/
1424  * ' ZD CHAINS NN,PT ECM: ',4f12.4/
1425  * ' HH CHAINS NN,PT ECM: ',4f12.4/
1426  * ' DI CHAINS NN,PT ECM: ',4f12.4/
1427  * ' ZZ CHAINS NN,PT ECM: ',4f12.4)
1428  ENDIF
1429 C
1430 C
1431 C
1432 C============================================================
1433  DO 33 i=1,20
1434  DO 35 j=1,50
1435  DO 335 jjjj=0,1
1436  335 CONTINUE
1437  pty(j,i) =pty(j,i) /(nhkkh1*dpt)
1438  emty(j,i) =emty(j,i) /(nhkkh1*dpt)
1439 35 CONTINUE
1440 33 CONTINUE
1441  DO 5136 j=1,50
1442  WRITE(6, 37)(pty(j,i),i=1,11)
1443 5137 FORMAT (11e11.3)
1444  37 FORMAT (11e11.3)
1445 5136 CONTINUE
1446  WRITE(6,3654)
1447  3654 FORMAT(' pt-Distribution')
1448  DO 3614 j=1,50
1449  WRITE(6, 37)ptp(j,1),(pty(j,i),i=1,10)
1450  3614 CONTINUE
1451  DO 5914 j=1,50
1452  WRITE(6, 37)ptp(j,1),(pty(j,i),i=11,20)
1453  5914 CONTINUE
1454  WRITE(6,4654)
1455  4654 FORMAT(' Et-Distribution')
1456  WRITE(6, 37)ptp(1,1),(emty(1,i),i=1,10)
1457  DO 4614 j=2,50
1458  WRITE(6, 37)ptp(j-1,1),(emty(j,i),i=1,10)
1459  WRITE(6, 37)ptp(j,1),(emty(j,i),i=1,10)
1460  4614 CONTINUE
1461  WRITE(6, 37)ptp(1,1),(emty(1,i),i=11,20)
1462  DO 6914 j=2,50
1463  WRITE(6, 37)ptp(j-1,1),(emty(j,i),i=11,20)
1464  WRITE(6, 37)ptp(j,1),(emty(j,i),i=11,20)
1465  6914 CONTINUE
1466  DO 514 j=1,50
1467  DO 312 i=1,20
1468  pty(j,i) =log10(abs(pty(j,i))+1.e-8)
1469  emty(j,i) =log10(abs(emty(j,i))+1.e-8)
1470 312 CONTINUE
1471  pty(j,10) =pty(j,11)
1472  emty(j,10) =emty(j,11)
1473 514 CONTINUE
1474  WRITE(6,3654)
1475  CALL plot(ptp,pty,1000,20,50,zero,dpt,-oneone,0.05d0)
1476 C CALL PLOT(PTP,PTY,1000,20,50,ZERO,DPT,ONEONE,0.05D0)
1477  WRITE(6,4654)
1478 
1479  CALL plot(ptp,emty,1000,20,50,zero,dpt,-6.0d0,0.10d0)
1480 C
1481  RETURN
1482  END
1483 C
1484 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1485 C
1486  SUBROUTINE histog(I)
1487 C
1488  RETURN
1489  END
1490 C----------------------------------------------------------------
1491  SUBROUTINE diseva(IOP,NHKKH1,PO,IGENER)
1492  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1493 *KEEP,HKKEVT.
1494 c INCLUDE (HKKEVT)
1495  parameter(nmxhkk= 89998)
1496 c PARAMETER (NMXHKK=25000)
1497  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
1498  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
1499  +(4,nmxhkk)
1500  COMMON /extevt/ idres(nmxhkk),idxres(nmxhkk),nobam(nmxhkk),
1501  & idbam(nmxhkk),idch(nmxhkk),npoint(10)
1502 C
1503 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
1504 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
1505 C THE POSITIONS OF THE PROJECTILE NUCLEONS
1506 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
1507 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
1508 C COMPLETELY CONSISTENT. THE TIMES IN THE
1509 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
1510 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
1511 C
1512 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
1513 C
1514 C NMXHKK: maximum numbers of entries (partons/particles) that can be
1515 C stored in the commonblock.
1516 C
1517 C NHKK: the actual number of entries stored in current event. These are
1518 C found in the first NHKK positions of the respective arrays below.
1519 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
1520 C entry.
1521 C
1522 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
1523 C = 0 : null entry.
1524 C = 1 : an existing entry, which has not decayed or fragmented.
1525 C This is the main class of entries which represents the
1526 C "final state" given by the generator.
1527 C = 2 : an entry which has decayed or fragmented and therefore
1528 C is not appearing in the final state, but is retained for
1529 C event history information.
1530 C = 3 : a documentation line, defined separately from the event
1531 C history. (incoming reacting
1532 C particles, etc.)
1533 C = 4 - 10 : undefined, but reserved for future standards.
1534 C = 11 - 20 : at the disposal of each model builder for constructs
1535 C specific to his program, but equivalent to a null line in the
1536 C context of any other program. One example is the cone defining
1537 C vector of HERWIG, another cluster or event axes of the JETSET
1538 C analysis routines.
1539 C = 21 - : at the disposal of users, in particular for event tracking
1540 C in the detector.
1541 C
1542 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
1543 C standard.
1544 C
1545 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
1546 C The value is 0 for initial entries.
1547 C
1548 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
1549 C one mother exist, in which case the value 0 is used. In cluster
1550 C fragmentation models, the two mothers would correspond to the q
1551 C and qbar which join to form a cluster. In string fragmentation,
1552 C the two mothers of a particle produced in the fragmentation would
1553 C be the two endpoints of the string (with the range in between
1554 C implied).
1555 C
1556 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
1557 C entry has not decayed, this is 0.
1558 C
1559 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
1560 C entry has not decayed, this is 0. It is assumed that the daughters
1561 C of a particle (or cluster or string) are stored sequentially, so
1562 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
1563 C daughters. Even in cases where only one daughter is defined (e.g.
1564 C K0 -> K0S) both values should be defined, to make for a uniform
1565 C approach in terms of loop constructions.
1566 C
1567 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
1568 C
1569 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
1570 C
1571 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
1572 C
1573 C PHKK(4,IHKK) : energy, in GeV.
1574 C
1575 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
1576 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
1577 C
1578 C VHKK(1,IHKK) : production vertex x position, in mm.
1579 C
1580 C VHKK(2,IHKK) : production vertex y position, in mm.
1581 C
1582 C VHKK(3,IHKK) : production vertex z position, in mm.
1583 C
1584 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
1585 C********************************************************************
1586 C------------------------------------------------------------------
1587 C--------------------------------evaporation-----------------------
1588  dimension aneva(4),pevap(50,2),xpevap(50,2),amevap(250),
1589  * xmevap(250)
1590  dimension anevap(4),amevpp(250)
1591  COMMON /final/ifinal
1592  COMMON /nomije/ ptmije(10),nnmije(10)
1593  COMMON /nomiju/ nnmiju(10)
1594  CHARACTER*8 aname
1595  COMMON /dpar/aname(210),aam(210),ga(210),tau(210),iich(210),
1596  +iibar(210),k1(210),k2(210)
1597 C------------------------------------------------------------------
1598 C------------------------------------------------------------------
1599  go to(1,2,3),iop
1600 C------------------------------------------------------------------
1601 C--------------------------------evaporation-----------------------
1602  1 CONTINUE
1603  dpeva=0.02d0
1604  dit=it/50+1
1605  IF(dit.LT.1.d0)dit=1
1606  IF (ifinal.EQ.0) THEN
1607  DO 3112 i=1,4
1608  aneva(i)=0
1609  anevap(i)=0
1610  3112 CONTINUE
1611  DO 3123 i=1,250
1612  amevap(i)=0
1613  amevpp(i)=0
1614  xmevap(i)=i
1615  3123 CONTINUE
1616  DO 3113 i=1,50
1617  DO 3113 ii=1,2
1618  pevap(i,ii)=0
1619  xpevap(i,ii)=i*dpeva
1620  3113 CONTINUE
1621  ENDIF
1622 C s(shower), g (grey), h(heavy) and b (black) particles
1623 C n- (neg prongs) sp8,sp5 (slow protons) p=0.15-0.8, 0.15-0.5
1624  ansho=0
1625  angre=0
1626  angre2=0
1627  anhea=0
1628  anbla=0
1629  anmin=0
1630  ansp8=0
1631  ansp5=0
1632  RETURN
1633 C------------------------------------------------------------------
1634 C------------------------------------------------------------------
1635  2 CONTINUE
1636 C------------------------------------------------------------------
1637 C--------------------------------evaporation-----------------------
1638 C WRITE(6,'(A)')' DISTR 2 '
1639  DO 1121 i=nhkkh1,nhkk
1640  IF (isthkk(i).EQ.1)THEN
1641  nrhkk=mcihad(idhkk(i))
1642  ichhkk=iich(nrhkk)
1643  pptt=sqrt(phkk(1,i)**2+phkk(2,i)**2+0.000001)
1644  ptt=pptt**2
1645  ptot=sqrt(pptt**2+phkk(3,i)**2+0.000001)
1646  IF(ichhkk.NE.0)THEN
1647  betp=ptot/phkk(4,i)
1648  IF(ichhkk.EQ.-1)anmin=anmin+1
1649  IF(betp.GE.0.7d0)ansho=ansho+1
1650  IF(betp.LE.0.7d0)anhea=anhea+1
1651  IF(betp.LE.0.2d0)anbla=anbla+1
1652  IF(betp.GE.0.2d0.AND.betp.LE.0.7d0)angre=angre+1
1653  IF(betp.GE.0.23d0.AND.betp.LE.0.7d0)angre2=angre2+1
1654  ENDIF
1655  ENDIF
1656  1121 CONTINUE
1657  IF (ifinal.EQ.0) THEN
1658  DO 3114 i=nhkkh1,nhkk
1659  IF(isthkk(i).EQ.-1) THEN
1660 C Neutrons
1661  IF (idhkk(i).EQ.2112)THEN
1662  pptt=(phkk(1,i)**2+phkk(2,i)**2+0.000001)
1663  ptot=sqrt(pptt+phkk(3,i)**2+0.000001)
1664  IF(nobam(i).EQ.2)THEN
1665  aneva(2)=aneva(2)+1.d0
1666  iptot=ptot/dpeva
1667  IF(iptot.LT.1)iptot=1
1668  IF(iptot.GT.50)iptot=50
1669  pevap(iptot,2)=pevap(iptot,2)+1.e0
1670  ELSEIF(nobam(i).EQ.1)THEN
1671  anevap(2)=anevap(2)+1.d0
1672  ENDIF
1673 C Protons
1674  ELSEIF(idhkk(i).EQ.2212)THEN
1675  pptt=(phkk(1,i)**2+phkk(2,i)**2+0.000001)
1676  ptot=sqrt(pptt+phkk(3,i)**2+0.000001)
1677  IF(nobam(i).EQ.2)THEN
1678  aneva(1)=aneva(1)+1.d0
1679  IF(ptot.GT.0.15d0.AND.ptot.LE.0.8d0)ansp8=ansp8+1
1680  IF(ptot.GT.0.15d0.AND.ptot.LE.0.5d0)ansp5=ansp5+1
1681  iptot=ptot/dpeva
1682  IF(iptot.LT.1)iptot=1
1683  IF(iptot.GT.50)iptot=50
1684  pevap(iptot,1)=pevap(iptot,1)+1.e0
1685  betp=ptot/phkk(4,i)
1686  IF(betp.GE.0.7d0)ansho=ansho+1
1687  IF(betp.LE.0.7d0)anhea=anhea+1
1688  IF(betp.LE.0.23d0)anbla=anbla+1
1689  IF(betp.GE.0.2d0.AND.betp.LE.0.7d0)angre=angre+1
1690 C IF(BETP.GE.0.23D0.AND.BETP.LE.0.7D0)
1691 C * ANGRE2=ANGRE2+1
1692  IF(ptot.GE.0.026d0.AND.ptot.LE.0.375d0)
1693  * angre2=angre2+1
1694  ELSEIF(nobam(i).EQ.1)THEN
1695  anevap(1)=anevap(1)+1.d0
1696  ENDIF
1697 C Photons
1698  ELSEIF(idhkk(i).EQ.22)THEN
1699  pptt=(phkk(1,i)**2+phkk(2,i)**2+0.000001)
1700  ptot=sqrt(pptt+phkk(3,i)**2+0.000001)
1701  IF(nobam(i).EQ.2)THEN
1702  aneva(3)=aneva(3)+1.d0
1703  ELSEIF(nobam(i).EQ.1)THEN
1704  anevap(3)=anevap(3)+1.d0
1705  ENDIF
1706 C transform deexcitation photons into
1707 C electron neutrinos (pdg=12, bamjet=5
1708 C to allow separate histograms of pi-0
1709 C and deexcitation photons
1710  isthkk(i)=1
1711  idhkk(i) =12
1712  ENDIF
1713  ENDIF
1714  IF(idhkk(i).EQ.80000)THEN
1715 C heavy fragments
1716  pptt=(phkk(1,i)**2+phkk(2,i)**2+0.000001)
1717  ptot=sqrt(pptt+phkk(3,i)**2+0.000001)
1718  IF(nobam(i).EQ.2)THEN
1719  aneva(4)=aneva(4)+1.d0
1720  idit=idres(i)
1721  IF(idit.EQ.0)go to 3115
1722  amevap(idit)=amevap(idit)+1.d0
1723  anhea=anhea+1
1724  anbla=anbla+1
1725  ELSEIF(nobam(i).EQ.1)THEN
1726  anevap(4)=anevap(4)+1.d0
1727  idit=idres(i)
1728  IF(idit.EQ.0)go to 3115
1729  amevpp(idit)=amevpp(idit)+1.d0
1730  ENDIF
1731  3115 CONTINUE
1732  ENDIF
1733  3114 CONTINUE
1734  ENDIF
1735  RETURN
1736 C------------------------------------------------------------------
1737 C------------------------------------------------------------------
1738  3 CONTINUE
1739 C------------------------------------------------------------------
1740 C--------------------------------evaporation-----------------------
1741  IF (ifinal.EQ.0) THEN
1742  WRITE(6,'(A)')' Target fragments'
1743  DO 3117 i=1,4
1744  aneva(i)=aneva(i)/nhkkh1
1745  3117 CONTINUE
1746  WRITE(6,'(A,F10.3)')' Number of evap. protons: ',aneva(1)
1747  WRITE(6,'(A,F10.3)')' Number of evap. neutrans: ',aneva(2)
1748  WRITE(6,'(A,F10.3)')' Number of ex. gammas: ',aneva(3)
1749  WRITE(6,'(A,F10.3)')' Number of heavy fragments:',aneva(4)
1750  WRITE(6,'(A)')' Projectile fragments'
1751  DO 3118 i=1,4
1752  anevap(i)=anevap(i)/nhkkh1
1753  3118 CONTINUE
1754  WRITE(6,'(A,F10.3)')
1755  * ' Number of evap. protons: ',anevap(1)
1756  WRITE(6,'(A,F10.3)')
1757  * ' Number of evap. neutrans: ',anevap(2)
1758  WRITE(6,'(A,F10.3)')
1759  * ' Number of ex. gammas: ',anevap(3)
1760  WRITE(6,'(A,F10.3)')
1761  * ' Number of heavy fragments:',anevap(4)
1762  ansho=ansho/nhkkh1
1763  angre=angre/nhkkh1
1764  angre2=angre2/nhkkh1
1765  anhea=anhea/nhkkh1
1766  anbla=anbla/nhkkh1
1767  anmin=anmin/nhkkh1
1768  ansp8=ansp8/nhkkh1
1769  ansp5=ansp5/nhkkh1
1770  WRITE(6,'(A,F10.3)')' Number of shower part.:',ansho
1771  WRITE(6,'(A,F10.3)')' Number of grey part.:',angre
1772  WRITE(6,'(A,F10.3)')' Number of grey23 part.:',angre2
1773  WRITE(6,'(A,F10.3)')' Number of heavy part.:',anhea
1774  WRITE(6,'(A,F10.3)')' Number of black part.:',anbla
1775  WRITE(6,'(A,F10.3)')' Number of neg.ch.part.:',anmin
1776  WRITE(6,'(A,F10.3)')' Number of slow 8 prot.:',ansp8
1777  WRITE(6,'(A,F10.3)')' Number of slow 5 prot.:',ansp5
1778  DO 3128 i=1,250
1779  amevap(i)=amevap(i)/nhkkh1
1780  amevpp(i)=amevpp(i)/nhkkh1
1781  WRITE(6,'(F12.4,2E15.5)')xmevap(i),amevap(i),amevpp(i)
1782  amevap(i)=log10(amevap(i)+1.d-9)
1783  amevpp(i)=log10(amevpp(i)+1.d-9)
1784  3128 CONTINUE
1785  WRITE(6,'(A)')' Heavy fragment spectrum'
1786  dit=it/50+1
1787  IF(dit.LT.1.d0)dit=1
1788  CALL plot(xmevap,amevap,it,1,it,0.,dit,-5.,0.1)
1789  dip=ip/50+1
1790  IF(dip.LT.1.d0)dip=1
1791  CALL plot(xmevap,amevpp,ip,1,ip,0.,dip,-5.,0.1)
1792  DO 3119 i=1,50
1793  DO 3119 ii=1,2
1794  pevap(i,ii)=log10(pevap(i,ii)/(nhkkh1*dpeva)+1.d-9)
1795  3119 CONTINUE
1796  ENDIF
1797  WRITE(6,'(A)')' Ev. p and n momentum spectrum'
1798  CALL plot(xpevap,pevap,100,2,50,0.,dpeva,-5.,0.1)
1799 C------------------------------------------------------------------
1800 C------------------------------------------------------------------
1801  RETURN
1802  END
c2_transformation::zero
static float_type zero(float_type)
utility function f(x)=0 useful in axis transforms
Definition: c2_function.hh:1203
z
Double_t z
Definition: plot.C:279
mcihad
function mcihad(MCIND)
Definition: dpm25nulib.f:364
d
Float_t d
Definition: plot.C:237
plomb
subroutine plomb(I, PP, CHAR, XF, ITIF, IJPROJ)
Definition: dpm25hist.f:1015
disres
subroutine disres(IOP, IJPROJ, PPN)
Definition: dpm25hist.f:1007
dispt
subroutine dispt(IOP, NHKKH1, PO)
Definition: dpm25hist.f:1025
a
G4double a
Definition: TRTMaterials.hh:39
nmxhkk
const int nmxhkk
Definition: G4DPMJET2_5Interface.hh:78
edensi
subroutine edensi(I, J)
Definition: dpm25hist.f:1011
distrc
subroutine distrc(IOP, NHKKH1, PO, IGENER)
Definition: dpm25hist.f:525
plot
subroutine plot(X, Y, N, M, MM, XO, DX, YO, DY)
Definition: dpm25nulib.f:176
histog
subroutine histog(I)
Definition: dpm25hist.f:1486
sewew
subroutine sewew(IOP, NHKKH1)
Definition: dpm25nuc4.f:33
distpa
subroutine distpa(IOP)
Definition: dpm25hist.f:1003
testem0.App.init
def init
Definition: testem0.py:56
distco
subroutine distco(IOP, IJPROJ, PPN, IDUMMY)
Definition: dpm25hist.f:999
diseva
subroutine diseva(IOP, NHKKH1, PO, IGENER)
Definition: dpm25hist.f:1491
HepGeom::Transform3D::dy
double dy() const
Definition: Transform3D.h:282
pt
TMarker * pt
Definition: egs.C:22
plombc
subroutine plombc(I, PP, CHAR, XF, ITIF, IJPROJ)
Definition: dpm25hist.f:1020
title
subroutine title(NA, NB, NCA, NCB)
Definition: dpm25nuc7.f:1744
c2_factory::log
static c2_log_p< float_type > & log()
make a *new object
Definition: c2_factory.hh:138
G4Abla::tau
G4double tau(G4double bet, G4double homega, G4double ef, G4double t)
Definition: G4Abla.cc:2586
c2_factory::sqrt
static c2_sqrt_p< float_type > & sqrt()
make a *new object
Definition: c2_factory.hh:142
distr
subroutine distr(IOP, NHKKH1, PO, IGENER)
Definition: dpm25hist.f:9
distrp
subroutine distrp(IOP, NHKKH1, PO)
Definition: dpm25hist.f:994
Generated on Sat Dec 14 2013 14:33:15 for Geant4_10 by   doxygen 1.8.5