--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/05/09 07:50:58 1.22 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/08/17 14:36:05 1.27 @@ -578,6 +578,9 @@ parameter (ndivx=30) + + +c$$$ print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy resxPAM = 0 resyPAM = 0 @@ -593,10 +596,16 @@ 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) @@ -604,6 +613,21 @@ nplx = npl(VIEW(icx)) 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 * -------------------------- @@ -664,7 +688,7 @@ elseif(PFAx.eq.'ETA2')then stripx = stripx + pfaeta2(icx,angx) - resxPAM = risx_eta2(abs(angx)) + resxPAM = risxeta2(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) @@ -672,26 +696,35 @@ elseif(PFAx.eq.'ETA3')then stripx = stripx + pfaeta3(icx,angx) - resxPAM = risx_eta3(abs(angx)) + resxPAM = risxeta3(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) +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 = risx_eta4(abs(angx)) + resxPAM = risxeta4(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) +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) - resxPAM = ris_eta(icx,angx) +c resxPAM = riseta(icx,angx) + resxPAM = riseta(viewx,angx) resxPAM = resxPAM*fbad_eta(icx,angx) - if(DEBUG.and.fbad_cog(2,icx).ne.1) - $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) +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 @@ -715,8 +748,9 @@ c$$$ print*,icx,' *** ',resxPAM endif - endif - + 10 endif + + * ----------------- * CLUSTER Y * ----------------- @@ -729,6 +763,20 @@ 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 print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' @@ -782,32 +830,41 @@ elseif(PFAy.eq.'ETA2')then stripy = stripy + pfaeta2(icy,angy) - resyPAM = risy_eta2(abs(angy)) + resyPAM = risyeta2(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) +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) - if(DEBUG.and.fbad_cog(3,icy).ne.1) - $ print*,'BAD icy >>> ',viewy,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) - if(DEBUG.and.fbad_cog(4,icy).ne.1) - $ print*,'BAD icy >>> ',viewy,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) - resyPAM = ris_eta(icy,angy) +c resyPAM = riseta(icy,angy) + resyPAM = riseta(viewy,angy) resyPAM = resyPAM*fbad_eta(icy,angy) - if(DEBUG.and.fbad_cog(2,icy).ne.1) - $ print*,'BAD icy >>> ',viewy,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.'ETAL')then + + 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 @@ -832,9 +889,9 @@ endif - endif + 20 endif -c print*,'## stripx,stripy ',stripx,stripy +c$$$ print*,'## stripx,stripy ',stripx,stripy c=========================================================== C COUPLE @@ -1050,6 +1107,185 @@ 100 continue end +************************************************************************ +* Call xyz_PAM subroutine with default PFA and fill the mini2 common. +* (done to be called from c/c++) +************************************************************************ + + subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy) + + include 'commontracker.f' + include 'level1.f' + include 'common_mini_2.f' + include 'common_xyzPAM.f' + include 'common_mech.f' + include 'calib.f' + +* flag to chose PFA +c$$$ character*10 PFA +c$$$ common/FINALPFA/PFA + + integer icx,icy !X-Y cluster ID + integer sensor + character*4 PFAx,PFAy !PFA to be used + real ax,ay !X-Y geometric angle + real bfx,bfy !X-Y b-field components + + ipx=0 + ipy=0 + +c$$$ PFAx = 'COG4'!PFA +c$$$ PFAy = 'COG4'!PFA + + if(icx.gt.nclstr1.or.icy.gt.nclstr1)then + print*,'xyzpam: ***WARNING*** clusters ',icx,icy + $ ,' does not exists (nclstr1=',nclstr1,')' + icx = -1*icx + icy = -1*icy + return + + endif + + call idtoc(pfaid,PFAx) + call idtoc(pfaid,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 + $ ,' does not belong to plane: ',ip + icx = -1*icx + return + endif + if( (nplanes-ipy+1).ne.ip )then + print*,'xyzpam: ***WARNING*** cluster ',icy + $ ,' does not belong to plane: ',ip + icy = -1*icy + return + endif + + call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) + + xgood(ip) = 1. + ygood(ip) = 1. + resx(ip) = resxPAM + resy(ip) = resyPAM + + xm(ip) = xPAM + ym(ip) = yPAM + zm(ip) = zPAM + xm_A(ip) = 0. + ym_A(ip) = 0. + xm_B(ip) = 0. + ym_B(ip) = 0. + +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 + icy = -1*icy + return + endif + + call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) + + xgood(ip) = 0. + ygood(ip) = 1. + resx(ip) = 1000. + resy(ip) = resyPAM + + xm(ip) = -100. + ym(ip) = -100. + zm(ip) = (zPAM_A+zPAM_B)/2. + xm_A(ip) = xPAM_A + ym_A(ip) = yPAM_A + xm_B(ip) = xPAM_B + ym_B(ip) = yPAM_B + +c zv(ip) = (zPAM_A+zPAM_B)/2. + + elseif(icx.ne.0.and.icy.eq.0)then + + ipx=npl(VIEW(icx)) +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 + resy(ip) = 1000. + + xm(ip) = -100. + ym(ip) = -100. + zm(ip) = (zPAM_A+zPAM_B)/2. + xm_A(ip) = xPAM_A + ym_A(ip) = yPAM_A + xm_B(ip) = xPAM_B + ym_B(ip) = yPAM_B + +c zv(ip) = (zPAM_A+zPAM_B)/2. + + else + + il = 2 + if(lad.ne.0)il=lad + is = 1 + if(sensor.ne.0)is=sensor +c print*,nplanes-ip+1,il,is + + xgood(ip) = 0. + ygood(ip) = 0. + resx(ip) = 1000. + resy(ip) = 1000. + + xm(ip) = -100. + ym(ip) = -100. + zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 + xm_A(ip) = 0. + ym_A(ip) = 0. + xm_B(ip) = 0. + ym_B(ip) = 0. + +c zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 + + endif + + if(DEBUG)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 ******************************************************************************** ******************************************************************************** @@ -1522,7 +1758,8 @@ * mask views with too many clusters do iv=1,nviews if( ncl_view(iv).gt. nclusterlimit)then - mask_view(iv) = 1 +c mask_view(iv) = 1 + mask_view(iv) = mask_view(iv) + 2**0 if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > ' $ ,nclusterlimit,' on view ', iv,' --> masked!' endif @@ -1667,8 +1904,10 @@ $ '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 @@ -1769,8 +2008,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 @@ -1823,7 +2064,8 @@ 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 @@ -1902,7 +2144,8 @@ 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 @@ -2136,7 +2379,8 @@ 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 @@ -2358,7 +2602,8 @@ 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 @@ -2691,7 +2936,8 @@ 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 @@ -3525,18 +3771,14 @@ c >>> is a couple cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id)) cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id)) + +c$$$ 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))) -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))) if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0) $ dedx_x(ip,ntr)=-dedx_x(ip,ntr) @@ -3544,21 +3786,22 @@ $ dedx_y(ip,ntr)=-dedx_y(ip,ntr) 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 + if(mod(VIEW(icl),2).eq.0)then cltrx(ip,ntr)=icl - nnnnn = npfastrips(icl,PFA,angx) - xbad(ip,ntr) = nbadstrips(nnnnn,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) elseif(mod(VIEW(icl),2).eq.1)then cltry(ip,ntr)=icl - nnnnn = npfastrips(icl,PFA,angy) - ybad(ip,ntr) = nbadstrips(nnnnn,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) endif + endif enddo @@ -3593,7 +3836,8 @@ 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 do icl=1,nclstr1