Completed first version of merge utility. Small additional changes and bugfixes.

This commit is contained in:
Yggdrasil 2026-04-28 18:18:31 +02:00
parent 01c9bf2018
commit c7b344530b
5 changed files with 401 additions and 191 deletions

View File

@ -19,8 +19,8 @@ TARGET1 = event2
TARGET2 = event2_merge TARGET2 = event2_merge
# Source files # Source files
SRCS1 = event2_03.f gbook.f SRCS1 = event2_03.f gbook.f utilities.f
SRCS2 = event2_merge.f gbook.f SRCS2 = event2_merge.f gbook.f utilities.f
# Executable paths # Executable paths
EXE1 = $(OUT)/$(TARGET1) EXE1 = $(OUT)/$(TARGET1)

View File

@ -346,188 +346,6 @@ C Close file.
END 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----------------------------------------------------------------------- C-----------------------------------------------------------------------
SUBROUTINE DEMOIN SUBROUTINE DEMOIN
IMPLICIT NONE IMPLICIT NONE

View File

@ -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 PROGRAM MAIN
CALL GREAD('./results/gtopdraw_THRUST_00000.dat') IMPLICIT NONE
CALL GTOPER( 3,1,1,0,100+3) INTEGER I
CALL GTOPER(10+3,0,0,0,110+3) 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 END PROGRAM
C----------------------------------------------------------------------- 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-----------------------------------------------------------------------

40
gbook.f
View File

@ -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 NOT 1., GIVE FACTORS BY WHICH THE ID1 AND ID2 BIN CONTENTS ARE
C MULTIPLIED BEFORE THE INDICATED OPERATION. (DIVISION WITH A C MULTIPLIED BEFORE THE INDICATED OPERATION. (DIVISION WITH A
C VANISHING BIN CONTENT WILL GIVE 0.) C VANISHING BIN CONTENT WILL GIVE 0.)
C OPER= 'A', 'S', 'L': FOR 'S' THE SQUARE ROOT OF THE CONTENT IN ID1 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) AND FOR 'L' THE 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 10-LOGARITHM IS TAKEN (A NONPOSITIVE BIN CONTENT IS BEFORE THAT
C REPLACED BY 0.8 TIMES THE SMALLEST POSITIVE BIN CONTENT). C REPLACED BY 0.8 TIMES THE SMALLEST POSITIVE BIN CONTENT).
C THEREAFTER, IN ALL THREE CASES, THE CONTENT IS MULTIPLIED BY F1 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 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 OPER= 'M': INTENDED FOR STATISTICAL ANALYSIS, BIN-BY-BIN MEAN AND
C STANDARD DEVIATION OF A VARIABLE, ASSUMING THAT ID1 CONTAINS C STANDARD DEVIATION OF A VARIABLE, ASSUMING THAT ID1 CONTAINS
C ACCUMULATED WEIGHTS, ID2 ACCUMULATED WEIGHT*VARIABLE AND 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)= 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)) &SQRT(MAX(A(IS3+IC)/A(IS1+IC)-A(IS2+IC)**2,ZERO))
180 A(IS1+IC)=F1*A(IS1+IC) 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 ENDIF
RETURN RETURN
END END
@ -931,6 +937,34 @@ C-----------------------------------------------------------------------
GOTO 200 GOTO 200
ENDIF ENDIF
END 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----------------------------------------------------------------------- C-----------------------------------------------------------------------
SUBROUTINE GREAD(FILE) SUBROUTINE GREAD(FILE)
IMPLICIT INTEGER (I-N) IMPLICIT INTEGER (I-N)

181
utilities.f Normal file
View File

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