/[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.18 by pam-fi, Fri Jun 1 15:01:19 2007 UTC revision 1.26 by mocchiut, Thu Jan 16 15:29:56 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    
28          double precision NX, NY ! EM GCC4.7
29  c------------------------------------------------------------------------  c------------------------------------------------------------------------
30  c     variables used in the tracking procedure (mini and its subroutines)  c     variables used in the tracking procedure (mini and its subroutines)
31  c  c
# Line 44  c      DATA XGOOD,YGOOD/nplanes*1.,nplan Line 47  c      DATA XGOOD,YGOOD/nplanes*1.,nplan
47  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
48  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
49        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
50        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"        DATA ALMIN/dinfneg,dinfneg,-1.,dinfneg,dinfneg/ !"
51    
52  c$$$      DIMENSION DAL(5)                    !increment of vector alfa  c$$$      DIMENSION DAL(5)                    !increment of vector alfa
53        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 72  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 75  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
75    
76        DIMENSION AL0(5)        DIMENSION AL0(5)
77        LOGICAL SUCCESS_NEW,SUCCESS_OLD        LOGICAL SUCCESS_NEW,SUCCESS_OLD
78    
79    c$$$      PRINT*,'==========' ! TEST
80    c$$$      PRINT*,'START MINI' ! TEST
81    c$$$      PRINT*,'==========' ! TEST
82    
83  *  *
84  *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)  *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)
85  *  *
# Line 96  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 104  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
104  *     ----------------------------------------------------------  *     ----------------------------------------------------------
105        AVRESX = RESXAV        AVRESX = RESXAV
106        AVRESY = RESYAV        AVRESY = RESYAV
107          NX = 0.0
108          NY = 0.0
109        DO IP=1,6        DO IP=1,6
110           IF( XGOOD(IP).EQ.1 )THEN           IF( XGOOD(IP).EQ.1 )THEN
111              NX=NX+1              NX=NX+1.0
112              AVRESX=AVRESX+RESX(IP)              AVRESX=AVRESX+RESX(IP)
113           ENDIF           ENDIF
          IF(NX.NE.0)AVRESX=AVRESX/NX  
114           IF( YGOOD(IP).EQ.1 )THEN           IF( YGOOD(IP).EQ.1 )THEN
115              NY=NY+1              NY=NY+1.0
116              AVRESY=AVRESY+RESY(IP)              AVRESY=AVRESY+RESY(IP)
117           ENDIF           ENDIF
          IF(NX.NE.0)AVRESY=AVRESY/NY          
118        ENDDO        ENDDO
119          IF(NX.NE.0.0)AVRESX=AVRESX/NX
120          IF(NY.NE.0.0)AVRESY=AVRESY/NY        
121    
122  *     ----------------------------------------------------------  *     ----------------------------------------------------------
123  *     define ALTOL(5) ---> tolerances on state vector  *     define ALTOL(5) ---> tolerances on state vector
# Line 134  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 144  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
144        CHI2=0        CHI2=0
145    
146        if(TRKDEBUG) print*,'guess: ',al        if(TRKDEBUG) print*,'guess: ',al
147        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
148    
149  *      *    
150  *     -----------------------  *     -----------------------
# Line 146  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 156  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
156  * **** Chi2+gaussian minimization  * **** Chi2+gaussian minimization
157  * -------------------------------  * -------------------------------
158    
159        IF(.NOT.STUDENT.OR.FIRSTSTEPS) THEN        IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN
160    
161           IF(ISTEP.GE.3) FIRSTSTEPS = .false.           IF(ISTEP.GE.3) FIRSTSTEPS = .false.
162    
# Line 244  c$$$            CHI2D(I)=CHI2D(I)*COST Line 254  c$$$            CHI2D(I)=CHI2D(I)*COST
254              ENDDO              ENDDO
255           ENDIF           ENDIF
256                
257           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
258    
259  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
260  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
# Line 296  c$$$         PRINT*,CHI2 Line 306  c$$$         PRINT*,CHI2
306           FC = CHI2           FC = CHI2
307           EC = 0.           EC = 0.
308    
309             ICOUNT = 0
310   100     CONTINUE   100     CONTINUE
311             ICOUNT = ICOUNT+1
312    
313           DO I=1,5           DO I=1,5
314              AL0(I)=AL(I)              AL0(I)=AL(I)
315           ENDDO           ENDDO
# Line 340  c$$$         PRINT*,E,CHI2_NEW Line 353  c$$$         PRINT*,E,CHI2_NEW
353              ENDIF              ENDIF
354  c$$$            E = BETA*E  c$$$            E = BETA*E
355           ENDIF           ENDIF
356             IF(ICOUNT.GT.20) GOTO 101
357           GOTO 100           GOTO 100
358    
359   101     CONTINUE   101     CONTINUE
# Line 407  c$$$     $        ISTEPMAX Line 421  c$$$     $        ISTEPMAX
421  *     ---------------------------------------------  *     ---------------------------------------------
422  *------------------------------------------------------------*  *------------------------------------------------------------*
423  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
424          IF(FACT.EQ.0)THEN
425             IFAIL=1
426             RETURN
427          ENDIF
428        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
429        ALTOL(1) = ALTOL(5)/DELETA1        ALTOL(1) = ALTOL(5)/DELETA1
430        ALTOL(2) = ALTOL(1)        ALTOL(2) = ALTOL(1)
# Line 543  c$$$               DAL(I)=0. Line 561  c$$$               DAL(I)=0.
561  *     ------------------------------------  *     ------------------------------------
562  *     Reduced chi^2  *     Reduced chi^2
563        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
   
564  c      print*,'mini2: chi2 ',chi2  c      print*,'mini2: chi2 ',chi2
565    
566   11   CONTINUE         11   CONTINUE      
567    
568        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,AL(5)
569    
570        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
571    
# Line 989  cpp      ENDDO       Line 1006  cpp      ENDDO      
1006  c$$$      print*,'POSXY (prima) ',vout  c$$$      print*,'POSXY (prima) ',vout
1007    
1008        DO I=1,nplanes        DO I=1,nplanes
1009  cpp         step=vout(3)-zv(i)  c$$$         ipass = 0 ! TEST
1010           step=vout(3)-zm(i)  c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST
1011    cPPP         step=vout(3)-zm(i)
1012    cPP         step=(zm(i)-vout(3))/VOUT(6)
1013   10      DO J=1,7   10      DO J=1,7
1014              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1015              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1016           ENDDO           ENDDO
1017    cPPP         step=vect(3)-zm(i)
1018             IF(VOUT(6).GE.0.) THEN
1019                IFAIL=1
1020                if(TRKVERBOSE)      
1021         $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1022                RETURN
1023             ENDIF
1024             step=(zm(i)-vect(3))/VOUT(6)
1025   11      continue   11      continue
1026           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1027    c$$$         ipass = ipass + 1 ! TEST
1028    c$$$         PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST
1029           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
1030              IFAIL=1              IFAIL=1
1031              if(TRKVERBOSE)              if(TRKVERBOSE)
# Line 1038  c$$$            if(.TRUE.)print*,'step', Line 1067  c$$$            if(.TRUE.)print*,'step',
1067              VOUT(7) = VOUT(7) * 0.997 !0.9968              VOUT(7) = VOUT(7) * 0.997 !0.9968
1068  *     -----------------------------------------------  *     -----------------------------------------------
1069           ENDIF           ENDIF
1070    c$$$         PRINT *,'TRACKING -> END' ! TEST
1071    
1072        ENDDO        ENDDO
1073    
# Line 1076  c$$$      print*,'POSXY (dopo) ',vout Line 1106  c$$$      print*,'POSXY (dopo) ',vout
1106           YM(IP) = -100.         !0.           YM(IP) = -100.         !0.
1107           XM_A(IP) = -100.         !0.           XM_A(IP) = -100.         !0.
1108           YM_A(IP) = -100.         !0.           YM_A(IP) = -100.         !0.
1109  c         ZM_A(IP) = 0           ZM_A(IP) = fitz(nplanes-ip+1) !init to mech. position
1110           XM_B(IP) =  -100.         !0.           XM_B(IP) =  -100.         !0.
1111           YM_B(IP) =  -100.         !0.           YM_B(IP) =  -100.         !0.
1112  c         ZM_B(IP) = 0           ZM_B(IP) = fitz(nplanes-ip+1) !init to mech. position
1113           RESX(IP) = 1000.       !3.d-4           RESX(IP) = 1000.       !3.d-4
1114           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
1115           XGOOD(IP) = 0           XGOOD(IP) = 0
# Line 1107  c         ZM_B(IP) = 0 Line 1137  c         ZM_B(IP) = 0
1137    
1138        subroutine guess()        subroutine guess()
1139    
1140  c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z) ! EM GCC4.7
1141    
1142        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
1143        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
1144                
1145        REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)        REAL*8 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES) ! EM GCC4.7
1146        REAL*4 CHI,XC,ZC,RADIUS        REAL*4 CHI,XC,ZC,RADIUS
1147  *     ----------------------------------------  *     ----------------------------------------
1148  *     Y view  *     Y view
# Line 1228  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1258  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1258        AL(2) = Y0        AL(2) = Y0
1259        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1260        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  
1261    
1262        AL(4)=0.        AL(4)=0.
1263        IF( AX.NE.0.OR.AY.NE.0. ) THEN        IF( AX.NE.0.OR.AY.NE.0. ) THEN
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.23