--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/08/28 13:25:50 1.30 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2008/11/25 14:41:38 1.36 @@ -688,6 +688,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 @@ -701,16 +704,16 @@ 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 +cc 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 ' @@ -727,7 +730,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 +750,14 @@ * -------------------------- * 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) - - 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 + call applypfa(PFAx,icx,angx,corr,res) + stripx = stripx + corr + resxPAM = res 10 endif - * ----------------- * CLUSTER Y @@ -876,7 +768,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 +792,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 - - 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 - + call applypfa(PFAy,icy,angy,corr,res) + stripy = stripy + corr + resyPAM = res 20 endif -c$$$ print*,'## stripx,stripy ',stripx,stripy +cc print*,'## stripx,stripy ',stripx,stripy c=========================================================== C COUPLE @@ -1024,11 +822,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 +859,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 +885,7 @@ nldx = nldy viewx = viewy + 1 - yi = acoordsi(stripy,viewy) + yi = dcoordsi(stripy,viewy) xi_A = edgeY_d - SiDimX/2 yi_A = yi @@ -1119,7 +916,7 @@ $ ' WARNING: false X strip: strip ',stripx endif endif - xi = acoordsi(stripx,viewx) + xi = dcoordsi(stripx,viewx) xi_A = xi yi_A = edgeX_d - SiDimY/2 @@ -1192,9 +989,9 @@ c in PAMELA reference system c------------------------------------------------------------------------ - xPAM = 0. - yPAM = 0. - zPAM = 0. + xPAM = 0.D0 + yPAM = 0.D0 + zPAM = 0.D0 xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4 yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4 @@ -1256,7 +1053,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 @@ -1303,10 +1100,10 @@ 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 @@ -1330,13 +1127,25 @@ resx(ip) = 1000. resy(ip) = resyPAM +cPP --- old --- +c$$$ xm(ip) = -100. +c$$$ ym(ip) = -100. +c$$$ zm(ip) = (zPAM_A+zPAM_B)/2. +c$$$ xm_A(ip) = xPAM_A +c$$$ ym_A(ip) = yPAM_A +c$$$ xm_B(ip) = xPAM_B +c$$$ ym_B(ip) = yPAM_B +cPP --- new --- xm(ip) = -100. ym(ip) = -100. - zm(ip) = (zPAM_A+zPAM_B)/2. + zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 xm_A(ip) = xPAM_A ym_A(ip) = yPAM_A + zm_A(ip) = zPAM_A xm_B(ip) = xPAM_B ym_B(ip) = yPAM_B + zm_B(ip) = zPAM_B +cPP ----------- c zv(ip) = (zPAM_A+zPAM_B)/2. @@ -1361,14 +1170,26 @@ resx(ip) = resxPAM resy(ip) = 1000. +cPP --- old --- +c$$$ xm(ip) = -100. +c$$$ ym(ip) = -100. +c$$$ zm(ip) = (zPAM_A+zPAM_B)/2. +c$$$ xm_A(ip) = xPAM_A +c$$$ ym_A(ip) = yPAM_A +c$$$ xm_B(ip) = xPAM_B +c$$$ ym_B(ip) = yPAM_B +cPP --- new --- xm(ip) = -100. ym(ip) = -100. - zm(ip) = (zPAM_A+zPAM_B)/2. + zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 xm_A(ip) = xPAM_A ym_A(ip) = yPAM_A + zm_A(ip) = zPAM_A xm_B(ip) = xPAM_B ym_B(ip) = yPAM_B - + zm_B(ip) = zPAM_B +cPP ----------- + c zv(ip) = (zPAM_A+zPAM_B)/2. else @@ -1426,7 +1247,7 @@ * ******************************************************************************** - real function distance_to(XPP,YPP) + real function distance_to(rXPP,rYPP) include 'common_xyzPAM.f' @@ -1435,9 +1256,14 @@ * ( 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. @@ -1564,17 +1390,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 +1415,15 @@ 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(((mod(int(stripx+0.5)-1,1024)+1).le.3) +c $ .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. +c endif + xi = dcoordsi(stripx,viewx) + yi = dcoordsi(stripy,viewy) + zi = 0.D0 c------------------------------------------------------------------------ c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c------------------------------------------------------------------------ @@ -3250,6 +3076,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) @@ -3351,8 +3179,8 @@ 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, @@ -3498,6 +3326,7 @@ *----- single clusters ----------------------------------------------- c print*,'## ncls(',ip,') ',ncls(ip) do ic=1,ncls(ip) !loop on single clusters +c print*,'-',ic,'-' 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) @@ -3586,9 +3415,13 @@ * ---------------------------- xm_A(nplanes-ip+1) = xmm_A ym_A(nplanes-ip+1) = ymm_A + zm_A(nplanes-ip+1) = zmm_A xm_B(nplanes-ip+1) = xmm_B ym_B(nplanes-ip+1) = ymm_B + zm_B(nplanes-ip+1) = zmm_B zm(nplanes-ip+1) = (zmm_A+zmm_B)/2. +c$$$ zm(nplanes-ip+1) = +c$$$ $ z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< * ---------------------------- @@ -3732,6 +3565,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 +3590,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 @@ -3924,12 +3767,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 @@ -3950,10 +3797,6 @@ cl_used(cltrx(ip,ntr)) = 1 !tag used clusters cl_used(cltry(ip,ntr)) = 1 !tag used clusters -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))) @@ -3963,23 +3806,51 @@ 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 + 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 @@ -4036,12 +3907,21 @@ ip=nplanes-npl(VIEW(icl))+1 if(cl_used(icl).eq.0)then !cluster not included in any track + +cc print*,' ic ',icl,' --- stored ' + 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) + +cc print*,icl,' >>>> ',sxbad(nclsx) + do is=1,2 c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.) c call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.) @@ -4059,6 +3939,11 @@ 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) + +cc print*,icl,' >>>> ',sybad(nclsy) + do is=1,2 c call xyz_PAM(0,icl,is,' ','COG1',0.,0.) c call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.) @@ -4082,24 +3967,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