*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** * this file contains all subroutines and functions * that are needed for position finding algorithms * * *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function pfa_eta(ic,angle) *-------------------------------------------------------------- * this function returns the position (in strip units) * it calls: * - pfa_eta2(ic,angle) * - pfa_eta3(ic,angle) * - pfa_eta4(ic,angle) * according to the angle *-------------------------------------------------------------- include 'commontracker.f' c include 'calib.f' include 'level1.f' pfa_eta = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view pfa_eta = pfa_eta2(ic,angle) else !X-view if(abs(angle).le.10.)then pfa_eta = pfa_eta2(ic,angle) elseif(abs(angle).gt.10..and.abs(angle).le.15.)then pfa_eta = pfa_eta3(ic,angle) elseif(abs(angle).gt.15.)then pfa_eta = pfa_eta4(ic,angle) endif endif 100 return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function ris_eta(ic,angle) *-------------------------------------------------------------- * this function returns the average spatial resolution * (in cm) for the ETA algorithm (function pfa_eta(ic,angle)) * it calls: * - risx_eta2(angle) * - risy_eta2(angle) * - risx_eta3(angle) * - risx_eta4(angle) * according to the angle *-------------------------------------------------------------- include 'commontracker.f' c include 'calib.f' include 'level1.f' c$$$ logical DEBUG c$$$ common/dbg/DEBUG ris_eta = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view ris_eta = risy_eta2(angle) if(abs(angle).gt.21.)ris_eta = risy_eta2(21.) else !X-view if(abs(angle).le.10.)then ris_eta = risx_eta2(angle) elseif(abs(angle).gt.10..and.abs(angle).le.15.)then ris_eta = risx_eta3(angle) elseif(abs(angle).gt.15..and.abs(angle).le.21.)then ris_eta = risx_eta4(angle) elseif(abs(angle).gt.21.)then ris_eta = risx_eta4(21.) endif endif c$$$ if(DEBUG)print*,'ris (ic ',ic,' ang',angle,')' c$$$ $ ,' -->',ris_eta 100 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 *------------------------------------------------------- include 'commontracker.f' include 'level1.f' * include 'calib.f' fbad_eta = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view fbad_eta = fbad_cog(2,ic) else !X-view if(abs(angle).le.10.)then fbad_eta = fbad_cog(2,ic) elseif(abs(angle).gt.10..and.abs(angle).le.15.)then fbad_eta = fbad_cog(3,ic) elseif(abs(angle).gt.15.)then fbad_eta = fbad_cog(4,ic) endif endif return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** c***************************************************** cccccc 02/02/2006 modified by Elena Vannuccini --> (1) c***************************************************** c real function pfa_eta2(cog2,view,lad,angle) real function pfa_eta2(ic,angle) !(1) *-------------------------------------------------------------- * 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 c logical DEBUG c common/dbg/DEBUG iview = VIEW(ic) !(1) lad = nld(MAXS(ic),VIEW(ic)) !(1) cog2 = cog(2,ic) !(1) pfa_eta2=cog2 * 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) $ print*,'pfa_eta2 *** warning *** angle out of range: ',angle if(angle.lt.angL(1))iang=1 if(angle.gt.angR(nangbin))iang=nangbin 98 continue !jump here if ok c$$$* find extremes of interpolation c$$$ iflag=0 c$$$* -------------------------------- c$$$ if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then c$$$c print*,'pfa_eta2 *** warning *** argument out of range: ',cog2 c$$$* goto 100 c$$$* ---------------------------------------------- c$$$* non salto piu`, ma scalo di 1 o -1 c$$$* nel caso si tratti di un cluster c$$$* in cui la strip con il segnale massimo non coincide c$$$* con la strip con il rapposto s/n massimo!!! c$$$* ---------------------------------------------- c$$$ if(cog2.lt.eta2(1,iang))then !temp c$$$ cog2=cog2+1. !temp c$$$ iflag=1 !temp c$$$ else !temp c$$$ cog2=cog2-1. !temp c$$$ iflag=-1 !temp c$$$ endif !temp c$$$c print*,'shifted >>> ',cog2 c$$$ endif iadd=0 10 continue if(cog2.lt.eta2(1,iang))then cog2 = cog2 + 1 iadd = iadd + 1 goto 10 endif 20 continue if(cog2.gt.eta2(netaval,iang))then cog2 = cog2 - 1 iadd = iadd - 1 goto 20 endif * -------------------------------- 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 pfa_eta2 = AA*cog2+BB pfa_eta2 = pfa_eta2 - iadd c$$$ if(iflag.eq.1)then c$$$ pfa_eta2=pfa_eta2-1. !temp c$$$ cog2=cog2-1. !temp c$$$ endif c$$$ if(iflag.eq.-1)then c$$$ pfa_eta2=pfa_eta2+1. !temp c$$$ cog2=cog2+1. !temp c$$$ endif if(DEBUG)print*,'ETA2 (ic ',ic,' ang',angle,')' $ ,cog2-iadd,' -->',pfa_eta2 100 return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** c***************************************************** cccccc 02/02/2006 modified by Elena Vannuccini --> (1) c***************************************************** c real function pfa_eta3(cog3,view,lad,angle) real function pfa_eta3(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 c logical DEBUG c common/dbg/DEBUG iview = VIEW(ic) !(1) lad = nld(MAXS(ic),VIEW(ic)) !(1) cog3 = cog(3,ic) !(1) pfa_eta3=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) $ print*,'pfa_eta3 *** warning *** angle out of range: ',angle if(angle.lt.angL(1))iang=1 if(angle.gt.angR(nangbin))iang=nangbin 98 continue !jump here if ok c$$$* find extremes of interpolation c$$$ iflag=0 c$$$* -------------------------------- c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then c$$$* ---------------------------------------------- c$$$* non salto piu`, ma scalo di 1 o -1 c$$$* nel caso si tratti di un cluster c$$$* in cui la strip con il segnale massimo non coincide c$$$* con la strip con il rapposto s/n massimo!!! c$$$* ---------------------------------------------- c$$$ if(cog2.lt.eta2(1,iang))then !temp c$$$ cog2=cog2+1. !temp c$$$ iflag=1 !temp c$$$ else !temp c$$$ cog2=cog2-1. !temp c$$$ iflag=-1 !temp c$$$ endif !temp c$$$c print*,'shifted >>> ',cog2 c$$$ endif iadd=0 10 continue if(cog3.lt.eta3(1,iang))then cog3 = cog3 + 1 iadd = iadd + 1 goto 10 endif 20 continue if(cog3.gt.eta3(netaval,iang))then cog3 = cog3 - 1 iadd = iadd - 1 goto 20 endif * -------------------------------- 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 pfa_eta3 = AA*cog3+BB pfa_eta3 = pfa_eta3 - iadd c$$$ if(iflag.eq.1)then c$$$ pfa_eta2=pfa_eta2-1. !temp c$$$ cog2=cog2-1. !temp c$$$ endif c$$$ if(iflag.eq.-1)then c$$$ pfa_eta2=pfa_eta2+1. !temp c$$$ cog2=cog2+1. !temp c$$$ endif if(DEBUG)print*,'ETA3 (ic ',ic,' ang',angle,')' $ ,cog3-iadd,' -->',pfa_eta3 100 return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** c***************************************************** cccccc 02/02/2006 modified by Elena Vannuccini --> (1) c***************************************************** c real function pfa_eta4(cog4,view,lad,angle) real function pfa_eta4(ic,angle) !(1) *-------------------------------------------------------------- * 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 c logical DEBUG c common/dbg/DEBUG iview = VIEW(ic) !(1) lad = nld(MAXS(ic),VIEW(ic)) !(1) cog4=cog(4,ic) !(1) pfa_eta4=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) $ print*,'pfa_eta4 *** warning *** angle out of range: ',angle if(angle.lt.angL(1))iang=1 if(angle.gt.angR(nangbin))iang=nangbin 98 continue !jump here if ok c$$$* find extremes of interpolation c$$$ iflag=0 c$$$* -------------------------------- c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then c$$$* ---------------------------------------------- c$$$* non salto piu`, ma scalo di 1 o -1 c$$$* nel caso si tratti di un cluster c$$$* in cui la strip con il segnale massimo non coincide c$$$* con la strip con il rapposto s/n massimo!!! c$$$* ---------------------------------------------- c$$$ if(cog2.lt.eta2(1,iang))then !temp c$$$ cog2=cog2+1. !temp c$$$ iflag=1 !temp c$$$ else !temp c$$$ cog2=cog2-1. !temp c$$$ iflag=-1 !temp c$$$ endif !temp c$$$c print*,'shifted >>> ',cog2 c$$$ endif iadd=0 10 continue if(cog4.lt.eta4(1,iang))then cog4 = cog4 + 1 iadd = iadd + 1 goto 10 endif 20 continue if(cog4.gt.eta4(netaval,iang))then cog4 = cog4 - 1 iadd = iadd - 1 goto 20 endif * -------------------------------- 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 pfa_eta4 = AA*cog4+BB pfa_eta4 = pfa_eta4 - iadd c$$$ if(iflag.eq.1)then c$$$ pfa_eta2=pfa_eta2-1. !temp c$$$ cog2=cog2-1. !temp c$$$ endif c$$$ if(iflag.eq.-1)then c$$$ pfa_eta2=pfa_eta2+1. !temp c$$$ cog2=cog2+1. !temp c$$$ endif if(DEBUG)print*,'ETA4 (ic ',ic,' ang',angle,')' $ ,cog4-iadd,' -->',pfa_eta4 100 return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function cog0(ncog,ic) *------------------------------------------------- * this function returns * * - the Center-Of-Gravity of the cluster IC * evaluated using NCOG strips, * calculated relative to MAXS(IC) * * - zero in case that not enough strips * have a positive signal * * NOTE: * This is the old definition, used by Straulino. * The new routine, according to Landi, * is COG(NCOG,IC) *------------------------------------------------- include 'commontracker.f' include 'level1.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)) ************************************************************ * COG computation ************************************************************ c print*,sl2,sl1,sc,sr1,sr2 COG = 0. if(sl1.gt.sr1.and.sl1.gt.0.)then if(ncog.eq.2.and.sl1.ne.0)then COG = -sl1/(sl1+sc) elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then COG = (sr1-sl1)/(sl1+sc+sr1) elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) else COG = 0. endif elseif(sl1.le.sr1.and.sr1.gt.0.)then if(ncog.eq.2.and.sr1.ne.0)then COG = sr1/(sc+sr1) elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then COG = (sr1-sl1)/(sl1+sc+sr1) elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) else COG = 0. endif endif COG0 = COG c print *,ncog,ic,cog,'/////////////' 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' c logical DEBUG c common/dbg/DEBUG if (ncog.gt.0) then * =========================== * ETA2 ETA3 ETA4 computation * =========================== * --> signal of the central strip sc = CLSIGNAL(INDMAX(ic)) !center * signal of adjacent strips sl1 = 0 !left 1 if( $ (INDMAX(ic)-1).ge.INDSTART(ic) $ ) $ sl1 = CLSIGNAL(INDMAX(ic)-1) sl2 = 0 !left 2 if( $ (INDMAX(ic)-2).ge.INDSTART(ic) $ ) $ sl2 = CLSIGNAL(INDMAX(ic)-2) 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 = 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 = CLSIGNAL(INDMAX(ic)+2) COG = 0. if(ncog.eq.2)then if(sl1.gt.sr1)then COG = -sl1/(sl1+sc) elseif(sl1.le.sr1)then COG = sr1/(sc+sr1) endif elseif(ncog.eq.3)then COG = (sr1-sl1)/(sl1+sc+sr1) elseif(ncog.eq.4)then if(sl2.gt.sr2)then COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) elseif(sl2.le.sr2)then COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) endif else print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG print*,' (NCOG must be <= 4)' 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 mu = 0 do i = istart,istop ipos = i-INDMAX(ic) ivk = nvk(MAXS(ic)+ipos) is = nst(MAXS(ic)+ipos) * print*,'******************',istart,istop,ipos * $ ,MAXS(ic)+ipos,iv,ivk,is cut = incut*SIGMA(iv,ivk,is) if(SIGMA(iv,ivk,is).ne.CLSIGMA(i)) $ print*,'cog(0,ic) --> hai fatto qualche cazzata' if(CLSIGNAL(i).ge.cut)then COG = COG + ipos*CLSIGNAL(i) mu = mu + 1 c print*,ipos,CLSIGNAL(i),incut,cut endif enddo if(DEDX(ic).le.0)then print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic) print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) print*,(CLSIGNAL(i),i=istart,istop) print*,'cog(0,ic) --> NOT EVALUATED ' else COG=COG/DEDX(ic) endif c if(DEBUG)print*,'COG (ic ',ic,' m',mu,')' c $ ,cog else COG=0 print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG print*,' (NCOG must be >= 0)' endif c print *,ncog,ic,cog,'/////////////' 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 f = 4. si = 8.4 incut=incuty else !X-view f = 6. si = 3.9 incut=incutx endif fbad_cog = 1. if (ncog.gt.0) then * --> signal of the central strip sc = CLSIGNAL(INDMAX(ic)) !center fsc = 1 if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)fsc=f * --> 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) if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f c else c fsl1 = 0 endif sl2 = 0 !left 2 fsl2 = 1 !left 2 if( $ (INDMAX(ic)-2).ge.INDSTART(ic) $ )then sl2 = CLSIGNAL(INDMAX(ic)-2) if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f c else c fsl2 = 0 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) if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f c else c fsr1 = 0 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) if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f c else c fsr2 = 0 endif ************************************************************ * COG computation ************************************************************ c print*,sl2,sl1,sc,sr1,sr2 COG = 0. if(ncog.eq.2)then if(sl1.gt.sr1)then COG = -sl1/(sl1+sc) fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2) fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2) elseif(sl1.le.sr1)then COG = sr1/(sc+sr1) fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2) fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2) endif elseif(ncog.eq.3)then COG = (sr1-sl1)/(sl1+sc+sr1) fbad_cog = $ (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2) fbad_cog = $ fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2) elseif(ncog.eq.4)then if(sl2.gt.sr2)then COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) fbad_cog = $ (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2 $ +fsc*(-COG)**2+fsr1*(1-COG)**2) fbad_cog = $ fbad_cog / ((-2-COG)**2+(-1-COG)**2 $ +(-COG)**2+(1-COG)**2) elseif(sl2.le.sr2)then COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) fbad_cog = $ (fsl1*(-1-COG)**2 $ +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2) fbad_cog = $ fbad_cog / ((-1-COG)**2 $ +(-COG)**2+(1-COG)**2+(2-COG)**2) endif else print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG print*,' (NCOG must be <= 4)' COG = 0. endif elseif(ncog.eq.0)then iv=VIEW(ic) istart=INDSTART(IC) istop=TOTCLLENGTH if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 COG=0. SNU=0. SDE=0. do i=istart,istop ipos=i-INDMAX(ic) il=nvk(MAXS(ic)+ipos) is=nst(MAXS(ic)+ipos) cut=incut*SIGMA(iv,il,is) if(CLSIGNAL(i).gt.cut)then COG = COG + ipos*CLSIGNAL(i) endif enddo COG=COG/DEDX(ic) do i=istart,istop ipos=i-INDMAX(ic) il=nvk(MAXS(ic)+ipos) is=nst(MAXS(ic)+ipos) cut=incut*SIGMA(iv,il,is) if(CLSIGNAL(i).gt.cut)then fs=1 if(BAD(iv,il,is).eq.0)fs=f SNU = SNU + fs*(ipos-COG)**2 SDE = SDE + (ipos-COG)**2 endif enddo if(SDE.ne.0)FBAD_COG=SNU/SDE 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 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 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risx_eta2(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) 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_eta2=HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risx_eta3(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) 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_eta3 = HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risx_eta4(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) 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_eta4=HQUADF* 1e-4 END *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** FUNCTION risy_eta2(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) 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_eta2=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) 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) 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 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***