source/magnet/magnetic-field.tar

common_B_inner.f0100640000076700007640000000271110252260040014041 0ustar  vannucciwizard*************************************************************************
*     
*     Common common_B_inner.f
*     
*     to be included in:
*     - ../magnet/read_B_inner.f
*     - ../magnet/inter_B.f
*     - ../magnet/inter_B_inner.f
*     
*************************************************************************

c      implicit double precision (a-h,o-z)


      parameter (nx=29, ny=23, nz=101) !number of measures along X, Y and Z axes

c     coordinates in m of the edges of the volume in which the field
c     is interpolated according to the inner maps
      parameter (edgexmin=-0.085,edgexmax=0.085
     $     ,edgeymin=-0.07,edgeymax=0.07
     $     ,edgezmin=-0.26,edgezmax=0.26)

************
c     first map
      real*8 px1(nx,3),py1(ny,3),pz1(nz,3) !coordinates of measure points:
c     e.g. py1(ny,1) = Y coordinates of Bx (=1) component of magnetic field

      real*8 b1(nx,ny,nz,3)     !magnetic field values:
c     e.g. b1(nx,ny,nz,2) = By (=2) component of magnetic field measured in (nx, ny, nz)

      real*8 px1max(3),px1min(3),py1max(3),py1min(3),pz1max(3),pz1min(3) !grid edges

      common/interpolation1/px1,py1,pz1,b1
     $     ,px1max,px1min,py1max,py1min,pz1max,pz1min

************
c     second map
      real*8 px2(nx,3),py2(ny,3),pz2(nz,3)
      real*8 b2(nx,ny,nz,3)
      real*8 px2max(3),px2min(3),py2max(3),py2min(3),pz2max(3),pz2min(3)

      common/interpolation2/px2,py2,pz2,b2
     $     ,px2max,px2min,py2max,py2min,pz2max,pz2min
common_B_outer.f0100640000076700007640000000303410252260012014062 0ustar  vannucciwizard*************************************************************************
*     
*     Common common_B_outer.f
*     
*     to be included in:
*     - ../magnet/read_B_outer.f
*     - ../magnet/inter_B.f
*     - ../magnet/inter_B_outer.f
*     
*************************************************************************

c      implicit double precision (a-h,o-z)


c     number of measures along X, Y and Z axes
      parameter (nox=13, noy=13, noz=4) 

c     coordinates in m of the edges of the volume in which the field
c     is interpolated according to the inner maps
c     UPPER VOLUME
      parameter (edgeuxmin=-0.18)
      parameter (edgeuxmax=0.18)
      parameter (edgeuymin=-0.18)
      parameter (edgeuymax=0.18)
      parameter (edgeuzmin=0.28)
      parameter (edgeuzmax=0.37)
c     LOWER VOLUME
      parameter (edgelxmin=edgeuxmin)
      parameter (edgelxmax=edgeuxmax)
      parameter (edgelymin=edgeuymin)
      parameter (edgelymax=edgeuymax)
      parameter (edgelzmin=-0.37)
      parameter (edgelzmax=-0.28)

************
c     MAGNETIC-FIELD MAP
      real*8 pox(nox,3),poy(noy,3),poz(noz,3) !coordinates of measure points:
c     e.g. py1(ny,1) = Y coordinates of Bx (=1) component of magnetic field

      real*8 bo(nox,noy,noz,3)     !magnetic field values:
c     e.g. b1(nx,ny,nz,2) = By (=2) component of magnetic field measured in (nx, ny, nz)

      real*8 poxmax(3),poxmin(3),poymax(3),poymin(3),pozmax(3),pozmin(3) 
c     grid edges

      common/interpolationo/pox,poy,poz,bo
     $     ,poxmax,poxmin,poymax,poymin,pozmax,pozmin

inter_B.f0100640000076700007640000001021710252303110012473 0ustar  vannucciwizard*************************************************************************
*     
*     Subroutine inter_B.f
*     
*     it computes the magnetic field in a chosen point x,y,z inside or 
*     outside the magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     if the point falls outside the interpolation volume, set the field to 0
*     
*     needs:
*     - ../common/common_B_inner.f common file for the inner magnetic field map
*     - ./inter_B_inner.f common file for the inner magnetic field map
*     - ../common/common_B_outer.f common file for the outer magnetic field map
*     - ./inter_B_outer.f common file for the outer magnetic field map
*     
*     to be called after ./read_B.f (magnetic field map reading subroutine)
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

      subroutine inter_B(x,y,z,res) !coordinates in m, magnetic field in T

      implicit double precision (a-h,o-z)
      include '../common/common_B_inner.f'
      include '../common/common_B_outer.f'


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

      real*8 x,y,z              !point of interest
      real*8 res(3)             !interpolated B components: res = (Bx, By, Bz)

      real*8 zl,zu
      real*8 resu(3),resl(3)

c------------------------------------------------------------------------
c     
c     set the field outside the interpolation volume to be 0
c     
c------------------------------------------------------------------------

      do ip=1,3
        res(ip)=0.
      enddo


c------------------------------------------------------------------------
c     
c     check if the point falls inside the interpolation volumes
c     
c------------------------------------------------------------------------

*     -----------------------
*     INNER MAP
*     -----------------------
      if(
     $          (x.ge.edgexmin).and.(x.le.edgexmax)
     $     .and.(y.ge.edgeymin).and.(y.le.edgeymax)
     $     .and.(z.ge.edgezmin).and.(z.le.edgezmax)
     $     ) then

        call inter_B_inner(x,y,z,res)
c        print*,'INNER - ',z,res

      endif

*     -----------------------
*     OUTER MAP
*     -----------------------
      if(
     $         ((x.ge.edgeuxmin).and.(x.le.edgeuxmax)
     $     .and.(y.ge.edgeuymin).and.(y.le.edgeuymax)
     $     .and.(z.ge.edgeuzmin).and.(z.le.edgeuzmax))
     $     .or.
     $         ((x.ge.edgelxmin).and.(x.le.edgelxmax)
     $     .and.(y.ge.edgelymin).and.(y.le.edgelymax)
     $     .and.(z.ge.edgelzmin).and.(z.le.edgelzmax))
     $     ) then
         
        call inter_B_outer(x,y,z,res)
c        res(2)=res(2)*10
c        print*,'OUTER - ',z,res

      endif

*     --------------------------------
*     GAP between INNER and OUTER MAPS
*     --------------------------------
      if(
     $          (x.gt.edgexmin).and.(x.lt.edgexmax)
     $     .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
     $     .and.(z.gt.edgezmax).and.(z.lt.edgeuzmin)
     $     )then
         
         zu = edgeuzmin
         zl = edgezmax
         call inter_B_inner(x,y,zl,resu)
         call inter_B_outer(x,y,zu,resl)

         do i=1,3
            res(i) = z  * ((resu(i)-resl(i))/(zu-zl))
     $           + resu(i)
     $           -   zu * ((resu(i)-resl(i))/(zu-zl))
         enddo
c        print*,'GAP U - ',z,res

      elseif(
     $          (x.gt.edgexmin).and.(x.lt.edgexmax)
     $     .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
     $     .and.(z.gt.edgelzmax).and.(z.lt.edgezmin)
     $     ) then
         
        
         zu = edgezmin
         zl = edgelzmax
         call inter_B_inner(x,y,zu,resu)
         call inter_B_outer(x,y,zl,resl)
         
         do i=1,3
            res(i) = z  * ((resu(i)-resl(i))/(zu-zl))
     $           + resu(i)
     $           -   zu * ((resu(i)-resl(i))/(zu-zl))
         enddo
c        print*,'GAP D - ',z,res

      endif

      return
      end


      include './inter_B_inner.f'
      include './inter_B_outer.f'
inter_B_inner.f0100640000076700007640000001662510252260077013715 0ustar  vannucciwizard*************************************************************************
*     
*     Subroutine inter_B_inner.f
*     
*     it computes the magnetic field in a chosen point x,y,z inside the
*     magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     the value is computed for two different inner maps and then averaged
*     
*     needs:
*     - ../common/common_B_inner.f
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

      subroutine inter_B_inner(x,y,z,res) !coordinates in m, magnetic field in T

      implicit double precision (a-h,o-z)
      include '../common/common_B_inner.f'


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

      real*8 x,y,z              !point of interpolation
      real*8 res(3)             !interpolated B components: res = (Bx, By, Bz)
      real*8 res1(3),res2(3)    !interpolated B components for the two maps

      integer ic                !index for B components:
                                ! ic=1 ---> Bx
                                ! ic=2 ---> By
                                ! ic=3 ---> Bz

      integer cube(3)           !vector of indexes identifying the cube
                                ! containing the point of interpolation
                                ! (see later...)
      
      real*8 xl,xh,yl,yh,zl,zh  !cube vertexes coordinates

      real*8 xr,yr,zr           !reduced variables (coordinates of the 
                                ! point of interpolation inside the cube)

      real*8 Bp(8)              !vector of values of B component 
                                ! being computed, on the eight cube vertexes


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

      do ic=1,3                 !loops on the three B components

c------------------------------------------------------------------------
c     
c     chooses the coordinates interval containing the input point
c     
c------------------------------------------------------------------------
c     e.g.:
c     
c     x1    x2    x3    x4    x5...
c     |-----|-+---|-----|-----|----
c     ~~~~~~~~x
c     
c     in this case the right interval is identified by indexes 2-3, so the 
c     value assigned to cube variable is 2

        cube(1)=INT((nx-1)*(x-px1min(ic))/(px1max(ic)-px1min(ic))) +1
        cube(2)=INT((ny-1)*(y-py1min(ic))/(py1max(ic)-py1min(ic))) +1
        cube(3)=INT((nz-1)*(z-pz1min(ic))/(pz1max(ic)-pz1min(ic))) +1
        
c------------------------------------------------------------------------
c     
c     if the point falls beyond the extremes of the grid...
c     
c------------------------------------------------------------------------
c     
c     ~~~~~~~~~~x1    x2    x3...
c     - - + - - |-----|-----|----
c     ~~~~x
c     
c     in the case cube = 1
c     
c     
c     ...nx-2  nx-1  nx
c     ----|-----|-----| - - - + - -
c     ~~~~~~~~~~~~~~~~~~~~~~~~x
c     
c     in this case cube = nx-1

        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nx) cube(1) = nx-1
        if (cube(2).ge.ny) cube(2) = ny-1
        if (cube(3).ge.nz) cube(3) = nz-1


c------------------------------------------------------------------------
c     
c     temporary variables definition for field computation
c     
c------------------------------------------------------------------------

        xl = px1(cube(1),ic)    !X coordinates of cube vertexes
        xh = px1(cube(1)+1,ic)

        yl = py1(cube(2),ic)    !Y coordinates of cube vertexes
        yh = py1(cube(2)+1,ic)

        zl = pz1(cube(3),ic)    !Z coordinates of cube vertexes
        zh = pz1(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)   !reduced variables
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = b1(cube(1)  ,cube(2)  ,cube(3)  ,ic) !ic-th component of B
        Bp(2) = b1(cube(1)+1,cube(2)  ,cube(3)  ,ic) ! on the eight cube
        Bp(3) = b1(cube(1)  ,cube(2)+1,cube(3)  ,ic) ! vertexes
        Bp(4) = b1(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = b1(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = b1(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = b1(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = b1(cube(1)+1,cube(2)+1,cube(3)+1,ic)

c------------------------------------------------------------------------
c     
c     computes interpolated ic-th component of B in (x,y,z)
c     
c------------------------------------------------------------------------

        res1(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr


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

c     second map is rotated by 180 degree along the Z axis. so change sign
c     of x and y coordinates and then change sign to Bx and By components
c     to obtain the correct result

      x=-x
      y=-y

      do ic=1,3

        cube(1)=INT((nx-1)*(x-px2min(ic))/(px2max(ic)-px2min(ic))) +1
        cube(2)=INT((ny-1)*(y-py2min(ic))/(py2max(ic)-py2min(ic))) +1
        cube(3)=INT((nz-1)*(z-pz2min(ic))/(pz2max(ic)-pz2min(ic))) +1
        
        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nx) cube(1) = nx-1
        if (cube(2).ge.ny) cube(2) = ny-1
        if (cube(3).ge.nz) cube(3) = nz-1

        xl = px2(cube(1),ic)
        xh = px2(cube(1)+1,ic)

        yl = py2(cube(2),ic)
        yh = py2(cube(2)+1,ic)

        zl = pz2(cube(3),ic)
        zh = pz2(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = b2(cube(1)  ,cube(2)  ,cube(3)  ,ic)
        Bp(2) = b2(cube(1)+1,cube(2)  ,cube(3)  ,ic)
        Bp(3) = b2(cube(1)  ,cube(2)+1,cube(3)  ,ic)
        Bp(4) = b2(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = b2(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = b2(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = b2(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = b2(cube(1)+1,cube(2)+1,cube(3)+1,ic)

        res2(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr

      enddo

c     change Bx and By components sign
      res2(1)=-res2(1)
      res2(2)=-res2(2)

c     change back the x and y coordinate signs
      x=-x
      y=-y


c------------------------------------------------------------------------
c     
c     average the two maps results
c     
c------------------------------------------------------------------------

      do ic=1,3
        res(ic)=(res1(ic)+res2(ic))/2
      enddo


      return
      end
inter_B_outer.f0100640000076700007640000001234110252267445013735 0ustar  vannucciwizard*************************************************************************
*     
*     Subroutine inter_B_outer.f
*     
*     it computes the magnetic field in a chosen point x,y,z OUTSIDE the
*     magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     the value is computed for the outer map
*     
*     needs:
*     - ../common/common_B_outer.f
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

      subroutine inter_B_outer(x,y,z,res) !coordinates in m, magnetic field in T

      implicit double precision (a-h,o-z)
      include '../common/common_B_outer.f'


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

      real*8 x,y,z              !point of interpolation
      real*8 res(3)             !interpolated B components: res = (Bx, By, Bz)
      real*8 zin

      integer ic                
c     !index for B components:
c     ! ic=1 ---> Bx
c     ! ic=2 ---> By
c     ! ic=3 ---> Bz
      
      integer cube(3)           
c     !vector of indexes identifying the cube
c     ! containing the point of interpolation
c     ! (see later...)
      
      real*8 xl,xh,yl,yh,zl,zh  !cube vertexes coordinates

      real*8 xr,yr,zr           
c     !reduced variables (coordinates of the 
c     ! point of interpolation inside the cube)

      real*8 Bp(8)              
c     !vector of values of B component 
c     ! being computed, on the eight cube vertexes


c     LOWER MAP
c     ---> up/down simmetry
      zin=z
      if(zin.le.edgelzmax)z=-z

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

      do ic=1,3                 !loops on the three B components

c------------------------------------------------------------------------
c     
c     chooses the coordinates interval containing the input point
c     
c------------------------------------------------------------------------
c     e.g.:
c     
c     x1    x2    x3    x4    x5...               xN
c     |-----|-+---|-----|-----|---- ... ----|-----|
c     ~~~~~~~~x
c     
c     in this case the right interval is identified by indexes 2-3, so the 
c     value assigned to cube variable is 2

        cube(1)=INT((nox-1)*(x-poxmin(ic))/(poxmax(ic)-poxmin(ic))) +1
        cube(2)=INT((noy-1)*(y-poymin(ic))/(poymax(ic)-poymin(ic))) +1
        cube(3)=INT((noz-1)*(z-pozmin(ic))/(pozmax(ic)-pozmin(ic))) +1
        
c------------------------------------------------------------------------
c     
c     if the point falls beyond the extremes of the grid...
c     
c------------------------------------------------------------------------
c     
c     ~~~~~~~~~~x1    x2    x3...
c     - - + - - |-----|-----|----
c     ~~~~x
c     
c     in the case cube = 1
c     
c     
c     ...nx-2  nx-1  nx
c     ----|-----|-----| - - - + - -
c     ~~~~~~~~~~~~~~~~~~~~~~~~x
c     
c     in this case cube = nx-1

        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nox) cube(1) = nox-1
        if (cube(2).ge.noy) cube(2) = noy-1
        if (cube(3).ge.noz) cube(3) = noz-1


c------------------------------------------------------------------------
c     
c     temporary variables definition for field computation
c     
c------------------------------------------------------------------------

        xl = pox(cube(1),ic)    !X coordinates of cube vertexes
        xh = pox(cube(1)+1,ic)

        yl = poy(cube(2),ic)    !Y coordinates of cube vertexes
        yh = poy(cube(2)+1,ic)

        zl = poz(cube(3),ic)    !Z coordinates of cube vertexes
        zh = poz(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)   !reduced variables
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = bo(cube(1)  ,cube(2)  ,cube(3)  ,ic) !ic-th component of B
        Bp(2) = bo(cube(1)+1,cube(2)  ,cube(3)  ,ic) ! on the eight cube
        Bp(3) = bo(cube(1)  ,cube(2)+1,cube(3)  ,ic) ! vertexes
        Bp(4) = bo(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = bo(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = bo(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = bo(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = bo(cube(1)+1,cube(2)+1,cube(3)+1,ic)

c------------------------------------------------------------------------
c     
c     computes interpolated ic-th component of B in (x,y,z)
c     
c------------------------------------------------------------------------

        res(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr


      enddo     

c     LOWER MAP
c     ---> up/down simmetry 
      if(zin.le.edgelzmax)then
         z=-z                   !back to initial ccoordinate
         res(3)=-res(3)         !invert BZ component
      endif

      return
      end
read_B.f0100640000076700007640000000156310252256652012313 0ustar  vannucciwizard*************************************************************************
*     
*     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)
*     
*************************************************************************

      subroutine read_B

c     call the subroutine which reads the maps of the measurements taken
c     inside the magnetic cavity
      call read_B_inner

c     call the subroutine which reads the maps of the measurements taken
c     outside the magnetic cavity
      call read_B_outer        

      return
      end

      include './read_B_inner.f'
      include './read_B_outer.f' 
read_B_inner.f0100640000076700007640000002564510252261245013507 0ustar  vannucciwizard*************************************************************************
*     
*     Subroutine read_B_inner.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_inner.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 read_B_inner

      implicit double precision (a-h,o-z)
      include '../common/common_B_inner.f'


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

      character*64 Bmap_file    !magnetic field file name
c      character*120 cmd1
c      character*120 cmd2
      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------------------------------------------------------------------------
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$$$      cmd1='cp $TRK_GRND/source/magnet/'
c$$$     $     //Bmap_file(1:LNBLNK(Bmap_file))//' .'         
c$$$      call system(cmd1)         
      
c     opens magnetic field map first file
      print *,'Opening file: ',Bmap_file
      call HROPEN
     $     (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
        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$$$      cmd1='cp $TRK_GRND/source/magnet/'
c$$$     $     //Bmap_file(1:LNBLNK(Bmap_file))//' .'         
c$$$      call system(cmd1)         

c     opens magnetic field map first file
      print *,'Opening file: ',Bmap_file
      call HROPEN
     $     (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

      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$$$      cmd2='rm -f '
c$$$     $     //Bmap_file(1:LNBLNK(Bmap_file))
c$$$      call system(cmd2)


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: '
     $     ,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
read_B_outer.f0100640000076700007640000001510610252261222013514 0ustar  vannucciwizard*************************************************************************
*     
*     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

      implicit double precision (a-h,o-z)
      include '../common/common_B_outer.f'


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------------------------------------------------------------------------
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
      print *,'Opening file: ',Bmap_file
      call HROPEN
     $     (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) = 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------------------------------------------------------------------------

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: '
     $     ,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
