diff --git a/event2_03.f b/event2_03.f index 604608e..58438bc 100644 --- a/event2_03.f +++ b/event2_03.f @@ -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) diff --git a/event2_merge.f b/event2_merge.f index 73761a3..908ad18 100644 --- a/event2_merge.f +++ b/event2_merge.f @@ -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) diff --git a/gbook.f b/gbook.f index 6f4676e..ae34bc8 100644 --- a/gbook.f +++ b/gbook.f @@ -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