C C MAIN program for HERA conditions; C needs to link files: C + js_sub.f < The essential Model C + vegas.f < Integrating routine C + pdf_grv.f < Structure functions C C###################################################################C C C C A program to simulate the production of J/psi mesons C C C C in the framework of Semihard Approach and/or Parton Model C C via gluon-gluon, photon-gluon or photon-photon fusion C C C C###################################################################C IMPLICIT DOUBLE PRECISION (A-G,P-Z) EXTERNAL FXN COMMON/FIGS/hsum,hfig(10) !weights accumulated in histograms COMMON/BVEG1/XL(10),XU(10),ACC COMMON/BVEGG/NDIM,NCALL,ITMX,NPRN COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC2,PSI COMMON/TYPE/ITYPE,MODEL COMMON/SEED/NUM COMMON/PAWC/HBOOK(1000000) CALL METRIC NUM=1 C C... Theoretical model, interaction type, and cms total energy MODEL = 1 ! 0 = Parton Model; ! 1 = Semihard Theory ITYPE = 2 ! 0 = proton*proton [glu+glu->Jpsi+gam] ! 1 = proton*proton [glu+glu->Jpsi+glu] ! 2 = lepton*proton [gam+glu->Jpsi+glu] ! 3 = lepton*lepton [gam+gam->Jpsi+gam] c CME = 7. ! sqrt(s) [GeV] HERMES c CME = 14. ! sqrt(s) [GeV] COMPASS c CME = 19. ! sqrt(s) [GeV] SMC c CME = 190. ! sqrt(s) [GeV] HERA CME = 300. ! sqrt(s) [GeV] HERA main C... Parameters and constants PI = 3.1415926 ! pi constant ALP = 1./137. ! electromagnetic constant ALS = 0.2 ! strong coupling constant XJ = 3.1 ! J/psi mass [GeV] PSI = 0.038/XJ ! J/psi wave function c XJ =9.5 !Upsilon mass [GeV] c PSI=0.4/XJ !Upsilon wave function XC = XJ/2. ! charmed quark mass QC2 =(2./3.)**2 ! charmed quark charge, squared XJ2 = XJ*XJ XC2 = XC*XC C C...Open output file. OPEN(UNIT=2,file='jpsi.outpt') C...Book histograms. CALL HLIMIT(1000000) CALL HROPEN(1,'BB','jpsi.histos','N',1024,ISTAT) hpt=8. CALL HBOOK1(102,'Pt(Jpsi)',48,0.,hpt,0.) CALL HBOOK1(112,'Pt(Jpsi)',48,0.,hpt,0.) CALL HBOOK1(8102,'Pt(Jpsi)',48,0.,hpt,0.) CALL HBOOK1(8112,'Pt(Jpsi)',48,0.,hpt,0.) hy =4. CALL HBOOK1(202,'Y(Jpsi)',40,-hy,hy,0.) CALL HBOOK1(212,'Y(Jpsi)',40,-hy,hy,0.) CALL HBOOK1(8202,'Y(Jpsi)',40,-hy,hy,0.) CALL HBOOK1(8212,'Y(Jpsi)',40,-hy,hy,0.) c CALL HBOOK1(302,'Z(Jpsi)',40,0.,1.,0.) CALL HBOOK1(312,'Z(Jpsi)',40,0.,1.,0.) CALL HBOOK1(8302,'Z(Jpsi)',40,0.,1.,0.) CALL HBOOK1(8312,'Z(Jpsi)',40,0.,1.,0.) hsum =0.0 DO 1 i=1,10 1 hfig(i)=0.0 C C... Integration limits for VEGAS NDIM=7 XL(1) =alog( 0.002) ! min transverse momentum of 1st parton XU(1) =alog( 5.000) ! max... XL(2) =alog( 0.002) ! min transverse momentum of 2nd parton XU(2) =alog(10.000) ! max... XL(3) =alog( 0.002) ! min transverse momentum of Jpsi XU(3) =alog( 8.000) ! max... IF(MODEL.EQ.0) XU(2)=alog(0.0021) RAPLIM= dlog(CME/XJ)*0.8 XL(4) =-RAPLIM ! min rapidity of Jpsi XU(4) = RAPLIM ! max... XL(5) =-RAPLIM*5. ! min rapidity of accompanying quantum XU(5) = RAPLIM*5. ! max... XL(6) = 0.D0 ! 1st parton azimuthal angle, units of 2*pi XU(6) = 1.D0 XL(7) = 0.D0 ! 2nd parton azimuthal angle, units of 2*pi XU(7) = 1.D0 C... Monte-Carlo integration by VEGAS NCALL = 10000 ! number of points per iteration ITMX = 10 ! number of iterations NPRN = 0 ! print/noprint statistics CALL VEGAS(FXN,AVGI,SD,CHI2A) C NCALL = 600000 ! number of points per iteration ITMX = 10 ! number of iterations NPRN = 1 ! print statistics CALL VEGAS1(FXN,AVGI,SD,CHI2A) C C... Write output: cross sections and histograms CALL HISTDO CALL HROUT(0,ICYCLE,' ') CALL HREND('BB') WRITE(2,102) hsum,hfig 102 format(//E12.3,//5E12.3,//5E12.3) WRITE(2,*) 'Last random number=',NUM CLOSE(UNIT=2) STOP END SUBROUTINE WRIOUT(hfxn) C*******************************************************************C C Filling output histograms C C*******************************************************************C IMPLICIT DOUBLE PRECISION (A-G,P-Z) DIMENSION jcut(10) !indicators for kinematical cuts COMMON/FIGS/hsum,hfig(10) !weights accumulated in histograms COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC,PSI COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/KINEM/x1,x2,g1T,g2T, pjT,Yj,Zj COMMON/OUTPT/CXPOS(-1:1),CXNEG(-1:1) COMMON/TYPE/ITYPE,MODEL tiny=1.D-5 C... Pseudurapidity of Jpsi C> th1=datan2(pjT+tiny,pj(1)) C> eta1=-dlog(tan(th1/2.)) C C___HERA cuts____________________________ DO 10 j=1,10 10 jcut(j)=1 IF(Zj.LE.0.05 .OR. Zj.GE.0.95) jcut(9)=0 IF(Zj.LE.0.40 .OR. Zj.GE.0.80) jcut(8)=0 C C...Histograms for presentation IF(jcut(9).EQ.0) GOTO 99 h0 = CXPOS(-1)+CXPOS(0)+CXPOS(1)+CXNEG(-1)+CXNEG(0)+CXNEG(1) h1 = CXPOS(-1) +CXPOS(1) h2 = CXNEG(-1) +CXNEG(1) h3 = CXPOS(-1) +CXPOS(1)+CXNEG(-1) +CXNEG(1) h4 = CXPOS(0) h5 = CXNEG(0) h6 = CXPOS(0) +CXNEG(0) h7 = CXPOS(-1)+CXPOS(0)+CXPOS(1) h8 = CXNEG(-1)+CXNEG(0)+CXNEG(1) h0=h0*hfxn h1=h1*hfxn h2=h2*hfxn h3=h3*hfxn h4=h4*hfxn h5=h5*hfxn h6=h6*hfxn h7=h7*hfxn h8=h8*hfxn hfig(1)=hfig(1)+h1 hfig(2)=hfig(2)+h2 hfig(3)=hfig(3)+h3 hfig(4)=hfig(4)+h4 hfig(5)=hfig(5)+h5 hfig(6)=hfig(6)+h6 hfig(7)=hfig(7)+h7 hfig(8)=hfig(8)+h8 hsum =hsum +h0 hpjT=sngl(pjT) np = 6 ! 6 bins per 1 GeV CALL HF1(102,hpjT,(h6*np)) CALL HF1(112,hpjT,(h0*np)) if(jcut(8).ne.0) CALL HF1(8102,hpjT,(h6*np)) if(jcut(8).ne.0) CALL HF1(8112,hpjT,(h0*np)) hYj= sngl(Yj) ny = 5 ! 5 bins per rapidity unit CALL HF1(202,hYj,(h6*ny)) CALL HF1(212,hYj,(h0*ny)) if(jcut(8).ne.0) CALL HF1(8202,hYj,(h6*ny)) if(jcut(8).ne.0) CALL HF1(8212,hYj,(h0*ny)) hZj= sngl(Zj) nz = 40 ! 40 bins per unit Z CALL HF1(302,hZj,(h6*nz)) CALL HF1(312,hZj,(h0*nz)) if(jcut(8).ne.0) CALL HF1(8302,hZj,(h6*nz)) if(jcut(8).ne.0) CALL HF1(8312,hZj,(h0*nz)) 99 CONTINUE RETURN END C###################################################################C C C C Routines to simulate the production of J/psi mesons C C C C in the framework of Semihard Approach and/or Parton Model C C via gluon-gluon, photon-gluon or photon-photon fusion C C C C###################################################################C DOUBLE PRECISION FUNCTION FXN(X,WGT) C*******************************************************************C C The integrand expression for VEGAS + kinematics C C*******************************************************************C IMPLICIT DOUBLE PRECISION (A-G,P-Z) DIMENSION X(10), SYMSYM(-1:1),ASYASY(-1:1) COMMON/BVEG1/XL(10),XU(10),ACC COMMON/BVEGG/NDIM,NCALL,ITMX,NPRN COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC,PSI COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/KINEM/x1,x2,g1T,g2T, pjT,yj,zj COMMON/OUTPT/CXPOS(-1:1),CXNEG(-1:1) COMMON/TYPE/ITYPE,MODEL COMMON/SEED/NUM FXN=0.D0 C... Random numbers within VEGAS g1T =dexp(X(1))!Transverse momentum of 1st parton g2T =dexp(X(2))!Transverse momentum of 2nd parton pjT =dexp(X(3))!Transverse momentum of J/psi yj = X(4) !Rapidity of J/psi y3 = X(5) !Rapidity of accompanying Gluon or Photon ph1 = X(6)*2.*PI !Azimuthal angle of 1st Gluon or Photon ph2 = X(7)*2.*PI !Azimuthal angle of 2nd Gluon or Photon C... Random numbers other than VEGAS C> RJ = RANDOM(NUM) C> phj = (2.*RJ)*PI !Azimuthal angle of J/psi phj = 0.D0 ph1 = ph1-phj ph2 = ph2-phj C C... Particle and Parton momenta, P(i) = (Beam, P_t, P_t, Energy) p1(1) = CME/2. p1(2) = 0.D0 p1(3) = 0.D0 p1(4) = CME/2. p2(1) =-CME/2. p2(2) = 0.D0 p2(3) = 0.D0 p2(4) = CME/2. g1(2) = g1T*dcos(ph1) g1(3) = g1T*dsin(ph1) g2(2) = g2T*dcos(ph2) g2(3) = g2T*dsin(ph2) pj(2) = pjT*dcos(phj) pj(3) = pjT*dsin(phj) tmj= dsqrt(XJ2 +pjT*pjT) g3(2) = g1(2) +g2(2) -pj(2) g3(3) = g1(3) +g2(3) -pj(3) tm3= dsqrt(g3(2)**2 +g3(3)**2) pj(1) = tmj*dsinh(yj) pj(4) = tmj*dcosh(yj) g3(1) = tm3*dsinh(y3) g3(4) = tm3*dcosh(y3) x1 = (tmj*dexp(yj) + tm3*dexp(y3))/CME x2 = (tmj*dexp(-yj)+ tm3*dexp(-y3))/CME IF(x1.GE..95D0 .OR. x2.GE..95D0) RETURN IF(x1.LE.1.D-4 .OR. x2.LE.1.D-4) RETURN y1 = 15.D0 y2 =-15.D0 IF(dabs(y1).GE.20.) RETURN IF(dabs(y2).GE.20.) RETURN DO 3 k=1,7 !Iterative solution arg1= CME*(1.-x1)-g2T*dexp(y2) arg2= CME*(1.-x2)-g1T*dexp(-y1) y1 = dlog(arg1/g1T) y2 =-dlog(arg2/g2T) 3 CONTINUE IF(y1.LE. 2.) RETURN IF(y2.GE.-2.) RETURN g1(1) = CME/2. -g1T*dsinh(y1) g1(4) = CME/2. -g1T*dcosh(y1) g2(1) =-CME/2. -g2T*dsinh(y2) g2(4) = CME/2. -g2T*dcosh(y2) ptest1=pj(1)+g3(1)-g1(1)-g2(1) ptest4=pj(4)+g3(4)-g1(4)-g2(4) IF(dabs(ptest1).GE. 1.D-5) RETURN IF(dabs(ptest4).GE. 1.D-5) RETURN zj = DOT(pj,p2)/DOT(g1,p2) IF(zj .LE. 0.D0) RETURN IF(zj .GE. 0.99) RETURN C C... Phase Space factor and integration Jacobian G1W = DOT(g1,g1) G2W = DOT(g2,g2) STOT= CME*CME SXX = XJ2 +2.*DOT(pj,g3) SXg = G2W +2.*DOT(p1,g2) IF(AF2(SXX,G1W,G2W).LE.0.D0) RETURN T1 = G1W +XJ2 -2.*DOT(g1,pj) T2 = G2W -2.*DOT(g2,g3) IF(GF(SXX,T1,0.D0,G1W,G2W,XJ2).GE.0.D0) RETURN IF(GF(SXX,T2,XJ2,G2W,G1W,0.D0).GE.0.D0) RETURN SXX = x1*x2*STOT SXg = x2*STOT IF(ITYPE.LE.1) DENOMR= STOT*SXX/PI/PI/4. IF(ITYPE.EQ.2) DENOMR= STOT*SXg*G1W**2 *(1.-x1)/PI IF(ITYPE.EQ.3) DENOMR= 4.*(STOT*G1W*G2W)**2 *(1.-x1)*(1.-x2) FACTOR= 2.*pjT/DENOMR *(g1T*g2T*pjT) c XJACOB= tmj*tm3*dsinh(dabs(yj-y3))/SXX IF(MODEL.EQ.0 .AND. ITYPE.LE.1) FACTOR=FACTOR*XJACOB IF(MODEL.EQ.0 .AND. ITYPE.EQ.2) FACTOR=FACTOR*XJACOB IF(MODEL.EQ.1 .AND. ITYPE.NE.3) FACTOR=FACTOR/10. C C... Couplings and other constants; averaging over spins and colors XNF= 3. ! Number of flavours Q2 = SXX/4.! Momentum scale c> Q2 = tmj*tmj ALS= ALPHAS(Q2,XNF) IF(ITYPE.EQ.0) FACTOR = FACTOR*PSI *ALP *ALS**2 /4./64. IF(ITYPE.EQ.1) FACTOR = FACTOR*PSI *ALS**3 /4./64. IF(ITYPE.EQ.2) FACTOR = FACTOR*PSI *ALP**2 *ALS**2 /4./8. IF(ITYPE.EQ.3) FACTOR = FACTOR*PSI *ALP**5 /4. C... Gluon and Photon Distribution Functions w1T = g1T*g1T w2T = g2T*g2T DEL1=0.999 DEL2=0.999 IF(MODEL.NE.0) GOTO 11 !This is Parton Model IF(ITYPE.LE.1) PDF1 = GRVg(x1,Q2,0)/x1 /(XU(1)-XL(1))/g1T IF(ITYPE.GE.2) PDF1 = 1.D0 IF(ITYPE.LE.2) PDF2 = GRVg(x2,Q2,0)/x2 /(XU(2)-XL(2))/g2T IF(ITYPE.EQ.3) PDF2 = 1.D0 IF(ITYPE.LE.1) DEL1 = x1*(1.-x1)**3 /0.12 IF(ITYPE.LE.2) DEL2 = x2*(1.-x2)**3 /0.12 chooseIF(ITYPE.LE.1) DEL1 = DELTAg(x1,2)/PDF1 chooseIF(ITYPE.LE.2) DEL2 = DELTAg(x2,2)/PDF2 chooseIF(ITYPE.LE.1) DEL1 = GRSVg(x1,2) chooseIF(ITYPE.LE.2) DEL2 = GRSVg(x2,2) chooseIF(ITYPE.LE.1) DEL1 = GSg(x1,1)/PDF1 chooseIF(ITYPE.LE.2) DEL2 = GSg(x2,1)/PDF2 11 IF(MODEL.NE.1) GOTO 12 !This is Semihard Theory IF(ITYPE.LE.1) PDF1 = BLUEML(x1,w1T,Q2)/x1 *2.*g1T IF(ITYPE.GE.2) PDF1 = 1.D0 IF(ITYPE.LE.2) PDF2 = BLUEML(x2,w2T,Q2)/x2 *2.*g2T IF(ITYPE.EQ.3) PDF2 = 1.D0 12 FACTOR =FACTOR*PDF1*PDF2 IF(PDF1.LT.0.D0) then c write(6,*) 'Negative PDF1=',PDF1 FACTOR=0.D0 endif IF(PDF2.LT.0.D0) then c write(6,*) 'Negative PDF2=',PDF2 FACTOR=0.D0 endif C C... Partonic cross section CALL XSEC(SYMSYM,ASYASY) DO 13 i=-1,1 CXPOS(i) = FACTOR*(SYMSYM(i) + DEL1*DEL2*ASYASY(i))/2. CXNEG(i) = FACTOR*(SYMSYM(i) - DEL1*DEL2*ASYASY(i))/2. IF(CXPOS(i).lt.0.D0) then cw write(6,*) 'Negative CXPOS(',i,') !!!!!!!!!!' CXPOS(i)=0.D0 endif IF(CXNEG(i).lt.0.D0) then cw write(6,*) 'Negative CXNEG(',i,') !!!!!!!!!!' CXNEG(i)=0.D0 endif 13 CONTINUE FXN = FACTOR*(SYMSYM(-1)+SYMSYM(0)+SYMSYM(1)) IF(FXN.LT.0.D0) then write(6,*) 'Negative total cross section !!!!!!!!!' cw write(6,*) 'x1,x2',x1,x2 cw write(6,*) 'y1,y2',y1,y2 cw write(6,*) 'yj,y3',yj,y3 cw write(6,*) 'g1T,g2T,pjT',g1T,g2T,pjT FXN=0.D0 endif C C... Filling histograms IF(NPRN.NE.1) RETURN hwgt = sngl(WGT)/FLOAT(ITMX) CALL WRIOUT(hwgt) RETURN END SUBROUTINE XSEC(SYMSYM,ASYASY) C*******************************************************************C C Partonic differential cross section C C*******************************************************************C IMPLICIT DOUBLE PRECISION (A-G,P-Z) DIMENSION SYMSYM(-1:1),ASYASY(-1:1) COMMON/TYPE/ITYPE,MODEL COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC,PSI COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/LOOPJ/AMPJ(4,4,-1:1,4) COMMON/LOOP1/SYM1(4,4),ASY1(4,4) COMMON/LOOP2/SYM2(4,4),ASY2(4,4) COMMON/GMUNU/DF(4,4),DC(4) C...Dot-products, polarizations CALL SCALAR CALL GAUGEJ IF(MODEL.EQ.0) CALL GAUGE0 IF(MODEL.EQ.1) CALL GAUGE1 C...Evaluating the Matrix Elements CALL FEYNJ DO 1 j=-1,1 SYMSYM(j)=0.D0 1 ASYASY(j)=0.D0 IF(ITYPE.EQ.0) COLOR= 8./3. *QC ! 0 =glu*glu -> Jpsi*gam IF(ITYPE.EQ.1) COLOR=10./9. ! 1 =glu*glu -> Jpsi*glu IF(ITYPE.EQ.2) COLOR= 8./3. *QC ! 2 =gam*glu -> Jpsi*glu IF(ITYPE.EQ.3) COLOR= 12. *QC**3 ! 3 =gam*gam -> Jpsi*gam DO 11 n1=1,4 DO 11 n2=1,4 DO 11 m1=1,4 DO 11 m2=1,4 SYMM=SYM1(n1,m1)*SYM2(n2,m2)*DC(n1)*DC(n2)*DC(m1)*DC(m2) ASYM=ASY1(n1,m1)*ASY2(n2,m2)*DC(n1)*DC(n2)*DC(m1)*DC(m2) DO 11 j=-1,1 DO 11 l=1,4 SYMSYM(j)=SYMSYM(j)-SYMM*AMPJ(n1,n2,j,l)*AMPJ(m1,m2,j,l)*DC(l) ASYASY(j)=ASYASY(j)-ASYM*AMPJ(n1,n2,j,l)*AMPJ(m1,m2,j,l)*DC(l) 11 CONTINUE DO 2 j=-1,1 SYMSYM(j)=SYMSYM(j)*COLOR/25.E26*1.E33 !Conversion to nanobarns ASYASY(j)=ASYASY(j)*COLOR/25.E26*1.E33 !Conversion to nanobarns 2 CONTINUE RETURN END SUBROUTINE SCALAR C*******************************************************************C C Dot-products of particle momenta C C*******************************************************************C IMPLICIT DOUBLE PRECISION (A-G,P-Z) COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/DOTPR/DG1P,DG2P,DG3P,DG12,DG13,DG23,G1W,G2W DG1P=DOT(g1,pj) DG2P=DOT(g2,pj) DG3P=DOT(g3,pj) DG12=DOT(g1,g2) DG13=DOT(g1,g3) DG23=DOT(g2,g3) G1W=DOT(g1,g1) G2W=DOT(g2,g2) RETURN END SUBROUTINE GAUGEJ C********************************************************************C C Polarization vector of outgoing J/psi C C********************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC,PSI COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/SPINJ/EJ(4,-1:1) QJ=dsqrt(pj(3)**2 +pj(2)**2 +pj(1)**2) Qt=dsqrt(pj(3)**2 +pj(2)**2) CT=pj(1)/QJ ST=Qt/QJ CP=pj(2)/Qt SP=pj(3)/Qt C... Scalar polarization EJ(1,0)= pj(1)*pj(4)/QJ/XJ EJ(2,0)= pj(2)*pj(4)/QJ/XJ EJ(3,0)= pj(3)*pj(4)/QJ/XJ EJ(4,0)= QJ/XJ C... Transverse polarization EJ(1,-1)=-ST EJ(2,-1)= CT*CP EJ(3,-1)= CT*SP EJ(4,-1)= 0.D0 C... Transverse polarization EJ(1,1)= 0.D0 EJ(2,1)=-SP EJ(3,1)= CP EJ(4,1)= 0.D0 RETURN END SUBROUTINE GAUGE0 C*******************************************************************C C Boson polarization matrix within conventional parton model C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) DIMENSION ps1(4),ps2(4),aux1(4),aux2(4),auxM(4),auxN(4) COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/KINEM/x1,x2,g1T,g2T, pjT,yj,zj COMMON/DOTPR/DG1P,DG2P,DG3P,DG12,DG13,DG23,G1W,G2W COMMON/LOOP1/SYM1(4,4),ASY1(4,4) COMMON/LOOP2/SYM2(4,4),ASY2(4,4) COMMON/GMUNU/DF(4,4),DC(4) COMMON/TYPE/ITYPE,MODEL DO 1 i=1,4 aux1(i)=0.D0 aux2(i)=0.D0 ps1(i)=2.*p1(i)-g1(i) 1 ps2(i)=2.*p2(i)-g2(i) DO 2 i=2,3 aux1(i)=g1(i)/g1T 2 aux2(i)=g2(i)/g2T IF(ITYPE-2) 10,20,30 10 DO 19 M=1,4 !...Gluon(1) * Gluon(2) DO 19 N=1,4 SYM1(M,N)=0.D0 SYM2(M,N)=0.D0 ASY1(M,N)=0.D0 ASY2(M,N)=0.D0 19 CONTINUE SYM1(2,2)= 1.D0 SYM1(3,3)= 1.D0 SYM2(2,2)= 1.D0 SYM2(3,3)= 1.D0 ASY1(2,3)=-1.D0 ASY1(3,2)= 1.D0 ASY2(2,3)=-1.D0 ASY2(3,2)= 1.D0 GOTO 40 20 DO 29 M=1,4 !..Photon(1) * Gluon(2) DO 29 N=1,4 SYM1(M,N)= 2.*(ps1(M)*ps1(N) -g1(M)*g1(N) +DF(M,N)*G1W) c SYM1(M,N)= 8.*aux1(M)*aux1(N)*g1T*g1T/x1/x1 SYM2(M,N)= 0.D0 ASY2(M,N)= 0.D0 DO 21 L=1,4 auxM(L)=0.D0 21 auxN(L)=0.D0 auxM(M)=DC(M) auxN(N)=DC(N) ASY1(M,N)= 2.*EPS(ps1,g1,auxM,auxN) 29 CONTINUE SYM2(2,2)= 1.D0 SYM2(3,3)= 1.D0 ASY2(2,3)=-1.D0 ASY2(3,2)= 1.D0 GOTO 40 30 DO 39 M=1,4 !.Photon(1) * Photon(2) DO 39 N=1,4 SYM1(M,N)= 2.*(ps1(M)*ps1(N) -g1(M)*g1(N) +DF(M,N)*G1W) SYM2(M,N)= 2.*(ps2(M)*ps2(N) -g2(M)*g2(N) +DF(M,N)*G2W) c SYM1(M,N)= 8.*aux1(M)*aux1(N)*g1T*g1T/x1/x1 c SYM2(M,N)= 8.*aux2(M)*aux2(N)*g2T*g2T/x2/x2 DO 31 L=1,4 auxM(L)=0.D0 31 auxN(L)=0.D0 auxM(M)=DC(M) auxN(N)=DC(N) ASY1(M,N)= 2.*EPS(ps1,g1,auxM,auxN) ASY2(M,N)= 2.*EPS(ps2,g2,auxM,auxN) 39 CONTINUE 40 CONTINUE RETURN !Comment out this line for gauge invariance tests DO 99 M=1,4 DO 99 N=1,4 SYM1(M,N)= g1(M)*g1(N) !For gauge invariance tests only SYM2(M,N)= g2(M)*g2(N) !For gauge invariance tests only ASY1(M,N)= 0.D0 ASY2(M,N)= 0.D0 99 CONTINUE RETURN END SUBROUTINE GAUGE1 C*******************************************************************C C Boson polarization matrix within semihard approach C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) DIMENSION ps1(4),ps2(4),aux1(4),aux2(4),auxM(4),auxN(4) COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/KINEM/x1,x2,g1T,g2T, pjT,yj,zj COMMON/DOTPR/DG1P,DG2P,DG3P,DG12,DG13,DG23,G1W,G2W COMMON/LOOP1/SYM1(4,4),ASY1(4,4) COMMON/LOOP2/SYM2(4,4),ASY2(4,4) COMMON/GMUNU/DF(4,4),DC(4) COMMON/TYPE/ITYPE,MODEL DO 1 i=1,4 aux1(i)=0.D0 aux2(i)=0.D0 ps1(i)=2.*p1(i)-g1(i) 1 ps2(i)=2.*p2(i)-g2(i) DO 2 i=2,3 aux1(i)=g1(i)/g1T 2 aux2(i)=g2(i)/g2T xgt1=x1*x1/g1T/g1T/4. xgt2=x2*x2/g2T/g2T/4. IF(ITYPE-2) 10,20,30 10 DO 19 M=1,4 !...Gluon(1) * Gluon(2) DO 19 N=1,4 SYM1(M,N)= 2.*(ps1(M)*ps1(N) -g1(M)*g1(N) +DF(M,N)*G1W)*xgt1 SYM2(M,N)= 2.*(ps2(M)*ps2(N) -g2(M)*g2(N) +DF(M,N)*G2W)*xgt2 c SYM1(M,N)= 2.*aux1(M)*aux1(N) c SYM2(M,N)= 2.*aux2(M)*aux2(N) DO 11 L=1,4 auxM(L)=0.D0 11 auxN(L)=0.D0 auxM(M)=DC(M) auxN(N)=DC(N) ASY1(M,N)= 2.*EPS(ps1,g1,auxM,auxN)*xgt1 ASY2(M,N)= 2.*EPS(ps2,g2,auxM,auxN)*xgt2 19 CONTINUE GOTO 40 20 DO 29 M=1,4 !..Photon(1) * Gluon(2) DO 29 N=1,4 SYM1(M,N)= 2.*(ps1(M)*ps1(N) -g1(M)*g1(N) +DF(M,N)*G1W) c SYM1(M,N)= 2.*aux1(M)*aux1(N)/xgt1 SYM2(M,N)= 2.*(ps2(M)*ps2(N) -g2(M)*g2(N) +DF(M,N)*G2W)*xgt2 c SYM2(M,N)= 2.*aux2(M)*aux2(N) DO 21 L=1,4 auxM(L)=0.D0 21 auxN(L)=0.D0 auxM(M)=DC(M) auxN(N)=DC(N) ASY1(M,N)= 2.*EPS(ps1,g1,auxM,auxN) ASY2(M,N)= 2.*EPS(ps2,g2,auxM,auxN)*xgt2 29 CONTINUE GOTO 40 30 DO 39 M=1,4 !.Photon(1) * Photon(2) DO 39 N=1,4 SYM1(M,N)= 2.*(ps1(M)*ps1(N) -g1(M)*g1(N) +DF(M,N)*G1W) SYM2(M,N)= 2.*(ps2(M)*ps2(N) -g2(M)*g2(N) +DF(M,N)*G2W) c SYM1(M,N)= 2.*aux1(M)*aux1(N)/xgt1 c SYM2(M,N)= 2.*aux2(M)*aux2(N)/xgt2 DO 31 L=1,4 auxM(L)=0.D0 31 auxN(L)=0.D0 auxM(M)=DC(M) auxN(N)=DC(N) ASY1(M,N)= 2.*EPS(ps1,g1,auxM,auxN) ASY2(M,N)= 2.*EPS(ps2,g2,auxM,auxN) 39 CONTINUE 40 CONTINUE RETURN !Comment out this line for gauge invariance tests DO 99 M=1,4 DO 99 N=1,4 SYM1(M,N)= g1(M)*g1(N) !For gauge invariance tests only SYM2(M,N)= g2(M)*g2(N) !For gauge invariance tests only ASY1(M,N)= 0.D0 ASY2(M,N)= 0.D0 99 CONTINUE RETURN END SUBROUTINE FEYNJ C*******************************************************************C C Matrix Elements for all Feynman diagrams C C*******************************************************************C IMPLICIT DOUBLE PRECISION (A-G,P-Z) COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC,PSI COMMON/MOMEN/p1(4),p2(4),g1(4),g2(4),pj(4),g3(4) COMMON/DOTPR/DG1P,DG2P,DG3P,DG12,DG13,DG23,W1,W2 COMMON/LOOPJ/AMPJ(4,4,-1:1,4) COMMON/GMUNU/DF(4,4),DC(4) COMMON/SPINJ/EJ(4,-1:1) C...Denominators W3=0.D0 !real final gluon Z1=(W1-DG1P)*(W3+DG3P) Z2=(W1-DG1P)*(W2-DG2P) Z3=(W2-DG2P)*(W3+DG3P) C...Fill main array DO 1 j=-1,1 DO 1 n1=1,4 DO 1 n2=1,4 DO 1 n3=1,4 AMPJ(n1,n2,j,n3)=0.D0 1 CONTINUE if(dabs(Z1).LT.1.D-3) return if(dabs(Z2).LT.1.D-3) return if(dabs(Z3).LT.1.D-3) return DO 10 n1=1,4 DO 10 n2=1,4 DO 10 n3=1,4 DO 10 nj=1,4 C...Diagram =<< 1 >>= TERM1= - DF(n1,n2)*DF(n3,nj)*DG13 - DF(n1,n2)*g1(n3)*g1(nj) + + DF(n1,n2)*g1(n3)*g3(nj) - DF(n1,n2)*g1(nj)*g2(n3) + DF(n1,n3)* + DF(n2,nj)*DG13 - DF(n1,n3)*g1(n2)*g3(nj) - DF(n1,n3)*g1(nj)* + g3(n2) - DF(n1,nj)*DF(n2,n3)*DG13 + DF(n1,nj)*g1(n2)*g1(n3) + + DF(n1,nj)*g1(n2)*g2(n3) + DF(n1,nj)*g1(n3)*g3(n2) + DF(n2,n3)* + g1(nj)*g3(n1) + DF(n2,n3)*g2(n1)*g3(nj) - DF(n2,n3)*g3(n1)* + g3(nj) - DF(n2,nj)*g1(n3)*g2(n1) - DF(n2,nj)*g2(n1)*g2(n3) + + DF(n2,nj)*g2(n3)*g3(n1) + DF(n3,nj)*g1(n2)*g3(n1) - DF(n3,nj)* + g2(n1)*g3(n2) + DF(n3,nj)*g3(n1)*g3(n2) C...Diagram =<< 2 >>= TERM2= - DF(n1,n2)*DF(n3,nj)*DG12 + DF(n1,n2)*g1(n3)*g2(nj) + + DF(n1,n2)*g1(nj)*g2(n3) + DF(n1,n3)*DF(n2,nj)*DG12 - DF(n1,n3)* + g1(n2)*g1(nj) - DF(n1,n3)*g1(n2)*g2(nj) + DF(n1,n3)*g1(nj)* + g3(n2) + DF(n1,nj)*DF(n2,n3)*DG12 + DF(n1,nj)*g1(n2)*g1(n3) - + DF(n1,nj)*g1(n2)*g2(n3) - DF(n1,nj)*g1(n3)*g3(n2) - DF(n2,n3)* + g1(nj)*g2(n1) - DF(n2,n3)*g2(n1)*g2(nj) + DF(n2,n3)*g2(nj)* + g3(n1) - DF(n2,nj)*g1(n3)*g2(n1) + DF(n2,nj)*g2(n1)*g2(n3) - + DF(n2,nj)*g2(n3)*g3(n1) + DF(n3,nj)*g1(n2)*g3(n1) + DF(n3,nj)* + g2(n1)*g3(n2) - DF(n3,nj)*g3(n1)*g3(n2) C...Diagram =<< 3 >>= TERM3= - DF(n1,n2)*DF(n3,nj)*DG23 - DF(n1,n2)*g1(n3)*g2(nj) - + DF(n1,n2)*g2(n3)*g2(nj) + DF(n1,n2)*g2(n3)*g3(nj) - DF(n1,n3)* + DF(n2,nj)*DG23 + DF(n1,n3)*g1(n2)*g3(nj) + DF(n1,n3)*g2(nj)* + g3(n2) - DF(n1,n3)*g3(n2)*g3(nj) + DF(n1,nj)*DF(n2,n3)*DG23 - + DF(n1,nj)*g1(n2)*g1(n3) - DF(n1,nj)*g1(n2)*g2(n3) + DF(n1,nj)* + g1(n3)*g3(n2) - DF(n2,n3)*g2(n1)*g3(nj) - DF(n2,n3)*g2(nj)* + g3(n1) + DF(n2,nj)*g1(n3)*g2(n1) + DF(n2,nj)*g2(n1)*g2(n3) + + DF(n2,nj)*g2(n3)*g3(n1) - DF(n3,nj)*g1(n2)*g3(n1) + DF(n3,nj)* + g2(n1)*g3(n2) + DF(n3,nj)*g3(n1)*g3(n2) DO 10 j=-1,1 AMPJ(n1,n2,j,n3) = AMPJ(n1,n2,j,n3)+ + 4.*XC*(TERM1/Z1 +TERM2/Z2 +TERM3/Z3)*EJ(nj,j)*DC(nj) 10 CONTINUE RETURN END DOUBLE PRECISION FUNCTION ALPHAS(Q2,XNF) C*******************************************************************C C Strong coupling constant C C*******************************************************************C IMPLICIT DOUBLE PRECISION (A-G,P-Z) ALPHAS=0.D0 ALAM2 =(0.2)**2 Q2RUN =dabs(Q2) IF(Q2RUN.LE.ALAM2) RETURN ALPHA =12.*3.1415926/(33.-2.*XNF)/dlog(Q2RUN/ALAM2) IF(ALPHA.GE.1.D0) ALPHA=1.D0 ALPHAS=ALPHA RETURN END DOUBLE PRECISION FUNCTION AF2(X,Y,Z) IMPLICIT DOUBLE PRECISION (A-G,P-Z) AF2 = X*X + Y*Y + Z*Z - 2.*X*Y - 2.*X*Z - 2.*Y*Z RETURN END DOUBLE PRECISION FUNCTION AF(X,Y,Z) IMPLICIT DOUBLE PRECISION (A-G,P-Z) AF = dsqrt(X*X + Y*Y + Z*Z - 2.*X*Y - 2.*X*Z - 2.*Y*Z) RETURN END DOUBLE PRECISION FUNCTION GF(X,Y,Z,U,V,W) IMPLICIT DOUBLE PRECISION (A-G,P-Z) GF = X*Z*W + X*U*V + Y*Z*V + Y*U*W . -X*Y*(Z+U+V+W-X-Y) -Z*U*(X+Y+V+W-Z-U) -V*W*(X+Y+Z+U-V-W) RETURN END SUBROUTINE METRIC IMPLICIT DOUBLE PRECISION(A-G,P-Z) COMMON/GMUNU/DF(4,4),DC(4) DO 1 I=1,4 DO 1 J=1,4 DF(I,J)= 0.D0 1 CONTINUE !P(i) =(Beam, P_t, P_t, Energy) DO 2 I=1,3 DF(I,I)=-1.D0 2 DC(I) =-1.D0 DF(4,4)= 1.D0 DC(4) = 1.D0 RETURN END DOUBLE PRECISION FUNCTION DOT(xx,yy) IMPLICIT DOUBLE PRECISION (A-G,P-Z) DIMENSION xx(4),yy(4) !P(i) = (Beam, P_t, P_t, Energy) DOT= -xx(1)*yy(1) -xx(2)*yy(2) -xx(3)*yy(3) +xx(4)*yy(4) RETURN END DOUBLE PRECISION FUNCTION EPS(q1,q2,q3,q4) IMPLICIT DOUBLE PRECISION (A-G,P-Z) DIMENSION q1(4),q2(4),q3(4),q4(4) G1 = q1(2)*q2(3)*q3(4) +q1(3)*q2(4)*q3(2) +q1(4)*q2(2)*q3(3) . -q1(4)*q2(3)*q3(2) -q1(2)*q2(4)*q3(3) -q1(3)*q2(2)*q3(4) G2 = q1(1)*q2(3)*q3(4) +q1(3)*q2(4)*q3(1) +q1(4)*q2(1)*q3(3) . -q1(4)*q2(3)*q3(1) -q1(1)*q2(4)*q3(3) -q1(3)*q2(1)*q3(4) G3 = q1(1)*q2(2)*q3(4) +q1(2)*q2(4)*q3(1) +q1(4)*q2(1)*q3(2) . -q1(4)*q2(2)*q3(1) -q1(1)*q2(4)*q3(2) -q1(2)*q2(1)*q3(4) G4 = q1(1)*q2(2)*q3(3) +q1(2)*q2(3)*q3(1) +q1(3)*q2(1)*q3(2) . -q1(3)*q2(2)*q3(1) -q1(1)*q2(3)*q3(2) -q1(2)*q2(1)*q3(3) EPS= -G1*q4(1) +G2*q4(2) -G3*q4(3) +G4*q4(4) RETURN END DOUBLE PRECISION FUNCTION BLUEML(x,g2,Q2) C*******************************************************************C C Gluon distribution by the method of J.Bluemlein C C BLUEML = x*gluon(x,kt2,Q2) DESY 94-121 C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) EXTERNAL BJ0,BI0 COMMON/CONVOL/cx,cg2,cQ2 BLUEML=0.D0 cx = x cg2= g2 cQ2= dabs(Q2) acc= 2.D-3 !Precision of integration IF(cg2.GE.cQ2) BLUEML=DGAUSS(BI0, x, 1.D0, acc) IF(cg2.LT.cQ2) BLUEML=DGAUSS(BJ0, x, 1.D0, acc) BLUEML=BLUEML/cg2 RETURN END DOUBLE PRECISION FUNCTION BJ0(eta) IMPLICIT DOUBLE PRECISION(A-G,P-Z) COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC,PSI COMMON/CONVOL/cx,cg2,cQ2 LO=0 !Leading order GRV Alsb= ALS*3./PI argj= 2.*dsqrt(-Alsb*dlog(eta)*dabs(dlog(cg2/cQ2))) argg= cx/eta BJ0 = Alsb/eta * DBESJ0(argj) * GRVg(argg,cQ2,LO) RETURN END DOUBLE PRECISION FUNCTION BI0(eta) IMPLICIT DOUBLE PRECISION(A-G,P-Z) COMMON/CONST/CME,PI,ALP,ALS, XC,XC2,XJ,XJ2, QC,PSI COMMON/CONVOL/cx,cg2,cQ2 LO=0 !Leading order GRV Alsb= ALS*3./PI argj= 2.*dsqrt(-Alsb*dlog(eta)*dabs(dlog(cg2/cQ2))) argg= cx/eta BI0 = Alsb/eta * DBESI0(argj) * GRVg(argg,cQ2,LO) RETURN END C###################################################################C C Parton distributions C C M.Gluck & E.Reya & A.Vogt: DO-TH 94/24; DESY 94-206 C C###################################################################C C C DOUBLE PRECISION FUNCTION GRVg(X,Q2,L) C*******************************************************************C C Gluon distribution function: GRVg = x*gluon(x,Q2) C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) C... L=0 >Leading Order; L=1 >NLO (MS bar); L=2 >NLO (DIS) QQ2=DABS(Q2) IF(QQ2.LT.0.4) QQ2=0.4D0 T0=DLOG(.23D0/(0.232**2)) T =DLOG(QQ2/(0.232**2)) S =DLOG(T/T0) S2=S**2 S3=S2*S SR=DSQRT(S) IF(L-1) 1,2,3 C... Leading Order set 1 AL = 0.524 BE = 1.088 AA = 1.742 -0.930*S BB = -0.399*S2 A = 7.486 -2.185*S B = 16.69 -22.74*S +5.779*S2 C =-25.59 +29.71*S -7.296*S2 D = 2.792 +2.215*S +0.422*S2 -0.104*S3 E = 0.807 +2.005*S EP = 3.841 +0.316*S GOTO 4 C... Next-to-Leading Order set, MS bar scheme 2 AL = 1.014 BE = 1.738 AA = 1.724 +0.157*S BB = 0.800 +1.016*S A = 7.517 -2.547*S B = 34.09 +17.47*S -52.21*DSQRT(S) C = 4.039 +1.491*S D = 3.404 +0.830*S E =-1.112 +3.438*S -0.302*S2 EP = 3.256 +0.436*S GOTO 4 C... Next-to-Leading order set, DIS scheme 3 AL = 1.258 BE = 1.846 AA = 2.423 BB = 2.427 +1.311*S -0.153*S2 A = 25.09 -7.935*S B =-14.84 +72.18*S -124.3*SR C = 590.3 -173.8*S D = 5.196 +1.857*S E =-1.648 +3.988*S -0.432*S2 EP = 3.232 -0.542*S 4 GRVg =(X**AA *(A +B*X +C*X**2)*(-DLOG(X))**BB . +S**AL *DEXP(-E +DSQRT(-EP*S**BE*DLOG(X))))*(1.-X)**D RETURN END DOUBLE PRECISION FUNCTION EPA(X,Q2) C*******************************************************************C C Equivalent Photon Approximation (Weiszaecker & Williams) C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) EPA =0.D0 IF(dabs(Q2).LE.1.D-6 .OR. X.LE.0.D0 .OR. X.GE.1.D0) RETURN EPA = (1.+(1.-X)**2)/DABS(Q2) RETURN END DOUBLE PRECISION FUNCTION WWA(X) C*******************************************************************C C Weiszaecker & Williams Approximation C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) WWA =0.D0 IF(X.LE.0.D0 .OR. X.GE.1.D0) RETURN WWA = (1.+(1.-X)**2) RETURN END C###################################################################C C Polarized Parton distributions C C###################################################################C C DOUBLE PRECISION FUNCTION DELTAg(X,L) C*******************************************************************C C Gluon distribution: DELTAg = x*(gluon+(x) - gluon-(x)) C C PL B378, 255 (1996) C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) IF(L-2) 1,2,3 C..."AB" set 1 AA= 0.1 A = -0.47 B = 2.6 E = 1.52/.876 GOTO 4 C..."AR" set 2 AA= 0.9 A = 0.03 B = 7.0 E = 1.11/.127 GOTO 4 C..."OS" set 3 AA= -1.25 A = 0.35 B = 4.0 E = 0.99/.712E-1 4 DELTAg = X**A *(1.-X)**B *(1.+AA*X) *E*X RETURN END DOUBLE PRECISION FUNCTION GRSVg(X,L) C*******************************************************************C C Gluon asymmetry: GRSVg = (gluon+(x) - gluon-(x)) / gluon(x) C C M.Gluck, E.Reya, M.Stratmann, W.Vogelsang: PRD53, 4775 (1996) C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) C L=1 Standard NLO C 2 Standard LO C 3 Valence NLO C 4 Valence LO IF(L.ge.3) GOTO 3 C..."Standard" scenario 1 IF(L.ge.2) GOTO 2 A = 1.1 B = 4.3 W =10.47 GOTO 5 2 A = 0.76 B = 6.84 W =11.80 GOTO 5 C..."Valence" scenario 3 IF(L.ge.4) GOTO 4 A = 1.05 B = 4.44 W = 9.87 GOTO 5 4 A = 0.60 B = 5.93 W = 7.40 5 GRSVg = X**A *(1.-X)**B *W RETURN END DOUBLE PRECISION FUNCTION GSg(X,L) C*******************************************************************C C Gluon distribution: GSg = x * (gluon+(x) - gluon-(x)) C C T.Gehrmann, W.J.Stirling: PRD53, 6100 (1996) C C*******************************************************************C IMPLICIT DOUBLE PRECISION(A-G,P-Z) C L=1,2,3 LO sets A,B,C C 4,5,6 NLO sets A,B,C DIMENSION A(6),B(6),C(6),R(6),W(6),D(6) DATA A /0.520, 0.524, 0.456, 0.724, 0.670, 0.425/ DATA B /9.45, 6.87, 8.72, 5.71, 5.34, 11.05/ DATA C / 0., 1., 0., 0., 1., 0. / DATA R / 0., -2., -3., 0., -2., -3. / DATA W / .19, .19, .019, .171, .163, .102 / DATA D /.0236, .0144,-.00164,.0312, .0130,.00336/ rx=dsqrt(X) GSg = X**A(L) *(1.-X)**B(L) *(1. +C(L)*X +R(L)*rx) *W(L)/D(L) RETURN END SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A) C==================================================================== C SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N C - BY G.P. LEPAGE SEPT 1976/(REV)APR 1978 C==================================================================== IMPLICIT DOUBLE PRECISION(A-G,P-Z) EXTERNAL FXN DIMENSION QRAN(10) COMMON/BVEG1/XL(10),XU(10),ACC COMMON/BVEGG/NDIM,NCALL,ITMX,NPRN COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI COMMON/BVE22/NDO,IT DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10) 1 ,KG(10),IA(10) DATA NDMX/50/,ALPH/1.5D0/,QNE/1.D0/,MDS/1/ C NDO=1 DO 1 J=1,NDIM 1 XI(1,J)=QNE C ENTRY VEGAS1(FXN,AVGI,SD,CHI2A) C - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID IT=0 SI=0.D0 SI2=SI SWGT=SI SCHI=SI C ENTRY VEGAS2(FXN,AVGI,SD,CHI2A) C - NO INITIALIZATION ND=NDMX NG=1 IF(MDS.EQ.0) GO TO 2 NG=INT((DFLOAT(NCALL)/2.D0)**(1.D0/DFLOAT(NDIM))) MDS=1 IF((2*NG-NDMX).LT.0) GO TO 2 MDS=-1 NPG=NG/NDMX+1 ND=NG/NPG NG=NPG*ND 2 K=NG**NDIM NPG=NCALL/K IF(NPG.LT.2) NPG=2 CALLS=NPG*K DXG=QNE/NG DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-QNE) XND=ND NDM=ND-1 DXG=DXG*XND XJAC=QNE/CALLS DO 3 J=1,NDIM DX(J)=XU(J)-XL(J) 3 XJAC=XJAC*DX(J) C C REBIN PRESERVING BIN DENSITY C IF(ND.EQ.NDO) GO TO 8 RC=NDO/XND DO 7 J=1,NDIM K=0 XN=0.D0 DR=XN I=K 4 K=K+1 DR=DR+QNE XO=XN XN=XI(K,J) 5 IF(RC.GT.DR) GO TO 4 I=I+1 DR=DR-RC XIN(I)=XN-(XN-XO)*DR IF(I.LT.NDM) GO TO 5 DO 6 I=1,NDM 6 XI(I,J)=XIN(I) 7 XI(ND,J)=QNE NDO=ND C 8 IF(NPRN.NE.0) WRITE(6,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND 1 ,(XL(J),XU(J),J=1,NDIM) C ENTRY VEGAS3(FXN,AVGI,SD,CHI2A) C - MAIN INTEGRATION LOOP 9 IT=IT+1 TI=0.D0 TSI=TI DO 10 J=1,NDIM KG(J)=1 DO 10 I=1,ND D(I,J)=TI 10 DI(I,J)=TI C 11 FB=0.D0 F2B=FB K=0 12 K=K+1 CALL ARAN9(QRAN,NDIM) WGT=XJAC DO 15 J=1,NDIM XN=(DFLOAT(KG(J))-QRAN(J))*DXG+QNE IA(J)=XN IF(IA(J).GT.1) GO TO 13 XO=XI(IA(J),J) RC=(XN-IA(J))*XO GO TO 14 13 XO=XI(IA(J),J)-XI(IA(J)-1,J) RC=XI(IA(J)-1,J)+(XN-IA(J))*XO 14 X(J)=XL(J)+RC*DX(J) 15 WGT=WGT*XO*XND C F=WGT F=F*FXN(X,WGT) F2=F*F FB=FB+F F2B=F2B+F2 DO 16 J=1,NDIM DI(IA(J),J)=DI(IA(J),J)+F 16 IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2 IF(K.LT.NPG) GO TO 12 C IF(FB.lt.1.d-20) FB = 0.D0 C write(6,*) 'TI = ',ti,tsi,fb,fb2 F2B=DSQRT(F2B*NPG) F2B=(F2B-FB)*(F2B+FB) TI=TI+FB TSI=TSI+F2B IF(MDS.GE.0) GO TO 18 DO 17 J=1,NDIM 17 D(IA(J),J)=D(IA(J),J)+F2B 18 K=NDIM 19 KG(K)=MOD(KG(K),NG)+1 IF(KG(K).NE.1) GO TO 11 K=K-1 IF(K.GT.0) GO TO 19 C C FINAL RESULTS FOR THIS ITERATION C TSI=TSI*DV2G TI2=TI*TI WGT=TI2/TSI c write(6,*) 'wgt ',wgt,ti2,tsi,ti SI=SI+TI*WGT SI2=SI2+TI2 SWGT=SWGT+WGT SCHI=SCHI+TI2*WGT AVGI=SI/SWGT SD=SWGT*IT/SI2 CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(DFLOAT(IT)-.999D0) SD=DSQRT(QNE/SD) C IF(NPRN.EQ.0) GO TO 21 TSI=DSQRT(TSI) WRITE(6,201) IT,TI,TSI,AVGI,SD,CHI2A c WRITE(6,*) IT,TI,TSI,AVGI,SD,CHI2A IF(NPRN.GE.0) GO TO 21 DO 20 J=1,NDIM 20 WRITE(6,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND) C C REFINE GRID C 21 DO 23 J=1,NDIM XO=D(1,J) XN=D(2,J) D(1,J)=(XO+XN)/2.D0 DT(J)=D(1,J) DO 22 I=2,NDM D(I,J)=XO+XN XO=XN XN=D(I+1,J) D(I,J)=(D(I,J)+XN)/3.D0 22 DT(J)=DT(J)+D(I,J) D(ND,J)=(XN+XO)/2.D0 23 DT(J)=DT(J)+D(ND,J) C DO 28 J=1,NDIM RC=0.D0 DO 24 I=1,ND R(I)=0. IF(D(I,J).LE.0.) GO TO 24 XO=DT(J)/D(I,J) R(I)=((XO-QNE)/XO/DLOG(XO))**ALPH 24 RC=RC+R(I) RC=RC/XND K=0 XN=0.D0 DR=XN I=K 25 K=K+1 DR=DR+R(K) XO=XN XN=XI(K,J) 26 IF(RC.GT.DR) GO TO 25 I=I+1 DR=DR-RC XIN(I)=XN-(XN-XO)*DR/R(K) IF(I.LT.NDM) GO TO 26 DO 27 I=1,NDM 27 XI(I,J)=XIN(I) 28 XI(ND,J)=QNE C IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9 200 FORMAT('0INPUT PARAMETERS FOR VEGAS: NDIM=',I3,' NCALL=',F8.0 1 /28X,' IT=',I5,' ITMX=',I5/28X,' ACC=',G9.3 2 /28X,' MDS=',I3,' ND=',I4/28X,' (XL,XU)=', 3 (T40,'( ',G12.6,' , ',G12.6,' )')) 201 FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3, 1 ': INTEGRAL =',G14.8/21X,'STD DEV =',G10.4 / 2 ' ACCUMULATED RESULTS: INTEGRAL =',G14.8 / 3 24X,'STD DEV =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4) 202 FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,' DELT I ', 1 2X,' CONV''CE ',11X,'X',7X,' DELT I ',2X,' CONV''CE ' 2 ,11X,'X',7X,' DELT I ',2X,' CONV''CE ' / 2 (' ',3G12.4,5X,3G12.4,5X,3G12.4)) RETURN END SUBROUTINE SAVE(NDIM) IMPLICIT DOUBLE PRECISION(A-G,P-Z) COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI COMMON/BVE22/NDO,IT C C STORES VEGAS DATA (UNIT 7) FOR LATER RE-INITIALIZATION C WRITE(7,200) NDO,IT,SI,SI2,SWGT,SCHI, 1 ((XI(I,J),I=1,NDO),J=1,NDIM) RETURN ENTRY RESTR(NDIM) C C ENTERS INITIALIZATION DATA FOR VEGAS C READ(7,200) NDO,IT,SI,SI2,SWGT,SCHI, 1 ((XI(I,J),I=1,NDO),J=1,NDIM) 200 FORMAT(2I8,4Z16/(5Z16)) RETURN END SUBROUTINE ARAN9(QRAN,NDIM) IMPLICIT DOUBLE PRECISION(A-G,P-Z) COMMON/SEED/NUM DIMENSION QRAN(10) DO 1 I=1,NDIM QRAN(I)=RANDOM(NUM) 1 CONTINUE RETURN END FUNCTION RANDOM(SEED) * ---------------------------------------------------------------- * Ref.: K. Park and K.W. Miller, Comm. of the ACM 31 (1988) p.1192 * Use seed = 1 as first value. * ---------------------------------------------------------------- IMPLICIT INTEGER(A-Z) DOUBLE PRECISION MINV,RANDOM SAVE PARAMETER(M=2147483647,A=16807,Q=127773,R=2836) PARAMETER(MINV=0.46566128752458d-09) HI = SEED/Q LO = MOD(SEED,Q) SEED = A*LO - R*HI IF(SEED.LE.0) SEED = SEED + M RANDOM = SEED*MINV RETURN END