--- DarthVader/TrackerLevel2/src/F77/analysisflight.f 2007/02/16 14:56:02 1.10 +++ DarthVader/TrackerLevel2/src/F77/analysisflight.f 2009/12/30 14:22:22 1.21 @@ -6,7 +6,7 @@ * - create LEVEL2 ntuple * ************************************************************************* - subroutine analysisflight(fin) + subroutine analysisflight include 'commontracker.f' include 'level1.f' @@ -19,9 +19,7 @@ * input flag * -* 0 - PFA = ETA -* 1 - PFA = COG4 - integer fin +c integer fin * flag to chose PFA character*10 PFA @@ -33,6 +31,11 @@ external npl external acoordsi,coordsi,nld,coord,dcoord + integer hitplanex(6) + integer hitplaney(6) + + logical BAD_NUCLEUS + ************************************************************ ************************************************************ ************************************************************ @@ -42,48 +45,43 @@ ************************************************************ ************************************************************ ************************************************************ - PFA='ETA' - if(fin.eq.0)PFA='ETA' - if(fin.eq.2)PFA='ETA2' - if(fin.eq.3)PFA='ETA3' - if(fin.eq.4)PFA='ETA4' - if(fin.eq.10)PFA='COG' - if(fin.eq.11)PFA='COG1' - if(fin.eq.12)PFA='COG2' - if(fin.eq.13)PFA='COG3' - if(fin.eq.14)PFA='COG4' - -c$$$ PFA='ETA' -c$$$ DEBUG=.true. - if(DEBUG)PRINT*,'P.F.A. --> ',fin,PFA + call idtoc(pfaid,PFA) + - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'----------------------------------' -c 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 -*------------------------------------------------------ - call init_level2 - call init_hough -*------------------------------------------------------ + 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 -*------------------------------------------------------ -c$$$ if(nclstr1.gt.nclstrmax_level2)then -c$$$ goto 8800 !fill nt-uple and go to next event -c$$$ 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 + + * /////////////////////////////////////////////// * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ @@ -93,16 +91,214 @@ 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. +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 * ********************************************************** @@ -111,8 +307,7 @@ call fill_level2_siglets - if(DEBUG)then - + if(DEBUG.EQ.1)then print*,'' print*,'DONE!' print*,'' @@ -127,17 +322,16 @@ 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 +c 8800 continue + continue * /////////////////////////////////////////////// * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ - 100 continue +c 100 continue + continue return end