--- DarthVader/TrackerLevel2/src/F77/analysisflight.f 2007/08/20 16:07:16 1.16 +++ DarthVader/TrackerLevel2/src/F77/analysisflight.f 2008/12/05 08:26:47 1.17 @@ -31,6 +31,11 @@ external npl external acoordsi,coordsi,nld,coord,dcoord + integer hitplanex(6) + integer hitplaney(6) + + logical BAD_NUCLEUS + ************************************************************ ************************************************************ ************************************************************ @@ -40,25 +45,9 @@ ************************************************************ ************************************************************ ************************************************************ -c$$$ TRACKMODE = 0 -c$$$ FACT = 100. -c$$$ ISTEPMIN = 3 call idtoc(pfaid,PFA) -c$$$ PFA='COG4' -c$$$ if(pfaid.eq.0)PFA='ETA' -c$$$ if(pfaid.eq.2)PFA='ETA2' -c$$$ if(pfaid.eq.3)PFA='ETA3' -c$$$ if(pfaid.eq.4)PFA='ETA4' -c$$$ if(pfaid.eq.10)PFA='COG' -c$$$ if(pfaid.eq.11)PFA='COG1' -c$$$ if(pfaid.eq.12)PFA='COG2' -c$$$ if(pfaid.eq.13)PFA='COG3' -c$$$ if(pfaid.eq.14)PFA='COG4' -*********************************************************** - -c if(DEBUG.EQ.1)PRINT*,'P.F.A. --> ',fin,PFA if(DEBUG.EQ.1)then print*,'----------------------------------' @@ -70,27 +59,25 @@ endif -*------------------------------------------------------ - call init_level2 - call init_hough -*------------------------------------------------------ - -*------------------------------------------------------ -* cut on maximum number of clusters -*------------------------------------------------------ -c$$$ if(nclstr1.gt.nclstrmax_level2)then -c$$$ goto 8800 !fill nt-uple and go to next event -c$$$ 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 - -c$$$ if(DEBUG.EQ.1)then -c$$$ print*,'----------------------------------' -c$$$ print*,iev,' ** ',nev2 -c$$$ endif + dedx_x_min = dedx_x_min_mip + dedx_y_min = dedx_y_min_mip +*------------------------------------------------------ + call init_level2 + call init_hough +*------------------------------------------------------ + 888 continue + + * /////////////////////////////////////////////// * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ @@ -100,15 +87,213 @@ goto 880 !fill ntp and go to next event endif - call fill_hough +* ------------------------------------- +* 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. + RECOVER_SINGLETS=.true. + ntrk_old=ntrk + + if(DEBUG.EQ.1) +cc if(.true.) + $ print*,'***** NEW SEARCH ****' + goto 888 !back to track finding + endif 880 continue @@ -119,7 +304,6 @@ call fill_level2_siglets if(DEBUG.EQ.1)then - print*,'' print*,'DONE!' print*,'' @@ -134,9 +318,6 @@ do iv = 1,nviews ngood = ngood + good1(iv) enddo -c$$$ if(ngood.ne.0)print*,'* WARNING * Event ' -c$$$ $ ,':LEVEL2 event status: ' -c$$$ $ ,(good2(i),i=1,nviews) 8800 continue