--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/08/28 13:25:50 1.30 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2009/12/30 14:22:23 1.44 @@ -20,8 +20,6 @@ include 'calib.f' include 'level2.f' - - c print*,'======================================================' c$$$ do ic=1,NCLSTR1 c$$$ if(.false. @@ -74,12 +72,11 @@ *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- - call cl_to_couples(iflag) if(iflag.eq.1)then !bad event goto 880 !go to next event endif - + if(ncp_tot.eq.0)goto 880 !go to next event *----------------------------------------------------- *----------------------------------------------------- * HOUGH TRASFORM @@ -113,6 +110,7 @@ if(iflag.eq.1)then !bad event goto 880 !go to next event endif + if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event *------------------------------------------------------------------------------- @@ -163,37 +161,65 @@ 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 - goto 878 - endif + 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 + + if(nclouds_yz.eq.0)goto 880 !go to next event - if(planehit.eq.3) goto 881 + +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 - 881 continue -* if there is at least three planes on the Y view decreases cuts on X view - if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and. - $ nplxz_min.ne.y_min_start)then - nptxz_min=x_min_step - nplxz_min=x_min_start-x_min_step - goto 879 - endif +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 @@ -243,15 +269,20 @@ * *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- - 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 @@ -272,7 +303,7 @@ * 1st) decreasing n.points * 2nd) increasing chi**2 * ------------------------------------------------------- - rchi2best=1000000000. + rchi2best=1000000000. ndofbest=0 do i=1,ntracks ndof=0 @@ -295,24 +326,6 @@ endif enddo -c$$$ rchi2best=1000000000. -c$$$ ndofbest=0 !(1) -c$$$ do i=1,ntracks -c$$$ if(RCHI2_STORE(i).lt.rchi2best.and. -c$$$ $ RCHI2_STORE(i).gt.0)then -c$$$ ndof=0 !(1) -c$$$ do ii=1,nplanes !(1) -c$$$ ndof=ndof !(1) -c$$$ $ +int(xgood_store(ii,i)) !(1) -c$$$ $ +int(ygood_store(ii,i)) !(1) -c$$$ enddo !(1) -c$$$ if(ndof.ge.ndofbest)then !(1) -c$$$ ibest=i -c$$$ rchi2best=RCHI2_STORE(i) -c$$$ ndofbest=ndof !(1) -c$$$ endif !(1) -c$$$ endif -c$$$ enddo if(ibest.eq.0)goto 880 !>> no good candidates *------------------------------------------------------------------------------- @@ -350,7 +363,6 @@ do i=1,5 AL_GUESS(i)=AL(i) enddo -c print*,'## guess: ',al do i=1,5 AL(i)=dble(AL_STORE(i,icand)) @@ -365,20 +377,14 @@ if(VERBOSE.EQ.1)iprint=1 if(DEBUG.EQ.1)iprint=2 call mini2(jstep,ifail,iprint) - if(ifail.ne.0) then +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 *** (after refinement) ' $ ,iev - -c$$$ print*,'guess: ',(al_guess(i),i=1,5) -c$$$ print*,'previous: ',(al_store(i,icand),i=1,5) -c$$$ print*,'result: ',(al(i),i=1,5) -c$$$ print*,'xgood ',xgood -c$$$ print*,'ygood ',ygood -c$$$ print*,'----------------------------------------------' endif -c chi2=-chi2 endif if(DEBUG.EQ.1)then @@ -424,37 +430,14 @@ if(WARNING.EQ.1)print*,'** WARNING ** sensor=0' goto 122 !jump endif -c print*,'is1 ',is1 do ip=1,NPLANES is2 = SENSOR_STORE(ip,icand) !sensor -c print*,'is2 ',is2,' ip ',ip 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! ' -* now search for track-image among track candidates -c$$$ do i=1,ntracks -c$$$ iimage=i -c$$$ do ip=1,nplanes -c$$$ if( CP_STORE(nplanes-ip+1,icand).ne. -c$$$ $ -1*CP_STORE(nplanes-ip+1,i).and. -c$$$ $ CP_STORE(nplanes-ip+1,i).ne.0.and. -c$$$ $ CP_STORE(nplanes-ip+1,icand).ne.0 )iimage=0 -c$$$ print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i) -c$$$ $ ,CP_STORE(nplanes-ip+1,icand),iimage -c$$$ enddo -c$$$ if( iimage.ne.0.and. -c$$$c $ RCHI2_STORE(i).le.CHI2MAX.and. -c$$$c $ RCHI2_STORE(i).gt.0.and. -c$$$ $ .true.)then -c$$$ if(DEBUG.EQ.1)print*,'Track candidate ',iimage -c$$$ $ ,' >>> TRACK IMAGE >>> of' -c$$$ $ ,ibest -c$$$ goto 122 !image track found -c$$$ endif -c$$$ enddo * --------------------------------------------------------------- * take the candidate with the greatest number of matching couples * if more than one satisfies the requirement, @@ -486,8 +469,6 @@ endif endif -c$$$ print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i) -c$$$ $ ,CP_STORE(nplanes-ip+1,icand),ncp enddo 117 continue trackimage(i)=ncp !number of matching couples @@ -501,8 +482,9 @@ 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.1)then + if(ngood.gt.0)then * ------------------------------------------------------- * order track-candidates according to: * 1st) decreasing n.points @@ -532,27 +514,20 @@ endif endif enddo + + if(DEBUG.EQ.1)then + print*,'Track candidate ',iimage + $ ,' >>> TRACK IMAGE >>> of' + $ ,ibest + endif endif - if(DEBUG.EQ.1)print*,'Track candidate ',iimage - $ ,' >>> TRACK IMAGE >>> of' - $ ,ibest 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) - if(ntrk.eq.NTRKMAX)then if(verbose.eq.1) $ print*, @@ -562,43 +537,29 @@ 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.EQ.1)then - print*,'***** NEW SEARCH ****' - endif - goto 11111 !try new search - - endif -* ********************************************** - - 880 return end @@ -688,6 +649,9 @@ 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 @@ -696,21 +660,19 @@ parameter (ndivx=30) -c$$$ print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy 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. -c print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy + 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 ' @@ -727,7 +689,7 @@ viewx = VIEW(icx) nldx = nld(MAXS(icx),VIEW(icx)) nplx = npl(VIEW(icx)) - resxPAM = RESXAV +c resxPAM = RESXAV stripx = float(MAXS(icx)) if( @@ -747,125 +709,15 @@ * -------------------------- * magnetic-field corrections * -------------------------- - angtemp = ax - bfytemp = bfy -* ///////////////////////////////// -* AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!! -* *grvzkkjsdgjhhhgngbn###>:( -* ///////////////////////////////// -c if(nplx.eq.6) angtemp = -1. * ax -c if(nplx.eq.6) bfytemp = -1. * bfy - if(viewx.eq.12) angtemp = -1. * ax - if(viewx.eq.12) bfytemp = -1. * bfy - tgtemp = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001 - angx = 180.*atan(tgtemp)/acos(-1.) - stripx = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX -c$$$ print*,nplx,ax,bfy/10. -c$$$ print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX -c$$$ print*,'========================' -c$$$ if(bfy.ne.0.)print*,viewx,'-x- ' -c$$$ $ ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ -* -------------------------- - -c$$$ print*,'--- x-cl ---' -c$$$ istart = INDSTART(ICX) -c$$$ istop = TOTCLLENGTH -c$$$ if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1 -c$$$ print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) -c$$$ print*,(CLSIGNAL(i),i=istart,istop) -c$$$ print*,INDMAX(icx) -c$$$ print*,cog(4,icx) -c$$$ print*,fbad_cog(4,icx) + stripx = stripx + fieldcorr(viewx,bfy) + angx = effectiveangle(ax,viewx,bfy) + call applypfa(PFAx,icx,angx,corr,res) + stripx = stripx + corr + resxPAM = res - if(PFAx.eq.'COG1')then - - stripx = stripx - resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM - - elseif(PFAx.eq.'COG2')then - - stripx = stripx + cog(2,icx) - resxPAM = risx_cog(abs(angx))!TEMPORANEO - resxPAM = resxPAM*fbad_cog(2,icx) - - elseif(PFAx.eq.'COG3')then - - stripx = stripx + cog(3,icx) - resxPAM = risx_cog(abs(angx))!TEMPORANEO - resxPAM = resxPAM*fbad_cog(3,icx) - - elseif(PFAx.eq.'COG4')then - - stripx = stripx + cog(4,icx) - resxPAM = risx_cog(abs(angx))!TEMPORANEO - resxPAM = resxPAM*fbad_cog(4,icx) - - elseif(PFAx.eq.'ETA2')then - - stripx = stripx + pfaeta2(icx,angx) - resxPAM = risxeta2(abs(angx)) - resxPAM = resxPAM*fbad_cog(2,icx) -c$$$ if(DEBUG.EQ.1.and.fbad_cog(2,icx).ne.1) -c$$$ $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) - - elseif(PFAx.eq.'ETA3')then - - stripx = stripx + pfaeta3(icx,angx) - resxPAM = risxeta3(abs(angx)) - resxPAM = resxPAM*fbad_cog(3,icx) -c if(DEBUG.and.fbad_cog(3,icx).ne.1) -c $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx) - - elseif(PFAx.eq.'ETA4')then - - stripx = stripx + pfaeta4(icx,angx) - resxPAM = risxeta4(abs(angx)) - resxPAM = resxPAM*fbad_cog(4,icx) -c if(DEBUG.and.fbad_cog(4,icx).ne.1) -c $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx) - - elseif(PFAx.eq.'ETA')then - - stripx = stripx + pfaeta(icx,angx) -c resxPAM = riseta(icx,angx) - resxPAM = riseta(viewx,angx) - resxPAM = resxPAM*fbad_eta(icx,angx) -c if(DEBUG.and.fbad_cog(2,icx).ne.1) -c $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) - - elseif(PFAx.eq.'ETAL')then - - stripx = stripx + pfaetal(icx,angx) - resxPAM = riseta(viewx,angx) - resxPAM = resxPAM*fbad_eta(icx,angx) -c if(DEBUG.and.fbad_cog(2,icx).ne.1) -c $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) - - elseif(PFAx.eq.'COG')then - - stripx = stripx + cog(0,icx) - resxPAM = risx_cog(abs(angx)) - resxPAM = resxPAM*fbad_cog(0,icx) - - else - if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx - endif - - -* ====================================== -* temporary patch for saturated clusters -* ====================================== - if( nsatstrips(icx).gt.0 )then - stripx = stripx + cog(4,icx) - resxPAM = pitchX*1e-4/sqrt(12.) -cc cc=cog(4,icx) -c$$$ print*,icx,' *** ',cc -c$$$ print*,icx,' *** ',resxPAM - endif - - 10 endif - + 10 continue + endif * ----------------- * CLUSTER Y @@ -876,7 +728,7 @@ viewy = VIEW(icy) nldy = nld(MAXS(icy),VIEW(icy)) nply = npl(VIEW(icy)) - resyPAM = RESYAV +c resyPAM = RESYAV stripy = float(MAXS(icy)) if( @@ -900,114 +752,20 @@ endif goto 100 endif + * -------------------------- * magnetic-field corrections * -------------------------- - tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001 - angy = 180.*atan(tgtemp)/acos(-1.) - stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY -c$$$ if(bfx.ne.0.)print*,viewy,'-y- ' -c$$$ $ ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ -* -------------------------- + stripy = stripy + fieldcorr(viewy,bfx) + angy = effectiveangle(ay,viewy,bfx) -c$$$ print*,'--- y-cl ---' -c$$$ istart = INDSTART(ICY) -c$$$ istop = TOTCLLENGTH -c$$$ if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1 -c$$$ print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) -c$$$ print*,(CLSIGNAL(i),i=istart,istop) -c$$$ print*,INDMAX(icy) -c$$$ print*,cog(4,icy) -c$$$ print*,fbad_cog(4,icy) - - if(PFAy.eq.'COG1')then - - stripy = stripy - resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM - - elseif(PFAy.eq.'COG2')then - - stripy = stripy + cog(2,icy) - resyPAM = risy_cog(abs(angy))!TEMPORANEO - resyPAM = resyPAM*fbad_cog(2,icy) - - elseif(PFAy.eq.'COG3')then - - stripy = stripy + cog(3,icy) - resyPAM = risy_cog(abs(angy))!TEMPORANEO - resyPAM = resyPAM*fbad_cog(3,icy) - - elseif(PFAy.eq.'COG4')then - - stripy = stripy + cog(4,icy) - resyPAM = risy_cog(abs(angy))!TEMPORANEO - resyPAM = resyPAM*fbad_cog(4,icy) - - elseif(PFAy.eq.'ETA2')then - - stripy = stripy + pfaeta2(icy,angy) - resyPAM = risyeta2(abs(angy)) - resyPAM = resyPAM*fbad_cog(2,icy) -c if(DEBUG.and.fbad_cog(2,icy).ne.1) -c $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy) - - elseif(PFAy.eq.'ETA3')then - - stripy = stripy + pfaeta3(icy,angy) - resyPAM = resyPAM*fbad_cog(3,icy) -c if(DEBUG.and.fbad_cog(3,icy).ne.1) -c $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy) - - elseif(PFAy.eq.'ETA4')then - - stripy = stripy + pfaeta4(icy,angy) - resyPAM = resyPAM*fbad_cog(4,icy) -c if(DEBUG.and.fbad_cog(4,icy).ne.1) -c $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy) - - elseif(PFAy.eq.'ETA')then - - stripy = stripy + pfaeta(icy,angy) -c resyPAM = riseta(icy,angy) - resyPAM = riseta(viewy,angy) - resyPAM = resyPAM*fbad_eta(icy,angy) -c if(DEBUG.and.fbad_cog(2,icy).ne.1) -c $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy) - - elseif(PFAy.eq.'ETAL')then + call applypfa(PFAy,icy,angy,corr,res) + stripy = stripy + corr + resyPAM = res - stripy = stripy + pfaetal(icy,angy) - resyPAM = riseta(viewy,angy) - resyPAM = resyPAM*fbad_eta(icy,angy) -c if(DEBUG.and.fbad_cog(2,icy).ne.1) -c $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy) - - elseif(PFAy.eq.'COG')then - - stripy = stripy + cog(0,icy) - resyPAM = risy_cog(abs(angy)) - resyPAM = resyPAM*fbad_cog(0,icy) - - else - if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx - endif - - -* ====================================== -* temporary patch for saturated clusters -* ====================================== - if( nsatstrips(icy).gt.0 )then - stripy = stripy + cog(4,icy) - resyPAM = pitchY*1e-4/sqrt(12.) -cc cc=cog(4,icy) -c$$$ print*,icy,' *** ',cc -c$$$ print*,icy,' *** ',resyPAM - endif - - - 20 endif + 20 continue + endif -c$$$ print*,'## stripx,stripy ',stripx,stripy c=========================================================== C COUPLE @@ -1024,11 +782,10 @@ $ ' 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------------------------------------------------------------------------ @@ -1062,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. @@ -1088,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 @@ -1098,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=========================================================== @@ -1110,8 +867,6 @@ 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... if(DEBUG.EQ.1) then @@ -1119,7 +874,10 @@ $ ' 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 @@ -1135,8 +893,6 @@ yi_B = yi endif -c print*,'X-cl ',icx,stripx,' --> ',xi -c print*,yi_A,' <--> ',yi_B else if(DEBUG.EQ.1) then @@ -1182,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 @@ -1192,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 @@ -1205,7 +982,6 @@ 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 if(DEBUG.EQ.1) then @@ -1216,9 +992,6 @@ endif -c print*,'## xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM -c print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A -c print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B 100 continue end @@ -1256,7 +1029,7 @@ if(icx.gt.nclstr1.or.icy.gt.nclstr1)then print*,'xyzpam: ***WARNING*** clusters ',icx,icy - $ ,' does not exists (nclstr1=',nclstr1,')' + $ ,' do not exists (n.clusters=',nclstr1,')' icx = -1*icx icy = -1*icy return @@ -1266,19 +1039,11 @@ call idtoc(pfaid,PFAx) call idtoc(pfaid,PFAy) -cc print*,PFAx,PFAy - -c$$$ call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) - -c$$$ print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy if(icx.ne.0.and.icy.ne.0)then ipx=npl(VIEW(icx)) ipy=npl(VIEW(icy)) -c$$$ if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip ) -c$$$ $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy -c$$$ $ ,' does not belong to the correct plane: ',ip,ipx,ipy if( (nplanes-ipx+1).ne.ip )then print*,'xyzpam: ***WARNING*** cluster ',icx @@ -1303,19 +1068,16 @@ xm(ip) = xPAM ym(ip) = yPAM zm(ip) = zPAM - xm_A(ip) = 0. - ym_A(ip) = 0. - xm_B(ip) = 0. - ym_B(ip) = 0. + 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)) -c$$$ if((nplanes-ipy+1).ne.ip) -c$$$ $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy -c$$$ $ ,' does not belong to the correct plane: ',ip,ipx,ipy if( (nplanes-ipy+1).ne.ip )then print*,'xyzpam: ***WARNING*** cluster ',icy $ ,' does not belong to plane: ',ip @@ -1330,9 +1092,12 @@ resx(ip) = 1000. resy(ip) = resyPAM - xm(ip) = -100. - ym(ip) = -100. - zm(ip) = (zPAM_A+zPAM_B)/2. +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 @@ -1343,9 +1108,6 @@ elseif(icx.ne.0.and.icy.eq.0)then ipx=npl(VIEW(icx)) -c$$$ if((nplanes-ipx+1).ne.ip) -c$$$ $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy -c$$$ $ ,' does not belong to the correct plane: ',ip,ipx,ipy if( (nplanes-ipx+1).ne.ip )then print*,'xyzpam: ***WARNING*** cluster ',icx @@ -1361,9 +1123,12 @@ resx(ip) = resxPAM resy(ip) = 1000. - xm(ip) = -100. - ym(ip) = -100. - zm(ip) = (zPAM_A+zPAM_B)/2. +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 @@ -1377,7 +1142,6 @@ if(lad.ne.0)il=lad is = 1 if(sensor.ne.0)is=sensor -c print*,nplanes-ip+1,il,is xgood(ip) = 0. ygood(ip) = 0. @@ -1397,12 +1161,10 @@ endif if(DEBUG.EQ.1)then -c print*,'----------------------------- track coord' 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) -c$$$ print*,'-----------------------------' endif end @@ -1426,7 +1188,7 @@ * ******************************************************************************** - real function distance_to(XPP,YPP) + real function distance_to(rXPP,rYPP) include 'common_xyzPAM.f' @@ -1435,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. @@ -1484,16 +1251,13 @@ 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. @@ -1514,17 +1278,9 @@ c$$$ $ ((yPAM-YPP)/resyPAM)**2 distance=dsqrt(distance) -c$$$ print*,xPAM,yPAM,zPAM -c$$$ $ ,' --- distance_to --- ',xpp,ypp -c$$$ print*,' resolution ',resxPAM,resyPAM else -c print* -c $ ,' function distance_to ---> wrong usage!!!' -c print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM -c print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b ' -c $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b endif distance_to = sngl(distance) @@ -1564,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) @@ -1589,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... -c print*,'whichsensor: ', -c $ ' 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------------------------------------------------------------------------ @@ -1622,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. @@ -1748,7 +1496,6 @@ is_cp=0 if(id.lt.0)is_cp=1 if(id.gt.0)is_cp=2 -c if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!' return end @@ -1847,11 +1594,13 @@ integer iflag 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 @@ -1864,8 +1613,8 @@ 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 @@ -1893,6 +1642,8 @@ 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) * ---------------------------------------------------- @@ -1900,7 +1651,7 @@ * ---------------------------------------------------- * cut on charge (X VIEW) * ---------------------------------------------------- - if(sgnl(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 @@ -1941,6 +1692,8 @@ 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) @@ -1950,7 +1703,7 @@ * ---------------------------------------------------- * cut on charge (Y VIEW) * ---------------------------------------------------- - if(sgnl(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 @@ -2028,8 +1781,8 @@ $ ,'( ',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 + mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1 + mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1 goto 10 endif @@ -2049,7 +1802,6 @@ 10 continue enddo !end loop on clusters(X) - do icl=1,nclstr1 if(cl_single(icl).eq.1)then ip=npl(VIEW(icl)) @@ -2057,18 +1809,86 @@ cls(ip,ncls(ip))=icl endif enddo + +c 80 continue + continue if(DEBUG.EQ.1)then print*,'clusters ',nclstr1 print*,'good ',(cl_good(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 + + + if(.not.RECOVER_SINGLETS)goto 81 + +* //////////////////////////////////////////////// +* PATCH to recover events with less than 3 couples +* //////////////////////////////////////////////// +* loop over singlet and create "fake" couples +* (with clx=0 or cly=0) +* + + 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 +* ---------------------------------------------------- +* jump masked views +* ---------------------------------------------------- + if( mask_view(VIEW(icl)).ne.0 ) goto 21 + if(mod(VIEW(icl),2).eq.0)then !=== X views +* ---------------------------------------------------- +* cut on charge (X VIEW) +* ---------------------------------------------------- + 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!!! +* ---------------------------------------------------- + + else !=== Y views +* ---------------------------------------------------- +* cut on charge (Y VIEW) +* ---------------------------------------------------- + if(sgnl(icl).lt.dedx_y_min) goto 21 + + 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!!! +* ---------------------------------------------------- + + endif + endif !end "single" condition + 21 continue + enddo !end loop over clusters + + if(DEBUG.EQ.1) + $ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) + + + 81 continue - do ip=1,6 + ncp_tot=0 + do ip=1,NPLANES ncp_tot = ncp_tot + ncp_plane(ip) enddo + if(DEBUG.EQ.1) + $ print*,'n.couple tot: ',ncp_tot return end @@ -2131,8 +1951,6 @@ * -------------------------------------------- do ip=1,nplanes if(ncp_plane(ip).gt.ncouplelimit)then -c mask_view(nviewx(ip)) = 8 -c mask_view(nviewy(ip)) = 8 mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7 mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7 endif @@ -2143,27 +1961,38 @@ ntrpt=0 !number of triplets do ip1=1,(nplanes-1) !loop on planes - 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) 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=xPAM ym1=yPAM zm1=zPAM -c print*,'***',is1,xm1,ym1,zm1 do ip2=(ip1+1),nplanes !loop on planes - 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 @@ -2173,7 +2002,17 @@ xm2=xPAM ym2=yPAM zm2=zPAM - + +* --------------------------------------------------- +* 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) @@ -2183,15 +2022,19 @@ $ '** 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 +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) @@ -2200,31 +2043,37 @@ alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2) * y0 (cm) alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1 - + **** -----------------------------------------------**** **** 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(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 @@ -2232,6 +2081,23 @@ 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 @@ -2242,83 +2108,131 @@ xm3=xPAM ym3=yPAM zm3=zPAM + + * find the circle passing through the three points -c$$$ call tricircle(3,xp,zp,angp,resp,chi -c$$$ $ ,xc,zc,radius,iflag) - iflag = DEBUG + 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*ZINI+BX + + 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.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 + $ '** 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 +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-(ZINI-zc)**2) + alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) + 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-(ZINI-zc)**2) + alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) + 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 - 31 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 @@ -2390,9 +2304,6 @@ 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 @@ -2428,19 +2339,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 @@ -2456,13 +2359,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 @@ -2479,7 +2382,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 @@ -2489,9 +2394,7 @@ do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo -c print*,'>>>> ',ncpused,npt,nplused -c if(ncpused.lt.ncpyz_min)goto 2228 !next doublet - if(npt.lt.nptyz_min)goto 2228 !next doublet + if(nplused.lt.nplyz_min)goto 2228 !next doublet * ~~~~~~~~~~~~~~~~~ @@ -2523,11 +2426,9 @@ 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.EQ.1)then - print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points' print*,'- alfayz1 ',alfayz1_av(nclouds_yz) print*,'- alfayz2 ',alfayz2_av(nclouds_yz) @@ -2536,10 +2437,6 @@ print*,'cpcloud_yz ' $ ,(cpcloud_yz(nclouds_yz,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) endif * >>> NEW CLOUD <<< * ~~~~~~~~~~~~~~~~~ @@ -2554,9 +2451,7 @@ endif if(DEBUG.EQ.1)then - print*,'---------------------- ' print*,'Y-Z total clouds ',nclouds_yz - print*,' ' endif @@ -2619,8 +2514,6 @@ 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 @@ -2654,13 +2547,16 @@ 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) - + + * ------------------------------------------------------------------------ * FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE * ------------------------------------------------------------------------ @@ -2673,7 +2569,6 @@ $ .true.)istrimage=1 if(distance.lt.cutdistxz.or.istrimage.eq.1)then -c print*,idb1,idb2,distance,' cloud ',nclouds_yz 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 @@ -2692,13 +2587,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 @@ -2713,13 +2608,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 @@ -2729,8 +2625,6 @@ do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo -c if(ncpused.lt.ncpxz_min)goto 22288 !next triplet - if(npt.lt.nptxz_min)goto 22288 !next triplet if(nplused.lt.nplxz_min)goto 22288 !next triplet * ~~~~~~~~~~~~~~~~~ @@ -2764,8 +2658,7 @@ npt_tot=npt_tot+npt if(DEBUG.EQ.1)then - print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' - print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' + print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz) print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz) print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz) @@ -2774,10 +2667,6 @@ 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 <<< * ~~~~~~~~~~~~~~~~~ @@ -2791,9 +2680,7 @@ endif if(DEBUG.EQ.1)then - print*,'---------------------- ' print*,'X-Z total clouds ',nclouds_xz - print*,' ' endif @@ -2849,8 +2736,8 @@ 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 @@ -2859,28 +2746,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 + 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 -* 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 + + 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 @@ -2907,66 +2838,33 @@ do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo - - + + 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*,'Combination ',iyz,ixz - $ ,' db ',ptcloud_yz(iyz) - $ ,' tr ',ptcloud_xz(ixz) - $ ,' -----> # matching couples ',ncp_ok + 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 - -c if(nplused.lt.nplxz_min)goto 888 !next combination - if(nplused.lt.nplyz_min)goto 888 !next combination - if(ncp_ok.lt.ncpok)goto 888 !next combination - -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 <------- -cccc SBAGLIATO -c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz)) -c$$$ AL_INI(2) = dreal(alfayz1_av(iyz)) -c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz)) -c$$$ $ /dreal(alfaxz2_av(ixz))) -c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4)) -c$$$ AL_INI(3) = tath/sqrt(1+tath**2) -c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0. -cccc GIUSTO (ma si sua guess()) -c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz)) -c$$$ AL_INI(2) = dreal(alfayz1_av(iyz)) -c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4)) -c$$$ AL_INI(3) = tath/sqrt(1+tath**2) -c$$$ IF(alfaxz2_av(ixz).NE.0)THEN -c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz)) -c$$$ $ /dreal(alfaxz2_av(ixz))) -c$$$ ELSE -c$$$ AL_INI(4) = acos(-1.)/2 -c$$$ IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.) -c$$$ ENDIF -c$$$ IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4) -c$$$ AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs -c$$$ -c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0. -c$$$ -c$$$ if(AL_INI(5).gt.defmax)goto 888 !next cloud if(DEBUG.EQ.1)then - print*,'track candidate', ntracks+1 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)) @@ -2999,6 +2897,16 @@ 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 @@ -3047,13 +2955,42 @@ $ 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 @@ -3094,6 +3031,8 @@ * ********************************************************** 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 @@ -3120,7 +3059,7 @@ 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)) @@ -3139,11 +3078,19 @@ $ cp_match(ip,hit_plane(ip)) SENSOR_STORE(nplanes-ip+1,ntracks) $ = is_cp(cp_match(ip,hit_plane(ip))) - LADDER_STORE(nplanes-ip+1,ntracks) - $ = LADDER( + + 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 @@ -3177,18 +3124,9 @@ if(ntracks.eq.0)then iflag=1 - return +cc return endif -c$$$ if(DEBUG.EQ.1)then -c$$$ print*,'****** TRACK CANDIDATES ***********' -c$$$ print*,'# R. chi2 RIG' -c$$$ do i=1,ntracks -c$$$ print*,i,' --- ',rchi2_store(i),' --- ' -c$$$ $ ,1./abs(AL_STORE(5,i)) -c$$$ enddo -c$$$ print*,'***********************************' -c$$$ endif if(DEBUG.EQ.1)then print*,'****** TRACK CANDIDATES *****************' print*,'# R. chi2 RIG ndof' @@ -3250,6 +3188,8 @@ 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) @@ -3266,7 +3206,42 @@ * 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) @@ -3291,16 +3266,48 @@ $ bxyz(2) $ ) - 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 - - dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) - dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) + 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 * ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- @@ -3312,6 +3319,10 @@ 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 @@ -3334,8 +3345,7 @@ * =========================================== * STEP 1 >>>>>>> try to include a new couple * =========================================== -c if(DEBUG.EQ.1)print*,'>>>> try to include a new couple' - distmin=1000000. + distmin=100000000. xmm = 0. ymm = 0. zmm = 0. @@ -3348,11 +3358,12 @@ 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, @@ -3366,7 +3377,8 @@ distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI id=id_cp(ip,icp,ist) - if(DEBUG.EQ.1)print*,'( couple ',id + if(DEBUG.EQ.1) + $ print*,'( couple ',id $ ,' ) distance ',distance if(distance.lt.distmin)then xmm = xPAM @@ -3378,14 +3390,12 @@ idm = id dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) -c QUIQUI --> non devo moltiplicare per clinc?!?!?! - clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI - $ +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI + clincnewc=10*sqrt(rymm**2+rxmm**2 + $ +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) endif 1188 continue enddo !end loop on couples on plane icp -c if(distmin.le.clinc)then !QUIQUI - if(distmin.le.clincnewc)then !QUIQUI + if(distmin.le.clincnewc)then * ----------------------------------- xm(nplanes-ip+1) = xmm !<<< ym(nplanes-ip+1) = ymm !<<< @@ -3399,14 +3409,13 @@ * ----------------------------------- CP_STORE(nplanes-ip+1,ibest)=idm if(DEBUG.EQ.1)print*,'%%%% included couple ',idm - $ ,' (dist.= ',distmin,', cut ',clinc,' )' + $ ,' (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.EQ.1)print*,'>>>> try to include a new cluster' distmin=1000000. xmm_A = 0. !--------------------------- ymm_A = 0. ! init variables that @@ -3425,6 +3434,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 @@ -3443,7 +3453,8 @@ $ ) distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI - if(DEBUG.EQ.1)print*,'( cl-X ',icx + if(DEBUG.EQ.1) + $ print*,'( cl-X ',icx $ ,' in cp ',id,' ) distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A @@ -3476,7 +3487,8 @@ $ ) distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI - if(DEBUG.EQ.1)print*,'( cl-Y ',icy + if(DEBUG.EQ.1) + $ print*,'( cl-Y ',icy $ ,' in cp ',id,' ) distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A @@ -3496,10 +3508,8 @@ 11882 continue enddo !end loop on cluster inside couples *----- single clusters ----------------------------------------------- -c print*,'## ncls(',ip,') ',ncls(ip) 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 @@ -3521,10 +3531,10 @@ distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI - if(DEBUG.EQ.1)print*,'( cl-s ',icl + if(DEBUG.EQ.1) + $ print*,'( cl-s ',icl $ ,' ) distance ',distance if(distance.lt.distmin)then -c if(DEBUG.EQ.1)print*,'YES' xmm_A = xPAM_A ymm_A = yPAM_A zmm_A = zPAM_A @@ -3545,10 +3555,7 @@ endif 18882 continue enddo !end loop on single clusters -c print*,'## distmin ', distmin,' clinc ',clinc -c QUIQUI------------ -c anche qui: non ci vuole clinc??? if(iclm.ne.0)then if(mod(VIEW(iclm),2).eq.0)then clincnew= @@ -3559,30 +3566,28 @@ $ 10* $ sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2)) endif -c QUIQUI------------ - - if(distmin.le.clincnew)then !QUIQUI -c if(distmin.le.clinc)then !QUIQUI + + if(distmin.le.clincnew)then CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< * ---------------------------- -c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' if(mod(VIEW(iclm),2).eq.0)then XGOOD(nplanes-ip+1)=1. resx(nplanes-ip+1)=rxmm - if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm + if(DEBUG.EQ.1) + $ print*,'%%%% included X-cl ',iclm $ ,'( chi^2, ',RCHI2_STORE(ibest) $ ,', dist.= ',distmin - $ ,', cut ',clinc,' )' + $ ,', cut ',clincnew,' )' else YGOOD(nplanes-ip+1)=1. resy(nplanes-ip+1)=rymm - if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm + if(DEBUG.EQ.1) + $ print*,'%%%% included Y-cl ',iclm $ ,'( chi^2, ',RCHI2_STORE(ibest) $ ,', dist.= ', distmin - $ ,', cut ',clinc,' )' + $ ,', cut ',clincnew,' )' endif -c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' * ---------------------------- xm_A(nplanes-ip+1) = xmm_A ym_A(nplanes-ip+1) = ymm_A @@ -3603,6 +3608,7 @@ return end + *************************************************** * * * * @@ -3612,78 +3618,6 @@ * * ************************************************** * - subroutine clean_XYclouds(ibest,iflag) - - include 'commontracker.f' - include 'level1.f' - include 'common_momanhough.f' - include 'level2.f' - - if(DEBUG.EQ.1)print*,'clean_XYclouds:' - - 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)=ntrk !tag used clusters -c$$$ cl_used(icly)=ntrk !tag used clusters - elseif(icl.ne.0)then -c$$$ cl_used(icl)=ntrk !tag used clusters - 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.EQ.1)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 - - - @@ -3732,6 +3666,12 @@ DEDX_Y(IP,IT) = 0 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 @@ -3751,6 +3691,10 @@ ys(1,ip)=0 ys(2,ip)=0 sgnlys(ip)=0 + sxbad(ip)=0 + sybad(ip)=0 + multmaxsx(ip)=0 + multmaxsy(ip)=0 enddo end @@ -3873,6 +3817,8 @@ integer ssensor,sladder pig=acos(-1.) + + * ------------------------------------- chi2_nt(ntr) = sngl(chi2) nstep_nt(ntr) = nstep @@ -3920,18 +3866,24 @@ dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)/factor) dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)/factor) + +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) - if(ip.eq.6) ax = -1. * axv_nt(ip,ntr) - if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND) - tgtemp = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001 - angx = 180.*atan(tgtemp)/acos(-1.) - tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001 - angy = 180.*atan(tgtemp)/acos(-1.) +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) + -c print*,'* ',ip,bfx,bfy,angx,angy id = CP_STORE(ip,IDCAND) ! couple id icl = CLS_STORE(ip,IDCAND) @@ -3942,44 +3894,91 @@ if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align LS(IP,ntr) = ssensor+10*sladder - if(id.ne.0)then +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)) - cl_used(cltrx(ip,ntr)) = 1 !tag used clusters - cl_used(cltry(ip,ntr)) = 1 !tag used clusters + 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))) -c$$$ nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx) -c$$$ nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy) -c$$$ xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id))) -c$$$ ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id))) - xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id))) - ybad(ip,ntr)= nbadstrips(4,cly(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) - if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0) - $ dedx_y(ip,ntr)=-dedx_y(ip,ntr) + 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 - elseif(icl.ne.0)then +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 -c$$$ nnnnn = npfastrips(icl,PFA,angx) -c$$$ xbad(ip,ntr) = nbadstrips(nnnnn,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 -c$$$ nnnnn = npfastrips(icl,PFA,angy) -c$$$ ybad(ip,ntr) = nbadstrips(nnnnn,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 @@ -3992,15 +3991,12 @@ 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 -c$$$ print*,(xgood(i),i=1,6) -c$$$ print*,(ygood(i),i=1,6) -c$$$ print*,(ls(i,ntr),i=1,6) -c$$$ print*,(dedx_x(i,ntr),i=1,6) -c$$$ print*,(dedx_y(i,ntr),i=1,6) -c$$$ print*,'-----------------------' - end subroutine fill_level2_siglets @@ -4036,40 +4032,40 @@ ip=nplanes-npl(VIEW(icl))+1 if(cl_used(icl).eq.0)then !cluster not included in any track + if(mod(VIEW(icl),2).eq.0)then !=== X views + nclsx = nclsx + 1 planex(nclsx) = ip 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.) c call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.) call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.) xs(is,nclsx) = (xPAM_A+xPAM_B)/2 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) = 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.) c call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.) call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.) ys(is,nclsy) = (yPAM_A+yPAM_B)/2 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 @@ -4082,24 +4078,7 @@ * associati ad una traccia, e permettere di salvare * solo questi nell'albero di uscita * -------------------------------------------------- - - -c$$$ print*,' cl ',icl,' --> ',cl_used(icl) -c$$$ -c$$$ if( cl_used(icl).ne.0 )then -c$$$ if( -c$$$ $ mod(VIEW(icl),2).eq.0.and. -c$$$ $ cltrx(ip,whichtrack(icl)).ne.icl ) -c$$$ $ print*,'**WARNING** cltrx(',ip,',',whichtrack(icl) -c$$$ $ ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl -c$$$ if( -c$$$ $ mod(VIEW(icl),2).eq.1.and. -c$$$ $ cltry(ip,whichtrack(icl)).ne.icl ) -c$$$ $ print*,'**WARNING** cltry(',ip,',',whichtrack(icl) -c$$$ $ ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl -c$$$ endif - - + enddo end @@ -4161,7 +4140,7 @@ alfayz2_av_nt(iyz)=alfayz2_av(iyz) nnn=nnn+ptcloud_yz(iyz) enddo - do ipt=1,nnn + do ipt=1,min(ndblt_max_nt,nnn) db_cloud_nt(ipt)=db_cloud(ipt) enddo endif @@ -4174,7 +4153,7 @@ alfaxz3_av_nt(ixz)=alfaxz3_av(ixz) nnn=nnn+ptcloud_xz(ixz) enddo - do ipt=1,nnn + do ipt=1,min(ntrpt_max_nt,nnn) tr_cloud_nt(ipt)=tr_cloud(ipt) enddo endif