C ============== PROGRAM QUPROG C ============== C--- this can be made into a subroutinep within another routine- c but you will need c the Common/Logical/Integer/DAta/ and Parameter statements c-- to be in yr routine... c and you will need the intialisation logic which follows it IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS COMMON/INPUT/alambda,flavor,qsct,qsdt,iord,mode COMMON/GRID/IQ0,IQC,IQB,NXGRI,NQGRI COMMON/QNOUGHT/Q0,ZM2,AMSR,chmq,chmq2,chmq24,Q2C,Q2B COMMON/CONTROL/READIN,HEAVY,NACT,IVFN COMMON/PARVAL/STVAL(45),PASVAL(45) COMMON/PARCOV/C(11),DUPP(11,11),DDOWN(11,11) COMMON/NORMS/UN,DN,GN,UBDBN dimension PVAL(45) LOGICAL DOERRS,READIN,HEAVY,VFN INTEGER NACT,IVFN,IORD,mode INTEGER IQ0,IQC,IQB,NXGRI,NQGRI INTEGER NFREE,NDIMS,NSYSF INTEGER I,J,II,K,L,IK c-- for zm set heavy=.false. and vfn=.false. c--for ff set heavy=.true. and vfn=.false. c--for tr set heavy=.false, and vfn=.true. c--set readin=.false. c--set doerrs=.true. only if you want full errors- it takes c--22*times as long! DATA DOERRS,READIN,HEAVY,VFN + /.TRUE.,.FALSE.,.FALSE.,.TRUE./ PARAMETER (NFREE=11) C--NFREE IS THE NUMBER OF FREE PARAMETERS, NDIMS IS THE TOTAL NUMBER C--MOST OF THESE ARE ZERO-AND ARE UNUSED- SEE STVAL ARRAY PARAMETER (NDIMS=45) INTEGER IFREE(NFREE) C-- THIS SPECIFIES WHICH OF THE PARAMETERS ARE FREE DATA IFREE/2,4,6,8,9,10,11,12,13,14,16/ DATA ZM/91.187/ DATA THETAW/0.2315/ DIMENSION PWGT(20),GUP(2),DGDPAR(NDIMS) DATA STVAL/0.5000,4.0065,0.0000,5.0402,0.5000,5.3359,0.0000, . 6.2124,0.6022,-0.2356,8.8592,6.7673,-0.2044, . 6.1582,0.0000,0.2652,0.118, . 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, . 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, . 0.0,0.0,0.0, . 0.0,0.0,0.0/ c----------------------------------------------------------------------------- c---------------------------------------------------------------------------- C-------INITIALISATION CAlls------------------------------------------------ c--try 3 way logic ffn/zm-vfn/rt-vfn IF(HEAVY)THEN IVFN=1 ELSE IVFN=0 ENDIF IF(VFN)THEN IVFN=IVFN+2 ELSE IVFN=IVFN ENDIF C IVFN=0 IS ZM-VFN, 1 IS FFN,2 IS RT-VFN, 3 IS NOT ALLOWED IF(IVFN.EQ.3)THEN WRITE(*,*)'IVFN=3 SO STOP',IVFN STOP ENDIF c--roberts/THORNE initialisation and qcdnum initialisations CALL INITIAL write(*,*)'IVFN=',IVFN c--read central fit PDF params if(ivfn.eq.0)then open(unit=05,file='parcen_zm.dat',status='old') elseif(ivfn.eq.1)then open(unit=05,file='parcen_ff.dat',status='old') elseif(ivfn.eq.2)then open(unit=05,file='parcen_tr.dat',status='old') endif read(05,57)(c(i), i=1,nfree) c read from parcen.dat 57 FORMAT(11(E12.4,1X)) c--and convert it to THE pval array c i.e fill in the other necessary fixed params DO I=1,NDIMS PVAL(I)=STVAL(I) ENDDO DO IK=1,NFREE I=IFREE(IK) PVAL(I)=C(IK) ENDDO c--first compute using only the central values of the main fit C-- do evolution over whole grid CALL GRCUTS(-1D0,-1D0,-1D0,-1D0) CALL EVOLVE(PVAL) alambda=15.400*PVAL(17)-1.5095 if(alambda.le.0.01)alambda=0.01 c-------------------------END OF initialisation-------------------------- C************************************************************ C here input x,q2 values required c********************************************************** Q2=20. c-- work out A and B coefficients for F2,XF3 NC call q2combs(q2) X=1.e-3 c---------------------------------------------------------------------------- c-- call the resullts routine C--THIS ROUTINE gives ALL QUANTITIES OF INTEREST for the current pval array CALL RESULTS(Q2,X,F2,FL,F2C,F2NC,FLNC,XF3NC + ,UQ,DQ,GQ,SQ,UB,DB,ST,CH,BT,DU + ,SIGNCEM,SIGCCEM,SIGNCEP,SIGCCEP) c----------------------------------------------------------------------------- c-- F2, FL are F2, FL for ep by gamma exchange, F2C is F2 charm c--F2NC,FLNC,XF3Nc are the ep structure functions for gamma and Z0 c--UQ,DQ,GQ,SQ are uvalence,dvalence, gluon and total sea PDfs c--UB,DB,ST,CH,BT,DU are ubar,dbar,sbar,cbar,bbar,and d/u C-SIGNCEM,SIGCCEM,SIGNCEP,SIGCCEP are the reduced cross-sections for C-neutral current e-, Charged current e-, Neutral current e+, c-Charged current e+, respectivelu C------------------------------------------------------------------------- C GLUON is GQ GCENT=GQ write(*,*)'glue central,x,q2',GCENT,x,q2 c-- at this point you have various structure functions, parton distributions c--and reduced cross-sectiosn of interest- c-- the rest of the code is to get errors on them- which takes c--2 *nfree= 22 times as long! c----------------------------------------------------------------------------- write(*,*)'calc for central values complete-now error bands' c---------------------------------------------------------------------------- ERRG=0.0 IF(DOERRS)THEN c --now get errors on all these from the 2*Neigenvector PDF sets c--there are input files for statistical (ie uncorrelated) c--and correlated systematic--and total errors c--this example is for total errors c open(unit=03,file='instat_**.dat',status='old') c open(unit=03,file='insys_**.dat',status='old') if(ivfn.eq.0)then open(unit=03,file='intot_zm.dat',status='old') elseif(ivfn.eq.1)then open(unit=03,file='intot_ff.dat',status='old') elseif(ivfn.eq.2)then open(unit=03,file='intot_tr.dat',status='old') endif 110 format(11(1x,f11.7)) 112 format(1x,'params up k=',i2,11(1x,f11.7)) 113 format(1x,'params down k=',i2,11(1x,f11.7)) DO IK=1,NFREE C-read up and down errors for each of nfree eigenvectors C-read ddup and ddown read(03,*)(dupp(i,ik),i=1,nfree) read(03,*)(ddown(i,ik),i=1,nfree) ENDDO c do ik=1,nfree c print 112, ik,(dupp(i,ik),i=1,nfree) c print 113, ik,(ddown(i,ik),i=1,nfree) c enddo C get the new PDF array for evolution c-fromm ddup or down and put it in pval C-first need to modify it to put in the fixed parameter values, just as we did c--for the central vaues DO I=1,NDIMS PVAL(I)=STVAL(I) ENDDO DO IK=1,NFREE c---assign j=1 to up errors and j=2 to down errors DO J=1,2 DO II=1,NFREE I=IFREE(II) IF(J.EQ.1)THEN PVAL(I)=DUPP(II,IK) ELSEIF(J.EQ.2)then pval(i)=DDOWN(II,IK) ENDIF ENDDO c--end FOR II LOOP c--check pval values are reasonable DO I=1,NDIMS PASVAL(I)=PVAL(I) ENDDO c--evolve within this j=1,2 loop for each of the nfree eigenvectors in turn CALL GRCUTS(-1D0,-1D0,-1D0,-1D0) CALL EVOLVE(PVAL) alambda=15.400*PVAL(17)-1.5095 if(alambda.le.0.01)alambda=0.01 c-- now ready to evaluate all the quantities again: once for each of the c--2 *Neigenvector PDF sets c-- in this routine the loop j=1,2 does the up and down shifts c--and the loop ik=1,nfree goes over the nfree eigenvectors c-- this is done for each x,q2 value of interest c--x,q2 have already been given shdnt need again c Q2=20. c CALL Q2COMBS(Q2) c x=1.e-3 CALL RESULTS(Q2,X,F2,FL,F2C,F2NC,FLNC,XF3NC + ,UQ,DQ,GQ,SQ,UB,DB,ST,CH,BT,DU + ,SIGNCEM,SIGCCEM,SIGNCEP,SIGCCEP) c--save this results for up and down GUP(J)=GQ c write(*,*)'pval end j,ik',j,ik,pval(ik) C--END 1,2 UP DOwN LOOP ENDDO C--STILL In IK=1,NFREE LOOP HERE c - now have up an down of each of the iK=1,nfree eigenvectOR shifts c-add up ((up-down)/2)**2 for each i and take sqrt for THE GLUON DGDPAR(IK) = (GUP(1)-GUP(2))**2. C--HERE IS THE IK ENDDO 89 ENDDO c--calculation of diagonal errors now c--note the dividing by two is done here- RGSKL=0.0 DO Ik=1,NFREE RGSKL=DGDPAR(IK) +RGSKL ENDDO ERRG=0.5*SQRT(RGSKL) C--THIS IS THE ENDIF FOR .DOERRS. ENDIF c------------------------------------------------------------------- write(*,*)'calculation of error bands complete-output' c------------------------------------------------------------------------ C--NOW WE'VE CONSTRUCTED THE QUANTITIES OF INTEREST ACROSS A WHOLE X,Q2 C--GRID we'll OUTPUT c--then the error output GUPDUM=GCENT+ERRG GDNDUM=GCENT-ERRG WRITE(*,*)'glu central,up, down', GCENT,GUPDUM,GDNDUM write(*,*)'x,q2',x,q2 998 CONTINUE 53 FORMAT(3(E12.4,1X)) STOP 999 END c------------------------------------------------------------------------------ SUBROUTINE INITIAL IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS COMMON/INPUT/alambda,flavor,qsct,qsdt,iord,mode COMMON/GRID/IQ0,IQC,IQB,NXGRI,NQGRI COMMON/QNOUGHT/Q0,ZM2,AMSR,chmq,chmq2,chmq24,Q2C,Q2B COMMON/CONTROL/READIN,HEAVY,NACT,IVFN COMMON/PARVAL/STVAL(45),PASVAL(45) LOGICAL READIN,HEAVY,VFN INTEGER NACT,IVFN,IORD INTEGER IQ0,IQC,IQB,NXGRI,NQGRI PARAMETER (NSTARTP=7) DIMENSION QSP(NSTARTP) DATA QSP/10.,20.,30.,40.,50.,80.,100./ DIMENSION PWGT(20) C-- Initialisation c-- stuff for dick roberts routines call wate96 iord=1 flavor=3 qsdt=7.29 qsct=74. mode=2 alambda=15.4*STVAL(17)-1.5095 if(alambda.le.0.01)alambda=0.01 c--qcdnum initialisation CALL QNINIT c--se thresholds Q0=7.0 ZM=91.187D0 ZM2=ZM*ZM ALPHAS=QNALFA(ZM2) Q2C=1.8225 Q2B=18.49 IF (Q0.LT.Q2C) THEN NACT=3 ELSE NACT=4 ENDIF c--this merely defines nact where we startevolution c--namely at q0 IF (HEAVY) NACT=3 CALL QNRSET('MCSTF',SQRT(Q2C)) CALL QNRSET('MBSTF',SQRT(Q2B)) CALL QNRSET('MCALF',SQRT(Q2C)) CALL QNRSET('MBALF',SQRT(Q2B)) IF (HEAVY) THEN CALL QTHRES(1D6,2D6) ELSE CALL QTHRES(Q2C,Q2B) ENDIF C-- x - Q2 grid DO I=1,NSTARTP CALL GRQINP(QSP(I),1) ENDDO CALL GRQINP(Q0,1) CALL GRQINP(Q2C,1) CALL GRQINP(Q2B,1) c qcdnum grid not my grid CALL GRXLIM(120,97D-8) CALL GRQLIM(61,29D-2,200D3) C-- Get final grid definitions and grid indices of Q0, Q2C and Q2B CALL GRGIVE(NXGRI,XMI,XMA,NQGRI,QMI,QMA) c WRITE(*,*)'NX,XL,XH,NQ,QL,QH',NXGRI,XMI,XMA,NQGRI,QMI,QMA IQ0 = IQFROMQ(Q0) IQC = IQFROMQ(Q2C) IQB = IQFROMQ(Q2B) C-- Allow for heavy weights IF (HEAVY) THEN CALL QNLSET('WTF2C',.TRUE.) CALL QNLSET('WTF2B',.TRUE.) CALL QNLSET('CLOWQ',.FALSE.) CALL QNLSET('WTFLC',.TRUE.) CALL QNLSET('WTFLB',.TRUE.) ENDIF C-- Compute weights and dump, or read in IF (READIN) THEN OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED', . STATUS='UNKNOWN') CALL QNREAD(24,ISTOP,IERR) ELSE CALL QNFILW(0,0) IF (HEAVY) THEN OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED', . STATUS='UNKNOWN') CALL QNDUMP(24) ENDIF ENDIF C-- Apply cuts to grid c--taking away the s cut at 600d0 CALL GRCUTS(-1D0,-1D0,-1D0,-1D0) C-- Choose renormalisation and factorisation scales CALL QNRSET('AAAR2',1D0) ! renormalisation CALL QNRSET('BBBR2',0D0) CALL QNRSET('AAM2L',1D0) ! factorisation (light) CALL QNRSET('BBM2L',0D0) CALL QNRSET('AAM2H',1D0) ! factorisation (heavy) CALL QNRSET('BBM2H',0D0) ZM=91.187D0 AS=STVAL(17) CALL QNRSET('ALFQ0',ZM*ZM) CALL QNRSET('ALFAS',AS) ZM2=ZM*ZM ALPHAS=QNALFA(ZM2) WRITE(*,*)'ALPHAS AT Mz2',ALPHAS C-- Book non-singlet distributions CALL QNBOOK(2,'UPLUS') CALL QNBOOK(3,'DPLUS') CALL QNBOOK(4,'SPLUS') CALL QNBOOK(5,'CPLUS') CALL QNBOOK(6,'BPLUS') CALL QNBOOK(7,'UPVAL') CALL QNBOOK(8,'DNVAL') C-- Book linear combinations for proton for f = 3,4,5 flavours C--FIRST THE WORLD SF DATA STUFF CALL VZERO(PWGT,20) PWGT(7)=1. PWGT(8)=1. CALL QNLINC(26,'VALENCE',3,PWGT) CALL QNLINC(26,'VALENCE',4,PWGT) CALL QNLINC(26,'VALENCE',5,PWGT) CALL VZERO(PWGT,20) PWGT(1) = 4./18. PWGT(2) = 4./9. PWGT(3) = 1./9. PWGT(4) = 1./9. CALL QNLINC(27,'PROTON',3,PWGT) PWGT(1) = 5./18. PWGT(5) = 4./9. CALL QNLINC(27,'PROTON',4,PWGT) PWGT(1) = 11./45. PWGT(6) = 1./9. CALL QNLINC(27,'PROTON',5,PWGT) PWGT(1) = 4./18. PWGT(2) = 1./9. PWGT(3) = 4./9. PWGT(4) = 1./9. CALL QNLINC(28,'NEUT',3,PWGT) PWGT(1) = 5./18. PWGT(5) = 4./9. CALL QNLINC(28,'NEUT',4,PWGT) PWGT(1) = 11./45. PWGT(6) = 1./9. CALL QNLINC(28,'NEUT',5,PWGT) CALL VZERO(PWGT,20) PWGT(1) = 4./18. PWGT(2) = 5./18. PWGT(3) = 5./18. PWGT(4) = 1./9. CALL QNLINC(29,'DEUTE',3,PWGT) PWGT(1) = 5./18. PWGT(5) = 4./9. CALL QNLINC(29,'DEUTE',4,PWGT) PWGT(1) = 11./45. PWGT(6) = 1./9. CALL QNLINC(29,'DEUTE',5,PWGT) CALL VZERO(PWGT,20) PWGT(1) = 1. PWGT(2) = 0. PWGT(3) = 0. PWGT(4) = 0. PWGT(5) = 0. PWGT(6) = 0. PWGT(7) = -1. PWGT(8) = -1. CALL QNLINC(16,'SEAQK',3,PWGT) CALL QNLINC(16,'SEAQK',4,PWGT) CALL QNLINC(16,'SEAQK',5,PWGT) C--THEN THE XSECN STUFF C--THESE COMBINATIONS ARE FOR UNPOLARISED LEPTON BEAM.. C--BUT ARE EASILY GENERALISED e.g. for e+ (1+P)* these combs c--for e- (1-P)* these combs , where P=1 for rh P=-1 for lh c--first cc combinations CALL VZERO(PWGT,20) PWGT(1)=0.5 PWGT(7)=-0.5 PWGT(8)=0.5 CALL QNLINC(12,'CCEP2',3,PWGT) CALL QNLINC(12,'CCEP2',4,PWGT) call qnlinc(12,'CCEP2',5,PWGT) PWGT(7) = 0.5 PWGT(8) = -0.5 CALL QNLINC(13,'CCEM2',3,PWGT) CALL QNLINC(13,'CCEM2',4,PWGT) CALL QNLINC(13,'CCEM2',5,PWGT) CALL VZERO(PWGT,20) c---these need cabibbo suppression factors correctly worked out with c--sigma..actually for ccep3 we are always above c threcshold cos2thc=0.95 sin2thc=0.05 PWGT(2) = -0.5 PWGT(3) = 0.5*cos2thc PWGT(4) = 0.5*sin2thc PWGT(1) = 0. PWGT(7) = +0.5 PWGT(8) = +0.5*cos2thc CALL QNLINC(14,'CCEP3',3,PWGT) PWGT(2) = -0.5 PWGT(3) = 0.5 PWGT(4) = 0.5 PWGT(1) = 0.0 PWGT(5) = -0.5 PWGT(7) = +0.5 PWGT(8) = +0.5 CALL QNLINC(14,'CCEP3',4,PWGT) c PWGT(6) = 0.5 c PWGT(1) = 0.5/5. c--leave 5 flavour the same as 4 flavour since b to t is cabibbo supressed c--rather strikingly CALL QNLINC(14,'CCEP3',5,PWGT) IF(HEAVY)THEN C--FOR HEAVY MAKE IT 4 FLAVOUR EVEN WHEN NACT=3, SO THAT ITS AT LEAST C--GOT THE CABIBBO UNSUPRESSIONS..OTHERWISe ITS JUST NOT RIGHT PWGT(2) = -0.5 PWGT(3) = 0.5 PWGT(4) = 0.5 PWGT(1) = 0.0 PWGT(5) = -0.5 PWGT(7) = +0.5 PWGT(8) = +0.5 CALL QNLINC(14,'CCEP3',3,PWGT) ENDIF CALL VZERO(PWGT,20) cos2thc=0.95 sin2thc=0.05 PWGT(2) = 0.5 PWGT(3) = -0.5*cos2thc PWGT(4) = -0.5*sin2thc PWGT(1) = 0. PWGT(7) = +0.5 PWGT(8) = 0.5*cos2thc CALL QNLINC(15,'CCEM3',3,PWGT) PWGT(2) = 0.5 PWGT(3) = -0.5 PWGT(4) = -0.5 PWGT(1) = 0.0 PWGT(7) = +0.5 PWGT(8) = +0.5 PWGT(1) = 0.0 PWGT(5) = 0.5 CALL QNLINC(15,'CCEM3',4,PWGT) C PWGT(6) = -0.5 C PWGT(1) = -0.5/NACT C--LEAVE 5 AS FOR 4 BECAUSE BBAR TO TBAR IS CABIBBO SUPRESSED CALL QNLINC(15,'CCEM3',5,PWGT) IF(HEAVY)THEN C--FOR HEAVY MAKE IT 4 FLAVOUR EVEN WHEN NACT=3, SO THAT ITS AT LEAST C--GOT THE CABIBBO UNSUPRESSIONS..OTHERWISe ITS JUST NOT RIGHT PWGT(2) = 0.5 PWGT(3) = -0.5 PWGT(4) = -0.5 PWGT(1) = 0.0 PWGT(5) = 0.5 PWGT(7) = +0.5 PWGT(8) = +0.5 CALL QNLINC(15,'CCEM3',3,PWGT) ENDIF c--define some quark pdfs CALL VZERO(PWGT,20) PWGT(2) = 0.5 PWGT(7) = -0.5 PWGT(1) = 0.5/3. CALL QNLINC(17,'UB',3,PWGT) PWGT(1) = 0.5/4. CALL QNLINC(17,'UB',4,PWGT) PWGT(1) = 0.5/5. CALL QNLINC(17,'UB',5,PWGT) CALL VZERO(PWGT,20) PWGT(4) = 0.5 PWGT(1) = 0.5/3. CALL QNLINC(18,'SB',3,PWGT) PWGT(1) = 0.5/4. CALL QNLINC(18,'SB',4,PWGT) PWGT(1) = 0.5/5. CALL QNLINC(18,'SB',5,PWGT) CALL VZERO(PWGT,20) CALL QNLINC(19,'CB',3,PWGT) PWGT(5) = 0.5 PWGT(1) = 0.5/4. CALL QNLINC(19,'CB',4,PWGT) PWGT(1) = 0.5/5. CALL QNLINC(19,'CB',5,PWGT) CALL VZERO(PWGT,20) PWGT(3) = 0.5 PWGT(8) = -0.5 PWGT(1) = 0.5/3. CALL QNLINC(20,'DB',3,PWGT) PWGT(1) = 0.5/4. CALL QNLINC(20,'DB',4,PWGT) PWGT(1) = 0.5/5. CALL QNLINC(20,'DB',5,PWGT) CALL VZERO(PWGT,20) CALL QNLINC(21,'BB',3,PWGT) CALL QNLINC(21,'BB',4,PWGT) PWGT(6) = 0.5 PWGT(1) = 0.5/5. CALL QNLINC(21,'BB',5,PWGT) c--- CALL QNRGET('MCSTF',chmq) chmq2=chmq*chmq chmq24=4.*chmq2 RETURN END C----------------------------------------------------------------------------- SUBROUTINE Q2COMBS(Q2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS COMMON/INPUT/alambda,flavor,qsct,qsdt,iord,mode COMMON/GRID/IQ0,IQC,IQB,NXGRI,NQGRI COMMON/QNOUGHT/Q0,ZM2,AMSR,chmq,chmq2,chmq24,Q2C,Q2B COMMON/CONTROL/READIN,HEAVY,NACT,IVFN COMMON/PARVAL/STVAL(45),PASVAL(45) LOGICAL READIN,HEAVY,VFN INTEGER NACT,IVFN,IORD INTEGER IQ0,IQC,IQB,NXGRI,NQGRI DIMENSION PWGT(20) C--NOW NEED TO KNOW Q2 THAT WE ARE INTERESTED IN DATA ZM/91.187/ DATA THETAW/0.2315/ c--code for ncem xsecn c--general stuff for ncem xsecn EL=-1. T3L=-0.5 AL=T3L VL=T3L-2.*EL*THETAW EU=2./3. T3U=+0.5 AU=T3U VU=T3U-2.*EU*THETAW ED=-1./3. T3D=-0.5 AD=T3D VD=T3D-2.*ED*THETAW PROP=Q2/(Q2+ZM*ZM) PZ=PROP/(4.*THETAW*(1.-THETAW)) ALU=EU*EU+PZ*2.*EL*EU*VU*(VL+AL)+ . PZ*PZ*(VL+AL)**2.*(VU*VU+AU*AU) BLU=PZ*2.*EL*EU*AU*(VL+AL) + . PZ*PZ*(VL+AL)**2.*2.*VU*AU ALD=ED*ED+PZ*2.*EL*ED*VD*(VL+AL)+ . PZ*PZ*(VL+AL)**2.*(VD*VD+AD*AD) BLD=PZ*2.*EL*ED*AD*(VL+AL) + . PZ*PZ*(VL+AL)**2.*2.*VD*AD ARU=EU*EU+PZ*2.*EL*EU*VU*(VL-AL)+ . PZ*PZ*(VL-AL)**2.*(VU*VU+AU*AU) BRU=-PZ*2.*EL*EU*AU*(VL-AL) - . PZ*PZ*(VL-AL)**2.*2.*VU*AU ARD=ED*ED+PZ*2.*EL*ED*VD*(VL-AL)+ . PZ*PZ*(VL-AL)**2.*(VD*VD+AD*AD) BRD=-PZ*2.*EL*ED*AD*(VL-AL) - . PZ*PZ*(VL-AL)**2.*2.*VD*AD C--NOW DEFINE COMBINATIONS OF LEFT AND RIGHT FOR UNPOLARISED C--AU+=(1-P)/2.*ALU +(1+P)/2.*ARU..SIMILARYLY FOR AD+,BU+BD+ AUM=0.5*(ALU+ARU) ADM=0.5*(ALD+ARD) BUM=0.5*(BLU+BRU) BDM=0.5*(BLD+BRD) C--NOW WE CAN DEFINE THE RIGHT COMBINATIONS TO GET F2,XF3 CALL VZERO(PWGT,20) PWGT(1)=(AUM+2.*ADM)/3. PWGT(2)=AUM PWGT(3)=ADM PWGT(4)=ADM CALL QNLINC(23,'NCEM2',3,PWGT) PWGT(1)=2.*(AUM+ADM)/4. PWGT(5)=AUM CALL QNLINC(23,'NCEM2',4,PWGT) PWGT(6)=ADM PWGT(1)=(2.*AUM+3.*ADM)/5. CALL QNLINC(23,'NCEM2',5,PWGT) CALL VZERO(PWGT,20) PWGT(7) = BUM PWGT(8) = BDM CALL QNLINC(25,'NCEM3',3,PWGT) CALL QNLINC(25,'NCEM3',4,PWGT) CALL QNLINC(25,'NCEM3',5,PWGT) c--code for nceP xsecn c--general stuff for nce xsecn EL=1. T3L=+0.5 AL=T3L VL=T3L-2.*EL*THETAW EU=2./3. T3U=+0.5 AU=T3U VU=T3U-2.*EU*THETAW ED=-1./3. T3D=-0.5 AD=T3D VD=T3D-2.*ED*THETAW ALU=EU*EU+PZ*2.*EL*EU*VU*(VL+AL)+ . PZ*PZ*(VL+AL)**2.*(VU*VU+AU*AU) BLU=PZ*2.*EL*EU*AU*(VL+AL) + . PZ*PZ*(VL+AL)**2.*2.*VU*AU ALD=ED*ED+PZ*2.*EL*ED*VD*(VL+AL)+ . PZ*PZ*(VL+AL)**2.*(VD*VD+AD*AD) BLD=PZ*2.*EL*ED*AD*(VL+AL) + . PZ*PZ*(VL+AL)**2.*2.*VD*AD ARU=EU*EU+PZ*2.*EL*EU*VU*(VL-AL)+ . PZ*PZ*(VL-AL)**2.*(VU*VU+AU*AU) BRU=-PZ*2.*EL*EU*AU*(VL-AL) - . PZ*PZ*(VL-AL)**2.*2.*VU*AU ARD=ED*ED+PZ*2.*EL*ED*VD*(VL-AL)+ . PZ*PZ*(VL-AL)**2.*(VD*VD+AD*AD) BRD=-PZ*2.*EL*ED*AD*(VL-AL) - . PZ*PZ*(VL-AL)**2.*2.*VD*AD C--NOW DEFINE COMBINATIONS OF LEFT AND RIGHT FOR UNPOLARISED C--AU+=(1+P)/2.*ALU +(1-P)/2.*ARU..SIMILARYLY FOR AD+,BU+BD+ AUP=0.5*(ALU+ARU) ADP=0.5*(ALD+ARD) BUP=0.5*(BLU+BRU) BDP=0.5*(BLD+BRD) C--NOW WE CAN DEFINE THE RIGHT COMBINATIONS TO GET F2,XF3 CALL VZERO(PWGT,20) PWGT(1)=(AUP+2.*ADP)/3. PWGT(2)=AUP PWGT(3)=ADP PWGT(4)=ADP CALL QNLINC(22,'NCEP2',3,PWGT) PWGT(1)=2.*(AUP+ADP)/4. PWGT(5)=AUP CALL QNLINC(22,'NCEP2',4,PWGT) PWGT(6)=ADP PWGT(1)=(2.*AUP+3.*ADP)/5. CALL QNLINC(22,'NCEP2',5,PWGT) CALL VZERO(PWGT,20) PWGT(7) = BUP PWGT(8) = BDP CALL QNLINC(24,'NCEP3',3,PWGT) CALL QNLINC(24,'NCEP3',4,PWGT) CALL QNLINC(24,'NCEP3',5,PWGT) c-- now these combinations are defined we can GET QUANTITIES OF INTEREST FOR EAC--CH X,Q2 RETURN END C ============== SUBROUTINE RESULTS(Q2,X,F2CENT,FLCENT,F2CCEN,F2NPEN,FLNCPN,XF3PNEN + ,UCENT,DCENT,GCENT,SCENT + ,UBCEN,DBCEN,STCEN,CHCEN,BTCEN,DUCEN + ,SIGNEN,SIGCEN,SIGPNEN,SIGPCEN) C ============== IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS COMMON/INPUT/alambda,flavor,qsct,qsdt,iord,mode COMMON/GRID/IQ0,IQC,IQB,NXGRI,NQGRI COMMON/QNOUGHT/Q0,ZM2,AMSR,chmq,chmq2,chmq24,Q2C,Q2B COMMON/CONTROL/READIN,HEAVY,NACT,IVFN COMMON/PARCOV/C(11),DUPP(11,11),DDOWN(11,11) COMMON/PARVAL/STVAL(45),PASVAL(45) COMMON/NORMS/UN,DN,GN,UBDBN DIMENSION TMPVAL(45) LOGICAL READIN,HEAVY,VFN INTEGER NACT,NTRIES,IVFN,IORD DATA ZM/91.187/ DATA THETAW/0.2315/ REAL GLU,UVAL,DVAL,RAT DIMENSION PWGT(20) IF(IVFN.EQ.1)THEN F2CENT=QSTFXQ('F2','PROTON',X,Q2,IFLAG)+ . QSTFXQ('F2C','PROTON',X,Q2,IFLAG)+ . QSTFXQ('F2B','PROTON',X,Q2,IFLAG) FLCENT=QSTFXQ('FL','PROTON',X,Q2,IFLAG)+ . QSTFXQ('FLC','PROTON',X,Q2,IFLAG)+ . QSTFXQ('FLB','PROTON',X,Q2,IFLAG) F2CCEN=QSTFXQ('F2C','PROTON',X,Q2,IFLAG) ELSEIF(IVFN.EQ.0)THEN F2CENT=QSTFXQ('F2','PROTON',X,Q2,IFLAG) FLCENT=QSTFXQ('FL','PROTON',X,Q2,IFLAG) F2CCEN=8./9.*QPDFXQ('CB',X,Q2,IFLAG) ELSEIF(IVFN.EQ.2)THEN f2cbit=0 c--here call the jacksmithstuff- c-- this is very time consuming even in the linking so is commented out c--only put it in if doing thorne roberts coefficient functions AND c--wanting specifically to look at f2charm c W2 = Q2*(1.-X)/X c--this is not quite w2 but corresponds to the jacsmith check c--on out of range c--mass of charm quark c if(w2.le.chmq24)go to 567 c if(Q2.lt.1.25.or.x.lt.1d-5)go to 567 c if(q2.le.chmq2)then c call js(X,Q2,chmq,fhad) c f2cbit=fhad c else c call js(X,chmq2,chmq,fhad) c f2cbit=fhad c endif c 567 continue c write(*,*)'w2,x,q2,fhad,f2cbit',w2,x,q2,fhad,f2cbit if(x.ge.0.00000112.and.x.le.0.92)then call sfun(x,q2,mode,f2p,flp,f1p,rp,f2n,fln,f1n,rn, xf2c,flc,f1c,f2b,flb,f1b) else f2p=0. flp=0. f2c=0. endif F2CENT=F2P FLCENT=FLP F2CCEN=F2C +f2cbit ENDIF UCENT=QPDFXQ('UPVAL',X,Q2,IFAIL) DCENT=QPDFXQ('DNVAL',X,Q2,IFAIL) GCENT=QPDFXQ('GLUON',X,Q2,IFAIL) UBCEN=QPDFXQ('UB',X,Q2,IFAIL) DBCEN=QPDFXQ('DB',X,Q2,IFAIL) STCEN=QPDFXQ('SB',X,Q2,IFAIL) IF(Q2.GE.Q2C)THEN CHCEN=QPDFXQ('CB',X,Q2,IFAIL) ELSE CHCEN=0.0 ENDIF IF(Q2.GE.Q2B)THEN BTCEN=QPDFXQ('BB',X,Q2,IFAIL) ELSE BTCEN=0.0 ENDIF UR=UCENT+2.*UBCEN DR=DCENT+2.*DBCEN IF(UR.GT.0.0)THEN DUCEN=DR/UR ELSE DUCEN=0.0 ENDIF SCENT=2.*(UBCEN+DBCEN+STCEN+CHCEN+BTCEN) c-- a new patch of code for cc red xsecn here c--notE ccem and ccep code will NOT work for fixed flavour number c--QCDNUM was NOT built for adding on the charmed quark contribution to CC C--PROCESSES IF(HEAVY)THEN XF3=0. F2=0. FL=0. ELSE XF3=QSTFXQ('XF3','CCEM3',X,Q2,IFLAG) F2=QSTFXQ('F2','CCEM2',X,Q2,IFLAG) FL=QSTFXQ('FL','CCEM2',X,Q2,IFLAG) ENDIF XF3CEN=XF3 FLCCMN=FL F2CMEN=F2 C---- and nc em red xsec is needed too XF3=QSTFXQ('XF3','NCEM3',X,Q2,IFLAG) XF3NEN=XF3 IF(IVFN.EQ.1)THEN F2=QSTFXQ('F2','NCEM2',X,Q2,IFLAG)+ . QSTFXQ('F2C','PROTON',X,Q2,IFLAG)+ . QSTFXQ('F2B','PROTON',X,Q2,IFLAG) FL=QSTFXQ('FL','NCEM2',X,Q2,IFLAG)+ . QSTFXQ('FLC','PROTON',X,Q2,IFLAG)+ . QSTFXQ('FLB','PROTON',X,Q2,IFLAG) ELSEIF(IVFN.EQ.0)THEN F2=QSTFXQ('F2','NCEM2',X,Q2,IFLAG) FL=QSTFXQ('FL','NCEM2',X,Q2,IFLAG) ELSEIF(IVFN.EQ.2)THEN F2Z=QSTFXQ('F2','NCEM2',X,Q2,IFLAG) FLZ=QSTFXQ('FL','NCEM2',X,Q2,IFLAG) F2GAM=QSTFXQ('F2','PROTON',X,Q2,IFLAG) FLGAM=QSTFXQ('FL','PROTON',X,Q2,IFLAG) if(x.ge.0.00000112.and.x.le.0.92)then call sfun(x,q2,mode,f2p,flp,f1p,rp,f2n,fln,f1n,rn, xf2c,flc,f1c,f2b,flb,f1b) else f2p=0. flp=0. endif IF(F2GAM.GT.0.0.AND.FLGAM.GT.0.0)THEN F2=F2P*F2Z/F2GAM FL=FLP*FLZ/FLGAM ELSE F2=0. FL=0. ENDIF c--ivfn endif ENDIF F2NMEN=F2 FLNCMN=FL SNOM=4.*920.*27.5 Y=Q2/(SNOM*X) c--put on kinematic limits for actually reduced cross-sections IF(Y.LT.1.0)THEN YP=1 + (1-Y)*(1-Y) YM=1 - (1-Y)*(1-Y) SIGCEN=(F2CMEN*YP - Y*Y*FLCCMN + YM*XF3CEN)*0.5 SIGNEN=(F2NMEN*YP - Y*Y*FLNCMN + YM*XF3NEN)/YP ELSE SIGCEN=0. SIGNEN=0. C--Y KINEMATIC REION ENDIF ENDIF c-- a new patch of code for cc eP red xsecn here c--note ccem and ccep code will NOT work for fixed flavour number c--QCDNUM was NOT built for adding on the charmed quark contribution to C CC PROCESSES IF(HEAVY)THEN XF3=0. F2=0. FL=0. ELSE XF3=QSTFXQ('XF3','CCEP3',X,Q2,IFLAG) F2=QSTFXQ('F2','CCEP2',X,Q2,IFLAG) FL=QSTFXQ('FL','CCEP2',X,Q2,IFLAG) ENDIF XF3PCEN=XF3 FLCCPN=FL F2CPEN=F2 C---- and nc eP red xsec is needed too XF3=QSTFXQ('XF3','NCEP3',X,Q2,IFLAG) XF3PNEN=XF3 IF(IVFN.EQ.1)THEN F2=QSTFXQ('F2','NCEP2',X,Q2,IFLAG)+ . QSTFXQ('F2C','PROTON',X,Q2,IFLAG)+ . QSTFXQ('F2B','PROTON',X,Q2,IFLAG) FL=QSTFXQ('FL','NCEP2',X,Q2,IFLAG)+ . QSTFXQ('FLC','PROTON',X,Q2,IFLAG)+ . QSTFXQ('FLB','PROTON',X,Q2,IFLAG) ELSEIF(IVFN.EQ.0)THEN F2=QSTFXQ('F2','NCEP2',X,Q2,IFLAG) FL=QSTFXQ('FL','NCEP2',X,Q2,IFLAG) ELSEIF(IVFN.EQ.2)THEN F2Z=QSTFXQ('F2','NCEP2',X,Q2,IFLAG) FLZ=QSTFXQ('FL','NCEP2',X,Q2,IFLAG) F2GAM=QSTFXQ('F2','PROTON',X,Q2,IFLAG) FLGAM=QSTFXQ('FL','PROTON',X,Q2,IFLAG) if(x.ge.0.00000112.and.x.le.0.92)then call sfun(x,q2,mode,f2p,flp,f1p,rp,f2n,fln,f1n,rn, xf2c,flc,f1c,f2b,flb,f1b) else f2p=0. flp=0. endif IF(F2GAM.GT.0.0.AND.FLGAM.GT.0.0)THEN F2=F2P*F2Z/F2GAM FL=FLP*FLZ/FLGAM ELSE F2=0. FL=0. ENDIF c--ivfn endif ENDIF F2NPEN=F2 FLNCPN=FL c--y kin limit C SNOM=4.*820.*27.5 SNOM=4.*920.*27.5 Y=Q2/(SNOM*X) IF(Y.LT.1.0)THEN YP=1 + (1-Y)*(1-Y) YM=1 - (1-Y)*(1-Y) SIGPCEN=(F2CPEN*YP - Y*Y*FLCCPN - YM*XF3PCEN)*0.5 SIGPNEN=(F2NPEN*YP - Y*Y*FLNCPN - YM*XF3PNEN)/YP ELSE SIGPCEN=0. SIGPNEN=0. ALEXCEN=0. C--Y KINEMATIC REION ENDIF ENDIF return end C------------------------------------------------------------------------ SUBROUTINE EVOLVE(XVAL) C------------------------------------------------------------------------ IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION XVAL(45) COMMON/GRID/IQ0,IQC,IQB,NXGRI,NQGRI COMMON/QNOUGHT/Q0,ZM2,AMSR,chmq,chmq2,chmq24,Q2C,Q2B COMMON/CONTROL/READIN,HEAVY,NACT,IVFN COMMON/NORMS/UN,DN,GN,UBDBN LOGICAL READIN,HEAVY INTEGER NACT,IVFN ALPHAS=QNALFA(ZM2) UA=XVAL(1) UB=XVAL(2) UE=XVAL(3) UC=XVAL(4) DA=UA DB=XVAL(6) DE=XVAL(7) DC=XVAL(8) SN=XVAL(9) SA=XVAL(10) SB=XVAL(11) SE=0. SC=XVAL(12) GA=XVAL(13) GB=XVAL(14) GE=0. GC=XVAL(15) DLN=XVAL(16) DLA=0.5 DLB=XVAL(11)+2. DLE=0. DLC=0. AS=XVAL(17) CALL QNRSET('ALFAS',AS) C-- Input quark distributions at Q2 = Q02 GeV2 UN=2./AREA(UA-1,UB,UE,UC) DN=1./AREA(DA-1,DB,DE,DC) c UBDBN=DLN/AREA(DLA-1,DLB,DLE,DLC) DO IX = 1,NXGRI X = XFROMIX(IX) UPVAL=UN*FF(X,UA,UB,UE,UC) DNVAL=DN*FF(X,DA,DB,DE,DC) SEA=SN*FF(X,SA,SB,SE,SC) GN=(1-UN*AREA(UA,UB,UE,UC)- . DN*AREA(DA,DB,DE,DC)- . SN*AREA(SA,SB,SE,SC))/AREA(GA,GB,GE,GC) GLUON=GN*FF(X,GA,GB,GE,GC) c UMD=UBDBN*FF(X,DLA,DLB,DLE,DLC) UMD=DLN*FF(X,DLA,DLB,DLE,DLC) SINGL=UPVAL+DNVAL+SEA CSEA=0.0 SSEA=0.2*SEA USEA=(SEA-SSEA-CSEA)/2.-UMD DSEA=(SEA-SSEA-CSEA)/2.+UMD UPLUS=UPVAL+USEA-1./NACT*SINGL DPLUS=DNVAL+DSEA-1./NACT*SINGL SPLUS=SSEA-1./NACT*SINGL CALL QNPSET('GLUON',IX,IQ0,GLUON) CALL QNPSET('SINGL',IX,IQ0,SINGL) CALL QNPSET('UPLUS',IX,IQ0,UPLUS) CALL QNPSET('DPLUS',IX,IQ0,DPLUS) CALL QNPSET('SPLUS',IX,IQ0,SPLUS) CALL QNPSET('UPVAL',IX,IQ0,UPVAL) CALL QNPSET('DNVAL',IX,IQ0,DNVAL) ENDDO C--THINGS ARE FINE FOR HEAVY SO DO IT IF (HEAVY) THEN C-- Evolve over whole grid CALL EVOLSG(IQ0,1,NQGRI) CALL EVPLUS('UPLUS',IQ0,1,NQGRI) CALL EVPLUS('DPLUS',IQ0,1,NQGRI) CALL EVPLUS('SPLUS',IQ0,1,NQGRI) CALL EVOLNM('UPVAL',IQ0,1,NQGRI) CALL EVOLNM('DNVAL',IQ0,1,NQGRI) ELSE C-- Evolve gluon, singlet and valence over whole grid CALL EVOLSG(IQ0,1,NQGRI) CALL EVOLNM('UPVAL',IQ0,1,NQGRI) CALL EVOLNM('DNVAL',IQ0,1,NQGRI) C-- Be more careful with plus distributions IF (NACT.EQ.3) THEN C--THINGS ARE ALSO FINE IF 1Q0 IS BELOW 1QC THEN CLEARLY CSEA=0. IS OK C--SITUATION CD BE 1