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