************************************************************************ * * 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... 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/-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 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 COMMON/TRKD/TRKDEBUG,TRKVERBOSE 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 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 * * ---------------------------------------------------------- * 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,1./AL(5) * * ----------------------- * START MINIMIZATION LOOP * ----------------------- 10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! 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 COST=1e-5 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(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) 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) 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,1./AL(5) *------------------------------------------------------------* * ---------------------------------------------------- * *------------------------------------------------------------* * 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 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 * 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(TRKVERBOSE)THEN CHI2=-9999. if(TRKVERBOSE) $ 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(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 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 DAL(I)=0. DO J=1,4 COV(I,J)=2.*COST*CHI2DD_R(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,1./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 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 cpp step=vout(3)-zv(i) step=vout(3)-zm(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 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 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. 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 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) c$$$ IF(AX.NE.0)THEN c$$$ AL(4)= atan(AY/AX) c$$$ ELSE c$$$ AL(4) = acos(-1.)/2 c$$$ IF(AY.LT.0)AL(4) = AL(4)+acos(-1.) c$$$ ENDIF c$$$ IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4) c$$$ AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys. c$$$ AL(4) = 0. c$$$ IF(AX.NE.0.AND.AY.NE.0)THEN c$$$ AL(4)= atan(AY/AX) c$$$ ELSEIF(AY.EQ.0)THEN c$$$ AL(4) = 0. c$$$ IF(AX.LT.0)AL(4) = AL(4)+acos(-1.) c$$$ ELSEIF(AX.EQ.0)THEN c$$$ AL(4) = acos(-1.)/2 c$$$ IF(AY.LT.0)AL(4) = AL(4)+acos(-1.) c$$$ ENDIF c$$$ IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4) c$$$ AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys. c$$$ AL(4)=0. c$$$ IF( AX.NE.0.OR.AY.NE.0. ) THEN c$$$ AL(4) = ASIN(AY/SQRT(AX**2+AY**2)) c$$$ IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4) c$$$ ENDIF 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