* RUNPACK file SUBROUTINE CAIM(CARD) * * CAIM is a SUBROUTINE for the format-free input of steering * data. * CHARACTER*(*) CARD LOGICAL START CHARACTER*16 KEYWRD, DLABEL * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== CHARACTER*18 OPTE(10) DATA OPTE/'SMOOTHX ','FXPOSITIVE','ZEROLEFT','ZERORIGHT', + ' ',' ',' ',' ', + ' ',' '/ SAVE START DATA START/.TRUE./ * ... IF(START) THEN START=.FALSE. * INITIALIZATION - RESET ETC LUNP=6 DO J=1,33 DO I=1,12 JTAB(I,J)=0 END DO END DO INWT=0 INIX=0 TITLE='UNFOLDING PROGRAM' DATLUM=0.0 XMCLUM=1.0 FACLUM=1.0 DO J=0,32 VTEXT(J)=' ' END DO DO J=1,10 LOPT(J)=0.0 END DO DO J=1,5 SOPTE(J)=' ' END DO DO J=1,NHU HUW(J)=1.0 HUF(J)=1.0 HUH(J)=0.0 HUS(J)=0.0 HUA(J)=0.0 END DO BINS(1)=0.0 WLIN=0.0 NTAB=0 NV=0 NBDEF=0 XLOW=0.0 XHIG=0.0 PERC=1.0 LENGD=LENGDP LPFLAG=-3 NOPT=0 LUSFUN=0 CALL SMTEXT(' ') CALL SMTEXT('RUN - Program for Regularized UNfolding') CALL SMTEXT(' ') CALL SMTEXT('The algorithm of the program is described in the + Reports:') CALL SMTEXT('V. Blobel, Unfolding methods in high-energy + physics experiments, DESY 84-118, 1984') CALL SMTEXT('Proceedings from the 1984 CERN School of + Computing, Aiguablava, Spain 9 - 22 September 1984, + CERN 85-02') CALL SMTEXT('Copyright V. Blobel (1996)') CALL SMTEXT(' ') CALL SMTEXT('RUN - User instructions') CALL SMTEXT('=======================') END IF * ...... IF(LPFLAG.GT.0) WRITE(*,101) CARD CALL CTRTXT(NK,CARD) DO IK=1,NK * interprete text from CARD (NK = number of keywords) IDMB=0 * get first keyword from statement CALL TVAR(IK,KEYWRD,NUM,IA,IB) C WRITE(*,*) IK,KEYWRD,NUM,'ik keywrd num' IF(KEYWRD.EQ.'DATA') THEN IDMB=1 ELSE IF(KEYWRD.EQ.'MONTECARLO') THEN IDMB=2 ELSE IF(KEYWRD.EQ.'BACKGROUND') THEN IDMB=3 END IF IF(IDMB.NE.0) THEN * Data or montecarlo or background IF(NTAB.LT.32) THEN NTAB=NTAB+1 JTAB(1,NTAB)=IDMB CALL TVAR(IK+1,DLABEL,NN,IA,IB) CALL CHTOHL(DLABEL(1:4),IHL,NHL) JTAB(2,NTAB)=IHL * integrated lumi is first number ... STAB(4,NTAB)=FVAR(2,1) IF(IDMB.EQ.1) THEN * data DATLUM=DATLUM+STAB(4,NTAB) C IF(STAB(4,NTAB).NE.0.0) THEN C STAB(4,NTAB)=1.0/STAB(4,NTAB) C ELSE C STAB(4,NTAB)=1.0 C END IF C IF(FACLUM.EQ.1.0) THEN C IF(STAB(4,NTAB).NE.1.0) FACLUM=STAB(4,NTAB) C END IF STAB(9,NTAB)=1.0 ELSE IF(IDMB.EQ.2) THEN * Montecarlo IF(STAB(4,NTAB).EQ.0.0) XMCLUM=0.0 IF(STAB(4,NTAB).NE.0.0.AND.XMCLUM.NE.0.0) THEN XMCLUM=STAB(4,NTAB) END IF C IF(STAB(4,NTAB).NE.0.0) THEN C STAB(4,NTAB)=1.0/STAB(4,NTAB) C ELSE C STAB(4,NTAB)=1.0 C END IF C STAB(9,NTAB)=STAB(4,NTAB) C IF(NUM.GE.2) THEN * second number in second keyword STAB(9,NTAB)=FVAR(2,2) IF(STAB(9,NTAB).EQ.0.0) STAB(9,NTAB)=1.0 C STAB(9,NTAB)=1.0/STAB(9,NTAB) C END IF END IF CALL TVAR(IK+2,DLABEL,NN,IA,IB) IF(DLABEL.EQ.'LOOK'.AND.NN.GE.1) THEN JTAB(11,NTAB)=IVAR(IK+2,1) JTAB(12,NTAB)=0 IF(NN.GE.2) JTAB(12,NTAB)=IVAR(IK+2,2) END IF GOTO 100 END IF ELSE IF(KEYWRD.EQ.'FRLOOK') THEN INWT=IVAR(IK,1) INIX=IVAR(IK,2) ELSE IF(KEYWRD.EQ.'PRINTUNIT') THEN LUNP=IVAR(IK,1) ELSE IF(KEYWRD.EQ.'NTUPLEUNITS') THEN LUN1=IVAR(IK,1) LUN2=IVAR(IK,2) IF(LUN1.EQ.0) LUN1=51 IF(LUN2.EQ.0) LUN2=52 CALL ABUNIT(LUN1,LUN2) ELSE IF(KEYWRD.EQ.'VARIABLE') THEN * index of variable JP=IVAR(IK,1) IF(JP.GE.0.AND.JP.LE.32) THEN CALL TVAR(0,VTEXT(JP),NUM,IA,IB) IF(JP.GT.0) THEN * number of bins for variable NDIM=IVAR(1,2) IF(NDIM.GT.0.AND.NV.LE.2) THEN NV=NV+1 NDV(NV)=JP NBV(NV)=NDIM XYFL(NV)=0.0 XYFH(NV)=0.0 IF(NUM.GE.4) THEN XYFL(NV)=AMIN1(FVAR(IK,3),FVAR(IK,4)) XYFH(NV)=AMAX1(FVAR(IK,3),FVAR(IK,4)) END IF END IF END IF END IF ELSE IF(KEYWRD.EQ.'TITLE') THEN CALL TVAR(0,TITLE,NUM,IA,IB) ELSE IF(KEYWRD.EQ.'KNOTS') THEN NKN=IVAR(IK,1) ELSE IF(KEYWRD.EQ.'NRDF') THEN NDF=IVAR(IK,1) ELSE IF(KEYWRD.EQ.'XLIMITS') THEN XLOW=FVAR(IK,1) XHIG=FVAR(IK,2) ELSE IF(KEYWRD.EQ.'XBINS') THEN * NBDEF=number of values -1 NBDEF=NUM-1 DO I=1,NBDEF+1 BINS(I)=FVAR(IK,I) END DO ELSE IF(KEYWRD.EQ.'CONSTRAINT') THEN IF(NCNS.LT.10) THEN NCNS=NCNS+1 FCNS(1,NCNS)=FVAR(IK,1) FCNS(2,NCNS)=FVAR(IK,2) END IF ELSE IF(KEYWRD.EQ.'LINWEIGHT') THEN WLIN=FVAR(IK,1) ELSE IF(KEYWRD.EQ.'USEFRACT') THEN PERC=FVAR(IK,1) ELSE IF(KEYWRD.EQ.'PRINT') THEN LPFLAG=-IABS(IVAR(IK,1)) ELSE IF(KEYWRD.EQ.'UNCOMD') THEN LENGD=IVAR(IK,1) ELSE DO I=1,10 IF(KEYWRD.EQ.OPTE(I)) THEN LOPT(I)=1 NOPT=NOPT+1 SOPTE(NOPT)=OPTE(I) END IF END DO END IF END DO IF(LPFLAG.LT.0) WRITE(*,101) CARD LPFLAG=IABS(LPFLAG) 100 RETURN 101 FORMAT(4X,A) END SUBROUTINE RDCAIM(LUN) * read text from unit LUn until data end or END or ENDQ CHARACTER*80 CARD, KEYWD*16 10 READ(LUN,101,END=100) CARD IF(CARD.EQ.'END'.OR.CARD.EQ.'ENDQ') GOTO 100 CALL CAIM(CARD) CALL TVAR(2,KEYWD,NN,IA,IB) IF(NN.GE.0) GOTO 10 CALL TVAR(1,KEYWD,NN,IA,IB) IF(NN.NE.0) GOTO 10 IF(KEYWD.NE.'END'.AND.KEYWD.NE.'ENDQ') GOTO 10 100 RETURN 101 FORMAT(A) END SUBROUTINE UPRPAR * Print unfolding parameter * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== CALL SMTEXT(' ') CALL SMTEXT('Unfolding parameter') CALL SMTEXT('-------------------') WRITE(LUNP,101) NKN,NDF,XLOW,XHIG WRITE(LUNP,102) NBDEF,(BINS(I),I=1,NBDEF+1) * text fro variables DO J=1,32 IF(VTEXT(J).NE.' ') THEN WRITE(LUNP,107) J,VTEXT(J) END IF END DO DO K=1,NV WRITE(LUNP,108) NDV(K),NBV(K),XYFL(K),XYFH(K) END DO DO I=1,NCNS WRITE(LUNP,103) FCNS(1,I),FCNS(2,I) END DO WRITE(LUNP,104) WLIN WRITE(LUNP,105) (SOPTE(I),I=1,NOPT) C BLA=USFUN(0.5*(XLOW+XHIG)) IF(LUSFUN.EQ.4711) THEN * standard USFUN WRITE(LUNP,109) 'Standard default version of USFUN used' ELSE * USER version of USFUN WRITE(LUNP,109) 'USER version of USFUN used' END IF IF(PERC.EQ.1.0) GOTO 100 WRITE(LUNP,106) PERC 100 RETURN 101 FORMAT( 2 3X,'Number of knots ',I6/ 3 3X,'Number of degress of freedom ',I6/ 4 3X,'X-range (Xlow,Xhig) ',2G12.4) 102 FORMAT( 1 3X,'Number of X-bins ',I6/ 2 3X,'Bin limits ',5G12.4/ 3( 3X,' ',5G12.4)) 103 FORMAT( 1 3X,'Constraint (X,F(X)) ',2G12.4) 104 FORMAT( 1 3X,'Weight for lin. regul. contribution ',G12.4) 105 FORMAT( 1 3X,'Options ',5(2X,A8)/ 2( 3X,' ',5(2X,A8)) ) 106 FORMAT(11X,'Test! Only fraction',F5.2,' of records used') 107 FORMAT( 1 3X,'Text for variable ',I3,2X,A50) 108 FORMAT( 1 3X,'Variable for unfolding fit Index =',I3,2X, 2 I3,' bins, limits = ',2G12.4) 109 FORMAT(3X,A) END SUBROUTINE RPASS1(NAME) * * first pass of unfolding * read one event from file A * allow user to modify event data * make statistic, determine min and max values * store accepted event on file b * at end: print statistic, min and max values * * COMMON/UNCOMA/NY,WT,X,Y(100) * NAMEH = label (type hollerith) * NY = number of y-values * WT = weight of event * X = x value * Y(*) = y values * * 10 CALL RPASS1(NAME) * * NAME is label from ntuple * IF(NAME.EQ.' ') GOTO 20 * ... * GOTO 10 * 20 CONTINUE * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== COMMON/CUSFUN/SPHUFS(NHU) REAL SEWT(3),SLWT(3),SQWT(3) CHARACTER*4 NAME,LNAME,NAMEUN CHARACTER*4 PTEXT(0:3)/'????','DATA','MOCA','BACK'/ LOGICAL START EXTERNAL USFUN SAVE START,SEWT,SLWT,SQWT,LNAME,LNAMEH,IOUTX,RANA,RANN,ICOUNT,IUNK DATA START/.TRUE./ * ... IF(.NOT.START) GOTO 10 * * initialization * START=.FALSE. ICOUNT=0 IUNK=0 NAMEUN=' ' IF(LPFLAG.GE.3) THEN CALL SMTEXT(' ') CALL SMTEXT('-> Subroutine RPASS1 starting ...') C CALL SMTEXT('RUN - Start of pass 1 - '//TITLE//'.') C CALL SMTEXT('=====================') CALL SMTEXT('The ntuples from unit A are read and analysed. + The user may modify the data within the MAIN program. ') CALL SMTEXT('A statistic on the events is collected and + printed. ') CALL SMTEXT('The min and max values of the variables are + determined from the DATA and MOCA events.') END IF BLA=USFUN(0.0) IF(LPFLAG.GE.2) THEN * print steering parameters CALL UPRPAR END IF ITAB =0 * number of bins for moca histogram NHA=4*(NKN-3) * reset experimental limits DO J=0,32 JXY(J) =0 XYMIN(J)=0.0 XYMAX(J)=0.0 END DO IOUTX=0 LNAME=' ' CALL CHTOHL(LNAME,LNAMEH,NHL) * data for pseudo random numbers RANA=SQRT(13.0) RANA=AMOD(RANA,1.0) RANN=0.0 * go to read first record GOTO 50 * * check user defined bin limits for fit variables * 10 DO K=1,NV IY=NDV(K) IF(XYFL(K).LT.XYFH(K)) THEN IF(Y(IY).LT.XYFL(K).OR.Y(IY).GT.XYFH(K)) WT=-ABS(WT) END IF END DO * check monte carlo X IF(IDB.EQ.2.AND.WT.GT.0.0) THEN C IF(Y(0).LT.XLOW.OR.Y(0).GT.XHIG) THEN IF(X .LT.XLOW.OR.X .GT.XHIG) THEN * outside! IOUTX=IOUTX+1 IF(IOUTX.LE.5) WRITE(*,*) X,' = X outside limits ',XLOW,XHIG END IF END IF C IF(ICOUNT.LE.100) WRITE(*,*) IDMB,WT,' is IDMB and WT' * * count event in table * IF(ITAB.NE.0) THEN IF(NAMEH.EQ.JTAB(2,ITAB)) GOTO 20 END IF DO ITAB=1,NTAB IF(NAMEH.EQ.JTAB(2,ITAB)) GOTO 20 END DO * add new row to table ITAB=NTAB+1 JTAB(2,ITAB)=NAMEH 20 IF(LPFLAG.GE.3) THEN IF(STAB(5,ITAB).LT.5.0) THEN IF(ITAB.EQ.1.AND.STAB(5,ITAB).EQ.0.0) WRITE(*,109) * print first 5 events from each type WRITE(*,108) NAMEH,NY,WT,X,(Y(I),I=1,NY) END IF END IF * add to table columns STAB(5,ITAB)=STAB(5,ITAB)+1.0 STAB(6,ITAB)=STAB(6,ITAB)+ABS(WT) IF(WT.GT.0.0) THEN STAB(7,ITAB)=STAB(7,ITAB)+WT STAB(8,ITAB)=STAB(8,ITAB)+WT*WT END IF * define flag IDMB=JTAB(1,ITAB) IF(IDMB.LE.0.OR.IDMB.GT.3) GOTO 50 * * make moca histogram of input x-distribution * (including ntuples with negative weight) IF(IDMB.EQ.2) THEN AWT=ABS(WT)*STAB(9,ITAB)/(XHIG-XLOW)*FLOAT(NHA) AWT=ABS(WT)/(XHIG-XLOW)*FLOAT(NHA) J=1.0+(X-XLOW)/(XHIG-XLOW)*FLOAT(NHA) IF(J.GE.1.AND.J.LE.NHA) THEN HUH(J)=HUH(J)+AWT HUS(J)=HUS(J)+AWT*AWT IF(WT.GT.0.0) HUA(J)=HUA(J)+AWT END IF END IF * * determine min and max values from data and Moca * IF(WT.GT.0.0.AND.(IDMB.EQ.1.OR.IDMB.EQ.2)) THEN IF(IDMB.EQ.1) JY=1 IF(IDMB.EQ.2) JY=0 DO J=JY,NY IF(JXY(J).EQ.0) THEN JXY(J)=1 XYMIN(J)=XY(J) XYMAX(J)=XY(J) ELSE XYMIN(J)=MIN(XYMIN(J),XY(J)) XYMAX(J)=MAX(XYMAX(J),XY(J)) END IF END DO END IF * * write to intermediate file * IF(ITAB.LE.NTAB.AND.WT.GT.0.0) THEN MY=NY CALL BSTORE(ITAB,REC,2+MY) END IF * * read next event from file ---------------------------------------- * 50 CALL AFETCH(NAMEH,REC,NR) ICOUNT=ICOUNT+1 C IF(MOD(ICOUNT,100).EQ.1) C + WRITE(*,199) ICOUNT,NAMEH,(REC(JJ),JJ=1,NR) 199 FORMAT(1X,I7,1X,A4,1X,4F8.3) IF(NR.LT.0) GOTO 60 * event (ntuple) read in IF(PERC.NE.1.0) THEN * eventually skip in case of reduced statistic RANN=AMOD(RANN+RANA,1.0) IF(RANN.GE.PERC) GOTO 50 END IF * NY=NR-2 IF(NAMEH.NE.LNAMEH) THEN * new name - convert to text string LNAMEH=NAMEH CALL HLTOCH(LNAMEH,1,LNAME,NCH) C WRITE(*,*) 'NEW NAME IS ',LNAME END IF NAME=LNAME C IF(ICOUNT.LE.100) WRITE(*,*) NAME,' is name ' * find corresponding table row IF(ITAB.NE.0) THEN IF(NAMEH.EQ.JTAB(2,ITAB)) GOTO 53 END IF DO ITAB=1,NTAB IF(NAMEH.EQ.JTAB(2,ITAB)) GOTO 53 END DO * skip event with unknown label IF(NAMEUN.NE.NAME) THEN IUNK=IUNK+1 IF(IUNK.LE.50) THEN WRITE(*,*) IUNK,'.th record with unknown label:',NAME NAMEUN=NAME END IF END IF GOTO 50 * table row found - define IDMB flag 53 IDMB=JTAB(1,ITAB) C IF(ICOUNT.LE.100) WRITE(*,*) ' Assigned IDMB =', IDMB GOTO 100 * * end-of-data = end of pass 1 -------------------------------------- * 60 NAME=' ' * CALL BUNIT(52) CALL BSTORE(ID,IDUM,-1) * * complete table * DO ITAB=1,NTAB * normalization factor for background IF(JTAB(1,ITAB).EQ.3) THEN IF(STAB(6,ITAB).NE.0.0) STAB( 9,ITAB)=STAB(4,ITAB)/STAB(6,ITAB) STAB(9,ITAB)=PERC*STAB(9,ITAB) END IF STAB(10,ITAB)=STAB(7,ITAB)*STAB(9,ITAB) END DO * * printout of table * IF(LPFLAG.GE.2) THEN CALL SMTEXT(' ') CALL SMTEXT('Event statistic') CALL SMTEXT('---------------') WRITE(LUNP,101) DO ITAB=1,NTAB IDMB=JTAB(1,ITAB) EFF=0.0 DIV=STAB(8,ITAB)*STAB(9,ITAB) IF(DIV.NE.0.0) EFF=STAB(7,ITAB)/DIV WRITE(LUNP,107) PTEXT(IDMB),JTAB(2,ITAB), + (STAB(I,ITAB),I=4,10),EFF END DO END IF DO I=1,3 SLWT(I)=0.0 SQWT(I)=0.0 SEWT(I)=0.0 END DO DO ITAB=1,NTAB I=JTAB(1,ITAB) IF(I.NE.0) THEN SLWT(I)=SLWT(I)+STAB( 7,ITAB)*STAB( 9,ITAB) SQWT(I)=SQWT(I)+STAB( 8,ITAB)*STAB( 9,ITAB)**2 END IF END DO IF(LPFLAG.GE.2) THEN WRITE(LUNP,103) DO I=1,3 IF(SQWT(I).NE.0.0) THEN EFF=SLWT(I)/SQWT(I) WRITE(LUNP,104) PTEXT(I),SLWT(I),EFF END IF END DO END IF IF(LPFLAG.GE.3) THEN CALL SMTEXT(' ') CALL SMTEXT('Explanation of table:') CALL SMTEXT('ARG = constant, specified in the user instruction') CALL SMTEXT('Nr of ev = total number of weighted events + = number of events') CALL SMTEXT('Sum |W| = sum of all absolute values of weights') CALL SMTEXT('Sum W = sum of all positive weights') CALL SMTEXT('Sum W**2 = sum of square of weights for positive + weights') CALL SMTEXT('Factor = 1.0 for DATA') CALL SMTEXT('Factor = Arg for MONTE CARLO') CALL SMTEXT('Factor = Ratio of ARG/SUM |W| for BACKGROUND. This + factor is applied to the background events') CALL SMTEXT('Fin.nr = final event number for unfolding. In + case of background it is the event number (eventually + after cuts), which is used in the fit-distribution') CALL SMTEXT('Ev-eff = event efficiency') END IF * check bin limits DO J=0,32 IF(JXY(J).NE.0) THEN IF(XYMIN(J).NE.XYMAX(J)) THEN * rounded limits for final histograms CALL ROUNDB(XYMIN(J),XYMAX(J),120,XYL(J),XYH(J)) * correct limits slightly C XYMIN(J)=AMAX1(XYMIN(J)*1.0001-0.0001*XYMAX(J),XYL(J)) C XYMAX(J)=AMIN1(XYMAX(J)*1.0001-0.0001*XYMIN(J),XYH(J)) ELSE JXY(J)=0 END IF END IF END DO * define missing limits for fit variables DO K=1,NV J=NDV(K) IF(XYFL(K).EQ.XYFH(K)) THEN XYFL(K)=XYMIN(J)*1.0001-0.0001*XYMAX(J) IF(XYFL(K).LT.0.0.AND.XYMIN(J).GE.0.0) XYFL(K)=0.0 XYFH(K)=XYMAX(J)*1.0001-0.0001*XYMIN(J) IF(XYFH(K).GT.0.0.AND.XYMAX(J).LE.0.0) XYFH(K)=0.0 END IF END DO * * print min max values * IF(LPFLAG.GE.2) THEN CALL SMTEXT(' ') CALL SMTEXT('Range of variables, exact and rounded') CALL SMTEXT('-------------------------------------') WRITE(LUNP,102) DO J=0,32 IF(JXY(J).NE.0) THEN IF(XYMIN(J).NE.XYMAX(J)) THEN WRITE(LUNP,106) J,XYMIN(J),XYMAX(J),XYL(J),XYH(J) END IF END IF END DO WRITE(LUNP,111) DO K=1,NV J=NDV(K) WRITE(LUNP,112) K,J,NDV(K),XYFL(K),XYFH(K) END DO END IF IF(LPFLAG.GE.2) THEN WRITE(LUNP,105) XLOW,XHIG IF(IOUTX.NE.0) THEN WRITE(LUNP,110) IOUTX END IF END IF * copy histogram of x values to HUF ... IHIS=0 DO I=1,NHA IF(HUS(I).NE.0.0) IHIS=1 HUF(I)=HUH(I) END DO * for option SMOOTHX smooth histogram before the fit ... IF(IHIS.EQ.1) THEN IF(LOPT(1).EQ.1) THEN CALL SPLFT(HUF,NHA) ELSE * fit of MOCA histogram -> spline coefficients CALL ESPQ(HUF,HUS,NHA,XLOW,XHIG,NKN,HUF,ASPMC) END IF * spline function definition from smoothed data CALL USPC(HUF,SPHUFS,NHA) END IF * * Normalization consideration =========================== * WRITE(LUNP,'(1X)') CALL SMTEXT('================================================') CALL SMTEXT('....NOTES on the normalization of the result....') *2345678901234567890123456789012346789012345678901234567890123456789012 WRITE(LUNP,'(1X)') CALL SMTEXT('The calculations are based on the identity') WRITE(LUNP,118) ' N = Lumi(integrated) * Int( f(x) dx )' WRITE(LUNP,118) 'or ' WRITE(LUNP,118) ' dN/dx = Lumi(integrated) * f(x) ' WRITE(LUNP,118) * GAMMA=1.0 * integrate function USFUN * average value in X region SAVF=AVFUN(XLOW,XHIG,USFUN) * integral over X region SINT=SAVF*(XHIG-XLOW) * determine number of Monta Carlo events (sum of weights) SMCI=0.0 SMCT=0.0 DO ITAB=1,NTAB IDMB=JTAB(1,ITAB) IF(IDMB.EQ.2) THEN * MOCA only SMCI=SMCI+STAB(7,ITAB) SMCT=SMCT+STAB(6,ITAB) END IF END DO * if Lexp not given - 1.0 assumed IF(DATLUM.EQ.0.0) DATLUM=1.0 * test function USFUN C BLA=USFUN(0.5*(XLOW+XHIG)) FLUMUS=XMCLUM IF(LUSFUN.EQ.4711) THEN * NO user supplied function IF(XMCLUM.EQ.0.0) THEN CALL SMTEXT('This is case 1:') CALL SMTEXT('The function XFUNCT is the default version, ') CALL SMTEXT('which is defined by the actual ') CALL SMTEXT('x-distribution of all MC events, ') CALL SMTEXT('including the rejected events with ') CALL SMTEXT('events with a negative weight. ') CALL SMTEXT('Note that these rejected events have to be ') CALL SMTEXT('included in the MC n-tuple file.' ) CALL SMTEXT('The default function XFUNCT is normalized ') CALL SMTEXT('to an integral of one.') * normalization of the function GAMMA=1.0/SINT DO I=1,NKN ASPMC(I)=GAMMA*ASPMC(I) END DO SINT=SINT*GAMMA SAVF=SAVF*GAMMA CALL SMTEXT('The MC Luminosity is not specified! ') CALL SMTEXT('The value is set to the sum of ') CALL SMTEXT('|weight|s. => see below.') FLUMUS=SMCT FACLUM=SMCT/DATLUM DATLUM=1.0 ELSE CALL SMTEXT('This is case 2:') CALL SMTEXT('The function XFUNCT is the default version, ') CALL SMTEXT('which is defined by the actual ') CALL SMTEXT('x-distribution of all MC events, ') CALL SMTEXT('including the rejected events with ') CALL SMTEXT('events with a negative weight. ') CALL SMTEXT('Note that these rejected events have to be ') CALL SMTEXT('included in the MC n-tuple file.' ) CALL SMTEXT('The MC Luminosity is specified.') CALL SMTEXT('Therefore the normalization of the function ') CALL SMTEXT('XFUNCT is determined to be consistent with ') CALL SMTEXT('the given MC Luminosity and the sum of the ') CALL SMTEXT('|weight|s => see below.') GAMMA=SMCT/(XMCLUM*SINT) SINT=SINT*GAMMA SAVF=SAVF*GAMMA FACLUM=1.0/DATLUM DATLUM=1.0 END IF *2345678901234567890123456789012346789012345678901234567890123456789012 ELSE * function is defined by user IF(XMCLUM.EQ.0.0) THEN CALL SMTEXT('This is case 3: ') CALL SMTEXT('The function XFUNCT is supplied by the user. ') CALL SMTEXT('The MC Luminosity is not specified and ') CALL SMTEXT('therefore calculated as the sum of all ') CALL SMTEXT('MC |weight|s, divided by the integral of ') CALL SMTEXT('the user function XFUNCT => see below') FLUMUS=SMCT/SINT FACLUM=SMCT/SINT/DATLUM DATLUM=1.0 ELSE CALL SMTEXT('This is case 4:') CALL SMTEXT('The function XFUNCT is supplied by the user.') CALL SMTEXT('The MC Luminosity is specified and thus ') CALL SMTEXT('everything is known. As a consistency check ') CALL SMTEXT('the overall acceptance can be calculated as ') CALL SMTEXT('the sum of all MC weights, divided by the ') CALL SMTEXT('product of MC Luminosity and the integral ') CALL SMTEXT('of the user function XFUNCT.') WRITE(LUNP,117) ' Mean acceptance= ', + SMCI/(FLUMUS*SINT) FACLUM=1.0 END IF END IF IF(DATLUM.NE.1.0) THEN FACLUM=FACLUM/DATLUM DATLUM=1.0 END IF * print-out WRITE(LUNP,114) SAVF,SINT WRITE(LUNP,115) SMCI,SMCT,SMCI/SMCT * IF(XMCLUM.NE.0.0) THEN * MC lumi specified - scale it SUML=0.0 DO ITAB=1,NTAB IDMB=JTAB(1,ITAB) IF(IDMB.EQ.2) THEN * MOCA only SUML=SUML+STAB(4,ITAB) END IF END DO FACLUM=FACLUM*SUML DATLUM=DATLUM*SUML END IF CALL SMTEXT(' ') CALL SMTEXT('================================================') CALL SMTEXT('end of pass 1 ') CALL SMTEXT('-- Subroutine RPASS1 ending.') 100 RETURN 101 FORMAT('0',9X,' Name Arg Nr of ev Sum |W 1| Sum W Sum W**2 Factor FIN.NR ev-ef 2f'/) 102 FORMAT('0',14X,'J',2(5X,'Min',10X,'Max',5X)) 103 FORMAT('0',49X,'Fin.nr ev-eff'/) 104 FORMAT(30X,A4,14X,G12.4,F10.3) 105 FORMAT(7X,'Xlow Xhig',2G13.5) 106 FORMAT(12X,I4,4G13.5) 107 FORMAT(2X,A4,8X,3X,A4,6X,7G13.5,F10.3) 108 FORMAT(1X,A4,I3,5G12.4/(8X,5G12.4)) 109 FORMAT(' '/' Content of first five records for each type'/ + ' Name NY WT X Y... '/ + ' ---- -- ------ ------ -------------- ') 110 FORMAT(1X,'Warning!',I6,' MC ntuple outside limits!') 111 FORMAT(1X,'Variables for unfolding fit'/ + 3X,'# j bins bin_limits..............') 112 FORMAT(1X,I3,I3,I8,2G12.4) C 113 FORMAT(1X,I3,'.th case of unknown name, name=',A4) 114 FORMAT(1X,' Properties of function USFUN(X)'/ + ' average value in XLOW...XHIG:',G12.4/ + ' integral XLOW...XHIG:',G12.4,' (int(fun))') 115 FORMAT(' Sum of MC weights, accepted events:',G12.4/ + ' all events:',G12.4,' (Sum_o_w)'/ + ' -> acceptance:',F08.4) 116 FORMAT(' MC-Luminosity:',G12.4/ + ' Data-Luminosity:',G12.4) 117 FORMAT(1X,A,G12.4) 118 FORMAT(10X,A) END SUBROUTINE RPASS2 * * Second pass of unfolding program * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== COMMON/RUNRES/XW(527),RSP(2082) * ------------------------------------------------------------------ * Internal common for matrices, vectors etc for up to 63 knots * and up to 30 data points PARAMETER (NKMAX=63,NDAMAX=30) COMMON/UINTER/UM(NKMAX*NKMAX),HE((NKMAX+NKMAX*NKMAX)/2), + GR(NKMAX),XX(NKMAX),AM(NKMAX),DE(NKMAX), + SCV(NKMAX), + SCM(NKMAX*NKMAX), + FUN,SF,SD,TAU * M=NKN-1 * UM = tranformation matrix * HE = hessian, replaced by covariance matrix * GR = gradient, replaced by step vector * XX = coefficient * AM = amplitudes M * DE = eigenvalues M * SC = scratch storage * ------------------------------------------------------------------ SAVE /RUNRES/,/UINTER/ REAL RH(11,2),SH(2),S(4),RSREL(2),RSNUM(2) COMMON/UNCOMD/H(LENGDP) INTEGER NH(LENGDP) EQUIVALENCE (H(1),NH(1)) CHARACTER*1 XMK * ... IF(LPFLAG.GE.3) THEN CALL SMTEXT(' ') CALL SMTEXT('-> Subroutine RPASS2 starting pass 2.1 ...') C CALL SMTEXT('RUN - Start of pass 2.1 - '//TITLE//'.') C CALL SMTEXT('=======================') END IF LENGD=LENGDP IPRINT=0 * * prepare weights for MC ... DO I=1,NHA HUW(I)=1.0 IF(HUH(I).NE.0.0.AND.HUF(I).NE.0.0) THEN HUW(I)=HUF(I)/HUH(I) END IF END DO * ... and fill array with user function DO I=1,NHA X=XLOW+(FLOAT(I)-0.5)*(XHIG-XLOW)/FLOAT(NHA) HUP(I)=USFUN(X) END DO * * determine total number of bins => = NPROD * NPROD=1 DO J=1,NV IF(NBV(J).NE.0) NPROD=NPROD*NBV(J) END DO IF(LPFLAG.GE.2) THEN WRITE(LUNP,104) LENGD,4*NPROD END IF * check length of common IF(4*NPROD.GT.LENGD) GOTO 90 * * double loop to fill histograms * first loop: fill data and moca histograms with all NPROD bins * * * second loop: fill all non-zero histogram bins * * number of words first_loop second_loop * NPROD * * * * ------------------------------------------------------------------ DO 50 IREP=1,2 * clear histogram space DO I=(IREP-1)*NPROD+1,LENGD H(I)=0.0 END DO * * input of next B-ntuple - define IDMB flag 11 CALL BFETCH(ITAB,REC,NR) IF(NR.LT.0) GOTO 45 NY=NR-2 IDMB=JTAB(1,ITAB) * * determine linear index "IN" of this event for a NV-dim. distrib. * IN=1 NN=1 DO K=1,NV IY=NDV(K) L=1.0+(Y(IY)-XYFL(K))/(XYFH(K)-XYFL(K))*FLOAT(NBV(K)) * entries outside the bin limits are ignored IF(L.LT.1.OR.L.GT.NBV(K)) GOTO 11 IN=IN+(L-1)*NN NN=NN*NBV(K) END DO * modify weight WT for Monte Carlo IF(IDMB.EQ.2) THEN L=1.0+(X-XLOW)/(XHIG-XLOW)*120.0 IF(L.GE.1.AND.L.LE.120) WT=WT*HUW(L) WT=WT*STAB(9,ITAB)*DATLUM IF(STAB(4,ITAB).NE.0.0) WT=WT/STAB(4,ITAB) END IF * * filling bins .... * IF(IREP.EQ.1) THEN * in first loop (IREP = 1) * fill complete nf-dimensional histograms * total space requirement is 5*nprod * 0 indices (reserved space) * 1 data histogram * 2 moca histogram (w) * 3 moca histogram (w**2) IF(IDMB.EQ.1) THEN * data in part 1 H(IN+1*NPROD)=H(IN+1*NPROD)+WT ELSE IF(IDMB.EQ.2) THEN * moca as wt, wt* in part 2 and 3 H(IN+2*NPROD)=H(IN+2*NPROD)+WT H(IN+3*NPROD)=H(IN+3*NPROD)+WT*WT END IF ELSE * --------------------------------------------------------------- * Second loop (IREP = 2) * compact histograms (only NBN non-zero bins) * space requirement NPROD+NBN*(NKN+2) * 0 indices (NPROD words) * 1 data histogram (NBN words) * 2 background histogram (NBN words) * 3... MC histograms (NBN*NKN words) for all knots * translate bin index IN into new (compact) bin index IN IN=NH(IN) IF(IN.EQ.0) GOTO 11 * ... but this should never happen for data or moca IF(IDMB.EQ.1) THEN * data H(IN+NPROD)=H(IN+NPROD)+WT ELSE IF(IDMB.EQ.2) THEN * moca - find four weights for four B-splines CALL ESPB(X,NKN,XLOW,XHIG,JNULL,S) II=NPROD+NBN*2+(IN-1)*NKN+JNULL DO J=1,4 H(II+J)=H(II+J)+WT*S(J) END DO ELSE IF(IDMB.EQ.3) THEN * background - reweight WT=WT*STAB(9,ITAB) H(IN+NPROD+NBN)=H(IN+NPROD+NBN)+WT END IF END IF GOTO 11 * * end-of-data, determine index-table for IREP=1 -------------------- * 45 IF(IREP.NE.1) GOTO 50 * reset counter DO J=1,2 RSREL(J)=0.0 RSNUM(J)=0.0 DO I=1,11 RH(I,J)=0.0 END DO END DO * loop on all bins and check content NBN=0 DO I=1,NPROD IF(H(NPROD+I).EQ.0.0.AND.H(2*NPROD+I).EQ.0.0) THEN * data and moca bin empty - pointer set to zero NH(I)=0 ELSE * non-empty bin - pointer defined NBN=NBN+1 NH(I)=NBN * bin index for relative error of data ... J=11 IF(H(NPROD+I).NE.0.0) THEN * relative error of data RELER=1.0/SQRT(ABS(H(NPROD+I))) J=1.0+10.0*RELER IF(J.GT.10) J=10 END IF * ... and update of table RH(J,1)=RH(J,1)+1.0 RELER=AMIN1(1.0,RELER) IF(J.NE.11) THEN RSREL(1)=RSREL(1)+RELER RSNUM(1)=RSNUM(1)+1.0 END IF * bin index for relative error of moca ... J=11 IF(H(2*NPROD+I).NE.0.0) THEN * relative error of MOCA RELER=SQRT(H(3*NPROD+I))/ABS(H(2*NPROD+I)) J=1.0+10.0*RELER IF(J.GT.10) J=10 END IF * ... and update for table RH(J,2)=RH(J,2)+1.0 RELER=AMIN1(1.0,RELER) IF(J.NE.11) THEN RSREL(2)=RSREL(2)+RELER RSNUM(2)=RSNUM(2)+1.0 END IF END IF END DO * * printout of table on statistical accuracy * IF(LPFLAG.GE.2) THEN WRITE(LUNP,101) NPROD+NBN*(NKN+2),NPROD,NBN END IF IF(NPROD+NBN*(NKN+2).GT.LENGD) GOTO 90 * sums for printout DO J=1,2 SH(J)=0.0 DO I=1,11 SH(J)=SH(J)+RH(I,J) END DO END DO IF(LPFLAG.GE.3) THEN CALL SMTEXT(' ') CALL SMTEXT('Statistical accuracy of bins in fit-distribution') CALL SMTEXT('------------------------------------------------') WRITE(LUNP,109) DO I=1,10 WRITE(LUNP,102) 10*I-10,10*I, + RH(I,1),100.0*RH(I,1)/SH(1), + RH(I,2),100.0*RH(I,2)/SH(2) END DO WRITE(6,103) RH(11,1),100.0*RH(11,1)/SH(1), + RH(11,2),100.0*RH(11,2)/SH(2) RELDAT=0.0 RELMOC=0.0 IF(RSNUM(1).NE.0.0) RELDAT=RSREL(1)/RSNUM(1) IF(RSNUM(2).NE.0.0) RELMOC=RSREL(2)/RSNUM(2) WRITE(LUNP,105) 100.0*RELDAT,100.0*RELMOC C CALL SMTEXT('PAGE') C CALL SMTEXT('RUN - Start of pass 2.2 - '//TITLE//'.') CALL SMTEXT('-> Subroutine RPASS2 starting pass 2.2 ...') END IF 50 IPRINT=0 * * unfolding fit * CALL RUNFIT * * * ... IF(LPFLAG.GE.2) THEN CALL RSOLVP(1) WRITE(6,111) SNDF=0.0 SUMB=0.0 DO I=1,NKN-1 IF(I.GT.2.AND.DE(I).LE.0.0) GOTO 51 DFAC=1.0/(1.0+TAU*DE(I)) AMP2=AM(I)*AM(I) AMPD=AM(I)*DFAC SNDF=SNDF+DFAC TD=TAU*DE(I) SSDF=AMPD**2*(TD*(2.0+TD))*3.0/(3.0+TAU*DE(I)) SUMB=SUMB+SSDF XMK=' ' IF(AMP2.LT.9.00) XMK='X' WRITE(6,112) I,DE(I),AMP2,XMK,AM(I),AMPD,DFAC,SNDF,SSDF,SUMB END DO 51 WRITE(6,115) TAUIN=0.0 IF(TAU.GT.0.0) TAUIN=1.0/TAU WRITE(6,118) TAU,TAUIN END IF * calculate binned result NB=0 CALL RURSLT(NB,XW(2),LPFLAG,XW) * IF(NBDEF.NE.0) THEN CALL RURSLT(NBDEF,BINS,LPFLAG,XW) END IF IF(LPFLAG.GE.3) THEN CALL SMTEXT('-- Subroutine RPASS2 ending.') C CALL SMTEXT('End of pass 2') END IF GOTO 100 * error_exit_error_exit_ error_exit_ error_exit_ error_exit_ error_e 90 IF(LPFLAG.GE.1) THEN CALL SMTEXT('Stop - COMMON/UNCOMD/ too short for distribution. + Either reduce number of bins for fit-distributions or increase + length of COMMON') END IF STOP 100 RETURN 101 FORMAT(' Nr of words required in pass 2.2 ',I7/ + ' Nr of bins in distribution ',I7/ + ' Nr of non-empty bins ',I7) 102 FORMAT(9X,I2,' -',I3,10X,2(F7.0,F12.2,9X)) 103 FORMAT(12X,'empty bins ',2(F7.0,F12.2,9X)) 104 FORMAT('0Lenght of COMMON/UNCOMD/ ',I7/ + ' Nr of words required in pass 2.1 ',I7) 105 FORMAT('0Mean relative error in percent', + 5X,F7.2,20X,F7.2) 109 FORMAT('0',30X,' DATA MOCA'/ + ' REL. ERROR IN PERCENT ',2(' NR OF BINS PERCENT OF BINS')) 111 FORMAT(//'0 I EIGENVALUE AMPL**2 AMPL 1 D-AMPL FAC SUM NOISE SUM'/) 112 FORMAT(1X,I5,G17.5,F15.4,1X,A1,2X,2F12.4,8X,2F8.4,8X,2F8.4) 115 FORMAT(87X,'NDF',17X,7H'NOISE') 118 FORMAT(' TAU =',G15.5,' 1/TAU =',G11.4/) END SUBROUTINE RPASS3(NAME) * * final pass of unfolding program / new version * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== COMMON/RUNRES/XW(527),RSP(2082) * ------------------------------------------------------------------ * Internal common for matrices, vectors etc for up to 63 knots * and up to 30 data points PARAMETER (NKMAX=63,NDAMAX=30) COMMON/UINTER/UM(NKMAX*NKMAX),HE((NKMAX+NKMAX*NKMAX)/2), + GR(NKMAX),XX(NKMAX),AM(NKMAX),DE(NKMAX), + SCV(NKMAX), + SCM(NKMAX*NKMAX), + FUN,SF,SD,TAU * M=NKN-1 * UM = tranformation matrix * HE = hessian, replaced by covariance matrix * GR = gradient, replaced by step vector * XX = coefficient * AM = amplitudes M * DE = eigenvalues M * SC = scratch storage * ------------------------------------------------------------------ SAVE /RUNRES/,/UINTER/ C REAL BX(2),BY(2) LOGICAL START CHARACTER*4 NAME,LNAME SAVE START DATA START/.TRUE./ * ... IF(.NOT.START) GOTO 10 START=.FALSE. CALL SMTEXT('-> Subroutine RPASS3 starting ...') C CALL SMTEXT('Start of pass 3') C CALL SMTEXT('===============') * reset histogram buffer DO K=1,32 DO J=1,3 DO I=1,120 HIST(I,J,K)=0.0 END DO END DO END DO LNAME=' ' * convert character string to hollerith CALL CHTOHL(LNAME,LNAMEH,NHL) * fetch next record ------------------------------------------------ 10 CALL BFETCH(ITAB,REC,NR) IF(NR.LT.0) GOTO 40 NY=NR-2 * set flag from identifier IDMB=JTAB(1,ITAB) * reweighting IF(IDMB.EQ.2) THEN * monte carlo - modify weight by unfolding fit result L=1.0+(X-XLOW)/(XHIG-XLOW)*FLOAT(NHA) IF(L.GE.1.AND.L.LE.NHA) WT=WT*HUW(L) WT=WT*STAB(9,ITAB)*DATLUM * weight from unfolding fit NISP=RSP(1)+0.5 WT=WT*ESPF(X,RSP(4),NISP,RSP(2),RSP(3)) IF(STAB(4,ITAB).NE.0.0) WT=WT/STAB(4,ITAB) ELSE IF(IDMB.EQ.3) THEN * background WT=WT*STAB(9,ITAB) END IF * loop on all measured quantities DO J=1,NY IF(JXY(J).NE.0) THEN C IF(XYL(J).NE.XYH(J)) THEN * bin number ... L=1.0+(Y(J)-XYL(J))/(XYH(J)-XYL(J))*120.0 IF(L.GE.1.AND.L.LE.120) THEN * ... and add to histogram HIST(L,IDMB,J)=HIST(L,IDMB,J)+WT END IF END IF END DO * reconstruct name for outside plots NAMEH=JTAB(2,ITAB) IF(NAMEH.NE.LNAMEH) THEN LNAMEH=NAMEH * convert hollerith to character CALL HLTOCH(LNAMEH,1,LNAME,NCH) END IF NAME=LNAME GOTO 100 * ... after end-of-data -------------------------------------------- 40 DO J=1,32 IF(JXY(J).NE.0) THEN * add background[3] to moca[2] DO I=1,120 HIST(I,2,J)=HIST(I,2,J)+HIST(I,3,J) END DO END IF END DO NAME=' ' CALL SMTEXT('-- Subroutine RPASS3 ending.') 100 RETURN END SUBROUTINE RUNFIT * * new '95 version * * input data * XLOW, XHIG = lower and upper limit of x * NKN = number of knots * NDF = requested number of data points * = 0 automatic * > 0 fixed * < 0 use abs(ndf) for starting iterations * NCNS = number of constraints * FCNS(2,10) = pairs x/function(x) * ILEFT/IRIGHT * WLIN = relative weight * LUNP * NBDEF = number of bins * BINS(31) = bin limits * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== COMMON/RUNRES/XW(527),RSP(2082) * ------------------------------------------------------------------ * Internal common for matrices, vectors etc for up to 63 knots * and up to 30 data points PARAMETER (NKMAX=63,NDAMAX=30) COMMON/UINTER/UM(NKMAX*NKMAX),HE((NKMAX+NKMAX*NKMAX)/2), + GR(NKMAX),XX(NKMAX),AM(NKMAX),DE(NKMAX), + SCV(NKMAX), + SCM(NKMAX*NKMAX), + FUN,SF,SD,TAU * M=NKN-1 * UM = tranformation matrix * HE = hessian, replaced by covariance matrix * GR = gradient, replaced by step vector * XX = coefficient * AM = amplitudes M * DE = eigenvalues M * SC = scratch storage * ------------------------------------------------------------------ SAVE /RUNRES/,/UINTER/ REAL SC(4000) * common for MINST-package COMMON/CMINST/DFUNST,DPARST,SECDST,IPRNST,XLIMST,NEVAST SAVE /CMINST/ * ... * defaults and user options ILEFT =1 IRIGHT=1 IF(LOPT(3).NE.0) ILEFT =0 IF(LOPT(4).NE.0) IRIGHT=0 IFXPOS=1 IF(LOPT(2).NE.0) IFXPOS=0 LG=0 NDATA=NDF * M=NKN-1 * UM = TRANFORMATION MATRIX * HE = HESSIAN, REPLACED BY COVARIANCE MATRIX * GR = GRADIENT, REPLACED BY STEP VECTOR * XX = COEFFICIENT * AM = AMPLITUDES M * DE = EIGENVALUES M * SC = SCRATCH STORAGE * commentcommentcommentcommentcommentcommentcommentcommentcomment * print report CALL SMTEXT(' ') CALL SMTEXT('Unfolding fit starting') CALL SMTEXT('The fit uses the one-, two- or three-dimensional + distribution of events. The Monte-Carlo distribution plus + background distribution is fitted to the data distribution. ') CALL SMTEXT('In the first iteration the least squares method + is used, to determine initial values of the fit + parameter. ') CALL SMTEXT('The following iterations are based on the + poisson statistic of data events and use the maximum + likelihood method.') CALL SMTEXT('The function value is a chi square for the + first iteration and similar to a chi suare for the + following iterations') * commentcommentcommentcommentcommentcommentcommentcommentcomment * ... * reset coefficients WRITE(*,*) 'NKN= ',NKN write(*,*) 'NCNS = ',NCNS DO I=1,NKN XX(I)=0.0 END DO IT=0 * define number of knots and limits ... 1 CALL RSOLVE(NKN,XLOW,XHIG,ILEFT,IRIGHT,WLIN) * specify constraints IF(IFXPOS.EQ.0) THEN CALL RSOLFC(0.0,0.0) END IF * specify all contraint value/function pairs DO L=1,NCNS CALL RSOLFC(FCNS(1,L),FCNS(2,L)/USFUN(FCNS(1,L))/FACLUM) write(*,*) ' constraint ',L, + FCNS(1,l),fcns(2,l)/usfun(fcns(1,l))/FACLUM END DO CALL RSOLVP(0) * start next iteration 10 IT=IT+1 * calculate FUN,GR(.),HE(.) (function, gradient, Hessian) CALL RUNFUN(IT) * define regularization matrix SC ... CALL RSOLRM(SC) C IF(MARKE.NE.4711) CALL SMPRT(SC,NKN,'Regularization matrix') C MARKE=4711 * calculate contribution from regularization term XCX without tau CALL SMAVA(SC,XX,XCX,NKN,1) * * calculation of: * CALL RSOLRT(UM,HE,GR,XX,AM,DE,TAU,XDFR,NKN-1,NIKN) * * determine number of degrees of freedom * IF(IT.LT.6) THEN IF(NDF.GT.0) THEN NDFA=NDF ELSE IF(NDF.LT.0.AND.IT.LE.2) THEN NDFA=-NDF ELSE NDFA=XDFR+1.0 END IF END IF NDFA=MIN(NDFA,30) * * * CALL RSOLEQ(UM,HE,GR,XX,AM,DE,TAU,FMAX,IVETO,FLOAT(NDFA)) FXT=FUN+0.5*TAU*XCX * initial values set to step vector IF(IT.EQ.1) THEN * first iteration WRITE(LUNP,104) IT AVR=0.0 IVR=0 DO I=1,NKN IF(GR(I).LT.0.0) IVR=1 AVR=AVR+GR(I) XX(I)=GR(I) END DO IF(IVR.EQ.0) GOTO 10 AVR=AVR/FLOAT(NKN) DO I=1,NKN XX(I)=AVR END DO GOTO 10 END IF * not first iteration - check convergence 30 IF(IT.GT.3.AND.IVETO.EQ.0.AND. + (AMAX1(TEST,TESTL).LT.0.02.OR.TEST.EQ.0.0)) GOTO 70 IF(IT.EQ.12) THEN WRITE(6,106) GOTO 71 END IF * minimization along step vector, restricted to keep nonnegative * content in each knot interval IPR=2 WRITE(LUNP,102) IT,TAU,NDFA * start-init-MINST for minimum search * search in positive direction only CALL MINSTI(0.0,FMAX,IPR) C IPRNST=IPR C XLIMST=FMAX C NEVAST=0 * end-init-MINST * minimization along line GR(.) in MINST FAC=0.0 50 CALL RUNFUN(-IT) CALL RSOLRM(SC) CALL SMAVA(SC,XX,XCX,NKN,1) FXT=FUN+0.5*TAU*XCX CALL MINST(FXT,FAC,NC) DO I=1,NKN XX(I)=XX(I)+FAC*GR(I) END DO IF(NC.LT.0) GOTO 50 * MINST convergence or non-convergence TESTL=TEST TEST =DFUNST CALL MINSTE(TEST,BLA,BLA) * ... next iteration GOTO 10 * * convergence reached ---------------------------------------------- * 70 WRITE(LUNP,103) XDFR 71 WRITE(6,105) SD,SF,NIKN CALL SMTEXT('Unfolding - fit ending') * get covariance matrix CALL RSOLVS(UM,HE) CALL RUNFUN(-IT) NSYM=(NKN*NKN+NKN)/2 * store result: number of knots, limits, ... RSP(1)=NKN RSP(2)=XLOW RSP(3)=XHIG DO I=1,NKN * spline coefficients ... RSP(3+I)=XX(I) END DO DO I=1,NSYM * ... and covariance matrix RSP(3+NKN+I)=HE(I) END DO NWD=3+NKN+(NKN*NKN+NKN)/2 100 RETURN 101 FORMAT('1',12X,'I T E R A T I V E S O L U T I O N'/) 102 FORMAT('0Iteration',I3, 2 ' Tau =',G12.4,' for NDF = ',I3/) 103 FORMAT('0',12X,'C O N V E R G E N C E ', 1 12X,'Recommended NDF from noise analysis is:', 2 ' NDF > ',F7.2/) 104 FORMAT('0Iteration',I3,' least squares step to define initial', 1 ' values'/) 105 FORMAT('0Nr of events Data',F12.1,' Fit',F12.1/ + ' Nr of indep. amplitudes reduced to',I4) 106 FORMAT('0No convergence after 12 iterations - continue without') 107 FORMAT(2X,A) 108 FORMAT(2X,5E14.6) END SUBROUTINE RUNFUN(IARG) * * calculate function value, gradient and hessian * for least squares sum (abs(arg)<3), or * for maximum likelihood sum (abs(arg)>2) * for arg<0 calculate only the function value * * FUN = function value * GR(.) = gradient * HE(.) = Hessian, matrix of second derivatives * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== COMMON/UNCOMD/H(LENGDP) COMMON/RUNRES/XW(527),RSP(2082) * ------------------------------------------------------------------ * Internal common for matrices, vectors etc for up to 63 knots * and up to 30 data points PARAMETER (NKMAX=63,NDAMAX=30) COMMON/UINTER/UM(NKMAX*NKMAX),HE((NKMAX+NKMAX*NKMAX)/2), + GR(NKMAX),XX(NKMAX),AM(NKMAX),DE(NKMAX), + SCV(NKMAX), + SCM(NKMAX*NKMAX), + FUN,SF,SD,TAU * M=NKN-1 * UM = tranformation matrix * HE = hessian, replaced by covariance matrix * GR = gradient, replaced by step vector * XX = coefficient * AM = amplitudes M * DE = eigenvalues M * SC = scratch storage * ------------------------------------------------------------------ SAVE /RUNRES/,/UINTER/ * INPUT * NKN,NBN,H(NPROD+L),H(NPROD+NBN+L),H(NPROD+NBN+NBN+L) DOUBLE PRECISION CH * ... C WRITE(*,*) 'RUNFUN called with argument',IARG * clear scalars CH=0.0 SF=0.0 SD=0.0 IF(IARG.GT.0) THEN * clear gradient, Hessian IJ=0 DO I=1,NKN * gradient GR(I)=0.0 DO J=1,I IJ=IJ+1 * Hessian HE(IJ)=0.0 END DO END DO END IF * loop over all bins LKN=NPROD+NBN+NBN C WRITE(*,*) 'LKN,NPROD,NBN',LKN,NPROD,NBN DO L=1,NBN * data and ... D=H(NPROD+L) SD=SD+D * ... and fitted events (background plus Monte Carlo) F=H(NPROD+NBN+L)+VXY(XX,H(LKN+1),NKN) SF=SF+F IF(D.NE.0.0.OR.F.NE.0.0) THEN * data events and/or fitted events in bin IF(IABS(IARG).LE.1.OR.F.LT.1.0) THEN * least squares sum WT=1.0/AMAX1(1.0,D) * function sum CH=CH+WT*(F-D)**2 IF(IARG.GT.0) THEN * calculate ... IJ=0 DO I=1,NKN * ... gradient ... GR(I)=GR(I)+WT*(D-F)*H(LKN+I) DO J=1,I IJ=IJ+1 * ... and hessian HE(IJ)=HE(IJ)+WT*H(LKN+I)*H(LKN+J) END DO END DO END IF ELSE * Poisson (maximum likelihood) sum IF(D.EQ.0.0) THEN CH=CH+2.0*F-1.0 ELSE * constants (D) added, to have likelihood fcn like chisqr CH=CH-2.0*(D*ALOG(F/D)-(F-D)) END IF IF(IARG.GT.0) THEN * calculate ... IJ=0 DO I=1,NKN * ... gradient ... GR(I)=GR(I)+(D/F-1.0)*H(LKN+I) DO J=1,I IJ=IJ+1 * ... hessian HE(IJ)=HE(IJ)+D*H(LKN+I)/F*H(LKN+J)/F END DO END DO END IF END IF END IF LKN=LKN+NKN END DO FUN=CH C WRITE(*,*) 'FUN',FUN C WRITE(*,*) 'GR ',(GR(II),II=1,NKN) C WRITE(*,*) 'XX ',(XX(II),II=1,NKN) 100 RETURN END FUNCTION USFUN(XX) * * This function has to return the value of the function f(x) for the * true (undisturbed) argument value x, which has been used to create * the Monte Carlo events. * * The default values (essentially) returns a function value * according to the x-distribution of the Monte Carlo events, * assuming, that no generated events are lost. * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== COMMON/CUSFUN/SPHUFS(NHU) LOGICAL START,USERF,DEFIN EXTERNAL XFUNCT DATA START/.TRUE./ * ... IF(START) THEN START=.FALSE. DEFIN=.TRUE. TEST=AVFUN(XLOW,XHIG,XFUNCT) IF(TEST.EQ.0.0) THEN * XFUNCT seems to be dummy USERF=.FALSE. LUSFUN=4711 ELSE * XFUNCT is a users function USERF=.TRUE. END IF GOTO 100 END IF IF(USERF) THEN USFUN=XFUNCT(XX) ELSE IF(DEFIN) THEN DEFIN=.FALSE. XXA=(FLOAT(NHA+NHA-1)*XLOW+XHIG)/FLOAT(NHA+NHA) XXB=(XLOW+FLOAT(NHA+NHA-1)*XHIG)/FLOAT(NHA+NHA) END IF USFUN=GAMMA*USPF(XX,SPHUFS,NHA,XXA,XXB) END IF 100 RETURN END SUBROUTINE ANACOV(N,V,VV,U,X) * * analysis of covariance matrix * REAL V(1) REAL VV(*),U(N,N),X(*),Y(60),HANA(50) * ... * dim of VV as of V. X and Y are vectors of length N. CALL SMTEXT(' ') CALL SMTEXT('Analysis of covariance matrix') CALL SMTEXT('-----------------------------') CALL SMTEXT('The size of the correlations between different 1 data points is analysed. Deviations are generated 2000 times 2 by monte carlo, according to the true covariance matrix. ') CALL SMTEXT(' From the deviations a chi**2 value is calculated, 1 ignoring correlations, i.e. using the diagonal elements only. ') CALL SMTEXT('From the chi**2 values the probability is 2 calculated and shown below in a histogram. ') CALL SMTEXT('If the distribution of probability is almost 1 flat, the correlations can be ignored. ') CALL SMTEXT('This means, that curves from theoretical predictions 1 can be compared directly with the data points, taking into 2 account only the (diagonal) errors shown in the figure.') * transform vv to matrix u for Monte Carlo generation NSYM=(N*N+N)/2 DO I=1,NSYM VV(I)=V(I) END DO CALL SMHHL(VV,N,Y,U,-1) * transpose matrix DO I=1,N DO J=1,I IF(I.NE.J) THEN UIJ =U(I,J) U(I,J)=U(J,I) U(J,I)=UIJ END IF END DO END DO WRITE(6,104) WRITE(6,102) (Y(I),I=1,N) COND=Y(1)/Y(N) WRITE(6,103) COND DMIN=0.0 DO I=1,N IF(Y(I).GT.0.0) YMIN=Y(I) F=SQRT(AMAX1(DMIN,Y(I))) DO J=1,N U(I,J)=F*U(I,J) END DO END DO * * NRANS random data sets * NRANS=5000 * prepare histogram of probabilities DO I=1,50 HANA(I)=0.0 END DO * loop on NRANS cases SUM1=0.0 SUM2=0.0 DO L=1,NRANS * start with X = N independent Gaussians DO I=1,N X(I)=RNGA() END DO * transform X to Y CALL GMAB(U,X,Y,N,N,1) * sum pseudo chi square SUM=0.0 IJ=0 DO I=1,N IJ=IJ+I SUM=SUM+Y(I)**2/V(IJ) END DO * ... and probability from chi square PR=PROB(SUM,N) * ... into histogram IB=1.0+PR*50.0 IB=MAX(1,MIN(IB,50)) HANA(IB)=HANA(IB)+1 * calculate mean value of chisquare/ndf SUM=SUM/FLOAT(N) SUM1=SUM1+PR SUM2=SUM2+SUM END DO SUM1=SUM1/FLOAT(NRANS) SUM2=SUM2/FLOAT(NRANS) * print histogram CALL VPRINT(HANA,50,0.0,1.0,' Histogram HANA') 100 RETURN 102 FORMAT(1X,10G12.4) 103 FORMAT('0Condition number is',G12.4/) 104 FORMAT('0Eigenvalues of covariance matrix') END SUBROUTINE MBYFUN(NB,XC,RC,VC,FACTOR) * multiply all data by average function value of function USFUN(x) REAL XC(*),RC(*),VC(*) EXTERNAL USFUN DO I=1,NB FAC=FACTOR*AVFUN(XC(I),XC(I+1),USFUN) RC(I)=FAC*RC(I) IJ=(I*I-I)/2 DO J=1,NB IF(J.LE.I) IJ=IJ+1 VC(IJ)=FAC*VC(IJ) IF(I.EQ.J) VC(IJ)=FAC*VC(IJ) IF(J.GE.I) IJ=IJ+J END DO END DO END SUBROUTINE RURSLT(NB,BINLIM,IPRINT,WW) * replaces RSOLVR * User callable routine * RUN result in x bins for NB bins in bin limits * BINLIM(1)...BINLIM(NB+1) stored in WW * * NB = 0 automatic bin limits * NB < 0 equidistant bins * NB > 0 bins limits in BINLIM * * IPRINT = 0 result calculation only * = 1 compact printout of result array * = 2 printout of result * = 3 printer plot of result * = 3 analysis of covariance matrix and printout * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== COMMON/RUNRES/XW(527),RSP(2082) * ------------------------------------------------------------------ * Internal common for matrices, vectors etc for up to 63 knots * and up to 30 data points PARAMETER (NKMAX=63,NDAMAX=30) COMMON/UINTER/UM(NKMAX*NKMAX),HE((NKMAX+NKMAX*NKMAX)/2), + GR(NKMAX),XX(NKMAX),AM(NKMAX),DE(NKMAX), + SCV(NKMAX), + SCM(NKMAX*NKMAX), + FUN,SF,SD,TAU * M=NKN-1 * UM = tranformation matrix * HE = hessian, replaced by covariance matrix * GR = gradient, replaced by step vector * XX = coefficient * AM = amplitudes M * DE = eigenvalues M * SC = scratch storage * ------------------------------------------------------------------ SAVE /RUNRES/,/UINTER/ REAL BINLIM(*),WW(*) COMMON/SPCOM/I0,N0,S(100) * ... reset output array ... IF(IPRINT.GT.0)CALL SMTEXT(' ') IF(IPRINT.GT.0)CALL SMTEXT(' -> Subroutine RURLST starting ...') C WRITE(*,*) 'RURSLT called with NB,IPRINT= ',NB,IPRINT NN=NB XLOW=RSP(2) XHIG=RSP(3) IF(NN.EQ.0) THEN * Automatic bin definition NN=NDF IF(NN.LE.0) NN=NDFA IF(NN.EQ.0) NN=12 * bins from next (NN+1)-th) orthogonal function CALL RSOLVB(NN+1,UM,NN,WW(2)) * check failure ... IF(NN.GT.0) THEN IF(WW(2).EQ.WW(2+NN)) NN=-NN ELSE NN=NDF IF(NN.LE.0) NN=NDFA NN=-NN END IF ELSE IF(NN.GT.0) THEN * copy from argument DO I=1,NN+1 WW(1+I)=BINLIM(I) END DO END IF IF(NN.LT.0) THEN * Equidistant bins NN=-NN * ... calculate equidistant bin limits WW(2)=XLOW DO I=1,NN WW(2+I)=XLOW+FLOAT(I)*(XHIG-XLOW)/FLOAT(NN) END DO WW(2+NB)=XHIG END IF * ... reset output array ... DO I=NN+3,2+NN+NN+(NN+NN+NN)/2 WW(I)=0.0 END DO * ... and insert bin number WW(1)=NN * get result from internal spline array RSP NKN=RSP(1)+0.5 IF(NKN.LE.0) GOTO 100 * correct outer bin boundaries IF(WW( 2).LT.XLOW) WW( 2)=XLOW IF(WW(2+NN).GT.XHIG) WW(2+NN)=XHIG * define matrix for transformation to bin contents IJ=0 * set_up transformation matrix SCMA(.) DO J=1,NN * integration between bin boundaries EI=ESPI(WW(J+1),WW(J+2),RSP(4),NKN,XLOW,XHIG) DO I=1,NKN IJ=IJ+1 * matrix element (S(I) is spine weight determined in ESPI) SCM(IJ)=S(I)/(WW(J+2)-WW(J+1)) END DO END DO * transform to bin content RC ... CALL GMAB(SCM,RSP(4),WW(NN+3),NN,NKN,1) * .. and transform covariance matrix VC CALL SMAVA(RSP(4+NKN),SCM,WW(NN+NN+3),NKN,NN) * finally scale data by user function USFUN(X) CALL MBYFUN(NN,WW(2),WW(NN+3),WW(NN+NN+3),FACLUM) IF(IPRINT.LE.0) GOTO 100 * ========================================== * print-out of result NN=WW(1)+0.5 NWD=NN+NN+(NN*NN+NN)/2+2 IF(IPRINT.GE.2) THEN CALL BGTEXT('PAGE') CALL BGTEXT(' - R E S U L T - ') * * print result * II=0 SUM1=0.0 SUM2=0.0 WRITE(6,106) DO I=1,NN II=II+I DY=SQRT(WW(NN+NN+2+II)) YY=WW(NN+2+I) XE=WW(1+I) XR=WW(2+I) YINT=YY*(XR-XE) EVEF=0.0 IF(DY.GT.0.0) EVEF=(YY/DY)**2 SUM1=SUM1+YINT SUM2=SUM2+EVEF WRITE(6,104) I,XE,XR,YY,DY,YINT,SUM1,EVEF,SUM2 END DO CALL SMPRV(WW(NN+3),WW(NN+NN+3),NN) END IF * IF(IPRINT.GE.1) THEN WRITE(*,107) 'Result using compact format (2X,5E14.6)' WRITE(*,107) TITLE WRITE(*,108) (WW(I),I=1,NWD) WRITE(*,108) END IF * IF(IPRINT.GE.3) THEN * * Printer plot of result * CALL FPXY(' ',0.5*(XLOW+XHIG),0.0) DO I=1,101 XF=0.01*(FLOAT(101-I)*XLOW+FLOAT(I-1)*XHIG) TX=ESPF(XF,RSP(4),NKN,XLOW,XHIG)*USFUN(XF) CALL FPNXY('*',1,XF,FACLUM*TX) END DO II=0 DO I=1,NN II=II+I DY=SQRT(WW(NN+NN+2+II)) YY=WW(NN+2+I) XE=WW(1+I) XR=WW(2+I) CALL FPXDYD('O',0.5*(XE+XR),0.5*(XR-XE),YY,DY) END DO CALL FPC('O = DATA POINT = INTEGRAL OVER BIN') CALL FPL END IF * IF(IPRINT.GE.3) THEN * printout of analysis of covariance matrix CALL SMTEXT('PAGE') CALL ANACOV(NN,WW(NN+NN+3),SCM,HE,GR) END IF IF(IPRINT.GT.0) CALL SMTEXT(' -- Subroutine RURLST ending.') IF(IPRINT.GT.0) CALL SMTEXT(' ') 100 RETURN 104 FORMAT(1X,I5,4G12.4,3X,2G12.4,2X,2G12.4) 106 FORMAT('0 I --- R A N G E --- AV.F(X) ERROR A 1V.F(X)*BIN SUM EFF.EV. SUM'/) 107 FORMAT(2X,A) 108 FORMAT(2X,5E14.6) END FUNCTION AVFUN(X1,X2,FUN) * Calculate average function value of function fun(x) in (x1,x2) * Simpson integration rule is used. * Relative precision is at least 0.001. EXTERNAL FUN DOUBLE PRECISION SUM * ... AVFUN=0.0 N=16 DO L=1,8 SUM=0.0 * form simpson sum DO I=0,N X=(FLOAT(N-I)*X1+FLOAT(I)*X2)/FLOAT(N) F=FUN(X) IF(I.EQ.0.OR.I.EQ.N) THEN SUM=SUM+F ELSE IF(MOD(I,2).EQ.0) THEN SUM=SUM+2.0*F ELSE SUM=SUM+4.0*F ENDIF END DO AVFUL=AVFUN AVFUN=SUM/DFLOAT(3*N) * estimate precision by comparison with previous value IF(L.GE.2) THEN IF(ABS(AVFUN-AVFUL).LE.0.0005*(AVFUN+AVFUL)) RETURN END IF N=N+N END DO END * RUSOLM file SUBROUTINE RSOLVE(NKN,XLOW,XHIG,ILA,IRA,WLIN) * * Define all necessary parameters of regularization mathematics * * ------------------------------------------------------------------ * M=NKN-1 * UM = tranformation matrix * HE = hessian, replaced by covariance matrix * GR = gradient, replaced by step vector * XX = coefficient * AM = amplitudes M * DE = eigenvalues M * SCV = scratch storage for vector * SCM = scratch storage for matrix * DGH = eigenvalues of Hessian H * DGM = eigenvalues of Hessian H, modified * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ * ... * copy arguments N = NKN ! number of knots XA = XLOW ! limit left XB = XHIG ! limit right * (re)set flags NC = 0 ! MC = 0 ! number of constraint pairs ICFL = 0 ! flag (=1 means non-negative solution) ILF = ILA ! flag for left limit (=1 no constraint) IRF = IRA ! flag for right limit (=1 no constraint) WLN = WLIN ! relative weight for linear weight * define parameter for conditions: * IL = 0 function zero outside left, otherwise nonzero * IR = 0 function zero outside right, otherwise nonzero * WLN= relative weight for linear regularization term * determine precision of computer DO IP=1,20 EPSMAP=1.0+10.0**(-IP) IF(EPSMAP.EQ.1.0) GOTO 10 END DO IP=20 10 DO JP=1,100 EPSMA=0.5*FLOAT(JP)*10.0**(-IP) EPSMAP=1.0+EPSMA IF(EPSMAP.NE.1.0) GOTO 20 END DO * default precision 20 PREC=100.0 RETURN END SUBROUTINE RSOLFC(XV,FV) * * Add flags for constraint * * function(XV) = FV forced * * special case XV = FV = 0 means: function non-negative forced * * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ * ... IF(XV.EQ.0.0.AND.VL.EQ.0.0) THEN * set flag to one ICFL=1 ELSE IF(MC.LT.10) THEN * store new function value ? MC=MC+1 NC=MC DX=(XB-XA)/FLOAT(N-3) IA=(XV-XA)/DX+1.0 * force 1 =< IA =< N-3 IA=MIN0(MAX0(1,IA),N-3) IC(MC)=IA CN(MC)=FV*24.0 END IF RETURN END SUBROUTINE RSOLRM(SC) * * CALL RSOLVM(SC) * -- * * Define regularization matrix * REAL SC(*) * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ REAL CC(6),CB(6),CD(6) * Constants for regularization matrix * second derivative DATA CC/70.0,40.0,10.0,0.0,30.0,15.0/ * first derivative DATA CB/74.0,40.0,6.0,12.0,22.0,-7.0/ * ... IA=ILF IB=IRF * factor EPS is linear weight EPS=WLN * * square of second derivative -------------------------------------- * IJ=0 DO I=1,N IR=N+1-I DO J=1,I IJ=IJ+1 IF(I.EQ.J) THEN * diagonal SC(IJ)=80.0 IF(IA.NE.0.AND.I .LE.3) SC(IJ)=SC(IJ)-CC(I ) IF(IB.NE.0.AND.IR.LE.3) SC(IJ)=SC(IJ)-CC(IR) ELSE IF(I.EQ.J+1) THEN * off-diagonal +1 SC(IJ)=-45.0 IF(IA.NE.0.AND.I .LE.3) SC(IJ)=SC(IJ)+CC(3+I ) IF(IB.NE.0.AND.IR.LE.2) SC(IJ)=SC(IJ)+CC(4+IR) ELSE IF(I.EQ.J+3) THEN * off-diagonal +3 SC(IJ)=5.0 ELSE * all remaining elements are zero SC(IJ)=0.0 END IF END DO END DO * * square of first derivative --------------------------------------- * EPS1=80.0*EPS EPS2=15.0*EPS EPS3=24.0*EPS * scale first derivative by EPS DO I=1,6 CD(I)=EPS*CB(I) END DO IJ=0 DO I=1,N IR=N+1-I DO J=1,I IJ=IJ+1 IF(I.EQ.J) THEN * diagonal SC(IJ)=SC(IJ)+EPS1 IF(IA.NE.0.AND.I .LE.3) SC(IJ)=SC(IJ)-CD(I ) IF(IB.NE.0.AND.IR.LE.3) SC(IJ)=SC(IJ)-CD(IR) ELSE IF(I.EQ.J+1) THEN * off-diagonal +1 SC(IJ)=SC(IJ)-EPS2 IF(IA.NE.0.AND.I .LE.3) SC(IJ)=SC(IJ)+CD(3+I ) IF(IB.NE.0.AND.IR.LE.2) SC(IJ)=SC(IJ)+CD(4+IR) ELSE IF(I.EQ.J+2) THEN * off-diagonal +2 SC(IJ)=SC(IJ)-EPS3 IF(IA.NE.0.AND.I.LE.3) SC(IJ)=SC(IJ)+CD(4) IF(IB.NE.0.AND.I.EQ.N) SC(IJ)=SC(IJ)+CD(4) ELSE IF(I.EQ.J+3) THEN * off-diagonal +3 SC(IJ)=SC(IJ)-EPS END IF END DO END DO RETURN END SUBROUTINE RSOLRT(U,H,G,X,A,D,TAU,XDFR,MLIM,NIKN) * * - - * CALL RSOLRT(U,H,G,X,A,D,TAU,XDFR,MLIM) * H(.) = Hessian * G(.) = gradient * * * ---- - - * CALL RSOLV1(M,MLIM,U,H,G,X,A,SC) * - - - | * +> N WORDS * * determine K'th orthogonal function, K = 1...M * (alternative to RSOLV0, if M is known before) * The K'th function is first defined as the legendre polynomial of * order (K-1). It is orthonormalized with metric H by subtraction of * the components of all previous orthogonal functions. The spline * coefficients are stored in the K'th row of the matrix U. The * amplitude for the K'th function is calculated by * T * A(K) = G U (ROW K) * * transformation of matrix H to diagonal * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ REAL X(*),G(*),H(*),U(*),A(*),D(*) * ... * ------------------------------------------------------------------ * Transform Hessian to diagonal matrix D with transf. matrix U * * copy Hessian H to scratch storage DO I=1,(N*N+N)/2 SCMA(I)=H(I) END DO C CALL SMPR(SCMA,N,'before SMHHL') * transformation of H to diagonal form CALL SMHHL(SCMA,N,DGH,U,-1) C CALL SMPR(SCMA,N,'after SMHHL') * ... to DGH(1) ... DGH(N) = diagonal values, decreasing in size * U(.) = transformation matrix * * eventually reduce size of matrix, ignoring all negative (!) * eigenvalues and the smallest positive eigenvalue * define M to ignore all negative eigenvalues and ... * ... and the smallest positive eigenvalue, not more than MLIM * analyse eigenvalue spectrum NPOS=0 DO I=1,N DGM(I)=DGH(I) IF(DGH(I).GT.0.0) THEN * index of positive eigenvalue NPOS=I END IF END DO C DGLOW=AMAX1(DGH(1)*(PREC*EPSMA)**2,DGH(NPOS-1)) C DO I=1,N * modify very small eigenvalues C DGM(I)=DGH(I) C DGM(I)=AMAX1(0.0,DGH(I))+DGLOW C END DO * divide eigenvectors by square root of eigenvalue C M=N * test: M=NPOS-1 C write(*,*) 'npos ',npos M=NPOS-1 * number of independent linear combinations NIKN=M IJ=0 DO J=1,M FACT=1.0/SQRT(DGM(J)) * test: FACT=0 for all negative and the smallest positive EV IF(J.GE.NPOS) FACT=0.0 * end of test modification ===================xumod========= DO I=1,N IJ=IJ+1 U(IJ)=U(IJ)*FACT END DO END DO * ------------------------------------------------------------------ * complete transformation to system, which is diagonal in reg. matr. * matrix U and amplitudes are redefined * * define regularization matrix CALL RSOLRM(SCMA) * transform SCMA using matrix U -> SCMB CALL SMAVA(SCMA,U,SCMB,N,M) * * now SCMB contains the matrix C1 of the paper * * spectral decomposition of C1 (increasing order of eigenvalues) CALL SMHHL(SCMB,M,D,SCMA,+1) * D(.) = diagonal matrix (= matrix S of the paper) * SMCA(.) = rotation matrix * * insert exact zeros as eigenvalues, if appropriate IF(ILF.EQ.1.AND.IRF.EQ.1) THEN D(1)=0.0 IF(WLN.EQ.0.0) D(2)=0.0 END IF DO I=1,M C write(*,*) i,'D(i)=',D(I) END DO * form complete transformation matrix SCMB and store in U CALL GMAB(SCMA,U,SCMB,M,M,N) IJ=0 DO J=1,M DO I=1,N IJ=IJ+1 U(IJ)=SCMB(IJ) END DO END DO * now in array U the matrix * U2T * D**(-1/2) * U1T * of the paper is stored * transform gradient CALL GMAB(U,G,SCV,M,N,1) * transform previous solution vector CALL SMVA(H,X,SCMA,N,1) CALL GMAB(U,SCMA,A,M,N,1) * add both to obtain normalized amplitudes DO I=1,M A(I)=A(I)+SCV(I) END DO * Now the vector A contains the vector * a_bar_prime = unregularized solution is stored * * ------------------------------------------------------------------ * Determine estimate of TAU by following method: * * For a given value of TAU the difference * A(i)**2 - (damped A(i))**2 * can be considered as the random noise, which on average should be * 1. Thus one could determine the value of TAU, for which the random * noise contributions are 1 on average. Sometimes there are larger * values of amplitudes at a high index, which disturb the algorithm. * Therefore the contributions are weighted by * 3/(3+TAU*D(i)) * which is nearly one in the region of significant amplitudes and * also in the transition region, but small in the noise region. * * NOISE = Ai**2 - (DAMPED Ai)**2 * = Ai**2 * (1 - 1/(1+TAU*Di)**2) * = Ai**2 * (TAU*Di)*(2+TAU*Di)/(1+TAU*Di)**2 * expected is NDF = SUM of noise * stop SUM if 0.001 < tau D < 10.0 DO J=3,M IF(D(J).GT.0.0) THEN * try this value for TAU TAU=1.0/D(J) FUN=0.0 IA=0 IB=0 DO I=1,M TD=TAU*D(I) IF(TD.GT.0.001.AND.TD.LT.10.0) THEN IF(IA.EQ.0) IA=I WT=3.0/(3.0+TD) FB=A(I)**2*(2.0+TD)*TD/(1.0+TD)**2-1.0 FUN=FUN+WT*FB IB=I END IF C WRITE(*,222) J,TAU,I,WT,FB,FUN,TD 222 FORMAT(' J,TAU,I,WT,FB,FUN ',I3,G12.3,I3,4G12.4) END DO IF(FUN.LT.0.0) GOTO 10 END IF END DO * * FUN is negative the first time. Use the current value of TAU as * estimate and improve TAU by iterative method. * * 10 IT=0 C WRITE(*,*) FUN,TAU,IA,IB,' IA IB' 20 IT=IT+1 FUN =0.0 DFUN=0.0 DO I=IA,IB TD=TAU*D(I) WT=3.0/(3.0+TD) FB=A(I)**2*(2.0+TD)*TD/(1.0+TD)**2-1.0 * function and ... FUN=FUN +WT*FB DA=-D(I)*WT/(3.0+TD) DB=A(I)**2*2.0*D(I)/(1.0+TD)**3 * ... and derivative DFUN=DFUN+WT*DB+DA*FB C WRITE(*,111) 'RSOLRT',I,DFUN,FA,FB,DA,TAU,D(I),A(I) 111 FORMAT(1X,A,I3,8G12.4) END DO C write(*,*) 'fun dfun ',fun,dfun IF(DFUN.EQ.0.0) THEN DTAU=0.0 ELSE DTAU=-FUN/DFUN END IF TAU =TAU+DTAU C WRITE(*,*) 'ITERATION ',IT,TAU,DTAU IF(ABS(DTAU).GT.TAU*1.0E-5.AND.IT.LE.10) GOTO 20 * convergence reached * Determine nr of degrees of freedom for the current value of TAU XDFR=0.0 DO I=1,M XDFR=XDFR+1.0/(1.0+TAU*D(I)) END DO C WRITE(*,*) '----------> XDFR = ',XDFR RETURN END SUBROUTINE RSOLEQ(U,H,G,X,A,D,TAU,FMAX,IVETO,XDF) * * determine tau and solve equations, observing constraints * REAL U(*),H(*),G(*),X(*),A(*),D(*) * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ REAL B1(100),B2(100),A2(100),AS(100) * ... * regularization parameter tau is determined from the identity * * M 1 * XDF = Sum ---------- * I=1 1+TAU*D(I) TAU=0.0 IF(XDF.LE.0.0.OR.XDF.GE.FLOAT(M)) GOTO 20 * loop to find starting value for tau DO J=2,M-1 IF(D(J).NE.0.0) THEN TAU=1.0/D(J) FS=-XDF DO I=1,M FS=FS+1.0/(1.0+TAU*D(I)) END DO IF(FS.GT.0.0) THEN * Newton refinement of tau IT=0 10 IT=IT+1 FS=-XDF FP=0.0 DO I=1,M FS=FS+1.0/(1.0+TAU*D(I)) FP=FP+D(I)/(1.0+TAU*D(I))**2 END DO DTAU=+FS/FP DELTA=ABS(DTAU/TAU) TAU=TAU+DTAU C WRITE(*,*) 'RSOLEQ TAU DTAU ',TAU,DTAU IF(DELTA.GT.1.0E-5.AND.IT.LT.10) GOTO 10 * tau accurately enough GOTO 20 END IF END IF END DO TAU=0.0 * * Calculate expected change of function value test and step vector g * using regularization parameter tau by the formulas * * * 1 * G = U A' - X with A' = ---------- A(I) * I 1+TAU*D(I) * * T * test= 0.5 * g * (gradient - tau C x) C = regularization matrix * 20 IR =0 IVETO=0 * define inverse matrix 30 IJ=0 DO I=1,M A2(I)=A(I) DO J=1,I IJ=IJ+1 SCMA(IJ)=0.0 END DO * diagonal element not zero SCMA(IJ)=1.0/(1.0+TAU*D(I)) END DO * same value of M ME=M * add known constraints DO L=1,NC IA=IC(L) DO J=1,M IJ=(J-1)*N+IA-1 IF(L.GT.MC) THEN * positivity constraint B1(J)=U(IJ+1) ELSE * function value constraint B1(J)=U(IJ+1)+11.0*(U(IJ+2)+U(IJ+3))+U(IJ+4) END IF END DO * update inverse matrix B1(ME+1)=0.0 CALL SMVA(SCMA,B1,B2,ME,1) EL=-1.0/VXY(B1,B2,ME) IJ=0 DO I=1,ME DO J=1,I IJ=IJ+1 SCMA(IJ)=SCMA(IJ)+EL*B2(I)*B2(J) END DO END DO DO I=1,ME IJ=IJ+1 SCMA(IJ)=-EL*B2(I) END DO SCMA(IJ+1)=EL ME=ME+1 IF(L.GT.MC) THEN * positivity constraint A2(ME)=0.0 ELSE * function value constraint A2(ME)=CN(L) END IF END DO * solve system by multiplication ... CALL SMVA(SCMA,A2,AS,ME,1) * ... and back by rotation to G CALL GMAB(AS,U,G,1,M,N) IF(IR.NE.0) GOTO 90 * maximum step FMAX=5.0 IF(ICFL.EQ.0) GOTO 90 * check negative content or step to negative content INEG=0 FNEG=0.0 DO 40 I=1,N DO L=MC+1,NC IF(IC(L).EQ.I) GOTO 40 END DO FN=X(I) FD=G(I)-FN IF(FD.NE.0.0) THEN FCA=-FN/FD IF(FCA.LT.0.1.AND.FD.LT.0.0) THEN IF(INEG.EQ.0.OR.FD.LT.FNEG) THEN * set flag INEG to index with largest negative value INEG=I FNEG=FD C WRITE(*,*) ' NEGATIVE>> :',INEG,FNEG END IF ELSE IF(FCA.GT.0.0) THEN * determine FMAX as maximum allowed step vector FMAX=AMIN1(FMAX,FCA) END IF END IF 40 CONTINUE IF(INEG.NE.0) THEN * add zero constraint and ... IF(NC.GE.110) GOTO 90 NC=NC+1 IVETO=1 * ... store index which could go negative ... IC(NC)=INEG * ... and repeat setting up equations GOTO 30 END IF * check, wether one constraint can be removed (step positive) IM =0 ASM=0.0 DO L=MC+1,NC * check Lagrangian IF(IR.NE.0.OR.AS(L+M).GT.ASM) THEN IR=L ASM=AS(L+M) END IF END DO IF(IR.EQ.0.OR.ASM.LT.0.0) GOTO 90 * remove constraint ... DO L=IR+1,NC IC(L-1)=IC(L) END DO NC=NC-1 IVETO=1 * ... and repeat setting up equations GOTO 30 * finally calculation of step as difference 90 DO I=1,N G(I)=G(I)-X(I) END DO RETURN END SUBROUTINE RSOLVS(U,H) * * calculate covariance matrix H by the formula * * T 1 * H = U V' U with V' = --------------- * II (1+TAU*D(I))**2 * REAL U(*),H(*) * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ * ... * SCMA still keeps inverse * square of symmetric matrix (M-by-M) CALL SMQU(SCMA,SCMB,M) * transpose of matrix CALL GMTR(U,SCMA,M,N) * product CALL SMAVA(SCMB,SCMA,H,M,N) RETURN END SUBROUTINE RSOLVB(K,U,NB,XC) * * Determine limits XC(1)...XC(NB+1) from extremas of K'th function * - - * CALL RSOLVB(K,U,NB,XC) * -- -- * REAL U(*),XC(*) * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ * ... XC(1)=XA * left edge FL=ESPF(XA,U(1+(K-1)*N),N,XA,XB) L=1 XL=XA * next argument with first derivative = 0.0 10 XNEW=ESPX(XL,0.0,U(1+(K-1)*N),N,XA,XB,1) * ... and function value at this argument FNEW=ESPF(XNEW, U(1+(K-1)*N),N,XA,XB) IF(XNEW.LE.XB) THEN IF(FL*FNEW.LT.0.0) THEN * different sign of derivative L=L+1 XC(L)=XNEW FL=FNEW END IF XL=XNEW GOTO 10 END IF * right edge FNEW=ESPF(XB ,U(1+(K-1)*N),N,XA,XB) IF(FL*FNEW.GT.0.0) L=L-1 XC(L+1)=XB * final number of bins (NB) +1 is number of bin limits NB=L RETURN END SUBROUTINE RSOLVP(MODE) * * print parameter from internal common * * internal common with parameters and scratch arrays --------------- PARAMETER (NKMAX=63,NDAMAX=30) COMMON/CSOLVE/N,M,XA,XB,NC,MC,ICFL,ILF,IRF,WLN,EPSMA,EPSMAP,PREC, + IC(110),CN(10),SCV(NKMAX), + SCMA(NKMAX*NKMAX),SCMB(NKMAX*NKMAX), + DGH(NKMAX),DGM(NKMAX) * ------------------------------------------------------------------ * ... IF(MODE.LE.0) THEN * print initial parameters WRITE(*,101) N,XA,XB,NC,MC IF(ILF.NE.1) THEN WRITE(*,102) ' Result forced to zero at left limit' END IF IF(IRF.NE.1) THEN WRITE(*,102) ' Result forced to zero at right limit' END IF WRITE(*,102) ' Relative linear weight',WLN WRITE(*,102) ' Relative machine precision',EPSMA WRITE(*,102) ' Precision parameter',PREC END IF RETURN 101 FORMAT(' Number of knots = number of spline functions ',I6/ + ' left X limit ',G12.6/ + ' right X limit ',G12.6/ + ' NC, MC=',2I6) 102 FORMAT(1X,A40,1X,G9.2) C 104 FORMAT(10X,I3,2X,2(G13.4,G13.3,3X)) END * RUNPRT file SUBROUTINE RUNPRT * * interface between RUN and printer graphics * * UNCOMA: ntuple data accesible to user ---------------------------- COMMON/UNCOMA/IDMB,NY,WT,X,Y(100) REAL REC(102),XY(0:100) EQUIVALENCE (REC(1),WT),(XY(0),X) * IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background * NY = number of measured variables * WT = weight * X = true variable (defined only in Monte Carlo) = XY(0) * Y(j) = measured variables, j = 1 ... NY = XY(j) * REC(.) = array with: WT, X, Y(.) (by equivalence) * ================================================================== * UNCOMB: internal variables and arrays ---------------------------- PARAMETER (LENGDP=50000) PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3)) COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU), + HUP(NHU),ASPMC(NKNMAX), + NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG, + XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD, + NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN, + WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3), + NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM, + HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32) REAL STAB(12,33) EQUIVALENCE (JTAB(1,1),STAB(1,1)) * * ITAB = index of last ntuple * NTAB = number of different ntuple types * J/STAB( 1,.) = 1 data =2 moca =3 bg * J/STAB( 2,.) = hollerith text * J/STAB( 3,.) = * J/STAB( 4,.) = Lumi Lumi evtnr (first argument * J/STAB( 5,.) = * J/STAB( 6,.) = * J/STAB( 7,.) = * J/STAB( 8,.) = * J/STAB( 9,.) = 1.0 wfac2 (second argument) * J/STAB(10,.) = * J/STAB(11,.) = ifg for LOOK ntuples * J/STAB(12,.) = inr for LOOK ntuples * * INWT = index weight * INIX = index X * * * NV = number of fit variables (from text input) * NDV(j) = index of fit variable - " - * NBV(j) = number of bins - " - * * XLOW, XHIG = lower and upper limit of x (from text input) * mandatory, not modified by program * * XMIN(j),XMAX(j) = min and max value of Y(j), determined from * ntuples with positive weight (except background) * JXY(j) = 0 if XMIN/XMAX undefined * JXY(j) = 1 if XMIN/XMAX defined * * XYL(j),XYH(j) = limits of histogram, rounded for 120 bins, * for final histograms * * XYFL(K),XYFH(K) = limits as used by unfolding fit * * * Where are bin limits used: * * for unfolding fit, usually with small number of bins: XYFL/XYFH * * for final histograms, with 120 bins: XYL/XYH * * * * * HIST(.,.,.) = histogram of meaured variables, filled in PASS 3 * ================================================================== * UNCOMC: character strings----------------------------------------- CHARACTER*50 TITLE,VTEXT CHARACTER*8 SOPTE COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10) * TITLE = text as comment for problem * VTEXT(j) = text as comment for variable j * SOPTE(i) = selected options * ------------------------------------------------------------------ SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/ * ================================================================== REAL HUT(120),BX(2),BY(2) * ... CALL SMTEXT(' ') CALL SMTEXT('-> Subroutine RUNPRT starting ...') * fill HUT with dependence of user function USFUN DO I=1,120 X=XLOW+(FLOAT(I)-0.5)*(XHIG-XLOW)/120.0 HUT(I)=USFUN(X) END DO CALL SMTEXT(' ') CALL SMTEXT('X-dependence of function USFUN') CALL VPRINT(HUT,120,XLOW,XHIG,'X-dependenxe of USFUN') CALL FPL * * histograms of variables after pass3 * DO J=1,10 IF(JXY(J).NE.0) THEN * determine bins from data CALL HPRB(HIST(1,1,J),120) * plot background by B BY(1)=0.0 DO I=1,120 BY(2)=HIST(I,3,J) IF(BY(2).NE.0.0) THEN BX(1)=XYL(J)+(FLOAT(I)-0.5)*(XYH(J)-XYL(J))/120.0 BX(2)=BX(1) CALL FPNXY('B',2,BX,BY) END IF END DO * plot expected distribution (= MOCA + background) CALL HPRC(HIST(1,2,J),120,XYL(J),XYH(J)) * plot data CALL HPRH(HIST(1,1,J),HIST(1,1,J),120,XYL(J),XYH(J)) CALL FPC(VTEXT(J)) CALL FPC('DATA WITH MOCA-CURVE B=BACKGROUND') CALL FPL * calculate difference DO I=1,120 HIST(I,3,J)=HIST(I,1,J)-HIST(I,2,J) END DO CALL HPRH(HIST(1,3,J),HIST(1,1,J),120,XYL(J),XYH(J)) CALL FPC(VTEXT(J)) CALL FPNXY('.',1,XYL(J),0.0) CALL FPNXY('.',1,XYH(J),0.0) CALL FPC('DIFFERENCE DATA MINUS FIT') CALL FPL END IF END DO CALL SMTEXT('-> Subroutine RUNPRT ending.') CALL SMTEXT(' ') 100 RETURN END