--- DarthVader/TrackerLevel2/src/F77/functionspfa.f 2008/03/05 17:00:20 1.23 +++ DarthVader/TrackerLevel2/src/F77/functionspfa.f 2008/03/22 08:32:50 1.24 @@ -8,6 +8,9 @@ * * integer function npfastrips(ic,angle) * +* ----------------------------------------------------------------- +* p.f.a. +* ----------------------------------------------------------------- * real function pfaeta(ic,angle) * real function pfaetal(ic,angle) * real function pfaeta2(ic,angle) @@ -15,17 +18,32 @@ * real function pfaeta4(ic,angle) * real function cog(ncog,ic) * +* ----------------------------------------------------------------- +* risoluzione spaziale media, stimata dalla simulazione (samuele) +* ----------------------------------------------------------------- +* FUNCTION risxeta2(angle) +* FUNCTION risxeta3(angle) +* FUNCTION risxeta4(angle) +* FUNCTION risyeta2(angle) +* FUNCTION risy_cog(angle) +* FUNCTION risx_cog(angle) +* real function riseta(iview,angle) +* ----------------------------------------------------------------- +* fattore moltiplicativo per tenere conto della dipendenza della +* risoluzione dal rumore delle strip +* ----------------------------------------------------------------- * real function fbad_cog(ncog,ic) * real function fbad_eta(ic,angle) * -* real function riseta(iview,angle) -* FUNCTION risxeta2(x) -* FUNCTION risxeta3(x) -* FUNCTION risxeta4(x) -* FUNCTION risyeta2(x) -* FUNCTION risy_cog(x) -* FUNCTION risx_cog(x) +* ----------------------------------------------------------------- +* NUOVO APPROCCIO PER LA STIMA DELLA RISOLUZIONE +* ----------------------------------------------------------------- +* real function riscogtheor(ncog,ic) +* real function risetatheor(ncog,ic,angle) * +* ----------------------------------------------------------------- +* correzione landi +* ----------------------------------------------------------------- * real function pfacorr(ic,angle) * * real function effectiveangle(ang,iview,bbb) @@ -552,6 +570,13 @@ * resolution. * It calls the function FBAD_COG(NCOG,IC), * accordingto the angle +* +* >>> cosi` non e` corretto!! +* >>> l'errore sulla coordinata eta si ottiene moltiplicando +* >>> l'errore sulla coordinata cog per la derivata della +* >>> distribuzione eta... pur sapendolo l'ho sempre ignorato... +* >>> deve essere modificato!!!! +* *------------------------------------------------------- include 'commontracker.f' @@ -589,7 +614,7 @@ end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** - real function pfaeta2(ic,angle) !(1) + real function pfaeta2(ic,angle) *-------------------------------------------------------------- * this function returns * @@ -608,10 +633,10 @@ real cog2,angle integer iview,lad - iview = VIEW(ic) - lad = nld(MAXS(ic),VIEW(ic)) - cog2 = cog(2,ic) - pfaeta2=cog2 + iview = VIEW(ic) + lad = nld(MAXS(ic),VIEW(ic)) + cog2 = cog(2,ic) + pfaeta2 = cog2 * ---------------- * find angular bin @@ -1311,116 +1336,351 @@ end -c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** -c$$$ real function fbad_cog0(ncog,ic) -c$$$*------------------------------------------------------- -c$$$* this function returns a factor that takes into -c$$$* account deterioration of the spatial resolution -c$$$* in the case BAD strips are included in the cluster. -c$$$* This factor should multiply the nominal spatial -c$$$* resolution. -c$$$* -c$$$* NB!!! -c$$$* (this is the old version. It consider only the two -c$$$* strips with the greatest signal. The new one is -c$$$* fbad_cog(ncog,ic) ) -c$$$* -c$$$*------------------------------------------------------- -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'level1.f' -c$$$ include 'calib.f' -c$$$ -c$$$* --> signal of the central strip -c$$$ sc = CLSIGNAL(INDMAX(ic)) !center -c$$$ -c$$$* signal of adjacent strips -c$$$* --> left -c$$$ sl1 = 0 !left 1 -c$$$ if( -c$$$ $ (INDMAX(ic)-1).ge.INDSTART(ic) -c$$$ $ ) -c$$$ $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1)) -c$$$ -c$$$ sl2 = 0 !left 2 -c$$$ if( -c$$$ $ (INDMAX(ic)-2).ge.INDSTART(ic) -c$$$ $ ) -c$$$ $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2)) -c$$$ -c$$$* --> right -c$$$ sr1 = 0 !right 1 -c$$$ if( -c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) -c$$$ $ .or. -c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) -c$$$ $ ) -c$$$ $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1)) -c$$$ -c$$$ sr2 = 0 !right 2 -c$$$ if( -c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) -c$$$ $ .or. -c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) -c$$$ $ ) -c$$$ $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2)) -c$$$ -c$$$ -c$$$ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view -c$$$ f = 4. -c$$$ si = 8.4 -c$$$ else !X-view -c$$$ f = 6. -c$$$ si = 3.9 -c$$$ endif -c$$$ -c$$$ fbad_cog = 1. -c$$$ f0 = 1 -c$$$ f1 = 1 -c$$$ f2 = 1 -c$$$ f3 = 1 -c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then -c$$$ -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f -c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f -c$$$ -c$$$ if(ncog.eq.2.and.sl1.ne.0)then -c$$$ fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.) -c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then -c$$$ fbad_cog = 1. -c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then -c$$$ fbad_cog = 1. -c$$$ else -c$$$ fbad_cog = 1. -c$$$ endif -c$$$ -c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then -c$$$ -c$$$ -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f -c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f -c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f -c$$$ -c$$$ if(ncog.eq.2.and.sr1.ne.0)then -c$$$ fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.) -c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then -c$$$ fbad_cog = 1. -c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then -c$$$ fbad_cog = 1. -c$$$ else -c$$$ fbad_cog = 1. -c$$$ endif -c$$$ -c$$$ endif -c$$$ -c$$$ fbad_cog0 = sqrt(fbad_cog) -c$$$ -c$$$ return -c$$$ end -c$$$ -c$$$ -c$$$ +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** + + real function riscogtheor(ncog,ic) +*------------------------------------------------------- +* +* this function returns the expected resolution +* obtained by propagating the strip noise +* to the center-of-gravity coordinate +* +* ncog = n.strip used in the coordinate evaluation +* (ncog=0 => all strips above threshold) +* +*------------------------------------------------------- + + include 'commontracker.f' + include 'level1.f' + include 'calib.f' + + if(mod(int(VIEW(ic)),2).eq.1)then !Y-view + incut = incuty + pitch = pitchY / 1.e4 + else !X-view + incut = incutx + pitch = pitchX / 1.e4 + endif + + func = 100000. + stot = 0. + + if (ncog.gt.0) then + +* --> signal of the central strip + sc = CLSIGNAL(INDMAX(ic)) !center + fsc = clsigma(INDMAX(ic)) +* --> signal of adjacent strips + sl1 = 0 !left 1 + fsl1 = 1 !left 1 + if( + $ (INDMAX(ic)-1).ge.INDSTART(ic) + $ )then + sl1 = CLSIGNAL(INDMAX(ic)-1) + fsl1 = clsigma(INDMAX(ic)-1) + endif + + sl2 = 0 !left 2 + fsl2 = 1 !left 2 + if( + $ (INDMAX(ic)-2).ge.INDSTART(ic) + $ )then + sl2 = CLSIGNAL(INDMAX(ic)-2) + fsl2 = clsigma(INDMAX(ic)-2) + endif + sr1 = 0 !right 1 + fsr1 = 1 !right 1 + if( + $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) + $ .or. + $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) + $ )then + sr1 = CLSIGNAL(INDMAX(ic)+1) + fsr1 = clsigma(INDMAX(ic)+1) + endif + sr2 = 0 !right 2 + fsr2 = 1 !right 2 + if( + $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) + $ .or. + $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) + $ )then + sr2 = CLSIGNAL(INDMAX(ic)+2) + fsr2 = clsigma(INDMAX(ic)+2) + endif + + + +************************************************************ +* COG2-3-4 computation +************************************************************ + +c print*,sl2,sl1,sc,sr1,sr2 + + vCOG = cog(ncog,ic)!0. + + if(ncog.eq.1)then + func = 1./12. + stot = 1. + elseif(ncog.eq.2)then + if(sl1.gt.sr1)then + func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2) + stot = sl1+sc + elseif(sl1.le.sr1)then + func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) + stot = sc+sr1 + endif + elseif(ncog.eq.3)then + func = + $ (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) + stot = sl1+sc+sr1 + elseif(ncog.eq.4)then + if(sl2.gt.sr2)then + func = + $ (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2 + $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) + stot = sl2+sl1+sc+sr1 + elseif(sl2.le.sr2)then + func = + $ (fsl1*(-1-vCOG)**2 + $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2) + stot = sl2+sl1+sc+sr1 + endif + else + print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG + $ ,' not implemented' + endif + + elseif(ncog.eq.0)then +* ========================= +* COG computation +* ========================= + + vCOG = cog(0,ic) + + iv = VIEW(ic) + istart = INDSTART(IC) + istop = TOTCLLENGTH + if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1 +ccc SGN = 0. + SNU = 0. +ccc SDE = 0. + + do i=INDMAX(IC),istart,-1 + ipos = i-INDMAX(ic) + cut = incut*CLSIGMA(i) + if(CLSIGNAL(i).gt.cut)then + fs = clsigma(i) + SNU = SNU + fs*(ipos-vCOG)**2 + stot = stot + CLSIGNAL(i) + else + goto 10 + endif + enddo + 10 continue + do i=INDMAX(IC)+1,istop + ipos = i-INDMAX(ic) + cut = incut*CLSIGMA(i) + if(CLSIGNAL(i).gt.cut)then + fs = clsigma(i) + SNU = SNU + fs*(ipos-vCOG)**2 + stot = stot + CLSIGNAL(i) + else + goto 20 + endif + enddo + 20 continue + if(SDE.ne.0)then + FUNC=SNU + else + + endif + + else + + FUNC=0 + print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG + print*,' (NCOG must be >= 0)' + + + endif + + + if(stot.gt.0..and.func.gt.0.)then + func = sqrt(func) + func = pitch * func/stot + endif + + riscogtheor = func + + return + end + + + +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** + + real function risetatheor(ncog,ic,angle) +*------------------------------------------------------- +* +* this function returns the expected resolution +* obtained by propagating the strip noise +* to the coordinate evaluated with non-linear eta-function +* +* ncog = n.strip used in the coordinate evaluation +* (ncog=0 => ncog=2,3,4 according to angle) +* +*------------------------------------------------------- + + include 'commontracker.f' + include 'level1.f' + include 'calib.f' + + + func = 1. + + iview = VIEW(ic) + lad = nld(MAXS(ic),VIEW(ic)) + +* ------------------------------------------------ +* number of strip to be used (in case of ncog = 0) +* ------------------------------------------------ + + inoeta = 0 + + if(mod(int(iview),2).eq.1)then !Y-view + + pitch = pitchY / 1.e4 + + if(ncog.eq.0)then + if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then + ncog = 2 + elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then + ncog = 3 + elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then + ncog = 4 + else + ncog = 4 + inoeta = 1 + endif + endif + + else !X-view + + pitch = pitchX / 1.e4 + + if(ncog.eq.0)then + if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then + ncog = 2 + elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then + ncog = 3 + elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then + ncog = 4 + else + ncog = 4 + inoeta = 1 + endif + endif + + endif + + func = riscogtheor(ncog,ic) + + risetatheor = func + + if(inoeta.eq.1)return ! no eta correction is applied --> exit + if(ncog.lt.1.or.ncog.gt.4)return + +* ---------------- +* find angular bin +* ---------------- +* (in futuro possiamo pensare di interpolare anche sull'angolo) + do iang=1,nangbin + if(angL(iang).lt.angle.and.angR(iang).ge.angle)then + iangle=iang + goto 98 + endif + enddo + if(DEBUG.EQ.1)print* + $ ,'risetatheor *** warning *** angle out of range: ',angle + if(angle.le.angL(1))iang=1 + if(angle.ge.angR(nangbin))iang=nangbin + 98 continue !jump here if ok + +* ------------- +* within +/-0.5 +* ------------- + + vcog = cog(ncog,ic) + + etamin = eta2(1,iang) + etamax = eta2(netaval,iang) + + iaddmax=10 + iadd=0 + 10 continue + if(vcog.lt.etamin)then + vcog = vcog + 1 + iadd = iadd + 1 + if(iadd>iaddmax)goto 111 + goto 10 + endif + 20 continue + if(vcog.gt.etamax)then + vcog = vcog - 1 + iadd = iadd - 1 + if(iadd<-1*iaddmax)goto 111 + goto 20 + endif + goto 1111 + 111 continue + if(DEBUG.eq.1) + $ print*,'risetatheor *** warning *** anomalous cluster' + if(DEBUG.eq.1) + $ print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)' + vcog=0 + 1111 continue + +* ------------------------------------------------ +* interpolation +* ------------------------------------------------ + + + ibin = netaval + do i=2,netaval + if(ncog.eq.2)eta=eta2(i,iang) + if(ncog.eq.3)eta=eta3(i,iang) + if(ncog.eq.4)eta=eta4(i,iang) + if(eta.ge.vcog)then + ibin = i + goto 99 + endif + enddo + 99 continue + + if(ncog.eq.2)then + x1 = eta2(ibin-1,iang) + x2 = eta2(ibin,iang) + y1 = feta2(ibin-1,iview,lad,iang) + y2 = feta2(ibin,iview,lad,iang) + elseif(ncog.eq.3)then + x1 = eta3(ibin-1,iang) + x2 = eta3(ibin,iang) + y1 = feta3(ibin-1,iview,lad,iang) + y2 = feta3(ibin,iview,lad,iang) + elseif(ncog.eq.4)then + x1 = eta4(ibin-1,iang) + x2 = eta4(ibin,iang) + y1 = feta4(ibin-1,iview,lad,iang) + y2 = feta4(ibin,iview,lad,iang) + endif + + func = func * (y2-y1)/(x2-x1) + + risetatheor = func + + return + end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** @@ -1948,8 +2208,8 @@ pfacorr = fcorr(iview,lad,iang) - if(DEBUG.eq.1)print*,'CORR (ic ',ic,' ang',angle,') -->',pfacorr - + if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr + 100 return end