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