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