C 06/05/87 710230854 MEMBER NAME LEPTO52 (LUND) M FORTRAN77 C######################################################################C C C C L E P T O C C C C version 5.2, July 31, 1987 C C C C The Lund Monte Carlo for Deep Inelastic Lepton-Nucleon Scattering C C C C Author: Gunnar Ingelman, DESY theory group, C C Notkestrasse 85, D-2000 Hamburg 52, FRG C C EARN/BITNET: T00ING @ DHHDESY3 C C Contribution on parton cascades from C C M. Bengtsson, T. Sjoestrand, Lund University C C C C Documentation: short manual in T00ING.LUND(INFOL52) C C Source code in T00ING.LUND(LEPTO52) C C Please report any errors or suggestions for impromevents. C C C C######################################################################C BLOCK DATA LEPTOD C...Give sensible default values to switches and parameters. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LFLMIX/ CABIBO(4,4) COMMON /LGRID/ NXX,NWW,XX(15),WW(15),PQG(15,15,3),PQQB(15,15,2), &QGMAX(15,15,3),QQBMAX(15,15,2),YCUT(15,15),XTOT(15,15),NP COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),COMFAC COMMON /PYPARA/ IPY(80),PYPAR(80),PYVAR(80) COMMON /PYMINU/ XKIN(4),UKIN(4),WKIN(4),AIN(4),BIN(4), &MAXFIN,RELUP,RELERR,RELER2,FCNMAX COMMON /PYMINC/ NAMKIN(4),NAM(30) CHARACTER*10 NAMKIN,NAM C...LEPTOU: Cuts, basic switches and parameters. DATA CUT/0.001,1.,0.,1.,4.,1.E+08,5.,1.E+08,1.,1.E+08,1.,1.E+08, &0.,3.1416/ DATA LST/4*1,3,3*1,4,1,4,4,4,2*1,0,2,0,3,21*0/ DATA PARL/2*1.,.44,.75,.226,0.,4.,0.0025,2.,.3,.01,4.,.1,.44, &.01,7.2973E-03,1.16637E-05,0.0711,12*0./ C...Internally used variables. DATA PARI/40*0./ DATA QC/.66667,2*-.33333,.66667,-.33333,.66667,-.33333,.66667/ DATA CABIBO/.95,.05,2*0.,.05,.948,.002,2*0.,.002,.998,4*0.,1./ DATA NWW,NXX/2*15/ DATA PQG,PQQB,QGMAX,QQBMAX/2250*0./,YCUT/225*0./,XTOT/225*0./ DATA XKIN/1.,2.,3.,4./,UKIN,WKIN,AIN,BIN/16*0./,MAXFIN/2000/ DATA RELUP,RELERR,RELER2/0.1,0.05,0.05/ DATA NAMKIN/' x',' Q**2 (y)',' ',' '/ DATA IPY/ 1 0, 0, 2, 2, 6, 1, 1, 6, 3, 1, 2 3, 1, 1, 2, 1, 1, 4, 1, 1, 1, 3 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 4 1, 2, 1, 1, 30, 33, 1, 1, 7, 0, 5 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6 0, 0, 0, 1, 100, 0, 0, 0, 0, 0, 7 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/ DATA (PYPAR(I),I=1,40)/ 1 7.299E-03, 2.290E-01, 2.000E-01, 2.500E-01, 4.000E+00, 1 1.000E+00, 4.400E-01, 4.400E-01, 7.500E-02, 0.000E+00, 2 2.000E+00, 2.000E+00, 1.000E+00, 0.000E+00, 3.000E+00, 2 1.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.000E+00, 3 2.500E-01, 1.000E+00, 2.000E+00, 1.000E-03, 1.000E+00, 3 1.000E+00, 1.000E+00, -2.000E-02, -1.000E-02, 0.000E+00, 4 0.000E+00, 1.600E+00, 0.500E+00, 0.200E+00, 3.894E-01, 4 1.000E+00, 3.300E-01, 6.600E-01, 0.000E+00, 1.000E+00/ DATA (PYPAR(I),I=41,80)/ 5 2.260E+00, 1.000E+04, 1.000E-04, 0.000E+00, 0.000E+00, 5 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 6 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 6 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 7 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 7 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 8 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 8 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/ DATA PYVAR/80*0./ END C ********************************************************************** SUBROUTINE LINIT(LFILE,LEPIN,PLX,PLY,PLZ,PPX,PPY,PPZ,INTER) C...Initialize for an incoming lepton (type LEPIN, momentum PLX,PLY,PLZ) C...and target nucleon (momentum PPX,PPY,PPZ) to interact via INTER. C...Find maximum of differential cross section, calculate QCD event C...probabilities or read them from logical file LFILE (if >0). C...Numerical integration to obtain total cross-section. COMMON /LINTEG/ NTOT,NPASS COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LGRID/ NXX,NWW,XX(15),WW(15),PQG(15,15,3),PQQB(15,15,2), &QGMAX(15,15,3),QQBMAX(15,15,2),YCUT(15,15),XTOT(15,15),NP COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),COMFAC COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT1/ MST(40),PAR(80) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR COMMON /PYMINU/ XKIN(4),UKIN(4),WKIN(4),AIN(4),BIN(4), &MAXFIN,RELUP,RELERR,RELER2,FCNMAX COMMON /PYMINC/ NAMKIN(4),NAM(30) COMMON /PYPARA/ IPY(80),PYPAR(80),PYVAR(80) CHARACTER*10 NAMKIN,NAM DIMENSION LSTW(40),PARLW(30) DOUBLE PRECISION DTHETA,DPHI,DBETA,DE,DPZ,DPT EXTERNAL DCROSS,DLOWER,DUPPER DATA PI/3.1415927/ IF(LST(8).LT.2) THEN C...Fragmentation parameters suitable for use without parton cascades, C...default values in JETSET 6.3 matched with cascade. PAR(12)=0.4 PAR(23)=1.1 PAR(25)=1.1 PAR(31)=1. PAR(32)=0.7 PAR(33)=1. PAR(34)=0.7 PAR(35)=0.5 PAR(36)=0.7 ELSE C...Reset PYTHIA parameters from LEPTO parameters. IF(LST(8).EQ.3.OR.LST(8).EQ.5) IPY(13)=0 IF(LST(8).EQ.4.OR.LST(8).EQ.5) IPY(14)=0 IPY(8)=LST(12) IPY(77)=LST(17) PYPAR(7)=PARL(3) ENDIF IF(LST(18).NE.1) THEN C...W, Z masses from theta-Weinberg, Fermi constant GF and rad. corr. PMAS(3)=SQRT(PI*PARL(16)/(SQRT(2.)*PARL(17)*PARL(5)* & (1.-PARL(18)))) PMAS(2)=PMAS(3)/SQRT(1.-PARL(5)) ENDIF C...Couplings between Z0 and left/right-handed leptons and quarks. ZL(1,1)=-.5+PARL(5) ZL(1,2)=PARL(5) ZL(2,1)=ZL(1,2) ZL(2,2)=ZL(1,1) ZL(1,3)=0.5 ZL(2,3)=0. ZL(1,4)=0. ZL(2,4)=0.5 DO 10 IFL=1,8 ZQ(1,IFL)=SIGN(0.5,QC(IFL))-QC(IFL)*PARL(5) 10 ZQ(2,IFL)=-QC(IFL)*PARL(5) C...Set initial state. LST(1)=INTER KSAVE(1)=LEPIN KSAVE(2)=41 K(1,1)=40000 K(2,1)=40000 K(1,2)=KSAVE(1) K(2,2)=KSAVE(2) P(1,1)=PLX P(1,2)=PLY P(1,3)=PLZ P(1,5)=ULMASS(0,LEPIN) P(1,4)=SQRT(P(1,1)**2+P(1,2)**2+P(1,3)**2+P(1,5)**2) P(2,1)=PPX P(2,2)=PPY P(2,3)=PPZ P(2,5)=ULMASS(0,41) P(2,4)=SQRT(P(2,1)**2+P(2,2)**2+P(2,3)**2+P(2,5)**2) N=2 C...Save momentum vectors of incoming particles; for lepton to be used C...when making cuts on the angle of outgoing lepton w.r.t. incoming C...lepton in lab system, for nucleon to be used when varying a lepton C...beam energy. DO 20 I=1,2 DO 20 J=1,5 20 PLAB(I,J)=P(I,J) C...Transform to cms of incoming particles, lepton along +z axis. DO 30 J=1,3 30 DBETA(1,J)=DBLE(P(1,J)+P(2,J))/DBLE(P(1,4)+P(2,4)) CALL DUROBO(0.D0,0.D0,-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(P(1,1),P(1,2)) CALL DUROBO(0.D0,-DPHI(1),0.D0,0.D0,0.D0) DTHETA(1)=ULANGL(P(1,3),P(1,1)) CALL DUROBO(-DTHETA(1),0.D0,0.D0,0.D0,0.D0) C...Save momentum vectors of incoming particles in their cm frame, C...to be used when initializing each new event to be simulated. DO 40 I=1,2 DO 40 J=1,5 40 PCMS(I,J)=P(I,J) C...CMS energy squared (neglecting initial state particle masses). PARL(21)=2.*(P(1,4)*P(2,4)-P(1,3)*P(2,3)) WRITE(6,1000) LEPIN,PLX,PLY,PLZ,PARL(1),PARL(2),PPX,PPY,PPZ,INTER, &SQRT(PARL(21)) C...Effective limits on kinematic variables x, y, Q**2. PM2=P(2,5)**2 S=PARL(21) CUT(6)=MIN(CUT(6),S) CUT(8)=MIN(CUT(8),S) CUT(10)=MIN(CUT(10),S/2.) XMIN=MAX(CUT(1),CUT(5)/(S*CUT(4)),CUT(5)/(2.*P(2,5)*CUT(10)), &1.-(CUT(8)-PM2)/(S*MAX(CUT(3),1.E-10))) YMIN=MAX(CUT(3),CUT(5)/(S*CUT(2)),(CUT(7)-PM2)/(S*(1.-CUT(1)))) Q2MIN=MAX(CUT(5),S*CUT(1)*CUT(3),S*CUT(3)-CUT(8)+PM2, &2.*P(2,5)*CUT(9)*CUT(1)) XMAX=MIN(CUT(2),CUT(6)/(S*MAX(CUT(3),1.E-10)), &1.-(CUT(7)-PM2)/(S*CUT(4)),CUT(6)/(2.*P(2,5)*MAX(CUT(9),1.E-10))) YMAX=MIN(CUT(4),CUT(6)/(S*MAX(CUT(1),1.E-10)), &(CUT(8)-PM2)/(S*MAX(1.-CUT(2),1.E-06))) Q2MAX=MIN(CUT(6),S*CUT(2)*CUT(4),S*CUT(4)-CUT(7)+PM2, &2.*P(2,5)*CUT(10)*CUT(2)) WRITE(6,1050) CUT,XMIN,XMAX,YMIN,YMAX,Q2MIN,Q2MAX IF(XMAX.LT.XMIN.OR.YMAX.LT.YMIN.OR.Q2MAX.LT.Q2MIN) THEN WRITE(6,1100) STOP ENDIF PARI(11)=(PARL(1)-PARL(2))/PARL(1) KSAVE(4)=LEPIN ILEP=1 IF(LEPIN.LT.0) ILEP=2 INU=0 IF(IABS(LEPIN).EQ.8.OR.IABS(LEPIN).EQ.10) INU=1 IF(INU.EQ.1) THEN C...Set full polarisation for incoming neutrino. PARL(6)=-1. IF(LEPIN.LT.0) PARL(6)=1. ENDIF IF(LST(1).EQ.1.AND.INU.EQ.0) THEN C...Electromagnetic interaction. LST(23)=1 KSAVE(3)=1 IG=1 IZ=0 ELSEIF(LST(1).EQ.2) THEN C...Weak charged current, only one helicity state contributes. IF(KSAVE(1).LT.0.AND.PARL(6).LT.-0.99 & .OR.KSAVE(1).GT.0.AND.PARL(6).GT.0.99) THEN WRITE(6,1150) LEPIN,PARL(6) STOP ENDIF LST(23)=1+ILEP IF(MOD(IABS(LEPIN),2).EQ.0) THEN KSAVE(3)=ISIGN(3,LEPIN) KSAVE(4)=ISIGN(IABS(LEPIN)-1,LEPIN) ELSE KSAVE(3)=ISIGN(3,-LEPIN) KSAVE(4)=ISIGN(IABS(LEPIN)+1,LEPIN) ENDIF ELSEIF(LST(1).EQ.3.OR.(LST(1).EQ.4.AND.INU.EQ.1)) THEN C...Weak neutral current. LST(23)=3+ILEP KSAVE(3)=2 IG=0 IZ=1 ELSEIF(LST(1).EQ.4.AND.INU.EQ.0) THEN C...Neutral current, electromagnetic and weak with interference. LST(23)=6 KSAVE(3)=5 IG=1 IZ=1 ELSE WRITE(6,1200) INTER,LEPIN STOP ENDIF C...Set process dependent optimization parameters. DO 50 I=1,4 OPTX(I)=0. OPTY(I)=0. 50 OPTQ2(I)=0. LST(31)=1 IF(LST(1).EQ.1) THEN LST(31)=0 OPTX(2)=1. OPTY(1)=1. OPTQ2(3)=1. ELSEIF(LST(1).EQ.4) THEN LST(31)=0 OPTX(1)=0.1 OPTX(2)=1. OPTY(1)=1. OPTQ2(1)=0.5 OPTQ2(2)=0.5 OPTQ2(3)=1. ELSE OPTX(1)=1. OPTY(1)=1. OPTQ2(1)=1. ENDIF C...Initialize Monte Carlo estimate of cross section. PARL(24)=0. PARI(27)=0. PARI(28)=0. PARI(29)=0. PARI(30)=0. PARI(32)=0. IF(LST(1).EQ.2) THEN C...Constant factor GF**2/pi for CC, transformation to picobarn. PARI(31)=PARL(17)**2/PI*0.39E+09 ELSE C...Constant factor 2*pi*alpha**2 for NC, transformation to picobarn. PARI(31)=2.*PI*PARL(16)**2*0.39E+09 ENDIF WRITE(6,1250) (I,LST(I),LST(I+10),PARL(I),PARL(I+10),I=1,10) PARL(23)=0. IF(LST(10).EQ.1) THEN C...Get integrated cross-section through adaptive Gaussian integration. WRITE(6,1300) ACCUR=PARL(15) IT=0 100 IT=IT+1 NTOT=0 NPASS=0 EPS=ACCUR CALL PYTIME(TI1) CALL GADAP2(XMIN,XMAX,DLOWER,DUPPER,DCROSS,EPS,PARL(23)) CALL PYTIME(TI2) WRITE(6,1310) IT,NTOT,NPASS,TI2-TI1,PARL(23) IF(PARL(23).GT.1.) THEN WRITE(6,1320) ACCUR ELSE WRITE(6,1330) ACCUR,ACCUR/PARL(23),PARL(15) ACCUR=PARL(23)*PARL(15) IF(IT.LT.2) GOTO 100 ENDIF ENDIF IF(LST(2).EQ.1) THEN C...Find max value of differential cross section for rejection. UKIN(1)=(XMAX+XMIN)/2. WKIN(1)=0.8*(XMAX-XMIN)/2. AIN(1)=XMIN BIN(1)=XMAX IF(LST(31).EQ.0) THEN UKIN(2)=(Q2MAX+Q2MIN)/2. WKIN(2)=0.8*(Q2MAX-Q2MIN)/2. AIN(2)=Q2MIN BIN(2)=Q2MAX NAMKIN(2)=' Q**2' ELSE UKIN(2)=(YMAX+YMIN)/2. WKIN(2)=0.8*(YMAX-YMIN)/2. AIN(2)=YMIN BIN(2)=YMAX NAMKIN(2)=' y' ENDIF C...Maximum obtained by minimizing -(diff. x-section). CALL PYTIME(TI1) CALL PYMINN CALL PYTIME(TI2) PARI(LST(23))=FCNMAX*1.1 WRITE(6,1400) PARI(LST(23)),TI2-TI1 ENDIF IF(LFILE.GT.0) THEN C...Read QCD weights from file. READ(LFILE) LSTW,PARLW,NXX,NWW,NP,XX,WW IPMAX=2 IF(LSTW(16).NE.0) IPMAX=3 READ(LFILE) (((PQG(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,NP), & (((PQQB(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,NP), & (((QGMAX(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,IPMAX), & (((QQBMAX(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,MIN(2,IPMAX)), & YCUT IF(NP.NE.1) READ(LFILE) XTOT CLOSE(LFILE) C...Check current parameter values against those used when C...calculating weights. IF(LST(12).NE.LSTW(12).OR.LST(13).NE.LSTW(13) & .OR.LST(15).NE.LSTW(15).OR.LST(16).NE.LSTW(16) & .OR.LST(23).NE.LSTW(23) & .OR.ABS(PARL(1)-PARLW(1)).GT.0.1.OR.ABS(PARL(2)-PARLW(2)).GT.0.1 & .OR.ABS(PARL(5)-PARLW(5)).GT.0.01 & .OR.ABS(PARL(6)-PARLW(6)).GT.0.1) THEN WRITE(6,1500) LST(12),LSTW(12),LST(13),LSTW(13),LST(15), & LSTW(15),LST(16),LSTW(16),LST(23),LSTW(23),PARL(1),PARLW(1), & PARL(2),PARLW(2),PARL(5),PARLW(5),PARL(6),PARLW(6) STOP ENDIF ELSEIF(LST(8).EQ.1) THEN C...Calculate weights if 1st order QCD requested. CALL PYTIME(TI1) CALL LWEITS(LFILE) CALL PYTIME(TI2) WRITE(6,1510) TI2-TI1 ENDIF C...Reset counter to zero for Monte Carlo estimate of cross section. PARI(28)=0. LST(20)=0 RETURN 1000 FORMAT('1',//,5X,'THE LUND MONTE CARLO FOR DEEP INELASTIC LEPTON-' &,'NUCLEON SCATTERING',/,5X,65('='),//, &25X,'LEPTO version 5.2, July 31, 1987',//, &' Lepton: type =',I3,5X,'momentum (px,py,pz) =',3F8.1, &' GeV/c',//,' Target: A, Z =',2F3.0,2X, &'momentum (px,py,pz) =',3F8.1,' GeV/c',//, &' Interaction :',I3,//, &' cms energy =',G12.4,/) 1050 FORMAT(/,' User applied cuts: ', & F12.4,' < x < ',F12.4, &/,20X,F12.4,' < y < ',F12.4, &/,20X,G12.4,' < Q**2 < ',G12.4, &/,20X,G12.4,' < W**2 < ',G12.4, &/,20X,G12.4,' < nu < ',G12.4, &/,20X,G12.4,' < E'' < ',G12.4, &/,20X,G12.4,' < theta < ',G12.4,/, &/, ' Effective ranges: ', & F12.4,' < x < ',F12.4, &/,20X,F12.4,' < y < ',F12.4, &/,20X,G12.4,' < Q**2 < ',G12.4) 1100 FORMAT(' Effective max value smaller that effective min ', &'value. Execution stopped !',/) 1150 FORMAT(/,' Weak charged current cross section zero for specified', &' lepton helicity; LEPIN, PARL(6) =',I3,F5.2, &/,' Execution stopped!',/) 1200 FORMAT(' Unrecognized interaction in LINIT call; INTER =',I5, &' for lepton LEPIN =',I5,/,' Execution stopped !',/) 1250 FORMAT(/,' Parameter values:',//,9X,'I',4X,'LST(I)',1X, &'LST(I+10)',8X,'PARL(I)',5X,'PARL(I+10)', &/,5X,55('-'),10(/,3I10,2G15.4),/) 1300 FORMAT(/,' Adaptive Gaussian integration of cross section', & /,' ----------------------------------------------'/,) 1310 FORMAT(' Iteration #',I3, &/,' Total # of function evaluations =',I8, &/,' # of non-zero function evaluations =',I8, &/,' Time for integration =',F8.3,' seconds.', &/,' ===> Cross-section (in region allowed by cuts) =',G10.2,' pb') 1320 FORMAT(' Required relative error = ',G10.2,/) 1330 FORMAT(' Effective absolute error = ',G10.2,/, & ' Effective relative error = ',G10.2,/, & ' Required relative error = ',G10.2,/) 1400 FORMAT(' Max of differential cross section (for weighting) =', &E12.4,/,' obtained in ',F5.1,' seconds.',/) 1500 FORMAT(//,' Warning: current parameter values do not match ', &'with those used when calculating weights.',//,15X, &'current value value for weights',/, &/,' LST(12) ',I12,10X,I12, &/,' LST(13) ',I12,10X,I12, &/,' LST(15) ',I12,10X,I12, &/,' LST(16) ',I12,10X,I12, &/,' LST(23) ',I12,10X,I12, &/,' PARL(1) ',E12.4,10X,E12.4, &/,' PARL(2) ',E12.4,10X,E12.4, &/,' PARL(5) ',E12.4,10X,E12.4, &/,' PARL(6) ',E12.4,10X,E12.4, &/,' Execution stopped !',/) 1510 FORMAT(/,' Time for calculating QCD weights =',F5.1,' seconds',/) END C ********************************************************************** SUBROUTINE LEPTO C...Administer the generation of an event. C...Note that if the error flag LST(21) is not zero, no event has been C...generated. COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /GITEST/ PEPCMS(4,5),TLAB(2,5),KQ(20),KQG(20),KQQ(20) COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT1/ MST(40),PAR(80) COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR DOUBLE PRECISION DTHETA,DPHI,DBETA,DE,DPZ,DPT,DETOT DIMENSION SPQ(17) 1 LST(21)=0 DO 5 I=1,4 5 K(I,2)=KSAVE(I) N=2 IF(LST(16).NE.0.AND.LST(2).GT.0) THEN C...Lepton energy is allowed to vary from event to event, C...momentum vector taken from P(1,1), P(1,2) and P(1,3). DO 10 I=1,2 DO 10 J=1,5 10 P(I,J)=PLAB(I,J) P(1,5)=ULMASS(0,K(1,2)) P(1,4)=SQRT(P(1,1)**2+P(1,2)**2+P(1,3)**2+P(1,5)**2) DO 20 J=1,5 20 PLAB(1,J)=P(1,J) C...Transform to cms of incoming particles, lepton along +z axis. DO 30 J=1,3 30 DBETA(1,J)=DBLE(P(1,J)+P(2,J))/DBLE(P(1,4)+P(2,4)) CALL DUROBO(0.D0,0.D0,-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(P(1,1),P(1,2)) CALL DUROBO(0.D0,-DPHI(1),0.D0,0.D0,0.D0) DTHETA(1)=ULANGL(P(1,3),P(1,1)) CALL DUROBO(-DTHETA(1),0.D0,0.D0,0.D0,0.D0) PARL(21)=2.*(P(1,4)*P(2,4)-P(1,3)*P(2,3)) ELSE C...Initial state momenta fixed from LINIT call. DO 40 I=1,2 DO 40 J=1,5 40 P(I,J)=PCMS(I,J) ENDIF CALL LEPTOX C...Return if error or if no event to be generated. IF(LST(21).NE.0.OR.LST(2).LE.0.OR.LST(3).EQ.0) RETURN PARI(29)=PARI(29)+1. C...Scattered lepton and exchanged virtual boson added to event record. P(4,1)=SNGL(DPT) P(4,2)=0. P(4,3)=SNGL(DPZ) P(4,5)=ULMASS(0,K(4,2)) P(3,1)=-P(4,1) P(3,2)=0. P(3,3)=P(1,3)-P(4,3) P(3,4)=P(1,4)-P(4,4) P(3,5)=-SQRT(Q2) K(1,1)=40000 K(2,1)=40000 K(3,1)=50000 K(4,1)=0 C...Inactivate scattered lepton in event record if desired. IF(LST(4).EQ.0) K(4,1)=40000 N=4 C...Prepare for parton cascade. IF(LST(8).GE.2) CALL CASCAD(0) CGITEST DO 19 I=1,4 DO 19 J=1,5 IF(I.LE.2) TLAB(I,J)=PLAB(I,J) 19 PEPCMS(I,J)=P(I,J) C...Transform to hadronic cms. Boost parameters in double precision. DETOT=DBLE(P(1,4))-DE+DBLE(P(2,4)) DBETA(2,1)=-DPT/DETOT DBETA(2,2)=0.D+0 DBETA(2,3)=(DBLE(P(1,3))-DPZ+DBLE(P(2,3)))/DETOT CALL DUROBO(0.D0,0.D0,-DBETA(2,1),-DBETA(2,2),-DBETA(2,3)) DPHI(2)=0. DTHETA(2)=ULANGL(P(3,3),P(3,1)) CALL DUROBO(-DTHETA(2),0.D0,0.D0,0.D0,0.D0) C...Save momentum of exchanged boson (used in subroutine LFRAME). DO 120 J=1,5 120 PB(J)=P(3,J) MST(1)=N+1 LST(26)=MST(1) PARL(25)=DALPS(Q2) IF(LST(8).NE.1) THEN C...QPM model without QCD corrections (cascade applied later). 200 CALL LQEV IF(LST(21).NE.0) GOTO 200 ELSE DO 210 I=1,17 210 SPQ(I)=PQ(I) C...Probabilities for hard, first order QCD events. CALL LQCDP(QG,QQB) 300 SRLU=RLU(0) IF(SRLU.GT.QQB+QG) THEN DO 310 I=1,17 310 PQ(I)=SPQ(I) CALL LQEV ELSEIF(SRLU.GT.QQB) THEN CALL LQGEV ELSE CALL LQQBEV ENDIF IF(LST(21).NE.0) GOTO 300 ENDIF IF(LST(8).LT.2) THEN C...No parton cascade, introduce primordial kt. IF(PARL(3).GT.1.E-03) THEN CALL LPRIKT(PARL(3),PT,PHI) CALL DUROBO(0.D0,DBLE(-PHI),0.D0,0.D0,0.D0) CALL DUROBO(DBLE(ATAN(2.*PT/SQRT(W2))),DBLE(PHI),0.D0,0.D0,0.D0) ENDIF MST(1)=0 C...Check system against fragmentation cuts. MST(26)=0 CALL LUPREP IF(MST(26).NE.0) THEN WRITE(6,*) ' Error in LUPREP! MST(25), MST(26) = ',MST(25), & MST(26) LST(21)=11 ENDIF ELSE MST(1)=0 C...Include parton cascades, target remnant and primordial kt. CALL CASCAD(1) ENDIF IF(LST(21).NE.0) GOTO 1 PARI(30)=PARI(30)+1. IF(LST(1).EQ.2) PARL(24)=PARL(24)*PARI(30)/PARI(29) IF(LST(7).EQ.1) CALL LUEXEC LST(28)=1 LST(29)=0 PHIR=6.2832*RLU(0) IF(LST(5).GE.2) CALL LFRAME(LST(5),LST(6)) RETURN END C ********************************************************************** SUBROUTINE LEPTOX C...Select type of process and choose kinematical variables (x,y or x,Q* C...according to the relevant total cross section. COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),COMFAC COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) DIMENSION XPQ(-6:6),PQH(17,2),PNT(2,2) DO 10 IH=1,2 DO 5 I=1,2 5 PNT(I,IH)=0. DO 6 I=1,8 EWQC(1,IH,I)=0. 6 EWQC(2,IH,I)=0. DO 10 I=1,17 10 PQH(I,IH)=0. DO 20 I=1,17 20 PQ(I)=0. LST(21)=0 IF(LST(2).NE.1) GOTO 110 100 PARI(28)=PARI(28)+1 101 CONTINUE C...Choose x according to the distribution C...hx(x) = a + b/x + c/x**2 + d/x**3. In detail C...hq=OPTX(1)/(XMAX-XMIN) + 1/ln(XMAX/XMIN)*OPTX(2)/X C... +XMIN*XMAX/(XMAX-XMIN)*OPTX(3)/X**2 C... +2*(XMIN*XMAX)**2/(XMAX**2-XMIN**2)*OPTX(4)/X**3 WHICH=(OPTX(1)+OPTX(2)+OPTX(3)+OPTX(4))*RLU(0) IF(WHICH.LE.OPTX(1)) THEN X=XMIN+RLU(0)*(XMAX-XMIN) ELSEIF(WHICH.LE.(OPTX(1)+OPTX(2))) THEN X=XMIN*(XMAX/XMIN)**RLU(0) ELSEIF(WHICH.LE.(OPTX(1)+OPTX(2)+OPTX(3))) THEN X=XMIN*XMAX/(XMAX+RLU(0)*(XMIN-XMAX)) ELSE X=SQRT((XMIN*XMAX)**2/(XMAX**2+RLU(0)*(XMIN**2-XMAX**2))) ENDIF IF(LST(31).EQ.0) THEN C...Choose Q**2 according to the distribution C...hq(Q2) = a + b/(Q2) + c/(Q2)**2 + d/(Q2)**3. In detail C...hq=OPTQ2(1)/(Q2MAX-Q2MIN) + 1/ln(Q2MAX/Q2MIN)*OPTQ2(2)/Q2 C... +Q2MIN*Q2MAX/(Q2MAX-Q2MIN)*OPTQ2(3)/Q2**2 C... +2*(Q2MIN*Q2MAX)**2/(Q2MAX**2-Q2MIN**2)*OPTQ2(4)/Q2**3 Q2LOW=MAX(Q2MIN,X*YMIN*PARL(21)) Q2UPP=MIN(Q2MAX,X*YMAX*PARL(21)) IF(Q2UPP.LT.Q2LOW) GOTO 101 WHICH=(OPTQ2(1)+OPTQ2(2)+OPTQ2(3)+OPTQ2(4))*RLU(0) IF(WHICH.LE.OPTQ2(1)) THEN Q2=Q2LOW+RLU(0)*(Q2UPP-Q2LOW) ELSEIF(WHICH.LE.(OPTQ2(1)+OPTQ2(2))) THEN Q2=Q2LOW*(Q2UPP/Q2LOW)**RLU(0) ELSEIF(WHICH.LE.(OPTQ2(1)+OPTQ2(2)+OPTQ2(3))) THEN Q2=Q2LOW*Q2UPP/(Q2UPP+RLU(0)*(Q2LOW-Q2UPP)) ELSE Q2=SQRT((Q2LOW*Q2UPP)**2/(Q2UPP**2+RLU(0)*(Q2LOW**2-Q2UPP**2))) ENDIF Y=Q2/(PARL(21)*X) IF(Y.LT.YMIN.OR.Y.GT.YMAX) GOTO 100 ELSE C...Choose y according to the distribution C...hy(y) = a + b/y + c/y**2 + d/y**3. In detail C...hy=OPTY(1)/(YMAX-YMIN) + 1/ln(YMAX/YMIN)*OPTY(2)/Y C... +YMIN*YMAX/(YMAX-YMIN)*OPTY(3)/Y**2 C... +2*(YMIN*YMAX)**2/(YMAX**2-YMIN**2)*OPTY(4)/Y**3 WHICH=(OPTY(1)+OPTY(2)+OPTY(3)+OPTY(4))*RLU(0) IF(WHICH.LE.OPTY(1)) THEN Y=YMIN+RLU(0)*(YMAX-YMIN) ELSEIF(WHICH.LE.(OPTY(1)+OPTY(2))) THEN Y=YMIN*(YMAX/YMIN)**RLU(0) ELSEIF(WHICH.LE.(OPTY(1)+OPTY(2)+OPTY(3))) THEN Y=YMIN*YMAX/(YMAX+RLU(0)*(YMIN-YMAX)) ELSE Y=SQRT((YMIN*YMAX)**2/(YMAX**2+RLU(0)*(YMIN**2-YMAX**2))) ENDIF Q2=X*Y*PARL(21) IF(Q2.LT.Q2MIN.OR.Q2.GT.Q2MAX) GOTO 100 ENDIF 110 IF(LKINEM(LST(2)).NE.0) THEN IF(LST(2).EQ.1) GOTO 100 LST(21)=1 RETURN ENDIF PARI(24)=(1.+(1.-Y)**2)/2. PARI(25)=1.-Y PARI(26)=(1.-(1.-Y)**2)/2. CALL LNUCL(X,Q2,XPQ) C...Lepton helicity state, only one contributes in some cases. IH=1 IF(PARL(6).GT.+0.99) IH=2 200 LST(30)=SIGN(1.,IH-1.5) PQH(17,IH)=0. PNT(1,IH)=0. PNT(2,IH)=0. IF(LST(1).EQ.2) THEN C...Zero CC cross-section for one helicity state. IF(KSAVE(1).LT.0.AND.IH.EQ.1 & .OR.KSAVE(1).GT.0.AND.IH.EQ.2) GOTO 240 YQ=PARI(24)-LST(30)*PARI(26) YQB=PARI(24)+LST(30)*PARI(26) IF(PARI(11).GT.1.E-06) THEN PNT(1,IH)=(1.-PARI(11))*(PARI(12)+PARI(13))*YQ PNT(2,IH)=PARI(11)*(PARI(12)+PARI(13))*YQ ENDIF DO 220 I=1,LST(12) IF(K(3,2)*QC(I).LT.0) THEN PQH(I,IH)=XPQ(I)*YQ ELSE PQH(I+LST(12),IH)=XPQ(-I)*YQB ENDIF 220 CONTINUE ELSE C...Neutral current, electromagnetic and weak with interference. GFQ2=Q2/(PMAS(2)**2+Q2)/PARL(5)/(1.-PARL(5)) II=3-IH ZLEP=ZL(IH,ILEP+2*INU) DO 230 I=1,MAX(LST(12),LST(13)) A=(-QC(I)*IG+IZ*GFQ2*ZLEP*ZQ(IH,I))**2 B=(-QC(I)*IG+IZ*GFQ2*ZLEP*ZQ(II,I))**2 C...Save helicity-dependent electroweak quark couplings for later use. EWQC(1,IH,I)=A EWQC(2,IH,I)=B IF(I.GT.LST(12)) GOTO 230 FYQ=(A+B)*PARI(24)+(A-B)*PARI(26) PQH(I,IH)=XPQ(I)*FYQ IF(I.LE.2.AND.PARI(11).GT.1.E-06) THEN PNT(1,IH)=PNT(1,IH)+(1.-PARI(11))*PARI(11+I)*FYQ PNT(2,IH)=PNT(2,IH)+PARI(11)*PARI(14-I)*FYQ ENDIF PQH(I+LST(12),IH)=XPQ(-I)*((A+B)*PARI(24)-(A-B)*PARI(26)) 230 CONTINUE ENDIF 240 CONTINUE DO 300 I=1,LST(12) 300 PQH(17,IH)=PQH(17,IH)+PQH(I,IH)+PQH(I+LST(12),IH) IF(ABS(PARL(6)).LT.0.99.AND.IH.EQ.1) THEN IH=2 GOTO 200 ENDIF DO 310 I=1,17 310 PQ(I)=(1.-PARL(6))/2.*PQH(I,1)+(1.+PARL(6))/2.*PQH(I,2) C...Common factor for matrix elements. IF(LST(31).EQ.0) THEN IF(LST(1).EQ.2) THEN COMFAC=1./X/(1.+Q2/PMAS(3)**2)**2 ELSE COMFAC=1./X/Q2**2 ENDIF ELSE IF(LST(1).EQ.2) THEN COMFAC=1./(1.+Q2/PMAS(3)**2)**2*PARL(21) ELSE COMFAC=1./Q2**2*PARL(21) ENDIF ENDIF IF(LST(2).LE.-2) RETURN HX=OPTX(1)/(XMAX-XMIN) + 1./ALOG(XMAX/XMIN)*OPTX(2)/X &+XMIN*XMAX/(XMAX-XMIN)*OPTX(3)/X**2 &+2*(XMIN*XMAX)**2/(XMAX**2-XMIN**2)*OPTX(4)/X**3 XFACT=OPTX(1)+OPTX(2)+OPTX(3)+OPTX(4) IF(LST(31).EQ.0) THEN Q2LOW=MAX(Q2MIN,X*YMIN*PARL(21)) Q2UPP=MIN(Q2MAX,X*YMAX*PARL(21)) HQ2=OPTQ2(1)/(Q2UPP-Q2LOW) & +1./ALOG(Q2UPP/Q2LOW)*OPTQ2(2)/Q2 & +Q2LOW*Q2UPP/(Q2UPP-Q2LOW)*OPTQ2(3)/Q2**2 & +2*(Q2LOW*Q2UPP)**2/(Q2UPP**2-Q2LOW**2)*OPTQ2(4)/Q2**3 Q2FACT=OPTQ2(1)+OPTQ2(2)+OPTQ2(3)+OPTQ2(4) COMFAC=COMFAC*XFACT*Q2FACT/HX/HQ2 ELSE HY=OPTY(1)/(YMAX-YMIN)+1./ALOG(YMAX/YMIN)*OPTY(2)/Y & +YMIN*YMAX/(YMAX-YMIN)*OPTY(3)/Y**2 & +2*(YMIN*YMAX)**2/(YMAX**2-YMIN**2)*OPTY(4)/Y**3 YFACT=OPTY(1)+OPTY(2)+OPTY(3)+OPTY(4) COMFAC=COMFAC*XFACT*YFACT/HX/HY ENDIF IF(LST(2).LE.0) RETURN SIGL=(1.-PARL(6))/2.*PQH(17,1) SIGR=(1.+PARL(6))/2.*PQH(17,2) SIGMA=SIGL+SIGR IF(LST(2).EQ.1) THEN C...When chosing (x,y), reject according to maximum of "cross-section", C...update cross-section estimate. PARI(27)=PARI(27)+SIGMA*COMFAC VIOL=SIGMA*COMFAC/PARI(LST(23)) IF(VIOL.GT.PARI(32)) THEN PARI(32)=VIOL IF(PARI(32).GT.1.) THEN PARI(LST(23))=PARI(LST(23))*PARI(32) WRITE(6,1300) PARI(32),INT(PARI(30)+1),PARI(LST(23)), & X,Y,Q2,W2 PARI(32)=1. ENDIF ENDIF IF(VIOL.LT.RLU(0)) GOTO 100 PARL(24)=PARI(31)*PARI(27)/PARI(28) ENDIF IF(ABS(PARL(6)).LT.0.99) THEN C...Choose helicity of incoming lepton. IH=1 IF(RLU(0)*SIGMA.GT.SIGL) IH=2 ENDIF LST(30)=SIGN(1.,IH-1.5) C...Choose target nucleon, proton or neutron. LST(22)=1 IF(PARI(11).GT.1.E-06) THEN IF(RLU(0).LT.(PARI(11)*(PQH(17,IH)-PNT(1,IH)-PNT(2,IH))+ & PNT(2,IH))/PQH(17,IH)) LST(22)=2 ENDIF K(2,2)=40+LST(22) RETURN 1300 FORMAT(' Warning: maximum violated by a factor ',F7.3, &' in event ',I7,/,' maximum increased by this factor to ',E12.3, &/,' Point of violation: x, y, Q**2, W**2 = ',4G10.3) END C ********************************************************************** FUNCTION LKINEM(L) C...Calculate kinematical variables and reject (optionally) if outside C...required limits. COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT1/ MST(40),PAR(80) COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR DOUBLE PRECISION DTHETA,DPHI,DBETA,DE,DPZ,DPT LKINEM=1 IF(L.EQ.-3) THEN C...x,W known from LWEITS, no cuts applied. U=(W2-P(2,5)**2)/(2.*P(2,5)*(1.-X)) Q2=2.*P(2,5)*U*X Y=Q2/(PARL(21)*X) GOTO 200 ENDIF C...x,y given. PARL(22)=Y*PARL(21) Q2=X*PARL(22) U=PARL(22)/(2.*P(2,5)) W2=PARL(22)*(1.-X)+P(2,5)**2 DE=DBLE(P(1,4))*(1.-DBLE(Y))+DBLE(X)*DBLE(Y)*DBLE(ABS(P(2,3))) P(4,4)=SNGL(DE) DPZ=DE-DBLE(X)*DBLE(Y)*(DBLE(P(2,4))+DBLE(ABS(P(2,3)))) DPT=DE*DE-DPZ*DPZ IF(SNGL(DPT).LT.0.) RETURN DPT=DSQRT(DPT) IF(L.EQ.3) GOTO 200 IF(X.LT.CUT(1).OR.X.GT.CUT(2)) RETURN IF(Y.LT.CUT(3).OR.Y.GT.CUT(4)) RETURN IF(Q2.LT.CUT(5).OR.Q2.GT.CUT(6)) RETURN IF(W2.LT.CUT(7).OR.W2.GT.CUT(8)) RETURN IF(U.LT.CUT(9).OR.U.GT.CUT(10)) RETURN C...Transform scattered lepton back to lab system to make cut C...in energy and angle (defined as space angle to incoming lepton). K(6,1)=0 K(6,2)=KSAVE(4) P(6,1)=SNGL(DPT) P(6,2)=0. P(6,3)=SNGL(DPZ) P(6,4)=SNGL(DE) P(6,5)=ULMASS(0,KSAVE(4)) N=6 MST(1)=6 CALL DUROBO(DTHETA(1),DPHI(1),0.D0,0.D0,0.D0) CALL DUROBO(0.D0,0.D0,DBETA(1,1),DBETA(1,2),DBETA(1,3)) MST(1)=0 N=2 IF(P(6,4).LT.CUT(11).OR.P(6,4).GT.CUT(12)) RETURN THETAL=ACOS((PLAB(1,1)*P(6,1)+PLAB(1,2)*P(6,2)+PLAB(1,3)*P(6,3))/ &SQRT(PLAB(1,1)**2+PLAB(1,2)**2+PLAB(1,3)**2)/ &SQRT(P(6,1)**2+P(6,2)**2+P(6,3)**2)) IF(THETAL.LT.CUT(13).OR.THETAL.GT.CUT(14)) RETURN 200 LKINEM=0 RETURN END C ********************************************************************** SUBROUTINE LQCDP(QG,QQB) C...Probabilities for hard QCD events, for given x and W, are obtained C...by making a linear interpolation of the values on the x-W grid. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LGRID/ NXX,NWW,XX(15),WW(15),PQG(15,15,3),PQQB(15,15,2), &QGMAX(15,15,3),QQBMAX(15,15,2),YCUT(15,15),XTOT(15,15),NP DATA NOUT,NABOVE/2*0/,NWARN/10/ QG=0. QQB=0. W=SQRT(W2) C...If W is very small or x close to 1, set QCD weights to zero. IF(WW(1).LT.6..AND.W.LT.WW(1)) RETURN IF(X.GT.XX(NXX)) RETURN XP=X IF(X.LT.XX(1).OR.X.GT.XX(NXX).OR. &W.LT.WW(1).OR.W.GT.WW(NWW)) THEN C...x and/or W ourside limits of grid, write warning for C...first NWARN cases. NOUT=NOUT+1 IF(NOUT.LE.NWARN) WRITE(6,1000) X,W,NWARN IF(X.LT.XX(1)) XP=XX(1) IF(X.GT.XX(NXX)) XP=XX(NXX) IF(W.LT.WW(1)) W=WW(1) IF(W.GT.WW(NWW)) W=WW(NWW) ENDIF IH=1 IF(LST(30).EQ.1) IH=2 IX=0 100 IX=IX+1 IF(XP.GT.XX(IX+1)) GOTO 100 IW=0 200 IW=IW+1 IF(W.GT.WW(IW+1)) GOTO 200 WD=(W-WW(IW))/(WW(IW+1)-WW(IW)) XD=(XP-XX(IX))/(XX(IX+1)-XX(IX)) DO 500 IP=1,NP X1P=(PQG(IX+1,IW,IP)-PQG(IX,IW,IP))*XD+PQG(IX,IW,IP) X2P=(PQG(IX+1,IW+1,IP)-PQG(IX,IW+1,IP))*XD+PQG(IX,IW+1,IP) QGIP=(X2P-X1P)*WD+X1P IF(NP.EQ.1) THEN QG=QGIP PARI(15)=MAX(QGMAX(IX,IW,IH),QGMAX(IX+1,IW+1,IH), & QGMAX(IX+1,IW,IH),QGMAX(IX,IW+1,IH)) ELSE QG=QG+PARI(23+IP)*QGIP PARI(14+IP)=MAX(QGMAX(IX,IW,IP),QGMAX(IX+1,IW+1,IP), & QGMAX(IX+1,IW,IP),QGMAX(IX,IW+1,IP)) ENDIF IF(IP.EQ.3) GOTO 500 X1P=(PQQB(IX+1,IW,IP)-PQQB(IX,IW,IP))*XD+PQQB(IX,IW,IP) X2P=(PQQB(IX+1,IW+1,IP)-PQQB(IX,IW+1,IP))*XD+PQQB(IX,IW+1,IP) QQBIP=(X2P-X1P)*WD+X1P IF(NP.EQ.1) THEN QQB=QQBIP PARI(18)=MAX(QQBMAX(IX,IW,IH),QQBMAX(IX+1,IW+1,IH), & QQBMAX(IX+1,IW,IH),QQBMAX(IX,IW+1,IH)) ELSE QQB=QQB+PARI(23+IP)*QQBIP PARI(17+IP)=MAX(QQBMAX(IX,IW,IP),QQBMAX(IX+1,IW+1,IP), & QQBMAX(IX+1,IW,IP),QQBMAX(IX,IW+1,IP)) ENDIF 500 CONTINUE IF(NP.NE.1) THEN C...Get total x-section from interpolation to be used for normalization. X1P=(XTOT(IX+1,IW)-XTOT(IX,IW))*XD+XTOT(IX,IW) X2P=(XTOT(IX+1,IW+1)-XTOT(IX,IW+1))*XD+XTOT(IX,IW+1) PQ17=(X2P-X1P)*WD+X1P QG=QG/PQ17 QQB=QQB/PQ17 ENDIF C...Get value of y-cut, ie minimum scaled invariant mass squared. PARL(27)=MAX(YCUT(IX,IW),YCUT(IX+1,IW+1), &YCUT(IX+1,IW),YCUT(IX,IW+1)) C...Include alpha-strong in weight. QG=QG*PARL(25) QQB=QQB*PARL(25) IF(QG+QQB.GT.1) THEN C...Sum of QCD event probabilities larger than unity, rescale to unity C...and print warning for first NWARN cases. NABOVE=NABOVE+1 IF(NABOVE.LE.NWARN) WRITE(6,1100) QG,QQB,X,W,NWARN QGQQB=QG+QQB QG=QG/QGQQB QQB=QQB/QGQQB ENDIF 1000 FORMAT(' WARNING: x (=',F7.4,') or W (=',F6.1,') outside grid,', &' weight on limit of grid used.',/,' Note, only first',I5, &' warnings printed.',/) 1100 FORMAT(' WARNING: Sum of QCD probabilities larger than unity!', &' QG, QQB =',2F8.4,' at x, W =',F8.4,F6.1,/, &' Rescaled to unit sum. Note, only first',I5, &' warnings printed.',/) RETURN END C ********************************************************************** SUBROUTINE LQEV C...Generate an ordinary 2-jet event, q-event. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT1/ MST(40),PAR(80) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) LST(24)=1 W=SQRT(W2) C...Choose flavour of scattered quark and target remnant. 200 CALL LFLAV(IFL,IFLR) IF(LST(21).NE.0) GOTO 200 IF(LST(14).EQ.0.OR.IFLR.GT.10.OR.LST(8).GE.2) THEN C...Check if energy in jet system is enough for fragmentation. IF(W.LT.ULMASS(2,IFL)+ULMASS(2,IFLR)+PAR(22)) GOTO 200 CALL LU2JET(MST(1),IFL,IFLR,W) ELSE C...Target remnant is not a simple diquark, special treatment needed. AMIFL=ULMASS(2,IFL) IF(W.LT.AMIFL+ULMASS(2,IFLR)+0.9+PAR(22)) GOTO 200 IFLRO=IFLR NREMH=0 300 NREMH=NREMH+1 IF(NREMH.GT.100) GOTO 999 CALL LREMH(IFLRO,IFLR,K2,XT) AMK2=ULMASS(1,K2) AMIFLR=ULMASS(2,IFLR) C...Give balancing pt to IFLQ and IFLQQ. CALL LPRIKT(PARL(14),PT,PHI) PT2=PT**2 TM2K2=AMK2**2+PT2 EK2=.5*(XT*W+TM2K2/XT/W) PZK2=-.5*(XT*W-TM2K2/XT/W) EPZ=W-TM2K2/XT/W WT=(1.-XT)*W*EPZ-PT2 C...Check if energy in jet system is enough for fragmentation. IF(WT.LT.(AMIFL+AMIFLR+PAR(22))**2) GOTO 300 WT=SQRT(WT+PT2) TMIFLR=AMIFLR**2+PT2 EIFL=.5*(WT+(AMIFL**2-TMIFLR)/WT) EIFLR=.5*(WT+(-AMIFL**2+TMIFLR)/WT) THER=ULANGL(-SQRT(EIFLR**2-TMIFLR),PT) C...Form jet system. CALL LU1JET(-MST(1),IFL,0,0,EIFL,0.,0.) CALL LU1JET(MST(1)+1,IFLR,0,0,EIFLR,THER,PHI) CALL DUROBO(0.D0,0.D0,0.D0,0.D0, & (DBLE(EPZ)-(1.D0-DBLE(XT))*DBLE(W))/ & (DBLE(EPZ)+(1.D0-DBLE(XT))*DBLE(W))) THEK2=ULANGL(PZK2,PT) C...Add formed "target" particle. MST(9)=1 P(MST(1)+2,5)=AMK2 CALL LUPART(MST(1)+2,K2,EK2,THEK2,PHI+3.1415927) MST(9)=0 ENDIF LST(21)=0 RETURN 999 LST(21)=1 RETURN END C ********************************************************************** SUBROUTINE LQGEV C...Generate quark-gluon jet event, choose xp and zp according to QCD C...matrix elements and apply cuts for soft and collinear gluons. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT1/ MST(40),PAR(80) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) LST(24)=2 W=SQRT(W2) J1=MST(1) J2=MST(1)+1 J3=MST(1)+2 J4=MST(1)+3 CALL LXP(XP,IFAIL) IF(IFAIL.NE.0) GOTO 999 C...Choose flavour of scattered quark and target remnant. 200 CALL LFLAV(IFL,IFLR) IF(LST(21).NE.0) RETURN CALL LZP(XP,ZP,IFAIL) IF(IFAIL.NE.0) GOTO 999 AMIFL=ULMASS(2,IFL) AMIFLR=ULMASS(2,IFLR) C CALL HF2(88,XP,ZP,1.) IF(LST(14).EQ.0.OR.IFLR.GT.10) THEN IF(W.LT.AMIFL+AMIFLR+PAR(22)) GOTO 999 IF(LQMCUT(XP,ZP,AMIFL,0.,AMIFLR).NE.0) GOTO 999 CALL LU3JET(J1,IFL,IFLR,W,PARI(21),PARI(23)) CALL DUROBO(DBLE(ACOS(-P(J3,3)/SQRT(P(J3,3)**2+P(J3,1)**2))), & 0.D0,0.D0,0.D0,0.D0) ELSE C...Target remnant is not a simple diquark, special treatment needed. IF(W.LT.AMIFL+AMIFLR+1.+PAR(22)) GOTO 999 IF(LQMCUT(XP,ZP,AMIFL,0.,1.).NE.0) GOTO 999 IFLRO=IFLR NREMH=0 300 NREMH=NREMH+1 IF(NREMH.GT.100) GOTO 999 CALL LREMH(IFLRO,IFLR,K2,XT) AMK2=ULMASS(1,K2) AMIFLR=ULMASS(2,IFLR) P(J1,5)=AMIFL P(J2,5)=0. CALL LPRIKT(PARL(14),PT,PHI) PT2=PT**2 TM2K2=AMK2**2+PT2 TMIFLR=AMIFLR**2+PT2 P(J3,5)=SQRT(TM2K2/XT+TMIFLR/(1.-XT)) IF(LQMCUT(XP,ZP,AMIFL,0.,P(J3,5)).NE.0) GOTO 300 MST(9)=1 CALL LU3JET(J1,IFL,IFLR,W,PARI(21),PARI(23)) MST(9)=0 CALL DUROBO(DBLE(ACOS(-P(J3,3)/SQRT(P(J3,3)**2+P(J3,1)**2))), & 0.D0,0.D0,0.D0,0.D0) EPZ=P(J3,4)-P(J3,3) P(J3,1)=PT*COS(PHI) P(J3,2)=PT*SIN(PHI) P(J3,3)=-0.5*((1.-XT)*EPZ-TMIFLR/(1.-XT)/EPZ) P(J3,4)= 0.5*((1.-XT)*EPZ+TMIFLR/(1.-XT)/EPZ) P(J3,5)=AMIFLR P(J4,1)=-P(J3,1) P(J4,2)=-P(J3,2) P(J4,3)=-0.5*(XT*EPZ-TM2K2/XT/EPZ) P(J4,4)= 0.5*(XT*EPZ+TM2K2/XT/EPZ) P(J4,5)=AMK2 K(J4,1)=0 K(J4,2)=K2 N=J4 IF((P(J3,4)+P(J2,4)/2.)**2-(P(J3,1)+P(J2,1)/2.)**2-P(J3,2)**2 & -(P(J3,3)+P(J2,3)/2.)**2.LT.(AMIFLR+2.5*PAR(22))**2) GOTO 300 ENDIF CALL LAZIMU(XP,ZP) C CALL HF2(89,XP,ZP,1.) LST(21)=0 RETURN 999 LST(21)=1 RETURN END C ********************************************************************** SUBROUTINE LQQBEV C...Generate boson-gluon fusion event, choose xp and zp according to C...QCD matrix elements and apply cuts for softness and collinearness. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT1/ MST(40),PAR(80) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) LST(24)=3 W=SQRT(W2) J1=MST(1) J2=MST(1)+1 J3=MST(1)+2 J4=MST(1)+3 CALL LXP(XP,IFAIL) IF(IFAIL.NE.0) GOTO 999 C...Choose flavour of produced quark-antiquark pair. 200 CALL LFLAV(IFL1,IFL3) IF(LST(21).NE.0) RETURN CALL LZP(XP,ZP,IFAIL) IF(IFAIL.NE.0) GOTO 999 IFL1A=IABS(IFL1) IFL3A=IABS(IFL3) AMIFL1=ULMASS(2,IFL1) AMIFL3=ULMASS(2,IFL3) C CALL HF2(90,XP,ZP,1.) IF(LST(14).EQ.0) THEN C...If baryon production from target remnant is neglected the C...target remnant is approximated by a gluon. IF(W.LT.AMIFL1+AMIFL3+PAR(22)) GOTO 999 IF(LQMCUT(XP,ZP,AMIFL1,0.,AMIFL3).NE.0) GOTO 999 CALL LU3JET(J1,IFL1,IFL3,W,PARI(21),PARI(23)) CALL DUROBO(DBLE(-ACOS(-P(J2,3)/SQRT(P(J2,3)**2+P(J2,1)**2))), & 0.D0,0.D0,0.D0,0.D0) GOTO 400 ENDIF IF(W.LT.AMIFL1+AMIFL3+0.9+2.*PAR(22)) GOTO 999 IF(LQMCUT(XP,ZP,AMIFL1,1.,AMIFL3).NE.0) GOTO 999 P(J1,5)=AMIFL1 P(J3,5)=AMIFL3 C...Choose target valence quark/diquark to form jet system with C...produced antiquark/quark. IFLR2=INT(1.+LST(22)/3.+RLU(0)) IF(IFLR2.EQ.LST(22)) THEN IFLR1=12 IF(RLU(0).GT.PARL(4)) IFLR1=21 ELSE IFLR1=10*LST(22)+LST(22) ENDIF AMR1=ULMASS(2,IFLR1) AMR2=ULMASS(2,IFLR2) NREMH=0 310 NREMH=NREMH+1 IF(NREMH.GT.100) GOTO 999 XT=1.-SQRT(RLU(0)) CALL LPRIKT(PARL(14),PT,PHI) PT2=PT**2 TM2R1=AMR1**2+PT2 TM2R2=AMR2**2+PT2 P(J2,5)=SQRT(TM2R1/(1.-XT)+TM2R2/XT) IF(LQMCUT(XP,ZP,AMIFL1,P(J2,5),AMIFL3).NE.0) GOTO 310 MST(9)=1 CALL LU3JET(J1,IFL1,IFL3,W,PARI(21),PARI(23)) MST(9)=0 CALL DUROBO(DBLE(-ACOS(-P(J2,3)/SQRT(P(J2,3)**2+P(J2,1)**2))), &0.D0,0.D0,0.D0,0.D0) EPZ=P(J2,4)-P(J2,3) IF(IFL1.GT.0) THEN IR1=J2 IR2=J4 ELSE IR1=J4 IR2=J2 ENDIF P(IR1,1)=PT*COS(PHI) P(IR1,2)=PT*SIN(PHI) P(IR1,3)=-0.5*((1.-XT)*EPZ-TM2R1/(1.-XT)/EPZ) P(IR1,4)= 0.5*((1.-XT)*EPZ+TM2R1/(1.-XT)/EPZ) P(IR1,5)=AMR1 P(IR2,1)=-P(IR1,1) P(IR2,2)=-P(IR1,2) P(IR2,3)=-0.5*(XT*EPZ-TM2R2/XT/EPZ) P(IR2,4)= 0.5*(XT*EPZ+TM2R2/XT/EPZ) P(IR2,5)=AMR2 K(IR1,1)=0 K(IR1,2)=IFLR1+500 K(IR2,1)=0 K(IR2,2)=IFLR2+500 K(J3,1)=10000 N=J4 IF((P(J1,4)+P(J2,4))**2-(P(J1,1)+P(J2,1))**2-(P(J1,3)+P(J2,3))**2 & -P(J2,2)**2.LT.(P(J1,5)+P(J2,5)+PAR(22))**2) GOTO 310 IF((P(J3,4)+P(J4,4))**2-(P(J3,1)+P(J4,1))**2-(P(J3,3)+P(J4,3))**2 & -P(J4,2)**2.LT.(P(J3,5)+P(J4,5)+PAR(22))**2) GOTO 310 400 CONTINUE CALL LAZIMU(XP,ZP) C CALL HF2(91,XP,ZP,1.) LST(21)=0 RETURN 999 LST(21)=1 RETURN END C ********************************************************************** SUBROUTINE LXP(XP,IFAIL) C...Choose value of XP according to QCD matrix elements weighted by C...structure functions. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) DOUBLE PRECISION DXPMAX IFAIL=1 XPMIN=DBLE(X)/(1.D0-2.D0*(1.D0-DBLE(X))*DBLE(PARL(27))) DXPMAX=DBLE(X)/(DBLE(X)+(1.D0-DBLE(X))*DBLE(PARL(27))) XPMAX=SNGL(DXPMAX) IF(XPMIN.GE.XPMAX) RETURN AP=1.-XPMIN BP=(1.D0-DXPMAX)/AP IF(LST(24).EQ.2) THEN QXPMAX=PARI(15) IF(LST(16).NE.0) QXPMAX=PARI(24)*PARI(15)+PARI(25)*PARI(16)+ & PARI(26)*PARI(17) ELSE QXPMAX=PARI(18) IF(LST(16).NE.0) QXPMAX=PARI(24)*PARI(18)+PARI(25)*PARI(19) ENDIF C...Safety factor on max value. QXPMAX=QXPMAX*1.05 LOOP=0 100 LOOP=LOOP+1 IF(LOOP.GT.1000) RETURN XP=1.-AP*BP**RLU(0) XPWEIT=DSIGMA(XP)/QXPMAX C CALL HFILL(10,XPWEIT,(XP-XPMIN)/(XPMAX-XPMIN),1.) C CALL HF1(80+(LST(24)-2)*2,XPWEIT,1.) IF(XPWEIT.LT.RLU(0)) GOTO 100 C CALL HF1(81+(LST(24)-2)*2,XPWEIT,1.) IFAIL=0 RETURN END C ********************************************************************** SUBROUTINE LZP(XP,ZP,IFAIL) C...Choose value of ZP according to QCD matrix elements. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) DATA C1,C2/0.2122066,0.0795775/ FQG(DZ,DX,DA,DB,DC)=DA*(DZ**2+DX**2)/(1.-DX)+2.*DA*DX*DZ*(1.-DZ) &+2.*DA*(1.-DZ)+4.*DB*DX*DZ*(1.-DZ)+DC*(DZ**2+DX**2)/(1.-DX)+ &2.*DC*(DX+DZ)*(1.-DZ) FQQ(DZ,DX,DA,DB,DC,DD,DE)=DA*DD*(DZ**2+(1.-DZ)**2)+DB*DE*DZ* &(1.-DZ)+DC*DD*(2.*DZ-1.) IFAIL=1 IH=1 IF(LST(30).EQ.1) IH=2 ZPMIN=(1.-X)*XP/(XP-X)*PARL(27) IF(ZPMIN.GE.0.5) RETURN ZPMAX=1.-ZPMIN I=IABS(LST(25)) AP=1.-ZPMIN BP=ZPMIN/AP IF(LST(1).EQ.2) THEN A=PARI(24) B=PARI(25) CSIGN=-LST(30)*ISIGN(1,LST(25))*PARI(26) ELSE A=(EWQC(1,IH,I)+EWQC(2,IH,I))*PARI(24) B=(EWQC(1,IH,I)+EWQC(2,IH,I))*PARI(25) C=(EWQC(1,IH,I)-EWQC(2,IH,I))*PARI(26) CSIGN=-C*LST(30)*ISIGN(1,LST(25)) ENDIF IF(LST(24).EQ.2) THEN DZPMAX=MAX(FQG(ZPMIN,XP,A,B,CSIGN),FQG(ZPMAX,XP,A,B,CSIGN)) AA=2.*(A+CSIGN)/(1.-XP)-4.*A*XP-8.*B*XP-4.*CSIGN IF(ABS(AA).GT.1.E-20) THEN BB=2.*A*(XP-1.)+4.*B*XP+2.*CSIGN*(1.-XP) Z1=-BB/AA IF(Z1.GT.ZPMIN.AND.Z1.LT.ZPMAX) THEN DZPMAX=MAX(DZPMAX,FQG(Z1,XP,A,B,CSIGN)) ENDIF ENDIF DZPMAX=DZPMAX*C1*1.05 ELSEIF(LST(24).EQ.3) THEN CP=1./BP**2 D=XP**2+(1.-XP)**2 E=8.*XP*(1-XP) DZPMAX=MAX(FQQ(ZPMIN,XP,A,B,CSIGN,D,E), & FQQ(ZPMAX,XP,A,B,CSIGN,D,E)) AA=4.*A*D-2.*B*E IF(ABS(AA).GT.1.E-20) THEN BB=B*E-2.*A*D+2.*CSIGN*D Z1=-BB/AA IF(Z1.GT.ZPMIN.AND.Z1.LT.ZPMAX) THEN DZPMAX=MAX(DZPMAX,FQQ(Z1,XP,A,B,CSIGN,D,E)) ENDIF ENDIF DZPMAX=DZPMAX*C2*1.05 ENDIF IPART=LST(24)-1 LOOP=0 100 LOOP=LOOP+1 IF(LOOP.GT.1000) RETURN IF(LST(24).EQ.2) THEN ZP=1.-AP*BP**RLU(0) SZP=1.-ZP ELSEIF(LST(24).EQ.3) THEN DP=BP*CP**RLU(0) ZP=DP/(1.+DP) SZP=ZP*(1.-ZP) ENDIF ZPWEIT=SZP*(A*DQCD(0,IPART,1,XP,ZP,0.)+B*DQCD(0,IPART,2,XP,ZP,0.) &+CSIGN*DQCD(0,IPART,3,XP,ZP,0.))/DZPMAX C CALL HF1(84+(LST(24)-2)*2,ZPWEIT,1.) IF(ZPWEIT.LT.RLU(0)) GOTO 100 C CALL HF1(85+(LST(24)-2)*2,ZPWEIT,1.) IFAIL=0 RETURN END C ********************************************************************** SUBROUTINE LAZIMU(XP,ZP) C...Choose azimuthal angle (PHI) according to QCD matrix elements. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) J=LST(24)-1 SGN=SIGN(1.,2.5-LST(24)) IFL=LST(25) I=IABS(IFL) IH=1 IF(LST(30).EQ.1) IH=2 IF(LST(1).EQ.2) THEN A=PARI(24)*DQCD(0,J,1,XP,ZP,Y)+PARI(25)*DQCD(0,J,2,XP,ZP,Y) & -LST(30)*ISIGN(1,IFL)*PARI(26)*DQCD(0,J,3,XP,ZP,Y) B=DQCD(1,J,1,XP,ZP,Y) & +SGN*LST(30)*ISIGN(1,IFL)*DQCD(1,J,3,XP,ZP,Y) C=DQCD(2,J,1,XP,ZP,Y) ELSE A=(EWQC(1,IH,I)+EWQC(2,IH,I))*(PARI(24)*DQCD(0,J,1,XP,ZP,Y)+ & PARI(25)*DQCD(0,J,2,XP,ZP,Y)) & -LST(30)*ISIGN(1,IFL)*(EWQC(1,IH,I)-EWQC(2,IH,I)) & *PARI(26)*DQCD(0,J,3,XP,ZP,Y) B=(EWQC(1,IH,I)+EWQC(2,IH,I))*DQCD(1,J,1,XP,ZP,Y) & +SGN*LST(30)*ISIGN(1,IFL)*(EWQC(1,IH,I)-EWQC(2,IH,I)) & *DQCD(1,J,3,XP,ZP,Y) C=(EWQC(1,IH,I)+EWQC(2,IH,I))*DQCD(2,J,1,XP,ZP,Y) ENDIF PHIMAX=ABS(A)+ABS(B)+ABS(C) 100 PHI=6.2832*RLU(0) IF(A+B*COS(PHI)+C*COS(2.*PHI).LT.RLU(0)*PHIMAX) GOTO 100 CALL DUROBO(0.D0,DBLE(PHI),0.D0,0.D0,0.D0) C CALL HF1(95+J,PHI*180./3.1415927,1.) RETURN END C ********************************************************************** FUNCTION DSIGMA(XP) C...Differential cross section for first order QCD processes. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /LUJETS/ N,K(2000,2),P(2000,5) DIMENSION XPQ(-6:6),PQH(17,2) DSIGMA=0. DO 10 I=1,17 PQH(I,1)=0. PQH(I,2)=0. 10 PQ(I)=0. AMU=ULMASS(2,1) IF(LST(20).EQ.0.OR.LST(16).EQ.0) THEN IL=1 IU=3 IF(LST(1).EQ.1.OR.LST(24).EQ.3) IU=2 ELSE IL=LST(20) IU=LST(20) ENDIF XI=X/XP ZPMIN=(1.-X)*XP/(XP-X)*PARL(27) IF(ZPMIN.GE.0.5) RETURN ZPMAX=1.D0-DBLE(ZPMIN) CALL LNUCL(XI,Q2,XPQ) IF(LST(24).EQ.3) GOTO 3000 C...Gluon bremsstrahlung process, i.e. qg-event. 2000 DO 2500 IP=IL,IU SIG=DQCDI(1,IP,XP,ZPMIN,ZPMAX) SGN=SIGN(1.,5.-2.*IP) DO 2300 IH=1,2 IF(IH.EQ.1) THEN IF(PARL(6).GT.0.99) GOTO 2300 IF(LST(20).EQ.0.AND.LST(30).NE.-1) GOTO 2300 ELSEIF(IH.EQ.2) THEN IF(PARL(6).LT.-0.99) GOTO 2300 IF(LST(20).EQ.0.AND.LST(30).NE.1) GOTO 2300 ENDIF IF(LST(20).NE.0) LST(30)=SIGN(1.,IH-1.5) IF(LST(1).NE.2) THEN DO 2100 I=1,LST(12) WQ=XPQ(I)*SIG*(EWQC(1,IH,I)+SGN*EWQC(2,IH,I)) WQB=XPQ(-I)*SIG*SGN*(EWQC(1,IH,I)+SGN*EWQC(2,IH,I)) C...Include y-dependence. WQ=WQ*PARI(23+IP) WQB=WQB*PARI(23+IP) PQH(I,IH)=PQH(I,IH)+WQ PQH(I+LST(12),IH)=PQH(I+LST(12),IH)+WQB PQH(17,IH)=PQH(17,IH)+WQ+WQB 2100 CONTINUE ELSEIF(LST(1).EQ.2) THEN C...Zero CC cross-section for one helicity state. IF(KSAVE(1).LT.0.AND.IH.EQ.1 & .OR.KSAVE(1).GT.0.AND.IH.EQ.2) GOTO 2300 IF(IP.EQ.3) THEN TQ=-LST(30) TQB=-TQ ELSE TQ=1. TQB=1. ENDIF DO 2200 I=1,LST(12) T1=-K(3,2)*QC(I) IF(T1.GT.0) THEN WQ=XPQ(I)*SIG*TQ WQB=0. ELSE WQB=XPQ(-I)*SIG*TQB WQ=0. ENDIF C...Include y-dependence. WQ=WQ*PARI(23+IP) WQB=WQB*PARI(23+IP) PQH(I,IH)=PQH(I,IH)+WQ PQH(I+LST(12),IH)=PQH(I+LST(12),IH)+WQB PQH(17,IH)=PQH(17,IH)+WQ+WQB 2200 CONTINUE ENDIF 2300 CONTINUE 2500 CONTINUE DO 2600 I=1,17 2600 PQ(I)=(1.-PARL(6))/2.*PQH(I,1)+(1.+PARL(6))/2.*PQH(I,2) IF(LST(20).EQ.0) THEN C...Simulation: cross section for chosen helicity state only. IH=1 IF(LST(30).EQ.1) IH=2 DSIGMA=PQH(17,IH) ELSE C...Integration: normalize and include alpha_s*1/(1-xp) DSIGMA=PQ(17)/PARI(20)*PARL(25)/(1.-XP) IF(LST(16).EQ.0) THEN C...Fixed beam energy, max of dsigma/dxp for L- and R-handed lepton. IF(PQH(17,1).GT.PARI(15)) PARI(15)=PQH(17,1) IF(PQH(17,2).GT.PARI(16)) PARI(16)=PQH(17,2) ELSE C...Variable beam energy, max of dsigma/dxp for S, T, I contributions. IF(PQ(17)/PARI(23+LST(20)).GT.PARI(14+LST(20))) & PARI(14+LST(20))=PQ(17)/PARI(23+LST(20)) ENDIF ENDIF RETURN C...Boson-gluon fusion, i.e. qq-event. 3000 S13=Q2*(1.-XP)/XP IF(S13.LT.4.*AMU**2) RETURN DO 3500 IP=IL,IU SIG=XPQ(0)*DQCDI(2,IP,XP,ZPMIN,ZPMAX) DO 3300 IH=1,2 IF(IH.EQ.1) THEN IF(PARL(6).GT.0.99) GOTO 3300 IF(LST(20).EQ.0.AND.LST(30).NE.-1) GOTO 3300 ELSEIF(IH.EQ.2) THEN IF(PARL(6).LT.-0.99) GOTO 3300 IF(LST(20).EQ.0.AND.LST(30).NE.1) GOTO 3300 ENDIF IF(LST(20).NE.0) LST(30)=SIGN(1.,IH-1.5) IF(LST(1).NE.2) THEN DO 3100 I=1,LST(13) IF(S13.LT.4.*ULMASS(2,I)**2) GOTO 3100 WQ=SIG/2.*(EWQC(1,IH,I)+EWQC(2,IH,I)) WQB=WQ C...Include y-dependence. WQ=WQ*PARI(23+IP) WQB=WQB*PARI(23+IP) PQH(I,IH)=PQH(I,IH)+WQ PQH(I+LST(13),IH)=PQH(I+LST(13),IH)+WQB PQH(17,IH)=PQH(17,IH)+WQ+WQB 3100 CONTINUE ELSEIF(LST(1).EQ.2) THEN C...Zero CC cross-section for one helicity state. IF(KSAVE(1).LT.0.AND.IH.EQ.1 & .OR.KSAVE(1).GT.0.AND.IH.EQ.2) GOTO 3300 DO 3200 I=1,LST(13) IF(S13.LT.(AMU+ULMASS(2,I))**2) GOTO 3200 IF(K(3,2)*QC(I).LT.0) THEN WQ=SIG WQB=0. ELSE WQB=SIG WQ=0. ENDIF C...Include y-dependence. WQ=WQ*PARI(23+IP) WQB=WQB*PARI(23+IP) PQH(I,IH)=PQH(I,IH)+WQ PQH(I+LST(13),IH)=PQH(I+LST(13),IH)+WQB PQH(17,IH)=PQH(17,IH)+WQ+WQB 3200 CONTINUE ENDIF 3300 CONTINUE 3500 CONTINUE DO 3600 I=1,17 3600 PQ(I)=(1.-PARL(6))/2.*PQH(I,1)+(1.+PARL(6))/2.*PQH(I,2) IF(LST(20).EQ.0) THEN C...Simulation: cross section for chosen helicity state only. IH=1 IF(LST(30).EQ.1) IH=2 DSIGMA=PQH(17,IH) ELSE C...Integration: normalize and include alpha_s*1/(1-xp) DSIGMA=PQ(17)/PARI(20)*PARL(25)/(1.-XP) IF(LST(16).EQ.0) THEN C...Fixed beam energy, max of dsigma/dxp for L- and R-handed lepton. IF(PQH(17,1).GT.PARI(18)) PARI(18)=PQH(17,1) IF(PQH(17,2).GT.PARI(19)) PARI(19)=PQH(17,2) ELSE C...Variable beam energy, max of dsigma/dxp for S, T, I contributions. IF(PQ(17)/PARI(23+LST(20)).GT.PARI(17+LST(20))) & PARI(17+LST(20))=PQ(17)/PARI(23+LST(20)) ENDIF ENDIF RETURN END C ********************************************************************** FUNCTION LQMCUT(XP,ZP,AM1,AM2,AM3) C...Apply cuts, if necessary, on the event configuration C...obtained from QCD matrix elements. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) IF(LST(24).EQ.2) THEN S12=Q2*(1.-XP)/XP S23=Q2*(XP-X)*(1.-ZP)/X/XP+AM2**2+AM3**2 S13=Q2*(XP-X)*ZP/X/XP+AM1**2+AM3**2 ELSEIF(LST(24).EQ.3) THEN S13=Q2*(1.-XP)/XP S23=Q2*(XP-X)*(1.-ZP)/X/XP+AM2**2+AM3**2 S12=Q2*(XP-X)*ZP/X/XP+AM1**2+AM2**2 IF(S13.LT.(AM1+AM3)**2) GOTO 900 ENDIF W=SQRT(W2) X1=1.-(S23-AM1**2)/W2 X3=1.-(S12-AM3**2)/W2 X2=2.-X1-X3 C CALL HF2(92+(LST(24)-2)*2,X1,X3,1.) PARI(21)=X1 PARI(22)=X2 PARI(23)=X3 IF(X1.GT.1..OR.X2.GT.1..OR.X3.GT.1.) GOTO 900 IF(X1*W/2..LT.AM1.OR.X2*W/2..LT.AM2.OR.X3*W/2..LT.AM3) GOTO 900 PA1=SQRT((0.5*X1*W)**2-AM1**2) PA2=SQRT((0.5*X2*W)**2-AM2**2) PA3=SQRT((0.5*X3*W)**2-AM3**2) IF(ABS((PA3**2-PA1**2-PA2**2)/(2.*PA1*PA2)).GE.1.) GOTO 900 IF(ABS((PA2**2-PA1**2-PA3**2)/(2.*PA1*PA3)).GE.1.) GOTO 900 C CALL HF2(93+(LST(24)-2)*2,X1,X3,1.) LQMCUT=0 RETURN 900 LQMCUT=1 RETURN END C ********************************************************************** FUNCTION DQCD(ICOSFI,IPART,IP,XP,ZP,Y) C...First order QCD matrix elements from R.D. Peccei and R. Ruckl: C...Nucl. Phys. B162 (1980) 125 C...Constants C1 to C5 are resp. 2/3/pi, 1/4/pi, 4/3/pi, 1/2/pi, 1/pi DATA C1,C2,C3,C4,C5/0.2122066,0.0795775,0.4244132,0.1591549, & 0.3183099/ IF(ICOSFI.EQ.0) THEN IF(IPART.EQ.1) THEN IF(IP.EQ.1) THEN DQCD=C1*((ZP**2+XP**2)/(1.-XP)/(1.-ZP)+2.*(XP*ZP+1.)) ELSEIF(IP.EQ.2) THEN DQCD=C1*4.*XP*ZP ELSEIF(IP.EQ.3) THEN DQCD=C1*((ZP**2+XP**2)/(1.-XP)/(1.-ZP)+2.*(XP+ZP)) ELSE WRITE(6,1000) ICOSFI,IPART,IP ENDIF ELSEIF(IPART.EQ.2) THEN IF(IP.EQ.1) THEN DQCD=C2*(XP**2+(1.-XP)**2)*(ZP**2+(1.-ZP)**2)/(1.-ZP)/ZP ELSEIF(IP.EQ.2) THEN DQCD=C2*8.*XP*(1.-XP) ELSEIF(IP.EQ.3) THEN DQCD=C2*(XP**2+(1.-XP)**2)*(ZP-(1.-ZP))/(1.-ZP)/ZP ELSE WRITE(6,1000) ICOSFI,IPART,IP ENDIF ELSE WRITE(6,1000) ICOSFI,IPART,IP ENDIF ELSEIF(ICOSFI.EQ.1) THEN IF(IPART.EQ.1) THEN IF(IP.EQ.1) THEN DQCD=C3*Y*SQRT((1.-Y)*XP*ZP/(1.-XP)/(1.-ZP))* & (1.-2./Y)*(1.-ZP-XP+2.*XP*ZP) ELSEIF(IP.EQ.3) THEN DQCD=C3*Y*SQRT((1.-Y)*XP*ZP/(1.-XP)/(1.-ZP))* & (1.-XP-ZP) ELSE WRITE(6,1000) ICOSFI,IPART,IP ENDIF ELSEIF(IPART.EQ.2) THEN IF(IP.EQ.1) THEN DQCD=C4*Y*SQRT((1.-Y)*XP*(1.-XP)/ZP/(1.-ZP))* & (1.-2./Y)*(1.-2.*ZP)*(1.-2.*XP) ELSEIF(IP.EQ.3) THEN DQCD=C4*Y*SQRT((1.-Y)*XP*(1.-XP)/ZP/(1.-ZP))* & (1.-2.*XP) ELSE WRITE(6,1000) ICOSFI,IPART,IP ENDIF ENDIF ELSEIF(ICOSFI.EQ.2) THEN IF(IPART.EQ.1) THEN DQCD=C3*(1.-Y)*XP*ZP ELSEIF(IPART.EQ.2) THEN DQCD=C5*(1.-Y)*XP*(1.-XP) ELSE WRITE(6,1000) ICOSFI,IPART,IP ENDIF ELSE WRITE(6,1000) ICOSFI,IPART,IP ENDIF RETURN 1000 FORMAT(' Error in routine DQCD !! ', &' ICOSFI, IPART, IP = ',3I10) END C ********************************************************************** FUNCTION DQCDI(IPART,IP,XP,ZPMIN,ZPMAX) C...First order QCD matrix elements as in function DQCD but analytically C...integrated over ZP from ZPMIN to ZPMAX, also a factor 1/(1-XP) is C...factored out (since XP chosen randomly according to 1/(1-XP) distr.) DATA C1,C2/0.2122066,0.0795775/ IF(IPART.EQ.1) THEN IF(IP.EQ.1) THEN ZLOG=ALOG(ZPMAX/ZPMIN) DQCDI=C1*(XP**2*ZLOG+ZPMIN-ZPMAX+(ZPMIN**2-ZPMAX**2)/2.+ZLOG+ & XP*(1.-XP)*(ZPMAX**2-ZPMIN**2)+2.*(1.-XP)*(ZPMAX-ZPMIN)) ELSEIF(IP.EQ.2) THEN DQCDI=C1*2.*XP*(1.-XP)*(ZPMAX**2-ZPMIN**2) ELSEIF(IP.EQ.3) THEN ZLOG=ALOG(ZPMAX/ZPMIN) DQCDI=C1*(XP**2*ZLOG+ZPMIN-ZPMAX+(ZPMIN**2-ZPMAX**2)/2.+ZLOG+ & 2.*XP*(1.-XP)*(ZPMAX-ZPMIN)+(1.-XP)*(ZPMAX**2-ZPMIN**2)) ELSE WRITE(6,1000) IPART,IP ENDIF ELSEIF(IPART.EQ.2) THEN IF(IP.EQ.1) THEN DQCDI=C2*(1.-XP)*(XP**2+(1.-XP)**2)*(2.*(ZPMIN-ZPMAX)+ & ALOG(ZPMAX**2/ZPMIN**2)) ELSEIF(IP.EQ.2) THEN DQCDI=C2*8.*XP*(1.-XP)**2*(ZPMAX-ZPMIN) ELSEIF(IP.EQ.3) THEN DQCDI=0. ELSE WRITE(6,1000) IPART,IP ENDIF ELSE WRITE(6,1000) IPART,IP ENDIF RETURN 1000 FORMAT(' Error in routine DQCDI !! ', &' IPART, IP = ',2I10) END C ******************************************************************** FUNCTION DALPS(Q2A) C...Alpha-strong for hard QCD effects. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) C...Number of active flavours. NF=2 DO 100 IFL=3,LST(9) 100 IF(Q2A.GT.4.*ULMASS(2,IFL)**2) NF=NF+1 DALPS=12.*3.1416/((33.-2.*NF)*ALOG(MAX(Q2A,PARL(7))/PARL(10)**2)) RETURN END C ********************************************************************** SUBROUTINE LFLAV(IFL,IFLR) C...Choose flavour of struck quark and the C...corresponding flavour of the target remnant jet. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) COMMON /LFLMIX/ CABIBO(4,4) LST(21)=0 LST(27)=0 IF(LST(24).EQ.3) THEN NFL=LST(13) ELSE NFL=LST(12) ENDIF 20 R=RLU(0)*PQ(17) PSUB=0. DO 30 I=1,2*NFL IFL=I PSUB=PSUB+PQ(I) IF(R.LE.PSUB) GOTO 40 30 CONTINUE 40 CONTINUE IF(IFL.GT.NFL) IFL=NFL-IFL LST(25)=IFL IFLR=-IFL IF(LST(1).EQ.2) THEN C...Weak charged current, change the flavour of the struck C...quark using generalized Cabibbo mixing matrix. IFLA=IABS(IFL) J1=(IFLA+1)/2 IF(IFLA.LE.2) IFLA=3-IFLA M1=MOD(IFLA,2) M2=MOD(IFLA+1,2) R=RLU(0) PSUB=0. DO 100 J=1,4 J2=J PSUB=PSUB+CABIBO(M1*J2+M2*J1,M2*J2+M1*J1) IF(R.LT.PSUB) GOTO 200 100 CONTINUE 200 IFL=2*J2-M2 IF(J2.EQ.1) IFL=2*MOD(IFL,2)+MOD(IFL+1,2) IF(LST(25).LT.0) IFL=-IFL ENDIF IFLA=IABS(IFL) IFLRA=IABS(IFLR) IF(IFLA.GE.LST(11).OR.IFLRA.GE.LST(11)) THEN C...Threshold function for heavy quarks of flavour IFLA and IFLRA. IF(1.-(.938+PMAS(100+IFLA)+PMAS(100+IFLRA)+2.*PMAS(101))**2/W2 & .LT.RLU(0)) GOTO(20,999,999) LST(24) ENDIF IF(LST(24).EQ.3) RETURN C...With LST(14)=0/1(default) baryon production from the target remnant C...is excluded/included. IF(LST(14).EQ.0) RETURN IF(IFLR.EQ.-1) THEN IF(LST(22).EQ.1) THEN IFLR=12 IF(RLU(0).GT.PARL(4)) IFLR=21 ELSE IFLR=22 ENDIF ELSEIF(IFLR.EQ.-2) THEN IF(LST(22).EQ.1) THEN IFLR=11 ELSE IFLR=12 IF(RLU(0).GT.PARL(4)) IFLR=21 ENDIF ENDIF RETURN 999 LST(21)=1 RETURN END C ********************************************************************** SUBROUTINE LREMH(IFLRO,IFLR,K2,XT) C...Gives flavour and energy-momentum fraction XT for the particle C...to be produced out of the target remnant when that is not a C...simple diquark. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) 100 IFLQ=INT(1.+LST(22)/3.+RLU(0)) IF(IFLQ.EQ.LST(22)) THEN IFLQQ=12 IF(RLU(0).GT.PARL(4)) IFLQQ =21 ELSE IFLQQ=10*LST(22)+LST(22) ENDIF XT=1.-SQRT(RLU(0)) IF(IFLRO.GT.0) THEN CALL LUIFLD(IFLQQ,IFLQQ/10,IFLRO,IDUM,K2) IF(K2.EQ.0) GOTO 100 IFLR=IFLQ XT=1.-XT ELSE CALL LUIFLD(IFLQ,0,IFLRO,IDUM,K2) IFLR=IFLQQ ENDIF LST(27)=1 RETURN END C ********************************************************************** SUBROUTINE LPRIKT(S,PT,PHI) C...Size (PT) and azimuthal angle (PHI) of primordial kt according C...to a Gaussian distribution. PT=S*SQRT(-ALOG(RLU(0))) PHI=6.2832*RLU(0) RETURN END C ********************************************************************** SUBROUTINE LFRAME(IFR,IPH) C...Make transformation from hadronic CM frame to lab frame. COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR DOUBLE PRECISION DTHETA,DPHI,DBETA,DE,DPZ,DPT INTEGER IFR,IPH,IFRAME,IPHI IFRAME=IFR IPHI=IPH IF(IFRAME.LT.1.OR.IFRAME.GT.4.OR.IPHI.LT.0.OR.IPHI.GT.1) GOTO 999 IF(IFRAME.EQ.1) IPHI=0 N=N+1 DO 5 J=1,5 5 P(N,J)=PB(J) 10 CONTINUE IF(IPHI.NE.LST(29)) THEN IFRAME=2 ELSE IFRAME=IFR ENDIF IF((IFRAME.EQ.LST(28)).AND.(IPHI.EQ.LST(29))) THEN DO 15 J=1,5 15 PB(J)=P(N,J) N=N-1 RETURN ENDIF GOTO(100,200,300,400), LST(28) GOTO 999 100 IF(IFRAME.GE.2) THEN CALL DUROBO(DTHETA(2),DPHI(2),0.D0,0.D0,0.D0) CALL DUROBO(0.D0,0.D0,DBETA(2,1),DBETA(2,2),DBETA(2,3)) LST(28)=2 ELSE GOTO 999 ENDIF GOTO 10 200 IF(IPHI.NE.LST(29)) THEN CALL DUROBO(0.D0,DBLE(SIGN(PHIR,FLOAT(IPHI-LST(29)))), & 0.D0,0.D0,0.D0) LST(29)=IPHI ENDIF IF(IFRAME.EQ.1) THEN CALL DUROBO(0.D0,0.D0,-DBETA(2,1),-DBETA(2,2),-DBETA(2,3)) CALL DUROBO(-DTHETA(2),0.D0,0.D0,0.D0,0.D0) LST(28)=1 ELSEIF(IFRAME.GE.3) THEN CALL DUROBO(DTHETA(1),DPHI(1),0.D0,0.D0,0.D0) CALL DUROBO(0.D0,0.D0,DBETA(1,1),DBETA(1,2),DBETA(1,3)) LST(28)=3 ENDIF GOTO 10 300 IF(IFRAME.LE.2) THEN CALL DUROBO(0.D0,0.D0,-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) CALL DUROBO(0.D0,-DPHI(1),0.D0,0.D0,0.D0) CALL DUROBO(-DTHETA(1),0.D0,0.D0,0.D0,0.D0) LST(28)=2 ELSEIF(IFRAME.EQ.4) THEN THEBOS=PLU(N,13) PHIBOS=PLU(N,15) CALL DUROBO(0.D0,DBLE(-PHIBOS),0.D0,0.D0,0.D0) CALL DUROBO(DBLE(-THEBOS),0.D0,0.D0,0.D0,0.D0) LST(28)=4 ENDIF GOTO 10 400 IF(IFRAME.LE.3) THEN CALL DUROBO(DBLE(THEBOS),DBLE(PHIBOS),0.D0,0.D0,0.D0) LST(28)=3 ENDIF GOTO 10 999 WRITE(6,1000) IFRAME,IPHI,LST(28),LST(29) 1000 FORMAT(' BAD VARIABLES IN SUBROUTINE LFRAME: IFRAME,IPHI,', &'LST(28),LST(29) =',4I5) RETURN END C ******************************************************************** SUBROUTINE DUROBO(THE,PHI,BEX,BEY,BEZ) C...Rotation and boost in double precision, modification of LUROBO. COMMON/LUJETS/N,K(2000,2),P(2000,5) COMMON/LUDAT1/MST(40),PAR(80) DOUBLE PRECISION THE,PHI,BEX,BEY,BEZ,DGA,DBEP,DGABEP DOUBLE PRECISION ROT(3,3),PV(3),DP(4) IMAX=N IF(MST(2).GT.0) IMAX=MST(2) IF(SNGL(THE**2+PHI**2).GT.1E-20) THEN C...Rotate (typically from z axis to direction theta,phi) ROT(1,1)=DCOS(THE)*DCOS(PHI) ROT(1,2)=-DSIN(PHI) ROT(1,3)=DSIN(THE)*DCOS(PHI) ROT(2,1)=DCOS(THE)*DSIN(PHI) ROT(2,2)=DCOS(PHI) ROT(2,3)=DSIN(THE)*DSIN(PHI) ROT(3,1)=-DSIN(THE) ROT(3,2)=0. ROT(3,3)=DCOS(THE) DO 120 I=MAX(1,MST(1)),IMAX IF(MOD(K(I,1)/10000,10).GE.6) GOTO 120 DO 100 J=1,3 100 PV(J)=P(I,J) DO 110 J=1,3 110 P(I,J)=ROT(J,1)*PV(1)+ROT(J,2)*PV(2)+ROT(J,3)*PV(3) 120 CONTINUE ENDIF IF(SNGL(BEX**2+BEY**2+BEZ**2).GT.1E-20) THEN C...Lorentz boost (typically from rest to momentum/energy=beta) DGA=1D0/DSQRT(1D0-BEX**2-BEY**2-BEZ**2) DO 140 I=MAX(1,MST(1)),IMAX IF(MOD(K(I,1)/10000,10).GE.6) GOTO 140 DO 130 J=1,4 130 DP(J)=P(I,J) DBEP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3) DGABEP=DGA*(DGA*DBEP/(1D0+DGA)+DP(4)) P(I,1)=DP(1)+DGABEP*BEX P(I,2)=DP(2)+DGABEP*BEY P(I,3)=DP(3)+DGABEP*BEZ P(I,4)=DGA*(DP(4)+DBEP) 140 CONTINUE ENDIF RETURN END C ********************************************************************** SUBROUTINE LWBB(ENU) C...Gives energy (ENU) of a (anti-)neutrino chosen from a simple C...parametrization of a wide band beam. DATA EMEAN,SLOPE,EMIN,EMAX/30.,0.02,12.,300./ A1=1./(EMEAN-12.) A2=EXP(EMEAN*SLOPE) 100 ENU=300.*RLU(0) IF(ENU.LT.EMEAN)THEN E=A1*(ENU-12.) ELSE E=A2*EXP(-ENU*SLOPE) ENDIF IF(ENU.LT.EMIN.OR.ENU.GT.EMAX) GOTO 100 IF(E.LT.RLU(0)) GOTO 100 RETURN END C ********************************************************************** SUBROUTINE LWEITS(LFILE) C...Integrates the QCD matrix elements to obtain probabilities for C...qg- and qq-events as a function of (x,W). Also finds various C...maximum values to be used for the QCD simulation. Results stored C...in common LGRID and optionally written to logical file LFILE. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),COMFAC COMMON /LGRID/ NXX,NWW,XX(15),WW(15),PQG(15,15,3),PQQB(15,15,2), &QGMAX(15,15,3),QQBMAX(15,15,2),YCUT(15,15),XTOT(15,15),NP DIMENSION WWI(15,4),XXI(15,4) DATA WWI/5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,17.5,20.,22.5,25., &5.,7.5,10.,12.5,15.,17.5,20.,22.5,25.,27.5,30.,32.5,35.,40.,45., &5.,10.,20.,30.,50.,75.,100.,125.,150.,175.,200.,225.,250.,300., &350.,5.,10.,20.,35.,60.,100.,150.,225.,350.,500.,700.,1000., &1400.,1900.,2500./ DATA XXI/.01,.02,.04,.06,.08,.1,.125,.15,.2,.25,.3,.45,.6,.75,.99, &.01,.02,.04,.06,.08,.1,.125,.15,.2,.25,.3,.45,.6,.75,.99, &.001,.005,.01,.02,.04,.06,.08,.1,.125,.15,.2,.3,.5,.75,.99, &.001,.005,.01,.02,.04,.06,.08,.1,.125,.15,.2,.3,.5,.75,.99/ EXTERNAL DSIGMA LST2=LST(2) LST(2)=-3 WMAX=SQRT(PARL(21)) WRITE(6,1000) PARL(11),LST(13),LST(9),PARL(10), &PARL(8),PARL(9),PARL(12),PARL(13) IF(LST(16).EQ.0) THEN NP=1 IPMAX=2 WRITE(6,1010) ELSE NP=3 IF(LST(1).EQ.1) NP=2 IPMAX=3 WRITE(6,1020) ENDIF IF(LST(19).GE.1.AND.LST(19).LE.4) THEN C...Grid taken from data in arrays WWI, XXI. DO 10 IW=1,NWW 10 WW(IW)=WWI(IW,LST(19)) DO 20 IX=1,NXX 20 XX(IX)=XXI(IX,LST(19)) ELSE C...Grid specified by user. READ(5,*) NWW,NXX READ(5,*) (WW(IW),IW=1,NWW) READ(5,*) (XX(IX),IX=1,NXX) IF(XX(NXX).GT..99) XX(NXX)=.99 ENDIF WRITE(6,1030) LST(19),NWW,NXX,WW,XX IF(WMAX.GT.WW(NWW)) THEN WRITE(6,1040) WMAX,WW(NWW) STOP ENDIF WRITE(6,1100) LW=0 DO 500 IW=1,NWW W=WW(IW) IF(LW.GT.0) GOTO 600 IF(W.GT.WMAX) LW=LW+1 W2=W**2 LX=0 DO 400 IX=1,NXX X=XX(IX) IF(LX.GT.0) GOTO 500 IF(X.GT.1.-W2/WMAX**2) LX=LX+1 CALL LEPTO PQCOM=PARI(31)*PQ(17)*COMFAC PARL(25)=DALPS(Q2) PARI(20)=PQ(17) XTOT(IX,IW)=PQ(17) PARL(27)=MAX(PARL(9)**2/W2,PARL(8)) IYCUT=0 100 IYCUT=IYCUT+1 RQG=0. RQQB=0. XPMIN=DBLE(X)/(1.D0-2.D0*(1.D0-DBLE(X))*DBLE(PARL(27))) XPMAX=DBLE(X)/(DBLE(X)+(1.D0-DBLE(X))*DBLE(PARL(27))) IF(XPMIN.GE.XPMAX) GOTO 210 DO 200 IP=1,NP IF(LST(16).EQ.0) THEN PARI(15)=0. PARI(16)=0. PARI(18)=0. PARI(19)=0. ELSE PARI(14+IP)=0. IF(IP.LE.2) PARI(17+IP)=0. ENDIF LST(20)=IP LST(24)=2 EPS=PARL(11) CALL GADAP(XPMIN,XPMAX,DSIGMA,EPS,RESULT) RQG=RQG+RESULT PQG(IX,IW,IP)=RESULT/PARL(25) IF(LST(16).EQ.0) THEN QGMAX(IX,IW,1)=PARI(15) QGMAX(IX,IW,2)=PARI(16) ELSE PQG(IX,IW,IP)=RESULT*PARI(20)/PARI(23+IP)/PARL(25) QGMAX(IX,IW,IP)=PARI(14+IP) ENDIF IF(IP.EQ.3) GOTO 200 LST(24)=3 EPS=PARL(11) CALL GADAP(XPMIN,XPMAX,DSIGMA,EPS,RESULT) RQQB=RQQB+RESULT PQQB(IX,IW,IP)=RESULT/PARL(25) IF(LST(16).EQ.0) THEN QQBMAX(IX,IW,1)=PARI(18) QQBMAX(IX,IW,2)=PARI(19) ELSE PQQB(IX,IW,IP)=RESULT*PARI(20)/PARI(23+IP)/PARL(25) QGMAX(IX,IW,IP)=PARI(17+IP) ENDIF 200 CONTINUE 210 CONTINUE RQ=1.-RQG-RQQB IF(RQ.LT.0.) THEN C...QCD probabilities > 1, increase cutoff. YCLOW=PARL(27) POT=SQRT(1./(RQG+RQQB)) PARL(27)=(1./PARL(12)+0.01)*(PARL(12)*PARL(27))**POT GOTO 100 ELSEIF(IYCUT.GT.1.AND.RQ.GT.PARL(13)) THEN C...Cutoff increased too much, try lower. PARL(27)=(PARL(27)+YCLOW)/2. GOTO 100 ENDIF YCUT(IX,IW)=PARL(27) WRITE(6,1200) W,X,Y,Q2,PARL(25),PQCOM,PARL(27),IYCUT, &RQ,RQG,RQQB,(QGMAX(IX,IW,IP),IP=1,IPMAX), &(QQBMAX(IX,IW,IP),IP=1,MIN(2,IPMAX)) 400 CONTINUE 500 CONTINUE 600 CONTINUE LST(2)=LST2 IF(LFILE.LT.0) THEN C...Write results on logical file number IABS(LFILE) WRITE(IABS(LFILE)) LST,PARL,NXX,NWW,NP,XX,WW WRITE(IABS(LFILE))(((PQG(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,NP), & (((PQQB(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,NP), & (((QGMAX(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,IPMAX), & (((QQBMAX(IX,IW,IP),IX=1,NXX),IW=1,NWW),IP=1,MIN(2,IPMAX)), & YCUT IF(NP.NE.1) WRITE(IABS(LFILE)) XTOT CLOSE(IABS(LFILE)) ENDIF RETURN 1000 FORMAT('1',/,5X,'Integration of 1st order QCD matrix elements', & /,5X,'============================================', &/,' for gluon radiation (qg-event) and boson-gluon fusion ', &'(qq-event) probability.', &//,' Required precision in integration, PARL(11) =',F8.4, &//,' Heaviest flavour produced in boson-gluon fusion, LST(13) =', &I5,//,' Alpha-strong parameters: max flavour, LST(9) =',I3, &' QCD lambda, PARL(10) =',F6.3,' GeV', &//,' Cuts on matrix elements:', &/,' PARL(8), PARL(9), PARL(12), PARL(13) =',4F8.4,/) 1010 FORMAT(' Lepton energy not allowed to vary in simulation.',/) 1020 FORMAT(' Lepton energy allowed to vary in simulation, ',/, &' y in table below calculated assuming max energy.',/) 1030 FORMAT(' Grid choice, LST(19) =',I3,5X,'# grid points in W, X =', &2I5,/,' W-values in array WW:',/,15F8.1, &/,' x-values in array XX:',/,15F8.4,/) 1040 FORMAT(' Max W outside grid, execution stopped ! Wmax, grid-max =' &,2F12.1) 1100 FORMAT(//,5X,'W',6X,'x',6X,'y',5X,'Q**2',1X,'alpha',2X,'dsigma', &8X,'cut',' it',2X,'q-event',1X,'qg-event', &1X,'qq-event',' max of matrix elements qg & qq (L,R or T,S,I)', &/,1X,129(1H-),/) 1200 FORMAT(F6.1,2F7.3,F9.1,F6.2,G12.3,F7.3,I3,3F9.4,5E9.2) END C ********************************************************************** SUBROUTINE LPRWTS(NSTEP) C...Prints probabilities for q-, qg- and qqbar-events using the present C...QCD weights stored in common block LGRID. COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LGRID/ NXX,NWW,XX(15),WW(15),PQG(15,15,3),PQQB(15,15,2), &QGMAX(15,15,3),QQBMAX(15,15,2),YCUT(15,15),XTOT(15,15),NP COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ WMAX=SQRT(PARL(21)) WRITE(6,1000) PARL(11),LST(13),LST(9),PARL(10), &PARL(8),PARL(9),PARL(12),PARL(13) IF(NP.EQ.1) THEN WRITE(6,1010) ELSE WRITE(6,1020) ENDIF WRITE(6,1030) LST(19),NWW,NXX,WW,XX IF(WMAX.GT.WW(NWW)) WRITE(6,1040) WMAX,WW(NWW) WRITE(6,1100) DO 10 I=1,2 DO 10 J=1,5 10 P(I,J)=PCMS(I,J) LW=0 DO 500 IW=1,NWW,MAX(1,NSTEP) W=WW(IW) IF(LW.GT.0) GOTO 600 IF(W.GT.WMAX) LW=LW+1 W2=W**2 LX=0 DO 400 IX=1,NXX,MAX(1,NSTEP) X=XX(IX) IF(LX.GT.0) GOTO 500 U=(W2-P(2,5)**2)/(2.*P(2,5)*(1.-X)) Q2=2.*P(2,5)*U*X Y=Q2/(PARL(21)*X) PARI(24)=(1.+(1.-Y)**2)/2. PARI(25)=1.-Y PARI(26)=(1.-(1.-Y)**2)/2. PARL(25)=DALPS(Q2) IF(Y.GT.1.) LX=LX+1 RQG=0. RQQB=0. DO 200 IP=1,NP IF(NP.EQ.1) THEN RQG=PQG(IX,IW,IP) RQQB=PQQB(IX,IW,IP) ELSE RQG=RQG+PQG(IX,IW,IP)*PARI(23+IP)/XTOT(IX,IW) IF(IP.LT.3) RQQB=RQQB+PQQB(IX,IW,IP)*PARI(23+IP)/XTOT(IX,IW) ENDIF 200 CONTINUE C...Include alpha-strong in weight. RQG=RQG*PARL(25) RQQB=RQQB*PARL(25) RQ=1.-RQG-RQQB WRITE(6,1200) W,X,Y,Q2,PARL(25),YCUT(IX,IW),RQ,RQG,RQQB 400 CONTINUE 500 CONTINUE 600 CONTINUE RETURN 1000 FORMAT('1',/,5X,'SUMMARY OF QCD MATRIX ELEMENT INTEGRATION', & /,5X,'-----------------------------------------',//, &/,' for gluon radiation (qg-event) and boson-gluon fusion ', &'(qq-event) probability.', &//,' Required precision in integration, PARL(11) =',F8.4, &//,' Heaviest flavour produced in boson-gluon fusion, LST(13) =', &I5,//,' Alpha-strong parameters: max flavour, LST(9) =',I3, &' QCD lambda, PARL(10) =',F6.3,' GeV', &//,' Cuts on matrix elements:', &/,' PARL(8), PARL(9), PARL(12), PARL(13) =',4F8.4,/) 1010 FORMAT(' Lepton energy not allowed to vary in simulation.',/) 1020 FORMAT(' Lepton energy allowed to vary in simulation, ',/, &' y in table below calculated assuming max energy.',/) 1030 FORMAT(' Grid choice, LST(19) =',I3,5X,'# grid points in W, X =', &2I5,/,' W-values in array WW:',/,15F8.1, &/,' x-values in array XX:',/,15F8.4,/) 1040 FORMAT(' Max W outside grid, execution stopped ! Wmax, grid-max =' &,2F12.1) 1100 FORMAT(//,5X,'W',6X,'x',6X,'y',7X,'Q**2',1X,'alpha', &4X,'cut',2X,'q-event',1X,'qg-event',1X,'qq-event', &/,1X,70(1H-),/) 1200 FORMAT(F6.1,2F7.3,F11.1,F6.2,F7.3,3F9.3) END C ********************************************************************** SUBROUTINE PYSIGM(NPAR,DERIV,DIFSIG,XF,IFLAG) C...Calculates the negative of the differential cross-section. C...In the generation procedure the maximum of the differential cross- C...section is needed for weighting purposes. This maximum is found by C...minimizing the negative differential cross-section using the MINUIT C...routines which are then calling this routine. C...More precisly, only the part of the cross-section formula which is C...needed for the weighting procedure is included here. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),COMFAC COMMON /LUJETS/ N,K(2000,2),P(2000,5) DIMENSION DERIV(30),XF(30) DUMMY=NPAR+DERIV(1) IF(IFLAG.EQ.1) NCALLS=0 IF(IFLAG.EQ.2) WRITE(6,1000) DIFSIG=1.E+10 NCALLS=NCALLS+1 X=XF(1) IF(X.LT.XMIN.OR.X.GT.XMAX) RETURN IF(LST(31).EQ.0) THEN Q2=XF(2) Y=Q2/(PARL(21)*X) ELSE Y=XF(2) Q2=Y*X*PARL(21) ENDIF IF(Y.LT.YMIN.OR.Y.GT.YMAX) RETURN Q2LOW=MAX(Q2MIN,X*YMIN*PARL(21)) Q2UPP=MIN(Q2MAX,X*YMAX*PARL(21)) IF(Q2.LT.Q2LOW.OR.Q2.GT.Q2UPP) RETURN LST2=LST(2) LST(2)=-1 CALL LEPTO LST(2)=LST2 IF(LST(21).NE.0) RETURN DIFSIG=-PQ(17)*COMFAC IF(IFLAG.EQ.3) WRITE(6,1100) NCALLS,DIFSIG,X,Y,Q2,W2 RETURN 1000 FORMAT(' WARNING ! IFLAG = 2 in call to PYSIGM, which does not ' &,'calculate derivatives.') 1100 FORMAT(/,5X,'Terminating entry in PYSIGM after ',I5,' calls.',/, &5X,'Best estimate of minimum found to be ',E12.4,/, &5X,'located at x, y, Q**2, W**2 = ',4G10.3,/) END C ********************************************************************** FUNCTION DLOWER(V1) C...Lower limit on second variable, y or Q**2, depending on first C...variable x=V1. Used for integrating differential cross-section. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ IF(LST(31).EQ.0) THEN DLOWER=MAX(Q2MIN,PARL(21)*V1*YMIN,PARL(21)*YMIN-CUT(8)+0.938**2, & 2.*0.938*CUT(9)*V1) ELSE DLOWER=MAX(YMIN,Q2MIN/(PARL(21)*V1), & (CUT(7)-0.938**2)/(PARL(21)*(1.-V1))) ENDIF RETURN END C ********************************************************************** FUNCTION DUPPER(V1) C...Upper limit on second variable, y or Q**2, depending on first C...variable x=V1. Used for integrating differential cross-section. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ IF(LST(31).EQ.0) THEN DUPPER=MIN(Q2MAX,PARL(21)*V1*YMAX,PARL(21)*YMAX-CUT(7)+0.938**2, & 2.*0.938*CUT(10)*V1) ELSE DUPPER=MIN(YMAX,Q2MAX/(PARL(21)*MAX(V1,1.E-10)), & (CUT(8)-0.938**2)/(PARL(21)*MAX(1.-V1,1.E-10))) ENDIF RETURN END C ********************************************************************** FUNCTION DCROSS(V1,V2) C...Differential cross-section, dsigma/dxdQ**2 or dsigma/dxdy. C...Can be used for numerical integration etc. C...Note, non-zero result only for region defined by cuts through CUT. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LINTRL/ PLAB(2,5),PCMS(2,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX, &Q2MIN,Q2MAX,ILEP,INU,IG,IZ COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),COMFAC COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) COMMON /LINTEG/ NTOT,NPASS DCROSS=0. NTOT=NTOT+1 C...Variable V1 is x, variable V2 is either Q**2 or y. X=V1 IF(X.LT.XMIN.OR.X.GT.XMAX) RETURN IF(LST(31).EQ.0) THEN Q2=V2 Y=Q2/(PARL(21)*X) ELSE Y=V2 Q2=Y*X*PARL(21) ENDIF IF(Y.LT.YMIN.OR.Y.GT.YMAX) RETURN Q2LOW=MAX(Q2MIN,X*YMIN*PARL(21)) Q2UPP=MIN(Q2MAX,X*YMAX*PARL(21)) IF(Q2.LT.Q2LOW.OR.Q2.GT.Q2UPP) RETURN LST2=LST(2) LST(2)=-2 CALL LEPTO LST(2)=LST2 IF(LST(21).NE.0) RETURN NPASS=NPASS+1 DCROSS=PARI(31)*PQ(17)*COMFAC RETURN END C ******************************************************************** SUBROUTINE LNUCL(X,Q2,XPQ) C...Get structure function per nucleon for a proton/neutron mixture C...according to defined nucleus. COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) DIMENSION XPQ(-6:6) CALL PYSTFU(41,X,Q2,XPQ) IF(PARI(11).LE.1.E-06) RETURN C...For nuclear target, mix u- and d-valence distributions. XUV=XPQ(1)-XPQ(-1) XDV=XPQ(2)-XPQ(-2) XPQ(1)=(1.-PARI(11))*XUV+PARI(11)*XDV + XPQ(-1) XPQ(2)=(1.-PARI(11))*XDV+PARI(11)*XUV + XPQ(-2) C...Save valence u- and d-probabilities. PARI(12)=XUV PARI(13)=XDV RETURN END C ********************************************************************** SUBROUTINE PYSTFU(KF,X,Q2,XPQ) C...Gives proton and pi+ parton structure functions according C...to a few different parametrizations. C...Note that what is coded is x times the probability distribution, C...i.e. xq(x,Q2) etc. COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),XP,YP,W2P,Q2P,UP COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17) COMMON /LUDAT2/ KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) DIMENSION XPQ(-6:6),XQ(6),TX(6),TT(6),TS(6),NEHLQ(8,2), &CEHLQ(6,6,2,8,2),CDO(3,6,5,2),COW(3,5,4,2) C...The following data lines are coefficients needed in the C...Eichten, Hinchliffe, Lane, Quigg proton structure function C...parametrizations, see below. C...Powers of 1-x in different cases DATA NEHLQ/3,4,7,5,7,7,7,7,3,4,7,6,7,7,7,7/ C...Expansion coefficients for up valence quark distribution DATA (((CEHLQ(IX,IT,NX,1,1),IX=1,6),IT=1,6),NX=1,2)/ 1 7.677E-01,-2.087E-01,-3.303E-01,-2.517E-02,-1.570E-02,-1.000E-04, 2-5.326E-01,-2.661E-01, 3.201E-01, 1.192E-01, 2.434E-02, 7.620E-03, 3 2.162E-01, 1.881E-01,-8.375E-02,-6.515E-02,-1.743E-02,-5.040E-03, 4-9.211E-02,-9.952E-02, 1.373E-02, 2.506E-02, 8.770E-03, 2.550E-03, 5 3.670E-02, 4.409E-02, 9.600E-04,-7.960E-03,-3.420E-03,-1.050E-03, 6-1.549E-02,-2.026E-02,-3.060E-03, 2.220E-03, 1.240E-03, 4.100E-04, 1 2.395E-01, 2.905E-01, 9.778E-02, 2.149E-02, 3.440E-03, 5.000E-04, 2 1.751E-02,-6.090E-03,-2.687E-02,-1.916E-02,-7.970E-03,-2.750E-03, 3-5.760E-03,-5.040E-03, 1.080E-03, 2.490E-03, 1.530E-03, 7.500E-04, 4 1.740E-03, 1.960E-03, 3.000E-04,-3.400E-04,-2.900E-04,-1.800E-04, 5-5.300E-04,-6.400E-04,-1.700E-04, 4.000E-05, 6.000E-05, 4.000E-05, 6 1.700E-04, 2.200E-04, 8.000E-05, 1.000E-05,-1.000E-05,-1.000E-05/ DATA (((CEHLQ(IX,IT,NX,1,2),IX=1,6),IT=1,6),NX=1,2)/ 1 7.237E-01,-2.189E-01,-2.995E-01,-1.909E-02,-1.477E-02, 2.500E-04, 2-5.314E-01,-2.425E-01, 3.283E-01, 1.119E-01, 2.223E-02, 7.070E-03, 3 2.289E-01, 1.890E-01,-9.859E-02,-6.900E-02,-1.747E-02,-5.080E-03, 4-1.041E-01,-1.084E-01, 2.108E-02, 2.975E-02, 9.830E-03, 2.830E-03, 5 4.394E-02, 5.116E-02,-1.410E-03,-1.055E-02,-4.230E-03,-1.270E-03, 6-1.991E-02,-2.539E-02,-2.780E-03, 3.430E-03, 1.720E-03, 5.500E-04, 1 2.410E-01, 2.884E-01, 9.369E-02, 1.900E-02, 2.530E-03, 2.400E-04, 2 1.765E-02,-9.220E-03,-3.037E-02,-2.085E-02,-8.440E-03,-2.810E-03, 3-6.450E-03,-5.260E-03, 1.720E-03, 3.110E-03, 1.830E-03, 8.700E-04, 4 2.120E-03, 2.320E-03, 2.600E-04,-4.900E-04,-3.900E-04,-2.300E-04, 5-6.900E-04,-8.200E-04,-2.000E-04, 7.000E-05, 9.000E-05, 6.000E-05, 6 2.400E-04, 3.100E-04, 1.100E-04, 0.000E+00,-2.000E-05,-2.000E-05/ C...Expansion coefficients for down valence quark distribution DATA (((CEHLQ(IX,IT,NX,2,1),IX=1,6),IT=1,6),NX=1,2)/ 1 3.813E-01,-8.090E-02,-1.634E-01,-2.185E-02,-8.430E-03,-6.200E-04, 2-2.948E-01,-1.435E-01, 1.665E-01, 6.638E-02, 1.473E-02, 4.080E-03, 3 1.252E-01, 1.042E-01,-4.722E-02,-3.683E-02,-1.038E-02,-2.860E-03, 4-5.478E-02,-5.678E-02, 8.900E-03, 1.484E-02, 5.340E-03, 1.520E-03, 5 2.220E-02, 2.567E-02,-3.000E-05,-4.970E-03,-2.160E-03,-6.500E-04, 6-9.530E-03,-1.204E-02,-1.510E-03, 1.510E-03, 8.300E-04, 2.700E-04, 1 1.261E-01, 1.354E-01, 3.958E-02, 8.240E-03, 1.660E-03, 4.500E-04, 2 3.890E-03,-1.159E-02,-1.625E-02,-9.610E-03,-3.710E-03,-1.260E-03, 3-1.910E-03,-5.600E-04, 1.590E-03, 1.590E-03, 8.400E-04, 3.900E-04, 4 6.400E-04, 4.900E-04,-1.500E-04,-2.900E-04,-1.800E-04,-1.000E-04, 5-2.000E-04,-1.900E-04, 0.000E+00, 6.000E-05, 4.000E-05, 3.000E-05, 6 7.000E-05, 8.000E-05, 2.000E-05,-1.000E-05,-1.000E-05,-1.000E-05/ DATA (((CEHLQ(IX,IT,NX,2,2),IX=1,6),IT=1,6),NX=1,2)/ 1 3.578E-01,-8.622E-02,-1.480E-01,-1.840E-02,-7.820E-03,-4.500E-04, 2-2.925E-01,-1.304E-01, 1.696E-01, 6.243E-02, 1.353E-02, 3.750E-03, 3 1.318E-01, 1.041E-01,-5.486E-02,-3.872E-02,-1.038E-02,-2.850E-03, 4-6.162E-02,-6.143E-02, 1.303E-02, 1.740E-02, 5.940E-03, 1.670E-03, 5 2.643E-02, 2.957E-02,-1.490E-03,-6.450E-03,-2.630E-03,-7.700E-04, 6-1.218E-02,-1.497E-02,-1.260E-03, 2.240E-03, 1.120E-03, 3.500E-04, 1 1.263E-01, 1.334E-01, 3.732E-02, 7.070E-03, 1.260E-03, 3.400E-04, 2 3.660E-03,-1.357E-02,-1.795E-02,-1.031E-02,-3.880E-03,-1.280E-03, 3-2.100E-03,-3.600E-04, 2.050E-03, 1.920E-03, 9.800E-04, 4.400E-04, 4 7.700E-04, 5.400E-04,-2.400E-04,-3.900E-04,-2.400E-04,-1.300E-04, 5-2.600E-04,-2.300E-04, 2.000E-05, 9.000E-05, 6.000E-05, 4.000E-05, 6 9.000E-05, 1.000E-04, 2.000E-05,-2.000E-05,-2.000E-05,-1.000E-05/ C...Expansion coefficients for up and down sea quark distributions DATA (((CEHLQ(IX,IT,NX,3,1),IX=1,6),IT=1,6),NX=1,2)/ 1 6.870E-02,-6.861E-02, 2.973E-02,-5.400E-03, 3.780E-03,-9.700E-04, 2-1.802E-02, 1.400E-04, 6.490E-03,-8.540E-03, 1.220E-03,-1.750E-03, 3-4.650E-03, 1.480E-03,-5.930E-03, 6.000E-04,-1.030E-03,-8.000E-05, 4 6.440E-03, 2.570E-03, 2.830E-03, 1.150E-03, 7.100E-04, 3.300E-04, 5-3.930E-03,-2.540E-03,-1.160E-03,-7.700E-04,-3.600E-04,-1.900E-04, 6 2.340E-03, 1.930E-03, 5.300E-04, 3.700E-04, 1.600E-04, 9.000E-05, 1 1.014E+00,-1.106E+00, 3.374E-01,-7.444E-02, 8.850E-03,-8.700E-04, 2 9.233E-01,-1.285E+00, 4.475E-01,-9.786E-02, 1.419E-02,-1.120E-03, 3 4.888E-02,-1.271E-01, 8.606E-02,-2.608E-02, 4.780E-03,-6.000E-04, 4-2.691E-02, 4.887E-02,-1.771E-02, 1.620E-03, 2.500E-04,-6.000E-05, 5 7.040E-03,-1.113E-02, 1.590E-03, 7.000E-04,-2.000E-04, 0.000E+00, 6-1.710E-03, 2.290E-03, 3.800E-04,-3.500E-04, 4.000E-05, 1.000E-05/ DATA (((CEHLQ(IX,IT,NX,3,2),IX=1,6),IT=1,6),NX=1,2)/ 1 1.008E-01,-7.100E-02, 1.973E-02,-5.710E-03, 2.930E-03,-9.900E-04, 2-5.271E-02,-1.823E-02, 1.792E-02,-6.580E-03, 1.750E-03,-1.550E-03, 3 1.220E-02, 1.763E-02,-8.690E-03,-8.800E-04,-1.160E-03,-2.100E-04, 4-1.190E-03,-7.180E-03, 2.360E-03, 1.890E-03, 7.700E-04, 4.100E-04, 5-9.100E-04, 2.040E-03,-3.100E-04,-1.050E-03,-4.000E-04,-2.400E-04, 6 1.190E-03,-1.700E-04,-2.000E-04, 4.200E-04, 1.700E-04, 1.000E-04, 1 1.081E+00,-1.189E+00, 3.868E-01,-8.617E-02, 1.115E-02,-1.180E-03, 2 9.917E-01,-1.396E+00, 4.998E-01,-1.159E-01, 1.674E-02,-1.720E-03, 3 5.099E-02,-1.338E-01, 9.173E-02,-2.885E-02, 5.890E-03,-6.500E-04, 4-3.178E-02, 5.703E-02,-2.070E-02, 2.440E-03, 1.100E-04,-9.000E-05, 5 8.970E-03,-1.392E-02, 2.050E-03, 6.500E-04,-2.300E-04, 2.000E-05, 6-2.340E-03, 3.010E-03, 5.000E-04,-3.900E-04, 6.000E-05, 1.000E-05/ C...Expansion coefficients for gluon distribution DATA (((CEHLQ(IX,IT,NX,4,1),IX=1,6),IT=1,6),NX=1,2)/ 1 9.482E-01,-9.578E-01, 1.009E-01,-1.051E-01, 3.456E-02,-3.054E-02, 2-9.627E-01, 5.379E-01, 3.368E-01,-9.525E-02, 1.488E-02,-2.051E-02, 3 4.300E-01,-8.306E-02,-3.372E-01, 4.902E-02,-9.160E-03, 1.041E-02, 4-1.925E-01,-1.790E-02, 2.183E-01, 7.490E-03, 4.140E-03,-1.860E-03, 5 8.183E-02, 1.926E-02,-1.072E-01,-1.944E-02,-2.770E-03,-5.200E-04, 6-3.884E-02,-1.234E-02, 5.410E-02, 1.879E-02, 3.350E-03, 1.040E-03, 1 2.948E+01,-3.902E+01, 1.464E+01,-3.335E+00, 5.054E-01,-5.915E-02, 2 2.559E+01,-3.955E+01, 1.661E+01,-4.299E+00, 6.904E-01,-8.243E-02, 3-1.663E+00, 1.176E+00, 1.118E+00,-7.099E-01, 1.948E-01,-2.404E-02, 4-2.168E-01, 8.170E-01,-7.169E-01, 1.851E-01,-1.924E-02,-3.250E-03, 5 2.088E-01,-4.355E-01, 2.239E-01,-2.446E-02,-3.620E-03, 1.910E-03, 6-9.097E-02, 1.601E-01,-5.681E-02,-2.500E-03, 2.580E-03,-4.700E-04/ DATA (((CEHLQ(IX,IT,NX,4,2),IX=1,6),IT=1,6),NX=1,2)/ 1 2.367E+00, 4.453E-01, 3.660E-01, 9.467E-02, 1.341E-01, 1.661E-02, 2-3.170E+00,-1.795E+00, 3.313E-02,-2.874E-01,-9.827E-02,-7.119E-02, 3 1.823E+00, 1.457E+00,-2.465E-01, 3.739E-02, 6.090E-03, 1.814E-02, 4-1.033E+00,-9.827E-01, 2.136E-01, 1.169E-01, 5.001E-02, 1.684E-02, 5 5.133E-01, 5.259E-01,-1.173E-01,-1.139E-01,-4.988E-02,-2.021E-02, 6-2.881E-01,-3.145E-01, 5.667E-02, 9.161E-02, 4.568E-02, 1.951E-02, 1 3.036E+01,-4.062E+01, 1.578E+01,-3.699E+00, 6.020E-01,-7.031E-02, 2 2.700E+01,-4.167E+01, 1.770E+01,-4.804E+00, 7.862E-01,-1.060E-01, 3-1.909E+00, 1.357E+00, 1.127E+00,-7.181E-01, 2.232E-01,-2.481E-02, 4-2.488E-01, 9.781E-01,-8.127E-01, 2.094E-01,-2.997E-02,-4.710E-03, 5 2.506E-01,-5.427E-01, 2.672E-01,-3.103E-02,-1.800E-03, 2.870E-03, 6-1.128E-01, 2.087E-01,-6.972E-02,-2.480E-03, 2.630E-03,-8.400E-04/ C...Expansion coefficients for strange sea quark distribution DATA (((CEHLQ(IX,IT,NX,5,1),IX=1,6),IT=1,6),NX=1,2)/ 1 4.968E-02,-4.173E-02, 2.102E-02,-3.270E-03, 3.240E-03,-6.700E-04, 2-6.150E-03,-1.294E-02, 6.740E-03,-6.890E-03, 9.000E-04,-1.510E-03, 3-8.580E-03, 5.050E-03,-4.900E-03,-1.600E-04,-9.400E-04,-1.500E-04, 4 7.840E-03, 1.510E-03, 2.220E-03, 1.400E-03, 7.000E-04, 3.500E-04, 5-4.410E-03,-2.220E-03,-8.900E-04,-8.500E-04,-3.600E-04,-2.000E-04, 6 2.520E-03, 1.840E-03, 4.100E-04, 3.900E-04, 1.600E-04, 9.000E-05, 1 9.235E-01,-1.085E+00, 3.464E-01,-7.210E-02, 9.140E-03,-9.100E-04, 2 9.315E-01,-1.274E+00, 4.512E-01,-9.775E-02, 1.380E-02,-1.310E-03, 3 4.739E-02,-1.296E-01, 8.482E-02,-2.642E-02, 4.760E-03,-5.700E-04, 4-2.653E-02, 4.953E-02,-1.735E-02, 1.750E-03, 2.800E-04,-6.000E-05, 5 6.940E-03,-1.132E-02, 1.480E-03, 6.500E-04,-2.100E-04, 0.000E+00, 6-1.680E-03, 2.340E-03, 4.200E-04,-3.400E-04, 5.000E-05, 1.000E-05/ DATA (((CEHLQ(IX,IT,NX,5,2),IX=1,6),IT=1,6),NX=1,2)/ 1 6.478E-02,-4.537E-02, 1.643E-02,-3.490E-03, 2.710E-03,-6.700E-04, 2-2.223E-02,-2.126E-02, 1.247E-02,-6.290E-03, 1.120E-03,-1.440E-03, 3-1.340E-03, 1.362E-02,-6.130E-03,-7.900E-04,-9.000E-04,-2.000E-04, 4 5.080E-03,-3.610E-03, 1.700E-03, 1.830E-03, 6.800E-04, 4.000E-04, 5-3.580E-03, 6.000E-05,-2.600E-04,-1.050E-03,-3.800E-04,-2.300E-04, 6 2.420E-03, 9.300E-04,-1.000E-04, 4.500E-04, 1.700E-04, 1.100E-04, 1 9.868E-01,-1.171E+00, 3.940E-01,-8.459E-02, 1.124E-02,-1.250E-03, 2 1.001E+00,-1.383E+00, 5.044E-01,-1.152E-01, 1.658E-02,-1.830E-03, 3 4.928E-02,-1.368E-01, 9.021E-02,-2.935E-02, 5.800E-03,-6.600E-04, 4-3.133E-02, 5.785E-02,-2.023E-02, 2.630E-03, 1.600E-04,-8.000E-05, 5 8.840E-03,-1.416E-02, 1.900E-03, 5.800E-04,-2.500E-04, 1.000E-05, 6-2.300E-03, 3.080E-03, 5.500E-04,-3.700E-04, 7.000E-05, 1.000E-05/ C...Expansion coefficients for charm sea quark distribution DATA (((CEHLQ(IX,IT,NX,6,1),IX=1,6),IT=1,6),NX=1,2)/ 1 9.270E-03,-1.817E-02, 9.590E-03,-6.390E-03, 1.690E-03,-1.540E-03, 2 5.710E-03,-1.188E-02, 6.090E-03,-4.650E-03, 1.240E-03,-1.310E-03, 3-3.960E-03, 7.100E-03,-3.590E-03, 1.840E-03,-3.900E-04, 3.400E-04, 4 1.120E-03,-1.960E-03, 1.120E-03,-4.800E-04, 1.000E-04,-4.000E-05, 5 4.000E-05,-3.000E-05,-1.800E-04, 9.000E-05,-5.000E-05,-2.000E-05, 6-4.200E-04, 7.300E-04,-1.600E-04, 5.000E-05, 5.000E-05, 5.000E-05, 1 8.098E-01,-1.042E+00, 3.398E-01,-6.824E-02, 8.760E-03,-9.000E-04, 2 8.961E-01,-1.217E+00, 4.339E-01,-9.287E-02, 1.304E-02,-1.290E-03, 3 3.058E-02,-1.040E-01, 7.604E-02,-2.415E-02, 4.600E-03,-5.000E-04, 4-2.451E-02, 4.432E-02,-1.651E-02, 1.430E-03, 1.200E-04,-1.000E-04, 5 1.122E-02,-1.457E-02, 2.680E-03, 5.800E-04,-1.200E-04, 3.000E-05, 6-7.730E-03, 7.330E-03,-7.600E-04,-2.400E-04, 1.000E-05, 0.000E+00/ DATA (((CEHLQ(IX,IT,NX,6,2),IX=1,6),IT=1,6),NX=1,2)/ 1 9.980E-03,-1.945E-02, 1.055E-02,-6.870E-03, 1.860E-03,-1.560E-03, 2 5.700E-03,-1.203E-02, 6.250E-03,-4.860E-03, 1.310E-03,-1.370E-03, 3-4.490E-03, 7.990E-03,-4.170E-03, 2.050E-03,-4.400E-04, 3.300E-04, 4 1.470E-03,-2.480E-03, 1.460E-03,-5.700E-04, 1.200E-04,-1.000E-05, 5-9.000E-05, 1.500E-04,-3.200E-04, 1.200E-04,-6.000E-05,-4.000E-05, 6-4.200E-04, 7.600E-04,-1.400E-04, 4.000E-05, 7.000E-05, 5.000E-05, 1 8.698E-01,-1.131E+00, 3.836E-01,-8.111E-02, 1.048E-02,-1.300E-03, 2 9.626E-01,-1.321E+00, 4.854E-01,-1.091E-01, 1.583E-02,-1.700E-03, 3 3.057E-02,-1.088E-01, 8.022E-02,-2.676E-02, 5.590E-03,-5.600E-04, 4-2.845E-02, 5.164E-02,-1.918E-02, 2.210E-03,-4.000E-05,-1.500E-04, 5 1.311E-02,-1.751E-02, 3.310E-03, 5.100E-04,-1.200E-04, 5.000E-05, 6-8.590E-03, 8.380E-03,-9.200E-04,-2.600E-04, 1.000E-05,-1.000E-05/ C...Expansion coefficients for bottom sea quark distribution DATA (((CEHLQ(IX,IT,NX,7,1),IX=1,6),IT=1,6),NX=1,2)/ 1 9.010E-03,-1.401E-02, 7.150E-03,-4.130E-03, 1.260E-03,-1.040E-03, 2 6.280E-03,-9.320E-03, 4.780E-03,-2.890E-03, 9.100E-04,-8.200E-04, 3-2.930E-03, 4.090E-03,-1.890E-03, 7.600E-04,-2.300E-04, 1.400E-04, 4 3.900E-04,-1.200E-03, 4.400E-04,-2.500E-04, 2.000E-05,-2.000E-05, 5 2.600E-04, 1.400E-04,-8.000E-05, 1.000E-04, 1.000E-05, 1.000E-05, 6-2.600E-04, 3.200E-04, 1.000E-05,-1.000E-05, 1.000E-05,-1.000E-05, 1 8.029E-01,-1.075E+00, 3.792E-01,-7.843E-02, 1.007E-02,-1.090E-03, 2 7.903E-01,-1.099E+00, 4.153E-01,-9.301E-02, 1.317E-02,-1.410E-03, 3-1.704E-02,-1.130E-02, 2.882E-02,-1.341E-02, 3.040E-03,-3.600E-04, 4-7.200E-04, 7.230E-03,-5.160E-03, 1.080E-03,-5.000E-05,-4.000E-05, 5 3.050E-03,-4.610E-03, 1.660E-03,-1.300E-04,-1.000E-05, 1.000E-05, 6-4.360E-03, 5.230E-03,-1.610E-03, 2.000E-04,-2.000E-05, 0.000E+00/ DATA (((CEHLQ(IX,IT,NX,7,2),IX=1,6),IT=1,6),NX=1,2)/ 1 8.980E-03,-1.459E-02, 7.510E-03,-4.410E-03, 1.310E-03,-1.070E-03, 2 5.970E-03,-9.440E-03, 4.800E-03,-3.020E-03, 9.100E-04,-8.500E-04, 3-3.050E-03, 4.440E-03,-2.100E-03, 8.500E-04,-2.400E-04, 1.400E-04, 4 5.300E-04,-1.300E-03, 5.600E-04,-2.700E-04, 3.000E-05,-2.000E-05, 5 2.000E-04, 1.400E-04,-1.100E-04, 1.000E-04, 0.000E+00, 0.000E+00, 6-2.600E-04, 3.200E-04, 0.000E+00,-3.000E-05, 1.000E-05,-1.000E-05, 1 8.672E-01,-1.174E+00, 4.265E-01,-9.252E-02, 1.244E-02,-1.460E-03, 2 8.500E-01,-1.194E+00, 4.630E-01,-1.083E-01, 1.614E-02,-1.830E-03, 3-2.241E-02,-5.630E-03, 2.815E-02,-1.425E-02, 3.520E-03,-4.300E-04, 4-7.300E-04, 8.030E-03,-5.780E-03, 1.380E-03,-1.300E-04,-4.000E-05, 5 3.460E-03,-5.380E-03, 1.960E-03,-2.100E-04, 1.000E-05, 1.000E-05, 6-4.850E-03, 5.950E-03,-1.890E-03, 2.600E-04,-3.000E-05, 0.000E+00/ C...Expansion coefficients for top sea quark distribution DATA (((CEHLQ(IX,IT,NX,8,1),IX=1,6),IT=1,6),NX=1,2)/ 1 4.410E-03,-7.480E-03, 3.770E-03,-2.580E-03, 7.300E-04,-7.100E-04, 2 3.840E-03,-6.050E-03, 3.030E-03,-2.030E-03, 5.800E-04,-5.900E-04, 3-8.800E-04, 1.660E-03,-7.500E-04, 4.700E-04,-1.000E-04, 1.000E-04, 4-8.000E-05,-1.500E-04, 1.200E-04,-9.000E-05, 3.000E-05, 0.000E+00, 5 1.300E-04,-2.200E-04,-2.000E-05,-2.000E-05,-2.000E-05,-2.000E-05, 6-7.000E-05, 1.900E-04,-4.000E-05, 2.000E-05, 0.000E+00, 0.000E+00, 1 6.623E-01,-9.248E-01, 3.519E-01,-7.930E-02, 1.110E-02,-1.180E-03, 2 6.380E-01,-9.062E-01, 3.582E-01,-8.479E-02, 1.265E-02,-1.390E-03, 3-2.581E-02, 2.125E-02, 4.190E-03,-4.980E-03, 1.490E-03,-2.100E-04, 4 7.100E-04, 5.300E-04,-1.270E-03, 3.900E-04,-5.000E-05,-1.000E-05, 5 3.850E-03,-5.060E-03, 1.860E-03,-3.500E-04, 4.000E-05, 0.000E+00, 6-3.530E-03, 4.460E-03,-1.500E-03, 2.700E-04,-3.000E-05, 0.000E+00/ DATA (((CEHLQ(IX,IT,NX,8,2),IX=1,6),IT=1,6),NX=1,2)/ 1 4.260E-03,-7.530E-03, 3.830E-03,-2.680E-03, 7.600E-04,-7.300E-04, 2 3.640E-03,-6.050E-03, 3.030E-03,-2.090E-03, 5.900E-04,-6.000E-04, 3-9.200E-04, 1.710E-03,-8.200E-04, 5.000E-04,-1.200E-04, 1.000E-04, 4-5.000E-05,-1.600E-04, 1.300E-04,-9.000E-05, 3.000E-05, 0.000E+00, 5 1.300E-04,-2.100E-04,-1.000E-05,-2.000E-05,-2.000E-05,-1.000E-05, 6-8.000E-05, 1.800E-04,-5.000E-05, 2.000E-05, 0.000E+00, 0.000E+00, 1 7.146E-01,-1.007E+00, 3.932E-01,-9.246E-02, 1.366E-02,-1.540E-03, 2 6.856E-01,-9.828E-01, 3.977E-01,-9.795E-02, 1.540E-02,-1.790E-03, 3-3.053E-02, 2.758E-02, 2.150E-03,-4.880E-03, 1.640E-03,-2.500E-04, 4 9.200E-04, 4.200E-04,-1.340E-03, 4.600E-04,-8.000E-05,-1.000E-05, 5 4.230E-03,-5.660E-03, 2.140E-03,-4.300E-04, 6.000E-05, 0.000E+00, 6-3.890E-03, 5.000E-03,-1.740E-03, 3.300E-04,-4.000E-05, 0.000E+00/ C...The following data lines are coefficients needed in the C...Duke, Owens proton structure function parametrizations, see below. C...Expansion coefficients for (up+down) valence quark distribution DATA ((CDO(IP,IS,1,1),IS=1,6),IP=1,3)/ 1 4.190E-01, 3.460E+00, 4.400E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2 4.000E-03, 7.240E-01,-4.860E+00, 0.000E+00, 0.000E+00, 0.000E+00, 3-7.000E-03,-6.600E-02, 1.330E+00, 0.000E+00, 0.000E+00, 0.000E+00/ DATA ((CDO(IP,IS,1,2),IS=1,6),IP=1,3)/ 1 3.740E-01, 3.330E+00, 6.030E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2 1.400E-02, 7.530E-01,-6.220E+00, 0.000E+00, 0.000E+00, 0.000E+00, 3 0.000E+00,-7.600E-02, 1.560E+00, 0.000E+00, 0.000E+00, 0.000E+00/ C...Expansion coefficients for down valence quark distribution DATA ((CDO(IP,IS,2,1),IS=1,6),IP=1,3)/ 1 7.630E-01, 4.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2-2.370E-01, 6.270E-01,-4.210E-01, 0.000E+00, 0.000E+00, 0.000E+00, 3 2.600E-02,-1.900E-02, 3.300E-02, 0.000E+00, 0.000E+00, 0.000E+00/ DATA ((CDO(IP,IS,2,2),IS=1,6),IP=1,3)/ 1 7.610E-01, 3.830E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2-2.320E-01, 6.270E-01,-4.180E-01, 0.000E+00, 0.000E+00, 0.000E+00, 3 2.300E-02,-1.900E-02, 3.600E-02, 0.000E+00, 0.000E+00, 0.000E+00/ C...Expansion coefficients for (up+down+strange) sea quark distribution DATA ((CDO(IP,IS,3,1),IS=1,6),IP=1,3)/ 1 1.265E+00, 0.000E+00, 8.050E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2-1.132E+00,-3.720E-01, 1.590E+00, 6.310E+00,-1.050E+01, 1.470E+01, 3 2.930E-01,-2.900E-02,-1.530E-01,-2.730E-01,-3.170E+00, 9.800E+00/ DATA ((CDO(IP,IS,3,2),IS=1,6),IP=1,3)/ 1 1.670E+00, 0.000E+00, 9.150E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2-1.920E+00,-2.730E-01, 5.300E-01, 1.570E+01,-1.010E+02, 2.230E+02, 3 5.820E-01,-1.640E-01,-7.630E-01,-2.830E+00, 4.470E+01,-1.170E+02/ C...Expansion coefficients for charm sea quark distribution DATA ((CDO(IP,IS,4,1),IS=1,6),IP=1,3)/ 1 0.000E+00,-3.600E-02, 6.350E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2 1.350E-01,-2.220E-01, 3.260E+00,-3.030E+00, 1.740E+01,-1.790E+01, 3-7.500E-02,-5.800E-02,-9.090E-01, 1.500E+00,-1.130E+01, 1.560E+01/ DATA ((CDO(IP,IS,4,2),IS=1,6),IP=1,3)/ 1 0.000E+00,-1.200E-01, 3.510E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2 6.700E-02,-2.330E-01, 3.660E+00,-4.740E-01, 9.500E+00,-1.660E+01, 3-3.100E-02,-2.300E-02,-4.530E-01, 3.580E-01,-5.430E+00, 1.550E+01/ C...Expansion coefficients for gluon distribution DATA ((CDO(IP,IS,5,1),IS=1,6),IP=1,3)/ 1 1.560E+00, 0.000E+00, 6.000E+00, 9.000E+00, 0.000E+00, 0.000E+00, 2-1.710E+00,-9.490E-01, 1.440E+00,-7.190E+00,-1.650E+01, 1.530E+01, 3 6.380E-01, 3.250E-01,-1.050E+00, 2.550E-01, 1.090E+01,-1.010E+01/ DATA ((CDO(IP,IS,5,2),IS=1,6),IP=1,3)/ 1 8.790E-01, 0.000E+00, 4.000E+00, 9.000E+00, 0.000E+00, 0.000E+00, 2-9.710E-01,-1.160E+00, 1.230E+00,-5.640E+00,-7.540E+00,-5.960E-01, 3 4.340E-01, 4.760E-01,-2.540E-01,-8.170E-01, 5.500E+00, 1.260E-01/ C...The following data lines are coefficients needed in the C...Owens pion structure function parametrizations, see below. C...Expansion coefficients for up and down valence quark distributions DATA ((COW(IP,IS,1,1),IS=1,5),IP=1,3)/ 1 4.0000E-01, 7.0000E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2 -6.2120E-02, 6.4780E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 3 -7.1090E-03, 1.3350E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00/ DATA ((COW(IP,IS,1,2),IS=1,5),IP=1,3)/ 1 4.0000E-01, 6.2800E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2 -5.9090E-02, 6.4360E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 3 -6.5240E-03, 1.4510E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00/ C...Expansion coefficients for gluon distribution DATA ((COW(IP,IS,2,1),IS=1,5),IP=1,3)/ 1 8.8800E-01, 0.0000E+00, 3.1100E+00, 6.0000E+00, 0.0000E+00, 2 -1.8020E+00, -1.5760E+00, -1.3170E-01, 2.8010E+00, -1.7280E+01, 3 1.8120E+00, 1.2000E+00, 5.0680E-01, -1.2160E+01, 2.0490E+01/ DATA ((COW(IP,IS,2,2),IS=1,5),IP=1,3)/ 1 7.9400E-01, 0.0000E+00, 2.8900E+00, 6.0000E+00, 0.0000E+00, 2 -9.1440E-01, -1.2370E+00, 5.9660E-01, -3.6710E+00, -8.1910E+00, 3 5.9660E-01, 6.5820E-01, -2.5500E-01, -2.3040E+00, 7.7580E+00/ C...Expansion coefficients for (up+down+strange) quark sea distribution DATA ((COW(IP,IS,3,1),IS=1,5),IP=1,3)/ 1 9.0000E-01, 0.0000E+00, 5.0000E+00, 0.0000E+00, 0.0000E+00, 2 -2.4280E-01, -2.1200E-01, 8.6730E-01, 1.2660E+00, 2.3820E+00, 3 1.3860E-01, 3.6710E-03, 4.7470E-02, -2.2150E+00, 3.4820E-01/ DATA ((COW(IP,IS,3,2),IS=1,5),IP=1,3)/ 1 9.0000E-01, 0.0000E+00, 5.0000E+00, 0.0000E+00, 0.0000E+00, 2 -1.4170E-01, -1.6970E-01, -2.4740E+00, -2.5340E+00, 5.6210E-01, 3 -1.7400E-01, -9.6230E-02, 1.5750E+00, 1.3780E+00, -2.7010E-01/ C...Expansion coefficients for charm quark sea distribution DATA ((COW(IP,IS,4,1),IS=1,5),IP=1,3)/ 1 0.0000E+00, -2.2120E-02, 2.8940E+00, 0.0000E+00, 0.0000E+00, 2 7.9280E-02, -3.7850E-01, 9.4330E+00, 5.2480E+00, 8.3880E+00, 3 -6.1340E-02, -1.0880E-01, -1.0852E+01, -7.1870E+00, -1.1610E+01/ DATA ((COW(IP,IS,4,2),IS=1,5),IP=1,3)/ 1 0.0000E+00, -8.8200E-02, 1.9240E+00, 0.0000E+00, 0.0000E+00, 2 6.2290E-02, -2.8920E-01, 2.4240E-01, -4.4630E+00, -8.3670E-01, 3 -4.0990E-02, -1.0820E-01, 2.0360E+00, 5.2090E+00, -4.8400E-02/ C...Euler's beta function, requires ordinary gamma function EULBET(X,Y)=GAMMA(X)*GAMMA(Y)/GAMMA(X+Y) C...Reset structure functions, check x and hadron flavour ALAM=0. DO 100 IFL=-6,6 100 XPQ(IFL)=0. IF(X.LT.0..OR.X.GT.1.) THEN WRITE(6,1000) X RETURN ENDIF KFA=IABS(KF) IF(KFA.NE.17.AND.KFA.NE.41.AND.KFA.NE.42) THEN WRITE(6,1100) KF RETURN ENDIF IF(KFA.EQ.17) GOTO 190 C...Parametrizations of the proton structure functions. IF(LST(15).EQ.0) THEN C...Simple scaling structure functions, valence quarks and gluon only. XPQ(0)=3.*(1.-X)**5 XPQ(2)=(1.-X)**3*(1.274+.589*(1.-X)-1.675*(1.-X)**2) XPQ(1)=2.*XPQ(2) ELSEIF(LST(15).EQ.1.OR.LST(15).EQ.2) THEN C...Proton structure functions from Eichten, Hinchliffe, Lane, Quigg: C...Rev. Mod. Phys. 56 (1984) 579; as revised in Eichten, Hinchliffe, C...Lane, Quigg: Fermilab-PUB-86/75T (1986) C...Allowed variable range : 5 GeV2 < Q2 < 1E8 GeV2, 1E-4 < x < 1 C...Determine set, lamdba and x and t expansion variables NSET=LST(15) IF(NSET.EQ.1) ALAM=0.2 IF(NSET.EQ.2) ALAM=0.29 T=ALOG(Q2/ALAM**2) TMIN=ALOG(5./ALAM**2) TMAX=ALOG(1E8/ALAM**2) VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN))) NX=1 IF(X.LE.0.1) NX=2 IF(NX.EQ.1) VX=(2.*X-1.1)/0.9 IF(NX.EQ.2) VX=MAX(-1.,(2.*ALOG(X)+11.51293)/6.90776) C...Chebyshev polynomials for x and t expansion TX(1)=1. TX(2)=VX TX(3)=2.*VX**2-1. TX(4)=4.*VX**3-3.*VX TX(5)=8.*VX**4-8.*VX**2+1. TX(6)=16.*VX**5-20.*VX**3+5.*VX TT(1)=1. TT(2)=VT TT(3)=2.*VT**2-1. TT(4)=4.*VT**3-3.*VT TT(5)=8.*VT**4-8.*VT**2+1. TT(6)=16.*VT**5-20.*VT**3+5.*VT C...Calculate structure functions DO 120 IFL=1,6 XQSUM=0. DO 110 IT=1,6 DO 110 IX=1,6 110 XQSUM=XQSUM+CEHLQ(IX,IT,NX,IFL,NSET)*TX(IX)*TT(IT) 120 XQ(IFL)=XQSUM*(1.-X)**NEHLQ(IFL,NSET) C...Put into output array XPQ(0)=XQ(4) XPQ(1)=XQ(1)+XQ(3) XPQ(2)=XQ(2)+XQ(3) XPQ(3)=XQ(5) XPQ(4)=XQ(6) XPQ(-1)=XQ(3) XPQ(-2)=XQ(3) XPQ(-3)=XQ(5) XPQ(-4)=XQ(6) C...Special expansion for bottom (threshold effects) IF(LST(12).GE.5) THEN IF(NSET.EQ.1) TMIN=8.1905 IF(NSET.EQ.2) TMIN=7.4474 IF(T.LE.TMIN) GOTO 140 VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN))) TT(1)=1. TT(2)=VT TT(3)=2.*VT**2-1. TT(4)=4.*VT**3-3.*VT TT(5)=8.*VT**4-8.*VT**2+1. TT(6)=16.*VT**5-20.*VT**3+5.*VT XQSUM=0. DO 130 IT=1,6 DO 130 IX=1,6 130 XQSUM=XQSUM+CEHLQ(IX,IT,NX,7,NSET)*TX(IX)*TT(IT) XPQ(5)=XQSUM*(1.-X)**NEHLQ(7,NSET) XPQ(-5)=XPQ(5) 140 CONTINUE ENDIF C...Special expansion for top (threshold effects) IF(LST(12).GE.6) THEN IF(NSET.EQ.1) TMIN=11.5528 IF(NSET.EQ.2) TMIN=10.8097 TMIN=TMIN+2.*ALOG(PMAS(106)/30.) TMAX=TMAX+2.*ALOG(PMAS(106)/30.) IF(T.LE.TMIN) GOTO 160 VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN))) TT(1)=1. TT(2)=VT TT(3)=2.*VT**2-1. TT(4)=4.*VT**3-3.*VT TT(5)=8.*VT**4-8.*VT**2+1. TT(6)=16.*VT**5-20.*VT**3+5.*VT XQSUM=0. DO 150 IT=1,6 DO 150 IX=1,6 150 XQSUM=XQSUM+CEHLQ(IX,IT,NX,8,NSET)*TX(IX)*TT(IT) XPQ(6)=XQSUM*(1.-X)**NEHLQ(8,NSET) XPQ(-6)=XPQ(6) 160 CONTINUE ENDIF ELSEIF(LST(15).EQ.3.OR.LST(15).EQ.4) THEN C...Proton structure functions from Duke, Owens : C...Phys. Rev. D30 (1984) 49 C...Allowed variable range : 4 GeV2 < Q2 < approx 1E6 GeV2 C...Determine set, lambda and s expansion parameter NSET=LST(15)-2 IF(NSET.EQ.1) ALAM=0.2 IF(NSET.EQ.2) ALAM=0.4 SD=ALOG(ALOG(MAX(Q2,4.)/ALAM**2)/ALOG(4./ALAM**2)) C...Calculate structure functions DO 180 IFL=1,5 DO 170 IS=1,6 170 TS(IS)=CDO(1,IS,IFL,NSET)+CDO(2,IS,IFL,NSET)*SD+ & CDO(3,IS,IFL,NSET)*SD**2 IF(IFL.LE.2) THEN XQ(IFL)=X**TS(1)*(1.-X)**TS(2)*(1.+TS(3)*X)/(EULBET(TS(1), & TS(2)+1.)*(1.+TS(3)*TS(1)/(TS(1)+TS(2)+1.))) ELSE XQ(IFL)=TS(1)*X**TS(2)*(1.-X)**TS(3)*(1.+TS(4)*X+TS(5)*X**2+ & TS(6)*X**3) ENDIF 180 CONTINUE C...Put into output arrays XPQ(0)=XQ(5) XPQ(1)=3.*XQ(1)-XQ(2)+XQ(3)/6. XPQ(2)=XQ(2)+XQ(3)/6. XPQ(3)=XQ(3)/6. XPQ(4)=XQ(4) XPQ(-1)=XQ(3)/6. XPQ(-2)=XQ(3)/6. XPQ(-3)=XQ(3)/6. XPQ(-4)=XQ(4) ELSEIF(LST(15).EQ.5) THEN C...Proton structure functions from Gluck, Hoffmann, Reya : C...Z. Phys. C13 (1982) 119. C...Allowed variable range : 4 GeV2 < Q2 < 4E4 GeV2, 0.01 < x < 0.8 ALAM=0.4 SD=ALOG(ALOG(MAX(Q2,4.)/ALAM**2)/ALOG(4./ALAM**2)) C...Valence u-quark distribution. A=.421-.0412*SD C=2.-.6223*SD**.8 D=3.37+.4319*SD XU=2.*C*X**A*(1.-X**C)**D/EULBET(D+1.,A/C) C...Valence d-quark distribution. A=.364-.0368*SD C=2.-.5414*SD**.8 D=5.09+.3463*SD XD=C*X**A*(1.-X**C)**D/EULBET(D+1.,A/C) C...xub=xu=xubar=xd=xdbar symmetric sea distribution. A=.25+.088*SD**1.3 B=.8128*SD-2.003*SD**1.8+.0831*SD**2 C=3.97*SD D=7.+1.666*SD E=.2487*SD**2.5 F=27.8+59.68*SD XUB=A*(1.+B*X+C*X**2)*(1.-X)**D + E*EXP(-F*X) C...xs=xs=xsbar sea distribution. A=.0625+.1132*SD**1.3 B=12.64*SD-51.70*SD**1.8+38.02*SD**2 C=4.448*SD D=7.+1.562*SD E=.3081*SD**2.5 F=47.24+67.91*SD XS=A*(1.+B*X+C*X**2)*(1.-X)**D + E*EXP(-F*X) C...Charm sea is zero. C...Gluon distribution. A=.9243+2.51*SD**0.5 B=8.558-9.227*SD**0.3-.655*SD**1.5 C=53.57-68.78*SD**0.3+19.3*SD D=6.+1.454*SD E=11.29*SD**2 F=41.24+50.71*SD XG=A*(1.+B*X+C*X**2)*(1.-X)**D + E*EXP(-F*X) C...Put into output array XPQ(0)=XG XPQ(1)=XU+XUB XPQ(2)=XD+XUB XPQ(3)=XS XPQ(-1)=XUB XPQ(-2)=XUB XPQ(-3)=XS ELSE WRITE(6,1200) LST(15) ENDIF C...Isospin conjugation for neutron IF(KFA.EQ.42) THEN XPS=XPQ(1) XPQ(1)=XPQ(2) XPQ(2)=XPS ENDIF GOTO 400 C...Parametrizations of the pion (pi+) structure functions. 190 IF(LST(15).EQ.0) THEN C...Simple scaling structure functions, valence quarks and gluon only. XPQ(0)=2.*(1.-X)**3 XPQ(1)=SQRT(X)*(1.-X) XPQ(-2)=XPQ(1) ELSEIF(LST(15).GE.1.AND.LST(15).LE.5) THEN C...Pion structure functions from Owens : C...Phys. Rev. D30 (1984) 943 C...Allowed variable range : 4 GeV2 < Q2 < approx 2000 GeV2 C...Determine set, lambda and s expansion variable NSET=LST(15) IF(LST(15).GE.3) NSET=MIN(LST(15)-2,2) IF(NSET.EQ.1) ALAM=0.2 IF(NSET.EQ.2) ALAM=0.4 SD=ALOG(ALOG(MAX(Q2,4.)/ALAM**2)/ALOG(4./ALAM**2)) C...Calculate structure functions DO 210 IFL=1,4 DO 200 IS=1,5 200 TS(IS)=COW(1,IS,IFL,NSET)+COW(2,IS,IFL,NSET)*SD+ & COW(3,IS,IFL,NSET)*SD**2 IF(IFL.EQ.1) THEN XQ(IFL)=X**TS(1)*(1.-X)**TS(2)/EULBET(TS(1),TS(2)+1.) ELSE XQ(IFL)=TS(1)*X**TS(2)*(1.-X)**TS(3)*(1.+TS(4)*X+TS(5)*X**2) ENDIF 210 CONTINUE C...Put into output arrays XPQ(0)=XQ(2) XPQ(1)=XQ(1)+XQ(3)/6. XPQ(2)=XQ(3)/6. XPQ(3)=XQ(3)/6. XPQ(4)=XQ(4) XPQ(-1)=XQ(3)/6. XPQ(-2)=XQ(1)+XQ(3)/6. XPQ(-3)=XQ(3)/6. XPQ(-4)=XQ(4) ELSE WRITE(6,1200) LST(15) ENDIF C...Charge conjugation for antiparticle, positivity, max flavour 400 IF(KF.LT.0) THEN DO 410 IFL=1,4 XPS=XPQ(IFL) XPQ(IFL)=XPQ(-IFL) 410 XPQ(-IFL)=XPS ENDIF DO 420 IFL=-6,6 XPQ(IFL)=MAX(0.,XPQ(IFL)) 420 IF(IABS(IFL).GT.LST(12)) XPQ(IFL)=0. PARL(26)=ALAM C...Formats for error printouts 1000 FORMAT(1X,'ERROR: X VALUE OUTSIDE PHYSICAL RANGE; X =',1P,E12.3) 1100 FORMAT(1X,'ERROR: ILLEGAL PARTICLE CODE FOR STRUCTURE FUNCTION;', &1X,'KF =',I5) 1200 FORMAT(1X,'ERROR: ILLEGAL VALUE OF PARAMETER LST(15) IN PYSTFU;', &1X,'LST(15) =',I5) RETURN END C....................................................................... C C PURPOSE - INTEGRATE A FUNCTION F(X) C METHOD - ADAPTIVE GAUSSIAN C USAGE - CALL GADAP(A0,B0,F,EPS,SUM) C PARAMETERS A0 - LOWER LIMIT (INPUT,REAL) C B0 - UPPER LIMIT (INPUT,REAL) C F - FUNCTION F(X) TO BE INTEGRATED. MUST BE C SUPPLIED BY THE USER. (INPUT,REAL FUNCTION) C EPS - DESIRED RELATIVE ACCURACY. IF SUM IS SMALL EPS C WILL BE ABSOLUTE ACCURACY INSTEAD. (INPUT,REAL) C SUM - CALCULATED VALUE FOR THE INTEGRAL (OUTPUT,REAL) C PRECISION - SINGLE C REQ'D PROG'S - F C AUTHOR - T. JOHANSSON, LUND UNIV. COMPUTER CENTER, 1973 C REFERENCE(S) - THE AUSTRALIAN COMPUTER JOURNAL,3 P.126 AUG. -71 C C....................................................................... SUBROUTINE GADAP(A0,B0,F,EPS,SUM) COMMON/GADAP1/ NUM,IFU EXTERNAL F DIMENSION A(300),B(300),F1(300),F2(300),F3(300),S(300),N(300) 1 FORMAT(16H GADAP:I TOO BIG) DSUM(F1F,F2F,F3F,AA,BB)=5./18.*(BB-AA)*(F1F+1.6*F2F+F3F) IF(EPS.LT.1.0E-8) EPS=1.0E-8 RED=1.3 L=1 I=1 SUM=0. C=SQRT(15.)/5. A(1)=A0 B(1)=B0 F1(1)=F(0.5*(1+C)*A0+0.5*(1-C)*B0) F2(1)=F(0.5*(A0+B0)) F3(1)=F(0.5*(1-C)*A0+0.5*(1+C)*B0) IFU=3 S(1)= DSUM(F1(1),F2(1),F3(1),A0,B0) 100 CONTINUE L=L+1 N(L)=3 EPS=EPS*RED A(I+1)=A(I)+C*(B(I)-A(I)) B(I+1)=B(I) A(I+2)=A(I)+B(I)-A(I+1) B(I+2)=A(I+1) A(I+3)=A(I) B(I+3)=A(I+2) W1=A(I)+(B(I)-A(I))/5. U2=2.*W1-(A(I)+A(I+2))/2. F1(I+1)=F(A(I)+B(I)-W1) F2(I+1)=F3(I) F3(I+1)=F(B(I)-A(I+2)+W1) F1(I+2)=F(U2) F2(I+2)=F2(I) F3(I+2)=F(B(I+2)+A(I+2)-U2) F1(I+3)=F(A(I)+A(I+2)-W1) F2(I+3)=F1(I) F3(I+3)=F(W1) IFU=IFU+6 IF(IFU.GT.5000) GOTO 130 S(I+1)= DSUM(F1(I+1),F2(I+1),F3(I+1),A(I+1),B(I+1)) S(I+2)= DSUM(F1(I+2),F2(I+2),F3(I+2),A(I+2),B(I+2)) S(I+3)= DSUM(F1(I+3),F2(I+3),F3(I+3),A(I+3),B(I+3)) SS=S(I+1)+S(I+2)+S(I+3) I=I+3 IF(I.GT.300)GOTO 120 SOLD=S(I-3) IF(ABS(SOLD-SS).GT.EPS*(1.+ABS(SS))/2.) GOTO 100 SUM=SUM+SS I=I-4 N(L)=0 L=L-1 110 CONTINUE IF(L.EQ.1) GOTO 130 N(L)=N(L)-1 EPS=EPS/RED IF(N(L).NE.0) GOTO 100 I=I-1 L=L-1 GOTO 110 120 WRITE(6,1) 130 RETURN END C....................................................................... C C PURPOSE - INTEGRATE A FUNCTION F(X,Y) OF TWO VARIABLES C METHOD - ADAPTIVE GAUSSIAN IN BOTH DIRECTIONS C USAGE - CALL GADAP2(A0,B0,FL,FU,F,EPS,SUM) C PARAMETERS A0 - LOWER X-LIMIT (INPUT,REAL) C B0 - UPPER X-LIMIT (INPUT,REAL) C FL - USER SUPPLIED FUNCTION FL(X) GIVING THE LOWER C Y-LIMIT FOR A GIVEN X-VALUE C (INPUT,REAL FUNCTION) C FU - USER SUPPLIED FUNCTION FU(X) GIVING THE UPPER C Y-LIMIT FOR A GIVEN X-VALUE C (INPUT,REAL FUNCTION) C F - USER SUPPLIED FUNCTION F(X,Y) TO BE INTEGRATED C (INPUT,REAL FUNCTION) C EPS - DESIRED ACCURACY (INPUT,REAL) C SUM - CALCULATED VALUE FOR THE INTEGRAL (OUTPUT,REAL) C PRECISION - SINGLE C REQ'D PROG'S - FL,FU,F,FGADAP C AUTHOR - THOMAS JOHANSSON, LDC,1973 C C....................................................................... SUBROUTINE GADAP2(A0,B0,FL,FU,F,EPS,SUM) COMMON/GADAP1/ NUM,IFU EXTERNAL F,FL,FU DIMENSION A(300),B(300),F1(300),F2(300),F3(300),S(300),N(300) 1 FORMAT(16H GADAP:I TOO BIG) DSUM(F1F,F2F,F3F,AA,BB)=5./18.*(BB-AA)*(F1F+1.6*F2F+F3F) IF(EPS.LT.1.0E-8) EPS=1.0E-8 RED=1.4 L=1 I=1 SUM=0. C=SQRT(15.)/5. A(1)=A0 B(1)=B0 X=0.5*(1+C)*A0+0.5*(1-C)*B0 AY=FL(X) BY=FU(X) F1(1)=FGADAP(X,AY,BY,F,EPS) X=0.5*(A0+B0) AY=FL(X) BY=FU(X) F2(1)=FGADAP(X,AY,BY,F,EPS) X=0.5*(1-C)*A0+0.5*(1+C)*B0 AY=FL(X) BY=FU(X) F3(1)=FGADAP(X,AY,BY,F,EPS) IFU=3 S(1)= DSUM(F1(1),F2(1),F3(1),A0,B0) 100 CONTINUE L=L+1 N(L)=3 EPS=EPS*RED A(I+1)=A(I)+C*(B(I)-A(I)) B(I+1)=B(I) A(I+2)=A(I)+B(I)-A(I+1) B(I+2)=A(I+1) A(I+3)=A(I) B(I+3)=A(I+2) W1=A(I)+(B(I)-A(I))/5. U2=2.*W1-(A(I)+A(I+2))/2. X=A(I)+B(I)-W1 AY=FL(X) BY=FU(X) F1(I+1)=FGADAP(X,AY,BY,F,EPS) F2(I+1)=F3(I) X=B(I)-A(I+2)+W1 AY=FL(X) BY=FU(X) F3(I+1)=FGADAP(X,AY,BY,F,EPS) X=U2 AY=FL(X) BY=FU(X) F1(I+2)=FGADAP(X,AY,BY,F,EPS) F2(I+2)=F2(I) X=B(I+2)+A(I+2)-U2 AY=FL(X) BY=FU(X) F3(I+2)=FGADAP(X,AY,BY,F,EPS) X=A(I)+A(I+2)-W1 AY=FL(X) BY=FU(X) F1(I+3)=FGADAP(X,AY,BY,F,EPS) F2(I+3)=F1(I) X=W1 AY=FL(X) BY=FU(X) F3(I+3)=FGADAP(X,AY,BY,F,EPS) IFU=IFU+6 IF(IFU.GT.5000) GOTO 130 S(I+1)= DSUM(F1(I+1),F2(I+1),F3(I+1),A(I+1),B(I+1)) S(I+2)= DSUM(F1(I+2),F2(I+2),F3(I+2),A(I+2),B(I+2)) S(I+3)= DSUM(F1(I+3),F2(I+3),F3(I+3),A(I+3),B(I+3)) SS=S(I+1)+S(I+2)+S(I+3) I=I+3 IF(I.GT.300)GOTO 120 SOLD=S(I-3) IF(ABS(SOLD-SS).GT.EPS*(1.+ABS(SS))/2.) GOTO 100 SUM=SUM+SS I=I-4 N(L)=0 L=L-1 110 CONTINUE IF(L.EQ.1) GOTO 130 N(L)=N(L)-1 EPS=EPS/RED IF(N(L).NE.0) GOTO 100 I=I-1 L=L-1 GOTO 110 120 WRITE(6,1) 130 RETURN END FUNCTION FGADAP(X,A0,B0,F,EPS) COMMON/GADAP1/ NUM,IFU EXTERNAL F DIMENSION A(300),B(300),F1(300),F2(300),F3(300),S(300),N(300) 1 FORMAT(16H GADAP:I TOO BIG) DSUM(F1F,F2F,F3F,AA,BB)=5./18.*(BB-AA)*(F1F+1.6*F2F+F3F) IF(EPS.LT.1.0E-8) EPS=1.0E-8 RED=1.4 L=1 I=1 SUM=0. C=SQRT(15.)/5. A(1)=A0 B(1)=B0 F1(1)=F(X,0.5*(1+C)*A0+0.5*(1-C)*B0) F2(1)=F(X,0.5*(A0+B0)) F3(1)=F(X,0.5*(1-C)*A0+0.5*(1+C)*B0) IFU=3 S(1)= DSUM(F1(1),F2(1),F3(1),A0,B0) 100 CONTINUE L=L+1 N(L)=3 EPS=EPS*RED A(I+1)=A(I)+C*(B(I)-A(I)) B(I+1)=B(I) A(I+2)=A(I)+B(I)-A(I+1) B(I+2)=A(I+1) A(I+3)=A(I) B(I+3)=A(I+2) W1=A(I)+(B(I)-A(I))/5. U2=2.*W1-(A(I)+A(I+2))/2. F1(I+1)=F(X,A(I)+B(I)-W1) F2(I+1)=F3(I) F3(I+1)=F(X,B(I)-A(I+2)+W1) F1(I+2)=F(X,U2) F2(I+2)=F2(I) F3(I+2)=F(X,B(I+2)+A(I+2)-U2) F1(I+3)=F(X,A(I)+A(I+2)-W1) F2(I+3)=F1(I) F3(I+3)=F(X,W1) IFU=IFU+6 IF(IFU.GT.5000) GOTO 130 S(I+1)= DSUM(F1(I+1),F2(I+1),F3(I+1),A(I+1),B(I+1)) S(I+2)= DSUM(F1(I+2),F2(I+2),F3(I+2),A(I+2),B(I+2)) S(I+3)= DSUM(F1(I+3),F2(I+3),F3(I+3),A(I+3),B(I+3)) SS=S(I+1)+S(I+2)+S(I+3) I=I+3 IF(I.GT.300)GOTO 120 SOLD=S(I-3) IF(ABS(SOLD-SS).GT.EPS*(1.+ABS(SS))/2.) GOTO 100 SUM=SUM+SS I=I-4 N(L)=0 L=L-1 110 CONTINUE IF(L.EQ.1) GOTO 130 N(L)=N(L)-1 EPS=EPS/RED IF(N(L).NE.0) GOTO 100 I=I-1 L=L-1 GOTO 110 120 WRITE(6,1) 130 FGADAP=SUM EPS=EPS/RED RETURN END C*********************************************************************** SUBROUTINE PYTIME(TIME) C...Interface routine to transfer a call to some machine-dependent C...routine to get the elapsed time during the minimization step. C...This time information is, however, by no means essential. TIME=0. C...As an example, at CERN the TIMEX routine can be used: CALL TIMEX(TIME) RETURN END C*********************************************************************** C...The following routines are slightly modified minimization routines C...from the MINUIT program package. SUBROUTINE PYCMND C...This is the MINUIT routine COMAND. CC GETS IFORMATION FROM /PYMINU/ AND TAKES APPROPRIATE ACTION, CC EITHER DIRECTLY BY SKIPPING TO THE CORRESPONDING CODE IN CC PYCMND, OR BY SETTING UP A CALL TO A SUBROUTINE CC COMMON /PYMINU/ XKIN(4),UKIN(4),WKIN(4),AIN(4),BIN(4), &MAXFIN,RELUP,RELERR,RELER2,FCNMAX COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC FVAL3 = 2.0*AMIN+1.0 C . . . . . . . . . . ERROR DEF WORD7(1)=RELUP*ABS(AMIN) UP = WORD7(1) IF (UP .LE. 0.) UP = 1.0 IF (ISW(2) .GE. 1) CALL PYMPRI(1,AMIN) WORD7(1)=MAXFIN WORD7(2)=RELERR*UP NFCNMX = WORD7(1) + 0.5 IF (NFCNMX .LE. 0) NFCNMX = 1000 EPSI = WORD7(2) IF (EPSI .LE. 0.) EPSI = 0.1 * UP NEWMIN = 0 ITAUR = 0 ISW(1) = 0 CALL PYSIMP IF(ABS(DIRIN(1)).LE.ABS(EPSMAC*X(1)).AND. & ABS(DIRIN(2)).LE.ABS(EPSMAC*X(2))) THEN WRITE(6,2100) GOTO 500 ENDIF WORD7(1)=MAXFIN RELERR=RELER2*RELERR WORD7(2)=RELERR*UP NFCNMX = WORD7(1) + 0.5 IF (NFCNMX .LE. 0) NFCNMX = 1000 EPSI = WORD7(2) IF (EPSI .LE. 0.) EPSI = 0.1 * UP CALL PYSIMP 500 FCNMAX=ABS(AMIN) IF(ISW(1).GE.1) THEN WRITE(6,2200) FCNMAX=FCNMAX*1.25 ENDIF FMAX=ABS(AMIN) C . . . . . . . . . . END, EXIT WORD7(1)=0. 1100 IT = WORD7(1) + 0.5 IF (FVAL3 .EQ. AMIN .OR. IT .GT. 0) RETURN IFLAG = 3 CALL PYSIGM(NPAR,GIN,F,U,IFLAG) NFCN = NFCN + 1 IF(ABS(F).GT.FMAX) WRITE(6,2300) F RETURN 2100 FORMAT(' Warning! Stepsizes are less than machine accuracy ', &'times variable values. No further minimization attempted.') 2200 FORMAT(' Warning! Simplex minimization has not converged ', &'properly. Returned maximum increased by a factor 1.25.') 2300 FORMAT(' Warning from PYCMND: function at minimum, ',E12.4, &' is smaller than stored minimum.') END SUBROUTINE PYINTO(PINT) C...This is the MINUIT routine INTOEX. CC TRANSFORMS FROM INTERNAL COORDINATES (PINT) TO EXTERNAL CC PARAMETERS (U). THE MINIMIZING ROUTINES WHICH WORK IN CC INTERNAL COORDINATES CALL THIS ROUTINE BEFORE CALLING FCN. COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC DIMENSION PINT(2) DO 100 I= 1, NU J = LCORSP(I) IF ( J ) 100,100,50 50 CONTINUE IF (LCODE(I) .EQ. 1) GO TO 80 AL = ALIM(I) U(I) = AL + 0.5 *(SIN(PINT(J)) +1.0) * (BLIM(I) -AL) GO TO 100 80 U(I) = PINT(J) 100 CONTINUE RETURN END SUBROUTINE PYMIDA C...This is the MINUIT routine MIDATA. CC GETS PARAMETERS FROM /PYMINU/ AND /PYMINC/ CC AND SETS UP THE STARTING PARAMETER LISTS. CC CONTROL THEN PASSES TO PYCMND FOR READING THE COMMAND "CARDS". CC COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC COMMON /PYMINU/ XKIN(4),UKIN(4),WKIN(4),AIN(4),BIN(4), &MAXFIN,RELUP,RELERR,RELER2,FCNMAX COMMON /PYMINC/ NAMKIN(4),NAM(30) CHARACTER*10 NAMKIN,NAM,NAMK,BLANK CHARACTER XTITLE*60 DATA BLANK/' '/ DATA XTITLE/' FIND MINIMUM OF -(DIFFERENTIAL CROSS SECTION)'/ DATA MNINIT /0/ C . INITIALIZE NEW DATA BLOCK . . IF (MNINIT .EQ. 0) NBLOCK=0 MNINIT = 1 NBLOCK = NBLOCK + 1 VERSN = 11.79 WRITE (ISYSWR,1004) MAXINT,MAXEXT,VERSN,NBLOCK WRITE (ISYSWR,1005) DO 50 I= 1, 7 50 ISW(I) = 0 SIGMA = 0. CALL PYTIME(TIME) WRITE (ISYSWR,1110) XTITLE,TIME,EPSMAC WRITE (ISYSWR,1005) NPFIX = 0 NINT = 0 NU = 0 NPAR = 0 IFATAL = 0 WRITE (ISYSWR,1005) DO 100 I= 1, MAXEXT U(I) = 0.0 NAM(I) = BLANK ERP(I) = 0.0 ERN(I) = 0.0 LCODE(I) = 0 100 LCORSP (I) = 0 UP = 1.0 ISW(5) = 1 IUNIT = ISYSRD C . . . READ PARAMETER CARDS . . ENTRY PYMID2 DO 200 I= 1, 200 XK=XKIN(I) NAMK=NAMKIN(I) UK=UKIN(I) WK=WKIN(I) A=AIN(I) B=BIN(I) K = XK + 0.1 NU = MAX0(NU,K) IF (K .LE. 0) GO TO 250 IF (K .LE. MAXEXT) GO TO 115 IFATAL = IFATAL + 1 WRITE (ISYSWR,1009) K,MAXEXT WRITE (ISYSWR,1002) K,NAMK,UK,WK,A,B GO TO 200 115 CONTINUE IF(NAM(K).EQ.BLANK) GO TO 120 C PREVIOUSLY DEFINED PARAMETER IS BEING REDEFINED WRITE(ISYSWR,1007) IF(WERR(K).GT..0) NINT=NINT-1 120 CONTINUE NAM(K) = NAMK U(K) = UK WERR(K) = WK IF (WK .GT. 0.0) GO TO 122 C . . . FIXED PARAMETER . . . . WRITE (ISYSWR, 1002) K,NAMK,UK LCODE(K) = 0 GO TO 160 C . . . VARIABLE PARAMETER . . . 122 WRITE (ISYSWR, 1002) K,NAMK,UK,WK,A,B NINT = NINT + 1 ISW(2) = 0 IF (A) 140,130,140 130 IF (B) 140,135,140 135 LCODE(K) = 1 GO TO 160 140 IF (B-A) 145,142,150 142 IFATAL = IFATAL + 1 WRITE (ISYSWR,1010) GO TO 150 145 SAV = B B = A A = SAV WRITE (ISYSWR,1003) 150 ALIM(K) = A BLIM(K) = B LCODE(K) = 4 IF ((B-U(K))*(U(K)-A)) 153,155,160 153 IFATAL = IFATAL + 1 WRITE (ISYSWR,1011) GO TO 160 155 WRITE (ISYSWR,1006) 160 CONTINUE 200 CONTINUE IFATAL = IFATAL + 1 WRITE (ISYSWR,1012) C . . . END PARAMETER CARDS C . . . STOP IF FATAL ERROR 250 WRITE (ISYSWR,1005) IF (NINT .LE. MAXINT) GO TO 253 WRITE (ISYSWR,1008) NINT,MAXINT IFATAL = IFATAL + 1 253 IF (IFATAL .LE. 0) GO TO 280 WRITE (ISYSWR,1013) IFATAL STOP C CALCULATE STEP SIZES DIRIN 280 NPAR = 0 DO 300 K= 1, NU IF (LCODE(K) .LE. 0) GO TO 300 NPAR = NPAR + 1 LCORSP(K) = NPAR SAV = U(K) X(NPAR) = PYPINT(SAV,K) XT(NPAR) = X(NPAR) SAV2 = SAV + WERR(K) VPLU = PYPINT(SAV2,K) - X(NPAR) SAV2 = SAV - WERR(K) VMINU = PYPINT(SAV2,K) - X(NPAR) DIRIN(NPAR) = 0.5 * (ABS(VPLU) +ABS(VMINU)) G2(NPAR) = 2.0 / DIRIN(NPAR)**2 GSTEP(NPAR) = DIRIN(NPAR) IF (LCODE(K) .GT. 1) GSTEP(NPAR) = -GSTEP(NPAR) 300 CONTINUE SIGMA = 1.0E10 IUNIT = ISYSRD RETURN C... THE FORMAT BELOW IS MACHINE-DEPENDENT. (A10) , (A6,4X) , ETC. 1002 FORMAT (I10,2X,A10,2X,2G12.6,2X,2G12.6) 1003 FORMAT(' WARNING - ABOVE LIMITS HAVE BEEN REVERSED.') 1004 FORMAT (1H1/42X,21(1H*)/42X,21H* D506 MINUIT */42X, 112H* DIMENSIONS, I3, 1H/, I3, 2H */ 42X, 1'* MODIFICATION OF *',/,42X, 111H* VERSION ,F6.2,4H */42X,16H* DATA BLOCK NO. ,I3,2H *) 1005 FORMAT (4X,96(1H*)) 1006 FORMAT(' WARNING - ABOVE PARAMETER IS AT LIMIT ') 1007 FORMAT(' WARNING ******* - PARAMETER REQUESTED ON FOLLOWING', 1' CARD HAS ALREADY APPEARED. PREVIOUS VALUES IGNORED.') 1008 FORMAT('0 TOO MANY VARIABLE PARAMETERS. YOU REQUEST',I5/, +' THIS VERSION OF MINUIT IS ONLY DIMENSIONED FOR',I4//) 1009 FORMAT('0FATAL ERROR. PARAMETER NUMBER',I11,' GREATER THAN ', +'ALLOWED MAXIMUM',I4) 1010 FORMAT(' FATAL ERROR. UPPER AND LOWER LIMITS ARE EQUAL.') 1011 FORMAT(' FATAL ERROR. PARAMETER OUTSIDE LIMITS',/) 1012 FORMAT('0FATAL ERROR. MORE THAN 200 PARAMETER CARDS',/) 1013 FORMAT(/I5,' FATAL ERRORS ON PARAMETER CARDS. ABORT.',//) 1110 FORMAT(5X,A60,5X,'TIME',F8.3,' SECONDS',/,70X,'MACH. PREC.=', &E10.2) END SUBROUTINE PYMINN C...This is the MINUIT routine MINNEW. CC THIS IS THE MAIN PROGRAM, DISGUISED AS A SUBROUTINE FOR CC REASONS OF COMPATIBILITY BETWEEN SYSTEMS. IT INITIALIZES CC SOME CONSTANTS IN COMMON (INCLUDING THE LOGICAL I/O UNIT NOS.) CC THEN VERIFIES THAT FCN GIVES THE SAME VALUE WHEN CALLED CC TWICE WITH THE SAME ARGUMENTS, AND PASSES CONTROL TO PYCMND. CC COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC C UNIT NUMBERS FOR CARD READER, PRINTER, PUNCH C ISYSRD = 5 ISYSWR = 6 ISYSPU = 7 MAXINT=15 MAXEXT=30 C DETERMINE MACHINE ACCURACY EPSMAC EPSMAC = 0.5 DO 33 I= 1, 100 EPSMAC = EPSMAC * 0.5 IF ((1.0+EPSMAC) .EQ. 1.0) GO TO 35 33 CONTINUE EPSMAC = 1.0E-6 35 EPSMAC = 2.0 * EPSMAC C . . . . . . . . . 110 CONTINUE NFCN = 1 CALL PYMIDA CALL PYINTO(X) WRITE (ISYSWR,120) 120 FORMAT (/,'0FIRST ENTRY TO FCN ') CALL PYSIGM(NPAR,GIN,AMIN,U,1) CALL PYSIGM(NPAR,GIN,AMIN,U,4) CALL PYMPRI(1,AMIN) CALL PYSIGM(NPAR,GIN,F ,U,4) IF (F .NE. AMIN) GO TO 160 NFCN = 3 CALL PYCMND RETURN 160 CONTINUE WRITE (ISYSWR,880) AMIN, F STOP 880 FORMAT('0FOR THE ABOVE VALUES OF THE PARAMETERS, FCN IS TIME-', +'DEPENDENT',/,'0F = ',E22.14,' FOR FIRST CALL',/,' F =',E22.14, +' FOR SECOND') END SUBROUTINE PYMPRI (IKODE,FVAL) C...This is the MINUIT routine MPRINT. CC PRINTS THE VALUES OF THE PARAMETERS AT THE TIME OF THE CALL. CC ALSO PRINTS OTHER RELEVANT INFORMATION SUCH AS FUNCTION VALUE, CC ESTIMATED DISTANCE TO MINIMUM, PARAMETER ERRORS, STEP SIZES. CC ACCORDING TO THE VALUE OF IKODE,THE PRINTOUT IS LONG FORMAT, CC SHORT FORMAT, OR MINOS FORMAT (0,1,2) CC COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC COMMON /PYMINC/ NAMKIN(4),NAM(30) CHARACTER*10 NAMKIN,NAM C . GET TIME AND PRINT HEADINGS . CALL PYTIME(TI) WRITE (ISYSWR,1000) E = SIGMA KOUNT = 0 C . . . LOOP OVER PARAMETERS . . DO 200 I= 1, NU IF(NAM(I).EQ.' ') GOTO 200 20 L = LCORSP(I) IF (L .EQ. 0) GO TO 55 C VARIABLE PARAMETER. CALCULATE EXTERNAL ERROR IF V EXISTS IF (ISW(2) .LT. 1) GO TO 27 DX = SQRT(ABS(V(L,L)*UP)) IF (LCODE(I) .LE. 1) GO TO 26 AL = ALIM(I) BA = BLIM(I) - AL DU1 = AL + 0.5 *(SIN(X(L)+DX) +1.0) * BA - U(I) DU2 = AL + 0.5 *(SIN(X(L)-DX) +1.0) * BA - U(I) IF (DX .GT. 1.0) DU1 = BA DX = 0.5 * (ABS(DU1) + ABS(DU2)) 26 WERR(I) = DX 27 X1 = X(L) X2 = DIRIN(L) IF (IKODE .LT. 2) GO TO 29 X1 = ERP(I) X2 = ERN(I) 29 IF (KOUNT) 30,30,40 30 KOUNT = 1 WRITE (ISYSWR,1001) FVAL,NFCN,TI,E, L,I,NAM(I),U(I),WERR(I),X1,X2 GO TO 45 40 WRITE (ISYSWR,1002) L,I,NAM(I),U(I),WERR(I),X1,X2 45 IF (LCODE(I) .LE. 1) GO TO 200 IF (ABS(COS(X(L))) .LT. 0.001) WRITE (ISYSWR,1004) GO TO 200 C FIXED PARAMETER. PRINT ONLY IF IKODE .GT.0 55 IF (IKODE .EQ. 0) GO TO 200 IF (KOUNT) 60,60,70 60 KOUNT = 1 WRITE (ISYSWR,1001) FVAL,NFCN,TI,E, L,I,NAM(I),U(I) GO TO 200 70 WRITE (ISYSWR,1003) I,NAM(I),U(I) 200 CONTINUE IF (IKODE.GE.1 .AND.ISW(2).GE.1) WRITE (ISYSWR,1005) UP RETURN 1000 FORMAT(/ 4X,'FCN VALUE',5X,'CALLS',4X,'TIME',4X,' EDM ',4X , +'INT.EXT. PARAMETER VALUE ERROR INTERN.VALUE ', +'INT.STEP SIZE') 1001 FORMAT(E15.7,I7,F9.2,E11.2,I6,I4,1X,A10,4E14.5) 1002 FORMAT(1H ,41X,I6,I4,1X,A10,4E14.5) 1003 FORMAT(1H ,47X ,I4,1X,A10,4E14.5) 1004 FORMAT(1H ,52X ,'WARNING - - ABOVE PARAMETER IS AT LIMIT.') 1005 FORMAT(/45X,'ERRORS CORRESPOND TO FUNCTION CHANGE OF ',E12.4) END FUNCTION PYPINT(PEXTI,I) C...This is the MINUIT routine PINTF. CC CALCULATES THE INTERNAL PARAMETER VALUE PYPINT CORRESPONDING CC TO THE EXTERNAL VALUE PEXTI FOR PARAMETER I. CC COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC DATA BIG, SMALL / 1.570796326795 , -1.570796326795 / IGO = LCODE(I) GO TO (100,200,300,400),IGO C-- IGO = 1 MEANS NO LIMITS 100 PYPINT = PEXTI GO TO 800 200 CONTINUE 300 CONTINUE C-- IGO = 4 MEANS THERE ARE TWO LIMITS 400 ALIMI = ALIM(I) BLIMI = BLIM(I) IF (PEXTI-ALIMI) 440,500,460 440 A = SMALL 450 PYPINT = A PEXTI = ALIMI + 0.5* (BLIMI-ALIMI) *(SIN(A) +1.0) LIMSET=1 WRITE (ISYSWR,241) I GO TO 800 460 IF (BLIMI-PEXTI) 470,520,480 470 A = BIG GO TO 450 480 YY=2.0*(PEXTI-ALIMI)/(BLIMI-ALIMI) - 1.0 PYPINT = ATAN(YY/SQRT(1.0- YY**2) ) GO TO 800 500 PYPINT = SMALL GO TO 800 520 PYPINT = BIG 800 RETURN 241 FORMAT(' WARNING - VARIABLE',I3,' HAS BEEN BROUGHT BACK IN', +'SIDE LIMITS BY PYPINT.') END SUBROUTINE PYRAZZ(YNEW,PNEW) C...This is the MINUIT routine RAZZIA. CC CALLED ONLY BY SIMPLEX (AND IMPROV) TO ADD A NEW POINT CC AND REMOVE AN OLD ONE FROM THE CURRENT SIMPLEX, AND GET THE CC ESTIMATED DISTANCE TO MINIMUM. CC COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC DIMENSION PNEW(15) DO 10 I=1,NPAR 10 P(I,JH)=PNEW(I) Y(JH)=YNEW IF(YNEW.GE.AMIN) GO TO 18 DO 15 I=1,NPAR 15 X(I)=PNEW(I) CALL PYINTO(X) AMIN=YNEW JL=JH 18 CONTINUE JH=1 NPARP1=NPAR+1 20 DO 25 J=2,NPARP1 IF (Y(J) .GT. Y(JH)) JH = J 25 CONTINUE SIGMA = Y(JH) - Y(JL) IF (SIGMA .LE. 0.) GO TO 45 US = 1.0/SIGMA DO 35 I= 1, NPAR PBIG = P(I,1) PLIT = PBIG DO 30 J= 2, NPARP1 IF (P(I,J) .GT. PBIG) PBIG = P(I,J) IF (P(I,J) .LT. PLIT) PLIT = P(I,J) 30 CONTINUE DIRIN(I) = PBIG - PLIT IF (ITAUR .LT. 1 ) V(I,I) = 0.5*(V(I,I) +US*DIRIN(I)**2) 35 CONTINUE 40 RETURN 45 WRITE (ISYSWR, 1000) NPAR GO TO 40 1000 FORMAT('0***** FUNCTION VALUE DOES NOT SEEM TO DEPEND ON ANY ', +'OF THE',I3,' VARIABLE PARAMETERS',/15X ,'VERIFY THAT STEP SIZES', +' ARE BIG ENOUGH AND CHECK FCN LOGIC.',/1X,81(1H*)/1X,81(1H*)//) END SUBROUTINE PYSIMP C...This is the MINUIT routine SIMPLEX. CC PERFORMS A MINIMIZATION USING THE SIMPLEX METHOD OF NELDER CC AND MEAD (REF. -- COMP. J. 7,308 (1965)). COMMON /PYMINU/ XKIN(4),UKIN(4),WKIN(4),AIN(4),BIN(4), &MAXFIN,RELUP,RELERR,RELER2,FCNMAX COMMON 1/PYMINE/ ERP(30) ,ERN(30) 2/PYPARI/ X(15) ,XT(15) ,DIRIN(15) ,MAXINT ,NPAR 3/PYPARE/ U(30) ,WERR(30) ,MAXEXT ,NU 4/PYLIMI/ ALIM(30) ,BLIM(30) ,LCODE(30) ,LCORSP(30) ,LIMSET 5/PYVARI/ V(15,15) 7/PYFIX / IPFIX(15),XS(15) ,XTS(15) ,DIRINS(15) ,NPFIX 7/PYFIX2/ GRDS(15) ,G2S(15) ,GSTEPS(15),ABERFS(15) C/PYCASC/ Y(16) ,JH ,JL F/PYDERI/ GIN(30) ,GRD(15) ,G2(15) ,GSTEP(15) ,ABERF(15) G/PYSIMV/ P(15,16) ,PSTAR(15),PSTST(15) ,PBAR(15) ,PRHO(15) J/PYVART/ VT(15,15) COMMON 6/PYUNIT/ ISYSRD ,ISYSWR ,ISYSPU 8/PYTITL/ TITLE(13),DATE(2) ,ISW(7) ,NBLOCK 9/PYCONV/ EPSI ,APSI ,VTEST ,NSTEPQ ,NFCN ,NFCNMX A/PYCARD/ CWORD ,CWORD2 ,CWORD3 ,WORD7(7) B/PYMINI/ AMIN ,UP ,NEWMIN ,ITAUR ,SIGMA,EPSMAC DATA ALPHA,BETA,GAMMA,RHOMIN,RHOMAX / 1.0, 0.5, 2.0, 4.0, 8.0/ IF (NPAR .LE. 0) RETURN NPFN=NFCN NPARP1=NPAR+1 RHO1 = 1.0 + ALPHA RHO2 = RHO1 + ALPHA*GAMMA WG = 1.0/FLOAT(NPAR) IFLAG=4 WRITE(ISYSWR,100) EPSI DO 2 I= 1, NPAR IF (ISW(2) .GE. 1) DIRIN(I) = SQRT(V(I,I)*UP) IF (ABS(DIRIN(I)) .LT. 1.0E-10*ABS(X(I))) DIRIN(I)=1.0E-8*X(I) IF(ITAUR.LT. 1) V(I,I) = DIRIN(I)**2/UP 2 CONTINUE IF (ITAUR .LT. 1) ISW(2) = 1 C** CHOOSE THE INITIAL SIMPLEX USING SINGLE-PARAMETER SEARCHES 1 CONTINUE YNPP1 = AMIN JL = NPARP1 Y(NPARP1) = AMIN ABSMIN = AMIN DO 10 I= 1, NPAR AMING = AMIN PBAR(I) = X(I) BESTX = X(I) KG = 0 NS = 0 NF = 0 4 X(I) = BESTX + DIRIN(I) CALL PYINTO(X) CALL PYSIGM(NPAR,GIN, F, U, 4) NFCN = NFCN + 1 IF (F .LE. AMING) GO TO 6 C FAILURE IF (KG .EQ. 1) GO TO 8 KG = -1 NF = NF + 1 DIRIN(I) = DIRIN(I) * (-0.4) IF (NF .LT. 3) GO TO 4 NS = 6 C SUCCESS 6 BESTX = X(I) DIRIN(I) = DIRIN(I) * 3.0 AMING = F KG = 1 NS = NS + 1 IF (NS .LT. 6) GO TO 4 C LOCAL MINIMUM FOUND IN ITH DIRECTION 8 Y(I) = AMING IF (AMING .LT. ABSMIN) JL = I IF (AMING .LT. ABSMIN) ABSMIN = AMING X(I) = BESTX DO 9 K= 1, NPAR 9 P(K,I) = X(K) 10 CONTINUE JH = NPARP1 AMIN=Y(JL) CALL PYRAZZ(YNPP1,PBAR) DO 20 I= 1, NPAR 20 X(I) = P(I,JL) CALL PYINTO(X) DO 30 I=1,NPAR 30 IF(ABS(DIRIN(I)).LE.ABS(EPSMAC*X(I))) DIRIN(I)=4.*EPSMAC*X(I) IF (ISW(5) .GE. 1) CALL PYMPRI(0,AMIN) SIGMA = SIGMA * 10. SIG2 = SIGMA NCYCL=0 C . . . . . START MAIN LOOP 50 CONTINUE C...Change in SIMPLX; error redefined for second call to PYSIMP. UP=RELUP*ABS(AMIN) EPSI=RELERR*UP IF (SIG2 .LT. EPSI .AND. SIGMA.LT.EPSI) GO TO 76 SIG2 = SIGMA IF ((NFCN-NPFN) .GT. NFCNMX) GO TO 78 C CALCULATE NEW POINT * BY REFLECTION DO 60 I= 1, NPAR PB = 0. DO 59 J= 1, NPARP1 59 PB = PB + WG * P(I,J) PBAR(I) = PB - WG * P(I,JH) 60 PSTAR(I)=(1.+ALPHA)*PBAR(I)-ALPHA*P(I,JH) CALL PYINTO(PSTAR) CALL PYSIGM(NPAR,GIN,YSTAR,U,4) NFCN=NFCN+1 IF(YSTAR.GE.AMIN) GO TO 70 C POINT * BETTER THAN JL, CALCULATE NEW POINT ** DO 61 I=1,NPAR 61 PSTST(I)=GAMMA*PSTAR(I)+(1.-GAMMA)*PBAR(I) CALL PYINTO(PSTST) CALL PYSIGM(NPAR,GIN,YSTST,U,4) NFCN=NFCN+1 C TRY A PARABOLA THROUGH PH, PSTAR, PSTST. MIN = PRHO Y1 = (YSTAR-Y(JH)) * RHO2 Y2 = (YSTST-Y(JH)) * RHO1 RHO = 0.5 * (RHO2*Y1 -RHO1*Y2) / (Y1 -Y2) IF (RHO .LT. RHOMIN) GO TO 66 IF (RHO .GT. RHOMAX) RHO = RHOMAX DO 64 I= 1, NPAR 64 PRHO(I) = RHO*PBAR(I) + (1.0-RHO)*P(I,JH) CALL PYINTO(PRHO) CALL PYSIGM(NPAR,GIN,YRHO, U,4) NFCN = NFCN + 1 IF (YRHO .LT. Y(JL) .AND. YRHO .LT. YSTST) GO TO 65 IF (YSTST .LT. Y(JL)) GO TO 67 IF (YRHO .GT. Y(JL)) GO TO 66 C ACCEPT MINIMUM POINT OF PARABOLA, PRHO 65 CALL PYRAZZ (YRHO,PRHO) GO TO 68 66 IF (YSTST .LT. Y(JL)) GO TO 67 CALL PYRAZZ(YSTAR,PSTAR) GO TO 68 67 CALL PYRAZZ(YSTST,PSTST) 68 NCYCL=NCYCL+1 IF (ISW(5) .LT. 2) GO TO 50 IF (ISW(5) .GE. 3 .OR. MOD(NCYCL, 10) .EQ. 0) CALL PYMPRI(0,AMIN) GO TO 50 C POINT * IS NOT AS GOOD AS JL 70 IF (YSTAR .GE. Y(JH)) GO TO 73 JHOLD = JH CALL PYRAZZ(YSTAR,PSTAR) IF (JHOLD .NE. JH) GO TO 50 C CALCULATE NEW POINT ** 73 DO 74 I=1,NPAR 74 PSTST(I)=BETA*P(I,JH)+(1.-BETA)*PBAR(I) CALL PYINTO (PSTST) CALL PYSIGM(NPAR,GIN,YSTST,U,4) NFCN=NFCN+1 IF(YSTST.GT.Y(JH)) GO TO 1 C POINT ** IS BETTER THAN JH IF (YSTST .LT. AMIN) GO TO 67 CALL PYRAZZ(YSTST,PSTST) GO TO 50 C . . . . . . END MAIN LOOP 76 WRITE(ISYSWR,120) GO TO 80 78 WRITE(ISYSWR,130) ISW(1) = 1 80 DO 82 I=1,NPAR PB = 0. DO 81 J=1,NPARP1 81 PB = PB + WG * P(I,J) 82 PBAR(I) = PB - WG * P(I,JH) CALL PYINTO(PBAR) CALL PYSIGM(NPAR,GIN,YPBAR,U,IFLAG) NFCN=NFCN+1 IF (YPBAR .LT. AMIN) CALL PYRAZZ(YPBAR,PBAR) CALL PYINTO(X) IF (NFCNMX+NPFN-NFCN .LT. 3*NPAR) GO TO 90 IF (SIGMA .GT. 2.0*EPSI) GO TO 1 90 CALL PYMPRI(1-ITAUR, AMIN) RETURN 100 FORMAT(' START SIMPLEX MINIMIZATION ',8X ,'CON', +'VERGENCE CRITERION -- ESTIMATED DISTANCE TO MINIMUM (EDM) .LT.', +E10.2 ) 120 FORMAT(1H ,'SIMPLEX MINIMIZATION HAS CONVERGED') 130 FORMAT(1H ,'SIMPLEX TERMINATES WITHOUT CONVERGENCE') END C*********************************************************************** C C The following routines for parton cascades were made together with C M. Bengtsson and T. Sjostrand, Lund Univ. They are modifications C of routines in PYTHIA 4.8, LU TP 87-3. For details of parton showers C in DIS see M. Bengtsson, T. Sjostrand, LU TP 87-10 and M. Bengtsson, C G. Ingelman, T. Sjostrand, DESY preprint in preparation. C C ********************************************************************** SUBROUTINE CASCAD(ICALL) COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U COMMON /LUJETS/ N,K(2000,2),P(2000,5) COMMON /LUDAT1/ MST(40),PAR(80) COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),DE,DPZ,DPT,PB(5),PHIR COMMON /PYPARA/ IPY(80),PYPAR(80),PYVAR(80) COMMON /PYPROC/ ISUB,KFL(3,2),XPY(2),SH,TH,UH,Q2PY,XSEC(0:40) COMMON /PYINT1/ XQPY(2,-6:6),DSIG(-6:6,-6:6,5),FSIG(10,10,3) DOUBLE PRECISION DTHETA,DPHI,DBETA,DE,DPZ,DPT DOUBLE PRECISION DPQ2,DPB(3),DPA(3),DCTHET DIMENSION KS(9,2),PS(9,5),ROBO(5),XPQ(-6:6) IF(ICALL.EQ.0) THEN C...Initialize cascade for each event, save event record in overall cms. MG=1 KS(1,1)=40000 KS(2,1)=40000 KS(3,1)=0 KS(4,1)=0 KS(5,1)=40003 KS(7,1)=40005 KS(1,2)=K(1,2) KS(2,2)=K(2,2) KS(3,2)=0 KS(4,2)=0 KS(5,2)=K(3,2) KS(7,2)=K(4,2) DO 10 J=1,5 PS(1,J)=P(1,J) PS(2,J)=P(2,J) PS(3,J)=0. PS(4,J)=0. PS(5,J)=P(3,J) 10 PS(7,J)=P(4,J) RETURN ENDIF C...Apply parton cascade on QPM event. C...Save incoming and outgoing quark as well as scattered lepton. KS(6,1)=40004 KS(6,2)=ISIGN(500+IABS(LST(25)),LST(25)) KS(8,1)=40006 KS(8,2)=K(5,2) KS(9,1)=80005 KS(9,2)=K(4,2) DO 110 J=1,5 PS(6,J)=0. PS(8,J)=P(5,J) 110 PS(9,J)=P(4,J) XR=X DPQ2=DBLE(Q2) PMA1=0. PS(6,5)=PMA1 PMA2=PS(8,5) DPB(1)=0.5D0*(DPQ2*(1D0/XR-1D0)+DBLE(PS(1,5))**2- &ULMASS(0,IABS(KS(7,2)))**2)/(PS(1,4)+PS(2,4)) DPB(2)=DSQRT(DPB(1)**2+DPQ2) DCTHET=(DBLE(PS(2,4))*DPB(1)-DPQ2/(2D0*XR))/(DBLE(PS(2,3))* &DPB(2)) DPA(1)=(DPB(2)*DCTHET)**2-DPB(1)**2 DPA(2)=DPQ2-DBLE(PMA1)**2+DBLE(PMA2)**2 PS(6,4)=-(DPA(2)*DPB(1)-DPB(2)*DCTHET*DSQRT(DPA(2)**2+4D0* &DBLE(PMA1)**2*DPA(1)))/(2D0*DPA(1)) PS(6,3)=-SQRT((PS(6,4)+PMA1)*(PS(6,4)-PMA1)) C...Partons with colour information in hadronic cms frame. DO 120 I=10,26 K(I,1)=80000 K(I,2)=0 DO 120 J=1,5 120 P(I,J)=0. NS=20 K(NS+1,1)=40003 K(NS+1,2)=K(3,2) K(NS+2,1)=70021 K(NS+3,2)=KS(6,2) DO 130 J=1,5 130 P(NS+1,J)=P(3,J) K(NS+3,1)=2+20000 IF(MG.EQ.1) THEN PW2=W2 DPA(3)=DSQRT(DPA(2)**2+4D0*DPQ2*DBLE(PMA1)**2) DPB(1)=(1D0/DBLE(XR)-2D0)*DPQ2/(2D0*SQRT(PW2)) DPB(2)=DSQRT(DPB(1)**2+DPQ2) P(NS+3,4)=(DPB(2)*DPA(3)-DPB(1)*DPA(2))/(2D0*DPQ2) P(NS+3,3)=-P(NS+3,4) ENDIF P(NS+3,5)=0. K(NS+4,1)=70003+NS P(NS+4,3)=NS+5 P(NS+4,4)=NS+5 K(NS+5,1)=8 K(NS+5,2)=KS(8,2) K(NS+6,1)=70005+NS DO 140 J=1,4 IF(MG.EQ.1) THEN P(NS+5,J)=P(NS+1,J)+P(NS+3,J) ELSE P(NS+5,J)=P(5,J) P(NS+3,J)=P(NS+5,J)-P(NS+1,J) ENDIF 140 CONTINUE P(NS+5,5)=PMA2 P(NS+6,1)=NS+3 P(NS+6,2)=NS+3 N=NS+6 C...Copy saved record in overall cms to line 1 through 9. DO 150 I=1,9 K(I,1)=KS(I,1) K(I,2)=KS(I,2) DO 150 J=1,5 150 P(I,J)=PS(I,J) C...Scale for bremsstrahlung etc. Q2PY=Q2 IPY(40)=8 IPY(47)=N C...Save quantities for later use. XPY(1)=1. XPY(2)=XR CALL PYSTFU(K(2,2),XR,Q2,XPQ) DO 160 IFL=-6,6 160 XQPY(2,IFL)=XPQ(IFL) IF(LST(1).EQ.1) THEN ISUB=39 IPY(11)=1 ELSEIF(LST(1).EQ.3) THEN ISUB=39 IPY(11)=2 ELSEIF(LST(1).EQ.4) THEN ISUB=39 IPY(11)=3 ELSEIF(LST(1).EQ.2) THEN ISUB=40 ENDIF IF(ISUB.EQ.39.AND.IPY(11).EQ.1) THEN KFL(2,1)=1 ELSEIF(ISUB.EQ.39.AND.IPY(11).EQ.2) THEN KFL(2,1)=2 ELSEIF(ISUB.EQ.39.AND.IPY(11).EQ.3) THEN KFL(2,1)=5 ELSEIF(ISUB.EQ.40) THEN KFL(2,1)=-3 ENDIF IFL1=MOD(K(6,2),500) IFL2=MOD(K(8,2),500) KFL(2,2)=ISIGN(500+IABS(IFL1),2*IFL1+1) KFL(1,1)=KFL(2,1) KFL(1,2)=KFL(2,2) IF(ISUB.EQ.39) KFL(3,1)=K(1,2) IF(ISUB.EQ.40) KFL(3,1)=K(1,2)+ISIGN(1,K(1,2)) KFL(3,2)=ISIGN(500+IABS(IFL2),2*IFL2+1) PYVAR(2)=(P(1,4)+P(2,4))**2 PYVAR(1)=SQRT(PYVAR(2)) PYVAR(3)=P(1,5) PYVAR(4)=P(2,5) PYVAR(5)=P(1,3) IPY(41)=K(1,2) IPY(42)=K(2,2) IPY(48)=0 C...Generate timelike parton shower (if required) IF(IPY(13).EQ.1) THEN IF(LST(17).EQ.1) THEN QMAX=MIN(SQRT(PYPAR(25)*Q2),P(25,4)) ELSE QMAX=MIN(SQRT(PYPAR(25)*W2),P(25,4)) ENDIF CALL LUSHOW(25,0,QMAX) ENDIF IT=25 IF(N.GE.27) IT=27 NS=N C...Generate spacelike parton shower (if required) IPU1=0 IPU2=23 IF(XPY(2)*(1.+(P(IT,5)**2+PYPAR(22))/P(21,5)**2).GT.0.999) THEN LST(21)=47 RETURN ENDIF IF(IPY(14).GE.1) THEN CALL PYSSPA(IPU1,IPU2) ELSEIF(IPY(13).GE.1) THEN I1=NS+1 I2=NS+3 DO 230 J=1,5 P(NS+1,J)=0. P(NS+2,J)=0. P(NS+3,J)=0. P(NS+4,J)=0. 230 P(I1,J)=P(21,J) K(NS+1,1)=20000 K(NS+3,1)=20000 K(NS+2,1)=70001+NS K(NS+4,1)=70003+NS K(I1,2)=KFL(2,1) K(I2,2)=KFL(2,2) P(I2,3)=(P(IT,5)**2+Q2)*(P(21,4)-P(21,3))/(2.*Q2) P(I2,4)=-P(I2,3) P(I2+1,3)=23 P(I2+1,4)=23 P(24,1)=I2 P(24,2)=I2 IPU1=0 IPU2=I2 N=N+4 ENDIF C...Rotate and boost outgoing parton shower IF(N.GT.30) THEN K(N+1,1)=0 DO 210 J=1,4 210 P(N+1,J)=P(NS+1,J)+P(NS+3,J) IF(P(N+1,4).LE.1.01*P(IT,5)) THEN LST(21)=50 RETURN ENDIF MST(1)=25 MST(2)=NS ROBO(1)=ULANGL(P(IT,3),SQRT(P(IT,1)**2+P(IT,2)**2)) ROBO(2)=ULANGL(P(IT,1),P(IT,2)) CALL LUROBO(0.,-ROBO(2),0.,0.,0.) CALL LUROBO(-ROBO(1),0.,0.,0.,0.) ROBO(5)=-(P(IT,3)*P(IT,4)-P(N+1,4)*SQRT(P(N+1,4)**2- & P(IT,4)**2+P(IT,3)**2))/(P(IT,3)**2+P(N+1,4)**2) CALL LUROBO(0.,0.,0.,0.,ROBO(5)) ROBO(1)=ULANGL(P(N+1,3),SQRT(P(N+1,1)**2+P(N+1,2)**2)) ROBO(2)=ULANGL(P(N+1,1),P(N+1,2)) CALL LUROBO(ROBO(1),ROBO(2),0.,0.,0.) MST(1)=0 MST(2)=0 ENDIF Q2PY=Q2 C...Hadron remnant and primordial kt IPY(47)=N CALL PYREMN(IPU1,IPU2) IF(IPY(48).EQ.1) THEN LST(21)=48 RETURN ENDIF C...Transform line 1,2 and 5-8 to hadronic cms frame. MST(1)=0 MST(2)=8 K(3,1)=K(3,1)+40000 K(4,1)=K(4,1)+40000 CALL DUROBO(0.D0,0.D0,-DBETA(2,1),-DBETA(2,2),-DBETA(2,3)) CALL DUROBO(-DTHETA(2),0.D0,0.D0,0.D0,0.D0) MST(2)=0 K(3,1)=K(3,1)-40000 K(4,1)=K(4,1)-40000 C...Activate line with scattered lepton. K(9,1)=5 IF(LST(4).EQ.0) K(9,1)=40005 C...Deactivate exchanged boson in line 21. K(21,1)=MOD(K(21,1),10000)+80000 C...Rearrange partons along strings MST(21)=2-IPY(34) IF(IPY(34).EQ.0) MST(21)=3 MST(26)=0 CALL LUPREP IF(MST(26).NE.0) THEN WRITE(6,*) ' Error in LUPREP! MST(25), MST(26) = ',MST(25), & MST(26) LST(21)=49 RETURN ENDIF RETURN END C*********************************************************************** SUBROUTINE PYSSPA(IPU1,IPU2) C...NEW X REDEFINITION C...GENERATES SPACELIKE PARTON SHOWERS COMMON/LUJETS/N,K(2000,2),P(2000,5) COMMON/LUDAT1/MST(40),PAR(80) COMMON/LUDAT2/KTYP(120),PMAS(120),PWID(60),KFR(80),CFR(40) COMMON/PYPARA/IPY(80),PYPAR(80),PYVAR(80) COMMON/PYPROC/ISUB,KFL(3,2),X(2),SH,TH,UH,Q2,XSEC(0:40) COMMON/PYINT1/XQ(2,-6:6),DSIG(-6:6,-6:6,5),FSIG(10,10,3) DIMENSION IFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVS(2),ROBO(5), &XFS(2,-6:6),XFA(-6:6),XFB(-6:6),WTAP(-6:6),WTSF(-6:6) DOUBLE PRECISION DQ2(3),DSH,DSHZ,DSHR,DPLCM,DPC(3),DPD(4),DMS, &DMSMA,DPT2,DPB(4),DBE1(4),DBE2(4),DBEP,DGABEP,DPQ(4),DPQS(2), &DM2,DQ2B C-GI &DQ23,DPH(4),DM2,DQ2B,DQM2 C...COMMON CONSTANTS, SET UP INITIAL VALUES ILEP=0 IF(IPU1.EQ.0) ILEP=1 IF(IPU2.EQ.0) ILEP=2 Q2E=Q2 CGI IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.3) Q2E=Q2E/PYPAR(26) IF(ISUB.EQ.27) Q2E=PMAS(2)**2 IF(ISUB.EQ.28) Q2E=PMAS(3)**2 TMAX=ALOG(PYPAR(26)*PYPAR(27)*Q2E/PYPAR(21)**2) IF(ILEP.GE.1) THEN SH=P(25,5)**2 IF(N.GE.27) SH=P(27,5)**2 Q2E=Q2 IF(IPY(77).EQ.2) Q2E=Q2*(1.-X(3-ILEP))/X(3-ILEP) Q2E=MAX(PYPAR(21)**2,MIN(Q2E,(0.95/X(3-ILEP)-1.)*Q2-SH, & Q2/2.+SH)) TMAX=ALOG(Q2E/PYPAR(21)**2) ENDIF IF(PYPAR(26)*Q2E.LT.MAX(PYPAR(22),2.*PYPAR(21)**2).OR. &TMAX.LT.0.2) RETURN IF(ILEP.EQ.0) XE0=2.*PYPAR(23)/PYVAR(1) B0=(33.-2.*IPY(8))/6. NS=N MST(2)=0 100 N=NS IF(ILEP.GE.1) THEN NQ=IPU2-2 IF(ILEP.EQ.2) NQ=IPU1+2 DPQS(1)=DBLE(P(NQ,3)) DPQS(2)=DBLE(P(NQ,4)) XBMIN=X(3-ILEP)*MAX(0.5,SH/Q2) CALL PYSTFU(IPY(43-ILEP),XBMIN,Q2,XFB) DO 110 IFL=-6,6 110 XQ(3-ILEP,IFL)=XFB(IFL) ENDIF DO 120 JT=1,2 IFLS(JT)=MOD(KFL(2,JT),500) IFLS(JT+2)=IFLS(JT) XS(JT)=X(JT) ZS(JT)=1. IF(ILEP.EQ.0) Q2S(JT)=PYPAR(26)*Q2E TEVS(JT)=TMAX DO 120 IFL=-6,6 120 XFS(JT,IFL)=XQ(JT,IFL) IF(ILEP.GE.1) THEN Q2S(ILEP)=P(NQ,5)**2 DQ2(ILEP)=Q2S(ILEP) Q2S(3-ILEP)=Q2E ENDIF DSH=SH IHFC=0 IHFX=0 C...PICK UP LEG WITH HIGHEST VIRTUALITY 130 N=N+2 JT=1 IF((N.GT.NS+2.AND.Q2S(2).GT.Q2S(1).AND.ILEP.EQ.0).OR.ILEP.EQ.1) &JT=2 JR=3-JT IFLB=IFLS(JT) XB=XS(JT) IF(ILEP.GE.1.AND.N.EQ.NS+2) XB=XS(JT)*MAX(SH/Q2,0.5) DO 140 IFL=-6,6 140 XFB(IFL)=XFS(JT,IFL) Q2B=Q2S(JT) TEVB=TEVS(JT) IF(IPY(14).GE.9.AND.N.GT.NS+4) THEN Q2B=0.5*(1./ZS(JT)+1.)*Q2S(JT)+0.5*(1./ZS(JT)-1.)*(Q2S(3-JT)- & SNGL(DSH)+SQRT((SNGL(DSH)+Q2S(1)+Q2S(2))**2+8.*Q2S(1)*Q2S(2)* & ZS(JT)/(1.-ZS(JT)))) TEVB=ALOG(PYPAR(27)*Q2B/PYPAR(21)**2) ENDIF IF(ILEP.EQ.0) THEN DSHR=2.*DSQRT(DSH) DSHZ=DSH/DBLE(ZS(JT)) ELSEIF(ILEP.GE.1) THEN DSHZ=DSH IF(N.GT.NS+4) DSHZ=(DSH+DQ2(JR)-DQ2(JT))/ZS(JT)-DQ2(JR)+ & PYPAR(22) DPD(2)=DSHZ+DQ2(JR)+DBLE(PYPAR(22)) C...CHECK IF QUARK PAIR CREATION ONLY POSSIBILITY IF(DQ2(JR).LT.4.*ULMASS(2,IABS(IFLB))**2) THEN DM2=ULMASS(2,IABS(IFLB))**2 DPC(1)=DQ2(JR)*(DBLE(PYPAR(22))+DM2)**2 DPC(2)=DPD(2)*(DPD(2)-2D0*PYPAR(22))*(PYPAR(22)+DM2) DPC(3)=PYPAR(22)*(DPD(2)-2D0*PYPAR(22))**2 XE0=1D0-(DPC(2)-DSQRT(DPC(2)**2-4D0*DPC(1)*DPC(3)))/ & (2D0*DPC(1)) ELSE XE0=1D0-(DPD(2)-2D0*DBLE(PYPAR(22)))*(DPD(2)-DSQRT(DPD(2)**2- & 4D0*DQ2(JR)*DBLE(PYPAR(22))))/(2D0*DQ2(JR)*DBLE(PYPAR(22))) ENDIF ENDIF 145 XE=MAX(XE0,XB*(1./(1.-PYPAR(24))-1.)) IF(XB+XE.GE.0.999) THEN Q2B=0. GOTO 210 ENDIF C...CALCULATE ALTARELLI-PARISI AND STRUCTURE FUNCTION WEIGHTS DO 150 IFL=-6,6 WTAP(IFL)=0. 150 WTSF(IFL)=0. IF(IFLB.EQ.0) THEN WTAPQ=16.*(1.-SQRT(XB+XE))/(3.*SQRT(XB)) DO 160 IFL=-IPY(8),IPY(8) IF(IFL.EQ.0) WTAP(IFL)=6.*ALOG((1.-XB)/XE) 160 IF(IFL.NE.0) WTAP(IFL)=WTAPQ ELSE WTAP(0)=0.5*XB*(1./(XB+XE)-1.) WTAP(IFLB)=8.*ALOG((1.-XB)*(XB+XE)/XE)/3. ENDIF 170 WTSUM=0. IF(IHFC.EQ.0) THEN DO 180 IFL=-IPY(8),IPY(8) WTSF(IFL)=XFB(IFL)/MAX(1E-10,XFB(IFLB)) 180 WTSUM=WTSUM+WTAP(IFL)*WTSF(IFL) IF(IABS(IFLB).GE.4.AND.WTSUM.GT.1E3) THEN IHFX=1 DO 185 IFL=-IPY(8),IPY(8) 185 WTSF(IFL)=WTSF(IFL)*1E3/WTSUM WTSUM=1E3 ENDIF ENDIF C...CHOOSE NEW T AND FLAVOUR 190 IF(IPY(14).LE.6.OR.IPY(14).GE.9) THEN TEVB=TEVB*RLU(0)**(B0/MAX(0.0001,WTSUM)) ELSE TEVB=TEVB*RLU(0)**(B0/MAX(0.0001,5.*WTSUM)) ENDIF Q2REF=PYPAR(21)**2*EXP(TEVB)/PYPAR(27) Q2B=Q2REF/PYPAR(27) DQ2B=Q2B IF(ILEP.GE.1) THEN DSHZ=DSH IF(N.GT.NS+4) DSHZ=(DSH+DQ2(JR)-DQ2(JT))/DBLE(ZS(JT))-DQ2(JR)+ & DQ2B ENDIF IF(Q2B.LT.PYPAR(22)) THEN Q2B=0. ELSE WTRAN=RLU(0)*WTSUM IFLA=-IPY(8)-1 200 IFLA=IFLA+1 WTRAN=WTRAN-WTAP(IFLA)*WTSF(IFLA) IF(IFLA.LT.IPY(8).AND.WTRAN.GT.0.) GOTO 200 C...CHOOSE Z VALUE AND CORRECTIVE WEIGHT IF(IFLB.EQ.0.AND.IFLA.EQ.0) THEN Z=1./(1.+((1.-XB)/XB)*(XE/(1.-XB))**RLU(0)) WTZ=(1.-Z*(1.-Z))**2 ELSEIF(IFLB.EQ.0) THEN Z=XB/(1.-RLU(0)*(1.-SQRT(XB+XE)))**2 WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z) ELSEIF(IFLA.EQ.0) THEN Z=XB*(1.+RLU(0)*(1./(XB+XE)-1.)) WTZ=1.-2.*Z*(1.-Z) ELSE Z=1.-(1.-XB)*(XE/((XB+XE)*(1.-XB)))**RLU(0) WTZ=0.5*(1.+Z**2) ENDIF C...REWEIGHT FIRST LEG BECAUSE OF MODIFIED XB OR CHECK PHASE SPACE IF(ILEP.GE.1.AND.N.EQ.NS+2) THEN XBNEW=X(JT)*(1.+(DSH-Q2B)/DQ2(JR)) IF(XBNEW.GT.MIN(Z,0.999)) GOTO 190 XB=XBNEW ENDIF C...SUM UP SOFT GLUON EMISSION AS EFFECTIVE Z SHIFT IF(IPY(15).GE.1) THEN RSOFT=6. IF(IFLB.NE.0) RSOFT=8./3. Z=Z*(TEVB/TEVS(JT))**(RSOFT*XE/((XB+XE)*B0)) IF(Z.LE.XB) GOTO 190 ENDIF C...CHECK IF HEAVY FLAVOUR BELOW THRESHOLD IHFT=0 IF(ILEP.GE.1.AND.IABS(IFLB).GE.4.AND.(XFB(IFLB).LT.1E-10.OR. & Q2B.LT.5.*ULMASS(2,IABS(IFLB))**2)) THEN IHFT=1 IFLA=0 ENDIF C...FOR LEPTOPRODUCTION, CHECK Z AGAINST NEW LIMIT IF(ILEP.GE.1) THEN DPD(2)=DSHZ+DQ2(JR)+DQ2B DM2=ULMASS(2,IABS(IFLA-IFLB))**2 DPC(1)=DQ2(JR)*(DQ2B+DM2)**2 DPC(2)=DPD(2)*(DPD(2)-2D0*DQ2B)*(DQ2B+DM2) DPC(3)=DQ2B*(DPD(2)-2D0*DQ2B)**2 ZU=(DPC(2)-DSQRT(DPC(2)**2-4D0*DPC(1)*DPC(3)))/(2D0*DPC(1)) IF(Z.GE.ZU) GOTO 190 ENDIF C...OPTION WITH EVOLUTION IN KT2=(1-Z)Q2: IF(IPY(14).GE.5.AND.IPY(14).LE.6.AND.N.LE.NS+4) THEN C...CHECK THAT (Q2)LAST BRANCHING < (Q2)HARD IF(Q2B/(1.-Z).GT.PYPAR(26)*Q2) GOTO 190 ELSEIF(IPY(14).GE.3.AND.IPY(14).LE.6.AND.N.GE.NS+6) THEN C...CHECK THAT Z,Q2 COMBINATION IS KINEMATICALLY ALLOWED Q2MAX=0.5*(1./ZS(JT)+1.)*DQ2(JT)+0.5*(1./ZS(JT)-1.)* & (DQ2(3-JT)-DSH+SQRT((DSH+DQ2(1)+DQ2(2))**2+8.*DQ2(1)*DQ2(2)* & ZS(JT)/(1.-ZS(JT)))) IF(Q2B/(1.-Z).GE.Q2MAX) GOTO 190 ELSEIF(IPY(14).EQ.7.OR.IPY(14).EQ.8) THEN C...OPTION WITH ALPHAS((1-Z)Q2): DEMAND KT2 > CUTOFF, REWEIGHT IF((1.-Z)*Q2B.LT.PYPAR(22)) GOTO 190 ALPRAT=TEVB/(TEVB+ALOG(1.-Z)) IF(ALPRAT.LT.5.*RLU(0)) GOTO 190 IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5. ENDIF C...WEIGHTING WITH NEW STRUCTURE FUNCTIONS CALL PYSTFU(IPY(40+JT),XB,Q2REF,XFB) XA=XB/Z CALL PYSTFU(IPY(40+JT),XA,Q2REF,XFA) IF(IHFT.EQ.1.OR.IHFX.EQ.1) THEN IF(XFA(IFLA).LT.1E-10) IHFC=1 GOTO 210 ELSEIF(XFB(IFLB).LT.1E-20) THEN GOTO 100 ENDIF IF(WTZ*XFA(IFLA)/XFB(IFLB).LT.RLU(0)*WTSF(IFLA)) THEN IF(ILEP.GE.1.AND.N.EQ.NS+2) GOTO 145 GOTO 170 ENDIF ENDIF 210 IF(N.EQ.NS+4-2*MIN(1,ILEP)) THEN C...DEFINE TWO HARD SCATTERERS IN THEIR CM-FRAME DQ2(JT)=Q2B IF(IPY(14).GE.3.AND.IPY(14).LE.6) DQ2(JT)=Q2B/(1.-Z) IF(ILEP.EQ.0) THEN DPLCM=DSQRT((DSH+DQ2(1)+DQ2(2))**2-4.*DQ2(1)*DQ2(2))/DSHR DO 220 JR=1,2 I=NS+2*JR-1 IPO=19+2*JR K(I,1)=20000 K(I,2)=ISIGN(500+IABS(IFLS(JR+2)),2*IFLS(JR+2)+1) P(I,1)=0. P(I,2)=0. P(I,3)=DPLCM*(-1)**(JR+1) P(I,4)=(DSH+DQ2(3-JR)-DQ2(JR))/DSHR P(I,5)=-SQRT(SNGL(DQ2(JR))) K(I+1,1)=70000+I K(I+1,2)=K(IPO+1,2) P(I+1,1)=0. P(I+1,2)=0. P(I+1,3)=IPO P(I+1,4)=IPO P(I+1,5)=0. P(IPO+1,1)=I P(IPO+1,2)=I 220 CONTINUE ELSE C..LEPTOPRODUCTION EVENTS: BOSON AND HADRON REST FRAME I1=NS+2*ILEP-1 I2=NS-2*ILEP+5 DO 230 J=1,5 P(NS+1,J)=0. P(NS+2,J)=0. P(NS+3,J)=0. P(NS+4,J)=0. 230 P(I1,J)=P(NQ,J) K(NS+1,1)=20000 K(NS+3,1)=20000 K(NS+2,1)=70001+NS K(NS+4,1)=70003+NS K(I1,2)=KFL(2,ILEP) K(I2,2)=KFL(2,3-ILEP) DPD(1)=DSH+DQ2(1)+DQ2(2) DPD(3)=(3-2*ILEP)*DSQRT(DPD(1)**2-4D0*DQ2(1)*DQ2(2)) P(I2,3)=(DPQS(2)*DPD(3)-DPQS(1)*DPD(1))/ & (2D0*DQ2(JR)) P(I2,4)=(DPQS(1)*DPD(3)-DPQS(2)*DPD(1))/ & (2D0*DQ2(JR)) P(I2,5)=-SQRT(SNGL(DQ2(3-ILEP))) P(I2+1,3)=MAX(IPU1,IPU2) P(I2+1,4)=MAX(IPU1,IPU2) P(26-2*ILEP,1)=I2 P(26-2*ILEP,2)=I2 N=N+2 ENDIF ELSEIF(N.GT.NS+4) THEN C...FIND MAXIMUM ALLOWED MASS OF TIMELIKE PARTON DQ2(3)=Q2B IF(IPY(14).GE.3.AND.IPY(14).LE.6) DQ2(3)=Q2B/(1.-Z) DPC(1)=P(IS(1),4) DPC(2)=P(IS(2),4) DPC(3)=0.5*(ABS(P(IS(1),3))+ABS(P(IS(2),3))) DPD(1)=DSH+DQ2(JR)+DQ2(JT) DPD(2)=DSHZ+DQ2(JR)+DQ2(3) DPD(3)=DSQRT(DPD(1)**2-4.*DQ2(JR)*DQ2(JT)) DPD(4)=DSQRT(DPD(2)**2-4.*DQ2(JR)*DQ2(3)) IKIN=0 IF((Q2S(JR).GE.0.5*PYPAR(22).AND.DPD(1)-DPD(3).GE.1D-10*DPD(1)) & .OR.ILEP.GE.1) IKIN=1 IF(IKIN.EQ.0) DMSMA=(DQ2(JT)/DBLE(ZS(JT))-DQ2(3))*(DSH/ & (DSH+DQ2(JT))-DSH/(DSHZ+DQ2(3))) IF(IKIN.EQ.1) DMSMA=(DPD(1)*DPD(2)-DPD(3)*DPD(4))/(2.* & DQ2(JR))-DQ2(JT)-DQ2(3) C...GENERATE TIMELIKE PARTON SHOWER (IF REQUIRED) IT=N-1 K(IT,1)=0 K(IT,2)=ISIGN(500+IABS(IFLB-IFLS(JT+2)),2*(IFLB-IFLS(JT+2))+1) P(IT,5)=ULMASS(0,K(IT,2)) IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100 P(IT,2)=0. K(IT+1,1)=70000+IT K(IT+1,2)=K(IS(JT)+1,2) DO 240 J=1,5 240 P(IT+1,J)=0. IF(MOD(IPY(14),2).EQ.0) THEN P(IT,1)=0. IF(ILEP.EQ.0) P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR IF(ILEP.GE.1) P(IT,4)=0.5*(P(IS(JT),3)*DPD(2)+ & DPQS(1)*(DQ2(JT)+DQ2(3)+P(IT,5)**2))/(P(IS(JT),3)*DPQS(2)- & P(IS(JT),4)*DPQS(1))-DPC(JT) P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2) CALL LUSHOW(IT,0,SQRT(MIN(SNGL(DMSMA),PYPAR(25)*Q2))) IF(N.GE.IT+2) P(IT,5)=P(IT+2,5) ENDIF C...RECONSTRUCT KINEMATICS OF BRANCHING: TIMELIKE PARTON SHOWER DMS=P(IT,5)**2 IF(IKIN.EQ.0.AND.ILEP.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/ & (DSH+DQ2(JT)) IF(IKIN.EQ.1.AND.ILEP.EQ.0) DPT2=(DMSMA-DMS)*(0.5*DPD(1)* & DPD(2)+0.5*DPD(3)*DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/ & (4.*DSH*DPC(3)**2) IF(IKIN.EQ.1.AND.ILEP.GE.1) DPT2=(DMSMA-DMS)*(0.5*DPD(1)* & DPD(2)+0.5*DPD(3)*DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/ & DPD(3)**2 IF(DPT2.LT.0.) GOTO 100 K(IT,1)=K(IT,1)+N+1 P(IT,1)=SQRT(SNGL(DPT2)) IF(ILEP.EQ.0) THEN DPB(1)=(0.5*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/ & DSHR)/DPC(3)-DPC(3) P(IT,3)=DPB(1)*(-1)**(JT+1) P(IT,4)=(DSHZ-DSH-DMS)/DSHR ELSE DPC(3)=DQ2(JT)+DQ2(3)+DMS DPB(2)=DPQS(2)*DBLE(P(IS(JT),3))-DPQS(1)*DPC(JT) DPB(1)=0.5D0*(DPC(JT)*DPD(2)+DPQS(2)*DPC(3))/DPB(2)- & DBLE(P(IS(JT),3)) P(IT,3)=DPB(1) P(IT,4)=0.5D0*(DBLE(P(IS(JT),3))*DPD(2)+ & DPQS(1)*DPC(3))/DPB(2)-DPC(JT) ENDIF IF(N.GE.IT+2) THEN MST(1)=IT+2 DPB(1)=DSQRT(DPB(1)**2+DPT2) DPB(2)=DSQRT(DPB(1)**2+DMS) DPB(3)=P(IT+2,3) DPB(4)=DSQRT(DPB(3)**2+DMS) BEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)* & DPB(1)) CALL LUROBO(0.,0.,0.,0.,BEZ) THE=ULANGL(P(IT,3),P(IT,1)) CALL LUROBO(THE,0.,0.,0.,0.) ENDIF C...RECONSTRUCT KINEMATICS OF BRANCHING: SPACELIKE PARTON K(N+1,1)=20000 K(N+1,2)=ISIGN(500+IABS(IFLB),2*IFLB+1) P(N+1,1)=P(IT,1) P(N+1,2)=0. P(N+1,3)=P(IT,3)+P(IS(JT),3) P(N+1,4)=P(IT,4)+P(IS(JT),4) P(N+1,5)=-SQRT(SNGL(DQ2(3))) K(N+2,1)=70001+N K(N+2,2)=K(IS(JT)+1,2) DO 250 J=1,5 250 P(N+2,J)=0. C...DEFINE COLOUR FLOW OF BRANCHING K(IS(JT),1)=20001+N ID1=IT IF((K(N+1,2).GE.501.AND.K(ID1,2).GE.501).OR.(K(N+1,2).LT.0.AND. & K(ID1,2).EQ.500).OR.(K(N+1,2).EQ.500.AND.K(ID1,2).EQ.500.AND. & RLU(0).GT.0.5).OR.(K(N+1,2).EQ.500.AND.K(ID1,2).LT.0)) & ID1=IS(JT) ID2=IT+IS(JT)-ID1 P(N+2,3)=ID1 P(N+2,4)=ID2 P(ID1+1,1)=N+1 P(ID1+1,2)=ID2 P(ID2+1,1)=ID1 P(ID2+1,2)=N+1 N=N+2 C...BOOST TO NEW CM-FRAME MST(1)=NS+1 IF(ILEP.EQ.0) THEN CALL LUROBO(0.,0.,-(P(N-1,1)+P(IS(JR),1))/(P(N-1,4)+P(IS(JR), & 4)),0.,-(P(N-1,3)+P(IS(JR),3))/(P(N-1,4)+P(IS(JR),4))) IR=N-1+(JT-1)*(IS(1)-N+1) CALL LUROBO(-ULANGL(P(IR,3),P(IR,1)),PAR(72)*RLU(0),0.,0.,0.) ELSE C...REORIENTATE EVENT WITHOUT CHANGING THE BOSON FOUR MOMENTUM DO 260 J=1,4 260 DPQ(J)=P(NQ,J) DBE1(4)=DPQ(4)+DBLE(P(N-1,4)) DO 270 J=1,3,2 270 DBE1(J)=-(DPQ(J)+DBLE(P(N-1,J)))/DBE1(4) DBE1(4)=1D0/DSQRT(1D0-DBE1(1)**2-DBE1(3)**2) DBEP=DBE1(1)*DPQ(1)+DBE1(3)*DPQ(3) DGABEP=DBE1(4)*(DBE1(4)*DBEP/(1D0+DBE1(4))+DPQ(4)) DO 280 J=1,3,2 280 DPQ(J)=DPQ(J)+DGABEP*DBE1(J) DPQ(4)=DBE1(4)*(DPQ(4)+DBEP) DPC(1)=DSQRT(DPQ(1)**2+DPQ(3)**2) DBE2(4)=-(DPQ(4)*DPC(1)-DPQS(2)*DSQRT(DPQS(2)**2+DPC(1)**2- & DPQ(4)**2))/(DPC(1)**2+DPQS(2)**2) THE=ULANGL(SNGL(DPQ(3)),SNGL(DPQ(1))) DBE2(1)=DBE2(4)*DSIN(DBLE(THE)) DBE2(3)=DBE2(4)*DCOS(DBLE(THE)) DBE2(4)=1D0/(1D0-DBE2(1)**2-DBE2(3)*2) C...CONSTRUCT THE COMBINED BOOST DPB(1)=DBE1(4)**2*DBE2(4)/(1D0+DBE1(4)) DPB(2)=DBE1(1)*DBE2(1)+DBE1(3)*DBE2(3) DPB(3)=DBE1(4)*DBE2(4)*(1D0+DPB(2)) DO 290 JB=1,3,2 290 ROBO(JB+2)=(DBE1(4)*DBE2(4)*DBE1(JB)+DBE2(4)*DBE2(JB)+ & DPB(1)*DBE1(JB)*DPB(2))/DPB(3) CALL LUROBO(0.,0.,ROBO(3),0.,ROBO(5)) IF(ILEP.EQ.1) THE=ULANGL(P(NS+1,3),P(NS+1,1)) IF(ILEP.EQ.2) THE=PAR(71)+ULANGL(P(NS+3,3),P(NS+3,1)) CALL LUROBO(-THE,PAR(72)*RLU(0),0.,0.,0.) ENDIF MST(1)=0 ENDIF C...SAVE QUANTITIES, LOOP BACK IS(JT)=N-1 IF(ILEP.EQ.2.AND.N.EQ.NS+4) IS(JT)=N-3 Q2S(JT)=Q2B DQ2(JT)=Q2B IF(IPY(14).GE.3.AND.IPY(14).LE.6) DQ2(JT)=Q2B/(1.-Z) DSH=DSHZ IF(Q2B.GE.0.5*PYPAR(22)) THEN IFLS(JT+2)=IFLS(JT) IFLS(JT)=IFLA XS(JT)=XA ZS(JT)=Z DO 300 IFL=-6,6 300 XFS(JT,IFL)=XFA(IFL) TEVS(JT)=TEVB ELSE IF(JT.EQ.1) IPU1=N-1 IF(JT.EQ.2) IPU2=N-1 ENDIF IF(MAX(IABS(1-ILEP)*Q2S(1),MIN(1,2-ILEP)*Q2S(2)).GE.0.5*PYPAR(22) &.OR.N.LE.NS+2) GOTO 130 IF(ILEP.EQ.0) THEN C...BOOST HARD SCATTERING PARTONS TO FRAME OF SHOWER INITIATORS DO 310 J=1,3 310 ROBO(J+2)=(P(NS+1,J)+P(NS+3,J))/(P(NS+1,4)+P(NS+3,4)) DO 320 J=1,5 320 P(N+2,J)=P(NS+1,J) MST(1)=N+2 MST(2)=N+2 CALL LUROBO(0.,0.,-ROBO(3),-ROBO(4),-ROBO(5)) ROBO(2)=ULANGL(P(N+2,1),P(N+2,2)) ROBO(1)=ULANGL(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2)) MST(1)=4 MST(2)=NS CALL LUROBO(ROBO(1),ROBO(2),ROBO(3),ROBO(4),ROBO(5)) MST(1)=0 MST(2)=0 ENDIF C...STORE USER INFORMATION K(21,1)=20001+NS K(23,1)=20003+NS DO 330 JT=1,2 KFL(1,JT)=ISIGN(500+IABS(IFLS(JT)),2*IFLS(JT)+1) 330 PYVAR(30+JT)=XS(JT) RETURN END C*********************************************************************** SUBROUTINE PYREMN(IPU1,IPU2) C...ADDS ON TARGET REMNANTS (ONE OR TWO FROM EACH SIDE) AND C...INCLUDES PRIMORDIAL KT. COMMON/LUJETS/N,K(2000,2),P(2000,5) COMMON/LUDAT1/MST(40),PAR(80) COMMON/PYPARA/IPY(80),PYPAR(80),PYVAR(80) COMMON/PYPROC/ISUB,KFL(3,2),X(2),SH,TH,UH,Q2,XSEC(0:40) DIMENSION KFLCH(2),KFLSP(2),CHI(2),PMS(6),IS(2),ROBO(5) C...FIND EVENT TYPE, SET POINTERS IF(IPU1.EQ.0.AND.IPU2.EQ.0) RETURN ILEP=0 IF(IPU1.EQ.0) ILEP=1 IF(IPU2.EQ.0) ILEP=2 IF(ISUB.EQ.7) ILEP=-1 IF(ILEP.EQ.1) IQ=21 IF(ILEP.EQ.2) IQ=23 IP=MAX(IPU1,IPU2) NS=N C...DEFINE INITIAL PARTONS, INCLUDING PRIMORDIAL KT 100 DO 120 I=3,4 IF(I.EQ.3) IPU=IPU1 IF(I.EQ.4) IPU=IPU2 K(I,1)=40000+I-2 DO 110 J=1,5 110 P(I,J)=0. IF(ISUB.EQ.7) THEN K(I,2)=500 SHS=0. ELSEIF(IPU.NE.0) THEN K(I,2)=K(IPU,2) P(I,5)=P(IPU,5) CALL PYPRKT(P(I,1),P(I,2)) PMS(I-2)=P(I,5)**2+P(I,1)**2+P(I,2)**2 ELSE K(I,2)=K(IQ,2) P(I,5)=-SQRT(Q2) PMS(I-2)=-Q2 SHS=(1.-X(5-I))*Q2/X(5-I)+PYVAR(7-I)**2 ENDIF 120 CONTINUE C...KINEMATICS CONSTRUCTION FOR INITIAL PARTONS IF(ILEP.EQ.0) SHS=PYVAR(31)*PYVAR(32)*PYVAR(2)+ &(P(3,1)+P(4,1))**2+(P(3,2)+P(4,2))**2 SHR=SQRT(MAX(0.,SHS)) IF(ILEP.EQ.0) THEN IF((SHS-PMS(1)-PMS(2))**2-4.*PMS(1)*PMS(2).LE.0.) GOTO 100 P(3,4)=0.5*(SHR+(PMS(1)-PMS(2))/SHR) P(3,3)=SQRT(MAX(0.,P(3,4)**2-PMS(1))) P(4,4)=SHR-P(3,4) P(4,3)=-P(3,3) ELSEIF(ILEP.EQ.1) THEN P(3,4)=P(IQ,4) P(3,3)=P(IQ,3) P(4,4)=P(IP,4) P(4,3)=P(IP,3) ELSEIF(ILEP.EQ.2) THEN P(3,4)=P(IP,4) P(3,3)=P(IP,3) P(4,4)=P(IQ,4) P(4,3)=P(IQ,3) ENDIF C...TRANSFORM PARTONS TO OVERALL CM-FRAME (NOT FOR LEPTOPRODUCTION) IF(ILEP.EQ.0) THEN MST(1)=3 MST(2)=4 ROBO(3)=(P(3,1)+P(4,1))/SHR ROBO(4)=(P(3,2)+P(4,2))/SHR CALL LUROBO(0.,0.,-ROBO(3),-ROBO(4),0.) ROBO(2)=ULANGL(P(3,1),P(3,2)) CALL LUROBO(0.,-ROBO(2),0.,0.,0.) ROBO(1)=ULANGL(P(3,3),P(3,1)) CALL LUROBO(-ROBO(1),0.,0.,0.,0.) MST(2)=MAX(IPY(47),IPU1,IPU2) CALL LUROBO(ROBO(1),ROBO(2),ROBO(3),ROBO(4),0.) ROBO(5)=MAX(-0.999999,MIN(0.999999,(PYVAR(31)-PYVAR(32))/ & (PYVAR(31)+PYVAR(32)))) CALL LUROBO(0.,0.,0.,0.,ROBO(5)) MST(1)=0 MST(2)=0 ENDIF C...CHECK INVARIANT MASS OF REMNANT SYSTEM: C...HADRONIC EVENTS OR LEPTOPRODUCTION IF(ILEP.LE.0) THEN IF(IPY(12).LE.0.OR.ISUB.EQ.7) PYVAR(33)=0. IF(IPY(12).LE.0.OR.ISUB.EQ.7) PYVAR(34)=0. PEH=P(3,4)+P(4,4)+0.5*PYVAR(1)*(PYVAR(33)+PYVAR(34)) PZH=P(3,3)+P(4,3)+0.5*PYVAR(1)*(PYVAR(33)-PYVAR(34)) SHH=(PYVAR(1)-PEH)**2-(P(3,1)+P(4,1))**2-(P(3,2)+P(4,2))**2- & PZH**2 PMMIN=P(1,5)+P(2,5)+ULMASS(0,K(3,2))+ULMASS(0,K(4,2)) IF(SHR.GE.PYVAR(1).OR.SHH.LE.(PMMIN+PYPAR(12))**2) THEN IPY(48)=1 RETURN ENDIF SHR=SQRT(SHH+(P(3,1)+P(4,1))**2+(P(3,2)+P(4,2))**2) ELSE PEI=P(IQ,4)+P(IP,4) PZI=P(IQ,3)+P(IP,3) PMS(ILEP)=MAX(0.,PEI**2-PZI**2+P(5-ILEP,1)**2+P(5-ILEP,2)**2) PMMIN=P(3-ILEP,5)+ULMASS(0,K(5-ILEP,2))+SQRT(PMS(ILEP)) IF(SHR.LE.PMMIN+PYPAR(12)) THEN IPY(48)=1 RETURN ENDIF ENDIF C...SUBDIVIDE REMNANT IF NECESSARY, STORE FIRST PARTON 130 I=NS-1 DO 160 JT=1,2 IF(JT.EQ.ILEP) GOTO 160 IF(JT.EQ.1) IPU=IPU1 IF(JT.EQ.2) IPU=IPU2 CALL PYSPLI(IPY(40+JT),KFL(1,JT),KFLCH(JT),KFLSP(JT)) I=I+2 IS(JT)=I K(I,1)=0 K(I,2)=KFLSP(JT) P(I,5)=ULMASS(0,K(I,2)) C...FIRST PARTON COLOUR CONNECTIONS AND TRANSVERSE MASS K(I+1,1)=70000+I K(I+1,2)=1000 IF(IPY(34).GE.1) K(I+1,2)=1000+JT DO 140 J=1,5 140 P(I+1,J)=0. IFLS=(3-ISIGN(1,KFLSP(JT)*(510-IABS(KFLSP(JT)))))/2 P(I+1,IFLS+2)=IPU P(IPU+1,3-IFLS)=I IF(KFLCH(JT).EQ.0) THEN P(I,1)=-P(JT+2,1) P(I,2)=-P(JT+2,2) PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2 ELSE C...WHEN EXTRA REMNANT PARTON OR HADRON: FIND RELATIVE PT, STORE CALL LUPTDI(1,P(I,1),P(I,2)) PMS(JT+2)=P(I,5)**2+P(I,1)**2+P(I,2)**2 I=I+2 K(I,1)=0 K(I,2)=KFLCH(JT) P(I,5)=ULMASS(0,K(I,2)) P(I,1)=-P(JT+2,1)-P(I-2,1) P(I,2)=-P(JT+2,2)-P(I-2,2) PMS(JT+4)=P(I,5)**2+P(I,1)**2+P(I,2)**2 C...RELATIVE DISTRIBUTION OF ENERGY; EXTRA PARTON COLOUR CONNECTION CALL PYCHID(IPY(40+JT),KFLCH(JT),CHI(JT)) PMS(JT)=PMS(JT+4)/CHI(JT)+PMS(JT+2)/(1.-CHI(JT)) K(I+1,1)=70000+I K(I+1,2)=1000 IF(IPY(34).GE.1) K(I+1,2)=1000+JT DO 150 J=1,5 150 P(I+1,J)=0. IF(IABS(KFLCH(JT)).GE.500) THEN IFLS=(3-ISIGN(1,KFLCH(JT)*(510-IABS(KFLCH(JT)))))/2 P(I+1,IFLS+2)=IPU P(IPU+1,3-IFLS)=I ELSE IF(IPY(34).GE.1) K(I,1)=JT ENDIF ENDIF 160 CONTINUE IF(SHR.LE.SQRT(PMS(1))+SQRT(PMS(2))) GOTO 130 N=I+1 C...RECONSTRUCT KINEMATICS OF REMNANTS DO 170 JT=1,2 IF(JT.EQ.ILEP) GOTO 170 PE=0.5*(SHR+(PMS(JT)-PMS(3-JT))/SHR) PZ=SQRT(PE**2-PMS(JT)) IF(KFLCH(JT).EQ.0) THEN P(IS(JT),4)=PE P(IS(JT),3)=PZ*(-1)**(JT-1) ELSE PW1=CHI(JT)*(PE+PZ) P(IS(JT)+2,4)=0.5*(PW1+PMS(JT+4)/PW1) P(IS(JT)+2,3)=0.5*(PW1-PMS(JT+4)/PW1)*(-1)**(JT-1) P(IS(JT),4)=PE-P(IS(JT)+2,4) P(IS(JT),3)=PZ*(-1)**(JT-1)-P(IS(JT)+2,3) ENDIF 170 CONTINUE C...HADRONIC EVENTS: BOOST REMNANTS TO CORRECT LONGITUDINAL FRAME IF(ILEP.LE.0) THEN MST(1)=NS+1 CALL LUROBO(0.,0.,0.,0.,-PZH/(PYVAR(1)-PEH)) MST(1)=0 C...LEPTOPRODUCTION EVENTS: BOOST COLLIDING SUBSYSTEM ELSE MST(1)=21 MST(2)=MAX(IP,IPY(47)) PEF=SHR-PE PZF=PZ*(-1)**(ILEP-1) PT2=P(5-ILEP,1)**2+P(5-ILEP,2)**2 PHIPT=ULANGL(P(5-ILEP,1),P(5-ILEP,2)) CALL LUROBO(0.,-PHIPT,0.,0.,0.) RQP=P(IQ,3)*(PT2+PEI**2)-P(IQ,4)*PEI*PZI SINTH=P(IQ,4)*SQRT(PT2*(PT2+PEI**2)/(RQP**2+PT2* & P(IQ,4)**2*PZI**2))*SIGN(1.,-RQP) CALL LUROBO(ASIN(SINTH),0.,0.,0.,0.) BETAX=(-PEI*PZI*SINTH+SQRT(PT2*(PT2+PEI**2-(PZI*SINTH)**2)))/ & (PT2+PEI**2) CALL LUROBO(0.,0.,BETAX,0.,0.) CALL LUROBO(0.,PHIPT,0.,0.,0.) PEM=P(IQ,4)+P(IP,4) PZM=P(IQ,3)+P(IP,3) BETAZ=(-PEM*PZM+PZF*SQRT(PZF**2+PEM**2-PZM**2))/(PZF**2+PEM**2) CALL LUROBO(0.,0.,0.,0.,BETAZ) MST(1)=3 MST(2)=4 CALL LUROBO(ASIN(SINTH),0.,BETAX,0.,0.) CALL LUROBO(0.,PHIPT,0.,0.,BETAZ) MST(1)=0 MST(2)=0 ENDIF RETURN END C*********************************************************************** SUBROUTINE PYCHID(KPART,KFL,CHI) C...ENERGY DISTRIBUTION AMONG FRAGMENTS OF HADRON REMNANT. COMMON/PYPARA/IPY(80),PYPAR(80),PYVAR(80) C...ENERGY DISTRIBUTION FOR PARTICLE INTO TWO JETS IF(IABS(KFL).GE.500) THEN IF(IABS(KPART).LE.40) CHIK=PYPAR(13) IF(IABS(KPART).GT.40) CHIK=PYPAR(15) IF(IPY(17).LE.1) THEN IF(IABS(KPART).LE.40) CHI=RLU(0) IF(IABS(KPART).GT.40) CHI=1.-SQRT(RLU(0)) ELSEIF(IPY(17).EQ.2) THEN CHI=1.-RLU(0)**(1./(1.+CHIK)) ELSEIF(IPY(17).EQ.3) THEN CUT=2.*0.3/PYVAR(1) 100 CHI=RLU(0)**2 IF((CHI**2/(CHI**2+CUT**2))**0.25*(1.-CHI)**CHIK.LT. & RLU(0)) GOTO 100 ELSE CUT=2.*0.3/PYVAR(1) CUTR=(1.+SQRT(1.+CUT**2))/CUT 110 CHIR=CUT*CUTR**RLU(0) CHI=(CHIR**2-CUT**2)/(2.*CHIR) IF((1.-CHI)**CHIK.LT.RLU(0)) GOTO 110 ENDIF C...ENERGY DISTRIBUTION FOR PARTICLE INTO JET PLUS PARTICLE ELSE IF(IPY(17).LE.1) THEN IF(IABS(KPART).LE.40) CHI=RLU(0) IF(IABS(KPART).GT.40) CHI=1.-SQRT(RLU(0)) ELSE IF(IABS(KPART).LE.40) CHI=1.-RLU(0)**(1./(1.+PYPAR(14))) IF(IABS(KPART).GT.40) CHI=1.-RLU(0)**(1./(1.+PYPAR(16))) ENDIF IF(IABS(KFL).GT.40) CHI=1.-CHI ENDIF RETURN END C*********************************************************************** SUBROUTINE PYPRKT(PTX,PTY) C...GIVES PRIMORDIAL KT TO REACTING PARTONS. COMMON/PYPARA/IPY(80),PYPAR(80),PYVAR(80) COMMON/LUDAT1/MST(40),PAR(80) C...NO PRIMORDIAL KT OR CHOSEN ACCORDING TO TRUNCATED GAUSSIAN OR C...EXPONENTIAL 100 IF(IPY(16).LE.0) THEN PT=0. ELSEIF(IPY(16).EQ.1) THEN PT=PYPAR(7)*SQRT(-ALOG(RLU(0))) ELSE RPT1=RLU(0) RPT2=RLU(0) PT=-PYPAR(8)*ALOG(RPT1*RPT2) ENDIF IF(PT.GT.PYPAR(11)) GOTO 100 PHI=PAR(72)*RLU(0) PTX=PT*COS(PHI) PTY=PT*SIN(PHI) RETURN END C*********************************************************************** SUBROUTINE PYSPLI(KPART,KFLIN,KFLCH,KFLSP) C...IN CASE OF A HADRON REMNANT WHICH IS MORE COMPLICATED THAN JUST A C...QUARK OR A DIQUARK, SPLIT IT INTO TWO (PARTONS OR HADRON + PARTON). IFLIN=KFLIN-ISIGN(500,KFLIN) KSIGN=ISIGN(1,KPART) IFL=KFLIN*KSIGN KFLCH=0 IDUM=0 IF(IABS(KPART).EQ.17) THEN C...DECOMPOSE PI+ (PI-). IF(IFL.EQ.501) THEN C...VALENCE U (UBAR) REMOVED. KFLSP=-502*KSIGN ELSEIF(IFL.EQ.-502) THEN C...VALENCE D (DBAR) REMOVED. KFLSP=501*KSIGN ELSEIF(KFLIN.EQ.500) THEN C...GLUON REMOVED. R=2.*RLU(0) IF(R.LT.1.) THEN KFLCH=501*KSIGN KFLSP=-502*KSIGN ELSE KFLCH=-502*KSIGN KFLSP=501*KSIGN ENDIF ELSEIF(IFL.GT.500.AND.IFL.NE.501) THEN C...SEA QUARK (ANTIQUARK) REMOVED. CALL LUIFLD(-IFLIN,0,KSIGN,IDUM,KFLCH) KFLSP=-502*KSIGN ELSEIF(IFL.LT.-500.AND.IFL.NE.-502) THEN C...SEA ANTIQUARK (QUARK) REMOVED. CALL LUIFLD(-IFLIN,0,-2*KSIGN,IDUM,KFLCH) KFLSP=501*KSIGN ENDIF ELSEIF(IABS(KPART).EQ.41) THEN C...DECOMPOSE PROTON (ANTIPROTON). IF(IFL.EQ.501) THEN C...VALENCE U (UBAR) REMOVED. R=4.*RLU(0) IF(R.LT.3.) THEN KFLSP=512*KSIGN ELSE KFLSP=521*KSIGN ENDIF ELSEIF(IFL.EQ.502) THEN C...VALENCE D (DBAR) REMOVED. KFLSP=511*KSIGN ELSEIF(KFLIN.EQ.500) THEN C...GLUON REMOVED. R=6.*RLU(0) IF(R.LT.3.) THEN KFLCH=501*KSIGN KFLSP=512*KSIGN ELSEIF(R.LT.4.) THEN KFLCH=501*KSIGN KFLSP=521*KSIGN ELSE KFLCH=502*KSIGN KFLSP=511*KSIGN ENDIF ELSEIF(IFL.GT.502) THEN C...SEA QUARK (ANTIQUARK) REMOVED. R=6*RLU(0) IF(R.LT.3.) THEN CALL LUIFLD(-IFLIN,0,KSIGN,IDUM,KFLCH) KFLSP=512*KSIGN ELSEIF(R.LT.4.) THEN CALL LUIFLD(-IFLIN,0,KSIGN,IDUM,KFLCH) KFLSP=521*KSIGN ELSE CALL LUIFLD(-IFLIN,0,2*KSIGN,IDUM,KFLCH) KFLSP=511*KSIGN ENDIF ELSEIF(IFL.LT.-500) THEN C...SEA ANTIQUARK (QUARK) REMOVED. 100 R=6*RLU(0) IF(R.LT.3.) THEN CALL LUIFLD(12*KSIGN,KSIGN,-IFLIN,IDUM,KFLCH) KFLSP=501*KSIGN ELSEIF(R.LT.4.) THEN CALL LUIFLD(21*KSIGN,2*KSIGN,-IFLIN,IDUM,KFLCH) KFLSP=501*KSIGN ELSE CALL LUIFLD(11*KSIGN,KSIGN,-IFLIN,IDUM,KFLCH) KFLSP=502*KSIGN ENDIF IF(KFLCH.EQ.0) GOTO 100 ENDIF ELSEIF(IABS(KPART).EQ.42) THEN C...DECOMPOSE NEUTRON (ANTINEUTRON). IF(IFL.EQ.501) THEN C...VALENCE U (UBAR) REMOVED. KFLSP=522*KSIGN ELSEIF(IFL.EQ.502) THEN C...VALENCE D (DBAR) REMOVED. R=4.*RLU(0) IF(R.LT.3.) THEN KFLSP=512*KSIGN ELSE KFLSP=521*KSIGN ENDIF ELSEIF(KFLIN.EQ.500) THEN C...GLUON REMOVED. R=6.*RLU(0) IF(R.LT.2.) THEN KFLCH=501*KSIGN KFLSP=522*KSIGN ELSEIF(R.LT.5.) THEN KFLCH=502*KSIGN KFLSP=512*KSIGN ELSE KFLCH=502*KSIGN KFLSP=521*KSIGN ENDIF ELSEIF(IFL.GT.502) THEN C...SEA QUARK (ANTIQUARK) REMOVED. R=6*RLU(0) IF(R.LT.2.) THEN CALL LUIFLD(-IFLIN,0,KSIGN,IDUM,KFLCH) KFLSP=522*KSIGN ELSEIF(R.LT.5.) THEN CALL LUIFLD(-IFLIN,0,2*KSIGN,IDUM,KFLCH) KFLSP=512*KSIGN ELSE CALL LUIFLD(-IFLIN,0,2*KSIGN,IDUM,KFLCH) KFLSP=521*KSIGN ENDIF ELSEIF(IFL.LT.-500) THEN C...SEA ANTIQUARK (QUARK) REMOVED. 110 R=6*RLU(0) IF(R.LT.2.) THEN CALL LUIFLD(22*KSIGN,KSIGN,-IFLIN,IDUM,KFLCH) KFLSP=501*KSIGN ELSEIF(R.LT.5.) THEN CALL LUIFLD(12*KSIGN,KSIGN,-IFLIN,IDUM,KFLCH) KFLSP=502*KSIGN ELSE CALL LUIFLD(21*KSIGN,2*KSIGN,-IFLIN,IDUM,KFLCH) KFLSP=502*KSIGN ENDIF IF(KFLCH.EQ.0) GOTO 110 ENDIF ENDIF RETURN END