*************************************************************************
*
* 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 readB(cpath)
include 'common_c2f.f'
LOGICAL DEBUG
LOGICAL VERBOSE
LOGICAL WARNING
COMMON/DBG/DEBUG,VERBOSE,WARNING
SAVE/DBG/
character*80 cpath
character*80 ppath
b_error=0
la=80
do i=1,80
if(cpath(i:i).eq.'/')la=i
enddo
ppath=cpath(1:la)
b_path = ppath
b_pathlen = la
c FUNZIONAAAAAAAAAAAAAAA!!!!!!!!!!!!!!!!!!
c if(b_debug.eq.1)print*,'Field loaded: ',b_loaded
if(DEBUG)print*,'Field loaded: ',b_loaded
if(b_loaded.eq.0)then
c call the subroutine which reads the maps of the measurements taken
c inside the magnetic cavity
call readBinner(ppath)
if(b_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(b_error.eq.1)return
b_loaded = 1
endif
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'
LOGICAL DEBUG
LOGICAL VERBOSE
LOGICAL WARNING
COMMON/DBG/DEBUG,VERBOSE,WARNING
SAVE/DBG/
character*80 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)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)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
b_error = 1
c if(b_debug.eq.1)
if(DEBUG)
$ 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
b_error = 1
c if(b_debug.eq.1)
if(DEBUG)
$ 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
LOGICAL DEBUG
LOGICAL VERBOSE
LOGICAL WARNING
COMMON/DBG/DEBUG,VERBOSE,WARNING
SAVE/DBG/
character*80 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)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
b_error = 1
c if(b_debug.eq.1)
if(DEBUG)
$ 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
b_error = 1
c if(b_debug.eq.1)
if(DEBUG)
$ 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