/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/mini.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/mini.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by pam-fi, Tue Nov 7 15:55:11 2006 UTC revision 1.13 by pam-fi, Mon Feb 19 16:28:39 2007 UTC
# Line 36  c--------------------------------------- Line 36  c---------------------------------------
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,1.,dinf,dinf/ !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
47        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
48    
# Line 119  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 121  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
121  *      *    
122        ISTEP=0                   !num. steps to minimize chi^2        ISTEP=0                   !num. steps to minimize chi^2
123        JFAIL=0                   !error flag        JFAIL=0                   !error flag
124          CHI2=0
125    
126  c     elena--------        if(TRKDEBUG) print*,'guess: ',al
       CHI2   = 0  
 c     elena--------  
127        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
128    
129  *      *    
# Line 140  c     elena-------- Line 141  c     elena--------
141           RETURN           RETURN
142        ENDIF        ENDIF
143                
144        COST=1e-7        COST=1e-5
145        DO I=1,5        DO I=1,5
146           DO J=1,5           DO J=1,5
147              CHI2DD(I,J)=CHI2DD(I,J)*COST              CHI2DD(I,J)=CHI2DD(I,J)*COST
# Line 247  c     elena-------- Line 248  c     elena--------
248  *     check number of steps:  *     check number of steps:
249  *------------------------------------------------------------*  *------------------------------------------------------------*
250        IF(ISTEP.ge.ISTEPMAX) then        IF(ISTEP.ge.ISTEPMAX) then
251           IFAIL=1  c$$$         IFAIL=1
252           if(TRKVERBOSE)  c$$$         if(TRKVERBOSE)
253       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',  c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
254       $        ISTEPMAX  c$$$     $        ISTEPMAX
255           goto 11           goto 11
256        endif        endif
257  *------------------------------------------------------------*  *------------------------------------------------------------*
# Line 272  c$$$         if(TRKVERBOSE)print*,i,' -- Line 273  c$$$         if(TRKVERBOSE)print*,i,' --
273  c$$$      ENDDO  c$$$      ENDDO
274  c$$$      print*,'chi2 -- ',DCHI2  c$$$      print*,'chi2 -- ',DCHI2
275    
276          IF(ISTEP.LT.3) GOTO 10 ! ***PP***
277        DO I=1,5        DO I=1,5
278           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
279        ENDDO        ENDDO
# Line 466  c$$$      ENDDO Line 468  c$$$      ENDDO
468       +                                       *( XGOOD(I)*(1-YGOOD(I)) )       +                                       *( XGOOD(I)*(1-YGOOD(I)) )
469       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
470       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
471    c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
472    c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
473    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
474    c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
475    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
476    c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
477    c$$$         print*,XV(I),XM(I),XGOOD(I)
478    c$$$         print*,YV(I),YM(I),YGOOD(I)
479        ENDDO        ENDDO
480    c$$$      print*,'CHISQ ',chi2
481  *     ------------------------------------------------  *     ------------------------------------------------
482  *      *    
483  *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)      *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)    
# Line 641  c       Line 652  c      
652        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
653        IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))        IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
654        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
655    
656    c$$$      print*,'POSXY (prima) ',vout
657    
658        DO I=1,nplanes        DO I=1,nplanes
659           step=vout(3)-zv(i)           step=vout(3)-zv(i)
660   10      DO J=1,7   10      DO J=1,7
# Line 674  c$$$            if(.TRUE.)print*,'step', Line 688  c$$$            if(.TRUE.)print*,'step',
688              GOTO 11              GOTO 11
689           ENDIF           ENDIF
690    
691    
692  *     -----------------------------------------------  *     -----------------------------------------------
693  *        evaluate track coordinates  *        evaluate track coordinates
694   100     XV(I)=VOUT(1)   100     XV(I)=VOUT(1)
# Line 683  c$$$            if(.TRUE.)print*,'step', Line 698  c$$$            if(.TRUE.)print*,'step',
698           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
699  *     -----------------------------------------------  *     -----------------------------------------------
700    
701             IF(TRACKMODE.EQ.1) THEN
702    *     -----------------------------------------------
703    *        change of energy by bremsstrahlung for electrons
704                VOUT(7) = VOUT(7) * 0.997 !0.9968
705    *     -----------------------------------------------
706             ENDIF
707    
708        ENDDO        ENDDO
709    
710    c$$$      print*,'POSXY (dopo) ',vout
711    
712    
713        RETURN        RETURN
714        END        END
715    
# Line 795  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 820  c      IMPLICIT DOUBLE PRECISION (A-H,O-
820              ZP(NP)=ZM(I)              ZP(NP)=ZM(I)
821           ENDIF           ENDIF
822        ENDDO        ENDDO
823          IFLAG=0                   !no debug mode
824        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
825  c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI  c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
826        IF(IFLAG.NE.0)GOTO 10 !straigth fit        IF(IFLAG.NE.0)GOTO 10 !straigth fit
827          if(CHI.gt.100)GOTO 10 !straigth fit
828        ARG = RADIUS**2-(ZINI-ZC)**2        ARG = RADIUS**2-(ZINI-ZC)**2
829        IF(ARG.LT.0)GOTO 10       !straigth fit        IF(ARG.LT.0)GOTO 10       !straigth fit
830        DC = SQRT(ARG)              DC = SQRT(ARG)      
# Line 806  c      print*,' circle: ',XC,ZC,RADIUS,' Line 833  c      print*,' circle: ',XC,ZC,RADIUS,'
833        AX = -(ZINI-ZC)/DC        AX = -(ZINI-ZC)/DC
834        DEF=100./(RADIUS*0.3*0.43)        DEF=100./(RADIUS*0.3*0.43)
835        IF(XC.GT.0)DEF=-DEF        IF(XC.GT.0)DEF=-DEF
836          
837          IF(ABS(X0).GT.30)THEN
838    c$$$         PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
839    c$$$     $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
840             GOTO 10       !straigth fit
841          ENDIF
842        GOTO 20                   !guess is ok        GOTO 20                   !guess is ok
843    
844  *     ----------------------------------------  *     ----------------------------------------
# Line 847  c      print*,' circle: ',XC,ZC,RADIUS,' Line 880  c      print*,' circle: ',XC,ZC,RADIUS,'
880        AL(2) = Y0        AL(2) = Y0
881        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
882        AL(3) = tath/sqrt(1+tath**2)        AL(3) = tath/sqrt(1+tath**2)
883        IF(AX.NE.0)THEN  c$$$      IF(AX.NE.0)THEN
884           AL(4)= atan(AY/AX)  c$$$         AL(4)= atan(AY/AX)
885        ELSE  c$$$      ELSE
886           AL(4) = acos(-1.)/2  c$$$         AL(4) = acos(-1.)/2
887           IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)  c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
888    c$$$  ENDIF
889    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
890    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
891    
892    c$$$      AL(4) =  0.
893    c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN
894    c$$$         AL(4)= atan(AY/AX)
895    c$$$      ELSEIF(AY.EQ.0)THEN
896    c$$$         AL(4) = 0.
897    c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)        
898    c$$$      ELSEIF(AX.EQ.0)THEN
899    c$$$         AL(4) = acos(-1.)/2
900    c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
901    c$$$      ENDIF
902    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
903    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
904    
905    c$$$      AL(4)=0.
906    c$$$      IF( AX.NE.0.OR.AY.NE.0. ) THEN
907    c$$$         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
908    c$$$         IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)
909    c$$$      ENDIF
910    
911          AL(4)=0.
912          IF( AX.NE.0.OR.AY.NE.0. ) THEN
913             AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
914             IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
915             IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
916        ENDIF        ENDIF
917        IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)        IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
918        AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs        IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)
919    
920        AL(5) = DEF        AL(5) = DEF
921    
922  c      print*,' guess: ',(al(i),i=1,5)  c      print*,' guess: ',(al(i),i=1,5)
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.23