Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dpm25nuc1g4.f
Go to the documentation of this file.
1  BLOCK DATA blkd41
2  IMPLICIT DOUBLE PRECISION (a-h,o-z)
3  SAVE
4 C
5 *KEEP,PANAME.
6 C------------------
7 C
8 C /PANAME/ CONTAINS PARTICLE NAMES
9 C BTYPE = LITERAL NAME OF THE PARTICLE
10 C
11  CHARACTER*8 btype
12  COMMON /paname/ btype(30)
13 *KEEP,NUCIMP.
14  COMMON /nucimp/ prmom(5,248),tamom(5,248), prmfep,prmfen,tamfep,
15  +tamfen, prefep,prefen,taefep,taefen, prepot(210),taepot(210),
16  +prebin,taebin,fermod,etacou
17 *KEEP,DPRIN.
18  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
19 *KEEP,DROPPT.
20  LOGICAL intpt,fermp,ihadss,ihadsv,ihadvs,ihadvv,ihada, ipadis,
21  +ishmal,lpauli
22  COMMON /droppt/ intpt,fermp,ihadss,ihadsv,ihadvs,ihadvv,ihada,
23  +ipadis,ishmal,lpauli
24 *KEEP,HADTHR.
25  COMMON /hadthr/ ehadth,inthad
26 *KEEP,NSHMAK.
27  COMMON /nshmak/ nnshma,npshma,ntshma,nshmac,nshma2
28 *KEEP,REJEC.
29  COMMON /rejec/ irco1,irco2,irco3,irco4,irco5, irss11,irss12,
30  +irss13,irss14, irsv11,irsv12,irsv13,irsv14, irvs11,irvs12,irvs13,
31  +irvs14, irvv11,irvv12,irvv13,irvv14
32 *KEEP,DSHM.
33  COMMON /dshm/ rash,rbsh,bmax,bstep,sigsh,rosh,gsh,
34  * bsite(0:1,200),nstatb,nsiteb
35 *KEEP,DAMP.
36  COMPLEX*16 ca,ci
37  COMMON /damp/ ca,ci,ga
38 *KEEP,INTMX.
39  parameter(intmx=2488,intmd=252)
40 *KEND.
41  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
42  COMMON /nomije/ ptmije(10),nnmije(10)
43  DATA ptmije /5.d0,7.d0,9.d0,11.d0,13.d0,15.d0,17.d0
44  +,19.d0,21.d0,23.d0 /
45 *
46 *---rejection counters
47  DATA irco1,irco2,irco3,irco4,irco5 /5*0/
48  DATA irss11,irss12,irss13,irss14,irsv11,irsv12,irsv13,irsv14 /8*0/
49  DATA irvs11,irvs12,irvs13,irvs14,irvv11,irvv12,irvv13,irvv14 /8*0/
50 C-------------------
51  DATA inthad /0/
52 *---predefinition of nuclear potentials, average binding energies, ...
53  DATA prepot /210*0.0/
54  DATA taepot /210*0.0/
55  DATA taebin,prebin,fermod /2*0.0d0,0.6d0/
56 *---internal particle names
57  DATA btype /'PROTON ' , 'APROTON ' , 'ELECTRON' , 'POSITRON' ,
58  +'NEUTRIE ' , 'ANEUTRIE' , 'PHOTON ' , 'NEUTRON ' , 'ANEUTRON' ,
59  +'MUON+ ' , 'MUON- ' , 'KAONLONG' , 'PION+ ' , 'PION- ' ,
60  +'KAON+ ' , 'KAON- ' , 'LAMBDA ' , 'ALAMBDA ' , 'KAONSHRT' ,
61  +'SIGMA- ' , 'SIGMA+ ' , 'SIGMAZER' , 'PIZERO ' , 'KAONZERO' ,
62  +'AKAONZER' , 'RESERVED' , 'BLANK ' , 'BLANK ' , 'BLANK ' ,
63  +'BLANK ' /
64 *---print options
65  DATA ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr /0, 0, 0, -1, 0,
66  +0, 0, 0/
67 C-------------------
68  DATA intpt, fermp, ihadss,ihadsv,ihadvs,ihadvv, ihada /.true.,
69  +.true., 4*.false., .true./
70  DATA ipadis, ishmal, lpauli /.false., .false., .true./
71 C----------------------------------
72  DATA nshmac /0/
73  DATA nshma2 /0/
74 *---parameters for Glauber initialization / calculation
75  DATA nstatb, nsiteb /2000, 200/
76  DATA ci /(1.0,0.0)/
77 *---parameters for combination of q-aq chains to color ropes
78 C DATA LCOMBI /.FALSE./, NCUTOF /100/
79  DATA isingd,idiftp,ioudif,iflagd /0,0,0,0/
80 *
81  END
82 ************************************************************************
83 ************************************************************************
84 *
85  SUBROUTINE dttest(CODEWD,WHAT,SDUM)
86 *
87 * -- not for normal user --
88 * contains input options unrecognized by DTPREP and
89 * performs special initialisations or tasks for program devoloping
90 *
91 C COMMON fully commented in DTPREP
92 *
93 * **********************************************************************
94 * * DESCRIPTION OF THE COMMON BLOCK(S), VARIABLE(S) AND DECLARATIONS *
95 * **********************************************************************
96 *
97 *
98 * /USER/ contains the parameters, expected to be modified by normal user
99 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100  IMPLICIT DOUBLE PRECISION(a-h,o-z)
101  SAVE
102  CHARACTER*80 title
103  CHARACTER*8 projty,targty
104 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
105 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
106  COMMON /user1/title,projty,targty
107  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
108 *
109 * /COLLE/ contains the input specifying the MC. run
110 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111  COMMON /colle/nevhad,nvers,ihadrz,nfile
112 *
113 * /COLLIS/ contains the input specifying the considered event
114 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115  common/collis/s,ijproj,ijtar,ptthr,ptthr2,iophrd,ijprlu,ijtalu
116 *
117 *
118 * /BOOKLT/ contains the final particle names and PPDB-numbers
119 * of 30 final particle types
120 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121  CHARACTER*8 btype
122  common/booklt/btype(30),nbook(30)
123 *
124 * /POLMN/ stores arrays describing probabilities of parton
125 * configurations
126 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
128  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
129 *
130 * /POMTYP/ contains parameters determining X-sections
131 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
133 *
134 * various smaller commons
135 * in alphabetical order
136 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137  COMMON /dropjj/dropjt,dropva
138  COMMON /gluspl/nugluu,nsgluu
139  COMMON /ptlarg/xsmax
140  COMMON /ptsamp/ isampt
141  COMMON /stars/istar2,istar3
142  COMMON /strufu/istrum,istrut
143  COMMON /popcor/pdb,ajsdef
144 *
145 * ********************************************************************
146 * declarations outside of commons
147 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148  CHARACTER*8 codewd,sdum
149  dimension what(6)
150 * ********************************************************************
151 *
152 *
153 * **********************************************************************
154 * * Print the warning that no normal codeword *
155 * **********************************************************************
156 *
157  WRITE(6,9)
158  9 FORMAT( ' special code word was used ')
159 *
160 *
161 * The following additional CODEWD options exist at the moment:
162 * RANDOMIZ SIGMAPOM PARTEV SELHARD
163 * GLUSPLIT
164 *
165 * The cards marked with )+ have to be followed by data cards of
166 * special format
167 *
168 *
169 *
170 * **********************************************************************
171 * * parse stored imput card *
172 * **********************************************************************
173 *
174 *
175 *
176 * *********************************************************************
177 * input card: CODEWD = RANDOMIZE
178 * Sets the SEED for the random number generators
179 *
180 * WHAT(1) = 1: gets testrun otherwise: reset
181 * WHAT(2,3,4,5) = giving the SEED for the random number generators.
182 * Default ISEED1/2/3/4=12,34,56,78
183 * Note. It is advisable to use only the seeds given by the
184 * program in earlier runs. Otherwise the number sequence might
185 * have defects in its randomness.
186 C Since 1999 RANDOMIZE initializes the RM48 generator
187 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188 *
189 C IF(CODEWD.EQ.'RANDOMIZ') THEN
190 *
191 * Reinitialize random generator
192 C IF(WHAT(2).NE.0.D0) THEN
193 C ISEED1=WHAT(2)
194 C ISEED2=WHAT(3)
195 C ISEED3=WHAT(4)
196 C ISEED4=WHAT(5)
197 C CALL RNDMST(ISEED1,ISEED2,ISEED3,ISEED4)
198 C ENDIF
199 * Test random generator and test
200 C IF(WHAT(1).EQ.1) CALL RNDMTE(1)
201 *
202 C GO TO 1
203 *
204 *
205 * *********************************************************************
206 * input card: CODEWD = SIGMAPOM
207 * Defines the options for the calculation of the u
208 * tarized hard and soft multi-pomeron cross sectio
209 * and demands a testrun
210 *
211 * WHAT(1) = ITEST testrun for ITEST = 1
212 * WHAT(2) = ISIG characterizing X sections, see SIGSHD
213 * default: 10
214 * Only ISIG=10 kept in dpmjet-II.5
215 *
216 * WHAT(3) = characterizes the method to calculate
217 * the cut pomeron X-section SIGMA(Lsoft,Mhard,Ntrp) and
218 * how to attribute them to strings
219 C!!!!!!!!!!!!!!!!!!!
220 C!!!!!!!!!!!!!!!!!!! THE FOLLOWING DISTRIBUTIONS WITH 2 CHANNEL
221 C!!!!!!!!!!!!!!!!!!! EIKONAL + H. MASS DIFFRACTION + HARD SCATTERING
222 C!!!!!!!!!!!!!!!!!!!
223 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
224 C
225 C =482: SEE PRBLM2
226 C
227 * DEFAULT: 482
228 * for all cases
229 * WHAT(4) = LMAX maximal considered number soft Pomerons
230 * WHAT(5) = MMAX maximal considered number hard Pomerons
231 * WHAT(6) = NMAX maximal considered number trippel Pomerons
232 *
233 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234 *
235 C ELSEIF(CODEWD.EQ.'SIGMAPOM') THEN
236  IF(codewd.EQ.'SIGMAPOM') THEN
237 *
238  itest = what(1)
239  IF(what(2).NE.0.) isig =int(what(2))
240  IF(what(3).NE.0.) ipim =int(what(3))
241  IF(ipim.GT.10) THEN
242  icon = ipim /10
243  ipim = ipim-10*icon
244  ENDIF
245  IF(ipim.EQ.1) THEN
246  difel = what(4)
247  difnu = what(5)
248  lmax =int(what(6))
249  mmax = lmax
250  ELSE
251  lmax =int(what(4))
252  mmax =int(what(5))
253  nmax =int(what(6))
254  ENDIF
255  IF (itest.EQ.1)CALL pomdi
256  go to 1
257 *
258 * ********************************************************************
259 * input card: CODEWD = GLUSPLIT
260 * Prevents the splitting of Gluons
261 *
262 * WHAT(1)=NUGLUU Default: 1.
263 * =1. only one jet in hard gluon scattering
264 * WHAT(2)=NSGLUU Default: 0.
265 * =0. two jets in soft sea gluons
266 * =1. only one jet in soft sea gluons
267 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268 *
269  ELSEIF(codewd.EQ.'GLUSPLIT') THEN
270 *
271  nugluu = what(1)
272  nsgluu = what(2)
273  go to 1
274 *
275 * ********************************************************************
276 *
277 * *********************************************************************
278 * input card: CODEWD = PARTEV
279 * defines the parton level collision events
280 * the X's PT's and flavors and
281 * demands a test run
282 *
283 * WHAT(1) = 1.: testrun ; other values: no testrun
284 * WHAT(2) = number of events is NPEV default: 30
285 * WHAT(3) = version of PARTEV is NVERS default: 1
286 * NVERS=1 all hard partons considered to be gluons
287 * SOFT X-VALUEA BY REJECTION
288 * NVERS=2 all hard partons considered to be gluons
289 * SOFT X-VALUEA BY AURENCHE -MAIRE METHOD
290 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291 *
292  ELSEIF(codewd.EQ.'PARTEV ') THEN
293 *
294  itest = int(what(1))
295  IF (what(2).EQ.0.d0)npev=30
296  IF (what(2).NE.0.d0)npev=int(what(2))
297  IF (what(3).EQ.0.d0)nvers=1
298  IF (what(3).NE.0.d0)nvers=int(what(3))
299 * first initialize NSOFT-NHARD event selection
300 * corresponing to choosen one options
301  IF(itest.EQ.1) THEN
302  IF(ipim.EQ.2)CALL prblm2(cmener)
303 * initialize hard scattering
304  CALL jtdtu(0)
305  CALL samppt(0,pt)
306 C CALL PARTEV(NPEV)
307  CALL timdat
308  ENDIF
309  go to 1
310 *
311 * *********************************************************************
312 * input card: CODEWD = SELHARD
313 * defines the selection of X's,PT's and
314 * flavors for hard scattering
315 *
316 * WHAT(2) = IOPHRD selects the model default: 2
317 * WHAT(4) = DROPJT=10:DROP FIRST HARD JET PAIR DEFAULT: 0
318 * THIS OPTION SHOULD ALLOW TO SIMULATE DRELL-YAN
319 * OR W OR Z PRODUCTION EVENTS
320 * WHAT(5) = PTTHR threshold PT for hard scattering default: 3.
321 * WHAT(6) = PTTHR2 threshold PT for FIRST HARD scattering default: 3.
322 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
323 *
324  ELSEIF(codewd.EQ.'SELHARD ') THEN
325 *
326  IF(what(2).NE.0.d0)iophrd=int(what(2))
327  IF(what(4).NE.0.d0)dropjt=what(4)
328  IF(what(6).NE.0.d0)ptthr2=what(6)
329  IF(what(5).NE.0.d0)THEN
330  ptthr=what(5)
331  IF(cmener.LT.2000.0d0.AND.isig.EQ.3)ptthr=what(5)
332  IF (cmener.GE.2000.0d0.AND.isig.EQ.3)
333  * ptthr=0.25*log(cmener/2000.)+2.
334  IF(ptthr2.LT.ptthr)ptthr2=ptthr
335  IF(istrut.EQ.1)THEN
336  ptthr=2.1+0.15*(log10(cmener/50.))**3
337  ptthr2=ptthr
338  ELSEIF(istrut.EQ.2)THEN
339  ptthr=2.5+0.12*(log10(cmener/50.))**3
340  ptthr2=ptthr
341  ENDIF
342  WRITE(6,1244)ptthr
343  1244 FORMAT (' THRESHOLD PT FOR HARD SCATTERING PTTHR=',f12.2)
344  ENDIF
345  go to 1
346 *
347 * *********************************************************************
348 * input card: CODEWD = XSLAPT
349 * calculates inclusive large PT cross sections
350 * and testrun of the large pt and minijet sampling
351 * in file laptabr
352 *
353 * WHAT(I) = not used at this level
354 * special DATA CARDS are required
355 * see description at beginning of LAPTABR
356 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ----
357 *
358  ELSEIF(codewd.EQ.'XSLAPT ') THEN
359 *
360  CALL timdat
361  CALL laptab
362  CALL timdat
363  go to 1
364 *
365 * *********************************************************************
366 * parameter card: CODEWD = SAMPT
367 * defines the options of soft pt sampling in
368 * subroutine SAMPPT
369 * WHAT(1) = number defining the option of soft pt sampling
370 * (default: 4 )
371 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ----
372 *
373  ELSEIF(codewd.EQ.'SAMPT ') THEN
374 *
375  isampt = int( what(1) )
376  IF( isampt.LT.0 .OR. isampt.GT.4 ) isampt=0
377  go to 1
378 *
379 * *********************************************************************
380 * ending special parsing the code word of "input card"
381 *
382  ENDIF
383 *
384 * 1) not recognized cards
385 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ----
386 * a warning will be issued in DTPREP
387  RETURN
388 *
389 * 2) recognized cards
390 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ----
391 * action was done,
392 * CODEWD is set to a value ignored in DTPREP
393 *
394  1 codewd='-zzzzzzz'
395  RETURN
396  END
397 *
398 ************************************************************************
399 ************************************************************************
400 *
401  BLOCK DATA bookle
402 *
403 * *********************************************************************
404 * /BOOKLT/
405 *
406 * neeeded in the following routines: DTUMAIN
407 *
408 * description
409 * /BOOKLT/ contains the final particle names and PPDB-numbers
410 * BTYPE = literal name of the particle
411 * NBOOK = the number of the particle
412 * proposed in the particle data booklet (90)
413 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414  IMPLICIT DOUBLE PRECISION(a-h,o-z)
415  SAVE
416  CHARACTER*8 btype
417  common/booklt/btype(30),nbook(30)
418 *
419  DATA btype /'PROTON ' , 'APROTON ' , 'ELECTRON' ,
420  1 'POSITRON' , 'NEUTRIE ' , 'ANEUTRIE' ,
421  2 'PHOTON ' , 'NEUTRON ' , 'ANEUTRON' ,
422  3 'MUON+ ' , 'MUON- ' , 'KAONLONG' ,
423  4 'PION+ ' , 'PION- ' , 'KAON+ ' ,
424  5 'KAON- ' , 'LAMBDA ' , 'ALAMBDA ' ,
425  6 'KAONSHRT' , 'SIGMA- ' , 'SIGMA+ ' ,
426  7 'SIGMAZER' , 'PIZERO ' , 'KAONZERO' ,
427  9 'AKAONZER' , ' ' , ' ' ,
428  z ' ' , ' ' , ' ' /
429 *
430 *
431  DATA nbook / 2212 , -2212 , 11 ,
432  1 -11 , 14 , -14 ,
433  2 22 , 2112 , -2112 ,
434  3 -13 , 13 , 130 ,
435  4 211 , -211 , 321 ,
436  5 -321 , 3122 , -3122 ,
437  6 310 , 3114 , 3224 ,
438  7 3214 , 111 , 311 ,
439  9 -311 , 0 , 0 ,
440  z 0 , 0 , 0 /
441 *
442  END
443 
444 
445 C______________________________________________________________________
446  SUBROUTINE samppt(MODE,PT)
447 * pt for partons at the end of soft chains
448 * this pt is sampled from the distribution
449 * exp(-b*pt^2) with pt=0..ptcut
450 * MODE = 0 - initialization to determine parameter b from total soft
451 * and differential hard cross section
452 * MODE = 1 - sample pt
453 * MODE = 2 - PLOT PT DISTRIBUTION SAMPLED
454 *
455  IMPLICIT DOUBLE PRECISION(a-h,o-z)
456  SAVE
457  parameter( zero=0.d0, one=1.d0)
458  parameter( alfa=0.56268d-01, beta=0.17173d+03 )
459  parameter( acc = 0.0001d0 )
460  COMMON /xsecpt/ ptcut,sigs,dsigh
461  COMMON /sigma / sigsof,bs,zsof,sighar,fill(7)
462  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
463 C repl. COMMON/COLLIS/ECM,S,IJPROJ,IJTAR,PTTHR,IOPHRD,IJ1LU,IJ2LU,PTTHR2
464  common/collis/s,ijproj,ijtar,ptthr,ptthr2,iophrd,ijprlu,ijtalu
465  CHARACTER*80 title
466  CHARACTER*8 projty,targty
467 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
468 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
469  COMMON /user1/title,projty,targty
470  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
471 *
472  common/ptsamp/ isampt
473  dimension pptt(50),dpptt(50)
474  DATA ecm0 /0.1d0/
475 C to keep identical commons for patchy ect
476  ecm=cmener
477  ptthr=2.5+0.12*(log10(cmener/50.))**3
478  ptcut=ptthr
479  CALL sigshd(ecm)
480  IF ( mode.EQ.0 ) THEN
481  DO 201 ii=1,50
482  pptt(ii)=ii*ptcut/50.
483  dpptt(ii)=0.
484  201 CONTINUE
485  sigs = 0.15*sigsof
486  IF(ecm.LT.1000.)THEN
487  aacucu=0.85*(ecm-400.)/600.
488  sigs=(1.-aacucu)*sigsof
489  ENDIF
490 C*************************************************************
491 C
492 C OPTIONS FOR SOFT PT SAMPLING
493 C
494 CWRITE(6,'(A,4E12.4)')' SAMPPT:ECM,PTCUT,SIGS,DSIGH',
495 C * ECM,PTCUT,SIGS,DSIGH
496  IF(ecm0.NE.ecm)THEN
497 C WRITE(6,'(A,5E12.4)')' SAMPPT:ECM,PTCUT,SIGS,DSIGH,SIGHAR',
498 C * ECM,PTCUT,SIGS,DSIGH,SIGHAR
499 C WRITE(6,5559)PTCUT,SIGSOF,SIGHAR,ISAMPT
500  5559 FORMAT(' SAMPPT:PTCUT,SIDSOF.SIGHATD,ISAMPT:',3e12.3,i5)
501  ENDIF
502  IF( isampt.EQ.0 ) THEN
503  c = dsigh/(2.*sigs*ptcut)
504  b = bsofpt(acc,c,ptcut)
505  IF( ioutpa.GE.1 )WRITE(6,9010)mode,isampt,ptcut,sigs,dsigh,b
506  * ,c,sigsof,sighar,rmin
507  ELSEIF( isampt.EQ.1 ) THEN
508  eb = ecm/beta
509  c = alfa*log(eb)
510  b = bsofpt(acc,c,ptcut)
511  IF( ioutpa.GE.1 )WRITE(6,9010)mode,isampt,ptcut,sigs,dsigh,b
512  * ,c,sigsof,sighar,rmin
513  ELSEIF( isampt.EQ.2 ) THEN
514  b=-6.
515  IF( ioutpa.GE.1 )WRITE(6,9010)mode,isampt,ptcut,sigs,dsigh,b
516  * ,c,sigsof,sighar,rmin
517  ELSEIF( isampt.EQ.3 ) THEN
518  b=1.e-06
519  IF( ioutpa.GE.1 )WRITE(6,9010)mode,isampt,ptcut,sigs,dsigh,b
520  * ,c,sigsof,sighar,rmin
521  ELSEIF( isampt.EQ.4)THEN
522  aaaa=ptcut**2*(sigsof+sighar)
523  IF (aaaa.LE.0.00001d0) THEN
524  aaaa=abs(aaaa)+0.0002
525 C WRITE(6,5559)PTCUT,SIGSOF,SIGHAR
526 C5559 FORMAT(' SAMPPT:PTCUT,SIDSOF.SIGHATD:',3E12.3)
527  ENDIF
528  c=sighar/aaaa
529  b = 0.5*bsofpt(acc,c,ptcut)
530  IF( ioutpa.GE.1 )WRITE(6,9010)mode,isampt,ptcut,sigs,dsigh,b
531  * ,c,sigsof,sighar,rmin
532  ENDIF
533  ecm0=ecm
534 C*************************************************************
535  rmin = exp(b*ptcut**2)
536 C
537 C IOUTPO=IOUTPA
538 C IOUTPA=1
539  IF( ioutpa.GE.1 )WRITE(6,9010)mode,isampt,ptcut,sigs,dsigh,b
540  * ,c,sigsof,sighar,rmin
541 9010 FORMAT(' SAMPPT MODE,ISAMPT,PTCUT,SIGS,DSIGH,B,C,SIGSOF',
542  *' SIGHAR,RMIN ',
543  * 2i2,f5.2,7e13.6)
544 C IOUTPA=IOUTPO
545  ELSEIF ( mode.EQ.1 ) THEN
546  IF( ioutpa.GE.1 )WRITE(6,9010)mode,isampt,ptcut,sigs,dsigh,b
547  * ,c,sigsof,sighar,rmin
548  ptt =log(1.0-rndm(v)*(1.0-rmin))/(b+0.00001d0)
549  pt=sqrt(ptt)
550  iipt=pt*50./ptcut+1.
551  iipt=min( iipt,50 )
552  dpptt(iipt)=dpptt(iipt)+1./(pt+0.000001d0)
553 C WRITE(6,111)MODE,PTT,PT,B,RMIN
554 C 111 FORMAT ('SAMPPT: MODE,PTT,PT,B RMIN',I5,4E15.8)
555  ELSEIF(mode.EQ.2)THEN
556  DO 202 ii=1,50
557  dpptt(ii)=log10(1.e-8+dpptt(ii))
558  202 CONTINUE
559  IF(iouxev.GE.-1)THEN
560  WRITE (6,203)
561  203 FORMAT(' PT DISTRIBUTION OF SOFT PARTONS AS SAMPLED IN BSOFPT')
562  CALL plot(pptt,dpptt,50,1,50,zero,ptcut/50.d0,zero,0.05d0*one)
563  ENDIF
564  ENDIF
565  RETURN
566  END
567 C****************************************************************8****
568  real
569  * * 8
570  * FUNCTION bsofpt(ACC,CC,PPTCUT)
571  IMPLICIT DOUBLE PRECISION(a-h,o-z)
572  SAVE
573  LOGICAL succes
574  COMMON /bsoff1/c,ptcut
575  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
576  EXTERNAL bsofc1,bsof1
577  dimension x(50),y(50)
578  c=cc
579  ptcut=pptcut
580  DO 100 i=1,50
581  x(i)=-2.5+i*0.1
582  y(i)=bsof1(x(i))
583  100 CONTINUE
584 C CALL PLOT (X,Y,50,1,50,-2.5D0,0.1D0,-1.D0,0.02D0)
585  IF(c.LT.1.d-10) THEN
586  bb=-30.
587  go to 999
588  ENDIF
589 C IF (C.GT.1.) THEN
590  kkkk=0
591  jjjj=0
592  b1=c+3.
593  b2=0.0001
594  go to 300
595  400 CONTINUE
596  kkkk=kkkk+1
597 C ENDIF
598 C IF (C.LT.1.)THEN
599  b1=-0.00001
600  b2=-3.
601 C ENDIF
602  300 CONTINUE
603  CALL zbrac(bsof1,b1,b2,succes)
604  IF (.NOT.succes)THEN
605  IF (kkkk.EQ.0)go to 400
606  jjjj=1
607  ENDIF
608  IF(iouxev.GE.1)WRITE(6,111)b1,b2
609  111 FORMAT(2f10.4)
610  IF (succes)THEN
611  bb=rtsafe(bsofc1,b1,b2,acc)
612  ENDIF
613  IF (jjjj.EQ.1)bb=0.
614  999 CONTINUE
615  bsofpt=bb
616  RETURN
617  END
618  SUBROUTINE bsofc1(B,F,DF)
619  IMPLICIT DOUBLE PRECISION(a-h,o-z)
620  SAVE
621  COMMON /bsoff1/c,ptcut
622  aaa=exp(b*ptcut**2)
623  f=c*(aaa-1.)-b*aaa
624  df=c*ptcut**2*aaa-aaa
625  * -b*ptcut**2*aaa
626  RETURN
627  END
628 *
629 *******************************************************************
630 *
631  real
632  * * 8
633  * FUNCTION bsof1(B)
634  IMPLICIT DOUBLE PRECISION(a-h,o-z)
635  SAVE
636  COMMON /bsoff1/c,ptcut
637  qqq=b*ptcut**2
638  aaa=0.
639  IF(qqq.GT.-60.) THEN
640  aaa=exp(b*ptcut**2)
641  ENDIF
642  bsof1=c*(aaa-1.)-b*aaa
643 C WRITE(6,10)B,PTCUT,BSOF1,C
644 C 10 FORMAT (4E15.4)
645  RETURN
646  END
647 *
648 *******************************************************************************
649  real
650  * * 8
651  * FUNCTION rtsafe(FUNCD,X1,X2,XACC)
652  IMPLICIT DOUBLE PRECISION(a-h,o-z)
653  SAVE
654  parameter(maxit=200,itepri=0)
655  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
656  CALL funcd(x1,fl,df)
657 C IF(IOUXEV.GE.0.AND.ITEPRI.EQ.1) WRITE(6,9999) FL,DF
658  CALL funcd(x2,fh,df)
659 C IF(IOUXEV.GE.0.AND.ITEPRI.EQ.1) WRITE(6,9999) FH,DF
660  IF(fl*fh.GE.0.) pause 'ROOT MUST BE BRACKETED'
661  IF(fl.LT.0.)THEN
662  xl=x1
663  xh=x2
664  ELSE
665  xh=x1
666  xl=x2
667  swap=fl
668  fl=fh
669  fh=swap
670  ENDIF
671  rtsafe=.5*(x1+x2)
672  dxold=abs(x2-x1)
673  dx=dxold
674  CALL funcd(rtsafe,f,df)
675 C IF(IOUXEV.GE.0.AND.ITEPRI.EQ.1) WRITE(6,9998) RTSAFE,F,DF
676 C IF(IOUXEV.GE.0.AND.ITEPRI.EQ.1) WRITE(6,9996)
677  DO 11 j=1,maxit
678 C IF(IOUXEV.GE.0.AND.ITEPRI.EQ.1)
679 C * WRITE(6,9997) RTSAFE,XH,XL,DXOLD,F,DF
680  vr1 = var( rtsafe,xh,df,f )
681  vr2 = var( rtsafe,xl,df,f )
682 C IF(IOUXEV.GE.0.AND.ITEPRI.EQ.1) WRITE(6,9995) VR1,VR2
683 C IF(((RTSAFE-XH)*DF-F)*((RTSAFE-XL)*DF-F).GE.0.
684 C * .OR. ABS(2.*F).GT.ABS(DXOLD*DF) ) THEN
685  IF( vr1*vr2 .GE. 0.
686  * .OR. abs(2.*f).GT.abs(dxold*df) ) THEN
687  dxold=dx
688  dx=0.5*(xh-xl)
689  rtsafe=xl+dx
690  IF(xl.EQ.rtsafe)RETURN
691  ELSE
692  dxold=dx
693 C IF(IOUXEV.GE.0.AND.ITEPRI.EQ.1) WRITE(6,9999) F,DF
694  dx=f/df
695  temp=rtsafe
696  rtsafe=rtsafe-dx
697  IF(temp.EQ.rtsafe)RETURN
698  ENDIF
699  IF(abs(dx).LT.xacc) RETURN
700  CALL funcd(rtsafe,f,df)
701  IF(f.LT.0.) THEN
702  xl=rtsafe
703  fl=f
704  ELSE
705  xh=rtsafe
706  fh=f
707  ENDIF
708 11 CONTINUE
709  pause 'RTSAFE EXCEEDING MAXIMUM ITERATIONS'
710  RETURN
711 9995 FORMAT(' VR1,VR2:',2e12.5)
712 9996 FORMAT(' RTSAFE,XH,XL,DXOLD,F,DF IN LOOP 11 J=1,MAXIT')
713 9997 FORMAT(3x,6e10.3)
714 9998 FORMAT(' RTSAFE: RTSAFE,F,DF =',3e12.5)
715 9999 FORMAT(' RTSAFE: F,DF =',2e12.5)
716  END
717 *
718 *****************************************************************
719  real
720  * * 8
721  * FUNCTION var(A,B,C,D)
722  IMPLICIT DOUBLE PRECISION(a-h,o-z)
723  SAVE
724  parameter( ambmax = 1.0d+38, epsi = 1.2d-38, one=1.d0 )
725  amb = a - b
726  siab= sign(one,amb)
727  abl = abs(amb)
728  abl = log10( abl + epsi )
729  sicc= sign(one, c )
730  ccl = abs( c )
731  ccl = log10( ccl + epsi )
732  rcheck=abl + ccl
733  IF( rcheck .LE. 38.d0 ) THEN
734  var = amb*c-d
735  ELSE
736  var = ambmax*siab*sicc - d
737  ENDIF
738  IF( var .GT. 1.0d+18 ) var = 1.0e+18
739  IF( var .LT. -1.0d+18 ) var = -1.0e+18
740  RETURN
741  END
742 C
743  SUBROUTINE zbrac(FUNC,X1,X2,SUCCES)
744  IMPLICIT DOUBLE PRECISION(a-h,o-z)
745  SAVE
746  EXTERNAL func
747  parameter(factor=1.6d0,ntry=50)
748  LOGICAL succes
749  IF(x1.EQ.x2)pause 'You have to guess an initial range'
750  f1=func(x1)
751  f2=func(x2)
752  succes=.true.
753  DO 11 j=1,ntry
754  IF(f1*f2.LT.0.d0)RETURN
755  IF(abs(f1).LT.abs(f2))THEN
756  x1=x1+factor*(x1-x2)
757  f1=func(x1)
758  ELSE
759  x2=x2+factor*(x2-x1)
760  f2=func(x2)
761  ENDIF
762 11 CONTINUE
763  succes=.false.
764  RETURN
765  END
766 *
767 C
768 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
769 C
770 C
771 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
772 C
773 c subroutine hibldr
774 c COMMON /DREAC/ umo(296),plabf(296),siin(296),wk(5184),
775 c * nrk(2,268),nure(30,2)
776 c read(2,1)umo,plabf,siin,wk
777 c 1 format (5e16.7)
778 c read(2,2)nrk,nure
779 c 2 format (8i10)
780 c return
781 c end
782 C
783 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
784 
785 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
786 C
787  SUBROUTINE checkf(EPN,PPN,IREJ,IORIG)
788  IMPLICIT DOUBLE PRECISION (a-h,o-z)
789  SAVE
790 *KEEP,HKKEVT.
791 c INCLUDE (HKKEVT)
792  parameter(amuamu=0.93149432d0)
793  parameter(nmxhkk= 89998)
794 c PARAMETER (NMXHKK=25000)
795  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
796  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
797  +(4,nmxhkk)
798 C
799 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
800 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
801 C THE POSITIONS OF THE PROJECTILE NUCLEONS
802 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
803 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
804 C COMPLETELY CONSISTENT. THE TIMES IN THE
805 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
806 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
807 C
808 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
809 C
810 C NMXHKK: maximum numbers of entries (partons/particles) that can be
811 C stored in the commonblock.
812 C
813 C NHKK: the actual number of entries stored in current event. These are
814 C found in the first NHKK positions of the respective arrays below.
815 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
816 C entry.
817 C
818 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
819 C = 0 : null entry.
820 C = 1 : an existing entry, which has not decayed or fragmented.
821 C This is the main class of entries which represents the
822 C "final state" given by the generator.
823 C = 2 : an entry which has decayed or fragmented and therefore
824 C is not appearing in the final state, but is retained for
825 C event history information.
826 C = 3 : a documentation line, defined separately from the event
827 C history. (incoming reacting
828 C particles, etc.)
829 C = 4 - 10 : undefined, but reserved for future standards.
830 C = 11 - 20 : at the disposal of each model builder for constructs
831 C specific to his program, but equivalent to a null line in the
832 C context of any other program. One example is the cone defining
833 C vector of HERWIG, another cluster or event axes of the JETSET
834 C analysis routines.
835 C = 21 - : at the disposal of users, in particular for event tracking
836 C in the detector.
837 C
838 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
839 C standard.
840 C
841 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
842 C The value is 0 for initial entries.
843 C
844 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
845 C one mother exist, in which case the value 0 is used. In cluster
846 C fragmentation models, the two mothers would correspond to the q
847 C and qbar which join to form a cluster. In string fragmentation,
848 C the two mothers of a particle produced in the fragmentation would
849 C be the two endpoints of the string (with the range in between
850 C implied).
851 C
852 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
853 C entry has not decayed, this is 0.
854 C
855 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
856 C entry has not decayed, this is 0. It is assumed that the daughters
857 C of a particle (or cluster or string) are stored sequentially, so
858 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
859 C daughters. Even in cases where only one daughter is defined (e.g.
860 C K0 -> K0S) both values should be defined, to make for a uniform
861 C approach in terms of loop constructions.
862 C
863 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
864 C
865 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
866 C
867 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
868 C
869 C PHKK(4,IHKK) : energy, in GeV.
870 C
871 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
872 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
873 C
874 C VHKK(1,IHKK) : production vertex x position, in mm.
875 C
876 C VHKK(2,IHKK) : production vertex y position, in mm.
877 C
878 C VHKK(3,IHKK) : production vertex z position, in mm.
879 C
880 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
881 C********************************************************************
882 *KEEP,DELP.
883  COMMON /delp/ delpx,delpy,delpz,delpe
884 *KEEP,TANUIN.
885  COMMON /tanuin/ tasuma,tasubi,tabi,tamasu,tama,taimma
886 *KEEP,NUCC.
887  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
888 *KEEP,DPRIN.
889  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
890 *KEND.
891  DATA icheck/0/
892  help=log10(epn)
893  phelp=0.
894  IF(help.GT.5.d0)phelp=help-5.
895  pthelp=12.+phelp*5.
896  irej=0
897  irejj=0
898  px=0.
899  py=0.
900  pz=0.
901  pe=0.
902  eext=0.
903  eexp=0.
904  eee1=0.
905  eeem1=0.
906  ee1001=0.
907  DO 10 i=1,nhkk
908 C Projectiles
909  IF(isthkk(i).EQ.11.OR.isthkk(i).EQ.13) THEN
910  gam=epn/phkk(5,i)
911  bgam=ppn/phkk(5,i)
912  px=px - phkk(1,i)
913  py=py - phkk(2,i)
914  pz=pz - gam*phkk(3,i) - bgam*phkk(4,i)
915  pe=pe - gam*phkk(4,i) - bgam*phkk(3,i)
916  ENDIF
917 C Target
918  IF(isthkk(i).EQ.12.OR.isthkk(i).EQ.14) THEN
919  px=px - phkk(1,i)
920  py=py - phkk(2,i)
921  pz=pz - phkk(3,i)
922  pe=pe - phkk(4,i)
923  ENDIF
924 * sum final state momenta
925  IF(isthkk(i).EQ.1) THEN
926  px=px + phkk(1,i)
927  py=py + phkk(2,i)
928  pz=pz + phkk(3,i)
929  pe=pe + phkk(4,i)
930  ENDIF
931 C noninteracting Projectiles
932  IF(isthkk(i).EQ.13.AND.jdahkk(2,i).EQ.0) THEN
933  gam=epn/phkk(5,i)
934  bgam=ppn/phkk(5,i)
935  px=px + phkk(1,i)
936  py=py + phkk(2,i)
937  pz=pz + gam*phkk(3,i) + bgam*phkk(4,i)
938  pe=pe + gam*phkk(4,i) + bgam*phkk(3,i)
939  ENDIF
940  IF(isthkk(i).EQ.14.AND.jdahkk(2,i).EQ.0) THEN
941 C noninteracting Targets
942  px=px + phkk(1,i)
943  py=py + phkk(2,i)
944  pz=pz + phkk(3,i)
945  pe=pe + phkk(4,i)
946  ENDIF
947  IF(isthkk(i).EQ.16) THEN
948  imo=jmohkk(1,i)
949  px=px + phkk(1,i)
950  py=py + phkk(2,i)
951  pz=pz + phkk(3,i)
952  pe=pe + phkk(4,i)
953  eext=eext + phkk(4,i) - phkk(4,imo)
954  ENDIF
955  IF(isthkk(i).EQ.15) THEN
956  imo=jmohkk(1,i)
957  gam=epn/phkk(5,i)
958  bgam=ppn/phkk(5,i)
959  px=px + phkk(1,i)
960  py=py + phkk(2,i)
961  pz=pz + gam*phkk(3,i) + bgam*phkk(4,i)
962  pe=pe + gam*phkk(4,i) + bgam*phkk(3,i)
963  eext=eext + phkk(4,i) - phkk(4,imo)
964  ENDIF
965  IF(isthkk(i).EQ.1) THEN
966  eee1=eee1+phkk(4,i)
967  ENDIF
968  IF(isthkk(i).EQ.-1) THEN
969  eeem1=eeem1+phkk(4,i)
970  ENDIF
971  IF(isthkk(i).EQ.1001) THEN
972  ee1001=ee1001+phkk(4,i)
973  ENDIF
974  10 CONTINUE
975  eee=eee1+eeem1+ee1001
976  epnto=epn*ip
977  aeee=eee/epnto
978  eeee=eee-epn*ip
979  aeeee=eee/epn-ip
980  aeee1=eee1/epnto
981  aeeem1=eeem1/epnto
982  aee101=ee1001/epnto
983  aip=1
984  ait=it
985  aitz=itz
986  aip=aip+(ait*amuamu+1.d-3*energy(ait,aitz))/epnto
987  delle=abs(aip-aeee)
988  elle=delle*epnto
989  tole=0.030
990 C TOLE=0.012
991  IF(it.LE.50)THEN
992  IF(it.EQ.ip)tole=0.02
993 C IF(IT.EQ.IP)TOLE=0.05
994  ENDIF
995  IF(delle.GE.tole)irej=1
996  IF(irej.EQ.1)THEN
997  icheck=icheck+1
998  IF(icheck.LE.100)THEN
999  WRITE(6,'(A,I5,E10.3,5F10.4)')
1000  * ' IP,EPN,AEEE,AEEEE,AEEE1,AEEEM1,AEE101:',
1001  * ip,epn,aeee,aeeee,aeee1,aeeem1,aee101
1002  WRITE(6,'(A,I5,E10.3,7E12.4)')
1003  * ' IP,EPN,EEE,EEEE,EEE1,EEEM1,EE1001,DELLE,ELLE:',
1004  * ip,epn,eee,eeee,eee1,eeem1,ee1001,delle,elle
1005  ENDIF
1006  ENDIF
1007 C PX=PX + DELPX
1008 C PY=PY + DELPY
1009 C PZ=PZ + DELPZ
1010 C PE=PE + DELPE
1011 C IF(IPRI.GT.1) THEN
1012 C IF (ABS(PX).GT.PTHELP.OR. ABS(PY).GT.PTHELP.OR.
1013 C * ABS(PZ)/EPN.GT.0.025*IP.
1014 C + OR. ABS(PE)/EPN.GT.0.025*IP) THEN
1015 C IREJJ=1
1016 C ICHECK=ICHECK+1
1017 C IF(ICHECK.LE.500.AND.IREJJ.EQ.1)THEN
1018 C WRITE(6,1000) PX,PY,PZ,PE,EEXT,EEXP,DELPX,DELPY,DELPZ,DELPE,IORIG
1019 C ENDIF
1020  1000 FORMAT(' CHECKF: PX,PY,PZ,PE,EEXT,EEXP',2f7.3,2e12.3,2f7.3
1021  * / 8x,' DELPX/Y/Z/E',4f7.3,i10,' IORIG')
1022 C WRITE(6,'(8X,A,6F8.3)') ' TASUMA,TASUBI,TABI,TAMASU,TAMA,TAIMMA',
1023 C +TASUMA,TASUBI,TABI,TAMASU,TAMA,TAIMMA
1024 
1025  IF(ipri.GT.1) THEN
1026  DO 20 i=1,nhkk
1027  IF(isthkk(i).EQ.11) THEN
1028  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1029  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1030  + (vhkk(khkk,i),khkk=1,4)
1031 
1032  ENDIF
1033  IF(isthkk(i).EQ.12) THEN
1034  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1035  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1036  + (vhkk(khkk,i),khkk=1,4)
1037 
1038  ENDIF
1039  IF(isthkk(i).EQ.1) THEN
1040  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1041  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1042  + (vhkk(khkk,i),khkk=1,4)
1043 
1044  1010 FORMAT (i6,i4,5i6,9(1pe10.2))
1045  ENDIF
1046  IF(isthkk(i).EQ.16) THEN
1047  imo=jmohkk(1,i)
1048  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1049  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1050  + (vhkk(khkk,i),khkk=1,4)
1051 
1052  ENDIF
1053  20 CONTINUE
1054  ENDIF
1055 C ENDIF
1056  RETURN
1057  END
1058 C
1059 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1060 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1061 C
1062  SUBROUTINE checkn(EPN,PPN,IREJ,IORIG)
1063  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1064  SAVE
1065 *KEEP,HKKEVT.
1066 c INCLUDE (HKKEVT)
1067  parameter(amuamu=0.93149432d0)
1068  parameter(nmxhkk= 89998)
1069 c PARAMETER (NMXHKK=25000)
1070  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
1071  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
1072  +(4,nmxhkk)
1073 C
1074 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
1075 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
1076 C THE POSITIONS OF THE PROJECTILE NUCLEONS
1077 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
1078 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
1079 C COMPLETELY CONSISTENT. THE TIMES IN THE
1080 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
1081 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
1082 C
1083 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
1084 C
1085 C NMXHKK: maximum numbers of entries (partons/particles) that can be
1086 C stored in the commonblock.
1087 C
1088 C NHKK: the actual number of entries stored in current event. These are
1089 C found in the first NHKK positions of the respective arrays below.
1090 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
1091 C entry.
1092 C
1093 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
1094 C = 0 : null entry.
1095 C = 1 : an existing entry, which has not decayed or fragmented.
1096 C This is the main class of entries which represents the
1097 C "final state" given by the generator.
1098 C = 2 : an entry which has decayed or fragmented and therefore
1099 C is not appearing in the final state, but is retained for
1100 C event history information.
1101 C = 3 : a documentation line, defined separately from the event
1102 C history. (incoming reacting
1103 C particles, etc.)
1104 C = 4 - 10 : undefined, but reserved for future standards.
1105 C = 11 - 20 : at the disposal of each model builder for constructs
1106 C specific to his program, but equivalent to a null line in the
1107 C context of any other program. One example is the cone defining
1108 C vector of HERWIG, another cluster or event axes of the JETSET
1109 C analysis routines.
1110 C = 21 - : at the disposal of users, in particular for event tracking
1111 C in the detector.
1112 C
1113 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
1114 C standard.
1115 C
1116 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
1117 C The value is 0 for initial entries.
1118 C
1119 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
1120 C one mother exist, in which case the value 0 is used. In cluster
1121 C fragmentation models, the two mothers would correspond to the q
1122 C and qbar which join to form a cluster. In string fragmentation,
1123 C the two mothers of a particle produced in the fragmentation would
1124 C be the two endpoints of the string (with the range in between
1125 C implied).
1126 C
1127 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
1128 C entry has not decayed, this is 0.
1129 C
1130 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
1131 C entry has not decayed, this is 0. It is assumed that the daughters
1132 C of a particle (or cluster or string) are stored sequentially, so
1133 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
1134 C daughters. Even in cases where only one daughter is defined (e.g.
1135 C K0 -> K0S) both values should be defined, to make for a uniform
1136 C approach in terms of loop constructions.
1137 C
1138 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
1139 C
1140 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
1141 C
1142 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
1143 C
1144 C PHKK(4,IHKK) : energy, in GeV.
1145 C
1146 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
1147 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
1148 C
1149 C VHKK(1,IHKK) : production vertex x position, in mm.
1150 C
1151 C VHKK(2,IHKK) : production vertex y position, in mm.
1152 C
1153 C VHKK(3,IHKK) : production vertex z position, in mm.
1154 C
1155 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
1156 C********************************************************************
1157 *KEEP,DELP.
1158  COMMON /delp/ delpx,delpy,delpz,delpe
1159 *KEEP,TANUIN.
1160  COMMON /tanuin/ tasuma,tasubi,tabi,tamasu,tama,taimma
1161 *KEEP,NUCC.
1162  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
1163 *KEEP,DPRIN.
1164  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1165 *KEND.
1166  DATA icheck/0/
1167  help=log10(epn)
1168  phelp=0.
1169  IF(help.GT.5.d0)phelp=help-5.d0
1170  pthelp=12.d0+phelp*5.d0
1171  irej=0
1172  irejj=0
1173  px=0.d0
1174  py=0.d0
1175  pz=0.d0
1176  pe=0.d0
1177  eext=0.d0
1178  eexp=0.d0
1179  eee1=0.d0
1180  eeem1=0.d0
1181  ee1001=0.d0
1182  pz1=0.d0
1183  pzm1=0.d0
1184  pz1001=0.d0
1185  px1=0.d0
1186  pxm1=0.0
1187  px1001=0.d0
1188  DO 10 i=1,nhkk
1189 C Projectiles
1190  IF(isthkk(i).EQ.11.OR.isthkk(i).EQ.13) THEN
1191  gam=epn/phkk(5,i)
1192  bgam=ppn/phkk(5,i)
1193  px=px - phkk(1,i)
1194  py=py - phkk(2,i)
1195  pz=pz - gam*phkk(3,i) - bgam*phkk(4,i)
1196  pe=pe - gam*phkk(4,i) - bgam*phkk(3,i)
1197  ENDIF
1198 C Target
1199  IF(isthkk(i).EQ.12.OR.isthkk(i).EQ.14) THEN
1200  px=px - phkk(1,i)
1201  py=py - phkk(2,i)
1202  pz=pz - phkk(3,i)
1203  pe=pe - phkk(4,i)
1204  ENDIF
1205 * sum final state momenta
1206  IF(isthkk(i).EQ.1) THEN
1207  px=px + phkk(1,i)
1208  py=py + phkk(2,i)
1209  pz=pz + phkk(3,i)
1210  pe=pe + phkk(4,i)
1211  ENDIF
1212 C noninteracting Projectiles
1213  IF(isthkk(i).EQ.13.AND.jdahkk(2,i).EQ.0) THEN
1214  gam=epn/phkk(5,i)
1215  bgam=ppn/phkk(5,i)
1216  px=px + phkk(1,i)
1217  py=py + phkk(2,i)
1218  pz=pz + gam*phkk(3,i) + bgam*phkk(4,i)
1219  pe=pe + gam*phkk(4,i) + bgam*phkk(3,i)
1220  ENDIF
1221  IF(isthkk(i).EQ.14.AND.jdahkk(2,i).EQ.0) THEN
1222 C noninteracting Targets
1223  px=px + phkk(1,i)
1224  py=py + phkk(2,i)
1225  pz=pz + phkk(3,i)
1226  pe=pe + phkk(4,i)
1227  ENDIF
1228  IF(isthkk(i).EQ.16) THEN
1229  imo=jmohkk(1,i)
1230  px=px + phkk(1,i)
1231  py=py + phkk(2,i)
1232  pz=pz + phkk(3,i)
1233  pe=pe + phkk(4,i)
1234  eext=eext + phkk(4,i) - phkk(4,imo)
1235  ENDIF
1236  IF(isthkk(i).EQ.15) THEN
1237  imo=jmohkk(1,i)
1238  gam=epn/phkk(5,i)
1239  bgam=ppn/phkk(5,i)
1240  px=px + phkk(1,i)
1241  py=py + phkk(2,i)
1242  pz=pz + gam*phkk(3,i) + bgam*phkk(4,i)
1243  pe=pe + gam*phkk(4,i) + bgam*phkk(3,i)
1244  eext=eext + phkk(4,i) - phkk(4,imo)
1245  ENDIF
1246  IF(isthkk(i).EQ.1) THEN
1247  eee1=eee1+phkk(4,i)
1248  pz1=pz1+phkk(3,i)
1249  px1=px1+phkk(1,i)
1250  ENDIF
1251  IF(isthkk(i).EQ.-1) THEN
1252  eeem1=eeem1+phkk(4,i)
1253  pzm1=pzm1+phkk(3,i)
1254  pxm1=pxm1+phkk(1,i)
1255  ENDIF
1256  IF(isthkk(i).EQ.1001) THEN
1257  ee1001=ee1001+phkk(4,i)
1258  pz1001=pz1001+phkk(3,i)
1259  px1001=px1001+phkk(1,i)
1260  ENDIF
1261  10 CONTINUE
1262  eee=eee1+eeem1+ee1001
1263  pzpz=pz1+pzm1+pz1001
1264  pxpx=px1+pxm1+px1001
1265 C--------------------------------------------------------------
1266 C
1267 C patch to correct pz of residual nuclei
1268 C
1269 C--------------------------------------------------------------
1270  delpz=ppn-pzpz
1271  pzpz=pzpz+delpz
1272  pz1001=pz1001+delpz
1273  ee1001=0.d0
1274  DO 101 i=1,nhkk
1275  IF(isthkk(i).EQ.1001) THEN
1276  phkk(3,i)=phkk(3,i)+delpz
1277  phkk(4,i)=sqrt(phkk(1,i)**2+phkk(2,i)**2+phkk(3,i)**2
1278  * +phkk(5,i)**2)
1279  ee1001=ee1001+phkk(4,i)
1280  ENDIF
1281  101 CONTINUE
1282  eee=eee1+eeem1+ee1001
1283 C--------------------------------------------------------------
1284  aip=1
1285  epnto=epn*aip
1286  eeee=eee-epn*aip
1287  aip=1
1288  ait=it
1289  aitz=itz
1290  bip=epn+(ait*amuamu+1.d-3*energy(ait,aitz))
1291  bmi=1.d-3*energy(ait,aitz)
1292  delle=abs(bip-eee)
1293 C TOLE=EPN/450000.D0
1294 C TOLE=EPN/2500.D0
1295  tole=0.16d0
1296  IF(delle.GE.tole)irej=1
1297  IF(irej.EQ.1)THEN
1298  icheck=icheck+1
1299  IF(icheck.LE.20)THEN
1300  WRITE(6,'(A,I5,E10.3,4F10.4)')
1301  * ' IP,EPN,PXPX,PX1,PXM1,PX1001:',
1302  * ip,epn,pxpx,px1,pxm1,px1001
1303  WRITE(6,'(A,I5,E10.3,6F10.4)')
1304  * ' IP,PPN,PZPZ,PZ1,PZM1,PZ1001,BIP,BMI:',
1305  * ip,ppn,pzpz,pz1,pzm1,pz1001,bip,bmi
1306  WRITE(6,'(A,I5,E10.3,5E12.4)')
1307  * ' IP,EPN,EEE,EEE1,EEEM1,EE1001,DELLE:',
1308  * ip,epn,eee,eee1,eeem1,ee1001,delle
1309  ENDIF
1310  ENDIF
1311 C PX=PX + DELPX
1312 C PY=PY + DELPY
1313 C PZ=PZ + DELPZ
1314 C PE=PE + DELPE
1315 C IF(IPRI.GT.1) THEN
1316 C IF (ABS(PX).GT.PTHELP.OR. ABS(PY).GT.PTHELP.OR.
1317 C * ABS(PZ)/EPN.GT.0.025D0*IP.
1318 C + OR. ABS(PE)/EPN.GT.0.025D0*IP) THEN
1319 C IREJJ=1
1320 C ICHECK=ICHECK+1
1321 C IF(ICHECK.LE.500.AND.IREJJ.EQ.1)THEN
1322 C WRITE(6,1000) PX,PY,PZ,PE,EEXT,EEXP,DELPX,DELPY,DELPZ,DELPE,IORIG
1323 C ENDIF
1324  1000 FORMAT(' CHECKF: PX,PY,PZ,PE,EEXT,EEXP',2f7.3,2e12.3,2f7.3
1325  * / 8x,' DELPX/Y/Z/E',4f7.3,i10,' IORIG')
1326 C WRITE(6,'(8X,A,6F8.3)') ' TASUMA,TASUBI,TABI,TAMASU,TAMA,TAIMMA',
1327 C +TASUMA,TASUBI,TABI,TAMASU,TAMA,TAIMMA
1328 
1329  IF(ipri.GT.1) THEN
1330  DO 20 i=1,nhkk
1331  IF(isthkk(i).EQ.11) THEN
1332  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1333  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1334  + (vhkk(khkk,i),khkk=1,4)
1335 
1336  ENDIF
1337  IF(isthkk(i).EQ.12) THEN
1338  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1339  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1340  + (vhkk(khkk,i),khkk=1,4)
1341 
1342  ENDIF
1343  IF(isthkk(i).EQ.1) THEN
1344  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1345  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1346  + (vhkk(khkk,i),khkk=1,4)
1347 
1348  1010 FORMAT (i6,i4,5i6,9(1pe10.2))
1349  ENDIF
1350  IF(isthkk(i).EQ.16) THEN
1351  imo=jmohkk(1,i)
1352  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1353  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1354  + (vhkk(khkk,i),khkk=1,4)
1355 
1356  ENDIF
1357  20 CONTINUE
1358  ENDIF
1359 C ENDIF
1360  RETURN
1361  END
1362 C
1363 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1364 C
1365  SUBROUTINE checko(EPN,PPN,IREJ,IORIG)
1366  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1367  SAVE
1368 *KEEP,HKKEVT.
1369 c INCLUDE (HKKEVT)
1370  parameter(amuamu=0.93149432d0)
1371  parameter(nmxhkk= 89998)
1372 c PARAMETER (NMXHKK=25000)
1373  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
1374  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
1375  +(4,nmxhkk)
1376 C
1377 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
1378 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
1379 C THE POSITIONS OF THE PROJECTILE NUCLEONS
1380 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
1381 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
1382 C COMPLETELY CONSISTENT. THE TIMES IN THE
1383 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
1384 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
1385 C
1386 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
1387 C
1388 C NMXHKK: maximum numbers of entries (partons/particles) that can be
1389 C stored in the commonblock.
1390 C
1391 C NHKK: the actual number of entries stored in current event. These are
1392 C found in the first NHKK positions of the respective arrays below.
1393 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
1394 C entry.
1395 C
1396 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
1397 C = 0 : null entry.
1398 C = 1 : an existing entry, which has not decayed or fragmented.
1399 C This is the main class of entries which represents the
1400 C "final state" given by the generator.
1401 C = 2 : an entry which has decayed or fragmented and therefore
1402 C is not appearing in the final state, but is retained for
1403 C event history information.
1404 C = 3 : a documentation line, defined separately from the event
1405 C history. (incoming reacting
1406 C particles, etc.)
1407 C = 4 - 10 : undefined, but reserved for future standards.
1408 C = 11 - 20 : at the disposal of each model builder for constructs
1409 C specific to his program, but equivalent to a null line in the
1410 C context of any other program. One example is the cone defining
1411 C vector of HERWIG, another cluster or event axes of the JETSET
1412 C analysis routines.
1413 C = 21 - : at the disposal of users, in particular for event tracking
1414 C in the detector.
1415 C
1416 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
1417 C standard.
1418 C
1419 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
1420 C The value is 0 for initial entries.
1421 C
1422 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
1423 C one mother exist, in which case the value 0 is used. In cluster
1424 C fragmentation models, the two mothers would correspond to the q
1425 C and qbar which join to form a cluster. In string fragmentation,
1426 C the two mothers of a particle produced in the fragmentation would
1427 C be the two endpoints of the string (with the range in between
1428 C implied).
1429 C
1430 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
1431 C entry has not decayed, this is 0.
1432 C
1433 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
1434 C entry has not decayed, this is 0. It is assumed that the daughters
1435 C of a particle (or cluster or string) are stored sequentially, so
1436 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
1437 C daughters. Even in cases where only one daughter is defined (e.g.
1438 C K0 -> K0S) both values should be defined, to make for a uniform
1439 C approach in terms of loop constructions.
1440 C
1441 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
1442 C
1443 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
1444 C
1445 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
1446 C
1447 C PHKK(4,IHKK) : energy, in GeV.
1448 C
1449 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
1450 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
1451 C
1452 C VHKK(1,IHKK) : production vertex x position, in mm.
1453 C
1454 C VHKK(2,IHKK) : production vertex y position, in mm.
1455 C
1456 C VHKK(3,IHKK) : production vertex z position, in mm.
1457 C
1458 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
1459 C********************************************************************
1460  COMMON /zentra/ icentr
1461 *KEEP,DELP.
1462  COMMON /delp/ delpx,delpy,delpz,delpe
1463 *KEEP,TANUIN.
1464  COMMON /tanuin/ tasuma,tasubi,tabi,tamasu,tama,taimma
1465 *KEEP,NUCC.
1466  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
1467 *KEEP,DPRIN.
1468  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1469 *KEND.
1470  DATA icheck/0/
1471  help=log10(epn)
1472  phelp=0.
1473  IF(help.GT.5.d0)phelp=help-5.
1474  pthelp=12.+phelp*5.
1475  irej=0
1476  irejj=0
1477  px=0.
1478  py=0.
1479  pz=0.
1480  pe=0.
1481  eext=0.
1482  eexp=0.
1483  eee1=0.
1484  eeem1=0.
1485  ee1001=0.
1486  DO 10 i=1,nhkk
1487 C Projectiles
1488  IF(isthkk(i).EQ.11.OR.isthkk(i).EQ.13) THEN
1489  gam=epn/phkk(5,i)
1490  bgam=ppn/phkk(5,i)
1491  px=px - phkk(1,i)
1492  py=py - phkk(2,i)
1493  pz=pz - gam*phkk(3,i) - bgam*phkk(4,i)
1494  pe=pe - gam*phkk(4,i) - bgam*phkk(3,i)
1495  ENDIF
1496 C Target
1497  IF(isthkk(i).EQ.12.OR.isthkk(i).EQ.14) THEN
1498  px=px - phkk(1,i)
1499  py=py - phkk(2,i)
1500  pz=pz - phkk(3,i)
1501  pe=pe - phkk(4,i)
1502  ENDIF
1503 * sum final state momenta
1504  IF(isthkk(i).EQ.1) THEN
1505  px=px + phkk(1,i)
1506  py=py + phkk(2,i)
1507  pz=pz + phkk(3,i)
1508  pe=pe + phkk(4,i)
1509  ENDIF
1510 C noninteracting Projectiles
1511  IF(isthkk(i).EQ.13.AND.jdahkk(2,i).EQ.0) THEN
1512  gam=epn/phkk(5,i)
1513  bgam=ppn/phkk(5,i)
1514  px=px + phkk(1,i)
1515  py=py + phkk(2,i)
1516  pz=pz + gam*phkk(3,i) + bgam*phkk(4,i)
1517  pe=pe + gam*phkk(4,i) + bgam*phkk(3,i)
1518  ENDIF
1519  IF(isthkk(i).EQ.14.AND.jdahkk(2,i).EQ.0) THEN
1520 C noninteracting Targets
1521  px=px + phkk(1,i)
1522  py=py + phkk(2,i)
1523  pz=pz + phkk(3,i)
1524  pe=pe + phkk(4,i)
1525  ENDIF
1526  IF(isthkk(i).EQ.16) THEN
1527  imo=jmohkk(1,i)
1528  px=px + phkk(1,i)
1529  py=py + phkk(2,i)
1530  pz=pz + phkk(3,i)
1531  pe=pe + phkk(4,i)
1532  eext=eext + phkk(4,i) - phkk(4,imo)
1533  ENDIF
1534  IF(isthkk(i).EQ.15) THEN
1535  imo=jmohkk(1,i)
1536  gam=epn/phkk(5,i)
1537  bgam=ppn/phkk(5,i)
1538  px=px + phkk(1,i)
1539  py=py + phkk(2,i)
1540  pz=pz + gam*phkk(3,i) + bgam*phkk(4,i)
1541  pe=pe + gam*phkk(4,i) + bgam*phkk(3,i)
1542  eext=eext + phkk(4,i) - phkk(4,imo)
1543  ENDIF
1544  IF(isthkk(i).EQ.1) THEN
1545  eee1=eee1+phkk(3,i)
1546  ENDIF
1547  10 CONTINUE
1548  eee=eee1
1549  epnto=ppn*ip
1550  aeee=eee/epnto
1551  eeee=eee-ppn*ip
1552  aeeee=eee/ppn-ip
1553  aeee1=eee1/epnto
1554  aip=1
1555  ait=it
1556  aitz=itz
1557  aip=aip
1558  delle=abs(aip-aeee)
1559  elle=delle*epnto
1560 C IF(DELLE.GE.0.025)IREJ=1
1561 C IF(IREJ.EQ.1)
1562 C * WRITE(6,'(A,I5,E10.3,3F10.4)')
1563 C *' IP,EPN,AEEE,AEEEE,AEEE1:',
1564 C * IP,EPN,AEEE,AEEEE,AEEE1
1565 C IF(IREJ.EQ.1)
1566 C * WRITE(6,'(A,I5,E10.3,5F10.4)')
1567 C *' IP,EPN,EEE,EEEE,EEE1,DELLE,ELLE:',
1568 C * IP,EPN,EEE,EEEE,EEE1,DELLE,ELLE
1569  px=px + delpx
1570  py=py + delpy
1571  pz=pz + delpz
1572  pe=pe + delpe
1573  tole=0.025d0*ip
1574  IF(ip.EQ.it.AND.it.GT.1)tole=0.05d0*ip
1575 C IF(ICENTR.EQ.1)TOLE=TOLE*2.
1576  IF(epn.LE.5.d0)tole=3.d0*tole
1577 C IF(IPRI.GT.1) THEN
1578  IF (abs(px).GT.pthelp.OR. abs(py).GT.pthelp.OR.
1579  * abs(pz)/epn.GT.tole.
1580  + or. abs(pe)/epn.GT.tole) THEN
1581  irej=1
1582  icheck=icheck+1
1583  IF(icheck.LE.50.AND.irej.EQ.1)THEN
1584  WRITE(6,1000) px,py,pz,pe,eext,eexp,delpx,delpy,delpz,delpe,iorig
1585  ENDIF
1586  1000 FORMAT(' CHECKO: PX,PY,PZ,PE,EEXT,EEXP',2f7.3,2e12.3,2f7.3
1587  * / 8x,' DELPX/Y/Z/E',4f7.3,i10,' IORIG')
1588 C WRITE(6,'(8X,A,6F8.3)') ' TASUMA,TASUBI,TABI,TAMASU,TAMA,TAIMMA',
1589 C +TASUMA,TASUBI,TABI,TAMASU,TAMA,TAIMMA
1590  IF(ipri.GE.1)THEN
1591  WRITE(6,1000) px,py,pz,pe,eext,eexp,delpx,delpy,delpz,delpe,iorig
1592  ENDIF
1593  IF(ipri.GT.1) THEN
1594  DO 20 i=1,nhkk
1595  IF(isthkk(i).EQ.11) THEN
1596  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1597  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1598  + (vhkk(khkk,i),khkk=1,4)
1599 
1600  ENDIF
1601  IF(isthkk(i).EQ.12) THEN
1602  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1603  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1604  + (vhkk(khkk,i),khkk=1,4)
1605 
1606  ENDIF
1607  IF(isthkk(i).EQ.1) THEN
1608  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1609  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1610  + (vhkk(khkk,i),khkk=1,4)
1611 
1612  1010 FORMAT (i6,i4,5i6,9(1pe10.2))
1613  ENDIF
1614  IF(isthkk(i).EQ.16) THEN
1615  imo=jmohkk(1,i)
1616  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1617  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1618  + (vhkk(khkk,i),khkk=1,4)
1619 
1620  ENDIF
1621  20 CONTINUE
1622  ENDIF
1623  ENDIF
1624  IF(ipri.GE.1)THEN
1625  WRITE(6,1000) px,py,pz,pe,eext,eexp,delpx,delpy,delpz,delpe,iorig
1626  ENDIF
1627  RETURN
1628  END
1629 C
1630  SUBROUTINE checke(EPN,PPN)
1631  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1632  SAVE
1633 *KEEP,HKKEVT.
1634 c INCLUDE (HKKEVT)
1635  parameter(nmxhkk= 89998)
1636 c PARAMETER (NMXHKK=25000)
1637  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
1638  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
1639  +(4,nmxhkk)
1640 C
1641 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
1642 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
1643 C THE POSITIONS OF THE PROJECTILE NUCLEONS
1644 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
1645 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
1646 C COMPLETELY CONSISTENT. THE TIMES IN THE
1647 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
1648 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
1649 C
1650 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
1651 C
1652 C NMXHKK: maximum numbers of entries (partons/particles) that can be
1653 C stored in the commonblock.
1654 C
1655 C NHKK: the actual number of entries stored in current event. These are
1656 C found in the first NHKK positions of the respective arrays below.
1657 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
1658 C entry.
1659 C
1660 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
1661 C = 0 : null entry.
1662 C = 1 : an existing entry, which has not decayed or fragmented.
1663 C This is the main class of entries which represents the
1664 C "final state" given by the generator.
1665 C = 2 : an entry which has decayed or fragmented and therefore
1666 C is not appearing in the final state, but is retained for
1667 C event history information.
1668 C = 3 : a documentation line, defined separately from the event
1669 C history. (incoming reacting
1670 C particles, etc.)
1671 C = 4 - 10 : undefined, but reserved for future standards.
1672 C = 11 - 20 : at the disposal of each model builder for constructs
1673 C specific to his program, but equivalent to a null line in the
1674 C context of any other program. One example is the cone defining
1675 C vector of HERWIG, another cluster or event axes of the JETSET
1676 C analysis routines.
1677 C = 21 - : at the disposal of users, in particular for event tracking
1678 C in the detector.
1679 C
1680 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
1681 C standard.
1682 C
1683 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
1684 C The value is 0 for initial entries.
1685 C
1686 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
1687 C one mother exist, in which case the value 0 is used. In cluster
1688 C fragmentation models, the two mothers would correspond to the q
1689 C and qbar which join to form a cluster. In string fragmentation,
1690 C the two mothers of a particle produced in the fragmentation would
1691 C be the two endpoints of the string (with the range in between
1692 C implied).
1693 C
1694 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
1695 C entry has not decayed, this is 0.
1696 C
1697 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
1698 C entry has not decayed, this is 0. It is assumed that the daughters
1699 C of a particle (or cluster or string) are stored sequentially, so
1700 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
1701 C daughters. Even in cases where only one daughter is defined (e.g.
1702 C K0 -> K0S) both values should be defined, to make for a uniform
1703 C approach in terms of loop constructions.
1704 C
1705 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
1706 C
1707 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
1708 C
1709 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
1710 C
1711 C PHKK(4,IHKK) : energy, in GeV.
1712 C
1713 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
1714 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
1715 C
1716 C VHKK(1,IHKK) : production vertex x position, in mm.
1717 C
1718 C VHKK(2,IHKK) : production vertex y position, in mm.
1719 C
1720 C VHKK(3,IHKK) : production vertex z position, in mm.
1721 C
1722 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
1723 C********************************************************************
1724 *KEEP,DELP.
1725  COMMON /delp/ delpx,delpy,delpz,delpe
1726 *KEEP,TANUIN.
1727  COMMON /tanuin/ tasuma,tasubi,tabi,tamasu,tama,taimma
1728 *KEEP,DPRIN.
1729  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1730 *KEND.
1731  px=0.
1732  py=0.
1733  pz=0.
1734  pe=0.
1735  eext=0.
1736  eexp=0.
1737  DO 10 i=1,nhkk
1738 C Projectiles
1739  IF(isthkk(i).EQ.11.OR.isthkk(i).EQ.13) THEN
1740  gam=epn/phkk(5,i)
1741  bgam=ppn/phkk(5,i)
1742  px=px - phkk(1,i)
1743  py=py - phkk(2,i)
1744  pz=pz - gam*phkk(3,i) - bgam*phkk(4,i)
1745  pe=pe - gam*phkk(4,i) - bgam*phkk(3,i)
1746  ENDIF
1747 C Target
1748  IF(isthkk(i).EQ.12.OR.isthkk(i).EQ.14) THEN
1749  px=px - phkk(1,i)
1750  py=py - phkk(2,i)
1751  pz=pz - phkk(3,i)
1752  pe=pe - phkk(4,i)
1753  ENDIF
1754 * sum final state momenta
1755  IF(isthkk(i).EQ.1) THEN
1756  px=px + phkk(1,i)
1757  py=py + phkk(2,i)
1758  pz=pz + phkk(3,i)
1759  pe=pe + phkk(4,i)
1760  ENDIF
1761 C noninteracting Projectiles
1762  IF(isthkk(i).EQ.13.AND.jdahkk(1,i).EQ.0) THEN
1763  gam=epn/phkk(5,i)
1764  bgam=ppn/phkk(5,i)
1765  px=px + phkk(1,i)
1766  py=py + phkk(2,i)
1767  pz=pz + gam*phkk(3,i) + bgam*phkk(4,i)
1768  pe=pe + gam*phkk(4,i) + bgam*phkk(3,i)
1769  ENDIF
1770  IF(isthkk(i).EQ.14.AND.jdahkk(1,i).EQ.0) THEN
1771 C noninteracting Targets
1772  px=px + phkk(1,i)
1773  py=py + phkk(2,i)
1774  pz=pz + phkk(3,i)
1775  pe=pe + phkk(4,i)
1776  ENDIF
1777  IF(isthkk(i).EQ.16) THEN
1778  imo=jmohkk(1,i)
1779  px=px + phkk(1,i)
1780  py=py + phkk(2,i)
1781  pz=pz + phkk(3,i)
1782  pe=pe + phkk(4,i)
1783  eext=eext + phkk(4,i) - phkk(4,imo)
1784  ENDIF
1785  IF(isthkk(i).EQ.15) THEN
1786  imo=jmohkk(1,i)
1787  px=px + phkk(1,i)
1788  py=py + phkk(2,i)
1789  pz=pz + phkk(3,i)
1790  pe=pe + phkk(4,i)
1791  eext=eext + phkk(4,i) - phkk(4,imo)
1792  ENDIF
1793  10 CONTINUE
1794  px=px + delpx
1795  py=py + delpy
1796  pz=pz + delpz
1797  pe=pe + delpe
1798  WRITE(6,1000) px,py,pz,pe,eext,eexp,delpx,delpy,delpz,delpe
1799  1000 FORMAT(' CHECKE: PX,PY,PZ,PE,EEXT,EEXP',6f7.3/ 8x,' DELPX/Y/Z/E',4
1800  +f7.3)
1801  WRITE(6,'(8X,A,6F8.3)') ' TASUMA,TASUBI,TABI,TAMASU,TAMA,TAIMMA',
1802  +tasuma,tasubi,tabi,tamasu,tama,taimma
1803  IF(ipri.GT.1) THEN
1804  IF (abs(px).GT.0.004.OR. abs(py).GT.0.004.OR. abs(pz).GT.0.004.
1805  + or. abs(pe).GT.0.004) THEN
1806 
1807 
1808  DO 20 i=1,nhkk
1809  IF(isthkk(i).EQ.11) THEN
1810  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1811  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1812  + (vhkk(khkk,i),khkk=1,4)
1813 
1814  ENDIF
1815  IF(isthkk(i).EQ.12) THEN
1816  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1817  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1818  + (vhkk(khkk,i),khkk=1,4)
1819 
1820  ENDIF
1821  IF(isthkk(i).EQ.1) THEN
1822  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1823  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1824  + (vhkk(khkk,i),khkk=1,4)
1825 
1826  1010 FORMAT (i6,i4,5i6,9(1pe10.2))
1827  ENDIF
1828  IF(isthkk(i).EQ.16) THEN
1829  imo=jmohkk(1,i)
1830  WRITE(6,1010) i,isthkk(i),idhkk(i),jmohkk(1,i),jmohkk
1831  + (2,i), jdahkk(1,i),jdahkk(2,i),(phkk(khkk,i),khkk=1,5),
1832  + (vhkk(khkk,i),khkk=1,4)
1833 
1834  ENDIF
1835  20 CONTINUE
1836  ENDIF
1837  ENDIF
1838  RETURN
1839  END
1840 C######################################################################
1841 C######################################################################
1842 C######################################################################
1843 *
1844 C######################################################################
1845 C######################################################################
1846 C######################################################################
1847 C
1848 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1849 C
1850 C######################################################################
1851 C######################################################################
1852 C######################################################################
1853 C######################################################################
1854 C######################################################################
1855 C######################################################################
1856 C
1857 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1858 C
1859  DOUBLE PRECISION FUNCTION ebind(IA,IZ)
1860  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1861  SAVE
1862 C***
1863 C Binding energy for nuclei with mass number IA
1864 C and atomic number IZ
1865 C from Shirokov & Yudin, Yad. Fizika, Nauka, Moskva 1972
1866 C***
1867  DATA a1,a2,a3,a4,a5 /0.01575, 0.0178, 0.000710, 0.0237, 0.034/
1868 C
1869 C WRITE (6,'(A,2I5)')' EBIND IA,IZ ',IA,IZ
1870  IF(ia.LE.1.OR.iz.EQ.0)THEN
1871  ebind=0
1872  RETURN
1873  ENDIF
1874  aa=ia
1875  ebind = a1*aa - a2*aa**0. 666667- a3*iz*iz*aa**(-0.333333) - a4
1876  +*(ia-2*iz)**2/aa
1877  IF (mod(ia,2).EQ.1) THEN
1878  ia5=0
1879  ELSEIF (mod(iz,2).EQ.1) THEN
1880  ia5=1
1881  ELSE
1882  ia5=-1
1883  ENDIF
1884  ebind=ebind - ia5*a5*aa**(-0.75)
1885  RETURN
1886  END
1887 C
1888 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1889 C
1890  SUBROUTINE defaul(EPN,PPN)
1891  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1892  SAVE
1893 *---set default values for some parameters
1894 *KEEP,HKKEVT.
1895 c INCLUDE (HKKEVT)
1896  parameter(nmxhkk= 89998)
1897 c PARAMETER (NMXHKK=25000)
1898  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
1899  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
1900  +(4,nmxhkk)
1901 C
1902 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
1903 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
1904 C THE POSITIONS OF THE PROJECTILE NUCLEONS
1905 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
1906 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
1907 C COMPLETELY CONSISTENT. THE TIMES IN THE
1908 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
1909 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
1910 C
1911 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
1912 C
1913 C NMXHKK: maximum numbers of entries (partons/particles) that can be
1914 C stored in the commonblock.
1915 C
1916 C NHKK: the actual number of entries stored in current event. These are
1917 C found in the first NHKK positions of the respective arrays below.
1918 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
1919 C entry.
1920 C
1921 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
1922 C = 0 : null entry.
1923 C = 1 : an existing entry, which has not decayed or fragmented.
1924 C This is the main class of entries which represents the
1925 C "final state" given by the generator.
1926 C = 2 : an entry which has decayed or fragmented and therefore
1927 C is not appearing in the final state, but is retained for
1928 C event history information.
1929 C = 3 : a documentation line, defined separately from the event
1930 C history. (incoming reacting
1931 C particles, etc.)
1932 C = 4 - 10 : undefined, but reserved for future standards.
1933 C = 11 - 20 : at the disposal of each model builder for constructs
1934 C specific to his program, but equivalent to a null line in the
1935 C context of any other program. One example is the cone defining
1936 C vector of HERWIG, another cluster or event axes of the JETSET
1937 C analysis routines.
1938 C = 21 - : at the disposal of users, in particular for event tracking
1939 C in the detector.
1940 C
1941 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
1942 C standard.
1943 C
1944 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
1945 C The value is 0 for initial entries.
1946 C
1947 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
1948 C one mother exist, in which case the value 0 is used. In cluster
1949 C fragmentation models, the two mothers would correspond to the q
1950 C and qbar which join to form a cluster. In string fragmentation,
1951 C the two mothers of a particle produced in the fragmentation would
1952 C be the two endpoints of the string (with the range in between
1953 C implied).
1954 C
1955 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
1956 C entry has not decayed, this is 0.
1957 C
1958 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
1959 C entry has not decayed, this is 0. It is assumed that the daughters
1960 C of a particle (or cluster or string) are stored sequentially, so
1961 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
1962 C daughters. Even in cases where only one daughter is defined (e.g.
1963 C K0 -> K0S) both values should be defined, to make for a uniform
1964 C approach in terms of loop constructions.
1965 C
1966 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
1967 C
1968 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
1969 C
1970 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
1971 C
1972 C PHKK(4,IHKK) : energy, in GeV.
1973 C
1974 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
1975 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
1976 C
1977 C VHKK(1,IHKK) : production vertex x position, in mm.
1978 C
1979 C VHKK(2,IHKK) : production vertex y position, in mm.
1980 C
1981 C VHKK(3,IHKK) : production vertex z position, in mm.
1982 C
1983 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
1984 C********************************************************************
1985 *KEEP,DPAR.
1986 C /DPAR/ CONTAINS PARTICLE PROPERTIES
1987 C ANAME = LITERAL NAME OF THE PARTICLE
1988 C AAM = PARTICLE MASS IN GEV
1989 C GA = DECAY WIDTH
1990 C TAU = LIFE TIME OF INSTABLE PARTICLES
1991 C IICH = ELECTRIC CHARGE OF THE PARTICLE
1992 C IIBAR = BARYON NUMBER
1993 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
1994 C
1995  CHARACTER*8 aname
1996  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
1997  +iibar(210),k1(210),k2(210)
1998 C------------------
1999 *KEEP,FACTMO.
2000  COMMON /factmo/ ifacto
2001 *KEEP,TAUFO.
2002  COMMON /taufo/ taufor,ktauge,itauve,incmod
2003 *KEEP,HADTHR.
2004  COMMON /hadthr/ ehadth,inthad
2005 *KEEP,NUCC.
2006 C COMMON /NUCCC/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
2007 C COMMON /NUCC/ JT,JTZ,JP,JPZ,JJPROJ,JBPROJ,JJTARG,JBTARG
2008  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
2009  COMMON /nuccc/ jt,jtz,jp,jpz,jjproj,jbproj,jjtarg,jbtarg
2010 *KEEP,ZENTRA.
2011  COMMON /zentra/ icentr
2012 *KEEP,CMHICO.
2013  COMMON /cmhico/ cmhis
2014 *KEEP,RESONA.
2015  COMMON /resona/ ireso
2016 *KEEP,XSEADI.
2017  COMMON /xseadi/ xseacu,unon,unom,unosea, cvq,cdq,csea,ssmima,
2018  +ssmimq,vvmthr
2019 *KEEP,COULO.
2020  common/coulo/icoul
2021 *KEEP,EDENS.
2022  common/edens/ieden
2023 *KEEP,PROJK.
2024  COMMON /projk/ iprojk
2025 *KEEP,NUCIMP.
2026  COMMON /nucimp/ prmom(5,248),tamom(5,248), prmfep,prmfen,tamfep,
2027  +tamfen, prefep,prefen,taefep,taefen, prepot(210),taepot(210),
2028  +prebin,taebin,fermod,etacou
2029 *KEND.
2030  COMMON /recom/irecom
2031 C---------------------
2032 *---minimum bias interactions
2033  icentr=0
2034 *---threshold for the use of HADRIN in the primary hadron-nucleon collision
2035  ehadth=5.
2036 *---lab energy and momentum of the projectile: pion+
2037  epn=200.
2038  ijproj=13
2039  jjproj=13
2040  ppn=sqrt((epn-aam(ijproj))*(epn+aam(ijproj)))
2041  ibproj=iibar(ijproj)
2042  ip=1
2043  ipz=1
2044  jbproj=iibar(ijproj)
2045  jp=1
2046  jpz=1
2047 *---copper target
2048  it=14
2049  itz=7
2050  jt=14
2051  jtz=7
2052 *---formation zone intranuclear cascade
2053  taufor=105.
2054  ktauge=0
2055  itauve=1
2056 *---inclusion of Coulomb potential for hA interactions
2057  icoul=1
2058  icoull=1
2059 *---cascade within projectile switched off
2060  iprojk=1
2061 *---nucleus independent meson potential
2062  potmes=0.002
2063  taepot(13)=potmes
2064  taepot(14)=potmes
2065  taepot(15)=potmes
2066  taepot(16)=potmes
2067  taepot(23)=potmes
2068  taepot(24)=potmes
2069  taepot(25)=potmes
2070 *---definition of soft quark distributions
2071  xseacu=0.05
2072  unon=1.11d0
2073  unom=1.11d0
2074  unosea=5.0d0
2075 *---cutoff parameters for x-sampling
2076  cvq=2.0d0
2077  cdq=2.0d0
2078  csea=0.3d0
2079  ssmima=0.90d0
2080  ssmimq=ssmima**2
2081 *---
2082  ireso=0
2083  cmhis=0.d0
2084  ieden=0
2085  ifacto=0
2086 C---Chain recombination
2087 C IRECOM=1
2088  RETURN
2089  END
2090 C
2091 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2092 C
2093 C
2094 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2095 C
2096  SUBROUTINE hadhad(EPN,PPN,NHKKH1,IHTAWW,ITTA,IREJFO)
2097  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2098  SAVE
2099 C*** AT LOW ENERGIES EPN.LE.EHADTH HADRIN IS USED WITH ONE INTERACTION
2100 C*** IN THE NUCLEUS INSTEAD OF THE DUAL PARTON MODEL (GLAUBER CASCADE)
2101 C*** ONLY FOR HADRON-NUCLEUS COLLISIONS!
2102 C
2103 *KEEP,HKKEVT.
2104 c INCLUDE (HKKEVT)
2105  parameter(nmxhkk= 89998)
2106 c PARAMETER (NMXHKK=25000)
2107  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
2108  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
2109  +(4,nmxhkk)
2110 C
2111 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
2112 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
2113 C THE POSITIONS OF THE PROJECTILE NUCLEONS
2114 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
2115 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
2116 C COMPLETELY CONSISTENT. THE TIMES IN THE
2117 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
2118 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
2119 C
2120 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
2121 C
2122 C NMXHKK: maximum numbers of entries (partons/particles) that can be
2123 C stored in the commonblock.
2124 C
2125 C NHKK: the actual number of entries stored in current event. These are
2126 C found in the first NHKK positions of the respective arrays below.
2127 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
2128 C entry.
2129 C
2130 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
2131 C = 0 : null entry.
2132 C = 1 : an existing entry, which has not decayed or fragmented.
2133 C This is the main class of entries which represents the
2134 C "final state" given by the generator.
2135 C = 2 : an entry which has decayed or fragmented and therefore
2136 C is not appearing in the final state, but is retained for
2137 C event history information.
2138 C = 3 : a documentation line, defined separately from the event
2139 C history. (incoming reacting
2140 C particles, etc.)
2141 C = 4 - 10 : undefined, but reserved for future standards.
2142 C = 11 - 20 : at the disposal of each model builder for constructs
2143 C specific to his program, but equivalent to a null line in the
2144 C context of any other program. One example is the cone defining
2145 C vector of HERWIG, another cluster or event axes of the JETSET
2146 C analysis routines.
2147 C = 21 - : at the disposal of users, in particular for event tracking
2148 C in the detector.
2149 C
2150 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
2151 C standard.
2152 C
2153 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
2154 C The value is 0 for initial entries.
2155 C
2156 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
2157 C one mother exist, in which case the value 0 is used. In cluster
2158 C fragmentation models, the two mothers would correspond to the q
2159 C and qbar which join to form a cluster. In string fragmentation,
2160 C the two mothers of a particle produced in the fragmentation would
2161 C be the two endpoints of the string (with the range in between
2162 C implied).
2163 C
2164 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
2165 C entry has not decayed, this is 0.
2166 C
2167 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
2168 C entry has not decayed, this is 0. It is assumed that the daughters
2169 C of a particle (or cluster or string) are stored sequentially, so
2170 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
2171 C daughters. Even in cases where only one daughter is defined (e.g.
2172 C K0 -> K0S) both values should be defined, to make for a uniform
2173 C approach in terms of loop constructions.
2174 C
2175 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
2176 C
2177 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
2178 C
2179 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
2180 C
2181 C PHKK(4,IHKK) : energy, in GeV.
2182 C
2183 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
2184 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
2185 C
2186 C VHKK(1,IHKK) : production vertex x position, in mm.
2187 C
2188 C VHKK(2,IHKK) : production vertex y position, in mm.
2189 C
2190 C VHKK(3,IHKK) : production vertex z position, in mm.
2191 C
2192 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
2193 C********************************************************************
2194 *KEEP,NUCC.
2195  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
2196 *KEEP,DPAR.
2197 C /DPAR/ CONTAINS PARTICLE PROPERTIES
2198 C ANAME = LITERAL NAME OF THE PARTICLE
2199 C AAM = PARTICLE MASS IN GEV
2200 C GA = DECAY WIDTH
2201 C TAU = LIFE TIME OF INSTABLE PARTICLES
2202 C IICH = ELECTRIC CHARGE OF THE PARTICLE
2203 C IIBAR = BARYON NUMBER
2204 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
2205 C
2206  CHARACTER*8 aname
2207  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
2208  +iibar(210),k1(210),k2(210)
2209 C------------------
2210 *KEEP,DPRIN.
2211  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
2212 *KEEP,DFINLS.
2213  parameter(maxfin=10)
2214  COMMON /dfinls/ itrh(maxfin),cxrh(maxfin),cyrh(maxfin), czrh
2215  +(maxfin),elrh(maxfin),plrh(maxfin),irh
2216 *KEEP,NUCIMP.
2217  COMMON /nucimp/ prmom(5,248),tamom(5,248), prmfep,prmfen,tamfep,
2218  +tamfen, prefep,prefen,taefep,taefen, prepot(210),taepot(210),
2219  +prebin,taebin,fermod,etacou
2220 *KEND.
2221 C------------------------------------------------------------------
2222 C PPN=SQRT((EPN-AAM(IJPROJ))*(EPN+AAM(IJPROJ)))
2223 c ipaupr=2
2224 c ipri=3
2225  irejfo=0
2226  IF(ipri.GE.2) WRITE(6,1001) ijproj,ijproj,ppn,epn,cccxp,cccyp,
2227  +ccczp,ihtaww,itta,ieline
2228  1001 FORMAT(' HADHAD 1:',
2229  +' IJPROJ,IJPROJ,PPN,EPN,CCCXP,CCCYP,CCCZP,IHTAWW,ITTA,IELINE'/ 2
2230  +i3,2e12.3,3f7.3,3i4)
2231  cccxp=0.
2232  cccyp=0.
2233  ccczp=1.
2234  ieline=0
2235  CALL sihnin(ijproj,itta,ppn,sight)
2236  CALL sihnel(ijproj,itta,ppn,sighte)
2237  sigtot=sight + sighte
2238  IF (sigtot*rndm(bb).LE.sighte)ieline=1
2239  IF(ipri.GE.2) WRITE(6,1000) ijproj,ijproj,ppn,epn,cccxp,cccyp,
2240  +ccczp,ihtaww,itta,ieline
2241  1000 FORMAT(' HADHAD 2 nach si...:',
2242  +' IJPROJ,IJPROJ,PPN,EPN,CCCXP,CCCYP,CCCZP,IHTAWW,ITTA,IELINE'/ 2
2243  +i3,2e12.3,3f7.3,3i4)
2244  ihadha=0
2245  12 CONTINUE
2246  ihadha=ihadha+1
2247  IF(ipri.GE.2) WRITE(6,1012) ijproj,ijproj,ppn,epn,cccxp,cccyp,
2248  +ccczp,ihtaww,itta,ieline
2249  1012 FORMAT(' HADHAD 12 loop:',
2250  +' IJPROJ,IJPROJ,PPN,EPN,CCCXP,CCCYP,CCCZP,IHTAWW,ITTA,IELINE'/ 2
2251  +i3,2e12.3,3f7.3,3i4)
2252  IF(ipri.GE.3) THEN
2253  do 1212 ii=1,irh
2254  WRITE(6,'(I3,5(1PE12.4),I5/3X,5(1PE12.4))') ii,elrh(ii),plrh
2255  + (ii),cxrh(ii),cyrh(ii),czrh(ii),itrh(ii), (phkk(jjj,nhkk),jjj
2256  + =1,5)
2257  1212 continue
2258  ENDIF
2259 * repeated entry if Pauli blocking was active
2260  CALL fhad(ijproj,ijproj,ppn,epn,cccxp,cccyp,ccczp, ihtaww,itta,
2261  +ieline,irejfh)
2262  IF(irejfh.EQ.1)THEN
2263  irejfo=1
2264  IF(ipri.GE.3)
2265  + WRITE(6,'(A)')' exit from hadhad with irejfo=1 '
2266  RETURN
2267  ENDIF
2268 *
2269 * require Pauli blocking for final state nucleons
2270 *
2271  IF (ihadha.LT.3)THEN
2272  DO 11 ii=1,irh
2273  itsec=itrh(ii)
2274  IF(itsec.EQ.1.AND.elrh(ii).LE.taefep+aam(itsec)) goto 12
2275  IF(itsec.EQ.8.AND.elrh(ii).LE.taefen+aam(itsec)) goto 12
2276  IF(iibar(itsec).NE.1.AND.elrh(ii)-aam(itsec)
2277  + .LE.taepot(itsec)) goto 12
2278  11 CONTINUE
2279  ENDIF
2280  nhkkh1=nhkk
2281 C
2282  IF (ipri.GE.2) WRITE (6,1010)irh,nhkkh1,ihtaww,itta
2283  1010 FORMAT (' HADHAD IRH,NHKKH1,IHTAWW,ITTA = ',4i5)
2284 C
2285  IF(ipri.GE.3) THEN
2286  WRITE(6,'(A/5X,A)')
2287  + ' HADHAD - PARTICLE TRANSFER FROM /FINLSP/ INTO /HKKEVT/',
2288  + ' II, ELRH, PLRH, CXRH, CYRH, CZRH / PHKK(1-5)'
2289  ENDIF
2290 C
2291  isthkk(1)=11
2292  DO 10 ii=1,irh
2293 C IF( (ITSEC.EQ.1.AND.ELRH(II).GE.TAEFEP+AAM(ITSEC)) .OR.
2294 C + (ITSEC.EQ.8.AND.ELRH(II).GE.TAEFEN+AAM(ITSEC)) )THEN
2295  nhkk=nhkk+1
2296  IF (nhkk.EQ.nmxhkk)THEN
2297  WRITE (6,'(A,2I5)') .EQ.' HADHAD:NHKKNMXHKK ',nhkk,nmxhkk
2298  RETURN
2299  ENDIF
2300  itsec=itrh(ii)
2301  idhkk(nhkk)=mpdgha(itsec)
2302  jmohkk(1,nhkk)=1
2303  jmohkk(2,nhkk)=ihtaww
2304  jdahkk(1,nhkk)=0
2305  jdahkk(2,nhkk)=0
2306  phkk(1,nhkk)=plrh(ii)*cxrh(ii)
2307  phkk(2,nhkk)=plrh(ii)*cyrh(ii)
2308  phkk(3,nhkk)=plrh(ii)*czrh(ii)
2309  phkk(4,nhkk)=elrh(ii)
2310  IF(phkk(4,nhkk)-aam(itsec).LE.taepot(itsec).
2311  + and.iibar(itsec).EQ.1)THEN
2312  isthkk(nhkk)=16
2313  ELSE
2314  isthkk(nhkk)=1
2315  ENDIF
2316  phkk(5,nhkk)=aam(itrh(ii))
2317 C
2318  IF(ipri.GE.3) THEN
2319  WRITE(6,'(I3,5(1PE12.4),I5/3X,5(1PE12.4),I5)')
2320  + ii,elrh(ii),plrh
2321  + (ii),cxrh(ii),cyrh(ii),czrh(ii),itrh(ii), (phkk(jjj,nhkk),jjj
2322  + =1,5),irejfo
2323  ENDIF
2324  vhkk(1,nhkk)=vhkk(1,ihtaww)
2325  vhkk(2,nhkk)=vhkk(2,ihtaww)
2326  vhkk(3,nhkk)=vhkk(3,ihtaww)
2327  vhkk(4,nhkk)=vhkk(4,1)
2328 C ENDIF
2329  10 CONTINUE
2330  jdahkk(1,1)=nhkkh1+1
2331  jdahkk(2,1)=nhkk
2332  jdahkk(1,ihtaww)=nhkkh1+1
2333  jdahkk(2,ihtaww)=nhkk
2334 c ipaupr=0
2335 c ipri=0
2336  IF(ipri.GE.3)
2337  + WRITE(6,'(A)')' exit from hadhad with irejfo=0 '
2338  RETURN
2339  END
2340  SUBROUTINE chebch(IREJ,NHKKH1)
2341  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2342  SAVE
2343 *KEEP,HKKEVT.
2344 c INCLUDE (HKKEVT)
2345  parameter(nmxhkk= 89998)
2346 c PARAMETER (NMXHKK=25000)
2347  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
2348  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
2349  +(4,nmxhkk)
2350  COMMON /extevt/ idres(nmxhkk),idxres(nmxhkk),nobam(nmxhkk),
2351  & idbam(nmxhkk),idch(nmxhkk),npoint(10)
2352 C
2353 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
2354 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
2355 C THE POSITIONS OF THE PROJECTILE NUCLEONS
2356 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
2357 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
2358 C COMPLETELY CONSISTENT. THE TIMES IN THE
2359 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
2360 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
2361 C
2362 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
2363 C
2364 C NMXHKK: maximum numbers of entries (partons/particles) that can be
2365 C stored in the commonblock.
2366 C
2367 C NHKK: the actual number of entries stored in current event. These are
2368 C found in the first NHKK positions of the respective arrays below.
2369 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
2370 C entry.
2371 C
2372 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
2373 C = 0 : null entry.
2374 C = 1 : an existing entry, which has not decayed or fragmented.
2375 C This is the main class of entries which represents the
2376 C "final state" given by the generator.
2377 C = 2 : an entry which has decayed or fragmented and therefore
2378 C is not appearing in the final state, but is retained for
2379 C event history information.
2380 C = 3 : a documentation line, defined separately from the event
2381 C history. (incoming reacting
2382 C particles, etc.)
2383 C = 4 - 10 : undefined, but reserved for future standards.
2384 C = 11 - 20 : at the disposal of each model builder for constructs
2385 C specific to his program, but equivalent to a null line in the
2386 C context of any other program. One example is the cone defining
2387 C vector of HERWIG, another cluster or event axes of the JETSET
2388 C analysis routines.
2389 C = 21 - : at the disposal of users, in particular for event tracking
2390 C in the detector.
2391 C
2392 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
2393 C standard.
2394 C
2395 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
2396 C The value is 0 for initial entries.
2397 C
2398 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
2399 C one mother exist, in which case the value 0 is used. In cluster
2400 C fragmentation models, the two mothers would correspond to the q
2401 C and qbar which join to form a cluster. In string fragmentation,
2402 C the two mothers of a particle produced in the fragmentation would
2403 C be the two endpoints of the string (with the range in between
2404 C implied).
2405 C
2406 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
2407 C entry has not decayed, this is 0.
2408 C
2409 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
2410 C entry has not decayed, this is 0. It is assumed that the daughters
2411 C of a particle (or cluster or string) are stored sequentially, so
2412 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
2413 C daughters. Even in cases where only one daughter is defined (e.g.
2414 C K0 -> K0S) both values should be defined, to make for a uniform
2415 C approach in terms of loop constructions.
2416 C
2417 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
2418 C
2419 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
2420 C
2421 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
2422 C
2423 C PHKK(4,IHKK) : energy, in GeV.
2424 C
2425 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
2426 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
2427 C
2428 C VHKK(1,IHKK) : production vertex x position, in mm.
2429 C
2430 C VHKK(2,IHKK) : production vertex y position, in mm.
2431 C
2432 C VHKK(3,IHKK) : production vertex z position, in mm.
2433 C
2434 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
2435 C********************************************************************
2436 *KEEP,DPAR.
2437 C /DPAR/ CONTAINS PARTICLE PROPERTIES
2438 C ANAME = LITERAL NAME OF THE PARTICLE
2439 C AAM = PARTICLE MASS IN GEV
2440 C GA = DECAY WIDTH
2441 C TAU = LIFE TIME OF INSTABLE PARTICLES
2442 C IICH = ELECTRIC CHARGE OF THE PARTICLE
2443 C IIBAR = BARYON NUMBER
2444 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
2445 C
2446  CHARACTER*8 aname
2447  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
2448  +iibar(210),k1(210),k2(210)
2449 C------------------
2450 *KEEP,NUCC.
2451  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
2452 *KEEP,DPRIN.
2453  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
2454 *KEND.
2455  COMMON /chabai/chargi,barnui
2456  COMMON /evappp/ievap
2457 C-----------------------------------------------------------------------
2458  DATA ievl /0/
2459 C----------------------------------------------------------------------
2460  zero=0
2461  oneone=1
2462  twotwo=2
2463  nhad=0
2464  nip=ip
2465  aip=ip
2466  ievl=ievl+1
2467  chaeve=0.
2468  baeve=0.
2469  IF(ievap.EQ.0)THEN
2470  DO 1171 i=1,nhkkh1
2471  IF (isthkk(i).EQ.13)THEN
2472  baeve=baeve+1
2473  IF(idhkk(i).EQ.2212)chaeve=chaeve+1.d0
2474  ENDIF
2475  IF (isthkk(i).EQ.14)THEN
2476  baeve=baeve+1
2477  IF(idhkk(i).EQ.2212)chaeve=chaeve+1.d0
2478  ENDIF
2479  1171 CONTINUE
2480  DO 521 i=nhkkh1,nhkk
2481  IF (isthkk(i).EQ.1.OR.isthkk(i).EQ.15.OR.isthkk(i).EQ.16)THEN
2482  nhad=nhad+1
2483  nrhkk=mcihad(idhkk(i))
2484  IF (nrhkk.LE.0.OR.nrhkk.GT.410)THEN
2485  WRITE(6,1389)nrhkk,i,idhkk(i),nhkkh1,nhkk
2486  1389 FORMAT (' distr: NRHKK ERROR ',5i10)
2487  nrhkk=1
2488  ENDIF
2489  ichhkk=iich(nrhkk)
2490  ibhkk=iibar(nrhkk)
2491  chaeve=chaeve+ichhkk
2492  baeve=baeve+ibhkk
2493  ENDIF
2494  521 CONTINUE
2495  ELSEIF(ievap.EQ.1)THEN
2496  DO 1521 i=1,nhkk
2497  IF (isthkk(i).EQ.1)THEN
2498  nrhkk=mcihad(idhkk(i))
2499  ichhkk=iich(nrhkk)
2500  ibhkk=iibar(nrhkk)
2501  chaeve=chaeve+ichhkk
2502  baeve=baeve+ibhkk
2503  ENDIF
2504  1521 CONTINUE
2505 C WRITE(6,'(A,2F12.1)')' after isthkk=1 ',CHAEVE,BAEVE
2506  DO 2521 i=1,nhkk
2507  IF (isthkk(i).EQ.-1)THEN
2508  IF(idhkk(i).EQ.2112)THEN
2509  baeve=baeve+1
2510 C WRITE(6,'(A,2F12.1)')' evap isthkk=-1',CHAEVE,BAEVE
2511  ENDIF
2512  IF(idhkk(i).EQ.2212)THEN
2513  chaeve=chaeve+1
2514  baeve=baeve+1
2515 C WRITE(6,'(A,2F12.1)')' evap isthkk=-1',CHAEVE,BAEVE
2516  ENDIF
2517  ENDIF
2518  IF((idhkk(i).EQ.80000).AND.(isthkk(i).NE.1000))THEN
2519  chaeve=chaeve+idxres(i)
2520  baeve=baeve+idres(i)
2521 C WRITE(6,'(A,2F12.1,2I5)')' h.f.',CHAEVE,BAEVE,IDXRES(I),IDRES(I)
2522  ENDIF
2523  2521 CONTINUE
2524  ENDIF
2525  IF(ievl.LE.10)WRITE(6,'(2A,4F10.2)')' Event charge and B-number',
2526  * '=',chaeve,baeve,chargi,barnui
2527  IF(chaeve-chargi.NE.0.d0.OR.baeve-barnui.NE.0.d0)THEN
2528 C DO 775 JJJ=1,200
2529  IF(ievl.LE.1000)WRITE(6,'(2A,4F10.2)')'Event charge and B-numb',
2530  *'(violated) =',chaeve,baeve,chargi,barnui
2531 C 775 CONTINUE
2532  irej=1
2533  ENDIF
2534  RETURN
2535  END
2536 C****************************************************************
2537 C
2538  SUBROUTINE parpt(IFL,PT1,PT2,IPT,NEVT)
2539 C Plot parton pt distribution
2540 C
2541  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2542  SAVE
2543  dimension pt(50,10),ypt(50,10)
2544  go to(1,2,3),ifl
2545  1 CONTINUE
2546  dpt=0.1
2547  DO 10 i=1,10
2548  DO 10 j=1,50
2549  pt(j,i)=j*dpt-dpt/2.
2550  ypt(j,i)=1.d-50
2551  10 CONTINUE
2552  RETURN
2553  2 CONTINUE
2554  ipt1=pt1/dpt+1.
2555  ipt2=pt2/dpt+1.
2556  IF(ipt1.GT.50)ipt1=50
2557  IF(ipt2.GT.50)ipt2=50
2558  ypt(ipt1,ipt)=ypt(ipt1,ipt)+1.
2559  ypt(ipt2,ipt)=ypt(ipt2,ipt)+1.
2560  ypt(ipt1,10)=ypt(ipt1,10)+1.
2561  ypt(ipt2,10)=ypt(ipt2,10)+1.
2562  RETURN
2563  3 CONTINUE
2564  DO 30 i=1,10
2565  DO 30 j=1,50
2566  ypt(j,i)=ypt(j,i)/nevt
2567  ypt(j,i)=log10(ypt(j,i)+1.d-18)
2568  30 CONTINUE
2569 C WRITE(6,*)' Parton pt distribution,vv=1,vsr=+2,sv=3,ss=4,zz=5,
2570 C * hh=6,10=all'
2571 C CALL PLOT(PT,YPT,500,10,50,0.D0,DPT,-5.D0,0.1D0)
2572  RETURN
2573  END
2574 *
2575 *
2576 *
2577 *===hkkfil=============================================================*
2578 *
2579  SUBROUTINE hkkfil(IST,ID,M1,M2,PX,PY,PZ,E,NHKKAU,KORMO,ICALL)
2580 
2581  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2582  SAVE
2583  parameter(lout=6,llook=9)
2584  parameter(tiny10=1.0d-10,tiny4=1.0d-3)
2585 
2586 *KEEP,HKKEVT.
2587 c INCLUDE (HKKEVT)
2588  parameter(nmxhkk= 89998)
2589 c PARAMETER (NMXHKK=25000)
2590  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
2591  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
2592  +(4,nmxhkk)
2593 C
2594 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
2595 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
2596 C THE POSITIONS OF THE PROJECTILE NUCLEONS
2597 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
2598 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
2599 C COMPLETELY CONSISTENT. THE TIMES IN THE
2600 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
2601 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
2602 C
2603 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
2604 C
2605 C NMXHKK: maximum numbers of entries (partons/particles) that can be
2606 C stored in the commonblock.
2607 C
2608 C NHKK: the actual number of entries stored in current event. These are
2609 C found in the first NHKK positions of the respective arrays below.
2610 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
2611 C entry.
2612 C
2613 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
2614 C = 0 : null entry.
2615 C = 1 : an existing entry, which has not decayed or fragmented.
2616 C This is the main class of entries which represents the
2617 C "final state" given by the generator.
2618 C = 2 : an entry which has decayed or fragmented and therefore
2619 C is not appearing in the final state, but is retained for
2620 C event history information.
2621 C = 3 : a documentation line, defined separately from the event
2622 C history. (incoming reacting
2623 C particles, etc.)
2624 C = 4 - 10 : undefined, but reserved for future standards.
2625 C = 11 - 20 : at the disposal of each model builder for constructs
2626 C specific to his program, but equivalent to a null line in the
2627 C context of any other program. One example is the cone defining
2628 C vector of HERWIG, another cluster or event axes of the JETSET
2629 C analysis routines.
2630 C = 21 - : at the disposal of users, in particular for event tracking
2631 C in the detector.
2632 C
2633 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
2634 C standard.
2635 C
2636 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
2637 C The value is 0 for initial entries.
2638 C
2639 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
2640 C one mother exist, in which case the value 0 is used. In cluster
2641 C fragmentation models, the two mothers would correspond to the q
2642 C and qbar which join to form a cluster. In string fragmentation,
2643 C the two mothers of a particle produced in the fragmentation would
2644 C be the two endpoints of the string (with the range in between
2645 C implied).
2646 C
2647 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
2648 C entry has not decayed, this is 0.
2649 C
2650 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
2651 C entry has not decayed, this is 0. It is assumed that the daughters
2652 C of a particle (or cluster or string) are stored sequentially, so
2653 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
2654 C daughters. Even in cases where only one daughter is defined (e.g.
2655 C K0 -> K0S) both values should be defined, to make for a uniform
2656 C approach in terms of loop constructions.
2657 C
2658 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
2659 C
2660 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
2661 C
2662 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
2663 C
2664 C PHKK(4,IHKK) : energy, in GeV.
2665 C
2666 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
2667 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
2668 C
2669 C VHKK(1,IHKK) : production vertex x position, in mm.
2670 C
2671 C VHKK(2,IHKK) : production vertex y position, in mm.
2672 C
2673 C VHKK(3,IHKK) : production vertex z position, in mm.
2674 C
2675 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
2676 C********************************************************************
2677  COMMON /nncms/ gacms,bgcms,umo,pcm,eproj,pproj
2678  COMMON /trafop/galab,bglab,blab
2679  COMMON /projk/ iprojk
2680  COMMON /ndon/ndone
2681 
2682 C IF (MODE.GT.100) THEN
2683 C WRITE(LOUT,'(1X,A,I5,A,I5)')
2684 C & 'HKKFIL: reset NHKK = ',NHKK,' to NHKK =',NHKK-MODE+100
2685 C NHKK = NHKK-MODE+100
2686 C RETURN
2687 C ENDIF
2688  mo1 = m1
2689  mo2 = m2
2690  nhkk = nhkk+1
2691  IF (nhkk.GT.nmxhkk) THEN
2692  WRITE(lout,1000) nhkk
2693  1000 FORMAT(1x,'HKKFIL: NHKK exeeds NMXHKK = ',i7,
2694  & '! program execution stopped..')
2695  stop
2696  ENDIF
2697  IF (m1.LT.0) mo1 = nhkk+m1
2698  IF (m2.LT.0) mo2 = nhkk+m2
2699  isthkk(nhkk) = ist
2700  idhkk(nhkk) = id
2701  IF(kormo.EQ.999)THEN
2702  jmohkk(1,nhkk) = mo1
2703  jmohkk(2,nhkk) = mo2
2704  ELSE
2705  jmohkk(1,nhkk)=nhkkau+kormo-1
2706  jmohkk(2,nhkk)=0
2707  ENDIF
2708  IF(nhkk.LE.jmohkk(1,nhkk))THEN
2709 C SUBROUTINE HKKFIL(IST,ID,M1,M2,PX,PY,PZ,E,NHKKAU,KORMO,ICALL)
2710  WRITE(6,*)' HKKFIL(IST,ID,M1,M2,PX,PY,PZ,E,NHKKAU,KORMO)',
2711  * nhkk,ist,id,m1,m2,px,py,pz,e,nhkkau,kormo,icall,jmohkk(1,nhkk)
2712  ENDIF
2713  jdahkk(1,nhkk) = 0
2714  jdahkk(2,nhkk) = 0
2715  IF (mo1.GT.0) THEN
2716  IF (jdahkk(1,mo1).NE.0) THEN
2717  jdahkk(2,mo1) = nhkk
2718  ELSE
2719  jdahkk(1,mo1) = nhkk
2720  ENDIF
2721  jdahkk(1,mo1)=nhkkau
2722  ENDIF
2723  IF (mo2.GT.0) THEN
2724  IF (jdahkk(1,mo2).NE.0) THEN
2725  jdahkk(2,mo2) = nhkk
2726  ELSE
2727  jdahkk(1,mo2) = nhkk
2728  ENDIF
2729  jdahkk(1,mo2) = nhkkau
2730  ENDIF
2731  phkk(1,nhkk) = px
2732  phkk(2,nhkk) = py
2733  phkk(3,nhkk) = pz
2734  phkk(4,nhkk) = e
2735  phkk(5,nhkk) = phkk(4,nhkk)**2-phkk(1,nhkk)**2-
2736  & phkk(2,nhkk)**2-phkk(3,nhkk)**2
2737  IF ((phkk(5,nhkk).LT.0.0d0).AND.(abs(phkk(5,nhkk)).GT.tiny4))
2738  & WRITE(lout,'(1X,A,G10.3)')
2739  & 'HKKFIL: negative mass**2 ',phkk(5,nhkk)
2740  phkk(5,nhkk) = sqrt(abs(phkk(5,nhkk)))
2741  IF (ist.EQ.88888.OR.ist.EQ.88887.OR.ist.EQ.88889) THEN
2742 * special treatment for chains:
2743 * position of chain in Lab = pos. of target nucleon
2744 * time of chain-creation in Lab = time of passage of projectile
2745 * nucleus at pos. of taget nucleus
2746  DO 1 i=1,3
2747  vhkk(i,nhkk) = vhkk(i,mo2)
2748  1 CONTINUE
2749  vhkk(4,nhkk) = vhkk(3,mo2)/blab-vhkk(3,mo1)/bglab
2750  ELSE
2751  IF(mo1.GE.1)THEN
2752  DO 2 i=1,4
2753  vhkk(i,nhkk) = vhkk(i,mo1)
2754  IF (iprojk.EQ.1) THEN
2755  whkk(i,nhkk) = whkk(i,mo1)
2756  ENDIF
2757  2 CONTINUE
2758  ENDIF
2759  ENDIF
2760 
2761  RETURN
2762  END
2763 
2764 C*********************************************************************
2765 C*********************************************************************
2766 
2767  SUBROUTINE jsparr
2768  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2769  SAVE
2770  INTEGER pycomp
2771 
2772 C...Purpose: to give program heading, or list an event, or particle
2773 C...data, or current parameter values.
2774  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
2775  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
2776  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
2777  common/pydat3/mdcy(500,3),mdme(4000,2),brat(4000),kfdp(4000,5)
2778  SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/
2779  CHARACTER chap*16,chan*16,chad(5)*16
2780  dimension kbamdp(5)
2781 
2782 C...List parton/particle data table. Check whether to be listed.
2783  WRITE(mstu(11),6800)
2784  mstj24=mstj(24)
2785  mstj(24)=0
2786  kfmax=20883
2787  IF(mstu(2).NE.0) kfmax=mstu(2)
2788 C KF = PDG Particle number
2789  DO 220 kf=100,kfmax
2790 C KC = Lund particle number
2791  kc=pycomp(kf)
2792  IF(kc.EQ.0) goto 220
2793  IF(mstu(14).EQ.0.AND.kf.GT.100.AND.kc.LE.100) goto 220
2794  IF(mstu(14).GT.0.AND.kf.GT.100.AND.max(mod(kf/1000,10),
2795  & mod(kf/100,10)).GT.mstu(14)) goto 220
2796 C BAMJET particle number
2797  kbam=mcihad(kf)
2798 C BAMJET antiparticle number
2799  kabam=mcihad(-kf)
2800 
2801 C...Find Lund particle name and mass. Print information.
2802  CALL pyname(kf,chap)
2803  IF(kf.LE.100.AND.chap.EQ.' '.AND.mdcy(kc,2).EQ.0) goto 220
2804 C Lund Antiparticle Name
2805  CALL pyname(-kf,chan)
2806  pm=pymass(kf)
2807  idc1=mdcy(kc,2)
2808  idc2=mdcy(kc,2)+mdcy(kc,3)-1
2809  WRITE(mstu(11),6900)kbam,
2810  & kf,kc,idc1,idc2,chap,chan,kchg(kc,1),kchg(kc,2),
2811  & kchg(kc,3),pm,pmas(kc,2),pmas(kc,3),pmas(kc,4),mdcy(kc,1)
2812  WRITE(26,6900)kbam,
2813  & kf,kc,idc1,idc2,chap,chan,kchg(kc,1),kchg(kc,2),
2814  & kchg(kc,3),pm,pmas(kc,2),pmas(kc,3),pmas(kc,4),mdcy(kc,1)
2815 
2816 C...Particle decay: channel number, branching ration, matrix element,
2817 C...decay products.
2818  IF(kf.GT.100.AND.kc.LE.100) goto 220
2819  DO 210 idc=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1
2820  DO 200 j=1,5
2821 C Lund names of decay products
2822  CALL pyname(kfdp(idc,j),chad(j))
2823 C Bamjet numbers of decay products
2824  kbamdp(j)=mcihad(kfdp(idc,j))
2825  IF(kbamdp(j).EQ.26)kbamdp(j)=0
2826  200 CONTINUE
2827  WRITE(26,7001) idc,mdme(idc,1),mdme(idc,2),brat(idc),
2828  & (kbamdp(j),j=1,5)
2829  210 WRITE(mstu(11),7000) idc,mdme(idc,1),mdme(idc,2),brat(idc),
2830  & (chad(j),j=1,5)
2831 C The same for the antiparticle, if it exists
2832  IF(kabam.NE.410)THEN
2833  WRITE(mstu(11),6900)kabam,
2834  & -kf,-kc,idc1,idc2,chan,chap,-kchg(kc,1),kchg(kc,2),
2835  & kchg(kc,3),pm,pmas(kc,2),pmas(kc,3),pmas(kc,4),mdcy(kc,1)
2836  WRITE(26,6900)kabam,
2837  & -kf,-kc,idc1,idc2,chan,chap,-kchg(kc,1),kchg(kc,2),
2838  & kchg(kc,3),pm,pmas(kc,2),pmas(kc,3),pmas(kc,4),mdcy(kc,1)
2839  DO 211 idc=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1
2840  DO 201 j=1,5
2841 C KC = Lund particle number
2842  kcdp=pycomp(kfdp(idc,j))
2843  IF(kcdp.LE.0.OR.kcdp.GT.500)THEN
2844 C WRITE(MSTU(11),'(A,I10)')' KCDP= ',KCDP
2845  kcdp=1
2846  ENDIF
2847 C Bamjet numbers of decay products
2848  kfdpm=-kfdp(idc,j)
2849  IF(kchg(kcdp,3).EQ.0)kfdpm=kfdp(idc,j)
2850  kbamdp(j)=mcihad(kfdpm)
2851  IF(kbamdp(j).EQ.26)kbamdp(j)=0
2852 C Lund names of decay products
2853  CALL pyname(kfdpm,chad(j))
2854  201 CONTINUE
2855  WRITE(26,7001) idc,mdme(idc,1),mdme(idc,2),brat(idc),
2856  & (kbamdp(j),j=1,5)
2857  211 WRITE(mstu(11),7000) idc,mdme(idc,1),mdme(idc,2),brat(idc),
2858  & (chad(j),j=1,5)
2859  ENDIF
2860  220 CONTINUE
2861  mstj(24)=mstj24
2862 
2863 
2864 C...Format statements for output on unit MSTU(11) (by default 6).
2865  6800 FORMAT(///30x,'Particle/parton data table'//1x,'BAM',
2866  &1x,'ABAM',1x,'KF',1x,'KC',1x,'DCF',1x,'DCL',1x,
2867  &'particle',8x,'antiparticle',6x,'chg col anti',8x,'mass',7x,
2868  &'width',7x,'w-cut',5x,'lifetime',1x,'decay'/11x,'IDC',1x,'on/off',
2869  &1x,'ME',3x,'Br.rat.',4x,'decay products')
2870  6900 FORMAT(/1x,i4,i6,i4,2i5,a16,a16,3i3,1x,f12.5,2(1x,f11.5),
2871  &1x,f12.5,1x,i2)
2872  7000 FORMAT(10x,i4,2x,i3,2x,i3,2x,f8.5,4x,5a16)
2873  7001 FORMAT(10x,i4,2x,i3,2x,i3,2x,f8.5,4x,5i5)
2874 
2875  RETURN
2876  END
2877 
2878 C*********************************************************************
2879 
Generated on Sat May 25 2013 14:32:29 for Geant4 by   doxygen 1.8.4