/[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.27 by mocchiut, Fri Jan 17 12:56:51 2014 UTC
# Line 23  c      logical DEBUG Line 23  c      logical DEBUG
23  c      common/dbg/DEBUG  c      common/dbg/DEBUG
24    
25        parameter (dinf=1.d15)      !just a huge number...        parameter (dinf=1.d15)      !just a huge number...
26          parameter (dinfneg=-dinf)   ! just a huge negative number...
27  c------------------------------------------------------------------------  c------------------------------------------------------------------------
28  c     variables used in the tracking procedure (mini and its subroutines)  c     variables used in the tracking procedure (mini and its subroutines)
29  c  c
# Line 44  c      DATA XGOOD,YGOOD/nplanes*1.,nplan Line 45  c      DATA XGOOD,YGOOD/nplanes*1.,nplan
45  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
46  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
47        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
48        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"        DATA ALMIN/dinfneg,dinfneg,-1.,dinfneg,dinfneg/ !"
49    
50  c$$$      DIMENSION DAL(5)                    !increment of vector alfa  c$$$      DIMENSION DAL(5)                    !increment of vector alfa
51        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 67  c--------------------------------------- Line 68  c---------------------------------------
68                    
69  c      LOGICAL TRKDEBUG,TRKVERBOSE  c      LOGICAL TRKDEBUG,TRKVERBOSE
70  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
71        LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT        LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT
72        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
73    
74        DIMENSION AL0(5)        DIMENSION AL0(5)
75        LOGICAL SUCCESS_NEW,SUCCESS_OLD        LOGICAL SUCCESS_NEW,SUCCESS_OLD
76    
77    c$$$      PRINT*,'==========' ! TEST
78    c$$$      PRINT*,'START MINI' ! TEST
79    c$$$      PRINT*,'==========' ! TEST
80    
81  *  *
82  *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)  *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)
83  *  *
84        STUDENT = .false.        STUDENT = .false.
85          FIRSTSTEPS = .true.
86          FIRSTSTUDENT = .true.
87        IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.        IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.
88    
89        IF(IPRINT.EQ.1) THEN        IF(IPRINT.EQ.1) THEN
# Line 94  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 102  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
102  *     ----------------------------------------------------------  *     ----------------------------------------------------------
103        AVRESX = RESXAV        AVRESX = RESXAV
104        AVRESY = RESYAV        AVRESY = RESYAV
105          NX = 0 !EM GCC4.7
106          NY = 0 !EM GCC4.7
107        DO IP=1,6        DO IP=1,6
108           IF( XGOOD(IP).EQ.1 )THEN           IF( XGOOD(IP).EQ.1 )THEN
109              NX=NX+1              NX=NX+1!EM GCC4.7
110              AVRESX=AVRESX+RESX(IP)              AVRESX=AVRESX+RESX(IP)
111           ENDIF           ENDIF
          IF(NX.NE.0)AVRESX=AVRESX/NX  
112           IF( YGOOD(IP).EQ.1 )THEN           IF( YGOOD(IP).EQ.1 )THEN
113              NY=NY+1              NY=NY+1!EM GCC4.7
114              AVRESY=AVRESY+RESY(IP)              AVRESY=AVRESY+RESY(IP)
115           ENDIF           ENDIF
          IF(NX.NE.0)AVRESY=AVRESY/NY          
116        ENDDO        ENDDO
117          IF(NX.NE.0.0)AVRESX=AVRESX/NX
118          IF(NY.NE.0.0)AVRESY=AVRESY/NY        
119    
120  *     ----------------------------------------------------------  *     ----------------------------------------------------------
121  *     define ALTOL(5) ---> tolerances on state vector  *     define ALTOL(5) ---> tolerances on state vector
# Line 132  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 142  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
142        CHI2=0        CHI2=0
143    
144        if(TRKDEBUG) print*,'guess: ',al        if(TRKDEBUG) print*,'guess: ',al
145        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
146    
147  *      *    
148  *     -----------------------  *     -----------------------
# Line 144  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 154  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
154  * **** Chi2+gaussian minimization  * **** Chi2+gaussian minimization
155  * -------------------------------  * -------------------------------
156    
157        IF(.NOT.STUDENT) THEN        IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN
158    
159             IF(ISTEP.GE.3) FIRSTSTEPS = .false.
160    
161           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
162           IF(JFAIL.NE.0) THEN           IF(JFAIL.NE.0) THEN
# Line 240  c$$$            CHI2D(I)=CHI2D(I)*COST Line 252  c$$$            CHI2D(I)=CHI2D(I)*COST
252              ENDDO              ENDDO
253           ENDIF           ENDIF
254                
255           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
256    
257  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
258  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
# Line 252  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOL Line 264  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOL
264  * **** Likelihood+Student minimization  * **** Likelihood+Student minimization
265  * -------------------------------  * -------------------------------
266    
267        IF(STUDENT) THEN        IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN
268    
269             IF(FIRSTSTUDENT) THEN
270                FIRSTSTUDENT = .false.
271                ISTEP = 1
272             ENDIF
273    
274           CALL CHISQSTT(1,JFAIL)           CALL CHISQSTT(1,JFAIL)
275           DO I=1,5                                         DO I=1,5                              
276              DAL(I)=0.                        DAL(I)=0.          
# Line 286  c$$$         PRINT*,CHI2 Line 304  c$$$         PRINT*,CHI2
304           FC = CHI2           FC = CHI2
305           EC = 0.           EC = 0.
306    
307             ICOUNT = 0
308   100     CONTINUE   100     CONTINUE
309             ICOUNT = ICOUNT+1
310    
311           DO I=1,5           DO I=1,5
312              AL0(I)=AL(I)              AL0(I)=AL(I)
313           ENDDO           ENDDO
# Line 330  c$$$         PRINT*,E,CHI2_NEW Line 351  c$$$         PRINT*,E,CHI2_NEW
351              ENDIF              ENDIF
352  c$$$            E = BETA*E  c$$$            E = BETA*E
353           ENDIF           ENDIF
354             IF(ICOUNT.GT.20) GOTO 101
355           GOTO 100           GOTO 100
356    
357   101     CONTINUE   101     CONTINUE
# Line 397  c$$$     $        ISTEPMAX Line 419  c$$$     $        ISTEPMAX
419  *     ---------------------------------------------  *     ---------------------------------------------
420  *------------------------------------------------------------*  *------------------------------------------------------------*
421  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
422          IF(FACT.EQ.0)THEN
423             IFAIL=1
424             RETURN
425          ENDIF
426        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
427        ALTOL(1) = ALTOL(5)/DELETA1        ALTOL(1) = ALTOL(5)/DELETA1
428        ALTOL(2) = ALTOL(1)        ALTOL(2) = ALTOL(1)
# Line 533  c$$$               DAL(I)=0. Line 559  c$$$               DAL(I)=0.
559  *     ------------------------------------  *     ------------------------------------
560  *     Reduced chi^2  *     Reduced chi^2
561        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
   
562  c      print*,'mini2: chi2 ',chi2  c      print*,'mini2: chi2 ',chi2
563    
564   11   CONTINUE         11   CONTINUE      
565    
566        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,AL(5)
567    
568        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
569    
# Line 979  cpp      ENDDO       Line 1004  cpp      ENDDO      
1004  c$$$      print*,'POSXY (prima) ',vout  c$$$      print*,'POSXY (prima) ',vout
1005    
1006        DO I=1,nplanes        DO I=1,nplanes
1007  cpp         step=vout(3)-zv(i)  c$$$         ipass = 0 ! TEST
1008           step=vout(3)-zm(i)  c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST
1009    cPPP         step=vout(3)-zm(i)
1010    cPP         step=(zm(i)-vout(3))/VOUT(6)
1011   10      DO J=1,7   10      DO J=1,7
1012              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1013              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1014           ENDDO           ENDDO
1015    cPPP         step=vect(3)-zm(i)
1016             IF(VOUT(6).GE.0.) THEN
1017                IFAIL=1
1018                if(TRKVERBOSE)      
1019         $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1020                RETURN
1021             ENDIF
1022             step=(zm(i)-vect(3))/VOUT(6)
1023   11      continue   11      continue
1024           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1025    c$$$         ipass = ipass + 1 ! TEST
1026    c$$$         PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST
1027           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
1028              IFAIL=1              IFAIL=1
1029              if(TRKVERBOSE)              if(TRKVERBOSE)
# Line 1028  c$$$            if(.TRUE.)print*,'step', Line 1065  c$$$            if(.TRUE.)print*,'step',
1065              VOUT(7) = VOUT(7) * 0.997 !0.9968              VOUT(7) = VOUT(7) * 0.997 !0.9968
1066  *     -----------------------------------------------  *     -----------------------------------------------
1067           ENDIF           ENDIF
1068    c$$$         PRINT *,'TRACKING -> END' ! TEST
1069    
1070        ENDDO        ENDDO
1071    
# Line 1066  c$$$      print*,'POSXY (dopo) ',vout Line 1104  c$$$      print*,'POSXY (dopo) ',vout
1104           YM(IP) = -100.         !0.           YM(IP) = -100.         !0.
1105           XM_A(IP) = -100.         !0.           XM_A(IP) = -100.         !0.
1106           YM_A(IP) = -100.         !0.           YM_A(IP) = -100.         !0.
1107  c         ZM_A(IP) = 0           ZM_A(IP) = fitz(nplanes-ip+1) !init to mech. position
1108           XM_B(IP) =  -100.         !0.           XM_B(IP) =  -100.         !0.
1109           YM_B(IP) =  -100.         !0.           YM_B(IP) =  -100.         !0.
1110  c         ZM_B(IP) = 0           ZM_B(IP) = fitz(nplanes-ip+1) !init to mech. position
1111           RESX(IP) = 1000.       !3.d-4           RESX(IP) = 1000.       !3.d-4
1112           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
1113           XGOOD(IP) = 0           XGOOD(IP) = 0
# Line 1117  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 1155  c      IMPLICIT DOUBLE PRECISION (A-H,O-
1155        S1=0.        S1=0.
1156        DO I=1,nplanes        DO I=1,nplanes
1157           IF(YGOOD(I).EQ.1)THEN           IF(YGOOD(I).EQ.1)THEN
1158              YY = YM(I)              YY = REAL(YM(I))!EM GCC4.7
1159              IF(XGOOD(I).EQ.0)THEN              IF(XGOOD(I).EQ.0)THEN
1160                 YY = (YM_A(I) + YM_B(I))/2                 YY = REAL((YM_A(I) + YM_B(I))/2.)!EM GCC4.7
1161              ENDIF              ENDIF
1162              SZZ=SZZ+ZM(I)*ZM(I)              SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7
1163              SZY=SZY+ZM(I)*YY              SZY=SZY+REAL(ZM(I)*YY)!EM GCC4.7
1164              SSY=SSY+YY              SSY=SSY+YY
1165              SZ=SZ+ZM(I)              SZ=SZ+REAL(ZM(I))!EM GCC4.7
1166              S1=S1+1.              S1=S1+1.
1167           ENDIF           ENDIF
1168        ENDDO        ENDDO
1169        DET=SZZ*S1-SZ*SZ        DET=SZZ*S1-SZ*SZ
1170        AY=(SZY*S1-SZ*SSY)/DET        AY=(SZY*S1-SZ*SSY)/DET
1171        BY=(SZZ*SSY-SZY*SZ)/DET        BY=(SZZ*SSY-SZY*SZ)/DET
1172        Y0 = AY*ZINI+BY        Y0 = REAL(AY*ZINI+BY)!EM GCC4.7
1173  *     ----------------------------------------  *     ----------------------------------------
1174  *     X view  *     X view
1175  *     ----------------------------------------  *     ----------------------------------------
# Line 1141  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 1179  c      IMPLICIT DOUBLE PRECISION (A-H,O-
1179        NP=0        NP=0
1180        DO I=1,nplanes        DO I=1,nplanes
1181           IF(XGOOD(I).EQ.1)THEN           IF(XGOOD(I).EQ.1)THEN
1182              XX = XM(I)              XX = REAL(XM(I))!EM GCC4.7
1183              IF(YGOOD(I).EQ.0)THEN              IF(YGOOD(I).EQ.0)THEN
1184                 XX = (XM_A(I) + XM_B(I))/2                 XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7
1185              ENDIF              ENDIF
1186              NP=NP+1              NP=NP+1
1187              XP(NP)=XX              XP(NP)=XX
1188              ZP(NP)=ZM(I)              ZP(NP)=REAL(ZM(I))!EM GCC4.7
1189           ENDIF           ENDIF
1190        ENDDO        ENDDO
1191        IFLAG=0                   !no debug mode        IFLAG=0                   !no debug mode
# Line 1161  c$$$      print*,' XP ',(rp(i),i=1,np) Line 1199  c$$$      print*,' XP ',(rp(i),i=1,np)
1199    
1200        IF(IFLAG.NE.0)GOTO 10 !straigth fit        IF(IFLAG.NE.0)GOTO 10 !straigth fit
1201  c      if(CHI.gt.100)GOTO 10 !straigth fit  c      if(CHI.gt.100)GOTO 10 !straigth fit
1202        ARG = RADIUS**2-(ZINI-ZC)**2        ARG = REAL(RADIUS**2-(ZINI-ZC)**2)!EM GCC4.7
1203        IF(ARG.LT.0)GOTO 10       !straigth fit        IF(ARG.LT.0)GOTO 10       !straigth fit
1204        DC = SQRT(ARG)              DC = SQRT(ARG)      
1205        IF(XC.GT.0)DC=-DC        IF(XC.GT.0)DC=-DC
1206        X0=XC+DC        X0=XC+DC
1207        AX = -(ZINI-ZC)/DC        AX = REAL(-(ZINI-ZC)/DC)!EM GCC4.7
1208        DEF=100./(RADIUS*0.3*0.43)        DEF=100./(RADIUS*0.3*0.43)
1209        IF(XC.GT.0)DEF=-DEF        IF(XC.GT.0)DEF=-DEF
1210                
# Line 1192  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1230  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1230        S1=0.        S1=0.
1231        DO I=1,nplanes        DO I=1,nplanes
1232           IF(XGOOD(I).EQ.1)THEN           IF(XGOOD(I).EQ.1)THEN
1233              XX = XM(I)              XX = REAL(XM(I))!EM GCC4.7
1234              IF(YGOOD(I).EQ.0)THEN              IF(YGOOD(I).EQ.0)THEN
1235                 XX = (XM_A(I) + XM_B(I))/2                 XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7
1236              ENDIF              ENDIF
1237              SZZ=SZZ+ZM(I)*ZM(I)              SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7
1238              SZX=SZX+ZM(I)*XX              SZX=SZX+REAL(ZM(I)*XX)!EM GCC4.7
1239              SSX=SSX+XX              SSX=SSX+XX
1240              SZ=SZ+ZM(I)              SZ=SZ+REAL(ZM(I))!EM GCC4.7
1241              S1=S1+1.              S1=S1+1.
1242           ENDIF           ENDIF
1243        ENDDO        ENDDO
# Line 1207  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1245  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1245        AX=(SZX*S1-SZ*SSX)/DET        AX=(SZX*S1-SZ*SSX)/DET
1246        BX=(SZZ*SSX-SZX*SZ)/DET        BX=(SZZ*SSX-SZX*SZ)/DET
1247        DEF = 0        DEF = 0
1248        X0  = AX*ZINI+BX        X0  = REAL(AX*ZINI+BX)!EM GCC4.7
1249    
1250   20   CONTINUE   20   CONTINUE
1251  *     ----------------------------------------  *     ----------------------------------------
# Line 1218  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1256  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1256        AL(2) = Y0        AL(2) = Y0
1257        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1258        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  
1259    
1260        AL(4)=0.        AL(4)=0.
1261        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.27
  ViewVC Help
Powered by ViewVC 1.1.23