LHAPDF/lhapdf5.5.1/src/wrapzeus.f

   1 ! -*- F90 -*-
   2
   3
   4       subroutine ZEUSevolve(x,Q,pdf)
   5       implicit double precision (a-h,o-z)
   6       include 'parmsetup.inc'
   7       character*64 gridname
   8       character*16 name(nmxset)
   9       integer nmem(nmxset),ndef(nmxset),mem
  10       common/NAME/name,nmem,ndef,mem
  11       integer nset,iset,isetlast
  12       data isetlast/-1/
  13       integer Eorder
  14       real*8 mc,mc2,mb,mb2,mt,mt2
  15       real*8 f(-6:6),pdf(-6:6)
  16       integer qnerr
  17       parameter(nstartp=7)
  18       DIMENSION QSP(NSTARTP)
  19       DATA QSP/10.,20.,30.,40.,50.,80.,100./
  20       real*8 xval(45)
  21       logical HEAVY,VFN
  22       real*8 pwgt(20)
  23       save
  24 !
  25       call getnset(iset)
  26 !      print *,'iset=',iset,' now calling get_pdfqcd'
  27       if(iset.ne.isetlast) then
  28         call get_pdfqcd(iset)
  29         isetlast = iset
  30       endif
  31 !
  32       Q2=Q*Q
  33       UCENT=QPDFXQ('UPVAL',X,Q2,IFAIL)
  34       DCENT=QPDFXQ('DNVAL',X,Q2,IFAIL)
  35       GCENT=QPDFXQ('GLUON',X,Q2,IFAIL)
  36       UBCEN=QPDFXQ('UB',X,Q2,IFAIL)
  37       DBCEN=QPDFXQ('DB',X,Q2,IFAIL)
  38       STCEN=QPDFXQ('SB',X,Q2,IFAIL)
  39       IF(Q2.GE.Q2C)THEN
  40         CHCEN=QPDFXQ('CB',X,Q2,IFAIL)
  41       ELSE
  42         CHCEN=0.0
  43       ENDIF
  44       IF(Q2.GE.Q2B)THEN
  45         BTCEN=QPDFXQ('BB',X,Q2,IFAIL)
  46       ELSE
  47         BTCEN=0.0
  48       ENDIF
  49       pdf(1)=dcent+dbcen
  50       pdf(2)=ucent+ubcen
  51       pdf(3)=stcen
  52       pdf(4)=chcen
  53       pdf(5)=btcen
  54       pdf(6)=0d0
  55       pdf(0)=gcent
  56       pdf(-1)=dbcen
  57       pdf(-2)=ubcen
  58       pdf(-3)=stcen
  59       pdf(-4)=chcen
  60       pdf(-5)=btcen
  61       pdf(-6)=0d0
  62 !
  63       return
  64 !
  65 !=======================================================================
  66       entry ZEUSalfa(alfas,Q)
  67       Q2=Q*Q
  68       nf=6
  69       if (Q2.lt.mt2) nf=5
  70       if (Q2.lt.mb2) nf=4
  71       if (Q2.lt.mc2) nf=3
  72       alfas=QALFAS(Q2,Qlam,nf,iflag)
  73 !      print *,q2,alfas,qlam,nf,iflag
  74       return
  75 !
  76 !=======================================================================
  77       entry ZEUSread(nset)
  78          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
  79       return
  80 !
  81 !=======================================================================
  82       entry ZEUSinit(nset,Eorder,Q2fit)
  83 !      print *,name(nset)
  84       if(name(nset).eq.'QCDNUM_ZEUS_TR') then
  85          HEAVY = .FALSE.
  86          VFN = .TRUE.
  87       else if(name(nset).eq.'QCDNUM_ZEUS_FF') then
  88          HEAVY = .TRUE.
  89          VFN = .FALSE.
  90       else if(name(nset).eq.'QCDNUM_ZEUS_ZM') then
  91          HEAVY = .FALSE.
  92          VFN = .FALSE.
  93       else
  94         print *,'name/scheme not recognized'
  95         stop 1
  96       endif
  97 !--try 3 way logic ffn/zm-vfn/rt-vfn
  98       IF(HEAVY)THEN
  99       IVFN=1
 100       ELSE
 101       IVFN=0
 102       ENDIF
 103       IF(VFN)THEN
 104       IVFN=IVFN+2
 105       ELSE
 106       IVFN=IVFN
 107       ENDIF
 108 ! IVFN=0 IS ZM-VFN, 1 IS FFN,2 IS RT-VFN, 3 IS NOT ALLOWED
 109
 110       IF(IVFN.EQ.3)THEN
 111       WRITE(*,*)'IVFN=3 SO STOP',IVFN
 112       STOP
 113       ENDIF
 114
 115
 116
 117 !--qcdnum initialisation
 118       CALL QNINIT
 119 !--se thresholds
 120       Q0=Q2fit
 121       ZM=91.187D0
 122       ZM2=ZM*ZM
 123       ALPHAS=QNALFA(ZM2)
 124
 125       call getQmassM(nset,4,mc)
 126       mc2=mc*mc
 127       call getQmassM(nset,5,mb)
 128       mb2=mb*mb
 129       call getQmassM(nset,6,mt)
 130       mt2=mt*mt
 131
 132 !      Q2C=1.8225
 133       Q2C=mc2
 134 !      Q2B=18.49
 135       Q2B=mb2
 136 !      print  *,q2c,q2b
 137
 138       IF (Q0.LT.Q2C) THEN
 139         NACT=3
 140       ELSE
 141         NACT=4
 142       ENDIF
 143 !--this merely defines nact where we startevolution
 144 !--namely at q0
 145       IF (HEAVY) NACT=3
 146
 147       CALL QNRSET('MCSTF',SQRT(Q2C))
 148       CALL QNRSET('MBSTF',SQRT(Q2B))
 149       CALL QNRSET('MCALF',SQRT(Q2C))
 150       CALL QNRSET('MBALF',SQRT(Q2B))
 151
 152       IF (HEAVY) THEN
 153         CALL QTHRES(1D10,2D10)
 154 !        CALL QTHRES(1D6,2D6)
 155       ELSE
 156         CALL QTHRES(Q2C,Q2B)
 157       ENDIF
 158
 159       DO I=1,NSTARTP
 160         CALL GRQINP(QSP(I),1)
 161       ENDDO
 162       CALL GRQINP(Q0,1)
 163       CALL GRQINP(Q2C,1)
 164       CALL GRQINP(Q2B,1)
 165 !      qcdnum grid not my grid
 166
 167 !      CALL GRXLIM(120,97D-8)
 168       CALL GRXLIM(nx,xmin)
 169
 170 !      CALL GRQLIM(61,29D-2,200D3)
 171       CALL GRQLIM(nq,qmin,qmax)
 172
 173 !--   Get final grid definitions and grid indices of Q0, Q2C and Q2B
 174
 175       CALL GRGIVE(NXGRI,XMI,XMA,NQGRI,QMI,QMA)
 176 !      WRITE(*,*)'NX,XL,XH,NQ,QL,QH',NXGRI,XMI,XMA,NQGRI,QMI,QMA
 177       IQ0 = IQFROMQ(Q0)
 178       IQC = IQFROMQ(Q2C)
 179       IQB = IQFROMQ(Q2B)
 180 !--   Allow for heavy weights
 181
 182       IF (HEAVY) THEN
 183         CALL QNLSET('WTF2C',.TRUE.)
 184         CALL QNLSET('WTF2B',.TRUE.)
 185         CALL QNLSET('CLOWQ',.FALSE.)
 186         CALL QNLSET('WTFLC',.TRUE.)
 187         CALL QNLSET('WTFLB',.TRUE.)
 188       ENDIF
 189
 190 !--   Compute weights and dump, or read in
 191 !
 192 !      IF (READIN) THEN
 193 !        OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
 194 !     .                                  STATUS='UNKNOWN')
 195 !        CALL QNREAD(24,ISTOP,IERR)
 196 !      ELSE
 197 !        CALL QNFILW(0,0)
 198 !        IF (HEAVY) THEN
 199 !          OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
 200 !     .                                    STATUS='UNKNOWN')
 201 !          CALL QNDUMP(24)
 202 !        ENDIF
 203 !      ENDIF
 204
 205
 206       if (index(gridname,'none').eq.1) then
 207          call qnfilw(0,0)
 208       else
 209          qnerr=-1
 210          open(unit=2,status='old',file=gridname,                        &
 211      &        form='unformatted',err=1)
 212          call QNREAD(2,1,qnerr)
 213     1    close(2)
 214          if (qnerr.ne.0) then
 215             write(*,*) 'Grid file problem: ',gridname
 216             if (qnerr.lt.0) then
 217                write(*,*) 'Grid file does not exist'
 218                write(*,*) 'Calculating and creating grid file'
 219                call qnfilw(0,0)
 220                open(unit=2,status='unknown',file=gridname,              &
 221      &              form='unformatted')
 222                call QNDUMP(2)
 223                close(2)
 224             else
 225                write(*,*) 'Existing grid file is inconsistent'
 226                if (qnerr.eq.1)                                          &
 227      &              write(*,*) 'Defined grid different'
 228                if (qnerr.eq.2)                                          &
 229      &              write(*,*) 'Heavy quark weight table different'
 230                if (qnerr.eq.3)                                          &
 231      &              write(*,*) 'Charm mass different'
 232                if (qnerr.eq.4)                                          &
 233      &              write(*,*) 'Bottom mass different'
 234                stop
 235             endif
 236          endif
 237       endif
 238
 239 !--   Apply cuts to grid
 240 !--taking away the s cut at 600d0
 241       CALL GRCUTS(-1D0,-1D0,-1D0,-1D0)
 242
 243
 244
 245
 246 !--   Choose renormalisation and factorisation scales
 247
 248                                 ! renormalisation
 249       CALL QNRSET('AAAR2',1D0)
 250       CALL QNRSET('BBBR2',0D0)
 251                                 ! factorisation (light)
 252       CALL QNRSET('AAM2L',1D0)
 253       CALL QNRSET('BBM2L',0D0)
 254                                 ! factorisation (heavy)
 255       CALL QNRSET('AAM2H',1D0)
 256       CALL QNRSET('BBM2H',0D0)
 257
 258 !       ZM=91.187D0
 259       imem=0
 260 !      print *,imem
 261 ! -- only need call to listPDF here to get alphas
 262       call listPDF(nset,imem,xval)
 263 !      print *,xval
 264         AS=XVAL(1)
 265 !        AS=0.118d0
 266       CALL QNRSET('ALFQ0',ZM*ZM)
 267       CALL QNRSET('ALFAS',AS)
 268
 269 !      ZM2=ZM*ZM
 270       ALPHAS=QNALFA(ZM2)
 271 !      WRITE(*,*)'ALPHAS AT Mz2',ALPHAS
 272
 273 !--   Book non-singlet distributions
 274
 275       CALL QNBOOK(2,'UPLUS')
 276       CALL QNBOOK(3,'DPLUS')
 277       CALL QNBOOK(4,'SPLUS')
 278       CALL QNBOOK(5,'CPLUS')
 279       CALL QNBOOK(6,'BPLUS')
 280       CALL QNBOOK(7,'UPVAL')
 281       CALL QNBOOK(8,'DNVAL')
 282
 283 !--   Book linear combinations for proton for f = 3,4,5 flavours
 284
 285 !--define some quark pdfs
 286          CALL dVZERO(PWGT,20)
 287         PWGT(2) = 0.5
 288
 289         PWGT(7) = -0.5
 290
 291         PWGT(1) = 0.5/3.
 292         CALL QNLINC(17,'UB',3,PWGT)
 293         PWGT(1) = 0.5/4.
 294         CALL QNLINC(17,'UB',4,PWGT)
 295         PWGT(1) = 0.5/5.
 296         CALL QNLINC(17,'UB',5,PWGT)
 297         CALL dVZERO(PWGT,20)
 298
 299         PWGT(4) = 0.5
 300         PWGT(1) = 0.5/3.
 301         CALL QNLINC(18,'SB',3,PWGT)
 302         PWGT(1) = 0.5/4.
 303         CALL QNLINC(18,'SB',4,PWGT)
 304         PWGT(1) = 0.5/5.
 305         CALL QNLINC(18,'SB',5,PWGT)
 306          CALL dVZERO(PWGT,20)
 307         CALL QNLINC(19,'CB',3,PWGT)
 308         PWGT(5) = 0.5
 309         PWGT(1) = 0.5/4.
 310
 311         CALL QNLINC(19,'CB',4,PWGT)
 312         PWGT(1) = 0.5/5.
 313         CALL QNLINC(19,'CB',5,PWGT)
 314         CALL dVZERO(PWGT,20)
 315         PWGT(3) = 0.5
 316
 317         PWGT(8) = -0.5
 318
 319         PWGT(1) = 0.5/3.
 320         CALL QNLINC(20,'DB',3,PWGT)
 321         PWGT(1) = 0.5/4.
 322         CALL QNLINC(20,'DB',4,PWGT)
 323         PWGT(1) = 0.5/5.
 324         CALL QNLINC(20,'DB',5,PWGT)
 325          CALL dVZERO(PWGT,20)
 326         CALL QNLINC(21,'BB',3,PWGT)
 327        CALL QNLINC(21,'BB',4,PWGT)
 328         PWGT(6) = 0.5
 329
 330         PWGT(1) = 0.5/5.
 331         CALL QNLINC(21,'BB',5,PWGT)
 332 !---
 333
 334       return
 335 !
 336 !=======================================================================
 337       entry ZEUSpdf(nset)
 338 !      ZM = 91.187d0
 339 !      zm2 = zm*zm
 340
 341 !      ALPHAS=QNALFA(ZM2)
 342
 343
 344 !      imem=mem
 345       call getnmem(nset,imem)
 346       call listPDF(nset,imem,xval)
 347 !      print *,nset,imem
 348 !      print *,xval
 349 !      print *,imem,xval
 350
 351       UA=XVAL(3)
 352       UB=XVAL(4)
 353       UE=0.0d0
 354       UC=XVAL(5)
 355
 356       DA=XVAL(7)
 357       DB=XVAL(8)
 358       DE=0.0d0
 359       DC=XVAL(9)
 360
 361       GA=XVAL(11)
 362       GB=XVAL(12)
 363       GE=0.0d0
 364       GC=XVAL(13)
 365
 366       SN=XVAL(14)
 367       SA=XVAL(15)
 368       SB=XVAL(16)
 369       SE=0.0d0
 370       SC=XVAL(17)
 371
 372       DLN=XVAL(18)
 373       DLA=XVAL(19)
 374 !      DLB=XVAL(20)
 375       DLB=XVAL(16)+2.0d0
 376       DLE=0.0d0
 377       DLC=XVAL(21)
 378       AS=XVAL(1)
 379        CALL QNRSET('ALFAS',AS)
 380 !--   Input quark distributions at Q2 = Q02 GeV2
 381 !-- WRITE IN ACTUAL VALUES TO SAVE USING dgamma
 382 !       UN=2./AREA(UA-1,UB,UE,UC)
 383 !       DN=1./AREA(DA-1,DB,DE,DC)
 384       UN = XVAL(2)
 385       DN=XVAL(6)
 386 !       UBDBN=DLN/AREA(DLA-1,DLB,DLE,DLC)
 387 !      call grgive(nxgri,xxmin,xxmax,nqgri,qqmin,qqmax)
 388       nxgri=nx
 389       DO IX = 1,NXGRI
 390         xX = XFROMIX(IX)
 391         UPVAL=UN*FF_LHA(Xx,UA,UB,UE,UC)
 392
 393
 394         DNVAL=DN*FF_LHA(Xx,DA,DB,DE,DC)
 395
 396
 397         SEA=SN*FF_LHA(Xx,SA,SB,SE,SC)
 398
 399 !        GN=(1-UN*AREA(UA,UB,UE,UC)-
 400 !     .        DN*AREA(DA,DB,DE,DC)-
 401 !     .        SN*AREA(SA,SB,SE,SC))/AREA(GA,GB,GE,GC)
 402         GN=XVAL(10)
 403         GLUON=GN*FF_LHA(Xx,GA,GB,GE,GC)
 404
 405
 406
 407 !        UMD=UBDBN*FF_LHA(X,DLA,DLB,DLE,DLC)
 408         UMD=DLN*FF_LHA(Xx,DLA,DLB,DLE,DLC)
 409         SINGL=UPVAL+DNVAL+SEA
 410
 411 !        print *,un,dn,sn,gn,dln
 412
 413         CSEA=0.0
 414         SSEA=0.2*SEA
 415         USEA=(SEA-SSEA-CSEA)/2.-UMD
 416         DSEA=(SEA-SSEA-CSEA)/2.+UMD
 417         UPLUS=UPVAL+USEA-1./NACT*SINGL
 418         DPLUS=DNVAL+DSEA-1./NACT*SINGL
 419         SPLUS=SSEA-1./NACT*SINGL
 420
 421
 422         CALL QNPSET('GLUON',IX,IQ0,GLUON)
 423         CALL QNPSET('SINGL',IX,IQ0,SINGL)
 424         CALL QNPSET('UPLUS',IX,IQ0,UPLUS)
 425         CALL QNPSET('DPLUS',IX,IQ0,DPLUS)
 426         CALL QNPSET('SPLUS',IX,IQ0,SPLUS)
 427         CALL QNPSET('UPVAL',IX,IQ0,UPVAL)
 428         CALL QNPSET('DNVAL',IX,IQ0,DNVAL)
 429       ENDDO
 430 !--THINGS ARE FINE FOR HEAVY SO DO IT
 431       IF (HEAVY) THEN
 432
 433 !--     Evolve over whole grid
 434
 435         CALL EVOLSG(IQ0,1,NQGRI)
 436         CALL EVPLUS('UPLUS',IQ0,1,NQGRI)
 437         CALL EVPLUS('DPLUS',IQ0,1,NQGRI)
 438         CALL EVPLUS('SPLUS',IQ0,1,NQGRI)
 439         CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
 440         CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
 441
 442       ELSE
 443
 444 !--     Evolve gluon, singlet and valence over whole grid
 445
 446         CALL EVOLSG(IQ0,1,NQGRI)
 447         CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
 448         CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
 449
 450 !--     Be more careful with plus distributions
 451
 452         IF (NACT.EQ.3) THEN
 453 !--THINGS ARE ALSO FINE IF 1Q0 IS BELOW 1QC THEN CLEARLY CSEA=0. IS OK
 454 !--SITUATION CD BE 1<Q0<Q2C<Q2B ETC
 455 !--GO DOWN QO TO 1 UP Q0 TO  Q2C
 456           CALL EVPLUS('UPLUS',IQ0,1,IQC)
 457           CALL EVPLUS('DPLUS',IQ0,1,IQC)
 458           CALL EVPLUS('SPLUS',IQ0,1,IQC)
 459 !--DEAL WITH CHARM THRESH
 460           FACTOR=1./3.-1./4.
 461           CALL QADDSI('UPLUS',IQC,FACTOR)
 462           CALL QADDSI('DPLUS',IQC,FACTOR)
 463           CALL QADDSI('SPLUS',IQC,FACTOR)
 464           CALL QNPNUL('CPLUS')
 465           FACTOR=-1/4.
 466           CALL QADDSI('CPLUS',IQC,FACTOR)
 467 !--GO UP Q2C TO Q2B
 468           CALL EVPLUS('UPLUS',IQC,IQC,IQB)
 469           CALL EVPLUS('DPLUS',IQC,IQC,IQB)
 470           CALL EVPLUS('SPLUS',IQC,IQC,IQB)
 471           CALL EVPLUS('CPLUS',IQC,IQC,IQB)
 472
 473         ELSEIF(NACT.EQ.4)THEN
 474 !--1<1QC<1Q0<1QB<ETC
 475 !--NOW WE NEED A RETHINK OF THE INITIAL CONDITIONS
 476 !--FIRST DEAL WITH CPLUS TRUNING ON AT IQC
 477 !--AND GOING UP TO IQB (MIDDLE AGAIN)
 478           CALL QNPNUL('CPLUS')
 479           FACTOR=-1/4.
 480           CALL QADDSI('CPLUS',IQC,FACTOR)
 481           CALL EVPLUS('CPLUS',IQC,IQC,IQB)
 482 !--REDIFINE THE PLUS DUSTNS TAKING INTO ACCOUNT CPLUS
 483           CALL QNPNUL('SPLUS')
 484           CALL QNPNUL('UPLUS')
 485           CALL QNPNUL('DPLUS')
 486 !--NOW YOU NEED A DO LOOP AGAIN
 487           DO IX = 1,NXGRI
 488           Xx = XFROMIX(IX)
 489           UPVAL=UN*FF_LHA(Xx,UA,UB,UE,UC)
 490
 491           DNVAL=DN*FF_LHA(Xx,DA,DB,DE,DC)
 492           SEA=SN*FF_LHA(Xx,SA,SB,SE,SC)
 493           SINGL=UPVAL+DNVAL+SEA
 494           UMD=DLN*FF_LHA(Xx,DLA,DLB,DLE,DLC)
 495           SSEA=0.2*SEA
 496           CSEA=QPDFIJ('CPLUS',IX,IQ0,IFL) + 1./NACT*SINGL
 497           USEA=(SEA-SSEA-CSEA)/2.-UMD
 498           DSEA=(SEA-SSEA-CSEA)/2.+UMD
 499           UPLUS=UPVAL+USEA-1./NACT*SINGL
 500           DPLUS=DNVAL+DSEA-1./NACT*SINGL
 501           SPLUS=SSEA-1./NACT*SINGL
 502           CALL QNPSET('UPLUS',IX,IQ0,UPLUS)
 503           CALL QNPSET('DPLUS',IX,IQ0,DPLUS)
 504           CALL QNPSET('SPLUS',IX,IQ0,SPLUS)
 505           ENDDO
 506 !-NOW DO MIDDLE FOR UPLUS,DPLUS,SPLUS
 507           CALL EVPLUS('UPLUS',IQ0,IQC,IQB)
 508           CALL EVPLUS('DPLUS',IQ0,IQC,IQB)
 509           CALL EVPLUS('SPLUS',IQ0,IQC,IQB)
 510 !--THEN GO DOWN IQC TO 1
 511              IF(IQC.GT.1)THEN
 512              FACTOR=1./4.-1./3.
 513              CALL QADDSI('UPLUS',IQC,FACTOR)
 514              CALL QADDSI('DPLUS',IQC,FACTOR)
 515              CALL QADDSI('SPLUS',IQC,FACTOR)
 516              CALL EVPLUS('UPLUS',IQC,1,IQC)
 517              CALL EVPLUS('DPLUS',IQC,1,IQC)
 518
 519              CALL EVPLUS('SPLUS',IQC,1,IQC)
 520              ENDIF
 521
 522         ENDIF
 523 !--THEN DEAL WITH B THRESHOLD FOR ALL4
 524         FACTOR=1./4.-1./5.
 525         CALL QADDSI('UPLUS',IQB,FACTOR)
 526         CALL QADDSI('DPLUS',IQB,FACTOR)
 527         CALL QADDSI('SPLUS',IQB,FACTOR)
 528         CALL QADDSI('CPLUS',IQB,FACTOR)
 529 !--THEN GO UP
 530         IF(IQB.LT.NQGRI)THEN
 531        CALL EVPLUS('UPLUS',IQB,IQB,NQGRI)
 532         CALL EVPLUS('DPLUS',IQB,IQB,NQGRI)
 533         CALL EVPLUS('SPLUS',IQB,IQB,NQGRI)
 534         CALL EVPLUS('CPLUS',IQB,IQB,NQGRI)
 535         ENDIF
 536 !C--THEN DEAL WITH B TURNING ON AT IQB AND GO UP
 537         CALL QNPNUL('BPLUS')
 538         FACTOR=-1./5.
 539        CALL QADDSI('BPLUS',IQB,FACTOR)
 540         CALL EVPLUS('BPLUS',IQB,IQB,NQGRI)
 541       ENDIF
 542 !
 543       call getnset(iset)
 544       call save_pdfqcd(iset)
 545
 546       return
 547 !
 548       END
 549 !-----------------------------------------------------------------------
 550 !      DOUBLE PRECISION FUNCTION AREA(A,B,E,C)
 551 !-----------------------------------------------------------------------
 552 !
 553 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 554 !
 555 !      IF (A.LE.-0.99.OR.B.LE.-0.99) THEN
 556 !        AREA=1D6
 557 !        RETURN
 558 !      ENDIF
 559 !
 560 !      AR1=(DGAMMA_LHA(A+1)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2)
 561 !      AR2=E*(DGAMMA_LHA(A+1.5)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2.5)
 562 !      AR3=C*(DGAMMA_LHA(A+2)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+3)
 563 !      AREA=AR1+AR2+AR3
 564 !
 565 !      IF (AREA.LE.1D-6) AREA=1D-6
 566 !
 567 !      RETURN
 568 !      END
 569 !
 570 !-----------------------------------------------------------------------
 571       DOUBLE PRECISION FUNCTION FF_LHA(X,A,B,E,C)
 572 !-----------------------------------------------------------------------
 573
 574       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 575
 576       FF_LHA=X**A*(1D0-X)**B*(1+E*SQRT(X)+C*X)
 577
 578       RETURN
 579       END
 580 !-----------------------------------------------------------------------
 581       SUBROUTINE DVZERO (A,N)
 582 !
 583 ! CERN PROGLIB# F121    VZERO           .VERSION KERNFOR  4.40  940929
 584 ! ORIG. 01/07/71, modif. 24/05/87 to set integer zero
 585 !                 modif. 25/05/94 to depend on QINTZERO
 586 !
 587       implicit real*8 (a-h,o-z)
 588       DIMENSION A(*)
 589       IF (N.LE.0)  RETURN
 590       DO 9 I= 1,N
 591     9 A(I)= 0d0
 592       RETURN
 593       END