1 subroutine ZEUSevolve(x,Q,pdf)
2 implicit double precision (a-h,o-z)
3 include 'parmsetup.inc'
5 character*16 name(nmxset)
6 integer nmem(nmxset),ndef(nmxset),mem
7 common/NAME/name,nmem,ndef,mem
8 integer nset,iset,isetlast
11 real*8 mc,mc2,mb,mb2,mt,mt2
12 real*8 f(-6:6),pdf(-6:6)
15 DIMENSION QSP(NSTARTP)
16 DATA QSP/10.,20.,30.,40.,50.,80.,100./
23 c print *,'iset=',iset,' now calling get_pdfqcd'
24 if(iset.ne.isetlast) then
30 UCENT=QPDFXQ('UPVAL',X,Q2,IFAIL)
31 DCENT=QPDFXQ('DNVAL',X,Q2,IFAIL)
32 GCENT=QPDFXQ('GLUON',X,Q2,IFAIL)
33 UBCEN=QPDFXQ('UB',X,Q2,IFAIL)
34 DBCEN=QPDFXQ('DB',X,Q2,IFAIL)
35 STCEN=QPDFXQ('SB',X,Q2,IFAIL)
37 CHCEN=QPDFXQ('CB',X,Q2,IFAIL)
42 BTCEN=QPDFXQ('BB',X,Q2,IFAIL)
62 c===========================================================================
63 entry ZEUSalfa(alfas,Q)
69 alfas=QALFAS(Q2,Qlam,nf,iflag)
70 c print *,q2,alfas,qlam,nf,iflag
73 c===========================================================================
75 read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
78 c===========================================================================
79 entry ZEUSinit(nset,Eorder,Q2fit)
81 if(name(nset).eq.'QCDNUM_ZEUS_TR') then
84 else if(name(nset).eq.'QCDNUM_ZEUS_FF') then
87 else if(name(nset).eq.'QCDNUM_ZEUS_ZM') then
91 print *,'name/scheme not recognized'
94 c--try 3 way logic ffn/zm-vfn/rt-vfn
105 C IVFN=0 IS ZM-VFN, 1 IS FFN,2 IS RT-VFN, 3 IS NOT ALLOWED
108 WRITE(*,*)'IVFN=3 SO STOP',IVFN
114 c--qcdnum initialisation
122 call getQmassM(nset,4,mc)
124 call getQmassM(nset,5,mb)
126 call getQmassM(nset,6,mt)
140 c--this merely defines nact where we startevolution
144 CALL QNRSET('MCSTF',SQRT(Q2C))
145 CALL QNRSET('MBSTF',SQRT(Q2B))
146 CALL QNRSET('MCALF',SQRT(Q2C))
147 CALL QNRSET('MBALF',SQRT(Q2B))
150 CALL QTHRES(1D10,2D10)
151 c CALL QTHRES(1D6,2D6)
157 CALL GRQINP(QSP(I),1)
162 c qcdnum grid not my grid
164 c CALL GRXLIM(120,97D-8)
167 c CALL GRQLIM(61,29D-2,200D3)
168 CALL GRQLIM(nq,qmin,qmax)
170 C-- Get final grid definitions and grid indices of Q0, Q2C and Q2B
172 CALL GRGIVE(NXGRI,XMI,XMA,NQGRI,QMI,QMA)
173 c WRITE(*,*)'NX,XL,XH,NQ,QL,QH',NXGRI,XMI,XMA,NQGRI,QMI,QMA
177 C-- Allow for heavy weights
180 CALL QNLSET('WTF2C',.TRUE.)
181 CALL QNLSET('WTF2B',.TRUE.)
182 CALL QNLSET('CLOWQ',.FALSE.)
183 CALL QNLSET('WTFLC',.TRUE.)
184 CALL QNLSET('WTFLB',.TRUE.)
187 C-- Compute weights and dump, or read in
190 c OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
191 c . STATUS='UNKNOWN')
192 c CALL QNREAD(24,ISTOP,IERR)
196 c OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
197 c . STATUS='UNKNOWN')
203 if (index(gridname,'none').eq.1) then
207 open(unit=2,status='old',file=gridname,
208 . form='unformatted',err=1)
209 call QNREAD(2,1,qnerr)
212 write(*,*) 'Grid file problem: ',gridname
214 write(*,*) 'Grid file does not exist'
215 write(*,*) 'Calculating and creating grid file'
217 open(unit=2,status='unknown',file=gridname,
218 . form='unformatted')
222 write(*,*) 'Existing grid file is inconsistent'
224 . write(*,*) 'Defined grid different'
226 . write(*,*) 'Heavy quark weight table different'
228 . write(*,*) 'Charm mass different'
230 . write(*,*) 'Bottom mass different'
236 C-- Apply cuts to grid
237 c--taking away the s cut at 600d0
238 CALL GRCUTS(-1D0,-1D0,-1D0,-1D0)
243 C-- Choose renormalisation and factorisation scales
245 CALL QNRSET('AAAR2',1D0) ! renormalisation
246 CALL QNRSET('BBBR2',0D0)
247 CALL QNRSET('AAM2L',1D0) ! factorisation (light)
248 CALL QNRSET('BBM2L',0D0)
249 CALL QNRSET('AAM2H',1D0) ! factorisation (heavy)
250 CALL QNRSET('BBM2H',0D0)
255 c -- only need call to listPDF here to get alphas
256 call listPDF(nset,imem,xval)
260 CALL QNRSET('ALFQ0',ZM*ZM)
261 CALL QNRSET('ALFAS',AS)
265 c WRITE(*,*)'ALPHAS AT Mz2',ALPHAS
267 C-- Book non-singlet distributions
269 CALL QNBOOK(2,'UPLUS')
270 CALL QNBOOK(3,'DPLUS')
271 CALL QNBOOK(4,'SPLUS')
272 CALL QNBOOK(5,'CPLUS')
273 CALL QNBOOK(6,'BPLUS')
274 CALL QNBOOK(7,'UPVAL')
275 CALL QNBOOK(8,'DNVAL')
277 C-- Book linear combinations for proton for f = 3,4,5 flavours
279 c--define some quark pdfs
286 CALL QNLINC(17,'UB',3,PWGT)
288 CALL QNLINC(17,'UB',4,PWGT)
290 CALL QNLINC(17,'UB',5,PWGT)
295 CALL QNLINC(18,'SB',3,PWGT)
297 CALL QNLINC(18,'SB',4,PWGT)
299 CALL QNLINC(18,'SB',5,PWGT)
301 CALL QNLINC(19,'CB',3,PWGT)
305 CALL QNLINC(19,'CB',4,PWGT)
307 CALL QNLINC(19,'CB',5,PWGT)
314 CALL QNLINC(20,'DB',3,PWGT)
316 CALL QNLINC(20,'DB',4,PWGT)
318 CALL QNLINC(20,'DB',5,PWGT)
320 CALL QNLINC(21,'BB',3,PWGT)
321 CALL QNLINC(21,'BB',4,PWGT)
325 CALL QNLINC(21,'BB',5,PWGT)
330 c==========================================================================
339 call getnmem(nset,imem)
340 call listPDF(nset,imem,xval)
373 CALL QNRSET('ALFAS',AS)
374 C-- Input quark distributions at Q2 = Q02 GeV2
375 C-- WRITE IN ACTUAL VALUES TO SAVE USING dgamma
376 C UN=2./AREA(UA-1,UB,UE,UC)
377 C DN=1./AREA(DA-1,DB,DE,DC)
380 c UBDBN=DLN/AREA(DLA-1,DLB,DLE,DLC)
381 c call grgive(nxgri,xxmin,xxmax,nqgri,qqmin,qqmax)
385 UPVAL=UN*FF_LHA(Xx,UA,UB,UE,UC)
388 DNVAL=DN*FF_LHA(Xx,DA,DB,DE,DC)
391 SEA=SN*FF_LHA(Xx,SA,SB,SE,SC)
393 C GN=(1-UN*AREA(UA,UB,UE,UC)-
394 C . DN*AREA(DA,DB,DE,DC)-
395 C . SN*AREA(SA,SB,SE,SC))/AREA(GA,GB,GE,GC)
397 GLUON=GN*FF_LHA(Xx,GA,GB,GE,GC)
401 c UMD=UBDBN*FF_LHA(X,DLA,DLB,DLE,DLC)
402 UMD=DLN*FF_LHA(Xx,DLA,DLB,DLE,DLC)
403 SINGL=UPVAL+DNVAL+SEA
405 c print *,un,dn,sn,gn,dln
409 USEA=(SEA-SSEA-CSEA)/2.-UMD
410 DSEA=(SEA-SSEA-CSEA)/2.+UMD
411 UPLUS=UPVAL+USEA-1./NACT*SINGL
412 DPLUS=DNVAL+DSEA-1./NACT*SINGL
413 SPLUS=SSEA-1./NACT*SINGL
416 CALL QNPSET('GLUON',IX,IQ0,GLUON)
417 CALL QNPSET('SINGL',IX,IQ0,SINGL)
418 CALL QNPSET('UPLUS',IX,IQ0,UPLUS)
419 CALL QNPSET('DPLUS',IX,IQ0,DPLUS)
420 CALL QNPSET('SPLUS',IX,IQ0,SPLUS)
421 CALL QNPSET('UPVAL',IX,IQ0,UPVAL)
422 CALL QNPSET('DNVAL',IX,IQ0,DNVAL)
424 C--THINGS ARE FINE FOR HEAVY SO DO IT
427 C-- Evolve over whole grid
429 CALL EVOLSG(IQ0,1,NQGRI)
430 CALL EVPLUS('UPLUS',IQ0,1,NQGRI)
431 CALL EVPLUS('DPLUS',IQ0,1,NQGRI)
432 CALL EVPLUS('SPLUS',IQ0,1,NQGRI)
433 CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
434 CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
438 C-- Evolve gluon, singlet and valence over whole grid
440 CALL EVOLSG(IQ0,1,NQGRI)
441 CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
442 CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
444 C-- Be more careful with plus distributions
447 C--THINGS ARE ALSO FINE IF 1Q0 IS BELOW 1QC THEN CLEARLY CSEA=0. IS OK
448 C--SITUATION CD BE 1<Q0<Q2C<Q2B ETC
449 C--GO DOWN QO TO 1 UP Q0 TO Q2C
450 CALL EVPLUS('UPLUS',IQ0,1,IQC)
451 CALL EVPLUS('DPLUS',IQ0,1,IQC)
452 CALL EVPLUS('SPLUS',IQ0,1,IQC)
453 C--DEAL WITH CHARM THRESH
455 CALL QADDSI('UPLUS',IQC,FACTOR)
456 CALL QADDSI('DPLUS',IQC,FACTOR)
457 CALL QADDSI('SPLUS',IQC,FACTOR)
460 CALL QADDSI('CPLUS',IQC,FACTOR)
462 CALL EVPLUS('UPLUS',IQC,IQC,IQB)
463 CALL EVPLUS('DPLUS',IQC,IQC,IQB)
464 CALL EVPLUS('SPLUS',IQC,IQC,IQB)
465 CALL EVPLUS('CPLUS',IQC,IQC,IQB)
467 ELSEIF(NACT.EQ.4)THEN
469 C--NOW WE NEED A RETHINK OF THE INITIAL CONDITIONS
470 C--FIRST DEAL WITH CPLUS TRUNING ON AT IQC
471 C--AND GOING UP TO IQB (MIDDLE AGAIN)
474 CALL QADDSI('CPLUS',IQC,FACTOR)
475 CALL EVPLUS('CPLUS',IQC,IQC,IQB)
476 C--REDIFINE THE PLUS DUSTNS TAKING INTO ACCOUNT CPLUS
480 C--NOW YOU NEED A DO LOOP AGAIN
483 UPVAL=UN*FF_LHA(Xx,UA,UB,UE,UC)
485 DNVAL=DN*FF_LHA(Xx,DA,DB,DE,DC)
486 SEA=SN*FF_LHA(Xx,SA,SB,SE,SC)
487 SINGL=UPVAL+DNVAL+SEA
488 UMD=DLN*FF_LHA(Xx,DLA,DLB,DLE,DLC)
490 CSEA=QPDFIJ('CPLUS',IX,IQ0,IFL) + 1./NACT*SINGL
491 USEA=(SEA-SSEA-CSEA)/2.-UMD
492 DSEA=(SEA-SSEA-CSEA)/2.+UMD
493 UPLUS=UPVAL+USEA-1./NACT*SINGL
494 DPLUS=DNVAL+DSEA-1./NACT*SINGL
495 SPLUS=SSEA-1./NACT*SINGL
496 CALL QNPSET('UPLUS',IX,IQ0,UPLUS)
497 CALL QNPSET('DPLUS',IX,IQ0,DPLUS)
498 CALL QNPSET('SPLUS',IX,IQ0,SPLUS)
500 C-NOW DO MIDDLE FOR UPLUS,DPLUS,SPLUS
501 CALL EVPLUS('UPLUS',IQ0,IQC,IQB)
502 CALL EVPLUS('DPLUS',IQ0,IQC,IQB)
503 CALL EVPLUS('SPLUS',IQ0,IQC,IQB)
504 C--THEN GO DOWN IQC TO 1
507 CALL QADDSI('UPLUS',IQC,FACTOR)
508 CALL QADDSI('DPLUS',IQC,FACTOR)
509 CALL QADDSI('SPLUS',IQC,FACTOR)
510 CALL EVPLUS('UPLUS',IQC,1,IQC)
511 CALL EVPLUS('DPLUS',IQC,1,IQC)
513 CALL EVPLUS('SPLUS',IQC,1,IQC)
517 C--THEN DEAL WITH B THRESHOLD FOR ALL4
519 CALL QADDSI('UPLUS',IQB,FACTOR)
520 CALL QADDSI('DPLUS',IQB,FACTOR)
521 CALL QADDSI('SPLUS',IQB,FACTOR)
522 CALL QADDSI('CPLUS',IQB,FACTOR)
525 CALL EVPLUS('UPLUS',IQB,IQB,NQGRI)
526 CALL EVPLUS('DPLUS',IQB,IQB,NQGRI)
527 CALL EVPLUS('SPLUS',IQB,IQB,NQGRI)
528 CALL EVPLUS('CPLUS',IQB,IQB,NQGRI)
530 cC--THEN DEAL WITH B TURNING ON AT IQB AND GO UP
533 CALL QADDSI('BPLUS',IQB,FACTOR)
534 CALL EVPLUS('BPLUS',IQB,IQB,NQGRI)
538 call save_pdfqcd(iset)
543 C------------------------------------------------------------------------
544 c DOUBLE PRECISION FUNCTION AREA(A,B,E,C)
545 C------------------------------------------------------------------------
547 c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
549 c IF (A.LE.-0.99.OR.B.LE.-0.99) THEN
554 c AR1=(DGAMMA_LHA(A+1)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2)
555 c AR2=E*(DGAMMA_LHA(A+1.5)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2.5)
556 c AR3=C*(DGAMMA_LHA(A+2)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+3)
559 c IF (AREA.LE.1D-6) AREA=1D-6
564 C------------------------------------------------------------------------
565 DOUBLE PRECISION FUNCTION FF_LHA(X,A,B,E,C)
566 C------------------------------------------------------------------------
568 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
570 FF_LHA=X**A*(1D0-X)**B*(1+E*SQRT(X)+C*X)
574 C------------------------------------------------------------------------
575 SUBROUTINE DVZERO (A,N)
577 C CERN PROGLIB# F121 VZERO .VERSION KERNFOR 4.40 940929
578 C ORIG. 01/07/71, modif. 24/05/87 to set integer zero
579 C modif. 25/05/94 to depend on QINTZERO
581 implicit real*8 (a-h,o-z)