/[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.17 by pam-fi, Thu May 24 13:29:09 2007 UTC revision 1.23 by pam-fi, Fri Dec 5 08:26:47 2008 UTC
# Line 67  c--------------------------------------- Line 67  c---------------------------------------
67                    
68  c      LOGICAL TRKDEBUG,TRKVERBOSE  c      LOGICAL TRKDEBUG,TRKVERBOSE
69  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
70        LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT        LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT
71        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
72    
73        DIMENSION AL0(5)        DIMENSION AL0(5)
74        LOGICAL SUCCESS_NEW,SUCCESS_OLD        LOGICAL SUCCESS_NEW,SUCCESS_OLD
75    
76    c$$$      PRINT*,'==========' ! TEST
77    c$$$      PRINT*,'START MINI' ! TEST
78    c$$$      PRINT*,'==========' ! TEST
79    
80  *  *
81  *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)  *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)
82  *  *
83        STUDENT = .false.        STUDENT = .false.
84          FIRSTSTEPS = .true.
85          FIRSTSTUDENT = .true.
86        IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.        IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.
87    
88        IF(IPRINT.EQ.1) THEN        IF(IPRINT.EQ.1) THEN
# Line 94  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 101  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
101  *     ----------------------------------------------------------  *     ----------------------------------------------------------
102        AVRESX = RESXAV        AVRESX = RESXAV
103        AVRESY = RESYAV        AVRESY = RESYAV
104          NX = 0.0
105          NY = 0.0
106        DO IP=1,6        DO IP=1,6
107           IF( XGOOD(IP).EQ.1 )THEN           IF( XGOOD(IP).EQ.1 )THEN
108              NX=NX+1              NX=NX+1.0
109              AVRESX=AVRESX+RESX(IP)              AVRESX=AVRESX+RESX(IP)
110           ENDIF           ENDIF
          IF(NX.NE.0)AVRESX=AVRESX/NX  
111           IF( YGOOD(IP).EQ.1 )THEN           IF( YGOOD(IP).EQ.1 )THEN
112              NY=NY+1              NY=NY+1.0
113              AVRESY=AVRESY+RESY(IP)              AVRESY=AVRESY+RESY(IP)
114           ENDIF           ENDIF
          IF(NX.NE.0)AVRESY=AVRESY/NY          
115        ENDDO        ENDDO
116          IF(NX.NE.0.0)AVRESX=AVRESX/NX
117          IF(NY.NE.0.0)AVRESY=AVRESY/NY        
118    
119  *     ----------------------------------------------------------  *     ----------------------------------------------------------
120  *     define ALTOL(5) ---> tolerances on state vector  *     define ALTOL(5) ---> tolerances on state vector
# Line 144  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 153  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
153  * **** Chi2+gaussian minimization  * **** Chi2+gaussian minimization
154  * -------------------------------  * -------------------------------
155    
156        IF(.NOT.STUDENT) THEN        IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN
157    
158             IF(ISTEP.GE.3) FIRSTSTEPS = .false.
159    
160           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
161           IF(JFAIL.NE.0) THEN           IF(JFAIL.NE.0) THEN
# Line 252  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOL Line 263  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOL
263  * **** Likelihood+Student minimization  * **** Likelihood+Student minimization
264  * -------------------------------  * -------------------------------
265    
266        IF(STUDENT) THEN        IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN
267    
268             IF(FIRSTSTUDENT) THEN
269                FIRSTSTUDENT = .false.
270                ISTEP = 1
271             ENDIF
272    
273           CALL CHISQSTT(1,JFAIL)           CALL CHISQSTT(1,JFAIL)
274           DO I=1,5                                         DO I=1,5                              
275              DAL(I)=0.                        DAL(I)=0.          
# Line 286  c$$$         PRINT*,CHI2 Line 303  c$$$         PRINT*,CHI2
303           FC = CHI2           FC = CHI2
304           EC = 0.           EC = 0.
305    
306             ICOUNT = 0
307   100     CONTINUE   100     CONTINUE
308             ICOUNT = ICOUNT+1
309    
310           DO I=1,5           DO I=1,5
311              AL0(I)=AL(I)              AL0(I)=AL(I)
312           ENDDO           ENDDO
# Line 330  c$$$         PRINT*,E,CHI2_NEW Line 350  c$$$         PRINT*,E,CHI2_NEW
350              ENDIF              ENDIF
351  c$$$            E = BETA*E  c$$$            E = BETA*E
352           ENDIF           ENDIF
353             IF(ICOUNT.GT.20) GOTO 101
354           GOTO 100           GOTO 100
355    
356   101     CONTINUE   101     CONTINUE
# Line 979  cpp      ENDDO       Line 1000  cpp      ENDDO      
1000  c$$$      print*,'POSXY (prima) ',vout  c$$$      print*,'POSXY (prima) ',vout
1001    
1002        DO I=1,nplanes        DO I=1,nplanes
1003  cpp         step=vout(3)-zv(i)  c$$$         ipass = 0 ! TEST
1004           step=vout(3)-zm(i)  c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST
1005    cPPP         step=vout(3)-zm(i)
1006    cPP         step=(zm(i)-vout(3))/VOUT(6)
1007   10      DO J=1,7   10      DO J=1,7
1008              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1009              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1010           ENDDO           ENDDO
1011    cPPP         step=vect(3)-zm(i)
1012             IF(VOUT(6).GE.0.) THEN
1013                IFAIL=1
1014                if(TRKVERBOSE)      
1015         $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1016                RETURN
1017             ENDIF
1018             step=(zm(i)-vect(3))/VOUT(6)
1019   11      continue   11      continue
1020           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1021    c$$$         ipass = ipass + 1 ! TEST
1022    c$$$         PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST
1023           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
1024              IFAIL=1              IFAIL=1
1025              if(TRKVERBOSE)              if(TRKVERBOSE)
# Line 1028  c$$$            if(.TRUE.)print*,'step', Line 1061  c$$$            if(.TRUE.)print*,'step',
1061              VOUT(7) = VOUT(7) * 0.997 !0.9968              VOUT(7) = VOUT(7) * 0.997 !0.9968
1062  *     -----------------------------------------------  *     -----------------------------------------------
1063           ENDIF           ENDIF
1064    c$$$         PRINT *,'TRACKING -> END' ! TEST
1065    
1066        ENDDO        ENDDO
1067    
# Line 1066  c$$$      print*,'POSXY (dopo) ',vout Line 1100  c$$$      print*,'POSXY (dopo) ',vout
1100           YM(IP) = -100.         !0.           YM(IP) = -100.         !0.
1101           XM_A(IP) = -100.         !0.           XM_A(IP) = -100.         !0.
1102           YM_A(IP) = -100.         !0.           YM_A(IP) = -100.         !0.
1103  c         ZM_A(IP) = 0           ZM_A(IP) = fitz(nplanes-ip+1) !init to mech. position
1104           XM_B(IP) =  -100.         !0.           XM_B(IP) =  -100.         !0.
1105           YM_B(IP) =  -100.         !0.           YM_B(IP) =  -100.         !0.
1106  c         ZM_B(IP) = 0           ZM_B(IP) = fitz(nplanes-ip+1) !init to mech. position
1107           RESX(IP) = 1000.       !3.d-4           RESX(IP) = 1000.       !3.d-4
1108           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
1109           XGOOD(IP) = 0           XGOOD(IP) = 0
# Line 1218  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1252  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1252        AL(2) = Y0        AL(2) = Y0
1253        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1254        AL(3) = tath/sqrt(1+tath**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  
1255    
1256        AL(4)=0.        AL(4)=0.
1257        IF( AX.NE.0.OR.AY.NE.0. ) THEN        IF( AX.NE.0.OR.AY.NE.0. ) THEN
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23