--- DarthVader/TrackerLevel2/src/F77/mini.f 2006/10/26 16:22:37 1.3 +++ DarthVader/TrackerLevel2/src/F77/mini.f 2006/11/07 15:55:11 1.4 @@ -22,7 +22,7 @@ c logical DEBUG c common/dbg/DEBUG - parameter (inf=1.e8) !just a huge number... + parameter (dinf=1.d15) !just a huge number... c------------------------------------------------------------------------ c variables used in the tracking procedure (mini and its subroutines) c @@ -36,20 +36,21 @@ c DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking DATA STEPAL/5*1.d-7/ !alpha vector step - DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization + DATA ISTEPMAX/120/ !maximum number of steps in the chi^2 minimization DATA TOLL/1.d-8/ !tolerance in reaching the next plane during * !the tracking procedure DATA STEPMAX/100./ !maximum number of steps in the trackin gprocess -c DATA ALMAX/inf,inf,inf,inf,0.25e2/ !limits on alpha vector components -c DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !" - DATA ALMAX/inf,inf,1.,inf,0.25e2/ !limits on alpha vector components - DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !" - + DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components + DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" DIMENSION DAL(5) !increment of vector alfa DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2 - + +c elena-------- + REAL*8 AVRESX,AVRESY +c elena-------- + INTEGER IFLAG c-------------------------------------------------------- c IFLAG =1 ---- chi2 derivatives computed by using @@ -62,33 +63,68 @@ c-------------------------------------------------------- DATA IFLAG/2/ - LOGICAL TRKDEBUG - COMMON/TRKD/TRKDEBUG +c LOGICAL TRKDEBUG,TRKVERBOSE +c COMMON/TRKD/TRKDEBUG,TRKVERBOSE + LOGICAL TRKDEBUG,TRKVERBOSE + COMMON/TRKD/TRKDEBUG,TRKVERBOSE IF(IPRINT.EQ.1) THEN - TRKDEBUG = .TRUE. + TRKVERBOSE = .TRUE. + TRKDEBUG = .FALSE. + ELSEIF(IPRINT.EQ.2)THEN + TRKVERBOSE = .TRUE. + TRKDEBUG = .TRUE. ELSE - TRKDEBUG = .FALSE. + TRKVERBOSE = .FALSE. + TRKDEBUG = .FALSE. ENDIF * ---------------------------------------------------------- +* evaluate average spatial resolution +* ---------------------------------------------------------- + AVRESX = RESXAV + AVRESY = RESYAV + DO IP=1,6 + IF( XGOOD(IP).EQ.1 )THEN + NX=NX+1 + AVRESX=AVRESX+RESX(IP) + ENDIF + IF(NX.NE.0)AVRESX=AVRESX/NX + IF( YGOOD(IP).EQ.1 )THEN + NY=NY+1 + AVRESY=AVRESY+RESY(IP) + ENDIF + IF(NX.NE.0)AVRESY=AVRESY/NY + ENDDO + +* ---------------------------------------------------------- * define ALTOL(5) ---> tolerances on state vector * * ---------------------------------------------------------- - FACT=10. !scale factor to define - !tolerance on alfa - ALTOL(1)=RESX(1)/FACT !al(1) = x - ALTOL(2)=RESY(1)/FACT !al(2) = y - ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta) - $ +RESY(1)**2)/44.51/FACT - ALTOL(4)=ALTOL(3) !al(4)=phi +* changed in order to evaluate energy-dependent +* tolerances on all 5 parameters + FACT=100. !scale factor to define tolerance on alfa c deflection error (see PDG) - DELETA1=0.01*RESX(1)/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) - DELETA2=0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) + DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) + DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) +c$$$ ALTOL(1) = AVRESX/FACT !al(1) = x +c$$$ ALTOL(2) = AVRESY/FACT !al(2) = y +c$$$ ALTOL(3) = DSQRT(AVRESX**2 !al(3)=sin(theta) +c$$$ $ +AVRESY**2)/44.51/FACT +c$$$ ALTOL(4) = ALTOL(3) !al(4)=phi +c deflection error (see PDG) +c$$$ DELETA1 = 0.01*AVRESX/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) +c$$$ DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) * ---------------------------------------------------------- * ISTEP=0 !num. steps to minimize chi^2 JFAIL=0 !error flag + +c elena-------- + CHI2 = 0 +c elena-------- + if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) + * * ----------------------- * START MINIMIZATION LOOP @@ -99,15 +135,11 @@ IF(JFAIL.NE.0) THEN IFAIL=1 CHI2=-9999. - if(TRKDEBUG) + if(TRKVERBOSE) $ PRINT *,'*** ERROR in mini *** wrong CHISQ' RETURN ENDIF -* CHI2_P=CHI2 - -c print*,'@@@@@ ',istep,' - ',al - COST=1e-7 DO I=1,5 DO J=1,5 @@ -123,7 +155,7 @@ *------------------------------------------------------------* CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant IF(IFA.NE.0) THEN !not positive-defined - if(TRKDEBUG)then + if(TRKVERBOSE)then PRINT *, $ '*** ERROR in mini ***'// $ 'on matrix inversion (not pos-def)' @@ -137,7 +169,7 @@ CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion * ******************************************* * find new value of AL-pha -* +* ******************************************* DO I=1,5 DAL(I)=0. DO J=1,5 @@ -161,7 +193,7 @@ ENDDO CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) IF(IFA.NE.0) THEN - if(TRKDEBUG)then + if(TRKVERBOSE)then PRINT *, $ '*** ERROR in mini ***'// $ 'on matrix inversion (not pos-def)' @@ -173,6 +205,9 @@ RETURN ENDIF CALL DSFINV(4,CHI2DD_R,4) +* ******************************************* +* find new value of AL-pha +* ******************************************* DO I=1,4 DAL(I)=0. DO J=1,4 @@ -185,14 +220,17 @@ AL(I)=AL(I)+DAL(I) ENDDO ENDIF + + if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) + *------------------------------------------------------------* * ---------------------------------------------------- * *------------------------------------------------------------* -* ******************************************* * check parameter bounds: +*------------------------------------------------------------* DO I=1,5 IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN - if(TRKDEBUG)then + if(TRKVERBOSE)then PRINT*,' *** WARNING in mini *** ' PRINT*,'MINI_2 ==> AL(',I,') out of range' PRINT*,' value: ',AL(I), @@ -205,21 +243,35 @@ RETURN ENDIF ENDDO - +*------------------------------------------------------------* * check number of steps: - IF(ISTEP.gt.ISTEPMAX) then +*------------------------------------------------------------* + IF(ISTEP.ge.ISTEPMAX) then IFAIL=1 - if(TRKDEBUG) + if(TRKVERBOSE) $ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=', $ ISTEPMAX goto 11 endif +*------------------------------------------------------------* * --------------------------------------------- * evaluate deflection tolerance on the basis of * estimated deflection * --------------------------------------------- - ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT +*------------------------------------------------------------* +c$$$ ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT + ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT + ALTOL(1) = ALTOL(5)/DELETA1 + ALTOL(2) = ALTOL(1) + ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51 + ALTOL(4) = ALTOL(3) + *---- check tolerances: +c$$$ DO I=1,5 +c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! +c$$$ ENDDO +c$$$ print*,'chi2 -- ',DCHI2 + DO I=1,5 IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! ENDDO @@ -229,9 +281,9 @@ CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives IF(JFAIL.NE.0) THEN IFAIL=1 - if(TRKDEBUG)THEN + if(TRKVERBOSE)THEN CHI2=-9999. - if(TRKDEBUG) + if(TRKVERBOSE) $ PRINT *,'*** ERROR in mini *** wrong CHISQ' ENDIF RETURN @@ -246,7 +298,7 @@ IF(PFIXED.EQ.0.) THEN CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant IF(IFA.NE.0) THEN !not positive-defined - if(TRKDEBUG)then + if(TRKVERBOSE)then PRINT *, $ '*** ERROR in mini ***'// $ 'on matrix inversion (not pos-def)' @@ -273,7 +325,7 @@ ENDDO CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) IF(IFA.NE.0) THEN - if(TRKDEBUG)then + if(TRKVERBOSE)then PRINT *, $ '*** ERROR in mini ***'// $ 'on matrix inversion (not pos-def)' @@ -306,13 +358,18 @@ if(pfixed.ne.0.) ndof=ndof-4 ! ***PP*** if(ndof.le.0.) then ndof = 1 - if(TRKDEBUG) + if(TRKVERBOSE) $ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' endif + + if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5) + * ------------------------------------ * Reduced chi^2 CHI2 = CHI2/dble(ndof) +c print*,'mini2: chi2 ',chi2 + 11 CONTINUE NSTEP=ISTEP ! ***PP*** @@ -343,8 +400,10 @@ $ ,XV0(nplanes),YV0(nplanes) DIMENSION AL_P(5) - LOGICAL TRKDEBUG - COMMON/TRKD/TRKDEBUG +c LOGICAL TRKVERBOSE +c COMMON/TRKD/TRKVERBOSE + LOGICAL TRKDEBUG,TRKVERBOSE + COMMON/TRKD/TRKDEBUG,TRKVERBOSE * * chi^2 computation * @@ -354,7 +413,7 @@ JFAIL=0 !error flag CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes IF(JFAIL.NE.0) THEN - IF(TRKDEBUG) + IF(TRKVERBOSE) $ PRINT *,'CHISQ ==> error from trk routine POSXYZ !!' IFAIL=1 RETURN @@ -426,7 +485,7 @@ JFAIL=0 CALL POSXYZ(AL_P,JFAIL) IF(JFAIL.NE.0) THEN - IF(TRKDEBUG) + IF(TRKVERBOSE) *23456789012345678901234567890123456789012345678901234567890123456789012 $ PRINT *,'CHISQ ==> error from trk routine POSXYZ' IFAIL=1 @@ -440,7 +499,7 @@ JFAIL=0 CALL POSXYZ(AL_P,JFAIL) IF(JFAIL.NE.0) THEN - IF(TRKDEBUG) + IF(TRKVERBOSE) $ PRINT *,'CHISQ ==> error from trk routine POSXYZ' IFAIL=1 RETURN @@ -472,7 +531,7 @@ COSTHE=DSQRT(1.-AL(3)**2) IF(COSTHE.EQ.0.) THEN - IF(TRKDEBUG)PRINT *,'=== WARNING ===> COSTHE=0' + IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0' IFAIL=1 RETURN ENDIF @@ -559,8 +618,10 @@ include 'commontracker.f' !tracker general common include 'common_mini_2.f' !common for the tracking procedure - LOGICAL TRKDEBUG - COMMON/TRKD/TRKDEBUG +c LOGICAL TRKVERBOSE +c COMMON/TRKD/TRKVERBOSE + LOGICAL TRKDEBUG,TRKVERBOSE + COMMON/TRKD/TRKDEBUG,TRKVERBOSE c DIMENSION AL_P(5) * @@ -590,12 +651,16 @@ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) IF(VOUT(3).GT.VECT(3)) THEN IFAIL=1 - if(TRKDEBUG) + if(TRKVERBOSE) $ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' - if(.TRUE.)print*,'charge',charge - if(.TRUE.)print*,'vect',vect - if(.TRUE.)print*,'vout',vout - if(.TRUE.)print*,'step',step +c$$$ if(.TRUE.)print*,'charge',charge +c$$$ if(.TRUE.)print*,'vect',vect +c$$$ if(.TRUE.)print*,'vout',vout +c$$$ if(.TRUE.)print*,'step',step + if(TRKVERBOSE)print*,'charge',charge + if(TRKVERBOSE)print*,'vect',vect + if(TRKVERBOSE)print*,'vout',vout + if(TRKVERBOSE)print*,'step',step RETURN ENDIF Z=VOUT(3) @@ -664,3 +729,134 @@ return end + + +*************************************************** +* * +* * +* * +* * +* * +* * +************************************************** + + subroutine guess() + +c IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + include 'commontracker.f' !tracker general common + include 'common_mini_2.f' !common for the tracking procedure + + REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES) + REAL*4 CHI,XC,ZC,RADIUS +* ---------------------------------------- +* Y view +* ---------------------------------------- +* ---------------------------------------- +* initial guess with a straigth line +* ---------------------------------------- + SZZ=0. + SZY=0. + SSY=0. + SZ=0. + S1=0. + DO I=1,nplanes + IF(YGOOD(I).EQ.1)THEN + YY = YM(I) + IF(XGOOD(I).EQ.0)THEN + YY = (YM_A(I) + YM_B(I))/2 + ENDIF + SZZ=SZZ+ZM(I)*ZM(I) + SZY=SZY+ZM(I)*YY + SSY=SSY+YY + SZ=SZ+ZM(I) + S1=S1+1. + ENDIF + ENDDO + DET=SZZ*S1-SZ*SZ + AY=(SZY*S1-SZ*SSY)/DET + BY=(SZZ*SSY-SZY*SZ)/DET + Y0 = AY*ZINI+BY +* ---------------------------------------- +* X view +* ---------------------------------------- +* ---------------------------------------- +* 1) initial guess with a circle +* ---------------------------------------- + NP=0 + DO I=1,nplanes + IF(XGOOD(I).EQ.1)THEN + XX = XM(I) + IF(YGOOD(I).EQ.0)THEN + XX = (XM_A(I) + XM_B(I))/2 + ENDIF + NP=NP+1 + XP(NP)=XX + ZP(NP)=ZM(I) + ENDIF + ENDDO + CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) +c print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI + IF(IFLAG.NE.0)GOTO 10 !straigth fit + ARG = RADIUS**2-(ZINI-ZC)**2 + IF(ARG.LT.0)GOTO 10 !straigth fit + DC = SQRT(ARG) + IF(XC.GT.0)DC=-DC + X0=XC+DC + AX = -(ZINI-ZC)/DC + DEF=100./(RADIUS*0.3*0.43) + IF(XC.GT.0)DEF=-DEF + GOTO 20 !guess is ok + +* ---------------------------------------- +* 2) initial guess with a straigth line +* - if circle does not intersect reference plane +* - if bad chi**2 +* ---------------------------------------- + 10 CONTINUE + SZZ=0. + SZX=0. + SSX=0. + SZ=0. + S1=0. + DO I=1,nplanes + IF(XGOOD(I).EQ.1)THEN + XX = XM(I) + IF(YGOOD(I).EQ.0)THEN + XX = (XM_A(I) + XM_B(I))/2 + ENDIF + SZZ=SZZ+ZM(I)*ZM(I) + SZX=SZX+ZM(I)*XX + SSX=SSX+XX + SZ=SZ+ZM(I) + S1=S1+1. + ENDIF + ENDDO + DET=SZZ*S1-SZ*SZ + AX=(SZX*S1-SZ*SSX)/DET + BX=(SZZ*SSX-SZX*SZ)/DET + DEF = 0 + X0 = AX*ZINI+BX + + 20 CONTINUE +* ---------------------------------------- +* guess +* ---------------------------------------- + + AL(1) = X0 + AL(2) = Y0 + tath = sqrt(AY**2+AX**2) + AL(3) = tath/sqrt(1+tath**2) + IF(AX.NE.0)THEN + AL(4)= atan(AY/AX) + ELSE + AL(4) = acos(-1.)/2 + IF(AY.LT.0)AL(4) = AL(4)+acos(-1.) + ENDIF + IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4) + AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs + AL(5) = DEF + +c print*,' guess: ',(al(i),i=1,5) + + end