************************************************************************ * * 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 MINI2(ISTEP,IFAIL,IPRINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) include 'commontracker.f' !tracker general common include 'common_mini_2.f' !common for the tracking procedure c logical DEBUG c common/dbg/DEBUG parameter (dinf=1.d15) !just a huge number... parameter (dinfneg=-dinf) ! just a huge negative 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/ !!! ***PP*** to be changed !z coordinate of the reference plane 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 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/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components c DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components DATA ALMIN/dinfneg,dinfneg,-1.,dinfneg,dinfneg/ !" c$$$ 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 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/ c LOGICAL TRKDEBUG,TRKVERBOSE c COMMON/TRKD/TRKDEBUG,TRKVERBOSE LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT COMMON/TRKD/TRKDEBUG,TRKVERBOSE DIMENSION AL0(5) LOGICAL SUCCESS_NEW,SUCCESS_OLD c$$$ PRINT*,'==========' ! TEST c$$$ PRINT*,'START MINI' ! TEST c$$$ PRINT*,'==========' ! TEST * * define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student) * STUDENT = .false. FIRSTSTEPS = .true. FIRSTSTUDENT = .true. IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true. IF(IPRINT.EQ.1) THEN TRKVERBOSE = .TRUE. TRKDEBUG = .FALSE. ELSEIF(IPRINT.EQ.2)THEN TRKVERBOSE = .TRUE. TRKDEBUG = .TRUE. ELSE TRKVERBOSE = .FALSE. TRKDEBUG = .FALSE. ENDIF * ---------------------------------------------------------- * evaluate average spatial resolution * ---------------------------------------------------------- AVRESX = RESXAV AVRESY = RESYAV NX = 0.0 NY = 0.0 DO IP=1,6 IF( XGOOD(IP).EQ.1 )THEN NX=NX+1.0 AVRESX=AVRESX+RESX(IP) ENDIF IF( YGOOD(IP).EQ.1 )THEN NY=NY+1.0 AVRESY=AVRESY+RESY(IP) ENDIF ENDDO IF(NX.NE.0.0)AVRESX=AVRESX/NX IF(NY.NE.0.0)AVRESY=AVRESY/NY * ---------------------------------------------------------- * define ALTOL(5) ---> tolerances on state vector * * ---------------------------------------------------------- * changed in order to evaluate energy-dependent * tolerances on all 5 parameters cPP FACT=1.0e10 !scale factor to define tolerance on alfa c deflection error (see PDG) 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 CHI2=0 if(TRKDEBUG) print*,'guess: ',al if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5) * * ----------------------- * START MINIMIZATION LOOP * ----------------------- 10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! * ------------------------------- * **** Chi2+gaussian minimization * ------------------------------- IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN IF(ISTEP.GE.3) FIRSTSTEPS = .false. CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives IF(JFAIL.NE.0) THEN IFAIL=1 CHI2=-9999. if(TRKVERBOSE) $ PRINT *,'*** ERROR in mini *** wrong CHISQ' RETURN ENDIF c COST=1e-5 COST=1. DO I=1,5 IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2 ENDDO DO I=1,5 DO J=1,5 CHI2DD(I,J)=CHI2DD(I,J)*COST ENDDO c$$$ 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(TRKVERBOSE)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) *COST 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(TRKVERBOSE)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) * ******************************************* * find new value of AL-pha * ******************************************* DO I=1,4 DAL(I)=0. DO J=1,4 DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) *COST 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 if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5) c$$$ PRINT*,'DAL ',(DAL(K),K=1,5) c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5) ENDIF * ------------------------------- * **** Likelihood+Student minimization * ------------------------------- IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN IF(FIRSTSTUDENT) THEN FIRSTSTUDENT = .false. ISTEP = 1 ENDIF CALL CHISQSTT(1,JFAIL) DO I=1,5 DAL(I)=0. DO J=1,5 DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) ENDDO ENDDO DO I=1,5 DO j=1,5 COV(I,J) = 2.*CHI2DD(I,J) ENDDO ENDDO CHI2TOLL = 1.E-3 ALPHA = 3.0 BETA = -0.4 E=1. EA=1. EB=1. EC=1. FA=1. FB=1. FC=1. SUCCESS_OLD = .FALSE. SUCCESS_NEW = .FALSE. CALL CHISQSTT(0,JFAIL) c$$$ PRINT*,CHI2 CHI2_NEW = CHI2 FC = CHI2 EC = 0. ICOUNT = 0 100 CONTINUE ICOUNT = ICOUNT+1 DO I=1,5 AL0(I)=AL(I) ENDDO DO I=1,5 AL(I)=AL(I)+E*DAL(I) ENDDO CALL CHISQSTT(0,JFAIL) CHI2_OLD = CHI2_NEW CHI2_NEW = CHI2 FA = FB FB = FC FC = CHI2 EA = EB EB = EC EC = E c$$$ PRINT*,E,CHI2_NEW IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101 SUCCESS_OLD = SUCCESS_NEW SUCCESS_NEW = .TRUE. E = E*ALPHA ELSE ! failure SUCCESS_OLD = SUCCESS_NEW SUCCESS_NEW = .FALSE. CHI2_NEW = CHI2_OLD DO I=1,5 AL(I)=AL0(I) ENDDO IF(SUCCESS_OLD) THEN DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA) IF(DENOM.NE.0.) THEN E = EB - 0.5*( (EB-EA)**2*(FB-FC) $ - (EB-EC)**2*(FB-FA) ) / DENOM ELSE E = BETA*E ENDIF ELSE E = BETA*E ENDIF c$$$ E = BETA*E ENDIF IF(ICOUNT.GT.20) GOTO 101 GOTO 100 101 CONTINUE DO I=1,5 DAL(I)=E*DAL(I) ENDDO c$$$ print*,' ' c$$$ PRINT*,'DAL ',(DAL(K),K=1,5) c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5) c$$$ print*,'==== CHI2 ====' c$$$ print*,chi2 c$$$ print*,'==== CHI2d ====' c$$$ print*,(chi2d(i),i=1,5) c$$$ print*,'==== CHI2dd ====' c$$$ do j=1,5 c$$$ print*,(chi2dd(j,i),i=1,5) c$$$ enddo c$$$ print*,'================' c$$$ print*,' ' *========= FIN QUI ============= ENDIF *------------------------------------------------------------* * ---------------------------------------------------- * *------------------------------------------------------------* * check parameter bounds: *------------------------------------------------------------* DO I=1,5 IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN if(TRKVERBOSE)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 *------------------------------------------------------------* * check number of steps: *------------------------------------------------------------* IF(ISTEP.ge.ISTEPMAX) then c$$$ IFAIL=1 c$$$ if(TRKVERBOSE) c$$$ $ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=', c$$$ $ ISTEPMAX goto 11 endif *------------------------------------------------------------* * --------------------------------------------- * evaluate deflection tolerance on the basis of * estimated deflection * --------------------------------------------- *------------------------------------------------------------* c$$$ ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT IF(FACT.EQ.0)THEN IFAIL=1 RETURN ENDIF 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) c$$$ print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step! *---- check tolerances: c$$$ DO I=1,5 c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! c$$$ ENDDO c$$$ print*,'chi2 -- ',DCHI2 IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP*** DO I=1,5 IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! ENDDO ***************************** * final estimate of chi^2 ***************************** * ------------------------------- * **** Chi2+gaussian minimization * ------------------------------- IF(.NOT.STUDENT) THEN JFAIL=0 !error flag CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives IF(JFAIL.NE.0) THEN IFAIL=1 if(TRKVERBOSE)THEN CHI2=-9999. if(TRKVERBOSE) $ PRINT *,'*** ERROR in mini *** wrong CHISQ' ENDIF RETURN ENDIF c COST=1e-7 COST=1. DO I=1,5 IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2 ENDDO DO I=1,5 DO J=1,5 CHI2DD(I,J)=CHI2DD(I,J)*COST ENDDO 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(TRKVERBOSE)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 c$$$ 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(TRKVERBOSE)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 c$$$ DAL(I)=0. DO J=1,4 COV(I,J)=2.*COST*CHI2DD_R(I,J) ENDDO ENDDO ENDIF ENDIF * ------------------------------- * **** Likelihood+student minimization * ------------------------------- IF(STUDENT) THEN CALL CHISQSTT(1,JFAIL) DO I=1,5 DO j=1,5 COV(I,J) = 2.*CHI2DD(I,J) ENDDO ENDDO ENDIF ***************************** * ------------------------------------ * Number of Degree Of Freedom ndof=0 do ip=1,nplanes ndof=ndof $ +int(xgood(ip)) $ +int(ygood(ip)) enddo if(pfixed.eq.0.) ndof=ndof-5 ! ***PP*** if(pfixed.ne.0.) ndof=ndof-4 ! ***PP*** if(ndof.le.0.) then ndof = 1 if(TRKVERBOSE) $ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' endif * ------------------------------------ * Reduced chi^2 CHI2 = CHI2/dble(ndof) c print*,'mini2: chi2 ',chi2 11 CONTINUE if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,AL(5) NSTEP=ISTEP ! ***PP*** c$$$ print*,'>>>>> NSTEP = ',NSTEP 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 'commontracker.f' !tracker general common include '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) c LOGICAL TRKVERBOSE c COMMON/TRKD/TRKVERBOSE LOGICAL TRKDEBUG,TRKVERBOSE COMMON/TRKD/TRKDEBUG,TRKVERBOSE * * 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 IF(TRKVERBOSE) $ PRINT *,'CHISQ ==> error from trk 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) ) c$$$ print*,(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) ) c$$$ print*,(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) ) c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 c$$$ + *( XGOOD(I)*(1-YGOOD(I)) ) c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 c$$$ + *( (1-XGOOD(I))*YGOOD(I) ) c$$$ print*,XV(I),XM(I),XGOOD(I) c$$$ print*,YV(I),YM(I),YGOOD(I) ENDDO c$$$ print*,'CHISQ ',chi2 * ------------------------------------------------ * * 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 IF(TRKVERBOSE) *23456789012345678901234567890123456789012345678901234567890123456789012 $ PRINT *,'CHISQ ==> error from trk 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 IF(TRKVERBOSE) $ PRINT *,'CHISQ ==> error from trk 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 IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0' IFAIL=1 RETURN 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 Likelihodd+Student 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 CHISQSTT(IFLAG,JFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) include 'commontracker.f' !tracker general common include 'common_mini_2.f' !common for the tracking procedure LOGICAL TRKDEBUG,TRKVERBOSE COMMON/TRKD/TRKDEBUG,TRKVERBOSE DIMENSION AL_P(5) DIMENSION VECTEMP(5) c$$$ DIMENSION U(5) ! BFGS 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 IF(TRKVERBOSE) $ PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!' IFAIL=1 RETURN ENDIF 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 IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0' IFAIL=1 RETURN 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 IF(IFLAG.EQ.0) THEN ! function calulation 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 TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/ $ (TAILX(I)*RESX(I)**2) ) TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/ $ (TAILY(I)*RESY(I)**2) ) CHI2=CHI2 $ +(TAILX(I)+1.0)*TERMX *( XGOOD(I) ) $ +(TAILY(I)+1.0)*TERMY *( YGOOD(I) ) ENDDO ENDIF IF(IFLAG.EQ.1) THEN ! derivative calulation DO I=1,5 CHI2DOLD(I)=CHI2D(I) ENDDO DO J=1,5 CHI2D(J)=0. DO I=1,nplanes CHI2D(J)=CHI2D(J) $ +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/ $ (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)* $ DXDAL(I,J) *XGOOD(I) $ +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/ $ (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)* $ DYDAL(I,J) *YGOOD(I) ENDDO ENDDO DO K=1,5 VECTEMP(K)=0. DO M=1,5 VECTEMP(K) = VECTEMP(K) + $ COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M)) ENDDO ENDDO DOWN1 = 0. DO K=1,5 DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K)) ENDDO IF(DOWN1.EQ.0.) THEN PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0' IFAIL=1 RETURN ENDIF DOWN2 = 0. DO K=1,5 DO M=1,5 DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K) ENDDO ENDDO IF(DOWN2.EQ.0.) THEN PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0' IFAIL=1 RETURN ENDIF c$$$ DO K=1,5 ! BFGS c$$$ U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2 c$$$ ENDDO DO I=1,5 DO J=1,5 CHI2DD(I,J) = COV(I,J)/2. $ +DAL(I)*DAL(J)/DOWN1 $ -VECTEMP(I)*VECTEMP(J)/DOWN2 c$$$ $ +DOWN2*U(I)*U(J) ! BFGS ENDDO ENDDO ENDIF 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 'commontracker.f' !tracker general common include 'common_mini_2.f' !common for the tracking procedure c LOGICAL TRKVERBOSE c COMMON/TRKD/TRKVERBOSE LOGICAL TRKDEBUG,TRKVERBOSE COMMON/TRKD/TRKDEBUG,TRKVERBOSE c DIMENSION AL_P(5) * cpp DO I=1,nplanes cpp ZV(I)=ZM(I) ! cpp 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 c$$$ print*,'POSXY (prima) ',vout DO I=1,nplanes c$$$ ipass = 0 ! TEST c$$$ PRINT *,'TRACKING -> START PLANE: ',I ! TEST cPPP step=vout(3)-zm(i) cPP step=(zm(i)-vout(3))/VOUT(6) 10 DO J=1,7 VECT(J)=VOUT(J) VECTINI(J)=VOUT(J) ENDDO cPPP step=vect(3)-zm(i) IF(VOUT(6).GE.0.) THEN IFAIL=1 if(TRKVERBOSE) $ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' RETURN ENDIF step=(zm(i)-vect(3))/VOUT(6) 11 continue CALL GRKUTA(CHARGE,STEP,VECT,VOUT) c$$$ ipass = ipass + 1 ! TEST c$$$ PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST IF(VOUT(3).GT.VECT(3)) THEN IFAIL=1 if(TRKVERBOSE) $ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' 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) 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.) * ----------------------------------------------- IF(TRACKMODE.EQ.1) THEN * ----------------------------------------------- * change of energy by bremsstrahlung for electrons VOUT(7) = VOUT(7) * 0.997 !0.9968 * ----------------------------------------------- ENDIF c$$$ PRINT *,'TRACKING -> END' ! TEST ENDDO c$$$ print*,'POSXY (dopo) ',vout RETURN END * ********************************************************** * Some initialization routines * ********************************************************** * ---------------------------------------------------------- * Routine to initialize COMMON/TRACK/ * subroutine track_init IMPLICIT DOUBLE PRECISION (A-H,O-Z) include 'commontracker.f' !tracker general common include 'common_mini_2.f' !common for the tracking procedure include '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. ZM_A(IP) = fitz(nplanes-ip+1) !init to mech. position XM_B(IP) = -100. !0. YM_B(IP) = -100. !0. ZM_B(IP) = fitz(nplanes-ip+1) !init to mech. position RESX(IP) = 1000. !3.d-4 RESY(IP) = 1000. !12.d-4 XGOOD(IP) = 0 YGOOD(IP) = 0 DEDXTRK_X(IP) = 0 DEDXTRK_Y(IP) = 0 AXV(IP) = 0 AYV(IP) = 0 XV(IP) = -100 YV(IP) = -100 enddo 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 IFLAG=0 !no debug mode CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) c$$$ print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG c$$$ print*,' XP ',(xp(i),i=1,np) c$$$ print*,' ZP ',(zp(i),i=1,np) c$$$ print*,' AP ',(ap(i),i=1,np) c$$$ print*,' XP ',(rp(i),i=1,np) IF(IFLAG.NE.0)GOTO 10 !straigth fit c if(CHI.gt.100)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 IF(ABS(X0).GT.30)THEN c$$$ PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS c$$$ $ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF GOTO 10 !straigth fit ENDIF 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) AL(4)=0. IF( AX.NE.0.OR.AY.NE.0. ) THEN AL(4) = ASIN(AY/SQRT(AX**2+AY**2)) IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4) IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4) ENDIF IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0) IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0) AL(5) = DEF c print*,' guess: ',(al(i),i=1,5) end