*************************************************************************
*     Program analysis.f
*
*     - reads cluster information (LEVEL1, reduction.f output ntuple)
*     - perform track identification and fit
*     - create LEVEL2 ntuple
*
*************************************************************************
      program momanhough

      include '../common/commontracker.f'
      include '../common/common_momanhough.f'
c      include '../common/common_level2debug.f'
      include '../common/common_mech.f'
      include '../common/common_xyzPAM.f'
      include '../common/common_mini_2.f'
      include '../common/calib.f'
      include '../common/level1.f'
      include '../common/level2.f'

      include '../common/momanhough_init.f'

*     flag to set debug mode
      logical DEBUG
      common/dbg/DEBUG
      logical DBG_FILLED
      data DBG_FILLED/.false./

*     flag to chose PFA
      character*10 PFA
      common/FINALPFA/PFA


c      parameter (inf=1.e8)      !just a huge number...
        
*     external functions
      external npl
      external acoordsi,coordsi,nld,coord,dcoord
c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------

      character*24 processing_date

      parameter (lun_data_level1=71) !data file id number
      parameter (lun_data_level2=72) !data file id number
      parameter (lun_data_calib=74) !data file id number

      character*74 data_file    !data file name
      character*74 data_dir    !data file name
      character*74 data_file_calib
      character*74 data_file_level1
      character*74 data_file_level2
# ifndef TEST2003
      parameter(ncalibmax=50)
      character*40 file_calib(ncalibmax)
      parameter(lun_calib_list=66) !calibration list file id
      integer which_calib_last
# endif

      integer minevent          !first event to be analysed
      logical FIMAGE            !


      COMMON/QUEST/IQUEST(100)  
c     !permette di ottenere ntuple funzionanti nonostante
c     ! il messaggio dei 64K di RZOUT...!???

c*****************************************************
cccccc 11/9/2005 modified by david fedele
c      swcode=202
cccccc 12/10/2005 modified by Elena Vannuccini
      swcode=300
c****************************************************
c------------------------------------------------------------------------
c     
c     HBOOK initialization
c     
c------------------------------------------------------------------------
     
      call HLIMIT(NWPAWC)


c------------------------------------------------------------------------
c     
c     reads input informations
c     
c------------------------------------------------------------------------
      call fdate(processing_date)
      write(*,101)
     $     processing_date
 101  format(/
     $      ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $      ,'*                                                 *',/
     $      ,'* ANALYSIS                                        *',/     
     $      ,'*                                                 *',/
     $      ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $     ,a24,/
     $     )

 111  format(a)
      print*,'Data directory:'      
      read(*,111)data_dir
      print*,data_dir
      print*,'File identifier: (DATE_NUM)'    
      read(*,*)data_file     
      print*,data_file
      minevent=1
      print*,'Maximum number of events to be analized:' !20000
      read(*,*) ntotev
      print*,ntotev    
      print*,'Position-finding algorythm: (GOG2,ETA2)'
      read(*,*)PFA
      PFA=PFA(1:lnblnk(PFA))
      print*,PFA
      print*,'DEBUG mode: (T, F)'
 11   format(l1)
      read(*,11)DEBUG
      print*,DEBUG
      print*,'---------------------------------------------------'


****** INITIALIZATIONS *************************************

*     1) read charge-correlation parameters
      print*,' '
      print*,'- read charge-correlation parameters'
      print*,' '
      call readchargeparam

*     1) read mip parameters
      print*,' '
      print*,'- read mip parameters'
      print*,' '
      call readmipparam

*     2) read z coordinates of the planes
      print*,' '
      print*,'- read z coordinates of the planes'
      print*,' '
      call mech_sensor          !reads sensors centres coordinates
      do ip=1,nplanes
        fitz(ip)=z_mech_sensor(ip,1,1)*0.1 !cm
*     gets planes mechanical z positions
*     (in mm) and sets them in micrometers
      enddo                     


*     3) read eta PFA parameters
      print*,' '
      print*,'- read P.F.A. parameters'
      print*,' '
      call readetaparam
      
      print*,' '
      print*,'- First guess P.F.A. >>>> ',PFAdef
      print*,' '
      print*,'- Final       P.F.A. >>>> ',PFA
      print*,' '

*     4) read magnetic field map      
      print*,' '
      print*,'- read magnetic field map'
      print*,' '
c      call read_B(5)            !legge il nome da STI!!.. temporaneo
c      call read_B_2maps(5)            !legge il nome da STI!!.. temporaneo
      call read_B

*     5) read allignment parameters
      print*,' '
      print*,'- read aligment parameters'
      print*,' '
      call readalignparam
 
************************************************************





c------------------------------------------------------------------------
c     
c     LEVEL 2 ntuple booking
c     
c------------------------------------------------------------------------


 503  format(a,'DW_',a,'_level2.rz')
      write(data_file_level2,503)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))

      print*,''
      print*,'------------------------------------'
      print*,'  Creating LEVEL2 rz file'
      print*,'  ',data_file_level2
      print*,'------------------------------------'


C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     largest RZ file: IQUEST(10) records x LREC words x 4 byte
C     with the following settings: 65000 x 4096 x 4 = 1G
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      IQUEST(10)=65000          
c     !permette di ottenere ntuple funzionanti nonostante
c     ! il messaggio dei 64K di RZOUT...!???
      call HROPEN(lun_data_level2,
     $     'LEVEL2',data_file_level2,'QNP',4096,istat) !opens rz  
      if(istat.ne.0) goto 19      
      call book_level2

      if(DEBUG)then
         print*,'(creating DEBUG nt-uple and histos)'
         call book_debug
      endif
      


c------------------------------------------------------------------------
c     
c     fills bad variables with online DSP bad strips from calib histograms
c     for bad strip exclusion
c     
c------------------------------------------------------------------------
        
# ifdef TEST2003
c------------------------------------------------------------------------
c     
c     opens calib file
c     
c------------------------------------------------------------------------

 505  format(a,'output_',a,'_calib.rz')
      write(data_file_calib,505)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
         
      print*,' '
      print*,'OPENING CALIB FILE...', data_file_calib
      
      IQUEST(10)=65000
      call HROPEN(lun_data_calib,'IN',data_file_calib,'QP',4096,istat)
      if(istat.ne.0) goto 17
      call HCDIR('//IN',' ')
      
      call HRIN(0,9999,0)       !puts histograms in memory
      
      print*,' '

      print*,' '
      print*,'READING PEDESTAL, SIGMA AND BADSTRIP HISTOGRAMS...'
      print*,' '
      
      call fillpedsig
      
      call HREND('IN')
      close (lun_data_calib)  
# else

      print*,' '
      print*,'OPENING CALIBRATION-LIST FILE:'
 501  format(a,'DW_',a,'_calib.txt')
      write(data_file_calib,501)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
      print*,data_file_calib
      open(lun_calib_list,
     $     FILE=data_file_calib(1:LNBLNK(data_file_calib))
     $     ,STATUS='UNKNOWN'
     $     ,IOSTAT=iostat
     $     )

 113  format(i5,'    ',a25)

      do i=1,ncalibmax

         read(lun_calib_list,113,IOSTAT=iostat) n_cal_list
     $        ,file_calib(i)(1:LNBLNK(file_calib(i)))
         if(iostat.ne.0)then
            ncal=i-1
            goto 2000
         endif
         print*,n_cal_list,' - '
     $        ,file_calib(i)(1:LNBLNK(file_calib(i)))

      enddo

 2000 close(lun_calib_list)
      which_calib_last=0

# endif 
      
c------------------------------------------------------------------------
c     
c     opens level1 file
c     
c------------------------------------------------------------------------

# ifdef TEST2003
 504  format(a,'output_',a,'_level1.rz')
# else
 504  format(a,'DW_',a,'_level1.rz')
# endif
      write(data_file_level1,504)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
      print*,''
      print*,'OPENING LEVEL1 FILE:'
      print*,data_file_level1
      IQUEST(10)=65000          
c     !permette di ottenere ntuple funzionanti nonostante
c     ! il messaggio dei 64K di RZOUT...!???
      call HROPEN(lun_data_level1,
     $     'LEVEL1',data_file_level1,'QP',4096,istat) !opens rz  
      if(istat.ne.0) goto 19      
      call HRIN(ntp_level1,9999,20)       
      call access_level1
c      call HPRNTU(ntp_level1+20)  
      call HNOENT(ntp_level1+20,iemax0)
      print*,'  events',iemax0


************************************************************
************************************************************
************************************************************
*
*     start track analysis
*
************************************************************
************************************************************
************************************************************
      maxevent=minevent+ntotev-1
      do iev = minevent,MIN(iemax0,maxevent) !loop on events

         call HCDIR('//LEVEL1',' ')
         call HGNT(ntp_level1+20,iev,ierr) !reads an event
         if(ierr.ne.0) goto 21
*------------------------------------------------------
*        LEVEL2 N-TUPLE INITIALIZATIONS
         call init_level2
         if(DEBUG)call init_level2_debug
         if(DEBUG)DBG_FILLED=.false.

         if(.not.good1)then
            goto 8800           !fill nt-uple and go to next event
         endif
*------------------------------------------------------


# ifndef TEST2003
      
         if(which_calib.ne.which_calib_last.and.
     $        which_calib.ne.0)then
            
            data_file_calib=
     $           data_dir(1:LNBLNK(data_dir))//
     $           file_calib(which_calib)
     $           (1:LNBLNK(file_calib(which_calib)))
c     print*,data_file_calib
            print*,''
            print*,
     $           '@ event ',nev2
     $           ,'   (CPU pkt N.',pkt_num1,')'
            print*,'--> ',data_file_calib            
            IQUEST(10)=65000          
            call HROPEN(lun_data_calib,
     $           'CALIB',data_file_calib,'QP',4096,istat) !opens
            if(istat.ne.0) goto 19      
            call HRIN(0,9999,0)       
            call fillpedsig
            do iview=1,nviews
               call HDELET(id_hi_bad+iview)
               call HDELET(id_hi_ped+iview)
               call HDELET(id_hi_sig+iview)
            enddo
            call HREND('CALIB')
            close(lun_data_calib)            
            which_calib_last=which_calib           
            
         elseif(which_calib.eq.0)then
            
            nocalib=nocalib+1
            good2=.false.
            goto 8800 !fill nt-uple and go to next event
            
         endif

# endif

      
*------------------------------------------------------
*     cut on maximum number of clusters
*------------------------------------------------------
         if(nclstr1.gt.nclstrmax_level2)then
            goto 8800           !fill nt-uple and go to next event
         endif

         do i=1,nclstr1
            cl_used(i)=0        !init mask of clusters associated to a track
         enddo

         if(DEBUG)then
            print*,'----------------------------------'
            print*,iev,'   ** ',nev2
         endif

*-------------------------------------------------------------------------------
*     STEP 1
*-------------------------------------------------------------------------------
*     X-Y cluster association
*     
*     Clusters are associated to form COUPLES
*     Clusters not associated in any couple are called SINGLETS
*
*     Track identification (Hough transform) and fitting is first done on couples. 
*     Hence singlets are possibly added to the track.
*    
*     Variables assigned by the routine "cl_to_couples" are those in the 
*     common blocks:
*     - common/clusters/cl_good
*     - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2
*     - common/singlets/ncls,cls,cl_single
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

         iflag=0
         call cl_to_couples(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*     selezione di tracce pulite per diagnostica
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c$$$         if(DEBUG)then
c$$$            do ip=1,nplanes 
c$$$               if(ncp_plane(ip).ne.1)good2=.false.
c$$$            enddo 
c$$$c            if(good2.eq.0)goto 100!next event 
c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event 
c$$$         endif
         
 

*-----------------------------------------------------
*-----------------------------------------------------
*     HOUGH TRASFORM
*-----------------------------------------------------
*-----------------------------------------------------


*-------------------------------------------------------------------------------
*     STEP 2
*-------------------------------------------------------------------------------
*     
*     Association of couples to form
*     - DOUBLETS in YZ view
*     - TRIPLETS in XZ view
*    
*     Variables assigned by the routine "cp_to_doubtrip" are those in the 
*     common blocks:
*     - common/hough_param/
*     $     alfayz1,  !Y0
*     $     alfayz2,  !tg theta-yz
*     $     alfaxz1,  !X0
*     $     alfaxz2,  !tg theta-xz
*     $     alfaxz3   !1/r
*     - common/doublets/ndblt,cpyz1,cpyz2
*     - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

         iflag=0
         call cp_to_doubtrip(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif


*-------------------------------------------------------------------------------
*     STEP 3
*-------------------------------------------------------------------------------
*     
*     Classification of doublets and triplets to form CLOUDS,  
*     according to distance in parameter space. 
*     
*     cloud = cluster of points (doublets/triplets) in parameter space
*
*     
*    
*     Variables assigned by the routine "doub_to_YZcloud" are those in the 
*     common blocks:
*     - common/clouds_yz/                   
*     $     nclouds_yz                        
*     $     ,alfayz1_av,alfayz2_av           
*     $     ,ptcloud_yz,db_cloud,cpcloud_yz
*
*     Variables assigned by the routine "trip_to_XZcloud" are those in the 
*     common blocks:
*      common/clouds_xz/                   
*     $      nclouds_xz     xz2_av,alfaxz3_av 
*     $     ,ptcloud_xz,tr_cloud,cpcloud_xz           
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

         iflag=0
         call doub_to_YZcloud(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif
         iflag=0
         call trip_to_XZcloud(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif

c*****************************************************
cccccc 01/12/2005 modified by elena
         if(DEBUG)then
            call fill_level2_clouds
            call HCDIR('//LEVEL2',' ')
            call HFNT(ntp_level2+1) !fill DEBUG nt-uple
            DBG_FILLED=.true.
         endif
c*****************************************************

*-------------------------------------------------------------------------------
*     STEP 4   (ITERATED until any other physical track isn't found)
*-------------------------------------------------------------------------------
*     
*     YZ and XZ clouds are combined in order to obtain the initial guess
*     of the candidate-track parameters.
*     A minimum number of matching couples between YZ and XZ clouds is required.
*
*     A TRACK CANDIDATE is defined by 
*     - the couples resulting from the INTERSECTION of the two clouds, and 
*     - the associated track parameters (evaluated by performing a zero-order  
*       track fitting) 
*
*     The NTRACKS candidate-track parameters are stored in common block:
*    
*     - common/track_candidates/NTRACKS,AL_STORE
*     $     ,XV_STORE,YV_STORE,ZV_STORE
*     $     ,XM_STORE,YM_STORE,ZM_STORE
*     $     ,RESX_STORE,RESY_STORE
*     $     ,AXV_STORE,AYV_STORE
*     $     ,XGOOD_STORE,YGOOD_STORE
*     $     ,CP_STORE,RCHI2_STORE
*
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
         ntrk=0                 !counter of identified physical tracks

11111    continue               !<<<<<<< come here when performing a new search

         iflag=0
         call clouds_to_ctrack(iflag)
         if(iflag.eq.1)then     !no candidate tracks found
            goto 880            !fill ntp and go to next event   
         endif

         FIMAGE=.false.         !processing best track (not track image)
         ibest=0                !best track among candidates
         iimage=0               !track image
*     ------------- select the best track -------------
         rchi2best=1000000000.
         do i=1,ntracks
            if(RCHI2_STORE(i).lt.rchi2best.and.
     $         RCHI2_STORE(i).gt.0)then
               ibest=i
               rchi2best=RCHI2_STORE(i)
            endif
         enddo
         if(ibest.eq.0)goto 880 !>> no good candidates 
*-------------------------------------------------------------------------------     
*     The best track candidate (ibest) is selected and a new fitting is performed. 
*     Previous to this, the track is refined by:
*     - possibly adding new COUPLES or SINGLETS from the missing planes
*     - evaluating the coordinates with improved PFAs 
*       ( angle-dependent ETA algorithms )
*-------------------------------------------------------------------------------

 1212    continue               !<<<<< come here to fit track-image

         if(.not.FIMAGE)then    !processing best candidate 
            icand=ibest            
         else                   !processing image
            icand=iimage
            iimage=0
         endif
         if(icand.eq.0)then
            print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
     $           ,ibest,iimage
            return
         endif

*     *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
         call refine_track(icand)
*     *-*-*-*-*-*-*-*-*-*-*-*-*-*-*

*     **********************************************************
*     ************************** FIT *** FIT *** FIT *** FIT ***
*     **********************************************************
         do i=1,5
            AL(i)=dble(AL_STORE(i,icand))
         enddo
         ifail=0                !error flag in chi2 computation
         jstep=0                !# minimization steps

         call mini_2(jstep,ifail)
         if(ifail.ne.0) then
            if(DEBUG)then
               print *,
     $              '*** MINIMIZATION FAILURE *** (mini_2) '
     $              ,iev
            endif
            chi2=-chi2 
         endif 
         
         if(DEBUG)then
            print*,'----------------------------- improved track coord'
22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
            do ip=1,6
               write(*,22222)ip,zm(ip),xm(ip),ym(ip)
     $              ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
     $              ,xgood(ip),ygood(ip),resx(ip),resy(ip)
            enddo
         endif

c         rchi2=chi2/dble(ndof)
         if(DEBUG)then
            print*,' ' 
            print*,'****** SELECTED TRACK *************'
            print*,'#         R. chi2        RIG'
            print*,' --- ',chi2,' --- '
     $           ,1./abs(AL(5)) 
            print*,'***********************************'
         endif
*     **********************************************************
*     ************************** FIT *** FIT *** FIT *** FIT ***
*     **********************************************************


*     ------------- search if the track has an IMAGE -------------
*     ------------- (also this is stored )           -------------
         if(FIMAGE)goto 122     !>>> jump! (this is already an image)
*     now search for track-image, by comparing couples IDs
         do i=1,ntracks
            iimage=i
            do ip=1,nplanes
               if(     CP_STORE(nplanes-ip+1,icand).ne.
     $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0
            enddo
            if(  iimage.ne.0.and.
c     $           RCHI2_STORE(i).le.CHI2MAX.and.
c     $           RCHI2_STORE(i).gt.0.and.
     $           .true.)then 
               if(DEBUG)print*,'Track candidate ',iimage
     $              ,' >>> TRACK IMAGE >>> of'
     $              ,ibest
               goto 122         !image track found
            endif
         enddo
 122     continue

*     --- and store the results --------------------------------
         ntrk = ntrk + 1                   !counter of found tracks
         if(.not.FIMAGE
     $        .and.iimage.eq.0) image(ntrk)= 0
         if(.not.FIMAGE
     $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
         if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous

         call fill_level2_tracks(ntrk)     !==> good2=.true.
c         print*,'++++++++++ iimage,fimage,ntrk,image '
c     $        ,iimage,fimage,ntrk,image(ntrk)

         if(ntrk.eq.NTRKMAX)then
            if(DEBUG)
     $           print*,
     $           '** warning ** number of identified '// 
     $           'tracks exceeds vector dimension '
     $           ,'( ',NTRKMAX,' )'
cc            good2=.false.
            goto 880            !fill ntp and go to next event
         endif
         if(iimage.ne.0)then
            FIMAGE=.true.       !
            goto 1212           !>>> fit image-track
         endif

*     --- then remove selected clusters (ibest+iimage) from clouds ----
         call clean_XYclouds(ibest)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif

*     **********************************************************
*     condition to start a new search
*     **********************************************************
         ixznew=0
         do ixz=1,nclouds_xz
            if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1
         enddo
         iyznew=0
         do iyz=1,nclouds_yz
            if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1
         enddo
         
         if(ixznew.ne.0.and.
     $      iyznew.ne.0.and.
     $        rchi2best.le.CHI2MAX.and.
c     $        rchi2best.lt.15..and.
     $        .true.)then
            if(DEBUG)then
               print*,'***** NEW SEARCH ****'
            endif
            goto 11111          !try new search
            
         endif
*     **********************************************

*     >>>>>>>>>>>>>>>>>>> NT-UPLE filling <<<<<<<<<<<<<<<<<<<<
*       +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +
*      / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ /
*     +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   

 880     continue

*     **********************************************************
*     stores info about clusters not associated with any track
*     **********************************************************

c*****************************************************
cccccc 27/9/2005 modified by david fedele
*     count #cluster per plane not associated to any track
c$$$         do icl=1,nclstr1
c$$$            if(cl_used(icl).eq.0)then !cluster not included in any track
c$$$               ip=nplanes-npl(VIEW(icl))+1
c$$$               if(mod(VIEW(icl),2).eq.0)nclsx(ip)=nclsx(ip)+1
c$$$               if(mod(VIEW(icl),2).eq.1)nclsy(ip)=nclsy(ip)+1
c$$$            endif
c$$$c            print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
c$$$         enddo
c**********************************************************************

         call fill_level2_siglets

c         print*,'****** ',iev,' - ',ntrk,nclsx,nclsy
c         print*,'****** ',iev,' - ',ntrk,' --- ',(BdL(i),i=1,ntrk)

         if(DEBUG)then

            print*,''
            print*,'DONE!'
            print*,''
            print*,'* summary *'
            print*,'tracks    ',ntrk
            print*,'cl used   ',(cl_used(i),i=1,nclstr1)
c*****************************************************
c$$$cccccc 27/9/2005 modified by david fedele
c$$$            print*,'cl unused (x-y)'
c$$$            do ip=1,nplanes
c$$$               print*,ip,' << ',nclsx(ip),nclsy(ip)
c$$$            enddo
c******************************************************
            print*,''
            print*,''
         endif

 8800    continue

         call HCDIR('//LEVEL2',' ')
         call HFNT(ntp_level2)  !fill LEVEL2 nt-uple         
         if(.not.DBG_FILLED.and.DEBUG)
     $        print*,'@@@@@@@@@@@@',ntp_level2+1,nev2_nt
         if(.not.DBG_FILLED.and.DEBUG)call HFNT(ntp_level2+1)

 100     continue
      enddo                     !end loop on events

c------------------------------------------------------------------------
c     
c     no error exit
c     
c------------------------------------------------------------------------

c$$$      print*,' '
c$$$      print*,'REDUCTION SUCCESSFULLY COMPLETED'
c$$$      print*,' '
c$$$      print*,' '
      
      goto 9000                 !happy ending


c------------------------------------------------------------------------
c     
c     data file opening error
c     
c------------------------------------------------------------------------
 19   continue
      
      print*,' '
      print*,'ERROR OPENING DATA FILE: ',data_file
      print*,' '
      print*,' '
      
      goto 9000                 !the end

c------------------------------------------------------------------------
c     
c     level1 ntuple event reading error
c     
c------------------------------------------------------------------------
      
 21   continue
      
      print*,' '
      print*,'ERROR WHILE READING LEVEL1 NTUPLE, AT EVENT
     $     : ',iev
      print*,' '
      print*,' '
      
      goto 9000                 !the end

c------------------------------------------------------------------------
c     
c     calib file opening error
c     
c------------------------------------------------------------------------
      
 17   continue
      
      print*,' '
      print*,'preanalysis: ERROR OPENING INPUT FILE: ',data_file_calib
      print*,' '
      print*,' '
      
      goto 9000                 !the end

c------------------------------------------------------------------------
c     
c     closes files and exits
c     
c------------------------------------------------------------------------

 9000 continue
      
      if(DEBUG)call HPRNTU(ntp_level2+1)
      call HPRNTU(ntp_level2)
      print*,''
      call HCDIR('//LEVEL2',' ')
      call HROUT(ntp_level2,ICYCLE,'T')
      print *,'- Stored LEVEL2 nt-uple (',nev2,' entries )'
      if(DEBUG)call HROUT(ntp_level2+1,ICYCLE,'T')
      if(DEBUG)print *,'- Stored DEBUG nt-uple '
      call HREND('level2')
      close(lun_data_level2)

      call HCDIR('//LEVEL1',' ')
      call HREND('level1')
      close(lun_data_level1)

      stop
      end


************************************************************


# include "momanhough-subroutines.F"


