*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
* 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)
*
* 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)
*
* 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
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
* 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
100 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
100 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
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
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) !(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
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.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*,'pfaeta2 *** 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
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
100 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)
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.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
pfaeta3 = AA*cog3+BB
pfaeta3 = pfaeta3 - 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*,'ETA3 (ic ',ic,' ang',angle,')'
$ ,cog3-iadd,' -->',pfaeta3
100 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.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
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
100 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,'/////////////'
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
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$$$
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
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.lt.angL(1))iang=1
if(angle.gt.angR(nangbin))iang=nangbin
98 continue !jump here if ok
pfacorr = fcorr(iview,lad,iang)
if(DEBUG.eq.1)print*,'CORR (ic ',ic,' ang',angle,') -->',pfacorr
100 return
end