--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2008/12/23 11:28:36 1.38 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2014/01/16 15:29:50 1.45 @@ -271,7 +271,8 @@ *------------------------------------------------------------------------------- ccc ntrk=0 !counter of identified physical tracks -11111 continue !<<<<<<< come here when performing a new search +c11111 continue !<<<<<<< come here when performing a new search + continue !<<<<<<< come here when performing a new search if(nclouds_xz.eq.0)goto 880 !go to next event if(nclouds_yz.eq.0)goto 880 !go to next event @@ -302,7 +303,7 @@ * 1st) decreasing n.points * 2nd) increasing chi**2 * ------------------------------------------------------- - rchi2best=1000000000. + rchi2best=1000000000. ndofbest=0 do i=1,ntracks ndof=0 @@ -527,14 +528,6 @@ * --- and store the results -------------------------------- - ntrk = ntrk + 1 !counter of found tracks - if(.not.FIMAGE - $ .and.iimage.eq.0) image(ntrk)= 0 - if(.not.FIMAGE - $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next - if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous - call fill_level2_tracks(ntrk) !==> good2=.true. - if(ntrk.eq.NTRKMAX)then if(verbose.eq.1) $ print*, @@ -544,6 +537,24 @@ cc good2=.false. goto 880 !fill ntp and go to next event endif + + ntrk = ntrk + 1 !counter of found tracks + if(.not.FIMAGE + $ .and.iimage.eq.0) image(ntrk)= 0 + if(.not.FIMAGE + $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next + if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous + call fill_level2_tracks(ntrk) !==> good2=.true. + +c$$$ if(ntrk.eq.NTRKMAX)then +c$$$ if(verbose.eq.1) +c$$$ $ print*, +c$$$ $ '** warning ** number of identified '// +c$$$ $ 'tracks exceeds vector dimension ' +c$$$ $ ,'( ',NTRKMAX,' )' +c$$$cc good2=.false. +c$$$ goto 880 !fill ntp and go to next event +c$$$ endif if(iimage.ne.0)then FIMAGE=.true. ! goto 1212 !>>> fit image-track @@ -705,7 +716,8 @@ stripx = stripx + corr resxPAM = res - 10 endif + 10 continue + endif * ----------------- * CLUSTER Y @@ -751,7 +763,8 @@ stripy = stripy + corr resyPAM = res - 20 endif + 20 continue + endif c=========================================================== @@ -1366,8 +1379,8 @@ iv1=iside iv2=mod(iside,4)+1 * straight line passing trhough two consecutive vertexes - AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2)) - BB = yvv(iv1) - AA*xvv(iv1) + AA = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7 + BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7 * point along the straight line closer to the track xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2) yoo = AA*xoo + BB @@ -1379,8 +1392,8 @@ iv1=iside iv2=mod(iside,4)+1 * straight line passing trhough two consecutive vertexes - AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2)) - BB = xvv(iv1) - AA*yvv(iv1) + AA = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7 + BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7 * point along the straight line closer to the track yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2) xoo = AA*yoo + BB @@ -1581,7 +1594,8 @@ integer iflag integer badseed,badclx,badcly - + + iflag = iflag if(DEBUG.EQ.1)print*,'cl_to_couples:' cc if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80 @@ -1767,8 +1781,8 @@ $ ,'( ',ncouplemax,' ) --> masked!' 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 + mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1 + mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1 goto 10 endif @@ -1796,7 +1810,8 @@ endif enddo - 80 continue +c 80 continue + continue if(DEBUG.EQ.1)then @@ -1961,9 +1976,9 @@ c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',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 + xm1=REAL(xPAM) !EM GCC4.7 + ym1=REAL(yPAM) !EM GCC4.7 + zm1=REAL(zPAM) !EM GCC4.7 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2 c$$$ print*,'(2) ip ',ip2 @@ -1984,9 +1999,9 @@ c $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.) - xm2=xPAM - ym2=yPAM - zm2=zPAM + xm2=REAL(xPAM) !EM GCC4.7 + ym2=REAL(yPAM) !EM GCC4.7 + zm2=REAL(zPAM) !EM GCC4.7 * --------------------------------------------------- * both couples must have a y-cluster @@ -2027,7 +2042,7 @@ * tg(th_yz) alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2) * y0 (cm) - alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1 + alfayz1(ndblt)=alfayz2(ndblt)*(REAL(zini)-zm1)+ym1 ! EM GCC4.7 zm1, ym1 and alfayz1/2 are REAL **** -----------------------------------------------**** **** reject non phisical couples **** @@ -2090,9 +2105,9 @@ call xyz_PAM $ (icx3,icy3,is3,PFAdef,PFAdef $ ,0.,0.,0.,0.) - xm3=xPAM - ym3=yPAM - zm3=zPAM + xm3=REAL(xPAM) !EM GCC4.7 + ym3=REAL(yPAM) !EM GCC4.7 + zm3=REAL(zPAM) !EM GCC4.7 * find the circle passing through the three points @@ -2124,12 +2139,12 @@ SZX=SZX+ZP(I)*XX SSX=SSX+XX SZ=SZ+ZP(I) - S1=S1+1. + S1=S1+1. ENDDO DET=SZZ*S1-SZ*SZ AX=(SZX*S1-SZ*SSX)/DET BX=(SZZ*SSX-SZX*SZ)/DET - X0 = AX*ZINI+BX + X0 = AX*REAL(ZINI)+BX ! EM GCC4.7 endif @@ -2168,14 +2183,16 @@ if(radius.ne.0.and.xc.lt.0)then *************POSITIVE DEFLECTION - alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2) - alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) - alfaxz3(ntrpt) = 1/radius + alfaxz1(ntrpt) = xc+sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 + alfaxz2(ntrpt) = (REAL(ZINI)-zc)/ + $ sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 + alfaxz3(ntrpt) = 1/radius else if(radius.ne.0.and.xc.ge.0)then *************NEGATIVE DEFLECTION - alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2) - alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) - alfaxz3(ntrpt) = -1/radius + alfaxz1(ntrpt) = xc-sqrt(radius**2-(REAL(ZINI)-zc)**2) + alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/ + $ sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 + alfaxz3(ntrpt) = -1/radius else if(radius.eq.0)then *************straight fit alfaxz1(ntrpt) = X0 @@ -2210,12 +2227,14 @@ enddo !end loop on planes - COPPIA 3 31 continue - 1 enddo !end loop on COPPIA 2 +c 1 enddo !end loop on COPPIA 2 + enddo !end loop on COPPIA 2 enddo !end loop on sensors - COPPIA 2 20 continue enddo !end loop on planes - COPPIA 2 - 11 continue +c 11 continue + continue enddo !end loop on COPPIA1 enddo !end loop on sensors - COPPIA 1 10 continue @@ -2322,7 +2341,7 @@ * doublet distance in parameter space distance= $ ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2 - $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2 + $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2 distance = sqrt(distance) if(distance.lt.cutdistyz)then @@ -2347,7 +2366,8 @@ 1118 continue enddo !end loop (2) on DOUBLETS - 1188 continue +c 1188 continue + continue enddo !end loop on... bo? nptloop=npv @@ -2535,7 +2555,7 @@ * solo i due parametri spaziali per il momemnto distance= $ ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2 - $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2 + $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2 distance = sqrt(distance) @@ -2574,7 +2594,8 @@ 11188 continue enddo !end loop (2) on TRIPLETS -11888 continue +c11888 continue + continue enddo !end loop on... bo? nptloop=npv @@ -2639,7 +2660,7 @@ npt_tot=npt_tot+npt if(DEBUG.EQ.1)then - print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' + print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz) print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz) print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz) @@ -2942,7 +2963,7 @@ xm(nplanes-ip+1)=xPAM ym(nplanes-ip+1)=yPAM zm(nplanes-ip+1)=zPAM - resx(nplanes-ip+1)=resxPAM + resx(nplanes-ip+1)=resxPAM resy(nplanes-ip+1)=resyPAM if(DEBUG.EQ.1)print*,'(X,Y)' $ ,nplanes-ip+1,xPAM,yPAM @@ -2953,7 +2974,7 @@ ym_B(nplanes-ip+1) = yPAM_B zm(nplanes-ip+1) $ = (zPAM_A+zPAM_B)/2. - resx(nplanes-ip+1) = resxPAM + resx(nplanes-ip+1) = resxPAM resy(nplanes-ip+1) = resyPAM if(icx.eq.0.and.icy.gt.0)then xgood(nplanes-ip+1)=0. @@ -3040,7 +3061,7 @@ XV_STORE(ip,ntracks)=sngl(xv(ip)) YV_STORE(ip,ntracks)=sngl(yv(ip)) - ZV_STORE(ip,ntracks)=sngl(zv(ip)) + ZV_STORE(ip,ntracks)=sngl(zv(ip)) XM_STORE(ip,ntracks)=sngl(xm(ip)) YM_STORE(ip,ntracks)=sngl(ym(ip)) ZM_STORE(ip,ntracks)=sngl(zm(ip)) @@ -3085,7 +3106,7 @@ enddo enddo - RCHI2_STORE(ntracks)=chi2 + RCHI2_STORE(ntracks)=REAL(chi2) * -------------------------------- * STORE candidate TRACK INFO - end @@ -3153,6 +3174,12 @@ character*10 PFA common/FINALPFA/PFA + double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7 + double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7 + double precision zmm,zmm_A,zmm_B !EM GCC4.7 + double precision clincnewc !EM GCC4.7 + double precision clincnew !EM GCC4.7 + real k(6) DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/ @@ -3371,8 +3398,8 @@ idm = id dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) - clincnewc=10*sqrt(rymm**2+rxmm**2 - $ +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) + clincnewc=10.*dsqrt(rymm**2+rxmm**2 + $ +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))) ! EM GCC4.7 endif 1188 continue enddo !end loop on couples on plane icp @@ -3540,12 +3567,14 @@ if(iclm.ne.0)then if(mod(VIEW(iclm),2).eq.0)then clincnew= - $ 20* - $ sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1)) + $ 20.* !EM GCC4.7 + $ dsqrt(rxmm**2 + + $ DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1)) else if(mod(VIEW(iclm),2).ne.0)then clincnew= - $ 10* - $ sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2)) + $ 10.* !EM GCC4.7 + $ dsqrt(rymm**2 + + $ DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2)) endif if(distmin.le.clincnew)then @@ -3794,7 +3823,7 @@ character*10 PFA common/FINALPFA/PFA - real sinth,phi,pig + real sinth,phi,pig, npig ! EM GCC4.7 integer ssensor,sladder pig=acos(-1.) @@ -3804,8 +3833,8 @@ chi2_nt(ntr) = sngl(chi2) nstep_nt(ntr) = nstep * ------------------------------------- - phi = al(4) - sinth = al(3) + phi = REAL(al(4)) + sinth = REAL(al(3)) if(sinth.lt.0)then sinth = -sinth phi = phi + pig @@ -3875,7 +3904,21 @@ if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align LS(IP,ntr) = ssensor+10*sladder - if(id.ne.0)then +c if(id.ne.0)then +CCCCCC(10 novembre 2009) PATCH X NUCLEI +C non ho capito perche', ma durante il ritracciamento dei nuclei +C (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile +C che non viene reinizializzata correttamente e i cluster esclusi +C dal fit risultano ancora inclusi... +C + cltrx(ip,ntr) = 0 + cltry(ip,ntr) = 0 +c$$$ if( +c$$$ $ xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1 +c$$$ $ .and. +c$$$ $ id.ne.0)then + if(id.ne.0)then !patch 30/12/09 + c >>> is a couple cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id)) cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id)) @@ -3884,7 +3927,7 @@ cl_used(cltrx(ip,ntr)) = 1 !tag used clusters - xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id))) + xbad(ip,ntr)= nbadstrips(4,clx(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) @@ -3914,7 +3957,9 @@ ypu(ip,ntr) = corr endif - elseif(icl.ne.0)then +c$$$ elseif(icl.ne.0)then + endif !patch 30/12/09 + if(icl.ne.0)then !patch 30/12/09 cl_used(icl) = 1 !tag used clusters @@ -4012,8 +4057,8 @@ do is=1,2 c call xyz_PAM(icl,0,is,'COG1',' ',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 + call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.) + xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7 enddo else !=== Y views nclsy = nclsy + 1 @@ -4028,8 +4073,8 @@ do is=1,2 c call xyz_PAM(0,icl,is,' ','COG1',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 + call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.) + ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7 enddo endif endif @@ -4105,7 +4150,7 @@ alfayz2_av_nt(iyz)=alfayz2_av(iyz) nnn=nnn+ptcloud_yz(iyz) enddo - do ipt=1,nnn + do ipt=1,min(ndblt_max_nt,nnn) db_cloud_nt(ipt)=db_cloud(ipt) enddo endif @@ -4118,7 +4163,7 @@ alfaxz3_av_nt(ixz)=alfaxz3_av(ixz) nnn=nnn+ptcloud_xz(ixz) enddo - do ipt=1,nnn + do ipt=1,min(ntrpt_max_nt,nnn) tr_cloud_nt(ipt)=tr_cloud(ipt) enddo endif