!******************************************************************************
!
!.....Calculation of the Q^2 distribution and total cross
!     section of the process
!
!       \nu_\ell + 12C --> \ell^- + X
!
!     in the CCQE, resonance production and DIS channels.
!
!     ilept = 1,2,3 for muon, electron or tau neutrino
!     ig selects the vector form factors (see below)
!
!     Omar Benhar
!     Last revised: Roma, September 2018
!
!******************************************************************************
!
      include 'd2sigma_Ar.f'
!
      implicit none
!
      integer istop,nnn,nmax,nevent,n_gauss
!
      parameter( nnn = 500000 , nmax = 500 , nevent = 50 )
      parameter( n_gauss =  128 )
!
      real*8 event(nevent),x(nnn),y(nnn),z(nnn),omega(n_gauss)
     $      ,T_lept(n_gauss),Q2(nmax),distQ2(nmax),F_int(n_gauss)
     $      ,weight(n_gauss)
!
      real*8 zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept,EA,ER,E_min,eps
     $      ,pi,E_nu,time_0,time_1,time_2,h_Q2,E_lept,p_lept
     $      ,omega_min,omega_max,dsigma,aux,JD,cos_theta_b,small
     $      ,k_F,sigma_tot,h_Q2_1,h_Q2_2,Q2_min,Q2_max,Q2_mid
     $      ,dsigma_QE,dsigma_RES,dsigma_DIS,xmpi,gg,amu
!
      integer nsamples,ig,ilept,nestr,iev,indev,i,j,k
     $       ,ip,np,n_Q2,n_Q2_1,n_Q2_2,n_last,n_omega
!
      common / config1 / nsamples
      common / config2 / x,y,z
!
      common / constants / zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept
     $                    ,EA,ER,E_min,eps,k_F
!
      common / flags / ig,ilept
c
      open (3,file='Q2dist.in')
      open (7,file='evgen_argon_n.tape')
      open (8,file='Q2dist.out')
      open (9,file='sigma_tot.out')
c
      call cpu_time(time_0)
c
c.....set constants [energies are in GeV]
c
      pi = 4.*atan(1.)
      zero = 0.
      small = 1.e-4
      amu = 930.4940954
      xmp = .938272081  ! proton mass [GeV]
      xmn = .939565413  ! neutron mass [GeV]
      xmpi = .13957018  ! pion mass [GeV]
      xm = .5*(xmp + xmn)   ! mass of isoscalar nucleon [GeV]
      xmA = 39.962      ! mass of 40Ar nucleus [GeV]
      xmAm1 = 38.964*amu    ! mass of 39Ar nucleus [GeV]
c     xmAm1 = 38.968*amu    ! mass of 39Cl nucleus [GeV]
      xmR_0 = xmA - xmn
      EA = 18.*xmp + 22.*xmn - xmA  ! carbon binding energy [GeV]
      ER = 18.*xmp + 21.*xmn - xmAm1    ! chlorine binding energy [GeV]
      E_min = EA - ER
      k_F = .251    ! Ar Fermi momentum. TAken from Ankowski & Sobczyk
c
      n_omega = 128
c     n_omega = 128
      omega_min = zero
c
c.....read ig, to select the vector form factors,
c     ilept, to select neutrino flavour and
c     and the width of the delta-function representation [GeV]
c
c     ig = 1  Dipole parametrization
c     ig = 2  Kelly's parametrization [PRC 70, 068282 (2004)]
c     ig = 3  BBBA parametrization [NPB Proc. Suppl., 159, 127 (2006)]
c
c     ilept = 1 muon neutrino
c     ilept = 2 electron neutrino
c     ilept = 3 tau neutrino
c
      read(3,*)ig
      read(3,*)ilept
      read(3,*)eps
c
c.....set mass of charged lepton
c
      if(ilept.eq.1)then
      xm_lept = 105.6583745/1000.               ! muon mass [GeV]
      else if(ilept.eq.2)then
      xm_lept = .5109989461/1000.               ! electron mass [GeV]
      else if(ilept.eq.3)then
      xm_lept = 1776.86/1000.                   ! tau mass [GeV]
      end if
c
      write(6,*)'               '
      if(ilept.eq.1)then
      write(6,*)' *** muon neutrino beam '
      else if(ilept.eq.2)then
      write(6,*)' *** electron neutrino beam '
      else if(ilept.eq.3)then
      write(6,*)' *** tau neutrino beam '
      end if
      write(6,*)'     charged lepton mass = ',xm_lept,' GeV '
      write(6,*)'               '
      if(ig.eq.1)then
      write(6,*)'     dipole fit of vector form factors'
      else if(ig.eq.2)then
      write(6,*)'     Kelly fit of vector form factors '
      else if(ig.eq.3)then
      write(6,*)'     BBBA fit of vector form factors  '
      end if
      write(6,*)'               '
      write(6,*)'     width of delta-function [GeV] = ',eps
      write(6,*)'               '
c
c.....read total number of samples from unit 7
c
      read(7,*)nestr
      write(6,*)'     # of samples ',nestr
      write(6,*)'               '
c
c.....read sampled momentum, energy and polar angle from unit 7
c     [energy and momentum in GeV]
c
      do iev = 1,nestr
      read(7,9000)indev,(event(j),j=1,3)
      x(iev) = event(1)
      y(iev) = event(2)
      z(iev) = event(3)
      end do
c
c.....read number of samples to be used in the calculation of the cross
c     section (cannot exceed nestr) from unit 5
c
      read(3,*)nsamples
c
      if(nsamples.gt.nestr)then
      write(6,*)' >>> Error: nsamples > nestr '
      stop
      end if
c
c.....read neutrino energy [GeV]
c
      read(3,*)E_nu
      write(6,*)' *** neutrino energy     = ',E_nu,' GeV '
      write(6,*)'                                        '
c
      call cpu_time(time_1)
c
      write(6,*)'                                          '
      write(6,*)' *** elapsed time ',time_1-time_0,' sec   '
c
      write(6,*)'                                           '
      write(6,*)' *** start calculation of Q^2-distribution '
      write(6,*)'                                           '
c
c.....set Q2 grid [GeV^2]
c
c.....low-Q^2 region
c
c     electron or muon neutrinos: 0 = Q2_min < Q2 < Q2_mid = 40*0.015 = 0.6 GeV^2
c     tau neutrino: 0 = Q2_min < Q2 < Q2_mid = 100*0.02 = 2.0 GeV^2
c
      if(ilept.eq.1.or.ilept.eq.2)then
         Q2_min = zero
         n_Q2_1 = 40
         h_Q2_1 = 0.015
      else if(ilept.eq.3)then
         Q2_min = zero
         n_Q2_1 = 100
         h_Q2_1 = 0.020
      end if
c
!     write(6,*)E_nu,xm_lept,Q2_min
c
c.....high-Q^2 region: Q2_mid < Q2 < Q2_max
c
      Q2_mid = Q2_min + n_Q2_1*h_Q2_1
c
c.....high Q^2 region Q2_mid < Q2 < Q2_max
c
c.....Q2_max modified on April 25, 2018
c
      Q2_max = 4.*E_nu**2 - xm_lept**2
c
      if(Q2_max.gt.5.)then
      Q2_max = 5.
      end if
c
      n_Q2_2 = 100
      h_Q2_2 = (Q2_max - Q2_mid)/n_Q2_2
      if(h_Q2_2.le.h_Q2_1)then
      h_Q2_2 = h_Q2_1
      end if
      n_Q2 = n_Q2_1 + n_Q2_2
c
      write(6,*)' >>>>> Q2_min h_Q2_1 n_Q2_1 ',Q2_min,h_Q2_1,n_Q2_1
      write(6,*)' >>>>> Q2_mid Q2_max h_Q2_2 n_Q2_2 '
     $                 ,Q2_mid,Q2_max,h_Q2_2,n_Q2_2
c
      do j = 1,n_Q2
      if(j.le.n_Q2_1)then
      Q2(j) = Q2_min + (j - 0.5)*h_Q2_1
      else
      Q2(j) = Q2(n_Q2_1) + (j - n_Q2_1 - 0.5)*h_Q2_2
      end if
      end do
c
c     write(56,6556)(j,Q2(j),j=1,n_Q2)
c6556 format(i5,e15.5)
c     write(6,*)' enter istop...'
c     read(5,*)istop
c     if(istop.gt.0)stop
c
c.....calculation of Q2 distribution [10^-38 cm^2/GeV^2]
c
      do j = 1,n_Q2             ! start Q2 loop
      distQ2(j) = zero
c
c.....set omega grid [GeV]
c
      gg = ( Q2(j) + xm_lept**2 )/2./E_nu
      omega_max = E_nu - gg/2. - xm_lept**2/2./gg
c     write(6,*)' +++++ ',E_nu,Q2(j),omega_max
c     if(j.eq.1)then
c     n_omega = 4
c     end if
      call dsetgauss(omega_min,omega_max,n_omega,omega,weight)
c
c     if(j.eq.1)then
c     write(6,6226)(k,omega(k),k=1,n_omega)
c6226 format(i5,e15.5)
c     end if
c
         do k = 1,n_omega       ! start omega loop
c
         T_lept(k) = E_nu - omega(k) - xm_lept
         E_lept = T_lept(k) + xm_lept
         p_lept = sqrt(E_lept**2 - xm_lept**2)
         JD = 2.*E_nu*p_lept
         cos_theta_b = ( 2.*E_nu*E_lept - Q2(j) - xm_lept**2 )/JD
c
            if(cos_theta_b.ge.-1..and.cos_theta_b.le.1.)then
c           write(18,1818)Q2(j),omega(k)
c    $                   ,dsigma(E_nu,cos_theta_b,omega(k))
c1818       format(3e25.12)
            F_int(k) = 2.*pi*dsigma(E_nu,cos_theta_b,omega(k))/JD
            else
c           write(28,1818)Q2(j),T_lept(k)
c    $           ,dsigma(E_nu,cos_theta_b,omega(k))
            F_int(k) = zero
            end if
c
         distQ2(j) = distQ2(j)
     $             + F_int(k)*weight(k)
c
         end do             ! end omega loop
c
      write(6,*)j,Q2(j),distQ2(j)
      if(j.gt.n_Q2_1.and.j.lt.n_Q2.and.distQ2(j).lt.small)then
      n_last = j-1
      go to 1000
      else
      n_last = j
      endif
c
      end do                ! end Q2 loop
c
 1000 continue
c
c.....print Q2 distribution
c
c     write(6,6767)(Q2(j),distQ2(j),j=1,n_last)
      write(8,6767)(Q2(j),distQ2(j),j=1,n_last)
c
c.....calculation of total cross section [10^-38 cm^2]
c
      sigma_tot = zero
      do i = 1,n_last
c     do i = 2,n_last
      if(i.le.n_Q2_1)then
      h_Q2 = h_Q2_1
      else
      h_Q2 = h_Q2_2
      end if
      sigma_tot = sigma_tot + distQ2(i)*h_Q2
      end do
c
      write(6,*)'                                        '
      write(6,*)' *** E_nu [Gev] sigma_tot [10^-38 cm^2] '
     $          ,E_nu,sigma_tot
      write(6,*)'     sigma_tot/E_nu [10^-38 cm^2/GeV] ',sigma_tot/E_nu
      write(9,7000)E_nu,sigma_tot,sigma_tot/E_nu
c
      call cpu_time(time_2)
c
      write(6,*)'                                          '
      write(6,*)' *** elapsed time ',time_2-time_1,' sec '
c
      stop
c
 9000 format(i8,3e15.5)
 1212 format(i8,2e15.5)
 6767 format(2e16.5)
 7000 format(3e16.5)
c
      end
c
c******************************************************************************
c
c     Returns abscissas and weights for gaussian integration
c     works with n=4,6,8,10,16,20,24,32,64,128
c
      subroutine dsetgauss(a,b,n,x,w)
      implicit none
      integer i,n,nm,kiplus,kiminus
      real*8 a,b,gauss4(4),gauss6(6),gauss8(8),work(256)
      real*8 gauss10(10),gauss16(16),gauss20(20)
      real*8 gauss24(24),gauss32(32),gauss64(64),gauss128(128)
      real*8 x(128),w(128),xuni(128),wuni(128)
      data gauss4/0.339981043584856d0,0.652145154862546d0,
     #0.861136311594053d0,0.347854845137454d0/
      data gauss6/0.238619186083197d0,0.467913934572691,
     #0.661209386466265d0,0.360761573048139,
     #0.932469514203152d0,0.171324492379170d0/
      data gauss8/0.183434642495650d0,0.362683783378362d0,
     #0.525532409916329d0,0.313706645877887d0,
     #0.796666477413627d0,0.222381034453374d0,
     #0.960289856497536d0,0.101228536290376d0/
      data gauss10/0.148874338981631d0,0.295524224714753d0,
     #0.433395394129247d0,0.269266719309996d0,
     #0.679409568299024d0,0.219086362515982d0,
     #0.865063366688985d0,0.149451349150581d0,
     #0.973906528517172d0,0.066671344308688d0/
      data gauss16/0.095012509837637d0,0.189450610455068d0,
     #0.281603550779258d0,0.182603415044923d0,
     #0.458016777657227d0,0.169156519395002d0,
     #0.617876244402643d0,0.149595988816576d0,
     #0.755404408355003d0,0.124628971255533d0,
     #0.865631202387831d0,0.095158511682492d0,
     #0.944575023073232d0,0.062253523938647d0,
     #0.989400934991649d0,0.027152459411754d0/
      data gauss20/0.076526521133497d0,0.152753387130725d0,
     #0.227785851141645d0,0.149172986472603d0,
     #0.373706088715419d0,0.142096109318382d0,
     #0.510867001950827d0,0.131688638449176d0,
     #0.636053680726515d0,0.118194531961518d0,
     #0.746331906460150d0,0.101930119817240d0,
     #0.839116971822218d0,0.083276741576704d0,
     #0.912234428251325d0,0.062672048334109d0,
     #0.963971927277913d0,0.040601429800386d0,
     #0.993128599185094d0,0.017614007139152d0/
      data gauss24/0.064056892862605d0,0.127938195346752d0,
     #0.191118867473616d0,0.125837456346828d0,
     #0.315042679696163d0,0.121670472927803d0,
     #0.433793507626045d0,0.115505668053725d0,
     #0.545421471388839d0,0.107444270115965d0,
     #0.648093651936975d0,0.097618652104113d0,
     #0.740124191578554d0,0.086190161531953d0,
     #0.820001985973902d0,0.073346481411080d0,
     #0.886415527004401d0,0.059298584915436d0,
     #0.938274552002732d0,0.044277438817419d0,
     #0.974728555971309d0,0.028531388628933d0,
     #0.995187219997021d0,0.012341229799987d0/
      data gauss32/0.048307665687738d0,0.096540088514727d0,
     #0.144471961582796d0,0.095638720079274d0,
     #0.239287362252137d0,0.093844399080804d0,
     #0.331868602282127d0,0.091173878695763d0,
     #0.421351276130635d0,0.087652093004403d0,
     #0.506899908932229d0,0.083311924226946d0,
     #0.587715757240762d0,0.078193895787070d0,
     #0.663044266930215d0,0.072345794108848d0,
     #0.732182118740289d0,0.065822222776361d0,
     #0.794483795967942d0,0.058684093478535d0,
     #0.849367613732569d0,0.050998059262376d0,
     #0.896321155766052d0,0.042835898022226d0,
     #0.934906075937739d0,0.034273862913021d0,
     #0.964762255587506d0,0.025392065309262d0,
     #0.985611511545268d0,0.016274394730905d0,
     #0.997263861849481d0,0.007018610009470d0/
      data gauss64/0.02435029266342d0,0.04869095700914d0,
     &0.07299312178780d0,0.04857546744150d0,
     &0.12146281929612d0,0.04834476223477d0,
     &0.16964442042399d0,0.04799938859631d0,
     &0.21742364374001d0,0.04754016571443d0,
     &0.26468716220877d0,0.04696818281529d0,
     &0.31132287199021d0,0.04628479657947d0,
     &0.35722015833767d0,0.04549162792405d0,
     &0.40227015796399d0,0.04459055815801d0,
     &0.44636601725346d0,0.04358372451999d0,
     &0.48940314570705d0,0.04247351510908d0,
     &0.53127946401989d0,0.04126256322055d0,
     &0.57189564620263d0,0.03995374110007d0,
     &0.61115535517239d0,0.03855015317862d0,
     &0.64896547125466d0,0.03705512854024d0,
     &0.68523631305423d0,0.03547221325688d0,
     &0.71988185017161d0,0.03380516183714d0,
     &0.75281990726053d0,0.03205792835485d0,
     &0.78397235894334d0,0.03023465707240d0,
     &0.81326531512280d0,0.02833967261426d0,
     &0.84062929625258d0,0.02637746971505d0,
     &0.86599939815409d0,0.02435270256871d0,
     &0.88931544599511d0,0.02227017380838d0,
     &0.91052213707850d0,0.02013482315353d0,
     &0.92956917213194d0,0.01795171577570d0,
     &0.94641137485840d0,0.01572603047602d0,
     &0.96100879965205d0,0.01346304789672d0,
     &0.97332682778991d0,0.01116813946013d0,
     &0.98333625388463d0,0.00884675982634d0,
     &0.99101337147674d0,0.00650445796881d0,
     &0.99634011677196d0,0.00414703325659d0,
     &0.99930504173577d0,0.00178327963660d0/
      data gauss128/0.01222369896062d0,0.02444618019626d0,
     &0.03666379096873d0,0.02443156909785d0,
     &0.06108196960414d0,0.02440235563385d0,
     &0.08546364050452d0,0.02435855726469d0,
     &0.10979423112764d0,0.02430020016797d0,
     &0.13405919946119d0,0.02422731922281d0,
     &0.15824404271422d0,0.02413995798902d0,
     &0.18233430598534d0,0.02403816868102d0,
     &0.20631559090208d0,0.02392201213669d0,
     &0.23017356422666d0,0.02379155778099d0,
     &0.25389396642269d0,0.02364688358442d0,
     &0.27746262017790d0,0.02348807601650d0,
     &0.30086543887768d0,0.02331522999401d0,
     &0.32408843502441d0,0.02312844882432d0,
     &0.34711772859764d0,0.02292784414359d0,
     &0.36993955534986d0,0.02271353585011d0,
     &0.39254027503327d0,0.02248565203258d0,
     &0.41490637955227d0,0.02224432889359d0,
     &0.43702450103710d0,0.02198971066819d0,
     &0.45888141983355d0,0.02172194953771d0,
     &0.48046407240417d0,0.02144120553878d0,
     &0.50175955913614d0,0.02114764646770d0,
     &0.52275515205118d0,0.02084144778010d0,
     &0.54343830241281d0,0.02052279248617d0,
     &0.56379664822662d0,0.02019187104117d0,
     &0.58381802162876d0,0.01984888123167d0,
     &0.60349045615855d0,0.01949402805730d0,
     &0.62280219391058d0,0.01912752360826d0,
     &0.64174169256231d0,0.01874958693853d0,
     &0.66029763227265d0,0.01836044393492d0,
     &0.67845892244772d0,0.01796032718214d0,
     &0.69621470836951d0,0.01754947582371d0,
     &0.71355437768359d0,0.01712813541906d0,
     &0.73046756674191d0,0.01669655779679d0,
     &0.74694416679706d0,0.01625500090410d0,
     &0.76297433004409d0,0.01580372865266d0,
     &0.77854847550641d0,0.01534301076088d0,
     &0.79365729476219d0,0.01487312259269d0,
     &0.80829175750791d0,0.01439434499294d0,
     &0.82244311695564d0,0.01390696411961d0,
     &0.83610291506091d0,0.01341127127272d0,
     &0.84926298757797d0,0.01290756272029d0,
     &0.86191546893955d0,0.01239613952122d0,
     &0.87405279695803d0,0.01187730734541d0,
     &0.88566771734540d0,0.01135137629108d0,
     &0.89675328804916d0,0.01081866069944d0,
     &0.90730288340176d0,0.01027947896690d0,
     &0.91731019808096d0,0.00973415335479d0,
     &0.92676925087895d0,0.00918300987166d0,
     &0.93567438827792d0,0.00862637779862d0,
     &0.94402028783022d0,0.00806458989049d0,
     &0.95180196134126d0,0.00749798192563d0,
     &0.95901475785370d0,0.00692689256690d0,
     &0.96565436643197d0,0.00635166316171d0,
     &0.97171681874714d0,0.00577263754287d0,
     &0.97719849146391d0,0.00519016183268d0,
     &0.98209610843572d0,0.00460458425670d0,
     &0.98640674272459d0,0.00401625498374d0,
     &0.99012781849173d0,0.00342552604091d0,
     &0.99325711290021d0,0.00283275147146d0,
     &0.99579275853498d0,0.00223828843096d0,
     &0.99773324862551d0,0.00164250301862d0,
     &0.99907745997738d0,0.00104581267832d0,
     &0.99982488794713d0,0.00044938068609d0/
c
      nm=n/2

      if(n.eq.4)then
         do i=1,n
            work(i)=gauss4(i)
         enddo
      else if(n.eq.6)then
         do i=1,n
            work(i)=gauss6(i)
         enddo
      else if(n.eq.8)then
         do i=1,n
            work(i)=gauss8(i)
         enddo
      else if(n.eq.10)then
         do i=1,n
            work(i)=gauss10(i)
         enddo
      else if(n.eq.16)then
         do i=1,n
            work(i)=gauss16(i)
         enddo
      else if(n.eq.20)then
         do i=1,n
            work(i)=gauss20(i)
         enddo
      else if(n.eq.24)then
         do i=1,n
            work(i)=gauss24(i)
         enddo
      else if(n.eq.32)then
         do i=1,n
            work(i)=gauss32(i)
         enddo
      else if(n.eq.64)then
         do i=1,n
            work(i)=gauss64(i)
         enddo
      else if(n.eq.128)then
         do i=1,n
            work(i)=gauss128(i)
         enddo
      else
        write(6,*)
     $  ' >>> Error: wrong number of points for gaussian integration '
         stop
      endif
c
      do 10 i=1,nm
      kiplus=nm+i
      kiminus=nm+1-i
      xuni(kiplus)=work(2*i-1)
      xuni(kiminus)=-work(2*i-1)
      wuni(kiplus)=work(2*i)
      wuni(kiminus)=work(2*i)
 10   continue
c
      do i=1,n
      x(i)=(b-a)/2.d0*xuni(i)+(b+a)/2.
      w(i)=(b-a)/2.d0*wuni(i)
      end do
c
      return
      end
c
c******************************************************************************
c