*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
* this file contains all subroutines and functions
* that are needed for position finding algorithms
*
*
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
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'
c include 'calib.f'
include 'level1.f'
pfaeta = 0
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
pfaeta = pfaeta2(ic,angle)
else !X-view
if(abs(angle).le.10.)then
pfaeta = pfaeta2(ic,angle)
elseif(abs(angle).gt.10..and.abs(angle).le.15.)then
pfaeta = pfaeta3(ic,angle)
elseif(abs(angle).gt.15.)then
pfaeta = pfaeta4(ic,angle)
endif
endif
c print*,'pfaeta ',pfaeta, angle
100 return
end
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
real function ris_eta(ic,angle)
*--------------------------------------------------------------
* this function returns the average spatial resolution
* (in cm) for the ETA algorithm (function pfaeta(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 pfaeta2(cog2,view,lad,angle)
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
c logical DEBUG
c common/dbg/DEBUG
c print*,'## pfaeta2 ',ic,angle
iview = VIEW(ic) !(1)
lad = nld(MAXS(ic),VIEW(ic)) !(1)
cog2 = cog(2,ic) !(1)
pfaeta2=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*,'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)print*,'ETA2 (ic ',ic,' ang',angle,')'
$ ,cog2-iadd,' -->',pfaeta2
100 return
end
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c*****************************************************
cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
c*****************************************************
c real function pfaeta3(cog3,view,lad,angle)
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
c logical DEBUG
c common/dbg/DEBUG
c print*,'## pfaeta3 ',ic,angle
iview = VIEW(ic) !(1)
lad = nld(MAXS(ic),VIEW(ic)) !(1)
cog3 = cog(3,ic) !(1)
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)
$ 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)print*,'ETA3 (ic ',ic,' ang',angle,')'
$ ,cog3-iadd,' -->',pfaeta3
100 return
end
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c*****************************************************
cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
c*****************************************************
c real function pfaeta4(cog4,view,lad,angle)
real function pfaeta4(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
c print*,'## pfaeta4 ',ic,angle
iview = VIEW(ic) !(1)
lad = nld(MAXS(ic),VIEW(ic)) !(1)
cog4=cog(4,ic) !(1)
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)
$ 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)print*,'ETA4 (ic ',ic,' ang',angle,')'
$ ,cog4-iadd,' -->',pfaeta4
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.
c print*,'## ',sl2,sl1,sc,sr1,sr2
if(ncog.eq.1)then
COG = 0.
elseif(ncog.eq.2)then
if(sl1.gt.sr1)then
if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)
elseif(sl1.le.sr1)then
if((sc+sr1).ne.0)COG = sr1/(sc+sr1)
endif
elseif(ncog.eq.3)then
if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)
elseif(ncog.eq.4)then
if(sl2.gt.sr2)then
if((sl2+sl1+sc+sr1).ne.0)
$ COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
elseif(sl2.le.sr2)then
if((sl2+sl1+sc+sr1).ne.0)
$ COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
endif
else
print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
print*,' (NCOG must be 0-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)
cc ivk = nvk(MAXS(ic)+ipos)
cc is = nst(MAXS(ic)+ipos)
* print*,'******************',istart,istop,ipos
* $ ,MAXS(ic)+ipos,iv,ivk,is
cc cut = incut*SIGMA(iv,ivk,is)
cut = incut*CLSIGMA(i)
c if(SIGMA(iv,ivk,is).ne.CLSIGMA(i))
c $ 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 *,'## 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
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
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***