--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/10/02 13:12:48 1.6 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2014/01/16 15:29:50 1.45 @@ -12,19 +12,47 @@ 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' include 'level2.f' - include 'momanhough_init.f' +c print*,'======================================================' +c$$$ do ic=1,NCLSTR1 +c$$$ if(.false. +c$$$ $ .or.nsatstrips(ic).gt.0 +c$$$c $ .or.nbadstrips(0,ic).gt.0 +c$$$c $ .or.nbadstrips(4,ic).gt.0 +c$$$c $ .or.nbadstrips(3,ic).gt.0 +c$$$ $ .or..false.)then +c$$$ print*,'--- cl-',ic,' ------------------------' +c$$$ istart = INDSTART(IC) +c$$$ istop = TOTCLLENGTH +c$$$ if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 +c$$$ print*,'ADC ',(CLADC(i),i=istart,istop) +c$$$ print*,'s/n ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) +c$$$ print*,'sgnl ',(CLSIGNAL(i),i=istart,istop) +c$$$ print*,'strip ',(i-INDMAX(ic),i=istart,istop) +c$$$ print*,'view ',VIEW(ic) +c$$$ print*,'maxs ',MAXS(ic) +c$$$ print*,'COG4 ',cog(4,ic) +c$$$ ff = fbad_cog(4,ic) +c$$$ print*,'fbad ',ff +c$$$ print*,(CLBAD(i),i=istart,istop) +c$$$ bb=nbadstrips(0,ic) +c$$$ print*,'#BAD (tot)',bb +c$$$ bb=nbadstrips(4,ic) +c$$$ print*,'#BAD (4)',bb +c$$$ bb=nbadstrips(3,ic) +c$$$ print*,'#BAD (3)',bb +c$$$ ss=nsatstrips(ic) +c$$$ print*,'#saturated ',ss +c$$$ endif +c$$$ enddo -c logical DEBUG -c common/dbg/DEBUG - *------------------------------------------------------------------------------- * STEP 1 *------------------------------------------------------------------------------- @@ -44,25 +72,11 @@ *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- -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 - - - + if(ncp_tot.eq.0)goto 880 !go to next event *----------------------------------------------------- *----------------------------------------------------- * HOUGH TRASFORM @@ -91,11 +105,12 @@ *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- -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 + if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event *------------------------------------------------------------------------------- @@ -123,18 +138,89 @@ * $ ,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 + + cutdistyz=cutystart + cutdistxz=cutxstart -c iflag=0 + 878 continue call doub_to_YZcloud(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event + endif +* ------------------------------------------------ +* first try to release the tolerance +* ------------------------------------------------ + 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 + if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz + goto 878 + endif +* ------------------------------------------------ +* hence reduce the minimum number of plane +* ------------------------------------------------ + if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then + nplyz_min=nplyz_min-1 + if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min + goto 878 endif -c iflag=0 + + if(nclouds_yz.eq.0)goto 880 !go to next event + + +ccc 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 - +* ------------------------------------------------ +* first try to release the tolerance +* ------------------------------------------------ + if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then + cutdistxz=cutdistxz+cutxstep + if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz + goto 879 + endif +* ------------------------------------------------ +* hence reduce the minimum number of plane +* ------------------------------------------------ + if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then + nplxz_min=nplxz_min-1 + if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min + goto 879 + endif + + if(nclouds_xz.eq.0)goto 880 !go to next event + + +c$$$ 881 continue +c$$$* if there is at least three planes on the Y view decreases cuts on X view +c$$$ if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and. +c$$$ $ nplxz_min.ne.y_min_start)then +c$$$ nptxz_min=x_min_step +c$$$ nplxz_min=x_min_start-x_min_step +c$$$ goto 879 +c$$$ endif + 880 return end @@ -144,20 +230,19 @@ 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 trackimage(NTRACKSMAX) + real*8 AL_GUESS(5) *------------------------------------------------------------------------------- * STEP 4 (ITERATED until any other physical track isn't found) @@ -184,28 +269,64 @@ * *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- - ntrk=0 !counter of identified physical tracks +ccc ntrk=0 !counter of identified physical tracks -11111 continue !<<<<<<< come here when performing a new search +c11111 continue !<<<<<<< come here when performing a new search + continue !<<<<<<< come here when performing a new search + + if(nclouds_xz.eq.0)goto 880 !go to next event + if(nclouds_yz.eq.0)goto 880 !go to next event 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 + if(ntracks.eq.0)goto 880 !go to next event FIMAGE=.false. !processing best track (not track image) 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 + +* ------------------------------------------------------- +* order track-candidates according to: +* 1st) decreasing n.points +* 2nd) increasing chi**2 +* ------------------------------------------------------- rchi2best=1000000000. + ndofbest=0 do i=1,ntracks - if(RCHI2_STORE(i).lt.rchi2best.and. - $ RCHI2_STORE(i).gt.0)then + ndof=0 + do ii=1,nplanes + ndof=ndof + $ +int(xgood_store(ii,i)) + $ +int(ygood_store(ii,i)) + enddo + if(ndof.gt.ndofbest)then + ibest=i + rchi2best=RCHI2_STORE(i) + ndofbest=ndof + elseif(ndof.eq.ndofbest)then + if(RCHI2_STORE(i).lt.rchi2best.and. + $ RCHI2_STORE(i).gt.0)then ibest=i rchi2best=RCHI2_STORE(i) - endif + ndofbest=ndof + endif + endif enddo + + if(ibest.eq.0)goto 880 !>> no good candidates *------------------------------------------------------------------------------- * The best track candidate (ibest) is selected and a new fitting is performed. @@ -224,8 +345,10 @@ iimage=0 endif if(icand.eq.0)then - print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand - $ ,ibest,iimage + if(VERBOSE.EQ.1)then + print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand + $ ,ibest,iimage + endif return endif @@ -236,24 +359,35 @@ * ********************************************************** * ************************** 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) - if(ifail.ne.0) then - if(DEBUG)then + iprint=0 +c if(DEBUG.EQ.1)iprint=1 + if(VERBOSE.EQ.1)iprint=1 + if(DEBUG.EQ.1)iprint=2 + call mini2(jstep,ifail,iprint) +cc if(ifail.ne.0) then + if(ifail.ne.0.or.CHI2.ne.CHI2) then !new + if(CHI2.ne.CHI2)CHI2=-9999. !new + if(VERBOSE.EQ.1)then print *, - $ '*** MINIMIZATION FAILURE *** (mini_2) ' + $ '*** MINIMIZATION FAILURE *** (after refinement) ' $ ,iev endif - chi2=-chi2 endif - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'----------------------------- improved track coord' 22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5) do ip=1,6 @@ -264,7 +398,7 @@ endif c rchi2=chi2/dble(ndof) - if(DEBUG)then + if(DEBUG.EQ.1)then print*,' ' print*,'****** SELECTED TRACK *************' print*,'# R. chi2 RIG' @@ -280,38 +414,122 @@ * ------------- 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 + +* ----------------------------------------------------- +* first check if the track is ambiguous +* ----------------------------------------------------- +* (modified on august 2007 by ElenaV) + is1=0 + do ip=1,NPLANES + if(SENSOR_STORE(ip,icand).ne.0)then + is1=SENSOR_STORE(ip,icand) + if(ip.eq.6)is1=3-is1 !last plane inverted + endif + enddo + if(is1.eq.0)then + if(WARNING.EQ.1)print*,'** WARNING ** sensor=0' + goto 122 !jump + endif + do ip=1,NPLANES + is2 = SENSOR_STORE(ip,icand) !sensor + if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted + if( + $ (is1.ne.is2.and.is2.ne.0) + $ )goto 122 !jump (this track cannot have an image) + enddo + if(DEBUG.eq.1)print*,' >>> ambiguous track! ' +* --------------------------------------------------------------- +* take the candidate with the greatest number of matching couples +* if more than one satisfies the requirement, +* choose the one with more points and lower chi2 +* --------------------------------------------------------------- +* count the number of matching couples + ncpmax = 0 + iimage = 0 !id of candidate with better matching do i=1,ntracks - iimage=i + ncp=0 do ip=1,nplanes - if( CP_STORE(nplanes-ip+1,icand).ne. - $ -1*CP_STORE(nplanes-ip+1,i) )iimage=0 + if(CP_STORE(nplanes-ip+1,icand).ne.0)then + if( + $ CP_STORE(nplanes-ip+1,i).ne.0 + $ .and. + $ CP_STORE(nplanes-ip+1,icand).eq. + $ -1*CP_STORE(nplanes-ip+1,i) + $ )then + ncp=ncp+1 !ok + elseif( + $ CP_STORE(nplanes-ip+1,i).ne.0 + $ .and. + $ CP_STORE(nplanes-ip+1,icand).ne. + $ -1*CP_STORE(nplanes-ip+1,i) + $ )then + ncp = 0 + goto 117 !it is not an image candidate + else + + endif + endif 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 + 117 continue + trackimage(i)=ncp !number of matching couples + if(ncp>ncpmax)then + ncpmax=ncp + iimage=i + endif + enddo +* check if there are more than one image candidates + ngood=0 + do i=1,ntracks + if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1 + enddo + if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood +* if there are, choose the best one + if(ngood.gt.0)then +* ------------------------------------------------------- +* order track-candidates according to: +* 1st) decreasing n.points +* 2nd) increasing chi**2 +* ------------------------------------------------------- + rchi2best=1000000000. + ndofbest=0 + do i=1,ntracks + if( trackimage(i).eq.ncpmax )then + ndof=0 + do ii=1,nplanes + ndof=ndof + $ +int(xgood_store(ii,i)) + $ +int(ygood_store(ii,i)) + enddo + if(ndof.gt.ndofbest)then + iimage=i + rchi2best=RCHI2_STORE(i) + ndofbest=ndof + elseif(ndof.eq.ndofbest)then + if(RCHI2_STORE(i).lt.rchi2best.and. + $ RCHI2_STORE(i).gt.0)then + iimage=i + rchi2best=RCHI2_STORE(i) + ndofbest=ndof + endif + endif + endif + enddo + + if(DEBUG.EQ.1)then + print*,'Track candidate ',iimage $ ,' >>> TRACK IMAGE >>> of' $ ,ibest - goto 122 !image track found endif - enddo + + endif + + 122 continue -* --- 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) +* --- and store the results -------------------------------- if(ntrk.eq.NTRKMAX)then - if(verbose) + if(verbose.eq.1) $ print*, $ '** warning ** number of identified '// $ 'tracks exceeds vector dimension ' @@ -319,220 +537,34 @@ cc good2=.false. goto 880 !fill ntp and go to next event endif + + 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$$$ if(ntrk.eq.NTRKMAX)then +c$$$ if(verbose.eq.1) +c$$$ $ print*, +c$$$ $ '** warning ** number of identified '// +c$$$ $ 'tracks exceeds vector dimension ' +c$$$ $ ,'( ',NTRKMAX,' )' +c$$$cc good2=.false. +c$$$ goto 880 !fill ntp and go to next event +c$$$ 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 - -* ********************************************************** -* 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 'commontracker.f' -c$$$ include '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 'commontracker.f' -c$$$ include '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 'commontracker.f' -c$$$ include '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$$$ - ************************************************************ ************************************************************ @@ -557,6 +589,8 @@ * PFAy - Position Finding Algorithm in y (COG2,ETA2,...) * angx - Projected angle in x * angy - Projected angle in y +* bfx - x-component of magnetci field +* bfy - y-component of magnetic field * * --------- COUPLES ------------------------------------------------------- * The couple defines a point in the space. @@ -595,178 +629,144 @@ * * - subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy) + subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) -c***************************************************** -c 07/10/2005 modified by elena vannuccini --> (1) -c 01/02/2006 modified by elena vannuccini --> (2) -c 02/02/2006 modified by Elena Vannuccini --> (3) -c (implemented new p.f.a.) -c 03/02/2006 modified by Elena Vannuccini --> (4) -c (implemented variable resolution) -c***************************************************** include 'commontracker.f' - include 'calib.f' include 'level1.f' + include 'calib.f' include 'common_align.f' include 'common_mech.f' include 'common_xyzPAM.f' - include 'common_resxy.f' - -c logical DEBUG -c common/dbg/DEBUG integer icx,icy !X-Y cluster ID integer sensor integer viewx,viewy character*4 PFAx,PFAy !PFA to be used - real angx,angy !X-Y angle + real ax,ay !X-Y geometric angle + real angx,angy !X-Y effective angle + real bfx,bfy !X-Y b-field components real stripx,stripy + double precision xi,yi,zi + double precision xi_A,yi_A,zi_A + double precision xi_B,yi_B,zi_B double precision xrt,yrt,zrt double precision xrt_A,yrt_A,zrt_A double precision xrt_B,yrt_B,zrt_B -c double precision xi,yi,zi -c double precision xi_A,yi_A,zi_A -c double precision xi_B,yi_B,zi_B parameter (ndivx=30) + + resxPAM = 0 resyPAM = 0 - xPAM = 0. - yPAM = 0. - zPAM = 0. - xPAM_A = 0. - yPAM_A = 0. - zPAM_A = 0. - xPAM_B = 0. - yPAM_B = 0. - zPAM_B = 0. + xPAM = 0.D0 + yPAM = 0.D0 + zPAM = 0.D0 + xPAM_A = 0.D0 + yPAM_A = 0.D0 + zPAM_A = 0.D0 + xPAM_B = 0.D0 + yPAM_B = 0.D0 + zPAM_B = 0.D0 + + if(sensor.lt.1.or.sensor.gt.2)then + print*,'xyz_PAM ***ERROR*** wrong input ' + print*,'sensor ',sensor + icx=0 + icy=0 + endif * ----------------- * CLUSTER X -* ----------------- - +* ----------------- if(icx.ne.0)then - viewx = VIEW(icx) - nldx = nld(MAXS(icx),VIEW(icx)) - nplx = npl(VIEW(icx)) - resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! - - stripx = float(MAXS(icx)) - if(PFAx.eq.'COG1')then !(1) - stripx = stripx !(1) - resxPAM = resxPAM !(1) - elseif(PFAx.eq.'COG2')then - stripx = stripx + cog(2,icx) - 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 stripx = stripx + etacorr - stripx = stripx + pfa_eta2(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) - 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) - 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) - 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) -c resxPAM = resxPAM*fbad_cog(2,icx) !(3)TEMPORANEO - resxPAM = resxPAM*fbad_eta(icx,angx) !(3)(4) - elseif(PFAx.eq.'COG')then !(2) - stripx = stripx + cog(0,icx) !(2) - resxPAM = risx_cog(angx) ! (4) - resxPAM = resxPAM*fbad_cog(0,icx)!(2) - else - print*,'*** Non valid p.f.a. (x) --> ',PFAx + + viewx = VIEW(icx) + nldx = nld(MAXS(icx),VIEW(icx)) + nplx = npl(VIEW(icx)) +c resxPAM = RESXAV + stripx = float(MAXS(icx)) + + if( + $ viewx.lt.1.or. + $ viewx.gt.12.or. + $ nldx.lt.1.or. + $ nldx.gt.3.or. + $ stripx.lt.1.or. + $ stripx.gt.3072.or. + $ .false.)then + print*,'xyz_PAM ***ERROR*** wrong input ' + print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx + icx = 0 + goto 10 endif +* -------------------------- +* magnetic-field corrections +* -------------------------- + stripx = stripx + fieldcorr(viewx,bfy) + angx = effectiveangle(ax,viewx,bfy) + + call applypfa(PFAx,icx,angx,corr,res) + stripx = stripx + corr + resxPAM = res + + 10 continue endif -c if(icy.eq.0.and.icx.ne.0) -c $ print*,PFAx,icx,angx,stripx,resxPAM,'***' - + * ----------------- * CLUSTER Y * ----------------- if(icy.ne.0)then + viewy = VIEW(icy) nldy = nld(MAXS(icy),VIEW(icy)) nply = npl(VIEW(icy)) - resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! +c resyPAM = RESYAV + stripy = float(MAXS(icy)) + if( + $ viewy.lt.1.or. + $ viewy.gt.12.or. + $ nldy.lt.1.or. + $ nldy.gt.3.or. + $ stripy.lt.1.or. + $ stripy.gt.3072.or. + $ .false.)then + print*,'xyz_PAM ***ERROR*** wrong input ' + print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy + icy = 0 + goto 20 + endif if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then - print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' - $ ,icx,icy + if(DEBUG.EQ.1) then + print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' + $ ,icx,icy + endif goto 100 endif - - stripy = float(MAXS(icy)) - if(PFAy.eq.'COG1')then !(1) - stripy = stripy !(1) - resyPAM = resyPAM !(1) - elseif(PFAy.eq.'COG2')then - stripy = stripy + cog(2,icy) - 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 stripy = stripy + etacorr - stripy = stripy + pfa_eta2(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) - 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) - 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) - resyPAM = ris_eta(icy,angy) ! (4) -c resyPAM = resyPAM*fbad_cog(2,icy) !(3)TEMPORANEO - resyPAM = resyPAM*fbad_eta(icy,angy) ! (4) - if(DEBUG.and.fbad_cog(2,icy).ne.1) !(3) - $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3) - elseif(PFAy.eq.'COG')then - stripy = stripy + cog(0,icy) - resyPAM = risy_cog(angy) ! (4) -c resyPAM = ris_eta(icy,angy) ! (4) - resyPAM = resyPAM*fbad_cog(0,icy) - else - print*,'*** Non valid p.f.a. (x) --> ',PFAx - endif +* -------------------------- +* magnetic-field corrections +* -------------------------- + stripy = stripy + fieldcorr(viewy,bfx) + angy = effectiveangle(ay,viewy,bfx) + + call applypfa(PFAy,icy,angy,corr,res) + stripy = stripy + corr + resyPAM = res + + 20 continue endif - + c=========================================================== C COUPLE C=========================================================== @@ -777,14 +777,15 @@ c------------------------------------------------------------------------ if(((mod(int(stripx+0.5)-1,1024)+1).le.3) $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips... - print*,'xyz_PAM (couple):', - $ ' WARNING: false X strip: strip ',stripx + if(DEBUG.EQ.1) then + print*,'xyz_PAM (couple):', + $ ' WARNING: false X strip: strip ',stripx + endif endif - xi = acoordsi(stripx,viewx) - yi = acoordsi(stripy,viewy) - zi = 0. + xi = dcoordsi(stripx,viewx) + yi = dcoordsi(stripy,viewy) + zi = 0.D0 - c------------------------------------------------------------------------ c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c------------------------------------------------------------------------ @@ -818,13 +819,13 @@ yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4 zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4 - xPAM_A = 0. - yPAM_A = 0. - zPAM_A = 0. - - xPAM_B = 0. - yPAM_B = 0. - zPAM_B = 0. + xPAM_A = 0.D0 + yPAM_A = 0.D0 + zPAM_A = 0.D0 + + xPAM_B = 0.D0 + yPAM_B = 0.D0 + zPAM_B = 0.D0 elseif( $ (icx.ne.0.and.icy.eq.0).or. @@ -844,7 +845,9 @@ nldx = nldy viewx = viewy + 1 - yi = acoordsi(stripy,viewy) + xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center + yi = dcoordsi(stripy,viewy) + zi = 0.D0 xi_A = edgeY_d - SiDimX/2 yi_A = yi @@ -854,8 +857,6 @@ yi_B = yi zi_B = 0. -c print*,'Y-cl ',icy,stripy,' --> ',yi -c print*,xi_A,' <--> ',xi_B elseif(icx.ne.0)then c=========================================================== @@ -866,14 +867,17 @@ nldy = nldx viewy = viewx - 1 -c print*,'X-singlet ',icx,nplx,nldx,viewx,stripx -c if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips... if(((mod(int(stripx+0.5)-1,1024)+1).le.3) $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips... - print*,'xyz_PAM (X-singlet):', - $ ' WARNING: false X strip: strip ',stripx + if(DEBUG.EQ.1) then + print*,'xyz_PAM (X-singlet):', + $ ' WARNING: false X strip: strip ',stripx + endif endif - xi = acoordsi(stripx,viewx) + + xi = dcoordsi(stripx,viewx) + yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center + zi = 0.D0 xi_A = xi yi_A = edgeX_d - SiDimY/2 @@ -889,14 +893,13 @@ yi_B = yi endif -c print*,'X-cl ',icx,stripx,' --> ',xi -c print*,yi_A,' <--> ',yi_B else - - print *,'routine xyz_PAM ---> not properly used !!!' - print *,'icx = ',icx - print *,'icy = ',icy + if(DEBUG.EQ.1) then + print *,'routine xyz_PAM ---> not properly used !!!' + print *,'icx = ',icx + print *,'icy = ',icy + endif goto 100 endif @@ -935,6 +938,24 @@ $ + zi_B $ + dz(nplx,nldx,sensor) + + + xrt = xi + $ - omega(nplx,nldx,sensor)*yi + $ + gamma(nplx,nldx,sensor)*zi + $ + dx(nplx,nldx,sensor) + + yrt = omega(nplx,nldx,sensor)*xi + $ + yi + $ - beta(nplx,nldx,sensor)*zi + $ + dy(nplx,nldx,sensor) + + zrt = -gamma(nplx,nldx,sensor)*xi + $ + beta(nplx,nldx,sensor)*yi + $ + zi + $ + dz(nplx,nldx,sensor) + + c xrt = xi c yrt = yi @@ -945,9 +966,12 @@ c in PAMELA reference system c------------------------------------------------------------------------ - xPAM = 0. - yPAM = 0. - zPAM = 0. + xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4 + yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4 + zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4 +c$$$ xPAM = 0.D0 +c$$$ yPAM = 0.D0 +c$$$ zPAM = 0.D0 xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4 yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4 @@ -958,19 +982,191 @@ zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4 -c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')' else - - print *,'routine xyz_PAM ---> not properly used !!!' - print *,'icx = ',icx - print *,'icy = ',icy - + if(DEBUG.EQ.1) then + print *,'routine xyz_PAM ---> not properly used !!!' + print *,'icx = ',icx + print *,'icy = ',icy + endif endif + + 100 continue end +************************************************************************ +* Call xyz_PAM subroutine with default PFA and fill the mini2 common. +* (done to be called from c/c++) +************************************************************************ + + subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy) + + include 'commontracker.f' + include 'level1.f' + include 'common_mini_2.f' + include 'common_xyzPAM.f' + include 'common_mech.f' + include 'calib.f' + +* flag to chose PFA +c$$$ character*10 PFA +c$$$ common/FINALPFA/PFA + + integer icx,icy !X-Y cluster ID + integer sensor + character*4 PFAx,PFAy !PFA to be used + real ax,ay !X-Y geometric angle + real bfx,bfy !X-Y b-field components + + ipx=0 + ipy=0 + +c$$$ PFAx = 'COG4'!PFA +c$$$ PFAy = 'COG4'!PFA + + + if(icx.gt.nclstr1.or.icy.gt.nclstr1)then + print*,'xyzpam: ***WARNING*** clusters ',icx,icy + $ ,' do not exists (n.clusters=',nclstr1,')' + icx = -1*icx + icy = -1*icy + return + + endif + + call idtoc(pfaid,PFAx) + call idtoc(pfaid,PFAy) + + + if(icx.ne.0.and.icy.ne.0)then + + ipx=npl(VIEW(icx)) + ipy=npl(VIEW(icy)) + + if( (nplanes-ipx+1).ne.ip )then + print*,'xyzpam: ***WARNING*** cluster ',icx + $ ,' does not belong to plane: ',ip + icx = -1*icx + return + endif + if( (nplanes-ipy+1).ne.ip )then + print*,'xyzpam: ***WARNING*** cluster ',icy + $ ,' does not belong to plane: ',ip + icy = -1*icy + return + endif + + call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) + + xgood(ip) = 1. + ygood(ip) = 1. + resx(ip) = resxPAM + resy(ip) = resyPAM + + xm(ip) = xPAM + ym(ip) = yPAM + zm(ip) = zPAM + xm_A(ip) = 0.D0 + ym_A(ip) = 0.D0 + xm_B(ip) = 0.D0 + ym_B(ip) = 0.D0 + +c zv(ip) = zPAM + + elseif(icx.eq.0.and.icy.ne.0)then + + ipy=npl(VIEW(icy)) + if( (nplanes-ipy+1).ne.ip )then + print*,'xyzpam: ***WARNING*** cluster ',icy + $ ,' does not belong to plane: ',ip + icy = -1*icy + return + endif + + call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) + + xgood(ip) = 0. + ygood(ip) = 1. + resx(ip) = 1000. + resy(ip) = resyPAM + +c$$$ xm(ip) = -100. +c$$$ ym(ip) = -100. +c$$$ zm(ip) = (zPAM_A+zPAM_B)/2. + xm(ip) = xPAM + ym(ip) = yPAM + zm(ip) = zPAM + xm_A(ip) = xPAM_A + ym_A(ip) = yPAM_A + xm_B(ip) = xPAM_B + ym_B(ip) = yPAM_B + +c zv(ip) = (zPAM_A+zPAM_B)/2. + + elseif(icx.ne.0.and.icy.eq.0)then + + ipx=npl(VIEW(icx)) + + if( (nplanes-ipx+1).ne.ip )then + print*,'xyzpam: ***WARNING*** cluster ',icx + $ ,' does not belong to plane: ',ip + icx = -1*icx + return + endif + + call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) + + xgood(ip) = 1. + ygood(ip) = 0. + resx(ip) = resxPAM + resy(ip) = 1000. + +c$$$ xm(ip) = -100. +c$$$ ym(ip) = -100. +c$$$ zm(ip) = (zPAM_A+zPAM_B)/2. + xm(ip) = xPAM + ym(ip) = yPAM + zm(ip) = zPAM + xm_A(ip) = xPAM_A + ym_A(ip) = yPAM_A + xm_B(ip) = xPAM_B + ym_B(ip) = yPAM_B + +c zv(ip) = (zPAM_A+zPAM_B)/2. + + else + + il = 2 + if(lad.ne.0)il=lad + is = 1 + if(sensor.ne.0)is=sensor + + xgood(ip) = 0. + ygood(ip) = 0. + resx(ip) = 1000. + resy(ip) = 1000. + + xm(ip) = -100. + ym(ip) = -100. + zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 + xm_A(ip) = 0. + ym_A(ip) = 0. + xm_B(ip) = 0. + ym_B(ip) = 0. + +c zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 + + endif + + if(DEBUG.EQ.1)then +22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5) + 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) + endif + end ******************************************************************************** ******************************************************************************** @@ -992,7 +1188,7 @@ * ******************************************************************************** - real function distance_to(XPP,YPP) + real function distance_to(rXPP,rYPP) include 'common_xyzPAM.f' @@ -1001,14 +1197,19 @@ * ( i.e. distance/resolution ) * ----------------------------------- + real rXPP,rYPP + double precision XPP,YPP double precision distance,RE double precision BETA,ALFA,xmi,ymi + XPP=DBLE(rXPP) + YPP=DBLE(rYPP) + * ---------------------- if ( - + xPAM.eq.0.and. - + yPAM.eq.0.and. - + zPAM.eq.0.and. +c$$$ + xPAM.eq.0.and. +c$$$ + yPAM.eq.0.and. +c$$$ + zPAM.eq.0.and. + xPAM_A.ne.0.and. + yPAM_A.ne.0.and. + zPAM_A.ne.0.and. @@ -1046,19 +1247,17 @@ endif distance= - $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2 + $ ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI +cc $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2 distance=dsqrt(distance) -c$$$ print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b -c$$$ $ ,' --- distance_to --- ',xpp,ypp -c$$$ print*,' resolution ',re * ---------------------- elseif( - + xPAM.ne.0.and. - + yPAM.ne.0.and. - + zPAM.ne.0.and. +c$$$ + xPAM.ne.0.and. +c$$$ + yPAM.ne.0.and. +c$$$ + zPAM.ne.0.and. + xPAM_A.eq.0.and. + yPAM_A.eq.0.and. + zPAM_A.eq.0.and. @@ -1071,22 +1270,17 @@ * ---------------------- distance= - $ ((xPAM-XPP)/resxPAM)**2 - $ + - $ ((yPAM-YPP)/resyPAM)**2 + $ ((xPAM-XPP))**2 !QUIQUI + $ + + $ ((yPAM-YPP))**2 +c$$$ $ ((xPAM-XPP)/resxPAM)**2 +c$$$ $ + +c$$$ $ ((yPAM-YPP)/resyPAM)**2 distance=dsqrt(distance) -c$$$ print*,xPAM,yPAM,zPAM -c$$$ $ ,' --- distance_to --- ',xpp,ypp -c$$$ print*,' resolution ',resxPAM,resyPAM else - print* - $ ,' function distance_to ---> wrong usage!!!' - print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM - print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b ' - $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b endif distance_to = sngl(distance) @@ -1126,17 +1320,17 @@ data c1/1.,0.,0.,1./ data c2/1.,-1.,-1.,1./ data c3/1.,1.,0.,0./ - real*8 yvvv,xvvv + double precision yvvv,xvvv double precision xi,yi,zi double precision xrt,yrt,zrt real AA,BB - real yvv(4),xvv(4) + double precision yvv(4),xvv(4) * tollerance to consider the track inside the sensitive area real ptoll data ptoll/150./ !um - external nviewx,nviewy,acoordsi,dcoord + external nviewx,nviewy,dcoordsi,dcoord nplpt = nplPAM !plane viewx = nviewx(nplpt) @@ -1151,15 +1345,9 @@ c------------------------------------------------------------------------ c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame c------------------------------------------------------------------------ - if(((mod(int(stripx+0.5)-1,1024)+1).le.3) - $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips... -c if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips... - print*,'whichsensor: ', - $ ' WARNING: false X strip: strip ',stripx - endif - xi = acoordsi(stripx,viewx) - yi = acoordsi(stripy,viewy) - zi = 0. + xi = dcoordsi(stripx,viewx) + yi = dcoordsi(stripy,viewy) + zi = 0.D0 c------------------------------------------------------------------------ c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c------------------------------------------------------------------------ @@ -1184,8 +1372,6 @@ yvv(iv)=sngl(yvvv) xvv(iv)=sngl(xvvv) -c print*,'LADDER ',il,' SENSOR ',is,' vertexes >> ' -c $ ,iv,xvv(iv),yvv(iv) enddo !end loop on sensor vertexes dtot=0. @@ -1193,8 +1379,8 @@ iv1=iside iv2=mod(iside,4)+1 * straight line passing trhough two consecutive vertexes - AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2)) - BB = yvv(iv1) - AA*xvv(iv1) + AA = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7 + BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7 * point along the straight line closer to the track xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2) yoo = AA*xoo + BB @@ -1206,8 +1392,8 @@ iv1=iside iv2=mod(iside,4)+1 * straight line passing trhough two consecutive vertexes - AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2)) - BB = xvv(iv1) - AA*yvv(iv1) + AA = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7 + BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7 * point along the straight line closer to the track yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2) xoo = AA*yoo + BB @@ -1284,6 +1470,7 @@ * it returns the plane number * include 'commontracker.f' + include 'level1.f' c include 'common_analysis.f' include 'common_momanhough.f' @@ -1309,7 +1496,6 @@ is_cp=0 if(id.lt.0)is_cp=1 if(id.gt.0)is_cp=2 - if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!' return end @@ -1321,6 +1507,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 +1536,7 @@ * positive if sensor =2 * include 'commontracker.f' + include 'level1.f' c include 'calib.f' c include 'level1.f' c include 'common_analysis.f' @@ -1378,274 +1566,6 @@ ************************************************************************* ************************************************************************* ************************************************************************* -c$$$ subroutine book_debug -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'common_momanhough.f' -c$$$ include 'common_level2debug.f' -c$$$ -c$$$ character*35 block1,block2,block3!,block4 -c$$$ $ ,block5!,block6 -c$$$ -c$$$* * * * * * * * * * * * * * * * * * * * * * * * -c$$$* HOUGH TRANSFORM PARAMETERS -c$$$ -c$$$ call HBOOK2(1003 -c$$$ $ ,'y vs tg thyz' -c$$$ $ ,300,-1.,1. !x -c$$$ $ ,3000,-70.,70.,0.) !y -c$$$ -c$$$ call HBOOK1(1004 -c$$$ $ ,'Dy' -c$$$ $ ,100,0.,2.,0.) -c$$$ -c$$$ call HBOOK1(1005 -c$$$ $ ,'D thyz' -c$$$ $ ,100,0.,.05,0.) -c$$$ -c$$$ -c$$$ -c$$$* DEBUG ntuple: -c$$$ call HBNT(ntp_level2+1,'LEVEL2',' ') -c$$$ -c$$$ call HBNAME(ntp_level2+1,'EVENT',good2_nt, -c$$$ $ 'GOOD2:L,NEV2:I') -c$$$ -c$$$ 411 format('NDBLT:I::[0,',I5,']') -c$$$ write(block1,411) ndblt_max_nt -c$$$ call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt, -c$$$ $ block1//' -c$$$ $ ,ALFAYZ1(NDBLT):R -c$$$ $ ,ALFAYZ2(NDBLT):R -c$$$ $ ,DB_CLOUD(NDBLT):I -c$$$ $ ') -c$$$ -c$$$ 412 format('NTRPT:I::[0,',I5,']') -c$$$ write(block2,412) ntrpt_max_nt -c$$$ call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt, -c$$$ $ block2//' -c$$$ $ ,ALFAXZ1(NTRPT):R -c$$$ $ ,ALFAXZ2(NTRPT):R -c$$$ $ ,ALFAXZ3(NTRPT):R -c$$$ $ ,TR_CLOUD(NTRPT):I -c$$$ $ ') -c$$$ -c$$$ -c$$$ 413 format('NCLOUDS_YZ:I::[0,',I4,']') -c$$$c$$$ 414 format('DB_CLOUD(',I4,'):I') -c$$$ write(block3,413) ncloyz_max -c$$$c$$$ write(block4,414) ndblt_max_nt -c$$$ call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ, -c$$$ $ block3//' -c$$$ $ ,ALFAYZ1_AV(NCLOUDS_YZ):R -c$$$ $ ,ALFAYZ2_AV(NCLOUDS_YZ):R -c$$$ $ ,PTCLOUD_YZ(NCLOUDS_YZ):I' -c$$$c$$$ $ ,'//block4 -c$$$ $ ) -c$$$ -c$$$ 415 format('NCLOUDS_XZ:I::[0,',I4,']') -c$$$c$$$ 416 format('TR_CLOUD(',I5,'):I') -c$$$ write(block5,415) ncloxz_max -c$$$c$$$ write(block6,416) ntrpt_max_nt -c$$$ call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ, -c$$$ $ block5//' -c$$$ $ ,ALFAXZ1_AV(NCLOUDS_XZ):R -c$$$ $ ,ALFAXZ2_AV(NCLOUDS_XZ):R -c$$$ $ ,ALFAXZ3_AV(NCLOUDS_XZ):R -c$$$ $ ,PTCLOUD_XZ(NCLOUDS_XZ):I' -c$$$c$$$ $ ,'//block6 -c$$$ $ ) -c$$$ -c$$$ -c$$$ return -c$$$ end -***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** -* -* -* -* -* -* -* -* -* -***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** -c$$$ subroutine book_level2 -c$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2) -c$$$c***************************************************** -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'common_momanhough.f' -c$$$ include 'level2.f' -c$$$ -c$$$ character*35 block1,block2 -c$$$ -c$$$c print*,'__________ booking LEVEL2 n-tuple __________' -c$$$ -c$$$* LEVEL1 ntuple: -c$$$ call HBNT(ntp_level2,'LEVEL2',' ') -c$$$ -c$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$c call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I') -c$$$cccccc 06/10/2005 modified by elena vannuccini -c$$$c call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I -c$$$c $ ,WHIC_CALIB:I,SWCODE:I') -c$$$ call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I -c$$$ $ ,WHICH_CALIB:I,SWCODE:I,CRC(12):L') -c$$$c********************************************************* -c$$$ -c$$$ -c$$$c$$$# ifndef TEST2003 -c$$$c$$$ -c$$$c$$$ call HBNAME(ntp_level2,'CPU',pkt_type -c$$$c$$$ $ ,'PKT_TYPE:I::[0,50] -c$$$c$$$ $ ,PKT_NUM:I -c$$$c$$$ $ ,OBT:I'// -c$$$c$$$c******************************************************** -c$$$c$$$cccccc 11/9/2005 modified by david fedele -c$$$c$$$c $ ,WHICH_CALIB:I::[0,50]') -c$$$c$$$ $ ',CPU_CRC:L') -c$$$c$$$c******************************************************** -c$$$c$$$ -c$$$c$$$# endif -c$$$ -c$$$ 417 format('NTRK:I::[0,',I4,']') -c$$$ 418 format(',IMAGE(NTRK):I::[0,',I4,']') -c$$$ write(block1,417)NTRKMAX -c$$$ write(block2,418)NTRKMAX -c$$$ call HBNAME(ntp_level2,'TRACKS',NTRK, -c$$$ $ block1// -c$$$ $ block2//' -c$$$ $ ,XM(6,NTRK):R -c$$$ $ ,YM(6,NTRK):R -c$$$ $ ,ZM(6,NTRK):R -c$$$ $ ,RESX(6,NTRK):R -c$$$ $ ,RESY(6,NTRK):R -c$$$ $ ,AL(5,NTRK):R -c$$$ $ ,COVAL(5,5,NTRK):R -c$$$ $ ,CHI2(NTRK):R -c$$$ $ ,XGOOD(6,NTRK):I::[0,1] -c$$$ $ ,YGOOD(6,NTRK):I::[0,1] -c$$$ $ ,XV(6,NTRK):R -c$$$ $ ,YV(6,NTRK):R -c$$$ $ ,ZV(6,NTRK):R -c$$$ $ ,AXV(6,NTRK):R -c$$$ $ ,AYV(6,NTRK):R'// -c$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$c $ ,DEDXP(6,NTRK):R'// -c$$$c $ ') -c$$$ $ ',DEDX_X(6,NTRK):R -c$$$ $ ,DEDX_Y(6,NTRK):R'// -c$$$c**************************************************** -c$$$cccccc 06/10/2005 modified by elena vannuccini -c$$$c $ ,CRC(12):L -c$$$ $ ',BdL(NTRK):R' -c$$$ $ ) -c$$$c**************************************************** -c$$$ -c$$$ -c$$$ call HBNAME(ntp_level2,'SINGLETX',nclsx, -c$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$c $ 'NCLSX(6):I,NCLSY(6):I') -c$$$ $ 'NCLSX:I::[0,500],PLANEX(NCLSX):I -c$$$ $ ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2) -c$$$c $ ,XS(NCLSX):R,SGNLXS(NCLSX):R') !(2) -c$$$ call HBNAME(ntp_level2,'SINGLETY',nclsy, -c$$$ $ 'NCLSY:I::[0,500],PLANEY(NCLSY):I -c$$$ $ ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2) -c$$$c $ ,YS(NCLSY):R,SGNLYS(NCLSY):R') !(2) -c$$$ return -c$$$ end - -c$$$ subroutine fill_level2_clouds -c$$$c***************************************************** -c$$$c 29/11/2005 created by elena vannuccini -c$$$c***************************************************** -c$$$ -c$$$* ------------------------------------------------------- -c$$$* This routine fills the variables related to the hough -c$$$* transform, for the debig n-tuple -c$$$* ------------------------------------------------------- -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'common_momanhough.f' -c$$$ include 'common_level2debug.f' -c$$$ include 'level2.f' -c$$$ -c$$$ good2_nt=.true.!good2 -c$$$c nev2_nt=nev2 -c$$$ -c$$$ if(.false. -c$$$ $ .or.ntrpt.gt.ntrpt_max_nt -c$$$ $ .or.ndblt.gt.ndblt_max_nt -c$$$ $ .or.NCLOUDS_XZ.gt.ncloxz_max -c$$$ $ .or.NCLOUDS_yZ.gt.ncloyz_max -c$$$ $ )then -c$$$ good2_nt=.false. -c$$$ ntrpt_nt=0 -c$$$ ndblt_nt=0 -c$$$ NCLOUDS_XZ_nt=0 -c$$$ NCLOUDS_YZ_nt=0 -c$$$ else -c$$$ ndblt_nt=ndblt -c$$$ ntrpt_nt=ntrpt -c$$$ if(ndblt.ne.0)then -c$$$ do id=1,ndblt -c$$$ alfayz1_nt(id)=alfayz1(id) !Y0 -c$$$ alfayz2_nt(id)=alfayz2(id) !tg theta-yz -c$$$c db_cloud_nt(id)=db_cloud(id) -c$$$ enddo -c$$$ endif -c$$$ if(ndblt.ne.0)then -c$$$ do it=1,ntrpt -c$$$ alfaxz1_nt(it)=alfaxz1(it) !X0 -c$$$ alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz -c$$$ alfaxz3_nt(it)=alfaxz3(it) !1/r -c$$$c tr_cloud_nt(it)=tr_cloud(it) -c$$$ enddo -c$$$ endif -c$$$ nclouds_yz_nt=nclouds_yz -c$$$ nclouds_xz_nt=nclouds_xz -c$$$ if(nclouds_yz.ne.0)then -c$$$ nnn=0 -c$$$ do iyz=1,nclouds_yz -c$$$ ptcloud_yz_nt(iyz)=ptcloud_yz(iyz) -c$$$ alfayz1_av_nt(iyz)=alfayz1_av(iyz) -c$$$ alfayz2_av_nt(iyz)=alfayz2_av(iyz) -c$$$ nnn=nnn+ptcloud_yz(iyz) -c$$$ enddo -c$$$ do ipt=1,nnn -c$$$ db_cloud_nt(ipt)=db_cloud(ipt) -c$$$ enddo -c$$$c print*,'#### ntupla #### ptcloud_yz ' -c$$$c $ ,(ptcloud_yz(i),i=1,nclouds_yz) -c$$$c print*,'#### ntupla #### db_cloud ' -c$$$c $ ,(db_cloud(i),i=1,nnn) -c$$$ endif -c$$$ if(nclouds_xz.ne.0)then -c$$$ nnn=0 -c$$$ do ixz=1,nclouds_xz -c$$$ ptcloud_xz_nt(ixz)=ptcloud_xz(ixz) -c$$$ alfaxz1_av_nt(ixz)=alfaxz1_av(ixz) -c$$$ alfaxz2_av_nt(ixz)=alfaxz2_av(ixz) -c$$$ alfaxz3_av_nt(ixz)=alfaxz3_av(ixz) -c$$$ nnn=nnn+ptcloud_xz(ixz) -c$$$ enddo -c$$$ do ipt=1,nnn -c$$$ tr_cloud_nt(ipt)=tr_cloud(ipt) -c$$$ enddo -c$$$c print*,'#### ntupla #### ptcloud_xz ' -c$$$c $ ,(ptcloud_xz(i),i=1,nclouds_xz) -c$$$c print*,'#### ntupla #### tr_cloud ' -c$$$c $ ,(tr_cloud(i),i=1,nnn) -c$$$ endif -c$$$ endif -c$$$ end *************************************************** @@ -1660,13 +1580,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 * -------------- @@ -1677,8 +1595,12 @@ integer badseed,badclx,badcly + iflag = iflag + if(DEBUG.EQ.1)print*,'cl_to_couples:' + +cc if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80 + * init variables - ncp_tot=0 do ip=1,nplanes do ico=1,ncouplemax clx(ip,ico)=0 @@ -1691,19 +1613,45 @@ ncls(ip)=0 enddo do icl=1,nclstrmax_level2 - cl_single(icl)=1 - cl_good(icl)=0 + cl_single(icl) = 1 !all are single + cl_good(icl) = 0 !all are bad + 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. nclusterlimit)then +c mask_view(iv) = 1 + mask_view(iv) = mask_view(iv) + 2**0 + if(DEBUG.EQ.1) + $ print*,' * WARNING * cl_to_couple: n.clusters > ' + $ ,nclusterlimit,' 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 + if(cl_used(icx).ne.0)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 + if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then cl_single(icx)=0 goto 10 endif @@ -1744,11 +1692,18 @@ do icy=1,nclstr1 !loop on cluster (Y) if(mod(VIEW(icy),2).eq.0)goto 20 + + if(cl_used(icx).ne.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 + if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then cl_single(icy)=0 goto 20 endif @@ -1794,24 +1749,21 @@ * charge correlation * (modified to be applied only below saturation... obviously) -* ------------------------------------------------------------- -* >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<< -* ------------------------------------------------------------- - if( .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx) + if( .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx) $ .and. - $ .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx) + $ .not.(sgnl(icy).lt.chmipy.and.sgnl(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=(sgnl(icy) + $ -kch(nplx,nldx)*sgnl(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) + sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx) $ -kch(nplx,nldx)*cch(nplx,nldx)) sss=sss/sqrt(kch(nplx,nldx)**2+1) cut = chcut * (16 + sss/50.) @@ -1820,48 +1772,36 @@ goto 20 !charge not consistent endif endif - -* ------------------> COUPLE <------------------ -* 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 - - if(ncp_plane(nplx).eq.ncouplemax)then - if(verbose)print*, + + if(ncp_plane(nplx).gt.ncouplemax)then + if(verbose.eq.1)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 + $ ,'( ',ncouplemax,' ) --> masked!' +c mask_view(nviewx(nplx)) = 2 +c mask_view(nviewy(nply)) = 2 + mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1 + mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1 + goto 10 endif +* ------------------> COUPLE <------------------ 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 * ---------------------------------------------- + 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)) @@ -1869,254 +1809,89 @@ cls(ip,ncls(ip))=icl endif enddo + +c 80 continue + continue - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'clusters ',nclstr1 print*,'good ',(cl_good(i),i=1,nclstr1) - print*,'singles ',(cl_single(i),i=1,nclstr1) + print*,'used ',(cl_used(i),i=1,nclstr1) + print*,'singlets',(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(verbose)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 - -*************************************************** -* * -* * -* * -* * -* * -* * -************************************************** - 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 + + if(.not.RECOVER_SINGLETS)goto 81 - integer badseed,badcl +* //////////////////////////////////////////////// +* PATCH to recover events with less than 3 couples +* //////////////////////////////////////////////// +* loop over singlet and create "fake" couples +* (with clx=0 or cly=0) +* -* 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 - + if(DEBUG.EQ.1) + $ print*,'>>> Recover singlets ' + $ ,'(creates fake couples) <<<' + do icl=1,nclstr1 + if( + $ cl_single(icl).eq.1.and. + $ cl_used(icl).eq.0.and. + $ .true.)then * ---------------------------------------------------- -* 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 +* jump masked views * ---------------------------------------------------- - - 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 - + if( mask_view(VIEW(icl)).ne.0 ) goto 21 + if(mod(VIEW(icl),2).eq.0)then !=== X views * ---------------------------------------------------- -* 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 +* cut on charge (X VIEW) * ---------------------------------------------------- - - - 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 -* =========================================================== + if(sgnl(icl).lt.dedx_x_min) goto 21 + nplx=npl(VIEW(icl)) +* ------------------> (FAKE) COUPLE <----------------- + ncp_plane(nplx) = ncp_plane(nplx) + 1 + clx(nplx,ncp_plane(nplx))=icl + cly(nplx,ncp_plane(nplx))=0 +c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!! +* ---------------------------------------------------- -* ------------------> COUPLE <------------------ -* 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 + else !=== Y views +* ---------------------------------------------------- +* cut on charge (Y VIEW) +* ---------------------------------------------------- + if(sgnl(icl).lt.dedx_y_min) goto 21 - 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 - iflag=1 - return - endif + nply=npl(VIEW(icl)) +* ------------------> (FAKE) COUPLE <----------------- + ncp_plane(nply) = ncp_plane(nply) + 1 + clx(nply,ncp_plane(nply))=0 + cly(nply,ncp_plane(nply))=icl +c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!! +* ---------------------------------------------------- - 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 -* ---------------------------------------------- + endif + endif !end "single" condition + 21 continue + enddo !end loop over clusters - 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 + if(DEBUG.EQ.1) + $ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) + + + 81 continue - do ip=1,6 - ncp_tot=ncp_tot+ncp_plane(ip) + ncp_tot=0 + do ip=1,NPLANES + 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(verbose)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 + if(DEBUG.EQ.1) + $ print*,'n.couple tot: ',ncp_tot return end - *************************************************** * * @@ -2128,20 +1903,14 @@ ************************************************** subroutine cp_to_doubtrip(iflag) -c***************************************************** -c 02/02/2006 modified by Elena Vannuccini --> (1) -c***************************************************** include 'commontracker.f' + include 'level1.f' include 'common_momanhough.f' - include 'momanhough_init.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'calib.f' - include 'level1.f' -c logical DEBUG -c common/dbg/DEBUG * output flag * -------------- @@ -2173,51 +1942,99 @@ real xc,zc,radius * ----------------------------- + if(DEBUG.EQ.1)print*,'cp_to_doubtrip:' + +* -------------------------------------------- +* put a limit to the maximum number of couples +* per plane, in order to apply hough transform +* (couples recovered during track refinement) +* -------------------------------------------- + do ip=1,nplanes + if(ncp_plane(ip).gt.ncouplelimit)then + mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7 + mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7 + endif + enddo ndblt=0 !number of doublets ntrpt=0 !number of triplets do ip1=1,(nplanes-1) !loop on planes - COPPIA 1 - do is1=1,2 !loop on sensors - COPPIA 1 - +c$$$ print*,'(1) ip ',ip1 +c$$$ $ ,mask_view(nviewx(ip1)) +c$$$ $ ,mask_view(nviewy(ip1)) + if( mask_view(nviewx(ip1)).ne.0 .or. + $ mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane + do is1=1,2 !loop on sensors - COPPIA 1 do icp1=1,ncp_plane(ip1) !loop on COPPIA 1 icx1=clx(ip1,icp1) icy1=cly(ip1,icp1) + +c$$$ print*,'(1) ip ',ip1,' icp ',icp1 + c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1) - call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1) - xm1=xPAM - ym1=yPAM - zm1=zPAM -c print*,'***',is1,xm1,ym1,zm1 +c call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1) + call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.) + xm1=REAL(xPAM) !EM GCC4.7 + ym1=REAL(yPAM) !EM GCC4.7 + zm1=REAL(zPAM) !EM GCC4.7 + do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2 - do is2=1,2 !loop on sensors -ndblt COPPIA 2 - +c$$$ print*,'(2) ip ',ip2 +c$$$ $ ,mask_view(nviewx(ip2)) +c$$$ $ ,mask_view(nviewy(ip2)) + if( mask_view(nviewx(ip2)).ne.0 .or. + $ mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane + do is2=1,2 !loop on sensors -ndblt COPPIA 2 do icp2=1,ncp_plane(ip2) !loop on COPPIA 2 icx2=clx(ip2,icp2) icy2=cly(ip2,icp2) + +c$$$ print*,'(2) ip ',ip2,' icp ',icp2 + c call xyz_PAM c $ (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1) +c call xyz_PAM +c $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM - $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1) - xm2=xPAM - ym2=yPAM - zm2=zPAM - + $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.) + xm2=REAL(xPAM) !EM GCC4.7 + ym2=REAL(yPAM) !EM GCC4.7 + zm2=REAL(zPAM) !EM GCC4.7 + +* --------------------------------------------------- +* both couples must have a y-cluster +* (condition necessary when in RECOVER_SINGLETS mode) +* --------------------------------------------------- + if(icy1.eq.0.or.icy2.eq.0)goto 111 + + if(cl_used(icy1).ne.0)goto 111 + if(cl_used(icy2).ne.0)goto 111 + + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on Y VIEW * (2 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ndblt.eq.ndblt_max)then - if(verbose)print*, + if(verbose.eq.1)print*, $ '** warning ** number of identified '// $ 'doublets exceeds vector dimention ' $ ,'( ',ndblt_max,' )' -c good2=.false. -c goto 880 !fill ntp and go to next event +c good2=.false. +c goto 880 !fill ntp and go to next event + do iv=1,12 +c mask_view(iv) = 3 + mask_view(iv) = mask_view(iv)+ 2**2 + enddo iflag=1 return endif + + +ccc print*,' ',icp1,icp2 + ndblt = ndblt + 1 * store doublet info cpyz1(ndblt)=id_cp(ip1,icp1,is1) @@ -2225,116 +2042,205 @@ * tg(th_yz) alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2) * y0 (cm) - alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1 - + alfayz1(ndblt)=alfayz2(ndblt)*(REAL(zini)-zm1)+ym1 ! EM GCC4.7 zm1, ym1 and alfayz1/2 are REAL + **** -----------------------------------------------**** **** reject non phisical couples **** **** -----------------------------------------------**** + if(SECOND_SEARCH)goto 111 if( $ abs(alfayz2(ndblt)).gt.alfyz2_max $ .or. $ abs(alfayz1(ndblt)).gt.alfyz1_max $ )ndblt = ndblt-1 -c$$$ if(iev.eq.33)then -c$$$ print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2 -c$$$ $ ,' || ',icx1,icy1,icx2,icy2 -c$$$ $ ,' || ',xm1,ym1,xm2,ym2 -c$$$ $ ,' || ',alfayz2(ndblt),alfayz1(ndblt) -c$$$ endif -c$$$ + + 111 continue * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on Y VIEW - end * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples + if(icx1.ne.0)then + if(cl_used(icx1).ne.0)goto 31 + endif + if(icx2.ne.0)then + if(cl_used(icx2).ne.0)goto 31 + endif + + if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples + do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3 +c$$$ print*,'(3) ip ',ip3 +c$$$ $ ,mask_view(nviewx(ip3)) +c$$$ $ ,mask_view(nviewy(ip3)) + if( mask_view(nviewx(ip3)).ne.0 .or. + $ mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane do is3=1,2 !loop on sensors - COPPIA 3 do icp3=1,ncp_plane(ip3) !loop on COPPIA 3 icx3=clx(ip3,icp3) icy3=cly(ip3,icp3) + +c$$$ print*,'(3) ip ',ip3,' icp ',icp3 + +* --------------------------------------------------- +* all three couples must have a x-cluster +* (condition necessary when in RECOVER_SINGLETS mode) +* --------------------------------------------------- + if( + $ icx1.eq.0.or. + $ icx2.eq.0.or. + $ icx3.eq.0.or. + $ .false.)goto 29 + + if(cl_used(icx1).ne.0)goto 29 + if(cl_used(icx2).ne.0)goto 29 + if(cl_used(icx3).ne.0)goto 29 + c call xyz_PAM c $ (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1) +c call xyz_PAM +c $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM - $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1) - xm3=xPAM - ym3=yPAM - zm3=zPAM + $ (icx3,icy3,is3,PFAdef,PFAdef + $ ,0.,0.,0.,0.) + xm3=REAL(xPAM) !EM GCC4.7 + ym3=REAL(yPAM) !EM GCC4.7 + zm3=REAL(zPAM) !EM GCC4.7 + + * find the circle passing through the three points + iflag_t = DEBUG call tricircle(3,xp,zp,angp,resp,chi - $ ,xc,zc,radius,iflag) -c print*,xc,zc,radius + $ ,xc,zc,radius,iflag_t) * the circle must intersect the reference plane - if( -c $ (xc.le.-1.*xclimit.or. -c $ xc.ge.xclimit).and. - $ radius**2.ge.(ZINI-zc)**2.and. - $ iflag.eq.0.and. - $ .true.)then - +cc if(iflag.ne.0)goto 29 + if(iflag_t.ne.0)then +* if tricircle fails, evaluate a straight line + if(DEBUG.eq.1) + $ print*,'TRICIRCLE failure' + $ ,' >>> straight line' + radius=0. + xc=0. + yc=0. + + SZZ=0. + SZX=0. + SSX=0. + SZ=0. + S1=0. + X0=0. + Ax=0. + BX=0. + DO I=1,3 + XX = XP(I) + SZZ=SZZ+ZP(I)*ZP(I) + SZX=SZX+ZP(I)*XX + SSX=SSX+XX + SZ=SZ+ZP(I) + S1=S1+1. + ENDDO + DET=SZZ*S1-SZ*SZ + AX=(SZX*S1-SZ*SSX)/DET + BX=(SZZ*SSX-SZX*SZ)/DET + X0 = AX*REAL(ZINI)+BX ! EM GCC4.7 + + endif + + if( .not.SECOND_SEARCH.and. + $ radius**2.lt.(ZINI-zc)**2)goto 29 + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on X VIEW * (3 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ntrpt.eq.ntrpt_max)then - 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 + if(verbose.eq.1)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 +c mask_view(iv) = 4 + mask_view(iv) = + $ mask_view(iv)+ 2**3 + enddo iflag=1 return endif + +ccc print*,' ',icp1,icp2,icp3 + ntrpt = ntrpt +1 * store triplet info cpxz1(ntrpt)=id_cp(ip1,icp1,is1) cpxz2(ntrpt)=id_cp(ip2,icp2,is2) cpxz3(ntrpt)=id_cp(ip3,icp3,is3) - if(xc.lt.0)then + if(radius.ne.0.and.xc.lt.0)then *************POSITIVE DEFLECTION - alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2) - alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) - alfaxz3(ntrpt) = 1/radius - else + alfaxz1(ntrpt) = xc+sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 + alfaxz2(ntrpt) = (REAL(ZINI)-zc)/ + $ sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 + alfaxz3(ntrpt) = 1/radius + else if(radius.ne.0.and.xc.ge.0)then *************NEGATIVE DEFLECTION - alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2) - alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) - alfaxz3(ntrpt) = -1/radius - endif - + alfaxz1(ntrpt) = xc-sqrt(radius**2-(REAL(ZINI)-zc)**2) + alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/ + $ sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 + alfaxz3(ntrpt) = -1/radius + else if(radius.eq.0)then +*************straight fit + alfaxz1(ntrpt) = X0 + alfaxz2(ntrpt) = AX + alfaxz3(ntrpt) = 0. + endif + +c$$$ print*,'alfaxz1 ', alfaxz1(ntrpt) +c$$$ print*,'alfaxz2 ', alfaxz2(ntrpt) +c$$$ print*,'alfaxz3 ', alfaxz3(ntrpt) + **** -----------------------------------------------**** **** reject non phisical triplets **** **** -----------------------------------------------**** - if( - $ abs(alfaxz2(ntrpt)).gt.alfxz2_max - $ .or. - $ abs(alfaxz1(ntrpt)).gt.alfxz1_max - $ )ntrpt = ntrpt-1 - - -c print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt) + if(SECOND_SEARCH)goto 29 + if( + $ abs(alfaxz2(ntrpt)).gt. + $ alfxz2_max + $ .or. + $ abs(alfaxz1(ntrpt)).gt. + $ alfxz1_max + $ )ntrpt = ntrpt-1 + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on X VIEW - end * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - endif + + 29 continue enddo !end loop on COPPIA 3 enddo !end loop on sensors - COPPIA 3 + 30 continue enddo !end loop on planes - COPPIA 3 - 30 continue - - 1 enddo !end loop on COPPIA 2 + + 31 continue +c 1 enddo !end loop on COPPIA 2 + enddo !end loop on COPPIA 2 enddo !end loop on sensors - COPPIA 2 + 20 continue enddo !end loop on planes - COPPIA 2 - + +c 11 continue + continue enddo !end loop on COPPIA1 enddo !end loop on sensors - COPPIA 1 + 10 continue enddo !end loop on planes - COPPIA 1 - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'--- doublets ',ndblt print*,'--- triplets ',ntrpt endif @@ -2359,11 +2265,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 * -------------- @@ -2382,6 +2287,7 @@ integer cp_useds1(ncouplemaxtot) ! sensor 1 integer cp_useds2(ncouplemaxtot) ! sensor 2 + if(DEBUG.EQ.1)print*,'doub_to_YZcloud:' *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * classification of DOUBLETS @@ -2395,12 +2301,11 @@ 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 -c print*,'--------------' -c print*,'** ',idb1,' **' - do icp=1,ncp_tot cp_useds1(icp)=0 !init cp_useds2(icp)=0 !init @@ -2436,19 +2341,11 @@ * doublet distance in parameter space distance= $ ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2 - $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2 + $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2 distance = sqrt(distance) -c$$$ if(iev.eq.33)then -c$$$ if(distance.lt.100) -c$$$ $ print*,'********* ',idb1,idbref,idb2,distance -c$$$ if(distance.lt.100) -c$$$ $ print*,'********* ',alfayz1(idbref),alfayz1(idb2) -c$$$ $ ,alfayz2(idbref),alfayz2(idb2) -c$$$ endif if(distance.lt.cutdistyz)then -c print*,idb1,idb2,distance,' cloud ',nclouds_yz if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1 if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1 if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1 @@ -2464,13 +2361,13 @@ temp1 = temp1 + alfayz1(idb2) temp2 = temp2 + alfayz2(idb2) -c print*,'* idbref,idb2 ',idbref,idb2 endif 1118 continue enddo !end loop (2) on DOUBLETS - 1188 continue +c 1188 continue + continue enddo !end loop on... bo? nptloop=npv @@ -2487,7 +2384,9 @@ enddo ncpused=0 do icp=1,ncp_tot - if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then + if( + $ (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and. + $ .true.)then ncpused=ncpused+1 ip=ip_cp(icp) hit_plane(ip)=1 @@ -2497,21 +2396,23 @@ do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo -c print*,'>>>> ',ncpused,npt,nplused - 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 * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< if(nclouds_yz.ge.ncloyz_max)then - if(verbose)print*, + if(verbose.eq.1)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 +c mask_view(iv) = 5 + mask_view(iv) = mask_view(iv) + 2**4 + enddo iflag=1 return endif @@ -2527,21 +2428,17 @@ c ptcloud_yz_nt(nclouds_yz)=npt do ipt=1,npt db_cloud(npt_tot+ipt) = db_all(ipt) -c print*,'>> ',ipt,db_all(ipt) enddo npt_tot=npt_tot+npt - if(DEBUG)then - print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' + if(DEBUG.EQ.1)then print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points' - print*,'- alfayz1 ',alfayz1_av(nclouds_yz) - print*,'- alfayz2 ',alfayz2_av(nclouds_yz) - print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) - print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) - print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) -c$$$ print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = ' -c$$$ $ ,ptcloud_yz(nclouds_yz) -c$$$ print*,'nt-uple: db_cloud(...) = ' -c$$$ $ ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot) + print*,'- alfayz1 ',alfayz1_av(nclouds_yz) + print*,'- alfayz2 ',alfayz2_av(nclouds_yz) + print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) + print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) + print*,'cpcloud_yz ' + $ ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot) + print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) endif * >>> NEW CLOUD <<< * ~~~~~~~~~~~~~~~~~ @@ -2549,10 +2446,14 @@ enddo !end loop (1) on DOUBLETS - if(DEBUG)then - print*,'---------------------- ' + if(nloop.lt.nstepy)then + cutdistyz = cutdistyz+cutystep + nloop = nloop+1 + goto 90 + endif + + if(DEBUG.EQ.1)then print*,'Y-Z total clouds ',nclouds_yz - print*,' ' endif @@ -2575,11 +2476,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 * -------------- @@ -2599,6 +2499,8 @@ integer cp_useds1(ncouplemaxtot) ! sensor 1 integer cp_useds2(ncouplemaxtot) ! sensor 2 + if(DEBUG.EQ.1)print*,'trip_to_XZcloud:' + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * classification of TRIPLETS * according to distance in parameter space @@ -2610,10 +2512,10 @@ 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*,'--------------' -c print*,'** ',itr1,' **' do icp=1,ncp_tot cp_useds1(icp)=0 @@ -2647,15 +2549,28 @@ do itr2=1,ntrpt !loop (2) on TRIPLETS if(itr2.eq.itr1)goto 11188 !next triplet if(tr_used(itr2).eq.1)goto 11188 !next triplet + + * triplet distance in parameter space * solo i due parametri spaziali per il momemnto distance= $ ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2 - $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2 + $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2 distance = sqrt(distance) - - if(distance.lt.cutdistxz)then -c print*,idb1,idb2,distance,' cloud ',nclouds_yz + + +* ------------------------------------------------------------------------ +* FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE +* ------------------------------------------------------------------------ +* (added in august 2007) + istrimage=0 + if( + $ abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and. + $ abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and. + $ abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and. + $ .true.)istrimage=1 + + if(distance.lt.cutdistxz.or.istrimage.eq.1)then if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1 if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1 if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1 @@ -2674,13 +2589,13 @@ temp1 = temp1 + alfaxz1(itr2) temp2 = temp2 + alfaxz2(itr2) temp3 = temp3 + alfaxz3(itr2) -c print*,'* itrref,itr2 ',itrref,itr2,distance endif 11188 continue enddo !end loop (2) on TRIPLETS -11888 continue +c11888 continue + continue enddo !end loop on... bo? nptloop=npv @@ -2695,13 +2610,14 @@ * 1bis) * 2) it is not already stored * ------------------------------------------ -c print*,'check cp_used' do ip=1,nplanes hit_plane(ip)=0 enddo ncpused=0 do icp=1,ncp_tot - if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then + if( + $ (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and. + $ .true.)then ncpused=ncpused+1 ip=ip_cp(icp) hit_plane(ip)=1 @@ -2711,19 +2627,21 @@ do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo - 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 + if(nplused.lt.nplxz_min)goto 22288 !next triplet * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< if(nclouds_xz.ge.ncloxz_max)then - if(verbose)print*, + if(verbose.eq.1)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 +c mask_view(iv) = 6 + mask_view(iv) = mask_view(iv) + 2**5 + enddo iflag=1 return endif @@ -2741,29 +2659,30 @@ enddo npt_tot=npt_tot+npt - if(DEBUG)then - print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' - print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' - print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz) - print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz) - print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz) - print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) - print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) + if(DEBUG.EQ.1)then + print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' + print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz) + print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz) + print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz) + print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) + print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) + print*,'cpcloud_xz ' + $ ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot) print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) -c$$$ print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = ' -c$$$ $ ,ptcloud_xz(nclouds_xz) -c$$$ print*,'nt-uple: tr_cloud(...) = ' -c$$$ $ ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot) endif * >>> NEW CLOUD <<< * ~~~~~~~~~~~~~~~~~ 22288 continue enddo !end loop (1) on DOUBLETS - - if(DEBUG)then - print*,'---------------------- ' + + if(nloop.lt.nstepx)then + cutdistxz=cutdistxz+cutxstep + nloop=nloop+1 + goto 91 + endif + + if(DEBUG.EQ.1)then print*,'X-Z total clouds ',nclouds_xz - print*,' ' endif @@ -2781,19 +2700,15 @@ ************************************************** subroutine clouds_to_ctrack(iflag) -c***************************************************** -c 02/02/2006 modified by Elena Vannuccini --> (1) -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 logical DEBUG -c common/dbg/DEBUG + * output flag * -------------- @@ -2809,23 +2724,22 @@ * ----------------------------------------------------------- * 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 real fitz(nplanes) !z coordinates of the planes in cm + if(DEBUG.EQ.1)print*,'clouds_to_ctrack:' ntracks=0 !counter of track candidates - do iyz=1,nclouds_yz !loop on YZ couds - do ixz=1,nclouds_xz !loop on XZ couds + do iyz=1,nclouds_yz !loop on YZ clouds + do ixz=1,nclouds_xz !loop on XZ clouds * -------------------------------------------------- * check of consistency of the clouds @@ -2834,27 +2748,72 @@ * of the two clouds * -------------------------------------------------- do ip=1,nplanes - hit_plane(ip)=0 - ncp_match(ip)=0 + hit_plane(ip)=0 !n.matching couples (REAL couples, not SINGLETS) + ncp_match(ip)=0 !n.matching couples per plane do icpp=1,ncouplemax cp_match(ip,icpp)=0 !init couple list enddo enddo - ncp_ok=0 - do icp=1,ncp_tot !loop on couples -* get info on - cpintersec(icp)=min( - $ cpcloud_yz(iyz,icp), - $ cpcloud_xz(ixz,icp)) - if( - $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or. - $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or. - $ .false.)cpintersec(icp)=0 + ncp_ok=0 !count n.matching-couples + ncpx_ok=0 !count n.matching-couples with x cluster + ncpy_ok=0 !count n.matching-couples with y cluster + + + do icp=1,ncp_tot !loop over couples + + if(.not.RECOVER_SINGLETS)then +* ------------------------------------------------------ +* if NOT in RECOVER_SINGLETS mode, take the intersection +* between xz yz clouds +* ------------------------------------------------------ + cpintersec(icp)=min( + $ cpcloud_yz(iyz,icp), + $ cpcloud_xz(ixz,icp)) +* cpintersec is >0 if yz and xz clouds contain the same image of couple icp +* ------------------------------------------------------ +* discard the couple if the sensor is in conflict +* ------------------------------------------------------ + if( + $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or. + $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or. + $ .false.)cpintersec(icp)=0 + else +* ------------------------------------------------------ +* if RECOVER_SINGLETS take the union +* (otherwise the fake couples formed by singlets would be +* discarded) +* ------------------------------------------------------ + cpintersec(icp)=max( + $ cpcloud_yz(iyz,icp), + $ cpcloud_xz(ixz,icp)) +c$$$ if(cpcloud_yz(iyz,icp).gt.0) +c$$$ $ cpintersec(icp)=cpcloud_yz(iyz,icp) +* cpintersec is >0 if either yz or xz cloud contains the couple icp + endif + +c$$$ print*,icp,ip_cp(icp),' -- ',cpintersec(icp) + if(cpintersec(icp).ne.0)then - ncp_ok=ncp_ok+1 ip=ip_cp(icp) hit_plane(ip)=1 +c$$$ if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0) +c$$$ $ ncp_ok=ncp_ok+1 +c$$$ if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0) +c$$$ $ ncpx_ok=ncpx_ok+1 +c$$$ if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0) +c$$$ $ ncpy_ok=ncpy_ok+1 + + if( cpcloud_yz(iyz,icp).gt.0.and. + $ cpcloud_xz(ixz,icp).gt.0) + $ ncp_ok=ncp_ok+1 + if( cpcloud_yz(iyz,icp).gt.0.and. + $ cpcloud_xz(ixz,icp).eq.0) + $ ncpy_ok=ncpy_ok+1 + if( cpcloud_yz(iyz,icp).eq.0.and. + $ cpcloud_xz(ixz,icp).gt.0) + $ ncpx_ok=ncpx_ok+1 + if(cpintersec(icp).eq.1)then * 1) only the couple image in sensor 1 matches id=-icp @@ -2881,52 +2840,33 @@ do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo - - if(nplused.lt.nplxz_min)goto 888 !next doublet - if(ncp_ok.lt.ncpok)goto 888 !next cloud - - if(DEBUG)then - print*,'Combination ',iyz,ixz - $ ,' db ',ptcloud_yz(iyz) - $ ,' tr ',ptcloud_xz(ixz) - $ ,' -----> # matching couples ',ncp_ok - endif -c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~' -c$$$ print*,'Configurazione cluster XZ' -c$$$ print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1)) -c$$$ print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1)) -c$$$ print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1)) -c$$$ print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1)) -c$$$ print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1)) -c$$$ print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1)) -c$$$ print*,'Configurazione cluster YZ' -c$$$ print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1)) -c$$$ print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1)) -c$$$ print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1)) -c$$$ print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1)) -c$$$ print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1)) -c$$$ print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1)) -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) - - if(DEBUG)then + + if(nplused.lt.3)goto 888 !next combination +ccc if(nplused.lt.nplxz_min)goto 888 !next combination +ccc if(nplused.lt.nplyz_min)goto 888 !next combination +* ----------------------------------------------------------- +* if in RECOVER_SINGLET mode, the two clouds must have +* at least ONE intersecting real couple +* ----------------------------------------------------------- + if(ncp_ok.lt.1)goto 888 !next combination + + if(DEBUG.EQ.1)then + print*,'////////////////////////////' + print*,'Cloud combination (Y,X): ',iyz,ixz + print*,' db ',ptcloud_yz(iyz) + print*,' tr ',ptcloud_xz(ixz) + print*,' -----> # matching couples ',ncp_ok + print*,' -----> # fake couples (X)',ncpx_ok + print*,' -----> # fake couples (Y)',ncpy_ok + do icp=1,ncp_tot + print*,'cp ',icp,' >' + $ ,' x ',cpcloud_xz(ixz,icp) + $ ,' y ',cpcloud_yz(iyz,icp) + $ ,' ==> ',cpintersec(icp) + enddo + endif + + if(DEBUG.EQ.1)then print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6)) print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5)) print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4)) @@ -2959,10 +2899,45 @@ hit_plane(6)=icp6 if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6 - + if(DEBUG.eq.1) + $ print*,'combination: ' + $ ,cp_match(6,icp1) + $ ,cp_match(5,icp2) + $ ,cp_match(4,icp3) + $ ,cp_match(3,icp4) + $ ,cp_match(2,icp5) + $ ,cp_match(1,icp6) + + +* --------------------------------------- +* check if this group of couples has been +* already fitted +* --------------------------------------- + do ica=1,ntracks + isthesame=1 + do ip=1,NPLANES + if(hit_plane(ip).ne.0)then + if( CP_STORE(nplanes-ip+1,ica) + $ .ne. + $ cp_match(ip,hit_plane(ip)) ) + $ isthesame=0 + else + if( CP_STORE(nplanes-ip+1,ica) + $ .ne. + $ 0 ) + $ isthesame=0 + endif + enddo + if(isthesame.eq.1)then + if(DEBUG.eq.1) + $ print*,'(already fitted)' + goto 666 !jump to next combination + endif + enddo + call track_init !init TRACK common - do ip=1,nplanes !loop on planes + do ip=1,nplanes !loop on planes (bottom to top) if(hit_plane(ip).ne.0)then id=cp_match(ip,hit_plane(ip)) is=is_cp(id) @@ -2976,68 +2951,117 @@ * ************************* c call xyz_PAM(icx,icy,is, c $ 'COG2','COG2',0.,0.) +c call xyz_PAM(icx,icy,is, !(1) +c $ PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM(icx,icy,is, !(1) - $ PFAdef,PFAdef,0.,0.) !(1) + $ PFAdef,PFAdef,0.,0.,0.,0.) * ************************* * ----------------------------- - xgood(nplanes-ip+1)=1. - ygood(nplanes-ip+1)=1. - xm(nplanes-ip+1)=xPAM - ym(nplanes-ip+1)=yPAM - zm(nplanes-ip+1)=zPAM - resx(nplanes-ip+1)=resxPAM - resy(nplanes-ip+1)=resyPAM + if(icx.gt.0.and.icy.gt.0)then + xgood(nplanes-ip+1)=1. + ygood(nplanes-ip+1)=1. + xm(nplanes-ip+1)=xPAM + ym(nplanes-ip+1)=yPAM + zm(nplanes-ip+1)=zPAM + resx(nplanes-ip+1)=resxPAM + resy(nplanes-ip+1)=resyPAM + if(DEBUG.EQ.1)print*,'(X,Y)' + $ ,nplanes-ip+1,xPAM,yPAM + else + xm_A(nplanes-ip+1) = xPAM_A + ym_A(nplanes-ip+1) = yPAM_A + xm_B(nplanes-ip+1) = xPAM_B + ym_B(nplanes-ip+1) = yPAM_B + zm(nplanes-ip+1) + $ = (zPAM_A+zPAM_B)/2. + resx(nplanes-ip+1) = resxPAM + resy(nplanes-ip+1) = resyPAM + if(icx.eq.0.and.icy.gt.0)then + xgood(nplanes-ip+1)=0. + ygood(nplanes-ip+1)=1. + resx(nplanes-ip+1) = 1000. + if(DEBUG.EQ.1)print*,'( Y)' + $ ,nplanes-ip+1,xPAM,yPAM + elseif(icx.gt.0.and.icy.eq.0)then + xgood(nplanes-ip+1)=1. + ygood(nplanes-ip+1)=0. + if(DEBUG.EQ.1)print*,'(X )' + $ ,nplanes-ip+1,xPAM,yPAM + resy(nplanes-ip+1) = 1000. + else + print*,'both icx=0 and icy=0' + $ ,' ==> IMPOSSIBLE!!' + endif + endif * ----------------------------- endif enddo !end loop on planes * ********************************************************** * ************************** 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.EQ.1)iprint=1 + if(DEBUG.EQ.1)iprint=2 + call mini2(jstep,ifail,iprint) if(ifail.ne.0) then - if(DEBUG)then + if(DEBUG.EQ.1)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 *** * ********************************************************** if(chi2.le.0.)goto 666 + if(chi2.ge.1.e08)goto 666 !OPTIMIZATION + if(chi2.ne.chi2)goto 666 !OPTIMIZATION * -------------------------- * STORE candidate TRACK INFO * -------------------------- if(ntracks.eq.NTRACKSMAX)then - if(verbose)print*, + if(verbose.eq.1)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 +c mask_view(iv) = 7 + mask_view(iv) = mask_view(iv) + 2**6 + enddo iflag=1 return endif ntracks = ntracks + 1 -c$$$ ndof=0 - do ip=1,nplanes -c$$$ ndof=ndof -c$$$ $ +int(xgood(ip)) -c$$$ $ +int(ygood(ip)) + do ip=1,nplanes !top to bottom + XV_STORE(ip,ntracks)=sngl(xv(ip)) YV_STORE(ip,ntracks)=sngl(yv(ip)) - ZV_STORE(ip,ntracks)=sngl(zv(ip)) + ZV_STORE(ip,ntracks)=sngl(zv(ip)) XM_STORE(ip,ntracks)=sngl(xm(ip)) YM_STORE(ip,ntracks)=sngl(ym(ip)) ZM_STORE(ip,ntracks)=sngl(zm(ip)) @@ -3050,23 +3074,39 @@ AYV_STORE(ip,ntracks)=sngl(ayv(ip)) XGOOD_STORE(ip,ntracks)=sngl(xgood(ip)) YGOOD_STORE(ip,ntracks)=sngl(ygood(ip)) +* NB! hit_plane is defined from bottom to top if(hit_plane(ip).ne.0)then CP_STORE(nplanes-ip+1,ntracks)= $ cp_match(ip,hit_plane(ip)) + SENSOR_STORE(nplanes-ip+1,ntracks) + $ = is_cp(cp_match(ip,hit_plane(ip))) + + icl= + $ clx(ip,icp_cp( + $ cp_match(ip,hit_plane(ip) + $ ))); + if(icl.eq.0) + $ icl= + $ cly(ip,icp_cp( + $ cp_match(ip,hit_plane(ip) + $ ))); + + LADDER_STORE(nplanes-ip+1,ntracks) + $ = LADDER(icl); else CP_STORE(nplanes-ip+1,ntracks)=0 + SENSOR_STORE(nplanes-ip+1,ntracks)=0 + LADDER_STORE(nplanes-ip+1,ntracks)=0 endif - CLS_STORE(nplanes-ip+1,ntracks)=0 + BX_STORE(ip,ntracks)=0!I dont need it now + BY_STORE(ip,ntracks)=0!I dont need it now + CLS_STORE(ip,ntracks)=0 do i=1,5 AL_STORE(i,ntracks)=sngl(AL(i)) enddo enddo -c$$$ * Number of Degree Of Freedom -c$$$ ndof=ndof-5 -c$$$ * reduced chi^2 -c$$$ rchi2=chi2/dble(ndof) - RCHI2_STORE(ntracks)=chi2 + RCHI2_STORE(ntracks)=REAL(chi2) * -------------------------------- * STORE candidate TRACK INFO - end @@ -3086,17 +3126,23 @@ if(ntracks.eq.0)then iflag=1 - return +cc return endif - if(DEBUG)then - print*,'****** TRACK CANDIDATES ***********' - print*,'# R. chi2 RIG' - do i=1,ntracks - print*,i,' --- ',rchi2_store(i),' --- ' - $ ,1./abs(AL_STORE(5,i)) - enddo - print*,'***********************************' + if(DEBUG.EQ.1)then + print*,'****** TRACK CANDIDATES *****************' + print*,'# R. chi2 RIG ndof' + do i=1,ntracks + 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) + print*,i,' --- ',rchi2_store(i),' --- ' + $ ,1./abs(AL_STORE(5,i)),' --- ',ndof + enddo + print*,'*****************************************' endif @@ -3115,28 +3161,33 @@ subroutine refine_track(ibest) -c****************************************************** -cccccc 06/10/2005 modified by elena vannuccini ---> (1) -cccccc 31/01/2006 modified by elena vannuccini ---> (2) -cccccc 12/08/2006 modified by elena vannucicni ---> (3) -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' include 'calib.f' -c logical DEBUG -c common/dbg/DEBUG - * flag to chose PFA character*10 PFA common/FINALPFA/PFA + double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7 + double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7 + double precision zmm,zmm_A,zmm_B !EM GCC4.7 + double precision clincnewc !EM GCC4.7 + double precision clincnew !EM GCC4.7 + + real k(6) + DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/ + + real xp,yp,zp + real xyzp(3),bxyz(3) + equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3)) + + if(DEBUG.EQ.1)print*,'refine_track:' * ================================================= * new estimate of positions using ETA algorithm * and @@ -3145,12 +3196,60 @@ call track_init do ip=1,nplanes !loop on planes + if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... ' + + xP=XV_STORE(nplanes-ip+1,ibest) + yP=YV_STORE(nplanes-ip+1,ibest) + zP=ZV_STORE(nplanes-ip+1,ibest) + call gufld(xyzp,bxyz) + BX_STORE(nplanes-ip+1,ibest)=bxyz(1) + BY_STORE(nplanes-ip+1,ibest)=bxyz(2) +c$$$ bxyz(1)=0 +c$$$ bxyz(2)=0 +c$$$ bxyz(3)=0 +* ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * 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. +* ||||||||||||||||||||||||||||||||||||||||||||||||| +c$$$ if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and. +c$$$ $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then +c$$$ +c$$$ id=CP_STORE(nplanes-ip+1,ibest) +c$$$ +c$$$ is=is_cp(id) +c$$$ icp=icp_cp(id) +c$$$ if(ip_cp(id).ne.ip) +c$$$ $ print*,'OKKIO!!' +c$$$ $ ,'id ',id,is,icp +c$$$ $ ,ip_cp(id),ip +c$$$ icx=clx(ip,icp) +c$$$ icy=cly(ip,icp) +c$$$c call xyz_PAM(icx,icy,is, +c$$$c $ PFA,PFA, +c$$$c $ AXV_STORE(nplanes-ip+1,ibest), +c$$$c $ AYV_STORE(nplanes-ip+1,ibest)) +c$$$ call xyz_PAM(icx,icy,is, +c$$$ $ PFA,PFA, +c$$$ $ AXV_STORE(nplanes-ip+1,ibest), +c$$$ $ AYV_STORE(nplanes-ip+1,ibest), +c$$$ $ bxyz(1), +c$$$ $ bxyz(2) +c$$$ $ ) +c$$$ +c$$$ xm(nplanes-ip+1) = xPAM +c$$$ ym(nplanes-ip+1) = yPAM +c$$$ zm(nplanes-ip+1) = zPAM +c$$$ xgood(nplanes-ip+1) = 1 +c$$$ ygood(nplanes-ip+1) = 1 +c$$$ resx(nplanes-ip+1) = resxPAM +c$$$ resy(nplanes-ip+1) = resyPAM +c$$$ +c$$$ dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) +c$$$ dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) + if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or. $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then id=CP_STORE(nplanes-ip+1,ibest) @@ -3163,47 +3262,87 @@ $ ,ip_cp(id),ip icx=clx(ip,icp) icy=cly(ip,icp) +c call xyz_PAM(icx,icy,is, +c $ PFA,PFA, +c $ AXV_STORE(nplanes-ip+1,ibest), +c $ AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(icx,icy,is, -c $ 'ETA2','ETA2', $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest), - $ AYV_STORE(nplanes-ip+1,ibest)) -c$$$ call xyz_PAM(icx,icy,is, -c$$$ $ 'COG2','COG2', -c$$$ $ 0., -c$$$ $ 0.) - xm(nplanes-ip+1) = xPAM - ym(nplanes-ip+1) = yPAM - zm(nplanes-ip+1) = zPAM - xgood(nplanes-ip+1) = 1 - ygood(nplanes-ip+1) = 1 - resx(nplanes-ip+1) = resxPAM - resy(nplanes-ip+1) = resyPAM - -c dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1) - 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) + $ AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) + + if(icx.gt.0.and.icy.gt.0)then + xm(nplanes-ip+1) = xPAM + ym(nplanes-ip+1) = yPAM + zm(nplanes-ip+1) = zPAM + xm_A(nplanes-ip+1) = 0. + ym_A(nplanes-ip+1) = 0. + xm_B(nplanes-ip+1) = 0. + ym_B(nplanes-ip+1) = 0. + xgood(nplanes-ip+1) = 1 + ygood(nplanes-ip+1) = 1 + resx(nplanes-ip+1) = resxPAM + resy(nplanes-ip+1) = resyPAM + dedxtrk_x(nplanes-ip+1)= + $ sgnl(icx)/mip(VIEW(icx),LADDER(icx)) + dedxtrk_y(nplanes-ip+1)= + $ sgnl(icy)/mip(VIEW(icy),LADDER(icy)) + else + xm(nplanes-ip+1) = 0. + ym(nplanes-ip+1) = 0. + zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2. + xm_A(nplanes-ip+1) = xPAM_A + ym_A(nplanes-ip+1) = yPAM_A + xm_B(nplanes-ip+1) = xPAM_B + ym_B(nplanes-ip+1) = yPAM_B + xgood(nplanes-ip+1) = 0 + ygood(nplanes-ip+1) = 0 + resx(nplanes-ip+1) = 1000.!resxPAM + resy(nplanes-ip+1) = 1000.!resyPAM + dedxtrk_x(nplanes-ip+1)= 0 + dedxtrk_y(nplanes-ip+1)= 0 + if(icx.gt.0)then + xgood(nplanes-ip+1) = 1 + resx(nplanes-ip+1) = resxPAM + dedxtrk_x(nplanes-ip+1)= + $ sgnl(icx)/mip(VIEW(icx),LADDER(icx)) + elseif(icy.gt.0)then + ygood(nplanes-ip+1) = 1 + resy(nplanes-ip+1) = resyPAM + dedxtrk_y(nplanes-ip+1)= + $ sgnl(icy)/mip(VIEW(icy),LADDER(icy)) + endif + endif +* ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * If the plane has NOT been already included, * it tries to include a COUPLE or a single cluster * ------------------------------------------------- +* ||||||||||||||||||||||||||||||||||||||||||||||||| else xgood(nplanes-ip+1)=0 ygood(nplanes-ip+1)=0 + + CP_STORE(nplanes-ip+1,ibest)=0 !re-init + CLS_STORE(nplanes-ip+1,ibest)=0 + * -------------------------------------------------------------- * determine which ladder and sensor are intersected by the track - xP=XV_STORE(nplanes-ip+1,ibest) - yP=YV_STORE(nplanes-ip+1,ibest) - zP=ZV_STORE(nplanes-ip+1,ibest) call whichsensor(ip,xP,yP,nldt,ist) * if the track hit the plane in a dead area, go to the next plane if(nldt.eq.0.or.ist.eq.0)goto 133 + + SENSOR_STORE(nplanes-ip+1,IBEST)=ist + LADDER_STORE(nplanes-ip+1,IBEST)=nldt * -------------------------------------------------------------- - if(DEBUG)then + if(DEBUG.EQ.1)then print*, $ '------ Plane ',ip,' intersected on LADDER ',nldt $ ,' SENSOR ',ist @@ -3214,8 +3353,7 @@ * =========================================== * STEP 1 >>>>>>> try to include a new couple * =========================================== -c if(DEBUG)print*,'>>>> try to include a new couple' - distmin=1000000. + distmin=100000000. xmm = 0. ymm = 0. zmm = 0. @@ -3228,23 +3366,28 @@ do icp=1,ncp_plane(ip) !loop on couples on plane icp icx=clx(ip,icp) icy=cly(ip,icp) + if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next if(LADDER(icx).ne.nldt.or. !If the ladder number does not match c $ cl_used(icx).eq.1.or. !or the X cluster is already used c $ cl_used(icy).eq.1.or. !or the Y cluster is already used - $ cl_used(icx).ne.0.or. !or the X cluster is already used !(3) - $ cl_used(icy).ne.0.or. !or the Y cluster is already used !(3) + $ cl_used(icx).ne.0.or. !or the X cluster is already used + $ cl_used(icy).ne.0.or. !or the Y cluster is already used $ .false.)goto 1188 !then jump to next couple. * call xyz_PAM(icx,icy,ist, $ PFA,PFA, -c $ 'ETA2','ETA2', $ AXV_STORE(nplanes-ip+1,ibest), - $ AYV_STORE(nplanes-ip+1,ibest)) + $ AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) distance = distance_to(XP,YP) +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI id=id_cp(ip,icp,ist) - if(DEBUG)print*,'( couple ',id - $ ,' ) normalized distance ',distance + if(DEBUG.EQ.1) + $ print*,'( couple ',id + $ ,' ) distance ',distance if(distance.lt.distmin)then xmm = xPAM ymm = yPAM @@ -3253,35 +3396,34 @@ rymm = resyPAM distmin = distance idm = id -c dedxmm = (dedx(icx)+dedx(icy))/2. !(1) - dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) - dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) + dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) + dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) + clincnewc=10.*dsqrt(rymm**2+rxmm**2 + $ +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))) ! EM GCC4.7 endif 1188 continue enddo !end loop on couples on plane icp - if(distmin.le.clinc)then + if(distmin.le.clincnewc)then * ----------------------------------- - xm(nplanes-ip+1) = xmm !<<< - ym(nplanes-ip+1) = ymm !<<< - zm(nplanes-ip+1) = zmm !<<< - xgood(nplanes-ip+1) = 1 !<<< - ygood(nplanes-ip+1) = 1 !<<< - resx(nplanes-ip+1)=rxmm !<<< - resy(nplanes-ip+1)=rymm !<<< -c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1) - dedxtrk_x(nplanes-ip+1) = dedxmmx !(1) - dedxtrk_y(nplanes-ip+1) = dedxmmy !(1) + xm(nplanes-ip+1) = xmm !<<< + ym(nplanes-ip+1) = ymm !<<< + zm(nplanes-ip+1) = zmm !<<< + xgood(nplanes-ip+1) = 1 !<<< + ygood(nplanes-ip+1) = 1 !<<< + resx(nplanes-ip+1)=rxmm !<<< + resy(nplanes-ip+1)=rymm !<<< + dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< + dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< * ----------------------------------- CP_STORE(nplanes-ip+1,ibest)=idm - if(DEBUG)print*,'%%%% included couple ',idm - $ ,' (norm.dist.= ',distmin,', cut ',clinc,' )' + if(DEBUG.EQ.1)print*,'%%%% included couple ',idm + $ ,' (dist.= ',distmin,', cut ',clincnewc,' )' goto 133 !next plane endif * ================================================ * STEP 2 >>>>>>> try to include a single cluster * either from a couple or single * ================================================ -c if(DEBUG)print*,'>>>> try to include a new cluster' distmin=1000000. xmm_A = 0. !--------------------------- ymm_A = 0. ! init variables that @@ -3300,6 +3442,7 @@ do icp=1,ncp_plane(ip) !loop on cluster inside couples icx=clx(ip,icp) icy=cly(ip,icp) + if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next id=id_cp(ip,icp,ist) if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match * !jump to the next couple @@ -3307,14 +3450,20 @@ c if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used !(3) * !jump to the Y cluster +c call xyz_PAM(icx,0,ist, +c $ PFA,PFA, +c $ AXV_STORE(nplanes-ip+1,ibest),0.) call xyz_PAM(icx,0,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ AXV_STORE(nplanes-ip+1,ibest),0.) + $ AXV_STORE(nplanes-ip+1,ibest),0., + $ bxyz(1), + $ bxyz(2) + $ ) distance = distance_to(XP,YP) -c if(DEBUG)print*,'normalized distance ',distance - if(DEBUG)print*,'( cl-X ',icx - $ ,' in cp ',id,' ) normalized distance ',distance +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI + if(DEBUG.EQ.1) + $ print*,'( cl-X ',icx + $ ,' in cp ',id,' ) distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A @@ -3326,8 +3475,8 @@ rymm = resyPAM distmin = distance iclm = icx -c dedxmm = dedx(icx) !(1) - dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) +c dedxmm = sgnl(icx) !(1) + dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = 0. !(1) endif 11881 continue @@ -3335,13 +3484,20 @@ c if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3) * !jump to the next couple +c call xyz_PAM(0,icy,ist, +c $ PFA,PFA, +c $ 0.,AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(0,icy,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ 0.,AYV_STORE(nplanes-ip+1,ibest)) + $ 0.,AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) distance = distance_to(XP,YP) - if(DEBUG)print*,'( cl-Y ',icy - $ ,' in cp ',id,' ) normalized distance ',distance +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI + if(DEBUG.EQ.1) + $ print*,'( cl-Y ',icy + $ ,' in cp ',id,' ) distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A @@ -3353,34 +3509,39 @@ rymm = resyPAM distmin = distance iclm = icy -c dedxmm = dedx(icy) !(1) +c dedxmm = sgnl(icy) !(1) dedxmmx = 0. !(1) - dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) + dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) endif 11882 continue enddo !end loop on cluster inside couples -*----- single clusters ----------------------------------------------- +*----- single clusters ----------------------------------------------- do ic=1,ncls(ip) !loop on single clusters icl=cls(ip,ic) -c if(cl_used(icl).eq.1.or. !if the cluster is already used if(cl_used(icl).ne.0.or. !if the cluster is already used !(3) $ LADDER(icl).ne.nldt.or. !or the ladder number does not match $ .false.)goto 18882 !jump to the next singlet if(mod(VIEW(icl),2).eq.0)then!<---- X view call xyz_PAM(icl,0,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ AXV_STORE(nplanes-ip+1,ibest),0.) + $ AXV_STORE(nplanes-ip+1,ibest),0., + $ bxyz(1), + $ bxyz(2) + $ ) else !<---- Y view call xyz_PAM(0,icl,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ 0.,AYV_STORE(nplanes-ip+1,ibest)) + $ 0.,AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) endif distance = distance_to(XP,YP) - if(DEBUG)print*,'( cl-s ',icl - $ ,' ) normalized distance ',distance +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI + if(DEBUG.EQ.1) + $ print*,'( cl-s ',icl + $ ,' ) distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A @@ -3392,43 +3553,61 @@ rymm = resyPAM distmin = distance iclm = icl -c dedxmm = dedx(icl) !(1) if(mod(VIEW(icl),2).eq.0)then !<---- X view - dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) - dedxmmy = 0. !(1) + dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) + dedxmmy = 0. else !<---- Y view - dedxmmx = 0. !(1) - dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) + dedxmmx = 0. + dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) endif endif 18882 continue enddo !end loop on single clusters - if(distmin.le.clinc)then - - CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< -* ---------------------------- + if(iclm.ne.0)then 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,' )' - else - YGOOD(nplanes-ip+1)=1. - resy(nplanes-ip+1)=rymm - if(DEBUG)print*,'%%%% included Y-cl ',iclm - $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )' + clincnew= + $ 20.* !EM GCC4.7 + $ dsqrt(rxmm**2 + + $ DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1)) + else if(mod(VIEW(iclm),2).ne.0)then + clincnew= + $ 10.* !EM GCC4.7 + $ dsqrt(rymm**2 + + $ DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2)) + endif + + if(distmin.le.clincnew)then + + CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< +* ---------------------------- + if(mod(VIEW(iclm),2).eq.0)then + XGOOD(nplanes-ip+1)=1. + resx(nplanes-ip+1)=rxmm + if(DEBUG.EQ.1) + $ print*,'%%%% included X-cl ',iclm + $ ,'( chi^2, ',RCHI2_STORE(ibest) + $ ,', dist.= ',distmin + $ ,', cut ',clincnew,' )' + else + YGOOD(nplanes-ip+1)=1. + resy(nplanes-ip+1)=rymm + if(DEBUG.EQ.1) + $ print*,'%%%% included Y-cl ',iclm + $ ,'( chi^2, ',RCHI2_STORE(ibest) + $ ,', dist.= ', distmin + $ ,', cut ',clincnew,' )' + endif +* ---------------------------- + xm_A(nplanes-ip+1) = xmm_A + ym_A(nplanes-ip+1) = ymm_A + xm_B(nplanes-ip+1) = xmm_B + ym_B(nplanes-ip+1) = ymm_B + zm(nplanes-ip+1) = (zmm_A+zmm_B)/2. + dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< + dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< +* ---------------------------- endif -* ---------------------------- - xm_A(nplanes-ip+1) = xmm_A - ym_A(nplanes-ip+1) = ymm_A - xm_B(nplanes-ip+1) = xmm_B - ym_B(nplanes-ip+1) = ymm_B - zm(nplanes-ip+1) = (zmm_A+zmm_B)/2. -c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1) - dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1) - dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1) -* ---------------------------- endif endif 133 continue @@ -3439,6 +3618,7 @@ return end + *************************************************** * * * * @@ -3447,187 +3627,7 @@ * * * * ************************************************** -cccccc 12/08/2006 modified by elena ---> (1) * - subroutine clean_XYclouds(ibest,iflag) - - include 'commontracker.f' - include 'common_momanhough.f' - 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 - - id=CP_STORE(nplanes-ip+1,ibest) - icl=CLS_STORE(nplanes-ip+1,ibest) - if(id.ne.0.or.icl.ne.0)then - if(id.ne.0)then - iclx=clx(ip,icp_cp(id)) - icly=cly(ip,icp_cp(id)) -c cl_used(iclx)=1 !tag used clusters -c cl_used(icly)=1 !tag used clusters - cl_used(iclx)=ntrk !tag used clusters !(1) - cl_used(icly)=ntrk !tag used clusters !(1) - elseif(icl.ne.0)then -c cl_used(icl)=1 !tag used clusters - cl_used(icl)=ntrk !tag used clusters !1) - endif - -c if(DEBUG)then -c print*,ip,' <<< ',id -c endif -* ----------------------------- -* remove the couple from clouds -* remove also vitual couples containing the -* selected clusters -* ----------------------------- - do icp=1,ncp_plane(ip) - if( - $ clx(ip,icp).eq.iclx - $ .or. - $ clx(ip,icp).eq.icl - $ .or. - $ cly(ip,icp).eq.icly - $ .or. - $ cly(ip,icp).eq.icl - $ )then - id=id_cp(ip,icp,1) - if(DEBUG)then - print*,ip,' <<< cp ',id - $ ,' ( cl-x ' - $ ,clx(ip,icp) - $ ,' cl-y ' - $ ,cly(ip,icp),' ) --> removed' - endif -* ----------------------------- -* remove the couple from clouds - do iyz=1,nclouds_yz - if(cpcloud_yz(iyz,abs(id)).ne.0)then - ptcloud_yz(iyz)=ptcloud_yz(iyz)-1 - cpcloud_yz(iyz,abs(id))=0 - endif - enddo - do ixz=1,nclouds_xz - if(cpcloud_xz(ixz,abs(id)).ne.0)then - ptcloud_xz(ixz)=ptcloud_xz(ixz)-1 - cpcloud_xz(ixz,abs(id))=0 - endif - enddo -* ----------------------------- - endif - enddo - - endif - enddo !end loop on planes - - return - end - - - - -c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** -c$$$ real function fbad_cog(ncog,ic) -c$$$ -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'level1.f' -c$$$ include 'calib.f' -c$$$ -c$$$* --> signal of the central strip -c$$$ sc = CLSIGNAL(INDMAX(ic)) !center -c$$$ -c$$$* signal of adjacent strips -c$$$* --> left -c$$$ sl1 = 0 !left 1 -c$$$ if( -c$$$ $ (INDMAX(ic)-1).ge.INDSTART(ic) -c$$$ $ ) -c$$$ $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1)) -c$$$ -c$$$ sl2 = 0 !left 2 -c$$$ if( -c$$$ $ (INDMAX(ic)-2).ge.INDSTART(ic) -c$$$ $ ) -c$$$ $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2)) -c$$$ -c$$$* --> right -c$$$ sr1 = 0 !right 1 -c$$$ if( -c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) -c$$$ $ .or. -c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) -c$$$ $ ) -c$$$ $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1)) -c$$$ -c$$$ sr2 = 0 !right 2 -c$$$ if( -c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) -c$$$ $ .or. -c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) -c$$$ $ ) -c$$$ $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2)) -c$$$ -c$$$ -c$$$ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view -c$$$ f = 4. -c$$$ si = 8.4 -c$$$ else !X-view -c$$$ f = 6. -c$$$ si = 3.9 -c$$$ endif -c$$$ -c$$$ fbad_cog = 1. -c$$$ f0 = 1 -c$$$ f1 = 1 -c$$$ f2 = 1 -c$$$ f3 = 1 -c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then -c$$$ -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f -c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f -c$$$ -c$$$ if(ncog.eq.2.and.sl1.ne.0)then -c$$$ fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.) -c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then -c$$$ fbad_cog = 1. -c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then -c$$$ fbad_cog = 1. -c$$$ else -c$$$ fbad_cog = 1. -c$$$ endif -c$$$ -c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then -c$$$ -c$$$ -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f -c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f -c$$$ -c$$$ if(ncog.eq.2.and.sr1.ne.0)then -c$$$ fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.) -c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then -c$$$ fbad_cog = 1. -c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then -c$$$ fbad_cog = 1. -c$$$ else -c$$$ fbad_cog = 1. -c$$$ endif -c$$$ -c$$$ endif -c$$$ -c$$$ fbad_cog = sqrt(fbad_cog) -c$$$ -c$$$ return -c$$$ end -c$$$ @@ -3635,63 +3635,53 @@ 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' +* --------------------------------- +* variables initialized from level1 +* --------------------------------- do i=1,nviews good2(i)=good1(i) + do j=1,nva1_view + vkflag(i,j)=1 + if(cnnev(i,j).le.0)then + vkflag(i,j)=cnnev(i,j) + endif + enddo 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***************************************************** - +* ---------------- +* level2 variables +* ---------------- 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 ZM_nt(IP,IT) = 0 RESX_nt(IP,IT) = 0 RESY_nt(IP,IT) = 0 + TAILX_nt(IP,IT) = 0 + TAILY_nt(IP,IT) = 0 + XBAD(IP,IT) = 0 + YBAD(IP,IT) = 0 XGOOD_nt(IP,IT) = 0 YGOOD_nt(IP,IT) = 0 -c***************************************************** -cccccc 11/9/2005 modified by david fedele + LS(IP,IT) = 0 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 + multmaxx(ip,it) = 0 + seedx(ip,it) = 0 + xpu(ip,it) = 0 + multmaxy(ip,it) = 0 + seedy(ip,it) = 0 + ypu(ip,it) = 0 enddo do ipa=1,5 AL_nt(IPA,IT) = 0 @@ -3700,25 +3690,22 @@ 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 + sxbad(ip)=0 + sybad(ip)=0 + multmaxsx(ip)=0 + multmaxsy(ip)=0 enddo -c******************************************************* end @@ -3733,6 +3720,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,ntrpt_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,ntrpt_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) * ------------------------------------------------------- @@ -3744,37 +3815,43 @@ include 'commontracker.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) + include 'calib.f' + + character*10 PFA + common/FINALPFA/PFA + + real sinth,phi,pig, npig ! EM GCC4.7 + integer ssensor,sladder pig=acos(-1.) -c good2=1!.true. - chi2_nt(ntr) = sngl(chi2) - nstep_nt(ntr) = 0!nstep - 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) -***************************************************** + +* ------------------------------------- + chi2_nt(ntr) = sngl(chi2) + nstep_nt(ntr) = nstep +* ------------------------------------- + phi = REAL(al(4)) + sinth = REAL(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 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 xgood_nt(ip,ntr) = int(xgood(ip)) ygood_nt(ip,ntr) = int(ygood(ip)) @@ -3783,44 +3860,156 @@ zm_nt(ip,ntr) = sngl(zm(ip)) RESX_nt(IP,ntr) = sngl(resx(ip)) RESY_nt(IP,ntr) = sngl(resy(ip)) + TAILX_nt(IP,ntr) = 0. + TAILY_nt(IP,ntr) = 0. xv_nt(ip,ntr) = sngl(xv(ip)) yv_nt(ip,ntr) = sngl(yv(ip)) 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) + ayv_nt(ip,ntr) = sngl(ayv(ip)) + + factor = sqrt( + $ tan( acos(-1.) * sngl(axv(ip)) /180. )**2 + + $ tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 + + $ 1. ) + + dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)/factor) + dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)/factor) - id = CP_STORE(ip,IDCAND) + +ccc print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr) + + ax = axv_nt(ip,ntr) + ay = ayv_nt(ip,ntr) + bfx = BX_STORE(ip,IDCAND) + bfy = BY_STORE(ip,IDCAND) +c$$$ if(ip.eq.6) ax = -1. * axv_nt(ip,ntr) +c$$$ if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND) +c$$$ tgtemp = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001 +c$$$ angx = 180.*atan(tgtemp)/acos(-1.) +c$$$ tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001 +c$$$ angy = 180.*atan(tgtemp)/acos(-1.) + + angx = effectiveangle(ax,2*ip,bfy) + angy = effectiveangle(ay,2*ip-1,bfx) + + + + id = CP_STORE(ip,IDCAND) ! couple id icl = CLS_STORE(ip,IDCAND) - if(id.ne.0)then + ssensor = -1 + sladder = -1 + ssensor = SENSOR_STORE(ip,IDCAND) + sladder = LADDER_STORE(ip,IDCAND) + if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align + LS(IP,ntr) = ssensor+10*sladder + +c if(id.ne.0)then +CCCCCC(10 novembre 2009) PATCH X NUCLEI +C non ho capito perche', ma durante il ritracciamento dei nuclei +C (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile +C che non viene reinizializzata correttamente e i cluster esclusi +C dal fit risultano ancora inclusi... +C + cltrx(ip,ntr) = 0 + cltry(ip,ntr) = 0 +c$$$ if( +c$$$ $ xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1 +c$$$ $ .and. +c$$$ $ id.ne.0)then + if(id.ne.0)then !patch 30/12/09 + +c >>> is a couple cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id)) cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id)) -c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr) - elseif(icl.ne.0)then - if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl - if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl -c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr) + + if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then + + cl_used(cltrx(ip,ntr)) = 1 !tag used clusters + + xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id))) + + if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0) + $ dedx_x(ip,ntr)=-dedx_x(ip,ntr) + + multmaxx(ip,ntr) = maxs(cltrx(ip,ntr)) + $ +10000*mult(cltrx(ip,ntr)) + seedx(ip,ntr) = clsignal(indmax(cltrx(ip,ntr))) + $ /clsigma(indmax(cltrx(ip,ntr))) + call applypfa(PFA,cltrx(ip,ntr),angx,corr,res) + xpu(ip,ntr) = corr + + endif + if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then + + cl_used(cltry(ip,ntr)) = 1 !tag used clusters + + ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id))) + + if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0) + $ dedx_y(ip,ntr)=-dedx_y(ip,ntr) + + multmaxy(ip,ntr) = maxs(cltry(ip,ntr)) + $ +10000*mult(cltry(ip,ntr)) + seedy(ip,ntr) = clsignal(indmax(cltry(ip,ntr))) + $ /clsigma(indmax(cltry(ip,ntr))) + call applypfa(PFA,cltry(ip,ntr),angy,corr,res) + ypu(ip,ntr) = corr + endif + +c$$$ elseif(icl.ne.0)then + endif !patch 30/12/09 + if(icl.ne.0)then !patch 30/12/09 + + cl_used(icl) = 1 !tag used clusters + + if(mod(VIEW(icl),2).eq.0)then + cltrx(ip,ntr)=icl + xbad(ip,ntr) = nbadstrips(4,icl) + + if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr) + + multmaxx(ip,ntr) = maxs(cltrx(ip,ntr)) + $ +10000*mult(cltrx(ip,ntr)) + seedx(ip,ntr) = clsignal(indmax(cltrx(ip,ntr))) + $ /clsigma(indmax(cltrx(ip,ntr))) + call applypfa(PFA,cltrx(ip,ntr),angx,corr,res) + xpu(ip,ntr) = corr + + elseif(mod(VIEW(icl),2).eq.1)then + cltry(ip,ntr)=icl + ybad(ip,ntr) = nbadstrips(4,icl) + + if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr) + + multmaxy(ip,ntr) = maxs(cltry(ip,ntr)) + $ +10000*mult(cltry(ip,ntr)) + seedy(ip,ntr) = clsignal(indmax(cltry(ip,ntr))) + $ /clsigma(indmax(cltry(ip,ntr))) + call applypfa(PFA,cltry(ip,ntr),angy,corr,res) + ypu(ip,ntr) = corr + + endif + 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) + if(DEBUG.eq.1)then + print*,'> STORING TRACK ',ntr + print*,'clusters: ' + do ip=1,6 + print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr) + enddo + print*,'dedx: ' + do ip=1,6 + print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr) + enddo + endif 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 @@ -3829,62 +4018,155 @@ * ------------------------------------------------------- include 'commontracker.f' - include 'level1.f' - include 'level2.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 -c good2=1!.true. nclsx = 0 nclsy = 0 + do iv = 1,nviews +c if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv) + good2(iv) = good2(iv) + mask_view(iv)*2**8 + enddo + + if(DEBUG.eq.1)then + print*,'> STORING SINGLETS ' + endif + do icl=1,nclstr1 + + ip=nplanes-npl(VIEW(icl))+1 + if(cl_used(icl).eq.0)then !cluster not included in any track - ip=nplanes-npl(VIEW(icl))+1 + if(mod(VIEW(icl),2).eq.0)then !=== X views + nclsx = nclsx + 1 planex(nclsx) = ip - sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2) + sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) + if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx) clsx(nclsx) = icl + sxbad(nclsx) = nbadstrips(1,icl) + multmaxsx(nclsx) = maxs(icl)+10000*mult(icl) + + do is=1,2 c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.) - call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.) - xs(is,nclsx) = (xPAM_A+xPAM_B)/2 +c call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.) + call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.) + xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7 enddo -c$$$ print*,'nclsx ',nclsx -c$$$ print*,'planex(nclsx) ',planex(nclsx) -c$$$ print*,'sgnlxs(nclsx) ',sgnlxs(nclsx) -c$$$ print*,'xs(1,nclsx) ',xs(1,nclsx) -c$$$ print*,'xs(2,nclsx) ',xs(2,nclsx) else !=== Y views nclsy = nclsy + 1 planey(nclsy) = ip - sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2) + sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) + if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy) clsy(nclsy) = icl + sybad(nclsy) = nbadstrips(1,icl) + multmaxsy(nclsy) = maxs(icl)+10000*mult(icl) + + do is=1,2 c call xyz_PAM(0,icl,is,' ','COG1',0.,0.) - call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.) - ys(is,nclsy) = (yPAM_A+yPAM_B)/2 +c call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.) + call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.) + ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7 enddo -c$$$ print*,'nclsy ',nclsy -c$$$ print*,'planey(nclsy) ',planey(nclsy) -c$$$ print*,'sgnlys(nclsy) ',sgnlys(nclsy) -c$$$ print*,'ys(1,nclsy) ',ys(1,nclsy) -c$$$ print*,'ys(2,nclsy) ',ys(2,nclsy) endif endif -c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip) ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO whichtrack(icl) = cl_used(icl) - +* -------------------------------------------------- +* per non perdere la testa... +* whichtrack(icl) e` una variabile del common level1 +* che serve solo per sapere quali cluster sono stati +* associati ad una traccia, e permettere di salvare +* solo questi nell'albero di uscita +* -------------------------------------------------- + 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,min(ndblt_max_nt,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,min(ntrpt_max_nt,nnn) + tr_cloud_nt(ipt)=tr_cloud(ipt) + enddo + endif + endif + end +