*************************************************************************
* 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
do i=1,nclstr1
cl_used(i)=0 !init mask of clusters associated to a track
enddo
dedx_x_min = dedx_x_min_mip
dedx_y_min = dedx_y_min_mip
*------------------------------------------------------
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
* ////////////////////////////////////////////////
* 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
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
8800 continue
* ///////////////////////////////////////////////
* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
100 continue
return
end
************************************************************