                                                                                
*     RUNLOOK file                                                              
                                                                                
                                                                                
                                                                                
                                                                                
      SUBROUTINE FRLOOK                                                         
*                                                                               
*     Copy of LOOK vectors to RUN ntuples                                       
*                                                                               
      CHARACTER*2 CHNUM,RECNUM*8,RECNU2*8, LABEL*4, MSG*80                      
      LOGICAL EXS,OPN                                                           
*     UNCOMA: ntuple data accesible to user ----------------------------        
      COMMON/UNCOMA/IDMB,NY,WT,X,Y(100)                                         
      REAL REC(102),XY(0:100)                                                   
      EQUIVALENCE (REC(1),WT),(XY(0),X)                                         
*     IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background           
*     NY   = number of measured variables                                       
*     WT   = weight                                                             
*     X    = true variable (defined only in Monte Carlo) = XY(0)                
*     Y(j) = measured variables, j = 1 ... NY            = XY(j)                
*     REC(.) = array with: WT, X, Y(.) (by equivalence)                         
*     ==================================================================        
*     UNCOMB: internal variables and arrays ----------------------------        
      PARAMETER (LENGDP=50000)                                                  
      PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3))                                    
      COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU),               
     +     HUP(NHU),ASPMC(NKNMAX),                                              
     +     NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG,                
     +     XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD,                             
     +     NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN,               
     +    WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3),                  
     +    NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM,              
     +    HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32)                           
      REAL STAB(12,33)                                                          
      EQUIVALENCE (JTAB(1,1),STAB(1,1))                                         
*                                                                               
*     ITAB         = index of last ntuple                                       
*     NTAB         = number of different ntuple types                           
*     J/STAB( 1,.) = 1 data  =2 moca  =3 bg                                     
*     J/STAB( 2,.) = hollerith text                                             
*     J/STAB( 3,.) =                                                            
*     J/STAB( 4,.) = Lumi  Lumi  evtnr    (first argument                       
*     J/STAB( 5,.) =                                                            
*     J/STAB( 6,.) =                                                            
*     J/STAB( 7,.) =                                                            
*     J/STAB( 8,.) =                                                            
*     J/STAB( 9,.) = 1.0   wfac2           (second argument)                    
*     J/STAB(10,.) =                                                            
*     J/STAB(11,.) = ifg for LOOK ntuples                                       
*     J/STAB(12,.) = inr for LOOK ntuples                                       
*                                                                               
*     INWT = index weight                                                       
*     INIX = index X                                                            
*                                                                               
*                                                                               
*     NV = number of fit variables                (from text input)             
*     NDV(j) = index of fit variable                   - " -                    
*     NBV(j) = number of bins                          - " -                    
*                                                                               
*     XLOW, XHIG = lower and upper limit of x     (from text input)             
*        mandatory, not modified by program                                     
*                                                                               
*     XMIN(j),XMAX(j) = min and max value of Y(j), determined from              
*        ntuples with positive weight (except background)                       
*     JXY(j) = 0 if XMIN/XMAX undefined                                         
*     JXY(j) = 1 if XMIN/XMAX defined                                           
*                                                                               
*     XYL(j),XYH(j) = limits of histogram, rounded for 120 bins,                
*                     for final histograms                                      
*                                                                               
*     XYFL(K),XYFH(K) = limits as used by unfolding fit                         
*                                                                               
*                                                                               
*     Where are bin limits used:                                                
*                                                                               
*     for unfolding fit, usually with small number of bins: XYFL/XYFH           
*                                                                               
*     for final histograms, with 120 bins: XYL/XYH                              
*                                                                               
*                                                                               
*                                                                               
*                                                                               
*     HIST(.,.,.) = histogram of meaured variables, filled in PASS 3            
*     ==================================================================        
*     UNCOMC: character strings-----------------------------------------        
      CHARACTER*50 TITLE,VTEXT                                                  
      CHARACTER*8 SOPTE                                                         
      COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10)                                 
*     TITLE    =  text as comment for problem                                   
*     VTEXT(j) =  text as comment for variable j                                
*     SOPTE(i) =  selected options                                              
*     ------------------------------------------------------------------        
      SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/                                           
*     ==================================================================        
      COMMON/ABTUPL/LUND(2)                                                     
*     ...                                                                       
      CALL SMTEXT('-> FRLOOK: starting ...')                                    
*     check LOOK vectors                                                        
      NTOT=0                                                                    
      DO ITAB=1,NTAB                                                            
       IF(JTAB(11,ITAB).NE.0) THEN                                              
          IFG=JTAB(11,ITAB)                                                     
          INR=JTAB(12,ITAB)                                                     
          CALL GPARM('VEC',IFG,INR,NVEC,NENT,NWORDS)                            
C         WRITE(*,*) 'VEC',IFG,INR,NVEC,NENT,NWORDS                             
          IF(NVEC.GT.0) THEN                                                    
             NTOT=NTOT+NWORDS+4*NENT                                            
          END IF                                                                
       END IF                                                                   
      END DO                                                                    
      IF(NTOT.EQ.0) GOTO 100                                                    
*     calculate number of records                                               
      NTOT=1+(NTOT*4-1)/4316                                                    
*     define defalt units                                                       
      CALL ABUNIT(51,52)                                                        
*     check files A and B and use FPACK open eventually                         
      DO ILUN=1,2                                                               
       LUNIT=LUND(ILUN)                                                         
C      WRITE(*,*) 'UNIT =',LUNIT                                                
       IF(LUNIT.NE.0) THEN                                                      
*         inquire for given unit                                                
          INQUIRE (UNIT=LUNIT,IOSTAT=IOS,EXIST=EXS,OPENED=OPN)                  
          IF(IOS.NE.0) THEN                                                     
             WRITE(*,*) 'Inquire failed for unit',LUNIT,                        
     +                  ' Error code=',IOS                                      
             STOP                                                               
          END IF                                                                
C         WRITE(*,*) EXS,OPN,' EXS,OPN'                                         
          IF(.NOT.OPN) THEN                                                     
             WRITE(CHNUM ,'(I2)') LUNIT                                         
             WRITE(RECNUM,'(I8)') NTOT                                          
             WRITE(RECNU2,'(I8)') 1+NTOT/5                                      
C            WRITE(*,*) 'OPEN UNIT='//CHNUM//',RECL=4320,ACTION                 
C    +       = READWRITE,ACCESS=SEQ,STATUS=SCR,FORM=BINARY,                     
C    +       NREC='//RECNUM//',NREC2='//RECNU2                                  
             CALL FPARM('OPEN UNIT='//CHNUM//',RECL=4320,ACTION                 
     +       = READWRITE,ACCESS=SEQ,STATUS=SCR,FORM=BINARY,                     
     +       NREC='//RECNUM//',NREC2='//RECNU2)                                 
             CALL FERMES(MSG,0)                                                 
             WRITE(*,*) MSG                                                     
             CALL FERMES(MSG,1)                                                 
             WRITE(*,*) MSG                                                     
          END IF                                                                
       END IF                                                                   
      END DO                                                                    
C     WRITE(*,*) INIX,INWT,' INIX,INWT'                                         
*     transfer LOOK ntuples to RUN                                              
      DO ITAB=1,NTAB                                                            
       IF(JTAB(11,ITAB).NE.0) THEN                                              
          IFG=JTAB(11,ITAB)                                                     
          INR=JTAB(12,ITAB)                                                     
          CALL HLTOCH(JTAB(2,ITAB),1,LABEL,NCH)                                 
          CALL GPARM('VEC',IFG,INR,NVEC,NENT,NWORDS)                            
          IF(NVEC.GT.0) THEN                                                    
             CALL RWDATA('VEC',IFG,INR)                                         
             X=0.0                                                              
   10        CALL GDATA('VEC',IFG,INR,NR,Y,NVEC)                                
             IF(NR.NE.0) THEN                                                   
                WT=Y(INWT)                                                      
                X =Y(INIX)                                                      
                CALL UEVENT(LABEL,NVEC,WT,X,Y)                                  
                GOTO 10                                                         
             END IF                                                             
             CALL PURGDF('VEC',IFG,INR)                                         
          END IF                                                                
       END IF                                                                   
      END DO                                                                    
      CALL UEVENT(LABEL,-1,WT,X,Y)                                              
      CALL SMTEXT('-- FRLOOK: transfer of ntuples finished')                    
  100 RETURN                                                                    
      END                                                                       
                                                                                
      SUBROUTINE RULOOK                                                         
*                                                                               
*     interface between RUN and LOOP for graphical display                      
*                                                                               
*     UNCOMA: ntuple data accesible to user ----------------------------        
      COMMON/UNCOMA/IDMB,NY,WT,X,Y(100)                                         
      REAL REC(102),XY(0:100)                                                   
      EQUIVALENCE (REC(1),WT),(XY(0),X)                                         
*     IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background           
*     NY   = number of measured variables                                       
*     WT   = weight                                                             
*     X    = true variable (defined only in Monte Carlo) = XY(0)                
*     Y(j) = measured variables, j = 1 ... NY            = XY(j)                
*     REC(.) = array with: WT, X, Y(.) (by equivalence)                         
*     ==================================================================        
*     UNCOMB: internal variables and arrays ----------------------------        
      PARAMETER (LENGDP=50000)                                                  
      PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3))                                    
      COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU),               
     +     HUP(NHU),ASPMC(NKNMAX),                                              
     +     NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG,                
     +     XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD,                             
     +     NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN,               
     +    WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3),                  
     +    NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM,              
     +    HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32)                           
      REAL STAB(12,33)                                                          
      EQUIVALENCE (JTAB(1,1),STAB(1,1))                                         
*                                                                               
*     ITAB         = index of last ntuple                                       
*     NTAB         = number of different ntuple types                           
*     J/STAB( 1,.) = 1 data  =2 moca  =3 bg                                     
*     J/STAB( 2,.) = hollerith text                                             
*     J/STAB( 3,.) =                                                            
*     J/STAB( 4,.) = Lumi  Lumi  evtnr    (first argument                       
*     J/STAB( 5,.) =                                                            
*     J/STAB( 6,.) =                                                            
*     J/STAB( 7,.) =                                                            
*     J/STAB( 8,.) =                                                            
*     J/STAB( 9,.) = 1.0   wfac2           (second argument)                    
*     J/STAB(10,.) =                                                            
*     J/STAB(11,.) = ifg for LOOK ntuples                                       
*     J/STAB(12,.) = inr for LOOK ntuples                                       
*                                                                               
*     INWT = index weight                                                       
*     INIX = index X                                                            
*                                                                               
*                                                                               
*     NV = number of fit variables                (from text input)             
*     NDV(j) = index of fit variable                   - " -                    
*     NBV(j) = number of bins                          - " -                    
*                                                                               
*     XLOW, XHIG = lower and upper limit of x     (from text input)             
*        mandatory, not modified by program                                     
*                                                                               
*     XMIN(j),XMAX(j) = min and max value of Y(j), determined from              
*        ntuples with positive weight (except background)                       
*     JXY(j) = 0 if XMIN/XMAX undefined                                         
*     JXY(j) = 1 if XMIN/XMAX defined                                           
*                                                                               
*     XYL(j),XYH(j) = limits of histogram, rounded for 120 bins,                
*                     for final histograms                                      
*                                                                               
*     XYFL(K),XYFH(K) = limits as used by unfolding fit                         
*                                                                               
*                                                                               
*     Where are bin limits used:                                                
*                                                                               
*     for unfolding fit, usually with small number of bins: XYFL/XYFH           
*                                                                               
*     for final histograms, with 120 bins: XYL/XYH                              
*                                                                               
*                                                                               
*                                                                               
*                                                                               
*     HIST(.,.,.) = histogram of meaured variables, filled in PASS 3            
*     ==================================================================        
*     UNCOMC: character strings-----------------------------------------        
      CHARACTER*50 TITLE,VTEXT                                                  
      CHARACTER*8 SOPTE                                                         
      COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10)                                 
*     TITLE    =  text as comment for problem                                   
*     VTEXT(j) =  text as comment for variable j                                
*     SOPTE(i) =  selected options                                              
*     ------------------------------------------------------------------        
      SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/                                           
*     ==================================================================        
      COMMON/RUNRES/XW(527),RSP(2082)                                           
*     ------------------------------------------------------------------        
*     Internal common for matrices, vectors etc for up to 63 knots              
*     and up to 30 data points                                                  
      PARAMETER (NKMAX=63,NDAMAX=30)                                            
      COMMON/UINTER/UM(NKMAX*NKMAX),HE((NKMAX+NKMAX*NKMAX)/2),                  
     +              GR(NKMAX),XX(NKMAX),AM(NKMAX),DE(NKMAX),                    
     +              SCV(NKMAX),                                                 
     +              SCM(NKMAX*NKMAX),                                           
     +              FUN,SF,SD,TAU                                               
*     M=NKN-1                                                                   
*     UM = tranformation matrix                                                 
*     HE = hessian, replaced by covariance matrix                               
*     GR = gradient, replaced by step vector                                    
*     XX = coefficient                                                          
*     AM = amplitudes     M                                                     
*     DE = eigenvalues    M                                                     
*     SC = scratch storage                                                      
*     ------------------------------------------------------------------        
      SAVE /RUNRES/,/UINTER/                                                    
      COMMON/UNCOMD/H(LENGDP)                                                   
      INTEGER NH(10000)                                                         
      EQUIVALENCE (H(1),NH(1))                                                  
      REAL RH(11,2),SH(2),S(4),RSREL(2),RSNUM(2)                                
      REAL BX(2),BY(2)                                                          
      REAL HUT(120),HUQ(120)                                                    
      CHARACTER*7  COLOR(7), CHNUM*3                                            
      DOUBLE PRECISION SUMDB,SUMMC                                              
      DATA COLOR/'WHITE','RED','GREEN','BLUE','YELLOW','MAGENTA','CYAN'/        
*     ...                                                                       
      CALL SMTEXT(' ')                                                          
      CALL SMTEXT('-> Subroutine RULOOK starting ...')                          
C     CALL SMTEXT('RUN - Start of RULOOK   - '//TITLE//'.')                     
C     CALL SMTEXT('=====================')                                      
*                                                                               
*     ------------------------------------------------------------------        
*     plots of group 1: Result of unfolding and related distribution            
*                                                                               
*                                                                               
*     Data points from unfolding                                                
      NN=XW(1)+0.5                                                              
      II=0                                                                      
      DO I=1,NN                                                                 
       II=II+I                                                                  
       DY=SQRT(XW(NN+NN+2+II))                                                  
       YY=XW(NN+2+I)                                                            
       XE=XW(1+I)                                                               
       XR=XW(2+I)                                                               
       CALL SPDATA(101,0,0.5*(XE+XR),YY,0.5*(XR-XE),DY)                         
       CALL SPDATA(102,0,0.5*(XE+XR),YY,0.5*(XR-XE),DY)                         
      END DO                                                                    
      CALL STEXT(101,4,'Solution of unfolding')                                 
*     Unfolding result as curve                                                 
      NKN=RSP(1)+0.5                                                            
      XLOW=RSP(2)                                                               
      XHIG=RSP(3)                                                               
      K=XW(1)+1.5                                                               
      DO I=1,201                                                                
       XF=0.005*(FLOAT(201-I)*XLOW+FLOAT(I-1)*XHIG)                             
       FU=ESPF(XF,RSP(4),NKN,XLOW,XHIG)                                         
       TX=FU*USFUN(XF)                                                          
       CALL SPCURV(102,1,XF,FACLUM*TX)                                          
       TT=ESPF(XF,UM(1+(K-1)*NKN),NKN,XLOW,XHIG)*USFUN(XF)                      
       CALL SPCURV(103,1,XF,FACLUM*TX   )                                       
*     ... plus N+1'th orthogonal function                                       
       CALL SPCURV(103,2,XF,FACLUM*(TX+TT))                                     
       CALL SPCURV(103,3,XF,FACLUM*(TX-TT))                                     
*      unfolding result                                                         
                                                                                
      END DO                                                                    
      CALL SATTR('XY',102,1,'FULL BLUE CURVE')                                  
      CALL SATTR('XY',103,1,'FULL BLUE CURVE')                                  
      CALL SATTR('XY',103,2,'FULL RED  CURVE')                                  
      CALL SATTR('XY',103,3,'FULL RED  CURVE')                                  
      CALL STEXT(103,4,'Total solution +- next orthogonal fcn')                 
                                                                                
*     figure 104: function USFUN                                                
      DO I=1,201                                                                
       XF=0.005*(FLOAT(201-I)*XLOW+FLOAT(I-1)*XHIG)                             
       TX=USFUN(XF)                                                             
       CALL SPCURV(104,1,XF,TX)                                                 
      END DO                                                                    
      CALL SATTR('XY',104,1,'FULL GREEN CURVE')                                 
      CALL STEXT(104,4,'User function USFUN(x)')                                
                                                                                
                                                                                
*     fill array HUT with user function USFUN ... and plot                      
C     DO I=1,NHA                                                                
C      X=XLOW+(FLOAT(I)-0.5)*(XHIG-XLOW)/FLOAT(NHA)                             
C      HUT(I)=USFUN(X)                                                          
C     END DO                                                                    
                                                                                
      DO I=1,201                                                                
       XF=0.005*(FLOAT(201-I)*XLOW+FLOAT(I-1)*XHIG)                             
       TX=ESPF(XF,ASPMC,NKN,XLOW,XHIG)                                          
       CALL SPCURV(105,1,XF,TX)                                                 
       FU=ESPF(XF,RSP(4),NKN,XLOW,XHIG)                                         
       CALL SPCURV(106,1,XF,FACLUM*FU)                                          
      END DO                                                                    
C     CALL LKHIST(105,3,HUH,NHA,XLOW,XHIG)                                      
C     CALL LKHIST(105,4,HUA,NHA,XLOW,XHIG)                                      
      CALL SPHIST(105,1,HUH,NHA,XLOW,XHIG)                                      
      CALL SPHIST(105,2,HUA,NHA,XLOW,XHIG)                                      
      CALL STEXT(105,4,'MOCA input distribution and spline fit')                
      CALL STEXT(106,4,'Unfolding distribution with USFUN     ')                
                                                                                
*     figure 107: weight as a function of x                                     
      CALL SPHIST(107,1,HUW,NHA,XLOW,XHIG)                                      
*                                                                               
*     ------------------------------------------------------------------        
*     Plots of group 2: Orthogonal functions                                    
*                                                                               
      DO K=1,NKN                                                                
       CALL ESPLK(200+K,1,UM(1+(K-1)*NKN),NKN,XLOW,XHIG,0)                      
       CALL SATTR('XY',200+K,1,'FULL BLUE CURVE')                               
       WRITE(CHNUM,'(I3)') K                                                    
       CALL STEXT(200+K,5,CHNUM//'. Orthogonal combination')                    
      END DO                                                                    
*                                                                               
*     ------------------------------------------------------------------        
*     Plots of group 3: Fit variables and B-splines                             
*                                                                               
*     curve of damping factor *10.0                                             
      DO I=1,NKN-1                                                              
       CALL SXY(301,1,FLOAT(I),10.0/(1.0+TAU*DE(I)))                            
      END DO                                                                    
      CALL STEXT(301,1,'Index I            ')                                   
      CALL STEXT(301,2,'Amplitude A(I)     ')                                   
      CALL STEXT(301,4,'Full and damped amplitudes ')                           
      CALL STEXT(301,5,'curve = 10 * damping factor')                           
      CALL SATTR('XY',301,1,'FULL RED  LINE')                                   
*     amplitudes                                                                
      WID=0.3                                                                   
      DO I=1,NKN-1                                                              
       FLI=I                                                                    
       AA=ABS(AM(I))                                                            
       IF(AA.GT.12.0) AA=12.0                                                   
       CALL SXY(301,2,FLI-WID,0.0)                                              
       CALL SXY(301,2,FLI-WID,AA)                                               
       IF(ABS(AM(I)).GT.12.0) AA=AA-1.0                                         
       CALL SXY(301,2,FLI+WID,AA)                                               
       CALL SXY(301,2,FLI+WID,0.0)                                              
       CALL SXY(301,3,FLI,AA)                                                   
      END DO                                                                    
      CALL SATTR('XY',301,2,'FULL BLUE  LINE')                                  
      CALL SATTR('XY',301,3,'TRIA BLUE  SYMB')                                  
                                                                                
      DO I=1,NKN-1                                                              
       FLI=I                                                                    
       AD=ABS(AM(I))/(1.0+TAU*DE(I))                                            
       AA=ABS(AD)                                                               
       IF(AA.GT.12.0) AA=12.0                                                   
C      CALL SXY(301,4,FLI,AA)                                                   
       CALL SXY(301,5,FLI-WID,0.0)                                              
       CALL SXY(301,5,FLI-WID,AA)                                               
       IF(AD.GT.12.0) AA=AA-1.0                                                 
       CALL SXY(301,5,FLI+WID,AA)                                               
       CALL SXY(301,5,FLI+WID,0.0)                                              
      END DO                                                                    
C     CALL SATTR('XY',301,4,'GREEN CIRCLE SYMF')                                
      CALL SATTR('XY',301,5,'FULL GREEN LINE')                                  
*     damped amplitudes                                                         
      DO I=1,NKN-1                                                              
       FLI=I                                                                    
       TD=TAU*DE(I)                                                             
       XNOIS=ABS(AM(I))*SQRT(TD*(2.0+TD))/(1+TD)                                
       IF(XNOIS.LT.11.9) CALL SXY(301,6,FLI,XNOIS)                              
      END DO                                                                    
      CALL SATTR('XY',301,6,'RED  CIRCLE SYMF  ')                               
*     ------------------------------------------------------------------        
*     Plots of group 4: MOCA distributions and B-spline contributions           
*                                                                               
*     sum data + background ...                                                 
      SUMDB=0.0                                                                 
      DO I=1,NBN                                                                
       SUMDB=SUMDB+H(NPROD+I)+H(NPROD+NBN+I)                                    
      END DO                                                                    
*     ... sum Monte Carlo ...                                                   
      SUMMC=0.0                                                                 
      DO I=1,NBN*NKN                                                            
       SUMMC=SUMMC+H(NPROD+NBN+NBN+I)                                           
      END DO                                                                    
*     ... and determine ratio for normalization of plot                         
      RATIO=SUMDB/SUMMC                                                         
*                                                                               
      NI=1                                                                      
      NJ=1                                                                      
      NK=1                                                                      
      DO L=1,NV                                                                 
       IF(L.EQ.1) NI=NBV(L)                                                     
       IF(L.EQ.2) NJ=NBV(L)                                                     
       IF(L.EQ.3) NK=NBV(L)                                                     
      END DO                                                                    
*     do loop all fit variables                                                 
      DO L=1,NV                                                                 
C      WRITE(*,*) L,NV,NBV(L),'L,NV,NBV(L)'                                     
*      reset total histogram for variable L                                     
       DO I=1,NBV(L)                                                            
        HUQ(I)=0.0                                                              
       END DO                                                                   
      DO IKN=1,NKN                                                              
*           IKN = knot (Bspline) index                                          
       IF(L.EQ.1) THEN                                                          
          DO IJK=1,NBV(L)                                                       
           HUT(IJK)=RATIO*AIJK(IKN,NI,NJ,NK,IJK,-1,-1)                          
           HUQ(IJK)=HUQ(IJK)+HUT(IJK)                                           
          END DO                                                                
       ELSE IF(L.EQ.2) THEN                                                     
          DO IJK=1,NBV(L)                                                       
           HUT(IJK)=RATIO*AIJK(IKN,NI,NJ,NK,-1,IJK,-1)                          
           HUQ(IJK)=HUQ(IJK)+HUT(IJK)                                           
          END DO                                                                
       ELSE IF(L.EQ.3) THEN                                                     
          DO IJK=1,NBV(L)                                                       
           HUT(IJK)=RATIO*AIJK(IKN,NI,NJ,NK,-1,-1,IJK)                          
           HUQ(IJK)=HUQ(IJK)+HUT(IJK)                                           
          END DO                                                                
       END IF                                                                   
*      histogram number 411...414, 421...424, 431...434                         
       IHISTR=400+L*10+MOD(IKN-1,4)+1                                           
*      colour index 1...7                                                       
       ICOLOR=(IKN+3)/4                                                         
       ICOLOR=MOD(ICOLOR-1,7)+1                                                 
C      WRITE(*,*) L,IKN,IHISTR,ICOLOR,' is L IKN, IHISTR,ICOLOR'                
*      contribution to B-spline IKN                                             
C      CALL LKHIST(IHISTR,IKN,HUT,NBV(L), XYL(NDV(L)),XYH(NDV(L)))              
C     WRITE(*,*) 'con ',IHISTR,IKN                                              
C      WRITE(*,222) (HUT(MM),MM=1,NBV(L))                                       
  222  FORMAT(1X,10G12.3)                                                       
C      CALL LKHIST(IHISTR,IKN,HUT,NBV(L),XYFL(L),XYFH(L))                       
       CALL SPHIST(IHISTR,IKN,HUT,NBV(L),XYFL(L),XYFH(L))                       
       CALL SATTR('HS',IHISTR,IKN,'HIST NORMAL'//COLOR(ICOLOR))                 
      END DO                                                                    
*     and sum = total histogram for each figure                                 
      DO J =1,4                                                                 
C      CALL LKHIST(400+L*10+J,0,HUQ,NBV(L), XYL(NDV(L)),XYH(NDV(L)))            
C     WRITE(*,*) 'sum ',400+L*10+J                                              
C      WRITE(*,222) (HUQ(MM),MM=1,NBV(L))                                       
CC     CALL LKHIST(400+L*10+J,0,HUQ,NBV(L),XYFL(L),XYFH(L))                     
       CALL SPHIST(400+L*10+J,0,HUQ,NBV(L),XYFL(L),XYFH(L))                     
       CALL SATTR('HS',400+L*10+J,0,'HIST NORMAL WHITE')                        
      END DO                                                                    
       CALL STEXT(400+L*10+1,5,'Sum and B-splines 1, 5, 9 ...')                 
       CALL STEXT(400+L*10+2,5,'Sum and B-splines 2, 6, 10 ...')                
       CALL STEXT(400+L*10+3,5,'Sum and B-splines 3, 7, 11 ...')                
       CALL STEXT(400+L*10+4,5,'Sum and B-splines 4, 8, 12 ...')                
      END DO                                                                    
*                                                                               
*     Plots of groups 5 to 8 -------------------------------------------        
*                                                                               
      DO J=1,10                                                                 
       IF(JXY(J).NE.0) THEN                                                     
*         event numbers in bins of constant width                               
*         transfer data histogram to LOOK                                       
C         CALL LKHIST(500+J,0,HIST(1,1,J),120,XYL(J),XYH(J))                    
          CALL SPHIST(500+J,0,HIST(1,1,J),120,XYL(J),XYH(J))                    
*         copy sum of MOCA and background ...                                   
          DO I=1,120                                                            
           HSM(I)=HIST(I,2,J)                                                   
          END DO                                                                
*         ... smooth it ...                                                     
          CALL SPLFT(HSM,120)                                                   
*         ... and add histogram to LOOK as curve                                
C         CALL LKHIST(500+J,1,HSM,120,XYL(J),XYH(J))                            
          CALL SPHIST(500+J,1,HSM,120,XYL(J),XYH(J))                            
          CALL SATTR('HS',500+J,1,'CURVE FULL'//COLOR(ICOLOR))                  
*                                                                               
*         determine bins from data                                              
C         CALL GRB(HIST(1,1,J),120)                                             
*         plot expected distribution (= MOCA + background)                      
C         CALL GRC(600+J,1,HIST(1,2,J),120,XYL(J),XYH(J))                       
          CALL SPCURH(600+J,1,HIST(1,2,J),120,XYL(J),XYH(J))                    
*         plot data                                                             
C         CALL GRH(600+J,1,HIST(1,1,J),HIST(1,1,J),120,XYL(J),XYH(J))           
*         histogram condensed with determination of new binning                 
          CALL SPHISC(600+J,1,                                                  
     +             HIST(1,1,J),HIST(1,1,J),120,XYL(J),XYH(J),0)                 
*         calculate difference                                                  
          DO I=1,120                                                            
           HSM(I)=HIST(I,1,J)-HIST(I,2,J)                                       
          END DO                                                                
C         CALL LKHIST(700+J,0,HSM,120,XYL(J),XYH(J))                            
          CALL SPHIST(700+J,0,HSM,120,XYL(J),XYH(J))                            
*         histogram condensed                                                   
          CALL SPHISC(800+J,1,HSM,HIST(I,1,J),120,XYL(J),XYH(J),1)              
       END IF                                                                   
      END DO                                                                    
      CALL SMTEXT('-- Subroutine RULOOK ending.')                               
C     CALL SMTEXT('End of RULOOK ')                                             
      CALL SMTEXT(' ')                                                          
  100 RETURN                                                                    
      END                                                                       
                                                                                
      FUNCTION AIJK(IKN,NI,NJ,NK,LI,LJ,LK)                                      
*                                                                               
*     Extraction of projections from three-dimensional array                    
*                                                                               
*                                                                               
*     SUM =AIJK (IKN,NI,NJ,NK,LI,LJ,LK)                                         
*                                                                               
*                                                                               
*                                                                               
*     IKN = knot index                                                          
*     LI = -1             sum over I=1  . . NI                                  
*     LI =  I             index                                                 
*                                                                               
*     LJ, LK respectively                                                       
*                                                                               
      DOUBLE PRECISION SUM                                                      
*     UNCOMA: ntuple data accesible to user ----------------------------        
      COMMON/UNCOMA/IDMB,NY,WT,X,Y(100)                                         
      REAL REC(102),XY(0:100)                                                   
      EQUIVALENCE (REC(1),WT),(XY(0),X)                                         
*     IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background           
*     NY   = number of measured variables                                       
*     WT   = weight                                                             
*     X    = true variable (defined only in Monte Carlo) = XY(0)                
*     Y(j) = measured variables, j = 1 ... NY            = XY(j)                
*     REC(.) = array with: WT, X, Y(.) (by equivalence)                         
*     ==================================================================        
*     UNCOMB: internal variables and arrays ----------------------------        
      PARAMETER (LENGDP=50000)                                                  
      PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3))                                    
      COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU),               
     +     HUP(NHU),ASPMC(NKNMAX),                                              
     +     NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG,                
     +     XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD,                             
     +     NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN,               
     +    WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3),                  
     +    NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM,              
     +    HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32)                           
      REAL STAB(12,33)                                                          
      EQUIVALENCE (JTAB(1,1),STAB(1,1))                                         
*                                                                               
*     ITAB         = index of last ntuple                                       
*     NTAB         = number of different ntuple types                           
*     J/STAB( 1,.) = 1 data  =2 moca  =3 bg                                     
*     J/STAB( 2,.) = hollerith text                                             
*     J/STAB( 3,.) =                                                            
*     J/STAB( 4,.) = Lumi  Lumi  evtnr    (first argument                       
*     J/STAB( 5,.) =                                                            
*     J/STAB( 6,.) =                                                            
*     J/STAB( 7,.) =                                                            
*     J/STAB( 8,.) =                                                            
*     J/STAB( 9,.) = 1.0   wfac2           (second argument)                    
*     J/STAB(10,.) =                                                            
*     J/STAB(11,.) = ifg for LOOK ntuples                                       
*     J/STAB(12,.) = inr for LOOK ntuples                                       
*                                                                               
*     INWT = index weight                                                       
*     INIX = index X                                                            
*                                                                               
*                                                                               
*     NV = number of fit variables                (from text input)             
*     NDV(j) = index of fit variable                   - " -                    
*     NBV(j) = number of bins                          - " -                    
*                                                                               
*     XLOW, XHIG = lower and upper limit of x     (from text input)             
*        mandatory, not modified by program                                     
*                                                                               
*     XMIN(j),XMAX(j) = min and max value of Y(j), determined from              
*        ntuples with positive weight (except background)                       
*     JXY(j) = 0 if XMIN/XMAX undefined                                         
*     JXY(j) = 1 if XMIN/XMAX defined                                           
*                                                                               
*     XYL(j),XYH(j) = limits of histogram, rounded for 120 bins,                
*                     for final histograms                                      
*                                                                               
*     XYFL(K),XYFH(K) = limits as used by unfolding fit                         
*                                                                               
*                                                                               
*     Where are bin limits used:                                                
*                                                                               
*     for unfolding fit, usually with small number of bins: XYFL/XYFH           
*                                                                               
*     for final histograms, with 120 bins: XYL/XYH                              
*                                                                               
*                                                                               
*                                                                               
*                                                                               
*     HIST(.,.,.) = histogram of meaured variables, filled in PASS 3            
*     ==================================================================        
*     UNCOMC: character strings-----------------------------------------        
      CHARACTER*50 TITLE,VTEXT                                                  
      CHARACTER*8 SOPTE                                                         
      COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10)                                 
*     TITLE    =  text as comment for problem                                   
*     VTEXT(j) =  text as comment for variable j                                
*     SOPTE(i) =  selected options                                              
*     ------------------------------------------------------------------        
      SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/                                           
*     ==================================================================        
      COMMON/UNCOMD/H(LENGDP)                                                   
      INTEGER NH(LENGDP)                                                        
      EQUIVALENCE (H(1),NH(1))                                                  
*     ...                                                                       
C     WRITE(*,111) IKN,NI,NJ,NK,LI,LJ,LK,'IKN,NI,NJ,NK,LI,LJ,LK'                
  111 FORMAT(1X,7I4,1X,A)                                                       
      KS=1                                                                      
      JS=1                                                                      
      IS=1                                                                      
*                                                                               
*     loop                                                                      
*                                                                               
      SUM=0.0D0                                                                 
      K=LK                                                                      
      IF(LK.GT.0) GOTO 20                                                       
      K=KS-1                                                                    
   15 K=K+1                                                                     
                                                                                
   20 J=LJ                                                                      
      IF(LJ.GT.0) GOTO 30                                                       
      J=JS-1                                                                    
   25 J=J+1                                                                     
                                                                                
   30 I=LI                                                                      
      IF(LI.GT.0) GOTO 40                                                       
      I=IS-1                                                                    
   35 I=I+1                                                                     
*                                                                               
   40 IJK=I+(J-1)*NI+(K-1)*NI*NJ                                                
      JN=NH(IJK)                                                                
C     WRITE(*,112) I,J,K,IJK,'I,J,K,IJK',JN                                     
  112 FORMAT(1X,3I4,I7,1X,A)                                                    
      IF(JN.NE.0) SUM=SUM+H(NPROD+NBN*2+(JN-1)*NKN+IKN)                         
*                                                                               
      IF(LI.LT.0.AND.I.LT.NI) GOTO 35                                           
      IF(LJ.LT.0.AND.J.LT.NJ) GOTO 25                                           
      IF(LK.LT.0.AND.K.LT.NK) GOTO 15                                           
*     end loop                                                                  
      AIJK=SUM                                                                  
      SU=SUM                                                                    
C     WRITE(*,101) SUM,IKN,NI,NJ,NK,LI,LJ,LK                                    
*                                                                               
  101 FORMAT(1X,G12.3,7I7,'AIJK,IKN,NI,NJ,NK,LI,LJ,LK')                         
      END                                                                       
      SUBROUTINE ESPLK(IFG,INR,A,N,XA,XB,ND)                                    
C                                                                               
C     SUBROUTINE ESPL                                                           
C     ---------------                                                           
C     PLOT SPLINE FUNCTION OR DERIVATIVE OR INTEGRAL                            
C                                                                               
C     CALL ESPL(A,N,XA,XB,ND)      A(1)...A(N) = B-SPLINE COEFFICIENTS          
C                                  XA, XB = LOWER AND UPPER LIMIT               
C     ND = 0        FUNCTION                                                    
C        = 1, 2, 3  ND'TH DERIVATIVE                                            
C        = -1       INTEGRAL WITH LOWER LIMIT XA                                
C                                                                               
C     THE SUBPROGRAM FPL IS USED. THE USER HAS TO CALL FPL AFTER THE            
C     CALL OF ESPL.                                                             
C                                                                               
C     IF N IS LARGER THAN 200, THE USER HAS TO INCLUDE A COMMON                 
C     COMMON/SPCOM/I0,N0,S(N).                                                  
C                                                                               
      REAL A(N)                                                                 
                                                                                
      IF(ND.GT.3.OR.ND.LT.(-1)) GOTO 100                                        
      K=FLOAT(120-N)/FLOAT(N-1)+0.5                                             
      IF(K.LT.1) K=1                                                            
      M=N+(N-1)*K                                                               
C     WRITE(*,*) 'ESPLK M N K ',M,N,K                                           
      DO I=1,M                                                                  
       X=(FLOAT(M-I)*XA+FLOAT(I-1)*XB)/FLOAT(M-1)                               
       IF(ND.EQ.0) THEN                                                         
          F=ESPF(X,A,N,XA,XB)                                                   
       ELSE IF(ND.GT.0) THEN                                                    
          F=ESPD(X,A,N,XA,XB,ND)                                                
       ELSE                                                                     
          F=ESPI(XA,X,A,N,XA,XB)                                                
       END IF                                                                   
       CALL SPCURV(IFG,INR,X,F)                                                 
      END DO                                                                    
      IF(ND.NE.0) GOTO 100                                                      
C     M=N-2                                                                     
C     DO I=1,M                                                                  
C      X=(FLOAT(M-I)*XA+FLOAT(I-1)*XB)/FLOAT(M-1)                               
C      F=ESPF(X,A,N,XA,XB)                                                      
C      CALL SPDATA(IFG,INR+1,X,F,0.0,0.0)                                       
C     END DO                                                                    
  100 RETURN                                                                    
      END                                                                       
                                                                                
      SUBROUTINE TSLOOK                                                         
*                                                                               
*     LOOK test histograms                                                      
*                                                                               
*     common for TEST histograms ---------------------------------------        
      PARAMETER (NBINS=80,NDIM=10+NBINS)                                        
      COMMON/CRTEST/FTI(101),                                                   
     + DATAXT(NDIM),DATAXA(NDIM),DATAYT(NDIM),DATAYA(NDIM),DATAYB(NDIM),        
     + RANDXT(NDIM),RANDXA(NDIM),RANDYT(NDIM),RANDYA(NDIM),BACKYT(NDIM),        
     + NDATA,NMOCA,NDACK,NBACK                                                  
*     end of common for TEST histograms --------------------------------        
     +                                                                          
*     ...                                                                       
*     data histograms X                                                         
      CALL AHISTL(901,1,DATAXT)                                                 
      CALL AHISTL(901,2,DATAXA)                                                 
      CALL  STEXT(901,4,'True data X, total and accepted')                      
      CALL  STEXT(901,1,'True X')                                               
*     data histograms Y                                                         
      CALL AHISTL(902,1,DATAYT)                                                 
      CALL AHISTL(902,2,DATAYA)                                                 
      CALL AHISTL(902,4,DATAYB)                                                 
      CALL  STEXT(902,4,'Data Y, total, accepted, background')                  
      CALL  STEXT(902,1,'Y(1) = measured')                                      
                                                                                
*     moca histograms X                                                         
      CALL AHISTL(903,1,RANDXT)                                                 
      CALL AHISTL(903,2,RANDXA)                                                 
      CALL  STEXT(903,4,'True moca X, total and accepted')                      
      CALL  STEXT(902,1,'True X')                                               
*     moca histograms Y                                                         
      CALL AHISTL(904,1,RANDYT)                                                 
      CALL AHISTL(904,2,RANDYA)                                                 
      CALL  STEXT(904,4,'True moca Y, total and accepted')                      
      CALL  STEXT(904,1,'Y(1) = measured')                                      
                                                                                
*     back histograms Y                                                         
      CALL AHISTL(905,1,BACKYT)                                                 
      CALL  STEXT(905,4,'Background total')                                     
      CALL  STEXT(905,1,'Y(1) = measured')                                      
*     plot histograms                                                           
                                                                                
*     integrated luminosities for data ...                                      
      DLUMI=FLOAT(NDATA)        ! assuming total cross section of 1             
*                                                                               
*     plot true function as figure 901                                          
      XA=0.0                                                                    
      XB=2.0                                                                    
      FACTOR=DLUMI*2.0/FLOAT(NBINS)                                             
      DO I=1,101                                                                
*      x value                                                                  
       X=0.01*(XA*FLOAT(101-I)+XB*FLOAT(I-1))                                   
*      function plot                                                            
       CALL SPCURV(901,0,X,FACTOR*FTI(I))                                       
      END DO                                                                    
*                                                                               
*     plot true function as figure 101                                          
      XA=0.0                                                                    
      XB=2.0                                                                    
      FACTOR=DLUMI*2.0/FLOAT(NBINS)                                             
      DO I=1,101                                                                
*      x value                                                                  
       X=0.01*(XA*FLOAT(101-I)+XB*FLOAT(I-1))                                   
*      function plot of normalized function                                     
       CALL SPCURV(101,0,X,FTI(I))                                              
      END DO                                                                    
                                                                                
  100 RETURN                                                                    
      END                                                                       
                                                                                
                                                                                
                                                                                
      SUBROUTINE AHISTG(IFG,ISTYLE,A)                                           
*     histogram to graphic (LOOK) as data points                                
      REAL A(*)                                                                 
      NB=A(4)                                                                   
      IF(A(5).EQ.0.0.OR.NB.EQ.0) RETURN                                         
      DX=0.5*(A(3)-A(2))/FLOAT(NB)                                              
      DO I=1,NB                                                                 
       X=A(2)+FLOAT(2*I)*DX-DX                                                  
       Y=A(10+I)                                                                
       DY=0.0                                                                   
       IF(A(5).EQ.A(8)) DY=SQRT(Y)                                              
       CALL SPDATA(IFG,ISTYLE,X,Y,DX,DY)                                        
      END DO                                                                    
      END                                                                       
                                                                                
      SUBROUTINE AHISTL(IFG,ISTYLE,A)                                           
*     histogram to graphic (LOOK) as data points                                
      REAL A(*)                                                                 
      NB=A(4)                                                                   
      IF(A(5).EQ.0.0.OR.NB.EQ.0) RETURN                                         
      CALL SPHIST(IFG,ISTYLE,A(11),NB,A(2),A(3))                                
*        ISTYLE.EQ.0  NORMAL      WHITE                                         
*        ISTYLE.EQ.1  LHAT   HIST RED                                           
*        ISTYLE.EQ.2  RHAT   HIST GREEN                                         
*        ISTYLE.EQ.3  SOLID  HIST BLUE                                          
      END                                                                       
