*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** * this file contains all subroutines and functions * that are needed for position finding algorithms: * * subroutine idtoc(ipfa,cpfa) * * subroutine applypfa(PFAtt,ic,ang,corr,res) * * integer function npfastrips(ic,angle) * * ----------------------------------------------------------------- * p.f.a. * ----------------------------------------------------------------- * 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) * * ----------------------------------------------------------------- * 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) * * ----------------------------------------------------------------- * 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) * 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 * *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** subroutine idtoc(ipfa,cpfa) integer ipfa character*10 cpfa CPFA='COG4' if(ipfa.eq.0)CPFA='ETA' if(ipfa.eq.2)CPFA='ETA2' if(ipfa.eq.3)CPFA='ETA3' if(ipfa.eq.4)CPFA='ETA4' if(ipfa.eq.5)CPFA='ETAL' if(ipfa.eq.10)CPFA='COG' if(ipfa.eq.11)CPFA='COG1' if(ipfa.eq.12)CPFA='COG2' if(ipfa.eq.13)CPFA='COG3' 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 cc tgtemp = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!! tgtemp = tan(angx*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 c corr = cog(4,ic) corr = digsat(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 c corr = cog(4,ic) corr = digsat(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 * to evaluate the position of a cluster, according to the p.f.a. *-------------------------------------------------------------- include 'commontracker.f' include 'level1.f' include 'calib.f' character*4 usedPFA call idtoc(pfaid,usedPFA) npfastrips=-1 if(usedPFA.eq.'COG1')npfastrips=1 if(usedPFA.eq.'COG2')npfastrips=2 if(usedPFA.eq.'COG3')npfastrips=3 if(usedPFA.eq.'COG4')npfastrips=4 if(usedPFA.eq.'ETA2')npfastrips=2 if(usedPFA.eq.'ETA3')npfastrips=3 if(usedPFA.eq.'ETA4')npfastrips=4 * ---------------------------------------------------------------- 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 npfastrips=2 elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then npfastrips=3 elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then npfastrips=4 else npfastrips=4 !COG4 endif else !X-view if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then npfastrips=2 elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then npfastrips=3 elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then npfastrips=4 else npfastrips=4 !COG4 endif endif endif * ---------------------------------------------------------------- if(usedPFA.eq.'COG')then 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*,pfaid,usedPFA,angle,npfastrips return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function pfaeta(ic,angle) *-------------------------------------------------------------- * this function returns the position (in strip units) * it calls: * - pfaeta2(ic,angle) * - pfaeta3(ic,angle) * - pfaeta4(ic,angle) * according to the angle *-------------------------------------------------------------- include 'commontracker.f' include 'level1.f' include 'calib.f' pfaeta = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view 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).lt.e3tay )then pfaeta = pfaeta3(ic,angle) elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then pfaeta = pfaeta4(ic,angle) else pfaeta = cog(4,ic) endif else !X-view if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then pfaeta = pfaeta2(ic,angle) elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then pfaeta = pfaeta3(ic,angle) elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then pfaeta = pfaeta4(ic,angle) else pfaeta = cog(4,ic) endif endif c 100 return return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function pfaetal(ic,angle) *-------------------------------------------------------------- * this function returns the position (in strip units) * it calls: * - pfaeta2(ic,angle)+pfcorr(ic,angle) * - pfaeta3(ic,angle)+pfcorr(ic,angle) * - pfaeta4(ic,angle)+pfcorr(ic,angle) * according to the angle *-------------------------------------------------------------- include 'commontracker.f' include 'level1.f' include 'calib.f' pfaetal = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view 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).lt.e3tay )then pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle) elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle) else pfaetal = cog(4,ic) endif else !X-view 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).lt.e3tax )then pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle) elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle) else pfaetal = cog(4,ic) endif endif c 100 return return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** c real function riseta(ic,angle) real function riseta(iview,angle) *-------------------------------------------------------------- * this function returns the average spatial resolution * (in cm) for the ETA algorithm (function pfaeta(ic,angle)) * it calls: * - risxeta2(angle) * - risyeta2(angle) * - risxeta3(angle) * - risxeta4(angle) * according to the angle *-------------------------------------------------------------- include 'commontracker.f' include 'level1.f' include 'calib.f' riseta = 0 c if(mod(int(VIEW(ic)),2).eq.1)then !Y-view if(mod(iview,2).eq.1)then !Y-view if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then riseta = risyeta2(angle) elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then riseta = risy_cog(angle) !ATTENZIONE!! elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then riseta = risy_cog(angle) !ATTENZIONE!! else riseta = risy_cog(angle) endif else !X-view if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then riseta = risxeta2(angle) elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then riseta = risxeta3(angle) elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then riseta = risxeta4(angle) else riseta = risx_cog(angle) endif endif c 100 return return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function fbad_eta(ic,angle) *------------------------------------------------------- * 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. * 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' include 'level1.f' include 'calib.f' fbad_eta = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then fbad_eta = fbad_cog(2,ic) elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then fbad_eta = fbad_cog(3,ic) elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then fbad_eta = fbad_cog(4,ic) else fbad_eta = fbad_cog(4,ic) endif else !X-view if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then fbad_eta = fbad_cog(2,ic) elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then fbad_eta = fbad_cog(3,ic) elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then fbad_eta = fbad_cog(4,ic) else fbad_eta = fbad_cog(4,ic) endif endif return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function pfaeta2(ic,angle) *-------------------------------------------------------------- * this function returns * * - the position (in strip units) * corrected according to the ETA2 Position Finding Algorithm. * The function performs an interpolation of FETA2%ETA2 * * - if the angle is out of range, the calibration parameters * of the lowest or higher bin are used * *-------------------------------------------------------------- include 'commontracker.f' include 'calib.f' include 'level1.f' real cog2,angle integer iview,lad iview = VIEW(ic) lad = nld(MAXS(ic),VIEW(ic)) cog2 = cog(2,ic) pfaeta2 = cog2 * ---------------- * 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*,'pfaeta2 *** 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 * ------------- iaddmax=10 iadd=0 10 continue if(cog2.lt.eta2(1,iang))then cog2 = cog2 + 1 iadd = iadd + 1 if(iadd>iaddmax)goto 111 goto 10 endif 20 continue if(cog2.gt.eta2(netaval,iang))then cog2 = cog2 - 1 iadd = iadd - 1 if(iadd<-1*iaddmax)goto 111 goto 20 endif goto 1111 111 continue if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster' if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)' cog2=0 1111 continue * -------------------------------- c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 do i=2,netaval if(eta2(i,iang).gt.cog2)then x1 = eta2(i-1,iang) x2 = eta2(i,iang) y1 = feta2(i-1,iview,lad,iang) y2 = feta2(i,iview,lad,iang) c print*,'*****',i,view,lad,iang c print*,'-----',x1,x2,y1,y2 goto 99 endif enddo 99 continue AA=(y2-y1)/(x2-x1) BB=y1-AA*x1 pfaeta2 = AA*cog2+BB pfaeta2 = pfaeta2 - iadd c$$$ if(iflag.eq.1)then c$$$ pfaeta2=pfaeta2-1. !temp c$$$ cog2=cog2-1. !temp c$$$ endif c$$$ if(iflag.eq.-1)then c$$$ pfaeta2=pfaeta2+1. !temp c$$$ cog2=cog2+1. !temp c$$$ endif if(DEBUG.EQ.1)print*,'ETA2 (ic ',ic,' ang',angle,')' $ ,cog2-iadd,' -->',pfaeta2 c 100 return return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function pfaeta3(ic,angle) !(1) *-------------------------------------------------------------- * this function returns * * - the position (in strip units) * corrected according to the ETA3 Position Finding Algorithm. * The function performs an interpolation of FETA3%ETA3 * * - if the angle is out of range, the calibration parameters * of the lowest or higher bin are used * *-------------------------------------------------------------- include 'commontracker.f' include 'calib.f' include 'level1.f' real cog3,angle integer iview,lad iview = VIEW(ic) lad = nld(MAXS(ic),VIEW(ic)) cog3 = cog(3,ic) cc = cog3 cog3 = cc pfaeta3=cog3 * ---------------- * find angular bin * ---------------- * (in futuro possiamo pensare di interpolare anche sull'angolo) do iang=1,nangbin c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle if(angL(iang).lt.angle.and.angR(iang).ge.angle)then iangle=iang goto 98 endif enddo if(DEBUG.EQ.1) $ print*,'pfaeta3 *** 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 * ------------- iaddmax=10 iadd=0 10 continue if(cog3.lt.eta3(1,iang))then cog3 = cog3 + 1. iadd = iadd + 1 if(iadd>iaddmax) goto 111 goto 10 endif 20 continue if(cog3.gt.eta3(netaval,iang))then cog3 = cog3 - 1. iadd = iadd - 1 if(iadd<-1*iaddmax) goto 111 goto 20 endif goto 1111 111 continue if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster' if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)' cog3=0 1111 continue * -------------------------------- c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 do i=2,netaval if(eta3(i,iang).gt.cog3)then x1 = eta3(i-1,iang) x2 = eta3(i,iang) y1 = feta3(i-1,iview,lad,iang) y2 = feta3(i,iview,lad,iang) c print*,'*****',i,view,lad,iang c print*,'-----',x1,x2,y1,y2 goto 99 endif enddo 99 continue AA=(y2-y1)/(x2-x1) BB=y1-AA*x1 pfaeta3 = AA*cog3+BB pfaeta3 = pfaeta3 - iadd if(DEBUG.EQ.1)print*,'ETA3 (ic ',ic,' ang',angle,')' $ ,cog3-iadd,' -->',pfaeta3 c 100 return return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function pfaeta4(ic,angle) *-------------------------------------------------------------- * this function returns * * - the position (in strip units) * corrected according to the ETA4 Position Finding Algorithm. * The function performs an interpolation of FETA3%ETA3 * * - if the angle is out of range, the calibration parameters * of the lowest or higher bin are used * *-------------------------------------------------------------- include 'commontracker.f' include 'calib.f' include 'level1.f' real cog4,angle integer iview,lad iview = VIEW(ic) lad = nld(MAXS(ic),VIEW(ic)) cog4=cog(4,ic) pfaeta4=cog4 * ---------------- * find angular bin * ---------------- * (in futuro possiamo pensare di interpolare anche sull'angolo) do iang=1,nangbin c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle if(angL(iang).lt.angle.and.angR(iang).ge.angle)then iangle=iang goto 98 endif enddo if(DEBUG.EQ.1) $ print*,'pfaeta4 *** 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 * ------------- iaddmax=10 iadd=0 10 continue if(cog4.lt.eta4(1,iang))then cog4 = cog4 + 1 iadd = iadd + 1 if(iadd>iaddmax)goto 111 goto 10 endif 20 continue if(cog4.gt.eta4(netaval,iang))then cog4 = cog4 - 1 iadd = iadd - 1 if(iadd<-1*iaddmax)goto 111 goto 20 endif goto 1111 111 continue if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster' if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)' cog4=0 1111 continue * -------------------------------- c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 do i=2,netaval if(eta4(i,iang).gt.cog4)then x1 = eta4(i-1,iang) x2 = eta4(i,iang) y1 = feta4(i-1,iview,lad,iang) y2 = feta4(i,iview,lad,iang) c print*,'*****',i,view,lad,iang c print*,'-----',x1,x2,y1,y2 goto 99 endif enddo 99 continue AA=(y2-y1)/(x2-x1) BB=y1-AA*x1 pfaeta4 = AA*cog4+BB pfaeta4 = pfaeta4 - iadd c$$$ if(iflag.eq.1)then c$$$ pfaeta2=pfaeta2-1. !temp c$$$ cog2=cog2-1. !temp c$$$ endif c$$$ if(iflag.eq.-1)then c$$$ pfaeta2=pfaeta2+1. !temp c$$$ cog2=cog2+1. !temp c$$$ endif if(DEBUG.EQ.1)print*,'ETA4 (ic ',ic,' ang',angle,')' $ ,cog4-iadd,' -->',pfaeta4 c 100 return return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function digsat(ic) *------------------------------------------------- * * *------------------------------------------------- include 'commontracker.f' include 'calib.f' include 'level1.f' integer nsat real pitchsat nsat = 0 pitchsat = 0. iv=VIEW(ic) istart = INDSTART(IC) istop = TOTCLLENGTH if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 do i = INDMAX(IC),istart,-1 if( (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx) $ .or. $ (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then nsat = nsat + 1 pitchsat = pitchsat + i - INDMAX(IC) else goto 10 endif enddo 10 continue do i = INDMAX(IC)+1,istop if( (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx) $ .or. $ (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then nsat = nsat + 1 pitchsat = pitchsat + i - INDMAX(IC) else goto 20 endif enddo 20 continue digsat = 0 if (nsat.gt.0) digsat = pitchsat / nsat return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function cog(ncog,ic) *------------------------------------------------- * this function returns * * - if NCOG=0, the Center-Of-Gravity of the * cluster IC, relative to MAXS(IC), according to * the cluster multiplicity * * - if NCOG>0, the Center-Of-Gravity of the cluster IC * evaluated using NCOG strips, even if they have a * negative signal (according to Landi) * *------------------------------------------------- include 'commontracker.f' include 'calib.f' include 'level1.f' if (ncog.gt.0) then * =========================== * ETA2 ETA3 ETA4 computation * =========================== * --> signal of the central strip sc = CLSIGNAL(INDMAX(ic)) !center * signal of adjacent strips sl1 = -9999. !left 1 if( $ (INDMAX(ic)-1).ge.INDSTART(ic) $ ) $ sl1 = CLSIGNAL(INDMAX(ic)-1) sl2 = -9999. !left 2 if( $ (INDMAX(ic)-2).ge.INDSTART(ic) $ ) $ sl2 = CLSIGNAL(INDMAX(ic)-2) sr1 = -9999. !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 = CLSIGNAL(INDMAX(ic)+1) sr2 = -9999. !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 = CLSIGNAL(INDMAX(ic)+2) COG = 0. c print *,'## ',sl2,sl1,sc,sr1,sr2 c ============================================================== if(ncog.eq.1)then COG = 0. 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 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1) $ .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) endif c if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic) c $ ,' : ',sl2,sl1,sc,sr1,sr2 c ============================================================== elseif(ncog.eq.3)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(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+sss).ne.0) $ COG = (COG-2*sl2)/(sl2+sss) elseif(sl2.lt.sr2)then 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+sss).ne.0 ) $ cog = (cog-2*sl2)/(sl2+sss) if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2) $ .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' COG = 0. endif c print*,'NCOG ',ncog,ic,' @@@ ',sl1,sc,sr1,' @@@ ',cog elseif(ncog.eq.0)then * ========================= * COG computation * ========================= 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 COG = 0 SGN = 0. mu = 0 c print*,'-------' do i = INDMAX(IC),istart,-1 ipos = i-INDMAX(ic) cut = incut*CLSIGMA(i) if(CLSIGNAL(i).ge.cut)then COG = COG + ipos*CLSIGNAL(i) SGN = SGN + CLSIGNAL(i) mu = mu + 1 c print*,ipos,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).ge.cut)then COG = COG + ipos*CLSIGNAL(i) SGN = SGN + CLSIGNAL(i) mu = mu + 1 c print*,ipos,CLSIGNAL(i) else goto 20 endif enddo 20 continue if(SGN.le.0)then print*,'cog(0,ic) --> ic, dedx ',ic,SGN print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) print*,(CLSIGNAL(i),i=istart,istop) c print*,'cog(0,ic) --> NOT EVALUATED ' else COG=COG/SGN endif c print*,'-------' else COG=0 print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG print*,' (NCOG must be >= 0)' endif c print *,'## cog ',ncog,ic,cog,'/////////////' if(COG.lt.-0.75.or.COG.gt.+0.75)then if(DEBUG.eq.1) $ print*,'cog *** warning *** anomalous cluster ??? --> ' if(DEBUG.eq.1) $ print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG endif return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function fbad_cog(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. * *------------------------------------------------------- include 'commontracker.f' include 'level1.f' include 'calib.f' if(mod(int(VIEW(ic)),2).eq.1)then !Y-view si = 8.4 !average good-strip noise f = 4. !average bad-strip noise: f*si incut=incuty else !X-view si = 3.9 !average good-strip noise f = 6. !average bad-strip noise: f*si incut=incutx endif fbad_cog = 1. if (ncog.gt.0) then * --> signal of the central strip sc = CLSIGNAL(INDMAX(ic)) !center fsc = 1 c if( CLBAD(INDMAX(ic)).eq.0 )fsc=f fsc = clsigma(INDMAX(ic))/si * --> 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) c if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f fsl1 = clsigma(INDMAX(ic)-1)/si endif sl2 = 0 !left 2 fsl2 = 1 !left 2 if( $ (INDMAX(ic)-2).ge.INDSTART(ic) $ )then sl2 = CLSIGNAL(INDMAX(ic)-2) c if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f fsl2 = clsigma(INDMAX(ic)-2)/si 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) c if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f fsr1 = clsigma(INDMAX(ic)+1)/si 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) c if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f fsr2 = clsigma(INDMAX(ic)+2)/si endif ************************************************************ * COG2-3-4 computation ************************************************************ c print*,sl2,sl1,sc,sr1,sr2 vCOG = cog(ncog,ic)!0. if(ncog.eq.2)then if(sl1.gt.sr1)then c COG = -sl1/(sl1+sc) fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2) fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2) elseif(sl1.le.sr1)then c COG = sr1/(sc+sr1) fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2) endif elseif(ncog.eq.3)then c COG = (sr1-sl1)/(sl1+sc+sr1) fbad_cog = $ (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) fbad_cog = $ fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2) elseif(ncog.eq.4)then if(sl2.gt.sr2)then c COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) fbad_cog = $ (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2 $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) fbad_cog = $ fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2 $ +(-vCOG)**2+(1-vCOG)**2) elseif(sl2.le.sr2)then c COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) fbad_cog = $ (fsl1*(-1-vCOG)**2 $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2) fbad_cog = $ fbad_cog / ((-1-vCOG)**2 $ +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2) endif else print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG print*,' (NCOG must be <= 4)' c COG = 0. 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 SGN = 0. SNU = 0. 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)/si SNU = SNU + fs*(ipos-vCOG)**2 SDE = SDE + (ipos-vCOG)**2 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)/si SNU = SNU + fs*(ipos-vCOG)**2 SDE = SDE + (ipos-vCOG)**2 else goto 20 endif enddo 20 continue if(SDE.ne.0)then FBAD_COG=SNU/SDE else endif else FBAD_COG=0 print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG print*,' (NCOG must be >= 0)' endif fbad_cog = sqrt(fbad_cog) return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** 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 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risxeta2(x) DOUBLE PRECISION V( 1) INTEGER NPAR, NDIM, IMQFUN, I, J DOUBLE PRECISION HQDJ, VV, VCONST DOUBLE PRECISION SIGVMI( 1), SIGVT( 1) DOUBLE PRECISION SIGV( 18, 1) DOUBLE PRECISION SIGDEL( 18) DOUBLE PRECISION SIGA( 18) DATA NPAR, NDIM, IMQFUN / 18, 1, 1/ DATA VCONST / 0.000000000000 / DATA SIGVMI / -20.50000000000 / DATA SIGVT / 41.00000000000 / DATA SIGV / 0.6097560748458E-01 +, 0.1097560971975 +, 0.1341463327408 +, 0.1829268187284 +, 0.2317073047161 +, 0.4268292486668 +, 0.4756097495556 +, 0.4999999701977 +, 0.5243902206421 +, 0.5731707215309 +, 0.7682926654816 +, 0.8170731663704 +, 0.8658536076546 +, 0.8902438879013 +, 0.9390243291855 +, 0.000000000000 +, 1.000000000000 +, 0.3658536374569 +/ DATA SIGDEL / 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.1999999994950E-05 +, 0.1999999994950E-05 +, 0.9756097197533E-01 +/ DATA SIGA / 51.65899502118 +, -150.4733247841 +, 143.0468613786 +, -16.56096738997 +, 5.149319798083 +, 21.57149712673 +, -39.46652322782 +, 47.13181632948 +, -32.93197883680 +, 16.38645317092 +, 1.453688482992 +, -10.00547244421 +, 131.3517670587 +, -140.6351538257 +, 49.05515749582 +, -23.00028974788 +, -22.58470403729 +, -3.824682486418 +/ V(1)= abs(x) if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR HQDJ = 0. DO 10 I = 1, NDIM VV = (V (I) - SIGVMI (I)) / SIGVT (I) HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2 10 CONTINUE HQDJ = HQDJ + SIGDEL (J) ** 2 HQDJ = SQRT (HQDJ) HQUADF = HQUADF + SIGA (J) * HQDJ 20 CONTINUE IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) risxeta2=HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risxeta3(x) DOUBLE PRECISION V( 1) INTEGER NPAR, NDIM, IMQFUN, I, J DOUBLE PRECISION HQDJ, VV, VCONST DOUBLE PRECISION SIGVMI( 1), SIGVT( 1) DOUBLE PRECISION SIGV( 18, 1) DOUBLE PRECISION SIGDEL( 18) DOUBLE PRECISION SIGA( 18) DATA NPAR, NDIM, IMQFUN / 18, 1, 1/ DATA VCONST / 0.000000000000 / DATA SIGVMI / -20.50000000000 / DATA SIGVT / 41.00000000000 / DATA SIGV / 0.6097560748458E-01 +, 0.1097560971975 +, 0.1341463327408 +, 0.1829268187284 +, 0.2317073047161 +, 0.4756097495556 +, 0.4999999701977 +, 0.5243902206421 +, 0.7682926654816 +, 0.8170731663704 +, 0.8658536076546 +, 0.8902438879013 +, 0.9390243291855 +, 0.000000000000 +, 1.000000000000 +, 0.3658536374569 +, 0.4146341383457 +, 0.6097560524940 +/ DATA SIGDEL / 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.1999999994950E-05 +, 0.1999999994950E-05 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +/ DATA SIGA / 55.18284054458 +, -160.3358431242 +, 144.6939185763 +, -20.45200854118 +, 5.223570087108 +,-0.4171476953945 +, -27.67911907462 +, 17.70327157495 +, -1.867165491707 +, -8.884458169181 +, 124.3526608791 +, -143.3309398345 +, 50.80345027122 +, -16.44454904415 +, -15.73785568450 +, -22.71810502561 +, 36.86170101430 +, 2.437918198452 +/ V(1) = abs(x) if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR HQDJ = 0. DO 10 I = 1, NDIM VV = (V (I) - SIGVMI (I)) / SIGVT (I) HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2 10 CONTINUE HQDJ = HQDJ + SIGDEL (J) ** 2 HQDJ = SQRT (HQDJ) HQUADF = HQUADF + SIGA (J) * HQDJ 20 CONTINUE IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) risxeta3 = HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risxeta4(x) DOUBLE PRECISION V( 1) INTEGER NPAR, NDIM, IMQFUN, I, J DOUBLE PRECISION HQDJ, VV, VCONST DOUBLE PRECISION SIGVMI( 1), SIGVT( 1) DOUBLE PRECISION SIGV( 18, 1) DOUBLE PRECISION SIGDEL( 18) DOUBLE PRECISION SIGA( 18) DATA NPAR, NDIM, IMQFUN / 18, 1, 1/ DATA VCONST / 0.000000000000 / DATA SIGVMI / -20.50000000000 / DATA SIGVT / 41.00000000000 / DATA SIGV / 0.3658536449075E-01 +, 0.6097560748458E-01 +, 0.1097560971975 +, 0.1341463327408 +, 0.4756097495556 +, 0.5243902206421 +, 0.8658536076546 +, 0.8902438879013 +, 0.9390243291855 +, 0.9634146094322 +, 0.000000000000 +, 1.000000000000 +, 0.3658536374569 +, 0.4146341383457 +, 0.6097560524940 +, 0.6585365533829 +, 0.7560975551605 +, 0.2439024299383 +/ DATA SIGDEL / 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.1999999994950E-05 +, 0.1999999994950E-05 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.1951219439507 +/ DATA SIGA / -43.61551887895 +, 57.88466995373 +, -92.04113299504 +, 74.08166649890 +, -9.768686062558 +, -4.304496875334 +, 72.62237333937 +, -91.21920840618 +, 56.75519978630 +, -43.21115751243 +, 12.79984505413 +, 12.10074868595 +, -6.238587250860 +, 23.43447356326 +, 17.98221401495 +, -7.980332610975 +, -3.426733307051 +, -8.683439558751 +/ V(1)=abs(x) if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR HQDJ = 0. DO 10 I = 1, NDIM VV = (V (I) - SIGVMI (I)) / SIGVT (I) HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2 10 CONTINUE HQDJ = HQDJ + SIGDEL (J) ** 2 HQDJ = SQRT (HQDJ) HQUADF = HQUADF + SIGA (J) * HQDJ 20 CONTINUE IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) risxeta4=HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risyeta2(x) DOUBLE PRECISION V( 1) INTEGER NPAR, NDIM, IMQFUN, I, J DOUBLE PRECISION HQDJ, VV, VCONST DOUBLE PRECISION SIGVMI( 1), SIGVT( 1) DOUBLE PRECISION SIGV( 12, 1) DOUBLE PRECISION SIGDEL( 12) DOUBLE PRECISION SIGA( 12) DATA NPAR, NDIM, IMQFUN / 12, 1, 1/ DATA VCONST / 0.000000000000 / DATA SIGVMI / -20.50000000000 / DATA SIGVT / 41.00000000000 / DATA SIGV / 0.1585365831852 +, 0.4024389982224 +, 0.4756097495556 +, 0.5243902206421 +, 0.5975609421730 +, 0.8414633870125 +, 0.000000000000 +, 1.000000000000 +, 0.2682926654816 +, 0.3170731663704 +, 0.7073170542717 +, 0.7560975551605 +/ DATA SIGDEL / 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.1999999994950E-05 +, 0.1999999994950E-05 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +/ DATA SIGA / 14.57433603529 +, -15.93532436156 +, -13.24628335221 +, -14.31193855410 +, -12.67339684488 +, 18.19876051780 +, -5.270493486725 +, -5.107670990828 +, -9.553262933901 +, 43.34150727448 +, 55.91366786432 +, -29.38037318563 +/ v(1)= abs(x) if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR HQDJ = 0. DO 10 I = 1, NDIM VV = (V (I) - SIGVMI (I)) / SIGVT (I) HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2 10 CONTINUE HQDJ = HQDJ + SIGDEL (J) ** 2 HQDJ = SQRT (HQDJ) HQUADF = HQUADF + SIGA (J) * HQDJ 20 CONTINUE IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) risyeta2=HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risy_cog(x) DOUBLE PRECISION V( 1) INTEGER NPAR, NDIM, IMQFUN, I, J DOUBLE PRECISION HQDJ, VV, VCONST DOUBLE PRECISION SIGVMI( 1), SIGVT( 1) DOUBLE PRECISION SIGV( 10, 1) DOUBLE PRECISION SIGDEL( 10) DOUBLE PRECISION SIGA( 10) DATA NPAR, NDIM, IMQFUN / 10, 1, 1/ DATA VCONST / 0.000000000000 / DATA SIGVMI / -20.50000000000 / DATA SIGVT / 41.00000000000 / DATA SIGV / 0.1585365831852 +, 0.8414633870125 +, 0.000000000000 +, 1.000000000000 +, 0.4634146094322 +, 0.5121951103210 +, 0.5609756112099 +, 0.6585365533829 +, 0.7073170542717 +, 0.3414633870125 +/ DATA SIGDEL / 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.1999999994950E-05 +, 0.1999999994950E-05 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.9756097197533E-01 +, 0.1951219439507 +/ DATA SIGA / 23.73833445988 +, 24.10182100013 +, 1.865894323190 +, 1.706006262931 +, -1.075607857202 +, -22.11489493403 +, 1.663100707801 +, 4.089852595440 +, -4.314993873697 +, -2.174479487744 +/ V(1)=abs(x) if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR HQDJ = 0. DO 10 I = 1, NDIM VV = (V (I) - SIGVMI (I)) / SIGVT (I) HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2 10 CONTINUE HQDJ = HQDJ + SIGDEL (J) ** 2 HQDJ = SQRT (HQDJ) HQUADF = HQUADF + SIGA (J) * HQDJ 20 CONTINUE IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) risy_cog=HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risx_cog(x) DOUBLE PRECISION V( 1) INTEGER NPAR, NDIM, IMQFUN, I, J DOUBLE PRECISION HQDJ, VV, VCONST DOUBLE PRECISION SIGVMI( 1), SIGVT( 1) DOUBLE PRECISION SIGV( 15, 1) DOUBLE PRECISION SIGDEL( 15) DOUBLE PRECISION SIGA( 15) DATA NPAR, NDIM, IMQFUN / 15, 1, 1/ DATA VCONST / 0.000000000000 / DATA SIGVMI / -20.50000000000 / DATA SIGVT / 41.00000000000 / DATA SIGV / 0.6097560748458E-01 +, 0.8536584675312E-01 +, 0.1341463327408 +, 0.2317073047161 +, 0.2804878056049 +, 0.3780487775803 +, 0.6219512224197 +, 0.7195121645927 +, 0.7682926654816 +, 0.8658536076546 +, 0.9146341085434 +, 0.9390243291855 +, 0.000000000000 +, 1.000000000000 +, 0.5121951103210 +/ DATA SIGDEL / 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.4878048598766E-01 +, 0.1999999994950E-05 +, 0.1999999994950E-05 +, 0.9756097197533E-01 +/ DATA SIGA / 31.95672945139 +, -34.23286209245 +, -6.298459168211 +, 10.98847700545 +,-0.3052213535054 +, 13.10517991464 +, 15.60290821679 +, -1.956118448507 +, 12.41453816720 +, -7.354056408553 +, -32.32512668778 +, 30.61116178966 +, 1.418505329236 +, 1.583492573619 +, -18.48799977042 +/ V(1)=abs(x) if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR HQDJ = 0. DO 10 I = 1, NDIM VV = (V (I) - SIGVMI (I)) / SIGVT (I) HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2 10 CONTINUE HQDJ = HQDJ + SIGDEL (J) ** 2 HQDJ = SQRT (HQDJ) HQUADF = HQUADF + SIGA (J) * HQDJ 20 CONTINUE IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) risx_cog = HQUADF * 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function pfacorr(ic,angle) *-------------------------------------------------------------- * this function returns the landi correction for this cluster *-------------------------------------------------------------- include 'commontracker.f' include 'calib.f' include 'level1.f' real angle integer iview,lad iview = VIEW(ic) lad = nld(MAXS(ic),VIEW(ic)) * 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*,'pfacorr *** 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 pfacorr = fcorr(iview,lad,iang) if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr c 100 return return end