************************************************************************* * 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 ************************************************************