*************************************************************************
*     
*     Subroutine read_B_outer.f
*     
*     it reads from rz files the two magnetic field maps taken inside the
*     spectrometer cavity and fills the variables in common_B_inner.f
*     
*     needs:
*     - ../common/common_B_outer.f common file for the outer magnetic field map
*     - .rz map files in ./ containing coordinates of measured points, Bx, By
*     and Bz components + errors
*     
*     output variables: (see common_B_outer.f)
*     - pxo(nx,3)        
*     - pyo(ny,3)
*     - pzo(nz,3)
*     - bo(nx,ny,nz,3)
*     
*************************************************************************

      subroutine read_B_outer(idata_dir)
      include './common_B_outer.for'
c      implicit double precision (a-h,o-z)
      CHARACTER*250 idata_dir
      integer lastnb
c      integer lung



c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------

      character*264 Bmap_file    !magnetic field file name
      parameter (lun_Bmap_file=66) !magnetic field map file id number
      parameter (ntpl_Bmap=20)  !ntuple identifier

      REAL PFX(3),FX,DFX       !Bx field component coordinates in m, value and error in T
      REAL PFY(3),FY,DFY
      REAL PFZ(3),FZ,DFZ
      INTEGER INDEX(3)          !point index
      COMMON /ntplvars/ INDEX,PFX,FX,DFX,PFY,FY,DFY,PFZ,FZ,DFZ
c      save/ntplvars/
c
      parameter(NWPAWC=200000)
      common/pawc/hmem(NWPAWC)
c      save/pawc/
      common/quest/iquest(100)
c      save/quest/
c
c      call hlimit(nwpawc)

c------------------------------------------------------------------------
c     
c     *** FIRST MAP ***
c     
c------------------------------------------------------------------------

c------------------------------------------------------------------------
c     
c     initialization and map file opening
c     
c------------------------------------------------------------------------

c      print*,' '
c      print*,' '
c      idata_dir='/wizard3/pamela/sw/calib/External_top_map_n4_150402.rz'
c      lastnb = length(idata_dir)
c      lung = lastnb+2
      lastnb = length(idata_dir)-1
c      idata_dir(lastnb:lastnb+29)='/External_top_map_n4_150402.rz'
c      lung = lastnb+29
      Bmap_file = idata_dir(1:lastnb)//'/External_top_map_n4_150402.rz'
c      goto 9000
c      Bmap_file='measure_n3_290302.rz'

c$$$      cmd1='cp $TRK_GRND/source/magnet/'
c$$$     $     //Bmap_file(1:LNBLNK(Bmap_file))//' .'         
c$$$      call system(cmd1)         
      
c     opens magnetic field map first file
c      print *,'Opening file: -',idata_dir(1:lung),'- '
      print *,'Opening file: ',Bmap_file(1:lastnb+30)

c      Bmap_file='External_top_map_n4_150402.rz'
      
c     opens magnetic field map first file
c      print *,'Opening file: ',Bmap_file
c
      call HROPEN
     $     (lun_Bmap_file,'Bmap',Bmap_file,'P',1024,istat) 
c     $     (lun_Bmap_file,'Bmap',idata_dir(1:lung),'P',1024,istat) 
c     $     (lun_Bmap_file,'Bmap',idata_dir(1:lung),' ',1024,istat) 
c     $     (lun_Bmap_file,'Bmap','./bin-aux/'//Bmap_file,'P',1024,istat) 
      if(istat.ne.0) goto 21


      call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory

c      call HPRNTU(ntpl_Bmap)
      call HBNAME(ntpl_Bmap,' ',0,'$CLEAR')
      call HBNAME(ntpl_Bmap,'INDEX',index,'$SET')
      call HBNAME(ntpl_Bmap,'BX',pfx,'$SET')
      call HBNAME(ntpl_Bmap,'BY',pfy,'$SET')
      call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET')


c------------------------------------------------------------------------
c     
c     reads events and fills variables
c     
c------------------------------------------------------------------------

      call HNOENT(ntpl_Bmap,iemax) !number of events

c     initializes measurement grid edges
      do ic=1,3
        poxmax(ic)=0.
        poxmin(ic)=0.
        poymax(ic)=0.
        poymin(ic)=0.
        pozmax(ic)=0.
        pozmin(ic)=0.
      enddo


      do iev=1,iemax            !event loop

        call HGNT(ntpl_Bmap,iev,ierr) !reads event
        if(ierr.ne.0) goto 22

c     the output consists of matrices for coordinates, B components values 
c     and errors:
c     e.g. px1(4,2) = X coordinate of the point with index = 4 along X, 
c     in which By (=2) component has been measured
c     e.g. b1(3,23,4,1) = Bx (=1) component value, measured in the point with
c     indexes = 3,23,4 along X, Y and Z

c     Bx component
        pox(index(1),1) = dble(pfx(1))
        if(pox(index(1),1).lt.poxmin(1)) poxmin(1)=pox(index(1),1)
        if(pox(index(1),1).gt.poxmax(1)) poxmax(1)=pox(index(1),1)
        poy(index(2),1) = dble(pfx(2))
        if(poy(index(2),1).lt.poymin(1)) poymin(1)=poy(index(2),1)
        if(poy(index(2),1).gt.poymax(1)) poymax(1)=poy(index(2),1)
        poz(index(3),1) = dble(pfx(3))
        if(poz(index(3),1).lt.pozmin(1)) pozmin(1)=poz(index(3),1)
        if(poz(index(3),1).gt.pozmax(1)) pozmax(1)=poz(index(3),1)

        bo(index(1),index(2),index(3),1) = dble(fx)


c     By component
        pox(index(1),2) = dble(pfy(1))
        if(pox(index(1),2).lt.poxmin(2)) poxmin(2)=pox(index(1),2)
        if(pox(index(1),2).gt.poxmax(2)) poxmax(2)=pox(index(1),2)
        poy(index(2),2) = dble(pfy(2))
        if(poy(index(2),2).lt.poymin(2)) poymin(2)=poy(index(2),2)
        if(poy(index(2),2).gt.poymax(2)) poymax(2)=poy(index(2),2)
        poz(index(3),2) = dble(pfy(3))
        if(poz(index(3),2).lt.pozmin(2)) pozmin(2)=poz(index(3),2)
        if(poz(index(3),2).gt.pozmax(2)) pozmax(2)=poz(index(3),2)

        bo(index(1),index(2),index(3),2) = dble(fy)


c     Bz component
        pox(index(1),3) = dble(pfz(1))
        if(pox(index(1),3).lt.poxmin(3)) poxmin(3)=pox(index(1),3)
        if(pox(index(1),3).gt.poxmax(3)) poxmax(3)=pox(index(1),3)
        poy(index(2),3) = dble(pfz(2))
        if(poy(index(2),3).lt.poymin(3)) poymin(3)=poy(index(2),3)
        if(poy(index(2),3).gt.poymax(3)) poymax(3)=poy(index(2),3)
        poz(index(3),3) = dble(pfz(3))
        if(poz(index(3),3).lt.pozmin(3)) pozmin(3)=poz(index(3),3)
        if(poz(index(3),3).gt.pozmax(3)) pozmax(3)=poz(index(3),3)

        bo(index(1),index(2),index(3),3) = dble(fz)

      enddo


c------------------------------------------------------------------------
c     
c     closes files
c     
c------------------------------------------------------------------------
      call HCDIR('//Bmap',' ')
      call HREND('Bmap')
      close(lun_Bmap_file)


c------------------------------------------------------------------------
c     
c     no error exit
c     
c------------------------------------------------------------------------

c      print*,' '
c      print*,'MAGNETIC FIELD SUCCESSFULLY READ'
c      print*,' '
c      print*,' '
     
      goto 9000                 !happy ending

c------------------------------------------------------------------------
c     
c     magnetic field map file opening error
c     
c------------------------------------------------------------------------
      
 21   continue
      
      print*,' '
      print*,'read_B_inner: ERROR OPENING MAGNETIC FIELD MAP FILE: '
c     $     ,idata_dir(1:lung)
     $     ,Bmap_file
      print*,' '
      print*,' '
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     ntuple event reading error
c     
c------------------------------------------------------------------------
      
 22   continue
      
      print*,' '
      print*,'read_B_inner: ERROR WHILE READING NTUPLE, AT EVENT
     $     : ',iev
      print*,' '
      print*,' '
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     exit
c     
c------------------------------------------------------------------------

 9000 continue

      return
      end
