+TITLE. DIPSI 2.4 960430 C +PATCH, DIPCO. +KEEP, DIPCO. IMPLICIT NONE COMMON/PESI/NWT,SUMWT,SUMWT2 COMMON/MYCARD/JEVE,QSQLOW,QSQUP,NUTO,NTPFLAG,EMC, + JMESON,JDKLEP,ICRXGX,IQ2EVO,IFORFA, + FORFAS,KEWGEN,IMASGE,YMIN,YMAX,YGEN, + IQ2BIN,BIPT,PTMIN,PTMAX,FALPHA_S, + XMALIM(2),USRGLU,ETA COMMON/RIEMPI/AMASS,XL,PSI(5),ESE(5),ESP(5),PHILE, + EMU2(5),EMU3(5),VECCM(5),GAMCM(5), + PSICM(5),PHIPRO,PHIMU2,PHIMU3,EBE(5), + EBP(5),QSQ,Y,AVU,PT2CM,SIGMA_GAMMA, + SIGMA_MU,WEIGHT,XBAR,Z,THELE,THEMU2,T, + THEMU3,THEPRO,GAM(5),EGA1(5),EGA2(5), + EPI0(5),SWGA2,PTESP,XBJORK,R,COPHI2, + SETHE2,WTGAMP,SWGA,TT,WSQ,PPSICM,EPSICM, + PGACM,EGACM,TTT,MCCOST,MCPHI,MCPSI, + MCPHIL COMMON/STAL/IJK,ISTIL(11) COMMON/CINEMA/S,PI,AMP,AMP2,AMU,DMU2,AMEL,DMEL2,AMPI, + PSIMAS,PSIMAS2,AMUMA INTEGER NWT, JEVE, NUTO, NTPFLAG, EMC, IJK INTEGER ISTIL, JMESON, JDKLEP, ICRXGX INTEGER IQ2EVO, IFORFA, KEWGEN, IMASGE INTEGER YGEN, IQ2BIN, USRGLU DOUBLE PRECISION ESE, ESP, PHILE, PSI, AMASS, XL DOUBLE PRECISION EMU2, EMU3, PHIPRO, PHIMU2, PHIMU3 DOUBLE PRECISION EBE, EBP, GAM, EGA1, EGA2, EPI0 DOUBLE PRECISION VECCM, GAMCM, PSICM, PGACM REAL XBJORK, PTESP, AVU, Y, QSQ, WEIGHT REAL PT2CM, XBAR, Z, SIGMA_MU, THEPRO REAL THEMU2, R, WTGAMP, SWGA2, THEMU3 REAL THELE, SIGMA_GAMMA, SUMWT, SUMWT2 REAL ETA3, S, PI, AMP, AMP2, AMUMA, AMU REAL AMEL, DMEL2, DMU2, PSIMAS, PSIMAS2 REAL QSQLOW, QSQUP, ETAE, ETAP, T, ETA2 REAL COPHI2, SETHE2, FORFAS, YMIN, YMAX REAL WSQ, AMPI, SWGA, EBEAM, PBEAM, BIPT REAL PTMIN, PTMAX, FALPHA_S, PPSICM, ETA REAL EPSICM, EGACM, TT, TTT REAL MCCOST, MCPHI, MCPSI, XMALIM, MCPHIL +KEEP, DIPSE. COMMON/SEMI/ISEED INTEGER ISEED +PATCH,DIPSI. +DECK,MAIN. C C------------------------------------------------------------------------------ C C ***** Program DIPSI ***** C ***** DIFFRACTIVE VECTOR MESON PRODUCTION GENERATOR ***** C M. Arneodo, L. Lamberti, M. Ryskin, submitted to Comp. Phys. Comm. C C------------------------------------------------------------------------------ C +CDE, DIPCO. COMMON/PAWC/JUNK(200000) REAL JUNK, MEANUT INTEGER KERR CALL HLIMIT(200000) CALL DIPSINI(KERR) DO 111 IJK = 1, NUTO CALL DIPSIGEN CALL TIMTEM(MEANUT) 111 CONTINUE C CALL DIPSOUT(MEANUT) STOP END +PATCH,DIINI. +DECK,DIINI. SUBROUTINE DIPSINI(KERR) C C------------------------------------------------------------------------------ C C *** INITIALISES DIPSI *** C C------------------------------------------------------------------------------ C +CDE, DIPCO. +CDE, DIPSE. CHARACTER DUC*12 INTEGER KERR REAL UT,SPACE LOGICAL MASLOG COMMON/CFREAD/SPACE(1000) MASLOG = .FALSE. KERR = 0 WRITE(6,999) C C *** READ RANDOM NUMBER SEEDS FROM EXTERNAL FILE C READ(9,*,ERR=10) ISEED REWIND(9) 10 IF (ISEED.EQ.0) THEN WRITE(19,*) ' ' WRITE(19,*) 'Dipsi warning: could not find file DIPSI.RND' WRITE(19,*) ' random number seed = 0' WRITE(19,*) ' ' KERR = 1 ENDIF C C *** READ USER CARDS FROM EXTERNAL FILE, SEE APPENDICES B AND C *** C CALL FFINIT(1000) CALL FFSET('SIZE',6) CALL FFSET('LINP',8) C C *** EBEAM: LEPTON BEAM MOMENTUM (NEGATIVE...) *** C CALL FFKEY('EBEAM ',EBEAM,1,'REAL') C C *** PBEAM: PROTON BEAM MOMENTUM *** C CALL FFKEY('PBEAM ',PBEAM,1,'REAL') C C *** NTPFLAG = 1 : FILL N-TUPLES *** C CALL FFKEY('NTPFLA',NTPFLAG,1,'INTE') C C *** EMC KINEMATICS FLAG: EMC = 1 muon beam, otherwise electrons *** C CALL FFKEY('EMC ',EMC,1,'INTE') C C *** NUTO = TOTAL NUMBER OF GENERATED EVENTS *** C CALL FFKEY('NUTO ',NUTO,1,'INTE') CALL FFKEY('JEVE ',JEVE,1,'INTE') C C *** Y GENERATION PARAMETERS *** C CALL FFKEY('YMIN ',YMIN,1,'REAL') CALL FFKEY('YMAX ',YMAX,1,'REAL') CALL FFKEY('YGEN ',YGEN,1,'INTE') C C *** QSQ GENERATION PARAMETERS *** C CALL FFKEY('QSQLOW',QSQLOW,1,'REAL') CALL FFKEY('QSQUP ',QSQUP,1,'REAL') CALL FFKEY('KEWGEN',KEWGEN,1,'INTE') C C *** VECTOR MESON CODE *** C CALL FFKEY('JMESON',JMESON,1,'INTE') C C *** DECAY MODE *** C CALL FFKEY('JDKLEP',JDKLEP,1,'INTE') C C *** MASS GENERATION SPECTRUM AND LIMITS *** C CALL FFKEY('IMASGE',IMASGE,1,'INTE') CALL FFKEY('MASMIN',XMALIM(1),1,'REAL') CALL FFKEY('MASMAX',XMALIM(2),1,'REAL') C C *** GLUON DISTRIBUTION CARD *** C CALL FFKEY('ICRXGX',ICRXGX,1,'INTE') CALL FFKEY('USRGLU',USRGLU,1,'INTE') C C *** GLUON DISTRIBUTION Q2 VALUE *** C CALL FFKEY('IQ2BIN',IQ2BIN,1,'INTE') C C *** GLUON DISTRIBUTION EVOLUTION *** C CALL FFKEY('IQ2EVO',IQ2EVO,1,'INTE') C C *** FORM FACTOR CARD *** C CALL FFKEY('IFORFA',IFORFA,1,'INTE') C C *** EXPONENTIAL FORM FACTOR SLOPE *** C CALL FFKEY('FORFAS',FORFAS,1,'REAL') C C *** PT2 GENERATION SLOPE *** C CALL FFKEY('BIPT ',BIPT,1,'REAL') CALL FFKEY('PTMIN ',PTMIN,1,'REAL') CALL FFKEY('PTMAX ',PTMAX,1,'REAL') C C *** ALPHA STRONG *** C CALL FFKEY('ALPHAS',FALPHA_S,1,'REAL') C C *** ETA (NORMALISATION CORRECTION FOR RELATIVISTIC C *** WAVE FUNCTION) *** C CALL FFKEY('ETA',ETA,1,'REAL') CALL FFGO C C *** SET SOME DEFAULTS IF VALUES READ IN MAKE NO SENSE C IF (EBEAM.GE.0) EBEAM = -27.5 IF (PBEAM.LT.0) PBEAM = 820. IF (NUTO.LE.0) NUTO = 10 IF (YMIN.LE.0.) YMIN = 1.E-3 IF (YMAX.GT.1.) YMAX = 1. IF (YMAX.LE.YMIN) THEN WRITE(6,*) ' ' WRITE(6,*) ' YMIN > YMAX. Set to default values.' WRITE(6,*) ' ' YMIN = 1.E-3 YMAX = 1. ENDIF IF (QSQLOW.EQ.0.) QSQLOW = 1.E-7 IF (QSQUP .EQ.0.) QSQUP = 1.E+2 IF (QSQUP.LE.QSQLOW) THEN WRITE(6,*) ' ' WRITE(6,*) ' QSQLOW > QSQUP. Set to default values.' WRITE(6,*) ' ' QSQLOW = 1.E-7 QSQUP = 1.E+2 ENDIF IF (PTMIN.LT.0.) PTMIN = 0. IF (PTMAX.LE.PTMIN) THEN PTMIN = 0. PTMAX = 1000. ENDIF IF (XMALIM(2).LE.XMALIM(1)) THEN WRITE(6,*) ' ' WRITE(6,*) ' MASMIN > MASMAX. Set to default values.' WRITE(6,*) ' ' MASLOG = .TRUE. ENDIF IF (BIPT.LE.0) BIPT = 3. IF ((IQ2BIN.LE.1) .OR. (IQ2BIN.GT.14)) IQ2BIN = 7 IF (JMESON.EQ.1) THEN PSIMAS = 0.7683 IF (MASLOG) THEN XMALIM(1) = 2.*0.1395675 XMALIM(2) = 1.5 ENDIF ELSE IF (JMESON.EQ.2) THEN PSIMAS = 1.019412 IF (MASLOG) THEN XMALIM(1) = 2.*0.497672 XMALIM(2) = 1.2 ENDIF ELSE IF (JMESON.EQ.3) THEN PSIMAS = 9.46032 ELSE IF (JMESON.EQ.4) THEN PSIMAS = 0.78195 IF (MASLOG) THEN XMALIM(1) = 0.732 XMALIM(2) = 0.832 ENDIF ELSE IF (JMESON.EQ.10) THEN PSIMAS = 3.686 ELSE IF (JMESON.EQ.11) THEN PSIMAS = 1.451 IF (MASLOG) THEN XMALIM(1) = 1.0 XMALIM(2) = 2.4 ENDIF ELSE IF (JMESON.EQ.21) THEN PSIMAS = 1.712 IF (MASLOG) THEN XMALIM(1) = 1.0 XMALIM(2) = 2.4 ENDIF ELSE PSIMAS = 3.0969 ENDIF C C *** ETA: CORRECTION TO NORMALISATION OF CROSS SECTION C *** TO TAKE INTO ACCOUNT RELATIVISTIC VS C *** NON-RELATIVISTIC VECTOR MESON WAVEFUNCTIONS *** C IF (ETA.LE.0.) THEN ETA=1. ENDIF WRITE(6,*) 'DIPSI cards: ' WRITE(6,*) ' ' WRITE(6,*) 'DIPSI-EBEAM : ', EBEAM WRITE(6,*) 'DIPSI-PBEAM : ', PBEAM WRITE(6,*) 'DIPSI-NTPFLAG : ', NTPFLAG WRITE(6,*) 'DIPSI-EMC : ', EMC WRITE(6,*) 'DIPSI-NUTO : ', NUTO WRITE(6,*) 'DIPSI-JEVE : ', JEVE WRITE(6,*) 'DIPSI-YMIN : ', YMIN WRITE(6,*) 'DIPSI-YMAX : ', YMAX WRITE(6,*) 'DIPSI-YGEN : ', YGEN WRITE(6,*) 'DIPSI-QSQLOW : ', QSQLOW WRITE(6,*) 'DIPSI-QSQUP : ', QSQUP WRITE(6,*) 'DIPSI-KEWGEN : ', KEWGEN WRITE(6,*) 'DIPSI-JMESON : ', JMESON WRITE(6,*) 'DIPSI-JDKLEP : ', JDKLEP WRITE(6,*) 'DIPSI-IMASGE : ', IMASGE WRITE(6,*) 'DIPSI-MASMIN : ', XMALIM(1) WRITE(6,*) 'DIPSI-MASMAX : ', XMALIM(2) WRITE(6,*) 'DIPSI-USRGLU : ', USRGLU WRITE(6,*) 'DIPSI-ICRXGX : ', ICRXGX WRITE(6,*) 'DIPSI-IQ2BIN : ', IQ2BIN WRITE(6,*) 'DIPSI-IQ2EVO : ', IQ2EVO WRITE(6,*) 'DIPSI-IFORFA : ', IFORFA WRITE(6,*) 'DIPSI-FORFAS : ', FORFAS WRITE(6,*) 'DIPSI-BIPT : ', BIPT WRITE(6,*) 'DIPSI-PTMIN : ', PTMIN WRITE(6,*) 'DIPSI-PTMAX : ', PTMAX WRITE(6,*) 'DIPSI-ALPHAS : ', FALPHA_S WRITE(6,*) 'DIPSI-ETA : ', ETA WRITE(6,*) ' ' WRITE(18,*) 'DIPSI cards: ' WRITE(18,*) ' ' WRITE(18,*) 'DIPSI-EBEAM : ', EBEAM WRITE(18,*) 'DIPSI-PBEAM : ', PBEAM WRITE(18,*) 'DIPSI-NTPFLAG : ', NTPFLAG WRITE(18,*) 'DIPSI-EMC : ', EMC WRITE(18,*) 'DIPSI-NUTO : ', NUTO WRITE(18,*) 'DIPSI-JEVE : ', JEVE WRITE(18,*) 'DIPSI-YMIN : ', YMIN WRITE(18,*) 'DIPSI-YMAX : ', YMAX WRITE(18,*) 'DIPSI-YGEN : ', YGEN WRITE(18,*) 'DIPSI-QSQLOW : ', QSQLOW WRITE(18,*) 'DIPSI-QSQUP : ', QSQUP WRITE(18,*) 'DIPSI-KEWGEN : ', KEWGEN WRITE(18,*) 'DIPSI-JMESON : ', JMESON WRITE(18,*) 'DIPSI-JDKLEP : ', JDKLEP WRITE(18,*) 'DIPSI-IMASGE : ', IMASGE WRITE(18,*) 'DIPSI-MASMIN : ', XMALIM(1) WRITE(18,*) 'DIPSI-MASMAX : ', XMALIM(2) WRITE(18,*) 'DIPSI-USRGLU : ', USRGLU WRITE(18,*) 'DIPSI-ICRXGX : ', ICRXGX WRITE(18,*) 'DIPSI-IQ2BIN : ', IQ2BIN WRITE(18,*) 'DIPSI-IQ2EVO : ', IQ2EVO WRITE(18,*) 'DIPSI-IFORFA : ', IFORFA WRITE(18,*) 'DIPSI-FORFAS : ', FORFAS WRITE(18,*) 'DIPSI-BIPT : ', BIPT WRITE(18,*) 'DIPSI-PTMIN : ', PTMIN WRITE(18,*) 'DIPSI-PTMAX : ', PTMAX WRITE(18,*) 'DIPSI-ALPHAS : ', FALPHA_S WRITE(18,*) 'DIPSI-ETA : ', ETA CALL LUBOOK() C C *** CONSTANTS *** C PI = 3.141592654 AMPI = 0.1349739 AMP = 0.938272 AMP2 = AMP*AMP AMUMA = 0.105658387 AMEL = 0.000511 DMEL2 = AMEL*AMEL EBP(1) = 0. EBP(2) = 0. EBP(3) = DBLE(PBEAM) EBP(5) = .938272 EBP(4) = SQRT(EBP(3)**2 + EBP(5)**2) EBE(1) = 0. EBE(2) = 0. EBE(3) = DBLE(EBEAM) EBE(5) = DBLE(AMEL) EBE(4) = SQRT(EBE(3)**2 + EBE(5)**2) IF (EMC.EQ.1) THEN EBE(5) = DBLE(AMUMA) EBE(4) = SQRT(EBE(3)**2 + EBE(5)**2) ENDIF S = EBE(5)**2 + AMP2 + 2*(EBE(4)*EBP(4) - EBE(3)*EBP(3)) 999 FORMAT(/,10X, + '****************************************************',/, + 10X,'* *',/, + 10X,'* *',/, + 10X,'* Generator DIPSI 2.4 *',/, + 10X,'* *',/, + 10X,'* *',/, + 10X,'* Elastic vector meson production *',/, + 10X,'* *',/, + 10X,'* *',/, + 10X,'* M. Arneodo, L. Lamberti and M. Ryskin *',/, + 10X,'* *',/, + 10X,'* *',/, + 10X,'* Comments and questions to: LAMBL@VXDESY.DESY.DE *',/, + 10X,'* *',/, + 10X,'* *',/, + 10X,'****************************************************',/) IF (QSQUP.GT.1000.) THEN KERR = 1 ENDIF CALL TIMED(UT) RETURN END +DECK,NTUP. SUBROUTINE NTPINIT(NTPNAM) C C------------------------------------------------------------------------------ C C *** N-TUPLE INITIALISATION *** C C------------------------------------------------------------------------------ C IMPLICIT NONE CHARACTER NTPNAM(70)*8 NTPNAM(1) = 'Q2' NTPNAM(2) = 'Y' NTPNAM(3) = 'NU' NTPNAM(4) = 'PT2CM' NTPNAM(5) = 'WSQ' NTPNAM(6) = 'Z' NTPNAM(7) = 'T' NTPNAM(8) = 'XL' NTPNAM(9) = 'PT' NTPNAM(10) = 'XBAR' NTPNAM(11) = 'Q2BAR' NTPNAM(12) = 'HCOSTH' NTPNAM(13) = 'HPHI' NTPNAM(14) = 'HPHIC' NTPNAM(15) = 'HPSI' NTPNAM(16) = 'WEIGHT' NTPNAM(17) = 'WTGAMP' NTPNAM(18) = 'EBE1' NTPNAM(19) = 'EBE2' NTPNAM(20) = 'EBE3' NTPNAM(21) = 'EBE4' NTPNAM(22) = 'EBP1' NTPNAM(23) = 'EBP2' NTPNAM(24) = 'EBP3' NTPNAM(25) = 'EBP4' NTPNAM(26) = 'ESE1' NTPNAM(27) = 'ESE2' NTPNAM(28) = 'ESE3' NTPNAM(29) = 'ESE4' NTPNAM(30) = 'ESE5' NTPNAM(31) = 'ESP1' NTPNAM(32) = 'ESP2' NTPNAM(33) = 'ESP3' NTPNAM(34) = 'ESP4' NTPNAM(35) = 'ESP5' NTPNAM(36) = 'GAM1' NTPNAM(37) = 'GAM2' NTPNAM(38) = 'GAM3' NTPNAM(39) = 'GAM4' NTPNAM(40) = 'GAM5' NTPNAM(41) = 'VEC1' NTPNAM(42) = 'VEC2' NTPNAM(43) = 'VEC3' NTPNAM(44) = 'VEC4' NTPNAM(45) = 'VEC5' NTPNAM(46) = 'MUP1' NTPNAM(47) = 'MUP2' NTPNAM(48) = 'MUP3' NTPNAM(49) = 'MUP4' NTPNAM(50) = 'MUP5' NTPNAM(51) = 'MUM1' NTPNAM(52) = 'MUM2' NTPNAM(53) = 'MUM3' NTPNAM(54) = 'MUM4' NTPNAM(55) = 'MUM5' NTPNAM(56) = 'PIZ1' NTPNAM(57) = 'PIZ2' NTPNAM(58) = 'PIZ3' NTPNAM(59) = 'PIZ4' NTPNAM(60) = 'PIZ5' NTPNAM(61) = 'PH1' NTPNAM(62) = 'PH2' NTPNAM(63) = 'PH3' NTPNAM(64) = 'PH4' NTPNAM(65) = 'PH5' NTPNAM(66) = 'PHO1' NTPNAM(67) = 'PHO2' NTPNAM(68) = 'PHO3' NTPNAM(69) = 'PHO4' NTPNAM(70) = 'PHO5' RETURN END +DECK,LUBO. SUBROUTINE LUBOOK() C C------------------------------------------------------------------------------ C C *** BOOKING OF HISTOGRAMS AND N-TUPLES *** C C------------------------------------------------------------------------------ C +CDE, DIPCO. CHARACTER NTPNAM(70)*8 REAL HQ2UP INTEGER ISTAT IF (NTPFLAG.EQ.1) THEN CALL NTPINIT(NTPNAM) +SELF, IF=VMS. OPEN(41,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=1024, + STATUS='UNKNOWN') CALL HRFILE(41,'ENTUPLE','N') +SELF, IF=UNIX. CALL HROPEN(41,'ENTUPLE','fort.41','N',1024,ISTAT) +SELF. CALL RZQUOT(99999) CALL HBOOKN(666,'DIFFRACTIVE JPSI EVENTS',70,'ENTUPLE', + 70000,NTPNAM) ENDIF CALL HBOOK1(1001,'Q2 GENERATED',100,-7.,3.,0.) IF (EMC.EQ.1) THEN CALL HBOOK1(1002,'VU GENERATED',100,0.,300.,0.) ELSE CALL HBOOK1(1002,'VU GENERATED',100,0.,50000.,0.) ENDIF CALL HBOOK1(1003,'Y GENERATED',100,0.,1.,0.) CALL HBOOK1(1004,'PT2 GENERATED',100,0.,5.,0.) CALL HBOOK1(1005,'LOG10(WEIGHT)',100,-20.,10.,0.) CALL HBOOK1(1006,'LOG10(SIGMA GAMMA)',100,-20.,10.,0.) CALL HBOOK1(1007,'LOG10(SIGMA MU)',100,-20.,10.,0.) CALL HBOOK1(1008,'XBJ GENERATED',100,0.,1.,0.) C CALL HBOOK1(1101,'SIGMA VS Q2',100,-7.,3.,0.) CALL HIDOPT(1101,'LOGY') IF (EMC.EQ.1) THEN CALL HBOOK1(1102,'SIGMA VS VU',100,0.,300.,0.) ELSE CALL HBOOK1(1102,'SIGMA VS VU',100,0.,50000.,0.) ENDIF CALL HBOOK1(1103,'SIGMA VS Y',100,0.,1.,0.) CALL HBOOK1(1104,'SIGMA VS PT2',100,0.,3.,0.) CALL HIDOPT(1104,'LOGY') CALL HBOOK1(11044,'SIGMA VS PT2, z.gt.0.9',100,0.,5.,0.) CALL HIDOPT(11044,'LOGY') CALL HBOOK1(1105,'SIGMA VS Z',100,0.5,1.,0.) CALL HIDOPT(1105,'LOGY') CALL HBOOK1(11055,'SIGMA VS Z',10,0.5,1.,0.) CALL HIDOPT(11055,'LOGY') CALL HBOOK1(1107,'SIGMA VS X GLUE',100,0.,0.1,0.) CALL HIDOPT(1107,'LOGY') CALL HBOOK1(301,'SIGMA VS P(1) PSI',100,-5.,5.,0.) CALL HBOOK1(302,'SIGMA VS P(2) PSI',100,-5.,5.,0.) IF (EMC.EQ.1) THEN CALL HBOOK1(303,'SIGMA VS P(3) PSI',100,-300.,0.,0.) CALL HBOOK1(304,'SIGMA VS E PSI',100,0.,300.,0.) ELSE CALL HBOOK1(303,'SIGMA VS P(3) PSI',100,-40.,200.,0.) CALL HBOOK1(304,'SIGMA VS E PSI',100,0.,200.,0.) ENDIF CALL HBOOK1(2112,'SIGMA VS P(1) MU 1',100,-10.,10.,0.) CALL HBOOK1(2113,'SIGMA VS P(2) MU 1',100,-10.,10.,0.) CALL HBOOK1(2115,'SIGMA VS P(1) MU 2',100,-10.,10.,0.) CALL HBOOK1(2116,'SIGMA VS P(2) MU 2',100,-10.,10.,0.) IF (EMC.EQ.1) THEN CALL HBOOK1(2114,'SIGMA VS P(3) MU 1',100,-300.,0.,0.) CALL HBOOK1(2108,'SIGMA VS E MU 1',100,0.,300.,0.) CALL HBOOK1(2117,'SIGMA VS P(3) MU 2',100,-300.,0.,0.) CALL HBOOK1(2109,'SIGMA VS E MU 2',100,0.,300.,0.) ELSE CALL HBOOK1(2114,'SIGMA VS P(3) MU 1',100,-40.,200.,0.) CALL HBOOK1(2108,'SIGMA VS E MU 1',100,0.,200.,0.) CALL HBOOK1(2117,'SIGMA VS P(3) MU 2',100,-40.,200.,0.) CALL HBOOK1(2109,'SIGMA VS E MU 2',100,0.,200.,0.) ENDIF C IF (JMESON.EQ.1) THEN CALL HBOOK1(1018,'MASS GENERATED',100,0.,1.500,0.) CALL HBOOK1(1118,'INVARIANT MASS 2 PIONS, not weighted', + 100,0.,1.500,0.) CALL HBOOK1(2118,'INVARIANT MASS 2 PIONS',100,0.,1.500,0.) ELSE IF (JMESON.EQ.2) THEN CALL HBOOK1(1018,'MASS GENERATED',100,1.,1.040,0.) CALL HBOOK1(1118,'INVARIANT MASS 2 KAONS not weighted', + 100,1.,1.040,0.) CALL HBOOK1(2118,'INVARIANT MASS 2 KAONS',100,1.,1.040,0.) ELSE IF (JMESON.EQ.4) THEN CALL HBOOK1(1018,'MASS GENERATED',100,0.70,0.84,0.) CALL HBOOK1(1118,'INVARIANT MASS 3 PIONS not weighted', + 100,1.,1.040,0.) CALL HBOOK1(2118,'INVARIANT MASS 3 PIONS',100,0.70,0.84,0.) ELSE IF (JMESON.EQ.3) THEN CALL HBOOK1(1118,'INVARIANT MASS 2 LEPTONS not weighted', + 100,9.459,9.461,0.) CALL HBOOK1(2118,'INVARIANT MASS 2 LEPTONS',100,9.459,9.461,0.) ELSE IF (JMESON.EQ.10) THEN CALL HBOOK1(1118,'INVARIANT MASS 2 LEPTONS not weighted', + 100,3.685,3.687,0.) CALL HBOOK1(2118,'INVARIANT MASS 2 LEPTONS',100,3.685,3.687,0.) ELSE IF (JMESON.EQ.11) THEN CALL HBOOK1(1118,'INVARIANT MASS 2 LEPTONS not weighted', + 100,1.0,2.4,0.) CALL HBOOK1(2118,'INVARIANT MASS 2 LEPTONS',100,1.0,2.4,0.) ELSE IF (JMESON.EQ.21) THEN CALL HBOOK1(1118,'INVARIANT MASS 2 LEPTONS not weighted', + 100,1.0,2.4,0.) CALL HBOOK1(2118,'INVARIANT MASS 2 LEPTONS',100,1.0,2.4,0.) ELSE CALL HBOOK1(1118,'INVARIANT MASS 2 LEPTONS not weighted', + 100,3.096,3.098,0.) CALL HBOOK1(2118,'INVARIANT MASS 2 LEPTONS',100,3.096,3.098,0.) ENDIF C CALL HBOOK1(2125,'SIGMA VS SCATT E MOMENTUM 1',100,-10.,10.,0.) CALL HBOOK1(2126,'SIGMA VS SCATT E MOMENTUM 2',100,-10.,10.,0.) CALL HBOOK1(2129,'SIGMA VS SCATT P MOMENTUM 1',100,-10.,10.,0.) CALL HBOOK1(2130,'SIGMA VS SCATT P MOMENTUM 2',100,-10.,10.,0.) C IF (EMC.EQ.1) THEN CALL HBOOK1(2127,'SIGMA VS SCATT E MOMENTUM 3',100,-300.,0.,0.) CALL HBOOK1(2128,'SIGMA VS SCATT E MOMENTUM 4',100,0.,300.,0.) CALL HBOOK1(2131,'SIGMA VS SCATT P MOMENTUM 3',100,-40.,0.,0.) CALL HBOOK1(2132,'SIGMA VS SCATT P MOMENTUM 4',100,0.,40.,0.) ELSE CALL HBOOK1(2127,'SIGMA VS SCATT E MOMENTUM 3',100,-30.,5.,0.) CALL HBOOK1(2128,'SIGMA VS SCATT E MOMENTUM 4',100,0.,30.,0.) CALL HBOOK1(2131,'SIGMA VS SCATT P MOMENTUM 3',100,600.,850.,0.) CALL HBOOK1(2132,'SIGMA VS SCATT P MOMENTUM 4',100,600.,850.,0.) ENDIF C CALL HBOOK1(99,'SIGMA VS PTESP',100,0.,10.0,0.) CALL HBOOK1(100,'SIGMA VS XL',100,0.9,1.0,0.) IF (EMC.EQ.1) THEN CALL HBOOK1(401,'SIGMA VS ETA E',100,-9.,2.,0.) CALL HBOOK1(402,'SIGMA VS ETA P',100,-9.,2.,0.) CALL HBOOK1(403,'SIGMA VS ETA MU2',100,-9.,2.,0.) CALL HBOOK1(404,'SIGMA VS ETA MU3',100,-9.,2.,0.) ELSE CALL HBOOK1(401,'SIGMA VS ETA E',100,-12.,10.,0.) CALL HBOOK1(402,'SIGMA VS ETA P',100,-12.,10.,0.) CALL HBOOK1(403,'SIGMA VS ETA MU2',100,-12.,10.,0.) CALL HBOOK1(404,'SIGMA VS ETA MU3',100,-12.,10.,0.) ENDIF CALL HBOOK1(101,'SIGMA VS THETA E ',90,90.,180.,0.) CALL HBOOK1(102,'SIGMA VS THETA MU1',180,0.,180.,0.) CALL HBOOK1(103,'SIGMA VS THETA MU2',180,0.,180.,0.) CALL HBOOK1(104,'SIGMA VS THETA P ',100,0.,1.,0.) C CALL HBOOK1(105,'SIGMA VS PHI E ',180,0.,360.,0.) CALL HBOOK1(106,'SIGMA VS PHI MU1',180,0.,360.,0.) CALL HBOOK1(107,'SIGMA VS PHI MU2',180,0.,360.,0.) CALL HBOOK1(108,'SIGMA VS PHI P ',180,0.,360.,0.) C CALL HBOOK2(2101,'LOG10(Q2) VS VU WEIGHTED', + 100,0.,50000.,100,-7.,2.0,0.) CALL HBOOK2(201,'THETA MU 1 VS THETA MU 2 WEIGHTED', + 100,0.,180.,100,0.,180.,0.) RETURN END +PATCH,DIOUT. +DECK,DIOUT. SUBROUTINE DIPSOUT(MEANUT) C C------------------------------------------------------------------------------ C C *** TERMINATES DIPSI *** C C------------------------------------------------------------------------------ C +CDE, DIPCO. +CDE, DIPSE. INTEGER IJL, ID2, IT2, ID3, IT3, IT4 REAL CROSS, SIGSIG, CROSGA, SIGGAM, MEANUT CHARACTER VECTOR*7, DECPA*16, YSPEC*4, QSPEC*6, GLUO*10 WRITE(9,*) ISEED CROSS = SUMWT/NUTO SIGSIG = SQRT(SUMWT2)/NUTO CROSGA = SWGA/NUTO SIGGAM = SQRT(SWGA2)/NUTO WRITE(6,1003) NUTO WRITE(6,1008) NWT, SUMWT, SUMWT2 WRITE(18,1006) WRITE(18,*) '****************************************************' WRITE(18,1006) WRITE(18,1008) NWT, SUMWT, SUMWT2 WRITE(18,1006) WRITE(18,*) '****************************************************' WRITE(18,1006) WRITE(18,1004) CROSS, SIGSIG WRITE(18,1005) CROSGA, SIGGAM WRITE(18,1006) WRITE(18,*) '****************************************************' WRITE(18,1006) WRITE(18,1007) MEANUT WRITE(18,1006) WRITE(18,*) '****************************************************' WRITE(18,1006) C DO 1001 IJL = 1, 11 WRITE(18,1002) IJL, ISTIL(IJL) 1001 CONTINUE CALL DATIME(ID2,IT2) ID3 = 19000000 + ID2 IT3 = INT(IT2/100) IT4 = MOD(IT2,100) IF (JMESON.EQ.0) VECTOR = ' JPSI' IF (JMESON.EQ.10) VECTOR = 'PSIprim' IF (JMESON.EQ.1) VECTOR = ' RHO' IF (JMESON.EQ.11) VECTOR = 'RHOprim' IF (JMESON.EQ.21) VECTOR = 'RHOseco' IF (JMESON.EQ.2) VECTOR = ' PHI' IF (JMESON.EQ.3) VECTOR = 'UPSILON' IF (JMESON.EQ.4) VECTOR = ' OMEGA' IF (JDKLEP.EQ.0) DECPA = 'mu+ mu- ' IF (JDKLEP.EQ.1) DECPA = 'e+ e- ' IF (JDKLEP.EQ.2) DECPA = 'pi+ pi- ' IF (JDKLEP.EQ.3) DECPA = 'k+ k- ' IF (JDKLEP.EQ.4) DECPA = 'k0s k0l ' IF (JDKLEP.EQ.5) DECPA = 'pi+ pi- pi0 ' IF (JDKLEP.EQ.10) DECPA = 'pi+ pi- mu+ mu- ' IF (JDKLEP.EQ.11) DECPA = 'pi+ pi- e+ e- ' IF (JDKLEP.EQ.12) DECPA = 'pi+ pi- pi+ pi- ' IF (JDKLEP.EQ.15) DECPA = 'pi0 pi0 mu+ mu- ' IF (JDKLEP.EQ.16) DECPA = 'pi0 pi0 e+ e- ' IF (JDKLEP.EQ.17) DECPA = 'pi0 pi0 pi+ pi- ' IF (YGEN.EQ.0) YSPEC = '1/Y ' IF (YGEN.EQ.1) YSPEC = 'flat' IF (KEWGEN.EQ.0) QSPEC = '1/Q**2' IF (KEWGEN.EQ.1) QSPEC = '1/Q**4' IF (KEWGEN.EQ.2) QSPEC = 'flat' IF (ICRXGX.EQ.0) GLUO = '3(1-x)**5 ' IF (ICRXGX.GT.0) GLUO = 'NMC distr.' IF (USRGLU.EQ.1) GLUO = 'User def.' WRITE(6,1000) VECTOR, DECPA, BIPT, YSPEC, QSPEC, GLUO, ETA, ID3, + IT3, IT4, MEANUT, NUTO, SUMWT, SUMWT2, QSQLOW, QSQUP, + YMIN, YMAX, PTMIN, PTMAX, CROSS, SIGSIG, CROSGA, SIGGAM WRITE(18,1000) VECTOR, DECPA, BIPT, YSPEC, QSPEC, GLUO, ETA, ID3, + IT3, IT4, MEANUT, NUTO, SUMWT, SUMWT2, QSQLOW, QSQUP, + YMIN, YMAX, PTMIN, PTMAX, CROSS, SIGSIG, CROSGA, SIGGAM 1000 FORMAT(/, + 10X,'****************************************************',/, + 10X,'* *',/, + 10X,'* *',/, + 10X,'* DIPSI generator output *',/, + 10X,'* *',/, + 10X,'* *',/, + 10X,'****************************************************',/, + 10X,'* *',/, + 10X,'* Selected process: ',A7,' --> ',A16,'*', /, + 10X,'* *',/, + 10X,'* Generation parameters *',/, + 10X,'* t slope: ',F6.2,' *',/, + 10X,'* y spectrum: ',A4,' *',/, + 10X,'* Q2 spectrum: ',A6,' *',/, + 10X,'* *',/, + 10X,'* Input gluon distribution *',/, + 10X,'* ',A10,' *', /, + 10X,'* *',/, + 10X,'* eta= ',F6.2,' *',/, + 10X,'* *',/, + 10X,'* Sample produced: ',I8,'; time: ',I2,':',I2,9X,' *',/, + 10X,'* *',/, + 10X,'****************************************************',/, + 10X,'* *',/, + 10X,'* Time to generate one event is ',F10.6, ' seconds *',/, + 10X,'* *',/, + 10X,'****************************************************',/, + 10X,'* *',/, + 10X,'* Generated events and weights *',/, + 10X,'* *',/, + 10X,'* Generated events: ',I14,X,'*', /, + 10X,'* Sum of generated weights: ',E14.8,X,'*', /, + 10X,'* Sum of generated weights squared: ',E14.8,X,'*', /, + 10X,'* *',/, + 10X,'****************************************************',/, + 10X,'* *',/, + 10X,'* Computed cross sections *',/, + 10X,'* *',/, + 10X,'* integrated in the range *',/, + 10X,'* *',/, + 10X,'*',12X,E9.3,' < Q2 < ',E9.3,12X,'*', /, + 10X,'*',15X,F6.4,' < y < ',F6.4,16X,'*', /, + 10X,'*',11X,F9.4,' < pt2 < ',F9.4,12X,'*', /, + 10X,'* *',/, + 10X,'* sigma_ep = (',E13.6,' +/- ',E13.6,') nb *', /, + 10X,'* sigma_gp = (',E13.6,' +/- ',E13.6,') nb *', /, + 10X,'* *',/, + 10X,'****************************************************',/, + 10X,'* *',/, + 10X,'* Questions to: LAMBL@VXDESY.DESY.DE *',/, + 10X,'* *',/, + 10X,'****************************************************') CALL LUBEND(NTPFLAG) 1007 FORMAT(1X,' TIME TO GENERATE ONE EVENT IS ',F13.6,' SECONDS') 1003 FORMAT(1X,' GENERATED EVENTS: ',I10) 1002 FORMAT(1X,' ISTIL(',I2,') = ',I10) 1004 FORMAT(1X,' TOTAL EP CROSS SECTION: ',E13.6,' +/- ',E13.6) 1005 FORMAT(1X,' TOTAL GP CROSS SECTION: ',E13.6,' +/- ',E13.6) 1008 FORMAT(1X,' NWT = ',I10,' SUMWT = ',E15.5,' SUMWT2 = ',E15.5) 1006 FORMAT(/) RETURN END +DECK,BEND. SUBROUTINE LUBEND(NTPFLAG) C C------------------------------------------------------------------------------ C C *** CLOSING OF HBOOK UNFORMATTED AND ASCII FILES *** C C------------------------------------------------------------------------------ C IMPLICIT NONE INTEGER NTPFLAG, ICY, ICW, ISTAT ICY = 1 ICW = 1 IF (NTPFLAG.EQ.1) THEN CALL HCDIR('//ENTUPLE',' ') CALL HROUT(666,ICY,' ') CALL HREND('ENTUPLE') CLOSE(41) ENDIF +SELF, IF=VMS. OPEN(40,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=1024, + STATUS='UNKNOWN') CALL HRFILE(40,'TIPTAP','N') +SELF, IF=UNIX. CALL HROPEN(40,'TIPTAP','fort.40','N',1024,ISTAT) +SELF. CALL HROUT(0,ICW,' ') CALL HREND('TIPTAP') CLOSE(40) CALL HISTDO RETURN END +PATCH,DIPSIGEN. +DECK,DIGEN. SUBROUTINE DIPSIGEN C C------------------------------------------------------------------------------ C C *** GENERATES ELASTIC VECTOR MESON EVENTS ACCORDING TO THE MODEL C *** OF M.RYSKIN (LUND PREPRINT LUTP 92-12, MARCH 1992 AND Z. PHYS. C57 C *** (1993) 89) C C *** BASED ON EMC AND NMC INELASTIC J/PSI GENERATORS BY T. SLOAN, C *** C. MARIOTTI AND M. DE JONG C C------------------------------------------------------------------------------ C +CDE, DIPCO. +CDE, DIPSE. DOUBLE PRECISION PROCM(5), DAMU, BI, CMS(5), DPH, DCT, DCP, DST DOUBLE PRECISION DQ2MIN, DQ2MAX, BETAX, BETAY, BETAZ, PHIS, DSP DOUBLE PRECISION P(5), AMASS2, B(5), GAMMA, PESC, CTH, STH, DTH DOUBLE PRECISION WTE, WTP, WTM, WT3, ETE, ETP, ETM, ET3, A(5) DOUBLE PRECISION EGR(5), EBR(5), PBR(5), SPPR(5), PSPR(5) DOUBLE PRECISION EBPCM(5), ULANGL, DCTPTR, DCTMTR, DPHITR DOUBLE PRECISION DCTPLO, DCTMLO, DPHILO REAL RANN, P2PRO, QSQMIN, QSQMAX, COMET, PFPT2, AWM REAL PFY, PFQSQ, P2MU2, P2MU3, YY, YG, Q2FLUX, PAWW, PBWW REAL ANUM, ADEN, PT2MIN, PT2MAX, Z1, P2ELE, TWIPHI, AWN REAL PHICM, SINTH, THCM, RHOMAS, SIGRHO, RAN1, RAN2, F REAL FMAX, SCALE, SODING, PARSOE(6), FACQSQ, PARBREI(2) REAL BREIT, EPSILON, WIGNER, PARWIG(3), FREQ, AWT, AWL INTEGER JJ, IROTA LOGICAL POL FREQ = 1000. ISTIL(1)=ISTIL(1)+1 IF (INT(ISTIL(1)/FREQ).EQ.(ISTIL(1)/FREQ)) THEN WRITE(6,6) ISTIL(1) ENDIF 6 FORMAT(X,'DIPSI: generating event ',I8,'.') C C *** SET DECAY PARTICLE MASSES FOR TWO BODY DECAYS *** C IF ((JMESON.EQ.1) .OR. (JMESON.EQ.11) .OR. (JMESON.EQ.21)) THEN AMU = 0.1395675 IF (XMALIM(1).LT.(2.*AMU)) XMALIM(1) = 2.*AMU IF (XMALIM(2).LE.(2.*AMU)) XMALIM(2) = 1.500 ELSE IF (JMESON.EQ.2) THEN IF (JDKLEP.EQ.4) THEN AMU = 0.497671 IF (XMALIM(1).LT.(2.*AMU)) XMALIM(1) = 2.*AMU IF (XMALIM(2).LE.(2.*AMU)) XMALIM(2) = 1.200 ELSE IF (JDKLEP.EQ.5) THEN AMU = 0. IF (XMALIM(1).LT.0.420) XMALIM(1) = 0.420 IF (XMALIM(2).LE.0.420) XMALIM(2) = 1.200 ELSE AMU = 0.493646 IF (XMALIM(1).LT.(2.*AMU)) XMALIM(1) = 2.*AMU IF (XMALIM(2).LE.(2.*AMU)) XMALIM(2) = 1.200 ENDIF ELSE IF (JMESON.EQ.4) THEN IF (JDKLEP.EQ.2) THEN AMU = 0.1395675 IF (XMALIM(1).LT.(2.*AMU)) XMALIM(1) = 2.*AMU IF (XMALIM(2).LE.(2.*AMU)) XMALIM(2) = 0.9 ELSE AMU = 0. IF (XMALIM(1).LT.0.420) XMALIM(1) = 0.420 IF (XMALIM(2).LE.0.420) XMALIM(2) = 0.900 ENDIF ELSE IF (JDKLEP.EQ.1) THEN AMU = 0.000511 ELSE IF (JDKLEP.EQ.5) THEN AMU = 0. ELSE AMU = 0.105658387 ENDIF ENDIF C C ---------------------------------------------------- C >>>>>>>> LOOP OVER KINEMATIC VARIABLES: <<<<<<<<<< C ---------------------------------------------------- C 10 CONTINUE C ISTIL(2)=ISTIL(2)+1 NWT=NWT+1 IF (YGEN.EQ.1) THEN C *** GENERATE Y FLAT: *** RANN = RAN(ISEED) Y = YMIN + RANN*(YMAX-YMIN) PFY = 1./(YMAX - YMIN) ELSE C *** GENERATE Y AS 1/Y: *** RANN = RAN(ISEED) Y = YMIN*(YMAX/YMIN)**RANN PFY = (1./Y)*(1./ALOG(YMAX/YMIN)) ENDIF C C *** NU *** C AVU = (EBE(4) * (EBP(3) + EBP(4))/AMP) * Y C C *** GENERATE PSIMAS IN CASE OF RHO OR PHI *** C IF (JMESON.EQ.1) THEN PSIMAS = 0.7683 SIGRHO = 0.1491 612 CONTINUE RAN1 = RAN(ISEED) RAN2 = RAN(ISEED) IF (IMASGE.EQ.1) THEN C *** FLAT/GAUSSIAN *** RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 C RHOMAS = PSIMAS + SIN(2*PI*RAN1)*SIGRHO*SQRT(-2*ALOG(RAN2)) IF (RHOMAS.LE.(2.*AMU)) GOTO 612 ELSE IF (IMASGE.EQ.2) THEN C *** SOEDING *** FMAX = 80. RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARSOE(1) = 9.86 PARSOE(2) = PSIMAS PARSOE(3) = SIGRHO PARSOE(4) = 3.47 PARSOE(5) = 0.0 PARSOE(6) = AMU F = SODING(RHOMAS,PARSOE) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 612 ELSE IF (IMASGE.EQ.3) THEN C *** RELATIVISTIC BREIT-WIGNER *** FMAX = 6.80 RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARWIG(1) = PSIMAS PARWIG(2) = SIGRHO PARWIG(3) = AMU F = WIGNER(RHOMAS,PARWIG) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 612 ELSE C *** BREIT-WIGNER *** FMAX = 180. RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARBREI(1) = PSIMAS PARBREI(2) = SIGRHO F = BREIT(RHOMAS,PARBREI) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 612 ENDIF ISTIL(3)=ISTIL(3)+1 PSIMAS = RHOMAS ELSE IF (JMESON.EQ.2) THEN PSIMAS = 1.019412 SIGRHO = 0.00441 613 CONTINUE RAN1 = RAN(ISEED) RAN2 = RAN(ISEED) IF (IMASGE.EQ.1) THEN C *** FLAT/GAUSSIAN *** RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 C RHOMAS = PSIMAS + SIN(2*PI*RAN1)*SIGRHO*SQRT(-2*ALOG(RAN2)) IF (JDKLEP.EQ.5) THEN IF (RHOMAS.LE..420) GOTO 613 ELSE IF (RHOMAS.LE.(2.*AMU)) GOTO 613 ENDIF ELSE IF ((IMASGE.EQ.3) .AND. (JDKLEP.NE.5)) THEN C *** RELATIVISTIC BREIT-WIGNER *** FMAX = 225. RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARWIG(1) = PSIMAS PARWIG(2) = SIGRHO PARWIG(3) = AMU F = WIGNER(RHOMAS,PARWIG) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 613 ELSE C *** BREIT-WIGNER *** FMAX = 200000. RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARBREI(1) = PSIMAS PARBREI(2) = SIGRHO F = BREIT(RHOMAS,PARBREI) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 613 ENDIF ISTIL(3)=ISTIL(3)+1 PSIMAS = RHOMAS ELSE IF (JMESON.EQ.4) THEN PSIMAS = 0.78195 SIGRHO = 0.00843 614 CONTINUE RAN1 = RAN(ISEED) RAN2 = RAN(ISEED) IF (IMASGE.EQ.1) THEN C *** FLAT/GAUSSIAN *** RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 C RHOMAS = PSIMAS + SIN(2*PI*RAN1)*SIGRHO*SQRT(-2*ALOG(RAN2)) IF (JDKLEP.EQ.5) THEN IF (RHOMAS.LE..420) GOTO 614 ELSE IF (RHOMAS.LE.(2.*AMU)) GOTO 614 ENDIF ELSE IF ((IMASGE.EQ.3) .AND. (JDKLEP.NE.5)) THEN C *** RELATIVISTIC BREIT-WIGNER *** FMAX = 120. RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARWIG(1) = PSIMAS PARWIG(2) = SIGRHO PARWIG(3) = AMU F = WIGNER(RHOMAS,PARWIG) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 614 ELSE C *** BREIT-WIGNER *** FMAX = 56500. RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARBREI(1) = PSIMAS PARBREI(2) = SIGRHO F = BREIT(RHOMAS,PARBREI) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 614 ENDIF ISTIL(3)=ISTIL(3)+1 PSIMAS = RHOMAS ELSE IF (JMESON.EQ.11) THEN PSIMAS = 1.451 SIGRHO = 0.310 615 CONTINUE RAN1 = RAN(ISEED) RAN2 = RAN(ISEED) IF (IMASGE.EQ.1) THEN C *** FLAT/GAUSSIAN *** RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 C RHOMAS = PSIMAS + SIN(2*PI*RAN1)*SIGRHO*SQRT(-2*ALOG(RAN2)) IF (JDKLEP.EQ.12) THEN IF (RHOMAS.LE.1.00) GOTO 615 ELSE IF (RHOMAS.LE.(2.*AMU)) GOTO 615 ENDIF ELSE IF ((IMASGE.EQ.3) .AND. (JDKLEP.NE.5)) THEN C *** RELATIVISTIC BREIT-WIGNER *** FMAX = 3.25 RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARWIG(1) = PSIMAS PARWIG(2) = SIGRHO PARWIG(3) = AMU F = WIGNER(RHOMAS,PARWIG) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 615 ELSE C *** BREIT-WIGNER *** FMAX = 41.7 RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARBREI(1) = PSIMAS PARBREI(2) = SIGRHO F = BREIT(RHOMAS,PARBREI) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 615 ENDIF ISTIL(3)=ISTIL(3)+1 PSIMAS = RHOMAS ELSE IF (JMESON.EQ.21) THEN PSIMAS = 1.712 SIGRHO = 0.213 616 CONTINUE RAN1 = RAN(ISEED) RAN2 = RAN(ISEED) IF (IMASGE.EQ.1) THEN C *** FLAT/GAUSSIAN *** RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 C RHOMAS = PSIMAS + SIN(2*PI*RAN1)*SIGRHO*SQRT(-2*ALOG(RAN2)) IF (JDKLEP.EQ.12) THEN IF (RHOMAS.LE.1.0) GOTO 616 ELSE IF (RHOMAS.LE.(2.*AMU)) GOTO 616 ENDIF ELSE IF ((IMASGE.EQ.3) .AND. (JDKLEP.NE.5)) THEN C *** RELATIVISTIC BREIT-WIGNER *** FMAX = 4.72 RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARWIG(1) = PSIMAS PARWIG(2) = SIGRHO PARWIG(3) = AMU F = WIGNER(RHOMAS,PARWIG) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 616 ELSE C *** BREIT-WIGNER *** FMAX = 90. RHOMAS = XMALIM(1) + (XMALIM(2)-XMALIM(1))*RAN1 PARBREI(1) = PSIMAS PARBREI(2) = SIGRHO F = BREIT(RHOMAS,PARBREI) SCALE = FMAX*RAN2 IF (SCALE.GT.F) GOTO 616 ENDIF ISTIL(3)=ISTIL(3)+1 PSIMAS = RHOMAS ENDIF PSIMAS2 = PSIMAS*PSIMAS DQ2MIN = EBE(5)**2*Y**2/(1-Y) Q2FLUX = SNGL(DQ2MIN) DQ2MAX = AMP * (2.0D0*AVU +AMP) - (PSIMAS+AMP)**2 IF ( DQ2MIN .LT. QSQLOW*1.0D0 ) DQ2MIN = QSQLOW*1.0D0 IF ( DQ2MAX .GT. QSQUP *1.0D0 ) DQ2MAX = QSQUP *1.0D0 C QSQMIN = SNGL(DQ2MIN) QSQMAX = SNGL(DQ2MAX) C IF (QSQMAX.LE.QSQMIN) GOTO 10 ISTIL(4)=ISTIL(4)+1 C C *** QSQ GENERATION *** C IF (KEWGEN.EQ.2) THEN C *** GENERATE QSQ FLAT *** QSQ = QSQMIN + (QSQMAX-QSQMIN)*RAN(ISEED) PFQSQ = 1./(QSQMAX-QSQMIN) ELSE IF (KEWGEN.EQ.1) THEN C *** GENERATE QSQ LIKE 1/QSQ**2 *** FACQSQ = 1./QSQMIN - 1./QSQMAX QSQ = 1./(1./QSQMIN - FACQSQ*RAN(ISEED)) PFQSQ = 1./(FACQSQ*QSQ*QSQ) ELSE C *** GENERATE QSQ LIKE 1/QSQ *** QSQ = QSQMIN*(QSQMAX/QSQMIN)**RAN(ISEED) PFQSQ = (1./QSQ)*(1./ALOG(QSQMAX/QSQMIN)) ENDIF C C *** BJORKEN_X CUT: *** C XBJORK = QSQ / (2.0*AVU*AMP) IF (XBJORK.GT.1.0) GOTO 10 ISTIL(5)=ISTIL(5)+1 C C *** MASS HADRONIC FINAL STATE NOT BELOW J/PSI MASS: *** C WSQ = (2.0*AMP*AVU)+(AMP*AMP)-QSQ IF (WSQ.LT.((PSIMAS+AMP)**2)) GOTO 10 ISTIL(6)=ISTIL(6)+1 C C *** SCATTERED ELECTRON 5-VECTOR. *** C ESE(4) = QSQ/(4*EBE(4)) + EBE(4)*(1-Y) ESE(5) = EBE(5) C CTH = - ( 1.0 - (QSQ/(2*EBE(4)*ESE(4))) ) IF (ABS(CTH).GT.1.0) GOTO 10 ISTIL(7)=ISTIL(7)+1 STH = SQRT(1.0-(CTH*CTH)) C IF ((ESE(4)*ESE(4))-ESE(5)*ESE(5).LT.0.) GOTO 10 ISTIL(8)=ISTIL(8)+1 PESC = SQRT((ESE(4)*ESE(4))-ESE(5)*ESE(5)) C C *** AZIMUTHAL DISTRIBUTION OF SCATTERED ELECTRON: *** C RANN = RAN(ISEED) PHIS = 2.0*PI*RANN ESE(1) = PESC*STH*COS(PHIS) ESE(2) = PESC*STH*SIN(PHIS) ESE(3) = PESC*CTH IF ((-2*ESE(5)**2 + 2*EBE(4)*ESE(4) - 2*ABS(PESC*EBE(3))).GT.QSQ) + GOTO 10 ISTIL(9)=ISTIL(9)+1 C C *** VIRTUAL PHOTON 5 VECTOR *** C CALL SUB5(EBE,ESE,GAM) C C *** BOOST TO PROTON REST FRAME *** C BETAX = EBP(1)/EBP(4) BETAY = EBP(2)/EBP(4) BETAZ = EBP(3)/EBP(4) GAMMA = EBP(4)/EBP(5) DO JJ = 1, 5 A(JJ) = EBE(JJ) B(JJ) = EBP(JJ) P(JJ) = GAM(JJ) ENDDO CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,A) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,B) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,P) DO JJ = 1, 5 EBR(JJ) = A(JJ) PBR(JJ) = B(JJ) EGR(JJ) = P(JJ) ENDDO DO JJ = 1, 4 CMS(JJ) = EBP(JJ) + GAM(JJ) ENDDO CMS(5) = DBLE(SQRT(WSQ)) BETAX = CMS(1)/CMS(4) BETAY = CMS(2)/CMS(4) BETAZ = CMS(3)/CMS(4) GAMMA = CMS(4)/CMS(5) DO JJ = 1, 5 B(JJ) = GAM(JJ) P(JJ) = EBP(JJ) ENDDO CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,B) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,P) DO JJ = 1, 5 GAMCM(JJ) = B(JJ) EBPCM(JJ) = P(JJ) ENDDO C C *** ROTATE TO A SYSTEM WHERE THE PHOTON IS ALIGNED WITH THE Z AXIS *** C DPH = ULANGL(EBPCM(1),EBPCM(2)) CALL RUOTA(0.D0,-DPH,GAMCM) CALL RUOTA(0.D0,-DPH,EBPCM) DTH = ULANGL(EBPCM(3),EBPCM(1)) CALL RUOTA(-DTH,0.D0,GAMCM) CALL RUOTA(-DTH,0.D0,EBPCM) C C *** WORK OUT MOMENTUM OF J/PSI IN CM *** C *** USE E_J/PSI + E_PROTON = W C *** N.B. P_J/PSI = -P_PROTON IN CM C ANUM = 4.*PSIMAS2*AMP2-(PSIMAS2+AMP2-WSQ)**2 ADEN = -4.*WSQ PPSICM = SQRT(ANUM/ADEN) EPSICM = SQRT(PPSICM**2+PSIMAS2) PGACM = SQRT((WSQ-AMP2+QSQ)**2+4.*QSQ*AMP2)/(2.*SQRT(WSQ)) C PT2MIN = PTMIN IF (PTMAX.GT.PPSICM) THEN PT2MAX = PPSICM**2 ELSE PT2MAX = PTMAX ENDIF BI = DBLE(BIPT) C C *** GENERATE PT2CM AS EXP(-BI*PT2CM) *** C COMET = DEXP(-BI*PT2MIN) - DEXP(-BI*PT2MAX) RANN = RAN(ISEED) Z1 = RANN*COMET PT2CM = -(1.0/BI)*DLOG(DEXP(-BI*PT2MIN)-Z1) PFPT2 = BI*EXP(-BI*PT2CM)/COMET C C *** GENERATE PHI OF J/PSI AROUND VIRTUAL PHOTON DIRECTION, FLAT IN CM *** C RANN = RAN(ISEED) PHICM = RANN*2.0*PI C C *** DETERMINE THE COMPONENTS OF THE J/PSI 4-MOMENTUM IN CM *** C SINTH = SQRT(PT2CM)/PPSICM IF(SINTH.GT.1.0) GO TO 10 THCM = ASIN(SINTH) PSICM(1) = SQRT(PT2CM)*COS(PHICM) PSICM(2) = SQRT(PT2CM)*SIN(PHICM) PSICM(3) = -PPSICM*COS(THCM) PSICM(4) = SQRT(PPSICM**2+PSIMAS2) PSICM(5) = PSIMAS C C *** COMPONENTS OF THE PROTON 4-MOMENTUM IN CM *** C PROCM(1) = -PSICM(1) PROCM(2) = -PSICM(2) PROCM(3) = -PSICM(3) PROCM(4) = SQRT(PPSICM**2+AMP2) PROCM(5) = DBLE(AMP) CALL RUOTA(DTH,0.D0,PROCM) CALL RUOTA(0.D0,DPH,PROCM) CALL RUOTA(DTH,0.D0,PSICM) CALL RUOTA(0.D0,DPH,PSICM) CALL RUOTA(DTH,0.D0,EBPCM) CALL RUOTA(0.D0,DPH,EBPCM) C C *** BOOST FROM CM TO LAB FRAME *** C DO 41 JJ = 1, 5 P(JJ) = PSICM(JJ) B(JJ) = PROCM(JJ) 41 CONTINUE BETAX = CMS(1)/CMS(4) BETAY = CMS(2)/CMS(4) BETAZ = CMS(3)/CMS(4) GAMMA = CMS(4)/CMS(5) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,B) DO 42 JJ = 1, 5 PSI(JJ) = P(JJ) ESP(JJ) = B(JJ) 42 CONTINUE P2PRO = SNGL(DSQRT(ESP(1)**2 + ESP(2)**2 + ESP(3)**2)) IF (EBP(3).NE.0.) XL = DBLE(P2PRO)/EBP(3) PTESP = SQRT(ESP(1)**2+ESP(2)**2) TTT = -2. * AMP2 + 2.*EBP(4)*ESP(4) - 2.*EBP(3)*ESP(3) TT = QSQ - PSIMAS2 + 2.*GAM(4)*PSI(4) - 2.*GAM(3)*PSI(3) + - 2.*GAM(2)*PSI(2) - 2.*GAM(1)*PSI(1) C C *** COMPUTE CROSS SECTION *** C CALL JPRYSKIN() IF (XBAR.GT.1.) GOTO 10 ISTIL(10)=ISTIL(10)+1 IF (SIGMA_MU.LE.0.) GOTO 10 ISTIL(11)=ISTIL(11)+1 WEIGHT = SIGMA_MU/(PFQSQ*PFY*PFPT2) WTGAMP = SIGMA_GAMMA/(PFPT2) C C *** Z = FRACTION OF VU CARRIED BY J/PSI (IN PROTON REST...) *** C Z = (PSI(4)*EBP(4)-PSI(3)*EBP(3))/(GAM(4)*EBP(4)-GAM(3)*EBP(3)) C C *** TWIPHI ANGLE BETWEEN ELECTRON SCATTERING PLANE AND GAMMA *** C *** POLARISATION PLANE *** C 11 RANN = RAN(ISEED) PAWW = 1./(1.+((2.-Y)/Y)**2) PBWW = ((2.-Y)/Y)**2/(1.+((2.-Y)/Y)**2) IF (RANN.LT.PAWW) THEN TWIPHI = PI/2. POL = .FALSE. ELSE TWIPHI = 0. POL = .TRUE. ENDIF C C *** BEGIN OF DECAY SECTION *** C EPSILON = (1-Y)/((1-Y+Y**2/2.) - (1-Y)*Q2FLUX/QSQ) IF ((JDKLEP.EQ.5) .AND. (JMESON.NE.1)) THEN CALL TREDK(PSI,EMU2,EMU3,EPI0,EGA1,EGA2,JMESON,JDKLEP) CALL HELOMEGA() SETHE2 = 1. - MCCOST**2 COPHI2 = COS(MCPSI+TWIPHI)*COS(MCPSI+TWIPHI) AWT = 3. * SETHE2*COPHI2 AWL = 3. * (1 - SETHE2) YY = (AWL*R*(1.-Y) + AWT*((1.-Y+Y**2/2.) - (1.-Y)*Q2FLUX/QSQ)) YG = (AWL*EPSILON*R + AWT) ELSE IF ((JDKLEP.GE.10) .AND. (JMESON.GE.10)) THEN CALL TREDK(PSI,EMU2,EMU3,EPI0,EGA1,EGA2,JMESON,JDKLEP) CALL HELOMEGA() AWT = 1. AWL = 1. YY = (AWL*R*(1.-Y) + AWT*((1.-Y+Y**2/2.) - (1.-Y)*Q2FLUX/QSQ)) YG = (AWL*EPSILON*R + AWT) ELSE DMU2 = AMU*AMU DAMU = DBLE(AMU) CALL JDK(DAMU,PSI,EMU2,EMU3) CALL DPHELI() SETHE2 = 1. - MCCOST**2 COPHI2 = COS(MCPSI+TWIPHI)*COS(MCPSI+TWIPHI) C C *** READJUST WEIGHTS IN ORDER TO TAKE INTO ACCOUNT CROSS SECTION *** C *** AND DECAY ANGULAR DISTRIBUTION *** C IF ((JMESON.EQ.0) .OR. (JMESON.EQ.3) .OR. (JMESON.EQ.10)) + THEN AWT = 3./2. * (1. - SETHE2*COPHI2) AWL = 3./2. * SETHE2 YY = (AWL*R*(1.-Y) + AWT*((1.-Y+Y**2/2.)-(1.-Y)*Q2FLUX/QSQ)) YG = (AWL*EPSILON*R + AWT) ELSE IF ((JMESON.EQ.1) .OR. (JMESON.EQ.2) .OR. (JMESON.EQ.4) + .OR. (JMESON.EQ.11) .OR. (JMESON.EQ.21)) THEN AWT = 3. * SETHE2*COPHI2 AWL = 3. * (1 - SETHE2) YY = (AWL*R*(1.-Y) + AWT*((1.-Y+Y**2/2.)-(1.-Y)*Q2FLUX/QSQ)) YG = (AWL*EPSILON*R + AWT) ENDIF ENDIF WEIGHT = WEIGHT * YY WTGAMP = WTGAMP * YG C C *** END OF DECAY SECTION *** C *** C *** WRITE DOWN SOME NUMBERS FOR CHECKING *** C IF (JDKLEP.LT.5) THEN AMASS2 = 2.0*DBLE(DMU2) + + 2.0*(EMU2(4)*EMU3(4)-EMU2(1)*EMU3(1) - + EMU2(2)*EMU3(2)-EMU2(3)*EMU3(3)) ELSE DAMU = 0.1395675 AMASS2 = 2.*DAMU**2 + EPI0(5)**2 + + 2.0*(EMU2(4)*EMU3(4)-EMU2(1)*EMU3(1) - + EMU2(2)*EMU3(2)-EMU2(3)*EMU3(3)) + + 2.0*(EMU2(4)*EPI0(4)-EMU2(1)*EPI0(1) - + EMU2(2)*EPI0(2)-EMU2(3)*EPI0(3)) + + 2.0*(EPI0(4)*EMU3(4)-EPI0(1)*EMU3(1) - + EPI0(2)*EMU3(2)-EPI0(3)*EMU3(3)) ENDIF IF (AMASS2.LE.0.) GOTO 10 AMASS = SQRT(AMASS2) C C *** FINAL STATE PARTICLE RAPIDITY *** C WTP = ESP(3) WTM = EMU2(3) WT3 = EMU3(3) WTE = ESE(3) ETP = ESP(4) ETM = EMU2(4) ET3 = EMU3(4) ETE = ESE(4) IF (((ETE+WTE).LE.0.) .OR. ((ETP+WTP).LE.0.) .OR. + ((ETM+WTM).LE.0.) .OR. ((ET3+WT3).LE.0.) ) THEN WRITE(19,*) ' ' WRITE(19,*) + 'Dipsi warning: Could not compute rapidities for event:', IJK WRITE(19,*) IJK, WTE, ETE WRITE(19,*) IJK, WTP, ETP WRITE(19,*) IJK, WTM, ETM WRITE(19,*) IJK, WT3, ET3 WRITE(19,*) ' ' GOTO 10 ENDIF C C *** FINAL STATE PARTICLE ANGLES *** C P2MU2 = SQRT(EMU2(1)**2 + EMU2(2)**2 + EMU2(3)**2) P2MU3 = SQRT(EMU3(1)**2 + EMU3(2)**2 + EMU3(3)**2) P2ELE = SQRT(ESE(1)**2 + ESE(2)**2 + ESE(3)**2) THEPRO = (ACOS(SNGL(WTP)/P2PRO)) * (180./PI) THEMU2 = (ACOS(SNGL(WTM)/P2MU2)) * (180./PI) THEMU3 = (ACOS(SNGL(WT3)/P2MU3)) * (180./PI) THELE = (ACOS(SNGL(WTE)/P2ELE)) * (180./PI) CALL DEPHI(ESP,PHIPRO) CALL DEPHI(EMU2,PHIMU2) CALL DEPHI(EMU3,PHIMU3) CALL DEPHI(ESE,PHILE) PHIPRO = PHIPRO * (180./PI) PHIMU2 = PHIMU2 * (180./PI) PHIMU3 = PHIMU3 * (180./PI) PHILE = PHILE * (180./PI) IF (IJK.EQ.JEVE) THEN WRITE(17,12) WRITE(17,13) IJK, QSQ, Y, AVU, XBJORK, WSQ WRITE(17,14) WRITE(17,13) IJK, AMASS, (1.-Z), XL, PTESP, TT WRITE(17,18) WRITE(17,16) WRITE(17,17) EBE(1), EBE(2), EBE(3), EBE(4), EBE(5) WRITE(17,17) EBP(1), EBP(2), EBP(3), EBP(4), EBP(5) WRITE(17,17) ESE(1), ESE(2), ESE(3), ESE(4), ESE(5) WRITE(17,17) ESP(1), ESP(2), ESP(3), ESP(4), ESP(5) WRITE(17,17) GAM(1), GAM(2), GAM(3), GAM(4), GAM(5) WRITE(17,17) PSI(1), PSI(2), PSI(3), PSI(4), PSI(5) WRITE(17,17) EMU2(1), EMU2(2), EMU2(3), EMU2(4), EMU2(5) WRITE(17,17) EMU3(1), EMU3(2), EMU3(3), EMU3(4), EMU3(5) IF ((JDKLEP.EQ.5) .OR. (JDKLEP.GE.10)) THEN WRITE(17,17) EPI0(1), EPI0(2), EPI0(3), EPI0(4), EPI0(5) WRITE(17,17) EGA1(1), EGA1(2), EGA1(3), EGA1(4), EGA1(5) WRITE(17,17) EGA2(1), EGA2(2), EGA2(3), EGA2(4), EGA2(5) ENDIF WRITE(17,18) WRITE(17,15) WEIGHT, WTGAMP ENDIF C C *** SUM UP THE WEIGHTS *** C SUMWT = SUMWT + WEIGHT SWGA = SWGA + WTGAMP IF (WEIGHT.GE.1E-16) SUMWT2 = SUMWT2 + WEIGHT**2 IF (WTGAMP.GE.1E-16) SWGA2 = SWGA2 + WTGAMP**2 CALL LUFIL() 12 FORMAT(5X,'Event',7X,'Q2',11X,'y',12X,'Nu',11X,'Xbj',10X,'Wsq') 13 FORMAT(2X,I8,5(2X,E11.5)) 14 FORMAT(17X,'Mass',9X,'1-z',10X,'xl',11X,'pt',11X,'t') 15 FORMAT(3X,'ep weight: ',E10.5,'; gp weight: ',E10.5,'.') 16 FORMAT(6X,'Px',13X,'Py',13X,'Pz',13X,'E',14X,'m') 17 FORMAT(5F15.9) 18 FORMAT(/) RETURN END +DECK,WIGNER. FUNCTION WIGNER(RHOMAS,PARWIG) C C------------------------------------------------------------------------ C C *** RELATIVISTIC BREIT-WIGNER MASS DISTRIBUTION *** C C------------------------------------------------------------------------ C IMPLICIT NONE REAL RHOMAS, PARWIG(3), WIGNER, AMPI, QQ, Q0, WIDT, WID0 REAL VECMAS INTEGER ISECOND ISECOND = 0 VECMAS = PARWIG(1) AMPI = PARWIG(3) WID0 = PARWIG(2) Q0 = SQRT(VECMAS**2 - 4. * AMPI**2) QQ = SQRT(RHOMAS**2 - 4. * AMPI**2) WIDT = WID0 * (QQ/Q0)**3 * (VECMAS/RHOMAS) WIGNER = (RHOMAS * VECMAS * WIDT) / / ((RHOMAS-VECMAS)**2 + VECMAS**2 * WIDT**2) ISECOND = 1 RETURN END +DECK,BREIT. FUNCTION BREIT(RHOMAS,PARBREI) C C------------------------------------------------------------------------ C C *** NON-RELATIVISTIC BREIT-WIGNER MASS DISTRIBUTION *** C C------------------------------------------------------------------------ C IMPLICIT NONE REAL RHOMAS, PARBREI(2), BREIT INTEGER ISECOND ISECOND = 0 BREIT = 1./((RHOMAS-PARBREI(1))**2 + PARBREI(2)**2/4.) ISECOND = 1 RETURN END +DECK,SOED. FUNCTION SODING(RHOMAS,PARSOE) C C------------------------------------------------------------------------ C C *** MASS DISTRIBUTION A LA SOEDING *** C C------------------------------------------------------------------------ C IMPLICIT NONE REAL RHOMAS, PAR(5), PIMA, Q0, QSQ, RAT, G, D2M, SODING REAL PARSOE(6), THRES INTEGER IFIRST PAR(1) = PARSOE(1) PAR(2) = PARSOE(2) PAR(3) = PARSOE(3) PAR(4) = PARSOE(4) PAR(5) = PARSOE(5) PIMA = PARSOE(6) THRES = 2.*PIMA IFIRST = 0 IF (RHOMAS.LT.THRES) THEN SODING = 0. ELSE Q0 = SQRT((PAR(2)/2.)**2 - PIMA**2) QSQ = RHOMAS**2/4. - PIMA**2 RAT = SQRT(QSQ)/Q0 G = 2. * PAR(3) * RAT**3 * PAR(2)/RHOMAS D2M = (PAR(2) + RHOMAS)*(PAR(2) - RHOMAS) SODING = PAR(1) * (PAR(2)*RHOMAS*G/(D2M**2+PAR(2)**2*G**2)) + + PAR(4)*D2M/(D2M**2+PAR(2)**2*G**2) + PAR(5) IF (SODING.LT.0.) SODING = 0. ENDIF IFIRST = 1 RETURN END +DECK,JPRYS. SUBROUTINE JPRYSKIN() C C------------------------------------------------------------------------------ C C *** CROSS SECTION CALCULATION (DDD / DT DQ2 DY) C *** FOR THE PROCESS E(MU)+P-->E(MU) P VECTOR MESON C *** CF. SECTION 2 OF THE PAPER C C------------------------------------------------------------------------------ C +CDE, DIPCO. REAL HBARC2, ALPHA, BRJEE, MASSA2, MASSA3, ALPHA_S, DTDPT2 REAL HBARC2A, BRJEEA, Q2BAR, MASSA, AACUT, FF, LAMBDA REAL CONST, AKIN, XGX, FLUX, BKIN, FKIN, QP C C *** CONSTANTS *** C HBARC2 = 0.389379E+06 ! GeV2 * nanobarn ALPHA = 1/137.036 LAMBDA = 0.2 HBARC2A = 0.38937966 C C *** WIDTH OF VECTOR MESON --> e+ e- (5.36 KEV FOR PSI) *** C IF (JMESON.EQ.1) THEN BRJEE = 6.7700E-06 BRJEEA = 6.7700 MASSA = 0.7683 ELSE IF (JMESON.EQ.2) THEN BRJEE = 1.3700E-06 BRJEEA = 1.3700 MASSA = 1.0194 ELSE IF (JMESON.EQ.3) THEN BRJEE = 1.3400E-06 BRJEEA = 1.3400 MASSA = PSIMAS ELSE IF (JMESON.EQ.4) THEN BRJEE = 0.600E-06 BRJEEA = 0.600 MASSA = 0.78195 ELSE IF (JMESON.EQ.10) THEN BRJEE = 2.1400E-06 BRJEEA = 2.1400 MASSA = 3.0969 ELSE IF (JMESON.EQ.11) THEN BRJEE = 2.1400E-06 ! not available on PR D50 (?) BRJEEA = 2.1400 ! not available on PR D50 (?) MASSA = 0.7683 ELSE IF (JMESON.EQ.21) THEN BRJEE = 2.1400E-06 ! not available on PR D50 (?) BRJEEA = 2.1400 ! not available on PR D50 (?) MASSA = 0.7683 ELSE BRJEE = 5.3600E-06 BRJEEA = 5.3600 MASSA = PSIMAS ENDIF MASSA2 = MASSA*MASSA MASSA3 = MASSA*MASSA2 C C *** KINEMATIC VARIABLES *** C Q2BAR = 0.25*(QSQ+MASSA2+PT2CM) XBAR = 4.0*Q2BAR/WSQ QP = Q2BAR - PT2CM/4. IF (XBAR.GT.1.) RETURN TTT = -2. * AMP2 + 2.*EBP(4)*ESP(4) - 2.*EBP(3)*ESP(3) T = TTT C C *** ALPHA STRONG *** C IF ((FALPHA_S.LE.0.) .OR. (FALPHA_S.GE.1.)) THEN ALPHA_S = (12. * PI)/(25. * ALOG(Q2BAR/LAMBDA**2)) IF (ALPHA_S.GT.0.7) ALPHA_S = 0.7 ELSE ALPHA_S = FALPHA_S ENDIF C C *** THE CROSS SECTION *** C CONST = HBARC2A * ALPHA_S**2 * BRJEEA * MASSA3 * PI**3 + / ( 3. * ALPHA ) AACUT = 0.5 IF (PT2CM.LE.AACUT) THEN FKIN = (4.*QP + AACUT)/(PT2CM + AACUT) ELSE FKIN = ((PT2CM + AACUT)*(PT2CM/2. + 2*QP)**2)/ / ((4*QP + AACUT)*PT2CM**2) ENDIF BKIN = 8.*Q2BAR/AACUT AKIN = ALOG(FKIN)/ALOG(BKIN)/(2.*Q2BAR*(2.*Q2BAR-PT2CM)) C C *** DETERMINE THE GLUON DISTRIBUTION *** C IF (USRGLU.EQ.1) THEN CALL USRGLS(XBAR,Q2BAR,XGX) ELSE CALL GLUONS(ICRXGX,IQ2EVO,IQ2BIN,XBAR,Q2BAR,XGX) ENDIF C C *** FORM FACTOR *** C IF (IFORFA.EQ.1) THEN FF = EXP(-FORFAS * T) ELSE FF = 1.0/(1.0 + T/0.71)**2 ENDIF C C *** T TO PT2 JACOBIAN *** C DTDPT2 = PGACM/SQRT(PPSICM**2-PT2CM) C SIGMA_GAMMA = FF**2 * CONST * XGX**2 * AKIN**2 * ETA**2 SIGMA_GAMMA = SIGMA_GAMMA * DTDPT2 C C *** R = SIGMA_L / SIGMA_T C R = QSQ/MASSA2 C C *** (# OF PHOTONS)/DY DQ2 *** C FLUX = ALPHA/(PI*Y*QSQ) C C *** D SIGMA(MU P-->J/PSI P)/DY DQ2 DT *** C SIGMA_MU = FLUX * SIGMA_GAMMA RETURN END +DECK,USRGLS. SUBROUTINE USRGLS(XBAR,Q2BAR,XGX) C C------------------------------------------------------------------------------ C C *** USER ROUTINE FOR GLUON DISTRIBUTION DETERMINATION AT GIVEN C *** XBAR AND Q2BAR VALUES. C C *** INPUT: XBAR, Q2BAR C *** OUTPUT: XGX C C *** SEE APPENDIX B C------------------------------------------------------------------------------ C IMPLICIT NONE REAL XBAR, Q2BAR, XGX, SCALE INTEGER LLL DOUBLE PRECISION U, D, S, C, B, T, SEA, GLUE, VAL(20) CHARACTER PARM(20)*20 C C *** MY FAVOURITE DISTRIBUTION... *** C XGX = 1. C C *** ... OR WHATEVER YOU LIKE... *** C XGX = 3. * (1-XBAR)**5. C C *** AS AN EXAMPLE USE PDFLIB *** C LLL = LLL + 1 IF (LLL.EQ.1) THEN PARM(1) = 'Nptype' PARM(2) = 'Ngroup' PARM(3) = 'Nset' VAL(1) = 1. ! proton VAL(2) = 3. ! 3=MRS, 5=GRV VAL(3) = 37. ! 31=mrsd-p, 28=mrsd-, 37=mrsa CALL PDFSET(PARM,VAL) ENDIF SCALE = SQRT(Q2BAR) CALL STRUCTF(DBLE(XBAR),DBLE(SCALE),U,D,SEA,S,C,B,T,GLUE) XGX = SNGL(GLUE) RETURN END +DECK,GLUONS. SUBROUTINE GLUONS(ICRXGX,IQ2EVO,IQ2BIN,XBAR,Q2BAR,XGX) C C------------------------------------------------------------------------------ C C *** COMPUTES GLUON DISTRIBUTION FOR GIVEN XBAR AND Q2BAR C *** SEE APPENDIX B C C------------------------------------------------------------------------------ C IMPLICIT NONE INTEGER ICRXGX, IQ2EVO, IQ2BIN REAL XBAR, Q2BAR, XGX REAL FAU, GLUMA(14,7), XGXLO, AXBAR, XZERO, QZERO, AQ2 REAL C1, C2, C3, CC, CET, WWW, FAC, XGXUP, FAD, AXG REAL Q2LOW, Q2UP, CC_O, CET_O, C1_O, C2_O, C3_O, BXG INTEGER IQW, INF QZERO = 2. XZERO = 1.E-6 Q2UP = 0. Q2LOW = 0. CC = 0. CC_O = 0. CET = 0. CET_O = 0. C1 = 0. C1_O = 0. C2 = 0. C2_O = 0. C3 = 0. C3_O = 0. C C *** GLUON DISTRIBUTION FOR 14 VALUES OF Q2 *** C *** THE FIRST INDEX RUNS OVER Q2 VALUES, THE SECOND ONE *** C *** GIVES THE CORRESPONDING COEFFICIENTS *** C GLUMA(1,1) = 1. ! Coeff. # GLUMA(1,2) = 1.25000 ! Q2 value GLUMA(1,3) = 3.88171 ! CC coeff. GLUMA(1,4) = 6.09535 ! CET coeff. GLUMA(1,5) = -1.49500 ! C1 coeff. GLUMA(1,6) = 3.56485 ! C2 coeff. GLUMA(1,7) = -3.28218 ! C3 coeff. GLUMA(2,1) = 2. GLUMA(2,2) = 1.75000 GLUMA(2,3) = 3.89190 GLUMA(2,4) = 6.40695 GLUMA(2,5) = -1.15749 GLUMA(2,6) = 2.84132 GLUMA(2,7) = -2.62854 GLUMA(3,1) = 3. GLUMA(3,2) = 2.50000 GLUMA(3,3) = 3.88090 GLUMA(3,4) = 6.71073 GLUMA(3,5) = -0.84483 GLUMA(3,6) = 2.15951 GLUMA(3,7) = -1.99870 GLUMA(4,1) = 4. GLUMA(4,2) = 3.50000 GLUMA(4,3) = 3.84898 GLUMA(4,4) = 6.96252 GLUMA(4,5) = -0.58866 GLUMA(4,6) = 1.58156 GLUMA(4,7) = -1.44337 GLUMA(5,1) = 5. GLUMA(5,2) = 4.50000 GLUMA(5,3) = 3.83377 GLUMA(5,4) = 7.14367 GLUMA(5,5) = -0.39835 GLUMA(5,6) = 1.11460 GLUMA(5,7) = -0.99992 GLUMA(6,1) = 6. GLUMA(6,2) = 5.50000 GLUMA(6,3) = 3.83322 GLUMA(6,4) = 7.30859 GLUMA(6,5) = -0.37582 GLUMA(6,6) = 1.14878 GLUMA(6,7) = -0.95404 GLUMA(7,1) = 7. GLUMA(7,2) = 7.00000 GLUMA(7,3) = 3.77872 GLUMA(7,4) = 7.42547 GLUMA(7,5) = -0.10379 GLUMA(7,6) = 0.38301 GLUMA(7,7) = -0.28221 GLUMA(8,1) = 8. GLUMA(8,2) = 9.00000 GLUMA(8,3) = 3.80739 GLUMA(8,4) = 7.63949 GLUMA(8,5) = -0.11761 GLUMA(8,6) = 0.41827 GLUMA(8,7) = -0.18860 GLUMA(9,1) = 9. GLUMA(9,2) = 11.50000 GLUMA(9,3) = 3.70374 GLUMA(9,4) = 7.70437 GLUMA(9,5) = 0.21066 GLUMA(9,6) = -0.42927 GLUMA(9,7) = 0.51813 GLUMA(10,1) = 10. GLUMA(10,2) = 15.00000 GLUMA(10,3) = 3.72441 GLUMA(10,4) = 7.90967 GLUMA(10,5) = 0.19062 GLUMA(10,6) = -0.39415 GLUMA(10,7) = 0.61469 GLUMA(11,1) = 11. GLUMA(11,2) = 20.00000 GLUMA(11,3) = 3.60713 GLUMA(11,4) = 7.97344 GLUMA(11,5) = 0.53833 GLUMA(11,6) = -1.30814 GLUMA(11,7) = 1.39402 GLUMA(12,1) = 12. GLUMA(12,2) = 27.00000 GLUMA(12,3) = 3.62692 GLUMA(12,4) = 8.18696 GLUMA(12,5) = 0.51083 GLUMA(12,6) = -1.27250 GLUMA(12,7) = 1.50023 GLUMA(13,1) = 13. GLUMA(13,2) = 36.00000 GLUMA(13,3) = 3.64793 GLUMA(13,4) = 8.39098 GLUMA(13,5) = 0.47587 GLUMA(13,6) = -1.24645 GLUMA(13,7) = 1.61343 GLUMA(14,1) = 14. GLUMA(14,2) = 48.00000 GLUMA(14,3) = 3.47188 GLUMA(14,4) = 8.37209 GLUMA(14,5) = 1.00476 GLUMA(14,6) = -2.64519 GLUMA(14,7) = 2.74528 C C *** SATURATE GLUON DISTRIBUTION *** C IF ((IQ2EVO.GT.1) .AND. (Q2BAR.LE.QZERO)) THEN AQ2 = Q2BAR Q2BAR = QZERO ENDIF IF ((IQ2EVO.EQ.3) .AND. (XBAR.LE.XZERO)) THEN AXBAR = XBAR XBAR = XZERO ENDIF C C *** TRY DIFFERENT XGX STRUCTURE FUNCTIONS *** C IF (ICRXGX.EQ.1) THEN C *** LOWER NMC *** C1 = 0.205 C2 = 0.068 C3 = -0.029 CC = 3.733 CET = 6.391 ELSE IF (ICRXGX.EQ.2) THEN C *** CENTRAL NMC *** IF (IQ2EVO.EQ.0) THEN C1 = GLUMA(IQ2BIN,5) C2 = GLUMA(IQ2BIN,6) C3 = GLUMA(IQ2BIN,7) CC = GLUMA(IQ2BIN,3) CET = GLUMA(IQ2BIN,4) ELSE IF (IQ2EVO.GT.0) THEN C *** Q2 EVOLUTION *** C INF = 2 DO 77 IQW = INF, 14 IF ((IQW-INF) .NE. 0) THEN Q2LOW = GLUMA(IQW-1,2) CC_O = GLUMA(IQW-1,3) CET_O = GLUMA(IQW-1,4) C1_O = GLUMA(IQW-1,5) C2_O = GLUMA(IQW-1,6) C3_O = GLUMA(IQW-1,7) ENDIF Q2UP = GLUMA(IQW,2) CC = GLUMA(IQW,3) CET = GLUMA(IQW,4) C1 = GLUMA(IQW,5) C2 = GLUMA(IQW,6) C3 = GLUMA(IQW,7) IF (Q2UP.GT.Q2BAR) GOTO 78 77 CONTINUE 78 CONTINUE ENDIF ELSE IF (ICRXGX.EQ.3) THEN C *** Upper NMC *** C1 = -1.515 C2 = 2.679 C3 = -1.669 CC = 5.332 CET = 10.610 ELSE C *** DEFAULT *** C1 = 0. C2 = 0. C3 = 0. CC = 3. CET = 5. ENDIF WWW = 0.1 * ALOG(1. + EXP(10. - 100.*XBAR)) IF ((ICRXGX.EQ.2) .AND. (IQ2EVO.GT.0)) THEN FAD = (1. + C1_O*WWW + C2_O*WWW**2 + C3_O*WWW**3) XGXLO = CC_O * ((1.-XBAR)**CET_O) * FAD FAU = (1. + C1*WWW + C2*WWW**2 + C3*WWW**3) XGXUP = CC * ((1.-XBAR)**CET) * FAU C C *** LINEAR INTERPOLATION BETWEEN TWO POINTS *** C *** xgx = axg * q2bar + bxg *** C AXG = (XGXUP - XGXLO)/(Q2UP - Q2LOW) BXG = XGXUP - AXG * Q2UP XGX = AXG * Q2BAR + BXG IF ((IQ2EVO.EQ.3) .AND. (Q2BAR.LE.QZERO)) THEN XGX = XGX * AQ2/Q2BAR ENDIF ELSE FAC = (1. + C1*WWW + C2*WWW**2 + C3*WWW**3) XGX = CC * ((1.-XBAR)**CET) * FAC ENDIF IF ((IQ2EVO.GT.1) .AND. (Q2BAR.LE.QZERO)) THEN Q2BAR = AQ2 ENDIF IF ((IQ2EVO.EQ.3) .AND. (XBAR.LE.XZERO)) THEN XBAR = AXBAR ENDIF RETURN END +DECK,TREDK. SUBROUTINE TREDK(OMEG,P1,P2,P0,P3,P4,JMESON,JDKLEP) C C--------------------------------------------------------------------------- C C *** LET OMEG DECAY INTO PI a AND Pi* AND THEN Pi* INTO Pi b AND Pi c. *** C C *** INPUT: OMEG (PRIMARY MESON) C *** OUTPUT: P1 (pi+), P2 (pi-), P0 (UNSTABLE MESON), P3 AND P4 (PARTICLES C FROM SECONDARY MESON DECAY) C C---------------------------------------------------------------------------- C IMPLICIT NONE +CDE,DIPSE. DOUBLE PRECISION P0(5), PS(5), P1(5), P2(5), BETAX, BETAY, BETAZ DOUBLE PRECISION PI, AMO, DMO, APAI0, DPAI0, APAI, DPAI, GAMMA, RL DOUBLE PRECISION APAIS, DPAIS, THETA_O, PHI_O, THETA_PAI, PHI_PAI DOUBLE PRECISION PMOD2, PMOD0, EPAI0, P3(5), P4(5), E2, MA, AML DOUBLE PRECISION OMEG(5), MB, MC, PEND, B(3,5), P(3,5) DOUBLE PRECISION FOUR(3,3), AMASS, WT1, DDDD, QMAX REAL EZERO, AMPIS, WHICH, RAN, PBR1, PBR2, RM, RW, RANN REAL WMAX, WT2, XW231, XW230, XW131, XW130, XW121, XW120 DOUBLE PRECISION WT12, WT13, WT23 DOUBLE PRECISION XM12, XM13, XM23 INTEGER JJ, II, IJL, JMESON, JDKLEP, ICL, ITO ICL = 0 RM = 0.7683 RW = 0.1491 PI = 3.141592654 APAI = 0.1395675 AMO = OMEG(5) ITO = ITO + 1 10 CONTINUE ICL = ICL + 1 C *** CHOOSE SECONDARY MESON MASS *** IF (JDKLEP.EQ.5) APAI0 = 0.1349739 IF ((JMESON.EQ.10) .AND. (JDKLEP.GE.10)) THEN APAI0 = 3.0969 IF (JDKLEP.GE.15) APAI = 0.1349739 ENDIF IF ((JMESON.EQ.11) .OR. (JMESON.EQ.21)) THEN 612 CONTINUE PBR1 = RAN(ISEED) - 0.5 PBR2 = RAN(ISEED) - 0.5 IF ((PBR1**2+PBR2**2).GT.0.25) GOTO 612 APAI0 = 0.5 * RW * PBR1/PBR2 + RM IF (JDKLEP.EQ.17) APAI = 0.1349739 IF (APAI0.LT.(2.*APAI)) GOTO 612 IF ((APAI0+2.*APAI).GE.OMEG(5)) GOTO 612 ENDIF DPAI0 = APAI0*APAI0 DPAI = APAI*APAI DMO = AMO*AMO WHICH = RAN(ISEED) IF (WHICH.LE.1./3.) THEN MA = APAI0 MB = APAI MC = APAI ELSE IF ((WHICH.GT.1./3.) .AND. (WHICH.LE.2./3.)) THEN MA = APAI MB = APAI0 MC = APAI ELSE IF (WHICH.GT.2./3.) THEN MA = APAI MB = APAI MC = APAI0 ENDIF C *** COMPUTE PI A MOMENTUM END POINT FROM DECAYING PARTICLE MASS *** AMPIS = MB + MC EZERO = (DMO + MA**2 - AMPIS**2)/(2.*AMO) PEND = DSQRT(EZERO**2 - MA**2) C *** GENERATE PI a MOMENTUM *** RANN = RAN(ISEED) PMOD0 = PEND*DBLE(RANN) EPAI0 = DSQRT(PMOD0**2 + MA**2) IF (EPAI0.GE.EZERO) GOTO 10 C *** GENERATE AZIMUTAL AND POLAR DISTRIBUTION *** RANN = RAN(ISEED) THETA_O = DACOS(1. - 2.*DBLE(RANN)) RANN = RAN(ISEED) PHI_O = 2.*PI*DBLE(RANN) C *** VECTOR MESON REST FRAME *** APAIS = DSQRT(DMO + MA**2 - 2.*AMO*EPAI0) DPAIS = APAIS**2 IF (APAIS.LE.(MB+MC)) GOTO 10 B(1,1) = PMOD0*DSIN(THETA_O)*DCOS(PHI_O) B(1,2) = PMOD0*DSIN(THETA_O)*DSIN(PHI_O) B(1,3) = PMOD0*DCOS(THETA_O) B(1,4) = EPAI0 B(1,5) = MA PS(1) = -B(1,1) PS(2) = -B(1,2) PS(3) = -B(1,3) PS(4) = DSQRT(PMOD0**2 + DPAIS) PS(5) = APAIS C *** PI* DECAY INTO PI b AND PI c C *** GENERATE AZIMUTAL AND POLAR DISTRIBUTION *** RANN = RAN(ISEED) THETA_PAI = DACOS(1. - 2.*DBLE(RANN)) RANN = RAN(ISEED) PHI_PAI = 2.*PI*DBLE(RANN) C *** PI* REST FRAME *** E2 = (DPAIS + MB**2 - MC**2)/(2.*APAIS) PMOD2 = DSQRT(E2**2-MB**2) B(2,1) = PMOD2*DSIN(THETA_PAI)*DCOS(PHI_PAI) B(2,2) = PMOD2*DSIN(THETA_PAI)*DSIN(PHI_PAI) B(2,3) = PMOD2*DCOS(THETA_PAI) B(2,4) = E2 B(2,5) = MB B(3,1) = -B(2,1) B(3,2) = -B(2,2) B(3,3) = -B(2,3) B(3,4) = DSQRT(PMOD2**2 + MC**2) B(3,5) = MC C *** BOOST PI b AND PI c TO VECTOR MESON REST FRAME *** BETAX = PS(1)/PS(4) BETAY = PS(2)/PS(4) BETAZ = PS(3)/PS(4) GAMMA = PS(4)/PS(5) IF (WHICH.LE.1./3.) THEN DO IJL = 1, 5 P0(IJL) = B(1,IJL) P1(IJL) = B(2,IJL) P2(IJL) = B(3,IJL) ENDDO CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P1) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P2) ELSE IF ((WHICH.GT.1./3.) .AND. (WHICH.LE.2./3.)) THEN DO IJL = 1, 5 P0(IJL) = B(2,IJL) P1(IJL) = B(3,IJL) P2(IJL) = B(1,IJL) ENDDO CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P0) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P1) ELSE IF (WHICH.GT.2./3.) THEN DO IJL = 1, 5 P0(IJL) = B(3,IJL) P1(IJL) = B(1,IJL) P2(IJL) = B(2,IJL) ENDDO CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P2) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P0) ENDIF DO JJ = 1, 5 P(1,JJ) = P0(JJ) P(2,JJ) = P1(JJ) P(3,JJ) = P2(JJ) ENDDO 11 CONTINUE XM12 = DSQRT(P(1,5)**2 + P(2,5)**2 + 2.*( + P(1,4)*P(2,4) - P(1,1)*P(2,1) - P(1,2)*P(2,2) - + P(1,3)*P(2,3))) XM23 = DSQRT(P(3,5)**2 + P(2,5)**2 + 2.*( + P(3,4)*P(2,4) - P(3,1)*P(2,1) - P(3,2)*P(2,2) - + P(3,3)*P(2,3))) XM13 = DSQRT(P(1,5)**2 + P(3,5)**2 + 2.*( + P(1,4)*P(3,4) - P(1,1)*P(3,1) - P(1,2)*P(3,2) - + P(1,3)*P(3,3))) XW121 = SQRT(XM12**2-(P(1,5)+P(2,5))**2)/2 XW120 = SQRT((AMO**2-(XM12+P(3,5))**2)*(AMO**2-(XM12-P(3,5))**2)) / /(2.*AMO) XW231 = SQRT(XM23**2-(P(3,5)+P(2,5))**2)/2 XW230 = SQRT((AMO**2-(XM23+P(1,5))**2)*(AMO**2-(XM23-P(1,5))**2)) / /(2.*AMO) XW131 = SQRT(XM13**2-(P(1,5)+P(3,5))**2)/2 XW130 = SQRT((AMO**2-(XM13+P(2,5))**2)*(AMO**2-(XM13-P(2,5))**2)) / /(2.*AMO) WT12 = (XW121*XW120)**3 WT23 = (XW231*XW230)**3 WT13 = (XW131*XW130)**3 IF (JMESON.EQ.4) THEN WMAX = 2.E-4 ELSE IF (JMESON.EQ.3) THEN WMAX = 750. ELSE IF (JMESON.EQ.2) THEN WMAX = 2.E-3 ELSE IF ((JMESON.EQ.10) .AND. (JDKLEP.GE.10)) THEN WMAX = 2.25E-3 ELSE IF ((JMESON.EQ.10) .AND. (JDKLEP.EQ.5)) THEN WMAX = 2.25 ELSE IF (JMESON.EQ.11) THEN WMAX = 0.001 ELSE IF (JMESON.EQ.21) THEN WMAX = 0.001 ELSE WMAX = 0.8 ENDIF RANN = RAN(ISEED) IF (WT23.LT.(RANN*WMAX)) GOTO 10 IF ((JMESON.EQ.10) .AND. (JDKLEP.GE.10)) RANN = DBLE(RAN(ISEED)) IF (WT12.LT.(RANN*WMAX)) GOTO 10 IF ((JMESON.EQ.10) .AND. (JDKLEP.GE.10)) RANN = DBLE(RAN(ISEED)) IF (WT13.LT.(RANN*WMAX)) GOTO 10 C *** THEN LET THE SECONDARY MESON DECAY INTO TWO PARTICLES *** IF (JDKLEP.EQ.5) AML = 0.D0 IF (JMESON.EQ.10) THEN IF (JDKLEP.EQ.5) THEN AML = 0.D0 ELSE IF ((JDKLEP.EQ.11) .OR. (JDKLEP.EQ.16)) THEN AML = 0.000511D0 ELSE AML = 0.105658389D0 ENDIF ELSE IF ((JMESON.EQ.11) .OR. (JMESON.EQ.21)) THEN AML = 0.1395675D0 ENDIF CALL JDK(AML,P0,P3,P4) C *** NOW BOOST ALL PARTICLES BACK TO THE LAB FRAME *** BETAX = OMEG(1)/OMEG(4) BETAY = OMEG(2)/OMEG(4) BETAZ = OMEG(3)/OMEG(4) GAMMA = OMEG(4)/OMEG(5) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P0) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P1) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P2) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P3) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P4) RETURN END +DECK,JDK. SUBROUTINE JDK(AMU,PPI,P1,P2) C C------------------------------------------------------------------------------ C C GENERATES TWO-BODY DECAY DECAY OF VECTOR MESON C C------------------------------------------------------------------------------ C IMPLICIT NONE +CDE,DIPSE. DOUBLE PRECISION P1(5), P2(5), PPI(5), BETAX, GAMMA, PHI DOUBLE PRECISION CPJPSI, SPJPSI, COSP, SINP, COST, SINT DOUBLE PRECISION PHDK, W11, W22, W33, E2, E3, P, PLEN, PI DOUBLE PRECISION AMU, THDK, CTJPSI, STJPSI, BETAY, BETAZ REAL RANN INTEGER ITT PI = 3.141592654 C C *** GENERATE AZIMUTHAL AND POLAR DISTRIBUTION DECAY MUON_MOMENTA *** C RANN = RAN(ISEED) CTJPSI = DBLE(2.0*RANN-1.0) THDK = DACOS(CTJPSI) STJPSI = DSIN(THDK) C C *** AZIMUTHAL ANGLE FLAT: *** C RANN = RAN(ISEED) PHI = 2.0*PI*(DBLE(RANN-0.5)) CPJPSI = DCOS(PHI) SPJPSI = DSIN(PHI) PHDK = PHI C C *** REST FRAME OF VECTOR MESON *** C W11 = PPI(5)*PPI(5) W22 = AMU*AMU W33 = AMU*AMU E3 = (W11-W22+W33)/(2.0*PPI(5)) E2 = PPI(5)-E3 P = SQRT(E2*E2-W22) C P1(1) = P*STJPSI*CPJPSI P1(2) = P*STJPSI*SPJPSI P1(3) = P*CTJPSI P1(4) = E2 P1(5) = AMU C P2(1) = -P1(1) P2(2) = -P1(2) P2(3) = -P1(3) P2(4) = E3 P2(5) = AMU C C *** BOOST TO LAB FRAME *** C BETAX = PPI(1)/PPI(4) BETAY = PPI(2)/PPI(4) BETAZ = PPI(3)/PPI(4) GAMMA = PPI(4)/PPI(5) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P1) CALL LORENZ(BETAX,BETAY,BETAZ,GAMMA,P2) RETURN END +DECK,LORE. SUBROUTINE LORENZ(BETAX,BETAY,BETAZ,GAMMA,P) C C------------------------------------------------------------------------------ C C *** LORENZ TRANSFORMS VECTOR P ALONG 3 AXIS. *** C C------------------------------------------------------------------------------ C IMPLICIT NONE DOUBLE PRECISION P(5), AX, AY, BETAX, BETAY, BETAZ, GAMMA AX = BETAX*P(1) + BETAY*P(2) + BETAZ*P(3) AY = GAMMA*(GAMMA/(1.D0 + GAMMA)*AX + P(4)) P(1) = P(1) + (BETAX*AY) P(2) = P(2) + (BETAY*AY) P(3) = P(3) + (BETAZ*AY) P(4) = GAMMA*(P(4) + AX) RETURN END +DECK,555. SUBROUTINE SUB5(P1,P2,P3) C C------------------------------------------------------------------------------ C C *** SUBTRACTS 2 5-VECTORS P1-P2=P3 *** C C------------------------------------------------------------------------------ C IMPLICIT NONE DOUBLE PRECISION P1(5), P2(5), P3(5) REAL AL, VMOD, AMSQ, V(5) INTEGER I DO 1 I = 1, 4 P3(I) = P1(I) - P2(I) V(I) = SNGL(P3(I)) 1 CONTINUE AL=VMOD(V,3) AMSQ=(P3(4)*P3(4))-(AL*AL) P3(5)=SQRT(ABS(AMSQ)) IF(AMSQ.LT.0.0) P3(5)= -P3(5) RETURN END +DECK,PHIPHI. SUBROUTINE DEPHI(P,PHI) C C------------------------------------------------------------------------------ C C *** COMPUTES PHI FROM ATAN2(P(2),P(1)) *** C C------------------------------------------------------------------------------ C IMPLICIT NONE DOUBLE PRECISION P(2), PHI REAL PI PI = 3.141592654 IF ((P(1).EQ.0.).AND.(P(2).EQ.0.)) THEN PHI = 0. ELSE PHI = ATAN2(P(2),P(1)) IF (P(2).GE.0.) PHI = PHI IF (P(2).LT.0.) PHI = PHI + 2.*PI ENDIF RETURN END +DECK,LUFI. SUBROUTINE LUFIL() C C------------------------------------------------------------------------------ C C *** FILLING OF HBOOK HISTOGRAMS AND N-TUPLES *** C C------------------------------------------------------------------------------ C +CDE, DIPCO. REAL NTPVAR(70), Q2BAR Q2BAR = XBAR*WSQ/4. NTPVAR(1) = QSQ NTPVAR(2) = Y NTPVAR(3) = AVU NTPVAR(4) = PT2CM NTPVAR(5) = WSQ NTPVAR(6) = Z NTPVAR(7) = TT NTPVAR(8) = SNGL(XL) NTPVAR(9) = PTESP NTPVAR(10) = XBAR NTPVAR(11) = Q2BAR NTPVAR(12) = MCCOST NTPVAR(13) = MCPHI NTPVAR(14) = MCPHIL NTPVAR(15) = MCPSI NTPVAR(16) = WEIGHT NTPVAR(17) = WTGAMP NTPVAR(18) = SNGL(EBE(1)) NTPVAR(19) = SNGL(EBE(2)) NTPVAR(20) = SNGL(EBE(3)) NTPVAR(21) = SNGL(EBE(4)) NTPVAR(22) = SNGL(EBP(1)) NTPVAR(23) = SNGL(EBP(2)) NTPVAR(24) = SNGL(EBP(3)) NTPVAR(25) = SNGL(EBP(4)) NTPVAR(26) = SNGL(ESE(1)) NTPVAR(27) = SNGL(ESE(2)) NTPVAR(28) = SNGL(ESE(3)) NTPVAR(29) = SNGL(ESE(4)) NTPVAR(30) = SNGL(ESE(5)) NTPVAR(31) = SNGL(ESP(1)) NTPVAR(32) = SNGL(ESP(2)) NTPVAR(33) = SNGL(ESP(3)) NTPVAR(34) = SNGL(ESP(4)) NTPVAR(35) = SNGL(ESP(5)) NTPVAR(36) = SNGL(GAM(1)) NTPVAR(37) = SNGL(GAM(2)) NTPVAR(38) = SNGL(GAM(3)) NTPVAR(39) = SNGL(GAM(4)) NTPVAR(40) = SNGL(GAM(5)) NTPVAR(41) = SNGL(PSI(1)) NTPVAR(42) = SNGL(PSI(2)) NTPVAR(43) = SNGL(PSI(3)) NTPVAR(44) = SNGL(PSI(4)) NTPVAR(45) = SNGL(PSI(5)) NTPVAR(46) = SNGL(EMU2(1)) NTPVAR(47) = SNGL(EMU2(2)) NTPVAR(48) = SNGL(EMU2(3)) NTPVAR(49) = SNGL(EMU2(4)) NTPVAR(50) = SNGL(EMU2(5)) NTPVAR(51) = SNGL(EMU3(1)) NTPVAR(52) = SNGL(EMU3(2)) NTPVAR(53) = SNGL(EMU3(3)) NTPVAR(54) = SNGL(EMU3(4)) NTPVAR(55) = SNGL(EMU3(5)) NTPVAR(56) = SNGL(EPI0(1)) NTPVAR(57) = SNGL(EPI0(2)) NTPVAR(58) = SNGL(EPI0(3)) NTPVAR(59) = SNGL(EPI0(4)) NTPVAR(60) = SNGL(EPI0(5)) NTPVAR(61) = SNGL(EGA1(1)) NTPVAR(62) = SNGL(EGA1(2)) NTPVAR(63) = SNGL(EGA1(3)) NTPVAR(64) = SNGL(EGA1(4)) NTPVAR(65) = SNGL(EGA1(5)) NTPVAR(66) = SNGL(EGA2(1)) NTPVAR(67) = SNGL(EGA2(2)) NTPVAR(68) = SNGL(EGA2(3)) NTPVAR(69) = SNGL(EGA2(4)) NTPVAR(70) = SNGL(EGA2(5)) IF (NTPFLAG.EQ.1) THEN CALL HCDIR('//ENTUPLE',' ') CALL HFN(666,NTPVAR) ENDIF CALL HFILL(1001,alog10(QSQ),0.,1.) CALL HFILL(1002,AVU,0.,1.) CALL HFILL(1003,Y,0.,1.) CALL HFILL(1004,PT2CM,0.,1.) CALL HFILL(1008,XBJORK,0.,1.) CALL HFILL(1118,SNGL(AMASS),0.,1.) IF (WEIGHT.GT.0.) THEN CALL HFILL(1005,ALOG10(WEIGHT),0.,1.) CALL HFILL(1006,ALOG10(SIGMA_GAMMA),0.,1.) CALL HFILL(1007,ALOG10(SIGMA_MU),0.,1.) ENDIF C CALL HFILL(1101,alog10(QSQ),0.,WEIGHT) CALL HFILL(1102,AVU,0.,WEIGHT) CALL HFILL(1103,Y,0.,WEIGHT) CALL HFILL(1104,PT2CM,0.,WEIGHT) IF(Z.GT.0.9) CALL HFILL(11044,PT2CM,0.,WEIGHT) CALL HFILL(1105,Z,0.,WEIGHT) CALL HFILL(11055,Z,0.,WEIGHT) CALL HFILL(1107,XBAR,0.,WEIGHT) C CALL HFILL(2108,SNGL(EMU2(4)),0.,WEIGHT) CALL HFILL(2109,SNGL(EMU3(4)),0.,WEIGHT) CALL HFILL(2112,SNGL(EMU2(1)),0.,WEIGHT) CALL HFILL(2113,SNGL(EMU2(2)),0.,WEIGHT) CALL HFILL(2114,SNGL(EMU2(3)),0.,WEIGHT) CALL HFILL(2115,SNGL(EMU3(1)),0.,WEIGHT) CALL HFILL(2116,SNGL(EMU3(2)),0.,WEIGHT) CALL HFILL(2117,SNGL(EMU3(3)),0.,WEIGHT) IF (PSI(5).LT.3.) CALL HFILL(1018,sngl(PSI(5)),0.,1.) CALL HFILL(2118,SNGL(AMASS),0.,WEIGHT) CALL HFILL(2125,SNGL(ESE(1)),0.,WEIGHT) CALL HFILL(2126,SNGL(ESE(2)),0.,WEIGHT) CALL HFILL(2127,SNGL(ESE(3)),0.,WEIGHT) CALL HFILL(2128,SNGL(ESE(4)),0.,WEIGHT) CALL HFILL(2129,SNGL(ESP(1)),0.,WEIGHT) CALL HFILL(2130,SNGL(ESP(2)),0.,WEIGHT) CALL HFILL(2131,SNGL(ESP(3)),0.,WEIGHT) CALL HFILL(2132,SNGL(ESP(4)),0.,WEIGHT) CALL HFILL(301,SNGL(PSI(1)),0.,WEIGHT) CALL HFILL(302,SNGL(PSI(2)),0.,WEIGHT) CALL HFILL(303,SNGL(PSI(3)),0.,WEIGHT) CALL HFILL(304,SNGL(PSI(4)),0.,WEIGHT) CALL HFILL(401,ETAE,0.,WEIGHT) CALL HFILL(402,ETAP,0.,WEIGHT) CALL HFILL(403,ETA2,0.,WEIGHT) CALL HFILL(404,ETA3,0.,WEIGHT) CALL HFILL(100,SNGL(XL),0.,WEIGHT) CALL HFILL(99,PTESP,0.,WEIGHT) CALL HFILL(101,THELE ,0.,WEIGHT) CALL HFILL(102,THEMU2,0.,WEIGHT) CALL HFILL(103,THEMU3,0.,WEIGHT) CALL HFILL(104,THEPRO,0.,WEIGHT) CALL HFILL(105,SNGL(PHILE),0.,WEIGHT) CALL HFILL(106,SNGL(PHIMU2),0.,WEIGHT) CALL HFILL(107,SNGL(PHIMU3),0.,WEIGHT) CALL HFILL(108,SNGL(PHIPRO),0.,WEIGHT) C CALL HFILL(2101,AVU,ALOG10(QSQ),WEIGHT) CALL HFILL(201,THEMU2,THEMU3,WEIGHT) RETURN END +DECK, TIMT. SUBROUTINE TIMTEM(MEANUT) C C---------------------------------------------------------------------- C C *** TIMING *** C C---------------------------------------------------------------------- C REAL UT, MEANUT, SUMUT INTEGER IT IT = IT + 1 CALL TIMED(UT) SUMUT = SUMUT + UT MEANUT = SUMUT/FLOAT(IT) RETURN END +DECK,ULANGL. DOUBLE PRECISION FUNCTION ULANGL(X,Y) C C----------------------------------------------------------------------- C C *** RECONSTRUCT ANGLE OF VECTOR WRT X AXIS FROM X AND Y COORDINATES *** C C----------------------------------------------------------------------- C IMPLICIT NONE C DOUBLE PRECISION PIGREC, X, Y, R PIGREC = 4.0D0*ATAN(1.0D0) ULANGL = 0.D0 R = SQRT(X**2+Y**2) IF(R.LT.1E-20) RETURN IF(ABS(X)/R.LT.0.8) THEN ULANGL = SIGN(ACOS(X/R),Y) ELSE ULANGL = ASIN(Y/R) IF ((X.LT.0.) .AND. (ULANGL.GE.0.)) THEN ULANGL = PIGREC-ULANGL ELSE IF (X.LT.0.) THEN ULANGL = -PIGREC-ULANGL ENDIF ENDIF RETURN END +DECK, RUOTA. SUBROUTINE RUOTA(DTH,DPH,P) C C------------------------------------------------------------------------- C C *** REFERENCE FRAME ROTATION *** C C------------------------------------------------------------------------- C IMPLICIT NONE DOUBLE PRECISION DROT(3,3), P(3), P8(3), DTH, DPH INTEGER JJ IF ((DTH**2+DPH**2) .GT. (1E-20)) THEN DROT(1,1) = DCOS(DTH)*DCOS(DPH) DROT(1,2) =-DSIN(DPH) DROT(1,3) = DSIN(DTH)*DCOS(DPH) DROT(2,1) = DCOS(DTH)*DSIN(DPH) DROT(2,2) = DCOS(DPH) DROT(2,3) = DSIN(DTH)*DSIN(DPH) DROT(3,1) =-DSIN(DTH) DROT(3,2) = 0.D+0 DROT(3,3) = DCOS(DTH) DO JJ = 1, 3 P8(JJ) = P(JJ) ENDDO DO JJ = 1, 3 P(JJ) = DROT(JJ,1)*P8(1) + DROT(JJ,2)*P8(2) + + DROT(JJ,3)*P8(3) ENDDO ENDIF RETURN END +DECK, HELI. SUBROUTINE DPHELI() C----------------------------------------------------------------------- C C *** TWO BODY DECAY ANGULAR DISTRIBUTION *** C C----------------------------------------------------------------------- +CDE, DIPCO. DOUBLE PRECISION DCTHETA, PGR(4), VSCAT(3), DPHIS, PRHO(4) DOUBLE PRECISION VPROD(3), VSCPR(3), PCM(4), PX(4), PXR(4) DOUBLE PRECISION ZH(3), YH(3), XH(3), PPIR(4), PVGUN(3), VM DOUBLE PRECISION SPHIS, CPHIS, SPHIL, CPHIL, DPHIL, DPSIS DOUBLE PRECISION BETAX, BETAY, BETAZ, GAMMA, VV(4), a(5) DOUBLE PRECISION b(5), c(5) INTEGER JJ C C *** ANGLES ARE DEFINED AS IN THE NMC RHO PAPER, NUCL. PHYS. B429 C *** (1994) 503. THE REFERENCE FOR THIS DEFINITION IS GIVEN IN C *** P. JOOS ET AL, NUCL. PHYS. B113 (1976) 53. *** C DO JJ = 1, 4 prho(JJ) = EMU2(JJ) + EMU3(JJ) PCM(JJ) = GAM(JJ) + EBP(JJ) PX(JJ) = PCM(JJ) - PRHO(JJ) ENDDO C C *** BOOST TO PROTON REST FRAME *** C BETAX = EBP(1)/EBP(4) BETAY = EBP(2)/EBP(4) BETAZ = EBP(3)/EBP(4) GAMMA = EBP(4)/EBP(5) DO JJ = 1, 4 A(JJ) = PRHO(JJ) B(JJ) = GAM(JJ) C(JJ) = EBE(JJ) ENDDO A(5) = DBLE(PSIMAS) B(5) = GAM(5) C(5) = EBE(5) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,A) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,B) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,C) C C *** BOOST SCATTERED PROTON TO VECTOR MESON REST FRAME *** C BETAX = PRHO(1)/PRHO(4) BETAY = PRHO(2)/PRHO(4) BETAZ = PRHO(3)/PRHO(4) GAMMA = PRHO(4)/PSIMAS DO JJ = 1, 4 VV(JJ) = PX(JJ) ENDDO CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,VV) DO JJ = 1, 4 PXR(JJ) = VV(JJ) ENDDO C C *** DEFINE PHOTON FOUR MOMENTUM *** C DO JJ = 1, 4 VV(JJ) = GAM(JJ) ENDDO CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,VV) DO JJ = 1, 4 PGR(JJ) = VV(JJ) ENDDO VM = DSQRT(PGR(1)**2+PGR(2)**2+PGR(3)**2) DO JJ = 1, 3 YH(JJ) = PGR(JJ)/VM ENDDO C C *** DEFINE QUANTISATION AXIS (OPPOSITE TO PROTON DIRECTION IN VM REST C *** FRAME) *** C DO JJ = 1, 3 ZH(JJ) = -PXR(JJ) ENDDO VM = DSQRT(ZH(1)**2+ZH(2)**2+ZH(3)**2) DO JJ = 1, 3 ZH(JJ) = ZH(JJ)/VM ENDDO C C *** DEFINE POSITIVE PION FOUR MOMENTUM *** C DO JJ = 1, 4 VV(JJ) = EMU2(JJ) ENDDO CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,VV) DO JJ = 1, 4 PPIR(JJ) = VV(JJ) ENDDO VM = DSQRT(PPIR(1)**2+PPIR(2)**2+PPIR(3)**2) DO JJ = 1, 3 PPIR(JJ) = PPIR(JJ)/VM ENDDO C C *** COSINE THETA *** C DCTHETA = PPIR(1)*ZH(1) + PPIR(2)*ZH(2) + PPIR(3)*ZH(3) C C *** VSCAT: VECTOR PERPENDICULAR TO PRODUCTION PLANE *** C CALL DCROSS(YH,ZH,VSCAT) VM = DSQRT(VSCAT(1)**2+VSCAT(2)**2+VSCAT(3)**2) DO JJ = 1, 3 VSCAT(JJ) = VSCAT(JJ)/VM ENDDO C C *** VPROD: VECTOR PERPENDICULAR TO DECAY PLANE *** C CALL DCROSS(ZH,PPIR,VPROD) VM = DSQRT(VPROD(1)**2+VPROD(2)**2+VPROD(3)**2) DO JJ = 1, 3 VPROD(JJ) = VPROD(JJ)/VM ENDDO C C *** VSCPR: VECTOR PERPENDICULAR TO THE RHO DIRECTION WITHIN THE C PRODUCTION PLANE *** CALL DCROSS(VSCAT,ZH,VSCPR) VM = DSQRT(VSCPR(1)**2+VSCPR(2)**2+VSCPR(3)**2) DO JJ = 1, 3 VSCPR(JJ) = VSCPR(JJ)/VM ENDDO C *** SMALL PHI *** CPHIS = VSCAT(1)*VPROD(1)+VSCAT(2)*VPROD(2)+VSCAT(3)*VPROD(3) SPHIS = = -1.*(VSCPR(1)*VPROD(1)+VSCPR(2)*VPROD(2)+VSCPR(3)*VPROD(3)) IF ((CPHIS.EQ.0.) .AND. (SPHIS.EQ.0.)) THEN DPHIS = 0. ELSE DPHIS = DATAN2(SPHIS,CPHIS) ENDIF IF (DPHIS.LT.0.0) DPHIS = DPHIS + 2 * 3.14159265359d0 C *** VSCAT: VECTOR PERPENDICULAR TO SCATTERING PLANE *** CALL DCROSS(B,C,VSCAT) VM = DSQRT(VSCAT(1)**2+VSCAT(2)**2+VSCAT(3)**2) DO JJ = 1, 3 VSCAT(JJ) = VSCAT(JJ)/VM ENDDO C *** VPROD: VECTOR PERPENDICULAR TO PRODUCTION PLANE *** CALL DCROSS(B,A,VPROD) VM = DSQRT(VPROD(1)**2+VPROD(2)**2+VPROD(3)**2) DO JJ = 1, 3 VPROD(JJ) = VPROD(JJ)/VM ENDDO CALL DCROSS(VPROD,VSCAT,VSCPR) VM = DSQRT(B(1)**2+B(2)**2+B(3)**2) DO JJ = 1, 3 PVGUN(JJ) = B(JJ)/VM ENDDO C *** CAPITAL PHI *** SPHIL = VSCPR(1)*PVGUN(1) + VSCPR(2)*PVGUN(2) + VSCPR(3)*PVGUN(3) CPHIL = VSCAT(1)*VPROD(1) + VSCAT(2)*VPROD(2) + VSCAT(3)*VPROD(3) IF ((CPHIL.EQ.0.) .AND. (SPHIL.EQ.0.)) THEN DPHIL = 0. ELSE DPHIL = DATAN2(SPHIL,CPHIL) ENDIF IF (dPHIL.LT.0.0) DPHIL = DPHIL + 2 * 3.14159265359d0 c *** PSI *** DPSIS = DPHIS - DPHIL IF (DPSIS.LT.0.0) DPSIS = DPSIS + 2 * 3.14159265359d0 MCCOST = SNGL(DCTHETA) MCPHI = SNGL(DPHIS) MCPHIL = SNGL(DPHIL) MCPSI = SNGL(DPSIS) RETURN END SUBROUTINE HELOMEGA() C C------------------------------------------------------------------------- C C *** HELICITY ANGLES FOR THREE-BODY DECAY *** C C------------------------------------------------------------------------- C +CDE, DIPCO. REAL PPI1(4), PPI2(4), PG(4), PRHO(4) REAL PPR(4), XMRHO, THEL, PHIL, PCM(4), PX(4), AXIS(3) REAL PXR(4), ZH(3), YH(3), XH(3), PPIR1(4), PPIR2(4) REAL XPRO, YPRO, VDOT INTEGER I, JJ DOUBLE PRECISION BETAX, BETAY, BETAZ, GAMMA DOUBLE PRECISION A(5), B(5), C(5) DOUBLE PRECISION VPROD(3), VSCAT(3), VSCPR(3), PVGUN(3), VM DOUBLE PRECISION SPHIL, CPHIL, DPHIL, DPSIS DO I = 1, 4 PPI1(I) = SNGL(EMU2(I)) PPI2(I) = SNGL(EMU3(I)) PG(I) = SNGL(GAM(I)) PRHO(I) = SNGL(PSI(I)) PPR(I) = SNGL(EBP(I)) ENDDO XMRHO = PSIMAS DO I = 1, 5 A(I) = PSI(I) B(I) = GAM(I) C(I) = EBE(I) ENDDO C *** BOOST TO PROTON REST FRAME *** BETAX = EBP(1)/EBP(4) BETAY = EBP(2)/EBP(4) BETAZ = EBP(3)/EBP(4) GAMMA = EBP(4)/EBP(5) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,A) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,B) CALL LORENZ(-BETAX,-BETAY,-BETAZ,GAMMA,C) CALL VADD(PG,PPR,PCM,4) CALL VSUB(PCM,PRHO,PX,4) CALL CROSS(PG,PRHO,AXIS) CALL VUNIT(AXIS,AXIS,3) CALL LORENF(XMRHO,PRHO,PX,PXR) CALL VSCALE(PXR,-1.,ZH,3) CALL VUNIT(ZH,ZH,3) CALL UCOPY(AXIS,YH,3) CALL CROSS(YH,ZH,XH) CALL LORENF(XMRHO,PRHO,PPI1,PPIR1) CALL LORENF(XMRHO,PRHO,PPI2,PPIR2) CALL CROSS(PPIR1,PPIR2,AXIS) CALL VUNIT(AXIS,AXIS,3) C *** THETA AND SMALL PHI *** THEL=VDOT(AXIS,ZH,3) C PHIL=ATAN2(PPIR(2),PPIR(1)) XPRO=VDOT(AXIS,XH,3) YPRO=VDOT(AXIS,YH,3) PHIL=XPRO/SQRT(XPRO**2+YPRO**2) PHIL=ACOS(PHIL) IF(YPRO.LT.0.) PHIL=6.2832-PHIL C *** TOWARDS CAPITAL PHI AND PSI C *** VSCAT: VECTOR PERPENDICULAR TO SCATTERING PLANE *** CALL DCROSS(B,C,VSCAT) VM = DSQRT(VSCAT(1)**2+VSCAT(2)**2+VSCAT(3)**2) DO JJ = 1, 3 VSCAT(JJ) = VSCAT(JJ)/VM ENDDO C *** VPROD: VECTOR PERPENDICULAR TO PRODUCTION PLANE *** CALL DCROSS(B,A,VPROD) VM = DSQRT(VPROD(1)**2+VPROD(2)**2+VPROD(3)**2) DO JJ = 1, 3 VPROD(JJ) = VPROD(JJ)/VM ENDDO CALL DCROSS(VPROD,VSCAT,VSCPR) VM = DSQRT(B(1)**2+B(2)**2+B(3)**2) DO JJ = 1, 3 PVGUN(JJ) = B(JJ)/VM ENDDO C *** CAPITAL PHI *** SPHIL = VSCPR(1)*PVGUN(1) + VSCPR(2)*PVGUN(2) + VSCPR(3)*PVGUN(3) CPHIL = VSCAT(1)*VPROD(1) + VSCAT(2)*VPROD(2) + VSCAT(3)*VPROD(3) IF ((CPHIL.EQ.0.) .AND. (SPHIL.EQ.0.)) THEN DPHIL = 0. ELSE DPHIL = DATAN2(SPHIL,CPHIL) ENDIF IF (DPHIL.LT.0.0) DPHIL = DPHIL + 2 * 3.14159265359D0 C *** PSI *** DPSIS = PHIL - DPHIL IF (DPSIS.LT.0.0) DPSIS = DPSIS + 2 * 3.14159265359d0 MCCOST = THEL MCPHI = PHIL MCPHIL = SNGL(DPHIL) MCPSI = SNGL(DPSIS) RETURN END +DECK, DCROSS. SUBROUTINE DCROSS(A,B,PV) C C------------------------------------------------------------------------- C C *** VECTOR PRODUCT *** C C------------------------------------------------------------------------- C IMPLICIT NONE DOUBLE PRECISION A(3), B(3), PV(3) PV(1) = A(2)*B(3) - A(3)*B(2) PV(2) = A(3)*B(1) - A(1)*B(3) PV(3) = A(1)*B(2) - A(2)*B(1) RETURN END