************************************************************************ * * subroutine to evaluate the vector alfa (AL) * which minimizes CHI^2 * * - modified from mini.f in order to call differente chi^2 routine. * The new one includes also single clusters: in this case * the residual is defined as the distance between the track and the * segment AB associated to the single cluster. * * ************************************************************************ SUBROUTINE MINI_2(ISTEP,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) include '../common/commontracker.f' !tracker general common include '../common/common_mini_2.f' !common for the tracking procedure logical DEBUG common/dbg/DEBUG parameter (inf=1.e8) !just a huge number... c------------------------------------------------------------------------ c variables used in the tracking procedure (mini and its subroutines) c c N.B.: in mini & C. (and in the following block of variables too) c the plane ordering is reversed in respect of normal 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 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 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/ !" DIMENSION DAL(5) !increment of vector alfa INTEGER IFLAG c-------------------------------------------------------- c IFLAG =1 ---- chi2 derivatives computed by using c incremental ratios and posxyz.f c IFLAG =2 ---- the approximation of Golden is used c (see chisq.f) c c NB: the two metods gives equivalent results BUT c method 2 is faster!! c-------------------------------------------------------- DATA IFLAG/2/ * ---------------------------------------------------------- * 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 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) * ---------------------------------------------------------- * 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 c print*,'@@@@@ ',istep,' - ',al 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 *------------------------------------------------------------* * 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 * ******************************************* * 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 !* * ******************************************* * 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** ' PRINT*,'MINI_2 ==> AL(',I,') out of range' PRINT*,' value: ',AL(I), $ ' limits: ',ALMIN(I),ALMAX(I) print*,'istep ',istep endif 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 goto 11 endif * --------------------------------------------- * evaluate deflection tolerance on the basis of * estimated deflection * --------------------------------------------- ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT *---- check tolerances: DO I=1,5 IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! ENDDO * ------------------------------------ * Number of Degree Of Freedom ndof=0 do ip=1,nplanes ndof=ndof $ +int(xgood(ip)) $ +int(ygood(ip)) enddo ndof=ndof-5 * ------------------------------------ * Reduced chi^2 CHI2 = CHI2/dble(ndof) 11 CONTINUE 101 CONTINUE c print*,'END MINI' RETURN END ****************************************************************************** * * routine to compute chi^2 and its derivatives * * * (modified in respect to the previous one in order to include * single clusters. In this case the residual is evaluated by * calculating the distance between the track intersection and the * segment AB associated to the single cluster) * ****************************************************************************** SUBROUTINE CHISQ(IFLAG,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) include '../common/commontracker.f' !tracker general common include '../common/common_mini_2.f' !common for the tracking procedure DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes) $ ,XV0(nplanes),YV0(nplanes) DIMENSION AL_P(5) * * chi^2 computation * DO I=1,5 AL_P(I)=AL(I) ENDDO 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 !!' IFAIL=1 RETURN ENDIF DO I=1,nplanes XV0(I)=XV(I) YV0(I)=YV(I) ENDDO * ------------------------------------------------ c$$$ CHI2=0. c$$$ DO I=1,nplanes c$$$ CHI2=CHI2 c$$$ + +(XV(I)-XM(I))**2/RESX(i)**2 *XGOOD(I)*YGOOD(I) c$$$ + +(YV(I)-YM(I))**2/RESY(i)**2 *YGOOD(I)*XGOOD(I) c$$$ ENDDO * --------------------------------------------------------- * For planes with only a X or Y-cl included, instead of * a X-Y couple, the residual for chi^2 calculation is * evaluated by finding the point x-y, along the segment AB, * closest to the track. * The X or Y coordinate, respectivelly for X and Y-cl, is * then assigned to XM or YM, which is then considered the * measured position of the cluster. * --------------------------------------------------------- CHI2=0. DO I=1,nplanes IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I)) ALFA = XM_A(I) - BETA * YM_A(I) YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2) if(YM(I).lt.dmin1(YM_A(I),YM_B(I))) $ YM(I)=dmin1(YM_A(I),YM_B(I)) if(YM(I).gt.dmax1(YM_A(I),YM_B(I))) $ YM(I)=dmax1(YM_A(I),YM_B(I)) XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I)) ALFA = YM_A(I) - BETA * XM_A(I) XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2) if(XM(I).lt.dmin1(XM_A(I),XM_B(I))) $ XM(I)=dmin1(XM_A(I),XM_B(I)) if(XM(I).gt.dmax1(XM_A(I),XM_B(I))) $ XM(I)=dmax1(XM_A(I),XM_B(I)) YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates ENDIF CHI2=CHI2 + +(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) ) + +(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) ) + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 + *( XGOOD(I)*(1-YGOOD(I)) ) + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 + *( (1-XGOOD(I))*YGOOD(I) ) ENDDO * ------------------------------------------------ * * calculation of derivatives (dX/dAL_fa and dY/dAL_fa) * * ////////////////////////////////////////////////// * METHOD 1 -- incremental ratios * ////////////////////////////////////////////////// IF(IFLAG.EQ.1) THEN DO J=1,5 DO JJ=1,5 AL_P(JJ)=AL(JJ) ENDDO AL_P(J)=AL_P(J)+STEPAL(J)/2. JFAIL=0 CALL POSXYZ(AL_P,JFAIL) IF(JFAIL.NE.0) THEN PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' IFAIL=1 RETURN ENDIF DO I=1,nplanes XV2(I)=XV(I) YV2(I)=YV(I) ENDDO AL_P(J)=AL_P(J)-STEPAL(J) JFAIL=0 CALL POSXYZ(AL_P,JFAIL) IF(JFAIL.NE.0) THEN PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' IFAIL=1 RETURN ENDIF DO I=1,nplanes XV1(I)=XV(I) YV1(I)=YV(I) ENDDO DO I=1,nplanes DXDAL(I,J)=(XV2(I)-XV1(I))/STEPAL(J) DYDAL(I,J)=(YV2(I)-YV1(I))/STEPAL(J) ENDDO ENDDO ENDIF * ////////////////////////////////////////////////// * METHOD 2 -- Bob Golden * ////////////////////////////////////////////////// IF(IFLAG.EQ.2) THEN DO I=1,nplanes DXDAL(I,1)=1. DYDAL(I,1)=0. DXDAL(I,2)=0. DYDAL(I,2)=1. COSTHE=DSQRT(1.-AL(3)**2) IF(COSTHE.EQ.0.) THEN PRINT *,'=== WARNING ===> COSTHE=0' STOP ENDIF DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3 DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE IF(AL(5).NE.0.) THEN DXDAL(I,5)= + (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I)) + *DCOS(AL(4))))/AL(5) DYDAL(I,5)= + (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I)) + *DSIN(AL(4))))/AL(5) ELSE DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 ) DYDAL(I,5)=0. ENDIF ENDDO ENDIF * * x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x * >>> CHI2D evaluation * DO J=1,5 CHI2D(J)=0. DO I=1,nplanes CHI2D(J)=CHI2D(J) + +2.*(XV0(I)-XM(I))/RESX(i)**2*DXDAL(I,J) *XGOOD(I) + +2.*(YV0(I)-YM(I))/RESY(i)**2*DYDAL(I,J) *YGOOD(I) ENDDO ENDDO * * >>> CHI2DD evaluation * DO I=1,5 DO J=1,5 CHI2DD(I,J)=0. DO K=1,nplanes CHI2DD(I,J)=CHI2DD(I,J) + +2.*DXDAL(K,I)*DXDAL(K,J)/RESX(k)**2 *XGOOD(K) + +2.*DYDAL(K,I)*DYDAL(K,J)/RESY(k)**2 *YGOOD(K) ENDDO ENDDO ENDDO * x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x RETURN END ***************************************************************** * * Routine to compute the track intersection points * on the tracking-system planes, given the track parameters * * The routine is based on GRKUTA, which computes the * trajectory of a charged particle in a magnetic field * by solving the equatins of motion with Runge-Kuta method. * * Variables that have to be assigned when the subroutine * is called are: * * ZM(1,NPLANES) ----> z coordinates of the planes * AL_P(1,5) ----> track-parameter vector * * ----------------------------------------------------------- * NB !!! * The routine works properly only if the * planes are numbered in descending order starting from the * reference plane (ZINI) * ----------------------------------------------------------- * ***************************************************************** SUBROUTINE POSXYZ(AL_P,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) include '../common/commontracker.f' !tracker general common include '../common/common_mini_2.f' !common for the tracking procedure c DIMENSION AL_P(5) * DO I=1,nplanes ZV(I)=ZM(I) ! ENDDO * * set parameters for GRKUTA * IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5)) IF(AL_P(5).EQ.0) CHARGE=1. VOUT(1)=AL_P(1) VOUT(2)=AL_P(2) VOUT(3)=ZINI ! DBLE(Z0)-DBLE(ZSPEC) VOUT(4)=AL_P(3)*DCOS(AL_P(4)) VOUT(5)=AL_P(3)*DSIN(AL_P(4)) VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 DO I=1,nplanes step=vout(3)-zv(i) 10 DO J=1,7 VECT(J)=VOUT(J) VECTINI(J)=VOUT(J) ENDDO 11 continue CALL GRKUTA(CHARGE,STEP,VECT,VOUT) IF(VOUT(3).GT.VECT(3)) THEN IFAIL=1 PRINT *,'posxy (grkuta): WARNING ===> backward track!!' print*,'charge',charge print*,'vect',vect print*,'vout',vout print*,'step',step RETURN ENDIF Z=VOUT(3) IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100 IF(Z.GT.ZM(I)+TOLL) GOTO 10 IF(Z.LE.ZM(I)-TOLL) THEN STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3)) DO J=1,7 VECT(J)=VECTINI(J) ENDDO GOTO 11 ENDIF * ----------------------------------------------- * evaluate track coordinates 100 XV(I)=VOUT(1) YV(I)=VOUT(2) ZV(I)=VOUT(3) AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.) AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) * ----------------------------------------------- ENDDO RETURN END * ********************************************************** * Some initialization routines * ********************************************************** * ---------------------------------------------------------- * Routine to initialize COMMON/TRACK/ * subroutine track_init IMPLICIT DOUBLE PRECISION (A-H,O-Z) include '../common/commontracker.f' !tracker general common include '../common/common_mini_2.f' !common for the tracking procedure include '../common/common_mech.f' do i=1,5 AL(i) = 0. enddo do ip=1,NPLANES ZM(IP) = fitz(nplanes-ip+1) !init to mech. position XM(IP) = -100. !0. YM(IP) = -100. !0. XM_A(IP) = -100. !0. YM_A(IP) = -100. !0. c ZM_A(IP) = 0 XM_B(IP) = -100. !0. YM_B(IP) = -100. !0. c ZM_B(IP) = 0 RESX(IP) = 1000. !3.d-4 RESY(IP) = 1000. !12.d-4 XGOOD(IP) = 0 YGOOD(IP) = 0 enddo return end