--- DarthVader/TrackerLevel2/src/F77/analysisflight.f 2006/05/30 16:30:37 1.2 +++ DarthVader/TrackerLevel2/src/F77/analysisflight.f 2009/12/30 14:22:22 1.21 @@ -9,14 +9,18 @@ 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 'level1.f' include 'level2.f' +* input flag +* +c integer fin + * flag to chose PFA character*10 PFA common/FINALPFA/PFA @@ -27,6 +31,11 @@ external npl external acoordsi,coordsi,nld,coord,dcoord + integer hitplanex(6) + integer hitplaney(6) + + logical BAD_NUCLEUS + ************************************************************ ************************************************************ ************************************************************ @@ -36,38 +45,43 @@ ************************************************************ ************************************************************ ************************************************************ - PFA='ETA' - if(DEBUG)then + call idtoc(pfaid,PFA) + + + if(DEBUG.EQ.1)then print*,'----------------------------------' - print*,'START... ',good1,nclstr1,nclstrmax_level2 + print*,'Settings: ' + print*,'PFA ',pfaid,' ---> ',PFA + print*,'tracking mode ',trackmode + print*,'fit-tolerance factor ',fact + print*,'minimum n.step ',istepmin + endif -*------------------------------------------------------ -* LEVEL2 N-TUPLE INITIALIZATIONS - call init_level2 -c if(.not.good1)then - if(good1.eq.0)then - goto 8800 !fill nt-uple and go to next event - 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 ! + + 887 continue + *------------------------------------------------------ -* 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 + call init_level2 + call init_hough +*------------------------------------------------------ + + 888 continue + + * /////////////////////////////////////////////// * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ @@ -77,14 +91,214 @@ 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) - if(iflag.eq.1)then !bad event - goto 880 !fill ntp and go to next event - endif + +* //////////////////////////////////////////////// +* 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 * ********************************************************** @@ -93,8 +307,7 @@ call fill_level2_siglets - if(DEBUG)then - + if(DEBUG.EQ.1)then print*,'' print*,'DONE!' print*,'' @@ -104,14 +317,21 @@ print*,'' print*,'' endif + + ngood = 0 + do iv = 1,nviews + ngood = ngood + good1(iv) + enddo - 8800 continue +c 8800 continue + continue * /////////////////////////////////////////////// * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ - 100 continue +c 100 continue + continue return end