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