SUBROUTINE GPRBI ************************************************************************* * * Subroutine gprbi.f (DERIVED FROM read_B_inner routine from * tracker software analysis) * CALLED BY gpdat.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: * - .rz map files in containing coordinates of measured points, Bx, By * and Bz components + errors * * output variables: (see ....) * - px#(nx,3) with #=1,2 for the 2 maps * - py#(ny,3) * - pz#(nz,3) * - b#(nx,ny,nz,3) * ************************************************************************* IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "gpfield.inc" c------------------------------------------------------------------------ c c local variables c c------------------------------------------------------------------------ parameter (ntpl_Bmap=20) !ntuple identifier REAL PFX(3),FFX,DFX, !Bx field component coordinates in m, value and error in T $ PFY(3),FFY,DFY $ ,PFZ(3),FFZ,DFZ INTEGER INDEX(3) !point index COMMON /PAWCR4/ INDEX,PFX,FFX,DFX,PFY,FFY,DFY,PFZ,FFZ,DFZ c------------------------------------------------------------------------ c c *** FIRST MAP *** c c------------------------------------------------------------------------ c------------------------------------------------------------------------ c c initialization and map file opening c c------------------------------------------------------------------------ CALL HCDIR('//FIELD1',' ') call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory 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) = DBLE(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) = DBLE(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) = DBLE(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) = DBLE(ffx) c By component px1(index(1),2) = DBLE(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) = DBLE(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) = DBLE(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) = DBLE(ffy) c Bz component px1(index(1),3) = DBLE(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) = DBLE(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) = DBLE(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) = DBLE(ffz) enddo c------------------------------------------------------------------------ c c closes files c c------------------------------------------------------------------------ call HREND('FIELD1') c------------------------------------------------------------------------ c c *** SECOND MAP *** c c------------------------------------------------------------------------ CALL HCDIR('//FIELD2',' ') call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory 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) = DBLE(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) = DBLE(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) = DBLE(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) = DBLE(ffx) c By component px2(index(1),2) = DBLE(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) = DBLE(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) = DBLE(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) = DBLE(ffy) c Bz component px2(index(1),3) = DBLE(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) = DBLE(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) = DBLE(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) = DBLE(ffz) enddo c------------------------------------------------------------------------ c c closes files c c------------------------------------------------------------------------ call HREND('FIELD2') c------------------------------------------------------------------------ c c no error exit c c------------------------------------------------------------------------ print*,' ' print*,'MAGNETIC FIELD SUCCESSFULLY READ' print*,' ' print*,' ' goto 9000 !happy ending c------------------------------------------------------------------------ c c magnetic field map file opening error c c------------------------------------------------------------------------ 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