Extra header added to the list
[u/mrichter/AliRoot.git] / DPMJET / phojet1.12-35c4.f
1 C***********************************************************************
2 C
3 C
4 C
5 C                       PHOJET version 1.12
6 C                       -------------------
7 C
8 C
9 C    ($Revision: 1.12.1.35 $, $Date: 2000/06/25 21:59:19 $)
10 C
11 C
12 C    Authors: Ralph Engel
13 C             (ralph.engel@fzk.de)
14 C
15 C             Johannes Ranft
16 C             (johannes.ranft@cern.ch)
17 C
18 C             Stefan Roesler
19 C             (Stefan.Roesler@cern.ch)
20 C
21 C
22 C    For the latest version and documentation check
23 C       http://www-ik.fzk.de/~engel/phojet.html
24 C
25 C
26 C    Bug reports, questions, complaints are welcome
27 C    (please send a mail to ralph.engel@fzk.de).
28 C
29 C
30 C    Note that the code is available with several interfaces to
31 C    Lund fragmentation programs (JETSET7.x, 1.x and a double
32 C    precision JETSET version). This file is the code with
33 C
34
35 C                interface to PYTHIA 6.1 (or higher)
36
37 C     for usage in DPMJET 3.x (Lund common block dimensions increased)
38
39 C
40 C***********************************************************************
41 C
42 C
43 C             List of subroutines and functions
44 C             ---------------------------------
45 C
46 C
47 C  main event simulation routines
48 C
49 C      PHO_EVENT
50 C      PHO_PARTON
51 C      PHO_POSPOM
52 C
53 C      PHO_STDPAR
54 C      PHO_POMSCA
55 C
56 C
57 C  user steering interface
58 C
59 C      PHO_SETMDL
60 C      PHO_PRESEL
61 C
62 C
63 C  experimental setup / photon flux calculation
64 C
65 C      PHO_FIXLAB
66 C      PHO_FIXCOL
67 C      PHO_GPHERA
68 C      PHO_GGEPEM
69 C      PHO_WGEPEM
70 C      PHO_GGBLSR
71 C      PHO_GGBEAM
72 C      PHO_GGHIOF
73 C      PHO_GGHIOG
74 C      PHO_GGFLCL
75 C      PHO_GGFLCR
76 C      PHO_GGFAUX
77 C      PHO_GGFNUC
78 C      PHO_GHHIOF
79 C      PHO_GHHIAS
80 C
81 C
82 C  initialization
83 C
84 C      PHO_INIT
85 C      PHO_DATINI
86 C      PHO_PARDAT
87 C      PHO_MCINI
88 C
89 C      PHO_EVEINI
90 C
91 C      PHO_HARINI
92 C      PHO_FRAINI
93 C
94 C      PHO_FITPAR
95 C
96 C
97 C  cross section calculation
98 C
99 C      PHO_CSINT
100 C
101 C      PHO_XSECT
102 C      PHO_BORNCS
103 C      PHO_HARXTO
104 C
105 C      PHO_DSIGDT
106 C
107 C      PHO_TRIREG
108 C      PHO_LOOREG
109 C      PHO_TRXPOM
110 C
111 C      PHO_EIKON
112 C      PHO_CHAN2A
113 C
114 C      PHO_SCALES
115 C
116 C
117 C  multiple interaction structure
118 C
119 C      PHO_IMPAMP
120 C      PHO_PRBDIS
121 C      PHO_SAMPRO
122 C      PHO_SAMPRB
123 C
124 C
125 C  hadron / photon remnant treatment, soft x selection
126 C
127 C      PHO_HARREM
128 C      PHO_PARREM
129 C
130 C      PHO_HADSP2
131 C      PHO_HADSP3
132 C      PHO_SOFTXX
133 C      PHO_SELSXR
134 C      PHO_SELSX2
135 C      PHO_SELSXS
136 C      PHO_SELSXI
137 C
138 C      PHO_VALFLA
139 C      PHO_REGFLA
140 C      PHO_SEAFLA
141 C      PHO_FLAUX
142 C      PHO_BETAF
143 C      IPHO_DIQU
144 C
145 C
146 C  primordial kt and soft parton pt
147 C
148 C      PHO_PRIMKT
149 C      PHO_PARTPT
150 C      PHO_SOFTPT
151 C      PHO_SELPT
152 C
153 C      PHO_CONN0
154 C      PHO_CONN1
155 C
156 C
157 C  simulation of hard scattering, initial state radiation
158 C
159 C      PHO_HARCOL
160 C      PHO_SELCOL
161 C      PHO_HARCOR
162 C
163 C      PHO_HARDIR
164 C      PHO_HARX12
165 C      PHO_HARDX1
166 C      PHO_HARKIN
167 C      PHO_HARWGH
168 C      PHO_HARSCA
169 C      PHO_HARFAC
170 C      PHO_HARWGX
171 C      PHO_HARWGI
172 C      PHO_HARINT
173 C      PHO_HARMCI
174 C
175 C      PHO_HARXR3
176 C      PHO_HARXR2
177 C      PHO_HARXD2
178 C      PHO_HARXPT
179 C      PHO_HARISR
180 C      PHO_HARZSP
181 C
182 C      PHO_PTCUT
183 C      PHO_ALPHAE
184 C      PHO_ALPHAS
185 C
186 C
187 C  diffraction dissociation
188 C
189 C      PHO_DIFDIS
190 C      PHO_DIFPRO
191 C      PHO_DIFPAR
192 C      PHO_QELAST
193 C      PHO_CDIFF
194 C      PHO_DFWRAP
195 C
196 C      PHO_SAMASS
197 C      PHO_DSIGDM
198 C      PHO_DFMASS
199 C
200 C      PHO_SDECAY
201 C      PHO_SDECY2
202 C      PHO_SDECY3
203 C
204 C      PHO_DIFSLP
205 C      PHO_DIFKIN
206 C      PHO_VECRES
207 C      PHO_DIFRES
208 C
209 C      PHO_REGPAR
210 C
211 C      PHO_PECMS
212 C      PHO_SETPAR
213 C
214 C
215 C  fragmentation, treatment of low-mass strings
216 C
217 C      PHO_STRING
218 C      PHO_STRFRA
219 C
220 C      PHO_ID2STR
221 C      PHO_MCHECK
222 C      PHO_POMCOR
223 C      PHO_MASCOR
224 C      PHO_PARCOR
225 C
226 C      PHO_GLU2QU
227 C      PHO_GLUSPL
228 C
229 C      PHO_DQMASS
230 C      PHO_BAMASS
231 C      PHO_MEMASS
232 C
233 C
234 C  particle code tables, particle numbering conversion
235 C
236 C      PHO_PNAME
237 C      PHO_PMASS
238 C      IPHO_CHR3
239 C      IPHO_BAR3
240 C
241 C      IPHO_ANTI
242 C
243 C      IPHO_PDG2ID
244 C      IPHO_ID2PDG
245 C      IPHO_LU2PDG
246 C      IPHO_PDG2LU
247 C
248 C      IPHO_CNV1
249 C      PHO_HACODE
250 C
251 C
252 C
253 C  Lorentz transformations, rotations and mass adjustment
254 C
255 C      PHO_ALTRA
256 C      PHO_LTRANS
257 C      PHO_TRANS
258 C      PHO_TRANI
259 C
260 C      PHO_MKSLTR
261 C      PHO_GETLTR
262 C
263 C      PHO_LTRHEP
264 C
265 C      PHO_MSHELL
266 C      PHO_MASSAD
267 C
268 C
269 C  program debugging and internal cross-checks
270 C
271 C      PHO_PREVNT
272 C      PHO_PRSTRG
273 C      PHO_CHECK
274 C
275 C      PHO_TRACE
276 C
277 C      PHO_REJSTA
278 C
279 C      PHO_ABORT
280 C
281 C
282 C  cross section fitting
283 C
284 C      PHO_FITMAI
285 C      PHO_FITINP
286 C      PHO_FITDAT
287 C      PHO_FITOUT
288 C      PHO_FITAMP
289 C      PHO_FITTST
290 C      PHO_FITMSQ
291 C      PHO_FITVD1
292 C      PHO_FITCN1
293 C      PHO_FITINI
294 C
295 C
296 C  cross section parametrizations
297 C
298 C      PHO_HADCSL
299 C      PHO_ALLM97
300 C      PHO_CSDIFF
301 C
302
303 C
304 C  random numbers
305 C
306
307 C      DPMJET random number generator DT_RNDM used
308
309 C
310 C      PHO_SFECFE
311 C      PHO_RNDBET
312 C      PHO_RNDGAM
313 C
314 C
315 C  auxiliary routines / numerical methods
316 C
317 C      PHO_GAUSET
318 C      PHO_GAUDAT
319 C
320 C      pho_samp1d
321 C
322 C      PHO_DZEROX
323 C      PHO_EXPINT
324 C      PHO_BESSJ0
325 C      PHO_BESSI0
326 C      pho_ExpBessI0
327 C      PHO_BESSI1
328 C      PHO_BESSK0
329 C      PHO_BESSK1
330 C
331 C      PHO_XLAM
332 C
333 C      PHO_SWAPD
334 C      PHO_SWAPI
335 C
336 C
337 C  parton density parametrization management / interface
338 C
339 C      PHO_PDF
340 C
341 C      PHO_SETPDF
342 C      PHO_GETPDF
343 C      PHO_ACTPDF
344 C
345 C      PHO_QPMPDF
346 C
347 C      PHO_PDFTST
348 C
349 C
350 C  parton density parametrizations from other authors
351 C
352 C      PHO_DOR98LO
353 C      PHO_DOR98SC
354 C      PHO_DOR94LO
355 C      PHO_DOR94HO
356 C      PHO_DOR94DI
357 C      PHO_DOR92LO
358 C      PHO_DOR92HO
359 C      PHO_DORPLO
360 C      PHO_DORPHO
361 C      PHO_DORGLO
362 C      PHO_DORGHO
363 C      PHO_DORGH0
364 C      PHO_DOR94FV
365 C      PHO_DOR94FW
366 C      PHO_DOR94FS
367 C      PHO_DOR92FV
368 C      PHO_DOR92FW
369 C      PHO_DOR92FS
370 C      PHO_DORFVP
371 C      PHO_DORFGP
372 C      PHO_DORFQP
373 C      PHO_DORGF
374 C      PHO_DORGFS
375 C      PHO_grsf1
376 C      PHO_grsf2
377 C
378 C      PHO_CKMTPA
379 C      PHO_CKMTPD
380 C      PHO_CKMTPO
381 C      PHO_CKMTFV
382 C
383 C      PHO_DBFINT
384 C
385 C      PHO_SASGAM
386 C      PHO_SASVMD
387 C      PHO_SASANO
388 C      PHO_SASBEH
389 C      PHO_SASDIR
390 C
391 C      PHO_PHGAL
392 C      PHVAL
393 C
394 C
395 C***********************************************************************
396
397 *$ CREATE PHO_INIT.FOR
398 *COPY PHO_INIT
399 CDECK  ID>, PHO_INIT
400       SUBROUTINE PHO_INIT(LINP,LOUT,IREJ)
401 C***********************************************************************
402 C
403 C     main subroutine to configure and manage PHOJET calculations
404 C
405 C     input:  LINP       input unit to read from
406 C                        -1 to skip reading of input file
407 C             LOUT       output unit to write to
408 C
409 C     output: IREJ       0  success
410 C                        1  failure
411 C
412 C***********************************************************************
413       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
414       SAVE
415
416 C  input/output channels
417       INTEGER LI,LO
418       COMMON /POINOU/ LI,LO
419 C  event debugging information
420       INTEGER NMAXD
421       PARAMETER (NMAXD=100)
422       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
423      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
424       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
425      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
426 C  model switches and parameters
427       CHARACTER*8 MDLNA
428       INTEGER ISWMDL,IPAMDL
429       DOUBLE PRECISION PARMDL
430       COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
431 C  general process information
432       INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
433       COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
434
435 C  global event kinematics and particle IDs
436       INTEGER IFPAP,IFPAB
437       DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
438       COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
439 C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
440       INTEGER IDEQP,IDEQB,IHFLD,IHFLS
441       DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
442       COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
443      &                IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
444 C  integration precision for hard cross sections (obsolete)
445       INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
446       COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
447 C  some hadron information, will be deleted in future versions
448       INTEGER NFS
449       DOUBLE PRECISION QMASS,BET,PCOUDI,PNORM,VALPRG
450       COMMON /POHDRN/ QMASS(6),BET,PCOUDI,PNORM,VALPRG(2),NFS
451 C  obsolete cut-off information
452       DOUBLE PRECISION PTCUT,PTANO,FPS,FPH,PSOMIN,XSOMIN
453       COMMON /POCUT1/ PTCUT(4),PTANO(4),FPS(4),FPH(4),PSOMIN,XSOMIN
454 C  photon flux kinematics and cuts
455       DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
456      &                 YMIN1,YMAX1,YMIN2,YMAX2,
457      &                 Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
458      &                 THMIN1,THMAX1,THMIN2,THMAX2
459       INTEGER          ITAG1,ITAG2
460       COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
461      &                YMIN1,YMAX1,YMIN2,YMAX2,
462      &                Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
463      &                THMIN1,THMAX1,THMIN2,THMAX2,
464      &                ITAG1,ITAG2
465 C  cut probability distribution
466       INTEGER IEETA1,IIMAX,KKMAX
467       PARAMETER( IEETA1=20, IIMAX=20, KKMAX=20 )
468       INTEGER IEEMAX,IMAX,KMAX
469       REAL PROB
470       DOUBLE PRECISION EPTAB
471       COMMON /POPROB/ PROB(4,IEETA1,0:IIMAX,0:KKMAX),EPTAB(4,IEETA1),
472      &                IEEMAX,IMAX,KMAX
473 C  event weights and generated cross section
474       INTEGER IPOWGC,ISWCUT,IVWGHT
475       DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
476       COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
477      &                IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
478 C  names of hard scattering processes
479       INTEGER Max_pro_1
480       PARAMETER ( Max_pro_1 = 16 )
481       CHARACTER*18 PROC
482       COMMON /POHPRO/ PROC(0:Max_pro_1)
483 C  hard cross sections and MC selection weights
484       INTEGER Max_pro_2
485       PARAMETER ( Max_pro_2 = 16 )
486       INTEGER IHa_last,IHb_last,MH_pro_on,MH_tried,
487      &  MH_acc_1,MH_acc_2
488       DOUBLE PRECISION Hfac,HWgx,HSig,Hdpt,HEcm_last,HQ2a_last,HQ2b_last
489       COMMON /POHRCS/ Hfac(-1:Max_pro_2),HWgx(-1:Max_pro_2),
490      &  HSig(-1:Max_pro_2),Hdpt(-1:Max_pro_2),
491      &  HEcm_last,HQ2a_last,HQ2b_last,IHa_last,IHb_last,
492      &  MH_pro_on(-1:Max_pro_2,0:4),MH_tried(-1:Max_pro_2,0:4),
493      &  MH_acc_1(-1:Max_pro_2,0:4),MH_acc_2(-1:Max_pro_2,0:4)
494
495       INTEGER MSTU,MSTJ
496       DOUBLE PRECISION PARU,PARJ
497       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
498
499       INTEGER KCHG
500       DOUBLE PRECISION  PMAS,PARF,VCKM
501       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
502
503       INTEGER MDCY,MDME,KFDP
504       DOUBLE PRECISION  BRAT
505       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
506
507       INTEGER PYCOMP
508
509       DIMENSION ITMP(0:11)
510       CHARACTER*10 CNAME
511       CHARACTER*70 NUMBER,FILENA
512
513  14   FORMAT(A10,A69)
514  15   FORMAT(A12)
515
516 C  define input/output units
517       IF(LINP.GE.0) THEN
518         LI = LINP
519       ELSE
520         LI = 5
521       ENDIF
522       LO = LOUT
523
524       IREJ = 0
525
526       WRITE(LO,*)
527       WRITE(LO,*) ' ==================================================='
528       WRITE(LO,*) '                                                    '
529       WRITE(LO,*) '      ----      PHOJET version 1.12      ----      '
530       WRITE(LO,*) '                                                    '
531       WRITE(LO,*) ' ==================================================='
532       WRITE(LO,*) '     Authors: Ralph Engel      (FZ Karlsruhe)'
533       WRITE(LO,*) '              Johannes Ranft   (Siegen Univ.)'
534       WRITE(LO,*) '              Stefan Roesler   (CERN)'
535       WRITE(LO,*) ' ---------------------------------------------------'
536       WRITE(LO,*) '   Manual, updates, and further information:'
537       WRITE(LO,*) '    http://www-ik.fzk.de/~engel/phojet.html'
538       WRITE(LO,*) ' ---------------------------------------------------'
539       WRITE(LO,*) '    please send suggestions / bug reports etc. to:'
540       WRITE(LO,*) '             ralph.engel@fzk.de'
541       WRITE(LO,*) ' ==================================================='
542       WRITE(LO,*) '   $Date: 2000/06/25 21:59:19 $'
543       WRITE(LO,*) '   $Revision: 1.12.1.35 $'
544
545       WRITE(LO,*) '   (code version with interface to PYTHIA 6.x)'
546
547       WRITE(LO,*) '   (code version for usage in DPMJET 3.x)'
548
549       WRITE(LO,*) ' ==================================================='
550       WRITE(LO,*)
551
552 C  standard initializations
553       CALL PHO_DATINI
554       CALL PHO_PARDAT
555       DUM = PHO_PMASS(0,-1)
556
557 C  initialize standard PDFs
558 C  proton
559       CALL PHO_SETPDF(2212,IDUM,5,6,0,0,-1)
560       CALL PHO_SETPDF(-2212,IDUM,5,6,0,0,-1)
561 C  neutron
562       CALL PHO_SETPDF(2112,IDUM,5,6,0,0,-1)
563       CALL PHO_SETPDF(-2112,IDUM,5,6,0,0,-1)
564 C  photon
565       CALL PHO_SETPDF(22,IDUM,5,3,0,0,-1)
566 C  pomeron
567       CALL PHO_SETPDF(990,IDUM,4,0,0,0,-1)
568 C  pions
569       CALL PHO_SETPDF(211,IDUM,5,2,0,0,-1)
570       CALL PHO_SETPDF(-211,IDUM,5,2,0,0,-1)
571       CALL PHO_SETPDF(111,IDUM,5,2,0,0,-1)
572 C  kaons
573       CALL PHO_SETPDF(321,IDUM,5,2,0,0,-1)
574       CALL PHO_SETPDF(-321,IDUM,5,2,0,0,-1)
575       CALL PHO_SETPDF(130,IDUM,5,2,0,0,-1)
576       CALL PHO_SETPDF(310,IDUM,5,2,0,0,-1)
577
578 C  nothing to be done
579       IF(LINP.LT.0) RETURN
580
581 C  main loop to read input cards
582  1200 CONTINUE
583         READ(LINP,14,END=1300) CNAME,NUMBER
584         IF(CNAME.EQ.'ENDINPUT  ') THEN
585           GOTO 1300
586         ELSE IF(CNAME.EQ.'STOP      ') THEN
587           WRITE(LO,*) 'STOP'
588           STOP
589         ELSE IF(CNAME.EQ.'COMMENT   ') THEN
590           WRITE(LO,'(1X,A10,A69)') 'COMMENT   ',NUMBER
591         ELSE IF(CNAME(1:1).EQ.'*') THEN
592           WRITE(LO,'(1X,A10,A69)') CNAME,NUMBER
593         ELSE IF(CNAME.EQ.'PTCUT     ') THEN
594           READ(NUMBER,*) PARMDL(36),PARMDL(37),PARMDL(38),PARMDL(39)
595           WRITE(LO,*) 'PTCUT     ',PARMDL(36),PARMDL(37),
596      &      PARMDL(38),PARMDL(39)
597         ELSE IF(CNAME.EQ.'PROCESS   ') THEN
598           READ(NUMBER,*) (IPRON(KK,1),KK=1,8)
599           WRITE(LO,*) 'PROCESS   ',(IPRON(KK,1),KK=1,8)
600         ELSE IF(CNAME.EQ.'DIFF-PROC ') THEN
601           READ(NUMBER,*) (ITMP(KK),KK=0,11)
602           WRITE(LO,*) 'DIFF-PROC ',(ITMP(KK),KK=0,8)
603           DO 112 KK=1,8
604             IPRON(KK,ITMP(0)) = ITMP(KK)
605  112      CONTINUE
606         ELSE IF(CNAME.EQ.'SUBPROCESS') THEN
607           READ(NUMBER,*) IMPRO,IP,ION
608           WRITE(LO,*) 'SUBPROCESS',IMPRO,IP,ION
609           MH_pro_on(IMPRO,IP) = ION
610         ELSE IF(CNAME.EQ.'PARTICLE1 ') THEN
611           READ(NUMBER,*) IDPDG,PVIR
612           IHFLS(1) = 1
613           XPSUB = 1.D0
614           CALL PHO_SETPAR(1,IDPDG,0,PVIR)
615           WRITE(LO,*) 'PARTICLE1  ',IDPDG,PVIR
616         ELSE IF(CNAME.EQ.'PARTICLE2 ') THEN
617           READ(NUMBER,*) IDPDG,PVIR
618           IHFLS(2) = 1
619           XTSUB = 1.D0
620           CALL PHO_SETPAR(2,IDPDG,0,PVIR)
621           WRITE(LO,*) 'PARTICLE2  ',IDPDG,PVIR
622         ELSE IF(CNAME.EQ.'REMNANT1  ') THEN
623           READ(NUMBER,*) IDPDG,IFL1,IFL2,IVAL,XSUB
624           IHFLS(1) = IVAL
625           IHFLD(1,1) = IFL1
626           IHFLD(1,2) = IFL2
627           XPSUB = XSUB
628           PVIR = 0.D0
629           CALL PHO_SETPAR(1,IDPDG,-1,PVIR)
630           WRITE(LO,*) 'REMNANT1   ',IDPDG,IFL1,IFL2,IVAL,XSUB
631         ELSE IF(CNAME.EQ.'REMNANT2  ') THEN
632           READ(NUMBER,*) IDPDG,IFL1,IFL2,IVAL,XSUB
633           IHFLS(2) = IVAL
634           IHFLD(2,1) = IFL1
635           IHFLD(2,2) = IFL2
636           XTSUB = XSUB
637           PVIR = 0.D0
638           CALL PHO_SETPAR(2,IDPDG,-1,PVIR)
639           WRITE(LO,*) 'REMNANT2   ',IDPDG,IFL1,IFL2,IVAL,XSUB
640         ELSE IF(CNAME.EQ.'PDF       ') THEN
641           READ(NUMBER,*) IDPDG,IPAR,ISET,IEXT
642           WRITE(LO,*) 'PDF        ',IDPDG,IPAR,ISET,IEXT
643           CALL PHO_SETPDF(IDPDG,IDUM,IPAR,ISET,IEXT,0,-1)
644         ELSE IF(CNAME.EQ.'SETMODEL  ') THEN
645           READ(NUMBER,*) I,IVAL
646           WRITE(LO,*) 'SETMODEL   ',I,IVAL
647           CALL PHO_SETMDL(I,IVAL,1)
648         ELSE IF(CNAME.EQ.'SETPARAM  ') THEN
649           READ(NUMBER,*) I,PARNEW
650           WRITE(LO,*) 'SETPARAM   ',I,PARNEW
651           PARMDL(I) = PARNEW
652         ELSE IF(CNAME.EQ.'DEBUG     ') THEN
653           READ(NUMBER,*) IDEBF,IDEBN,IDLEV
654           WRITE(LO,*) 'DEBUG      ',IDEBF,IDEBN,IDLEV
655           CALL PHO_TRACE(IDEBF,IDEBN,IDLEV)
656         ELSE IF(CNAME.EQ.'TRACE     ') THEN
657           READ(NUMBER,*) IDEBF,IDLEV
658           WRITE(LO,*) 'TRACE      ',IDEBF,IDLEV
659           IDEB(IDEBF) = IDLEV
660         ELSE IF(CNAME.EQ.'SETICUT   ') THEN
661           READ(NUMBER,*) I,ICUT
662           WRITE(LO,*) 'SETICUT    ',I,ICUT
663           ISWCUT(I) = ICUT
664         ELSE IF(CNAME.EQ.'SETFCUT   ') THEN
665           READ(NUMBER,*) I,PARNEW
666           WRITE(LO,*) 'SETFCUT    ',I,PARNEW
667           HSWCUT(I) = PARNEW
668         ELSE IF(CNAME.EQ.'LUND-MSTU ') THEN
669           READ(NUMBER,*) I,IVAL
670           WRITE(LO,*) 'LUND-MSTU  ',I,IVAL
671           MSTU(I) = IVAL
672         ELSE IF(CNAME.EQ.'LUND-MSTJ ') THEN
673           READ(NUMBER,*) I,IVAL
674           WRITE(LO,*) 'LUND-MSTJ  ',I,IVAL
675           MSTJ(I) = IVAL
676         ELSE IF(CNAME.EQ.'LUND-PARJ ') THEN
677           READ(NUMBER,*) I,EE
678           WRITE(LO,*) 'LUND-PARJ  ',I,EE
679           PARJ(I) = REAL(EE)
680         ELSE IF(CNAME.EQ.'LUND-PARU ') THEN
681           READ(NUMBER,*) I,EE
682           WRITE(LO,*) 'LUND-PARU  ',I,EE
683           PARU(I) = REAL(EE)
684         ELSE IF(CNAME.EQ.'LUND-DECAY') THEN
685           READ(NUMBER,*) ID,ION
686           WRITE(LO,*) 'LUND-DECAY ',ID,ION
687
688           KC=PYCOMP(ID)
689
690           MDCY(KC,1) = ION
691         ELSE IF(CNAME.EQ.'PSOFTMIN  ') THEN
692           READ(NUMBER,*) PSOMIN
693           WRITE(LO,*) 'PSOFTMIN   ',PSOMIN
694         ELSE IF(CNAME.EQ.'INTPREC   ') THEN
695           READ(NUMBER,*) NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
696           WRITE(LO,*) 'INTPREC    ',NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
697
698 C  PDF test utility
699         ELSE IF(CNAME.EQ.'PDFTEST   ') THEN
700           READ(NUMBER,*) IDPDG,SCALE2,PVIRT2
701           PVIRT2 = ABS(PVIRT2)
702           WRITE(LO,*) 'PDFTEST   ',IDPDG,' ',SCALE2,' ',PVIRT2
703           CALL PHO_PDFTST(IDPDG,SCALE2,PVIRT2)
704
705 C  mass cut on gamma-gamma or gamma-hadron system
706         ELSE IF(CNAME.EQ.'ECMS-CUT  ') THEN
707           READ(NUMBER,*) ECMIN,ECMAX
708           WRITE(LO,*) 'ECMS-CUT  ',ECMIN,ECMAX
709
710 C  beam lepton (anti-)tagging system
711         ELSE IF(CNAME.EQ.'TAG-METHOD') THEN
712           READ(NUMBER,*) ITAG1,ITAG2
713           WRITE(LO,*) 'TAG-METHOD',ITAG1,ITAG2
714         ELSE IF(CNAME.EQ.'E-TAG1    ') THEN
715           READ(NUMBER,*)
716      &      EEMIN1,YMIN1,YMAX1,Q2MIN1,Q2MAX1,THMIN1,THMAX1
717           WRITE(LO,*) 'E-TAG1    ',EEMIN1,YMIN1,YMAX1,
718      &      Q2MIN1,Q2MAX1,THMIN1,THMAX1
719         ELSE IF(CNAME.EQ.'E-TAG2    ') THEN
720           READ(NUMBER,*)
721      &      EEMIN2,YMIN2,YMAX2,Q2MIN2,Q2MAX2,THMIN2,THMAX2
722           WRITE(LO,*) 'E-TAG2    ',EEMIN2,YMIN2,YMAX2,
723      &      Q2MIN2,Q2MAX2,THMIN2,THMAX2
724
725 C  sampling of gamma-p events in ep (HERA)
726         ELSE IF(    (CNAME.EQ.'WW-HERA   ')
727      &          .OR.(CNAME.EQ.'GP-HERA   ')) THEN
728           READ(NUMBER,*) EE1,EE2,NEV
729           WRITE(LO,*) 'GP-HERA   ',EE1,EE2,NEV
730           IF(YMAX2.LT.0.D0) THEN
731             WRITE(LO,*) ' PHO_INIT:ERROR:ELECTRON TAGGER NOT SET'
732           ELSE
733             CALL PHO_GPHERA(NEV,EE1,EE2)
734             KEVENT = 0
735           ENDIF
736
737 C  sampling of gamma-gamma events in e+e- (LEP)
738         ELSE IF(    (CNAME.EQ.'GG-EPEM   ')
739      &          .OR.(CNAME.EQ.'WW-EPEM   ')) THEN
740           READ(NUMBER,*) EE1,EE2,NEV
741           WRITE(LO,*) 'GG-EPEM   ',EE1,EE2,NEV
742           IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
743             WRITE(LO,*) ' PHO_INIT:ERROR:ELECTRON TAGGERS NOT SET'
744           ELSE
745             CALL PHO_GGEPEM(-1,EE1,EE2)
746             CALL PHO_GGEPEM(NEV,EE1,EE2)
747             CALL PHO_GGEPEM(-2,sig_tot,sig_gg)
748             KEVENT = 0
749           ENDIF
750
751 C  sampling of gamma-gamma in heavy-ion collisions
752         ELSE IF(CNAME.EQ.'GG-HION-F ') THEN
753           READ(NUMBER,*) EE,NA,NZ,NEV
754           WRITE(LO,*) 'GG-HION-F ',EE,NA,NZ,NEV
755           IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
756             WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
757           ELSE
758             CALL PHO_GGHIOF(NEV,EE,NA,NZ)
759             KEVENT = 0
760           ENDIF
761         ELSE IF(CNAME.EQ.'GG-HION-G ') THEN
762           READ(NUMBER,*) EE,NA,NZ,NEV
763           WRITE(LO,*) 'GG-HION-G ',EE,NA,NZ,NEV
764           IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
765             WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
766           ELSE
767             CALL PHO_GGHIOG(NEV,EE,NA,NZ)
768             KEVENT = 0
769           ENDIF
770
771 C  sampling of gamma-hadron events in heavy ion collisions
772         ELSE IF(CNAME.EQ.'GH-HION-F ') THEN
773           READ(NUMBER,*) EE,NA,NZ,NEV
774           WRITE(LO,*) 'GH-HION-F ',EE,NA,NZ,NEV
775           IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
776             WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
777           ELSE
778             CALL PHO_GHHIOF(NEV,EE,NA,NZ)
779             KEVENT = 0
780           ENDIF
781
782 C  sampling of hadron-gamma events in hadron - heavy ion collisions
783         ELSE IF(CNAME.EQ.'HG-HIAS-F ') THEN
784           READ(NUMBER,*) EP,EE,NA,NZ,NEV
785           WRITE(LO,*) 'HG-HIAS-F ',EP,EE,NA,NZ,NEV
786           IF(YMAX2.LT.0.D0) THEN
787             WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
788           ELSE
789             CALL PHO_GHHIAS(NEV,EP,EE,NA,NZ)
790             KEVENT = 0
791           ENDIF
792
793 C  sampling of photoproduction events e+e-, backscattered laser
794         ELSE IF(CNAME.EQ.'BLASER    ') THEN
795           READ(NUMBER,*) EE1,EE2,Pl_lam_1,Pl_lam_2,X_1,X_2,rho,A,NEV
796           WRITE(LO,*) 'BLASER    ',EE1,EE2,
797      &      Pl_lam_1,Pl_lam_2,X_1,X_2,rho,A,NEV
798           CALL PHO_GGBLSR(NEV,EE1,EE2,Pl_lam_1,Pl_lam_2,X_1,X_2,rho,A)
799           KEVENT = 0
800
801 C  sampling of photoproduction events beamstrahlung
802         ELSE IF(CNAME.EQ.'BEAMST    ') THEN
803           READ(NUMBER,*) EE1,YPSI,SIGX,SIGY,SIGZ,AEB,NEV
804           WRITE(LO,*) 'BEAMST    ',EE1,YPSI,SIGX,SIGY,SIGZ,AEB,NEV
805           IF(YMAX1.LT.0.D0) THEN
806             WRITE(LO,*) ' PHO_INIT:ERROR:ELECTRON TAGGER 1 NOT SET'
807           ELSE
808             CALL PHO_GGBEAM(NEV,EE1,YPSI,SIGX,SIGY,SIGZ,AEB)
809             KEVENT = 0
810           ENDIF
811
812 C  fixed-energy events in LAB system of particle 2
813         ELSE IF(CNAME.EQ.'EVENT-LAB ') THEN
814           READ(NUMBER,*) PLAB,NEV
815           WRITE(LO,*) 'EVENT-LAB ',PLAB,NEV
816           CALL PHO_FIXLAB(PLAB,NEV)
817           KEVENT = 0
818
819 C  fixed-energy events in CM system
820         ELSE IF(CNAME.EQ.'EVENT-CMS ') THEN
821           READ(NUMBER,*) ECM,NEV
822           WRITE(LO,*) 'EVENT-CMS ',ECM,NEV
823           PMASS1 = PHO_PMASS(IFPAB(1),0)-SQRT(PVIRT(1))
824           PMASS2 = PHO_PMASS(IFPAB(2),0)-SQRT(PVIRT(2))
825           CALL PHO_PECMS(1,PMASS1,PMASS2,ECM,PCM,EE)
826           E1 = EE
827           E2 = ECM-EE
828           THETA = 0.D0
829           PHI   = 0.D0
830           CALL PHO_FIXCOL(E1,E2,THETA,PHI,NEV)
831           KEVENT = 0
832
833 C  fixed-energy events for collider setup with crossing angle
834         ELSE IF(CNAME.EQ.'EVENT-COLL') THEN
835           READ(NUMBER,*) E1,E2,THETA,PHI,NEV
836           WRITE(LO,*) 'EVENT-COLL',E1,E2,THETA,PHI,NEV
837           CALL PHO_FIXCOL(E1,E2,THETA,PHI,NEV)
838           KEVENT = 0
839
840 C  unknown data card
841         ELSE
842           WRITE(LO,*) 'PHO_INIT: unknown data card: ',CNAME,NUMBER
843         ENDIF
844
845       GOTO 1200
846  1300 CONTINUE
847       WRITE(LO,*) ' RETURN'
848
849       END
850
851 *$ CREATE PHO_SETMDL.FOR
852 *COPY PHO_SETMDL
853 CDECK  ID>, PHO_SETMDL
854       SUBROUTINE PHO_SETMDL(INDX,IVAL,IMODE)
855 C**********************************************************************
856 C
857 C     set model switches
858 C
859 C     input:  INDX       model parameter number
860 C                        (positive: ISWMDL, negative: IPAMDL)
861 C             IVAL       new value
862 C             IMODE      -1  print value of parameter INDX
863 C                        1   set new value
864 C                        -2  print current settings
865 C
866 C**********************************************************************
867       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
868       SAVE
869
870 C  input/output channels
871       INTEGER LI,LO
872       COMMON /POINOU/ LI,LO
873 C  model switches and parameters
874       CHARACTER*8 MDLNA
875       INTEGER ISWMDL,IPAMDL
876       DOUBLE PRECISION PARMDL
877       COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
878
879       IF(IMODE.EQ.-2) THEN
880         WRITE(LO,'(/1X,A,/1X,A,/)') 'PHO_SETMDL: current settings',
881      &                             '----------------------------'
882         DO 100 I=1,48,3
883           IF(ISWMDL(I).EQ.-9999) GOTO 200
884           IF(ISWMDL(I+1).EQ.-9999) THEN
885             WRITE(LO,'(5X,I3,A1,A,I6)') I,':',MDLNA(I),ISWMDL(I)
886             GOTO 200
887           ELSE IF(ISWMDL(I+2).EQ.-9999) THEN
888             WRITE(LO,'(2(5X,I3,A1,A,I6))') I,':',MDLNA(I),ISWMDL(I),
889      &        I+1,':',MDLNA(I+1),ISWMDL(I+1)
890             GOTO 200
891           ELSE
892             WRITE(LO,'(3(5X,I3,A1,A,I6))')
893      &        (I+K,':',MDLNA(I+K),ISWMDL(I+K),K=0,2)
894           ENDIF
895  100    CONTINUE
896  200    CONTINUE
897       ELSE IF(IMODE.EQ.-1) THEN
898         WRITE(LO,'(1X,A,1X,A,I6)')
899      &    'PHO_SETMDL:',MDLNA(INDX),ISWMDL(INDX)
900       ELSE IF(IMODE.EQ.1) THEN
901         IF(INDX.GT.0) THEN
902           IF(ISWMDL(INDX).NE.IVAL) THEN
903             WRITE(LO,'(1X,A,I4,1X,A,2I6)')
904      &        'PHO_SETMDL:ISWMDL(OLD/NEW):',
905      &        INDX,MDLNA(INDX),ISWMDL(INDX),IVAL
906             ISWMDL(INDX) = IVAL
907           ENDIF
908         ELSE IF(INDX.LT.0) THEN
909           IF(IPAMDL(-INDX).NE.IVAL) THEN
910             WRITE(LO,'(1X,A,I4,1X,2I6)') 'PHO_SETMDL:IPAMDL(OLD/NEW):',
911      &        -INDX,IPAMDL(-INDX),IVAL
912             IPAMDL(-INDX) = IVAL
913           ENDIF
914         ENDIF
915       ELSE
916         WRITE(LO,'(/1X,A,I6)')
917      &    'PHO_SETMDL:ERROR: unsupported mode',IMODE
918       ENDIF
919       END
920
921 *$ CREATE PHO_DATINI.FOR
922 *COPY PHO_DATINI
923 CDECK  ID>, PHO_DATINI
924       SUBROUTINE PHO_DATINI
925 C*********************************************************************
926 C
927 C     initialization of variables and switches
928 C
929 C*********************************************************************
930       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
931       SAVE
932
933 C  input/output channels
934       INTEGER LI,LO
935       COMMON /POINOU/ LI,LO
936 C  some constants
937       DOUBLE PRECISION PI,PI2,PI4,GEV2MB,Q_ch,Q_ch2,Q_ch4
938       COMMON /POCONS/ PI,PI2,PI4,GEV2MB,
939      &  Q_ch(-6:6),Q_ch2(-6:6),Q_ch4(-6:6)
940 C  event debugging information
941       INTEGER NMAXD
942       PARAMETER (NMAXD=100)
943       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
944      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
945       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
946      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
947 C  event weights and generated cross section
948       INTEGER IPOWGC,ISWCUT,IVWGHT
949       DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
950       COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
951      &                IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
952 C  scale parameters for parton model calculations
953       INTEGER NQQAL,NQQALI,NQQALF,NQQPD
954       DOUBLE PRECISION AQQAL,AQQALI,AQQALF,AQQPD
955       COMMON /POHSCL/ AQQAL,AQQALI,AQQALF,AQQPD,
956      &                NQQAL,NQQALI,NQQALF,NQQPD
957 C  integration precision for hard cross sections (obsolete)
958       INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
959       COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
960 C  hard scattering parameters used for most recent hard interaction
961       INTEGER NFbeta,NF
962       DOUBLE PRECISION ALQCD2,BQCD
963       COMMON /POHAPA/ ALQCD2(3,4),BQCD(4),NFbeta,NF
964 C  cut probability distribution
965       INTEGER IEETA1,IIMAX,KKMAX
966       PARAMETER( IEETA1=20, IIMAX=20, KKMAX=20 )
967       INTEGER IEEMAX,IMAX,KMAX
968       REAL PROB
969       DOUBLE PRECISION EPTAB
970       COMMON /POPROB/ PROB(4,IEETA1,0:IIMAX,0:KKMAX),EPTAB(4,IEETA1),
971      &                IEEMAX,IMAX,KMAX
972 C  gamma-lepton or gamma-hadron vertex information
973       INTEGER IGHEL,IDPSRC,IDBSRC
974       DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
975      &                 RADSRC,AMSRC,GAMSRC
976       COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
977      &                GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
978      &                IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
979 C  photon flux kinematics and cuts
980       DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
981      &                 YMIN1,YMAX1,YMIN2,YMAX2,
982      &                 Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
983      &                 THMIN1,THMAX1,THMIN2,THMAX2
984       INTEGER          ITAG1,ITAG2
985       COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
986      &                YMIN1,YMAX1,YMIN2,YMAX2,
987      &                Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
988      &                THMIN1,THMAX1,THMIN2,THMAX2,
989      &                ITAG1,ITAG2
990 C  obsolete cut-off information
991       DOUBLE PRECISION PTCUT,PTANO,FPS,FPH,PSOMIN,XSOMIN
992       COMMON /POCUT1/ PTCUT(4),PTANO(4),FPS(4),FPH(4),PSOMIN,XSOMIN
993 C  global event kinematics and particle IDs
994       INTEGER IFPAP,IFPAB
995       DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
996       COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
997 C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
998       INTEGER IDEQP,IDEQB,IHFLD,IHFLS
999       DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
1000       COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
1001      &                IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
1002 C  some hadron information, will be deleted in future versions
1003       INTEGER NFS
1004       DOUBLE PRECISION QMASS,BET,PCOUDI,PNORM,VALPRG
1005       COMMON /POHDRN/ QMASS(6),BET,PCOUDI,PNORM,VALPRG(2),NFS
1006 C  model switches and parameters
1007       CHARACTER*8 MDLNA
1008       INTEGER ISWMDL,IPAMDL
1009       DOUBLE PRECISION PARMDL
1010       COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
1011 C  general process information
1012       INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
1013       COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
1014 C  parameters of the "simple" Vector Dominance Model
1015       DOUBLE PRECISION VMAS,GAMM,RMIN,RMAX,VMSL,VMFA
1016       COMMON /POSVDM/ VMAS(4),GAMM(4),RMIN(4),RMAX(4),VMSL(4),VMFA(4)
1017 C  parameters for DGLAP backward evolution in ISR
1018       INTEGER NFSISR
1019       DOUBLE PRECISION Q2MISR,PMISR,ZMISR,AL2ISR
1020       COMMON /PODGL1/ Q2MISR(2),PMISR(2),ZMISR(2),AL2ISR(2),NFSISR
1021 C  particles created by initial state evolution
1022       INTEGER MXISR1,MXISR2
1023       PARAMETER ( MXISR1 = 150, MXISR2 = 50 )
1024       INTEGER IFLISR,IPOISR,IMXISR
1025       DOUBLE PRECISION PHISR
1026       COMMON /POPISR/ IFLISR(2,MXISR1),PHISR(2,4,MXISR1),
1027      &                IPOISR(2,2,MXISR2),IMXISR(2)
1028 C  names of hard scattering processes
1029       INTEGER Max_pro_1
1030       PARAMETER ( Max_pro_1 = 16 )
1031       CHARACTER*18 PROC
1032       COMMON /POHPRO/ PROC(0:Max_pro_1)
1033 C  hard cross sections and MC selection weights
1034       INTEGER Max_pro_2
1035       PARAMETER ( Max_pro_2 = 16 )
1036       INTEGER IHa_last,IHb_last,MH_pro_on,MH_tried,
1037      &  MH_acc_1,MH_acc_2
1038       DOUBLE PRECISION Hfac,HWgx,HSig,Hdpt,HEcm_last,HQ2a_last,HQ2b_last
1039       COMMON /POHRCS/ Hfac(-1:Max_pro_2),HWgx(-1:Max_pro_2),
1040      &  HSig(-1:Max_pro_2),Hdpt(-1:Max_pro_2),
1041      &  HEcm_last,HQ2a_last,HQ2b_last,IHa_last,IHb_last,
1042      &  MH_pro_on(-1:Max_pro_2,0:4),MH_tried(-1:Max_pro_2,0:4),
1043      &  MH_acc_1(-1:Max_pro_2,0:4),MH_acc_2(-1:Max_pro_2,0:4)
1044 C  interpolation tables for hard cross section and MC selection weights
1045       INTEGER Max_tab_E,Max_tab_Q2,Max_pro_tab
1046       PARAMETER ( Max_tab_E = 20, Max_tab_Q2 = 10, Max_pro_tab = 16 )
1047       INTEGER IH_Q2a_up,IH_Q2b_up,IH_Ecm_up
1048       DOUBLE PRECISION Hfac_tab,HWgx_tab,HSig_tab,Hdpt_tab,
1049      &  HQ2a_tab,HQ2b_tab,HEcm_tab
1050       COMMON /POHTAB/
1051      &  Hfac_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1052      &  HWgx_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1053      &  HSig_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1054      &  Hdpt_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1055      &  HQ2a_tab(1:Max_tab_Q2,0:4),HQ2b_tab(1:Max_tab_Q2,0:4),
1056      &  HEcm_tab(1:Max_tab_E,0:4),
1057      &  IH_Q2a_up(0:4),IH_Q2b_up(0:4),IH_Ecm_up(0:4)
1058
1059 C  initialize /POCONS/
1060       PI   = ATAN(1.D0)*4.D0
1061       PI2  = 2.D0*PI
1062       PI4  = 2.D0*PI2
1063 C  GeV**-2 --> millibarn (multiply by GEV2MB to get mb as units)
1064       GEV2MB = 0.389365D0
1065 C  precalculate quark charges
1066       do i=1,6
1067         Q_ch(i) = dble(2-3*mod(i,2))/3.D0
1068         Q_ch(-i) = -Q_ch(i)
1069
1070         Q_ch2(i) = Q_ch(i)**2
1071         Q_ch2(-i) = Q_ch2(i)
1072
1073         Q_ch4(i) = Q_ch2(i)**2
1074         Q_ch4(-i) = Q_ch4(i)
1075       enddo
1076       Q_ch(0)  = 0.D0
1077       Q_ch2(0) = 0.D0
1078       Q_ch4(0) = 0.D0
1079
1080 C  initialize /GLOCMS/
1081       ECM    = 50.D0
1082       PMASS(1) = 0.D0
1083       PVIRT(1) = 0.D0
1084       PMASS(2) = 0.D0
1085       PVIRT(2) = 0.D0
1086       IFPAP(1) = 22
1087       IFPAP(2) = 22
1088 C  initialize /HADVAL/
1089       IHFLD(1,1) = 0
1090       IHFLD(1,2) = 0
1091       IHFLD(2,1) = 0
1092       IHFLD(2,2) = 0
1093       IHFLS(1) = 1
1094       IHFLS(2) = 1
1095 C  initialize /MODELS/
1096       ISWMDL(1)  = 3
1097       MDLNA(1)  = 'AMPL MOD'
1098       ISWMDL(2)  = 1
1099       MDLNA(2)  = 'MIN-BIAS'
1100       ISWMDL(3)  = 1
1101       MDLNA(3)  = 'PTS DISH'
1102       ISWMDL(4)  = 1
1103       MDLNA(4)  = 'PTS DISP'
1104       ISWMDL(5)  = 2
1105       MDLNA(5)  = 'PTS ASSI'
1106       ISWMDL(6)  = 3
1107       MDLNA(6)  = 'HADRONIZ'
1108       ISWMDL(7)  = 2
1109       MDLNA(7)  = 'MASS COR'
1110       ISWMDL(8)  = 3
1111       MDLNA(8)  = 'PAR SHOW'
1112       ISWMDL(9)  = 0
1113       MDLNA(9)  = 'GLU SPLI'
1114       ISWMDL(10) = 2
1115       MDLNA(10) = 'VIRT PHO'
1116       ISWMDL(11) = 0
1117       MDLNA(11) = 'LARGE NC'
1118       ISWMDL(12) = 0
1119       MDLNA(12) = 'LIPA POM'
1120       ISWMDL(13) = 1
1121       MDLNA(13) = 'QELAS VM'
1122       ISWMDL(14) = 2
1123       MDLNA(14) = 'ENHA GRA'
1124       ISWMDL(15) = 4
1125       MDLNA(15) = 'MULT SCA'
1126       ISWMDL(16) = 4
1127       MDLNA(16) = 'MULT DIF'
1128       ISWMDL(17) = 4
1129       MDLNA(17) = 'MULT CDF'
1130       ISWMDL(18) = 0
1131       MDLNA(18) = 'BALAN PT'
1132       ISWMDL(19) = 1
1133       MDLNA(19) = 'POMV FLA'
1134       ISWMDL(20) = 0
1135       MDLNA(20) = 'SEA  FLA'
1136       ISWMDL(21) = 2
1137       MDLNA(21) = 'SPIN DEC'
1138       ISWMDL(22) = 1
1139       MDLNA(22) = 'DIF.MASS'
1140       ISWMDL(23) = 1
1141       MDLNA(23) = 'DIFF RES'
1142       ISWMDL(24) = 0
1143       MDLNA(24) = 'PTS HPOM'
1144       ISWMDL(25) = 0
1145       MDLNA(25) = 'POM CORR'
1146       ISWMDL(26) = 1
1147       MDLNA(26) = 'OVERLAP '
1148       ISWMDL(27) = 0
1149       MDLNA(27) = 'MUL R/AN'
1150       ISWMDL(28) = 1
1151       MDLNA(28) = 'SUR PROB'
1152       ISWMDL(29) = 1
1153       MDLNA(29) = 'PRIMO KT'
1154       ISWMDL(30) = 0
1155       MDLNA(30) = 'DIFF. CS'
1156       ISWMDL(31) = -9999
1157 C  mass-independent sea flavour ratios (for low-mass strings)
1158       PARMDL(1)  = 0.425D0
1159       PARMDL(2)  = 0.425D0
1160       PARMDL(3)  = 0.15D0
1161       PARMDL(4)  = 0.D0
1162       PARMDL(5)  = 0.D0
1163       PARMDL(6)  = 0.D0
1164 C  suppression by energy momentum conservation
1165       PARMDL(8)  = 9.D0
1166       PARMDL(9)  = 7.D0
1167 C  VDM factors
1168       PARMDL(10) = 0.866D0
1169       PARMDL(11) = 0.288D0
1170       PARMDL(12) = 0.288D0
1171       PARMDL(13) = 0.288D0
1172       PARMDL(14) = 0.866D0
1173       PARMDL(15) = 0.288D0
1174       PARMDL(16) = 0.288D0
1175       PARMDL(17) = 0.288D0
1176       PARMDL(18) = 0.D0
1177 C  lower energy limit for initialization
1178       PARMDL(19) = 5.D0
1179 C  soft pt for hard scattering remnants
1180       PARMDL(20) = 5.D0
1181 C  low energy beta of soft pt distribution 1
1182       PARMDL(21) = 4.5D0
1183 C  high energy beta of soft pt distribution 1
1184       PARMDL(22) = 3.0D0
1185 C  low energy beta of soft pt distribution 0
1186       PARMDL(23) = 2.5D0
1187 C  high energy beta of soft pt distribution 0
1188       PARMDL(24) = 0.4D0
1189 C  effective quark mass in photon wave function
1190       PARMDL(25) = 0.2D0
1191 C  normalization of unevolved Pomeron PDFs
1192       PARMDL(26) = 0.3D0
1193 C  effective VDM parameters for Q**2 dependence of cross section
1194       PARMDL(27) = 0.65D0
1195       PARMDL(28) = 0.08D0
1196       PARMDL(29) = 0.05D0
1197       PARMDL(30) = 0.22D0
1198       PARMDL(31) = 0.589824D0
1199       PARMDL(32) = 0.609961D0
1200       PARMDL(33) = 1.038361D0
1201       PARMDL(34) = 1.96D0
1202 C  Q**2 suppression of multiple interactions
1203       PARMDL(35) = 0.59D0
1204 C  pt cutoff defaults
1205       PARMDL(36) = 2.5D0
1206       PARMDL(37) = 2.5D0
1207       PARMDL(38) = 2.5D0
1208       PARMDL(39) = 2.5D0
1209 C  enhancement factor for diffractive cross sections
1210       PARMDL(40) = 1.D0
1211       PARMDL(41) = 1.D0
1212       PARMDL(42) = 1.D0
1213 C  mass in soft pt distribution
1214       PARMDL(43) = 0.D0
1215 C  maximum of x allowed for leading particle
1216       PARMDL(44) = 0.9D0
1217 C  max. mass sampled in diffraction
1218       PARMDL(45) = sqrt(0.4D0)
1219 C  mass threshold in diffraction (2pi mass)
1220       PARMDL(46) = 0.3D0
1221 C  regularization of slope parameter in diffraction
1222       PARMDL(47) = 4.D0
1223 C  renormalized intercept for enhanced graphs
1224       PARMDL(48) = 1.08D0
1225 C  coherence constraint for diff. cross sections
1226       PARMDL(49) = sqrt(0.05D0)
1227 C  exponents of x distributions
1228 C  baryon
1229       PARMDL(50) = 1.5D0
1230       PARMDL(51) = -0.5D0
1231       PARMDL(52) = -0.99D0
1232       PARMDL(53) = -0.99D0
1233 C  meson (non-strangeness part)
1234       PARMDL(54) = -0.5D0
1235       PARMDL(55) = -0.5D0
1236       PARMDL(56) = -0.99D0
1237       PARMDL(57) = -0.99D0
1238 C  meson (strangeness part)
1239       PARMDL(58) = -0.2D0
1240       PARMDL(59) = -0.2D0
1241       PARMDL(60) = -0.99D0
1242       PARMDL(61) = -0.99D0
1243 C  particle remnant (no valence quarks)
1244       PARMDL(62) = -0.5D0
1245       PARMDL(63) = -0.5D0
1246       PARMDL(64) = -0.99D0
1247       PARMDL(65) = -0.99D0
1248 C  ratio beetween triple-pomeron/reggeon couplings grrp/gppp
1249       PARMDL(66) = 10.D0
1250 C  ratio beetween triple-pomeron/reggeon couplings gppr/gppp
1251       PARMDL(67) = 10.D0
1252 C  min. abs(t) in diffraction
1253       PARMDL(68) = 0.D0
1254 C  max. abs(t) in diffraction
1255       PARMDL(69) = 10.D0
1256 C  min. mass for elastic pomerons in central diffraction
1257       PARMDL(70) = 2.D0
1258 C  min. mass of diffractive blob in central diffraction
1259       PARMDL(71) = 2.D0
1260 C  min. Feynman x cut in central diffraction
1261       PARMDL(72) = 0.D0
1262 C  direct pomeron coupling
1263       PARMDL(74) = 0.D0
1264 C  relative deviation allowed for energy-momentum conservation
1265 C  energy-momentum relative deviation
1266       PARMDL(75) = 0.01D0
1267 C  transverse momentum deviation
1268       PARMDL(76) = 0.01D0
1269 C  couplings for unitarization in diffraction
1270 C  non-unitarized pomeron coupling (sqrt(mb))
1271       PARMDL(77)  = 3.D0
1272 C  rescaling factor for pomeron PDF
1273       PARMDL(78)  = 3.D0
1274 C  coupling probabilities
1275       PARMDL(79)  = 1.D0
1276       PARMDL(80)  = 0.D0
1277 C  scales to calculate alpha-s of matrix element
1278       PARMDL(81) = 1.D0
1279       PARMDL(82) = 1.D0
1280       PARMDL(83) = 1.D0
1281 C  scales to calculate alpha-s of initial state radiation
1282       PARMDL(84) = 1.D0
1283       PARMDL(85) = 1.D0
1284       PARMDL(86) = 1.D0
1285 C  scales to calculate alpha-s of final state radiation
1286       PARMDL(87) = 1.D0
1287       PARMDL(88) = 1.D0
1288       PARMDL(89) = 1.D0
1289 C  scales to calculate PDFs
1290       PARMDL(90) = 1.D0
1291       PARMDL(91) = 1.D0
1292       PARMDL(92) = 1.D0
1293 C  scale for ISR starting virtuality
1294       PARMDL(93) = 1.D0
1295 C  min. virtuality to generate time-like showers in ISR
1296       PARMDL(94) = 2.D0
1297 C  factor to scale the max. allowed time-like parton shower virtuality
1298       PARMDL(95) = 4.D0
1299 C  max. transverse momentum for primordial kt
1300       PARMDL(100) = 2.D0
1301 C  weight factors for pt-distribution
1302       PARMDL(101) = 2.D0
1303       PARMDL(102) = 2.D0
1304       PARMDL(103) = 4.D0
1305       PARMDL(104) = 2.D0
1306       PARMDL(105) = 6.D0
1307       PARMDL(106) = 4.D0
1308 C
1309 *     PARMDL(110-125)  reserved for hard scattering
1310 C  currently chosen scales for hard scattering
1311       DO 10 I=1,16
1312         PARMDL(109+I) = 0.D0
1313  10   CONTINUE
1314 C  virtuality cutoff in initial state evolution
1315       PARMDL(126) = PARMDL(36)**2
1316       PARMDL(127) = PARMDL(37)**2
1317       PARMDL(128) = PARMDL(38)**2
1318       PARMDL(129) = PARMDL(39)**2
1319 C  virtuality cutoff for direct contribution to photon PDF
1320       PARMDL(130) = 1.D30
1321       PARMDL(131) = 1.D30
1322       PARMDL(132) = 1.D30
1323       PARMDL(133) = 1.D30
1324 C  fraction of events without popcorn
1325       PARMDL(134) = -1.D0
1326 C  fraction of diquarks with spin 1 (relative to sum of spin 1 and 0)
1327       PARMDL(135) = 0.5D0
1328 C  soft color re-connection (fraction)
1329 C  g g final state
1330       PARMDL(140) = 1.D0/64.D0
1331 C  g q final state
1332       PARMDL(141) = 1.D0/24.D0
1333 C  q q final state
1334       PARMDL(142) = 1.D0/9.D0
1335 C  effective scale in Drees-Godbole like suppresion in photon PDF
1336       PARMDL(144) = 0.766D0**2
1337 C  QCD scales (if PDF scales are not used, 4 active flavours)
1338       PARMDL(145) = 0.2D0**2
1339       PARMDL(146) = 0.2D0**2
1340       PARMDL(147) = 0.2D0**2
1341 C  threshold scales for variable flavour calculation (GeV**2)
1342       PARMDL(148) = 1.5D0**2
1343       PARMDL(149) = 4.5D0**2
1344       PARMDL(150) = 175.D0**2
1345 C  constituent quark masses
1346       PARMDL(151) = 0.3D0
1347       PARMDL(152) = 0.3D0
1348       PARMDL(153) = 0.5D0
1349       PARMDL(154) = 1.6D0
1350       PARMDL(155) = 5.D0
1351       PARMDL(156) = 174.D0
1352 C  min. masses of valence quark
1353       PARMDL(157) = 0.3D0
1354 C  min. masses of valence diquark
1355       PARMDL(158) = 0.8D0
1356 C  min. mass of sea quark
1357       PARMDL(159) = 0.D0
1358 C  suppression of strange quarks as photon valences
1359       PARMDL(160) = 0.2D0
1360 C  min. masses for strings (used in PHO_SOFTXX)
1361       PARMDL(161) = 1.D0
1362       PARMDL(162) = 1.D0
1363       PARMDL(163) = 1.D0
1364       PARMDL(164) = 1.D0
1365 C  min. momentum fraction for soft processes
1366       PARMDL(165) = 0.3D0
1367 C  min. phase space for x-sampling
1368       PARMDL(166) = 0.135D0
1369 C  Ross-Stodolsky exponent
1370       PARMDL(170) = 4.2D0
1371 C  cutoff on photon-pomeron invariant mass in hadron-hadron collisions
1372       PARMDL(175) = 2.D0
1373
1374 **sr
1375 *  extra factor multiplying difference between Goulianos and PHOJET-
1376 *  diff. cross sections
1377       PARMDL(200) = 0.6D0
1378 **
1379
1380 C  complex amplitudes, eikonal functions
1381       IPAMDL(1)  = 0
1382 C  allow for Reggeon cuts
1383       IPAMDL(2)  = 1
1384 C  decay of hadron resonances in diffraction (0 iso, 1 trans, 2 long)
1385       IPAMDL(3)  = 0
1386 C  polarization of photon resonances (0 none, 1 trans, 2 long)
1387       IPAMDL(4)  = 1
1388 C  pt of valence partons
1389       IPAMDL(5)  = 1
1390 C  pt of hard scattering remnant
1391       IPAMDL(6)  = 2
1392 C  running cutoff for hard scattering
1393       IPAMDL(7)  = 1
1394 C  intercept used for the calculation of enhanced graphs
1395       IPAMDL(8)  = 1
1396 C  effective slope of hard scattering amplitde
1397       IPAMDL(9)  = 1
1398 C  mass dependence of slope parameters
1399       IPAMDL(10) = 0
1400 C  lepton-photon vertex 1
1401       IPAMDL(11) = 0
1402 C  lepton-photon vertex 2
1403       IPAMDL(12) = 0
1404 C  call by DPMJET
1405       IPAMDL(13) = 0
1406 C  method to sample x distributions
1407       IPAMDL(14) = 3
1408 C  energy-momentum check
1409       IPAMDL(15) = 1
1410 C  phase space correction for DPMJET interface
1411       IPAMDL(16) = 1
1412 C  fragment strings from projectile/target/central diff. separately
1413       IPAMDL(17) = 1
1414 C  method to construct strings for hard interactions
1415       IPAMDL(18) = 1
1416 C  method to construct strings for soft sea (pomeron cuts)
1417       IPAMDL(19) = 0
1418 C  method to construct strings in pomeron interactions
1419       IPAMDL(20) = 0
1420 C  soft color re-connection
1421       IPAMDL(21) = 0
1422 C  resummation of triple- and loop-Pomeron
1423       IPAMDL(24) = 1
1424 C  resummation of X iterated triple-Pomeron
1425       IPAMDL(25) = 1
1426 C  dimension of interpolation table for weights in hard scattering
1427       IPAMDL(30) = Max_tab_E
1428 C  dimension of interpolation table for pomeron cut distribution
1429       IPAMDL(31) = IEETA1
1430 C  number of cut soft pomerons (restriction by field dimension)
1431       IPAMDL(32) = IIMAX
1432 C  number of cut hard pomerons (restriction by field dimension)
1433       IPAMDL(33) = KKMAX
1434 C  tau pair production in direct photon-photon collisions
1435       IPAMDL(64) = 0
1436 C  currently chosen scales for hard scattering
1437 C  ATTENTION:   IPAMDL(65-80)  reserved for hard scattering!
1438       DO 15 I=1,16
1439         IPAMDL(64+I) = -99999
1440  15   CONTINUE
1441 C  scales to calculate alpha-s of matrix element
1442       IPAMDL(81) = 1
1443       IPAMDL(82) = 1
1444       IPAMDL(83) = 1
1445 C  scales to calculate alpha-s of initial state radiation
1446       IPAMDL(84) = 1
1447       IPAMDL(85) = 1
1448       IPAMDL(86) = 1
1449 C  scales to calculate alpha-s of final state radiation
1450       IPAMDL(87) = 1
1451       IPAMDL(88) = 1
1452       IPAMDL(89) = 1
1453 C  scales to calculate PDFs
1454       IPAMDL(90) = 1
1455       IPAMDL(91) = 1
1456       IPAMDL(92) = 1
1457 C  where to get the parameter sets from
1458       IPAMDL(99) = 1
1459 C  program PHO_ABORT for fatal errors (simulation of division by zero)
1460       IPAMDL(100) = 0
1461 C  initial state parton showers for all / hardest interaction(s)
1462       IPAMDL(101) = 1
1463 C  final state parton showers for all / hardest interaction(s)
1464       IPAMDL(102) = 1
1465 C  initial virtuality for ISR generation
1466       IPAMDL(109) = 1
1467 C  qqbar-gamma coupling in initial state showers
1468       IPAMDL(110) = 1
1469 C  generation of time-like showers during ISR
1470       IPAMDL(111) = 1
1471 C  reweighting of multiple soft contributions for virtual photons
1472       IPAMDL(114) = 1
1473 C  reweighting / use photon virtuality in photon PDF calculations
1474       IPAMDL(115) = 0
1475 C  use full QPM model incl. interference terms (direct part in gam-gam)
1476       IPAMDL(116) = 0
1477 C  matching sigma_tot to F2 as given by parton density at high Q2
1478       IPAMDL(117) = 1
1479 C  use virtuality of target in F2 calculations (two-gamma only)
1480       IPAMDL(118) = 1
1481 C  calculation of alpha_em
1482       IPAMDL(120) = 1
1483 C  strict pt cutoff for gamma-gamma events
1484       IPAMDL(121) = 0
1485 C  photon virtuality sampled in photon flux approximations
1486       IPAMDL(174) = 1
1487 C  photon-pomeron: 0,1,2: both,left,right photon emission
1488       IPAMDL(175) = 0
1489 C  keep full history information in PHOJET-JETSET interface
1490       IPAMDL(178) = 1
1491 C  max. number of conservation law violations allowed in one run
1492       IPAMDL(179) = 20
1493 C  selection of soft X values
1494 C  max. iteration number in PHO_SELSXS
1495       IPAMDL(180) = 50
1496 C  max. iteration number in PHO_SELSXR
1497       IPAMDL(181) = 200
1498 C  max. iteration number in PHO_SELSX2
1499       IPAMDL(182) = 100
1500 C  max. iteration number in PHO_SELSXI
1501       IPAMDL(183) = 50
1502
1503 C  initialize /PROBAB/
1504       IEEMAX = IEETA1
1505       IMAX   = IIMAX
1506       KMAX   = KKMAX
1507
1508       DO 20 I=1,30
1509         PARMDL(300+I) = -100000.D0
1510  20   CONTINUE
1511 C  initialize /POHDRN/
1512       QMASS(1) =  PARMDL(151)
1513       QMASS(2) =  PARMDL(152)
1514       QMASS(3) =  PARMDL(153)
1515       QMASS(4) =  PARMDL(154)
1516       QMASS(5) =  PARMDL(155)
1517       QMASS(6) =  PARMDL(156)
1518       BET      = 8.D0
1519       PCOUDI   = 0.D0
1520       VALPRG(1) = 1.D0
1521       VALPRG(2) = 1.D0
1522 C  number of light flavours (quarks treated as massless)
1523       NFS      = 4
1524 C  initialize /POCUT1/
1525       PTCUT(1) = PARMDL(36)
1526       PTCUT(2) = PARMDL(37)
1527       PTCUT(3) = PARMDL(38)
1528       PTCUT(4) = PARMDL(39)
1529       PSOMIN = 0.D0
1530       XSOMIN = 0.D0
1531 C  initialize /POHAPA/
1532       NFbeta  = 4
1533       NF      = 4
1534       BQCD(1) = PI4/(11.D0-(2.D0/3.D0)*3)
1535       BQCD(2) = PI4/(11.D0-(2.D0/3.D0)*4)
1536       BQCD(3) = PI4/(11.D0-(2.D0/3.D0)*5)
1537       BQCD(4) = PI4/(11.D0-(2.D0/3.D0)*6)
1538 C  initialize /POGAUP/
1539       NGAUP1 = 12
1540       NGAUP2 = 12
1541       NGAUET = 16
1542       NGAUIN = 12
1543       NGAUSO = 96
1544 C  initialize //
1545       DO 30 I=1,100
1546         IDEB(I) = 0
1547  30   CONTINUE
1548 C  initialize /PROCES/
1549       DO 35 I=1,11
1550         IPRON(I,1) = 1
1551  35   CONTINUE
1552
1553 C  DPMJET default: no elastic scattering
1554       IPRON(2,1) = 0
1555
1556       DO 36 K=2,4
1557         DO 37 I=2,11
1558           IPRON(I,K) = 0
1559  37     CONTINUE
1560         IPRON(1,K) = 1
1561         IPRON(8,K) = 1
1562  36   CONTINUE
1563 C  initialize /POSVDM/
1564       TWOPIM = 0.28D0
1565       RMIN(1) = 0.285D0
1566       RMIN(2) = 0.45D0
1567       RMIN(3) = 1.D0
1568       RMIN(4) = TWOPIM
1569       VMAS(1) = 0.770D0
1570       VMAS(2) = 0.787D0
1571       VMAS(3) = 1.02D0
1572       VMAS(4) = TWOPIM
1573       GAMM(1) = 0.155D0
1574       GAMM(2) = 0.01D0
1575       GAMM(3) = 0.0045D0
1576       GAMM(4) = 1.D0
1577       RMAX(1) = VMAS(1)+TWOPIM
1578       RMAX(2) = VMAS(2)+TWOPIM
1579       RMAX(3) = VMAS(3)+TWOPIM
1580       RMAX(4) = VMAS(1)+TWOPIM
1581       VMSL(1) = 11.D0
1582       VMSL(2) = 10.D0
1583       VMSL(3) = 6.D0
1584       VMSL(4) = 4.D0
1585       VMFA(1) = 0.0033D0
1586       VMFA(2) = 0.00036D0
1587       VMFA(3) = 0.0002D0
1588       VMFA(4) = 0.0002D0
1589 C  initialize /PODGL1/
1590       Q2MISR(1) = PARMDL(36)**2
1591       Q2MISR(2) = PARMDL(36)**2
1592       PMISR(1) = 1.D0
1593       PMISR(2) = 1.D0
1594       ZMISR(1) = 0.001D0
1595       ZMISR(2) = 0.001D0
1596       AL2ISR(1) = 0.046D0
1597       AL2ISR(2) = 0.046D0
1598       NFSISR  = 4
1599 C  initialize /POPISR/
1600       DO 40 I=1,50
1601         IPOISR(1,2,I) = 0
1602         IPOISR(2,2,I) = 0
1603  40   CONTINUE
1604 C  initialize /POHPRO/
1605       PROC(0) = 'sum over processes'
1606       PROC(1) = 'G  +G  --> G  +G  '
1607       PROC(2) = 'Q  +QB --> G  +G  '
1608       PROC(3) = 'G  +Q  --> G  +Q  '
1609       PROC(4) = 'G  +G  --> Q  +QB '
1610       PROC(5) = 'Q  +QB --> Q  +QB '
1611       PROC(6) = 'Q  +QB --> QP +QBP'
1612       PROC(7) = 'Q  +Q  --> Q  +Q  '
1613       PROC(8) = 'Q  +QP --> Q  +QP '
1614       PROC(9) = 'resolved processes'
1615       PROC(10) = 'gam+Q  --> G  +Q  '
1616       PROC(11) = 'gam+G  --> Q  +QB '
1617       PROC(12) = 'Q  +gam--> G  +Q  '
1618       PROC(13) = 'G  +gam--> Q  +QB '
1619       PROC(14) = 'gam+gam--> Q  +QB '
1620       PROC(15) = 'direct processes  '
1621       PROC(16) = 'gam+gam--> l+ +l- '
1622
1623 C  initialize /POHRCS/
1624       do M=1,Max_pro_2
1625         HWgx(M) = 0.D0
1626         HSig(M) = 0.D0
1627         Hdpt(M) = 0.D0
1628       enddo
1629       DO I=0,4
1630         DO M=-1,Max_pro_2
1631 C  switch all hard subprocesses on
1632           MH_pro_on(M,I) = 1
1633 C  reset all counters
1634           MH_tried(M,I) = 0
1635           MH_acc_1(M,I) = 0
1636           MH_acc_2(M,I) = 0
1637         ENDDO
1638         MH_pro_on(16,I) = 0
1639       ENDDO
1640
1641 C  initialize /POHTAB/
1642       do I=0,4
1643         IH_Ecm_up(I) = 0
1644         IH_Q2a_up(I) = 0
1645         IH_Q2b_up(I) = 0
1646         HEcm_tab(1,I) = 0.D0
1647       enddo
1648       HEcm_last = 0.D0
1649       IHa_last = 0.D0
1650       IHb_last = 0.D0
1651
1652 C  initialize /POFSRC/
1653       IGHEL(1) = -1
1654       IGHEL(2) = -1
1655 C  initialize /LEPCUT/
1656       ECMIN = 5.D0
1657       ECMAX = 1.D+30
1658       EEMIN1 = 1.D0
1659       EEMIN2 = 1.D0
1660       YMAX1 = -1.D0
1661       YMAX2 = -1.D0
1662       THMIN1 = 0.D0
1663       THMAX1 = PI
1664       THMIN2 = 0.D0
1665       THMAX2 = PI
1666       ITAG1 = 1
1667       ITAG2 = 1
1668 C  initialize /POWGHT/
1669       DO 70 I=1,20
1670         HSWCUT(I) = 0.D0
1671         ISWCUT(I) = 0
1672  70   CONTINUE
1673       EVWGHT(1) = 1.D0
1674       IVWGHT(1) = 0
1675       SIGGEN(1) = 0.D0
1676       SIGGEN(2) = 0.D0
1677       SIGGEN(3) = 0.D0
1678       SIGGEN(4) = 0.D0
1679
1680       END
1681
1682 *$ CREATE PHO_PARDAT.FOR
1683 *COPY PHO_PARDAT
1684 CDECK  ID>, PHO_PARDAT
1685       SUBROUTINE PHO_PARDAT
1686 C***********************************************************************
1687 C
1688 C     particle data (based on 1996 PDG naming scheme and data tables)
1689 C
1690 C***********************************************************************
1691
1692       IMPLICIT NONE
1693
1694       SAVE
1695
1696 C  input/output channels
1697       INTEGER LI,LO
1698       COMMON /POINOU/ LI,LO
1699 C  event debugging information
1700       INTEGER NMAXD
1701       PARAMETER (NMAXD=100)
1702       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
1703      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
1704       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
1705      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
1706 C  particle ID translation table
1707       integer         ID_pdg_list,ID_list,ID_pdg_max
1708       character*12    name_list
1709       COMMON /POPAR1/ ID_pdg_list(300),ID_list(577),name_list(300),
1710      &                ID_pdg_max
1711 C  general particle data
1712       double precision xm_list,tau_list,gam_list,
1713      &  xm_psm2_list,xm_vem2_list,xm_b82_list,xm_b102_list,
1714      &  xm_bb82_list,xm_bb102_list
1715       integer          ich3_list,iba3_list,iq_list,
1716      &                 id_psm_list,id_vem_list,id_b8_list,id_b10_list
1717       COMMON /POPAR2/ xm_list(300),tau_list(300),gam_list(300),
1718      &  xm_psm2_list(6,6),xm_vem2_list(6,6),
1719      &  xm_b82_list(6,6,6),xm_b102_list(6,6,6),
1720      &  xm_bb82_list(6,6,6,6),xm_bb102_list(6,6,6,6),
1721      &  ich3_list(300),iba3_list(300),iq_list(3,300),
1722      &  id_psm_list(6,6),id_vem_list(6,6),
1723      &  id_b8_list(6,6,6),id_b10_list(6,6,6)
1724 C  particle decay data
1725       double precision wg_sec_list
1726       integer          idec_list,isec_list
1727       COMMON /POPAR3/ wg_sec_list(500),idec_list(3,300),
1728      &  isec_list(3,500)
1729
1730 C  external functions
1731
1732       integer ipho_pdg2id
1733       double precision pho_pmass
1734
1735 C  local variables for storing data tables
1736
1737       integer number,ich3,iba3,iq_linear,idec_linear,isec_linear,
1738      &  id_psm_linear,id_vem_linear,id_b8_linear,id_b10_linear
1739
1740       dimension number(300),ich3(300),iba3(300),iq_linear(900),
1741      &  idec_linear(900),isec_linear(900),id_psm_linear(36),
1742      &  id_vem_linear(36),id_b8_linear(216),id_b10_linear(216)
1743
1744       double precision xmass,gamma,wg_chan
1745       dimension xmass(300),gamma(300),wg_chan(300)
1746
1747       character*12 name
1748       dimension name(300)
1749
1750       integer i,i1,i2,ii,j,jj,k,l,ichan,i_tab_max,K8,K10,L8,L10
1751       double precision AM1,AM2,AM2P,AM2V,AM82,AM102,AMM
1752
1753       integer itmp
1754
1755       DATA i_tab_max /260/
1756
1757       DATA (number(K),K=    1,  171) /
1758      &     1,     2,     3,     4,     5,     6,  1103,  2101,  2103,
1759      &  2203,  3101,  3103,  3201,  3203,  3303,  4101,  4103,  4201,
1760      &  4203,  4301,  4303,  4403,    81,    82,    90,    91,    92,
1761      &   110,   990,    21,    22,    24,    23,    11,    13,    15,
1762      &    12,    14,    16,   211,   111,   221,   113,   213,   223,
1763      &   331, 10221, 10111, 10211,   333, 10223, 10113, 10213, 20113,
1764      & 20213,   225, 20223, 20221, 20111, 20211,   115,   215, 30223,
1765      & 50223, 40113, 40213, 50221,   335, 60223,   227, 10115, 10215,
1766      & 10333,   117,   217, 30113, 30213, 60221,   337, 20225,   229,
1767      & 30225, 40225,   321,   311,   310,   130,   323,   313, 10313,
1768      & 10323, 20313, 20323, 30313, 30323, 10311, 10321,   325,   315,
1769      & 40313, 40323, 10315, 10325,   317,   327, 20315, 20325,   319,
1770      &   329,   411,   421,   423,   413, 10423,   425,   415,   431,
1771      &   433, 10433,   521,   511,   513,   523,   531,   441,   443,
1772      & 10441, 10443,   445, 20443, 30443, 40443, 50443, 60443,   553,
1773      &   551, 10553,   555, 20553, 10551, 70553, 10555, 30553, 40553,
1774      & 50553, 60553,  2212,  2112, 12112, 12212,  1214,  2124, 22112,
1775      & 22212, 32112, 32212,  2116,  2216, 12116, 12216, 21214, 22124,
1776      & 42112, 42212, 31214, 32124,  1218,  2128,  1114,  2114,  2214/
1777       DATA (number(K),K=  172,  260) /
1778      &  2224, 31114, 32114, 32214, 32224,  1112,  1212,  2122,  2222,
1779      & 11114, 12114, 12214, 12224,  1116,  1216,  2126,  2226, 21112,
1780      & 21212, 22122, 22222, 21114, 22114, 22214, 22224, 11116, 11216,
1781      & 12126, 12226,  1118,  2118,  2218,  2228,  3122, 13122,  3124,
1782      & 23122, 33122, 13124, 43122, 53122,  3126, 13126, 23124,  3128,
1783      & 23126,  3222,  3212,  3112,  3224,  3214,  3114, 13112, 13212,
1784      & 13222, 13114, 13214, 13224, 23112, 23212, 23222,  3116,  3216,
1785      &  3226, 13116, 13216, 13226, 23114, 23214, 23224,  3118,  3218,
1786      &  3228,  3322,  3312,  3324,  3314, 13314, 13324,  3334,  4122,
1787      & 14122,  4222,  4212,  4112,  4232,  4132,  4332,  5122/
1788       DATA (name(K),K=    1,   76) /
1789      &'d           ','u           ','s           ','c           ',
1790      &'b           ','t           ','(dd)_1      ','(ud)_0      ',
1791      &'(ud)_1      ','(uu)_1      ','(sd)_0      ','(sd)_1      ',
1792      &'(su)_0      ','(su)_1      ','(ss)_1      ','(cd)_0      ',
1793      &'(cd)_1      ','(cu)_0      ','(cu)_1      ','(cs)_0      ',
1794      &'(cs)_1      ','(cc)_1      ','remnant 1   ','remnant 2   ',
1795      &'string      ','mod. string ','coll. string','reggeon     ',
1796      &'pomeron     ','gluon       ','gamma       ','W           ',
1797      &'Z           ','e           ','mu          ','tau         ',
1798      &'nu(e)       ','nu(mu)      ','nu(tau)     ','pi          ',
1799      &'pi          ','eta         ','rho(770)    ','rho(770)    ',
1800      &'ome(782)    ','etap(958)   ','f(0)(980)   ','a(0)(980)   ',
1801      &'a(0)(980)   ','phi(1020)   ','h(1)(1170)  ','b(1)(1235)  ',
1802      &'b(1)(1235)  ','a(1)(1260)  ','a(1)(1260)  ','f(2)(1270)  ',
1803      &'f(1)(1285)  ','eta(1295)   ','pi(1300)    ','pi(1300)    ',
1804      &'a(2)(1320)  ','a(2)(1320)  ','f(1)(1420)  ','ome(1420)   ',
1805      &'rho(1450)   ','rho(1450)   ','f(0)(1500)  ','f(2)p(1525) ',
1806      &'ome(1600)   ','ome(3)(1670)','pi(2)(1670) ','pi(2)(1670) ',
1807      &'phi(1680)   ','rho(3)(1690)','rho(3)(1690)','rho(1700)   '/
1808       DATA (name(K),K=   77,  152) /
1809      &'rho(1700)   ','f(J)(1710)  ','phi(3)(1850)','f(2)(2010)  ',
1810      &'f(4)(2050)  ','f(2)(2300)  ','f(2)(2340)  ','K           ',
1811      &'K           ','K(S)        ','K(L)        ','K*(892)     ',
1812      &'K*(892)     ','K(1)(1270)  ','K(1)(1270)  ','K(1)(1400)  ',
1813      &'K(1)(1400)  ','K*(1410)    ','K*(1410)    ','K(0)*(1430) ',
1814      &'K(0)*(1430) ','K(2)*(1430) ','K(2)*(1430) ','K*(1680)    ',
1815      &'K*(1680)    ','K(2)(1770)  ','K(2)(1770)  ','K(3)*(1780) ',
1816      &'K(3)*(1780) ','K(2)(1820)  ','K(2)(1820)  ','K(4)*(2045) ',
1817      &'K(4)*(2045) ','D           ','D           ','D*(2007)    ',
1818      &'D*(2010)    ','D(1)(2420)  ','D(2)*(2460) ','D(2)*(2460) ',
1819      &'D(s)        ','D(s)*       ','D(s1)(2536) ','B           ',
1820      &'B           ','B*          ','B*          ','B(s)        ',
1821      &'eta(c)(1S)  ','J/psi(1S)   ','chi(c0)(1P) ','chi(c1)(1P) ',
1822      &'chi(c2)(1P) ','psi(2S)     ','psi(3770)   ','psi(4040)   ',
1823      &'psi(4160)   ','psi(4415)   ','Ups(1S)     ','chi(b0)(1P) ',
1824      &'chi(b1)(1P) ','chi(b2)(1P) ','Ups(2S)     ','chi(b0)(2P) ',
1825      &'chi(b1)(2P) ','chi(b2)(2P) ','Ups(3S)     ','Ups(4S)     ',
1826      &'Ups(10860)  ','Ups(11020)  ','p           ','n           ',
1827      &'N(1440)     ','N(1440)     ','N(1520)     ','N(1520)     '/
1828       DATA (name(K),K=  153,  228) /
1829      &'N(1535)     ','N(1535)     ','N(1650)     ','N(1650)     ',
1830      &'N(1675)     ','N(1675)     ','N(1680)     ','N(1680)     ',
1831      &'N(1700)     ','N(1700)     ','N(1710)     ','N(1710)     ',
1832      &'N(1720)     ','N(1720)     ','N(2190)     ','N(2190)     ',
1833      &'Del(1232)   ','Del(1232)   ','Del(1232)   ','Del(1232)   ',
1834      &'Del(1600)   ','Del(1600)   ','Del(1600)   ','Del(1600)   ',
1835      &'Del(1620)   ','Del(1620)   ','Del(1620)   ','Del(1620)   ',
1836      &'Del(1700)   ','Del(1700)   ','Del(1700)   ','Del(1700)   ',
1837      &'Del(1905)   ','Del(1905)   ','Del(1905)   ','Del(1905)   ',
1838      &'Del(1910)   ','Del(1910)   ','Del(1910)   ','Del(1910)   ',
1839      &'Del(1920)   ','Del(1920)   ','Del(1920)   ','Del(1920)   ',
1840      &'Del(1930)   ','Del(1930)   ','Del(1930)   ','Del(1930)   ',
1841      &'Del(1950)   ','Del(1950)   ','Del(1950)   ','Del(1950)   ',
1842      &'Lambda      ','Lam(1405)   ','Lam(1520)   ','Lam(1600)   ',
1843      &'Lam(1670)   ','Lam(1690)   ','Lam(1800)   ','Lam(1810)   ',
1844      &'Lam(1820)   ','Lam(1830)   ','Lam(1890)   ','Lam(2100)   ',
1845      &'Lam(2110)   ','Sigma       ','Sigma       ','Sigma       ',
1846      &'Sig(1385)   ','Sig(1385)   ','Sig(1385)   ','Sig(1660)   ',
1847      &'Sig(1660)   ','Sig(1660)   ','Sig(1670)   ','Sig(1670)   '/
1848       DATA (name(K),K=  229,  260) /
1849      &'Sig(1670)   ','Sig(1750)   ','Sig(1750)   ','Sig(1750)   ',
1850      &'Sig(1775)   ','Sig(1775)   ','Sig(1775)   ','Sig(1915)   ',
1851      &'Sig(1915)   ','Sig(1915)   ','Sig(1940)   ','Sig(1940)   ',
1852      &'Sig(1940)   ','Sig(2030)   ','Sig(2030)   ','Sig(2030)   ',
1853      &'Xi          ','Xi          ','Xi(1530)    ','Xi(1530)    ',
1854      &'Xi(1820)    ','Xi(1820)    ','Omega       ','Lam(c)      ',
1855      &'Lam(c)(2593)','Sig(c)(2455)','Sig(c)(2455)','Sig(c)(2455)',
1856      &'Xi(c)       ','Xi(c)       ','Ome(c)      ','Lam(b)      '/
1857       DATA (ich3(K),K=    1,  260) /
1858      &-1, 2,-1, 2,-1, 2,-2, 1, 1, 4,-2,-2, 1, 1,-2, 1, 1, 4, 4, 1, 1, 4,
1859      & 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0,-3,-3,-3, 0, 0, 0, 3, 0, 0, 0, 3,
1860      & 0, 0, 0, 0, 3, 0, 0, 0, 3, 0, 3, 0, 0, 0, 0, 3, 0, 3, 0, 0, 0, 3,
1861      & 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 3, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 3,
1862      & 0, 0, 3, 0, 3, 0, 3, 0, 3, 3, 0, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 3,
1863      & 0, 0, 3, 0, 0, 3, 3, 3, 3, 3, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1864      & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 3, 0, 3,
1865      & 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3,-3, 0, 3, 6,-3, 0, 3, 6,
1866      &-3, 0, 3, 6,-3, 0, 3, 6,-3, 0, 3, 6,-3, 0, 3, 6,-3, 0, 3, 6,-3, 0,
1867      & 3, 6,-3, 0, 3, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0,-3,
1868      & 3, 0,-3,-3, 0, 3,-3, 0, 3,-3, 0, 3,-3, 0, 3,-3, 0, 3,-3, 0, 3,-3,
1869      & 0, 3, 0,-3, 0,-3,-3, 0,-3, 3, 3, 6, 3, 0, 3, 0, 0, 0/
1870       DATA (iba3(K),K=    1,  260) /
1871      &1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,
1872      &0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1873      &0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1874      &0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1875      &0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
1876      &3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
1877      &3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
1878      &3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3/
1879       DATA (iq_linear(K),K=    1,  418) /
1880      & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 2,
1881      & 1, 0, 2, 1, 0, 2, 2, 0, 3, 1, 0, 3, 1, 0, 3, 2, 0, 3, 2, 0, 3, 3,
1882      & 0, 4, 1, 0, 4, 1, 0, 4, 2, 0, 4, 2, 0, 4, 3, 0, 4, 3, 0, 4, 4, 0,
1883      & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1884      & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1885      & 0, 0, 0, 0, 0, 0, 0, 2,-1, 0, 1,-1, 0, 2,-2, 0, 1,-1, 0, 2,-1, 0,
1886      & 2,-2, 0, 3,-3, 0, 2,-2, 0, 1,-1, 0, 2,-1, 0, 3,-3, 0, 2,-2, 0, 1,
1887      &-1, 0, 2,-1, 0, 1,-1, 0, 2, 1, 0, 2,-2, 0, 2,-2, 0, 2,-2, 0, 1,-1,
1888      & 0, 2,-1, 0, 1,-1, 0, 2,-1, 0, 2,-2, 0, 2,-2, 0, 1,-1, 0, 2,-1, 0,
1889      & 2,-2, 0, 3,-3, 0, 2,-2, 0, 2,-2, 0, 1,-1, 0, 2,-1, 0, 3,-3, 0, 1,
1890      &-1, 0, 2,-1, 0, 1,-1, 0, 2,-1, 0, 2,-2, 0, 3,-3, 0, 2,-2, 0, 2,-2,
1891      & 0, 2,-2, 0, 2,-2, 0, 2,-3, 0, 1,-3, 0, 1,-3, 0, 3,-1, 0, 2,-3, 0,
1892      & 1,-3, 0, 1,-3, 0, 2,-3, 0, 1,-3, 0, 2,-3, 0, 1,-3, 0, 2,-3, 0, 1,
1893      &-3, 0, 2,-3, 0, 2,-3, 0, 1,-3, 0, 1,-3, 0, 2,-3, 0, 1,-3, 0, 2,-3,
1894      & 0, 1,-3, 0, 2,-3, 0, 1,-3, 0, 2,-3, 0, 1,-3, 0, 2,-3, 0, 4,-1, 0,
1895      & 4,-2, 0, 4,-2, 0, 4,-1, 0, 4,-2, 0, 4,-2, 0, 4,-1, 0, 4,-3, 0, 4,
1896      &-3, 0, 4,-3, 0, 2,-5, 0, 1,-5, 0, 1,-5, 0, 2,-5, 0, 3,-5, 0, 4,-4,
1897      & 0, 4,-4, 0, 4,-4, 0, 4,-4, 0, 4,-4, 0, 4,-4, 0, 4,-4, 0, 4,-4, 0,
1898      & 4,-4, 0, 4,-4, 0, 5,-5, 0, 5,-5, 0, 5,-5, 0, 5,-5, 0, 5,-5, 0, 5/
1899       DATA (iq_linear(K),K=  419,  780) /
1900      &-5, 0, 5,-5, 0, 5,-5, 0, 5,-5, 0, 5,-5, 0, 5,-5, 0, 5,-5, 0, 2, 2,
1901      & 1, 2, 1, 1, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 2, 1,
1902      & 2, 1, 1, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2,
1903      & 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1,
1904      & 1, 2, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 2, 2,
1905      & 1, 1, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 2,
1906      & 2, 2, 1, 1, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 1, 2, 1,
1907      & 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 1, 1, 1, 1, 2, 1,
1908      & 2, 1, 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 2, 2, 3, 1, 2, 3,
1909      & 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1,
1910      & 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 2, 2, 3, 2, 1, 3, 1, 1, 3,
1911      & 3, 2, 2, 3, 2, 1, 3, 1, 1, 3, 1, 1, 3, 2, 1, 3, 2, 2, 3, 1, 1, 3,
1912      & 2, 1, 3, 2, 2, 3, 1, 1, 3, 2, 1, 3, 2, 2, 3, 1, 1, 3, 2, 1, 3, 2,
1913      & 2, 3, 1, 1, 3, 2, 1, 3, 2, 2, 3, 1, 1, 3, 2, 1, 3, 2, 2, 3, 1, 1,
1914      & 3, 2, 1, 3, 2, 2, 2, 3, 3, 1, 3, 3, 2, 3, 3, 1, 3, 3, 3, 3, 1, 3,
1915      & 3, 2, 3, 3, 3, 2, 1, 4, 4, 1, 2, 2, 2, 4, 2, 1, 2, 1, 1, 4, 3, 2,
1916      & 2, 3, 1, 2, 3, 3, 4, 5, 1, 2/
1917       DATA (xmass(K),K=    1,  114) /
1918      &3.0000E-01,3.0000E-01,3.5000E-01,1.4500E+00,4.5000E+00,1.7400E+02,
1919      &7.7133E-01,5.7933E-01,7.7133E-01,7.7133E-01,8.0473E-01,9.2953E-01,
1920      &8.0473E-01,9.2953E-01,1.0936E+00,1.9691E+00,2.0081E+00,1.9691E+00,
1921      &2.0081E+00,2.1543E+00,2.1797E+00,3.2753E+00,0.0000E+00,0.0000E+00,
1922      &0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,
1923      &0.0000E+00,8.0410E+01,9.1187E+01,5.1100E-04,1.0566E-01,1.7771E+00,
1924      &0.0000E+00,0.0000E+00,0.0000E+00,1.3957E-01,1.3498E-01,5.4730E-01,
1925      &7.7000E-01,7.7000E-01,7.8194E-01,9.5778E-01,9.8000E-01,9.8340E-01,
1926      &9.8340E-01,1.0194E+00,1.1700E+00,1.2295E+00,1.2295E+00,1.2300E+00,
1927      &1.2300E+00,1.2750E+00,1.2819E+00,1.2970E+00,1.3000E+00,1.3000E+00,
1928      &1.3181E+00,1.3181E+00,1.4262E+00,1.4190E+00,1.4650E+00,1.4650E+00,
1929      &1.5000E+00,1.5250E+00,1.6490E+00,1.6670E+00,1.6700E+00,1.6700E+00,
1930      &1.6800E+00,1.6910E+00,1.6910E+00,1.7000E+00,1.7000E+00,1.7120E+00,
1931      &1.8540E+00,2.0100E+00,2.0440E+00,2.2970E+00,2.3400E+00,4.9368E-01,
1932      &4.9767E-01,4.9767E-01,4.9767E-01,8.9166E-01,8.9610E-01,1.2720E+00,
1933      &1.2720E+00,1.4020E+00,1.4020E+00,1.4140E+00,1.4140E+00,1.4290E+00,
1934      &1.4290E+00,1.4256E+00,1.4324E+00,1.7170E+00,1.7170E+00,1.7730E+00,
1935      &1.7730E+00,1.7760E+00,1.7760E+00,1.8160E+00,1.8160E+00,2.0450E+00,
1936      &2.0450E+00,1.8693E+00,1.8646E+00,2.0067E+00,2.0100E+00,2.4222E+00/
1937       DATA (xmass(K),K=  115,  228) /
1938      &2.4589E+00,2.4590E+00,1.9685E+00,2.1124E+00,2.5353E+00,5.2789E+00,
1939      &5.2792E+00,5.3249E+00,5.3249E+00,5.3693E+00,2.9798E+00,3.0969E+00,
1940      &3.4173E+00,3.5105E+00,3.5562E+00,3.6860E+00,3.7699E+00,4.0400E+00,
1941      &4.1590E+00,4.4150E+00,9.4604E+00,9.8598E+00,9.8919E+00,9.9132E+00,
1942      &1.0023E+01,1.0232E+01,1.0255E+01,1.0268E+01,1.0355E+01,1.0580E+01,
1943      &1.0865E+01,1.1019E+01,9.3827E-01,9.3957E-01,1.4400E+00,1.4400E+00,
1944      &1.5200E+00,1.5200E+00,1.5350E+00,1.5350E+00,1.6500E+00,1.6500E+00,
1945      &1.6750E+00,1.6750E+00,1.6800E+00,1.6800E+00,1.7000E+00,1.7000E+00,
1946      &1.7100E+00,1.7100E+00,1.7200E+00,1.7200E+00,2.1900E+00,2.1900E+00,
1947      &1.2320E+00,1.2320E+00,1.2320E+00,1.2320E+00,1.6000E+00,1.6000E+00,
1948      &1.6000E+00,1.6000E+00,1.6200E+00,1.6200E+00,1.6200E+00,1.6200E+00,
1949      &1.7000E+00,1.7000E+00,1.7000E+00,1.7000E+00,1.9050E+00,1.9050E+00,
1950      &1.9050E+00,1.9050E+00,1.9100E+00,1.9100E+00,1.9100E+00,1.9100E+00,
1951      &1.9200E+00,1.9200E+00,1.9200E+00,1.9200E+00,1.9300E+00,1.9300E+00,
1952      &1.9300E+00,1.9300E+00,1.9500E+00,1.9500E+00,1.9500E+00,1.9500E+00,
1953      &1.1157E+00,1.4070E+00,1.5195E+00,1.6000E+00,1.6700E+00,1.6900E+00,
1954      &1.8000E+00,1.8100E+00,1.8200E+00,1.8300E+00,1.8900E+00,2.1000E+00,
1955      &2.1100E+00,1.1894E+00,1.1926E+00,1.1974E+00,1.3828E+00,1.3837E+00,
1956      &1.3872E+00,1.6600E+00,1.6600E+00,1.6600E+00,1.6700E+00,1.6700E+00/
1957       DATA (xmass(K),K=  229,  260) /
1958      &1.6700E+00,1.7500E+00,1.7500E+00,1.7500E+00,1.7750E+00,1.7750E+00,
1959      &1.7750E+00,1.9150E+00,1.9150E+00,1.9150E+00,1.9400E+00,1.9400E+00,
1960      &1.9400E+00,2.0300E+00,2.0300E+00,2.0300E+00,1.3149E+00,1.3213E+00,
1961      &1.5318E+00,1.5350E+00,1.8230E+00,1.8230E+00,1.6724E+00,2.2849E+00,
1962      &2.5939E+00,2.4528E+00,2.4536E+00,2.4522E+00,2.4656E+00,2.4703E+00,
1963      &2.7040E+00,5.6240E+00/
1964       DATA (gamma(K),K=    1,  114) /
1965      &8.0000E-01,8.0000E-01,8.0000E-01,8.0000E-01,8.0000E-01,8.0000E-01,
1966      &8.0000E-01,8.0000E-01,8.0000E-01,8.0000E-01,8.0000E-01,8.0000E-01,
1967      &8.0000E-01,8.0000E-01,8.0000E-01,0.0000E+00,0.0000E+00,0.0000E+00,
1968      &0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,
1969      &0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,0.0000E+00,
1970      &0.0000E+00,2.0600E+00,2.4900E+00,0.0000E+00,2.9959E-19,2.2700E-12,
1971      &0.0000E+00,0.0000E+00,0.0000E+00,2.5284E-17,7.8000E-09,1.1800E-06,
1972      &1.5070E-01,1.5070E-01,8.4100E-03,2.0300E-04,0.0000E+00,0.0000E+00,
1973      &0.0000E+00,4.4300E-03,3.6000E-01,1.4200E-01,1.4200E-01,0.0000E+00,
1974      &0.0000E+00,1.8550E-01,2.4000E-02,5.3000E-02,0.0000E+00,0.0000E+00,
1975      &1.0700E-01,1.0700E-01,5.5000E-02,1.7000E-01,3.1000E-01,3.1000E-01,
1976      &1.1200E-01,7.6000E-02,2.2000E-01,1.6800E-01,2.5800E-01,2.5800E-01,
1977      &1.5000E-01,1.6000E-01,1.6000E-01,2.4000E-01,2.4000E-01,1.3300E-01,
1978      &8.7000E-02,2.0000E-01,2.0800E-01,1.5000E-01,3.2000E-01,5.3140E-17,
1979      &0.0000E+00,7.3730E-15,1.2730E-17,5.0800E-02,5.0500E-02,9.0000E-02,
1980      &9.0000E-02,1.7400E-01,1.7400E-01,2.3200E-01,2.3200E-01,2.8700E-01,
1981      &2.8700E-01,9.8500E-02,1.0900E-01,3.2000E-01,3.2000E-01,1.8600E-01,
1982      &1.8600E-01,1.5900E-01,1.5900E-01,2.7600E-01,2.7600E-01,1.9800E-01,
1983      &1.9800E-01,6.2300E-13,1.5860E-12,5.0000E-03,2.0000E-03,1.8900E-02/
1984       DATA (gamma(K),K=  115,  228) /
1985      &2.3000E-02,2.5000E-02,1.4100E-12,2.0000E-03,0.0000E+00,3.9900E-13,
1986      &4.2200E-13,0.0000E+00,0.0000E+00,4.2700E-13,1.3200E-02,8.7000E-05,
1987      &1.4000E-02,8.8000E-04,2.0000E-03,2.7700E-04,2.3600E-02,5.2000E-02,
1988      &7.8000E-02,4.3000E-02,5.2500E-05,0.0000E+00,0.0000E+00,0.0000E+00,
1989      &4.4000E-05,0.0000E+00,0.0000E+00,0.0000E+00,2.6300E-05,1.0000E-02,
1990      &1.1000E-01,7.9000E-02,0.0000E+00,7.4240E-28,3.5000E-01,3.5000E-01,
1991      &1.2000E-01,1.2000E-01,1.5000E-01,1.5000E-01,1.5000E-01,1.5000E-01,
1992      &1.5000E-01,1.5000E-01,1.3000E-01,1.3000E-01,1.0000E-01,1.0000E-01,
1993      &1.0000E-01,1.0000E-01,1.5000E-01,1.5000E-01,4.5000E-01,4.5000E-01,
1994      &1.2000E-01,1.2000E-01,1.2000E-01,1.2000E-01,3.5000E-01,3.5000E-01,
1995      &3.5000E-01,3.5000E-01,1.5000E-01,1.5000E-01,1.5000E-01,1.5000E-01,
1996      &3.0000E-01,3.0000E-01,3.0000E-01,3.0000E-01,3.5000E-01,3.5000E-01,
1997      &3.5000E-01,3.5000E-01,2.5000E-01,2.5000E-01,2.5000E-01,2.5000E-01,
1998      &2.0000E-01,2.0000E-01,2.0000E-01,2.0000E-01,3.5000E-01,3.5000E-01,
1999      &3.5000E-01,3.5000E-01,3.0000E-01,3.0000E-01,3.0000E-01,3.0000E-01,
2000      &2.5010E-15,5.0000E-02,1.5600E-02,1.5000E-01,3.5000E-02,6.0000E-02,
2001      &3.0000E-01,1.5000E-01,8.0000E-02,9.5000E-02,1.0000E-01,2.0000E-01,
2002      &2.0000E-01,8.2400E-15,8.9000E-06,4.4500E-15,3.5800E-02,3.6000E-02,
2003      &3.9400E-02,1.0000E-01,1.0000E-01,1.0000E-01,6.0000E-02,6.0000E-02/
2004       DATA (gamma(K),K=  229,  260) /
2005      &6.0000E-02,9.0000E-02,9.0000E-02,9.0000E-02,1.2000E-01,1.2000E-01,
2006      &1.2000E-01,1.2000E-01,1.2000E-01,1.2000E-01,2.2000E-01,2.2000E-01,
2007      &2.2000E-01,1.8000E-01,1.8000E-01,1.8000E-01,2.2700E-15,4.0200E-15,
2008      &9.1000E-03,9.9000E-03,2.4000E-02,2.4000E-02,8.0100E-15,3.1900E-12,
2009      &3.6000E-03,0.0000E+00,0.0000E+00,0.0000E+00,1.8600E-12,6.7000E-12,
2010      &1.0200E-11,5.3100E-13/
2011       DATA (idec_linear(K),K=    1,  304) /
2012      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2013      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2014      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2015      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2016      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2017      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2018      &  0,  0,  0,  0,  0,  0,  3,  1,  1,  2,  2,  6,  0,  0,  0,  0,
2019      &  0,  0,  0,  0,  0,  3,  7,  7,  3,  8,  9,  1, 10, 14,  1, 15,
2020      & 16,  1, 17, 17,  1, 18, 20,  1, 21, 24,  0,  0,  0,  0,  0,  0,
2021      &  0,  0,  0,  1, 25, 29,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2022      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2023      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1, 30, 32,
2024      &  1, 33, 34,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1, 35, 37,  0,
2025      &  0,  0,  0,  0,  0,  0,  0,  0,  1, 38, 39,  0,  0,  0,  0,  0,
2026      &  0,  1, 40, 40,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2027      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  3, 41, 46,  0,  0,  0,  3,
2028      & 47, 48,  3, 49, 52,  1, 53, 54,  1, 55, 56,  1, 57, 58,  1, 59,
2029      & 60,  0,  0,  0,  0,  0,  0,  1, 61, 68,  1, 69, 76,  0,  0,  0,
2030      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0/
2031       DATA (idec_linear(K),K=  305,  608) /
2032      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2033      &  0,  0,  0,  0,  0,  0,  0,  2, 77, 78,  2, 79, 82,  1, 83, 84,
2034      &  1, 85, 87,  0,  0,  0,  0,  0,  0,  0,  0,  0,  2, 88, 90,  1,
2035      & 91, 92,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2036      &  0,  0,  0,  0,  2, 93, 95,  1, 96, 98,  0,  0,  0,  0,  0,  0,
2037      &  0,  0,  0,  1, 99,101,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2038      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2039      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2040      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  3,102,102,  1,103,112,  1,
2041      &113,122,  0,  0,  0,  0,  0,  0,  1,123,129,  1,130,136,  0,  0,
2042      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2043      &  0,  0,  0,  0,  0,  0,  1,137,144,  1,145,152,  0,  0,  0,  0,
2044      &  0,  0,  0,  0,  0,  0,  0,  0,  1,153,153,  1,154,155,  1,156,
2045      &157,  1,158,158,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2046      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,159,162,  1,
2047      &163,169,  1,170,176,  1,177,180,  0,  0,  0,  0,  0,  0,  0,  0,
2048      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2049      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2050      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0/
2051       DATA (idec_linear(K),K=  609,  780) /
2052      &  0,  0,  0,  0,  3,181,182,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2053      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2054      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  3,183,184,  3,185,
2055      &185,  3,186,186,  1,187,189,  1,190,192,  1,193,194,  0,  0,  0,
2056      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2057      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,195,203,  0,  0,
2058      &  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
2059      &  0,  0,  0,  0,  0,  0,  1,204,216,  0,  0,  0,  3,217,217,  3,
2060      &218,218,  1,219,220,  1,221,222,  0,  0,  0,  0,  0,  0,  2,223,
2061      &225,  2,226,239,  0,  0,  0,  2,240,240,  2,241,241,  2,242,242,
2062      &  2,243,246,  2,247,251,  2,252,255,  0,  0,  0/
2063       DATA (isec_linear(K),K=    1,  152) /
2064      &     11,     12,    -12,     13,    -14,     16,     11,    -12,
2065      &     16,   -213,     16,      0,   -211,     16,      0,   -323,
2066      &     16,      0,    -13,     12,      0,     22,     22,      0,
2067      &     22,    -11,     11,     22,     22,      0,    111,     22,
2068      &     22,    111,    111,    111,    211,   -211,    111,    211,
2069      &   -211,     22,    211,   -211,      0,    111,    111,      0,
2070      &    211,    111,      0,    211,   -211,    111,    211,   -211,
2071      &      0,    111,     22,      0,    221,    211,   -211,    221,
2072      &    111,    111,    211,   -211,     22,     22,     22,      0,
2073      &    321,   -321,      0,    130,    310,      0,    113,    111,
2074      &      0,    211,   -211,    111,    221,     22,      0,    113,
2075      &    111,      0,   -213,    211,      0,    213,   -211,      0,
2076      &    211,   -211,      0,    111,    111,      0,    113,    111,
2077      &      0,   -213,    211,      0,    213,   -211,      0,    311,
2078      &   -313,      0,   -311,    313,      0,    113,    211,   -211,
2079      &    -13,     12,      0,    211,    111,      0,    211,    211,
2080      &   -211,    211,    111,    111,    -13,    111,     12,    -11,
2081      &    111,     12,    211,   -211,      0,    111,    111,      0,
2082      &    111,    111,    111,    211,   -211,    111,    211,     13/
2083       DATA (isec_linear(K),K=  153,  304) /
2084      &     12,    211,     11,     12,    321,    111,      0,    311,
2085      &    211,      0,    311,    111,      0,    321,   -211,      0,
2086      &    311,    111,      0,    321,   -211,      0,    321,    111,
2087      &      0,    311,    211,      0,    311,    111,      0,    321,
2088      &   -211,      0,    313,    111,      0,    323,   -211,      0,
2089      &    311,    113,      0,    321,   -213,      0,    311,    223,
2090      &      0,    311,    221,      0,    321,    111,      0,    311,
2091      &    211,      0,    323,    111,      0,    313,    211,      0,
2092      &    321,    113,      0,    311,    213,      0,    321,    223,
2093      &      0,    321,    221,      0,   -321,    211,    211,   -311,
2094      &    211,      0,   -321,    211,      0,   -321,    211,    111,
2095      &    311,    211,   -211,    311,    111,      0,    421,    111,
2096      &      0,    421,     22,      0,    421,    211,      0,    411,
2097      &    111,      0,    411,     22,      0,    221,    211,      0,
2098      &    321,   -321,    321,    321,   -311,      0,    431,     22,
2099      &      0,    431,     22,      0,    111,    111,      0,    211,
2100      &   -211,      0,     22,     22,      0,    -11,     11,      0,
2101      &    -13,     13,      0,    211,   -211,    111,    443,    211,
2102      &   -211,    443,    111,    111,    443,    221,      0,   2212/
2103       DATA (isec_linear(K),K=  305,  456) /
2104      &     11,     12,   2112,    111,      0,   2212,   -211,      0,
2105      &   2112,    111,    111,   2112,    211,   -211,   1114,    211,
2106      &      0,   2114,    111,      0,   2214,   -211,      0,   2112,
2107      &    113,      0,   2212,   -213,      0,   2112,    221,      0,
2108      &   2212,    111,      0,   2112,    211,      0,   2212,    111,
2109      &    111,   2212,    211,   -211,   2224,   -211,      0,   2214,
2110      &    111,      0,   2114,    211,      0,   2212,    113,      0,
2111      &   2112,    213,      0,   2212,    221,      0,   2212,   -211,
2112      &      0,   2112,    111,      0,   2214,   -211,      0,   2114,
2113      &    111,      0,   1114,    211,      0,   2212,   -213,      0,
2114      &   2112,    113,      0,   2212,    111,      0,   2112,    211,
2115      &      0,   2224,   -211,      0,   2214,    111,      0,   2114,
2116      &    211,      0,   2212,    113,      0,   2112,    213,      0,
2117      &   2212,   -211,      0,   2112,    111,      0,   2212,   -213,
2118      &      0,   2112,    113,      0,   3122,    311,      0,   3212,
2119      &    311,      0,   3112,    321,      0,   2112,    221,      0,
2120      &   2212,    111,      0,   2112,    211,      0,   2212,    113,
2121      &      0,   2112,    213,      0,   3122,    321,      0,   3222,
2122      &    311,      0,   3212,    321,      0,   2212,    221,      0/
2123       DATA (isec_linear(K),K=  457,  608) /
2124      &   2112,   -211,      0,   2212,   -211,      0,   2112,    111,
2125      &      0,   2212,    111,      0,   2112,    211,      0,   2212,
2126      &    211,      0,   2112,   -211,      0,   2114,   -211,      0,
2127      &   1114,    111,      0,   2112,   -213,      0,   2212,   -211,
2128      &      0,   2112,    111,      0,   2214,   -211,      0,   2114,
2129      &    111,      0,   1114,    211,      0,   2212,   -213,      0,
2130      &   2112,    113,      0,   2212,    111,      0,   2112,    211,
2131      &      0,   2224,   -211,      0,   2214,    111,      0,   2114,
2132      &    211,      0,   2212,    113,      0,   2112,    213,      0,
2133      &   2212,    211,      0,   2224,    111,      0,   2214,    211,
2134      &      0,   2212,    213,      0,   2212,   -211,      0,   2112,
2135      &    111,      0,   2212,    111,      0,   2112,    211,      0,
2136      &   3122,     22,      0,   2112,   -211,      0,   3122,    211,
2137      &      0,   3212,    211,      0,   3222,    111,      0,   3122,
2138      &    111,      0,   3222,   -211,      0,   3112,    211,      0,
2139      &   3122,   -211,      0,   3212,   -211,      0,   2112,   -311,
2140      &      0,   2212,   -321,      0,   3222,   -211,      0,   3212,
2141      &    111,      0,   3112,    211,      0,   3122,    221,      0,
2142      &   3224,   -211,      0,   3114,    211,      0,   3214,    111/
2143       DATA (isec_linear(K),K=  609,  760) /
2144      &      0,   2112,   -311,      0,   2212,   -321,      0,   3122,
2145      &    111,      0,   3122,    223,      0,   3122,    113,      0,
2146      &   3222,   -213,      0,   3112,    213,      0,   3212,    113,
2147      &      0,   3122,    221,      0,   3212,    221,      0,   3222,
2148      &   -211,      0,   3112,    211,      0,   3212,    111,      0,
2149      &   3122,    111,      0,   3122,   -211,      0,   3322,    111,
2150      &      0,   3312,    211,      0,   3322,   -211,      0,   3312,
2151      &    111,      0,   3322,   -211,      0,   3312,    111,      0,
2152      &   3122,   -321,      0,   3222,    221,      0,   3222,    331,
2153      &      0,   2212,   -311,      0,   3322,    321,      0,   3224,
2154      &    221,      0,   2214,    331,      0,   2224,   -321,      0,
2155      &   3122,    213,      0,   3212,    213,      0,   3222,    113,
2156      &      0,   3222,    223,      0,   2212,   -313,      0,   2214,
2157      &   -313,      0,   2224,   -323,      0,   4122,    211,      0,
2158      &   4122,    111,      0,   4122,   -211,      0,   3222,   -311,
2159      &      0,   3322,    211,      0,   3222,   -313,      0,   3322,
2160      &    213,      0,   3212,   -313,      0,   3222,   -323,      0,
2161      &   3322,    223,      0,   3312,    213,      0,   3214,   -313,
2162      &      0,   3322,   -311,      0,   3322,    313,      0,   3334/
2163       DATA (isec_linear(K),K=  761,  765) /
2164      &    213,      0,   3334,    211,      0/
2165       DATA (wg_chan(K),K=    1,  114) /
2166      &1.0000E+00,2.8000E-01,2.8000E-01,3.5000E-01,7.0000E-02,2.0000E-02,
2167      &1.0000E+00,9.9000E-01,1.0000E-02,3.8000E-01,3.0000E-02,3.0000E-01,
2168      &2.4000E-01,5.0000E-02,1.0000E+00,0.0000E+00,1.0000E+00,8.8800E-01,
2169      &2.5000E-02,8.7000E-02,4.8000E-01,2.4000E-01,2.6000E-01,2.0000E-02,
2170      &4.9100E-01,3.4400E-01,1.2900E-01,2.4000E-02,1.2000E-02,4.0000E-01,
2171      &3.0000E-01,3.0000E-01,6.0000E-01,4.0000E-01,4.0000E-01,3.0000E-01,
2172      &3.0000E-01,5.0000E-01,5.0000E-01,1.0000E+00,6.4000E-01,2.1000E-01,
2173      &6.0000E-02,2.0000E-02,3.0000E-02,4.0000E-02,6.9000E-01,3.1000E-01,
2174      &2.1000E-01,1.2000E-01,2.7000E-01,4.0000E-01,3.3000E-01,6.7000E-01,
2175      &3.3000E-01,6.7000E-01,3.3000E-01,6.7000E-01,3.3000E-01,6.7000E-01,
2176      &1.9000E-01,3.8000E-01,9.0000E-02,2.0000E-01,3.0000E-02,4.0000E-02,
2177      &5.0000E-02,2.0000E-02,1.9000E-01,3.8000E-01,9.0000E-02,2.0000E-01,
2178      &3.0000E-02,4.0000E-02,5.0000E-02,2.0000E-02,7.0000E-01,3.0000E-01,
2179      &1.0000E-01,5.0000E-01,1.6000E-01,2.4000E-01,5.5000E-01,4.5000E-01,
2180      &6.8000E-01,3.0000E-01,2.0000E-02,3.0000E-01,4.0000E-01,3.0000E-01,
2181      &9.0000E-01,1.0000E-01,4.9000E-01,4.9000E-01,2.0000E-02,1.0000E-01,
2182      &1.0000E-01,8.0000E-01,6.0000E-01,3.0000E-01,1.0000E-01,1.0000E+00,
2183      &1.5000E-01,3.5000E-01,7.0000E-02,1.8000E-01,1.1000E-01,6.0000E-02,
2184      &3.0000E-02,1.0000E-02,3.0000E-02,1.0000E-02,1.5000E-01,3.5000E-01/
2185       DATA (wg_chan(K),K=  115,  228) /
2186      &7.0000E-02,1.8000E-01,1.1000E-01,6.0000E-02,3.0000E-02,1.0000E-02,
2187      &3.0000E-02,1.0000E-02,3.7000E-01,1.8000E-01,4.0000E-02,8.0000E-02,
2188      &1.3000E-01,1.3000E-01,7.0000E-02,1.8000E-01,3.7000E-01,1.3000E-01,
2189      &8.0000E-02,4.0000E-02,7.0000E-02,1.3000E-01,1.3000E-01,7.0000E-02,
2190      &4.7000E-01,2.3000E-01,5.0000E-02,1.0000E-02,2.0000E-02,2.0000E-02,
2191      &7.0000E-02,1.3000E-01,2.3000E-01,4.7000E-01,5.0000E-02,2.0000E-02,
2192      &1.0000E-02,2.0000E-02,1.0000E+00,3.3000E-01,6.7000E-01,6.7000E-01,
2193      &3.3000E-01,1.0000E+00,2.5000E-01,1.8000E-01,2.7000E-01,3.0000E-01,
2194      &8.0000E-02,1.7000E-01,2.4000E-01,3.0000E-02,1.8000E-01,1.0000E-01,
2195      &2.0000E-01,1.7000E-01,8.0000E-02,1.8000E-01,3.0000E-02,2.4000E-01,
2196      &2.0000E-01,1.0000E-01,2.5000E-01,2.7000E-01,1.8000E-01,3.0000E-01,
2197      &6.4000E-01,3.6000E-01,5.2000E-01,4.8000E-01,1.0000E+00,1.0000E+00,
2198      &8.8000E-01,6.0000E-02,6.0000E-02,8.8000E-01,6.0000E-02,6.0000E-02,
2199      &8.8000E-01,1.2000E-01,1.9000E-01,1.9000E-01,1.6000E-01,1.6000E-01,
2200      &1.7000E-01,3.0000E-02,3.0000E-02,3.0000E-02,4.0000E-02,1.0000E-01,
2201      &1.0000E-01,2.0000E-01,1.2000E-01,1.0000E-01,4.0000E-02,4.0000E-02,
2202      &5.0000E-02,7.5000E-02,7.5000E-02,3.0000E-02,3.0000E-02,4.0000E-02,
2203      &1.0000E+00,1.0000E+00,3.3000E-01,6.7000E-01,6.7000E-01,3.3000E-01,
2204      &2.5000E-01,2.5000E-01,5.0000E-01,2.0000E-02,3.0000E-02,7.0000E-02/
2205       DATA (wg_chan(K),K=  229,  255) /
2206      &2.0000E-02,2.0000E-02,4.0000E-02,1.3000E-01,7.0000E-02,6.0000E-02,
2207      &6.0000E-02,2.0000E-01,1.4000E-01,4.0000E-02,1.0000E-01,1.0000E+00,
2208      &1.0000E+00,1.0000E+00,2.5000E-01,3.0000E-02,3.0000E-01,4.2000E-01,
2209      &2.2000E-01,3.5000E-01,1.9000E-01,1.6000E-01,8.0000E-02,3.7000E-01,
2210      &2.0000E-01,3.6000E-01,7.0000E-02/
2211       DATA (id_psm_linear(K),K=    1,   36) /
2212      &    111,    211,   -311,    411,      0,      0,   -211,    111,
2213      &   -321,    421,      0,      0,    311,    321,    221,    431,
2214      &      0,      0,   -411,   -421,   -431,    441,      0,      0,
2215      &      0,      0,      0,      0,      0,      0,      0,      0,
2216      &      0,      0,      0,      0/
2217       DATA (id_vem_linear(K),K=    1,   36) /
2218      &    113,    213,   -313,    413,      0,      0,   -213,    113,
2219      &   -323,    423,      0,      0,    313,    323,    333,    433,
2220      &      0,      0,   -413,   -423,   -433,  20443,      0,      0,
2221      &      0,      0,      0,      0,      0,      0,      0,      0,
2222      &      0,      0,      0,      0/
2223       DATA (id_b8_linear(K),K=    1,  171) /
2224      &  1114,  2112,  3112,  4112,     0,     0,  2112,  2212,  3212,
2225      &  4122,     0,     0,  3112,  3212,  3312,  4132,     0,     0,
2226      &  4112,  4122,  4132,  4412,     0,     0,     0,     0,     0,
2227      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2228      &  2112,  2212,  3212,  4122,     0,     0,  2212,  2224,  3222,
2229      &  4222,     0,     0,  3212,  3222,  3322,  4232,     0,     0,
2230      &  4122,  4222,  4232,  4422,     0,     0,     0,     0,     0,
2231      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2232      &  3112,  3212,  3312,  4132,     0,     0,  3212,  3222,  3322,
2233      &  4232,     0,     0,  3312,  3322,  3334,  4332,     0,     0,
2234      &  4132,  4232,  4332,  4432,     0,     0,     0,     0,     0,
2235      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2236      &  4112,  4122,  4132,  4412,     0,     0,  4122,  4222,  4232,
2237      &  4422,     0,     0,  4132,  4232,  4332,  4432,     0,     0,
2238      &  4412,  4422,  4432,  4444,     0,     0,     0,     0,     0,
2239      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2240      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2241      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2242      &     0,     0,     0,     0,     0,     0,     0,     0,     0/
2243       DATA (id_b8_linear(K),K=  172,  216) /
2244      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2245      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2246      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2247      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2248      &     0,     0,     0,     0,     0,     0,     0,     0,     0/
2249       DATA (id_b10_linear(K),K=    1,  171) /
2250      &  1114,  2114,  3114,  4114,     0,     0,  2114,  2214,  3214,
2251      &  4214,     0,     0,  3114,  3214,  3314,  4314,     0,     0,
2252      &  4114,  4214,  4314,  4414,     0,     0,     0,     0,     0,
2253      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2254      &  2114,  2214,  3214,  4214,     0,     0,  2214,  2224,  3224,
2255      &  4224,     0,     0,  3214,  3224,  3324,  4324,     0,     0,
2256      &  4214,  4224,  4324,  4424,     0,     0,     0,     0,     0,
2257      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2258      &  3114,  3214,  3314,  4314,     0,     0,  3214,  3224,  3324,
2259      &  4324,     0,     0,  3314,  3324,  3334,  4334,     0,     0,
2260      &  4314,  4324,  4334,  4434,     0,     0,     0,     0,     0,
2261      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2262      &  4114,  4214,  4314,  4414,     0,     0,  4214,  4224,  4324,
2263      &  4424,     0,     0,  4314,  4324,  4334,  4434,     0,     0,
2264      &  4414,  4424,  4434,  4444,     0,     0,     0,     0,     0,
2265      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2266      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2267      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2268      &     0,     0,     0,     0,     0,     0,     0,     0,     0/
2269       DATA (id_b10_linear(K),K=  172,  216) /
2270      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2271      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2272      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2273      &     0,     0,     0,     0,     0,     0,     0,     0,     0,
2274      &     0,     0,     0,     0,     0,     0,     0,     0,     0/
2275
2276       ID_pdg_max = i_tab_max
2277
2278 C  copy from local to global variables
2279       do i=1,i_tab_max
2280         ID_pdg_list(i) = number(i)
2281         name_list(i)   = name(i)
2282         xm_list(i)     = xmass(i)
2283         gam_list(i)    = gamma(i)
2284         ich3_list(i)   = ich3(i)
2285         iba3_list(i)   = iba3(i)
2286         do j=1,3
2287           iq_list(j,i)   = iq_linear(3*(i-1)+j)
2288           idec_list(j,i) = idec_linear(3*(i-1)+j)
2289         enddo
2290       enddo
2291
2292 C  initialize hash table
2293       call pho_cpcini(ID_pdg_max,ID_pdg_list,ID_list)
2294
2295       itmp = IDEB(71)
2296       IDEB(71) = -1
2297
2298 C  quark index table for mesons
2299       do i=1,6
2300         do j=1,6
2301           id_psm_list(i,j) = ipho_pdg2id(id_psm_linear(6*(j-1)+i))
2302           id_vem_list(i,j) = ipho_pdg2id(id_vem_linear(6*(j-1)+i))
2303         enddo
2304       enddo
2305
2306 C  quark index table for baryons
2307       do i=1,6
2308         do j=1,6
2309           do k=1,6
2310             id_b8_list(i,j,k)  =
2311      &        ipho_pdg2id(id_b8_linear(36*(k-1)+6*(j-1)+i))
2312             id_b10_list(i,j,k) =
2313      &        ipho_pdg2id(id_b10_linear(36*(k-1)+6*(j-1)+i))
2314           enddo
2315         enddo
2316       enddo
2317
2318       IDEB(71) = itmp
2319
2320 C  copy secondary particles
2321 C  (translate PDG-ID to CPC and sort according to CPC)
2322       ichan = 0
2323       do i=1,i_tab_max
2324         if(idec_list(1,i).ne.0) then
2325           do j=idec_list(2,i),idec_list(3,i)
2326             ichan = ichan+1
2327             wg_sec_list(ichan) = wg_chan(j)
2328             do k=1,3
2329               if(isec_linear(3*(j-1)+k).ne.0) then
2330                 isec_list(k,ichan) = ipho_pdg2id(isec_linear(3*(j-1)+k))
2331               else
2332                 isec_list(k,ichan) = 0
2333               endif
2334             enddo
2335           enddo
2336         endif
2337       enddo
2338
2339 C  add two-pion background (low-mass photon dissociation)
2340       i = ipho_pdg2id(92)
2341       ichan = ichan+1
2342       idec_list(1,i) = 1
2343       idec_list(2,i) = ichan
2344       idec_list(3,i) = ichan
2345       wg_sec_list(ichan) = 1.D0
2346       isec_list(1,ichan) = ipho_pdg2id(211)
2347       isec_list(2,ichan) = ipho_pdg2id(-211)
2348       isec_list(3,ichan) = 0
2349
2350 C  min. mass limits for strings: q-qbar
2351       do i=1,6
2352         do j=1,6
2353           AM2P = 1000.D0
2354           AM2V = 1000.D0
2355           do k=1,3
2356 C  pseudo-scalar mesons
2357             i1 = iabs(id_psm_list(i,k))
2358             if(i1.ne.0) then
2359               AM1 = xm_list(i1)
2360             else
2361               AM1 = pho_pmass(i,3)+pho_pmass(k,3)
2362             endif
2363             i2 = iabs(id_psm_list(k,j))
2364             if(i2.ne.0) then
2365               AM2 = xm_list(i2)
2366             else
2367               AM2 = pho_pmass(k,3)+pho_pmass(j,3)
2368             endif
2369             AM2P = MIN(AM2P,AM1+AM2)
2370 C  vector mesons
2371             i1 = iabs(id_vem_list(i,k))
2372             if(i1.ne.0) then
2373               AM1 = xm_list(i1)
2374             else
2375               AM1 = pho_pmass(i,3)+pho_pmass(k,3)
2376             endif
2377             i2 = iabs(id_vem_list(k,j))
2378             if(i2.ne.0) then
2379               AM2 = xm_list(i2)
2380             else
2381               AM2 = pho_pmass(k,3)+pho_pmass(j,3)
2382             endif
2383             AM2V = MIN(AM2V,AM1+AM2)
2384           enddo
2385           xm_psm2_list(i,j) = AM2P
2386           xm_vem2_list(i,j) = AM2V
2387         enddo
2388       enddo
2389
2390 C  min. mass limits for strings: qq-q
2391       do i=1,6
2392         do j=1,6
2393           do k=1,6
2394             AM82  = 1000.D0
2395             AM102 = 1000.D0
2396             do l=1,3
2397 C  pseudo-scalar meson
2398               i1 = iabs(id_psm_list(k,l))
2399               if(i1.ne.0) then
2400                 AM1 = xm_list(i1)
2401               else
2402                 AM1 = pho_pmass(i,3)+pho_pmass(k,3)
2403               endif
2404 C  vector meson
2405               i2 = iabs(id_vem_list(k,l))
2406               if(i2.ne.0) then
2407                 AM2 = xm_list(i2)
2408               else
2409                 AM2 = pho_pmass(i,3)+pho_pmass(k,3)
2410               endif
2411 C  octet baryon
2412               AMM = min(AM1,AM2)
2413               K8  = id_b8_list(i,j,l)
2414               if(K8.ne.0) then
2415                 AM1 = xm_list(K8)
2416               else
2417                 AM1 = pho_pmass(i,3)+pho_pmass(j,3)+pho_pmass(l,3)
2418               endif
2419               AM82  = MIN(AM82, AM1 + AMM)
2420 C  decuplet baryon
2421               K10 = id_b10_list(i,j,l)
2422               if(K10.ne.0) then
2423                 AM2 = xm_list(K10)
2424               else
2425                 AM2 = pho_pmass(i,3)+pho_pmass(j,3)+pho_pmass(l,3)
2426               endif
2427               AM102 = MIN(AM102, AM2 + AMM)
2428             enddo
2429             xm_b82_list(i,j,k)  = AM82
2430             xm_b102_list(i,j,k) = AM102
2431           enddo
2432         enddo
2433       enddo
2434
2435 C  min. mass limits for strings: qq-qbarqbar
2436       do i=1,6
2437         do j=1,6
2438           do ii=1,6
2439             do jj=1,6
2440               AM82  = 1000.D0
2441               AM102 = 1000.D0
2442               do l=1,3
2443 C  octet baryons
2444                 K8  = id_b8_list(i,j,l)
2445                 if(K8.ne.0) then
2446                   AM1 = xm_list(K8)
2447                 else
2448                   AM1 = pho_pmass(i,3)+pho_pmass(j,3)+pho_pmass(l,3)
2449                 endif
2450                 L8  = id_b8_list(ii,jj,l)
2451                 if(L8.ne.0) then
2452                   AM2 = xm_list(L8)
2453                 else
2454                   AM2 = pho_pmass(ii,3)+pho_pmass(jj,3)+pho_pmass(l,3)
2455                 endif
2456                 AM82  = MIN(AM82, AM1+AM2)
2457 C  decuplet baryons
2458                 K10 = id_b10_list(i,j,l)
2459                 if(K10.ne.0) then
2460                   AM1 = xm_list(K10)
2461                 else
2462                   AM1 = pho_pmass(i,3)+pho_pmass(j,3)+pho_pmass(l,3)
2463                 endif
2464                 L10 = id_b10_list(ii,jj,l)
2465                 if(L10.ne.0) then
2466                   AM2 = xm_list(L10)
2467                 else
2468                   AM2 = pho_pmass(ii,3)+pho_pmass(jj,3)+pho_pmass(l,3)
2469                 endif
2470                 AM102 = MIN(AM102, AM1+AM2)
2471               enddo
2472               xm_bb82_list(i,j,ii,jj)  = AM82
2473               xm_bb102_list(i,j,ii,jj) = AM102
2474             enddo
2475           enddo
2476         enddo
2477       enddo
2478
2479       END
2480
2481 *$ CREATE PHO_PRESEL.FOR
2482 *COPY PHO_PRESEL
2483 CDECK  ID>, PHO_PRESEL
2484       SUBROUTINE PHO_PRESEL(MODE,IREJ)
2485 C**********************************************************************
2486 C
2487 C     user specific function to pre-select events during generation
2488 C
2489 C     input:   MODE  5  electron and photon kinematics
2490 C                   10  process and number of cut Pomerons
2491 C                   15  partons without construction of strings
2492 C                   20  partons assigned to strings
2493 C                   25  after fragmentation, complete final state
2494 C
2495 C     output:  IREJ  0  event accepted
2496 C                   50  event rejected
2497 C
2498 C**********************************************************************
2499       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2500       SAVE
2501
2502 C  input/output channels
2503       INTEGER LI,LO
2504       COMMON /POINOU/ LI,LO
2505 C  event debugging information
2506       INTEGER NMAXD
2507       PARAMETER (NMAXD=100)
2508       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2509      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2510       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2511      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2512
2513 C  standard particle data interface
2514       INTEGER NMXHEP
2515
2516       PARAMETER (NMXHEP=4000)
2517
2518       INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
2519       DOUBLE PRECISION PHEP,VHEP
2520       COMMON /POEVT1/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
2521      &                JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),
2522      &                VHEP(4,NMXHEP)
2523 C  extension to standard particle data interface (PHOJET specific)
2524       INTEGER IMPART,IPHIST,ICOLOR
2525       COMMON /POEVT2/ IMPART(NMXHEP),IPHIST(2,NMXHEP),ICOLOR(2,NMXHEP)
2526
2527 C  global event kinematics and particle IDs
2528       INTEGER IFPAP,IFPAB
2529       DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
2530       COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
2531 C  gamma-lepton or gamma-hadron vertex information
2532       INTEGER IGHEL,IDPSRC,IDBSRC
2533       DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
2534      &                 RADSRC,AMSRC,GAMSRC
2535       COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
2536      &                GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
2537      &                IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
2538 C  hard scattering data
2539       INTEGER MSCAHD
2540       PARAMETER ( MSCAHD = 50 )
2541       INTEGER LSCAHD,LSC1HD,LSIDX,
2542      &        NINHD,N0INHD,NIVAL,N0IVAL,NOUTHD,NBRAHD,NPROHD
2543       DOUBLE PRECISION PPH,PTHD,ETAHD,Q2SCA,PDFVA,XHD,VHD,X0HD
2544       COMMON /POHSLT/ LSCAHD,LSC1HD,LSIDX(MSCAHD),
2545      &                PPH(8*MSCAHD,2),PTHD(MSCAHD),ETAHD(MSCAHD,2),
2546      &                Q2SCA(MSCAHD,2),PDFVA(MSCAHD,2),
2547      &                XHD(MSCAHD,2),VHD(MSCAHD),X0HD(MSCAHD,2),
2548      &                NINHD(MSCAHD,2),N0INHD(MSCAHD,2),
2549      &                NIVAL(MSCAHD,2),N0IVAL(MSCAHD,2),
2550      &                NOUTHD(MSCAHD,2),NBRAHD(MSCAHD,2),NPROHD(MSCAHD)
2551 C  event weights and generated cross section
2552       INTEGER IPOWGC,ISWCUT,IVWGHT
2553       DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2554       COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2555      &                IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2556
2557       IREJ = 0
2558
2559 *     XBJ = GQ2(2)/(GGECM**2+GQ2(2))
2560 *     IF(XBJ.LT.0.002D0) IREJ = 1
2561
2562       END
2563
2564 *$ CREATE PHO_FIXCOL.FOR
2565 *COPY PHO_FIXCOL
2566 CDECK  ID>, PHO_FIXCOL
2567       SUBROUTINE PHO_FIXCOL(E1,E2,THETA,PHI,NEV)
2568 C**********************************************************************
2569 C
2570 C     interface to call PHOJET (fixed energy run) with
2571 C     collider kinematics
2572 C
2573 C     equivalen photon approximation to get photon flux
2574 C
2575 C     input:     NEV     number of events to generate
2576 C                THETA   azimuthal angle (micro radians)
2577 C                PHI     beam crossing angle
2578 C                        (with respect to x, in degrees)
2579 C                E1      energy of particle 1 (+z direction, GeV)
2580 C                E2      energy of particle 2 (-z direction, GeV)
2581 C
2582 C     note: particle types have to be specified before
2583 C           with PHO_SETPAR
2584 C
2585 C**********************************************************************
2586       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2587       SAVE
2588
2589       PARAMETER(TWOPI=6.283185307D0,BOG=TWOPI/360.0D0)
2590
2591 C  input/output channels
2592       INTEGER LI,LO
2593       COMMON /POINOU/ LI,LO
2594 C  event debugging information
2595       INTEGER NMAXD
2596       PARAMETER (NMAXD=100)
2597       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2598      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2599       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2600      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2601 C  general process information
2602       INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
2603       COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
2604 C  global event kinematics and particle IDs
2605       INTEGER IFPAP,IFPAB
2606       DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
2607       COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
2608 C  model switches and parameters
2609       CHARACTER*8 MDLNA
2610       INTEGER ISWMDL,IPAMDL
2611       DOUBLE PRECISION PARMDL
2612       COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
2613 C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
2614       INTEGER IDEQP,IDEQB,IHFLD,IHFLS
2615       DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
2616       COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
2617      &                IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
2618 C  integration precision for hard cross sections (obsolete)
2619       INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2620       COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2621 C  event weights and generated cross section
2622       INTEGER IPOWGC,ISWCUT,IVWGHT
2623       DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2624       COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2625      &                IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2626
2627       DIMENSION P1(4),P2(4)
2628
2629 C  remnant initialization (only needed for DPMJET)
2630       ISAVP1 = IFPAP(1)
2631       ISAVB1 = IFPAB(1)
2632       IF(IFPAP(1).EQ.81) THEN
2633         IFPAP(1) = IDEQP(1)
2634         IFPAB(1) = IDEQB(1)
2635       ENDIF
2636       ISAVP2 = IFPAP(2)
2637       ISAVB2 = IFPAB(2)
2638       IF(IFPAP(2).EQ.82) THEN
2639         IFPAP(2) = IDEQP(2)
2640         IFPAB(2) = IDEQB(2)
2641       ENDIF
2642       PMASS1 = PHO_PMASS(IFPAB(1),0)-SQRT(PVIRT(1))
2643       PMASS2 = PHO_PMASS(IFPAB(2),0)-SQRT(PVIRT(2))
2644       PP1 = SQRT(E1**2-PMASS1**2)
2645       PP2 = SQRT(E2**2-PMASS2**2)
2646 C  beam crossing angle
2647       TH = 1.D-6*THETA/2.D0
2648       PH = PHI*BOG
2649       P1(1) = PP1*SIN(TH)*COS(PH)
2650       P1(2) = PP1*SIN(TH)*SIN(PH)
2651       P1(3) = PP1*COS(TH)
2652       P1(4) = E1
2653       P2(1) = PP2*SIN(TH)*COS(PH)
2654       P2(2) = PP2*SIN(TH)*SIN(PH)
2655       P2(3) = -PP2*COS(TH)
2656       P2(4) = E2
2657       CALL PHO_EVENT(-1,P1,P2,SIGMAX,IREJ)
2658       IFPAP(1) = ISAVP1
2659       IFPAB(1) = ISAVB1
2660       IFPAP(2) = ISAVP2
2661       IFPAB(2) = ISAVB2
2662       ITRY = 0
2663       CALL PHO_PHIST(-1,SIGMAX)
2664       CALL PHO_LHIST(-1,SIGMAX)
2665 C  test of DPMJET interface (default is IPAMDL(13)=0)
2666       if(IPAMDL(13).gt.0) then
2667         MODE = IPAMDL(13)
2668         IPAMDL(13) = 0
2669       else
2670         MODE = 1
2671       endif
2672 C  main generation loop
2673       DO 50 I=1,NEV
2674  55     CONTINUE
2675         ITRY = ITRY+1
2676         CALL PHO_EVENT(MODE,P1,P2,SIGCUR,IREJ)
2677         IF(IREJ.NE.0) GOTO 55
2678         CALL PHO_PHIST(1,HSWGHT(0))
2679         CALL PHO_LHIST(1,HSWGHT(0))
2680  50   CONTINUE
2681
2682       IF(NEV.GT.0) THEN
2683         SIGMAX = SIGMAX*DBLE(NEV)/DBLE(ITRY)
2684         WRITE(LO,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
2685      &  '=========================================================',
2686      &  ' *****   simulated cross section: ',SIGMAX,' mb  *****',
2687      &  '========================================================='
2688         CALL PHO_EVENT(-2,P1,P2,SIGCUR,IREJ)
2689         CALL PHO_PHIST(-2,SIGMAX)
2690         CALL PHO_LHIST(-2,SIGMAX)
2691       ELSE
2692         WRITE(LO,'(1X,A,I5)') 'POFCOL: no events simulated',NEV
2693       ENDIF
2694
2695       END
2696
2697 *$ CREATE PHO_FIXLAB.FOR
2698 *COPY PHO_FIXLAB
2699 CDECK  ID>, PHO_FIXLAB
2700       SUBROUTINE PHO_FIXLAB(PLAB,NEV)
2701 C**********************************************************************
2702 C
2703 C     interface to call PHOJET (fixed energy run) with
2704 C     LAB kinematics (second particle as target)
2705 C
2706 C     equivalent photon approximation to get photon flux
2707 C
2708 C     input:     NEV     number of events to generate
2709 C                PLAB    LAB momentum of particle 1
2710 C
2711 C     note: particle types have to be specified before
2712 C           with PHO_SETPAR
2713 C
2714 C**********************************************************************
2715       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2716       SAVE
2717
2718 C  input/output channels
2719       INTEGER LI,LO
2720       COMMON /POINOU/ LI,LO
2721 C  event debugging information
2722       INTEGER NMAXD
2723       PARAMETER (NMAXD=100)
2724       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2725      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2726       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2727      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2728 C  general process information
2729       INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
2730       COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
2731 C  global event kinematics and particle IDs
2732       INTEGER IFPAP,IFPAB
2733       DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
2734       COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
2735 C  model switches and parameters
2736       CHARACTER*8 MDLNA
2737       INTEGER ISWMDL,IPAMDL
2738       DOUBLE PRECISION PARMDL
2739       COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
2740 C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
2741       INTEGER IDEQP,IDEQB,IHFLD,IHFLS
2742       DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
2743       COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
2744      &                IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
2745 C  integration precision for hard cross sections (obsolete)
2746       INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2747       COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2748 C  event weights and generated cross section
2749       INTEGER IPOWGC,ISWCUT,IVWGHT
2750       DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2751       COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2752      &                IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2753
2754       DIMENSION P1(4),P2(4)
2755
2756 C  remnant initialization (only needed for DPMJET)
2757       SPCM = PLAB
2758       ISAVP1 = IFPAP(1)
2759       ISAVB1 = IFPAB(1)
2760       IF(IFPAP(1).EQ.81) THEN
2761         IFPAP(1) = IDEQP(1)
2762         IFPAB(1) = IDEQB(1)
2763       ENDIF
2764       ISAVP2 = IFPAP(2)
2765       ISAVB2 = IFPAB(2)
2766       IF(IFPAP(2).EQ.82) THEN
2767         IFPAP(2) = IDEQP(2)
2768         IFPAB(2) = IDEQB(2)
2769       ENDIF
2770 C  get momenta in LAB system
2771       PMASS1 = PHO_PMASS(IFPAB(1),0)**2-PVIRT(1)
2772       PMASS2 = PHO_PMASS(IFPAB(2),0)**2-PVIRT(2)
2773       IF(PMASS2.LT.0.1D0) THEN
2774         WRITE(LO,'(/1X,2A,2I7)') 'PHO_FIXLAB:ERROR: ',
2775      &    'no LAB system possible',IFPAB(1),IFPAB(2)
2776       ELSE
2777         P1(1) = 0.D0
2778         P1(2) = 0.D0
2779         P1(3) = PLAB
2780         P1(4) = SQRT(PMASS1+PLAB**2)
2781         P2(1) = 0.D0
2782         P2(2) = 0.D0
2783         P2(3) = 0.D0
2784         P2(4) = SQRT(PMASS2)
2785         CALL PHO_EVENT(-1,P1,P2,SIGMAX,IREJ)
2786         IFPAP(1) = ISAVP1
2787         IFPAB(1) = ISAVB1
2788         IFPAP(2) = ISAVP2
2789         IFPAB(2) = ISAVB2
2790         ITRY = 0
2791         CALL PHO_PHIST(-1,SIGMAX)
2792         CALL PHO_LHIST(-1,SIGMAX)
2793 C  event generation loop
2794         DO 40 I=1,NEV
2795  45       CONTINUE
2796           ITRY = ITRY+1
2797           CALL PHO_EVENT(1,P1,P2,SIGCUR,IREJ)
2798           IF(IREJ.NE.0) GOTO 45
2799           CALL PHO_LHIST(1,HSWGHT(0))
2800
2801           CALL PHO_PHIST(10,HSWGHT(0))
2802
2803  40     CONTINUE
2804         IF(NEV.GT.0) THEN
2805           SIGMAX = SIGMAX*DBLE(NEV)/DBLE(ITRY)
2806           WRITE(LO,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
2807      &    '=========================================================',
2808      &    ' *****   simulated cross section: ',SIGMAX,' mb  *****',
2809      &    '========================================================='
2810           CALL PHO_EVENT(-2,P1,P2,SIGCUR,IREJ)
2811           CALL PHO_PHIST(-2,SIGMAX)
2812           CALL PHO_LHIST(-2,SIGMAX)
2813         ELSE
2814           WRITE(LO,'(1X,A,I5)')
2815      &      'PHO_FIXLAB: no events simulated',NEV
2816         ENDIF
2817       ENDIF
2818
2819       END
2820
2821 *$ CREATE PHO_GPHERA.FOR
2822 *COPY PHO_GPHERA
2823 CDECK  ID>, PHO_GPHERA
2824       SUBROUTINE PHO_GPHERA(NEVENT,EE1,EE2)
2825 C**********************************************************************
2826 C
2827 C     interface to call PHOJET (variable energy run) with
2828 C     HERA kinematics, photon as particle 2
2829 C
2830 C     equivalent photon approximation to get photon flux
2831 C
2832 C     input:     NEVENT  number of events to generate
2833 C                EE1     proton energy (LAB system)
2834 C                EE2     electron energy (LAB system)
2835 C             from /POFCUT/:
2836 C                YMIN2    lower limit of Y
2837 C                        (energy fraction taken by photon from electron)
2838 C                YMAX2    upper limit of Y
2839 C                Q2MIN2   lower limit of photon virtuality
2840 C                Q2MAX2   upper limit of photon virtuality
2841 C
2842 C**********************************************************************
2843       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2844       SAVE
2845
2846       PARAMETER ( DEPS = 1.D-10,
2847      &            PI   = 3.14159265359D0 )
2848
2849 C  input/output channels
2850       INTEGER LI,LO
2851       COMMON /POINOU/ LI,LO
2852 C  event debugging information
2853       INTEGER NMAXD
2854       PARAMETER (NMAXD=100)
2855       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2856      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2857       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2858      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2859 C  model switches and parameters
2860       CHARACTER*8 MDLNA
2861       INTEGER ISWMDL,IPAMDL
2862       DOUBLE PRECISION PARMDL
2863       COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
2864 C  photon flux kinematics and cuts
2865       DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
2866      &                 YMIN1,YMAX1,YMIN2,YMAX2,
2867      &                 Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
2868      &                 THMIN1,THMAX1,THMIN2,THMAX2
2869       INTEGER          ITAG1,ITAG2
2870       COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
2871      &                YMIN1,YMAX1,YMIN2,YMAX2,
2872      &                Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
2873      &                THMIN1,THMAX1,THMIN2,THMAX2,
2874      &                ITAG1,ITAG2
2875 C  gamma-lepton or gamma-hadron vertex information
2876       INTEGER IGHEL,IDPSRC,IDBSRC
2877       DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
2878      &                 RADSRC,AMSRC,GAMSRC
2879       COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
2880      &                GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
2881      &                IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
2882 C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
2883       INTEGER IDEQP,IDEQB,IHFLD,IHFLS
2884       DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
2885       COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
2886      &                IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
2887 C  event weights and generated cross section
2888       INTEGER IPOWGC,ISWCUT,IVWGHT
2889       DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2890       COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2891      &                IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2892
2893       DIMENSION P1(4),P2(4)
2894
2895       WRITE(LO,'(//1X,A,I10)') 'PHO_GPHERA: events to process',NEVENT
2896 C  assign particle momenta according to HERA kinematics
2897 C  proton data
2898       PROM = PHO_PMASS(2212,1)
2899       PROM2 = PROM**2
2900       IDPSRC(1) = 0
2901       IDBSRC(1) = 0
2902 C  electron data
2903       ELEM = 0.512D-03
2904       ELEM2 = ELEM**2
2905       AMSRC(2) = ELEM
2906       IDPSRC(2) = 11
2907       IDBSRC(2) = ipho_pdg2id(11)
2908 C
2909       Q2MIN = Q2MIN2
2910       Q2MAX = Q2MAX2
2911 C
2912       XIMAX = LOG(YMAX2)
2913       XIMIN = LOG(YMIN2)
2914       XIDEL = XIMAX-XIMIN
2915 C
2916       IF(Q2MIN.GT.ELEM2*YMIN2**2/(1.D0-YMIN2))
2917      &  WRITE(LO,'(/1X,A,1P2E11.4)')
2918      &  'PHO_GPHERA: lower Q2 cutoff larger than kin. limit:',
2919      &  Q2MIN,ELEM2*YMIN2**2/(1.D0-YMIN2)
2920 C
2921       Max_tab = 50
2922       DELLY = LOG(YMAX2/YMIN2)/DBLE(Max_tab-1)
2923       FLUXT = 0.D0
2924       FLUXL = 0.D0
2925       IF(IDEB(30).GE.1) WRITE(LO,'(1X,A,I5)')
2926      &  'PHO_GPHERA: table of photon flux (trans/long)',Max_tab
2927       DO 100 I=1,Max_tab
2928         Y = EXP(XIMIN+DELLY*DBLE(I-1))
2929         Q2LOW = MAX(Q2MIN,ELEM2*Y**2/(1.D0-Y))
2930         FFT = ((1.D0+(1.D0-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
2931      &         -2.D0*ELEM2*Y*(1.D0/Q2LOW-1.D0/Q2MAX))/(2.D0*PI*137.D0)
2932         FFL = 2.D0*(1.D0-Y)/Y*LOG(Q2MAX/Q2LOW)/(2.D0*PI*137.D0)
2933         FLUXT = FLUXT + Y*FFT
2934         FLUXL = FLUXL + Y*FFL
2935         IF(IDEB(30).GE.1) WRITE(LO,'(5X,1P3E14.4)') Y,FFT,FFL
2936  100  CONTINUE
2937       FLUXT = FLUXT*DELLY
2938       FLUXL = FLUXL*DELLY
2939       IF(IDEB(30).GE.1) WRITE(LO,'(1X,A,1P2E12.4)')
2940      &  'PHO_GPHERA: integrated flux (trans./long.):',FLUXT,FLUXL
2941 C
2942       AY = 0.D0
2943       AY2 = 0.D0
2944       YY = YMIN2
2945       Q2LOW = MAX(Q2MIN,ELEM2*YY**2/(1.D0-YY))
2946       WGMAX = (1.D0+(1.D0-YY)**2)*LOG(Q2MAX/Q2LOW)
2947      &        -2.D0*ELEM2*YY*(1.D0/Q2LOW-1.D0/Q2MAX)*YY
2948       IF(ISWMDL(10).GE.2) WGMAX = WGMAX+2.D0*(1.D0-YY)*LOG(Q2MAX/Q2LOW)
2949 C
2950 C  initialization of PHOJET at upper energy limit
2951 C  proton momentum
2952       P1(1) = 0.D0
2953       P1(2) = 0.D0
2954       P1(3) = SQRT(EE1**2-PROM2+DEPS)
2955       P1(4) = EE1
2956 C  photon momentum
2957       EGAM = YMAX2*EE2
2958       P2(1) = 0.D0
2959       P2(2) = 0.D0
2960       P2(3) = -EGAM
2961       P2(4) = EGAM
2962 C  sum of both photon polarizations
2963       IGHEL(2) = -1
2964 C
2965       CALL PHO_SETPAR(1,2212,0,0.D0)
2966       CALL PHO_SETPAR(2,22,0,0.D0)
2967       CALL PHO_EVENT(-1,P1,P2,SIGMAX,IREJ)
2968       CALL PHO_PHIST(-1,SIGMAX)
2969       CALL PHO_LHIST(-1,SIGMAX)
2970 C
2971 C  generation of events, flux calculation
2972
2973       ECMIN2 = ECMIN**2
2974       ECMAX2 = ECMAX**2
2975       AY = 0.D0
2976       AY2 = 0.D0
2977       Q22MIN = 1.D30
2978       Q22AVE = 0.D0
2979       Q22AV2 = 0.D0
2980       Q22MAX = 0.D0
2981       AN2MIN = 1.D30
2982       AN2MAX = 0.D0
2983       YY2MIN = 1.D30
2984       YY2MAX = 0.D0
2985       NITER = NEVENT
2986       ITRY = 0
2987       ITRW = 0
2988       DO 200 I=1,NITER
2989  150    CONTINUE
2990 C  sample y
2991         ITRY = ITRY+1
2992  175    CONTINUE
2993           ITRW = ITRW+1
2994           YY = EXP(XIDEL*DT_RNDM(AY)+XIMIN)
2995           IF(ISWMDL(10).GE.2) THEN
2996             YEFF = 1.D0+(1.D0-YY)**2+2.D0*(1.D0-YY)
2997           ELSE
2998             YEFF = 1.D0+(1.D0-YY)**2
2999           ENDIF
3000           Q2LOW = MAX(Q2MIN,ELEM2*YY**2/(1.D0-YY))
3001           Q2LOG = LOG(Q2MAX/Q2LOW)
3002           WGH = YEFF*Q2LOG-2.D0*ELEM2*YY**2*(1.D0/Q2LOW-1.D0/Q2MAX)
3003           IF(WGMAX.LT.WGH) THEN
3004             WRITE(LO,'(1X,A,3E12.5)')
3005      &        'PHO_GPHERA: inconsistent weight:',YY,WGMAX,WGH
3006           ENDIF
3007         IF(DT_RNDM(AY2)*WGMAX.GT.WGH) GOTO 175
3008 C  sample Q2
3009         IF(IPAMDL(174).EQ.1) THEN
3010  185      CONTINUE
3011             Q2 = Q2LOW*EXP(Q2LOG*DT_RNDM(YY))
3012             WEIGHT = (YEFF-2.D0*ELEM2*YY**2/Q2)/YEFF
3013           IF(WEIGHT.LT.DT_RNDM(Q2)) GOTO 185
3014         ELSE
3015           Q2 = Q2LOW
3016         ENDIF
3017 C
3018
3019 C  incoming electron
3020         PINI(1,2) = 0.D0
3021         PINI(2,2) = 0.D0
3022         PINI(3,2) = -EE2
3023         PINI(4,2) = EE2
3024         PINI(5,2) = 0.D0
3025 C  outgoing electron
3026         YQ2 = SQRT((1.D0-YY)*Q2)
3027         Q2E = Q2/(4.D0*EE2)
3028         E1Y = EE2*(1.D0-YY)
3029         CALL PHO_SFECFE(SIF,COF)
3030         PFIN(1,2) = YQ2*COF
3031         PFIN(2,2) = YQ2*SIF
3032         PFIN(3,2) = -E1Y+Q2E
3033         PFIN(4,2) = E1Y+Q2E
3034         PFIN(5,2) = 0.D0
3035 C  set /POFSRC/
3036         GYY(2) = YY
3037         GQ2(2) = Q2
3038 C  polar angle
3039         PFTHE(2) = ACOS(PFIN(3,2)/PFIN(4,2))
3040 C  electron tagger
3041         IF(PFIN(4,2).GT.EEMIN2) THEN
3042           IF((PFTHE(2).LT.THMIN2).OR.(PFTHE(2).GT.THMAX2)) GOTO 175
3043         ENDIF
3044 C  azimuthal angle
3045         PFPHI(2) = ATAN2(COF,SIF)
3046 C  photon momentum
3047         P2(1) = -PFIN(1,2)
3048         P2(2) = -PFIN(2,2)
3049         P2(3) = PINI(3,2)-PFIN(3,2)
3050         P2(4) = PINI(4,2)-PFIN(4,2)
3051 C  proton momentum
3052         P1(1) = 0.D0
3053         P1(2) = 0.D0
3054         P1(3) = SQRT(EE1**2-PROM2)
3055         P1(4) = EE1
3056 C  ECMS cut
3057         GGECM = (P1(4)+P2(4))**2-(P1(1)+P2(1))**2
3058      &         -(P1(2)+P2(2))**2-(P1(3)+P2(3))**2
3059         IF((GGECM.LT.ECMIN2).OR.(GGECM.GT.ECMAX2)) GOTO 175
3060         GGECM = SQRT(GGECM)
3061 C
3062         PGAM(1,2) = P2(1)
3063         PGAM(2,2) = P2(2)
3064         PGAM(3,2) = P2(3)
3065         PGAM(4,2) = P2(4)
3066         PGAM(5,2) = -SQRT(Q2)
3067 C  photon helicity
3068         IF(ISWMDL(10).GE.2) THEN
3069           WGH  = YEFF-2.D0*ELEM2*YY**2/Q2
3070           WGHL = 2.D0*(1-YY)
3071           IF(DT_RNDM(YY).GE.WGHL/WGH) THEN
3072             IGHEL(2) = 1
3073           ELSE
3074             IGHEL(2) = 0
3075           ENDIF
3076         ELSE
3077           IGHEL(2) = -1
3078         ENDIF
3079 C  user cuts
3080         CALL PHO_PRESEL(5,IREJ)
3081         IF(IREJ.NE.0) GOTO 175
3082 C  event generation
3083         CALL PHO_EVENT(1,P1,P2,SIGCUR,IREJ)
3084         IF(IREJ.NE.0) GOTO 150
3085
3086 C  statistics
3087         AY = AY+YY
3088         AY2 = AY2+YY*YY
3089         YY2MIN = MIN(YY2MIN,YY)
3090         YY2MAX = MAX(YY2MAX,YY)
3091         Q22MIN = MIN(Q22MIN,Q2)
3092         Q22MAX = MAX(Q22MAX,Q2)
3093         Q22AVE = Q22AVE+Q2
3094         Q22AV2 = Q22AV2+Q2*Q2
3095         AN2MIN = MIN(AN2MIN,PFTHE(2))
3096         AN2MAX = MAX(AN2MAX,PFTHE(2))
3097 C  histograms
3098         CALL PHO_PHIST(1,HSWGHT(0))
3099         CALL PHO_LHIST(1,HSWGHT(0))
3100  200  CONTINUE
3101 C
3102       WGY = WGMAX*DBLE(ITRY)/DBLE(ITRW)/(137.D0*2.D0*PI)
3103       WGY = WGY*LOG(YMAX2/YMIN2)
3104       AY  = AY/DBLE(NITER)
3105       AY2 = AY2/DBLE(NITER)
3106       DAY = SQRT((AY2-AY**2)/DBLE(NITER))
3107       Q22AVE = Q22AVE/DBLE(NITER)
3108       Q22AV2 = Q22AV2/DBLE(NITER)
3109       Q22AV2 = SQRT((Q22AV2-Q22AVE**2)/DBLE(NITER))
3110       WEIGHT = WGY*SIGMAX*DBLE(NITER)/DBLE(ITRY)
3111 C  output of histograms
3112       WRITE(LO,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
3113      &'=========================================================',
3114      &' *****   simulated cross section: ',WEIGHT,' mb  *****',
3115      &'========================================================='
3116       WRITE(LO,'(//1X,A,3I10)')
3117      &  'PHO_GPHERA:SUMMARY:NITER,ITRY,ITRW',NITER,ITRY,ITRW
3118       WRITE(LO,'(1X,A,1P2E12.4)') 'EFFECTIVE WEIGHT (FLUX,TOTAL)',
3119      &  WGY,WEIGHT
3120       WRITE(LO,'(1X,A,1P2E12.4)') 'AVERAGE Y,DY                 ',AY,DAY
3121       WRITE(LO,'(1X,A,1P2E12.4)') 'SAMPLED Y RANGE PHOTON       ',
3122      &  YY2MIN,YY2MAX
3123       WRITE(LO,'(1X,A,1P2E12.4)') 'AVERAGE Q2,DQ2               ',
3124      &  Q22AVE,Q22AV2
3125       WRITE(LO,'(1X,A,1P2E12.4)') 'SAMPLED Q2 RANGE PHOTON      ',
3126      &  Q22MIN,Q22MAX
3127       WRITE(LO,'(1X,A,1P4E12.4)') 'SAMPLED THETA RANGE ELECTRON ',
3128      &  AN2MIN,AN2MAX,PI-AN2MAX,PI-AN2MIN
3129 C
3130       CALL PHO_EVENT(-2,P1,P2,WEIGHT,IREJ)
3131       IF(NITER.GT.1) THEN
3132         CALL PHO_PHIST(-2,WEIGHT)
3133         CALL PHO_LHIST(-2,WEIGHT)
3134       ELSE
3135         WRITE(LO,'(1X,A,I4)') 'PHO_GPHERA:NO OUTPUT OF HISTOGRAMS',NITER
3136       ENDIF
3137
3138       END
3139
3140 *$ CREATE PHO_GGEPEM.FOR
3141 *COPY PHO_GGEPEM
3142 CDECK  ID>, PHO_GGEPEM
3143       SUBROUTINE PHO_GGEPEM(NEVENT,EE1,EE2)
3144 C**********************************************************************
3145 C
3146 C     interface to call PHOJET (variable energy run) for
3147 C     gamma-gamma collisions on e+e- collider
3148 C
3149 C     fully differential equivalent (improved) photon approximation
3150 C     to get photon flux
3151 C
3152 C     input:     EE1     LAB system energy of electron/positron 1
3153 C                EE2     LAB system energy of electron/positron 2
3154 C                NEVENT  >0  number of events to generate
3155 C                        -1   initialization
3156 C                        -2   final call (cross section calculation)
3157 C            from /LEPCUT/:
3158 C                YMIN1   lower limit of Y1
3159 C                        (energy fraction taken by photon from electron)
3160 C                YMAX1   upper limit of Y1
3161 C                Q2MIN1  lower limit of photon virtuality
3162 C                Q2MAX1  upper limit of photon virtuality
3163 C                THMIN1  lower limit of scattered electron
3164 C                THMAX1  upper limit of scattered electron
3165 C                YMIN2   lower limit of Y2
3166 C                        (energy fraction taken by photon from electron)
3167 C                YMAX2   upper limit of Y2
3168 C                Q2MIN2  lower limit of photon virtuality
3169 C                Q2MAX2  upper limit of photon virtuality
3170 C                THMIN2  lower limit of scattered electron
3171 C                THMAX2  upper limit of scattered electron
3172 C
3173 C     output:    after final call with NEVENT=-2
3174 C                EE1     e+ e- cross section (mb)
3175 C                EE2     gamma-gamma cross section (mb)
3176 C
3177 C**********************************************************************
3178
3179       IMPLICIT NONE
3180
3181       SAVE
3182
3183       DOUBLE PRECISION EE1,EE2
3184       INTEGER NEVENT
3185
3186 C  input/output channels
3187       INTEGER LI,LO
3188       COMMON /POINOU/ LI,LO
3189 C  event debugging information
3190       INTEGER NMAXD
3191       PARAMETER (NMAXD=100)
3192       INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
3193      &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3194       COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
3195      &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3196 C  model switches and parameters
3197       CHARACTER*8 MDLNA
3198       INTEGER ISWMDL,IPAMDL
3199       DOUBLE PRECISION PARMDL
3200       COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
3201 C  some constants
3202       DOUBLE PRECISION PI,PI2,PI4,GEV2MB,Q_ch,Q_ch2,Q_ch4
3203       COMMON /POCONS/ PI,PI2,PI4,GEV2MB,
3204      &  Q_ch(-6:6),Q_ch2(-6:6),Q_ch4(-6:6)
3205 C  photon flux kinematics and cuts
3206       DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
3207      &                 YMIN1,YMAX1,YMIN2,YMAX2,
3208      &                 Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
3209      &                 THMIN1,THMAX1,THMIN2,THMAX2
3210       INTEGER          ITAG1,ITAG2
3211       COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
3212      &                YMIN1,YMAX1,YMIN2,YMAX2,
3213      &                Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
3214      &                THMIN1,THMAX1,THMIN2,THMAX2,
3215      &                ITAG1,ITAG2
3216 C  gamma-lepton or gamma-hadron vertex information
3217       INTEGER IGHEL,IDPSRC,IDBSRC
3218       DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
3219      &                 RADSRC,AMSRC,GAMSRC
3220       COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
3221      &                GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
3222      &                IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
3223 C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
3224       INTEGER IDEQP,IDEQB,IHFLD,IHFLS
3225       DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
3226       COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
3227      &                IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
3228 C  event weights and generated cross section
3229       INTEGER IPOWGC,ISWCUT,IVWGHT
3230       DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
3231       COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
3232      &                IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
3233
3234 C  external functions
3235       DOUBLE PRECISION DT_RNDM
3236
3237 C  local variables
3238       DOUBLE PRECISION AN1MAX,AN1MIN,AN2MAX,AN2MIN,AY1,AY2,AYS1,AYS2,
3239      &  COF1,COF2,CPFTHE,DAY1,DAY2,DELLY,DITRY,DITRW,
3240      &  ECFRAC,ECMAX2,ECMIN2,EGAM,ELEM,ELEM2,FFL,FFT,FLUXL,FLUXT,
3241      &  FLXAPP,FLXQPM,GGECM2,P1,P2,PP,PT,PT2,Q21AV2,Q21AVE,Q21MAX,
3242      &  Q21MIN,Q22AV2,Q22AVE,Q22MAX,Q22MIN,Q2LOG1,Q2LOG2,Q2LOW1,
3243      &  Q2LOW2,Q2P1,Q2P2,SIF1,SIF2,SIGCUR,SIGMAX,THMAC1,
3244      &  THMAC2,THMIC1,THMIC2,WEIGHT,WG,WGFX,WGH,WGHAPP,WGHL,WGHQPM,
3245      &  WGMAX,WGY,X1DEL,X1MAX,X1MIN,X2DEL,X2MAX,X2MIN,Y1,Y2,YEFF1,YEFF2,
3246      &  YMI,YY1MAX,YY1MIN,YY2MAX,YY2MIN
3247
3248       INTEGER I,IHEAC1,IHEAC2,IHETRY,IREJ,ITRW_low,ITRW_high,ITRY_low,
3249      &  ITRY_high,K,Max_tab,NITER,ITG1,ITG2
3250
3251       DIMENSION P1(4),P2(4),IHETRY(4),IHEAC1(4),IHEAC2(4)
3252       integer ipho_pdg2id
3253
3254 C  initialization of event generation
3255
3256       if(NEVENT.eq.-1) then
3257
3258         DO 10 I=1,4
3259           IHETRY(I) = 0
3260           IHEAC1(I) = 0
3261           IHEAC2(I) = 0
3262  10     CONTINUE
3263
3264         WRITE(LO,'(//1X,A)') 'PHO_GGEPEM: initialization'
3265
3266 C  electron data
3267         ELEM = 0.512D-03
3268         ELEM2 = ELEM**2
3269         AMSRC(1) = ELEM
3270         AMSRC(2) = ELEM
3271 C  lepton numbers
3272         IDPSRC(1) = 11
3273         IDPSRC(2) = -11
3274         IDBSRC(1) = ipho_pdg2id(11)
3275         IDBSRC(2) = ipho_pdg2id(-11)
3276
3277 C  check/update kinematic limitations
3278
3279         Ymi = min(Ymax1,1.D0-ELEM/EE1)
3280         if(Ymi.lt.Ymax1) then
3281           WRITE(LO,'(/1X,A,2E12.5)')
3282      &      'PHO_GGEPEM: Ymax1 decreased (old/new)',Ymax1,Ymi
3283           Ymax1 = YMI
3284         endif
3285         Ymi = min(Ymax2,1.D0-ELEM/EE2)
3286         if(Ymi.lt.Ymax2) then
3287           WRITE(LO,'(/1X,A,2E12.5)')
3288      &      'PHO_GGEPEM: Ymax2 decreased (old/new)',Ymax2,Ymi
3289           Ymax2 = YMI
3290         endif
3291
3292         YMI = ECMIN**2/(4.D0*EE1*EE2*YMAX2)
3293         IF(YMIN1.LT.YMI) THEN
3294           WRITE(LO,'(/1X,A,2E12.5)')
3295      &      'PHO_GGEPEM: Ymin1 increased (old/new)',YMIN1,YMI
3296           YMIN1 = YMI
3297         ELSE IF(YMIN1.GT.YMI) THEN
3298           WRITE(LO,'(/1X,A,/1X,A,E12.5,A,E12.5)')
3299      &      'PHO_GGEPEM:','ECM-CUT&