--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/08/31 14:56:52 1.31 +++ 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 ' @@ -802,7 +805,7 @@ 20 endif -c$$$ print*,'## stripx,stripy ',stripx,stripy +cc print*,'## stripx,stripy ',stripx,stripy c=========================================================== C COUPLE @@ -819,9 +822,9 @@ $ ' 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 @@ -856,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. @@ -882,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 @@ -913,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 @@ -986,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 @@ -1050,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 @@ -1097,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 @@ -1124,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. @@ -1155,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 @@ -1220,7 +1247,7 @@ * ******************************************************************************** - real function distance_to(XPP,YPP) + real function distance_to(rXPP,rYPP) include 'common_xyzPAM.f' @@ -1229,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. @@ -1358,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) @@ -1383,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------------------------------------------------------------------------ @@ -3044,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) @@ -3145,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, @@ -3292,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) @@ -3380,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 !<<< * ---------------------------- @@ -3551,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 @@ -3779,11 +3822,11 @@ elseif(icl.ne.0)then + cl_used(icl) = 1 !tag used clusters if(mod(VIEW(icl),2).eq.0)then cltrx(ip,ntr)=icl - xbad(ip,ntr) = nbadstrips(4,icl) if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr) @@ -3797,7 +3840,6 @@ elseif(mod(VIEW(icl),2).eq.1)then cltry(ip,ntr)=icl - ybad(ip,ntr) = nbadstrips(4,icl) if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr) @@ -3865,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.) @@ -3888,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.) @@ -3911,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