--- DarthVader/TrackerLevel2/src/F77/functionspfa.f 2007/08/27 12:57:15 1.20 +++ DarthVader/TrackerLevel2/src/F77/functionspfa.f 2007/08/31 14:56:52 1.21 @@ -4,6 +4,8 @@ * * subroutine idtoc(ipfa,cpfa) * +* subroutine applypfa(PFAtt,ic,ang,corr,res) +* * integer function npfastrips(ic,angle) * * real function pfaeta(ic,angle) @@ -26,6 +28,9 @@ * * real function pfacorr(ic,angle) * +* real function effectiveangle(ang,iview,bbb) +* real function fieldcorr(iview,bbb) +* * NB - The angle is the "effective angle", which is relative * to the sensor and it takes into account the magnetic field * @@ -49,9 +54,257 @@ if(ipfa.eq.14)CPFA='COG4' end +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** + real function effectiveangle(ang,iview,bbb) + + include 'commontracker.f' + + effectiveangle = 0. + + if(mod(iview,2).eq.0)then +c ================================================= +c X view +c ================================================= +c here bbb is the y component of the m.field + angx = ang + by = bbb + if(iview.eq.12) angx = -1. * ang + if(iview.eq.12) by = -1. * bbb + tgtemp = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 + + elseif(mod(iview,2).eq.1)then +c ================================================= +c Y view +c ================================================= +c here bbb is the x component of the m.filed + angy = ang + bx = bbb + tgtemp = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001 + + endif + effectiveangle = 180.*atan(tgtemp)/acos(-1.) + + return + end +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** + real function fieldcorr(iview,bbb) + + include 'commontracker.f' + + fieldcorr = 0. + + if(mod(iview,2).eq.0)then + +c ================================================= +c X view +c ================================================= +c here bbb is the y component of the m.field + by = bbb + if(iview.eq.12) by = -1. * bbb + fieldcorr = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX + + elseif(mod(iview,2).eq.1)then +c ================================================= +c Y view +c ================================================= +c here bbb is the x component of the m.filed + bx = bbb + fieldcorr = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY + + endif + + return + end +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** + + subroutine applypfa(PFAtt,ic,ang,corr,res) +*--------------------------------------------------------------- +* this subroutine calculate the coordinate of cluster ic (in +* strip units), relative to the strip with the maximum signal, +* and its spatial resolution (in cm), applying PFAtt. +* ang is the effective angle, relative to the sensor +*--------------------------------------------------------------- + + character*4 PFAtt + include 'commontracker.f' + include 'level1.f' + + corr = 0 + res = 0 + + if(ic.le.0)return + + iview = VIEW(ic) + + if(mod(iview,2).eq.0)then +c ================================================= +c X view +c ================================================= + + res = RESXAV + + if(PFAtt.eq.'COG1')then + + corr = 0 + res = 1e-4*pitchX/sqrt(12.)!!res + + elseif(PFAtt.eq.'COG2')then + + corr = cog(2,ic) + res = risx_cog(abs(ang))!TEMPORANEO + res = res*fbad_cog(2,ic) + + elseif(PFAtt.eq.'COG3')then + + corr = cog(3,ic) + res = risx_cog(abs(ang))!TEMPORANEO + res = res*fbad_cog(3,ic) + + elseif(PFAtt.eq.'COG4')then + + corr = cog(4,ic) + res = risx_cog(abs(ang))!TEMPORANEO + res = res*fbad_cog(4,ic) + + elseif(PFAtt.eq.'ETA2')then + + corr = pfaeta2(ic,ang) + res = risxeta2(abs(ang)) + res = res*fbad_cog(2,ic) + + elseif(PFAtt.eq.'ETA3')then + + corr = pfaeta3(ic,ang) + res = risxeta3(abs(ang)) + res = res*fbad_cog(3,ic) + + elseif(PFAtt.eq.'ETA4')then + + corr = pfaeta4(ic,ang) + res = risxeta4(abs(ang)) + res = res*fbad_cog(4,ic) + + elseif(PFAtt.eq.'ETA')then + corr = pfaeta(ic,ang) +c res = riseta(ic,ang) + res = riseta(iview,ang) + res = res*fbad_eta(ic,ang) + elseif(PFAtt.eq.'ETAL')then + corr = pfaetal(ic,ang) + res = riseta(iview,ang) + res = res*fbad_eta(ic,ang) + + elseif(PFAtt.eq.'COG')then + + corr = cog(0,ic) + res = risx_cog(abs(ang)) + res = res*fbad_cog(0,ic) + + else + if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt + endif + + +* ====================================== +* temporary patch for saturated clusters +* ====================================== + if( nsatstrips(ic).gt.0 )then + corr = cog(4,ic) + res = pitchX*1e-4/sqrt(12.) +cc cc=cog(4,ic) +c$$$ print*,ic,' *** ',cc +c$$$ print*,ic,' *** ',res + endif + + + elseif(mod(iview,2).eq.1)then +c ================================================= +c Y view +c ================================================= + + res = RESYAV + + if(PFAtt.eq.'COG1')then + + corr = 0 + res = 1e-4*pitchY/sqrt(12.)!res + + elseif(PFAtt.eq.'COG2')then + + corr = cog(2,ic) + res = risy_cog(abs(ang))!TEMPORANEO + res = res*fbad_cog(2,ic) + + elseif(PFAtt.eq.'COG3')then + + corr = cog(3,ic) + res = risy_cog(abs(ang))!TEMPORANEO + res = res*fbad_cog(3,ic) + + elseif(PFAtt.eq.'COG4')then + + corr = cog(4,ic) + res = risy_cog(abs(ang))!TEMPORANEO + res = res*fbad_cog(4,ic) + + elseif(PFAtt.eq.'ETA2')then + + corr = pfaeta2(ic,ang) + res = risyeta2(abs(ang)) + res = res*fbad_cog(2,ic) + + elseif(PFAtt.eq.'ETA3')then + + corr = pfaeta3(ic,ang) + res = res*fbad_cog(3,ic) + + elseif(PFAtt.eq.'ETA4')then + + corr = pfaeta4(ic,ang) + res = res*fbad_cog(4,ic) + + elseif(PFAtt.eq.'ETA')then + + corr = pfaeta(ic,ang) +c res = riseta(ic,ang) + res = riseta(iview,ang) + res = res*fbad_eta(ic,ang) + + elseif(PFAtt.eq.'ETAL')then + + corr = pfaetal(ic,ang) + res = riseta(iview,ang) + res = res*fbad_eta(ic,ang) + + elseif(PFAtt.eq.'COG')then + + corr = cog(0,ic) + res = risy_cog(abs(ang)) + res = res*fbad_cog(0,ic) + + else + if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt + endif + + +* ====================================== +* temporary patch for saturated clusters +* ====================================== + if( nsatstrips(ic).gt.0 )then + corr = cog(4,ic) + res = pitchY*1e-4/sqrt(12.) +cc cc=cog(4,ic) +c$$$ print*,ic,' *** ',cc +c$$$ print*,ic,' *** ',res + endif + + endif + end + +*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** integer function npfastrips(ic,angle) *-------------------------------------------------------------- * thid function returns the number of strips used @@ -749,43 +1002,72 @@ c ============================================================== if(ncog.eq.1)then COG = 0. - if(sr1.gt.sc)cog=1. !NEW - if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW + if(sr1.gt.sc)cog=1. + if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. c ============================================================== elseif(ncog.eq.2)then + COG = 0. if(sl1.gt.sr1)then if((sl1+sc).ne.0)COG = -sl1/(sl1+sc) elseif(sl1.lt.sr1)then if((sc+sr1).ne.0)COG = sr1/(sc+sr1) - elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW + elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1) - $ .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW + $ .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1) - $ .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW + $ .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) endif c if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic) c $ ,' : ',sl2,sl1,sc,sr1,sr2 c ============================================================== elseif(ncog.eq.3)then - if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1) -c if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic) + COG = 0 + sss = sc + if( sl1.ne.-9999. )COG = COG-sl1 + if( sl1.ne.-9999. )sss = sss+sl1 + if( sr1.ne.-9999. )COG = COG+sr1 + if( sr1.ne.-9999. )sss = sss+sr1 + if(sss.ne.0)COG=COG/sss + +c if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1) +c if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic) c $ ,' : ',sl2,sl1,sc,sr1,sr2 c ============================================================== elseif(ncog.eq.4)then + + COG = 0 + sss = sc + if( sl1.ne.-9999. )COG = COG-sl1 + if( sl1.ne.-9999. )sss = sss+sl1 + if( sr1.ne.-9999. )COG = COG+sr1 + if( sr1.ne.-9999. )sss = sss+sr1 if(sl2.gt.sr2)then - if((sl2+sl1+sc+sr1).ne.0) - $ COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) + if((sl2+sss).ne.0) + $ COG = (COG-2*sl2)/(sl2+sss) elseif(sl2.lt.sr2)then - if((sr2+sl1+sc+sr1).ne.0) - $ COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1) - elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW + if((sr2+sss).ne.0) + $ COG = (2*sr2+COG)/(sr2+sss) + elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2) - $ .and.(sl2+sl1+sc+sr1).ne.0 ) - $ cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW + $ .and.(sl2+sss).ne.0 ) + $ cog = (cog-2*sl2)/(sl2+sss) if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2) - $ .and.(sr2+sl1+sc+sr1).ne.0 ) - $ cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1) !NEW + $ .and.(sr2+sss).ne.0 ) + $ cog = (2*sr2+cog)/(sr2+sss) endif +c ============================================================== + elseif(ncog.eq.5)then + COG = 0 + sss = sc + if( sl1.ne.-9999. )COG = COG-sl1 + if( sl1.ne.-9999. )sss = sss+sl1 + if( sr1.ne.-9999. )COG = COG+sr1 + if( sr1.ne.-9999. )sss = sss+sr1 + if( sl2.ne.-9999. )COG = COG-2*sl2 + if( sl2.ne.-9999. )sss = sss+sl2 + if( sr2.ne.-9999. )COG = COG+2*sr2 + if( sr2.ne.-9999. )sss = sss+sr2 + if(sss.ne.0)COG=COG/sss else print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG $ ,' not implemented'