*$ CREATE GPLEVBIN.FOR *COPY GPLEVBIN * *=== Gplevbin =========================================================* * PROGRAM GPLEVBIN * *----------------------------------------------------------------------* * * * Copyright (C) 1990-2022 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * GnuPlot LEVel BINning setup: * * This code is a derivative of Pawlevbin, which in turn was a * * derivative going back to 1990 of a code preparing the binning * * data to be plotted with the level plotting routine of the * * Mizer graphics package * * * * Authors: Alfredo Ferrari & Paola Sala * * * * * * Created on Sometimes 1990 by Alfredo Ferrari & Paola Sala * * INFN - Milan * * * * Last change on 05-Feb-22 by Alfredo Ferrari * * Private * * * *----------------------------------------------------------------------* * PARAMETER ( NM = 50000 ) PARAMETER ( NMMIN = 5000 ) PARAMETER ( MXZTA = 3000 ) PARAMETER ( MXRAD = 3000 ) PARAMETER ( MXXXXX = MXZTA * MXRAD ) PARAMETER ( EPS = 0.25 ) PARAMETER ( EPS2 = 0.0625) CHARACTER YES*1,FILE1*250,FILE3*250,FILE4*250, & CARD*80,TITPLO*80,FILCMP*80 * CHARACTER CHELP*20 DIMENSION FLUMAT(MXXXXX), FLEMAT(MXXXXX) DIMENSION XX(NM),YY(NM),XXX(NM),YYY(NM) DOUBLE PRECISION XPOINT (NM), YPOINT (NM), ZPOINT (NM), WHAT (6) DOUBLE PRECISION WX0, WY0, WZ0, WX1, WY1, WZ1, WXX, WXY, WXZ, & WYX, WYY, WYZ, WXAXLN, WYAXLN COMMON /LGCURL/NUM,IVET(MXXXXX) LOGICAL LCMPFL, LCMPOU, LINSID, LINPRE, LSMPOI, LGEOPL LOGICAL LRHOFL, LROTA , LAXSWP, LGAXSW LOGICAL LRPHI DOUBLE PRECISION X, Y, XM, DDX, DXBIN, DYBIN CHARACTER SDUM*8 * CALL CMSPPR CALL RTMTXI NUM = MXXXXX FACTX = 1.E+00 FACTY = 1.E+00 GGXM = 2047.0 GGYM = 2047.0 WX0 = 0.D+00 WY0 = 0.D+00 WZ0 = 0.D+00 WXX = 0.D+00 WXY = 0.D+00 WXZ = 0.D+00 WYX = 0.D+00 WYY = 0.D+00 WYZ = 0.D+00 LSMPOI=.FALSE. LGEOPL=.FALSE. LCMPFL=.FALSE. LCMPOU=.FALSE. LROTA =.FALSE. WRITE(*,*)' PLOTGEOM file (def. PLOTGEOM.STORE):' READ(*,'(A50)')FILE1 IF ( FILE1 .EQ. ' ' ) FILE1 = 'PLOTGEOM.STORE' WRITE (*,'(A)') ' Swap plotgeom/bin axis (def=n)?' READ (*,'(A)') YES LAXSWP = YES .EQ. 'Y' .OR. YES .EQ. 'y' LGAXSW = LAXSWP IF ( LAXSWP ) THEN WRITE (*,'(A)') &' Swap: both geo and bin axis (def=0), geo only (1) bin only (2)?' READ (*,'(A)') YES IFSWAP = 0 IF ( YES .NE. ' ' ) READ (YES,*) IFSWAP IF ( IFSWAP .LT. 0 .OR. IFSWAP .GT. 2 ) IFSWAP = 0 LAXSWP = MOD (IFSWAP,2) .EQ. 0 LGAXSW = IFSWAP .LE. 1 END IF * OPEN ( UNIT=99, FILE=FILE1,FORM='UNFORMATTED', & STATUS='OLD', ERR=133 ) * LGEOPL=.TRUE. * New for (geometry) rotations: WRITE (*,*)' Rotate geometry ? (y,n), def n' READ (*,'(A)')YES IF ( YES .EQ. 'y' .OR. YES .EQ. 'Y' ) THEN LROTA =.TRUE. WRITE(*,'(A)')' Enter the 6 rot-defi what''s' READ (*,*) (WHAT(I),I=1,6) SDUM=' ' CALL SETROT (WHAT,SDUM) CALL PRNROT KROTAT = 1 ELSE KROTAT = 0 LROTA=.FALSE. END IF * This part is not reuired when using Flair, still required for Paw: * WRITE(*,*)' Type the "X" axis expansion factor (def=1):' * READ (*,'(A)') CHELP * IF ( CHELP .NE. ' ' ) READ (CHELP,*) FACTX * WRITE(*,*) FACTX * WRITE(*,*)' Type the "Y" axis expansion factor (def=1):' * READ (*,'(A)') CHELP * IF ( CHELP .NE. ' ' ) READ (CHELP,*) FACTY * WRITE(*,*) FACTY OPEN (UNIT=99,FILE=FILE1,FORM='UNFORMATTED', 1 STATUS='OLD') READ(99)TITPLO IF ( TITPLO .EQ. ' ***COMPRESSED***COMPRESSED*** ' ) THEN WRITE (*,*) ' !!! Compressed input file !!! ' LCMPFL = .TRUE. LCMPOU = .FALSE. READ (99) TITPLO ELSE LCMPFL = .FALSE. WRITE (*,*)' (Possible) compressed output file (def='' ''=no)' READ (*,'(A)')FILCMP IF ( FILCMP .NE. ' ' ) THEN LCMPOU = .TRUE. * Start_VAX_seq * OPEN ( UNIT = 98, FILE=FILCMP, STATUS='NEW', * & FORM ='UNFORMATTED' ) * End_VAX_seq * Start_UNIX_seq OPEN ( UNIT = 98, FILE=FILCMP, STATUS='UNKNOWN', & FORM ='UNFORMATTED' ) * End_UNIX_seq WRITE(98)' ***COMPRESSED***COMPRESSED*** ' WRITE(98) TITPLO ELSE LCMPOU=.FALSE. END IF END IF WRITE(*,'(A80)')TITPLO READ(99)X0,Y0,Z0,X1,Y1,Z1,XU,YU,ZU,XV,YV,ZV,XAXLEN,YAXLEN IF (LCMPOU) &WRITE(98)X0,Y0,Z0,X1,Y1,Z1,XU,YU,ZU,XV,YV,ZV,XAXLEN,YAXLEN * X0,Y0,Z0 is the left-bottom corner, X1,Y1,Z1 is the right-top one * Xu,Yu,Zu is the bottom-top axis versor (Y axis just to make * confusion...), the length is Yaxlen * Xv,Yv,Zv is the left-right axis versor (X axis), the length is Xaxlen * Origins of the R (X) and T (Y) axis: T00=X0*XU+Y0*YU+Z0*ZU R00=X0*XV+Y0*YV+Z0*ZV EPSY = 0.99999E+00 * YAXLEN / 2000.E+00 EPSX = 0.99999E+00 * XAXLEN / 2000.E+00 DXEPS = 0.13E+00 * EPSX DYEPS = 0.13E+00 * EPSY *+++++ Debug: WRITE(*,*)' R00,T00',R00,T00 *----- WX0=X0 WY0=Y0 WZ0=Z0 WX1=X1 WY1=Y1 WZ1=Z1 WYX=XU WYY=YU WYZ=ZU WXX=XV WXY=YV WXZ=ZV WXAXLN =XAXLEN WYAXLN =YAXLEN * X0,Y0,Z0 is now the left-top corner: X0 =WX0+WYAXLN*WYX Y0 =WY0+WYAXLN*WYY Z0 =WZ0+WYAXLN*WYZ * X1,Y1,Z1 is now the right-bottom corner: X1 =WX0+WXAXLN*WXX Y1 =WY0+WXAXLN*WXY Z1 =WZ0+WXAXLN*WXZ * Xu,Yu,Zu is now the top-bottom axis versor: XU =-XU YU =-YU ZU =-ZU IF ( LROTA ) THEN XPOINT (1) = X0 YPOINT (1) = Y0 ZPOINT (1) = Z0 XPOINT (2) = X1 YPOINT (2) = Y1 ZPOINT (2) = Z1 NPOINT = 2 CALL DOTRSF ( NPOINT, XPOINT, YPOINT, ZPOINT, KROTAT ) X0 = XPOINT (1) Y0 = YPOINT (1) Z0 = ZPOINT (1) X1 = XPOINT (2) Y1 = YPOINT (2) Z1 = ZPOINT (2) XPOINT (1) = XU YPOINT (1) = YU ZPOINT (1) = ZU XPOINT (2) = XV YPOINT (2) = YV ZPOINT (2) = ZV NPOINT = 2 CALL DORTNO ( NPOINT, XPOINT, YPOINT, ZPOINT, KROTAT ) XU = XPOINT (1) YU = YPOINT (1) ZU = ZPOINT (1) XV = XPOINT (2) YV = YPOINT (2) ZV = ZPOINT (2) END IF *+++++ Debug: WRITE(*,*)X0,Y0,Z0 WRITE(*,*)X1,Y1,Z1 WRITE(*,*)XU,YU,ZU WRITE(*,*)XV,YV,ZV *----- * End reshaffling of input positions/directions: T0=-(X0*XU+Y0*YU+Z0*ZU) T1=-(X1*XU+Y1*YU+Z1*ZU) R0=X0*XV+Y0*YV+Z0*ZV R1=X1*XV+Y1*YV+Z1*ZV XUMIN = MIN (R0,R1) XUMAX = MAX (R0,R1) YUMIN = MIN (T0,T1) YUMAX = MAX (T0,T1) IF ( LGAXSW ) THEN WRITE (*,*) YUMIN, YUMAX WRITE (*,*) ' X_min - X_max (def. provided)' READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) YUMIN,YUMAX WRITE (*,*) XUMIN, XUMAX WRITE (*,*) ' Y_min - Y_max (def. provided)' READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) XUMIN,XUMAX ELSE WRITE (*,*) XUMIN, XUMAX WRITE (*,*) ' X_min - X_max (def. provided)' READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) XUMIN,XUMAX WRITE (*,*) YUMIN, YUMAX WRITE (*,*) ' Y_min - Y_max (def. provided)' READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) YUMIN,YUMAX END IF DX = XUMAX-XUMIN DY = YUMAX-YUMIN *++++++++++ XWMIN=XUMIN-0.05E+00*DX YWMIN=YUMIN-0.05E+00*DY DX = XUMAX - XWMIN DY = YUMAX - YWMIN *---------- G0XM=200. G0YM=200. GMM=1947. IF (FACTX*DX .GE. FACTY*DY) THEN GXM=GMM GYM=DY*FACTY/DX/FACTX*(GXM-G0XM)+G0XM ELSE GYM=GMM GXM=DX*FACTX/DY/FACTY*(GYM-G0YM)+G0YM END IF G1XM = G0XM + (XUMIN-XWMIN) * (GXM-G0XM) / (XUMAX-XWMIN) G1YM = G0YM + (YUMIN-YWMIN) * (GYM-G0YM) / (YUMAX-YWMIN) IF ( LAXSWP .NEQV. LGAXSW ) THEN ZTALOW = YUMIN ZTAHGH = YUMAX RADLOW = XUMIN RADHGH = XUMAX ELSE ZTALOW = XUMIN ZTAHGH = XUMAX RADLOW = YUMIN RADHGH = YUMAX END IF ZMIN = 0.E+00 ZMAX = 0.E+00 133 CONTINUE DO 50 ID=1,100 IPEN=1 WRITE(*,*)' Bin file:' READ(*,'(A)')FILE3 IF (FILE3 .EQ. ' ' .AND. ID .GT. 1 ) GO TO 60 IF (FILE3 .EQ. ' ') GO TO 55 WRITE(*,*) & ' Density file (name) or density value (pos.) or norm. (neg.):' READ(*,'(A)')FILE4 OPEN (UNIT=11,FILE=FILE3,STATUS='OLD',FORM='UNFORMATTED' * Start_VAX_seq * & ,READONLY * End_VAX_seq & ) IF (FILE4 .EQ. ' ') THEN LRHOFL = .FALSE. RHOFIX = 1.E+00 ELSE RHOFIX =-1.E+00 LRHOFL = .TRUE. OPEN (UNIT=12,FILE=FILE4,STATUS='OLD',FORM='UNFORMATTED', * Start_VAX_seq * & READONLY, * End_VAX_seq & ERR=53) GO TO 54 53 CONTINUE READ(FILE4,*)RHOFIX IF ( RHOFIX .LT. 0.E+00 ) THEN LRHOFL = .FALSE. RHOFIX =-RHOFIX END IF 54 CONTINUE END IF RHOTH=1.E-30 WRITE(*,*)' Threshold density (def. provided):' READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) RHOTH IF (RHOTH .LE. 0.E+00) RHOTH=5.E-3 WRITE (*,*) ' Rho_threshold: ', RHOTH WRITE (*,*) ' Which binning?' READ (*,*) NBG * NBG=1 CALL UBDOS2(FLUMAT,FLEMAT,NZTABN,NRADBN,ZTALOW,ZTAHGH,LGEOPL, & RADLOW,RADHGH,RHOTH,MXXXXX,NBG,LRHOFL,RHOFIX,LRPHI) CLOSE (UNIT=11) IF (LRHOFL) CLOSE (UNIT=12) * Start_VAX_seq * OPEN ( UNIT=93, FILE='gplevh.dat', STATUS ='NEW' ) * End_VAX_seq * Start_UNIX_seq OPEN ( UNIT=93, FILE='gplevh.dat', STATUS ='UNKNOWN' ) * End_UNIX_seq * DXBIN = (DBLE(ZTAHGH)-DBLE(ZTALOW))/DBLE(NZTABN) ZMIN = +1.E+30 ZMAX = -1.E+30 * +-------------------------------------------------------------------* * | 2D histogram IF (NRADBN .GT. 1) THEN DYBIN = (DBLE(RADHGH)-DBLE(RADLOW))/DBLE(NRADBN) * If LRPHI find multiples of lines to dump for simulating a R-Phi plot IF ( LRPHI .AND. NZTABN.LT.32) THEN MULT = INT(32.D+00 / DBLE(NZTABN)) ! +/- mod/2 from the 32 DDX = DXBIN / DBLE(MULT) ELSE MULT = 1 DDX = 0.D+00 ENDIF DO IE=1,NZTABN X=DBLE(ZTALOW)+(DBLE(IE)-1.0D+00) * DXBIN DO IM=1,MULT IF ( IM .EQ. 1 ) THEN XM = X ELSE XM = X + DBLE(IM-1)*DDX END IF DO IA=1,NRADBN Y = DBLE(RADLOW)+(DBLE(IA)-1.0D+00)*DYBIN IDOS = (IA-1)*NZTABN + IE IF ( LAXSWP ) THEN WRITE (93,'(1P,2E18.10,2E13.5)') & Y, XM, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 ELSE WRITE (93,'(1P,2E18.10,2E13.5)') & XM, Y, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 END IF IF(ABS(FLUMAT(IDOS)).LT.ZMIN.AND.ABS(FLUMAT(IDOS)) & .GT.2.E-38) ZMIN=ABS(FLUMAT(IDOS)) IF(ABS(FLUMAT(IDOS)).GT.ZMAX) ZMAX=ABS(FLUMAT(IDOS)) END DO IDOS = (NRADBN-1)*NZTABN + IE IF ( LAXSWP ) THEN WRITE (93,'(1P,2E18.10,2E13.5)') & Y+DYBIN, XM, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 ELSE WRITE (93,'(1P,2E18.10,2E13.5)') & XM, Y+DYBIN, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 END IF WRITE (93,*) END DO END DO X = X + DXBIN DO IA=1,NRADBN Y = DBLE(RADLOW)+(DBLE(IA)-1.0D+00)*DYBIN IDOS = (IA-1)*NZTABN + NZTABN IF ( LAXSWP ) THEN WRITE (93,'(1P,2E18.10,2E13.5)') & Y, X, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 ELSE WRITE (93,'(1P,2E18.10,2E13.5)') & X, Y, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 END IF END DO IDOS = (NRADBN-1)*NZTABN + NZTABN IF ( LAXSWP ) THEN WRITE (93,'(1P,2E18.10,2E13.5)') & Y+DYBIN, X, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 ELSE WRITE (93,'(1P,2E18.10,2E13.5)') & X, Y+DYBIN, FLUMAT(IDOS), FLEMAT(IDOS)*100.0 END IF * | * +-------------------------------------------------------------------* * | 1D histogram ELSE WRITE (93,'(A)') '# Detector n: 1D Projection' DO IE=1,NZTABN X = DBLE(ZTALOW)+(DBLE(IE)-1.0D+00)*DXBIN WRITE (93,'(1P,2E18.10,2E13.5)') & X, X+DXBIN, FLUMAT(IE), FLEMAT(IE)*100.0 IF(ABS(FLUMAT(IE)).LT.ZMIN.AND.ABS(FLUMAT(IE)) & .GT.2.E-38) ZMIN=ABS(FLUMAT(IE)) IF(ABS(FLUMAT(IE)).GT.ZMAX) ZMAX=ABS(FLUMAT(IE)) END DO WRITE (93,'(1P,2E18.10,2E13.5)') & X+DXBIN, X+DXBIN, FLUMAT(NZTABN), FLEMAT(NZTABN)*100.0 END IF * | * +-------------------------------------------------------------------* WRITE (*,*) ZMIN, ZMAX 55 CONTINUE CLOSE (UNIT = 93) IF (ID .EQ. 1) THEN * Start_VAX_seq * OPEN(UNIT=89,FILE='gplevh.npo',STATUS='NEW') * OPEN(UNIT=90,FILE='gplevh.poi',STATUS='NEW') * OPEN(UNIT=91,FILE='gplevh.lim',STATUS='NEW') * End_VAX_seq * Start_UNIX_seq OPEN(UNIT=89,FILE='gplevh.npo',STATUS='UNKNOWN') OPEN(UNIT=90,FILE='gplevh.poi',STATUS='UNKNOWN') OPEN(UNIT=91,FILE='gplevh.lim',STATUS='UNKNOWN') * End_UNIX_seq IF ( .NOT. LGEOPL ) THEN IF ( LAXSWP ) THEN YUMIN = ZTALOW XUMIN = RADLOW YUMAX = ZTAHGH XUMAX = RADHGH ELSE XUMIN = ZTALOW YUMIN = RADLOW XUMAX = ZTAHGH YUMAX = RADHGH END IF END IF IF ( LGAXSW ) THEN WRITE(91,*)YUMIN,XUMIN WRITE(91,*)YUMAX,XUMAX ELSE WRITE(91,*)XUMIN,YUMIN WRITE(91,*)XUMAX,YUMAX END IF WRITE(91,*)ZTALOW,RADLOW WRITE(91,*)ZTAHGH,RADHGH WRITE(91,*)ZMIN,ZMAX CLOSE(91) LSMPOI=.FALSE. * Plot now the geometry: IF ( LGEOPL ) THEN DO 3000 II=1,100 * Read the number of worms on tape: READ (99,ERR=3600,END=3600) NWORMS, NDUMMY IF ( NWORMS.EQ.-1 .AND. NDUMMY.EQ.-1 ) GO TO 3500 IF ( LCMPOU ) WRITE(98) NWORMS, NDUMMY DO 2000 IW = 1, NWORMS READ (99) IWORM,IDUMMY,NPOINT IF ( IWORM .NE. IW ) STOP 'IWORM-IW' READ (99) (XX(J),YY(J),J=1,NPOINT) IF (.NOT.LCMPFL) CALL SHRINK (XX,YY,NPOINT,DXEPS,DYEPS) IF ( LCMPOU ) THEN WRITE(98)IWORM,IDUMMY,NPOINT WRITE(98)(XX(J),YY(J),J=1,NPOINT) END IF DO 1500 IP=1,NPOINT XX (IP) = XX (IP) + R00 YY (IP) = YY (IP) + T00 1500 CONTINUE IF ( LROTA ) THEN DO J=1,NPOINT XPOINT (J)= WX0 + ( XX (J) - R00 ) * WXX & + ( YY (J) - T00 ) * WYX YPOINT (J)= WY0 + ( XX (J) - R00 ) * WXY & + ( YY (J) - T00 ) * WYY ZPOINT (J)= WZ0 + ( XX (J) - R00 ) * WXZ & + ( YY (J) - T00 ) * WYZ END DO CALL DOTRSF ( NPOINT, XPOINT, YPOINT, ZPOINT, KROTAT ) DO J=1,NPOINT XX (J) = XPOINT (J) * WXX + YPOINT (J) * WXY & + ZPOINT (J) * WXZ YY (J) = XPOINT (J) * WYX + YPOINT (J) * WYY & + ZPOINT (J) * WYZ END DO END IF IF ( XUMIN .LT. 0.E+00 ) THEN XZMIN = 1.000002E+00 * XUMIN ELSE XZMIN = 0.999998E+00 * XUMIN END IF IF ( YUMIN .LT. 0.E+00 ) THEN YZMIN = 1.000002E+00 * YUMIN ELSE YZMIN = 0.999998E+00 * YUMIN END IF IF ( XUMAX .GT. 0.E+00 ) THEN XZMAX = 1.000002E+00 * XUMAX ELSE XZMAX = 0.999998E+00 * XUMAX END IF IF ( YUMAX .GT. 0.E+00 ) THEN YZMAX = 1.000002E+00 * YUMAX ELSE YZMAX = 0.999998E+00 * YUMAX END IF NPPPP = 0 XXP = XX (1) YYP = YY (1) LINSID=XXP.GE.XZMIN.AND.XXP.LE.XZMAX.AND.YYP.GE.YZMIN & .AND.YYP.LE.YZMAX IF ( LINSID ) THEN NPPPP = NPPPP + 1 XXX (NPPPP) = XXP YYY (NPPPP) = YYP END IF DO 1600 IP=2,NPOINT XXO = XXP YYO = YYP LINPRE=LINSID XXP = XX(IP) YYP = YY(IP) LINSID=XXP.GE.XZMIN.AND.XXP.LE.XZMAX.AND.YYP.GE.YZMIN & .AND.YYP.LE.YZMAX IF ( LINPRE .AND. LINSID ) THEN NPPPP = NPPPP + 1 XXX (NPPPP) = XXP YYY (NPPPP) = YYP ELSE IF ( .NOT. LINPRE .AND. LINSID ) THEN DDXX = XXP - XXO DDX0 = 0.E+00 IF ( XXO .LT. XZMIN ) THEN DDX0 = XZMIN - XXO ELSE IF ( XXO .GT. XZMAX ) THEN DDX0 = XZMAX - XXO END IF DDYY = YYP - YYO DDY0 = 0.E+00 IF ( YYO .LT. YZMIN ) THEN DDY0 = YZMIN - YYO ELSE IF ( YYO .GT. YZMAX ) THEN DDY0 = YZMAX - YYO END IF IF ( ABS (DDXX) .GT. 1.E-10 .AND. ABS (DDYY) & .GT. 1.E-10 ) THEN DDYX = DDYY / DDXX * DDX0 ELSE IF ( ABS (DDXX) .GT. 1.E-10 ) THEN DDYX = 1.E+10 ELSE DDYX = 0.E+00 END IF IF ( ABS (DDYX) .GT. ABS (DDY0) ) THEN XXA = XXO + DDX0 YYA = YYO + DDYY / DDXX * DDX0 ELSE YYA = YYO + DDY0 XXA = XXO + DDXX / DDYY * DDY0 END IF NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA NPPPP = NPPPP + 1 XXX (NPPPP) = XXP YYY (NPPPP) = YYP ELSE IF ( LINPRE .AND. .NOT. LINSID ) THEN DDXX = XXO - XXP DDX1 = 0.E+00 IF ( XXP .LT. XZMIN ) THEN DDX1 = XZMIN - XXP ELSE IF ( XXP .GT. XZMAX ) THEN DDX1 = XZMAX - XXP END IF DDYY = YYO - YYP DDY1 = 0.E+00 IF ( YYP .LT. YZMIN ) THEN DDY1 = YZMIN - YYP ELSE IF ( YYP .GT. YZMAX ) THEN DDY1 = YZMAX - YYP END IF IF ( ABS (DDXX) .GT. 1.E-10 .AND. ABS (DDYY) & .GT. 1.E-10 ) THEN DDYX = DDYY / DDXX * DDX1 ELSE IF ( ABS (DDXX) .GT. 1.E-10 ) THEN DDYX = 1.E+10 ELSE DDYX = 0.E+00 END IF IF ( ABS (DDYX) .GT. ABS (DDY1) ) THEN XXA = XXP + DDX1 YYA = YYP + DDYY / DDXX * DDX1 ELSE YYA = YYP + DDY1 XXA = XXP + DDXX / DDYY * DDY1 END IF NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA ELSE IF ( .NOT. LINPRE .AND. .NOT. LINSID ) THEN DDXX = XXO - XXP DDX1 = 0.E+00 IF ( XXP .LT. XZMIN .AND. XXO .LT. XZMIN ) THEN ELSE IF ( XXP .GT. XZMAX .AND. XXO .GT. XZMAX )THEN ELSE IF ( YYP .LT. YZMIN .AND. YYO .LT. YZMIN )THEN ELSE IF ( YYP .GT. YZMAX .AND. YYO .GT. YZMAX )THEN ELSE DDYY = YYO - YYP IF ( ABS (DDXX) .LT. 1.E-07 * ABS (XZMAX-XZMIN) & ) THEN XXA = XXO YYA = MAX ( YZMIN, MIN ( YYP, YYO ) ) NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA XXA = XXP YYA = MIN ( YZMAX, MAX ( YYP, YYO ) ) NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA ELSE IF ( ABS (DDYY) .LT. 1.E-07 & * ABS (YZMAX-YZMIN) ) THEN YYA = YYO XXA = MAX ( XZMIN, MIN ( XXP, XXO ) ) NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA YYA = YYP XXA = MIN ( XZMAX, MAX ( XXP, XXO ) ) NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA ELSE IPPP = 0 EMME = DDYY / DDXX XXA = MAX ( XZMIN, MIN ( XXP, XXO ) ) YYA = EMME * ( XXA - XXP ) + YYP IF ( YYA .GE. YZMIN .AND. YYA .LE. YZMAX ) & THEN IPPP = IPPP + 1 NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA END IF XXA = MIN ( XZMAX, MAX ( XXP, XXO ) ) YYA = EMME * ( XXA - XXP ) + YYP IF ( YYA .GE. YZMIN .AND. YYA .LE. YZMAX ) & THEN IPPP = IPPP + 1 NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA END IF EMME = DDXX / DDYY YYA = MIN ( YZMAX, MAX ( YYP, YYO ) ) XXA = EMME * ( YYA - YYP ) + XXP IF ( XXA .GE. XZMIN .AND. XXA .LE. XZMAX & .AND. IPPP .LT. 2 ) THEN IPPP = IPPP + 1 NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA END IF YYA = MAX ( YZMIN, MIN ( YYP, YYO ) ) XXA = EMME * ( YYA - YYP ) + XXP IF ( XXA .GE. XZMIN .AND. XXA .LE. XZMAX & .AND. IPPP .LT. 2 ) THEN IPPP = IPPP + 1 NPPPP = NPPPP + 1 XXX (NPPPP) = XXA YYY (NPPPP) = YYA END IF END IF END IF END IF 1600 CONTINUE IF ( NPPPP .GT. 0 ) WRITE(89,*)NPPPP LSMPOI = LSMPOI .OR. NPPPP .GT. 0 IF ( LGAXSW ) THEN DO J=1,NPPPP WRITE(90,*)YYY(J),XXX(J) END DO ELSE DO J=1,NPPPP WRITE(90,*)XXX(J),YYY(J) END DO END IF WRITE(90,*) 2000 CONTINUE 3000 CONTINUE 3500 CONTINUE WRITE(*,*) 'Materials and magnetic field information found.' * Start_UNIX_seq OPEN ( UNIT = 92, FILE = 'plotgeom.his', STATUS='OLD', ERR=3550 ) CLOSE ( UNIT = 92, STATUS = 'DELETE' ) 3550 CONTINUE * End_UNIX_seq * Start_VAX_seq * OPEN ( UNIT=93, FILE='plotgeom.dat', STATUS = 'NEW' ) * End_VAX_seq * Start_UNIX_seq OPEN ( UNIT=93, FILE='plotgeom.dat', STATUS ='UNKNOWN' ) * End_UNIX_seq READ(99,ERR=3600,END=3600) IMXPFX, IMXPFY IF ( LGAXSW ) THEN XMIN=YUMIN XMAX=YUMAX YMIN=XUMIN YMAX=XUMAX ELSE XMIN=XUMIN XMAX=XUMAX YMIN=YUMIN YMAX=YUMAX END IF DXSTEP = (XMAX-XMIN)/FLOAT(IMXPFX) DYSTEP = (YMAX-YMIN)/FLOAT(IMXPFY) DO IY=1,IMXPFY Y=YMIN+(FLOAT(IY)-0.5)*DYSTEP DO IX=1,IMXPFX X=XMIN+(FLOAT(IX)-0.5)*DXSTEP READ(99) NREG,MAT,LAT,BX,BY,B IF ( LGAXSW ) THEN XN = Y YN = X ELSE XN = X YN = Y END IF WRITE (93,*) XN,YN,NREG,MAT,LAT,BX,BY,B END DO WRITE (93,*) END DO IC = 0 3600 CONTINUE END IF IF ( .NOT. LSMPOI ) THEN WRITE (89,*) 0 WRITE (90,*) 0.E+00, 0.E+00 END IF CLOSE (UNIT=99) IF ( LCMPOU ) CLOSE (UNIT=98) END IF IF (FILE3 .EQ. ' ') GO TO 60 * End plot of the geometry 50 CONTINUE 60 CONTINUE 9999 CONTINUE CLOSE(89) CLOSE(90) END *$ CREATE SHRINK.FOR *COPY SHRINK * *=== Shrink ===========================================================* * SUBROUTINE SHRINK (XX,YY,NPOINT,DXEPS,DYEPS) * *----------------------------------------------------------------------* * * * Copyright (C) 1990-2022 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * SHRINK user binnings : * * * * Authors: Alfredo Ferrari & Paola Sala * * * * * * Created on sometimes 1990 by Alfredo Ferrari & Paola Sala * * INFN - Milan * * * * Last change on 05-Feb-22 by Alfredo Ferrari * * Private * * * *----------------------------------------------------------------------* * DIMENSION XX(*), YY(*) IF ( NPOINT .LE. 3 ) RETURN XSCALE = DXEPS YSCALE = DYEPS DUEPS = DXEPS / XSCALE DVEPS = DYEPS / YSCALE DDEPS = SQRT ( DUEPS**2 + DYEPS**2 ) IPOINT = 1 IPP = 1 ULAST = XX (1) / XSCALE VLAST = YY (1) / YSCALE DO 1000 IP = 2, NPOINT - 1 UU = XX (IP) / XSCALE VV = YY (IP) / YSCALE DU = ABS ( UU - ULAST ) DV = ABS ( VV - VLAST ) IF ( ( DU .GE. DUEPS .OR. DV .GE. DVEPS ) .AND. IPOINT .EQ. 1 ) & THEN IPOINT = IPOINT + 1 XX (IPOINT) = XX (IP) YY (IPOINT) = YY (IP) ULSM1 = ULAST VLSM1 = VLAST ULAST = UU VLAST = VV IPP = IP IF ( ABS (ULAST-ULSM1) .LE. 1.0E-02*DUEPS ) THEN EMME = 1.E+31 ELSE IF ( ABS (VLAST-VLSM1) .LE. 1.0E-02*DVEPS ) THEN EMME = 0.E+00 ELSE EMME = ( VLAST - VLSM1 ) / ( ULAST - ULSM1 ) QU = VLSM1 - EMME * ULSM1 END IF * The line is : y = m x + q ELSE IF ( DU .GE. DUEPS .OR. DV .GE. DVEPS ) THEN * x=const line: IF ( ABS (EMME) .GT. 1.E+30 ) THEN DIST = ABS (UU-ULAST) ELSE IF ( ABS (EMME) .LT. 1.E-30 ) THEN DIST = ABS (VV-VLAST) ELSE DIST = ABS ( QU + EMME * UU - VV ) & / SQRT ( 1.E+00 + EMME**2 ) END IF IF ( DIST .GE. DDEPS ) THEN IF ( IPP .LT. IP - 1 ) THEN IPOINT = IPOINT + 1 XX (IPOINT) = XX (IP-1) YY (IPOINT) = YY (IP-1) ULSM1 = ULAST VLSM1 = VLAST ULAST = XX (IPOINT) / XSCALE VLAST = YY (IPOINT) / YSCALE IPP = IP - 1 END IF IPOINT = IPOINT + 1 XX (IPOINT) = XX (IP) YY (IPOINT) = YY (IP) ULSM1 = ULAST VLSM1 = VLAST ULAST = UU VLAST = VV IPP = IP IF ( ABS (ULAST-ULSM1) .LE. 1.0E-02*DUEPS ) THEN EMME = 1.E+31 ELSE IF ( ABS (VLAST-VLSM1) .LE. 1.0E-02*DVEPS ) THEN EMME = 0.E+00 ELSE EMME = ( VLAST - VLSM1 ) / ( ULAST - ULSM1 ) QU = VLSM1 - EMME * ULSM1 END IF ELSE IPP = IP ULAST = UU VLAST = VV XX (IPOINT) = XX (IP) YY (IPOINT) = YY (IP) * Do not recompute m,q: * IF ( ABS (ULAST-ULSM1) .LE. 1.0E-02*DUEPS ) THEN * EMME = 1.E+31 * ELSE IF ( ABS (VLAST-VLSM1) .LE. 1.0E-02*DVEPS ) THEN * EMME = 0.E+00 * ELSE * EMME = ( VLAST - VLSM1 ) / ( ULAST - ULSM1 ) * QU = VLSM1 - EMME * ULSM1 * END IF END IF END IF 1000 CONTINUE IPOINT = IPOINT + 1 XX (IPOINT) = XX (NPOINT) YY (IPOINT) = YY (NPOINT) NPOINT = IPOINT RETURN *=== End of subroutine Shrink =========================================* END *$ CREATE UBDOS2.FOR *COPY UBDOS2 * *=== Ubdos2 ===========================================================* * SUBROUTINE UBDOS2 (FLUMAT,FLEMAT,NZTABN,NRADBN,SZTLOW,SZTHGH, & LGEOPL,SRDLOW,SRDHGH,RHOTH,MXXXXX,NBG,LRHOFL, & RHOFIX, LRPHI) * *----------------------------------------------------------------------* * * * Copyright (C) 1990-2022 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * User Binning DOSe 2nd version : * * * * Author: Alfredo Ferrari & Paola Sala * * * * * * Created on sometimes 1990 by Alfredo Ferrari & Paola Sala * * INFN - Milan * * * * Last change on 05-Feb-22 by Alfredo Ferrari * * Private * * * *----------------------------------------------------------------------* * CHARACTER RUNTIT*80,RUNTIM*32,CHSTAT*10 * *=== usrbin ==========================================================* * *---------------------------------------------------------------------* * Module USRBIN: * * A. Ferrari & A. Fasso': user defined binnings * * Last change A. Ferrari 09-apr-1994 * * * * * * Up to MXUSBN user defined binnings are allowed * * lusbin = logical flag, .true. if at least 1 user defined * * binning is used * * lusevt = logical flag, .true. if at least 1 user defined * * binning event by event is used * * lustkb = logical flag, .true. if at least 1 user defined * * track-length binning is used. Track length bin- * * ning are recognized as binnings where accurate * * deposition along the track is requested but * * for generalized particles other than 208/211 * * nusrbn = number of user defined binnings used * * itusbn = type of binning: 0 = cartesian, .ne. 0 = RZ * * idusbn = distribution to be scored: the usual values are * * allowed. * * titusb = binning name * * ipusbn = logical unit to print the results on: formatted * * if > 0, unformatted if < 0 * * kbusbn = initial location in blank common of the consi- * * dered binning (real*8 address) * * nxbin = number of x (r for RZ) intervals * * nybin = number of y (1 for RZ) intervals * * nzbin = number of z intervals * * xlow/high = minimum and maximum x (r for RZ) * * ylow/high = minimum and maximum y (meaningless for RZ) * * zlow/high = minimum and maximum z * * dxusbn = x (r) bin width * * dyusbn = y bin width * * dzusbn = z bin width * * tcusbn = time cut-off (seconds) for this binning * * bkusbn = 1st Birk's law parameter for this binning * * (meaningful only for energy scoring) * * b2usbn = 2nd Birk's law parameter for this binning * * (meaningful only for energy scoring) * * levtbn = logical flag for binning to be printed at the * * end of each event only * * lntzer = logical flag for printing only non zero cells * * ltrkbn = logical flag for flagging track-length binnings * * * *---------------------------------------------------------------------* * PARAMETER ( MXUSBN = 100 ) LOGICAL LUSBIN, LEVTBN, LNTZER, LUSEVT, LUSTKB, LTRKBN CHARACTER*10 TITUSB COMMON /UUSRBN/ XLOW (MXUSBN), XHIGH (MXUSBN), YLOW (MXUSBN), & YHIGH (MXUSBN), ZLOW (MXUSBN), ZHIGH (MXUSBN), & DXUSBN(MXUSBN), DYUSBN(MXUSBN), DZUSBN(MXUSBN), & TCUSBN(MXUSBN), BKUSBN(MXUSBN), B2USBN(MXUSBN), & NXBIN (MXUSBN), NYBIN (MXUSBN), NZBIN (MXUSBN), & ITUSBN(MXUSBN), IDUSBN(MXUSBN), KBUSBN(MXUSBN), & IPUSBN(MXUSBN), LEVTBN(MXUSBN), LNTZER(MXUSBN), & LTRKBN(MXUSBN), NUSRBN, LUSBIN, LUSEVT, LUSTKB COMMON /USRCH/ TITUSB(MXUSBN) PARAMETER ( PIPIPI = 3.141592653589793D+00 ) PARAMETER ( MXZTA = 3000 ) PARAMETER ( MXRAD = 3000 ) PARAMETER ( MEMORY = 40000000 ) DIMENSION GMSTOR (MEMORY), GESTOR (MEMORY), ACCDIS (MXUSBN), & JB (MXUSBN), RHO (MEMORY), GMRZTA (MXRAD,MXZTA), & VLRZTA(MXRAD,MXZTA), FLUMAT(MXXXXX), FLEMAT(MXXXXX), & GERZTA (MXRAD,MXZTA) DOUBLE PRECISION ACCDIS, ACCUDS, DOSMAX, DOSMIN, DOSINT, & DOSECO, VOLUCO, ZTAZTA, GMRZTA, VLRZTA, ZTALOW, ZTAHGH, & DZTABN, RADLOW, RADHGH, DRADBN, XXXMIN, XXXMAX, & YYYMIN, YYYMAX, ZZZMIN, ZZZMAX, XCOLOW, XCOHGH, DXCOBN, & YCOLOW, YCOHGH, DYCOBN, ZCOLOW, ZCOHGH, DZCOBN, & WEIGXT, WEIGX2, WEIGYT, WEIGY2, WEIGZT, WEIGZ2, XCSXCS, YUIYUI, & PARVOL, VOLCHK, VOLPAR, VOLPA0, PARVHL, DOSERR, GERZTA LOGICAL LCYL, LFIRST, LRHOFL, LGEOPL, LRADBN, LSTAT, LRPHI CHARACTER CARD*80 PARAMETER ( GEVGRA = 1.60219E-07 ) * For the Atlas TP: 1E34 adopted as peak luminosity PARAMETER ( PLUMIN = 1.0 E+00 ) * PARAMETER ( PLUMIN = 1.0 E+34 ) ** PARAMETER ( PLUMIN = 1.7 E+34 ) PARAMETER ( ALUMIN = 1.0 E+00 ) * PARAMETER ( ALUMIN = 1.0 E+34 ) PARAMETER ( SIGIN = 1.0 E+00 ) * PARAMETER ( SIGIN = 8.0 E-26 ) PARAMETER ( TIMEYR = 1.0 E+00 ) * PARAMETER ( TIMEYR = 1.0 E+07 ) PARAMETER ( TIMKHZ = 1.0 E+00 ) * PARAMETER ( TIMKHZ = 1.0 E-03 ) * DO I = 1, MEMORY GMSTOR (I) = 0.D+00 GESTOR (I) = 0.D+00 RHO (I) = 0.E+00 END DO * IF ( .NOT. LRHOFL ) THEN IOPT = 1 ULUMIN = PLUMIN TIMUSD = TIMKHZ * ULUMIN = ALUMIN * TIMUSD = TIMEYR ELSE ULUMIN = ALUMIN TIMUSD = TIMEYR END IF SITILU = SIGIN*ULUMIN*TIMUSD LFIRST = .TRUE. IB0 = 0 LSTAT = .FALSE. IRECRD = 0 READ (11,ERR=120) RUNTIT,RUNTIM,WEIPRI,NCASE,MCASE,NBATCH IRECRD = IRECRD + 1 WRITE(*,*)RUNTIT WRITE(*,*)RUNTIM WRITE(*,*)WEIPRI WRITE(*,*)MCASE, NCASE WRITE(*,*)NBATCH GO TO 150 100 CONTINUE IRECRD = 0 REWIND 11 READ (11,ERR=100) RUNTIT,RUNTIM,WEIPRI,NCASE,NBATCH IRECRD = IRECRD + 1 MCASE = 0 WRITE(*,*)RUNTIT WRITE(*,*)RUNTIM WRITE(*,*)WEIPRI WRITE(*,*)NCASE WRITE(*,*)NBATCH GO TO 150 120 CONTINUE IRECRD = 0 REWIND 11 READ (11,ERR=100) RUNTIT,RUNTIM,WEIPRI,NCASE IRECRD = IRECRD + 1 NBATCH = 1 MCASE = 0 WRITE(*,*)RUNTIT WRITE(*,*)RUNTIM WRITE(*,*)WEIPRI WRITE(*,*)NCASE 150 CONTINUE ACCUTT = 0.D0 ACCUDS = 0.D0 KLAST = 0 LCYL = .FALSE. NBLAST = 0 * +-------------------------------------------------------------------* * | Read Binnings DO 1000 IB=IB0+1,1000 NB=IB READ (11,ERR=400,END=2000) MB, TITUSB(NB), ITUSBN(NB), & IDUSBN(NB), & XLOW(NB),XHIGH(NB),NXBIN(NB), & DXUSBN(NB),YLOW(NB), & YHIGH(NB),NYBIN(NB),DYUSBN(NB), & ZLOW(NB),ZHIGH(NB),NZBIN(NB), & DZUSBN(NB),LNTZER(NB),BKUSBN(NB), & B2USBN(NB),TCUSBN(NB) IRECRD = IRECRD + 1 GO TO 500 400 CONTINUE * BACKSPACE 11 REWIND 11 DO IR = 1, IRECRD READ (11) END DO READ (11,ERR=450,END=2000) MB, TITUSB(NB), ITUSBN(NB), & IDUSBN(NB), & XLOW(NB),XHIGH(NB),NXBIN(NB), & DXUSBN(NB),YLOW(NB), & YHIGH(NB),NYBIN(NB),DYUSBN(NB), & ZLOW(NB),ZHIGH(NB),NZBIN(NB), & DZUSBN(NB) IRECRD = IRECRD + 1 LNTZER(NB)=.FALSE. TCUSBN(NB)=1.E+38 BKUSBN(NB)=0.E+00 B2USBN(NB)=0.E+00 GO TO 500 450 CONTINUE * BACKSPACE 11 REWIND 11 DO IR = 1, IRECRD READ (11) END DO READ (11,ERR=460) CHSTAT,IJX IRECRD = IRECRD + 1 IF ( CHSTAT .EQ. 'STATISTICS' ) THEN LSTAT = .TRUE. GO TO 2000 END IF 460 CONTINUE STOP 'UNKNOWN RECORD TYPE' 500 CONTINUE NBLAST = NB WRITE(*,*)' NB',NB,' ',TITUSB(NB) * | +----------------------------------------------------------------* * | | Requested bin ask user for limits IF ( NB .EQ. NBG ) THEN IF ( MOD (ITUSBN(NB),10) .EQ. 1 .OR. & MOD (ITUSBN(NB),10) .EQ. 7 ) THEN LRADBN = .TRUE. ELSE LRADBN = .FALSE. ENDIF IRAXIS = 0 IZAXIS = 0 WRITE (*,*) XLOW (NB), XHIGH(NB), NXBIN (NB) IF ( LRADBN ) THEN WRITE (*,*)' Rmin, Rmax (def. provided)' ELSE WRITE (*,*)' Xmin, Xmax (def. provided)' ENDIF XCOLOW = XLOW (NB) XCOHGH = XHIGH(NB) READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) XCOLOW, XCOHGH IF ( LRADBN ) THEN WRITE (*,*)' NR (def. N_R_max)' ELSE WRITE (*,*)' Nx (def. N_x_max)' END IF NXCOBN = NXBIN (NB) READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) NXCOBN DXCOBN = ( XCOHGH - XCOLOW ) / DBLE (NXCOBN) IF ( NXCOBN .GT. 1 ) THEN IRAXIS = 1 RADLOW = XCOLOW RADHGH = XCOHGH NRADBN = NXCOBN DRADBN = DXCOBN END IF WRITE (*,*) YLOW (NB), YHIGH(NB), NYBIN (NB) IF ( LRADBN ) THEN WRITE (*,*)' PHImin, PHImax (def. provided)' ELSE WRITE (*,*)' Ymin, Ymax (def. provided)' ENDIF YCOLOW = YLOW (NB) YCOHGH = YHIGH(NB) READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) YCOLOW, YCOHGH IF ( LRADBN ) THEN WRITE (*,*)' NPHI (def. N_PHI_max)' ELSE WRITE (*,*)' Ny (def. N_y_max)' ENDIF NYCOBN = NYBIN (NB) READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) NYCOBN DYCOBN = ( YCOHGH - YCOLOW ) / DBLE (NYCOBN) IF ( NYCOBN .GT. 1 .AND. IRAXIS .LE. 0 ) THEN IRAXIS = 2 RADLOW = YCOLOW RADHGH = YCOHGH NRADBN = NYCOBN DRADBN = DYCOBN ELSE IF ( NYCOBN .GT. 1 .AND. IRAXIS .GT. 0 ) THEN IZAXIS = 2 ZTALOW = YCOLOW ZTAHGH = YCOHGH NZTABN = NYCOBN DZTABN = DYCOBN END IF WRITE (*,*) ZLOW (NB), ZHIGH(NB), NZBIN (NB) WRITE (*,*)' Zmin, Zmax (def. provided)' ZCOLOW = ZLOW (NB) ZCOHGH = ZHIGH(NB) READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) ZCOLOW, ZCOHGH WRITE (*,*)' Nz (def. N_z_max)' NZCOBN = NZBIN (NB) READ (*,'(A)') CARD IF ( CARD .NE. ' ' ) READ (CARD,*) NZCOBN DZCOBN = ( ZCOHGH - ZCOLOW ) / DBLE (NZCOBN) LRPHI = LRADBN .AND. NZCOBN .EQ. 1 IF ( NZCOBN .GT. 1 .AND. IZAXIS .LE. 0 ) THEN IZAXIS = 3 ZTALOW = ZCOLOW ZTAHGH = ZCOHGH NZTABN = NZCOBN DZTABN = DZCOBN ELSE IF ( NZCOBN .LE. 1 .AND. IZAXIS .GT. 0 ) THEN END IF * One dimensional scans * Z axis: IF (NXCOBN .EQ. 1 .AND. NYCOBN .EQ. 1) THEN IRAXIS = 0 IZAXIS = 3 ZTALOW = ZCOLOW ZTAHGH = ZCOHGH NZTABN = NZCOBN DZTABN = DZCOBN NRADBN = 1 * Y axis: ELSE IF (NXCOBN .EQ. 1 .AND. NZCOBN .EQ. 1) THEN IRAXIS = 0 IZAXIS = 2 ZTALOW = YCOLOW ZTAHGH = YCOHGH NZTABN = NYCOBN DZTABN = DYCOBN NRADBN = 1 * X axis: ELSE IF (NYCOBN .EQ. 1 .AND. NZCOBN .EQ. 1) THEN IRAXIS = 0 IZAXIS = 1 ZTALOW = XCOLOW ZTAHGH = XCOHGH NZTABN = NXCOBN DZTABN = DXCOBN NRADBN = 1 ELSE IF (NXCOBN .EQ. 1 .AND. NYCOBN .EQ. 1.AND. & NZCOBN .EQ. 1) THEN STOP ' ABSURD BINNING' END IF WRITE (*,*)' IRAXIS,IZAXIS',IRAXIS,IZAXIS WRITE (*,*)' ZTALOW,ZTAHGH,DZTABN',ZTALOW,ZTAHGH,DZTABN WRITE (*,*)' RADLOW,RADHGH,DRADBN',RADLOW,RADHGH,DRADBN DO IE=1,NZTABN DO IA=1,NRADBN VLRZTA (IA,IE) = 0.D+00 GMRZTA (IA,IE) = 0.D+00 GERZTA (IA,IE) = 0.D+00 END DO END DO END IF * | | * | +----------------------------------------------------------------* JB(IB) = NB ACCDIS(NB) = 0.D+00 K0 = KLAST + 1 K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1 KLAST = K1 KBUSBN(NB) = K0 * | Read data READ (11) (GMSTOR(J), J = K0, K1) IRECRD = IRECRD + 1 * | +----------------------------------------------------------------* * | | IF ( NB .EQ. NBG .AND. LFIRST ) THEN LFIRST =.FALSE. KK0 = 1 KK1 = K1 - K0 + 1 IF ( LRHOFL .AND. RHOFIX .LT. 0.E+00 ) THEN READ(12) READ(12) READ(12) (RHO(J), J = KK0, KK1) CLOSE (UNIT=12) ELSE DO J = KK0,KK1 RHO (J) = RHOFIX END DO END IF END IF * | | * | +----------------------------------------------------------------* 1000 CONTINUE * | * +-------------------------------------------------------------------* 2000 CONTINUE * +-------------------------------------------------------------------* * | Read statistics IF ( LSTAT ) THEN KLAST = 0 * | +----------------------------------------------------------------* * | | Read errors DO IB = 1, NBLAST K0 = KLAST + 1 K1 = NXBIN(IB) * NYBIN(IB) * NZBIN(IB) + K0 - 1 KLAST = K1 READ (11) (GESTOR(J), J = K0, K1) IRECRD = IRECRD + 1 END DO * | | * | +----------------------------------------------------------------* END IF * | * +-------------------------------------------------------------------* NB = NBG IF ( LRADBN ) THEN PARVO0 = 0.5D+00 * DYUSBN(NB) * DZUSBN(NB) * DXUSBN(NB)**2 PARVHL = 2.D+00 * XLOW(NB) / DXUSBN (NB) ELSE PARVOL = DXUSBN(NB) * DYUSBN(NB) * DZUSBN(NB) END IF * +-------------------------------------------------------------------* * | DO 7540 IX = 1,NXBIN(NB) * | +----------------------------------------------------------------* * | | IF ( LRADBN ) THEN PARVOL = PARVO0 * ( DBLE ( 2*IX - 1 ) + PARVHL ) IF ( PARVOL .LE. 0.0001 ) THEN WRITE(*,*)'PARVOL : PARVO0,PARVHL ,ix ' WRITE(*,*)PARVOL,PARVO0,PARVHL ,ix END IF END IF * | | * | +----------------------------------------------------------------* * | +----------------------------------------------------------------* * | | DO 7530 IY = 0,NYBIN(NB)-1 * | | +-------------------------------------------------------------* * | | | DO 7520 IZ = 0,NZBIN(NB)-1 IXYZ0 = ( IZ * NYBIN (NB) + IY ) * NXBIN (NB) & + IX IXYZ = IXYZ0 + KBUSBN (NB) - 1 XXX = XLOW(NB)+(FLOAT(IX)-0.5E+00)*DXUSBN(NB) YYY = YLOW(NB)+(FLOAT(IY)+0.5E+00)*DYUSBN(NB) ZZZ = ZLOW(NB)+(FLOAT(IZ)+0.5E+00)*DZUSBN(NB) * | | | +----------------------------------------------------------* * | | | | IF ( ABS (GMSTOR (IXYZ)) .GT. 0.E+00 ) THEN IF ( RHO (IXYZ0) .LE. 1.E-05 .AND. LRHOFL ) THEN WRITE (77,*)' RHO,IXYZ0,X,Y,Z,GMSTOR,NB', & RHO(IXYZ0),IXYZ0,XXX,YYY,ZZZ,GMSTOR(IXYZ),NB ACCDIS(NB) = ACCDIS(NB) + GMSTOR (IXYZ) * PARVOL GMSTOR (IXYZ) = 0.E+00 GESTOR (IXYZ) = 0.E+00 DOSECO = 0.D+00 DOSERR = 0.D+00 VOLUCO = 0.D+00 ELSE IF ( RHO (IXYZ0) .LE. RHOTH ) THEN ACCDIS(NB) = ACCDIS(NB) + GMSTOR (IXYZ) * PARVOL GMSTOR (IXYZ) = 0.E+00 DOSECO = 0.D+00 DOSERR = 0.D+00 VOLUCO = 0.D+00 ELSE GMSTOR (IXYZ) = GMSTOR (IXYZ) / RHO (IXYZ0) & * SITILU IF ( LRHOFL ) THEN GMSTOR (IXYZ) = GMSTOR (IXYZ) * GEVGRA END IF DOSECO = DBLE(GMSTOR(IXYZ)) * PARVOL & * RHO (IXYZ0) DOSERR = (DBLE(GESTOR(IXYZ)) * DOSECO)**2 IF ( ABS (DOSECO) .LT. 1.D-70 ) THEN & DOSECO = SIGN ( 1.D-70, DOSECO ) VOLUCO = DBLE (PARVOL*RHO(IXYZ0)) END IF * | | | | * | | | +----------------------------------------------------------* * | | | | ELSE IF ( RHO (IXYZ0) .GT. 1.E-05 ) THEN GMSTOR (IXYZ) = 0.E+00 GESTOR (IXYZ) = 0.E+00 DOSECO = 1.D-70 DOSERR = 0.D+00 VOLUCO = DBLE (PARVOL*RHO(IXYZ0)) ELSE GMSTOR (IXYZ) = 0.E+00 GESTOR (IXYZ) = 0.E+00 DOSECO = 0.D+00 DOSERR = 0.D+00 VOLUCO = 0.D+00 END IF END IF * | | | | * | | | +----------------------------------------------------------* * | | | +----------------------------------------------------------* * | | | | IF ( ABS (DOSECO) .GT. 0.D+00 ) THEN IF ( ABS (DOSECO) .LT. 1.D-60 ) DOSECO = 0.D+00 RHOVOL = 0.D+00 VOLCHK = 0.D+00 * LPOINT=1 LPOINT=3 * LPOINT=10 * MPOINT=1 MPOINT=3 * MPOINT=10 * NPOINT=1 NPOINT=3 * NPOINT=10 WEIGXT=1.D+00/DBLE(LPOINT) WEIGX2=0.5D+00*WEIGXT WEIGYT=1.D+00/DBLE(MPOINT) WEIGY2=0.5D+00*WEIGYT WEIGZT=1.D+00/DBLE(NPOINT) WEIGZ2=0.5D+00*WEIGZT IF ( LRADBN ) THEN VOLPA0 = WEIGZT * DZUSBN(NB) * WEIGY2 & * DYUSBN (NB) ELSE VOLPA0 = WEIGZT * WEIGXT * WEIGYT * PARVOL END IF * | | | | +-------------------------------------------------------* * | | | | | DO 7510 JX=1,LPOINT XXXMIN = XLOW(NB)+(DBLE(IX-1)+WEIGXT*DBLE(JX-1)) & * DXUSBN(NB) XXXMAX = XLOW(NB)+(DBLE(IX-1)+WEIGXT*DBLE(JX)) & * DXUSBN(NB) IF ( LRADBN ) THEN VOLPAR = VOLPA0 * ( XXXMAX - XXXMIN ) & * ( XXXMAX + XXXMIN ) ELSE VOLPAR = VOLPA0 END IF * | | | | | +----------------------------------------------------* * | | | | | | DO 7508 JY=1,MPOINT YYYMIN = YLOW(NB)+(DBLE(IY) & + WEIGYT*DBLE(JY-1)) * DYUSBN(NB) YYYMAX = YLOW(NB)+(DBLE(IY) & + WEIGYT*DBLE(JY)) * DYUSBN(NB) * | | | | | | +-------------------------------------------------* * | | | | | | | DO 7505 JZ=1,NPOINT ZZZMIN = ZLOW(NB)+(DBLE(IZ) & + WEIGZT*DBLE(JZ-1)) * DZUSBN(NB) ZZZMAX = ZLOW(NB)+(DBLE(IZ) & + WEIGZT*DBLE(JZ)) * DZUSBN(NB) XCSXCS = 0.5D+00 * ( XXXMIN + XXXMAX ) YUIYUI = 0.5D+00 * ( YYYMIN + YYYMAX ) ZTAZTA = 0.5D+00 * ( ZZZMIN + ZZZMAX ) VOLCHK = VOLCHK + VOLPAR IF ( XCSXCS .LT. XCOLOW .OR. XCSXCS .GE. & XCOHGH .OR. YUIYUI .LT. YCOLOW .OR. & YUIYUI .GE. YCOHGH .OR. ZTAZTA .LT. & ZCOLOW .OR. ZTAZTA .GE. ZCOHGH ) & GO TO 7505 IF ( IRAXIS .EQ. 0 ) THEN IRAD = 1 ELSE IF ( IRAXIS .EQ. 1 ) THEN IRAD = INT ( ( XCSXCS - RADLOW ) & / DRADBN ) + 1 ELSE IF ( IRAXIS .EQ. 2 ) THEN IRAD = INT ( ( YUIYUI - RADLOW ) & / DRADBN ) + 1 ELSE STOP ' IRAXIS??' END IF IF ( IZAXIS .EQ. 1 ) THEN IZTA = INT ( ( XCSXCS - ZTALOW ) & / DZTABN ) + 1 ELSE IF ( IZAXIS .EQ. 2 ) THEN IZTA = INT ( ( YUIYUI - ZTALOW ) & / DZTABN ) + 1 ELSE IF ( IZAXIS .EQ. 3 ) THEN IZTA = INT ( ( ZTAZTA - ZTALOW ) & / DZTABN ) + 1 ELSE STOP ' IZAXIS??' END IF IF ( IRAD .LE. 0 ) GO TO 7505 IF ( IZTA .LE. 0 ) GO TO 7505 IRAD = MIN ( IRAD, NRADBN ) IZTA = MIN ( IZTA, NZTABN ) GMRZTA (IRAD,IZTA) = GMRZTA (IRAD,IZTA) & + DOSECO * VOLPAR / PARVOL GERZTA (IRAD,IZTA) = GERZTA (IRAD,IZTA) & + DOSERR * VOLPAR / PARVOL VLRZTA (IRAD,IZTA) = VLRZTA (IRAD,IZTA) & + VOLUCO * VOLPAR / PARVOL 7505 CONTINUE * | | | | | | | * | | | | | | +-------------------------------------------------* 7508 CONTINUE * | | | | | | * | | | | | +----------------------------------------------------* 7510 CONTINUE * | | | | | * | | | | +-------------------------------------------------------* IF (ABS(PARVOL-VOLCHK) .GT. 1.D-05*PARVOL) & WRITE(*,*)' IX,IY,IZ,NB,PARVOL,VOLCHK', & IX,IY,IZ,NB,PARVOL,VOLCHK END IF * | | | | * | | | +----------------------------------------------------------* 7520 CONTINUE * | | | * | | +-------------------------------------------------------------* 7530 CONTINUE * | | * | +----------------------------------------------------------------* 7540 CONTINUE * | * +-------------------------------------------------------------------* DOSMAX = -1.D+70 DOSMIN = +1.D+70 DOSINT = 0.D+00 DO IE=1,NZTABN DO IA=1,NRADBN DOSINT = DOSINT + GMRZTA (IA,IE) IF ( ABS (GMRZTA (IA,IE)) .GT. 0.D+00 ) THEN GERZTA (IA,IE) = SQRT (GERZTA(IA,IE)) / GMRZTA (IA,IE) GMRZTA (IA,IE) = GMRZTA (IA,IE) / VLRZTA (IA,IE) END IF IDOS = (IA-1)*NZTABN + IE * Introduce here the ABS!! FLUMAT (IDOS) = GMRZTA (IA,IE) FLEMAT (IDOS) = GERZTA (IA,IE) IF ( DBLE (ABS(FLUMAT (IDOS))) .LT. 0.999999D+38 ) & DOSMAX = MAX ( DOSMAX, DBLE (ABS(FLUMAT (IDOS))) ) DOSMIN = MIN ( DOSMIN, DBLE (ABS(FLUMAT (IDOS))) ) END DO END DO WRITE (*,*) ' Dosmax, Dosmin: ', DOSMAX, DOSMIN WRITE (*,*) ' Binning integral: ', DOSINT DO IE=1,NZTABN DO IA=1,NRADBN IDOS = (IA-1)*NZTABN + IE IF ( DBLE (ABS(FLUMAT(IDOS))) .GT. 0.D+00 .AND. & DBLE (ABS(FLUMAT(IDOS))) .LT. 1.D-25*DOSMAX ) THEN WRITE (*,*)' *** Zeroed!',FLUMAT(IDOS),SNGL(DOSMAX), & IE,IA FLUMAT (IDOS) = 0.E+00 END IF END DO END DO * Start_VAX_seq * OPEN(UNIT=33,FILE='DOSLEV.DAT',STATUS='NEW') * End_VAX_seq * Start_UNIX_seq OPEN(UNIT=33,FILE='doslev.dat',STATUS='UNKNOWN') * End_UNIX_seq WRITE(33,8700)((FLUMAT((IA-1)*NZTABN+IE),IA=1,NRADBN),IE=1,NZTABN) 8700 FORMAT (1(5X,1P,5(1X,E11.4))) CLOSE(UNIT=33) IB0=IB-1 SRDLOW=RADLOW SRDHGH=RADHGH SZTLOW=ZTALOW SZTHGH=ZTAHGH RETURN *=== End of subroutine Ubdos2 =========================================* END