C
C     $Id: CKM_read_user.F,v 1.15 2001/03/07 18:43:22 vavilov Exp $
C     $Author: vavilov $
C
      SUBROUTINE CKM_read_user
C
C  This is the user routine called after reading an event from
C    a hit file.  The example here fills an ntuple with some event
C    variables.    EJR 12/25/00
C1
      IMPLICIT   NONE
C
C    These are the CKM common blocks from each sub detector
C 
#include "ckm_units.inc"
#include "ckm_files.inc"
#include "ckm_para.inc"
#include "umsdata.inc"
#include "richdata.inc"
#include "keatdata.inc"
#include "vvscom.inc"
#include "strawdata.inc"
#include "fvsdata.inc"
#include "mvsdata.inc"
c
#include "umsanapar.inc"
#include "umsana.inc"
#include "keatana.inc"
#include "dmsana.inc"
#include "richana.inc"

c
      integer i,j,k,nread,ndecay,naccept
      common /mycom/nread,ndecay,naccept

      real pkaon_mc,ppion_mc,pkgen_mc,ppigen_mc
      real pk_mc(3),ppi_mc(3),xk_mc(3),xpi_mc(3),xk0,yk0
      real*8 ppion,pnorm,ppi(3),xpi(3),pkaon,pk(3),xk(3)
      real*8 hitdiff(3),cross(3),abscross
      real pionfvs,muonfvs,xbiv,ybiv,rbiv
      logical first
      Integer ievent,ntrack,ippart(100),itype(100),it(100),ikaon,ipion
      Real lvert(100,3),lpvert(100,4),ubuf(100),xvtx,yvtx,zvtx
      real*8 theta,mm2,doca,theta_mc,mm2_mc,doca_mc
      real*8 mkaon,mpion,pa(0:3),pb(0:3)
      data mkaon/.493679/,mpion/.13957018/

      common /ckmevt/ievent
      common /ckmtrack/ntrack,ippart,itype,it,lvert,lpvert,ubuf

      real     data(84)
      integer idata(84)
      equivalence(idata,data)

      Logical LFirst,ldebug
      Data LFirst,ldebug/.true.,.false./
      Save LFirst,ldebug
c
c     Vertex stuff. DV.
c
      integer    vmaxtrk
      parameter (vmaxtrk=4)
      integer    vtrklist(vmaxtrk)
      data       vtrklist /vmaxtrk*1/
      save       vtrklist

      real       xyz(3), dxyz(vmaxtrk), x2xyz
      integer    ifl
c
c-------------------------------------------------------------
c  Debug printout if needed
c-------------------------------------------------------------
c
      if(mod(nread,100).eq.0)print *,'Event #',nread

      if(ldebug)then
         print *,'------------------'
         print *,'NREAD = ',nread
         print *,' ntrack =',ntrack
         do i=1,min(ntrack,4)
          print *,'   vertex =',(lvert(i,j),j=1,3)
          print *,'     mom  =',(lpvert(i,j),j=1,3)
          print *,'     part =',ippart(i)
         end do
         print *,'    nhit_wire=',nhit_wire
         print *,'    ntracks_krch,nhits_krch,ndigs_krch=',
     &                ntracks_krch,nhits_krch,ndigs_krch
         print *,'    ntracks_prch,nhits_prch,ndigs_prch=',
     &                ntracks_prch,nhits_prch,ndigs_prch
         print *,'    nhit_keat=',nhit_keat
         print *,'    nhit_straw=',nhit_straw
         print *,'    n_vvsdata=',n_vvsdata
         do i=1,n_vvsdata
          print *,'       #phi modules=',nhits_vvs(i)
         end do
         print *,'    nhit_fvs=',nhit_fvs
         print *,'    nhit_mvs=',nhit_mvs
      end if

c-------------------------------------------------------------
c   Get Monte Carlo values of momenta
c-------------------------------------------------------------
c
c  See if pion or muon hits the FVS
c      
      pionfvs   = 0.
      muonfvs = 0.
      do i=1,nhit_fvs
        if(ipart_fvs(i).eq.8) pionfvs = 1.
        if(ipart_fvs(i).eq.5) muonfvs = 1.
      end do
c
c    MC momenta of kaon and pion from track bank
c
      ikaon=0
      ipion=0
      pkaon_mc=0.
      ppion_mc=0.
      pkgen_mc=0.
      ppigen_mc=0.
      ppion=0
      pkaon=0
      
      do i=1,ntrack
        if(ippart(i).eq.11.or.ippart(i).eq.99.or.ippart(i).eq.98)then
           pkgen_mc=sqrt(lpvert(i,1)**2+lpvert(i,2)**2+lpvert(i,3)**2)
           ikaon=i
        end if
        if(ippart(i).eq.8)then
           ipion=i
           ppigen_mc=sqrt(lpvert(i,1)**2+lpvert(i,2)**2+lpvert(i,3)**2)
        end if
      end do
C
C    MC kaon momentum at last KEAT and pion momentum at first DMS
C
      do i=1,3
        pk_mc(i)=0.
        xk_mc(i)=-9999.
        pk(i)=0.
        xk(i)=-9999.
        ppi_mc(i)=0
        xpi_mc(i)=-9999.
        ppi(i)=0.
        xpi(i)=-9999.
      end do

      do i=1,nhit_keat
        if(itra_keat(i).eq.1.and.
     &     (ipart_keat(i).eq.98.or.ipart_keat(i).eq.99
     &      .or.ipart_keat(i).eq.11) )then
          pk_mc(1)=vect_keat(i,4)*vect_keat(i,7)
          pk_mc(2)=vect_keat(i,5)*vect_keat(i,7)
          pk_mc(3)=sqrt(vect_keat(i,7)**2-pk_mc(1)**2-pk_mc(2)**2)
          pkaon_mc=sqrt(pk_mc(1)**2+pk_mc(2)**2+pk_mc(3)**2)
          xk_mc(1)=vect_keat(i,1)
          xk_mc(2)=vect_keat(i,2)
          xk_mc(3)=vect_keat(i,3)
c
          xbiv=xk_mc(1)+vect_keat(i,4)*550.
          ybiv=xk_mc(2)+vect_keat(i,5)*550.
          rbiv=sqrt(xbiv**2+ybiv**2)
        end if
      end do

      Lfirst=.true.
      do i=1,nhit_straw
       if(nhit_straw.gt.0.and.
     &      ipart_straw(i).eq.8.and.Lfirst)then
          ppi_mc(1)=vect_straw(i,4)*vect_straw(i,7)
          ppi_mc(2)=vect_straw(i,5)*vect_straw(i,7)
          ppi_mc(3)=sqrt(vect_straw(i,7)**2-ppi_mc(1)**2-
     &                   ppi_mc(2)**2)
          ppion_mc=sqrt(ppi_mc(1)**2+ppi_mc(2)**2+ppi_mc(3)**2)
          xpi_mc(1)=vect_straw(i,1)
          xpi_mc(2)=vect_straw(i,2)
          xpi_mc(3)=vect_straw(i,3)
          Lfirst=.false.
       end if
      end do
c
c-------------------------------------------------------------
c  Call fitting routines to get reconstructed values
c-------------------------------------------------------------
c
c   UMS reconstruction  
c
      nttrk=0
      call  umsana
c      if(nttrk.gt.0)then
c         pkaon=pmtmtk(4,1)
c         pnorm=sqrt(1+rxsltk(1)**2+rysltk(1)**2)
c         pk(1)=rxsltk(1)*pkaon/pnorm
c         pk(2)=rysltk(1)*pkaon/pnorm
c         pk(3)=pkaon/pnorm
c         xk(1)=rxintk(1)+rxsltk(1)*xk_mc(3)
c         xk(2)=ryintk(1)+rysltk(1)*xk_mc(3)
c         xk(3)=xk_mc(3)
c      end if

      nkeattrk=0
      call keatana
      if(nkeattrk.eq.1)then
         pkaon=keatptrk(1)
         pnorm=pkaon/sqrt(1+keatxusl(1)**2+keatyusl(1)**2)
         pk(1)=pnorm*keatxusl(1)
         pk(2)=pnorm*keatyusl(1)
         pk(3)=pnorm
         xk(1)=keatxint(1)+keatxusl(1)*xk_mc(3)
         xk(2)=keatyint(1)+keatyusl(1)*xk_mc(3)
         xk(3)=xk_mc(3)
      end if
c
c   DMS reconstruction
c
      ndmstrk=0
      call dmsana
      if(ndmstrk.eq.1)then
         ppion=dmsptrk(1)
         pnorm=ppion/sqrt(1+dmsxusl(1)**2+dmsyusl(1)**2)
         ppi(1)=pnorm*dmsxusl(1)
         ppi(2)=pnorm*dmsyusl(1)
         ppi(3)=pnorm
         xpi(1)=dmsxint(1)+dmsxusl(1)*xpi_mc(3)
         xpi(2)=dmsyint(1)+dmsyusl(1)*xpi_mc(3)
         xpi(3)=xpi_mc(3)
      end if
c
c   RICH reconstruction
c
      call richana
c
c-------------------------------------------------------------
c  Now calculate DOCA, THETA from tracking
c-------------------------------------------------------------
c
c   First from MC values
c
      if(pkaon_mc.gt.0.and.ppion_mc.gt.0)then
         hitdiff(1)=dble(xk_mc(1)-xpi_mc(1))
         hitdiff(2)=dble(xk_mc(2)-xpi_mc(2))
         hitdiff(3)=dble(xk_mc(3)-xpi_mc(3))
         cross(1)=ppi_mc(2)*pk_mc(3)-ppi_mc(3)*pk_mc(2)
         cross(2)=ppi_mc(3)*pk_mc(1)-ppi_mc(1)*pk_mc(3)
         cross(3)=ppi_mc(1)*pk_mc(2)-ppi_mc(2)*pk_mc(1)
         abscross=dsqrt(cross(1)**2+cross(2)**2+cross(3)**2)
         doca_mc=(hitdiff(1)*cross(1) + hitdiff(2)*cross(2) +
     &         hitdiff(3)*cross(3))/abscross
         theta_mc = acos((ppi_mc(1)*pk_mc(1)+ppi_mc(2)*pk_mc(2)+
     &              ppi_mc(3)*pk_mc(3))/pkaon_mc/ppion_mc)
      else
         doca_mc=-9999.
         theta_mc=-9999.
      end if
c
c   Then from reconstructed values
c
      if(pkaon.gt.0.and.ppion.gt.0)then
         hitdiff(1)=dble(xk(1)-xpi(1))
         hitdiff(2)=dble(xk(2)-xpi(2))
         hitdiff(3)=dble(xk(3)-xpi(3))
         cross(1)=ppi(2)*pk(3)-ppi(3)*pk(2)
         cross(2)=ppi(3)*pk(1)-ppi(1)*pk(3)
         cross(3)=ppi(1)*pk(2)-ppi(2)*pk(1)
         abscross=dsqrt(cross(1)**2+cross(2)**2+cross(3)**2)
         doca=(hitdiff(1)*cross(1) + hitdiff(2)*cross(2) +
     &         hitdiff(3)*cross(3))/abscross
         theta = acos((ppi(1)*pk(1)+ppi(2)*pk(2)+
     &                 ppi(3)*pk(3))/pkaon/ppion)
         call find_vtx(1, 1, vtrklist, xyz, dxyz, x2xyz, ifl)
         if(ifl.eq.0) then
            xvtx = xyz(1)
            yvtx = xyz(2)
            zvtx = xyz(3)
         else
            xvtx=-9999.
            yvtx=-9999.
            zvtx=-9999.
         endif
      else
         doca=-9999.
         theta=-9999.
         xvtx=-9999.
         yvtx=-9999.
         zvtx=-9999.
      end if
c
c Angle of decay from RICH 
c
      if(prch(2).gt.0)then
        thetarch = acos((prch(2)*rchangx(2)*pk(1) + 
     &                  prch(2)*rchangy(2)*pk(2) + 
     &                  prch(2)*
     &      sqrt(1-rchangx(2)**2-rchangy(2)**2)*pk(3))/pkaon/prch(2))
      else
        thetarch=-9999.
      end if
c
c  Now calculate missing mass square exactly
c
c
c            MC:
c
      pa(0)=dsqrt(pkaon_mc**2 + mkaon**2)
      do i=1,3
        pa(i)=pk_mc(i)
      end do
      pb(0)=dsqrt(ppion_mc**2 + mpion**2)
      do i=1,3
        pb(i)=ppi_mc(i)
      end do
      mm2_mc = (pa(0)-pb(0))**2 
     &     -(pa(1)-pb(1))**2
     &     -(pa(2)-pb(2))**2
     &     -(pa(3)-pb(3))**2
c
c            Tracker:
c
      pa(0)=dsqrt(pkaon**2 + mkaon**2)
      do i=1,3
        pa(i)=pk(i)
      end do
      pb(0)=dsqrt(ppion**2 + mpion**2)
      do i=1,3
        pb(i)=ppi(i)
      end do
      mm2 = (pa(0)-pb(0))**2 
     &     -(pa(1)-pb(1))**2
     &     -(pa(2)-pb(2))**2
     &     -(pa(3)-pb(3))**2
c
c          RICH:
c
      pa(0)=dsqrt(prch(1)**2 + mkaon**2)
      pnorm=prch(1)/dsqrt(1+rchangx(1)**2+rchangy(1)**2)
      pa(1)=pnorm*rchangx(1)
      pa(2)=pnorm*rchangy(1)
      pa(3)=pnorm
      pb(0)=dsqrt(prch(2)**2 + mpion**2)
      pnorm=prch(2)/dsqrt(1+rchangx(2)**2+rchangy(2)**2)
      pb(1)=pnorm*rchangx(2)
      pb(2)=pnorm*rchangy(2)
      pb(3)=pnorm
      mm2rch = (pa(0)-pb(0))**2 
     &        -(pa(1)-pb(1))**2
     &        -(pa(2)-pb(2))**2
     &        -(pa(3)-pb(3))**2


c
c-------------------------------------------------------------
C  Fill histos and Ntuple with information in common blocks 
C     For list of variables, see 
c
c         $CKMGEANT_DIR/include/ *data.inc  and  *ana.inc
c
c-------------------------------------------------------------
C
      call hfill(10,pkaon_mc,0.,1.)
      call hfill(20,ppion_mc,0.,1.)

      i=1
      data(i)=float(ievent)
      i=i+1
      idata(i) = ckm_nrndm(1)
      i = i + 1
      idata(i) = ckm_nrndm(2)
      i = i + 1
      data(i)=float(ntrack)

      i=i+1
      data(i)=ippart(1)
      i=i+1
      data(i)=lvert(1,1)
      i=i+1
      data(i)=lvert(1,2)
      i=i+1
      data(i)=lvert(1,3)
      i=i+1
      data(i)=lpvert(1,1)
      i=i+1
      data(i)=lpvert(1,2)
      i=i+1
      data(i)=lpvert(1,3)

      i=i+1
      data(i) =float(ndigs_krch)
      i=i+1
      data(i) =float(ndigs_prch)
      i=i+1
      data(i) =float(nhit_wire)
      i=i+1
      data(i) =float(nhit_straw)
      i=i+1
      data(i) =float(nhit_keat)
      i=i+1
      data(i)=float(n_vvsdata)
      i=i+1
      data(i)=float(nhit_fvs)
      i=i+1
      data(i)=float(nhit_mvs)
      i=i+1
      data(i)=pionfvs
      i=i+1
      data(i)=muonfvs
      if(ntrack.ge.3)then
        i=i+1
        data(i)=lvert(2,1)
        i=i+1
        data(i)=lvert(2,2)   
        i=i+1
        data(i)=lvert(2,3)
      else
        do j=1,3
         i=i+1
         data(i)=-9999
        end do
      end if
      i=i+1
      data(i)=xvtx
      i=i+1
      data(i)=yvtx
      i=i+1
      data(i)=zvtx
      i=i+1
      data(i)=x2xyz
      i=i+1
      data(i)=pkgen_mc
      i=i+1
      data(i)=pkaon_mc
      i=i+1
      data(i)=ppion_mc
      i=i+1
      data(i)=sngl(pkaon)
      i=i+1
      data(i)=sngl(ppion)
      do j=1,3
        i=i+1
        data(i)=pk_mc(j)
      end do
      do j=1,3
        i=i+1
        data(i)=pk(j)
      end do
      do j=1,3
        i=i+1
        data(i)=ppi_mc(j)
      end do
      do j=1,3
        i=i+1
        data(i)=ppi(j)
      end do
      i=i+1
      data(i)=xk_mc(1)
      i=i+1
      data(i)=xk_mc(2)
      i=i+1
      data(i)=xk(1)
      i=i+1
      data(i)=xk(2)
      i=i+1
      data(i)=rbiv
      i=i+1
      data(i)=xpi_mc(1)
      i=i+1
      data(i)=xpi_mc(2)
      i=i+1
      data(i)=xpi(1)
      i=i+1
      data(i)=xpi(2)
      i=i+1
      data(i)=float(nttrk)
      i=i+1
      data(i)=rchsqtk(1)
      i=i+1
      data(i)=rxintk(1)
      i=i+1
      data(i)=rxsltk(1)
      i=i+1
      data(i)=ryintk(1)
      i=i+1
      data(i)=rysltk(1)
      i=i+1
      data(i)=float(ndmstrk)
      i=i+1
      data(i)=dmschi2(1)
      if(ndmstrk.eq.0)data(i)=0
      i=i+1
      data(i)=dmsxint(1)
      i=i+1
      data(i)=dmsxusl(1)
      i=i+1
      data(i)=dmsyint(1)
      i=i+1
      data(i)=dmsyusl(1)
      i=i+1
      data(i)=sngl(theta_mc)
      i=i+1
      data(i)=sngl(doca_mc)
      i=i+1
      data(i)=sngl(mm2_mc)
      i=i+1
      data(i)=sngl(theta)
      i=i+1
      data(i)=sngl(doca)
      i=i+1
      data(i)=sngl(mm2)
      i=i+1
      data(i)=radius_dig(1)
      i=i+1
      data(i)=radius_dig(2)
      i=i+1
      data(i)=rchchi2(1)
      i=i+1
      data(i)=rchchi2(2)
      i=i+1
      data(i)=sngl(prch(1))
      i=i+1
      data(i)=sngl(prch(2))
      i=i+1
      data(i)=sngl(thetarch)
      i=i+1
      data(i)=sngl(mm2rch)
      i=i+1
      data(i)=rchangx(1)
      i=i+1
      data(i)=rchangy(1)
      i=i+1
      data(i)=rchangx(2)
      i=i+1
      data(i)=rchangy(2)

      call hfn(1,data)

c-------------------------------------------------------------
c  Count how many kaons decay within accepted decay area
c-------------------------------------------------------------
c
      if(zvtx.ge.1400.and.zvtx.le.4400)ndecay=ndecay+1
C
      If(  zvtx.ge.2000.and.zvtx.le.4000.and.
     &     ppion.ge.14.and.ppion.le.20.and.
     &     mm2.gt.0.and.mm2.lt..008.and.theta.gt..0025
     &     ) then
        naccept=naccept+1
      end if

      Return
      End
C
C  $Log: CKM_read_user.F,v $
c Revision 1.15  2001/03/07  18:43:22  vavilov
c Initial particle kinematics as well as vertex "chi2" were added to Ntuple.
c Variable names: 'ipst','xvst','yvst','zvst','pxst','pyst','pzst' and 'chi2vtx'
c
C  Revision 1.14  2001/03/07 17:08:02  kushnir
C  Added random seeds to ntuple
C
c Revision 1.13  2001/03/07  02:26:57  ramberg
c Calculate missing mass exactly instead of using approximation
c
C  Revision 1.12  2001/03/07 00:50:36  vavilov
C  {z,x,y}vtx are filled with a new vertexing code. x2xyz - a-la chi2 for a
C  vertex (actually sum of squares of distances from vertex to tracks)
C
C  Revision 1.11  2001/02/27 21:18:17  ramberg
C  Added umsana to the analysis and save its variables in ntuple.
C
C  Revision 1.10  2001/02/23 16:01:34  jurgen
C  change ifimx to nhits_vvs
C
C  Revision 1.9  2001/02/13 16:59:48  jurgen
C  fix comment (linux complains)
C
C  Revision 1.8  2001/02/12 05:40:26  ramberg
C  Added calls to routines for analysis and fill ntuple appropriately
C
C  Revision 1.7  2001/01/19 21:58:59  ramberg
C  Use single precision for a few quantities when using HFN
C
C  Revision 1.6  2001/01/18 15:19:42  jurgen
C  change hitfvs variable to real. change comment to please linux
C
C  Revision 1.5  2001/01/18 01:19:29  jurgen
C  fix problem with double end-if in debug at end
C
C  Revision 1.4  2001/01/14 17:17:51  ramberg
C  Took out sections that weren't general.  Fixed bug on if loop
C
C  Revision 1.3  2001/01/13 05:40:10  ramberg
C  new example ckm_read user routine
C
C  Revision 1.2  2001/01/05 20:12:21  jurgen
C  add cvs headers
C
C
