/[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.9 by pam-fi, Tue Nov 21 17:12:26 2006 UTC revision 1.16 by pam-fi, Mon May 14 11:03:06 2007 UTC
# Line 43  c      DATA XGOOD,YGOOD/nplanes*1.,nplan Line 43  c      DATA XGOOD,YGOOD/nplanes*1.,nplan
43    
44  c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components  c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
45  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"  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,1.,dinf,dinf/ !limits on alpha vector components
47        DATA ALMIN/-dinf,-dinf,-dinf,-dinf,-dinf/ !"        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
48    
49        DIMENSION DAL(5)                    !increment of vector alfa        DIMENSION DAL(5)                    !increment of vector alfa
50        DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2        DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2
# Line 105  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 105  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
105  *     ----------------------------------------------------------  *     ----------------------------------------------------------
106  *     changed in order to evaluate energy-dependent  *     changed in order to evaluate energy-dependent
107  *     tolerances on all 5 parameters  *     tolerances on all 5 parameters
108        FACT=100.                  !scale factor to define tolerance on alfa  cPP      FACT=1.0e10                  !scale factor to define tolerance on alfa
109  c     deflection error (see PDG)  c     deflection error (see PDG)
110        DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))        DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
111        DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)        DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
# Line 121  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        if(TRKDEBUG) print*,'guess: ',al        if(TRKDEBUG) print*,'guess: ',al
127        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
# Line 140  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 141  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
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 266  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE Line 267  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE
267        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
268        ALTOL(4) = ALTOL(3)                ALTOL(4) = ALTOL(3)        
269                
270    c$$$      print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step!
271    
272  *---- check tolerances:  *---- check tolerances:
273  c$$$      DO I=1,5  c$$$      DO I=1,5
274  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
275  c$$$      ENDDO  c$$$      ENDDO
276  c$$$      print*,'chi2 -- ',DCHI2  c$$$      print*,'chi2 -- ',DCHI2
277    
278          IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP***
279        DO I=1,5        DO I=1,5
280           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
281        ENDDO        ENDDO
# Line 362  c$$$      print*,'chi2 -- ',DCHI2 Line 366  c$$$      print*,'chi2 -- ',DCHI2
366       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
367        endif        endif
368    
       if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)  
   
369  *     ------------------------------------  *     ------------------------------------
370  *     Reduced chi^2  *     Reduced chi^2
371        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
# Line 372  c      print*,'mini2: chi2 ',chi2 Line 374  c      print*,'mini2: chi2 ',chi2
374    
375   11   CONTINUE         11   CONTINUE      
376    
377          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
378    
379        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
380    
381    c$$$      print*,'>>>>> NSTEP = ',NSTEP
382    
383        RETURN        RETURN
384        END                END        
385    
# Line 439  c$$$      ENDDO Line 445  c$$$      ENDDO
445  *     measured position of the cluster.  *     measured position of the cluster.
446  *     ---------------------------------------------------------  *     ---------------------------------------------------------
447        CHI2=0.        CHI2=0.
448        DO I=1,nplanes        DO I=1,nplanes        
449           IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl           IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
450              BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))              BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
451              ALFA = XM_A(I) - BETA * YM_A(I)              ALFA = XM_A(I) - BETA * YM_A(I)
# Line 466  c$$$      ENDDO Line 472  c$$$      ENDDO
472       +                                       *( XGOOD(I)*(1-YGOOD(I)) )       +                                       *( XGOOD(I)*(1-YGOOD(I)) )
473       +        +((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
474       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
475    c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
476    c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
477    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
478    c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
479    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
480    c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
481    c$$$         print*,XV(I),XM(I),XGOOD(I)
482    c$$$         print*,YV(I),YM(I),YGOOD(I)
483        ENDDO        ENDDO
484  c      print*,'CHISQ ',chi2  c$$$      print*,'CHISQ ',chi2
485  *     ------------------------------------------------  *     ------------------------------------------------
486  *      *    
487  *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)      *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)    
# Line 626  c      COMMON/TRKD/TRKVERBOSE Line 640  c      COMMON/TRKD/TRKVERBOSE
640  c        c      
641        DIMENSION AL_P(5)        DIMENSION AL_P(5)
642  *      *    
643        DO I=1,nplanes  cpp      DO I=1,nplanes
644           ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
645        ENDDO        cpp      ENDDO      
646  *      *    
647  *     set parameters for GRKUTA  *     set parameters for GRKUTA
648  *      *    
# Line 643  c       Line 657  c      
657        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))
658        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
659    
660  c$$$      print*,'POSXY ',vout  c$$$      print*,'POSXY (prima) ',vout
661    
662        DO I=1,nplanes        DO I=1,nplanes
663           step=vout(3)-zv(i)  cpp         step=vout(3)-zv(i)
664             step=vout(3)-zm(i)
665   10      DO J=1,7   10      DO J=1,7
666              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
667              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
# Line 678  c$$$            if(.TRUE.)print*,'step', Line 693  c$$$            if(.TRUE.)print*,'step',
693              GOTO 11              GOTO 11
694           ENDIF           ENDIF
695    
696    
697  *     -----------------------------------------------  *     -----------------------------------------------
698  *        evaluate track coordinates  *        evaluate track coordinates
699   100     XV(I)=VOUT(1)   100     XV(I)=VOUT(1)
# Line 687  c$$$            if(.TRUE.)print*,'step', Line 703  c$$$            if(.TRUE.)print*,'step',
703           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
704  *     -----------------------------------------------  *     -----------------------------------------------
705    
706             IF(TRACKMODE.EQ.1) THEN
707    *     -----------------------------------------------
708    *        change of energy by bremsstrahlung for electrons
709                VOUT(7) = VOUT(7) * 0.997 !0.9968
710    *     -----------------------------------------------
711             ENDIF
712    
713        ENDDO        ENDDO
714    
715    c$$$      print*,'POSXY (dopo) ',vout
716    
717    
718        RETURN        RETURN
719        END        END
720    
# Line 729  c         ZM_B(IP) = 0 Line 755  c         ZM_B(IP) = 0
755           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
756           XGOOD(IP) = 0           XGOOD(IP) = 0
757           YGOOD(IP) = 0           YGOOD(IP) = 0
758             DEDXTRK_X(IP) = 0
759             DEDXTRK_Y(IP) = 0
760             AXV(IP) = 0
761             AYV(IP) = 0
762             XV(IP) = -100
763             YV(IP) = -100
764        enddo        enddo
765    
766        return        return
# Line 801  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 833  c      IMPLICIT DOUBLE PRECISION (A-H,O-
833        ENDDO        ENDDO
834        IFLAG=0                   !no debug mode        IFLAG=0                   !no debug mode
835        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
836  c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG  
837    c$$$      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
838    c$$$      print*,' XP ',(xp(i),i=1,np)
839    c$$$      print*,' ZP ',(zp(i),i=1,np)
840    c$$$      print*,' AP ',(ap(i),i=1,np)
841    c$$$      print*,' XP ',(rp(i),i=1,np)
842    
843        IF(IFLAG.NE.0)GOTO 10 !straigth fit        IF(IFLAG.NE.0)GOTO 10 !straigth fit
844        if(CHI.gt.100)GOTO 10 !straigth fit  c      if(CHI.gt.100)GOTO 10 !straigth fit
845        ARG = RADIUS**2-(ZINI-ZC)**2        ARG = RADIUS**2-(ZINI-ZC)**2
846        IF(ARG.LT.0)GOTO 10       !straigth fit        IF(ARG.LT.0)GOTO 10       !straigth fit
847        DC = SQRT(ARG)              DC = SQRT(ARG)      
# Line 813  c      print*,' circle: ',XC,ZC,RADIUS,' Line 851  c      print*,' circle: ',XC,ZC,RADIUS,'
851        DEF=100./(RADIUS*0.3*0.43)        DEF=100./(RADIUS*0.3*0.43)
852        IF(XC.GT.0)DEF=-DEF        IF(XC.GT.0)DEF=-DEF
853                
854    
855    
856        IF(ABS(X0).GT.30)THEN        IF(ABS(X0).GT.30)THEN
857           PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS  c$$$         PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
858       $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF  c$$$     $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
859           GOTO 10       !straigth fit           GOTO 10       !straigth fit
860        ENDIF        ENDIF
861        GOTO 20                   !guess is ok        GOTO 20                   !guess is ok
# Line 859  c      print*,' circle: ',XC,ZC,RADIUS,' Line 899  c      print*,' circle: ',XC,ZC,RADIUS,'
899        AL(2) = Y0        AL(2) = Y0
900        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
901        AL(3) = tath/sqrt(1+tath**2)        AL(3) = tath/sqrt(1+tath**2)
902        IF(AX.NE.0)THEN  c$$$      IF(AX.NE.0)THEN
903           AL(4)= atan(AY/AX)  c$$$         AL(4)= atan(AY/AX)
904        ELSE  c$$$      ELSE
905           AL(4) = acos(-1.)/2  c$$$         AL(4) = acos(-1.)/2
906           IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)  c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
907    c$$$  ENDIF
908    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
909    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
910    
911    c$$$      AL(4) =  0.
912    c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN
913    c$$$         AL(4)= atan(AY/AX)
914    c$$$      ELSEIF(AY.EQ.0)THEN
915    c$$$         AL(4) = 0.
916    c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)        
917    c$$$      ELSEIF(AX.EQ.0)THEN
918    c$$$         AL(4) = acos(-1.)/2
919    c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
920    c$$$      ENDIF
921    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
922    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
923    
924    c$$$      AL(4)=0.
925    c$$$      IF( AX.NE.0.OR.AY.NE.0. ) THEN
926    c$$$         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
927    c$$$         IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)
928    c$$$      ENDIF
929    
930          AL(4)=0.
931          IF( AX.NE.0.OR.AY.NE.0. ) THEN
932             AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
933             IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
934             IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
935        ENDIF        ENDIF
936        IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)        IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
937        AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs        IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)
938    
939        AL(5) = DEF        AL(5) = DEF
940    
941  c      print*,' guess: ',(al(i),i=1,5)  c      print*,' guess: ',(al(i),i=1,5)
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.23