--- DarthVader/TrackerLevel2/src/F77/mini.f 2006/05/30 16:30:37 1.2 +++ DarthVader/TrackerLevel2/src/F77/mini.f 2006/10/26 16:22:37 1.3 @@ -12,7 +12,7 @@ ************************************************************************ - SUBROUTINE MINI_2(ISTEP,IFAIL) + SUBROUTINE MINI2(ISTEP,IFAIL,IPRINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) @@ -31,9 +31,9 @@ c ordering, but they maintain their Z coordinates. so plane number 1 is c the first one that a particle meets, and its Z coordinate is > 0 c------------------------------------------------------------------------ - DATA ZINI/23.5/ !z coordinate of the reference plane + DATA ZINI/23.5/ !!! ***PP*** to be changed !z coordinate of the reference plane - DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking +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 @@ -48,6 +48,7 @@ DIMENSION DAL(5) !increment of vector alfa + DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2 INTEGER IFLAG c-------------------------------------------------------- @@ -59,7 +60,16 @@ c NB: the two metods gives equivalent results BUT c method 2 is faster!! c-------------------------------------------------------- - DATA IFLAG/2/ + DATA IFLAG/2/ + + LOGICAL TRKDEBUG + COMMON/TRKD/TRKDEBUG + + IF(IPRINT.EQ.1) THEN + TRKDEBUG = .TRUE. + ELSE + TRKDEBUG = .FALSE. + ENDIF * ---------------------------------------------------------- * define ALTOL(5) ---> tolerances on state vector @@ -79,86 +89,129 @@ * ISTEP=0 !num. steps to minimize chi^2 JFAIL=0 !error flag - CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives - IF(JFAIL.NE.0) THEN - IFAIL=1 - if(DEBUG) - $ PRINT *,'mini_2 ===> error on CHISQ computation !!!' - RETURN - ENDIF * * ----------------------- * START MINIMIZATION LOOP * ----------------------- 10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! - CHI2_P=CHI2 + + CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives + IF(JFAIL.NE.0) THEN + IFAIL=1 + CHI2=-9999. + if(TRKDEBUG) + $ PRINT *,'*** ERROR in mini *** wrong CHISQ' + RETURN + ENDIF + +* CHI2_P=CHI2 c print*,'@@@@@ ',istep,' - ',al - cost=1e-7 + COST=1e-7 DO I=1,5 DO J=1,5 CHI2DD(I,J)=CHI2DD(I,J)*COST ENDDO CHI2D(I)=CHI2D(I)*COST ENDDO + + IF(PFIXED.EQ.0.) THEN + *------------------------------------------------------------* * track fitting with FREE deflection *------------------------------------------------------------* - CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant - IF(IFA.NE.0) THEN !not positive-defined - if(DEBUG)then - PRINT *, - $ 'MINI_HOUGH ==> '// - $ '** ERROR ** on matrix inversion (not positive-defined)!!!' - $ ,DET - endif - IFAIL=1 - RETURN - ENDIF - CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion + CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant + IF(IFA.NE.0) THEN !not positive-defined + if(TRKDEBUG)then + PRINT *, + $ '*** ERROR in mini ***'// + $ 'on matrix inversion (not pos-def)' + $ ,DET + endif + IF(CHI2.EQ.0) CHI2=-9999. + IF(CHI2.GT.0) CHI2=-CHI2 + IFAIL=1 + RETURN + ENDIF + 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 !* - DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) !* - COV(I,J)=CHI2DD(I,J) !* - ENDDO !* - ENDDO !* - DO I=1,5 !* - AL(I)=AL(I)+DAL(I) !* - ENDDO !* +* find new value of AL-pha +* + DO I=1,5 + DAL(I)=0. + DO J=1,5 + DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) + COV(I,J)=2.*COST*CHI2DD(I,J) + ENDDO + ENDDO + DO I=1,5 + AL(I)=AL(I)+DAL(I) + ENDDO +*------------------------------------------------------------* +* track fitting with FIXED deflection +*------------------------------------------------------------* + ELSE + AL(5)=1./PFIXED + DO I=1,4 + CHI2D_R(I)=CHI2D(I) + DO J=1,4 + CHI2DD_R(I,J)=CHI2DD(I,J) + ENDDO + ENDDO + CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) + IF(IFA.NE.0) THEN + if(TRKDEBUG)then + PRINT *, + $ '*** ERROR in mini ***'// + $ 'on matrix inversion (not pos-def)' + $ ,DET + endif + IF(CHI2.EQ.0) CHI2=-9999. + IF(CHI2.GT.0) CHI2=-CHI2 + IFAIL=1 + RETURN + ENDIF + CALL DSFINV(4,CHI2DD_R,4) + DO I=1,4 + DAL(I)=0. + DO J=1,4 + DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) + COV(I,J)=2.*COST*CHI2DD_R(I,J) + ENDDO + ENDDO + DAL(5)=0. + DO I=1,4 + AL(I)=AL(I)+DAL(I) + ENDDO + ENDIF +*------------------------------------------------------------* +* ---------------------------------------------------- * +*------------------------------------------------------------* * ******************************************* * check parameter bounds: DO I=1,5 IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN - if(DEBUG)then - PRINT*,' **WARNING** ' + if(TRKDEBUG)then + PRINT*,' *** WARNING in mini *** ' PRINT*,'MINI_2 ==> AL(',I,') out of range' PRINT*,' value: ',AL(I), $ ' limits: ',ALMIN(I),ALMAX(I) print*,'istep ',istep endif + IF(CHI2.EQ.0) CHI2=-9999. + IF(CHI2.GT.0) CHI2=-CHI2 IFAIL=1 RETURN ENDIF ENDDO -* new estimate of chi^2: - JFAIL=0 !error flag - CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives - IF(JFAIL.NE.0) THEN - IFAIL=1 - if(DEBUG) - $ PRINT *,'mini_2: ===> error on CHISQ computation !!!' - RETURN - ENDIF + * check number of steps: IF(ISTEP.gt.ISTEPMAX) then IFAIL=1 - if(DEBUG) - $ PRINT *,'mini_2: WARNING ===> ISTEP.GT.ISTEPMAX=',ISTEPMAX + if(TRKDEBUG) + $ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=', + $ ISTEPMAX goto 11 endif * --------------------------------------------- @@ -171,6 +224,75 @@ IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! ENDDO +* new estimate of chi^2: + JFAIL=0 !error flag + CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives + IF(JFAIL.NE.0) THEN + IFAIL=1 + if(TRKDEBUG)THEN + CHI2=-9999. + if(TRKDEBUG) + $ PRINT *,'*** ERROR in mini *** wrong CHISQ' + ENDIF + RETURN + ENDIF + COST=1e-7 + DO I=1,5 + DO J=1,5 + CHI2DD(I,J)=CHI2DD(I,J)*COST + ENDDO + CHI2D(I)=CHI2D(I)*COST + ENDDO + 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 + PRINT *, + $ '*** ERROR in mini ***'// + $ 'on matrix inversion (not pos-def)' + $ ,DET + endif + IF(CHI2.EQ.0) CHI2=-9999. + IF(CHI2.GT.0) CHI2=-CHI2 + IFAIL=1 + RETURN + ENDIF + CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion + DO I=1,5 + DAL(I)=0. + DO J=1,5 + COV(I,J)=2.*COST*CHI2DD(I,J) + ENDDO + ENDDO + ELSE + DO I=1,4 + CHI2D_R(I)=CHI2D(I) + DO J=1,4 + CHI2DD_R(I,J)=CHI2DD(I,J) + ENDDO + ENDDO + CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) + IF(IFA.NE.0) THEN + if(TRKDEBUG)then + PRINT *, + $ '*** ERROR in mini ***'// + $ 'on matrix inversion (not pos-def)' + $ ,DET + endif + IF(CHI2.EQ.0) CHI2=-9999. + IF(CHI2.GT.0) CHI2=-CHI2 + IFAIL=1 + RETURN + ENDIF + CALL DSFINV(4,CHI2DD_R,4) + DO I=1,4 + DAL(I)=0. + DO J=1,4 + COV(I,J)=2.*COST*CHI2DD_R(I,J) + ENDDO + ENDDO + ENDIF +***************************** * ------------------------------------ * Number of Degree Of Freedom @@ -180,17 +302,20 @@ $ +int(xgood(ip)) $ +int(ygood(ip)) enddo - ndof=ndof-5 + if(pfixed.eq.0.) ndof=ndof-5 ! ***PP*** + if(pfixed.ne.0.) ndof=ndof-4 ! ***PP*** + if(ndof.le.0.) then + ndof = 1 + if(TRKDEBUG) + $ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' + endif * ------------------------------------ * Reduced chi^2 CHI2 = CHI2/dble(ndof) - 11 CONTINUE - 101 CONTINUE - -c print*,'END MINI' + NSTEP=ISTEP ! ***PP*** RETURN END @@ -217,6 +342,9 @@ DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes) $ ,XV0(nplanes),YV0(nplanes) DIMENSION AL_P(5) + + LOGICAL TRKDEBUG + COMMON/TRKD/TRKDEBUG * * chi^2 computation * @@ -226,7 +354,8 @@ JFAIL=0 !error flag CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes IF(JFAIL.NE.0) THEN - PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' + IF(TRKDEBUG) + $ PRINT *,'CHISQ ==> error from trk routine POSXYZ !!' IFAIL=1 RETURN ENDIF @@ -297,7 +426,9 @@ JFAIL=0 CALL POSXYZ(AL_P,JFAIL) IF(JFAIL.NE.0) THEN - PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' + IF(TRKDEBUG) +*23456789012345678901234567890123456789012345678901234567890123456789012 + $ PRINT *,'CHISQ ==> error from trk routine POSXYZ' IFAIL=1 RETURN ENDIF @@ -309,7 +440,8 @@ JFAIL=0 CALL POSXYZ(AL_P,JFAIL) IF(JFAIL.NE.0) THEN - PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' + IF(TRKDEBUG) + $ PRINT *,'CHISQ ==> error from trk routine POSXYZ' IFAIL=1 RETURN ENDIF @@ -340,8 +472,9 @@ COSTHE=DSQRT(1.-AL(3)**2) IF(COSTHE.EQ.0.) THEN - PRINT *,'=== WARNING ===> COSTHE=0' - STOP + IF(TRKDEBUG)PRINT *,'=== WARNING ===> COSTHE=0' + IFAIL=1 + RETURN ENDIF DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 @@ -425,6 +558,9 @@ include 'commontracker.f' !tracker general common include 'common_mini_2.f' !common for the tracking procedure + + LOGICAL TRKDEBUG + COMMON/TRKD/TRKDEBUG c DIMENSION AL_P(5) * @@ -454,12 +590,12 @@ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) IF(VOUT(3).GT.VECT(3)) THEN IFAIL=1 - if(WARNING) + if(TRKDEBUG) $ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' - if(WARNING)print*,'charge',charge - if(WARNING)print*,'vect',vect - if(WARNING)print*,'vout',vout - if(WARNING)print*,'step',step + if(.TRUE.)print*,'charge',charge + if(.TRUE.)print*,'vect',vect + if(.TRUE.)print*,'vout',vout + if(.TRUE.)print*,'step',step RETURN ENDIF Z=VOUT(3)