!! A second example program using a tabulation. It differs from
!! tabulation_example.f90 in that it shows how to consistently carry
!! out initialisation and evolution from separate sub-programs, as
!! might for example be needed in a PDF fit.
!!
!! It outputs a subset of table 15 of hep-ph/0511119 and this output
!! should be identical to the contents of the file
!! tabulation_example.default_output
!!
!! NB: for the full functionality used in generating the HeraLHC and
!!     Les Houches comparison tables, see ../benchmarks/benchmarks.f90
!!     and carefully read the comments at the beginning. Subtleties
!!     exist in particular wrt the treatment of scales muF/=muR.
!!
!! NB: commented code shows usage with LHAPDF (e.g. for cross-checking 
!!     some public PDF set's evolution) -- to use this part of the code
!!     you will also need to link with LHAPDF
!!

!! This module holds the bits and pieces needed to store the grid,
!! the dglap_holder, the coupling and a table
!! 
module table_module
  use hoppet_v1
  !! if using LHAPDF, rename a couple of hoppet functions which
  !! would otherwise conflict with LHAPDF 
  !use hoppet_v1, EvolvePDF_hoppet => EvolvePDF, InitPDF_hoppet => InitPDF
  implicit none
  !! holds information about the grid (not strictly needed here, but handy)
  type(grid_def),  save :: grid, gdarray(4)
  !! holds the splitting functions
  type(dglap_holder), save :: dh
  !! hold the PDF tabulation in the module table_module
  type(pdf_table)       :: table
  !! hold the coupling
  type(running_coupling), save :: coupling
end module table_module


!----------------------------------------------------------------------
!! Either place the intial condition subroutine in a module and "use" that
!! module, otherwise define an explicit interface for the subroutine.
!!
!! Here we go for the first option
module initial_condition
  use hoppet_v1 ! 
  implicit none
contains 
  !======================================================================
  !! The dummy PDF suggested by Vogt as the initial condition for the 
  !! unpolarized evolution (as used in hep-ph/0511119).
  function unpolarized_dummy_pdf(xvals) result(pdf)
    real(dp), intent(in) :: xvals(:)
    real(dp)             :: pdf(size(xvals),ncompmin:ncompmax)
    real(dp) :: uv(size(xvals)), dv(size(xvals))
    real(dp) :: ubar(size(xvals)), dbar(size(xvals))
    !---------------------
    real(dp), parameter :: N_g = 1.7_dp, N_ls = 0.387975_dp
    real(dp), parameter :: N_uv=5.107200_dp, N_dv = 3.064320_dp
    real(dp), parameter :: N_db = half*N_ls
  
    pdf = zero
    ! clean method for labelling as PDF as being in the human representation
    ! (not actually needed after setting pdf=0
    call LabelPdfAsHuman(pdf)

    !-- remember that these are all xvals*q(xvals)
    uv = N_uv * xvals**0.8_dp * (1-xvals)**3
    dv = N_dv * xvals**0.8_dp * (1-xvals)**4
    dbar = N_db * xvals**(-0.1_dp) * (1-xvals)**6
    ubar = dbar * (1-xvals)

    ! labels iflv_g, etc., come from the hoppet_v1 module, inherited
    ! from the main program
    pdf(:, iflv_g) = N_g * xvals**(-0.1_dp) * (1-xvals)**5
    pdf(:,-iflv_s) = 0.2_dp*(dbar + ubar)
    pdf(:, iflv_s) = pdf(:,-iflv_s)
    pdf(:, iflv_u) = uv + ubar
    pdf(:,-iflv_u) = ubar
    pdf(:, iflv_d) = dv + dbar
    pdf(:,-iflv_d) = dbar
  end function unpolarized_dummy_pdf
end module initial_condition


!----------------------------------------------------------------------
!! the main program that does the work once the evolution and table 
!! have been initialised
Program main
  use table_module
  implicit none
  !! hold results at some x, Q
  real(dp) :: Q, pdf_at_xQ(-6:6)
  real(dp), parameter :: heralhc_xvals(9) = &
       & (/1e-5_dp,1e-4_dp,1e-3_dp,1e-2_dp,0.1_dp,0.3_dp,0.5_dp,0.7_dp,0.9_dp/)
  integer  :: ix, iloop

  ! do the common initialisation (including allocating a table)
  ! (call this just once)
  call tabulation_example_startup

  ! the following might be replaced by the "evolve and calculate
  ! chi^2" loop in a fit program; here we just loop over it once.
  ! NB: if you do plan to run a fit program, then you'd be advised
  !     to use the precalculated evolution (see below)
  do iloop = 1, 1
    ! evolve an initial condition to fill the table
    ! (do this for each initial condition)
    call tabulation_example_evolve_table
  
    ! get the value of the tabulation at some point
    Q = 100.0_dp
    write(6,'(a)')
    write(6,'(a,f8.3,a)') "           Evaluating PDFs at Q = ",Q," GeV"
    write(6,'(a5,2a12,a14,a10,a12)') "x",&
         & "u-ubar","d-dbar","2(ubr+dbr)","c+cbar","gluon"
    do ix = 1, size(heralhc_xvals)
       call EvalPdfTable_xQ(table,heralhc_xvals(ix),Q,pdf_at_xQ)
       write(6,'(es7.1,5es12.4)') heralhc_xvals(ix), &
            &  pdf_at_xQ(2)-pdf_at_xQ(-2), &
            &  pdf_at_xQ(1)-pdf_at_xQ(-1), &
            &  2*(pdf_at_xQ(-1)+pdf_at_xQ(-2)), &
            &  (pdf_at_xQ(-4)+pdf_at_xQ(4)), &
            &  pdf_at_xQ(0)
    end do
  end do

  !! deallocate things where relevant (right at end of program,
  !! not strictly needed)
  call tabulation_example_cleanup

End program main


!----------------------------------------------------------------------
!! do the initialisation
subroutine tabulation_example_startup
  use table_module
  implicit none
  real(dp) :: dy, ymax
  integer  :: order, nloop
  real(dp) :: quark_masses(4:6)

  ! set up parameters for grid
  order = -6
  ymax  = 12.0_dp
  dy    = 0.1_dp

  ! set up the grid itself -- we use 4 nested subgrids
  call InitGridDef(gdarray(4),dy/27.0_dp,0.2_dp, order=order)
  call InitGridDef(gdarray(3),dy/9.0_dp,0.5_dp, order=order)
  call InitGridDef(gdarray(2),dy/3.0_dp,2.0_dp, order=order)
  call InitGridDef(gdarray(1),dy,       ymax  ,order=order)
  call InitGridDef(grid,gdarray(1:4),locked=.true.)

  ! initialise the splitting-function holder
  nloop = 3
  call InitDglapHolder(grid,dh,factscheme=factscheme_MSbar,&
       &                      nloop=nloop,nflo=3,nfhi=6)
  write(6,'(a)') "Splitting functions initialised!"

  !--- if doing many evolutions with different couplings, the
  !    following lines should be placed inside
  !    tabulation_example_evolve_table, and the coupling should be
  !    deleted each time around once you're done with it.
  !
  ! allocate and initialise the running coupling with a given
  ! set of quark masses (NB: charm mass just above Q0).
  quark_masses(4:6) = (/1.414213563_dp, 4.5_dp, 175.0_dp/)
  call InitRunningCoupling(coupling,alfas=0.35_dp,Q=sqrt(2.0_dp),nloop=nloop,&
       &                   quark_masses = quark_masses)

  ! create the tables that will contain our copy of the user's pdf
  ! as well as the convolutions with the pdf.
  call AllocPdfTable(grid, table, Qmin=1.0_dp, Qmax=10000.0_dp, & 
       & dlnlnQ = dy/4.0_dp, freeze_at_Qmin=.true.)
  ! add information about the nf transitions to the table (improves
  ! interpolation quality)
  call AddNfInfoToPdfTable(table,coupling)

  ! "pre-evolve" so that subsequent evolutions are faster,
  ! specifying the initial scale for PDFs (here sqrt(2.0))
  !call PreEvolvePdfTable(table, sqrt(2.0_dp), dh, coupling)

end subroutine tabulation_example_startup

!----------------------------------------------------------------------
!!
subroutine tabulation_example_evolve_table
  use table_module
  use initial_condition
  implicit none
  !! hold the initial pdf
  real(dp), pointer :: pdf0(:,:)
  real(dp) :: Q0_pdf
  integer  :: nloop_ev
  !! define the interfaces for LHA pdf (by default not used)
  !! (NB: unfortunately this conflicts with an internal hoppet name,
  !! so make sure that you "redefine" the internal hoppet name, 
  !! as illustrated in the commented "use" line above:
  !! use hoppet_v1, EvolvePDF_hoppet => EvolvePDF, ...)
  ! interface
  !    subroutine EvolvePDF(x,Q,res)
  !      use types; implicit none
  !      real(dp), intent(in)  :: x,Q
  !      real(dp), intent(out) :: res(*)
  !    end subroutine EvolvePDF
  ! end interface


  ! initialise a PDF from the function below (must be contained,
  ! in a "used" module, or with an explicitly defined interface)
  call AllocPDF(grid, pdf0)

  ! get the PDF from the routine contained in the initial_condition
  ! module
  pdf0 = unpolarized_dummy_pdf(xValues(grid))
  Q0_pdf   = sqrt(2.0_dp)

  !! or get it from LHAPDF
  ! set up LHAPDF with cteq
  !call InitPDFsetByName("cteq61.LHgrid")
  !call InitPDF(0)
  ! allocate and set up the initial pdf from LHAPDF ...
  !Q0_pdf = 5.0_dp 
  !call InitPDF_LHAPDF(grid, pdf0, EvolvePDF, Q0_pdf)

  ! choose the number of loops (must not be higher than what's
  ! been initialised)
  nloop_ev = 3

  ! create the tabulation based on the evolution of pdf0 from scale Q0
  call EvolvePdfTable(table, Q0_pdf, pdf0, dh, coupling, nloop=nloop_ev)
  ! use the following if instead the table has been pre-evolved
  !call EvolvePdfTable(table,pdf0)
  write(6,'(a)') "Evolution done!"

  call Delete(pdf0)
end subroutine tabulation_example_evolve_table


!----------------------------------------------------------------------
!! deallocate things 
subroutine tabulation_example_cleanup
  use table_module
  implicit none
  ! some cleaning up (not strictly speaking needed, but illustrates
  ! how it's done)
  call Delete(table)
  call Delete(dh)
  call Delete(coupling)
  call Delete(grid)

end subroutine tabulation_example_cleanup