--- tracker/ground/source/analysis/momanhough-subroutines.F 2006/03/08 15:00:39 1.1.1.1 +++ tracker/ground/source/analysis/momanhough-subroutines.F 2006/03/20 19:43:33 1.2 @@ -1,173 +1,538 @@ ************************************************************ +* The following subroutines +* - track_finding >> hough transform +* - track_fitting >> bob golden fitting +* all the procedures to create LEVEL2 data, starting from LEVEL1 data. +* +* +* +* (This subroutine and all the dependent subroutines +* will be included in the flight software) +************************************************************ + subroutine track_finding(iflag) - subroutine readmipparam - include '../common/commontracker.f' + include '../common/common_momanhough.f' + include '../common/common_mech.f' + include '../common/common_xyzPAM.f' + include '../common/common_mini_2.f' include '../common/calib.f' + include '../common/level1.f' + include '../common/level2.f' - character*60 fname_param - 201 format('trk-LADDER',i1,'-mip.dat') - do ilad=1,nladders_view - write(fname_param,201)ilad - print *,'Opening file: ',fname_param - open(10, - $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) - $ ,STATUS='UNKNOWN' - $ ,IOSTAT=iostat - $ ) - if(iostat.ne.0)then - print*,'READMIPPARAM: *** Error in opening file ***' - return - endif - do iv=1,nviews - read(10,* - $ ,IOSTAT=iostat - $ )pip, - $ mip(int(pip),ilad) -c print*,ilad,iv,pip,mip(int(pip),ilad) - enddo - close(10) - enddo + include '../common/momanhough_init.f' + + logical DEBUG + common/dbg/DEBUG - return - end -*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** - subroutine readchargeparam +*------------------------------------------------------------------------------- +* STEP 1 +*------------------------------------------------------------------------------- +* X-Y cluster association +* +* Clusters are associated to form COUPLES +* Clusters not associated in any couple are called SINGLETS +* +* Track identification (Hough transform) and fitting is first done on couples. +* Hence singlets are possibly added to the track. +* +* Variables assigned by the routine "cl_to_couples" are those in the +* common blocks: +* - common/clusters/cl_good +* - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2 +* - common/singlets/ncls,cls,cl_single +*------------------------------------------------------------------------------- +*------------------------------------------------------------------------------- + +c iflag=0 + call cl_to_couples(iflag) + if(iflag.eq.1)then !bad event + goto 880 !fill ntp and go to next event + endif + +*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +* selezione di tracce pulite per diagnostica +*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c$$$ if(DEBUG)then +c$$$ do ip=1,nplanes +c$$$ if(ncp_plane(ip).ne.1)good2=.false. +c$$$ enddo +c$$$c if(good2.eq.0)goto 100!next event +c$$$c if(good2.eq.0)goto 880!fill ntp and go to next event +c$$$ endif + + + +*----------------------------------------------------- +*----------------------------------------------------- +* HOUGH TRASFORM +*----------------------------------------------------- +*----------------------------------------------------- + + +*------------------------------------------------------------------------------- +* STEP 2 +*------------------------------------------------------------------------------- +* +* Association of couples to form +* - DOUBLETS in YZ view +* - TRIPLETS in XZ view +* +* Variables assigned by the routine "cp_to_doubtrip" are those in the +* common blocks: +* - common/hough_param/ +* $ alfayz1, !Y0 +* $ alfayz2, !tg theta-yz +* $ alfaxz1, !X0 +* $ alfaxz2, !tg theta-xz +* $ alfaxz3 !1/r +* - common/doublets/ndblt,cpyz1,cpyz2 +* - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3 +*------------------------------------------------------------------------------- +*------------------------------------------------------------------------------- + +c iflag=0 + call cp_to_doubtrip(iflag) + if(iflag.eq.1)then !bad event + goto 880 !fill ntp and go to next event + endif + + +*------------------------------------------------------------------------------- +* STEP 3 +*------------------------------------------------------------------------------- +* +* Classification of doublets and triplets to form CLOUDS, +* according to distance in parameter space. +* +* cloud = cluster of points (doublets/triplets) in parameter space +* +* +* +* Variables assigned by the routine "doub_to_YZcloud" are those in the +* common blocks: +* - common/clouds_yz/ +* $ nclouds_yz +* $ ,alfayz1_av,alfayz2_av +* $ ,ptcloud_yz,db_cloud,cpcloud_yz +* +* Variables assigned by the routine "trip_to_XZcloud" are those in the +* common blocks: +* common/clouds_xz/ +* $ nclouds_xz xz2_av,alfaxz3_av +* $ ,ptcloud_xz,tr_cloud,cpcloud_xz +*------------------------------------------------------------------------------- +*------------------------------------------------------------------------------- + +c iflag=0 + call doub_to_YZcloud(iflag) + if(iflag.eq.1)then !bad event + goto 880 !fill ntp and go to next event + endif +c iflag=0 + call trip_to_XZcloud(iflag) + if(iflag.eq.1)then !bad event + goto 880 !fill ntp and go to next event + endif + 880 return + end + +************************************************************ + + subroutine track_fitting(iflag) + include '../common/commontracker.f' + include '../common/common_momanhough.f' + include '../common/common_mech.f' + include '../common/common_xyzPAM.f' + include '../common/common_mini_2.f' include '../common/calib.f' + include '../common/level1.f' + include '../common/level2.f' + + include '../common/momanhough_init.f' + + logical DEBUG + common/dbg/DEBUG + + logical FIMAGE ! - character*60 fname_param - 201 format('charge-l',i1,'.dat') - do ilad=1,nladders_view - write(fname_param,201)ilad - print *,'Opening file: ',fname_param - open(10, - $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) - $ ,STATUS='UNKNOWN' - $ ,IOSTAT=iostat - $ ) - if(iostat.ne.0)then - print*,'READCHARGEPARAM: *** Error in opening file ***' +*------------------------------------------------------------------------------- +* STEP 4 (ITERATED until any other physical track isn't found) +*------------------------------------------------------------------------------- +* +* YZ and XZ clouds are combined in order to obtain the initial guess +* of the candidate-track parameters. +* A minimum number of matching couples between YZ and XZ clouds is required. +* +* A TRACK CANDIDATE is defined by +* - the couples resulting from the INTERSECTION of the two clouds, and +* - the associated track parameters (evaluated by performing a zero-order +* track fitting) +* +* The NTRACKS candidate-track parameters are stored in common block: +* +* - common/track_candidates/NTRACKS,AL_STORE +* $ ,XV_STORE,YV_STORE,ZV_STORE +* $ ,XM_STORE,YM_STORE,ZM_STORE +* $ ,RESX_STORE,RESY_STORE +* $ ,AXV_STORE,AYV_STORE +* $ ,XGOOD_STORE,YGOOD_STORE +* $ ,CP_STORE,RCHI2_STORE +* +*------------------------------------------------------------------------------- +*------------------------------------------------------------------------------- + ntrk=0 !counter of identified physical tracks + +11111 continue !<<<<<<< come here when performing a new search + +c iflag=0 + call clouds_to_ctrack(iflag) + if(iflag.eq.1)then !no candidate tracks found + goto 880 !fill ntp and go to next event + endif + + FIMAGE=.false. !processing best track (not track image) + ibest=0 !best track among candidates + iimage=0 !track image +* ------------- select the best track ------------- + rchi2best=1000000000. + do i=1,ntracks + if(RCHI2_STORE(i).lt.rchi2best.and. + $ RCHI2_STORE(i).gt.0)then + ibest=i + rchi2best=RCHI2_STORE(i) + endif + enddo + if(ibest.eq.0)goto 880 !>> no good candidates +*------------------------------------------------------------------------------- +* The best track candidate (ibest) is selected and a new fitting is performed. +* Previous to this, the track is refined by: +* - possibly adding new COUPLES or SINGLETS from the missing planes +* - evaluating the coordinates with improved PFAs +* ( angle-dependent ETA algorithms ) +*------------------------------------------------------------------------------- + + 1212 continue !<<<<< come here to fit track-image + + if(.not.FIMAGE)then !processing best candidate + icand=ibest + else !processing image + icand=iimage + iimage=0 + endif + if(icand.eq.0)then + print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand + $ ,ibest,iimage return endif - do ip=1,nplanes - read(10,* - $ ,IOSTAT=iostat - $ )pip, - $ kch(ip,ilad),cch(ip,ilad),sch(ip,ilad) -c print*,ilad,ip,pip,kch(ip,ilad), -c $ cch(ip,ilad),sch(ip,ilad) - enddo - close(10) - enddo - return - end -*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** - subroutine readetaparam -* ----------------------------------------- -* read eta2,3,4 calibration parameters -* and fill variables: -* -* eta2(netabin,nladders_view,nviews) -* eta3(2*netabin,nladders_view,nviews) -* eta4(2*netabin,nladders_view,nviews) -* - include '../common/commontracker.f' - include '../common/calib.f' +* *-*-*-*-*-*-*-*-*-*-*-*-*-*-* + call refine_track(icand) +* *-*-*-*-*-*-*-*-*-*-*-*-*-*-* - character*40 fname_binning - character*40 fname_param -c character*120 cmd1 -c character*120 cmd2 +* ********************************************************** +* ************************** FIT *** FIT *** FIT *** FIT *** +* ********************************************************** + do i=1,5 + AL(i)=dble(AL_STORE(i,icand)) + enddo + ifail=0 !error flag in chi2 computation + jstep=0 !# minimization steps + call mini_2(jstep,ifail) + if(ifail.ne.0) then + if(DEBUG)then + print *, + $ '*** MINIMIZATION FAILURE *** (mini_2) ' + $ ,iev + endif + chi2=-chi2 + endif + + if(DEBUG)then + print*,'----------------------------- improved track coord' +22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5) + do ip=1,6 + write(*,22222)ip,zm(ip),xm(ip),ym(ip) + $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip) + $ ,xgood(ip),ygood(ip),resx(ip),resy(ip) + enddo + endif -******retrieve ANGULAR BINNING info - fname_binning='binning.dat' - print *,'Opening file: ',fname_binning - open(10, - $ FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning)) - $ ,STATUS='UNKNOWN' - $ ,IOSTAT=iostat - $ ) - if(iostat.ne.0)then - print*,'READETAPARAM: *** Error in opening file ***' - return - endif - print*,'---- ANGULAR BINNING ----' - print*,'Bin - angL - angR' - 101 format(i2,' ',f6.2,' ',f6.2) - do ibin=1,nangmax - read(10,* - $ ,IOSTAT=iostat - $ )xnn,angL(ibin),angR(ibin) - if(iostat.ne.0)goto 1000 - write(*,101)int(xnn),angL(ibin),angR(ibin) - enddo - 1000 nangbin=int(xnn) - close(10) - print*,'-------------------------' - +c rchi2=chi2/dble(ndof) + if(DEBUG)then + print*,' ' + print*,'****** SELECTED TRACK *************' + print*,'# R. chi2 RIG' + print*,' --- ',chi2,' --- ' + $ ,1./abs(AL(5)) + print*,'***********************************' + endif +* ********************************************************** +* ************************** FIT *** FIT *** FIT *** FIT *** +* ********************************************************** - do ieta=2,4 !loop on eta 2,3,4 -******retrieve correction parameters - 200 format(' Opening eta',i1,' files...') - write(*,200)ieta - - 201 format('eta',i1,'-bin',i1,'-l',i1,'.dat') - 202 format('eta',i1,'-bin',i2,'-l',i1,'.dat') - do iang=1,nangbin - do ilad=1,nladders_view - if(iang.lt.10)write(fname_param,201)ieta,iang,ilad - if(iang.ge.10)write(fname_param,202)ieta,iang,ilad -c print *,'Opening file: ',fname_param - open(10, - $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) - $ ,STATUS='UNKNOWN' - $ ,IOSTAT=iostat - $ ) - if(iostat.ne.0)then - print*,'READETAPARAM: *** Error in opening file ***' - return - endif - do ival=1,netavalmax - if(ieta.eq.2)read(10,* - $ ,IOSTAT=iostat - $ ) - $ eta2(ival,iang), - $ (feta2(ival,iv,ilad,iang),iv=1,nviews) - if(ieta.eq.3)read(10,* - $ ,IOSTAT=iostat - $ ) - $ eta3(ival,iang), - $ (feta3(ival,iv,ilad,iang),iv=1,nviews) - if(ieta.eq.4)read(10,* - $ ,IOSTAT=iostat - $ ) - $ eta4(ival,iang), - $ (feta4(ival,iv,ilad,iang),iv=1,nviews) - if(iostat.ne.0)then - netaval=ival-1 -c$$$ if(eta2(1,iang).ne.-eta2(netaval,iang)) -c$$$ $ print*,'**** ERROR on parameters !!! ****' - goto 2000 - endif - enddo - 2000 close(10) -* print*,'... done' +* ------------- search if the track has an IMAGE ------------- +* ------------- (also this is stored ) ------------- + if(FIMAGE)goto 122 !>>> jump! (this is already an image) +* now search for track-image, by comparing couples IDs + do i=1,ntracks + iimage=i + do ip=1,nplanes + if( CP_STORE(nplanes-ip+1,icand).ne. + $ -1*CP_STORE(nplanes-ip+1,i) )iimage=0 enddo + if( iimage.ne.0.and. +c $ RCHI2_STORE(i).le.CHI2MAX.and. +c $ RCHI2_STORE(i).gt.0.and. + $ .true.)then + if(DEBUG)print*,'Track candidate ',iimage + $ ,' >>> TRACK IMAGE >>> of' + $ ,ibest + goto 122 !image track found + endif enddo + 122 continue - enddo !end loop on eta 2,3,4 +* --- and store the results -------------------------------- + ntrk = ntrk + 1 !counter of found tracks + if(.not.FIMAGE + $ .and.iimage.eq.0) image(ntrk)= 0 + if(.not.FIMAGE + $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next + if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous + + call fill_level2_tracks(ntrk) !==> good2=.true. +c print*,'++++++++++ iimage,fimage,ntrk,image ' +c $ ,iimage,fimage,ntrk,image(ntrk) + + if(ntrk.eq.NTRKMAX)then + if(DEBUG) + $ print*, + $ '** warning ** number of identified '// + $ 'tracks exceeds vector dimension ' + $ ,'( ',NTRKMAX,' )' +cc good2=.false. + goto 880 !fill ntp and go to next event + endif + if(iimage.ne.0)then + FIMAGE=.true. ! + goto 1212 !>>> fit image-track + endif +* --- then remove selected clusters (ibest+iimage) from clouds ---- + call clean_XYclouds(ibest,iflag) + if(iflag.eq.1)then !bad event + goto 880 !fill ntp and go to next event + endif - return +* ********************************************************** +* condition to start a new search +* ********************************************************** + ixznew=0 + do ixz=1,nclouds_xz + if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1 + enddo + iyznew=0 + do iyz=1,nclouds_yz + if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1 + enddo + + if(ixznew.ne.0.and. + $ iyznew.ne.0.and. + $ rchi2best.le.CHI2MAX.and. +c $ rchi2best.lt.15..and. + $ .true.)then + if(DEBUG)then + print*,'***** NEW SEARCH ****' + endif + goto 11111 !try new search + + endif +* ********************************************** + + + + 880 return end + + +c$$$************************************************************ +c$$$ +c$$$ subroutine readmipparam +c$$$ +c$$$ include '../common/commontracker.f' +c$$$ include '../common/calib.f' +c$$$ +c$$$ character*60 fname_param +c$$$ 201 format('trk-LADDER',i1,'-mip.dat') +c$$$ do ilad=1,nladders_view +c$$$ write(fname_param,201)ilad +c$$$ print *,'Opening file: ',fname_param +c$$$ open(10, +c$$$ $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) +c$$$ $ ,STATUS='UNKNOWN' +c$$$ $ ,IOSTAT=iostat +c$$$ $ ) +c$$$ if(iostat.ne.0)then +c$$$ print*,'READMIPPARAM: *** Error in opening file ***' +c$$$ return +c$$$ endif +c$$$ do iv=1,nviews +c$$$ read(10,* +c$$$ $ ,IOSTAT=iostat +c$$$ $ )pip, +c$$$ $ mip(int(pip),ilad) +c$$$c print*,ilad,iv,pip,mip(int(pip),ilad) +c$$$ enddo +c$$$ close(10) +c$$$ enddo +c$$$ +c$$$ return +c$$$ end +c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** +c$$$ subroutine readchargeparam +c$$$ +c$$$ +c$$$ include '../common/commontracker.f' +c$$$ include '../common/calib.f' +c$$$ +c$$$ character*60 fname_param +c$$$ 201 format('charge-l',i1,'.dat') +c$$$ do ilad=1,nladders_view +c$$$ write(fname_param,201)ilad +c$$$ print *,'Opening file: ',fname_param +c$$$ open(10, +c$$$ $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) +c$$$ $ ,STATUS='UNKNOWN' +c$$$ $ ,IOSTAT=iostat +c$$$ $ ) +c$$$ if(iostat.ne.0)then +c$$$ print*,'READCHARGEPARAM: *** Error in opening file ***' +c$$$ return +c$$$ endif +c$$$ do ip=1,nplanes +c$$$ read(10,* +c$$$ $ ,IOSTAT=iostat +c$$$ $ )pip, +c$$$ $ kch(ip,ilad),cch(ip,ilad),sch(ip,ilad) +c$$$c print*,ilad,ip,pip,kch(ip,ilad), +c$$$c $ cch(ip,ilad),sch(ip,ilad) +c$$$ enddo +c$$$ close(10) +c$$$ enddo +c$$$ +c$$$ return +c$$$ end +c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** +c$$$ subroutine readetaparam +c$$$* ----------------------------------------- +c$$$* read eta2,3,4 calibration parameters +c$$$* and fill variables: +c$$$* +c$$$* eta2(netabin,nladders_view,nviews) +c$$$* eta3(2*netabin,nladders_view,nviews) +c$$$* eta4(2*netabin,nladders_view,nviews) +c$$$* +c$$$ include '../common/commontracker.f' +c$$$ include '../common/calib.f' +c$$$ +c$$$ character*40 fname_binning +c$$$ character*40 fname_param +c$$$c character*120 cmd1 +c$$$c character*120 cmd2 +c$$$ +c$$$ +c$$$******retrieve ANGULAR BINNING info +c$$$ fname_binning='binning.dat' +c$$$ print *,'Opening file: ',fname_binning +c$$$ open(10, +c$$$ $ FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning)) +c$$$ $ ,STATUS='UNKNOWN' +c$$$ $ ,IOSTAT=iostat +c$$$ $ ) +c$$$ if(iostat.ne.0)then +c$$$ print*,'READETAPARAM: *** Error in opening file ***' +c$$$ return +c$$$ endif +c$$$ print*,'---- ANGULAR BINNING ----' +c$$$ print*,'Bin - angL - angR' +c$$$ 101 format(i2,' ',f6.2,' ',f6.2) +c$$$ do ibin=1,nangmax +c$$$ read(10,* +c$$$ $ ,IOSTAT=iostat +c$$$ $ )xnn,angL(ibin),angR(ibin) +c$$$ if(iostat.ne.0)goto 1000 +c$$$ write(*,101)int(xnn),angL(ibin),angR(ibin) +c$$$ enddo +c$$$ 1000 nangbin=int(xnn) +c$$$ close(10) +c$$$ print*,'-------------------------' +c$$$ +c$$$ +c$$$ +c$$$ do ieta=2,4 !loop on eta 2,3,4 +c$$$******retrieve correction parameters +c$$$ 200 format(' Opening eta',i1,' files...') +c$$$ write(*,200)ieta +c$$$ +c$$$ 201 format('eta',i1,'-bin',i1,'-l',i1,'.dat') +c$$$ 202 format('eta',i1,'-bin',i2,'-l',i1,'.dat') +c$$$ do iang=1,nangbin +c$$$ do ilad=1,nladders_view +c$$$ if(iang.lt.10)write(fname_param,201)ieta,iang,ilad +c$$$ if(iang.ge.10)write(fname_param,202)ieta,iang,ilad +c$$$c print *,'Opening file: ',fname_param +c$$$ open(10, +c$$$ $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) +c$$$ $ ,STATUS='UNKNOWN' +c$$$ $ ,IOSTAT=iostat +c$$$ $ ) +c$$$ if(iostat.ne.0)then +c$$$ print*,'READETAPARAM: *** Error in opening file ***' +c$$$ return +c$$$ endif +c$$$ do ival=1,netavalmax +c$$$ if(ieta.eq.2)read(10,* +c$$$ $ ,IOSTAT=iostat +c$$$ $ ) +c$$$ $ eta2(ival,iang), +c$$$ $ (feta2(ival,iv,ilad,iang),iv=1,nviews) +c$$$ if(ieta.eq.3)read(10,* +c$$$ $ ,IOSTAT=iostat +c$$$ $ ) +c$$$ $ eta3(ival,iang), +c$$$ $ (feta3(ival,iv,ilad,iang),iv=1,nviews) +c$$$ if(ieta.eq.4)read(10,* +c$$$ $ ,IOSTAT=iostat +c$$$ $ ) +c$$$ $ eta4(ival,iang), +c$$$ $ (feta4(ival,iv,ilad,iang),iv=1,nviews) +c$$$ if(iostat.ne.0)then +c$$$ netaval=ival-1 +c$$$c$$$ if(eta2(1,iang).ne.-eta2(netaval,iang)) +c$$$c$$$ $ print*,'**** ERROR on parameters !!! ****' +c$$$ goto 2000 +c$$$ endif +c$$$ enddo +c$$$ 2000 close(10) +c$$$* print*,'... done' +c$$$ enddo +c$$$ enddo +c$$$ +c$$$ enddo !end loop on eta 2,3,4 +c$$$ +c$$$ +c$$$ return +c$$$ end +c$$$ + ************************************************************ ************************************************************