36 |
c 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 |
37 |
|
|
38 |
DATA STEPAL/5*1.d-7/ !alpha vector step |
DATA STEPAL/5*1.d-7/ !alpha vector step |
39 |
DATA ISTEPMAX/120/ !maximum number of steps in the chi^2 minimization |
DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization |
40 |
DATA TOLL/1.d-8/ !tolerance in reaching the next plane during |
DATA TOLL/1.d-8/ !tolerance in reaching the next plane during |
41 |
* !the tracking procedure |
* !the tracking procedure |
42 |
DATA STEPMAX/100./ !maximum number of steps in the trackin gprocess |
DATA STEPMAX/100./ !maximum number of steps in the trackin gprocess |
43 |
|
|
44 |
|
c DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components |
45 |
|
c DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" |
46 |
DATA ALMAX/dinf,dinf,dinf,dinf,dinf/ !limits on alpha vector components |
DATA ALMAX/dinf,dinf,dinf,dinf,dinf/ !limits on alpha vector components |
47 |
DATA ALMIN/-dinf,-dinf,-dinf,-dinf,-dinf/ !" |
DATA ALMIN/-dinf,-dinf,-dinf,-dinf,-dinf/ !" |
48 |
|
|
140 |
RETURN |
RETURN |
141 |
ENDIF |
ENDIF |
142 |
|
|
143 |
COST=1e-9 |
COST=1e-7 |
|
costfac=1.1 |
|
|
flagstep=0 |
|
|
|
|
|
7 continue |
|
|
|
|
|
IF(costfac.le.1.) THEN |
|
|
IFAIL=1 |
|
|
PRINT *,'=== WARNING ===> no matrix inversion ' |
|
|
RETURN |
|
|
ENDIF |
|
|
|
|
144 |
DO I=1,5 |
DO I=1,5 |
145 |
DO J=1,5 |
DO J=1,5 |
146 |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
154 |
* track fitting with FREE deflection |
* track fitting with FREE deflection |
155 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
156 |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
157 |
IF(IFA.NE.0.or.jfa.ne.0) THEN !not positive-defined |
IF(IFA.NE.0) THEN !not positive-defined |
158 |
if(ifa.eq.-1)then |
if(TRKVERBOSE)then |
|
if(TRKVERBOSE)then |
|
159 |
PRINT *, |
PRINT *, |
160 |
$ '*** ERROR in mini ***'// |
$ '*** ERROR in mini ***'// |
161 |
$ 'on matrix inversion (not pos-def)' |
$ 'on matrix inversion (not pos-def)' |
162 |
$ ,DET |
$ ,DET |
163 |
endif |
endif |
164 |
IF(CHI2.EQ.0) CHI2=-9999. |
IF(CHI2.EQ.0) CHI2=-9999. |
165 |
IF(CHI2.GT.0) CHI2=-CHI2 |
IF(CHI2.GT.0) CHI2=-CHI2 |
166 |
IFAIL=1 |
IFAIL=1 |
167 |
RETURN |
RETURN |
|
endif |
|
|
if(jfa.eq.-1)then |
|
|
if(flagstep.eq.-1.) costfac=(costfac-1)/2+1 |
|
|
cost=cost*costfact |
|
|
flagstep=1. |
|
|
goto 7 |
|
|
elseif(jfa.eq.1)then |
|
|
if(flagstep.eq.1.) costfac=(costfac-1)/2+1 |
|
|
cost=cost/costfac |
|
|
flagstep=-1. |
|
|
goto 7 |
|
|
endif |
|
168 |
ENDIF |
ENDIF |
169 |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
170 |
* ******************************************* |
* ******************************************* |
247 |
* check number of steps: |
* check number of steps: |
248 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
249 |
IF(ISTEP.ge.ISTEPMAX) then |
IF(ISTEP.ge.ISTEPMAX) then |
250 |
IFAIL=1 |
c$$$ IFAIL=1 |
251 |
if(TRKVERBOSE) |
c$$$ if(TRKVERBOSE) |
252 |
$ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=', |
c$$$ $ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=', |
253 |
$ ISTEPMAX |
c$$$ $ ISTEPMAX |
254 |
goto 11 |
goto 11 |
255 |
endif |
endif |
256 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
800 |
ENDIF |
ENDIF |
801 |
ENDDO |
ENDDO |
802 |
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
803 |
c print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI |
print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG |
804 |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
805 |
|
if(CHI.gt.100)GOTO 10 !straigth fit |
806 |
ARG = RADIUS**2-(ZINI-ZC)**2 |
ARG = RADIUS**2-(ZINI-ZC)**2 |
807 |
IF(ARG.LT.0)GOTO 10 !straigth fit |
IF(ARG.LT.0)GOTO 10 !straigth fit |
808 |
DC = SQRT(ARG) |
DC = SQRT(ARG) |
811 |
AX = -(ZINI-ZC)/DC |
AX = -(ZINI-ZC)/DC |
812 |
DEF=100./(RADIUS*0.3*0.43) |
DEF=100./(RADIUS*0.3*0.43) |
813 |
IF(XC.GT.0)DEF=-DEF |
IF(XC.GT.0)DEF=-DEF |
814 |
|
|
815 |
|
IF(ABS(X0).GT.30)THEN |
816 |
|
PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS |
817 |
|
$ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF |
818 |
|
GOTO 10 !straigth fit |
819 |
|
ENDIF |
820 |
GOTO 20 !guess is ok |
GOTO 20 !guess is ok |
821 |
|
|
822 |
* ---------------------------------------- |
* ---------------------------------------- |