--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/09/28 14:04:40 1.4 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/11/08 16:42:28 1.13 @@ -12,19 +12,17 @@ subroutine track_finding(iflag) 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' +c include 'level1.f' include 'level2.f' - include 'momanhough_init.f' +c include 'momanhough_init.f' -c logical DEBUG -c common/dbg/DEBUG - *------------------------------------------------------------------------------- * STEP 1 *------------------------------------------------------------------------------- @@ -47,22 +45,9 @@ c iflag=0 call cl_to_couples(iflag) if(iflag.eq.1)then !bad event - goto 880 !fill ntp and go to next event + goto 880 !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 @@ -94,7 +79,7 @@ c iflag=0 call cp_to_doubtrip(iflag) if(iflag.eq.1)then !bad event - goto 880 !fill ntp and go to next event + goto 880 !go to next event endif @@ -123,18 +108,61 @@ * $ ,ptcloud_xz,tr_cloud,cpcloud_xz *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- +* count number of hit planes + planehit=0 + do np=1,nplanes + if(ncp_plane(np).ne.0)then + planehit=planehit+1 + endif + enddo + if(planehit.lt.3) goto 880 ! exit + + nptxz_min=x_min_start + nplxz_min=x_min_start + + nptyz_min=y_min_start + nplyz_min=y_min_start -c iflag=0 + cutdistyz=cutystart + cutdistxz=cutxstart + + 878 continue 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 + endif + if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then + if(cutdistyz.lt.maxcuty/2)then + cutdistyz=cutdistyz+cutystep + else + cutdistyz=cutdistyz+(3*cutystep) + endif + goto 878 + endif + + if(planehit.eq.3) goto 881 + + 879 continue call trip_to_XZcloud(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif - + + if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then + cutdistxz=cutdistxz+cutxstep + goto 879 + endif + + + 881 continue +* if there is at least three planes on the Y view decreases cuts on X view + if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and. + $ nplxz_min.ne.y_min_start)then + nptxz_min=x_min_step + nplxz_min=x_min_start-x_min_step + goto 879 + endif + 880 return end @@ -144,20 +172,18 @@ subroutine track_fitting(iflag) 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' - include 'momanhough_init.f' +c include 'momanhough_init.f' -c logical DEBUG -c common/dbg/DEBUG - logical FIMAGE ! + real*8 AL_GUESS(5) *------------------------------------------------------------------------------- * STEP 4 (ITERATED until any other physical track isn't found) @@ -198,13 +224,33 @@ ibest=0 !best track among candidates iimage=0 !track image * ------------- select the best track ------------- +c$$$ rchi2best=1000000000. +c$$$ do i=1,ntracks +c$$$ if(RCHI2_STORE(i).lt.rchi2best.and. +c$$$ $ RCHI2_STORE(i).gt.0)then +c$$$ ibest=i +c$$$ rchi2best=RCHI2_STORE(i) +c$$$ endif +c$$$ enddo +c$$$ if(ibest.eq.0)goto 880 !>> no good candidates + rchi2best=1000000000. + ndofbest=0 !(1) do i=1,ntracks - if(RCHI2_STORE(i).lt.rchi2best.and. - $ RCHI2_STORE(i).gt.0)then + if(RCHI2_STORE(i).lt.rchi2best.and. + $ RCHI2_STORE(i).gt.0)then + ndof=0 !(1) + do ii=1,nplanes !(1) + ndof=ndof !(1) + $ +int(xgood_store(ii,i)) !(1) + $ +int(ygood_store(ii,i)) !(1) + enddo !(1) + if(ndof.ge.ndofbest)then !(1) ibest=i rchi2best=RCHI2_STORE(i) - endif + ndofbest=ndof !(1) + endif !(1) + endif enddo if(ibest.eq.0)goto 880 !>> no good candidates *------------------------------------------------------------------------------- @@ -236,21 +282,37 @@ * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** + call guess() + do i=1,5 + AL_GUESS(i)=AL(i) + enddo + do i=1,5 AL(i)=dble(AL_STORE(i,icand)) enddo + IDCAND = icand !fitted track-candidate ifail=0 !error flag in chi2 computation jstep=0 !# minimization steps - call mini_2(jstep,ifail) + iprint=0 +c if(DEBUG)iprint=1 + if(VERBOSE)iprint=1 + call mini2(jstep,ifail,iprint) if(ifail.ne.0) then - if(DEBUG)then + if(.true.)then print *, - $ '*** MINIMIZATION FAILURE *** (mini_2) ' + $ '*** MINIMIZATION FAILURE *** (after refinement) ' $ ,iev + +c$$$ print*,'guess: ',(al_guess(i),i=1,5) +c$$$ print*,'previous: ',(al_store(i,icand),i=1,5) +c$$$ print*,'result: ',(al(i),i=1,5) +c$$$ print*,'xgood ',xgood +c$$$ print*,'ygood ',ygood +c$$$ print*,'----------------------------------------------' endif - chi2=-chi2 +c chi2=-chi2 endif if(DEBUG)then @@ -311,7 +373,7 @@ c $ ,iimage,fimage,ntrk,image(ntrk) if(ntrk.eq.NTRKMAX)then - if(DEBUG) + if(verbose) $ print*, $ '** warning ** number of identified '// $ 'tracks exceeds vector dimension ' @@ -607,12 +669,13 @@ c***************************************************** include 'commontracker.f' - include 'calib.f' include 'level1.f' + include 'calib.f' +c include 'level1.f' include 'common_align.f' include 'common_mech.f' include 'common_xyzPAM.f' - include 'common_resxy.f' +c include 'common_resxy.f' c logical DEBUG c common/dbg/DEBUG @@ -667,27 +730,27 @@ resxPAM = resxPAM*fbad_cog(2,icx) elseif(PFAx.eq.'ETA2')then c cog2 = cog(2,icx) -c etacorr = pfa_eta2(cog2,viewx,nldx,angx) +c etacorr = pfaeta2(cog2,viewx,nldx,angx) c stripx = stripx + etacorr - stripx = stripx + pfa_eta2(icx,angx) !(3) + stripx = stripx + pfaeta2(icx,angx) !(3) resxPAM = risx_eta2(angx) ! (4) if(DEBUG.and.fbad_cog(2,icx).ne.1) $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) resxPAM = resxPAM*fbad_cog(2,icx) elseif(PFAx.eq.'ETA3')then !(3) - stripx = stripx + pfa_eta3(icx,angx) !(3) + stripx = stripx + pfaeta3(icx,angx) !(3) resxPAM = risx_eta3(angx) ! (4) if(DEBUG.and.fbad_cog(3,icx).ne.1) !(3) $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3) resxPAM = resxPAM*fbad_cog(3,icx) !(3) elseif(PFAx.eq.'ETA4')then !(3) - stripx = stripx + pfa_eta4(icx,angx) !(3) + stripx = stripx + pfaeta4(icx,angx) !(3) resxPAM = risx_eta4(angx) ! (4) if(DEBUG.and.fbad_cog(4,icx).ne.1) !(3) $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3) resxPAM = resxPAM*fbad_cog(4,icx) !(3) elseif(PFAx.eq.'ETA')then !(3) - stripx = stripx + pfa_eta(icx,angx) !(3) + stripx = stripx + pfaeta(icx,angx) !(3) resxPAM = ris_eta(icx,angx) ! (4) if(DEBUG.and.fbad_cog(2,icx).ne.1) !(3) $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3) @@ -731,25 +794,25 @@ resyPAM = resyPAM*fbad_cog(2,icy) elseif(PFAy.eq.'ETA2')then c cog2 = cog(2,icy) -c etacorr = pfa_eta2(cog2,viewy,nldy,angy) +c etacorr = pfaeta2(cog2,viewy,nldy,angy) c stripy = stripy + etacorr - stripy = stripy + pfa_eta2(icy,angy) !(3) + stripy = stripy + pfaeta2(icy,angy) !(3) resyPAM = risy_eta2(angy) ! (4) resyPAM = resyPAM*fbad_cog(2,icy) if(DEBUG.and.fbad_cog(2,icy).ne.1) $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy) elseif(PFAy.eq.'ETA3')then !(3) - stripy = stripy + pfa_eta3(icy,angy) !(3) + stripy = stripy + pfaeta3(icy,angy) !(3) resyPAM = resyPAM*fbad_cog(3,icy) !(3) if(DEBUG.and.fbad_cog(3,icy).ne.1) !(3) $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3) elseif(PFAy.eq.'ETA4')then !(3) - stripy = stripy + pfa_eta4(icy,angy) !(3) + stripy = stripy + pfaeta4(icy,angy) !(3) resyPAM = resyPAM*fbad_cog(4,icy) !(3) if(DEBUG.and.fbad_cog(4,icy).ne.1) !(3) $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3) elseif(PFAy.eq.'ETA')then !(3) - stripy = stripy + pfa_eta(icy,angy) !(3) + stripy = stripy + pfaeta(icy,angy) !(3) resyPAM = ris_eta(icy,angy) ! (4) c resyPAM = resyPAM*fbad_cog(2,icy) !(3)TEMPORANEO resyPAM = resyPAM*fbad_eta(icy,angy) ! (4) @@ -1284,6 +1347,7 @@ * it returns the plane number * include 'commontracker.f' + include 'level1.f' c include 'common_analysis.f' include 'common_momanhough.f' @@ -1321,6 +1385,7 @@ * it returns the id number ON THE PLANE * include 'commontracker.f' + include 'level1.f' c include 'common_analysis.f' include 'common_momanhough.f' @@ -1349,6 +1414,7 @@ * positive if sensor =2 * include 'commontracker.f' + include 'level1.f' c include 'calib.f' c include 'level1.f' c include 'common_analysis.f' @@ -1660,13 +1726,11 @@ subroutine cl_to_couples(iflag) include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' - include 'momanhough_init.f' +c include 'momanhough_init.f' include 'calib.f' - include 'level1.f' - -c logical DEBUG -c common/dbg/DEBUG +c include 'level1.f' * output flag * -------------- @@ -1675,7 +1739,7 @@ * -------------- integer iflag - integer badseed,badcl + integer badseed,badclx,badcly * init variables ncp_tot=0 @@ -1691,16 +1755,38 @@ ncls(ip)=0 enddo do icl=1,nclstrmax_level2 - cl_single(icl)=1 - cl_good(icl)=0 + cl_single(icl) = 1 + cl_good(icl) = 0 + enddo + do iv=1,nviews + ncl_view(iv) = 0 + mask_view(iv) = 0 !all included enddo +* count number of cluster per view + do icl=1,nclstr1 + ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1 + enddo +* mask views with too many clusters + do iv=1,nviews + if( ncl_view(iv).gt. nclustermax)then + mask_view(iv) = 1 + if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > ' + $ ,nclustermax,' on view ', iv,' --> masked!' + endif + enddo + + * start association ncouples=0 do icx=1,nclstr1 !loop on cluster (X) if(mod(VIEW(icx),2).eq.1)goto 10 * ---------------------------------------------------- +* jump masked views (X VIEW) +* ---------------------------------------------------- + if( mask_view(VIEW(icx)).ne.0 ) goto 10 +* ---------------------------------------------------- * cut on charge (X VIEW) * ---------------------------------------------------- if(dedx(icx).lt.dedx_x_min)then @@ -1717,7 +1803,7 @@ else ilast=TOTCLLENGTH endif - badcl=badseed + badclx=badseed do igood=-ngoodstr,ngoodstr ibad=1 if((INDMAX(icx)+igood).gt.ifirst.and. @@ -1727,7 +1813,7 @@ $ nvk(MAXS(icx)+igood), $ nst(MAXS(icx)+igood)) endif - badcl=badcl*ibad + badclx=badclx*ibad enddo * ---------------------------------------------------- * >>> eliminato il taglio sulle BAD <<< @@ -1746,6 +1832,11 @@ if(mod(VIEW(icy),2).eq.0)goto 20 * ---------------------------------------------------- +* jump masked views (Y VIEW) +* ---------------------------------------------------- + if( mask_view(VIEW(icy)).ne.0 ) goto 20 + +* ---------------------------------------------------- * cut on charge (Y VIEW) * ---------------------------------------------------- if(dedx(icy).lt.dedx_y_min)then @@ -1762,7 +1853,7 @@ else ilast=TOTCLLENGTH endif - badcl=badseed + badcly=badseed do igood=-ngoodstr,ngoodstr ibad=1 if((INDMAX(icy)+igood).gt.ifirst.and. @@ -1771,7 +1862,7 @@ $ ibad=BAD(VIEW(icy), $ nvk(MAXS(icy)+igood), $ nst(MAXS(icy)+igood)) - badcl=badcl*ibad + badcly=badcly*ibad enddo * ---------------------------------------------------- * >>> eliminato il taglio sulle BAD <<< @@ -1794,43 +1885,58 @@ * charge correlation * (modified to be applied only below saturation... obviously) -* ------------------------------------------------------------- -* >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<< -* ------------------------------------------------------------- -c$$$ if(dedx(icy).lt.chsaty.or.dedx(icx).lt.chsatx)then -c$$$ ddd=(dedx(icy) -c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) -c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1) -c$$$ cut=chcut*sch(nplx,nldx) -c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent -c$$$ endif + if( .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx) + $ .and. + $ .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx) + $ .and. + $ (badclx.eq.1.and.badcly.eq.1) + $ .and. + $ .true.)then + + ddd=(dedx(icy) + $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) + ddd=ddd/sqrt(kch(nplx,nldx)**2+1) + +c cut = chcut * sch(nplx,nldx) + + sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx) + $ -kch(nplx,nldx)*cch(nplx,nldx)) + sss=sss/sqrt(kch(nplx,nldx)**2+1) + cut = chcut * (16 + sss/50.) + + if(abs(ddd).gt.cut)then + goto 20 !charge not consistent + endif + endif * ------------------> COUPLE <------------------ * check to do not overflow vector dimentions - if(ncp_plane(nplx).gt.ncouplemax)then - if(DEBUG)print*, - $ ' ** warning ** number of identified'// - $ ' couples on plane ',nplx, - $ ' exceeds vector dimention'// - $ ' ( ',ncouplemax,' )' -c good2=.false. -c goto 880 !fill ntp and go to next event - iflag=1 - return - endif - -c$$$ if(ncp_plane(nplx).eq.ncouplemax)then +c$$$ if(ncp_plane(nplx).gt.ncouplemax)then c$$$ if(DEBUG)print*, -c$$$ $ '** warning ** number of identified '// -c$$$ $ 'couples on plane ',nplx, -c$$$ $ 'exceeds vector dimention ' -c$$$ $ ,'( ',ncouplemax,' )' +c$$$ $ ' ** warning ** number of identified'// +c$$$ $ ' couples on plane ',nplx, +c$$$ $ ' exceeds vector dimention'// +c$$$ $ ' ( ',ncouplemax,' )' c$$$c good2=.false. -c$$$c goto 880 !fill ntp and go to next event +c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif + if(ncp_plane(nplx).eq.ncouplemax)then + if(verbose)print*, + $ '** warning ** number of identified '// + $ 'couples on plane ',nplx, + $ 'exceeds vector dimention ' + $ ,'( ',ncouplemax,' )' +c good2=.false. +c goto 880 !fill ntp and go to next event + mask_view(nviewx(nplx)) = 2 + mask_view(nviewy(nply)) = 2 +c iflag=1 +c return + endif + ncp_plane(nplx) = ncp_plane(nplx) + 1 clx(nplx,ncp_plane(nplx))=icx cly(nply,ncp_plane(nplx))=icy @@ -1863,20 +1969,18 @@ endif do ip=1,6 - ncp_tot=ncp_tot+ncp_plane(ip) + ncp_tot = ncp_tot + ncp_plane(ip) enddo c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) - if(ncp_tot.gt.ncp_max)then - if(DEBUG)print*, - $ '** warning ** number of identified '// - $ 'couples exceeds upper limit for Hough tr. ' - $ ,'( ',ncp_max,' )' -c good2=.false. -c goto 880 !fill ntp and go to next event - iflag=1 - return - endif +c$$$ if(ncp_tot.gt.ncp_max)then +c$$$ if(verbose)print*, +c$$$ $ '** warning ** number of identified '// +c$$$ $ 'couples exceeds upper limit for Hough tr. ' +c$$$ $ ,'( ',ncp_max,' )' +c$$$ iflag=1 +c$$$ return +c$$$ endif return end @@ -1889,228 +1993,15 @@ * * * * ************************************************** - subroutine cl_to_couples_nocharge(iflag) - - include 'commontracker.f' - include 'common_momanhough.f' - include 'momanhough_init.f' - include 'calib.f' - include 'level1.f' - -c logical DEBUG -c common/dbg/DEBUG - -* output flag -* -------------- -* 0 = good event -* 1 = bad event -* -------------- - integer iflag - - integer badseed,badcl - -* init variables - ncp_tot=0 - do ip=1,nplanes - do ico=1,ncouplemax - clx(ip,ico)=0 - cly(ip,ico)=0 - enddo - ncp_plane(ip)=0 - do icl=1,nclstrmax_level2 - cls(ip,icl)=1 - enddo - ncls(ip)=0 - enddo - do icl=1,nclstrmax_level2 - cl_single(icl)=1 - cl_good(icl)=0 - enddo - -* start association - ncouples=0 - do icx=1,nclstr1 !loop on cluster (X) - if(mod(VIEW(icx),2).eq.1)goto 10 - -* ---------------------------------------------------- -* cut on charge (X VIEW) - if(dedx(icx).lt.dedx_x_min)then - cl_single(icx)=0 - goto 10 - endif -* cut BAD (X VIEW) - badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx))) - ifirst=INDSTART(icx) - if(icx.ne.nclstr1) then - ilast=INDSTART(icx+1)-1 - else - ilast=TOTCLLENGTH - endif - badcl=badseed - do igood=-ngoodstr,ngoodstr - ibad=1 - if((INDMAX(icx)+igood).gt.ifirst.and. - $ (INDMAX(icx)+igood).lt.ilast.and. - $ .true.)then - ibad=BAD(VIEW(icx), - $ nvk(MAXS(icx)+igood), - $ nst(MAXS(icx)+igood)) - endif - badcl=badcl*ibad - enddo - if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut - cl_single(icx)=0 !<<<<<<<<<<<<<< BAD cut - goto 10 !<<<<<<<<<<<<<< BAD cut - endif !<<<<<<<<<<<<<< BAD cut -* ---------------------------------------------------- - - cl_good(icx)=1 - nplx=npl(VIEW(icx)) - nldx=nld(MAXS(icx),VIEW(icx)) - - do icy=1,nclstr1 !loop on cluster (Y) - if(mod(VIEW(icy),2).eq.0)goto 20 - -* ---------------------------------------------------- -* cut on charge (Y VIEW) - if(dedx(icy).lt.dedx_y_min)then - cl_single(icy)=0 - goto 20 - endif -* cut BAD (Y VIEW) - badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy))) - ifirst=INDSTART(icy) - if(icy.ne.nclstr1) then - ilast=INDSTART(icy+1)-1 - else - ilast=TOTCLLENGTH - endif - badcl=badseed - do igood=-ngoodstr,ngoodstr - ibad=1 - if((INDMAX(icy)+igood).gt.ifirst.and. - $ (INDMAX(icy)+igood).lt.ilast.and. - $ .true.) - $ ibad=BAD(VIEW(icy), - $ nvk(MAXS(icy)+igood), - $ nst(MAXS(icy)+igood)) - badcl=badcl*ibad - enddo - if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut - cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut - goto 20 !<<<<<<<<<<<<<< BAD cut - endif !<<<<<<<<<<<<<< BAD cut -* ---------------------------------------------------- - - - cl_good(icy)=1 - nply=npl(VIEW(icy)) - nldy=nld(MAXS(icy),VIEW(icy)) - -* ---------------------------------------------- -* CONDITION TO FORM A COUPLE -* ---------------------------------------------- -* geometrical consistency (same plane and ladder) - if(nply.eq.nplx.and.nldy.eq.nldx)then -* charge correlation -* =========================================================== -* this version of the subroutine is used for the calibration -* thus charge-correlation selection is obviously removed -* =========================================================== -c$$$ ddd=(dedx(icy) -c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) -c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1) -c$$$ cut=chcut*sch(nplx,nldx) -c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent -* =========================================================== - - -* ------------------> COUPLE <------------------ -* check to do not overflow vector dimentions - if(ncp_plane(nplx).gt.ncouplemax)then - if(DEBUG)print*, - $ ' ** warning ** number of identified'// - $ ' couples on plane ',nplx, - $ ' exceeds vector dimention'// - $ ' ( ',ncouplemax,' )' -c good2=.false. -c goto 880 !fill ntp and go to next event - iflag=1 - return - endif - -c$$$ if(ncp_plane(nplx).eq.ncouplemax)then -c$$$ if(DEBUG)print*, -c$$$ $ '** warning ** number of identified '// -c$$$ $ 'couples on plane ',nplx, -c$$$ $ 'exceeds vector dimention ' -c$$$ $ ,'( ',ncouplemax,' )' -c$$$c good2=.false. -c$$$c goto 880 !fill ntp and go to next event -c$$$ iflag=1 -c$$$ return -c$$$ endif - - ncp_plane(nplx) = ncp_plane(nplx) + 1 - clx(nplx,ncp_plane(nplx))=icx - cly(nply,ncp_plane(nplx))=icy - cl_single(icx)=0 - cl_single(icy)=0 - endif -* ---------------------------------------------- - - 20 continue - enddo !end loop on clusters(Y) - - 10 continue - enddo !end loop on clusters(X) - - - do icl=1,nclstr1 - if(cl_single(icl).eq.1)then - ip=npl(VIEW(icl)) - ncls(ip)=ncls(ip)+1 - cls(ip,ncls(ip))=icl - endif - enddo - - - if(DEBUG)then - print*,'clusters ',nclstr1 - print*,'good ',(cl_good(i),i=1,nclstr1) - print*,'singles ',(cl_single(i),i=1,nclstr1) - print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) - endif - - do ip=1,6 - ncp_tot=ncp_tot+ncp_plane(ip) - enddo -c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) - - if(ncp_tot.gt.ncp_max)then - if(DEBUG)print*, - $ '** warning ** number of identified '// - $ 'couples exceeds upper limit for Hough tr. ' - $ ,'( ',ncp_max,' )' -c good2=.false. -c goto 880 !fill ntp and go to next event - iflag=1 - return - endif - - return - end - -c$$$ subroutine cl_to_couples_2(iflag) +c$$$ subroutine cl_to_couples_nocharge(iflag) c$$$ c$$$ include 'commontracker.f' +c$$$ include 'level1.f' c$$$ include 'common_momanhough.f' -c$$$ include 'momanhough_init.f' +c$$$c include 'momanhough_init.f' c$$$ include 'calib.f' -c$$$ include 'level1.f' +c$$$c include 'level1.f' c$$$ -c$$$ logical DEBUG -c$$$ common/dbg/DEBUG c$$$ c$$$* output flag c$$$* -------------- @@ -2170,11 +2061,10 @@ c$$$ endif c$$$ badcl=badcl*ibad c$$$ enddo -c$$$* print*,'icx ',icx,badcl -c$$$ if(badcl.eq.0)then -c$$$ cl_single(icx)=0 -c$$$ goto 10 -c$$$ endif +c$$$ if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut +c$$$ cl_single(icx)=0 !<<<<<<<<<<<<<< BAD cut +c$$$ goto 10 !<<<<<<<<<<<<<< BAD cut +c$$$ endif !<<<<<<<<<<<<<< BAD cut c$$$* ---------------------------------------------------- c$$$ c$$$ cl_good(icx)=1 @@ -2209,11 +2099,10 @@ c$$$ $ nst(MAXS(icy)+igood)) c$$$ badcl=badcl*ibad c$$$ enddo -c$$$* print*,'icy ',icy,badcl -c$$$ if(badcl.eq.0)then -c$$$ cl_single(icy)=0 -c$$$ goto 20 -c$$$ endif +c$$$ if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut +c$$$ cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut +c$$$ goto 20 !<<<<<<<<<<<<<< BAD cut +c$$$ endif !<<<<<<<<<<<<<< BAD cut c$$$* ---------------------------------------------------- c$$$ c$$$ @@ -2226,30 +2115,35 @@ c$$$* ---------------------------------------------- c$$$* geometrical consistency (same plane and ladder) c$$$ if(nply.eq.nplx.and.nldy.eq.nldx)then -c$$$ -c$$$c$$$* charge correlation +c$$$* charge correlation +c$$$* =========================================================== +c$$$* this version of the subroutine is used for the calibration +c$$$* thus charge-correlation selection is obviously removed +c$$$* =========================================================== c$$$c$$$ ddd=(dedx(icy) c$$$c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) c$$$c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1) c$$$c$$$ cut=chcut*sch(nplx,nldx) c$$$c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent +c$$$* =========================================================== +c$$$ c$$$ c$$$* ------------------> COUPLE <------------------ c$$$* check to do not overflow vector dimentions -c$$$ if(ncp_plane(nplx).gt.ncouplemax)then -c$$$ if(DEBUG)print*, -c$$$ $ ' ** warning ** number of identified'// -c$$$ $ ' couples on plane ',nplx, -c$$$ $ ' exceeds vector dimention'// -c$$$ $ ' ( ',ncouplemax,' )' -c$$$c good2=.false. -c$$$c goto 880 !fill ntp and go to next event -c$$$ iflag=1 -c$$$ return -c$$$ endif +c$$$c$$$ if(ncp_plane(nplx).gt.ncouplemax)then +c$$$c$$$ if(DEBUG)print*, +c$$$c$$$ $ ' ** warning ** number of identified'// +c$$$c$$$ $ ' couples on plane ',nplx, +c$$$c$$$ $ ' exceeds vector dimention'// +c$$$c$$$ $ ' ( ',ncouplemax,' )' +c$$$c$$$c good2=.false. +c$$$c$$$c goto 880 !fill ntp and go to next event +c$$$c$$$ iflag=1 +c$$$c$$$ return +c$$$c$$$ endif c$$$ c$$$ if(ncp_plane(nplx).eq.ncouplemax)then -c$$$ if(DEBUG)print*, +c$$$ if(verbose)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'couples on plane ',nplx, c$$$ $ 'exceeds vector dimention ' @@ -2265,7 +2159,6 @@ c$$$ cly(nply,ncp_plane(nplx))=icy c$$$ cl_single(icx)=0 c$$$ cl_single(icy)=0 -c$$$c print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy c$$$ endif c$$$* ---------------------------------------------- c$$$ @@ -2298,7 +2191,7 @@ c$$$c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) c$$$ c$$$ if(ncp_tot.gt.ncp_max)then -c$$$ if(DEBUG)print*, +c$$$ if(verbose)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'couples exceeds upper limit for Hough tr. ' c$$$ $ ,'( ',ncp_max,' )' @@ -2310,6 +2203,7 @@ c$$$ c$$$ return c$$$ end +c$$$ *************************************************** * * @@ -2326,15 +2220,14 @@ c***************************************************** include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' - include 'momanhough_init.f' +c include 'momanhough_init.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'calib.f' - include 'level1.f' +c include 'level1.f' -c logical DEBUG -c common/dbg/DEBUG * output flag * -------------- @@ -2402,12 +2295,15 @@ * (2 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ndblt.eq.ndblt_max)then - if(DEBUG)print*, + if(verbose)print*, $ '** warning ** number of identified '// $ 'doublets exceeds vector dimention ' $ ,'( ',ndblt_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event + do iv=1,12 + mask_view(iv) = 3 + enddo iflag=1 return endif @@ -2472,12 +2368,15 @@ * (3 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ntrpt.eq.ntrpt_max)then - if(DEBUG)print*, + if(verbose)print*, $ '** warning ** number of identified '// $ 'triplets exceeds vector dimention ' $ ,'( ',ntrpt_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event + do iv=1,nviews + mask_view(iv) = 4 + enddo iflag=1 return endif @@ -2552,11 +2451,10 @@ subroutine doub_to_YZcloud(iflag) include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' - include 'momanhough_init.f' +c include 'momanhough_init.f' -c logical DEBUG -c common/dbg/DEBUG * output flag * -------------- @@ -2588,6 +2486,8 @@ distance=0 nclouds_yz=0 !number of clouds npt_tot=0 + nloop=0 + 90 continue do idb1=1,ndblt !loop (1) on DOUBLETS if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud @@ -2691,7 +2591,7 @@ nplused=nplused+ hit_plane(ip) enddo c print*,'>>>> ',ncpused,npt,nplused - if(ncpused.lt.ncpyz_min)goto 2228 !next doublet +c if(ncpused.lt.ncpyz_min)goto 2228 !next doublet if(npt.lt.nptyz_min)goto 2228 !next doublet if(nplused.lt.nplyz_min)goto 2228 !next doublet @@ -2699,12 +2599,15 @@ * >>> NEW CLOUD <<< if(nclouds_yz.ge.ncloyz_max)then - if(DEBUG)print*, + if(verbose)print*, $ '** warning ** number of identified '// $ 'YZ clouds exceeds vector dimention ' $ ,'( ',ncloyz_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event + do iv=1,nviews + mask_view(iv) = 5 + enddo iflag=1 return endif @@ -2742,6 +2645,12 @@ enddo !end loop (1) on DOUBLETS + if(nloop.lt.nstepy)then + cutdistyz = cutdistyz+cutystep + nloop = nloop+1 + goto 90 + endif + if(DEBUG)then print*,'---------------------- ' print*,'Y-Z total clouds ',nclouds_yz @@ -2768,11 +2677,10 @@ subroutine trip_to_XZcloud(iflag) include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' - include 'momanhough_init.f' +c include 'momanhough_init.f' -c logical DEBUG -c common/dbg/DEBUG * output flag * -------------- @@ -2803,6 +2711,8 @@ distance=0 nclouds_xz=0 !number of clouds npt_tot=0 !total number of selected triplets + nloop=0 + 91 continue do itr1=1,ntrpt !loop (1) on TRIPLETS if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud c print*,'--------------' @@ -2904,19 +2814,22 @@ do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo - if(ncpused.lt.ncpxz_min)goto 22288 !next triplet +c if(ncpused.lt.ncpxz_min)goto 22288 !next triplet if(npt.lt.nptxz_min)goto 22288 !next triplet if(nplused.lt.nplxz_min)goto 22288 !next doublet * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< if(nclouds_xz.ge.ncloxz_max)then - if(DEBUG)print*, + if(verbose)print*, $ '** warning ** number of identified '// $ 'XZ clouds exceeds vector dimention ' $ ,'( ',ncloxz_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event + do iv=1,nviews + mask_view(iv) = 6 + enddo iflag=1 return endif @@ -2952,7 +2865,13 @@ * ~~~~~~~~~~~~~~~~~ 22288 continue enddo !end loop (1) on DOUBLETS - + + if(nloop.lt.nstepx)then + cutdistxz=cutdistxz+cutxstep + nloop=nloop+1 + goto 91 + endif + if(DEBUG)then print*,'---------------------- ' print*,'X-Z total clouds ',nclouds_xz @@ -2979,14 +2898,13 @@ c***************************************************** include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'common_mech.f' - include 'momanhough_init.f' +c include 'momanhough_init.f' -c logical DEBUG -c common/dbg/DEBUG * output flag * -------------- @@ -3002,14 +2920,14 @@ * ----------------------------------------------------------- * list of matching couples in the combination * between a XZ and YZ cloud - integer cp_match(nplanes,ncouplemax) + integer cp_match(nplanes,2*ncouplemax) integer ncp_match(nplanes) * ----------------------------------------------------------- integer hit_plane(nplanes) * ----------------------------------------------------------- * variables for track fitting double precision AL_INI(5) - double precision tath +c double precision tath * ----------------------------------------------------------- c real fitz(nplanes) !z coordinates of the planes in cm @@ -3102,23 +3020,33 @@ c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~' * -------> INITIAL GUESS <------- - AL_INI(1)=dreal(alfaxz1_av(ixz)) - AL_INI(2)=dreal(alfayz1_av(iyz)) - AL_INI(4)=datan(dreal(alfayz2_av(iyz)) - $ /dreal(alfaxz2_av(ixz))) - tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4)) - AL_INI(3)=tath/sqrt(1+tath**2) - AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0. - -c print*,'*******',AL_INI(5) - if(AL_INI(5).gt.defmax)goto 888 !next cloud - -c print*,'alfaxz2, alfayz2 ' -c $ ,alfaxz2_av(ixz),alfayz2_av(iyz) - -* -------> INITIAL GUESS <------- -c print*,'AL_INI ',(al_ini(i),i=1,5) - +cccc SBAGLIATO +c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz)) +c$$$ AL_INI(2) = dreal(alfayz1_av(iyz)) +c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz)) +c$$$ $ /dreal(alfaxz2_av(ixz))) +c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4)) +c$$$ AL_INI(3) = tath/sqrt(1+tath**2) +c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0. +cccc GIUSTO (ma si sua guess()) +c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz)) +c$$$ AL_INI(2) = dreal(alfayz1_av(iyz)) +c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4)) +c$$$ AL_INI(3) = tath/sqrt(1+tath**2) +c$$$ IF(alfaxz2_av(ixz).NE.0)THEN +c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz)) +c$$$ $ /dreal(alfaxz2_av(ixz))) +c$$$ ELSE +c$$$ AL_INI(4) = acos(-1.)/2 +c$$$ IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.) +c$$$ ENDIF +c$$$ IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4) +c$$$ AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs +c$$$ +c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0. +c$$$ +c$$$ if(AL_INI(5).gt.defmax)goto 888 !next cloud + if(DEBUG)then print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6)) print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5)) @@ -3186,19 +3114,34 @@ * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** +cccc scommentare se si usa al_ini della nuvola +c$$$ do i=1,5 +c$$$ AL(i)=AL_INI(i) +c$$$ enddo + call guess() do i=1,5 - AL(i)=AL_INI(i) + AL_INI(i)=AL(i) enddo ifail=0 !error flag in chi^2 computation jstep=0 !number of minimization steps - call mini_2(jstep,ifail) + iprint=0 +c if(DEBUG)iprint=1 + if(DEBUG)iprint=1 + call mini2(jstep,ifail,iprint) if(ifail.ne.0) then if(DEBUG)then print *, $ '*** MINIMIZATION FAILURE *** ' - $ //'(mini_2 in clouds_to_ctrack)' + $ //'(clouds_to_ctrack)' + print*,'initial guess: ' + + print*,'AL_INI(1) = ',AL_INI(1) + print*,'AL_INI(2) = ',AL_INI(2) + print*,'AL_INI(3) = ',AL_INI(3) + print*,'AL_INI(4) = ',AL_INI(4) + print*,'AL_INI(5) = ',AL_INI(5) endif - chi2=-chi2 +c chi2=-chi2 endif * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** @@ -3211,12 +3154,15 @@ * -------------------------- if(ntracks.eq.NTRACKSMAX)then - if(DEBUG)print*, + if(verbose)print*, $ '** warning ** number of candidate tracks '// $ ' exceeds vector dimension ' $ ,'( ',NTRACKSMAX,' )' c good2=.false. c goto 880 !fill ntp and go to next event + do iv=1,nviews + mask_view(iv) = 7 + enddo iflag=1 return endif @@ -3315,16 +3261,15 @@ c****************************************************** include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'common_mech.f' - include 'momanhough_init.f' - include 'level1.f' +c include 'momanhough_init.f' +c include 'level1.f' include 'calib.f' -c logical DEBUG -c common/dbg/DEBUG * flag to chose PFA character*10 PFA @@ -3338,11 +3283,13 @@ call track_init do ip=1,nplanes !loop on planes +* ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * If the plane has been already included, it just * computes again the coordinates of the x-y couple * using improved PFAs * ------------------------------------------------- +* ||||||||||||||||||||||||||||||||||||||||||||||||| if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and. $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then @@ -3377,10 +3324,12 @@ dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) +* ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * If the plane has NOT been already included, * it tries to include a COUPLE or a single cluster * ------------------------------------------------- +* ||||||||||||||||||||||||||||||||||||||||||||||||| else xgood(nplanes-ip+1)=0 @@ -3435,6 +3384,7 @@ $ AYV_STORE(nplanes-ip+1,ibest)) distance = distance_to(XP,YP) + distance = distance / RCHI2_STORE(ibest)!<<< MS id=id_cp(ip,icp,ist) if(DEBUG)print*,'( couple ',id $ ,' ) normalized distance ',distance @@ -3505,7 +3455,7 @@ $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest),0.) distance = distance_to(XP,YP) -c if(DEBUG)print*,'normalized distance ',distance + distance = distance / RCHI2_STORE(ibest)!<<< MS if(DEBUG)print*,'( cl-X ',icx $ ,' in cp ',id,' ) normalized distance ',distance if(distance.lt.distmin)then @@ -3533,6 +3483,7 @@ $ PFA,PFA, $ 0.,AYV_STORE(nplanes-ip+1,ibest)) distance = distance_to(XP,YP) + distance = distance / RCHI2_STORE(ibest)!<<< MS if(DEBUG)print*,'( cl-Y ',icy $ ,' in cp ',id,' ) normalized distance ',distance if(distance.lt.distmin)then @@ -3572,6 +3523,7 @@ endif distance = distance_to(XP,YP) + distance = distance / RCHI2_STORE(ibest)!<<< MS if(DEBUG)print*,'( cl-s ',icl $ ,' ) normalized distance ',distance if(distance.lt.distmin)then @@ -3601,17 +3553,25 @@ CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< * ---------------------------- +c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' if(mod(VIEW(iclm),2).eq.0)then XGOOD(nplanes-ip+1)=1. resx(nplanes-ip+1)=rxmm if(DEBUG)print*,'%%%% included X-cl ',iclm - $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )' +c if(.true.)print*,'%%%% included X-cl ',iclm + $ ,'( chi^2, ',RCHI2_STORE(ibest) + $ ,', norm.dist.= ',distmin + $ ,', cut ',clinc,' )' else YGOOD(nplanes-ip+1)=1. resy(nplanes-ip+1)=rymm if(DEBUG)print*,'%%%% included Y-cl ',iclm - $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )' +c if(.true.)print*,'%%%% included Y-cl ',iclm + $ ,'( chi^2, ',RCHI2_STORE(ibest) + $ ,', norm.dist.= ', distmin + $ ,', cut ',clinc,' )' endif +c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' * ---------------------------- xm_A(nplanes-ip+1) = xmm_A ym_A(nplanes-ip+1) = ymm_A @@ -3645,14 +3605,13 @@ subroutine clean_XYclouds(ibest,iflag) include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' - include 'momanhough_init.f' +c include 'momanhough_init.f' include 'level2.f' !(1) c include 'calib.f' c include 'level1.f' -c logical DEBUG -c common/dbg/DEBUG do ip=1,nplanes !loop on planes @@ -3828,47 +3787,21 @@ subroutine init_level2 -c***************************************************** -c 07/10/2005 modified by elena vannuccini --> (1) -c***************************************************** - include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' include 'level2.f' - include 'level1.f' +c include 'level1.f' do i=1,nviews good2(i)=good1(i) enddo -c good2 = 0!.false. -c$$$ nev2 = nev1 - -c$$$# ifndef TEST2003 -c$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$c pkt_type = pkt_type1 -c$$$c pkt_num = pkt_num1 -c$$$c obt = obt1 -c$$$c which_calib = which_calib1 -c$$$ swcode = 302 -c$$$ -c$$$ which_calib = which_calib1 -c$$$ pkt_type = pkt_type1 -c$$$ pkt_num = pkt_num1 -c$$$ obt = obt1 -c$$$ cpu_crc = cpu_crc1 -c$$$ do iv=1,12 -c$$$ crc(iv)=crc1(iv) -c$$$ enddo -c$$$# endif -c***************************************************** NTRK = 0 - do it=1,NTRKMAX!NTRACKSMAX + do it=1,NTRKMAX IMAGE(IT)=0 CHI2_nt(IT) = -100000. -c BdL(IT) = 0. do ip=1,nplanes XM_nt(IP,IT) = 0 YM_nt(IP,IT) = 0 @@ -3877,12 +3810,8 @@ RESY_nt(IP,IT) = 0 XGOOD_nt(IP,IT) = 0 YGOOD_nt(IP,IT) = 0 -c***************************************************** -cccccc 11/9/2005 modified by david fedele DEDX_X(IP,IT) = 0 DEDX_Y(IP,IT) = 0 -c****************************************************** -cccccc 17/8/2006 modified by elena CLTRX(IP,IT) = 0 CLTRY(IP,IT) = 0 enddo @@ -3893,25 +3822,18 @@ enddo enddo enddo - - -c***************************************************** -cccccc 11/9/2005 modified by david fedele nclsx=0 nclsy=0 do ip=1,NSINGMAX planex(ip)=0 -c xs(ip)=0 xs(1,ip)=0 xs(2,ip)=0 sgnlxs(ip)=0 planey(ip)=0 -c ys(ip)=0 ys(1,ip)=0 ys(2,ip)=0 sgnlys(ip)=0 enddo -c******************************************************* end @@ -3926,6 +3848,90 @@ ************************************************************ + subroutine init_hough + + include 'commontracker.f' + include 'level1.f' + include 'common_momanhough.f' + include 'common_hough.f' + include 'level2.f' + + ntrpt_nt=0 + ndblt_nt=0 + NCLOUDS_XZ_nt=0 + NCLOUDS_YZ_nt=0 + do idb=1,ndblt_max_nt + db_cloud_nt(idb)=0 + alfayz1_nt(idb)=0 + alfayz2_nt(idb)=0 + enddo + do itr=1,ntrpl_max_nt + tr_cloud_nt(itr)=0 + alfaxz1_nt(itr)=0 + alfaxz2_nt(itr)=0 + alfaxz3_nt(itr)=0 + enddo + do idb=1,ncloyz_max + ptcloud_yz_nt(idb)=0 + alfayz1_av_nt(idb)=0 + alfayz2_av_nt(idb)=0 + enddo + do itr=1,ncloxz_max + ptcloud_xz_nt(itr)=0 + alfaxz1_av_nt(itr)=0 + alfaxz2_av_nt(itr)=0 + alfaxz3_av_nt(itr)=0 + enddo + + ntrpt=0 + ndblt=0 + NCLOUDS_XZ=0 + NCLOUDS_YZ=0 + do idb=1,ndblt_max + db_cloud(idb)=0 + cpyz1(idb)=0 + cpyz2(idb)=0 + alfayz1(idb)=0 + alfayz2(idb)=0 + enddo + do itr=1,ntrpl_max + tr_cloud(itr)=0 + cpxz1(itr)=0 + cpxz2(itr)=0 + cpxz3(itr)=0 + alfaxz1(itr)=0 + alfaxz2(itr)=0 + alfaxz3(itr)=0 + enddo + do idb=1,ncloyz_max + ptcloud_yz(idb)=0 + alfayz1_av(idb)=0 + alfayz2_av(idb)=0 + do idbb=1,ncouplemaxtot + cpcloud_yz(idb,idbb)=0 + enddo + enddo + do itr=1,ncloxz_max + ptcloud_xz(itr)=0 + alfaxz1_av(itr)=0 + alfaxz2_av(itr)=0 + alfaxz3_av(itr)=0 + do itrr=1,ncouplemaxtot + cpcloud_xz(itr,itrr)=0 + enddo + enddo + end +************************************************************ +* +* +* +* +* +* +* +************************************************************ + + subroutine fill_level2_tracks(ntr) * ------------------------------------------------------- @@ -3936,36 +3942,35 @@ include 'commontracker.f' +c include 'level1.f' include 'level1.f' + include 'common_momanhough.f' include 'level2.f' include 'common_mini_2.f' - include 'common_momanhough.f' - real sinth,phi,pig !(4) + real sinth,phi,pig pig=acos(-1.) -c good2=1!.true. chi2_nt(ntr) = sngl(chi2) - nstep_nt(ntr) = 0!nstep + nstep_nt(ntr) = nstep + + phi = al(4) + sinth = al(3) + if(sinth.lt.0)then + sinth = -sinth + phi = phi + pig + endif + npig = aint(phi/(2*pig)) + phi = phi - npig*2*pig + if(phi.lt.0) + $ phi = phi + 2*pig + al(4) = phi + al(3) = sinth - phi = al(4) !(4) - sinth = al(3) !(4) - if(sinth.lt.0)then !(4) - sinth = -sinth !(4) - phi = phi + pig !(4) - endif !(4) - npig = aint(phi/(2*pig)) !(4) - phi = phi - npig*2*pig !(4) - if(phi.lt.0) !(4) - $ phi = phi + 2*pig !(4) - al(4) = phi !(4) - al(3) = sinth !(4) -***************************************************** do i=1,5 al_nt(i,ntr) = sngl(al(i)) do j=1,5 coval(i,j,ntr) = sngl(cov(i,j)) enddo -c print*,al_nt(i,ntr) enddo do ip=1,nplanes ! loop on planes @@ -3981,7 +3986,6 @@ zv_nt(ip,ntr) = sngl(zv(ip)) axv_nt(ip,ntr) = sngl(axv(ip)) ayv_nt(ip,ntr) = sngl(ayv(ip)) -c dedxp(ip,ntr) = sngl(dedxtrk(ip)) !(1) dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)) !(2) dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)) !(2) @@ -3998,22 +4002,11 @@ endif enddo -c call CalcBdL(100,xxxx,IFAIL) -c if(ifps(xxxx).eq.1)BdL(ntr) = xxxx -c$$$ print*,'xgood(ip,ntr) ',(xgood_nt(ip,ntr),ip=1,6) -c$$$ print*,'ygood(ip,ntr) ',(ygood_nt(ip,ntr),ip=1,6) -c$$$ print*,'dedx_x(ip,ntr) ',(dedx_x(ip,ntr),ip=1,6) -c$$$ print*,'dedx_y(ip,ntr) ',(dedx_y(ip,ntr),ip=1,6) end subroutine fill_level2_siglets -c***************************************************** -c 07/10/2005 created by elena vannuccini -c 31/01/2006 modified by elena vannuccini -* to convert adc to mip --> (2) -c***************************************************** * ------------------------------------------------------- * This routine fills the elements of the variables @@ -4022,10 +4015,11 @@ * ------------------------------------------------------- include 'commontracker.f' - include 'level1.f' - include 'level2.f' +c include 'level1.f' include 'calib.f' + include 'level1.f' include 'common_momanhough.f' + include 'level2.f' include 'common_xyzPAM.f' * count #cluster per plane not associated to any track @@ -4033,6 +4027,10 @@ nclsx = 0 nclsy = 0 + do iv = 1,nviews + if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv) + enddo + do icl=1,nclstr1 if(cl_used(icl).eq.0)then !cluster not included in any track ip=nplanes-npl(VIEW(icl))+1 @@ -4076,8 +4074,81 @@ enddo end +*************************************************** +* * +* * +* * +* * +* * +* * +************************************************** + subroutine fill_hough +* ------------------------------------------------------- +* This routine fills the variables related to the hough +* transform, for the debig n-tuple +* ------------------------------------------------------- + include 'commontracker.f' + include 'level1.f' + include 'common_momanhough.f' + include 'common_hough.f' + include 'level2.f' - + if(.false. + $ .or.ntrpt.gt.ntrpt_max_nt + $ .or.ndblt.gt.ndblt_max_nt + $ .or.NCLOUDS_XZ.gt.ncloxz_max + $ .or.NCLOUDS_yZ.gt.ncloyz_max + $ )then + ntrpt_nt=0 + ndblt_nt=0 + NCLOUDS_XZ_nt=0 + NCLOUDS_YZ_nt=0 + else + ndblt_nt=ndblt + ntrpt_nt=ntrpt + if(ndblt.ne.0)then + do id=1,ndblt + alfayz1_nt(id)=alfayz1(id) !Y0 + alfayz2_nt(id)=alfayz2(id) !tg theta-yz + enddo + endif + if(ndblt.ne.0)then + do it=1,ntrpt + alfaxz1_nt(it)=alfaxz1(it) !X0 + alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz + alfaxz3_nt(it)=alfaxz3(it) !1/r + enddo + endif + nclouds_yz_nt=nclouds_yz + nclouds_xz_nt=nclouds_xz + if(nclouds_yz.ne.0)then + nnn=0 + do iyz=1,nclouds_yz + ptcloud_yz_nt(iyz)=ptcloud_yz(iyz) + alfayz1_av_nt(iyz)=alfayz1_av(iyz) + alfayz2_av_nt(iyz)=alfayz2_av(iyz) + nnn=nnn+ptcloud_yz(iyz) + enddo + do ipt=1,nnn + db_cloud_nt(ipt)=db_cloud(ipt) + enddo + endif + if(nclouds_xz.ne.0)then + nnn=0 + do ixz=1,nclouds_xz + ptcloud_xz_nt(ixz)=ptcloud_xz(ixz) + alfaxz1_av_nt(ixz)=alfaxz1_av(ixz) + alfaxz2_av_nt(ixz)=alfaxz2_av(ixz) + alfaxz3_av_nt(ixz)=alfaxz3_av(ixz) + nnn=nnn+ptcloud_xz(ixz) + enddo + do ipt=1,nnn + tr_cloud_nt(ipt)=tr_cloud(ipt) + enddo + endif + endif + end +