subroutine pval
      vector vinp(7),veq(8),veqs(8),vhou(8),vhous(8)
      vector par(1),step(1),pmin(1),pmax(1),errpar(1)
      real xj(8),veqs2(8),veqsi2(8)

      iseed=123456
      seed=0.123456

      if(vinp(1).eq.1) then
        siginv=0.
        vigid=0.
        do j=1,8
          veqs2(j)=veqs(j)**2
          veqsi2(j)=1./veqs2(j)
          siginv=siginv+veqsi2(j)
          vigid=vigid+veq(j)/veqs2(j)
        enddo
        vigid=vigid/siginv
        sigerr=sqrt(1./siginv)
        chi2d=0.
        do j=1,8
          chi2d=chi2d+(veq(j)-vigid)**2/veqs2(j)
        enddo
        print*,' data:',vigid,sigerr,chi2d
        pvald=0
        isim=100000
        rsim=isim
        do i=1,isim
          vigi=0.
          do j=1,8
            call rannor(iseed,xg1)
            xj(j)=vigid+xg1*veqs(j)
            vigi=vigi+xj(j)/veqs2(j)
          enddo
          vigi=vigi/siginv
          call hfill(101,vigi,0.,1.)
          chi2=0.
          do j=1,8
            chi2=chi2+(xj(j)-vigi)**2/veqs2(j)
          enddo
          call hfill(100,chi2,0.,1.)
          if(chi2.ge.chi2d) pvald=pvald+1
        enddo
        pvalue=pvald/rsim
        print*,' p-value',pvalue
      endif

      return
c23456789012345678901234567890123456789012345678901234567890123456789012
      end