--- DarthVader/TrackerLevel2/src/F77/functionspfa.f 2007/02/16 14:56:02 1.8 +++ DarthVader/TrackerLevel2/src/F77/functionspfa.f 2007/04/27 10:39:58 1.9 @@ -1,3 +1,4 @@ + *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** * this file contains all subroutines and functions * that are needed for position finding algorithms @@ -6,6 +7,98 @@ *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** + 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) *-------------------------------------------------------------- @@ -17,29 +110,37 @@ * according to the angle *-------------------------------------------------------------- include 'commontracker.f' -c include 'calib.f' include 'level1.f' + include 'calib.f' pfaeta = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view - pfaeta = pfaeta2(ic,angle) - + 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).le.10.)then + if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then pfaeta = pfaeta2(ic,angle) - elseif(abs(angle).gt.10..and.abs(angle).le.15.)then + elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then pfaeta = pfaeta3(ic,angle) - elseif(abs(angle).gt.15.)then + elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then pfaeta = pfaeta4(ic,angle) - endif + else + pfaeta = cog(4,ic) + endif endif -c print*,'pfaeta ',pfaeta, angle - 100 return end @@ -56,37 +157,38 @@ * according to the angle *-------------------------------------------------------------- include 'commontracker.f' -c include 'calib.f' include 'level1.f' + include 'calib.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.) + + if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then + ris_eta = risy_eta2(angle) + elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then + ris_eta = risy_cog(angle) !ATTENZIONE!! + elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then + ris_eta = risy_cog(angle) !ATTENZIONE!! + else + ris_eta = risy_cog(angle) + endif else !X-view - if(abs(angle).le.10.)then + if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )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 + elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then + ris_eta = risx_eta3(angle) + elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then ris_eta = risx_eta4(angle) - elseif(abs(angle).gt.21.)then - ris_eta = risx_eta4(21.) - endif + else + ris_eta = risx_cog(angle) + endif endif -c$$$ if(DEBUG)print*,'ris (ic ',ic,' ang',angle,')' -c$$$ $ ,' -->',ris_eta - - 100 return end @@ -104,22 +206,32 @@ include 'commontracker.f' include 'level1.f' -* include 'calib.f' + include 'calib.f' fbad_eta = 0 if(mod(int(VIEW(ic)),2).eq.1)then !Y-view - fbad_eta = fbad_cog(2,ic) + 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).le.10.)then + if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then fbad_eta = fbad_cog(2,ic) - elseif(abs(angle).gt.10..and.abs(angle).le.15.)then + elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then fbad_eta = fbad_cog(3,ic) - elseif(abs(angle).gt.15.)then + elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then fbad_eta = fbad_cog(4,ic) - endif + else + fbad_eta = fbad_cog(4,ic) + endif endif @@ -127,10 +239,6 @@ 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 @@ -150,19 +258,14 @@ 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) + 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 -c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle if(angL(iang).lt.angle.and.angR(iang).ge.angle)then iangle=iang goto 98 @@ -252,10 +355,6 @@ 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 @@ -275,14 +374,10 @@ 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) + iview = VIEW(ic) + lad = nld(MAXS(ic),VIEW(ic)) + cog3 = cog(3,ic) pfaeta3=cog3 * find angular bin @@ -376,11 +471,7 @@ 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) + real function pfaeta4(ic,angle) *-------------------------------------------------------------- * this function returns * @@ -399,14 +490,10 @@ 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) + iview = VIEW(ic) + lad = nld(MAXS(ic),VIEW(ic)) + cog4=cog(4,ic) pfaeta4=cog4 * find angular bin @@ -618,8 +705,6 @@ include 'calib.f' include 'level1.f' -c logical DEBUG -c common/dbg/DEBUG if (ncog.gt.0) then @@ -696,28 +781,48 @@ 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 - do i = istart,istop +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 + 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 + print*,ipos,CLSIGNAL(i) + else + goto 20 endif enddo - if(SGNL(ic).le.0)then - print*,'cog(0,ic) --> ic, dedx ',ic,SGNL(ic) + 20 continue + if(SGN.le.0)then +c print*,'cog(0,ic) --> ic, dedx ',ic,SGN print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) print*,(CLSIGNAL(i),i=istart,istop) - print*,'cog(0,ic) --> NOT EVALUATED ' +c print*,'cog(0,ic) --> NOT EVALUATED ' else - COG=COG/SGNL(ic) + COG=COG/SGN endif +c print*,'-------' else @@ -734,6 +839,7 @@ end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** + real function fbad_cog(ncog,ic) *------------------------------------------------------- * this function returns a factor that takes into @@ -749,12 +855,12 @@ include 'calib.f' if(mod(int(VIEW(ic)),2).eq.1)then !Y-view - f = 4. - si = 8.4 + si = 8.4 !average good-strip noise + f = 4. !average bad-strip noise: f*si incut=incuty else !X-view - f = 6. - si = 3.9 + si = 3.9 !average good-strip noise + f = 6. !average bad-strip noise: f*si incut=incutx endif @@ -765,8 +871,8 @@ * --> signal of the central strip sc = CLSIGNAL(INDMAX(ic)) !center fsc = 1 -c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)fsc=f - if( CLBAD(INDMAX(ic)).eq.0 )fsc=f +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 @@ -774,10 +880,8 @@ $ (INDMAX(ic)-1).ge.INDSTART(ic) $ )then sl1 = CLSIGNAL(INDMAX(ic)-1) -c if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f - if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f -c else -c fsl1 = 0 +c if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f + fsl1 = clsigma(INDMAX(ic)-1)/si endif sl2 = 0 !left 2 @@ -786,10 +890,8 @@ $ (INDMAX(ic)-2).ge.INDSTART(ic) $ )then sl2 = CLSIGNAL(INDMAX(ic)-2) -c if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f - if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f -c else -c fsl2 = 0 +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 @@ -799,10 +901,8 @@ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) $ )then sr1 = CLSIGNAL(INDMAX(ic)+1) -c if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f - if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f -c else -c fsr1 = 0 +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 @@ -812,95 +912,134 @@ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) $ )then sr2 = CLSIGNAL(INDMAX(ic)+2) -c if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f - if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f -c else -c fsr2 = 0 +c if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f + fsr2 = clsigma(INDMAX(ic)+2)/si endif ************************************************************ -* COG computation +* COG2-3-4 computation ************************************************************ c print*,sl2,sl1,sc,sr1,sr2 - COG = 0. + vCOG = cog(ncog,ic)!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) +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 - COG = sr1/(sc+sr1) - fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2) - fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2) +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 - COG = (sr1-sl1)/(sl1+sc+sr1) +c COG = (sr1-sl1)/(sl1+sc+sr1) fbad_cog = - $ (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2) + $ (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) fbad_cog = - $ fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2) + $ fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2) elseif(ncog.eq.4)then if(sl2.gt.sr2)then - COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) +c 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) + $ (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2 + $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) fbad_cog = - $ fbad_cog / ((-2-COG)**2+(-1-COG)**2 - $ +(-COG)**2+(1-COG)**2) + $ fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2 + $ +(-vCOG)**2+(1-vCOG)**2) elseif(sl2.le.sr2)then - COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) +c 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) + $ (fsl1*(-1-vCOG)**2 + $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2) fbad_cog = - $ fbad_cog / ((-1-COG)**2 - $ +(-COG)**2+(1-COG)**2+(2-COG)**2) + $ 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)' - COG = 0. +c COG = 0. endif elseif(ncog.eq.0)then - +* ========================= +* COG computation +* ========================= - 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) + 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 - COG = COG + ipos*CLSIGNAL(i) - endif + fs = clsigma(i)/si + SNU = SNU + fs*(ipos-vCOG)**2 + SDE = SDE + (ipos-vCOG)**2 + else + goto 10 + endif enddo - COG=COG/SGNL(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) + 10 continue + do i=INDMAX(IC)+1,istop + ipos = i-INDMAX(ic) + cut = incut*CLSIGMA(i) 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 + fs = clsigma(i)/si + SNU = SNU + fs*(ipos-vCOG)**2 + SDE = SDE + (ipos-vCOG)**2 + else + goto 20 endif enddo - if(SDE.ne.0)FBAD_COG=SNU/SDE + 20 continue + if(SDE.ne.0)then + FBAD_COG=SNU/SDE + else + + endif else @@ -1103,6 +1242,7 @@ +/ V(1)= abs(x) + if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR @@ -1193,6 +1333,8 @@ +/ V(1) = abs(x) + if(V(1).gt.20.)V(1)=20. + HQUADF = 0. DO 20 J = 1, NPAR HQDJ = 0. @@ -1281,6 +1423,7 @@ +/ V(1)=abs(x) + if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR @@ -1352,6 +1495,7 @@ +/ v(1)= abs(x) + if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR @@ -1418,6 +1562,7 @@ +/ V(1)=abs(x) + if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR @@ -1498,6 +1643,7 @@ +/ V(1)=abs(x) + if(V(1).gt.20.)V(1)=20. HQUADF = 0. DO 20 J = 1, NPAR @@ -1515,4 +1661,3 @@ risx_cog = HQUADF * 1e-4 END -*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***