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

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'common_mech.f'
      include 'common_xyzPAM.f'
      include 'common_mini_2.f'
      include 'calib.f'
      include 'level2.f'

*     input flag
*
c      integer fin

*     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

      integer hitplanex(6)
      integer hitplaney(6)

      logical BAD_NUCLEUS

************************************************************
************************************************************
************************************************************
*
*     track analysis
*
************************************************************
************************************************************
************************************************************

      call idtoc(pfaid,PFA)
      

      if(DEBUG.EQ.1)then
         print*,'----------------------------------'
         print*,'Settings: '
         print*,'PFA                  ',pfaid,' ---> ',PFA
         print*,'tracking mode        ',trackmode
         print*,'fit-tolerance factor ',fact
         print*,'minimum n.step       ',istepmin
         
      endif

      RECOVER_SINGLETS = .false.
      RECOVER_NUCLEI   = .false.
      SECOND_SEARCH    = .false.
      
      ntrk_old = 0

      dedx_x_min = dedx_x_min_mip !mip search
      dedx_y_min = dedx_y_min_mip !

c 887  continue ! Emiliano, fix compilation warning (gcc 4.1), Label 887 defined but not used

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

      call init_level2
      call init_hough      
*------------------------------------------------------

 888  continue

     
*     ///////////////////////////////////////////////
*     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

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

*     ----------------------------------------------------
*     fill hough output variables only for the first searc
*     ----------------------------------------------------
      if(.not.SECOND_SEARCH)call fill_hough

      iflag=0
      call track_fitting(iflag)
      
     
*     ////////////////////////////////////////////////
*     PATCH to recover events with less than 3 couples
*     ////////////////////////////////////////////////     
*     if no tracks have been found 
*     and n.planes with at least one couple (planehit) 
*     is less than 3 try to recover applying the hough 
*     transform to singlet
      nplanehit=0
      do ip=1,NPLANES
         if(ncp_plane(ip).gt.0)nplanehit=nplanehit+1
      enddo
      if(DEBUG.eq.1)then
cc      if(.true.)then
         print*,''
         print*,'-----------------------------------'
         print*,'n.tracks     ',ntrk
         print*,'n.hit-planes ',nplanehit
      endif
      if(
     $     .not.RECOVER_SINGLETS.and. !recover only once
     $     ntrk.eq.0            .and. !no tracks has been found
     $     nplanehit.le.3       .and. !if less than 3 hit planes (n.real couples)
     $     .true.)then
         RECOVER_SINGLETS=.true.
cc         if(.true.)
         if(DEBUG.EQ.1)
     $        print*,' -----> ((( TRY TO RECOVER SINGLETS )))'
         goto 888               !back to track finding
      endif


*     @~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@
*
*     EV may 2014 ==> SKIP NUCLEI PATCH
*     (done differently...)
*
*     @~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@
      goto 7575


*     ////////////////////////////////////////////////
*     PATCH to recover nuclei
*     ////////////////////////////////////////////////
      itrk = ntrk               !check the last track first
      isimage = 0
      BAD_NUCLEUS=.false.
 150  continue

      dedxav  = 0.
      npt     = 0
      dedxmin = 1000000
      if(itrk.gt.0)then
         do ip=1,NPLANES
cc            print*,ip,' ** ',dedx_x(ip,itrk),dedx_y(ip,itrk)
            if(xgood_nt(ip,itrk).eq.1)then
               dedxav = dedxav + abs(dedx_x(ip,itrk))
               if(abs(dedx_x(ip,itrk)).lt.dedxmin)
     $              dedxmin=abs(dedx_x(ip,itrk))
               npt = npt + 1
            endif         
            if(ygood_nt(ip,itrk).eq.1)then 
               dedxav = dedxav + abs(dedx_y(ip,itrk))
               if(abs(dedx_y(ip,itrk)).lt.dedxmin)
     $              dedxmin=abs(dedx_y(ip,itrk))
               npt = npt + 1
            endif            
         enddo
         if(npt.gt.0)dedxav = dedxav/npt
         if(  dedxav.gt.dedx_min_nuclei )then
            if( (dedxav-dedxmin)/dedxav.gt.ddedx_min_nuclei )
     $           BAD_NUCLEUS=.true.
         endif
      endif
      if(DEBUG.eq.1)then
cc      if(.true.)then
         print*,''
         print*,'-----------------------------------'
         print*,'itrk          ',itrk
         print*,'dedx(average) ',dedxav,' cut ',dedx_min_nuclei
         print*,'dedx(minimum) ',dedxmin
         if(dedxav.gt.0)print*,'(av-min)/av',(dedxav-dedxmin)/dedxav
     $        ,' cut ',ddedx_min_nuclei
         print*,'BAD_NUCLEUS ? ',BAD_NUCLEUS
      endif

      itrk=0
      if(ntrk.gt.0)itrk = image(ntrk)              
      if(itrk.gt.0.and.isimage.eq.0)then      !hence check the image
         isimage=1         
         goto 150
      endif


      if(
     $     .not.RECOVER_NUCLEI.and. !recover only once
     $     (ntrk.eq.0.or.BAD_NUCLEUS).and. !no tracks found or bad nucleus
     $     .true.)then

         RECOVER_NUCLEI   = .true.
         dedx_x_min = dedx_x_min_nuclei
         dedx_y_min = dedx_y_min_nuclei
cc         if(.true.)
         if(DEBUG.EQ.1)
     $     print*,' -----> ((( TRY TO RECOVER NUCLEI )))'
*     ----------------------------------------------
*     in case of re-tracking, un-tag used clusters
*     ----------------------------------------------
         if(ntrk.ne.0)then
            do ip=1,NPLANES
               if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
               if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
            enddo
            ntrk = ntrk-1       !decrement track
            if(isimage.eq.1)then
               do ip=1,NPLANES
                  if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
                  if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
               enddo
               ntrk = ntrk-1    !decrement track image
            endif
         endif
         goto 888               !back to track finding
         
      endif
*     ----------------------------------------------
*     back to default parameters for the second-track search
*     ----------------------------------------------
      if(RECOVER_NUCLEI)then
         dedx_x_min = dedx_x_min_mip
         dedx_y_min = dedx_y_min_mip
         RECOVER_NUCLEI   = .false.     
      endif

*     @~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@
 7575 continue
*     @~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@

      if(ntrk.eq.0)goto 880

*     /////////////////////////////////////////////////////
*     PATCH to perform a more efficient second-track search
*     ///////////////////////////////////////////////////// 
*     -------------------------------------------------
*     n. hit planes and number of clusters
*     -------------------------------------------------
      do ip=1,NPLANES
         hitplanex(ip)=0
         hitplaney(ip)=0
         nclusters=0
      enddo
      do icl=1,nclstr1
         ip = npl(VIEW(icl))
         if(cl_good(icl).eq.0.or.cl_used(icl).eq.1)goto 9
         if(mod(VIEW(icl),2).eq.1)hitplaney(ip)=hitplaney(ip)+1 !YY         
         if(mod(VIEW(icl),2).eq.0)hitplanex(ip)=hitplanex(ip)+1 !XX         
         nclusters=nclusters+1
 9       continue
      enddo
      nhitplaney=0
      nhitplanex=0
      do ip=1,NPLANES
         if(hitplaney(ip).gt.0)nhitplaney=nhitplaney+1
         if(hitplanex(ip).gt.0)nhitplanex=nhitplanex+1
      enddo      
*     -------------------------------------------------
*     chi2 cut
*     -------------------------------------------------
      p0 = 1.111588e+00
      p1 = 1.707656e+00
      p2 = 1.489693e-01
      chi2m025 = p0 
     $     + abs(al_nt(5,ntrk))*p1 
     $     + al_nt(5,ntrk)*al_nt(5,ntrk)*p2
      chi2m025_0 = p0 + 0.07*p1 + 0.07*0.07*p2
cc      chi2m = (chi2m025-chi2m025_0+12**0.25)**4. !95%
      chi2m = (chi2m025-chi2m025_0+20.**0.25)**4. !???
*     -------------------------------------------------
*     image track
*     -------------------------------------------------
      ntrk_best = ntrk
      chi2_best = chi2_nt(ntrk)
      if(image(ntrk).ne.0)
     $     chi2_best = min(chi2_nt(ntrk),chi2_nt(image(ntrk)))
      if(DEBUG.eq.1)then
cc      if(.true.)then
         print*,''
         print*,'-----------------------------------'
         print*,'chi2_nt(',ntrk_best,') = ',chi2_best,' < ',chi2m,' ?'
         print*,'nhitplanex',nhitplanex
         print*,'nhitplaney',nhitplaney
         print*,'nclusters',nclusters
      endif
      if(
     $     ntrk.gt.ntrk_old.and. !one track must have been already found
     $     chi2_best.gt.0 .and.     !this track must not be garbage
     $     chi2_best.lt.chi2m .and. !"
     $     nhitplaney.ge.3.and. !a minimum number of hit planes is required
     $     nhitplanex.ge.3.and. !"
     $     nclusters.le.20.and. !limit on number of clusters to save time
     $     .true.)then

         SECOND_SEARCH=.true.
c         RECOVER_SINGLETS=.true.      !OPTIMIZATION  
         ntrk_old=ntrk

         if(DEBUG.EQ.1)
cc         if(.true.)
     $        print*,'***** NEW SEARCH ****'
         goto 888               !back to track finding

      endif

 880  continue

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

      call fill_level2_siglets
      
      if(DEBUG.EQ.1)then
         print*,''
         print*,'DONE!'
         print*,''
         print*,'* summary *'
         print*,'tracks    ',ntrk
         print*,'cl used   ',(cl_used(i),i=1,nclstr1)
         print*,''
         print*,''
      endif

      ngood = 0
      do iv = 1,nviews
         ngood = ngood + good1(iv)
      enddo
      
c 8800 continue
      continue
      
      
*     ///////////////////////////////////////////////
*     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

c 100  continue
      continue
      
      return
      end
      

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


