*************************************************************************
*     
*     Subroutine read_B.f
*     
*     it calls the subroutines which read the magnetic field maps for
*     the PAMELA spectrometer
*     
*     needs:
*     - ./read_B_inner.f inner map reading subroutine
*     - ./read_B_outer.f outer map reading subroutine
*     
*     to be called before ./inter_B.f (interpolation subroutine)
*     
*************************************************************************

c$$$      subroutine readB(cpath)
c$$$      include 'common_c2f.f'
c$$$
c$$$      LOGICAL DEBUG
c$$$      LOGICAL VERBOSE
c$$$      LOGICAL WARNING
c$$$      COMMON/DBG/DEBUG,VERBOSE,WARNING
c$$$      SAVE/DBG/
c$$$
c$$$      character*256 cpath
c$$$      character*256 ppath
c$$$
c$$$      
c$$$      b_error=0
c$$$      ppath=''
c$$$
c$$$      la=256
c$$$      do i=1,256
c$$$         if(cpath(i:i).eq.'/')la=i
c$$$      enddo
c$$$      ppath=cpath(1:la)
c$$$
c$$$      b_path    = ppath
c$$$      b_pathlen = la
c$$$
c$$$      if(DEBUG)print*,'Field loaded: ',b_loaded
c$$$      if(b_loaded.eq.0)then
c$$$
c$$$c     call the subroutine which reads the maps of the measurements taken
c$$$c     inside the magnetic cavity      
c$$$         call readBinner(ppath)
c$$$         if(b_error.eq.1)return
c$$$
c$$$c     call the subroutine which reads the maps of the measurements taken
c$$$c     outside the magnetic cavity
c$$$         call readBouter(ppath)   
c$$$         if(b_error.eq.1)return
c$$$
c$$$         b_loaded  = 1
c$$$      endif
c$$$
c$$$      return
c$$$      end
      subroutine readB
      include 'common_c2f.f'

c$$$      LOGICAL DEBUG
c$$$      LOGICAL VERBOSE
c$$$      LOGICAL WARNING
c$$$      COMMON/DBG/DEBUG,VERBOSE,WARNING
c$$$      SAVE/DBG/
      INTEGER DEBUG
      INTEGER VERBOSE
      INTEGER WARNING
      COMMON/dbg/DEBUG,VERBOSE,WARNING
      SAVE/dbg/

      

c$$$      character*256 cpath
      character*256 ppath

      
      c2f_error=0

      ppath=c2f_path(1:c2f_pathlen)


c      VERBOSE=.true.

c     call the subroutine which reads the maps of the measurements taken
c     inside the magnetic cavity      
      call readBinner(ppath)
      if(c2f_error.eq.1)return
      
c     call the subroutine which reads the maps of the measurements taken
c     outside the magnetic cavity
      call readBouter(ppath)   
      if(c2f_error.eq.1)return


      return
      end

************************************************************


*************************************************************************
*     
*     Subroutine readBinner
*     
*     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_B.f common file for the inner magnetic field map
*     - .rz map files in ./ containing coordinates of measured points, Bx, By
*     and Bz components + errors
*     
*     output variables: (see common_B_inner.f)
*     - px#(nx,3)         with #=1,2 for the 2 maps
*     - py#(ny,3)
*     - pz#(nz,3)
*     - b#(nx,ny,nz,3)
*     
*************************************************************************

      subroutine readBinner(path)

c     implicit double precision (a-h,o-z)
      include 'common_B.f'
      include 'common_c2f.f'

c$$$      LOGICAL DEBUG
c$$$      LOGICAL VERBOSE
c$$$      LOGICAL WARNING
c$$$      COMMON/DBG/DEBUG,VERBOSE,WARNING
c$$$      SAVE/DBG/
      INTEGER DEBUG
      INTEGER VERBOSE
      INTEGER WARNING
      COMMON/dbg/DEBUG,VERBOSE,WARNING
      SAVE/dbg/

      character*256 path
      
C
      REAL hmemor(10000000)
      integer Iquest(100)    
      COMMON /pawc/hmemor
      save /pawc/
C
      Common /QUEST/ Iquest
      save /quest/
     

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

      character*64 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
     $     PFY(3),FY,DFY
     $     ,PFZ(3),FZ,DFZ
      INTEGER INDEX(3)          !point index

      COMMON /PAWCR4/ INDEX,PFX,FX,DFX,PFY,FY,DFY,PFZ,FZ,DFZ

   
      CALL HLIMIT(10000000)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     largest RZ file: IQUEST(10) records x LREC words x 4 byte
C     with the following settings: 65000 x 4096 x 4 = 1G
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      
c     permette di ottenere ntuple funzionanti nonostante
c     il messaggio dei 64K di RZOUT
      Iquest(10) = 512000
c------------------------------------------------------------------------
c     
c     *** FIRST MAP ***
c     
c------------------------------------------------------------------------

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

c     print*,' '
c     print*,' '

      Bmap_file='measure_n3_290302.rz'
      
c     opens magnetic field map first file
c      if(b_debug.eq.1)print *                 
      if(VERBOSE.EQ.1)print *                 
     $     ,path(1:LNBLNK(path))//Bmap_file
      call HROPEN
     $     (lun_Bmap_file,'Bmap'
     $     ,path(1:LNBLNK(path))//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
        px1max(ic)=0.
        px1min(ic)=0.
        py1max(ic)=0.
        py1min(ic)=0.
        pz1max(ic)=0.
        pz1min(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
        px1(index(1),1) = pfx(1)
        if(px1(index(1),1).lt.px1min(1)) px1min(1)=px1(index(1),1)
        if(px1(index(1),1).gt.px1max(1)) px1max(1)=px1(index(1),1)
        py1(index(2),1) = pfx(2)
        if(py1(index(2),1).lt.py1min(1)) py1min(1)=py1(index(2),1)
        if(py1(index(2),1).gt.py1max(1)) py1max(1)=py1(index(2),1)
        pz1(index(3),1) = pfx(3)
        if(pz1(index(3),1).lt.pz1min(1)) pz1min(1)=pz1(index(3),1)
        if(pz1(index(3),1).gt.pz1max(1)) pz1max(1)=pz1(index(3),1)

        b1(index(1),index(2),index(3),1) = fx


c     By component
        px1(index(1),2) = pfy(1)
        if(px1(index(1),2).lt.px1min(2)) px1min(2)=px1(index(1),2)
        if(px1(index(1),2).gt.px1max(2)) px1max(2)=px1(index(1),2)
        py1(index(2),2) = pfy(2)
        if(py1(index(2),2).lt.py1min(2)) py1min(2)=py1(index(2),2)
        if(py1(index(2),2).gt.py1max(2)) py1max(2)=py1(index(2),2)
        pz1(index(3),2) = pfy(3)
        if(pz1(index(3),2).lt.pz1min(2)) pz1min(2)=pz1(index(3),2)
        if(pz1(index(3),2).gt.pz1max(2)) pz1max(2)=pz1(index(3),2)

        b1(index(1),index(2),index(3),2) = fy


c     Bz component
        px1(index(1),3) = pfz(1)
        if(px1(index(1),3).lt.px1min(3)) px1min(3)=px1(index(1),3)
        if(px1(index(1),3).gt.px1max(3)) px1max(3)=px1(index(1),3)
        py1(index(2),3) = pfz(2)
        if(py1(index(2),3).lt.py1min(3)) py1min(3)=py1(index(2),3)
        if(py1(index(2),3).gt.py1max(3)) py1max(3)=py1(index(2),3)
        pz1(index(3),3) = pfz(3)
        if(pz1(index(3),3).lt.pz1min(3)) pz1min(3)=pz1(index(3),3)
        if(pz1(index(3),3).gt.pz1max(3)) pz1max(3)=pz1(index(3),3)

        b1(index(1),index(2),index(3),3) = fz

      enddo

c------------------------------------------------------------------------
c     
c     closes files
c     
c------------------------------------------------------------------------

      call HREND('Bmap')
      close(lun_Bmap_file)
c$$$  cmd2='rm -f '
c$$$  $     //Bmap_file(1:LNBLNK(Bmap_file))
c$$$  call system(cmd2)
c$$$  


c------------------------------------------------------------------------
c     
c     *** SECOND MAP ***
c     
c------------------------------------------------------------------------

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

c     print*,' '
c     print*,' '

      Bmap_file='measure_n4_110402_corrected.rz'
c     opens magnetic field map first file
c      if(b_debug.eq.1)print *                   !,'Opening file: '
      if(VERBOSE.EQ.1)print *                   !,'Opening file: '
     $     ,path(1:LNBLNK(path))//Bmap_file
      call HROPEN
     $     (lun_Bmap_file,'Bmap'
     $     ,path(1:LNBLNK(path))//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      print*,'iemax ',iemax

      do ic=1,3                 !grid edges
        px2max(ic)=0.
        px2min(ic)=0.
        py2max(ic)=0.
        py2min(ic)=0.
        pz2max(ic)=0.
        pz2min(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. px(4,2) = X coordinate of the point with index = 4 along X, 
c     in which By (=2) component has been measured
c     e.g. b(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
        px2(index(1),1) = pfx(1)
        if(px2(index(1),1).lt.px2min(1)) px2min(1)=px2(index(1),1)
        if(px2(index(1),1).gt.px2max(1)) px2max(1)=px2(index(1),1)
        py2(index(2),1) = pfx(2)
        if(py2(index(2),1).lt.py2min(1)) py2min(1)=py2(index(2),1)
        if(py2(index(2),1).gt.py2max(1)) py2max(1)=py2(index(2),1)
        pz2(index(3),1) = pfx(3)
        if(pz2(index(3),1).lt.pz2min(1)) pz2min(1)=pz2(index(3),1)
        if(pz2(index(3),1).gt.pz2max(1)) pz2max(1)=pz2(index(3),1)

        b2(index(1),index(2),index(3),1) = fx


c     By component
        px2(index(1),2) = pfy(1)
        if(px2(index(1),2).lt.px2min(2)) px2min(2)=px2(index(1),2)
        if(px2(index(1),2).gt.px2max(2)) px2max(2)=px2(index(1),2)
        py2(index(2),2) = pfy(2)
        if(py2(index(2),2).lt.py2min(2)) py2min(2)=py2(index(2),2)
        if(py2(index(2),2).gt.py2max(2)) py2max(2)=py2(index(2),2)
        pz2(index(3),2) = pfy(3)
        if(pz2(index(3),2).lt.pz2min(2)) pz2min(2)=pz2(index(3),2)
        if(pz2(index(3),2).gt.pz2max(2)) pz2max(2)=pz2(index(3),2)

        b2(index(1),index(2),index(3),2) = fy


c     Bz component
        px2(index(1),3) = pfz(1)
        if(px2(index(1),3).lt.px2min(3)) px2min(3)=px2(index(1),3)
        if(px2(index(1),3).gt.px2max(3)) px2max(3)=px2(index(1),3)
        py2(index(2),3) = pfz(2)
        if(py2(index(2),3).lt.py2min(3)) py2min(3)=py2(index(2),3)
        if(py2(index(2),3).gt.py2max(3)) py2max(3)=py2(index(2),3)
        pz2(index(3),3) = pfz(3)
        if(pz2(index(3),3).lt.pz2min(3)) pz2min(3)=pz2(index(3),3)
        if(pz2(index(3),3).gt.pz2max(3)) pz2max(3)=pz2(index(3),3)

        b2(index(1),index(2),index(3),3) = fz

      enddo

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

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

      goto 9000                 !happy ending

c------------------------------------------------------------------------
c     
c     magnetic field map file opening error
c     
c------------------------------------------------------------------------
      
 21   continue
      
      c2f_error = 1
c      if(b_debug.eq.1)
      if(DEBUG.EQ.1)
     $     print*
     $     ,'read_B_inner: ERROR OPENING MAGNETIC FIELD MAP FILE: '
     $     ,Bmap_file
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     ntuple event reading error
c     
c------------------------------------------------------------------------
      
 22   continue
      
      c2f_error = 1
c      if(b_debug.eq.1.eq.1)
      if(DEBUG.EQ.1)
     $     print*,'read_B_inner: ERROR WHILE READING NTUPLE, at entry
     $     : ',iev
      
      goto 9000                 !the end


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

 9000 continue

      return
      end
*************************************************************************
*     
*     Subroutine readBouter
*     
*     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_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 readBouter(path)

c     implicit double precision (a-h,o-z)
      include 'common_B.f'
      include 'common_c2f.f'
C
c$$$      LOGICAL DEBUG
c$$$      LOGICAL VERBOSE
c$$$      LOGICAL WARNING
c$$$      COMMON/DBG/DEBUG,VERBOSE,WARNING
c$$$      SAVE/DBG/
      INTEGER DEBUG
      INTEGER VERBOSE
      INTEGER WARNING
      COMMON/dbg/DEBUG,VERBOSE,WARNING
      SAVE/dbg/

      character*256 path

      REAL hmemor(10000000)
      integer Iquest(100)    
      COMMON /pawc/hmemor
      save /pawc/
C
      Common /QUEST/ Iquest
      save /quest/
     


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

      character*64 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
     $     PFY(3),FY,DFY
     $     ,PFZ(3),FZ,DFZ
      INTEGER INDEX(3)          !point index

      COMMON /PAWCR4/ INDEX,PFX,FX,DFX,PFY,FY,DFY,PFZ,FZ,DFZ

   
c      print*,'Calling HLIMIT'
      CALL HLIMIT(10000000)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     largest RZ file: IQUEST(10) records x LREC words x 4 byte
C     with the following settings: 65000 x 4096 x 4 = 1G
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      
c     permette di ottenere ntuple funzionanti nonostante
c     il messaggio dei 64K di RZOUT
      Iquest(10) = 512000
c------------------------------------------------------------------------
c     
c     *** FIRST MAP ***
c     
c------------------------------------------------------------------------

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

c     print*,' '
c     print*,' '

      Bmap_file='External_top_map_n4_150402.rz'
      
c     opens magnetic field map first file
c      if(b_debug.eq.1)print *
      if(VERBOSE.EQ.1)print *
     $     ,path(1:LNBLNK(path))//Bmap_file
      call HROPEN
     $     (lun_Bmap_file,'Bmap'
     $     ,path(1:LNBLNK(path))//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) = 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) = 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) = 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) = fx


c     By component
        pox(index(1),2) = 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) = 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) = 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) = fy


c     Bz component
        pox(index(1),3) = 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) = 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) = 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) = fz

      enddo


c------------------------------------------------------------------------
c     
c     closes files
c     
c------------------------------------------------------------------------

      call HREND('Bmap')
      close(lun_Bmap_file)


c------------------------------------------------------------------------
c     
c     no error exit
c     
c------------------------------------------------------------------------
      
      goto 9000                 !happy ending

c------------------------------------------------------------------------
c     
c     magnetic field map file opening error
c     
c------------------------------------------------------------------------
      
 21   continue
      
      c2f_error = 1
c      if(b_debug.eq.1)
      if(DEBUG.EQ.1)
     $     print*
     $     ,'read_B_inner: ERROR OPENING MAGNETIC FIELD MAP FILE: '
     $     ,Bmap_file
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     ntuple evevent reading error
c     
c------------------------------------------------------------------------
      
 22   continue
      
      c2f_error  = 1
c      if(b_debug.eq.1)
      if(DEBUG.EQ.1)
     $     print*
     $     ,'read_B_inner: ERROR WHILE READING NTUPLE, at event
     $     : ',iev
      
      goto 9000                 !the end


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

 9000 continue

      return
      end

