--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/11/21 17:13:31 1.15 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/08/17 14:36:05 1.27 @@ -20,7 +20,40 @@ include 'calib.f' include 'level2.f' -c include 'momanhough_init.f' + + +c print*,'======================================================' +c$$$ do ic=1,NCLSTR1 +c$$$ if(.false. +c$$$ $ .or.nsatstrips(ic).gt.0 +c$$$c $ .or.nbadstrips(0,ic).gt.0 +c$$$c $ .or.nbadstrips(4,ic).gt.0 +c$$$c $ .or.nbadstrips(3,ic).gt.0 +c$$$ $ .or..false.)then +c$$$ print*,'--- cl-',ic,' ------------------------' +c$$$ istart = INDSTART(IC) +c$$$ istop = TOTCLLENGTH +c$$$ if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 +c$$$ print*,'ADC ',(CLADC(i),i=istart,istop) +c$$$ print*,'s/n ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) +c$$$ print*,'sgnl ',(CLSIGNAL(i),i=istart,istop) +c$$$ print*,'strip ',(i-INDMAX(ic),i=istart,istop) +c$$$ print*,'view ',VIEW(ic) +c$$$ print*,'maxs ',MAXS(ic) +c$$$ print*,'COG4 ',cog(4,ic) +c$$$ ff = fbad_cog(4,ic) +c$$$ print*,'fbad ',ff +c$$$ print*,(CLBAD(i),i=istart,istop) +c$$$ bb=nbadstrips(0,ic) +c$$$ print*,'#BAD (tot)',bb +c$$$ bb=nbadstrips(4,ic) +c$$$ print*,'#BAD (4)',bb +c$$$ bb=nbadstrips(3,ic) +c$$$ print*,'#BAD (3)',bb +c$$$ ss=nsatstrips(ic) +c$$$ print*,'#saturated ',ss +c$$$ endif +c$$$ enddo *------------------------------------------------------------------------------- * STEP 1 @@ -41,7 +74,7 @@ *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- -c iflag=0 + call cl_to_couples(iflag) if(iflag.eq.1)then !bad event goto 880 !go to next event @@ -75,7 +108,7 @@ *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- -c iflag=0 + call cp_to_doubtrip(iflag) if(iflag.eq.1)then !bad event goto 880 !go to next event @@ -239,25 +272,25 @@ * 2nd) increasing chi**2 * ------------------------------------------------------- rchi2best=1000000000. - ndofbest=0 !(1) + ndofbest=0 do i=1,ntracks - 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) + ndof=0 + do ii=1,nplanes + ndof=ndof + $ +int(xgood_store(ii,i)) + $ +int(ygood_store(ii,i)) + enddo + if(ndof.gt.ndofbest)then ibest=i rchi2best=RCHI2_STORE(i) - ndofbest=ndof !(1) - elseif(ndof.eq.ndofbest)then !(1) + ndofbest=ndof + elseif(ndof.eq.ndofbest)then if(RCHI2_STORE(i).lt.rchi2best.and. $ RCHI2_STORE(i).gt.0)then ibest=i rchi2best=RCHI2_STORE(i) - ndofbest=ndof !(1) - endif !(1) + ndofbest=ndof + endif endif enddo @@ -298,8 +331,10 @@ iimage=0 endif if(icand.eq.0)then - print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand - $ ,ibest,iimage + if(VERBOSE)then + print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand + $ ,ibest,iimage + endif return endif @@ -314,6 +349,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)) @@ -326,9 +362,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 @@ -475,6 +512,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. @@ -513,46 +552,35 @@ * * - 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' include 'calib.f' -c include 'level1.f' include 'common_align.f' include 'common_mech.f' include 'common_xyzPAM.f' -c include 'common_resxy.f' - -c logical DEBUG -c common/dbg/DEBUG integer icx,icy !X-Y cluster ID 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 double precision xrt,yrt,zrt double precision xrt_A,yrt_A,zrt_A double precision xrt_B,yrt_B,zrt_B -c double precision xi,yi,zi -c double precision xi_A,yi_A,zi_A -c double precision xi_B,yi_B,zi_B parameter (ndivx=30) + + +c$$$ print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy resxPAM = 0 resyPAM = 0 @@ -566,126 +594,305 @@ xPAM_B = 0. yPAM_B = 0. 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) - nldx = nld(MAXS(icx),VIEW(icx)) - nplx = npl(VIEW(icx)) - resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! - - stripx = float(MAXS(icx)) - if(PFAx.eq.'COG1')then !(1) - stripx = stripx !(1) - resxPAM = resxPAM !(1) + + viewx = VIEW(icx) + nldx = nld(MAXS(icx),VIEW(icx)) + 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 +* -------------------------- + 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) + + + if(PFAx.eq.'COG1')then + + stripx = stripx + resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM + elseif(PFAx.eq.'COG2')then - stripx = stripx + cog(2,icx) + + 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 -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 = 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) - 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 = 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 - print*,'*** Non valid p.f.a. (x) --> ',PFAx + if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx endif - endif -c if(icy.eq.0.and.icx.ne.0) -c $ print*,PFAx,icx,angx,stripx,resxPAM,'***' - + +* ====================================== +* temporary patch for saturated clusters +* ====================================== + if( nsatstrips(icx).gt.0 )then + stripx = stripx + cog(4,icx) + resxPAM = pitchX*1e-4/sqrt(12.) + cc=cog(4,icx) +c$$$ print*,icx,' *** ',cc +c$$$ print*,icx,' *** ',resxPAM + endif + + 10 endif + + * ----------------- * 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( + $ 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 - print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' - $ ,icx,icy + if(DEBUG) then + print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' + $ ,icx,icy + endif goto 100 endif - - stripy = float(MAXS(icy)) - if(PFAy.eq.'COG1')then !(1) - stripy = stripy !(1) - resyPAM = resyPAM !(1) +* -------------------------- +* 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 +* -------------------------- + +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) + + 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 -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) + + stripy = stripy + pfaeta2(icy,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) - elseif(PFAy.eq.'ETA3')then !(3) - stripy = stripy + pfaeta3(icy,angy) !(3) - resyPAM = resyPAM*fbad_cog(3,icy) !(3) - if(DEBUG.and.fbad_cog(3,icy).ne.1) !(3) - $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3) - elseif(PFAy.eq.'ETA4')then !(3) - stripy = stripy + pfaeta4(icy,angy) !(3) - resyPAM = resyPAM*fbad_cog(4,icy) !(3) - if(DEBUG.and.fbad_cog(4,icy).ne.1) !(3) - $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3) - elseif(PFAy.eq.'ETA')then !(3) - stripy = stripy + pfaeta(icy,angy) !(3) - resyPAM = ris_eta(icy,angy) ! (4) -c resyPAM = resyPAM*fbad_cog(2,icy) !(3)TEMPORANEO - resyPAM = resyPAM*fbad_eta(icy,angy) ! (4) - if(DEBUG.and.fbad_cog(2,icy).ne.1) !(3) - $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3) +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(angy) ! (4) -c resyPAM = ris_eta(icy,angy) ! (4) + + stripy = stripy + cog(0,icy) + resyPAM = risy_cog(abs(angy)) resyPAM = resyPAM*fbad_cog(0,icy) + else - print*,'*** Non valid p.f.a. (x) --> ',PFAx + if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx endif - endif - +* ====================================== +* temporary patch for saturated clusters +* ====================================== + if( nsatstrips(icy).gt.0 )then + stripy = stripy + cog(4,icy) + resyPAM = pitchY*1e-4/sqrt(12.) + cc=cog(4,icy) +c$$$ print*,icy,' *** ',cc +c$$$ print*,icy,' *** ',resyPAM + endif + + + 20 endif + +c$$$ print*,'## stripx,stripy ',stripx,stripy + c=========================================================== C COUPLE C=========================================================== @@ -696,8 +903,10 @@ 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... - print*,'xyz_PAM (couple):', - $ ' WARNING: false X strip: strip ',stripx + if(DEBUG) then + print*,'xyz_PAM (couple):', + $ ' WARNING: false X strip: strip ',stripx + endif endif xi = acoordsi(stripx,viewx) yi = acoordsi(stripy,viewy) @@ -789,8 +998,10 @@ c if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips... 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... - print*,'xyz_PAM (X-singlet):', - $ ' WARNING: false X strip: strip ',stripx + if(DEBUG) then + print*,'xyz_PAM (X-singlet):', + $ ' WARNING: false X strip: strip ',stripx + endif endif xi = acoordsi(stripx,viewx) @@ -812,10 +1023,11 @@ c print*,yi_A,' <--> ',yi_B else - - print *,'routine xyz_PAM ---> not properly used !!!' - print *,'icx = ',icx - print *,'icy = ',icy + if(DEBUG) then + print *,'routine xyz_PAM ---> not properly used !!!' + print *,'icx = ',icx + print *,'icy = ',icy + endif goto 100 endif @@ -880,16 +1092,200 @@ c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')' else - - print *,'routine xyz_PAM ---> not properly used !!!' - print *,'icx = ',icx - print *,'icy = ',icy - + if(DEBUG) then + print *,'routine xyz_PAM ---> not properly used !!!' + print *,'icx = ',icx + print *,'icy = ',icy + endif 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 +************************************************************************ +* 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 ******************************************************************************** ******************************************************************************** @@ -965,7 +1361,8 @@ endif distance= - $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2 + $ ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI +cc $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2 distance=dsqrt(distance) c$$$ print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b @@ -990,9 +1387,12 @@ * ---------------------- distance= - $ ((xPAM-XPP)/resxPAM)**2 - $ + - $ ((yPAM-YPP)/resyPAM)**2 + $ ((xPAM-XPP))**2 !QUIQUI + $ + + $ ((yPAM-YPP))**2 +c$$$ $ ((xPAM-XPP)/resxPAM)**2 +c$$$ $ + +c$$$ $ ((yPAM-YPP)/resyPAM)**2 distance=dsqrt(distance) c$$$ print*,xPAM,yPAM,zPAM @@ -1001,11 +1401,11 @@ else - print* - $ ,' function distance_to ---> wrong usage!!!' - print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM - print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b ' - $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b +c print* +c $ ,' function distance_to ---> wrong usage!!!' +c print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM +c print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b ' +c $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b endif distance_to = sngl(distance) @@ -1073,8 +1473,8 @@ 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((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips... - print*,'whichsensor: ', - $ ' WARNING: false X strip: strip ',stripx +c print*,'whichsensor: ', +c $ ' WARNING: false X strip: strip ',stripx endif xi = acoordsi(stripx,viewx) yi = acoordsi(stripy,viewy) @@ -1229,7 +1629,7 @@ is_cp=0 if(id.lt.0)is_cp=1 if(id.gt.0)is_cp=2 - if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!' +c if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!' return end @@ -1358,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 @@ -1377,7 +1778,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 @@ -1427,7 +1828,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 @@ -1473,21 +1874,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.) @@ -1503,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 @@ -1605,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 @@ -1622,7 +2027,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 @@ -1638,8 +2044,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 @@ -1656,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 @@ -1703,8 +2112,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 @@ -1732,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 @@ -1966,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 @@ -2176,7 +2590,7 @@ enddo c if(ncpused.lt.ncpxz_min)goto 22288 !next triplet if(npt.lt.nptxz_min)goto 22288 !next triplet - if(nplused.lt.nplxz_min)goto 22288 !next doublet + if(nplused.lt.nplxz_min)goto 22288 !next triplet * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< @@ -2188,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 @@ -2253,9 +2668,6 @@ ************************************************** subroutine clouds_to_ctrack(iflag) -c***************************************************** -c 02/02/2006 modified by Elena Vannuccini --> (1) -c***************************************************** include 'commontracker.f' include 'level1.f' @@ -2263,7 +2675,7 @@ include 'common_xyzPAM.f' include 'common_mini_2.f' include 'common_mech.f' -c include 'momanhough_init.f' + * output flag @@ -2458,8 +2870,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. @@ -2487,7 +2901,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 @@ -2522,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 @@ -2530,11 +2945,8 @@ ntracks = ntracks + 1 -c$$$ ndof=0 do ip=1,nplanes -c$$$ ndof=ndof -c$$$ $ +int(xgood(ip)) -c$$$ $ +int(ygood(ip)) + XV_STORE(ip,ntracks)=sngl(xv(ip)) YV_STORE(ip,ntracks)=sngl(yv(ip)) ZV_STORE(ip,ntracks)=sngl(zv(ip)) @@ -2553,19 +2965,26 @@ if(hit_plane(ip).ne.0)then CP_STORE(nplanes-ip+1,ntracks)= $ cp_match(ip,hit_plane(ip)) + SENSOR_STORE(nplanes-ip+1,ntracks) + $ = is_cp(cp_match(ip,hit_plane(ip))) + LADDER_STORE(nplanes-ip+1,ntracks) + $ = LADDER( + $ clx(ip,icp_cp( + $ cp_match(ip,hit_plane(ip) + $ )))); else CP_STORE(nplanes-ip+1,ntracks)=0 + SENSOR_STORE(nplanes-ip+1,ntracks)=0 + LADDER_STORE(nplanes-ip+1,ntracks)=0 endif + BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now + BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now CLS_STORE(nplanes-ip+1,ntracks)=0 do i=1,5 AL_STORE(i,ntracks)=sngl(AL(i)) enddo enddo -c$$$ * Number of Degree Of Freedom -c$$$ ndof=ndof-5 -c$$$ * reduced chi^2 -c$$$ rchi2=chi2/dble(ndof) RCHI2_STORE(ntracks)=chi2 * -------------------------------- @@ -2589,14 +3008,29 @@ return endif +c$$$ if(DEBUG)then +c$$$ print*,'****** TRACK CANDIDATES ***********' +c$$$ print*,'# R. chi2 RIG' +c$$$ do i=1,ntracks +c$$$ print*,i,' --- ',rchi2_store(i),' --- ' +c$$$ $ ,1./abs(AL_STORE(5,i)) +c$$$ enddo +c$$$ print*,'***********************************' +c$$$ endif if(DEBUG)then - print*,'****** TRACK CANDIDATES ***********' - print*,'# R. chi2 RIG' - do i=1,ntracks - print*,i,' --- ',rchi2_store(i),' --- ' - $ ,1./abs(AL_STORE(5,i)) - enddo - print*,'***********************************' + print*,'****** TRACK CANDIDATES *****************' + print*,'# R. chi2 RIG ndof' + do i=1,ntracks + 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) + print*,i,' --- ',rchi2_store(i),' --- ' + $ ,1./abs(AL_STORE(5,i)),' --- ',ndof + enddo + print*,'*****************************************' endif @@ -2615,11 +3049,6 @@ subroutine refine_track(ibest) -c****************************************************** -cccccc 06/10/2005 modified by elena vannuccini ---> (1) -cccccc 31/01/2006 modified by elena vannuccini ---> (2) -cccccc 12/08/2006 modified by elena vannucicni ---> (3) -c****************************************************** include 'commontracker.f' include 'level1.f' @@ -2627,15 +3056,19 @@ include 'common_xyzPAM.f' include 'common_mini_2.f' include 'common_mech.f' -c include 'momanhough_init.f' -c include 'level1.f' include 'calib.f' - * flag to chose PFA character*10 PFA common/FINALPFA/PFA + real k(6) + DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/ + + 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 @@ -2644,6 +3077,15 @@ 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) + BX_STORE(nplanes-ip+1,ibest)=bxyz(1) + BY_STORE(nplanes-ip+1,ibest)=bxyz(2) +c$$$ bxyz(1)=0 +c$$$ bxyz(2)=0 +c$$$ bxyz(3)=0 * ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * If the plane has been already included, it just @@ -2664,15 +3106,18 @@ $ ,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)) -c$$$ call xyz_PAM(icx,icy,is, -c$$$ $ 'COG2','COG2', -c$$$ $ 0., -c$$$ $ 0.) + $ AYV_STORE(nplanes-ip+1,ibest), + $ bxyz(1), + $ bxyz(2) + $ ) + xm(nplanes-ip+1) = xPAM ym(nplanes-ip+1) = yPAM zm(nplanes-ip+1) = zPAM @@ -2681,9 +3126,8 @@ 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) + dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) + dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) * ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- @@ -2698,12 +3142,12 @@ * -------------------------------------------------------------- * 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) 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 + + SENSOR_STORE(nplanes-ip+1,IBEST)=ist + LADDER_STORE(nplanes-ip+1,IBEST)=nldt * -------------------------------------------------------------- if(DEBUG)then @@ -2740,15 +3184,17 @@ * 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 +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI id=id_cp(ip,icp,ist) if(DEBUG)print*,'( couple ',id - $ ,' ) normalized distance ',distance + $ ,' ) distance ',distance if(distance.lt.distmin)then xmm = xPAM ymm = yPAM @@ -2757,28 +3203,30 @@ 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) + dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) + dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) +c QUIQUI --> non devo moltiplicare per clinc?!?!?! + clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI + $ +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI endif 1188 continue enddo !end loop on couples on plane icp - if(distmin.le.clinc)then +c if(distmin.le.clinc)then !QUIQUI + if(distmin.le.clincnewc)then !QUIQUI * ----------------------------------- - xm(nplanes-ip+1) = xmm !<<< - ym(nplanes-ip+1) = ymm !<<< - zm(nplanes-ip+1) = zmm !<<< - xgood(nplanes-ip+1) = 1 !<<< - ygood(nplanes-ip+1) = 1 !<<< - resx(nplanes-ip+1)=rxmm !<<< - resy(nplanes-ip+1)=rymm !<<< -c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1) - dedxtrk_x(nplanes-ip+1) = dedxmmx !(1) - dedxtrk_y(nplanes-ip+1) = dedxmmy !(1) + xm(nplanes-ip+1) = xmm !<<< + ym(nplanes-ip+1) = ymm !<<< + zm(nplanes-ip+1) = zmm !<<< + xgood(nplanes-ip+1) = 1 !<<< + ygood(nplanes-ip+1) = 1 !<<< + resx(nplanes-ip+1)=rxmm !<<< + resy(nplanes-ip+1)=rymm !<<< + dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< + dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< * ----------------------------------- CP_STORE(nplanes-ip+1,ibest)=idm if(DEBUG)print*,'%%%% included couple ',idm - $ ,' (norm.dist.= ',distmin,', cut ',clinc,' )' + $ ,' (dist.= ',distmin,', cut ',clinc,' )' goto 133 !next plane endif * ================================================ @@ -2811,14 +3259,19 @@ 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 +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI if(DEBUG)print*,'( cl-X ',icx - $ ,' in cp ',id,' ) normalized distance ',distance + $ ,' in cp ',id,' ) distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A @@ -2830,8 +3283,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 @@ -2839,14 +3292,19 @@ 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 +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI if(DEBUG)print*,'( cl-Y ',icy - $ ,' in cp ',id,' ) normalized distance ',distance + $ ,' in cp ',id,' ) distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A @@ -2858,13 +3316,14 @@ 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 if(cl_used(icl).eq.1.or. !if the cluster is already used @@ -2873,21 +3332,26 @@ $ .false.)goto 18882 !jump to the next singlet if(mod(VIEW(icl),2).eq.0)then!<---- X view 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 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 +c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI if(DEBUG)print*,'( cl-s ',icl - $ ,' ) normalized distance ',distance + $ ,' ) distance ',distance,'<',distmin,' ?' if(distance.lt.distmin)then + if(DEBUG)print*,'YES' xmm_A = xPAM_A ymm_A = yPAM_A zmm_A = zPAM_A @@ -2898,51 +3362,64 @@ rymm = resyPAM distmin = distance iclm = icl -c dedxmm = dedx(icl) !(1) if(mod(VIEW(icl),2).eq.0)then !<---- X view - dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) - dedxmmy = 0. !(1) + dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) + dedxmmy = 0. else !<---- Y view - dedxmmx = 0. !(1) - dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) + dedxmmx = 0. + dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) 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 !<<<< -* ---------------------------- -c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' +c QUIQUI------------ +c anche qui: non ci vuole clinc??? + if(iclm.ne.0)then if(mod(VIEW(iclm),2).eq.0)then - XGOOD(nplanes-ip+1)=1. - resx(nplanes-ip+1)=rxmm - if(DEBUG)print*,'%%%% included X-cl ',iclm -c if(.true.)print*,'%%%% included X-cl ',iclm - $ ,'( chi^2, ',RCHI2_STORE(ibest) - $ ,', norm.dist.= ',distmin - $ ,', cut ',clinc,' )' - else - YGOOD(nplanes-ip+1)=1. - resy(nplanes-ip+1)=rymm - if(DEBUG)print*,'%%%% included Y-cl ',iclm -c if(.true.)print*,'%%%% included Y-cl ',iclm - $ ,'( chi^2, ',RCHI2_STORE(ibest) - $ ,', norm.dist.= ', distmin - $ ,', cut ',clinc,' )' + clincnew= + $ 20* + $ sqrt(rxmm**2+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)) + endif +c QUIQUI------------ + + if(distmin.le.clincnew)then !QUIQUI +c if(distmin.le.clinc)then !QUIQUI + + CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< +* ---------------------------- +c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + if(mod(VIEW(iclm),2).eq.0)then + XGOOD(nplanes-ip+1)=1. + resx(nplanes-ip+1)=rxmm + if(DEBUG)print*,'%%%% included X-cl ',iclm + $ ,'( chi^2, ',RCHI2_STORE(ibest) + $ ,', dist.= ',distmin + $ ,', cut ',clinc,' )' + else + YGOOD(nplanes-ip+1)=1. + resy(nplanes-ip+1)=rymm + if(DEBUG)print*,'%%%% included Y-cl ',iclm + $ ,'( chi^2, ',RCHI2_STORE(ibest) + $ ,', dist.= ', distmin + $ ,', cut ',clinc,' )' + endif +c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' +* ---------------------------- + xm_A(nplanes-ip+1) = xmm_A + ym_A(nplanes-ip+1) = ymm_A + xm_B(nplanes-ip+1) = xmm_B + ym_B(nplanes-ip+1) = ymm_B + zm(nplanes-ip+1) = (zmm_A+zmm_B)/2. + dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< + dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< +* ---------------------------- endif -c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' -* ---------------------------- - xm_A(nplanes-ip+1) = xmm_A - ym_A(nplanes-ip+1) = ymm_A - xm_B(nplanes-ip+1) = xmm_B - ym_B(nplanes-ip+1) = ymm_B - zm(nplanes-ip+1) = (zmm_A+zmm_B)/2. -c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1) - dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1) - dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1) -* ---------------------------- endif endif 133 continue @@ -2961,19 +3438,13 @@ * * * * ************************************************** -cccccc 12/08/2006 modified by elena ---> (1) * subroutine clean_XYclouds(ibest,iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' -c include 'momanhough_init.f' - include 'level2.f' !(1) -c include 'calib.f' -c include 'level1.f' - - + include 'level2.f' do ip=1,nplanes !loop on planes @@ -2983,18 +3454,12 @@ if(id.ne.0)then iclx=clx(ip,icp_cp(id)) icly=cly(ip,icp_cp(id)) -c cl_used(iclx)=1 !tag used clusters -c cl_used(icly)=1 !tag used clusters - cl_used(iclx)=ntrk !tag used clusters !(1) - cl_used(icly)=ntrk !tag used clusters !(1) + cl_used(iclx)=ntrk !tag used clusters + cl_used(icly)=ntrk !tag used clusters elseif(icl.ne.0)then -c cl_used(icl)=1 !tag used clusters - cl_used(icl)=ntrk !tag used clusters !1) + cl_used(icl)=ntrk !tag used clusters endif -c if(DEBUG)then -c print*,ip,' <<< ',id -c endif * ----------------------------- * remove the couple from clouds * remove also vitual couples containing the @@ -3055,13 +3520,22 @@ include 'level1.f' include 'common_momanhough.f' include 'level2.f' -c include 'level1.f' +* --------------------------------- +* variables initialized from level1 +* --------------------------------- do i=1,nviews good2(i)=good1(i) + do j=1,nva1_view + vkflag(i,j)=1 + if(cnnev(i,j).le.0)then + vkflag(i,j)=cnnev(i,j) + endif + enddo enddo - - +* ---------------- +* level2 variables +* ---------------- NTRK = 0 do it=1,NTRKMAX IMAGE(IT)=0 @@ -3072,8 +3546,13 @@ ZM_nt(IP,IT) = 0 RESX_nt(IP,IT) = 0 RESY_nt(IP,IT) = 0 + TAILX_nt(IP,IT) = 0 + TAILY_nt(IP,IT) = 0 + XBAD(IP,IT) = 0 + YBAD(IP,IT) = 0 XGOOD_nt(IP,IT) = 0 YGOOD_nt(IP,IT) = 0 + LS(IP,IT) = 0 DEDX_X(IP,IT) = 0 DEDX_Y(IP,IT) = 0 CLTRX(IP,IT) = 0 @@ -3206,17 +3685,23 @@ include 'commontracker.f' -c include 'level1.f' include 'level1.f' include 'common_momanhough.f' include 'level2.f' include 'common_mini_2.f' - real sinth,phi,pig + include 'calib.f' + + character*10 PFA + common/FINALPFA/PFA + + real sinth,phi,pig + integer ssensor,sladder pig=acos(-1.) +* ------------------------------------- chi2_nt(ntr) = sngl(chi2) nstep_nt(ntr) = nstep - +* ------------------------------------- phi = al(4) sinth = al(3) if(sinth.lt.0)then @@ -3229,14 +3714,13 @@ $ phi = phi + 2*pig al(4) = phi al(3) = sinth - do i=1,5 al_nt(i,ntr) = sngl(al(i)) do j=1,5 coval(i,j,ntr) = sngl(cov(i,j)) enddo enddo - +* ------------------------------------- do ip=1,nplanes ! loop on planes xgood_nt(ip,ntr) = int(xgood(ip)) ygood_nt(ip,ntr) = int(ygood(ip)) @@ -3245,29 +3729,91 @@ zm_nt(ip,ntr) = sngl(zm(ip)) RESX_nt(IP,ntr) = sngl(resx(ip)) RESY_nt(IP,ntr) = sngl(resy(ip)) + TAILX_nt(IP,ntr) = 0. + TAILY_nt(IP,ntr) = 0. xv_nt(ip,ntr) = sngl(xv(ip)) 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( + $ tan( acos(-1.) * sngl(axv(ip)) /180. )**2 + + $ tan( 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) + ax = axv_nt(ip,ntr) + 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 print*,'* ',ip,bfx,bfy,angx,angy + + id = CP_STORE(ip,IDCAND) ! couple id icl = CLS_STORE(ip,IDCAND) + ssensor = -1 + sladder = -1 + ssensor = SENSOR_STORE(ip,IDCAND) + sladder = LADDER_STORE(ip,IDCAND) + 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 >>> is a couple cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id)) cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id)) -c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr) + +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))) + + + if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0) + $ dedx_x(ip,ntr)=-dedx_x(ip,ntr) + if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0) + $ dedx_y(ip,ntr)=-dedx_y(ip,ntr) + elseif(icl.ne.0)then - if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl - if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl -c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr) + + 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) + 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) + endif + endif enddo +c$$$ print*,(xgood(i),i=1,6) +c$$$ print*,(ygood(i),i=1,6) +c$$$ print*,(ls(i,ntr),i=1,6) +c$$$ print*,(dedx_x(i,ntr),i=1,6) +c$$$ print*,(dedx_y(i,ntr),i=1,6) +c$$$ print*,'-----------------------' + end subroutine fill_level2_siglets @@ -3279,7 +3825,6 @@ * ------------------------------------------------------- include 'commontracker.f' -c include 'level1.f' include 'calib.f' include 'level1.f' include 'common_momanhough.f' @@ -3287,12 +3832,12 @@ include 'common_xyzPAM.f' * count #cluster per plane not associated to any track -c good2=1!.true. nclsx = 0 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 @@ -3301,11 +3846,13 @@ 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)) + if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx) 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 @@ -3316,11 +3863,13 @@ 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)) + if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy) 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 @@ -3330,7 +3879,6 @@ c$$$ print*,'ys(2,nclsy) ',ys(2,nclsy) endif endif -c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip) ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO whichtrack(icl) = cl_used(icl)