subroutine idtoc(ipfa,cpfa) integer ipfa character*4 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.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 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** * this file contains all subroutines and functions * that are needed for position finding algorithms * * *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** integer function npfastrips(ic,PFA,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,PFA usedPFA=PFA npfastrips=0 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')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 c usedPFA='COG' 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 c usedPFA='COG' endif endif endif * ---------------------------------------------------------------- if(usedPFA.eq.'COG')then 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 mu = 0 do i = INDMAX(IC),istart,-1 ipos = i-INDMAX(ic) cut = incut*CLSIGMA(i) if(CLSIGNAL(i).ge.cut)then mu = mu + 1 print*,i,mu 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 mu = mu + 1 print*,i,mu else goto 20 endif enddo 20 continue npfastrips=mu endif * ---------------------------------------------------------------- c print*,pfastrips 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).le.e2tay )then pfaeta = pfaeta2(ic,angle) elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then pfaeta = pfaeta3(ic,angle) elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then pfaeta = pfaeta4(ic,angle) else pfaeta = cog(4,ic) endif else !X-view if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then pfaeta = pfaeta2(ic,angle) elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then pfaeta = pfaeta3(ic,angle) elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then pfaeta = pfaeta4(ic,angle) else pfaeta = 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 print*,'---- ',riseta,iview,angle 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) $ 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 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** 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) $ 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 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** 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) $ 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' 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((sr2+sl1+sc+sr1).ne.0) $ COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1) endif 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. c$$$ do i=INDMAX(IC),istart,-1 c$$$ ipos = i-INDMAX(ic) c$$$ cut = incut*CLSIGMA(i) c$$$ if(CLSIGNAL(i).gt.cut)then c$$$ COG = COG + ipos*CLSIGNAL(i) c$$$ SGN = SGN + CLSIGNAL(i) 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).gt.cut)then c$$$ COG = COG + ipos*CLSIGNAL(i) c$$$ SGN = SGN + CLSIGNAL(i) c$$$ else c$$$ goto 20 c$$$ endif c$$$ enddo c$$$ 20 continue c$$$ if(SGN.le.0)then c$$$ print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN c$$$ print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) c$$$ print*,(CLSIGNAL(i),i=istart,istop) c$$$ print*,'fbad_cog(0,ic) --> NOT EVALUATED ' c$$$ else c$$$ COG=COG/SGN c$$$ endif do i=INDMAX(IC),istart,-1 ipos = i-INDMAX(ic) cut = incut*CLSIGMA(i) if(CLSIGNAL(i).gt.cut)then fs = clsigma(i)/si SNU = SNU + fs*(ipos-vCOG)**2 SDE = SDE + (ipos-vCOG)**2 else goto 10 endif enddo 10 continue do i=INDMAX(IC)+1,istop ipos = i-INDMAX(ic) cut = incut*CLSIGMA(i) if(CLSIGNAL(i).gt.cut)then fs = clsigma(i)/si SNU = SNU + fs*(ipos-vCOG)**2 SDE = SDE + (ipos-vCOG)**2 else goto 20 endif enddo 20 continue if(SDE.ne.0)then FBAD_COG=SNU/SDE else endif else FBAD_COG=0 print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG print*,' (NCOG must be >= 0)' endif fbad_cog = sqrt(fbad_cog) return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** real function 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 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