--- DarthVader/TrackerLevel2/src/F77/functionspfa.f 2007/08/21 15:21:22 1.19 +++ DarthVader/TrackerLevel2/src/F77/functionspfa.f 2007/08/27 12:57:15 1.20 @@ -1,4 +1,35 @@ - +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** +* this file contains all subroutines and functions +* that are needed for position finding algorithms: +* +* subroutine idtoc(ipfa,cpfa) +* +* integer function npfastrips(ic,angle) +* +* real function pfaeta(ic,angle) +* real function pfaetal(ic,angle) +* real function pfaeta2(ic,angle) +* real function pfaeta3(ic,angle) +* real function pfaeta4(ic,angle) +* real function cog(ncog,ic) +* +* 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) +* +* real function pfacorr(ic,angle) +* +* NB - The angle is the "effective angle", which is relative +* to the sensor and it takes into account the magnetic field +* +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** subroutine idtoc(ipfa,cpfa) @@ -19,15 +50,9 @@ end -*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** -* this file contains all subroutines and functions -* that are needed for position finding algorithms -* -* -*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** - integer function npfastrips(ic,PFA,angle) + integer function npfastrips(ic,angle) *-------------------------------------------------------------- * thid function returns the number of strips used * to evaluate the position of a cluster, according to the p.f.a. @@ -36,12 +61,13 @@ include 'level1.f' include 'calib.f' - character*4 usedPFA,PFA + character*4 usedPFA + - usedPFA=PFA + call idtoc(pfaid,usedPFA) - npfastrips=0 + npfastrips=-1 if(usedPFA.eq.'COG1')npfastrips=1 if(usedPFA.eq.'COG2')npfastrips=2 @@ -51,7 +77,7 @@ if(usedPFA.eq.'ETA3')npfastrips=3 if(usedPFA.eq.'ETA4')npfastrips=4 * ---------------------------------------------------------------- - if(usedPFA.eq.'ETA')then + if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then c print*,VIEW(ic),angle if(mod(int(VIEW(ic)),2).eq.1)then !Y-view if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then @@ -61,8 +87,7 @@ elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then npfastrips=4 else - npfastrips=4 -c usedPFA='COG' + npfastrips=4 !COG4 endif else !X-view if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then @@ -72,49 +97,50 @@ elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then npfastrips=4 else - npfastrips=4 -c usedPFA='COG' + npfastrips=4 !COG4 endif endif endif * ---------------------------------------------------------------- if(usedPFA.eq.'COG')then - iv=VIEW(ic) - if(mod(iv,2).eq.1)incut=incuty - if(mod(iv,2).eq.0)incut=incutx - istart = INDSTART(IC) - istop = TOTCLLENGTH - if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 - mu = 0 - do i = INDMAX(IC),istart,-1 - ipos = i-INDMAX(ic) - cut = incut*CLSIGMA(i) - if(CLSIGNAL(i).ge.cut)then - mu = mu + 1 - print*,i,mu - else - goto 10 - endif - enddo - 10 continue - do i = INDMAX(IC)+1,istop - ipos = i-INDMAX(ic) - cut = incut*CLSIGMA(i) - if(CLSIGNAL(i).ge.cut)then - mu = mu + 1 - print*,i,mu - else - goto 20 - endif - enddo - 20 continue - npfastrips=mu + npfastrips=0 + +c$$$ iv=VIEW(ic) +c$$$ if(mod(iv,2).eq.1)incut=incuty +c$$$ if(mod(iv,2).eq.0)incut=incutx +c$$$ istart = INDSTART(IC) +c$$$ istop = TOTCLLENGTH +c$$$ if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 +c$$$ mu = 0 +c$$$ do i = INDMAX(IC),istart,-1 +c$$$ ipos = i-INDMAX(ic) +c$$$ cut = incut*CLSIGMA(i) +c$$$ if(CLSIGNAL(i).ge.cut)then +c$$$ mu = mu + 1 +c$$$ print*,i,mu +c$$$ else +c$$$ goto 10 +c$$$ endif +c$$$ enddo +c$$$ 10 continue +c$$$ do i = INDMAX(IC)+1,istop +c$$$ ipos = i-INDMAX(ic) +c$$$ cut = incut*CLSIGMA(i) +c$$$ if(CLSIGNAL(i).ge.cut)then +c$$$ mu = mu + 1 +c$$$ print*,i,mu +c$$$ else +c$$$ goto 20 +c$$$ endif +c$$$ enddo +c$$$ 20 continue +c$$$ npfastrips=mu endif * ---------------------------------------------------------------- -c print*,pfastrips +c print*,pfaid,usedPFA,angle,npfastrips return end @@ -137,12 +163,12 @@ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view - if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then + if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then pfaeta = pfaeta2(ic,angle) cc print*,pfaeta2(ic,angle) - elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then + elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then pfaeta = pfaeta3(ic,angle) - elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then + elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then pfaeta = pfaeta4(ic,angle) else pfaeta = cog(4,ic) @@ -150,11 +176,11 @@ else !X-view - if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then + if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then pfaeta = pfaeta2(ic,angle) - elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then + elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then pfaeta = pfaeta3(ic,angle) - elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then + elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then pfaeta = pfaeta4(ic,angle) else pfaeta = cog(4,ic) @@ -183,12 +209,12 @@ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view - if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then + if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle) cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle) - elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then + elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle) - elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then + elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle) else pfaetal = cog(4,ic) @@ -196,12 +222,12 @@ else !X-view - if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then + if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle) cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle) - elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then + elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle) - elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then + elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle) else pfaetal = cog(4,ic) @@ -659,103 +685,6 @@ *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** -c$$$ real function cog0(ncog,ic) -c$$$*------------------------------------------------- -c$$$* this function returns -c$$$* -c$$$* - the Center-Of-Gravity of the cluster IC -c$$$* evaluated using NCOG strips, -c$$$* calculated relative to MAXS(IC) -c$$$* -c$$$* - zero in case that not enough strips -c$$$* have a positive signal -c$$$* -c$$$* NOTE: -c$$$* This is the old definition, used by Straulino. -c$$$* The new routine, according to Landi, -c$$$* is COG(NCOG,IC) -c$$$*------------------------------------------------- -c$$$ -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'level1.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$$$* COG computation -c$$$************************************************************ -c$$$ -c$$$c print*,sl2,sl1,sc,sr1,sr2 -c$$$ -c$$$ COG = 0. -c$$$ -c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then -c$$$ -c$$$ if(ncog.eq.2.and.sl1.ne.0)then -c$$$ COG = -sl1/(sl1+sc) -c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then -c$$$ COG = (sr1-sl1)/(sl1+sc+sr1) -c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then -c$$$ COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) -c$$$ else -c$$$ COG = 0. -c$$$ endif -c$$$ -c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then -c$$$ -c$$$ if(ncog.eq.2.and.sr1.ne.0)then -c$$$ COG = sr1/(sc+sr1) -c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then -c$$$ COG = (sr1-sl1)/(sl1+sc+sr1) -c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then -c$$$ COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) -c$$$ else -c$$$ COG = 0. -c$$$ endif -c$$$ -c$$$ endif -c$$$ -c$$$ COG0 = COG -c$$$ -c$$$c print *,ncog,ic,cog,'/////////////' -c$$$ -c$$$ return -c$$$ end - -*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function cog(ncog,ic) *------------------------------------------------- * this function returns @@ -1072,36 +1001,6 @@ SGN = 0. SNU = 0. SDE = 0. -c$$$ do i=INDMAX(IC),istart,-1 -c$$$ ipos = i-INDMAX(ic) -c$$$ cut = incut*CLSIGMA(i) -c$$$ if(CLSIGNAL(i).gt.cut)then -c$$$ COG = COG + ipos*CLSIGNAL(i) -c$$$ SGN = SGN + CLSIGNAL(i) -c$$$ else -c$$$ goto 10 -c$$$ endif -c$$$ enddo -c$$$ 10 continue -c$$$ do i=INDMAX(IC)+1,istop -c$$$ ipos = i-INDMAX(ic) -c$$$ cut = incut*CLSIGMA(i) -c$$$ if(CLSIGNAL(i).gt.cut)then -c$$$ COG = COG + ipos*CLSIGNAL(i) -c$$$ SGN = SGN + CLSIGNAL(i) -c$$$ else -c$$$ goto 20 -c$$$ endif -c$$$ enddo -c$$$ 20 continue -c$$$ if(SGN.le.0)then -c$$$ print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN -c$$$ print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) -c$$$ print*,(CLSIGNAL(i),i=istart,istop) -c$$$ print*,'fbad_cog(0,ic) --> NOT EVALUATED ' -c$$$ else -c$$$ COG=COG/SGN -c$$$ endif do i=INDMAX(IC),istart,-1 ipos = i-INDMAX(ic) @@ -1149,116 +1048,116 @@ end -*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** - real function fbad_cog0(ncog,ic) -*------------------------------------------------------- -* this function returns a factor that takes into -* account deterioration of the spatial resolution -* in the case BAD strips are included in the cluster. -* This factor should multiply the nominal spatial -* resolution. -* -* NB!!! -* (this is the old version. It consider only the two -* strips with the greatest signal. The new one is -* fbad_cog(ncog,ic) ) -* -*------------------------------------------------------- - - include 'commontracker.f' - include 'level1.f' - include 'calib.f' - -* --> signal of the central strip - sc = CLSIGNAL(INDMAX(ic)) !center - -* signal of adjacent strips -* --> left - sl1 = 0 !left 1 - if( - $ (INDMAX(ic)-1).ge.INDSTART(ic) - $ ) - $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1)) - - sl2 = 0 !left 2 - if( - $ (INDMAX(ic)-2).ge.INDSTART(ic) - $ ) - $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2)) - -* --> right - sr1 = 0 !right 1 - if( - $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) - $ .or. - $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) - $ ) - $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1)) - - sr2 = 0 !right 2 - if( - $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) - $ .or. - $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) - $ ) - $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2)) - - - if(mod(int(VIEW(ic)),2).eq.1)then !Y-view - f = 4. - si = 8.4 - else !X-view - f = 6. - si = 3.9 - endif - - fbad_cog = 1. - f0 = 1 - f1 = 1 - f2 = 1 - f3 = 1 - if(sl1.gt.sr1.and.sl1.gt.0.)then - - if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f - if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f -c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f - - if(ncog.eq.2.and.sl1.ne.0)then - fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.) - elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then - fbad_cog = 1. - elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then - fbad_cog = 1. - else - fbad_cog = 1. - endif - - elseif(sl1.le.sr1.and.sr1.gt.0.)then - - - if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f - if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f -c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f - - if(ncog.eq.2.and.sr1.ne.0)then - fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.) - elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then - fbad_cog = 1. - elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then - fbad_cog = 1. - else - fbad_cog = 1. - endif - - endif - - fbad_cog0 = sqrt(fbad_cog) - - return - 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$$$ *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***