*
* $Id$
*
* $Log$
* Revision 1.1  1996/01/16 17:08:06  mclareni
* Initial revision
*
*
#include "hbook/pilot.h"
*CMZ :  4.22/11 23/08/94  14.05.57  by  Rene Brun
*-- 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)
 
#include "hbook/hcbook.inc"
#include "hbook/hcmcpm.inc"
 
C Passed to the routine..
       INTEGER IDDATA,IDMC(*),IDWT(*),NSRC,IERR,IBIN
       CHARACTER*(*) CHOPT
 
#if defined(CERNLIB_DOUBLE)
       DOUBLE PRECISION WJI,AJI,WASUM
#endif
#if !defined(CERNLIB_DOUBLE)
       REAL WJI,AJI,WASUM
#endif
 
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
 
       LOGICAL HEXIST
       EXTERNAL HEXIST
       DATA IFIRST/0/
 
       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=0
       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))
            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
 
