--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/05/14 11:03:06 1.23 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/08/31 14:56:52 1.31 @@ -215,6 +215,7 @@ c include 'momanhough_init.f' logical FIMAGE ! + real trackimage(NTRACKSMAX) real*8 AL_GUESS(5) *------------------------------------------------------------------------------- @@ -331,7 +332,7 @@ iimage=0 endif if(icand.eq.0)then - if(VERBOSE)then + if(VERBOSE.EQ.1)then print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand $ ,ibest,iimage endif @@ -360,12 +361,12 @@ jstep=0 !# minimization steps iprint=0 -c if(DEBUG)iprint=1 - if(VERBOSE)iprint=1 - if(DEBUG)iprint=2 +c if(DEBUG.EQ.1)iprint=1 + if(VERBOSE.EQ.1)iprint=1 + if(DEBUG.EQ.1)iprint=2 call mini2(jstep,ifail,iprint) if(ifail.ne.0) then - if(VERBOSE)then + if(VERBOSE.EQ.1)then print *, $ '*** MINIMIZATION FAILURE *** (after refinement) ' $ ,iev @@ -380,7 +381,7 @@ c 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 @@ -391,7 +392,7 @@ endif c rchi2=chi2/dble(ndof) - if(DEBUG)then + if(DEBUG.EQ.1)then print*,' ' print*,'****** SELECTED TRACK *************' print*,'# R. chi2 RIG' @@ -407,25 +408,140 @@ * ------------- 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 +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, +* 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 +c$$$ print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i) +c$$$ $ ,CP_STORE(nplanes-ip+1,icand),ncp enddo - if( iimage.ne.0.and. -c $ RCHI2_STORE(i).le.CHI2MAX.and. -c $ RCHI2_STORE(i).gt.0.and. - $ .true.)then - if(DEBUG)print*,'Track candidate ',iimage - $ ,' >>> TRACK IMAGE >>> of' - $ ,ibest - goto 122 !image track found + 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 there are, choose the best one + if(ngood.gt.1)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 + + 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 @@ -438,7 +554,7 @@ c $ ,iimage,fimage,ntrk,image(ntrk) if(ntrk.eq.NTRKMAX)then - if(verbose) + if(verbose.eq.1) $ print*, $ '** warning ** number of identified '// $ 'tracks exceeds vector dimension ' @@ -474,7 +590,7 @@ $ rchi2best.le.CHI2MAX.and. c $ rchi2best.lt.15..and. $ .true.)then - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'***** NEW SEARCH ****' endif goto 11111 !try new search @@ -596,130 +712,50 @@ zPAM_B = 0. c print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy + 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 +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 * -------------------------- - 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 = risx_eta2(abs(angx)) - resxPAM = resxPAM*fbad_cog(2,icx) - if(DEBUG.and.fbad_cog(2,icx).ne.1) - $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) - - elseif(PFAx.eq.'ETA3')then - - stripx = stripx + pfaeta3(icx,angx) - resxPAM = risx_eta3(abs(angx)) - resxPAM = resxPAM*fbad_cog(3,icx) - if(DEBUG.and.fbad_cog(3,icx).ne.1) - $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx) - - elseif(PFAx.eq.'ETA4')then - - stripx = stripx + pfaeta4(icx,angx) - resxPAM = risx_eta4(abs(angx)) - resxPAM = resxPAM*fbad_cog(4,icx) - if(DEBUG.and.fbad_cog(4,icx).ne.1) - $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx) - - elseif(PFAx.eq.'ETA')then - - stripx = stripx + pfaeta(icx,angx) - resxPAM = ris_eta(icx,angx) - resxPAM = resxPAM*fbad_eta(icx,angx) - if(DEBUG.and.fbad_cog(2,icx).ne.1) - $ 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) 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=cog(4,icx) -c$$$ print*,icx,' *** ',cc -c$$$ print*,icx,' *** ',resxPAM - endif - - endif - + 10 endif + * ----------------- * CLUSTER Y * ----------------- @@ -729,113 +765,42 @@ viewy = VIEW(icy) nldy = nld(MAXS(icy),VIEW(icy)) nply = npl(VIEW(icy)) - resyPAM = RESYAV +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 - if(DEBUG) then + if(DEBUG.EQ.1) then print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' $ ,icx,icy 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 = risy_eta2(abs(angy)) - 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 - - stripy = stripy + pfaeta3(icy,angy) - resyPAM = resyPAM*fbad_cog(3,icy) - if(DEBUG.and.fbad_cog(3,icy).ne.1) - $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy) + call applypfa(PFAy,icy,angy,corr,res) + stripy = stripy + corr + resyPAM = res - elseif(PFAy.eq.'ETA4')then - - stripy = stripy + pfaeta4(icy,angy) - resyPAM = resyPAM*fbad_cog(4,icy) - if(DEBUG.and.fbad_cog(4,icy).ne.1) - $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy) - - elseif(PFAy.eq.'ETA')then - - stripy = stripy + pfaeta(icy,angy) - resyPAM = ris_eta(icy,angy) - resyPAM = resyPAM*fbad_eta(icy,angy) - if(DEBUG.and.fbad_cog(2,icy).ne.1) - $ 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) 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=cog(4,icy) -c$$$ print*,icy,' *** ',cc -c$$$ print*,icy,' *** ',resyPAM - endif - - - endif + 20 endif c$$$ print*,'## stripx,stripy ',stripx,stripy @@ -849,7 +814,7 @@ 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... - if(DEBUG) then + if(DEBUG.EQ.1) then print*,'xyz_PAM (couple):', $ ' WARNING: false X strip: strip ',stripx endif @@ -858,7 +823,6 @@ yi = acoordsi(stripy,viewy) zi = 0. - c------------------------------------------------------------------------ c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c------------------------------------------------------------------------ @@ -944,7 +908,7 @@ 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) then + if(DEBUG.EQ.1) then print*,'xyz_PAM (X-singlet):', $ ' WARNING: false X strip: strip ',stripx endif @@ -969,7 +933,7 @@ c print*,yi_A,' <--> ',yi_B else - if(DEBUG) then + if(DEBUG.EQ.1) then print *,'routine xyz_PAM ---> not properly used !!!' print *,'icx = ',icx print *,'icy = ',icy @@ -1038,7 +1002,7 @@ c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')' else - if(DEBUG) then + if(DEBUG.EQ.1) then print *,'routine xyz_PAM ---> not properly used !!!' print *,'icx = ',icx print *,'icy = ',icy @@ -1082,11 +1046,23 @@ c$$$ PFAx = 'COG4'!PFA c$$$ PFAy = 'COG4'!PFA + + + if(icx.gt.nclstr1.or.icy.gt.nclstr1)then + print*,'xyzpam: ***WARNING*** clusters ',icx,icy + $ ,' does not exists (nclstr1=',nclstr1,')' + icx = -1*icx + icy = -1*icy + return + + endif call idtoc(pfaid,PFAx) call idtoc(pfaid,PFAy) - call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) +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 @@ -1094,9 +1070,24 @@ ipx=npl(VIEW(icx)) ipy=npl(VIEW(icy)) - if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip ) - $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy - $ ,' does not belong to the correct plane: ',ip,ipx,ipy +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 + $ ,' 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. @@ -1116,10 +1107,18 @@ elseif(icx.eq.0.and.icy.ne.0)then ipy=npl(VIEW(icy)) - if((nplanes-ipy+1).ne.ip) - $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy - $ ,' does not belong to the correct plane: ',ip,ipx,ipy +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 + 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. @@ -1138,10 +1137,19 @@ elseif(icx.ne.0.and.icy.eq.0)then ipx=npl(VIEW(icx)) - if((nplanes-ipx+1).ne.ip) - $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy - $ ,' does not belong to the correct plane: ',ip,ipx,ipy +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 + $ ,' 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 @@ -1182,7 +1190,7 @@ endif - if(DEBUG)then + 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) @@ -1633,6 +1641,8 @@ integer iflag integer badseed,badclx,badcly + + if(DEBUG.EQ.1)print*,'cl_to_couples:' * init variables ncp_tot=0 @@ -1663,9 +1673,11 @@ * mask views with too many clusters do iv=1,nviews if( ncl_view(iv).gt. nclusterlimit)then - mask_view(iv) = 1 - if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > ' - $ ,nclusterlimit,' on view ', iv,' --> masked!' +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 @@ -1803,13 +1815,15 @@ endif if(ncp_plane(nplx).gt.ncouplemax)then - if(verbose)print*, + if(verbose.eq.1)print*, $ '** warning ** number of identified '// $ 'couples on plane ',nplx, $ 'exceeds vector dimention ' $ ,'( ',ncouplemax,' ) --> masked!' - mask_view(nviewx(nplx)) = 2 - mask_view(nviewy(nply)) = 2 +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 @@ -1839,10 +1853,10 @@ enddo - 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*,'singlets',(cl_single(i),i=1,nclstr1) print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) endif @@ -1902,6 +1916,7 @@ real xc,zc,radius * ----------------------------- + if(DEBUG.EQ.1)print*,'cp_to_doubtrip:' * -------------------------------------------- * put a limit to the maximum number of couples @@ -1910,8 +1925,10 @@ * -------------------------------------------- do ip=1,nplanes if(ncp_plane(ip).gt.ncouplelimit)then - mask_view(nviewx(ip)) = 8 - mask_view(nviewy(ip)) = 8 +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 enddo @@ -1937,8 +1954,7 @@ do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2 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 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) @@ -1957,14 +1973,15 @@ * (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 do iv=1,12 - mask_view(iv) = 3 +c mask_view(iv) = 3 + mask_view(iv) = mask_view(iv)+ 2**2 enddo iflag=1 return @@ -2020,6 +2037,9 @@ 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 call tricircle(3,xp,zp,angp,resp,chi $ ,xc,zc,radius,iflag) c print*,xc,zc,radius @@ -2036,14 +2056,15 @@ * (3 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ntrpt.eq.ntrpt_max)then - if(verbose)print*, + 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 - mask_view(iv) = 4 +c mask_view(iv) = 4 + mask_view(iv)=mask_view(iv)+ 2**3 enddo iflag=1 return @@ -2097,7 +2118,7 @@ 10 continue enddo !end loop on planes - COPPIA 1 - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'--- doublets ',ndblt print*,'--- triplets ',ntrpt endif @@ -2144,6 +2165,7 @@ integer cp_useds1(ncouplemaxtot) ! sensor 1 integer cp_useds2(ncouplemaxtot) ! sensor 2 + if(DEBUG.EQ.1)print*,'doub_to_YZcloud:' *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * classification of DOUBLETS @@ -2270,14 +2292,15 @@ * >>> 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 - mask_view(iv) = 5 +c mask_view(iv) = 5 + mask_view(iv) = mask_view(iv) + 2**4 enddo iflag=1 return @@ -2297,14 +2320,16 @@ c print*,'>> ',ipt,db_all(ipt) enddo npt_tot=npt_tot+npt - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' 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) + 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) c$$$ print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = ' c$$$ $ ,ptcloud_yz(nclouds_yz) c$$$ print*,'nt-uple: db_cloud(...) = ' @@ -2322,7 +2347,7 @@ goto 90 endif - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'---------------------- ' print*,'Y-Z total clouds ',nclouds_yz print*,' ' @@ -2371,6 +2396,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 @@ -2428,7 +2455,18 @@ $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2 distance = sqrt(distance) - if(distance.lt.cutdistxz)then +* ------------------------------------------------------------------------ +* 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 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 @@ -2492,14 +2530,15 @@ * ~~~~~~~~~~~~~~~~~ * >>> 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 - mask_view(iv) = 6 +c mask_view(iv) = 6 + mask_view(iv) = mask_view(iv) + 2**5 enddo iflag=1 return @@ -2518,14 +2557,16 @@ enddo npt_tot=npt_tot+npt - if(DEBUG)then + if(DEBUG.EQ.1)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) + 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) @@ -2543,7 +2584,7 @@ goto 91 endif - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'---------------------- ' print*,'X-Z total clouds ',nclouds_xz print*,' ' @@ -2595,10 +2636,9 @@ * ----------------------------------------------------------- * variables for track fitting double precision AL_INI(5) -c 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 @@ -2620,7 +2660,7 @@ enddo enddo ncp_ok=0 - do icp=1,ncp_tot !loop on couples + do icp=1,ncp_tot !loop over couples * get info on cpintersec(icp)=min( $ cpcloud_yz(iyz,icp), @@ -2629,6 +2669,7 @@ $ (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 +* cpintersec is >0 if yz and xz clouds contain the same image of couple icp if(cpintersec(icp).ne.0)then ncp_ok=ncp_ok+1 @@ -2661,16 +2702,18 @@ nplused=nplused+ hit_plane(ip) enddo -c if(nplused.lt.nplxz_min)goto 888 !next doublet - if(nplused.lt.nplyz_min)goto 888 !next doublet - if(ncp_ok.lt.ncpok)goto 888 !next cloud - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'Combination ',iyz,ixz $ ,' db ',ptcloud_yz(iyz) $ ,' tr ',ptcloud_xz(ixz) $ ,' -----> # matching couples ',ncp_ok 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)) @@ -2716,7 +2759,8 @@ c$$$ c$$$ if(AL_INI(5).gt.defmax)goto 888 !next cloud - if(DEBUG)then + 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)) @@ -2749,10 +2793,35 @@ hit_plane(6)=icp6 if(ncp_match(6).eq.0)hit_plane(6)=0 !-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) @@ -2796,11 +2865,11 @@ ifail=0 !error flag in chi^2 computation jstep=0 !number of minimization steps iprint=0 -c if(DEBUG)iprint=1 - if(DEBUG)iprint=2 +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 *** ' $ //'(clouds_to_ctrack)' @@ -2825,14 +2894,15 @@ * -------------------------- 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 - mask_view(iv) = 7 +c mask_view(iv) = 7 + mask_view(iv) = mask_view(iv) + 2**6 enddo iflag=1 return @@ -2840,7 +2910,7 @@ ntracks = ntracks + 1 - do ip=1,nplanes + do ip=1,nplanes !top to bottom XV_STORE(ip,ntracks)=sngl(xv(ip)) YV_STORE(ip,ntracks)=sngl(yv(ip)) @@ -2857,6 +2927,7 @@ 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)) @@ -2872,9 +2943,9 @@ SENSOR_STORE(nplanes-ip+1,ntracks)=0 LADDER_STORE(nplanes-ip+1,ntracks)=0 endif - BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now - BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now - 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 @@ -2903,7 +2974,7 @@ return endif -c$$$ if(DEBUG)then +c$$$ if(DEBUG.EQ.1)then c$$$ print*,'****** TRACK CANDIDATES ***********' c$$$ print*,'# R. chi2 RIG' c$$$ do i=1,ntracks @@ -2912,7 +2983,7 @@ c$$$ enddo c$$$ print*,'***********************************' c$$$ endif - if(DEBUG)then + if(DEBUG.EQ.1)then print*,'****** TRACK CANDIDATES *****************' print*,'# R. chi2 RIG ndof' do i=1,ntracks @@ -2964,6 +3035,7 @@ 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 @@ -3045,7 +3117,7 @@ LADDER_STORE(nplanes-ip+1,IBEST)=nldt * -------------------------------------------------------------- - if(DEBUG)then + if(DEBUG.EQ.1)then print*, $ '------ Plane ',ip,' intersected on LADDER ',nldt $ ,' SENSOR ',ist @@ -3056,7 +3128,7 @@ * =========================================== * STEP 1 >>>>>>> try to include a new couple * =========================================== -c if(DEBUG)print*,'>>>> try to include a new couple' +c if(DEBUG.EQ.1)print*,'>>>> try to include a new couple' distmin=1000000. xmm = 0. ymm = 0. @@ -3088,7 +3160,7 @@ distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI id=id_cp(ip,icp,ist) - if(DEBUG)print*,'( couple ',id + if(DEBUG.EQ.1)print*,'( couple ',id $ ,' ) distance ',distance if(distance.lt.distmin)then xmm = xPAM @@ -3120,7 +3192,7 @@ dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< * ----------------------------------- CP_STORE(nplanes-ip+1,ibest)=idm - if(DEBUG)print*,'%%%% included couple ',idm + if(DEBUG.EQ.1)print*,'%%%% included couple ',idm $ ,' (dist.= ',distmin,', cut ',clinc,' )' goto 133 !next plane endif @@ -3128,7 +3200,7 @@ * STEP 2 >>>>>>> try to include a single cluster * either from a couple or single * ================================================ -c if(DEBUG)print*,'>>>> try to include a new cluster' +c if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster' distmin=1000000. xmm_A = 0. !--------------------------- ymm_A = 0. ! init variables that @@ -3165,7 +3237,7 @@ $ ) distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI - if(DEBUG)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 @@ -3198,7 +3270,7 @@ $ ) distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI - if(DEBUG)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 @@ -3243,10 +3315,10 @@ distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI - if(DEBUG)print*,'( cl-s ',icl - $ ,' ) distance ',distance,'<',distmin,' ?' + if(DEBUG.EQ.1)print*,'( cl-s ',icl + $ ,' ) distance ',distance if(distance.lt.distmin)then - if(DEBUG)print*,'YES' +c if(DEBUG.EQ.1)print*,'YES' xmm_A = xPAM_A ymm_A = yPAM_A zmm_A = zPAM_A @@ -3292,14 +3364,14 @@ 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 + if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm $ ,'( chi^2, ',RCHI2_STORE(ibest) $ ,', dist.= ',distmin $ ,', cut ',clinc,' )' else YGOOD(nplanes-ip+1)=1. resy(nplanes-ip+1)=rymm - if(DEBUG)print*,'%%%% included Y-cl ',iclm + if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm $ ,'( chi^2, ',RCHI2_STORE(ibest) $ ,', dist.= ', distmin $ ,', cut ',clinc,' )' @@ -3341,6 +3413,8 @@ 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) @@ -3349,10 +3423,10 @@ if(id.ne.0)then iclx=clx(ip,icp_cp(id)) icly=cly(ip,icp_cp(id)) - cl_used(iclx)=ntrk !tag used clusters - cl_used(icly)=ntrk !tag used clusters +c$$$ cl_used(iclx)=ntrk !tag used clusters +c$$$ cl_used(icly)=ntrk !tag used clusters elseif(icl.ne.0)then - cl_used(icl)=ntrk !tag used clusters +c$$$ cl_used(icl)=ntrk !tag used clusters endif * ----------------------------- @@ -3371,7 +3445,7 @@ $ cly(ip,icp).eq.icl $ )then id=id_cp(ip,icp,1) - if(DEBUG)then + if(DEBUG.EQ.1)then print*,ip,' <<< cp ',id $ ,' ( cl-x ' $ ,clx(ip,icp) @@ -3452,6 +3526,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 @@ -3631,12 +3711,12 @@ zv_nt(ip,ntr) = sngl(zv(ip)) axv_nt(ip,ntr) = sngl(axv(ip)) ayv_nt(ip,ntr) = sngl(ayv(ip)) -c l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + factor = sqrt( $ tan( acos(-1.) * sngl(axv(ip)) /180. )**2 + $ tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 + $ 1. ) -c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)/factor) dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)/factor) @@ -3644,12 +3724,16 @@ 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 @@ -3667,43 +3751,77 @@ cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id)) cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id)) -c$$$ if(is_cp(id).ne.ssensor) -c$$$ $ print*,'ERROR is sensor assignment (couple)' -c$$$ $ ,is_cp(id),ssensor -c$$$ if(LADDER(clx(nplanes-ip+1,icp_cp(id))).ne.sladder) -c$$$ $ print*,'ERROR is ladder assignment (couple)' -c$$$ $ ,LADDER(clx(nplanes-ip+1,icp_cp(id))),sladder - - nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx) - nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy) - xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id))) - ybad(ip,ntr)= nbadstrips(nnnny,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 + + 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) + 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 + + 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 + elseif(icl.ne.0)then -c >>> is a singlet -c$$$ if(LADDER(icl).ne.sladder) -c$$$ $ print*,'ERROR is ladder assignment (single)' -c$$$ $ ,LADDER(icl),sladder + + cl_used(icl) = 1 !tag used clusters + if(mod(VIEW(icl),2).eq.0)then cltrx(ip,ntr)=icl - nnnnn = npfastrips(icl,PFA,angx) - 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 - nnnnn = npfastrips(icl,PFA,angy) - 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 enddo + if(DEBUG.eq.1)then + print*,'> STORING TRACK ',ntr + print*,'clusters: ' + do ip=1,6 + print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr) + enddo + endif c$$$ print*,(xgood(i),i=1,6) c$$$ print*,(ygood(i),i=1,6) @@ -3734,12 +3852,19 @@ nclsy = 0 do iv = 1,nviews - if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv) +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 @@ -3779,6 +3904,30 @@ ***** 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 +* -------------------------------------------------- + + +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