gchiurato
/
EVENT2
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added log binning.

This commit is contained in:
Yggdrasil 2026-05-04 17:54:01 +02:00
parent c7b344530b
commit 3a4b0fd29f
3 changed files with 120 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

93
event2_03.f
View File

@ -21,10 +21,11 @@ C VERSION 0.2 - 15th March 1996 - removed need for Lorentz boosts
C VERSION 0.3 - 10th November 1997 - improved numerical convergence
C and safety against arithmetic errors
C - added CUTUP variable
C VERSION 0.4 - 23th April 2026 - separated hist ouput
C - added runcard input
C - allowed different bin sizes
C - seeds can now be externally set
C VERSION 0.4 - 5th May 2026 - separated hist ouput
C - added runcard input
C - allowed different bin sizes
C - seeds can now be externally set
C - added log-binning capability
C
C - NEV IS THE NUMBER OF EVENTS TO GENERATE
C - EM IS THE TOTAL CENTRE-OF-MASS ENERGY
@ -194,8 +195,8 @@ C---PRINT OPENING MESSAGE
WRITE (*,'(2A)') ' S.Catani & M.H.Seymour,',
$ ' Phys. Lett. B378 (1996) 287.'
WRITE (*,'(/A)') ' Written by Mike Seymour, January 1996'
WRITE (*,'(A)') ' Edited by Giorgio Chiurato, April 2026'
WRITE (*,'(A/)') ' Version 0.4, April 2026'
WRITE (*,'(A)') ' Edited by Giorgio Chiurato, April 2026'
WRITE (*,'(A/)') ' Version 0.4, May 2026'
END
C-----------------------------------------------------------------------
SUBROUTINE READCARD
@ -219,6 +220,7 @@ C-----------------------------------------------------------------------
INTEGER CB1,CB2,DB,TB,YB,EB
DOUBLE PRECISION CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH
DOUBLE PRECISION CBS1,CBS2,DBS,TBS,YBS,EBS
INTEGER LHI
DOUBLE PRECISION ISEEDIN(2)
@ -234,6 +236,7 @@ C-----------------------------------------------------------------------
COMMON /GBINS/CB1,CB2,DB,TB,YB,EB
COMMON /GBINS/CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH
COMMON /GBINS/CBS1,CBS2,DBS,TBS,YBS,EBS
COMMON /GBINS/LHI
COMMON /ID/ RID
common /runcard/keys,settings
@ -341,6 +344,8 @@ C TODO: need to add a proper info page and maybe change the input para
call readparm('EL ', EL , -1.0D0)
call readparm('EH ', EH , 1.00D0)
call readmode('LH ', LHI , 0)
C Close file.
close(9)
@ -354,13 +359,15 @@ C-----------------------------------------------------------------------
INTEGER CB1,CB2,DB,TB,YB,EB
DOUBLE PRECISION CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH
DOUBLE PRECISION CBS1,CBS2,DBS,TBS,YBS,EBS
LOGICAL LH
DOUBLE PRECISION CSUM,CSQR,BSUM,BSQR
COMMON /DEMCOM/CSUM(8),CSQR(8),BSUM(2,6),BSQR(2,6)
COMMON /GBINS/CB1,CB2,DB,TB,YB,EB
COMMON /GBINS/CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH
COMMON /GBINS/CBS1,CBS2,DBS,TBS,YBS,EBS
COMMON /GBINS/LH
COMMON /ID/ RID
CBS1=(CH1-CL1)/CB1
CBS2=(CH2-CL2)/CB2
DBS=(DH-DL)/DB
@ -378,36 +385,36 @@ C-----------------------------------------------------------------------
BSQR(2,I)=0
ENDIF
ENDDO
CALL GBOOK1( 1,'C-PARAMETER_'//RID,CB1,CL1,CH1)
CALL GBOOK1(101,'C-PARAMETER_'//RID,CB1,CL1,CH1)
CALL GBOOK1(201,'C-PARAMETER_'//RID,CB1,CL1,CH1)
CALL GBOOK1( 11,'C-PARAMETER_'//RID,CB2,CL2,CH2)
CALL GBOOK1(111,'C-PARAMETER_'//RID,CB2,CL2,CH2)
CALL GBOOK1(211,'C-PARAMETER_'//RID,CB2,CL2,CH2)
CALL GBOOK1( 2,'D-PARAMETER_'//RID,DB,DL,DH)
CALL GBOOK1(102,'D-PARAMETER_'//RID,DB,DL,DH)
CALL GBOOK1(202,'D-PARAMETER_'//RID,DB,DL,DH)
CALL GBOOK1( 12,'D-PARAMETER_'//RID,DB,DL,DH)
CALL GBOOK1(112,'D-PARAMETER_'//RID,DB,DL,DH)
CALL GBOOK1(212,'D-PARAMETER_'//RID,DB,DL,DH)
CALL GBOOK1( 3,'THRUST_'//RID,TB,TL,TH)
CALL GBOOK1(103,'THRUST_'//RID,TB,TL,TH)
CALL GBOOK1(203,'THRUST_'//RID,TB,TL,TH)
CALL GBOOK1( 13,'THRUST_'//RID,TB,TL,TH)
CALL GBOOK1(113,'THRUST_'//RID,TB,TL,TH)
CALL GBOOK1(213,'THRUST_'//RID,TB,TL,TH)
CALL GBOOK1( 4,'y3(JADE,P)_'//RID,YB,YL,YH)
CALL GBOOK1(104,'y3(JADE,P)_'//RID,YB,YL,YH)
CALL GBOOK1(204,'y3(JADE,P)_'//RID,YB,YL,YH)
CALL GBOOK1( 14,'y3(JADE,P)_'//RID,YB,YL,YH)
CALL GBOOK1(114,'y3(JADE,P)_'//RID,YB,YL,YH)
CALL GBOOK1(214,'y3(JADE,P)_'//RID,YB,YL,YH)
CALL GBOOK1( 5,'EEC_'//RID,EB,EL,EH)
CALL GBOOK1(105,'EEC_'//RID,EB,EL,EH)
CALL GBOOK1(205,'EEC_'//RID,EB,EL,EH)
CALL GBOOK1( 15,'EEC_'//RID,EB,EL,EH)
CALL GBOOK1(115,'EEC_'//RID,EB,EL,EH)
CALL GBOOK1(215,'EEC_'//RID,EB,EL,EH)
CALL GBOOK1( 1,'C-PARAMETER_'//RID,CB1,CL1,CH1,LH)
CALL GBOOK1(101,'C-PARAMETER_'//RID,CB1,CL1,CH1,LH)
CALL GBOOK1(201,'C-PARAMETER_'//RID,CB1,CL1,CH1,LH)
CALL GBOOK1( 11,'C-PARAMETER_'//RID,CB2,CL2,CH2,LH)
CALL GBOOK1(111,'C-PARAMETER_'//RID,CB2,CL2,CH2,LH)
CALL GBOOK1(211,'C-PARAMETER_'//RID,CB2,CL2,CH2,LH)
CALL GBOOK1( 2,'D-PARAMETER_'//RID,DB,DL,DH,LH)
CALL GBOOK1(102,'D-PARAMETER_'//RID,DB,DL,DH,LH)
CALL GBOOK1(202,'D-PARAMETER_'//RID,DB,DL,DH,LH)
CALL GBOOK1( 12,'D-PARAMETER_'//RID,DB,DL,DH,LH)
CALL GBOOK1(112,'D-PARAMETER_'//RID,DB,DL,DH,LH)
CALL GBOOK1(212,'D-PARAMETER_'//RID,DB,DL,DH,LH)
CALL GBOOK1( 3,'THRUST_'//RID,TB,TL,TH,LH)
CALL GBOOK1(103,'THRUST_'//RID,TB,TL,TH,LH)
CALL GBOOK1(203,'THRUST_'//RID,TB,TL,TH,LH)
CALL GBOOK1( 13,'THRUST_'//RID,TB,TL,TH,LH)
CALL GBOOK1(113,'THRUST_'//RID,TB,TL,TH,LH)
CALL GBOOK1(213,'THRUST_'//RID,TB,TL,TH,LH)
CALL GBOOK1( 4,'y3(JADE,P)_'//RID,YB,YL,YH,LH)
CALL GBOOK1(104,'y3(JADE,P)_'//RID,YB,YL,YH,LH)
CALL GBOOK1(204,'y3(JADE,P)_'//RID,YB,YL,YH,LH)
CALL GBOOK1( 14,'y3(JADE,P)_'//RID,YB,YL,YH,LH)
CALL GBOOK1(114,'y3(JADE,P)_'//RID,YB,YL,YH,LH)
CALL GBOOK1(214,'y3(JADE,P)_'//RID,YB,YL,YH,LH)
CALL GBOOK1( 5,'EEC_'//RID,EB,EL,EH,LH)
CALL GBOOK1(105,'EEC_'//RID,EB,EL,EH,LH)
CALL GBOOK1(205,'EEC_'//RID,EB,EL,EH,LH)
CALL GBOOK1( 15,'EEC_'//RID,EB,EL,EH,LH)
CALL GBOOK1(115,'EEC_'//RID,EB,EL,EH,LH)
CALL GBOOK1(215,'EEC_'//RID,EB,EL,EH,LH)
END
C-----------------------------------------------------------------------
SUBROUTINE DEMOUT(NEV)
@ -472,12 +479,14 @@ C-----------------------------------------------------------------------
INTEGER CB1,CB2,DB,TB,YB,EB
DOUBLE PRECISION CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH
DOUBLE PRECISION CBS1,CBS2,DBS,TBS,YBS,EBS
LOGICAL LH
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /DEMCOM/CSUM,CSQR,BSUM,BSQR
COMMON /GBINS/CB1,CB2,DB,TB,YB,EB
COMMON /GBINS/CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH
COMMON /GBINS/CBS1,CBS2,DBS,TBS,YBS,EBS
COMMON /GBINS/LH
DATA CS,BS/20*0/
C---FILL THE HISTOGRAMS AT THE END OF EACH EVENT
@ -530,7 +539,7 @@ C---CALCULATE THE C-PARAMETER
C=C-3*DOT(P,I,J)**2*OS/(E(I)*E(J))
ENDDO
ENDDO
CALL GFILL1(201+ORD,C,C*WEIGHT/CBS1)
CALL GFILLSC1(201+ORD,C,C*WEIGHT)
IF (ORD.EQ.0) CS(3)=CS(3)+WEIGHT*C
C---CALCULATE THE D-PARAMETER
IF (N.EQ.4) THEN
@ -546,7 +555,7 @@ C---CALCULATE THE D-PARAMETER
ELSE
D=0
ENDIF
CALL GFILL1(202+ORD,D,D*WEIGHT/DBS)
CALL GFILLSC1(202+ORD,D,D*WEIGHT)
IF (ORD.EQ.0) CS(4)=CS(4)+WEIGHT*D
C---CALCULATE THE THRUST
IF (N.EQ.4) THEN
@ -589,7 +598,7 @@ C---CALCULATE THE THRUST
T=MAX(E(1),E(2),E(3))*ORS*2
ENDIF
TAU = 1 - T
CALL GFILL1(203+ORD,TAU,TAU*WEIGHT/TBS)
CALL GFILLSC1(203+ORD,TAU,TAU*WEIGHT)
IF (ORD.EQ.0) CS(5)=CS(5)+TAU*WEIGHT
C---CALCULATE THE Y3 VALUE (USES P SCHEME FOR NO PARTICULAR REASON)
IF (N.EQ.4) THEN
@ -630,7 +639,7 @@ C---CALCULATE THE Y3 VALUE (USES P SCHEME FOR NO PARTICULAR REASON)
ELSE
Y3=1-T
ENDIF
CALL GFILL1(204+ORD,Y3,Y3*WEIGHT/YBS)
CALL GFILLSC1(204+ORD,Y3,Y3*WEIGHT)
IF (ORD.EQ.0) CS(6)=CS(6)+Y3*WEIGHT
C---CALCULATE THE ENERGY-ENERGY CORRELATION
DO I=2,N
@ -638,7 +647,7 @@ C---CALCULATE THE ENERGY-ENERGY CORRELATION
COSANG=1-DOT(P,I,J)*OE(I)*OE(J)
IF (ABS(COSANG).GE.1) COSANG=SIGN(1-1D-12,COSANG)
BEEC=(1-COSANG**2)*2*WEIGHT*E(I)*E(J)*OS
CALL GFILL1(205+ORD,COSANG,BEEC/EBS)
CALL GFILLSC1(205+ORD,COSANG,BEEC)
COSLIM=0.1D0
IF (ORD.EQ.0) THEN
IF (ABS(COSANG).LT.COSLIM)

14
event2_merge.f
View File

@ -55,13 +55,13 @@ C-----------------------------------------------------------------------
CHARACTER*256 LINE
CHARACTER*60 TITLE
INTEGER LENTIT
INTEGER ID
INTEGER NX
INTEGER ID, NX, LH
DOUBLE PRECISION XMIN, XMAX
NX = 0
XMIN = 0
XMAX = 0
LH = 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',
@ -70,16 +70,16 @@ C --- NEED TO USE UNIT 20 TO AVOID CONFLICT WITH UNIT 10 USED BY GREAD
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)
CALL GPROP1(ID, NX, XMIN, XMAX,LH)
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)
CALL GBOOK1(200+ID,TITLE,NX,XMIN,XMAX,LH)
CALL GBOOK1(300+ID,TITLE,NX,XMIN,XMAX,LH)
CALL GBOOK1(210+ID,TITLE,NX,XMIN,XMAX,LH)
CALL GBOOK1(310+ID,TITLE,NX,XMIN,XMAX,LH)
C --- ACCUMULATE WEIGHTS
CALL GOPERA(100+ID,'*',100+ID,100+ID,1D0,1D0)

69
gbook.f
View File

@ -155,7 +155,7 @@ C Returns year, month, day as integers
RETURN
END
C*********************************************************************
SUBROUTINE GBOOK1(ID,TITLE,NX,XL,XU)
SUBROUTINE GBOOK1(ID,TITLE,NX,XL,XU,ILOG)
IMPLICIT INTEGER (I-N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NSIZE=200000,NMAX=2000)
@ -172,7 +172,19 @@ C*********************************************************************
A(IS+1)=NX
A(IS+2)=XL
A(IS+3)=XU
A(IS+4)=(XU-XL)/NX
A(IS+6)=ILOG
IF (ILOG.EQ.0) THEN
C ------ LINEAR BINNING
A(IS+4)=(XU-XL)/NX
ELSE
C ------ LOG BINNING
IF (XL.LE.0D0.OR.XU.LE.0D0) THEN
A(IS+6)=0
A(IS+4)=(XU-XL)/NX
ELSE
A(IS+4) = (LOG(XU)-LOG(XL))/NX
ENDIF
ENDIF
A(IS+5)=1
CALL GRESET(ID)
TITFX=TITLE//' '
@ -235,10 +247,42 @@ C*********************************************************************
IF(X.GE.A(IS+3)) IOX=3
A(IS+12+IOX)=A(IS+12+IOX)+W
IF(IOX.NE.2) RETURN
IX=(X-A(IS+2))/A(IS+4)
IF(A(IS+6).EQ.1) THEN
IF (X.LE.0D0) RETURN
IX=(LOG(X)-LOG(A(IS+2)))/A(IS+4)
ELSE
IX=(X-A(IS+2))/A(IS+4)
ENDIF
A(IS+19+IX)=A(IS+19+IX)+W
RETURN
END
C*********************************************************************
SUBROUTINE GFILLSC1(ID,X,W)
IMPLICIT INTEGER (I-N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NSIZE=200000,NMAX=2000)
COMMON /GBOOK/ A(NSIZE)
IF (ID.GT.NMAX) RETURN
IS=A(ID+2)+0.5
A(IS+9)=A(IS+9)+1.
IOX=2
IF(X.LT.A(IS+2)) IOX=1
IF(X.GE.A(IS+3)) IOX=3
A(IS+12+IOX)=A(IS+12+IOX)+W
IF(IOX.NE.2) RETURN
IF(A(IS+6).EQ.1) THEN
IF (X.LE.0D0) RETURN
IX=(LOG(X)-LOG(A(IS+2)))/A(IS+4)
XI = EXP((IX)*A(IS+4))
C ------ MAYBE CACHE LAST FACTOR FOR EFFICIENCY?
DX = A(IS+2)*XI*(EXP(A(IS+4))-1)
ELSE
IX=(X-A(IS+2))/A(IS+4)
DX = A(IS+4)
ENDIF
A(IS+19+IX)=A(IS+19+IX)+W/DX
RETURN
END
C*********************************************************************
SUBROUTINE GFILL2(ID,X,Y,W)
IMPLICIT INTEGER (I-N)
@ -714,7 +758,11 @@ C-----------------------------------------------------------------------
50 IF (A(IS+18+IX).GT.0.0) YMIN=MIN(YMIN,A(IS+18+IX))
DIV=.FALSE.
DO 100 IX=1, NX
X(IX)=(IX-0.5)*A(IS+4)+A(IS+2)
IF (A(IS+6).EQ.1) THEN
X(IX)=EXP(DLOG(A(IS+2))+(IX-0.5)*A(IS+4))
ELSE
X(IX)=(IX-0.5)*A(IS+4)+A(IS+2)
ENDIF
Y(IX)=A(IS+18+IX)
E(IX)=0
IF (ISERR.NE.0) E(IX)=A(ISERR+18+IX)
@ -743,6 +791,9 @@ C SET UP PAGE...
ELSE
WRITE (21,*) 'SET LIMITS X',SNGL(A(IS+2)),SNGL(A(IS+3))
ENDIF
IF(A(IS+6).EQ.1) THEN
WRITE (21,*) 'SET X LOG BINS'
ENDIF
WRITE (21,*) 'SET PATT 0.02 0.08'
WRITE (21,*) 'SET ORDER X Y DY'
NTYPE=0
@ -939,7 +990,7 @@ C-----------------------------------------------------------------------
END
C-----------------------------------------------------------------------
C --- SIMPLE UTILITY TO OBTAIN HIST PROPERTIES
SUBROUTINE GPROP1(ID, NX, XMIN, XMAX)
SUBROUTINE GPROP1(ID, NX, XMIN, XMAX, LH)
IMPLICIT INTEGER (I-N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
@ -963,6 +1014,8 @@ C --- SIMPLE UTILITY TO OBTAIN HIST PROPERTIES
NX = NXMAX
ENDIF
LH = A(IS+6)
RETURN
END
C-----------------------------------------------------------------------
@ -1000,6 +1053,8 @@ C --- PROCESS KEYS
IF(KEYS(2).EQ.'LIMITS') THEN
READ(KEYS(4),*) XMIN
READ(KEYS(5),*) XMAX
ELSE IF(KEYS(2).EQ.'X'.AND.KEYS(3).EQ.'LOG') THEN
LH = 1
ENDIF
ELSE IF (KEYS(1).EQ.'TITLE') THEN
@ -1017,9 +1072,9 @@ C --- PROCESS KEYS
ELSE IF (KEYS(1).EQ.'DATA') THEN
INDATA = .TRUE.
CALL GBOOK1(ID,TITLE,NX,XMIN,XMAX)
CALL GBOOK1(ID,TITLE,NX,XMIN,XMAX, LH)
IF(IDERR.NE.0) THEN
CALL GBOOK1(IDERR,TITLE,NX,XMIN,XMAX)
CALL GBOOK1(IDERR,TITLE,NX,XMIN,XMAX, LH)
ENDIF
C --- DATA LINES