************************************************************************* * * 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