PROGRAM MAKE_PARAM IMPLICIT NONE INTEGER HMEMOR,hsize PARAMETER (HSIZE=2000000) COMMON/PAWC/HMEMOR(HSIZE) #include "make_param.inc" CALL HLIMIT(HSIZE) CALL READSETUP(1) ! READ SETUP AND PUT STEERING VARIABLES IN COMMON CALL NTOPEN C C DEFINE HISTOS + FILL THEM + BUILD CORRELATIONS C C LEVEL OF HISTOGRAM RANGE DETECTION "REFINEMENT" IS DETERMINED BY THE C FLAG: 1(POOR) -> 5(GOOD) C CALL DEFHIST(5) CALL NTCLOSE(1) PRINT*,'------------ Now all the final fits --------------' CALL FITPARAM(1) PRINT*,'--------------------- Done -----------------------' IF (NBIN1.GT.3.AND.NBIN2.GT.3.AND.SMOOTHFLAG.GT.0) THEN CALL SMOOTHING ELSEIF (SMOOTHFLAG.GT.0) THEN PRINT*,'No smoothing done !!!!' PRINT*,'The reason is that you have only 1-d histos or .le.' & //' 3x3 bins.' PRINT*,'You should do interpolation instead (on your own).' ELSE PRINT*,'No smoothing done !!!!' ENDIF IF (OUTPUT_FORM.EQ.1) CALL WRITE_INC ! OUTPUT PARAMETERS + CORRELATION MATRICES IF (OUTPUT_FORM.EQ.2) CALL OUTPARAM ! OUTPUT PARAMETERS + CORRELATION MATRICES 999 CONTINUE CALL NTCLOSE(2) END C====================================================================== C Readin setup and important variables C====================================================================== SUBROUTINE READSETUP(IFLAG) IMPLICIT NONE INTEGER IFLAG #include "make_param.inc" #include "abcfit_bmatrix.inc" LOGICAL LEXIST INTEGER INPUT_CARDS_UNIT,I,NUMCLASS,J,INIT1,PPOINT(MAXNUMOFPAR) PARAMETER (INIT1=MAXNUMOFPAR+1) REAL RTEMP1,RTEMP2 CHARACTER*70 INPUT_CARDS_FILE,DUMMY_CHAR DATA PARAMFLAG /INIT1*-1/ C CALL VZERO(URANGE,2*MAXNUMOFHIS) C C GET UNITS TO USE C INPUT_CARDS_UNIT=90 LEXIST=.TRUE. DO WHILE (LEXIST) INQUIRE(UNIT=INPUT_CARDS_UNIT,OPENED=LEXIST) INPUT_CARDS_UNIT=INPUT_CARDS_UNIT+1 ENDDO INPUT_CARDS_UNIT=INPUT_CARDS_UNIT-1 UNIT_IN=INPUT_CARDS_UNIT+1 LEXIST=.TRUE. DO WHILE (LEXIST) INQUIRE(UNIT=UNIT_IN,OPENED=LEXIST) UNIT_IN=UNIT_IN+1 ENDDO UNIT_IN=UNIT_IN-1 UNIT_CHECK=UNIT_IN+1 LEXIST=.TRUE. DO WHILE (LEXIST) INQUIRE(UNIT=UNIT_CHECK,OPENED=LEXIST) UNIT_CHECK=UNIT_CHECK+1 ENDDO UNIT_CHECK=UNIT_CHECK-1 UNIT_OUT=UNIT_CHECK+1 LEXIST=.TRUE. DO WHILE (LEXIST) INQUIRE(UNIT=UNIT_OUT,OPENED=LEXIST) UNIT_OUT=UNIT_OUT+1 ENDDO UNIT_OUT=UNIT_OUT-1 C C READ SETUP FROM FILE C INPUT_CARDS_FILE='make_param.in' INQUIRE (FILE=INPUT_CARDS_FILE,EXIST=LEXIST) IF (.NOT.LEXIST) THEN PRINT*,'Cannot find input cards-file: ',INPUT_CARDS_FILE STOP ENDIF PRINT*,'Opening cards-file: ', INPUT_CARDS_FILE OPEN (INPUT_CARDS_UNIT,FILE=INPUT_CARDS_FILE,ACCESS= & 'SEQUENTIAL',FORM='FORMATTED',STATUS='OLD') READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(A70)') INPUT_FILE READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I6)') OUTPUT_FORM READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(A70)') OUTPUT_FILE_CHECK READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I6)') BINLEVEL READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I6)') EVENTFLAG READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,*) EVENTCUT READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I2)') MATCHFLAG READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I6)') USENUMOFENTRIES READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I2)') NUMCLASS READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR DO I=1,NUMCLASS READ(INPUT_CARDS_UNIT,'(I6)') PPOINT(I) ENDDO READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR DO I=1,NUMCLASS READ(INPUT_CARDS_UNIT,'(I6)') PARAMFLAG(PPOINT(I)) ENDDO READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR DO I=1,NUMCLASS READ(INPUT_CARDS_UNIT,*) B1MIN(PPOINT(I)),B1MAX(PPOINT(I)) & ,B2MIN(PPOINT(I)),B2MAX(PPOINT(I)),USEPAR(1,PPOINT(I)) & ,USEPAR(2,PPOINT(I)) ENDDO READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(2I)') NBIN1,NBIN2 READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I)') FITFLAG READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I)') NBINFIT READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I)') FITFUNC READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,*) SCALEFRAC SCALEFRAC=1.0-SCALEFRAC READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,*) RANGESCALECUT READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(I2)') SMOOTHFLAG READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR DO I=1,NUMCLASS DO J=1,NPAR(PARAMFLAG(PPOINT(I))) READ(INPUT_CARDS_UNIT,'(2I)') & SMOOTHFUNC(J,PPOINT(I)),SMOOTHFUNCERR(J,PPOINT(I)) ENDDO ENDDO READ(INPUT_CARDS_UNIT,'(A50)') DUMMY_CHAR DO WHILE(.TRUE.) READ(INPUT_CARDS_UNIT,*,ERR=900) I,RTEMP1,RTEMP2 URANGE(1,I)=RTEMP1 URANGE(2,I)=RTEMP2 ENDDO 900 CONTINUE CLOSE (INPUT_CARDS_UNIT) 999 RETURN END C====================================================================== C DEFINE HISTOGRAM RANGES + FILL THEM + BUILD CORRELATIONS C====================================================================== SUBROUTINE DEFHIST(IFLAGIN) IMPLICIT NONE INTEGER IFLAGIN #include "make_param.inc" #include "abcfit_bmatrix.inc" LOGICAL HEXIST INTEGER IFLAG,H,I,J,K,L,M,N,IPAR(MAXNUMOFPAR),IHIST,ISTAT,IBIN & ,NPARAMS,IPARAMS,HNLEN,NBINOK REAL RMEAN,LOWER,UPPER,HI,HSTATI,ELOW,EUP,CTHLOW,CTHUP,HSUM,BINMAX DOUBLE PRECISION BOUND(5,4,MAXNUMOFPAR,0:400),SUM(2,25) DOUBLE PRECISION VAR1,VAR2,P(4,MAXNUMOFPAR),PREF(4,MAXNUMOFPAR) & ,Y0(25),Y1(25) CHARACTER*100 HISNAM EXTERNAL HEXIST EXTERNAL HI,HSTATI,HSUM IF (USENUMOFENTRIES.LE.0) USENUMOFENTRIES=NEVT NPARAMS=0 DO I=1,NJET IF (PID(I).GT.0.AND.PARAMFLAG(PID(I)).GE.0) NPARAMS=NPARAMS & +NPAR(PARAMFLAG(PID(I))) ENDDO CALL VZERO(BOUND,2*5*4*MAXNUMOFPAR*400) CALL VZERO(SUM,2*25*2) CALL VZERO(COR,NPARAMS*NPARAMS*2) NHIST=0 HISNAM='Parameter histogram' HNLEN=19 IFLAG=MIN(5,MAX(2,IFLAGIN)) DO H=1,IFLAG DO I=1,USENUMOFENTRIES CALL NTREAD(I,ISTAT) IF (ISTAT.NE.0) THEN PRINT*,'Error reading entry: ',I,' from inputfile' PRINT*,'Continuing with ',I-1,' entries!!!' USENUMOFENTRIES=I-1 GOTO 900 ENDIF IF (MATCHFLAG.GT.0) THEN IPARAMS=0 DO J=1,NJET IF (PID(J).GT.IPARAMS) IPARAMS=PID(J) ENDDO DO J=1,IPARAMS IF (PARAMFLAG(J).GE.0) CALL & MATCH(MATCHFLAG,J) ENDDO ENDIF IF (CUT.LT.EVENTCUT.OR.EVFLAG.NE.EVENTFLAG) GOTO 800 DO J=1,NJET DO K=1,4 IF (BINLEVEL.EQ.1) P(K,J)=DBLE(PTRUE(K,J)) IF (BINLEVEL.EQ.2) P(K,J)=DBLE(PRECO(K,J)) PREF(K,J)=DBLE(PTRUE(K,J)) ENDDO IPAR(J)=-1 IF (PID(J).GT.0) IPAR(J)=PARAMFLAG(PID(J)) ENDDO CALL VZERO(Y1,25*2) CALL DEFUVEC(IPAR,NJET,PRECO,ISTAT) DO J=1,NJET IPAR(J)=-1 ENDDO IF (ISTAT.NE.0) THEN PRINT*,'Error calculating parameterisation reference' PRINT*,'system for event: ',I,' -- skipping event!!' GOTO 800 ENDIF IPARAMS=1 DO J=1,NJET IF (PID(J).LE.0) GOTO 700 IF (PARAMFLAG(PID(J)).LT.0) GOTO 700 IPAR(J)=PARAMFLAG(PID(J)) CALL CVP2Y(NJET,IPAR,PREF,Y0,ISTAT) C PRINT*,'J,IPAR(J)',J,IPAR(J),Y0 IPAR(J)=-1 IF (ISTAT.NE.0) THEN c IF (H.EQ.IFLAG) THEN c PRINT*,'Error calculating parameterisation' c PRINT*,'for particle: ',J,' in event: ',I c & ,' -- skipping particle!!' c ENDIFprint y0 GOTO 800 ENDIF VAR1=P(4,J) VAR2=SQRT(P(1,J)**2+P(2,J)**2+P(3,J)**2) IF (VAR2.GT.0.0D0) THEN VAR2=P(3,J)/VAR2 ELSE PRINT*,'Warning: Jet with P=0!!!(evt=',I,')' VAR2=-2.0D0 ENDIF CALL GETBIN(PID(J),VAR1,VAR2,IBIN) IF (IBIN.GT.0) THEN DO K=1,NPAR(PARAMFLAG(PID(J))) IF (H.EQ.1) THEN BOUND(1,K,PID(J),IBIN)=BOUND(1,K,PID(J),IBIN)+1.0D0 IF (BOUND(2,K,PID(J),IBIN).EQ.0.0D0) BOUND(2,K,PID(J) !HOLDS MINIMUM & ,IBIN)=1.0D32 IF (BOUND(3,K,PID(J),IBIN).EQ.0.0D0) BOUND(3,K,PID(J) !HOLDS MAXIMUM & ,IBIN)=-1.0D32 IF (Y0(K).LT.BOUND(2,K,PID(J) & ,IBIN)) BOUND(2,K,PID(J),IBIN)=Y0(K) IF (Y0(K).GT.BOUND(3,K,PID(J) & ,IBIN)) BOUND(3,K,PID(J),IBIN)=Y0(K) BOUND(4,K,PID(J),IBIN)=BOUND(4,K,PID(J),IBIN)+Y0(K) BOUND(5,K,PID(J),IBIN)=BOUND(5,K,PID(J),IBIN)+Y0(K) & *Y0(K) ELSEIF (H.EQ.2) THEN N=0 DO L=2,PID(J) IF (PARAMFLAG(L-1).GE.0) N=N+NPAR(PARAMFLAG(L-1)) ENDDO IHIST=IBIN+((K-1)+N)*MAX(NBIN1,1)*MAX(NBIN2,1) IF (.NOT.HEXIST(IHIST)) THEN IF (BOUND(1,K,PID(J),IBIN).GE.3.0D0) THEN RMEAN=BOUND(4,K,PID(J),IBIN)/BOUND(1,K,PID(J),IBIN & ) UPPER=SQRT(REAL(BOUND(5,K,PID(J),IBIN)/BOUND(1,K & ,PID(J),IBIN))-RMEAN*RMEAN) IF (ABS(RMEAN-RANGESCALECUT*UPPER).LT & .ABS(REAL(BOUND(2,K,PID(J),IBIN)))) THEN LOWER=RMEAN-RANGESCALECUT*UPPER ELSE LOWER=REAL(BOUND(2,K,PID(J),IBIN)) ENDIF IF (ABS(RMEAN+RANGESCALECUT*UPPER).LT & .ABS(REAL(BOUND(3,K,PID(J),IBIN)))) THEN UPPER=RMEAN+RANGESCALECUT*UPPER ELSE UPPER=REAL(BOUND(3,K,PID(J),IBIN)) ENDIF UPPER=MIN(ABS(RMEAN-LOWER),ABS(RMEAN-UPPER)) LOWER=RMEAN-UPPER UPPER=RMEAN+UPPER RMEAN=SQRT((REAL(BOUND(5,K,PID(J),IBIN)/BOUND(1,K & ,PID(J),IBIN))-RMEAN*RMEAN)/REAL(BOUND(1,K,PID(J & ),IBIN))) ELSE LOWER=-10.0 UPPER=10.0 RMEAN=1.0 ENDIF IF (H.EQ.IFLAG) THEN CALL GETLIM(PID(J),VAR1,VAR2,ELOW,EUP,CTHLOW,CTHUP & ) CALL GENHISNAM(PID(J),K,ELOW,EUP,CTHLOW,CTHUP & ,HNLEN,HISNAM) NHIST=NHIST+1 IF (URANGE(IHIST,1).LT.URANGE(IHIST,2)) THEN LOWER=URANGE(IHIST,1) UPPER=URANGE(IHIST,2) PRINT* & ,'Using user specified range for histogram: ' & ,IHIST ENDIF CALL HBOOK1(IHIST,HISNAM(1:HNLEN),NBINFIT & ,LOWER,UPPER,0.0) ELSE L=INT((UPPER-LOWER)/RMEAN) CALL HBOOK1(IHIST,HISNAM(1:HNLEN),L & ,LOWER,UPPER,0.0) BOUND(5,K,PID(J),IBIN)=DBLE(L) ENDIF CALL HIDOPT(IHIST,'STAT') ENDIF CALL HFILL(IHIST,REAL(Y0(K)),0.0,1.0) ELSEIF (H.EQ.3) THEN N=0 DO L=2,PID(J) IF (PARAMFLAG(L-1).GE.0) N=N+NPAR(PARAMFLAG(L-1)) ENDDO IHIST=IBIN+((K-1)+N)*MAX(NBIN1,1)*MAX(NBIN2,1) IF (BOUND(1,K,PID(J),IBIN).GT.2.0D0) THEN UPPER=0.0 RMEAN=0.0 BINMAX=0.0 NBINOK=0 DO N=1,INT(BOUND(5,K,PID(J),IBIN)) IF (HI(IHIST,N).GT.0.0) NBINOK=NBINOK+1 IF (HI(IHIST,N).GT.UPPER) THEN UPPER=HI(IHIST,N) CALL HIX(IHIST,N,BINMAX) ENDIF ENDDO CALL HIX(IHIST,1,ELOW) CALL HIX(IHIST,2,EUP) ELOW=EUP-ELOW LOWER=0.0 L=1 N=INT(HSUM(IHIST)) DO WHILE((LOWER+HI(IHIST,L)/REAL(N)).LT.(SCALEFRAC/2 & .0).AND.L.LT.INT(BOUND(5,K,PID(J),IBIN))) LOWER=LOWER+HI(IHIST,L)/REAL(N) L=L+1 ENDDO CALL HIX(IHIST,L,LOWER) UPPER=0.0 M=1 DO WHILE ((UPPER+HI(IHIST,M)/REAL(N)).LT.(1.0 & -SCALEFRAC/2.0).AND.M.LT.INT(BOUND(5,K,PID(J),IBIN & ))) UPPER=UPPER+HI(IHIST,M)/REAL(N) M=M+1 ENDDO CALL HIX(IHIST,M,UPPER) RMEAN=HSTATI(IHIST,1,' ',1) IF (M.GT.L) THEN BOUND(3,K,PID(J),IBIN)=DBLE(RMEAN+MAX(ABS(RMEAN & -LOWER),ABS(UPPER-RMEAN))) ELSE BOUND(3,K,PID(J),IBIN)=DBLE(RMEAN+MAX(ABS(RMEAN & -LOWER),ABS(UPPER-RMEAN))) ENDIF UPPER=RMEAN+ELOW*REAL(NBINOK)/2.0 IF (ABS(UPPER).LT.ABS(REAL(BOUND(3,K,PID(J),IBIN)))) & THEN BOUND(3,K,PID(J),IBIN)=DBLE(UPPER) ENDIF BOUND(2,K,PID(J),IBIN)=DBLE(RMEAN-ABS(BOUND(3,K & ,PID(J),IBIN)-RMEAN)) CALL HDELET(IHIST) IF (H.EQ.IFLAG) THEN CALL GETLIM(PID(J),VAR1,VAR2,ELOW,EUP,CTHLOW,CTHUP & ) CALL GENHISNAM(PID(J),K,ELOW,EUP,CTHLOW,CTHUP & ,HNLEN,HISNAM) IF (URANGE(IHIST,1).LT.URANGE(IHIST,2)) THEN BOUND(2,K,PID(J),IBIN)=DBLE(URANGE(IHIST,1)) BOUND(3,K,PID(J),IBIN)=DBLE(URANGE(IHIST,2)) PRINT* & ,'Using user specified range for histogram: ' & ,IHIST ENDIF ENDIF CALL HBOOK1(IHIST,HISNAM(1:HNLEN),NBINFIT & ,REAL(BOUND(2,K,PID(J),IBIN)),REAL(BOUND(3,K,PID(J & ),IBIN)),0.0) NHIST=NHIST+1 CALL HIDOPT(IHIST,'STAT') BOUND(1,K,PID(J),IBIN)=0.0D0 BOUND(1,1,1,0)=0.0D0 ELSEIF (BOUND(1,K,PID(J),IBIN).GT.0.0D0) THEN nhist=nhist+1 call hreset(ihist,' ') BOUND(1,K,PID(J),IBIN)=0.0D0 ENDIF CALL HFILL(IHIST,REAL(Y0(K)),0.0,1.0) ELSEIF (H.EQ.4) THEN C C FIT ALL HISTO ONCE TO GET REASONABLE LIMITS C IF (BOUND(1,1,1,0).LE.0.0D0) THEN BOUND(1,1,1,0)=1.0D0 ! just a tag-flag PRINT*,'---- Doing 1st range-fits ----' CALL FITPARAM(2) PRINT*,'---- 1st range-fits done ----' DO L=1,NHIST IF (HEXIST(L)) CALL HDELET(L) ENDDO ENDIF N=0 DO L=2,PID(J) IF (PARAMFLAG(L-1).GE.0) N=N+NPAR(PARAMFLAG(L-1)) ENDDO IHIST=IBIN+((K-1)+N)*MAX(NBIN1,1)*MAX(NBIN2,1) IF (.NOT.HEXIST(IHIST)) THEN IF (FITOK(IHIST).LE.2) THEN BOUND(2,K,PID(J),IBIN)=DBLE(FITMEAN(IHIST,1 & )-RANGESCALECUT*FITSIGMA(IHIST,1)) BOUND(3,K,PID(J),IBIN)=DBLE(FITMEAN(IHIST,1 & )+RANGESCALECUT*FITSIGMA(IHIST,1)) ELSE BOUND(2,K,PID(J),IBIN)=-10.0D0 BOUND(3,K,PID(J),IBIN)=10.0D0 ENDIF IF (H.EQ.IFLAG) THEN CALL GETLIM(PID(J),VAR1,VAR2,ELOW,EUP,CTHLOW,CTHUP & ) CALL GENHISNAM(PID(J),K,ELOW,EUP,CTHLOW,CTHUP & ,HNLEN,HISNAM) IF (URANGE(IHIST,1).LT.URANGE(IHIST,2)) THEN BOUND(2,K,PID(J),IBIN)=DBLE(URANGE(IHIST,1)) BOUND(3,K,PID(J),IBIN)=DBLE(URANGE(IHIST,2)) PRINT* & ,'Using user specified range for histogram: ' & ,IHIST ENDIF ENDIF CALL HBOOK1(IHIST,HISNAM(1:HNLEN),NBINFIT & ,REAL(BOUND(2,K,PID(J),IBIN)),REAL(BOUND(3,K & ,PID(J),IBIN)),0.0) CALL HIDOPT(IHIST,'STAT') ENDIF CALL HFILL(IHIST,REAL(Y0(K)),0.0,1.0) ELSEIF (H.EQ.5) THEN C C REDO ALL BAD FIT TO GET REASONABLE LIMITS C IF (BOUND(1,1,1,0).EQ.1.0D0) THEN BOUND(1,1,1,0)=0.0D0 ! just a tag-flag PRINT*,'---- Doing 2nd range-fits ----' CALL FITPARAM(3) ! re-fit only failed fits PRINT*,'---- 2nd range-fits done ----' DO L=1,NHIST IF (HEXIST(L).AND.FITOK(L).LE.2) THEN FITMEAN(L,2)=HSTATI(L,1,' ',1) FITSIGMA(L,2)=HSTATI(L,2,' ',1) CALL HDELET(L) ELSEIF (HEXIST(L)) THEN CALL HRESET(L,' ') BOUND(2,K,PID(J),IBIN)=-10.0D0 BOUND(3,K,PID(J),IBIN)=10.0D0 ENDIF ENDDO ENDIF N=0 DO L=2,PID(J) IF (PARAMFLAG(L-1).GE.0) N=N+NPAR(PARAMFLAG(L-1)) ENDDO IHIST=IBIN+((K-1)+N)*MAX(NBIN1,1)*MAX(NBIN2,1) IF (.NOT.HEXIST(IHIST)) THEN IF (FITOK(IHIST).LE.2) THEN BOUND(2,K,PID(J),IBIN)=DBLE(FITMEAN(IHIST,1 & )-RANGESCALECUT*FITSIGMA(IHIST,1)) BOUND(3,K,PID(J),IBIN)=DBLE(FITMEAN(IHIST,1 & )+RANGESCALECUT*FITSIGMA(IHIST,1)) IF (ABS(DBLE(FITMEAN(IHIST,2)-RANGESCALECUT & *FITSIGMA(IHIST,2))).LT.1.5D0*ABS(BOUND(2,K & ,PID(J),IBIN))) THEN BOUND(2,K,PID(J),IBIN)=DBLE(FITMEAN(IHIST,2) & -RANGESCALECUT*FITSIGMA(IHIST,2)) BOUND(3,K,PID(J),IBIN)=DBLE(FITMEAN(IHIST,2) & +RANGESCALECUT*FITSIGMA(IHIST,2)) ENDIF ENDIF IF (H.EQ.IFLAG) THEN CALL GETLIM(PID(J),VAR1,VAR2,ELOW,EUP,CTHLOW,CTHUP & ) CALL GENHISNAM(PID(J),K,ELOW,EUP,CTHLOW,CTHUP & ,HNLEN,HISNAM) IF (URANGE(IHIST,1).LT.URANGE(IHIST,2)) THEN BOUND(2,K,PID(J),IBIN)=DBLE(URANGE(IHIST,1)) BOUND(3,K,PID(J),IBIN)=DBLE(URANGE(IHIST,2)) PRINT* & ,'Using user specified range for histogram: ' & ,IHIST ENDIF ENDIF CALL HBOOK1(IHIST,HISNAM(1:HNLEN),NBINFIT & ,REAL(BOUND(2,K,PID(J),IBIN)),REAL(BOUND(3,K & ,PID(J),IBIN)),0.0) CALL HIDOPT(IHIST,'STAT') ENDIF CALL HFILL(IHIST,REAL(Y0(K)),0.0,1.0) ENDIF ENDDO ENDIF DO K=1,NPAR(PARAMFLAG(PID(J))) Y1(IPARAMS-1+K)=Y0(K) ENDDO IPARAMS=IPARAMS+NPAR(PARAMFLAG(PID(J))) 700 CONTINUE ENDDO IF (H.EQ.IFLAG) THEN IPARAMS=1 DO L=1,NJET IF (PID(L).GT.0.AND.PARAMFLAG(PID(L)).GE.0.AND.(Y1(IPARAMS & )+Y1(IPARAMS+1)+Y1(IPARAMS+2)).NE.0.0D0) THEN DO K=0,NPAR(PARAMFLAG(PID(L)))-1 SUM(1,IPARAMS+K)=SUM(1,IPARAMS+K)+1.0D0 SUM(2,IPARAMS+K)=SUM(2,IPARAMS+K)+Y1(IPARAMS+K) ENDDO J=1 DO M=1,L IF (PID(M).GT.0.AND.PARAMFLAG(PID(M)).GE.0.AND.(Y1(J & )+Y1(J+1)+Y1(J+2)).NE.0.0D0) THEN DO K=0,NPAR(PARAMFLAG(PID(L)))-1 DO N=0,NPAR(PARAMFLAG(PID(M)))-1 COR(IPARAMS+K,J+N)=COR(IPARAMS+K,J+N)+Y1(IPARAMS & +K)*Y1(J+N) ENDDO ENDDO J=J+NPAR(PARAMFLAG(PID(M))) ENDIF ENDDO IPARAMS=IPARAMS+NPAR(PARAMFLAG(PID(L))) ENDIF ENDDO ENDIF 800 CONTINUE ENDDO 900 CONTINUE ENDDO C C Construct correlation matrix C IPARAMS=1 DO L=1,NJET IF (PID(L).GT.0.AND.PARAMFLAG(PID(L)).GE.0) THEN J=1 DO M=1,L IF (PID(M).GT.0.AND.PARAMFLAG(PID(M)).GE.0) THEN DO K=0,NPAR(PARAMFLAG(PID(L)))-1 DO N=0,NPAR(PARAMFLAG(PID(M)))-1 IF (IPARAMS+K.NE.J+N) THEN COR(IPARAMS+K,J+N)=(COR(IPARAMS+K,J+N)-SUM(2,IPARAMS & +K)*SUM(2,J+N)/SUM(1,IPARAMS+K))/SQRT(COR(IPARAMS & +K,IPARAMS+K)-SUM(2,IPARAMS+K)**2/SUM(1,IPARAMS+K) & )/SQRT(COR(J+N,J+N)-SUM(2,J+N)**2/SUM(1,J+N)) COR(J+N,IPARAMS+K)=COR(IPARAMS+K,J+N) ENDIF ENDDO ENDDO J=J+NPAR(PARAMFLAG(PID(M))) ENDIF ENDDO IPARAMS=IPARAMS+NPAR(PARAMFLAG(PID(L))) ENDIF ENDDO NPARAMS=0 DO I=1,NJET IF (PID(I).GT.0.AND.PARAMFLAG(PID(I)).GE.0) NPARAMS=NPARAMS & +NPAR(PARAMFLAG(PID(I))) ENDDO DO I=1,NPARAMS COR(I,I)=1.0D0 ENDDO 999 RETURN END SUBROUTINE GETBIN(CLFLAG,VAR1,VAR2,IBIN) IMPLICIT NONE INTEGER CLFLAG,IBIN DOUBLE PRECISION VAR1,VAR2 #include "make_param.inc" REAL RTEMP IBIN=0 IF (NBIN1.GT.0) THEN RTEMP = (B1MAX(CLFLAG)-B1MIN(CLFLAG))/REAL(NBIN1) IF (REAL(VAR1)-B1MIN(CLFLAG).GE.0.0.AND.REAL(VAR1)-B1MAX(CLFLAG) & .LE.0.0) IBIN = INT((REAL(VAR1)-B1MIN(CLFLAG))/RTEMP)+1 ENDIF IF (NBIN2.GT.0) THEN RTEMP = (B2MAX(CLFLAG)-B2MIN(CLFLAG))/REAL(NBIN2) IF (REAL(VAR2)-B2MIN(CLFLAG).GE.0.0.AND.REAL(VAR2)-B2MAX(CLFLAG) & .LE.0.0) THEN IF (IBIN.GT.0) THEN IBIN =IBIN+INT((REAL(VAR2)-B2MIN(CLFLAG))/RTEMP)*NBIN1 ELSE IF (NBIN1.EQ.0) IBIN = INT((REAL(VAR2)-B2MIN(CLFLAG))/RTEMP) & +1 ENDIF ELSE IBIN=0 ENDIF ENDIF 999 RETURN END C====================================================================== C Routine to fit all the histograms according to choice of fitting function C====================================================================== SUBROUTINE FITPARAM(IFLAGIN) IMPLICIT NONE INTEGER IFLAGIN #include "make_param.inc" LOGICAL HEXIST,WROUT INTEGER NY,NWT,LOC,IERFLAG,I,ISTAT,NPARI,NPARX,FITSTAT,SAVEFITFUNC & ,IFLAG REAL XMI,XMA,YMI,YMA,HSTATI,HI,RDLHOOD,HMAX,rtemp1(5),rtemp2(5) DOUBLE PRECISION DSTVAL(5),DMEAN_FIT,DMEAN_UNC & ,DSIGMA_FIT,DSIGMA_UNC,FCN,DBNDL(5),DBNDU(5),DSTEP(5) & ,DARGLIS(10),FMIN,DGAUSS,DLHOOD,FEDM,ERRDEF,FMEAN CHARACTER*10 CHNAM(5) CHARACTER*80 TITL DATA CHNAM /'mean ','sigma','dummy','const','slope'/ EXTERNAL FCN,HSTATI,HI,DGAUSS,DLHOOD,HEXIST,RDLHOOD,HMAX SAVE CHNAM IFLAG=IFLAGIN WROUT=.TRUE. SAVEFITFUNC=FITFUNC IF (IFLAG.GT.1) THEN WROUT=.FALSE. FITFUNC=1 ENDIF CALL MNINIT(5,6,7) DO HISTOFIT=1,NHIST IF (IFLAGIN.GT.2.AND.FITOK(HISTOFIT).NE.2) GOTO 900 FITMEAN(HISTOFIT,1)=0.0 FITSIGMA(HISTOFIT,1)=0.0 FITMEAN(HISTOFIT,2)=0.0 FITSIGMA(HISTOFIT,2)=0.0 NEVNT=0.0D0 IF (HEXIST(HISTOFIT)) NEVNT=DBLE(HSTATI(HISTOFIT,3,' ',1)) IF (.NOT.HEXIST(HISTOFIT)) GOTO 900 CALL HGIVE(HISTOFIT,TITL,OLI_NX,XMI,XMA,NY,YMI,YMA,NWT,LOC) IF (NEVNT.LE.2.D0) THEN IF (WROUT) THEN PRINT*,'***********************************************' PRINT*,'FIT FOR HISTOGRAM', HISTOFIT,' IS NOT POSSIBLE' PRINT*,'ONLY ',INT(NEVNT),' ENTRIES' PRINT*,'USE MORE EVENTS OR INCREASE SIZE OF THE BINS' PRINT*,'***********************************************' ENDIF FITMEAN(HISTOFIT,1)=(XMA-XMI)/2.0 FITSIGMA(HISTOFIT,1)=XMA-FITMEAN(HISTOFIT,1) FITMEAN(HISTOFIT,2)=MAX(10000.0,10.0*FITMEAN(HISTOFIT,1)) FITSIGMA(HISTOFIT,2)=MAX(10000.0,10.0*FITSIGMA(HISTOFIT,1)) FITOK(HISTOFIT)=3 GOTO 900 ELSE FITMEAN(HISTOFIT,1)=HSTATI(HISTOFIT,1,' ',0) FITSIGMA(HISTOFIT,1)=HSTATI(HISTOFIT,2,' ',0) FITMEAN(HISTOFIT,2)=FITSIGMA(HISTOFIT,1)*2./SQRT(REAL(NEVNT)) FITSIGMA(HISTOFIT,2)=FITSIGMA(HISTOFIT,1)*3./SQRT(REAL(NEVNT)) IF (NEVNT.LT.50.0D0) THEN IF (WROUT) THEN PRINT*,'***********************************************' PRINT*,'FIT FOR HISTOGRAM', HISTOFIT,' IS NOT POSSIBLE' PRINT*,'ONLY ',INT(NEVNT),' ENTRIES' PRINT*,'USE MEAN AND RMS FROM HISTOGRAM' PRINT*,'***********************************************' PRINT*,'MEAN=',FITMEAN(HISTOFIT,1), '+-',FITMEAN(HISTOFIT & ,2) PRINT*,'RMS =',FITSIGMA(HISTOFIT,1),'+-',FITSIGMA(HISTOFIT & ,2) PRINT*,'***********************************************' ENDIF FITOK(HISTOFIT)=1 GOTO 900 ENDIF DSTVAL(1) = DBLE(FITMEAN(HISTOFIT,1)) DSTVAL(2) = DBLE(FITSIGMA(HISTOFIT,1)) ENDIF NY=1 OLI_EPS=.001D0 OLI_A=DBLE(XMI) OLI_B=DBLE(XMA) C COMPUTE LIKELIHOOD NORM DMEAN=DSTVAL(1) DSIGMA=DSTVAL(2) DDUMMY=1.0D0 DCONST=0D0 IF (FITFUNC.EQ.2) DCONST=1D0 DSLOPE=0.0D0 OLI_WX = (XMA-XMI)/REAL(OLI_NX) DSTVAL(3) = DGAUSS(DLHOOD,OLI_A,OLI_B,OLI_EPS) DSTVAL(3) = NEVNT*DBLE(OLI_WX)/DSTVAL(3) DSTVAL(4) = 0D0 IF (FITFUNC.GE.2) DSTVAL(4) = 1D0 DSTVAL(5) = 0D0 C C SET THE PARAMETER VALUES C DBNDL(1) = -10D0 DBNDU(1) = 10D0 DBNDL(2) = 0D0 DBNDU(2) = 100D0 DBNDL(3) = 0D0 DBNDU(3) = 1000000D0 DBNDL(4) = 0D0 DBNDU(4) = 10000D0 DBNDL(5) = -1000D0 DBNDU(5) = 1000D0 DO I = 1,5 DSTEP(I) = 0.001D0 ENDDO 800 CONTINUE FITOK(HISTOFIT)=0 LLIM=-1.0E32 ULIM=1.0E32 IF (IFLAG.GT.2) THEN DSTVAL(2)=DBLE(OLI_WX) LLIM=REAL(DSTVAL(1))-MAX(3.0*OLI_WX,FITSIGMA(HISTOFIT,1)) ULIM=REAL(DSTVAL(1))+MAX(3.0*OLI_WX,FITSIGMA(HISTOFIT,1)) ENDIF DARGLIS(1) = -1.0D0 CALL MNEXCM (FCN,'SET PRINTOUT ',DARGLIS,1,IERFLAG & ) DO I = 1,3 CALL MNPARM (I,CHNAM(I),DSTVAL(I),DSTEP(I),DBNDL(I) & ,DBNDU(I),IERFLAG) IF (IERFLAG.NE.0) THEN IF (WROUT) PRINT*,'UNABLE TO DEFINE PARAMETER NUMBER',I & ,IERFLAG ENDIF ENDDO IF (FITFUNC.GE.2) THEN I=4 CALL MNPARM (I,CHNAM(I),DSTVAL(I),DSTEP(I),DBNDL(I) & ,DBNDU(I),IERFLAG) IF (IERFLAG.NE.0) THEN IF (WROUT) PRINT*,'UNABLE TO DEFINE PARAMETER NUMBER',I & ,IERFLAG ENDIF ENDIF IF (FITFUNC.EQ.3) THEN I=5 CALL MNPARM (I,CHNAM(I),DSTVAL(I),DSTEP(I),DBNDL(I) & ,DBNDU(I),IERFLAG) IF (IERFLAG.NE.0) THEN IF (WROUT) PRINT*,'UNABLE TO DEFINE PARAMETER NUMBER',I & ,IERFLAG ENDIF DARGLIS(1) = 4.0D0 CALL MNEXCM (FCN,'FIX ',DARGLIS,1 & ,IERFLAG) DARGLIS(1) = 5.0D0 CALL MNEXCM (FCN,'FIX ',DARGLIS,1 & ,IERFLAG) ENDIF DARGLIS(1) = 1.0D0 CALL MNEXCM (FCN,'CALL FCN ',DARGLIS,1,IERFLAG) DARGLIS(1) = 0.0D0 CALL MNEXCM (FCN,'MIGRAD ',DARGLIS,0,IERFLAG) IF (FITFUNC.EQ.3) THEN DARGLIS(1) = 4.0D0 CALL MNEXCM (FCN,'RELEASE ',DARGLIS,1 & ,IERFLAG) DARGLIS(1) = 0.0D0 CALL MNEXCM (FCN,'MIGRAD ',DARGLIS,0,IERFLAG & ) DARGLIS(1) = 5.0D0 CALL MNEXCM (FCN,'RELEASE ',DARGLIS,1 & ,IERFLAG) CALL MNEXCM (FCN,'MIGRAD ',DARGLIS,0,IERFLAG & ) ENDIF CALL MNSTAT (FMIN,FEDM,ERRDEF,NPARI,NPARX,ISTAT) IF (FITFUNC.EQ.3.AND.(ISTAT.NE.3.OR.FMIN.LE.0.000001D0)) THEN FITFUNC=2 CALL MNEXCM (FCN,'CLEAR ',0,0,IERFLAG) GOTO 800 ENDIF IF (IERFLAG.EQ.0.AND.ISTAT.EQ.3.AND.FMIN.GT.0.001D0) THEN C COPY THE FITTED VALUES OF THE MEAN AND SIGMA AND THEIR C THEIR ERRORS CALL MNPOUT(1,CHNAM(1),DMEAN_FIT,DMEAN_UNC,DARGLIS(9) & ,DARGLIS(10),ISTAT) CALL MNPOUT(2,CHNAM(2),DSIGMA_FIT,DSIGMA_UNC,DARGLIS(9) & ,DARGLIS(10),ISTAT) CALL MNPOUT(3,CHNAM(3),DARGLIS(6),DARGLIS(8),DARGLIS(9) & ,DARGLIS(10),ISTAT) FITMEAN(HISTOFIT,1) = REAL(DMEAN_FIT) FITMEAN(HISTOFIT,2) = REAL(DMEAN_UNC) FITSIGMA(HISTOFIT,1) = REAL(DSIGMA_FIT) FITSIGMA(HISTOFIT,2) = REAL(DSIGMA_UNC) FITHEIGHT(HISTOFIT,1) = REAL(DARGLIS(6)) FITHEIGHT(HISTOFIT,2) = REAL(DARGLIS(8)) FITCHI2(HISTOFIT)=REAL(FMIN)/REAL(OLI_NBIN-NPARI) IF (WROUT.OR.IFLAG.GT.2) THEN PRINT*,'Result for histogram ',HISTOFIT,' (fmin=',FMIN,'):' PRINT*,' MEAN=',FITMEAN(HISTOFIT,1), '+-' & ,FITMEAN(HISTOFIT,2) PRINT*,' RMS =',FITSIGMA(HISTOFIT,1),'+-' & ,FITSIGMA(HISTOFIT,2) ENDIF DMEAN = DMEAN_FIT DSIGMA = DSIGMA_FIT DDUMMY = DARGLIS(6) DCONST = 0.0D0 DSLOPE = 0.0D0 IF (FITFUNC.GE.2) THEN CALL MNPOUT(4,CHNAM(4),DCONST,DARGLIS(8) & ,DARGLIS(9),DARGLIS(10),ISTAT) FITHEIGHT(HISTOFIT,1) = FITHEIGHT(HISTOFIT,1)+REAL(DCONST) FITHEIGHT(HISTOFIT,2) = SQRT(REAL(DARGLIS(8))**2 & +FITHEIGHT(HISTOFIT,2)**2) ENDIF IF (FITFUNC.EQ.3) THEN CALL MNPOUT(5,CHNAM(5),DSLOPE,DARGLIS(8) & ,DARGLIS(9),DARGLIS(10),ISTAT) ENDIF IF (WROUT) THEN YMA=REAL(HISTOFIT) CALL CVR2C(YMA,0,I,TITL) CALL HBFUN1 (MAXNUMOFHIS+HISTOFIT & ,'Fit function for histogram '//TITL(1:I),MAX(100,4 & *NBINFIT),XMI,XMA,RDLHOOD) ENDIF CALL MNEXCM (FCN,'CLEAR ',0,0,IERFLAG) ELSE PRINT*,'Fit to histogram ',HISTOFIT,' FAILED! Using:' PRINT*,' MEAN=',FITMEAN(HISTOFIT,1), '+-' & ,FITMEAN(HISTOFIT,2) PRINT*,' RMS =',FITSIGMA(HISTOFIT,1),'+-' & ,FITSIGMA(HISTOFIT,2) FITOK(HISTOFIT)=2 ENDIF IF (SAVEFITFUNC.EQ.3.AND.FITFUNC.EQ.2) FITFUNC=3 IF (IFLAG.GT.1.AND.NY.GT.0.AND.(FITOK(HISTOFIT).EQ.2.OR & .(FITOK(HISTOFIT).EQ.0.AND.ABS(REAL(DDUMMY)-HMAX(HISTOFIT)).GT & .7.0*SQRT(HMAX(HISTOFIT)+FITHEIGHT(HISTOFIT,2)**2)))) THEN IF (IFLAG.GT.2) THEN PRINT*,'Congratulations -- You`ve got a REALLY nasty fit' PRINT*,'Last retry before giving up for histogram: ' & ,HISTOFIT ENDIF YMA=0.0 NY=0 DO I=1,OLI_NX IF (HI(HISTOFIT,I).GT.YMA) THEN YMA=HI(HISTOFIT,I) NY=I ENDIF ENDDO IF (FITOK(HISTOFIT).EQ.0) DSTVAL(2)=DBLE(OLI_WX) IF (NY.GT.0) THEN DSTVAL(1)=OLI_A+OLI_WX*(DBLE(NY)-0.5D0) NY=0 IFLAG=3 GOTO 800 ELSE PRINT*,'This histogram makes no sense -- giving up' ENDIF ENDIF 900 CONTINUE ENDDO IF (WROUT) CALL CHECKFIT FITFUNC=SAVEFITFUNC 999 RETURN END SUBROUTINE CHECKFIT IMPLICIT NONE #include "make_param.inc" LOGICAL HEXIST INTEGER H,I,J,IHIST,NOENT REAL FITCHI,CHIMIN,CHIMAX,HI,HSTATI,WARNINGLIMIT,CHICUT,RBIN,RTEMP & ,HMAX,CHICUT1,HISFLOWLIM,HSUM,HISCHECK(10) CHARACTER*10 NTNAME(10) DATA NTNAME /'HISID','TAIL','FITOUTPUT','HEIGHT','HISFLOW' & ,'FITFAIL',' ',' ',' ',' '/ EXTERNAL HI,HEXIST,HSTATI,HMAX,HSUM IHIST=0 CHIMIN=1.0E16 CHIMAX=-1.0E16 WARNINGLIMIT=0.05 HISFLOWLIM=0.15 IF (FITFLAG.EQ.1) THEN CHICUT=4.0 CHICUT1=5.0 ENDIF IF (FITFLAG.EQ.2) THEN CHICUT=20.0 CHICUT1=3.0 ENDIF PRINT*,'List of possible problematic histograms' PRINT*,'PLEASE CHECK!!!!!!' PRINT*,'-------------START-OF-LIST--------------' CALL HBOOKN(3*MAXNUMOFHIS,'His check',10,'//OUTPUT',1024,NTNAME) DO H=1,3 DO I=1,NHIST CALL VZERO(HISCHECK,10) IF (FITOK(I).EQ.0) THEN IF (H.EQ.1) THEN IF (FITCHI2(I).LT.CHIMIN) CHIMIN=FITCHI2(I) IF (FITCHI2(I).GT.CHIMAX) CHIMAX=FITCHI2(I) IHIST=IHIST+1 ELSEIF (H.EQ.2) THEN IF (.NOT.HEXIST(NHIST+1)) CALL HBOOK1(NHIST+1,'Chi2 check' & ,50,CHIMIN,CHIMAX,0.0) CALL HFILL(NHIST+1,FITCHI2(I),0.0,1.0/REAL(IHIST)) IF (I.EQ.NHIST) THEN ! FINAL HIST-VALUES HAS BEEN FILLED CHIMIN=0.0 J=1 DO WHILE((CHIMIN+HI(NHIST+1,J)).LT.WARNINGLIMIT) CHIMIN=CHIMIN+HI(NHIST+1,J) J=J+1 ENDDO CALL HIX(NHIST+1,J,CHIMIN) CHIMAX=0.0 J=1 DO WHILE(CHIMAX+HI(NHIST+1,J).LT.1.0-WARNINGLIMIT) CHIMAX=CHIMAX+HI(NHIST+1,J) J=J+1 ENDDO CALL HIX(NHIST+1,J+1,CHIMAX) ENDIF ELSEIF (H.EQ.3) THEN C C NOW CALCULATE A LOT OF HISTGRAM CHECKS C CALL HIX(I,1,RBIN) CALL HIX(I,2,RTEMP) RBIN=RTEMP-RBIN FITCHI=(FITMEAN(I,1)-HSTATI(I,1,' ',1))**2/(FITMEAN(I,2) & **2+HSTATI(I,2,' ',1)**2/HSTATI(I,3,' ',1))+(FITSIGMA(I & ,1)-HSTATI(I,2,' ',1))**2/ & (FITSIGMA(I,2)**2+HSTATI(I,2,' ',1)**2/HSTATI(I & ,3,' ',1)) CALL HNOENT(I,NOENT) IF (FITCHI2(I).GT.CHIMAX.OR.FITCHI/2.0.GT.CHICUT.OR & .ABS(HMAX(I)-HMAX(MAXNUMOFHIS+I))/SQRT(HMAX(I) & +FITHEIGHT(I,2)**2).GT.CHICUT1.OR.(REAL(NOENT)/HSUM(I)-1 & .0).GT.HISFLOWLIM)THEN HISCHECK(1)=REAL(I) PRINT'(1x,a,i4)','WARNING for histogram: ',I IF (FITCHI2(I).GT.CHIMAX) THEN PRINT'(1x,a,f6.2,a)' & ,' Fit-Chi2 lies in the UPPER ',100 & *WARNINGLIMIT,'% tail' HISCHECK(2)=1.0 ENDIF IF (FITCHI/2.0.GT.CHICUT) THEN PRINT'(1x,a,f6.2,a)' & ,' Fit-output is more than ',SQRT(CHICUT) & ,' sigmas away from RMS or MEAN of histogram.' HISCHECK(3)=1.0 ENDIF IF (ABS(HMAX(I)-HMAX(MAXNUMOFHIS+I)) & /SQRT(HMAX(I)+FITHEIGHT(I,2)**2).GT.CHICUT1) THEN PRINT'(1x,a,f6.2,a)' & ,' Fit Max height is more than ',CHICUT1 & ,' sigma away for histogram.' HISCHECK(4)=1.0 ENDIF IF ((REAL(NOENT)/HSUM(I)-1.0).GT.HISFLOWLIM) THEN PRINT'(1x,a,f6.2,a)',' More than ',100.0 & *HISFLOWLIM,'% under/overflow in histogram.' HISCHECK(5)=1.0 ENDIF ENDIF ENDIF ELSE PRINT'(1x,a,i4)','WARNING for histogram: ',I PRINT'(1x,a)',' Histogram fit failed!!!! ' HISCHECK(1)=REAL(I) HISCHECK(6)=1.0 ENDIF IF (HISCHECK(1).GT.0.0) CALL HFN(3*MAXNUMOFHIS,HISCHECK) ENDDO ENDDO PRINT*,'--------------END-OF-LIST---------------' 999 RETURN END C====================================================================== C Fit function C====================================================================== SUBROUTINE FCN (NPARAM,GRAD,FVAL,XVAL,IFLAG) IMPLICIT NONE INTEGER NPARAM,IFLAG DOUBLE PRECISION GRAD(*),FVAL,XVAL(*) #include "make_param.inc" INTEGER IBIN REAL HI,HIE DOUBLE PRECISION DL_CAL,PROB,XBIN,YBIN,FYBIN,DLHOOD EXTERNAL HI,HIE EXTERNAL DLHOOD C C TRANSFER VALUES TO GAUSS-FUNCTION C DMEAN = XVAL(1) DSIGMA = XVAL(2) DDUMMY = XVAL(3) DCONST = 0.0D0 DSLOPE = 0.0D0 IF (FITFUNC.GE.2) DCONST = XVAL(4) IF (FITFUNC.EQ.3) DSLOPE = XVAL(5) IF (DDUMMY .LE. 0D0) RETURN C C LOOP OVER ALL BINS C DL_CAL = 0D0 OLI_NBIN=0 DO IBIN=1,OLI_NX XBIN=OLI_A+DBLE(OLI_WX*(REAL(IBIN)-0.5)) YBIN=DBLE(HI(HISTOFIT,IBIN)) FYBIN=DBLE(HIE(HISTOFIT,IBIN)) C C SUM UP LIKELIHOOD C PROB = -10.0 IF (REAL(XBIN).GE.LLIM.AND.REAL(XBIN).LE.ULIM) PROB = & DLHOOD(XBIN) IF (FITFLAG.EQ.1) THEN IF (PROB.GT.0D0) THEN DL_CAL = DL_CAL+2D0*(PROB - YBIN*DLOG(PROB)) OLI_NBIN=OLI_NBIN+1 ENDIF ELSE IF (YBIN.GT.0.0.AND.FYBIN.GT.0.0.AND.PROB.GE.0.0D0) THEN DL_CAL=DL_CAL+((DBLE(YBIN)-PROB)/DBLE(FYBIN))**2 OLI_NBIN=OLI_NBIN+1 ENDIF ENDIF ENDDO FVAL = DL_CAL 999 RETURN END C====================================================================== C Misc. functions for the FIT-FUNCTION C====================================================================== DOUBLE PRECISION FUNCTION NUMBER5(X) IMPLICIT NONE DOUBLE PRECISION X #include "make_param.inc" DOUBLE PRECISION XDUMMY,XDUMMY1,XDUMMY2,XNEW,OLISGAUSS EXTERNAL OLISGAUSS XDUMMY=DSQRT(1D0/(27D0*DCUBE**3)+(X/2D0/DCUBE)**2) XDUMMY1=X/2D0/DCUBE+XDUMMY XDUMMY1=DSIGN(1D0,XDUMMY1)*(ABS(XDUMMY1))**(1D0/3D0) XDUMMY2=X/2D0/DCUBE-XDUMMY XDUMMY2=DSIGN(1D0,XDUMMY2)*(ABS(XDUMMY2))**(1D0/3D0) XNEW=XDUMMY1+XDUMMY2 NUMBER5 = OLISGAUSS(XNEW, DMEAN, DSIGMA) 999 RETURN END DOUBLE PRECISION FUNCTION DLHOOD(X) IMPLICIT NONE DOUBLE PRECISION X #include "make_param.inc" DOUBLE PRECISION OLISGAUSS EXTERNAL OLISGAUSS DLHOOD = DDUMMY*OLISGAUSS(X, DMEAN, DSIGMA)+DCONST+DSLOPE*X 999 RETURN END REAL FUNCTION RDLHOOD(X) IMPLICIT NONE REAL X #include "make_param.inc" DOUBLE PRECISION OLISGAUSS EXTERNAL OLISGAUSS RDLHOOD = REAL(DDUMMY*OLISGAUSS(DBLE(X), DMEAN, DSIGMA)+DCONST & +DSLOPE*X) 999 RETURN END DOUBLE PRECISION FUNCTION OLISGAUSS(X,XMEAN,SIGMA) IMPLICIT NONE DOUBLE PRECISION X,XMEAN,SIGMA DOUBLE PRECISION ARG OLISGAUSS = 0D0 IF (SIGMA.LE.0D0) RETURN ARG = .5D0*((X - XMEAN)/SIGMA)**2 OLISGAUSS = DEXP(-ARG) 999 RETURN END C====================================================================== C Smoothing subroutines C====================================================================== SUBROUTINE SMOOTHING IMPLICIT NONE #include "make_param.inc" #include "abcfit_bmatrix.inc" INTEGER I,J,K,N,PIDMAX,LSTR,LSTR1 REAL RTEMP CHARACTER*10 STR,STR1 CHARACTER*100 HISNAM PIDMAX=0 DO I=1,MAXNUMOFPAR IF (PID(I).GT.PIDMAX) PIDMAX=PID(I) ENDDO DO I=1,PIDMAX IF (PARAMFLAG(I).GE.0) THEN DO J=1,NPAR(PARAMFLAG(I)) RTEMP=REAL(I) CALL CVR2C(RTEMP,0,LSTR,STR) RTEMP=REAL(J) CALL CVR2C(RTEMP,0,LSTR1,STR1) IF (NBIN1.GT.0.AND.NBIN2.GT.0) THEN CALL HBOOK2(NHIST+I*20+J*4,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth vs E mean',NBIN1 & ,B1MIN(I),B1MAX(I),NBIN2,B2MIN(I),B2MAX(I),0.0) CALL HBOOK2(NHIST+I*20+J*4+1,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth vs E mean error' & ,NBIN1,B1MIN(I),B1MAX(I),NBIN2,B2MIN(I),B2MAX(I),0.0) CALL HBOOK2(NHIST+I*20+J*4+2,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth vs E sigma',NBIN1 & ,B1MIN(I),B1MAX(I),NBIN2,B2MIN(I),B2MAX(I),0.0) CALL HBOOK2(NHIST+I*20+J*4+3,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth vs E sigma error' & ,NBIN1,B1MIN(I),B1MAX(I),NBIN2,B2MIN(I),B2MAX(I),0.0) ENDIF IF (NBIN1.GT.0.AND.NBIN2.LE.0) THEN CALL HBOOK1(NHIST+20*I+J*4,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Energy mean',NBIN1 & ,B1MIN(I),B1MAX(I),0.0) CALL HBOOK1(NHIST+20*I+J*4+1,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Energy mean error',NBIN1 & ,B1MIN(I),B1MAX(I),0.0) CALL HBOOK1(NHIST+20*I+J*4+2,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Energy sigma',NBIN1 & ,B1MIN(I),B1MAX(I),0.0) CALL HBOOK1(NHIST+20*I+J*4+3,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Energy sigma error',NBIN1 & ,B1MIN(I),B1MAX(I),0.0) ENDIF IF (NBIN2.GT.0.AND.NBIN1.LE.0) THEN CALL HBOOK1(NHIST+20*I+J*4,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth mean',NBIN2,B2MIN(I) & ,B2MAX(I),0.0) CALL HBOOK1(NHIST+20*I+J*4+1,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth mean error',NBIN2 & ,B2MIN(I),B2MAX(I),0.0) CALL HBOOK1(NHIST+20*I+J*4+2,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth sigma',NBIN2,B2MIN(I) & ,B2MAX(I),0.0) CALL HBOOK1(NHIST+20*I+J*4+3,'PID='//STR(1:LSTR)// & '. PARNUM='//STR1(1:LSTR1)//'. Cth sigma error',NBIN2 & ,B2MIN(I),B2MAX(I),0.0) ENDIF DO K=0,3,2 RTEMP=REAL(NHIST+I*20+J*4+K) CALL CVR2C(RTEMP,0,LSTR,STR) CALL HCOPY (NHIST+I*20+J*4+K,NHIST+I*20+J*4+K+2 & *MAXNUMOFHIS,'Smoothed version of '//STR(1:LSTR)) ENDDO N=0 DO K=2,I IF(PARAMFLAG(K-1).GE.0) N=N+NPAR(PARAMFLAG(K-1)) ENDDO K=1+((J-1)+N)*MAX(NBIN1,1)*MAX(NBIN2,1) CALL HPAK(NHIST+I*20+J*4,FITMEAN(K,1)) CALL HPAK(NHIST+I*20+J*4+1,FITMEAN(K,2)) CALL HPAK(NHIST+I*20+J*4+2,FITSIGMA(K,1)) CALL HPAK(NHIST+I*20+J*4+3,FITSIGMA(K,2)) N=0 DO K=2,I IF(PARAMFLAG(K-1).GE.0) N=N+2*NPAR(PARAMFLAG(K-1)) ENDDO C SMOOTH Y0 CALL FITSMOOTH(SMOOTHFUNC(J,I),N+J,NHIST & +I*20+J*4) C SMOOTH DY0 CALL FITSMOOTH(SMOOTHFUNCERR(J,I),N+NPAR(PARAMFLAG(I))+J & ,NHIST+I*20+J*4+2) ENDDO ENDIF ENDDO 999 RETURN END SUBROUTINE FITSMOOTH(IFLAG,POINTER,IDHIST) IMPLICIT NONE INTEGER IFLAG,POINTER,IDHIST #include "make_param.inc" INTEGER NWT,LOC,I,J,IERFLAG,FITFCN,NPARI,NPARX,ISTAT REAL ASYMFUNCTIONET,SYMFUNCTIONET,XMA,YMA,HIJ,XTH,XE,ROLI_A,ROLI_B DOUBLE PRECISION BNDL(10),BNDU(10),STVAL(10),STEP(10),DARGLIS(10) & ,FMIN,FEDM,ERRDEF,DPAR,DPAR_ERR,BND1,BND2 CHARACTER*10 CHNAM(25) CHARACTER*80 TITL EXTERNAL ASYMFUNCTIONET,SYMFUNCTIONET,HIJ,FITFCN DATA CHNAM /'par(1)','par(2)', 'par(3)', 'par(4)', + 'par(5)','par(6)', 'par(3)', 'par(8)','par(9)', + 'par(10)','par(11)', 'par(12)', 'par(13)','par(14)', + 'par(15)','par(16)', 'par(13)', 'par(18)','par(19)', + 'par(20)','par(21)', 'par(22)', 'par(23)','par(24)', + 'par(25)'/ C C INITIALIZE VALUES C PRINT*,'==================================================' PRINT*,'==================================================' PRINT*,'= FITTING HISTOGRAM ',IDHIST,' =' PRINT*,'==================================================' PRINT*,'==================================================' C C GET THE NUMBER OF ENERGY AND THETA BINS (NENEERGY=;NTHETA=) C CALL HGIVE(IDHIST,TITL,OLI_NX,ROLI_A,XMA,OLI_NY,ROLI_B,YMA,NWT,LOC & ) OLI_A=DBLE(ROLI_A) OLI_B=DBLE(ROLI_B) OLI_WX = (XMA-OLI_A)/REAL(OLI_NX) OLI_WY = (YMA-OLI_B)/REAL(MAX(1,OLI_NY)) HISTOFIT=IDHIST PARPOINTER=POINTER C C LOOP OVER ALL ENERGY SLICES AND AND FIT THETA-DEPENDANCE C SMFUNCFLAG=IFLAG ITHETA=0 IF (NBIN2.LE.0) ITHETA=-1 ! TAG 1-D HISTO FOR ENERGY IF (NBIN1.GT.0) THEN DO IENERGY = 1,OLI_NX C SET THE MINUIT PARAMETER VALUES, LIKE LOWER AND UPPER BOUNDS (NONE) ETC DO I = 1,3 BNDL(I) = 0.0D0 BNDU(I) = 0.D0 STVAL(I) = 0.0D0 IF (I.EQ.1) THEN STEP(I) = 0.0D0 DO J=1,OLI_NY IF (HIJ(IDHIST,IENERGY,J).NE.0.0) STEP(I)=STEP(I)+1.0D0 STVAL(I)=STVAL(I)+DBLE(HIJ(IDHIST,IENERGY,J)) ENDDO IF (STEP(I).GT.0.0D0) STVAL(I)=STVAL(I)/STEP(I) ENDIF STEP(I) = 0.001D0 ENDDO CALL VZERO(SMOOTHPAR(6*(IENERGY-1)+1,1,POINTER),6) IF (STVAL(1).NE.0.0D0) THEN CALL MNINIT(5,6,7) DO I = 1,3 CALL MNPARM (I,CHNAM(I),STVAL(I),STEP(I), + BNDL(I),BNDU(I),IERFLAG) IF (IERFLAG.NE.0) PRINT* & ,'SMOOTH:UNABLE TO DEFINE PARAMETER NUMBER',I,IERFLAG ENDDO DARGLIS(1) = 1.0D0 CALL MNEXCM (FITFCN,'SET PRINTOUT ', + DARGLIS,1, IERFLAG) DARGLIS(1) = 2.0D0 CALL MNEXCM (FITFCN,'SET STRATEGY ', + DARGLIS,1,IERFLAG) DARGLIS(1) = 1.0D0 CALL MNEXCM (FITFCN,'CALL FCN ', + DARGLIS,1,IERFLAG) DARGLIS(1) = 0.0D0 CALL MNEXCM (FITFCN,'MIGRAD ', + DARGLIS,0,IERFLAG) CALL MNSTAT (FMIN,FEDM,ERRDEF,NPARI,NPARX,ISTAT) IF(IERFLAG.EQ.0.AND.ISTAT.EQ.3) THEN C C COPY THE FITTED VALUES OF IN SOME WEIRD COMMON C DO I = 1,3 CALL MNPOUT(I,CHNAM(I),DPAR,DPAR_ERR,BND1,BND2,ISTAT) SMOOTHPAR(6*(IENERGY-1)+(I-1)*2+1,1,POINTER)=REAL(DPAR) SMOOTHPAR(6*(IENERGY-1)+(I-1)*2+2,1,POINTER) & =REAL(DPAR_ERR) ENDDO ENDIF CALL MNEXCM (FITFCN,'CLEAR ', + 0,0,IERFLAG) ENDIF ENDDO ENDIF C C LOOP OVER ALL THETA SLICES AND AND FIT ENERGY-DEPENDENCE C IENERGY=0 IF (NBIN1.LE.0) IENERGY=-1 ! TAG 1-D HISTO FOR CTH IF (NBIN2.GT.0) THEN DO ITHETA=1,3 CALL VZERO(SMOOTHPAR(6*(ITHETA-1)+1,2,POINTER),6) CALL MNINIT(5,6,7) c set the MINUIT parameter values, like lower and upper bounds (none) etc DO I = 1,6 BNDL(I) = 0.D0 BNDU(I) = 0.D0 STEP(I) = 0.001 STVAL(I) = 0.D0 ENDDO c c copy vector and errors c IF (SMOOTHPAR((ITHETA-1)*2+2,1,POINTER).GT.0.0) THEN DO I = 1,6 CALL MNPARM (I,CHNAM(I),STVAL(I),STEP(I),BNDL(I),BNDU(I) & ,IERFLAG) IF (IERFLAG.NE.0) THEN PRINT*,'UNABLE TO DEFINE PARAMETER NUMBER',I,IERFLAG ENDIF ENDDO DARGLIS(1) = 1.0D0 CALL MNEXCM (FITFCN,'SET PRINTOUT ', + DARGLIS,1, IERFLAG) DARGLIS(1) = 2.0D0 CALL MNEXCM (FITFCN,'SET STRATEGY ', + DARGLIS,1,IERFLAG) DARGLIS(1) = 1.0D0 CALL MNEXCM (FITFCN,'CALL FCN ', + DARGLIS,1,IERFLAG) DARGLIS(1) = 0.0D0 CALL MNEXCM (FITFCN,'MIGRAD ', + DARGLIS,0,IERFLAG) CALL MNSTAT (FMIN,FEDM,ERRDEF,NPARI,NPARX,ISTAT) IF(IERFLAG.EQ.0.AND.ISTAT.GE.1) THEN C C COPY THE FITTED VALUES OF IN SOME WEIRD COMMON C DO I = 1,6 CALL MNPOUT(I,CHNAM(I),DPAR,DPAR_ERR,BND1,BND2,ISTAT) SMOOTHPAR(6*(ITHETA-1)+I,2,POINTER)=REAL(DPAR) ENDDO ENDIF CALL MNEXCM (FITFCN,'CLEAR ', + 0,0,IERFLAG) ENDIF ENDDO ENDIF DO I = 1,OLI_NX DO J = 1,MAX(1,OLI_NY) XE = OLI_A+OLI_WX*REAL(I-0.5) XTH = OLI_B+OLI_WY*REAL(J-0.5) IF (IENERGY.LT.0) THEN ! CHECK FOR 1-D HISTO IN CTH XTH=XE XE=0.0 ENDIF C to plot only cthe functions in each energy bin: C IF (SMFUNCFLAG.EQ.1) THEN C STVAL(1)=SMOOTHPAR(6*(I-1)+1,1,POINTER) C STVAL(2)=SMOOTHPAR(6*(I-1)+3,1,POINTER) C STVAL(3)=SMOOTHPAR(6*(I-1)+5,1,POINTER) C CALL HFILL(IDHIST+2*MAXNUMOFHIS,XE,XTH,REAL(STVAL(1)+STVAL(2 C & )*XTH**3+STVAL(3)*XTH**5)) C ELSE C STVAL(1)=SMOOTHPAR(6*(I-1)+1,1,POINTER) C STVAL(2)=SMOOTHPAR(6*(I-1)+3,1,POINTER) C STVAL(3)=SMOOTHPAR(6*(I-1)+5,1,POINTER) C CALL HFILL(IDHIST+2*MAXNUMOFHIS,XE,XTH,REAL(STVAL(1)+STVAL(2 C & )*XTH**6+STVAL(3)*XTH**8)) C ENDIF IF (SMFUNCFLAG.EQ.1) CALL HFILL(IDHIST+2*MAXNUMOFHIS,XE,XTH & ,ASYMFUNCTIONET(XTH,XE)) IF (SMFUNCFLAG.EQ.2) CALL HFILL(IDHIST+2*MAXNUMOFHIS,XE,XTH & ,SYMFUNCTIONET(XTH,XE)) ENDDO ENDDO 999 RETURN END SUBROUTINE FITFCN(NPAR,GRAD,FVAL,PAR,IFLAG) IMPLICIT NONE INTEGER NPAR,IFLAG DOUBLE PRECISION GRAD(*),FVAL,PAR(*) #include "make_param.inc" INTEGER I REAL XE,XTH,VAL,ERR,HIJ EXTERNAL HIJ C C TAG ON DIMENSION C FVAL=0.0D0 IF (IENERGY.GT.0.AND.ITHETA.EQ.0) THEN DO I = 1,OLI_NY XTH = OLI_B+OLI_WY*REAL(I-0.5) VAL = HIJ(HISTOFIT,IENERGY,I) ERR = HIJ(HISTOFIT+1,IENERGY,I) IF (ERR.GT.0.000001) THEN IF (SMFUNCFLAG.EQ.1) FVAL=FVAL+(VAL-PAR(1)-PAR(2)*XTH**3 & -PAR(3)*XTH**5)**2/ERR**2 IF (SMFUNCFLAG.EQ.2) FVAL=FVAL+(VAL-PAR(1)-PAR(2)*XTH**6 & -PAR(3)*XTH**8)**2/ERR**2 ENDIF ENDDO ELSEIF (ITHETA.GT.0.AND.IENERGY.EQ.0) THEN DO I = 1,OLI_NX XE = OLI_A+OLI_WX*REAL(I-0.5) IF (SMOOTHPAR(6*(I-1)+(ITHETA-1)*2+2,1,PARPOINTER).GT.0 & .000001) & FVAL=FVAL+( & DBLE(SMOOTHPAR(6*(I-1)+(ITHETA-1)*2+1,1,PARPOINTER))- & (PAR(1)+PAR(2)*XE+PAR(3)*XE**2+PAR(4)*XE**3+PAR(5)*XE**4 & +PAR(6)*XE**5) & )**2/DBLE(SMOOTHPAR(6*(I-1)+(ITHETA-1)*2+2,1,PARPOINTER))**2 ENDDO ENDIF 999 RETURN END REAL FUNCTION ASYMFUNCTIONET(TH,E) C SMOOTH FUNCTION FOR THETA AND ENERGY SMOOTHING IMPLICIT NONE REAL TH,E #include "make_param.inc" INTEGER I REAL E_PAR(6,3),TH_PAR(3) DO I = 1,6 E_PAR(I,1) = SMOOTHPAR(I,2,PARPOINTER) E_PAR(I,2) = SMOOTHPAR(6+I,2,PARPOINTER) E_PAR(I,3) = SMOOTHPAR(12+I,2,PARPOINTER) ENDDO DO I=1,3 TH_PAR(I) = E_PAR(1,I)+E_PAR(2,I)*E+E_PAR(3,I)*E**2+E_PAR(4,I)*E & **3+E_PAR(5,I)*E**4+E_PAR(6,I)*E**5 ENDDO ASYMFUNCTIONET = TH_PAR(1)+TH_PAR(2)*TH**3+TH_PAR(3)*TH**5 999 RETURN END REAL FUNCTION SYMFUNCTIONET(TH,E) C SMOOTH FUNCTION FOR THETA AND ENERGY SMOOTHING IMPLICIT NONE REAL TH,E #include "make_param.inc" INTEGER I REAL E_PAR(6,3),TH_PAR(3) DO I = 1,6 E_PAR(I,1) = SMOOTHPAR(I,2,PARPOINTER) E_PAR(I,2) = SMOOTHPAR(6+I,2,PARPOINTER) E_PAR(I,3) = SMOOTHPAR(12+I,2,PARPOINTER) ENDDO DO I=1,3 TH_PAR(I) = E_PAR(1,I)+E_PAR(2,I)*E+E_PAR(3,I)*E**2+E_PAR(4,I)*E & **3+E_PAR(5,I)*E**4+E_PAR(6,I)*E**5 ENDDO SYMFUNCTIONET = TH_PAR(1)+TH_PAR(2)*TH**6+TH_PAR(3)*TH**8 999 RETURN END C====================================================================== C Output parameter routine C====================================================================== SUBROUTINE OUTPARAM IMPLICIT NONE #include "make_param.inc" #include "abcfit_bmatrix.inc" INTEGER I,J,K,L,N,PIDMAX,NPARAMS,FLEN CHARACTER*10 NUMBERS CHARACTER*100 FILENAME DATA NUMBERS /'0123456789'/ C C CONSTRUCT FILENAME C FILENAME='aibi_evol_' FLEN=10 PIDMAX=0 DO I=1,MAXNUMOFPAR IF (PID(I).GT.PIDMAX) PIDMAX=PID(I) IF (PARAMFLAG(PID(I)).GE.0.AND.PID(I).GT.0) THEN FILENAME=FILENAME(1:FLEN)//NUMBERS(PARAMFLAG(PID(I)) & +1:PARAMFLAG(PID(I))+1) FLEN=FLEN+1 ENDIF ENDDO FILENAME=FILENAME(1:FLEN)//'_'//NUMBERS(ABS(BINLEVEL-3 & ):ABS(BINLEVEL-3))//'0x.dat' FLEN=FLEN+8 OPEN(UNIT_OUT,FORM='FORMATTED',STATUS='NEW',FILE & =FILENAME(1:FLEN)) C C FIRST WRITE 20 COMMENT LINES (TO BE EDITED BY THE USER AFTERWARDS) C DO I=1,20 IF (I.EQ.10) THEN WRITE(UNIT_OUT,'(a)',ERR=900) & '* Add your comments in the first 20 lines' ELSE WRITE(UNIT_OUT,'(a)',ERR=900) '*' ENDIF ENDDO WRITE(UNIT_OUT,'(a)',ERR=900) '*Number of particles required' & //' for this parameterisation file' WRITE(UNIT_OUT,'(i1)',ERR=900) NJET WRITE(UNIT_OUT,'(a)',ERR=900) '*Specify different'// & ' parameterisations given in this parameterisation file' WRITE(UNIT_OUT,'(6(i2,a),i2)',ERR=900) PARAMFLAG(1),',' & ,PARAMFLAG(2),',',PARAMFLAG(3),',',PARAMFLAG(4),',',PARAMFLAG(5) & ,',',PARAMFLAG(6),',',PARAMFLAG(7) WRITE(UNIT_OUT,'(a)',ERR=900) '*Specify the version of the ' & //'read-in structure of the rest of the file' C C IF ANY-ONE _EVER_ CHANGES THE STRUCTURE OF THE FILE THEN SET VERSION!!! C AND REMEMBER TO ADD READ-IN CODE TO AIBI_EVOL.F CODE.... C WRITE(UNIT_OUT,'(i1,a,i1)',ERR=900) 1,',',1 C C NOW WRITE FUNCTION DEFINITIONS AND PARAMETERS FOR Y0 AND DY0 C DO I=1,PIDMAX N=0 DO K=2,I IF (PARAMFLAG(K-1).GE.0) N=N+2*NPAR(PARAMFLAG(K-1)) ENDDO IF (PARAMFLAG(I).GE.0) THEN WRITE(UNIT_OUT,'(a)',ERR=900) & '*Parameterisation and its energy'//' and cos_theta ranges' WRITE(UNIT_OUT,'(i2,a,3(f6.2,a),f6.2)',ERR=900) PARAMFLAG(I), & ',',B1MIN(I),',',B1MAX(I),',',B2MIN(I),',',B2MAX(I) DO J=1,NPAR(PARAMFLAG(I)) WRITE(UNIT_OUT,'(a)',ERR=900) '*FUNCTYP,PARNUM' & //',START_PARTICLE,END_PARTICLE,#PARAMETERS' WRITE(UNIT_OUT,'(4(i2,a),i3)',ERR=900) SMOOTHFUNC(J,I),',',J & ,',',USEPAR(1,I),',',USEPAR(2,I),',',18 WRITE(UNIT_OUT,'(a)',ERR=900) '*PARAMETERS FOR THE FUNCTION' DO K=1,6 WRITE(UNIT_OUT,'(E12.6,A,E12.6,A,E12.6)',ERR=900) & SMOOTHPAR((K-1)*3+1,2,N+J),',',SMOOTHPAR((K-1)*3+2,2,N+J & ),',',SMOOTHPAR((K-1)*3+3,2,N+J) ENDDO WRITE(UNIT_OUT,'(a)',ERR=900) & '*FUNCTYP,PARNUM(ERR)' & //',START_PARTICLE,END_PARTICLE,#PARAMETERS' WRITE(UNIT_OUT,'(4(i2,a),i3)',ERR=900) SMOOTHFUNCERR(J,I) & ,',',J+NPAR(PARAMFLAG(I)),',',USEPAR(1,I),',',USEPAR(2,I), & ',',18 WRITE(UNIT_OUT,'(a)',ERR=900) '*PARAMETERS FOR THE FUNCTION' & //'(ERR)' DO K=1,6 WRITE(UNIT_OUT,'(E12.6,A,E12.6,A,E12.6)',ERR=900) & SMOOTHPAR((K-1)*3+1,2,N+NPAR(PARAMFLAG(I))+J),',' & ,SMOOTHPAR((K-1)*3+2,2,N+NPAR(PARAMFLAG(I))+J),',' & ,SMOOTHPAR((K-1)*3+3,2,N+NPAR(PARAMFLAG(I))+J) ENDDO ENDDO ENDIF ENDDO C C NOW WRITE DEFINITIONS AND PARAMETERS FOR CORRELATION MATRICES C WRITE(UNIT_OUT,'(a)',ERR=900) '*Particle_number_for_ordering,' & //'Nbin_energy,Nbin_theta for correlation matrices' C C ONLY ONE GLOBAL CORRELATION MATRIX... C WRITE(UNIT_OUT,'(I2,A,I2,A,I2)',ERR=900) 1,',',1,',',1 WRITE(UNIT_OUT,'(a)',ERR=900) '*Energy and cos_theta ranges' C C FOR GLOBAL CORRELATION WE USE ALL EVENTS -- CHOOSE WIDE RANGES C WRITE(UNIT_OUT,'(F8.2,A,F8.2,A,F8.2,A,F8.2)',ERR=900) 0.0,',',300 & .0,',',-1.0,',',1.0 WRITE(UNIT_OUT,'(a)',ERR=900) '*Correlation matrix...' NPARAMS=1 DO WHILE (COR(NPARAMS,NPARAMS).GT.0.0) NPARAMS=NPARAMS+1 ENDDO NPARAMS=NPARAMS-1 DO I=1,NPARAMS DO J=1,NPARAMS,3 IF (J+2.LE.NPARAMS) WRITE(UNIT_OUT,'(E12.6,A,E12.6,A,E12.6)' & ,ERR=900) COR(J,I),',',COR(J+1,I),',',COR(J+2,I) IF (J+2.EQ.NPARAMS+1) WRITE(UNIT_OUT,'(E12.6,A,E12.6)' & ,ERR=900) COR(J,I),',',COR(J+1,I) IF (J+2.EQ.NPARAMS+2) WRITE(UNIT_OUT,'(E12.6)' & ,ERR=900) COR(J,I) ENDDO ENDDO C C DONE C PRINT*,'File: ',FILENAME(1:FLEN),' has been successfully written' PRINT*,'Please check and edit/add information to the ' & //'first 20 comment lines!!!!' PRINT*,'AND replace x in the filename with your prefered ITEVOLN' & //' value for this event type!!!!' PRINT*,'Have fun....' CLOSE(UNIT_OUT) RETURN 900 CONTINUE CLOSE(UNIT_OUT) PRINT*,'ERROR writing File: ',FILENAME(1:FLEN) & ,' -- operation aborted!!!!' 999 RETURN END C====================================================================== C True<->Reco Matching subroutines C====================================================================== SUBROUTINE MATCH(MATCH_FLAG,PIDFLAG) IMPLICIT NONE C INPUT INTEGER MATCH_FLAG,PIDFLAG #include "make_param.inc" C INTERNAL INTEGER I,J,K,NPAR INTEGER IPERM(MAXNUMOFPAR),NPERM,IPID(MAXNUMOFPAR) & ,IORDER(MAXNUMOFPAR) REAL CHI,CHIMIN,PTEMP(4,MAXNUMOFPAR),VAL NPAR=0 DO I=1,MAXNUMOFPAR IF (PID(I).EQ.PIDFLAG) THEN NPAR=NPAR+1 IPID(NPAR)=I ENDIF ENDDO C NPAR.EQ.1 IS SIMPLE IF (NPAR.LE.1) RETURN C CALCULATE THE NUMBER OF PERMUTATIONS = NPAR!=NPERM NPERM=1 DO I=1,NPAR NPERM=NPERM*I ENDDO CHIMIN=10000000.0 IPERM(1)=0 ! set array values to 1,2,3,....NPAR IF (MATCH_FLAG.EQ.3) THEN NPERM=1 DO K=1,4 PTEMP(K,1)=0.0D0 PTEMP(K,2)=0.0D0 ENDDO ENDIF DO I=1,NPERM CALL PERMU(IPERM,NPAR) CHI=0.0 DO J=1,NPAR K=IPID(IPERM(J)) IF (MATCH_FLAG.EQ.1) CALL ANGLE_CAL(IPID(J),K,VAL) IF (MATCH_FLAG.EQ.2) CALL MASS_CAL(IPID(J),K,VAL) IF (MATCH_FLAG.EQ.3) THEN DO K=1,4 PTEMP(K,1)=PTEMP(K,1)+PRECO(K,IPID(J)) PTEMP(K,2)=PTEMP(K,2)+PTRUE(K,IPID(J)) ENDDO ENDIF CHI=CHI+VAL ENDDO IF (CHI.LE.CHIMIN) THEN CHIMIN=CHI DO K=1,NPAR IORDER(K)=IPERM(K) ENDDO ENDIF ENDDO IF (MATCH_FLAG.EQ.3) THEN DO I=2,NPAR PID(IPID(I))=0 ENDDO DO I=1,4 PRECO(I,IPID(1))=PTEMP(I,1) PTRUE(I,IPID(1))=PTEMP(I,2) ENDDO ELSE DO I=1,NPAR DO J=1,4 PTEMP(J,IPID(I))=PRECO(J,IPID(IORDER(I))) ENDDO ENDDO DO I=1,NPAR DO J=1,4 PRECO(J,IPID(I))=PTEMP(J,IPID(I)) ENDDO ENDDO ENDIF 999 RETURN END SUBROUTINE MASS_CAL(I1,I2,VAL) IMPLICIT NONE INTEGER I1,I2,I REAL VAL,PSUM(4),P #include "ntparam.inc" P=0.0 DO I=1,4 PSUM(I)=PTRUE(I,I1)+PRECO(I,I2) IF (I.LE.3) P=P+PSUM(I)**2 ENDDO VAL=SQRT(MAX(0.0,PSUM(4)**2-P)) RETURN END SUBROUTINE ANGLE_CAL(I1,I2,VAL) IMPLICIT NONE INTEGER I1,I2,I REAL VAL,P1,P2 #include "ntparam.inc" VAL=0.0 P1=0.0 P2=0.0 DO I=1,3 VAL=VAL+PTRUE(I,I1)*PRECO(I,I2) P1=P1+PTRUE(I,I1)**2 P2=P2+PRECO(I,I2)**2 ENDDO VAL=VAL/SQRT(P1)/SQRT(P2) IF (VAL.GE.0.0) VAL=ACOS(MIN(1.,VAL)) IF (VAL.LT.0.0) VAL=ACOS(MAX(-1.,VAL)) RETURN END C====================================================================== C Ntuple IO subroutines C====================================================================== SUBROUTINE NTOPEN IMPLICIT NONE #include "make_param.inc" INTEGER ISTAT C C OPEN INPUT FILE C CALL HROPEN(UNIT_IN,'INPUT',INPUT_FILE,' ',IREC,ISTAT) IF (ISTAT.NE.0) THEN PRINT*,'**********************************************' PRINT*,'**********************************************' PRINT*,'ERROR BY OPENING INPUT-FILE',INPUT_FILE PRINT*,'**********************************************' PRINT*,'**********************************************' STOP ENDIF PRINT*,'OPENED INPUT-FILE: ',INPUT_FILE CALL HRIN(ID_IN, 999999, 0) CALL HPRNT(ID_IN) CALL HBNAME(ID_IN,' ',0,'$CLEAR') CALL HBNAME(ID_IN,'CFPARM',NJET,'$SET') CALL HNOENT(ID_IN, NEVT) PRINT*,'NUMBER OF EVENTS ARE',NEVT C C OPEN OUTPUT FILE C CALL HROPEN(UNIT_CHECK,'OUTPUT',OUTPUT_FILE_CHECK,'N ',IREC & ,ISTAT) IF (ISTAT.NE.0) THEN PRINT*,'**********************************************' PRINT*,'**********************************************' PRINT*,'ERROR BY OPENING OUTPUT-CHECK-FILE',OUTPUT_FILE_CHECK PRINT*,'**********************************************' PRINT*,'**********************************************' STOP ENDIF PRINT*,'OPENED OUTPUT-CHECK-FILE: ',OUTPUT_FILE_CHECK 999 RETURN END SUBROUTINE NTREAD(I,IERR) IMPLICIT NONE INTEGER I,IERR #include "make_param.inc" CHARACTER*50 OLDDIR IERR=0 CALL HCDIR(OLDDIR,'R') CALL HCDIR('//INPUT',' ') CALL HGNTB(ID_IN,'CFPARM',I,IERR) CALL HCDIR(OLDDIR,' ') 999 RETURN END SUBROUTINE NTCLOSE(IFLAG) IMPLICIT NONE INTEGER IFLAG #include "make_param.inc" INTEGER ISTAT IF (IFLAG.EQ.1) THEN CALL HDELET(ID_IN) CALL HREND('INPUT') CLOSE(UNIT_IN) PRINT*,'INPUT_FILE ',INPUT_FILE,' IS NOW CLOSED' ENDIF IF (IFLAG.EQ.2) THEN ISTAT=0 CALL HCDIR('//OUTPUT',' ') CALL HROUT(0,ISTAT,' ') CALL HREND('OUTPUT') CLOSE(UNIT_CHECK) PRINT*,'OUTPUT-CHECK-FILE ',OUTPUT_FILE_CHECK,' IS NOW CLOSED' ENDIF 999 RETURN END C====================================================================== C Special misc. subroutines C====================================================================== SUBROUTINE CVR2C(NUM,NDEC,L,STR) IMPLICIT NONE INTEGER NDEC,L REAL NUM CHARACTER*(*) STR INTEGER I,J,K CHARACTER*11 NUMBER DATA NUMBER /'12345678900'/ SAVE NUMBER I=INT(NUM*10.0**NDEC) STR='0' L=1 IF (I.EQ.0) RETURN L=0 IF (I.LT.0) THEN I=ABS(I) STR='-' L=1 ENDIF IF (ABS(NUM).LT.1.0) THEN IF (L.LE.0) THEN STR='0' L=1 ELSE STR='-0' L=2 ENDIF ENDIF DO J=1+INT(LOG(REAL(I))/LOG(10.0)),1,-1 K=MOD(I/10**(J-1),10) IF (K.EQ.0) K=10 IF (J.EQ.NDEC.AND.NDEC.GT.0) THEN STR=STR(1:L)//'.' L=L+1 ENDIF IF (L.EQ.0) THEN STR=NUMBER(K:K) ELSE STR=STR(1:L)//NUMBER(K:K) ENDIF L=L+1 ENDDO 999 RETURN END SUBROUTINE GENHISNAM(PID,PARNUM,ELOW,EUP,CTHLOW,CTHUP,L,HISNAM) IMPLICIT NONE INTEGER PID,PARNUM,L REAL ELOW,EUP,CTHLOW,CTHUP CHARACTER*(*) HISNAM INTEGER I REAL RTEMP CHARACTER*10 STR RTEMP=REAL(PID) CALL CVR2C(RTEMP,0,I,STR) HISNAM=' '//STR(1:I)//' ' L=1+I+1 RTEMP=REAL(PARNUM) CALL CVR2C(RTEMP,0,I,STR) HISNAM=HISNAM(1:L)//STR(1:I)//' (' L=L+I+2 CALL CVR2C(ELOW,2,I,STR) HISNAM=HISNAM(1:L)//STR(1:I)//',' L=L+I+1 CALL CVR2C(EUP,2,I,STR) HISNAM=HISNAM(1:L)//STR(1:I)//')(' L=L+I+2 CALL CVR2C(CTHLOW,2,I,STR) HISNAM=HISNAM(1:L)//STR(1:I)//',' L=L+I+1 CALL CVR2C(CTHUP,2,I,STR) HISNAM=HISNAM(1:L)//STR(1:I)//')' L=L+I+1 999 RETURN END SUBROUTINE GETLIM(CLFLAG,VAR1,VAR2,ELOW,EUP,CTHLOW,CTHUP) IMPLICIT NONE INTEGER CLFLAG REAL ELOW,EUP,CTHLOW,CTHUP DOUBLE PRECISION VAR1,VAR2 #include "make_param.inc" REAL RTEMP ELOW=0.0 EUP=0.0 IF (NBIN1.GT.0) THEN RTEMP = (B1MAX(CLFLAG)-B1MIN(CLFLAG))/REAL(NBIN1) ELOW=B1MIN(CLFLAG)+REAL(MAX(0,INT((REAL(VAR1)-B1MIN(CLFLAG)) & /RTEMP)))*RTEMP EUP=B1MIN(CLFLAG)+REAL(MIN(NBIN1,MAX(0,INT((REAL(VAR1) & -B1MIN(CLFLAG))/RTEMP)))+1)*RTEMP ENDIF CTHLOW=0.0 CTHUP=0.0 IF (NBIN2.GT.0) THEN RTEMP = (B2MAX(CLFLAG)-B2MIN(CLFLAG))/REAL(NBIN2) CTHLOW=B2MIN(CLFLAG)+REAL(MAX(0,INT((REAL(VAR2)-B2MIN(CLFLAG)) & /RTEMP)))*RTEMP CTHUP=B2MIN(CLFLAG)+REAL(MIN(NBIN2,MAX(0,INT((REAL(VAR2) & -B2MIN(CLFLAG))/RTEMP))+1))*RTEMP ENDIF 999 RETURN END SUBROUTINE WRITE_INC IMPLICIT NONE #include "make_param.inc" #include "abcfit_bmatrix.inc" INTEGER I,J,K,L,N,PIDMAX,NPARAMS,NCLASS UNIT_OUT=UNIT_OUT+1 PIDMAX=0 NCLASS=0 DO I=1,MAXNUMOFPAR IF (PARAMFLAG(I).GE.0) NCLASS=NCLASS+1 IF (PID(I).GT.PIDMAX) PIDMAX=PID(I) ENDDO OPEN(UNIT_OUT,FORM='FORMATTED',STATUS='NEW',FILE='aibi_evol.inc' & ) C C FIRST WRITE 20 COMMENT LINES (TO BE EDITED BY THE USER AFTERWARDS) C DO I=1,20 IF (I.EQ.10) THEN WRITE(UNIT_OUT,'(A1)',ERR=900) & 'C Add your comments in the first 20 lines' ELSE WRITE(UNIT_OUT,'(A1)',ERR=900) 'C' ENDIF ENDDO C C NOW WRITE FUNCTION DEFINITIONS AND PARAMETERS FOR Y0 AND DY0 C WRITE(UNIT_OUT,'(A)',ERR=900) 'C' WRITE(UNIT_OUT,'(A)',ERR=900) & 'C PAR_SMOTH(NPAR_SMOOTH,NPAR_PARTICLE,NCLASS)' WRITE(UNIT_OUT,'(A)',ERR=900) & 'C WITH NPAR_SMOOTH = NUMBER OF PARAMETERS FOR '// & 'SMOOTHING FUNCTION' WRITE(UNIT_OUT,'(A)',ERR=900) & 'C NPAR_PARTICLE = NUMBER OF PARAMETERS FOR '// & 'EACH PARTICLE ' WRITE(UNIT_OUT,'(A)',ERR=900) & 'C NCLASS = NUMBER OF CLASSES FOR '// & 'PARAMETRISATION ' WRITE(UNIT_OUT,'(A)',ERR=900) 'C' NPARAMS=1 DO WHILE (COR(NPARAMS,NPARAMS).GT.0.0) NPARAMS=NPARAMS+1 ENDDO NPARAMS=NPARAMS-1 I=3 J=18 WRITE(UNIT_OUT, & '(A6,A18,A12,I2,A1,I1,A1,I1,A2,A12,I2,A1,I1,A1,I1,A1)')' ' & ,'DOUBLE PRECISION ','PAR_SMOOTH(',J,',',I,',',NCLASS,'),' & ,'EPAR_SMOOTH(',J,',',I,',',NCLASS,')' WRITE(UNIT_OUT & ,'(A6,A12,I2,A1,I2,A1)')' &',' ,COR(',NPARAMS,',' & ,NPARAMS,')' WRITE(UNIT_OUT & ,'(A6,A5,A6,I1,A1,A1,A8,I1,A1,A1,A6,I1,A1,A1,A8,I1,A1)')' ' & ,'REAL ','E_MIN(',PIDMAX,')',',','CTH_MIN(',PIDMAX,')',',' & ,'E_MAX(',PIDMAX,')',',','CTH_MAX(',PIDMAX,')' WRITE(UNIT_OUT & ,'(A6,A8,A8,I1,A1,I1,A2,A9,I1,A1,I1,A2,A6,I1,A2,A7,A7)' & ) ' ','INTEGER ','FUNCTYP(',I,',',PIDMAX,'),','FUNCTYPE(' & ,I,',',PIDMAX,'),','ITYPP(',PIDMAX,'),','EVETYP,','BINTRUE' WRITE(UNIT_OUT,'(A6,A7,A8,I1,A1)') & ' &','PIDMAX',',USEP(2,',PIDMAX,')' DO I=1,3 WRITE(UNIT_OUT,'(A1)',ERR=900) 'C' ENDDO WRITE(UNIT_OUT,'(A6,A7,I2)')' ','EVETYP=',EVENTFLAG WRITE(UNIT_OUT,'(A6,A7,I2)')' ','PIDMAX=',PIDMAX WRITE(UNIT_OUT,'(A6,A8,I2)')' ','BINTRUE=',BINLEVEL DO I=1,PIDMAX WRITE(UNIT_OUT,'(A20,I2)',ERR=900) 'C-------------CLASS=',I WRITE(UNIT_OUT,'(A6,A6,I2,A2,I2)')' ','ITYPP(',I,')=' & ,PARAMFLAG(I) WRITE(UNIT_OUT,'(A6,A6,I2,A2,F10.4)')' ','E_MIN(',I,')=' & ,B1MIN(I) WRITE(UNIT_OUT,'(A6,A6,I2,A2,F10.4)')' ','E_MAX(',I,')=' & ,B1MAX(I) WRITE(UNIT_OUT,'(A6,A8,I2,A2,F10.4)')' ','CTH_MIN(',I,')=' & ,B2MIN(I) WRITE(UNIT_OUT,'(A6,A8,I2,A2,F10.4)')' ','CTH_MAX(',I,')=' & ,B2MAX(I) WRITE(UNIT_OUT,'(A6,A7,I2,A2,I2)')' ','USEP(1,',I,')=' & ,USEPAR(1,I) WRITE(UNIT_OUT,'(A6,A7,I2,A2,I2)')' ','USEP(2,',I,')=' & ,USEPAR(2,I) IF (PARAMFLAG(I).GE.0) THEN DO J=1,NPAR(PARAMFLAG(I)) WRITE(UNIT_OUT,'(A6,A8,I1,A1,I1,A3,I2)')' ','FUNCTYP(',J & ,',',I,')=',SMOOTHFUNC(J,I) WRITE(UNIT_OUT,'(A6,A9,I1,A1,I1,A3,I2)')' ','FUNCTYPE(',J & ,',',I,')=',SMOOTHFUNCERR(J,I) ENDDO ENDIF ENDDO DO I=1,3 WRITE(UNIT_OUT,'(A1)',ERR=900) 'C' ENDDO L=0 WRITE(UNIT_OUT,'(A6,A20)') & ' ','DATA PAR_SMOOTH/' DO I=1,PIDMAX N=0 DO K=1,I IF (PARAMFLAG(K).GE.0) N=N+2*NPAR(PARAMFLAG(K)) ENDDO IF (PARAMFLAG(I).GE.0) THEN DO J=1,NPAR(PARAMFLAG(I)) C---------------------write out the parameters for the y0's DO K=1,6 L=L+1 IF (L.LE.NCLASS*NPAR(PARAMFLAG(I))*6-1) THEN WRITE(UNIT_OUT,'(A6,E12.6,A,E12.6,A,E12.6,A)',ERR & =900)' &',SMOOTHPAR((K-1)*3+1,2,N+J),',' & ,SMOOTHPAR((K-1)*3+2,2,N+J),',',SMOOTHPAR((K-1)*3 & +3,2,N+J),',' ELSE WRITE(UNIT_OUT,'(A6,E12.6,A,E12.6,A,E12.6)',ERR=900) & ' &',SMOOTHPAR((K-1)*3+1,2,N+J),',',SMOOTHPAR((K & -1)*3+2,2,N+J),',',SMOOTHPAR((K-1)*3+3,2,N+J) ENDIF ENDDO ENDDO ENDIF ENDDO WRITE(UNIT_OUT,'(A6,A1)') & ' &','/' C---------------------write out the parameters for the dy0's L=0 WRITE(UNIT_OUT,'(A6,A20)') & ' ','DATA EPAR_SMOOTH/' DO I=1,PIDMAX N=0 DO K=1,I IF (PARAMFLAG(K).GE.0) N=N+2*NPAR(PARAMFLAG(K)) ENDDO IF (PARAMFLAG(I).GE.0) THEN DO J=1,NPAR(PARAMFLAG(I)) DO K=1,6 L=L+1 IF (L.LE.NCLASS*NPAR(PARAMFLAG(I))*6-1) THEN WRITE(UNIT_OUT,'(A6,E12.6,A,E12.6,A,E12.6,A)',ERR & =900)' &',SMOOTHPAR((K-1)*3+1,2,N & +NPAR(PARAMFLAG(I))+J),',',SMOOTHPAR((K-1)*3+2,2,N & +NPAR(PARAMFLAG(I))+J),',',SMOOTHPAR((K-1)*3+3,2,N & +NPAR(PARAMFLAG(I))+J),',' ELSE WRITE(UNIT_OUT,'(A6,E12.6,A,E12.6,A,E12.6)',ERR=900) & ' &',SMOOTHPAR((K-1)*3+1,2,N+NPAR(PARAMFLAG(I))+J & ),',',SMOOTHPAR((K-1)*3+2,2,N+NPAR(PARAMFLAG(I))+J & ),',',SMOOTHPAR((K-1)*3+3,2,N+NPAR(PARAMFLAG(I))+J & ) ENDIF ENDDO ENDDO ENDIF ENDDO WRITE(UNIT_OUT,'(A6,A1)') & ' &','/' C C ONLY ONE GLOBAL CORRELATION MATRIX... C WRITE(UNIT_OUT,'(a)',ERR=900) '*Correlation matrix...' WRITE(UNIT_OUT,'(A6,A20)') & ' ','DATA COR/' L=0 DO I=1,NPARAMS DO J=1,NPARAMS L=L+1 IF (L.LE.NPARAMS*NPARAMS-1) THEN WRITE(UNIT_OUT & ,'(A6,E12.6,A1)',ERR=900) ' &',COR(J,I),',' ELSE WRITE(UNIT_OUT & ,'(A6,E12.6)',ERR=900) ' &',COR(J,I) ENDIF ENDDO ENDDO WRITE(UNIT_OUT,'(A6,A1)') & ' &','/' C C DONE C PRINT*,'File: aibi_evol.inc has been successfully written' PRINT*,'Please check and edit/add information to the ' & //'first 20 comment lines!!!!' PRINT*,'Have fun....' RETURN 900 CONTINUE CLOSE(UNIT_OUT) PRINT*,'ERROR writing File aibi_evol.inc -- operation aborted!!!!' 999 RETURN END