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