*CMZ :  4.20/13 19/10/93  15.35.31  by  Roger Barlow, Christine Beeston
*CMZ :  4.20/09 24/09/93  16.27.34  by  Roger Barlow, Christine Beeston
*-- Author :    Roger Barlow, Christine Beeston   24/09/93
C glitch caused by comparing reals corrected (I hope)  RJB 29/11/03
C      thanks to Peter Litchfield

      SUBROUTINE HADJUST(TI,I,P,KZERO,AKI)

C Subroutine which given a set of Pj (answer), solves for the ti
C (ti = 1 - di/fi, where di is the number of data events in bin i,
C and fi is the number of mc events predicted in bin i, given the Pj
C provided.  fi = Nd * sum over j of (Pj aji/Nj), where Nd = total data
C events, Nj is total number of mc events from source j, aji is number
C of mc events from source j in bin i.

*KEEP,HCMCPM.

C some maximum values - max number of bins and dimensions
       INTEGER NSRCMX,NSRCMN,NOPTS
       PARAMETER(NSRCMX=20,NSRCMN=2,NOPTS=7)

C Numbers of monte carlo and data events, links to data, MC and weight
C histograms, Number of mc sources, options, total number of bins.
C normalisation constants for weight histograms

       DOUBLE PRECISION BJ(NSRCMX)
       INTEGER NDATEV,NMCEV(NSRCMX),
     + NMCSRC,IOPT(NOPTS),NTOT
       COMMON/HINPUTS/NDATEV,NMCEV,NMCSRC,IOPT,NTOT,BJ
C Histogram IDs
       INTEGER IDD,IDM(NSRCMX),IDW(NSRCMX)
       COMMON/HMCIDS/IDD,IDM,IDW

*KEND.

       DOUBLE PRECISION TI,P(NSRCMX),WPMAX,TMIN,AKI,
     + STEP,FUNC,DERIV,D,DELTA,WJI,WKI

       INTEGER I,KZERO,IMAX,J,NPMAX,APMAXS


C set up largest step value
        STEP=0.2D0

C First case - di=0 -> ti=1
       IF(NINT(HI(IDD,I)).EQ.0)THEN
          TI=1.0D0
          KZERO=0
          RETURN
       ENDIF

C Find the largest pj (for case when one or more of the aji are zero)
       WPMAX=HI(IDW(1),I)*P(1)
       IMAX=1
       DO 10 J=2,NMCSRC
          WJI=HI(IDW(J),I)
          IF(WJI*P(J).GT.WPMAX)THEN
             WPMAX=WJI*P(J)
             IMAX=J
          ENDIF
10     CONTINUE
C count the number of sources having p=pmax, and the sum of their aji
       NPMAX=0
       APMAXS=0
       DO 20J=1,NMCSRC
          WJI=HI(IDW(J),I)
C glitch correction
          IF(WJI*P(J).GE.WPMAX.or.J.EQ.IMAX) THEN
C          IF(WJI*P(J).EQ.WPMAX)THEN
             NPMAX=NPMAX+1
             APMAXS=APMAXS+HI(IDM(J),I)
          ENDIF
C check that none of the p(j) are zero as this causes crashes
          IF(WJI*P(J).EQ.0.0D0)THEN
             WRITE(6,*)'ADJUST: P(',J,') =',P(J),' - zero
     + fractions cause crashes'
             WRITE(6,*)'ADJUST: dont set starting values
     + to zero - if you didnt then'
             WRITE(6,*)'ADJUST: try setting limits on the P(J)
     + (0<P(J)<1) before any SCAN or MIGRad calls'
          ENDIF
20     CONTINUE

       TMIN=-1./WPMAX


       KZERO=IMAX
       IF(APMAXS.EQ.0)THEN
          AKI=HI(IDD,I)/(1+WPMAX)
          WKI=HI(IDW(KZERO),I)
          DO 30 J=1,NMCSRC
             IF(J.NE.KZERO)THEN
                WJI=HI(IDW(J),I)
                IF(WKI*P(KZERO).NE.WJI*P(J))THEN
                   AKI=AKI-HI(IDM(J),I)*WJI*P(J)/
     +                     (WKI*P(KZERO)-WJI*P(J))
                ENDIF
             ENDIF
30        CONTINUE
          IF(AKI.GT.0)THEN
             AKI=AKI/NPMAX
             TI=TMIN
             RETURN
          ENDIF
       ENDIF
       KZERO=0

C Now zero cases have been sorted out, solve for the ti using
C Newton's method
       TI=0.

11     CONTINUE

C check starting value for ti is sensible, if not then zero is a good
C place to start.
       IF(TI.GE.1.0.OR.TI.LE.TMIN)THEN
          STEP=STEP/10.0
          TI=0.
       ENDIF
C value of function
       FUNC=-HI(IDD,I)/(1.0D0-TI)
       DERIV=FUNC/(1.0D0-TI)

       DO 40 J=1,NMCSRC
          WJI=HI(IDW(J),I)
          D=1.0D0/(TI+1.0D0/(WJI*P(J)))
          FUNC=FUNC+HI(IDM(J),I)*D
          DERIV=DERIV-HI(IDM(J),I)*D**2
40     CONTINUE

C check convergence - require either F small or step small
       IF(ABS(FUNC).GT.1.0D-12)THEN
          DELTA=-FUNC/DERIV
C limit the step size
          IF(DABS(DELTA).GT.STEP)THEN
             DELTA=DSIGN(STEP,DELTA)
          ENDIF
          TI=TI+DELTA
          IF(DABS(DELTA).GT.1.0D-13)THEN
             GOTO 11
          ENDIF
       ENDIF

       RETURN
       END

*CMZ :          12/06/95  17.50.13  by  Unknown
*CMZ :  4.20/09 24/09/93  16.55.42  by  Christine Beeston
*-- Author :    Roger Barlow, Christine Beeston   24/09/93

      SUBROUTINE HMCINI(IDDATA,IDMC,IDWT,NSRC,CHOPT,IERR)

* Initialisation routine for function HMCLNL, needs to be called
* each time a new set of histograms is introduced (generally once at the
* beginning of each fit).  Performs some error checking and sets up a dummy
* weight histogram if necessary.  Normalises the weight histograms (once
* only, so the same weight histogram may be used for 2 simultaneous
* fits and calculates the overall normalisation constants bj - these are
* stored with the histograms using HNORMA.  Ids of the weight histograms
* are returned in IDWT.

* Input parameters:
* ----------------
* IDDATA Data histogram identifier.
* IDMC Array of dimension NMCSRC containing Monte Carlo histogram
* identifiers.
* IDWT Array of dimension NMCSRC containing weight histogram
* identifiers. ('W' option only).
* NMCSRC Number of Monte Carlo sources.
* CHOPT...
* 'W' Use the weight histograms provided.  For non existent weight
* histograms, and if the W option is not requested, a dummy weight histogram
* in which all entries are 1 is booked.
* IERR - 0 if initialisation fine, 1 if parameters sent to HMCINI were
* not useable (e.g. number of sources 0 or 1)

*KEEP,HCBOOK.
      INTEGER     NWPAW,IXPAWC,IHDIV,IXHIGZ,IXKU,        LMAIN
      REAL                                       FENC   ,      HCV
      COMMON/PAWC/NWPAW,IXPAWC,IHDIV,IXHIGZ,IXKU,FENC(5),LMAIN,HCV(9989)
      INTEGER   IQ        ,LQ
      REAL            Q
      DIMENSION IQ(2),Q(2),LQ(8000)
      EQUIVALENCE (LQ(1),LMAIN),(IQ(1),LQ(9)),(Q(1),IQ(1))
      INTEGER       HVERSN,IHWORK,LHBOOK,LHPLOT,LGTIT,LHWORK,
     +LCDIR,LSDIR,LIDS,LTAB,LCID,LCONT,LSCAT,LPROX,LPROY,LSLIX,
     +LSLIY,LBANX,LBANY,LPRX,LPRY,LFIX,LLID,LR1,LR2,LNAME,LCHAR,LINT,
     +LREAL,LBLOK,LLBLK,LBUFM,LBUF,LTMPM,LTMP,LTMP1,LHPLIP,LHDUM,
     +LHFIT,LFUNC,LHFCO,LHFNA,LCIDN
      COMMON/HCBOOK/HVERSN,IHWORK,LHBOOK,LHPLOT,LGTIT,LHWORK,
     +LCDIR,LSDIR,LIDS,LTAB,LCID,LCONT,LSCAT,LPROX,LPROY,LSLIX,
     +LSLIY,LBANX,LBANY,LPRX,LPRY,LFIX,LLID,LR1,LR2,LNAME,LCHAR,LINT,
     +LREAL,LBLOK,LLBLK,LBUFM,LBUF,LTMPM,LTMP,LTMP1,LHPLIP,LHDUM(9),
     +LHFIT,LFUNC,LHFCO,LHFNA,LCIDN
*
      INTEGER   KNCX   ,KXMIN  ,KXMAX  ,KMIN1  ,KMAX1 ,KNORM  , KTIT1,
     +          KNCY   ,KYMIN  ,KYMAX  ,KMIN2  ,KMAX2 ,KSCAL2 , KTIT2,
     +          KNBIT  ,KNOENT ,KSTAT1 ,KNSDIR  ,KNRH ,
     +          KCON1  ,KCON2  ,KBITS  ,KNTOT
      PARAMETER(KNCX=3,KXMIN=4,KXMAX=5,KMIN1=7,KMAX1=8,KNORM=9,KTIT1=10,
     +          KNCY=7,KYMIN=8,KYMAX=9,KMIN2=6,KMAX2=10,KSCAL2=11,
     +          KTIT2=12,KNBIT=1,KNOENT=2,KSTAT1=3,KNSDIR=5,KNRH=6,
     +          KCON1=9,KCON2=3,KBITS=1,KNTOT=2)
*
*KEEP,HCBITS.
      INTEGER           I1,   I2,   I3,   I4,   I5,   I6,   I7,   I8,
     +                  I9,   I10,  I11,  I12,  I13,  I14,  I15,  I16,
     +I17,  I18,  I19,  I20,  I21,  I22,  I23,  I24,  I25,  I26,  I27,
     +I28,  I29,  I30,  I31,  I32,  I33,  I34,  I35,  I123, I230
      COMMON / HCBITS  / I1,   I2,   I3,   I4,   I5,   I6,   I7,   I8,
     +                  I9,   I10,  I11,  I12,  I13,  I14,  I15,  I16,
     +I17,  I18,  I19,  I20,  I21,  I22,  I23,  I24,  I25,  I26,  I27,
     +I28,  I29,  I30,  I31,  I32,  I33,  I34,  I35,  I123, I230
*
*KEEP,HCMCPM.

C some maximum values - max number of bins and dimensions
       INTEGER NSRCMX,NSRCMN,NOPTS
       PARAMETER(NSRCMX=20,NSRCMN=2,NOPTS=7)

C Numbers of monte carlo and data events, links to data, MC and weight
C histograms, Number of mc sources, options, total number of bins.
C normalisation constants for weight histograms

       DOUBLE PRECISION BJ(NSRCMX)
       INTEGER NDATEV,NMCEV(NSRCMX),
     + NMCSRC,IOPT(NOPTS),NTOT
       COMMON/HINPUTS/NDATEV,NMCEV,NMCSRC,IOPT,NTOT,BJ
C Histogram IDs
       INTEGER IDD,IDM(NSRCMX),IDW(NSRCMX)
       COMMON/HMCIDS/IDD,IDM,IDW

*KEND.

C Passed to the routine..
       INTEGER IDDATA,IDMC(*),IDWT(*),NSRC,IERR,IBIN
       CHARACTER*(*) CHOPT

       DOUBLE PRECISION WJI,AJI,WASUM

C Some variables for use with HGIVE
       CHARACTER*80 CHTITL
       INTEGER NY,NWT,LOC
       REAL XMI,XMA,YMI,YMA

       LOGICAL LWTHIS
       INTEGER IDWTH,NWTHIS,NCXD,NCXM,JSRC,JWT
       PARAMETER(IDWTH=12345)
       REAL DUMWT,CNORM,FACT
       PARAMETER(DUMWT=1.0)

       INTEGER IRET,INORM,IFIRST
       DATA IFIRST/0/

       LOGICAL HEXIST
       EXTERNAL HEXIST

       IF(IFIRST.EQ.0)THEN
          WRITE(6,*)' *********************************************'
          WRITE(6,*)' *** THERE IS A NEW VERSION OF HMCLNL      ***'
          WRITE(6,*)' *** YOU MUST NOW PASS THE HISTOGRAM IDS   ***'
          WRITE(6,*)' *** TO ALLOW FOR MULTIPLE SIMULTANEOUS    ***'
          WRITE(6,*)' *** FITS - HMCINI MUST ALSO BE CALLED     ***'
          WRITE(6,*)' *** FOR EACH SET OF HISTOGRAMS IF A       ***'
          WRITE(6,*)' *** SIMULTANEOUS FIT IS BEING PERFORMED   ***'
          WRITE(6,*)' *** THE NEW HMCLNL IS USED LIKE THIS:     ***'
          WRITE(6,*)' ***                                       ***'
          WRITE(6,*)' *** RLNL=HMCLNL(IDDATA,IDMC,IDWT,NSRC,PJ) ***'
          WRITE(6,*)' ***                                       ***'
          WRITE(6,*)' *** THE OLD VERSION WILL CONTINUE TO BE   ***'
          WRITE(6,*)' *** AVAILABLE AS HMCLNO(PJ) FOR A COUPLE  ***'
          WRITE(6,*)' *** OF MONTHS.                            ***'
          WRITE(6,*)' ***                                       ***'
          WRITE(6,*)' *** REPORT PROBLEMS TO BEESTON@CERNVM     ***'
          WRITE(6,*)' *********************************************'
          IFIRST=1
       ENDIF
C Initialise error flag
       IERR=0
C put number of sources in common
       NMCSRC=NSRC

       LWTHIS=.FALSE.

C Decode options
       CALL HUOPTC(CHOPT,'FLWSNPE',IOPT)


C Initial checks - first check number of MC histos doesnt exceed NSRCMX
       IF(NMCSRC.GT.NSRCMX)THEN
          CALL HBUG('Number of MC histograms exceeds maximum of 20',
     +'HMCINI',0)
          IERR=1
          GOTO 99
       ENDIF
       IF(NMCSRC.LT.NSRCMN)THEN
          CALL HBUG('Number of MC histograms less than minimum of 2',
     +'HMCINI',0)
          IERR=1
          GOTO 99
       ENDIF

C Next get links to histos, stop if any histo cant be found
       IF(HEXIST(IDDATA))THEN
          IDD=IDDATA
       ELSE
          CALL HBUG('Data histogram missing','HMCINI',IDDATA)
          IERR=1
          GOTO 99
       ENDIF
       DO 140 JSRC=1,NMCSRC
          IF(HEXIST(IDMC(JSRC)))THEN
             IDM(JSRC)=IDMC(JSRC)
          ELSE
             CALL HBUG('MC histogram missing','HMCINI',IDMC(JSRC))
             IERR=1
             GOTO 99
          ENDIF
140    CONTINUE
C check all histos have entries
       NDATEV=NINT(HSUM(IDD))
       IF(NDATEV.LE.0)THEN
          CALL HBUG('Data histogram has no entries - cant fit',
     +'HMCINI',IDD)
          IERR=1
          GOTO 99
       ENDIF
       DO 150 JSRC=1,NMCSRC
          NMCEV(JSRC)=NINT(HSUM(IDM(JSRC)))
          IF(NMCEV(JSRC).LE.0)THEN
             CALL HBUG('MC histogram has no entries - cant fit',
     +'HMCINI',IDM(JSRC))
             IERR=1
             GOTO 99
          ENDIF
150    CONTINUE


C check all histos have the same number of bins

       CALL HGIVE(IDD,CHTITL,NCXD,XMI,XMA,NY,YMI,YMA,NWT,LOC)
       DO 160 JSRC=1,NMCSRC
          CALL HGIVE(IDM(JSRC),CHTITL,NCXM,XMI,XMA,NY,YMI,YMA,NWT,LOC)
          IF(NCXM.NE.NCXD)THEN
             CALL HBUG(
     +'MC histo has different number of bins to data histo',
     +'HMCINI',IDM(JSRC))
             IERR=1
             GOTO 99
          ENDIF
160    CONTINUE
       NTOT=NCXD

C Weight histograms - check existence, if histograms not provided then
C use 12345+jwt (filled with 1).  If 'W' option but no histos provided, use
C 12345 for all but issue warning message
C ..first find an unused id for the dummy weight histogram
       JWT=0
165    IF(HEXIST(IDWTH+JWT))THEN
          JWT=JWT+1
          GOTO 165
       ENDIF

C now check existence of histograms if weight option chosen
       IF(IOPT(3).EQ.1)THEN
          NWTHIS=0
          DO 170 JSRC=1,NMCSRC
             IF(.NOT.HEXIST(IDWT(JSRC)))THEN
                CALL HBUG(
     +'Weight histo does not exist - using default',
     +'HMCINI',IDW(JSRC))
                IDW(JSRC)=IDWTH+JWT
                LWTHIS=.TRUE.
             ELSE
                IDW(JSRC)=IDWT(JSRC)
                NWTHIS=NWTHIS+1
             ENDIF
170       CONTINUE
       ELSE
          DO 180 JSRC=1,NMCSRC
             IDW(JSRC)=IDWTH+JWT
             LWTHIS=.TRUE.
180       CONTINUE
       ENDIF
C Put new wt histo ids in idwt to pass back to user
       DO 185 JSRC=1,NMCSRC
          IDWT(JSRC)=IDW(JSRC)
185    CONTINUE

C book + fill the default weights histogram.
       IF(LWTHIS)THEN
          CALL HBOOK1(IDWTH+JWT,'Dummy weight histo',NTOT,0.5,
     +    FLOAT(NTOT)+0.5,0.)
          DO 190 JBIN=1,NTOT
             CALL HFILL(IDWTH+JWT,FLOAT(JBIN),0.,DUMWT)
190       CONTINUE
       ENDIF

C Check all weights histos have right number of bins
       DO 200 JSRC=1,NMCSRC
          CALL HGIVE(IDW(JSRC),CHTITL,NCXM,XMI,XMA,NY,YMI,YMZ,NWT,LOC)
          IF(NCXM.NE.NTOT)THEN
             CALL HBUG(
     +'Weight histo has different number of bins to data histo',
     +'HMCINI',IDM(JSRC))
             IERR=1
             GOTO 99
          ENDIF
200    CONTINUE

C Now normalise weight histograms - first calculate overall normalisation
C factors
       DO 210 JSRC=1,NMCSRC
          IRET=3
          CALL HLOOP(IDW(JSRC),'HMCINI',IRET)
          IF(IRET.EQ.0)THEN
             CALL HBUG(
     +'Weight histogram does not exist ','HMCINI',IDW(JSRC))
             STOP
          ELSE
             INORM=JBIT(IQ(LCID+KBITS),18)
          ENDIF

* only normalise histograms once!
          IF(INORM.EQ.0)THEN
            WASUM=0.0D0
            DO 220 IBIN=1,NTOT
               WJI=HI(IDW(JSRC),IBIN)
               AJI=HI(IDM(JSRC),IBIN)
               WASUM=WASUM+WJI*AJI
220         CONTINUE
            BJ(JSRC)=WASUM/FLOAT(NMCEV(JSRC))

            WRITE(6,*)'****************************************'
            WRITE(6,*)'WASUM IS',WASUM
            WRITE(6,*)'NMCEV IS',NMCEV(JSRC)
            WRITE(6,*)'NORM FACTOR',BJ(JSRC)

            IF(BJ(JSRC).LE.0.0)THEN
               IERR=1
               CALL HBUG(
     +'Normalisation constant for weight histo negative or zero',
     +'HMCINI',IDW(JSRC))
               GOTO 99
            ENDIF
            CNORM=1.0/(BJ(JSRC))
            CALL HOPERA(IDW(JSRC),'+',IDW(JSRC),IDW(JSRC),CNORM,0.0)
* musnt pass double precision to hnorma
            FACT=BJ(JSRC)
C Pack normalisation constant bj in histogram using HNORMA
            CALL HNORMA(IDW(JSRC),FACT)

          ENDIF
210    CONTINUE

99     CONTINUE



       RETURN
       END

*CMZ :          12/04/94  12.48.47  by  Unknown
*CMZ :  4.21/06 24/11/93  16.37.33  by  Christine Beeston
*CMZ :  4.20/13 20/10/93  15.27.09  by  Christine Beeston
*CMZ :  4.20/10 07/10/93  15.29.59  by  Fons Rademakers
*CMZ :  4.20/09 24/09/93  16.30.17  by  Roger Barlow, Christine Beeston
*-- Author :    Roger Barlow, Christine Beeston   24/09/93

       SUBROUTINE HMCMLL(IDDATA,IDMC,IDWT,NSRC,CHOPT,
     +                   IFIXMC,FIXFRA,FLIM,START,STEP,UP,
     +                   ANSMIG,DANSMIG)

* Fits the given Monte Carlo distributions to the data distribution, using
* a binned maximum likelihood fit which includes the effect of both data and
* Monte Carlo statistics, and allows weights to be
* provided for each Monte Carlo distribution.  The data and Monte Carlo
* distributions must be presented in 1 dimensional histograms.
* The best estimate of the fraction of each Monte
* Carlo distribution present in the data distribution is returned, with an
* error estimate where required.

* Input parameters:
* -----------------
* IDDATA Data histogram identifier.
* IDMC Array of dimension NMCSRC containing Monte Carlo histogram
* identifiers.
* IDWT Array of dimension NMCSRC containing weight histogram
* identifiers. ('W' option only).
* NMCSRC Number of Monte Carlo sources.
* CHOPT - character options as follows  1 2 3 4 5 6 7
*                                       F L W S N P E
* 'F' Fix one or more of the fractions.
* 'L' Set limits on the fractions as given in FLIM.
* 'W' Use the weight histograms provided.  For non existent weight
* histograms, and if the W option is not requested, a dummy weight histogram
* in which all entries are 1 is booked.
* 'S' Scan the likelihood function with respect to each fit parameter,
* before and after the fit.  If the 'N' option is specified, the function
* will only be scanned once for each parameter.
* 'N' Do not perform the fit.
* 'P' Use the parameter start points and initial step sizes
* provided in START and STEP.  If negative values are provided the
* defaults are used.
* If the P option is
* not specified then the start point for each free parameter is
* 1.-sum of fixed fractions/nsrc-number of fixed fractions and the initial
* step size is 0.01.
* 'E' Perform a detailed error analysis using the MINUIT routines
* HESSE and MINOS.

* IFIXMC Array of dimension NMCSRC containing '1' if a
* parameter is to be fixed in the fit, '0' otherwise. ('F' option
* only).
* FRAC Array of dimension NMCSRC with the values at which
* parameters are to be fixed. ('F' option only).
* FLIM Array of dimension (NMCSRC,2) with the lower, then
* upper limits on the parameters. ('L' option only).
* START Array of dimension NMCSRC with the start values for
* each parameter. ('P' option only).
* STEP Array of dimension NMCSRC with initial step sizes for the
* parameters ('P' option only).
* UP 'UP' value for the error estimate ('E' option only).  See the Minuit
* manual for definition of UP.

* Output parameters:
* -----------------
* PAR Array of dimension NMCSRC with the final fitted values
* of the parameters.
* DPAR Array of dimension NMCSRC with the errors
* on the final fitted values of the parameters.

*KEEP,HCMCPM.

C some maximum values - max number of bins and dimensions
       INTEGER NSRCMX,NSRCMN,NOPTS
       PARAMETER(NSRCMX=20,NSRCMN=2,NOPTS=7)

C Numbers of monte carlo and data events, links to data, MC and weight
C histograms, Number of mc sources, options, total number of bins.
C normalisation constants for weight histograms

       DOUBLE PRECISION BJ(NSRCMX)
       INTEGER NDATEV,NMCEV(NSRCMX),
     + NMCSRC,IOPT(NOPTS),NTOT
       COMMON/HINPUTS/NDATEV,NMCEV,NMCSRC,IOPT,NTOT,BJ
C Histogram IDs
       INTEGER IDD,IDM(NSRCMX),IDW(NSRCMX)
       COMMON/HMCIDS/IDD,IDM,IDW

*KEND.

C Passed to the routine..
       INTEGER IDDATA,IDMC(*),IDWT(*),NSRC,IFIXMC(*)
       CHARACTER*(*) CHOPT
       DOUBLE PRECISION FLIM(2,*),START(*),STEP(*),UP,FIXFRA(*)

C Number of free MC sourceS, Sum of fixed MC fractions
       INTEGER MCFREE
       DOUBLE PRECISION FRACSM
       COMMON/HMCFIX/FRACSM,MCFREE
       DOUBLE PRECISION ANSMIG(NSRCMX),DANSMIG(NSRCMX)

       CHARACTER*10 PNAME(NSRCMX)

C Stuff for Minuit
       DOUBLE PRECISION PSTART(NSRCMX),PSTEP(NSRCMX)
       DOUBLE PRECISION ARGLIS(10),BLO,BHI,BND1,BND2,VAL,ERROR

       LOGICAL HEXIST

       INTEGER JSRC,IFAIL,IVARBL,IERR

       EXTERNAL HFCNMIN,HEXIST

       DATA PNAME/'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10',
     +  'P11','P12','P13','P14','P15','P16','P17','P18','P19','P20'/
       DATA DEFSTP/0.01D0/

C Initialisation - HMCINI checks existence and compatibility of
C Data, MC and Weight histograms.
       CALL HMCINI(IDDATA,IDMC,IDWT,NSRC,CHOPT,IERR)
       IF(IERR.EQ.1)THEN
          CALL HBUG('Returning from HMCMLL - HMCINI error',
     +    'HMCMLL',0)
          RETURN
       ENDIF

       CALL MNINIT(5,6,7)

C number of free mc sources
       MCFREE=0
       FRACSM=0.0D0
       DO 90 JSRC=1,NMCSRC
          IF(IOPT(1).EQ.0.OR.IFIXMC(JSRC).EQ.0)THEN
             MCFREE=MCFREE+1
          ELSE
             FRACSM=FRACSM+FIXFRA(JSRC)
          ENDIF
90     CONTINUE
       WRITE(6,1000)MCFREE,NMCSRC-MCFREE
1000   FORMAT(' HMCMLL: You have ',I2,' free fractions and ',
     + I2,' fixed')

C set up start values for parameter start points, limits, step sizes
       DO 100 JSRC =1,NMCSRC
          IF(IOPT(2).EQ.1)THEN
             BLO=FLIM(1,JSRC)
             BHI=FLIM(2,JSRC)
          ELSE
             BLO=0.0D0
             BHI=0.0D0
          ENDIF

          IF(IOPT(6).EQ.1)THEN
             IF(START(JSRC).LT.0.0)THEN
                WRITE(6,*)
     +' HMCMLL: Using default for start value source ',JSRC
                PSTART(JSRC)=(1.-FRACSM)/FLOAT(MCFREE)
             ELSE
                PSTART(JSRC)=START(JSRC)
             ENDIF
             IF(STEP(JSRC).LT.0.0)THEN
                WRITE(6,*)
     +' HMCMLL: Using default for step value source',JSRC
                PSTEP(JSRC)=DEFSTP
             ELSE
                PSTEP(JSRC)=STEP(JSRC)
             ENDIF
          ELSE
             PSTART(JSRC)=(1.-FRACSM)/FLOAT(MCFREE)
             PSTEP(JSRC)=DEFSTP
          ENDIF
          IF(IOPT(1).EQ.1.AND.IFIXMC(JSRC).EQ.1)THEN
             PSTART(JSRC)=FIXFRA(JSRC)
          ENDIF

C Initialise Minuit parameters
          CALL MNPARM(JSRC,PNAME(JSRC),PSTART(JSRC),PSTEP(JSRC),
     +    BLO,BHI,IFAIL)
          IF(IFAIL.NE.0)THEN
             CALL HBUG(
     +'Unable to define a parameter in MNPARM call','HMCMLL',0)
             STOP
          ENDIF
100    CONTINUE

C set title
       CALL MNSETI('*** HBOOK New ll maximisation')

C initialisation call to FCN
       ARGLIS(1)=1.
       CALL MNEXCM(HFCNMIN,'CALL FCN',ARGLIS,1,IFAIL)

C fix chosen fractions
       DO 110 JSRC=1,NMCSRC
          IF(IOPT(1).EQ.1.AND.IFIXMC(JSRC).EQ.1)THEN
             ARGLIS(1)=FLOAT(JSRC)
             CALL MNEXCM(HFCNMIN,'FIX',ARGLIS,1,IFAIL)
             WRITE(6,3000)JSRC,FIXFRA(JSRC)
          ENDIF
110    CONTINUE
3000   FORMAT(' HMCMLL: Fixed parameter number ',I2,' at ',F8.5)

C Set initial UP value to 0.5
       ARGLIS(1)=0.5
       CALL MNEXCM(HFCNMIN,'SET ERR',ARGLIS,1,IFAIL)

C scan log likelihood around start
       IF(IOPT(4).EQ.1)THEN
          CALL MNEXCM(HFCNMIN,'SCAN',ARGLIS,0,IFAIL)
       ENDIF

C maximise log - likelihood
       IF(IOPT(5).NE.1)THEN
          IFAIL=0
          CALL MNEXCM(HFCNMIN,'MIGRAD',ARGLIS,0,IFAIL)
          IF(IFAIL.NE.0)THEN
             WRITE(6,3500)IFAIL
          ENDIF
3500      FORMAT(' HMCMLL: MIGRAD ERROR FLAG ',I2)

C scan log likelihood around end point
          IF(IOPT(4).EQ.1)THEN
             CALL MNEXCM(HFCNMIN,'SCAN',ARGLIS,0,IFAIL)
          ENDIF
       ENDIF

C calculate minos errors
C First set UP value as appropriate for  parameters
       IF(IOPT(7).EQ.1)THEN
          IFAIL=0
          IF(UP.LE.0)THEN
             CALL HBUG(
     +'+Negative or zero UP value given - using 0.5',
     +       'HMCMLL',0)
             ARGLIS(1)=0.5
          ELSE
             ARGLIS(1)=UP
          ENDIF
          CALL MNEXCM(HFCNMIN,'SET ERR',ARGLIS,1,IFAIL)
          IF(IFAIL.EQ.0)THEN
             WRITE(6,4000)ARGLIS(1),MCFREE
          ENDIF
4000      FORMAT(' HMCMLL: SET UP VALUE TO ',F5.2,/,
     +           ' HMCMLL: FOR ',I2,' FREE PARAMETERS')

          IFAIL=0
          CALL MNEXCM(HFCNMIN,'MINOS',ARGLIS,0,IFAIL)
          IF(IFAIL.NE.0)THEN
             WRITE(6,5000)IFAIL
          ENDIF
5000      FORMAT(' HMCMLL: MINOS ERROR FLAG ',I2)
       ENDIF

C termination
       DO 130 JSRC=1,NMCSRC
          CALL MNPOUT(JSRC,PNAME(JSRC),VAL,ERROR,BND1,BND2,IVARBL)
          ANSMIG(JSRC)=VAL
          DANSMIG(JSRC)=ERROR
130    CONTINUE

       RETURN
       END

*CMZ :          26/04/94  22.25.05  by  Unknown
*CMZ :  4.21/00 23/10/93  17.02.51  by  Rene Brun
*CMZ :  4.20/09 24/09/93  16.31.04  by  Roger Barlow, Christine Beeston
*-- Author :    Roger Barlow, Christine Beeston   24/09/93

       SUBROUTINE HFCNMIN(NPAR,GRAD,FVAL,XVAL,IFLAG,FUTIL)

*KEEP,HCMCPM.

C some maximum values - max number of bins and dimensions
       INTEGER NSRCMX,NSRCMN,NOPTS
       PARAMETER(NSRCMX=20,NSRCMN=2,NOPTS=7)

C Numbers of monte carlo and data events, links to data, MC and weight
C histograms, Number of mc sources, options, total number of bins.
C normalisation constants for weight histograms

       DOUBLE PRECISION BJ(NSRCMX)
       INTEGER NDATEV,NMCEV(NSRCMX),
     + NMCSRC,IOPT(NOPTS),NTOT
       COMMON/HINPUTS/NDATEV,NMCEV,NMCSRC,IOPT,NTOT,BJ
C Histogram IDs
       INTEGER IDD,IDM(NSRCMX),IDW(NSRCMX)
       COMMON/HMCIDS/IDD,IDM,IDW

*KEND.

       INTEGER NPAR,IFLAG
       DOUBLE PRECISION GRAD(*),FVAL,XVAL(*),FUTIL
       DOUBLE PRECISION HMCLNL
       EXTERNAL HMCLNL

       IF(IFLAG.EQ.1)THEN
          RETURN
       ENDIF

       FVAL=-HMCLNL(IDD,IDM,IDW,NMCSRC,XVAL)

C termination code here.
       IF(IFLAG.EQ.3)THEN
          WRITE(6,*)'HFCNMIN: FCN termination for MINUIT'
          WRITE(6,*)'HFCNMIN: END VALUES OF PARAMETERS'
          WRITE(6,*)'HFCNMIN: ',(XVAL(JSRC),JSRC=1,NMCSRC)
       ENDIF

       RETURN
       END
*CMZ :          25/02/97  20.21.34  by  Unknown
*CMZ :  4.21/06 01/12/93  15.11.40  by  Christine Beeston
*CMZ :  4.20/13 19/10/93  15.33.26  by  Roger Barlow, Christine Beeston
*CMZ :  4.20/10 28/09/93  16.50.28  by  Rene Brun
*CMZ :  4.20/09 24/09/93  16.33.59  by  Roger Barlow, Christine Beeston
*-- Author :    Roger Barlow, Christine Beeston   24/09/93

       DOUBLE PRECISION FUNCTION HMCLNL(IDDATA,IDMC,IDWT,NSRC,PJ)

* HMCLNL returns the log likelihood (including effect of both data
* and Monte Carlo statistics) that the data distribution arose from a
* distribution given by combining the Monte Carlo distributions, weighted
* by the weights provided, using the fractions given in FRAC.
* HMCINI must be called before this function may be used.
*
* Input parameters:
* ----------------
* FRAC: Double precision array of dimension NMCSRC containing the
* fraction of each Monte Carlo distribution you wish to assume is in the data
* distribution, in order to calculate the log likelihood.

*KEEP,HCBOOK.
      INTEGER     NWPAW,IXPAWC,IHDIV,IXHIGZ,IXKU,        LMAIN
      REAL                                       FENC   ,      HCV
      COMMON/PAWC/NWPAW,IXPAWC,IHDIV,IXHIGZ,IXKU,FENC(5),LMAIN,HCV(9989)
      INTEGER   IQ        ,LQ
      REAL            Q
      DIMENSION IQ(2),Q(2),LQ(8000)
      EQUIVALENCE (LQ(1),LMAIN),(IQ(1),LQ(9)),(Q(1),IQ(1))
      INTEGER       HVERSN,IHWORK,LHBOOK,LHPLOT,LGTIT,LHWORK,
     +LCDIR,LSDIR,LIDS,LTAB,LCID,LCONT,LSCAT,LPROX,LPROY,LSLIX,
     +LSLIY,LBANX,LBANY,LPRX,LPRY,LFIX,LLID,LR1,LR2,LNAME,LCHAR,LINT,
     +LREAL,LBLOK,LLBLK,LBUFM,LBUF,LTMPM,LTMP,LTMP1,LHPLIP,LHDUM,
     +LHFIT,LFUNC,LHFCO,LHFNA,LCIDN
      COMMON/HCBOOK/HVERSN,IHWORK,LHBOOK,LHPLOT,LGTIT,LHWORK,
     +LCDIR,LSDIR,LIDS,LTAB,LCID,LCONT,LSCAT,LPROX,LPROY,LSLIX,
     +LSLIY,LBANX,LBANY,LPRX,LPRY,LFIX,LLID,LR1,LR2,LNAME,LCHAR,LINT,
     +LREAL,LBLOK,LLBLK,LBUFM,LBUF,LTMPM,LTMP,LTMP1,LHPLIP,LHDUM(9),
     +LHFIT,LFUNC,LHFCO,LHFNA,LCIDN
*
      INTEGER   KNCX   ,KXMIN  ,KXMAX  ,KMIN1  ,KMAX1 ,KNORM  , KTIT1,
     +          KNCY   ,KYMIN  ,KYMAX  ,KMIN2  ,KMAX2 ,KSCAL2 , KTIT2,
     +          KNBIT  ,KNOENT ,KSTAT1 ,KNSDIR  ,KNRH ,
     +          KCON1  ,KCON2  ,KBITS  ,KNTOT
      PARAMETER(KNCX=3,KXMIN=4,KXMAX=5,KMIN1=7,KMAX1=8,KNORM=9,KTIT1=10,
     +          KNCY=7,KYMIN=8,KYMAX=9,KMIN2=6,KMAX2=10,KSCAL2=11,
     +          KTIT2=12,KNBIT=1,KNOENT=2,KSTAT1=3,KNSDIR=5,KNRH=6,
     +          KCON1=9,KCON2=3,KBITS=1,KNTOT=2)
*
*KEEP,HCBITS.
      INTEGER           I1,   I2,   I3,   I4,   I5,   I6,   I7,   I8,
     +                  I9,   I10,  I11,  I12,  I13,  I14,  I15,  I16,
     +I17,  I18,  I19,  I20,  I21,  I22,  I23,  I24,  I25,  I26,  I27,
     +I28,  I29,  I30,  I31,  I32,  I33,  I34,  I35,  I123, I230
      COMMON / HCBITS  / I1,   I2,   I3,   I4,   I5,   I6,   I7,   I8,
     +                  I9,   I10,  I11,  I12,  I13,  I14,  I15,  I16,
     +I17,  I18,  I19,  I20,  I21,  I22,  I23,  I24,  I25,  I26,  I27,
     +I28,  I29,  I30,  I31,  I32,  I33,  I34,  I35,  I123, I230
*
*KEEP,HCMCPM.

C some maximum values - max number of bins and dimensions
       INTEGER NSRCMX,NSRCMN,NOPTS
       PARAMETER(NSRCMX=20,NSRCMN=2,NOPTS=7)

C Numbers of monte carlo and data events, links to data, MC and weight
C histograms, Number of mc sources, options, total number of bins.
C normalisation constants for weight histograms

       DOUBLE PRECISION BJ(NSRCMX)
       INTEGER NDATEV,NMCEV(NSRCMX),
     + NMCSRC,IOPT(NOPTS),NTOT
       COMMON/HINPUTS/NDATEV,NMCEV,NMCSRC,IOPT,NTOT,BJ
C Histogram IDs
       INTEGER IDD,IDM(NSRCMX),IDW(NSRCMX)
       COMMON/HMCIDS/IDD,IDM,IDW

*KEND.

       INTEGER IDDATA,IDMC(NSRCMX),IDWT(NSRCMX),NSRC

       DOUBLE PRECISION ANSWER(NSRCMX), X(NSRCMX), AKI, TI,
     +                  F,AJI,WJI,PJ(NSRCMX),PJSUM,BPJSUM

       INTEGER KZERO,I,J,IRET,IFIRST
       DATA IFIRST/0/

C Warning that this is new version!!
       IF(IFIRST.EQ.0)THEN
       WRITE(6,*)' ************************************************'
       WRITE(6,*)' *** YOU ARE USING THE NEW VERSION            ***'
       WRITE(6,*)' *** OF HMCLNL - HISTO IDS MUST BE PASSED     ***'
       WRITE(6,*)' *** TO ALLOW FOR >1 INDEPENDENT FIT          ***'
       WRITE(6,*)' *** HMCINI MUST BE CALLED ONCE FOR EACH SET  ***'
       WRITE(6,*)' *** OF HISTOGRAMS.  THEN LOG L IS GIVEN BY   ***'
       WRITE(6,*)' ***                                          ***'
       WRITE(6,*)' ***  RLNL = HMCLNL(IDDATA,IDMC,IDWT,NSRC,PJ) ***'
       WRITE(6,*)' ***                                          ***'
       WRITE(6,*)' ***  THE OLD VERSION OF HMCLNL IS STILL      ***'
       WRITE(6,*)' ***  AVAILABLE - HMCLNO(PJ).  THIS WILL      ***'
       WRITE(6,*)' ***  BE REMOVED IN A COUPLE OF MONTHS        ***'
       WRITE(6,*)' ***  REPORT PROBLEMS TO Roger.barlow@cern.ch ***'
       WRITE(6,*)' ************************************************'
       IFIRST=1
       ENDIF

C put ids in common
       IDD=IDDATA
       NMCSRC=NSRC
       DO 10 J=1,NMCSRC
          IDM(J)=IDMC(J)
          IDW(J)=IDWT(J)
10     CONTINUE

C retrieve histogram normalisation factors(bj) from weight histograms
       DO 20 J=1,NMCSRC
          BJ(J)=1.
          IRET=3
          CALL HLOOP(IDW(J),'HMCLNL',IRET)
          IF(IRET.EQ.0)THEN
             STOP
          ENDIF
          IF(I18.NE.0)THEN
             BJ(J)=Q(LCID+KNORM)
          ENDIF
20     CONTINUE
C get total numbers of events in histograms - ignoring under + overflows
       NDATEV=NINT(HSUM(IDD))
       DO 25 J=1,NMCSRC
          NMCEV(J)=NINT(HSUM(IDM(J)))
25     CONTINUE

C first convert Pj to Pj'
       PJSUM=0.0D0
       BPJSUM=0.0D0
       DO 30 J=1,NMCSRC
          PJSUM=PJSUM+PJ(J)
          BPJSUM=BPJSUM+BJ(J)*PJ(J)
30     CONTINUE
       DO 40 J=1,NMCSRC
          ANSWER(J)=BJ(J)*PJ(J)*PJSUM/BPJSUM
40     CONTINUE

C convert ANSWER to normalised fractions in X
       DO 50 J=1,NMCSRC
          X(J)=ANSWER(J)*NDATEV/NMCEV(J)
50     CONTINUE

       HMCLNL=0.0D0

       DO 80 I=1,NTOT

          CALL HADJUST(TI,I,X,KZERO,AKI)
          F=0.0
          DO 70 J=1,NMCSRC
             WJI=HI(IDW(J),I)
C now get AJI from the aji and the Ti
             IF(KZERO.NE.0.AND.X(J).EQ.X(KZERO))THEN
                AJI=AKI
             ELSE
                AJI=0
                IF(NINT(HI(IDM(J),I)).NE.0)THEN
                   AJI=HI(IDM(J),I)/(1.0D0+X(J)*WJI*TI)
                ENDIF
             ENDIF

             F=F+X(J)*WJI*AJI

             IF(NINT(HI(IDM(J),I)).NE.0)THEN
                IF(AJI.LE.0.0)THEN
                   CONTINUE
                ELSE
                   HMCLNL=HMCLNL+HI(IDM(J),I)*DLOG(AJI)
                ENDIF
             ENDIF

             HMCLNL=HMCLNL-AJI
70        CONTINUE

          IF(NINT(HI(IDD,I)).NE.0)THEN
             HMCLNL=HMCLNL+HI(IDD,I)*DLOG(F)
          ENDIF

          HMCLNL=HMCLNL-F

80     CONTINUE

       RETURN
       END
*CMZ :          26/04/94  21.56.24  by  Unknown
*CMZ :  4.21/06 01/12/93  15.11.40  by  Christine Beeston
*CMZ :  4.20/13 19/10/93  15.33.26  by  Roger Barlow, Christine Beeston
*CMZ :  4.20/10 28/09/93  16.50.28  by  Rene Brun
*CMZ :  4.20/09 24/09/93  16.33.59  by  Roger Barlow, Christine Beeston
*-- Author :    Roger Barlow, Christine Beeston   24/09/93

       DOUBLE PRECISION FUNCTION HMCLNO(PJ)

* HMCLNO returns the log likelihood (including effect of both data
* and Monte Carlo statistics) that the data distribution arose from a
* distribution given by combining the Monte Carlo distributions, weighted
* by the weights provided, using the fractions given in FRAC.
* HMCINI must be called before this function may be used.
*
* Input parameters:
* ----------------
* FRAC: Double precision array of dimension NMCSRC containing the
* fraction of each Monte Carlo distribution you wish to assume is in the data
* distribution, in order to calculate the log likelihood.

*KEEP,HCBOOK.
      INTEGER     NWPAW,IXPAWC,IHDIV,IXHIGZ,IXKU,        LMAIN
      REAL                                       FENC   ,      HCV
      COMMON/PAWC/NWPAW,IXPAWC,IHDIV,IXHIGZ,IXKU,FENC(5),LMAIN,HCV(9989)
      INTEGER   IQ        ,LQ
      REAL            Q
      DIMENSION IQ(2),Q(2),LQ(8000)
      EQUIVALENCE (LQ(1),LMAIN),(IQ(1),LQ(9)),(Q(1),IQ(1))
      INTEGER       HVERSN,IHWORK,LHBOOK,LHPLOT,LGTIT,LHWORK,
     +LCDIR,LSDIR,LIDS,LTAB,LCID,LCONT,LSCAT,LPROX,LPROY,LSLIX,
     +LSLIY,LBANX,LBANY,LPRX,LPRY,LFIX,LLID,LR1,LR2,LNAME,LCHAR,LINT,
     +LREAL,LBLOK,LLBLK,LBUFM,LBUF,LTMPM,LTMP,LTMP1,LHPLIP,LHDUM,
     +LHFIT,LFUNC,LHFCO,LHFNA,LCIDN
      COMMON/HCBOOK/HVERSN,IHWORK,LHBOOK,LHPLOT,LGTIT,LHWORK,
     +LCDIR,LSDIR,LIDS,LTAB,LCID,LCONT,LSCAT,LPROX,LPROY,LSLIX,
     +LSLIY,LBANX,LBANY,LPRX,LPRY,LFIX,LLID,LR1,LR2,LNAME,LCHAR,LINT,
     +LREAL,LBLOK,LLBLK,LBUFM,LBUF,LTMPM,LTMP,LTMP1,LHPLIP,LHDUM(9),
     +LHFIT,LFUNC,LHFCO,LHFNA,LCIDN
*
      INTEGER   KNCX   ,KXMIN  ,KXMAX  ,KMIN1  ,KMAX1 ,KNORM  , KTIT1,
     +          KNCY   ,KYMIN  ,KYMAX  ,KMIN2  ,KMAX2 ,KSCAL2 , KTIT2,
     +          KNBIT  ,KNOENT ,KSTAT1 ,KNSDIR  ,KNRH ,
     +          KCON1  ,KCON2  ,KBITS  ,KNTOT
      PARAMETER(KNCX=3,KXMIN=4,KXMAX=5,KMIN1=7,KMAX1=8,KNORM=9,KTIT1=10,
     +          KNCY=7,KYMIN=8,KYMAX=9,KMIN2=6,KMAX2=10,KSCAL2=11,
     +          KTIT2=12,KNBIT=1,KNOENT=2,KSTAT1=3,KNSDIR=5,KNRH=6,
     +          KCON1=9,KCON2=3,KBITS=1,KNTOT=2)
*
*KEEP,HCBITS.
      INTEGER           I1,   I2,   I3,   I4,   I5,   I6,   I7,   I8,
     +                  I9,   I10,  I11,  I12,  I13,  I14,  I15,  I16,
     +I17,  I18,  I19,  I20,  I21,  I22,  I23,  I24,  I25,  I26,  I27,
     +I28,  I29,  I30,  I31,  I32,  I33,  I34,  I35,  I123, I230
      COMMON / HCBITS  / I1,   I2,   I3,   I4,   I5,   I6,   I7,   I8,
     +                  I9,   I10,  I11,  I12,  I13,  I14,  I15,  I16,
     +I17,  I18,  I19,  I20,  I21,  I22,  I23,  I24,  I25,  I26,  I27,
     +I28,  I29,  I30,  I31,  I32,  I33,  I34,  I35,  I123, I230
*
*KEEP,HCMCPM.

C some maximum values - max number of bins and dimensions
       INTEGER NSRCMX,NSRCMN,NOPTS
       PARAMETER(NSRCMX=20,NSRCMN=2,NOPTS=7)

C Numbers of monte carlo and data events, links to data, MC and weight
C histograms, Number of mc sources, options, total number of bins.
C normalisation constants for weight histograms

       DOUBLE PRECISION BJ(NSRCMX)
       INTEGER NDATEV,NMCEV(NSRCMX),
     + NMCSRC,IOPT(NOPTS),NTOT
       COMMON/HINPUTS/NDATEV,NMCEV,NMCSRC,IOPT,NTOT,BJ
C Histogram IDs
       INTEGER IDD,IDM(NSRCMX),IDW(NSRCMX)
       COMMON/HMCIDS/IDD,IDM,IDW

*KEND.


       DOUBLE PRECISION ANSWER(NSRCMX), X(NSRCMX), AKI, TI,
     +                  F,AJI,WJI,PJ(NSRCMX),PJSUM,BPJSUM

       INTEGER KZERO,I,J,IRET,IFIRST
       DATA IFIRST/0/

       IF(IFIRST.EQ.0)THEN
       WRITE(6,*)' ************************************************'
       WRITE(6,*)' ***   THIS IS THE OLD VERSION OF             ***'
       WRITE(6,*)' ***   HMCLNL:  IT WILL BE DELETED IN         ***'
       WRITE(6,*)' ***   A COUPLE OF MONTHS.  THE NEW           ***'
       WRITE(6,*)' ***   VERSION IS NO SLOWER AND               ***'
       WRITE(6,*)' ***   ALLOWS FOR MULTIPLE SIMULTANEOUS       ***'
       WRITE(6,*)' ***   FITS.  THE HISTOGRAM IDS MUST BE       ***'
       WRITE(6,*)' ***   PASSED AS FOLLOWS:                     ***'
       WRITE(6,*)' ***                                          ***'
       WRITE(6,*)' ***   RLNL=HMCLNL(IDDATA,IDMC,IDWT,NSRC,PJ)  ***'
       WRITE(6,*)' ***                                          ***'
       WRITE(6,*)' *** PLEASE REPORT PROBLEMS TO BEESTON@CERNVM ***'
       WRITE(6,*)' ************************************************'
       IFIRST=1
       ENDIF

C retrieve histogram normalisation factors(bj) from weight histograms
       DO 20 J=1,NMCSRC
          BJ(J)=1.
          IRET=3
          CALL HLOOP(IDW(J),'HMCLNO',IRET)
          IF(IRET.EQ.0)THEN
             STOP
          ENDIF
          IF(I18.NE.0)THEN
             BJ(J)=Q(LCID+KNORM)
          ENDIF
20     CONTINUE
C get total numbers of events in histograms - ignoring under + overflows
       NDATEV=NINT(HSUM(IDD))
       DO 25 J=1,NMCSRC
          NMCEV(J)=NINT(HSUM(IDM(J)))
25     CONTINUE

C first convert Pj to Pj'
       PJSUM=0.0D0
       BPJSUM=0.0D0
       DO 30 J=1,NMCSRC
          PJSUM=PJSUM+PJ(J)
          BPJSUM=BPJSUM+BJ(J)*PJ(J)
30     CONTINUE

       DO 40 J=1,NMCSRC
          ANSWER(J)=BJ(J)*PJ(J)*PJSUM/BPJSUM
40     CONTINUE

C convert ANSWER to normalised fractions in X
       DO 50 J=1,NMCSRC
          X(J)=ANSWER(J)*NDATEV/NMCEV(J)
50     CONTINUE

       HMCLNO=0.0D0

       DO 80 I=1,NTOT

          CALL HADJUST(TI,I,X,KZERO,AKI)
          F=0.0
          DO 70 J=1,NMCSRC
             WJI=HI(IDW(J),I)
C now get AJI from the aji and the Ti
             IF(KZERO.NE.0.AND.X(J).EQ.X(KZERO))THEN
                AJI=AKI
             ELSE
                AJI=0
                IF(NINT(HI(IDM(J),I)).NE.0)THEN
                   AJI=HI(IDM(J),I)/(1.0D0+X(J)*WJI*TI)
                ENDIF
             ENDIF

             F=F+X(J)*WJI*AJI

             IF(NINT(HI(IDM(J),I)).NE.0)THEN
                IF(AJI.LE.0.0)THEN
                   CONTINUE
                ELSE
                   HMCLNO=HMCLNO+HI(IDM(J),I)*DLOG(AJI)
                ENDIF
             ENDIF

             HMCLNO=HMCLNO-AJI
70        CONTINUE

          IF(NINT(HI(IDD,I)).NE.0)THEN
             HMCLNO=HMCLNO+HI(IDD,I)*DLOG(F)
          ENDIF

          HMCLNO=HMCLNO-F

80     CONTINUE

       RETURN
       END