--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/11/21 14:00:40 1.14 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/02/16 14:56:02 1.18 @@ -233,24 +233,53 @@ c$$$ enddo c$$$ if(ibest.eq.0)goto 880 !>> no good candidates - rchi2best=1000000000. +* ------------------------------------------------------- +* order track-candidates according to: +* 1st) decreasing n.points +* 2nd) increasing chi**2 +* ------------------------------------------------------- + rchi2best=1000000000. ndofbest=0 !(1) do i=1,ntracks - if(RCHI2_STORE(i).lt.rchi2best.and. - $ RCHI2_STORE(i).gt.0)then - ndof=0 !(1) - do ii=1,nplanes !(1) - ndof=ndof !(1) - $ +int(xgood_store(ii,i)) !(1) - $ +int(ygood_store(ii,i)) !(1) - enddo !(1) - if(ndof.ge.ndofbest)then !(1) + ndof=0 !(1) + do ii=1,nplanes !(1) + ndof=ndof !(1) + $ +int(xgood_store(ii,i)) !(1) + $ +int(ygood_store(ii,i)) !(1) + enddo !(1) + if(ndof.gt.ndofbest)then !(1) + ibest=i + rchi2best=RCHI2_STORE(i) + ndofbest=ndof !(1) + elseif(ndof.eq.ndofbest)then !(1) + if(RCHI2_STORE(i).lt.rchi2best.and. + $ RCHI2_STORE(i).gt.0)then ibest=i rchi2best=RCHI2_STORE(i) ndofbest=ndof !(1) endif !(1) 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 *------------------------------------------------------------------------------- * The best track candidate (ibest) is selected and a new fitting is performed. @@ -285,6 +314,7 @@ 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)) @@ -297,9 +327,10 @@ iprint=0 c if(DEBUG)iprint=1 if(VERBOSE)iprint=1 + if(DEBUG)iprint=2 call mini2(jstep,ifail,iprint) if(ifail.ne.0) then - if(.true.)then + if(VERBOSE)then print *, $ '*** MINIMIZATION FAILURE *** (after refinement) ' $ ,iev @@ -446,6 +477,8 @@ * PFAy - Position Finding Algorithm in y (COG2,ETA2,...) * angx - Projected angle in x * angy - Projected angle in y +* bfx - x-component of magnetci field +* bfy - y-component of magnetic field * * --------- COUPLES ------------------------------------------------------- * The couple defines a point in the space. @@ -484,16 +517,8 @@ * * - subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy) + subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) -c***************************************************** -c 07/10/2005 modified by elena vannuccini --> (1) -c 01/02/2006 modified by elena vannuccini --> (2) -c 02/02/2006 modified by Elena Vannuccini --> (3) -c (implemented new p.f.a.) -c 03/02/2006 modified by Elena Vannuccini --> (4) -c (implemented variable resolution) -c***************************************************** include 'commontracker.f' include 'level1.f' @@ -511,7 +536,9 @@ integer sensor integer viewx,viewy character*4 PFAx,PFAy !PFA to be used - real angx,angy !X-Y angle + real ax,ay !X-Y geometric angle + real angx,angy !X-Y effective angle + real bfx,bfy !X-Y b-field components real stripx,stripy @@ -537,94 +564,126 @@ xPAM_B = 0. yPAM_B = 0. zPAM_B = 0. - +c print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy * ----------------- * CLUSTER X * ----------------- if(icx.ne.0)then + viewx = VIEW(icx) nldx = nld(MAXS(icx),VIEW(icx)) nplx = npl(VIEW(icx)) - resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! - + resxPAM = RESXAV stripx = float(MAXS(icx)) - if(PFAx.eq.'COG1')then !(1) - stripx = stripx !(1) - resxPAM = resxPAM !(1) +* -------------------------- +* magnetic-field corrections +* -------------------------- +c$$$ print*,nplx,ax,bfy/10. + angtemp = ax + bfytemp = bfy + if(nplx.eq.6) angtemp = -1. * ax + if(nplx.eq.6) 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*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX +c$$$ print*,'========================' +* -------------------------- + + if(PFAx.eq.'COG1')then + stripx = stripx + resxPAM = resxPAM elseif(PFAx.eq.'COG2')then stripx = stripx + cog(2,icx) resxPAM = resxPAM*fbad_cog(2,icx) + elseif(PFAx.eq.'COG3')then + stripx = stripx + cog(3,icx) + resxPAM = resxPAM*fbad_cog(3,icx) + elseif(PFAx.eq.'COG4')then +c print*,'COG4' + stripx = stripx + cog(4,icx) + resxPAM = resxPAM*fbad_cog(4,icx) elseif(PFAx.eq.'ETA2')then -c cog2 = cog(2,icx) -c etacorr = pfaeta2(cog2,viewx,nldx,angx) -c stripx = stripx + etacorr - stripx = stripx + pfaeta2(icx,angx) !(3) - resxPAM = risx_eta2(angx) ! (4) + stripx = stripx + pfaeta2(icx,angx) + resxPAM = risx_eta2(abs(angx)) if(DEBUG.and.fbad_cog(2,icx).ne.1) $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) resxPAM = resxPAM*fbad_cog(2,icx) - elseif(PFAx.eq.'ETA3')then !(3) - stripx = stripx + pfaeta3(icx,angx) !(3) - resxPAM = risx_eta3(angx) ! (4) - if(DEBUG.and.fbad_cog(3,icx).ne.1) !(3) - $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3) - resxPAM = resxPAM*fbad_cog(3,icx) !(3) - elseif(PFAx.eq.'ETA4')then !(3) - stripx = stripx + pfaeta4(icx,angx) !(3) - resxPAM = risx_eta4(angx) ! (4) - if(DEBUG.and.fbad_cog(4,icx).ne.1) !(3) - $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3) - resxPAM = resxPAM*fbad_cog(4,icx) !(3) - elseif(PFAx.eq.'ETA')then !(3) - stripx = stripx + pfaeta(icx,angx) !(3) - resxPAM = ris_eta(icx,angx) ! (4) - if(DEBUG.and.fbad_cog(2,icx).ne.1) !(3) - $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3) -c resxPAM = resxPAM*fbad_cog(2,icx) !(3)TEMPORANEO - resxPAM = resxPAM*fbad_eta(icx,angx) !(3)(4) - elseif(PFAx.eq.'COG')then !(2) - stripx = stripx + cog(0,icx) !(2) - resxPAM = risx_cog(angx) ! (4) - resxPAM = resxPAM*fbad_cog(0,icx)!(2) + elseif(PFAx.eq.'ETA3')then + stripx = stripx + pfaeta3(icx,angx) + resxPAM = risx_eta3(abs(angx)) + if(DEBUG.and.fbad_cog(3,icx).ne.1) + $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx) + resxPAM = resxPAM*fbad_cog(3,icx) + elseif(PFAx.eq.'ETA4')then + stripx = stripx + pfaeta4(icx,angx) + resxPAM = risx_eta4(abs(angx)) + if(DEBUG.and.fbad_cog(4,icx).ne.1) + $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx) + resxPAM = resxPAM*fbad_cog(4,icx) + elseif(PFAx.eq.'ETA')then +c print*,'ETA' + stripx = stripx + pfaeta(icx,angx) + resxPAM = ris_eta(icx,angx) + if(DEBUG.and.fbad_cog(2,icx).ne.1) + $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) + resxPAM = resxPAM*fbad_eta(icx,angx) + elseif(PFAx.eq.'COG')then + stripx = stripx + cog(0,icx) + resxPAM = risx_cog(abs(angx)) + resxPAM = resxPAM*fbad_cog(0,icx) else print*,'*** Non valid p.f.a. (x) --> ',PFAx endif +c print*,'%%%%%%%%%%%%' + endif -c if(icy.eq.0.and.icx.ne.0) -c $ print*,PFAx,icx,angx,stripx,resxPAM,'***' * ----------------- * CLUSTER Y * ----------------- if(icy.ne.0)then + viewy = VIEW(icy) nldy = nld(MAXS(icy),VIEW(icy)) nply = npl(VIEW(icy)) - resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! - + resyPAM = RESYAV + stripy = float(MAXS(icy)) if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' $ ,icx,icy 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 +* -------------------------- - stripy = float(MAXS(icy)) if(PFAy.eq.'COG1')then !(1) stripy = stripy !(1) resyPAM = resyPAM !(1) elseif(PFAy.eq.'COG2')then stripy = stripy + cog(2,icy) resyPAM = resyPAM*fbad_cog(2,icy) + elseif(PFAy.eq.'COG3')then + stripy = stripy + cog(3,icy) + resyPAM = resyPAM*fbad_cog(3,icy) + elseif(PFAy.eq.'COG4')then + stripy = stripy + cog(4,icy) + resyPAM = resyPAM*fbad_cog(4,icy) elseif(PFAy.eq.'ETA2')then c cog2 = cog(2,icy) c etacorr = pfaeta2(cog2,viewy,nldy,angy) c stripy = stripy + etacorr stripy = stripy + pfaeta2(icy,angy) !(3) - resyPAM = risy_eta2(angy) ! (4) + resyPAM = risy_eta2(abs(angy)) ! (4) resyPAM = resyPAM*fbad_cog(2,icy) if(DEBUG.and.fbad_cog(2,icy).ne.1) $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy) @@ -647,7 +706,7 @@ $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3) elseif(PFAy.eq.'COG')then stripy = stripy + cog(0,icy) - resyPAM = risy_cog(angy) ! (4) + resyPAM = risy_cog(abs(angy)) ! (4) c resyPAM = ris_eta(icy,angy) ! (4) resyPAM = resyPAM*fbad_cog(0,icy) else @@ -656,7 +715,8 @@ endif - +c print*,'## stripx,stripy ',stripx,stripy + c=========================================================== C COUPLE C=========================================================== @@ -858,6 +918,11 @@ 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 @@ -1348,7 +1413,7 @@ * ---------------------------------------------------- * cut on charge (X VIEW) * ---------------------------------------------------- - if(dedx(icx).lt.dedx_x_min)then + if(sgnl(icx).lt.dedx_x_min)then cl_single(icx)=0 goto 10 endif @@ -1398,7 +1463,7 @@ * ---------------------------------------------------- * cut on charge (Y VIEW) * ---------------------------------------------------- - if(dedx(icy).lt.dedx_y_min)then + if(sgnl(icy).lt.dedx_y_min)then cl_single(icy)=0 goto 20 endif @@ -1444,21 +1509,21 @@ * charge correlation * (modified to be applied only below saturation... obviously) - if( .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx) + if( .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx) $ .and. - $ .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx) + $ .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx) $ .and. $ (badclx.eq.1.and.badcly.eq.1) $ .and. $ .true.)then - ddd=(dedx(icy) - $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) + ddd=(sgnl(icy) + $ -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx)) ddd=ddd/sqrt(kch(nplx,nldx)**2+1) c cut = chcut * sch(nplx,nldx) - sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx) + sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx) $ -kch(nplx,nldx)*cch(nplx,nldx)) sss=sss/sqrt(kch(nplx,nldx)**2+1) cut = chcut * (16 + sss/50.) @@ -1467,15 +1532,8 @@ goto 20 !charge not consistent endif endif - -* ------------------> COUPLE <------------------ - ncp_plane(nplx) = ncp_plane(nplx) + 1 - clx(nplx,ncp_plane(nplx))=icx - cly(nply,ncp_plane(nplx))=icy - cl_single(icx)=0 - cl_single(icy)=0 - if(ncp_plane(nplx).eq.ncouplemax)then + if(ncp_plane(nplx).gt.ncouplemax)then if(verbose)print*, $ '** warning ** number of identified '// $ 'couples on plane ',nplx, @@ -1483,7 +1541,15 @@ $ ,'( ',ncouplemax,' ) --> masked!' mask_view(nviewx(nplx)) = 2 mask_view(nviewy(nply)) = 2 + goto 10 endif + +* ------------------> COUPLE <------------------ + ncp_plane(nplx) = ncp_plane(nplx) + 1 + clx(nplx,ncp_plane(nplx))=icx + cly(nply,ncp_plane(nplx))=icy + cl_single(icx)=0 + cl_single(icy)=0 * ---------------------------------------------- endif @@ -1592,7 +1658,8 @@ icx1=clx(ip1,icp1) icy1=cly(ip1,icp1) c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1) - call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,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 @@ -1608,8 +1675,10 @@ icy2=cly(ip2,icp2) c call xyz_PAM c $ (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1) +c call xyz_PAM +c $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM - $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1) + $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.) xm2=xPAM ym2=yPAM zm2=zPAM @@ -1673,8 +1742,11 @@ icy3=cly(ip3,icp3) c call xyz_PAM c $ (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1) +c call xyz_PAM +c $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM - $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1) + $ (icx3,icy3,is3,PFAdef,PFAdef + $ ,0.,0.,0.,0.) xm3=xPAM ym3=yPAM zm3=zPAM @@ -2323,7 +2395,8 @@ nplused=nplused+ hit_plane(ip) enddo - if(nplused.lt.nplxz_min)goto 888 !next doublet +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 @@ -2427,8 +2500,10 @@ * ************************* c call xyz_PAM(icx,icy,is, c $ 'COG2','COG2',0.,0.) +c call xyz_PAM(icx,icy,is, !(1) +c $ PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM(icx,icy,is, !(1) - $ PFAdef,PFAdef,0.,0.) !(1) + $ PFAdef,PFAdef,0.,0.,0.,0.) * ************************* * ----------------------------- xgood(nplanes-ip+1)=1. @@ -2456,7 +2531,7 @@ jstep=0 !number of minimization steps iprint=0 c if(DEBUG)iprint=1 - if(DEBUG)iprint=1 + if(DEBUG)iprint=2 call mini2(jstep,ifail,iprint) if(ifail.ne.0) then if(DEBUG)then @@ -2600,11 +2675,14 @@ c include 'level1.f' include 'calib.f' - * flag to chose PFA character*10 PFA common/FINALPFA/PFA + real xp,yp,zp + real xyzp(3),bxyz(3) + equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3)) + * ================================================= * new estimate of positions using ETA algorithm * and @@ -2613,6 +2691,13 @@ call track_init do ip=1,nplanes !loop on planes + xP=XV_STORE(nplanes-ip+1,ibest) + yP=YV_STORE(nplanes-ip+1,ibest) + zP=ZV_STORE(nplanes-ip+1,ibest) + call gufld(xyzp,bxyz) +c$$$ bxyz(1)=0 +c$$$ bxyz(2)=0 +c$$$ bxyz(3)=0 * ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * If the plane has been already included, it just @@ -2633,11 +2718,17 @@ $ ,ip_cp(id),ip icx=clx(ip,icp) icy=cly(ip,icp) +c call xyz_PAM(icx,icy,is, +c $ PFA,PFA, +c $ AXV_STORE(nplanes-ip+1,ibest), +c $ AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(icx,icy,is, -c $ 'ETA2','ETA2', $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest), - $ AYV_STORE(nplanes-ip+1,ibest)) + $ AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) c$$$ call xyz_PAM(icx,icy,is, c$$$ $ 'COG2','COG2', c$$$ $ 0., @@ -2650,9 +2741,9 @@ resx(nplanes-ip+1) = resxPAM resy(nplanes-ip+1) = resyPAM -c dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1) - dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) - dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) +c dedxtrk(nplanes-ip+1) = (sgnl(icx)+sgnl(icy))/2. !(1) + dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) + dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) * ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- @@ -2667,9 +2758,9 @@ * -------------------------------------------------------------- * determine which ladder and sensor are intersected by the track - xP=XV_STORE(nplanes-ip+1,ibest) - yP=YV_STORE(nplanes-ip+1,ibest) - zP=ZV_STORE(nplanes-ip+1,ibest) +c$$$ xP=XV_STORE(nplanes-ip+1,ibest) +c$$$ yP=YV_STORE(nplanes-ip+1,ibest) +c$$$ zP=ZV_STORE(nplanes-ip+1,ibest) call whichsensor(ip,xP,yP,nldt,ist) * if the track hit the plane in a dead area, go to the next plane if(nldt.eq.0.or.ist.eq.0)goto 133 @@ -2707,11 +2798,17 @@ $ cl_used(icy).ne.0.or. !or the Y cluster is already used !(3) $ .false.)goto 1188 !then jump to next couple. * +c call xyz_PAM(icx,icy,ist, +c $ PFA,PFA, +c $ AXV_STORE(nplanes-ip+1,ibest), +c $ AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(icx,icy,ist, $ PFA,PFA, -c $ 'ETA2','ETA2', $ AXV_STORE(nplanes-ip+1,ibest), - $ AYV_STORE(nplanes-ip+1,ibest)) + $ AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) distance = distance_to(XP,YP) distance = distance / RCHI2_STORE(ibest)!<<< MS @@ -2726,9 +2823,9 @@ rymm = resyPAM distmin = distance idm = id -c dedxmm = (dedx(icx)+dedx(icy))/2. !(1) - dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) - dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) +c dedxmm = (sgnl(icx)+sgnl(icy))/2. !(1) + dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) + dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) endif 1188 continue enddo !end loop on couples on plane icp @@ -2780,10 +2877,15 @@ c if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used !(3) * !jump to the Y cluster +c call xyz_PAM(icx,0,ist, +c $ PFA,PFA, +c $ AXV_STORE(nplanes-ip+1,ibest),0.) call xyz_PAM(icx,0,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ AXV_STORE(nplanes-ip+1,ibest),0.) + $ AXV_STORE(nplanes-ip+1,ibest),0., + $ bxyz(1), + $ bxyz(2) + $ ) distance = distance_to(XP,YP) distance = distance / RCHI2_STORE(ibest)!<<< MS if(DEBUG)print*,'( cl-X ',icx @@ -2799,8 +2901,8 @@ rymm = resyPAM distmin = distance iclm = icx -c dedxmm = dedx(icx) !(1) - dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) +c dedxmm = sgnl(icx) !(1) + dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = 0. !(1) endif 11881 continue @@ -2808,10 +2910,15 @@ c if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3) * !jump to the next couple +c call xyz_PAM(0,icy,ist, +c $ PFA,PFA, +c $ 0.,AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(0,icy,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ 0.,AYV_STORE(nplanes-ip+1,ibest)) + $ 0.,AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) distance = distance_to(XP,YP) distance = distance / RCHI2_STORE(ibest)!<<< MS if(DEBUG)print*,'( cl-Y ',icy @@ -2827,36 +2934,49 @@ rymm = resyPAM distmin = distance iclm = icy -c dedxmm = dedx(icy) !(1) +c dedxmm = sgnl(icy) !(1) dedxmmx = 0. !(1) - dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) + dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) endif 11882 continue enddo !end loop on cluster inside couples -*----- single clusters ----------------------------------------------- +*----- single clusters ----------------------------------------------- +c print*,'## ncls(',ip,') ',ncls(ip) do ic=1,ncls(ip) !loop on single clusters icl=cls(ip,ic) +c print*,'## ic ',ic,' ist ',ist 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 if(mod(VIEW(icl),2).eq.0)then!<---- X view +c call xyz_PAM(icl,0,ist, +c $ PFA,PFA, +c $ AXV_STORE(nplanes-ip+1,ibest),0.) call xyz_PAM(icl,0,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ AXV_STORE(nplanes-ip+1,ibest),0.) + $ AXV_STORE(nplanes-ip+1,ibest),0., + $ bxyz(1), + $ bxyz(2) + $ ) else !<---- Y view +c call xyz_PAM(0,icl,ist, +c $ PFA,PFA, +c $ 0.,AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(0,icl,ist, -c $ 'ETA2','ETA2', $ PFA,PFA, - $ 0.,AYV_STORE(nplanes-ip+1,ibest)) + $ 0.,AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) endif distance = distance_to(XP,YP) distance = distance / RCHI2_STORE(ibest)!<<< MS if(DEBUG)print*,'( cl-s ',icl - $ ,' ) normalized distance ',distance + $ ,' ) normalized distance ',distance,'<',distmin,' ?' if(distance.lt.distmin)then + if(DEBUG)print*,'YES' xmm_A = xPAM_A ymm_A = yPAM_A zmm_A = zPAM_A @@ -2867,18 +2987,18 @@ rymm = resyPAM distmin = distance iclm = icl -c dedxmm = dedx(icl) !(1) +c dedxmm = sgnl(icl) !(1) if(mod(VIEW(icl),2).eq.0)then !<---- X view - dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) + dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) dedxmmy = 0. !(1) else !<---- Y view dedxmmx = 0. !(1) - dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) + dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) endif endif 18882 continue enddo !end loop on single clusters - +c print*,'## distmin ', distmin,' clinc ',clinc if(distmin.le.clinc)then CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< @@ -3218,9 +3338,15 @@ yv_nt(ip,ntr) = sngl(yv(ip)) zv_nt(ip,ntr) = sngl(zv(ip)) axv_nt(ip,ntr) = sngl(axv(ip)) - ayv_nt(ip,ntr) = sngl(ayv(ip)) - dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)) !(2) - dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)) !(2) + ayv_nt(ip,ntr) = sngl(ayv(ip)) +c l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + factor = sqrt( + $ sin( acos(-1.) * sngl(axv(ip)) /180. )**2 + + $ sin( 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) id = CP_STORE(ip,IDCAND) icl = CLS_STORE(ip,IDCAND) @@ -3270,11 +3396,12 @@ if(mod(VIEW(icl),2).eq.0)then !=== X views nclsx = nclsx + 1 planex(nclsx) = ip - sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2) + sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2) clsx(nclsx) = icl do is=1,2 c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.) - call xyz_PAM(icl,0,is,PFAdef,' ',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 @@ -3285,11 +3412,12 @@ else !=== Y views nclsy = nclsy + 1 planey(nclsy) = ip - sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2) + sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2) clsy(nclsy) = icl do is=1,2 c call xyz_PAM(0,icl,is,' ','COG1',0.,0.) - call xyz_PAM(0,icl,is,' ',PFAdef,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