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)')