diff --git a/Makefile b/Makefile index 3285b70..a2102b7 100644 --- a/Makefile +++ b/Makefile @@ -19,8 +19,8 @@ TARGET1 = event2 TARGET2 = event2_merge # Source files -SRCS1 = event2_03.f gbook.f -SRCS2 = event2_merge.f gbook.f +SRCS1 = event2_03.f gbook.f utilities.f +SRCS2 = event2_merge.f gbook.f utilities.f # Executable paths EXE1 = $(OUT)/$(TARGET1) diff --git a/event2_03.f b/event2_03.f index 2b58af1..604608e 100644 --- a/event2_03.f +++ b/event2_03.f @@ -346,188 +346,6 @@ C Close file. END -c----------------------------------------------------------------------- -c Auxilliary subroutines. -c----------------------------------------------------------------------- - - subroutine getline(unit, line, stat) - implicit none - integer, intent(in) :: unit - integer, intent(out) :: stat - character(72), intent(out) :: line - integer :: size - integer :: i,j - integer :: stat2 - character(72) :: buffer - character(2) :: pattern -c List of characers where blanks after/before will be eliminated. - character(*), parameter :: killtrail = "=,>[*+" - character(*), parameter :: killlead = "=,>]*+" - -c Read the full line. - line = '' - do - read(unit, "(A)", iostat=stat) line - if (stat > 0) return - exit - end do - -c Replace all `tab` characters by a blank. - do - i = index(line, char(9)) - if (i.eq.0) exit - line(i:i) = " " - end do - -c Kill leading blanks. - line = trim(adjustl( line )) -c Kill possible comments. - i = index(line, "!") -c Kill trailing blanks. - if (i.gt.0) line = trim(adjustl(line(:i-1))) - -c Kill blanks before special characters. - do j=1,len(killlead) - pattern = ' ' // killlead(j:j) - do - i = index(line,pattern) - if (i.eq.0) exit - line = line(:i-1) // killtrail(j:j) // line(i+2:) - end do - end do - -c Kill blanks after special characters. - do j=1,len(killlead) - pattern = killlead(j:j) // ' ' - do - i = index(line,pattern) - if (i.eq.0) exit - line = line(:i-1) // killlead(j:j) // line(i+2:) - end do - end do - - return - end - -************************************************************************ - - subroutine readmode(cmode, var, def) - implicit none - character(12), intent(in) :: cmode - integer, intent(in) :: def - integer, intent(out) :: var - integer :: i, imode - character(12) :: keys(20),settings(20) -c Common blocks. - common/runcard/keys,settings - -c Try to find mode with name 'cmode' in settings. - imode = -1 - do i=1,20 - if (keys(i).eq.cmode)then - imode = i - exit - endif - end do -c If not found, set to default. - if (imode.lt.0) var = def -c Otherwise set to value present in settings. - if (imode.ge.0) call readInt(settings(imode), var) - - return - end - -************************************************************************ - - subroutine readparm(cparm, var, def) - implicit none - character(12), intent(in) :: cparm - real(8), intent(in) :: def - real(8), intent(out) :: var - integer :: i,iparm - character(12) :: keys(20),settings(20) -c Common blocks. - common/runcard/keys,settings - -c Try to find mode with name 'cparm' in settings. - iparm = -1 - do i=1,20 - if (keys(i).eq.cparm) iparm = i - end do -c If not found, set to default. - if (iparm.lt.0) var = def -c Otherwise set to value present in settings. - if (iparm.ge.0) read(settings(iparm),*) var - - return - end - -************************************************************************ - -c Auxiliary helper subroutine to read integers in different formats. - - subroutine readint(string,var) - implicit none - integer, intent(out) :: var - character(8), intent(in) :: string - integer :: iposk,iposm,ipose,iposd - real(8) :: helper - - iposk = index(string,'k') - if (iposk.eq.0) iposk = index(string,'K') - iposm = index(string,'m') - if (iposm.eq.0) iposm = index(string,'M') - ipose = index(string,'e') - if (ipose.eq.0) ipose = index(string,'E') - iposd = index(string,'d') - if (iposd.eq.0) iposd = index(string,'D') - - if (iposk.ne.0)then - read(string(1:iposk-1),'(I16)') var - var = 1000*var - elseif (iposm.ne.0)then - read(string(1:iposm-1),'(I16)') var - var = 1000000*var - elseif (ipose.ne.0 .or. iposd.ne.0)then - read(string,'(F16.0)') helper - var = helper - else - read(string, '(I16)') var - endif - - return - end - -************************************************************************ - -c Auxiliary helper subroutine to create output folder - - SUBROUTINE CREATEOUT(PATH) - IMPLICIT NONE - CHARACTER*40 PATH - - INTEGER L,EXITSTAT - CHARACTER*40 OPATH - COMMON /FPATH/ OPATH - - OPATH = PATH - L = LEN_TRIM(OPATH) - - IF (L == 0) THEN - OPATH = './' - ELSE IF (OPATH(L:L) /= '/') THEN - OPATH(L+1:L+1) = '/' - END IF - - CALL EXECUTE_COMMAND_LINE("mkdir -p " // TRIM(OPATH), - $ exitstat=EXITSTAT) - IF (EXITSTAT /= 0) THEN - PRINT *, "Error creating directory:", TRIM(OPATH) - STOP - END IF - - END - C----------------------------------------------------------------------- SUBROUTINE DEMOIN IMPLICIT NONE diff --git a/event2_merge.f b/event2_merge.f index 5fdec9b..73761a3 100644 --- a/event2_merge.f +++ b/event2_merge.f @@ -1,6 +1,183 @@ +C----------------------------------------------------------------------- +C --- HISTOGRAM MERGING UTILITY, VERSION 0.0 +C --- DESIGNED TO BE USED IN CONJUNCTION WITH EERAD2 V0.4 +C +C VERSION 0.0 - 28th April 2026 +C +C WRITTEN BY G. CHIURATO +C PROGRAM MAIN - CALL GREAD('./results/gtopdraw_THRUST_00000.dat') - CALL GTOPER( 3,1,1,0,100+3) - CALL GTOPER(10+3,0,0,0,110+3) + IMPLICIT NONE + INTEGER I + CHARACTER*256 IPATH, OPATH, LPATH + CHARACTER*30 OBS(5) + COMMON /CONST/OBS + +C --- TODO: ADD EXTERNALLY CONFIG. PATHS + IPATH = './results/' + OPATH = './merge/' + LPATH = TRIM(OPATH)//'lists/' + + CALL PRINTHEADER + CALL CREATEOUT(LPATH) + CALL CREATEOUT(OPATH) + CALL SCANDIR(IPATH, LPATH) + DO 10 I = 1,5 +C ------ TODO: FIXED IDS ARE NOT GREAT + CALL MERGE(OBS(I), I, OPATH, LPATH) + 10 CONTINUE END PROGRAM -C----------------------------------------------------------------------- \ No newline at end of file +C----------------------------------------------------------------------- + SUBROUTINE SCANDIR(IPATH, OPATH) + IMPLICIT NONE + INTEGER I + INTEGER STAT + CHARACTER*30 OBS(5), ESCOBS + CHARACTER*256 IPATH, OPATH + COMMON /CONST/OBS + DATA OBS/'C-PARAMETER','D-PARAMETER', + $ 'THRUST','y3(JADE,P)','EEC'/ + +C CREATE FILE LIST USING SYSTEM COMMAND + DO 10 I = 1, 5 + CALL BRACKETESC(OBS(I), ESCOBS) + CALL SYSTEM('ls '//TRIM(IPATH)//'gtopdraw_'//TRIM(ESCOBS)// + $ '_* > '//TRIM(OPATH)//'filelist_'//TRIM(ESCOBS)//'.txt', + $ STATUS=STAT) +10 CONTINUE + + END +C----------------------------------------------------------------------- + SUBROUTINE MERGE(OBS, ID, OPATH, LPATH) + IMPLICIT NONE + CHARACTER*30 OBS + CHARACTER*256 OPATH, LPATH + CHARACTER*256 LINE + CHARACTER*60 TITLE + INTEGER LENTIT + INTEGER ID + INTEGER NX + DOUBLE PRECISION XMIN, XMAX + + NX = 0 + XMIN = 0 + XMAX = 0 + +C --- NEED TO USE UNIT 20 TO AVOID CONFLICT WITH UNIT 10 USED BY GREAD + OPEN(UNIT=20, FILE=TRIM(LPATH)//'filelist_'//TRIM(OBS)//'.txt', + $ STATUS='OLD', ERR=30) + +C --- READ FIRST HIST AND SAVE BIN NUMBER AND TITLE + READ(20, '(A)', END=20) LINE + CALL GREAD(LINE) + CALL GPROP1(ID, NX, XMIN, XMAX) +C --- TODO: IN THIS WAY THE MERGED HIST HAS THE SAME NAME AS THE LAST +C OPENED HIST. NEED TO CHANGE THIS TO SOMETHING DIFFERENT + CALL GTITLE(ID, TITLE, LENTIT) + +C --- 2XX IS VALUE ACC., 3XX IS WEIGHT ACC. + CALL GBOOK1(200+ID,TITLE,NX,XMIN,XMAX) + CALL GBOOK1(300+ID,TITLE,NX,XMIN,XMAX) + CALL GBOOK1(210+ID,TITLE,NX,XMIN,XMAX) + CALL GBOOK1(310+ID,TITLE,NX,XMIN,XMAX) + +C --- ACCUMULATE WEIGHTS + CALL GOPERA(100+ID,'*',100+ID,100+ID,1D0,1D0) + CALL GOPERA(100+ID,'I', 0,100+ID,1D0,1D0) + CALL GOPERA(300+ID,'+',100+ID,300+ID,1D0,1D0) + CALL GOPERA(110+ID,'*',110+ID,110+ID,1D0,1D0) + CALL GOPERA(110+ID,'I', 0,110+ID,1D0,1D0) + CALL GOPERA(310+ID,'+',110+ID,310+ID,1D0,1D0) +C --- ACCUMULATE VALUES + CALL GOPERA( ID,'*',100+ID, ID,1D0,1D0) + CALL GOPERA(200+ID,'+', ID,200+ID,1D0,1D0) + CALL GOPERA( 10+ID,'*',110+ID, 10+ID,1D0,1D0) + CALL GOPERA(210+ID,'+', 10+ID,210+ID,1D0,1D0) + + 10 CONTINUE + CALL GRESET(ID) + CALL GRESET(100+ID) + READ(20, '(A)', END=20) LINE + CALL GREAD(LINE) +C --- ACCUMULATE WEIGHTS + CALL GOPERA(100+ID,'*',100+ID,100+ID,1D0,1D0) + CALL GOPERA(100+ID,'I', 0,100+ID,1D0,1D0) + CALL GOPERA(300+ID,'+',100+ID,300+ID,1D0,1D0) + CALL GOPERA(110+ID,'*',110+ID,110+ID,1D0,1D0) + CALL GOPERA(110+ID,'I', 0,110+ID,1D0,1D0) + CALL GOPERA(310+ID,'+',110+ID,310+ID,1D0,1D0) +C --- ACCUMULATE VALUES + CALL GOPERA( ID,'*',100+ID, ID,1D0,1D0) + CALL GOPERA(200+ID,'+', ID,200+ID,1D0,1D0) + CALL GOPERA( 10+ID,'*',110+ID, 10+ID,1D0,1D0) + CALL GOPERA(210+ID,'+', 10+ID,210+ID,1D0,1D0) + GOTO 10 + + 20 CONTINUE + CALL GRESET(ID) + CALL GRESET(100+ID) +C --- COMPUTE FINAL VALUES + CALL GOPERA(200+ID,'/',300+ID,ID,1D0,1D0) + CALL GOPERA(210+ID,'/',310+ID,10+ID,1D0,1D0) +C --- COMPUTE FINAL ERRORS + CALL GOPERA(300+ID,'I',0,300+ID,1D0,1D0) + CALL GOPERA(300+ID,'S',0,100+ID,1D0,0D0) + CALL GOPERA(310+ID,'I',0,310+ID,1D0,1D0) + CALL GOPERA(310+ID,'S',0,110+ID,1D0,0D0) +C --- WRITE OUT RESULTS + CALL GTOPER( ID,1,1,0,100+ID) + CALL GTOPER( 10+ID,0,0,0,110+ID) + + CLOSE(10) + RETURN + 30 PRINT *, 'Error opening filelist for observable '//TRIM(OBS) + STOP + END +C----------------------------------------------------------------------- + SUBROUTINE BRACKETESC(IN, OUT) +C Escapes ( and ) for Unix shell globbing + CHARACTER*(*) IN, OUT + INTEGER I, J, L + CHARACTER C + + OUT = ' ' + J = 1 + L = LEN(IN) + + DO 10 I = 1, L + C = IN(I:I) + + IF (C .EQ. '(') THEN + OUT(J:J) = '\' + J = J + 1 + OUT(J:J) = '(' + + ELSE IF (C .EQ. ')') THEN + OUT(J:J) = '\' + J = J + 1 + OUT(J:J) = ')' + + ELSE + OUT(J:J) = C + END IF + + J = J + 1 + 10 CONTINUE + + RETURN + END +C----------------------------------------------------------------------- + SUBROUTINE PRINTHEADER + IMPLICIT NONE +C --- PRINT OPENING MESSAGE + WRITE (*,'(/2A)') ' This is EVENT2 merge utility, ', + $ ' used to merge EVENT2 histogram from different runs.' + WRITE (*,'(A)') ' Results will be stored in merge folder.' + WRITE (*,'(2A)') ' If you use this program, on EVENT2 runs, ', + $ ' please reference:' + WRITE (*,'(2A)') ' S.Catani & M.H.Seymour,', + $ ' Phys. Lett. B378 (1996) 287.' + WRITE (*,'(/A)') ' Written by Giorgio Chiurato, April 2026' + WRITE (*,'(A/)') ' Version 0.0, April 2026' + END +C----------------------------------------------------------------------- \ No newline at end of file diff --git a/gbook.f b/gbook.f index 943d9a0..6f4676e 100644 --- a/gbook.f +++ b/gbook.f @@ -86,13 +86,14 @@ C CONTENTS IN ID1 AND ID2 AND PUT THE RESULT IN ID3. F1 AND F2, IF C NOT 1., GIVE FACTORS BY WHICH THE ID1 AND ID2 BIN CONTENTS ARE C MULTIPLIED BEFORE THE INDICATED OPERATION. (DIVISION WITH A C VANISHING BIN CONTENT WILL GIVE 0.) -C OPER= 'A', 'S', 'L': FOR 'S' THE SQUARE ROOT OF THE CONTENT IN ID1 -C IS TAKEN (RESULT 0 FOR NEGATIVE BIN CONTENTS) AND FOR 'L' THE +C OPER= 'A', 'S', 'L', 'I': FOR 'S' THE SQUARE ROOT OF THE CONTENT IN ID1 +C IS TAKEN (RESULT 0 FOR NEGATIVE BIN CONTENTS). FOR 'L' THE C 10-LOGARITHM IS TAKEN (A NONPOSITIVE BIN CONTENT IS BEFORE THAT C REPLACED BY 0.8 TIMES THE SMALLEST POSITIVE BIN CONTENT). C THEREAFTER, IN ALL THREE CASES, THE CONTENT IS MULTIPLIED BY F1 C AND ADDED WITH F2, AND THE RESULT IS PLACED IN ID3. THUS ID2 -C IS DUMMY IN THESE CASES. +C IS DUMMY IN THESE CASES. FOR 'I' THE INVERSE OF THE CONTENT IN +C ID1 IS COMPUTED. C OPER= 'M': INTENDED FOR STATISTICAL ANALYSIS, BIN-BY-BIN MEAN AND C STANDARD DEVIATION OF A VARIABLE, ASSUMING THAT ID1 CONTAINS C ACCUMULATED WEIGHTS, ID2 ACCUMULATED WEIGHT*VARIABLE AND @@ -324,6 +325,11 @@ C********************************************************************* IF(ID3.NE.0.AND.ABS(A(IS1+IC)).GT.1E-10) A(IS3+IC)= &SQRT(MAX(A(IS3+IC)/A(IS1+IC)-A(IS2+IC)**2,ZERO)) 180 A(IS1+IC)=F1*A(IS1+IC) + ELSEIF(OPER.EQ.'I') THEN + DO 190 IC=10,18+NC + FA2=F2*A(IS1+IC) + IF(ABS(FA2).LE.1E-10) A(IS3+IC)=0. + 190 IF(ABS(FA2).GT.1E-10) A(IS3+IC)=F1/FA2 ENDIF RETURN END @@ -931,6 +937,34 @@ C----------------------------------------------------------------------- GOTO 200 ENDIF END +C----------------------------------------------------------------------- +C --- SIMPLE UTILITY TO OBTAIN HIST PROPERTIES + SUBROUTINE GPROP1(ID, NX, XMIN, XMAX) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + PARAMETER (NXMAX=1000) + + IF (ID.GT.NMAX) THEN + NX = 0 + XMIN = 0 + XMAX = 0 + RETURN + ENDIF + + IS = NINT(A(ID+2)) + NX = NINT(A(IS+1)) + XMIN = SNGL(A(IS+2)) + XMAX = SNGL(A(IS+3)) + + IF (NX.GT.NXMAX) THEN + NX = NXMAX + ENDIF + + RETURN + END C----------------------------------------------------------------------- SUBROUTINE GREAD(FILE) IMPLICIT INTEGER (I-N) diff --git a/utilities.f b/utilities.f new file mode 100644 index 0000000..75cf7fa --- /dev/null +++ b/utilities.f @@ -0,0 +1,181 @@ +c----------------------------------------------------------------------- +c Auxilliary subroutines. +c----------------------------------------------------------------------- + + subroutine getline(unit, line, stat) + implicit none + integer, intent(in) :: unit + integer, intent(out) :: stat + character(72), intent(out) :: line + integer :: size + integer :: i,j + integer :: stat2 + character(72) :: buffer + character(2) :: pattern +c List of characers where blanks after/before will be eliminated. + character(*), parameter :: killtrail = "=,>[*+" + character(*), parameter :: killlead = "=,>]*+" + +c Read the full line. + line = '' + do + read(unit, "(A)", iostat=stat) line + if (stat > 0) return + exit + end do + +c Replace all `tab` characters by a blank. + do + i = index(line, char(9)) + if (i.eq.0) exit + line(i:i) = " " + end do + +c Kill leading blanks. + line = trim(adjustl( line )) +c Kill possible comments. + i = index(line, "!") +c Kill trailing blanks. + if (i.gt.0) line = trim(adjustl(line(:i-1))) + +c Kill blanks before special characters. + do j=1,len(killlead) + pattern = ' ' // killlead(j:j) + do + i = index(line,pattern) + if (i.eq.0) exit + line = line(:i-1) // killtrail(j:j) // line(i+2:) + end do + end do + +c Kill blanks after special characters. + do j=1,len(killlead) + pattern = killlead(j:j) // ' ' + do + i = index(line,pattern) + if (i.eq.0) exit + line = line(:i-1) // killlead(j:j) // line(i+2:) + end do + end do + + return + end + +************************************************************************ + + subroutine readmode(cmode, var, def) + implicit none + character(12), intent(in) :: cmode + integer, intent(in) :: def + integer, intent(out) :: var + integer :: i, imode + character(12) :: keys(20),settings(20) +c Common blocks. + common/runcard/keys,settings + +c Try to find mode with name 'cmode' in settings. + imode = -1 + do i=1,20 + if (keys(i).eq.cmode)then + imode = i + exit + endif + end do +c If not found, set to default. + if (imode.lt.0) var = def +c Otherwise set to value present in settings. + if (imode.ge.0) call readInt(settings(imode), var) + + return + end + +************************************************************************ + + subroutine readparm(cparm, var, def) + implicit none + character(12), intent(in) :: cparm + real(8), intent(in) :: def + real(8), intent(out) :: var + integer :: i,iparm + character(12) :: keys(20),settings(20) +c Common blocks. + common/runcard/keys,settings + +c Try to find mode with name 'cparm' in settings. + iparm = -1 + do i=1,20 + if (keys(i).eq.cparm) iparm = i + end do +c If not found, set to default. + if (iparm.lt.0) var = def +c Otherwise set to value present in settings. + if (iparm.ge.0) read(settings(iparm),*) var + + return + end + +************************************************************************ + +c Auxiliary helper subroutine to read integers in different formats. + + subroutine readint(string,var) + implicit none + integer, intent(out) :: var + character(8), intent(in) :: string + integer :: iposk,iposm,ipose,iposd + real(8) :: helper + + iposk = index(string,'k') + if (iposk.eq.0) iposk = index(string,'K') + iposm = index(string,'m') + if (iposm.eq.0) iposm = index(string,'M') + ipose = index(string,'e') + if (ipose.eq.0) ipose = index(string,'E') + iposd = index(string,'d') + if (iposd.eq.0) iposd = index(string,'D') + + if (iposk.ne.0)then + read(string(1:iposk-1),'(I16)') var + var = 1000*var + elseif (iposm.ne.0)then + read(string(1:iposm-1),'(I16)') var + var = 1000000*var + elseif (ipose.ne.0 .or. iposd.ne.0)then + read(string,'(F16.0)') helper + var = helper + else + read(string, '(I16)') var + endif + + return + end + +************************************************************************ + +c Auxiliary helper subroutine to create output folder + + SUBROUTINE CREATEOUT(PATH) + IMPLICIT NONE + CHARACTER*40 PATH + + INTEGER L,EXITSTAT + CHARACTER*40 OPATH + COMMON /FPATH/ OPATH + + OPATH = PATH + L = LEN_TRIM(OPATH) + + IF (L == 0) THEN + OPATH = './' + ELSE IF (OPATH(L:L) /= '/') THEN + OPATH(L+1:L+1) = '/' + END IF + + CALL EXECUTE_COMMAND_LINE("mkdir -p " // TRIM(OPATH), + $ exitstat=EXITSTAT) + IF (EXITSTAT /= 0) THEN + PRINT *, "Error creating directory:", TRIM(OPATH) + STOP + END IF + + END \ No newline at end of file