LHAPDF/lhapdf5.3.1/wrapzeus.f

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