Double check if SM is running added. Some redundant output removed from SM
[u/mrichter/AliRoot.git] / DPMJET / phojet1.12-35c4.f
CommitLineData
7b076c76 1C***********************************************************************
2C
3C
4C
5C PHOJET version 1.12
6C -------------------
7C
8C
9C ($Revision: 1.12.1.35 $, $Date: 2000/06/25 21:59:19 $)
10C
11C
12C Authors: Ralph Engel
13C (ralph.engel@fzk.de)
14C
15C Johannes Ranft
16C (johannes.ranft@cern.ch)
17C
18C Stefan Roesler
19C (Stefan.Roesler@cern.ch)
20C
21C
22C For the latest version and documentation check
23C http://www-ik.fzk.de/~engel/phojet.html
24C
25C
26C Bug reports, questions, complaints are welcome
27C (please send a mail to ralph.engel@fzk.de).
28C
29C
30C Note that the code is available with several interfaces to
31C Lund fragmentation programs (JETSET7.x, 1.x and a double
32C precision JETSET version). This file is the code with
33C
34
35C interface to PYTHIA 6.1 (or higher)
36
37C for usage in DPMJET 3.x (Lund common block dimensions increased)
38
39C
40C***********************************************************************
41C
42C
43C List of subroutines and functions
44C ---------------------------------
45C
46C
47C main event simulation routines
48C
49C PHO_EVENT
50C PHO_PARTON
51C PHO_POSPOM
52C
53C PHO_STDPAR
54C PHO_POMSCA
55C
56C
57C user steering interface
58C
59C PHO_SETMDL
60C PHO_PRESEL
61C
62C
63C experimental setup / photon flux calculation
64C
65C PHO_FIXLAB
66C PHO_FIXCOL
67C PHO_GPHERA
68C PHO_GGEPEM
69C PHO_WGEPEM
70C PHO_GGBLSR
71C PHO_GGBEAM
72C PHO_GGHIOF
73C PHO_GGHIOG
74C PHO_GGFLCL
75C PHO_GGFLCR
76C PHO_GGFAUX
77C PHO_GGFNUC
78C PHO_GHHIOF
79C PHO_GHHIAS
80C
81C
82C initialization
83C
84C PHO_INIT
85C PHO_DATINI
86C PHO_PARDAT
87C PHO_MCINI
88C
89C PHO_EVEINI
90C
91C PHO_HARINI
92C PHO_FRAINI
93C
94C PHO_FITPAR
95C
96C
97C cross section calculation
98C
99C PHO_CSINT
100C
101C PHO_XSECT
102C PHO_BORNCS
103C PHO_HARXTO
104C
105C PHO_DSIGDT
106C
107C PHO_TRIREG
108C PHO_LOOREG
109C PHO_TRXPOM
110C
111C PHO_EIKON
112C PHO_CHAN2A
113C
114C PHO_SCALES
115C
116C
117C multiple interaction structure
118C
119C PHO_IMPAMP
120C PHO_PRBDIS
121C PHO_SAMPRO
122C PHO_SAMPRB
123C
124C
125C hadron / photon remnant treatment, soft x selection
126C
127C PHO_HARREM
128C PHO_PARREM
129C
130C PHO_HADSP2
131C PHO_HADSP3
132C PHO_SOFTXX
133C PHO_SELSXR
134C PHO_SELSX2
135C PHO_SELSXS
136C PHO_SELSXI
137C
138C PHO_VALFLA
139C PHO_REGFLA
140C PHO_SEAFLA
141C PHO_FLAUX
142C PHO_BETAF
143C IPHO_DIQU
144C
145C
146C primordial kt and soft parton pt
147C
148C PHO_PRIMKT
149C PHO_PARTPT
150C PHO_SOFTPT
151C PHO_SELPT
152C
153C PHO_CONN0
154C PHO_CONN1
155C
156C
157C simulation of hard scattering, initial state radiation
158C
159C PHO_HARCOL
160C PHO_SELCOL
161C PHO_HARCOR
162C
163C PHO_HARDIR
164C PHO_HARX12
165C PHO_HARDX1
166C PHO_HARKIN
167C PHO_HARWGH
168C PHO_HARSCA
169C PHO_HARFAC
170C PHO_HARWGX
171C PHO_HARWGI
172C PHO_HARINT
173C PHO_HARMCI
174C
175C PHO_HARXR3
176C PHO_HARXR2
177C PHO_HARXD2
178C PHO_HARXPT
179C PHO_HARISR
180C PHO_HARZSP
181C
182C PHO_PTCUT
183C PHO_ALPHAE
184C PHO_ALPHAS
185C
186C
187C diffraction dissociation
188C
189C PHO_DIFDIS
190C PHO_DIFPRO
191C PHO_DIFPAR
192C PHO_QELAST
193C PHO_CDIFF
194C PHO_DFWRAP
195C
196C PHO_SAMASS
197C PHO_DSIGDM
198C PHO_DFMASS
199C
200C PHO_SDECAY
201C PHO_SDECY2
202C PHO_SDECY3
203C
204C PHO_DIFSLP
205C PHO_DIFKIN
206C PHO_VECRES
207C PHO_DIFRES
208C
209C PHO_REGPAR
210C
211C PHO_PECMS
212C PHO_SETPAR
213C
214C
215C fragmentation, treatment of low-mass strings
216C
217C PHO_STRING
218C PHO_STRFRA
219C
220C PHO_ID2STR
221C PHO_MCHECK
222C PHO_POMCOR
223C PHO_MASCOR
224C PHO_PARCOR
225C
226C PHO_GLU2QU
227C PHO_GLUSPL
228C
229C PHO_DQMASS
230C PHO_BAMASS
231C PHO_MEMASS
232C
233C
234C particle code tables, particle numbering conversion
235C
236C PHO_PNAME
237C PHO_PMASS
238C IPHO_CHR3
239C IPHO_BAR3
240C
241C IPHO_ANTI
242C
243C IPHO_PDG2ID
244C IPHO_ID2PDG
245C IPHO_LU2PDG
246C IPHO_PDG2LU
247C
248C IPHO_CNV1
249C PHO_HACODE
250C
251C
252C
253C Lorentz transformations, rotations and mass adjustment
254C
255C PHO_ALTRA
256C PHO_LTRANS
257C PHO_TRANS
258C PHO_TRANI
259C
260C PHO_MKSLTR
261C PHO_GETLTR
262C
263C PHO_LTRHEP
264C
265C PHO_MSHELL
266C PHO_MASSAD
267C
268C
269C program debugging and internal cross-checks
270C
271C PHO_PREVNT
272C PHO_PRSTRG
273C PHO_CHECK
274C
275C PHO_TRACE
276C
277C PHO_REJSTA
278C
279C PHO_ABORT
280C
281C
282C cross section fitting
283C
284C PHO_FITMAI
285C PHO_FITINP
286C PHO_FITDAT
287C PHO_FITOUT
288C PHO_FITAMP
289C PHO_FITTST
290C PHO_FITMSQ
291C PHO_FITVD1
292C PHO_FITCN1
293C PHO_FITINI
294C
295C
296C cross section parametrizations
297C
298C PHO_HADCSL
299C PHO_ALLM97
300C PHO_CSDIFF
301C
302
303C
304C random numbers
305C
306
307C DPMJET random number generator DT_RNDM used
308
309C
310C PHO_SFECFE
311C PHO_RNDBET
312C PHO_RNDGAM
313C
314C
315C auxiliary routines / numerical methods
316C
317C PHO_GAUSET
318C PHO_GAUDAT
319C
320C pho_samp1d
321C
322C PHO_DZEROX
323C PHO_EXPINT
324C PHO_BESSJ0
325C PHO_BESSI0
326C pho_ExpBessI0
327C PHO_BESSI1
328C PHO_BESSK0
329C PHO_BESSK1
330C
331C PHO_XLAM
332C
333C PHO_SWAPD
334C PHO_SWAPI
335C
336C
337C parton density parametrization management / interface
338C
339C PHO_PDF
340C
341C PHO_SETPDF
342C PHO_GETPDF
343C PHO_ACTPDF
344C
345C PHO_QPMPDF
346C
347C PHO_PDFTST
348C
349C
350C parton density parametrizations from other authors
351C
352C PHO_DOR98LO
353C PHO_DOR98SC
354C PHO_DOR94LO
355C PHO_DOR94HO
356C PHO_DOR94DI
357C PHO_DOR92LO
358C PHO_DOR92HO
359C PHO_DORPLO
360C PHO_DORPHO
361C PHO_DORGLO
362C PHO_DORGHO
363C PHO_DORGH0
364C PHO_DOR94FV
365C PHO_DOR94FW
366C PHO_DOR94FS
367C PHO_DOR92FV
368C PHO_DOR92FW
369C PHO_DOR92FS
370C PHO_DORFVP
371C PHO_DORFGP
372C PHO_DORFQP
373C PHO_DORGF
374C PHO_DORGFS
375C PHO_grsf1
376C PHO_grsf2
377C
378C PHO_CKMTPA
379C PHO_CKMTPD
380C PHO_CKMTPO
381C PHO_CKMTFV
382C
383C PHO_DBFINT
384C
385C PHO_SASGAM
386C PHO_SASVMD
387C PHO_SASANO
388C PHO_SASBEH
389C PHO_SASDIR
390C
391C PHO_PHGAL
392C PHVAL
393C
394C
395C***********************************************************************
396
397*$ CREATE PHO_INIT.FOR
398*COPY PHO_INIT
399CDECK ID>, PHO_INIT
400 SUBROUTINE PHO_INIT(LINP,LOUT,IREJ)
401C***********************************************************************
402C
403C main subroutine to configure and manage PHOJET calculations
404C
405C input: LINP input unit to read from
406C -1 to skip reading of input file
407C LOUT output unit to write to
408C
409C output: IREJ 0 success
410C 1 failure
411C
412C***********************************************************************
413 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
414 SAVE
415
416C input/output channels
417 INTEGER LI,LO
418 COMMON /POINOU/ LI,LO
419C event debugging information
420 INTEGER NMAXD
421 PARAMETER (NMAXD=100)
422 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
423 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
424 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
425 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
426C model switches and parameters
427 CHARACTER*8 MDLNA
428 INTEGER ISWMDL,IPAMDL
429 DOUBLE PRECISION PARMDL
430 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
431C general process information
432 INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
433 COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
434
435C global event kinematics and particle IDs
436 INTEGER IFPAP,IFPAB
437 DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
438 COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
439C nucleon-nucleus / nucleus-nucleus interface to DPMJET
440 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
441 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
442 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
443 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
444C integration precision for hard cross sections (obsolete)
445 INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
446 COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
447C some hadron information, will be deleted in future versions
448 INTEGER NFS
449 DOUBLE PRECISION QMASS,BET,PCOUDI,PNORM,VALPRG
450 COMMON /POHDRN/ QMASS(6),BET,PCOUDI,PNORM,VALPRG(2),NFS
451C obsolete cut-off information
452 DOUBLE PRECISION PTCUT,PTANO,FPS,FPH,PSOMIN,XSOMIN
453 COMMON /POCUT1/ PTCUT(4),PTANO(4),FPS(4),FPH(4),PSOMIN,XSOMIN
454C photon flux kinematics and cuts
455 DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
456 & YMIN1,YMAX1,YMIN2,YMAX2,
457 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
458 & THMIN1,THMAX1,THMIN2,THMAX2
459 INTEGER ITAG1,ITAG2
460 COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
461 & YMIN1,YMAX1,YMIN2,YMAX2,
462 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
463 & THMIN1,THMAX1,THMIN2,THMAX2,
464 & ITAG1,ITAG2
465C cut probability distribution
466 INTEGER IEETA1,IIMAX,KKMAX
467 PARAMETER( IEETA1=20, IIMAX=20, KKMAX=20 )
468 INTEGER IEEMAX,IMAX,KMAX
469 REAL PROB
470 DOUBLE PRECISION EPTAB
471 COMMON /POPROB/ PROB(4,IEETA1,0:IIMAX,0:KKMAX),EPTAB(4,IEETA1),
472 & IEEMAX,IMAX,KMAX
473C event weights and generated cross section
474 INTEGER IPOWGC,ISWCUT,IVWGHT
475 DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
476 COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
477 & IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
478C names of hard scattering processes
479 INTEGER Max_pro_1
480 PARAMETER ( Max_pro_1 = 16 )
481 CHARACTER*18 PROC
482 COMMON /POHPRO/ PROC(0:Max_pro_1)
483C hard cross sections and MC selection weights
484 INTEGER Max_pro_2
485 PARAMETER ( Max_pro_2 = 16 )
486 INTEGER IHa_last,IHb_last,MH_pro_on,MH_tried,
487 & MH_acc_1,MH_acc_2
488 DOUBLE PRECISION Hfac,HWgx,HSig,Hdpt,HEcm_last,HQ2a_last,HQ2b_last
489 COMMON /POHRCS/ Hfac(-1:Max_pro_2),HWgx(-1:Max_pro_2),
490 & HSig(-1:Max_pro_2),Hdpt(-1:Max_pro_2),
491 & HEcm_last,HQ2a_last,HQ2b_last,IHa_last,IHb_last,
492 & MH_pro_on(-1:Max_pro_2,0:4),MH_tried(-1:Max_pro_2,0:4),
493 & MH_acc_1(-1:Max_pro_2,0:4),MH_acc_2(-1:Max_pro_2,0:4)
494
495 INTEGER MSTU,MSTJ
496 DOUBLE PRECISION PARU,PARJ
497 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
498
499 INTEGER KCHG
500 DOUBLE PRECISION PMAS,PARF,VCKM
501 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
502
503 INTEGER MDCY,MDME,KFDP
504 DOUBLE PRECISION BRAT
a374771e 505 COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
7b076c76 506
507 INTEGER PYCOMP
508
509 DIMENSION ITMP(0:11)
510 CHARACTER*10 CNAME
511 CHARACTER*70 NUMBER,FILENA
512
513 14 FORMAT(A10,A69)
514 15 FORMAT(A12)
515
516C define input/output units
517 IF(LINP.GE.0) THEN
518 LI = LINP
519 ELSE
520 LI = 5
521 ENDIF
522 LO = LOUT
523
524 IREJ = 0
525
526 WRITE(LO,*)
527 WRITE(LO,*) ' ==================================================='
528 WRITE(LO,*) ' '
529 WRITE(LO,*) ' ---- PHOJET version 1.12 ---- '
530 WRITE(LO,*) ' '
531 WRITE(LO,*) ' ==================================================='
532 WRITE(LO,*) ' Authors: Ralph Engel (FZ Karlsruhe)'
533 WRITE(LO,*) ' Johannes Ranft (Siegen Univ.)'
534 WRITE(LO,*) ' Stefan Roesler (CERN)'
535 WRITE(LO,*) ' ---------------------------------------------------'
536 WRITE(LO,*) ' Manual, updates, and further information:'
537 WRITE(LO,*) ' http://www-ik.fzk.de/~engel/phojet.html'
538 WRITE(LO,*) ' ---------------------------------------------------'
539 WRITE(LO,*) ' please send suggestions / bug reports etc. to:'
540 WRITE(LO,*) ' ralph.engel@fzk.de'
541 WRITE(LO,*) ' ==================================================='
542 WRITE(LO,*) ' $Date: 2000/06/25 21:59:19 $'
543 WRITE(LO,*) ' $Revision: 1.12.1.35 $'
544
545 WRITE(LO,*) ' (code version with interface to PYTHIA 6.x)'
546
547 WRITE(LO,*) ' (code version for usage in DPMJET 3.x)'
548
549 WRITE(LO,*) ' ==================================================='
550 WRITE(LO,*)
551
552C standard initializations
553 CALL PHO_DATINI
554 CALL PHO_PARDAT
555 DUM = PHO_PMASS(0,-1)
556
557C initialize standard PDFs
558C proton
559 CALL PHO_SETPDF(2212,IDUM,5,6,0,0,-1)
560 CALL PHO_SETPDF(-2212,IDUM,5,6,0,0,-1)
561C neutron
562 CALL PHO_SETPDF(2112,IDUM,5,6,0,0,-1)
563 CALL PHO_SETPDF(-2112,IDUM,5,6,0,0,-1)
564C photon
565 CALL PHO_SETPDF(22,IDUM,5,3,0,0,-1)
566C pomeron
567 CALL PHO_SETPDF(990,IDUM,4,0,0,0,-1)
568C pions
569 CALL PHO_SETPDF(211,IDUM,5,2,0,0,-1)
570 CALL PHO_SETPDF(-211,IDUM,5,2,0,0,-1)
571 CALL PHO_SETPDF(111,IDUM,5,2,0,0,-1)
572C kaons
573 CALL PHO_SETPDF(321,IDUM,5,2,0,0,-1)
574 CALL PHO_SETPDF(-321,IDUM,5,2,0,0,-1)
575 CALL PHO_SETPDF(130,IDUM,5,2,0,0,-1)
576 CALL PHO_SETPDF(310,IDUM,5,2,0,0,-1)
577
578C nothing to be done
579 IF(LINP.LT.0) RETURN
580
581C main loop to read input cards
582 1200 CONTINUE
583 READ(LINP,14,END=1300) CNAME,NUMBER
584 IF(CNAME.EQ.'ENDINPUT ') THEN
585 GOTO 1300
586 ELSE IF(CNAME.EQ.'STOP ') THEN
587 WRITE(LO,*) 'STOP'
588 STOP
589 ELSE IF(CNAME.EQ.'COMMENT ') THEN
590 WRITE(LO,'(1X,A10,A69)') 'COMMENT ',NUMBER
591 ELSE IF(CNAME(1:1).EQ.'*') THEN
592 WRITE(LO,'(1X,A10,A69)') CNAME,NUMBER
593 ELSE IF(CNAME.EQ.'PTCUT ') THEN
594 READ(NUMBER,*) PARMDL(36),PARMDL(37),PARMDL(38),PARMDL(39)
595 WRITE(LO,*) 'PTCUT ',PARMDL(36),PARMDL(37),
596 & PARMDL(38),PARMDL(39)
597 ELSE IF(CNAME.EQ.'PROCESS ') THEN
598 READ(NUMBER,*) (IPRON(KK,1),KK=1,8)
599 WRITE(LO,*) 'PROCESS ',(IPRON(KK,1),KK=1,8)
600 ELSE IF(CNAME.EQ.'DIFF-PROC ') THEN
601 READ(NUMBER,*) (ITMP(KK),KK=0,11)
602 WRITE(LO,*) 'DIFF-PROC ',(ITMP(KK),KK=0,8)
603 DO 112 KK=1,8
604 IPRON(KK,ITMP(0)) = ITMP(KK)
605 112 CONTINUE
606 ELSE IF(CNAME.EQ.'SUBPROCESS') THEN
607 READ(NUMBER,*) IMPRO,IP,ION
608 WRITE(LO,*) 'SUBPROCESS',IMPRO,IP,ION
609 MH_pro_on(IMPRO,IP) = ION
610 ELSE IF(CNAME.EQ.'PARTICLE1 ') THEN
611 READ(NUMBER,*) IDPDG,PVIR
612 IHFLS(1) = 1
613 XPSUB = 1.D0
614 CALL PHO_SETPAR(1,IDPDG,0,PVIR)
615 WRITE(LO,*) 'PARTICLE1 ',IDPDG,PVIR
616 ELSE IF(CNAME.EQ.'PARTICLE2 ') THEN
617 READ(NUMBER,*) IDPDG,PVIR
618 IHFLS(2) = 1
619 XTSUB = 1.D0
620 CALL PHO_SETPAR(2,IDPDG,0,PVIR)
621 WRITE(LO,*) 'PARTICLE2 ',IDPDG,PVIR
622 ELSE IF(CNAME.EQ.'REMNANT1 ') THEN
623 READ(NUMBER,*) IDPDG,IFL1,IFL2,IVAL,XSUB
624 IHFLS(1) = IVAL
625 IHFLD(1,1) = IFL1
626 IHFLD(1,2) = IFL2
627 XPSUB = XSUB
628 PVIR = 0.D0
629 CALL PHO_SETPAR(1,IDPDG,-1,PVIR)
630 WRITE(LO,*) 'REMNANT1 ',IDPDG,IFL1,IFL2,IVAL,XSUB
631 ELSE IF(CNAME.EQ.'REMNANT2 ') THEN
632 READ(NUMBER,*) IDPDG,IFL1,IFL2,IVAL,XSUB
633 IHFLS(2) = IVAL
634 IHFLD(2,1) = IFL1
635 IHFLD(2,2) = IFL2
636 XTSUB = XSUB
637 PVIR = 0.D0
638 CALL PHO_SETPAR(2,IDPDG,-1,PVIR)
639 WRITE(LO,*) 'REMNANT2 ',IDPDG,IFL1,IFL2,IVAL,XSUB
640 ELSE IF(CNAME.EQ.'PDF ') THEN
641 READ(NUMBER,*) IDPDG,IPAR,ISET,IEXT
642 WRITE(LO,*) 'PDF ',IDPDG,IPAR,ISET,IEXT
643 CALL PHO_SETPDF(IDPDG,IDUM,IPAR,ISET,IEXT,0,-1)
644 ELSE IF(CNAME.EQ.'SETMODEL ') THEN
645 READ(NUMBER,*) I,IVAL
646 WRITE(LO,*) 'SETMODEL ',I,IVAL
647 CALL PHO_SETMDL(I,IVAL,1)
648 ELSE IF(CNAME.EQ.'SETPARAM ') THEN
649 READ(NUMBER,*) I,PARNEW
650 WRITE(LO,*) 'SETPARAM ',I,PARNEW
651 PARMDL(I) = PARNEW
652 ELSE IF(CNAME.EQ.'DEBUG ') THEN
653 READ(NUMBER,*) IDEBF,IDEBN,IDLEV
654 WRITE(LO,*) 'DEBUG ',IDEBF,IDEBN,IDLEV
655 CALL PHO_TRACE(IDEBF,IDEBN,IDLEV)
656 ELSE IF(CNAME.EQ.'TRACE ') THEN
657 READ(NUMBER,*) IDEBF,IDLEV
658 WRITE(LO,*) 'TRACE ',IDEBF,IDLEV
659 IDEB(IDEBF) = IDLEV
660 ELSE IF(CNAME.EQ.'SETICUT ') THEN
661 READ(NUMBER,*) I,ICUT
662 WRITE(LO,*) 'SETICUT ',I,ICUT
663 ISWCUT(I) = ICUT
664 ELSE IF(CNAME.EQ.'SETFCUT ') THEN
665 READ(NUMBER,*) I,PARNEW
666 WRITE(LO,*) 'SETFCUT ',I,PARNEW
667 HSWCUT(I) = PARNEW
668 ELSE IF(CNAME.EQ.'LUND-MSTU ') THEN
669 READ(NUMBER,*) I,IVAL
670 WRITE(LO,*) 'LUND-MSTU ',I,IVAL
671 MSTU(I) = IVAL
672 ELSE IF(CNAME.EQ.'LUND-MSTJ ') THEN
673 READ(NUMBER,*) I,IVAL
674 WRITE(LO,*) 'LUND-MSTJ ',I,IVAL
675 MSTJ(I) = IVAL
676 ELSE IF(CNAME.EQ.'LUND-PARJ ') THEN
677 READ(NUMBER,*) I,EE
678 WRITE(LO,*) 'LUND-PARJ ',I,EE
679 PARJ(I) = REAL(EE)
680 ELSE IF(CNAME.EQ.'LUND-PARU ') THEN
681 READ(NUMBER,*) I,EE
682 WRITE(LO,*) 'LUND-PARU ',I,EE
683 PARU(I) = REAL(EE)
684 ELSE IF(CNAME.EQ.'LUND-DECAY') THEN
685 READ(NUMBER,*) ID,ION
686 WRITE(LO,*) 'LUND-DECAY ',ID,ION
687
688 KC=PYCOMP(ID)
689
690 MDCY(KC,1) = ION
691 ELSE IF(CNAME.EQ.'PSOFTMIN ') THEN
692 READ(NUMBER,*) PSOMIN
693 WRITE(LO,*) 'PSOFTMIN ',PSOMIN
694 ELSE IF(CNAME.EQ.'INTPREC ') THEN
695 READ(NUMBER,*) NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
696 WRITE(LO,*) 'INTPREC ',NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
697
698C PDF test utility
699 ELSE IF(CNAME.EQ.'PDFTEST ') THEN
700 READ(NUMBER,*) IDPDG,SCALE2,PVIRT2
701 PVIRT2 = ABS(PVIRT2)
702 WRITE(LO,*) 'PDFTEST ',IDPDG,' ',SCALE2,' ',PVIRT2
703 CALL PHO_PDFTST(IDPDG,SCALE2,PVIRT2)
704
705C mass cut on gamma-gamma or gamma-hadron system
706 ELSE IF(CNAME.EQ.'ECMS-CUT ') THEN
707 READ(NUMBER,*) ECMIN,ECMAX
708 WRITE(LO,*) 'ECMS-CUT ',ECMIN,ECMAX
709
710C beam lepton (anti-)tagging system
711 ELSE IF(CNAME.EQ.'TAG-METHOD') THEN
712 READ(NUMBER,*) ITAG1,ITAG2
713 WRITE(LO,*) 'TAG-METHOD',ITAG1,ITAG2
714 ELSE IF(CNAME.EQ.'E-TAG1 ') THEN
715 READ(NUMBER,*)
716 & EEMIN1,YMIN1,YMAX1,Q2MIN1,Q2MAX1,THMIN1,THMAX1
717 WRITE(LO,*) 'E-TAG1 ',EEMIN1,YMIN1,YMAX1,
718 & Q2MIN1,Q2MAX1,THMIN1,THMAX1
719 ELSE IF(CNAME.EQ.'E-TAG2 ') THEN
720 READ(NUMBER,*)
721 & EEMIN2,YMIN2,YMAX2,Q2MIN2,Q2MAX2,THMIN2,THMAX2
722 WRITE(LO,*) 'E-TAG2 ',EEMIN2,YMIN2,YMAX2,
723 & Q2MIN2,Q2MAX2,THMIN2,THMAX2
724
725C sampling of gamma-p events in ep (HERA)
726 ELSE IF( (CNAME.EQ.'WW-HERA ')
727 & .OR.(CNAME.EQ.'GP-HERA ')) THEN
728 READ(NUMBER,*) EE1,EE2,NEV
729 WRITE(LO,*) 'GP-HERA ',EE1,EE2,NEV
730 IF(YMAX2.LT.0.D0) THEN
731 WRITE(LO,*) ' PHO_INIT:ERROR:ELECTRON TAGGER NOT SET'
732 ELSE
733 CALL PHO_GPHERA(NEV,EE1,EE2)
734 KEVENT = 0
735 ENDIF
736
737C sampling of gamma-gamma events in e+e- (LEP)
738 ELSE IF( (CNAME.EQ.'GG-EPEM ')
739 & .OR.(CNAME.EQ.'WW-EPEM ')) THEN
740 READ(NUMBER,*) EE1,EE2,NEV
741 WRITE(LO,*) 'GG-EPEM ',EE1,EE2,NEV
742 IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
743 WRITE(LO,*) ' PHO_INIT:ERROR:ELECTRON TAGGERS NOT SET'
744 ELSE
745 CALL PHO_GGEPEM(-1,EE1,EE2)
746 CALL PHO_GGEPEM(NEV,EE1,EE2)
747 CALL PHO_GGEPEM(-2,sig_tot,sig_gg)
748 KEVENT = 0
749 ENDIF
750
751C sampling of gamma-gamma in heavy-ion collisions
752 ELSE IF(CNAME.EQ.'GG-HION-F ') THEN
753 READ(NUMBER,*) EE,NA,NZ,NEV
754 WRITE(LO,*) 'GG-HION-F ',EE,NA,NZ,NEV
755 IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
756 WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
757 ELSE
758 CALL PHO_GGHIOF(NEV,EE,NA,NZ)
759 KEVENT = 0
760 ENDIF
761 ELSE IF(CNAME.EQ.'GG-HION-G ') THEN
762 READ(NUMBER,*) EE,NA,NZ,NEV
763 WRITE(LO,*) 'GG-HION-G ',EE,NA,NZ,NEV
764 IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
765 WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
766 ELSE
767 CALL PHO_GGHIOG(NEV,EE,NA,NZ)
768 KEVENT = 0
769 ENDIF
770
771C sampling of gamma-hadron events in heavy ion collisions
772 ELSE IF(CNAME.EQ.'GH-HION-F ') THEN
773 READ(NUMBER,*) EE,NA,NZ,NEV
774 WRITE(LO,*) 'GH-HION-F ',EE,NA,NZ,NEV
775 IF((YMAX1.LT.0.D0).OR.(YMAX2.LT.0.D0)) THEN
776 WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
777 ELSE
778 CALL PHO_GHHIOF(NEV,EE,NA,NZ)
779 KEVENT = 0
780 ENDIF
781
782C sampling of hadron-gamma events in hadron - heavy ion collisions
783 ELSE IF(CNAME.EQ.'HG-HIAS-F ') THEN
784 READ(NUMBER,*) EP,EE,NA,NZ,NEV
785 WRITE(LO,*) 'HG-HIAS-F ',EP,EE,NA,NZ,NEV
786 IF(YMAX2.LT.0.D0) THEN
787 WRITE(LO,*) ' PHO_INIT:ERROR:Y RANGE FOR PHOTONS NOT SET'
788 ELSE
789 CALL PHO_GHHIAS(NEV,EP,EE,NA,NZ)
790 KEVENT = 0
791 ENDIF
792
793C sampling of photoproduction events e+e-, backscattered laser
794 ELSE IF(CNAME.EQ.'BLASER ') THEN
795 READ(NUMBER,*) EE1,EE2,Pl_lam_1,Pl_lam_2,X_1,X_2,rho,A,NEV
796 WRITE(LO,*) 'BLASER ',EE1,EE2,
797 & Pl_lam_1,Pl_lam_2,X_1,X_2,rho,A,NEV
798 CALL PHO_GGBLSR(NEV,EE1,EE2,Pl_lam_1,Pl_lam_2,X_1,X_2,rho,A)
799 KEVENT = 0
800
801C sampling of photoproduction events beamstrahlung
802 ELSE IF(CNAME.EQ.'BEAMST ') THEN
803 READ(NUMBER,*) EE1,YPSI,SIGX,SIGY,SIGZ,AEB,NEV
804 WRITE(LO,*) 'BEAMST ',EE1,YPSI,SIGX,SIGY,SIGZ,AEB,NEV
805 IF(YMAX1.LT.0.D0) THEN
806 WRITE(LO,*) ' PHO_INIT:ERROR:ELECTRON TAGGER 1 NOT SET'
807 ELSE
808 CALL PHO_GGBEAM(NEV,EE1,YPSI,SIGX,SIGY,SIGZ,AEB)
809 KEVENT = 0
810 ENDIF
811
812C fixed-energy events in LAB system of particle 2
813 ELSE IF(CNAME.EQ.'EVENT-LAB ') THEN
814 READ(NUMBER,*) PLAB,NEV
815 WRITE(LO,*) 'EVENT-LAB ',PLAB,NEV
816 CALL PHO_FIXLAB(PLAB,NEV)
817 KEVENT = 0
818
819C fixed-energy events in CM system
820 ELSE IF(CNAME.EQ.'EVENT-CMS ') THEN
821 READ(NUMBER,*) ECM,NEV
822 WRITE(LO,*) 'EVENT-CMS ',ECM,NEV
823 PMASS1 = PHO_PMASS(IFPAB(1),0)-SQRT(PVIRT(1))
824 PMASS2 = PHO_PMASS(IFPAB(2),0)-SQRT(PVIRT(2))
825 CALL PHO_PECMS(1,PMASS1,PMASS2,ECM,PCM,EE)
826 E1 = EE
827 E2 = ECM-EE
828 THETA = 0.D0
829 PHI = 0.D0
830 CALL PHO_FIXCOL(E1,E2,THETA,PHI,NEV)
831 KEVENT = 0
832
833C fixed-energy events for collider setup with crossing angle
834 ELSE IF(CNAME.EQ.'EVENT-COLL') THEN
835 READ(NUMBER,*) E1,E2,THETA,PHI,NEV
836 WRITE(LO,*) 'EVENT-COLL',E1,E2,THETA,PHI,NEV
837 CALL PHO_FIXCOL(E1,E2,THETA,PHI,NEV)
838 KEVENT = 0
839
840C unknown data card
841 ELSE
842 WRITE(LO,*) 'PHO_INIT: unknown data card: ',CNAME,NUMBER
843 ENDIF
844
845 GOTO 1200
846 1300 CONTINUE
847 WRITE(LO,*) ' RETURN'
848
849 END
850
851*$ CREATE PHO_SETMDL.FOR
852*COPY PHO_SETMDL
853CDECK ID>, PHO_SETMDL
854 SUBROUTINE PHO_SETMDL(INDX,IVAL,IMODE)
855C**********************************************************************
856C
857C set model switches
858C
859C input: INDX model parameter number
860C (positive: ISWMDL, negative: IPAMDL)
861C IVAL new value
862C IMODE -1 print value of parameter INDX
863C 1 set new value
864C -2 print current settings
865C
866C**********************************************************************
867 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
868 SAVE
869
870C input/output channels
871 INTEGER LI,LO
872 COMMON /POINOU/ LI,LO
873C model switches and parameters
874 CHARACTER*8 MDLNA
875 INTEGER ISWMDL,IPAMDL
876 DOUBLE PRECISION PARMDL
877 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
878
879 IF(IMODE.EQ.-2) THEN
880 WRITE(LO,'(/1X,A,/1X,A,/)') 'PHO_SETMDL: current settings',
881 & '----------------------------'
882 DO 100 I=1,48,3
883 IF(ISWMDL(I).EQ.-9999) GOTO 200
884 IF(ISWMDL(I+1).EQ.-9999) THEN
885 WRITE(LO,'(5X,I3,A1,A,I6)') I,':',MDLNA(I),ISWMDL(I)
886 GOTO 200
887 ELSE IF(ISWMDL(I+2).EQ.-9999) THEN
888 WRITE(LO,'(2(5X,I3,A1,A,I6))') I,':',MDLNA(I),ISWMDL(I),
889 & I+1,':',MDLNA(I+1),ISWMDL(I+1)
890 GOTO 200
891 ELSE
892 WRITE(LO,'(3(5X,I3,A1,A,I6))')
893 & (I+K,':',MDLNA(I+K),ISWMDL(I+K),K=0,2)
894 ENDIF
895 100 CONTINUE
896 200 CONTINUE
897 ELSE IF(IMODE.EQ.-1) THEN
898 WRITE(LO,'(1X,A,1X,A,I6)')
899 & 'PHO_SETMDL:',MDLNA(INDX),ISWMDL(INDX)
900 ELSE IF(IMODE.EQ.1) THEN
901 IF(INDX.GT.0) THEN
902 IF(ISWMDL(INDX).NE.IVAL) THEN
903 WRITE(LO,'(1X,A,I4,1X,A,2I6)')
904 & 'PHO_SETMDL:ISWMDL(OLD/NEW):',
905 & INDX,MDLNA(INDX),ISWMDL(INDX),IVAL
906 ISWMDL(INDX) = IVAL
907 ENDIF
908 ELSE IF(INDX.LT.0) THEN
909 IF(IPAMDL(-INDX).NE.IVAL) THEN
910 WRITE(LO,'(1X,A,I4,1X,2I6)') 'PHO_SETMDL:IPAMDL(OLD/NEW):',
911 & -INDX,IPAMDL(-INDX),IVAL
912 IPAMDL(-INDX) = IVAL
913 ENDIF
914 ENDIF
915 ELSE
916 WRITE(LO,'(/1X,A,I6)')
917 & 'PHO_SETMDL:ERROR: unsupported mode',IMODE
918 ENDIF
919 END
920
921*$ CREATE PHO_DATINI.FOR
922*COPY PHO_DATINI
923CDECK ID>, PHO_DATINI
924 SUBROUTINE PHO_DATINI
925C*********************************************************************
926C
927C initialization of variables and switches
928C
929C*********************************************************************
930 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
931 SAVE
932
933C input/output channels
934 INTEGER LI,LO
935 COMMON /POINOU/ LI,LO
936C some constants
937 DOUBLE PRECISION PI,PI2,PI4,GEV2MB,Q_ch,Q_ch2,Q_ch4
938 COMMON /POCONS/ PI,PI2,PI4,GEV2MB,
939 & Q_ch(-6:6),Q_ch2(-6:6),Q_ch4(-6:6)
940C event debugging information
941 INTEGER NMAXD
942 PARAMETER (NMAXD=100)
943 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
944 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
945 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
946 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
947C event weights and generated cross section
948 INTEGER IPOWGC,ISWCUT,IVWGHT
949 DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
950 COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
951 & IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
952C scale parameters for parton model calculations
953 INTEGER NQQAL,NQQALI,NQQALF,NQQPD
954 DOUBLE PRECISION AQQAL,AQQALI,AQQALF,AQQPD
955 COMMON /POHSCL/ AQQAL,AQQALI,AQQALF,AQQPD,
956 & NQQAL,NQQALI,NQQALF,NQQPD
957C integration precision for hard cross sections (obsolete)
958 INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
959 COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
960C hard scattering parameters used for most recent hard interaction
961 INTEGER NFbeta,NF
962 DOUBLE PRECISION ALQCD2,BQCD
963 COMMON /POHAPA/ ALQCD2(3,4),BQCD(4),NFbeta,NF
964C cut probability distribution
965 INTEGER IEETA1,IIMAX,KKMAX
966 PARAMETER( IEETA1=20, IIMAX=20, KKMAX=20 )
967 INTEGER IEEMAX,IMAX,KMAX
968 REAL PROB
969 DOUBLE PRECISION EPTAB
970 COMMON /POPROB/ PROB(4,IEETA1,0:IIMAX,0:KKMAX),EPTAB(4,IEETA1),
971 & IEEMAX,IMAX,KMAX
972C gamma-lepton or gamma-hadron vertex information
973 INTEGER IGHEL,IDPSRC,IDBSRC
974 DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
975 & RADSRC,AMSRC,GAMSRC
976 COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
977 & GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
978 & IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
979C photon flux kinematics and cuts
980 DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
981 & YMIN1,YMAX1,YMIN2,YMAX2,
982 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
983 & THMIN1,THMAX1,THMIN2,THMAX2
984 INTEGER ITAG1,ITAG2
985 COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
986 & YMIN1,YMAX1,YMIN2,YMAX2,
987 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
988 & THMIN1,THMAX1,THMIN2,THMAX2,
989 & ITAG1,ITAG2
990C obsolete cut-off information
991 DOUBLE PRECISION PTCUT,PTANO,FPS,FPH,PSOMIN,XSOMIN
992 COMMON /POCUT1/ PTCUT(4),PTANO(4),FPS(4),FPH(4),PSOMIN,XSOMIN
993C global event kinematics and particle IDs
994 INTEGER IFPAP,IFPAB
995 DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
996 COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
997C nucleon-nucleus / nucleus-nucleus interface to DPMJET
998 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
999 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
1000 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
1001 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
1002C some hadron information, will be deleted in future versions
1003 INTEGER NFS
1004 DOUBLE PRECISION QMASS,BET,PCOUDI,PNORM,VALPRG
1005 COMMON /POHDRN/ QMASS(6),BET,PCOUDI,PNORM,VALPRG(2),NFS
1006C model switches and parameters
1007 CHARACTER*8 MDLNA
1008 INTEGER ISWMDL,IPAMDL
1009 DOUBLE PRECISION PARMDL
1010 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
1011C general process information
1012 INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
1013 COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
1014C parameters of the "simple" Vector Dominance Model
1015 DOUBLE PRECISION VMAS,GAMM,RMIN,RMAX,VMSL,VMFA
1016 COMMON /POSVDM/ VMAS(4),GAMM(4),RMIN(4),RMAX(4),VMSL(4),VMFA(4)
1017C parameters for DGLAP backward evolution in ISR
1018 INTEGER NFSISR
1019 DOUBLE PRECISION Q2MISR,PMISR,ZMISR,AL2ISR
1020 COMMON /PODGL1/ Q2MISR(2),PMISR(2),ZMISR(2),AL2ISR(2),NFSISR
1021C particles created by initial state evolution
1022 INTEGER MXISR1,MXISR2
1023 PARAMETER ( MXISR1 = 150, MXISR2 = 50 )
1024 INTEGER IFLISR,IPOISR,IMXISR
1025 DOUBLE PRECISION PHISR
1026 COMMON /POPISR/ IFLISR(2,MXISR1),PHISR(2,4,MXISR1),
1027 & IPOISR(2,2,MXISR2),IMXISR(2)
1028C names of hard scattering processes
1029 INTEGER Max_pro_1
1030 PARAMETER ( Max_pro_1 = 16 )
1031 CHARACTER*18 PROC
1032 COMMON /POHPRO/ PROC(0:Max_pro_1)
1033C hard cross sections and MC selection weights
1034 INTEGER Max_pro_2
1035 PARAMETER ( Max_pro_2 = 16 )
1036 INTEGER IHa_last,IHb_last,MH_pro_on,MH_tried,
1037 & MH_acc_1,MH_acc_2
1038 DOUBLE PRECISION Hfac,HWgx,HSig,Hdpt,HEcm_last,HQ2a_last,HQ2b_last
1039 COMMON /POHRCS/ Hfac(-1:Max_pro_2),HWgx(-1:Max_pro_2),
1040 & HSig(-1:Max_pro_2),Hdpt(-1:Max_pro_2),
1041 & HEcm_last,HQ2a_last,HQ2b_last,IHa_last,IHb_last,
1042 & MH_pro_on(-1:Max_pro_2,0:4),MH_tried(-1:Max_pro_2,0:4),
1043 & MH_acc_1(-1:Max_pro_2,0:4),MH_acc_2(-1:Max_pro_2,0:4)
1044C interpolation tables for hard cross section and MC selection weights
1045 INTEGER Max_tab_E,Max_tab_Q2,Max_pro_tab
1046 PARAMETER ( Max_tab_E = 20, Max_tab_Q2 = 10, Max_pro_tab = 16 )
1047 INTEGER IH_Q2a_up,IH_Q2b_up,IH_Ecm_up
1048 DOUBLE PRECISION Hfac_tab,HWgx_tab,HSig_tab,Hdpt_tab,
1049 & HQ2a_tab,HQ2b_tab,HEcm_tab
1050 COMMON /POHTAB/
1051 & Hfac_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1052 & HWgx_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1053 & HSig_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1054 & Hdpt_tab(-1:Max_pro_tab,Max_tab_E,Max_tab_Q2,Max_tab_Q2,0:4),
1055 & HQ2a_tab(1:Max_tab_Q2,0:4),HQ2b_tab(1:Max_tab_Q2,0:4),
1056 & HEcm_tab(1:Max_tab_E,0:4),
1057 & IH_Q2a_up(0:4),IH_Q2b_up(0:4),IH_Ecm_up(0:4)
1058
1059C initialize /POCONS/
1060 PI = ATAN(1.D0)*4.D0
1061 PI2 = 2.D0*PI
1062 PI4 = 2.D0*PI2
1063C GeV**-2 --> millibarn (multiply by GEV2MB to get mb as units)
1064 GEV2MB = 0.389365D0
1065C precalculate quark charges
1066 do i=1,6
1067 Q_ch(i) = dble(2-3*mod(i,2))/3.D0
1068 Q_ch(-i) = -Q_ch(i)
1069
1070 Q_ch2(i) = Q_ch(i)**2
1071 Q_ch2(-i) = Q_ch2(i)
1072
1073 Q_ch4(i) = Q_ch2(i)**2
1074 Q_ch4(-i) = Q_ch4(i)
1075 enddo
1076 Q_ch(0) = 0.D0
1077 Q_ch2(0) = 0.D0
1078 Q_ch4(0) = 0.D0
1079
1080C initialize /GLOCMS/
1081 ECM = 50.D0
1082 PMASS(1) = 0.D0
1083 PVIRT(1) = 0.D0
1084 PMASS(2) = 0.D0
1085 PVIRT(2) = 0.D0
1086 IFPAP(1) = 22
1087 IFPAP(2) = 22
1088C initialize /HADVAL/
1089 IHFLD(1,1) = 0
1090 IHFLD(1,2) = 0
1091 IHFLD(2,1) = 0
1092 IHFLD(2,2) = 0
1093 IHFLS(1) = 1
1094 IHFLS(2) = 1
1095C initialize /MODELS/
1096 ISWMDL(1) = 3
1097 MDLNA(1) = 'AMPL MOD'
1098 ISWMDL(2) = 1
1099 MDLNA(2) = 'MIN-BIAS'
1100 ISWMDL(3) = 1
1101 MDLNA(3) = 'PTS DISH'
1102 ISWMDL(4) = 1
1103 MDLNA(4) = 'PTS DISP'
1104 ISWMDL(5) = 2
1105 MDLNA(5) = 'PTS ASSI'
1106 ISWMDL(6) = 3
1107 MDLNA(6) = 'HADRONIZ'
1108 ISWMDL(7) = 2
1109 MDLNA(7) = 'MASS COR'
1110 ISWMDL(8) = 3
1111 MDLNA(8) = 'PAR SHOW'
1112 ISWMDL(9) = 0
1113 MDLNA(9) = 'GLU SPLI'
1114 ISWMDL(10) = 2
1115 MDLNA(10) = 'VIRT PHO'
1116 ISWMDL(11) = 0
1117 MDLNA(11) = 'LARGE NC'
1118 ISWMDL(12) = 0
1119 MDLNA(12) = 'LIPA POM'
1120 ISWMDL(13) = 1
1121 MDLNA(13) = 'QELAS VM'
1122 ISWMDL(14) = 2
1123 MDLNA(14) = 'ENHA GRA'
1124 ISWMDL(15) = 4
1125 MDLNA(15) = 'MULT SCA'
1126 ISWMDL(16) = 4
1127 MDLNA(16) = 'MULT DIF'
1128 ISWMDL(17) = 4
1129 MDLNA(17) = 'MULT CDF'
1130 ISWMDL(18) = 0
1131 MDLNA(18) = 'BALAN PT'
1132 ISWMDL(19) = 1
1133 MDLNA(19) = 'POMV FLA'
1134 ISWMDL(20) = 0
1135 MDLNA(20) = 'SEA FLA'
1136 ISWMDL(21) = 2
1137 MDLNA(21) = 'SPIN DEC'
1138 ISWMDL(22) = 1
1139 MDLNA(22) = 'DIF.MASS'
1140 ISWMDL(23) = 1
1141 MDLNA(23) = 'DIFF RES'
1142 ISWMDL(24) = 0
1143 MDLNA(24) = 'PTS HPOM'
1144 ISWMDL(25) = 0
1145 MDLNA(25) = 'POM CORR'
1146 ISWMDL(26) = 1
1147 MDLNA(26) = 'OVERLAP '
1148 ISWMDL(27) = 0
1149 MDLNA(27) = 'MUL R/AN'
1150 ISWMDL(28) = 1
1151 MDLNA(28) = 'SUR PROB'
1152 ISWMDL(29) = 1
1153 MDLNA(29) = 'PRIMO KT'
1154 ISWMDL(30) = 0
1155 MDLNA(30) = 'DIFF. CS'
1156 ISWMDL(31) = -9999
1157C mass-independent sea flavour ratios (for low-mass strings)
1158 PARMDL(1) = 0.425D0
1159 PARMDL(2) = 0.425D0
1160 PARMDL(3) = 0.15D0
1161 PARMDL(4) = 0.D0
1162 PARMDL(5) = 0.D0
1163 PARMDL(6) = 0.D0
1164C suppression by energy momentum conservation
1165 PARMDL(8) = 9.D0
1166 PARMDL(9) = 7.D0
1167C VDM factors
1168 PARMDL(10) = 0.866D0
1169 PARMDL(11) = 0.288D0
1170 PARMDL(12) = 0.288D0
1171 PARMDL(13) = 0.288D0
1172 PARMDL(14) = 0.866D0
1173 PARMDL(15) = 0.288D0
1174 PARMDL(16) = 0.288D0
1175 PARMDL(17) = 0.288D0
1176 PARMDL(18) = 0.D0
1177C lower energy limit for initialization
1178 PARMDL(19) = 5.D0
1179C soft pt for hard scattering remnants
1180 PARMDL(20) = 5.D0
1181C low energy beta of soft pt distribution 1
1182 PARMDL(21) = 4.5D0
1183C high energy beta of soft pt distribution 1
1184 PARMDL(22) = 3.0D0
1185C low energy beta of soft pt distribution 0
1186 PARMDL(23) = 2.5D0
1187C high energy beta of soft pt distribution 0
1188 PARMDL(24) = 0.4D0
1189C effective quark mass in photon wave function
1190 PARMDL(25) = 0.2D0
1191C normalization of unevolved Pomeron PDFs
1192 PARMDL(26) = 0.3D0
1193C effective VDM parameters for Q**2 dependence of cross section
1194 PARMDL(27) = 0.65D0
1195 PARMDL(28) = 0.08D0
1196 PARMDL(29) = 0.05D0
1197 PARMDL(30) = 0.22D0
1198 PARMDL(31) = 0.589824D0
1199 PARMDL(32) = 0.609961D0
1200 PARMDL(33) = 1.038361D0
1201 PARMDL(34) = 1.96D0
1202C Q**2 suppression of multiple interactions
1203 PARMDL(35) = 0.59D0
1204C pt cutoff defaults
1205 PARMDL(36) = 2.5D0
1206 PARMDL(37) = 2.5D0
1207 PARMDL(38) = 2.5D0
1208 PARMDL(39) = 2.5D0
1209C enhancement factor for diffractive cross sections
1210 PARMDL(40) = 1.D0
1211 PARMDL(41) = 1.D0
1212 PARMDL(42) = 1.D0
1213C mass in soft pt distribution
1214 PARMDL(43) = 0.D0
1215C maximum of x allowed for leading particle
1216 PARMDL(44) = 0.9D0
1217C max. mass sampled in diffraction
1218 PARMDL(45) = sqrt(0.4D0)
1219C mass threshold in diffraction (2pi mass)
1220 PARMDL(46) = 0.3D0
1221C regularization of slope parameter in diffraction
1222 PARMDL(47) = 4.D0
1223C renormalized intercept for enhanced graphs
1224 PARMDL(48) = 1.08D0
1225C coherence constraint for diff. cross sections
1226 PARMDL(49) = sqrt(0.05D0)
1227C exponents of x distributions
1228C baryon
1229 PARMDL(50) = 1.5D0
1230 PARMDL(51) = -0.5D0
1231 PARMDL(52) = -0.99D0
1232 PARMDL(53) = -0.99D0
1233C meson (non-strangeness part)
1234 PARMDL(54) = -0.5D0
1235 PARMDL(55) = -0.5D0
1236 PARMDL(56) = -0.99D0
1237 PARMDL(57) = -0.99D0
1238C meson (strangeness part)
1239 PARMDL(58) = -0.2D0
1240 PARMDL(59) = -0.2D0
1241 PARMDL(60) = -0.99D0
1242 PARMDL(61) = -0.99D0
1243C particle remnant (no valence quarks)
1244 PARMDL(62) = -0.5D0
1245 PARMDL(63) = -0.5D0
1246 PARMDL(64) = -0.99D0
1247 PARMDL(65) = -0.99D0
1248C ratio beetween triple-pomeron/reggeon couplings grrp/gppp
1249 PARMDL(66) = 10.D0
1250C ratio beetween triple-pomeron/reggeon couplings gppr/gppp
1251 PARMDL(67) = 10.D0
1252C min. abs(t) in diffraction
1253 PARMDL(68) = 0.D0
1254C max. abs(t) in diffraction
1255 PARMDL(69) = 10.D0
1256C min. mass for elastic pomerons in central diffraction
1257 PARMDL(70) = 2.D0
1258C min. mass of diffractive blob in central diffraction
1259 PARMDL(71) = 2.D0
1260C min. Feynman x cut in central diffraction
1261 PARMDL(72) = 0.D0
1262C direct pomeron coupling
1263 PARMDL(74) = 0.D0
1264C relative deviation allowed for energy-momentum conservation
1265C energy-momentum relative deviation
1266 PARMDL(75) = 0.01D0
1267C transverse momentum deviation
1268 PARMDL(76) = 0.01D0
1269C couplings for unitarization in diffraction
1270C non-unitarized pomeron coupling (sqrt(mb))
1271 PARMDL(77) = 3.D0
1272C rescaling factor for pomeron PDF
1273 PARMDL(78) = 3.D0
1274C coupling probabilities
1275 PARMDL(79) = 1.D0
1276 PARMDL(80) = 0.D0
1277C scales to calculate alpha-s of matrix element
1278 PARMDL(81) = 1.D0
1279 PARMDL(82) = 1.D0
1280 PARMDL(83) = 1.D0
1281C scales to calculate alpha-s of initial state radiation
1282 PARMDL(84) = 1.D0
1283 PARMDL(85) = 1.D0
1284 PARMDL(86) = 1.D0
1285C scales to calculate alpha-s of final state radiation
1286 PARMDL(87) = 1.D0
1287 PARMDL(88) = 1.D0
1288 PARMDL(89) = 1.D0
1289C scales to calculate PDFs
1290 PARMDL(90) = 1.D0
1291 PARMDL(91) = 1.D0
1292 PARMDL(92) = 1.D0
1293C scale for ISR starting virtuality
1294 PARMDL(93) = 1.D0
1295C min. virtuality to generate time-like showers in ISR
1296 PARMDL(94) = 2.D0
1297C factor to scale the max. allowed time-like parton shower virtuality
1298 PARMDL(95) = 4.D0
1299C max. transverse momentum for primordial kt
1300 PARMDL(100) = 2.D0
1301C weight factors for pt-distribution
1302 PARMDL(101) = 2.D0
1303 PARMDL(102) = 2.D0
1304 PARMDL(103) = 4.D0
1305 PARMDL(104) = 2.D0
1306 PARMDL(105) = 6.D0
1307 PARMDL(106) = 4.D0
1308C
1309* PARMDL(110-125) reserved for hard scattering
1310C currently chosen scales for hard scattering
1311 DO 10 I=1,16
1312 PARMDL(109+I) = 0.D0
1313 10 CONTINUE
1314C virtuality cutoff in initial state evolution
1315 PARMDL(126) = PARMDL(36)**2
1316 PARMDL(127) = PARMDL(37)**2
1317 PARMDL(128) = PARMDL(38)**2
1318 PARMDL(129) = PARMDL(39)**2
1319C virtuality cutoff for direct contribution to photon PDF
1320 PARMDL(130) = 1.D30
1321 PARMDL(131) = 1.D30
1322 PARMDL(132) = 1.D30
1323 PARMDL(133) = 1.D30
1324C fraction of events without popcorn
1325 PARMDL(134) = -1.D0
1326C fraction of diquarks with spin 1 (relative to sum of spin 1 and 0)
1327 PARMDL(135) = 0.5D0
1328C soft color re-connection (fraction)
1329C g g final state
1330 PARMDL(140) = 1.D0/64.D0
1331C g q final state
1332 PARMDL(141) = 1.D0/24.D0
1333C q q final state
1334 PARMDL(142) = 1.D0/9.D0
1335C effective scale in Drees-Godbole like suppresion in photon PDF
1336 PARMDL(144) = 0.766D0**2
1337C QCD scales (if PDF scales are not used, 4 active flavours)
1338 PARMDL(145) = 0.2D0**2
1339 PARMDL(146) = 0.2D0**2
1340 PARMDL(147) = 0.2D0**2
1341C threshold scales for variable flavour calculation (GeV**2)
1342 PARMDL(148) = 1.5D0**2
1343 PARMDL(149) = 4.5D0**2
1344 PARMDL(150) = 175.D0**2
1345C constituent quark masses
1346 PARMDL(151) = 0.3D0
1347 PARMDL(152) = 0.3D0
1348 PARMDL(153) = 0.5D0
1349 PARMDL(154) = 1.6D0
1350 PARMDL(155) = 5.D0
1351 PARMDL(156) = 174.D0
1352C min. masses of valence quark
1353 PARMDL(157) = 0.3D0
1354C min. masses of valence diquark
1355 PARMDL(158) = 0.8D0
1356C min. mass of sea quark
1357 PARMDL(159) = 0.D0
1358C suppression of strange quarks as photon valences
1359 PARMDL(160) = 0.2D0
1360C min. masses for strings (used in PHO_SOFTXX)
1361 PARMDL(161) = 1.D0
1362 PARMDL(162) = 1.D0
1363 PARMDL(163) = 1.D0
1364 PARMDL(164) = 1.D0
1365C min. momentum fraction for soft processes
1366 PARMDL(165) = 0.3D0
1367C min. phase space for x-sampling
1368 PARMDL(166) = 0.135D0
1369C Ross-Stodolsky exponent
1370 PARMDL(170) = 4.2D0
1371C cutoff on photon-pomeron invariant mass in hadron-hadron collisions
1372 PARMDL(175) = 2.D0
1373
1374**sr
1375* extra factor multiplying difference between Goulianos and PHOJET-
1376* diff. cross sections
1377 PARMDL(200) = 0.6D0
1378**
1379
1380C complex amplitudes, eikonal functions
1381 IPAMDL(1) = 0
1382C allow for Reggeon cuts
1383 IPAMDL(2) = 1
1384C decay of hadron resonances in diffraction (0 iso, 1 trans, 2 long)
1385 IPAMDL(3) = 0
1386C polarization of photon resonances (0 none, 1 trans, 2 long)
1387 IPAMDL(4) = 1
1388C pt of valence partons
1389 IPAMDL(5) = 1
1390C pt of hard scattering remnant
1391 IPAMDL(6) = 2
1392C running cutoff for hard scattering
1393 IPAMDL(7) = 1
1394C intercept used for the calculation of enhanced graphs
1395 IPAMDL(8) = 1
1396C effective slope of hard scattering amplitde
1397 IPAMDL(9) = 1
1398C mass dependence of slope parameters
1399 IPAMDL(10) = 0
1400C lepton-photon vertex 1
1401 IPAMDL(11) = 0
1402C lepton-photon vertex 2
1403 IPAMDL(12) = 0
1404C call by DPMJET
1405 IPAMDL(13) = 0
1406C method to sample x distributions
1407 IPAMDL(14) = 3
1408C energy-momentum check
1409 IPAMDL(15) = 1
1410C phase space correction for DPMJET interface
1411 IPAMDL(16) = 1
1412C fragment strings from projectile/target/central diff. separately
1413 IPAMDL(17) = 1
1414C method to construct strings for hard interactions
1415 IPAMDL(18) = 1
1416C method to construct strings for soft sea (pomeron cuts)
1417 IPAMDL(19) = 0
1418C method to construct strings in pomeron interactions
1419 IPAMDL(20) = 0
1420C soft color re-connection
1421 IPAMDL(21) = 0
1422C resummation of triple- and loop-Pomeron
1423 IPAMDL(24) = 1
1424C resummation of X iterated triple-Pomeron
1425 IPAMDL(25) = 1
1426C dimension of interpolation table for weights in hard scattering
1427 IPAMDL(30) = Max_tab_E
1428C dimension of interpolation table for pomeron cut distribution
1429 IPAMDL(31) = IEETA1
1430C number of cut soft pomerons (restriction by field dimension)
1431 IPAMDL(32) = IIMAX
1432C number of cut hard pomerons (restriction by field dimension)
1433 IPAMDL(33) = KKMAX
1434C tau pair production in direct photon-photon collisions
1435 IPAMDL(64) = 0
1436C currently chosen scales for hard scattering
1437C ATTENTION: IPAMDL(65-80) reserved for hard scattering!
1438 DO 15 I=1,16
1439 IPAMDL(64+I) = -99999
1440 15 CONTINUE
1441C scales to calculate alpha-s of matrix element
1442 IPAMDL(81) = 1
1443 IPAMDL(82) = 1
1444 IPAMDL(83) = 1
1445C scales to calculate alpha-s of initial state radiation
1446 IPAMDL(84) = 1
1447 IPAMDL(85) = 1
1448 IPAMDL(86) = 1
1449C scales to calculate alpha-s of final state radiation
1450 IPAMDL(87) = 1
1451 IPAMDL(88) = 1
1452 IPAMDL(89) = 1
1453C scales to calculate PDFs
1454 IPAMDL(90) = 1
1455 IPAMDL(91) = 1
1456 IPAMDL(92) = 1
1457C where to get the parameter sets from
1458 IPAMDL(99) = 1
1459C program PHO_ABORT for fatal errors (simulation of division by zero)
1460 IPAMDL(100) = 0
1461C initial state parton showers for all / hardest interaction(s)
1462 IPAMDL(101) = 1
1463C final state parton showers for all / hardest interaction(s)
1464 IPAMDL(102) = 1
1465C initial virtuality for ISR generation
1466 IPAMDL(109) = 1
1467C qqbar-gamma coupling in initial state showers
1468 IPAMDL(110) = 1
1469C generation of time-like showers during ISR
1470 IPAMDL(111) = 1
1471C reweighting of multiple soft contributions for virtual photons
1472 IPAMDL(114) = 1
1473C reweighting / use photon virtuality in photon PDF calculations
1474 IPAMDL(115) = 0
1475C use full QPM model incl. interference terms (direct part in gam-gam)
1476 IPAMDL(116) = 0
1477C matching sigma_tot to F2 as given by parton density at high Q2
1478 IPAMDL(117) = 1
1479C use virtuality of target in F2 calculations (two-gamma only)
1480 IPAMDL(118) = 1
1481C calculation of alpha_em
1482 IPAMDL(120) = 1
1483C strict pt cutoff for gamma-gamma events
1484 IPAMDL(121) = 0
1485C photon virtuality sampled in photon flux approximations
1486 IPAMDL(174) = 1
1487C photon-pomeron: 0,1,2: both,left,right photon emission
1488 IPAMDL(175) = 0
1489C keep full history information in PHOJET-JETSET interface
1490 IPAMDL(178) = 1
1491C max. number of conservation law violations allowed in one run
1492 IPAMDL(179) = 20
1493C selection of soft X values
1494C max. iteration number in PHO_SELSXS
1495 IPAMDL(180) = 50
1496C max. iteration number in PHO_SELSXR
1497 IPAMDL(181) = 200
1498C max. iteration number in PHO_SELSX2
1499 IPAMDL(182) = 100
1500C max. iteration number in PHO_SELSXI
1501 IPAMDL(183) = 50
1502
1503C initialize /PROBAB/
1504 IEEMAX = IEETA1
1505 IMAX = IIMAX
1506 KMAX = KKMAX
1507
1508 DO 20 I=1,30
1509 PARMDL(300+I) = -100000.D0
1510 20 CONTINUE
1511C initialize /POHDRN/
1512 QMASS(1) = PARMDL(151)
1513 QMASS(2) = PARMDL(152)
1514 QMASS(3) = PARMDL(153)
1515 QMASS(4) = PARMDL(154)
1516 QMASS(5) = PARMDL(155)
1517 QMASS(6) = PARMDL(156)
1518 BET = 8.D0
1519 PCOUDI = 0.D0
1520 VALPRG(1) = 1.D0
1521 VALPRG(2) = 1.D0
1522C number of light flavours (quarks treated as massless)
1523 NFS = 4
1524C initialize /POCUT1/
1525 PTCUT(1) = PARMDL(36)
1526 PTCUT(2) = PARMDL(37)
1527 PTCUT(3) = PARMDL(38)
1528 PTCUT(4) = PARMDL(39)
1529 PSOMIN = 0.D0
1530 XSOMIN = 0.D0
1531C initialize /POHAPA/
1532 NFbeta = 4
1533 NF = 4
1534 BQCD(1) = PI4/(11.D0-(2.D0/3.D0)*3)
1535 BQCD(2) = PI4/(11.D0-(2.D0/3.D0)*4)
1536 BQCD(3) = PI4/(11.D0-(2.D0/3.D0)*5)
1537 BQCD(4) = PI4/(11.D0-(2.D0/3.D0)*6)
1538C initialize /POGAUP/
1539 NGAUP1 = 12
1540 NGAUP2 = 12
1541 NGAUET = 16
1542 NGAUIN = 12
1543 NGAUSO = 96
1544C initialize //
1545 DO 30 I=1,100
1546 IDEB(I) = 0
1547 30 CONTINUE
1548C initialize /PROCES/
1549 DO 35 I=1,11
1550 IPRON(I,1) = 1
1551 35 CONTINUE
1552
1553C DPMJET default: no elastic scattering
1554 IPRON(2,1) = 0
1555
1556 DO 36 K=2,4
1557 DO 37 I=2,11
1558 IPRON(I,K) = 0
1559 37 CONTINUE
1560 IPRON(1,K) = 1
1561 IPRON(8,K) = 1
1562 36 CONTINUE
1563C initialize /POSVDM/
1564 TWOPIM = 0.28D0
1565 RMIN(1) = 0.285D0
1566 RMIN(2) = 0.45D0
1567 RMIN(3) = 1.D0
1568 RMIN(4) = TWOPIM
1569 VMAS(1) = 0.770D0
1570 VMAS(2) = 0.787D0
1571 VMAS(3) = 1.02D0
1572 VMAS(4) = TWOPIM
1573 GAMM(1) = 0.155D0
1574 GAMM(2) = 0.01D0
1575 GAMM(3) = 0.0045D0
1576 GAMM(4) = 1.D0
1577 RMAX(1) = VMAS(1)+TWOPIM
1578 RMAX(2) = VMAS(2)+TWOPIM
1579 RMAX(3) = VMAS(3)+TWOPIM
1580 RMAX(4) = VMAS(1)+TWOPIM
1581 VMSL(1) = 11.D0
1582 VMSL(2) = 10.D0
1583 VMSL(3) = 6.D0
1584 VMSL(4) = 4.D0
1585 VMFA(1) = 0.0033D0
1586 VMFA(2) = 0.00036D0
1587 VMFA(3) = 0.0002D0
1588 VMFA(4) = 0.0002D0
1589C initialize /PODGL1/
1590 Q2MISR(1) = PARMDL(36)**2
1591 Q2MISR(2) = PARMDL(36)**2
1592 PMISR(1) = 1.D0
1593 PMISR(2) = 1.D0
1594 ZMISR(1) = 0.001D0
1595 ZMISR(2) = 0.001D0
1596 AL2ISR(1) = 0.046D0
1597 AL2ISR(2) = 0.046D0
1598 NFSISR = 4
1599C initialize /POPISR/
1600 DO 40 I=1,50
1601 IPOISR(1,2,I) = 0
1602 IPOISR(2,2,I) = 0
1603 40 CONTINUE
1604C initialize /POHPRO/
1605 PROC(0) = 'sum over processes'
1606 PROC(1) = 'G +G --> G +G '
1607 PROC(2) = 'Q +QB --> G +G '
1608 PROC(3) = 'G +Q --> G +Q '
1609 PROC(4) = 'G +G --> Q +QB '
1610 PROC(5) = 'Q +QB --> Q +QB '
1611 PROC(6) = 'Q +QB --> QP +QBP'
1612 PROC(7) = 'Q +Q --> Q +Q '
1613 PROC(8) = 'Q +QP --> Q +QP '
1614 PROC(9) = 'resolved processes'
1615 PROC(10) = 'gam+Q --> G +Q '
1616 PROC(11) = 'gam+G --> Q +QB '
1617 PROC(12) = 'Q +gam--> G +Q '
1618 PROC(13) = 'G +gam--> Q +QB '
1619 PROC(14) = 'gam+gam--> Q +QB '
1620 PROC(15) = 'direct processes '
1621 PROC(16) = 'gam+gam--> l+ +l- '
1622
1623C initialize /POHRCS/
1624 do M=1,Max_pro_2
1625 HWgx(M) = 0.D0
1626 HSig(M) = 0.D0
1627 Hdpt(M) = 0.D0
1628 enddo
1629 DO I=0,4
1630 DO M=-1,Max_pro_2
1631C switch all hard subprocesses on
1632 MH_pro_on(M,I) = 1
1633C reset all counters
1634 MH_tried(M,I) = 0
1635 MH_acc_1(M,I) = 0
1636 MH_acc_2(M,I) = 0
1637 ENDDO
1638 MH_pro_on(16,I) = 0
1639 ENDDO
1640
1641C initialize /POHTAB/
1642 do I=0,4
1643 IH_Ecm_up(I) = 0
1644 IH_Q2a_up(I) = 0
1645 IH_Q2b_up(I) = 0
1646 HEcm_tab(1,I) = 0.D0
1647 enddo
1648 HEcm_last = 0.D0
1649 IHa_last = 0.D0
1650 IHb_last = 0.D0
1651
1652C initialize /POFSRC/
1653 IGHEL(1) = -1
1654 IGHEL(2) = -1
1655C initialize /LEPCUT/
1656 ECMIN = 5.D0
1657 ECMAX = 1.D+30
1658 EEMIN1 = 1.D0
1659 EEMIN2 = 1.D0
1660 YMAX1 = -1.D0
1661 YMAX2 = -1.D0
1662 THMIN1 = 0.D0
1663 THMAX1 = PI
1664 THMIN2 = 0.D0
1665 THMAX2 = PI
1666 ITAG1 = 1
1667 ITAG2 = 1
1668C initialize /POWGHT/
1669 DO 70 I=1,20
1670 HSWCUT(I) = 0.D0
1671 ISWCUT(I) = 0
1672 70 CONTINUE
1673 EVWGHT(1) = 1.D0
1674 IVWGHT(1) = 0
1675 SIGGEN(1) = 0.D0
1676 SIGGEN(2) = 0.D0
1677 SIGGEN(3) = 0.D0
1678 SIGGEN(4) = 0.D0
1679
1680 END
1681
1682*$ CREATE PHO_PARDAT.FOR
1683*COPY PHO_PARDAT
1684CDECK ID>, PHO_PARDAT
1685 SUBROUTINE PHO_PARDAT
1686C***********************************************************************
1687C
1688C particle data (based on 1996 PDG naming scheme and data tables)
1689C
1690C***********************************************************************
1691
1692 IMPLICIT NONE
1693
1694 SAVE
1695
1696C input/output channels
1697 INTEGER LI,LO
1698 COMMON /POINOU/ LI,LO
1699C 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
1706C 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
1711C 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)
1724C 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
1730C external functions
1731
1732 integer ipho_pdg2id
1733 double precision pho_pmass
1734
1735C 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
2278C 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
2292C 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
2298C 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
2306C 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
2320C copy secondary particles
2321C (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
2339C 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
2350C 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
2356C 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)
2370C 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
2390C 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
2397C 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
2404C 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
2411C 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)
2420C 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
2435C 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
2443C 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)
2457C decuplet baryons
2458 K10 = id_b10_list(i,j,l)
2459 if(K10.ne.0) then
2460 AM1 = xm_list(K10)
2461 else
2462 AM1 = pho_pmass(i,3)+pho_pmass(j,3)+pho_pmass(l,3)
2463 endif
2464 L10 = id_b10_list(ii,jj,l)
2465 if(L10.ne.0) then
2466 AM2 = xm_list(L10)
2467 else
2468 AM2 = pho_pmass(ii,3)+pho_pmass(jj,3)+pho_pmass(l,3)
2469 endif
2470 AM102 = MIN(AM102, AM1+AM2)
2471 enddo
2472 xm_bb82_list(i,j,ii,jj) = AM82
2473 xm_bb102_list(i,j,ii,jj) = AM102
2474 enddo
2475 enddo
2476 enddo
2477 enddo
2478
2479 END
2480
2481*$ CREATE PHO_PRESEL.FOR
2482*COPY PHO_PRESEL
2483CDECK ID>, PHO_PRESEL
2484 SUBROUTINE PHO_PRESEL(MODE,IREJ)
2485C**********************************************************************
2486C
2487C user specific function to pre-select events during generation
2488C
2489C input: MODE 5 electron and photon kinematics
2490C 10 process and number of cut Pomerons
2491C 15 partons without construction of strings
2492C 20 partons assigned to strings
2493C 25 after fragmentation, complete final state
2494C
2495C output: IREJ 0 event accepted
2496C 50 event rejected
2497C
2498C**********************************************************************
2499 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2500 SAVE
2501
2502C input/output channels
2503 INTEGER LI,LO
2504 COMMON /POINOU/ LI,LO
2505C event debugging information
2506 INTEGER NMAXD
2507 PARAMETER (NMAXD=100)
2508 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2509 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2510 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2511 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2512
2513C standard particle data interface
2514 INTEGER NMXHEP
2515
2516 PARAMETER (NMXHEP=4000)
2517
2518 INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
2519 DOUBLE PRECISION PHEP,VHEP
2520 COMMON /POEVT1/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
2521 & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),
2522 & VHEP(4,NMXHEP)
2523C extension to standard particle data interface (PHOJET specific)
2524 INTEGER IMPART,IPHIST,ICOLOR
2525 COMMON /POEVT2/ IMPART(NMXHEP),IPHIST(2,NMXHEP),ICOLOR(2,NMXHEP)
2526
2527C global event kinematics and particle IDs
2528 INTEGER IFPAP,IFPAB
2529 DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
2530 COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
2531C gamma-lepton or gamma-hadron vertex information
2532 INTEGER IGHEL,IDPSRC,IDBSRC
2533 DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
2534 & RADSRC,AMSRC,GAMSRC
2535 COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
2536 & GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
2537 & IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
2538C hard scattering data
2539 INTEGER MSCAHD
2540 PARAMETER ( MSCAHD = 50 )
2541 INTEGER LSCAHD,LSC1HD,LSIDX,
2542 & NINHD,N0INHD,NIVAL,N0IVAL,NOUTHD,NBRAHD,NPROHD
2543 DOUBLE PRECISION PPH,PTHD,ETAHD,Q2SCA,PDFVA,XHD,VHD,X0HD
2544 COMMON /POHSLT/ LSCAHD,LSC1HD,LSIDX(MSCAHD),
2545 & PPH(8*MSCAHD,2),PTHD(MSCAHD),ETAHD(MSCAHD,2),
2546 & Q2SCA(MSCAHD,2),PDFVA(MSCAHD,2),
2547 & XHD(MSCAHD,2),VHD(MSCAHD),X0HD(MSCAHD,2),
2548 & NINHD(MSCAHD,2),N0INHD(MSCAHD,2),
2549 & NIVAL(MSCAHD,2),N0IVAL(MSCAHD,2),
2550 & NOUTHD(MSCAHD,2),NBRAHD(MSCAHD,2),NPROHD(MSCAHD)
2551C event weights and generated cross section
2552 INTEGER IPOWGC,ISWCUT,IVWGHT
2553 DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2554 COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2555 & IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2556
2557 IREJ = 0
2558
2559* XBJ = GQ2(2)/(GGECM**2+GQ2(2))
2560* IF(XBJ.LT.0.002D0) IREJ = 1
2561
2562 END
2563
2564*$ CREATE PHO_FIXCOL.FOR
2565*COPY PHO_FIXCOL
2566CDECK ID>, PHO_FIXCOL
2567 SUBROUTINE PHO_FIXCOL(E1,E2,THETA,PHI,NEV)
2568C**********************************************************************
2569C
2570C interface to call PHOJET (fixed energy run) with
2571C collider kinematics
2572C
2573C equivalen photon approximation to get photon flux
2574C
2575C input: NEV number of events to generate
2576C THETA azimuthal angle (micro radians)
2577C PHI beam crossing angle
2578C (with respect to x, in degrees)
2579C E1 energy of particle 1 (+z direction, GeV)
2580C E2 energy of particle 2 (-z direction, GeV)
2581C
2582C note: particle types have to be specified before
2583C with PHO_SETPAR
2584C
2585C**********************************************************************
2586 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2587 SAVE
2588
2589 PARAMETER(TWOPI=6.283185307D0,BOG=TWOPI/360.0D0)
2590
2591C input/output channels
2592 INTEGER LI,LO
2593 COMMON /POINOU/ LI,LO
2594C event debugging information
2595 INTEGER NMAXD
2596 PARAMETER (NMAXD=100)
2597 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2598 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2599 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2600 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2601C general process information
2602 INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
2603 COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
2604C global event kinematics and particle IDs
2605 INTEGER IFPAP,IFPAB
2606 DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
2607 COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
2608C model switches and parameters
2609 CHARACTER*8 MDLNA
2610 INTEGER ISWMDL,IPAMDL
2611 DOUBLE PRECISION PARMDL
2612 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
2613C nucleon-nucleus / nucleus-nucleus interface to DPMJET
2614 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
2615 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
2616 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
2617 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
2618C integration precision for hard cross sections (obsolete)
2619 INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2620 COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2621C event weights and generated cross section
2622 INTEGER IPOWGC,ISWCUT,IVWGHT
2623 DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2624 COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2625 & IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2626
2627 DIMENSION P1(4),P2(4)
2628
2629C remnant initialization (only needed for DPMJET)
2630 ISAVP1 = IFPAP(1)
2631 ISAVB1 = IFPAB(1)
2632 IF(IFPAP(1).EQ.81) THEN
2633 IFPAP(1) = IDEQP(1)
2634 IFPAB(1) = IDEQB(1)
2635 ENDIF
2636 ISAVP2 = IFPAP(2)
2637 ISAVB2 = IFPAB(2)
2638 IF(IFPAP(2).EQ.82) THEN
2639 IFPAP(2) = IDEQP(2)
2640 IFPAB(2) = IDEQB(2)
2641 ENDIF
2642 PMASS1 = PHO_PMASS(IFPAB(1),0)-SQRT(PVIRT(1))
2643 PMASS2 = PHO_PMASS(IFPAB(2),0)-SQRT(PVIRT(2))
2644 PP1 = SQRT(E1**2-PMASS1**2)
2645 PP2 = SQRT(E2**2-PMASS2**2)
2646C beam crossing angle
2647 TH = 1.D-6*THETA/2.D0
2648 PH = PHI*BOG
2649 P1(1) = PP1*SIN(TH)*COS(PH)
2650 P1(2) = PP1*SIN(TH)*SIN(PH)
2651 P1(3) = PP1*COS(TH)
2652 P1(4) = E1
2653 P2(1) = PP2*SIN(TH)*COS(PH)
2654 P2(2) = PP2*SIN(TH)*SIN(PH)
2655 P2(3) = -PP2*COS(TH)
2656 P2(4) = E2
2657 CALL PHO_EVENT(-1,P1,P2,SIGMAX,IREJ)
2658 IFPAP(1) = ISAVP1
2659 IFPAB(1) = ISAVB1
2660 IFPAP(2) = ISAVP2
2661 IFPAB(2) = ISAVB2
2662 ITRY = 0
2663 CALL PHO_PHIST(-1,SIGMAX)
2664 CALL PHO_LHIST(-1,SIGMAX)
2665C test of DPMJET interface (default is IPAMDL(13)=0)
2666 if(IPAMDL(13).gt.0) then
2667 MODE = IPAMDL(13)
2668 IPAMDL(13) = 0
2669 else
2670 MODE = 1
2671 endif
2672C main generation loop
2673 DO 50 I=1,NEV
2674 55 CONTINUE
2675 ITRY = ITRY+1
2676 CALL PHO_EVENT(MODE,P1,P2,SIGCUR,IREJ)
2677 IF(IREJ.NE.0) GOTO 55
2678 CALL PHO_PHIST(1,HSWGHT(0))
2679 CALL PHO_LHIST(1,HSWGHT(0))
2680 50 CONTINUE
2681
2682 IF(NEV.GT.0) THEN
2683 SIGMAX = SIGMAX*DBLE(NEV)/DBLE(ITRY)
2684 WRITE(LO,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
2685 & '=========================================================',
2686 & ' ***** simulated cross section: ',SIGMAX,' mb *****',
2687 & '========================================================='
2688 CALL PHO_EVENT(-2,P1,P2,SIGCUR,IREJ)
2689 CALL PHO_PHIST(-2,SIGMAX)
2690 CALL PHO_LHIST(-2,SIGMAX)
2691 ELSE
2692 WRITE(LO,'(1X,A,I5)') 'POFCOL: no events simulated',NEV
2693 ENDIF
2694
2695 END
2696
2697*$ CREATE PHO_FIXLAB.FOR
2698*COPY PHO_FIXLAB
2699CDECK ID>, PHO_FIXLAB
2700 SUBROUTINE PHO_FIXLAB(PLAB,NEV)
2701C**********************************************************************
2702C
2703C interface to call PHOJET (fixed energy run) with
2704C LAB kinematics (second particle as target)
2705C
2706C equivalent photon approximation to get photon flux
2707C
2708C input: NEV number of events to generate
2709C PLAB LAB momentum of particle 1
2710C
2711C note: particle types have to be specified before
2712C with PHO_SETPAR
2713C
2714C**********************************************************************
2715 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2716 SAVE
2717
2718C input/output channels
2719 INTEGER LI,LO
2720 COMMON /POINOU/ LI,LO
2721C event debugging information
2722 INTEGER NMAXD
2723 PARAMETER (NMAXD=100)
2724 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2725 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2726 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2727 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2728C general process information
2729 INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
2730 COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
2731C global event kinematics and particle IDs
2732 INTEGER IFPAP,IFPAB
2733 DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
2734 COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
2735C model switches and parameters
2736 CHARACTER*8 MDLNA
2737 INTEGER ISWMDL,IPAMDL
2738 DOUBLE PRECISION PARMDL
2739 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
2740C nucleon-nucleus / nucleus-nucleus interface to DPMJET
2741 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
2742 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
2743 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
2744 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
2745C integration precision for hard cross sections (obsolete)
2746 INTEGER NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2747 COMMON /POGAUP/ NGAUP1,NGAUP2,NGAUET,NGAUIN,NGAUSO
2748C event weights and generated cross section
2749 INTEGER IPOWGC,ISWCUT,IVWGHT
2750 DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2751 COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2752 & IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2753
2754 DIMENSION P1(4),P2(4)
2755
2756C remnant initialization (only needed for DPMJET)
2757 SPCM = PLAB
2758 ISAVP1 = IFPAP(1)
2759 ISAVB1 = IFPAB(1)
2760 IF(IFPAP(1).EQ.81) THEN
2761 IFPAP(1) = IDEQP(1)
2762 IFPAB(1) = IDEQB(1)
2763 ENDIF
2764 ISAVP2 = IFPAP(2)
2765 ISAVB2 = IFPAB(2)
2766 IF(IFPAP(2).EQ.82) THEN
2767 IFPAP(2) = IDEQP(2)
2768 IFPAB(2) = IDEQB(2)
2769 ENDIF
2770C get momenta in LAB system
2771 PMASS1 = PHO_PMASS(IFPAB(1),0)**2-PVIRT(1)
2772 PMASS2 = PHO_PMASS(IFPAB(2),0)**2-PVIRT(2)
2773 IF(PMASS2.LT.0.1D0) THEN
2774 WRITE(LO,'(/1X,2A,2I7)') 'PHO_FIXLAB:ERROR: ',
2775 & 'no LAB system possible',IFPAB(1),IFPAB(2)
2776 ELSE
2777 P1(1) = 0.D0
2778 P1(2) = 0.D0
2779 P1(3) = PLAB
2780 P1(4) = SQRT(PMASS1+PLAB**2)
2781 P2(1) = 0.D0
2782 P2(2) = 0.D0
2783 P2(3) = 0.D0
2784 P2(4) = SQRT(PMASS2)
2785 CALL PHO_EVENT(-1,P1,P2,SIGMAX,IREJ)
2786 IFPAP(1) = ISAVP1
2787 IFPAB(1) = ISAVB1
2788 IFPAP(2) = ISAVP2
2789 IFPAB(2) = ISAVB2
2790 ITRY = 0
2791 CALL PHO_PHIST(-1,SIGMAX)
2792 CALL PHO_LHIST(-1,SIGMAX)
2793C event generation loop
2794 DO 40 I=1,NEV
2795 45 CONTINUE
2796 ITRY = ITRY+1
2797 CALL PHO_EVENT(1,P1,P2,SIGCUR,IREJ)
2798 IF(IREJ.NE.0) GOTO 45
2799 CALL PHO_LHIST(1,HSWGHT(0))
2800
2801 CALL PHO_PHIST(10,HSWGHT(0))
2802
2803 40 CONTINUE
2804 IF(NEV.GT.0) THEN
2805 SIGMAX = SIGMAX*DBLE(NEV)/DBLE(ITRY)
2806 WRITE(LO,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
2807 & '=========================================================',
2808 & ' ***** simulated cross section: ',SIGMAX,' mb *****',
2809 & '========================================================='
2810 CALL PHO_EVENT(-2,P1,P2,SIGCUR,IREJ)
2811 CALL PHO_PHIST(-2,SIGMAX)
2812 CALL PHO_LHIST(-2,SIGMAX)
2813 ELSE
2814 WRITE(LO,'(1X,A,I5)')
2815 & 'PHO_FIXLAB: no events simulated',NEV
2816 ENDIF
2817 ENDIF
2818
2819 END
2820
2821*$ CREATE PHO_GPHERA.FOR
2822*COPY PHO_GPHERA
2823CDECK ID>, PHO_GPHERA
2824 SUBROUTINE PHO_GPHERA(NEVENT,EE1,EE2)
2825C**********************************************************************
2826C
2827C interface to call PHOJET (variable energy run) with
2828C HERA kinematics, photon as particle 2
2829C
2830C equivalent photon approximation to get photon flux
2831C
2832C input: NEVENT number of events to generate
2833C EE1 proton energy (LAB system)
2834C EE2 electron energy (LAB system)
2835C from /POFCUT/:
2836C YMIN2 lower limit of Y
2837C (energy fraction taken by photon from electron)
2838C YMAX2 upper limit of Y
2839C Q2MIN2 lower limit of photon virtuality
2840C Q2MAX2 upper limit of photon virtuality
2841C
2842C**********************************************************************
2843 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2844 SAVE
2845
2846 PARAMETER ( DEPS = 1.D-10,
2847 & PI = 3.14159265359D0 )
2848
2849C input/output channels
2850 INTEGER LI,LO
2851 COMMON /POINOU/ LI,LO
2852C event debugging information
2853 INTEGER NMAXD
2854 PARAMETER (NMAXD=100)
2855 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
2856 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2857 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
2858 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
2859C model switches and parameters
2860 CHARACTER*8 MDLNA
2861 INTEGER ISWMDL,IPAMDL
2862 DOUBLE PRECISION PARMDL
2863 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
2864C photon flux kinematics and cuts
2865 DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
2866 & YMIN1,YMAX1,YMIN2,YMAX2,
2867 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
2868 & THMIN1,THMAX1,THMIN2,THMAX2
2869 INTEGER ITAG1,ITAG2
2870 COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
2871 & YMIN1,YMAX1,YMIN2,YMAX2,
2872 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
2873 & THMIN1,THMAX1,THMIN2,THMAX2,
2874 & ITAG1,ITAG2
2875C gamma-lepton or gamma-hadron vertex information
2876 INTEGER IGHEL,IDPSRC,IDBSRC
2877 DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
2878 & RADSRC,AMSRC,GAMSRC
2879 COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
2880 & GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
2881 & IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
2882C nucleon-nucleus / nucleus-nucleus interface to DPMJET
2883 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
2884 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
2885 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
2886 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
2887C event weights and generated cross section
2888 INTEGER IPOWGC,ISWCUT,IVWGHT
2889 DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
2890 COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
2891 & IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
2892
2893 DIMENSION P1(4),P2(4)
2894
2895 WRITE(LO,'(//1X,A,I10)') 'PHO_GPHERA: events to process',NEVENT
2896C assign particle momenta according to HERA kinematics
2897C proton data
2898 PROM = PHO_PMASS(2212,1)
2899 PROM2 = PROM**2
2900 IDPSRC(1) = 0
2901 IDBSRC(1) = 0
2902C electron data
2903 ELEM = 0.512D-03
2904 ELEM2 = ELEM**2
2905 AMSRC(2) = ELEM
2906 IDPSRC(2) = 11
2907 IDBSRC(2) = ipho_pdg2id(11)
2908C
2909 Q2MIN = Q2MIN2
2910 Q2MAX = Q2MAX2
2911C
2912 XIMAX = LOG(YMAX2)
2913 XIMIN = LOG(YMIN2)
2914 XIDEL = XIMAX-XIMIN
2915C
2916 IF(Q2MIN.GT.ELEM2*YMIN2**2/(1.D0-YMIN2))
2917 & WRITE(LO,'(/1X,A,1P2E11.4)')
2918 & 'PHO_GPHERA: lower Q2 cutoff larger than kin. limit:',
2919 & Q2MIN,ELEM2*YMIN2**2/(1.D0-YMIN2)
2920C
2921 Max_tab = 50
2922 DELLY = LOG(YMAX2/YMIN2)/DBLE(Max_tab-1)
2923 FLUXT = 0.D0
2924 FLUXL = 0.D0
2925 IF(IDEB(30).GE.1) WRITE(LO,'(1X,A,I5)')
2926 & 'PHO_GPHERA: table of photon flux (trans/long)',Max_tab
2927 DO 100 I=1,Max_tab
2928 Y = EXP(XIMIN+DELLY*DBLE(I-1))
2929 Q2LOW = MAX(Q2MIN,ELEM2*Y**2/(1.D0-Y))
2930 FFT = ((1.D0+(1.D0-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
2931 & -2.D0*ELEM2*Y*(1.D0/Q2LOW-1.D0/Q2MAX))/(2.D0*PI*137.D0)
2932 FFL = 2.D0*(1.D0-Y)/Y*LOG(Q2MAX/Q2LOW)/(2.D0*PI*137.D0)
2933 FLUXT = FLUXT + Y*FFT
2934 FLUXL = FLUXL + Y*FFL
2935 IF(IDEB(30).GE.1) WRITE(LO,'(5X,1P3E14.4)') Y,FFT,FFL
2936 100 CONTINUE
2937 FLUXT = FLUXT*DELLY
2938 FLUXL = FLUXL*DELLY
2939 IF(IDEB(30).GE.1) WRITE(LO,'(1X,A,1P2E12.4)')
2940 & 'PHO_GPHERA: integrated flux (trans./long.):',FLUXT,FLUXL
2941C
2942 AY = 0.D0
2943 AY2 = 0.D0
2944 YY = YMIN2
2945 Q2LOW = MAX(Q2MIN,ELEM2*YY**2/(1.D0-YY))
2946 WGMAX = (1.D0+(1.D0-YY)**2)*LOG(Q2MAX/Q2LOW)
2947 & -2.D0*ELEM2*YY*(1.D0/Q2LOW-1.D0/Q2MAX)*YY
2948 IF(ISWMDL(10).GE.2) WGMAX = WGMAX+2.D0*(1.D0-YY)*LOG(Q2MAX/Q2LOW)
2949C
2950C initialization of PHOJET at upper energy limit
2951C proton momentum
2952 P1(1) = 0.D0
2953 P1(2) = 0.D0
2954 P1(3) = SQRT(EE1**2-PROM2+DEPS)
2955 P1(4) = EE1
2956C photon momentum
2957 EGAM = YMAX2*EE2
2958 P2(1) = 0.D0
2959 P2(2) = 0.D0
2960 P2(3) = -EGAM
2961 P2(4) = EGAM
2962C sum of both photon polarizations
2963 IGHEL(2) = -1
2964C
2965 CALL PHO_SETPAR(1,2212,0,0.D0)
2966 CALL PHO_SETPAR(2,22,0,0.D0)
2967 CALL PHO_EVENT(-1,P1,P2,SIGMAX,IREJ)
2968 CALL PHO_PHIST(-1,SIGMAX)
2969 CALL PHO_LHIST(-1,SIGMAX)
2970C
2971C generation of events, flux calculation
2972
2973 ECMIN2 = ECMIN**2
2974 ECMAX2 = ECMAX**2
2975 AY = 0.D0
2976 AY2 = 0.D0
2977 Q22MIN = 1.D30
2978 Q22AVE = 0.D0
2979 Q22AV2 = 0.D0
2980 Q22MAX = 0.D0
2981 AN2MIN = 1.D30
2982 AN2MAX = 0.D0
2983 YY2MIN = 1.D30
2984 YY2MAX = 0.D0
2985 NITER = NEVENT
2986 ITRY = 0
2987 ITRW = 0
2988 DO 200 I=1,NITER
2989 150 CONTINUE
2990C sample y
2991 ITRY = ITRY+1
2992 175 CONTINUE
2993 ITRW = ITRW+1
2994 YY = EXP(XIDEL*DT_RNDM(AY)+XIMIN)
2995 IF(ISWMDL(10).GE.2) THEN
2996 YEFF = 1.D0+(1.D0-YY)**2+2.D0*(1.D0-YY)
2997 ELSE
2998 YEFF = 1.D0+(1.D0-YY)**2
2999 ENDIF
3000 Q2LOW = MAX(Q2MIN,ELEM2*YY**2/(1.D0-YY))
3001 Q2LOG = LOG(Q2MAX/Q2LOW)
3002 WGH = YEFF*Q2LOG-2.D0*ELEM2*YY**2*(1.D0/Q2LOW-1.D0/Q2MAX)
3003 IF(WGMAX.LT.WGH) THEN
3004 WRITE(LO,'(1X,A,3E12.5)')
3005 & 'PHO_GPHERA: inconsistent weight:',YY,WGMAX,WGH
3006 ENDIF
3007 IF(DT_RNDM(AY2)*WGMAX.GT.WGH) GOTO 175
3008C sample Q2
3009 IF(IPAMDL(174).EQ.1) THEN
3010 185 CONTINUE
3011 Q2 = Q2LOW*EXP(Q2LOG*DT_RNDM(YY))
3012 WEIGHT = (YEFF-2.D0*ELEM2*YY**2/Q2)/YEFF
3013 IF(WEIGHT.LT.DT_RNDM(Q2)) GOTO 185
3014 ELSE
3015 Q2 = Q2LOW
3016 ENDIF
3017C
3018
3019C incoming electron
3020 PINI(1,2) = 0.D0
3021 PINI(2,2) = 0.D0
3022 PINI(3,2) = -EE2
3023 PINI(4,2) = EE2
3024 PINI(5,2) = 0.D0
3025C outgoing electron
3026 YQ2 = SQRT((1.D0-YY)*Q2)
3027 Q2E = Q2/(4.D0*EE2)
3028 E1Y = EE2*(1.D0-YY)
3029 CALL PHO_SFECFE(SIF,COF)
3030 PFIN(1,2) = YQ2*COF
3031 PFIN(2,2) = YQ2*SIF
3032 PFIN(3,2) = -E1Y+Q2E
3033 PFIN(4,2) = E1Y+Q2E
3034 PFIN(5,2) = 0.D0
3035C set /POFSRC/
3036 GYY(2) = YY
3037 GQ2(2) = Q2
3038C polar angle
3039 PFTHE(2) = ACOS(PFIN(3,2)/PFIN(4,2))
3040C electron tagger
3041 IF(PFIN(4,2).GT.EEMIN2) THEN
3042 IF((PFTHE(2).LT.THMIN2).OR.(PFTHE(2).GT.THMAX2)) GOTO 175
3043 ENDIF
3044C azimuthal angle
3045 PFPHI(2) = ATAN2(COF,SIF)
3046C photon momentum
3047 P2(1) = -PFIN(1,2)
3048 P2(2) = -PFIN(2,2)
3049 P2(3) = PINI(3,2)-PFIN(3,2)
3050 P2(4) = PINI(4,2)-PFIN(4,2)
3051C proton momentum
3052 P1(1) = 0.D0
3053 P1(2) = 0.D0
3054 P1(3) = SQRT(EE1**2-PROM2)
3055 P1(4) = EE1
3056C ECMS cut
3057 GGECM = (P1(4)+P2(4))**2-(P1(1)+P2(1))**2
3058 & -(P1(2)+P2(2))**2-(P1(3)+P2(3))**2
3059 IF((GGECM.LT.ECMIN2).OR.(GGECM.GT.ECMAX2)) GOTO 175
3060 GGECM = SQRT(GGECM)
3061C
3062 PGAM(1,2) = P2(1)
3063 PGAM(2,2) = P2(2)
3064 PGAM(3,2) = P2(3)
3065 PGAM(4,2) = P2(4)
3066 PGAM(5,2) = -SQRT(Q2)
3067C photon helicity
3068 IF(ISWMDL(10).GE.2) THEN
3069 WGH = YEFF-2.D0*ELEM2*YY**2/Q2
3070 WGHL = 2.D0*(1-YY)
3071 IF(DT_RNDM(YY).GE.WGHL/WGH) THEN
3072 IGHEL(2) = 1
3073 ELSE
3074 IGHEL(2) = 0
3075 ENDIF
3076 ELSE
3077 IGHEL(2) = -1
3078 ENDIF
3079C user cuts
3080 CALL PHO_PRESEL(5,IREJ)
3081 IF(IREJ.NE.0) GOTO 175
3082C event generation
3083 CALL PHO_EVENT(1,P1,P2,SIGCUR,IREJ)
3084 IF(IREJ.NE.0) GOTO 150
3085
3086C statistics
3087 AY = AY+YY
3088 AY2 = AY2+YY*YY
3089 YY2MIN = MIN(YY2MIN,YY)
3090 YY2MAX = MAX(YY2MAX,YY)
3091 Q22MIN = MIN(Q22MIN,Q2)
3092 Q22MAX = MAX(Q22MAX,Q2)
3093 Q22AVE = Q22AVE+Q2
3094 Q22AV2 = Q22AV2+Q2*Q2
3095 AN2MIN = MIN(AN2MIN,PFTHE(2))
3096 AN2MAX = MAX(AN2MAX,PFTHE(2))
3097C histograms
3098 CALL PHO_PHIST(1,HSWGHT(0))
3099 CALL PHO_LHIST(1,HSWGHT(0))
3100 200 CONTINUE
3101C
3102 WGY = WGMAX*DBLE(ITRY)/DBLE(ITRW)/(137.D0*2.D0*PI)
3103 WGY = WGY*LOG(YMAX2/YMIN2)
3104 AY = AY/DBLE(NITER)
3105 AY2 = AY2/DBLE(NITER)
3106 DAY = SQRT((AY2-AY**2)/DBLE(NITER))
3107 Q22AVE = Q22AVE/DBLE(NITER)
3108 Q22AV2 = Q22AV2/DBLE(NITER)
3109 Q22AV2 = SQRT((Q22AV2-Q22AVE**2)/DBLE(NITER))
3110 WEIGHT = WGY*SIGMAX*DBLE(NITER)/DBLE(ITRY)
3111C output of histograms
3112 WRITE(LO,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
3113 &'=========================================================',
3114 &' ***** simulated cross section: ',WEIGHT,' mb *****',
3115 &'========================================================='
3116 WRITE(LO,'(//1X,A,3I10)')
3117 & 'PHO_GPHERA:SUMMARY:NITER,ITRY,ITRW',NITER,ITRY,ITRW
3118 WRITE(LO,'(1X,A,1P2E12.4)') 'EFFECTIVE WEIGHT (FLUX,TOTAL)',
3119 & WGY,WEIGHT
3120 WRITE(LO,'(1X,A,1P2E12.4)') 'AVERAGE Y,DY ',AY,DAY
3121 WRITE(LO,'(1X,A,1P2E12.4)') 'SAMPLED Y RANGE PHOTON ',
3122 & YY2MIN,YY2MAX
3123 WRITE(LO,'(1X,A,1P2E12.4)') 'AVERAGE Q2,DQ2 ',
3124 & Q22AVE,Q22AV2
3125 WRITE(LO,'(1X,A,1P2E12.4)') 'SAMPLED Q2 RANGE PHOTON ',
3126 & Q22MIN,Q22MAX
3127 WRITE(LO,'(1X,A,1P4E12.4)') 'SAMPLED THETA RANGE ELECTRON ',
3128 & AN2MIN,AN2MAX,PI-AN2MAX,PI-AN2MIN
3129C
3130 CALL PHO_EVENT(-2,P1,P2,WEIGHT,IREJ)
3131 IF(NITER.GT.1) THEN
3132 CALL PHO_PHIST(-2,WEIGHT)
3133 CALL PHO_LHIST(-2,WEIGHT)
3134 ELSE
3135 WRITE(LO,'(1X,A,I4)') 'PHO_GPHERA:NO OUTPUT OF HISTOGRAMS',NITER
3136 ENDIF
3137
3138 END
3139
3140*$ CREATE PHO_GGEPEM.FOR
3141*COPY PHO_GGEPEM
3142CDECK ID>, PHO_GGEPEM
3143 SUBROUTINE PHO_GGEPEM(NEVENT,EE1,EE2)
3144C**********************************************************************
3145C
3146C interface to call PHOJET (variable energy run) for
3147C gamma-gamma collisions on e+e- collider
3148C
3149C fully differential equivalent (improved) photon approximation
3150C to get photon flux
3151C
3152C input: EE1 LAB system energy of electron/positron 1
3153C EE2 LAB system energy of electron/positron 2
3154C NEVENT >0 number of events to generate
3155C -1 initialization
3156C -2 final call (cross section calculation)
3157C from /LEPCUT/:
3158C YMIN1 lower limit of Y1
3159C (energy fraction taken by photon from electron)
3160C YMAX1 upper limit of Y1
3161C Q2MIN1 lower limit of photon virtuality
3162C Q2MAX1 upper limit of photon virtuality
3163C THMIN1 lower limit of scattered electron
3164C THMAX1 upper limit of scattered electron
3165C YMIN2 lower limit of Y2
3166C (energy fraction taken by photon from electron)
3167C YMAX2 upper limit of Y2
3168C Q2MIN2 lower limit of photon virtuality
3169C Q2MAX2 upper limit of photon virtuality
3170C THMIN2 lower limit of scattered electron
3171C THMAX2 upper limit of scattered electron
3172C
3173C output: after final call with NEVENT=-2
3174C EE1 e+ e- cross section (mb)
3175C EE2 gamma-gamma cross section (mb)
3176C
3177C**********************************************************************
3178
3179 IMPLICIT NONE
3180
3181 SAVE
3182
3183 DOUBLE PRECISION EE1,EE2
3184 INTEGER NEVENT
3185
3186C input/output channels
3187 INTEGER LI,LO
3188 COMMON /POINOU/ LI,LO
3189C event debugging information
3190 INTEGER NMAXD
3191 PARAMETER (NMAXD=100)
3192 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
3193 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3194 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
3195 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3196C model switches and parameters
3197 CHARACTER*8 MDLNA
3198 INTEGER ISWMDL,IPAMDL
3199 DOUBLE PRECISION PARMDL
3200 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
3201C some constants
3202 DOUBLE PRECISION PI,PI2,PI4,GEV2MB,Q_ch,Q_ch2,Q_ch4
3203 COMMON /POCONS/ PI,PI2,PI4,GEV2MB,
3204 & Q_ch(-6:6),Q_ch2(-6:6),Q_ch4(-6:6)
3205C photon flux kinematics and cuts
3206 DOUBLE PRECISION ECMIN,ECMAX,EEMIN1,EEMIN2,
3207 & YMIN1,YMAX1,YMIN2,YMAX2,
3208 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
3209 & THMIN1,THMAX1,THMIN2,THMAX2
3210 INTEGER ITAG1,ITAG2
3211 COMMON /POFCUT/ ECMIN,ECMAX,EEMIN1,EEMIN2,
3212 & YMIN1,YMAX1,YMIN2,YMAX2,
3213 & Q2MIN1,Q2MAX1,Q2MIN2,Q2MAX2,
3214 & THMIN1,THMAX1,THMIN2,THMAX2,
3215 & ITAG1,ITAG2
3216C gamma-lepton or gamma-hadron vertex information
3217 INTEGER IGHEL,IDPSRC,IDBSRC
3218 DOUBLE PRECISION PINI,PFIN,PGAM,GYY,GQ2,GGECM,GAIMP,PFTHE,PFPHI,
3219 & RADSRC,AMSRC,GAMSRC
3220 COMMON /POFSRC/ PINI(5,2),PFIN(5,2),PGAM(5,2),IGHEL(2),
3221 & GYY(2),GQ2(2),GGECM,GAIMP(2),PFTHE(2),PFPHI(2),
3222 & IDPSRC(2),IDBSRC(2),RADSRC(2),AMSRC(2),GAMSRC(2)
3223C nucleon-nucleus / nucleus-nucleus interface to DPMJET
3224 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
3225 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
3226 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
3227 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
3228C event weights and generated cross section
3229 INTEGER IPOWGC,ISWCUT,IVWGHT
3230 DOUBLE PRECISION SIGGEN,HSWGHT,HSWCUT,EVWGHT
3231 COMMON /POWGHT/ SIGGEN(4),HSWGHT(0:10),HSWCUT(20),EVWGHT(0:10),
3232 & IPOWGC(0:10),ISWCUT(20),IVWGHT(0:10)
3233
3234C external functions
3235 DOUBLE PRECISION DT_RNDM
3236
3237C local variables
3238 DOUBLE PRECISION AN1MAX,AN1MIN,AN2MAX,AN2MIN,AY1,AY2,AYS1,AYS2,
3239 & COF1,COF2,CPFTHE,DAY1,DAY2,DELLY,DITRY,DITRW,
3240 & ECFRAC,ECMAX2,ECMIN2,EGAM,ELEM,ELEM2,FFL,FFT,FLUXL,FLUXT,
3241 & FLXAPP,FLXQPM,GGECM2,P1,P2,PP,PT,PT2,Q21AV2,Q21AVE,Q21MAX,
3242 & Q21MIN,Q22AV2,Q22AVE,Q22MAX,Q22MIN,Q2LOG1,Q2LOG2,Q2LOW1,
3243 & Q2LOW2,Q2P1,Q2P2,SIF1,SIF2,SIGCUR,SIGMAX,THMAC1,
3244 & THMAC2,THMIC1,THMIC2,WEIGHT,WG,WGFX,WGH,WGHAPP,WGHL,WGHQPM,
3245 & WGMAX,WGY,X1DEL,X1MAX,X1MIN,X2DEL,X2MAX,X2MIN,Y1,Y2,YEFF1,YEFF2,
3246 & YMI,YY1MAX,YY1MIN,YY2MAX,YY2MIN
3247
3248 INTEGER I,IHEAC1,IHEAC2,IHETRY,IREJ,ITRW_low,ITRW_high,ITRY_low,
3249 & ITRY_high,K,Max_tab,NITER,ITG1,ITG2
3250
3251 DIMENSION P1(4),P2(4),IHETRY(4),IHEAC1(4),IHEAC2(4)
3252 integer ipho_pdg2id
3253
3254C initialization of event generation
3255
3256 if(NEVENT.eq.-1) then
3257
3258 DO 10 I=1,4
3259 IHETRY(I) = 0
3260 IHEAC1(I) = 0
3261 IHEAC2(I) = 0
3262 10 CONTINUE
3263
3264 WRITE(LO,'(//1X,A)') 'PHO_GGEPEM: initialization'
3265
3266C electron data
3267 ELEM = 0.512D-03
3268 ELEM2 = ELEM**2
3269 AMSRC(1) = ELEM
3270 AMSRC(2) = ELEM
3271C lepton numbers
3272 IDPSRC(1) = 11
3273 IDPSRC(2) = -11
3274 IDBSRC(1) = ipho_pdg2id(11)
3275 IDBSRC(2) = ipho_pdg2id(-11)
3276
3277C check/update kinematic limitations
3278
3279 Ymi = min(Ymax1,1.D0-ELEM/EE1)
3280 if(Ymi.lt.Ymax1) then
3281 WRITE(LO,'(/1X,A,2E12.5)')
3282 & 'PHO_GGEPEM: Ymax1 decreased (old/new)',Ymax1,Ymi
3283 Ymax1 = YMI
3284 endif
3285 Ymi = min(Ymax2,1.D0-ELEM/EE2)
3286 if(Ymi.lt.Ymax2) then
3287 WRITE(LO,'(/1X,A,2E12.5)')
3288 & 'PHO_GGEPEM: Ymax2 decreased (old/new)',Ymax2,Ymi
3289 Ymax2 = YMI
3290 endif
3291
3292 YMI = ECMIN**2/(4.D0*EE1*EE2*YMAX2)
3293 IF(YMIN1.LT.YMI) THEN
3294 WRITE(LO,'(/1X,A,2E12.5)')
3295 & 'PHO_GGEPEM: Ymin1 increased (old/new)',YMIN1,YMI
3296 YMIN1 = YMI
3297 ELSE IF(YMIN1.GT.YMI) THEN
3298 WRITE(LO,'(/1X,A,/1X,A,E12.5,A,E12.5)')
3299 & 'PHO_GGEPEM:','ECM-CUT corresponds to YMIN1 of',YMI,
3300 & ' INSTEAD OF',YMIN1
3301 ENDIF
3302 YMI = ECMIN**2/(4.D0*EE1*EE2*YMAX1)
3303 IF(YMIN2.LT.YMI) THEN
3304 WRITE(LO,'(/1X,A,2E12.5)')
3305 & 'PHO_GGEPEM: Ymin2 increased (old/new)',YMIN2,YMI
3306 YMIN2 = YMI
3307 ELSE IF(YMIN2.GT.YMI) THEN
3308 WRITE(LO,'(/1X,A,/1X,A,E12.5,A,E12.5)')
3309 & 'PHO_GGEPEM:','ECM-CUT corresponds to YMIN2 of',YMI,
3310 & ' INSTEAD OF',YMIN2
3311 ENDIF
3312
3313C store COS of angular tagging range
3314 THMIC1 = COS(MAX(0.D0,THMIN1))
3315 THMAC1 = COS(MIN(THMAX1,PI))
3316 THMIC2 = COS(MAX(0.D0,THMIN2))
3317 THMAC2 = COS(MIN(THMAX2,PI))
3318
3319 X1MAX = LOG(YMAX1)
3320 X1MIN = LOG(YMIN1)
3321 X1DEL = X1MAX-X1MIN
3322 X2MAX = LOG(YMAX2)
3323 X2MIN = LOG(YMIN2)
3324 X2DEL = X2MAX-X2MIN
3325
3326C debug: integrated photon flux
3327
3328 if(IDEB(30).ge.1) then
3329 Max_tab = 50
3330 FLUXT = 0.D0
3331 FLUXL = 0.D0
3332 DELLY = LOG(YMAX1/YMIN1)/DBLE(Max_tab-1)
3333 IF(IDEB(30).GE.2) WRITE(LO,'(1X,2A,I5)') 'PHO_GGEPEM: ',
3334 & 'table of photon flux (trans/long side 1)',Max_tab
3335 do I=1,Max_tab
3336 Y1 = EXP(X1MIN+DELLY*DBLE(I-1))
3337 if((1.D0-Y1).gt.1.D-8) then
3338 Q2LOW1 = MAX(Q2MIN1,ELEM2*Y1*Y1/(1.D0-Y1))
3339 else
3340 Q2low1 = 2.D0*Q2max1
3341 endif
3342 if(Q2low1.lt.Q2max1) then
3343 FFT = ((1.D0+(1.D0-Y1)**2)/Y1*LOG(Q2MAX1/Q2LOW1)
3344 & -2.D0*ELEM2*Y1*(1.D0/Q2LOW1-1.D0/Q2MAX1))/(2.D0*PI*137.D0)
3345 FFL = 2.D0*(1.D0-Y1)*LOG(Q2MAX1/Q2LOW1)/(2.D0*PI*137.D0)
3346 else
3347 FFT = 0.D0
3348 FFL = 0.D0
3349 endif
3350 FLUXT = FLUXT + Y1*FFL
3351 FLUXL = FLUXL + Y1*FFT
3352 IF(IDEB(30).GE.2) WRITE(LO,'(5X,1P3E14.4)') Y1,FFT,FFL
3353 enddo
3354 FLUXT = FLUXT*DELLY
3355 FLUXL = FLUXL*DELLY
3356 WRITE(LO,'(1X,2A,1P2E12.4)') 'PHO_GGEPEM: ',
3357 & 'integrated flux (trans/long side 1):',FLUXT,FLUXL
3358 endif
3359
3360C maximum weight
3361
3362 Q2LOW1 = MAX(Q2MIN1,ELEM2*YMIN1**2/(1.D0-YMIN1))
3363 Q2LOW2 = MAX(Q2MIN2,ELEM2*YMIN2**2/(1.D0-YMIN2))
3364 Y1 = YMIN1
3365 Y2 = YMIN2
3366 IF(ISWMDL(10).GE.2) THEN
3367C long. and transversely polarized photons
3368 WGMAX = ((1.D0+(1.D0-Y1)**2+2.D0*(1.D0-Y1))*LOG(Q2MAX1/Q2LOW1)
3369 & -2.D0*ELEM2*Y1*(1.D0/Q2LOW1-1.D0/Q2MAX1)*Y1)
3370 & *((1.D0+(1.D0-Y2)**2+2.D0*(1.D0-Y2))*LOG(Q2MAX2/Q2LOW2)
3371 & -2.D0*ELEM2*Y2*(1.D0/Q2LOW2-1.D0/Q2MAX2)*Y2)
3372 ELSE
3373C transversely polarized photons only
3374 WGMAX = ((1.D0+(1.D0-Y1)**2)*LOG(Q2MAX1/Q2LOW1)
3375 & -2.D0*ELEM2*Y1*(1.D0/Q2LOW1-1.D0/Q2MAX1)*Y1)
3376 & *((1.D0+(1.D0-Y2)**2)*LOG(Q2MAX2/Q2LOW2)
3377 & -2.D0*ELEM2*Y2*(1.D0/Q2LOW2-1.D0/Q2MAX2)*Y2)
3378 ENDIF
3379
3380C initialize gamma-gamma event generator
3381
3382C photon 1
3383 EGAM = YMAX1*EE1
3384 P1(1) = 0.D0
3385 P1(2) = 0.D0
3386 P1(3) = SQRT(EGAM**2-Q2LOW1)
3387 P1(4) = EGAM
3388C photon 2
3389 EGAM = YMAX2*EE2
3390 P2(1) = 0.D0
3391 P2(2) = 0.D0
3392 P2(3) = -SQRT(EGAM**2-Q2LOW2)
3393 P2(4) = EGAM
3394C sum of helicities
3395 IGHEL(1) = -1
3396 IGHEL(2) = -1
3397
3398C set min. energy for interpolation tables
3399 parmdl(19) = min(parmdl(19),ecmin)
3400
3401C initialize event gneration
3402 CALL PHO_SETPAR(1,22,0,0.D0)
3403 CALL PHO_SETPAR(2,22,0,0.D0)
3404 CALL PHO_EVENT(-1,P1,P2,SIGMAX,IREJ)
3405 CALL PHO_PHIST(-1,SIGMAX)
3406 CALL PHO_LHIST(-1,SIGMAX)
3407
3408C generation of events, flux calculation
3409
3410 ECMIN2 = ECMIN**2
3411 ECMAX2 = ECMAX**2
3412 ECFRAC = ECMIN2/(4.D0*EE1*EE2)
3413 AY1 = 0.D0
3414 AY2 = 0.D0
3415 AYS1 = 0.D0
3416 AYS2 = 0.D0
3417 Q21MIN = 1.D30
3418 Q22MIN = 1.D30
3419 Q21MAX = 0.D0
3420 Q22MAX = 0.D0
3421 Q21AVE = 0.D0
3422 Q22AVE = 0.D0
3423 Q21AV2 = 0.D0
3424 Q22AV2 = 0.D0
3425 AN1MIN = 1.D30
3426 AN2MIN = 1.D30
3427 AN1MAX = 0.D0
3428 AN2MAX = 0.D0
3429 YY1MIN = 1.D30
3430 YY2MIN = 1.D30
3431 YY1MAX = 0.D0
3432 YY2MAX = 0.D0
3433 NITER = 0
3434 ITRY_low = 0
3435 ITRY_high = 0
3436 ITRW_low = 0
3437 ITRW_high = 0
3438
3439C generate NEVENT events (might be just 1 per call)
3440
3441 else if(NEVENT.gt.0) then
3442
3443 NITER = NITER+NEVENT
3444
3445 DO 200 I=1,NEVENT
3446
3447C sample y1, y2
3448 150 CONTINUE
3449 ITRY_low = ITRY_low+1
3450 if(ITRY_low.eq.1000000) then
3451 ITRY_low = 0
3452 ITRY_high = ITRY_high+1
3453 endif
3454
3455 175 CONTINUE
3456 ITRW_low = ITRW_low+1
3457 if(ITRW_low.eq.1000000) then
3458 ITRW_low = 0
3459 ITRW_high = ITRW_high+1
3460 endif
3461
3462 Y1 = EXP(X1DEL*DT_RNDM(AY1)+X1MIN)
3463 Y2 = EXP(X2DEL*DT_RNDM(AY2)+X2MIN)
3464 IF(Y1*Y2.LT.ECFRAC) GOTO 175
3465 IF(ISWMDL(10).GE.2) THEN
3466 YEFF1 = 1.D0+(1.D0-Y1)**2+2.D0*(1.D0-Y1)
3467 YEFF2 = 1.D0+(1.D0-Y2)**2+2.D0*(1.D0-Y2)
3468 ELSE
3469 YEFF1 = 1.D0+(1.D0-Y1)**2
3470 YEFF2 = 1.D0+(1.D0-Y2)**2
3471 ENDIF
3472
3473 Q2LOW1 = MAX(Q2MIN1,ELEM2*Y1**2/(1.D0-Y1))
3474 Q2LOW2 = MAX(Q2MIN2,ELEM2*Y2**2/(1.D0-Y2))
3475 Q2LOG1 = LOG(Q2MAX1/Q2LOW1)
3476 Q2LOG2 = LOG(Q2MAX2/Q2LOW2)
3477 WGH = (YEFF1*Q2LOG1
3478 & -2.D0*ELEM2*Y1*(1.D0/Q2LOW1-1.D0/Q2MAX1)*Y1)
3479 & *(YEFF2*Q2LOG2
3480 & -2.D0*ELEM2*Y2*(1.D0/Q2LOW2-1.D0/Q2MAX2)*Y2)
3481 IF(WGMAX.LT.WGH) THEN
3482 WRITE(LO,'(1X,A,4E12.5)')
3483 & 'PHO_GGEPEM: inconsistent weight:',Y1,Y2,WGMAX,WGH
3484 ENDIF
3485 IF(DT_RNDM(AYS1)*WGMAX.GT.WGH) GOTO 175
3486
3487C limit on Ecm_gg (app. cut, precise cut applied later)
3488 GGECM2 = 4.D0*Y1*Y2*EE1*EE2
3489 if(GGECM2.lt.ECMIN2) goto 175
3490
3491C sample Q2
3492 IF(IPAMDL(174).EQ.1) THEN
3493 185 CONTINUE
3494 Q2P1 = Q2LOW1*EXP(Q2LOG1*DT_RNDM(Y1))
3495 WEIGHT = (YEFF1-2.D0*(1.D0-Y1)*Q2LOW1/Q2P1)/YEFF1
3496 IF(WEIGHT.LT.DT_RNDM(Q2P1)) GOTO 185
3497 ELSE
3498 Q2P1 = Q2LOW1
3499 ENDIF
3500
3501 IF(IPAMDL(174).EQ.1) THEN
3502 186 CONTINUE
3503 Q2P2 = Q2LOW2*EXP(Q2LOG2*DT_RNDM(Y2))
3504 WEIGHT = (YEFF2-2.D0*(1.D0-Y2)*Q2LOW2/Q2P2)/YEFF2
3505 IF(WEIGHT.LT.DT_RNDM(Q2P2)) GOTO 186
3506 ELSE
3507 Q2P2 = Q2LOW2
3508 ENDIF
3509
3510 GYY(1) = Y1
3511 GQ2(1) = Q2P1
3512 GYY(2) = Y2
3513 GQ2(2) = Q2P2
3514
3515C incoming electron 1
3516 PINI(1,1) = 0.D0
3517 PINI(2,1) = 0.D0
3518 PINI(3,1) = EE1*(1.D0-0.5D0*ELEM2/EE1**2)
3519 PINI(4,1) = EE1
3520 PINI(5,1) = ELEM
3521C photon 1
3522 PP = (2.D0*EE1**2*Y1+Q2P1)/(2.D0*PINI(3,1))
3523 PT2 = (EE1**2*(Q2P1*(1.D0-Y1)-ELEM2*Y1**2)
3524 & -0.25D0*Q2P1**2-Q2P1*ELEM2)/PINI(3,1)**2
3525 IF(PT2.LT.0.D0) GOTO 175
3526 PT = SQRT(PT2)
3527 CALL PHO_SFECFE(SIF1,COF1)
3528 P1(1) = COF1*PT
3529 P1(2) = SIF1*PT
3530 P1(3) = PP
3531 P1(4) = EE1*Y1
3532C outgoing electron 1
3533 PFIN(1,1) = -P1(1)
3534 PFIN(2,1) = -P1(2)
3535 PFIN(3,1) = PINI(3,1)-P1(3)
3536 PFIN(4,1) = PINI(4,1)-P1(4)
3537 PFIN(5,1) = ELEM
3538C incoming electron 2
3539 PINI(1,2) = 0.D0
3540 PINI(2,2) = 0.D0
3541 PINI(3,2) = -EE2*(1.D0-0.5D0*ELEM2/EE2**2)
3542 PINI(4,2) = EE2
3543 PINI(5,2) = 0.D0
3544C photon 2
3545 PP = (2.D0*EE2**2*Y2+Q2P2)/(2.D0*PINI(3,2))
3546 PT2 = (EE2**2*(Q2P2*(1.D0-Y2)-ELEM2*Y2**2)
3547 & -0.25D0*Q2P2**2-Q2P2*ELEM2)/PINI(3,2)**2
3548 IF(PT2.LT.0.D0) GOTO 175
3549 PT = SQRT(PT2)
3550 CALL PHO_SFECFE(SIF2,COF2)
3551 P2(1) = COF2*PT
3552 P2(2) = SIF2*PT
3553 P2(3) = PP
3554 P2(4) = EE2*Y2
3555C outgoing electron 2
3556 PFIN(1,2) = -P2(1)
3557 PFIN(2,2) = -P2(2)
3558 PFIN(3,2) = PINI(3,2)-P2(3)
3559 PFIN(4,2) = PINI(4,2)-P2(4)
3560 PFIN(5,2) = ELEM
3561
3562C precise ECMS cut
3563
3564 GGECM2 = (P1(4)+P2(4))**2-(P1(1)+P2(1))**2
3565 & -(P1(2)+P2(2))**2-(P1(3)+P2(3))**2
3566 IF((GGECM2.LT.ECMIN2).OR.(GGECM2.GT.ECMAX2)) GOTO 175
3567 GGECM = SQRT(GGECM2)
3568
3569C beam lepton detector acceptance
3570
3571C lepton tagger 1
3572 CPFTHE = PFIN(3,1)/PFIN(4,1)
3573 ITG1 = 0
3574 IF(PFIN(4,1).GE.EEMIN1) THEN
3575 IF((CPFTHE.LE.THMIC1).AND.(CPFTHE.GE.THMAC1)) ITG1 = 1
3576 ENDIF
3577
3578C lepton tagger 2
3579 CPFTHE = PFIN(3,2)/PFIN(4,2)
3580 ITG2 = 0
3581 IF(PFIN(4,2).GE.EEMIN2) THEN
3582 IF((CPFTHE.LE.THMIC2).AND.(CPFTHE.GE.THMAC2)) ITG2 = 1
3583 ENDIF
3584
3585C beam lepton taggers
3586
3587C anti-tag
3588 IF((ITAG1.EQ.-1).AND.(ITG1.NE.0)) GOTO 175
3589 IF((ITAG2.EQ.-1).AND.(ITG2.NE.0)) GOTO 175
3590C tag
3591 IF((ITAG1.EQ.1).AND.(ITG1.EQ.0)) GOTO 175
3592 IF((ITAG2.EQ.1).AND.(ITG2.EQ.0)) GOTO 175
3593C single-tag inclusive
3594 IF((ITAG1.EQ.0).AND.(ITAG2.EQ.0).AND.(ITG1+ITG2.EQ.0))
3595 & GOTO 175
3596C single-tag/anti-tag
3597 IF((ITAG1.EQ.2).AND.(ITAG2.EQ.2).AND.(ITG1+ITG2.NE.1))
3598 & GOTO 175
3599
3600 PGAM(1,1) = P1(1)
3601 PGAM(2,1) = P1(2)
3602 PGAM(3,1) = P1(3)
3603 PGAM(4,1) = P1(4)
3604 PGAM(5,1) = -SQRT(Q2P1)
3605 PGAM(1,2) = P2(1)
3606 PGAM(2,2) = P2(2)
3607 PGAM(3,2) = P2(3)
3608 PGAM(4,2) = P2(4)
3609 PGAM(5,2) = -SQRT(Q2P2)
3610
3611C photon helicities
3612 IF(ISWMDL(10).GE.2) THEN
3613 WGH = YEFF1-2.D0*ELEM2*Y1**2/Q2P1
3614 WGHL = 2.D0*(1-Y1)
3615 IF(DT_RNDM(Y1).GT.WGHL/WGH) THEN
3616 IGHEL(1) = 1
3617 ELSE
3618 IGHEL(1) = 0
3619 ENDIF
3620 WGH = YEFF2-2.D0*ELEM2*Y2**2/Q2P2
3621 WGHL = 2.D0*(1-Y2)
3622 IF(DT_RNDM(Y2).GT.WGHL/WGH) THEN
3623 IGHEL(2) = 1
3624 ELSE
3625 IGHEL(2) = 0
3626 ENDIF
3627 K = 2*IGHEL(1)+IGHEL(2)+1
3628 IHETRY(K) = IHETRY(K)+1
3629 ELSE
3630 IGHEL(1) = -1
3631 IGHEL(2) = -1
3632 ENDIF
3633
3634C user cuts
3635 CALL PHO_PRESEL(5,IREJ)
3636 IF(IREJ.NE.0) GOTO 175
3637
3638 WGFX = 1.D0
3639C reweight according to LO photon emission diagrams (Budnev et al.)
3640 IF(IPAMDL(116).GE.1) THEN
3641 CALL PHO_WGEPEM(FLXAPP,FLXQPM,0)
3642 WGFX = FLXQPM/FLXAPP
3643 if(WGFX.gt.1.D0) then
3644 WRITE(LO,'(1x,a,/,5x,1p,5e11.4)')
3645 & ' PHO_GGEPEM: flux weight > 1 (y1/2,Q21/2,W)',
3646 & Y1,Y2,Q2P1,Q2P2,GGECM
3647 endif
3648 ENDIF
3649
3650C event generation
3651* IVWGHT(1) = 1
3652* EVWGHT(1) = MAX(WGFX,1.D0)
3653 CALL PHO_EVENT(1,P1,P2,SIGCUR,IREJ)
3654 IF(IREJ.NE.0) GOTO 150
3655 IF(ISWMDL(10).GE.2) THEN
3656 K = 2*IGHEL(1)+IGHEL(2)+1
3657 IHEAC1(K) = IHEAC1(K)+1
3658 ENDIF
3659
3660C reweight according to QPM model (e+e- collider only)
3661 IF((KHDIR.GT.0).AND.
3662 & (IPAMDL(116).GE.2).AND.(ISWMDL(10).GE.2)) THEN
3663 CALL PHO_WGEPEM(WGHAPP,WGHQPM,1)
3664 WG = WGHQPM/WGHAPP/MAX(1.D0,WGFX)
3665 IF(DT_RNDM(WG).GT.WG) GOTO 150
3666 ELSE IF(IPAMDL(116).GE.1) THEN
3667 IF(DT_RNDM(WG).GT.WGFX) GOTO 150
3668 ENDIF
3669
3670C polar angle
3671 PFTHE(1) = ACOS(PFIN(3,1)/PFIN(4,1))
3672 PFTHE(2) = ACOS(PFIN(3,2)/PFIN(4,2))
3673C azimuthal angle
3674 PFPHI(1) = ATAN2(COF1,SIF1)
3675 PFPHI(2) = ATAN2(COF2,SIF2)
3676
3677C statistics
3678 AY1 = AY1+Y1
3679 AYS1 = AYS1+Y1*Y1
3680 AY2 = AY2+Y2
3681 AYS2 = AYS2+Y2*Y2
3682 Q21MIN = MIN(Q21MIN,Q2P1)
3683 Q22MIN = MIN(Q22MIN,Q2P2)
3684 Q21MAX = MAX(Q21MAX,Q2P1)
3685 Q22MAX = MAX(Q22MAX,Q2P2)
3686 AN1MIN = MIN(AN1MIN,PFTHE(1))
3687 AN2MIN = MIN(AN2MIN,PFTHE(2))
3688 AN1MAX = MAX(AN1MAX,PFTHE(1))
3689 AN2MAX = MAX(AN2MAX,PFTHE(2))
3690 YY1MIN = MIN(YY1MIN,Y1)
3691 YY2MIN = MIN(YY2MIN,Y2)
3692 YY1MAX = MAX(YY1MAX,Y1)
3693 YY2MAX = MAX(YY2MAX,Y2)
3694 Q21AVE = Q21AVE+Q2P1
3695 Q22AVE = Q22AVE+Q2P2
3696 Q21AV2 = Q21AV2+Q2P1*Q2P1
3697 Q22AV2 = Q22AV2+Q2P2*Q2P2
3698 IF(ISWMDL(10).GE.2) THEN
3699 K = 2*IGHEL(1)+IGHEL(2)+1
3700 IHEAC2(K) = IHEAC2(K)+1
3701 ENDIF
3702
3703C external histograms
3704 CALL PHO_PHIST(1,HSWGHT(0))
3705 CALL PHO_LHIST(1,HSWGHT(0))
3706 200 CONTINUE
3707
3708C final cross section calculation and event generation summary
3709
3710 else if(NEVENT.eq.-2) then
3711
3712* EVWGHT(1) = 1.D0
3713* IVWGHT(1) = 0
3714 DITRY = dble(ITRY_high)*1.D+6+dble(ITRY_low)
3715 DITRW = dble(ITRW_high)*1.D+6+dble(ITRW_low)
3716 WGY = WGMAX*DITRY/DITRW/(137.D0*2.D0*PI)**2
3717 WGY = WGY*LOG(YMAX1/YMIN1)*LOG(YMAX2/YMIN2)
3718 AY1 = AY1/DBLE(NITER)
3719 AYS1 = AYS1/DBLE(NITER)
3720 DAY1 = SQRT((AYS1-AY1**2)/DBLE(NITER))
3721 AY2 = AY2/DBLE(NITER)
3722 AYS2 = AYS2/DBLE(NITER)
3723 DAY2 = SQRT((AYS2-AY2**2)/DBLE(NITER))
3724 Q21AVE = Q21AVE/DBLE(NITER)
3725 Q21AV2 = Q21AV2/DBLE(NITER)
3726 Q21AV2 = SQRT((Q21AV2-Q21AVE**2)/DBLE(NITER))
3727 Q22AVE = Q22AVE/DBLE(NITER)
3728 Q22AV2 = Q22AV2/DBLE(NITER)
3729 Q22AV2 = SQRT((Q22AV2-Q22AVE**2)/DBLE(NITER))
3730 WEIGHT = WGY*SIGMAX*DBLE(NITER)/DITRY
3731 EE1 = WEIGHT
3732 EE2 = SIGMAX*DBLE(NITER)/DITRY
3733
3734C output of statistics, histograms
3735 WRITE(LO,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
3736 & '=========================================================',
3737 & ' ***** simulated cross section: ',WEIGHT,' mb *****',
3738 & '========================================================='
3739 WRITE(LO,'(//1X,A,I10,1p,2e14.6)')
3740 & 'PHO_GGEPEM:summary: NITER,ITRY,ITRW',NITER,DITRY,DITRW
3741 WRITE(LO,'(1X,A,1P2E12.4)') 'effective weight (FLUX,TOTAL)',
3742 & WGY,WEIGHT
3743 WRITE(LO,'(1X,A,1P2E12.4)') 'average Y1,DY1 ',
3744 & AY1,DAY1
3745 WRITE(LO,'(1X,A,1P2E12.4)') 'average Y2,DY2 ',
3746 & AY2,DAY2
3747 WRITE(LO,'(1X,A,1P2E12.4)') 'sampled Y range photon 1 ',
3748 & YY1MIN,YY1MAX
3749 WRITE(LO,'(1X,A,1P2E12.4)') 'sampled Y range photon 2 ',
3750 & YY2MIN,YY2MAX
3751 WRITE(LO,'(1X,A,1P2E12.4)') 'average Q2,DQ2 photon 1 ',
3752 & Q21AVE,Q21AV2
3753 WRITE(LO,'(1X,A,1P2E12.4)') 'sampled Q2 range photon 1 ',
3754 & Q21MIN,Q21MAX
3755 WRITE(LO,'(1X,A,1P2E12.4)') 'average Q2,DQ2 photon 2 ',
3756 & Q22AVE,Q22AV2
3757 WRITE(LO,'(1X,A,1P2E12.4)') 'sampled Q2 range photon 2 ',
3758 & Q22MIN,Q22MAX
3759 WRITE(LO,'(1X,A,1P2E12.4)') 'sampled THETA range electron1',
3760 & AN1MIN,AN1MAX
3761 WRITE(LO,'(1X,A,1P4E12.4)') 'sampled THETA range electron2',
3762 & AN2MIN,AN2MAX,PI-AN2MAX,PI-AN2MIN